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