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