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