Update of calibration task and replace AliTRDrunCalib.C by AddTaskTRDCalib.C (Raphaelle)
[u/mrichter/AliRoot.git] / TRD / Macros / AliTRDConfigPID.C
1 //
2 // Configuration for the Physics Data Challenge 2006 modified to generate
3 // e,mu,pi,K,p for PID studies of the TRD.
4 // Per event, AliGenBox is used to generate 100 particles per species
5 // (50 positive and 50 negative ones). Barrel detectors only.
6 // s.masciocchi@gsi.de
7 //
8
9 // One can use the configuration macro in compiled mode by
10 // root [0] gSystem->Load("libgeant321");
11 // root [0] gSystem->SetIncludePath("-I$ROOTSYS/include -I$ALICE_ROOT/include\
12 //                   -I$ALICE_ROOT -I$ALICE/geant3/TGeant3");
13 // root [0] .x grun.C(1,"Config_PDC06.C++")
14
15 #if !defined(__CINT__) || defined(__MAKECINT__)
16 #include <Riostream.h>
17 #include <TPDGCode.h>
18 #include <TRandom.h>
19 #include <TDatime.h>
20 #include <TSystem.h>
21 #include <TVirtualMC.h>
22 #include <TGeant3TGeo.h>
23 #include "EVGEN/AliGenCocktail.h"
24 #include "EVGEN/AliGenParam.h"
25 #include "EVGEN/AliGenMUONlib.h"
26 #include "STEER/AliRunLoader.h"
27 #include "STEER/AliRun.h"
28 #include "STEER/AliConfig.h"
29 #include "PYTHIA6/AliDecayerPythia.h"
30 #include "PYTHIA6/AliGenPythia.h"
31 #include "STEER/AliMagF.h"
32 #include "STRUCT/AliBODY.h"
33 #include "STRUCT/AliMAG.h"
34 #include "STRUCT/AliABSOv0.h"
35 #include "STRUCT/AliDIPOv2.h"
36 #include "STRUCT/AliHALL.h"
37 #include "STRUCT/AliFRAMEv2.h"
38 #include "STRUCT/AliSHILv2.h"
39 #include "STRUCT/AliPIPEv0.h"
40 #include "ITS/AliITSgeom.h"
41 #include "ITS/AliITSvPPRasymmFMD.h"
42 #include "TPC/AliTPCv2.h"
43 #include "TOF/AliTOFv5T0.h"
44 #include "HMPID/AliHMPIDv1.h"
45 #include "ZDC/AliZDCv2.h"
46 #include "TRD/AliTRDv1.h"
47 #include "FMD/AliFMDv1.h"
48 #include "MUON/AliMUONv1.h"
49 #include "PHOS/AliPHOSv1.h"
50 #include "PMD/AliPMDv1.h"
51 #include "T0/AliT0v1.h"
52 #include "EMCAL/AliEMCALv2.h"
53 #include "ACORDE/AliACORDEv1.h"
54 #include "VZERO/AliVZEROv6.h"
55 #endif
56
57
58 enum PDC06Proc_t 
59 {
60 //--- Heavy Flavour Production ---
61   kCharmPbPb5500,  kCharmpPb8800,  kCharmpp14000,  kCharmpp14000wmi,
62   kD0PbPb5500,     kD0pPb8800,     kD0pp14000,
63   kDPlusPbPb5500,  kDPluspPb8800,  kDPluspp14000,
64   kBeautyPbPb5500, kBeautypPb8800, kBeautypp14000, kBeautypp14000wmi, 
65 // -- Pythia Mb
66   kPyMbNoHvq, kPyOmegaPlus, kPyOmegaMinus, kRunMax
67 }; 
68
69 const char * pprRunName[] = {
70   "kCharmPbPb5500",  "kCharmpPb8800",  "kCharmpp14000",  "kCharmpp14000wmi",
71   "kD0PbPb5500",     "kD0pPb8800",     "kD0pp14000",
72   "kDPlusPbPb5500",  "kDPluspPb8800",  "kDPluspp14000",
73   "kBeautyPbPb5500", "kBeautypPb8800", "kBeautypp14000", "kBeautypp14000wmi", 
74   "kPyMbNoHvq", "kPyOmegaPlus", "kPyOmegaMinus"
75 };
76
77
78 //--- Decay Mode ---
79 enum DecayHvFl_t 
80 {
81   kNature,  kHadr, kSemiEl, kSemiMu
82 };
83 //--- Rapidity Cut ---
84 enum YCut_t
85 {
86   kFull, kBarrel, kMuonArm
87 };
88
89 //--- Functions ---
90 AliGenPythia *PythiaHVQ(PDC06Proc_t proc);
91 AliGenerator *MbCocktail();
92 AliGenerator *PyMbTriggered(Int_t pdg);
93 void ProcessEnvironmentVars();
94 Float_t EtaToTheta(Float_t arg){
95   return (180./TMath::Pi())*2.*atan(exp(-arg));
96 }
97
98 // This part for configuration
99 static DecayHvFl_t   decHvFl   = kNature; 
100 static YCut_t        ycut      = kBarrel;
101 static AliMagF::BMap_t         mag       = k5kG; 
102 //========================//
103 // Set Random Number seed //
104 //========================//
105 TDatime dt;
106 static UInt_t seed    = dt.Get();
107
108 // nEvts = -1  : you get 1 QQbar pair and all the fragmentation and 
109 //               decay chain
110 // nEvts = N>0 : you get N charm / beauty Hadrons 
111 Int_t nEvts = -1; 
112 // stars = kTRUE : all heavy resonances and their decay stored
113 //       = kFALSE: only final heavy hadrons and their decays stored
114 Bool_t stars = kTRUE;
115
116 // To be used only with kCharmppMNRwmi and kBeautyppMNRwmi
117 // To get a "reasonable" agreement with MNR results, events have to be 
118 // generated with the minimum ptHard set to 2.76 GeV.
119 // To get a "perfect" agreement with MNR results, events have to be 
120 // generated in four ptHard bins with the following relative 
121 // normalizations:
122 //  CHARM
123 // 2.76-3 GeV: 25%
124 //    3-4 GeV: 40%
125 //    4-8 GeV: 29%
126 //     >8 GeV:  6%
127 //  BEAUTY
128 // 2.76-4 GeV:  5% 
129 //    4-6 GeV: 31%
130 //    6-8 GeV: 28%
131 //     >8 GeV: 36%
132 Float_t ptHardMin =  2.76;
133 Float_t ptHardMax = -1.;
134
135
136 // Comment line
137 static TString comment;
138
139 void Config()
140 {
141  
142
143   // Get settings from environment variables
144   ProcessEnvironmentVars();
145
146   gRandom->SetSeed(seed);
147   cerr<<"Seed for random number generation= "<<seed<<endl; 
148
149   // libraries required by geant321
150 #if defined(__CINT__)
151   gSystem->Load("libgeant321");
152 #endif
153
154   new TGeant3TGeo("C++ Interface to Geant3");
155
156   //=======================================================================
157   //  Create the output file
158
159    
160   AliRunLoader* rl=0x0;
161
162   cout<<"Config.C: Creating Run Loader ..."<<endl;
163   rl = AliRunLoader::Open("galice.root",
164                           AliConfig::GetDefaultEventFolderName(),
165                           "recreate");
166   if (rl == 0x0)
167     {
168       gAlice->Fatal("Config.C","Can not instatiate the Run Loader");
169       return;
170     }
171   rl->SetCompressionLevel(2);
172   rl->SetNumberOfEventsPerFile(1000);
173   gAlice->SetRunLoader(rl);
174
175   //
176   //=======================================================================
177   // ************* STEERING parameters FOR ALICE SIMULATION **************
178   // --- Specify event type to be tracked through the ALICE setup
179   // --- All positions are in cm, angles in degrees, and P and E in GeV
180
181
182     gMC->SetProcess("DCAY",1);
183     gMC->SetProcess("PAIR",1);
184     gMC->SetProcess("COMP",1);
185     gMC->SetProcess("PHOT",1);
186     gMC->SetProcess("PFIS",0);
187     gMC->SetProcess("DRAY",0);
188     gMC->SetProcess("ANNI",1);
189     gMC->SetProcess("BREM",1);
190     gMC->SetProcess("MUNU",1);
191     gMC->SetProcess("CKOV",1);
192     gMC->SetProcess("HADR",1);
193     gMC->SetProcess("LOSS",2);
194     gMC->SetProcess("MULS",1);
195     gMC->SetProcess("RAYL",1);
196
197     Float_t cut = 1.e-3;        // 1MeV cut by default
198     Float_t tofmax = 1.e10;
199
200     gMC->SetCut("CUTGAM", cut);
201     gMC->SetCut("CUTELE", cut);
202     gMC->SetCut("CUTNEU", cut);
203     gMC->SetCut("CUTHAD", cut);
204     gMC->SetCut("CUTMUO", cut);
205     gMC->SetCut("BCUTE",  cut); 
206     gMC->SetCut("BCUTM",  cut); 
207     gMC->SetCut("DCUTE",  cut); 
208     gMC->SetCut("DCUTM",  cut); 
209     gMC->SetCut("PPCUTM", cut);
210     gMC->SetCut("TOFMAX", tofmax); 
211
212
213
214
215   // Set External decayer //
216   //======================//
217   TVirtualMCDecayer* decayer = new AliDecayerPythia();
218   // DECAYS
219   //
220   switch(decHvFl) {
221   case kNature:
222     decayer->SetForceDecay(kAll);
223     break;
224   case kHadr:
225     decayer->SetForceDecay(kHadronicD);
226     break;
227   case kSemiEl:
228     decayer->SetForceDecay(kSemiElectronic);
229     break;
230   case kSemiMu:
231     decayer->SetForceDecay(kSemiMuonic);
232     break;
233   }
234   decayer->Init();
235   gMC->SetExternalDecayer(decayer);
236
237   //=========================//
238   // Generator Configuration //
239   //=========================//
240   
241   AliGenCocktail *gener  = new AliGenCocktail();
242   Double_t momen         = 0.6;
243   Int_t Npart            = 50;
244   // Set pseudorapidity range from -1 to 1.
245   Float_t thmin          = EtaToTheta(1.);   // theta min. <---> eta max
246   Float_t thmax          = EtaToTheta(-1.);  // theta max. <---> eta min
247   
248   // electron generator
249   AliGenBox *eminus= new AliGenBox(Npart);
250   eminus->SetPart(kElectron);
251   eminus->SetMomentumRange(momen,momen);
252   eminus->SetPhiRange(0,360);
253   eminus->SetThetaRange(thmin,thmax);
254   eminus->SetOrigin(0., 0., 0.);    // primary vertex position
255   eminus->SetSigma(0.,0.,0.);       // sigma (x,y,z) in cm, on position
256
257   // positron generator
258   AliGenBox *eplus= new AliGenBox(Npart);
259   eplus->SetPart(kPositron);
260   eplus->SetMomentumRange(momen,momen);
261   eplus->SetPhiRange(0,360);
262   eplus->SetThetaRange(thmin,thmax);
263   eplus->SetOrigin(0., 0., 0.);    // primary vertex position
264   eplus->SetSigma(0.,0.,0.);       // sigma (x,y,z) in cm, on position
265
266   // mu- generator
267   AliGenBox *muminus= new AliGenBox(Npart);
268   muminus->SetPart(kMuonMinus);
269   muminus->SetMomentumRange(momen,momen);
270   muminus->SetPhiRange(0,360);
271   muminus->SetThetaRange(thmin,thmax);
272   muminus->SetOrigin(0., 0., 0.);    // primary vertex position
273   muminus->SetSigma(0.,0.,0.);       // sigma (x,y,z) in cm, on position
274
275   // mu+ generator
276   AliGenBox *muplus= new AliGenBox(Npart);
277   muplus->SetPart(kMuonPlus);
278   muplus->SetMomentumRange(momen,momen);
279   muplus->SetPhiRange(0,360);
280   muplus->SetThetaRange(thmin,thmax);
281   muplus->SetOrigin(0., 0., 0.);    // primary vertex position
282   muplus->SetSigma(0.,0.,0.);       // sigma (x,y,z) in cm, on position
283
284   // pi- generator
285   AliGenBox *piminus= new AliGenBox(Npart);
286   piminus->SetPart(kPiMinus);
287   piminus->SetMomentumRange(momen,momen);
288   piminus->SetPhiRange(0,360);
289   piminus->SetThetaRange(thmin,thmax);
290   piminus->SetOrigin(0., 0., 0.);    // primary vertex position
291   piminus->SetSigma(0.,0.,0.);       // sigma (x,y,z) in cm, on position
292
293   // pi+ generator
294   AliGenBox *piplus= new AliGenBox(Npart);
295   piplus->SetPart(kPiPlus);
296   piplus->SetMomentumRange(momen,momen);
297   piplus->SetPhiRange(0,360);
298   piplus->SetThetaRange(thmin,thmax);
299   piplus->SetOrigin(0., 0., 0.);    // primary vertex position
300   piplus->SetSigma(0.,0.,0.);       // sigma (x,y,z) in cm, on position
301
302   // K- generator
303   AliGenBox *kminus= new AliGenBox(Npart);
304   kminus->SetPart(kKMinus);
305   kminus->SetMomentumRange(momen,momen);
306   kminus->SetPhiRange(0,360);
307   kminus->SetThetaRange(thmin,thmax);
308   kminus->SetOrigin(0., 0., 0.);    // primary vertex position
309   kminus->SetSigma(0.,0.,0.);       // sigma (x,y,z) in cm, on position
310
311   // K+ generator
312   AliGenBox *kplus= new AliGenBox(Npart);
313   kplus->SetPart(kKPlus);
314   kplus->SetMomentumRange(momen,momen);
315   kplus->SetPhiRange(0,360);
316   kplus->SetThetaRange(thmin,thmax);
317   kplus->SetOrigin(0., 0., 0.);    // primary vertex position
318   kplus->SetSigma(0.,0.,0.);       // sigma (x,y,z) in cm, on position
319
320   // p generator
321   AliGenBox *proton= new AliGenBox(Npart);
322   proton->SetPart(kProton);
323   proton->SetMomentumRange(momen,momen);
324   proton->SetPhiRange(0,360);
325   proton->SetThetaRange(thmin,thmax);
326   proton->SetOrigin(0., 0., 0.);    // primary vertex position
327   proton->SetSigma(0.,0.,0.);       // sigma (x,y,z) in cm, on position
328
329   // anti-p generator
330   AliGenBox *aproton= new AliGenBox(Npart);
331   aproton->SetPart(kProtonBar);
332   aproton->SetMomentumRange(momen,momen);
333   aproton->SetPhiRange(0,360);
334   aproton->SetThetaRange(thmin,thmax);
335   aproton->SetOrigin(0., 0., 0.);    // primary vertex position
336   aproton->SetSigma(0.,0.,0.);       // sigma (x,y,z) in cm, on position
337
338
339
340   gener->AddGenerator(eminus,"e-",1);
341   gener->AddGenerator(eplus,"e+",1);
342   gener->AddGenerator(muminus,"mu-",1);
343   gener->AddGenerator(muplus,"mu+",1);
344   gener->AddGenerator(piminus,"pi-",1);
345   gener->AddGenerator(piplus,"pi+",1);
346   gener->AddGenerator(kminus,"K-",1);
347   gener->AddGenerator(kplus,"K+",1);
348   gener->AddGenerator(proton,"p",1);
349   gener->AddGenerator(aproton,"anti-p",1);
350
351   gener->Init();
352
353   // FIELD
354   //    
355   if (mag == AliMagF::k2kG) {
356     comment = comment.Append(" | L3 field 0.2 T");
357   } else if (mag == AliMagF::k5kG) {
358     comment = comment.Append(" | L3 field 0.5 T");
359   }
360   printf("\n \n Comment: %s \n \n", comment.Data());
361     
362   AliMagF* field = new AliMagF("Maps","Maps", 1., 1., mag);
363   field->SetL3ConstField(0); //Using const. field in the barrel
364   TGeoGlobalMagField::Instance()->SetField(field);
365   rl->CdGAFile();
366
367   Int_t iABSO  = 1;
368   Int_t iACORDE   = 0;
369   Int_t iDIPO  = 0;
370   Int_t iEMCAL = 0;
371   Int_t iFMD   = 0;
372   Int_t iFRAME = 1;
373   Int_t iHALL  = 1;
374   Int_t iITS   = 1;
375   Int_t iMAG   = 1;
376   Int_t iMUON  = 0;
377   Int_t iPHOS  = 1;
378   Int_t iPIPE  = 1;
379   Int_t iPMD   = 0;
380   Int_t iHMPID  = 1;
381   Int_t iSHIL  = 0;
382   Int_t iT0 = 0;
383   Int_t iTOF   = 1;
384   Int_t iTPC   = 1;
385   Int_t iTRD   = 1;
386   Int_t iVZERO = 0;
387   Int_t iZDC   = 0;
388
389
390     //=================== Alice BODY parameters =============================
391     AliBODY *BODY = new AliBODY("BODY", "Alice envelop");
392
393
394     if (iMAG)
395     {
396         //=================== MAG parameters ============================
397         // --- Start with Magnet since detector layouts may be depending ---
398         // --- on the selected Magnet dimensions ---
399         AliMAG *MAG = new AliMAG("MAG", "Magnet");
400     }
401
402
403     if (iABSO)
404     {
405         //=================== ABSO parameters ============================
406         AliABSO *ABSO = new AliABSOv0("ABSO", "Muon Absorber");
407     }
408
409     if (iDIPO)
410     {
411         //=================== DIPO parameters ============================
412
413         AliDIPO *DIPO = new AliDIPOv2("DIPO", "Dipole version 2");
414     }
415
416     if (iHALL)
417     {
418         //=================== HALL parameters ============================
419
420         AliHALL *HALL = new AliHALL("HALL", "Alice Hall");
421     }
422
423
424     if (iFRAME)
425     {
426         //=================== FRAME parameters ============================
427
428         AliFRAMEv2 *FRAME = new AliFRAMEv2("FRAME", "Space Frame");
429         FRAME->SetHoles(0);
430     }
431
432     if (iSHIL)
433     {
434         //=================== SHIL parameters ============================
435
436         AliSHIL *SHIL = new AliSHILv2("SHIL", "Shielding Version 2");
437     }
438
439
440     if (iPIPE)
441     {
442         //=================== PIPE parameters ============================
443
444         AliPIPE *PIPE = new AliPIPEv0("PIPE", "Beam Pipe");
445     }
446  
447     if(iITS) {
448
449     //=================== ITS parameters ============================
450     //
451     // As the innermost detector in ALICE, the Inner Tracking System "impacts" on
452     // almost all other detectors. This involves the fact that the ITS geometry
453     // still has several options to be followed in parallel in order to determine
454     // the best set-up which minimizes the induced background. All the geometries
455     // available to date are described in the following. Read carefully the comments
456     // and use the default version (the only one uncommented) unless you are making
457     // comparisons and you know what you are doing. In this case just uncomment the
458     // ITS geometry you want to use and run Aliroot.
459     //
460     // Detailed geometries:         
461     //
462     //
463     //AliITS *ITS  = new AliITSv5symm("ITS","Updated ITS TDR detailed version with symmetric services");
464     //
465     //AliITS *ITS  = new AliITSv5asymm("ITS","Updates ITS TDR detailed version with asymmetric services");
466     //
467         AliITSvPPRasymmFMD *ITS  = new AliITSvPPRasymmFMD("ITS","New ITS PPR detailed version with asymmetric services");
468         ITS->SetMinorVersion(2);  // don't touch this parameter if you're not an ITS developer
469         ITS->SetReadDet(kTRUE);   // don't touch this parameter if you're not an ITS developer
470     //    ITS->SetWriteDet("$ALICE_ROOT/ITS/ITSgeometry_vPPRasymm2.det");  // don't touch this parameter if you're not an ITS developer
471         ITS->SetThicknessDet1(200.);   // detector thickness on layer 1 must be in the range [100,300]
472         ITS->SetThicknessDet2(200.);   // detector thickness on layer 2 must be in the range [100,300]
473         ITS->SetThicknessChip1(200.);  // chip thickness on layer 1 must be in the range [150,300]
474         ITS->SetThicknessChip2(200.);  // chip thickness on layer 2 must be in the range [150,300]
475         ITS->SetRails(0);            // 1 --> rails in ; 0 --> rails out
476         ITS->SetCoolingFluid(1);   // 1 --> water ; 0 --> freon
477
478     // Coarse geometries (warning: no hits are produced with these coarse geometries and they unuseful 
479     // for reconstruction !):
480     //                                                     
481     //
482     //AliITSvPPRcoarseasymm *ITS  = new AliITSvPPRcoarseasymm("ITS","New ITS PPR coarse version with asymmetric services");
483     //ITS->SetRails(0);                // 1 --> rails in ; 0 --> rails out
484     //ITS->SetSupportMaterial(0);      // 0 --> Copper ; 1 --> Aluminum ; 2 --> Carbon
485     //
486     //AliITS *ITS  = new AliITSvPPRcoarsesymm("ITS","New ITS PPR coarse version with symmetric services");
487     //ITS->SetRails(0);                // 1 --> rails in ; 0 --> rails out
488     //ITS->SetSupportMaterial(0);      // 0 --> Copper ; 1 --> Aluminum ; 2 --> Carbon
489     //                      
490     //
491     //
492     // Geant3 <-> EUCLID conversion
493     // ============================
494     //
495     // SetEUCLID is a flag to output (=1) or not to output (=0) both geometry and
496     // media to two ASCII files (called by default ITSgeometry.euc and
497     // ITSgeometry.tme) in a format understandable to the CAD system EUCLID.
498     // The default (=0) means that you dont want to use this facility.
499     //
500         ITS->SetEUCLID(0);  
501     }
502
503     if (iTPC)
504     {
505       //============================ TPC parameters =====================
506         AliTPC *TPC = new AliTPCv2("TPC", "Default");
507     }
508
509
510     if (iTOF) {
511         //=================== TOF parameters ============================
512         AliTOF *TOF = new AliTOFv5T0("TOF", "normal TOF");
513     }
514
515
516     if (iHMPID)
517     {
518         //=================== HMPID parameters ===========================
519         AliHMPID *HMPID = new AliHMPIDv1("HMPID", "normal HMPID");
520
521     }
522
523
524     if (iZDC)
525     {
526         //=================== ZDC parameters ============================
527
528         AliZDC *ZDC = new AliZDCv2("ZDC", "normal ZDC");
529     }
530
531     if (iTRD)
532     {
533         //=================== TRD parameters ============================
534
535         AliTRD *TRD = new AliTRDv1("TRD", "TRD slow simulator");
536         AliTRDgeometry *geoTRD = TRD->GetGeometry();
537  
538     }
539
540     if (iFMD)
541     {
542         //=================== FMD parameters ============================
543         AliFMD *FMD = new AliFMDv1("FMD", "normal FMD");
544    }
545
546     if (iMUON)
547     {
548         //=================== MUON parameters ===========================
549         // New MUONv1 version (geometry defined via builders)
550         AliMUON *MUON = new AliMUONv1("MUON", "default");
551     }
552     //=================== PHOS parameters ===========================
553
554     if (iPHOS)
555     {
556         AliPHOS *PHOS = new AliPHOSv1("PHOS", "IHEP");
557     }
558
559
560     if (iPMD)
561     {
562         //=================== PMD parameters ============================
563         AliPMD *PMD = new AliPMDv1("PMD", "normal PMD");
564     }
565
566     if (iT0)
567     {
568         //=================== T0 parameters ============================
569         AliT0 *T0 = new AliT0v1("T0", "T0 Detector");
570     }
571
572     if (iEMCAL)
573     {
574         //=================== EMCAL parameters ============================
575         AliEMCAL *EMCAL = new AliEMCALv2("EMCAL", "SHISH");
576     }
577
578      if (iACORDE)
579     {
580         //=================== ACORDE parameters ============================
581         AliACORDE *ACORDE = new AliACORDEv1("ACORDE", "normal ACORDE");
582     }
583
584      if (iVZERO)
585     {
586         //=================== ACORDE parameters ============================
587         AliVZERO *VZERO = new AliVZEROv6("VZERO", "normal VZERO");
588     }
589 }
590
591
592 void ProcessEnvironmentVars()
593 {
594     // Run type
595     if (gSystem->Getenv("CONFIG_RUN_TYPE")) {
596       for (Int_t iRun = 0; iRun < kRunMax; iRun++) {
597         if (strcmp(gSystem->Getenv("CONFIG_RUN_TYPE"), pprRunName[iRun])==0) {
598           proc = (PDC06Proc_t)iRun;
599           cout<<"Run type set to "<<pprRunName[iRun]<<endl;
600         }
601       }
602     }
603
604     // Random Number seed
605     if (gSystem->Getenv("CONFIG_SEED")) {
606       seed = atoi(gSystem->Getenv("CONFIG_SEED"));
607     }
608 }
609
610
611