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