Configuration in synch with ConfigPPR.C (E. Futo).
[u/mrichter/AliRoot.git] / TFluka / macro / FlukaConfig.C
1 static Int_t    eventsPerRun = 100;
2 enum PprRun_t
3 {
4     test50,
5     kParam_8000,   kParam_4000,  kParam_2000,
6     kHijing_cent1, kHijing_cent2,
7     kHijing_per1,  kHijing_per2, kHijing_per3, kHijing_per4,  kHijing_per5,
8     kHijing_jj25,  kHijing_jj50, kHijing_jj75, kHijing_jj100, kHijing_jj200,
9     kHijing_gj25,  kHijing_gj50, kHijing_gj75, kHijing_gj100, kHijing_gj200,
10     kHijing_pA, kPythia6, kPythia6Jets, kD0PbPb5500, kD_TRD, kB_TRD, kJpsi_TRD,
11     kU_TRD, kPyJJ, kPyGJ
12 };
13                                                                                 
14 enum PprGeo_t
15 {
16     kHoles, kNoHoles
17 };
18                                                                                 
19 enum PprRad_t
20 {
21     kGluonRadiation, kNoGluonRadiation
22 };
23                                                                                 
24 enum PprMag_t
25 {
26     k2kG, k4kG, k5kG
27 };
28                                                                                 
29                                                                                 
30 // This part for configuration
31 //static PprRun_t srun = test50;
32 static PprRun_t srun = kPythia6;
33 static PprGeo_t sgeo = kHoles;
34 static PprRad_t srad = kGluonRadiation;
35 static PprMag_t smag = k5kG;
36                                                                                 
37 // Comment line
38 static TString  comment;
39                                                                                 
40 // Functions
41 Float_t EtaToTheta(Float_t arg);
42 //AliGenerator* GeneratorFactory(PprRun_t srun);
43 //AliGenHijing* HijingStandard();
44
45 void Config()
46 {
47   cout << "==> Config.C..." << endl;
48   
49   // Set Random Number seed
50   gRandom->SetSeed(12345);
51   cout<<"Seed for random number generation= "<<gRandom->GetSeed()<<endl;
52
53   
54   
55   // libraries required by fluka21
56
57   if (!gSystem->Getenv("WITH_ROOT")) {
58      cout << "=== RUNNING TFluka with FLUGG ===\n";
59      char * gvmc = gSystem->ExpandPathName("$(G4VMC)/examples/macro/g4libs.C");
60      gROOT->LoadMacro(gvmc);
61      g4libs();
62
63      cout << "\t* Loading Flugg..." << endl;  
64      gSystem->Load("libFlugg");    
65   } else {
66      cout << "=== RUNNING TFluka with TGeo ===\n";
67      gSystem->Load("libGeom");
68   }     
69
70   cout << "\t* Loading TFluka..." << endl;  
71   gSystem->Load("libTFluka");    
72     
73   cout << "\t* Instantiating TFluka..." << endl;
74   new  TFluka("C++ Interface to Fluka", 3/*verbositylevel*/);
75   
76   AliRunLoader* rl=0x0;
77                                                                                 
78   cout<<"Config.C: Creating Run Loader ..."<<endl;
79   rl = AliRunLoader::Open("galice.root",
80                           AliConfig::fgkDefaultEventFolderName,
81                           "recreate");
82   if (rl == 0x0)
83     {
84       gAlice->Fatal("Config.C","Can not instatiate the Run Loader");
85       return;
86     }
87   rl->SetCompressionLevel(2);
88   rl->SetNumberOfEventsPerFile(3);
89   gAlice->SetRunLoader(rl);
90                                                                                 
91   //
92   // Set External decayer
93   AliDecayer *decayer = new AliDecayerPythia();
94                                                                                
95   decayer->SetForceDecay(kAll);
96   decayer->Init();
97   gMC->SetExternalDecayer(decayer);
98  
99   TFluka *fluka = (TFluka *) gMC;
100   fluka->SetCoreInputFileName("corealice.inp");
101   fluka->SetInputFileName("alice.inp");
102   //
103   //
104   //
105   // Physics process control
106
107   gMC->SetProcess("DCAY",1);
108   gMC->SetProcess("PAIR",1);
109   gMC->SetProcess("COMP",1);
110   gMC->SetProcess("PHOT",1);
111   gMC->SetProcess("PFIS",0);
112   gMC->SetProcess("DRAY",0);
113   gMC->SetProcess("ANNI",1);
114   gMC->SetProcess("BREM",1);
115   gMC->SetProcess("MUNU",1);
116   gMC->SetProcess("CKOV",1);
117   gMC->SetProcess("HADR",1);
118   gMC->SetProcess("LOSS",2);
119   gMC->SetProcess("MULS",1);
120   gMC->SetProcess("RAYL",1);
121                                                                                 
122   Float_t cut = 1.e-3;        // 1MeV cut by default
123   Float_t tofmax = 1.e10;
124                                                                                 
125   gMC->SetCut("CUTGAM", cut);
126   gMC->SetCut("CUTELE", cut);
127   gMC->SetCut("CUTNEU", cut);
128   gMC->SetCut("CUTHAD", cut);
129   gMC->SetCut("CUTMUO", cut);
130   gMC->SetCut("BCUTE",  cut);
131   gMC->SetCut("BCUTM",  cut);
132   gMC->SetCut("DCUTE",  cut);
133   gMC->SetCut("DCUTM",  cut);
134   gMC->SetCut("PPCUTM", cut);
135   gMC->SetCut("TOFMAX", tofmax);
136
137   //
138   //=======================================================================
139   // ************* STEERING parameters FOR ALICE SIMULATION **************
140   // --- Specify event type to be tracked through the ALICE setup
141   // --- All positions are in cm, angles in degrees, and P and E in GeV
142   if (gSystem->Getenv("CONFIG_NPARTICLES"))
143       int     nParticles = atoi(gSystem->Getenv("CONFIG_NPARTICLES"));
144   else
145       int     nParticles = 10;
146   
147   cout << "\t* Creating and configuring generator for " << nParticles 
148        << " particles..." << endl;
149   
150   AliGenHIJINGpara *gener = new AliGenHIJINGpara(nParticles);
151   
152   gener->SetMomentumRange(0, 999);
153   gener->SetPhiRange(0, 360);
154   // Set pseudorapidity range from -8 to 8.
155   Float_t thmin = EtaToTheta( 0.1);   // theta min. <---> eta max
156   Float_t thmax = EtaToTheta(-0.1);  // theta max. <---> eta min 
157   gener->SetThetaRange(thmin,thmax);
158   gener->SetOrigin(0, 0, 0);  //vertex position
159   gener->SetSigma(0, 0, 0);   //Sigma in (X,Y,Z) (cm) on IP position
160   gener->Init();
161   // 
162   // Activate this line if you want the vertex smearing to happen
163   // track by track
164   //
165
166
167   gAlice->SetDebug(10);
168      if (smag == k2kG) {
169         comment = comment.Append(" | L3 field 0.2 T");
170     } else if (smag == k4kG) {
171         comment = comment.Append(" | L3 field 0.4 T");
172     } else if (smag == k5kG) {
173         comment = comment.Append(" | L3 field 0.5 T");
174     }
175                                                                                 
176                                                                                 
177     if (srad == kGluonRadiation)
178     {
179         comment = comment.Append(" | Gluon Radiation On");
180                                                                                 
181     } else {
182         comment = comment.Append(" | Gluon Radiation Off");
183     }
184                                                                                 
185     if (sgeo == kHoles)
186     {
187         comment = comment.Append(" | Holes for PHOS/RICH");
188                                                                                 
189     } else {
190         comment = comment.Append(" | No holes for PHOS/RICH");
191     }
192                                                                                 
193     printf("\n \n Comment: %s \n \n", comment.Data());
194                                                                                 
195                                                                                 
196 // Field (L3 0.4 T)
197     AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 2, 1., 10., smag);
198     field->SetL3ConstField(0); //Using const. field in the barrel
199     rl->CdGAFile();
200     gAlice->SetField(field);
201  
202   Int_t   iABSO  = 0; //1
203   Int_t   iCRT   = 0; //Not good ?
204   Int_t   iDIPO  = 0; //1
205   Int_t   iFMD   = 0; //1
206   Int_t   iFRAME = 0; //1
207   Int_t   iHALL  = 0; //1
208   Int_t   iITS   = 0; //1
209   Int_t   iMAG   = 0; //1
210   Int_t   iMUON  = 0; //1. Not good (newFlagLttc=10000 is outside array bounds)
211   Int_t   iPHOS  = 0; //1
212   Int_t   iPIPE  = 0; //1
213   Int_t   iPMD   = 0; //Not good (too many regions)
214   Int_t   iRICH  = 1; //1. Not good (no tracking with FRAME)
215   Int_t   iSHIL  = 0; //1. Not good (no tracking) (it works alone)
216   Int_t   iSTART = 0; //1. Not good (no tracking) (it works alone)
217   Int_t   iTOF   = 0; //1. Not good (no tracking) (newFlagLttc=10000 is outside array bounds if alone)
218   Int_t   iTPC   = 1;
219   Int_t   iTRD   = 0; //1. Not good (no tracking) (Crash alone with FRAME)
220   Int_t   iZDC   = 0; //1. Needs SHIL and others
221   Int_t   iEMCAL = 0; //Not good (Crash)
222   Int_t   iVZERO = 0;
223  
224   cout << "\t* Creating the detectors ..." << endl;
225   //=================== Alice BODY parameters =============================
226     //=================== Alice BODY parameters =============================
227     AliBODY *BODY = new AliBODY("BODY", "Alice envelop");
228                                                                                 
229                                                                                 
230     if (iMAG)
231     {
232         //=================== MAG parameters ============================
233         // --- Start with Magnet since detector layouts may be depending ---
234         // --- on the selected Magnet dimensions ---
235         AliMAG *MAG = new AliMAG("MAG", "Magnet");
236     }
237                                                                                 
238                                                                                 
239     if (iABSO)
240     {
241         //=================== ABSO parameters ============================
242         AliABSO *ABSO = new AliABSOv0("ABSO", "Muon Absorber");
243     }
244                                                                                 
245     if (iDIPO)
246     {
247         //=================== DIPO parameters ============================
248                                                                                 
249         AliDIPO *DIPO = new AliDIPOv2("DIPO", "Dipole version 2");
250     }
251                                                                                 
252     if (iHALL)
253     {
254         //=================== HALL parameters ============================
255                                                                                 
256         AliHALL *HALL = new AliHALL("HALL", "Alice Hall");
257     }
258                                                                                 
259                                                                                 
260     if (iFRAME)
261     {
262         //=================== FRAME parameters ============================
263                                                                                 
264         AliFRAMEv2 *FRAME = new AliFRAMEv2("FRAME", "Space Frame");
265         if (sgeo == kHoles) {
266             FRAME->SetHoles(1);
267         } else {
268             FRAME->SetHoles(0);
269         }
270     }
271                                                                                 
272     if (iSHIL)
273     {
274         //=================== SHIL parameters ============================
275                                                                                 
276         AliSHIL *SHIL = new AliSHILv2("SHIL", "Shielding Version 2");
277     }
278                                                                                 
279                                                                                 
280     if (iPIPE)
281     {
282         //=================== PIPE parameters ============================
283                                                                                 
284         AliPIPE *PIPE = new AliPIPEv0("PIPE", "Beam Pipe");
285     }
286                                                                                 
287     if(iITS) {
288                                                                                 
289     //=================== ITS parameters ============================
290     //
291     // As the innermost detector in ALICE, the Inner Tracking System "impacts" on
292     // almost all other detectors. This involves the fact that the ITS geometry
293     // still has several options to be followed in parallel in order to determine
294     // the best set-up which minimizes the induced background. All the geometries
295     // available to date are described in the following. Read carefully the comments
296     // and use the default version (the only one uncommented) unless you are making
297     // comparisons and you know what you are doing. In this case just uncomment the
298     // ITS geometry you want to use and run Aliroot.
299     //
300     // Detailed geometries:
301     //
302     //
303     //AliITS *ITS  = new AliITSv5symm("ITS","Updated ITS TDR detailed version with symmetric services");
304     //
305     //AliITS *ITS  = new AliITSv5asymm("ITS","Updates ITS TDR detailed version with asymmetric services");
306     //
307         AliITSvPPRasymmFMD *ITS  = new AliITSvPPRasymmFMD("ITS","New ITS PPR detailed version with asymmetric services");
308         ITS->SetMinorVersion(2);  // don't touch this parameter if you're not an ITS developer
309         ITS->SetReadDet(kTRUE);   // don't touch this parameter if you're not an ITS developer
310     //    ITS->SetWriteDet("$ALICE_ROOT/ITS/ITSgeometry_vPPRasymm2.det");  // don't touch this parameter if you're not an ITS developer
311         ITS->SetThicknessDet1(200.);   // detector thickness on layer 1 must be in the range [100,300]
312         ITS->SetThicknessDet2(200.);   // detector thickness on layer 2 must be in the range [100,300]
313         ITS->SetThicknessChip1(200.);  // chip thickness on layer 1 must be in the range [150,300]
314         ITS->SetThicknessChip2(200.);  // chip thickness on layer 2 must be in the range [150,300]
315         ITS->SetRails(0);            // 1 --> rails in ; 0 --> rails out
316         ITS->SetCoolingFluid(1);   // 1 --> water ; 0 --> freon
317                                                                                 
318     // Coarse geometries (warning: no hits are produced with these coarse geometries and they unuseful
319     // for reconstruction !):
320     //
321     //
322     //AliITSvPPRcoarseasymm *ITS  = new AliITSvPPRcoarseasymm("ITS","New ITS PPR coarse version with asymmetric services");
323     //ITS->SetRails(0);                // 1 --> rails in ; 0 --> rails out
324     //ITS->SetSupportMaterial(0);      // 0 --> Copper ; 1 --> Aluminum ; 2 --> Carbon
325     //
326     //AliITS *ITS  = new AliITSvPPRcoarsesymm("ITS","New ITS PPR coarse version with symmetric services");
327     //ITS->SetRails(0);                // 1 --> rails in ; 0 --> rails out
328     //ITS->SetSupportMaterial(0);      // 0 --> Copper ; 1 --> Aluminum ; 2 --> Carbon
329     //
330     //
331     //
332     // Geant3 <-> EUCLID conversion
333     // ============================
334     //
335     // SetEUCLID is a flag to output (=1) or not to output (=0) both geometry and
336     // media to two ASCII files (called by default ITSgeometry.euc and
337     // ITSgeometry.tme) in a format understandable to the CAD system EUCLID.
338     // The default (=0) means that you dont want to use this facility.
339     //
340         ITS->SetEUCLID(0);
341     }
342                                                                                 
343     if (iTPC)
344     {
345         //============================ TPC parameters ================================
346         // --- This allows the user to specify sectors for the SLOW (TPC geometry 2)
347         // --- Simulator. SecAL (SecAU) <0 means that ALL lower (upper)
348         // --- sectors are specified, any value other than that requires at least one
349         // --- sector (lower or upper)to be specified!
350         // --- Reminder: sectors 1-24 are lower sectors (1-12 -> z>0, 13-24 -> z<0)
351         // ---           sectors 25-72 are the upper ones (25-48 -> z>0, 49-72 -> z<0)
352         // --- SecLows - number of lower sectors specified (up to 6)
353         // --- SecUps - number of upper sectors specified (up to 12)
354         // --- Sens - sensitive strips for the Slow Simulator !!!
355         // --- This does NOT work if all S or L-sectors are specified, i.e.
356         // --- if SecAL or SecAU < 0
357         //
358         //
359         //-----------------------------------------------------------------------------
360                                                                                 
361         //  gROOT->LoadMacro("SetTPCParam.C");
362         //  AliTPCParam *param = SetTPCParam();
363         AliTPC *TPC = new AliTPCv2("TPC", "Default");
364                                                                                 
365         // All sectors included
366         TPC->SetSecAL(-1);
367         TPC->SetSecAU(-1);
368                                                                                 
369     }
370                                                                                 
371                                                                                 
372     if (iTOF) {
373         //=================== TOF parameters ============================
374         AliTOF *TOF = new AliTOFv4T0("TOF", "normal TOF");
375     }
376                                                                                 
377                                                                                 
378     if (iRICH)
379     {
380         //=================== RICH parameters ===========================
381         AliRICH *RICH = new AliRICHv1("RICH", "normal RICH");
382                                                                                 
383     }
384                                                                                 
385                                                                                 
386     if (iZDC)
387     {
388         //=================== ZDC parameters ============================
389                                                                                 
390         AliZDC *ZDC = new AliZDCv2("ZDC", "normal ZDC");
391     }
392                                                                                 
393     if (iTRD)
394     {
395         //=================== TRD parameters ============================
396                                                                                 
397         AliTRD *TRD = new AliTRDv1("TRD", "TRD slow simulator");
398                                                                                 
399         // Select the gas mixture (0: 97% Xe + 3% isobutane, 1: 90% Xe + 10% CO2)
400         TRD->SetGasMix(1);
401         if (sgeo == kHoles) {
402             // With hole in front of PHOS
403             TRD->SetPHOShole();
404             // With hole in front of RICH
405             TRD->SetRICHhole();
406         }
407             // Switch on TR
408             AliTRDsim *TRDsim = TRD->CreateTR();
409     }
410                                                                                 
411     if (iFMD)
412     {
413         //=================== FMD parameters ============================
414         AliFMD *FMD = new AliFMDv1("FMD", "normal FMD");
415    }
416                                                                                 
417     if (iMUON)
418     {
419         //=================== MUON parameters ===========================
420                                                                                 
421         AliMUON *MUON = new AliMUONv1("MUON", "default");
422         MUON->AddGeometryBuilder(new AliMUONSt1GeometryBuilder(MUON));
423         MUON->AddGeometryBuilder(new AliMUONSt2GeometryBuilder(MUON));
424         MUON->AddGeometryBuilder(new AliMUONSlatGeometryBuilder(MUON));
425         MUON->AddGeometryBuilder(new AliMUONTriggerGeometryBuilder(MUON));
426     }
427     //=================== PHOS parameters ===========================
428                                                                                 
429     if (iPHOS)
430     {
431         AliPHOS *PHOS = new AliPHOSv1("PHOS", "IHEP");
432     }
433                                                                                 
434                                                                                 
435     if (iPMD)
436     {
437         //=================== PMD parameters ============================
438         AliPMD *PMD = new AliPMDv1("PMD", "normal PMD");
439     }
440                                                                                 
441     if (iSTART)
442     {
443         //=================== START parameters ============================
444         AliSTART *START = new AliSTARTv1("START", "START Detector");
445     }
446                                                                                 
447     if (iEMCAL)
448     {
449         //=================== EMCAL parameters ============================
450         AliEMCAL *EMCAL = new AliEMCALv1("EMCAL", "EMCAL_55_25");
451     }
452                                                                                 
453      if (iCRT)
454     {
455         //=================== CRT parameters ============================
456         AliCRT *CRT = new AliCRTv0("CRT", "normal ACORDE");
457     }
458                                                                                 
459      if (iVZERO)
460     {
461         //=================== CRT parameters ============================
462         AliVZERO *VZERO = new AliVZEROv3("VZERO", "normal VZERO");
463     }
464                                                                                 
465                                                                                 
466 }
467                                                                                 
468 Float_t EtaToTheta(Float_t arg){
469   return (180./TMath::Pi())*2.*atan(exp(-arg));
470 }
471