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