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