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