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