]> 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
Fix for report #67958: Fortrans open call should open explicitly for read-only
[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 //
1016 if (ks == 3 || ks == 13 || ks == 14) {
1017 TParticle* particle = 0;
1018 if (fStack) {
1019 particle = fStack->Particle(nt);
1020 } else {
33c3c91a 1021 particle = AliRunLoader::Instance()->Stack()->Particle(nt);
39d810c8 1022 }
1023// particle->SetFirstDaughter(fPythia->GetK(2, i));
1024// particle->SetLastDaughter(fPythia->GetK(3, i));
1025 }
1026
1027 KeepTrack(nt);
1028 pParent[i] = nt;
1029 SetHighWaterMark(nt);
1030
1031 } // select particle
1032 } // particle loop
1033
1034 delete[] pParent;
1035
1036 return 1;
1037}
1038
1039
1040void AliGenPythiaPlus::FinishRun()
1041{
1042// Print x-section summary
1043 fPythia->PrintStatistics();
1044
1045 if (fNev > 0.) {
1046 fQ /= fNev;
1047 fX1 /= fNev;
1048 fX2 /= fNev;
1049 }
1050
1051 printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
1052 printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
1053}
1054
1055void AliGenPythiaPlus::AdjustWeights() const
1056{
1057// Adjust the weights after generation of all events
1058//
1059 if (gAlice) {
1060 TParticle *part;
1061 Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
1062 for (Int_t i=0; i<ntrack; i++) {
1063 part= gAlice->GetMCApp()->Particle(i);
1064 part->SetWeight(part->GetWeight()*fKineBias);
1065 }
1066 }
1067}
1068
1069void AliGenPythiaPlus::SetNuclei(Int_t a1, Int_t a2)
1070{
1071// Treat protons as inside nuclei with mass numbers a1 and a2
1072
1073 fAProjectile = a1;
1074 fATarget = a2;
1075 fSetNuclei = kTRUE;
1076}
1077
1078
1079void AliGenPythiaPlus::MakeHeader()
1080{
1081//
1082// Make header for the simulated event
1083//
1084 if (gAlice) {
1085 if (gAlice->GetEvNumber()>=fDebugEventFirst &&
1086 gAlice->GetEvNumber()<=fDebugEventLast) fPythia->EventListing();
1087 }
1088
1089// Builds the event header, to be called after each event
1090 if (fHeader) delete fHeader;
1091 fHeader = new AliGenPythiaEventHeader("Pythia");
1092//
1093// Event type
1094 ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->ProcessCode());
1095//
1096// Number of trials
1097 ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
1098//
1099// Event Vertex
1100 fHeader->SetPrimaryVertex(fVertex);
1101
1102//
1103// Number of primaries
1104 fHeader->SetNProduced(fNprimaries);
1105//
1106// Jets that have triggered
1107
1108 if (fProcess == kPyJets)
1109 {
1110 Int_t ntrig, njet;
1111 Float_t jets[4][10];
1112 GetJets(njet, ntrig, jets);
1113
1114
1115 for (Int_t i = 0; i < ntrig; i++) {
1116 ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
1117 jets[3][i]);
1118 }
1119 }
1120//
1121// Copy relevant information from external header, if present.
1122//
1123 Float_t uqJet[4];
1124
1125 if (fRL) {
1126 AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
1127 for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
1128 {
1129 printf("Adding Jet %d %d \n", i, exHeader->NTriggerJets());
1130
1131
1132 exHeader->TriggerJet(i, uqJet);
1133 ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
1134 }
1135 }
1136//
1137// Store quenching parameters
1138//
1139 if (fQuench){
1140 Double_t z[4];
1141 Double_t xp, yp;
1142 if (fQuench == 1) {
1143 // Pythia::Quench()
1144 fPythia->GetQuenchingParameters(xp, yp, z);
1145 } else {
1146 // Pyquen
1147 Double_t r1 = PARIMP.rb1;
1148 Double_t r2 = PARIMP.rb2;
1149 Double_t b = PARIMP.b1;
1150 Double_t r = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
1151 Double_t phi = PARIMP.psib1;
1152 xp = r * TMath::Cos(phi);
1153 yp = r * TMath::Sin(phi);
1154
1155 }
1156 ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
1157 ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
1158 }
1159//
1160// Store pt^hard
1161 ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetPtHard());
1162//
1163// Pass header
1164//
1165 AddHeader(fHeader);
1166 fHeader = 0x0;
1167}
1168
1169Bool_t AliGenPythiaPlus::CheckTrigger(TParticle* jet1, TParticle* jet2)
1170{
1171// Check the kinematic trigger condition
1172//
1173 Double_t eta[2];
1174 eta[0] = jet1->Eta();
1175 eta[1] = jet2->Eta();
1176 Double_t phi[2];
1177 phi[0] = jet1->Phi();
1178 phi[1] = jet2->Phi();
1179 Int_t pdg[2];
1180 pdg[0] = jet1->GetPdgCode();
1181 pdg[1] = jet2->GetPdgCode();
1182 Bool_t triggered = kFALSE;
1183
1184 if (fProcess == kPyJets) {
1185 Int_t njets = 0;
1186 Int_t ntrig = 0;
1187 Float_t jets[4][10];
1188//
1189// Use Pythia clustering on parton level to determine jet axis
1190//
1191 GetJets(njets, ntrig, jets);
1192
1193 if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
1194//
1195 } else {
1196 Int_t ij = 0;
1197 Int_t ig = 1;
1198 if (pdg[0] == kGamma) {
1199 ij = 1;
1200 ig = 0;
1201 }
1202 //Check eta range first...
1203 if ((eta[ij] < fEtaMaxJet && eta[ij] > fEtaMinJet) &&
1204 (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
1205 {
1206 //Eta is okay, now check phi range
1207 if ((phi[ij] < fPhiMaxJet && phi[ij] > fPhiMinJet) &&
1208 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
1209 {
1210 triggered = kTRUE;
1211 }
1212 }
1213 }
1214 return triggered;
1215}
1216
1217
1218
1219Bool_t AliGenPythiaPlus::CheckKinematicsOnChild(){
1220//
1221//Checking Kinematics on Child (status code 1, particle code ?, kin cuts
1222//
1223 Bool_t checking = kFALSE;
1224 Int_t j, kcode, ks, km;
1225 Int_t nPartAcc = 0; //number of particles in the acceptance range
1226 Int_t numberOfAcceptedParticles = 1;
1227 if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
8507138f 1228 Int_t npart = fParticles.GetEntriesFast();
39d810c8 1229
1230 for (j = 0; j<npart; j++) {
8507138f 1231 TParticle * jparticle = (TParticle *) fParticles.At(j);
39d810c8 1232 kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
1233 ks = jparticle->GetStatusCode();
1234 km = jparticle->GetFirstMother();
1235
1236 if( (ks == 1) && (kcode == fPdgCodeParticleforAcceptanceCut) && (KinematicSelection(jparticle,1)) ){
1237 nPartAcc++;
1238 }
1239 if( numberOfAcceptedParticles <= nPartAcc){
1240 checking = kTRUE;
1241 break;
1242 }
1243 }
1244
1245 return checking;
1246}
1247
1248void AliGenPythiaPlus::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
1249{
1250//
1251// Calls the Pythia jet finding algorithm to find jets in the current event
1252//
1253//
1254//
1255// Save jets
1256//
1257// Run Jet Finder
1258 fPythia->Pycell(njets);
1259 Int_t i;
1260 for (i = 0; i < njets; i++) {
1261 Float_t px, py, pz, e;
1262 fPythia->GetJet(i, px, py, pz, e);
1263 jets[0][i] = px;
1264 jets[1][i] = py;
1265 jets[2][i] = pz;
1266 jets[3][i] = e;
1267 }
1268}
1269
1270
1271void AliGenPythiaPlus::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
1272{
1273//
1274// Calls the Pythia clustering algorithm to find jets in the current event
1275//
1276 nJets = 0;
1277 nJetsTrig = 0;
1278 if (fJetReconstruction == kCluster) {
1279//
1280// Configure cluster algorithm
1281//
1282// fPythia->SetPARU(43, 2.);
1283// fPythia->SetMSTU(41, 1);
1284//
1285// Call cluster algorithm
1286//
1287 fPythia->Pyclus(nJets);
1288//
1289// Loading jets from common block
1290//
1291 } else {
1292
1293//
1294// Run Jet Finder
1295 fPythia->Pycell(nJets);
1296 }
1297
1298 Int_t i;
1299 for (i = 0; i < nJets; i++) {
1300 Float_t px, py, pz, e;
1301 fPythia->GetJet(i, px, py, pz, e);
1302 Float_t pt = TMath::Sqrt(px * px + py * py);
1303 Float_t phi = TMath::Pi() + TMath::ATan2(-py, -px);
1304 Float_t theta = TMath::ATan2(pt,pz);
1305 Float_t et = e * TMath::Sin(theta);
1306 Float_t eta = -TMath::Log(TMath::Tan(theta / 2.));
1307 if (
1308 eta > fEtaMinJet && eta < fEtaMaxJet &&
1309 phi > fPhiMinJet && phi < fPhiMaxJet &&
1310 et > fEtMinJet && et < fEtMaxJet
1311 )
1312 {
1313 jets[0][nJetsTrig] = px;
1314 jets[1][nJetsTrig] = py;
1315 jets[2][nJetsTrig] = pz;
1316 jets[3][nJetsTrig] = e;
1317 nJetsTrig++;
1318 } else {
1319 }
1320 }
1321}
1322
1323void AliGenPythiaPlus::GetSubEventTime()
1324{
1325 // Calculates time of the next subevent
1326 fEventTime = 0.;
1327 if (fEventsTime) {
1328 TArrayF &array = *fEventsTime;
1329 fEventTime = array[fCurSubEvent++];
1330 }
1331 // printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
1332 return;
1333}
1334
1335
1336
1337
1338Bool_t AliGenPythiaPlus::IsInEMCAL(Float_t phi, Float_t eta)
1339{
1340 // Is particle in EMCAL acceptance?
1341 // phi in degrees, etamin=-etamax
1342 if(phi > fEMCALMinPhi && phi < fEMCALMaxPhi &&
1343 eta < fEMCALEta )
1344 return kTRUE;
1345 else
1346 return kFALSE;
1347}
1348
1349Bool_t AliGenPythiaPlus::IsInPHOS(Float_t phi, Float_t eta)
1350{
1351 // Is particle in PHOS acceptance?
1352 // Acceptance slightly larger considered.
1353 // phi in degrees, etamin=-etamax
1354 if(phi > fPHOSMinPhi && phi < fPHOSMaxPhi &&
1355 eta < fPHOSEta )
1356 return kTRUE;
1357 else
1358 return kFALSE;
1359}
1360
1361void AliGenPythiaPlus::RotatePhi(Int_t iphcand, Bool_t& okdd)
1362{
1363 //calculate the new position random between fPHOSMinPhi and fPHOSMaxPhi
1364 Double_t phiPHOSmin = TMath::Pi()*fPHOSMinPhi/180;
1365 Double_t phiPHOSmax = TMath::Pi()*fPHOSMaxPhi/180;
1366 Double_t phiPHOS = gRandom->Uniform(phiPHOSmin,phiPHOSmax);
1367
1368 //calculate deltaphi
8507138f 1369 TParticle* ph = (TParticle *) fParticles.At(iphcand);
39d810c8 1370 Double_t phphi = ph->Phi();
1371 Double_t deltaphi = phiPHOS - phphi;
1372
1373
1374
1375 //loop for all particles and produce the phi rotation
8507138f 1376 Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
39d810c8 1377 Double_t oldphi, newphi;
1378 Double_t newVx, newVy, R, Vz, time;
1379 Double_t newPx, newPy, pt, Pz, e;
1380 for(Int_t i=0; i< np; i++) {
8507138f 1381 TParticle* iparticle = (TParticle *) fParticles.At(i);
39d810c8 1382 oldphi = iparticle->Phi();
1383 newphi = oldphi + deltaphi;
1384 if(newphi < 0) newphi = 2*TMath::Pi() + newphi; // correct angle
1385 if(newphi > 2*TMath::Pi()) newphi = newphi - 2*TMath::Pi(); // correct angle
1386
1387 R = iparticle->R();
1388 newVx = R*TMath::Cos(newphi);
1389 newVy = R*TMath::Sin(newphi);
1390 Vz = iparticle->Vz(); // don't transform
1391 time = iparticle->T(); // don't transform
1392
1393 pt = iparticle->Pt();
1394 newPx = pt*TMath::Cos(newphi);
1395 newPy = pt*TMath::Sin(newphi);
1396 Pz = iparticle->Pz(); // don't transform
1397 e = iparticle->Energy(); // don't transform
1398
1399 // apply rotation
1400 iparticle->SetProductionVertex(newVx, newVy, Vz, time);
1401 iparticle->SetMomentum(newPx, newPy, Pz, e);
1402
1403 } //end particle loop
1404
1405 // now let's check that we put correctly the candidate photon in PHOS
1406 Float_t phi = ph->Phi()*180./TMath::Pi(); //Convert to degrees
1407 Float_t eta =TMath::Abs(ph->Eta());//in calos etamin=-etamax
1408 if(IsInPHOS(phi,eta))
1409 okdd = kTRUE;
1410}
1411
1412
1413#ifdef never
1414void AliGenPythiaPlus::Streamer(TBuffer &R__b)
1415{
1416 // Stream an object of class AliGenPythia.
1417
1418 if (R__b.IsReading()) {
1419 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1420 AliGenerator::Streamer(R__b);
1421 R__b >> (Int_t&)fProcess;
1422 R__b >> (Int_t&)fStrucFunc;
1423 R__b >> (Int_t&)fForceDecay;
1424 R__b >> fEnergyCMS;
1425 R__b >> fKineBias;
1426 R__b >> fTrials;
1427 fParentSelect.Streamer(R__b);
1428 fChildSelect.Streamer(R__b);
1429 R__b >> fXsection;
1430// (AliPythia::Instance())->Streamer(R__b);
1431 R__b >> fPtHardMin;
1432 R__b >> fPtHardMax;
1433// if (fDecayer) fDecayer->Streamer(R__b);
1434 } else {
1435 R__b.WriteVersion(AliGenPythiaPlus::IsA());
1436 AliGenerator::Streamer(R__b);
1437 R__b << (Int_t)fProcess;
1438 R__b << (Int_t)fStrucFunc;
1439 R__b << (Int_t)fForceDecay;
1440 R__b << fEnergyCMS;
1441 R__b << fKineBias;
1442 R__b << fTrials;
1443 fParentSelect.Streamer(R__b);
1444 fChildSelect.Streamer(R__b);
1445 R__b << fXsection;
1446// R__b << fPythia;
1447 R__b << fPtHardMin;
1448 R__b << fPtHardMax;
1449 // fDecayer->Streamer(R__b);
1450 }
1451}
1452#endif
1453
1454
1455