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