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