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