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