0cba2f78828f650722692a479a1f4e3fcbe5f6b1
[u/mrichter/AliRoot.git] / prod / LHC08d8 / Config.C
1 // One can use the configuration macro in compiled mode by\r
2 // root [0] gSystem->Load("libgeant321");
3 // root [0] gSystem->SetIncludePath("-I$ROOTSYS/include -I$ALICE_ROOT/include\
4 //                   -I$ALICE_ROOT -I$ALICE/geant3/TGeant3");
5 // root [0] .x grun.C(1,"Config_PDC07.C++")
6
7 #if !defined(__CINT__) || defined(__MAKECINT__)
8 #include <Riostream.h>
9 #include <TRandom.h>
10 #include <TDatime.h>
11 #include <TSystem.h>
12 #include <TVirtualMC.h>
13 #include <TGeant3TGeo.h>
14 #include "EVGEN/AliGenCocktail.h"
15 #include "EVGEN/AliGenParam.h"
16 #include "EVGEN/AliGenMUONlib.h"
17 #include "STEER/AliRunLoader.h"
18 #include "STEER/AliRun.h"
19 #include "STEER/AliConfig.h"
20 #include "PYTHIA6/AliDecayerPythia.h"
21 #include "PYTHIA6/AliGenPythia.h"
22 #include "STEER/AliMagWrapCheb.h"
23 #include "STRUCT/AliBODY.h"
24 #include "STRUCT/AliMAG.h"
25 #include "STRUCT/AliABSOv3.h"
26 #include "STRUCT/AliDIPOv3.h"
27 #include "STRUCT/AliHALLv3.h"
28 #include "STRUCT/AliFRAMEv2.h"
29 #include "STRUCT/AliSHILv3.h"
30 #include "STRUCT/AliPIPEv3.h"
31 #include "ITS/AliITSv11Hybrid.h"
32 #include "TPC/AliTPCv2.h"
33 #include "TOF/AliTOFv6T0.h"
34 #include "HMPID/AliHMPIDv3.h"
35 #include "ZDC/AliZDCv3.h"
36 #include "TRD/AliTRDv1.h"
37 #include "TRD/AliTRDgeometry.h"
38 #include "FMD/AliFMDv1.h"
39 #include "MUON/AliMUONv1.h"
40 #include "PHOS/AliPHOSv1.h"
41 #include "PHOS/AliPHOSSimParam.h"
42 #include "PMD/AliPMDv1.h"
43 #include "T0/AliT0v1.h"
44 #include "EMCAL/AliEMCALv2.h"
45 #include "ACORDE/AliACORDEv1.h"
46 #include "VZERO/AliVZEROv7.h"
47 #endif
48
49
50 enum PDC07Proc_t
51 {
52 //--- Heavy Flavour Production ---
53   kCharmPbPb5500,  kCharmpPb8800,  kCharmpp14000,  kCharmpp14000wmi,
54   kD0PbPb5500,     kD0pPb8800,     kD0pp14000,
55   kDPlusPbPb5500,  kDPluspPb8800,  kDPluspp14000,
56   kBeautyPbPb5500, kBeautypPb8800, kBeautypp14000, kBeautypp14000wmi,
57   // -- Pythia Mb
58   kPyMbNoHvq, kPyOmegaPlus, kPyOmegaMinus, kPyJetJet, kBeautyppwmiJet,
59   kPyGammaJetPHOS, kPyJetJetPHOS, kPyJetJetPHOSv2, kPyGammaBremsPHOS,
60   kPyGammaJetEMCAL, kPyJetJetEMCAL, kPyGammaBremsEMCAL, kRunMax
61 };
62
63 const char * pprRunName[] = {
64   "kCharmPbPb5500",  "kCharmpPb8800",  "kCharmpp14000",  "kCharmpp14000wmi",
65   "kD0PbPb5500",     "kD0pPb8800",     "kD0pp14000",
66   "kDPlusPbPb5500",  "kDPluspPb8800",  "kDPluspp14000",
67   "kBeautyPbPb5500", "kBeautypPb8800", "kBeautypp14000", "kBeautypp14000wmi",
68   "kPyMbNoHvq", "kPyOmegaPlus", "kPyOmegaMinus", "kPyJetJet", "kBeautyppwmiJet",
69   "kPyGammaJetPHOS", "kPyJetJetPHOS",  "kPyJetJetPHOSv2", "kPyGammaBremsPHOS",
70   "kPyGammaJetEMCAL", "kPyJetJetEMCAL", "kPyGammaBremsEMCAL"
71 };
72
73
74 //--- Decay Mode ---
75 enum DecayHvFl_t
76 {
77   kNature,  kHadr, kSemiEl, kSemiMu
78 };
79 //--- Magnetic Field ---
80 enum Mag_t
81 {
82   kNoField, k5kG, kFieldMax
83 };
84 //--- Trigger config ---
85 enum TrigConf_t
86 {
87     kDefaultPPTrig, kDefaultPbPbTrig
88 };
89
90 const char * TrigConfName[] = {
91     "p-p","Pb-Pb"
92 };
93
94 //--- Functions ---
95 class AliGenPythia ;
96 AliGenPythia *PythiaHVQ(PDC07Proc_t proc);
97 AliGenPythia *PythiaHard(PDC07Proc_t proc);
98 AliGenerator *MbCocktail();
99 AliGenerator *PyMbTriggered(Int_t pdg);
100 void ProcessEnvironmentVars();
101
102 // This part for configuration
103 static DecayHvFl_t   decHvFl     = kNature;
104 static Mag_t         mag         = k5kG;
105 static TrigConf_t    trig        = kDefaultPPTrig; // default pp trigger configuration
106 static Int_t         runNumber   = 0;
107 static Int_t         eventNumber = 0;
108 //========================//
109 // Set Random Number seed //
110 //========================//
111 TDatime dt;
112 static UInt_t seed    = dt.Get();
113
114 // nEvts = -1  : you get 1 QQbar pair and all the fragmentation and
115 //               decay chain
116 // nEvts = N>0 : you get N charm / beauty Hadrons
117 Int_t nEvts = -1;
118 // stars = kTRUE : all heavy resonances and their decay stored
119 //       = kFALSE: only final heavy hadrons and their decays stored
120 Bool_t stars = kTRUE;
121
122 // To be used only with kCharmpp14000wmi and kBeautypp14000wmi
123 // To get a "reasonable" agreement with MNR results, events have to be
124 // generated with the minimum ptHard set to 2.76 GeV.
125 // To get a "perfect" agreement with MNR results, events have to be
126 // generated in four ptHard bins with the following relative
127 // normalizations:
128 //  CHARM
129 // 2.76-3 GeV: 25%
130 //    3-4 GeV: 40%
131 //    4-8 GeV: 29%
132 //     >8 GeV:  6%
133 //  BEAUTY
134 // 2.76-4 GeV:  5%
135 //    4-6 GeV: 31%
136 //    6-8 GeV: 28%
137 //     >8 GeV: 36%
138
139
140 Float_t eCMS=5500;
141 PDC07Proc_t proc = kPyJetJet;
142 //PDC07Proc_t proc = kPyGammaJetEMCAL;
143 //PDC07Proc_t proc = kBeautyppwmiJet;
144 Int_t nPtHardGJBins = 11;
145 Int_t ptHardGJLo[] = {5 ,10,20,30,40,50,60,70,80,90 ,100};
146 Int_t ptHardGJHi[] = {10,20,30,40,50,60,70,80,90,100,-1};
147 Int_t nPtHardJJBins = 17;
148 Int_t ptHardJJLo[] = {2 ,12,16,20,24,29,35,41,50,60,72,86 ,104,124,149,179,215};
149 Int_t ptHardJJHi[] = {12,16,20,24,29,35,41,50,60,72,86,104,124,149,179,215,258};
150 Float_t ptHardMin =  10.;
151 Float_t ptHardMax =  20.;
152 Float_t ptGammaPi0Min = 0.;
153 Int_t iquenching = 0;
154 Float_t qhat = 50.; //for  iquenching equal 1 and 4
155 //Float_t medlength = 6.; // for iquenching equal to 4
156 Int_t year = 0; //ALICE geometry configuration
157 Bool_t pi0FragPhotonInDet = kFALSE ;
158
159 // Comment line
160 static TString comment;
161
162 void Config()
163 {
164
165         // Get settings from environment variables
166         ProcessEnvironmentVars();
167         
168         gSystem->Load("liblhapdf.so");      // Parton density functions
169         gSystem->Load("libEGPythia6.so");   // TGenerator interface
170         gSystem->Load("libpythia6.so");     // Pythia
171         gSystem->Load("libAliPythia6.so");  // ALICE specific implementations
172         
173         // libraries required by geant321
174 #if defined(__CINT__)
175         gSystem->Load("liblhapdf");
176         gSystem->Load("libEGPythia6");
177         gSystem->Load("libpythia6");
178         gSystem->Load("libAliPythia6");
179         gSystem->Load("libgeant321");
180 #endif
181         
182         new TGeant3TGeo("C++ Interface to Geant3");
183         
184         // Output every 100 tracks
185         ((TGeant3*)gMC)->SetSWIT(4,100);
186         
187         //=======================================================================
188         
189         // Run loader
190         AliRunLoader* rl=0x0;
191         rl = AliRunLoader::Open("galice.root",
192                                                         AliConfig::GetDefaultEventFolderName(),
193                                                         "recreate");
194         if (rl == 0x0)
195     {
196                 gAlice->Fatal("Config.C","Can not instatiate the Run Loader");
197                 return;
198     }
199         rl->SetCompressionLevel(2);
200         rl->SetNumberOfEventsPerFile(1000);
201         gAlice->SetRunLoader(rl);
202         
203         // Set the trigger configuration
204         gAlice->SetTriggerDescriptor(TrigConfName[trig]);
205         cout<<"Trigger configuration is set to  "<<TrigConfName[trig]<<endl;
206         
207         //
208         //=======================================================================
209         // ************* STEERING parameters FOR ALICE SIMULATION **************
210         // --- Specify event type to be tracked through the ALICE setup
211         // --- All positions are in cm, angles in degrees, and P and E in GeV
212         
213         
214     gMC->SetProcess("DCAY",1);
215     gMC->SetProcess("PAIR",1);
216     gMC->SetProcess("COMP",1);
217     gMC->SetProcess("PHOT",1);
218     gMC->SetProcess("PFIS",0);
219     gMC->SetProcess("DRAY",0);
220     gMC->SetProcess("ANNI",1);
221     gMC->SetProcess("BREM",1);
222     gMC->SetProcess("MUNU",1);
223     gMC->SetProcess("CKOV",1);
224     gMC->SetProcess("HADR",1);
225     gMC->SetProcess("LOSS",2);
226     gMC->SetProcess("MULS",1);
227     gMC->SetProcess("RAYL",1);
228         
229     Float_t cut = 1.e-3;        // 1MeV cut by default
230     Float_t tofmax = 1.e10;
231         
232     gMC->SetCut("CUTGAM", cut);
233     gMC->SetCut("CUTELE", cut);
234     gMC->SetCut("CUTNEU", cut);
235     gMC->SetCut("CUTHAD", cut);
236     gMC->SetCut("CUTMUO", cut);
237     gMC->SetCut("BCUTE",  cut);
238     gMC->SetCut("BCUTM",  cut);
239     gMC->SetCut("DCUTE",  cut);
240     gMC->SetCut("DCUTM",  cut);
241     gMC->SetCut("PPCUTM", cut);
242     gMC->SetCut("TOFMAX", tofmax);
243         
244         //======================//
245         // Set External decayer //
246         //======================//
247         TVirtualMCDecayer* decayer = new AliDecayerPythia();
248         // DECAYS
249         //
250         switch(decHvFl) {
251                 case kNature:
252                         decayer->SetForceDecay(kNeutralPion);
253                         break;
254                 case kHadr:
255                         decayer->SetForceDecay(kHadronicD);
256                         break;
257                 case kSemiEl:
258                         decayer->SetForceDecay(kSemiElectronic);
259                         break;
260                 case kSemiMu:
261                         decayer->SetForceDecay(kSemiMuonic);
262                         break;
263         }
264         decayer->Init();
265         gMC->SetExternalDecayer(decayer);
266         //if(proc == kPyJetJetPHOSv2) // in this case we decay the pi0
267         //      decayer->SetForceDecay(kNeutralPion);
268         
269         //=========================//
270         // Generator Configuration //
271         //=========================//
272         AliGenPythia* gener = 0x0;
273         
274         if (proc <=   kBeautypp14000wmi) {
275                 AliGenPythia *pythia = PythiaHVQ(proc);
276                 // FeedDown option
277                 pythia->SetFeedDownHigherFamily(kFALSE);
278                 // Stack filling option
279                 if(!stars) pythia->SetStackFillOpt(AliGenPythia::kParentSelection);
280                 // Set Count mode
281                 if(nEvts>0) pythia->SetCountMode(AliGenPythia::kCountParents);
282                 //
283                 // DECAYS
284                 //
285                 switch(decHvFl) {
286                         case kNature:
287                                 pythia->SetForceDecay(kNeutralPion);
288                                 break;
289                         case kHadr:
290                                 pythia->SetForceDecay(kHadronicD);
291                                 break;
292                         case kSemiEl:
293                                 pythia->SetForceDecay(kSemiElectronic);
294                                 break;
295                         case kSemiMu:
296                                 pythia->SetForceDecay(kSemiMuonic);
297                                 break;
298                 }
299                 //
300                 // GEOM & KINE CUTS
301                 //
302                 pythia->SetMomentumRange(0,99999999);
303                 pythia->SetPhiRange(0., 360.);
304                 pythia->SetThetaRange(0,180);
305                 switch(ycut) {
306                         case kFull:
307                                 pythia->SetYRange(-12,12);
308                                 break;
309                         case kBarrel:
310                                 pythia->SetYRange(-2,2);
311                                 break;
312                         case kMuonArm:
313                                 pythia->SetYRange(1,6);
314                                 break;
315                 }
316                 gener = pythia;
317         } else if (proc == kPyMbNoHvq) {
318                 gener = MbCocktail();
319         } else if (proc == kPyOmegaMinus) {
320                 gener = PyMbTriggered(3334);
321         } else if (proc == kPyOmegaPlus) {
322                 gener = PyMbTriggered(-3334);
323         } else if (proc <= kPyGammaBremsEMCAL) {
324                 AliGenPythia *pythia = PythiaHard(proc);
325                 // FeedDown option
326                 pythia->SetFeedDownHigherFamily(kFALSE);
327                 // Set Count mode
328                 if(nEvts>0) pythia->SetCountMode(AliGenPythia::kCountParents);
329                 
330                 //
331                 // GEOM & KINE CUTS
332                 //
333                 pythia->SetMomentumRange(0,99999999);
334 //              if(proc == kPyJetJetPHOSv2)
335 //                      pythia->SetForceDecay(kNeutralPion);
336 //              else
337 //                      pythia->SetForceDecay(kAll);
338
339                 pythia->SetForceDecay(kNeutralPion);
340                 pythia->SetPycellParameters(2., 274, 432, 0., 4., 5., 1.0);
341                 pythia->SetGluonRadiation(1,1); //ISR and FRS effects on.
342                 pythia->SetPtKick(3.); // set the intrinsic kt to about 5.3 GeV/c
343                 gener = pythia;
344         }
345         
346         
347         // PRIMARY VERTEX
348         
349         gener->SetOrigin(0., 0., 0.);    // vertex position
350         
351         // Size of the interaction diamond
352         // Longitudinal
353         Float_t           sigmaz  = 7.55 / TMath::Sqrt(2.); // [cm] 14  TeV (p-p), 5.5 TeV (Pb-Pb)
354         if(eCMS == 10000) sigmaz  = 5.4  / TMath::Sqrt(2.); // [cm] 10  TeV
355         if(eCMS == 900)   sigmaz  = 10.5 / TMath::Sqrt(2.); // [cm] 0.9 TeV
356         cout<<">>> SigmaZ x srqt(2.) "<< sigmaz * TMath::Sqrt(2.) << endl;
357         
358     // Transverse
359         Float_t betast  = 10;                 // beta* [m]
360         Float_t eps     = 3.75e-6;            // emittance [m]
361         Float_t gamma   = eCMS / 2.0 / 0.938272;   // relativistic gamma [1]
362         Float_t sigmaxy = TMath::Sqrt(eps * betast / gamma) / TMath::Sqrt(2.) * 100.;  // [cm]
363         printf("\n \n Diamond size x-y: %10.3e z: %10.3e\n \n", sigmaxy, sigmaz);
364         
365         gener->SetSigma(sigmaxy, sigmaxy, sigmaz);      // Sigma in (X,Y,Z) (cm) on IP position
366         gener->SetCutVertexZ(3.);        // Truncate at 3 sigma
367         gener->SetVertexSmear(kPerEvent);
368         
369         gener->Init();
370         
371         //Quenching
372         gener->SetQuench(iquenching);
373         if(iquenching == 1){ //quenching weights
374                 Float_t k = 6e5*(qhat/1.7) ; //qhat=1.7, k = 6e5, default  value
375                 AliPythia::Instance()->InitQuenching(0.,0.1,k,0,0.95,6);
376         }
377         else if(iquenching == 4){ //q-pythia
378                 //Not yet implemented
379                 gener->SetQhat(qhat); //   !transport coefficient (GeV2/fm)
380         }
381
382   // FIELD
383   //
384   AliMagWrapCheb* field = 0x0;
385   if (mag == kNoField) {
386     comment = comment.Append(" | L3 field 0.0 T");
387     field = new AliMagWrapCheb("Maps","Maps", 2, 0., 10., AliMagWrapCheb::k2kG);
388   } else if (mag == k5kG) {
389     comment = comment.Append(" | L3 field 0.5 T");
390     field = new AliMagWrapCheb("Maps","Maps", 2, 1., 10., AliMagWrapCheb::k5kG,
391                                kTRUE,"$(ALICE_ROOT)/data/maps/mfchebKGI_sym.root");
392   }
393   printf("\n \n Comment: %s \n \n", comment.Data());
394
395   rl->CdGAFile();
396   gAlice->SetField(field);
397
398         
399         Int_t iABSO   = 1;
400         Int_t iACORDE = 0;
401         Int_t iDIPO   = 1;
402         Int_t iEMCAL  = 1;
403         Int_t iFMD    = 1;
404         Int_t iFRAME  = 1;
405         Int_t iHALL   = 1;
406         Int_t iITS    = 1;
407         Int_t iMAG    = 1;
408         Int_t iMUON   = 1;
409         Int_t iPHOS   = 1;
410         Int_t iPIPE   = 1;
411         Int_t iPMD    = 1;
412         Int_t iHMPID  = 1;
413         Int_t iSHIL   = 1;
414         Int_t iT0     = 1;
415         Int_t iTOF    = 1;
416         Int_t iTPC    = 1;
417         Int_t iTRD    = 1;
418         Int_t iVZERO  = 1;
419         Int_t iZDC    = 1;
420         
421         
422     //=================== Alice BODY parameters =============================
423     AliBODY *BODY = new AliBODY("BODY", "Alice envelop");
424         
425         
426     if (iMAG)
427     {
428         //=================== MAG parameters ============================
429         // --- Start with Magnet since detector layouts may be depending ---
430         // --- on the selected Magnet dimensions ---
431         AliMAG *MAG = new AliMAG("MAG", "Magnet");
432     }
433         
434         
435     if (iABSO)
436     {
437         //=================== ABSO parameters ============================
438         AliABSO *ABSO = new AliABSOv3("ABSO", "Muon Absorber");
439     }
440         
441     if (iDIPO)
442     {
443         //=================== DIPO parameters ============================
444                 
445         AliDIPO *DIPO = new AliDIPOv3("DIPO", "Dipole version 3");
446     }
447         
448     if (iHALL)
449     {
450         //=================== HALL parameters ============================
451                 
452         AliHALL *HALL = new AliHALLv3("HALL", "Alice Hall");
453     }
454         
455         
456     if (iFRAME)
457     {
458         //=================== FRAME parameters ============================
459                 
460         AliFRAMEv2 *FRAME = new AliFRAMEv2("FRAME", "Space Frame");
461             FRAME->SetHoles(1);
462     }
463         
464     if (iSHIL)
465     {
466         //=================== SHIL parameters ============================
467                 
468         AliSHIL *SHIL = new AliSHILv3("SHIL", "Shielding Version 3");
469     }
470         
471         
472     if (iPIPE)
473     {
474         //=================== PIPE parameters ============================
475                 
476         AliPIPE *PIPE = new AliPIPEv3("PIPE", "Beam Pipe");
477     }
478         
479     if (iITS)
480         {
481                 //=================== ITS parameters ============================
482                 
483                 AliITS *ITS  = new AliITSv11Hybrid("ITS","ITS v11Hybrid");
484     }
485         
486     if (iTPC)
487     {
488                 //============================ TPC parameters =====================
489
490         AliTPC *TPC = new AliTPCv2("TPC", "Default");
491     }
492         
493         
494     if (iTOF) {
495         //=================== TOF parameters ============================
496
497                 AliTOF *TOF = new AliTOFv6T0("TOF", "normal TOF");
498         }
499         
500         
501     if (iHMPID)
502     {
503         //=================== HMPID parameters ===========================
504  
505                 AliHMPID *HMPID = new AliHMPIDv3("HMPID", "normal HMPID");  
506         }
507         
508         
509     if (iZDC)
510     {
511         //=================== ZDC parameters ============================
512                 
513         AliZDC *ZDC = new AliZDCv3("ZDC", "normal ZDC");
514     }
515         
516     if (iTRD)
517     {
518         //=================== TRD parameters ============================
519                 
520         AliTRD *TRD = new AliTRDv1("TRD", "TRD slow simulator");
521         AliTRDgeometry *geoTRD = TRD->GetGeometry();
522                 if(year == 2009){
523                         // Partial geometry: modules at 0,8,9,17
524                         // Partial geometry: modules at 1,7,10,16 expected for 2009
525                         // starting at 3h in positive direction
526                         
527                         //geoTRD->SetSMstatus(1,0);
528                         //geoTRD->SetSMstatus(2,0);
529                         //geoTRD->SetSMstatus(3,0);
530                         //geoTRD->SetSMstatus(4,0);
531                         //geoTRD->SetSMstatus(5,0);
532                         //geoTRD->SetSMstatus(6,0);
533                         //geoTRD->SetSMstatus(7,0);
534                         
535                         
536                         //geoTRD->SetSMstatus(10,0);
537                         //geoTRD->SetSMstatus(11,0);
538                         //geoTRD->SetSMstatus(12,0);
539                         //geoTRD->SetSMstatus(13,0);
540                         //geoTRD->SetSMstatus(14,0);
541                         //geoTRD->SetSMstatus(15,0);
542                         //geoTRD->SetSMstatus(16,0);
543         
544                         geoTRD->SetSMstatus(2,0);
545                         geoTRD->SetSMstatus(3,0);
546                         geoTRD->SetSMstatus(4,0);
547                         geoTRD->SetSMstatus(5,0);
548                         geoTRD->SetSMstatus(6,0);
549
550                                geoTRD->SetSMstatus(11,0);
551                                geoTRD->SetSMstatus(12,0);
552                                        geoTRD->SetSMstatus(13,0);
553                                        geoTRD->SetSMstatus(14,0);
554                                        geoTRD->SetSMstatus(15,0);
555         }
556     }
557         
558     if (iFMD)
559     {
560         //=================== FMD parameters ============================
561                 AliFMD *FMD = new AliFMDv1("FMD", "normal FMD");
562         }
563         
564     if (iMUON)
565     {
566         //=================== MUON parameters ===========================
567         // New MUONv1 version (geometry defined via builders)
568         AliMUON *MUON = new AliMUONv1("MUON", "default");
569     }
570     //=================== PHOS parameters ===========================
571         
572     if (iPHOS)
573     {
574         AliPHOS *PHOS = new AliPHOSv1("PHOS", "IHEP");
575     }
576         
577         
578     if (iPMD)
579     {
580         //=================== PMD parameters ============================
581         AliPMD *PMD = new AliPMDv1("PMD", "normal PMD");
582     }
583         
584     if (iT0)
585     {
586         //=================== T0 parameters ============================
587         AliT0 *T0 = new AliT0v1("T0", "T0 Detector");
588     }
589         
590     if (iEMCAL)
591     {
592         //=================== EMCAL parameters ============================
593         AliEMCAL *EMCAL = new AliEMCALv2("EMCAL", "EMCAL_COMPLETE");
594     }
595         
596         if (iACORDE)
597     {
598         //=================== ACORDE parameters ============================
599         AliACORDE *ACORDE = new AliACORDEv1("ACORDE", "normal ACORDE");
600     }
601         
602         if (iVZERO)
603     {
604         //=================== VZERO parameters ============================
605         AliVZERO *VZERO = new AliVZEROv7("VZERO", "normal VZERO");
606     }
607         
608         AliLog::Message(AliLog::kInfo, "End of Config", "Config.C", "Config.C", "Config()"," Config.C", __LINE__);
609 }
610
611 //           PYTHIA
612
613 AliGenPythia *PythiaHVQ(PDC07Proc_t proc) {
614 //*******************************************************************//
615 // Configuration file for charm / beauty generation with PYTHIA      //
616 //                                                                   //
617 // The parameters have been tuned in order to reproduce the inclusive//
618 // heavy quark pt distribution given by the NLO pQCD calculation by  //
619 // Mangano, Nason and Ridolfi.                                       //
620 //                                                                   //
621 // For details and for the NORMALIZATION of the yields see:          //
622 //   N.Carrer and A.Dainese,                                         //
623 //   "Charm and beauty production at the LHC",                       //
624 //   ALICE-INT-2003-019, [arXiv:hep-ph/0311225];                     //
625 //   PPR Chapter 6.6, CERN/LHCC 2005-030 (2005).                     //
626 //*******************************************************************//
627   AliGenPythia * gener = 0x0;
628
629   switch(proc) {
630   case kCharmPbPb5500:
631       comment = comment.Append(" Charm in Pb-Pb at 5.5 TeV");
632       gener = new AliGenPythia(nEvts);
633       gener->SetProcess(kPyCharmPbPbMNR);
634       gener->SetStrucFunc(kCTEQ4L);
635       gener->SetPtHard(2.1,-1.0);
636       gener->SetEnergyCMS(5500.);
637       gener->SetNuclei(208,208);
638       break;
639   case kCharmpPb8800:
640       comment = comment.Append(" Charm in p-Pb at 8.8 TeV");
641       gener = new AliGenPythia(nEvts);
642       gener->SetProcess(kPyCharmpPbMNR);
643       gener->SetStrucFunc(kCTEQ4L);
644       gener->SetPtHard(2.1,-1.0);
645       gener->SetEnergyCMS(8800.);
646       gener->SetProjectile("P",1,1);
647       gener->SetTarget("Pb",208,82);
648       break;
649   case kCharmpp14000:
650       comment = comment.Append(" Charm in pp at 14 TeV");
651       gener = new AliGenPythia(nEvts);
652       gener->SetProcess(kPyCharmppMNR);
653       gener->SetStrucFunc(kCTEQ4L);
654       gener->SetPtHard(2.1,-1.0);
655       gener->SetEnergyCMS(14000.);
656       break;
657   case kCharmpp14000wmi:
658       comment = comment.Append(" Charm in pp at 14 TeV with mult. interactions");
659       gener = new AliGenPythia(-1);
660       gener->SetProcess(kPyCharmppMNRwmi);
661       gener->SetStrucFunc(kCTEQ5L);
662       gener->SetPtHard(ptHardMin,ptHardMax);
663       gener->SetEnergyCMS(14000.);
664       break;
665   case kD0PbPb5500:
666       comment = comment.Append(" D0 in Pb-Pb at 5.5 TeV");
667       gener = new AliGenPythia(nEvts);
668       gener->SetProcess(kPyD0PbPbMNR);
669       gener->SetStrucFunc(kCTEQ4L);
670       gener->SetPtHard(2.1,-1.0);
671       gener->SetEnergyCMS(5500.);
672       gener->SetNuclei(208,208);
673       break;
674   case kD0pPb8800:
675       comment = comment.Append(" D0 in p-Pb at 8.8 TeV");
676       gener = new AliGenPythia(nEvts);
677       gener->SetProcess(kPyD0pPbMNR);
678       gener->SetStrucFunc(kCTEQ4L);
679       gener->SetPtHard(2.1,-1.0);
680       gener->SetEnergyCMS(8800.);
681       gener->SetProjectile("P",1,1);
682       gener->SetTarget("Pb",208,82);
683       break;
684   case kD0pp14000:
685       comment = comment.Append(" D0 in pp at 14 TeV");
686       gener = new AliGenPythia(nEvts);
687       gener->SetProcess(kPyD0ppMNR);
688       gener->SetStrucFunc(kCTEQ4L);
689       gener->SetPtHard(2.1,-1.0);
690       gener->SetEnergyCMS(14000.);
691       break;
692   case kDPlusPbPb5500:
693       comment = comment.Append(" DPlus in Pb-Pb at 5.5 TeV");
694       gener = new AliGenPythia(nEvts);
695       gener->SetProcess(kPyDPlusPbPbMNR);
696       gener->SetStrucFunc(kCTEQ4L);
697       gener->SetPtHard(2.1,-1.0);
698       gener->SetEnergyCMS(5500.);
699       gener->SetNuclei(208,208);
700       break;
701   case kDPluspPb8800:
702       comment = comment.Append(" DPlus in p-Pb at 8.8 TeV");
703       gener = new AliGenPythia(nEvts);
704       gener->SetProcess(kPyDPluspPbMNR);
705       gener->SetStrucFunc(kCTEQ4L);
706       gener->SetPtHard(2.1,-1.0);
707       gener->SetEnergyCMS(8800.);
708       gener->SetProjectile("P",1,1);
709       gener->SetTarget("Pb",208,82);
710       break;
711   case kDPluspp14000:
712       comment = comment.Append(" DPlus in pp at 14 TeV");
713       gener = new AliGenPythia(nEvts);
714       gener->SetProcess(kPyDPlusppMNR);
715       gener->SetStrucFunc(kCTEQ4L);
716       gener->SetPtHard(2.1,-1.0);
717       gener->SetEnergyCMS(14000.);
718       break;
719   case kBeautyPbPb5500:
720       comment = comment.Append(" Beauty in Pb-Pb at 5.5 TeV");
721       gener = new AliGenPythia(nEvts);
722       gener->SetProcess(kPyBeautyPbPbMNR);
723       gener->SetStrucFunc(kCTEQ4L);
724       gener->SetPtHard(2.75,-1.0);
725       gener->SetEnergyCMS(5500.);
726       gener->SetNuclei(208,208);
727       break;
728   case kBeautypPb8800:
729       comment = comment.Append(" Beauty in p-Pb at 8.8 TeV");
730       gener = new AliGenPythia(nEvts);
731       gener->SetProcess(kPyBeautypPbMNR);
732       gener->SetStrucFunc(kCTEQ4L);
733       gener->SetPtHard(2.75,-1.0);
734       gener->SetEnergyCMS(8800.);
735       gener->SetProjectile("P",1,1);
736       gener->SetTarget("Pb",208,82);
737       break;
738   case kBeautypp14000:
739       comment = comment.Append(" Beauty in pp at 14 TeV");
740       gener = new AliGenPythia(nEvts);
741       gener->SetProcess(kPyBeautyppMNR);
742       gener->SetStrucFunc(kCTEQ4L);
743       gener->SetPtHard(2.75,-1.0);
744       gener->SetEnergyCMS(14000.);
745       break;
746   case kBeautypp14000wmi:
747       comment = comment.Append(" Beauty in pp at 14 TeV with mult. interactions");
748       gener = new AliGenPythia(-1);
749       gener->SetProcess(kPyBeautyppMNRwmi);
750       gener->SetStrucFunc(kCTEQ5L);
751       gener->SetPtHard(ptHardMin,ptHardMax);
752       gener->SetEnergyCMS(14000.);
753       break;
754    }
755
756   return gener;
757 }
758
759 AliGenPythia *PythiaHard(PDC07Proc_t proc) {
760 //*******************************************************************//
761 // Configuration file for hard QCD processes generation with PYTHIA  //
762 //                                                                   //
763 //*******************************************************************//
764   AliGenPythia * gener = 0x0;
765
766   switch(proc) {
767
768   case kPyJetJet:
769     comment = comment.Append(" pp->jet + jet over at 14 TeV, no restriction");
770     AliGenPythia * gener = new AliGenPythia(nEvts);
771     gener->SetEnergyCMS(eCMS);//        Centre of mass energy
772     gener->SetProcess(kPyJets);//        Process type
773     gener->SetJetEtaRange(-1.5, 1.5);//  Final state kinematic cuts
774     gener->SetJetPhiRange(0., 360.);
775     gener->SetJetEtRange(10., 1000.);
776     gener->SetPtHard(ptHardMin, ptHardMax);// Pt transfer of the hard scattering
777     gener->SetStrucFunc(kCTEQ4L);
778     return gener;
779   case kBeautyppwmiJet:
780         comment = comment.Append(" Beauty Jets over EMCAL");
781         AliGenPythia *gener = new AliGenPythia(-1);
782         gener->SetProcess(kPyBeautyppMNRwmi);
783         gener->SetStrucFunc(kCTEQ4L);
784         //Jet specific stuff below
785         gener->SetPtHard(ptHardMin,ptHardMax);
786         gener->SetJetEtaRange(-1,1);
787         gener->SetJetPhiRange(60.,210.);
788         gener->SetJetEtRange(10.,1000.);
789         //jet specific stuff above
790         gener->SetEnergyCMS(eCMS);
791         gener->SetForceDecay(kSemiElectronic);
792         gener->SetYRange(-2,2);
793         gener->SetFeedDownHigherFamily(kFALSE);
794         gener->SetCountMode(AliGenPythia::kCountParents);               
795       break;
796   case kPyGammaJetPHOS:
797       comment = comment.Append(" pp->jet + gamma over PHOS");
798       gener = new AliGenPythia(nEvts);
799       gener->SetEnergyCMS(eCMS);
800       gener->SetProcess(kPyDirectGamma);
801       gener->SetStrucFunc(kCTEQ4L);
802       gener->SetPtHard(ptHardMin,ptHardMax);
803       //gener->SetYHard(-1.0,1.0);
804       gener->SetGammaEtaRange(-0.13,0.13);
805       gener->SetGammaPhiRange(219.,321.);//Over 5 modules +-2 deg
806       break;
807   case kPyJetJetPHOS:
808       comment = comment.Append(" pp->jet + jet over PHOS");
809       gener = new AliGenPythia(nEvts);
810       gener->SetEnergyCMS(eCMS);
811       gener->SetProcess(kPyJets);
812       gener->SetStrucFunc(kCTEQ4L);
813       gener->SetPtHard(ptHardMin,ptHardMax);
814       //gener->SetYHard(-1.0,1.0);
815       gener->SetJetEtaRange(-1.,1.);
816       gener->SetJetPhiRange(200.,340.);
817       gener->SetJetEtRange(10., 1000.);
818       gener->SetPi0InPHOS(pi0FragPhotonInDet);
819       gener->SetFragPhotonOrPi0MinPt(ptGammaPi0Min);
820
821       printf("\n \n Event generator: Minimum pT of particle in calorimeter %f \n \n", ptGammaPi0Min);
822       break;
823   case kPyGammaBremsPHOS:
824       comment = comment.Append(" pp->jet + jet+bremsphoton over PHOS at 14 TeV");
825       gener = new AliGenPythia(nEvts);
826       gener->SetEnergyCMS(eCMS);
827       gener->SetProcess(kPyJets);
828       gener->SetStrucFunc(kCTEQ4L);
829       gener->SetPtHard(ptHardMin,ptHardMax);
830       //gener->SetYHard(-1.0,1.0);
831       gener->SetJetEtaRange(-1.,1.);
832       gener->SetJetPhiRange(200.,340.);
833       gener->SetFragPhotonInPHOS(pi0FragPhotonInDet);
834       gener->SetFragPhotonOrPi0MinPt(ptGammaPi0Min);
835           printf("\n \n Event generator: Minimum pT of particle in calorimeter %f \n \n", ptGammaPi0Min);
836           break;
837    case kPyJetJetPHOSv2:
838         comment = comment.Append(" pp->jet + jet over PHOS version2 ");
839         gener = new AliGenPythia(nEvts);
840         gener->SetEnergyCMS(eCMS);
841         gener->SetProcess(kPyJets);
842         gener->SetStrucFunc(kCTEQ4L);
843         gener->SetPtHard(ptHardMin,ptHardMax);
844         //gener->SetYHard(-1.0,1.0);
845         gener->SetJetEtaRange(-1.,1.);
846         gener->SetJetPhiRange(200.,340.);
847         //gener->SetPi0InPHOS(kTRUE);
848         gener->SetPhotonInPHOSeta(pi0FragPhotonInDet);
849         gener->SetPhotonMinPt(ptGammaPi0Min);
850         gener->SetForceDecay(kAll);
851         break;
852   case kPyGammaJetEMCAL:
853     comment = comment.Append(" pp->jet + gamma over EMCAL at 14 TeV");
854     gener = new AliGenPythia(nEvts);
855     gener->SetEnergyCMS(eCMS);
856     gener->SetProcess(kPyDirectGamma);
857     gener->SetStrucFunc(kCTEQ4L);
858     gener->SetPtHard(ptHardMin,ptHardMax);
859     //gener->SetYHard(-1.0,1.0);
860     gener->SetGammaEtaRange(-0.701,0.701);
861     gener->SetGammaPhiRange(79.,191.);//Over 6 supermodules +-2 deg
862     break;
863   case kPyJetJetEMCAL:
864     comment = comment.Append(" pp->jet + jet over EMCAL at 14 TeV");
865     gener = new AliGenPythia(nEvts);
866     gener->SetEnergyCMS(eCMS);
867     gener->SetProcess(kPyJets);
868     gener->SetStrucFunc(kCTEQ4L);
869     gener->SetPtHard(ptHardMin,ptHardMax);
870     //gener->SetYHard(-1.0,1.0);
871     gener->SetJetEtaRange(-1,1);
872     gener->SetJetPhiRange(60.,210.);
873         gener->SetJetEtRange(10., 1000.);
874     gener->SetPi0InEMCAL(pi0FragPhotonInDet);
875     gener->SetFragPhotonOrPi0MinPt(ptGammaPi0Min);
876     printf("\n \n Event generator: Minimum pT of particle in calorimeter %f \n \n", ptGammaPi0Min);
877     break;
878   case kPyGammaBremsEMCAL:
879       comment = comment.Append(" pp->jet + jet+bremsphoton over EMCAL at 14 TeV");
880       gener = new AliGenPythia(nEvts);
881       gener->SetEnergyCMS(eCMS);
882       gener->SetProcess(kPyJets);
883       gener->SetStrucFunc(kCTEQ4L);
884       gener->SetPtHard(ptHardMin,ptHardMax);
885       //gener->SetYHard(-1.0,1.0);
886       gener->SetJetEtaRange(-1,1);
887       gener->SetJetPhiRange(60.,210.); 
888           gener->SetJetEtRange(10., 1000.);
889       gener->SetFragPhotonInEMCAL(pi0FragPhotonInDet);
890       gener->SetFragPhotonOrPi0MinPt(ptGammaPi0Min);
891
892       printf("\n \n Event generator: Minimum pT of particle in calorimeter %f \n \n", ptGammaPi0Min);
893       break;
894
895   }
896
897   return gener;
898 }
899
900 AliGenerator* MbCocktail()
901 {
902       comment = comment.Append(" pp at 14 TeV: Pythia low-pt, no heavy quarks + J/Psi from parameterisation");
903       AliGenCocktail * gener = new AliGenCocktail();
904       gener->UsePerEventRates();
905
906 //    Pythia
907
908       AliGenPythia* pythia = new AliGenPythia(-1);
909       pythia->SetMomentumRange(0, 999999.);
910       pythia->SetThetaRange(0., 180.);
911       pythia->SetYRange(-12.,12.);
912       pythia->SetPtRange(0,1000.);
913       pythia->SetProcess(kPyMb);
914       pythia->SetEnergyCMS(14000.);
915       pythia->SwitchHFOff();
916
917 //   J/Psi parameterisation
918
919       AliGenParam* jpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF scaled", "Jpsi");
920       jpsi->SetPtRange(0.,100.);
921       jpsi->SetYRange(-8., 8.);
922       jpsi->SetPhiRange(0., 360.);
923       jpsi->SetForceDecay(kAll);
924
925       gener->AddGenerator(pythia, "Pythia", 1.);
926       gener->AddGenerator(jpsi,   "J/Psi", 8.e-4);
927
928       return gener;
929 }
930
931 AliGenerator* PyMbTriggered(Int_t pdg)
932 {
933     AliGenPythia* pythia = new AliGenPythia(-1);
934     pythia->SetMomentumRange(0, 999999.);
935     pythia->SetThetaRange(0., 180.);
936     pythia->SetYRange(-12.,12.);
937     pythia->SetPtRange(0,1000.);
938     pythia->SetProcess(kPyMb);
939     pythia->SetEnergyCMS(14000.);
940     pythia->SetTriggerParticle(pdg, 0.9);
941     return pythia;
942 }
943
944 void ProcessEnvironmentVars()
945 {
946         cout << "######################################" << endl;
947         cout << "## Processing environment variables ##" << endl;
948         cout << "######################################" << endl;
949         
950     // Run Number
951     if (gSystem->Getenv("RUN")) {
952                 runNumber = atoi(gSystem->Getenv("RUN"));
953     }
954     //cout<<"Run number "<<runNumber<<endl;
955         
956     // Event Number
957     if (gSystem->Getenv("EVENT")) {
958                 eventNumber = atoi(gSystem->Getenv("EVENT"));
959     }
960     //cout<<"Event number "<<eventNumber<<endl;
961         
962         // Random Number seed
963     if (gSystem->Getenv("CONFIG_SEED")) {
964                 seed = atoi(gSystem->Getenv("CONFIG_SEED"));
965     }
966     else if(gSystem->Getenv("EVENT") && gSystem->Getenv("RUN")){
967                 seed = runNumber * 100000 + eventNumber;
968     }
969         
970     gRandom->SetSeed(seed);
971         
972     cout<<"////////////////////////////////////////////////////////////////////////////////////"<<endl;
973     cout<<"Seed for random number generation= "<< seed<<"  "<< gRandom->GetSeed()<<endl;
974     cout<<"////////////////////////////////////////////////////////////////////////////////////"<<endl;
975         
976     // Run type
977     if (gSystem->Getenv("DC_RUN_TYPE")) {
978                 cout<<"run type "<<gSystem->Getenv("DC_RUN_TYPE")<<endl;
979                 for (Int_t iRun = 0; iRun < kRunMax; iRun++) {
980                         if (strcmp(gSystem->Getenv("DC_RUN_TYPE"), pprRunName[iRun])==0) {
981                                 proc = (PDC07Proc_t)iRun;
982                         }
983                 }
984     }
985     else
986                 cout << "Environment variable DC_RUN_TYPE is not defined" << endl;
987     if (gSystem->Getenv("YEAR"))
988                 year = atoi(gSystem->Getenv("YEAR"));
989         if (gSystem->Getenv("ECMS"))
990                 eCMS = atof(gSystem->Getenv("ECMS"));
991         if (gSystem->Getenv("TRIGGER"))
992                 trig = atoi(gSystem->Getenv("TRIGGER"));
993         if ( gSystem->Getenv("PI0GAMMAINDET") )
994                 pi0FragPhotonInDet      = atoi(gSystem->Getenv("PI0GAMMAINDET"));
995     if (gSystem->Getenv("PTGAMMAPI0MIN"))
996                 ptGammaPi0Min = atof(gSystem->Getenv("PTGAMMAPI0MIN"));
997         if (gSystem->Getenv("QUENCHING"))
998                 iquenching = atof(gSystem->Getenv("QUENCHING"));
999         if (gSystem->Getenv("QHAT"))
1000                 qhat = atof(gSystem->Getenv("QHAT"));
1001         //if (gSystem->Getenv("MEDLENGHT"))
1002         //      medlength = atof(gSystem->Getenv("MEDLENGHT"));
1003         if (gSystem->Getenv("PTHARDMIN"))
1004                 ptHardMin = atof(gSystem->Getenv("PTHARDMIN"));
1005         if (gSystem->Getenv("PTHARDMAX"))
1006                 ptHardMax = atof(gSystem->Getenv("PTHARDMAX"));
1007                                         
1008     if(gSystem->Getenv("PTHARDMIN") && gSystem->Getenv("PTHARDMAX")){
1009         if( strncmp(pprRunName[proc],"kPyJetJet",9)==0 || strncmp(pprRunName[proc],"kBeautyppwmiJet",15)==0 ){
1010                         ptHardMin = ptHardJJLo[runNumber%nPtHardJJBins];
1011                         ptHardMax = ptHardJJHi[runNumber%nPtHardJJBins];
1012                 }
1013                 else if( strncmp(pprRunName[proc],"kPyGammaJet",11)==0){
1014                         ptHardMin = ptHardGJLo[runNumber%nPtHardGJBins];
1015                         ptHardMax = ptHardGJHi[runNumber%nPtHardGJBins];        
1016                 }       
1017         }
1018         
1019     cout<<">> Run type set to "<<pprRunName[proc]<<endl;
1020     cout<<">> Collision energy set to "<<eCMS <<endl;
1021         cout<<">> Collisions trigger type "<<trig<<" "<<TrigConfName[trig]<<endl;
1022     cout<<">> ptHard limits: "<<ptHardMin<<" to " <<ptHardMax<<" GeV"<<endl;
1023         if(pi0FragPhotonInDet) cout<<">> pt gamma/pi0 threshold "<< ptGammaPi0Min<<" GeV "<<endl;
1024     if(iquenching) cout<<">> quenching on? "<< iquenching<<" qhat "<<qhat//<<" medlength "<<medlength
1025         <<endl;
1026         cout<<">> Year geometry : "<<year<<endl;        
1027     cout << "######################################" << endl;
1028     cout << "######################################" << endl;
1029 }