Support files added
[u/mrichter/AliRoot.git] / MFT / Config.C
1 #if !defined(__CINT__) || defined(__MAKECINT__)
2 #include <Riostream.h>
3 #include <TRandom.h>
4 #include <TDatime.h>
5 #include <TSystem.h>
6 #include <TVirtualMC.h>
7 #include <TGeant3TGeo.h>
8 #include "AliRunLoader.h"
9 #include "AliRun.h"
10 #include "AliConfig.h"
11 #include "AliDecayerPythia.h"
12 #include "AliGenPythia.h"
13 #include "AliGenDPMjet.h"
14 #include "AliMagFCheb.h"
15 #include "AliBODY.h"
16 #include "AliMAG.h"
17 #include "AliABSOv3.h"
18 #include "AliDIPOv3.h"
19 #include "AliHALLv3.h"
20 #include "AliFRAMEv2.h"
21 #include "AliSHILv3.h"
22 #include "AliPIPEv3.h"
23 #include "AliPIPEv4.h"
24 #include "AliITSv11.h"
25 #include "AliTPCv2.h"
26 #include "AliTOFv6T0.h"
27 #include "AliHMPIDv3.h"
28 #include "AliZDCv3.h"
29 #include "AliTRDv1.h"
30 #include "AliTRDgeometry.h"
31 #include "AliFMDv1.h"
32 #include "AliMUONv1.h"
33 #include "AliPHOSv1.h"
34 #include "AliPHOSSimParam.h"
35 #include "AliPMDv1.h"
36 #include "AliT0v1.h"
37 #include "AliEMCALv2.h"
38 #include "AliACORDEv1.h"
39 #include "AliVZEROv7.h"
40 #include "AliMFT.h"
41 #endif
42
43 enum PDCProc_t {kGenBox,
44                 kGenMuonLMR,
45                 kGenCorrHF,
46                 kPythia6,
47                 kPythiaPerugia0, 
48                 kPythiaPerugia0Jpsi2e, 
49                 kPythiaPerugia0BtoJpsi2e, 
50                 kHijing, 
51                 kHijing2500,
52                 kHijing2500Cocktail};
53
54 const Char_t* pprRunName[] = {"kGenBox",
55                               "kGenMuonLMR",
56                               "kGenCorrHF",
57                               "kPythia6",
58                               "kPythiaPerugia0", 
59                               "kPythiaPerugia0Jpsi2e", 
60                               "kPythiaPerugia0BtoJpsi2e", 
61                               "kHijing", 
62                               "kHijing2500", 
63                               "kHijing2500Cocktail"};
64
65 enum Mag_t { kNoField, k5kG, kFieldMax };
66
67 const Char_t* pprField[] = { "kNoField", "k5kG", "kFieldMax" };
68
69 void LoadLibs();
70
71 // ----------------------- Generator, field, beam energy,... ------------------------------------------------------------
72 static PDCProc_t     proc     = kGenBox;
73 static PDCProc_t     signal   = kGenBox;    // only in case kHijing2500Cocktail is the proc
74 static Mag_t         mag      = k5kG;
75 static Float_t       energy   = 14000.; // energy in CMS
76 static Float_t       bMin     = 0.;
77 static Float_t       bMax =   = 5.; // 0-5 fm corresponds to around 0-10% (see https://twiki.cern.ch/twiki/bin/viewauth/ALICE/CentStudies#Tables_with_centrality_bins_for)
78 static Double_t      JpsiPol  = 0; // Jpsi polarisation
79 static Bool_t        JpsiHarderPt = kFALSE; // Jpsi harder pt spectrum (8.8 TeV)
80 // ----------------------------------------------------------------------------------------------------------------------
81
82 static TString comment;
83
84 //====================================================================================================================================================
85
86 void Config() {
87
88   //  AliLog::SetClassDebugLevel("AliMFT", 1);
89
90   LoadLibs();
91
92   new TGeant3TGeo("C++ Interface to Geant3");
93
94   // Create the output file
95
96   AliRunLoader* rl=0x0;
97
98   cout<<"Config.C: Creating Run Loader ..."<<endl;
99   rl = AliRunLoader::Open("galice.root", AliConfig::GetDefaultEventFolderName(), "recreate");
100   if (rl == 0x0) {
101     gAlice->Fatal("Config.C","Can not instatiate the Run Loader");
102     return;
103   }
104   rl->SetCompressionLevel(2);
105   rl->SetNumberOfEventsPerFile(1000);
106   gAlice->SetRunLoader(rl);
107
108   // ************* STEERING parameters FOR ALICE SIMULATION **************
109   // --- Specify event type to be tracked through the ALICE setup
110   // --- All positions are in cm, angles in degrees, and P and E in GeV
111   
112   gMC->SetProcess("DCAY",1);
113   gMC->SetProcess("PAIR",1);
114   gMC->SetProcess("COMP",1);
115   gMC->SetProcess("PHOT",1);
116   gMC->SetProcess("PFIS",0);
117   gMC->SetProcess("DRAY",0);
118   gMC->SetProcess("ANNI",1);
119   gMC->SetProcess("BREM",1);
120   gMC->SetProcess("MUNU",1);
121   gMC->SetProcess("CKOV",1);
122   gMC->SetProcess("HADR",1);
123   gMC->SetProcess("LOSS",2);
124   gMC->SetProcess("MULS",1);
125   gMC->SetProcess("RAYL",1);
126   
127   Float_t cut = 1.e-3;        // 1MeV cut by default
128   Float_t tofmax = 1.e10;
129   
130   gMC->SetCut("CUTGAM", cut);
131   gMC->SetCut("CUTELE", cut);
132   gMC->SetCut("CUTNEU", cut);
133   gMC->SetCut("CUTHAD", cut);
134   gMC->SetCut("CUTMUO", cut);
135   gMC->SetCut("BCUTE",  cut); 
136   gMC->SetCut("BCUTM",  cut); 
137   gMC->SetCut("DCUTE",  cut); 
138   gMC->SetCut("DCUTM",  cut); 
139   gMC->SetCut("PPCUTM", cut);
140   gMC->SetCut("TOFMAX", tofmax); 
141   
142   TVirtualMCDecayer *decayer = new AliDecayerPythia();
143   decayer->SetForceDecay(kAll);
144   decayer->Init();
145   gMC->SetExternalDecayer(decayer);
146   
147   // Generator
148   AliGenerator* gener = 0x0;
149   if (proc == kPythia6)                   gener = MbPythia(); 
150   else if (proc == kPythiaPerugia0)       gener = MbPythiaTunePerugia0();
151   else if (proc == kHijing)               gener = Hijing();     
152   else if (proc == kHijing2500)           gener = Hijing2500(); 
153   else if (proc == kHijing2500Cocktail)   gener = Hijing2500Cocktail();
154   else if (proc == kGenBox)               gener = GenBox();
155   else if (proc == kGenMuonLMR)           gener = GenMuonLMR();
156   else if (proc == kGenCorrHF)            gener = GenCorrHF();
157
158   // Size of the interaction diamond
159   Float_t sigmaz  = 5.4 / TMath::Sqrt(2.);     // [cm]
160   Float_t betast  = 3.5;                       // beta* [m]
161   Float_t eps     = 3.75e-6;                   // emittance [m]
162   Float_t gamma   = energy / 2.0 / 0.938272;   // relativistic gamma [1]
163   Float_t sigmaxy = TMath::Sqrt(eps * betast / gamma) / TMath::Sqrt(2.) * 100.;  // [cm]
164
165   printf("\n \n Diamond size x-y: %10.3e z: %10.3e\n \n", sigmaxy, sigmaz);
166     
167 //   gener->SetSigma(sigmaxy, sigmaxy, sigmaz);      // Sigma in (X,Y,Z) (cm) on IP position
168 //   gener->SetVertexSmear(kPerEvent);
169   gener->Init();
170
171   printf("\n \n Comment: %s \n \n", comment.Data());
172
173   //  TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", -1., -1., AliMagF::k5kG, AliMagF::kBeamTypeAA, 2750.));
174   TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", -1., -1., AliMagF::k5kG, AliMagF::kBeamTypepp, 7000.));
175
176   rl->CdGAFile();
177   
178   // Detector Setup
179
180   Int_t iABSO  = 1;
181   Int_t iDIPO  = 1;
182   Int_t iHALL  = 1;
183   Int_t iMUON  = 1;
184   Int_t iPIPE  = 1;
185   Int_t iSHIL  = 1;
186   Int_t iT0    = 0;
187   Int_t iVZERO = 1;
188   Int_t iMFT   = 1;
189   Int_t iACORDE= 0;
190   Int_t iEMCAL = 0;
191   Int_t iFMD   = 0;
192   Int_t iFRAME = 0;
193   Int_t iITS   = 0;
194   Int_t iMAG   = 0;
195   Int_t iPHOS  = 0;
196   Int_t iPMD   = 0;
197   Int_t iHMPID = 0;
198   Int_t iTOF   = 0;
199   Int_t iTPC   = 0;
200   Int_t iTRD   = 0;
201   Int_t iZDC   = 0;
202   
203
204   AliBODY *BODY = new AliBODY("BODY", "Alice envelop");
205
206   if (iMAG)       AliMAG    *MAG    = new AliMAG("MAG", "Magnet");
207   if (iABSO)      AliABSO   *ABSO   = new AliABSOv3("ABSO", "Muon Absorber");
208   if (iDIPO)      AliDIPO   *DIPO   = new AliDIPOv3("DIPO", "Dipole version 3");
209   if (iHALL)      AliHALL   *HALL   = new AliHALLv3("HALL", "Alice Hall");
210   if (iSHIL)      AliSHIL   *SHIL   = new AliSHILv3("SHIL", "Shielding Version 3");
211   if (iITS)       AliITS    *ITS    = new AliITSv11("ITS","ITS v11");
212   if (iTPC)       AliTPC    *TPC    = new AliTPCv2("TPC", "Default");
213   if (iTOF)       AliTOF    *TOF    = new AliTOFv6T0("TOF", "normal TOF");
214   if (iHMPID)     AliHMPID  *HMPID  = new AliHMPIDv3("HMPID", "normal HMPID");
215   if (iFMD)       AliFMD    *FMD    = new AliFMDv1("FMD", "normal FMD");
216   if (iPHOS)      AliPHOS   *PHOS   = new AliPHOSv1("PHOS", "noCPV_Modules123");
217   if (iPMD)       AliPMD    *PMD    = new AliPMDv1("PMD", "normal PMD");
218   if (iT0)        AliT0     *T0     = new AliT0v1("T0", "T0 Detector");
219   if (iEMCAL)     AliEMCAL  *EMCAL  = new AliEMCALv2("EMCAL", "EMCAL_FIRSTYEARV1");
220   if (iACORDE)    AliACORDE *ACORDE = new AliACORDEv1("ACORDE", "normal ACORDE");
221   if (iVZERO)     AliVZERO  *VZERO  = new AliVZEROv7("VZERO", "normal VZERO");
222   if (iFRAME) {
223     AliFRAMEv2 *FRAME = new AliFRAMEv2("FRAME", "Space Frame");
224     FRAME->SetHoles(1);
225   }
226   if (iPIPE) {
227     //    AliPIPE *PIPE = new AliPIPEv3("PIPE", "Beam Pipe");
228     AliPIPE *PIPE = new AliPIPEv4("PIPE", "Beam Pipe",0.,1.98,0.08,8.0,-50., 3, -57.4, 2.0, 4.0);  // cylindre with new adaptator pipe starting at -59.7cm
229   }
230   if (iZDC) {
231     AliZDC *ZDC = new AliZDCv3("ZDC", "normal ZDC");
232     ZDC->SetSpectatorsTrack();  
233     ZDC->SetLumiLength(0.);
234   }
235   if (iTRD) {
236     AliTRD *TRD = new AliTRDv1("TRD", "TRD slow simulator");
237   }
238   if (iMUON) {
239     AliMUON *MUON = new AliMUONv1("MUON", "default");
240     MUON->SetTriggerEffCells(1); // not needed if raw masks 
241     Char_t* digitstore="AliMUONDigitStoreV2S";    
242     MUON->SetDigitStoreClassName(digitstore);
243   }
244   if (iMFT) {
245     AliMFT *MFT = new AliMFT("MFT", "normal MFT");
246     MFT->SetNSlices(1);
247     MFT->SetChargeDispersion(20.e-4);
248     MFT->SetNStepForChargeDispersion(4);
249   }
250
251   TIter next(gAlice->Modules());
252   AliModule *detector;
253   printf("gAlice->Modules:\n");
254   while((detector = (AliModule*)next())) printf("%s\n",detector->GetName());
255
256 }
257
258 //====================================================================================================================================================
259
260 AliGenerator* GenBox() {
261
262   AliGenBox *gener = new AliGenBox(2);
263   gener->SetMomentumRange(10.0, 10.1);
264   gener->SetPhiRange(0., 360.);         
265   gener->SetThetaRange(171.0,178.0);
266   gener->SetPart(kMuonMinus);           // Muons
267   gener->SetOrigin(0., 0., 0.);         // vertex position
268   
269   return gener;
270   
271 }
272
273 //====================================================================================================================================================
274
275 AliGenerator* GenMuonLMR() {
276   
277   AliGenMUONLMR *gener = new AliGenMUONLMR();
278   gener->SetMomentumRange(0,999);
279   gener->SetPtRange(0,100.);
280   gener->SetPhiRange(0., 360.);
281   gener->SetThetaRange(0., 180.);
282   gener->SetChildThetaRange(170.0,179.0);
283   gener->SetYRange(-4.5, -2.0);
284   gener->SetOrigin(0.0, 0.0, 0.0);             // vertex position
285   gener->SetSigma(0.0, 0.0, 0.0);              // vertex position smearing
286   gener->SetVertexSmear(kPerEvent);
287   enum {kEta2Body, kEtaDalitz, kRho2Body, kOmega2Body, kOmegaDalitz, kPhi2Body, kEtaPrimeDalitz, kPionLMR, kKaonLMR}; 
288   gener->GenerateSingleProcess(kPhi2Body, 500);
289   gener->SetCutOnChild(1);
290
291   return gener;
292
293 }
294
295 //====================================================================================================================================================
296
297 AliGenerator* GenCorrHF() {
298   
299   AliGenCorrHF *gener = new AliGenCorrHF(1, 4, 14);  // for charm, 1 pair per event
300   // AliGenCorrHF *gener = new AliGenCorrHF(1, 5, 14);  // for beauty, 1 pair per event
301   
302   gener->SetMomentumRange(0,9999);
303   gener->SetCutOnChild(1);          // 1/0 means cuts on children enable/disable
304   gener->SetChildThetaRange(170.0,179.0);
305   gener->SetOrigin(0,0,0);          //vertex position    
306   gener->SetForceDecay(kSemiMuonic);
307   gener->SetTrackingFlag(1);
308   gener->Init();
309   
310   return gener;
311   
312 }
313
314 //====================================================================================================================================================
315
316 AliGenerator* MbPythia() {
317   
318   comment = comment.Append(" pp: Pythia low-pt");
319   
320   //    Pythia
321   AliGenPythia* pythia = new AliGenPythia(-1); 
322   pythia->SetMomentumRange(0, 999999.);
323   //  pythia->SetThetaRange(0., 180.);
324   //  pythia->SetChildYRange(-12.,0.);
325   //  pythia->SetPtRange(0,1000.);
326   //  pythia->SetCutOnChild(1);
327   pythia->SetProcess(kPyMb);
328   pythia->SetEnergyCMS(energy);
329   pythia->SetForceDecay(kSemiMuonic);
330   
331   return pythia;
332 }
333
334 //====================================================================================================================================================
335
336 AliGenerator* MbPythiaTunePerugia0() {
337   
338   comment = comment.Append(" pp: Pythia low-pt (Perugia0)");
339   
340   //    Pythia
341   AliGenPythia* pythia = new AliGenPythia(-1); 
342   pythia->SetMomentumRange(0, 999999.);
343   pythia->SetThetaRange(0., 180.);
344   pythia->SetYRange(-12.,12.);
345   pythia->SetPtRange(0,1000.);
346   pythia->SetProcess(kPyMb);
347   pythia->SetEnergyCMS(energy);
348   //    Tune
349   //    320     Perugia 0
350   pythia->SetTune(320); 
351   pythia->UseNewMultipleInteractionsScenario();
352   
353   return pythia;
354 }
355
356 //====================================================================================================================================================
357
358 AliGenerator* MbPythiaTunePerugia0Jpsi2e() {
359   
360   comment = comment.Append("Jpsi forced to dielectrons");
361   AliGenParam *jpsi=0x0;
362   if(JpsiHarderPt) jpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp 8.8", "Jpsi");  // 8.8 TeV
363   else jpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp 7", "Jpsi");  // 7 TeV
364   jpsi->SetPtRange(0.,999.);
365   jpsi->SetYRange(-1.0, 1.0);
366   jpsi->SetPhiRange(0.,360.);
367   jpsi->SetForceDecay(kDiElectron);
368   return jpsi;
369
370 }
371
372 //====================================================================================================================================================
373
374 AliGenerator* MbPythiaTunePerugia0BtoJpsi2e() {
375
376   comment = comment.Append(" pp: Pythia (Perugia0) BtoJpsi (1 bbbar per event, 1 b-hadron in |y|<2, 1 J/psi in |y|<2");
377   
378   //    Pythia
379   AliGenPythia* pythia = new AliGenPythia(-1);
380   pythia->SetMomentumRange(0, 999999.);
381   pythia->SetThetaRange(0., 180.);
382   pythia->SetYRange(-2.,2.);
383   pythia->SetPtRange(0,1000.);
384   pythia->SetProcess(kPyBeautyppMNRwmi);
385   pythia->SetEnergyCMS(energy);
386   //    Tune
387   //    320     Perugia 0
388   pythia->SetTune(320);
389   pythia->UseNewMultipleInteractionsScenario();
390   //
391   //    decays
392   pythia->SetCutOnChild(1);
393   pythia->SetPdgCodeParticleforAcceptanceCut(443);
394   pythia->SetChildYRange(-2,2);
395   pythia->SetChildPtRange(0,10000.);
396   //
397   //    decays
398   pythia->SetForceDecay(kBJpsiDiElectron);
399   
400   return pythia;
401 }
402
403 //====================================================================================================================================================
404
405 AliGenerator* HijingParam() {
406
407   AliGenHIJINGpara *gener = new AliGenHIJINGpara(2000);
408   gener->SetMomentumRange(0,999);
409   gener->SetPhiRange(0,360);
410   gener->SetThetaRange(171, 179);
411   gener->SetOrigin(0,0,0);        // vertex position
412   gener->SetSigma(0,0,0);         // Sigma in (X,Y,Z) (cm) on IP position
413   gGener = gener;
414   
415   return gener;
416
417 }
418
419 //====================================================================================================================================================
420
421 AliGenerator* Hijing() {
422   
423   AliGenHijing *gener = new AliGenHijing(-1);
424   // centre of mass energy 
425   gener->SetEnergyCMS(energy);
426   gener->SetImpactParameterRange(bMin, bMax);   
427   // reference frame
428   gener->SetReferenceFrame("CMS");
429   // projectile
430   gener->SetProjectile("A", 208, 82);
431   gener->SetTarget    ("A", 208, 82);
432   // tell hijing to keep the full parent child chain
433   gener->KeepFullEvent();
434   // enable jet quenching
435   gener->SetJetQuenching(1);
436   // enable shadowing
437   gener->SetShadowing(1);
438   // Don't track spectators
439   gener->SetSpectators(0);
440   // kinematic selection
441   gener->SetSelectAll(0);
442   return gener;
443 }
444
445 //====================================================================================================================================================
446
447 AliGenerator* Hijing2500() {
448   
449   AliGenHijing *gener = (AliGenHijing*) Hijing();
450   gener->SetJetQuenching(0);    
451   gener->SetPtHardMin (4.5);
452   return gener;
453   
454 }
455
456 //====================================================================================================================================================
457
458 AliGenerator* Hijing2500Cocktail() {
459   
460   comment = comment.Append(" PbPb: Hjing2500 at 5.5 + muon signals");
461
462   AliGenCocktail *cocktail = new AliGenCocktail();
463   cocktail->SetProjectile("A", 208, 82);
464   cocktail->SetTarget    ("A", 208, 82);
465   cocktail->SetEnergyCMS(energy);
466
467   // 1 Hijing event: provides underlying event and collision geometry 
468   AliGenHijing *hijing = Hijing2500();
469   cocktail->AddGenerator(hijing,"hijing",1);
470
471   // generator for the custom signal
472   AliGenerator* signalGen = 0x0;      
473   switch (signal) {
474   case 0:
475     signalGen = GenBox();
476     break;
477   case 1:
478     signalGen = GenMuonLMR(); 
479     break;
480   case 2:
481     signalGen = GenCorrHF();
482     break;
483   default:
484     signalGen = GenBox();
485     break;
486   }
487   cocktail->AddGenerator(signalGen, "signal", 1);
488
489   cocktail->SetTrackingFlag(1);
490
491   return cocktail;
492
493 }
494
495 //====================================================================================================================================================
496
497 void LoadLibs() {
498
499 #if defined(__CINT__)
500
501   gSystem->Load("liblhapdf");      // Parton density functions
502   gSystem->Load("libEGPythia6");   // TGenerator interface
503   if (proc == kPythia6) {
504     gSystem->Load("libpythia6");        // Pythia 6.2
505     gSystem->Load("libAliPythia6");     // ALICE specific implementations
506   } 
507   else {
508     gSystem->Load("libpythia6.4.21");   // Pythia 6.4
509     gSystem->Load("libAliPythia6");     // ALICE specific implementations       
510   }
511   
512   if (proc == kHijing || proc == kHijing2500 || proc == kHijing2500Cocktail) {
513     gSystem->Load("libhijing"); 
514     gSystem->Load("libTHijing");
515   } 
516   
517   gSystem->Load("libgeant321");
518   
519 #endif
520
521 }
522
523 //====================================================================================================================================================
524