0d191a7a8759fa2e5f3ddafc5c6a5359c824d27f
[u/mrichter/AliRoot.git] / test / generators / therminator / Config.C
1 //
2 // Configuration for the first physics production 2008
3 //
4
5 // One can use the configuration macro in compiled mode by
6 // root [0] gSystem->Load("libgeant321");
7 // root [0] gSystem->SetIncludePath("-I$ROOTSYS/include -I$ALICE_ROOT/include\
8 //                   -I$ALICE_ROOT -I$ALICE/geant3/TGeant3");
9 // root [0] .x grun.C(1,"Config.C++")
10
11 #if !defined(__CINT__) || defined(__MAKECINT__)
12 #include <Riostream.h>
13 #include <TRandom.h>
14 #include <TDatime.h>
15 #include <TSystem.h>
16 #include <TVirtualMC.h>
17 #include <TGeant3TGeo.h>
18 #include "STEER/AliRunLoader.h"
19 #include "STEER/AliRun.h"
20 #include "STEER/AliConfig.h"
21 #include "PYTHIA6/AliDecayerPythia.h"
22 #include "PYTHIA6/AliGenPythia.h"
23 #include "TDPMjet/AliGenDPMjet.h"
24 #include "STEER/AliMagFCheb.h"
25 #include "STRUCT/AliBODY.h"
26 #include "STRUCT/AliMAG.h"
27 #include "STRUCT/AliABSOv3.h"
28 #include "STRUCT/AliDIPOv3.h"
29 #include "STRUCT/AliHALLv3.h"
30 #include "STRUCT/AliFRAMEv2.h"
31 #include "STRUCT/AliSHILv3.h"
32 #include "STRUCT/AliPIPEv3.h"
33 #include "ITS/AliITSv11Hybrid.h"
34 #include "TPC/AliTPCv2.h"
35 #include "TOF/AliTOFv6T0.h"
36 #include "HMPID/AliHMPIDv3.h"
37 #include "ZDC/AliZDCv3.h"
38 #include "TRD/AliTRDv1.h"
39 #include "TRD/AliTRDgeometry.h"
40 #include "FMD/AliFMDv1.h"
41 #include "MUON/AliMUONv1.h"
42 #include "PHOS/AliPHOSv1.h"
43 #include "PHOS/AliPHOSSimParam.h"
44 #include "PMD/AliPMDv1.h"
45 #include "T0/AliT0v1.h"
46 #include "EMCAL/AliEMCALv2.h"
47 #include "ACORDE/AliACORDEv1.h"
48 #include "VZERO/AliVZEROv7.h"
49 #endif
50
51 enum PDC06Proc_t
52 {
53   kPythia6, kPhojet, kTherminator, kRunMax
54 };
55
56 const char  * pprRunName[] =
57 {
58   "kPythia6", "kPhojet", "kTherminator"
59 };
60
61 enum Mag_t
62 {
63   kNoField, k5kG, kFieldMax
64 };
65
66 const char * pprField[] =
67 {
68 "kNoField", "k5kG"
69 };
70
71 enum PprTrigConf_t
72 {
73 kDefaultPPTrig, kDefaultPbPbTrig
74 };
75
76 const char * pprTrigConfName[] =
77 {
78 "p-p","Pb-Pb"
79 };
80
81 static PprTrigConf_t strig = kDefaultPbPbTrig;// default PP trigger configuration
82
83 //--- Functions ---
84 class AliGenPythia;
85 AliGenerator *MbPythia();
86 AliGenerator *MbPhojet();
87 AliGenerator *Therminator();
88 void ProcessEnvironmentVars();
89
90 // Geterator, field, beam energy
91 static PDC06Proc_t proc = kTherminator;
92 static Mag_t mag = k5kG;
93 static Float_t energy = 10000; // energy in CMS
94 //======================= =//
95 // Set Random Number seed //
96 //========================//
97 TDatime dt;
98 static UInt_t seed = dt.Get();
99
100 // Comment line
101 static TString comment;
102
103 void Config()
104 {
105
106   // Get settings from environment variables
107   ProcessEnvironmentVars();
108
109   gRandom->SetSeed(seed);
110   cerr<<"Seed for random number generation= "<<seed<<endl;
111
112   // Libraries required by geant321
113 #if defined(__CINT__)
114   gSystem->Load("liblhapdf"); // Parton density functions
115   gSystem->Load("libEGPythia6"); // TGenerator interface
116   gSystem->Load("libpythia6"); // Pythia
117   gSystem->Load("libAliPythia6"); // ALICE specific implementations
118   gSystem->Load("libgeant321");
119 #endif
120
121   new TGeant3TGeo("C++ Interface to Geant3");
122
123   //=======================================================================
124   //  Create the output file
125
126
127   AliRunLoader* rl=0x0;
128
129   cout<<"Config.C: Creating Run Loader ..."<<endl;
130   rl = AliRunLoader::Open("galice.root",
131       AliConfig::GetDefaultEventFolderName(),
132       "recreate");
133   if (rl == 0x0)
134   {
135     gAlice->Fatal("Config.C","Can not instatiate the Run Loader");
136     return;
137   }
138   rl->SetCompressionLevel(2);
139   rl->SetNumberOfEventsPerFile(1000);
140   gAlice->SetRunLoader(rl);
141
142   // Set the trigger configuration: proton-proton
143   gAlice->SetTriggerDescriptor(pprTrigConfName[strig]);
144   cout<<"Trigger configuration is set to  "<<pprTrigConfName[strig]<<endl;
145
146   //
147   //=======================================================================
148   // ************* STEERING parameters FOR ALICE SIMULATION **************
149   // --- Specify event type to be tracked through the ALICE setup
150   // --- All positions are in cm, angles in degrees, and P and E in GeV
151
152
153   gMC->SetProcess("DCAY",1);
154   gMC->SetProcess("PAIR",1);
155   gMC->SetProcess("COMP",1);
156   gMC->SetProcess("PHOT",1);
157   gMC->SetProcess("PFIS",0);
158   gMC->SetProcess("DRAY",0);
159   gMC->SetProcess("ANNI",1);
160   gMC->SetProcess("BREM",1);
161   gMC->SetProcess("MUNU",1);
162   gMC->SetProcess("CKOV",1);
163   gMC->SetProcess("HADR",1);
164   gMC->SetProcess("LOSS",2);
165   gMC->SetProcess("MULS",1);
166   gMC->SetProcess("RAYL",1);
167
168   Float_t cut = 1.e-3; // 1MeV cut by default
169   Float_t tofmax = 1.e10;
170
171   gMC->SetCut("CUTGAM", cut);
172   gMC->SetCut("CUTELE", cut);
173   gMC->SetCut("CUTNEU", cut);
174   gMC->SetCut("CUTHAD", cut);
175   gMC->SetCut("CUTMUO", cut);
176   gMC->SetCut("BCUTE", cut);
177   gMC->SetCut("BCUTM", cut);
178   gMC->SetCut("DCUTE", cut);
179   gMC->SetCut("DCUTM", cut);
180   gMC->SetCut("PPCUTM", cut);
181   gMC->SetCut("TOFMAX", tofmax);
182
183   //======================//
184   // Set External decayer //
185   //======================//
186   TVirtualMCDecayer* decayer = new AliDecayerPythia();
187   decayer->SetForceDecay(kAll);
188   decayer->Init();
189   gMC->SetExternalDecayer(decayer);
190
191   //=========================//
192   // Generator Configuration //
193   //=========================//
194   AliGenerator* gener = 0x0;
195
196   switch (proc)
197   {
198     case kPythia6:
199     gener = MbPythia();
200     break;
201     case kPhojet:
202     gener = MbPhojet();
203     break;
204     case kTherminator:
205     gener = Therminator();
206     break;
207   }
208
209   // PRIMARY VERTEX
210   //
211   gener->SetOrigin(0., 0., 0.); // vertex position
212   //
213   // Size of the interaction diamond
214   // Longitudinal
215   Float_t sigmaz = 5.4 / TMath::Sqrt(2.); // [cm]
216   if (energy == 900)
217   sigmaz = 10.5 / TMath::Sqrt(2.); // [cm]
218   //
219   // Transverse
220   Float_t betast = 10; // beta* [m]
221   Float_t eps = 3.75e-6; // emittance [m]
222   Float_t gamma = energy / 2.0 / 0.938272; // relativistic gamma [1]
223   Float_t sigmaxy = TMath::Sqrt(eps * betast / gamma) / TMath::Sqrt(2.) * 100.; // [cm]
224   printf("\n \n Diamond size x-y: %10.3e z: %10.3e\n \n", sigmaxy, sigmaz);
225
226   gener->SetSigma(sigmaxy, sigmaxy, sigmaz); // Sigma in (X,Y,Z) (cm) on IP position
227   gener->SetCutVertexZ(3.); // Truncate at 3 sigma
228   gener->SetVertexSmear(kPerEvent);
229
230   gener->Init();
231
232   // FIELD
233   AliMagF* field = 0x0;
234
235   if (mag == kNoField)
236   {
237     comment = comment.Append(" | L3 field 0.0 T");
238     field = new AliMagF("Maps","Maps", -1., -1., AliMagF::k2kG);
239   }
240   else if (mag == k5kG)
241   {
242     comment = comment.Append(" | L3 field 0.5 T");
243     field = new AliMagF("Maps","Maps", -1., -1., AliMagF::k5kG);
244   }
245   printf("\n \n Comment: %s \n \n", comment.Data());
246
247   TGeoGlobalMagField::Instance()->SetField(field);
248
249   rl->CdGAFile();
250
251   Int_t iABSO = 1;
252   Int_t iACORDE= 0;
253   Int_t iDIPO = 1;
254   Int_t iEMCAL = 0;
255   Int_t iFMD = 1;
256   Int_t iFRAME = 1;
257   Int_t iHALL = 1;
258   Int_t iITS = 1;
259   Int_t iMAG = 1;
260   Int_t iMUON = 1;
261   Int_t iPHOS = 1;
262   Int_t iPIPE = 1;
263   Int_t iPMD = 0;
264   Int_t iHMPID = 1;
265   Int_t iSHIL = 1;
266   Int_t iT0 = 1;
267   Int_t iTOF = 1;
268   Int_t iTPC = 1;
269   Int_t iTRD = 1;
270   Int_t iVZERO = 1;
271   Int_t iZDC = 1;
272
273   //=================== Alice BODY parameters =============================
274   AliBODY *BODY = new AliBODY("BODY", "Alice envelop");
275
276   if (iMAG)
277   {
278     //=================== MAG parameters ============================
279     // --- Start with Magnet since detector layouts may be depending ---
280     // --- on the selected Magnet dimensions ---
281     AliMAG *MAG = new AliMAG("MAG", "Magnet");
282   }
283
284   if (iABSO)
285   {
286     //=================== ABSO parameters ============================
287     AliABSO *ABSO = new AliABSOv3("ABSO", "Muon Absorber");
288   }
289
290   if (iDIPO)
291   {
292     //=================== DIPO parameters ============================
293
294     AliDIPO *DIPO = new AliDIPOv3("DIPO", "Dipole version 3");
295   }
296
297   if (iHALL)
298   {
299     //=================== HALL parameters ============================
300
301     AliHALL *HALL = new AliHALLv3("HALL", "Alice Hall");
302   }
303
304   if (iFRAME)
305   {
306     //=================== FRAME parameters ============================
307
308     AliFRAMEv2 *FRAME = new AliFRAMEv2("FRAME", "Space Frame");
309     FRAME->SetHoles(1);
310   }
311
312   if (iSHIL)
313   {
314     //=================== SHIL parameters ============================
315
316     AliSHIL *SHIL = new AliSHILv3("SHIL", "Shielding Version 3");
317   }
318
319   if (iPIPE)
320   {
321     //=================== PIPE parameters ============================
322
323     AliPIPE *PIPE = new AliPIPEv3("PIPE", "Beam Pipe");
324   }
325
326   if (iITS)
327   {
328     //=================== ITS parameters ============================
329
330     AliITS *ITS = new AliITSv11Hybrid("ITS","ITS v11Hybrid");
331   }
332
333   if (iTPC)
334   {
335     //============================ TPC parameters =====================
336
337     AliTPC *TPC = new AliTPCv2("TPC", "Default");
338   }
339
340   if (iTOF)
341   {
342     //=================== TOF parameters ============================
343
344     AliTOF *TOF = new AliTOFv6T0("TOF", "normal TOF");
345   }
346
347   if (iHMPID)
348   {
349     //=================== HMPID parameters ===========================
350
351     AliHMPID *HMPID = new AliHMPIDv3("HMPID", "normal HMPID");
352
353   }
354
355   if (iZDC)
356   {
357     //=================== ZDC parameters ============================
358
359     AliZDC *ZDC = new AliZDCv3("ZDC", "normal ZDC");
360   }
361
362   if (iTRD)
363   {
364     //=================== TRD parameters ============================
365
366     AliTRD *TRD = new AliTRDv1("TRD", "TRD slow simulator");
367     AliTRDgeometry *geoTRD = TRD->GetGeometry();
368     // Partial geometry: modules at 0,1,7,8,9,10,17
369     // starting at 3h in positive direction
370     geoTRD->SetSMstatus(2,0);
371     geoTRD->SetSMstatus(3,0);
372     geoTRD->SetSMstatus(4,0);
373     geoTRD->SetSMstatus(5,0);
374     geoTRD->SetSMstatus(6,0);
375     geoTRD->SetSMstatus(11,0);
376     geoTRD->SetSMstatus(12,0);
377     geoTRD->SetSMstatus(13,0);
378     geoTRD->SetSMstatus(14,0);
379     geoTRD->SetSMstatus(15,0);
380     geoTRD->SetSMstatus(16,0);
381   }
382
383   if (iFMD)
384   {
385     //=================== FMD parameters ============================
386
387     AliFMD *FMD = new AliFMDv1("FMD", "normal FMD");
388   }
389
390   if (iMUON)
391   {
392     //=================== MUON parameters ===========================
393     // New MUONv1 version (geometry defined via builders)
394
395     AliMUON *MUON = new AliMUONv1("MUON", "default");
396   }
397
398   if (iPHOS)
399   {
400     //=================== PHOS parameters ===========================
401
402     AliPHOS *PHOS = new AliPHOSv1("PHOS", "IHEP");
403     //Set simulation parameters different from the default ones.
404     AliPHOSSimParam* simEmc = AliPHOSSimParam::GetInstance();
405
406     // APD noise of warm (+20C) PHOS:
407     // a2 = a1*(Y1/Y2)*(M1/M2), where a1 = 0.012 is APD noise at -25C,
408     // Y1 = 4.3 photo-electrons/MeV, Y2 = 1.7 p.e/MeV - light yields at -25C and +20C,
409     // M1 = 50, M2 = 50 - APD gain factors chosen for t1 = -25C and t2 = +20C,
410     // Y = MeanLightYield*APDEfficiency.
411
412     Float_t apdNoise = 0.012*2.5;
413     simEmc->SetAPDNoise(apdNoise);
414
415     //Raw Light Yield at +20C
416     simEmc->SetMeanLightYield(18800);
417
418     //ADC channel width at +18C.
419     simEmc->SetADCchannelW(0.0125);
420   }
421
422   if (iPMD)
423   {
424     //=================== PMD parameters ============================
425
426     AliPMD *PMD = new AliPMDv1("PMD", "normal PMD");
427       }
428
429       if (iT0)
430 {
431   //=================== T0 parameters ============================
432 AliT0  *T0 = new AliT0v1("T0", "T0 Detector");
433 }
434
435 if (iEMCAL)
436 {
437   //=================== EMCAL parameters ============================
438
439   AliEMCAL *EMCAL = new AliEMCALv2("EMCAL", "EMCAL_COMPLETE");
440 }
441
442 if (iACORDE)
443 {
444   //  =================== ACORDE parameters ============================
445
446 AliACORDE *ACORDE = new AliACORDEv1("ACORDE", "normal ACORDE");
447 }
448
449 if (iVZERO)
450 {
451 //=================== ACORDE parameters ============================
452
453 AliVZERO *VZERO = new AliVZEROv7("VZERO", "normal VZERO");
454 }
455 }
456 //
457 //           PYTHIA
458 //
459
460 AliGenerator* MbPythia()
461 {
462 comment = comment.Append(" pp at 14 TeV: Pythia low-pt");
463 //
464 //    Pythia
465 AliGenPythia* pythia = new AliGenPythia(-1);
466 pythia->SetMomentumRange(0, 999999.);
467 pythia->SetThetaRange(0., 180.);
468 pythia->SetYRange(-12.,12.);
469 pythia->SetPtRange(0,1000.);
470 pythia->SetProcess(kPyMb);
471 pythia->SetEnergyCMS(energy);
472
473 return pythia;
474 }
475
476 AliGenerator* MbPhojet()
477 {
478 comment = comment.Append(" pp at 14 TeV: Phojet low-pt");
479 //
480 //    DPMJET
481 #if defined(__CINT__)
482 gSystem->Load("libdpmjet"); // Parton density functions
483 gSystem->Load("libTDPMjet"); // Parton density functions
484 #endif
485 AliGenDPMjet* dpmjet = new AliGenDPMjet(-1);
486 dpmjet->SetMomentumRange(0, 999999.);
487 dpmjet->SetThetaRange(0., 180.);
488 dpmjet->SetYRange(-12.,12.);
489 dpmjet->SetPtRange(0,1000.);
490 dpmjet->SetProcess(kDpmMb);
491 dpmjet->SetEnergyCMS(energy);
492
493 return dpmjet;
494 }
495
496 AliGenerator* Therminator()
497 {
498 comment = comment.Append(" pp at 14 TeV: Therminator");
499
500 #if defined(__CINT__)
501 gSystem->Load("libTTherminator"); // Therminator library
502 #endif
503 AliGenTherminator *genther = new AliGenTherminator();
504 genther->SetFileName("event.out");
505 genther->SetEventNumberInFile(1);
506 genther->SetTemperature(0.145);
507 genther->SetMiuI(-0.0009);
508 genther->SetMiuS(0.000);
509 genther->SetMiuB(0.0008);
510 genther->SetAlfaRange(8.0);
511 genther->SetRapRange(4.0);
512 genther->SetRhoMax(7.74);
513 genther->SetTau(9.74);
514 genther->SetModel("Lhyquid3D");
515 genther->SetLhyquidSet("LHC500C2030");
516
517 return genther;
518 }
519
520 void ProcessEnvironmentVars()
521 {
522 // Run type
523 if (gSystem->Getenv("CONFIG_RUN_TYPE"))
524 {
525   for (Int_t iRun = 0; iRun < kRunMax; iRun++)
526   {
527     if (strcmp(gSystem->Getenv("CONFIG_RUN_TYPE"), pprRunName[iRun])==0)
528     {
529       proc = (PDC06Proc_t)iRun;
530       cout<<"Run type set to "<<pprRunName[iRun]<<endl;
531     }
532   }
533 }
534
535 // Field
536 if (gSystem->Getenv("CONFIG_FIELD"))
537 {
538   for (Int_t iField = 0; iField < kFieldMax; iField++)
539   {
540     if (strcmp(gSystem->Getenv("CONFIG_FIELD"), pprField[iField])==0)
541     {
542       mag = (Mag_t)iField;
543       cout<<"Field set to "<<pprField[iField]<<endl;
544     }
545   }
546 }
547
548 // Energy
549 if (gSystem->Getenv("CONFIG_ENERGY"))
550 {
551   energy = atoi(gSystem->Getenv("CONFIG_ENERGY"));
552   cout<<"Energy set to "<<energy<<" GeV"<<endl;
553 }
554
555 // Random Number seed
556 if (gSystem->Getenv("CONFIG_SEED"))
557 {
558   seed = atoi(gSystem->Getenv("CONFIG_SEED"));
559 }
560 }