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