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