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

/**************************************************************************
 * Copyright(c) 2009-2014, 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 light nuclei for event generators
// such as PYTHIA and PHOJET
//
// Light nuclei are generated whenever a cluster of nucleons is found
// within a sphere of radius p0 (coalescence momentum), i.e. have the
// same momentum.
//
// By default it starts with He4 nuclei which are the most stable,
// then He3 nuclei, tritons and finally deuterons. It can also generate
// a single nucleus species by disabling the others.
//
// Sample code for PYTHIA:
//
//    AliGenLightNuclei* gener = new AliGenLightNuclei();
//
//    AliGenPythia* pythia = new AliGenPythia(-1);
//    pythia->SetCollisionSystem("p+", "p+");
//    pythia->SetEnergyCMS(7000);
//
//    gener->UsePerEventRates();
//    gener->AddGenerator(pythia, "PYTHIA", 1);
//    gener->SetCoalescenceMomentum(0.100); // default (GeV/c)
//
//////////////////////////////////////////////////////////////////////

//
// Author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
//

#include "TMath.h"
#include "TPDGCode.h"
#include "TMCProcess.h"
#include "TList.h"
#include "TVector3.h"
#include "TParticle.h"
#include "AliStack.h"
#include "AliRun.h"
#include "AliLog.h"
#include "AliMC.h"
#include "AliGenLightNuclei.h"

ClassImp(AliGenLightNuclei)

AliGenLightNuclei::AliGenLightNuclei()
:AliGenCocktail()
,fP0(0.100)
,fGenDeuterons(kTRUE)
,fGenTritons(kTRUE)
,fGenHe3Nuclei(kTRUE)
,fGenHe4Nuclei(kTRUE)
{
//
// default constructor
//
}

AliGenLightNuclei::~AliGenLightNuclei()
{
//
// default destructor
//
}

void AliGenLightNuclei::Generate()
{
//
// delegate the particle generation to the cocktail
// and modify the stack adding the light nuclei
//
	AliGenCocktail::Generate();
	
	// find the nucleons and anti-nucleons
	
	TList* protons      = new TList();
	TList* neutrons     = new TList();
	TList* antiprotons  = new TList();
	TList* antineutrons = new TList();
	
	for (Int_t i=0; i < fStack->GetNprimary(); ++i)
	{
		TParticle* iParticle = fStack->Particle(i);
		
		if(iParticle->GetStatusCode() != 1) continue;
		
		switch(iParticle->GetPdgCode())
		{
			case kProton:
				protons->Add(iParticle);
				break;
			case kNeutron:
				neutrons->Add(iParticle);
				break;
			case kProtonBar:
				antiprotons->Add(iParticle);
				break;
			case kNeutronBar:
				antineutrons->Add(iParticle);
				break;
			default:
				break;
		}
	}
	
	// do not delete content
	protons->SetOwner(kFALSE);
	neutrons->SetOwner(kFALSE);
	antiprotons->SetOwner(kFALSE);
	antineutrons->SetOwner(kFALSE);
	
	// first try with He4 nuclei which are the most stable
	
	if(fGenHe4Nuclei)
	{
		this->GenerateHe4Nuclei(protons, neutrons);
		this->GenerateHe4Nuclei(antiprotons, antineutrons);
	}
	
	// then He3 nuclei
	
	if(fGenHe3Nuclei)
	{
		this->GenerateHe3Nuclei(protons, neutrons);
		this->GenerateHe3Nuclei(antiprotons, antineutrons);
	}
	
	// then tritons
	
	if(fGenTritons)
	{
		this->GenerateTritons(protons, neutrons);
		this->GenerateTritons(antiprotons, antineutrons);
	}
	
	// and finally deuterons
	
	if(fGenDeuterons)
	{
		this->GenerateDeuterons(protons, neutrons);
		this->GenerateDeuterons(antiprotons, antineutrons);
	}
	
	protons->Clear("nodelete");
	neutrons->Clear("nodelete");
	antiprotons->Clear("nodelete");
	antineutrons->Clear("nodelete");
	
	delete protons;
	delete neutrons;
	delete antiprotons;
	delete antineutrons;
}

Bool_t AliGenLightNuclei::Coalescence(const TParticle* n1, const TParticle* n2) const
{
//
// returns true if the nucleons are inside of an sphere of radius p0
// (assume the nucleons are in the same place e.g. PYTHIA, PHOJET,...)
//
	Double_t deltaP = this->GetPcm(n1->Px(), n1->Py(), n1->Pz(), n1->GetMass(), 
	                               n2->Px(), n2->Py(), n2->Pz(), n2->GetMass());
	
	return ( deltaP < fP0);
}

TParticle* AliGenLightNuclei::FindPartner(const TParticle* n0, const TList* nucleons, const TParticle* nx) const
{
//
// find the first nucleon partner within a sphere of radius p0
// centered at n0 and exclude nucleon nx
//
	TIter n_next(nucleons);
	while(TParticle* n1 = dynamic_cast<TParticle*>( n_next()) )
	{
		if(n1 == 0) continue;
		if(n1 == nx) continue;
		if(n1->GetStatusCode() == kCluster) continue;
		if( !this->Coalescence(n0, n1) ) continue;
		return n1;
	}
	
	return 0;
}

Int_t AliGenLightNuclei::GenerateDeuterons(const TList* protons, const TList* neutrons)
{
//
// a deuteron is generated from a pair of p-n nucleons
// (the center of the sphere is one of the nucleons)
//
	Int_t npart = 0;
	
	TIter p_next(protons);
	while(TParticle* n0 = dynamic_cast<TParticle*>(p_next()) )
	{
		if(n0 == 0) continue;
		if(n0->GetStatusCode() == kCluster) continue;
		
		TParticle* partner = this->FindPartner(n0, neutrons, 0);
		
		if(partner == 0) continue;
		
		this->PushDeuteron(n0, partner);
		
		++npart;
	}
	
	return npart;
}

Int_t AliGenLightNuclei::GenerateTritons(const TList* protons, const TList* neutrons)
{
//
// a triton is generated from a cluster of p-n-n nucleons with same momentum
// (triangular configuration)
//
	Int_t npart = 0;
	
	TIter p_next(protons);
	while(TParticle* n0 = dynamic_cast<TParticle*>(p_next()) )
	{
		if(n0 == 0) continue;
		if(n0->GetStatusCode() == kCluster) continue;
		
		TParticle* partner1 = this->FindPartner(n0, neutrons, 0);
		
		if(partner1 == 0) continue;
		
		TParticle* partner2 = this->FindPartner(n0, neutrons, partner1);
		
		if(partner2 == 0) continue;
		
		// check that the partners coalesce between themselves
		
		if(!this->Coalescence(partner1, partner2)) continue;
		
		this->PushTriton(n0, partner1, partner2);
		
		++npart;
	}
	
	return npart;
}

Int_t AliGenLightNuclei::GenerateHe3Nuclei(const TList* protons, const TList* neutrons)
{
//
// a He3 nucleus is generated from a cluster of p-n-p nucleons with same momentum
// (triangular configuration)
//
	Int_t npart = 0;
	
	TIter p_next(protons); // center of the sphere
	while(TParticle* n0 = dynamic_cast<TParticle*>(p_next()) )
	{
		if(n0 == 0) continue;
		if(n0->GetStatusCode() == kCluster) continue;
		
		TParticle* partner1 = this->FindPartner(n0, neutrons, 0);
		
		if(partner1 == 0) continue;
		
		TParticle* partner2 = this->FindPartner(n0, protons, n0);
		
		if(partner2 == 0) continue;
		
		// check that the partners coalesce between themselves
		
		if(!this->Coalescence(partner1, partner2)) continue;
		
		this->PushHe3Nucleus(n0, partner1, partner2);
		
		++npart;
	}
	
	return npart;
}

Int_t AliGenLightNuclei::GenerateHe4Nuclei(const TList* protons, const TList* neutrons)
{
//
// a He4 nucleus is generated from a cluster of p-n-p-n nucleons with same momentum
// (tetrahedron configuration)
//
	Int_t npart = 0;
	
	TIter p_next(protons); // center of the sphere
	while(TParticle* n0 = dynamic_cast<TParticle*>(p_next()) )
	{
		if(n0 == 0) continue;
		if(n0->GetStatusCode() == kCluster) continue;
		
		TParticle* partner1 = this->FindPartner(n0, neutrons, 0);
		
		if(partner1 == 0) continue;
		
		TParticle* partner2 = this->FindPartner(n0, protons, n0);
		
		if(partner2 == 0) continue;
		
		TParticle* partner3 = this->FindPartner(n0, neutrons, partner1);
		
		if(partner3 == 0) continue;
		
		// check that the partners coalesce between themselves
		
		if(!this->Coalescence(partner1, partner2)) continue;
		if(!this->Coalescence(partner1, partner3)) continue;
		if(!this->Coalescence(partner2, partner3)) continue;
		
		this->PushHe4Nucleus(n0, partner1, partner2, partner3 );
		
		++npart;
	}
	
	return npart;
}

void AliGenLightNuclei::PushDeuteron(TParticle* parent1, TParticle* parent2)
{
//
// push a deuteron to the particle stack
//
	Int_t pdg = ( parent1->GetPdgCode() > 0 ) ? kDeuteron : kAntiDeuteron;
	this->PushNucleus(pdg, 1.87561282, parent1, parent2);
}

void AliGenLightNuclei::PushTriton(TParticle* parent1, TParticle* parent2, TParticle* parent3)
{
//
// push a triton to the particle stack
//
	Int_t pdg = ( parent1->GetPdgCode() > 0 ) ? kTriton : kAntiTriton;
	this->PushNucleus(pdg, 2.80925, parent1, parent2, parent3);
}

void AliGenLightNuclei::PushHe3Nucleus(TParticle* parent1, TParticle* parent2, TParticle* parent3)
{
//
// push a He3 nucleus to the particle stack
//
	Int_t pdg = ( parent1->GetPdgCode() > 0 ) ? kHe3Nucleus : kAntiHe3Nucleus;
	this->PushNucleus(pdg, 2.80923, parent1, parent2, parent3);
}

void AliGenLightNuclei::PushHe4Nucleus(TParticle* parent1, TParticle* parent2, TParticle* parent3, TParticle* parent4)
{
//
// push a He4 nucleus to the particle stack
//
	Int_t pdg = ( parent1->GetPdgCode() > 0 ) ? kAlpha : kAntiAlpha;
	this->PushNucleus(pdg, 3.727417, parent1, parent2, parent3, parent4);
}

void AliGenLightNuclei::PushNucleus(Int_t pdg, Double_t mass, TParticle* parent1, TParticle* parent2, TParticle* parent3, TParticle* parent4)
{
//
// push a nucleus to the stack and tag the parents with the kCluster status code
//
	Int_t ntrk;
	
	// momentum
	TVector3 p1(parent1->Px(), parent1->Py(), parent1->Pz());
	TVector3 p2(parent2->Px(), parent2->Py(), parent2->Pz());
	TVector3 p3(0, 0, 0);
	TVector3 p4(0, 0, 0);
	if(parent3 != 0) p3.SetXYZ(parent3->Px(), parent3->Py(), parent3->Pz());
	if(parent4 != 0) p4.SetXYZ(parent4->Px(), parent4->Py(), parent4->Pz());
	
	// momentum
	TVector3 pN = p1+p2+p3+p4;
	
	// E^2 = p^2 + m^2
	Double_t energy = TMath::Sqrt(pN.Mag2() + mass*mass);
	
	// production vertex same as the parent1'
	TVector3 vN(parent1->Vx(), parent1->Vy(), parent1->Vz());
	
	Double_t weight = 1;
	Int_t is = 1; // final state particle
	
	// add a new nucleus to current event stack
	fStack->PushTrack(1, -1, pdg,
	                 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);
	if(parent3 != 0) parent3->SetStatusCode(kCluster);
	if(parent4 != 0) parent4->SetStatusCode(kCluster);
	
	// set no transport for the parents
	parent1->SetBit(kDoneBit);
	parent2->SetBit(kDoneBit);
	if(parent3 != 0) parent3->SetBit(kDoneBit);
	if(parent4 != 0) parent4->SetBit(kDoneBit);
}

Double_t AliGenLightNuclei::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 AliGenLightNuclei::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 for 2 particles
//
	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));
}
Commit	Line	Data
a36106b2	1	/**************************************************************************
3fdac9ed	2	* Copyright(c) 2009-2014, ALICE Experiment at CERN, All rights reserved. *
a36106b2	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	//
	18	// Afterburner to generate light nuclei for event generators
	19	// such as PYTHIA and PHOJET
	20	//
	21	// Light nuclei are generated whenever a cluster of nucleons is found
	22	// within a sphere of radius p0 (coalescence momentum), i.e. have the
	23	// same momentum.
	24	//
	25	// By default it starts with He4 nuclei which are the most stable,
3fdac9ed	26	// then He3 nuclei, tritons and finally deuterons. It can also generate
a36106b2	27	// a single nucleus species by disabling the others.
	28	//
	29	// Sample code for PYTHIA:
	30	//
	31	// AliGenLightNuclei* gener = new AliGenLightNuclei();
	32	//
	33	// AliGenPythia* pythia = new AliGenPythia(-1);
	34	// pythia->SetCollisionSystem("p+", "p+");
	35	// pythia->SetEnergyCMS(7000);
	36	//
	37	// gener->UsePerEventRates();
	38	// gener->AddGenerator(pythia, "PYTHIA", 1);
3fdac9ed	39	// gener->SetCoalescenceMomentum(0.100); // default (GeV/c)
a36106b2	40	//
	41	//////////////////////////////////////////////////////////////////////
	42
	43	//
	44	// Author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
	45	//
	46
	47	#include "TMath.h"
	48	#include "TPDGCode.h"
	49	#include "TMCProcess.h"
	50	#include "TList.h"
	51	#include "TVector3.h"
	52	#include "TParticle.h"
	53	#include "AliStack.h"
	54	#include "AliRun.h"
	55	#include "AliLog.h"
	56	#include "AliMC.h"
	57	#include "AliGenLightNuclei.h"
	58
	59	ClassImp(AliGenLightNuclei)
	60
	61	AliGenLightNuclei::AliGenLightNuclei()
	62	:AliGenCocktail()
3fdac9ed	63	,fP0(0.100)
a36106b2	64	,fGenDeuterons(kTRUE)
	65	,fGenTritons(kTRUE)
	66	,fGenHe3Nuclei(kTRUE)
	67	,fGenHe4Nuclei(kTRUE)
	68	{
	69	//
	70	// default constructor
	71	//
	72	}
	73
	74	AliGenLightNuclei::~AliGenLightNuclei()
	75	{
	76	//
	77	// default destructor
	78	//
	79	}
	80
	81	void AliGenLightNuclei::Generate()
	82	{
	83	//
	84	// delegate the particle generation to the cocktail
	85	// and modify the stack adding the light nuclei
	86	//
	87	AliGenCocktail::Generate();
	88
	89	// find the nucleons and anti-nucleons
	90
	91	TList* protons = new TList();
	92	TList* neutrons = new TList();
	93	TList* antiprotons = new TList();
	94	TList* antineutrons = new TList();
	95
	96	for (Int_t i=0; i < fStack->GetNprimary(); ++i)
	97	{
	98	TParticle* iParticle = fStack->Particle(i);
	99
	100	if(iParticle->GetStatusCode() != 1) continue;
	101
	102	switch(iParticle->GetPdgCode())
	103	{
	104	case kProton:
	105	protons->Add(iParticle);
	106	break;
	107	case kNeutron:
	108	neutrons->Add(iParticle);
	109	break;
	110	case kProtonBar:
	111	antiprotons->Add(iParticle);
	112	break;
	113	case kNeutronBar:
	114	antineutrons->Add(iParticle);
	115	break;
	116	default:
	117	break;
	118	}
	119	}
	120
	121	// do not delete content
	122	protons->SetOwner(kFALSE);
	123	neutrons->SetOwner(kFALSE);
	124	antiprotons->SetOwner(kFALSE);
	125	antineutrons->SetOwner(kFALSE);
	126
	127	// first try with He4 nuclei which are the most stable
128
129	if(fGenHe4Nuclei)
130	{
131	this->GenerateHe4Nuclei(protons, neutrons);
132	this->GenerateHe4Nuclei(antiprotons, antineutrons);
133	}
134
135	// then He3 nuclei
136
137	if(fGenHe3Nuclei)
138	{
139	this->GenerateHe3Nuclei(protons, neutrons);
140	this->GenerateHe3Nuclei(antiprotons, antineutrons);
141	}
142
143	// then tritons
144
145	if(fGenTritons)
146	{
147	this->GenerateTritons(protons, neutrons);
148	this->GenerateTritons(antiprotons, antineutrons);
149	}
150
151	// and finally deuterons
152
153	if(fGenDeuterons)
154	{
155	this->GenerateDeuterons(protons, neutrons);
156	this->GenerateDeuterons(antiprotons, antineutrons);
157	}
158
159	protons->Clear("nodelete");
160	neutrons->Clear("nodelete");
161	antiprotons->Clear("nodelete");
162	antineutrons->Clear("nodelete");
163
164	delete protons;
165	delete neutrons;
166	delete antiprotons;
167	delete antineutrons;
168	}
169
170	Bool_t AliGenLightNuclei::Coalescence(const TParticle* n1, const TParticle* n2) const
171	{
172	//
173	// returns true if the nucleons are inside of an sphere of radius p0
174	// (assume the nucleons are in the same place e.g. PYTHIA, PHOJET,...)
175	//
3fdac9ed	176	Double_t deltaP = this->GetPcm(n1->Px(), n1->Py(), n1->Pz(), n1->GetMass(),
3fdac9ed	177	n2->Px(), n2->Py(), n2->Pz(), n2->GetMass());
a36106b2	178
	179	return ( deltaP < fP0);
	180	}
	181
	182	TParticle* AliGenLightNuclei::FindPartner(const TParticle* n0, const TList* nucleons, const TParticle* nx) const
	183	{
	184	//
	185	// find the first nucleon partner within a sphere of radius p0
	186	// centered at n0 and exclude nucleon nx
	187	//
	188	TIter n_next(nucleons);
	189	while(TParticle* n1 = dynamic_cast<TParticle*>( n_next()) )
	190	{
	191	if(n1 == 0) continue;
	192	if(n1 == nx) continue;
	193	if(n1->GetStatusCode() == kCluster) continue;
	194	if( !this->Coalescence(n0, n1) ) continue;
	195	return n1;
	196	}
	197
	198	return 0;
	199	}
	200
	201	Int_t AliGenLightNuclei::GenerateDeuterons(const TList* protons, const TList* neutrons)
	202	{
	203	//
	204	// a deuteron is generated from a pair of p-n nucleons
	205	// (the center of the sphere is one of the nucleons)
	206	//
	207	Int_t npart = 0;
	208
	209	TIter p_next(protons);
	210	while(TParticle* n0 = dynamic_cast<TParticle*>(p_next()) )
	211	{
	212	if(n0 == 0) continue;
	213	if(n0->GetStatusCode() == kCluster) continue;
	214
	215	TParticle* partner = this->FindPartner(n0, neutrons, 0);
	216
	217	if(partner == 0) continue;
	218
	219	this->PushDeuteron(n0, partner);
	220
	221	++npart;
	222	}
	223
	224	return npart;
	225	}
	226
	227	Int_t AliGenLightNuclei::GenerateTritons(const TList* protons, const TList* neutrons)
	228	{
	229	//
	230	// a triton is generated from a cluster of p-n-n nucleons with same momentum
	231	// (triangular configuration)
	232	//
	233	Int_t npart = 0;
	234
	235	TIter p_next(protons);
	236	while(TParticle* n0 = dynamic_cast<TParticle*>(p_next()) )
	237	{
	238	if(n0 == 0) continue;
	239	if(n0->GetStatusCode() == kCluster) continue;
	240
	241	TParticle* partner1 = this->FindPartner(n0, neutrons, 0);
242
243	if(partner1 == 0) continue;
244
245	TParticle* partner2 = this->FindPartner(n0, neutrons, partner1);
246
247	if(partner2 == 0) continue;
248
249	// check that the partners coalesce between themselves
250
251	if(!this->Coalescence(partner1, partner2)) continue;
252
253	this->PushTriton(n0, partner1, partner2);
254
255	++npart;
256	}
257
258	return npart;
259	}
260
261	Int_t AliGenLightNuclei::GenerateHe3Nuclei(const TList* protons, const TList* neutrons)
262	{
263	//
264	// a He3 nucleus is generated from a cluster of p-n-p nucleons with same momentum
265	// (triangular configuration)
266	//
267	Int_t npart = 0;
268
269	TIter p_next(protons); // center of the sphere
270	while(TParticle* n0 = dynamic_cast<TParticle*>(p_next()) )
271	{
272	if(n0 == 0) continue;
273	if(n0->GetStatusCode() == kCluster) continue;
274
275	TParticle* partner1 = this->FindPartner(n0, neutrons, 0);
276
277	if(partner1 == 0) continue;
278
279	TParticle* partner2 = this->FindPartner(n0, protons, n0);
280
281	if(partner2 == 0) continue;
282
283	// check that the partners coalesce between themselves
284
285	if(!this->Coalescence(partner1, partner2)) continue;
286
287	this->PushHe3Nucleus(n0, partner1, partner2);
288
289	++npart;
290	}
291
292	return npart;
293	}
294
295	Int_t AliGenLightNuclei::GenerateHe4Nuclei(const TList* protons, const TList* neutrons)
296	{
297	//
298	// a He4 nucleus is generated from a cluster of p-n-p-n nucleons with same momentum
299	// (tetrahedron configuration)
300	//
301	Int_t npart = 0;
302
303	TIter p_next(protons); // center of the sphere
304	while(TParticle* n0 = dynamic_cast<TParticle*>(p_next()) )
305	{
306	if(n0 == 0) continue;
307	if(n0->GetStatusCode() == kCluster) continue;
308
309	TParticle* partner1 = this->FindPartner(n0, neutrons, 0);
310
311	if(partner1 == 0) continue;
312
313	TParticle* partner2 = this->FindPartner(n0, protons, n0);
314
315	if(partner2 == 0) continue;
316
317	TParticle* partner3 = this->FindPartner(n0, neutrons, partner1);
318
319	if(partner3 == 0) continue;
320
321	// check that the partners coalesce between themselves
322
323	if(!this->Coalescence(partner1, partner2)) continue;
324	if(!this->Coalescence(partner1, partner3)) continue;
325	if(!this->Coalescence(partner2, partner3)) continue;
326
327	this->PushHe4Nucleus(n0, partner1, partner2, partner3 );
328
329	++npart;
330	}
331
332	return npart;
333	}
334
335	void AliGenLightNuclei::PushDeuteron(TParticle* parent1, TParticle* parent2)
336	{
337	//
338	// push a deuteron to the particle stack
339	//
340	Int_t pdg = ( parent1->GetPdgCode() > 0 ) ? kDeuteron : kAntiDeuteron;
341	this->PushNucleus(pdg, 1.87561282, parent1, parent2);
342	}
343
344	void AliGenLightNuclei::PushTriton(TParticle* parent1, TParticle* parent2, TParticle* parent3)
345	{
346	//
347	// push a triton to the particle stack
348	//
349	Int_t pdg = ( parent1->GetPdgCode() > 0 ) ? kTriton : kAntiTriton;
350	this->PushNucleus(pdg, 2.80925, parent1, parent2, parent3);
351	}
352
353	void AliGenLightNuclei::PushHe3Nucleus(TParticle* parent1, TParticle* parent2, TParticle* parent3)
354	{
355	//
356	// push a He3 nucleus to the particle stack
357	//
358	Int_t pdg = ( parent1->GetPdgCode() > 0 ) ? kHe3Nucleus : kAntiHe3Nucleus;
359	this->PushNucleus(pdg, 2.80923, parent1, parent2, parent3);
360	}
361
362	void AliGenLightNuclei::PushHe4Nucleus(TParticle* parent1, TParticle* parent2, TParticle* parent3, TParticle* parent4)
363	{
364	//
365	// push a He4 nucleus to the particle stack
366	//
367	Int_t pdg = ( parent1->GetPdgCode() > 0 ) ? kAlpha : kAntiAlpha;
368	this->PushNucleus(pdg, 3.727417, parent1, parent2, parent3, parent4);
369	}
370
371	void AliGenLightNuclei::PushNucleus(Int_t pdg, Double_t mass, TParticle* parent1, TParticle* parent2, TParticle* parent3, TParticle* parent4)
372	{
373	//
374	// push a nucleus to the stack and tag the parents with the kCluster status code
375	//
376	Int_t ntrk;
377
378	// momentum
379	TVector3 p1(parent1->Px(), parent1->Py(), parent1->Pz());
380	TVector3 p2(parent2->Px(), parent2->Py(), parent2->Pz());
381	TVector3 p3(0, 0, 0);
382	TVector3 p4(0, 0, 0);
383	if(parent3 != 0) p3.SetXYZ(parent3->Px(), parent3->Py(), parent3->Pz());
384	if(parent4 != 0) p4.SetXYZ(parent4->Px(), parent4->Py(), parent4->Pz());
385
386	// momentum
387	TVector3 pN = p1+p2+p3+p4;
388
389	// E^2 = p^2 + m^2
390	Double_t energy = TMath::Sqrt(pN.Mag2() + mass*mass);
391
392	// production vertex same as the parent1'
393	TVector3 vN(parent1->Vx(), parent1->Vy(), parent1->Vz());
394
395	Double_t weight = 1;
396	Int_t is = 1; // final state particle
397
398	// add a new nucleus to current event stack
399	fStack->PushTrack(1, -1, pdg,
400	pN.X(), pN.Y(), pN.Z(), energy,
401	vN.X(), vN.Y(), vN.Z(), parent1->T(),
402	0., 0., 0., kPNCapture, ntrk, weight, is);
403
404	// change the status code of the parents
405	parent1->SetStatusCode(kCluster);
406	parent2->SetStatusCode(kCluster);
407	if(parent3 != 0) parent3->SetStatusCode(kCluster);
408	if(parent4 != 0) parent4->SetStatusCode(kCluster);
409
410	// set no transport for the parents
411	parent1->SetBit(kDoneBit);
412	parent2->SetBit(kDoneBit);
413	if(parent3 != 0) parent3->SetBit(kDoneBit);
414	if(parent4 != 0) parent4->SetBit(kDoneBit);
415	}
416
417	Double_t AliGenLightNuclei::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 AliGenLightNuclei::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 for 2 particles
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	}