]> 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
Coverity fix.
[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
500 Float_t polar[3] = {0,0,0};
501 Float_t origin[3] = {0,0,0};
502 Float_t p[4];
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
751 Float_t tof = kconv*iparticle->T();
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 }
859 if(!ok)
860 return 0;
861 }
862
863
864 // Select events with a photon pt > min pt going to PHOS eta acceptance or exactly PHOS eta phi
865 if ((fProcess == kPyJets || fProcess == kPyDirectGamma) && fPhotonInCalo && (fCheckPHOSeta || fCheckPHOS)){
866
867 Bool_t okd = kFALSE;
868
869 Int_t pdg = 22;
870 Int_t iphcand = -1;
871 for (i=0; i< np; i++) {
8507138f 872 TParticle* iparticle = (TParticle *) fParticles.At(i);
39d810c8 873 Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
874 Float_t eta =TMath::Abs(iparticle->Eta());//in calos etamin=-etamax
875
876 if(iparticle->GetStatusCode() == 1
877 && iparticle->GetPdgCode() == pdg
878 && iparticle->Pt() > fPhotonMinPt
879 && eta < fPHOSEta){
880
881 // first check if the photon is in PHOS phi
882 if(IsInPHOS(phi,eta)){
883 okd = kTRUE;
884 break;
885 }
886 if(fCheckPHOSeta) iphcand = i; // candiate photon to rotate in phi
887
888 }
889 }
890
891 if(!okd && iphcand != -1) // execute rotation in phi
892 RotatePhi(iphcand,okd);
893
894 if(!okd)
895 return 0;
896 }
897
898 if (fTriggerParticle) {
899 Bool_t triggered = kFALSE;
900 for (i = 0; i < np; i++) {
8507138f 901 TParticle * iparticle = (TParticle *) fParticles.At(i);
39d810c8 902 kf = CheckPDGCode(iparticle->GetPdgCode());
903 if (kf != fTriggerParticle) continue;
904 if (iparticle->Pt() == 0.) continue;
905 if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
906 triggered = kTRUE;
907 break;
908 }
909 if (!triggered) {
910 delete [] pParent;
911 return 0;
912 }
913 }
914
915
916 // Check if there is a ccbar or bbbar pair with at least one of the two
917 // in fYMin < y < fYMax
918 if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi) {
9ccc3d0e 919 TParticle *partCheck;
920 TParticle *mother;
39d810c8 921 Bool_t theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
9ccc3d0e 922 Bool_t theChild=kFALSE;
923 Float_t y;
924 Int_t pdg,mpdg,mpdgUpperFamily;
39d810c8 925 for(i=0; i<np; i++) {
9ccc3d0e 926 partCheck = (TParticle*)fParticles.At(i);
927 pdg = partCheck->GetPdgCode();
928 if(TMath::Abs(pdg) == fFlavorSelect) { // quark
929 if(pdg>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
930 y = 0.5*TMath::Log((partCheck->Energy()+partCheck->Pz()+1.e-13)/
931 (partCheck->Energy()-partCheck->Pz()+1.e-13));
932 if(y>fYMin && y<fYMax) inYcut=kTRUE;
933 }
934
935 if(fCutOnChild && TMath::Abs(pdg) == fPdgCodeParticleforAcceptanceCut) {
936 Int_t mi = partCheck->GetFirstMother() - 1;
937 if(mi<0) continue;
938 mother = (TParticle*)fParticles.At(mi);
939 mpdg=TMath::Abs(mother->GetPdgCode());
75fde694 940 mpdgUpperFamily=(mpdg>1000 ? mpdg+1000 : mpdg+100); // keep e from c from b
9ccc3d0e 941 if ( ParentSelected(mpdg) ||
942 (fFlavorSelect==5 && ParentSelected(mpdgUpperFamily))) {
943 if (KinematicSelection(partCheck,1)) {
944 theChild=kTRUE;
945 }
946 }
947 }
39d810c8 948 }
9ccc3d0e 949 if (!theQ || !theQbar || !inYcut) { // one of the c/b conditions not satisfied
39d810c8 950 delete[] pParent;
951 return 0;
952 }
9ccc3d0e 953 if (fCutOnChild && !theChild) { // one of the child conditions not satisfied
954 delete[] pParent;
955 return 0;
956 }
957
39d810c8 958 }
959
960 //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
961 if ( (fProcess == kPyW ||
962 fProcess == kPyZ ||
963 fProcess == kPyMbDefault ||
964 fProcess == kPyMb ||
04081a91 965 fProcess == kPyMbWithDirectPhoton ||
39d810c8 966 fProcess == kPyMbNonDiffr)
967 && (fCutOnChild == 1) ) {
968 if ( !CheckKinematicsOnChild() ) {
969 delete[] pParent;
970 return 0;
971 }
972 }
973
974
975 for (i = 0; i < np; i++) {
976 Int_t trackIt = 0;
8507138f 977 TParticle * iparticle = (TParticle *) fParticles.At(i);
39d810c8 978 kf = CheckPDGCode(iparticle->GetPdgCode());
979 Int_t ks = iparticle->GetStatusCode();
980 Int_t km = iparticle->GetFirstMother();
981 if ((ks == 1 && kf!=0 && KinematicSelection(iparticle, 0)) ||
982 (ks != 1) ||
983 (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
984 nc++;
985 if (ks == 1) trackIt = 1;
986
987 Int_t ipa = (fPythia->Version() == 6) ? (iparticle->GetFirstMother() - 1) :(iparticle->GetFirstMother()) ;
988 iparent = (ipa > -1) ? pParent[ipa] : -1;
989 if (ipa >= np) fPythia->EventListing();
990
991//
992// store track information
993 p[0] = iparticle->Px();
994 p[1] = iparticle->Py();
995 p[2] = iparticle->Pz();
996 p[3] = iparticle->Energy();
997
998
999 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
1000 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
1001 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
1002
1003 Float_t tof = fEventTime + kconv * iparticle->T();
1004
1005 PushTrack(fTrackIt*trackIt, iparent, kf,
1006 p[0], p[1], p[2], p[3],
1007 origin[0], origin[1], origin[2], tof,
1008 polar[0], polar[1], polar[2],
1009 kPPrimary, nt, 1., ks);
1010 fNprimaries++;
1011
1012
1013 //
1014 // Special Treatment to store color-flow
1015 //
3e54fb8c 1016 /*
39d810c8 1017 if (ks == 3 || ks == 13 || ks == 14) {
1018 TParticle* particle = 0;
1019 if (fStack) {
1020 particle = fStack->Particle(nt);
1021 } else {
33c3c91a 1022 particle = AliRunLoader::Instance()->Stack()->Particle(nt);
39d810c8 1023 }
3e54fb8c 1024 particle->SetFirstDaughter(fPythia->GetK(2, i));
1025 particle->SetLastDaughter(fPythia->GetK(3, i));
39d810c8 1026 }
3e54fb8c 1027 */
39d810c8 1028 KeepTrack(nt);
1029 pParent[i] = nt;
1030 SetHighWaterMark(nt);
1031
1032 } // select particle
1033 } // particle loop
1034
1035 delete[] pParent;
1036
1037 return 1;
1038}
1039
1040
1041void AliGenPythiaPlus::FinishRun()
1042{
1043// Print x-section summary
1044 fPythia->PrintStatistics();
1045
1046 if (fNev > 0.) {
1047 fQ /= fNev;
1048 fX1 /= fNev;
1049 fX2 /= fNev;
1050 }
1051
1052 printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
1053 printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
1054}
1055
1056void AliGenPythiaPlus::AdjustWeights() const
1057{
1058// Adjust the weights after generation of all events
1059//
1060 if (gAlice) {
1061 TParticle *part;
1062 Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
1063 for (Int_t i=0; i<ntrack; i++) {
1064 part= gAlice->GetMCApp()->Particle(i);
1065 part->SetWeight(part->GetWeight()*fKineBias);
1066 }
1067 }
1068}
1069
1070void AliGenPythiaPlus::SetNuclei(Int_t a1, Int_t a2)
1071{
1072// Treat protons as inside nuclei with mass numbers a1 and a2
1073
1074 fAProjectile = a1;
1075 fATarget = a2;
1076 fSetNuclei = kTRUE;
1077}
1078
1079
1080void AliGenPythiaPlus::MakeHeader()
1081{
1082//
1083// Make header for the simulated event
1084//
1085 if (gAlice) {
1086 if (gAlice->GetEvNumber()>=fDebugEventFirst &&
1087 gAlice->GetEvNumber()<=fDebugEventLast) fPythia->EventListing();
1088 }
1089
1090// Builds the event header, to be called after each event
1091 if (fHeader) delete fHeader;
1092 fHeader = new AliGenPythiaEventHeader("Pythia");
1093//
1094// Event type
1095 ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->ProcessCode());
1096//
1097// Number of trials
1098 ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
1099//
1100// Event Vertex
1101 fHeader->SetPrimaryVertex(fVertex);
1102
1103//
1104// Number of primaries
1105 fHeader->SetNProduced(fNprimaries);
1106//
1107// Jets that have triggered
1108
1109 if (fProcess == kPyJets)
1110 {
1111 Int_t ntrig, njet;
1112 Float_t jets[4][10];
1113 GetJets(njet, ntrig, jets);
1114
1115
1116 for (Int_t i = 0; i < ntrig; i++) {
1117 ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
1118 jets[3][i]);
1119 }
1120 }
1121//
1122// Copy relevant information from external header, if present.
1123//
1124 Float_t uqJet[4];
1125
1126 if (fRL) {
1127 AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
1128 for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
1129 {
1130 printf("Adding Jet %d %d \n", i, exHeader->NTriggerJets());
1131
1132
1133 exHeader->TriggerJet(i, uqJet);
1134 ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
1135 }
1136 }
1137//
1138// Store quenching parameters
1139//
1140 if (fQuench){
3e54fb8c 1141 Double_t z[4] = {0.};
1142 Double_t xp = 0.;
1143 Double_t yp = 0.;
39d810c8 1144 if (fQuench == 1) {
1145 // Pythia::Quench()
1146 fPythia->GetQuenchingParameters(xp, yp, z);
1147 } else {
1148 // Pyquen
1149 Double_t r1 = PARIMP.rb1;
1150 Double_t r2 = PARIMP.rb2;
1151 Double_t b = PARIMP.b1;
1152 Double_t r = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
1153 Double_t phi = PARIMP.psib1;
1154 xp = r * TMath::Cos(phi);
1155 yp = r * TMath::Sin(phi);
1156
1157 }
1158 ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
1159 ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
1160 }
1161//
1162// Store pt^hard
1163 ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetPtHard());
1164//
1165// Pass header
1166//
1167 AddHeader(fHeader);
1168 fHeader = 0x0;
1169}
1170
1171Bool_t AliGenPythiaPlus::CheckTrigger(TParticle* jet1, TParticle* jet2)
1172{
1173// Check the kinematic trigger condition
1174//
1175 Double_t eta[2];
1176 eta[0] = jet1->Eta();
1177 eta[1] = jet2->Eta();
1178 Double_t phi[2];
1179 phi[0] = jet1->Phi();
1180 phi[1] = jet2->Phi();
1181 Int_t pdg[2];
1182 pdg[0] = jet1->GetPdgCode();
1183 pdg[1] = jet2->GetPdgCode();
1184 Bool_t triggered = kFALSE;
1185
1186 if (fProcess == kPyJets) {
1187 Int_t njets = 0;
1188 Int_t ntrig = 0;
1189 Float_t jets[4][10];
1190//
1191// Use Pythia clustering on parton level to determine jet axis
1192//
1193 GetJets(njets, ntrig, jets);
1194
1195 if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
1196//
1197 } else {
1198 Int_t ij = 0;
1199 Int_t ig = 1;
1200 if (pdg[0] == kGamma) {
1201 ij = 1;
1202 ig = 0;
1203 }
1204 //Check eta range first...
1205 if ((eta[ij] < fEtaMaxJet && eta[ij] > fEtaMinJet) &&
1206 (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
1207 {
1208 //Eta is okay, now check phi range
1209 if ((phi[ij] < fPhiMaxJet && phi[ij] > fPhiMinJet) &&
1210 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
1211 {
1212 triggered = kTRUE;
1213 }
1214 }
1215 }
1216 return triggered;
1217}
1218
1219
1220
1221Bool_t AliGenPythiaPlus::CheckKinematicsOnChild(){
1222//
1223//Checking Kinematics on Child (status code 1, particle code ?, kin cuts
1224//
1225 Bool_t checking = kFALSE;
1226 Int_t j, kcode, ks, km;
1227 Int_t nPartAcc = 0; //number of particles in the acceptance range
1228 Int_t numberOfAcceptedParticles = 1;
1229 if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
8507138f 1230 Int_t npart = fParticles.GetEntriesFast();
39d810c8 1231
1232 for (j = 0; j<npart; j++) {
8507138f 1233 TParticle * jparticle = (TParticle *) fParticles.At(j);
39d810c8 1234 kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
1235 ks = jparticle->GetStatusCode();
1236 km = jparticle->GetFirstMother();
1237
1238 if( (ks == 1) && (kcode == fPdgCodeParticleforAcceptanceCut) && (KinematicSelection(jparticle,1)) ){
1239 nPartAcc++;
1240 }
1241 if( numberOfAcceptedParticles <= nPartAcc){
1242 checking = kTRUE;
1243 break;
1244 }
1245 }
1246
1247 return checking;
1248}
1249
1250void AliGenPythiaPlus::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
1251{
1252//
1253// Calls the Pythia jet finding algorithm to find jets in the current event
1254//
1255//
1256//
1257// Save jets
1258//
1259// Run Jet Finder
1260 fPythia->Pycell(njets);
1261 Int_t i;
1262 for (i = 0; i < njets; i++) {
1263 Float_t px, py, pz, e;
1264 fPythia->GetJet(i, px, py, pz, e);
1265 jets[0][i] = px;
1266 jets[1][i] = py;
1267 jets[2][i] = pz;
1268 jets[3][i] = e;
1269 }
1270}
1271
1272
1273void AliGenPythiaPlus::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
1274{
1275//
1276// Calls the Pythia clustering algorithm to find jets in the current event
1277//
1278 nJets = 0;
1279 nJetsTrig = 0;
1280 if (fJetReconstruction == kCluster) {
1281//
1282// Configure cluster algorithm
1283//
1284// fPythia->SetPARU(43, 2.);
1285// fPythia->SetMSTU(41, 1);
1286//
1287// Call cluster algorithm
1288//
1289 fPythia->Pyclus(nJets);
1290//
1291// Loading jets from common block
1292//
1293 } else {
1294
1295//
1296// Run Jet Finder
1297 fPythia->Pycell(nJets);
1298 }
1299
1300 Int_t i;
1301 for (i = 0; i < nJets; i++) {
1302 Float_t px, py, pz, e;
1303 fPythia->GetJet(i, px, py, pz, e);
1304 Float_t pt = TMath::Sqrt(px * px + py * py);
1305 Float_t phi = TMath::Pi() + TMath::ATan2(-py, -px);
1306 Float_t theta = TMath::ATan2(pt,pz);
1307 Float_t et = e * TMath::Sin(theta);
1308 Float_t eta = -TMath::Log(TMath::Tan(theta / 2.));
1309 if (
1310 eta > fEtaMinJet && eta < fEtaMaxJet &&
1311 phi > fPhiMinJet && phi < fPhiMaxJet &&
1312 et > fEtMinJet && et < fEtMaxJet
1313 )
1314 {
1315 jets[0][nJetsTrig] = px;
1316 jets[1][nJetsTrig] = py;
1317 jets[2][nJetsTrig] = pz;
1318 jets[3][nJetsTrig] = e;
1319 nJetsTrig++;
1320 } else {
1321 }
1322 }
1323}
1324
1325void AliGenPythiaPlus::GetSubEventTime()
1326{
1327 // Calculates time of the next subevent
1328 fEventTime = 0.;
1329 if (fEventsTime) {
1330 TArrayF &array = *fEventsTime;
1331 fEventTime = array[fCurSubEvent++];
1332 }
1333 // printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
1334 return;
1335}
1336
1337
1338
1339
1340Bool_t AliGenPythiaPlus::IsInEMCAL(Float_t phi, Float_t eta)
1341{
1342 // Is particle in EMCAL acceptance?
1343 // phi in degrees, etamin=-etamax
1344 if(phi > fEMCALMinPhi && phi < fEMCALMaxPhi &&
1345 eta < fEMCALEta )
1346 return kTRUE;
1347 else
1348 return kFALSE;
1349}
1350
1351Bool_t AliGenPythiaPlus::IsInPHOS(Float_t phi, Float_t eta)
1352{
1353 // Is particle in PHOS acceptance?
1354 // Acceptance slightly larger considered.
1355 // phi in degrees, etamin=-etamax
1356 if(phi > fPHOSMinPhi && phi < fPHOSMaxPhi &&
1357 eta < fPHOSEta )
1358 return kTRUE;
1359 else
1360 return kFALSE;
1361}
1362
1363void AliGenPythiaPlus::RotatePhi(Int_t iphcand, Bool_t& okdd)
1364{
1365 //calculate the new position random between fPHOSMinPhi and fPHOSMaxPhi
1366 Double_t phiPHOSmin = TMath::Pi()*fPHOSMinPhi/180;
1367 Double_t phiPHOSmax = TMath::Pi()*fPHOSMaxPhi/180;
1368 Double_t phiPHOS = gRandom->Uniform(phiPHOSmin,phiPHOSmax);
1369
1370 //calculate deltaphi
8507138f 1371 TParticle* ph = (TParticle *) fParticles.At(iphcand);
39d810c8 1372 Double_t phphi = ph->Phi();
1373 Double_t deltaphi = phiPHOS - phphi;
1374
1375
1376
1377 //loop for all particles and produce the phi rotation
8507138f 1378 Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
39d810c8 1379 Double_t oldphi, newphi;
1380 Double_t newVx, newVy, R, Vz, time;
1381 Double_t newPx, newPy, pt, Pz, e;
1382 for(Int_t i=0; i< np; i++) {
8507138f 1383 TParticle* iparticle = (TParticle *) fParticles.At(i);
39d810c8 1384 oldphi = iparticle->Phi();
1385 newphi = oldphi + deltaphi;
1386 if(newphi < 0) newphi = 2*TMath::Pi() + newphi; // correct angle
1387 if(newphi > 2*TMath::Pi()) newphi = newphi - 2*TMath::Pi(); // correct angle
1388
1389 R = iparticle->R();
1390 newVx = R*TMath::Cos(newphi);
1391 newVy = R*TMath::Sin(newphi);
1392 Vz = iparticle->Vz(); // don't transform
1393 time = iparticle->T(); // don't transform
1394
1395 pt = iparticle->Pt();
1396 newPx = pt*TMath::Cos(newphi);
1397 newPy = pt*TMath::Sin(newphi);
1398 Pz = iparticle->Pz(); // don't transform
1399 e = iparticle->Energy(); // don't transform
1400
1401 // apply rotation
1402 iparticle->SetProductionVertex(newVx, newVy, Vz, time);
1403 iparticle->SetMomentum(newPx, newPy, Pz, e);
1404
1405 } //end particle loop
1406
1407 // now let's check that we put correctly the candidate photon in PHOS
1408 Float_t phi = ph->Phi()*180./TMath::Pi(); //Convert to degrees
1409 Float_t eta =TMath::Abs(ph->Eta());//in calos etamin=-etamax
1410 if(IsInPHOS(phi,eta))
1411 okdd = kTRUE;
1412}
1413
1414
1415#ifdef never
1416void AliGenPythiaPlus::Streamer(TBuffer &R__b)
1417{
1418 // Stream an object of class AliGenPythia.
1419
1420 if (R__b.IsReading()) {
1421 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1422 AliGenerator::Streamer(R__b);
1423 R__b >> (Int_t&)fProcess;
1424 R__b >> (Int_t&)fStrucFunc;
1425 R__b >> (Int_t&)fForceDecay;
1426 R__b >> fEnergyCMS;
1427 R__b >> fKineBias;
1428 R__b >> fTrials;
1429 fParentSelect.Streamer(R__b);
1430 fChildSelect.Streamer(R__b);
1431 R__b >> fXsection;
1432// (AliPythia::Instance())->Streamer(R__b);
1433 R__b >> fPtHardMin;
1434 R__b >> fPtHardMax;
1435// if (fDecayer) fDecayer->Streamer(R__b);
1436 } else {
1437 R__b.WriteVersion(AliGenPythiaPlus::IsA());
1438 AliGenerator::Streamer(R__b);
1439 R__b << (Int_t)fProcess;
1440 R__b << (Int_t)fStrucFunc;
1441 R__b << (Int_t)fForceDecay;
1442 R__b << fEnergyCMS;
1443 R__b << fKineBias;
1444 R__b << fTrials;
1445 fParentSelect.Streamer(R__b);
1446 fChildSelect.Streamer(R__b);
1447 R__b << fXsection;
1448// R__b << fPythia;
1449 R__b << fPtHardMin;
1450 R__b << fPtHardMax;
1451 // fDecayer->Streamer(R__b);
1452 }
1453}
1454#endif
1455
1456
1457