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

/**************************************************************************\r
 * Copyright(c) 2009-2010, ALICE Experiment at CERN, All rights reserved. *\r
 *                                                                        *\r
 * Author: The ALICE Off-line Project.                                    *\r
 * Contributors are mentioned in the code where appropriate.              *\r
 *                                                                        *\r
 * Permission to use, copy, modify and distribute this software and its   *\r
 * documentation strictly for non-commercial purposes is hereby granted   *\r
 * without fee, provided that the above copyright notice appears in all   *\r
 * copies and that both the copyright notice and this permission notice   *\r
 * appear in the supporting documentation. The authors make no claims     *\r
 * about the suitability of this software for any purpose. It is          *\r
 * provided "as is" without express or implied warranty.                  *\r
 **************************************************************************/\r
\r
//\r
// A very simple model based on the article "Physics Letters B325 (1994) 7-12".\r
// The only addition is to model the freeze-out at which the\r
// light nuclei can form for the energies of the LHC, since PYTHIA does not give\r
// any detailed spatial description of the generated particles at the level of a\r
// few fermis.\r
//\r
// There are two options to place the nucleons: a thermal picture where all\r
// nucleons are placed randomly and homogeneously in an spherical volume of a\r
// few fermis, and expansion one where they are projected on its surface.\r
//\r
// A (anti)deuteron will form if there is a pair of (anti)proton-(anti)neutron\r
// with momentum difference less than ~ 300MeV and relative distance less than a\r
// ~ 2.1fm. Only 3/4 of this clusters are accepted due to the deuteron spin.\r
//\r
// Author: Eulogio Serradilla <eulogio.serradilla@ciemat.es>,\r
//         Arturo Menchaca <menchaca@fisica.unam.mx>\r
//\r
\r
#include "Riostream.h"\r
#include "TParticle.h"\r
#include "AliStack.h"\r
#include "AliGenDeuteron.h"\r
#include "AliGenCocktailAfterBurner.h"\r
#include "TMath.h"\r
#include "TMCProcess.h"\r
#include "TList.h"\r
#include "TVector3.h"\r
#include "AliMC.h"\r
#include "TArrayF.h"\r
#include "AliGenCocktailEventHeader.h"\r
\r
ClassImp(AliGenDeuteron)\r
\r
AliGenDeuteron::AliGenDeuteron(Int_t sign, Double_t pmax, Double_t rmax, Double_t rsrc, Int_t model):\r
fDeuteronMass(1.87561282), //pdg\r
fPmax(pmax), // GeV/c\r
fRmax(rmax*1.e-13), //cm\r
fRsrc(rsrc*1.e-13), // cm\r
fSpinProb(0.75),\r
fModel(model),\r
fSign(1)\r
{\r
	fSign = sign > 0 ? 1:-1;\r
}\r
\r
AliGenDeuteron::~AliGenDeuteron()\r
{\r
}\r
\r
void AliGenDeuteron::Init()\r
{\r
	// Standard AliGenerator Initializer\r
}\r
\r
void AliGenDeuteron::Generate()\r
{\r
// Modify the stack of each event, adding the new particles at the end and\r
// not transporting the parents.\r
\r
	Info("Generate","freeze-out model : %d (0 expansion, 1 thermal)",fModel);\r
	Info("Generate","relative momentum: %g GeV/c",fPmax);\r
	Info("Generate","relative distance: %g cm",fRmax);\r
	Info("Generate","source radius    : %g cm",fRsrc);\r
	Info("Generate","spin probability : %g ",fSpinProb);\r
	Info("Generate","sign             : %d ",fSign);\r
	\r
	Int_t ntr;\r
	\r
	// Get the cocktail generator\r
	AliGenCocktailAfterBurner* gener = (AliGenCocktailAfterBurner*)gAlice->GetMCApp()->Generator();\r
	\r
	// Loop over events\r
	for(Int_t ns = 0; ns < gener->GetNumberOfEvents(); ++ns)\r
	{\r
		gener->SetActiveEventNumber(ns);\r
		\r
		AliStack* stack = gener->GetStack(ns); // Stack of event ns\r
		if(!stack)\r
		{\r
			Info("Generate", "no stack, exiting");\r
			return;\r
		}\r
		\r
		// find emerging nucleons from the collision of current event\r
		\r
		TList* protons = new TList();\r
		protons->SetOwner(kFALSE);\r
		TList* neutrons = new TList();\r
		neutrons->SetOwner(kFALSE);\r
		\r
		// primary vertex of current event\r
		TArrayF primVtx;\r
		gener->GetActiveEventHeader()->PrimaryVertex(primVtx);\r
		TVector3 r0(primVtx[0],primVtx[1],primVtx[2]);\r
		\r
		for (Int_t i=0; i<stack->GetNtrack(); ++i)\r
		{\r
			TParticle* iParticle = stack->Particle(i);\r
			\r
			if(iParticle->TestBit(kDoneBit)) continue;\r
			// select only nucleons within the freeze-out volume\r
			TVector3 r(iParticle->Vx(),iParticle->Vy(),iParticle->Vz());\r
			if((r-r0).Mag()>fRsrc) continue;\r
			\r
			Int_t pdgCode = iParticle->GetPdgCode();\r
			if(pdgCode == fSign*2212)// (anti)proton\r
			{\r
				FixProductionVertex(iParticle);\r
				protons->Add(iParticle);\r
			}\r
			else if(pdgCode == fSign*2112) // (anti)neutron\r
			{\r
				FixProductionVertex(iParticle);\r
				neutrons->Add(iParticle);\r
			}\r
		}\r
		\r
		// coalescence conditions\r
		// A.J. Baltz et al., Phys. lett B 325(1994)7\r
		\r
		TIter p_next(protons);\r
		while(TParticle* n0 = (TParticle*) p_next())\r
		{\r
			TParticle* iProton = 0;\r
			TParticle* jNeutron = 0;\r
			\r
			// Select the first nucleon\r
			iProton = n0;\r
			\r
			TVector3 v0(n0->Vx(), n0->Vy(), n0->Vz());\r
			TVector3 p0(n0->Px(), n0->Py(), n0->Pz()), p1(0,0,0);\r
			\r
			// See if there is a neibourgh...\r
			TIter n_next(neutrons);\r
			while(TParticle* n1 = (TParticle*) n_next() )\r
			{\r
				TVector3 v(n1->Vx(), n1->Vy(), n1->Vz());\r
				TVector3 p(n1->Px(), n1->Py(), n1->Pz());\r
				\r
				if((p-p0).Mag() > fPmax) continue;\r
				if((v-v0).Mag() > fRmax) continue;\r
				\r
				jNeutron = n1;\r
				p1 = p;\r
				break;\r
			}\r
			\r
			if(jNeutron == 0) continue; // with next proton\r
			\r
			if(Rndm() > fSpinProb) continue;\r
			\r
			// neutron captured!\r
			\r
			TVector3 pN = p0+p1;\r
			TVector3 vN(iProton->Vx(), iProton->Vy(), iProton->Vz()); // production vertex = p = n = collision vertex\r
			\r
			// E^2 = p^2 + m^2\r
			Double_t EN = TMath::Sqrt(pN.Mag2()+fDeuteronMass*fDeuteronMass);\r
			\r
			// Add a new (anti)deuteron to the current event stack\r
			stack->PushTrack(1, stack->Particles()->IndexOf(iProton), fSign*1000010020,\r
				         pN.X(), pN.Y(), pN.Z(), EN,\r
				         vN.X(), vN.Y(), vN.Z(), iProton->T(),\r
				         0., 0., 0., kPNCapture, ntr,1.,0);\r
			\r
			//Info("Generate","neutron capture (NEW DEUTERON)");\r
			\r
			// Set bit not to transport for the nucleons\r
			iProton->SetBit(kDoneBit);\r
			jNeutron->SetBit(kDoneBit);\r
			\r
			// Remove from the list\r
			neutrons->Remove(jNeutron);\r
		}\r
		\r
		protons->Clear("nodelete");\r
		neutrons->Clear("nodelete");\r
		\r
		delete protons;\r
		delete neutrons;\r
	}\r
	\r
	Info("Generate","Coalescence afterburner: DONE");\r
}\r
\r
// create the freeze-out nucleon distribution around the collision vertex\r
void AliGenDeuteron::FixProductionVertex(TParticle* i)\r
{\r
	Double_t x,y,z;\r
	\r
	if(fModel == kThermal) // homogeneous volume\r
	{\r
		// random (r,theta,phi)\r
		Double_t r = fRsrc*TMath::Power(Rndm(),1./3.);\r
		Double_t theta = TMath::ACos(2.*Rndm()-1.);\r
		Double_t phi = 2.*TMath::Pi()*Rndm();\r
	\r
		// transform coordenates\r
		x = r*TMath::Sin(theta)*TMath::Cos(phi);\r
		y = r*TMath::Sin(theta)*TMath::Sin(phi);\r
		z = r*TMath::Cos(theta);\r
	}\r
	else // projection into the surface of an sphere\r
	{\r
		x = fRsrc*TMath::Sin(i->Theta())*TMath::Cos(i->Phi());\r
		y = fRsrc*TMath::Sin(i->Theta())*TMath::Sin(i->Phi());\r
		z = fRsrc*TMath::Cos(i->Theta());\r
	}\r
	\r
	// assume nucleons at the freeze-out have the same production vertex\r
	i->SetProductionVertex(i->Vx()+x, i->Vy()+y, i->Vz()+z, i->T());\r
}\r
Commit	Line	Data
e3b21c34	1	/**************************************************************************\r
	2	* Copyright(c) 2009-2010, ALICE Experiment at CERN, All rights reserved. *\r
	3	* *\r
	4	* Author: The ALICE Off-line Project. *\r
	5	* Contributors are mentioned in the code where appropriate. *\r
	6	* *\r
	7	* Permission to use, copy, modify and distribute this software and its *\r
	8	* documentation strictly for non-commercial purposes is hereby granted *\r
	9	* without fee, provided that the above copyright notice appears in all *\r
	10	* copies and that both the copyright notice and this permission notice *\r
	11	* appear in the supporting documentation. The authors make no claims *\r
	12	* about the suitability of this software for any purpose. It is *\r
	13	* provided "as is" without express or implied warranty. *\r
	14	**************************************************************************/\r
	15	\r
	16	//\r
	17	// A very simple model based on the article "Physics Letters B325 (1994) 7-12".\r
	18	// The only addition is to model the freeze-out at which the\r
	19	// light nuclei can form for the energies of the LHC, since PYTHIA does not give\r
	20	// any detailed spatial description of the generated particles at the level of a\r
	21	// few fermis.\r
	22	//\r
	23	// There are two options to place the nucleons: a thermal picture where all\r
	24	// nucleons are placed randomly and homogeneously in an spherical volume of a\r
	25	// few fermis, and expansion one where they are projected on its surface.\r
	26	//\r
	27	// A (anti)deuteron will form if there is a pair of (anti)proton-(anti)neutron\r
	28	// with momentum difference less than ~ 300MeV and relative distance less than a\r
	29	// ~ 2.1fm. Only 3/4 of this clusters are accepted due to the deuteron spin.\r
	30	//\r
	31	// Author: Eulogio Serradilla <eulogio.serradilla@ciemat.es>,\r
	32	// Arturo Menchaca <menchaca@fisica.unam.mx>\r
	33	//\r
	34	\r
539d3694	35	#include "Riostream.h"\r
	36	#include "TParticle.h"\r
	37	#include "AliStack.h"\r
	38	#include "AliGenDeuteron.h"\r
	39	#include "AliGenCocktailAfterBurner.h"\r
	40	#include "TMath.h"\r
	41	#include "TMCProcess.h"\r
	42	#include "TList.h"\r
	43	#include "TVector3.h"\r
539d3694	44	#include "AliMC.h"\r
e3b21c34	45	#include "TArrayF.h"\r
e3b21c34	46	#include "AliGenCocktailEventHeader.h"\r
539d3694	47	\r
	48	ClassImp(AliGenDeuteron)\r
	49	\r
539d3694	50	AliGenDeuteron::AliGenDeuteron(Int_t sign, Double_t pmax, Double_t rmax, Double_t rsrc, Int_t model):\r
	51	fDeuteronMass(1.87561282), //pdg\r
	52	fPmax(pmax), // GeV/c\r
	53	fRmax(rmax*1.e-13), //cm\r
	54	fRsrc(rsrc*1.e-13), // cm\r
	55	fSpinProb(0.75),\r
	56	fModel(model),\r
	57	fSign(1)\r
	58	{\r
	59	fSign = sign > 0 ? 1:-1;\r
	60	}\r
	61	\r
	62	AliGenDeuteron::~AliGenDeuteron()\r
	63	{\r
	64	}\r
	65	\r
	66	void AliGenDeuteron::Init()\r
	67	{\r
	68	// Standard AliGenerator Initializer\r
	69	}\r
	70	\r
	71	void AliGenDeuteron::Generate()\r
	72	{\r
	73	// Modify the stack of each event, adding the new particles at the end and\r
	74	// not transporting the parents.\r
	75	\r
	76	Info("Generate","freeze-out model : %d (0 expansion, 1 thermal)",fModel);\r
	77	Info("Generate","relative momentum: %g GeV/c",fPmax);\r
	78	Info("Generate","relative distance: %g cm",fRmax);\r
	79	Info("Generate","source radius : %g cm",fRsrc);\r
	80	Info("Generate","spin probability : %g ",fSpinProb);\r
	81	Info("Generate","sign : %d ",fSign);\r
	82	\r
539d3694	83	Int_t ntr;\r
	84	\r
	85	// Get the cocktail generator\r
	86	AliGenCocktailAfterBurner* gener = (AliGenCocktailAfterBurner*)gAlice->GetMCApp()->Generator();\r
	87	\r
	88	// Loop over events\r
	89	for(Int_t ns = 0; ns < gener->GetNumberOfEvents(); ++ns)\r
	90	{\r
	91	gener->SetActiveEventNumber(ns);\r
	92	\r
	93	AliStack* stack = gener->GetStack(ns); // Stack of event ns\r
	94	if(!stack)\r
	95	{\r
	96	Info("Generate", "no stack, exiting");\r
	97	return;\r
	98	}\r
	99	\r
	100	// find emerging nucleons from the collision of current event\r
	101	\r
	102	TList* protons = new TList();\r
	103	protons->SetOwner(kFALSE);\r
	104	TList* neutrons = new TList();\r
	105	neutrons->SetOwner(kFALSE);\r
	106	\r
e3b21c34	107	// primary vertex of current event\r
	108	TArrayF primVtx;\r
	109	gener->GetActiveEventHeader()->PrimaryVertex(primVtx);\r
	110	TVector3 r0(primVtx[0],primVtx[1],primVtx[2]);\r
539d3694	111	\r
	112	for (Int_t i=0; i<stack->GetNtrack(); ++i)\r
	113	{\r
	114	TParticle* iParticle = stack->Particle(i);\r
	115	\r
	116	if(iParticle->TestBit(kDoneBit)) continue;\r
	117	// select only nucleons within the freeze-out volume\r
	118	TVector3 r(iParticle->Vx(),iParticle->Vy(),iParticle->Vz());\r
	119	if((r-r0).Mag()>fRsrc) continue;\r
	120	\r
	121	Int_t pdgCode = iParticle->GetPdgCode();\r
	122	if(pdgCode == fSign*2212)// (anti)proton\r
	123	{\r
e3b21c34	124	FixProductionVertex(iParticle);\r
539d3694	125	protons->Add(iParticle);\r
	126	}\r
	127	else if(pdgCode == fSign*2112) // (anti)neutron\r
	128	{\r
e3b21c34	129	FixProductionVertex(iParticle);\r
539d3694	130	neutrons->Add(iParticle);\r
	131	}\r
	132	}\r
	133	\r
	134	// coalescence conditions\r
	135	// A.J. Baltz et al., Phys. lett B 325(1994)7\r
	136	\r
	137	TIter p_next(protons);\r
	138	while(TParticle* n0 = (TParticle*) p_next())\r
	139	{\r
	140	TParticle* iProton = 0;\r
	141	TParticle* jNeutron = 0;\r
	142	\r
	143	// Select the first nucleon\r
	144	iProton = n0;\r
	145	\r
	146	TVector3 v0(n0->Vx(), n0->Vy(), n0->Vz());\r
	147	TVector3 p0(n0->Px(), n0->Py(), n0->Pz()), p1(0,0,0);\r
	148	\r
	149	// See if there is a neibourgh...\r
	150	TIter n_next(neutrons);\r
	151	while(TParticle* n1 = (TParticle*) n_next() )\r
	152	{\r
	153	TVector3 v(n1->Vx(), n1->Vy(), n1->Vz());\r
	154	TVector3 p(n1->Px(), n1->Py(), n1->Pz());\r
	155	\r
	156	if((p-p0).Mag() > fPmax) continue;\r
	157	if((v-v0).Mag() > fRmax) continue;\r
	158	\r
	159	jNeutron = n1;\r
	160	p1 = p;\r
	161	break;\r
	162	}\r
	163	\r
	164	if(jNeutron == 0) continue; // with next proton\r
	165	\r
e3b21c34	166	if(Rndm() > fSpinProb) continue;\r
539d3694	167	\r
	168	// neutron captured!\r
	169	\r
	170	TVector3 pN = p0+p1;\r
	171	TVector3 vN(iProton->Vx(), iProton->Vy(), iProton->Vz()); // production vertex = p = n = collision vertex\r
	172	\r
	173	// E^2 = p^2 + m^2\r
	174	Double_t EN = TMath::Sqrt(pN.Mag2()+fDeuteronMass*fDeuteronMass);\r
	175	\r
	176	// Add a new (anti)deuteron to the current event stack\r
	177	stack->PushTrack(1, stack->Particles()->IndexOf(iProton), fSign*1000010020,\r
	178	pN.X(), pN.Y(), pN.Z(), EN,\r
	179	vN.X(), vN.Y(), vN.Z(), iProton->T(),\r
	180	0., 0., 0., kPNCapture, ntr,1.,0);\r
	181	\r
	182	//Info("Generate","neutron capture (NEW DEUTERON)");\r
	183	\r
	184	// Set bit not to transport for the nucleons\r
	185	iProton->SetBit(kDoneBit);\r
	186	jNeutron->SetBit(kDoneBit);\r
	187	\r
	188	// Remove from the list\r
	189	neutrons->Remove(jNeutron);\r
	190	}\r
	191	\r
	192	protons->Clear("nodelete");\r
	193	neutrons->Clear("nodelete");\r
	194	\r
	195	delete protons;\r
	196	delete neutrons;\r
	197	}\r
	198	\r
	199	Info("Generate","Coalescence afterburner: DONE");\r
	200	}\r
	201	\r
	202	// create the freeze-out nucleon distribution around the collision vertex\r
e3b21c34	203	void AliGenDeuteron::FixProductionVertex(TParticle* i)\r
539d3694	204	{\r
	205	Double_t x,y,z;\r
	206	\r
	207	if(fModel == kThermal) // homogeneous volume\r
	208	{\r
	209	// random (r,theta,phi)\r
e3b21c34	210	Double_t r = fRsrc*TMath::Power(Rndm(),1./3.);\r
	211	Double_t theta = TMath::ACos(2.*Rndm()-1.);\r
	212	Double_t phi = 2.TMath::Pi()Rndm();\r
539d3694	213	\r
	214	// transform coordenates\r
	215	x = rTMath::Sin(theta)TMath::Cos(phi);\r
	216	y = rTMath::Sin(theta)TMath::Sin(phi);\r
	217	z = r*TMath::Cos(theta);\r
	218	}\r
	219	else // projection into the surface of an sphere\r
	220	{\r
	221	x = fRsrcTMath::Sin(i->Theta())TMath::Cos(i->Phi());\r
	222	y = fRsrcTMath::Sin(i->Theta())TMath::Sin(i->Phi());\r
	223	z = fRsrc*TMath::Cos(i->Theta());\r
	224	}\r
	225	\r
e3b21c34	226	// assume nucleons at the freeze-out have the same production vertex\r
539d3694	227	i->SetProductionVertex(i->Vx()+x, i->Vy()+y, i->Vz()+z, i->T());\r
539d3694	228	}\r