39d810c8 |
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 | /* $Id$ */ |
17 | |
18 | // |
19 | // Generator using the TPythia interface (via AliPythia) |
20 | // to generate pp collisions. |
21 | // Using SetNuclei() also nuclear modifications to the structure functions |
22 | // can be taken into account. This makes, of course, only sense for the |
23 | // generation of the products of hard processes (heavy flavor, jets ...) |
24 | // |
25 | // andreas.morsch@cern.ch |
26 | // |
27 | |
28 | #include <TClonesArray.h> |
29 | #include <TDatabasePDG.h> |
30 | #include <TParticle.h> |
31 | #include <TPDGCode.h> |
32 | #include <TSystem.h> |
33 | #include <TTree.h> |
34 | #include "AliConst.h" |
35 | #include "AliDecayerPythia.h" |
36 | #include "AliGenPythiaPlus.h" |
37 | #include "AliHeader.h" |
38 | #include "AliGenPythiaEventHeader.h" |
39 | #include "AliPythiaBase.h" |
40 | #include "AliPythiaRndm.h" |
41 | #include "AliRun.h" |
42 | #include "AliStack.h" |
43 | #include "AliRunLoader.h" |
44 | #include "AliMC.h" |
45 | #include "PyquenCommon.h" |
46 | |
47 | ClassImp(AliGenPythiaPlus) |
48 | |
49 | |
50 | AliGenPythiaPlus::AliGenPythiaPlus(): |
51 | AliGenMC(), |
52 | fPythia(0), |
53 | fProcess(kPyCharm), |
54 | fStrucFunc(kCTEQ5L), |
55 | fEnergyCMS(5500.), |
56 | fKineBias(0.), |
57 | fTrials(0), |
58 | fTrialsRun(0), |
59 | fQ(0.), |
60 | fX1(0.), |
61 | fX2(0.), |
62 | fEventTime(0.), |
63 | fInteractionRate(0.), |
64 | fTimeWindow(0.), |
65 | fCurSubEvent(0), |
66 | fEventsTime(0), |
67 | fNev(0), |
68 | fFlavorSelect(0), |
69 | fXsection(0.), |
70 | fPtHardMin(0.), |
71 | fPtHardMax(1.e4), |
72 | fYHardMin(-1.e10), |
73 | fYHardMax(1.e10), |
74 | fGinit(1), |
75 | fGfinal(1), |
76 | fHadronisation(1), |
77 | fNpartons(0), |
78 | fReadFromFile(0), |
79 | fQuench(0), |
80 | fPtKick(1.), |
81 | fFullEvent(kTRUE), |
82 | fDecayer(new AliDecayerPythia()), |
83 | fDebugEventFirst(-1), |
84 | fDebugEventLast(-1), |
85 | fEtMinJet(0.), |
86 | fEtMaxJet(1.e4), |
87 | fEtaMinJet(-20.), |
88 | fEtaMaxJet(20.), |
89 | fPhiMinJet(0.), |
90 | fPhiMaxJet(2.* TMath::Pi()), |
91 | fJetReconstruction(kCell), |
92 | fEtaMinGamma(-20.), |
93 | fEtaMaxGamma(20.), |
94 | fPhiMinGamma(0.), |
95 | fPhiMaxGamma(2. * TMath::Pi()), |
96 | fPycellEtaMax(2.), |
97 | fPycellNEta(274), |
98 | fPycellNPhi(432), |
99 | fPycellThreshold(0.), |
100 | fPycellEtSeed(4.), |
101 | fPycellMinEtJet(10.), |
102 | fPycellMaxRadius(1.), |
103 | fStackFillOpt(kFlavorSelection), |
104 | fFeedDownOpt(kTRUE), |
105 | fFragmentation(kTRUE), |
106 | fSetNuclei(kFALSE), |
107 | fNewMIS(kFALSE), |
108 | fHFoff(kFALSE), |
109 | fTriggerParticle(0), |
110 | fTriggerEta(0.9), |
111 | fCountMode(kCountAll), |
112 | fHeader(0), |
113 | fRL(0), |
114 | fFileName(0), |
115 | fFragPhotonInCalo(kFALSE), |
116 | fPi0InCalo(kFALSE) , |
117 | fPhotonInCalo(kFALSE), |
118 | fCheckEMCAL(kFALSE), |
119 | fCheckPHOS(kFALSE), |
120 | fCheckPHOSeta(kFALSE), |
121 | fFragPhotonOrPi0MinPt(0), |
122 | fPhotonMinPt(0), |
123 | fPHOSMinPhi(219.), |
124 | fPHOSMaxPhi(321.), |
125 | fPHOSEta(0.13), |
126 | fEMCALMinPhi(79.), |
127 | fEMCALMaxPhi(191.), |
128 | fEMCALEta(0.71) |
129 | |
130 | { |
131 | // Default Constructor |
132 | SetNuclei(0,0); |
133 | if (!AliPythiaRndm::GetPythiaRandom()) |
134 | AliPythiaRndm::SetPythiaRandom(GetRandom()); |
135 | } |
136 | |
137 | AliGenPythiaPlus::AliGenPythiaPlus(AliPythiaBase* pythia) |
138 | :AliGenMC(-1), |
139 | fPythia(pythia), |
140 | fProcess(kPyCharm), |
141 | fStrucFunc(kCTEQ5L), |
142 | fEnergyCMS(5500.), |
143 | fKineBias(0.), |
144 | fTrials(0), |
145 | fTrialsRun(0), |
146 | fQ(0.), |
147 | fX1(0.), |
148 | fX2(0.), |
149 | fEventTime(0.), |
150 | fInteractionRate(0.), |
151 | fTimeWindow(0.), |
152 | fCurSubEvent(0), |
153 | fEventsTime(0), |
154 | fNev(0), |
155 | fFlavorSelect(0), |
156 | fXsection(0.), |
157 | fPtHardMin(0.), |
158 | fPtHardMax(1.e4), |
159 | fYHardMin(-1.e10), |
160 | fYHardMax(1.e10), |
161 | fGinit(kTRUE), |
162 | fGfinal(kTRUE), |
163 | fHadronisation(kTRUE), |
164 | fNpartons(0), |
165 | fReadFromFile(kFALSE), |
166 | fQuench(kFALSE), |
167 | fPtKick(1.), |
168 | fFullEvent(kTRUE), |
169 | fDecayer(new AliDecayerPythia()), |
170 | fDebugEventFirst(-1), |
171 | fDebugEventLast(-1), |
172 | fEtMinJet(0.), |
173 | fEtMaxJet(1.e4), |
174 | fEtaMinJet(-20.), |
175 | fEtaMaxJet(20.), |
176 | fPhiMinJet(0.), |
177 | fPhiMaxJet(2.* TMath::Pi()), |
178 | fJetReconstruction(kCell), |
179 | fEtaMinGamma(-20.), |
180 | fEtaMaxGamma(20.), |
181 | fPhiMinGamma(0.), |
182 | fPhiMaxGamma(2. * TMath::Pi()), |
183 | fPycellEtaMax(2.), |
184 | fPycellNEta(274), |
185 | fPycellNPhi(432), |
186 | fPycellThreshold(0.), |
187 | fPycellEtSeed(4.), |
188 | fPycellMinEtJet(10.), |
189 | fPycellMaxRadius(1.), |
190 | fStackFillOpt(kFlavorSelection), |
191 | fFeedDownOpt(kTRUE), |
192 | fFragmentation(kTRUE), |
193 | fSetNuclei(kFALSE), |
194 | fNewMIS(kFALSE), |
195 | fHFoff(kFALSE), |
196 | fTriggerParticle(0), |
197 | fTriggerEta(0.9), |
198 | fCountMode(kCountAll), |
199 | fHeader(0), |
200 | fRL(0), |
201 | fFileName(0), |
202 | fFragPhotonInCalo(kFALSE), |
203 | fPi0InCalo(kFALSE) , |
204 | fPhotonInCalo(kFALSE), |
205 | fCheckEMCAL(kFALSE), |
206 | fCheckPHOS(kFALSE), |
207 | fCheckPHOSeta(kFALSE), |
208 | fFragPhotonOrPi0MinPt(0), |
209 | fPhotonMinPt(0), |
210 | fPHOSMinPhi(219.), |
211 | fPHOSMaxPhi(321.), |
212 | fPHOSEta(0.13), |
213 | fEMCALMinPhi(79.), |
214 | fEMCALMaxPhi(191.), |
215 | fEMCALEta(0.71) |
216 | { |
217 | // default charm production at 5. 5 TeV |
218 | // semimuonic decay |
219 | // structure function GRVHO |
220 | // |
221 | fName = "Pythia"; |
222 | fTitle= "Particle Generator using PYTHIA"; |
223 | SetForceDecay(); |
224 | // Set random number generator |
225 | if (!AliPythiaRndm::GetPythiaRandom()) |
226 | AliPythiaRndm::SetPythiaRandom(GetRandom()); |
227 | fParticles = new TClonesArray("TParticle",1000); |
228 | SetNuclei(0,0); |
229 | } |
230 | |
231 | AliGenPythiaPlus::~AliGenPythiaPlus() |
232 | { |
233 | // Destructor |
234 | if(fEventsTime) delete fEventsTime; |
235 | } |
236 | |
237 | void AliGenPythiaPlus::SetInteractionRate(Float_t rate,Float_t timewindow) |
238 | { |
239 | // Generate pileup using user specified rate |
240 | fInteractionRate = rate; |
241 | fTimeWindow = timewindow; |
242 | GeneratePileup(); |
243 | } |
244 | |
245 | void AliGenPythiaPlus::GeneratePileup() |
246 | { |
247 | // Generate sub events time for pileup |
248 | fEventsTime = 0; |
249 | if(fInteractionRate == 0.) { |
250 | Warning("GeneratePileup","Zero interaction specified. Skipping pileup generation.\n"); |
251 | return; |
252 | } |
253 | |
254 | Int_t npart = NumberParticles(); |
255 | if(npart < 0) { |
256 | Warning("GeneratePileup","Negative number of particles. Skipping pileup generation.\n"); |
257 | return; |
258 | } |
259 | |
260 | if(fEventsTime) delete fEventsTime; |
261 | fEventsTime = new TArrayF(npart); |
262 | TArrayF &array = *fEventsTime; |
263 | for(Int_t ipart = 0; ipart < npart; ipart++) |
264 | array[ipart] = 0.; |
265 | |
266 | Float_t eventtime = 0.; |
267 | while(1) |
268 | { |
269 | eventtime += (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate); |
270 | if(eventtime > fTimeWindow) break; |
271 | array.Set(array.GetSize()+1); |
272 | array[array.GetSize()-1] = eventtime; |
273 | } |
274 | |
275 | eventtime = 0.; |
276 | while(1) |
277 | { |
278 | eventtime -= (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate); |
279 | if(TMath::Abs(eventtime) > fTimeWindow) break; |
280 | array.Set(array.GetSize()+1); |
281 | array[array.GetSize()-1] = eventtime; |
282 | } |
283 | |
284 | SetNumberParticles(fEventsTime->GetSize()); |
285 | } |
286 | |
287 | void AliGenPythiaPlus::SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi, |
288 | Float_t thresh, Float_t etseed, Float_t minet, Float_t r) |
289 | { |
290 | // Set pycell parameters |
291 | fPycellEtaMax = etamax; |
292 | fPycellNEta = neta; |
293 | fPycellNPhi = nphi; |
294 | fPycellThreshold = thresh; |
295 | fPycellEtSeed = etseed; |
296 | fPycellMinEtJet = minet; |
297 | fPycellMaxRadius = r; |
298 | } |
299 | |
300 | |
301 | |
302 | void AliGenPythiaPlus::SetEventListRange(Int_t eventFirst, Int_t eventLast) |
303 | { |
304 | // Set a range of event numbers, for which a table |
305 | // of generated particle will be printed |
306 | fDebugEventFirst = eventFirst; |
307 | fDebugEventLast = eventLast; |
308 | if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst; |
309 | } |
310 | |
311 | void AliGenPythiaPlus::Init() |
312 | { |
313 | // Initialisation |
314 | |
315 | // SetMC(AliPythia::Instance()); |
316 | // fPythia=(AliPythia*) fMCEvGen; |
317 | |
318 | // |
319 | fParentWeight=1./Float_t(fNpart); |
320 | // |
321 | |
322 | |
323 | fPythia->SetPtHardRange(fPtHardMin, fPtHardMax); |
324 | fPythia->SetYHardRange(fYHardMin, fYHardMax); |
325 | |
326 | if (fAProjectile > 0 && fATarget > 0) fPythia->SetNuclei(fAProjectile, fATarget); |
327 | // Fragmentation? |
328 | if (fFragmentation) { |
329 | fPythia->SetFragmentation(1); |
330 | } else { |
331 | fPythia->SetFragmentation(0); |
332 | } |
333 | |
334 | |
335 | // initial state radiation |
336 | fPythia->SetInitialAndFinalStateRadiation(fGinit, fGfinal); |
337 | |
338 | // pt - kick |
339 | fPythia->SetIntrinsicKt(fPtKick); |
340 | |
341 | if (fReadFromFile) { |
342 | fRL = AliRunLoader::Open(fFileName, "Partons"); |
343 | fRL->LoadKinematics(); |
344 | fRL->LoadHeader(); |
345 | } else { |
346 | fRL = 0x0; |
347 | } |
348 | // |
349 | fPythia->ProcInit(fProcess, fEnergyCMS, fStrucFunc); |
350 | // Forward Paramters to the AliPythia object |
351 | fDecayer->SetForceDecay(fForceDecay); |
352 | // Switch off Heavy Flavors on request |
353 | if (fHFoff) { |
354 | fPythia->SwitchHFOff(); |
355 | // Switch off g->QQbar splitting in decay table |
356 | ((AliDecayerPythia*) fDecayer)->HeavyFlavourOff(); |
357 | } |
358 | |
359 | fDecayer->Init(); |
360 | |
361 | |
362 | // Parent and Children Selection |
363 | switch (fProcess) |
364 | { |
365 | case kPyOldUEQ2ordered: |
366 | case kPyOldUEQ2ordered2: |
367 | case kPyOldPopcorn: |
368 | break; |
369 | case kPyCharm: |
370 | case kPyCharmUnforced: |
371 | case kPyCharmPbPbMNR: |
372 | case kPyCharmpPbMNR: |
373 | case kPyCharmppMNR: |
374 | case kPyCharmppMNRwmi: |
375 | fParentSelect[0] = 411; |
376 | fParentSelect[1] = 421; |
377 | fParentSelect[2] = 431; |
378 | fParentSelect[3] = 4122; |
379 | fFlavorSelect = 4; |
380 | break; |
381 | case kPyD0PbPbMNR: |
382 | case kPyD0pPbMNR: |
383 | case kPyD0ppMNR: |
384 | fParentSelect[0] = 421; |
385 | fFlavorSelect = 4; |
386 | break; |
387 | case kPyDPlusPbPbMNR: |
388 | case kPyDPluspPbMNR: |
389 | case kPyDPlusppMNR: |
390 | fParentSelect[0] = 411; |
391 | fFlavorSelect = 4; |
392 | break; |
393 | case kPyDPlusStrangePbPbMNR: |
394 | case kPyDPlusStrangepPbMNR: |
395 | case kPyDPlusStrangeppMNR: |
396 | fParentSelect[0] = 431; |
397 | fFlavorSelect = 4; |
398 | break; |
399 | case kPyBeauty: |
400 | case kPyBeautyPbPbMNR: |
401 | case kPyBeautypPbMNR: |
402 | case kPyBeautyppMNR: |
403 | case kPyBeautyppMNRwmi: |
404 | fParentSelect[0]= 511; |
405 | fParentSelect[1]= 521; |
406 | fParentSelect[2]= 531; |
407 | fParentSelect[3]= 5122; |
408 | fParentSelect[4]= 5132; |
409 | fParentSelect[5]= 5232; |
410 | fParentSelect[6]= 5332; |
411 | fFlavorSelect = 5; |
412 | break; |
413 | case kPyBeautyUnforced: |
414 | fParentSelect[0] = 511; |
415 | fParentSelect[1] = 521; |
416 | fParentSelect[2] = 531; |
417 | fParentSelect[3] = 5122; |
418 | fParentSelect[4] = 5132; |
419 | fParentSelect[5] = 5232; |
420 | fParentSelect[6] = 5332; |
421 | fFlavorSelect = 5; |
422 | break; |
423 | case kPyJpsiChi: |
424 | case kPyJpsi: |
425 | fParentSelect[0] = 443; |
426 | break; |
427 | case kPyMbDefault: |
428 | case kPyMb: |
429 | case kPyMbNonDiffr: |
430 | case kPyMbMSEL1: |
431 | case kPyJets: |
432 | case kPyDirectGamma: |
433 | case kPyLhwgMb: |
434 | break; |
435 | case kPyW: |
436 | case kPyZ: |
437 | break; |
438 | } |
439 | // |
440 | // |
441 | // JetFinder for Trigger |
442 | // |
443 | // Configure detector (EMCAL like) |
444 | // |
445 | fPythia->SetPycellParameters(fPycellEtaMax,fPycellNEta, fPycellNPhi, |
446 | fPycellThreshold, fPycellEtSeed, |
447 | fPycellMinEtJet, fPycellMaxRadius); |
448 | // |
449 | // This counts the total number of calls to Pyevnt() per run. |
450 | fTrialsRun = 0; |
451 | fQ = 0.; |
452 | fX1 = 0.; |
453 | fX2 = 0.; |
454 | fNev = 0 ; |
455 | // |
456 | // |
457 | // |
458 | AliGenMC::Init(); |
459 | // |
460 | // |
461 | // |
462 | if (fSetNuclei) { |
463 | fDyBoost = 0; |
464 | Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n"); |
465 | } |
466 | |
467 | if (fQuench) { |
468 | fPythia->InitQuenching(0., 0.1, 0.6e6, 0, 0.97, 30); |
469 | } |
470 | |
471 | // fPythia->SetPARJ(200, 0.0); |
472 | |
473 | // if (fQuench == 3) { |
474 | // // Nestor's change of the splittings |
475 | // fPythia->SetPARJ(200, 0.8); |
476 | // fPythia->SetMSTJ(41, 1); // QCD radiation only |
477 | // fPythia->SetMSTJ(42, 2); // angular ordering |
478 | // fPythia->SetMSTJ(44, 2); // option to run alpha_s |
479 | // fPythia->SetMSTJ(47, 0); // No correction back to hard scattering element |
480 | // fPythia->SetMSTJ(50, 0); // No coherence in first branching |
481 | // fPythia->SetPARJ(82, 1.); // Cut off for parton showers |
482 | // } |
483 | } |
484 | |
485 | void AliGenPythiaPlus::Generate() |
486 | { |
487 | // Generate one event |
488 | |
489 | fDecayer->ForceDecay(); |
490 | |
491 | Float_t polar[3] = {0,0,0}; |
492 | Float_t origin[3] = {0,0,0}; |
493 | Float_t p[4]; |
494 | // converts from mm/c to s |
495 | const Float_t kconv=0.001/2.999792458e8; |
496 | // |
497 | Int_t nt=0; |
498 | Int_t jev=0; |
499 | Int_t j, kf; |
500 | fTrials=0; |
501 | fEventTime = 0.; |
502 | |
503 | |
504 | |
505 | // Set collision vertex position |
506 | if (fVertexSmear == kPerEvent) Vertex(); |
507 | |
508 | // event loop |
509 | while(1) |
510 | { |
511 | // |
512 | // Produce event |
513 | // |
514 | // |
515 | // Switch hadronisation off |
516 | // |
517 | // fPythia->SwitchHadronisationOff(); |
518 | // |
519 | // Either produce new event or read partons from file |
520 | // |
521 | if (!fReadFromFile) { |
522 | if (!fNewMIS) { |
523 | fPythia->GenerateEvent(); |
524 | } else { |
525 | fPythia->GenerateMIEvent(); |
526 | } |
527 | fNpartons = fPythia->GetNumberOfParticles(); |
528 | } else { |
529 | printf("Loading Event %d\n",AliRunLoader::GetRunLoader()->GetEventNumber()); |
530 | fRL->GetEvent(AliRunLoader::GetRunLoader()->GetEventNumber()); |
531 | fPythia->SetNumberOfParticles(0); |
532 | fPythia->LoadEvent(fRL->Stack(), 0 , 1); |
533 | fPythia->EditEventList(21); |
534 | } |
535 | |
536 | // |
537 | // Run quenching routine |
538 | // |
539 | if (fQuench == 1) { |
540 | fPythia->Quench(); |
541 | } else if (fQuench == 2){ |
542 | fPythia->Pyquen(208., 0, 0.); |
543 | } else if (fQuench == 3) { |
544 | // Quenching is via multiplicative correction of the splittings |
545 | } |
546 | |
547 | // |
548 | // Switch hadronisation on |
549 | // |
550 | // fPythia->SwitchHadronisationOn(); |
551 | // |
552 | // .. and perform hadronisation |
553 | // printf("Calling hadronisation %d\n", fPythia->GetN()); |
554 | // fPythia->HadronizeEvent(); |
555 | fTrials++; |
556 | fPythia->GetParticles(fParticles); |
557 | Boost(); |
558 | // |
559 | // |
560 | // |
561 | Int_t i; |
562 | |
563 | fNprimaries = 0; |
564 | Int_t np = fParticles->GetEntriesFast(); |
565 | |
566 | if (np == 0) continue; |
567 | // |
568 | |
569 | // |
570 | Int_t* pParent = new Int_t[np]; |
571 | Int_t* pSelected = new Int_t[np]; |
572 | Int_t* trackIt = new Int_t[np]; |
573 | for (i = 0; i < np; i++) { |
574 | pParent[i] = -1; |
575 | pSelected[i] = 0; |
576 | trackIt[i] = 0; |
577 | } |
578 | |
579 | Int_t nc = 0; // Total n. of selected particles |
580 | Int_t nParents = 0; // Selected parents |
581 | Int_t nTkbles = 0; // Trackable particles |
582 | if (fProcess != kPyMbDefault && |
583 | fProcess != kPyMb && |
584 | fProcess != kPyJets && |
585 | fProcess != kPyDirectGamma && |
586 | fProcess != kPyMbNonDiffr && |
587 | fProcess != kPyMbMSEL1 && |
588 | fProcess != kPyW && |
589 | fProcess != kPyZ && |
590 | fProcess != kPyCharmppMNRwmi && |
591 | fProcess != kPyBeautyppMNRwmi) { |
592 | |
593 | for (i = 0; i < np; i++) { |
594 | TParticle* iparticle = (TParticle *) fParticles->At(i); |
595 | Int_t ks = iparticle->GetStatusCode(); |
596 | kf = CheckPDGCode(iparticle->GetPdgCode()); |
597 | // No initial state partons |
598 | if (ks==21) continue; |
599 | // |
600 | // Heavy Flavor Selection |
601 | // |
602 | // quark ? |
603 | kf = TMath::Abs(kf); |
604 | Int_t kfl = kf; |
605 | // Resonance |
606 | |
607 | if (kfl > 100000) kfl %= 100000; |
608 | if (kfl > 10000) kfl %= 10000; |
609 | // meson ? |
610 | if (kfl > 10) kfl/=100; |
611 | // baryon |
612 | if (kfl > 10) kfl/=10; |
613 | Int_t ipa = (fPythia->Version() == 6) ? (iparticle->GetFirstMother() - 1) :(iparticle->GetFirstMother()) ; |
614 | Int_t kfMo = 0; |
615 | // |
616 | // Establish mother daughter relation between heavy quarks and mesons |
617 | // |
618 | if (kf >= fFlavorSelect && kf <= 6) { |
619 | Int_t idau = (fPythia->Version() == 6) ? (iparticle->GetFirstDaughter() - 1) :(iparticle->GetFirstDaughter()); |
620 | if (idau > -1) { |
621 | TParticle* daughter = (TParticle *) fParticles->At(idau); |
622 | Int_t pdgD = daughter->GetPdgCode(); |
623 | if (pdgD == 91 || pdgD == 92) { |
624 | Int_t jmin = (fPythia->Version() == 6) ? (daughter->GetFirstDaughter() - 1) : (daughter->GetFirstDaughter()); |
625 | Int_t jmax = (fPythia->Version() == 6) ? (daughter->GetLastDaughter() - 1) : (daughter->GetLastDaughter()); |
626 | |
627 | for (Int_t j = jmin; j <= jmax; j++) |
628 | ((TParticle *) fParticles->At(j))->SetFirstMother(i+1); |
629 | } // is string or cluster |
630 | } // has daughter |
631 | } // heavy quark |
632 | |
633 | |
634 | if (ipa > -1) { |
635 | TParticle * mother = (TParticle *) fParticles->At(ipa); |
636 | kfMo = TMath::Abs(mother->GetPdgCode()); |
637 | } |
638 | |
639 | // What to keep in Stack? |
640 | Bool_t flavorOK = kFALSE; |
641 | Bool_t selectOK = kFALSE; |
642 | if (fFeedDownOpt) { |
643 | if (kfl >= fFlavorSelect) flavorOK = kTRUE; |
644 | } else { |
645 | if (kfl > fFlavorSelect) { |
646 | nc = -1; |
647 | break; |
648 | } |
649 | if (kfl == fFlavorSelect) flavorOK = kTRUE; |
650 | } |
651 | switch (fStackFillOpt) { |
652 | case kFlavorSelection: |
653 | selectOK = kTRUE; |
654 | break; |
655 | case kParentSelection: |
656 | if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE; |
657 | break; |
658 | } |
659 | if (flavorOK && selectOK) { |
660 | // |
661 | // Heavy flavor hadron or quark |
662 | // |
663 | // Kinematic seletion on final state heavy flavor mesons |
664 | if (ParentSelected(kf) && !KinematicSelection(iparticle, 0)) |
665 | { |
666 | continue; |
667 | } |
668 | pSelected[i] = 1; |
669 | if (ParentSelected(kf)) ++nParents; // Update parent count |
670 | // printf("\n particle (HF) %d %d %d", i, pSelected[i], kf); |
671 | } else { |
672 | // Kinematic seletion on decay products |
673 | if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf) |
674 | && !KinematicSelection(iparticle, 1)) |
675 | { |
676 | continue; |
677 | } |
678 | // |
679 | // Decay products |
680 | // Select if mother was selected and is not tracked |
681 | |
682 | if (pSelected[ipa] && |
683 | !trackIt[ipa] && // mother will be tracked ? |
684 | kfMo != 5 && // mother is b-quark, don't store fragments |
685 | kfMo != 4 && // mother is c-quark, don't store fragments |
686 | kf != 92) // don't store string |
687 | { |
688 | // |
689 | // Semi-stable or de-selected: diselect decay products: |
690 | // |
691 | // |
692 | if (pSelected[i] == -1 || fDecayer->GetLifetime(kf) > fMaxLifeTime) |
693 | { |
694 | Int_t ipF = iparticle->GetFirstDaughter(); |
695 | Int_t ipL = iparticle->GetLastDaughter(); |
696 | if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1; |
697 | } |
698 | // printf("\n particle (decay) %d %d %d", i, pSelected[i], kf); |
699 | pSelected[i] = (pSelected[i] == -1) ? 0 : 1; |
700 | } |
701 | } |
702 | if (pSelected[i] == -1) pSelected[i] = 0; |
703 | if (!pSelected[i]) continue; |
704 | // Count quarks only if you did not include fragmentation |
705 | if (fFragmentation && kf <= 10) continue; |
706 | |
707 | nc++; |
708 | // Decision on tracking |
709 | trackIt[i] = 0; |
710 | // |
711 | // Track final state particle |
712 | if (ks == 1) trackIt[i] = 1; |
713 | // Track semi-stable particles |
714 | if ((ks == 1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime)) trackIt[i] = 1; |
715 | // Track particles selected by process if undecayed. |
716 | if (fForceDecay == kNoDecay) { |
717 | if (ParentSelected(kf)) trackIt[i] = 1; |
718 | } else { |
719 | if (ParentSelected(kf)) trackIt[i] = 0; |
720 | } |
721 | if (trackIt[i] == 1) ++nTkbles; // Update trackable counter |
722 | // |
723 | // |
724 | |
725 | } // particle selection loop |
726 | if (nc > 0) { |
727 | for (i = 0; i < np; i++) { |
728 | if (!pSelected[i]) continue; |
729 | TParticle * iparticle = (TParticle *) fParticles->At(i); |
730 | kf = CheckPDGCode(iparticle->GetPdgCode()); |
731 | Int_t ks = iparticle->GetStatusCode(); |
732 | p[0] = iparticle->Px(); |
733 | p[1] = iparticle->Py(); |
734 | p[2] = iparticle->Pz(); |
735 | p[3] = iparticle->Energy(); |
736 | |
737 | origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm] |
738 | origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm] |
739 | origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm] |
740 | |
741 | Float_t tof = kconv*iparticle->T(); |
742 | Int_t ipa = (fPythia->Version() == 6) ? (iparticle->GetFirstMother() - 1) :(iparticle->GetFirstMother()) ; |
743 | Int_t iparent = (ipa > -1) ? pParent[ipa] : -1; |
744 | |
745 | PushTrack(fTrackIt*trackIt[i], iparent, kf, |
746 | p[0], p[1], p[2], p[3], |
747 | origin[0], origin[1], origin[2], tof, |
748 | polar[0], polar[1], polar[2], |
749 | kPPrimary, nt, 1., ks); |
750 | pParent[i] = nt; |
751 | KeepTrack(nt); |
752 | fNprimaries++; |
753 | } // PushTrack loop |
754 | } |
755 | } else { |
756 | nc = GenerateMB(); |
757 | } // mb ? |
758 | |
759 | GetSubEventTime(); |
760 | |
761 | delete[] pParent; |
762 | delete[] pSelected; |
763 | delete[] trackIt; |
764 | |
765 | if (nc > 0) { |
766 | switch (fCountMode) { |
767 | case kCountAll: |
768 | // printf(" Count all \n"); |
769 | jev += nc; |
770 | break; |
771 | case kCountParents: |
772 | // printf(" Count parents \n"); |
773 | jev += nParents; |
774 | break; |
775 | case kCountTrackables: |
776 | // printf(" Count trackable \n"); |
777 | jev += nTkbles; |
778 | break; |
779 | } |
780 | if (jev >= fNpart || fNpart == -1) { |
781 | fKineBias=Float_t(fNpart)/Float_t(fTrials); |
782 | fPythia->GetXandQ(fQ, fX1, fX2); |
783 | fTrialsRun += fTrials; |
784 | fNev++; |
785 | MakeHeader(); |
786 | break; |
787 | } |
788 | } |
789 | } // event loop |
790 | SetHighWaterMark(nt); |
791 | // Adjust weight due to kinematic selection |
792 | // AdjustWeights(); |
793 | // get cross-section |
794 | fXsection = fPythia->GetXSection(); |
795 | } |
796 | |
797 | Int_t AliGenPythiaPlus::GenerateMB() |
798 | { |
799 | // |
800 | // Min Bias selection and other global selections |
801 | // |
802 | Int_t i, kf, nt, iparent; |
803 | Int_t nc = 0; |
804 | Float_t p[4]; |
805 | Float_t polar[3] = {0,0,0}; |
806 | Float_t origin[3] = {0,0,0}; |
807 | // converts from mm/c to s |
808 | const Float_t kconv = 0.001 / 2.999792458e8; |
809 | |
810 | Int_t np = (fHadronisation) ? fParticles->GetEntriesFast() : fNpartons; |
811 | |
812 | Int_t* pParent = new Int_t[np]; |
813 | for (i=0; i< np; i++) pParent[i] = -1; |
814 | if (fProcess == kPyJets || fProcess == kPyDirectGamma) { |
815 | TParticle* jet1 = (TParticle *) fParticles->At(6); |
816 | TParticle* jet2 = (TParticle *) fParticles->At(7); |
817 | if (!CheckTrigger(jet1, jet2)) { |
818 | delete [] pParent; |
819 | return 0; |
820 | } |
821 | } |
822 | |
823 | // Select jets with fragmentation photon or pi0 going to PHOS or EMCAL |
824 | if (fProcess == kPyJets && (fFragPhotonInCalo || fPi0InCalo) ) { |
825 | |
826 | Bool_t ok = kFALSE; |
827 | |
828 | Int_t pdg = 0; |
829 | if (fFragPhotonInCalo) pdg = 22 ; // Photon |
830 | else if (fPi0InCalo) pdg = 111 ; // Pi0 |
831 | |
832 | for (i=0; i< np; i++) { |
833 | TParticle* iparticle = (TParticle *) fParticles->At(i); |
834 | if(iparticle->GetStatusCode()==1 && iparticle->GetPdgCode()==pdg && |
835 | iparticle->Pt() > fFragPhotonOrPi0MinPt){ |
836 | Int_t imother = (fPythia->Version() == 6) ? (iparticle->GetFirstMother() - 1) :(iparticle->GetFirstMother()) ; |
837 | TParticle* pmother = (TParticle *) fParticles->At(imother); |
838 | if(pdg == 111 || |
839 | (pdg == 22 && pmother->GetStatusCode() != 11))//No photon from hadron decay |
840 | { |
841 | Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees |
842 | Float_t eta =TMath::Abs(iparticle->Eta());//in calos etamin=-etamax |
843 | if((fCheckEMCAL && IsInEMCAL(phi,eta)) || |
844 | (fCheckPHOS && IsInPHOS(phi,eta)) ) |
845 | ok =kTRUE; |
846 | } |
847 | } |
848 | } |
849 | if(!ok) |
850 | return 0; |
851 | } |
852 | |
853 | |
854 | // Select events with a photon pt > min pt going to PHOS eta acceptance or exactly PHOS eta phi |
855 | if ((fProcess == kPyJets || fProcess == kPyDirectGamma) && fPhotonInCalo && (fCheckPHOSeta || fCheckPHOS)){ |
856 | |
857 | Bool_t okd = kFALSE; |
858 | |
859 | Int_t pdg = 22; |
860 | Int_t iphcand = -1; |
861 | for (i=0; i< np; i++) { |
862 | TParticle* iparticle = (TParticle *) fParticles->At(i); |
863 | Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees |
864 | Float_t eta =TMath::Abs(iparticle->Eta());//in calos etamin=-etamax |
865 | |
866 | if(iparticle->GetStatusCode() == 1 |
867 | && iparticle->GetPdgCode() == pdg |
868 | && iparticle->Pt() > fPhotonMinPt |
869 | && eta < fPHOSEta){ |
870 | |
871 | // first check if the photon is in PHOS phi |
872 | if(IsInPHOS(phi,eta)){ |
873 | okd = kTRUE; |
874 | break; |
875 | } |
876 | if(fCheckPHOSeta) iphcand = i; // candiate photon to rotate in phi |
877 | |
878 | } |
879 | } |
880 | |
881 | if(!okd && iphcand != -1) // execute rotation in phi |
882 | RotatePhi(iphcand,okd); |
883 | |
884 | if(!okd) |
885 | return 0; |
886 | } |
887 | |
888 | if (fTriggerParticle) { |
889 | Bool_t triggered = kFALSE; |
890 | for (i = 0; i < np; i++) { |
891 | TParticle * iparticle = (TParticle *) fParticles->At(i); |
892 | kf = CheckPDGCode(iparticle->GetPdgCode()); |
893 | if (kf != fTriggerParticle) continue; |
894 | if (iparticle->Pt() == 0.) continue; |
895 | if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue; |
896 | triggered = kTRUE; |
897 | break; |
898 | } |
899 | if (!triggered) { |
900 | delete [] pParent; |
901 | return 0; |
902 | } |
903 | } |
904 | |
905 | |
906 | // Check if there is a ccbar or bbbar pair with at least one of the two |
907 | // in fYMin < y < fYMax |
908 | if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi) { |
909 | TParticle *hvq; |
910 | Bool_t theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE; |
911 | Float_t yQ; |
912 | Int_t pdgQ; |
913 | for(i=0; i<np; i++) { |
914 | hvq = (TParticle*)fParticles->At(i); |
915 | pdgQ = hvq->GetPdgCode(); |
916 | if(TMath::Abs(pdgQ) != fFlavorSelect) continue; |
917 | if(pdgQ>0) { theQ=kTRUE; } else { theQbar=kTRUE; } |
918 | yQ = 0.5*TMath::Log((hvq->Energy()+hvq->Pz()+1.e-13)/ |
919 | (hvq->Energy()-hvq->Pz()+1.e-13)); |
920 | if(yQ>fYMin && yQ<fYMax) inYcut=kTRUE; |
921 | } |
922 | if (!theQ || !theQbar || !inYcut) { |
923 | delete[] pParent; |
924 | return 0; |
925 | } |
926 | } |
927 | |
928 | //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff |
929 | if ( (fProcess == kPyW || |
930 | fProcess == kPyZ || |
931 | fProcess == kPyMbDefault || |
932 | fProcess == kPyMb || |
933 | fProcess == kPyMbNonDiffr) |
934 | && (fCutOnChild == 1) ) { |
935 | if ( !CheckKinematicsOnChild() ) { |
936 | delete[] pParent; |
937 | return 0; |
938 | } |
939 | } |
940 | |
941 | |
942 | for (i = 0; i < np; i++) { |
943 | Int_t trackIt = 0; |
944 | TParticle * iparticle = (TParticle *) fParticles->At(i); |
945 | kf = CheckPDGCode(iparticle->GetPdgCode()); |
946 | Int_t ks = iparticle->GetStatusCode(); |
947 | Int_t km = iparticle->GetFirstMother(); |
948 | if ((ks == 1 && kf!=0 && KinematicSelection(iparticle, 0)) || |
949 | (ks != 1) || |
950 | (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) { |
951 | nc++; |
952 | if (ks == 1) trackIt = 1; |
953 | |
954 | Int_t ipa = (fPythia->Version() == 6) ? (iparticle->GetFirstMother() - 1) :(iparticle->GetFirstMother()) ; |
955 | iparent = (ipa > -1) ? pParent[ipa] : -1; |
956 | if (ipa >= np) fPythia->EventListing(); |
957 | |
958 | // |
959 | // store track information |
960 | p[0] = iparticle->Px(); |
961 | p[1] = iparticle->Py(); |
962 | p[2] = iparticle->Pz(); |
963 | p[3] = iparticle->Energy(); |
964 | |
965 | |
966 | origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm] |
967 | origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm] |
968 | origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm] |
969 | |
970 | Float_t tof = fEventTime + kconv * iparticle->T(); |
971 | |
972 | PushTrack(fTrackIt*trackIt, iparent, kf, |
973 | p[0], p[1], p[2], p[3], |
974 | origin[0], origin[1], origin[2], tof, |
975 | polar[0], polar[1], polar[2], |
976 | kPPrimary, nt, 1., ks); |
977 | fNprimaries++; |
978 | |
979 | |
980 | // |
981 | // Special Treatment to store color-flow |
982 | // |
983 | if (ks == 3 || ks == 13 || ks == 14) { |
984 | TParticle* particle = 0; |
985 | if (fStack) { |
986 | particle = fStack->Particle(nt); |
987 | } else { |
988 | particle = gAlice->Stack()->Particle(nt); |
989 | } |
990 | // particle->SetFirstDaughter(fPythia->GetK(2, i)); |
991 | // particle->SetLastDaughter(fPythia->GetK(3, i)); |
992 | } |
993 | |
994 | KeepTrack(nt); |
995 | pParent[i] = nt; |
996 | SetHighWaterMark(nt); |
997 | |
998 | } // select particle |
999 | } // particle loop |
1000 | |
1001 | delete[] pParent; |
1002 | |
1003 | return 1; |
1004 | } |
1005 | |
1006 | |
1007 | void AliGenPythiaPlus::FinishRun() |
1008 | { |
1009 | // Print x-section summary |
1010 | fPythia->PrintStatistics(); |
1011 | |
1012 | if (fNev > 0.) { |
1013 | fQ /= fNev; |
1014 | fX1 /= fNev; |
1015 | fX2 /= fNev; |
1016 | } |
1017 | |
1018 | printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun); |
1019 | printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2); |
1020 | } |
1021 | |
1022 | void AliGenPythiaPlus::AdjustWeights() const |
1023 | { |
1024 | // Adjust the weights after generation of all events |
1025 | // |
1026 | if (gAlice) { |
1027 | TParticle *part; |
1028 | Int_t ntrack=gAlice->GetMCApp()->GetNtrack(); |
1029 | for (Int_t i=0; i<ntrack; i++) { |
1030 | part= gAlice->GetMCApp()->Particle(i); |
1031 | part->SetWeight(part->GetWeight()*fKineBias); |
1032 | } |
1033 | } |
1034 | } |
1035 | |
1036 | void AliGenPythiaPlus::SetNuclei(Int_t a1, Int_t a2) |
1037 | { |
1038 | // Treat protons as inside nuclei with mass numbers a1 and a2 |
1039 | |
1040 | fAProjectile = a1; |
1041 | fATarget = a2; |
1042 | fSetNuclei = kTRUE; |
1043 | } |
1044 | |
1045 | |
1046 | void AliGenPythiaPlus::MakeHeader() |
1047 | { |
1048 | // |
1049 | // Make header for the simulated event |
1050 | // |
1051 | if (gAlice) { |
1052 | if (gAlice->GetEvNumber()>=fDebugEventFirst && |
1053 | gAlice->GetEvNumber()<=fDebugEventLast) fPythia->EventListing(); |
1054 | } |
1055 | |
1056 | // Builds the event header, to be called after each event |
1057 | if (fHeader) delete fHeader; |
1058 | fHeader = new AliGenPythiaEventHeader("Pythia"); |
1059 | // |
1060 | // Event type |
1061 | ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->ProcessCode()); |
1062 | // |
1063 | // Number of trials |
1064 | ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials); |
1065 | // |
1066 | // Event Vertex |
1067 | fHeader->SetPrimaryVertex(fVertex); |
1068 | |
1069 | // |
1070 | // Number of primaries |
1071 | fHeader->SetNProduced(fNprimaries); |
1072 | // |
1073 | // Jets that have triggered |
1074 | |
1075 | if (fProcess == kPyJets) |
1076 | { |
1077 | Int_t ntrig, njet; |
1078 | Float_t jets[4][10]; |
1079 | GetJets(njet, ntrig, jets); |
1080 | |
1081 | |
1082 | for (Int_t i = 0; i < ntrig; i++) { |
1083 | ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i], |
1084 | jets[3][i]); |
1085 | } |
1086 | } |
1087 | // |
1088 | // Copy relevant information from external header, if present. |
1089 | // |
1090 | Float_t uqJet[4]; |
1091 | |
1092 | if (fRL) { |
1093 | AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader()); |
1094 | for (Int_t i = 0; i < exHeader->NTriggerJets(); i++) |
1095 | { |
1096 | printf("Adding Jet %d %d \n", i, exHeader->NTriggerJets()); |
1097 | |
1098 | |
1099 | exHeader->TriggerJet(i, uqJet); |
1100 | ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]); |
1101 | } |
1102 | } |
1103 | // |
1104 | // Store quenching parameters |
1105 | // |
1106 | if (fQuench){ |
1107 | Double_t z[4]; |
1108 | Double_t xp, yp; |
1109 | if (fQuench == 1) { |
1110 | // Pythia::Quench() |
1111 | fPythia->GetQuenchingParameters(xp, yp, z); |
1112 | } else { |
1113 | // Pyquen |
1114 | Double_t r1 = PARIMP.rb1; |
1115 | Double_t r2 = PARIMP.rb2; |
1116 | Double_t b = PARIMP.b1; |
1117 | Double_t r = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b); |
1118 | Double_t phi = PARIMP.psib1; |
1119 | xp = r * TMath::Cos(phi); |
1120 | yp = r * TMath::Sin(phi); |
1121 | |
1122 | } |
1123 | ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp); |
1124 | ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z); |
1125 | } |
1126 | // |
1127 | // Store pt^hard |
1128 | ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetPtHard()); |
1129 | // |
1130 | // Pass header |
1131 | // |
1132 | AddHeader(fHeader); |
1133 | fHeader = 0x0; |
1134 | } |
1135 | |
1136 | Bool_t AliGenPythiaPlus::CheckTrigger(TParticle* jet1, TParticle* jet2) |
1137 | { |
1138 | // Check the kinematic trigger condition |
1139 | // |
1140 | Double_t eta[2]; |
1141 | eta[0] = jet1->Eta(); |
1142 | eta[1] = jet2->Eta(); |
1143 | Double_t phi[2]; |
1144 | phi[0] = jet1->Phi(); |
1145 | phi[1] = jet2->Phi(); |
1146 | Int_t pdg[2]; |
1147 | pdg[0] = jet1->GetPdgCode(); |
1148 | pdg[1] = jet2->GetPdgCode(); |
1149 | Bool_t triggered = kFALSE; |
1150 | |
1151 | if (fProcess == kPyJets) { |
1152 | Int_t njets = 0; |
1153 | Int_t ntrig = 0; |
1154 | Float_t jets[4][10]; |
1155 | // |
1156 | // Use Pythia clustering on parton level to determine jet axis |
1157 | // |
1158 | GetJets(njets, ntrig, jets); |
1159 | |
1160 | if (ntrig || fEtMinJet == 0.) triggered = kTRUE; |
1161 | // |
1162 | } else { |
1163 | Int_t ij = 0; |
1164 | Int_t ig = 1; |
1165 | if (pdg[0] == kGamma) { |
1166 | ij = 1; |
1167 | ig = 0; |
1168 | } |
1169 | //Check eta range first... |
1170 | if ((eta[ij] < fEtaMaxJet && eta[ij] > fEtaMinJet) && |
1171 | (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma)) |
1172 | { |
1173 | //Eta is okay, now check phi range |
1174 | if ((phi[ij] < fPhiMaxJet && phi[ij] > fPhiMinJet) && |
1175 | (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma)) |
1176 | { |
1177 | triggered = kTRUE; |
1178 | } |
1179 | } |
1180 | } |
1181 | return triggered; |
1182 | } |
1183 | |
1184 | |
1185 | |
1186 | Bool_t AliGenPythiaPlus::CheckKinematicsOnChild(){ |
1187 | // |
1188 | //Checking Kinematics on Child (status code 1, particle code ?, kin cuts |
1189 | // |
1190 | Bool_t checking = kFALSE; |
1191 | Int_t j, kcode, ks, km; |
1192 | Int_t nPartAcc = 0; //number of particles in the acceptance range |
1193 | Int_t numberOfAcceptedParticles = 1; |
1194 | if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; } |
1195 | Int_t npart = fParticles->GetEntriesFast(); |
1196 | |
1197 | for (j = 0; j<npart; j++) { |
1198 | TParticle * jparticle = (TParticle *) fParticles->At(j); |
1199 | kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) ); |
1200 | ks = jparticle->GetStatusCode(); |
1201 | km = jparticle->GetFirstMother(); |
1202 | |
1203 | if( (ks == 1) && (kcode == fPdgCodeParticleforAcceptanceCut) && (KinematicSelection(jparticle,1)) ){ |
1204 | nPartAcc++; |
1205 | } |
1206 | if( numberOfAcceptedParticles <= nPartAcc){ |
1207 | checking = kTRUE; |
1208 | break; |
1209 | } |
1210 | } |
1211 | |
1212 | return checking; |
1213 | } |
1214 | |
1215 | void AliGenPythiaPlus::RecJetsUA1(Int_t& njets, Float_t jets [4][50]) |
1216 | { |
1217 | // |
1218 | // Calls the Pythia jet finding algorithm to find jets in the current event |
1219 | // |
1220 | // |
1221 | // |
1222 | // Save jets |
1223 | // |
1224 | // Run Jet Finder |
1225 | fPythia->Pycell(njets); |
1226 | Int_t i; |
1227 | for (i = 0; i < njets; i++) { |
1228 | Float_t px, py, pz, e; |
1229 | fPythia->GetJet(i, px, py, pz, e); |
1230 | jets[0][i] = px; |
1231 | jets[1][i] = py; |
1232 | jets[2][i] = pz; |
1233 | jets[3][i] = e; |
1234 | } |
1235 | } |
1236 | |
1237 | |
1238 | void AliGenPythiaPlus::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10]) |
1239 | { |
1240 | // |
1241 | // Calls the Pythia clustering algorithm to find jets in the current event |
1242 | // |
1243 | nJets = 0; |
1244 | nJetsTrig = 0; |
1245 | if (fJetReconstruction == kCluster) { |
1246 | // |
1247 | // Configure cluster algorithm |
1248 | // |
1249 | // fPythia->SetPARU(43, 2.); |
1250 | // fPythia->SetMSTU(41, 1); |
1251 | // |
1252 | // Call cluster algorithm |
1253 | // |
1254 | fPythia->Pyclus(nJets); |
1255 | // |
1256 | // Loading jets from common block |
1257 | // |
1258 | } else { |
1259 | |
1260 | // |
1261 | // Run Jet Finder |
1262 | fPythia->Pycell(nJets); |
1263 | } |
1264 | |
1265 | Int_t i; |
1266 | for (i = 0; i < nJets; i++) { |
1267 | Float_t px, py, pz, e; |
1268 | fPythia->GetJet(i, px, py, pz, e); |
1269 | Float_t pt = TMath::Sqrt(px * px + py * py); |
1270 | Float_t phi = TMath::Pi() + TMath::ATan2(-py, -px); |
1271 | Float_t theta = TMath::ATan2(pt,pz); |
1272 | Float_t et = e * TMath::Sin(theta); |
1273 | Float_t eta = -TMath::Log(TMath::Tan(theta / 2.)); |
1274 | if ( |
1275 | eta > fEtaMinJet && eta < fEtaMaxJet && |
1276 | phi > fPhiMinJet && phi < fPhiMaxJet && |
1277 | et > fEtMinJet && et < fEtMaxJet |
1278 | ) |
1279 | { |
1280 | jets[0][nJetsTrig] = px; |
1281 | jets[1][nJetsTrig] = py; |
1282 | jets[2][nJetsTrig] = pz; |
1283 | jets[3][nJetsTrig] = e; |
1284 | nJetsTrig++; |
1285 | } else { |
1286 | } |
1287 | } |
1288 | } |
1289 | |
1290 | void AliGenPythiaPlus::GetSubEventTime() |
1291 | { |
1292 | // Calculates time of the next subevent |
1293 | fEventTime = 0.; |
1294 | if (fEventsTime) { |
1295 | TArrayF &array = *fEventsTime; |
1296 | fEventTime = array[fCurSubEvent++]; |
1297 | } |
1298 | // printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime); |
1299 | return; |
1300 | } |
1301 | |
1302 | |
1303 | |
1304 | |
1305 | Bool_t AliGenPythiaPlus::IsInEMCAL(Float_t phi, Float_t eta) |
1306 | { |
1307 | // Is particle in EMCAL acceptance? |
1308 | // phi in degrees, etamin=-etamax |
1309 | if(phi > fEMCALMinPhi && phi < fEMCALMaxPhi && |
1310 | eta < fEMCALEta ) |
1311 | return kTRUE; |
1312 | else |
1313 | return kFALSE; |
1314 | } |
1315 | |
1316 | Bool_t AliGenPythiaPlus::IsInPHOS(Float_t phi, Float_t eta) |
1317 | { |
1318 | // Is particle in PHOS acceptance? |
1319 | // Acceptance slightly larger considered. |
1320 | // phi in degrees, etamin=-etamax |
1321 | if(phi > fPHOSMinPhi && phi < fPHOSMaxPhi && |
1322 | eta < fPHOSEta ) |
1323 | return kTRUE; |
1324 | else |
1325 | return kFALSE; |
1326 | } |
1327 | |
1328 | void AliGenPythiaPlus::RotatePhi(Int_t iphcand, Bool_t& okdd) |
1329 | { |
1330 | //calculate the new position random between fPHOSMinPhi and fPHOSMaxPhi |
1331 | Double_t phiPHOSmin = TMath::Pi()*fPHOSMinPhi/180; |
1332 | Double_t phiPHOSmax = TMath::Pi()*fPHOSMaxPhi/180; |
1333 | Double_t phiPHOS = gRandom->Uniform(phiPHOSmin,phiPHOSmax); |
1334 | |
1335 | //calculate deltaphi |
1336 | TParticle* ph = (TParticle *) fParticles->At(iphcand); |
1337 | Double_t phphi = ph->Phi(); |
1338 | Double_t deltaphi = phiPHOS - phphi; |
1339 | |
1340 | |
1341 | |
1342 | //loop for all particles and produce the phi rotation |
1343 | Int_t np = (fHadronisation) ? fParticles->GetEntriesFast() : fNpartons; |
1344 | Double_t oldphi, newphi; |
1345 | Double_t newVx, newVy, R, Vz, time; |
1346 | Double_t newPx, newPy, pt, Pz, e; |
1347 | for(Int_t i=0; i< np; i++) { |
1348 | TParticle* iparticle = (TParticle *) fParticles->At(i); |
1349 | oldphi = iparticle->Phi(); |
1350 | newphi = oldphi + deltaphi; |
1351 | if(newphi < 0) newphi = 2*TMath::Pi() + newphi; // correct angle |
1352 | if(newphi > 2*TMath::Pi()) newphi = newphi - 2*TMath::Pi(); // correct angle |
1353 | |
1354 | R = iparticle->R(); |
1355 | newVx = R*TMath::Cos(newphi); |
1356 | newVy = R*TMath::Sin(newphi); |
1357 | Vz = iparticle->Vz(); // don't transform |
1358 | time = iparticle->T(); // don't transform |
1359 | |
1360 | pt = iparticle->Pt(); |
1361 | newPx = pt*TMath::Cos(newphi); |
1362 | newPy = pt*TMath::Sin(newphi); |
1363 | Pz = iparticle->Pz(); // don't transform |
1364 | e = iparticle->Energy(); // don't transform |
1365 | |
1366 | // apply rotation |
1367 | iparticle->SetProductionVertex(newVx, newVy, Vz, time); |
1368 | iparticle->SetMomentum(newPx, newPy, Pz, e); |
1369 | |
1370 | } //end particle loop |
1371 | |
1372 | // now let's check that we put correctly the candidate photon in PHOS |
1373 | Float_t phi = ph->Phi()*180./TMath::Pi(); //Convert to degrees |
1374 | Float_t eta =TMath::Abs(ph->Eta());//in calos etamin=-etamax |
1375 | if(IsInPHOS(phi,eta)) |
1376 | okdd = kTRUE; |
1377 | } |
1378 | |
1379 | |
1380 | #ifdef never |
1381 | void AliGenPythiaPlus::Streamer(TBuffer &R__b) |
1382 | { |
1383 | // Stream an object of class AliGenPythia. |
1384 | |
1385 | if (R__b.IsReading()) { |
1386 | Version_t R__v = R__b.ReadVersion(); if (R__v) { } |
1387 | AliGenerator::Streamer(R__b); |
1388 | R__b >> (Int_t&)fProcess; |
1389 | R__b >> (Int_t&)fStrucFunc; |
1390 | R__b >> (Int_t&)fForceDecay; |
1391 | R__b >> fEnergyCMS; |
1392 | R__b >> fKineBias; |
1393 | R__b >> fTrials; |
1394 | fParentSelect.Streamer(R__b); |
1395 | fChildSelect.Streamer(R__b); |
1396 | R__b >> fXsection; |
1397 | // (AliPythia::Instance())->Streamer(R__b); |
1398 | R__b >> fPtHardMin; |
1399 | R__b >> fPtHardMax; |
1400 | // if (fDecayer) fDecayer->Streamer(R__b); |
1401 | } else { |
1402 | R__b.WriteVersion(AliGenPythiaPlus::IsA()); |
1403 | AliGenerator::Streamer(R__b); |
1404 | R__b << (Int_t)fProcess; |
1405 | R__b << (Int_t)fStrucFunc; |
1406 | R__b << (Int_t)fForceDecay; |
1407 | R__b << fEnergyCMS; |
1408 | R__b << fKineBias; |
1409 | R__b << fTrials; |
1410 | fParentSelect.Streamer(R__b); |
1411 | fChildSelect.Streamer(R__b); |
1412 | R__b << fXsection; |
1413 | // R__b << fPythia; |
1414 | R__b << fPtHardMin; |
1415 | R__b << fPtHardMax; |
1416 | // fDecayer->Streamer(R__b); |
1417 | } |
1418 | } |
1419 | #endif |
1420 | |
1421 | |
1422 | |