[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 "AliRun.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
    // Fix for the production vertex\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
4a33c50d	45	#include "AliRun.h"\r
e3b21c34	46	#include "TArrayF.h"\r
e3b21c34	47	#include "AliGenCocktailEventHeader.h"\r
539d3694	48	\r
	49	ClassImp(AliGenDeuteron)\r
	50	\r
539d3694	51	AliGenDeuteron::AliGenDeuteron(Int_t sign, Double_t pmax, Double_t rmax, Double_t rsrc, Int_t model):\r
	52	fDeuteronMass(1.87561282), //pdg\r
	53	fPmax(pmax), // GeV/c\r
	54	fRmax(rmax*1.e-13), //cm\r
	55	fRsrc(rsrc*1.e-13), // cm\r
	56	fSpinProb(0.75),\r
	57	fModel(model),\r
	58	fSign(1)\r
	59	{\r
	60	fSign = sign > 0 ? 1:-1;\r
	61	}\r
	62	\r
	63	AliGenDeuteron::~AliGenDeuteron()\r
	64	{\r
	65	}\r
	66	\r
	67	void AliGenDeuteron::Init()\r
	68	{\r
	69	// Standard AliGenerator Initializer\r
	70	}\r
	71	\r
	72	void AliGenDeuteron::Generate()\r
	73	{\r
	74	// Modify the stack of each event, adding the new particles at the end and\r
	75	// not transporting the parents.\r
	76	\r
	77	Info("Generate","freeze-out model : %d (0 expansion, 1 thermal)",fModel);\r
	78	Info("Generate","relative momentum: %g GeV/c",fPmax);\r
	79	Info("Generate","relative distance: %g cm",fRmax);\r
	80	Info("Generate","source radius : %g cm",fRsrc);\r
	81	Info("Generate","spin probability : %g ",fSpinProb);\r
	82	Info("Generate","sign : %d ",fSign);\r
	83	\r
539d3694	84	Int_t ntr;\r
	85	\r
	86	// Get the cocktail generator\r
	87	AliGenCocktailAfterBurner* gener = (AliGenCocktailAfterBurner*)gAlice->GetMCApp()->Generator();\r
	88	\r
	89	// Loop over events\r
	90	for(Int_t ns = 0; ns < gener->GetNumberOfEvents(); ++ns)\r
	91	{\r
	92	gener->SetActiveEventNumber(ns);\r
	93	\r
	94	AliStack* stack = gener->GetStack(ns); // Stack of event ns\r
	95	if(!stack)\r
	96	{\r
	97	Info("Generate", "no stack, exiting");\r
	98	return;\r
	99	}\r
	100	\r
	101	// find emerging nucleons from the collision of current event\r
	102	\r
	103	TList* protons = new TList();\r
	104	protons->SetOwner(kFALSE);\r
	105	TList* neutrons = new TList();\r
	106	neutrons->SetOwner(kFALSE);\r
	107	\r
e3b21c34	108	// primary vertex of current event\r
	109	TArrayF primVtx;\r
	110	gener->GetActiveEventHeader()->PrimaryVertex(primVtx);\r
	111	TVector3 r0(primVtx[0],primVtx[1],primVtx[2]);\r
539d3694	112	\r
	113	for (Int_t i=0; i<stack->GetNtrack(); ++i)\r
	114	{\r
	115	TParticle* iParticle = stack->Particle(i);\r
	116	\r
	117	if(iParticle->TestBit(kDoneBit)) continue;\r
	118	// select only nucleons within the freeze-out volume\r
	119	TVector3 r(iParticle->Vx(),iParticle->Vy(),iParticle->Vz());\r
	120	if((r-r0).Mag()>fRsrc) continue;\r
	121	\r
	122	Int_t pdgCode = iParticle->GetPdgCode();\r
	123	if(pdgCode == fSign*2212)// (anti)proton\r
	124	{\r
e3b21c34	125	FixProductionVertex(iParticle);\r
539d3694	126	protons->Add(iParticle);\r
	127	}\r
	128	else if(pdgCode == fSign*2112) // (anti)neutron\r
	129	{\r
e3b21c34	130	FixProductionVertex(iParticle);\r
539d3694	131	neutrons->Add(iParticle);\r
	132	}\r
	133	}\r
	134	\r
	135	// coalescence conditions\r
	136	// A.J. Baltz et al., Phys. lett B 325(1994)7\r
	137	\r
	138	TIter p_next(protons);\r
	139	while(TParticle* n0 = (TParticle*) p_next())\r
	140	{\r
	141	TParticle* iProton = 0;\r
	142	TParticle* jNeutron = 0;\r
	143	\r
	144	// Select the first nucleon\r
	145	iProton = n0;\r
	146	\r
	147	TVector3 v0(n0->Vx(), n0->Vy(), n0->Vz());\r
	148	TVector3 p0(n0->Px(), n0->Py(), n0->Pz()), p1(0,0,0);\r
	149	\r
	150	// See if there is a neibourgh...\r
	151	TIter n_next(neutrons);\r
	152	while(TParticle* n1 = (TParticle*) n_next() )\r
	153	{\r
	154	TVector3 v(n1->Vx(), n1->Vy(), n1->Vz());\r
	155	TVector3 p(n1->Px(), n1->Py(), n1->Pz());\r
	156	\r
	157	if((p-p0).Mag() > fPmax) continue;\r
	158	if((v-v0).Mag() > fRmax) continue;\r
	159	\r
	160	jNeutron = n1;\r
	161	p1 = p;\r
	162	break;\r
	163	}\r
	164	\r
	165	if(jNeutron == 0) continue; // with next proton\r
	166	\r
e3b21c34	167	if(Rndm() > fSpinProb) continue;\r
539d3694	168	\r
	169	// neutron captured!\r
	170	\r
	171	TVector3 pN = p0+p1;\r
	172	TVector3 vN(iProton->Vx(), iProton->Vy(), iProton->Vz()); // production vertex = p = n = collision vertex\r
	173	\r
	174	// E^2 = p^2 + m^2\r
	175	Double_t EN = TMath::Sqrt(pN.Mag2()+fDeuteronMass*fDeuteronMass);\r
	176	\r
	177	// Add a new (anti)deuteron to the current event stack\r
	178	stack->PushTrack(1, stack->Particles()->IndexOf(iProton), fSign*1000010020,\r
	179	pN.X(), pN.Y(), pN.Z(), EN,\r
	180	vN.X(), vN.Y(), vN.Z(), iProton->T(),\r
	181	0., 0., 0., kPNCapture, ntr,1.,0);\r
	182	\r
	183	//Info("Generate","neutron capture (NEW DEUTERON)");\r
	184	\r
	185	// Set bit not to transport for the nucleons\r
	186	iProton->SetBit(kDoneBit);\r
	187	jNeutron->SetBit(kDoneBit);\r
	188	\r
	189	// Remove from the list\r
	190	neutrons->Remove(jNeutron);\r
	191	}\r
	192	\r
	193	protons->Clear("nodelete");\r
	194	neutrons->Clear("nodelete");\r
	195	\r
	196	delete protons;\r
	197	delete neutrons;\r
	198	}\r
	199	\r
	200	Info("Generate","Coalescence afterburner: DONE");\r
	201	}\r
	202	\r
	203	// create the freeze-out nucleon distribution around the collision vertex\r
e3b21c34	204	void AliGenDeuteron::FixProductionVertex(TParticle* i)\r
539d3694	205	{\r
4a33c50d	206	// Fix for the production vertex\r
539d3694	207	Double_t x,y,z;\r
	208	\r
	209	if(fModel == kThermal) // homogeneous volume\r
	210	{\r
	211	// random (r,theta,phi)\r
e3b21c34	212	Double_t r = fRsrc*TMath::Power(Rndm(),1./3.);\r
	213	Double_t theta = TMath::ACos(2.*Rndm()-1.);\r
	214	Double_t phi = 2.TMath::Pi()Rndm();\r
539d3694	215	\r
	216	// transform coordenates\r
	217	x = rTMath::Sin(theta)TMath::Cos(phi);\r
	218	y = rTMath::Sin(theta)TMath::Sin(phi);\r
	219	z = r*TMath::Cos(theta);\r
	220	}\r
	221	else // projection into the surface of an sphere\r
	222	{\r
	223	x = fRsrcTMath::Sin(i->Theta())TMath::Cos(i->Phi());\r
	224	y = fRsrcTMath::Sin(i->Theta())TMath::Sin(i->Phi());\r
	225	z = fRsrc*TMath::Cos(i->Theta());\r
	226	}\r
	227	\r
e3b21c34	228	// assume nucleons at the freeze-out have the same production vertex\r
539d3694	229	i->SetProductionVertex(i->Vx()+x, i->Vy()+y, i->Vz()+z, i->T());\r
539d3694	230	}\r