AliGenHerwig added
[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, 
32                 param1, param2, param3, param4, 
33                 cocktail, fluka, halo, ntuple, scan};
34
35 gentype_t gentype = herwig;
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 TGeant3("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::fgkDefaultEventFolderName,
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(13);             //GEANT particle type
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(5);              //GEANT particle type
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(5); 
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 param1:
294 //*******************************************************
295 // Example for J/psi  Production from  Parameterisation 
296 // using default library (AliMUONlib)                                       
297 //*******************************************************
298       {
299         AliGenParam *gener =
300           new AliGenParam(ntracks, AliGenMUONlib::kUpsilon);
301         gener->SetMomentumRange(0,999);
302         gener->SetPtRange(0,999);     
303         gener->SetPhiRange(0. , 360.);
304         gener->SetYRange(2.5,4);
305         gener->SetCutOnChild(1);
306         gener->SetChildThetaRange(2,9);
307         gener->SetOrigin(0,0,0);        //vertex position
308         gener->SetSigma(0,0,5.3);       //Sigma in (X,Y,Z) (cm) on IP position
309         gener->SetForceDecay(kDiMuon);
310         gener->SetTrackingFlag(0);
311         gGener = gener;
312       }
313       break;
314
315     case param2:
316 //*******************************************************
317 // Example for Omega  Production from  Parameterisation 
318 // specifying library.                                       
319 //*******************************************************
320       {
321         AliGenParam *gener = new AliGenParam(1000,new AliGenPHOSlib(), 
322                                              AliGenPHOSlib::kOmega);
323         gener->SetWeighting(kNonAnalog);
324         gener->SetForceDecay(kNoDecay);
325         gener->SetPtRange(0,100);
326         gener->SetThetaRange(45,135);
327         gener->SetTrackingFlag(0);
328         gGener = gener;
329       }
330       break;
331
332     case param3:
333 //*******************************************************
334 // Example for Upsilon  Production from  Parameterisation 
335 // specifying library.                                       
336 // GSI style
337 //*******************************************************
338       {
339         AliGenParam *gener = new AliGenParam(1000,new AliGenGSIlib(), 
340                                              AliGenGSIlib::kUpsilon, "MUON");
341         gener->SetMomentumRange(0,999);
342         gener->SetPtRange(0,999);     
343         gener->SetPhiRange(0., 360.);
344         gener->SetYRange(2.5,4);
345         gener->SetCutOnChild(1);
346         gener->SetChildThetaRange(2,9);
347         gener->SetOrigin(0,0,0);        //vertex position
348         gener->SetSigma(0,0,5.3);       //Sigma in (X,Y,Z) (cm) on IP position
349         gener->SetForceDecay(kDiMuon);
350         gener->SetTrackingFlag(0);
351         gGener = gener;
352       }
353       break;
354      
355     case param4:
356 //*******************************************************
357 // Example for Omega  Production from  Parameterisation 
358 // specifying library.
359 // The alternative way.                                       
360 //*******************************************************
361       {
362         AliGenLib* Lib=new AliGenPHOSlib();
363         Int_t iOmega = AliGenPHOSlib::kOmega;
364         AliGenParam *gener = new AliGenParam(50, iOmega,            
365                                              Lib->GetPt(iOmega, ""),
366                                              Lib->GetY (iOmega, ""),
367                                              Lib->GetIp(iOmega, ""));
368         gener->SetPtRange(0,999);     
369         gener->SetWeighting(kNonAnalog);
370         gener->SetForceDecay(kNoDecay);
371         gener->SetTrackingFlag(0);
372         gGener = gener;
373       }
374       break;
375       
376     case fluka:
377 //*******************************************************
378 // Example for a FLUKA Boundary Source                  *
379 //*******************************************************
380       {
381         AliGenFLUKAsource *gener = new AliGenFLUKAsource(-1);
382         gener->SetFileName("$(ALICE_ROOT)/data/all32.root"); 
383         gener->SetPartFlag(9);
384         gener->SetAgeMax(1.e-5);
385 //  31.7 events     
386         gener->SetFraction(0.0315);     
387 //     gener->SetFraction(0.75*0.0315);     
388         rl->CdGAFile();
389 //     gener->SetPartFlag(10);
390         gener->SetMomentumRange(0,999);
391         gener->SetPhiRange(0.,360.);
392         gener->SetThetaRange(0., 180.); 
393         gener->SetAgeMax(1.e-5);
394      
395 //  31.7 events     
396 //     gener->SetFraction(0.0315);     
397         gGener = gener;
398       }
399       break;
400
401     case ntuple:
402 //*******************************************************
403 // Example for reading from a external file                  *
404 //*******************************************************
405       {
406         AliGenExtFile *gener = new AliGenExtFile(-1); 
407         gener->SetVertexSmear(kPerEvent); 
408         gener->SetTrackingFlag(1);
409         
410         AliGenReaderTreeK * reader = new AliGenReaderTreeK();
411         reader->SetFileName("$(ALICE_ROOT)/data/dtujet93.root");
412         gener->SetReader(reader);
413         gGener = gener;
414       }
415       break;
416
417     case halo:
418 //*******************************************************
419 // Example for Tunnel Halo Source                       *
420 //*******************************************************
421       {
422         AliGenHalo *gener = new AliGenHalo(ntracks); 
423         gener->SetFileName("/h1/morsch/marsip/marsip5.mu");
424         gGener = gener;
425       }
426       break;
427       
428     case cocktail:
429 //*******************************************************
430 // Example for a Cocktail                               *
431 //*******************************************************
432       {
433         AliGenCocktail *gener = new AliGenCocktail(); 
434
435         gener->SetPhiRange(0,360);
436         gener->SetYRange(2.5,4);
437         gener->SetThetaRange(2,9);
438         gener->SetPtRange(0,10);
439         gener->SetOrigin(0,0,0);        //vertex position
440         gener->SetSigma(0,0,0);         //Sigma in (X,Y,Z) (cm) on IP position
441         gener->SetMomentumRange(0,999);
442
443         AliGenParam *jpsi = new AliGenParam(1,AliGenMUONlib::kJpsi);
444         jpsi->SetForceDecay(kDiMuon);
445         jpsi->SetCutOnChild(1);
446
447      
448         AliGenFLUKAsource *bg = new AliGenFLUKAsource(-1);
449         bg->AddFile("$(ALICE_ROOT)/data/all32.root"); 
450         rl->CdGAFile();
451         bg->SetPartFlag(9);
452         bg->SetAgeMax(1.e-5);
453 //  31.7 events     
454 //     gener->SetFraction(0.0315);     
455         bg->SetFraction(0.01*0.0315);     
456       
457         gener->AddGenerator(jpsi,"J/Psi", 1);
458         gener->AddGenerator(bg,"Background",1);
459
460         gGener = gener;
461       }
462       break;
463     }
464  
465 // Activate this line if you want the vertex smearing to happen
466 // track by track
467 //
468 // gener->SetVertexSmear(kPerTrack); 
469
470   gGener->Init();
471
472   gAlice->SetField(-999,2);    //Specify maximum magnetic field in Tesla (neg. ==> default field)
473
474   Int_t iMAG=1;
475   rl->CdGAFile();
476
477 //=================== Alice BODY parameters =============================
478   AliBODY *BODY = new AliBODY("BODY","Alice envelop");
479
480
481   if(iMAG) {
482 //=================== MAG parameters ============================
483 // --- Start with Magnet since detector layouts may be depending ---
484 // --- on the selected Magnet dimensions ---
485     AliMAG *MAG  = new AliMAG("MAG","Magnet");
486   }
487 }