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