Process kPyMbAtlasTuneMC09 added
[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 &&
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;
2f405d65 894 if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi) {
8507138f 895 TParticle* jet1 = (TParticle *) fParticles.At(6);
896 TParticle* jet2 = (TParticle *) fParticles.At(7);
234f6d32 897 if (!CheckTrigger(jet1, jet2)) {
898 delete [] pParent;
899 return 0;
900 }
8d2cd130 901 }
e0e89f40 902
9fd8e520 903 // Select jets with fragmentation photon or pi0 going to PHOS or EMCAL
904 if (fProcess == kPyJets && (fFragPhotonInCalo || fPi0InCalo) ) {
ec2c406e 905
906 Bool_t ok = kFALSE;
907
908 Int_t pdg = 0;
9fd8e520 909 if (fFragPhotonInCalo) pdg = 22 ; // Photon
82e0bdff 910 else if (fPi0InCalo) pdg = 111 ; // Pi0
ec2c406e 911
912 for (i=0; i< np; i++) {
8507138f 913 TParticle* iparticle = (TParticle *) fParticles.At(i);
ec2c406e 914 if(iparticle->GetStatusCode()==1 && iparticle->GetPdgCode()==pdg &&
9fd8e520 915 iparticle->Pt() > fFragPhotonOrPi0MinPt){
562cbbcf 916 Int_t imother = iparticle->GetFirstMother() - 1;
8507138f 917 TParticle* pmother = (TParticle *) fParticles.At(imother);
9fd8e520 918 if(pdg == 111 ||
82e0bdff 919 (pdg == 22 && pmother->GetStatusCode() != 11)) //No photon from hadron decay
9fd8e520 920 {
921 Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
82e0bdff 922 Float_t eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
9fd8e520 923 if((fCheckEMCAL && IsInEMCAL(phi,eta)) ||
924 (fCheckPHOS && IsInPHOS(phi,eta)) )
925 ok =kTRUE;
926 }
ec2c406e 927 }
928 }
929 if(!ok)
930 return 0;
931 }
9dfe63b3 932
933 // Select beauty jets with electron in EMCAL
934 if (fProcess == kPyBeautyJets && fEleInEMCAL) {
935
936 Bool_t ok = kFALSE;
937
938 Int_t pdg = 11; //electron
939
70574ff8 940 Float_t pt = 0.;
941 Float_t eta = 0.;
942 Float_t phi = 0.;
9dfe63b3 943 for (i=0; i< np; i++) {
944 TParticle* iparticle = (TParticle *) fParticles.At(i);
945 if(iparticle->GetStatusCode()==1 && TMath::Abs(iparticle->GetPdgCode())==pdg &&
946 iparticle->Pt() > fElectronMinPt){
947 pt = iparticle->Pt();
948 phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
949 eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
950 if(IsInEMCAL(phi,eta))
951 ok =kTRUE;
952 }
953 }
954 if(!ok)
955 return 0;
956 AliDebug(5,Form("Found an electron jet (pt,eta,phi) = (%f,%f,%f)",pt,eta,phi));
957 }
ec2c406e 958
700b9416 959 // Check for minimum multiplicity
960 if (fTriggerMultiplicity > 0) {
961 Int_t multiplicity = 0;
962 for (i = 0; i < np; i++) {
8507138f 963 TParticle * iparticle = (TParticle *) fParticles.At(i);
700b9416 964
965 Int_t statusCode = iparticle->GetStatusCode();
966
967 // Initial state particle
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 }
988
38112f3f 989 Printf("Triggered on event with multiplicity of %d >= %d", multiplicity, fTriggerMultiplicity);
700b9416 990 }
90a236ce 991
992 // Select events with a photon pt > min pt going to PHOS eta acceptance or exactly PHOS eta phi
993 if ((fProcess == kPyJets || fProcess == kPyDirectGamma) && fPhotonInCalo && (fCheckPHOSeta || fCheckPHOS)){
994
995 Bool_t okd = kFALSE;
996
997 Int_t pdg = 22;
998 Int_t iphcand = -1;
999 for (i=0; i< np; i++) {
8507138f 1000 TParticle* iparticle = (TParticle *) fParticles.At(i);
90a236ce 1001 Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
1002 Float_t eta =TMath::Abs(iparticle->Eta());//in calos etamin=-etamax
1003
1004 if(iparticle->GetStatusCode() == 1
1005 && iparticle->GetPdgCode() == pdg
1006 && iparticle->Pt() > fPhotonMinPt
04a26b0a 1007 && eta < fPHOSEta){
90a236ce 1008
1009 // first check if the photon is in PHOS phi
1010 if(IsInPHOS(phi,eta)){
1011 okd = kTRUE;
1012 break;
1013 }
1014 if(fCheckPHOSeta) iphcand = i; // candiate photon to rotate in phi
1015
1016 }
1017 }
1018
1019 if(!okd && iphcand != -1) // execute rotation in phi
1020 RotatePhi(iphcand,okd);
1021
1022 if(!okd)
1023 return 0;
1024 }
1025
7ce1321b 1026 if (fTriggerParticle) {
1027 Bool_t triggered = kFALSE;
1028 for (i = 0; i < np; i++) {
8507138f 1029 TParticle * iparticle = (TParticle *) fParticles.At(i);
7ce1321b 1030 kf = CheckPDGCode(iparticle->GetPdgCode());
38db0ee6 1031 if (kf != fTriggerParticle) continue;
7ce1321b 1032 if (iparticle->Pt() == 0.) continue;
1033 if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
1034 triggered = kTRUE;
1035 break;
1036 }
234f6d32 1037 if (!triggered) {
1038 delete [] pParent;
1039 return 0;
1040 }
7ce1321b 1041 }
e0e89f40 1042
1043
1044 // Check if there is a ccbar or bbbar pair with at least one of the two
1045 // in fYMin < y < fYMax
2f405d65 1046
9dfe63b3 1047 if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi || fProcess == kPyBeautyJets) {
9ccc3d0e 1048 TParticle *partCheck;
1049 TParticle *mother;
e0e89f40 1050 Bool_t theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
9ccc3d0e 1051 Bool_t theChild=kFALSE;
1052 Float_t y;
1053 Int_t pdg,mpdg,mpdgUpperFamily;
e0e89f40 1054 for(i=0; i<np; i++) {
9ccc3d0e 1055 partCheck = (TParticle*)fParticles.At(i);
1056 pdg = partCheck->GetPdgCode();
1057 if(TMath::Abs(pdg) == fFlavorSelect) { // quark
1058 if(pdg>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
1059 y = 0.5*TMath::Log((partCheck->Energy()+partCheck->Pz()+1.e-13)/
1060 (partCheck->Energy()-partCheck->Pz()+1.e-13));
1061 if(y>fYMin && y<fYMax) inYcut=kTRUE;
1062 }
1063 if(fCutOnChild && TMath::Abs(pdg) == fPdgCodeParticleforAcceptanceCut) {
1064 Int_t mi = partCheck->GetFirstMother() - 1;
1065 if(mi<0) continue;
1066 mother = (TParticle*)fParticles.At(mi);
1067 mpdg=TMath::Abs(mother->GetPdgCode());
75fde694 1068 mpdgUpperFamily=(mpdg>1000 ? mpdg+1000 : mpdg+100); // keep e from c from b
9ccc3d0e 1069 if ( ParentSelected(mpdg) ||
1070 (fFlavorSelect==5 && ParentSelected(mpdgUpperFamily))) {
1071 if (KinematicSelection(partCheck,1)) {
1072 theChild=kTRUE;
1073 }
1074 }
1075 }
e0e89f40 1076 }
9ccc3d0e 1077 if (!theQ || !theQbar || !inYcut) { // one of the c/b conditions not satisfied
234f6d32 1078 delete[] pParent;
e0e89f40 1079 return 0;
1080 }
9ccc3d0e 1081 if (fCutOnChild && !theChild) { // one of the child conditions not satisfied
1082 delete[] pParent;
1083 return 0;
1084 }
1085
e0e89f40 1086 }
aea21c57 1087
0f6ee828 1088 //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
f529e69b 1089 if ( (fProcess == kPyW ||
1090 fProcess == kPyZ ||
1091 fProcess == kPyMbDefault ||
1092 fProcess == kPyMb ||
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