EIC code displayed by LXR |
|
|
|
|||
File indexing completed on 2024-06-26 07:06:20
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
| [ Source navigation ] | [ Diff markup ] | [ Identifier search ] | [ general search ] |
|
This page was automatically generated by the 2.3.5 LXR engine. The LXR team |
|