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