--- /dev/null
+//\r
+// Configuration for the Physics Data Challenge 2007\r
+//\r
+//\r
+// Assuming inel = 70 mb (PPR v.1, p.64)\r
+//\r
+// 84.98% of MSEL=0 events (including diffractive) with \r
+// QQbar switched off (these events will include injected J/psi) => kPyMbNoHvq\r
+// \r
+// 14.14% of MSEL=1 events with ccbar (in 4 subsamples) => kCharmpp14000wmi \r
+// bin 1 25% (3.535%): 2.76 < pthard < 3 GeV/c \r
+// bin 2 40% (5.656%): 3 < pthard < 4 GeV/c \r
+// bin 3 29% (4.101%): 4 < pthard < 8 GeV/c \r
+// bin 4 6% (0.848%): pthard > 8 GeV/c \r
+//\r
+// 0.73% of MSEL=1 events with bbbar (in 4 subsamples) => kBeautypp14000wmi\r
+// bin 1 5% (0.037%): 2.76 < pthard < 4 GeV/c\r
+// bin 2 31% (0.226%): 4 < pthard < 6 GeV/c\r
+// bin 3 28% (0.204%): 6 < pthard < 8 GeV/c\r
+// bin 4 36% (0.263%): pthard >8 GeV/c\r
+//\r
+// 0.075% of MSEL=0 events with QQbar switched off and 1 Omega- => kPyOmegaMinus\r
+// 0.075% of MSEL=0 events with QQbar switched off and 1 OmegaBar+ => kPyOmegaPlus\r
+\r
+// One can use the configuration macro in compiled mode by\r
+// root [0] gSystem->Load("libgeant321");\r
+// root [0] gSystem->SetIncludePath("-I$ROOTSYS/include -I$ALICE_ROOT/include\\r
+// -I$ALICE_ROOT -I$ALICE/geant3/TGeant3");\r
+// root [0] .x grun.C(1,"Config_PDC07.C++")\r
+\r
+class AliGenPythia;\r
+\r
+#if !defined(__CINT__) || defined(__MAKECINT__)\r
+#include <Riostream.h>\r
+#include <TRandom.h>\r
+#include <TDatime.h>\r
+#include <TSystem.h>\r
+#include <TVirtualMC.h>\r
+#include <TGeant3TGeo.h>\r
+#include "EVGEN/AliGenCocktail.h"\r
+#include "EVGEN/AliGenParam.h"\r
+#include "EVGEN/AliGenMUONlib.h"\r
+#include "STEER/AliRunLoader.h"\r
+#include "STEER/AliRun.h"\r
+#include "STEER/AliConfig.h"\r
+#include "PYTHIA6/AliDecayerPythia.h"\r
+#include "PYTHIA6/AliGenPythia.h"\r
+#include "STEER/AliMagFMaps.h"\r
+#include "STRUCT/AliBODY.h"\r
+#include "STRUCT/AliMAG.h"\r
+#include "STRUCT/AliABSOv3.h"\r
+#include "STRUCT/AliDIPOv3.h"\r
+#include "STRUCT/AliHALLv3.h"\r
+#include "STRUCT/AliFRAMEv2.h"\r
+#include "STRUCT/AliSHILv3.h"\r
+#include "STRUCT/AliPIPEv3.h"\r
+#include "ITS/AliITSgeom.h"\r
+#include "ITS/AliITSv11Hybrid.h"\r
+#include "TPC/AliTPCv2.h"\r
+#include "TOF/AliTOFv6T0.h"\r
+#include "HMPID/AliHMPIDv3.h"\r
+#include "HMPID/AliHMPIDv2.h"\r
+#include "ZDC/AliZDCv3.h"\r
+#include "TRD/AliTRDv1.h"\r
+#include "FMD/AliFMDv1.h"\r
+#include "MUON/AliMUONv1.h"\r
+#include "PHOS/AliPHOSv1.h"\r
+#include "PMD/AliPMDv1.h"\r
+#include "T0/AliT0v1.h"\r
+#include "EMCAL/AliEMCALv2.h"\r
+#include "ACORDE/AliACORDEv1.h"\r
+#include "VZERO/AliVZEROv7.h"\r
+#endif\r
+\r
+\r
+ // libraries required by geant321\r
+\r
+\r
+\r
+enum PDC07Proc_t \r
+{\r
+//--- Heavy Flavour Production ---\r
+ kCharmPbPb5500, kCharmpPb8800, kCharmpp14000, kCharmpp14000wmi,\r
+ kD0PbPb5500, kD0pPb8800, kD0pp14000,\r
+ kDPlusPbPb5500, kDPluspPb8800, kDPluspp14000,\r
+ kBeautyPbPb5500, kBeautypPb8800, kBeautypp14000, kBeautypp14000wmi, \r
+// -- Pythia Mb\r
+ kPyMbNoHvq, kPyOmegaPlus, kPyOmegaMinus, kPyJetJet, kPyJetJetFixed,\r
+ kPyGammaJetPHOS, kPyJetJetPHOS, kPyJetJetPHOSv2, kPyGammaBremsPHOS, \r
+ kPyGammaJetEMCAL, kPyJetJetEMCAL, kPyGammaBremsEMCAL, kRunMax\r
+};\r
+\r
+const char * pprRunName[] = {\r
+ "kCharmPbPb5500", "kCharmpPb8800", "kCharmpp14000", "kCharmpp14000wmi",\r
+ "kD0PbPb5500", "kD0pPb8800", "kD0pp14000",\r
+ "kDPlusPbPb5500", "kDPluspPb8800", "kDPluspp14000",\r
+ "kBeautyPbPb5500", "kBeautypPb8800", "kBeautypp14000", "kBeautypp14000wmi", \r
+ "kPyMbNoHvq", "kPyOmegaPlus", "kPyOmegaMinus", "kPyJetJet","kPyJetJetFixed",\r
+ "kPyGammaJetPHOS", "kPyJetJetPHOS", "kPyJetJetPHOSv2", "kPyGammaBremsPHOS",\r
+ "kPyGammaJetEMCAL", "kPyJetJetEMCAL", "kPyGammaBremsEMCAL"\r
+};\r
+\r
+\r
+//--- Decay Mode ---\r
+enum DecayHvFl_t \r
+{\r
+ kNature, kHadr, kSemiEl, kSemiMu\r
+};\r
+//--- Magnetic Field ---\r
+enum Mag_t\r
+{\r
+ kNoField, k5kG\r
+};\r
+\r
+//--- Trigger config ---\r
+enum TrigConf_t\r
+{\r
+ kDefaultPPTrig, kDefaultPbPbTrig\r
+};\r
+\r
+const char * TrigConfName[] = {\r
+ "p-p","Pb-Pb"\r
+};\r
+\r
+Float_t eCMS=14000;\r
+PDC07Proc_t proc = kPyJetJetEMCAL;\r
+Int_t year = 2009;\r
+\r
+\r
+enum PprGeo_t\r
+ {\r
+ kHoles, kNoHoles\r
+ };\r
+\r
+static PprGeo_t geo = kHoles;\r
+\r
+//--- Functions ---\r
+AliGenPythia *PythiaHVQ(PDC07Proc_t proc);\r
+\r
+AliGenPythia* Blubb(Int_t proc) {\r
+ //*******************************************************************//\r
+ // Configuration file for hard QCD processes generation with PYTHIA //\r
+ // //\r
+ //*******************************************************************//\r
+ AliGenPythia * gener = 0x0;\r
+\r
+ switch(proc) {\r
+ \r
+ case kPyJetJet:\r
+ comment = comment.Append(Form("pp->jet + jet at %3.0f TeV, pt hard %3.0f - %3.0f",eCMS,ptHardMin,ptHardMax));\r
+ AliGenPythia * gener = new AliGenPythia(nEvts);\r
+ gener->SetEnergyCMS(eCMS);// Centre of mass energy\r
+ gener->SetProcess(kPyJets);// Process type\r
+ gener->SetJetEtaRange(-1.5, 1.5);// Final state kinematic cuts\r
+ gener->SetJetPhiRange(0., 360.);\r
+ gener->SetJetEtRange(10., 1000.);\r
+ gener->SetPtHard(ptHardMin, ptHardMax);// Pt transfer of the hard scattering\r
+ gener->SetStrucFunc(kCTEQ5L);\r
+ \r
+ return gener;\r
+ \r
+ case kPyGammaJetPHOS:\r
+ comment = comment.Append(" pp->jet + gamma over PHOS");\r
+ gener = new AliGenPythia(nEvts);\r
+ gener->SetEnergyCMS(eCMS);\r
+ gener->SetProcess(kPyDirectGamma);\r
+ gener->SetStrucFunc(kCTEQ4L);\r
+ gener->SetPtHard(ptHardMin,ptHardMax);\r
+ //gener->SetYHard(-1.0,1.0);\r
+ gener->SetGammaEtaRange(-0.13,0.13);\r
+ gener->SetGammaPhiRange(218.,322.);//Over 5 modules +-2 deg\r
+ break;\r
+ case kPyJetJetPHOS:\r
+ comment = comment.Append(" pp->jet + jet over PHOS");\r
+ gener = new AliGenPythia(nEvts);\r
+ gener->SetEnergyCMS(eCMS);\r
+ gener->SetProcess(kPyJets);\r
+ gener->SetStrucFunc(kCTEQ4L);\r
+ gener->SetPtHard(ptHardMin,ptHardMax);\r
+ //gener->SetYHard(-1.0,1.0);\r
+ gener->SetJetEtaRange(-1.,1.);\r
+ gener->SetJetPhiRange(200.,340.); \r
+ gener->SetPi0InPHOS(kTRUE);\r
+ gener->SetFragPhotonOrPi0MinPt(ptGammaPi0Min);\r
+ \r
+ printf("\n \n Event generator: Minimum pT of particle in calorimeter %f \n \n", ptGammaPi0Min);\r
+ break;\r
+ case kPyGammaBremsPHOS:\r
+ comment = comment.Append(" pp->jet + jet+bremsphoton over PHOS at 14 TeV");\r
+ gener = new AliGenPythia(nEvts);\r
+ gener->SetEnergyCMS(eCMS);\r
+ gener->SetProcess(kPyJets);\r
+ gener->SetStrucFunc(kCTEQ4L);\r
+ gener->SetPtHard(ptHardMin,ptHardMax);\r
+ //gener->SetYHard(-1.0,1.0);\r
+ gener->SetJetEtaRange(-1.,1.);\r
+ gener->SetJetPhiRange(200.,340.); \r
+ gener->SetFragPhotonInPHOS(kTRUE);\r
+ gener->SetFragPhotonOrPi0MinPt(ptGammaPi0Min);\r
+ printf("\n \n Event generator: Minimum pT of particle in calorimeter %f \n \n", ptGammaPi0Min);\r
+ break;\r
+ case kPyJetJetPHOSv2:\r
+ comment = comment.Append(" pp->jet + jet over PHOS version2 ");\r
+ gener = new AliGenPythia(nEvts);\r
+ gener->SetEnergyCMS(eCMS);\r
+ gener->SetProcess(kPyJets);\r
+ gener->SetStrucFunc(kCTEQ4L);\r
+ gener->SetPtHard(ptHardMin,ptHardMax);\r
+ //gener->SetYHard(-1.0,1.0);\r
+ gener->SetJetEtaRange(-1.,1.);\r
+ gener->SetJetPhiRange(200.,340.); \r
+ //gener->SetPi0InPHOS(kTRUE);\r
+ gener->SetPhotonInPHOSeta(kTRUE);\r
+ gener->SetPhotonMinPt(ptGammaPi0Min);\r
+ gener->SetForceDecay(kAll);\r
+ break;\r
+ case kPyGammaJetEMCAL:\r
+ comment = comment.Append(" pp->jet + gamma over EMCAL at 14 TeV");\r
+ gener = new AliGenPythia(nEvts);\r
+ gener->SetEnergyCMS(eCMS);\r
+ gener->SetProcess(kPyDirectGamma);\r
+ gener->SetStrucFunc(kCTEQ4L);\r
+ gener->SetPtHard(ptHardMin,ptHardMax);\r
+ //gener->SetYHard(-1.0,1.0);\r
+ gener->SetGammaEtaRange(-0.71,0.71);\r
+ gener->SetGammaPhiRange(78.,192.);//Over 6 supermodules +-2 deg \r
+ break;\r
+ case kPyJetJetEMCAL:\r
+ comment = comment.Append(" pp->jet + jet over EMCAL at 14 TeV");\r
+ gener = new AliGenPythia(nEvts);\r
+ gener->SetEnergyCMS(eCMS);\r
+ gener->SetProcess(kPyJets);\r
+ gener->SetStrucFunc(kCTEQ4L);\r
+ gener->SetPtHard(ptHardMin,ptHardMax);\r
+ //gener->SetYHard(-1.0,1.0);\r
+ gener->SetJetEtaRange(-1,1);\r
+ gener->SetJetPhiRange(60.,210.); \r
+ gener->SetPi0InEMCAL(kTRUE);\r
+ gener->SetFragPhotonOrPi0MinPt(ptGammaPi0Min);\r
+ printf("\n \n Event generator: Minimum pT of particle in calorimeter %f \n \n", ptGammaPi0Min);\r
+ break;\r
+ case kPyGammaBremsEMCAL:\r
+ comment = comment.Append(" pp->jet + jet+bremsphoton over EMCAL at 14 TeV");\r
+ gener = new AliGenPythia(nEvts);\r
+ gener->SetEnergyCMS(eCMS);\r
+ gener->SetProcess(kPyJets);\r
+ gener->SetStrucFunc(kCTEQ4L);\r
+ gener->SetPtHard(ptHardMin,ptHardMax);\r
+ //gener->SetYHard(-1.0,1.0);\r
+ gener->SetJetEtaRange(-1,1);\r
+ gener->SetJetPhiRange(60.,210.); //Over 2 uncovered PHOS modules\r
+ gener->SetFragPhotonInEMCAL(kTRUE);\r
+ gener->SetFragPhotonOrPi0MinPt(ptGammaPi0Min);\r
+ \r
+ printf("\n \n Event generator: Minimum pT of particle in calorimeter %f \n \n", ptGammaPi0Min);\r
+ break;\r
+ case kPyJetJetFixed:\r
+ comment = comment.Append(" pp->jet + jet+jet in Fixed bin for testing");\r
+ gener = new AliGenPythia(nEvts);\r
+ gener->SetEnergyCMS(14000);// Centre of mass energy\r
+ gener->SetProcess(kPyJets);// Process type\r
+ gener->SetJetEtaRange(-1.1, 1.1);// Final state kinematic cuts\r
+ gener->SetJetPhiRange(0., 360.);\r
+ gener->SetJetEtRange(95., 105.);\r
+ gener->SetPtHard(80, 1000);// Pt transfer of the hard scattering\r
+ gener->SetStrucFunc(kCTEQ5L);\r
+ gener->SetGluonRadiation(1,1); \r
+ \r
+\r
+ printf("\n \n Fixed Parameter \n");\r
+ break;\r
+\r
+ }\r
+\r
+ return gener;\r
+}\r
+\r
+AliGenerator *MbCocktail();\r
+AliGenerator *PyMbTriggered(Int_t pdg);\r
+void ProcessEnvironmentVars();\r
+\r
+// This part for configuration\r
+static DecayHvFl_t decHvFl = kNature; \r
+static Mag_t mag = k5kG; \r
+static TrigConf_t trig = kDefaultPPTrig; // default pp trigger configuration\r
+static Int_t runNumber= 0; \r
+static Int_t eventNumber= 0; \r
+//========================//\r
+// Set Random Number seed //\r
+//========================//\r
+TDatime dt;\r
+static UInt_t seed = dt.Get();\r
+\r
+// nEvts = -1 : you get 1 QQbar pair and all the fragmentation and \r
+// decay chain\r
+// nEvts = N>0 : you get N charm / beauty Hadrons \r
+Int_t nEvts = -1; \r
+// stars = kTRUE : all heavy resonances and their decay stored\r
+// = kFALSE: only final heavy hadrons and their decays stored\r
+Bool_t stars = kTRUE;\r
+\r
+// To be used only with kCharmpp14000wmi and kBeautypp14000wmi\r
+// To get a "reasonable" agreement with MNR results, events have to be \r
+// generated with the minimum ptHard set to 2.76 GeV.\r
+// To get a "perfect" agreement with MNR results, events have to be \r
+// generated in four ptHard bins with the following relative \r
+// normalizations:\r
+// CHARM\r
+// 2.76-3 GeV: 25%\r
+// 3-4 GeV: 40%\r
+// 4-8 GeV: 29%\r
+// >8 GeV: 6%\r
+// BEAUTY\r
+// 2.76-4 GeV: 5% \r
+// 4-6 GeV: 31%\r
+// 6-8 GeV: 28%\r
+// >8 GeV: 36%\r
+\r
+Float_t ptHardMin = 10.;\r
+Float_t ptHardMax = 20.;\r
+Float_t ptGammaPi0Min = 1.;\r
+Int_t iquenching = 0;\r
+Float_t qhat = 20.;\r
+\r
+// Comment line\r
+static TString comment;\r
+\r
+void Config()\r
+{\r
+ \r
+\r
+ // Get settings from environment variables\r
+ ProcessEnvironmentVars();\r
+\r
+#if defined(__CINT__)\r
+ gSystem->Load("liblhapdf");\r
+ gSystem->Load("libEGPythia6");\r
+ gSystem->Load("libpythia6");\r
+ gSystem->Load("libAliPythia6");\r
+ gSystem->Load("libgeant321");\r
+#endif\r
+\r
+ new TGeant3TGeo("C++ Interface to Geant3");\r
+\r
+ // Output every 100 tracks\r
+ ((TGeant3*)gMC)->SetSWIT(4,100);\r
+\r
+ //=======================================================================\r
+ // Run loader \r
+ AliRunLoader* rl=0x0;\r
+ rl = AliRunLoader::Open("galice.root",\r
+ AliConfig::GetDefaultEventFolderName(),\r
+ "recreate");\r
+ if (rl == 0x0)\r
+ {\r
+ gAlice->Fatal("Config.C","Can not instatiate the Run Loader");\r
+ return;\r
+ }\r
+ rl->SetCompressionLevel(2);\r
+ rl->SetNumberOfEventsPerFile(1000);\r
+ gAlice->SetRunLoader(rl);\r
+\r
+ // Run number\r
+ //gAlice->SetRunNumber(runNumber);\r
+ \r
+ // Set the trigger configuration\r
+ gAlice->SetTriggerDescriptor(TrigConfName[trig]);\r
+ cout<<"Trigger configuration is set to "<<TrigConfName[trig]<<endl;\r
+\r
+ //\r
+ //=======================================================================\r
+ // ************* STEERING parameters FOR ALICE SIMULATION **************\r
+ // --- Specify event type to be tracked through the ALICE setup\r
+ // --- All positions are in cm, angles in degrees, and P and E in GeV\r
+\r
+\r
+ gMC->SetProcess("DCAY",1);\r
+ gMC->SetProcess("PAIR",1);\r
+ gMC->SetProcess("COMP",1);\r
+ gMC->SetProcess("PHOT",1);\r
+ gMC->SetProcess("PFIS",0);\r
+ gMC->SetProcess("DRAY",0);\r
+ gMC->SetProcess("ANNI",1);\r
+ gMC->SetProcess("BREM",1);\r
+ gMC->SetProcess("MUNU",1);\r
+ gMC->SetProcess("CKOV",1);\r
+ gMC->SetProcess("HADR",1);\r
+ gMC->SetProcess("LOSS",2);\r
+ gMC->SetProcess("MULS",1);\r
+ gMC->SetProcess("RAYL",1);\r
+\r
+ Float_t cut = 1.e-3; // 1MeV cut by default\r
+ Float_t tofmax = 1.e10;\r
+\r
+ gMC->SetCut("CUTGAM", cut);\r
+ gMC->SetCut("CUTELE", cut);\r
+ gMC->SetCut("CUTNEU", cut);\r
+ gMC->SetCut("CUTHAD", cut);\r
+ gMC->SetCut("CUTMUO", cut);\r
+ gMC->SetCut("BCUTE", cut); \r
+ gMC->SetCut("BCUTM", cut); \r
+ gMC->SetCut("DCUTE", cut); \r
+ gMC->SetCut("DCUTM", cut); \r
+ gMC->SetCut("PPCUTM", cut);\r
+ gMC->SetCut("TOFMAX", tofmax); \r
+\r
+ //======================//\r
+ // Set External decayer //\r
+ //======================//\r
+ TVirtualMCDecayer* decayer = new AliDecayerPythia();\r
+ // DECAYS\r
+ //\r
+ switch(decHvFl) {\r
+ case kNature:\r
+ decayer->SetForceDecay(kAll);\r
+ break;\r
+ case kHadr:\r
+ decayer->SetForceDecay(kHadronicD);\r
+ break;\r
+ case kSemiEl:\r
+ decayer->SetForceDecay(kSemiElectronic);\r
+ break;\r
+ case kSemiMu:\r
+ decayer->SetForceDecay(kSemiMuonic);\r
+ break;\r
+ }\r
+ decayer->Init();\r
+ gMC->SetExternalDecayer(decayer);\r
+ if(proc == kPyJetJetPHOSv2) // in this case we decay the pi0\r
+ decayer->SetForceDecay(kNeutralPion);\r
+\r
+ //=========================//\r
+ // Generator Configuration //\r
+ //=========================//\r
+ AliGenPythia* gener = 0x0;\r
+ \r
+ if (proc <= kBeautypp14000wmi) {\r
+ AliGenPythia *pythia = PythiaHVQ(proc);\r
+ // FeedDown option\r
+ pythia->SetFeedDownHigherFamily(kFALSE);\r
+ // Stack filling option\r
+ if(!stars) pythia->SetStackFillOpt(AliGenPythia::kParentSelection);\r
+ // Set Count mode\r
+ if(nEvts>0) pythia->SetCountMode(AliGenPythia::kCountParents);\r
+ //\r
+ // DECAYS\r
+ // \r
+ switch(decHvFl) {\r
+ case kNature:\r
+ pythia->SetForceDecay(kAll);\r
+ break;\r
+ case kHadr:\r
+ pythia->SetForceDecay(kHadronicD);\r
+ break;\r
+ case kSemiEl:\r
+ pythia->SetForceDecay(kSemiElectronic);\r
+ break;\r
+ case kSemiMu:\r
+ pythia->SetForceDecay(kSemiMuonic);\r
+ break;\r
+ }\r
+ //\r
+ // GEOM & KINE CUTS\r
+ //\r
+ pythia->SetMomentumRange(0,99999999);\r
+ pythia->SetPhiRange(0., 360.);\r
+ pythia->SetThetaRange(0,180);\r
+ switch(ycut) {\r
+ case kFull:\r
+ pythia->SetYRange(-999,999);\r
+ break;\r
+ case kBarrel:\r
+ pythia->SetYRange(-2,2);\r
+ break;\r
+ case kMuonArm:\r
+ pythia->SetYRange(1,6);\r
+ break;\r
+ }\r
+ gener = pythia;\r
+ } else if (proc == kPyMbNoHvq) {\r
+ gener = MbCocktail();\r
+ } else if (proc == kPyOmegaMinus) {\r
+ gener = PyMbTriggered(3334);\r
+ } else if (proc == kPyOmegaPlus) {\r
+ gener = PyMbTriggered(-3334);\r
+ } else if (proc <= kPyGammaBremsEMCAL) {\r
+ \r
+ AliGenPythia *pythia = Blubb((Int_t) proc ); \r
+ \r
+ // FeedDown option\r
+ pythia->SetFeedDownHigherFamily(kFALSE);\r
+ // Set Count mode\r
+ if(nEvts>0) pythia->SetCountMode(AliGenPythia::kCountParents);\r
+ \r
+ //\r
+ // GEOM & KINE CUTS\r
+ //\r
+ pythia->SetMomentumRange(0,99999999);\r
+ // pythia->SetJetEtaRange(-1.5, 1.5);// Final state kinematic cuts\r
+ // pythia->SetJetPhiRange(0., 360.);\r
+ // pythia->SetThetaRange(45,135);\r
+ \r
+ if(proc == kPyJetJetPHOSv2)\r
+ pythia->SetForceDecay(kNeutralPion);\r
+ else \r
+ pythia->SetForceDecay(kAll);\r
+ \r
+ pythia->SetPycellParameters(2., 274, 432, 0., 4., 5., 1.0);\r
+ pythia->SetPtKick(3.); // set the intrinsic kt to 5 GeV/c\r
+ pythia->SetGluonRadiation(1,1);\r
+ gener = pythia;\r
+ }\r
+ \r
+\r
+ // PRIMARY VERTEX\r
+\r
+ gener->SetOrigin(0., 0., 0.); // vertex position\r
+\r
+ // Size of the interaction diamond\r
+ // Longitudinal\r
+ Float_t sigmaz = 7.55 / TMath::Sqrt(2.); // [cm]\r
+\r
+ // Transverse\r
+ Float_t betast = 10; // beta* [m]\r
+ Float_t eps = 3.75e-6; // emittance [m]\r
+ Float_t gamma = 7000. / 0.938272; // relativistic gamma [1]\r
+ Float_t sigmaxy = TMath::Sqrt(eps * betast / gamma) / TMath::Sqrt(2.) * 100.; // [cm]\r
+ printf("\n \n Diamond size x-y: %10.3e z: %10.3e\n \n", sigmaxy, sigmaz);\r
+ \r
+ gener->SetSigma(sigmaxy, sigmaxy, sigmaz); // Sigma in (X,Y,Z) (cm) on IP position\r
+ gener->SetCutVertexZ(3.); // Truncate at 3 sigma\r
+ gener->SetVertexSmear(kPerEvent);\r
+\r
+ gener->Init();\r
+\r
+ //Quenching\r
+ gener->SetQuench(iquenching);\r
+ if(iquenching == 1){\r
+ Float_t k = 6e5*(qhat/1.7) ; //qhat=1.7, k = 6e5, default value \r
+ AliPythia::Instance()->InitQuenching(0.,0.1,k,0,0.95,6);\r
+\r
+ }\r
+ // FIELD\r
+ \r
+ printf("\n \n Comment: %s \n \n", comment.Data());\r
+ if (mag == kNoField) {\r
+ comment = comment.Append(" | L3 field 0.0 T");\r
+ field = new AliMagWrapCheb("Maps","Maps", 2, 0., 10., AliMagWrapCheb::k2kG);\r
+ } else if (mag == k5kG) {\r
+ comment = comment.Append(" | L3 field 0.5 T");\r
+ // field = new AliMagWrapCheb("Maps","Maps", 2, 1., 10., AliMagWrapCheb::k5kG);\r
+ field = new AliMagWrapCheb("Maps","Maps", 2, 1., 10., AliMagWrapCheb::k5kG,kTRUE,"$(ALICE_ROOT)/data/maps/mfchebKGI_sym.root");\r
+ }\r
+\r
+ printf("\n \n Comment: %s \n \n", comment.Data());\r
+ rl->CdGAFile();\r
+ gAlice->SetField(field);\r
+\r
+\r
+\r
+ Int_t iABSO = 1;\r
+ Int_t iACORDE = 0;\r
+ Int_t iDIPO = 1;\r
+ Int_t iEMCAL = 1;\r
+ Int_t iFMD = 1;\r
+ Int_t iFRAME = 1;\r
+ Int_t iHALL = 1;\r
+ Int_t iITS = 1;\r
+ Int_t iMAG = 1;\r
+ Int_t iMUON = 1;\r
+ Int_t iPHOS = 1;\r
+ Int_t iPIPE = 1;\r
+ Int_t iPMD = 0 ;\r
+ Int_t iHMPID = 1;\r
+ Int_t iSHIL = 1;\r
+ Int_t iT0 = 1;\r
+ Int_t iTOF = 1;\r
+ Int_t iTPC = 1;\r
+ Int_t iTRD = 1;\r
+ Int_t iVZERO = 1;\r
+ Int_t iZDC = 1;\r
+ \r
+\r
+ //=================== Alice BODY parameters =============================\r
+ AliBODY *BODY = new AliBODY("BODY", "Alice envelop");\r
+\r
+\r
+ if (iMAG)\r
+ {\r
+ //=================== MAG parameters ============================\r
+ // --- Start with Magnet since detector layouts may be depending ---\r
+ // --- on the selected Magnet dimensions ---\r
+ AliMAG *MAG = new AliMAG("MAG", "Magnet");\r
+ }\r
+\r
+\r
+ if (iABSO)\r
+ {\r
+ //=================== ABSO parameters ============================\r
+ AliABSO *ABSO = new AliABSOv3("ABSO", "Muon Absorber");\r
+ }\r
+\r
+ if (iDIPO)\r
+ {\r
+ //=================== DIPO parameters ============================\r
+\r
+ AliDIPO *DIPO = new AliDIPOv3("DIPO", "Dipole version 3");\r
+ }\r
+\r
+ if (iHALL)\r
+ {\r
+ //=================== HALL parameters ============================\r
+\r
+ AliHALL *HALL = new AliHALLv3("HALL", "Alice Hall");\r
+ }\r
+\r
+\r
+ if (iFRAME)\r
+ {\r
+ //=================== FRAME parameters ============================\r
+\r
+ AliFRAMEv2 *FRAME = new AliFRAMEv2("FRAME", "Space Frame");\r
+ if (geo == kHoles) FRAME->SetHoles(1);\r
+ else FRAME->SetHoles(0);\r
+\r
+ }\r
+\r
+ if (iSHIL)\r
+ {\r
+ //=================== SHIL parameters ============================\r
+\r
+ AliSHIL *SHIL = new AliSHILv3("SHIL", "Shielding Version 3");\r
+ }\r
+\r
+\r
+ if (iPIPE)\r
+ {\r
+ //=================== PIPE parameters ============================\r
+\r
+ AliPIPE *PIPE = new AliPIPEv3("PIPE", "Beam Pipe");\r
+ }\r
+ \r
+ if(iITS) {\r
+\r
+ //=================== ITS parameters ============================\r
+ //\r
+ AliITS *ITS = new AliITSv11Hybrid("ITS","ITS v11Hybrid");\r
+ // AliITSvPPRasymmFMD *ITS = new AliITSvPPRasymmFMD("ITS","New ITS PPR detailed version with asymmetric services");\r
+ }\r
+\r
+ if (iTPC)\r
+ {\r
+ //============================ TPC parameters =====================\r
+ AliTPC *TPC = new AliTPCv2("TPC", "Default");\r
+ }\r
+\r
+\r
+ if (iTOF) {\r
+ //=================== TOF parameters ============================\r
+ AliTOF *TOF = new AliTOFv6T0("TOF", "normal TOF");\r
+ }\r
+\r
+\r
+ if (iHMPID)\r
+ {\r
+ //=================== HMPID parameters ===========================\r
+ // AliHMPID *HMPID = new AliHMPIDv2("HMPID", "normal HMPID");\r
+ AliHMPID *HMPID = new AliHMPIDv3("HMPID", "normal HMPID");\r
+\r
+ }\r
+\r
+\r
+ if (iZDC)\r
+ {\r
+ //=================== ZDC parameters ============================\r
+\r
+ AliZDC *ZDC = new AliZDCv3("ZDC", "normal ZDC");\r
+ }\r
+\r
+ if (iTRD)\r
+ {\r
+ //=================== TRD parameters ============================\r
+\r
+ AliTRD *TRD = new AliTRDv1("TRD", "TRD slow simulator");\r
+ AliTRDgeometry *geoTRD = TRD->GetGeometry();\r
+ // Partial geometry: modules at 0,8,9,17\r
+ // Partial geometry: modules at 1,7,10,16 expected for 2009\r
+ // starting at 3h in positive direction\r
+\r
+\r
+ geoTRD->SetSMstatus(2,0);\r
+ geoTRD->SetSMstatus(3,0);\r
+ geoTRD->SetSMstatus(4,0);\r
+ geoTRD->SetSMstatus(5,0);\r
+ geoTRD->SetSMstatus(6,0);\r
+ geoTRD->SetSMstatus(11,0);\r
+ geoTRD->SetSMstatus(12,0);\r
+ geoTRD->SetSMstatus(13,0);\r
+ geoTRD->SetSMstatus(14,0);\r
+ geoTRD->SetSMstatus(15,0);\r
+\r
+\r
+\r
+\r
+ /*\r
+ // Partial geometry: modules at 2,3,4,6,11,12,14,15\r
+ // starting at 6h in positive direction\r
+ //Hole on SM 13 and 14 for PHOS\r
+ geoTRD->SetSMstatus(13,0);\r
+ geoTRD->SetSMstatus(14,0);\r
+ geoTRD->SetSMstatus(15,0);\r
+ */\r
+\r
+ }\r
+\r
+ if (iFMD)\r
+ {\r
+ //=================== FMD parameters ============================\r
+ AliFMD *FMD = new AliFMDv1("FMD", "normal FMD");\r
+ }\r
+\r
+ if (iMUON)\r
+ {\r
+ //=================== MUON parameters ===========================\r
+ // New MUONv1 version (geometry defined via builders)\r
+ AliMUON *MUON = new AliMUONv1("MUON", "default");\r
+ }\r
+ //=================== PHOS parameters ===========================\r
+\r
+ if (iPHOS)\r
+ {\r
+ AliPHOS *PHOS = new AliPHOSv1("PHOS", "IHEP");\r
+ // Default for cold phos...\r
+ }\r
+\r
+\r
+ if (iPMD)\r
+ {\r
+ //=================== PMD parameters ============================\r
+ AliPMD *PMD = new AliPMDv1("PMD", "normal PMD");\r
+ }\r
+\r
+ if (iT0)\r
+ {\r
+ //=================== T0 parameters ============================\r
+ AliT0 *T0 = new AliT0v1("T0", "T0 Detector");\r
+ }\r
+\r
+ if (iEMCAL)\r
+ {\r
+ //=================== EMCAL parameters ============================\r
+ AliEMCAL *EMCAL = new AliEMCALv2("EMCAL", "EMCAL_COMPLETE");\r
+ }\r
+\r
+ if (iACORDE)\r
+ {\r
+ //=================== CRT parameters ============================\r
+ AliACORDE *ACORDE = new AliACORDEv1("CRT", "normal ACORDE");\r
+ }\r
+\r
+ if (iVZERO)\r
+ {\r
+ //=================== CRT parameters ============================\r
+ AliVZERO *VZERO = new AliVZEROv7("VZERO", "normal VZERO");\r
+ }\r
+}\r
+\r
+// PYTHIA\r
+\r
+AliGenPythia *PythiaHVQ(PDC07Proc_t proc) {\r
+//*******************************************************************//\r
+// Configuration file for charm / beauty generation with PYTHIA //\r
+// //\r
+// The parameters have been tuned in order to reproduce the inclusive//\r
+// heavy quark pt distribution given by the NLO pQCD calculation by //\r
+// Mangano, Nason and Ridolfi. //\r
+// //\r
+// For details and for the NORMALIZATION of the yields see: //\r
+// N.Carrer and A.Dainese, //\r
+// "Charm and beauty production at the LHC", //\r
+// ALICE-INT-2003-019, [arXiv:hep-ph/0311225]; //\r
+// PPR Chapter 6.6, CERN/LHCC 2005-030 (2005). //\r
+//*******************************************************************//\r
+ AliGenPythia * gener = 0x0;\r
+\r
+ switch(proc) {\r
+ case kCharmPbPb5500:\r
+ comment = comment.Append(" Charm in Pb-Pb at 5.5 TeV");\r
+ gener = new AliGenPythia(nEvts);\r
+ gener->SetProcess(kPyCharmPbPbMNR);\r
+ gener->SetStrucFunc(kCTEQ4L);\r
+ gener->SetPtHard(2.1,-1.0);\r
+ gener->SetEnergyCMS(5500.);\r
+ gener->SetNuclei(208,208);\r
+ break;\r
+ case kCharmpPb8800:\r
+ comment = comment.Append(" Charm in p-Pb at 8.8 TeV");\r
+ gener = new AliGenPythia(nEvts);\r
+ gener->SetProcess(kPyCharmpPbMNR);\r
+ gener->SetStrucFunc(kCTEQ4L);\r
+ gener->SetPtHard(2.1,-1.0);\r
+ gener->SetEnergyCMS(8800.);\r
+ gener->SetProjectile("P",1,1);\r
+ gener->SetTarget("Pb",208,82);\r
+ break;\r
+ case kCharmpp14000:\r
+ comment = comment.Append(" Charm in pp at 14 TeV");\r
+ gener = new AliGenPythia(nEvts);\r
+ gener->SetProcess(kPyCharmppMNR);\r
+ gener->SetStrucFunc(kCTEQ4L);\r
+ gener->SetPtHard(2.1,-1.0);\r
+ gener->SetEnergyCMS(14000.);\r
+ break;\r
+ case kCharmpp14000wmi:\r
+ comment = comment.Append(" Charm in pp at 14 TeV with mult. interactions");\r
+ gener = new AliGenPythia(-1);\r
+ gener->SetProcess(kPyCharmppMNRwmi);\r
+ gener->SetStrucFunc(kCTEQ5L);\r
+ gener->SetPtHard(ptHardMin,ptHardMax);\r
+ gener->SetEnergyCMS(14000.);\r
+ break;\r
+ case kD0PbPb5500:\r
+ comment = comment.Append(" D0 in Pb-Pb at 5.5 TeV");\r
+ gener = new AliGenPythia(nEvts);\r
+ gener->SetProcess(kPyD0PbPbMNR);\r
+ gener->SetStrucFunc(kCTEQ4L);\r
+ gener->SetPtHard(2.1,-1.0);\r
+ gener->SetEnergyCMS(5500.);\r
+ gener->SetNuclei(208,208);\r
+ break;\r
+ case kD0pPb8800:\r
+ comment = comment.Append(" D0 in p-Pb at 8.8 TeV");\r
+ gener = new AliGenPythia(nEvts);\r
+ gener->SetProcess(kPyD0pPbMNR);\r
+ gener->SetStrucFunc(kCTEQ4L);\r
+ gener->SetPtHard(2.1,-1.0);\r
+ gener->SetEnergyCMS(8800.);\r
+ gener->SetProjectile("P",1,1);\r
+ gener->SetTarget("Pb",208,82);\r
+ break;\r
+ case kD0pp14000:\r
+ comment = comment.Append(" D0 in pp at 14 TeV");\r
+ gener = new AliGenPythia(nEvts);\r
+ gener->SetProcess(kPyD0ppMNR);\r
+ gener->SetStrucFunc(kCTEQ4L);\r
+ gener->SetPtHard(2.1,-1.0);\r
+ gener->SetEnergyCMS(14000.);\r
+ break;\r
+ case kDPlusPbPb5500:\r
+ comment = comment.Append(" DPlus in Pb-Pb at 5.5 TeV");\r
+ gener = new AliGenPythia(nEvts);\r
+ gener->SetProcess(kPyDPlusPbPbMNR);\r
+ gener->SetStrucFunc(kCTEQ4L);\r
+ gener->SetPtHard(2.1,-1.0);\r
+ gener->SetEnergyCMS(5500.);\r
+ gener->SetNuclei(208,208);\r
+ break;\r
+ case kDPluspPb8800:\r
+ comment = comment.Append(" DPlus in p-Pb at 8.8 TeV");\r
+ gener = new AliGenPythia(nEvts);\r
+ gener->SetProcess(kPyDPluspPbMNR);\r
+ gener->SetStrucFunc(kCTEQ4L);\r
+ gener->SetPtHard(2.1,-1.0);\r
+ gener->SetEnergyCMS(8800.);\r
+ gener->SetProjectile("P",1,1);\r
+ gener->SetTarget("Pb",208,82);\r
+ break;\r
+ case kDPluspp14000:\r
+ comment = comment.Append(" DPlus in pp at 14 TeV");\r
+ gener = new AliGenPythia(nEvts);\r
+ gener->SetProcess(kPyDPlusppMNR);\r
+ gener->SetStrucFunc(kCTEQ4L);\r
+ gener->SetPtHard(2.1,-1.0);\r
+ gener->SetEnergyCMS(14000.);\r
+ break;\r
+ case kBeautyPbPb5500:\r
+ comment = comment.Append(" Beauty in Pb-Pb at 5.5 TeV");\r
+ gener = new AliGenPythia(nEvts);\r
+ gener->SetProcess(kPyBeautyPbPbMNR);\r
+ gener->SetStrucFunc(kCTEQ4L);\r
+ gener->SetPtHard(2.75,-1.0);\r
+ gener->SetEnergyCMS(5500.);\r
+ gener->SetNuclei(208,208);\r
+ break;\r
+ case kBeautypPb8800:\r
+ comment = comment.Append(" Beauty in p-Pb at 8.8 TeV");\r
+ gener = new AliGenPythia(nEvts);\r
+ gener->SetProcess(kPyBeautypPbMNR);\r
+ gener->SetStrucFunc(kCTEQ4L);\r
+ gener->SetPtHard(2.75,-1.0);\r
+ gener->SetEnergyCMS(8800.);\r
+ gener->SetProjectile("P",1,1);\r
+ gener->SetTarget("Pb",208,82);\r
+ break;\r
+ case kBeautypp14000:\r
+ comment = comment.Append(" Beauty in pp at 14 TeV");\r
+ gener = new AliGenPythia(nEvts);\r
+ gener->SetProcess(kPyBeautyppMNR);\r
+ gener->SetStrucFunc(kCTEQ4L);\r
+ gener->SetPtHard(2.75,-1.0);\r
+ gener->SetEnergyCMS(14000.);\r
+ break;\r
+ case kBeautypp14000wmi:\r
+ comment = comment.Append(" Beauty in pp at 14 TeV with mult. interactions");\r
+ gener = new AliGenPythia(-1);\r
+ gener->SetProcess(kPyBeautyppMNRwmi);\r
+ gener->SetStrucFunc(kCTEQ5L);\r
+ gener->SetPtHard(ptHardMin,ptHardMax);\r
+ gener->SetEnergyCMS(14000.);\r
+ break;\r
+ }\r
+\r
+ return gener;\r
+}\r
+\r
+AliGenerator* MbCocktail()\r
+{\r
+ comment = comment.Append(" pp at 14 TeV: Pythia low-pt, no heavy quarks + J/Psi from parameterisation");\r
+ AliGenCocktail * gener = new AliGenCocktail();\r
+ gener->UsePerEventRates();\r
+ \r
+// Pythia\r
+\r
+ AliGenPythia* pythia = new AliGenPythia(-1); \r
+ pythia->SetMomentumRange(0, 999999.);\r
+ pythia->SetThetaRange(0., 180.);\r
+ pythia->SetYRange(-12.,12.);\r
+ pythia->SetPtRange(0,1000.);\r
+ pythia->SetProcess(kPyMb);\r
+ pythia->SetEnergyCMS(14000.);\r
+ pythia->SwitchHFOff();\r
+ \r
+// J/Psi parameterisation\r
+\r
+ AliGenParam* jpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF scaled", "Jpsi");\r
+ jpsi->SetPtRange(0.,100.);\r
+ jpsi->SetYRange(-8., 8.);\r
+ jpsi->SetPhiRange(0., 360.);\r
+ jpsi->SetForceDecay(kAll);\r
+\r
+ gener->AddGenerator(pythia, "Pythia", 1.);\r
+ gener->AddGenerator(jpsi, "J/Psi", 8.e-4); \r
+ \r
+ return gener;\r
+}\r
+\r
+AliGenerator* PyMbTriggered(Int_t pdg)\r
+{\r
+ AliGenPythia* pythia = new AliGenPythia(-1); \r
+ pythia->SetMomentumRange(0, 999999.);\r
+ pythia->SetThetaRange(0., 180.);\r
+ pythia->SetYRange(-12.,12.);\r
+ pythia->SetPtRange(0,1000.);\r
+ pythia->SetProcess(kPyMb);\r
+ pythia->SetEnergyCMS(14000.);\r
+ pythia->SetTriggerParticle(pdg, 0.9);\r
+ return pythia;\r
+}\r
+\r
+void ProcessEnvironmentVars()\r
+{ \r
+ cout << "######################################" << endl;\r
+ cout << "## Processing environment variables ##" << endl;\r
+ cout << "######################################" << endl;\r
+ \r
+ // Run Number\r
+ if (gSystem->Getenv("DC_RUN")) {\r
+ runNumber = atoi(gSystem->Getenv("DC_RUN"));\r
+ }\r
+ //cout<<"Run number "<<runNumber<<endl;\r
+\r
+ // Event Number\r
+ if (gSystem->Getenv("DC_EVENT")) {\r
+ eventNumber = atoi(gSystem->Getenv("DC_EVENT"));\r
+ }\r
+ //cout<<"Event number "<<eventNumber<<endl;\r
+\r
+ // Random Number seed\r
+ if (gSystem->Getenv("CONFIG_SEED")) {\r
+ seed = atoi(gSystem->Getenv("CONFIG_SEED"));\r
+ }\r
+ else if(gSystem->Getenv("DC_EVENT") && gSystem->Getenv("DC_RUN")){\r
+ seed = runNumber * 100000 + eventNumber;\r
+ }\r
+\r
+ gRandom->SetSeed(seed);\r
+ \r
+ cout<<"////////////////////////////////////////////////////////////////////////////////////"<<endl;\r
+ cout<<"Seed for random number generation= "<< seed<<" "<< gRandom->GetSeed()<<endl; \r
+ cout<<"////////////////////////////////////////////////////////////////////////////////////"<<endl;\r
+\r
+ // Run type\r
+ if (gSystem->Getenv("DC_RUN_TYPE")) {\r
+ cout<<"run type "<<gSystem->Getenv("DC_RUN_TYPE")<<endl;\r
+ for (Int_t iRun = 0; iRun < kRunMax; iRun++) {\r
+ if (strcmp(gSystem->Getenv("DC_RUN_TYPE"), pprRunName[iRun])==0) {\r
+ proc = (PDC07Proc_t)iRun;\r
+ }\r
+ }\r
+ }\r
+ else \r
+ cout << "Environment variable DC_RUN_TYPE is not defined" << endl;\r
+\r
+ if (gSystem->Getenv("ECMS"))\r
+ eCMS = atof(gSystem->Getenv("ECMS"));\r
+ if (gSystem->Getenv("PTHARDMIN"))\r
+ ptHardMin = atof(gSystem->Getenv("PTHARDMIN"));\r
+ if (gSystem->Getenv("PTHARDMAX"))\r
+ ptHardMax = atof(gSystem->Getenv("PTHARDMAX"));\r
+ if (gSystem->Getenv("PTGAMMAPI0MIN"))\r
+ ptGammaPi0Min = atof(gSystem->Getenv("PTGAMMAPI0MIN"));\r
+ if (gSystem->Getenv("QUENCHING"))\r
+ iquenching = atof(gSystem->Getenv("QUENCHING"));\r
+ if (gSystem->Getenv("QHAT"))\r
+ qhat = atof(gSystem->Getenv("QHAT"));\r
+\r
+ cout<<">> Run type set to "<<pprRunName[proc]<<endl;\r
+ cout<<">> Collision energy set to "<<eCMS <<endl;\r
+ cout<<">> ptHard limits: "<<ptHardMin<<" to " <<ptHardMax<<" GeV"<<endl;\r
+ cout<<">> pt gamma/pi0 threshold "<< ptGammaPi0Min<<" GeV "<<endl;\r
+ cout<<">> quenching on? "<< iquenching<<" qhat "<<qhat<<endl;\r
+\r
+ cout << "######################################" << endl;\r
+ cout << "######################################" << endl;\r
+}\r