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