]>
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(), |
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( p1x*p1x + p1y*p1y + p1z*p1z + m1*m1); | |
423 | Double_t E2 = TMath::Sqrt( p2x*p2x + p2y*p2y + p2z*p2z + m2*m2); | |
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 | } |