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