efficiency and pt Resolution v1 (M. Verweij)
[u/mrichter/AliRoot.git] / MFT / Config.C
CommitLineData
c95698be 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
43enum PDCProc_t {kGenBox,
44 kGenMuonLMR,
6e8e4c4a 45 kGenPionKaon,
c95698be 46 kGenCorrHF,
47 kPythia6,
48 kPythiaPerugia0,
49 kPythiaPerugia0Jpsi2e,
50 kPythiaPerugia0BtoJpsi2e,
51 kHijing,
52 kHijing2500,
53 kHijing2500Cocktail};
54
55const Char_t* pprRunName[] = {"kGenBox",
56 "kGenMuonLMR",
6e8e4c4a 57 "kGenPionKaon",
c95698be 58 "kGenCorrHF",
59 "kPythia6",
60 "kPythiaPerugia0",
61 "kPythiaPerugia0Jpsi2e",
62 "kPythiaPerugia0BtoJpsi2e",
63 "kHijing",
64 "kHijing2500",
65 "kHijing2500Cocktail"};
66
67enum Mag_t { kNoField, k5kG, kFieldMax };
68
69const Char_t* pprField[] = { "kNoField", "k5kG", "kFieldMax" };
70
71void LoadLibs();
72
73// ----------------------- Generator, field, beam energy,... ------------------------------------------------------------
3a11a1cf 74static PDCProc_t proc = kGenBox;
c95698be 75static PDCProc_t signal = kGenBox; // only in case kHijing2500Cocktail is the proc
76static Mag_t mag = k5kG;
77static Float_t energy = 14000.; // energy in CMS
78static Float_t bMin = 0.;
79static 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)
80static Double_t JpsiPol = 0; // Jpsi polarisation
81static Bool_t JpsiHarderPt = kFALSE; // Jpsi harder pt spectrum (8.8 TeV)
82// ----------------------------------------------------------------------------------------------------------------------
83
84static TString comment;
85
86//====================================================================================================================================================
87
88void 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();
6e8e4c4a 159 else if (proc == kGenPionKaon) gener = GenParamPionKaon();
c95698be 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");
6e8e4c4a 231 AliPIPE *PIPE = new AliPIPEv4("PIPE", "Beam Pipe", 1.98, 0.08);
c95698be 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
263AliGenerator* 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
278AliGenerator* 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
6e8e4c4a 300AliGenerator* 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
c95698be 316AliGenerator* 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
335AliGenerator* 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
355AliGenerator* 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
377AliGenerator* 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
393AliGenerator* 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
424AliGenerator* 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
440AliGenerator* 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
466AliGenerator* Hijing2500() {
467
468 AliGenHijing *gener = (AliGenHijing*) Hijing();
469 gener->SetJetQuenching(0);
470 gener->SetPtHardMin (4.5);
471 return gener;
472
473}
474
475//====================================================================================================================================================
476
477AliGenerator* 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
516void 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