git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Some code duplication in Pythia configuration removed.
[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
28#include <TDatabasePDG.h>
29#include <TParticle.h>
30#include <TPDGCode.h>
31#include <TSystem.h>
32#include <TTree.h>
8d2cd130 33#include "AliConst.h"
34#include "AliDecayerPythia.h"
35#include "AliGenPythia.h"
5fa4b20b 36#include "AliHeader.h"
8d2cd130 37#include "AliGenPythiaEventHeader.h"
38#include "AliPythia.h"
7cdba479 39#include "AliPythiaRndm.h"
8d2cd130 40#include "AliRun.h"
7ea3ea5b 41#include "AliStack.h"
42#include "AliRunLoader.h"
5d12ce38 43#include "AliMC.h"
7c21f297 44#include "pyquenCommon.h"
8d2cd130 45
014a9521 46ClassImp(AliGenPythia)
8d2cd130 47
48AliGenPythia::AliGenPythia()
49 :AliGenMC()
50{
51// Default Constructor
52 fParticles = 0;
53 fPythia = 0;
e5c87a3d 54 fHeader = 0;
5fffd0fe 55 fReadFromFile = 0;
f913ec4f 56 fEventTime = 0.;
8ec89232 57 fInteractionRate = 0.;
9d7108a7 58 fTimeWindow = 0.;
59 fEventsTime = 0;
60 fCurSubEvent = 0;
8d2cd130 61 fDecayer = new AliDecayerPythia();
62 SetEventListRange();
63 SetJetPhiRange();
64 SetJetEtaRange();
65 SetJetEtRange();
66 SetGammaPhiRange();
67 SetGammaEtaRange();
7ce1321b 68 SetTriggerParticle();
8d2cd130 69 SetPtKick();
7ea3ea5b 70 SetQuench();
5fa4b20b 71 SetHadronisation();
1d568bc2 72 fSetNuclei = kFALSE;
beac474c 73 fNewMIS = kFALSE;
74 fHFoff = kFALSE;
7cdba479 75 if (!AliPythiaRndm::GetPythiaRandom())
76 AliPythiaRndm::SetPythiaRandom(GetRandom());
8d2cd130 77}
78
79AliGenPythia::AliGenPythia(Int_t npart)
80 :AliGenMC(npart)
81{
82// default charm production at 5. 5 TeV
83// semimuonic decay
84// structure function GRVHO
85//
86 fName = "Pythia";
87 fTitle= "Particle Generator using PYTHIA";
88 fXsection = 0.;
5fffd0fe 89 fReadFromFile = 0;
f913ec4f 90 fEventTime = 0.;
8ec89232 91 fInteractionRate = 0.;
9d7108a7 92 fTimeWindow = 0.;
93 fEventsTime = 0;
94 fCurSubEvent = 0;
8d2cd130 95 SetProcess();
96 SetStrucFunc();
97 SetForceDecay();
98 SetPtHard();
99 SetYHard();
100 SetEnergyCMS();
101 fDecayer = new AliDecayerPythia();
102 // Set random number generator
7cdba479 103 if (!AliPythiaRndm::GetPythiaRandom())
104 AliPythiaRndm::SetPythiaRandom(GetRandom());
8d2cd130 105 fFlavorSelect = 0;
106 // Produced particles
107 fParticles = new TClonesArray("TParticle",1000);
e5c87a3d 108 fHeader = 0;
8d2cd130 109 SetEventListRange();
110 SetJetPhiRange();
111 SetJetEtaRange();
112 SetJetEtRange();
113 SetGammaPhiRange();
114 SetGammaEtaRange();
7ce1321b 115 SetTriggerParticle();
8d2cd130 116 SetJetReconstructionMode();
7ea3ea5b 117 SetQuench();
5fa4b20b 118 SetHadronisation();
8d2cd130 119 SetPtKick();
120 // Options determining what to keep in the stack (Heavy flavour generation)
121 fStackFillOpt = kFlavorSelection; // Keep particle with selected flavor
122 fFeedDownOpt = kTRUE; // allow feed down from higher family
123 // Fragmentation on/off
124 fFragmentation = kTRUE;
125 // Default counting mode
126 fCountMode = kCountAll;
592f8307 127 // Pycel
128 SetPycellParameters();
1d568bc2 129 fSetNuclei = kFALSE;
beac474c 130 fNewMIS = kFALSE;
131 fHFoff = kFALSE;
8d2cd130 132}
133
134AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)
f936b8ea 135 :AliGenMC(Pythia)
8d2cd130 136{
137// copy constructor
138 Pythia.Copy(*this);
139}
140
141AliGenPythia::~AliGenPythia()
142{
143// Destructor
9d7108a7 144 if(fEventsTime) delete fEventsTime;
145}
146
147void AliGenPythia::SetInteractionRate(Float_t rate,Float_t timewindow)
148{
149// Generate pileup using user specified rate
150 fInteractionRate = rate;
151 fTimeWindow = timewindow;
152 GeneratePileup();
153}
154
155void AliGenPythia::GeneratePileup()
156{
157// Generate sub events time for pileup
158 fEventsTime = 0;
159 if(fInteractionRate == 0.) {
160 Warning("GeneratePileup","Zero interaction specified. Skipping pileup generation.\n");
161 return;
162 }
163
164 Int_t npart = NumberParticles();
165 if(npart < 0) {
166 Warning("GeneratePileup","Negative number of particles. Skipping pileup generation.\n");
167 return;
168 }
169
170 if(fEventsTime) delete fEventsTime;
171 fEventsTime = new TArrayF(npart);
172 TArrayF &array = *fEventsTime;
173 for(Int_t ipart = 0; ipart < npart; ipart++)
174 array[ipart] = 0.;
175
176 Float_t eventtime = 0.;
177 while(1)
178 {
179 eventtime += (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
180 if(eventtime > fTimeWindow) break;
181 array.Set(array.GetSize()+1);
182 array[array.GetSize()-1] = eventtime;
183 }
184
185 eventtime = 0.;
186 while(1)
187 {
188 eventtime -= (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
189 if(TMath::Abs(eventtime) > fTimeWindow) break;
190 array.Set(array.GetSize()+1);
191 array[array.GetSize()-1] = eventtime;
192 }
193
194 SetNumberParticles(fEventsTime->GetSize());
8d2cd130 195}
196
592f8307 197void AliGenPythia::SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
198 Float_t thresh, Float_t etseed, Float_t minet, Float_t r)
199{
200// Set pycell parameters
201 fPycellEtaMax = etamax;
202 fPycellNEta = neta;
203 fPycellNPhi = nphi;
204 fPycellThreshold = thresh;
205 fPycellEtSeed = etseed;
206 fPycellMinEtJet = minet;
207 fPycellMaxRadius = r;
208}
209
210
211
8d2cd130 212void AliGenPythia::SetEventListRange(Int_t eventFirst, Int_t eventLast)
213{
214 // Set a range of event numbers, for which a table
215 // of generated particle will be printed
216 fDebugEventFirst = eventFirst;
217 fDebugEventLast = eventLast;
218 if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst;
219}
220
221void AliGenPythia::Init()
222{
223// Initialisation
224
225 SetMC(AliPythia::Instance());
b88f5cea 226 fPythia=(AliPythia*) fMCEvGen;
e2bddf81 227
8d2cd130 228//
229 fParentWeight=1./Float_t(fNpart);
230//
231// Forward Paramters to the AliPythia object
232 fDecayer->SetForceDecay(fForceDecay);
233 fDecayer->Init();
234
235
236 fPythia->SetCKIN(3,fPtHardMin);
237 fPythia->SetCKIN(4,fPtHardMax);
238 fPythia->SetCKIN(7,fYHardMin);
239 fPythia->SetCKIN(8,fYHardMax);
240
1a626d4e 241 if (fAProjectile > 0 && fATarget > 0) fPythia->SetNuclei(fAProjectile, fATarget);
8d2cd130 242 // Fragmentation?
243 if (fFragmentation) {
244 fPythia->SetMSTP(111,1);
245 } else {
246 fPythia->SetMSTP(111,0);
247 }
248
249
250// initial state radiation
251 fPythia->SetMSTP(61,fGinit);
252// final state radiation
253 fPythia->SetMSTP(71,fGfinal);
254// pt - kick
255 if (fPtKick > 0.) {
256 fPythia->SetMSTP(91,1);
257 fPythia->SetPARP(91,fPtKick);
258 } else {
259 fPythia->SetMSTP(91,0);
260 }
261
5fa4b20b 262
263 if (fReadFromFile) {
264 fRL = AliRunLoader::Open(fFileName, "Partons");
265 fRL->LoadKinematics();
266 fRL->LoadHeader();
267 } else {
268 fRL = 0x0;
269 }
beac474c 270// Switch off Heavy Flavors on request
271 if (fHFoff) {
272 fPythia->SetMSTP(58, 3);
273 fPythia->SetMSTJ(45, 3);
274 for (Int_t i = 156; i <= 160; i++) fPythia->SetMDME(i, 1, 0);
275 }
8d2cd130 276 //
277 fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
278
279// Parent and Children Selection
280 switch (fProcess)
281 {
65f2626c 282 case kPyOldUEQ2ordered:
283 case kPyOldUEQ2ordered2:
284 case kPyOldPopcorn:
285 break;
8d2cd130 286 case kPyCharm:
287 case kPyCharmUnforced:
adf4d898 288 case kPyCharmPbPbMNR:
aabc7187 289 case kPyCharmppMNR:
290 case kPyCharmpPbMNR:
8d2cd130 291 fParentSelect[0] = 411;
292 fParentSelect[1] = 421;
293 fParentSelect[2] = 431;
294 fParentSelect[3] = 4122;
295 fFlavorSelect = 4;
296 break;
adf4d898 297 case kPyD0PbPbMNR:
298 case kPyD0pPbMNR:
299 case kPyD0ppMNR:
8d2cd130 300 fParentSelect[0] = 421;
301 fFlavorSelect = 4;
302 break;
90d7b703 303 case kPyDPlusPbPbMNR:
304 case kPyDPluspPbMNR:
305 case kPyDPlusppMNR:
306 fParentSelect[0] = 411;
307 fFlavorSelect = 4;
308 break;
8d2cd130 309 case kPyBeauty:
adf4d898 310 case kPyBeautyPbPbMNR:
311 case kPyBeautypPbMNR:
312 case kPyBeautyppMNR:
8d2cd130 313 fParentSelect[0]= 511;
314 fParentSelect[1]= 521;
315 fParentSelect[2]= 531;
316 fParentSelect[3]= 5122;
317 fParentSelect[4]= 5132;
318 fParentSelect[5]= 5232;
319 fParentSelect[6]= 5332;
320 fFlavorSelect = 5;
321 break;
322 case kPyBeautyUnforced:
323 fParentSelect[0] = 511;
324 fParentSelect[1] = 521;
325 fParentSelect[2] = 531;
326 fParentSelect[3] = 5122;
327 fParentSelect[4] = 5132;
328 fParentSelect[5] = 5232;
329 fParentSelect[6] = 5332;
330 fFlavorSelect = 5;
331 break;
332 case kPyJpsiChi:
333 case kPyJpsi:
334 fParentSelect[0] = 443;
335 break;
336 case kPyMb:
337 case kPyMbNonDiffr:
338 case kPyJets:
339 case kPyDirectGamma:
340 break;
589380c6 341 case kPyW:
0f6ee828 342 case kPyZ:
589380c6 343 break;
8d2cd130 344 }
345//
592f8307 346//
347// JetFinder for Trigger
348//
349// Configure detector (EMCAL like)
350//
351 fPythia->SetPARU(51, fPycellEtaMax);
352 fPythia->SetMSTU(51, fPycellNEta);
353 fPythia->SetMSTU(52, fPycellNPhi);
354//
355// Configure Jet Finder
356//
357 fPythia->SetPARU(58, fPycellThreshold);
358 fPythia->SetPARU(52, fPycellEtSeed);
359 fPythia->SetPARU(53, fPycellMinEtJet);
360 fPythia->SetPARU(54, fPycellMaxRadius);
361 fPythia->SetMSTU(54, 2);
362//
8d2cd130 363// This counts the total number of calls to Pyevnt() per run.
364 fTrialsRun = 0;
365 fQ = 0.;
366 fX1 = 0.;
367 fX2 = 0.;
368 fNev = 0 ;
369//
1d568bc2 370//
371//
8d2cd130 372 AliGenMC::Init();
1d568bc2 373//
374//
375//
376 if (fSetNuclei) {
377 fDyBoost = 0;
378 Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n");
379 }
32d6ef7d 380
381 if (fQuench) {
5fa4b20b 382 fPythia->InitQuenching(0., 0.1, 0.6e6, 0);
32d6ef7d 383 }
384
8d2cd130 385}
386
387void AliGenPythia::Generate()
388{
389// Generate one event
e5c87a3d 390
8d2cd130 391 fDecayer->ForceDecay();
392
393 Float_t polar[3] = {0,0,0};
394 Float_t origin[3] = {0,0,0};
a920faf9 395 Float_t p[4];
8d2cd130 396// converts from mm/c to s
397 const Float_t kconv=0.001/2.999792458e8;
398//
399 Int_t nt=0;
400 Int_t jev=0;
401 Int_t j, kf;
402 fTrials=0;
f913ec4f 403 fEventTime = 0.;
404
2590ccf9 405
8d2cd130 406
407 // Set collision vertex position
2590ccf9 408 if (fVertexSmear == kPerEvent) Vertex();
409
8d2cd130 410// event loop
411 while(1)
412 {
32d6ef7d 413//
5fa4b20b 414// Produce event
32d6ef7d 415//
32d6ef7d 416//
5fa4b20b 417// Switch hadronisation off
418//
419 fPythia->SetMSTJ(1, 0);
32d6ef7d 420//
5fa4b20b 421// Either produce new event or read partons from file
422//
423 if (!fReadFromFile) {
beac474c 424 if (!fNewMIS) {
425 fPythia->Pyevnt();
426 } else {
427 fPythia->Pyevnw();
428 }
5fa4b20b 429 fNpartons = fPythia->GetN();
430 } else {
431 printf("Loading Event %d\n",AliRunLoader::GetRunLoader()->GetEventNumber());
432 fRL->GetEvent(AliRunLoader::GetRunLoader()->GetEventNumber());
433 fPythia->SetN(0);
434 LoadEvent(fRL->Stack(), 0 , 1);
435 fPythia->Pyedit(21);
436 }
437
32d6ef7d 438//
439// Run quenching routine
440//
5fa4b20b 441 if (fQuench == 1) {
442 fPythia->Quench();
443 } else if (fQuench == 2){
444 fPythia->Pyquen(208., 0, 0.);
445 }
32d6ef7d 446//
5fa4b20b 447// Switch hadronisation on
32d6ef7d 448//
aea21c57 449 fPythia->SetMSTJ(1, 1);
5fa4b20b 450//
451// .. and perform hadronisation
aea21c57 452// printf("Calling hadronisation %d\n", fPythia->GetN());
5fa4b20b 453 fPythia->Pyexec();
8d2cd130 454 fTrials++;
8d2cd130 455 fPythia->ImportParticles(fParticles,"All");
1d568bc2 456 Boost();
8d2cd130 457//
458//
459//
460 Int_t i;
461
5fa4b20b 462
8d2cd130 463 Int_t np = fParticles->GetEntriesFast();
5fa4b20b 464
7c21f297 465 if (np == 0) continue;
8d2cd130 466//
2590ccf9 467
8d2cd130 468//
469 Int_t* pParent = new Int_t[np];
470 Int_t* pSelected = new Int_t[np];
471 Int_t* trackIt = new Int_t[np];
5fa4b20b 472 for (i = 0; i < np; i++) {
8d2cd130 473 pParent[i] = -1;
474 pSelected[i] = 0;
475 trackIt[i] = 0;
476 }
477
478 Int_t nc = 0; // Total n. of selected particles
479 Int_t nParents = 0; // Selected parents
480 Int_t nTkbles = 0; // Trackable particles
481 if (fProcess != kPyMb && fProcess != kPyJets &&
482 fProcess != kPyDirectGamma &&
589380c6 483 fProcess != kPyMbNonDiffr &&
0f6ee828 484 fProcess != kPyW && fProcess != kPyZ ) {
8d2cd130 485
5fa4b20b 486 for (i = 0; i < np; i++) {
2590ccf9 487 TParticle* iparticle = (TParticle *) fParticles->At(i);
8d2cd130 488 Int_t ks = iparticle->GetStatusCode();
489 kf = CheckPDGCode(iparticle->GetPdgCode());
490// No initial state partons
491 if (ks==21) continue;
492//
493// Heavy Flavor Selection
494//
495 // quark ?
496 kf = TMath::Abs(kf);
497 Int_t kfl = kf;
9ff6c04c 498 // Resonance
f913ec4f 499
9ff6c04c 500 if (kfl > 100000) kfl %= 100000;
183a5ca9 501 if (kfl > 10000) kfl %= 10000;
8d2cd130 502 // meson ?
503 if (kfl > 10) kfl/=100;
504 // baryon
505 if (kfl > 10) kfl/=10;
8d2cd130 506 Int_t ipa = iparticle->GetFirstMother()-1;
507 Int_t kfMo = 0;
f913ec4f 508//
509// Establish mother daughter relation between heavy quarks and mesons
510//
511 if (kf >= fFlavorSelect && kf <= 6) {
512 Int_t idau = iparticle->GetFirstDaughter() - 1;
513 if (idau > -1) {
514 TParticle* daughter = (TParticle *) fParticles->At(idau);
515 Int_t pdgD = daughter->GetPdgCode();
516 if (pdgD == 91 || pdgD == 92) {
517 Int_t jmin = daughter->GetFirstDaughter() - 1;
518 Int_t jmax = daughter->GetLastDaughter() - 1;
519 for (Int_t j = jmin; j <= jmax; j++)
520 ((TParticle *) fParticles->At(j))->SetFirstMother(i+1);
521 } // is string or cluster
522 } // has daughter
523 } // heavy quark
8d2cd130 524
f913ec4f 525
8d2cd130 526 if (ipa > -1) {
527 TParticle * mother = (TParticle *) fParticles->At(ipa);
528 kfMo = TMath::Abs(mother->GetPdgCode());
529 }
f913ec4f 530
8d2cd130 531 // What to keep in Stack?
532 Bool_t flavorOK = kFALSE;
533 Bool_t selectOK = kFALSE;
534 if (fFeedDownOpt) {
32d6ef7d 535 if (kfl >= fFlavorSelect) flavorOK = kTRUE;
8d2cd130 536 } else {
32d6ef7d 537 if (kfl > fFlavorSelect) {
538 nc = -1;
539 break;
540 }
541 if (kfl == fFlavorSelect) flavorOK = kTRUE;
8d2cd130 542 }
543 switch (fStackFillOpt) {
544 case kFlavorSelection:
32d6ef7d 545 selectOK = kTRUE;
546 break;
8d2cd130 547 case kParentSelection:
32d6ef7d 548 if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
549 break;
8d2cd130 550 }
551 if (flavorOK && selectOK) {
552//
553// Heavy flavor hadron or quark
554//
555// Kinematic seletion on final state heavy flavor mesons
556 if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
557 {
9ff6c04c 558 continue;
8d2cd130 559 }
560 pSelected[i] = 1;
561 if (ParentSelected(kf)) ++nParents; // Update parent count
562// printf("\n particle (HF) %d %d %d", i, pSelected[i], kf);
563 } else {
564// Kinematic seletion on decay products
565 if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf)
9ff6c04c 566 && !KinematicSelection(iparticle, 1))
8d2cd130 567 {
9ff6c04c 568 continue;
8d2cd130 569 }
570//
571// Decay products
572// Select if mother was selected and is not tracked
573
574 if (pSelected[ipa] &&
575 !trackIt[ipa] && // mother will be tracked ?
576 kfMo != 5 && // mother is b-quark, don't store fragments
577 kfMo != 4 && // mother is c-quark, don't store fragments
578 kf != 92) // don't store string
579 {
580//
581// Semi-stable or de-selected: diselect decay products:
582//
583//
584 if (pSelected[i] == -1 || fDecayer->GetLifetime(kf) > fMaxLifeTime)
585 {
586 Int_t ipF = iparticle->GetFirstDaughter();
587 Int_t ipL = iparticle->GetLastDaughter();
588 if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
589 }
590// printf("\n particle (decay) %d %d %d", i, pSelected[i], kf);
591 pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
592 }
593 }
594 if (pSelected[i] == -1) pSelected[i] = 0;
595 if (!pSelected[i]) continue;
596 // Count quarks only if you did not include fragmentation
597 if (fFragmentation && kf <= 10) continue;
9ff6c04c 598
8d2cd130 599 nc++;
600// Decision on tracking
601 trackIt[i] = 0;
602//
603// Track final state particle
604 if (ks == 1) trackIt[i] = 1;
605// Track semi-stable particles
d25cfd65 606 if ((ks == 1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime)) trackIt[i] = 1;
8d2cd130 607// Track particles selected by process if undecayed.
608 if (fForceDecay == kNoDecay) {
609 if (ParentSelected(kf)) trackIt[i] = 1;
610 } else {
611 if (ParentSelected(kf)) trackIt[i] = 0;
612 }
613 if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
614//
615//
616
617 } // particle selection loop
618 if (nc > 0) {
619 for (i = 0; i<np; i++) {
620 if (!pSelected[i]) continue;
621 TParticle * iparticle = (TParticle *) fParticles->At(i);
622 kf = CheckPDGCode(iparticle->GetPdgCode());
623 Int_t ks = iparticle->GetStatusCode();
624 p[0] = iparticle->Px();
625 p[1] = iparticle->Py();
626 p[2] = iparticle->Pz();
a920faf9 627 p[3] = iparticle->Energy();
628
2590ccf9 629 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
630 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
631 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
632
8d2cd130 633 Float_t tof = kconv*iparticle->T();
634 Int_t ipa = iparticle->GetFirstMother()-1;
635 Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
a920faf9 636
637 PushTrack(fTrackIt*trackIt[i], iparent, kf,
638 p[0], p[1], p[2], p[3],
639 origin[0], origin[1], origin[2], tof,
640 polar[0], polar[1], polar[2],
641 kPPrimary, nt, 1., ks);
8d2cd130 642 pParent[i] = nt;
643 KeepTrack(nt);
642f15cf 644 } // PushTrack loop
8d2cd130 645 }
646 } else {
647 nc = GenerateMB();
648 } // mb ?
f913ec4f 649
650 GetSubEventTime();
8d2cd130 651
652 if (pParent) delete[] pParent;
653 if (pSelected) delete[] pSelected;
654 if (trackIt) delete[] trackIt;
655
656 if (nc > 0) {
657 switch (fCountMode) {
658 case kCountAll:
659 // printf(" Count all \n");
660 jev += nc;
661 break;
662 case kCountParents:
663 // printf(" Count parents \n");
664 jev += nParents;
665 break;
666 case kCountTrackables:
667 // printf(" Count trackable \n");
668 jev += nTkbles;
669 break;
670 }
671 if (jev >= fNpart || fNpart == -1) {
672 fKineBias=Float_t(fNpart)/Float_t(fTrials);
8d2cd130 673 fQ += fPythia->GetVINT(51);
674 fX1 += fPythia->GetVINT(41);
675 fX2 += fPythia->GetVINT(42);
676 fTrialsRun += fTrials;
677 fNev++;
678 MakeHeader();
679 break;
680 }
681 }
682 } // event loop
683 SetHighWaterMark(nt);
684// adjust weight due to kinematic selection
b88f5cea 685// AdjustWeights();
8d2cd130 686// get cross-section
687 fXsection=fPythia->GetPARI(1);
688}
689
690Int_t AliGenPythia::GenerateMB()
691{
692//
693// Min Bias selection and other global selections
694//
695 Int_t i, kf, nt, iparent;
696 Int_t nc = 0;
bf950da8 697 Float_t p[4];
8d2cd130 698 Float_t polar[3] = {0,0,0};
699 Float_t origin[3] = {0,0,0};
700// converts from mm/c to s
701 const Float_t kconv=0.001/2.999792458e8;
702
5fa4b20b 703 Int_t np = (fHadronisation) ? fParticles->GetEntriesFast() : fNpartons;
aea21c57 704
5fa4b20b 705
8d2cd130 706 Int_t* pParent = new Int_t[np];
707 for (i=0; i< np; i++) pParent[i] = -1;
708 if (fProcess == kPyJets || fProcess == kPyDirectGamma) {
709 TParticle* jet1 = (TParticle *) fParticles->At(6);
710 TParticle* jet2 = (TParticle *) fParticles->At(7);
aea21c57 711 if (!CheckTrigger(jet1, jet2)) return 0;
8d2cd130 712 }
713
7ce1321b 714 if (fTriggerParticle) {
715 Bool_t triggered = kFALSE;
716 for (i = 0; i < np; i++) {
717 TParticle * iparticle = (TParticle *) fParticles->At(i);
718 kf = CheckPDGCode(iparticle->GetPdgCode());
719 if (TMath::Abs(kf) != fTriggerParticle) continue;
720 if (iparticle->Pt() == 0.) continue;
721 if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
722 triggered = kTRUE;
723 break;
724 }
725 if (!triggered) return 0;
726 }
727
aea21c57 728
0f6ee828 729 //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
730 if ( (fProcess == kPyW || fProcess == kPyZ || fProcess == kPyMb || fProcess == kPyMbNonDiffr)
731 && (fCutOnChild == 1) ) {
732 if ( !CheckKinematicsOnChild() ) {
733 if (pParent) delete[] pParent;
734 return 0;
735 }
aea21c57 736 }
737
738
7ce1321b 739
f913ec4f 740 for (i = 0; i < np; i++) {
8d2cd130 741 Int_t trackIt = 0;
742 TParticle * iparticle = (TParticle *) fParticles->At(i);
743 kf = CheckPDGCode(iparticle->GetPdgCode());
744 Int_t ks = iparticle->GetStatusCode();
745 Int_t km = iparticle->GetFirstMother();
746 if ((ks == 1 && kf!=0 && KinematicSelection(iparticle, 0)) ||
747 (ks != 1) ||
748 (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
749 nc++;
750 if (ks == 1) trackIt = 1;
751 Int_t ipa = iparticle->GetFirstMother()-1;
752
753 iparent = (ipa > -1) ? pParent[ipa] : -1;
754
755//
756// store track information
757 p[0] = iparticle->Px();
758 p[1] = iparticle->Py();
759 p[2] = iparticle->Pz();
a920faf9 760 p[3] = iparticle->Energy();
1406f599 761
a920faf9 762
2590ccf9 763 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
764 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
765 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
766
f913ec4f 767 Float_t tof = fEventTime + kconv * iparticle->T();
a920faf9 768
769 PushTrack(fTrackIt*trackIt, iparent, kf,
770 p[0], p[1], p[2], p[3],
771 origin[0], origin[1], origin[2], tof,
772 polar[0], polar[1], polar[2],
773 kPPrimary, nt, 1., ks);
5fa4b20b 774 //
775 // Special Treatment to store color-flow
776 //
777 if (ks == 3 || ks == 13 || ks == 14) {
778 TParticle* particle = 0;
779 if (fStack) {
780 particle = fStack->Particle(nt);
781 } else {
782 particle = gAlice->Stack()->Particle(nt);
783 }
784 particle->SetFirstDaughter(fPythia->GetK(2, i));
785 particle->SetLastDaughter(fPythia->GetK(3, i));
786 }
787
8d2cd130 788 KeepTrack(nt);
789 pParent[i] = nt;
f913ec4f 790 SetHighWaterMark(nt);
791
8d2cd130 792 } // select particle
793 } // particle loop
794
795 if (pParent) delete[] pParent;
f913ec4f 796 return 1;
8d2cd130 797}
798
799
800void AliGenPythia::FinishRun()
801{
802// Print x-section summary
803 fPythia->Pystat(1);
2779fc64 804
805 if (fNev > 0.) {
806 fQ /= fNev;
807 fX1 /= fNev;
808 fX2 /= fNev;
809 }
810
8d2cd130 811 printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
812 printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
8d2cd130 813}
814
815void AliGenPythia::AdjustWeights()
816{
817// Adjust the weights after generation of all events
818//
e2bddf81 819 if (gAlice) {
820 TParticle *part;
821 Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
822 for (Int_t i=0; i<ntrack; i++) {
823 part= gAlice->GetMCApp()->Particle(i);
824 part->SetWeight(part->GetWeight()*fKineBias);
825 }
8d2cd130 826 }
827}
828
829void AliGenPythia::SetNuclei(Int_t a1, Int_t a2)
830{
831// Treat protons as inside nuclei with mass numbers a1 and a2
1d568bc2 832
1a626d4e 833 fAProjectile = a1;
834 fATarget = a2;
1d568bc2 835 fSetNuclei = kTRUE;
8d2cd130 836}
837
838
839void AliGenPythia::MakeHeader()
840{
183a5ca9 841 if (gAlice) {
842 if (gAlice->GetEvNumber()>=fDebugEventFirst &&
f913ec4f 843 gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(2);
183a5ca9 844 }
845
8d2cd130 846// Builds the event header, to be called after each event
e5c87a3d 847 if (fHeader) delete fHeader;
848 fHeader = new AliGenPythiaEventHeader("Pythia");
8d2cd130 849//
850// Event type
e5c87a3d 851 ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
8d2cd130 852//
853// Number of trials
e5c87a3d 854 ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
8d2cd130 855//
856// Event Vertex
d25cfd65 857 fHeader->SetPrimaryVertex(fVertex);
8d2cd130 858//
859// Jets that have triggered
f913ec4f 860
8d2cd130 861 if (fProcess == kPyJets)
862 {
863 Int_t ntrig, njet;
864 Float_t jets[4][10];
865 GetJets(njet, ntrig, jets);
9ff6c04c 866
8d2cd130 867
868 for (Int_t i = 0; i < ntrig; i++) {
e5c87a3d 869 ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
8d2cd130 870 jets[3][i]);
871 }
872 }
5fa4b20b 873//
874// Copy relevant information from external header, if present.
875//
876 Float_t uqJet[4];
877
878 if (fRL) {
879 AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
880 for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
881 {
882 printf("Adding Jet %d %d \n", i, exHeader->NTriggerJets());
883
884
885 exHeader->TriggerJet(i, uqJet);
886 ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
887 }
888 }
889//
890// Store quenching parameters
891//
892 if (fQuench){
893 Double_t z[4];
894 Double_t xp, yp;
7c21f297 895 if (fQuench == 1) {
896 // Pythia::Quench()
897 fPythia->GetQuenchingParameters(xp, yp, z);
898 } else {
899 // Pyquen
900 Double_t r1 = PARIMP.rb1;
901 Double_t r2 = PARIMP.rb2;
902 Double_t b = PARIMP.b1;
903 Double_t r = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
904 Double_t phi = PARIMP.psib1;
905 xp = r * TMath::Cos(phi);
906 yp = r * TMath::Sin(phi);
907
908 }
909 ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
910 ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
911 }
beac474c 912//
913// Store pt^hard
914 ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
5fa4b20b 915//
cf57b268 916// Pass header
5fa4b20b 917//
cf57b268 918 AddHeader(fHeader);
8d2cd130 919}
cf57b268 920
921void AliGenPythia::AddHeader(AliGenEventHeader* header)
922{
923 // Add header to container or runloader
924 if (fContainer) {
925 fContainer->AddHeader(header);
926 } else {
927 AliRunLoader::GetRunLoader()->GetHeader()->SetGenEventHeader(header);
928 }
929}
930
8d2cd130 931
932Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
933{
934// Check the kinematic trigger condition
935//
936 Double_t eta[2];
937 eta[0] = jet1->Eta();
938 eta[1] = jet2->Eta();
939 Double_t phi[2];
940 phi[0] = jet1->Phi();
941 phi[1] = jet2->Phi();
942 Int_t pdg[2];
943 pdg[0] = jet1->GetPdgCode();
944 pdg[1] = jet2->GetPdgCode();
945 Bool_t triggered = kFALSE;
946
947 if (fProcess == kPyJets) {
948 Int_t njets = 0;
949 Int_t ntrig = 0;
950 Float_t jets[4][10];
951//
952// Use Pythia clustering on parton level to determine jet axis
953//
954 GetJets(njets, ntrig, jets);
955
956 if (ntrig) triggered = kTRUE;
957//
958 } else {
959 Int_t ij = 0;
960 Int_t ig = 1;
961 if (pdg[0] == kGamma) {
962 ij = 1;
963 ig = 0;
964 }
965 //Check eta range first...
966 if ((eta[ij] < fEtaMaxJet && eta[ij] > fEtaMinJet) &&
967 (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
968 {
969 //Eta is okay, now check phi range
970 if ((phi[ij] < fPhiMaxJet && phi[ij] > fPhiMinJet) &&
971 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
972 {
973 triggered = kTRUE;
974 }
975 }
976 }
977 return triggered;
978}
aea21c57 979
980
981//Checking Kinematics on Child (status code 1, particle code ?, kin cuts
982Bool_t AliGenPythia::CheckKinematicsOnChild(){
983
984 Bool_t checking = kFALSE;
985 Int_t j, kcode, ks, km;
986 Int_t nPartAcc = 0; //number of particles in the acceptance range
987 Int_t numberOfAcceptedParticles = 1;
988 if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
989 Int_t npart = fParticles->GetEntriesFast();
990
0f6ee828 991 for (j = 0; j<npart; j++) {
aea21c57 992 TParticle * jparticle = (TParticle *) fParticles->At(j);
993 kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
994 ks = jparticle->GetStatusCode();
995 km = jparticle->GetFirstMother();
996
997 if( (ks == 1) && (kcode == fPdgCodeParticleforAcceptanceCut) && (KinematicSelection(jparticle,1)) ){
998 nPartAcc++;
999 }
0f6ee828 1000 if( numberOfAcceptedParticles <= nPartAcc){
1001 checking = kTRUE;
1002 break;
1003 }
aea21c57 1004 }
0f6ee828 1005
aea21c57 1006 return checking;
aea21c57 1007}
1008
8d2cd130 1009
1010AliGenPythia& AliGenPythia::operator=(const AliGenPythia& rhs)
1011{
1012// Assignment operator
014a9521 1013 rhs.Copy(*this);
8d2cd130 1014 return *this;
1015}
1016
5fa4b20b 1017void AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
8d2cd130 1018{
1019//
1020// Load event into Pythia Common Block
1021//
5fa4b20b 1022
32d6ef7d 1023 Int_t npart = stack -> GetNprimary();
1024 Int_t n0 = 0;
1025
1026 if (!flag) {
1027 (fPythia->GetPyjets())->N = npart;
1028 } else {
1029 n0 = (fPythia->GetPyjets())->N;
1030 (fPythia->GetPyjets())->N = n0 + npart;
1031 }
1032
1033
8d2cd130 1034 for (Int_t part = 0; part < npart; part++) {
32d6ef7d 1035 TParticle *MPart = stack->Particle(part);
1036
5fa4b20b 1037 Int_t kf = MPart->GetPdgCode();
1038 Int_t ks = MPart->GetStatusCode();
1039 Int_t idf = MPart->GetFirstDaughter();
1040 Int_t idl = MPart->GetLastDaughter();
1041
1042 if (reHadr) {
1043 if (ks == 11 || ks == 12) {
1044 ks -= 10;
1045 idf = -1;
1046 idl = -1;
1047 }
1048 }
32d6ef7d 1049
8d2cd130 1050 Float_t px = MPart->Px();
1051 Float_t py = MPart->Py();
1052 Float_t pz = MPart->Pz();
1053 Float_t e = MPart->Energy();
a920faf9 1054 Float_t m = MPart->GetCalcMass();
8d2cd130 1055
1056
32d6ef7d 1057 (fPythia->GetPyjets())->P[0][part+n0] = px;
1058 (fPythia->GetPyjets())->P[1][part+n0] = py;
1059 (fPythia->GetPyjets())->P[2][part+n0] = pz;
1060 (fPythia->GetPyjets())->P[3][part+n0] = e;
1061 (fPythia->GetPyjets())->P[4][part+n0] = m;
8d2cd130 1062
32d6ef7d 1063 (fPythia->GetPyjets())->K[1][part+n0] = kf;
1064 (fPythia->GetPyjets())->K[0][part+n0] = ks;
5fa4b20b 1065 (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1066 (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1067 (fPythia->GetPyjets())->K[2][part+n0] = MPart->GetFirstMother() + 1;
8d2cd130 1068 }
1069}
1070
5fa4b20b 1071
014a9521 1072void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
8d2cd130 1073{
1074//
1075// Calls the Pythia jet finding algorithm to find jets in the current event
1076//
1077//
8d2cd130 1078//
1079// Save jets
1080 Int_t n = fPythia->GetN();
1081
1082//
1083// Run Jet Finder
1084 fPythia->Pycell(njets);
1085 Int_t i;
1086 for (i = 0; i < njets; i++) {
1087 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
1088 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
1089 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
1090 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
1091
1092 jets[0][i] = px;
1093 jets[1][i] = py;
1094 jets[2][i] = pz;
1095 jets[3][i] = e;
1096 }
1097}
1098
1099
1100
1101void AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
1102{
1103//
1104// Calls the Pythia clustering algorithm to find jets in the current event
1105//
1106 Int_t n = fPythia->GetN();
1107 nJets = 0;
1108 nJetsTrig = 0;
1109 if (fJetReconstruction == kCluster) {
1110//
1111// Configure cluster algorithm
1112//
1113 fPythia->SetPARU(43, 2.);
1114 fPythia->SetMSTU(41, 1);
1115//
1116// Call cluster algorithm
1117//
1118 fPythia->Pyclus(nJets);
1119//
1120// Loading jets from common block
1121//
1122 } else {
592f8307 1123
8d2cd130 1124//
1125// Run Jet Finder
1126 fPythia->Pycell(nJets);
1127 }
1128
1129 Int_t i;
1130 for (i = 0; i < nJets; i++) {
1131 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
1132 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
1133 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
1134 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
1135 Float_t pt = TMath::Sqrt(px * px + py * py);
a920faf9 1136 Float_t phi = TMath::Pi() + TMath::ATan2(-py, -px);
8d2cd130 1137 Float_t theta = TMath::ATan2(pt,pz);
1138 Float_t et = e * TMath::Sin(theta);
1139 Float_t eta = -TMath::Log(TMath::Tan(theta / 2.));
8d2cd130 1140 if (
1141 eta > fEtaMinJet && eta < fEtaMaxJet &&
675eb105 1142 phi > fPhiMinJet && phi < fPhiMaxJet &&
8d2cd130 1143 et > fEtMinJet && et < fEtMaxJet
1144 )
1145 {
1146 jets[0][nJetsTrig] = px;
1147 jets[1][nJetsTrig] = py;
1148 jets[2][nJetsTrig] = pz;
1149 jets[3][nJetsTrig] = e;
1150 nJetsTrig++;
5fa4b20b 1151// printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
8d2cd130 1152 } else {
1153// printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1154 }
1155 }
1156}
1157
f913ec4f 1158void AliGenPythia::GetSubEventTime()
1159{
1160 // Calculates time of the next subevent
9d7108a7 1161 fEventTime = 0.;
1162 if (fEventsTime) {
1163 TArrayF &array = *fEventsTime;
1164 fEventTime = array[fCurSubEvent++];
1165 }
1166 // printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
1167 return;
f913ec4f 1168}
8d2cd130 1169
1170#ifdef never
1171void AliGenPythia::Streamer(TBuffer &R__b)
1172{
1173 // Stream an object of class AliGenPythia.
1174
1175 if (R__b.IsReading()) {
1176 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1177 AliGenerator::Streamer(R__b);
1178 R__b >> (Int_t&)fProcess;
1179 R__b >> (Int_t&)fStrucFunc;
1180 R__b >> (Int_t&)fForceDecay;
1181 R__b >> fEnergyCMS;
1182 R__b >> fKineBias;
1183 R__b >> fTrials;
1184 fParentSelect.Streamer(R__b);
1185 fChildSelect.Streamer(R__b);
1186 R__b >> fXsection;
1187// (AliPythia::Instance())->Streamer(R__b);
1188 R__b >> fPtHardMin;
1189 R__b >> fPtHardMax;
1190// if (fDecayer) fDecayer->Streamer(R__b);
1191 } else {
1192 R__b.WriteVersion(AliGenPythia::IsA());
1193 AliGenerator::Streamer(R__b);
1194 R__b << (Int_t)fProcess;
1195 R__b << (Int_t)fStrucFunc;
1196 R__b << (Int_t)fForceDecay;
1197 R__b << fEnergyCMS;
1198 R__b << fKineBias;
1199 R__b << fTrials;
1200 fParentSelect.Streamer(R__b);
1201 fChildSelect.Streamer(R__b);
1202 R__b << fXsection;
1203// R__b << fPythia;
1204 R__b << fPtHardMin;
1205 R__b << fPtHardMax;
1206 // fDecayer->Streamer(R__b);
1207 }
1208}
1209#endif
1210
90d7b703 1211
589380c6 1212