--- /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:: 157 $: revision of last commit
+// $Author:: odjuvsla $: author of last commit
+// $Date:: 2013-10-06 16:17:54 +0200 #$: date of last commit
+//
+// Description:
+// calculates n-body phase space (constant matrix element) using various algorithms
+//
+// the n-body decay is split up into (n - 2) successive 2-body decays
+// each 2-body decay is considered in its own center-of-mass frame thereby
+// separating the mass from the (trivial) angular dependence
+//
+// the event is boosted into the same frame in which the n-body system is
+// given
+//
+// based on:
+// GENBOD (CERNLIB W515), see F. James, "Monte Carlo Phase Space", CERN 68-15 (1968)
+// NUPHAZ, see M. M. Block, "Monte Carlo phase space evaluation", Comp. Phys. Commun. 69, 459 (1992)
+// S. U. Chung, "Spin Formalism", CERN Yellow Report
+// S. U. Chung et. al., "Diffractive Dissociation for COMPASS"
+//
+// index convention:
+// - all vectors have the same size (= number of decay daughters)
+// - index i corresponds to the respective value in the (i + 1)-body system: effective mass M, break-up momentum, angles
+// - thus some vector elements are not used like breakupMom[0], theta[0], phi[0], ...
+// this overhead is negligible compared to the ease of notation
+//
+// the following graph illustrates how the n-body decay is decomposed into a sequence of two-body decays
+//
+// n-body ... 3-body 2-body single daughter
+//
+// m[n - 1] m[2] m[1]
+// ^ ^ ^
+// | | |
+// | | |
+// M[n - 1] --> ... --> M[2] --> M[1] --> M [0] = m[0]
+// theta[n - 1] ... theta[2] theta[1] theta[0] = 0 (not used)
+// phi [n - 1] ... phi [2] phi [1] phi [0] = 0 (not used)
+// mSum [n - 1] ... mSum [2] mSum [1] mSum [0] = m[0]
+// = sum_0^(n - 1) m[i] = m[2] + m[1] + m[0] = m[1] + m[0]
+// breakUpMom[n - 1] ... breakUpMom[2] breakUpMom[1] breakUpMom[0] = 0 (not used)
+// = q(M[n - 1], m[n - 1], M[n - 2]) = q(M[2], m[2], M[1]) = q(M[1], m[1], m[0])
+//
+//
+///////////////////////////////////////////////////////////////////////////
+
+
+#ifndef NBODYPHASESPACEGEN_H
+#define NBODYPHASESPACEGEN_H
+
+
+#include <iostream>
+#include <vector>
+
+#include "reportingUtils.h"
+#include "lorentzvector.h"
+#include "randomgenerator.h"
+#include "starlightconstants.h"
+
+
+// small helper functions
+// calculates factorial
+inline
+unsigned int
+factorial(const unsigned int n)
+{
+ unsigned int fac = 1;
+ for (unsigned int i = 1; i <= n; ++i)
+ fac *= i;
+ return fac;
+}
+
+
+// computes breakup momentum of 2-body decay
+inline
+double
+breakupMomentum(const double M, // mass of mother particle
+ const double m1, // mass of daughter particle 1
+ const double m2) // mass of daughter particle 2
+{
+ if (M < m1 + m2)
+ return 0;
+ return sqrt((M - m1 - m2) * (M + m1 + m2) * (M - m1 + m2) * (M + m1 - m2)) / (2 * M);
+}
+
+
+class nBodyPhaseSpaceGen {
+
+public:
+
+ nBodyPhaseSpaceGen();
+ virtual ~nBodyPhaseSpaceGen();
+
+ // generator setup
+ /// sets decay constants and prepares internal variables
+ bool setDecay(const std::vector<double>& daughterMasses); // daughter particle masses
+ bool setDecay(const unsigned int nmbOfDaughters, // number of daughter particles
+ const double* daughterMasses); // array of daughter particle masses
+
+ // random generator
+ double random () { return randyInstance.Rndom(); } ///< returns number from internal random generator
+
+ // high-level generator interface
+ /// generates full event with certain n-body mass and momentum and returns event weight
+ double generateDecay (const lorentzVector& nBody); // Lorentz vector of n-body system in lab frame
+ /// \brief 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 events are needed
+ bool generateDecayAccepted(const lorentzVector& nBody, // Lorentz vector of n-body system in lab frame
+ const double maxWeight = 0); // if positive, given value is used as maximum weight, otherwise _maxWeight
+
+ void setMaxWeight (const double maxWeight) { _maxWeight = maxWeight; } ///< sets maximum weight used for hit-miss MC
+ double maxWeight () const { return _maxWeight; } ///< returns maximum weight used for hit-miss MC
+ double normalization () const { return _norm; } ///< returns normalization used in weight calculation
+ double eventWeight () const { return _weight; } ///< returns weight of generated event
+ double maxWeightObserved () const { return _maxWeightObserved; } ///< returns maximum observed weight since instantiation
+ void resetMaxWeightObserved() { _maxWeightObserved = 0; } ///< sets maximum observed weight back to zero
+
+ /// estimates maximum weight for given n-body mass
+ double estimateMaxWeight(const double nBodyMass, // sic!
+ const unsigned int nmbOfIterations = 10000); // number of generated events
+
+ /// \brief applies event weight in form of hit-miss MC
+ /// assumes that event weight has been already calculated by calcWeight()
+ /// if maxWeight > 0 value is used as maximum weight, otherwise _maxWeight value is used
+ inline bool eventAccepted(const double maxWeight = 0);
+
+ //----------------------------------------------------------------------------
+ // trivial accessors
+ const lorentzVector& daughter (const int index) const { return _daughters[index]; } ///< returns Lorentz vector of daughter at index
+ const std::vector<lorentzVector>& daughters () const { return _daughters; } ///< returns Lorentz vectors of all daughters
+ unsigned int nmbOfDaughters () const { return _n; } ///< returns number of daughters
+ double daughterMass (const int index) const { return _m[index]; } ///< returns invariant mass of daughter at index
+ double intermediateMass(const int index) const { return _M[index]; } ///< returns intermediate mass of (index + 1)-body system
+ double breakupMom (const int index) const { return _breakupMom[index]; } ///< returns breakup momentum in (index + 1)-body RF
+ double cosTheta (const int index) const { return _cosTheta[index]; } ///< returns polar angle in (index + 1)-body RF
+ double phi (const int index) const { return _phi[index]; } ///< returns azimuth in (index + 1)-body RF
+
+
+ std::ostream& print(std::ostream& out = std::cout) const; ///< prints generator status
+ friend std::ostream& operator << (std::ostream& out,
+ const nBodyPhaseSpaceGen& gen)
+ { return gen.print(out); }
+
+private:
+
+ //----------------------------------------------------------------------------
+ // low-level generator interface
+ /// randomly choses the (n - 2) effective masses of the respective (i + 1)-body systems
+ void pickMasses(const double nBodyMass); // total energy of n-body system in its RF
+
+ /// \brief computes event weight and breakup momenta
+ /// operates on vector of intermediate two-body masses prepared by pickMasses()
+ double calcWeight();
+
+ /// randomly choses the (n - 1) polar and (n - 1) azimuthal angles in the respective (i + 1)-body RFs
+ inline void pickAngles();
+
+ /// \brief calculates full event kinematics from the effective masses of the (i + 1)-body systems and the Lorentz vector of the decaying system
+ /// uses the break-up momenta calculated by calcWeight() and angles from pickAngles()
+ void calcEventKinematics(const lorentzVector& nBody); // Lorentz vector of n-body system in lab frame
+
+ // external parameters
+ std::vector<double> _m; ///< masses of daughter particles
+
+ // internal variables
+ unsigned int _n; ///< number of daughter particles
+ std::vector<double> _M; ///< effective masses of (i + 1)-body systems
+ std::vector<double> _cosTheta; ///< cosine of polar angle of the 2-body decay of the (i + 1)-body system
+ std::vector<double> _phi; ///< azimuthal angle of the 2-body decay of the (i + 1)-body system
+ std::vector<double> _mSum; ///< sums of daughter particle masses
+ std::vector<double> _breakupMom; ///< breakup momenta for the two-body decays: (i + 1)-body --> daughter_(i + 1) + i-body
+ std::vector<lorentzVector> _daughters; ///< Lorentz vectors of the daughter particles
+ double _norm; ///< normalization value
+ double _weight; ///< phase space weight of generated event
+ double _maxWeightObserved; ///< maximum event weight calculated processing the input data
+ double _maxWeight; ///< maximum weight used to weight events in hit-miss MC
+
+
+};
+
+
+inline
+void
+nBodyPhaseSpaceGen::pickAngles()
+{
+ for (unsigned int i = 1; i < _n; ++i) { // loop over 2- to n-bodies
+ _cosTheta[i] = 2 * random() - 1; // range [-1, 1]
+ _phi[i] = starlightConstants::twoPi * random(); // range [ 0, 2 pi]
+ }
+}
+
+
+inline
+bool
+nBodyPhaseSpaceGen::eventAccepted(const double maxWeight) // if maxWeight > 0, given value is used as maximum weight, otherwise _maxWeight
+{
+ const double max = (maxWeight <= 0) ? _maxWeight : maxWeight;
+ if (max <= 0) {
+ printWarn << "maximum weight = " << max << " does not make sense. rejecting event." << std::endl;
+ return false;
+ }
+ if ((_weight / max) > random())
+ return true;
+ return false;
+}
+
+
+#endif // NBODYPHASESPACEGEN_H
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff