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