From 1363673db3ae34b87d8b71dfe3b8077be4b171b3 Mon Sep 17 00:00:00 2001 From: Axel Naumann Date: Thu, 4 Mar 2021 16:41:29 +0100 Subject: [PATCH 1/3] =?UTF-8?q?[physics]=20Add=20T{Gen{Foam}}Decay.=20By?= =?UTF-8?q?=20Rados=C5=82aw=20Kycia!?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- math/physics/CMakeLists.txt | 10 ++ math/physics/inc/TDecay.h | 118 +++++++++++++ math/physics/inc/TGenDecay.h | 126 ++++++++++++++ math/physics/inc/TGenFoamDecay.h | 201 ++++++++++++++++++++++ math/physics/inc/TGenInterface.h | 94 +++++++++++ math/physics/src/TDecay.cxx | 218 ++++++++++++++++++++++++ math/physics/src/TGenDecay.cxx | 109 ++++++++++++ math/physics/src/TGenFoamDecay.cxx | 180 ++++++++++++++++++++ math/physics/test/CMakeLists.txt | 11 ++ math/physics/test/tTDecay.cxx | 240 +++++++++++++++++++++++++++ math/physics/test/tTGenDecay.cxx | 211 +++++++++++++++++++++++ math/physics/test/tTGenFoamDecay.cxx | 211 +++++++++++++++++++++++ 12 files changed, 1729 insertions(+) create mode 100644 math/physics/inc/TDecay.h create mode 100644 math/physics/inc/TGenDecay.h create mode 100644 math/physics/inc/TGenFoamDecay.h create mode 100644 math/physics/inc/TGenInterface.h create mode 100644 math/physics/src/TDecay.cxx create mode 100644 math/physics/src/TGenDecay.cxx create mode 100644 math/physics/src/TGenFoamDecay.cxx create mode 100644 math/physics/test/CMakeLists.txt create mode 100644 math/physics/test/tTDecay.cxx create mode 100644 math/physics/test/tTGenDecay.cxx create mode 100644 math/physics/test/tTGenFoamDecay.cxx diff --git a/math/physics/CMakeLists.txt b/math/physics/CMakeLists.txt index c18b861f7281c..a03aba5095c04 100644 --- a/math/physics/CMakeLists.txt +++ b/math/physics/CMakeLists.txt @@ -11,6 +11,10 @@ ROOT_STANDARD_LIBRARY_PACKAGE(Physics HEADERS TFeldmanCousins.h + TDecay.h + TGenDecay.h + TGenFoamDecay.h + TGenInterface.h TGenPhaseSpace.h TLorentzRotation.h TLorentzVector.h @@ -22,6 +26,9 @@ ROOT_STANDARD_LIBRARY_PACKAGE(Physics TVector3.h SOURCES src/TFeldmanCousins.cxx + src/TDecay.cxx + src/TGenDecay.cxx + src/TGenFoamDecay.cxx src/TGenPhaseSpace.cxx src/TLorentzRotation.cxx src/TLorentzVector.cxx @@ -32,8 +39,11 @@ ROOT_STANDARD_LIBRARY_PACKAGE(Physics src/TVector2.cxx src/TVector3.cxx DEPENDENCIES + Foam Matrix MathCore DICTIONARY_OPTIONS -writeEmptyRootPCM ) + +ROOT_ADD_TEST_SUBDIRECTORY(test) diff --git a/math/physics/inc/TDecay.h b/math/physics/inc/TDecay.h new file mode 100644 index 0000000000000..4e2ce55e14d1e --- /dev/null +++ b/math/physics/inc/TDecay.h @@ -0,0 +1,118 @@ +/************************************************************************* +* Copyright (C) 1995-2021, Rene Brun and Fons Rademakers. * +* All rights reserved. * +* * +* For the licensing terms see $ROOTSYS/LICENSE. * +* For the list of contributors see $ROOTSYS/README/CREDITS. * +*************************************************************************/ + +/*! + @class TDecay + + @brief The generator decays a central blob of a given 4-momentum into particles of specific masses. Original GENBOD algorithm with w few modifications. + +### Authors + +Radosław Kycia (kycia.radoslaw@gmail.com), Piotr Lebiedowicz, Antoni Szczurek + +People who helped in the development of the project: Jacek Turnau, Janusz Chwastowski, Rafał Staszewski, Maciej Trzebiński. + +### Description + +This is a helper class that can be base to construct MC generators that simulate decay. It returns weight equals LIPS (Lorentz Invariant Phase Space): + +\f$ dLIPS = (2\pi)^4 \delta^{(4)}(P-\sum_{i=1}^{n}p_{i}) \prod_{i=1}^{n} \frac{d^{3}p_{i}}{(2\pi)^3 2E_{i}}\f$ + +The generator can be used as a building block in more complicated decays. + +The generator returns weighted events. + +The class is adapted from the TGenPhaseSpace ROOT package. + +This class is not recommended for non-advanced users. If you want a robust and adaptive MC generator and integrator, use TGenFoamDecay instead. + +#### The scheme of use + +1. Prepare 4-momenta of decaying particle P and mass[nt] array of final products. +2. Initialize generator: SetDecay(); +3. Generate decay: Generate(); The class requires random numbers given to Generate(). They are used to make an event. +4. For each of particles enumerated by 0..nt-1 get their 4-momentum using GetDecay(): pfi = GetDecay(i); +5. Repeat 3 and 4 for another decay. + + +### Further details + +- R.A. Kycia, P. Lebiedowicz, A. Szczurek, 'Decay: A Monte Carlo library for the decay of a particle with ROOT compatibility', arXiv:2011.14750 [hep-ph] +- R. A. Kycia, J. Chwastowski, R. Staszewski, and J. Turnau, 'GenEx: A simple generator structure for exclusive processes in high energy collisions', Commun. Comput. Phys.24no. 3, (2018) 860, arXiv:1411.6035 [hep-ph] +- R. A. Kycia, J. Turnau, J. J. Chwastowski, R. Staszewski, and M. Trzebinski, 'The adaptiveMonte Carlo toolbox for phase space integration and generation', Commun. Comput. Phys.25no. 5, (2019) 1547, arXiv:1711.06087 [hep-ph] +- R. Hagedorn, 'Relativistic Kinematics: A guide to the kinematic problems of high-energy physics'. W.A. Benjamin, New York, Amsterdam, 1964 +- H.Pilkhun, The Interactions of Hadrons North-Holland 1967 + +*/ + + +#ifndef ROOT_TDecay +#define ROOT_TDecay + +#include + +#include // std::cin, std::cout +#include // std::queue +#include + +using namespace std; + + +class TDecay : public TObject { +private: + Int_t fNt; // number of decay particles + Double_t fMass[18]; // masses of particles + Double_t fBeta[3]; // betas of decaying particle + Double_t fTeCmTm; // total energy in the C.M. minus the total mass + TLorentzVector fDecPro[18]; //kinematics of the generated particles + + Int_t kMAXP; //max number of particles (relict from TGenPhaseSpace) + + Double_t PDK(Double_t a, Double_t b, Double_t c); + + ///factorial function + int factorial(int n); + +public: + + /// constructor + TDecay(): fNt(0), fMass(), fBeta(), fTeCmTm(0.) {}; + + /// copy constructor + TDecay(const TDecay &gen); + + /// desctructor + virtual ~TDecay() {}; + + /// assignment + TDecay& operator=(const TDecay &gen); + + /// Sets up configuration of decay + /// @param P 4-momentum of decaying particle + /// @param nt number of products of decay + /// @param mass mass matrix of products of decay mass[nt] + /// @returns kTRUE - the decay is permitted by kinematics; kFALSE - the decay is forbidden by kinematics + /// @warning This should be first method to call since it sets up decay configuration. + Bool_t SetDecay(TLorentzVector &P, Int_t nt, const Double_t *mass); + + /// Generate a single decay + /// @param rnd a queue of 3*nt-4 random numbers from external source. + Double_t Generate(std::queue &rnd); + + /// @returns 4-vector of n-th product of decay + /// @param n number of final particle to get from range 1...nt + /// @warning You should call Generate() in first place. + TLorentzVector *GetDecay(Int_t n); + + /// @returns number of final particles + Int_t GetNt() const { return fNt;} + +}; + +#endif + diff --git a/math/physics/inc/TGenDecay.h b/math/physics/inc/TGenDecay.h new file mode 100644 index 0000000000000..c69fb235058c0 --- /dev/null +++ b/math/physics/inc/TGenDecay.h @@ -0,0 +1,126 @@ +/************************************************************************* +* Copyright (C) 1995-2021, Rene Brun and Fons Rademakers. * +* All rights reserved. * +* * +* For the licensing terms see $ROOTSYS/LICENSE. * +* For the list of contributors see $ROOTSYS/README/CREDITS. * +*************************************************************************/ + +/*! + @class TGenDecay + +@brief Generates spherical decay and returns the weight of event equals to LIPS (Lorentz Invariant Phase Space). The generator decays a central blob of a given 4-momentum into particles of specific masses. + +### Authors + +Radosław Kycia (kycia.radoslaw@gmail.com), Piotr Lebiedowicz, Antoni Szczurek + +People who helped in the development of the project: Jacek Turnau, Janusz Chwastowski, Rafał Staszewski, Maciej Trzebiński. + + +### Details + +Generator implements TGenInterface interface and works with it is along with this interface. + +The general way of working with the generator is as follows: +1. Prepare 4-momenta of decaying particle P and mass[nt] array of final products. +2. Initialize generator: SetDecay(); +3. Generate decay: Generate(); +4. For each of particles enumerated by 0..nt-1 get their 4-momentum using GetDecay(): pfi = GetDecay(i); +5. Repeat 3 and 4 for another decay. + +Generate() method returns the LIPS: + +\f$ dLIPS = (2\pi)^4 \delta^{(4)}(P-\sum_{i=1}^{n}p_{i}) \prod_{i=1}^{n} \frac{d^{3}p_{i}}{(2\pi)^3 2E_{i}}\f$ + +Comparing to TDecay, it contains a random number generator. + +The generator returns weighted events. + +This class is not recommended for non-advanced users. If you want a robust and adaptive MC generator and integrator, use TGenFoamDecay instead. + + +### Further details + +- R.A. Kycia, P. Lebiedowicz, A. Szczurek, 'Decay: A Monte Carlo library for the decay of a particle with ROOT compatibility', arXiv:2011.14750 [hep-ph] +- R. A. Kycia, J. Chwastowski, R. Staszewski, and J. Turnau, 'GenEx: A simple generator structure for exclusive processes in high energy collisions', Commun. Comput. Phys.24no. 3, (2018) 860, arXiv:1411.6035 [hep-ph] +- R. A. Kycia, J. Turnau, J. J. Chwastowski, R. Staszewski, and M. Trzebinski, 'The adaptiveMonte Carlo toolbox for phase space integration and generation', Commun. Comput. Phys.25no. 5, (2019) 1547, arXiv:1711.06087 [hep-ph] +- R. Hagedorn, 'Relativistic Kinematics: A guide to the kinematic problems of high-energy physics'. W.A. Benjamin, New York, Amsterdam, 1964 +- H.Pilkhun, The Interactions of Hadrons North-Holland 1967 +People who helped in the development of the project: + +*/ + +#ifndef ROOT_TGenDecay +#define ROOT_TGenDecay + +#include +#include + +#include +#include + +#include // std::cin, std::cout +#include // std::queue +#include + + +using namespace std; + + +class TGenDecay : public TObject, public TGenInterface { +private: + Int_t fNt; // number of decay particles + Double_t fMass[18]; // masses of particles + Double_t fBeta[3]; // betas of decaying particle + Double_t fTeCmTm; // total energy in the C.M. minus the total mass + TLorentzVector fDecPro[18]; //kinematics of the generated particles + Int_t seed; //seed for pseudorandom number generator + + Int_t kMAXP; //max number of particles (relict from TGenPhaseSpace) + + TDecay _decay; //decay engine + TRandom3 _pseRan; //pseudorandom numbers + + queue rndQueue; //queue for random numbers + +public: + + /// constructor + TGenDecay(): fNt(0), fMass(), fBeta(), fTeCmTm(0.), seed(4357) {} + /// copy constructor + TGenDecay(const TGenDecay &gen); + /// desctuctor + virtual ~TGenDecay() {} + /// assignment + TGenDecay& operator=(const TGenDecay &gen); + + /// Sets up configuration of decay + /// @param P 4-momentum of decaying particle (Momentum, Energy units are Gev/C, GeV) + /// @param nt number of products of decay + /// @param mass mass matrix of products of decay mass[nt] + /// @returns kTRUE - the decay is permitted by kinematics; kFALSE - the decay is forbidden by kinematics + /// @warning This should be first method to call since it sets up decay configuration. + /// @warning The method also initialize FOAM. + Bool_t SetDecay(TLorentzVector &P, Int_t nt, const Double_t *mass); + + /// Generate a single decay + Double_t Generate(void); + + /// Collect 4-vector of products of decay + /// @param n number of final particle to get from range 1...nt + /// @warning You should call Generate() in first place. + TLorentzVector *GetDecay(Int_t n); + + + /// @returns 4-vector of n-th product of decay + /// @param n number of final particle to get from range 1...nt + Int_t GetNt() const { return fNt;} + + ///sets seed for pseudorandom number generator + void setSeed( Int_t seed ) {seed = seed; _pseRan.SetSeed(seed);}; + +}; + +#endif + diff --git a/math/physics/inc/TGenFoamDecay.h b/math/physics/inc/TGenFoamDecay.h new file mode 100644 index 0000000000000..f5251a4635950 --- /dev/null +++ b/math/physics/inc/TGenFoamDecay.h @@ -0,0 +1,201 @@ +/************************************************************************* +* Copyright (C) 1995-2021, Rene Brun and Fons Rademakers. * +* All rights reserved. * +* * +* For the licensing terms see $ROOTSYS/LICENSE. * +* For the list of contributors see $ROOTSYS/README/CREDITS. * +*************************************************************************/ + +/*! + @class TGenFoamDecay + + + @brief Adaptive MC generator of decay and integrator over LIPS. The generator decays a central blob of a given 4-momentum into particles of specific masses. + +### Authors + +Radosław Kycia (kycia.radoslaw@gmail.com), Piotr Lebiedowicz, Antoni Szczurek + +People who helped in development of the project: Jacek Turnau, Janusz Chwastowski, Rafał Staszewski, Maciej Trzebiński. + +### Details + +This class is a multipurpose MC generator for decay, which also integrates functions over phase space. It uses FOAM adaptive integrator. + +To use this class, make your own generator class that inherits from TGenFoamDecay and redefine Integrand() method, filling it with an integrand function that will be integrated over LIPS (Lorent Invariant Phase Space) and possibly kinematics cuts (return 0.0). If you do not have any specific integrand, then the best choice is to make Integrand() return 1.0. + +The volume element of LIPS is: + +\f$ dLIPS = (2\pi)^4 \delta^{(4)}(P-\sum_{i=1}^{n}p_{i}) \prod_{i=1}^{n} \frac{d^{3}p_{i}}{(2\pi)^3 2E_{i}}\f$ + +Generator implements interface TGenInterface and work with it is along with this interface. + +The general way of working with the generator is as follows: + +1. Create your generator class derived from TGenFoamDecay and implement Integrand() method - return 1.0 if you do not have any specific integrand. +2. Prepare 4-momenta of decaying particle P and mass[nt] array of final products. +3. Initialize generator: SetDecay(); +4. Generate decay: Generate(); +5. For each of particles enumerated by 0..nt-1 get their 4-momentum using GetDecay(): pfi = GetDecay(i); +6. Repeat 3 and 4 for another decay. + +If you prepared the integrand to be integrated over LIPS, you can print the integral final value by calling the Finalize() method or getting the numerical value by GetIntegMC(). + +### Troubleshooting and suggestions + +- By default, the weight of the event is 1.0; however, see OptRej parameter. Sometimes FOAM, even for OptRej=1 weight are not precisely 1.0, so you can remove them or use weights, e.g., in histograms. +- If you do not have an integrand function (e.g., matrix element), then the optimal choice for Integrand() method is to return 1.0 (uniform distribution in the phase space). +- If the integrand has 'small' support in phase space and you get some nonsense result/errors, you probably should increase (say 10x or more) nSampl and nCells parameters of FOAM. It will increase the probability that adaptive MC will have a chance to spot support of the function. However, the exploration phase will last longer. If FOAM still cannot spot the support of the integral, then you can think of using a special Monte Carlo generator for non-spherical decay. +- FOAM has many various setting that can be tested when you deal with demanding integrand functions, see TFoam documentation or the Foam paper: + * S. Jadach, Foam: A General-Purpose Cellular Monte Carlo Event Generator, Comput.Phys.Commun. 152 (2003) 55-100 +- Kinematic cuts can be effectively implemented in Integrand() method - just simply return 0.0 when kinematic cuts are met. +- Foam produces a lot of diagnostic output. It is good to read it at the beginning of implementation. However, in production code or when the generator is embedded in more extensive software, it is better to turn them off. To suppress them set Chat to 0. + + +### Further details + +- R.A. Kycia, P. Lebiedowicz, A. Szczurek, 'Decay: A Monte Carlo library for the decay of a particle with ROOT compatibility', arXiv:2011.14750 [hep-ph] +- R. A. Kycia, J. Chwastowski, R. Staszewski, and J. Turnau, 'GenEx: A simple generator structure for exclusive processes in high energy collisions', Commun. Comput. Phys.24no. 3, (2018) 860, arXiv:1411.6035 [hep-ph] +- R. A. Kycia, J. Turnau, J. J. Chwastowski, R. Staszewski, and M. Trzebinski, 'The adaptiveMonte Carlo toolbox for phase space integration and generation', Commun. Comput. Phys.25no. 5, (2019) 1547, arXiv:1711.06087 [hep-ph] +- R. Hagedorn, 'Relativistic Kinematics: A guide to the kinematic problems of high-energy physics'. W.A. Benjamin, New York, Amsterdam, 1964 +- H.Pilkhun, The Interactions of Hadrons North-Holland 1967 +- S. Jadach, Foam: A General-Purpose Cellular Monte Carlo Event Generator, Comput.Phys.Commun. 152 (2003) 55-100 + +*/ + +#ifndef ROOT_TGenFoamDecay +#define ROOT_TGenFoamDecay + + + +#include + +#include +#include +#include +#include + +#include +#include // std::cin, std::cout +#include // std::queue +#include + +#include + +using namespace std; + + +class TGenFoamDecay : public TFoamIntegrand, public TGenInterface { +private: + Int_t fNt; // number of decay particles + Double_t fMass[18]; // masses of particles + Double_t fBeta[3]; // betas of decaying particle + Double_t fTeCmTm; // total energy in the C.M. minus the total mass + TLorentzVector fDecPro[18]; //kinematics of the generated particles + Int_t seed; //seed for pseudorandom generator + + Int_t kMAXP; //max number of particles (relict from TGenPhaseSpace) + + TDecay _decay; //decay engine + TFoam* _foam; //adaptive integrator + TRandom3 _pseRan; // pseudorandom number generator + + //FOAM parameters + Int_t nCells; // Number of Cells + Int_t nSampl; // Number of MC events per cell in build-up + Int_t nBin; // Number of bins in build-up + Int_t OptRej; // Wted events for OptRej=0; wt=1 for OptRej=1 (default) + Int_t OptDrive; // (D=2) Option, type of Drive =0,1,2 for TrueVol,Sigma,WtMax + Int_t EvPerBin; // Maximum events (equiv.) per bin in buid-up + Int_t Chat; // Chat level + + + + /// @returns weight of the process + /// @param nDim dimension of integration + /// @param Xarg vector of probablilistic variables from [0;1] of dim nDim + /// @warning it is required by Foam integrator + Double_t Density(int nDim, Double_t *Xarg); + + + +public: + + /// constructor + TGenFoamDecay(): fNt(0), fMass(), fBeta(), fTeCmTm(0.), seed(4357), nCells(1000), nSampl(1000), nBin(8), OptRej(1), OptDrive(2), EvPerBin(25), Chat(1) {_foam = new TFoam("FoamX");}; + /// copy constructor + TGenFoamDecay(const TGenFoamDecay &gen); + /// desctructor + virtual ~TGenFoamDecay() { delete _foam;} + /// assignemt + TGenFoamDecay& operator=(const TGenFoamDecay &gen); + + /// Sets up configuration of decay + /// @param P 4-momentum of decaying particle (Momentum, Energy units are Gev/C, GeV) + /// @param nt number of products of decay + /// @param mass mass matrix of products of decay mass[nt] + /// @returns kTRUE - the decay is permitted by kinematics; kFALSE - the decay is forbidden by kinematics + /// @warning This should be first method to call since it sets up decay configuration. + /// @warning The method also initialize FOAM. + Bool_t SetDecay(TLorentzVector &P, Int_t nt, const Double_t *mass); + + /// Generate a single decay + Double_t Generate(void); + + /// Collect 4-vector of products of decay + /// @param n number of final particle to get from range 1...nt + /// @warning You should call Generate() in first place. + TLorentzVector *GetDecay(Int_t n); + + /// @returns 4-vector of n-th product of decay + /// @param n number of final particle to get from range 1...nt + Int_t GetNt() const { return fNt;}; + + /// @returns the function under integral over LIPS (Lorentz Invariant Phase Space) + /// @param fNt number of outgoing particles + /// @param pf array of TLorentzVectors of outgoing particles + /// @warning It is set to 1.0. You should to redefine it (in your derived class, after inheritance) when you want to use full adaptation features. + virtual Double_t Integrand( int fNt, TLorentzVector* pf ); + + ///sets seed for pseudorandom number generator + virtual void setSeed( Int_t seed ) {seed = seed; _pseRan.SetSeed(seed);}; + + ///finalize Foam printing out MC integral and statistics + virtual void Finalize( void ); + + ///@returns integral +- error of MC integration of function in Integrand() over LIPS + virtual void GetIntegMC(Double_t & inetgral, Double_t & error); + + ///Sets FOAM number of Cells + /// @warning It should be done before Generate() + void SetnCells(Int_t nc) {nCells = nc;}; + + ///Sets FOAM number of MC events per cell in build-up + /// @warning It should be done before Generate() + void SetnSampl(Int_t ns) {nSampl = ns;}; + + ///Sets FOAM number of bins in build-up + /// @warning It should be done before Generate() + void SetnBin(Int_t nb) {nBin = nb;}; + + ///Sets FOAM Weigh events for OptRej=0; wt=1 for OptRej=1 (default) + /// @warning It determines if the events will be weighted of not (weight=1). + /// @warning It should be done before Generate() + void SetOptRej(Int_t OptR) {OptRej = OptR;}; + + ///Sets FOAM (D=2) option, type of Drive =0,1,2 for TrueVol,Sigma,WtMax + /// @warning It should be done before Generate() + void SetOptDrive(Int_t OptD) {OptDrive = OptD;}; + + ///Sets FOAM maximum events (equiv.) per bin in buid-up + /// @warning It should be done before Generate() + void SetEvPerBin(Int_t Ev) {EvPerBin = Ev;}; + + ///Sets FOAM chat level + /// @warning It should be done before Generate() + void SetChat(Int_t Ch) {Chat = Ch;}; + +}; + +#endif + diff --git a/math/physics/inc/TGenInterface.h b/math/physics/inc/TGenInterface.h new file mode 100644 index 0000000000000..632f592585d92 --- /dev/null +++ b/math/physics/inc/TGenInterface.h @@ -0,0 +1,94 @@ +/************************************************************************* +* Copyright (C) 1995-2021, Rene Brun and Fons Rademakers. * +* All rights reserved. * +* * +* For the licensing terms see $ROOTSYS/LICENSE. * +* For the list of contributors see $ROOTSYS/README/CREDITS. * +*************************************************************************/ + +/*! + @class TGenInterface + + @brief Minimal interface implemented by MC generators for the decay of a central blob of a given 4-momentum into particles of specific masses. + +### Authors + +Radosław Kycia (kycia.radoslaw@gmail.com), Piotr Lebiedowicz, Antoni Szczurek + +People who helped in development of the project: Jacek Turnau, Janusz Chwastowski, Rafał Staszewski, Maciej Trzebiński. + +### Details +The general concept of working with MC generators is as follows: + +1. Prepare 4-momenta of decaying particle P and mass[nt] array of final products. +2. Initialize generator: SetDecay(); +3. Generate decay: Generate(); +4. For each of particles enumerated by 0..nt-1 get their 4-momentum using GetDecay(): pfi = GetDecay(i); +5. Repeat 3 and 4 for another decay. + +@warning The class is an interface for integration of some distribution over phase space of decaying particle. In the end, you would probably like to normalize distribution properly. + +The interface is implemented in: +- TGenFoamDecay - adaptive Monte Carlo generator and integrator; Probably contains all you need. +- TGenDecay - non-adaptive Monte Carlo generator; returns proper LIPS (Lorentz Invariant Phase Space) weight for an event. +- TDecay - the base class that returns LISP for an event; Can be used in more advanced MC generators as a part that makes decay. +- TGenPhaseSpace - legacy MC generator. + +### Further details + +- R.A. Kycia, P. Lebiedowicz, A. Szczurek, 'Decay: A Monte Carlo library for the decay of a particle with ROOT compatibility', arXiv:2011.14750 [hep-ph] +- R. A. Kycia, J. Chwastowski, R. Staszewski, and J. Turnau, 'GenEx: A simple generator structure for exclusive processes in high energy collisions', Commun. Comput. Phys.24no. 3, (2018) 860, arXiv:1411.6035 [hep-ph] +- R. A. Kycia, J. Turnau, J. J. Chwastowski, R. Staszewski, and M. Trzebinski, 'The adaptiveMonte Carlo toolbox for phase space integration and generation', Commun. Comput. Phys.25no. 5, (2019) 1547, arXiv:1711.06087 [hep-ph] +- R. Hagedorn, 'Relativistic Kinematics: A guide to the kinematic problems of high-energy physics'. W.A. Benjamin, New York, Amsterdam, 1964 +- H.Pilkhun, The Interactions of Hadrons North-Holland 1967 + +*/ + + + +#ifndef ROOT_TGenInterface +#define ROOT_TGenInterface + +#include + + +using namespace std; + + +class TGenInterface { +private: + + Int_t fNt; // number of decay particles + +public: + + ///constructor + TGenInterface(): fNt(0) {}; + + ///destructor + virtual ~TGenInterface() {}; + + /// Sets up configuration of decay + /// @param P 4-momentum of decaying particle + /// @param nt number of products of decay + /// @param mass mass matrix of products of decay mass[nt] + /// @returns kTRUE - the decay is permitted by kinematics; kFALSE - the decay is forbidden by kinematics + /// @warning This should be first method to call since it sets up decay configuration. + virtual Bool_t SetDecay(TLorentzVector &P, Int_t nt, const Double_t *mass)=0; + + /// Generate a single decay + virtual Double_t Generate(void)=0; + + /// @returns 4-vector of n-th product of decay + /// @param n number of final particle to get from range 1...nt + /// @warning You should call Generate() in first place. + virtual TLorentzVector *GetDecay(Int_t n)=0; + + /// @returns number of final particles + virtual Int_t GetNt()const { return fNt;}; + + +}; + +#endif + diff --git a/math/physics/src/TDecay.cxx b/math/physics/src/TDecay.cxx new file mode 100644 index 0000000000000..d3c6eae47d854 --- /dev/null +++ b/math/physics/src/TDecay.cxx @@ -0,0 +1,218 @@ +/************************************************************************* +* Copyright (C) 1995-2021, Rene Brun and Fons Rademakers. * +* All rights reserved. * +* * +* For the licensing terms see $ROOTSYS/LICENSE. * +* For the list of contributors see $ROOTSYS/README/CREDITS. * +*************************************************************************/ + +#include "TDecay.h" + + +//const Int_t kMAXP = 18; + +//_____________________________________________________________________________________ +Double_t TDecay::PDK(Double_t a, Double_t b, Double_t c) +{ + //the PDK function + Double_t x = (a-b-c)*(a+b+c)*(a-b+c)*(a+b-c); + x = TMath::Sqrt(x)/(2*a); + return x; +} + +//_____________________________________________________________________________________ +Int_t DoubleMax(const void *a, const void *b) +{ + //special max function + Double_t aa = * ((Double_t *) a); + Double_t bb = * ((Double_t *) b); + if (aa > bb) return 1; + if (aa < bb) return -1; + return 0; + +} + +//__________________________________________________________________________________________________ +TDecay::TDecay(const TDecay &gen) : TObject(gen) +{ + //copy constructor + fNt = gen.fNt; + fTeCmTm = gen.fTeCmTm; + fBeta[0] = gen.fBeta[0]; + fBeta[1] = gen.fBeta[1]; + fBeta[2] = gen.fBeta[2]; + for (Int_t i=0;i &rnd) +{ + + Double_t rno[kMAXP]; + rno[0] = 0; + Int_t n; + + //check if phase space dim is ok + assert( (unsigned int) (3*fNt - 4) == rnd.size() ); + + double weight = 1.0; // additional weight + + if (fNt>2) + { + for (n=1; n compute the weight of the current event + // + Double_t wt = 1.0; + Double_t pd[kMAXP]; + for (n=0; nPx(); + Double_t z = v->Pz(); + v->SetPz( cY*z - sY*x ); + v->SetPx( sY*z + cY*x); // rotation around Y + x = v->Px(); + Double_t y = v->Py(); + v->SetPx( cZ*x - sZ*y ); + v->SetPy( sZ*x + cZ*y ); // rotation around Z + } + + if (i == (fNt-1)) break; + + Double_t beta = pd[i] / sqrt(pd[i]*pd[i] + invMas[i]*invMas[i]); + for (j=0; j<=i; j++) fDecPro[j].Boost(0,0,beta); + i++; + } + + + // + //---> final boost of all particles + // + for (n=0;n return the weigth of event + // + return wt; +} + +//__________________________________________________________________________________ +TLorentzVector *TDecay::GetDecay(Int_t n) +{ + //return Lorentz vector corresponding to decay n + if (n>fNt) return 0; + return fDecPro+n; +} + +//_____________________________________________________________________________________ +Bool_t TDecay::SetDecay(TLorentzVector &P, Int_t nt, const Double_t *mass) +{ + + kMAXP = nt; + + Int_t n; + fNt = nt; + if (fNt<2 || fNt>18) return kFALSE; // no more then 18 particle + + + // + fTeCmTm = P.Mag(); // total energy in C.M. minus the sum of the masses + for (n=0;n save the betas of the decaying particle + // + if (P.Beta()) { + Double_t w = P.Beta()/P.Rho(); + fBeta[0] = P(0)*w; + fBeta[1] = P(1)*w; + fBeta[2] = P(2)*w; + } + else fBeta[0]=fBeta[1]=fBeta[2]=0; + + return kTRUE; +} diff --git a/math/physics/src/TGenDecay.cxx b/math/physics/src/TGenDecay.cxx new file mode 100644 index 0000000000000..37afea2efb8e3 --- /dev/null +++ b/math/physics/src/TGenDecay.cxx @@ -0,0 +1,109 @@ +/************************************************************************* +* Copyright (C) 1995-2021, Rene Brun and Fons Rademakers. * +* All rights reserved. * +* * +* For the licensing terms see $ROOTSYS/LICENSE. * +* For the list of contributors see $ROOTSYS/README/CREDITS. * +*************************************************************************/ + +#include "TGenDecay.h" + + +//const Int_t kMAXP = 18; + + +//__________________________________________________________________________________________________ +TGenDecay::TGenDecay(const TGenDecay &gen) : TObject(gen) +{ + //copy constructor + fNt = gen.fNt; + fTeCmTm = gen.fTeCmTm; + fBeta[0] = gen.fBeta[0]; + fBeta[1] = gen.fBeta[1]; + fBeta[2] = gen.fBeta[2]; + _decay = gen._decay; + _pseRan = gen._pseRan; + for (Int_t i=0;ifNt) return 0; + + return _decay.GetDecay(n); +} + +//_____________________________________________________________________________________ +Bool_t TGenDecay::SetDecay(TLorentzVector &P, Int_t nt, const Double_t *mass) +{ + + kMAXP = nt; + + Int_t n; + fNt = nt; + if (fNt<2 || fNt>18) return kFALSE; // no more then 18 particle + + // + fTeCmTm = P.Mag(); // total energy in C.M. minus the sum of the masses + for (n=0;n rndQueue; + + //put rnd numbers into queue + for( int i = 0; i < 3*fNt-4; i++) + { + rndQueue.push( Xarg[i] ); + } + + //make decay and take d(LIPS) + double wtdecay = _decay.Generate( rndQueue ); + + //get out particles + TLorentzVector pf[fNt]; + for( int i = 0; i < fNt; i++ ) + { + pf[i] = *(_decay.GetDecay( i )); + } + + //calculate integrand + double integrand = Integrand( fNt, pf ); + + return wtdecay * integrand; +} + +//__________________________________________________________________________________________________ +Double_t TGenFoamDecay::Integrand( int fNt, TLorentzVector * pf ) +{ + return 1.0; //default and probably overloaded for matrix element +} + +//__________________________________________________________________________________________________ +TGenFoamDecay::TGenFoamDecay(const TGenFoamDecay &gen) +{ + //copy constructor + fNt = gen.fNt; + fTeCmTm = gen.fTeCmTm; + fBeta[0] = gen.fBeta[0]; + fBeta[1] = gen.fBeta[1]; + fBeta[2] = gen.fBeta[2]; + _decay = gen._decay; + _foam = gen._foam; + _pseRan = gen._pseRan; + for (Int_t i=0;iMakeEvent(); + + return _foam->GetMCwt( ); + +} + +//__________________________________________________________________________________ +TLorentzVector *TGenFoamDecay::GetDecay(Int_t n) +{ + + if (n>fNt) return 0; + + //return Lorentz vector corresponding to decay of n-th particle + return _decay.GetDecay( n ); +} + +//_____________________________________________________________________________________ +Bool_t TGenFoamDecay::SetDecay(TLorentzVector &P, Int_t nt, + const Double_t *mass) +{ + + kMAXP = nt; + + Int_t n; + fNt = nt; + if (fNt<2 || fNt>18) return kFALSE; // no more then 18 particle + + + fTeCmTm = P.Mag(); // total energy in C.M. minus the sum of the masses + for (n=0;n 0 ) + { + cout<<"***** Foam version "<< _foam->GetVersion() <<" *****"<SetkDim( 3*fNt-4); // Mandatory!!! + _foam->SetnCells( nCells); // optional + _foam->SetnSampl( nSampl); // optional + _foam->SetnBin( nBin); // optional + _foam->SetOptRej( OptRej); // optional + _foam->SetOptDrive( OptDrive); // optional + _foam->SetEvPerBin( EvPerBin); // optional + _foam->SetChat( Chat); // optional + //=============================== + _foam->SetRho(this); + _foam->SetPseRan(&_pseRan); + + // Initialize simulator + _foam->Initialize(); + + return kTRUE; +} + +//_____________________________________________________________________________________ +void TGenFoamDecay::Finalize( void ) +{ + Double_t MCresult,MCerror; + Double_t eps = 0.0005; + Double_t Effic, WtMax, AveWt, Sigma; + Double_t IntNorm, Errel; + _foam->Finalize( IntNorm, Errel); // final printout + _foam->GetIntegMC( MCresult, MCerror); // get MC intnegral + _foam->GetWtParams(eps, AveWt, WtMax, Sigma); // get MC wt parameters + long nCalls=_foam->GetnCalls(); + Effic=0; if(WtMax>0) Effic=AveWt/WtMax; + cout << "================================================================" << endl; + cout << " MCresult= " << MCresult << " +- " << MCerror << " RelErr= "<< MCerror/MCresult << endl; + cout << " Dispersion/= " << Sigma/AveWt << endl; + cout << " /WtMax= " << Effic <<", for epsilon = "< GetIntegMC( integral, error); +} diff --git a/math/physics/test/CMakeLists.txt b/math/physics/test/CMakeLists.txt new file mode 100644 index 0000000000000..53787242e4b21 --- /dev/null +++ b/math/physics/test/CMakeLists.txt @@ -0,0 +1,11 @@ +# Copyright (C) 1995-2019, Rene Brun and Fons Rademakers. +# All rights reserved. +# +# For the licensing terms see $ROOTSYS/LICENSE. +# For the list of contributors see $ROOTSYS/README/CREDITS. + +set(Libraries Physics EG) + +ROOT_ADD_GTEST(testTDecay tTDecay.cxx LIBRARIES ${Libraries}) +ROOT_ADD_GTEST(testTGenDecay tTGenDecay.cxx LIBRARIES ${Libraries}) +ROOT_ADD_GTEST(testTGenFoamDecay tTGenFoamDecay.cxx LIBRARIES ${Libraries}) diff --git a/math/physics/test/tTDecay.cxx b/math/physics/test/tTDecay.cxx new file mode 100644 index 0000000000000..0f4da6a42affe --- /dev/null +++ b/math/physics/test/tTDecay.cxx @@ -0,0 +1,240 @@ +/************************************************************************* + * Copyright (C) 1995-2021, Rene Brun and Fons Rademakers. * + * All rights reserved. * + * * + * For the licensing terms see $ROOTSYS/LICENSE. * + * For the list of contributors see $ROOTSYS/README/CREDITS. * + *************************************************************************/ + +/* +Unit tests for TDecay + +Author: Radosław Kycia (kycia.radoslaw@gmail.com) + +*/ + +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include "TGenDecay.h" + +#include + +using namespace std; + +// decay into 2 particles in final state +TEST(PhaseSpace2Dim, Integrand1) +{ + + // Properites of particles: + TDatabasePDG *PDGDatabese = TDatabasePDG::Instance(); + + // make an object of your generator + queue rndQueue; // queue for random numbers + TDecay generator; // decay engine + TRandom3 pseRan; // pseudorandom numbers + + // set up initial particle/blob that decays: + + // center of mass energy + const double mmu = 0.1057; // muon rest mass GeV + + // CM 4-momentum - blob that initiates decay + TLorentzVector pbCM; + + // set blob CM energy to mass of the muon with zero momentum - rest frame of particle + pbCM.SetPxPyPzE(0.0, 0.0, 0.0, mmu); + + // set up outgoing particles configuration: + + // number of outgoing particles + const int Nop = 2; + + // PDGID outgoing paricles (masses are taken from PDG table) + // full list: https://pdg.lbl.gov/2007/reviews/montecarlorpp.pdf + int idOut[Nop] = {11, 11}; // e, e + + // out particles array for storing their 4-vectors + TLorentzVector pf[Nop]; + + // masses of products - from PDGDatabase + double mass[Nop]; + + // fill in the table with masses of particles using their codes + for (int i = 0; i < Nop; i++) { + mass[i] = PDGDatabese->GetParticle(idOut[i])->Mass(); + } + + // set Generator to decay configuration + generator.SetDecay(pbCM, Nop, mass); + + // phase space element for decay + double wtdecay = 0.0; + + // helper variables for computations of MC integral and its error + double sumIntegrand = 0.0; + double sumIntegrand2 = 0.0; + + long NevTot = 1e6; // Total number of events + + // GENRATION LOOP + for (long loop = 0; loop < NevTot; loop++) { + + // clear queue - in case there was a previous run + while (!rndQueue.empty()) { + rndQueue.pop(); + } + + // number of degrees of freedom in the decay + int nDim = 3 * Nop - 4; + + // put rnd numbers into queue + for (int i = 0; i < nDim; i++) { + rndQueue.push(pseRan.Rndm()); + } + + // make a decay + wtdecay = generator.Generate(rndQueue); + + // get out particles into pf[] array - order as in mass[] array + for (int i = 0; i < Nop; i++) { + pf[i] = *(generator.GetDecay(i)); + } + + // integrand and MC integration + + double integrand = 1.0; + + // update MC sums + sumIntegrand += wtdecay * integrand; + sumIntegrand2 += pow(wtdecay * integrand, 2); + } + // END of GENRATION LOOP + + // caclulate MC integla and its error + double integral = sumIntegrand / (double)NevTot; + // double error = sqrt(abs( pow(integral,2) - sumIntegrand2/(double) NevTot ))/sqrt(NevTot); + + // Theoretical prediction //Theoretical prediction + double th_integral = M_PI * + sqrt((pow(mmu, 2) - pow(mass[0] + mass[1], 2)) * (pow(mmu, 2) - pow(mass[0] - mass[1], 2))) / + (2.0 * pow(mmu, 2) * pow(2.0 * M_PI, 2)); + + ASSERT_NEAR(integral, th_integral, 1e-12); +}; + +// decay into 3 particles in final state - 4-Fermi theory of muon decay +TEST(PhaseSpace3Dim, Fermi4ModelOfDecay) +{ + + // Properites of particles: + TDatabasePDG *PDGDatabese = TDatabasePDG::Instance(); + + // make an object of your generator + queue rndQueue; // queue for random numbers + TDecay generator; // decay engine + TRandom3 pseRan; // pseudorandom numbers + + // set up initial particle/blob that decays: + + // center of mass energy + const double mmu = 0.1057; // muon rest mass GeV + + // CM 4-momentum - blob that initiates decay + TLorentzVector pbCM; + + // set blob CM energy to mass of the muon with zero momentum - rest frame of particle + pbCM.SetPxPyPzE(0.0, 0.0, 0.0, mmu); + + // set up outgoing particles configuration: + + // number of outgoing particles + const int Nop = 3; + + // PDGID outgoing paricles (masses are taken from PDG table) + // full list: https://pdg.lbl.gov/2007/reviews/montecarlorpp.pdf + int idOut[Nop] = {11, 14, 12}; // e, \nu_\mi, \nu_e + + // out particles array for storing their 4-vectors + TLorentzVector pf[Nop]; + + // masses of products - from PDGDatabase + double mass[Nop]; + + // fill in the table with masses of particles using their codes + for (int i = 0; i < Nop; i++) { + mass[i] = PDGDatabese->GetParticle(idOut[i])->Mass(); + } + + // set Generator to decay configuration + generator.SetDecay(pbCM, Nop, mass); + + // phase space element for decay + double wtdecay = 0.0; + + // helper variables for computations of MC integral and its error + double sumIntegrand = 0.0; + double sumIntegrand2 = 0.0; + + long NevTot = 1e6; // Total number of events + + // GENRATION LOOP + for (long loop = 0; loop < NevTot; loop++) { + // Generate event + + // clear queue - in case there was a previous run + while (!rndQueue.empty()) { + rndQueue.pop(); + } + + // number of degrees of freedom in the decay + int nDim = 3 * Nop - 4; + + // put rnd numbers into queue + for (int i = 0; i < nDim; i++) { + rndQueue.push(pseRan.Rndm()); + } + + // make a decay + wtdecay = generator.Generate(rndQueue); + + // get out particles into pf[] array - order as in mass[] array + for (int i = 0; i < Nop; i++) { + pf[i] = *(generator.GetDecay(i)); + } + + // integrand and MC integration + + // Here put your function to integrate over LIPS (Lorentz Invariant Phase Space) + // You can use: pf[0]....pf[Nop-1] - 4 momenta of outgoing particles + const double G = 1.166e-5; + double M2 = 32.0 * G * G * (mmu * mmu - 2.0 * mmu * pf[2].E()) * mmu * pf[2].E(); + double coeff = 1.0 / (2.0 * mmu); + + double integrand = M2 * coeff; + + // update MC sums + sumIntegrand += wtdecay * integrand; + sumIntegrand2 += pow(wtdecay * integrand, 2); + } + // END of GENRATION LOOP + + // caclulate MC integla and its error + double integral = sumIntegrand / (double)NevTot; + double error = sqrt(abs(pow(integral, 2) - sumIntegrand2 / (double)NevTot)) / sqrt(NevTot); + + // Theoretical prediction from (see Problem 5.3 in M.D. Schwatz, 'Quantum Field Theory and The Standard Model', + // Cambridge 2014) + const double G = 1.166e-5; + double th_integral = (G * G * pow(mmu, 5)) / (192.0 * pow(M_PI, 3)); + + ASSERT_TRUE(abs(integral - th_integral) < error); +}; diff --git a/math/physics/test/tTGenDecay.cxx b/math/physics/test/tTGenDecay.cxx new file mode 100644 index 0000000000000..cc7e9792fda29 --- /dev/null +++ b/math/physics/test/tTGenDecay.cxx @@ -0,0 +1,211 @@ +/************************************************************************* + * Copyright (C) 1995-2021, Rene Brun and Fons Rademakers. * + * All rights reserved. * + * * + * For the licensing terms see $ROOTSYS/LICENSE. * + * For the list of contributors see $ROOTSYS/README/CREDITS. * + *************************************************************************/ +/* +Unit tests for TGenDecay + +Author: Radosław Kycia (kycia.radoslaw@gmail.com) +*/ + +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include "TGenDecay.h" + +#include + +using namespace std; + +// decay into 2 particles in final state +TEST(PhaseSpace2Dim, Integrand1) +{ + + // Properites of particles: + TDatabasePDG *PDGDatabese = TDatabasePDG::Instance(); + + // make an object of your generator + TGenDecay generator; + + // set up initial particle/blob that decays: + + // center of mass energy + const double mmu = 0.1057; // muon rest mass GeV + + // CM 4-momentum - blob that initiates decay + TLorentzVector pbCM; + + // set blob CM energy to mass of the muon with zero momentum - rest frame of particle + pbCM.SetPxPyPzE(0.0, 0.0, 0.0, mmu); + + // set up outgoing particles configuration: + + // number of outgoing particles + const int Nop = 2; + + // PDGID outgoing paricles (masses are taken from PDG table) + // full list: https://pdg.lbl.gov/2007/reviews/montecarlorpp.pdf + int idOut[Nop] = {11, 11}; // e, e + + // out particles array for storing their 4-vectors + TLorentzVector pf[Nop]; + + // masses of products - from PDGDatabase + double mass[Nop]; + + // fill in the table with masses of particles using their codes + for (int i = 0; i < Nop; i++) { + mass[i] = PDGDatabese->GetParticle(idOut[i])->Mass(); + } + + // set Generator to decay configuration + generator.SetDecay(pbCM, Nop, mass); + + // phase space element for decay + double wtdecay = 0.0; + + // helper variables for computations of MC integral and its error + double sumIntegrand = 0.0; + double sumIntegrand2 = 0.0; + + long NevTot = 1e6; // Total number of events + + // GENRATION LOOP + for (long loop = 0; loop < NevTot; loop++) { + + // Generate event + + // make decay + wtdecay = generator.Generate(); + + // get out particles into pf[] array - order as in mass[] array + for (int i = 0; i < Nop; i++) { + pf[i] = *(generator.GetDecay(i)); + } + + // integrand and MC integration + + double integrand = 1.0; + + // update MC sums + sumIntegrand += wtdecay * integrand; + sumIntegrand2 += pow(wtdecay * integrand, 2); + } + // END of GENRATION LOOP + + // caclulate MC integla and its error + double integral = sumIntegrand / (double)NevTot; + // double error = sqrt(abs( pow(integral,2) - sumIntegrand2/(double) NevTot ))/sqrt(NevTot); + + // Theoretical prediction //Theoretical prediction + double th_integral = M_PI * + sqrt((pow(mmu, 2) - pow(mass[0] + mass[1], 2)) * (pow(mmu, 2) - pow(mass[0] - mass[1], 2))) / + (2.0 * pow(mmu, 2) * pow(2.0 * M_PI, 2)); + + ASSERT_NEAR(integral, th_integral, 1e-12); +}; + +// decay into 3 particles in final state - 4-Fermi theory of muon decay +TEST(PhaseSpace3Dim, Fermi4ModelOfDecay) +{ + + // Properites of particles: + TDatabasePDG *PDGDatabese = TDatabasePDG::Instance(); + + // make an object of your generator + TGenDecay generator; + + // set up initial particle/blob that decays: + + // center of mass energy + const double mmu = 0.1057; // muon rest mass GeV + + // CM 4-momentum - blob that initiates decay + TLorentzVector pbCM; + + // set blob CM energy to mass of the muon with zero momentum - rest frame of particle + pbCM.SetPxPyPzE(0.0, 0.0, 0.0, mmu); + + // set up outgoing particles configuration: + + // number of outgoing particles + const int Nop = 3; + + // PDGID outgoing paricles (masses are taken from PDG table) + // full list: https://pdg.lbl.gov/2007/reviews/montecarlorpp.pdf + int idOut[Nop] = {11, 14, 12}; // e, \nu_\mi, \nu_e + + // out particles array for storing their 4-vectors + TLorentzVector pf[Nop]; + + // masses of products - from PDGDatabase + double mass[Nop]; + + // fill in the table with masses of particles using their codes + for (int i = 0; i < Nop; i++) { + mass[i] = PDGDatabese->GetParticle(idOut[i])->Mass(); + } + + // set Generator to decay configuration + generator.SetDecay(pbCM, Nop, mass); + + // phase space element for decay + double wtdecay = 0.0; + + // helper variables for computations of MC integral and its error + double sumIntegrand = 0.0; + double sumIntegrand2 = 0.0; + + long NevTot = 1e6; // Total number of events + + // GENRATION LOOP + for (long loop = 0; loop < NevTot; loop++) { + + // Generate event + + // make decay + wtdecay = generator.Generate(); + + // get out particles into pf[] array - order as in mass[] array + for (int i = 0; i < Nop; i++) { + pf[i] = *(generator.GetDecay(i)); + } + + // integrand and MC integration + + // Here put your function to integrate over LIPS (Lorentz Invariant Phase Space) + // You can use: pf[0]....pf[Nop-1] - 4 momenta of outgoing particles + const double G = 1.166e-5; + double M2 = 32.0 * G * G * (mmu * mmu - 2.0 * mmu * pf[2].E()) * mmu * pf[2].E(); + double coeff = 1.0 / (2.0 * mmu); + + double integrand = M2 * coeff; + + // update MC sums + sumIntegrand += wtdecay * integrand; + sumIntegrand2 += pow(wtdecay * integrand, 2); + } + // END of GENRATION LOOP + + // caclulate MC integla and its error + double integral = sumIntegrand / (double)NevTot; + double error = sqrt(abs(pow(integral, 2) - sumIntegrand2 / (double)NevTot)) / sqrt(NevTot); + + // Theoretical prediction from (see Problem 5.3 in M.D. Schwatz, 'Quantum Field Theory and The Standard Model', + // Cambridge 2014) + const double G = 1.166e-5; + double th_integral = (G * G * pow(mmu, 5)) / (192.0 * pow(M_PI, 3)); + + ASSERT_TRUE(abs(integral - th_integral) < error); +}; diff --git a/math/physics/test/tTGenFoamDecay.cxx b/math/physics/test/tTGenFoamDecay.cxx new file mode 100644 index 0000000000000..ab8b0b355c4f6 --- /dev/null +++ b/math/physics/test/tTGenFoamDecay.cxx @@ -0,0 +1,211 @@ +/************************************************************************* + * Copyright (C) 1995-2021, Rene Brun and Fons Rademakers. * + * All rights reserved. * + * * + * For the licensing terms see $ROOTSYS/LICENSE. * + * For the list of contributors see $ROOTSYS/README/CREDITS. * + *************************************************************************/ +/* +Unit tests for TGenFoamDecay +Author: Radosław Kycia (kycia.radoslaw@gmail.com) +*/ + +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include "TGenFoamDecay.h" + +#include + +using namespace std; + +// decay into 2 particles in final state +TEST(PhaseSpace2Dim, Integrand1) +{ + + // Properites of particles: + TDatabasePDG *PDGDatabese = TDatabasePDG::Instance(); + + // Create your own generator class derived from TGenFoamDecay + class Generator : public TGenFoamDecay { + Double_t Integrand(int /*nt*/, TLorentzVector * /*pf*/) + { + + // Here put your function to integrate over LIPS (Lorentz Invariant Phase Space) + // You can use: pf[0]....pf[Nop-1] - 4 momenta of outgoing particles + return 1.0; + }; + }; + + // make an object of your generator + Generator generator; + + // silent mode + generator.SetChat(0); + + // set up initial particle/blob that decays: + + // center of mass energy + const double mmu = 0.1057; // muon rest mass GeV + + // CM 4-momentum - blob that initiates decay + TLorentzVector pbCM; + + // set blob CM energy to mass of the muon with zero momentum - rest frame of particle + pbCM.SetPxPyPzE(0.0, 0.0, 0.0, mmu); + + // set up outgoing particles configuration: + + // number of outgoing particles + const int Nop = 2; + + // PDGID outgoing paricles (masses are taken from PDG table) + // full list: https://pdg.lbl.gov/2007/reviews/montecarlorpp.pdf + int idOut[Nop] = {11, 11}; // e, e + + // out particles array for storing their 4-vectors + TLorentzVector pf[Nop]; + + // masses of products - from PDGDatabase + double mass[Nop]; + + // fill in the table with masses of particles using their codes + for (int i = 0; i < Nop; i++) { + mass[i] = PDGDatabese->GetParticle(idOut[i])->Mass(); + } + + // set Generator to decay configuration + generator.SetDecay(pbCM, Nop, mass); + + long NevTot = 10e4; // Total number of events to generate + + // GENRATION LOOP + for (long loop = 0; loop < NevTot; loop++) { + + // Generate event + + // make decay + generator.Generate(); + + // get out particles into pf[] array - order as in mass[] array + for (int i = 0; i < Nop; i++) { + pf[i] = *(generator.GetDecay(i)); + } + } + // END of GENRATION LOOP + + // generator.Finalize(); + + // get integral and error from Generator + Double_t integral, error; + generator.GetIntegMC(integral, error); + + // Theoretical prediction //Theoretical prediction + double th_integral = M_PI * + sqrt((pow(mmu, 2) - pow(mass[0] + mass[1], 2)) * (pow(mmu, 2) - pow(mass[0] - mass[1], 2))) / + (2.0 * pow(mmu, 2) * pow(2.0 * M_PI, 2)); + + ASSERT_NEAR(integral, th_integral, 1e-15); +}; + +// decay into 3 particles in final state - 4-Fermi theory of muon decay +TEST(PhaseSpace3Dim, Fermi4ModelOfDecay) +{ + + // Properites of particles: + TDatabasePDG *PDGDatabese = TDatabasePDG::Instance(); + + // Create your own generator class derived from TGenFoamDecay + class Generator : public TGenFoamDecay { + Double_t Integrand(int /*nt*/, TLorentzVector *pf) + { + + // Here put your function to integrate over LIPS (Lorentz Invariant Phase Space) + // You can use: pf[0]....pf[Nop-1] - 4 momenta of outgoing particles + const double G = 1.166e-5; + const double mmu = 0.1057; + double M2 = 32.0 * G * G * (mmu * mmu - 2.0 * mmu * pf[2].E()) * mmu * pf[2].E(); + double coeff = 1.0 / (2.0 * mmu); + + return M2 * coeff; + }; + }; + + // make an object of your generator + Generator generator; + + // silent mode + generator.SetChat(0); + + // set up initial particle/blob that decays: + + // center of mass energy + const double mmu = 0.1057; // muon rest mass GeV + + // CM 4-momentum - blob that initiates decay + TLorentzVector pbCM; + + // set blob CM energy to mass of the muon with zero momentum - rest frame of particle + pbCM.SetPxPyPzE(0.0, 0.0, 0.0, mmu); + + // set up outgoing particles configuration: + + // number of outgoing particles + const int Nop = 3; + + // PDGID outgoing paricles (masses are taken from PDG table) + // full list: https://pdg.lbl.gov/2007/reviews/montecarlorpp.pdf + int idOut[Nop] = {11, 14, 12}; // e, \nu_\mi, \nu_e + + // out particles array for storing their 4-vectors + TLorentzVector pf[Nop]; + + // masses of products - from PDGDatabase + double mass[Nop]; + + // fill in the table with masses of particles using their codes + for (int i = 0; i < Nop; i++) { + mass[i] = PDGDatabese->GetParticle(idOut[i])->Mass(); + } + + // set Generator to decay configuration + generator.SetDecay(pbCM, Nop, mass); + + long NevTot = 10e4; // Total number of events to generate + + // GENRATION LOOP + for (long loop = 0; loop < NevTot; loop++) { + + // Generate event + + // make decay + generator.Generate(); + + // get out particles into pf[] array - order as in mass[] array + for (int i = 0; i < Nop; i++) { + pf[i] = *(generator.GetDecay(i)); + } + } + // END of GENRATION LOOP + + // generator.Finalize(); + + // get integral and error from Generator + Double_t integral, error; + generator.GetIntegMC(integral, error); + + // Theoretical prediction from (see Problem 5.3 in M.D. Schwatz, 'Quantum Field Theory and The Standard Model', + // Cambridge 2014) + const double G = 1.166e-5; + double th_integral = (G * G * pow(mmu, 5)) / (192.0 * pow(M_PI, 3)); + + ASSERT_TRUE(abs(integral - th_integral) < error); +}; From 3541197bb6345f2823577e5b4cdeae4a6c2103d5 Mon Sep 17 00:00:00 2001 From: Axel Naumann Date: Thu, 4 Mar 2021 18:14:23 +0100 Subject: [PATCH 2/3] [physics] Hide static / local DoubleMax symbol. --- math/physics/src/TGenPhaseSpace.cxx | 2 ++ 1 file changed, 2 insertions(+) diff --git a/math/physics/src/TGenPhaseSpace.cxx b/math/physics/src/TGenPhaseSpace.cxx index 6f62862d667cc..07caa6c4ac7b4 100644 --- a/math/physics/src/TGenPhaseSpace.cxx +++ b/math/physics/src/TGenPhaseSpace.cxx @@ -38,6 +38,7 @@ Double_t TGenPhaseSpace::PDK(Double_t a, Double_t b, Double_t c) return x; } +namespace { //////////////////////////////////////////////////////////////////////////////// /// Special max function @@ -50,6 +51,7 @@ Int_t DoubleMax(const void *a, const void *b) return 0; } +} //////////////////////////////////////////////////////////////////////////////// /// Copy constructor From 37c205e1fd6d1e34aeed5b9b3a758df4fe23a033 Mon Sep 17 00:00:00 2001 From: Axel Naumann Date: Thu, 4 Mar 2021 18:10:15 +0100 Subject: [PATCH 3/3] [physics] Coding style for T{Gen{Foam}}Decay. --- math/physics/inc/TDecay.h | 99 +++++------ math/physics/inc/TGenDecay.h | 101 +++++------ math/physics/inc/TGenFoamDecay.h | 264 ++++++++++++++++------------- math/physics/src/TDecay.cxx | 235 +++++++++++++------------ math/physics/src/TGenDecay.cxx | 110 ++++++------ math/physics/src/TGenFoamDecay.cxx | 226 ++++++++++++------------ 6 files changed, 525 insertions(+), 510 deletions(-) diff --git a/math/physics/inc/TDecay.h b/math/physics/inc/TDecay.h index 4e2ce55e14d1e..9cb7f6d86c497 100644 --- a/math/physics/inc/TDecay.h +++ b/math/physics/inc/TDecay.h @@ -1,35 +1,39 @@ /************************************************************************* -* Copyright (C) 1995-2021, Rene Brun and Fons Rademakers. * -* All rights reserved. * -* * -* For the licensing terms see $ROOTSYS/LICENSE. * -* For the list of contributors see $ROOTSYS/README/CREDITS. * -*************************************************************************/ + * Copyright (C) 1995-2021, Rene Brun and Fons Rademakers. * + * All rights reserved. * + * * + * For the licensing terms see $ROOTSYS/LICENSE. * + * For the list of contributors see $ROOTSYS/README/CREDITS. * + *************************************************************************/ /*! @class TDecay - @brief The generator decays a central blob of a given 4-momentum into particles of specific masses. Original GENBOD algorithm with w few modifications. - + @brief The generator decays a central blob of a given 4-momentum into particles of specific masses. Original GENBOD +algorithm with w few modifications. + ### Authors - + Radosław Kycia (kycia.radoslaw@gmail.com), Piotr Lebiedowicz, Antoni Szczurek -People who helped in the development of the project: Jacek Turnau, Janusz Chwastowski, Rafał Staszewski, Maciej Trzebiński. +People who helped in the development of the project: Jacek Turnau, Janusz Chwastowski, Rafał Staszewski, Maciej +Trzebiński. ### Description -This is a helper class that can be base to construct MC generators that simulate decay. It returns weight equals LIPS (Lorentz Invariant Phase Space): - +This is a helper class that can be base to construct MC generators that simulate decay. It returns weight equals LIPS +(Lorentz Invariant Phase Space): + \f$ dLIPS = (2\pi)^4 \delta^{(4)}(P-\sum_{i=1}^{n}p_{i}) \prod_{i=1}^{n} \frac{d^{3}p_{i}}{(2\pi)^3 2E_{i}}\f$ -The generator can be used as a building block in more complicated decays. +The generator can be used as a building block in more complicated decays. The generator returns weighted events. The class is adapted from the TGenPhaseSpace ROOT package. -This class is not recommended for non-advanced users. If you want a robust and adaptive MC generator and integrator, use TGenFoamDecay instead. +This class is not recommended for non-advanced users. If you want a robust and adaptive MC generator and integrator, use +TGenFoamDecay instead. #### The scheme of use @@ -39,58 +43,59 @@ This class is not recommended for non-advanced users. If you want a robust and a 4. For each of particles enumerated by 0..nt-1 get their 4-momentum using GetDecay(): pfi = GetDecay(i); 5. Repeat 3 and 4 for another decay. - + ### Further details -- R.A. Kycia, P. Lebiedowicz, A. Szczurek, 'Decay: A Monte Carlo library for the decay of a particle with ROOT compatibility', arXiv:2011.14750 [hep-ph] -- R. A. Kycia, J. Chwastowski, R. Staszewski, and J. Turnau, 'GenEx: A simple generator structure for exclusive processes in high energy collisions', Commun. Comput. Phys.24no. 3, (2018) 860, arXiv:1411.6035 [hep-ph] -- R. A. Kycia, J. Turnau, J. J. Chwastowski, R. Staszewski, and M. Trzebinski, 'The adaptiveMonte Carlo toolbox for phase space integration and generation', Commun. Comput. Phys.25no. 5, (2019) 1547, arXiv:1711.06087 [hep-ph] -- R. Hagedorn, 'Relativistic Kinematics: A guide to the kinematic problems of high-energy physics'. W.A. Benjamin, New York, Amsterdam, 1964 -- H.Pilkhun, The Interactions of Hadrons North-Holland 1967 - -*/ +- R.A. Kycia, P. Lebiedowicz, A. Szczurek, 'Decay: A Monte Carlo library for the decay of a particle with ROOT +compatibility', arXiv:2011.14750 [hep-ph] +- R. A. Kycia, J. Chwastowski, R. Staszewski, and J. Turnau, 'GenEx: A simple generator structure for exclusive +processes in high energy collisions', Commun. Comput. Phys.24no. 3, (2018) 860, arXiv:1411.6035 [hep-ph] +- R. A. Kycia, J. Turnau, J. J. Chwastowski, R. Staszewski, and M. Trzebinski, 'The adaptiveMonte Carlo toolbox for +phase space integration and generation', Commun. Comput. Phys.25no. 5, (2019) 1547, arXiv:1711.06087 [hep-ph] +- R. Hagedorn, 'Relativistic Kinematics: A guide to the kinematic problems of high-energy physics'. W.A. Benjamin, New +York, Amsterdam, 1964 +- H.Pilkhun, The Interactions of Hadrons North-Holland 1967 +*/ #ifndef ROOT_TDecay #define ROOT_TDecay #include -#include // std::cin, std::cout -#include // std::queue +#include // std::cin, std::cout +#include // std::queue #include using namespace std; - class TDecay : public TObject { -private: - Int_t fNt; // number of decay particles - Double_t fMass[18]; // masses of particles - Double_t fBeta[3]; // betas of decaying particle - Double_t fTeCmTm; // total energy in the C.M. minus the total mass - TLorentzVector fDecPro[18]; //kinematics of the generated particles +private: + Int_t fNt; // number of decay particles + Double_t fMass[18]; // masses of particles + Double_t fBeta[3]; // betas of decaying particle + Double_t fTeCmTm; // total energy in the C.M. minus the total mass + TLorentzVector fDecPro[18]; // kinematics of the generated particles - Int_t kMAXP; //max number of particles (relict from TGenPhaseSpace) + Int_t kMAXP; // max number of particles (relict from TGenPhaseSpace) Double_t PDK(Double_t a, Double_t b, Double_t c); - - ///factorial function - int factorial(int n); -public: + /// factorial function + int factorial(int n); +public: /// constructor - TDecay(): fNt(0), fMass(), fBeta(), fTeCmTm(0.) {}; - + TDecay() : fNt(0), fMass(), fBeta(), fTeCmTm(0.){}; + /// copy constructor TDecay(const TDecay &gen); - + /// desctructor - virtual ~TDecay() {}; - + virtual ~TDecay(){}; + /// assignment - TDecay& operator=(const TDecay &gen); + TDecay &operator=(const TDecay &gen); /// Sets up configuration of decay /// @param P 4-momentum of decaying particle @@ -99,20 +104,18 @@ class TDecay : public TObject { /// @returns kTRUE - the decay is permitted by kinematics; kFALSE - the decay is forbidden by kinematics /// @warning This should be first method to call since it sets up decay configuration. Bool_t SetDecay(TLorentzVector &P, Int_t nt, const Double_t *mass); - + /// Generate a single decay /// @param rnd a queue of 3*nt-4 random numbers from external source. Double_t Generate(std::queue &rnd); - + /// @returns 4-vector of n-th product of decay /// @param n number of final particle to get from range 1...nt /// @warning You should call Generate() in first place. - TLorentzVector *GetDecay(Int_t n); + TLorentzVector *GetDecay(Int_t n); /// @returns number of final particles - Int_t GetNt() const { return fNt;} - + Int_t GetNt() const { return fNt; } }; #endif - diff --git a/math/physics/inc/TGenDecay.h b/math/physics/inc/TGenDecay.h index c69fb235058c0..cf349754b1c66 100644 --- a/math/physics/inc/TGenDecay.h +++ b/math/physics/inc/TGenDecay.h @@ -1,26 +1,28 @@ /************************************************************************* -* Copyright (C) 1995-2021, Rene Brun and Fons Rademakers. * -* All rights reserved. * -* * -* For the licensing terms see $ROOTSYS/LICENSE. * -* For the list of contributors see $ROOTSYS/README/CREDITS. * -*************************************************************************/ + * Copyright (C) 1995-2021, Rene Brun and Fons Rademakers. * + * All rights reserved. * + * * + * For the licensing terms see $ROOTSYS/LICENSE. * + * For the list of contributors see $ROOTSYS/README/CREDITS. * + *************************************************************************/ /*! @class TGenDecay -@brief Generates spherical decay and returns the weight of event equals to LIPS (Lorentz Invariant Phase Space). The generator decays a central blob of a given 4-momentum into particles of specific masses. +@brief Generates spherical decay and returns the weight of event equals to LIPS (Lorentz Invariant Phase Space). The +generator decays a central blob of a given 4-momentum into particles of specific masses. ### Authors - + Radosław Kycia (kycia.radoslaw@gmail.com), Piotr Lebiedowicz, Antoni Szczurek -People who helped in the development of the project: Jacek Turnau, Janusz Chwastowski, Rafał Staszewski, Maciej Trzebiński. +People who helped in the development of the project: Jacek Turnau, Janusz Chwastowski, Rafał Staszewski, Maciej +Trzebiński. ### Details -Generator implements TGenInterface interface and works with it is along with this interface. +Generator implements TGenInterface interface and works with it is along with this interface. The general way of working with the generator is as follows: 1. Prepare 4-momenta of decaying particle P and mass[nt] array of final products. @@ -37,18 +39,23 @@ Comparing to TDecay, it contains a random number generator. The generator returns weighted events. -This class is not recommended for non-advanced users. If you want a robust and adaptive MC generator and integrator, use TGenFoamDecay instead. +This class is not recommended for non-advanced users. If you want a robust and adaptive MC generator and integrator, use +TGenFoamDecay instead. + - ### Further details -- R.A. Kycia, P. Lebiedowicz, A. Szczurek, 'Decay: A Monte Carlo library for the decay of a particle with ROOT compatibility', arXiv:2011.14750 [hep-ph] -- R. A. Kycia, J. Chwastowski, R. Staszewski, and J. Turnau, 'GenEx: A simple generator structure for exclusive processes in high energy collisions', Commun. Comput. Phys.24no. 3, (2018) 860, arXiv:1411.6035 [hep-ph] -- R. A. Kycia, J. Turnau, J. J. Chwastowski, R. Staszewski, and M. Trzebinski, 'The adaptiveMonte Carlo toolbox for phase space integration and generation', Commun. Comput. Phys.25no. 5, (2019) 1547, arXiv:1711.06087 [hep-ph] -- R. Hagedorn, 'Relativistic Kinematics: A guide to the kinematic problems of high-energy physics'. W.A. Benjamin, New York, Amsterdam, 1964 -- H.Pilkhun, The Interactions of Hadrons North-Holland 1967 +- R.A. Kycia, P. Lebiedowicz, A. Szczurek, 'Decay: A Monte Carlo library for the decay of a particle with ROOT +compatibility', arXiv:2011.14750 [hep-ph] +- R. A. Kycia, J. Chwastowski, R. Staszewski, and J. Turnau, 'GenEx: A simple generator structure for exclusive +processes in high energy collisions', Commun. Comput. Phys.24no. 3, (2018) 860, arXiv:1411.6035 [hep-ph] +- R. A. Kycia, J. Turnau, J. J. Chwastowski, R. Staszewski, and M. Trzebinski, 'The adaptiveMonte Carlo toolbox for +phase space integration and generation', Commun. Comput. Phys.25no. 5, (2019) 1547, arXiv:1711.06087 [hep-ph] +- R. Hagedorn, 'Relativistic Kinematics: A guide to the kinematic problems of high-energy physics'. W.A. Benjamin, New +York, Amsterdam, 1964 +- H.Pilkhun, The Interactions of Hadrons North-Holland 1967 People who helped in the development of the project: - + */ #ifndef ROOT_TGenDecay @@ -60,67 +67,65 @@ People who helped in the development of the project: #include #include -#include // std::cin, std::cout -#include // std::queue +#include // std::cin, std::cout +#include // std::queue #include - using namespace std; - class TGenDecay : public TObject, public TGenInterface { -private: - Int_t fNt; // number of decay particles - Double_t fMass[18]; // masses of particles - Double_t fBeta[3]; // betas of decaying particle - Double_t fTeCmTm; // total energy in the C.M. minus the total mass - TLorentzVector fDecPro[18]; //kinematics of the generated particles - Int_t seed; //seed for pseudorandom number generator +private: + Int_t fNt; // number of decay particles + Double_t fMass[18]; // masses of particles + Double_t fBeta[3]; // betas of decaying particle + Double_t fTeCmTm; // total energy in the C.M. minus the total mass + TLorentzVector fDecPro[18]; // kinematics of the generated particles + Int_t fSeed; // seed for pseudorandom number generator - Int_t kMAXP; //max number of particles (relict from TGenPhaseSpace) + Int_t kMAXP; // max number of particles (relict from TGenPhaseSpace) - TDecay _decay; //decay engine - TRandom3 _pseRan; //pseudorandom numbers - - queue rndQueue; //queue for random numbers + TDecay fDecay; // decay engine + TRandom3 fPseRan; // pseudorandom numbers -public: + queue fRndQueue; // queue for random numbers +public: /// constructor - TGenDecay(): fNt(0), fMass(), fBeta(), fTeCmTm(0.), seed(4357) {} + TGenDecay() : fNt(0), fMass(), fBeta(), fTeCmTm(0.), fSeed(4357) {} /// copy constructor TGenDecay(const TGenDecay &gen); /// desctuctor virtual ~TGenDecay() {} /// assignment - TGenDecay& operator=(const TGenDecay &gen); + TGenDecay &operator=(const TGenDecay &gen); /// Sets up configuration of decay /// @param P 4-momentum of decaying particle (Momentum, Energy units are Gev/C, GeV) /// @param nt number of products of decay /// @param mass mass matrix of products of decay mass[nt] /// @returns kTRUE - the decay is permitted by kinematics; kFALSE - the decay is forbidden by kinematics - /// @warning This should be first method to call since it sets up decay configuration. + /// @warning This should be first method to call since it sets up decay configuration. /// @warning The method also initialize FOAM. Bool_t SetDecay(TLorentzVector &P, Int_t nt, const Double_t *mass); - + /// Generate a single decay Double_t Generate(void); - + /// Collect 4-vector of products of decay /// @param n number of final particle to get from range 1...nt /// @warning You should call Generate() in first place. - TLorentzVector *GetDecay(Int_t n); - + TLorentzVector *GetDecay(Int_t n); /// @returns 4-vector of n-th product of decay /// @param n number of final particle to get from range 1...nt - Int_t GetNt() const { return fNt;} - - ///sets seed for pseudorandom number generator - void setSeed( Int_t seed ) {seed = seed; _pseRan.SetSeed(seed);}; - + Int_t GetNt() const { return fNt; } + + /// sets seed for pseudorandom number generator + void setSeed(Int_t seed) + { + fSeed = seed; + fPseRan.SetSeed(seed); + }; }; #endif - diff --git a/math/physics/inc/TGenFoamDecay.h b/math/physics/inc/TGenFoamDecay.h index f5251a4635950..24cf86f34c453 100644 --- a/math/physics/inc/TGenFoamDecay.h +++ b/math/physics/inc/TGenFoamDecay.h @@ -1,64 +1,84 @@ /************************************************************************* -* Copyright (C) 1995-2021, Rene Brun and Fons Rademakers. * -* All rights reserved. * -* * -* For the licensing terms see $ROOTSYS/LICENSE. * -* For the list of contributors see $ROOTSYS/README/CREDITS. * -*************************************************************************/ + * Copyright (C) 1995-2021, Rene Brun and Fons Rademakers. * + * All rights reserved. * + * * + * For the licensing terms see $ROOTSYS/LICENSE. * + * For the list of contributors see $ROOTSYS/README/CREDITS. * + *************************************************************************/ /*! @class TGenFoamDecay - @brief Adaptive MC generator of decay and integrator over LIPS. The generator decays a central blob of a given 4-momentum into particles of specific masses. - + @brief Adaptive MC generator of decay and integrator over LIPS. The generator decays a central blob of a given +4-momentum into particles of specific masses. + ### Authors - + Radosław Kycia (kycia.radoslaw@gmail.com), Piotr Lebiedowicz, Antoni Szczurek People who helped in development of the project: Jacek Turnau, Janusz Chwastowski, Rafał Staszewski, Maciej Trzebiński. ### Details -This class is a multipurpose MC generator for decay, which also integrates functions over phase space. It uses FOAM adaptive integrator. +This class is a multipurpose MC generator for decay, which also integrates functions over phase space. It uses FOAM +adaptive integrator. -To use this class, make your own generator class that inherits from TGenFoamDecay and redefine Integrand() method, filling it with an integrand function that will be integrated over LIPS (Lorent Invariant Phase Space) and possibly kinematics cuts (return 0.0). If you do not have any specific integrand, then the best choice is to make Integrand() return 1.0. +To use this class, make your own generator class that inherits from TGenFoamDecay and redefine Integrand() method, +filling it with an integrand function that will be integrated over LIPS (Lorent Invariant Phase Space) and possibly +kinematics cuts (return 0.0). If you do not have any specific integrand, then the best choice is to make Integrand() +return 1.0. The volume element of LIPS is: - + \f$ dLIPS = (2\pi)^4 \delta^{(4)}(P-\sum_{i=1}^{n}p_{i}) \prod_{i=1}^{n} \frac{d^{3}p_{i}}{(2\pi)^3 2E_{i}}\f$ - + Generator implements interface TGenInterface and work with it is along with this interface. - + The general way of working with the generator is as follows: - -1. Create your generator class derived from TGenFoamDecay and implement Integrand() method - return 1.0 if you do not have any specific integrand. + +1. Create your generator class derived from TGenFoamDecay and implement Integrand() method - return 1.0 if you do not +have any specific integrand. 2. Prepare 4-momenta of decaying particle P and mass[nt] array of final products. 3. Initialize generator: SetDecay(); 4. Generate decay: Generate(); 5. For each of particles enumerated by 0..nt-1 get their 4-momentum using GetDecay(): pfi = GetDecay(i); 6. Repeat 3 and 4 for another decay. - -If you prepared the integrand to be integrated over LIPS, you can print the integral final value by calling the Finalize() method or getting the numerical value by GetIntegMC(). + +If you prepared the integrand to be integrated over LIPS, you can print the integral final value by calling the +Finalize() method or getting the numerical value by GetIntegMC(). ### Troubleshooting and suggestions -- By default, the weight of the event is 1.0; however, see OptRej parameter. Sometimes FOAM, even for OptRej=1 weight are not precisely 1.0, so you can remove them or use weights, e.g., in histograms. -- If you do not have an integrand function (e.g., matrix element), then the optimal choice for Integrand() method is to return 1.0 (uniform distribution in the phase space). -- If the integrand has 'small' support in phase space and you get some nonsense result/errors, you probably should increase (say 10x or more) nSampl and nCells parameters of FOAM. It will increase the probability that adaptive MC will have a chance to spot support of the function. However, the exploration phase will last longer. If FOAM still cannot spot the support of the integral, then you can think of using a special Monte Carlo generator for non-spherical decay. -- FOAM has many various setting that can be tested when you deal with demanding integrand functions, see TFoam documentation or the Foam paper: +- By default, the weight of the event is 1.0; however, see OptRej parameter. Sometimes FOAM, even for OptRej=1 weight +are not precisely 1.0, so you can remove them or use weights, e.g., in histograms. +- If you do not have an integrand function (e.g., matrix element), then the optimal choice for Integrand() method is to +return 1.0 (uniform distribution in the phase space). +- If the integrand has 'small' support in phase space and you get some nonsense result/errors, you probably should +increase (say 10x or more) nSampl and nCells parameters of FOAM. It will increase the probability that adaptive MC will +have a chance to spot support of the function. However, the exploration phase will last longer. If FOAM still cannot +spot the support of the integral, then you can think of using a special Monte Carlo generator for non-spherical decay. +- FOAM has many various setting that can be tested when you deal with demanding integrand functions, see TFoam +documentation or the Foam paper: * S. Jadach, Foam: A General-Purpose Cellular Monte Carlo Event Generator, Comput.Phys.Commun. 152 (2003) 55-100 -- Kinematic cuts can be effectively implemented in Integrand() method - just simply return 0.0 when kinematic cuts are met. -- Foam produces a lot of diagnostic output. It is good to read it at the beginning of implementation. However, in production code or when the generator is embedded in more extensive software, it is better to turn them off. To suppress them set Chat to 0. +- Kinematic cuts can be effectively implemented in Integrand() method - just simply return 0.0 when kinematic cuts are +met. +- Foam produces a lot of diagnostic output. It is good to read it at the beginning of implementation. However, in +production code or when the generator is embedded in more extensive software, it is better to turn them off. To suppress +them set Chat to 0. ### Further details -- R.A. Kycia, P. Lebiedowicz, A. Szczurek, 'Decay: A Monte Carlo library for the decay of a particle with ROOT compatibility', arXiv:2011.14750 [hep-ph] -- R. A. Kycia, J. Chwastowski, R. Staszewski, and J. Turnau, 'GenEx: A simple generator structure for exclusive processes in high energy collisions', Commun. Comput. Phys.24no. 3, (2018) 860, arXiv:1411.6035 [hep-ph] -- R. A. Kycia, J. Turnau, J. J. Chwastowski, R. Staszewski, and M. Trzebinski, 'The adaptiveMonte Carlo toolbox for phase space integration and generation', Commun. Comput. Phys.25no. 5, (2019) 1547, arXiv:1711.06087 [hep-ph] -- R. Hagedorn, 'Relativistic Kinematics: A guide to the kinematic problems of high-energy physics'. W.A. Benjamin, New York, Amsterdam, 1964 -- H.Pilkhun, The Interactions of Hadrons North-Holland 1967 +- R.A. Kycia, P. Lebiedowicz, A. Szczurek, 'Decay: A Monte Carlo library for the decay of a particle with ROOT +compatibility', arXiv:2011.14750 [hep-ph] +- R. A. Kycia, J. Chwastowski, R. Staszewski, and J. Turnau, 'GenEx: A simple generator structure for exclusive +processes in high energy collisions', Commun. Comput. Phys.24no. 3, (2018) 860, arXiv:1411.6035 [hep-ph] +- R. A. Kycia, J. Turnau, J. J. Chwastowski, R. Staszewski, and M. Trzebinski, 'The adaptiveMonte Carlo toolbox for +phase space integration and generation', Commun. Comput. Phys.25no. 5, (2019) 1547, arXiv:1711.06087 [hep-ph] +- R. Hagedorn, 'Relativistic Kinematics: A guide to the kinematic problems of high-energy physics'. W.A. Benjamin, New +York, Amsterdam, 1964 +- H.Pilkhun, The Interactions of Hadrons North-Holland 1967 - S. Jadach, Foam: A General-Purpose Cellular Monte Carlo Event Generator, Comput.Phys.Commun. 152 (2003) 55-100 */ @@ -66,8 +86,6 @@ If you prepared the integrand to be integrated over LIPS, you can print the inte #ifndef ROOT_TGenFoamDecay #define ROOT_TGenFoamDecay - - #include #include @@ -76,126 +94,128 @@ If you prepared the integrand to be integrated over LIPS, you can print the inte #include #include -#include // std::cin, std::cout -#include // std::queue +#include // std::cin, std::cout +#include // std::queue #include #include using namespace std; - -class TGenFoamDecay : public TFoamIntegrand, public TGenInterface { -private: - Int_t fNt; // number of decay particles - Double_t fMass[18]; // masses of particles - Double_t fBeta[3]; // betas of decaying particle - Double_t fTeCmTm; // total energy in the C.M. minus the total mass - TLorentzVector fDecPro[18]; //kinematics of the generated particles - Int_t seed; //seed for pseudorandom generator - - Int_t kMAXP; //max number of particles (relict from TGenPhaseSpace) - - TDecay _decay; //decay engine - TFoam* _foam; //adaptive integrator - TRandom3 _pseRan; // pseudorandom number generator - - //FOAM parameters - Int_t nCells; // Number of Cells - Int_t nSampl; // Number of MC events per cell in build-up - Int_t nBin; // Number of bins in build-up - Int_t OptRej; // Wted events for OptRej=0; wt=1 for OptRej=1 (default) - Int_t OptDrive; // (D=2) Option, type of Drive =0,1,2 for TrueVol,Sigma,WtMax - Int_t EvPerBin; // Maximum events (equiv.) per bin in buid-up - Int_t Chat; // Chat level - - - - /// @returns weight of the process - /// @param nDim dimension of integration - /// @param Xarg vector of probablilistic variables from [0;1] of dim nDim - /// @warning it is required by Foam integrator - Double_t Density(int nDim, Double_t *Xarg); - - +class TGenFoamDecay : public TFoamIntegrand, public TGenInterface { +private: + Int_t fNt; // number of decay particles + Double_t fMass[18]; // masses of particles + Double_t fBeta[3]; // betas of decaying particle + Double_t fTeCmTm; // total energy in the C.M. minus the total mass + TLorentzVector fDecPro[18]; // kinematics of the generated particles + Int_t fSeed; // seed for pseudorandom generator + + Int_t kMAXP; // max number of particles (relict from TGenPhaseSpace) + + TDecay fDecay; // decay engine + TFoam *fFoam; // adaptive integrator + TRandom3 fPseRan; // pseudorandom number generator + + // FOAM parameters + Int_t fNCells; // Number of Cells + Int_t fNSampl; // Number of MC events per cell in build-up + Int_t fNBin; // Number of bins in build-up + Int_t fOptRej; // Wted events for OptRej=0; wt=1 for OptRej=1 (default) + Int_t fOptDrive; // (D=2) Option, type of Drive =0,1,2 for TrueVol,Sigma,WtMax + Int_t fEvPerBin; // Maximum events (equiv.) per bin in buid-up + Int_t fChat; // Chat level + + /// @returns weight of the process + /// @param nDim dimension of integration + /// @param Xarg vector of probablilistic variables from [0;1] of dim nDim + /// @warning it is required by Foam integrator + Double_t Density(int nDim, Double_t *Xarg); public: - /// constructor - TGenFoamDecay(): fNt(0), fMass(), fBeta(), fTeCmTm(0.), seed(4357), nCells(1000), nSampl(1000), nBin(8), OptRej(1), OptDrive(2), EvPerBin(25), Chat(1) {_foam = new TFoam("FoamX");}; + TGenFoamDecay() + : fNt(0), fMass(), fBeta(), fTeCmTm(0.), fSeed(4357), fNCells(1000), fNSampl(1000), fNBin(8), fOptRej(1), fOptDrive(2), + fEvPerBin(25), fChat(1) + { + fFoam = new TFoam("FoamX"); + }; /// copy constructor TGenFoamDecay(const TGenFoamDecay &gen); /// desctructor - virtual ~TGenFoamDecay() { delete _foam;} + virtual ~TGenFoamDecay() { delete fFoam; } /// assignemt - TGenFoamDecay& operator=(const TGenFoamDecay &gen); + TGenFoamDecay &operator=(const TGenFoamDecay &gen); /// Sets up configuration of decay /// @param P 4-momentum of decaying particle (Momentum, Energy units are Gev/C, GeV) /// @param nt number of products of decay /// @param mass mass matrix of products of decay mass[nt] /// @returns kTRUE - the decay is permitted by kinematics; kFALSE - the decay is forbidden by kinematics - /// @warning This should be first method to call since it sets up decay configuration. + /// @warning This should be first method to call since it sets up decay configuration. /// @warning The method also initialize FOAM. Bool_t SetDecay(TLorentzVector &P, Int_t nt, const Double_t *mass); - + /// Generate a single decay Double_t Generate(void); - + /// Collect 4-vector of products of decay /// @param n number of final particle to get from range 1...nt /// @warning You should call Generate() in first place. - TLorentzVector *GetDecay(Int_t n); + TLorentzVector *GetDecay(Int_t n); /// @returns 4-vector of n-th product of decay /// @param n number of final particle to get from range 1...nt - Int_t GetNt() const { return fNt;}; - - /// @returns the function under integral over LIPS (Lorentz Invariant Phase Space) - /// @param fNt number of outgoing particles - /// @param pf array of TLorentzVectors of outgoing particles - /// @warning It is set to 1.0. You should to redefine it (in your derived class, after inheritance) when you want to use full adaptation features. - virtual Double_t Integrand( int fNt, TLorentzVector* pf ); - - ///sets seed for pseudorandom number generator - virtual void setSeed( Int_t seed ) {seed = seed; _pseRan.SetSeed(seed);}; - - ///finalize Foam printing out MC integral and statistics - virtual void Finalize( void ); - - ///@returns integral +- error of MC integration of function in Integrand() over LIPS - virtual void GetIntegMC(Double_t & inetgral, Double_t & error); - - ///Sets FOAM number of Cells - /// @warning It should be done before Generate() - void SetnCells(Int_t nc) {nCells = nc;}; - - ///Sets FOAM number of MC events per cell in build-up - /// @warning It should be done before Generate() - void SetnSampl(Int_t ns) {nSampl = ns;}; - - ///Sets FOAM number of bins in build-up - /// @warning It should be done before Generate() - void SetnBin(Int_t nb) {nBin = nb;}; - - ///Sets FOAM Weigh events for OptRej=0; wt=1 for OptRej=1 (default) - /// @warning It determines if the events will be weighted of not (weight=1). - /// @warning It should be done before Generate() - void SetOptRej(Int_t OptR) {OptRej = OptR;}; - - ///Sets FOAM (D=2) option, type of Drive =0,1,2 for TrueVol,Sigma,WtMax - /// @warning It should be done before Generate() - void SetOptDrive(Int_t OptD) {OptDrive = OptD;}; - - ///Sets FOAM maximum events (equiv.) per bin in buid-up - /// @warning It should be done before Generate() - void SetEvPerBin(Int_t Ev) {EvPerBin = Ev;}; - - ///Sets FOAM chat level - /// @warning It should be done before Generate() - void SetChat(Int_t Ch) {Chat = Ch;}; - + Int_t GetNt() const { return fNt; }; + + /// @returns the function under integral over LIPS (Lorentz Invariant Phase Space) + /// @param fNt number of outgoing particles + /// @param pf array of TLorentzVectors of outgoing particles + /// @warning It is set to 1.0. You should to redefine it (in your derived class, after inheritance) when you want to + /// use full adaptation features. + virtual Double_t Integrand(int fNt, TLorentzVector *pf); + + /// sets seed for pseudorandom number generator + virtual void setSeed(Int_t seed) + { + fSeed = seed; + fPseRan.SetSeed(seed); + }; + + /// finalize Foam printing out MC integral and statistics + virtual void Finalize(void); + + ///@returns integral +- error of MC integration of function in Integrand() over LIPS + virtual void GetIntegMC(Double_t &inetgral, Double_t &error); + + /// Sets FOAM number of Cells + /// @warning It should be done before Generate() + void SetnCells(Int_t nc) { fNCells = nc; }; + + /// Sets FOAM number of MC events per cell in build-up + /// @warning It should be done before Generate() + void SetnSampl(Int_t ns) { fNSampl = ns; }; + + /// Sets FOAM number of bins in build-up + /// @warning It should be done before Generate() + void SetnBin(Int_t nb) { fNBin = nb; }; + + /// Sets FOAM Weigh events for OptRej=0; wt=1 for OptRej=1 (default) + /// @warning It determines if the events will be weighted of not (weight=1). + /// @warning It should be done before Generate() + void SetOptRej(Int_t OptR) { fOptRej = OptR; }; + + /// Sets FOAM (D=2) option, type of Drive =0,1,2 for TrueVol,Sigma,WtMax + /// @warning It should be done before Generate() + void SetOptDrive(Int_t OptD) { fOptDrive = OptD; }; + + /// Sets FOAM maximum events (equiv.) per bin in buid-up + /// @warning It should be done before Generate() + void SetEvPerBin(Int_t Ev) { fEvPerBin = Ev; }; + + /// Sets FOAM chat level + /// @warning It should be done before Generate() + void SetChat(Int_t Ch) { fChat = Ch; }; }; #endif - diff --git a/math/physics/src/TDecay.cxx b/math/physics/src/TDecay.cxx index d3c6eae47d854..bff6617312762 100644 --- a/math/physics/src/TDecay.cxx +++ b/math/physics/src/TDecay.cxx @@ -1,80 +1,79 @@ /************************************************************************* -* Copyright (C) 1995-2021, Rene Brun and Fons Rademakers. * -* All rights reserved. * -* * -* For the licensing terms see $ROOTSYS/LICENSE. * -* For the list of contributors see $ROOTSYS/README/CREDITS. * -*************************************************************************/ + * Copyright (C) 1995-2021, Rene Brun and Fons Rademakers. * + * All rights reserved. * + * * + * For the licensing terms see $ROOTSYS/LICENSE. * + * For the list of contributors see $ROOTSYS/README/CREDITS. * + *************************************************************************/ #include "TDecay.h" - -//const Int_t kMAXP = 18; +// const Int_t kMAXP = 18; //_____________________________________________________________________________________ -Double_t TDecay::PDK(Double_t a, Double_t b, Double_t c) +Double_t TDecay::PDK(Double_t a, Double_t b, Double_t c) { - //the PDK function - Double_t x = (a-b-c)*(a+b+c)*(a-b+c)*(a+b-c); - x = TMath::Sqrt(x)/(2*a); + // the PDK function + Double_t x = (a - b - c) * (a + b + c) * (a - b + c) * (a + b - c); + x = TMath::Sqrt(x) / (2 * a); return x; } +namespace { //_____________________________________________________________________________________ -Int_t DoubleMax(const void *a, const void *b) +Int_t DoubleMax(const void *a, const void *b) { - //special max function - Double_t aa = * ((Double_t *) a); - Double_t bb = * ((Double_t *) b); - if (aa > bb) return 1; - if (aa < bb) return -1; + // special max function + Double_t aa = *((Double_t *)a); + Double_t bb = *((Double_t *)b); + if (aa > bb) + return 1; + if (aa < bb) + return -1; return 0; - } +} // namespace //__________________________________________________________________________________________________ TDecay::TDecay(const TDecay &gen) : TObject(gen) { - //copy constructor - fNt = gen.fNt; - fTeCmTm = gen.fTeCmTm; + // copy constructor + fNt = gen.fNt; + fTeCmTm = gen.fTeCmTm; fBeta[0] = gen.fBeta[0]; fBeta[1] = gen.fBeta[1]; fBeta[2] = gen.fBeta[2]; - for (Int_t i=0;i &rnd) { @@ -82,36 +81,32 @@ Double_t TDecay::Generate(std::queue &rnd) Double_t rno[kMAXP]; rno[0] = 0; Int_t n; - - //check if phase space dim is ok - assert( (unsigned int) (3*fNt - 4) == rnd.size() ); - - double weight = 1.0; // additional weight - - if (fNt>2) - { - for (n=1; n 2) { + for (n = 1; n < fNt - 1; n++) { + rno[n] = rnd.front(); // fNt-2 random numbers + rnd.pop(); + } + + qsort(rno + 1, fNt - 2, sizeof(Double_t), DoubleMax); // sort them + + weight /= double(factorial(fNt - 2)); } - rno[fNt-1] = 1; - - - Double_t invMas[kMAXP], sum=0; - for (n=0; n &rnd) // Double_t wt = 1.0; Double_t pd[kMAXP]; - for (n=0; nPx(); Double_t z = v->Pz(); - v->SetPz( cY*z - sY*x ); - v->SetPx( sY*z + cY*x); // rotation around Y + v->SetPz(cY * z - sY * x); + v->SetPx(sY * z + cY * x); // rotation around Y x = v->Px(); Double_t y = v->Py(); - v->SetPx( cZ*x - sZ*y ); - v->SetPy( sZ*x + cZ*y ); // rotation around Z + v->SetPx(cZ * x - sZ * y); + v->SetPy(sZ * x + cZ * y); // rotation around Z } - if (i == (fNt-1)) break; + if (i == (fNt - 1)) + break; - Double_t beta = pd[i] / sqrt(pd[i]*pd[i] + invMas[i]*invMas[i]); - for (j=0; j<=i; j++) fDecPro[j].Boost(0,0,beta); + Double_t beta = pd[i] / sqrt(pd[i] * pd[i] + invMas[i] * invMas[i]); + for (j = 0; j <= i; j++) + fDecPro[j].Boost(0, 0, beta); i++; } - // - //---> final boost of all particles + //---> final boost of all particles // - for (n=0;n return the weigth of event @@ -176,43 +172,44 @@ Double_t TDecay::Generate(std::queue &rnd) } //__________________________________________________________________________________ -TLorentzVector *TDecay::GetDecay(Int_t n) -{ - //return Lorentz vector corresponding to decay n - if (n>fNt) return 0; - return fDecPro+n; +TLorentzVector *TDecay::GetDecay(Int_t n) +{ + // return Lorentz vector corresponding to decay n + if (n > fNt) + return 0; + return fDecPro + n; } //_____________________________________________________________________________________ -Bool_t TDecay::SetDecay(TLorentzVector &P, Int_t nt, const Double_t *mass) +Bool_t TDecay::SetDecay(TLorentzVector &P, Int_t nt, const Double_t *mass) { kMAXP = nt; - + Int_t n; fNt = nt; - if (fNt<2 || fNt>18) return kFALSE; // no more then 18 particle - + if (fNt < 2 || fNt > 18) + return kFALSE; // no more then 18 particle // - fTeCmTm = P.Mag(); // total energy in C.M. minus the sum of the masses - for (n=0;n save the betas of the decaying particle // if (P.Beta()) { - Double_t w = P.Beta()/P.Rho(); - fBeta[0] = P(0)*w; - fBeta[1] = P(1)*w; - fBeta[2] = P(2)*w; - } - else fBeta[0]=fBeta[1]=fBeta[2]=0; - - return kTRUE; + Double_t w = P.Beta() / P.Rho(); + fBeta[0] = P(0) * w; + fBeta[1] = P(1) * w; + fBeta[2] = P(2) * w; + } else + fBeta[0] = fBeta[1] = fBeta[2] = 0; + + return kTRUE; } diff --git a/math/physics/src/TGenDecay.cxx b/math/physics/src/TGenDecay.cxx index 37afea2efb8e3..bf9f9ddba84a4 100644 --- a/math/physics/src/TGenDecay.cxx +++ b/math/physics/src/TGenDecay.cxx @@ -1,109 +1,105 @@ /************************************************************************* -* Copyright (C) 1995-2021, Rene Brun and Fons Rademakers. * -* All rights reserved. * -* * -* For the licensing terms see $ROOTSYS/LICENSE. * -* For the list of contributors see $ROOTSYS/README/CREDITS. * -*************************************************************************/ + * Copyright (C) 1995-2021, Rene Brun and Fons Rademakers. * + * All rights reserved. * + * * + * For the licensing terms see $ROOTSYS/LICENSE. * + * For the list of contributors see $ROOTSYS/README/CREDITS. * + *************************************************************************/ #include "TGenDecay.h" - -//const Int_t kMAXP = 18; - +// const Int_t kMAXP = 18; //__________________________________________________________________________________________________ TGenDecay::TGenDecay(const TGenDecay &gen) : TObject(gen) { - //copy constructor - fNt = gen.fNt; - fTeCmTm = gen.fTeCmTm; + // copy constructor + fNt = gen.fNt; + fTeCmTm = gen.fTeCmTm; fBeta[0] = gen.fBeta[0]; fBeta[1] = gen.fBeta[1]; fBeta[2] = gen.fBeta[2]; - _decay = gen._decay; - _pseRan = gen._pseRan; - for (Int_t i=0;ifNt) return 0; - - return _decay.GetDecay(n); +TLorentzVector *TGenDecay::GetDecay(Int_t n) +{ + // return Lorentz vector corresponding to decay n + if (n > fNt) + return 0; + + return fDecay.GetDecay(n); } //_____________________________________________________________________________________ -Bool_t TGenDecay::SetDecay(TLorentzVector &P, Int_t nt, const Double_t *mass) +Bool_t TGenDecay::SetDecay(TLorentzVector &P, Int_t nt, const Double_t *mass) { kMAXP = nt; Int_t n; fNt = nt; - if (fNt<2 || fNt>18) return kFALSE; // no more then 18 particle + if (fNt < 2 || fNt > 18) + return kFALSE; // no more then 18 particle // - fTeCmTm = P.Mag(); // total energy in C.M. minus the sum of the masses - for (n=0;n rndQueue; + // put rnd numbers into queue + for (int i = 0; i < 3 * fNt - 4; i++) { + rndQueue.push(Xarg[i]); + } + // make decay and take d(LIPS) + double wtdecay = fDecay.Generate(rndQueue); -Double_t TGenFoamDecay::Density(int nDim, Double_t *Xarg) -{ - - //queue for random numbers - queue rndQueue; - - //put rnd numbers into queue - for( int i = 0; i < 3*fNt-4; i++) - { - rndQueue.push( Xarg[i] ); - } - - //make decay and take d(LIPS) - double wtdecay = _decay.Generate( rndQueue ); - - //get out particles - TLorentzVector pf[fNt]; - for( int i = 0; i < fNt; i++ ) - { - pf[i] = *(_decay.GetDecay( i )); - } - - //calculate integrand - double integrand = Integrand( fNt, pf ); - - return wtdecay * integrand; + // get out particles + TLorentzVector pf[fNt]; + for (int i = 0; i < fNt; i++) { + pf[i] = *(fDecay.GetDecay(i)); + } + + // calculate integrand + double integrand = Integrand(fNt, pf); + + return wtdecay * integrand; } //__________________________________________________________________________________________________ -Double_t TGenFoamDecay::Integrand( int fNt, TLorentzVector * pf ) +Double_t TGenFoamDecay::Integrand(int /*nt*/, TLorentzVector * /*pf*/) { - return 1.0; //default and probably overloaded for matrix element -} + return 1.0; // default and probably overloaded for matrix element +} //__________________________________________________________________________________________________ TGenFoamDecay::TGenFoamDecay(const TGenFoamDecay &gen) { - //copy constructor - fNt = gen.fNt; - fTeCmTm = gen.fTeCmTm; + // copy constructor + fNt = gen.fNt; + fTeCmTm = gen.fTeCmTm; fBeta[0] = gen.fBeta[0]; fBeta[1] = gen.fBeta[1]; fBeta[2] = gen.fBeta[2]; - _decay = gen._decay; - _foam = gen._foam; - _pseRan = gen._pseRan; - for (Int_t i=0;iMakeEvent(); - - return _foam->GetMCwt( ); + fFoam->MakeEvent(); + return fFoam->GetMCwt(); } //__________________________________________________________________________________ -TLorentzVector *TGenFoamDecay::GetDecay(Int_t n) -{ - - if (n>fNt) return 0; - - //return Lorentz vector corresponding to decay of n-th particle - return _decay.GetDecay( n ); +TLorentzVector *TGenFoamDecay::GetDecay(Int_t n) +{ + + if (n > fNt) + return 0; + + // return Lorentz vector corresponding to decay of n-th particle + return fDecay.GetDecay(n); } //_____________________________________________________________________________________ -Bool_t TGenFoamDecay::SetDecay(TLorentzVector &P, Int_t nt, - const Double_t *mass) +Bool_t TGenFoamDecay::SetDecay(TLorentzVector &P, Int_t nt, const Double_t *mass) { kMAXP = nt; Int_t n; fNt = nt; - if (fNt<2 || fNt>18) return kFALSE; // no more then 18 particle - + if (fNt < 2 || fNt > 18) + return kFALSE; // no more then 18 particle - fTeCmTm = P.Mag(); // total energy in C.M. minus the sum of the masses - for (n=0;n 0 ) - { - cout<<"***** Foam version "<< _foam->GetVersion() <<" *****"<SetkDim( 3*fNt-4); // Mandatory!!! - _foam->SetnCells( nCells); // optional - _foam->SetnSampl( nSampl); // optional - _foam->SetnBin( nBin); // optional - _foam->SetOptRej( OptRej); // optional - _foam->SetOptDrive( OptDrive); // optional - _foam->SetEvPerBin( EvPerBin); // optional - _foam->SetChat( Chat); // optional - //=============================== - _foam->SetRho(this); - _foam->SetPseRan(&_pseRan); - - // Initialize simulator - _foam->Initialize(); - - return kTRUE; + if (fTeCmTm <= 0) + return kFALSE; // not enough energy for this decay + + fDecay.SetDecay(P, fNt, fMass); // set decay to TDecay + + // initialize FOAM + //========================================================= + if (fChat > 0) { + cout << "***** Foam version " << fFoam->GetVersion() << " *****" << endl; + } + fFoam->SetkDim(3 * fNt - 4); // Mandatory!!! + fFoam->SetnCells(fNCells); // optional + fFoam->SetnSampl(fNSampl); // optional + fFoam->SetnBin(fNBin); // optional + fFoam->SetOptRej(fOptRej); // optional + fFoam->SetOptDrive(fOptDrive); // optional + fFoam->SetEvPerBin(fEvPerBin); // optional + fFoam->SetChat(fChat); // optional + //=============================== + fFoam->SetRho(this); + fFoam->SetPseRan(&fPseRan); + + // Initialize simulator + fFoam->Initialize(); + + return kTRUE; } //_____________________________________________________________________________________ -void TGenFoamDecay::Finalize( void ) +void TGenFoamDecay::Finalize(void) { - Double_t MCresult,MCerror; - Double_t eps = 0.0005; - Double_t Effic, WtMax, AveWt, Sigma; - Double_t IntNorm, Errel; - _foam->Finalize( IntNorm, Errel); // final printout - _foam->GetIntegMC( MCresult, MCerror); // get MC intnegral - _foam->GetWtParams(eps, AveWt, WtMax, Sigma); // get MC wt parameters - long nCalls=_foam->GetnCalls(); - Effic=0; if(WtMax>0) Effic=AveWt/WtMax; - cout << "================================================================" << endl; - cout << " MCresult= " << MCresult << " +- " << MCerror << " RelErr= "<< MCerror/MCresult << endl; - cout << " Dispersion/= " << Sigma/AveWt << endl; - cout << " /WtMax= " << Effic <<", for epsilon = "<Finalize(IntNorm, Errel); // final printout + fFoam->GetIntegMC(MCresult, MCerror); // get MC intnegral + fFoam->GetWtParams(eps, AveWt, WtMax, Sigma); // get MC wt parameters + long nCalls = fFoam->GetnCalls(); + Effic = 0; + if (WtMax > 0) + Effic = AveWt / WtMax; + cout << "================================================================" << endl; + cout << " MCresult= " << MCresult << " +- " << MCerror << " RelErr= " << MCerror / MCresult << endl; + cout << " Dispersion/= " << Sigma / AveWt << endl; + cout << " /WtMax= " << Effic << ", for epsilon = " << eps << endl; + cout << " nCalls (initialization only) = " << nCalls << endl; + cout << "================================================================" << endl; } //_____________________________________________________________________________________ -void TGenFoamDecay::GetIntegMC(Double_t & integral, Double_t & error) +void TGenFoamDecay::GetIntegMC(Double_t &integral, Double_t &error) { - _foam-> GetIntegMC( integral, error); + fFoam->GetIntegMC(integral, error); }