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