Coverity fix.
[u/mrichter/AliRoot.git] / EVGEN / AliGenDeuteron.cxx
CommitLineData
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
46#include "AliGenCocktailEventHeader.h"\r
539d3694 47\r
48ClassImp(AliGenDeuteron)\r
49\r
539d3694 50AliGenDeuteron::AliGenDeuteron(Int_t sign, Double_t pmax, Double_t rmax, Double_t rsrc, Int_t model):\r
51fDeuteronMass(1.87561282), //pdg\r
52fPmax(pmax), // GeV/c\r
53fRmax(rmax*1.e-13), //cm\r
54fRsrc(rsrc*1.e-13), // cm\r
55fSpinProb(0.75),\r
56fModel(model),\r
57fSign(1)\r
58{\r
59 fSign = sign > 0 ? 1:-1;\r
60}\r
61\r
62AliGenDeuteron::~AliGenDeuteron()\r
63{\r
64}\r
65\r
66void AliGenDeuteron::Init()\r
67{\r
68 // Standard AliGenerator Initializer\r
69}\r
70\r
71void 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 203void 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 = r*TMath::Sin(theta)*TMath::Cos(phi);\r
216 y = r*TMath::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 = fRsrc*TMath::Sin(i->Theta())*TMath::Cos(i->Phi());\r
222 y = fRsrc*TMath::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
228}\r