[u/mrichter/AliRoot.git] / STARLIGHT / starlight / src / nBodyPhaseSpaceGen.cpp

///////////////////////////////////////////////////////////////////////////
//
//    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:: 159                         $: revision of last commit
// $Author:: odjuvsla                 $: author of last commit
// $Date:: 2013-10-06 16:17:58 +0200 #$: 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;
}
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:: 159 $: revision of last commit
	24	// $Author:: odjuvsla $: author of last commit
	25	// $Date:: 2013-10-06 16:17:58 +0200 #$: date of last commit
	26	//
	27	// Description:
	28	// see nBodyPhaseSpaceGen.h
	29	//
	30	//
	31	///////////////////////////////////////////////////////////////////////////
	32
	33
	34	#include <algorithm>
	35
	36	#include "nBodyPhaseSpaceGen.h"
	37
	38
	39	using namespace std;
	40	using namespace starlightConstants;
	41
	42
	43	nBodyPhaseSpaceGen::nBodyPhaseSpaceGen()
	44	: _n (0),
	45	_norm (0),
	46	_weight (0),
	47	_maxWeightObserved(0),
	48	_maxWeight (0)
	49	{ }
	50
	51
	52	nBodyPhaseSpaceGen::~nBodyPhaseSpaceGen()
	53	{ }
	54
	55
	56	// sets decay constants and prepares internal variables
	57	bool
	58	nBodyPhaseSpaceGen::setDecay(const vector<double>& daughterMasses) // array of daughter particle masses
	59	{
	60	_n = daughterMasses.size();
	61	if (_n < 2) {
	62	printWarn << "number of daughters = " << _n << " does not make sense." << endl;
	63	return false;
	64	}
65	// copy daughter masses
66	_m.clear();
67	_m = daughterMasses;
68	// prepare effective mass vector
69	_M.clear();
70	_M.resize(_n, 0);
71	_M[0] = _m[0];
72	// prepare angle vectors
73	_cosTheta.clear();
74	_cosTheta.resize(_n, 0);
75	_phi.clear();
76	_phi.resize(_n, 0);
77	// calculate daughter mass sums
78	_mSum.clear();
79	_mSum.resize(_n, 0);
80	_mSum[0] = _m[0];
81	for (unsigned int i = 1; i < _n; ++i)
82	_mSum[i] = _mSum[i - 1] + _m[i];
83	// prepare breakup momentum vector
84	_breakupMom.clear();
85	_breakupMom.resize(_n, 0);
86	// prepare vector for daughter Lorentz vectors
87	_daughters.clear();
88	_daughters.resize(_n, lorentzVector(0, 0, 0, 0));
89	// calculate normalization
90	_norm = 1 / (2 * pow(twoPi, 2 * (int)_n - 3) * factorial(_n - 2));
91	resetMaxWeightObserved();
92	return true;
93	}
94
95
96	// set decay constants and prepare internal variables
97	bool
45d54d9a	98	nBodyPhaseSpaceGen::setDecay(const unsigned int nmbOfDaughters, // number of daughter particles
45d54d9a	99	const double* daughterMasses) // array of daughter particle masses
da32329d AM	100	{
da32329d AM	101	vector <double> m;
45d54d9a	102	m.resize(nmbOfDaughters, 0);
45d54d9a	103	for (unsigned int i = 0; i < nmbOfDaughters; ++i)
da32329d AM	104	m[i] = daughterMasses[i];
	105	return setDecay(m);
	106	}
	107
	108
	109	// generates event with certain n-body mass and momentum and returns event weigth
	110	// general purpose function
	111	double
	112	nBodyPhaseSpaceGen::generateDecay(const lorentzVector& nBody) // Lorentz vector of n-body system in lab frame
	113	{
	114	const double nBodyMass = nBody.M();
	115	if (_n < 2) {
	116	printWarn << "number of daughter particles = " << _n << " is smaller than 2. "
	117	<< "weight is set to 0." << endl;
	118	_weight = 0;
	119	} else if (nBodyMass < _mSum[_n - 1]) {
	120	printWarn << "n-body mass = " << nBodyMass << " is smaller than sum of daughter masses = "
	121	<< _mSum[_n - 1] << ". weight is set to 0." << endl;
	122	_weight = 0;
	123	} else {
	124	pickMasses(nBodyMass);
	125	calcWeight();
	126	pickAngles();
	127	calcEventKinematics(nBody);
	128	}
	129	return _weight;
	130	}
	131
	132
	133	// generates full event with certain n-body mass and momentum only, when event is accepted (return value = true)
	134	// this function is more efficient, if only weighted evens are needed
	135	bool
45d54d9a	136	nBodyPhaseSpaceGen::generateDecayAccepted(const lorentzVector& nBody, // Lorentz vector of n-body system in lab frame
45d54d9a	137	const double maxWeight) // if positive, given value is used as maximum weight, otherwise _maxWeight
da32329d AM	138	{
	139	const double nBodyMass = nBody.M();
	140	if (_n < 2) {
	141	printWarn << "number of daughter particles = " << _n << " is smaller than 2. "
	142	<< "no event generated." << endl;
	143	return false;
	144	} else if (nBodyMass < _mSum[_n - 1]) {
	145	printWarn << "n-body mass = " << nBodyMass << " is smaller than sum of daughter masses = "
	146	<< _mSum[_n - 1] << ". no event generated." << endl;
	147	return false;
	148	}
	149	pickMasses(nBodyMass);
	150	calcWeight();
45d54d9a	151	if (!eventAccepted(maxWeight))
da32329d AM	152	return false;
	153	pickAngles();
	154	calcEventKinematics(nBody);
	155	return true;
	156	}
	157
	158
	159	// randomly choses the (n - 2) effective masses of the respective (i + 1)-body systems
	160	void
	161	nBodyPhaseSpaceGen::pickMasses(const double nBodyMass) // total energy of the system in its RF
	162	{
	163	_M[_n - 1] = nBodyMass;
	164	// create vector of sorted random values
	165	vector<double> r(_n - 2, 0); // (n - 2) values needed for 2- through (n - 1)-body systems
	166	for (unsigned int i = 0; i < (_n - 2); ++i)
	167	r[i] = random();
	168	sort(r.begin(), r.end());
	169	// set effective masses of (intermediate) two-body decays
	170	const double massInterval = nBodyMass - _mSum[_n - 1]; // kinematically allowed mass interval
	171	for (unsigned int i = 1; i < (_n - 1); ++i) // loop over intermediate 2- to (n - 1)-bodies
	172	_M[i] = _mSum[i] + r[i - 1] * massInterval; // _mSum[i] is minimum effective mass
	173	}
	174
	175
	176	// computes event weight (= integrand value) and breakup momenta
	177	// uses vector of intermediate two-body masses prepared by pickMasses()
	178	double
	179	nBodyPhaseSpaceGen::calcWeight()
	180	{
	181	for (unsigned int i = 1; i < _n; ++i) // loop over 2- to n-bodies
	182	_breakupMom[i] = breakupMomentum(_M[i], _M[i - 1], _m[i]);
	183	double momProd = 1; // product of breakup momenta
	184	for (unsigned int i = 1; i < _n; ++i) // loop over 2- to n-bodies
	185	momProd *= _breakupMom[i];
	186	const double massInterval = _M[_n - 1] - _mSum[_n - 1]; // kinematically allowed mass interval
	187	_weight = _norm * pow(massInterval, (int)_n - 2) * momProd / _M[_n - 1];
	188	if (_weight > _maxWeightObserved)
	189	_maxWeightObserved = _weight;
	190	if (std::isnan(_weight))
	191	printWarn << "weight = " << _weight << endl;
	192	return _weight;
	193	}
	194
	195
	196	// calculates complete event from the effective masses of the (i + 1)-body
	197	// systems, the Lorentz vector of the decaying system, and the decay angles
	198	// uses the break-up momenta calculated by calcWeight()
	199	void
	200	nBodyPhaseSpaceGen::calcEventKinematics(const lorentzVector& nBody) // Lorentz vector of n-body system in lab frame
	201	{
	202	// build event starting in n-body RF going down to 2-body RF
	203	// is more efficicient than Raubold-Lynch method, since it requitres only one rotation and boost per daughter
	204	lorentzVector P = nBody; // Lorentz of (i + 1)-body system in lab frame
	205	for (unsigned int i = _n - 1; i >= 1; --i) { // loop from n-body down to 2-body
	206	// construct Lorentz vector of daughter _m[i] in (i + 1)-body RF
45d54d9a	207	const double sinTheta = sqrt(1 - _cosTheta[i] * _cosTheta[i]);
	208	const double pT = _breakupMom[i] * sinTheta;
	209	lorentzVector& daughter = _daughters[i];
	210	daughter.SetPxPyPzE(pT * cos(_phi[i]),
	211	pT * sin(_phi[i]),
	212	_breakupMom[i] * _cosTheta[i],
	213	sqrt(_m[i] * _m[i] + _breakupMom[i] * _breakupMom[i]));
da32329d	214	// boost daughter into lab frame
45d54d9a	215	daughter.Boost(P.BoostVector());
da32329d	216	// calculate Lorentz vector of i-body system in lab frame
45d54d9a	217	P -= daughter;
da32329d AM	218	}
	219	// set last daughter
	220	_daughters[0] = P;
	221	}
	222
	223
	224	// calculates maximum weight for given n-body mass
	225	double
	226	nBodyPhaseSpaceGen::estimateMaxWeight(const double nBodyMass, // sic!
	227	const unsigned int nmbOfIterations) // number of generated events
	228	{
45d54d9a	229	double maxWeight = 0;
da32329d AM	230	for (unsigned int i = 0; i < nmbOfIterations; ++i) {
	231	pickMasses(nBodyMass);
	232	calcWeight();
45d54d9a	233	maxWeight = max(_weight, maxWeight);
da32329d	234	}
45d54d9a	235	return maxWeight;
da32329d AM	236	}
	237
	238
	239	ostream&
	240	nBodyPhaseSpaceGen::print(ostream& out) const
	241	{
	242	out << "nBodyPhaseSpaceGen parameters:" << endl
	243	<< " number of daughter particles ............... " << _n << endl
	244	<< " masses of the daughter particles ........... " << _m << endl
	245	<< " sums of daughter particle masses ........... " << _mSum << endl
	246	<< " effective masses of (i + 1)-body systems ... " << _M << endl
	247	<< " cos(polar angle) in (i + 1)-body systems ... " << _cosTheta << endl
	248	<< " azimuth in (i + 1)-body systems ............ " << _phi << endl
	249	<< " breakup momenta in (i + 1)-body systems .... " << _breakupMom << endl
	250	<< " normalization value ........................ " << _norm << endl
	251	<< " weight of generated event .................. " << _weight << endl
	252	<< " maximum weight used in hit-miss MC ......... " << _maxWeight << endl
	253	<< " maximum weight since instantiation ......... " << _maxWeightObserved << endl
	254	<< " daughter four-momenta:" << endl;
	255	for (unsigned int i = 0; i < _n; ++i)
	256	out << " daughter " << i << ": " << _daughters[i] << endl;
	257	return out;
	258	}