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