Generator for Pythia6 and Pythia8 event generation.
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+/* $Id$ */
+
+//
+// Generator using the TPythia interface (via AliPythia)
+// to generate pp collisions.
+// Using SetNuclei() also nuclear modifications to the structure functions
+// can be taken into account. This makes, of course, only sense for the
+// generation of the products of hard processes (heavy flavor, jets ...)
+//
+// andreas.morsch@cern.ch
+//
+
+#include <TClonesArray.h>
+#include <TDatabasePDG.h>
+#include <TParticle.h>
+#include <TPDGCode.h>
+#include <TSystem.h>
+#include <TTree.h>
+#include "AliConst.h"
+#include "AliDecayerPythia.h"
+#include "AliGenPythiaPlus.h"
+#include "AliHeader.h"
+#include "AliGenPythiaEventHeader.h"
+#include "AliPythiaBase.h"
+#include "AliPythiaRndm.h"
+#include "AliRun.h"
+#include "AliStack.h"
+#include "AliRunLoader.h"
+#include "AliMC.h"
+#include "PyquenCommon.h"
+
+ClassImp(AliGenPythiaPlus)
+
+
+AliGenPythiaPlus::AliGenPythiaPlus():
+ AliGenMC(),
+ fPythia(0),
+ fProcess(kPyCharm),
+ fStrucFunc(kCTEQ5L),
+ fEnergyCMS(5500.),
+ fKineBias(0.),
+ fTrials(0),
+ fTrialsRun(0),
+ fQ(0.),
+ fX1(0.),
+ fX2(0.),
+ fEventTime(0.),
+ fInteractionRate(0.),
+ fTimeWindow(0.),
+ fCurSubEvent(0),
+ fEventsTime(0),
+ fNev(0),
+ fFlavorSelect(0),
+ fXsection(0.),
+ fPtHardMin(0.),
+ fPtHardMax(1.e4),
+ fYHardMin(-1.e10),
+ fYHardMax(1.e10),
+ fGinit(1),
+ fGfinal(1),
+ fHadronisation(1),
+ fNpartons(0),
+ fReadFromFile(0),
+ fQuench(0),
+ fPtKick(1.),
+ fFullEvent(kTRUE),
+ fDecayer(new AliDecayerPythia()),
+ fDebugEventFirst(-1),
+ fDebugEventLast(-1),
+ fEtMinJet(0.),
+ fEtMaxJet(1.e4),
+ fEtaMinJet(-20.),
+ fEtaMaxJet(20.),
+ fPhiMinJet(0.),
+ fPhiMaxJet(2.* TMath::Pi()),
+ fJetReconstruction(kCell),
+ fEtaMinGamma(-20.),
+ fEtaMaxGamma(20.),
+ fPhiMinGamma(0.),
+ fPhiMaxGamma(2. * TMath::Pi()),
+ fPycellEtaMax(2.),
+ fPycellNEta(274),
+ fPycellNPhi(432),
+ fPycellThreshold(0.),
+ fPycellEtSeed(4.),
+ fPycellMinEtJet(10.),
+ fPycellMaxRadius(1.),
+ fStackFillOpt(kFlavorSelection),
+ fFeedDownOpt(kTRUE),
+ fFragmentation(kTRUE),
+ fSetNuclei(kFALSE),
+ fNewMIS(kFALSE),
+ fHFoff(kFALSE),
+ fTriggerParticle(0),
+ fTriggerEta(0.9),
+ fCountMode(kCountAll),
+ fHeader(0),
+ fRL(0),
+ fFileName(0),
+ fFragPhotonInCalo(kFALSE),
+ fPi0InCalo(kFALSE) ,
+ fPhotonInCalo(kFALSE),
+ fCheckEMCAL(kFALSE),
+ fCheckPHOS(kFALSE),
+ fCheckPHOSeta(kFALSE),
+ fFragPhotonOrPi0MinPt(0),
+ fPhotonMinPt(0),
+ fPHOSMinPhi(219.),
+ fPHOSMaxPhi(321.),
+ fPHOSEta(0.13),
+ fEMCALMinPhi(79.),
+ fEMCALMaxPhi(191.),
+ fEMCALEta(0.71)
+
+{
+// Default Constructor
+ SetNuclei(0,0);
+ if (!AliPythiaRndm::GetPythiaRandom())
+ AliPythiaRndm::SetPythiaRandom(GetRandom());
+}
+
+AliGenPythiaPlus::AliGenPythiaPlus(AliPythiaBase* pythia)
+ :AliGenMC(-1),
+ fPythia(pythia),
+ fProcess(kPyCharm),
+ fStrucFunc(kCTEQ5L),
+ fEnergyCMS(5500.),
+ fKineBias(0.),
+ fTrials(0),
+ fTrialsRun(0),
+ fQ(0.),
+ fX1(0.),
+ fX2(0.),
+ fEventTime(0.),
+ fInteractionRate(0.),
+ fTimeWindow(0.),
+ fCurSubEvent(0),
+ fEventsTime(0),
+ fNev(0),
+ fFlavorSelect(0),
+ fXsection(0.),
+ fPtHardMin(0.),
+ fPtHardMax(1.e4),
+ fYHardMin(-1.e10),
+ fYHardMax(1.e10),
+ fGinit(kTRUE),
+ fGfinal(kTRUE),
+ fHadronisation(kTRUE),
+ fNpartons(0),
+ fReadFromFile(kFALSE),
+ fQuench(kFALSE),
+ fPtKick(1.),
+ fFullEvent(kTRUE),
+ fDecayer(new AliDecayerPythia()),
+ fDebugEventFirst(-1),
+ fDebugEventLast(-1),
+ fEtMinJet(0.),
+ fEtMaxJet(1.e4),
+ fEtaMinJet(-20.),
+ fEtaMaxJet(20.),
+ fPhiMinJet(0.),
+ fPhiMaxJet(2.* TMath::Pi()),
+ fJetReconstruction(kCell),
+ fEtaMinGamma(-20.),
+ fEtaMaxGamma(20.),
+ fPhiMinGamma(0.),
+ fPhiMaxGamma(2. * TMath::Pi()),
+ fPycellEtaMax(2.),
+ fPycellNEta(274),
+ fPycellNPhi(432),
+ fPycellThreshold(0.),
+ fPycellEtSeed(4.),
+ fPycellMinEtJet(10.),
+ fPycellMaxRadius(1.),
+ fStackFillOpt(kFlavorSelection),
+ fFeedDownOpt(kTRUE),
+ fFragmentation(kTRUE),
+ fSetNuclei(kFALSE),
+ fNewMIS(kFALSE),
+ fHFoff(kFALSE),
+ fTriggerParticle(0),
+ fTriggerEta(0.9),
+ fCountMode(kCountAll),
+ fHeader(0),
+ fRL(0),
+ fFileName(0),
+ fFragPhotonInCalo(kFALSE),
+ fPi0InCalo(kFALSE) ,
+ fPhotonInCalo(kFALSE),
+ fCheckEMCAL(kFALSE),
+ fCheckPHOS(kFALSE),
+ fCheckPHOSeta(kFALSE),
+ fFragPhotonOrPi0MinPt(0),
+ fPhotonMinPt(0),
+ fPHOSMinPhi(219.),
+ fPHOSMaxPhi(321.),
+ fPHOSEta(0.13),
+ fEMCALMinPhi(79.),
+ fEMCALMaxPhi(191.),
+ fEMCALEta(0.71)
+{
+// default charm production at 5. 5 TeV
+// semimuonic decay
+// structure function GRVHO
+//
+ fName = "Pythia";
+ fTitle= "Particle Generator using PYTHIA";
+ SetForceDecay();
+ // Set random number generator
+ if (!AliPythiaRndm::GetPythiaRandom())
+ AliPythiaRndm::SetPythiaRandom(GetRandom());
+ fParticles = new TClonesArray("TParticle",1000);
+ SetNuclei(0,0);
+ }
+
+AliGenPythiaPlus::~AliGenPythiaPlus()
+{
+// Destructor
+ if(fEventsTime) delete fEventsTime;
+}
+
+void AliGenPythiaPlus::SetInteractionRate(Float_t rate,Float_t timewindow)
+{
+// Generate pileup using user specified rate
+ fInteractionRate = rate;
+ fTimeWindow = timewindow;
+ GeneratePileup();
+}
+
+void AliGenPythiaPlus::GeneratePileup()
+{
+// Generate sub events time for pileup
+ fEventsTime = 0;
+ if(fInteractionRate == 0.) {
+ Warning("GeneratePileup","Zero interaction specified. Skipping pileup generation.\n");
+ return;
+ }
+
+ Int_t npart = NumberParticles();
+ if(npart < 0) {
+ Warning("GeneratePileup","Negative number of particles. Skipping pileup generation.\n");
+ return;
+ }
+
+ if(fEventsTime) delete fEventsTime;
+ fEventsTime = new TArrayF(npart);
+ TArrayF &array = *fEventsTime;
+ for(Int_t ipart = 0; ipart < npart; ipart++)
+ array[ipart] = 0.;
+
+ Float_t eventtime = 0.;
+ while(1)
+ {
+ eventtime += (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
+ if(eventtime > fTimeWindow) break;
+ array.Set(array.GetSize()+1);
+ array[array.GetSize()-1] = eventtime;
+ }
+
+ eventtime = 0.;
+ while(1)
+ {
+ eventtime -= (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
+ if(TMath::Abs(eventtime) > fTimeWindow) break;
+ array.Set(array.GetSize()+1);
+ array[array.GetSize()-1] = eventtime;
+ }
+
+ SetNumberParticles(fEventsTime->GetSize());
+}
+
+void AliGenPythiaPlus::SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
+ Float_t thresh, Float_t etseed, Float_t minet, Float_t r)
+{
+// Set pycell parameters
+ fPycellEtaMax = etamax;
+ fPycellNEta = neta;
+ fPycellNPhi = nphi;
+ fPycellThreshold = thresh;
+ fPycellEtSeed = etseed;
+ fPycellMinEtJet = minet;
+ fPycellMaxRadius = r;
+}
+
+
+
+void AliGenPythiaPlus::SetEventListRange(Int_t eventFirst, Int_t eventLast)
+{
+ // Set a range of event numbers, for which a table
+ // of generated particle will be printed
+ fDebugEventFirst = eventFirst;
+ fDebugEventLast = eventLast;
+ if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst;
+}
+
+void AliGenPythiaPlus::Init()
+{
+// Initialisation
+
+// SetMC(AliPythia::Instance());
+// fPythia=(AliPythia*) fMCEvGen;
+
+//
+ fParentWeight=1./Float_t(fNpart);
+//
+
+
+ fPythia->SetPtHardRange(fPtHardMin, fPtHardMax);
+ fPythia->SetYHardRange(fYHardMin, fYHardMax);
+
+ if (fAProjectile > 0 && fATarget > 0) fPythia->SetNuclei(fAProjectile, fATarget);
+ // Fragmentation?
+ if (fFragmentation) {
+ fPythia->SetFragmentation(1);
+ } else {
+ fPythia->SetFragmentation(0);
+ }
+
+
+// initial state radiation
+ fPythia->SetInitialAndFinalStateRadiation(fGinit, fGfinal);
+
+// pt - kick
+ fPythia->SetIntrinsicKt(fPtKick);
+
+ if (fReadFromFile) {
+ fRL = AliRunLoader::Open(fFileName, "Partons");
+ fRL->LoadKinematics();
+ fRL->LoadHeader();
+ } else {
+ fRL = 0x0;
+ }
+ //
+ fPythia->ProcInit(fProcess, fEnergyCMS, fStrucFunc);
+ // Forward Paramters to the AliPythia object
+ fDecayer->SetForceDecay(fForceDecay);
+// Switch off Heavy Flavors on request
+ if (fHFoff) {
+ fPythia->SwitchHFOff();
+ // Switch off g->QQbar splitting in decay table
+ ((AliDecayerPythia*) fDecayer)->HeavyFlavourOff();
+ }
+
+ fDecayer->Init();
+
+
+// Parent and Children Selection
+ switch (fProcess)
+ {
+ case kPyOldUEQ2ordered:
+ case kPyOldUEQ2ordered2:
+ case kPyOldPopcorn:
+ break;
+ case kPyCharm:
+ case kPyCharmUnforced:
+ case kPyCharmPbPbMNR:
+ case kPyCharmpPbMNR:
+ case kPyCharmppMNR:
+ case kPyCharmppMNRwmi:
+ fParentSelect[0] = 411;
+ fParentSelect[1] = 421;
+ fParentSelect[2] = 431;
+ fParentSelect[3] = 4122;
+ fFlavorSelect = 4;
+ break;
+ case kPyD0PbPbMNR:
+ case kPyD0pPbMNR:
+ case kPyD0ppMNR:
+ fParentSelect[0] = 421;
+ fFlavorSelect = 4;
+ break;
+ case kPyDPlusPbPbMNR:
+ case kPyDPluspPbMNR:
+ case kPyDPlusppMNR:
+ fParentSelect[0] = 411;
+ fFlavorSelect = 4;
+ break;
+ case kPyDPlusStrangePbPbMNR:
+ case kPyDPlusStrangepPbMNR:
+ case kPyDPlusStrangeppMNR:
+ fParentSelect[0] = 431;
+ fFlavorSelect = 4;
+ break;
+ case kPyBeauty:
+ case kPyBeautyPbPbMNR:
+ case kPyBeautypPbMNR:
+ case kPyBeautyppMNR:
+ case kPyBeautyppMNRwmi:
+ fParentSelect[0]= 511;
+ fParentSelect[1]= 521;
+ fParentSelect[2]= 531;
+ fParentSelect[3]= 5122;
+ fParentSelect[4]= 5132;
+ fParentSelect[5]= 5232;
+ fParentSelect[6]= 5332;
+ fFlavorSelect = 5;
+ break;
+ case kPyBeautyUnforced:
+ fParentSelect[0] = 511;
+ fParentSelect[1] = 521;
+ fParentSelect[2] = 531;
+ fParentSelect[3] = 5122;
+ fParentSelect[4] = 5132;
+ fParentSelect[5] = 5232;
+ fParentSelect[6] = 5332;
+ fFlavorSelect = 5;
+ break;
+ case kPyJpsiChi:
+ case kPyJpsi:
+ fParentSelect[0] = 443;
+ break;
+ case kPyMbDefault:
+ case kPyMb:
+ case kPyMbNonDiffr:
+ case kPyMbMSEL1:
+ case kPyJets:
+ case kPyDirectGamma:
+ case kPyLhwgMb:
+ break;
+ case kPyW:
+ case kPyZ:
+ break;
+ }
+//
+//
+// JetFinder for Trigger
+//
+// Configure detector (EMCAL like)
+//
+ fPythia->SetPycellParameters(fPycellEtaMax,fPycellNEta, fPycellNPhi,
+ fPycellThreshold, fPycellEtSeed,
+ fPycellMinEtJet, fPycellMaxRadius);
+//
+// This counts the total number of calls to Pyevnt() per run.
+ fTrialsRun = 0;
+ fQ = 0.;
+ fX1 = 0.;
+ fX2 = 0.;
+ fNev = 0 ;
+//
+//
+//
+ AliGenMC::Init();
+//
+//
+//
+ if (fSetNuclei) {
+ fDyBoost = 0;
+ Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n");
+ }
+
+ if (fQuench) {
+ fPythia->InitQuenching(0., 0.1, 0.6e6, 0, 0.97, 30);
+ }
+
+// fPythia->SetPARJ(200, 0.0);
+
+// if (fQuench == 3) {
+// // Nestor's change of the splittings
+// fPythia->SetPARJ(200, 0.8);
+// fPythia->SetMSTJ(41, 1); // QCD radiation only
+// fPythia->SetMSTJ(42, 2); // angular ordering
+// fPythia->SetMSTJ(44, 2); // option to run alpha_s
+// fPythia->SetMSTJ(47, 0); // No correction back to hard scattering element
+// fPythia->SetMSTJ(50, 0); // No coherence in first branching
+// fPythia->SetPARJ(82, 1.); // Cut off for parton showers
+// }
+}
+
+void AliGenPythiaPlus::Generate()
+{
+// Generate one event
+
+ fDecayer->ForceDecay();
+
+ Float_t polar[3] = {0,0,0};
+ Float_t origin[3] = {0,0,0};
+ Float_t p[4];
+// converts from mm/c to s
+ const Float_t kconv=0.001/2.999792458e8;
+//
+ Int_t nt=0;
+ Int_t jev=0;
+ Int_t j, kf;
+ fTrials=0;
+ fEventTime = 0.;
+
+
+
+ // Set collision vertex position
+ if (fVertexSmear == kPerEvent) Vertex();
+
+// event loop
+ while(1)
+ {
+//
+// Produce event
+//
+//
+// Switch hadronisation off
+//
+// fPythia->SwitchHadronisationOff();
+//
+// Either produce new event or read partons from file
+//
+ if (!fReadFromFile) {
+ if (!fNewMIS) {
+ fPythia->GenerateEvent();
+ } else {
+ fPythia->GenerateMIEvent();
+ }
+ fNpartons = fPythia->GetNumberOfParticles();
+ } else {
+ printf("Loading Event %d\n",AliRunLoader::GetRunLoader()->GetEventNumber());
+ fRL->GetEvent(AliRunLoader::GetRunLoader()->GetEventNumber());
+ fPythia->SetNumberOfParticles(0);
+ fPythia->LoadEvent(fRL->Stack(), 0 , 1);
+ fPythia->EditEventList(21);
+ }
+
+//
+// Run quenching routine
+//
+ if (fQuench == 1) {
+ fPythia->Quench();
+ } else if (fQuench == 2){
+ fPythia->Pyquen(208., 0, 0.);
+ } else if (fQuench == 3) {
+ // Quenching is via multiplicative correction of the splittings
+ }
+
+//
+// Switch hadronisation on
+//
+// fPythia->SwitchHadronisationOn();
+//
+// .. and perform hadronisation
+// printf("Calling hadronisation %d\n", fPythia->GetN());
+// fPythia->HadronizeEvent();
+ fTrials++;
+ fPythia->GetParticles(fParticles);
+ Boost();
+//
+//
+//
+ Int_t i;
+
+ fNprimaries = 0;
+ Int_t np = fParticles->GetEntriesFast();
+
+ if (np == 0) continue;
+//
+
+//
+ Int_t* pParent = new Int_t[np];
+ Int_t* pSelected = new Int_t[np];
+ Int_t* trackIt = new Int_t[np];
+ for (i = 0; i < np; i++) {
+ pParent[i] = -1;
+ pSelected[i] = 0;
+ trackIt[i] = 0;
+ }
+
+ Int_t nc = 0; // Total n. of selected particles
+ Int_t nParents = 0; // Selected parents
+ Int_t nTkbles = 0; // Trackable particles
+ if (fProcess != kPyMbDefault &&
+ fProcess != kPyMb &&
+ fProcess != kPyJets &&
+ fProcess != kPyDirectGamma &&
+ fProcess != kPyMbNonDiffr &&
+ fProcess != kPyMbMSEL1 &&
+ fProcess != kPyW &&
+ fProcess != kPyZ &&
+ fProcess != kPyCharmppMNRwmi &&
+ fProcess != kPyBeautyppMNRwmi) {
+
+ for (i = 0; i < np; i++) {
+ TParticle* iparticle = (TParticle *) fParticles->At(i);
+ Int_t ks = iparticle->GetStatusCode();
+ kf = CheckPDGCode(iparticle->GetPdgCode());
+// No initial state partons
+ if (ks==21) continue;
+//
+// Heavy Flavor Selection
+//
+ // quark ?
+ kf = TMath::Abs(kf);
+ Int_t kfl = kf;
+ // Resonance
+
+ if (kfl > 100000) kfl %= 100000;
+ if (kfl > 10000) kfl %= 10000;
+ // meson ?
+ if (kfl > 10) kfl/=100;
+ // baryon
+ if (kfl > 10) kfl/=10;
+ Int_t ipa = (fPythia->Version() == 6) ? (iparticle->GetFirstMother() - 1) :(iparticle->GetFirstMother()) ;
+ Int_t kfMo = 0;
+//
+// Establish mother daughter relation between heavy quarks and mesons
+//
+ if (kf >= fFlavorSelect && kf <= 6) {
+ Int_t idau = (fPythia->Version() == 6) ? (iparticle->GetFirstDaughter() - 1) :(iparticle->GetFirstDaughter());
+ if (idau > -1) {
+ TParticle* daughter = (TParticle *) fParticles->At(idau);
+ Int_t pdgD = daughter->GetPdgCode();
+ if (pdgD == 91 || pdgD == 92) {
+ Int_t jmin = (fPythia->Version() == 6) ? (daughter->GetFirstDaughter() - 1) : (daughter->GetFirstDaughter());
+ Int_t jmax = (fPythia->Version() == 6) ? (daughter->GetLastDaughter() - 1) : (daughter->GetLastDaughter());
+
+ for (Int_t j = jmin; j <= jmax; j++)
+ ((TParticle *) fParticles->At(j))->SetFirstMother(i+1);
+ } // is string or cluster
+ } // has daughter
+ } // heavy quark
+
+
+ if (ipa > -1) {
+ TParticle * mother = (TParticle *) fParticles->At(ipa);
+ kfMo = TMath::Abs(mother->GetPdgCode());
+ }
+
+ // What to keep in Stack?
+ Bool_t flavorOK = kFALSE;
+ Bool_t selectOK = kFALSE;
+ if (fFeedDownOpt) {
+ if (kfl >= fFlavorSelect) flavorOK = kTRUE;
+ } else {
+ if (kfl > fFlavorSelect) {
+ nc = -1;
+ break;
+ }
+ if (kfl == fFlavorSelect) flavorOK = kTRUE;
+ }
+ switch (fStackFillOpt) {
+ case kFlavorSelection:
+ selectOK = kTRUE;
+ break;
+ case kParentSelection:
+ if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
+ break;
+ }
+ if (flavorOK && selectOK) {
+//
+// Heavy flavor hadron or quark
+//
+// Kinematic seletion on final state heavy flavor mesons
+ if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
+ {
+ continue;
+ }
+ pSelected[i] = 1;
+ if (ParentSelected(kf)) ++nParents; // Update parent count
+// printf("\n particle (HF) %d %d %d", i, pSelected[i], kf);
+ } else {
+// Kinematic seletion on decay products
+ if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf)
+ && !KinematicSelection(iparticle, 1))
+ {
+ continue;
+ }
+//
+// Decay products
+// Select if mother was selected and is not tracked
+
+ if (pSelected[ipa] &&
+ !trackIt[ipa] && // mother will be tracked ?
+ kfMo != 5 && // mother is b-quark, don't store fragments
+ kfMo != 4 && // mother is c-quark, don't store fragments
+ kf != 92) // don't store string
+ {
+//
+// Semi-stable or de-selected: diselect decay products:
+//
+//
+ if (pSelected[i] == -1 || fDecayer->GetLifetime(kf) > fMaxLifeTime)
+ {
+ Int_t ipF = iparticle->GetFirstDaughter();
+ Int_t ipL = iparticle->GetLastDaughter();
+ if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
+ }
+// printf("\n particle (decay) %d %d %d", i, pSelected[i], kf);
+ pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
+ }
+ }
+ if (pSelected[i] == -1) pSelected[i] = 0;
+ if (!pSelected[i]) continue;
+ // Count quarks only if you did not include fragmentation
+ if (fFragmentation && kf <= 10) continue;
+
+ nc++;
+// Decision on tracking
+ trackIt[i] = 0;
+//
+// Track final state particle
+ if (ks == 1) trackIt[i] = 1;
+// Track semi-stable particles
+ if ((ks == 1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime)) trackIt[i] = 1;
+// Track particles selected by process if undecayed.
+ if (fForceDecay == kNoDecay) {
+ if (ParentSelected(kf)) trackIt[i] = 1;
+ } else {
+ if (ParentSelected(kf)) trackIt[i] = 0;
+ }
+ if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
+//
+//
+
+ } // particle selection loop
+ if (nc > 0) {
+ for (i = 0; i < np; i++) {
+ if (!pSelected[i]) continue;
+ TParticle * iparticle = (TParticle *) fParticles->At(i);
+ kf = CheckPDGCode(iparticle->GetPdgCode());
+ Int_t ks = iparticle->GetStatusCode();
+ p[0] = iparticle->Px();
+ p[1] = iparticle->Py();
+ p[2] = iparticle->Pz();
+ p[3] = iparticle->Energy();
+
+ origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
+ origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
+ origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
+
+ Float_t tof = kconv*iparticle->T();
+ Int_t ipa = (fPythia->Version() == 6) ? (iparticle->GetFirstMother() - 1) :(iparticle->GetFirstMother()) ;
+ Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
+
+ PushTrack(fTrackIt*trackIt[i], iparent, kf,
+ p[0], p[1], p[2], p[3],
+ origin[0], origin[1], origin[2], tof,
+ polar[0], polar[1], polar[2],
+ kPPrimary, nt, 1., ks);
+ pParent[i] = nt;
+ KeepTrack(nt);
+ fNprimaries++;
+ } // PushTrack loop
+ }
+ } else {
+ nc = GenerateMB();
+ } // mb ?
+
+ GetSubEventTime();
+
+ delete[] pParent;
+ delete[] pSelected;
+ delete[] trackIt;
+
+ if (nc > 0) {
+ switch (fCountMode) {
+ case kCountAll:
+ // printf(" Count all \n");
+ jev += nc;
+ break;
+ case kCountParents:
+ // printf(" Count parents \n");
+ jev += nParents;
+ break;
+ case kCountTrackables:
+ // printf(" Count trackable \n");
+ jev += nTkbles;
+ break;
+ }
+ if (jev >= fNpart || fNpart == -1) {
+ fKineBias=Float_t(fNpart)/Float_t(fTrials);
+ fPythia->GetXandQ(fQ, fX1, fX2);
+ fTrialsRun += fTrials;
+ fNev++;
+ MakeHeader();
+ break;
+ }
+ }
+ } // event loop
+ SetHighWaterMark(nt);
+// Adjust weight due to kinematic selection
+// AdjustWeights();
+// get cross-section
+ fXsection = fPythia->GetXSection();
+}
+
+Int_t AliGenPythiaPlus::GenerateMB()
+{
+//
+// Min Bias selection and other global selections
+//
+ Int_t i, kf, nt, iparent;
+ Int_t nc = 0;
+ Float_t p[4];
+ Float_t polar[3] = {0,0,0};
+ Float_t origin[3] = {0,0,0};
+// converts from mm/c to s
+ const Float_t kconv = 0.001 / 2.999792458e8;
+
+ Int_t np = (fHadronisation) ? fParticles->GetEntriesFast() : fNpartons;
+
+ Int_t* pParent = new Int_t[np];
+ for (i=0; i< np; i++) pParent[i] = -1;
+ if (fProcess == kPyJets || fProcess == kPyDirectGamma) {
+ TParticle* jet1 = (TParticle *) fParticles->At(6);
+ TParticle* jet2 = (TParticle *) fParticles->At(7);
+ if (!CheckTrigger(jet1, jet2)) {
+ delete [] pParent;
+ return 0;
+ }
+ }
+
+ // Select jets with fragmentation photon or pi0 going to PHOS or EMCAL
+ if (fProcess == kPyJets && (fFragPhotonInCalo || fPi0InCalo) ) {
+
+ Bool_t ok = kFALSE;
+
+ Int_t pdg = 0;
+ if (fFragPhotonInCalo) pdg = 22 ; // Photon
+ else if (fPi0InCalo) pdg = 111 ; // Pi0
+
+ for (i=0; i< np; i++) {
+ TParticle* iparticle = (TParticle *) fParticles->At(i);
+ if(iparticle->GetStatusCode()==1 && iparticle->GetPdgCode()==pdg &&
+ iparticle->Pt() > fFragPhotonOrPi0MinPt){
+ Int_t imother = (fPythia->Version() == 6) ? (iparticle->GetFirstMother() - 1) :(iparticle->GetFirstMother()) ;
+ TParticle* pmother = (TParticle *) fParticles->At(imother);
+ if(pdg == 111 ||
+ (pdg == 22 && pmother->GetStatusCode() != 11))//No photon from hadron decay
+ {
+ Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
+ Float_t eta =TMath::Abs(iparticle->Eta());//in calos etamin=-etamax
+ if((fCheckEMCAL && IsInEMCAL(phi,eta)) ||
+ (fCheckPHOS && IsInPHOS(phi,eta)) )
+ ok =kTRUE;
+ }
+ }
+ }
+ if(!ok)
+ return 0;
+ }
+
+
+ // Select events with a photon pt > min pt going to PHOS eta acceptance or exactly PHOS eta phi
+ if ((fProcess == kPyJets || fProcess == kPyDirectGamma) && fPhotonInCalo && (fCheckPHOSeta || fCheckPHOS)){
+
+ Bool_t okd = kFALSE;
+
+ Int_t pdg = 22;
+ Int_t iphcand = -1;
+ for (i=0; i< np; i++) {
+ TParticle* iparticle = (TParticle *) fParticles->At(i);
+ Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
+ Float_t eta =TMath::Abs(iparticle->Eta());//in calos etamin=-etamax
+
+ if(iparticle->GetStatusCode() == 1
+ && iparticle->GetPdgCode() == pdg
+ && iparticle->Pt() > fPhotonMinPt
+ && eta < fPHOSEta){
+
+ // first check if the photon is in PHOS phi
+ if(IsInPHOS(phi,eta)){
+ okd = kTRUE;
+ break;
+ }
+ if(fCheckPHOSeta) iphcand = i; // candiate photon to rotate in phi
+
+ }
+ }
+
+ if(!okd && iphcand != -1) // execute rotation in phi
+ RotatePhi(iphcand,okd);
+
+ if(!okd)
+ return 0;
+ }
+
+ if (fTriggerParticle) {
+ Bool_t triggered = kFALSE;
+ for (i = 0; i < np; i++) {
+ TParticle * iparticle = (TParticle *) fParticles->At(i);
+ kf = CheckPDGCode(iparticle->GetPdgCode());
+ if (kf != fTriggerParticle) continue;
+ if (iparticle->Pt() == 0.) continue;
+ if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
+ triggered = kTRUE;
+ break;
+ }
+ if (!triggered) {
+ delete [] pParent;
+ return 0;
+ }
+ }
+
+
+ // Check if there is a ccbar or bbbar pair with at least one of the two
+ // in fYMin < y < fYMax
+ if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi) {
+ TParticle *hvq;
+ Bool_t theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
+ Float_t yQ;
+ Int_t pdgQ;
+ for(i=0; i<np; i++) {
+ hvq = (TParticle*)fParticles->At(i);
+ pdgQ = hvq->GetPdgCode();
+ if(TMath::Abs(pdgQ) != fFlavorSelect) continue;
+ if(pdgQ>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
+ yQ = 0.5*TMath::Log((hvq->Energy()+hvq->Pz()+1.e-13)/
+ (hvq->Energy()-hvq->Pz()+1.e-13));
+ if(yQ>fYMin && yQ<fYMax) inYcut=kTRUE;
+ }
+ if (!theQ || !theQbar || !inYcut) {
+ delete[] pParent;
+ return 0;
+ }
+ }
+
+ //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
+ if ( (fProcess == kPyW ||
+ fProcess == kPyZ ||
+ fProcess == kPyMbDefault ||
+ fProcess == kPyMb ||
+ fProcess == kPyMbNonDiffr)
+ && (fCutOnChild == 1) ) {
+ if ( !CheckKinematicsOnChild() ) {
+ delete[] pParent;
+ return 0;
+ }
+ }
+
+
+ for (i = 0; i < np; i++) {
+ Int_t trackIt = 0;
+ TParticle * iparticle = (TParticle *) fParticles->At(i);
+ kf = CheckPDGCode(iparticle->GetPdgCode());
+ Int_t ks = iparticle->GetStatusCode();
+ Int_t km = iparticle->GetFirstMother();
+ if ((ks == 1 && kf!=0 && KinematicSelection(iparticle, 0)) ||
+ (ks != 1) ||
+ (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
+ nc++;
+ if (ks == 1) trackIt = 1;
+
+ Int_t ipa = (fPythia->Version() == 6) ? (iparticle->GetFirstMother() - 1) :(iparticle->GetFirstMother()) ;
+ iparent = (ipa > -1) ? pParent[ipa] : -1;
+ if (ipa >= np) fPythia->EventListing();
+
+//
+// store track information
+ p[0] = iparticle->Px();
+ p[1] = iparticle->Py();
+ p[2] = iparticle->Pz();
+ p[3] = iparticle->Energy();
+
+
+ origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
+ origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
+ origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
+
+ Float_t tof = fEventTime + kconv * iparticle->T();
+
+ PushTrack(fTrackIt*trackIt, iparent, kf,
+ p[0], p[1], p[2], p[3],
+ origin[0], origin[1], origin[2], tof,
+ polar[0], polar[1], polar[2],
+ kPPrimary, nt, 1., ks);
+ fNprimaries++;
+
+
+ //
+ // Special Treatment to store color-flow
+ //
+ if (ks == 3 || ks == 13 || ks == 14) {
+ TParticle* particle = 0;
+ if (fStack) {
+ particle = fStack->Particle(nt);
+ } else {
+ particle = gAlice->Stack()->Particle(nt);
+ }
+// particle->SetFirstDaughter(fPythia->GetK(2, i));
+// particle->SetLastDaughter(fPythia->GetK(3, i));
+ }
+
+ KeepTrack(nt);
+ pParent[i] = nt;
+ SetHighWaterMark(nt);
+
+ } // select particle
+ } // particle loop
+
+ delete[] pParent;
+
+ return 1;
+}
+
+
+void AliGenPythiaPlus::FinishRun()
+{
+// Print x-section summary
+ fPythia->PrintStatistics();
+
+ if (fNev > 0.) {
+ fQ /= fNev;
+ fX1 /= fNev;
+ fX2 /= fNev;
+ }
+
+ printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
+ printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
+}
+
+void AliGenPythiaPlus::AdjustWeights() const
+{
+// Adjust the weights after generation of all events
+//
+ if (gAlice) {
+ TParticle *part;
+ Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
+ for (Int_t i=0; i<ntrack; i++) {
+ part= gAlice->GetMCApp()->Particle(i);
+ part->SetWeight(part->GetWeight()*fKineBias);
+ }
+ }
+}
+
+void AliGenPythiaPlus::SetNuclei(Int_t a1, Int_t a2)
+{
+// Treat protons as inside nuclei with mass numbers a1 and a2
+
+ fAProjectile = a1;
+ fATarget = a2;
+ fSetNuclei = kTRUE;
+}
+
+
+void AliGenPythiaPlus::MakeHeader()
+{
+//
+// Make header for the simulated event
+//
+ if (gAlice) {
+ if (gAlice->GetEvNumber()>=fDebugEventFirst &&
+ gAlice->GetEvNumber()<=fDebugEventLast) fPythia->EventListing();
+ }
+
+// Builds the event header, to be called after each event
+ if (fHeader) delete fHeader;
+ fHeader = new AliGenPythiaEventHeader("Pythia");
+//
+// Event type
+ ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->ProcessCode());
+//
+// Number of trials
+ ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
+//
+// Event Vertex
+ fHeader->SetPrimaryVertex(fVertex);
+
+//
+// Number of primaries
+ fHeader->SetNProduced(fNprimaries);
+//
+// Jets that have triggered
+
+ if (fProcess == kPyJets)
+ {
+ Int_t ntrig, njet;
+ Float_t jets[4][10];
+ GetJets(njet, ntrig, jets);
+
+
+ for (Int_t i = 0; i < ntrig; i++) {
+ ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
+ jets[3][i]);
+ }
+ }
+//
+// Copy relevant information from external header, if present.
+//
+ Float_t uqJet[4];
+
+ if (fRL) {
+ AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
+ for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
+ {
+ printf("Adding Jet %d %d \n", i, exHeader->NTriggerJets());
+
+
+ exHeader->TriggerJet(i, uqJet);
+ ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
+ }
+ }
+//
+// Store quenching parameters
+//
+ if (fQuench){
+ Double_t z[4];
+ Double_t xp, yp;
+ if (fQuench == 1) {
+ // Pythia::Quench()
+ fPythia->GetQuenchingParameters(xp, yp, z);
+ } else {
+ // Pyquen
+ Double_t r1 = PARIMP.rb1;
+ Double_t r2 = PARIMP.rb2;
+ Double_t b = PARIMP.b1;
+ Double_t r = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
+ Double_t phi = PARIMP.psib1;
+ xp = r * TMath::Cos(phi);
+ yp = r * TMath::Sin(phi);
+
+ }
+ ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
+ ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
+ }
+//
+// Store pt^hard
+ ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetPtHard());
+//
+// Pass header
+//
+ AddHeader(fHeader);
+ fHeader = 0x0;
+}
+
+Bool_t AliGenPythiaPlus::CheckTrigger(TParticle* jet1, TParticle* jet2)
+{
+// Check the kinematic trigger condition
+//
+ Double_t eta[2];
+ eta[0] = jet1->Eta();
+ eta[1] = jet2->Eta();
+ Double_t phi[2];
+ phi[0] = jet1->Phi();
+ phi[1] = jet2->Phi();
+ Int_t pdg[2];
+ pdg[0] = jet1->GetPdgCode();
+ pdg[1] = jet2->GetPdgCode();
+ Bool_t triggered = kFALSE;
+
+ if (fProcess == kPyJets) {
+ Int_t njets = 0;
+ Int_t ntrig = 0;
+ Float_t jets[4][10];
+//
+// Use Pythia clustering on parton level to determine jet axis
+//
+ GetJets(njets, ntrig, jets);
+
+ if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
+//
+ } else {
+ Int_t ij = 0;
+ Int_t ig = 1;
+ if (pdg[0] == kGamma) {
+ ij = 1;
+ ig = 0;
+ }
+ //Check eta range first...
+ if ((eta[ij] < fEtaMaxJet && eta[ij] > fEtaMinJet) &&
+ (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
+ {
+ //Eta is okay, now check phi range
+ if ((phi[ij] < fPhiMaxJet && phi[ij] > fPhiMinJet) &&
+ (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
+ {
+ triggered = kTRUE;
+ }
+ }
+ }
+ return triggered;
+}
+
+
+
+Bool_t AliGenPythiaPlus::CheckKinematicsOnChild(){
+//
+//Checking Kinematics on Child (status code 1, particle code ?, kin cuts
+//
+ Bool_t checking = kFALSE;
+ Int_t j, kcode, ks, km;
+ Int_t nPartAcc = 0; //number of particles in the acceptance range
+ Int_t numberOfAcceptedParticles = 1;
+ if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
+ Int_t npart = fParticles->GetEntriesFast();
+
+ for (j = 0; j<npart; j++) {
+ TParticle * jparticle = (TParticle *) fParticles->At(j);
+ kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
+ ks = jparticle->GetStatusCode();
+ km = jparticle->GetFirstMother();
+
+ if( (ks == 1) && (kcode == fPdgCodeParticleforAcceptanceCut) && (KinematicSelection(jparticle,1)) ){
+ nPartAcc++;
+ }
+ if( numberOfAcceptedParticles <= nPartAcc){
+ checking = kTRUE;
+ break;
+ }
+ }
+
+ return checking;
+}
+
+void AliGenPythiaPlus::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
+{
+//
+// Calls the Pythia jet finding algorithm to find jets in the current event
+//
+//
+//
+// Save jets
+//
+// Run Jet Finder
+ fPythia->Pycell(njets);
+ Int_t i;
+ for (i = 0; i < njets; i++) {
+ Float_t px, py, pz, e;
+ fPythia->GetJet(i, px, py, pz, e);
+ jets[0][i] = px;
+ jets[1][i] = py;
+ jets[2][i] = pz;
+ jets[3][i] = e;
+ }
+}
+
+
+void AliGenPythiaPlus::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
+{
+//
+// Calls the Pythia clustering algorithm to find jets in the current event
+//
+ nJets = 0;
+ nJetsTrig = 0;
+ if (fJetReconstruction == kCluster) {
+//
+// Configure cluster algorithm
+//
+// fPythia->SetPARU(43, 2.);
+// fPythia->SetMSTU(41, 1);
+//
+// Call cluster algorithm
+//
+ fPythia->Pyclus(nJets);
+//
+// Loading jets from common block
+//
+ } else {
+
+//
+// Run Jet Finder
+ fPythia->Pycell(nJets);
+ }
+
+ Int_t i;
+ for (i = 0; i < nJets; i++) {
+ Float_t px, py, pz, e;
+ fPythia->GetJet(i, px, py, pz, e);
+ Float_t pt = TMath::Sqrt(px * px + py * py);
+ Float_t phi = TMath::Pi() + TMath::ATan2(-py, -px);
+ Float_t theta = TMath::ATan2(pt,pz);
+ Float_t et = e * TMath::Sin(theta);
+ Float_t eta = -TMath::Log(TMath::Tan(theta / 2.));
+ if (
+ eta > fEtaMinJet && eta < fEtaMaxJet &&
+ phi > fPhiMinJet && phi < fPhiMaxJet &&
+ et > fEtMinJet && et < fEtMaxJet
+ )
+ {
+ jets[0][nJetsTrig] = px;
+ jets[1][nJetsTrig] = py;
+ jets[2][nJetsTrig] = pz;
+ jets[3][nJetsTrig] = e;
+ nJetsTrig++;
+ } else {
+ }
+ }
+}
+
+void AliGenPythiaPlus::GetSubEventTime()
+{
+ // Calculates time of the next subevent
+ fEventTime = 0.;
+ if (fEventsTime) {
+ TArrayF &array = *fEventsTime;
+ fEventTime = array[fCurSubEvent++];
+ }
+ // printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
+ return;
+}
+
+
+
+
+Bool_t AliGenPythiaPlus::IsInEMCAL(Float_t phi, Float_t eta)
+{
+ // Is particle in EMCAL acceptance?
+ // phi in degrees, etamin=-etamax
+ if(phi > fEMCALMinPhi && phi < fEMCALMaxPhi &&
+ eta < fEMCALEta )
+ return kTRUE;
+ else
+ return kFALSE;
+}
+
+Bool_t AliGenPythiaPlus::IsInPHOS(Float_t phi, Float_t eta)
+{
+ // Is particle in PHOS acceptance?
+ // Acceptance slightly larger considered.
+ // phi in degrees, etamin=-etamax
+ if(phi > fPHOSMinPhi && phi < fPHOSMaxPhi &&
+ eta < fPHOSEta )
+ return kTRUE;
+ else
+ return kFALSE;
+}
+
+void AliGenPythiaPlus::RotatePhi(Int_t iphcand, Bool_t& okdd)
+{
+ //calculate the new position random between fPHOSMinPhi and fPHOSMaxPhi
+ Double_t phiPHOSmin = TMath::Pi()*fPHOSMinPhi/180;
+ Double_t phiPHOSmax = TMath::Pi()*fPHOSMaxPhi/180;
+ Double_t phiPHOS = gRandom->Uniform(phiPHOSmin,phiPHOSmax);
+
+ //calculate deltaphi
+ TParticle* ph = (TParticle *) fParticles->At(iphcand);
+ Double_t phphi = ph->Phi();
+ Double_t deltaphi = phiPHOS - phphi;
+
+
+
+ //loop for all particles and produce the phi rotation
+ Int_t np = (fHadronisation) ? fParticles->GetEntriesFast() : fNpartons;
+ Double_t oldphi, newphi;
+ Double_t newVx, newVy, R, Vz, time;
+ Double_t newPx, newPy, pt, Pz, e;
+ for(Int_t i=0; i< np; i++) {
+ TParticle* iparticle = (TParticle *) fParticles->At(i);
+ oldphi = iparticle->Phi();
+ newphi = oldphi + deltaphi;
+ if(newphi < 0) newphi = 2*TMath::Pi() + newphi; // correct angle
+ if(newphi > 2*TMath::Pi()) newphi = newphi - 2*TMath::Pi(); // correct angle
+
+ R = iparticle->R();
+ newVx = R*TMath::Cos(newphi);
+ newVy = R*TMath::Sin(newphi);
+ Vz = iparticle->Vz(); // don't transform
+ time = iparticle->T(); // don't transform
+
+ pt = iparticle->Pt();
+ newPx = pt*TMath::Cos(newphi);
+ newPy = pt*TMath::Sin(newphi);
+ Pz = iparticle->Pz(); // don't transform
+ e = iparticle->Energy(); // don't transform
+
+ // apply rotation
+ iparticle->SetProductionVertex(newVx, newVy, Vz, time);
+ iparticle->SetMomentum(newPx, newPy, Pz, e);
+
+ } //end particle loop
+
+ // now let's check that we put correctly the candidate photon in PHOS
+ Float_t phi = ph->Phi()*180./TMath::Pi(); //Convert to degrees
+ Float_t eta =TMath::Abs(ph->Eta());//in calos etamin=-etamax
+ if(IsInPHOS(phi,eta))
+ okdd = kTRUE;
+}
+
+
+#ifdef never
+void AliGenPythiaPlus::Streamer(TBuffer &R__b)
+{
+ // Stream an object of class AliGenPythia.
+
+ if (R__b.IsReading()) {
+ Version_t R__v = R__b.ReadVersion(); if (R__v) { }
+ AliGenerator::Streamer(R__b);
+ R__b >> (Int_t&)fProcess;
+ R__b >> (Int_t&)fStrucFunc;
+ R__b >> (Int_t&)fForceDecay;
+ R__b >> fEnergyCMS;
+ R__b >> fKineBias;
+ R__b >> fTrials;
+ fParentSelect.Streamer(R__b);
+ fChildSelect.Streamer(R__b);
+ R__b >> fXsection;
+// (AliPythia::Instance())->Streamer(R__b);
+ R__b >> fPtHardMin;
+ R__b >> fPtHardMax;
+// if (fDecayer) fDecayer->Streamer(R__b);
+ } else {
+ R__b.WriteVersion(AliGenPythiaPlus::IsA());
+ AliGenerator::Streamer(R__b);
+ R__b << (Int_t)fProcess;
+ R__b << (Int_t)fStrucFunc;
+ R__b << (Int_t)fForceDecay;
+ R__b << fEnergyCMS;
+ R__b << fKineBias;
+ R__b << fTrials;
+ fParentSelect.Streamer(R__b);
+ fChildSelect.Streamer(R__b);
+ R__b << fXsection;
+// R__b << fPythia;
+ R__b << fPtHardMin;
+ R__b << fPtHardMax;
+ // fDecayer->Streamer(R__b);
+ }
+}
+#endif
+
+
+
--- /dev/null
+#ifndef ALIGENPYTHIAPLUS_H
+#define ALIGENPYTHIAPLUS_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+
+/* $Id$ */
+
+//
+// Generator using the TPythia interface (via AliPythia)
+// to generate pp collisions.
+// Using SetNuclei() also nuclear modifications to the structure functions
+// can be taken into account. This makes, of course, only sense for the
+// generation of the products of hard processes (heavy flavor, jets ...)
+//
+// andreas.morsch@cern.ch
+//
+
+#include "AliGenMC.h"
+#include "AliPythia.h"
+
+class AliPythiaBase;
+class TParticle;
+class AliGenPythiaEventHeader;
+class AliGenEventHeader;
+class AliStack;
+class AliRunLoader;
+
+class AliGenPythiaPlus : public AliGenMC
+{
+ public:
+
+ typedef enum {kFlavorSelection, kParentSelection} StackFillOpt_t;
+ typedef enum {kCountAll, kCountParents, kCountTrackables} CountMode_t;
+ typedef enum {kCluster, kCell} JetRecMode_t;
+
+ AliGenPythiaPlus();
+ AliGenPythiaPlus(AliPythiaBase* pythia);
+ virtual ~AliGenPythiaPlus();
+ virtual void Generate();
+ virtual void Init();
+ // Range of events to be printed
+ virtual void SetEventListRange(Int_t eventFirst=-1, Int_t eventLast=-1);
+ // Select process type
+ virtual void SetProcess(Process_t proc = kPyCharm) {fProcess = proc;}
+ // Select structure function
+ virtual void SetStrucFunc(StrucFunc_t func = kCTEQ5L) {fStrucFunc = func;}
+ // Select pt of hard scattering
+ virtual void SetPtHard(Float_t ptmin = 0, Float_t ptmax = 1.e10)
+ {fPtHardMin = ptmin; fPtHardMax = ptmax; }
+ // y of hard scattering
+ virtual void SetYHard(Float_t ymin = -1.e10, Float_t ymax = 1.e10)
+ {fYHardMin = ymin; fYHardMax = ymax; }
+ // Set initial and final state gluon radiation
+ virtual void SetGluonRadiation(Int_t iIn, Int_t iFin)
+ {fGinit = iIn; fGfinal = iFin;}
+ // Intrinsic kT
+ virtual void SetPtKick(Float_t kt = 1.)
+ {fPtKick = kt;}
+ // Use the Pythia 6.3 new multiple interations scenario
+ virtual void UseNewMultipleInteractionsScenario() {fNewMIS = kTRUE;}
+ // Switch off heavy flavors
+ virtual void SwitchHFOff() {fHFoff = kTRUE;}
+ // Set centre of mass energy
+ virtual void SetEnergyCMS(Float_t energy = 5500) {fEnergyCMS = energy;}
+ // Treat protons as inside nuclei with mass numbers a1 and a2
+ virtual void SetNuclei(Int_t a1, Int_t a2);
+ //
+ // Trigger options
+ //
+ // Energy range for jet trigger
+ virtual void SetJetEtRange(Float_t etmin = 0., Float_t etmax = 1.e4)
+ {fEtMinJet = etmin; fEtMaxJet = etmax;}
+ // Eta range for jet trigger
+ virtual void SetJetEtaRange(Float_t etamin = -20., Float_t etamax = 20.)
+ {fEtaMinJet = etamin; fEtaMaxJet = etamax;}
+ // Phi range for jet trigger
+ virtual void SetJetPhiRange(Float_t phimin = 0., Float_t phimax = 360.)
+ {fPhiMinJet = TMath::Pi()*phimin/180.; fPhiMaxJet = TMath::Pi()*phimax/180.;}
+ // Jet reconstruction mode; default is cone algorithm
+ virtual void SetJetReconstructionMode(Int_t mode = kCell) {fJetReconstruction = mode;}
+ // Eta range for gamma trigger
+ virtual void SetGammaEtaRange(Float_t etamin = -20., Float_t etamax = 20.)
+ {fEtaMinGamma = etamin; fEtaMaxGamma = etamax;}
+ // Phi range for gamma trigger
+ virtual void SetGammaPhiRange(Float_t phimin = 0., Float_t phimax = 360.)
+ {fPhiMinGamma = TMath::Pi()*phimin/180.; fPhiMaxGamma = TMath::Pi()*phimax/180.;}
+ // Select jets with fragmentation photon or pi0 going to PHOS or EMCAL
+ virtual void SetFragPhotonInCalo(Bool_t b) {fFragPhotonInCalo = b;}
+ virtual void SetPi0InCalo (Bool_t b) {fPi0InCalo = b;}
+ virtual void SetPhotonInCalo(Bool_t b) {fPhotonInCalo = b;}
+ virtual void SetCheckPHOS (Bool_t b) {fCheckPHOS = b;}
+ virtual void SetCheckEMCAL(Bool_t b) {fCheckEMCAL = b;}
+ virtual void SetFragPhotonInEMCAL(Bool_t b) {fCheckEMCAL = b; fFragPhotonInCalo = b;}
+ virtual void SetFragPhotonInPHOS(Bool_t b) {fCheckPHOS = b; fFragPhotonInCalo = b;}
+ virtual void SetPi0InEMCAL(Bool_t b) {fCheckEMCAL = b; fPi0InCalo = b;}
+ virtual void SetPi0InPHOS(Bool_t b) {fCheckPHOS = b; fPi0InCalo = b;}
+ virtual void SetPhotonInEMCAL(Bool_t b) {fCheckEMCAL = b; fPhotonInCalo = b;}
+ virtual void SetPhotonInPHOS(Bool_t b) {fCheckPHOS = b; fPhotonInCalo = b;}
+ virtual void SetPhotonInPHOSeta(Bool_t b) {fCheckPHOSeta = b; fPhotonInCalo = b;}
+ virtual void SetFragPhotonOrPi0MinPt(Float_t pt) {fFragPhotonOrPi0MinPt = pt;}
+ virtual void SetPhotonMinPt(Float_t pt) {fPhotonMinPt = pt;}
+ // Trigger and rotate event
+ void RotatePhi(Int_t iphcand, Bool_t& okdd);
+ // Trigger on a single particle
+ virtual void SetTriggerParticle(Int_t particle = 0, Float_t etamax = 0.9)
+ {fTriggerParticle = particle; fTriggerEta = etamax;}
+ //
+ // Heavy flavor options
+ //
+ // Set option for feed down from higher family
+ virtual void SetFeedDownHigherFamily(Bool_t opt) {
+ fFeedDownOpt = opt;
+ }
+ // Set option for selecting particles kept in stack according to flavor
+ // or to parent selection
+ virtual void SetStackFillOpt(StackFillOpt_t opt) {
+ fStackFillOpt = opt;
+ }
+ // Set fragmentation option
+ virtual void SetFragmentation(Bool_t opt) {
+ fFragmentation = opt;
+ }
+ // Set counting mode
+ virtual void SetCountMode(CountMode_t mode) {
+ fCountMode = mode;
+ }
+ //
+ // Quenching
+ //
+ // Set quenching mode 0 = no, 1 = AM, 2 = IL
+ virtual void SetQuench(Int_t flag = 0) {fQuench = flag;}
+ virtual void SetHadronisation(Int_t flag = 1) {fHadronisation = flag;}
+ virtual void SetReadFromFile(const Text_t *filname) {fFileName = filname; fReadFromFile = 1;}
+
+ //
+ // Pile-up
+ //
+ // Get interaction rate for pileup studies
+ virtual void SetInteractionRate(Float_t rate,Float_t timewindow = 90.e-6);
+ virtual Float_t GetInteractionRate() const {return fInteractionRate;}
+ // Get cross section of process
+ virtual Float_t GetXsection() const {return fXsection;}
+ // Get triggered jets
+ void GetJets(Int_t& njets, Int_t& ntrig, Float_t jets[4][10]);
+ void RecJetsUA1(Int_t& njets, Float_t jets[4][50]);
+ void SetPycellParameters(Float_t etamax = 2., Int_t neta = 274, Int_t nphi = 432,
+ Float_t thresh = 0., Float_t etseed = 4.,
+ Float_t minet = 10., Float_t r = 1.);
+
+ // Getters
+ virtual Process_t GetProcess() const {return fProcess;}
+ virtual StrucFunc_t GetStrucFunc() const {return fStrucFunc;}
+ virtual void GetPtHard(Float_t& ptmin, Float_t& ptmax) const
+ {ptmin = fPtHardMin; ptmax = fPtHardMax;}
+ virtual Float_t GetEnergyCMS() const {return fEnergyCMS;}
+ virtual void GetNuclei(Int_t& a1, Int_t& a2) const
+ {a1 = fAProjectile; a2 = fATarget;}
+ virtual void GetJetEtRange(Float_t& etamin, Float_t& etamax) const
+ {etamin = fEtaMinJet; etamax = fEtaMaxJet;}
+ virtual void GetJetPhiRange(Float_t& phimin, Float_t& phimax) const
+ {phimin = fPhiMinJet*180./TMath::Pi(); phimax = fPhiMaxJet*180/TMath::Pi();}
+ virtual void GetGammaEtaRange(Float_t& etamin, Float_t& etamax) const
+ {etamin = fEtaMinGamma; etamax = fEtaMaxGamma;}
+ virtual void GetGammaPhiRange(Float_t& phimin, Float_t& phimax) const
+ {phimin = fPhiMinGamma*180./TMath::Pi(); phimax = fPhiMaxGamma*180./TMath::Pi();}
+ //
+ Bool_t IsInEMCAL(Float_t phi, Float_t eta);
+ Bool_t IsInPHOS(Float_t phi, Float_t eta);
+ //
+ virtual void FinishRun();
+ Bool_t CheckTrigger(TParticle* jet1, TParticle* jet2);
+ //Used in some processes to selected child properties
+ Bool_t CheckKinematicsOnChild();
+ void GetSubEventTime();
+
+ protected:
+ // adjust the weight from kinematic cuts
+ void AdjustWeights() const;
+ Int_t GenerateMB();
+ void MakeHeader();
+ void GeneratePileup();
+ AliPythiaBase *fPythia; //!Pythia
+ Process_t fProcess; //Process type
+ StrucFunc_t fStrucFunc; //Structure Function
+ Float_t fEnergyCMS; //Centre of mass energy
+ Float_t fKineBias; //!Bias from kinematic selection
+ Int_t fTrials; //!Number of trials for current event
+ Int_t fTrialsRun; //!Number of trials for run
+ Float_t fQ; //Mean Q
+ Float_t fX1; //Mean x1
+ Float_t fX2; //Mean x2
+ Float_t fEventTime; //Time of the subevent
+ Float_t fInteractionRate; //Interaction rate (set by user)
+ Float_t fTimeWindow; //Time window for pileup events (set by user)
+ Int_t fCurSubEvent; //Index of the current sub-event
+ TArrayF *fEventsTime; //Subevents time for pileup
+ Int_t fNev; //Number of events
+ Int_t fFlavorSelect; //Heavy Flavor Selection
+ Float_t fXsection; //Cross-section
+ Float_t fPtHardMin; //lower pT-hard cut
+ Float_t fPtHardMax; //higher pT-hard cut
+ Float_t fYHardMin; //lower y-hard cut
+ Float_t fYHardMax; //higher y-hard cut
+ Int_t fGinit; //initial state gluon radiation
+ Int_t fGfinal; //final state gluon radiation
+ Int_t fHadronisation; //hadronisation
+ Int_t fNpartons; //Number of partons before hadronisation
+ Int_t fReadFromFile; //read partons from file
+ Int_t fQuench; //Flag for quenching
+ Float_t fPtKick; //Transverse momentum kick
+ Bool_t fFullEvent; //!Write Full event if true
+ AliDecayer *fDecayer; //!Pointer to the decayer instance
+ Int_t fDebugEventFirst; //!First event to debug
+ Int_t fDebugEventLast; //!Last event to debug
+ Float_t fEtMinJet; //Minimum et of triggered Jet
+ Float_t fEtMaxJet; //Maximum et of triggered Jet
+ Float_t fEtaMinJet; //Minimum eta of triggered Jet
+ Float_t fEtaMaxJet; //Maximum eta of triggered Jet
+ Float_t fPhiMinJet; //Minimum phi of triggered Jet
+ Float_t fPhiMaxJet; //Maximum phi of triggered Jet
+ Int_t fJetReconstruction; //Jet Reconstruction mode
+ Float_t fEtaMinGamma; // Minimum eta of triggered gamma
+ Float_t fEtaMaxGamma; // Maximum eta of triggered gamma
+ Float_t fPhiMinGamma; // Minimum phi of triggered gamma
+ Float_t fPhiMaxGamma; // Maximum phi of triggered gamma
+ Float_t fPycellEtaMax; // Max. eta for Pycell
+ Int_t fPycellNEta; // Number of eta bins for Pycell
+ Int_t fPycellNPhi; // Number of phi bins for Pycell
+ Float_t fPycellThreshold; // Pycell threshold
+ Float_t fPycellEtSeed; // Pycell seed
+ Float_t fPycellMinEtJet; // Pycell min. jet et
+ Float_t fPycellMaxRadius; // Pycell cone radius
+ StackFillOpt_t fStackFillOpt; // Stack filling with all particles with
+ // that flavour or only with selected
+ // parents and their decays
+ Bool_t fFeedDownOpt; // Option to set feed down from higher
+ // quark families (e.g. b->c)
+ Bool_t fFragmentation; // Option to activate fragmentation by Pythia
+ Bool_t fSetNuclei; // Flag indicating that SetNuclei has been called
+ Bool_t fNewMIS; // Flag for the new multipple interactions scenario
+ Bool_t fHFoff; // Flag for switching heafy flavor production off
+ Int_t fTriggerParticle; // Trigger on this particle ...
+ Float_t fTriggerEta; // .. within |eta| < fTriggerEta
+ CountMode_t fCountMode; // Options for counting when the event will be finished.
+ // fCountMode = kCountAll --> All particles that end up in the
+ // stack are counted
+ // fCountMode = kCountParents --> Only selected parents are counted
+ // fCountMode = kCountTrackabless --> Only particles flagged for tracking
+ // are counted
+ //
+ //
+
+ AliGenPythiaEventHeader* fHeader; //! Event header
+ AliRunLoader* fRL; //! Run Loader
+ const Text_t* fFileName; //! Name of file to read from
+
+
+ Bool_t fFragPhotonInCalo; // Option to ask for Fragmentation Photon in calorimeters acceptance
+ Bool_t fPi0InCalo; // Option to ask for Pi0 in calorimeters acceptance
+ Bool_t fPhotonInCalo; // Option to ask for Decay Photon in calorimeter acceptance
+ Bool_t fCheckEMCAL; // Option to ask for FragPhoton or Pi0 in calorimeters EMCAL acceptance
+ Bool_t fCheckPHOS; // Option to ask for FragPhoton or Pi0 in calorimeters PHOS acceptance
+ Bool_t fCheckPHOSeta; // Option to ask for PHOS eta acceptance
+ Float_t fFragPhotonOrPi0MinPt; // Minimum momentum of Fragmentation Photon or Pi0
+ Float_t fPhotonMinPt; // Minimum momentum of Photon
+ //Calorimeters eta-phi acceptance
+ Float_t fPHOSMinPhi; // Minimum phi PHOS
+ Float_t fPHOSMaxPhi; // Maximum phi PHOS
+ Float_t fPHOSEta; // Minimum eta PHOS
+ Float_t fEMCALMinPhi; // Minimum phi EMCAL
+ Float_t fEMCALMaxPhi; // Maximum phi EMCAL
+ Float_t fEMCALEta; // Maximum eta EMCAL
+
+ private:
+ AliGenPythiaPlus(const AliGenPythiaPlus &Pythia);
+ AliGenPythiaPlus & operator=(const AliGenPythiaPlus & rhs);
+
+ ClassDef(AliGenPythiaPlus, 1) // AliGenerator interface to Pythia
+};
+#endif
+
+
+
+
+
--- /dev/null
+
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+/* $Id: AliPythia.cxx,v 1.40 2007/10/09 08:43:24 morsch Exp $ */
+
+#include "AliPythia6.h"
+#include "AliStack.h"
+#include "AliPythiaRndm.h"
+#include "AliFastGlauber.h"
+#include "AliQuenchingWeights.h"
+
+#include "TVector3.h"
+#include "TParticle.h"
+#include "PyquenCommon.h"
+
+ClassImp(AliPythia6)
+
+#ifndef WIN32
+# define pyclus pyclus_
+# define pycell pycell_
+# define pyshow pyshow_
+# define pyrobo pyrobo_
+# define pyquen pyquen_
+# define pyevnw pyevnw_
+# define type_of_call
+#else
+# define pyclus PYCLUS
+# define pycell PYCELL
+# define pyrobo PYROBO
+# define pyquen PYQUEN
+# define pyevnw PYEVNW
+# define type_of_call _stdcall
+#endif
+
+extern "C" void type_of_call pyclus(Int_t & );
+extern "C" void type_of_call pycell(Int_t & );
+extern "C" void type_of_call pyshow(Int_t &, Int_t &, Double_t &);
+extern "C" void type_of_call pyrobo(Int_t &, Int_t &, Double_t &, Double_t &, Double_t &, Double_t &, Double_t &);
+extern "C" void type_of_call pyquen(Double_t &, Int_t &, Double_t &);
+extern "C" void type_of_call pyevnw();
+
+
+//_____________________________________________________________________________
+
+AliPythia6* AliPythia6::fgAliPythia=NULL;
+
+AliPythia6::AliPythia6():
+ TPythia6(),
+ AliPythiaBase(),
+ fProcess(kPyMb),
+ fEcms(0.),
+ fStrucFunc(kCTEQ5L),
+ fXJet(0.),
+ fYJet(0.),
+ fNGmax(30),
+ fZmax(0.97),
+ fGlauber(0),
+ fQuenchingWeights(0)
+{
+// Default Constructor
+//
+// Set random number
+ if (!AliPythiaRndm::GetPythiaRandom())
+ AliPythiaRndm::SetPythiaRandom(GetRandom());
+ fGlauber = 0;
+ fQuenchingWeights = 0;
+}
+
+AliPythia6::AliPythia6(const AliPythia6& pythia):
+ TPythia6(),
+ AliPythiaBase(),
+ fProcess(kPyMb),
+ fEcms(0.),
+ fStrucFunc(kCTEQ5L),
+ fXJet(0.),
+ fYJet(0.),
+ fNGmax(30),
+ fZmax(0.97),
+ fGlauber(0),
+ fQuenchingWeights(0)
+{
+ // Copy Constructor
+ pythia.Copy(*this);
+}
+
+void AliPythia6::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfunc)
+{
+// Initialise the process to generate
+ if (!AliPythiaRndm::GetPythiaRandom())
+ AliPythiaRndm::SetPythiaRandom(GetRandom());
+
+ fProcess = process;
+ fEcms = energy;
+ fStrucFunc = strucfunc;
+//...Switch off decay of pi0, K0S, Lambda, Sigma+-, Xi0-, Omega-.
+ SetMDCY(Pycomp(111) ,1,0);
+ SetMDCY(Pycomp(310) ,1,0);
+ SetMDCY(Pycomp(3122),1,0);
+ SetMDCY(Pycomp(3112),1,0);
+ SetMDCY(Pycomp(3212),1,0);
+ SetMDCY(Pycomp(3222),1,0);
+ SetMDCY(Pycomp(3312),1,0);
+ SetMDCY(Pycomp(3322),1,0);
+ SetMDCY(Pycomp(3334),1,0);
+ // Select structure function
+ SetMSTP(52,2);
+ SetMSTP(51,AliStructFuncType::PDFsetIndex(strucfunc));
+ // Particles produced in string fragmentation point directly to either of the two endpoints
+ // of the string (depending in the side they were generated from).
+ SetMSTU(16,2);
+
+//
+// Pythia initialisation for selected processes//
+//
+// Make MSEL clean
+//
+ for (Int_t i=1; i<= 200; i++) {
+ SetMSUB(i,0);
+ }
+// select charm production
+ switch (process)
+ {
+ case kPyOldUEQ2ordered: //Old underlying events with Q2 ordered QCD processes
+// Multiple interactions on.
+ SetMSTP(81,1);
+// Double Gaussian matter distribution.
+ SetMSTP(82,4);
+ SetPARP(83,0.5);
+ SetPARP(84,0.4);
+// pT0.
+ SetPARP(82,2.0);
+// Reference energy for pT0 and energy rescaling pace.
+ SetPARP(89,1800);
+ SetPARP(90,0.25);
+// String drawing almost completely minimizes string length.
+ SetPARP(85,0.9);
+ SetPARP(86,0.95);
+// ISR and FSR activity.
+ SetPARP(67,4);
+ SetPARP(71,4);
+// Lambda_FSR scale.
+ SetPARJ(81,0.29);
+ break;
+ case kPyOldUEQ2ordered2:
+// Old underlying events with Q2 ordered QCD processes
+// Multiple interactions on.
+ SetMSTP(81,1);
+// Double Gaussian matter distribution.
+ SetMSTP(82,4);
+ SetPARP(83,0.5);
+ SetPARP(84,0.4);
+// pT0.
+ SetPARP(82,2.0);
+// Reference energy for pT0 and energy rescaling pace.
+ SetPARP(89,1800);
+ SetPARP(90,0.16); // here is the difference with kPyOldUEQ2ordered
+// String drawing almost completely minimizes string length.
+ SetPARP(85,0.9);
+ SetPARP(86,0.95);
+// ISR and FSR activity.
+ SetPARP(67,4);
+ SetPARP(71,4);
+// Lambda_FSR scale.
+ SetPARJ(81,0.29);
+ break;
+ case kPyOldPopcorn:
+// Old production mechanism: Old Popcorn
+ SetMSEL(1);
+ SetMSTJ(12,3);
+// (D=2) Like MSTJ(12)=2 but added prod ofthe 1er rank baryon
+ SetMSTP(88,2);
+// (D=1)see can be used to form baryons (BARYON JUNCTION)
+ SetMSTJ(1,1);
+ AtlasTuning();
+ break;
+ case kPyCharm:
+ SetMSEL(4);
+// heavy quark masses
+
+ SetPMAS(4,1,1.2);
+//
+// primordial pT
+ SetMSTP(91,1);
+ SetPARP(91,1.);
+ SetPARP(93,5.);
+//
+ break;
+ case kPyBeauty:
+ SetMSEL(5);
+ SetPMAS(5,1,4.75);
+ break;
+ case kPyJpsi:
+ SetMSEL(0);
+// gg->J/Psi g
+ SetMSUB(86,1);
+ break;
+ case kPyJpsiChi:
+ SetMSEL(0);
+// gg->J/Psi g
+ SetMSUB(86,1);
+// gg-> chi_0c g
+ SetMSUB(87,1);
+// gg-> chi_1c g
+ SetMSUB(88,1);
+// gg-> chi_2c g
+ SetMSUB(89,1);
+ break;
+ case kPyCharmUnforced:
+ SetMSEL(0);
+// gq->qg
+ SetMSUB(28,1);
+// gg->qq
+ SetMSUB(53,1);
+// gg->gg
+ SetMSUB(68,1);
+ break;
+ case kPyBeautyUnforced:
+ SetMSEL(0);
+// gq->qg
+ SetMSUB(28,1);
+// gg->qq
+ SetMSUB(53,1);
+// gg->gg
+ SetMSUB(68,1);
+ break;
+ case kPyMb:
+// Minimum Bias pp-Collisions
+//
+//
+// select Pythia min. bias model
+ SetMSEL(0);
+ SetMSUB(92,1); // single diffraction AB-->XB
+ SetMSUB(93,1); // single diffraction AB-->AX
+ SetMSUB(94,1); // double diffraction
+ SetMSUB(95,1); // low pt production
+
+ AtlasTuning();
+ break;
+ case kPyMbDefault:
+// Minimum Bias pp-Collisions
+//
+//
+// select Pythia min. bias model
+ SetMSEL(0);
+ SetMSUB(92,1); // single diffraction AB-->XB
+ SetMSUB(93,1); // single diffraction AB-->AX
+ SetMSUB(94,1); // double diffraction
+ SetMSUB(95,1); // low pt production
+
+ break;
+ case kPyLhwgMb:
+// Les Houches Working Group 05 Minimum Bias pp-Collisions: hep-ph/0604120
+// -> Pythia 6.3 or above is needed
+//
+ SetMSEL(0);
+ SetMSUB(92,1); // single diffraction AB-->XB
+ SetMSUB(93,1); // single diffraction AB-->AX
+ SetMSUB(94,1); // double diffraction
+ SetMSUB(95,1); // low pt production
+ SetMSTP(51,AliStructFuncType::PDFsetIndex(kCTEQ6ll));
+ SetMSTP(52,2);
+ SetMSTP(68,1);
+ SetMSTP(70,2);
+ SetMSTP(81,1); // Multiple Interactions ON
+ SetMSTP(82,4); // Double Gaussian Model
+ SetMSTP(88,1);
+
+ SetPARP(82,2.3); // [GeV] PT_min at Ref. energy
+ SetPARP(83,0.5); // Core density in proton matter distribution (def.value)
+ SetPARP(84,0.5); // Core radius
+ SetPARP(85,0.9); // Regulates gluon prod. mechanism
+ SetPARP(90,0.2); // 2*epsilon (exponent in power law)
+
+ break;
+ case kPyMbNonDiffr:
+// Minimum Bias pp-Collisions
+//
+//
+// select Pythia min. bias model
+ SetMSEL(0);
+ SetMSUB(95,1); // low pt production
+
+ AtlasTuning();
+ break;
+ case kPyMbMSEL1:
+ ConfigHeavyFlavor();
+// Intrinsic <kT^2>
+ SetMSTP(91,1);// Width (1=gaussian) primordial kT dist. inside hadrons
+ SetPARP(91,1.); // <kT^2> = PARP(91,1.)^2
+ SetPARP(93,5.); // Upper cut-off
+// Set Q-quark mass
+ SetPMAS(4,1,1.2); // Charm quark mass
+ SetPMAS(5,1,4.78); // Beauty quark mass
+ SetPARP(71,4.); // Defaut value
+// Atlas Tuning
+ AtlasTuning();
+ break;
+ case kPyJets:
+//
+// QCD Jets
+//
+ SetMSEL(1);
+ // Pythia Tune A (CDF)
+ //
+ SetPARP(67,2.5); // Regulates Initial State Radiation (value from best fit to D0 dijet analysis)
+ SetMSTP(82,4); // Double Gaussian Model
+ SetPARP(82,2.0); // [GeV] PT_min at Ref. energy
+ SetPARP(84,0.4); // Core radius
+ SetPARP(85,0.90) ; // Regulates gluon prod. mechanism
+ SetPARP(86,0.95); // Regulates gluon prod. mechanism
+ SetPARP(89,1800.); // [GeV] Ref. energy
+ SetPARP(90,0.25); // 2*epsilon (exponent in power law)
+ break;
+ case kPyDirectGamma:
+ SetMSEL(10);
+ break;
+ case kPyCharmPbPbMNR:
+ case kPyD0PbPbMNR:
+ case kPyDPlusPbPbMNR:
+ case kPyDPlusStrangePbPbMNR:
+ // Tuning of Pythia parameters aimed to get a resonable agreement
+ // between with the NLO calculation by Mangano, Nason, Ridolfi for the
+ // c-cbar single inclusive and double differential distributions.
+ // This parameter settings are meant to work with Pb-Pb collisions
+ // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
+ // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
+ // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
+ ConfigHeavyFlavor();
+ // Intrinsic <kT>
+ SetMSTP(91,1);
+ SetPARP(91,1.304);
+ SetPARP(93,6.52);
+ // Set c-quark mass
+ SetPMAS(4,1,1.2);
+ break;
+ case kPyCharmpPbMNR:
+ case kPyD0pPbMNR:
+ case kPyDPluspPbMNR:
+ case kPyDPlusStrangepPbMNR:
+ // Tuning of Pythia parameters aimed to get a resonable agreement
+ // between with the NLO calculation by Mangano, Nason, Ridolfi for the
+ // c-cbar single inclusive and double differential distributions.
+ // This parameter settings are meant to work with p-Pb collisions
+ // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
+ // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
+ // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
+ ConfigHeavyFlavor();
+ // Intrinsic <kT>
+ SetMSTP(91,1);
+ SetPARP(91,1.16);
+ SetPARP(93,5.8);
+
+ // Set c-quark mass
+ SetPMAS(4,1,1.2);
+ break;
+ case kPyCharmppMNR:
+ case kPyD0ppMNR:
+ case kPyDPlusppMNR:
+ case kPyDPlusStrangeppMNR:
+ // Tuning of Pythia parameters aimed to get a resonable agreement
+ // between with the NLO calculation by Mangano, Nason, Ridolfi for the
+ // c-cbar single inclusive and double differential distributions.
+ // This parameter settings are meant to work with pp collisions
+ // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
+ // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
+ // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
+ ConfigHeavyFlavor();
+ // Intrinsic <kT^2>
+ SetMSTP(91,1);
+ SetPARP(91,1.);
+ SetPARP(93,5.);
+
+ // Set c-quark mass
+ SetPMAS(4,1,1.2);
+ break;
+ case kPyCharmppMNRwmi:
+ // Tuning of Pythia parameters aimed to get a resonable agreement
+ // between with the NLO calculation by Mangano, Nason, Ridolfi for the
+ // c-cbar single inclusive and double differential distributions.
+ // This parameter settings are meant to work with pp collisions
+ // and with kCTEQ5L PDFs.
+ // Added multiple interactions according to ATLAS tune settings.
+ // To get a "reasonable" agreement with MNR results, events have to be
+ // generated with the minimum ptHard (AliGenPythia::SetPtHard)
+ // set to 2.76 GeV.
+ // To get a "perfect" agreement with MNR results, events have to be
+ // generated in four ptHard bins with the following relative
+ // normalizations:
+ // 2.76-3 GeV: 25%
+ // 3-4 GeV: 40%
+ // 4-8 GeV: 29%
+ // >8 GeV: 6%
+ ConfigHeavyFlavor();
+ // Intrinsic <kT^2>
+ SetMSTP(91,1);
+ SetPARP(91,1.);
+ SetPARP(93,5.);
+
+ // Set c-quark mass
+ SetPMAS(4,1,1.2);
+ AtlasTuning();
+ break;
+ case kPyBeautyPbPbMNR:
+ // Tuning of Pythia parameters aimed to get a resonable agreement
+ // between with the NLO calculation by Mangano, Nason, Ridolfi for the
+ // b-bbar single inclusive and double differential distributions.
+ // This parameter settings are meant to work with Pb-Pb collisions
+ // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
+ // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
+ // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
+ ConfigHeavyFlavor();
+ // QCD scales
+ SetPARP(67,1.0);
+ SetPARP(71,1.0);
+ // Intrinsic <kT>
+ SetMSTP(91,1);
+ SetPARP(91,2.035);
+ SetPARP(93,10.17);
+ // Set b-quark mass
+ SetPMAS(5,1,4.75);
+ break;
+ case kPyBeautypPbMNR:
+ // Tuning of Pythia parameters aimed to get a resonable agreement
+ // between with the NLO calculation by Mangano, Nason, Ridolfi for the
+ // b-bbar single inclusive and double differential distributions.
+ // This parameter settings are meant to work with p-Pb collisions
+ // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
+ // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
+ // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
+ ConfigHeavyFlavor();
+ // QCD scales
+ SetPARP(67,1.0);
+ SetPARP(71,1.0);
+ // Intrinsic <kT>
+ SetMSTP(91,1);
+ SetPARP(91,1.60);
+ SetPARP(93,8.00);
+ // Set b-quark mass
+ SetPMAS(5,1,4.75);
+ break;
+ case kPyBeautyppMNR:
+ // Tuning of Pythia parameters aimed to get a resonable agreement
+ // between with the NLO calculation by Mangano, Nason, Ridolfi for the
+ // b-bbar single inclusive and double differential distributions.
+ // This parameter settings are meant to work with pp collisions
+ // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
+ // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
+ // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
+ ConfigHeavyFlavor();
+ // QCD scales
+ SetPARP(67,1.0);
+ SetPARP(71,1.0);
+
+ // Intrinsic <kT>
+ SetMSTP(91,1);
+ SetPARP(91,1.);
+ SetPARP(93,5.);
+
+ // Set b-quark mass
+ SetPMAS(5,1,4.75);
+ break;
+ case kPyBeautyppMNRwmi:
+ // Tuning of Pythia parameters aimed to get a resonable agreement
+ // between with the NLO calculation by Mangano, Nason, Ridolfi for the
+ // b-bbar single inclusive and double differential distributions.
+ // This parameter settings are meant to work with pp collisions
+ // and with kCTEQ5L PDFs.
+ // Added multiple interactions according to ATLAS tune settings.
+ // To get a "reasonable" agreement with MNR results, events have to be
+ // generated with the minimum ptHard (AliGenPythia::SetPtHard)
+ // set to 2.76 GeV.
+ // To get a "perfect" agreement with MNR results, events have to be
+ // generated in four ptHard bins with the following relative
+ // normalizations:
+ // 2.76-4 GeV: 5%
+ // 4-6 GeV: 31%
+ // 6-8 GeV: 28%
+ // >8 GeV: 36%
+ ConfigHeavyFlavor();
+ // QCD scales
+ SetPARP(67,1.0);
+ SetPARP(71,1.0);
+
+ // Intrinsic <kT>
+ SetMSTP(91,1);
+ SetPARP(91,1.);
+ SetPARP(93,5.);
+
+ // Set b-quark mass
+ SetPMAS(5,1,4.75);
+
+ AtlasTuning();
+ break;
+ case kPyW:
+
+ //Inclusive production of W+/-
+ SetMSEL(0);
+ //f fbar -> W+
+ SetMSUB(2,1);
+ // //f fbar -> g W+
+ // SetMSUB(16,1);
+ // //f fbar -> gamma W+
+ // SetMSUB(20,1);
+ // //f g -> f W+
+ // SetMSUB(31,1);
+ // //f gamma -> f W+
+ // SetMSUB(36,1);
+
+ // Initial/final parton shower on (Pythia default)
+ // With parton showers on we are generating "W inclusive process"
+ SetMSTP(61,1); //Initial QCD & QED showers on
+ SetMSTP(71,1); //Final QCD & QED showers on
+
+ break;
+
+ case kPyZ:
+
+ //Inclusive production of Z
+ SetMSEL(0);
+ //f fbar -> Z/gamma
+ SetMSUB(1,1);
+
+ // // f fbar -> g Z/gamma
+ // SetMSUB(15,1);
+ // // f fbar -> gamma Z/gamma
+ // SetMSUB(19,1);
+ // // f g -> f Z/gamma
+ // SetMSUB(30,1);
+ // // f gamma -> f Z/gamma
+ // SetMSUB(35,1);
+
+ //only Z included, not gamma
+ SetMSTP(43,2);
+
+ // Initial/final parton shower on (Pythia default)
+ // With parton showers on we are generating "Z inclusive process"
+ SetMSTP(61,1); //Initial QCD & QED showers on
+ SetMSTP(71,1); //Final QCD & QED showers on
+
+ break;
+
+ }
+//
+// Initialize PYTHIA
+ SetMSTP(41,1); // all resonance decays switched on
+ Initialize("CMS","p","p",fEcms);
+
+}
+
+Int_t AliPythia6::CheckedLuComp(Int_t kf)
+{
+// Check Lund particle code (for debugging)
+ Int_t kc=Pycomp(kf);
+ return kc;
+}
+
+void AliPythia6::SetNuclei(Int_t a1, Int_t a2)
+{
+// Treat protons as inside nuclei with mass numbers a1 and a2
+// The MSTP array in the PYPARS common block is used to enable and
+// select the nuclear structure functions.
+// MSTP(52) : (D=1) choice of proton and nuclear structure-function library
+// =1: internal PYTHIA acording to MSTP(51)
+// =2: PDFLIB proton s.f., with MSTP(51) = 1000xNGROUP+NSET
+// If the following mass number both not equal zero, nuclear corrections of the stf are used.
+// MSTP(192) : Mass number of nucleus side 1
+// MSTP(193) : Mass number of nucleus side 2
+ SetMSTP(52,2);
+ SetMSTP(192, a1);
+ SetMSTP(193, a2);
+}
+
+
+AliPythia6* AliPythia6::Instance()
+{
+// Set random number generator
+ if (fgAliPythia) {
+ return fgAliPythia;
+ } else {
+ fgAliPythia = new AliPythia6();
+ return fgAliPythia;
+ }
+}
+
+void AliPythia6::PrintParticles()
+{
+// Print list of particl properties
+ Int_t np = 0;
+ char* name = new char[16];
+ for (Int_t kf=0; kf<1000000; kf++) {
+ for (Int_t c = 1; c > -2; c-=2) {
+ Int_t kc = Pycomp(c*kf);
+ if (kc) {
+ Float_t mass = GetPMAS(kc,1);
+ Float_t width = GetPMAS(kc,2);
+ Float_t tau = GetPMAS(kc,4);
+
+ Pyname(kf,name);
+
+ np++;
+
+ printf("\n mass, width, tau: %6d %s %10.3f %10.3e %10.3e",
+ c*kf, name, mass, width, tau);
+ }
+ }
+ }
+ printf("\n Number of particles %d \n \n", np);
+}
+
+void AliPythia6::ResetDecayTable()
+{
+// Set default values for pythia decay switches
+ Int_t i;
+ for (i = 1; i < 501; i++) SetMDCY(i,1,fDefMDCY[i]);
+ for (i = 1; i < 2001; i++) SetMDME(i,1,fDefMDME[i]);
+}
+
+void AliPythia6::SetDecayTable()
+{
+// Set default values for pythia decay switches
+//
+ Int_t i;
+ for (i = 1; i < 501; i++) fDefMDCY[i] = GetMDCY(i,1);
+ for (i = 1; i < 2001; i++) fDefMDME[i] = GetMDME(i,1);
+}
+
+void AliPythia6::Pyclus(Int_t& njet)
+{
+// Call Pythia clustering algorithm
+//
+ pyclus(njet);
+}
+
+void AliPythia6::Pycell(Int_t& njet)
+{
+// Call Pythia jet reconstruction algorithm
+//
+ pycell(njet);
+}
+
+void AliPythia6::GetJet(Int_t i, Float_t& px, Float_t& py, Float_t& pz, Float_t& e)
+{
+ // Get jet number i
+ Int_t n = GetN();
+ px = GetPyjets()->P[0][n+i];
+ py = GetPyjets()->P[1][n+i];
+ pz = GetPyjets()->P[2][n+i];
+ e = GetPyjets()->P[3][n+i];
+}
+
+void AliPythia6::Pyshow(Int_t ip1, Int_t ip2, Double_t qmax)
+{
+// Call Pythia showering
+//
+ pyshow(ip1, ip2, qmax);
+}
+
+void AliPythia6::Pyrobo(Int_t imi, Int_t ima, Double_t the, Double_t phi, Double_t bex, Double_t bey, Double_t bez)
+{
+ pyrobo(imi, ima, the, phi, bex, bey, bez);
+}
+
+
+
+void AliPythia6::InitQuenching(Float_t cMin, Float_t cMax, Float_t k, Int_t iECMethod, Float_t zmax, Int_t ngmax)
+{
+// Initializes
+// (1) The quenching model using quenching weights according to C. Salgado and U. Wiedemann
+// (2) The nuclear geometry using the Glauber Model
+//
+
+ fGlauber = new AliFastGlauber();
+ fGlauber->Init(2);
+ fGlauber->SetCentralityClass(cMin, cMax);
+
+ fQuenchingWeights = new AliQuenchingWeights();
+ fQuenchingWeights->InitMult();
+ fQuenchingWeights->SetK(k);
+ fQuenchingWeights->SetECMethod(AliQuenchingWeights::kECMethod(iECMethod));
+ fNGmax = ngmax;
+ fZmax = zmax;
+
+}
+
+
+void AliPythia6::Quench()
+{
+//
+//
+// Simple Jet Quenching routine:
+// =============================
+// The jet formed by all final state partons radiated by the parton created
+// in the hard collisions is quenched by a factor (1-z) using light cone variables in
+// the initial parton reference frame:
+// (E + p_z)new = (1-z) (E + p_z)old
+//
+//
+//
+//
+// The lost momentum is first balanced by one gluon with virtuality > 0.
+// Subsequently the gluon splits to yield two gluons with E = p.
+//
+//
+//
+ static Float_t eMean = 0.;
+ static Int_t icall = 0;
+
+ Double_t p0[4][5];
+ Double_t p1[4][5];
+ Double_t p2[4][5];
+ Int_t klast[4] = {-1, -1, -1, -1};
+
+ Int_t numpart = fPyjets->N;
+ Double_t px = 0., py = 0., pz = 0., e = 0., m = 0., p = 0., pt = 0., theta = 0., phi = 0.;
+ Double_t pxq[4], pyq[4], pzq[4], eq[4], yq[4], mq[4], pq[4], phiq[4], thetaq[4], ptq[4];
+ Bool_t quenched[4];
+ Double_t wjtKick[4];
+ Int_t nGluon[4];
+ Int_t qPdg[4];
+ Int_t imo, kst, pdg;
+
+//
+// Sore information about Primary partons
+//
+// j =
+// 0, 1 partons from hard scattering
+// 2, 3 partons from initial state radiation
+//
+ for (Int_t i = 2; i <= 7; i++) {
+ Int_t j = 0;
+ // Skip gluons that participate in hard scattering
+ if (i == 4 || i == 5) continue;
+ // Gluons from hard Scattering
+ if (i == 6 || i == 7) {
+ j = i - 6;
+ pxq[j] = fPyjets->P[0][i];
+ pyq[j] = fPyjets->P[1][i];
+ pzq[j] = fPyjets->P[2][i];
+ eq[j] = fPyjets->P[3][i];
+ mq[j] = fPyjets->P[4][i];
+ } else {
+ // Gluons from initial state radiation
+ //
+ // Obtain 4-momentum vector from difference between original parton and parton after gluon
+ // radiation. Energy is calculated independently because initial state radition does not
+ // conserve strictly momentum and energy for each partonic system independently.
+ //
+ // Not very clean. Should be improved !
+ //
+ //
+ j = i;
+ pxq[j] = fPyjets->P[0][i] - fPyjets->P[0][i+2];
+ pyq[j] = fPyjets->P[1][i] - fPyjets->P[1][i+2];
+ pzq[j] = fPyjets->P[2][i] - fPyjets->P[2][i+2];
+ mq[j] = fPyjets->P[4][i];
+ eq[j] = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j] + pzq[j] * pzq[j] + mq[j] * mq[j]);
+ }
+//
+// Calculate some kinematic variables
+//
+ yq[j] = 0.5 * TMath::Log((eq[j] + pzq[j] + 1.e-14) / (eq[j] - pzq[j] + 1.e-14));
+ pq[j] = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j] + pzq[j] * pzq[j]);
+ phiq[j] = TMath::Pi()+TMath::ATan2(-pyq[j], -pxq[j]);
+ ptq[j] = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j]);
+ thetaq[j] = TMath::ATan2(ptq[j], pzq[j]);
+ qPdg[j] = fPyjets->K[1][i];
+ }
+
+ Double_t int0[4];
+ Double_t int1[4];
+
+ fGlauber->GetI0I1ForPythiaAndXY(4, phiq, int0, int1, fXJet, fYJet, 15.);
+
+ for (Int_t j = 0; j < 4; j++) {
+ //
+ // Quench only central jets and with E > 10.
+ //
+
+
+ Int_t itype = (qPdg[j] == 21) ? 2 : 1;
+ Double_t eloss = fQuenchingWeights->GetELossRandomKFast(itype, int0[j], int1[j], eq[j]);
+
+ if (TMath::Abs(yq[j]) > 2.5 || eq[j] < 10.) {
+ fZQuench[j] = 0.;
+ } else {
+ if (eq[j] > 40. && TMath::Abs(yq[j]) < 0.5) {
+ icall ++;
+ eMean += eloss;
+ }
+ //
+ // Extra pt
+ Double_t l = fQuenchingWeights->CalcLk(int0[j], int1[j]);
+ wjtKick[j] = TMath::Sqrt(l * fQuenchingWeights->CalcQk(int0[j], int1[j]));
+ //
+ // Fractional energy loss
+ fZQuench[j] = eloss / eq[j];
+ //
+ // Avoid complete loss
+ //
+ if (fZQuench[j] == 1.) fZQuench[j] = fZmax;
+ //
+ // Some debug printing
+
+
+// printf("Initial parton # %3d, Type %3d Energy %10.3f Phi %10.3f Length %10.3f Loss %10.3f Kick %10.3f Mean: %10.3f %10.3f\n",
+// j, itype, eq[j], phiq[j], l, eloss, wjtKick[j], eMean / Float_t(icall+1), yq[j]);
+
+// fZQuench[j] = 0.8;
+// while (fZQuench[j] >= 0.95) fZQuench[j] = gRandom->Exp(0.2);
+ }
+
+ quenched[j] = (fZQuench[j] > 0.01);
+ } // primary partons
+
+
+
+ Double_t pNew[1000][4];
+ Int_t kNew[1000];
+ Int_t icount = 0;
+ Double_t zquench[4];
+
+//
+// System Loop
+ for (Int_t isys = 0; isys < 4; isys++) {
+// Skip to next system if not quenched.
+ if (!quenched[isys]) continue;
+
+ nGluon[isys] = 1 + Int_t(fZQuench[isys] / (1. - fZQuench[isys]));
+ if (nGluon[isys] > fNGmax) nGluon[isys] = fNGmax;
+ zquench[isys] = 1. - TMath::Power(1. - fZQuench[isys], 1./Double_t(nGluon[isys]));
+ wjtKick[isys] = wjtKick[isys] / TMath::Sqrt(Double_t(nGluon[isys]));
+
+
+
+ Int_t igMin = -1;
+ Int_t igMax = -1;
+ Double_t pg[4] = {0., 0., 0., 0.};
+
+//
+// Loop on radiation events
+
+ for (Int_t iglu = 0; iglu < nGluon[isys]; iglu++) {
+ while (1) {
+ icount = 0;
+ for (Int_t k = 0; k < 4; k++)
+ {
+ p0[isys][k] = 0.;
+ p1[isys][k] = 0.;
+ p2[isys][k] = 0.;
+ }
+// Loop over partons
+ for (Int_t i = 0; i < numpart; i++)
+ {
+ imo = fPyjets->K[2][i];
+ kst = fPyjets->K[0][i];
+ pdg = fPyjets->K[1][i];
+
+
+
+// Quarks and gluons only
+ if (pdg != 21 && TMath::Abs(pdg) > 6) continue;
+// Particles from hard scattering only
+
+ if (imo > 8 && imo < 1000) imo = fPyjets->K[2][imo - 1];
+ Int_t imom = imo % 1000;
+ if ((isys == 0 || isys == 1) && ((imom != (isys + 7)))) continue;
+ if ((isys == 2 || isys == 3) && ((imom != (isys + 1)))) continue;
+
+
+// Skip comment lines
+ if (kst != 1 && kst != 2) continue;
+//
+// Parton kinematic
+ px = fPyjets->P[0][i];
+ py = fPyjets->P[1][i];
+ pz = fPyjets->P[2][i];
+ e = fPyjets->P[3][i];
+ m = fPyjets->P[4][i];
+ pt = TMath::Sqrt(px * px + py * py);
+ p = TMath::Sqrt(px * px + py * py + pz * pz);
+ phi = TMath::Pi() + TMath::ATan2(-py, -px);
+ theta = TMath::ATan2(pt, pz);
+
+//
+// Save 4-momentum sum for balancing
+ Int_t index = isys;
+
+ p0[index][0] += px;
+ p0[index][1] += py;
+ p0[index][2] += pz;
+ p0[index][3] += e;
+
+ klast[index] = i;
+
+//
+// Fractional energy loss
+ Double_t z = zquench[index];
+
+
+// Don't fully quench radiated gluons
+//
+ if (imo > 1000) {
+// This small factor makes sure that the gluons are not too close in phase space to avoid recombination
+//
+
+ z = 0.02;
+ }
+// printf("z: %d %f\n", imo, z);
+
+
+//
+
+ //
+ //
+ // Transform into frame in which initial parton is along z-axis
+ //
+ TVector3 v(px, py, pz);
+ v.RotateZ(-phiq[index]); v.RotateY(-thetaq[index]);
+ Double_t pxs = v.X(); Double_t pys = v.Y(); Double_t pl = v.Z();
+
+ Double_t jt = TMath::Sqrt(pxs * pxs + pys * pys);
+ Double_t mt2 = jt * jt + m * m;
+ Double_t zmax = 1.;
+ //
+ // Kinematic limit on z
+ //
+ if (m > 0.) zmax = 1. - m / TMath::Sqrt(m * m + jt * jt);
+ //
+ // Change light-cone kinematics rel. to initial parton
+ //
+ Double_t eppzOld = e + pl;
+ Double_t empzOld = e - pl;
+
+ Double_t eppzNew = (1. - z) * eppzOld;
+ Double_t empzNew = empzOld - mt2 * z / eppzOld;
+ Double_t eNew = 0.5 * (eppzNew + empzNew);
+ Double_t plNew = 0.5 * (eppzNew - empzNew);
+
+ Double_t jtNew;
+ //
+ // if mt very small (or sometimes even < 0 for numerical reasons) set it to 0
+ Double_t mt2New = eppzNew * empzNew;
+ if (mt2New < 1.e-8) mt2New = 0.;
+ if (z < zmax) {
+ if (m * m > mt2New) {
+ //
+ // This should not happen
+ //
+ Fatal("Quench()", "This should never happen %e %e %e!", m, eppzNew, empzNew);
+ jtNew = 0;
+ } else {
+ jtNew = TMath::Sqrt(mt2New - m * m);
+ }
+ } else {
+ // If pT is to small (probably a leading massive particle) we scale only the energy
+ // This can cause negative masses of the radiated gluon
+ // Let's hope for the best ...
+ jtNew = jt;
+ eNew = TMath::Sqrt(plNew * plNew + mt2);
+
+ }
+ //
+ // Calculate new px, py
+ //
+ Double_t pxNew = jtNew / jt * pxs;
+ Double_t pyNew = jtNew / jt * pys;
+
+// Double_t dpx = pxs - pxNew;
+// Double_t dpy = pys - pyNew;
+// Double_t dpz = pl - plNew;
+// Double_t de = e - eNew;
+// Double_t dmass2 = de * de - dpx * dpx - dpy * dpy - dpz * dpz;
+// printf("New mass (1) %e %e %e %e %e %e %e \n", dmass2, jt, jtNew, pl, plNew, e, eNew);
+// printf("New mass (2) %e %e \n", pxNew, pyNew);
+ //
+ // Rotate back
+ //
+ TVector3 w(pxNew, pyNew, plNew);
+ w.RotateY(thetaq[index]); w.RotateZ(phiq[index]);
+ pxNew = w.X(); pyNew = w.Y(); plNew = w.Z();
+
+ p1[index][0] += pxNew;
+ p1[index][1] += pyNew;
+ p1[index][2] += plNew;
+ p1[index][3] += eNew;
+ //
+ // Updated 4-momentum vectors
+ //
+ pNew[icount][0] = pxNew;
+ pNew[icount][1] = pyNew;
+ pNew[icount][2] = plNew;
+ pNew[icount][3] = eNew;
+ kNew[icount] = i;
+ icount++;
+ } // parton loop
+ //
+ // Check if there was phase-space for quenching
+ //
+
+ if (icount == 0) quenched[isys] = kFALSE;
+ if (!quenched[isys]) break;
+
+ for (Int_t j = 0; j < 4; j++)
+ {
+ p2[isys][j] = p0[isys][j] - p1[isys][j];
+ }
+ p2[isys][4] = p2[isys][3] * p2[isys][3] - p2[isys][0] * p2[isys][0] - p2[isys][1] * p2[isys][1] - p2[isys][2] * p2[isys][2];
+ if (p2[isys][4] > 0.) {
+ p2[isys][4] = TMath::Sqrt(p2[isys][4]);
+ break;
+ } else {
+ printf("Warning negative mass squared in system %d %f ! \n", isys, zquench[isys]);
+ printf("4-Momentum: %10.3e %10.3e %10.3e %10.3e %10.3e \n", p2[isys][0], p2[isys][1], p2[isys][2], p2[isys][3], p2[isys][4]);
+ if (p2[isys][4] < -0.01) {
+ printf("Negative mass squared !\n");
+ // Here we have to put the gluon back to mass shell
+ // This will lead to a small energy imbalance
+ p2[isys][4] = 0.;
+ p2[isys][3] = TMath::Sqrt(p2[isys][0] * p2[isys][0] + p2[isys][1] * p2[isys][1] + p2[isys][2] * p2[isys][2]);
+ break;
+ } else {
+ p2[isys][4] = 0.;
+ break;
+ }
+ }
+ /*
+ zHeavy *= 0.98;
+ printf("zHeavy lowered to %f\n", zHeavy);
+ if (zHeavy < 0.01) {
+ printf("No success ! \n");
+ icount = 0;
+ quenched[isys] = kFALSE;
+ break;
+ }
+ */
+ } // iteration on z (while)
+
+// Update event record
+ for (Int_t k = 0; k < icount; k++) {
+// printf("%6d %6d %10.3e %10.3e %10.3e %10.3e\n", k, kNew[k], pNew[k][0],pNew[k][1], pNew[k][2], pNew[k][3] );
+ fPyjets->P[0][kNew[k]] = pNew[k][0];
+ fPyjets->P[1][kNew[k]] = pNew[k][1];
+ fPyjets->P[2][kNew[k]] = pNew[k][2];
+ fPyjets->P[3][kNew[k]] = pNew[k][3];
+ }
+ //
+ // Add the gluons
+ //
+ Int_t ish = 0;
+ Int_t iGlu;
+ if (!quenched[isys]) continue;
+//
+// Last parton from shower i
+ Int_t in = klast[isys];
+//
+// Continue if no parton in shower i selected
+ if (in == -1) continue;
+//
+// If this is the second initial parton and it is behind the first move pointer by previous ish
+ if (isys == 1 && klast[1] > klast[0]) in += ish;
+//
+// Starting index
+
+// jmin = in - 1;
+// How many additional gluons will be generated
+ ish = 1;
+ if (p2[isys][4] > 0.05) ish = 2;
+//
+// Position of gluons
+ iGlu = numpart;
+ if (iglu == 0) igMin = iGlu;
+ igMax = iGlu;
+ numpart += ish;
+ (fPyjets->N) += ish;
+
+ if (ish == 1) {
+ fPyjets->P[0][iGlu] = p2[isys][0];
+ fPyjets->P[1][iGlu] = p2[isys][1];
+ fPyjets->P[2][iGlu] = p2[isys][2];
+ fPyjets->P[3][iGlu] = p2[isys][3];
+ fPyjets->P[4][iGlu] = p2[isys][4];
+
+ fPyjets->K[0][iGlu] = 1;
+ if (iglu == nGluon[isys] - 1) fPyjets->K[0][iGlu] = 1;
+ fPyjets->K[1][iGlu] = 21;
+ fPyjets->K[2][iGlu] = fPyjets->K[2][in] + 1000;
+ fPyjets->K[3][iGlu] = -1;
+ fPyjets->K[4][iGlu] = -1;
+
+ pg[0] += p2[isys][0];
+ pg[1] += p2[isys][1];
+ pg[2] += p2[isys][2];
+ pg[3] += p2[isys][3];
+ } else {
+ //
+ // Split gluon in rest frame.
+ //
+ Double_t bx = p2[isys][0] / p2[isys][3];
+ Double_t by = p2[isys][1] / p2[isys][3];
+ Double_t bz = p2[isys][2] / p2[isys][3];
+ Double_t pst = p2[isys][4] / 2.;
+ //
+ // Isotropic decay ????
+ Double_t cost = 2. * gRandom->Rndm() - 1.;
+ Double_t sint = TMath::Sqrt(1. - cost * cost);
+ Double_t phi = 2. * TMath::Pi() * gRandom->Rndm();
+
+ Double_t pz1 = pst * cost;
+ Double_t pz2 = -pst * cost;
+ Double_t pt1 = pst * sint;
+ Double_t pt2 = -pst * sint;
+ Double_t px1 = pt1 * TMath::Cos(phi);
+ Double_t py1 = pt1 * TMath::Sin(phi);
+ Double_t px2 = pt2 * TMath::Cos(phi);
+ Double_t py2 = pt2 * TMath::Sin(phi);
+
+ fPyjets->P[0][iGlu] = px1;
+ fPyjets->P[1][iGlu] = py1;
+ fPyjets->P[2][iGlu] = pz1;
+ fPyjets->P[3][iGlu] = pst;
+ fPyjets->P[4][iGlu] = 0.;
+
+ fPyjets->K[0][iGlu] = 1 ;
+ fPyjets->K[1][iGlu] = 21;
+ fPyjets->K[2][iGlu] = fPyjets->K[2][in] + 1000;
+ fPyjets->K[3][iGlu] = -1;
+ fPyjets->K[4][iGlu] = -1;
+
+ fPyjets->P[0][iGlu+1] = px2;
+ fPyjets->P[1][iGlu+1] = py2;
+ fPyjets->P[2][iGlu+1] = pz2;
+ fPyjets->P[3][iGlu+1] = pst;
+ fPyjets->P[4][iGlu+1] = 0.;
+
+ fPyjets->K[0][iGlu+1] = 1;
+ if (iglu == nGluon[isys] - 1) fPyjets->K[0][iGlu+1] = 1;
+ fPyjets->K[1][iGlu+1] = 21;
+ fPyjets->K[2][iGlu+1] = fPyjets->K[2][in] + 1000;
+ fPyjets->K[3][iGlu+1] = -1;
+ fPyjets->K[4][iGlu+1] = -1;
+ SetMSTU(1,0);
+ SetMSTU(2,0);
+ //
+ // Boost back
+ //
+ Pyrobo(iGlu + 1, iGlu + 2, 0., 0., bx, by, bz);
+ }
+/*
+ for (Int_t ig = iGlu; ig < iGlu+ish; ig++) {
+ Double_t px, py, pz;
+ px = fPyjets->P[0][ig];
+ py = fPyjets->P[1][ig];
+ pz = fPyjets->P[2][ig];
+ TVector3 v(px, py, pz);
+ v.RotateZ(-phiq[isys]);
+ v.RotateY(-thetaq[isys]);
+ Double_t pxs = v.X(); Double_t pys = v.Y(); Double_t pzs = v.Z();
+ Double_t r = AliPythiaRndm::GetPythiaRandom()->Rndm();
+ Double_t jtKick = 0.3 * TMath::Sqrt(-TMath::Log(r));
+ if (ish == 2) jtKick = wjtKick[i] * TMath::Sqrt(-TMath::Log(r)) / TMath::Sqrt(2.);
+ Double_t phiKick = 2. * TMath::Pi() * AliPythiaRndm::GetPythiaRandom()->Rndm();
+ pxs += jtKick * TMath::Cos(phiKick);
+ pys += jtKick * TMath::Sin(phiKick);
+ TVector3 w(pxs, pys, pzs);
+ w.RotateY(thetaq[isys]);
+ w.RotateZ(phiq[isys]);
+ fPyjets->P[0][ig] = w.X();
+ fPyjets->P[1][ig] = w.Y();
+ fPyjets->P[2][ig] = w.Z();
+ fPyjets->P[2][ig] = w.Mag();
+ }
+*/
+ } // kGluon
+
+
+ // Check energy conservation
+ Double_t pxs = 0.;
+ Double_t pys = 0.;
+ Double_t pzs = 0.;
+ Double_t es = 14000.;
+
+ for (Int_t i = 0; i < numpart; i++)
+ {
+ kst = fPyjets->K[0][i];
+ if (kst != 1 && kst != 2) continue;
+ pxs += fPyjets->P[0][i];
+ pys += fPyjets->P[1][i];
+ pzs += fPyjets->P[2][i];
+ es -= fPyjets->P[3][i];
+ }
+ if (TMath::Abs(pxs) > 1.e-2 ||
+ TMath::Abs(pys) > 1.e-2 ||
+ TMath::Abs(pzs) > 1.e-1) {
+ printf("%e %e %e %e\n", pxs, pys, pzs, es);
+// Fatal("Quench()", "4-Momentum non-conservation");
+ }
+
+ } // end quenching loop (systems)
+// Clean-up
+ for (Int_t i = 0; i < numpart; i++)
+ {
+ imo = fPyjets->K[2][i];
+ if (imo > 1000) {
+ fPyjets->K[2][i] = fPyjets->K[2][i] % 1000;
+ }
+ }
+// this->Pylist(1);
+} // end quench
+
+
+void AliPythia6::Pyquen(Double_t a, Int_t ibf, Double_t b)
+{
+ // Igor Lokthine's quenching routine
+ // http://lokhtin.web.cern.ch/lokhtin/pyquen/pyquen.txt
+
+ pyquen(a, ibf, b);
+}
+
+void AliPythia6::SetPyquenParameters(Double_t t0, Double_t tau0, Int_t nf, Int_t iengl, Int_t iangl)
+{
+ // Set the parameters for the PYQUEN package.
+ // See comments in PyquenCommon.h
+
+
+ PYQPAR.t0 = t0;
+ PYQPAR.tau0 = tau0;
+ PYQPAR.nf = nf;
+ PYQPAR.iengl = iengl;
+ PYQPAR.iangl = iangl;
+}
+
+void AliPythia6::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
+{
+//
+// Load event into Pythia Common Block
+//
+
+ Int_t npart = stack -> GetNprimary();
+ Int_t n0 = 0;
+
+ if (!flag) {
+ GetPyjets()->N = npart;
+ } else {
+ n0 = GetPyjets()->N;
+ GetPyjets()->N = n0 + npart;
+ }
+
+
+ for (Int_t part = 0; part < npart; part++) {
+ TParticle *mPart = stack->Particle(part);
+
+ Int_t kf = mPart->GetPdgCode();
+ Int_t ks = mPart->GetStatusCode();
+ Int_t idf = mPart->GetFirstDaughter();
+ Int_t idl = mPart->GetLastDaughter();
+
+ if (reHadr) {
+ if (ks == 11 || ks == 12) {
+ ks -= 10;
+ idf = -1;
+ idl = -1;
+ }
+ }
+
+ Float_t px = mPart->Px();
+ Float_t py = mPart->Py();
+ Float_t pz = mPart->Pz();
+ Float_t e = mPart->Energy();
+ Float_t m = mPart->GetCalcMass();
+
+
+ (GetPyjets())->P[0][part+n0] = px;
+ (GetPyjets())->P[1][part+n0] = py;
+ (GetPyjets())->P[2][part+n0] = pz;
+ (GetPyjets())->P[3][part+n0] = e;
+ (GetPyjets())->P[4][part+n0] = m;
+
+ (GetPyjets())->K[1][part+n0] = kf;
+ (GetPyjets())->K[0][part+n0] = ks;
+ (GetPyjets())->K[3][part+n0] = idf + 1;
+ (GetPyjets())->K[4][part+n0] = idl + 1;
+ (GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
+ }
+}
+
+
+void AliPythia6::Pyevnw()
+{
+ // New multiple interaction scenario
+ pyevnw();
+}
+
+void AliPythia6::GetQuenchingParameters(Double_t& xp, Double_t& yp, Double_t z[4])
+{
+ // Return event specific quenching parameters
+ xp = fXJet;
+ yp = fYJet;
+ for (Int_t i = 0; i < 4; i++) z[i] = fZQuench[i];
+
+}
+
+void AliPythia6::ConfigHeavyFlavor()
+{
+ //
+ // Default configuration for Heavy Flavor production
+ //
+ // All QCD processes
+ //
+ SetMSEL(1);
+
+ // No multiple interactions
+ SetMSTP(81,0);
+ SetPARP(81, 0.);
+ SetPARP(82, 0.);
+ // Initial/final parton shower on (Pythia default)
+ SetMSTP(61,1);
+ SetMSTP(71,1);
+
+ // 2nd order alpha_s
+ SetMSTP(2,2);
+
+ // QCD scales
+ SetMSTP(32,2);
+ SetPARP(34,1.0);
+}
+
+void AliPythia6::AtlasTuning()
+{
+ //
+ // Configuration for the ATLAS tuning
+ SetMSTP(51,AliStructFuncType::PDFsetIndex(kCTEQ5L));
+ SetMSTP(81,1); // Multiple Interactions ON
+ SetMSTP(82,4); // Double Gaussian Model
+ SetPARP(81,1.9); // Min. pt for multiple interactions (default in 6.2-14)
+ SetPARP(82,1.8); // [GeV] PT_min at Ref. energy
+ SetPARP(89,1000.); // [GeV] Ref. energy
+ SetPARP(90,0.16); // 2*epsilon (exponent in power law)
+ SetPARP(83,0.5); // Core density in proton matter distribution (def.value)
+ SetPARP(84,0.5); // Core radius
+ SetPARP(85,0.33); // Regulates gluon prod. mechanism
+ SetPARP(86,0.66); // Regulates gluon prod. mechanism
+ SetPARP(67,1); // Regulates Initial State Radiation
+}
+
+void AliPythia6::SetPtHardRange(Float_t ptmin, Float_t ptmax)
+{
+ // Set the pt hard range
+ SetCKIN(3, ptmin);
+ SetCKIN(4, ptmax);
+}
+
+void AliPythia6::SetYHardRange(Float_t ymin, Float_t ymax)
+{
+ // Set the y hard range
+ SetCKIN(7, ymin);
+ SetCKIN(8, ymax);
+}
+
+
+void AliPythia6::SetFragmentation(Int_t flag)
+{
+ // Switch fragmentation on/off
+ SetMSTP(111, flag);
+}
+
+void AliPythia6::SetInitialAndFinalStateRadiation(Int_t flag1, Int_t flag2)
+{
+// initial state radiation
+ SetMSTP(61, flag1);
+// final state radiation
+ SetMSTP(71, flag2);
+}
+
+void AliPythia6::SetIntrinsicKt(Float_t kt)
+{
+ // Set the inreinsic kt
+ if (kt > 0.) {
+ SetMSTP(91,1);
+ SetPARP(91,kt);
+ SetPARP(93, 4. * kt);
+ } else {
+ SetMSTP(91,0);
+ }
+}
+
+void AliPythia6::SwitchHFOff()
+{
+ // Switch off heavy flavor
+ // Maximum number of quark flavours used in pdf
+ SetMSTP(58, 3);
+ // Maximum number of flavors that can be used in showers
+ SetMSTJ(45, 3);
+}
+
+void AliPythia6::SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
+ Float_t thresh, Float_t etseed, Float_t minet, Float_t r)
+{
+// Set pycell parameters
+ SetPARU(51, etamax);
+ SetMSTU(51, neta);
+ SetMSTU(52, nphi);
+ SetPARU(58, thresh);
+ SetPARU(52, etseed);
+ SetPARU(53, minet);
+ SetPARU(54, r);
+ SetMSTU(54, 2);
+}
+
+void AliPythia6::ModifiedSplitting()
+{
+ // Modified splitting probability as a model for quenching
+ SetPARJ(200, 0.8);
+ SetMSTJ(41, 1); // QCD radiation only
+ SetMSTJ(42, 2); // angular ordering
+ SetMSTJ(44, 2); // option to run alpha_s
+ SetMSTJ(47, 0); // No correction back to hard scattering element
+ SetMSTJ(50, 0); // No coherence in first branching
+ SetPARJ(82, 1.); // Cut off for parton showers
+}
+
+void AliPythia6::SwitchHadronisationOff()
+{
+ // Switch off hadronisarion
+ SetMSTJ(1, 0);
+}
+
+void AliPythia6::SwitchHadronisationOn()
+{
+ // Switch on hadronisarion
+ SetMSTJ(1, 1);
+}
+
+
+void AliPythia6::GetXandQ(Float_t& x1, Float_t& x2, Float_t& q)
+{
+ // Get x1, x2 and Q for this event
+
+ q = GetVINT(51);
+ x1 = GetVINT(41);
+ x2 = GetVINT(42);
+}
+
+Float_t AliPythia6::GetXSection()
+{
+ // Get the total cross-section
+ return (GetPARI(1));
+}
+
+Float_t AliPythia6::GetPtHard()
+{
+ // Get the pT hard for this event
+ return GetVINT(47);
+}
+
+Int_t AliPythia6::ProcessCode()
+{
+ // Get the subprocess code
+ return GetMSTI(1);
+}
+
+void AliPythia6::PrintStatistics()
+{
+ // End of run statistics
+ Pystat(1);
+}
+
+void AliPythia6::EventListing()
+{
+ // End of run statistics
+ Pylist(2);
+}
+
+AliPythia6& AliPythia6::operator=(const AliPythia6& rhs)
+{
+// Assignment operator
+ rhs.Copy(*this);
+ return *this;
+}
+
+ void AliPythia6::Copy(TObject&) const
+{
+ //
+ // Copy
+ //
+ Fatal("Copy","Not implemented!\n");
+}
--- /dev/null
+#ifndef ALIPYTHIA6_H
+#define ALIPYTHIA6_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+/* $Id: AliPythia.h,v 1.22 2007/10/09 08:43:24 morsch Exp $ */
+
+#include <TPythia6.h>
+#include "AliPythiaBase.h"
+
+
+class AliFastGlauber;
+class AliQuenchingWeights;
+class AliStack;
+
+class AliPythia6 : public TPythia6, public AliPythiaBase
+{
+
+ public:
+ virtual ~AliPythia6(){;}
+ virtual Int_t Version() {return (6);}
+ // convert to compressed code and print result (for debugging only)
+ virtual Int_t CheckedLuComp(Int_t kf);
+ // Pythia initialisation for selected processes
+ virtual void ProcInit
+ (Process_t process, Float_t energy, StrucFunc_t strucfunc);
+ virtual void GenerateEvent() {Pyevnt();}
+ virtual void GenerateMIEvent() {Pyevnw();}
+ virtual void HadronizeEvent() {Pyexec();}
+ virtual Int_t GetNumberOfParticles() {return GetN();}
+ virtual void SetNumberOfParticles(Int_t i) {SetN(i);}
+ virtual void EditEventList(Int_t i) {Pyedit(i);}
+ virtual void PrintStatistics();
+ virtual void EventListing();
+ virtual Int_t GetParticles(TClonesArray *particles) {return ImportParticles(particles, "All");}
+ // Treat protons as inside nuclei
+ virtual void SetNuclei(Int_t a1, Int_t a2);
+ // Print particle properties
+ virtual void PrintParticles();
+ // Reset the decay table
+ virtual void ResetDecayTable();
+ //
+ // Common Physics Configuration
+ virtual void SetPtHardRange(Float_t ptmin, Float_t ptmax);
+ virtual void SetYHardRange(Float_t ymin, Float_t ymax);
+ virtual void SetFragmentation(Int_t flag);
+ virtual void SetInitialAndFinalStateRadiation(Int_t flag1, Int_t flag2);
+ virtual void SetIntrinsicKt(Float_t kt);
+ virtual void SwitchHFOff();
+ virtual void SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
+ Float_t thresh, Float_t etseed, Float_t minet, Float_t r);
+ virtual void ModifiedSplitting();
+ virtual void SwitchHadronisationOff();
+ virtual void SwitchHadronisationOn();
+ //
+ // Common Getters
+ virtual void GetXandQ(Float_t& x1, Float_t& x2, Float_t& q);
+ virtual Float_t GetXSection();
+ virtual Int_t ProcessCode();
+ virtual Float_t GetPtHard();
+ //
+ //
+ virtual void SetDecayTable();
+ virtual void Pyevnw();
+ virtual void Pycell(Int_t& nclus);
+ virtual void Pyclus(Int_t& nclus);
+ virtual void GetJet(Int_t i, Float_t& px, Float_t& py, Float_t& pz, Float_t& e);
+ virtual void LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr);
+ virtual void Pyshow(Int_t ip1, Int_t ip2, Double_t qmax);
+ virtual void Pyrobo(Int_t imi, Int_t ima, Double_t the, Double_t phi, Double_t bex, Double_t bey, Double_t bez);
+ virtual void InitQuenching(Float_t bmin, Float_t bmax, Float_t k, Int_t iECMethod, Float_t zmax = 0.97, Int_t ngmax = 30);
+ virtual void SetPyquenParameters(Double_t t0, Double_t tau0, Int_t nf, Int_t iengl, Int_t iangl);
+ virtual void Pyquen(Double_t a, Int_t ibf, Double_t b);
+ virtual void GetQuenchingParameters(Double_t& xp, Double_t& yp, Double_t z[4]);
+ // return instance of the singleton
+ static AliPythia6* Instance();
+ virtual void Quench();
+ // Assignment Operator
+ AliPythia6 & operator=(const AliPythia6 & rhs);
+ void Copy(TObject&) const;
+ protected:
+ Process_t fProcess; // Process type
+ Float_t fEcms; // Centre of mass energy
+ StrucFunc_t fStrucFunc; // Structure function
+ Int_t fDefMDCY[501]; // ! Default decay switches per particle
+ Int_t fDefMDME[2001]; // ! Default decay switches per mode
+ Double_t fZQuench[4]; // ! Quenching fractions for this even
+ Double_t fXJet; // ! Jet production point X
+ Double_t fYJet; // ! Jet production point Y
+ Int_t fNGmax; // Maximum number of radiated gluons in quenching
+ Float_t fZmax; // Maximum energy loss in quenching
+ AliFastGlauber* fGlauber; // ! The Glauber model
+ AliQuenchingWeights* fQuenchingWeights; // ! The Quenching Weights model
+ static AliPythia6* fgAliPythia; // Pointer to single instance
+ private:
+ AliPythia6();
+ AliPythia6(const AliPythia6& pythia);
+ void ConfigHeavyFlavor();
+ void AtlasTuning();
+ ClassDef(AliPythia6,1) //ALICE UI to PYTHIA
+};
+
+#endif
+
+
+
+
+
#pragma link off all functions;
#pragma link C++ enum Process_t;
+#pragma link C++ class AliPythiaBase+;
#pragma link C++ class AliPythia+;
+#pragma link C++ class AliPythia6+;
#pragma link C++ class AliGenPythia+;
+#pragma link C++ class AliGenPythiaPlus+;
#pragma link C++ class AliDecayerPythia+;
-#pragma link C++ class AliGenPythiaJets+;
#endif
--- /dev/null
+
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+/* $Id$ */
+
+#include "AliPythiaBase.h"
+
+ClassImp(AliPythiaBase)
+
+AliPythiaBase::AliPythiaBase()
+{
+ // Default constructor
+}
+
+
+
--- /dev/null
+#ifndef ALIPYTHIABASE_H
+#define ALIPYTHIABASE_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+/* $Id$ */
+#include "AliRndm.h"
+#include <TObject.h>
+#include "AliStructFuncType.h"
+#include "PythiaProcesses.h"
+
+class AliFastGlauber;
+class AliQuenchingWeights;
+class AliStack;
+class TClonesArray;
+
+class AliPythiaBase : public AliRndm
+{
+
+ public:
+ AliPythiaBase();
+ virtual ~AliPythiaBase(){;}
+ virtual Int_t Version() = 0;
+ // Convert to compressed code and print result (for debugging only)
+ virtual Int_t CheckedLuComp(Int_t kf) = 0;
+ // Pythia initialisation for selected processes
+ virtual void ProcInit (Process_t process, Float_t energy, StrucFunc_t strucfunc) = 0;
+ virtual void GenerateEvent() = 0;
+ virtual void GenerateMIEvent() = 0;
+ virtual Int_t GetNumberOfParticles() = 0;
+ virtual void SetNumberOfParticles(Int_t i) = 0;
+ virtual void EditEventList(Int_t i) = 0;
+ virtual void HadronizeEvent() = 0;
+ virtual Int_t GetParticles(TClonesArray *particles) = 0;
+ virtual void PrintStatistics() = 0;
+ virtual void EventListing() = 0;
+ // Treat protons as inside nuclei
+ virtual void SetNuclei(Int_t a1, Int_t a2) = 0;
+ // Print particle properties
+ virtual void PrintParticles() = 0;
+ // Reset the decay table
+ virtual void ResetDecayTable() = 0;
+ //
+ // Common Physics Configuration
+ virtual void SetPtHardRange(Float_t ptmin, Float_t ptmax) = 0;
+ virtual void SetYHardRange(Float_t ymin, Float_t ymax) = 0;
+ virtual void SetFragmentation(Int_t flag) = 0;
+ virtual void SetInitialAndFinalStateRadiation(Int_t flag1, Int_t flag2) = 0;
+ virtual void SetIntrinsicKt(Float_t kt) = 0;
+ virtual void SwitchHFOff() = 0;
+ virtual void SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
+ Float_t thresh, Float_t etseed, Float_t minet, Float_t r) = 0;
+ virtual void ModifiedSplitting() = 0;
+ virtual void SwitchHadronisationOff() = 0;
+ virtual void SwitchHadronisationOn() = 0;
+ //
+ // Common Getters
+ virtual void GetXandQ(Float_t& x1, Float_t& x2, Float_t& q) = 0;
+ virtual Float_t GetXSection() = 0;
+ virtual Int_t ProcessCode() = 0;
+ virtual Float_t GetPtHard() = 0;
+ //
+ //
+ virtual void SetDecayTable() = 0;
+ virtual void Pycell(Int_t& nclus) = 0;
+ virtual void Pyclus(Int_t& nclus) = 0;
+ virtual void GetJet(Int_t i, Float_t& px, Float_t& py, Float_t& pz, Float_t& e) = 0;
+ virtual void LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr) = 0;
+ virtual void Pyshow(Int_t ip1, Int_t ip2, Double_t qmax) = 0;
+ virtual void Pyrobo(Int_t imi, Int_t ima, Double_t the, Double_t phi, Double_t bex, Double_t bey, Double_t bez) = 0;
+ virtual void InitQuenching(Float_t bmin, Float_t bmax, Float_t k, Int_t iECMethod, Float_t zmax, Int_t ngmax) = 0;
+ virtual void Pyquen(Double_t a, Int_t ibf, Double_t b) = 0;
+ virtual void GetQuenchingParameters(Double_t& xp, Double_t& yp, Double_t z[4]) = 0;
+ // return instance of the singleton
+ virtual void Quench() = 0;
+ virtual void ConfigHeavyFlavor() = 0;
+ virtual void AtlasTuning() = 0;
+ ClassDef(AliPythiaBase, 1) //ALICE UI to PYTHIA
+};
+
+#endif
+
+
+
+
+
-#-*- Mode: Makefile -*-
-
-SRCS= AliPythia.cxx AliGenPythia.cxx AliDecayerPythia.cxx \
-AliGenPythiaJets.cxx
+SRCS= AliPythiaBase.cxx AliPythia.cxx AliPythia6.cxx AliGenPythia.cxx AliGenPythiaPlus.cxx AliDecayerPythia.cxx
HDRS= $(SRCS:.cxx=.h)
DHDR:=AliPythia6LinkDef.h
-EXPORT:=AliPythia.h AliDecayerPythia.h
+EXPORT:=AliPythiaBase.h AliPythia.h AliDecayerPythia.h
EINCLUDE:=FASTSIM
-
-ifeq (win32gcc,$(ALICE_TARGET))
-PACKSOFLAGS:= $(SOFLAGS) -L$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET) \
- -lEVGEN -lSTEER -lpythia6 -lSTEERBase -lFASTSIM \
- -L$(shell root-config --libdir) -lVMC -lEG -lEGPythia6
-endif
git://git.uio.no / u/mrichter/AliRoot.git / commitdiff