]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PYTHIA6/AliGenPythia.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fixing coding violations. Method AliITSQADataMakerRec::AreEqual defined static
[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;
6c43eccb 587 fPythia->Qpygin0();
cd07c39b 588 }
32d6ef7d 589//
5fa4b20b 590// Either produce new event or read partons from file
591//
592 if (!fReadFromFile) {
beac474c 593 if (!fNewMIS) {
594 fPythia->Pyevnt();
595 } else {
596 fPythia->Pyevnw();
597 }
5fa4b20b 598 fNpartons = fPythia->GetN();
599 } else {
33c3c91a 600 printf("Loading Event %d\n",AliRunLoader::Instance()->GetEventNumber());
601 fRL->GetEvent(AliRunLoader::Instance()->GetEventNumber());
5fa4b20b 602 fPythia->SetN(0);
603 LoadEvent(fRL->Stack(), 0 , 1);
604 fPythia->Pyedit(21);
605 }
606
32d6ef7d 607//
608// Run quenching routine
609//
5fa4b20b 610 if (fQuench == 1) {
611 fPythia->Quench();
612 } else if (fQuench == 2){
613 fPythia->Pyquen(208., 0, 0.);
b976f7d8 614 } else if (fQuench == 3) {
615 // Quenching is via multiplicative correction of the splittings
5fa4b20b 616 }
b976f7d8 617
32d6ef7d 618//
5fa4b20b 619// Switch hadronisation on
32d6ef7d 620//
20e47f08 621 if (fHadronisation) {
622 fPythia->SetMSTJ(1, 1);
5fa4b20b 623//
624// .. and perform hadronisation
aea21c57 625// printf("Calling hadronisation %d\n", fPythia->GetN());
20e47f08 626 fPythia->Pyexec();
627 }
8d2cd130 628 fTrials++;
8507138f 629 fPythia->ImportParticles(&fParticles,"All");
df275cfa 630 if (TMath::Abs(fDyBoost) > 1.e-4) Boost();
8d2cd130 631//
632//
633//
634 Int_t i;
635
dbd64db6 636 fNprimaries = 0;
8507138f 637 Int_t np = fParticles.GetEntriesFast();
5fa4b20b 638
7c21f297 639 if (np == 0) continue;
8d2cd130 640//
2590ccf9 641
8d2cd130 642//
643 Int_t* pParent = new Int_t[np];
644 Int_t* pSelected = new Int_t[np];
645 Int_t* trackIt = new Int_t[np];
5fa4b20b 646 for (i = 0; i < np; i++) {
8d2cd130 647 pParent[i] = -1;
648 pSelected[i] = 0;
649 trackIt[i] = 0;
650 }
651
652 Int_t nc = 0; // Total n. of selected particles
653 Int_t nParents = 0; // Selected parents
654 Int_t nTkbles = 0; // Trackable particles
f529e69b 655 if (fProcess != kPyMbDefault &&
656 fProcess != kPyMb &&
6d2ec66d 657 fProcess != kPyMbAtlasTuneMC09 &&
04081a91 658 fProcess != kPyMbWithDirectPhoton &&
f529e69b 659 fProcess != kPyJets &&
8d2cd130 660 fProcess != kPyDirectGamma &&
589380c6 661 fProcess != kPyMbNonDiffr &&
d7de4a67 662 fProcess != kPyMbMSEL1 &&
f529e69b 663 fProcess != kPyW &&
664 fProcess != kPyZ &&
665 fProcess != kPyCharmppMNRwmi &&
9dfe63b3 666 fProcess != kPyBeautyppMNRwmi &&
667 fProcess != kPyBeautyJets) {
8d2cd130 668
5fa4b20b 669 for (i = 0; i < np; i++) {
8507138f 670 TParticle* iparticle = (TParticle *) fParticles.At(i);
8d2cd130 671 Int_t ks = iparticle->GetStatusCode();
672 kf = CheckPDGCode(iparticle->GetPdgCode());
673// No initial state partons
674 if (ks==21) continue;
675//
676// Heavy Flavor Selection
677//
678 // quark ?
679 kf = TMath::Abs(kf);
680 Int_t kfl = kf;
9ff6c04c 681 // Resonance
f913ec4f 682
9ff6c04c 683 if (kfl > 100000) kfl %= 100000;
183a5ca9 684 if (kfl > 10000) kfl %= 10000;
8d2cd130 685 // meson ?
686 if (kfl > 10) kfl/=100;
687 // baryon
688 if (kfl > 10) kfl/=10;
8d2cd130 689 Int_t ipa = iparticle->GetFirstMother()-1;
690 Int_t kfMo = 0;
f913ec4f 691//
692// Establish mother daughter relation between heavy quarks and mesons
693//
694 if (kf >= fFlavorSelect && kf <= 6) {
695 Int_t idau = iparticle->GetFirstDaughter() - 1;
696 if (idau > -1) {
8507138f 697 TParticle* daughter = (TParticle *) fParticles.At(idau);
f913ec4f 698 Int_t pdgD = daughter->GetPdgCode();
699 if (pdgD == 91 || pdgD == 92) {
700 Int_t jmin = daughter->GetFirstDaughter() - 1;
701 Int_t jmax = daughter->GetLastDaughter() - 1;
2ab330c9 702 for (Int_t jp = jmin; jp <= jmax; jp++)
703 ((TParticle *) fParticles.At(jp))->SetFirstMother(i+1);
f913ec4f 704 } // is string or cluster
705 } // has daughter
706 } // heavy quark
8d2cd130 707
f913ec4f 708
8d2cd130 709 if (ipa > -1) {
8507138f 710 TParticle * mother = (TParticle *) fParticles.At(ipa);
8d2cd130 711 kfMo = TMath::Abs(mother->GetPdgCode());
712 }
f913ec4f 713
8d2cd130 714 // What to keep in Stack?
715 Bool_t flavorOK = kFALSE;
716 Bool_t selectOK = kFALSE;
717 if (fFeedDownOpt) {
32d6ef7d 718 if (kfl >= fFlavorSelect) flavorOK = kTRUE;
8d2cd130 719 } else {
32d6ef7d 720 if (kfl > fFlavorSelect) {
721 nc = -1;
722 break;
723 }
724 if (kfl == fFlavorSelect) flavorOK = kTRUE;
8d2cd130 725 }
726 switch (fStackFillOpt) {
727 case kFlavorSelection:
32d6ef7d 728 selectOK = kTRUE;
729 break;
8d2cd130 730 case kParentSelection:
32d6ef7d 731 if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
732 break;
8d2cd130 733 }
734 if (flavorOK && selectOK) {
735//
736// Heavy flavor hadron or quark
737//
738// Kinematic seletion on final state heavy flavor mesons
739 if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
740 {
9ff6c04c 741 continue;
8d2cd130 742 }
743 pSelected[i] = 1;
744 if (ParentSelected(kf)) ++nParents; // Update parent count
745// printf("\n particle (HF) %d %d %d", i, pSelected[i], kf);
746 } else {
747// Kinematic seletion on decay products
748 if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf)
9ff6c04c 749 && !KinematicSelection(iparticle, 1))
8d2cd130 750 {
9ff6c04c 751 continue;
8d2cd130 752 }
753//
754// Decay products
755// Select if mother was selected and is not tracked
756
757 if (pSelected[ipa] &&
758 !trackIt[ipa] && // mother will be tracked ?
759 kfMo != 5 && // mother is b-quark, don't store fragments
760 kfMo != 4 && // mother is c-quark, don't store fragments
761 kf != 92) // don't store string
762 {
763//
764// Semi-stable or de-selected: diselect decay products:
765//
766//
767 if (pSelected[i] == -1 || fDecayer->GetLifetime(kf) > fMaxLifeTime)
768 {
769 Int_t ipF = iparticle->GetFirstDaughter();
770 Int_t ipL = iparticle->GetLastDaughter();
771 if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
772 }
773// printf("\n particle (decay) %d %d %d", i, pSelected[i], kf);
774 pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
775 }
776 }
777 if (pSelected[i] == -1) pSelected[i] = 0;
778 if (!pSelected[i]) continue;
779 // Count quarks only if you did not include fragmentation
780 if (fFragmentation && kf <= 10) continue;
9ff6c04c 781
8d2cd130 782 nc++;
783// Decision on tracking
784 trackIt[i] = 0;
785//
786// Track final state particle
787 if (ks == 1) trackIt[i] = 1;
788// Track semi-stable particles
d25cfd65 789 if ((ks == 1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime)) trackIt[i] = 1;
8d2cd130 790// Track particles selected by process if undecayed.
791 if (fForceDecay == kNoDecay) {
792 if (ParentSelected(kf)) trackIt[i] = 1;
793 } else {
794 if (ParentSelected(kf)) trackIt[i] = 0;
795 }
796 if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
797//
798//
799
800 } // particle selection loop
801 if (nc > 0) {
802 for (i = 0; i<np; i++) {
803 if (!pSelected[i]) continue;
8507138f 804 TParticle * iparticle = (TParticle *) fParticles.At(i);
8d2cd130 805 kf = CheckPDGCode(iparticle->GetPdgCode());
806 Int_t ks = iparticle->GetStatusCode();
807 p[0] = iparticle->Px();
808 p[1] = iparticle->Py();
809 p[2] = iparticle->Pz();
a920faf9 810 p[3] = iparticle->Energy();
811
2590ccf9 812 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
813 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
814 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
815
8d2cd130 816 Float_t tof = kconv*iparticle->T();
817 Int_t ipa = iparticle->GetFirstMother()-1;
818 Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
a920faf9 819
820 PushTrack(fTrackIt*trackIt[i], iparent, kf,
821 p[0], p[1], p[2], p[3],
822 origin[0], origin[1], origin[2], tof,
823 polar[0], polar[1], polar[2],
824 kPPrimary, nt, 1., ks);
8d2cd130 825 pParent[i] = nt;
dbd64db6 826 KeepTrack(nt);
827 fNprimaries++;
642f15cf 828 } // PushTrack loop
8d2cd130 829 }
830 } else {
831 nc = GenerateMB();
832 } // mb ?
f913ec4f 833
834 GetSubEventTime();
8d2cd130 835
234f6d32 836 delete[] pParent;
837 delete[] pSelected;
838 delete[] trackIt;
8d2cd130 839
840 if (nc > 0) {
841 switch (fCountMode) {
842 case kCountAll:
843 // printf(" Count all \n");
844 jev += nc;
845 break;
846 case kCountParents:
847 // printf(" Count parents \n");
848 jev += nParents;
849 break;
850 case kCountTrackables:
851 // printf(" Count trackable \n");
852 jev += nTkbles;
853 break;
854 }
855 if (jev >= fNpart || fNpart == -1) {
856 fKineBias=Float_t(fNpart)/Float_t(fTrials);
e0e89f40 857
8d2cd130 858 fQ += fPythia->GetVINT(51);
859 fX1 += fPythia->GetVINT(41);
860 fX2 += fPythia->GetVINT(42);
861 fTrialsRun += fTrials;
862 fNev++;
863 MakeHeader();
864 break;
865 }
866 }
867 } // event loop
868 SetHighWaterMark(nt);
869// adjust weight due to kinematic selection
b88f5cea 870// AdjustWeights();
8d2cd130 871// get cross-section
872 fXsection=fPythia->GetPARI(1);
873}
874
875Int_t AliGenPythia::GenerateMB()
876{
877//
878// Min Bias selection and other global selections
879//
880 Int_t i, kf, nt, iparent;
881 Int_t nc = 0;
bf950da8 882 Float_t p[4];
8d2cd130 883 Float_t polar[3] = {0,0,0};
884 Float_t origin[3] = {0,0,0};
885// converts from mm/c to s
886 const Float_t kconv=0.001/2.999792458e8;
887
e0e89f40 888
889
8507138f 890 Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
aea21c57 891
5fa4b20b 892
e0e89f40 893
8d2cd130 894 Int_t* pParent = new Int_t[np];
895 for (i=0; i< np; i++) pParent[i] = -1;
2f405d65 896 if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi) {
8507138f 897 TParticle* jet1 = (TParticle *) fParticles.At(6);
898 TParticle* jet2 = (TParticle *) fParticles.At(7);
234f6d32 899 if (!CheckTrigger(jet1, jet2)) {
900 delete [] pParent;
901 return 0;
902 }
8d2cd130 903 }
e0e89f40 904
9fd8e520 905 // Select jets with fragmentation photon or pi0 going to PHOS or EMCAL
906 if (fProcess == kPyJets && (fFragPhotonInCalo || fPi0InCalo) ) {
ec2c406e 907
908 Bool_t ok = kFALSE;
909
910 Int_t pdg = 0;
9fd8e520 911 if (fFragPhotonInCalo) pdg = 22 ; // Photon
6d2ec66d 912 else if (fPi0InCalo) pdg = 111 ; // Pi0
ec2c406e 913
914 for (i=0; i< np; i++) {
8507138f 915 TParticle* iparticle = (TParticle *) fParticles.At(i);
ec2c406e 916 if(iparticle->GetStatusCode()==1 && iparticle->GetPdgCode()==pdg &&
9fd8e520 917 iparticle->Pt() > fFragPhotonOrPi0MinPt){
562cbbcf 918 Int_t imother = iparticle->GetFirstMother() - 1;
8507138f 919 TParticle* pmother = (TParticle *) fParticles.At(imother);
9fd8e520 920 if(pdg == 111 ||
82e0bdff 921 (pdg == 22 && pmother->GetStatusCode() != 11)) //No photon from hadron decay
9fd8e520 922 {
923 Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
82e0bdff 924 Float_t eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
9fd8e520 925 if((fCheckEMCAL && IsInEMCAL(phi,eta)) ||
926 (fCheckPHOS && IsInPHOS(phi,eta)) )
927 ok =kTRUE;
928 }
ec2c406e 929 }
930 }
931 if(!ok)
932 return 0;
933 }
9dfe63b3 934
935 // Select beauty jets with electron in EMCAL
936 if (fProcess == kPyBeautyJets && fEleInEMCAL) {
937
938 Bool_t ok = kFALSE;
939
940 Int_t pdg = 11; //electron
941
70574ff8 942 Float_t pt = 0.;
943 Float_t eta = 0.;
944 Float_t phi = 0.;
9dfe63b3 945 for (i=0; i< np; i++) {
946 TParticle* iparticle = (TParticle *) fParticles.At(i);
947 if(iparticle->GetStatusCode()==1 && TMath::Abs(iparticle->GetPdgCode())==pdg &&
948 iparticle->Pt() > fElectronMinPt){
949 pt = iparticle->Pt();
950 phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
951 eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
952 if(IsInEMCAL(phi,eta))
953 ok =kTRUE;
954 }
955 }
956 if(!ok)
957 return 0;
958 AliDebug(5,Form("Found an electron jet (pt,eta,phi) = (%f,%f,%f)",pt,eta,phi));
959 }
700b9416 960 // Check for minimum multiplicity
961 if (fTriggerMultiplicity > 0) {
962 Int_t multiplicity = 0;
963 for (i = 0; i < np; i++) {
8507138f 964 TParticle * iparticle = (TParticle *) fParticles.At(i);
700b9416 965
966 Int_t statusCode = iparticle->GetStatusCode();
967
968 // Initial state particle
e296848e 969 if (statusCode != 1)
700b9416 970 continue;
38112f3f 971 // eta cut
700b9416 972 if (fTriggerMultiplicityEta > 0 && TMath::Abs(iparticle->Eta()) > fTriggerMultiplicityEta)
973 continue;
38112f3f 974 // pt cut
975 if (iparticle->Pt() < fTriggerMultiplicityPtMin)
976 continue;
977
700b9416 978 TParticlePDG* pdgPart = iparticle->GetPDG();
979 if (pdgPart && pdgPart->Charge() == 0)
980 continue;
981
982 ++multiplicity;
983 }
984
985 if (multiplicity < fTriggerMultiplicity) {
986 delete [] pParent;
987 return 0;
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 ||
6d2ec66d 1093 fProcess == kPyMbAtlasTuneMC09 ||
04081a91 1094 fProcess == kPyMbWithDirectPhoton ||
f529e69b 1095 fProcess == kPyMbNonDiffr)
0f6ee828 1096 && (fCutOnChild == 1) ) {
1097 if ( !CheckKinematicsOnChild() ) {
234f6d32 1098 delete[] pParent;
0f6ee828 1099 return 0;
1100 }
aea21c57 1101 }
1102
1103
f913ec4f 1104 for (i = 0; i < np; i++) {
8d2cd130 1105 Int_t trackIt = 0;
8507138f 1106 TParticle * iparticle = (TParticle *) fParticles.At(i);
8d2cd130 1107 kf = CheckPDGCode(iparticle->GetPdgCode());
1108 Int_t ks = iparticle->GetStatusCode();
1109 Int_t km = iparticle->GetFirstMother();
1110 if ((ks == 1 && kf!=0 && KinematicSelection(iparticle, 0)) ||
1111 (ks != 1) ||
9dfe63b3 1112 ((fProcess == kPyJets || fProcess == kPyBeautyJets) && ks == 21 && km == 0 && i>1)) {
8d2cd130 1113 nc++;
1114 if (ks == 1) trackIt = 1;
1115 Int_t ipa = iparticle->GetFirstMother()-1;
1116
1117 iparent = (ipa > -1) ? pParent[ipa] : -1;
1118
1119//
1120// store track information
1121 p[0] = iparticle->Px();
1122 p[1] = iparticle->Py();
1123 p[2] = iparticle->Pz();
a920faf9 1124 p[3] = iparticle->Energy();
1406f599 1125
a920faf9 1126
2590ccf9 1127 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
1128 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
1129 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
1130
f913ec4f 1131 Float_t tof = fEventTime + kconv * iparticle->T();
a920faf9 1132
1133 PushTrack(fTrackIt*trackIt, iparent, kf,
1134 p[0], p[1], p[2], p[3],
1135 origin[0], origin[1], origin[2], tof,
1136 polar[0], polar[1], polar[2],
1137 kPPrimary, nt, 1., ks);
dbd64db6 1138 fNprimaries++;
8d2cd130 1139 KeepTrack(nt);
1140 pParent[i] = nt;
f913ec4f 1141 SetHighWaterMark(nt);
1142
8d2cd130 1143 } // select particle
1144 } // particle loop
1145
234f6d32 1146 delete[] pParent;
e0e89f40 1147
f913ec4f 1148 return 1;
8d2cd130 1149}
1150
1151
1152void AliGenPythia::FinishRun()
1153{
1154// Print x-section summary
1155 fPythia->Pystat(1);
2779fc64 1156
1157 if (fNev > 0.) {
1158 fQ /= fNev;
1159 fX1 /= fNev;
1160 fX2 /= fNev;
1161 }
1162
8d2cd130 1163 printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
1164 printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
8d2cd130 1165}
1166
7184e472 1167void AliGenPythia::AdjustWeights() const
8d2cd130 1168{
1169// Adjust the weights after generation of all events
1170//
e2bddf81 1171 if (gAlice) {
1172 TParticle *part;
1173 Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
1174 for (Int_t i=0; i<ntrack; i++) {
1175 part= gAlice->GetMCApp()->Particle(i);
1176 part->SetWeight(part->GetWeight()*fKineBias);
1177 }
8d2cd130 1178 }
1179}
1180
20e47f08 1181void AliGenPythia::SetNuclei(Int_t a1, Int_t a2, Int_t pdfset)
8d2cd130 1182{
1183// Treat protons as inside nuclei with mass numbers a1 and a2
1d568bc2 1184
1a626d4e 1185 fAProjectile = a1;
1186 fATarget = a2;
20e47f08 1187 fNucPdf = pdfset; // 0 EKS98 1 EPS08
1d568bc2 1188 fSetNuclei = kTRUE;
8d2cd130 1189}
1190
1191
1192void AliGenPythia::MakeHeader()
1193{
7184e472 1194//
1195// Make header for the simulated event
1196//
183a5ca9 1197 if (gAlice) {
1198 if (gAlice->GetEvNumber()>=fDebugEventFirst &&
f913ec4f 1199 gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(2);
183a5ca9 1200 }
1201
8d2cd130 1202// Builds the event header, to be called after each event
e5c87a3d 1203 if (fHeader) delete fHeader;
1204 fHeader = new AliGenPythiaEventHeader("Pythia");
8d2cd130 1205//
1206// Event type
e5c87a3d 1207 ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
8d2cd130 1208//
1209// Number of trials
e5c87a3d 1210 ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
8d2cd130 1211//
1212// Event Vertex
d25cfd65 1213 fHeader->SetPrimaryVertex(fVertex);
d07f0af2 1214 fHeader->SetInteractionTime(fEventTime);
dbd64db6 1215//
1216// Number of primaries
1217 fHeader->SetNProduced(fNprimaries);
8d2cd130 1218//
1219// Jets that have triggered
f913ec4f 1220
9dfe63b3 1221 //Need to store jets for b-jet studies too!
2f405d65 1222 if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi)
8d2cd130 1223 {
1224 Int_t ntrig, njet;
1225 Float_t jets[4][10];
1226 GetJets(njet, ntrig, jets);
9ff6c04c 1227
8d2cd130 1228
1229 for (Int_t i = 0; i < ntrig; i++) {
e5c87a3d 1230 ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
8d2cd130 1231 jets[3][i]);
1232 }
1233 }
5fa4b20b 1234//
1235// Copy relevant information from external header, if present.
1236//
1237 Float_t uqJet[4];
1238
1239 if (fRL) {
1240 AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
1241 for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
1242 {
1243 printf("Adding Jet %d %d \n", i, exHeader->NTriggerJets());
1244
1245
1246 exHeader->TriggerJet(i, uqJet);
1247 ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
1248 }
1249 }
1250//
1251// Store quenching parameters
1252//
1253 if (fQuench){
1254 Double_t z[4];
1255 Double_t xp, yp;
7c21f297 1256 if (fQuench == 1) {
1257 // Pythia::Quench()
1258 fPythia->GetQuenchingParameters(xp, yp, z);
1bab4b79 1259 } else if (fQuench == 2){
7c21f297 1260 // Pyquen
1261 Double_t r1 = PARIMP.rb1;
1262 Double_t r2 = PARIMP.rb2;
1263 Double_t b = PARIMP.b1;
1264 Double_t r = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
1265 Double_t phi = PARIMP.psib1;
1266 xp = r * TMath::Cos(phi);
1267 yp = r * TMath::Sin(phi);
1268
1bab4b79 1269 } else if (fQuench == 4) {
1270 // QPythia
5831e75f 1271 Double_t xy[2];
e9719084 1272 Double_t i0i1[2];
1273 AliFastGlauber::Instance()->GetSavedXY(xy);
1274 AliFastGlauber::Instance()->GetSavedI0I1(i0i1);
5831e75f 1275 xp = xy[0];
1276 yp = xy[1];
e6fe9b82 1277 ((AliGenPythiaEventHeader*) fHeader)->SetImpactParameter(fImpact);
7c21f297 1278 }
1bab4b79 1279
7c21f297 1280 ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
1281 ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
1bab4b79 1282 }
beac474c 1283//
1284// Store pt^hard
1285 ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
5fa4b20b 1286//
cf57b268 1287// Pass header
5fa4b20b 1288//
cf57b268 1289 AddHeader(fHeader);
4c4eac97 1290 fHeader = 0x0;
8d2cd130 1291}
cf57b268 1292
8d2cd130 1293Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
1294{
1295// Check the kinematic trigger condition
1296//
1297 Double_t eta[2];
1298 eta[0] = jet1->Eta();
1299 eta[1] = jet2->Eta();
1300 Double_t phi[2];
1301 phi[0] = jet1->Phi();
1302 phi[1] = jet2->Phi();
1303 Int_t pdg[2];
1304 pdg[0] = jet1->GetPdgCode();
1305 pdg[1] = jet2->GetPdgCode();
1306 Bool_t triggered = kFALSE;
1307
2f405d65 1308 if (fProcess == kPyJets || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi) {
8d2cd130 1309 Int_t njets = 0;
1310 Int_t ntrig = 0;
1311 Float_t jets[4][10];
1312//
1313// Use Pythia clustering on parton level to determine jet axis
1314//
1315 GetJets(njets, ntrig, jets);
1316
76d6ba9a 1317 if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
8d2cd130 1318//
1319 } else {
1320 Int_t ij = 0;
1321 Int_t ig = 1;
1322 if (pdg[0] == kGamma) {
1323 ij = 1;
1324 ig = 0;
1325 }
1326 //Check eta range first...
1327 if ((eta[ij] < fEtaMaxJet && eta[ij] > fEtaMinJet) &&
1328 (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
1329 {
1330 //Eta is okay, now check phi range
1331 if ((phi[ij] < fPhiMaxJet && phi[ij] > fPhiMinJet) &&
1332 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
1333 {
1334 triggered = kTRUE;
1335 }
1336 }
1337 }
1338 return triggered;
1339}
aea21c57 1340
1341
aea21c57 1342
7184e472 1343Bool_t AliGenPythia::CheckKinematicsOnChild(){
1344//
1345//Checking Kinematics on Child (status code 1, particle code ?, kin cuts
1346//
aea21c57 1347 Bool_t checking = kFALSE;
1348 Int_t j, kcode, ks, km;
1349 Int_t nPartAcc = 0; //number of particles in the acceptance range
1350 Int_t numberOfAcceptedParticles = 1;
1351 if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
8507138f 1352 Int_t npart = fParticles.GetEntriesFast();
aea21c57 1353
0f6ee828 1354 for (j = 0; j<npart; j++) {
8507138f 1355 TParticle * jparticle = (TParticle *) fParticles.At(j);
aea21c57 1356 kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
1357 ks = jparticle->GetStatusCode();
1358 km = jparticle->GetFirstMother();
1359
1360 if( (ks == 1) && (kcode == fPdgCodeParticleforAcceptanceCut) && (KinematicSelection(jparticle,1)) ){
1361 nPartAcc++;
1362 }
0f6ee828 1363 if( numberOfAcceptedParticles <= nPartAcc){
1364 checking = kTRUE;
1365 break;
1366 }
aea21c57 1367 }
0f6ee828 1368
aea21c57 1369 return checking;
aea21c57 1370}
1371
5fa4b20b 1372void AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
8d2cd130 1373{
1058d9df 1374 //
1375 // Load event into Pythia Common Block
1376 //
1377
1378 Int_t npart = stack -> GetNprimary();
1379 Int_t n0 = 0;
1380
1381 if (!flag) {
1382 (fPythia->GetPyjets())->N = npart;
1383 } else {
1384 n0 = (fPythia->GetPyjets())->N;
1385 (fPythia->GetPyjets())->N = n0 + npart;
1386 }
1387
1388
1389 for (Int_t part = 0; part < npart; part++) {
1390 TParticle *mPart = stack->Particle(part);
32d6ef7d 1391
1058d9df 1392 Int_t kf = mPart->GetPdgCode();
1393 Int_t ks = mPart->GetStatusCode();
1394 Int_t idf = mPart->GetFirstDaughter();
1395 Int_t idl = mPart->GetLastDaughter();
1396
1397 if (reHadr) {
1398 if (ks == 11 || ks == 12) {
1399 ks -= 10;
1400 idf = -1;
1401 idl = -1;
1402 }
32d6ef7d 1403 }
1404
1058d9df 1405 Float_t px = mPart->Px();
1406 Float_t py = mPart->Py();
1407 Float_t pz = mPart->Pz();
1408 Float_t e = mPart->Energy();
1409 Float_t m = mPart->GetCalcMass();
32d6ef7d 1410
1058d9df 1411
1412 (fPythia->GetPyjets())->P[0][part+n0] = px;
1413 (fPythia->GetPyjets())->P[1][part+n0] = py;
1414 (fPythia->GetPyjets())->P[2][part+n0] = pz;
1415 (fPythia->GetPyjets())->P[3][part+n0] = e;
1416 (fPythia->GetPyjets())->P[4][part+n0] = m;
1417
1418 (fPythia->GetPyjets())->K[1][part+n0] = kf;
1419 (fPythia->GetPyjets())->K[0][part+n0] = ks;
1420 (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1421 (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1422 (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1423 }
1424}
1425
1426void AliGenPythia::LoadEvent(TObjArray* stack, Int_t flag, Int_t reHadr)
1427{
1428 //
1429 // Load event into Pythia Common Block
1430 //
1431
1432 Int_t npart = stack -> GetEntries();
1433 Int_t n0 = 0;
1434
1435 if (!flag) {
1436 (fPythia->GetPyjets())->N = npart;
1437 } else {
1438 n0 = (fPythia->GetPyjets())->N;
1439 (fPythia->GetPyjets())->N = n0 + npart;
1440 }
1441
1442
1443 for (Int_t part = 0; part < npart; part++) {
1444 TParticle *mPart = dynamic_cast<TParticle *>(stack->At(part));
1445 Int_t kf = mPart->GetPdgCode();
1446 Int_t ks = mPart->GetStatusCode();
1447 Int_t idf = mPart->GetFirstDaughter();
1448 Int_t idl = mPart->GetLastDaughter();
1449
1450 if (reHadr) {
5fa4b20b 1451 if (ks == 11 || ks == 12) {
1058d9df 1452 ks -= 10;
1453 idf = -1;
1454 idl = -1;
5fa4b20b 1455 }
8d2cd130 1456 }
1058d9df 1457
1458 Float_t px = mPart->Px();
1459 Float_t py = mPart->Py();
1460 Float_t pz = mPart->Pz();
1461 Float_t e = mPart->Energy();
1462 Float_t m = mPart->GetCalcMass();
1463
1464
1465 (fPythia->GetPyjets())->P[0][part+n0] = px;
1466 (fPythia->GetPyjets())->P[1][part+n0] = py;
1467 (fPythia->GetPyjets())->P[2][part+n0] = pz;
1468 (fPythia->GetPyjets())->P[3][part+n0] = e;
1469 (fPythia->GetPyjets())->P[4][part+n0] = m;
1470
1471 (fPythia->GetPyjets())->K[1][part+n0] = kf;
1472 (fPythia->GetPyjets())->K[0][part+n0] = ks;
1473 (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1474 (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1475 (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1476 }
8d2cd130 1477}
1478
5fa4b20b 1479
014a9521 1480void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
8d2cd130 1481{
1482//
1483// Calls the Pythia jet finding algorithm to find jets in the current event
1484//
1485//
8d2cd130 1486//
1487// Save jets
1488 Int_t n = fPythia->GetN();
1489
1490//
1491// Run Jet Finder
1492 fPythia->Pycell(njets);
1493 Int_t i;
1494 for (i = 0; i < njets; i++) {
1495 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
1496 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
1497 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
1498 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
1499
1500 jets[0][i] = px;
1501 jets[1][i] = py;
1502 jets[2][i] = pz;
1503 jets[3][i] = e;
1504 }
1505}
1506
1507
1508
1509void AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
1510{
1511//
1512// Calls the Pythia clustering algorithm to find jets in the current event
1513//
1514 Int_t n = fPythia->GetN();
1515 nJets = 0;
1516 nJetsTrig = 0;
1517 if (fJetReconstruction == kCluster) {
1518//
1519// Configure cluster algorithm
1520//
1521 fPythia->SetPARU(43, 2.);
1522 fPythia->SetMSTU(41, 1);
1523//
1524// Call cluster algorithm
1525//
1526 fPythia->Pyclus(nJets);
1527//
1528// Loading jets from common block
1529//
1530 } else {
592f8307 1531
8d2cd130 1532//
1533// Run Jet Finder
1534 fPythia->Pycell(nJets);
1535 }
1536
1537 Int_t i;
1538 for (i = 0; i < nJets; i++) {
1539 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
1540 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
1541 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
1542 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
1543 Float_t pt = TMath::Sqrt(px * px + py * py);
a920faf9 1544 Float_t phi = TMath::Pi() + TMath::ATan2(-py, -px);
8d2cd130 1545 Float_t theta = TMath::ATan2(pt,pz);
1546 Float_t et = e * TMath::Sin(theta);
1547 Float_t eta = -TMath::Log(TMath::Tan(theta / 2.));
8d2cd130 1548 if (
1549 eta > fEtaMinJet && eta < fEtaMaxJet &&
675eb105 1550 phi > fPhiMinJet && phi < fPhiMaxJet &&
8d2cd130 1551 et > fEtMinJet && et < fEtMaxJet
1552 )
1553 {
1554 jets[0][nJetsTrig] = px;
1555 jets[1][nJetsTrig] = py;
1556 jets[2][nJetsTrig] = pz;
1557 jets[3][nJetsTrig] = e;
1558 nJetsTrig++;
5fa4b20b 1559// printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
8d2cd130 1560 } else {
1561// printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1562 }
1563 }
1564}
1565
f913ec4f 1566void AliGenPythia::GetSubEventTime()
1567{
1568 // Calculates time of the next subevent
9d7108a7 1569 fEventTime = 0.;
1570 if (fEventsTime) {
1571 TArrayF &array = *fEventsTime;
1572 fEventTime = array[fCurSubEvent++];
1573 }
1574 // printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
1575 return;
f913ec4f 1576}
8d2cd130 1577
ec2c406e 1578Bool_t AliGenPythia::IsInEMCAL(Float_t phi, Float_t eta)
1579{
1580 // Is particle in EMCAL acceptance?
ec2c406e 1581 // phi in degrees, etamin=-etamax
9fd8e520 1582 if(phi > fEMCALMinPhi && phi < fEMCALMaxPhi &&
1583 eta < fEMCALEta )
ec2c406e 1584 return kTRUE;
1585 else
1586 return kFALSE;
1587}
1588
1589Bool_t AliGenPythia::IsInPHOS(Float_t phi, Float_t eta)
1590{
1591 // Is particle in PHOS acceptance?
1592 // Acceptance slightly larger considered.
1593 // phi in degrees, etamin=-etamax
9fd8e520 1594 if(phi > fPHOSMinPhi && phi < fPHOSMaxPhi &&
1595 eta < fPHOSEta )
ec2c406e 1596 return kTRUE;
1597 else
1598 return kFALSE;
1599}
1600
90a236ce 1601void AliGenPythia::RotatePhi(Int_t iphcand, Bool_t& okdd)
1602{
1603 //calculate the new position random between fPHOSMinPhi and fPHOSMaxPhi
1604 Double_t phiPHOSmin = TMath::Pi()*fPHOSMinPhi/180;
1605 Double_t phiPHOSmax = TMath::Pi()*fPHOSMaxPhi/180;
1606 Double_t phiPHOS = gRandom->Uniform(phiPHOSmin,phiPHOSmax);
1607
1608 //calculate deltaphi
8507138f 1609 TParticle* ph = (TParticle *) fParticles.At(iphcand);
90a236ce 1610 Double_t phphi = ph->Phi();
1611 Double_t deltaphi = phiPHOS - phphi;
1612
1613
1614
1615 //loop for all particles and produce the phi rotation
8507138f 1616 Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
90a236ce 1617 Double_t oldphi, newphi;
1618 Double_t newVx, newVy, R, Vz, time;
1619 Double_t newPx, newPy, pt, Pz, e;
1620 for(Int_t i=0; i< np; i++) {
8507138f 1621 TParticle* iparticle = (TParticle *) fParticles.At(i);
90a236ce 1622 oldphi = iparticle->Phi();
1623 newphi = oldphi + deltaphi;
1624 if(newphi < 0) newphi = 2*TMath::Pi() + newphi; // correct angle
1625 if(newphi > 2*TMath::Pi()) newphi = newphi - 2*TMath::Pi(); // correct angle
1626
1627 R = iparticle->R();
1628 newVx = R*TMath::Cos(newphi);
1629 newVy = R*TMath::Sin(newphi);
1630 Vz = iparticle->Vz(); // don't transform
1631 time = iparticle->T(); // don't transform
1632
1633 pt = iparticle->Pt();
1634 newPx = pt*TMath::Cos(newphi);
1635 newPy = pt*TMath::Sin(newphi);
1636 Pz = iparticle->Pz(); // don't transform
1637 e = iparticle->Energy(); // don't transform
1638
1639 // apply rotation
1640 iparticle->SetProductionVertex(newVx, newVy, Vz, time);
1641 iparticle->SetMomentum(newPx, newPy, Pz, e);
1642
1643 } //end particle loop
1644
1645 // now let's check that we put correctly the candidate photon in PHOS
1646 Float_t phi = ph->Phi()*180./TMath::Pi(); //Convert to degrees
1647 Float_t eta =TMath::Abs(ph->Eta());//in calos etamin=-etamax
1648 if(IsInPHOS(phi,eta))
1649 okdd = kTRUE;
1650}
ec2c406e 1651
1652
8d2cd130 1653#ifdef never
1654void AliGenPythia::Streamer(TBuffer &R__b)
1655{
1656 // Stream an object of class AliGenPythia.
1657
1658 if (R__b.IsReading()) {
1659 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1660 AliGenerator::Streamer(R__b);
1661 R__b >> (Int_t&)fProcess;
1662 R__b >> (Int_t&)fStrucFunc;
1663 R__b >> (Int_t&)fForceDecay;
1664 R__b >> fEnergyCMS;
1665 R__b >> fKineBias;
1666 R__b >> fTrials;
1667 fParentSelect.Streamer(R__b);
1668 fChildSelect.Streamer(R__b);
1669 R__b >> fXsection;
1670// (AliPythia::Instance())->Streamer(R__b);
1671 R__b >> fPtHardMin;
1672 R__b >> fPtHardMax;
1673// if (fDecayer) fDecayer->Streamer(R__b);
1674 } else {
1675 R__b.WriteVersion(AliGenPythia::IsA());
1676 AliGenerator::Streamer(R__b);
1677 R__b << (Int_t)fProcess;
1678 R__b << (Int_t)fStrucFunc;
1679 R__b << (Int_t)fForceDecay;
1680 R__b << fEnergyCMS;
1681 R__b << fKineBias;
1682 R__b << fTrials;
1683 fParentSelect.Streamer(R__b);
1684 fChildSelect.Streamer(R__b);
1685 R__b << fXsection;
1686// R__b << fPythia;
1687 R__b << fPtHardMin;
1688 R__b << fPtHardMax;
1689 // fDecayer->Streamer(R__b);
1690 }
1691}
1692#endif
1693
90d7b703 1694
589380c6 1695