This example shows how to use the Hydra's phase space Monte Carlo algorithms to generate a sample of B0 -> J/psi K pi and plot the Dalitz plot.
#ifndef PHSP_UNWEIGHTING_FUNCTOR_INL_
#define PHSP_UNWEIGHTING_FUNCTOR_INL_
#include <iostream>
#include <assert.h>
#include <time.h>
#include <vector>
#include <array>
#include <chrono>
#include <tclap/CmdLine.h>
#ifdef _ROOT_AVAILABLE_
#include <TROOT.h>
#include <TH1D.h>
#include <TH2D.h>
#include <TApplication.h>
#include <TCanvas.h>
#include <TStyle.h>
#endif //_ROOT_AVAILABLE_
{
try {
TCLAP::CmdLine cmd("Command line arguments for PHSP B0 -> J/psi K pi", '=');
TCLAP::ValueArg<size_t>
NArg(
"n",
"nevents",
"Number of events to generate. Default is [ 10e6 ].",
true, 10e6, "unsigned long");
}
catch (TCLAP::ArgException &e) {
std::cerr << "error: " << e.error() << " for arg " << e.argId()
<< std::endl;
}
#ifdef _ROOT_AVAILABLE_
TH2D Dalitz_W("Dalitz_W",
"Weighted Sample;"
"M^{2}(A B) [GeV^{2}/c^{4}];"
"M^{2}(B C) [GeV^{2}/c^{4}]",
TH2D Dalitz_U("Dalitz_U",
"Unweighted Sample;"
"M^{2}(A B) [GeV^{2}/c^{4}];"
"M^{2}(B C) [GeV^{2}/c^{4}]",
#endif
});
return resonance( (b+c).
mass());
} );
{
auto start = std::chrono::high_resolution_clock::now();
auto end = std::chrono::high_resolution_clock::now();
std::cout << std::endl;
std::cout << std::endl;
std::cout << "----------------- Device ----------------"<< std::endl;
std::cout << "| P -> A B C" << std::endl;
std::cout <<
"| Number of events :"<<
nentries << std::endl;
std::cout <<
"| Time (ms) :"<<
elapsed.count() << std::endl;
std::cout << "-----------------------------------------"<< std::endl;
auto dalitz_weights = Events | Events.GetEventWeightFunctor();
{100,100},
dalitz_variables, dalitz_weights);
{100,100},
dalitz_unweighted);
#ifdef _ROOT_AVAILABLE_
for(size_t i=0; i< 100; i++){
for(size_t j=0; j< 100; j++){
Dalitz_W.SetBinContent(i+1, j+1, Hist_Dalitz_W.GetBinContent({i,j}) );
Dalitz_U.SetBinContent(i+1, j+1, Hist_Dalitz_U.GetBinContent({i,j}) );
}
}
#endif
}
#ifdef _ROOT_AVAILABLE_
TApplication *m_app=new TApplication("myapp",0,0);
TCanvas canvas_d1("canvas_d1", "Phase-space weigted sample", 500, 500);
Dalitz_W.Draw("colz");
TCanvas canvas_d2("canvas_d2", "Phase-space unweigted sample", 500, 500);
Dalitz_U.Draw("colz");
m_app->Run();
#endif
return 0;
}
#endif