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