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