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