a5ce006b39b7e26d81c5fb146b7ea980068c6b3c
[u/mrichter/AliRoot.git] / test / vmctest / production / Config.C
1 // $Id$
2 //
3 // Configuration for the Geant4 production 2010
4 // By E. Sicking, CERN
5 //
6
7 #if !defined(__CINT__) || defined(__MAKECINT__)
8 #include <Riostream.h>
9 #include <TRandom.h>  
10 #include <TDatime.h>
11 #include <TSystem.h>
12 #include <TVirtualMC.h>
13 #include <TGeant3TGeo.h>
14 #include "STEER/AliRunLoader.h"
15 #include "STEER/AliRun.h"
16 #include "STEER/AliConfig.h"
17 #include "PYTHIA6/AliDecayerPythia.h"
18 #include "PYTHIA6/AliGenPythia.h"
19 #include "TDPMjet/AliGenDPMjet.h"
20 #include "STEER/AliMagFCheb.h"
21 #include "STRUCT/AliBODY.h"
22 #include "STRUCT/AliMAG.h"
23 #include "STRUCT/AliABSOv3.h"
24 #include "STRUCT/AliDIPOv3.h"
25 #include "STRUCT/AliHALLv3.h"
26 #include "STRUCT/AliFRAMEv2.h"
27 #include "STRUCT/AliSHILv3.h"
28 #include "STRUCT/AliPIPEv3.h"
29 #include "ITS/AliITSv11.h"
30 #include "TPC/AliTPCv2.h"
31 #include "TOF/AliTOFv6T0.h"
32 #include "HMPID/AliHMPIDv3.h"
33 #include "ZDC/AliZDCv3.h"
34 #include "TRD/AliTRDv1.h"
35 #include "TRD/AliTRDgeometry.h"
36 #include "FMD/AliFMDv1.h"
37 #include "MUON/AliMUONv1.h"
38 #include "PHOS/AliPHOSv1.h"
39 #include "PHOS/AliPHOSSimParam.h"
40 #include "PMD/AliPMDv1.h"
41 #include "T0/AliT0v1.h"
42 #include "EMCAL/AliEMCALv2.h"
43 #include "ACORDE/AliACORDEv1.h"
44 #include "VZERO/AliVZEROv7.h"
45 #endif
46
47
48 enum PDC06Proc_t 
49   {
50     kPythia6, kPythia6D6T, kPythia6ATLAS, kPythia6ATLAS_Flat, kPythiaPerugia0, kPhojet, kRunMax
51   };
52
53 const char * pprRunName[] = {
54   "kPythia6", "kPythia6D6T", "kPythia6ATLAS", "kPythia6ATLAS_Flat", "kPythiaPerugia0", "kPhojet" 
55 };
56
57 enum Mag_t
58   {
59     kNoField, k5kG, kFieldMax
60   };
61
62 const char * pprField[] = {
63   "kNoField", "k5kG"
64 };
65
66 enum PhysicsList_t 
67   {
68     QGSP_BERT_EMV,         CHIPS,         QGSP_BERT_CHIPS,         QGSP_FTFP_BERT,         FTFP_BERT,         FTFP_BERT_EMV,
69     QGSP_BERT_EMV_OPTICAL, CHIPS_OPTICAL, QGSP_BERT_CHIPS_OPTICAL, QGSP_FTFP_BERT_OPTICAL, FTFP_BERT_OPTICAL, FTFP_BERT_EMV_OPTICAL, kListMax
70   };
71
72 const char * physicsListName[] = {
73   "QGSP_BERT_EMV",         "CHIPS",         "QGSP_BERT_CHIPS",         "QGSP_FTFP_BERT",         "FTFP_BERT",         "FTFP_BERT_EMV",
74   "QGSP_BERT_EMV_OPTICAL", "CHIPS_OPTICAL", "QGSP_BERT_CHIPS_OPTICAL", "QGSP_FTFP_BERT_OPTICAL", "FTFP_BERT_OPTICAL", "FTFP_BERT_EMV_OPTICAL",
75 };
76
77 enum PprTrigConf_t
78   {
79     kDefaultPPTrig, kDefaultPbPbTrig
80   };
81
82 const char * pprTrigConfName[] = {
83   "p-p","Pb-Pb"
84 };
85
86 //--- Functions ---
87 class AliGenPythia;
88 AliGenerator *MbPythia();
89 AliGenerator *MbPythiaTuneD6T();
90 AliGenerator *MbPhojet();
91 void ProcessEnvironmentVars();
92
93 // Geterator, field, beam energy
94 static PDC06Proc_t   proc     = kPhojet;
95 static Mag_t         mag      = k5kG;
96 static Float_t       energy   = 10000; // energy in CMS
97 static Int_t         runNumber = 0;
98 static PprTrigConf_t strig = kDefaultPPTrig; // default pp trigger configuration
99 static PhysicsList_t physicslist  = QGSP_BERT_EMV;
100
101 //========================//
102 // Set Random Number seed //
103 //========================//
104 TDatime dt;
105 static UInt_t seed    = dt.Get();
106
107 // Comment line
108 static TString comment;
109
110 void Config()
111 {
112     
113
114   // Get settings from environment variables
115   ProcessEnvironmentVars();
116
117   gRandom->SetSeed(seed);
118   cerr<<"Seed for random number generation= "<<seed<<endl; 
119
120   // Libraries required by geant321
121 #if defined(__CINT__)
122   gSystem->Load("liblhapdf");      // Parton density functions
123   gSystem->Load("libEGPythia6");   // TGenerator interface
124   if (proc == kPythia6 || proc == kPhojet) {
125     gSystem->Load("libpythia6");        // Pythia 6.2
126   } else {
127     gSystem->Load("libpythia6.4.21");   // Pythia 6.4
128   }
129   gSystem->Load("libAliPythia6");  // ALICE specific implementations
130   // gSystem->Load("libgeant321");
131 #endif
132
133   // new TGeant3TGeo("C++ Interface to Geant3");
134
135   //=======================================================================
136   //  Create the output file
137
138    
139   AliRunLoader* rl=0x0;
140
141   cout<<"Config.C: Creating Run Loader ..."<<endl;
142   rl = AliRunLoader::Open("galice.root",
143                           AliConfig::GetDefaultEventFolderName(),
144                           "recreate");
145   if (rl == 0x0)
146     {
147       gAlice->Fatal("Config.C","Can not instatiate the Run Loader");
148       return;
149     }
150   rl->SetCompressionLevel(2);
151   rl->SetNumberOfEventsPerFile(1000);
152   gAlice->SetRunLoader(rl);
153   // gAlice->SetGeometryFromFile("geometry.root");
154   // gAlice->SetGeometryFromCDB();
155   
156   // Set the trigger configuration: proton-proton
157
158   AliSimulation::Instance()->SetTriggerConfig(pprTrigConfName[strig]);
159   cout <<"Trigger configuration is set to  "<<pprTrigConfName[strig]<<endl;
160
161
162   rl->CdGAFile();
163   
164   Int_t iABSO  = 1;
165   Int_t iACORDE= 0;
166   Int_t iDIPO  = 1;
167   Int_t iEMCAL = 1;
168   Int_t iFMD   = 1;
169   Int_t iFRAME = 1;
170   Int_t iHALL  = 1;
171   Int_t iITS   = 1;
172   Int_t iMAG   = 1;
173   Int_t iMUON  = 1;
174   Int_t iPHOS  = 1;
175   Int_t iPIPE  = 1;
176   Int_t iPMD   = 1;
177   Int_t iHMPID = 1;
178   Int_t iSHIL  = 1;
179   Int_t iT0    = 1;
180   Int_t iTOF   = 1;
181   Int_t iTPC   = 1;
182   Int_t iTRD   = 1;
183   Int_t iVZERO = 1;
184   Int_t iZDC   = 1;
185   
186
187   //=================== Alice BODY parameters =============================
188   AliBODY *BODY = new AliBODY("BODY", "Alice envelop");
189
190
191   if (iMAG)
192     {
193       //=================== MAG parameters ============================
194       // --- Start with Magnet since detector layouts may be depending ---
195       // --- on the selected Magnet dimensions ---
196       AliMAG *MAG = new AliMAG("MAG", "Magnet");
197     }
198
199
200   if (iABSO)
201     {
202       //=================== ABSO parameters ============================
203       AliABSO *ABSO = new AliABSOv3("ABSO", "Muon Absorber");
204     }
205
206   if (iDIPO)
207     {
208       //=================== DIPO parameters ============================
209
210       AliDIPO *DIPO = new AliDIPOv3("DIPO", "Dipole version 3");
211     }
212
213   if (iHALL)
214     {
215       //=================== HALL parameters ============================
216
217       AliHALL *HALL = new AliHALLv3("HALL", "Alice Hall");
218     }
219
220
221   if (iFRAME)
222     {
223       //=================== FRAME parameters ============================
224
225       AliFRAMEv2 *FRAME = new AliFRAMEv2("FRAME", "Space Frame");
226       FRAME->SetHoles(1);
227     }
228
229   if (iSHIL)
230     {
231       //=================== SHIL parameters ============================
232
233       AliSHIL *SHIL = new AliSHILv3("SHIL", "Shielding Version 3");
234     }
235
236
237   if (iPIPE)
238     {
239       //=================== PIPE parameters ============================
240
241       AliPIPE *PIPE = new AliPIPEv3("PIPE", "Beam Pipe");
242     }
243  
244   if (iITS)
245     {
246       //=================== ITS parameters ============================
247
248       AliITS *ITS  = new AliITSv11("ITS","ITS v11");
249     }
250
251   if (iTPC)
252     {
253       //============================ TPC parameters =====================
254
255       AliTPC *TPC = new AliTPCv2("TPC", "Default");
256       TPC->SetPrimaryIonisation();// not used with Geant3
257     }
258
259
260   if (iTOF) {
261     //=================== TOF parameters ============================
262
263     AliTOF *TOF = new AliTOFv6T0("TOF", "normal TOF");
264   }
265
266
267   if (iHMPID)
268     {
269       //=================== HMPID parameters ===========================
270
271       AliHMPID *HMPID = new AliHMPIDv3("HMPID", "normal HMPID");
272
273     }
274
275
276   if (iZDC)
277     {
278       //=================== ZDC parameters ============================
279         
280       AliZDC *ZDC = new AliZDCv3("ZDC", "normal ZDC");
281       //Collimators aperture
282       ZDC->SetVCollSideCAperture(0.85);
283       ZDC->SetVCollSideCCentre(0.);
284       ZDC->SetVCollSideAAperture(0.75);
285       ZDC->SetVCollSideACentre(0.);
286       //Detector position
287       ZDC->SetYZNC(1.6);
288       ZDC->SetYZNA(1.6);
289       ZDC->SetYZPC(1.6);
290       ZDC->SetYZPA(1.6);
291     }
292
293   if (iTRD)
294     {
295       //=================== TRD parameters ============================
296
297       AliTRD *TRD = new AliTRDv1("TRD", "TRD slow simulator");
298       AliTRDgeometry *geoTRD = TRD->GetGeometry();
299       // Partial geometry: modules at 0,1,7,8,9,16,17
300       // starting at 3h in positive direction
301       geoTRD->SetSMstatus(2,0);
302       geoTRD->SetSMstatus(3,0);
303       geoTRD->SetSMstatus(4,0);
304       geoTRD->SetSMstatus(5,0);
305       geoTRD->SetSMstatus(6,0);
306       geoTRD->SetSMstatus(11,0);
307       geoTRD->SetSMstatus(12,0);
308       geoTRD->SetSMstatus(13,0);
309       geoTRD->SetSMstatus(14,0);
310       geoTRD->SetSMstatus(15,0);
311       geoTRD->SetSMstatus(16,0);
312     }
313
314   if (iFMD)
315     {
316       //=================== FMD parameters ============================
317
318       AliFMD *FMD = new AliFMDv1("FMD", "normal FMD");
319     }
320
321   if (iMUON)
322     {
323       //=================== MUON parameters ===========================
324       // New MUONv1 version (geometry defined via builders)
325       AliMUON *MUON = new AliMUONv1("MUON", "default");
326       // activate trigger efficiency by cells
327       MUON->SetTriggerEffCells(1);
328     }
329
330   if (iPHOS)
331     {
332       //=================== PHOS parameters ===========================
333
334       AliPHOS *PHOS = new AliPHOSv1("PHOS", "noCPV_Modules123");
335
336     }
337
338
339   if (iPMD)
340     {
341       //=================== PMD parameters ============================
342
343       AliPMD *PMD = new AliPMDv1("PMD", "normal PMD");
344     }
345
346   if (iT0)
347     {
348       //=================== T0 parameters ============================
349       AliT0 *T0 = new AliT0v1("T0", "T0 Detector");
350     }
351
352   if (iEMCAL)
353     {
354       //=================== EMCAL parameters ============================
355
356       AliEMCAL *EMCAL = new AliEMCALv2("EMCAL", "EMCAL_FIRSTYEARV1");
357     }
358
359   if (iACORDE)
360     {
361       //=================== ACORDE parameters ============================
362
363       AliACORDE *ACORDE = new AliACORDEv1("ACORDE", "normal ACORDE");
364     }
365
366   if (iVZERO)
367     {
368       //=================== ACORDE parameters ============================
369
370       AliVZERO *VZERO = new AliVZEROv7("VZERO", "normal VZERO");
371     }
372
373
374   // Load Geant4 + Geant4 VMC libraries
375   //
376   if (gClassTable->GetID("TGeant4") == -1) {
377     TString g4libsMacro = "$G4INSTALL/macro/g4libs.C";
378     //TString g4libsMacro = "$ALICE/geant4_vmc/examples/macro/g4libs.C";
379     //Load Geant4 libraries 
380     if (!gInterpreter->IsLoaded(g4libsMacro.Data())) {
381       gROOT->LoadMacro(g4libsMacro.Data());
382       gInterpreter->ProcessLine("g4libs()");
383     }
384   }    
385
386
387   // Create Geant4 VMC
388   //  
389   TGeant4 *geant4 = 0;
390   if ( ! gMC ) {
391     TG4RunConfiguration* runConfiguration=0x0;
392     for (Int_t iList = 0; iList < kListMax; iList++) {
393       if(iList<kListMax/2){
394         if(physicslist == iList){
395           runConfiguration = 
396             new TG4RunConfiguration("geomRoot", 
397                                     physicsListName[iList], 
398                                     "specialCuts+stackPopper+stepLimiter",
399                                     true);
400         }
401       }
402       else if(iList>=kListMax/2){//add "optical" PL to HadronPhysicsList
403         if(physicslist == iList){
404           runConfiguration = 
405             new TG4RunConfiguration("geomRoot", 
406                                     Form("%s+optical",physicsListName[iList-kListMax/2]), 
407                                     "specialCuts+stackPopper+stepLimiter",
408                                     true);
409         }
410       }
411     }
412     geant4 = new TGeant4("TGeant4", "The Geant4 Monte Carlo", runConfiguration);
413     cout << "Geant4 has been created." << endl;
414   } 
415   else {
416     cout << "Monte Carlo has been already created." << endl;
417   }  
418   
419   
420   
421   // Customization of Geant4 VMC
422   //
423   geant4->ProcessGeantCommand("/mcVerbose/all 1");  
424   geant4->ProcessGeantCommand("/mcVerbose/geometryManager 1");  
425   geant4->ProcessGeantCommand("/mcVerbose/opGeometryManager 1");  
426   geant4->ProcessGeantCommand("/mcTracking/loopVerbose 1");     
427   geant4->ProcessGeantCommand("/mcPhysics/rangeCuts 0.01 mm"); 
428   // for Geant4 <= 9.4.p03
429   //geant4->ProcessGeantCommand("/mcPhysics/selectOpProcess Scintillation");
430   //geant4->ProcessGeantCommand("/mcPhysics/setOpProcessActivation false");
431   // for Geant4 >= 9.5
432   geant4->ProcessGeantCommand("/optics_engine/selectOpProcess Scintillation");
433   geant4->ProcessGeantCommand("/optics_engine/setOpProcessUse false");
434   geant4->ProcessGeantCommand("/optics_engine/selectOpProcess OpWLS");
435   geant4->ProcessGeantCommand("/optics_engine/setOpProcessUse false");
436   geant4->ProcessGeantCommand("/optics_engine/selectOpProcess OpMieHG");
437   geant4->ProcessGeantCommand("/optics_engine/setOpProcessUse false");
438   
439   geant4->ProcessGeantCommand("/mcVerbose/composedPhysicsList 2");  
440   geant4->ProcessGeantCommand("/mcTracking/skipNeutrino true");
441   // geant4->ProcessGeantCommand("/mcDet/setMaxStepInLowDensityMaterials 1 cm");
442
443
444   //
445   //=======================================================================
446   // ************* STEERING parameters FOR ALICE SIMULATION **************
447   // --- Specify event type to be tracked through the ALICE setup
448   // --- All positions are in cm, angles in degrees, and P and E in GeV
449
450
451   gMC->SetProcess("DCAY",1);
452   gMC->SetProcess("PAIR",1);
453   gMC->SetProcess("COMP",1);
454   gMC->SetProcess("PHOT",1);
455   gMC->SetProcess("PFIS",0);
456   gMC->SetProcess("DRAY",0);
457   gMC->SetProcess("ANNI",1);
458   gMC->SetProcess("BREM",1);
459   gMC->SetProcess("MUNU",1);
460   gMC->SetProcess("CKOV",1);
461   gMC->SetProcess("HADR",1);
462   gMC->SetProcess("LOSS",2);
463   gMC->SetProcess("MULS",1);
464   gMC->SetProcess("RAYL",1);
465
466   Float_t cut = 1.e-3;        // 1MeV cut by default
467   Float_t tofmax = 1.e10;
468
469   gMC->SetCut("CUTGAM", cut);
470   gMC->SetCut("CUTELE", cut);
471   gMC->SetCut("CUTNEU", cut);
472   gMC->SetCut("CUTHAD", cut);
473   gMC->SetCut("CUTMUO", cut);
474   gMC->SetCut("BCUTE",  cut); 
475   gMC->SetCut("BCUTM",  cut); 
476   gMC->SetCut("DCUTE",  cut); 
477   gMC->SetCut("DCUTM",  cut); 
478   gMC->SetCut("PPCUTM", cut);
479   gMC->SetCut("TOFMAX", tofmax); 
480
481
482
483
484   //======================//
485   // Set External decayer //
486   //======================//
487   TVirtualMCDecayer* decayer = new AliDecayerPythia();
488   decayer->SetForceDecay(kAll);
489   decayer->Init();
490   gMC->SetExternalDecayer(decayer);
491
492   //=========================//
493   // Generator Configuration //
494   //=========================//
495   AliGenerator* gener = 0x0;
496   
497   if (proc == kPythia6) {
498     gener = MbPythia();
499   } else if (proc == kPythia6D6T) {
500     gener = MbPythiaTuneD6T();
501   } else if (proc == kPythia6ATLAS) {
502     gener = MbPythiaTuneATLAS();
503   } else if (proc == kPythiaPerugia0) {
504     gener = MbPythiaTunePerugia0();
505   } else if (proc == kPythia6ATLAS_Flat) {
506     gener = MbPythiaTuneATLAS_Flat();
507   } else if (proc == kPhojet) {
508     gener = MbPhojet();
509   }
510   
511   
512   //
513   //
514   // Size of the interaction diamond
515   // Longitudinal
516   Float_t sigmaz  = 5.4 / TMath::Sqrt(2.); // [cm]
517   if (energy == 900)
518     //sigmaz  = 10.5 / TMath::Sqrt(2.); // [cm]
519     //sigmaz = 3.7;
520     if (energy == 7000)
521       sigmaz  = 6.3 / TMath::Sqrt(2.); // [cm]
522   
523   //
524   // Transverse
525
526   // beta*
527   Float_t betast                  = 10.0;  // beta* [m]
528   if (runNumber >= 117048) betast =  2.0;
529   if (runNumber >  122375) betast =  3.5;  // starting with fill 1179
530   //    
531   //
532   Float_t eps     = 5.0e-6;            // emittance [m]
533   Float_t gamma   = energy / 2.0 / 0.938272;  // relativistic gamma [1]
534   Float_t sigmaxy = TMath::Sqrt(eps * betast / gamma) / TMath::Sqrt(2.) * 100.;  // [cm]
535   printf("\n \n Diamond size x-y: %10.3e z: %10.3e\n \n", sigmaxy, sigmaz);
536     
537   gener->SetSigma(sigmaxy, sigmaxy, sigmaz);      // Sigma in (X,Y,Z) (cm) on IP position
538   gener->SetVertexSmear(kPerEvent);
539   gener->Init();
540
541   printf("\n \n Comment: %s \n \n", comment.Data());
542
543
544 }
545 //
546 //           PYTHIA
547 //
548
549 AliGenerator* MbPythia()
550 {
551   comment = comment.Append(" pp: Pythia low-pt");
552   //
553   //    Pythia
554   AliGenPythia* pythia = new AliGenPythia(-1); 
555   pythia->SetMomentumRange(0, 999999.);
556   pythia->SetThetaRange(0., 180.);
557   pythia->SetYRange(-12.,12.);
558   pythia->SetPtRange(0,1000.);
559   pythia->SetProcess(kPyMb);
560   pythia->SetEnergyCMS(energy);
561       
562   return pythia;
563 }
564
565 AliGenerator* MbPythiaTuneD6T()
566 {
567   comment = comment.Append(" pp: Pythia low-pt");
568   //
569   //    Pythia
570   AliGenPythia* pythia = new AliGenPythia(-1); 
571   pythia->SetMomentumRange(0, 999999.);
572   pythia->SetThetaRange(0., 180.);
573   pythia->SetYRange(-12.,12.);
574   pythia->SetPtRange(0,1000.);
575   pythia->SetProcess(kPyMb);
576   pythia->SetEnergyCMS(energy);
577   //    Tune
578   //    109     D6T : Rick Field's CDF Tune D6T (NB: needs CTEQ6L pdfs externally)
579   pythia->SetTune(109); // F I X 
580   pythia->SetStrucFunc(kCTEQ6l);
581   //
582   return pythia;
583 }
584
585 AliGenerator* MbPythiaTunePerugia0()
586 {
587   comment = comment.Append(" pp: Pythia low-pt (Perugia0)");
588   //
589   //    Pythia
590   AliGenPythia* pythia = new AliGenPythia(-1); 
591   pythia->SetMomentumRange(0, 999999.);
592   pythia->SetThetaRange(0., 180.);
593   pythia->SetYRange(-12.,12.);
594   pythia->SetPtRange(0,1000.);
595   pythia->SetProcess(kPyMb);
596   pythia->SetEnergyCMS(energy);
597   //    Tune
598   //    320     Perugia 0
599   pythia->SetTune(320); 
600   pythia->UseNewMultipleInteractionsScenario();
601   pythia->SetCrossingAngle(0,0.000280);
602
603   //
604   return pythia;
605 }
606
607
608 AliGenerator* MbPythiaTuneATLAS()
609 {
610   comment = comment.Append(" pp: Pythia low-pt");
611   //
612   //    Pythia
613   AliGenPythia* pythia = new AliGenPythia(-1); 
614   pythia->SetMomentumRange(0, 999999.);
615   pythia->SetThetaRange(0., 180.);
616   pythia->SetYRange(-12.,12.);
617   pythia->SetPtRange(0,1000.);
618   pythia->SetProcess(kPyMb);
619   pythia->SetEnergyCMS(energy);
620   //    Tune
621   //    C   306 ATLAS-CSC: Arthur Moraes' (new) ATLAS tune (needs CTEQ6L externally)
622   pythia->SetTune(306);
623   pythia->SetStrucFunc(kCTEQ6l);
624   //
625   return pythia;
626 }
627
628 AliGenerator* MbPythiaTuneATLAS_Flat()
629 {
630   AliGenPythia* pythia = MbPythiaTuneATLAS();
631       
632   comment = comment.Append("; flat multiplicity distribution");
633       
634   // set high multiplicity trigger
635   // this weight achieves a flat multiplicity distribution
636   Double_t cont[] =
637     {0, 
638      0.234836, 0.103484, 0.00984802, 0.0199906, 0.0260018, 0.0208481, 0.0101477, 0.00146998, -0.00681513, -0.0114928,
639      -0.0113352, -0.0102012, -0.00895238, -0.00534961, -0.00261648, -0.000819048, 0.00130816, 0.00177978, 0.00373838, 0.00566255,
640      0.00628156, 0.00687458, 0.00668017, 0.00702917, 0.00810848, 0.00876167, 0.00832783, 0.00848518, 0.0107709, 0.0106849,
641      0.00933702, 0.00912525, 0.0106553, 0.0102785, 0.0101756, 0.010962, 0.00957103, 0.00970448, 0.0117133, 0.012271,
642      0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113,
643      0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113,
644      0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113,
645      0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113,
646      0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113,
647      0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113,
648      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
649      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
650      0};
651
652   Double_t err[] =
653     {0, 
654      0.00168216, 0.000743548, 0.00103125, 0.00108605, 0.00117101, 0.00124577, 0.00129119, 0.00128341, 0.00121421, 0.00112431,
655      0.00100588, 0.000895078, 0.000790314, 0.000711673, 0.000634547, 0.000589133, 0.000542763, 0.000509552, 0.000487375, 0.000468906, 
656      0.000460196, 0.000455806, 0.00044843, 0.000449317, 0.00045007, 0.000458016, 0.000461167, 0.000474742, 0.00050161, 0.00051637, 
657      0.000542336, 0.000558854, 0.000599169, 0.000611982, 0.000663855, 0.000696322, 0.000722825, 0.000771686, 0.000838023, 0.000908317, 
658      0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003,
659      0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003,
660      0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003,
661      0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003,
662      0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003,
663      0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003,
664      0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
665      0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
666      0};
667
668   TH1F *weight = new TH1F("newweight","newweight",120,-0.5,119.5);
669
670   weight->SetContent(cont);
671   weight->SetError(err);
672         
673   Int_t limit = weight->GetRandom();
674   pythia->SetTriggerChargedMultiplicity(limit, 1.4);
675       
676   comment = comment.Append(Form("; multiplicity threshold set to %d in |eta| < 1.4", limit));
677
678   return pythia;
679 }
680
681 AliGenerator* MbPhojet()
682 {
683   comment = comment.Append(" pp: Pythia low-pt");
684   //
685   //    DPMJET
686 #if defined(__CINT__)
687   gSystem->Load("libdpmjet");      // Parton density functions
688   gSystem->Load("libTDPMjet");      // Parton density functions
689 #endif
690   AliGenDPMjet* dpmjet = new AliGenDPMjet(-1); 
691   dpmjet->SetMomentumRange(0, 999999.);
692   dpmjet->SetThetaRange(0., 180.);
693   dpmjet->SetYRange(-12.,12.);
694   dpmjet->SetPtRange(0,1000.);
695   dpmjet->SetProcess(kDpmMb);
696   dpmjet->SetEnergyCMS(energy);
697   dpmjet->SetCrossingAngle(0,0.000280);
698   return dpmjet;
699 }
700
701 void ProcessEnvironmentVars()
702 {
703   // Run type
704   if (gSystem->Getenv("CONFIG_RUN_TYPE")) {
705     for (Int_t iRun = 0; iRun < kRunMax; iRun++) {
706       if (strcmp(gSystem->Getenv("CONFIG_RUN_TYPE"), pprRunName[iRun])==0) {
707         proc = (PDC06Proc_t)iRun;
708         cout<<"Run type set to "<<pprRunName[iRun]<<endl;
709       }
710     }
711   }
712
713   // Field
714   if (gSystem->Getenv("CONFIG_FIELD")) {
715     for (Int_t iField = 0; iField < kFieldMax; iField++) {
716       if (strcmp(gSystem->Getenv("CONFIG_FIELD"), pprField[iField])==0) {
717         mag = (Mag_t)iField;
718         cout<<"Field set to "<<pprField[iField]<<endl;
719       }
720     }
721   }
722
723   // Energy
724   if (gSystem->Getenv("CONFIG_ENERGY")) {
725     energy = atoi(gSystem->Getenv("CONFIG_ENERGY"));
726     cout<<"Energy set to "<<energy<<" GeV"<<endl;
727   }
728
729   // Random Number seed
730   if (gSystem->Getenv("CONFIG_SEED")) {
731     seed = atoi(gSystem->Getenv("CONFIG_SEED"));
732   }
733
734   // Run number
735   if (gSystem->Getenv("DC_RUN")) {
736     runNumber = atoi(gSystem->Getenv("DC_RUN"));
737   }
738
739   // Physics lists
740   if (gSystem->Getenv("CONFIG_PHYSICSLIST")) {
741     for (Int_t iList = 0; iList < kListMax; iList++) {
742       if (strcmp(gSystem->Getenv("CONFIG_PHYSICSLIST"), physicsListName[iList])==0){
743         physicslist = (PhysicsList_t)iList;
744         cout<<"Physics list set to "<< physicsListName[iList]<<endl;
745       }
746     }
747   }
748
749 }