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