Line data Source code
1 : /**************************************************************************
2 : * Copyright(c) 2009-2015, 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 : // whose momenta in their CM frame is less than p0 (coalescence momentum).
23 : //
24 : // Sample code for PYTHIA:
25 : //
26 : // AliGenLightNuclei* gener = new AliGenLightNuclei();
27 : //
28 : // AliGenPythia* pythia = new AliGenPythia(-1);
29 : // pythia->SetCollisionSystem("p+", "p+");
30 : // pythia->SetEnergyCMS(7000);
31 : //
32 : // gener->UsePerEventRates();
33 : // gener->AddGenerator(pythia, "PYTHIA", 1);
34 : // gener->SetNucleusPdgCode(AliGenLightNuclei::kDeuteron); // default
35 : // gener->SetCoalescenceMomentum(0.100); // default (GeV/c)
36 : // gener->SetSpinProbability(0.75); // default: 1
37 : //
38 : //////////////////////////////////////////////////////////////////////
39 :
40 : //
41 : // Author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
42 : //
43 :
44 : #include "TMath.h"
45 : #include "TPDGCode.h"
46 : #include "TMCProcess.h"
47 : #include "TList.h"
48 : #include "TVector3.h"
49 : #include "TParticle.h"
50 : #include "AliStack.h"
51 : #include "AliRun.h"
52 : #include "AliLog.h"
53 : #include "AliMC.h"
54 : #include "TMath.h"
55 : #include "TLorentzVector.h"
56 : #include "AliGenLightNuclei.h"
57 :
58 6 : ClassImp(AliGenLightNuclei)
59 :
60 : AliGenLightNuclei::AliGenLightNuclei()
61 0 : :AliGenCocktail()
62 0 : ,fPdg(kDeuteron)
63 0 : ,fP0(0.100)
64 0 : ,fSpinProb(1)
65 0 : {
66 : //
67 : // default constructor
68 : //
69 0 : }
70 :
71 0 : AliGenLightNuclei::~AliGenLightNuclei()
72 0 : {
73 : //
74 : // default destructor
75 : //
76 0 : }
77 :
78 : void AliGenLightNuclei::Generate()
79 : {
80 : //
81 : // delegate the particle generation to the cocktail
82 : // and modify the stack adding the light nuclei
83 : //
84 0 : AliGenCocktail::Generate();
85 :
86 : // find the nucleons and anti-nucleons
87 0 : TList* protons = new TList();
88 0 : TList* neutrons = new TList();
89 0 : TList* lambdas = new TList();
90 0 : TList* antiprotons = new TList();
91 0 : TList* antineutrons = new TList();
92 0 : TList* antilambdas = new TList();
93 :
94 0 : for (Int_t i=0; i < fStack->GetNprimary(); ++i)
95 : {
96 0 : TParticle* iParticle = fStack->Particle(i);
97 :
98 0 : if(iParticle->GetStatusCode() != 1) continue;
99 :
100 0 : switch(iParticle->GetPdgCode())
101 : {
102 : case kProton:
103 0 : protons->Add(iParticle);
104 0 : break;
105 : case kNeutron:
106 0 : neutrons->Add(iParticle);
107 0 : break;
108 : case kLambda0:
109 0 : lambdas->Add(iParticle);
110 0 : break;
111 : case kProtonBar:
112 0 : antiprotons->Add(iParticle);
113 0 : break;
114 : case kNeutronBar:
115 0 : antineutrons->Add(iParticle);
116 0 : break;
117 : case kLambda0Bar:
118 0 : antilambdas->Add(iParticle);
119 0 : break;
120 : default:
121 : break;
122 : }
123 0 : }
124 :
125 : // do not delete content
126 0 : protons->SetOwner(kFALSE);
127 0 : neutrons->SetOwner(kFALSE);
128 0 : lambdas->SetOwner(kFALSE);
129 0 : antiprotons->SetOwner(kFALSE);
130 0 : antineutrons->SetOwner(kFALSE);
131 0 : antilambdas->SetOwner(kFALSE);
132 :
133 0 : if(TMath::Abs(fPdg) == kDeuteron)
134 : {
135 0 : this->GenerateNuclei(kDeuteron, 1.87561282, protons, neutrons);
136 0 : this->GenerateNuclei(-kDeuteron, 1.87561282, antiprotons, antineutrons);
137 0 : }
138 0 : else if(TMath::Abs(fPdg) == kTriton)
139 : {
140 0 : this->GenerateNuclei(kTriton, 2.80925, protons, neutrons, neutrons);
141 0 : this->GenerateNuclei(-kTriton, 2.80925, antiprotons, antineutrons, antineutrons);
142 0 : }
143 0 : else if(TMath::Abs(fPdg) == kHyperTriton )
144 : {
145 0 : this->GenerateNuclei(kHyperTriton, 2.99131, protons, neutrons, lambdas);
146 0 : this->GenerateNuclei(-kHyperTriton, 2.99131, antiprotons, antineutrons, antilambdas);
147 0 : }
148 0 : else if(TMath::Abs(fPdg) == kHe3Nucleus)
149 : {
150 0 : this->GenerateNuclei(kHe3Nucleus, 2.80923, protons, neutrons, protons);
151 0 : this->GenerateNuclei(-kHe3Nucleus, 2.80923, antiprotons, antineutrons, antiprotons);
152 0 : }
153 0 : else if(TMath::Abs(fPdg) == kAlpha )
154 : {
155 0 : this->GenerateNuclei(kAlpha, 3.727417, protons, neutrons, protons, neutrons);
156 0 : this->GenerateNuclei(-kAlpha, 3.727417, antiprotons, antineutrons, antiprotons, antineutrons);
157 0 : }
158 : else
159 : {
160 0 : AliFatal(Form("Unknown nucleus PDG: %d", fPdg));
161 : }
162 :
163 0 : protons->Clear("nodelete");
164 0 : neutrons->Clear("nodelete");
165 0 : lambdas->Clear("nodelete");
166 0 : antiprotons->Clear("nodelete");
167 0 : antineutrons->Clear("nodelete");
168 0 : antilambdas->Clear("nodelete");
169 :
170 0 : delete protons;
171 0 : delete neutrons;
172 0 : delete lambdas;
173 0 : delete antiprotons;
174 0 : delete antineutrons;
175 0 : delete antilambdas;
176 0 : }
177 :
178 : Bool_t AliGenLightNuclei::Coalescence(const TLorentzVector& p1, const TLorentzVector& p2) const
179 : {
180 : //
181 : // returns true if the 2 nucleons have momentum < p0 in their CM frame
182 : //
183 0 : TLorentzVector p1cm(p1);
184 0 : TLorentzVector p2cm(p2);
185 :
186 0 : TLorentzVector p = p1 + p2;
187 0 : TVector3 b = -p.BoostVector();
188 :
189 0 : p1cm.Boost(b);
190 0 : p2cm.Boost(b);
191 :
192 0 : return (p1cm.Vect().Mag() < fP0) && (p2cm.Vect().Mag() < fP0) && (Rndm() <= fSpinProb);
193 0 : }
194 :
195 : Bool_t AliGenLightNuclei::Coalescence(const TLorentzVector& p1, const TLorentzVector& p2, const TLorentzVector& p3) const
196 : {
197 : //
198 : // returns true if the 3 nucleons have momentum < p0 in their CM frame
199 : //
200 0 : TLorentzVector p1cm(p1);
201 0 : TLorentzVector p2cm(p2);
202 0 : TLorentzVector p3cm(p3);
203 :
204 0 : TLorentzVector p = p1 + p2 + p3;
205 0 : TVector3 b = -p.BoostVector();
206 :
207 0 : p1cm.Boost(b);
208 0 : p2cm.Boost(b);
209 0 : p3cm.Boost(b);
210 :
211 0 : return (p1cm.Vect().Mag() < fP0) && (p2cm.Vect().Mag() < fP0) && (p3cm.Vect().Mag() < fP0) && (Rndm() <= fSpinProb);
212 0 : }
213 :
214 : Bool_t AliGenLightNuclei::Coalescence(const TLorentzVector& p1, const TLorentzVector& p2, const TLorentzVector& p3, const TLorentzVector& p4) const
215 : {
216 : //
217 : // returns true if the 4 nucleons have momentum < p0 in their CM frame
218 : //
219 0 : TLorentzVector p1cm(p1);
220 0 : TLorentzVector p2cm(p2);
221 0 : TLorentzVector p3cm(p3);
222 0 : TLorentzVector p4cm(p4);
223 :
224 0 : TLorentzVector p = p1 + p2 + p3 + p4;
225 0 : TVector3 b = -p.BoostVector();
226 :
227 0 : p1cm.Boost(b);
228 0 : p2cm.Boost(b);
229 0 : p3cm.Boost(b);
230 0 : p4cm.Boost(b);
231 :
232 0 : return (p1cm.Vect().Mag() < fP0) && (p2cm.Vect().Mag() < fP0) && (p3cm.Vect().Mag() < fP0) && (p4cm.Vect().Mag() < fP0) && (Rndm() <= fSpinProb);
233 0 : }
234 :
235 : Int_t AliGenLightNuclei::GenerateNuclei(Int_t pdg, Double_t mass, const TList* l1, const TList* l2)
236 : {
237 : //
238 : // a nucleus with A=2 is generated from the first n1-n2 nucleon cluster
239 : // that fulfill the coalescence condition
240 : //
241 : Int_t npart = 0;
242 :
243 0 : TIter n1iter(l1);
244 0 : while(TParticle* n1 = dynamic_cast<TParticle*>(n1iter()))
245 : {
246 0 : if(n1 == 0) continue;
247 0 : if(n1->GetStatusCode() == kCluster) continue;
248 :
249 0 : TLorentzVector p1(n1->Px(), n1->Py(), n1->Pz(), n1->Energy());
250 :
251 0 : TIter n2iter(l2);
252 0 : if(l2 == l1) n2iter = n1iter;
253 :
254 0 : while(TParticle* n2 = dynamic_cast<TParticle*>( n2iter()))
255 : {
256 0 : if(n2 == 0) continue;
257 0 : if(n2 == n1) continue;
258 0 : if(n2->GetStatusCode() == kCluster) continue;
259 :
260 0 : TLorentzVector p2(n2->Px(), n2->Py(), n2->Pz(), n2->Energy());
261 :
262 0 : if(!this->Coalescence(p1, p2)) continue;
263 :
264 0 : this->PushNucleus(pdg, mass, n1, n2);
265 :
266 0 : ++npart;
267 :
268 0 : break;
269 0 : }
270 0 : }
271 :
272 : return npart;
273 0 : }
274 :
275 : Int_t AliGenLightNuclei::GenerateNuclei(Int_t pdg, Double_t mass, const TList* l1, const TList* l2, const TList* l3)
276 : {
277 : //
278 : // a nucleus with A=3 is generated from the first n1-n2-n3 nucleon cluster
279 : // that fulfill the coalescence condition
280 : //
281 : Int_t npart = 0;
282 :
283 0 : TIter n1iter(l1);
284 0 : while(TParticle* n1 = dynamic_cast<TParticle*>(n1iter()))
285 : {
286 0 : if(n1 == 0) continue;
287 0 : if(n1->GetStatusCode() == kCluster) continue;
288 :
289 0 : TLorentzVector p1(n1->Px(), n1->Py(), n1->Pz(), n1->Energy());
290 :
291 0 : TIter n2iter(l2);
292 0 : if(l2 == l1) n2iter = n1iter;
293 :
294 0 : while(TParticle* n2 = dynamic_cast<TParticle*>( n2iter()) )
295 : {
296 0 : if(n2 == 0) continue;
297 0 : if(n2 == n1) continue;
298 0 : if(n2->GetStatusCode() == kCluster) continue;
299 :
300 0 : TLorentzVector p2(n2->Px(), n2->Py(), n2->Pz(), n2->Energy());
301 :
302 0 : TIter n3iter(l3);
303 0 : if(l3 == l1) n3iter = n1iter;
304 0 : if(l3 == l2) n3iter = n2iter;
305 :
306 0 : while(TParticle* n3 = dynamic_cast<TParticle*>( n3iter()) )
307 : {
308 0 : if(n3 == 0) continue;
309 0 : if(n3 == n1) continue;
310 0 : if(n3 == n2) continue;
311 0 : if(n3->GetStatusCode() == kCluster) continue;
312 :
313 0 : TLorentzVector p3(n3->Px(), n3->Py(), n3->Pz(), n3->Energy());
314 :
315 0 : if(!this->Coalescence(p1, p2, p3)) continue;
316 :
317 0 : this->PushNucleus(pdg, mass, n1, n2, n3);
318 :
319 0 : ++npart;
320 :
321 0 : break;
322 0 : }
323 :
324 0 : if(n2->GetStatusCode() == kCluster) break;
325 0 : }
326 0 : }
327 :
328 : return npart;
329 0 : }
330 :
331 : Int_t AliGenLightNuclei::GenerateNuclei(Int_t pdg, Double_t mass, const TList* l1, const TList* l2, const TList* l3, const TList* l4)
332 : {
333 : //
334 : // a nucleus with A=4 is generated from the first n1-n2-n3-n4 nucleon cluster
335 : // that fulfill the coalescence condition
336 : //
337 : Int_t npart = 0;
338 :
339 0 : TIter n1iter(l1);
340 0 : while(TParticle* n1 = dynamic_cast<TParticle*>(n1iter()))
341 : {
342 0 : if(n1 == 0) continue;
343 0 : if(n1->GetStatusCode() == kCluster) continue;
344 :
345 0 : TLorentzVector p1(n1->Px(), n1->Py(), n1->Pz(), n1->Energy());
346 :
347 0 : TIter n2iter(l2);
348 0 : if(l2 == l1) n2iter = n1iter;
349 :
350 0 : while(TParticle* n2 = dynamic_cast<TParticle*>( n2iter()))
351 : {
352 0 : if(n2 == 0) continue;
353 0 : if(n2->GetStatusCode() == kCluster) continue;
354 :
355 0 : TLorentzVector p2(n2->Px(), n2->Py(), n2->Pz(), n2->Energy());
356 :
357 0 : TIter n3iter(l3);
358 0 : if(l3 == l1) n3iter = n1iter;
359 0 : if(l3 == l2) n3iter = n2iter;
360 :
361 0 : while(TParticle* n3 = dynamic_cast<TParticle*>( n3iter()))
362 : {
363 0 : if(n3 == 0) continue;
364 0 : if(n3 == n1) continue;
365 0 : if(n3 == n2) continue;
366 0 : if(n3->GetStatusCode() == kCluster) continue;
367 :
368 0 : TLorentzVector p3(n3->Px(), n3->Py(), n3->Pz(), n3->Energy());
369 :
370 0 : TIter n4iter(l4);
371 0 : if(l4 == l1) n4iter = n1iter;
372 0 : if(l4 == l2) n4iter = n2iter;
373 0 : if(l4 == l3) n4iter = n3iter;
374 :
375 0 : while(TParticle* n4 = dynamic_cast<TParticle*>( n4iter()))
376 : {
377 0 : if(n4 == 0) continue;
378 0 : if(n4 == n1) continue;
379 0 : if(n4 == n2) continue;
380 0 : if(n4 == n3) continue;
381 0 : if(n4->GetStatusCode() == kCluster) continue;
382 :
383 0 : TLorentzVector p4(n4->Px(), n4->Py(), n4->Pz(), n4->Energy());
384 :
385 0 : if(!this->Coalescence(p1, p2, p3, p4)) continue;
386 :
387 0 : this->PushNucleus(pdg, mass, n1, n2, n3, n4);
388 :
389 0 : ++npart;
390 :
391 0 : break;
392 0 : }
393 :
394 0 : if(n3->GetStatusCode() == kCluster) break;
395 0 : }
396 :
397 0 : if(n2->GetStatusCode() == kCluster) break;
398 0 : }
399 0 : }
400 :
401 : return npart;
402 0 : }
403 :
404 : void AliGenLightNuclei::PushNucleus(Int_t pdg, Double_t mass, TParticle* parent1, TParticle* parent2, TParticle* parent3, TParticle* parent4)
405 : {
406 : //
407 : // push a nucleus to the stack and tag the parents with the kCluster status code
408 : //
409 0 : Int_t ntrk;
410 :
411 : // momentum
412 0 : TVector3 p1(parent1->Px(), parent1->Py(), parent1->Pz());
413 0 : TVector3 p2(parent2->Px(), parent2->Py(), parent2->Pz());
414 0 : TVector3 p3(0, 0, 0);
415 0 : TVector3 p4(0, 0, 0);
416 0 : if(parent3 != 0) p3.SetXYZ(parent3->Px(), parent3->Py(), parent3->Pz());
417 0 : if(parent4 != 0) p4.SetXYZ(parent4->Px(), parent4->Py(), parent4->Pz());
418 :
419 : // momentum
420 0 : TVector3 pN = p1+p2+p3+p4;
421 :
422 : // E^2 = p^2 + m^2
423 0 : Double_t energy = TMath::Sqrt(pN.Mag2() + mass*mass);
424 :
425 : // production vertex same as the parent1'
426 0 : TVector3 vN(parent1->Vx(), parent1->Vy(), parent1->Vz());
427 :
428 : Double_t weight = 1;
429 : Int_t is = 1; // final state particle
430 :
431 : // add a new nucleus to current event stack
432 0 : fStack->PushTrack(1, -1, pdg,
433 0 : pN.X(), pN.Y(), pN.Z(), energy,
434 0 : vN.X(), vN.Y(), vN.Z(), parent1->T(),
435 : 0., 0., 0., kPNCapture, ntrk, weight, is);
436 :
437 : // change the status code of the parents
438 0 : parent1->SetStatusCode(kCluster);
439 0 : parent2->SetStatusCode(kCluster);
440 0 : if(parent3 != 0) parent3->SetStatusCode(kCluster);
441 0 : if(parent4 != 0) parent4->SetStatusCode(kCluster);
442 :
443 : // set no transport for the parents
444 0 : parent1->SetBit(kDoneBit);
445 0 : parent2->SetBit(kDoneBit);
446 0 : if(parent3 != 0) parent3->SetBit(kDoneBit);
447 0 : if(parent4 != 0) parent4->SetBit(kDoneBit);
448 0 : }
|
