STARLIGHT/starlight/include/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 //
  20 ///////////////////////////////////////////////////////////////////////////
  21 //
  22 // File and Version Information:
  23 // $Rev:: 157                         $: revision of last commit
  24 // $Author:: odjuvsla                 $: author of last commit
  25 // $Date:: 2013-10-06 16:17:54 +0200 #$: 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 //
  66 ///////////////////////////////////////////////////////////////////////////
  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();
 113         virtual ~nBodyPhaseSpaceGen();
 114
 115         // generator setup
 116         /// sets decay constants and prepares internal variables
 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 randyInstance.Rndom(); }  ///< returns number from internal random generator
 123
 124         // high-level generator interface
 125         /// generates full event with certain n-body mass and momentum and returns event weight
 126         double generateDecay        (const lorentzVector& nBody);          // Lorentz vector of n-body system in lab frame
 127         /// \brief generates full event with certain n-body mass and momentum only when event is accepted (return value = true)
 128         /// this function is more efficient, if only weighted events are needed
 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_;    }  ///< sets maximum weight used for hit-miss MC
 133         double maxWeight             () const                 { return _maxWeight;         }  ///< returns maximum weight used for hit-miss MC
 134         double normalization         () const                 { return _norm;              }  ///< returns normalization used in weight calculation
 135         double eventWeight           () const                 { return _weight;            }  ///< returns weight of generated event
 136         double maxWeightObserved     () const                 { return _maxWeightObserved; }  ///< returns maximum observed weight since instantiation
 137         void   resetMaxWeightObserved()                       { _maxWeightObserved = 0;    }  ///< sets maximum observed weight back to zero
 138
 139         /// estimates maximum weight for given n-body mass
 140         double estimateMaxWeight(const double       nBodyMass,                 // sic!
 141                                  const unsigned int nmbOfIterations = 10000);  // number of generated events
 142
 143         /// \brief applies event weight in form of hit-miss MC
 144         /// assumes that event weight has been already calculated by calcWeight()
 145         /// if maxWeight > 0 value is used as maximum weight, otherwise _maxWeight value is used
 146         inline bool eventAccepted(const double maxWeight = 0);
 147
 148         //----------------------------------------------------------------------------
 149         // trivial accessors
 150         const lorentzVector&              daughter        (const int index) const { return _daughters[index];  }  ///< returns Lorentz vector of daughter at index
 151         const std::vector<lorentzVector>& daughters       ()                const { return _daughters;         }  ///< returns Lorentz vectors of all daughters
 152         unsigned int                      nmbOfDaughters  ()                const { return _n;                 }  ///< returns number of daughters
 153         double                            daughterMass    (const int index) const { return _m[index];          }  ///< returns invariant mass of daughter at index
 154         double                            intermediateMass(const int index) const { return _M[index];          }  ///< returns intermediate mass of (index + 1)-body system
 155         double                            breakupMom      (const int index) const { return _breakupMom[index]; }  ///< returns breakup momentum in (index + 1)-body RF
 156         double                            cosTheta        (const int index) const { return _cosTheta[index];   }  ///< returns polar angle in (index + 1)-body RF
 157         double                            phi             (const int index) const { return _phi[index];        }  ///< returns azimuth in (index + 1)-body RF
 158
 159
 160         std::ostream& print(std::ostream& out = std::cout) const;  ///< prints generator status
 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
 169         /// randomly choses the (n - 2) effective masses of the respective (i + 1)-body systems
 170         void pickMasses(const double nBodyMass);  // total energy of n-body system in its RF
 171
 172         /// \brief computes event weight and breakup momenta
 173         /// operates on vector of intermediate two-body masses prepared by pickMasses()
 174         double calcWeight();
 175
 176         /// randomly choses the (n - 1) polar and (n - 1) azimuthal angles in the respective (i + 1)-body RFs
 177         inline void pickAngles();
 178
 179         /// \brief calculates full event kinematics from the effective masses of the (i + 1)-body systems and the Lorentz vector of the decaying system
 180         /// uses the break-up momenta calculated by calcWeight() and angles from pickAngles()
 181         void calcEventKinematics(const lorentzVector& nBody);  // Lorentz vector of n-body system in lab frame
 182
 183         // external parameters
 184         std::vector<double> _m;  ///< masses of daughter particles
 185
 186         // internal variables
 187         unsigned int               _n;                  ///< number of daughter particles
 188         std::vector<double>        _M;                  ///< effective masses of (i + 1)-body systems
 189         std::vector<double>        _cosTheta;           ///< cosine of polar angle of the 2-body decay of the (i + 1)-body system
 190         std::vector<double>        _phi;                ///< azimuthal angle of the 2-body decay of the (i + 1)-body system
 191         std::vector<double>        _mSum;               ///< sums of daughter particle masses
 192         std::vector<double>        _breakupMom;         ///< breakup momenta for the two-body decays: (i + 1)-body --> daughter_(i + 1) + i-body
 193         std::vector<lorentzVector> _daughters;          ///< Lorentz vectors of the daughter particles
 194         double                     _norm;               ///< normalization value
 195         double                     _weight;             ///< phase space weight of generated event
 196         double                     _maxWeightObserved;  ///< maximum event weight calculated processing the input data
 197         double                     _maxWeight;          ///< maximum weight used to weight events in hit-miss MC
 198
 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 maxWeightArg)  // if maxWeightArg > 0, given value is used as maximum weight, otherwise _maxWeight
 217 {
 218         const double max = (maxWeightArg <= 0) ? _maxWeight : maxWeightArg;
 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