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