Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2024-06-17 07:07:00

0001 ///////////////////////////////////////////////////////////////////////////
0002 //
0003 //    Copyright 2010
0004 //
0005 //    This file is part of starlight.
0006 //
0007 //    starlight is free software: you can redistribute it and/or modify
0008 //    it under the terms of the GNU General Public License as published by
0009 //    the Free Software Foundation, either version 3 of the License, or
0010 //    (at your option) any later version.
0011 //    
0012 //    starlight is distributed in the hope that it will be useful,
0013 //    but WITHOUT ANY WARRANTY; without even the implied warranty of
0014 //    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0015 //    GNU General Public License for more details.
0016 //    
0017 //    You should have received a copy of the GNU General Public License
0018 //    along with starlight. If not, see <http://www.gnu.org/licenses/>.
0019 //
0020 ///////////////////////////////////////////////////////////////////////////
0021 //
0022 // File and Version Information:
0023 // $Rev:: 211                         $: revision of last commit
0024 // $Author:: butter                   $: author of last commit
0025 // $Date:: 2015-08-10 03:05:09 +0100 #$: date of last commit
0026 //
0027 // Description:
0028 //    calculates n-body phase space (constant matrix element) using various algorithms
0029 //    
0030 //    the n-body decay is split up into (n - 2) successive 2-body decays
0031 //    each 2-body decay is considered in its own center-of-mass frame thereby
0032 //    separating the mass from the (trivial) angular dependence
0033 //    
0034 //    the event is boosted into the same frame in which the n-body system is
0035 //    given
0036 //    
0037 //    based on:
0038 //    GENBOD (CERNLIB W515), see F. James, "Monte Carlo Phase Space", CERN 68-15 (1968)
0039 //    NUPHAZ, see M. M. Block, "Monte Carlo phase space evaluation", Comp. Phys. Commun. 69, 459 (1992)
0040 //    S. U. Chung, "Spin Formalism", CERN Yellow Report
0041 //    S. U. Chung et. al., "Diffractive Dissociation for COMPASS"
0042 //    
0043 //    index convention:
0044 //    - all vectors have the same size (= number of decay daughters)
0045 //    - index i corresponds to the respective value in the (i + 1)-body system: effective mass M, break-up momentum, angles
0046 //    - thus some vector elements are not used like breakupMom[0], theta[0], phi[0], ...
0047 //      this overhead is negligible compared to the ease of notation
0048 //    
0049 //    the following graph illustrates how the n-body decay is decomposed into a sequence of two-body decays
0050 //    
0051 //    n-body       ...                   3-body                 2-body                 single daughter
0052 //    
0053 //    m[n - 1]                           m[2]                   m[1]
0054 //     ^                                  ^                      ^
0055 //     |                                  |                      |
0056 //     |                                  |                      |
0057 //    M[n - 1] --> ... -->               M[2] -->               M[1] -->               M    [0] = m[0]
0058 //    theta[n - 1] ...                   theta[2]               theta[1]               theta[0] = 0 (not used)
0059 //    phi  [n - 1] ...                   phi  [2]               phi  [1]               phi  [0] = 0 (not used)
0060 //    mSum [n - 1] ...                   mSum [2]               mSum [1]               mSum [0] = m[0]
0061 //    = sum_0^(n - 1) m[i]               = m[2] + m[1] + m[0]   = m[1] + m[0]
0062 //    breakUpMom[n - 1] ...              breakUpMom[2]          breakUpMom[1]          breakUpMom[0] = 0 (not used)
0063 //    = q(M[n - 1], m[n - 1], M[n - 2])  = q(M[2], m[2], M[1])  = q(M[1], m[1], m[0])
0064 //
0065 //
0066 ///////////////////////////////////////////////////////////////////////////
0067 
0068 
0069 #ifndef NBODYPHASESPACEGEN_H
0070 #define NBODYPHASESPACEGEN_H
0071 
0072 
0073 #include <iostream>
0074 #include <vector>
0075 
0076 #include "reportingUtils.h"
0077 #include "lorentzvector.h"
0078 #include "randomgenerator.h"
0079 #include "starlightconstants.h"
0080 
0081 
0082 // small helper functions
0083 // calculates factorial
0084 inline
0085 unsigned int
0086 factorial(const unsigned int n)
0087 {
0088     unsigned int fac = 1;
0089     for (unsigned int i = 1; i <= n; ++i)
0090         fac *= i;
0091     return fac;
0092 }
0093 
0094 
0095 // computes breakup momentum of 2-body decay
0096 inline
0097 double
0098 breakupMomentum(const double M,   // mass of mother particle
0099                 const double m1,  // mass of daughter particle 1
0100                 const double m2)  // mass of daughter particle 2
0101 {
0102   if (M < m1 + m2)
0103     return 0;
0104   return sqrt((M - m1 - m2) * (M + m1 + m2) * (M - m1 + m2) * (M + m1 - m2)) / (2 * M);
0105 }
0106 
0107 
0108 class nBodyPhaseSpaceGen {
0109 
0110 public:
0111 
0112     nBodyPhaseSpaceGen(const randomGenerator& randy);
0113     virtual ~nBodyPhaseSpaceGen();
0114   
0115     // generator setup
0116     /// sets decay constants and prepares internal variables
0117     bool setDecay(const std::vector<double>& daughterMasses);  // daughter particle masses
0118     bool setDecay(const unsigned int         nmbOfDaughters,   // number of daughter particles
0119                   const double*              daughterMasses);  // array of daughter particle masses
0120   
0121     // random generator
0122     double random ()                        { return _randy.Rndom(); }  ///< returns number from internal random generator
0123 
0124     // high-level generator interface
0125     /// generates full event with certain n-body mass and momentum and returns event weight
0126     double generateDecay        (const lorentzVector& nBody);          // Lorentz vector of n-body system in lab frame
0127     /// \brief generates full event with certain n-body mass and momentum only when event is accepted (return value = true)
0128     /// this function is more efficient, if only weighted events are needed
0129     bool   generateDecayAccepted(const lorentzVector& nBody,           // Lorentz vector of n-body system in lab frame
0130                                  const double         maxWeight = 0);  // if positive, given value is used as maximum weight, otherwise _maxWeight
0131 
0132     void   setMaxWeight          (const double maxWeight) { _maxWeight = maxWeight;    }  ///< sets maximum weight used for hit-miss MC
0133     double maxWeight             () const                 { return _maxWeight;         }  ///< returns maximum weight used for hit-miss MC
0134     double normalization         () const                 { return _norm;              }  ///< returns normalization used in weight calculation
0135     double eventWeight           () const                 { return _weight;            }  ///< returns weight of generated event
0136     double maxWeightObserved     () const                 { return _maxWeightObserved; }  ///< returns maximum observed weight since instantiation
0137     void   resetMaxWeightObserved()                       { _maxWeightObserved = 0;    }  ///< sets maximum observed weight back to zero
0138 
0139     /// estimates maximum weight for given n-body mass
0140     double estimateMaxWeight(const double       nBodyMass,                 // sic!
0141                              const unsigned int nmbOfIterations = 10000);  // number of generated events
0142 
0143     /// \brief applies event weight in form of hit-miss MC
0144     /// assumes that event weight has been already calculated by calcWeight()
0145     /// if maxWeight > 0 value is used as maximum weight, otherwise _maxWeight value is used
0146     inline bool eventAccepted(const double maxWeight = 0);
0147 
0148     //----------------------------------------------------------------------------
0149     // trivial accessors
0150     const lorentzVector&              daughter        (const int index) const { return _daughters[index];  }  ///< returns Lorentz vector of daughter at index
0151     const std::vector<lorentzVector>& daughters       ()                const { return _daughters;         }  ///< returns Lorentz vectors of all daughters
0152     unsigned int                      nmbOfDaughters  ()                const { return _n;                 }  ///< returns number of daughters
0153     double                            daughterMass    (const int index) const { return _m[index];          }  ///< returns invariant mass of daughter at index
0154     double                            intermediateMass(const int index) const { return _M[index];          }  ///< returns intermediate mass of (index + 1)-body system
0155     double                            breakupMom      (const int index) const { return _breakupMom[index]; }  ///< returns breakup momentum in (index + 1)-body RF
0156     double                            cosTheta        (const int index) const { return _cosTheta[index];   }  ///< returns polar angle in (index + 1)-body RF
0157     double                            phi             (const int index) const { return _phi[index];        }  ///< returns azimuth in (index + 1)-body RF
0158 
0159 
0160     std::ostream& print(std::ostream& out = std::cout) const;  ///< prints generator status
0161     friend std::ostream& operator << (std::ostream&             out,
0162                                       const nBodyPhaseSpaceGen& gen)
0163     { return gen.print(out); }
0164 
0165 private:
0166 
0167     //----------------------------------------------------------------------------
0168     // low-level generator interface
0169     /// randomly choses the (n - 2) effective masses of the respective (i + 1)-body systems
0170     void pickMasses(const double nBodyMass);  // total energy of n-body system in its RF
0171 
0172     /// \brief computes event weight and breakup momenta
0173     /// operates on vector of intermediate two-body masses prepared by pickMasses()
0174     double calcWeight();
0175 
0176     /// randomly choses the (n - 1) polar and (n - 1) azimuthal angles in the respective (i + 1)-body RFs
0177     inline void pickAngles();
0178         
0179     /// \brief calculates full event kinematics from the effective masses of the (i + 1)-body systems and the Lorentz vector of the decaying system
0180     /// uses the break-up momenta calculated by calcWeight() and angles from pickAngles()
0181     void calcEventKinematics(const lorentzVector& nBody);  // Lorentz vector of n-body system in lab frame
0182 
0183     // external parameters
0184     std::vector<double> _m;  ///< masses of daughter particles
0185 
0186     // internal variables
0187     unsigned int               _n;                  ///< number of daughter particles
0188     std::vector<double>        _M;                  ///< effective masses of (i + 1)-body systems
0189     std::vector<double>        _cosTheta;           ///< cosine of polar angle of the 2-body decay of the (i + 1)-body system
0190     std::vector<double>        _phi;                ///< azimuthal angle of the 2-body decay of the (i + 1)-body system
0191     std::vector<double>        _mSum;               ///< sums of daughter particle masses
0192     std::vector<double>        _breakupMom;         ///< breakup momenta for the two-body decays: (i + 1)-body --> daughter_(i + 1) + i-body
0193     std::vector<lorentzVector> _daughters;          ///< Lorentz vectors of the daughter particles
0194     double                     _norm;               ///< normalization value
0195     double                     _weight;             ///< phase space weight of generated event
0196     double                     _maxWeightObserved;  ///< maximum event weight calculated processing the input data
0197     double                     _maxWeight;          ///< maximum weight used to weight events in hit-miss MC
0198     randomGenerator        _randy;
0199 
0200 };
0201 
0202 
0203 inline
0204 void
0205 nBodyPhaseSpaceGen::pickAngles()
0206 {
0207     for (unsigned int i = 1; i < _n; ++i) {  // loop over 2- to n-bodies
0208         _cosTheta[i] = 2 * random() - 1;  // range [-1,    1]
0209         _phi[i]      = starlightConstants::twoPi * random();  // range [ 0, 2 pi]
0210     }
0211 }
0212 
0213 
0214 inline 
0215 bool
0216 nBodyPhaseSpaceGen::eventAccepted(const double maxWeight)  // if maxWeight > 0, given value is used as maximum weight, otherwise _maxWeight
0217 {
0218     const double max = (maxWeight <= 0) ? _maxWeight : maxWeight;
0219     if (max <= 0) {
0220         printWarn << "maximum weight = " << max << " does not make sense. rejecting event." << std::endl;
0221         return false;
0222     }
0223     if ((_weight / max) > random())
0224         return true;
0225     return false;
0226 }
0227 
0228 
0229 #endif  // NBODYPHASESPACEGEN_H