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