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