[u/mrichter/AliRoot.git] / EVGEN / AliGenDeuteron.cxx

/**************************************************************************
 * Copyright(c) 2009-2010, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

//////////////////////////////////////////////////////////////////////////
// Afterburner to generate (anti)deuterons simulating coalescence of
// (anti)nucleons as in A. J. Baltz et al., Phys. lett B 325(1994)7.
// If the nucleon generator does not provide a spatial description of 
// the source the afterburner can provide one.
//
// There two options for the source: a thermal picture where all nucleons
// are placed randomly and homogeneously in an spherical volume, or
// an expansion one where they are projected onto its surface.
//
// An (anti)deuteron will form if there is a pair of (anti)proton-(anti)neutron
// with momentum difference less than ~ 300MeV and relative distance less than
// ~ 2.1fm. Only 3/4 of these clusters are accepted due to the deuteron spin.
//
// When there are more than one pair fulfilling the coalescence conditions,
// the selected pair can be the one with the first partner, with
// the lowest relative momentum, the lowest relative distance or both.
//////////////////////////////////////////////////////////////////////////

// Author: Eulogio Serradilla <eulogio.serradilla@ciemat.es>,
//         Arturo Menchaca <menchaca@fisica.unam.mx>
//

#include "Riostream.h"
#include "TParticle.h"
#include "AliStack.h"
#include "AliRun.h"
#include "TMath.h"
#include "TMCProcess.h"
#include "TList.h"
#include "TVector3.h"
#include "AliMC.h"
#include "TArrayF.h"
#include "AliCollisionGeometry.h"
#include "AliGenCocktailEventHeader.h"
#include "AliGenCocktailEntry.h"
#include "AliGenCocktailAfterBurner.h"
#include "AliGenDeuteron.h"
#include "AliLog.h"

ClassImp(AliGenDeuteron)

AliGenDeuteron::AliGenDeuteron(Int_t sign, Double_t pmax, Double_t rmax, Int_t cluster)
 :fSign(1)
 ,fPmax(pmax)
 ,fRmax(rmax)
 ,fSpinProb(1)
 ,fClusterType(cluster)
 ,fModel(0)
 ,fTimeLength(2.5)
 ,fB(0)
 ,fR(0)
 ,fPsiR(0)
 ,fCurStack(0)
{
//
// constructor
//
	fSign = sign > 0 ? 1:-1;
	
}

AliGenDeuteron::~AliGenDeuteron()
{
//
// Default destructor
//
}

void AliGenDeuteron::Init()
{
//
// Standard AliGenerator initializer
//
}

void AliGenDeuteron::Generate()
{
//
// Look for coalescence of (anti)nucleons in the nucleon source
// provided by the generator cocktail
//
	AliInfo(Form("sign: %d, Pmax: %g GeV/c, Rmax: %g fm, cluster: %d",fSign, fPmax, fRmax, fClusterType));
	if(fModel!=kNone) AliInfo(Form("model: %d, time: %g fm/c", fModel, fTimeLength));
	
	AliGenCocktailAfterBurner* gener = (AliGenCocktailAfterBurner*)gAlice->GetMCApp()->Generator();
	
	// find nuclear radius, only for first generator and projectile=target
	Bool_t collisionGeometry=0;
	if(fModel != kNone && gener->FirstGenerator()->Generator()->ProvidesCollisionGeometry())
	{
		TString name;
		Int_t ap, zp, at, zt;
		gener->FirstGenerator()->Generator()->GetProjectile(name,ap,zp);
		gener->FirstGenerator()->Generator()->GetTarget(name,at,zt);
		if(ap != at) AliWarning("projectile and target have different size");
		fR = 1.29*TMath::Power(at, 1./3.);
		collisionGeometry = 1;
	}
	
	for(Int_t ns = 0; ns < gener->GetNumberOfEvents(); ++ns)
	{
		gener->SetActiveEventNumber(ns);
		
		if(fModel != kNone && collisionGeometry)
		{
			fPsiR = (gener->GetCollisionGeometry(ns))->ReactionPlaneAngle();
			fB = (gener->GetCollisionGeometry(ns))->ImpactParameter();
			
			if(fB >= 2.*fR) continue; // no collision
		}
		
		fCurStack = gener->GetStack(ns);
		if(!fCurStack)
		{
			AliWarning("no event stack");
			return;
		}
		
		TList* protons = new TList();
		protons->SetOwner(kFALSE);
		TList* neutrons = new TList();
		neutrons->SetOwner(kFALSE);
		
		// particles produced by the generator
		for (Int_t i=0; i < fCurStack->GetNprimary(); ++i)
		{
			TParticle* iParticle = fCurStack->Particle(i);
			
			if(iParticle->GetStatusCode() != 1) continue;
			
			Int_t pdgCode = iParticle->GetPdgCode();
			if(pdgCode == fSign*2212)// (anti)proton
			{
				FixProductionVertex(iParticle);
				protons->Add(iParticle);
			}
			else if(pdgCode == fSign*2112) // (anti)neutron
			{
				FixProductionVertex(iParticle);
				neutrons->Add(iParticle);
			}
		}
		
		// look for clusters
		if(fClusterType==kFirstPartner)
		{
			FirstPartner(protons, neutrons);
		}
		else
		{
			WeightMatrix(protons, neutrons);
		}
		
		protons->Clear("nodelete");
		neutrons->Clear("nodelete");
		
		delete protons;
		delete neutrons;
	}
	
	AliInfo("DONE");
}

Double_t AliGenDeuteron::GetCoalescenceProbability(const TParticle* nucleon1, const TParticle* nucleon2) const
{
//
// Coalescence conditions as in
// A. J. Baltz et al., Phys. lett B 325(1994)7
//
// returns < 0 if coalescence is not possible
// otherwise returns a coalescence probability
//
	const Double_t kProtonMass  = 0.938272013;
	const Double_t kNeutronMass = 0.939565378;
	
	TVector3 v1(nucleon1->Vx(), nucleon1->Vy(), nucleon1->Vz());
	TVector3 p1(nucleon1->Px(), nucleon1->Py(), nucleon1->Pz());
	
	TVector3 v2(nucleon2->Vx(), nucleon2->Vy(), nucleon2->Vz());
	TVector3 p2(nucleon2->Px(), nucleon2->Py(), nucleon2->Pz());
	
	Double_t deltaP = 2.*this->GetPcm(p1, kProtonMass, p2, kNeutronMass); // relative momentum in CM frame
	if( deltaP >= fPmax) return -1.;
	
	Double_t deltaR = (v2-v1).Mag();       // relative distance (cm)
	if(deltaR >= fRmax*1.e-13) return -1.;
	
	if(Rndm() > fSpinProb)    return -1.;  // spin
	
	if(fClusterType == kLowestMomentum) return 1. - deltaP/fPmax;
	if(fClusterType == kLowestDistance) return 1. - 1.e+13*deltaR/fRmax;
	
	return 1. - 1.e+13*(deltaP*deltaR)/(fRmax*fPmax);
}

void AliGenDeuteron::FirstPartner(const TList* protons, TList* neutrons)
{
//
// Clusters are made with the first nucleon pair that fulfill
// the coalescence conditions, starting with the protons
//
	TIter p_next(protons);
	while(TParticle* n0 = (TParticle*) p_next())
	{
		TParticle* partner = 0;
		TIter n_next(neutrons);
		while(TParticle* n1 = (TParticle*) n_next() )
		{
			if(GetCoalescenceProbability(n0, n1) < 0 ) continue; // with next neutron
			partner = n1;
			break;
		}
		
		if(partner == 0) continue; // with next proton
		
		PushDeuteron(n0, partner);
		
		// Remove from the list for the next iteration
		neutrons->Remove(partner);
	}
}

void AliGenDeuteron::WeightMatrix(const TList* protons, const TList* neutrons)
{
//
// Build all possible nucleon pairs with their own probability
// and select only those with the highest coalescence probability
//
	Int_t nMaxPairs = protons->GetSize()*neutrons->GetSize();
	
	TParticle** cProton = new TParticle*[nMaxPairs];
	TParticle** cNeutron = new TParticle*[nMaxPairs];
	Double_t*  cWeight = new Double_t[nMaxPairs];
	
	// build all pairs with probability > 0
	Int_t cIdx = -1;
	TIter p_next(protons);
	while(TParticle* n1 = (TParticle*) p_next())
	{
		TIter n_next(neutrons);
		while(TParticle* n2 = (TParticle*) n_next() )
		{
			Double_t weight = this->GetCoalescenceProbability(n1,n2);
			if(weight < 0) continue;
			++cIdx;
			cProton[cIdx]  = n1;
			cNeutron[cIdx] = n2;
			cWeight[cIdx]  = weight;
		}
		n_next.Reset();
	}
	p_next.Reset();
	
	Int_t nPairs = cIdx + 1;
	
	// find the interacting pairs:
	// remove repeated pairs and select only
	// the pair with the highest coalescence probability
	
	Int_t nMaxIntPair = TMath::Min(protons->GetSize(), neutrons->GetSize());
	
	TParticle** iProton = new TParticle*[nMaxIntPair];
	TParticle** iNeutron = new TParticle*[nMaxIntPair];
	Double_t* iWeight = new Double_t[nMaxIntPair];
	
	Int_t iIdx = -1;
	while(kTRUE)
	{
		Int_t j = -1;
		Double_t wMax = 0;
		for(Int_t i=0; i < nPairs; ++i)
		{
			if(cWeight[i] > wMax)
			{
				wMax=cWeight[i];
				j = i;
			}
		}
		
		if(j == -1 ) break; // end
		
		// Save the interacting pair
		++iIdx;
		iProton[iIdx]  = cProton[j];
		iNeutron[iIdx] = cNeutron[j];
		iWeight[iIdx]  = cWeight[j];
		
		// invalidate all combinations with these pairs for the next iteration
		for(Int_t i=0; i < nPairs; ++i)
		{
			if(cProton[i] == iProton[iIdx]) cWeight[i] = -1.;
			if(cNeutron[i] == iNeutron[iIdx]) cWeight[i] = -1.;
		}
	}
	
	Int_t nIntPairs = iIdx + 1;
	
	delete[] cProton;
	delete[] cNeutron;
	delete[] cWeight;
	
	// Add the (anti)deuterons to the current event stack
	for(Int_t i=0; i<nIntPairs; ++i)
	{
		TParticle* n1 = iProton[i];
		TParticle* n2 = iNeutron[i];
		PushDeuteron(n1,n2);
	}
	
	delete[] iProton;
	delete[] iNeutron;
	delete[] iWeight;
}

void AliGenDeuteron::PushDeuteron(TParticle* parent1, TParticle* parent2)
{
//
// Create an (anti)deuteron from parent1 and parent2,
// add to the current stack and set kDoneBit for the parents
//
	const Double_t kDeuteronMass = 1.87561282;
	const Int_t kDeuteronPdg = 1000010020;
	
	// momentum
	TVector3 p1(parent1->Px(), parent1->Py(), parent1->Pz());
	TVector3 p2(parent2->Px(), parent2->Py(), parent2->Pz());
	TVector3 pN = p1+p2;
	
	// production vertex same as the parent1's
	TVector3 vN(parent1->Vx(), parent1->Vy(), parent1->Vz());
	
	// E^2 = p^2 + m^2
	Double_t energy = TMath::Sqrt(pN.Mag2() + kDeuteronMass*kDeuteronMass);
	
	Int_t ntrk = 0;
	Double_t weight = 1;
	Int_t is = 1; // final state particle
	
	// Add a new (anti)deuteron to current event stack
	fCurStack->PushTrack(1, -1, fSign*kDeuteronPdg,
	                 pN.X(), pN.Y(), pN.Z(), energy,
	                 vN.X(), vN.Y(), vN.Z(), parent1->T(),
	                 0., 0., 0., kPNCapture, ntrk, weight, is);
	
	// change the status code of the parents
	parent1->SetStatusCode(kCluster);
	parent2->SetStatusCode(kCluster);
	
	// Set kDoneBit for the parents
	parent1->SetBit(kDoneBit);
	parent2->SetBit(kDoneBit);
}

void AliGenDeuteron::FixProductionVertex(TParticle* i)
{
//
// Replace current generator nucleon spatial distribution
// with a custom distribution according to the selected model
//
	if(fModel == kNone || fModel > kExpansion) return;
	
	// semi-axis from collision geometry (fm)
	Double_t a = fTimeLength + TMath::Sqrt(fR*fR - fB*fB/4.);
	Double_t b = fTimeLength + fR - fB/2.;
	Double_t c = fTimeLength;
	
	Double_t xx = 0;
	Double_t yy = 0;
	Double_t zz = 0;
	
	if(fModel == kThermal)
	{
		// uniformly ditributed in the volume on an ellipsoid
		// random (r,theta,phi) unit sphere
		Double_t r = TMath::Power(Rndm(),1./3.);
		Double_t theta = TMath::ACos(2.*Rndm()-1.);
		Double_t phi = 2.*TMath::Pi()*Rndm();
		
		// transform coordenates
		xx = a*r*TMath::Sin(theta)*TMath::Cos(phi);
		yy = b*r*TMath::Sin(theta)*TMath::Sin(phi);
		zz = c*r*TMath::Cos(theta);
	}
	else if(fModel == kExpansion)
	{
		// project into the surface of an ellipsoid
		xx = a*TMath::Sin(i->Theta())*TMath::Cos(i->Phi());
		yy = b*TMath::Sin(i->Theta())*TMath::Sin(i->Phi());
		zz = c*TMath::Cos(i->Theta());
	}
	
	// rotate by the reaction plane angle
	Double_t x = xx*TMath::Cos(fPsiR)+yy*TMath::Sin(fPsiR);
	Double_t y = -xx*TMath::Sin(fPsiR)+yy*TMath::Cos(fPsiR);
	Double_t z = zz;
	
	// translate by the production vertex (cm)
	i->SetProductionVertex(i->Vx() + 1.e-13*x, i->Vy() + 1.e-13*y, i->Vz() + 1.e-13*z, i->T());
}

Double_t AliGenDeuteron::GetS(Double_t p1x, Double_t p1y, Double_t p1z, Double_t m1, Double_t p2x, Double_t p2y, Double_t p2z, Double_t m2) const
{
//
// square of the center of mass energy
//
	Double_t E1 = TMath::Sqrt( p1x*p1x + p1y*p1y + p1z*p1z + m1*m1);
	Double_t E2 = TMath::Sqrt( p2x*p2x + p2y*p2y + p2z*p2z + m2*m2);
	
	return (E1+E2)*(E1+E2) - ((p1x+p2x)*(p1x+p2x) + (p1y+p2y)*(p1y+p2y) + (p1z+p2z)*(p1z+p2z));
}

Double_t AliGenDeuteron::GetPcm(Double_t p1x, Double_t p1y, Double_t p1z, Double_t m1, Double_t p2x, Double_t p2y, Double_t p2z, Double_t m2) const
{
//
// momentum in the CM frame
//
	Double_t s = this->GetS(p1x, p1y, p1z, m1, p2x, p2y, p2z, m2);
	
	return TMath::Sqrt( (s-(m1-m2)*(m1-m2))*(s-(m1+m2)*(m1+m2)) )/(2.*TMath::Sqrt(s));
}

Double_t AliGenDeuteron::GetPcm(const TVector3& p1, Double_t m1, const TVector3& p2, Double_t m2) const
{
//
// momentum in the CM frame
//
	return this->GetPcm(p1.Px(),p1.Py(),p1.Pz(),m1,p2.Px(),p2.Py(),p2.Pz(),m2);
}
Commit	Line	Data
be66cf70	1	/**************************************************************************
	2	* Copyright(c) 2009-2010, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	//////////////////////////////////////////////////////////////////////////
	17	// Afterburner to generate (anti)deuterons simulating coalescence of
	18	// (anti)nucleons as in A. J. Baltz et al., Phys. lett B 325(1994)7.
	19	// If the nucleon generator does not provide a spatial description of
	20	// the source the afterburner can provide one.
	21	//
	22	// There two options for the source: a thermal picture where all nucleons
	23	// are placed randomly and homogeneously in an spherical volume, or
	24	// an expansion one where they are projected onto its surface.
	25	//
	26	// An (anti)deuteron will form if there is a pair of (anti)proton-(anti)neutron
	27	// with momentum difference less than ~ 300MeV and relative distance less than
	28	// ~ 2.1fm. Only 3/4 of these clusters are accepted due to the deuteron spin.
	29	//
	30	// When there are more than one pair fulfilling the coalescence conditions,
	31	// the selected pair can be the one with the first partner, with
	32	// the lowest relative momentum, the lowest relative distance or both.
	33	//////////////////////////////////////////////////////////////////////////
	34
	35	// Author: Eulogio Serradilla <eulogio.serradilla@ciemat.es>,
	36	// Arturo Menchaca <menchaca@fisica.unam.mx>
	37	//
	38
	39	#include "Riostream.h"
	40	#include "TParticle.h"
be66cf70	41	#include "AliStack.h"
95145d47	42	#include "AliRun.h"
be66cf70	43	#include "TMath.h"
	44	#include "TMCProcess.h"
	45	#include "TList.h"
	46	#include "TVector3.h"
	47	#include "AliMC.h"
	48	#include "TArrayF.h"
	49	#include "AliCollisionGeometry.h"
	50	#include "AliGenCocktailEventHeader.h"
	51	#include "AliGenCocktailEntry.h"
	52	#include "AliGenCocktailAfterBurner.h"
	53	#include "AliGenDeuteron.h"
c93255fe	54	#include "AliLog.h"
be66cf70	55
	56	ClassImp(AliGenDeuteron)
	57
	58	AliGenDeuteron::AliGenDeuteron(Int_t sign, Double_t pmax, Double_t rmax, Int_t cluster)
	59	:fSign(1)
	60	,fPmax(pmax)
	61	,fRmax(rmax)
95145d47	62	,fSpinProb(1)
be66cf70	63	,fClusterType(cluster)
	64	,fModel(0)
	65	,fTimeLength(2.5)
	66	,fB(0)
	67	,fR(0)
	68	,fPsiR(0)
	69	,fCurStack(0)
be66cf70	70	{
	71	//
	72	// constructor
	73	//
	74	fSign = sign > 0 ? 1:-1;
	75
	76	}
	77
	78	AliGenDeuteron::~AliGenDeuteron()
	79	{
	80	//
	81	// Default destructor
	82	//
	83	}
	84
	85	void AliGenDeuteron::Init()
	86	{
	87	//
	88	// Standard AliGenerator initializer
	89	//
	90	}
	91
	92	void AliGenDeuteron::Generate()
	93	{
	94	//
	95	// Look for coalescence of (anti)nucleons in the nucleon source
	96	// provided by the generator cocktail
	97	//
	98	AliInfo(Form("sign: %d, Pmax: %g GeV/c, Rmax: %g fm, cluster: %d",fSign, fPmax, fRmax, fClusterType));
	99	if(fModel!=kNone) AliInfo(Form("model: %d, time: %g fm/c", fModel, fTimeLength));
	100
	101	AliGenCocktailAfterBurner* gener = (AliGenCocktailAfterBurner*)gAlice->GetMCApp()->Generator();
	102
	103	// find nuclear radius, only for first generator and projectile=target
	104	Bool_t collisionGeometry=0;
	105	if(fModel != kNone && gener->FirstGenerator()->Generator()->ProvidesCollisionGeometry())
	106	{
	107	TString name;
	108	Int_t ap, zp, at, zt;
	109	gener->FirstGenerator()->Generator()->GetProjectile(name,ap,zp);
	110	gener->FirstGenerator()->Generator()->GetTarget(name,at,zt);
	111	if(ap != at) AliWarning("projectile and target have different size");
	112	fR = 1.29*TMath::Power(at, 1./3.);
	113	collisionGeometry = 1;
	114	}
	115
	116	for(Int_t ns = 0; ns < gener->GetNumberOfEvents(); ++ns)
	117	{
	118	gener->SetActiveEventNumber(ns);
	119
	120	if(fModel != kNone && collisionGeometry)
	121	{
	122	fPsiR = (gener->GetCollisionGeometry(ns))->ReactionPlaneAngle();
	123	fB = (gener->GetCollisionGeometry(ns))->ImpactParameter();
	124
	125	if(fB >= 2.*fR) continue; // no collision
	126	}
	127
be66cf70	128	fCurStack = gener->GetStack(ns);
	129	if(!fCurStack)
	130	{
	131	AliWarning("no event stack");
	132	return;
	133	}
	134
	135	TList* protons = new TList();
	136	protons->SetOwner(kFALSE);
	137	TList* neutrons = new TList();
	138	neutrons->SetOwner(kFALSE);
	139
95145d47	140	// particles produced by the generator
95145d47	141	for (Int_t i=0; i < fCurStack->GetNprimary(); ++i)
be66cf70	142	{
	143	TParticle* iParticle = fCurStack->Particle(i);
	144
95145d47	145	if(iParticle->GetStatusCode() != 1) continue;
be66cf70	146
	147	Int_t pdgCode = iParticle->GetPdgCode();
	148	if(pdgCode == fSign*2212)// (anti)proton
	149	{
	150	FixProductionVertex(iParticle);
	151	protons->Add(iParticle);
	152	}
	153	else if(pdgCode == fSign*2112) // (anti)neutron
	154	{
	155	FixProductionVertex(iParticle);
	156	neutrons->Add(iParticle);
	157	}
	158	}
	159
	160	// look for clusters
	161	if(fClusterType==kFirstPartner)
	162	{
	163	FirstPartner(protons, neutrons);
	164	}
	165	else
	166	{
	167	WeightMatrix(protons, neutrons);
	168	}
	169
	170	protons->Clear("nodelete");
	171	neutrons->Clear("nodelete");
	172
	173	delete protons;
	174	delete neutrons;
	175	}
	176
	177	AliInfo("DONE");
	178	}
	179
	180	Double_t AliGenDeuteron::GetCoalescenceProbability(const TParticle* nucleon1, const TParticle* nucleon2) const
	181	{
	182	//
	183	// Coalescence conditions as in
	184	// A. J. Baltz et al., Phys. lett B 325(1994)7
	185	//
	186	// returns < 0 if coalescence is not possible
	187	// otherwise returns a coalescence probability
	188	//
95145d47	189	const Double_t kProtonMass = 0.938272013;
	190	const Double_t kNeutronMass = 0.939565378;
	191
be66cf70	192	TVector3 v1(nucleon1->Vx(), nucleon1->Vy(), nucleon1->Vz());
	193	TVector3 p1(nucleon1->Px(), nucleon1->Py(), nucleon1->Pz());
	194
	195	TVector3 v2(nucleon2->Vx(), nucleon2->Vy(), nucleon2->Vz());
	196	TVector3 p2(nucleon2->Px(), nucleon2->Py(), nucleon2->Pz());
	197
95145d47	198	Double_t deltaP = 2.*this->GetPcm(p1, kProtonMass, p2, kNeutronMass); // relative momentum in CM frame
95145d47	199	if( deltaP >= fPmax) return -1.;
be66cf70	200
	201	Double_t deltaR = (v2-v1).Mag(); // relative distance (cm)
	202	if(deltaR >= fRmax*1.e-13) return -1.;
	203
	204	if(Rndm() > fSpinProb) return -1.; // spin
	205
	206	if(fClusterType == kLowestMomentum) return 1. - deltaP/fPmax;
	207	if(fClusterType == kLowestDistance) return 1. - 1.e+13*deltaR/fRmax;
	208
	209	return 1. - 1.e+13(deltaPdeltaR)/(fRmax*fPmax);
	210	}
	211
	212	void AliGenDeuteron::FirstPartner(const TList* protons, TList* neutrons)
	213	{
	214	//
	215	// Clusters are made with the first nucleon pair that fulfill
	216	// the coalescence conditions, starting with the protons
	217	//
	218	TIter p_next(protons);
	219	while(TParticle* n0 = (TParticle*) p_next())
	220	{
	221	TParticle* partner = 0;
	222	TIter n_next(neutrons);
	223	while(TParticle* n1 = (TParticle*) n_next() )
	224	{
	225	if(GetCoalescenceProbability(n0, n1) < 0 ) continue; // with next neutron
	226	partner = n1;
	227	break;
	228	}
	229
	230	if(partner == 0) continue; // with next proton
	231
	232	PushDeuteron(n0, partner);
	233
	234	// Remove from the list for the next iteration
	235	neutrons->Remove(partner);
	236	}
	237	}
	238
	239	void AliGenDeuteron::WeightMatrix(const TList* protons, const TList* neutrons)
	240	{
	241	//
	242	// Build all possible nucleon pairs with their own probability
	243	// and select only those with the highest coalescence probability
	244	//
	245	Int_t nMaxPairs = protons->GetSize()*neutrons->GetSize();
	246
	247	TParticle** cProton = new TParticle*[nMaxPairs];
	248	TParticle** cNeutron = new TParticle*[nMaxPairs];
	249	Double_t* cWeight = new Double_t[nMaxPairs];
	250
	251	// build all pairs with probability > 0
	252	Int_t cIdx = -1;
	253	TIter p_next(protons);
	254	while(TParticle* n1 = (TParticle*) p_next())
	255	{
	256	TIter n_next(neutrons);
	257	while(TParticle* n2 = (TParticle*) n_next() )
	258	{
	259	Double_t weight = this->GetCoalescenceProbability(n1,n2);
	260	if(weight < 0) continue;
	261	++cIdx;
	262	cProton[cIdx] = n1;
	263	cNeutron[cIdx] = n2;
264	cWeight[cIdx] = weight;
265	}
266	n_next.Reset();
267	}
268	p_next.Reset();
269
270	Int_t nPairs = cIdx + 1;
271
272	// find the interacting pairs:
273	// remove repeated pairs and select only
274	// the pair with the highest coalescence probability
275
276	Int_t nMaxIntPair = TMath::Min(protons->GetSize(), neutrons->GetSize());
277
278	TParticle** iProton = new TParticle*[nMaxIntPair];
279	TParticle** iNeutron = new TParticle*[nMaxIntPair];
280	Double_t* iWeight = new Double_t[nMaxIntPair];
281
282	Int_t iIdx = -1;
95145d47	283	while(kTRUE)
be66cf70	284	{
	285	Int_t j = -1;
	286	Double_t wMax = 0;
	287	for(Int_t i=0; i < nPairs; ++i)
	288	{
	289	if(cWeight[i] > wMax)
	290	{
	291	wMax=cWeight[i];
	292	j = i;
	293	}
	294	}
	295
	296	if(j == -1 ) break; // end
	297
	298	// Save the interacting pair
	299	++iIdx;
	300	iProton[iIdx] = cProton[j];
	301	iNeutron[iIdx] = cNeutron[j];
	302	iWeight[iIdx] = cWeight[j];
	303
	304	// invalidate all combinations with these pairs for the next iteration
	305	for(Int_t i=0; i < nPairs; ++i)
	306	{
	307	if(cProton[i] == iProton[iIdx]) cWeight[i] = -1.;
	308	if(cNeutron[i] == iNeutron[iIdx]) cWeight[i] = -1.;
	309	}
	310	}
	311
	312	Int_t nIntPairs = iIdx + 1;
	313
	314	delete[] cProton;
	315	delete[] cNeutron;
	316	delete[] cWeight;
	317
	318	// Add the (anti)deuterons to the current event stack
	319	for(Int_t i=0; i<nIntPairs; ++i)
	320	{
	321	TParticle* n1 = iProton[i];
	322	TParticle* n2 = iNeutron[i];
	323	PushDeuteron(n1,n2);
	324	}
	325
	326	delete[] iProton;
	327	delete[] iNeutron;
	328	delete[] iWeight;
	329	}
	330
	331	void AliGenDeuteron::PushDeuteron(TParticle* parent1, TParticle* parent2)
	332	{
	333	//
	334	// Create an (anti)deuteron from parent1 and parent2,
	335	// add to the current stack and set kDoneBit for the parents
	336	//
	337	const Double_t kDeuteronMass = 1.87561282;
	338	const Int_t kDeuteronPdg = 1000010020;
	339
	340	// momentum
	341	TVector3 p1(parent1->Px(), parent1->Py(), parent1->Pz());
	342	TVector3 p2(parent2->Px(), parent2->Py(), parent2->Pz());
	343	TVector3 pN = p1+p2;
	344
	345	// production vertex same as the parent1's
	346	TVector3 vN(parent1->Vx(), parent1->Vy(), parent1->Vz());
	347
348	// E^2 = p^2 + m^2
349	Double_t energy = TMath::Sqrt(pN.Mag2() + kDeuteronMass*kDeuteronMass);
350
3d87465c	351	Int_t ntrk = 0;
	352	Double_t weight = 1;
	353	Int_t is = 1; // final state particle
	354
be66cf70	355	// Add a new (anti)deuteron to current event stack
0f52fba4	356	fCurStack->PushTrack(1, -1, fSign*kDeuteronPdg,
be66cf70	357	pN.X(), pN.Y(), pN.Z(), energy,
be66cf70	358	vN.X(), vN.Y(), vN.Z(), parent1->T(),
3d87465c	359	0., 0., 0., kPNCapture, ntrk, weight, is);
	360
	361	// change the status code of the parents
	362	parent1->SetStatusCode(kCluster);
	363	parent2->SetStatusCode(kCluster);
be66cf70	364
	365	// Set kDoneBit for the parents
	366	parent1->SetBit(kDoneBit);
	367	parent2->SetBit(kDoneBit);
	368	}
	369
	370	void AliGenDeuteron::FixProductionVertex(TParticle* i)
	371	{
	372	//
	373	// Replace current generator nucleon spatial distribution
	374	// with a custom distribution according to the selected model
	375	//
	376	if(fModel == kNone \|\| fModel > kExpansion) return;
	377
	378	// semi-axis from collision geometry (fm)
	379	Double_t a = fTimeLength + TMath::Sqrt(fRfR - fBfB/4.);
	380	Double_t b = fTimeLength + fR - fB/2.;
	381	Double_t c = fTimeLength;
	382
	383	Double_t xx = 0;
	384	Double_t yy = 0;
	385	Double_t zz = 0;
	386
	387	if(fModel == kThermal)
	388	{
	389	// uniformly ditributed in the volume on an ellipsoid
	390	// random (r,theta,phi) unit sphere
	391	Double_t r = TMath::Power(Rndm(),1./3.);
	392	Double_t theta = TMath::ACos(2.*Rndm()-1.);
	393	Double_t phi = 2.TMath::Pi()Rndm();
	394
	395	// transform coordenates
	396	xx = arTMath::Sin(theta)*TMath::Cos(phi);
	397	yy = brTMath::Sin(theta)*TMath::Sin(phi);
	398	zz = crTMath::Cos(theta);
	399	}
	400	else if(fModel == kExpansion)
	401	{
	402	// project into the surface of an ellipsoid
	403	xx = aTMath::Sin(i->Theta())TMath::Cos(i->Phi());
	404	yy = bTMath::Sin(i->Theta())TMath::Sin(i->Phi());
	405	zz = c*TMath::Cos(i->Theta());
	406	}
	407
	408	// rotate by the reaction plane angle
	409	Double_t x = xxTMath::Cos(fPsiR)+yyTMath::Sin(fPsiR);
	410	Double_t y = -xxTMath::Sin(fPsiR)+yyTMath::Cos(fPsiR);
	411	Double_t z = zz;
	412
	413	// translate by the production vertex (cm)
	414	i->SetProductionVertex(i->Vx() + 1.e-13x, i->Vy() + 1.e-13y, i->Vz() + 1.e-13*z, i->T());
	415	}
	416
95145d47	417	Double_t AliGenDeuteron::GetS(Double_t p1x, Double_t p1y, Double_t p1z, Double_t m1, Double_t p2x, Double_t p2y, Double_t p2z, Double_t m2) const
	418	{
	419	//
	420	// square of the center of mass energy
	421	//
	422	Double_t E1 = TMath::Sqrt( p1xp1x + p1yp1y + p1zp1z + m1m1);
	423	Double_t E2 = TMath::Sqrt( p2xp2x + p2yp2y + p2zp2z + m2m2);
	424
	425	return (E1+E2)(E1+E2) - ((p1x+p2x)(p1x+p2x) + (p1y+p2y)(p1y+p2y) + (p1z+p2z)(p1z+p2z));
	426	}
	427
	428	Double_t AliGenDeuteron::GetPcm(Double_t p1x, Double_t p1y, Double_t p1z, Double_t m1, Double_t p2x, Double_t p2y, Double_t p2z, Double_t m2) const
	429	{
	430	//
	431	// momentum in the CM frame
	432	//
	433	Double_t s = this->GetS(p1x, p1y, p1z, m1, p2x, p2y, p2z, m2);
	434
	435	return TMath::Sqrt( (s-(m1-m2)(m1-m2))(s-(m1+m2)(m1+m2)) )/(2.TMath::Sqrt(s));
	436	}
	437
	438	Double_t AliGenDeuteron::GetPcm(const TVector3& p1, Double_t m1, const TVector3& p2, Double_t m2) const
	439	{
	440	//
	441	// momentum in the CM frame
	442	//
	443	return this->GetPcm(p1.Px(),p1.Py(),p1.Pz(),m1,p2.Px(),p2.Py(),p2.Pz(),m2);
	444	}
	445