SetSeed implementation
[u/mrichter/AliRoot.git] / EVGEN / AliGenLightNuclei.cxx
1 /**************************************************************************
2  * Copyright(c) 2009-2014, ALICE Experiment at CERN, All rights reserved. *
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,
26 // then He3 nuclei, tritons and finally deuterons. It can also generate
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);
39 //    gener->SetCoalescenceMomentum(0.100); // default (GeV/c)
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()
63 ,fP0(0.100)
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 //
176         Double_t deltaP = this->GetPcm(n1->Px(), n1->Py(), n1->Pz(), n1->GetMass(), 
177                                        n2->Px(), n2->Py(), n2->Pz(), n2->GetMass());
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 }