PDC'06 configurations for charm and beauty added.
[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();
e0e89f40 71 SetTriggerParticle();
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();
115 SetJetReconstructionMode();
7ea3ea5b 116 SetQuench();
5fa4b20b 117 SetHadronisation();
8d2cd130 118 SetPtKick();
e0e89f40 119 SetTriggerParticle();
8d2cd130 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 kPyCharmpPbMNR:
e0e89f40 290 case kPyCharmppMNR:
291 case kPyCharmppMNRwmi:
8d2cd130 292 fParentSelect[0] = 411;
293 fParentSelect[1] = 421;
294 fParentSelect[2] = 431;
295 fParentSelect[3] = 4122;
296 fFlavorSelect = 4;
297 break;
adf4d898 298 case kPyD0PbPbMNR:
299 case kPyD0pPbMNR:
300 case kPyD0ppMNR:
8d2cd130 301 fParentSelect[0] = 421;
302 fFlavorSelect = 4;
303 break;
90d7b703 304 case kPyDPlusPbPbMNR:
305 case kPyDPluspPbMNR:
306 case kPyDPlusppMNR:
307 fParentSelect[0] = 411;
308 fFlavorSelect = 4;
309 break;
e0e89f40 310 case kPyDPlusStrangePbPbMNR:
311 case kPyDPlusStrangepPbMNR:
312 case kPyDPlusStrangeppMNR:
313 fParentSelect[0] = 431;
314 fFlavorSelect = 4;
315 break;
8d2cd130 316 case kPyBeauty:
adf4d898 317 case kPyBeautyPbPbMNR:
318 case kPyBeautypPbMNR:
319 case kPyBeautyppMNR:
e0e89f40 320 case kPyBeautyppMNRwmi:
8d2cd130 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 kPyBeautyUnforced:
331 fParentSelect[0] = 511;
332 fParentSelect[1] = 521;
333 fParentSelect[2] = 531;
334 fParentSelect[3] = 5122;
335 fParentSelect[4] = 5132;
336 fParentSelect[5] = 5232;
337 fParentSelect[6] = 5332;
338 fFlavorSelect = 5;
339 break;
340 case kPyJpsiChi:
341 case kPyJpsi:
342 fParentSelect[0] = 443;
343 break;
344 case kPyMb:
345 case kPyMbNonDiffr:
346 case kPyJets:
347 case kPyDirectGamma:
348 break;
589380c6 349 case kPyW:
0f6ee828 350 case kPyZ:
589380c6 351 break;
8d2cd130 352 }
353//
592f8307 354//
355// JetFinder for Trigger
356//
357// Configure detector (EMCAL like)
358//
359 fPythia->SetPARU(51, fPycellEtaMax);
360 fPythia->SetMSTU(51, fPycellNEta);
361 fPythia->SetMSTU(52, fPycellNPhi);
362//
363// Configure Jet Finder
364//
365 fPythia->SetPARU(58, fPycellThreshold);
366 fPythia->SetPARU(52, fPycellEtSeed);
367 fPythia->SetPARU(53, fPycellMinEtJet);
368 fPythia->SetPARU(54, fPycellMaxRadius);
369 fPythia->SetMSTU(54, 2);
370//
8d2cd130 371// This counts the total number of calls to Pyevnt() per run.
372 fTrialsRun = 0;
373 fQ = 0.;
374 fX1 = 0.;
375 fX2 = 0.;
376 fNev = 0 ;
377//
1d568bc2 378//
379//
8d2cd130 380 AliGenMC::Init();
1d568bc2 381//
382//
383//
384 if (fSetNuclei) {
385 fDyBoost = 0;
386 Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n");
387 }
32d6ef7d 388
389 if (fQuench) {
5fa4b20b 390 fPythia->InitQuenching(0., 0.1, 0.6e6, 0);
32d6ef7d 391 }
392
8d2cd130 393}
394
395void AliGenPythia::Generate()
396{
397// Generate one event
e5c87a3d 398
8d2cd130 399 fDecayer->ForceDecay();
400
401 Float_t polar[3] = {0,0,0};
402 Float_t origin[3] = {0,0,0};
a920faf9 403 Float_t p[4];
8d2cd130 404// converts from mm/c to s
405 const Float_t kconv=0.001/2.999792458e8;
406//
407 Int_t nt=0;
408 Int_t jev=0;
409 Int_t j, kf;
410 fTrials=0;
f913ec4f 411 fEventTime = 0.;
412
2590ccf9 413
8d2cd130 414
415 // Set collision vertex position
2590ccf9 416 if (fVertexSmear == kPerEvent) Vertex();
417
8d2cd130 418// event loop
419 while(1)
420 {
32d6ef7d 421//
5fa4b20b 422// Produce event
32d6ef7d 423//
32d6ef7d 424//
5fa4b20b 425// Switch hadronisation off
426//
427 fPythia->SetMSTJ(1, 0);
32d6ef7d 428//
5fa4b20b 429// Either produce new event or read partons from file
430//
431 if (!fReadFromFile) {
beac474c 432 if (!fNewMIS) {
433 fPythia->Pyevnt();
434 } else {
435 fPythia->Pyevnw();
436 }
5fa4b20b 437 fNpartons = fPythia->GetN();
438 } else {
439 printf("Loading Event %d\n",AliRunLoader::GetRunLoader()->GetEventNumber());
440 fRL->GetEvent(AliRunLoader::GetRunLoader()->GetEventNumber());
441 fPythia->SetN(0);
442 LoadEvent(fRL->Stack(), 0 , 1);
443 fPythia->Pyedit(21);
444 }
445
32d6ef7d 446//
447// Run quenching routine
448//
5fa4b20b 449 if (fQuench == 1) {
450 fPythia->Quench();
451 } else if (fQuench == 2){
452 fPythia->Pyquen(208., 0, 0.);
453 }
32d6ef7d 454//
5fa4b20b 455// Switch hadronisation on
32d6ef7d 456//
aea21c57 457 fPythia->SetMSTJ(1, 1);
5fa4b20b 458//
459// .. and perform hadronisation
aea21c57 460// printf("Calling hadronisation %d\n", fPythia->GetN());
5fa4b20b 461 fPythia->Pyexec();
8d2cd130 462 fTrials++;
8d2cd130 463 fPythia->ImportParticles(fParticles,"All");
1d568bc2 464 Boost();
8d2cd130 465//
466//
467//
468 Int_t i;
469
5fa4b20b 470
8d2cd130 471 Int_t np = fParticles->GetEntriesFast();
5fa4b20b 472
7c21f297 473 if (np == 0) continue;
8d2cd130 474//
2590ccf9 475
8d2cd130 476//
477 Int_t* pParent = new Int_t[np];
478 Int_t* pSelected = new Int_t[np];
479 Int_t* trackIt = new Int_t[np];
5fa4b20b 480 for (i = 0; i < np; i++) {
8d2cd130 481 pParent[i] = -1;
482 pSelected[i] = 0;
483 trackIt[i] = 0;
484 }
485
486 Int_t nc = 0; // Total n. of selected particles
487 Int_t nParents = 0; // Selected parents
488 Int_t nTkbles = 0; // Trackable particles
489 if (fProcess != kPyMb && fProcess != kPyJets &&
490 fProcess != kPyDirectGamma &&
589380c6 491 fProcess != kPyMbNonDiffr &&
e0e89f40 492 fProcess != kPyW && fProcess != kPyZ &&
493 fProcess != kPyCharmppMNRwmi && fProcess != kPyBeautyppMNRwmi) {
8d2cd130 494
5fa4b20b 495 for (i = 0; i < np; i++) {
2590ccf9 496 TParticle* iparticle = (TParticle *) fParticles->At(i);
8d2cd130 497 Int_t ks = iparticle->GetStatusCode();
498 kf = CheckPDGCode(iparticle->GetPdgCode());
499// No initial state partons
500 if (ks==21) continue;
501//
502// Heavy Flavor Selection
503//
504 // quark ?
505 kf = TMath::Abs(kf);
506 Int_t kfl = kf;
9ff6c04c 507 // Resonance
f913ec4f 508
9ff6c04c 509 if (kfl > 100000) kfl %= 100000;
183a5ca9 510 if (kfl > 10000) kfl %= 10000;
8d2cd130 511 // meson ?
512 if (kfl > 10) kfl/=100;
513 // baryon
514 if (kfl > 10) kfl/=10;
8d2cd130 515 Int_t ipa = iparticle->GetFirstMother()-1;
516 Int_t kfMo = 0;
f913ec4f 517//
518// Establish mother daughter relation between heavy quarks and mesons
519//
520 if (kf >= fFlavorSelect && kf <= 6) {
521 Int_t idau = iparticle->GetFirstDaughter() - 1;
522 if (idau > -1) {
523 TParticle* daughter = (TParticle *) fParticles->At(idau);
524 Int_t pdgD = daughter->GetPdgCode();
525 if (pdgD == 91 || pdgD == 92) {
526 Int_t jmin = daughter->GetFirstDaughter() - 1;
527 Int_t jmax = daughter->GetLastDaughter() - 1;
528 for (Int_t j = jmin; j <= jmax; j++)
529 ((TParticle *) fParticles->At(j))->SetFirstMother(i+1);
530 } // is string or cluster
531 } // has daughter
532 } // heavy quark
8d2cd130 533
f913ec4f 534
8d2cd130 535 if (ipa > -1) {
536 TParticle * mother = (TParticle *) fParticles->At(ipa);
537 kfMo = TMath::Abs(mother->GetPdgCode());
538 }
f913ec4f 539
8d2cd130 540 // What to keep in Stack?
541 Bool_t flavorOK = kFALSE;
542 Bool_t selectOK = kFALSE;
543 if (fFeedDownOpt) {
32d6ef7d 544 if (kfl >= fFlavorSelect) flavorOK = kTRUE;
8d2cd130 545 } else {
32d6ef7d 546 if (kfl > fFlavorSelect) {
547 nc = -1;
548 break;
549 }
550 if (kfl == fFlavorSelect) flavorOK = kTRUE;
8d2cd130 551 }
552 switch (fStackFillOpt) {
553 case kFlavorSelection:
32d6ef7d 554 selectOK = kTRUE;
555 break;
8d2cd130 556 case kParentSelection:
32d6ef7d 557 if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
558 break;
8d2cd130 559 }
560 if (flavorOK && selectOK) {
561//
562// Heavy flavor hadron or quark
563//
564// Kinematic seletion on final state heavy flavor mesons
565 if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
566 {
9ff6c04c 567 continue;
8d2cd130 568 }
569 pSelected[i] = 1;
570 if (ParentSelected(kf)) ++nParents; // Update parent count
571// printf("\n particle (HF) %d %d %d", i, pSelected[i], kf);
572 } else {
573// Kinematic seletion on decay products
574 if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf)
9ff6c04c 575 && !KinematicSelection(iparticle, 1))
8d2cd130 576 {
9ff6c04c 577 continue;
8d2cd130 578 }
579//
580// Decay products
581// Select if mother was selected and is not tracked
582
583 if (pSelected[ipa] &&
584 !trackIt[ipa] && // mother will be tracked ?
585 kfMo != 5 && // mother is b-quark, don't store fragments
586 kfMo != 4 && // mother is c-quark, don't store fragments
587 kf != 92) // don't store string
588 {
589//
590// Semi-stable or de-selected: diselect decay products:
591//
592//
593 if (pSelected[i] == -1 || fDecayer->GetLifetime(kf) > fMaxLifeTime)
594 {
595 Int_t ipF = iparticle->GetFirstDaughter();
596 Int_t ipL = iparticle->GetLastDaughter();
597 if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
598 }
599// printf("\n particle (decay) %d %d %d", i, pSelected[i], kf);
600 pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
601 }
602 }
603 if (pSelected[i] == -1) pSelected[i] = 0;
604 if (!pSelected[i]) continue;
605 // Count quarks only if you did not include fragmentation
606 if (fFragmentation && kf <= 10) continue;
9ff6c04c 607
8d2cd130 608 nc++;
609// Decision on tracking
610 trackIt[i] = 0;
611//
612// Track final state particle
613 if (ks == 1) trackIt[i] = 1;
614// Track semi-stable particles
d25cfd65 615 if ((ks == 1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime)) trackIt[i] = 1;
8d2cd130 616// Track particles selected by process if undecayed.
617 if (fForceDecay == kNoDecay) {
618 if (ParentSelected(kf)) trackIt[i] = 1;
619 } else {
620 if (ParentSelected(kf)) trackIt[i] = 0;
621 }
622 if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
623//
624//
625
626 } // particle selection loop
627 if (nc > 0) {
628 for (i = 0; i<np; i++) {
629 if (!pSelected[i]) continue;
630 TParticle * iparticle = (TParticle *) fParticles->At(i);
631 kf = CheckPDGCode(iparticle->GetPdgCode());
632 Int_t ks = iparticle->GetStatusCode();
633 p[0] = iparticle->Px();
634 p[1] = iparticle->Py();
635 p[2] = iparticle->Pz();
a920faf9 636 p[3] = iparticle->Energy();
637
2590ccf9 638 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
639 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
640 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
641
8d2cd130 642 Float_t tof = kconv*iparticle->T();
643 Int_t ipa = iparticle->GetFirstMother()-1;
644 Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
a920faf9 645
646 PushTrack(fTrackIt*trackIt[i], iparent, kf,
647 p[0], p[1], p[2], p[3],
648 origin[0], origin[1], origin[2], tof,
649 polar[0], polar[1], polar[2],
650 kPPrimary, nt, 1., ks);
8d2cd130 651 pParent[i] = nt;
652 KeepTrack(nt);
642f15cf 653 } // PushTrack loop
8d2cd130 654 }
655 } else {
656 nc = GenerateMB();
657 } // mb ?
f913ec4f 658
659 GetSubEventTime();
8d2cd130 660
661 if (pParent) delete[] pParent;
662 if (pSelected) delete[] pSelected;
663 if (trackIt) delete[] trackIt;
664
665 if (nc > 0) {
666 switch (fCountMode) {
667 case kCountAll:
668 // printf(" Count all \n");
669 jev += nc;
670 break;
671 case kCountParents:
672 // printf(" Count parents \n");
673 jev += nParents;
674 break;
675 case kCountTrackables:
676 // printf(" Count trackable \n");
677 jev += nTkbles;
678 break;
679 }
680 if (jev >= fNpart || fNpart == -1) {
681 fKineBias=Float_t(fNpart)/Float_t(fTrials);
e0e89f40 682
8d2cd130 683 fQ += fPythia->GetVINT(51);
684 fX1 += fPythia->GetVINT(41);
685 fX2 += fPythia->GetVINT(42);
686 fTrialsRun += fTrials;
687 fNev++;
688 MakeHeader();
689 break;
690 }
691 }
692 } // event loop
693 SetHighWaterMark(nt);
694// adjust weight due to kinematic selection
b88f5cea 695// AdjustWeights();
8d2cd130 696// get cross-section
697 fXsection=fPythia->GetPARI(1);
698}
699
700Int_t AliGenPythia::GenerateMB()
701{
702//
703// Min Bias selection and other global selections
704//
705 Int_t i, kf, nt, iparent;
706 Int_t nc = 0;
bf950da8 707 Float_t p[4];
8d2cd130 708 Float_t polar[3] = {0,0,0};
709 Float_t origin[3] = {0,0,0};
710// converts from mm/c to s
711 const Float_t kconv=0.001/2.999792458e8;
712
e0e89f40 713
714
5fa4b20b 715 Int_t np = (fHadronisation) ? fParticles->GetEntriesFast() : fNpartons;
aea21c57 716
5fa4b20b 717
e0e89f40 718
8d2cd130 719 Int_t* pParent = new Int_t[np];
720 for (i=0; i< np; i++) pParent[i] = -1;
721 if (fProcess == kPyJets || fProcess == kPyDirectGamma) {
722 TParticle* jet1 = (TParticle *) fParticles->At(6);
723 TParticle* jet2 = (TParticle *) fParticles->At(7);
aea21c57 724 if (!CheckTrigger(jet1, jet2)) return 0;
8d2cd130 725 }
e0e89f40 726
7ce1321b 727 if (fTriggerParticle) {
728 Bool_t triggered = kFALSE;
729 for (i = 0; i < np; i++) {
730 TParticle * iparticle = (TParticle *) fParticles->At(i);
731 kf = CheckPDGCode(iparticle->GetPdgCode());
732 if (TMath::Abs(kf) != fTriggerParticle) continue;
733 if (iparticle->Pt() == 0.) continue;
734 if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
735 triggered = kTRUE;
736 break;
737 }
738 if (!triggered) return 0;
739 }
e0e89f40 740
741
742 // Check if there is a ccbar or bbbar pair with at least one of the two
743 // in fYMin < y < fYMax
744 if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi) {
745 TParticle *hvq;
746 Bool_t theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
747 Float_t yQ;
748 Int_t pdgQ;
749 for(i=0; i<np; i++) {
750 hvq = (TParticle*)fParticles->At(i);
751 pdgQ = hvq->GetPdgCode();
752 if(TMath::Abs(pdgQ) != fFlavorSelect) continue;
753 if(pdgQ>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
754 yQ = 0.5*TMath::Log((hvq->Energy()+hvq->Pz()+1.e-13)/
755 (hvq->Energy()-hvq->Pz()+1.e-13));
756 if(yQ>fYMin && yQ<fYMax) inYcut=kTRUE;
757 }
758 if (!theQ || !theQbar || !inYcut) {
759 if (pParent) delete[] pParent;
760 return 0;
761 }
762 }
aea21c57 763
0f6ee828 764 //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
765 if ( (fProcess == kPyW || fProcess == kPyZ || fProcess == kPyMb || fProcess == kPyMbNonDiffr)
766 && (fCutOnChild == 1) ) {
767 if ( !CheckKinematicsOnChild() ) {
768 if (pParent) delete[] pParent;
769 return 0;
770 }
aea21c57 771 }
772
773
f913ec4f 774 for (i = 0; i < np; i++) {
8d2cd130 775 Int_t trackIt = 0;
776 TParticle * iparticle = (TParticle *) fParticles->At(i);
777 kf = CheckPDGCode(iparticle->GetPdgCode());
778 Int_t ks = iparticle->GetStatusCode();
779 Int_t km = iparticle->GetFirstMother();
780 if ((ks == 1 && kf!=0 && KinematicSelection(iparticle, 0)) ||
781 (ks != 1) ||
782 (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
783 nc++;
784 if (ks == 1) trackIt = 1;
785 Int_t ipa = iparticle->GetFirstMother()-1;
786
787 iparent = (ipa > -1) ? pParent[ipa] : -1;
788
789//
790// store track information
791 p[0] = iparticle->Px();
792 p[1] = iparticle->Py();
793 p[2] = iparticle->Pz();
a920faf9 794 p[3] = iparticle->Energy();
1406f599 795
a920faf9 796
2590ccf9 797 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
798 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
799 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
800
f913ec4f 801 Float_t tof = fEventTime + kconv * iparticle->T();
a920faf9 802
803 PushTrack(fTrackIt*trackIt, iparent, kf,
804 p[0], p[1], p[2], p[3],
805 origin[0], origin[1], origin[2], tof,
806 polar[0], polar[1], polar[2],
807 kPPrimary, nt, 1., ks);
5fa4b20b 808 //
809 // Special Treatment to store color-flow
810 //
811 if (ks == 3 || ks == 13 || ks == 14) {
812 TParticle* particle = 0;
813 if (fStack) {
814 particle = fStack->Particle(nt);
815 } else {
816 particle = gAlice->Stack()->Particle(nt);
817 }
818 particle->SetFirstDaughter(fPythia->GetK(2, i));
819 particle->SetLastDaughter(fPythia->GetK(3, i));
820 }
821
8d2cd130 822 KeepTrack(nt);
823 pParent[i] = nt;
f913ec4f 824 SetHighWaterMark(nt);
825
8d2cd130 826 } // select particle
827 } // particle loop
828
829 if (pParent) delete[] pParent;
e0e89f40 830
f913ec4f 831 return 1;
8d2cd130 832}
833
834
835void AliGenPythia::FinishRun()
836{
837// Print x-section summary
838 fPythia->Pystat(1);
2779fc64 839
840 if (fNev > 0.) {
841 fQ /= fNev;
842 fX1 /= fNev;
843 fX2 /= fNev;
844 }
845
8d2cd130 846 printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
847 printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
8d2cd130 848}
849
850void AliGenPythia::AdjustWeights()
851{
852// Adjust the weights after generation of all events
853//
e2bddf81 854 if (gAlice) {
855 TParticle *part;
856 Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
857 for (Int_t i=0; i<ntrack; i++) {
858 part= gAlice->GetMCApp()->Particle(i);
859 part->SetWeight(part->GetWeight()*fKineBias);
860 }
8d2cd130 861 }
862}
863
864void AliGenPythia::SetNuclei(Int_t a1, Int_t a2)
865{
866// Treat protons as inside nuclei with mass numbers a1 and a2
1d568bc2 867
1a626d4e 868 fAProjectile = a1;
869 fATarget = a2;
1d568bc2 870 fSetNuclei = kTRUE;
8d2cd130 871}
872
873
874void AliGenPythia::MakeHeader()
875{
183a5ca9 876 if (gAlice) {
877 if (gAlice->GetEvNumber()>=fDebugEventFirst &&
f913ec4f 878 gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(2);
183a5ca9 879 }
880
8d2cd130 881// Builds the event header, to be called after each event
e5c87a3d 882 if (fHeader) delete fHeader;
883 fHeader = new AliGenPythiaEventHeader("Pythia");
8d2cd130 884//
885// Event type
e5c87a3d 886 ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
8d2cd130 887//
888// Number of trials
e5c87a3d 889 ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
8d2cd130 890//
891// Event Vertex
d25cfd65 892 fHeader->SetPrimaryVertex(fVertex);
8d2cd130 893//
894// Jets that have triggered
f913ec4f 895
8d2cd130 896 if (fProcess == kPyJets)
897 {
898 Int_t ntrig, njet;
899 Float_t jets[4][10];
900 GetJets(njet, ntrig, jets);
9ff6c04c 901
8d2cd130 902
903 for (Int_t i = 0; i < ntrig; i++) {
e5c87a3d 904 ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
8d2cd130 905 jets[3][i]);
906 }
907 }
5fa4b20b 908//
909// Copy relevant information from external header, if present.
910//
911 Float_t uqJet[4];
912
913 if (fRL) {
914 AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
915 for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
916 {
917 printf("Adding Jet %d %d \n", i, exHeader->NTriggerJets());
918
919
920 exHeader->TriggerJet(i, uqJet);
921 ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
922 }
923 }
924//
925// Store quenching parameters
926//
927 if (fQuench){
928 Double_t z[4];
929 Double_t xp, yp;
7c21f297 930 if (fQuench == 1) {
931 // Pythia::Quench()
932 fPythia->GetQuenchingParameters(xp, yp, z);
933 } else {
934 // Pyquen
935 Double_t r1 = PARIMP.rb1;
936 Double_t r2 = PARIMP.rb2;
937 Double_t b = PARIMP.b1;
938 Double_t r = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
939 Double_t phi = PARIMP.psib1;
940 xp = r * TMath::Cos(phi);
941 yp = r * TMath::Sin(phi);
942
943 }
944 ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
945 ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
946 }
beac474c 947//
948// Store pt^hard
949 ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
5fa4b20b 950//
cf57b268 951// Pass header
5fa4b20b 952//
cf57b268 953 AddHeader(fHeader);
8d2cd130 954}
cf57b268 955
956void AliGenPythia::AddHeader(AliGenEventHeader* header)
957{
958 // Add header to container or runloader
959 if (fContainer) {
960 fContainer->AddHeader(header);
961 } else {
962 AliRunLoader::GetRunLoader()->GetHeader()->SetGenEventHeader(header);
963 }
964}
965
8d2cd130 966
967Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
968{
969// Check the kinematic trigger condition
970//
971 Double_t eta[2];
972 eta[0] = jet1->Eta();
973 eta[1] = jet2->Eta();
974 Double_t phi[2];
975 phi[0] = jet1->Phi();
976 phi[1] = jet2->Phi();
977 Int_t pdg[2];
978 pdg[0] = jet1->GetPdgCode();
979 pdg[1] = jet2->GetPdgCode();
980 Bool_t triggered = kFALSE;
981
982 if (fProcess == kPyJets) {
983 Int_t njets = 0;
984 Int_t ntrig = 0;
985 Float_t jets[4][10];
986//
987// Use Pythia clustering on parton level to determine jet axis
988//
989 GetJets(njets, ntrig, jets);
990
991 if (ntrig) triggered = kTRUE;
992//
993 } else {
994 Int_t ij = 0;
995 Int_t ig = 1;
996 if (pdg[0] == kGamma) {
997 ij = 1;
998 ig = 0;
999 }
1000 //Check eta range first...
1001 if ((eta[ij] < fEtaMaxJet && eta[ij] > fEtaMinJet) &&
1002 (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
1003 {
1004 //Eta is okay, now check phi range
1005 if ((phi[ij] < fPhiMaxJet && phi[ij] > fPhiMinJet) &&
1006 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
1007 {
1008 triggered = kTRUE;
1009 }
1010 }
1011 }
1012 return triggered;
1013}
aea21c57 1014
1015
1016//Checking Kinematics on Child (status code 1, particle code ?, kin cuts
1017Bool_t AliGenPythia::CheckKinematicsOnChild(){
1018
1019 Bool_t checking = kFALSE;
1020 Int_t j, kcode, ks, km;
1021 Int_t nPartAcc = 0; //number of particles in the acceptance range
1022 Int_t numberOfAcceptedParticles = 1;
1023 if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
1024 Int_t npart = fParticles->GetEntriesFast();
1025
0f6ee828 1026 for (j = 0; j<npart; j++) {
aea21c57 1027 TParticle * jparticle = (TParticle *) fParticles->At(j);
1028 kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
1029 ks = jparticle->GetStatusCode();
1030 km = jparticle->GetFirstMother();
1031
1032 if( (ks == 1) && (kcode == fPdgCodeParticleforAcceptanceCut) && (KinematicSelection(jparticle,1)) ){
1033 nPartAcc++;
1034 }
0f6ee828 1035 if( numberOfAcceptedParticles <= nPartAcc){
1036 checking = kTRUE;
1037 break;
1038 }
aea21c57 1039 }
0f6ee828 1040
aea21c57 1041 return checking;
aea21c57 1042}
1043
8d2cd130 1044
1045AliGenPythia& AliGenPythia::operator=(const AliGenPythia& rhs)
1046{
1047// Assignment operator
014a9521 1048 rhs.Copy(*this);
8d2cd130 1049 return *this;
1050}
1051
5fa4b20b 1052void AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
8d2cd130 1053{
1054//
1055// Load event into Pythia Common Block
1056//
5fa4b20b 1057
32d6ef7d 1058 Int_t npart = stack -> GetNprimary();
1059 Int_t n0 = 0;
1060
1061 if (!flag) {
1062 (fPythia->GetPyjets())->N = npart;
1063 } else {
1064 n0 = (fPythia->GetPyjets())->N;
1065 (fPythia->GetPyjets())->N = n0 + npart;
1066 }
1067
1068
8d2cd130 1069 for (Int_t part = 0; part < npart; part++) {
32d6ef7d 1070 TParticle *MPart = stack->Particle(part);
1071
5fa4b20b 1072 Int_t kf = MPart->GetPdgCode();
1073 Int_t ks = MPart->GetStatusCode();
1074 Int_t idf = MPart->GetFirstDaughter();
1075 Int_t idl = MPart->GetLastDaughter();
1076
1077 if (reHadr) {
1078 if (ks == 11 || ks == 12) {
1079 ks -= 10;
1080 idf = -1;
1081 idl = -1;
1082 }
1083 }
32d6ef7d 1084
8d2cd130 1085 Float_t px = MPart->Px();
1086 Float_t py = MPart->Py();
1087 Float_t pz = MPart->Pz();
1088 Float_t e = MPart->Energy();
a920faf9 1089 Float_t m = MPart->GetCalcMass();
8d2cd130 1090
1091
32d6ef7d 1092 (fPythia->GetPyjets())->P[0][part+n0] = px;
1093 (fPythia->GetPyjets())->P[1][part+n0] = py;
1094 (fPythia->GetPyjets())->P[2][part+n0] = pz;
1095 (fPythia->GetPyjets())->P[3][part+n0] = e;
1096 (fPythia->GetPyjets())->P[4][part+n0] = m;
8d2cd130 1097
32d6ef7d 1098 (fPythia->GetPyjets())->K[1][part+n0] = kf;
1099 (fPythia->GetPyjets())->K[0][part+n0] = ks;
5fa4b20b 1100 (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1101 (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1102 (fPythia->GetPyjets())->K[2][part+n0] = MPart->GetFirstMother() + 1;
8d2cd130 1103 }
1104}
1105
5fa4b20b 1106
014a9521 1107void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
8d2cd130 1108{
1109//
1110// Calls the Pythia jet finding algorithm to find jets in the current event
1111//
1112//
8d2cd130 1113//
1114// Save jets
1115 Int_t n = fPythia->GetN();
1116
1117//
1118// Run Jet Finder
1119 fPythia->Pycell(njets);
1120 Int_t i;
1121 for (i = 0; i < njets; i++) {
1122 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
1123 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
1124 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
1125 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
1126
1127 jets[0][i] = px;
1128 jets[1][i] = py;
1129 jets[2][i] = pz;
1130 jets[3][i] = e;
1131 }
1132}
1133
1134
1135
1136void AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
1137{
1138//
1139// Calls the Pythia clustering algorithm to find jets in the current event
1140//
1141 Int_t n = fPythia->GetN();
1142 nJets = 0;
1143 nJetsTrig = 0;
1144 if (fJetReconstruction == kCluster) {
1145//
1146// Configure cluster algorithm
1147//
1148 fPythia->SetPARU(43, 2.);
1149 fPythia->SetMSTU(41, 1);
1150//
1151// Call cluster algorithm
1152//
1153 fPythia->Pyclus(nJets);
1154//
1155// Loading jets from common block
1156//
1157 } else {
592f8307 1158
8d2cd130 1159//
1160// Run Jet Finder
1161 fPythia->Pycell(nJets);
1162 }
1163
1164 Int_t i;
1165 for (i = 0; i < nJets; i++) {
1166 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
1167 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
1168 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
1169 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
1170 Float_t pt = TMath::Sqrt(px * px + py * py);
a920faf9 1171 Float_t phi = TMath::Pi() + TMath::ATan2(-py, -px);
8d2cd130 1172 Float_t theta = TMath::ATan2(pt,pz);
1173 Float_t et = e * TMath::Sin(theta);
1174 Float_t eta = -TMath::Log(TMath::Tan(theta / 2.));
8d2cd130 1175 if (
1176 eta > fEtaMinJet && eta < fEtaMaxJet &&
675eb105 1177 phi > fPhiMinJet && phi < fPhiMaxJet &&
8d2cd130 1178 et > fEtMinJet && et < fEtMaxJet
1179 )
1180 {
1181 jets[0][nJetsTrig] = px;
1182 jets[1][nJetsTrig] = py;
1183 jets[2][nJetsTrig] = pz;
1184 jets[3][nJetsTrig] = e;
1185 nJetsTrig++;
5fa4b20b 1186// printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
8d2cd130 1187 } else {
1188// printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1189 }
1190 }
1191}
1192
f913ec4f 1193void AliGenPythia::GetSubEventTime()
1194{
1195 // Calculates time of the next subevent
9d7108a7 1196 fEventTime = 0.;
1197 if (fEventsTime) {
1198 TArrayF &array = *fEventsTime;
1199 fEventTime = array[fCurSubEvent++];
1200 }
1201 // printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
1202 return;
f913ec4f 1203}
8d2cd130 1204
1205#ifdef never
1206void AliGenPythia::Streamer(TBuffer &R__b)
1207{
1208 // Stream an object of class AliGenPythia.
1209
1210 if (R__b.IsReading()) {
1211 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1212 AliGenerator::Streamer(R__b);
1213 R__b >> (Int_t&)fProcess;
1214 R__b >> (Int_t&)fStrucFunc;
1215 R__b >> (Int_t&)fForceDecay;
1216 R__b >> fEnergyCMS;
1217 R__b >> fKineBias;
1218 R__b >> fTrials;
1219 fParentSelect.Streamer(R__b);
1220 fChildSelect.Streamer(R__b);
1221 R__b >> fXsection;
1222// (AliPythia::Instance())->Streamer(R__b);
1223 R__b >> fPtHardMin;
1224 R__b >> fPtHardMax;
1225// if (fDecayer) fDecayer->Streamer(R__b);
1226 } else {
1227 R__b.WriteVersion(AliGenPythia::IsA());
1228 AliGenerator::Streamer(R__b);
1229 R__b << (Int_t)fProcess;
1230 R__b << (Int_t)fStrucFunc;
1231 R__b << (Int_t)fForceDecay;
1232 R__b << fEnergyCMS;
1233 R__b << fKineBias;
1234 R__b << fTrials;
1235 fParentSelect.Streamer(R__b);
1236 fChildSelect.Streamer(R__b);
1237 R__b << fXsection;
1238// R__b << fPythia;
1239 R__b << fPtHardMin;
1240 R__b << fPtHardMax;
1241 // fDecayer->Streamer(R__b);
1242 }
1243}
1244#endif
1245
90d7b703 1246
589380c6 1247