]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PYTHIA6/AliGenPythia.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
bugfix: propagate CF output to ESD converter in order to fill the MC information
[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
943 Float_t pt, eta, phi;
944 for (i=0; i< np; i++) {
945 TParticle* iparticle = (TParticle *) fParticles.At(i);
946 if(iparticle->GetStatusCode()==1 && TMath::Abs(iparticle->GetPdgCode())==pdg &&
947 iparticle->Pt() > fElectronMinPt){
948 pt = iparticle->Pt();
949 phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
950 eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
951 if(IsInEMCAL(phi,eta))
952 ok =kTRUE;
953 }
954 }
955 if(!ok)
956 return 0;
957 AliDebug(5,Form("Found an electron jet (pt,eta,phi) = (%f,%f,%f)",pt,eta,phi));
958 }
ec2c406e 959
700b9416 960 // Check for minimum multiplicity
961 if (fTriggerMultiplicity > 0) {
962 Int_t multiplicity = 0;
963 for (i = 0; i < np; i++) {
8507138f 964 TParticle * iparticle = (TParticle *) fParticles.At(i);
700b9416 965
966 Int_t statusCode = iparticle->GetStatusCode();
967
968 // Initial state particle
969 if (statusCode > 20)
970 continue;
971
972 // skip quarks and gluons
973 Int_t pdgCode = TMath::Abs(iparticle->GetPdgCode());
974 if (pdgCode <= 10 || pdgCode == 21)
975 continue;
976
977 if (fTriggerMultiplicityEta > 0 && TMath::Abs(iparticle->Eta()) > fTriggerMultiplicityEta)
978 continue;
979
980 TParticlePDG* pdgPart = iparticle->GetPDG();
981 if (pdgPart && pdgPart->Charge() == 0)
982 continue;
983
984 ++multiplicity;
985 }
986
987 if (multiplicity < fTriggerMultiplicity) {
988 delete [] pParent;
989 return 0;
990 }
991
992 Printf("Triggered on event with multiplicity of %d > %d", multiplicity, fTriggerMultiplicity);
993 }
90a236ce 994
995 // Select events with a photon pt > min pt going to PHOS eta acceptance or exactly PHOS eta phi
996 if ((fProcess == kPyJets || fProcess == kPyDirectGamma) && fPhotonInCalo && (fCheckPHOSeta || fCheckPHOS)){
997
998 Bool_t okd = kFALSE;
999
1000 Int_t pdg = 22;
1001 Int_t iphcand = -1;
1002 for (i=0; i< np; i++) {
8507138f 1003 TParticle* iparticle = (TParticle *) fParticles.At(i);
90a236ce 1004 Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
1005 Float_t eta =TMath::Abs(iparticle->Eta());//in calos etamin=-etamax
1006
1007 if(iparticle->GetStatusCode() == 1
1008 && iparticle->GetPdgCode() == pdg
1009 && iparticle->Pt() > fPhotonMinPt
04a26b0a 1010 && eta < fPHOSEta){
90a236ce 1011
1012 // first check if the photon is in PHOS phi
1013 if(IsInPHOS(phi,eta)){
1014 okd = kTRUE;
1015 break;
1016 }
1017 if(fCheckPHOSeta) iphcand = i; // candiate photon to rotate in phi
1018
1019 }
1020 }
1021
1022 if(!okd && iphcand != -1) // execute rotation in phi
1023 RotatePhi(iphcand,okd);
1024
1025 if(!okd)
1026 return 0;
1027 }
1028
7ce1321b 1029 if (fTriggerParticle) {
1030 Bool_t triggered = kFALSE;
1031 for (i = 0; i < np; i++) {
8507138f 1032 TParticle * iparticle = (TParticle *) fParticles.At(i);
7ce1321b 1033 kf = CheckPDGCode(iparticle->GetPdgCode());
38db0ee6 1034 if (kf != fTriggerParticle) continue;
7ce1321b 1035 if (iparticle->Pt() == 0.) continue;
1036 if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
1037 triggered = kTRUE;
1038 break;
1039 }
234f6d32 1040 if (!triggered) {
1041 delete [] pParent;
1042 return 0;
1043 }
7ce1321b 1044 }
e0e89f40 1045
1046
1047 // Check if there is a ccbar or bbbar pair with at least one of the two
1048 // in fYMin < y < fYMax
9dfe63b3 1049 //
1050 // TO BE CHECKED: Should we require this for beauty jets?
1051 //
1052 if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi || fProcess == kPyBeautyJets) {
9ccc3d0e 1053 TParticle *partCheck;
1054 TParticle *mother;
e0e89f40 1055 Bool_t theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
9ccc3d0e 1056 Bool_t theChild=kFALSE;
1057 Float_t y;
1058 Int_t pdg,mpdg,mpdgUpperFamily;
e0e89f40 1059 for(i=0; i<np; i++) {
9ccc3d0e 1060 partCheck = (TParticle*)fParticles.At(i);
1061 pdg = partCheck->GetPdgCode();
1062 if(TMath::Abs(pdg) == fFlavorSelect) { // quark
1063 if(pdg>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
1064 y = 0.5*TMath::Log((partCheck->Energy()+partCheck->Pz()+1.e-13)/
1065 (partCheck->Energy()-partCheck->Pz()+1.e-13));
1066 if(y>fYMin && y<fYMax) inYcut=kTRUE;
1067 }
1068 if(fCutOnChild && TMath::Abs(pdg) == fPdgCodeParticleforAcceptanceCut) {
1069 Int_t mi = partCheck->GetFirstMother() - 1;
1070 if(mi<0) continue;
1071 mother = (TParticle*)fParticles.At(mi);
1072 mpdg=TMath::Abs(mother->GetPdgCode());
75fde694 1073 mpdgUpperFamily=(mpdg>1000 ? mpdg+1000 : mpdg+100); // keep e from c from b
9ccc3d0e 1074 if ( ParentSelected(mpdg) ||
1075 (fFlavorSelect==5 && ParentSelected(mpdgUpperFamily))) {
1076 if (KinematicSelection(partCheck,1)) {
1077 theChild=kTRUE;
1078 }
1079 }
1080 }
e0e89f40 1081 }
9ccc3d0e 1082 if (!theQ || !theQbar || !inYcut) { // one of the c/b conditions not satisfied
234f6d32 1083 delete[] pParent;
e0e89f40 1084 return 0;
1085 }
9ccc3d0e 1086 if (fCutOnChild && !theChild) { // one of the child conditions not satisfied
1087 delete[] pParent;
1088 return 0;
1089 }
1090
e0e89f40 1091 }
aea21c57 1092
0f6ee828 1093 //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
f529e69b 1094 if ( (fProcess == kPyW ||
1095 fProcess == kPyZ ||
1096 fProcess == kPyMbDefault ||
1097 fProcess == kPyMb ||
04081a91 1098 fProcess == kPyMbWithDirectPhoton ||
f529e69b 1099 fProcess == kPyMbNonDiffr)
0f6ee828 1100 && (fCutOnChild == 1) ) {
1101 if ( !CheckKinematicsOnChild() ) {
234f6d32 1102 delete[] pParent;
0f6ee828 1103 return 0;
1104 }
aea21c57 1105 }
1106
1107
f913ec4f 1108 for (i = 0; i < np; i++) {
8d2cd130 1109 Int_t trackIt = 0;
8507138f 1110 TParticle * iparticle = (TParticle *) fParticles.At(i);
8d2cd130 1111 kf = CheckPDGCode(iparticle->GetPdgCode());
1112 Int_t ks = iparticle->GetStatusCode();
1113 Int_t km = iparticle->GetFirstMother();
1114 if ((ks == 1 && kf!=0 && KinematicSelection(iparticle, 0)) ||
1115 (ks != 1) ||
9dfe63b3 1116 //
1117 //TO BE CHECKED: Should we require this for beauty jets?
1118 //
1119 ((fProcess == kPyJets || fProcess == kPyBeautyJets) && ks == 21 && km == 0 && i>1)) {
8d2cd130 1120 nc++;
1121 if (ks == 1) trackIt = 1;
1122 Int_t ipa = iparticle->GetFirstMother()-1;
1123
1124 iparent = (ipa > -1) ? pParent[ipa] : -1;
1125
1126//
1127// store track information
1128 p[0] = iparticle->Px();
1129 p[1] = iparticle->Py();
1130 p[2] = iparticle->Pz();
a920faf9 1131 p[3] = iparticle->Energy();
1406f599 1132
a920faf9 1133
2590ccf9 1134 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
1135 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
1136 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
1137
f913ec4f 1138 Float_t tof = fEventTime + kconv * iparticle->T();
a920faf9 1139
1140 PushTrack(fTrackIt*trackIt, iparent, kf,
1141 p[0], p[1], p[2], p[3],
1142 origin[0], origin[1], origin[2], tof,
1143 polar[0], polar[1], polar[2],
1144 kPPrimary, nt, 1., ks);
dbd64db6 1145 fNprimaries++;
8d2cd130 1146 KeepTrack(nt);
1147 pParent[i] = nt;
f913ec4f 1148 SetHighWaterMark(nt);
1149
8d2cd130 1150 } // select particle
1151 } // particle loop
1152
234f6d32 1153 delete[] pParent;
e0e89f40 1154
f913ec4f 1155 return 1;
8d2cd130 1156}
1157
1158
1159void AliGenPythia::FinishRun()
1160{
1161// Print x-section summary
1162 fPythia->Pystat(1);
2779fc64 1163
1164 if (fNev > 0.) {
1165 fQ /= fNev;
1166 fX1 /= fNev;
1167 fX2 /= fNev;
1168 }
1169
8d2cd130 1170 printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
1171 printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
8d2cd130 1172}
1173
7184e472 1174void AliGenPythia::AdjustWeights() const
8d2cd130 1175{
1176// Adjust the weights after generation of all events
1177//
e2bddf81 1178 if (gAlice) {
1179 TParticle *part;
1180 Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
1181 for (Int_t i=0; i<ntrack; i++) {
1182 part= gAlice->GetMCApp()->Particle(i);
1183 part->SetWeight(part->GetWeight()*fKineBias);
1184 }
8d2cd130 1185 }
1186}
1187
20e47f08 1188void AliGenPythia::SetNuclei(Int_t a1, Int_t a2, Int_t pdfset)
8d2cd130 1189{
1190// Treat protons as inside nuclei with mass numbers a1 and a2
1d568bc2 1191
1a626d4e 1192 fAProjectile = a1;
1193 fATarget = a2;
20e47f08 1194 fNucPdf = pdfset; // 0 EKS98 1 EPS08
1d568bc2 1195 fSetNuclei = kTRUE;
8d2cd130 1196}
1197
1198
1199void AliGenPythia::MakeHeader()
1200{
7184e472 1201//
1202// Make header for the simulated event
1203//
183a5ca9 1204 if (gAlice) {
1205 if (gAlice->GetEvNumber()>=fDebugEventFirst &&
f913ec4f 1206 gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(2);
183a5ca9 1207 }
1208
8d2cd130 1209// Builds the event header, to be called after each event
e5c87a3d 1210 if (fHeader) delete fHeader;
1211 fHeader = new AliGenPythiaEventHeader("Pythia");
8d2cd130 1212//
1213// Event type
e5c87a3d 1214 ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
8d2cd130 1215//
1216// Number of trials
e5c87a3d 1217 ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
8d2cd130 1218//
1219// Event Vertex
d25cfd65 1220 fHeader->SetPrimaryVertex(fVertex);
dbd64db6 1221
1222//
1223// Number of primaries
1224 fHeader->SetNProduced(fNprimaries);
8d2cd130 1225//
1226// Jets that have triggered
f913ec4f 1227
9dfe63b3 1228 //Need to store jets for b-jet studies too!
1229 if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets)
8d2cd130 1230 {
1231 Int_t ntrig, njet;
1232 Float_t jets[4][10];
1233 GetJets(njet, ntrig, jets);
9ff6c04c 1234
8d2cd130 1235
1236 for (Int_t i = 0; i < ntrig; i++) {
e5c87a3d 1237 ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
8d2cd130 1238 jets[3][i]);
1239 }
1240 }
5fa4b20b 1241//
1242// Copy relevant information from external header, if present.
1243//
1244 Float_t uqJet[4];
1245
1246 if (fRL) {
1247 AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
1248 for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
1249 {
1250 printf("Adding Jet %d %d \n", i, exHeader->NTriggerJets());
1251
1252
1253 exHeader->TriggerJet(i, uqJet);
1254 ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
1255 }
1256 }
1257//
1258// Store quenching parameters
1259//
1260 if (fQuench){
1261 Double_t z[4];
1262 Double_t xp, yp;
7c21f297 1263 if (fQuench == 1) {
1264 // Pythia::Quench()
1265 fPythia->GetQuenchingParameters(xp, yp, z);
1bab4b79 1266 } else if (fQuench == 2){
7c21f297 1267 // Pyquen
1268 Double_t r1 = PARIMP.rb1;
1269 Double_t r2 = PARIMP.rb2;
1270 Double_t b = PARIMP.b1;
1271 Double_t r = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
1272 Double_t phi = PARIMP.psib1;
1273 xp = r * TMath::Cos(phi);
1274 yp = r * TMath::Sin(phi);
1275
1bab4b79 1276 } else if (fQuench == 4) {
1277 // QPythia
5831e75f 1278 Double_t xy[2];
e9719084 1279 Double_t i0i1[2];
1280 AliFastGlauber::Instance()->GetSavedXY(xy);
1281 AliFastGlauber::Instance()->GetSavedI0I1(i0i1);
5831e75f 1282 xp = xy[0];
1283 yp = xy[1];
e9719084 1284 ((AliGenPythiaEventHeader*) fHeader)->SetInMediumLength(2. * i0i1[1] / i0i1[0]);
7c21f297 1285 }
1bab4b79 1286
7c21f297 1287 ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
1288 ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
1bab4b79 1289 }
beac474c 1290//
1291// Store pt^hard
1292 ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
5fa4b20b 1293//
cf57b268 1294// Pass header
5fa4b20b 1295//
cf57b268 1296 AddHeader(fHeader);
4c4eac97 1297 fHeader = 0x0;
8d2cd130 1298}
cf57b268 1299
8d2cd130 1300Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
1301{
1302// Check the kinematic trigger condition
1303//
1304 Double_t eta[2];
1305 eta[0] = jet1->Eta();
1306 eta[1] = jet2->Eta();
1307 Double_t phi[2];
1308 phi[0] = jet1->Phi();
1309 phi[1] = jet2->Phi();
1310 Int_t pdg[2];
1311 pdg[0] = jet1->GetPdgCode();
1312 pdg[1] = jet2->GetPdgCode();
1313 Bool_t triggered = kFALSE;
1314
9dfe63b3 1315 //
1316 //TO BE CHECKED: If we call this method for kPyBeautyJets, we need it here
1317 //
1318 if (fProcess == kPyJets || fProcess == kPyBeautyJets) {
8d2cd130 1319 Int_t njets = 0;
1320 Int_t ntrig = 0;
1321 Float_t jets[4][10];
1322//
1323// Use Pythia clustering on parton level to determine jet axis
1324//
1325 GetJets(njets, ntrig, jets);
1326
76d6ba9a 1327 if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
8d2cd130 1328//
1329 } else {
1330 Int_t ij = 0;
1331 Int_t ig = 1;
1332 if (pdg[0] == kGamma) {
1333 ij = 1;
1334 ig = 0;
1335 }
1336 //Check eta range first...
1337 if ((eta[ij] < fEtaMaxJet && eta[ij] > fEtaMinJet) &&
1338 (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
1339 {
1340 //Eta is okay, now check phi range
1341 if ((phi[ij] < fPhiMaxJet && phi[ij] > fPhiMinJet) &&
1342 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
1343 {
1344 triggered = kTRUE;
1345 }
1346 }
1347 }
1348 return triggered;
1349}
aea21c57 1350
1351
aea21c57 1352
7184e472 1353Bool_t AliGenPythia::CheckKinematicsOnChild(){
1354//
1355//Checking Kinematics on Child (status code 1, particle code ?, kin cuts
1356//
aea21c57 1357 Bool_t checking = kFALSE;
1358 Int_t j, kcode, ks, km;
1359 Int_t nPartAcc = 0; //number of particles in the acceptance range
1360 Int_t numberOfAcceptedParticles = 1;
1361 if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
8507138f 1362 Int_t npart = fParticles.GetEntriesFast();
aea21c57 1363
0f6ee828 1364 for (j = 0; j<npart; j++) {
8507138f 1365 TParticle * jparticle = (TParticle *) fParticles.At(j);
aea21c57 1366 kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
1367 ks = jparticle->GetStatusCode();
1368 km = jparticle->GetFirstMother();
1369
1370 if( (ks == 1) && (kcode == fPdgCodeParticleforAcceptanceCut) && (KinematicSelection(jparticle,1)) ){
1371 nPartAcc++;
1372 }
0f6ee828 1373 if( numberOfAcceptedParticles <= nPartAcc){
1374 checking = kTRUE;
1375 break;
1376 }
aea21c57 1377 }
0f6ee828 1378
aea21c57 1379 return checking;
aea21c57 1380}
1381
5fa4b20b 1382void AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
8d2cd130 1383{
1058d9df 1384 //
1385 // Load event into Pythia Common Block
1386 //
1387
1388 Int_t npart = stack -> GetNprimary();
1389 Int_t n0 = 0;
1390
1391 if (!flag) {
1392 (fPythia->GetPyjets())->N = npart;
1393 } else {
1394 n0 = (fPythia->GetPyjets())->N;
1395 (fPythia->GetPyjets())->N = n0 + npart;
1396 }
1397
1398
1399 for (Int_t part = 0; part < npart; part++) {
1400 TParticle *mPart = stack->Particle(part);
32d6ef7d 1401
1058d9df 1402 Int_t kf = mPart->GetPdgCode();
1403 Int_t ks = mPart->GetStatusCode();
1404 Int_t idf = mPart->GetFirstDaughter();
1405 Int_t idl = mPart->GetLastDaughter();
1406
1407 if (reHadr) {
1408 if (ks == 11 || ks == 12) {
1409 ks -= 10;
1410 idf = -1;
1411 idl = -1;
1412 }
32d6ef7d 1413 }
1414
1058d9df 1415 Float_t px = mPart->Px();
1416 Float_t py = mPart->Py();
1417 Float_t pz = mPart->Pz();
1418 Float_t e = mPart->Energy();
1419 Float_t m = mPart->GetCalcMass();
32d6ef7d 1420
1058d9df 1421
1422 (fPythia->GetPyjets())->P[0][part+n0] = px;
1423 (fPythia->GetPyjets())->P[1][part+n0] = py;
1424 (fPythia->GetPyjets())->P[2][part+n0] = pz;
1425 (fPythia->GetPyjets())->P[3][part+n0] = e;
1426 (fPythia->GetPyjets())->P[4][part+n0] = m;
1427
1428 (fPythia->GetPyjets())->K[1][part+n0] = kf;
1429 (fPythia->GetPyjets())->K[0][part+n0] = ks;
1430 (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1431 (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1432 (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1433 }
1434}
1435
1436void AliGenPythia::LoadEvent(TObjArray* stack, Int_t flag, Int_t reHadr)
1437{
1438 //
1439 // Load event into Pythia Common Block
1440 //
1441
1442 Int_t npart = stack -> GetEntries();
1443 Int_t n0 = 0;
1444
1445 if (!flag) {
1446 (fPythia->GetPyjets())->N = npart;
1447 } else {
1448 n0 = (fPythia->GetPyjets())->N;
1449 (fPythia->GetPyjets())->N = n0 + npart;
1450 }
1451
1452
1453 for (Int_t part = 0; part < npart; part++) {
1454 TParticle *mPart = dynamic_cast<TParticle *>(stack->At(part));
1455 Int_t kf = mPart->GetPdgCode();
1456 Int_t ks = mPart->GetStatusCode();
1457 Int_t idf = mPart->GetFirstDaughter();
1458 Int_t idl = mPart->GetLastDaughter();
1459
1460 if (reHadr) {
5fa4b20b 1461 if (ks == 11 || ks == 12) {
1058d9df 1462 ks -= 10;
1463 idf = -1;
1464 idl = -1;
5fa4b20b 1465 }
8d2cd130 1466 }
1058d9df 1467
1468 Float_t px = mPart->Px();
1469 Float_t py = mPart->Py();
1470 Float_t pz = mPart->Pz();
1471 Float_t e = mPart->Energy();
1472 Float_t m = mPart->GetCalcMass();
1473
1474
1475 (fPythia->GetPyjets())->P[0][part+n0] = px;
1476 (fPythia->GetPyjets())->P[1][part+n0] = py;
1477 (fPythia->GetPyjets())->P[2][part+n0] = pz;
1478 (fPythia->GetPyjets())->P[3][part+n0] = e;
1479 (fPythia->GetPyjets())->P[4][part+n0] = m;
1480
1481 (fPythia->GetPyjets())->K[1][part+n0] = kf;
1482 (fPythia->GetPyjets())->K[0][part+n0] = ks;
1483 (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1484 (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1485 (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1486 }
8d2cd130 1487}
1488
5fa4b20b 1489
014a9521 1490void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
8d2cd130 1491{
1492//
1493// Calls the Pythia jet finding algorithm to find jets in the current event
1494//
1495//
8d2cd130 1496//
1497// Save jets
1498 Int_t n = fPythia->GetN();
1499
1500//
1501// Run Jet Finder
1502 fPythia->Pycell(njets);
1503 Int_t i;
1504 for (i = 0; i < njets; i++) {
1505 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
1506 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
1507 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
1508 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
1509
1510 jets[0][i] = px;
1511 jets[1][i] = py;
1512 jets[2][i] = pz;
1513 jets[3][i] = e;
1514 }
1515}
1516
1517
1518
1519void AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
1520{
1521//
1522// Calls the Pythia clustering algorithm to find jets in the current event
1523//
1524 Int_t n = fPythia->GetN();
1525 nJets = 0;
1526 nJetsTrig = 0;
1527 if (fJetReconstruction == kCluster) {
1528//
1529// Configure cluster algorithm
1530//
1531 fPythia->SetPARU(43, 2.);
1532 fPythia->SetMSTU(41, 1);
1533//
1534// Call cluster algorithm
1535//
1536 fPythia->Pyclus(nJets);
1537//
1538// Loading jets from common block
1539//
1540 } else {
592f8307 1541
8d2cd130 1542//
1543// Run Jet Finder
1544 fPythia->Pycell(nJets);
1545 }
1546
1547 Int_t i;
1548 for (i = 0; i < nJets; i++) {
1549 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
1550 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
1551 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
1552 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
1553 Float_t pt = TMath::Sqrt(px * px + py * py);
a920faf9 1554 Float_t phi = TMath::Pi() + TMath::ATan2(-py, -px);
8d2cd130 1555 Float_t theta = TMath::ATan2(pt,pz);
1556 Float_t et = e * TMath::Sin(theta);
1557 Float_t eta = -TMath::Log(TMath::Tan(theta / 2.));
8d2cd130 1558 if (
1559 eta > fEtaMinJet && eta < fEtaMaxJet &&
675eb105 1560 phi > fPhiMinJet && phi < fPhiMaxJet &&
8d2cd130 1561 et > fEtMinJet && et < fEtMaxJet
1562 )
1563 {
1564 jets[0][nJetsTrig] = px;
1565 jets[1][nJetsTrig] = py;
1566 jets[2][nJetsTrig] = pz;
1567 jets[3][nJetsTrig] = e;
1568 nJetsTrig++;
5fa4b20b 1569// printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
8d2cd130 1570 } else {
1571// printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1572 }
1573 }
1574}
1575
f913ec4f 1576void AliGenPythia::GetSubEventTime()
1577{
1578 // Calculates time of the next subevent
9d7108a7 1579 fEventTime = 0.;
1580 if (fEventsTime) {
1581 TArrayF &array = *fEventsTime;
1582 fEventTime = array[fCurSubEvent++];
1583 }
1584 // printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
1585 return;
f913ec4f 1586}
8d2cd130 1587
ec2c406e 1588Bool_t AliGenPythia::IsInEMCAL(Float_t phi, Float_t eta)
1589{
1590 // Is particle in EMCAL acceptance?
ec2c406e 1591 // phi in degrees, etamin=-etamax
9fd8e520 1592 if(phi > fEMCALMinPhi && phi < fEMCALMaxPhi &&
1593 eta < fEMCALEta )
ec2c406e 1594 return kTRUE;
1595 else
1596 return kFALSE;
1597}
1598
1599Bool_t AliGenPythia::IsInPHOS(Float_t phi, Float_t eta)
1600{
1601 // Is particle in PHOS acceptance?
1602 // Acceptance slightly larger considered.
1603 // phi in degrees, etamin=-etamax
9fd8e520 1604 if(phi > fPHOSMinPhi && phi < fPHOSMaxPhi &&
1605 eta < fPHOSEta )
ec2c406e 1606 return kTRUE;
1607 else
1608 return kFALSE;
1609}
1610
90a236ce 1611void AliGenPythia::RotatePhi(Int_t iphcand, Bool_t& okdd)
1612{
1613 //calculate the new position random between fPHOSMinPhi and fPHOSMaxPhi
1614 Double_t phiPHOSmin = TMath::Pi()*fPHOSMinPhi/180;
1615 Double_t phiPHOSmax = TMath::Pi()*fPHOSMaxPhi/180;
1616 Double_t phiPHOS = gRandom->Uniform(phiPHOSmin,phiPHOSmax);
1617
1618 //calculate deltaphi
8507138f 1619 TParticle* ph = (TParticle *) fParticles.At(iphcand);
90a236ce 1620 Double_t phphi = ph->Phi();
1621 Double_t deltaphi = phiPHOS - phphi;
1622
1623
1624
1625 //loop for all particles and produce the phi rotation
8507138f 1626 Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
90a236ce 1627 Double_t oldphi, newphi;
1628 Double_t newVx, newVy, R, Vz, time;
1629 Double_t newPx, newPy, pt, Pz, e;
1630 for(Int_t i=0; i< np; i++) {
8507138f 1631 TParticle* iparticle = (TParticle *) fParticles.At(i);
90a236ce 1632 oldphi = iparticle->Phi();
1633 newphi = oldphi + deltaphi;
1634 if(newphi < 0) newphi = 2*TMath::Pi() + newphi; // correct angle
1635 if(newphi > 2*TMath::Pi()) newphi = newphi - 2*TMath::Pi(); // correct angle
1636
1637 R = iparticle->R();
1638 newVx = R*TMath::Cos(newphi);
1639 newVy = R*TMath::Sin(newphi);
1640 Vz = iparticle->Vz(); // don't transform
1641 time = iparticle->T(); // don't transform
1642
1643 pt = iparticle->Pt();
1644 newPx = pt*TMath::Cos(newphi);
1645 newPy = pt*TMath::Sin(newphi);
1646 Pz = iparticle->Pz(); // don't transform
1647 e = iparticle->Energy(); // don't transform
1648
1649 // apply rotation
1650 iparticle->SetProductionVertex(newVx, newVy, Vz, time);
1651 iparticle->SetMomentum(newPx, newPy, Pz, e);
1652
1653 } //end particle loop
1654
1655 // now let's check that we put correctly the candidate photon in PHOS
1656 Float_t phi = ph->Phi()*180./TMath::Pi(); //Convert to degrees
1657 Float_t eta =TMath::Abs(ph->Eta());//in calos etamin=-etamax
1658 if(IsInPHOS(phi,eta))
1659 okdd = kTRUE;
1660}
ec2c406e 1661
1662
8d2cd130 1663#ifdef never
1664void AliGenPythia::Streamer(TBuffer &R__b)
1665{
1666 // Stream an object of class AliGenPythia.
1667
1668 if (R__b.IsReading()) {
1669 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1670 AliGenerator::Streamer(R__b);
1671 R__b >> (Int_t&)fProcess;
1672 R__b >> (Int_t&)fStrucFunc;
1673 R__b >> (Int_t&)fForceDecay;
1674 R__b >> fEnergyCMS;
1675 R__b >> fKineBias;
1676 R__b >> fTrials;
1677 fParentSelect.Streamer(R__b);
1678 fChildSelect.Streamer(R__b);
1679 R__b >> fXsection;
1680// (AliPythia::Instance())->Streamer(R__b);
1681 R__b >> fPtHardMin;
1682 R__b >> fPtHardMax;
1683// if (fDecayer) fDecayer->Streamer(R__b);
1684 } else {
1685 R__b.WriteVersion(AliGenPythia::IsA());
1686 AliGenerator::Streamer(R__b);
1687 R__b << (Int_t)fProcess;
1688 R__b << (Int_t)fStrucFunc;
1689 R__b << (Int_t)fForceDecay;
1690 R__b << fEnergyCMS;
1691 R__b << fKineBias;
1692 R__b << fTrials;
1693 fParentSelect.Streamer(R__b);
1694 fChildSelect.Streamer(R__b);
1695 R__b << fXsection;
1696// R__b << fPythia;
1697 R__b << fPtHardMin;
1698 R__b << fPtHardMax;
1699 // fDecayer->Streamer(R__b);
1700 }
1701}
1702#endif
1703
90d7b703 1704
589380c6 1705