Hydra  4.0.1
A header-only templated C++ framework to perform data analysis on massively parallel platforms.
phsp_reweighting.inl

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.

/*----------------------------------------------------------------------------
*
* Copyright (C) 2016 - 2023 Antonio Augusto Alves Junior
*
* This file is part of Hydra Data Analysis Framework.
*
* Hydra is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Hydra is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Hydra. If not, see <http://www.gnu.org/licenses/>.
*
*---------------------------------------------------------------------------*/
/*
* phsp_basic.inl
*
* Created on: Jul 7, 2017
* Author: Antonio Augusto Alves Junior
*/
#ifndef PHSP_REWEIGHTING_INL_
#define PHSP_REWEIGHTING_INL_
/**
* \example phsp_reweighting.inl
* 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.
*/
/*---------------------------------
* std
* ---------------------------------
*/
#include <iostream>
#include <assert.h>
#include <time.h>
#include <array>
#include <chrono>
/*---------------------------------
* command line arguments
*---------------------------------
*/
#include <tclap/CmdLine.h>
/*---------------------------------
* Include hydra classes and
* algorithms for
*--------------------------------
*/
#include <hydra/Types.h>
#include <hydra/Vector4R.h>
#include <hydra/Function.h>
#include <hydra/Lambda.h>
#include <hydra/Tuple.h>
#include <hydra/Decays.h>
#include <hydra/Range.h>
/*-------------------------------------
* Include classes from ROOT to fill
* and draw histograms and plots.
*-------------------------------------
*/
#ifdef _ROOT_AVAILABLE_
#include <TROOT.h>
#include <TH1D.h>
#include <TH2D.h>
#include <TApplication.h>
#include <TCanvas.h>
#include <TColor.h>
#include <TString.h>
#include <TStyle.h>
#endif //_ROOT_AVAILABLE_
// Daughter particles
declarg(A, hydra::Vector4R)
declarg(B, hydra::Vector4R)
declarg(C, hydra::Vector4R)
//---------------------------
using namespace hydra::arguments;
int main(int argv, char** argc)
{
size_t nentries = 0; // number of events to generate, to be get from command line
double P_mass = 5.27955;
double A_mass = 3.0969;
double B_mass = 0.493677;
double C_mass = 0.13957061;
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");
cmd.add(NArg);
// Parse the argv array.
cmd.parse(argv, argc);
// Get the value parsed by each arg.
nentries = NArg.getValue();
}
catch (TCLAP::ArgException &e) {
std::cerr << "error: " << e.error() << " for arg " << e.argId()
<< std::endl;
}
#ifdef _ROOT_AVAILABLE_
//
TH2D Dalitz_Flat("Dalitz_Flat",
"PHSP;"
"M^{2}(AB) [GeV^{2}/c^{4}];"
"M^{2}(BC) [GeV^{2}/c^{4}]",
100, pow(A_mass + B_mass,2), pow(P_mass - C_mass,2),
100, pow(B_mass + C_mass,2), pow(P_mass - A_mass,2));
//
TH2D Dalitz_Resonant("Dalitz_Resonant",
"Resonance;"
"M^{2}(AB) [GeV^{2}/c^{4}];"
"M^{2}(BC) [GeV^{2}/c^{4}]",
100, pow(A_mass + B_mass,2), pow(P_mass - C_mass,2),
100, pow(B_mass + C_mass,2), pow(P_mass - A_mass,2));
#endif
hydra::Vector4R Parent(P_mass, 0.0, 0.0, 0.0);
double masses[3]{A_mass, B_mass, C_mass };
// Create PhaseSpace object for B0-> K pi J/psi
return hydra::make_tuple( (a + b).mass2(), (b + c).mass2());
});
auto breit_wigner = hydra::wrap_lambda( [ resonance ] __hydra_dual__ ( B b, C c ){
return resonance( (b+c).mass());
} );
//device
{
//allocate memory to hold the final states particles
auto start = std::chrono::high_resolution_clock::now();
//generate the final state particles
phsp.Generate(Parent, Events);
auto end = std::chrono::high_resolution_clock::now();
std::chrono::duration<double, std::milli> elapsed = end - start;
//output
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;
//the power of lazyness
auto dalitz_variables = Events | dalitz_calculator ;
auto dalitz_flat_weights = Events | Events.GetEventWeightFunctor();
auto dalitz_resonant_weights = Events | Events.GetEventWeightFunctor(breit_wigner);
auto Hist_Dalitz_Flat = hydra::make_dense_histogram<double,2>( hydra::device::sys,
{100,100},
{pow(A_mass + B_mass,2), pow(B_mass + C_mass,2)},
{pow(P_mass - C_mass,2), pow(P_mass - A_mass,2)},
dalitz_variables, dalitz_flat_weights);
auto Hist_Dalitz_Resonant = hydra::make_dense_histogram<double,2>( hydra::device::sys,
{100,100},
{pow(A_mass + B_mass,2), pow(B_mass + C_mass,2)},
{pow(P_mass - C_mass,2), pow(P_mass - A_mass,2)},
dalitz_variables, dalitz_resonant_weights);
#ifdef _ROOT_AVAILABLE_
//bring events to CPU memory space
//hydra::Events<3, hydra::host::sys_t > Events_h(Events_d);
for(size_t i=0; i< 100; i++){
for(size_t j=0; j< 100; j++){
Dalitz_Flat.SetBinContent(i+1, j+1, Hist_Dalitz_Flat.GetBinContent({i,j}) );
Dalitz_Resonant.SetBinContent(i+1, j+1, Hist_Dalitz_Resonant.GetBinContent({i,j}) );
}
}
#endif
}
#ifdef _ROOT_AVAILABLE_
TApplication *m_app=new TApplication("myapp",0,0);
TCanvas canvas_d1("canvas_d1", "Phase-space Device", 500, 500);
Dalitz_Flat.Draw("colz");
TCanvas canvas_d2("canvas_d2", "Phase-space Device", 500, 500);
Dalitz_Resonant.Draw("colz");
m_app->Run();
#endif
return 0;
}
#endif /* PHSP_REWEIGHTING_INL_ */