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