f74ff0157e4d893056e59dbf2f68370f41821ead
[u/mrichter/AliRoot.git] / FlukaConfig.C
1 static Int_t    eventsPerRun = 100;
2 enum PprGeo_t
3 {
4     kHoles, kNoHoles
5 };
6                                                                                 
7 enum PprRad_t
8 {
9     kGluonRadiation, kNoGluonRadiation
10 };
11                                                                                 
12 enum PprMag_t
13 {
14     k2kG, k4kG, k5kG
15 };
16                                                                                 
17                                                                                 
18 // This part for configuration
19 static PprGeo_t sgeo = kHoles;
20 static PprRad_t srad = kGluonRadiation;
21 static PprMag_t smag = k5kG;
22                                                                                 
23 // Comment line
24 static TString  comment;
25                                                                                 
26 // Functions
27 Float_t EtaToTheta(Float_t arg);
28
29
30 void Config()
31 {
32   cout << "==> Config.C..." << endl;
33   
34   // Set Random Number seed
35   gRandom->SetSeed(12345);
36   cout<<"Seed for random number generation= "<<gRandom->GetSeed()<<endl;
37
38   
39   
40   // libraries required by fluka21
41
42   Bool_t isFluka = kTRUE;
43   if (isFluka) {
44     gSystem->Load("libGeom");
45     cout << "\t* Loading TFluka..." << endl;  
46     gSystem->Load("libTFluka");    
47     
48     cout << "\t* Instantiating TFluka..." << endl;
49     new  TFluka("C++ Interface to Fluka", 0/*verbositylevel*/);
50   }
51   else {
52     cout << "\t* Loading Geant3..." << endl;  
53     gSystem->Load("libgeant321");
54     
55     cout << "\t* Instantiating Geant3TGeo..." << endl;
56     new     TGeant3TGeo("C++ Interface to Geant3");
57   }
58   
59   AliRunLoader* rl=0x0;
60                                                                                 
61   cout<<"Config.C: Creating Run Loader ..."<<endl;
62   rl = AliRunLoader::Open("galice.root",
63                             AliConfig::GetDefaultEventFolderName(),
64                             "recreate");
65   if (rl == 0x0)
66     {
67       gAlice->Fatal("Config.C","Can not instatiate the Run Loader");
68       return;
69     }
70   rl->SetCompressionLevel(2);
71   rl->SetNumberOfEventsPerFile(3);
72   gAlice->SetRunLoader(rl);
73                                                                                 
74   //
75   // Set External decayer
76   AliDecayer *decayer = new AliDecayerPythia();
77                                                                                
78   decayer->SetForceDecay(kAll);
79   decayer->Init();
80   gMC->SetExternalDecayer(decayer);
81
82   //
83   //
84   //
85   // Physics process control
86
87   gMC->SetProcess("DCAY",1);
88   gMC->SetProcess("PAIR",1);
89   gMC->SetProcess("COMP",1);
90   gMC->SetProcess("PHOT",1);
91   gMC->SetProcess("PFIS",0);
92   gMC->SetProcess("DRAY",1);
93   gMC->SetProcess("ANNI",1);
94   gMC->SetProcess("BREM",1);
95   gMC->SetProcess("MUNU",1);
96   gMC->SetProcess("CKOV",1); 
97   gMC->SetProcess("HADR",1);
98   gMC->SetProcess("LOSS",2);
99   gMC->SetProcess("MULS",1);
100   gMC->SetProcess("RAYL",1);
101                                                                                 
102   Float_t cut = 1.e-3;        // 1MeV cut by default
103   Float_t tofmax = 1.e10;
104                                                                                 
105   gMC->SetCut("CUTGAM", cut);
106   gMC->SetCut("CUTELE", cut);
107   gMC->SetCut("CUTNEU", cut);
108   gMC->SetCut("CUTHAD", cut);
109   gMC->SetCut("CUTMUO", cut);
110   gMC->SetCut("BCUTE",  cut);
111   gMC->SetCut("BCUTM",  cut);
112   gMC->SetCut("DCUTE",  cut);
113   gMC->SetCut("DCUTM",  cut);
114   gMC->SetCut("PPCUTM", cut);
115   gMC->SetCut("TOFMAX", tofmax);
116
117   //
118   //=======================================================================
119   // ************* STEERING parameters FOR ALICE SIMULATION **************
120   // --- Specify event type to be tracked through the ALICE setup
121   // --- All positions are in cm, angles in degrees, and P and E in GeV
122   if (gSystem->Getenv("CONFIG_NPARTICLES"))
123       int     nParticles = atoi(gSystem->Getenv("CONFIG_NPARTICLES"));
124   else
125       int     nParticles = 1000;
126   
127   cout << "\t* Creating and configuring generator for " << nParticles 
128        << " particles..." << endl;
129   
130   AliGenHIJINGpara *gener = new AliGenHIJINGpara(nParticles);
131   
132   gener->SetMomentumRange(0., 999);
133   gener->SetPhiRange(0, 360);
134   // Set pseudorapidity range from -3 to 3.
135   Float_t thmin = EtaToTheta( 3.);   // theta min. <---> eta max
136   Float_t thmax = EtaToTheta(-3.);   // theta max. <---> eta min 
137   gener->SetThetaRange(thmin,thmax);
138   gener->SetOrigin(0, 0, 0);  //vertex position
139   gener->SetSigma(0, 0, 0);   //Sigma in (X,Y,Z) (cm) on IP position
140   gener->Init();
141   // 
142   // Activate this line if you want the vertex smearing to happen
143   // track by track
144   //
145
146
147      if (smag == k2kG) {
148         comment = comment.Append(" | L3 field 0.2 T");
149     } else if (smag == k4kG) {
150         comment = comment.Append(" | L3 field 0.4 T");
151     } else if (smag == k5kG) {
152         comment = comment.Append(" | L3 field 0.5 T");
153     }
154                                                                                 
155                                                                                 
156     if (srad == kGluonRadiation)
157     {
158         comment = comment.Append(" | Gluon Radiation On");
159                                                                                 
160     } else {
161         comment = comment.Append(" | Gluon Radiation Off");
162     }
163                                                                                 
164     if (sgeo == kHoles)
165     {
166         comment = comment.Append(" | Holes for PHOS/HMPID");
167                                                                                 
168     } else {
169         comment = comment.Append(" | No holes for PHOS/HMPID");
170     }
171                                                                                 
172     printf("\n \n Comment: %s \n \n", comment.Data());
173                                                                                 
174                                                                                 
175 // Field (L3 0.4 T)
176     AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 2, 1., 10., smag);
177     field->SetL3ConstField(0); //Using const. field in the barrel
178     rl->CdGAFile();
179     gAlice->SetField(field);
180  
181   Int_t   iABSO  = 1; 
182   Int_t   iACORDE   = 1; 
183   Int_t   iDIPO  = 1; 
184   Int_t   iFMD   = 1; 
185   Int_t   iFRAME = 1; 
186   Int_t   iHALL  = 1; 
187   Int_t   iITS   = 1; 
188   Int_t   iMAG   = 1; 
189   Int_t   iMUON  = 1; 
190   Int_t   iPHOS  = 1; 
191   Int_t   iPIPE  = 1; 
192   Int_t   iPMD   = 1; 
193   Int_t   iHMPID  = 1; 
194   Int_t   iSHIL  = 1; 
195   Int_t   iT0 = 1; 
196   Int_t   iTOF   = 1; 
197   Int_t   iTPC   = 1;
198   Int_t   iTRD   = 1; 
199   Int_t   iZDC   = 1; 
200   Int_t   iEMCAL = 0; 
201   Int_t   iVZERO = 1;
202  
203   cout << "\t* Creating the detectors ..." << endl;
204   //=================== Alice BODY parameters =============================
205     //=================== Alice BODY parameters =============================
206     AliBODY *BODY = new AliBODY("BODY", "Alice envelop");
207                                                                                 
208
209     if (iMAG)
210     {
211         //=================== MAG parameters ============================
212         // --- Start with Magnet since detector layouts may be depending ---
213         // --- on the selected Magnet dimensions ---
214         AliMAG *MAG = new AliMAG("MAG", "Magnet");
215     }
216
217
218     if (iABSO)
219     {
220         //=================== ABSO parameters ============================
221         AliABSO *ABSO = new AliABSOv0("ABSO", "Muon Absorber");
222     }
223
224     if (iDIPO)
225     {
226         //=================== DIPO parameters ============================
227
228         AliDIPO *DIPO = new AliDIPOv2("DIPO", "Dipole version 2");
229     }
230
231     if (iHALL)
232     {
233         //=================== HALL parameters ============================
234
235         AliHALL *HALL = new AliHALL("HALL", "Alice Hall");
236     }
237
238
239     if (iFRAME)
240     {
241         //=================== FRAME parameters ============================
242
243         AliFRAMEv2 *FRAME = new AliFRAMEv2("FRAME", "Space Frame");
244     }
245
246     if (iSHIL)
247     {
248         //=================== SHIL parameters ============================
249
250         AliSHIL *SHIL = new AliSHILv2("SHIL", "Shielding Version 2");
251     }
252
253
254     if (iPIPE)
255     {
256         //=================== PIPE parameters ============================
257
258         AliPIPE *PIPE = new AliPIPEv0("PIPE", "Beam Pipe");
259     }
260  
261     if(iITS) {
262
263     //=================== ITS parameters ============================
264     //
265     // As the innermost detector in ALICE, the Inner Tracking System "impacts" on
266     // almost all other detectors. This involves the fact that the ITS geometry
267     // still has several options to be followed in parallel in order to determine
268     // the best set-up which minimizes the induced background. All the geometries
269     // available to date are described in the following. Read carefully the comments
270     // and use the default version (the only one uncommented) unless you are making
271     // comparisons and you know what you are doing. In this case just uncomment the
272     // ITS geometry you want to use and run Aliroot.
273     //
274     // Detailed geometries:         
275     //
276     //
277     //AliITS *ITS  = new AliITSv5symm("ITS","Updated ITS TDR detailed version with symmetric services");
278     //
279     //AliITS *ITS  = new AliITSv5asymm("ITS","Updates ITS TDR detailed version with asymmetric services");
280     //
281         AliITSvPPRasymmFMD *ITS  = new AliITSvPPRasymmFMD("ITS","New ITS PPR detailed version with asymmetric services");
282         ITS->SetMinorVersion(2);  // don't touch this parameter if you're not an ITS developer
283         ITS->SetReadDet(kFALSE);          // don't touch this parameter if you're not an ITS developer
284     //    ITS->SetWriteDet("$ALICE_ROOT/ITS/ITSgeometry_vPPRasymm2.det");  // don't touch this parameter if you're not an ITS developer
285         ITS->SetThicknessDet1(200.);   // detector thickness on layer 1 must be in the range [100,300]
286         ITS->SetThicknessDet2(200.);   // detector thickness on layer 2 must be in the range [100,300]
287         ITS->SetThicknessChip1(150.);  // chip thickness on layer 1 must be in the range [150,300]
288         ITS->SetThicknessChip2(150.);  // chip thickness on layer 2 must be in the range [150,300]
289         ITS->SetRails(0);            // 1 --> rails in ; 0 --> rails out
290         ITS->SetCoolingFluid(1);   // 1 --> water ; 0 --> freon
291
292     // Coarse geometries (warning: no hits are produced with these coarse geometries and they unuseful 
293     // for reconstruction !):
294     //                                                     
295     //
296     //AliITSvPPRcoarseasymm *ITS  = new AliITSvPPRcoarseasymm("ITS","New ITS PPR coarse version with asymmetric services");
297     //ITS->SetRails(0);                // 1 --> rails in ; 0 --> rails out
298     //ITS->SetSupportMaterial(0);      // 0 --> Copper ; 1 --> Aluminum ; 2 --> Carbon
299     //
300     //AliITS *ITS  = new AliITSvPPRcoarsesymm("ITS","New ITS PPR coarse version with symmetric services");
301     //ITS->SetRails(0);                // 1 --> rails in ; 0 --> rails out
302     //ITS->SetSupportMaterial(0);      // 0 --> Copper ; 1 --> Aluminum ; 2 --> Carbon
303     //                      
304     //
305     //
306     // Geant3 <-> EUCLID conversion
307     // ============================
308     //
309     // SetEUCLID is a flag to output (=1) or not to output (=0) both geometry and
310     // media to two ASCII files (called by default ITSgeometry.euc and
311     // ITSgeometry.tme) in a format understandable to the CAD system EUCLID.
312     // The default (=0) means that you dont want to use this facility.
313     //
314         ITS->SetEUCLID(0);  
315     }
316
317     if (iTPC)
318     {
319       //============================ TPC parameters =====================
320         AliTPC *TPC = new AliTPCv2("TPC", "Default");
321     }
322
323
324     if (iTOF) {
325         //=================== TOF parameters ============================
326         AliTOF *TOF = new AliTOFv5T0("TOF", "normal TOF");
327         // Partial geometry: modules at 2,3,4,6,7,11,12,14,15,16
328         // starting at 6h in positive direction
329         //      Int_t TOFSectors[18]={-1,-1,0,0,0,-1,0,0,-1,-1,-1,0,0,-1,0,0,0,0};
330         // Partial geometry: modules at 1,2,6,7,9,10,11,12,15,16,17
331         // (ALICE numbering convention)
332         Int_t TOFSectors[18]={-1,0,0,-1,-1,-1,0,0,-1,0,0,0,0,-1,-1,0,0,0};
333         TOF->SetTOFSectors(TOFSectors);
334     }
335
336
337     if (iHMPID)
338     {
339         //=================== HMPID parameters ===========================
340         AliHMPID *HMPID = new AliHMPIDv1("HMPID", "normal HMPID");
341
342     }
343
344
345     if (iZDC)
346     {
347         //=================== ZDC parameters ============================
348
349         AliZDC *ZDC = new AliZDCv2("ZDC", "normal ZDC");
350     }
351
352     if (iTRD)
353     {
354         //=================== TRD parameters ============================
355
356         AliTRD *TRD = new AliTRDv1("TRD", "TRD slow simulator");
357         AliTRDgeometry *geoTRD = TRD->GetGeometry();
358         // Partial geometry: modules at 2,3,4,6,11,12,14,15
359         // starting at 6h in positive direction
360         geoTRD->SetSMstatus(0,0);
361         geoTRD->SetSMstatus(1,0);
362         geoTRD->SetSMstatus(5,0);
363         geoTRD->SetSMstatus(7,0);
364         geoTRD->SetSMstatus(8,0);
365         geoTRD->SetSMstatus(9,0);
366         geoTRD->SetSMstatus(10,0);
367         geoTRD->SetSMstatus(13,0);
368         geoTRD->SetSMstatus(16,0);
369         geoTRD->SetSMstatus(17,0);
370     }
371
372     if (iFMD)
373     {
374         //=================== FMD parameters ============================
375         AliFMD *FMD = new AliFMDv1("FMD", "normal FMD");
376    }
377
378     if (iMUON)
379     {
380         //=================== MUON parameters ===========================
381         // New MUONv1 version (geometry defined via builders)
382         AliMUON *MUON = new AliMUONv1("MUON", "default");
383     }
384     //=================== PHOS parameters ===========================
385
386     if (iPHOS)
387     {
388         AliPHOS *PHOS = new AliPHOSv1("PHOS", "IHEP");
389     }
390
391
392     if (iPMD)
393     {
394         //=================== PMD parameters ============================
395         AliPMD *PMD = new AliPMDv1("PMD", "normal PMD");
396     }
397
398     if (iT0)
399     {
400         //=================== T0 parameters ============================
401         AliT0 *T0 = new AliT0v1("T0", "T0 Detector");
402     }
403
404     if (iEMCAL)
405     {
406         //=================== EMCAL parameters ============================
407         AliEMCAL *EMCAL = new AliEMCALv2("EMCAL", "SHISH_77_TRD1_2X2_FINAL_110DEG");
408     }
409
410      if (iACORDE)
411     {
412         //=================== ACORDE parameters ============================
413         AliACORDE *ACORDE = new AliACORDEv0("ACORDE", "normal ACORDE");
414     }
415
416      if (iVZERO)
417     {
418         //=================== ACORDE parameters ============================
419         AliVZERO *VZERO = new AliVZEROv7("VZERO", "normal VZERO");
420     }                                                                                
421 }
422                                                                                 
423 Float_t EtaToTheta(Float_t arg){
424   return (180./TMath::Pi())*2.*atan(exp(-arg));
425 }
426