1 #if !defined(__CINT__) || defined(__MAKECINT__)
6 #include <TVirtualMC.h>
7 #include <TGeant3TGeo.h>
8 #include "AliRunLoader.h"
10 #include "AliConfig.h"
11 #include "AliDecayerPythia.h"
12 #include "AliGenPythia.h"
13 #include "AliGenDPMjet.h"
14 #include "AliMagFCheb.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"
26 #include "AliTOFv6T0.h"
27 #include "AliHMPIDv3.h"
30 #include "AliTRDgeometry.h"
32 #include "AliMUONv1.h"
33 #include "AliPHOSv1.h"
34 #include "AliPHOSSimParam.h"
37 #include "AliEMCALv2.h"
38 #include "AliACORDEv1.h"
39 #include "AliVZEROv7.h"
43 // Git push Test: R. Tieulent
45 enum PDCProc_t {kGenBox,
52 kPythiaPerugia0Jpsi2mu,
53 kPythiaPerugia0BtoJpsi2mu,
58 const Char_t* pprRunName[] = {"kGenBox",
65 "kPythiaPerugia0Jpsi2mu",
66 "kPythiaPerugia0BtoJpsi2mu",
69 "kHijing2500Cocktail"};
71 enum Mag_t { kNoField, k5kG, kFieldMax };
73 const Char_t* pprField[] = { "kNoField", "k5kG", "kFieldMax" };
77 // ----------------------- Generator, field, beam energy,... ------------------------------------------------------------
78 static PDCProc_t proc = kGenBox;
79 static PDCProc_t signal = kGenBox; // only in case kHijing2500Cocktail is the proc
80 static Mag_t mag = k5kG;
81 static Float_t energy = 5500.; // energy in CMS
82 static Float_t bMin = 0.;
83 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)
84 static Double_t JpsiPol = 0; // Jpsi polarisation
85 static Bool_t JpsiHarderPt = kFALSE; // Jpsi harder pt spectrum (8.8 TeV)
86 // ----------------------------------------------------------------------------------------------------------------------
88 static TString comment;
90 //====================================================================================================================================================
94 // AliLog::SetClassDebugLevel("AliMFT", 1);
98 new TGeant3TGeo("C++ Interface to Geant3");
100 // Create the output file
102 AliRunLoader* rl=0x0;
104 printf("Config.C: Creating Run Loader ...");
105 rl = AliRunLoader::Open("galice.root", AliConfig::GetDefaultEventFolderName(), "recreate");
107 gAlice->Fatal("Config.C","Can not instatiate the Run Loader");
110 rl->SetCompressionLevel(2);
111 rl->SetNumberOfEventsPerFile(1000);
112 gAlice->SetRunLoader(rl);
114 // ************* STEERING parameters FOR ALICE SIMULATION **************
115 // --- Specify event type to be tracked through the ALICE setup
116 // --- All positions are in cm, angles in degrees, and P and E in GeV
118 gMC->SetProcess("DCAY",1);
119 gMC->SetProcess("PAIR",1);
120 gMC->SetProcess("COMP",1);
121 gMC->SetProcess("PHOT",1);
122 gMC->SetProcess("PFIS",0);
123 gMC->SetProcess("DRAY",0);
124 gMC->SetProcess("ANNI",1);
125 gMC->SetProcess("BREM",1);
126 gMC->SetProcess("MUNU",1);
127 gMC->SetProcess("CKOV",1);
128 gMC->SetProcess("HADR",1);
129 gMC->SetProcess("LOSS",2);
130 gMC->SetProcess("MULS",1);
131 gMC->SetProcess("RAYL",1);
133 Float_t cut = 1.e-3; // 1MeV cut by default
134 Float_t tofmax = 1.e10;
136 gMC->SetCut("CUTGAM", cut);
137 gMC->SetCut("CUTELE", cut);
138 gMC->SetCut("CUTNEU", cut);
139 gMC->SetCut("CUTHAD", cut);
140 gMC->SetCut("CUTMUO", cut);
141 gMC->SetCut("BCUTE", cut);
142 gMC->SetCut("BCUTM", cut);
143 gMC->SetCut("DCUTE", cut);
144 gMC->SetCut("DCUTM", cut);
145 gMC->SetCut("PPCUTM", cut);
146 gMC->SetCut("TOFMAX", tofmax);
148 TVirtualMCDecayer *decayer = new AliDecayerPythia();
149 decayer->SetForceDecay(kAll);
151 gMC->SetExternalDecayer(decayer);
154 AliGenerator* gener = 0x0;
155 if (proc == kPythia6) gener = MbPythia();
156 else if (proc == kPythiaPerugia0) gener = MbPythiaTunePerugia0();
157 else if (proc == kHijing) gener = Hijing();
158 else if (proc == kHijing2500) gener = Hijing2500();
159 else if (proc == kHijing2500Cocktail) gener = Hijing2500Cocktail();
160 else if (proc == kGenBox) gener = GenBox();
161 else if (proc == kGenMuonLMR) gener = GenMuonLMR();
162 else if (proc == kGenParamJpsi) gener = GenParamJpsi();
163 else if (proc == kGenCorrHF) gener = GenCorrHF();
164 else if (proc == kGenPionKaon) gener = GenParamPionKaon();
165 else if (proc == kPythiaPerugia0BtoJpsi2mu) gener = MbPythiaTunePerugia0BtoJpsi2mu();
167 // Size of the interaction diamond
168 Float_t sigmaz = 5.4 / TMath::Sqrt(2.); // [cm]
169 Float_t betast = 3.5; // beta* [m]
170 Float_t eps = 3.75e-6; // emittance [m]
171 Float_t gamma = energy / 2.0 / 0.938272; // relativistic gamma [1]
172 Float_t sigmaxy = TMath::Sqrt(eps * betast / gamma) / TMath::Sqrt(2.) * 100.; // [cm]
174 printf("\n \n Diamond size x-y: %10.3e z: %10.3e\n \n", sigmaxy, sigmaz);
176 // gener->SetSigma(sigmaxy, sigmaxy, sigmaz); // Sigma in (X,Y,Z) (cm) on IP position
177 // gener->SetVertexSmear(kPerEvent);
180 printf("\n \n Comment: %s \n \n", comment.Data());
182 // TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", -1., -1., AliMagF::k5kG, AliMagF::kBeamTypeAA, 2750.));
183 TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", -1., -1., AliMagF::k5kG, AliMagF::kBeamTypepp, 7000.));
213 AliBODY *BODY = new AliBODY("BODY", "Alice envelop");
215 if (iMAG) AliMAG *MAG = new AliMAG("MAG", "Magnet");
216 if (iABSO) AliABSO *ABSO = new AliABSOv3("ABSO", "Muon Absorber");
217 if (iDIPO) AliDIPO *DIPO = new AliDIPOv3("DIPO", "Dipole version 3");
218 if (iHALL) AliHALL *HALL = new AliHALLv3("HALL", "Alice Hall");
219 if (iSHIL) AliSHIL *SHIL = new AliSHILv3("SHIL", "Shielding Version 3");
220 if (iITS) AliITS *ITS = new AliITSv11("ITS","ITS v11");
221 if (iTPC) AliTPC *TPC = new AliTPCv2("TPC", "Default");
222 if (iTOF) AliTOF *TOF = new AliTOFv6T0("TOF", "normal TOF");
223 if (iHMPID) AliHMPID *HMPID = new AliHMPIDv3("HMPID", "normal HMPID");
224 if (iFMD) AliFMD *FMD = new AliFMDv1("FMD", "normal FMD");
225 if (iPHOS) AliPHOS *PHOS = new AliPHOSv1("PHOS", "noCPV_Modules123");
226 if (iPMD) AliPMD *PMD = new AliPMDv1("PMD", "normal PMD");
227 if (iT0) AliT0 *T0 = new AliT0v1("T0", "T0 Detector");
228 if (iEMCAL) AliEMCAL *EMCAL = new AliEMCALv2("EMCAL", "EMCAL_FIRSTYEARV1");
229 if (iACORDE) AliACORDE *ACORDE = new AliACORDEv1("ACORDE", "normal ACORDE");
230 if (iVZERO) AliVZERO *VZERO = new AliVZEROv7("VZERO", "normal VZERO");
232 AliFRAMEv2 *FRAME = new AliFRAMEv2("FRAME", "Space Frame");
236 // AliPIPE *PIPE = new AliPIPEv3("PIPE", "Beam Pipe");
237 AliPIPE *PIPE = new AliPIPEv4("PIPE", "Beam Pipe");
240 AliZDC *ZDC = new AliZDCv3("ZDC", "normal ZDC");
241 ZDC->SetSpectatorsTrack();
242 ZDC->SetLumiLength(0.);
245 AliTRD *TRD = new AliTRDv1("TRD", "TRD slow simulator");
248 AliMUON *MUON = new AliMUONv1("MUON", "default");
249 MUON->SetTriggerEffCells(1); // not needed if raw masks
250 Char_t* digitstore="AliMUONDigitStoreV2S";
251 MUON->SetDigitStoreClassName(digitstore);
254 AliMFT *MFT = new AliMFT("MFT", "normal MFT");
257 TIter next(gAlice->Modules());
259 printf("gAlice->Modules:\n");
260 while((detector = (AliModule*)next())) printf("%s\n",detector->GetName());
264 //====================================================================================================================================================
266 AliGenerator* GenBox() {
268 AliGenBox *gener = new AliGenBox(2);
269 gener->SetMomentumRange(10.0, 10.1);
270 gener->SetPhiRange(0., 360.);
271 gener->SetThetaRange(171.0,178.0);
272 gener->SetPart(kMuonMinus); // Muons
273 gener->SetOrigin(0., 0., 0.); // vertex position
279 //====================================================================================================================================================
281 AliGenerator* GenMuonLMR() {
283 AliGenMUONLMR *gener = new AliGenMUONLMR();
284 gener->SetMomentumRange(0,999);
285 gener->SetPtRange(0,100.);
286 gener->SetYRange(-4.5, -2.0);
287 gener->SetChildThetaRange(171.0,178.0);
288 gener->SetChildMomentumRange(4.0, 999.);
289 gener->SetOrigin(0.0, 0.0, 0.0); // vertex position
290 gener->SetSigma(0.0, 0.0, 0.0); // vertex position smearing
291 enum {kEta2Body, kEtaDalitz, kRho2Body, kOmega2Body, kOmegaDalitz, kPhi2Body, kEtaPrimeDalitz, kPionLMR, kKaonLMR};
292 gener->GenerateSingleProcess(kOmega2Body, 10);
293 gener->SetCutOnChild(1);
299 //====================================================================================================================================================
301 AliGenerator* GenParamJpsi() {
303 AliGenParam *gener = new AliGenParam(5, AliGenMUONlib::kJpsi);
304 gener->SetMomentumRange(0,999);
305 gener->SetPtRange(0,100.);
306 gener->SetYRange(-4.0, -2.5);
307 gener->SetPhiRange(0., 360.);
308 gener->SetChildThetaRange(171.0,177.0);
309 gener->SetChildMomentumRange(5.0, 999.);
310 gener->SetOrigin(0.0, 0.0, 0.0); // vertex position
311 gener->SetSigma(0.0, 0.0, 0.0); // Sigma in (X,Y,Z) (cm) on IP position
312 gener->SetForceDecay(kDiMuon);
313 gener->SetTrackingFlag(1);
314 gener->SetCutOnChild(1);
320 //====================================================================================================================================================
322 AliGenerator* GenParamPionKaon() {
324 AliGenParamPionsKaons *gener = new AliGenParamPionsKaons(100,"$ALICE_ROOT/MFT/PionKaonKinematics.root");
325 gener->SetPtRange(0, 5.);
326 gener->SetPhiRange(0., 360.);
327 gener->SetYRange(-10., 0.);
328 gener->SetOrigin(0.0, 0.0, 0.0); // vertex position
329 gener->SetSigma(0.0, 0.0, 0.0); // vertex position smearing
330 // gener->SetCutOnChild(1);
336 //====================================================================================================================================================
338 AliGenerator* GenCorrHF() {
340 AliGenCorrHF *gener = new AliGenCorrHF(1, 4, 6); // for charm, 1 pair per event
341 // AliGenCorrHF *gener = new AliGenCorrHF(1, 5, 6); // for beauty, 1 pair per event
343 gener->SetMomentumRange(0,9999);
344 gener->SetCutOnChild(1); // 1/0 means cuts on children enable/disable
345 gener->SetChildThetaRange(171.0,178.0);
346 gener->SetChildMomentumRange(4.0, 999.);
347 gener->SetOrigin(0,0,0); //vertex position
348 gener->SetForceDecay(kSemiMuonic);
349 gener->SetTrackingFlag(1);
356 //====================================================================================================================================================
358 AliGenerator* MbPythia() {
360 comment = comment.Append(" pp: Pythia low-pt");
363 AliGenPythia* pythia = new AliGenPythia(-1);
364 pythia->SetMomentumRange(0, 999999.);
365 // pythia->SetThetaRange(0., 180.);
366 // pythia->SetChildYRange(-12.,0.);
367 // pythia->SetPtRange(0,1000.);
368 // pythia->SetCutOnChild(1);
369 pythia->SetProcess(kPyMb);
370 pythia->SetEnergyCMS(energy);
371 pythia->SetForceDecay(kSemiMuonic);
376 //====================================================================================================================================================
378 AliGenerator* MbPythiaTunePerugia0() {
380 comment = comment.Append(" pp: Pythia low-pt (Perugia0)");
383 AliGenPythia* pythia = new AliGenPythia(-1);
384 pythia->SetMomentumRange(0, 999999.);
385 pythia->SetThetaRange(0., 180.);
386 pythia->SetYRange(-12.,12.);
387 pythia->SetPtRange(0,1000.);
388 pythia->SetProcess(kPyMb);
389 pythia->SetEnergyCMS(energy);
392 pythia->SetTune(320);
393 pythia->UseNewMultipleInteractionsScenario();
398 //====================================================================================================================================================
400 AliGenerator* MbPythiaTunePerugia0Jpsi2mu() {
402 comment = comment.Append("Jpsi forced to dimuons");
403 AliGenParam *jpsi=0x0;
404 if (JpsiHarderPt) jpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp 8.8", "Jpsi"); // 8.8 TeV
405 else jpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp 7", "Jpsi"); // 7 TeV
406 jpsi->SetPtRange(0.,999.);
407 jpsi->SetYRange(-1.0, 1.0);
408 jpsi->SetPhiRange(0.,360.);
409 jpsi->SetOrigin(0,0,0);
410 jpsi->SetForceDecay(kDiMuon);
415 //====================================================================================================================================================
417 AliGenerator* MbPythiaTunePerugia0BtoJpsi2mu() {
419 comment = comment.Append(" pp: Pythia (Perugia0) BtoJpsi (1 bbbar per event, 1 b-hadron, 1 J/psi");
422 AliGenPythia* pythia = new AliGenPythia(-1);
423 pythia->SetMomentumRange(0, 999999.);
424 pythia->SetYRange(-4.5, -2.0);
425 pythia->SetPtRange(0,1000.);
426 pythia->SetProcess(kPyBeautyppMNRwmi);
427 pythia->SetEnergyCMS(energy);
430 pythia->SetTune(320);
431 pythia->UseNewMultipleInteractionsScenario();
434 pythia->SetCutOnChild(1);
435 pythia->SetPdgCodeParticleforAcceptanceCut(443);
436 pythia->SetChildYRange(-4.5, -2.0);
437 pythia->SetChildPtRange(0,10000.);
440 pythia->SetForceDecay(kBJpsiDiMuon);
446 //====================================================================================================================================================
448 AliGenerator* HijingParam() {
450 AliGenHIJINGpara *gener = new AliGenHIJINGpara(2000);
451 gener->SetMomentumRange(0,999);
452 gener->SetPhiRange(0,360);
453 gener->SetThetaRange(171, 179);
454 gener->SetOrigin(0,0,0); // vertex position
455 gener->SetSigma(0,0,0); // Sigma in (X,Y,Z) (cm) on IP position
462 //====================================================================================================================================================
464 AliGenerator* Hijing() {
466 AliGenHijing *gener = new AliGenHijing(-1);
467 // centre of mass energy
468 gener->SetEnergyCMS(energy);
469 gener->SetImpactParameterRange(bMin, bMax);
471 gener->SetReferenceFrame("CMS");
473 gener->SetProjectile("A", 208, 82);
474 gener->SetTarget ("A", 208, 82);
475 // tell hijing to keep the full parent child chain
476 gener->KeepFullEvent();
477 // enable jet quenching
478 gener->SetJetQuenching(1);
480 gener->SetShadowing(1);
481 // Don't track spectators
482 gener->SetSpectators(0);
483 // kinematic selection
484 gener->SetSelectAll(0);
488 //====================================================================================================================================================
490 AliGenerator* Hijing2500() {
492 AliGenHijing *gener = (AliGenHijing*) Hijing();
493 gener->SetJetQuenching(0);
494 gener->SetPtHardMin (4.5);
499 //====================================================================================================================================================
501 AliGenerator* Hijing2500Cocktail() {
503 comment = comment.Append(" PbPb: Hjing2500 at 5.5 + muon signals");
505 AliGenCocktail *cocktail = new AliGenCocktail();
506 cocktail->SetProjectile("A", 208, 82);
507 cocktail->SetTarget ("A", 208, 82);
508 cocktail->SetEnergyCMS(energy);
510 // 1 Hijing event: provides underlying event and collision geometry
511 AliGenHijing *hijing = Hijing2500();
512 cocktail->AddGenerator(hijing,"hijing",1);
514 // generator for the custom signal
515 AliGenerator* signalGen = 0x0;
518 signalGen = GenBox();
521 signalGen = GenMuonLMR();
524 signalGen = GenCorrHF();
527 signalGen = GenBox();
530 cocktail->AddGenerator(signalGen, "signal", 1);
532 cocktail->SetTrackingFlag(1);
538 //====================================================================================================================================================
540 AliGenerator* QEDGeneratorPbPb() {
542 AliGenEpEmv1 *generQED = new AliGenEpEmv1();
543 generQED->SetEnergyCMS(energy);
544 generQED->SetProjectile("A", 208, 82);
545 generQED->SetTarget ("A", 208, 82);
546 generQED->SetYRange(-8.5, 0.);
547 generQED->SetPtRange(0.4e-3, 1.0); // Set pt limits (GeV) for e+-
548 generQED->SetDebug(0); // Set debug level (0 = silent)
552 //====================================================================================================================================================
556 #if defined(__CINT__)
558 gSystem->Load("liblhapdf"); // Parton density functions
559 gSystem->Load("libEGPythia6"); // TGenerator interface
560 if (proc == kPythia6) {
561 gSystem->Load("libpythia6"); // Pythia 6.2
562 gSystem->Load("libAliPythia6"); // ALICE specific implementations
565 gSystem->Load("libpythia6.4.21"); // Pythia 6.4
566 gSystem->Load("libAliPythia6"); // ALICE specific implementations
569 if (proc == kHijing || proc == kHijing2500 || proc == kHijing2500Cocktail) {
570 gSystem->Load("libhijing");
571 gSystem->Load("libTHijing");
574 gSystem->Load("libgeant321");
580 //====================================================================================================================================================