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