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