In DAs:
[u/mrichter/AliRoot.git] / macros / Config_gener.C
1 #if !defined(__CINT__) || defined(__MAKECINT__)
2 #include <Riostream.h>
3 #include <TRandom.h>
4 #include <TSystem.h>
5 #include <TVirtualMC.h>
6 #include <TGeant3.h>
7 #include "STEER/AliRunLoader.h"
8 #include "STEER/AliRun.h"
9 #include "STEER/AliConfig.h"
10 #include "PYTHIA6/AliDecayerPythia.h"
11 #include "EVGEN/AliGenCocktail.h"
12 #include "EVGEN/AliGenFixed.h"
13 #include "EVGEN/AliGenBox.h"
14 #include "EVGEN/AliGenScan.h"
15 #include "EVGEN/AliGenHIJINGpara.h"
16 #include "THijing/AliGenHijing.h"
17 #include "PYTHIA6/AliGenPythia.h"
18 #include "THerwig/AliGenHerwig.h"
19 #include "EVGEN/AliGenParam.h"
20 #include "EVGEN/AliGenMUONlib.h"
21 #include "EVGEN/AliGenPHOSlib.h"
22 #include "EVGEN/AliGenGSIlib.h"
23 #include "EVGEN/AliGenFLUKAsource.h"
24 #include "EVGEN/AliGenExtFile.h"
25 #include "EVGEN/AliGenHalo.h"
26 #include "EVGEN/AliGenReaderTreeK.h"
27 #include "STRUCT/AliBODY.h"
28 #include "STRUCT/AliMAG.h"
29 #endif
30
31 enum gentype_t {hijing, hijingParam, gun, box, pythia, herwig, mcatnlo, 
32                 param1, param2, param3, param4, 
33                 cocktail, fluka, halo, ntuple, scan};
34
35 gentype_t gentype = mcatnlo;
36
37 Int_t ntracks=1;
38
39 void Config()
40 {
41
42   // Set Random Number seed
43   gRandom->SetSeed(12345); //Set 0 to use the current time
44   cout<<"Seed for random number generation= "<<gRandom->GetSeed()<<endl; 
45
46
47   // libraries required by geant321
48 #if defined(__CINT__)
49   gSystem->Load("libgeant321");
50 #endif
51
52   new TGeant3TGeo("C++ Interface to Geant3");
53
54 //=======================================================================
55 //  Create the output file
56    
57   AliRunLoader* rl=0x0;
58
59   cout<<"Config.C: Creating Run Loader ..."<<endl;
60   rl = AliRunLoader::Open("galice.root",
61                           AliConfig::GetDefaultEventFolderName(),
62                           "recreate");
63   if (rl == 0x0)
64     {
65       gAlice->Fatal("Config.C","Can not instatiate the Run Loader");
66       return;
67     }
68   rl->SetCompressionLevel(2);
69   rl->SetNumberOfEventsPerFile(3);
70   gAlice->SetRunLoader(rl);
71
72 //
73 // Set External decayer
74   AliDecayer* decayer = new AliDecayerPythia();
75   decayer->SetForceDecay(kAll);
76   decayer->Init();
77   gMC->SetExternalDecayer(decayer);
78
79
80 //
81 //=======================================================================
82 // ************* STEERING parameters FOR ALICE SIMULATION **************
83 // --- Specify event type to be tracked through the ALICE setup
84 // --- All positions are in cm, angles in degrees, and P and E in GeV
85
86
87   gMC->SetProcess("DCAY",1);
88   gMC->SetProcess("PAIR",1);
89   gMC->SetProcess("COMP",1);
90   gMC->SetProcess("PHOT",1);
91   gMC->SetProcess("PFIS",0);
92   gMC->SetProcess("DRAY",0);
93   gMC->SetProcess("ANNI",1);
94   gMC->SetProcess("BREM",1);
95   gMC->SetProcess("MUNU",1);
96   gMC->SetProcess("CKOV",1);
97   gMC->SetProcess("HADR",1);
98   gMC->SetProcess("LOSS",2);
99   gMC->SetProcess("MULS",1);
100   gMC->SetProcess("RAYL",1);
101   
102   Float_t cut = 1.e-3;        // 1MeV cut by default
103   Float_t tofmax = 1.e10;
104   
105   gMC->SetCut("CUTGAM", cut);
106   gMC->SetCut("CUTELE", cut);
107   gMC->SetCut("CUTNEU", cut);
108   gMC->SetCut("CUTHAD", cut);
109   gMC->SetCut("CUTMUO", cut);
110   gMC->SetCut("BCUTE",  cut); 
111   gMC->SetCut("BCUTM",  cut); 
112   gMC->SetCut("DCUTE",  cut); 
113   gMC->SetCut("DCUTM",  cut); 
114   gMC->SetCut("PPCUTM", cut);
115   gMC->SetCut("TOFMAX", tofmax); 
116   
117
118   AliGenerator * gGener = 0x0;
119   switch(gentype)
120     {
121     case gun:
122 //*********************************************
123 // Example for Fixed Particle Gun             
124 //*********************************************
125       {
126         AliGenFixed *gener = new AliGenFixed(ntracks);
127         gener->SetMomentum(50);
128         gener->SetPhi(180.);
129         gener->SetTheta(5.);
130         gener->SetOrigin(0,0,0);        //vertex position
131         gener->SetPart(kNeutron);
132         gGener = gener;
133       }
134       break;
135     case box:  
136 //*********************************************
137 // Example for Moving Particle Gun            *
138 //*********************************************
139       {
140         AliGenBox *gener = new AliGenBox(ntracks);
141         gener->SetMomentumRange(3,4);
142         gener->SetPhiRange(0,360);
143         gener->SetThetaRange(90, 180. );
144         gener->SetOrigin(0,0,0);   
145         //vertex position
146         gener->SetSigma(0,0,0);         //Sigma in (X,Y,Z) (cm) on IP position
147         gener->SetPart(kMuonPlus);
148         gGener = gener;
149       }
150       break;
151     case scan:  
152 //*********************************************
153 // Scanning on a grid                         *
154 //*********************************************
155       {
156         AliGenScan *gener = new AliGenScan(-1);
157         gener->SetMomentumRange(4,4);
158         gener->SetPhiRange(0,360);
159         gener->SetThetaRange(9,9);
160         //vertex position
161         gener->SetSigma(6,6,0);         //Sigma in (X,Y,Z) (cm) on IP position
162         gener->SetPart(kMuonPlus); 
163         gener->SetRange(20, -100, 100, 20, -100, 100, 1, 500, 500);
164         gGener = gener;
165       }
166       break;
167      
168     case hijingParam:
169       {
170         AliGenHIJINGpara *gener = new AliGenHIJINGpara(ntracks);
171         gener->SetMomentumRange(0,999);
172         gener->SetPhiRange(0,360);
173         gener->SetThetaRange(2,10);
174         gener->SetOrigin(0,0,0);        //vertex position
175         gener->SetSigma(0,0,0);         //Sigma in (X,Y,Z) (cm) on IP position
176         gGener = gener;
177       }
178       break;
179     case hijing:
180       {
181         AliGenHijing *gener = new AliGenHijing(-1);
182 // centre of mass energy 
183         gener->SetEnergyCMS(5500);
184 // reference frame
185         gener->SetReferenceFrame("CMS     ");
186 // projectile
187         gener->SetProjectile("A       ", 208, 82);
188         gener->SetTarget    ("A       ", 208, 82);
189 // impact parameter range
190         gener->SetImpactParameterRange(0, 3.);
191 // evaluate cross section before run
192         gener->SetEvaluate(0);
193 // tell hijing to keep the full parent child chain
194         gener->KeepFullEvent();
195 // enable jet quenching
196         gener->SetJetQuenching(1);
197 // enable shadowing
198         gener->SetShadowing(1);
199 // neutral pion and heavy particle decays switched off
200         gener->SetDecaysOff(1);
201 // trigger
202         gener->SetTrigger(0);
203 // kinematic selection
204         gener->SetSelectAll(0);
205 // momentum range
206         gener->SetMomentumRange(0,999);
207 // phi range
208         gener->SetPhiRange(0.,360.);
209 // theta range 
210         gener->SetThetaRange(0,180.);
211 // select flavor (0: no, 4: charm+beauty, 5:beauty)
212         gener->SetFlavor(0);
213 //     
214         gener->SetOrigin(0., 0.0 ,0);
215         gener->SetSigma(0,0,5.3);
216         gener->SetVertexSmear(kPerEvent); 
217 // no tracking
218         gener->SetTrackingFlag(0);
219         gGener = gener;
220       }
221       break;
222      
223     case pythia:
224 //********************************************
225 // Example for Charm  Production with Pythia *
226 //********************************************
227       {
228         AliGenPythia *gener = new AliGenPythia(-1);
229 //   final state kinematic cuts
230         gener->SetMomentumRange(0,999);
231         gener->SetPhiRange(0. ,360.);
232         gener->SetThetaRange(0., 180.);
233         gener->SetYRange(-10,10);
234         gener->SetPtRange(0,100);
235 //   vertex position and smearing 
236         gener->SetOrigin(0,0,0);       // vertex position
237         gener->SetVertexSmear(kPerEvent);
238         gener->SetSigma(0,0,5.6);      // Sigma in (X,Y,Z) (cm) on IP position
239 //   Structure function. See the list in EVGEN/AliStructFuncType.h
240         gener->SetStrucFunc(kGRVHO);
241 // Select corection for nuclear structure functions
242 //     gener->SetNuclei(208,208);
243 //
244 //   Process type. See the list in PYTHIA6/AliPythia.h
245         gener->SetProcess(kPyBeauty);
246 //   
247 //   Pt transfer of the hard scattering
248         gener->SetPtHard(0.,5.);
249 //   Decay type (semielectronic, semimuonic, nodecay)
250         gener->SetForceDecay(kSemiElectronic);
251 //   Centre of mass energy 
252         gener->SetEnergyCMS(5500.);
253 //   No Tracking 
254         gener->SetTrackingFlag(0);
255         gGener = gener;
256       }
257       break;              
258
259     case herwig:
260 //********************************************
261 // Example for Charm  Production with Pythia *
262 //********************************************
263       {
264         AliGenHerwig *gener = new AliGenHerwig(-1);
265 //   final state kinematic cuts
266         gener->SetMomentumRange(0,7000);
267         gener->SetPhiRange(0. ,360.);
268         gener->SetThetaRange(0., 180.);
269         gener->SetYRange(-10,10);
270         gener->SetPtRange(0,7000);
271 //   vertex position and smearing 
272         gener->SetOrigin(0,0,0);       // vertex position
273         gener->SetVertexSmear(kPerEvent);
274         gener->SetSigma(0,0,5.6);      // Sigma in (X,Y,Z) (cm) on IP position
275 //   Beam momenta
276         gener->SetBeamMomenta(7000,7000);
277 //   Beams
278         gener->SetProjectile("P");
279         gener->SetTarget("P");
280 //   Structure function
281         gener->SetStrucFunc(kGRVHO);
282 //   Hard scatering
283         gener->SetPtHardMin(200);
284         gener->SetPtRMS(20);
285 //   Min bias
286         gener->SetProcess(8000);
287 //   No Tracking 
288         gener->SetTrackingFlag(0);
289         gGener = gener;
290       }
291       break;              
292
293     case mcatnlo:
294     {
295       AliGenHerwig *gener = new AliGenHerwig(-1);
296       gener->SetMomentumRange(0,7000);
297       gener->SetPhiRange(0. ,360.);
298       gener->SetThetaRange(0., 180.);
299       gener->SetYRange(-10,10);
300       gener->SetPtRange(0,7000);
301       gener->SetOrigin(0,0,0);       // vertex position
302       gener->SetVertexSmear(kPerEvent);
303       gener->SetSigma(0,0,5.6);      // Sigma in (X,Y,Z) (cm) on IP position
304       gener->SetBeamMomenta(7000,7000);
305       gener->SetProjectile("P");
306       gener->SetTarget("P");
307       gener->SetStrucFunc(kCTEQ5M);
308       gener->SetProcess(-1705);
309       gener->SetHardProcessFile("sb.events");
310       gener->SetEventListRange(0,1);
311       gener->SetTrackingFlag(0);
312       gGener = gener;
313     }
314     break;  
315       
316     case param1:
317 //*******************************************************
318 // Example for J/psi  Production from  Parameterisation 
319 // using default library (AliMUONlib)                                       
320 //*******************************************************
321       {
322         AliGenParam *gener =
323           new AliGenParam(ntracks, AliGenMUONlib::kUpsilon);
324         gener->SetMomentumRange(0,999);
325         gener->SetPtRange(0,999);     
326         gener->SetPhiRange(0. , 360.);
327         gener->SetYRange(2.5,4);
328         gener->SetCutOnChild(1);
329         gener->SetChildThetaRange(2,9);
330         gener->SetOrigin(0,0,0);        //vertex position
331         gener->SetSigma(0,0,5.3);       //Sigma in (X,Y,Z) (cm) on IP position
332         gener->SetForceDecay(kDiMuon);
333         gener->SetTrackingFlag(0);
334         gGener = gener;
335       }
336       break;
337
338     case param2:
339 //*******************************************************
340 // Example for Omega  Production from  Parameterisation 
341 // specifying library.                                       
342 //*******************************************************
343       {
344         AliGenParam *gener = new AliGenParam(1000,new AliGenPHOSlib(), 
345                                              AliGenPHOSlib::kOmega);
346         gener->SetWeighting(kNonAnalog);
347         gener->SetForceDecay(kNoDecay);
348         gener->SetPtRange(0,100);
349         gener->SetThetaRange(45,135);
350         gener->SetTrackingFlag(0);
351         gGener = gener;
352       }
353       break;
354
355     case param3:
356 //*******************************************************
357 // Example for Upsilon  Production from  Parameterisation 
358 // specifying library.                                       
359 // GSI style
360 //*******************************************************
361       {
362         AliGenParam *gener = new AliGenParam(1000,new AliGenGSIlib(), 
363                                              AliGenGSIlib::kUpsilon, "MUON");
364         gener->SetMomentumRange(0,999);
365         gener->SetPtRange(0,999);     
366         gener->SetPhiRange(0., 360.);
367         gener->SetYRange(2.5,4);
368         gener->SetCutOnChild(1);
369         gener->SetChildThetaRange(2,9);
370         gener->SetOrigin(0,0,0);        //vertex position
371         gener->SetSigma(0,0,5.3);       //Sigma in (X,Y,Z) (cm) on IP position
372         gener->SetForceDecay(kDiMuon);
373         gener->SetTrackingFlag(0);
374         gGener = gener;
375       }
376       break;
377      
378     case param4:
379 //*******************************************************
380 // Example for Omega  Production from  Parameterisation 
381 // specifying library.
382 // The alternative way.                                       
383 //*******************************************************
384       {
385         AliGenLib* Lib=new AliGenPHOSlib();
386         Int_t iOmega = AliGenPHOSlib::kOmega;
387         AliGenParam *gener = new AliGenParam(50, iOmega,            
388                                              Lib->GetPt(iOmega, ""),
389                                              Lib->GetY (iOmega, ""),
390                                              Lib->GetIp(iOmega, ""));
391         gener->SetPtRange(0,999);     
392         gener->SetWeighting(kNonAnalog);
393         gener->SetForceDecay(kNoDecay);
394         gener->SetTrackingFlag(0);
395         gGener = gener;
396       }
397       break;
398       
399     case fluka:
400 //*******************************************************
401 // Example for a FLUKA Boundary Source                  *
402 //*******************************************************
403       {
404         AliGenFLUKAsource *gener = new AliGenFLUKAsource(-1);
405         gener->SetFileName("$(ALICE_ROOT)/data/all32.root"); 
406         gener->SetPartFlag(9);
407         gener->SetAgeMax(1.e-5);
408 //  31.7 events     
409         gener->SetFraction(0.0315);     
410 //     gener->SetFraction(0.75*0.0315);     
411         rl->CdGAFile();
412 //     gener->SetPartFlag(10);
413         gener->SetMomentumRange(0,999);
414         gener->SetPhiRange(0.,360.);
415         gener->SetThetaRange(0., 180.); 
416         gener->SetAgeMax(1.e-5);
417      
418 //  31.7 events     
419 //     gener->SetFraction(0.0315);     
420         gGener = gener;
421       }
422       break;
423
424     case ntuple:
425 //*******************************************************
426 // Example for reading from a external file                  *
427 //*******************************************************
428       {
429         AliGenExtFile *gener = new AliGenExtFile(-1); 
430         gener->SetVertexSmear(kPerEvent); 
431         gener->SetTrackingFlag(1);
432         
433         AliGenReaderTreeK * reader = new AliGenReaderTreeK();
434         reader->SetFileName("$(ALICE_ROOT)/data/dtujet93.root");
435         gener->SetReader(reader);
436         gGener = gener;
437       }
438       break;
439
440     case halo:
441 //*******************************************************
442 // Example for Tunnel Halo Source                       *
443 //*******************************************************
444       {
445         AliGenHalo *gener = new AliGenHalo(ntracks); 
446         gener->SetFileName("/h1/morsch/marsip/marsip5.mu");
447         gGener = gener;
448       }
449       break;
450       
451     case cocktail:
452 //*******************************************************
453 // Example for a Cocktail                               *
454 //*******************************************************
455       {
456         AliGenCocktail *gener = new AliGenCocktail(); 
457
458         gener->SetPhiRange(0,360);
459         gener->SetYRange(2.5,4);
460         gener->SetThetaRange(2,9);
461         gener->SetPtRange(0,10);
462         gener->SetOrigin(0,0,0);        //vertex position
463         gener->SetSigma(0,0,0);         //Sigma in (X,Y,Z) (cm) on IP position
464         gener->SetMomentumRange(0,999);
465
466         AliGenParam *jpsi = new AliGenParam(1,AliGenMUONlib::kJpsi);
467         jpsi->SetForceDecay(kDiMuon);
468         jpsi->SetCutOnChild(1);
469
470      
471         AliGenFLUKAsource *bg = new AliGenFLUKAsource(-1);
472         bg->AddFile("$(ALICE_ROOT)/data/all32.root"); 
473         rl->CdGAFile();
474         bg->SetPartFlag(9);
475         bg->SetAgeMax(1.e-5);
476 //  31.7 events     
477 //     gener->SetFraction(0.0315);     
478         bg->SetFraction(0.01*0.0315);     
479       
480         gener->AddGenerator(jpsi,"J/Psi", 1);
481         gener->AddGenerator(bg,"Background",1);
482
483         gGener = gener;
484       }
485       break;
486     }
487  
488 // Activate this line if you want the vertex smearing to happen
489 // track by track
490 //
491 // gener->SetVertexSmear(kPerTrack); 
492
493   gGener->Init();
494
495   gAlice->SetField(-999,2);    //Specify maximum magnetic field in Tesla (neg. ==> default field)
496
497   Int_t iMAG=1;
498   rl->CdGAFile();
499
500 //=================== Alice BODY parameters =============================
501   AliBODY *BODY = new AliBODY("BODY","Alice envelop");
502
503
504   if(iMAG) {
505 //=================== MAG parameters ============================
506 // --- Start with Magnet since detector layouts may be depending ---
507 // --- on the selected Magnet dimensions ---
508     AliMAG *MAG  = new AliMAG("MAG","Magnet");
509   }
510 }