EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
nBodyPhaseSpaceGen.h
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file nBodyPhaseSpaceGen.h
1 
2 //
3 // Copyright 2010
4 //
5 // This file is part of starlight.
6 //
7 // starlight is free software: you can redistribute it and/or modify
8 // it under the terms of the GNU General Public License as published by
9 // the Free Software Foundation, either version 3 of the License, or
10 // (at your option) any later version.
11 //
12 // starlight is distributed in the hope that it will be useful,
13 // but WITHOUT ANY WARRANTY; without even the implied warranty of
14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 // GNU General Public License for more details.
16 //
17 // You should have received a copy of the GNU General Public License
18 // along with starlight. If not, see <http://www.gnu.org/licenses/>.
19 //
21 //
22 // File and Version Information:
23 // $Rev:: 211 $: revision of last commit
24 // $Author:: butter $: author of last commit
25 // $Date:: 2015-08-10 03:05:09 +0100 #$: date of last commit
26 //
27 // Description:
28 // calculates n-body phase space (constant matrix element) using various algorithms
29 //
30 // the n-body decay is split up into (n - 2) successive 2-body decays
31 // each 2-body decay is considered in its own center-of-mass frame thereby
32 // separating the mass from the (trivial) angular dependence
33 //
34 // the event is boosted into the same frame in which the n-body system is
35 // given
36 //
37 // based on:
38 // GENBOD (CERNLIB W515), see F. James, "Monte Carlo Phase Space", CERN 68-15 (1968)
39 // NUPHAZ, see M. M. Block, "Monte Carlo phase space evaluation", Comp. Phys. Commun. 69, 459 (1992)
40 // S. U. Chung, "Spin Formalism", CERN Yellow Report
41 // S. U. Chung et. al., "Diffractive Dissociation for COMPASS"
42 //
43 // index convention:
44 // - all vectors have the same size (= number of decay daughters)
45 // - index i corresponds to the respective value in the (i + 1)-body system: effective mass M, break-up momentum, angles
46 // - thus some vector elements are not used like breakupMom[0], theta[0], phi[0], ...
47 // this overhead is negligible compared to the ease of notation
48 //
49 // the following graph illustrates how the n-body decay is decomposed into a sequence of two-body decays
50 //
51 // n-body ... 3-body 2-body single daughter
52 //
53 // m[n - 1] m[2] m[1]
54 // ^ ^ ^
55 // | | |
56 // | | |
57 // M[n - 1] --> ... --> M[2] --> M[1] --> M [0] = m[0]
58 // theta[n - 1] ... theta[2] theta[1] theta[0] = 0 (not used)
59 // phi [n - 1] ... phi [2] phi [1] phi [0] = 0 (not used)
60 // mSum [n - 1] ... mSum [2] mSum [1] mSum [0] = m[0]
61 // = sum_0^(n - 1) m[i] = m[2] + m[1] + m[0] = m[1] + m[0]
62 // breakUpMom[n - 1] ... breakUpMom[2] breakUpMom[1] breakUpMom[0] = 0 (not used)
63 // = q(M[n - 1], m[n - 1], M[n - 2]) = q(M[2], m[2], M[1]) = q(M[1], m[1], m[0])
64 //
65 //
67 
68 
69 #ifndef NBODYPHASESPACEGEN_H
70 #define NBODYPHASESPACEGEN_H
71 
72 
73 #include <iostream>
74 #include <vector>
75 
76 #include "reportingUtils.h"
77 #include "lorentzvector.h"
78 #include "randomgenerator.h"
79 #include "starlightconstants.h"
80 
81 
82 // small helper functions
83 // calculates factorial
84 inline
85 unsigned int
86 factorial(const unsigned int n)
87 {
88  unsigned int fac = 1;
89  for (unsigned int i = 1; i <= n; ++i)
90  fac *= i;
91  return fac;
92 }
93 
94 
95 // computes breakup momentum of 2-body decay
96 inline
97 double
98 breakupMomentum(const double M, // mass of mother particle
99  const double m1, // mass of daughter particle 1
100  const double m2) // mass of daughter particle 2
101 {
102  if (M < m1 + m2)
103  return 0;
104  return sqrt((M - m1 - m2) * (M + m1 + m2) * (M - m1 + m2) * (M + m1 - m2)) / (2 * M);
105 }
106 
107 
108 class nBodyPhaseSpaceGen {
109 
110 public:
111 
112  nBodyPhaseSpaceGen(const randomGenerator& randy);
113  virtual ~nBodyPhaseSpaceGen();
114 
115  // generator setup
117  bool setDecay(const std::vector<double>& daughterMasses); // daughter particle masses
118  bool setDecay(const unsigned int nmbOfDaughters, // number of daughter particles
119  const double* daughterMasses); // array of daughter particle masses
120 
121  // random generator
122  double random () { return _randy.Rndom(); }
123 
124  // high-level generator interface
126  double generateDecay (const lorentzVector& nBody); // Lorentz vector of n-body system in lab frame
129  bool generateDecayAccepted(const lorentzVector& nBody, // Lorentz vector of n-body system in lab frame
130  const double maxWeight = 0); // if positive, given value is used as maximum weight, otherwise _maxWeight
131 
132  void setMaxWeight (const double maxWeight) { _maxWeight = maxWeight; }
133  double maxWeight () const { return _maxWeight; }
134  double normalization () const { return _norm; }
135  double eventWeight () const { return _weight; }
136  double maxWeightObserved () const { return _maxWeightObserved; }
137  void resetMaxWeightObserved() { _maxWeightObserved = 0; }
138 
140  double estimateMaxWeight(const double nBodyMass, // sic!
141  const unsigned int nmbOfIterations = 10000); // number of generated events
142 
146  inline bool eventAccepted(const double maxWeight = 0);
147 
148  //----------------------------------------------------------------------------
149  // trivial accessors
150  const lorentzVector& daughter (const int index) const { return _daughters[index]; }
151  const std::vector<lorentzVector>& daughters () const { return _daughters; }
152  unsigned int nmbOfDaughters () const { return _n; }
153  double daughterMass (const int index) const { return _m[index]; }
154  double intermediateMass(const int index) const { return _M[index]; }
155  double breakupMom (const int index) const { return _breakupMom[index]; }
156  double cosTheta (const int index) const { return _cosTheta[index]; }
157  double phi (const int index) const { return _phi[index]; }
158 
159 
160  std::ostream& print(std::ostream& out = std::cout) const;
161  friend std::ostream& operator << (std::ostream& out,
162  const nBodyPhaseSpaceGen& gen)
163  { return gen.print(out); }
164 
165 private:
166 
167  //----------------------------------------------------------------------------
168  // low-level generator interface
170  void pickMasses(const double nBodyMass); // total energy of n-body system in its RF
171 
174  double calcWeight();
175 
177  inline void pickAngles();
178 
181  void calcEventKinematics(const lorentzVector& nBody); // Lorentz vector of n-body system in lab frame
182 
183  // external parameters
184  std::vector<double> _m;
185 
186  // internal variables
187  unsigned int _n;
188  std::vector<double> _M;
189  std::vector<double> _cosTheta;
190  std::vector<double> _phi;
191  std::vector<double> _mSum;
192  std::vector<double> _breakupMom;
193  std::vector<lorentzVector> _daughters;
194  double _norm;
195  double _weight;
196  double _maxWeightObserved;
197  double _maxWeight;
198  randomGenerator _randy;
199 
200 };
201 
202 
203 inline
204 void
205 nBodyPhaseSpaceGen::pickAngles()
206 {
207  for (unsigned int i = 1; i < _n; ++i) { // loop over 2- to n-bodies
208  _cosTheta[i] = 2 * random() - 1; // range [-1, 1]
209  _phi[i] = starlightConstants::twoPi * random(); // range [ 0, 2 pi]
210  }
211 }
212 
213 
214 inline
215 bool
216 nBodyPhaseSpaceGen::eventAccepted(const double maxWeight) // if maxWeight > 0, given value is used as maximum weight, otherwise _maxWeight
217 {
218  const double max = (maxWeight <= 0) ? _maxWeight : maxWeight;
219  if (max <= 0) {
220  printWarn << "maximum weight = " << max << " does not make sense. rejecting event." << std::endl;
221  return false;
222  }
223  if ((_weight / max) > random())
224  return true;
225  return false;
226 }
227 
228 
229 #endif // NBODYPHASESPACEGEN_H