estarlight/src/nBodyPhaseSpaceGen.cpp

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 //     see nBodyPhaseSpaceGen.h
0029 //
0030 //
0031 ///////////////////////////////////////////////////////////////////////////
0032
0033
0034 #include <algorithm>
0035
0036 #include "nBodyPhaseSpaceGen.h"
0037
0038
0039 using namespace std;
0040 using namespace starlightConstants;
0041
0042
0043 nBodyPhaseSpaceGen::nBodyPhaseSpaceGen(const randomGenerator& randy)
0044     : _n                (0),
0045       _norm             (0),
0046       _weight           (0),
0047       _maxWeightObserved(0),
0048       _maxWeight        (0),
0049           _randy            (randy)
0050 { }
0051
0052
0053 nBodyPhaseSpaceGen::~nBodyPhaseSpaceGen()
0054 { }
0055
0056
0057 // sets decay constants and prepares internal variables
0058 bool
0059 nBodyPhaseSpaceGen::setDecay(const vector<double>& daughterMasses)  // array of daughter particle masses
0060 {
0061     _n = daughterMasses.size();
0062     if (_n < 2) {
0063         printWarn << "number of daughters = " << _n << " does not make sense." << endl;
0064         return false;
0065     }
0066     // copy daughter masses
0067     _m.clear();
0068     _m = daughterMasses;
0069     // prepare effective mass vector
0070     _M.clear();
0071     _M.resize(_n, 0);
0072     _M[0] = _m[0];
0073     // prepare angle vectors
0074     _cosTheta.clear();
0075     _cosTheta.resize(_n, 0);
0076     _phi.clear();
0077     _phi.resize(_n, 0);
0078     // calculate daughter mass sums
0079     _mSum.clear();
0080     _mSum.resize(_n, 0);
0081     _mSum[0] = _m[0];
0082     for (unsigned int i = 1; i < _n; ++i)
0083         _mSum[i] = _mSum[i - 1] + _m[i];
0084     // prepare breakup momentum vector
0085     _breakupMom.clear();
0086     _breakupMom.resize(_n, 0);
0087     // prepare vector for daughter Lorentz vectors
0088     _daughters.clear();
0089     _daughters.resize(_n, lorentzVector(0, 0, 0, 0));
0090     // calculate normalization
0091     _norm = 1 / (2 * pow(twoPi, 2 * (int)_n - 3) * factorial(_n - 2));
0092     resetMaxWeightObserved();
0093     return true;
0094 }
0095
0096
0097 // set decay constants and prepare internal variables
0098 bool
0099 nBodyPhaseSpaceGen::setDecay(const unsigned int nmbOfDaughters,  // number of daughter particles
0100                              const double*      daughterMasses)  // array of daughter particle masses
0101 {
0102     vector <double> m;
0103     m.resize(nmbOfDaughters, 0);
0104     for (unsigned int i = 0; i < nmbOfDaughters; ++i)
0105         m[i] = daughterMasses[i];
0106     return setDecay(m);
0107 }
0108
0109
0110 // generates event with certain n-body mass and momentum and returns event weigth
0111 // general purpose function
0112 double
0113 nBodyPhaseSpaceGen::generateDecay(const lorentzVector& nBody)  // Lorentz vector of n-body system in lab frame
0114 {
0115     const double nBodyMass = nBody.M();
0116     if (_n < 2) {
0117         printWarn << "number of daughter particles = " << _n << " is smaller than 2. "
0118                   << "weight is set to 0." << endl;
0119         _weight = 0;
0120     } else if (nBodyMass < _mSum[_n - 1]) {
0121         printWarn << "n-body mass = " << nBodyMass << " is smaller than sum of daughter masses = "
0122                   << _mSum[_n - 1] << ". weight is set to 0." << endl;
0123         _weight = 0;
0124     } else {
0125         pickMasses(nBodyMass);
0126         calcWeight();
0127         pickAngles();
0128         calcEventKinematics(nBody);
0129     }
0130     return _weight;
0131 }
0132
0133
0134 // generates full event with certain n-body mass and momentum only, when event is accepted (return value = true)
0135 // this function is more efficient, if only weighted evens are needed
0136 bool
0137 nBodyPhaseSpaceGen::generateDecayAccepted(const lorentzVector& nBody,      // Lorentz vector of n-body system in lab frame
0138                                           const double         maxWeight)  // if positive, given value is used as maximum weight, otherwise _maxWeight
0139 {
0140     const double nBodyMass = nBody.M();
0141     if (_n < 2) {
0142         printWarn << "number of daughter particles = " << _n << " is smaller than 2. "
0143                   << "no event generated." << endl;
0144         return false;
0145     } else if (nBodyMass < _mSum[_n - 1]) {
0146         printWarn << "n-body mass = " << nBodyMass << " is smaller than sum of daughter masses = "
0147                   << _mSum[_n - 1] << ". no event generated." << endl;
0148         return false;
0149     }
0150     pickMasses(nBodyMass);
0151     calcWeight();
0152     if (!eventAccepted(maxWeight))
0153         return false;
0154     pickAngles();
0155     calcEventKinematics(nBody);
0156     return true;
0157 }
0158
0159
0160 // randomly choses the (n - 2) effective masses of the respective (i + 1)-body systems
0161 void
0162 nBodyPhaseSpaceGen::pickMasses(const double nBodyMass)  // total energy of the system in its RF
0163 {
0164     _M[_n - 1] = nBodyMass;
0165     // create vector of sorted random values
0166     vector<double> r(_n - 2, 0);  // (n - 2) values needed for 2- through (n - 1)-body systems
0167     for (unsigned int i = 0; i < (_n - 2); ++i)
0168         r[i] = random();
0169     sort(r.begin(), r.end());
0170     // set effective masses of (intermediate) two-body decays
0171     const double massInterval = nBodyMass - _mSum[_n - 1];  // kinematically allowed mass interval
0172     for (unsigned int i = 1; i < (_n - 1); ++i)             // loop over intermediate 2- to (n - 1)-bodies
0173         _M[i] = _mSum[i] + r[i - 1] * massInterval;           // _mSum[i] is minimum effective mass
0174 }
0175
0176
0177 // computes event weight (= integrand value) and breakup momenta
0178 // uses vector of intermediate two-body masses prepared by pickMasses()
0179 double
0180 nBodyPhaseSpaceGen::calcWeight()
0181 {
0182     for (unsigned int i = 1; i < _n; ++i)  // loop over 2- to n-bodies
0183         _breakupMom[i] = breakupMomentum(_M[i], _M[i - 1], _m[i]);
0184     double momProd = 1;                    // product of breakup momenta
0185     for (unsigned int i = 1; i < _n; ++i)  // loop over 2- to n-bodies
0186         momProd *= _breakupMom[i];
0187     const double massInterval = _M[_n - 1] - _mSum[_n - 1];  // kinematically allowed mass interval
0188     _weight = _norm * pow(massInterval, (int)_n - 2) * momProd / _M[_n - 1];
0189     if (_weight > _maxWeightObserved)
0190         _maxWeightObserved = _weight;
0191     if (std::isnan(_weight))
0192         printWarn << "weight = " << _weight << endl;
0193     return _weight;
0194 }
0195
0196
0197 // calculates complete event from the effective masses of the (i + 1)-body
0198 // systems, the Lorentz vector of the decaying system, and the decay angles
0199 // uses the break-up momenta calculated by calcWeight()
0200 void
0201 nBodyPhaseSpaceGen::calcEventKinematics(const lorentzVector& nBody)  // Lorentz vector of n-body system in lab frame
0202 {
0203     // build event starting in n-body RF going down to 2-body RF
0204     // is more efficicient than Raubold-Lynch method, since it requitres only one rotation and boost per daughter
0205     lorentzVector P = nBody;  // Lorentz of (i + 1)-body system in lab frame
0206     for (unsigned int i = _n - 1; i >= 1; --i) {  // loop from n-body down to 2-body
0207         // construct Lorentz vector of daughter _m[i] in (i + 1)-body RF
0208         const double   sinTheta = sqrt(1 - _cosTheta[i] * _cosTheta[i]);
0209         const double   pT       = _breakupMom[i] * sinTheta;
0210         lorentzVector& daughter = _daughters[i];
0211         daughter.SetPxPyPzE(pT * cos(_phi[i]),
0212                             pT * sin(_phi[i]),
0213                             _breakupMom[i] * _cosTheta[i],
0214                             sqrt(_m[i] * _m[i] + _breakupMom[i] * _breakupMom[i]));
0215         // boost daughter into lab frame
0216         daughter.Boost(P.BoostVector());
0217         // calculate Lorentz vector of i-body system in lab frame
0218         P -= daughter;
0219     }
0220     // set last daughter
0221     _daughters[0] = P;
0222 }
0223
0224
0225 // calculates maximum weight for given n-body mass
0226 double
0227 nBodyPhaseSpaceGen::estimateMaxWeight(const double       nBodyMass,        // sic!
0228                                       const unsigned int nmbOfIterations)  // number of generated events
0229 {
0230     double maxWeight = 0;
0231     for (unsigned int i = 0; i < nmbOfIterations; ++i) {
0232         pickMasses(nBodyMass);
0233         calcWeight();
0234         maxWeight = max(_weight, maxWeight);
0235     }
0236     return maxWeight;
0237 }
0238
0239
0240 ostream&
0241 nBodyPhaseSpaceGen::print(ostream& out) const
0242 {
0243     out << "nBodyPhaseSpaceGen parameters:" << endl
0244         << "    number of daughter particles ............... " << _n                 << endl
0245         << "    masses of the daughter particles ........... " << _m                 << endl
0246         << "    sums of daughter particle masses ........... " << _mSum              << endl
0247         << "    effective masses of (i + 1)-body systems ... " << _M                 << endl
0248         << "    cos(polar angle) in (i + 1)-body systems ... " << _cosTheta          << endl
0249         << "    azimuth in (i + 1)-body systems ............ " << _phi               << endl
0250         << "    breakup momenta in (i + 1)-body systems .... " << _breakupMom        << endl
0251         << "    normalization value ........................ " << _norm              << endl
0252         << "    weight of generated event .................. " << _weight            << endl
0253         << "    maximum weight used in hit-miss MC ......... " << _maxWeight         << endl
0254         << "    maximum weight since instantiation ......... " << _maxWeightObserved << endl
0255         << "    daughter four-momenta:" << endl;
0256     for (unsigned int i = 0; i < _n; ++i)
0257         out << "        daughter " << i << ": " << _daughters[i] << endl;
0258     return out;
0259 }