Update timestamps for new AMANDA simulation (17/02/2015)
[u/mrichter/AliRoot.git] / EVGEN / AliGenLightNuclei.cxx
CommitLineData
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
59ClassImp(AliGenLightNuclei)
60
61AliGenLightNuclei::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
74AliGenLightNuclei::~AliGenLightNuclei()
75{
76//
77// default destructor
78//
79}
80
81void 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
170Bool_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
182TParticle* 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
201Int_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
227Int_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
261Int_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
295Int_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
335void 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
344void 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
353void 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
362void 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
371void 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
417Double_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
428Double_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}