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