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