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