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