1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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 ...)
25 // andreas.morsch@cern.ch
28 #include <TDatabasePDG.h>
29 #include <TParticle.h>
35 #include "AliDecayerPythia.h"
36 #include "AliGenPythia.h"
37 #include "AliGenPythiaEventHeader.h"
38 #include "AliPythia.h"
39 #include "AliPythiaRndm.h"
42 #include "AliRunLoader.h"
44 ClassImp(AliGenPythia)
46 AliGenPythia::AliGenPythia()
49 // Default Constructor
53 fDecayer = new AliDecayerPythia();
64 if (!AliPythiaRndm::GetPythiaRandom())
65 AliPythiaRndm::SetPythiaRandom(GetRandom());
68 AliGenPythia::AliGenPythia(Int_t npart)
71 // default charm production at 5. 5 TeV
73 // structure function GRVHO
76 fTitle= "Particle Generator using PYTHIA";
84 fDecayer = new AliDecayerPythia();
85 // Set random number generator
86 if (!AliPythiaRndm::GetPythiaRandom())
87 AliPythiaRndm::SetPythiaRandom(GetRandom());
90 fParticles = new TClonesArray("TParticle",1000);
98 SetJetReconstructionMode();
101 // Options determining what to keep in the stack (Heavy flavour generation)
102 fStackFillOpt = kFlavorSelection; // Keep particle with selected flavor
103 fFeedDownOpt = kTRUE; // allow feed down from higher family
104 // Fragmentation on/off
105 fFragmentation = kTRUE;
106 // Default counting mode
107 fCountMode = kCountAll;
109 SetPycellParameters();
113 AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)
120 AliGenPythia::~AliGenPythia()
125 void AliGenPythia::SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
126 Float_t thresh, Float_t etseed, Float_t minet, Float_t r)
128 // Set pycell parameters
129 fPycellEtaMax = etamax;
132 fPycellThreshold = thresh;
133 fPycellEtSeed = etseed;
134 fPycellMinEtJet = minet;
135 fPycellMaxRadius = r;
140 void AliGenPythia::SetEventListRange(Int_t eventFirst, Int_t eventLast)
142 // Set a range of event numbers, for which a table
143 // of generated particle will be printed
144 fDebugEventFirst = eventFirst;
145 fDebugEventLast = eventLast;
146 if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst;
149 void AliGenPythia::Init()
153 SetMC(AliPythia::Instance());
154 fPythia=(AliPythia*) fMCEvGen;
156 fParentWeight=1./Float_t(fNpart);
158 // Forward Paramters to the AliPythia object
159 fDecayer->SetForceDecay(fForceDecay);
163 fPythia->SetCKIN(3,fPtHardMin);
164 fPythia->SetCKIN(4,fPtHardMax);
165 fPythia->SetCKIN(7,fYHardMin);
166 fPythia->SetCKIN(8,fYHardMax);
168 if (fAProjectile > 0 && fATarget > 0) fPythia->SetNuclei(fAProjectile, fATarget);
170 if (fFragmentation) {
171 fPythia->SetMSTP(111,1);
173 fPythia->SetMSTP(111,0);
177 // initial state radiation
178 fPythia->SetMSTP(61,fGinit);
179 // final state radiation
180 fPythia->SetMSTP(71,fGfinal);
183 fPythia->SetMSTP(91,1);
184 fPythia->SetPARP(91,fPtKick);
186 fPythia->SetMSTP(91,0);
189 // fPythia->SetMSTJ(1,2);
191 fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
193 // Parent and Children Selection
197 case kPyCharmUnforced:
198 case kPyCharmPbPbMNR:
201 fParentSelect[0] = 411;
202 fParentSelect[1] = 421;
203 fParentSelect[2] = 431;
204 fParentSelect[3] = 4122;
210 fParentSelect[0] = 421;
214 case kPyBeautyPbPbMNR:
215 case kPyBeautypPbMNR:
217 fParentSelect[0]= 511;
218 fParentSelect[1]= 521;
219 fParentSelect[2]= 531;
220 fParentSelect[3]= 5122;
221 fParentSelect[4]= 5132;
222 fParentSelect[5]= 5232;
223 fParentSelect[6]= 5332;
226 case kPyBeautyUnforced:
227 fParentSelect[0] = 511;
228 fParentSelect[1] = 521;
229 fParentSelect[2] = 531;
230 fParentSelect[3] = 5122;
231 fParentSelect[4] = 5132;
232 fParentSelect[5] = 5232;
233 fParentSelect[6] = 5332;
238 fParentSelect[0] = 443;
248 // JetFinder for Trigger
250 // Configure detector (EMCAL like)
252 fPythia->SetPARU(51, fPycellEtaMax);
253 fPythia->SetMSTU(51, fPycellNEta);
254 fPythia->SetMSTU(52, fPycellNPhi);
256 // Configure Jet Finder
258 fPythia->SetPARU(58, fPycellThreshold);
259 fPythia->SetPARU(52, fPycellEtSeed);
260 fPythia->SetPARU(53, fPycellMinEtJet);
261 fPythia->SetPARU(54, fPycellMaxRadius);
262 fPythia->SetMSTU(54, 2);
264 // This counts the total number of calls to Pyevnt() per run.
279 Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n");
283 void AliGenPythia::Generate()
285 // Generate one event
287 fDecayer->ForceDecay();
289 Float_t polar[3] = {0,0,0};
290 Float_t origin[3] = {0,0,0};
292 // converts from mm/c to s
293 const Float_t kconv=0.001/2.999792458e8;
300 // Set collision vertex position
301 if(fVertexSmear==kPerEvent) {
302 fPythia->SetMSTP(151,1);
304 fPythia->SetPARP(151+j, fOsigma[j]*10.);
306 } else if (fVertexSmear==kPerTrack) {
307 fPythia->SetMSTP(151,0);
313 fPythia->SetMSTJ(1, 0);
320 fPythia->SetMSTJ(1, 1);
325 if (gAlice->GetEvNumber()>=fDebugEventFirst &&
326 gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(1);
329 fPythia->ImportParticles(fParticles,"All");
337 Int_t np = fParticles->GetEntriesFast();
338 if (np == 0 ) continue;
339 // Get event vertex and discard the event if the Z coord. is too big
340 TParticle *iparticle = (TParticle *) fParticles->At(0);
341 Float_t distz = iparticle->Vz()/10.;
342 if(TMath::Abs(distz)>fCutVertexZ*fOsigma[2]) continue;
344 fVertex[0] = iparticle->Vx()/10.+fOrigin.At(0);
345 fVertex[1] = iparticle->Vy()/10.+fOrigin.At(1);
346 fVertex[2] = iparticle->Vz()/10.+fOrigin.At(2);
348 Int_t* pParent = new Int_t[np];
349 Int_t* pSelected = new Int_t[np];
350 Int_t* trackIt = new Int_t[np];
351 for (i=0; i< np; i++) {
357 Int_t nc = 0; // Total n. of selected particles
358 Int_t nParents = 0; // Selected parents
359 Int_t nTkbles = 0; // Trackable particles
360 if (fProcess != kPyMb && fProcess != kPyJets &&
361 fProcess != kPyDirectGamma &&
362 fProcess != kPyMbNonDiffr) {
364 for (i = 0; i<np; i++) {
365 iparticle = (TParticle *) fParticles->At(i);
366 Int_t ks = iparticle->GetStatusCode();
367 kf = CheckPDGCode(iparticle->GetPdgCode());
368 // No initial state partons
369 if (ks==21) continue;
371 // Heavy Flavor Selection
377 if (kfl > 10) kfl/=100;
379 if (kfl > 10) kfl/=10;
380 if (kfl > 10) kfl/=10;
382 Int_t ipa = iparticle->GetFirstMother()-1;
386 TParticle * mother = (TParticle *) fParticles->At(ipa);
387 kfMo = TMath::Abs(mother->GetPdgCode());
389 // What to keep in Stack?
390 Bool_t flavorOK = kFALSE;
391 Bool_t selectOK = kFALSE;
393 if (kfl >= fFlavorSelect) flavorOK = kTRUE;
395 if (kfl > fFlavorSelect) {
399 if (kfl == fFlavorSelect) flavorOK = kTRUE;
401 switch (fStackFillOpt) {
402 case kFlavorSelection:
405 case kParentSelection:
406 if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
409 if (flavorOK && selectOK) {
411 // Heavy flavor hadron or quark
413 // Kinematic seletion on final state heavy flavor mesons
414 if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
419 if (ParentSelected(kf)) ++nParents; // Update parent count
420 // printf("\n particle (HF) %d %d %d", i, pSelected[i], kf);
422 // Kinematic seletion on decay products
423 if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf)
424 && !KinematicSelection(iparticle, 1))
430 // Select if mother was selected and is not tracked
432 if (pSelected[ipa] &&
433 !trackIt[ipa] && // mother will be tracked ?
434 kfMo != 5 && // mother is b-quark, don't store fragments
435 kfMo != 4 && // mother is c-quark, don't store fragments
436 kf != 92) // don't store string
439 // Semi-stable or de-selected: diselect decay products:
442 if (pSelected[i] == -1 || fDecayer->GetLifetime(kf) > fMaxLifeTime)
444 Int_t ipF = iparticle->GetFirstDaughter();
445 Int_t ipL = iparticle->GetLastDaughter();
446 if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
448 // printf("\n particle (decay) %d %d %d", i, pSelected[i], kf);
449 pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
452 if (pSelected[i] == -1) pSelected[i] = 0;
453 if (!pSelected[i]) continue;
454 // Count quarks only if you did not include fragmentation
455 if (fFragmentation && kf <= 10) continue;
457 // Decision on tracking
460 // Track final state particle
461 if (ks == 1) trackIt[i] = 1;
462 // Track semi-stable particles
463 if ((ks == 1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime)) trackIt[i] = 1;
464 // Track particles selected by process if undecayed.
465 if (fForceDecay == kNoDecay) {
466 if (ParentSelected(kf)) trackIt[i] = 1;
468 if (ParentSelected(kf)) trackIt[i] = 0;
470 if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
474 } // particle selection loop
476 for (i = 0; i<np; i++) {
477 if (!pSelected[i]) continue;
478 TParticle * iparticle = (TParticle *) fParticles->At(i);
479 kf = CheckPDGCode(iparticle->GetPdgCode());
480 Int_t ks = iparticle->GetStatusCode();
481 p[0] = iparticle->Px();
482 p[1] = iparticle->Py();
483 p[2] = iparticle->Pz();
484 origin[0] = fOrigin[0]+iparticle->Vx()/10.;
485 origin[1] = fOrigin[1]+iparticle->Vy()/10.;
486 origin[2] = fOrigin[2]+iparticle->Vz()/10.;
487 Float_t tof = kconv*iparticle->T();
488 Int_t ipa = iparticle->GetFirstMother()-1;
489 Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
490 PushTrack(fTrackIt*trackIt[i] ,
491 iparent, kf, p, origin, polar, tof, kPPrimary, nt, 1., ks);
500 if (pParent) delete[] pParent;
501 if (pSelected) delete[] pSelected;
502 if (trackIt) delete[] trackIt;
505 switch (fCountMode) {
507 // printf(" Count all \n");
511 // printf(" Count parents \n");
514 case kCountTrackables:
515 // printf(" Count trackable \n");
519 if (jev >= fNpart || fNpart == -1) {
520 fKineBias=Float_t(fNpart)/Float_t(fTrials);
521 printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
523 fQ += fPythia->GetVINT(51);
524 fX1 += fPythia->GetVINT(41);
525 fX2 += fPythia->GetVINT(42);
526 fTrialsRun += fTrials;
533 SetHighWaterMark(nt);
534 // adjust weight due to kinematic selection
537 fXsection=fPythia->GetPARI(1);
540 Int_t AliGenPythia::GenerateMB()
543 // Min Bias selection and other global selections
545 Int_t i, kf, nt, iparent;
548 Float_t polar[3] = {0,0,0};
549 Float_t origin[3] = {0,0,0};
550 // converts from mm/c to s
551 const Float_t kconv=0.001/2.999792458e8;
553 Int_t np = fParticles->GetEntriesFast();
554 Int_t* pParent = new Int_t[np];
555 for (i=0; i< np; i++) pParent[i] = -1;
556 if (fProcess == kPyJets || fProcess == kPyDirectGamma) {
557 TParticle* jet1 = (TParticle *) fParticles->At(6);
558 TParticle* jet2 = (TParticle *) fParticles->At(7);
559 if (!CheckTrigger(jet1, jet2)) return 0;
562 for (i = 0; i<np; i++) {
564 TParticle * iparticle = (TParticle *) fParticles->At(i);
565 kf = CheckPDGCode(iparticle->GetPdgCode());
566 Int_t ks = iparticle->GetStatusCode();
567 Int_t km = iparticle->GetFirstMother();
568 if ((ks == 1 && kf!=0 && KinematicSelection(iparticle, 0)) ||
570 (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
572 if (ks == 1) trackIt = 1;
573 Int_t ipa = iparticle->GetFirstMother()-1;
575 iparent = (ipa > -1) ? pParent[ipa] : -1;
578 // store track information
579 p[0] = iparticle->Px();
580 p[1] = iparticle->Py();
581 p[2] = iparticle->Pz();
582 origin[0] = fOrigin[0]+iparticle->Vx()/10.;
583 origin[1] = fOrigin[1]+iparticle->Vy()/10.;
584 origin[2] = fOrigin[2]+iparticle->Vz()/10.;
585 Float_t tof=kconv*iparticle->T();
586 PushTrack(fTrackIt*trackIt, iparent, kf, p, origin, polar,
587 tof, kPPrimary, nt, 1., ks);
593 if (pParent) delete[] pParent;
595 printf("\n I've put %i particles on the stack \n",nc);
600 void AliGenPythia::FinishRun()
602 // Print x-section summary
607 printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
608 printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
613 void AliGenPythia::AdjustWeights()
615 // Adjust the weights after generation of all events
618 Int_t ntrack=gAlice->GetNtrack();
619 for (Int_t i=0; i<ntrack; i++) {
620 part= gAlice->Particle(i);
621 part->SetWeight(part->GetWeight()*fKineBias);
625 void AliGenPythia::SetNuclei(Int_t a1, Int_t a2)
627 // Treat protons as inside nuclei with mass numbers a1 and a2
635 void AliGenPythia::MakeHeader()
637 // Builds the event header, to be called after each event
638 if (fHeader) delete fHeader;
639 fHeader = new AliGenPythiaEventHeader("Pythia");
642 ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
645 ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
648 fHeader->SetPrimaryVertex(fVertex);
650 // Jets that have triggered
651 if (fProcess == kPyJets)
655 GetJets(njet, ntrig, jets);
657 for (Int_t i = 0; i < ntrig; i++) {
658 ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
662 gAlice->SetGenEventHeader(fHeader);
666 Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
668 // Check the kinematic trigger condition
671 eta[0] = jet1->Eta();
672 eta[1] = jet2->Eta();
674 phi[0] = jet1->Phi();
675 phi[1] = jet2->Phi();
677 pdg[0] = jet1->GetPdgCode();
678 pdg[1] = jet2->GetPdgCode();
679 Bool_t triggered = kFALSE;
681 if (fProcess == kPyJets) {
686 // Use Pythia clustering on parton level to determine jet axis
688 GetJets(njets, ntrig, jets);
690 if (ntrig) triggered = kTRUE;
695 if (pdg[0] == kGamma) {
699 //Check eta range first...
700 if ((eta[ij] < fEtaMaxJet && eta[ij] > fEtaMinJet) &&
701 (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
703 //Eta is okay, now check phi range
704 if ((phi[ij] < fPhiMaxJet && phi[ij] > fPhiMinJet) &&
705 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
714 AliGenPythia& AliGenPythia::operator=(const AliGenPythia& rhs)
716 // Assignment operator
721 void AliGenPythia::LoadEvent()
724 // Load event into Pythia Common Block
727 AliRunLoader* rl = AliRunLoader::GetRunLoader();
728 Int_t npart = (rl->Stack())-> GetNprimary();
729 (fPythia->GetPyjets())->N = npart;
731 for (Int_t part = 0; part < npart; part++) {
732 TParticle *MPart = (rl->Stack())->Particle(part);
733 Int_t kf = MPart->GetPdgCode();
734 Int_t ks = MPart->GetStatusCode();
735 Float_t px = MPart->Px();
736 Float_t py = MPart->Py();
737 Float_t pz = MPart->Pz();
738 Float_t e = MPart->Energy();
739 Float_t p = TMath::Sqrt(px * px + py * py + pz * pz);
740 Float_t m = TMath::Sqrt(e * e - p * p);
743 (fPythia->GetPyjets())->P[0][part] = px;
744 (fPythia->GetPyjets())->P[1][part] = py;
745 (fPythia->GetPyjets())->P[2][part] = pz;
746 (fPythia->GetPyjets())->P[3][part] = e;
747 (fPythia->GetPyjets())->P[4][part] = m;
749 (fPythia->GetPyjets())->K[1][part] = kf;
750 (fPythia->GetPyjets())->K[0][part] = ks;
754 void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
757 // Calls the Pythia jet finding algorithm to find jets in the current event
762 Int_t n = fPythia->GetN();
766 fPythia->Pycell(njets);
768 for (i = 0; i < njets; i++) {
769 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
770 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
771 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
772 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
783 void AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
786 // Calls the Pythia clustering algorithm to find jets in the current event
788 Int_t n = fPythia->GetN();
791 if (fJetReconstruction == kCluster) {
793 // Configure cluster algorithm
795 fPythia->SetPARU(43, 2.);
796 fPythia->SetMSTU(41, 1);
798 // Call cluster algorithm
800 fPythia->Pyclus(nJets);
802 // Loading jets from common block
808 fPythia->Pycell(nJets);
812 for (i = 0; i < nJets; i++) {
813 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
814 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
815 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
816 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
817 Float_t pt = TMath::Sqrt(px * px + py * py);
818 Float_t phi = TMath::ATan2(py,px);
819 Float_t theta = TMath::ATan2(pt,pz);
820 Float_t et = e * TMath::Sin(theta);
821 Float_t eta = -TMath::Log(TMath::Tan(theta / 2.));
824 eta > fEtaMinJet && eta < fEtaMaxJet &&
825 phi > fPhiMinJet && phi < fPhiMaxJet &&
826 et > fEtMinJet && et < fEtMaxJet
829 jets[0][nJetsTrig] = px;
830 jets[1][nJetsTrig] = py;
831 jets[2][nJetsTrig] = pz;
832 jets[3][nJetsTrig] = e;
836 // printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
843 void AliGenPythia::Streamer(TBuffer &R__b)
845 // Stream an object of class AliGenPythia.
847 if (R__b.IsReading()) {
848 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
849 AliGenerator::Streamer(R__b);
850 R__b >> (Int_t&)fProcess;
851 R__b >> (Int_t&)fStrucFunc;
852 R__b >> (Int_t&)fForceDecay;
856 fParentSelect.Streamer(R__b);
857 fChildSelect.Streamer(R__b);
859 // (AliPythia::Instance())->Streamer(R__b);
862 // if (fDecayer) fDecayer->Streamer(R__b);
864 R__b.WriteVersion(AliGenPythia::IsA());
865 AliGenerator::Streamer(R__b);
866 R__b << (Int_t)fProcess;
867 R__b << (Int_t)fStrucFunc;
868 R__b << (Int_t)fForceDecay;
872 fParentSelect.Streamer(R__b);
873 fChildSelect.Streamer(R__b);
878 // fDecayer->Streamer(R__b);