d3e38bf37043548fa8a0100f7e0630d596d82c2e
[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / sim / Config.C
1 // -------------------------------------------------------------------
2 struct Setup
3 {
4   TString runType;    // Event generator chosen
5   UInt_t  seed;       // Random number seed - (env)
6   Float_t minB;       // Least imp.param. (env)
7   Float_t maxB;       // Largest imp.param. (env)  
8   Int_t   hftype;     // Heavy flavour type (random)
9   /** 
10    * Get the value of an environment variable as a float 
11    * 
12    * @param envName Enviroment variable name 
13    * @param def     Default value 
14    *
15    * @return As float, or default
16    */
17   static Float_t Env2Float(const char* envName, Float_t def) 
18   { 
19     TString val(gSystem->Getenv(envName));
20     if (val.IsNull()) return def;
21     return val.Atof();
22   }
23   /** 
24    * Get the value of an environment variable as a unsigned int 
25    * 
26    * @param envName Enviroment variable name 
27    * @param def     Default value 
28    *
29    * @return As unsigned int, or default
30    */
31   static UInt_t Env2UInt(const char* envName, UInt_t def)
32   {
33     TString val(gSystem->Getenv(envName));
34     if (val.IsNull()) return def;
35     return UInt_t(val.Atoll());
36   }
37   /** 
38    * Get the value of an environment variable as a int
39    * 
40    * @param envName Enviroment variable name 
41    * @param def     Default value 
42    *
43    * @return As int, or default
44    */
45   static UInt_t Env2Int(const char* envName, Int_t def)
46   {
47     TString val(gSystem->Getenv(envName));
48     if (val.IsNull()) return def;
49     return val.Atoi();
50   }
51
52   /** 
53    * Constructor - retrieves needed information from CDB manager and
54    * environment.
55    */
56   Setup() 
57     : runType(""),
58       seed(0),
59       minB(0), 
60       maxB(100),
61       hftype(-1)
62   {
63     TDatime now;
64     runType = gSystem->Getenv("CONFIG_RUN_TYPE");
65     seed    = Env2UInt("CONFIG_SEED", now.Get());
66     minB    = Env2Float("CONFIG_BMIN", 0);
67     maxB    = Env2Float("CONFIG_BMAX", 100);
68     if (runType[0] == 'k') runType.Remove(0,1);
69     runType.ToLower();
70
71     // gROOT->Macro("GetGRP.C");
72
73     Int_t pytR = gRandom->Rndm();
74     if      (pytR < 0.16) hftype = 0; 
75     else if (pytR < 0.32) hftype = 1; 
76     else if (pytR < 0.48) hftype = 2; 
77     else if (pytR < 0.64) hftype = 3; 
78     else if (pytR < 0.72) hftype = 4; 
79     else if (pytR < 0.80) hftype = 5; 
80     else                  hftype = 6; 
81
82     if (runType.IsNull() || runType == "default") DeduceRunType();
83
84     Print();
85   }
86   void Print()
87   {
88     Printf("=======================================================\n"
89            "           Set-up of the simulation\n");
90     grp->Print();
91     Printf("Run type: '%s'", runType.Data());
92     Printf("Seed:     %d", seed);
93     Printf("\n"
94            "=======================================================");
95   }
96   /** 
97    * Set the default generator based on the beam type 
98    *
99    * - p-p PYTHIA
100    * - p-A or A-p DPMJet
101    * - A-A Hijing 
102    */
103   void DeduceRunType()
104   {
105     if      (grp->IsPP())                runType = "pythia";
106     else if (grp->IsPA() || grp->IsAP()) runType = "dpmjet";
107     else if (grp->IsAA())                runType = "hijing";
108   }
109   void LoadGen() {
110     if (!gROOT->GetClass("AliStructFuncType")) 
111       gSystem->Load("liblhapdf");      // Parton density functions
112     if (!gROOT->GetClass("TPythia6"))
113       gSystem->Load("libEGPythia6")   // TGenerator interface
114     if (!runType.EqualTo("hydjet")) 
115       LoadPythia(false);
116   }
117   /** 
118    * Load the pythia libraries 
119    * 
120    * @param vers Optional version post-fix
121    */
122   void LoadPythia(Bool_t gen=true, const char* vers="6.4.21")
123   {
124     if (gen) LoadGen();
125     char m = vers[0];
126     if (gROOT->GetClass(Form("AliPythia6%c", m))) return;
127     gSystem->Load(Form("libpythia%s", vers));
128     gSystem->Load(Form("libAliPythia%c", m));
129   }
130   /** 
131    * Load HIJING libraries 
132    */
133   void LoadHijing() 
134   {
135     LoadPythia();
136     if (gROOT->GetClass("THijing")) return;
137     gSystem->Load("libhijing");
138     gSystem->Load("libTHijing");
139     AliPDG::AddParticlesToPdgDataBase();
140   }
141   /** 
142    * Load HydJet libraries 
143    */
144   void LoadHydjet()
145   {
146     if (gROOT->GetClass("TUHKMgen")) return;
147     gSystem->Load("libTUHKMgen");
148   }
149   /** 
150    * Load DPMJet libraries
151    */
152   void LoadDpmjet()
153   {
154     LoadPythia();
155     if (gROOT->GetClass("TDPMjet")) return;
156     gSystem->Load("libdpmjet");
157     gSystem->Load("libTDPMjet");
158   }
159   /** 
160    * Load AMPT libraries
161    */
162   void LoadAmpt() 
163   {
164     LoadPythia();
165     if (gROOT->GetClass("TAmpt")) return;
166     gSystem->Load("libampt");
167     gSystem->Load("libTAmpt");
168   }
169   /** 
170    * Make the generator 
171    * 
172    * @return Point to newly allocated generator or null 
173    */
174   AliGenerator* MakeGenerator()
175   {
176     Bool_t asym = grp->IsPA()||grp->IsAP();
177     TString& rt = runType;
178     AliGenerator* g = 0;
179     if      (rt.EndsWith("perugia0chadr"))     g=PythiaHF(0);
180     else if (rt.EndsWith("perugia0bchadr"))    g=PythiaHF(1);
181     else if (rt.EndsWith("perugia0cele"))      g=PythiaHF(2);
182     else if (rt.EndsWith("perugia0bele"))      g=PythiaHF(3);
183     else if (rt.EndsWith("perugia0jspi2e"))    g=PythiaHF(4);
184     else if (rt.EndsWith("perugia0btojspi2e")) g=PythiaHF(5);
185     else if (rt.BeginsWith("pythia"))          g=Pythia(rt);
186     else if (rt.BeginsWith("hijing2000hf"))    g=HFCocktail(rt);
187     else if (rt.BeginsWith("hijing2000"))      g=Hijing(asym, 
188                                                         false, 2.3);
189     else if (rt.BeginsWith("hijing"))          g=Hijing(asym, 
190                                                         grp->IsAA(), 0);
191     else if (rt.BeginsWith("ampthf"))          g=HFCocktail(rt);
192     else if (rt.BeginsWith("ampt"))            g=Ampt();
193     else if (rt.BeginsWith("dpmjet"))          g=Dpmjet();
194     else if (rt.BeginsWith("phojet"))          g=Dpmjet();
195     else if (rt.BeginsWith("hydjet"))          g=Hydjet();
196
197     if (g) g->SetVertexSmear(AliGenerator::kPerEvent);
198     else 
199       Fatal("", "Invalid run type \"%s\" specified", runType.Data());
200     return g;
201   }
202   TVirtualMCDecayer* MakeDecayer()
203   {
204     if (runType.BeginsWith("hydjet")) return 0;
205
206     LoadPythia();
207     TVirtualMCDecayer* decayer = new AliDecayerPythia();
208     if (runType.EqualTo("hijing2000hf") && hftype < 2) 
209       decayer->SetForceDecay(kHadronicD);
210     else 
211       decayer->SetForceDecay(kAll);
212     decayer->Init();
213     return decayer;
214   }
215
216   // === PYTHIA ========================================================
217   // Normal 
218   AliGenerator* Pythia(const TString & tune)
219   {
220     // Int_t kCTEQ6l = 8;
221     if (!grp->IsPP()) Fatal("Setup", "Pythia6 only works for pp");
222
223     TString t(tune);
224     t.ToUpper();
225     t.ReplaceAll("PYTHIA6", "");
226     t.ReplaceAll("PYTHIA", "");
227     Info("Setup", "Making Pythia6 event generator (tune: %s)", t.Data());
228
229     LoadPythia();
230     AliGenPythia* pythia = new AliGenPythia(-1); 
231     pythia->SetMomentumRange(0, 999999.);
232     pythia->SetThetaRange(0., 180.);
233     pythia->SetYRange(-12.,12.);
234     pythia->SetPtRange(0,1000.);
235     pythia->SetProcess(kPyMb);
236     pythia->SetEnergyCMS(grp->energy);
237
238     if (t == "D6T") {
239       //    Tune
240       //    109     D6T : Rick Field's CDF Tune D6T 
241       //                  (NB: needs CTEQ6L pdfs externally)
242       pythia->SetTune(109); // F I X 
243       pythia->SetStrucFunc(kCTEQ6l);
244     }
245     else if (t == "PERUGIA0") { 
246       //    Tune
247       //    320     Perugia 0
248       pythia->SetTune(320); 
249       pythia->UseNewMultipleInteractionsScenario();
250     }
251     else if (t == "ATLAS") {
252       //    Tune
253       //    C   306 ATLAS-CSC: Arthur Moraes' (new) ATLAS tune 
254       //                        (needs CTEQ6L externally)
255       pythia->SetTune(306);
256       pythia->SetStrucFunc(kCTEQ6l);
257     }
258     else if (t == "JETS") { 
259       pythia->SetProcess(kPyJets);
260       pythia->SetStrucFunc(kCTEQ6l);
261       pythia->SetJetEtaRange(-1.5, 1.5); 
262       pythia->SetJetEtRange(50., 800.);
263       pythia->SetPtHard(45., 1000.);
264       pythia->SetPycellParameters(2.2, 300, 432, 0., 4., 5., 0.7);
265     }
266     else if (t == "ATLAS_FLAT") {
267       // set high multiplicity trigger
268       // this weight achieves a flat multiplicity distribution
269       TH1 *weight = new TH1D("weight","weight",201,-0.5,200.5);
270       weight->SetBinContent(1,5.49443);
271       weight->SetBinContent(2,8.770816);
272       weight->SetBinContent(6,0.4568624);
273       weight->SetBinContent(7,0.2919915);
274       weight->SetBinContent(8,0.6674189);
275       weight->SetBinContent(9,0.364737);
276       weight->SetBinContent(10,0.8818444);
277       weight->SetBinContent(11,0.531885);
278       weight->SetBinContent(12,1.035197);
279       weight->SetBinContent(13,0.9394057);
280       weight->SetBinContent(14,0.9643193);
281       weight->SetBinContent(15,0.94543);
282       weight->SetBinContent(16,0.9426507);
283       weight->SetBinContent(17,0.9423649);
284       weight->SetBinContent(18,0.789456);
285       weight->SetBinContent(19,1.149026);
286       weight->SetBinContent(20,1.100491);
287       weight->SetBinContent(21,0.6350525);
288       weight->SetBinContent(22,1.351941);
289       weight->SetBinContent(23,0.03233504);
290       weight->SetBinContent(24,0.9574557);
291       weight->SetBinContent(25,0.868133);
292       weight->SetBinContent(26,1.030998);
293       weight->SetBinContent(27,1.08897);
294       weight->SetBinContent(28,1.251382);
295       weight->SetBinContent(29,0.1391099);
296       weight->SetBinContent(30,1.192876);
297       weight->SetBinContent(31,0.448944);
298       for (Int_t i = 32; i <= 201; i++) weight->SetBinContent(i,1);
299       weight->SetEntries(526);
300       
301       Int_t limit = weight->GetRandom();
302       pythia->SetTriggerChargedMultiplicity(limit, 1.4);
303     }
304     return pythia;
305   }
306   AliGenerator* PythiaHF(Int_t type, Bool_t harder=0) 
307   { 
308     LoadPythia();
309     if (type == 6) return Pythia("jets");
310     if (type == 4) { 
311       AliGenParam *jpsi =  AliGenParam(1, AliGenMUONlib::kJpsi,
312                                        (harder?"CDF pp 8.8":"CDF pp 7"),"Jpsi");
313       jpsi->SetPtRange(0.,999.);
314       jpsi->SetYRange(-1.0, 1.0);
315       jpsi->SetPhiRange(0.,360.);
316       jpsi->SetForceDecay(kDiElectron);
317       return jpsi;
318     }
319     AliGenPythia* pythia = static_cast<AliGenPythia*>(Pythia("PERUGIA0"));
320     switch (type) { 
321     case 0: // chadr
322       pythia->SetProcess(kPyCharmppMNRwmi);
323       pythia->SetForceDecay(kHadronicD);
324       break;
325     case 1: // bchadr
326       pythia->SetProcess(kPyBeautyppMNRwmi);
327       pythia->SetForceDecay(kHadronicD);
328       break;
329     case 2: // cele
330       pythia->SetProcess(kPyCharmppMNRwmi);
331       pythia->SetCutOnChild(1);
332       pythia->SetPdgCodeParticleforAcceptanceCut(11);
333       pythia->SetChildYRange(-1.2,1.2);
334       pythia->SetChildPtRange(0,10000.);
335       break;
336     case 3: // bele
337       pythia->SetProcess(kPyBeautyppMNRwmi);
338       pythia->SetCutOnChild(1);
339       pythia->SetPdgCodeParticleforAcceptanceCut(11);
340       pythia->SetChildYRange(-1.2,1.2);
341       pythia->SetChildPtRange(0,10000.);
342       break;
343     case 5:
344       pythia->SetProcess(kPyBeautyppMNRwmi);
345       pythia->SetCutOnChild(1);
346       pythia->SetPdgCodeParticleforAcceptanceCut(443);
347       pythia->SetChildYRange(-2,2);
348       pythia->SetChildPtRange(0,10000.);
349     }
350     return pythia;
351   }
352   /** 
353    * Make a Min-Bias AA, pA, or Ap Hijing generator 
354    * 
355    * @param slowN If true, make a cocktail with slow neutrons 
356    * 
357    * @return Generator 
358    */
359   AliGenerator* Hijing(Bool_t slowN=false, Bool_t quench=1, Float_t ptHard=0) 
360   {
361     LoadHijing();
362     AliGenHijing *gener = new AliGenHijing(-1);
363     // centre of mass energy 
364     gener->SetEnergyCMS(grp->energy);
365     gener->SetImpactParameterRange(minB, maxB); 
366     // reference frame
367     gener->SetReferenceFrame("CMS");
368     // projectil
369     gener->SetProjectile(grp->beam1.Name(), grp->beam1.a, grp->beam1.z);
370     gener->SetTarget    (grp->beam2.Name(), grp->beam2.a, grp->beam2.z);
371     // tell hijing to keep the full parent child chain
372     gener->KeepFullEvent();
373     // enable jet quenching
374     gener->SetJetQuenching(quench);
375     // enable shadowing
376     gener->SetShadowing(slowN);
377     // Don't track spectators
378     gener->SetSpectators(!slowN);
379     // 
380     if (ptHard > 0) hi->SetPtHardMin(ptHard);
381
382     // kinematic selection
383     gener->SetSelectAll(0);
384     // Boosted CMS 
385     gener->SetBoostLHC(grp->IsPA() || grp->IsAP());
386     // No need for cocktail 
387     if (!slowN || !grp->IsPA() || !grp->IsAP()) return gener;
388
389     AliGenCocktail* cocktail = new AliGenCocktail();
390     cocktail->SetProjectile(grp->beam1.Name(), grp->beam1.a, grp->beam1.z);
391     cocktail->SetTarget    (grp->beam2.Name(), grp->beam2.a, grp->beam2.z);
392     cocktail->SetEnergyCMS(grp->energy);
393
394     AliGenSlowNucleons*     gray  = new AliGenSlowNucleons(1);
395     AliCollisionGeometry*   coll  = gener->CollisionGeometry();
396     AliSlowNucleonModelExp* model = new AliSlowNucleonModelExp();
397     //  Not yet in the release...
398     //      model->SetSaturation(kTRUE);
399     gray->SetSlowNucleonModel(model);
400     gray->SetTarget(grp->beam2.a, grp->beam2.z);
401     gray->SetThetaDist(1);
402     gray->SetProtonDirection(grp->beam1.IsP() ? 1 : 2);
403     //      gray->SetDebug(1);
404     gray->SetNominalCmsEnergy(2*grp->beamEnergy);
405     gray->NeedsCollisionGeometry();
406     gray->SetCollisionGeometry(coll);
407
408     cocktail->AddGenerator(gener, "Hijing pPb", 1);
409     cocktail->AddGenerator(gray, "Gray Particles", 1);
410     
411     return cocktail;
412   }
413   /** 
414    * Make a DPMJet generator for AA, pA, or Ap. 
415    * 
416    * @param fragments If true, make fragments 
417    * 
418    * @return Generator 
419    */
420   AliGenerator* Dpmjet(Bool_t fragments=0)
421   {
422     LoadDpmjet();
423     AliGenDPMjet* dpmjet = new AliGenDPMjet(-1); 
424     dpmjet->SetEnergyCMS(grp->energy);
425     dpmjet->SetProjectile(grp->beam1.Name(), grp->beam1.a, grp->beam1.z);
426     dpmjet->SetTarget    (grp->beam2.Name(), grp->beam2.a, grp->beam2.z);
427     dpmjet->SetImpactParameterRange(minB, maxB);
428     dpmjet->SetProjectileBeamEnergy(grp->beam1.z*grp->beamEnergy/grp->beam1.a);
429     if (grp->IsAA()) { 
430       dpmjet->SetPi0Decay(0);
431     }
432     else if (grp->IsPA() || grp->IsAP()) { 
433       dpmjet->SetTriggerParticle(3312, 1.2, 2.0);
434       dpmjet->SetFragmentProd(fragments); // Alwas disabled 
435     }
436     else if (grp->IsPP()) { // PhoJet
437       dpmjet->SetMomentumRange(0, 999999.);
438       dpmjet->SetThetaRange(0., 180.);
439       dpmjet->SetYRange(-12.,12.);
440       dpmjet->SetPtRange(0,1000.);
441     }
442     return dpmjet;
443   }
444   /** 
445    * Make an AMPT generator for AA collisions 
446    * 
447    * @return Generator 
448    */
449   AliGenerator* Ampt()
450   {
451     LoadAmpt();
452     AliGenAmpt *genHi = new AliGenAmpt(-1);
453     genHi->SetEnergyCMS(grp->energy);
454     genHi->SetReferenceFrame("CMS");
455     genHi->SetProjectile(grp->beam1.Name(), grp->beam1.a, grp->beam1.z);
456     genHi->SetTarget    (grp->beam2.Name(), grp->beam2.a, grp->beam2.z);
457     genHi->SetPtHardMin (2);
458     genHi->SetImpactParameterRange(bMin,bMax);
459     // disable jet quenching
460     genHi->SetJetQuenching(0); 
461     // enable shadowing
462     genHi->SetShadowing(1);    
463     // neutral pion and heavy particle decays switched off
464     genHi->SetDecaysOff(1);
465     genHi->SetSpectators(0);   // track spectators 
466     genHi->KeepFullEvent();
467     genHi->SetSelectAll(0);
468     return genHi;
469   }
470   /** 
471    * Make an HydJet generator for A-A
472    * 
473    * @return Generator 
474    */
475   AliGenerator* Hydjet()
476   {
477     LoadHydjet();
478     AliGenUHKM *genHi = new AliGenUHKM(-1);
479     genHi->SetAllParametersLHC();
480     genHi->SetProjectile(grp->beam1.Name(), grp->beam1.a, grp->beam1.z);
481     genHi->SetTarget    (grp->beam2.Name(), grp->beam2.a, grp->beam2.z);
482     genHi->SetEcms(grp->energy);
483     genHi->SetEnergyCMS(grp->energy);
484     genHi->SetBmin(bMin);
485     genHi->SetBmax(bMax);
486     genHi->SetPyquenPtmin(9);
487     return genHi;
488   }
489   /** 
490    * Make a heavy flavour cocktail 
491    * 
492    * @param base Underlying event. 
493    * 
494    * @return Generator 
495    */
496   AliGeneator* HFCocktail(const TString& base) 
497   {
498     
499     AliGenCocktail *cocktail = new AliGenCocktail();
500     cocktail->SetProjectile(grp->beam1.Name(), grp->beam1.a, grp->beam1.z);
501     cocktail->SetTarget    (grp->beam2.Name(), grp->beam2.a, grp->beam2.z);
502     cocktail->SetEnergyCMS(grp->energy);
503     
504     // Add underlying event
505     if (base.BeginsWith("ampt", TString::kIgnoreCase)) { 
506       hi = Ampt();
507       cocktail->AddGenerator(hi,"ampt",1);
508     }
509     else { 
510       hi = Hijing(grp->IsPA() || grp->IsAP(), false, 2.3);
511       cocktail->AddGenerator(hi,"hijing",1);
512
513     }
514
515     // --- Default formula -------------------------------------------
516     TForumla* one = new TFormula("one", "1.");
517
518     // --- Pythia ----------------------------------------------------
519     AliGenerator* pythia = PythiaHF(hftype);
520     switch (hftype) { 
521     case 6: 
522       cocktail->AddGenerator(pythia, "pythiaJets", 1, one);
523       break;
524     defualt:
525       cocktail
526         ->AddGenerator(pythia, "pythiaHF", 1, 
527                        new TFormula("Signals", 
528                                     "20.*(x<5.)+80./3.*(1.-x/20.)*(x>5.)"));
529       break;
530     }
531     // --- Dummy -----------------------------------------------------
532     AliGenParam* param = 0;
533
534     // --- Phos stuff ------------------------------------------------
535     AliGenPHOSlib *plib = new AliGenPHOSlib();
536     Double_t lower[] = { 0, 3, 6, 9, 12, -1 };
537     Double_t *pLow   = lower;
538     for (Int_t i = 0; i < 5; i++) {
539       param = new AliGenParam(5, plib, AliGenPHOSlib::kPi0);
540       param->SetPhiRange(0., 360.) ;
541       param->SetYRange(-1.2, 1.2) ;
542       param->SetPtRange(lower[i], 30.) ;
543       cocktail->AddGenerator(param,Form("Pi0HagPt%d", i), 1., one);
544     }
545
546     // --- Jpsi->mu+ mu-----------------------------------------------
547     param = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp 3.94", "Jpsi");
548     param->SetPtRange(0.,999.);
549     param->SetYRange(-4.2, -2.3);
550     param->SetPhiRange(0.,360.);
551     param->SetForceDecay(kDiMuon);
552     param->SetCutOnChild(1);
553     param->SetChildPhiRange(0.,360.);
554     param->SetChildThetaRange(168.5,178.5);
555     cocktail->AddGenerator(param, "Jpsi2M", 1, one); 
556
557     // --- Chi_c -> J/Psi + gamma, J/Psi -> e+e- ---------------------
558     Float_t thmin  = (180./TMath::Pi())*2.*atan(exp(-1.2));  
559     Float_t thmax  = (180./TMath::Pi())*2.*atan(exp( 1.2));  
560     param  = new AliGenParam(1, AliGenMUONlib::kChic,"default","Chic");
561     param->SetMomentumRange(0, 999.);        // Wide cut on the momentum
562     param->SetPtRange(0, 100.);              // Wide cut on Pt
563     param->SetYRange(-2.5, 2.5);
564     param->SetCutOnChild(1);                 // Enable cuts on decay products
565     param->SetChildPhiRange(0., 360.);
566     // In the acceptance of the Central Barrel
567     param->SetChildThetaRange(thmin, thmax); 
568     // Chi_c -> J/Psi + gamma, J/Psi -> e+e-
569     param->SetForceDecay(kChiToJpsiGammaToElectronElectron); 
570     cocktail->AddGenerator(param, "Chi_c", 1, one); 
571
572     // --- Dummy -----------------------------------------------------
573     AliGenBox* box = 0;
574
575     // --- Some particles --------------------------------------------
576     Double_t    boxR   = gRandom->Integer(3)+1;
577     Int_t       boxP[] = { 310, 3122,   3312, 3322, 
578                            (boxR==1 ?3334: boxR==2 ?-3334: -3312) };
579     Int_t       boxN[] = {   1,    1,      3,    3, 1 }
580     const char* boxT[] = { "K0s", "Lambda", "Xi-", "Xi0", 
581                            (boxR==1 ? "Omega-": boxR==2 ? "Omega+": "Xi+") };
582     for (Int_t i = 0; i < 5; i++) {
583       box = new AliGenBox(boxN[i]);
584       box->SetPart(boxP[i]);
585       box->SetPtRange(0,13);
586       box->SetThetaRange(45, 135);
587       cocktail->AddGenerator(box, boxT[i], 1, one);
588     }
589     
590     // --- High pT charged particle ----------------------------------
591     TFormula* hptF = new TFormula("High Pt", 
592                                   "5.*(x<5.)+20./3.*(1.-x/20.)*(x > 5.)");
593     Int_t       hptP[] = { 211, 321, 2212 };
594     const char* hptT[] = { "pi", "K", "p" };
595     for (Int_t i = 0; i < 3; i++) { 
596       for (Int_t j = -1; j <= 1; j++) {
597         box->SetPart(j*hptP[i]);
598         box->SetPtRange(2., 50.);
599         box->SetThetaRange(thmin, thmax);
600         cocktail->AddGenerator(box, Form("%s%c",hptT[i],(j<0,'-','+')),1,hptF);
601       }
602     }
603     return cocktail;
604   }
605 };
606
607 void Config()
608 {
609   // --- Get settings from environment variables --------------------
610   Setup    s;
611
612   // ---- Seed random number generator -------------------------------
613   gRandom->SetSeed(s.seed);
614   std::cerr << "Seed for random number generation= " << s.seed << std::endl; 
615
616   //------------------------------------------------------------------
617   // 
618   // Geometry and tracking 
619   //
620   // --- Libraries required by geant321 ------------------------------
621   s.LoadGen();
622   gSystem->Load("libgeant321");
623   new TGeant3TGeo("C++ Interface to Geant3");
624
625   // -----------------------------------------------------------------
626   //  Create the output file
627   std::cout<< "Config.C: Creating Run Loader ..." << std::endl;
628   AliRunLoader* rl = AliRunLoader::Open("galice.root",
629                                         AliConfig::GetDefaultEventFolderName(),
630                                         "recreate");
631   if (!rl) Fatal("Config","Can not instatiate the Run Loader");
632
633   rl->SetCompressionLevel(2);
634   rl->SetNumberOfEventsPerFile(1000);
635   gAlice->SetRunLoader(rl);
636
637   // --- Trigger configuration ---------------------------------------
638   // AliSimulation::Instance()->SetTriggerConfig(grp->IsAA() ? "Pb-Pb" : "p-p");
639
640   //
641   //=======================================================================
642   // Steering parameters for ALICE simulation
643   // 
644   // --- Specify event type to be tracked through the ALICE setup
645   // --- All positions are in cm, angles in degrees, and P and E in GeV
646
647   // --- Process switches --------------------------------------------
648   gMC->SetProcess("DCAY",1);
649   gMC->SetProcess("PAIR",1);
650   gMC->SetProcess("COMP",1);
651   gMC->SetProcess("PHOT",1);
652   gMC->SetProcess("PFIS",0);
653   gMC->SetProcess("DRAY",0);
654   gMC->SetProcess("ANNI",1);
655   gMC->SetProcess("BREM",1);
656   gMC->SetProcess("MUNU",1);
657   gMC->SetProcess("CKOV",1);
658   gMC->SetProcess("HADR",1);
659   gMC->SetProcess("LOSS",2);
660   gMC->SetProcess("MULS",1);
661   gMC->SetProcess("RAYL",1);
662   
663   Float_t cut = 1.e-3;        // 1MeV cut by default
664   Float_t tofmax = 1.e10;
665
666   // --- Tracking cuts -----------------------------------------------
667   gMC->SetCut("CUTGAM", cut);
668   gMC->SetCut("CUTELE", cut);
669   gMC->SetCut("CUTNEU", cut);
670   gMC->SetCut("CUTHAD", cut);
671   gMC->SetCut("CUTMUO", cut);
672   gMC->SetCut("BCUTE",  cut); 
673   gMC->SetCut("BCUTM",  cut); 
674   gMC->SetCut("DCUTE",  cut); 
675   gMC->SetCut("DCUTM",  cut); 
676   gMC->SetCut("PPCUTM", cut);
677   gMC->SetCut("TOFMAX", tofmax); 
678   
679   // --- Set External decayer ----------------------------------------
680   TVirtualMCDecayer* decayer = s.MakeDecayer();
681   if (decayer) gMC->SetExternalDecayer(decayer);  
682
683   //------------------------------------------------------------------
684   // 
685   // Generator Configuration 
686   //
687   // --- Make the generator - this loads libraries 
688   AliGenerator* gener = s.MakeGenerator();
689   gener->Init();
690
691   // --- Go back to galice.root --------------------------------------
692   rl->CdGAFile();
693   
694   //=================== Alice BODY parameters =============================
695   AliBODY *BODY = new AliBODY("BODY", "Alice envelop");
696   
697   
698   if (detCfg->UseMAG())   new AliMAG("MAG", "Magnet");
699   if (detCfg->UseABSO())  new AliABSOv3("ABSO", "Muon Absorber");
700   if (detCfg->UseDIPO())  new AliDIPOv3("DIPO", "Dipole version 3");
701   if (detCfg->UseHALL())  new AliHALLv3("HALL", "Alice Hall");
702   if (detCfg->UseFRAME()) (new AliFRAMEv2("FRAME", "Space Frame"))->SetHoles(1);
703   if (detCfg->UseSHIL())  new AliSHILv3("SHIL", "Shielding Version 3");
704   if (detCfg->UsePIPE())  new AliPIPEv3("PIPE", "Beam Pipe");
705   if (detCfg->UseITS())   new AliITSv11("ITS","ITS v11");
706   // if (detCfg->UseITS())   new AliITSv11Hybrid("ITS","ITS v11Hybrid");
707   if (detCfg->UseTPC())   new AliTPCv2("TPC", "Default");
708   if (detCfg->UseTOF())   new AliTOFv6T0("TOF", "normal TOF");
709   if (detCfg->UseHMPID()) new AliHMPIDv3("HMPID", "normal HMPID");
710   if (detCfg->UseZDC()) {
711     AliZDC *ZDC = 0;
712     if (grp->period.EqualTo("LHC10h")) {
713       // Need to use older ZDC for PbPb 
714       ZDC = new AliZDCv3("ZDC", "normal ZDC");
715       ZDC->SetSpectatorsTrack();
716     }
717     else 
718       ZDC = new AliZDCv4("ZDC", "normal ZDC");
719     if (grp->Year() < 2011) { //?
720       // What are these? Do they need to be set properly? 
721       //Collimators aperture
722       ZDC->SetVCollSideCAperture(0.85);
723       ZDC->SetVCollSideCCentre(0.);
724       ZDC->SetVCollSideAAperture(0.75);
725       ZDC->SetVCollSideACentre(0.);
726       //Detector position
727       ZDC->SetYZNC(1.6);
728       ZDC->SetYZNA(1.6);
729       ZDC->SetYZPC(1.6);
730       ZDC->SetYZPA(1.6);
731     }
732     ZDC->SetLumiLength(0.);
733     if (grp->IsPA() || grp->IsAP()) { 
734       ZDC->SetpAsystem();
735       ZDC->SetBeamEnergy(82.*grp->beamEnergy/208.);
736     }
737   }
738   if (detCfg->UseTRD()) {
739     AliTRD *TRD = new AliTRDv1("TRD", "TRD slow simulator");
740     AliTRDgeometry *geoTRD = TRD->GetGeometry();
741     // Total of 18 super modules. We turn them all off by default 
742     for (Int_t i = 0; i < 18; i++) geoTRD->SetSMstatus(i, 0);
743
744     // '09-'10 had 7 super modules 
745     geoTRD->SetSMstatus( 0,1);
746     geoTRD->SetSMstatus( 1,1);
747     geoTRD->SetSMstatus( 7,1);
748     geoTRD->SetSMstatus( 8,1);
749     geoTRD->SetSMstatus( 9,1);
750     geoTRD->SetSMstatus(14,1);//?
751     geoTRD->SetSMstatus(17,1);
752
753     // In jan '11 3 more were added 
754     if (grp->Year() > 2010) { 
755       geoTRD->SetSMstatus(11, 1);
756       geoTRD->SetSMstatus(15, 1);
757       geoTRD->SetSMstatus(16, 1);//?
758     }
759
760     // In the 2012 shutdow 3 more were added 
761     if (grp->Year() > 2012) { 
762       geoTRD->SetSMstatus( 2,1);
763       geoTRD->SetSMstatus( 3,1);
764       geoTRD->SetSMstatus( 6,1);
765     }
766     if (grp->Year() > 2014) {
767       geoTRD->SetSMstatus( 4,1);
768       geoTRD->SetSMstatus( 5,1);
769       geoTRD->SetSMstatus(10,1);
770       geoTRD->SetSMstatus(12,1);
771       geoTRD->SetSMstatus(13,1);
772     }      
773   }
774   if (detCfg->UseFMD())    new AliFMDv1("FMD", "normal FMD");
775   if (detCfg->UseMUON()) {
776     AliMUON *MUON = new AliMUONv1("MUON", "default");
777     MUON->SetTriggerEffCells(1); // not needed if raw masks
778     MUON->SetTriggerResponseV1(2);
779   }
780   if (detCfg->UsePHOS())   new AliPHOSv1("PHOS", "noCPV_Modules123");
781   if (detCfg->UsePMD())    new AliPMDv1("PMD", "normal PMD");
782   if (detCfg->UseT0())     new AliT0v1("T0", "T0 Detector");
783   if (detCfg->UseEMCAL())  new AliEMCALv2("EMCAL", "EMCAL_COMPLETE12SMV1");
784   if (detCfg->UseACORDE()) new AliACORDEv1("ACORDE", "normal ACORDE");
785   if (detCfg->UseVZERO())  new AliVZEROv7("VZERO", "normal VZERO");
786 }
787
788
789
790
791 // 
792 // EOF
793 //