git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
TClonesArray* repaced by TClonesArray
[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),
55 fEnergyCMS(5500.),
56 fKineBias(0.),
57 fTrials(0),
58 fTrialsRun(0),
59 fQ(0.),
60 fX1(0.),
61 fX2(0.),
62 fEventTime(0.),
63 fInteractionRate(0.),
64 fTimeWindow(0.),
65 fCurSubEvent(0),
66 fEventsTime(0),
67 fNev(0),
68 fFlavorSelect(0),
69 fXsection(0.),
70 fPtHardMin(0.),
71 fPtHardMax(1.e4),
72 fYHardMin(-1.e10),
73 fYHardMax(1.e10),
74 fGinit(1),
75 fGfinal(1),
76 fHadronisation(1),
77 fNpartons(0),
78 fReadFromFile(0),
79 fQuench(0),
80 fPtKick(1.),
81 fFullEvent(kTRUE),
82 fDecayer(new AliDecayerPythia()),
83 fDebugEventFirst(-1),
84 fDebugEventLast(-1),
85 fEtMinJet(0.),
86 fEtMaxJet(1.e4),
87 fEtaMinJet(-20.),
88 fEtaMaxJet(20.),
89 fPhiMinJet(0.),
90 fPhiMaxJet(2.* TMath::Pi()),
91 fJetReconstruction(kCell),
92 fEtaMinGamma(-20.),
93 fEtaMaxGamma(20.),
94 fPhiMinGamma(0.),
95 fPhiMaxGamma(2. * TMath::Pi()),
96 fPycellEtaMax(2.),
97 fPycellNEta(274),
98 fPycellNPhi(432),
99 fPycellThreshold(0.),
100 fPycellEtSeed(4.),
101 fPycellMinEtJet(10.),
102 fPycellMaxRadius(1.),
103 fStackFillOpt(kFlavorSelection),
104 fFeedDownOpt(kTRUE),
105 fFragmentation(kTRUE),
106 fSetNuclei(kFALSE),
107 fNewMIS(kFALSE),
108 fHFoff(kFALSE),
109 fTriggerParticle(0),
110 fTriggerEta(0.9),
111 fCountMode(kCountAll),
112 fHeader(0),
113 fRL(0),
114 fFileName(0),
115 fFragPhotonInCalo(kFALSE),
116 fPi0InCalo(kFALSE) ,
117 fPhotonInCalo(kFALSE),
118 fCheckEMCAL(kFALSE),
119 fCheckPHOS(kFALSE),
120 fCheckPHOSeta(kFALSE),
121 fFragPhotonOrPi0MinPt(0),
122 fPhotonMinPt(0),
123 fPHOSMinPhi(219.),
124 fPHOSMaxPhi(321.),
125 fPHOSEta(0.13),
126 fEMCALMinPhi(79.),
127 fEMCALMaxPhi(191.),
128 fEMCALEta(0.71)
129
130{
131// Default Constructor
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),
142 fEnergyCMS(5500.),
143 fKineBias(0.),
144 fTrials(0),
145 fTrialsRun(0),
146 fQ(0.),
147 fX1(0.),
148 fX2(0.),
149 fEventTime(0.),
150 fInteractionRate(0.),
151 fTimeWindow(0.),
152 fCurSubEvent(0),
153 fEventsTime(0),
154 fNev(0),
155 fFlavorSelect(0),
156 fXsection(0.),
157 fPtHardMin(0.),
158 fPtHardMax(1.e4),
159 fYHardMin(-1.e10),
160 fYHardMax(1.e10),
161 fGinit(kTRUE),
162 fGfinal(kTRUE),
163 fHadronisation(kTRUE),
164 fNpartons(0),
165 fReadFromFile(kFALSE),
166 fQuench(kFALSE),
167 fPtKick(1.),
168 fFullEvent(kTRUE),
169 fDecayer(new AliDecayerPythia()),
170 fDebugEventFirst(-1),
171 fDebugEventLast(-1),
172 fEtMinJet(0.),
173 fEtMaxJet(1.e4),
174 fEtaMinJet(-20.),
175 fEtaMaxJet(20.),
176 fPhiMinJet(0.),
177 fPhiMaxJet(2.* TMath::Pi()),
178 fJetReconstruction(kCell),
179 fEtaMinGamma(-20.),
180 fEtaMaxGamma(20.),
181 fPhiMinGamma(0.),
182 fPhiMaxGamma(2. * TMath::Pi()),
183 fPycellEtaMax(2.),
184 fPycellNEta(274),
185 fPycellNPhi(432),
186 fPycellThreshold(0.),
187 fPycellEtSeed(4.),
188 fPycellMinEtJet(10.),
189 fPycellMaxRadius(1.),
190 fStackFillOpt(kFlavorSelection),
191 fFeedDownOpt(kTRUE),
192 fFragmentation(kTRUE),
193 fSetNuclei(kFALSE),
194 fNewMIS(kFALSE),
195 fHFoff(kFALSE),
196 fTriggerParticle(0),
197 fTriggerEta(0.9),
198 fCountMode(kCountAll),
199 fHeader(0),
200 fRL(0),
201 fFileName(0),
202 fFragPhotonInCalo(kFALSE),
203 fPi0InCalo(kFALSE) ,
204 fPhotonInCalo(kFALSE),
205 fCheckEMCAL(kFALSE),
206 fCheckPHOS(kFALSE),
207 fCheckPHOSeta(kFALSE),
208 fFragPhotonOrPi0MinPt(0),
209 fPhotonMinPt(0),
210 fPHOSMinPhi(219.),
211 fPHOSMaxPhi(321.),
212 fPHOSEta(0.13),
213 fEMCALMinPhi(79.),
214 fEMCALMaxPhi(191.),
215 fEMCALEta(0.71)
216{
217// default charm production at 5. 5 TeV
218// semimuonic decay
219// structure function GRVHO
220//
221 fName = "Pythia";
222 fTitle= "Particle Generator using PYTHIA";
223 SetForceDecay();
224 // Set random number generator
225 if (!AliPythiaRndm::GetPythiaRandom())
226 AliPythiaRndm::SetPythiaRandom(GetRandom());
227 fParticles = new TClonesArray("TParticle",1000);
228 SetNuclei(0,0);
229 }
230
231AliGenPythiaPlus::~AliGenPythiaPlus()
232{
233// Destructor
234 if(fEventsTime) delete fEventsTime;
235}
236
237void AliGenPythiaPlus::SetInteractionRate(Float_t rate,Float_t timewindow)
238{
239// Generate pileup using user specified rate
240 fInteractionRate = rate;
241 fTimeWindow = timewindow;
242 GeneratePileup();
243}
244
245void AliGenPythiaPlus::GeneratePileup()
246{
247// Generate sub events time for pileup
248 fEventsTime = 0;
249 if(fInteractionRate == 0.) {
250 Warning("GeneratePileup","Zero interaction specified. Skipping pileup generation.\n");
251 return;
252 }
253
254 Int_t npart = NumberParticles();
255 if(npart < 0) {
256 Warning("GeneratePileup","Negative number of particles. Skipping pileup generation.\n");
257 return;
258 }
259
260 if(fEventsTime) delete fEventsTime;
261 fEventsTime = new TArrayF(npart);
262 TArrayF &array = *fEventsTime;
263 for(Int_t ipart = 0; ipart < npart; ipart++)
264 array[ipart] = 0.;
265
266 Float_t eventtime = 0.;
267 while(1)
268 {
269 eventtime += (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
270 if(eventtime > fTimeWindow) break;
271 array.Set(array.GetSize()+1);
272 array[array.GetSize()-1] = eventtime;
273 }
274
275 eventtime = 0.;
276 while(1)
277 {
278 eventtime -= (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
279 if(TMath::Abs(eventtime) > fTimeWindow) break;
280 array.Set(array.GetSize()+1);
281 array[array.GetSize()-1] = eventtime;
282 }
283
284 SetNumberParticles(fEventsTime->GetSize());
285}
286
287void AliGenPythiaPlus::SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
288 Float_t thresh, Float_t etseed, Float_t minet, Float_t r)
289{
290// Set pycell parameters
291 fPycellEtaMax = etamax;
292 fPycellNEta = neta;
293 fPycellNPhi = nphi;
294 fPycellThreshold = thresh;
295 fPycellEtSeed = etseed;
296 fPycellMinEtJet = minet;
297 fPycellMaxRadius = r;
298}
299
300
301
302void AliGenPythiaPlus::SetEventListRange(Int_t eventFirst, Int_t eventLast)
303{
304 // Set a range of event numbers, for which a table
305 // of generated particle will be printed
306 fDebugEventFirst = eventFirst;
307 fDebugEventLast = eventLast;
308 if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst;
309}
310
311void AliGenPythiaPlus::Init()
312{
313// Initialisation
314
315// SetMC(AliPythia::Instance());
316// fPythia=(AliPythia*) fMCEvGen;
317
318//
319 fParentWeight=1./Float_t(fNpart);
320//
321
322
323 fPythia->SetPtHardRange(fPtHardMin, fPtHardMax);
324 fPythia->SetYHardRange(fYHardMin, fYHardMax);
325
326 if (fAProjectile > 0 && fATarget > 0) fPythia->SetNuclei(fAProjectile, fATarget);
327 // Fragmentation?
328 if (fFragmentation) {
329 fPythia->SetFragmentation(1);
330 } else {
331 fPythia->SetFragmentation(0);
332 }
333
334
335// initial state radiation
336 fPythia->SetInitialAndFinalStateRadiation(fGinit, fGfinal);
337
338// pt - kick
339 fPythia->SetIntrinsicKt(fPtKick);
340
341 if (fReadFromFile) {
342 fRL = AliRunLoader::Open(fFileName, "Partons");
343 fRL->LoadKinematics();
344 fRL->LoadHeader();
345 } else {
346 fRL = 0x0;
347 }
348 //
349 fPythia->ProcInit(fProcess, fEnergyCMS, fStrucFunc);
350 // Forward Paramters to the AliPythia object
351 fDecayer->SetForceDecay(fForceDecay);
352// Switch off Heavy Flavors on request
353 if (fHFoff) {
354 fPythia->SwitchHFOff();
355 // Switch off g->QQbar splitting in decay table
356 ((AliDecayerPythia*) fDecayer)->HeavyFlavourOff();
357 }
358
359 fDecayer->Init();
360
361
362// Parent and Children Selection
363 switch (fProcess)
364 {
365 case kPyOldUEQ2ordered:
366 case kPyOldUEQ2ordered2:
367 case kPyOldPopcorn:
368 break;
369 case kPyCharm:
370 case kPyCharmUnforced:
371 case kPyCharmPbPbMNR:
372 case kPyCharmpPbMNR:
373 case kPyCharmppMNR:
374 case kPyCharmppMNRwmi:
375 fParentSelect[0] = 411;
376 fParentSelect[1] = 421;
377 fParentSelect[2] = 431;
378 fParentSelect[3] = 4122;
379 fFlavorSelect = 4;
380 break;
381 case kPyD0PbPbMNR:
382 case kPyD0pPbMNR:
383 case kPyD0ppMNR:
384 fParentSelect[0] = 421;
385 fFlavorSelect = 4;
386 break;
387 case kPyDPlusPbPbMNR:
388 case kPyDPluspPbMNR:
389 case kPyDPlusppMNR:
390 fParentSelect[0] = 411;
391 fFlavorSelect = 4;
392 break;
393 case kPyDPlusStrangePbPbMNR:
394 case kPyDPlusStrangepPbMNR:
395 case kPyDPlusStrangeppMNR:
396 fParentSelect[0] = 431;
397 fFlavorSelect = 4;
398 break;
399 case kPyBeauty:
400 case kPyBeautyPbPbMNR:
401 case kPyBeautypPbMNR:
402 case kPyBeautyppMNR:
403 case kPyBeautyppMNRwmi:
404 fParentSelect[0]= 511;
405 fParentSelect[1]= 521;
406 fParentSelect[2]= 531;
407 fParentSelect[3]= 5122;
408 fParentSelect[4]= 5132;
409 fParentSelect[5]= 5232;
410 fParentSelect[6]= 5332;
411 fFlavorSelect = 5;
412 break;
413 case kPyBeautyUnforced:
414 fParentSelect[0] = 511;
415 fParentSelect[1] = 521;
416 fParentSelect[2] = 531;
417 fParentSelect[3] = 5122;
418 fParentSelect[4] = 5132;
419 fParentSelect[5] = 5232;
420 fParentSelect[6] = 5332;
421 fFlavorSelect = 5;
422 break;
423 case kPyJpsiChi:
424 case kPyJpsi:
425 fParentSelect[0] = 443;
426 break;
427 case kPyMbDefault:
428 case kPyMb:
429 case kPyMbNonDiffr:
430 case kPyMbMSEL1:
431 case kPyJets:
432 case kPyDirectGamma:
433 case kPyLhwgMb:
434 break;
435 case kPyW:
436 case kPyZ:
437 break;
438 }
439//
440//
441// JetFinder for Trigger
442//
443// Configure detector (EMCAL like)
444//
445 fPythia->SetPycellParameters(fPycellEtaMax,fPycellNEta, fPycellNPhi,
446 fPycellThreshold, fPycellEtSeed,
447 fPycellMinEtJet, fPycellMaxRadius);
448//
449// This counts the total number of calls to Pyevnt() per run.
450 fTrialsRun = 0;
451 fQ = 0.;
452 fX1 = 0.;
453 fX2 = 0.;
454 fNev = 0 ;
455//
456//
457//
458 AliGenMC::Init();
459//
460//
461//
462 if (fSetNuclei) {
463 fDyBoost = 0;
464 Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n");
465 }
466
467 if (fQuench) {
468 fPythia->InitQuenching(0., 0.1, 0.6e6, 0, 0.97, 30);
469 }
470
471// fPythia->SetPARJ(200, 0.0);
472
473// if (fQuench == 3) {
474// // Nestor's change of the splittings
475// fPythia->SetPARJ(200, 0.8);
476// fPythia->SetMSTJ(41, 1); // QCD radiation only
477// fPythia->SetMSTJ(42, 2); // angular ordering
478// fPythia->SetMSTJ(44, 2); // option to run alpha_s
479// fPythia->SetMSTJ(47, 0); // No correction back to hard scattering element
480// fPythia->SetMSTJ(50, 0); // No coherence in first branching
481// fPythia->SetPARJ(82, 1.); // Cut off for parton showers
482// }
483}
484
485void AliGenPythiaPlus::Generate()
486{
487// Generate one event
488
489 fDecayer->ForceDecay();
490
491 Float_t polar[3] = {0,0,0};
492 Float_t origin[3] = {0,0,0};
493 Float_t p[4];
494// converts from mm/c to s
495 const Float_t kconv=0.001/2.999792458e8;
496//
497 Int_t nt=0;
498 Int_t jev=0;
499 Int_t j, kf;
500 fTrials=0;
501 fEventTime = 0.;
502
503
504
505 // Set collision vertex position
506 if (fVertexSmear == kPerEvent) Vertex();
507
508// event loop
509 while(1)
510 {
511//
512// Produce event
513//
514//
515// Switch hadronisation off
516//
517// fPythia->SwitchHadronisationOff();
518//
519// Either produce new event or read partons from file
520//
521 if (!fReadFromFile) {
522 if (!fNewMIS) {
523 fPythia->GenerateEvent();
524 } else {
525 fPythia->GenerateMIEvent();
526 }
527 fNpartons = fPythia->GetNumberOfParticles();
528 } else {
529 printf("Loading Event %d\n",AliRunLoader::GetRunLoader()->GetEventNumber());
530 fRL->GetEvent(AliRunLoader::GetRunLoader()->GetEventNumber());
531 fPythia->SetNumberOfParticles(0);
532 fPythia->LoadEvent(fRL->Stack(), 0 , 1);
533 fPythia->EditEventList(21);
534 }
535
536//
537// Run quenching routine
538//
539 if (fQuench == 1) {
540 fPythia->Quench();
541 } else if (fQuench == 2){
542 fPythia->Pyquen(208., 0, 0.);
543 } else if (fQuench == 3) {
544 // Quenching is via multiplicative correction of the splittings
545 }
546
547//
548// Switch hadronisation on
549//
550// fPythia->SwitchHadronisationOn();
551//
552// .. and perform hadronisation
553// printf("Calling hadronisation %d\n", fPythia->GetN());
554// fPythia->HadronizeEvent();
555 fTrials++;
556 fPythia->GetParticles(fParticles);
557 Boost();
558//
559//
560//
561 Int_t i;
562
563 fNprimaries = 0;
564 Int_t np = fParticles->GetEntriesFast();
565
566 if (np == 0) continue;
567//
568
569//
570 Int_t* pParent = new Int_t[np];
571 Int_t* pSelected = new Int_t[np];
572 Int_t* trackIt = new Int_t[np];
573 for (i = 0; i < np; i++) {
574 pParent[i] = -1;
575 pSelected[i] = 0;
576 trackIt[i] = 0;
577 }
578
579 Int_t nc = 0; // Total n. of selected particles
580 Int_t nParents = 0; // Selected parents
581 Int_t nTkbles = 0; // Trackable particles
582 if (fProcess != kPyMbDefault &&
583 fProcess != kPyMb &&
584 fProcess != kPyJets &&
585 fProcess != kPyDirectGamma &&
586 fProcess != kPyMbNonDiffr &&
587 fProcess != kPyMbMSEL1 &&
588 fProcess != kPyW &&
589 fProcess != kPyZ &&
590 fProcess != kPyCharmppMNRwmi &&
591 fProcess != kPyBeautyppMNRwmi) {
592
593 for (i = 0; i < np; i++) {
594 TParticle* iparticle = (TParticle *) fParticles->At(i);
595 Int_t ks = iparticle->GetStatusCode();
596 kf = CheckPDGCode(iparticle->GetPdgCode());
597// No initial state partons
598 if (ks==21) continue;
599//
600// Heavy Flavor Selection
601//
602 // quark ?
603 kf = TMath::Abs(kf);
604 Int_t kfl = kf;
605 // Resonance
606
607 if (kfl > 100000) kfl %= 100000;
608 if (kfl > 10000) kfl %= 10000;
609 // meson ?
610 if (kfl > 10) kfl/=100;
611 // baryon
612 if (kfl > 10) kfl/=10;
613 Int_t ipa = (fPythia->Version() == 6) ? (iparticle->GetFirstMother() - 1) :(iparticle->GetFirstMother()) ;
614 Int_t kfMo = 0;
615//
616// Establish mother daughter relation between heavy quarks and mesons
617//
618 if (kf >= fFlavorSelect && kf <= 6) {
619 Int_t idau = (fPythia->Version() == 6) ? (iparticle->GetFirstDaughter() - 1) :(iparticle->GetFirstDaughter());
620 if (idau > -1) {
621 TParticle* daughter = (TParticle *) fParticles->At(idau);
622 Int_t pdgD = daughter->GetPdgCode();
623 if (pdgD == 91 || pdgD == 92) {
624 Int_t jmin = (fPythia->Version() == 6) ? (daughter->GetFirstDaughter() - 1) : (daughter->GetFirstDaughter());
625 Int_t jmax = (fPythia->Version() == 6) ? (daughter->GetLastDaughter() - 1) : (daughter->GetLastDaughter());
626
627 for (Int_t j = jmin; j <= jmax; j++)
628 ((TParticle *) fParticles->At(j))->SetFirstMother(i+1);
629 } // is string or cluster
630 } // has daughter
631 } // heavy quark
632
633
634 if (ipa > -1) {
635 TParticle * mother = (TParticle *) fParticles->At(ipa);
636 kfMo = TMath::Abs(mother->GetPdgCode());
637 }
638
639 // What to keep in Stack?
640 Bool_t flavorOK = kFALSE;
641 Bool_t selectOK = kFALSE;
642 if (fFeedDownOpt) {
643 if (kfl >= fFlavorSelect) flavorOK = kTRUE;
644 } else {
645 if (kfl > fFlavorSelect) {
646 nc = -1;
647 break;
648 }
649 if (kfl == fFlavorSelect) flavorOK = kTRUE;
650 }
651 switch (fStackFillOpt) {
652 case kFlavorSelection:
653 selectOK = kTRUE;
654 break;
655 case kParentSelection:
656 if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
657 break;
658 }
659 if (flavorOK && selectOK) {
660//
661// Heavy flavor hadron or quark
662//
663// Kinematic seletion on final state heavy flavor mesons
664 if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
665 {
666 continue;
667 }
668 pSelected[i] = 1;
669 if (ParentSelected(kf)) ++nParents; // Update parent count
670// printf("\n particle (HF) %d %d %d", i, pSelected[i], kf);
671 } else {
672// Kinematic seletion on decay products
673 if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf)
674 && !KinematicSelection(iparticle, 1))
675 {
676 continue;
677 }
678//
679// Decay products
680// Select if mother was selected and is not tracked
681
682 if (pSelected[ipa] &&
683 !trackIt[ipa] && // mother will be tracked ?
684 kfMo != 5 && // mother is b-quark, don't store fragments
685 kfMo != 4 && // mother is c-quark, don't store fragments
686 kf != 92) // don't store string
687 {
688//
689// Semi-stable or de-selected: diselect decay products:
690//
691//
692 if (pSelected[i] == -1 || fDecayer->GetLifetime(kf) > fMaxLifeTime)
693 {
694 Int_t ipF = iparticle->GetFirstDaughter();
695 Int_t ipL = iparticle->GetLastDaughter();
696 if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
697 }
698// printf("\n particle (decay) %d %d %d", i, pSelected[i], kf);
699 pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
700 }
701 }
702 if (pSelected[i] == -1) pSelected[i] = 0;
703 if (!pSelected[i]) continue;
704 // Count quarks only if you did not include fragmentation
705 if (fFragmentation && kf <= 10) continue;
706
707 nc++;
708// Decision on tracking
709 trackIt[i] = 0;
710//
711// Track final state particle
712 if (ks == 1) trackIt[i] = 1;
713// Track semi-stable particles
714 if ((ks == 1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime)) trackIt[i] = 1;
715// Track particles selected by process if undecayed.
716 if (fForceDecay == kNoDecay) {
717 if (ParentSelected(kf)) trackIt[i] = 1;
718 } else {
719 if (ParentSelected(kf)) trackIt[i] = 0;
720 }
721 if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
722//
723//
724
725 } // particle selection loop
726 if (nc > 0) {
727 for (i = 0; i < np; i++) {
728 if (!pSelected[i]) continue;
729 TParticle * iparticle = (TParticle *) fParticles->At(i);
730 kf = CheckPDGCode(iparticle->GetPdgCode());
731 Int_t ks = iparticle->GetStatusCode();
732 p[0] = iparticle->Px();
733 p[1] = iparticle->Py();
734 p[2] = iparticle->Pz();
735 p[3] = iparticle->Energy();
736
737 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
738 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
739 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
740
741 Float_t tof = kconv*iparticle->T();
742 Int_t ipa = (fPythia->Version() == 6) ? (iparticle->GetFirstMother() - 1) :(iparticle->GetFirstMother()) ;
743 Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
744
745 PushTrack(fTrackIt*trackIt[i], iparent, kf,
746 p[0], p[1], p[2], p[3],
747 origin[0], origin[1], origin[2], tof,
748 polar[0], polar[1], polar[2],
749 kPPrimary, nt, 1., ks);
750 pParent[i] = nt;
751 KeepTrack(nt);
752 fNprimaries++;
753 } // PushTrack loop
754 }
755 } else {
756 nc = GenerateMB();
757 } // mb ?
758
759 GetSubEventTime();
760
761 delete[] pParent;
762 delete[] pSelected;
763 delete[] trackIt;
764
765 if (nc > 0) {
766 switch (fCountMode) {
767 case kCountAll:
768 // printf(" Count all \n");
769 jev += nc;
770 break;
771 case kCountParents:
772 // printf(" Count parents \n");
773 jev += nParents;
774 break;
775 case kCountTrackables:
776 // printf(" Count trackable \n");
777 jev += nTkbles;
778 break;
779 }
780 if (jev >= fNpart || fNpart == -1) {
781 fKineBias=Float_t(fNpart)/Float_t(fTrials);
782 fPythia->GetXandQ(fQ, fX1, fX2);
783 fTrialsRun += fTrials;
784 fNev++;
785 MakeHeader();
786 break;
787 }
788 }
789 } // event loop
790 SetHighWaterMark(nt);
791// Adjust weight due to kinematic selection
792// AdjustWeights();
793// get cross-section
794 fXsection = fPythia->GetXSection();
795}
796
797Int_t AliGenPythiaPlus::GenerateMB()
798{
799//
800// Min Bias selection and other global selections
801//
802 Int_t i, kf, nt, iparent;
803 Int_t nc = 0;
804 Float_t p[4];
805 Float_t polar[3] = {0,0,0};
806 Float_t origin[3] = {0,0,0};
807// converts from mm/c to s
808 const Float_t kconv = 0.001 / 2.999792458e8;
809
810 Int_t np = (fHadronisation) ? fParticles->GetEntriesFast() : fNpartons;
811
812 Int_t* pParent = new Int_t[np];
813 for (i=0; i< np; i++) pParent[i] = -1;
814 if (fProcess == kPyJets || fProcess == kPyDirectGamma) {
815 TParticle* jet1 = (TParticle *) fParticles->At(6);
816 TParticle* jet2 = (TParticle *) fParticles->At(7);
817 if (!CheckTrigger(jet1, jet2)) {
818 delete [] pParent;
819 return 0;
820 }
821 }
822
823 // Select jets with fragmentation photon or pi0 going to PHOS or EMCAL
824 if (fProcess == kPyJets && (fFragPhotonInCalo || fPi0InCalo) ) {
825
826 Bool_t ok = kFALSE;
827
828 Int_t pdg = 0;
829 if (fFragPhotonInCalo) pdg = 22 ; // Photon
830 else if (fPi0InCalo) pdg = 111 ; // Pi0
831
832 for (i=0; i< np; i++) {
833 TParticle* iparticle = (TParticle *) fParticles->At(i);
834 if(iparticle->GetStatusCode()==1 && iparticle->GetPdgCode()==pdg &&
835 iparticle->Pt() > fFragPhotonOrPi0MinPt){
836 Int_t imother = (fPythia->Version() == 6) ? (iparticle->GetFirstMother() - 1) :(iparticle->GetFirstMother()) ;
837 TParticle* pmother = (TParticle *) fParticles->At(imother);
838 if(pdg == 111 ||
839 (pdg == 22 && pmother->GetStatusCode() != 11))//No photon from hadron decay
840 {
841 Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
842 Float_t eta =TMath::Abs(iparticle->Eta());//in calos etamin=-etamax
843 if((fCheckEMCAL && IsInEMCAL(phi,eta)) ||
844 (fCheckPHOS && IsInPHOS(phi,eta)) )
845 ok =kTRUE;
846 }
847 }
848 }
849 if(!ok)
850 return 0;
851 }
852
853
854 // Select events with a photon pt > min pt going to PHOS eta acceptance or exactly PHOS eta phi
855 if ((fProcess == kPyJets || fProcess == kPyDirectGamma) && fPhotonInCalo && (fCheckPHOSeta || fCheckPHOS)){
856
857 Bool_t okd = kFALSE;
858
859 Int_t pdg = 22;
860 Int_t iphcand = -1;
861 for (i=0; i< np; i++) {
862 TParticle* iparticle = (TParticle *) fParticles->At(i);
863 Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
864 Float_t eta =TMath::Abs(iparticle->Eta());//in calos etamin=-etamax
865
866 if(iparticle->GetStatusCode() == 1
867 && iparticle->GetPdgCode() == pdg
868 && iparticle->Pt() > fPhotonMinPt
869 && eta < fPHOSEta){
870
871 // first check if the photon is in PHOS phi
872 if(IsInPHOS(phi,eta)){
873 okd = kTRUE;
874 break;
875 }
876 if(fCheckPHOSeta) iphcand = i; // candiate photon to rotate in phi
877
878 }
879 }
880
881 if(!okd && iphcand != -1) // execute rotation in phi
882 RotatePhi(iphcand,okd);
883
884 if(!okd)
885 return 0;
886 }
887
888 if (fTriggerParticle) {
889 Bool_t triggered = kFALSE;
890 for (i = 0; i < np; i++) {
891 TParticle * iparticle = (TParticle *) fParticles->At(i);
892 kf = CheckPDGCode(iparticle->GetPdgCode());
893 if (kf != fTriggerParticle) continue;
894 if (iparticle->Pt() == 0.) continue;
895 if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
896 triggered = kTRUE;
897 break;
898 }
899 if (!triggered) {
900 delete [] pParent;
901 return 0;
902 }
903 }
904
905
906 // Check if there is a ccbar or bbbar pair with at least one of the two
907 // in fYMin < y < fYMax
908 if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi) {
909 TParticle *hvq;
910 Bool_t theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
911 Float_t yQ;
912 Int_t pdgQ;
913 for(i=0; i<np; i++) {
914 hvq = (TParticle*)fParticles->At(i);
915 pdgQ = hvq->GetPdgCode();
916 if(TMath::Abs(pdgQ) != fFlavorSelect) continue;
917 if(pdgQ>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
918 yQ = 0.5*TMath::Log((hvq->Energy()+hvq->Pz()+1.e-13)/
919 (hvq->Energy()-hvq->Pz()+1.e-13));
920 if(yQ>fYMin && yQ<fYMax) inYcut=kTRUE;
921 }
922 if (!theQ || !theQbar || !inYcut) {
923 delete[] pParent;
924 return 0;
925 }
926 }
927
928 //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
929 if ( (fProcess == kPyW ||
930 fProcess == kPyZ ||
931 fProcess == kPyMbDefault ||
932 fProcess == kPyMb ||
933 fProcess == kPyMbNonDiffr)
934 && (fCutOnChild == 1) ) {
935 if ( !CheckKinematicsOnChild() ) {
936 delete[] pParent;
937 return 0;
938 }
939 }
940
941
942 for (i = 0; i < np; i++) {
943 Int_t trackIt = 0;
944 TParticle * iparticle = (TParticle *) fParticles->At(i);
945 kf = CheckPDGCode(iparticle->GetPdgCode());
946 Int_t ks = iparticle->GetStatusCode();
947 Int_t km = iparticle->GetFirstMother();
948 if ((ks == 1 && kf!=0 && KinematicSelection(iparticle, 0)) ||
949 (ks != 1) ||
950 (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
951 nc++;
952 if (ks == 1) trackIt = 1;
953
954 Int_t ipa = (fPythia->Version() == 6) ? (iparticle->GetFirstMother() - 1) :(iparticle->GetFirstMother()) ;
955 iparent = (ipa > -1) ? pParent[ipa] : -1;
956 if (ipa >= np) fPythia->EventListing();
957
958//
959// store track information
960 p[0] = iparticle->Px();
961 p[1] = iparticle->Py();
962 p[2] = iparticle->Pz();
963 p[3] = iparticle->Energy();
964
965
966 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
967 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
968 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
969
970 Float_t tof = fEventTime + kconv * iparticle->T();
971
972 PushTrack(fTrackIt*trackIt, iparent, kf,
973 p[0], p[1], p[2], p[3],
974 origin[0], origin[1], origin[2], tof,
975 polar[0], polar[1], polar[2],
976 kPPrimary, nt, 1., ks);
977 fNprimaries++;
978
979
980 //
981 // Special Treatment to store color-flow
982 //
983 if (ks == 3 || ks == 13 || ks == 14) {
984 TParticle* particle = 0;
985 if (fStack) {
986 particle = fStack->Particle(nt);
987 } else {
988 particle = gAlice->Stack()->Particle(nt);
989 }
990// particle->SetFirstDaughter(fPythia->GetK(2, i));
991// particle->SetLastDaughter(fPythia->GetK(3, i));
992 }
993
994 KeepTrack(nt);
995 pParent[i] = nt;
996 SetHighWaterMark(nt);
997
998 } // select particle
999 } // particle loop
1000
1001 delete[] pParent;
1002
1003 return 1;
1004}
1005
1006
1007void AliGenPythiaPlus::FinishRun()
1008{
1009// Print x-section summary
1010 fPythia->PrintStatistics();
1011
1012 if (fNev > 0.) {
1013 fQ /= fNev;
1014 fX1 /= fNev;
1015 fX2 /= fNev;
1016 }
1017
1018 printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
1019 printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
1020}
1021
1022void AliGenPythiaPlus::AdjustWeights() const
1023{
1024// Adjust the weights after generation of all events
1025//
1026 if (gAlice) {
1027 TParticle *part;
1028 Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
1029 for (Int_t i=0; i<ntrack; i++) {
1030 part= gAlice->GetMCApp()->Particle(i);
1031 part->SetWeight(part->GetWeight()*fKineBias);
1032 }
1033 }
1034}
1035
1036void AliGenPythiaPlus::SetNuclei(Int_t a1, Int_t a2)
1037{
1038// Treat protons as inside nuclei with mass numbers a1 and a2
1039
1040 fAProjectile = a1;
1041 fATarget = a2;
1042 fSetNuclei = kTRUE;
1043}
1044
1045
1046void AliGenPythiaPlus::MakeHeader()
1047{
1048//
1049// Make header for the simulated event
1050//
1051 if (gAlice) {
1052 if (gAlice->GetEvNumber()>=fDebugEventFirst &&
1053 gAlice->GetEvNumber()<=fDebugEventLast) fPythia->EventListing();
1054 }
1055
1056// Builds the event header, to be called after each event
1057 if (fHeader) delete fHeader;
1058 fHeader = new AliGenPythiaEventHeader("Pythia");
1059//
1060// Event type
1061 ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->ProcessCode());
1062//
1063// Number of trials
1064 ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
1065//
1066// Event Vertex
1067 fHeader->SetPrimaryVertex(fVertex);
1068
1069//
1070// Number of primaries
1071 fHeader->SetNProduced(fNprimaries);
1072//
1073// Jets that have triggered
1074
1075 if (fProcess == kPyJets)
1076 {
1077 Int_t ntrig, njet;
1078 Float_t jets[4][10];
1079 GetJets(njet, ntrig, jets);
1080
1081
1082 for (Int_t i = 0; i < ntrig; i++) {
1083 ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
1084 jets[3][i]);
1085 }
1086 }
1087//
1088// Copy relevant information from external header, if present.
1089//
1090 Float_t uqJet[4];
1091
1092 if (fRL) {
1093 AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
1094 for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
1095 {
1096 printf("Adding Jet %d %d \n", i, exHeader->NTriggerJets());
1097
1098
1099 exHeader->TriggerJet(i, uqJet);
1100 ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
1101 }
1102 }
1103//
1104// Store quenching parameters
1105//
1106 if (fQuench){
1107 Double_t z[4];
1108 Double_t xp, yp;
1109 if (fQuench == 1) {
1110 // Pythia::Quench()
1111 fPythia->GetQuenchingParameters(xp, yp, z);
1112 } else {
1113 // Pyquen
1114 Double_t r1 = PARIMP.rb1;
1115 Double_t r2 = PARIMP.rb2;
1116 Double_t b = PARIMP.b1;
1117 Double_t r = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
1118 Double_t phi = PARIMP.psib1;
1119 xp = r * TMath::Cos(phi);
1120 yp = r * TMath::Sin(phi);
1121
1122 }
1123 ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
1124 ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
1125 }
1126//
1127// Store pt^hard
1128 ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetPtHard());
1129//
1130// Pass header
1131//
1132 AddHeader(fHeader);
1133 fHeader = 0x0;
1134}
1135
1136Bool_t AliGenPythiaPlus::CheckTrigger(TParticle* jet1, TParticle* jet2)
1137{
1138// Check the kinematic trigger condition
1139//
1140 Double_t eta[2];
1141 eta[0] = jet1->Eta();
1142 eta[1] = jet2->Eta();
1143 Double_t phi[2];
1144 phi[0] = jet1->Phi();
1145 phi[1] = jet2->Phi();
1146 Int_t pdg[2];
1147 pdg[0] = jet1->GetPdgCode();
1148 pdg[1] = jet2->GetPdgCode();
1149 Bool_t triggered = kFALSE;
1150
1151 if (fProcess == kPyJets) {
1152 Int_t njets = 0;
1153 Int_t ntrig = 0;
1154 Float_t jets[4][10];
1155//
1156// Use Pythia clustering on parton level to determine jet axis
1157//
1158 GetJets(njets, ntrig, jets);
1159
1160 if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
1161//
1162 } else {
1163 Int_t ij = 0;
1164 Int_t ig = 1;
1165 if (pdg[0] == kGamma) {
1166 ij = 1;
1167 ig = 0;
1168 }
1169 //Check eta range first...
1170 if ((eta[ij] < fEtaMaxJet && eta[ij] > fEtaMinJet) &&
1171 (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
1172 {
1173 //Eta is okay, now check phi range
1174 if ((phi[ij] < fPhiMaxJet && phi[ij] > fPhiMinJet) &&
1175 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
1176 {
1177 triggered = kTRUE;
1178 }
1179 }
1180 }
1181 return triggered;
1182}
1183
1184
1185
1186Bool_t AliGenPythiaPlus::CheckKinematicsOnChild(){
1187//
1188//Checking Kinematics on Child (status code 1, particle code ?, kin cuts
1189//
1190 Bool_t checking = kFALSE;
1191 Int_t j, kcode, ks, km;
1192 Int_t nPartAcc = 0; //number of particles in the acceptance range
1193 Int_t numberOfAcceptedParticles = 1;
1194 if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
1195 Int_t npart = fParticles->GetEntriesFast();
1196
1197 for (j = 0; j<npart; j++) {
1198 TParticle * jparticle = (TParticle *) fParticles->At(j);
1199 kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
1200 ks = jparticle->GetStatusCode();
1201 km = jparticle->GetFirstMother();
1202
1203 if( (ks == 1) && (kcode == fPdgCodeParticleforAcceptanceCut) && (KinematicSelection(jparticle,1)) ){
1204 nPartAcc++;
1205 }
1206 if( numberOfAcceptedParticles <= nPartAcc){
1207 checking = kTRUE;
1208 break;
1209 }
1210 }
1211
1212 return checking;
1213}
1214
1215void AliGenPythiaPlus::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
1216{
1217//
1218// Calls the Pythia jet finding algorithm to find jets in the current event
1219//
1220//
1221//
1222// Save jets
1223//
1224// Run Jet Finder
1225 fPythia->Pycell(njets);
1226 Int_t i;
1227 for (i = 0; i < njets; i++) {
1228 Float_t px, py, pz, e;
1229 fPythia->GetJet(i, px, py, pz, e);
1230 jets[0][i] = px;
1231 jets[1][i] = py;
1232 jets[2][i] = pz;
1233 jets[3][i] = e;
1234 }
1235}
1236
1237
1238void AliGenPythiaPlus::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
1239{
1240//
1241// Calls the Pythia clustering algorithm to find jets in the current event
1242//
1243 nJets = 0;
1244 nJetsTrig = 0;
1245 if (fJetReconstruction == kCluster) {
1246//
1247// Configure cluster algorithm
1248//
1249// fPythia->SetPARU(43, 2.);
1250// fPythia->SetMSTU(41, 1);
1251//
1252// Call cluster algorithm
1253//
1254 fPythia->Pyclus(nJets);
1255//
1256// Loading jets from common block
1257//
1258 } else {
1259
1260//
1261// Run Jet Finder
1262 fPythia->Pycell(nJets);
1263 }
1264
1265 Int_t i;
1266 for (i = 0; i < nJets; i++) {
1267 Float_t px, py, pz, e;
1268 fPythia->GetJet(i, px, py, pz, e);
1269 Float_t pt = TMath::Sqrt(px * px + py * py);
1270 Float_t phi = TMath::Pi() + TMath::ATan2(-py, -px);
1271 Float_t theta = TMath::ATan2(pt,pz);
1272 Float_t et = e * TMath::Sin(theta);
1273 Float_t eta = -TMath::Log(TMath::Tan(theta / 2.));
1274 if (
1275 eta > fEtaMinJet && eta < fEtaMaxJet &&
1276 phi > fPhiMinJet && phi < fPhiMaxJet &&
1277 et > fEtMinJet && et < fEtMaxJet
1278 )
1279 {
1280 jets[0][nJetsTrig] = px;
1281 jets[1][nJetsTrig] = py;
1282 jets[2][nJetsTrig] = pz;
1283 jets[3][nJetsTrig] = e;
1284 nJetsTrig++;
1285 } else {
1286 }
1287 }
1288}
1289
1290void AliGenPythiaPlus::GetSubEventTime()
1291{
1292 // Calculates time of the next subevent
1293 fEventTime = 0.;
1294 if (fEventsTime) {
1295 TArrayF &array = *fEventsTime;
1296 fEventTime = array[fCurSubEvent++];
1297 }
1298 // printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
1299 return;
1300}
1301
1302
1303
1304
1305Bool_t AliGenPythiaPlus::IsInEMCAL(Float_t phi, Float_t eta)
1306{
1307 // Is particle in EMCAL acceptance?
1308 // phi in degrees, etamin=-etamax
1309 if(phi > fEMCALMinPhi && phi < fEMCALMaxPhi &&
1310 eta < fEMCALEta )
1311 return kTRUE;
1312 else
1313 return kFALSE;
1314}
1315
1316Bool_t AliGenPythiaPlus::IsInPHOS(Float_t phi, Float_t eta)
1317{
1318 // Is particle in PHOS acceptance?
1319 // Acceptance slightly larger considered.
1320 // phi in degrees, etamin=-etamax
1321 if(phi > fPHOSMinPhi && phi < fPHOSMaxPhi &&
1322 eta < fPHOSEta )
1323 return kTRUE;
1324 else
1325 return kFALSE;
1326}
1327
1328void AliGenPythiaPlus::RotatePhi(Int_t iphcand, Bool_t& okdd)
1329{
1330 //calculate the new position random between fPHOSMinPhi and fPHOSMaxPhi
1331 Double_t phiPHOSmin = TMath::Pi()*fPHOSMinPhi/180;
1332 Double_t phiPHOSmax = TMath::Pi()*fPHOSMaxPhi/180;
1333 Double_t phiPHOS = gRandom->Uniform(phiPHOSmin,phiPHOSmax);
1334
1335 //calculate deltaphi
1336 TParticle* ph = (TParticle *) fParticles->At(iphcand);
1337 Double_t phphi = ph->Phi();
1338 Double_t deltaphi = phiPHOS - phphi;
1339
1340
1341
1342 //loop for all particles and produce the phi rotation
1343 Int_t np = (fHadronisation) ? fParticles->GetEntriesFast() : fNpartons;
1344 Double_t oldphi, newphi;
1345 Double_t newVx, newVy, R, Vz, time;
1346 Double_t newPx, newPy, pt, Pz, e;
1347 for(Int_t i=0; i< np; i++) {
1348 TParticle* iparticle = (TParticle *) fParticles->At(i);
1349 oldphi = iparticle->Phi();
1350 newphi = oldphi + deltaphi;
1351 if(newphi < 0) newphi = 2*TMath::Pi() + newphi; // correct angle
1352 if(newphi > 2*TMath::Pi()) newphi = newphi - 2*TMath::Pi(); // correct angle
1353
1354 R = iparticle->R();
1355 newVx = R*TMath::Cos(newphi);
1356 newVy = R*TMath::Sin(newphi);
1357 Vz = iparticle->Vz(); // don't transform
1358 time = iparticle->T(); // don't transform
1359
1360 pt = iparticle->Pt();
1361 newPx = pt*TMath::Cos(newphi);
1362 newPy = pt*TMath::Sin(newphi);
1363 Pz = iparticle->Pz(); // don't transform
1364 e = iparticle->Energy(); // don't transform
1365
1366 // apply rotation
1367 iparticle->SetProductionVertex(newVx, newVy, Vz, time);
1368 iparticle->SetMomentum(newPx, newPy, Pz, e);
1369
1370 } //end particle loop
1371
1372 // now let's check that we put correctly the candidate photon in PHOS
1373 Float_t phi = ph->Phi()*180./TMath::Pi(); //Convert to degrees
1374 Float_t eta =TMath::Abs(ph->Eta());//in calos etamin=-etamax
1375 if(IsInPHOS(phi,eta))
1376 okdd = kTRUE;
1377}
1378
1379
1380#ifdef never
1381void AliGenPythiaPlus::Streamer(TBuffer &R__b)
1382{
1383 // Stream an object of class AliGenPythia.
1384
1385 if (R__b.IsReading()) {
1386 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1387 AliGenerator::Streamer(R__b);
1388 R__b >> (Int_t&)fProcess;
1389 R__b >> (Int_t&)fStrucFunc;
1390 R__b >> (Int_t&)fForceDecay;
1391 R__b >> fEnergyCMS;
1392 R__b >> fKineBias;
1393 R__b >> fTrials;
1394 fParentSelect.Streamer(R__b);
1395 fChildSelect.Streamer(R__b);
1396 R__b >> fXsection;
1397// (AliPythia::Instance())->Streamer(R__b);
1398 R__b >> fPtHardMin;
1399 R__b >> fPtHardMax;
1400// if (fDecayer) fDecayer->Streamer(R__b);
1401 } else {
1402 R__b.WriteVersion(AliGenPythiaPlus::IsA());
1403 AliGenerator::Streamer(R__b);
1404 R__b << (Int_t)fProcess;
1405 R__b << (Int_t)fStrucFunc;
1406 R__b << (Int_t)fForceDecay;
1407 R__b << fEnergyCMS;
1408 R__b << fKineBias;
1409 R__b << fTrials;
1410 fParentSelect.Streamer(R__b);
1411 fChildSelect.Streamer(R__b);
1412 R__b << fXsection;
1413// R__b << fPythia;
1414 R__b << fPtHardMin;
1415 R__b << fPtHardMax;
1416 // fDecayer->Streamer(R__b);
1417 }
1418}
1419#endif
1420
1421
1422