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