4c039060 |
1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, 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 | $Log$ |
c0acde67 |
18 | Revision 1.55 2002/04/17 10:23:56 morsch |
19 | Coding Rule violations corrected. |
20 | |
6d4dd317 |
21 | Revision 1.54 2002/03/28 11:49:10 morsch |
22 | Pass status code in SetTrack. |
23 | |
5ec84200 |
24 | Revision 1.53 2002/03/25 14:51:13 morsch |
25 | New stack-fill and count options introduced (N. Carrer). |
26 | |
0d7c70a0 |
27 | Revision 1.51 2002/03/06 08:46:57 morsch |
28 | - Loop until np-1 |
29 | - delete dyn. alloc. arrays (N. Carrer) |
30 | |
b9817493 |
31 | Revision 1.50 2002/03/03 13:48:50 morsch |
32 | Option kPyCharmPbMNR added. Produce charm pairs in agreement with MNR |
33 | NLO calculations (Nicola Carrer). |
34 | |
fbd1348b |
35 | Revision 1.49 2002/02/08 16:50:50 morsch |
36 | Add name and title in constructor. |
37 | |
8b31bfac |
38 | Revision 1.48 2001/12/20 11:44:28 morsch |
39 | Add kinematic bias for direct gamma production. |
40 | |
b2bb450d |
41 | Revision 1.47 2001/12/19 14:45:00 morsch |
42 | Store number of trials in header. |
43 | |
fd8fb861 |
44 | Revision 1.46 2001/12/19 10:36:19 morsch |
45 | Add possibility if jet kinematic biasing. |
46 | |
988b118b |
47 | Revision 1.45 2001/11/28 08:06:52 morsch |
48 | Use fMaxLifeTime parameter. |
49 | |
8cdfd303 |
50 | Revision 1.44 2001/11/27 13:13:07 morsch |
51 | Maximum lifetime for long-lived particles to be put on the stack is parameter. |
52 | It can be set via SetMaximumLifetime(..). |
53 | |
47fc6bd5 |
54 | Revision 1.43 2001/10/21 18:35:56 hristov |
55 | Several pointers were set to zero in the default constructors to avoid memory management problems |
56 | |
2685bf00 |
57 | Revision 1.42 2001/10/15 08:21:55 morsch |
58 | Vertex truncation settings moved to AliGenMC. |
59 | |
4c07486d |
60 | Revision 1.41 2001/10/08 08:45:42 morsch |
61 | Possibility of vertex cut added. |
62 | |
a0b2b296 |
63 | Revision 1.40 2001/09/25 11:30:23 morsch |
64 | Pass event vertex to header. |
65 | |
c5d36e84 |
66 | Revision 1.39 2001/07/27 17:09:36 morsch |
67 | Use local SetTrack, KeepTrack and SetHighWaterMark methods |
68 | to delegate either to local stack or to stack owned by AliRun. |
69 | (Piotr Skowronski, A.M.) |
70 | |
a99cf51f |
71 | Revision 1.38 2001/07/13 10:58:54 morsch |
72 | - Some coded moved to AliGenMC |
73 | - Improved handling of secondary vertices. |
74 | |
fff02fee |
75 | Revision 1.37 2001/06/28 11:17:28 morsch |
76 | SetEventListRange setter added. Events in specified range are listed for |
77 | debugging. (Yuri Kharlov) |
78 | |
11dfaf8d |
79 | Revision 1.36 2001/03/30 07:05:49 morsch |
80 | Final print-out in finish run. |
81 | Write parton system for jet-production (preliminary solution). |
82 | |
9be6226f |
83 | Revision 1.35 2001/03/09 13:03:40 morsch |
84 | Process_t and Struc_Func_t moved to AliPythia.h |
85 | |
f1a48a38 |
86 | Revision 1.34 2001/02/14 15:50:40 hristov |
87 | The last particle in event marked using SetHighWaterMark |
88 | |
28adac24 |
89 | Revision 1.33 2001/01/30 09:23:12 hristov |
90 | Streamers removed (R.Brun) |
91 | |
a8a6107b |
92 | Revision 1.32 2001/01/26 19:55:51 hristov |
93 | Major upgrade of AliRoot code |
94 | |
2ab0c725 |
95 | Revision 1.31 2001/01/17 10:54:31 hristov |
96 | Better protection against FPE |
97 | |
0f9e52a3 |
98 | Revision 1.30 2000/12/18 08:55:35 morsch |
99 | Make AliPythia dependent generartors work with new scheme of random number generation |
100 | |
3356c022 |
101 | Revision 1.29 2000/12/04 11:22:03 morsch |
102 | Init of sRandom as in 1.15 |
103 | |
097ec7c1 |
104 | Revision 1.28 2000/12/02 11:41:39 morsch |
105 | Use SetRandom() to initialize random number generator in constructor. |
106 | |
1ed25e11 |
107 | Revision 1.27 2000/11/30 20:29:02 morsch |
108 | Initialise static variable sRandom in constructor: sRandom = fRandom; |
109 | |
ad1df918 |
110 | Revision 1.26 2000/11/30 07:12:50 alibrary |
111 | Introducing new Rndm and QA classes |
112 | |
65fb704d |
113 | Revision 1.25 2000/10/18 19:11:27 hristov |
114 | Division by zero fixed |
115 | |
919c2096 |
116 | Revision 1.24 2000/09/18 10:41:35 morsch |
117 | Add possibility to use nuclear structure functions from PDF library V8. |
118 | |
811826d8 |
119 | Revision 1.23 2000/09/14 14:05:40 morsch |
120 | dito |
121 | |
819f84ad |
122 | Revision 1.22 2000/09/14 14:02:22 morsch |
123 | - Correct conversion from mm to cm when passing particle vertex to MC. |
124 | - Correct handling of fForceDecay == all. |
125 | |
4451ef92 |
126 | Revision 1.21 2000/09/12 14:14:55 morsch |
127 | Call fDecayer->ForceDecay() at the beginning of Generate(). |
128 | |
00d6ce7d |
129 | Revision 1.20 2000/09/06 14:29:33 morsch |
130 | Use AliPythia for event generation an AliDecayPythia for decays. |
131 | Correct handling of "nodecay" option |
132 | |
18edb254 |
133 | Revision 1.19 2000/07/11 18:24:56 fca |
134 | Coding convention corrections + few minor bug fixes |
135 | |
aee8290b |
136 | Revision 1.18 2000/06/30 12:40:34 morsch |
137 | Pythia takes care of vertex smearing. Correct conversion from Pythia units (mm) to |
138 | Geant units (cm). |
139 | |
1512e357 |
140 | Revision 1.17 2000/06/09 20:34:07 morsch |
141 | All coding rule violations except RS3 corrected |
142 | |
f87cfe57 |
143 | Revision 1.16 2000/05/15 15:04:20 morsch |
144 | The full event is written for fNtrack = -1 |
145 | Coding rule violations corrected. |
146 | |
5ddeb374 |
147 | Revision 1.15 2000/04/26 10:14:24 morsch |
148 | Particles array has one entry more than pythia particle list. Upper bound of |
149 | particle loop changed to np-1 (R. Guernane, AM) |
150 | |
2a0e6f5b |
151 | Revision 1.14 2000/04/05 08:36:13 morsch |
152 | Check status code of particles in Pythia event |
153 | to avoid double counting as partonic state and final state particle. |
154 | |
23211d3e |
155 | Revision 1.13 1999/11/09 07:38:48 fca |
156 | Changes for compatibility with version 2.23 of ROOT |
157 | |
084c1b4a |
158 | Revision 1.12 1999/11/03 17:43:20 fca |
159 | New version from G.Martinez & A.Morsch |
160 | |
886b6f73 |
161 | Revision 1.11 1999/09/29 09:24:14 fca |
162 | Introduction of the Copyright and cvs Log |
4c039060 |
163 | */ |
164 | |
6d4dd317 |
165 | // |
166 | // Generator using the TPythia interface (via AliPythia) |
167 | // to generate pp collisions. |
168 | // Using SetNuclei() also nuclear modifications to the structure functions |
169 | // can be taken into account. This makes, of course, only sense for the |
170 | // generation of the products of hard processes (heavy flavor, jets ...) |
171 | // |
172 | // andreas.morsch@cern.ch |
173 | // |
174 | |
175 | |
fe4da5cc |
176 | #include "AliGenPythia.h" |
fff02fee |
177 | #include "AliGenPythiaEventHeader.h" |
18edb254 |
178 | #include "AliDecayerPythia.h" |
fe4da5cc |
179 | #include "AliRun.h" |
180 | #include "AliPythia.h" |
18edb254 |
181 | #include "AliPDG.h" |
1578254f |
182 | #include <TParticle.h> |
18edb254 |
183 | #include <TSystem.h> |
f87cfe57 |
184 | |
fe4da5cc |
185 | ClassImp(AliGenPythia) |
186 | |
187 | AliGenPythia::AliGenPythia() |
fff02fee |
188 | :AliGenMC() |
fe4da5cc |
189 | { |
18edb254 |
190 | // Default Constructor |
2685bf00 |
191 | fParticles = 0; |
192 | fPythia = 0; |
988b118b |
193 | fDecayer = new AliDecayerPythia(); |
11dfaf8d |
194 | SetEventListRange(); |
988b118b |
195 | SetJetPhiRange(); |
196 | SetJetEtaRange(); |
fe4da5cc |
197 | } |
198 | |
199 | AliGenPythia::AliGenPythia(Int_t npart) |
fff02fee |
200 | :AliGenMC(npart) |
fe4da5cc |
201 | { |
202 | // default charm production at 5. 5 TeV |
203 | // semimuonic decay |
204 | // structure function GRVHO |
205 | // |
8b31bfac |
206 | fName = "Pythia"; |
207 | fTitle= "Particle Generator using PYTHIA"; |
fe4da5cc |
208 | fXsection = 0.; |
811826d8 |
209 | fNucA1=0; |
210 | fNucA2=0; |
fe4da5cc |
211 | SetProcess(); |
212 | SetStrucFunc(); |
886b6f73 |
213 | SetForceDecay(); |
fe4da5cc |
214 | SetPtHard(); |
215 | SetEnergyCMS(); |
18edb254 |
216 | fDecayer = new AliDecayerPythia(); |
1ed25e11 |
217 | // Set random number generator |
097ec7c1 |
218 | sRandom=fRandom; |
fff02fee |
219 | fFlavorSelect = 0; |
220 | // Produced particles |
221 | fParticles = new TClonesArray("TParticle",1000); |
c5d36e84 |
222 | fEventVertex.Set(3); |
988b118b |
223 | SetEventListRange(); |
224 | SetJetPhiRange(); |
225 | SetJetEtaRange(); |
719b6b8a |
226 | // Options determining what to keep in the stack (Heavy flavour generation) |
227 | fStackFillOpt = kFlavorSelection; // Keep particle with selected flavor |
228 | fFeedDownOpt = kTRUE; // allow feed down from higher family |
229 | // Fragmentation on/off |
230 | fFragmentation = kTRUE; |
231 | // Default counting mode |
232 | fCountMode = kCountAll; |
fe4da5cc |
233 | } |
234 | |
f87cfe57 |
235 | AliGenPythia::AliGenPythia(const AliGenPythia & Pythia) |
236 | { |
237 | // copy constructor |
6d4dd317 |
238 | Pythia.Copy(*this); |
f87cfe57 |
239 | } |
240 | |
fe4da5cc |
241 | AliGenPythia::~AliGenPythia() |
242 | { |
5ddeb374 |
243 | // Destructor |
fe4da5cc |
244 | } |
245 | |
11dfaf8d |
246 | void AliGenPythia::SetEventListRange(Int_t eventFirst, Int_t eventLast) |
247 | { |
248 | // Set a range of event numbers, for which a table |
249 | // of generated particle will be printed |
250 | fDebugEventFirst = eventFirst; |
251 | fDebugEventLast = eventLast; |
252 | if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst; |
253 | } |
254 | |
fe4da5cc |
255 | void AliGenPythia::Init() |
256 | { |
5ddeb374 |
257 | // Initialisation |
18edb254 |
258 | SetMC(AliPythia::Instance()); |
fe4da5cc |
259 | fPythia=(AliPythia*) fgMCEvGen; |
260 | // |
261 | fParentWeight=1./Float_t(fNpart); |
262 | // |
263 | // Forward Paramters to the AliPythia object |
00d6ce7d |
264 | fDecayer->SetForceDecay(fForceDecay); |
18edb254 |
265 | fDecayer->Init(); |
00d6ce7d |
266 | |
18edb254 |
267 | |
fe4da5cc |
268 | fPythia->SetCKIN(3,fPtHardMin); |
269 | fPythia->SetCKIN(4,fPtHardMax); |
811826d8 |
270 | if (fNucA1 > 0 && fNucA2 > 0) fPythia->SetNuclei(fNucA1, fNucA2); |
719b6b8a |
271 | // Fragmentation? |
272 | if (fFragmentation) { |
273 | fPythia->SetMSTP(111,1); |
274 | } else { |
275 | fPythia->SetMSTP(111,0); |
276 | } |
277 | |
fe4da5cc |
278 | fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc); |
18edb254 |
279 | |
280 | // fPythia->Pylist(0); |
281 | // fPythia->Pystat(2); |
fe4da5cc |
282 | // Parent and Children Selection |
283 | switch (fProcess) |
284 | { |
f1a48a38 |
285 | case kPyCharm: |
f1a48a38 |
286 | case kPyCharmUnforced: |
fbd1348b |
287 | case kPyCharmPbMNR: |
b2bb450d |
288 | fParentSelect[0] = 411; |
289 | fParentSelect[1] = 421; |
290 | fParentSelect[2] = 431; |
291 | fParentSelect[3] = 4122; |
8b31bfac |
292 | fFlavorSelect = 4; |
fe4da5cc |
293 | break; |
719b6b8a |
294 | case kPyD0PbMNR: |
295 | fParentSelect[0] = 421; |
296 | fFlavorSelect = 4; |
297 | break; |
f1a48a38 |
298 | case kPyBeauty: |
c0acde67 |
299 | case kPyBeautyPbMNR: |
fff02fee |
300 | fParentSelect[0]= 511; |
301 | fParentSelect[1]= 521; |
302 | fParentSelect[2]= 531; |
303 | fParentSelect[3]= 5122; |
304 | fParentSelect[4]= 5132; |
305 | fParentSelect[5]= 5232; |
306 | fParentSelect[6]= 5332; |
307 | fFlavorSelect = 5; |
fe4da5cc |
308 | break; |
f1a48a38 |
309 | case kPyBeautyUnforced: |
fff02fee |
310 | fParentSelect[0] = 511; |
311 | fParentSelect[1] = 521; |
312 | fParentSelect[2] = 531; |
313 | fParentSelect[3] = 5122; |
314 | fParentSelect[4] = 5132; |
315 | fParentSelect[5] = 5232; |
316 | fParentSelect[6] = 5332; |
317 | fFlavorSelect = 5; |
fe4da5cc |
318 | break; |
f1a48a38 |
319 | case kPyJpsiChi: |
320 | case kPyJpsi: |
fff02fee |
321 | fParentSelect[0] = 443; |
fe4da5cc |
322 | break; |
f1a48a38 |
323 | case kPyMb: |
324 | case kPyJets: |
325 | case kPyDirectGamma: |
5be1fe76 |
326 | break; |
fe4da5cc |
327 | } |
fff02fee |
328 | AliGenMC::Init(); |
fe4da5cc |
329 | } |
330 | |
331 | void AliGenPythia::Generate() |
332 | { |
5ddeb374 |
333 | // Generate one event |
00d6ce7d |
334 | fDecayer->ForceDecay(); |
a0b2b296 |
335 | |
9be6226f |
336 | Float_t polar[3] = {0,0,0}; |
337 | Float_t origin[3] = {0,0,0}; |
fff02fee |
338 | Float_t p[3]; |
5be1fe76 |
339 | // converts from mm/c to s |
340 | const Float_t kconv=0.001/2.999792458e8; |
fe4da5cc |
341 | // |
fe4da5cc |
342 | Int_t nt=0; |
343 | Int_t jev=0; |
7921755b |
344 | Int_t j, kf; |
fe4da5cc |
345 | fTrials=0; |
fff02fee |
346 | |
fd8fb861 |
347 | // Set collision vertex position |
aee8290b |
348 | if(fVertexSmear==kPerEvent) { |
1512e357 |
349 | fPythia->SetMSTP(151,1); |
fe4da5cc |
350 | for (j=0;j<3;j++) { |
fff02fee |
351 | fPythia->SetPARP(151+j, fOsigma[j]*10.); |
fe4da5cc |
352 | } |
aee8290b |
353 | } else if (fVertexSmear==kPerTrack) { |
fe4da5cc |
354 | fPythia->SetMSTP(151,0); |
fe4da5cc |
355 | } |
fff02fee |
356 | // event loop |
fe4da5cc |
357 | while(1) |
358 | { |
084c1b4a |
359 | fPythia->Pyevnt(); |
11dfaf8d |
360 | if (gAlice->GetEvNumber()>=fDebugEventFirst && |
361 | gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(1); |
fe4da5cc |
362 | fTrials++; |
47fc6bd5 |
363 | |
fff02fee |
364 | fPythia->ImportParticles(fParticles,"All"); |
365 | |
366 | // |
367 | // |
368 | // |
369 | Int_t i; |
370 | |
371 | Int_t np = fParticles->GetEntriesFast(); |
a0b2b296 |
372 | if (np == 0 ) continue; |
373 | // Get event vertex and discard the event if the Z coord. is too big |
4c07486d |
374 | TParticle *iparticle = (TParticle *) fParticles->At(0); |
375 | Float_t distz = iparticle->Vz()/10.; |
376 | if(TMath::Abs(distz)>fCutVertexZ*fOsigma[2]) continue; |
a0b2b296 |
377 | // |
4c07486d |
378 | fEventVertex[0] = iparticle->Vx()/10.+fOrigin.At(0); |
379 | fEventVertex[1] = iparticle->Vy()/10.+fOrigin.At(1); |
380 | fEventVertex[2] = iparticle->Vz()/10.+fOrigin.At(2); |
a0b2b296 |
381 | // |
fff02fee |
382 | Int_t* pParent = new Int_t[np]; |
383 | Int_t* pSelected = new Int_t[np]; |
384 | Int_t* trackIt = new Int_t[np]; |
b9817493 |
385 | for (i=0; i< np; i++) { |
fff02fee |
386 | pParent[i] = -1; |
387 | pSelected[i] = 0; |
719b6b8a |
388 | trackIt[i] = 0; |
fff02fee |
389 | } |
719b6b8a |
390 | // printf("\n **************************************************%d\n",np); |
391 | Int_t nc = 0; // Total n. of selected particles |
392 | Int_t nParents = 0; // Selected parents |
393 | Int_t nTkbles = 0; // Trackable particles |
f1a48a38 |
394 | if (fProcess != kPyMb && fProcess != kPyJets && fProcess != kPyDirectGamma) { |
fd8fb861 |
395 | |
b9817493 |
396 | for (i = 0; i<np; i++) { |
c5d36e84 |
397 | iparticle = (TParticle *) fParticles->At(i); |
23211d3e |
398 | Int_t ks = iparticle->GetStatusCode(); |
a8228d85 |
399 | kf = CheckPDGCode(iparticle->GetPdgCode()); |
fff02fee |
400 | // No initial state partons |
2a0e6f5b |
401 | if (ks==21) continue; |
fe4da5cc |
402 | // |
fff02fee |
403 | // Heavy Flavor Selection |
5be1fe76 |
404 | // |
fff02fee |
405 | // quark ? |
406 | kf = TMath::Abs(kf); |
407 | Int_t kfl = kf; |
408 | // meson ? |
409 | if (kfl > 10) kfl/=100; |
410 | // baryon |
411 | if (kfl > 10) kfl/=10; |
412 | if (kfl > 10) kfl/=10; |
413 | |
414 | Int_t ipa = iparticle->GetFirstMother()-1; |
415 | Int_t kfMo = 0; |
416 | |
417 | if (ipa > -1) { |
418 | TParticle * mother = (TParticle *) fParticles->At(ipa); |
419 | kfMo = TMath::Abs(mother->GetPdgCode()); |
420 | } |
421 | // printf("\n particle (all) %d %d %d", i, pSelected[i], kf); |
719b6b8a |
422 | // What to keep in Stack? |
423 | Bool_t flavorOK = kFALSE; |
424 | Bool_t selectOK = kFALSE; |
425 | if (fFeedDownOpt) { |
426 | if (kfl >= fFlavorSelect) flavorOK = kTRUE; |
427 | } else { |
428 | if (kfl > fFlavorSelect) { |
429 | nc = -1; |
430 | break; |
431 | } |
432 | if (kfl == fFlavorSelect) flavorOK = kTRUE; |
433 | } |
434 | switch (fStackFillOpt) { |
435 | case kFlavorSelection: |
436 | selectOK = kTRUE; |
437 | break; |
438 | case kParentSelection: |
439 | if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE; |
440 | break; |
441 | } |
442 | if (flavorOK && selectOK) { |
5be1fe76 |
443 | // |
fff02fee |
444 | // Heavy flavor hadron or quark |
fe4da5cc |
445 | // |
fff02fee |
446 | // Kinematic seletion on final state heavy flavor mesons |
447 | if (ParentSelected(kf) && !KinematicSelection(iparticle, 0)) |
448 | { |
719b6b8a |
449 | continue; |
fff02fee |
450 | } |
451 | pSelected[i] = 1; |
719b6b8a |
452 | if (ParentSelected(kf)) ++nParents; // Update parent count |
fff02fee |
453 | // printf("\n particle (HF) %d %d %d", i, pSelected[i], kf); |
454 | } else { |
455 | // Kinematic seletion on decay products |
456 | if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf) |
457 | && !KinematicSelection(iparticle, 1)) |
458 | { |
719b6b8a |
459 | continue; |
fff02fee |
460 | } |
fe4da5cc |
461 | // |
fff02fee |
462 | // Decay products |
463 | // Select if mother was selected and is not tracked |
464 | |
465 | if (pSelected[ipa] && |
466 | !trackIt[ipa] && // mother will be tracked ? |
467 | kfMo != 5 && // mother is b-quark, don't store fragments |
468 | kfMo != 4 && // mother is c-quark, don't store fragments |
469 | kf != 92) // don't store string |
470 | { |
fe4da5cc |
471 | // |
fff02fee |
472 | // Semi-stable or de-selected: diselect decay products: |
fe4da5cc |
473 | // |
5be1fe76 |
474 | // |
47fc6bd5 |
475 | if (pSelected[i] == -1 || fDecayer->GetLifetime(kf) > fMaxLifeTime) |
fff02fee |
476 | { |
477 | Int_t ipF = iparticle->GetFirstDaughter(); |
478 | Int_t ipL = iparticle->GetLastDaughter(); |
479 | if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1; |
480 | } |
481 | // printf("\n particle (decay) %d %d %d", i, pSelected[i], kf); |
482 | pSelected[i] = (pSelected[i] == -1) ? 0 : 1; |
483 | } |
484 | } |
485 | if (pSelected[i] == -1) pSelected[i] = 0; |
486 | if (!pSelected[i]) continue; |
719b6b8a |
487 | // Count quarks only if you did not include fragmentation |
488 | if (fFragmentation && kf <= 10) continue; |
fff02fee |
489 | nc++; |
490 | // Decision on tracking |
491 | trackIt[i] = 0; |
492 | // |
493 | // Track final state particle |
494 | if (ks == 1) trackIt[i] = 1; |
495 | // Track semi-stable particles |
8cdfd303 |
496 | if ((ks ==1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime)) trackIt[i] = 1; |
fff02fee |
497 | // Track particles selected by process if undecayed. |
498 | if (fForceDecay == kNoDecay) { |
499 | if (ParentSelected(kf)) trackIt[i] = 1; |
500 | } else { |
501 | if (ParentSelected(kf)) trackIt[i] = 0; |
502 | } |
719b6b8a |
503 | if (trackIt[i] == 1) ++nTkbles; // Update trackable counter |
fff02fee |
504 | // |
505 | // |
506 | |
507 | } // particle selection loop |
719b6b8a |
508 | if (nc > 0) { |
b9817493 |
509 | for (i = 0; i<np; i++) { |
fff02fee |
510 | if (!pSelected[i]) continue; |
511 | TParticle * iparticle = (TParticle *) fParticles->At(i); |
512 | kf = CheckPDGCode(iparticle->GetPdgCode()); |
5ec84200 |
513 | Int_t ks = iparticle->GetStatusCode(); |
a0b2b296 |
514 | p[0] = iparticle->Px(); |
515 | p[1] = iparticle->Py(); |
516 | p[2] = iparticle->Pz(); |
517 | origin[0] = fOrigin[0]+iparticle->Vx()/10.; |
518 | origin[1] = fOrigin[1]+iparticle->Vy()/10.; |
519 | origin[2] = fOrigin[2]+iparticle->Vz()/10.; |
520 | Float_t tof = kconv*iparticle->T(); |
521 | Int_t ipa = iparticle->GetFirstMother()-1; |
fff02fee |
522 | Int_t iparent = (ipa > -1) ? pParent[ipa] : -1; |
a99cf51f |
523 | SetTrack(fTrackIt*trackIt[i] , |
5ec84200 |
524 | iparent, kf, p, origin, polar, tof, kPPrimary, nt, 1., ks); |
fff02fee |
525 | pParent[i] = nt; |
a99cf51f |
526 | KeepTrack(nt); |
fff02fee |
527 | } // SetTrack loop |
528 | } |
529 | } else { |
530 | nc = GenerateMB(); |
5be1fe76 |
531 | } // mb ? |
fff02fee |
532 | |
b9817493 |
533 | if (pParent) delete[] pParent; |
534 | if (pSelected) delete[] pSelected; |
535 | if (trackIt) delete[] trackIt; |
536 | |
fe4da5cc |
537 | if (nc > 0) { |
719b6b8a |
538 | switch (fCountMode) { |
539 | case kCountAll: |
540 | // printf(" Count all \n"); |
541 | jev += nc; |
542 | break; |
543 | case kCountParents: |
544 | // printf(" Count parents \n"); |
545 | jev += nParents; |
546 | break; |
547 | case kCountTrackables: |
548 | // printf(" Count trackable \n"); |
549 | jev += nTkbles; |
550 | break; |
551 | } |
5ddeb374 |
552 | if (jev >= fNpart || fNpart == -1) { |
fe4da5cc |
553 | fKineBias=Float_t(fNpart)/Float_t(fTrials); |
23211d3e |
554 | printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev); |
fd8fb861 |
555 | MakeHeader(); |
fe4da5cc |
556 | break; |
557 | } |
558 | } |
559 | } // event loop |
a99cf51f |
560 | SetHighWaterMark(nt); |
fe4da5cc |
561 | // adjust weight due to kinematic selection |
562 | AdjustWeights(); |
563 | // get cross-section |
564 | fXsection=fPythia->GetPARI(1); |
565 | } |
566 | |
fff02fee |
567 | Int_t AliGenPythia::GenerateMB() |
fe4da5cc |
568 | { |
5ec84200 |
569 | // |
570 | // Min Bias selection and other global selections |
571 | // |
fff02fee |
572 | Int_t i, kf, nt, iparent; |
573 | Int_t nc = 0; |
574 | Float_t p[3]; |
575 | Float_t polar[3] = {0,0,0}; |
576 | Float_t origin[3] = {0,0,0}; |
fff02fee |
577 | // converts from mm/c to s |
578 | const Float_t kconv=0.001/2.999792458e8; |
4451ef92 |
579 | |
fff02fee |
580 | Int_t np = fParticles->GetEntriesFast(); |
581 | Int_t* pParent = new Int_t[np]; |
b9817493 |
582 | for (i=0; i< np; i++) pParent[i] = -1; |
b2bb450d |
583 | if (fProcess == kPyJets || fProcess == kPyDirectGamma) { |
988b118b |
584 | TParticle* jet1 = (TParticle *) fParticles->At(6); |
585 | TParticle* jet2 = (TParticle *) fParticles->At(7); |
586 | if (!CheckTrigger(jet1, jet2)) return 0; |
587 | } |
588 | |
b9817493 |
589 | for (i = 0; i<np; i++) { |
fff02fee |
590 | Int_t trackIt = 0; |
591 | TParticle * iparticle = (TParticle *) fParticles->At(i); |
592 | kf = CheckPDGCode(iparticle->GetPdgCode()); |
593 | Int_t ks = iparticle->GetStatusCode(); |
594 | Int_t km = iparticle->GetFirstMother(); |
595 | // printf("\n Particle: %d %d %d", i, kf, ks); |
596 | |
597 | if ((ks == 1 && kf!=0 && KinematicSelection(iparticle, 0)) || |
598 | (ks != 1) || |
599 | (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) { |
600 | nc++; |
601 | if (ks == 1) trackIt = 1; |
602 | Int_t ipa = iparticle->GetFirstMother()-1; |
fe4da5cc |
603 | |
fff02fee |
604 | iparent = (ipa > -1) ? pParent[ipa] : -1; |
fe4da5cc |
605 | |
606 | // |
fff02fee |
607 | // store track information |
a0b2b296 |
608 | p[0] = iparticle->Px(); |
609 | p[1] = iparticle->Py(); |
610 | p[2] = iparticle->Pz(); |
611 | origin[0] = fOrigin[0]+iparticle->Vx()/10.; |
612 | origin[1] = fOrigin[1]+iparticle->Vy()/10.; |
613 | origin[2] = fOrigin[2]+iparticle->Vz()/10.; |
fff02fee |
614 | Float_t tof=kconv*iparticle->T(); |
a99cf51f |
615 | SetTrack(fTrackIt*trackIt, iparent, kf, p, origin, polar, |
5ec84200 |
616 | tof, kPPrimary, nt, 1., ks); |
a99cf51f |
617 | KeepTrack(nt); |
fff02fee |
618 | pParent[i] = nt; |
619 | } // select particle |
620 | } // particle loop |
621 | |
622 | if (pParent) delete[] pParent; |
fe4da5cc |
623 | |
fff02fee |
624 | printf("\n I've put %i particles on the stack \n",nc); |
fff02fee |
625 | return nc; |
626 | } |
fe4da5cc |
627 | |
fe4da5cc |
628 | |
fff02fee |
629 | void AliGenPythia::FinishRun() |
630 | { |
631 | // Print x-section summary |
632 | fPythia->Pystat(1); |
fe4da5cc |
633 | } |
fff02fee |
634 | |
fe4da5cc |
635 | void AliGenPythia::AdjustWeights() |
636 | { |
5ddeb374 |
637 | // Adjust the weights after generation of all events |
638 | // |
5ddeb374 |
639 | TParticle *part; |
fe4da5cc |
640 | Int_t ntrack=gAlice->GetNtrack(); |
641 | for (Int_t i=0; i<ntrack; i++) { |
2ab0c725 |
642 | part= gAlice->Particle(i); |
5ddeb374 |
643 | part->SetWeight(part->GetWeight()*fKineBias); |
fe4da5cc |
644 | } |
645 | } |
811826d8 |
646 | |
647 | void AliGenPythia::SetNuclei(Int_t a1, Int_t a2) |
648 | { |
649 | // Treat protons as inside nuclei with mass numbers a1 and a2 |
650 | fNucA1 = a1; |
651 | fNucA2 = a2; |
652 | } |
fff02fee |
653 | |
654 | |
6d4dd317 |
655 | void AliGenPythia::MakeHeader() const |
fff02fee |
656 | { |
657 | // Builds the event header, to be called after each event |
658 | AliGenEventHeader* header = new AliGenPythiaEventHeader("Pythia"); |
c5d36e84 |
659 | ((AliGenPythiaEventHeader*) header)->SetProcessType(fPythia->GetMSTI(1)); |
fd8fb861 |
660 | ((AliGenPythiaEventHeader*) header)->SetTrials(fTrials); |
c5d36e84 |
661 | header->SetPrimaryVertex(fEventVertex); |
662 | gAlice->SetGenEventHeader(header); |
fff02fee |
663 | } |
811826d8 |
664 | |
988b118b |
665 | |
5ec84200 |
666 | Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2) const |
988b118b |
667 | { |
668 | // Check the kinematic trigger condition |
669 | // |
b2bb450d |
670 | Double_t eta[2]; |
671 | eta[0] = jet1->Eta(); |
672 | eta[1] = jet2->Eta(); |
673 | Double_t phi[2]; |
674 | phi[0] = jet1->Phi(); |
675 | phi[1] = jet2->Phi(); |
676 | Int_t pdg[2]; |
677 | pdg[0] = jet1->GetPdgCode(); |
678 | pdg[1] = jet2->GetPdgCode(); |
988b118b |
679 | Bool_t triggered = kFALSE; |
b2bb450d |
680 | |
681 | if (fProcess == kPyJets) { |
682 | //Check eta range first... |
683 | if ((eta[0] < fEtaMaxJet && eta[0] > fEtaMinJet) || |
684 | (eta[1] < fEtaMaxJet && eta[1] > fEtaMinJet)) |
685 | { |
686 | //Eta is okay, now check phi range |
687 | if ((phi[0] < fPhiMaxJet && phi[0] > fPhiMinJet) || |
688 | (phi[1] < fPhiMaxJet && phi[1] > fPhiMinJet)) |
689 | { |
690 | triggered = kTRUE; |
691 | } |
692 | } |
693 | } else { |
694 | Int_t ij = 0; |
695 | Int_t ig = 1; |
696 | if (pdg[0] == kGamma) { |
697 | ij = 1; |
698 | ig = 0; |
699 | } |
700 | //Check eta range first... |
701 | if ((eta[ij] < fEtaMaxJet && eta[ij] > fEtaMinJet) && |
702 | (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma)) |
703 | { |
704 | //Eta is okay, now check phi range |
705 | if ((phi[ij] < fPhiMaxJet && phi[ij] > fPhiMinJet) && |
706 | (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma)) |
707 | { |
708 | triggered = kTRUE; |
709 | } |
710 | } |
988b118b |
711 | } |
b2bb450d |
712 | |
713 | |
988b118b |
714 | return triggered; |
715 | } |
811826d8 |
716 | |
f87cfe57 |
717 | AliGenPythia& AliGenPythia::operator=(const AliGenPythia& rhs) |
718 | { |
719 | // Assignment operator |
720 | return *this; |
721 | } |
fe4da5cc |
722 | |
723 | |
811826d8 |
724 | |
a8a6107b |
725 | #ifdef never |
18edb254 |
726 | void AliGenPythia::Streamer(TBuffer &R__b) |
727 | { |
728 | // Stream an object of class AliGenPythia. |
729 | |
730 | if (R__b.IsReading()) { |
731 | Version_t R__v = R__b.ReadVersion(); if (R__v) { } |
732 | AliGenerator::Streamer(R__b); |
733 | R__b >> (Int_t&)fProcess; |
734 | R__b >> (Int_t&)fStrucFunc; |
735 | R__b >> (Int_t&)fForceDecay; |
736 | R__b >> fEnergyCMS; |
737 | R__b >> fKineBias; |
738 | R__b >> fTrials; |
739 | fParentSelect.Streamer(R__b); |
740 | fChildSelect.Streamer(R__b); |
741 | R__b >> fXsection; |
742 | // (AliPythia::Instance())->Streamer(R__b); |
743 | R__b >> fPtHardMin; |
744 | R__b >> fPtHardMax; |
745 | // if (fDecayer) fDecayer->Streamer(R__b); |
746 | } else { |
747 | R__b.WriteVersion(AliGenPythia::IsA()); |
748 | AliGenerator::Streamer(R__b); |
749 | R__b << (Int_t)fProcess; |
750 | R__b << (Int_t)fStrucFunc; |
751 | R__b << (Int_t)fForceDecay; |
752 | R__b << fEnergyCMS; |
753 | R__b << fKineBias; |
754 | R__b << fTrials; |
755 | fParentSelect.Streamer(R__b); |
756 | fChildSelect.Streamer(R__b); |
757 | R__b << fXsection; |
758 | // R__b << fPythia; |
759 | R__b << fPtHardMin; |
760 | R__b << fPtHardMax; |
761 | // fDecayer->Streamer(R__b); |
762 | } |
763 | } |
a8a6107b |
764 | #endif |
18edb254 |
765 | |