+++ /dev/null
-///////////////////////////////////////////////////////////////////////////
-//
-// Copyright 2010
-//
-// This file is part of starlight.
-//
-// starlight is free software: you can redistribute it and/or modify
-// it under the terms of the GNU General Public License as published by
-// the Free Software Foundation, either version 3 of the License, or
-// (at your option) any later version.
-//
-// starlight is distributed in the hope that it will be useful,
-// but WITHOUT ANY WARRANTY; without even the implied warranty of
-// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-// GNU General Public License for more details.
-//
-// You should have received a copy of the GNU General Public License
-// along with starlight. If not, see <http://www.gnu.org/licenses/>.
-//
-///////////////////////////////////////////////////////////////////////////
-//
-// File and Version Information:
-// $Rev:: $: revision of last commit
-// $Author:: $: author of last commit
-// $Date:: $: date of last commit
-//
-// Description:
-// see nBodyPhaseSpaceGen.h
-//
-//
-///////////////////////////////////////////////////////////////////////////
-
-
-#include <algorithm>
-
-#include "nBodyPhaseSpaceGen.h"
-
-
-using namespace std;
-using namespace starlightConstants;
-
-
-nBodyPhaseSpaceGen::nBodyPhaseSpaceGen()
- : _n (0),
- _norm (0),
- _weight (0),
- _maxWeightObserved(0),
- _maxWeight (0)
-{ }
-
-
-nBodyPhaseSpaceGen::~nBodyPhaseSpaceGen()
-{ }
-
-
-// sets decay constants and prepares internal variables
-bool
-nBodyPhaseSpaceGen::setDecay(const vector<double>& daughterMasses) // array of daughter particle masses
-{
- _n = daughterMasses.size();
- if (_n < 2) {
- printWarn << "number of daughters = " << _n << " does not make sense." << endl;
- return false;
- }
- // copy daughter masses
- _m.clear();
- _m = daughterMasses;
- // prepare effective mass vector
- _M.clear();
- _M.resize(_n, 0);
- _M[0] = _m[0];
- // prepare angle vectors
- _cosTheta.clear();
- _cosTheta.resize(_n, 0);
- _phi.clear();
- _phi.resize(_n, 0);
- // calculate daughter mass sums
- _mSum.clear();
- _mSum.resize(_n, 0);
- _mSum[0] = _m[0];
- for (unsigned int i = 1; i < _n; ++i)
- _mSum[i] = _mSum[i - 1] + _m[i];
- // prepare breakup momentum vector
- _breakupMom.clear();
- _breakupMom.resize(_n, 0);
- // prepare vector for daughter Lorentz vectors
- _daughters.clear();
- _daughters.resize(_n, lorentzVector(0, 0, 0, 0));
- // calculate normalization
- _norm = 1 / (2 * pow(twoPi, 2 * (int)_n - 3) * factorial(_n - 2));
- resetMaxWeightObserved();
- return true;
-}
-
-
-// set decay constants and prepare internal variables
-bool
-nBodyPhaseSpaceGen::setDecay(const unsigned int nmbOfDaughters, // number of daughter particles
- const double* daughterMasses) // array of daughter particle masses
-{
- vector <double> m;
- m.resize(nmbOfDaughters, 0);
- for (unsigned int i = 0; i < nmbOfDaughters; ++i)
- m[i] = daughterMasses[i];
- return setDecay(m);
-}
-
-
-// generates event with certain n-body mass and momentum and returns event weigth
-// general purpose function
-double
-nBodyPhaseSpaceGen::generateDecay(const lorentzVector& nBody) // Lorentz vector of n-body system in lab frame
-{
- const double nBodyMass = nBody.M();
- if (_n < 2) {
- printWarn << "number of daughter particles = " << _n << " is smaller than 2. "
- << "weight is set to 0." << endl;
- _weight = 0;
- } else if (nBodyMass < _mSum[_n - 1]) {
- printWarn << "n-body mass = " << nBodyMass << " is smaller than sum of daughter masses = "
- << _mSum[_n - 1] << ". weight is set to 0." << endl;
- _weight = 0;
- } else {
- pickMasses(nBodyMass);
- calcWeight();
- pickAngles();
- calcEventKinematics(nBody);
- }
- return _weight;
-}
-
-
-// generates full event with certain n-body mass and momentum only, when event is accepted (return value = true)
-// this function is more efficient, if only weighted evens are needed
-bool
-nBodyPhaseSpaceGen::generateDecayAccepted(const lorentzVector& nBody, // Lorentz vector of n-body system in lab frame
- const double maxWeight) // if positive, given value is used as maximum weight, otherwise _maxWeight
-{
- const double nBodyMass = nBody.M();
- if (_n < 2) {
- printWarn << "number of daughter particles = " << _n << " is smaller than 2. "
- << "no event generated." << endl;
- return false;
- } else if (nBodyMass < _mSum[_n - 1]) {
- printWarn << "n-body mass = " << nBodyMass << " is smaller than sum of daughter masses = "
- << _mSum[_n - 1] << ". no event generated." << endl;
- return false;
- }
- pickMasses(nBodyMass);
- calcWeight();
- if (!eventAccepted(maxWeight))
- return false;
- pickAngles();
- calcEventKinematics(nBody);
- return true;
-}
-
-
-// randomly choses the (n - 2) effective masses of the respective (i + 1)-body systems
-void
-nBodyPhaseSpaceGen::pickMasses(const double nBodyMass) // total energy of the system in its RF
-{
- _M[_n - 1] = nBodyMass;
- // create vector of sorted random values
- vector<double> r(_n - 2, 0); // (n - 2) values needed for 2- through (n - 1)-body systems
- for (unsigned int i = 0; i < (_n - 2); ++i)
- r[i] = random();
- sort(r.begin(), r.end());
- // set effective masses of (intermediate) two-body decays
- const double massInterval = nBodyMass - _mSum[_n - 1]; // kinematically allowed mass interval
- for (unsigned int i = 1; i < (_n - 1); ++i) // loop over intermediate 2- to (n - 1)-bodies
- _M[i] = _mSum[i] + r[i - 1] * massInterval; // _mSum[i] is minimum effective mass
-}
-
-
-// computes event weight (= integrand value) and breakup momenta
-// uses vector of intermediate two-body masses prepared by pickMasses()
-double
-nBodyPhaseSpaceGen::calcWeight()
-{
- for (unsigned int i = 1; i < _n; ++i) // loop over 2- to n-bodies
- _breakupMom[i] = breakupMomentum(_M[i], _M[i - 1], _m[i]);
- double momProd = 1; // product of breakup momenta
- for (unsigned int i = 1; i < _n; ++i) // loop over 2- to n-bodies
- momProd *= _breakupMom[i];
- const double massInterval = _M[_n - 1] - _mSum[_n - 1]; // kinematically allowed mass interval
- _weight = _norm * pow(massInterval, (int)_n - 2) * momProd / _M[_n - 1];
- if (_weight > _maxWeightObserved)
- _maxWeightObserved = _weight;
- if (std::isnan(_weight))
- printWarn << "weight = " << _weight << endl;
- return _weight;
-}
-
-
-// calculates complete event from the effective masses of the (i + 1)-body
-// systems, the Lorentz vector of the decaying system, and the decay angles
-// uses the break-up momenta calculated by calcWeight()
-void
-nBodyPhaseSpaceGen::calcEventKinematics(const lorentzVector& nBody) // Lorentz vector of n-body system in lab frame
-{
- // build event starting in n-body RF going down to 2-body RF
- // is more efficicient than Raubold-Lynch method, since it requitres only one rotation and boost per daughter
- lorentzVector P = nBody; // Lorentz of (i + 1)-body system in lab frame
- for (unsigned int i = _n - 1; i >= 1; --i) { // loop from n-body down to 2-body
- // construct Lorentz vector of daughter _m[i] in (i + 1)-body RF
- const double sinTheta = sqrt(1 - _cosTheta[i] * _cosTheta[i]);
- const double pT = _breakupMom[i] * sinTheta;
- lorentzVector& daughter = _daughters[i];
- daughter.SetPxPyPzE(pT * cos(_phi[i]),
- pT * sin(_phi[i]),
- _breakupMom[i] * _cosTheta[i],
- sqrt(_m[i] * _m[i] + _breakupMom[i] * _breakupMom[i]));
- // boost daughter into lab frame
- daughter.Boost(P.BoostVector());
- // calculate Lorentz vector of i-body system in lab frame
- P -= daughter;
- }
- // set last daughter
- _daughters[0] = P;
-}
-
-
-// calculates maximum weight for given n-body mass
-double
-nBodyPhaseSpaceGen::estimateMaxWeight(const double nBodyMass, // sic!
- const unsigned int nmbOfIterations) // number of generated events
-{
- double maxWeight = 0;
- for (unsigned int i = 0; i < nmbOfIterations; ++i) {
- pickMasses(nBodyMass);
- calcWeight();
- maxWeight = max(_weight, maxWeight);
- }
- return maxWeight;
-}
-
-
-ostream&
-nBodyPhaseSpaceGen::print(ostream& out) const
-{
- out << "nBodyPhaseSpaceGen parameters:" << endl
- << " number of daughter particles ............... " << _n << endl
- << " masses of the daughter particles ........... " << _m << endl
- << " sums of daughter particle masses ........... " << _mSum << endl
- << " effective masses of (i + 1)-body systems ... " << _M << endl
- << " cos(polar angle) in (i + 1)-body systems ... " << _cosTheta << endl
- << " azimuth in (i + 1)-body systems ............ " << _phi << endl
- << " breakup momenta in (i + 1)-body systems .... " << _breakupMom << endl
- << " normalization value ........................ " << _norm << endl
- << " weight of generated event .................. " << _weight << endl
- << " maximum weight used in hit-miss MC ......... " << _maxWeight << endl
- << " maximum weight since instantiation ......... " << _maxWeightObserved << endl
- << " daughter four-momenta:" << endl;
- for (unsigned int i = 0; i < _n; ++i)
- out << " daughter " << i << ": " << _daughters[i] << endl;
- return out;
-}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff