[u/mrichter/AliRoot.git] / STARLIGHT / starlight / include / nBodyPhaseSpaceGen.h

///////////////////////////////////////////////////////////////////////////
//
//    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 maxWeightArg)  // if maxWeightArg > 0, given value is used as maximum weight, otherwise _maxWeight
{
	const double max = (maxWeightArg <= 0) ? _maxWeight : maxWeightArg;
	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
Commit	Line	Data
da32329d AM	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