Hopefully now Coverity will be happy
[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,
45 kGenCorrHF,
46 kPythia6,
47 kPythiaPerugia0,
48 kPythiaPerugia0Jpsi2e,
49 kPythiaPerugia0BtoJpsi2e,
50 kHijing,
51 kHijing2500,
52 kHijing2500Cocktail};
53
54const Char_t* pprRunName[] = {"kGenBox",
55 "kGenMuonLMR",
56 "kGenCorrHF",
57 "kPythia6",
58 "kPythiaPerugia0",
59 "kPythiaPerugia0Jpsi2e",
60 "kPythiaPerugia0BtoJpsi2e",
61 "kHijing",
62 "kHijing2500",
63 "kHijing2500Cocktail"};
64
65enum Mag_t { kNoField, k5kG, kFieldMax };
66
67const Char_t* pprField[] = { "kNoField", "k5kG", "kFieldMax" };
68
69void LoadLibs();
70
71// ----------------------- Generator, field, beam energy,... ------------------------------------------------------------
72static PDCProc_t proc = kGenBox;
73static PDCProc_t signal = kGenBox; // only in case kHijing2500Cocktail is the proc
74static Mag_t mag = k5kG;
75static Float_t energy = 14000.; // energy in CMS
76static Float_t bMin = 0.;
77static 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)
78static Double_t JpsiPol = 0; // Jpsi polarisation
79static Bool_t JpsiHarderPt = kFALSE; // Jpsi harder pt spectrum (8.8 TeV)
80// ----------------------------------------------------------------------------------------------------------------------
81
82static TString comment;
83
84//====================================================================================================================================================
85
86void 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");
f43b7000 228 AliPIPE *PIPE = new AliPIPEv4("PIPE", "Beam Pipe", 1.98, 0.08); // cylindre with new adaptator pipe
c95698be 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
260AliGenerator* 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
275AliGenerator* 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
297AliGenerator* 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
316AliGenerator* 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
336AliGenerator* 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
358AliGenerator* 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
374AliGenerator* 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
405AliGenerator* 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
421AliGenerator* 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
447AliGenerator* Hijing2500() {
448
449 AliGenHijing *gener = (AliGenHijing*) Hijing();
450 gener->SetJetQuenching(0);
451 gener->SetPtHardMin (4.5);
452 return gener;
453
454}
455
456//====================================================================================================================================================
457
458AliGenerator* 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
497void 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