git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
move production scripts out of test
[u/mrichter/AliRoot.git] / prod / LHC09a1 / Config.C
CommitLineData
afea418a 1//\r
2// Configuration for the Physics Data Challenge 2007\r
3//\r
4//\r
5// Assuming inel = 70 mb (PPR v.1, p.64)\r
6//\r
7// 84.98% of MSEL=0 events (including diffractive) with \r
8// QQbar switched off (these events will include injected J/psi) => kPyMbNoHvq\r
9// \r
10// 14.14% of MSEL=1 events with ccbar (in 4 subsamples) => kCharmpp14000wmi \r
11// bin 1 25% (3.535%): 2.76 < pthard < 3 GeV/c \r
12// bin 2 40% (5.656%): 3 < pthard < 4 GeV/c \r
13// bin 3 29% (4.101%): 4 < pthard < 8 GeV/c \r
14// bin 4 6% (0.848%): pthard > 8 GeV/c \r
15//\r
16// 0.73% of MSEL=1 events with bbbar (in 4 subsamples) => kBeautypp14000wmi\r
17// bin 1 5% (0.037%): 2.76 < pthard < 4 GeV/c\r
18// bin 2 31% (0.226%): 4 < pthard < 6 GeV/c\r
19// bin 3 28% (0.204%): 6 < pthard < 8 GeV/c\r
20// bin 4 36% (0.263%): pthard >8 GeV/c\r
21//\r
22// 0.075% of MSEL=0 events with QQbar switched off and 1 Omega- => kPyOmegaMinus\r
23// 0.075% of MSEL=0 events with QQbar switched off and 1 OmegaBar+ => kPyOmegaPlus\r
24\r
25// One can use the configuration macro in compiled mode by\r
26// root [0] gSystem->Load("libgeant321");\r
27// root [0] gSystem->SetIncludePath("-I$ROOTSYS/include -I$ALICE_ROOT/include\\r
28// -I$ALICE_ROOT -I$ALICE/geant3/TGeant3");\r
29// root [0] .x grun.C(1,"Config_PDC07.C++")\r
30\r
31class AliGenPythia;\r
32\r
33#if !defined(__CINT__) || defined(__MAKECINT__)\r
34#include <Riostream.h>\r
35#include <TRandom.h>\r
36#include <TDatime.h>\r
37#include <TSystem.h>\r
38#include <TVirtualMC.h>\r
39#include <TGeant3TGeo.h>\r
40#include "EVGEN/AliGenCocktail.h"\r
41#include "EVGEN/AliGenParam.h"\r
42#include "EVGEN/AliGenMUONlib.h"\r
43#include "STEER/AliRunLoader.h"\r
44#include "STEER/AliRun.h"\r
45#include "STEER/AliConfig.h"\r
46#include "PYTHIA6/AliDecayerPythia.h"\r
47#include "PYTHIA6/AliGenPythia.h"\r
99c7d495 48#include "STEER/AliMagF.h"\r
afea418a 49#include "STRUCT/AliBODY.h"\r
50#include "STRUCT/AliMAG.h"\r
51#include "STRUCT/AliABSOv3.h"\r
52#include "STRUCT/AliDIPOv3.h"\r
53#include "STRUCT/AliHALLv3.h"\r
54#include "STRUCT/AliFRAMEv2.h"\r
55#include "STRUCT/AliSHILv3.h"\r
56#include "STRUCT/AliPIPEv3.h"\r
57#include "ITS/AliITSgeom.h"\r
58#include "ITS/AliITSv11Hybrid.h"\r
59#include "TPC/AliTPCv2.h"\r
60#include "TOF/AliTOFv6T0.h"\r
61#include "HMPID/AliHMPIDv3.h"\r
62#include "HMPID/AliHMPIDv2.h"\r
63#include "ZDC/AliZDCv3.h"\r
64#include "TRD/AliTRDv1.h"\r
65#include "FMD/AliFMDv1.h"\r
66#include "MUON/AliMUONv1.h"\r
67#include "PHOS/AliPHOSv1.h"\r
68#include "PMD/AliPMDv1.h"\r
69#include "T0/AliT0v1.h"\r
70#include "EMCAL/AliEMCALv2.h"\r
71#include "ACORDE/AliACORDEv1.h"\r
72#include "VZERO/AliVZEROv7.h"\r
73#endif\r
74\r
75\r
76 // libraries required by geant321\r
77\r
78\r
79\r
80enum PDC07Proc_t \r
81{\r
82//--- Heavy Flavour Production ---\r
83 kCharmPbPb5500, kCharmpPb8800, kCharmpp14000, kCharmpp14000wmi,\r
84 kD0PbPb5500, kD0pPb8800, kD0pp14000,\r
85 kDPlusPbPb5500, kDPluspPb8800, kDPluspp14000,\r
86 kBeautyPbPb5500, kBeautypPb8800, kBeautypp14000, kBeautypp14000wmi, \r
87// -- Pythia Mb\r
88 kPyMbNoHvq, kPyOmegaPlus, kPyOmegaMinus, kPyJetJet, kPyJetJetFixed,\r
89 kPyGammaJetPHOS, kPyJetJetPHOS, kPyJetJetPHOSv2, kPyGammaBremsPHOS, \r
90 kPyGammaJetEMCAL, kPyJetJetEMCAL, kPyGammaBremsEMCAL, kRunMax\r
91};\r
92\r
93const char * pprRunName[] = {\r
94 "kCharmPbPb5500", "kCharmpPb8800", "kCharmpp14000", "kCharmpp14000wmi",\r
95 "kD0PbPb5500", "kD0pPb8800", "kD0pp14000",\r
96 "kDPlusPbPb5500", "kDPluspPb8800", "kDPluspp14000",\r
97 "kBeautyPbPb5500", "kBeautypPb8800", "kBeautypp14000", "kBeautypp14000wmi", \r
98 "kPyMbNoHvq", "kPyOmegaPlus", "kPyOmegaMinus", "kPyJetJet","kPyJetJetFixed",\r
99 "kPyGammaJetPHOS", "kPyJetJetPHOS", "kPyJetJetPHOSv2", "kPyGammaBremsPHOS",\r
100 "kPyGammaJetEMCAL", "kPyJetJetEMCAL", "kPyGammaBremsEMCAL"\r
101};\r
102\r
103\r
104//--- Decay Mode ---\r
105enum DecayHvFl_t \r
106{\r
107 kNature, kHadr, kSemiEl, kSemiMu\r
108};\r
109//--- Magnetic Field ---\r
110enum Mag_t\r
111{\r
112 kNoField, k5kG\r
113};\r
114\r
115//--- Trigger config ---\r
116enum TrigConf_t\r
117{\r
118 kDefaultPPTrig, kDefaultPbPbTrig\r
119};\r
120\r
121const char * TrigConfName[] = {\r
122 "p-p","Pb-Pb"\r
123};\r
124\r
125Float_t eCMS=14000;\r
126PDC07Proc_t proc = kPyJetJetEMCAL;\r
127Int_t year = 2009;\r
128\r
129\r
130enum PprGeo_t\r
131 {\r
132 kHoles, kNoHoles\r
133 };\r
134\r
135static PprGeo_t geo = kHoles;\r
136\r
137//--- Functions ---\r
138AliGenPythia *PythiaHVQ(PDC07Proc_t proc);\r
139\r
140AliGenPythia* Blubb(Int_t proc) {\r
141 //*******************************************************************//\r
142 // Configuration file for hard QCD processes generation with PYTHIA //\r
143 // //\r
144 //*******************************************************************//\r
145 AliGenPythia * gener = 0x0;\r
146\r
147 switch(proc) {\r
148 \r
149 case kPyJetJet:\r
150 comment = comment.Append(Form("pp->jet + jet at %3.0f TeV, pt hard %3.0f - %3.0f",eCMS,ptHardMin,ptHardMax));\r
151 AliGenPythia * gener = new AliGenPythia(nEvts);\r
152 gener->SetEnergyCMS(eCMS);// Centre of mass energy\r
153 gener->SetProcess(kPyJets);// Process type\r
154 gener->SetJetEtaRange(-1.5, 1.5);// Final state kinematic cuts\r
155 gener->SetJetPhiRange(0., 360.);\r
156 gener->SetJetEtRange(10., 1000.);\r
157 gener->SetPtHard(ptHardMin, ptHardMax);// Pt transfer of the hard scattering\r
158 gener->SetStrucFunc(kCTEQ5L);\r
159 \r
160 return gener;\r
161 \r
162 case kPyGammaJetPHOS:\r
163 comment = comment.Append(" pp->jet + gamma over PHOS");\r
164 gener = new AliGenPythia(nEvts);\r
165 gener->SetEnergyCMS(eCMS);\r
166 gener->SetProcess(kPyDirectGamma);\r
167 gener->SetStrucFunc(kCTEQ4L);\r
168 gener->SetPtHard(ptHardMin,ptHardMax);\r
169 //gener->SetYHard(-1.0,1.0);\r
170 gener->SetGammaEtaRange(-0.13,0.13);\r
171 gener->SetGammaPhiRange(218.,322.);//Over 5 modules +-2 deg\r
172 break;\r
173 case kPyJetJetPHOS:\r
174 comment = comment.Append(" pp->jet + jet over PHOS");\r
175 gener = new AliGenPythia(nEvts);\r
176 gener->SetEnergyCMS(eCMS);\r
177 gener->SetProcess(kPyJets);\r
178 gener->SetStrucFunc(kCTEQ4L);\r
179 gener->SetPtHard(ptHardMin,ptHardMax);\r
180 //gener->SetYHard(-1.0,1.0);\r
181 gener->SetJetEtaRange(-1.,1.);\r
182 gener->SetJetPhiRange(200.,340.); \r
183 gener->SetPi0InPHOS(kTRUE);\r
184 gener->SetFragPhotonOrPi0MinPt(ptGammaPi0Min);\r
185 \r
186 printf("\n \n Event generator: Minimum pT of particle in calorimeter %f \n \n", ptGammaPi0Min);\r
187 break;\r
188 case kPyGammaBremsPHOS:\r
189 comment = comment.Append(" pp->jet + jet+bremsphoton over PHOS at 14 TeV");\r
190 gener = new AliGenPythia(nEvts);\r
191 gener->SetEnergyCMS(eCMS);\r
192 gener->SetProcess(kPyJets);\r
193 gener->SetStrucFunc(kCTEQ4L);\r
194 gener->SetPtHard(ptHardMin,ptHardMax);\r
195 //gener->SetYHard(-1.0,1.0);\r
196 gener->SetJetEtaRange(-1.,1.);\r
197 gener->SetJetPhiRange(200.,340.); \r
198 gener->SetFragPhotonInPHOS(kTRUE);\r
199 gener->SetFragPhotonOrPi0MinPt(ptGammaPi0Min);\r
200 printf("\n \n Event generator: Minimum pT of particle in calorimeter %f \n \n", ptGammaPi0Min);\r
201 break;\r
202 case kPyJetJetPHOSv2:\r
203 comment = comment.Append(" pp->jet + jet over PHOS version2 ");\r
204 gener = new AliGenPythia(nEvts);\r
205 gener->SetEnergyCMS(eCMS);\r
206 gener->SetProcess(kPyJets);\r
207 gener->SetStrucFunc(kCTEQ4L);\r
208 gener->SetPtHard(ptHardMin,ptHardMax);\r
209 //gener->SetYHard(-1.0,1.0);\r
210 gener->SetJetEtaRange(-1.,1.);\r
211 gener->SetJetPhiRange(200.,340.); \r
212 //gener->SetPi0InPHOS(kTRUE);\r
213 gener->SetPhotonInPHOSeta(kTRUE);\r
214 gener->SetPhotonMinPt(ptGammaPi0Min);\r
215 gener->SetForceDecay(kAll);\r
216 break;\r
217 case kPyGammaJetEMCAL:\r
218 comment = comment.Append(" pp->jet + gamma over EMCAL at 14 TeV");\r
219 gener = new AliGenPythia(nEvts);\r
220 gener->SetEnergyCMS(eCMS);\r
221 gener->SetProcess(kPyDirectGamma);\r
222 gener->SetStrucFunc(kCTEQ4L);\r
223 gener->SetPtHard(ptHardMin,ptHardMax);\r
224 //gener->SetYHard(-1.0,1.0);\r
225 gener->SetGammaEtaRange(-0.71,0.71);\r
226 gener->SetGammaPhiRange(78.,192.);//Over 6 supermodules +-2 deg \r
227 break;\r
228 case kPyJetJetEMCAL:\r
229 comment = comment.Append(" pp->jet + jet over EMCAL at 14 TeV");\r
230 gener = new AliGenPythia(nEvts);\r
231 gener->SetEnergyCMS(eCMS);\r
232 gener->SetProcess(kPyJets);\r
233 gener->SetStrucFunc(kCTEQ4L);\r
234 gener->SetPtHard(ptHardMin,ptHardMax);\r
235 //gener->SetYHard(-1.0,1.0);\r
236 gener->SetJetEtaRange(-1,1);\r
237 gener->SetJetPhiRange(60.,210.); \r
238 gener->SetPi0InEMCAL(kTRUE);\r
239 gener->SetFragPhotonOrPi0MinPt(ptGammaPi0Min);\r
240 printf("\n \n Event generator: Minimum pT of particle in calorimeter %f \n \n", ptGammaPi0Min);\r
241 break;\r
242 case kPyGammaBremsEMCAL:\r
243 comment = comment.Append(" pp->jet + jet+bremsphoton over EMCAL at 14 TeV");\r
244 gener = new AliGenPythia(nEvts);\r
245 gener->SetEnergyCMS(eCMS);\r
246 gener->SetProcess(kPyJets);\r
247 gener->SetStrucFunc(kCTEQ4L);\r
248 gener->SetPtHard(ptHardMin,ptHardMax);\r
249 //gener->SetYHard(-1.0,1.0);\r
250 gener->SetJetEtaRange(-1,1);\r
251 gener->SetJetPhiRange(60.,210.); //Over 2 uncovered PHOS modules\r
252 gener->SetFragPhotonInEMCAL(kTRUE);\r
253 gener->SetFragPhotonOrPi0MinPt(ptGammaPi0Min);\r
254 \r
255 printf("\n \n Event generator: Minimum pT of particle in calorimeter %f \n \n", ptGammaPi0Min);\r
256 break;\r
257 case kPyJetJetFixed:\r
258 comment = comment.Append(" pp->jet + jet+jet in Fixed bin for testing");\r
259 gener = new AliGenPythia(nEvts);\r
260 gener->SetEnergyCMS(14000);// Centre of mass energy\r
261 gener->SetProcess(kPyJets);// Process type\r
262 gener->SetJetEtaRange(-1.1, 1.1);// Final state kinematic cuts\r
263 gener->SetJetPhiRange(0., 360.);\r
264 gener->SetJetEtRange(95., 105.);\r
265 gener->SetPtHard(80, 1000);// Pt transfer of the hard scattering\r
266 gener->SetStrucFunc(kCTEQ5L);\r
267 gener->SetGluonRadiation(1,1); \r
268 \r
269\r
270 printf("\n \n Fixed Parameter \n");\r
271 break;\r
272\r
273 }\r
274\r
275 return gener;\r
276}\r
277\r
278AliGenerator *MbCocktail();\r
279AliGenerator *PyMbTriggered(Int_t pdg);\r
280void ProcessEnvironmentVars();\r
281\r
282// This part for configuration\r
283static DecayHvFl_t decHvFl = kNature; \r
284static Mag_t mag = k5kG; \r
285static TrigConf_t trig = kDefaultPPTrig; // default pp trigger configuration\r
286static Int_t runNumber= 0; \r
287static Int_t eventNumber= 0; \r
288//========================//\r
289// Set Random Number seed //\r
290//========================//\r
291TDatime dt;\r
292static UInt_t seed = dt.Get();\r
293\r
294// nEvts = -1 : you get 1 QQbar pair and all the fragmentation and \r
295// decay chain\r
296// nEvts = N>0 : you get N charm / beauty Hadrons \r
297Int_t nEvts = -1; \r
298// stars = kTRUE : all heavy resonances and their decay stored\r
299// = kFALSE: only final heavy hadrons and their decays stored\r
300Bool_t stars = kTRUE;\r
301\r
302// To be used only with kCharmpp14000wmi and kBeautypp14000wmi\r
303// To get a "reasonable" agreement with MNR results, events have to be \r
304// generated with the minimum ptHard set to 2.76 GeV.\r
305// To get a "perfect" agreement with MNR results, events have to be \r
306// generated in four ptHard bins with the following relative \r
307// normalizations:\r
308// CHARM\r
309// 2.76-3 GeV: 25%\r
310// 3-4 GeV: 40%\r
311// 4-8 GeV: 29%\r
312// >8 GeV: 6%\r
313// BEAUTY\r
314// 2.76-4 GeV: 5% \r
315// 4-6 GeV: 31%\r
316// 6-8 GeV: 28%\r
317// >8 GeV: 36%\r
318\r
319Float_t ptHardMin = 10.;\r
320Float_t ptHardMax = 20.;\r
321Float_t ptGammaPi0Min = 1.;\r
322Int_t iquenching = 0;\r
323Float_t qhat = 20.;\r
324\r
325// Comment line\r
326static TString comment;\r
327\r
328void Config()\r
329{\r
330 \r
331\r
332 // Get settings from environment variables\r
333 ProcessEnvironmentVars();\r
334\r
335#if defined(__CINT__)\r
336 gSystem->Load("liblhapdf");\r
337 gSystem->Load("libEGPythia6");\r
338 gSystem->Load("libpythia6");\r
339 gSystem->Load("libAliPythia6");\r
340 gSystem->Load("libgeant321");\r
341#endif\r
342\r
343 new TGeant3TGeo("C++ Interface to Geant3");\r
344\r
345 // Output every 100 tracks\r
346 ((TGeant3*)gMC)->SetSWIT(4,100);\r
347\r
348 //=======================================================================\r
349 // Run loader \r
350 AliRunLoader* rl=0x0;\r
351 rl = AliRunLoader::Open("galice.root",\r
352 AliConfig::GetDefaultEventFolderName(),\r
353 "recreate");\r
354 if (rl == 0x0)\r
355 {\r
356 gAlice->Fatal("Config.C","Can not instatiate the Run Loader");\r
357 return;\r
358 }\r
359 rl->SetCompressionLevel(2);\r
360 rl->SetNumberOfEventsPerFile(1000);\r
361 gAlice->SetRunLoader(rl);\r
362\r
363 // Run number\r
364 //gAlice->SetRunNumber(runNumber);\r
365 \r
366 // Set the trigger configuration\r
367 gAlice->SetTriggerDescriptor(TrigConfName[trig]);\r
368 cout<<"Trigger configuration is set to "<<TrigConfName[trig]<<endl;\r
369\r
370 //\r
371 //=======================================================================\r
372 // ************* STEERING parameters FOR ALICE SIMULATION **************\r
373 // --- Specify event type to be tracked through the ALICE setup\r
374 // --- All positions are in cm, angles in degrees, and P and E in GeV\r
375\r
376\r
377 gMC->SetProcess("DCAY",1);\r
378 gMC->SetProcess("PAIR",1);\r
379 gMC->SetProcess("COMP",1);\r
380 gMC->SetProcess("PHOT",1);\r
381 gMC->SetProcess("PFIS",0);\r
382 gMC->SetProcess("DRAY",0);\r
383 gMC->SetProcess("ANNI",1);\r
384 gMC->SetProcess("BREM",1);\r
385 gMC->SetProcess("MUNU",1);\r
386 gMC->SetProcess("CKOV",1);\r
387 gMC->SetProcess("HADR",1);\r
388 gMC->SetProcess("LOSS",2);\r
389 gMC->SetProcess("MULS",1);\r
390 gMC->SetProcess("RAYL",1);\r
391\r
392 Float_t cut = 1.e-3; // 1MeV cut by default\r
393 Float_t tofmax = 1.e10;\r
394\r
395 gMC->SetCut("CUTGAM", cut);\r
396 gMC->SetCut("CUTELE", cut);\r
397 gMC->SetCut("CUTNEU", cut);\r
398 gMC->SetCut("CUTHAD", cut);\r
399 gMC->SetCut("CUTMUO", cut);\r
400 gMC->SetCut("BCUTE", cut); \r
401 gMC->SetCut("BCUTM", cut); \r
402 gMC->SetCut("DCUTE", cut); \r
403 gMC->SetCut("DCUTM", cut); \r
404 gMC->SetCut("PPCUTM", cut);\r
405 gMC->SetCut("TOFMAX", tofmax); \r
406\r
407 //======================//\r
408 // Set External decayer //\r
409 //======================//\r
410 TVirtualMCDecayer* decayer = new AliDecayerPythia();\r
411 // DECAYS\r
412 //\r
413 switch(decHvFl) {\r
414 case kNature:\r
415 decayer->SetForceDecay(kAll);\r
416 break;\r
417 case kHadr:\r
418 decayer->SetForceDecay(kHadronicD);\r
419 break;\r
420 case kSemiEl:\r
421 decayer->SetForceDecay(kSemiElectronic);\r
422 break;\r
423 case kSemiMu:\r
424 decayer->SetForceDecay(kSemiMuonic);\r
425 break;\r
426 }\r
427 decayer->Init();\r
428 gMC->SetExternalDecayer(decayer);\r
429 if(proc == kPyJetJetPHOSv2) // in this case we decay the pi0\r
430 decayer->SetForceDecay(kNeutralPion);\r
431\r
432 //=========================//\r
433 // Generator Configuration //\r
434 //=========================//\r
435 AliGenPythia* gener = 0x0;\r
436 \r
437 if (proc <= kBeautypp14000wmi) {\r
438 AliGenPythia *pythia = PythiaHVQ(proc);\r
439 // FeedDown option\r
440 pythia->SetFeedDownHigherFamily(kFALSE);\r
441 // Stack filling option\r
442 if(!stars) pythia->SetStackFillOpt(AliGenPythia::kParentSelection);\r
443 // Set Count mode\r
444 if(nEvts>0) pythia->SetCountMode(AliGenPythia::kCountParents);\r
445 //\r
446 // DECAYS\r
447 // \r
448 switch(decHvFl) {\r
449 case kNature:\r
450 pythia->SetForceDecay(kAll);\r
451 break;\r
452 case kHadr:\r
453 pythia->SetForceDecay(kHadronicD);\r
454 break;\r
455 case kSemiEl:\r
456 pythia->SetForceDecay(kSemiElectronic);\r
457 break;\r
458 case kSemiMu:\r
459 pythia->SetForceDecay(kSemiMuonic);\r
460 break;\r
461 }\r
462 //\r
463 // GEOM & KINE CUTS\r
464 //\r
465 pythia->SetMomentumRange(0,99999999);\r
466 pythia->SetPhiRange(0., 360.);\r
467 pythia->SetThetaRange(0,180);\r
468 switch(ycut) {\r
469 case kFull:\r
470 pythia->SetYRange(-999,999);\r
471 break;\r
472 case kBarrel:\r
473 pythia->SetYRange(-2,2);\r
474 break;\r
475 case kMuonArm:\r
476 pythia->SetYRange(1,6);\r
477 break;\r
478 }\r
479 gener = pythia;\r
480 } else if (proc == kPyMbNoHvq) {\r
481 gener = MbCocktail();\r
482 } else if (proc == kPyOmegaMinus) {\r
483 gener = PyMbTriggered(3334);\r
484 } else if (proc == kPyOmegaPlus) {\r
485 gener = PyMbTriggered(-3334);\r
486 } else if (proc <= kPyGammaBremsEMCAL) {\r
487 \r
488 AliGenPythia *pythia = Blubb((Int_t) proc ); \r
489 \r
490 // FeedDown option\r
491 pythia->SetFeedDownHigherFamily(kFALSE);\r
492 // Set Count mode\r
493 if(nEvts>0) pythia->SetCountMode(AliGenPythia::kCountParents);\r
494 \r
495 //\r
496 // GEOM & KINE CUTS\r
497 //\r
498 pythia->SetMomentumRange(0,99999999);\r
499 // pythia->SetJetEtaRange(-1.5, 1.5);// Final state kinematic cuts\r
500 // pythia->SetJetPhiRange(0., 360.);\r
501 // pythia->SetThetaRange(45,135);\r
502 \r
503 if(proc == kPyJetJetPHOSv2)\r
504 pythia->SetForceDecay(kNeutralPion);\r
505 else \r
506 pythia->SetForceDecay(kAll);\r
507 \r
508 pythia->SetPycellParameters(2., 274, 432, 0., 4., 5., 1.0);\r
509 pythia->SetPtKick(3.); // set the intrinsic kt to 5 GeV/c\r
510 pythia->SetGluonRadiation(1,1);\r
511 gener = pythia;\r
512 }\r
513 \r
514\r
515 // PRIMARY VERTEX\r
516\r
517 gener->SetOrigin(0., 0., 0.); // vertex position\r
518\r
519 // Size of the interaction diamond\r
520 // Longitudinal\r
521 Float_t sigmaz = 7.55 / TMath::Sqrt(2.); // [cm]\r
522\r
523 // Transverse\r
524 Float_t betast = 10; // beta* [m]\r
525 Float_t eps = 3.75e-6; // emittance [m]\r
526 Float_t gamma = 7000. / 0.938272; // relativistic gamma [1]\r
527 Float_t sigmaxy = TMath::Sqrt(eps * betast / gamma) / TMath::Sqrt(2.) * 100.; // [cm]\r
528 printf("\n \n Diamond size x-y: %10.3e z: %10.3e\n \n", sigmaxy, sigmaz);\r
529 \r
530 gener->SetSigma(sigmaxy, sigmaxy, sigmaz); // Sigma in (X,Y,Z) (cm) on IP position\r
531 gener->SetCutVertexZ(3.); // Truncate at 3 sigma\r
532 gener->SetVertexSmear(kPerEvent);\r
533\r
534 gener->Init();\r
535\r
536 //Quenching\r
537 gener->SetQuench(iquenching);\r
538 if(iquenching == 1){\r
539 Float_t k = 6e5*(qhat/1.7) ; //qhat=1.7, k = 6e5, default value \r
540 AliPythia::Instance()->InitQuenching(0.,0.1,k,0,0.95,6);\r
541\r
542 }\r
543 // FIELD\r
544 \r
99c7d495 545 AliMagF* field = 0x0;\r
afea418a 546 if (mag == kNoField) {\r
547 comment = comment.Append(" | L3 field 0.0 T");\r
4642ac4b 548 field = new AliMagF("Maps","Maps", 0., 0., AliMagF::k5kGUniform);\r
afea418a 549 } else if (mag == k5kG) {\r
550 comment = comment.Append(" | L3 field 0.5 T");\r
551 // field = new AliMagWrapCheb("Maps","Maps", 2, 1., 10., AliMagWrapCheb::k5kG);\r
4642ac4b 552 field = new AliMagF("Maps","Maps", -1., -1., AliMagF::k5kG);\r
afea418a 553 }\r
554\r
555 printf("\n \n Comment: %s \n \n", comment.Data());\r
99c7d495 556 TGeoGlobalMagField::Instance()->SetField(field);\r
afea418a 557 rl->CdGAFile();\r
afea418a 558\r
559\r
560\r
561 Int_t iABSO = 1;\r
562 Int_t iACORDE = 0;\r
563 Int_t iDIPO = 1;\r
564 Int_t iEMCAL = 1;\r
565 Int_t iFMD = 1;\r
566 Int_t iFRAME = 1;\r
567 Int_t iHALL = 1;\r
568 Int_t iITS = 1;\r
569 Int_t iMAG = 1;\r
570 Int_t iMUON = 1;\r
571 Int_t iPHOS = 1;\r
572 Int_t iPIPE = 1;\r
573 Int_t iPMD = 0 ;\r
574 Int_t iHMPID = 1;\r
575 Int_t iSHIL = 1;\r
576 Int_t iT0 = 1;\r
577 Int_t iTOF = 1;\r
578 Int_t iTPC = 1;\r
579 Int_t iTRD = 1;\r
580 Int_t iVZERO = 1;\r
581 Int_t iZDC = 1;\r
582 \r
583\r
584 //=================== Alice BODY parameters =============================\r
585 AliBODY *BODY = new AliBODY("BODY", "Alice envelop");\r
586\r
587\r
588 if (iMAG)\r
589 {\r
590 //=================== MAG parameters ============================\r
591 // --- Start with Magnet since detector layouts may be depending ---\r
592 // --- on the selected Magnet dimensions ---\r
593 AliMAG *MAG = new AliMAG("MAG", "Magnet");\r
594 }\r
595\r
596\r
597 if (iABSO)\r
598 {\r
599 //=================== ABSO parameters ============================\r
600 AliABSO *ABSO = new AliABSOv3("ABSO", "Muon Absorber");\r
601 }\r
602\r
603 if (iDIPO)\r
604 {\r
605 //=================== DIPO parameters ============================\r
606\r
607 AliDIPO *DIPO = new AliDIPOv3("DIPO", "Dipole version 3");\r
608 }\r
609\r
610 if (iHALL)\r
611 {\r
612 //=================== HALL parameters ============================\r
613\r
614 AliHALL *HALL = new AliHALLv3("HALL", "Alice Hall");\r
615 }\r
616\r
617\r
618 if (iFRAME)\r
619 {\r
620 //=================== FRAME parameters ============================\r
621\r
622 AliFRAMEv2 *FRAME = new AliFRAMEv2("FRAME", "Space Frame");\r
623 if (geo == kHoles) FRAME->SetHoles(1);\r
624 else FRAME->SetHoles(0);\r
625\r
626 }\r
627\r
628 if (iSHIL)\r
629 {\r
630 //=================== SHIL parameters ============================\r
631\r
632 AliSHIL *SHIL = new AliSHILv3("SHIL", "Shielding Version 3");\r
633 }\r
634\r
635\r
636 if (iPIPE)\r
637 {\r
638 //=================== PIPE parameters ============================\r
639\r
640 AliPIPE *PIPE = new AliPIPEv3("PIPE", "Beam Pipe");\r
641 }\r
642 \r
643 if(iITS) {\r
644\r
645 //=================== ITS parameters ============================\r
646 //\r
647 AliITS *ITS = new AliITSv11Hybrid("ITS","ITS v11Hybrid");\r
648 // AliITSvPPRasymmFMD *ITS = new AliITSvPPRasymmFMD("ITS","New ITS PPR detailed version with asymmetric services");\r
649 }\r
650\r
651 if (iTPC)\r
652 {\r
653 //============================ TPC parameters =====================\r
654 AliTPC *TPC = new AliTPCv2("TPC", "Default");\r
655 }\r
656\r
657\r
658 if (iTOF) {\r
659 //=================== TOF parameters ============================\r
660 AliTOF *TOF = new AliTOFv6T0("TOF", "normal TOF");\r
661 }\r
662\r
663\r
664 if (iHMPID)\r
665 {\r
666 //=================== HMPID parameters ===========================\r
667 // AliHMPID *HMPID = new AliHMPIDv2("HMPID", "normal HMPID");\r
668 AliHMPID *HMPID = new AliHMPIDv3("HMPID", "normal HMPID");\r
669\r
670 }\r
671\r
672\r
673 if (iZDC)\r
674 {\r
675 //=================== ZDC parameters ============================\r
676\r
677 AliZDC *ZDC = new AliZDCv3("ZDC", "normal ZDC");\r
678 }\r
679\r
680 if (iTRD)\r
681 {\r
682 //=================== TRD parameters ============================\r
683\r
684 AliTRD *TRD = new AliTRDv1("TRD", "TRD slow simulator");\r
685 AliTRDgeometry *geoTRD = TRD->GetGeometry();\r
686 // Partial geometry: modules at 0,8,9,17\r
687 // Partial geometry: modules at 1,7,10,16 expected for 2009\r
688 // starting at 3h in positive direction\r
689\r
690\r
691 geoTRD->SetSMstatus(2,0);\r
692 geoTRD->SetSMstatus(3,0);\r
693 geoTRD->SetSMstatus(4,0);\r
694 geoTRD->SetSMstatus(5,0);\r
695 geoTRD->SetSMstatus(6,0);\r
696 geoTRD->SetSMstatus(11,0);\r
697 geoTRD->SetSMstatus(12,0);\r
698 geoTRD->SetSMstatus(13,0);\r
699 geoTRD->SetSMstatus(14,0);\r
700 geoTRD->SetSMstatus(15,0);\r
701\r
702\r
703\r
704\r
705 /*\r
706 // Partial geometry: modules at 2,3,4,6,11,12,14,15\r
707 // starting at 6h in positive direction\r
708 //Hole on SM 13 and 14 for PHOS\r
709 geoTRD->SetSMstatus(13,0);\r
710 geoTRD->SetSMstatus(14,0);\r
711 geoTRD->SetSMstatus(15,0);\r
712 */\r
713\r
714 }\r
715\r
716 if (iFMD)\r
717 {\r
718 //=================== FMD parameters ============================\r
719 AliFMD *FMD = new AliFMDv1("FMD", "normal FMD");\r
720 }\r
721\r
722 if (iMUON)\r
723 {\r
724 //=================== MUON parameters ===========================\r
725 // New MUONv1 version (geometry defined via builders)\r
726 AliMUON *MUON = new AliMUONv1("MUON", "default");\r
727 }\r
728 //=================== PHOS parameters ===========================\r
729\r
730 if (iPHOS)\r
731 {\r
732 AliPHOS *PHOS = new AliPHOSv1("PHOS", "IHEP");\r
733 // Default for cold phos...\r
734 }\r
735\r
736\r
737 if (iPMD)\r
738 {\r
739 //=================== PMD parameters ============================\r
740 AliPMD *PMD = new AliPMDv1("PMD", "normal PMD");\r
741 }\r
742\r
743 if (iT0)\r
744 {\r
745 //=================== T0 parameters ============================\r
746 AliT0 *T0 = new AliT0v1("T0", "T0 Detector");\r
747 }\r
748\r
749 if (iEMCAL)\r
750 {\r
751 //=================== EMCAL parameters ============================\r
752 AliEMCAL *EMCAL = new AliEMCALv2("EMCAL", "EMCAL_COMPLETE");\r
753 }\r
754\r
755 if (iACORDE)\r
756 {\r
757 //=================== CRT parameters ============================\r
758 AliACORDE *ACORDE = new AliACORDEv1("CRT", "normal ACORDE");\r
759 }\r
760\r
761 if (iVZERO)\r
762 {\r
763 //=================== CRT parameters ============================\r
764 AliVZERO *VZERO = new AliVZEROv7("VZERO", "normal VZERO");\r
765 }\r
766}\r
767\r
768// PYTHIA\r
769\r
770AliGenPythia *PythiaHVQ(PDC07Proc_t proc) {\r
771//*******************************************************************//\r
772// Configuration file for charm / beauty generation with PYTHIA //\r
773// //\r
774// The parameters have been tuned in order to reproduce the inclusive//\r
775// heavy quark pt distribution given by the NLO pQCD calculation by //\r
776// Mangano, Nason and Ridolfi. //\r
777// //\r
778// For details and for the NORMALIZATION of the yields see: //\r
779// N.Carrer and A.Dainese, //\r
780// "Charm and beauty production at the LHC", //\r
781// ALICE-INT-2003-019, [arXiv:hep-ph/0311225]; //\r
782// PPR Chapter 6.6, CERN/LHCC 2005-030 (2005). //\r
783//*******************************************************************//\r
784 AliGenPythia * gener = 0x0;\r
785\r
786 switch(proc) {\r
787 case kCharmPbPb5500:\r
788 comment = comment.Append(" Charm in Pb-Pb at 5.5 TeV");\r
789 gener = new AliGenPythia(nEvts);\r
790 gener->SetProcess(kPyCharmPbPbMNR);\r
791 gener->SetStrucFunc(kCTEQ4L);\r
792 gener->SetPtHard(2.1,-1.0);\r
793 gener->SetEnergyCMS(5500.);\r
794 gener->SetNuclei(208,208);\r
795 break;\r
796 case kCharmpPb8800:\r
797 comment = comment.Append(" Charm in p-Pb at 8.8 TeV");\r
798 gener = new AliGenPythia(nEvts);\r
799 gener->SetProcess(kPyCharmpPbMNR);\r
800 gener->SetStrucFunc(kCTEQ4L);\r
801 gener->SetPtHard(2.1,-1.0);\r
802 gener->SetEnergyCMS(8800.);\r
803 gener->SetProjectile("P",1,1);\r
804 gener->SetTarget("Pb",208,82);\r
805 break;\r
806 case kCharmpp14000:\r
807 comment = comment.Append(" Charm in pp at 14 TeV");\r
808 gener = new AliGenPythia(nEvts);\r
809 gener->SetProcess(kPyCharmppMNR);\r
810 gener->SetStrucFunc(kCTEQ4L);\r
811 gener->SetPtHard(2.1,-1.0);\r
812 gener->SetEnergyCMS(14000.);\r
813 break;\r
814 case kCharmpp14000wmi:\r
815 comment = comment.Append(" Charm in pp at 14 TeV with mult. interactions");\r
816 gener = new AliGenPythia(-1);\r
817 gener->SetProcess(kPyCharmppMNRwmi);\r
818 gener->SetStrucFunc(kCTEQ5L);\r
819 gener->SetPtHard(ptHardMin,ptHardMax);\r
820 gener->SetEnergyCMS(14000.);\r
821 break;\r
822 case kD0PbPb5500:\r
823 comment = comment.Append(" D0 in Pb-Pb at 5.5 TeV");\r
824 gener = new AliGenPythia(nEvts);\r
825 gener->SetProcess(kPyD0PbPbMNR);\r
826 gener->SetStrucFunc(kCTEQ4L);\r
827 gener->SetPtHard(2.1,-1.0);\r
828 gener->SetEnergyCMS(5500.);\r
829 gener->SetNuclei(208,208);\r
830 break;\r
831 case kD0pPb8800:\r
832 comment = comment.Append(" D0 in p-Pb at 8.8 TeV");\r
833 gener = new AliGenPythia(nEvts);\r
834 gener->SetProcess(kPyD0pPbMNR);\r
835 gener->SetStrucFunc(kCTEQ4L);\r
836 gener->SetPtHard(2.1,-1.0);\r
837 gener->SetEnergyCMS(8800.);\r
838 gener->SetProjectile("P",1,1);\r
839 gener->SetTarget("Pb",208,82);\r
840 break;\r
841 case kD0pp14000:\r
842 comment = comment.Append(" D0 in pp at 14 TeV");\r
843 gener = new AliGenPythia(nEvts);\r
844 gener->SetProcess(kPyD0ppMNR);\r
845 gener->SetStrucFunc(kCTEQ4L);\r
846 gener->SetPtHard(2.1,-1.0);\r
847 gener->SetEnergyCMS(14000.);\r
848 break;\r
849 case kDPlusPbPb5500:\r
850 comment = comment.Append(" DPlus in Pb-Pb at 5.5 TeV");\r
851 gener = new AliGenPythia(nEvts);\r
852 gener->SetProcess(kPyDPlusPbPbMNR);\r
853 gener->SetStrucFunc(kCTEQ4L);\r
854 gener->SetPtHard(2.1,-1.0);\r
855 gener->SetEnergyCMS(5500.);\r
856 gener->SetNuclei(208,208);\r
857 break;\r
858 case kDPluspPb8800:\r
859 comment = comment.Append(" DPlus in p-Pb at 8.8 TeV");\r
860 gener = new AliGenPythia(nEvts);\r
861 gener->SetProcess(kPyDPluspPbMNR);\r
862 gener->SetStrucFunc(kCTEQ4L);\r
863 gener->SetPtHard(2.1,-1.0);\r
864 gener->SetEnergyCMS(8800.);\r
865 gener->SetProjectile("P",1,1);\r
866 gener->SetTarget("Pb",208,82);\r
867 break;\r
868 case kDPluspp14000:\r
869 comment = comment.Append(" DPlus in pp at 14 TeV");\r
870 gener = new AliGenPythia(nEvts);\r
871 gener->SetProcess(kPyDPlusppMNR);\r
872 gener->SetStrucFunc(kCTEQ4L);\r
873 gener->SetPtHard(2.1,-1.0);\r
874 gener->SetEnergyCMS(14000.);\r
875 break;\r
876 case kBeautyPbPb5500:\r
877 comment = comment.Append(" Beauty in Pb-Pb at 5.5 TeV");\r
878 gener = new AliGenPythia(nEvts);\r
879 gener->SetProcess(kPyBeautyPbPbMNR);\r
880 gener->SetStrucFunc(kCTEQ4L);\r
881 gener->SetPtHard(2.75,-1.0);\r
882 gener->SetEnergyCMS(5500.);\r
883 gener->SetNuclei(208,208);\r
884 break;\r
885 case kBeautypPb8800:\r
886 comment = comment.Append(" Beauty in p-Pb at 8.8 TeV");\r
887 gener = new AliGenPythia(nEvts);\r
888 gener->SetProcess(kPyBeautypPbMNR);\r
889 gener->SetStrucFunc(kCTEQ4L);\r
890 gener->SetPtHard(2.75,-1.0);\r
891 gener->SetEnergyCMS(8800.);\r
892 gener->SetProjectile("P",1,1);\r
893 gener->SetTarget("Pb",208,82);\r
894 break;\r
895 case kBeautypp14000:\r
896 comment = comment.Append(" Beauty in pp at 14 TeV");\r
897 gener = new AliGenPythia(nEvts);\r
898 gener->SetProcess(kPyBeautyppMNR);\r
899 gener->SetStrucFunc(kCTEQ4L);\r
900 gener->SetPtHard(2.75,-1.0);\r
901 gener->SetEnergyCMS(14000.);\r
902 break;\r
903 case kBeautypp14000wmi:\r
904 comment = comment.Append(" Beauty in pp at 14 TeV with mult. interactions");\r
905 gener = new AliGenPythia(-1);\r
906 gener->SetProcess(kPyBeautyppMNRwmi);\r
907 gener->SetStrucFunc(kCTEQ5L);\r
908 gener->SetPtHard(ptHardMin,ptHardMax);\r
909 gener->SetEnergyCMS(14000.);\r
910 break;\r
911 }\r
912\r
913 return gener;\r
914}\r
915\r
916AliGenerator* MbCocktail()\r
917{\r
918 comment = comment.Append(" pp at 14 TeV: Pythia low-pt, no heavy quarks + J/Psi from parameterisation");\r
919 AliGenCocktail * gener = new AliGenCocktail();\r
920 gener->UsePerEventRates();\r
921 \r
922// Pythia\r
923\r
924 AliGenPythia* pythia = new AliGenPythia(-1); \r
925 pythia->SetMomentumRange(0, 999999.);\r
926 pythia->SetThetaRange(0., 180.);\r
927 pythia->SetYRange(-12.,12.);\r
928 pythia->SetPtRange(0,1000.);\r
929 pythia->SetProcess(kPyMb);\r
930 pythia->SetEnergyCMS(14000.);\r
931 pythia->SwitchHFOff();\r
932 \r
933// J/Psi parameterisation\r
934\r
935 AliGenParam* jpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF scaled", "Jpsi");\r
936 jpsi->SetPtRange(0.,100.);\r
937 jpsi->SetYRange(-8., 8.);\r
938 jpsi->SetPhiRange(0., 360.);\r
939 jpsi->SetForceDecay(kAll);\r
940\r
941 gener->AddGenerator(pythia, "Pythia", 1.);\r
942 gener->AddGenerator(jpsi, "J/Psi", 8.e-4); \r
943 \r
944 return gener;\r
945}\r
946\r
947AliGenerator* PyMbTriggered(Int_t pdg)\r
948{\r
949 AliGenPythia* pythia = new AliGenPythia(-1); \r
950 pythia->SetMomentumRange(0, 999999.);\r
951 pythia->SetThetaRange(0., 180.);\r
952 pythia->SetYRange(-12.,12.);\r
953 pythia->SetPtRange(0,1000.);\r
954 pythia->SetProcess(kPyMb);\r
955 pythia->SetEnergyCMS(14000.);\r
956 pythia->SetTriggerParticle(pdg, 0.9);\r
957 return pythia;\r
958}\r
959\r
960void ProcessEnvironmentVars()\r
961{ \r
962 cout << "######################################" << endl;\r
963 cout << "## Processing environment variables ##" << endl;\r
964 cout << "######################################" << endl;\r
965 \r
966 // Run Number\r
967 if (gSystem->Getenv("DC_RUN")) {\r
968 runNumber = atoi(gSystem->Getenv("DC_RUN"));\r
969 }\r
970 //cout<<"Run number "<<runNumber<<endl;\r
971\r
972 // Event Number\r
973 if (gSystem->Getenv("DC_EVENT")) {\r
974 eventNumber = atoi(gSystem->Getenv("DC_EVENT"));\r
975 }\r
976 //cout<<"Event number "<<eventNumber<<endl;\r
977\r
978 // Random Number seed\r
979 if (gSystem->Getenv("CONFIG_SEED")) {\r
980 seed = atoi(gSystem->Getenv("CONFIG_SEED"));\r
981 }\r
982 else if(gSystem->Getenv("DC_EVENT") && gSystem->Getenv("DC_RUN")){\r
983 seed = runNumber * 100000 + eventNumber;\r
984 }\r
985\r
986 gRandom->SetSeed(seed);\r
987 \r
988 cout<<"////////////////////////////////////////////////////////////////////////////////////"<<endl;\r
989 cout<<"Seed for random number generation= "<< seed<<" "<< gRandom->GetSeed()<<endl; \r
990 cout<<"////////////////////////////////////////////////////////////////////////////////////"<<endl;\r
991\r
992 // Run type\r
993 if (gSystem->Getenv("DC_RUN_TYPE")) {\r
994 cout<<"run type "<<gSystem->Getenv("DC_RUN_TYPE")<<endl;\r
995 for (Int_t iRun = 0; iRun < kRunMax; iRun++) {\r
996 if (strcmp(gSystem->Getenv("DC_RUN_TYPE"), pprRunName[iRun])==0) {\r
997 proc = (PDC07Proc_t)iRun;\r
998 }\r
999 }\r
1000 }\r
1001 else \r
1002 cout << "Environment variable DC_RUN_TYPE is not defined" << endl;\r
1003\r
1004 if (gSystem->Getenv("ECMS"))\r
1005 eCMS = atof(gSystem->Getenv("ECMS"));\r
1006 if (gSystem->Getenv("PTHARDMIN"))\r
1007 ptHardMin = atof(gSystem->Getenv("PTHARDMIN"));\r
1008 if (gSystem->Getenv("PTHARDMAX"))\r
1009 ptHardMax = atof(gSystem->Getenv("PTHARDMAX"));\r
1010 if (gSystem->Getenv("PTGAMMAPI0MIN"))\r
1011 ptGammaPi0Min = atof(gSystem->Getenv("PTGAMMAPI0MIN"));\r
1012 if (gSystem->Getenv("QUENCHING"))\r
1013 iquenching = atof(gSystem->Getenv("QUENCHING"));\r
1014 if (gSystem->Getenv("QHAT"))\r
1015 qhat = atof(gSystem->Getenv("QHAT"));\r
1016\r
1017 cout<<">> Run type set to "<<pprRunName[proc]<<endl;\r
1018 cout<<">> Collision energy set to "<<eCMS <<endl;\r
1019 cout<<">> ptHard limits: "<<ptHardMin<<" to " <<ptHardMax<<" GeV"<<endl;\r
1020 cout<<">> pt gamma/pi0 threshold "<< ptGammaPi0Min<<" GeV "<<endl;\r
1021 cout<<">> quenching on? "<< iquenching<<" qhat "<<qhat<<endl;\r
1022\r
1023 cout << "######################################" << endl;\r
1024 cout << "######################################" << endl;\r
1025}\r