In vmctest/production:
[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/AliITSv11Hybrid.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,
69     QGSP_BERT_EMV_OPTICAL, CHIPS_OPTICAL, QGSP_BERT_CHIPS_OPTICAL, QGSP_FTFP_BERT_OPTICAL, kListMax
70   };
71
72 const char * physicsListName[] = {
73   "QGSP_BERT_EMV",         "CHIPS",         "QGSP_BERT_CHIPS",         "QGSP_FTFP_BERT",
74   "QGSP_BERT_EMV_OPTICAL", "CHIPS_OPTICAL", "QGSP_BERT_CHIPS_OPTICAL", "QGSP_FTFP_BERT_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 AliITSv11Hybrid("ITS","ITS v11Hybrid");
249       AliITS *ITS  = new AliITSv11("ITS","ITS v11");
250     }
251
252   if (iTPC)
253     {
254       //============================ TPC parameters =====================
255
256       AliTPC *TPC = new AliTPCv2("TPC", "Default");
257       TPC->SetPrimaryIonisation();// not used with Geant3
258     }
259
260
261   if (iTOF) {
262     //=================== TOF parameters ============================
263
264     AliTOF *TOF = new AliTOFv6T0("TOF", "normal TOF");
265   }
266
267
268   if (iHMPID)
269     {
270       //=================== HMPID parameters ===========================
271
272       AliHMPID *HMPID = new AliHMPIDv3("HMPID", "normal HMPID");
273
274     }
275
276
277   if (iZDC)
278     {
279       //=================== ZDC parameters ============================
280         
281       AliZDC *ZDC = new AliZDCv3("ZDC", "normal ZDC");
282       //Collimators aperture
283       ZDC->SetVCollSideCAperture(0.85);
284       ZDC->SetVCollSideCCentre(0.);
285       ZDC->SetVCollSideAAperture(0.75);
286       ZDC->SetVCollSideACentre(0.);
287       //Detector position
288       ZDC->SetYZNC(1.6);
289       ZDC->SetYZNA(1.6);
290       ZDC->SetYZPC(1.6);
291       ZDC->SetYZPA(1.6);
292     }
293
294   if (iTRD)
295     {
296       //=================== TRD parameters ============================
297
298       AliTRD *TRD = new AliTRDv1("TRD", "TRD slow simulator");
299       AliTRDgeometry *geoTRD = TRD->GetGeometry();
300       // Partial geometry: modules at 0,1,7,8,9,16,17
301       // starting at 3h in positive direction
302       geoTRD->SetSMstatus(2,0);
303       geoTRD->SetSMstatus(3,0);
304       geoTRD->SetSMstatus(4,0);
305       geoTRD->SetSMstatus(5,0);
306       geoTRD->SetSMstatus(6,0);
307       geoTRD->SetSMstatus(11,0);
308       geoTRD->SetSMstatus(12,0);
309       geoTRD->SetSMstatus(13,0);
310       geoTRD->SetSMstatus(14,0);
311       geoTRD->SetSMstatus(15,0);
312       geoTRD->SetSMstatus(16,0);
313     }
314
315   if (iFMD)
316     {
317       //=================== FMD parameters ============================
318
319       AliFMD *FMD = new AliFMDv1("FMD", "normal FMD");
320     }
321
322   if (iMUON)
323     {
324       //=================== MUON parameters ===========================
325       // New MUONv1 version (geometry defined via builders)
326       AliMUON *MUON = new AliMUONv1("MUON", "default");
327       // activate trigger efficiency by cells
328       MUON->SetTriggerEffCells(1);
329     }
330
331   if (iPHOS)
332     {
333       //=================== PHOS parameters ===========================
334
335       AliPHOS *PHOS = new AliPHOSv1("PHOS", "noCPV_Modules123");
336
337     }
338
339
340   if (iPMD)
341     {
342       //=================== PMD parameters ============================
343
344       AliPMD *PMD = new AliPMDv1("PMD", "normal PMD");
345     }
346
347   if (iT0)
348     {
349       //=================== T0 parameters ============================
350       AliT0 *T0 = new AliT0v1("T0", "T0 Detector");
351     }
352
353   if (iEMCAL)
354     {
355       //=================== EMCAL parameters ============================
356
357       AliEMCAL *EMCAL = new AliEMCALv2("EMCAL", "EMCAL_FIRSTYEARV1");
358     }
359
360   if (iACORDE)
361     {
362       //=================== ACORDE parameters ============================
363
364       AliACORDE *ACORDE = new AliACORDEv1("ACORDE", "normal ACORDE");
365     }
366
367   if (iVZERO)
368     {
369       //=================== ACORDE parameters ============================
370
371       AliVZERO *VZERO = new AliVZEROv7("VZERO", "normal VZERO");
372     }
373
374
375   // Load Geant4 + Geant4 VMC libraries
376   //
377   if (gClassTable->GetID("TGeant4") == -1) {
378     TString g4libsMacro = "$G4INSTALL/macro/g4libs.C";
379     //TString g4libsMacro = "$ALICE/geant4_vmc/examples/macro/g4libs.C";
380     //Load Geant4 libraries 
381     if (!gInterpreter->IsLoaded(g4libsMacro.Data())) {
382       gROOT->LoadMacro(g4libsMacro.Data());
383       gInterpreter->ProcessLine("g4libs()");
384     }
385   }    
386
387
388   // Create Geant4 VMC
389   //  
390   TGeant4 *geant4 = 0;
391   if ( ! gMC ) {
392     TG4RunConfiguration* runConfiguration=0x0;
393     for (Int_t iList = 0; iList < kListMax; iList++) {
394       if(iList<kListMax/2){
395         if(physicslist == iList){
396           runConfiguration = 
397             new TG4RunConfiguration("geomRoot", 
398                                     physicsListName[iList], 
399                                     "specialCuts+stackPopper+stepLimiter",
400                                     true);
401         }
402       }
403       else if(iList>=kListMax/2){//add "optical" PL to HadronPhysicsList
404         if(physicslist == iList){
405           runConfiguration = 
406             new TG4RunConfiguration("geomRoot", 
407                                     Form("%s+optical",physicsListName[iList-kListMax/2]), 
408                                     "specialCuts+stackPopper+stepLimiter",
409                                     true);
410         }
411       }
412     }
413     geant4 = new TGeant4("TGeant4", "The Geant4 Monte Carlo", runConfiguration);
414     cout << "Geant4 has been created." << endl;
415   } 
416   else {
417     cout << "Monte Carlo has been already created." << endl;
418   }  
419   
420   
421   
422   // Customization of Geant4 VMC
423   //
424   geant4->ProcessGeantCommand("/mcVerbose/all 1");  
425   geant4->ProcessGeantCommand("/mcVerbose/geometryManager 1");  
426   geant4->ProcessGeantCommand("/mcVerbose/opGeometryManager 1");  
427   geant4->ProcessGeantCommand("/mcTracking/loopVerbose 1");     
428   geant4->ProcessGeantCommand("/mcPhysics/rangeCuts 0.01 mm"); 
429   // for Geant4 <= 9.4.p03
430   //geant4->ProcessGeantCommand("/mcPhysics/selectOpProcess Scintillation");
431   //geant4->ProcessGeantCommand("/mcPhysics/setOpProcessActivation false");
432   // for Geant4 >= 9.5
433   geant4->ProcessGeantCommand("/optics_engine/selectOpProcess Scintillation");
434   geant4->ProcessGeantCommand("/optics_engine/setOpProcessUse false");
435   geant4->ProcessGeantCommand("/optics_engine/selectOpProcess OpWLS");
436   geant4->ProcessGeantCommand("/optics_engine/setOpProcessUse false");
437   geant4->ProcessGeantCommand("/optics_engine/selectOpProcess OpMieHG");
438   geant4->ProcessGeantCommand("/optics_engine/setOpProcessUse false");
439   
440   geant4->ProcessGeantCommand("/mcVerbose/composedPhysicsList 2");  
441   geant4->ProcessGeantCommand("/mcTracking/skipNeutrino true");
442   // geant4->ProcessGeantCommand("/mcDet/setMaxStepInLowDensityMaterials 1 cm");
443
444
445   //
446   //=======================================================================
447   // ************* STEERING parameters FOR ALICE SIMULATION **************
448   // --- Specify event type to be tracked through the ALICE setup
449   // --- All positions are in cm, angles in degrees, and P and E in GeV
450
451
452   gMC->SetProcess("DCAY",1);
453   gMC->SetProcess("PAIR",1);
454   gMC->SetProcess("COMP",1);
455   gMC->SetProcess("PHOT",1);
456   gMC->SetProcess("PFIS",0);
457   gMC->SetProcess("DRAY",0);
458   gMC->SetProcess("ANNI",1);
459   gMC->SetProcess("BREM",1);
460   gMC->SetProcess("MUNU",1);
461   gMC->SetProcess("CKOV",1);
462   gMC->SetProcess("HADR",1);
463   gMC->SetProcess("LOSS",2);
464   gMC->SetProcess("MULS",1);
465   gMC->SetProcess("RAYL",1);
466
467   Float_t cut = 1.e-3;        // 1MeV cut by default
468   Float_t tofmax = 1.e10;
469
470   gMC->SetCut("CUTGAM", cut);
471   gMC->SetCut("CUTELE", cut);
472   gMC->SetCut("CUTNEU", cut);
473   gMC->SetCut("CUTHAD", cut);
474   gMC->SetCut("CUTMUO", cut);
475   gMC->SetCut("BCUTE",  cut); 
476   gMC->SetCut("BCUTM",  cut); 
477   gMC->SetCut("DCUTE",  cut); 
478   gMC->SetCut("DCUTM",  cut); 
479   gMC->SetCut("PPCUTM", cut);
480   gMC->SetCut("TOFMAX", tofmax); 
481
482
483
484
485   //======================//
486   // Set External decayer //
487   //======================//
488   TVirtualMCDecayer* decayer = new AliDecayerPythia();
489   decayer->SetForceDecay(kAll);
490   decayer->Init();
491   gMC->SetExternalDecayer(decayer);
492
493   //=========================//
494   // Generator Configuration //
495   //=========================//
496   AliGenerator* gener = 0x0;
497   
498   if (proc == kPythia6) {
499     gener = MbPythia();
500   } else if (proc == kPythia6D6T) {
501     gener = MbPythiaTuneD6T();
502   } else if (proc == kPythia6ATLAS) {
503     gener = MbPythiaTuneATLAS();
504   } else if (proc == kPythiaPerugia0) {
505     gener = MbPythiaTunePerugia0();
506   } else if (proc == kPythia6ATLAS_Flat) {
507     gener = MbPythiaTuneATLAS_Flat();
508   } else if (proc == kPhojet) {
509     gener = MbPhojet();
510   }
511   
512   
513   //
514   //
515   // Size of the interaction diamond
516   // Longitudinal
517   Float_t sigmaz  = 5.4 / TMath::Sqrt(2.); // [cm]
518   if (energy == 900)
519     //sigmaz  = 10.5 / TMath::Sqrt(2.); // [cm]
520     //sigmaz = 3.7;
521     if (energy == 7000)
522       sigmaz  = 6.3 / TMath::Sqrt(2.); // [cm]
523   
524   //
525   // Transverse
526
527   // beta*
528   Float_t betast                  = 10.0;  // beta* [m]
529   if (runNumber >= 117048) betast =  2.0;
530   if (runNumber >  122375) betast =  3.5;  // starting with fill 1179
531   //    
532   //
533   Float_t eps     = 5.0e-6;            // emittance [m]
534   Float_t gamma   = energy / 2.0 / 0.938272;  // relativistic gamma [1]
535   Float_t sigmaxy = TMath::Sqrt(eps * betast / gamma) / TMath::Sqrt(2.) * 100.;  // [cm]
536   printf("\n \n Diamond size x-y: %10.3e z: %10.3e\n \n", sigmaxy, sigmaz);
537     
538   gener->SetSigma(sigmaxy, sigmaxy, sigmaz);      // Sigma in (X,Y,Z) (cm) on IP position
539   gener->SetVertexSmear(kPerEvent);
540   gener->Init();
541
542   printf("\n \n Comment: %s \n \n", comment.Data());
543
544
545 }
546 //
547 //           PYTHIA
548 //
549
550 AliGenerator* MbPythia()
551 {
552   comment = comment.Append(" pp: Pythia low-pt");
553   //
554   //    Pythia
555   AliGenPythia* pythia = new AliGenPythia(-1); 
556   pythia->SetMomentumRange(0, 999999.);
557   pythia->SetThetaRange(0., 180.);
558   pythia->SetYRange(-12.,12.);
559   pythia->SetPtRange(0,1000.);
560   pythia->SetProcess(kPyMb);
561   pythia->SetEnergyCMS(energy);
562       
563   return pythia;
564 }
565
566 AliGenerator* MbPythiaTuneD6T()
567 {
568   comment = comment.Append(" pp: Pythia low-pt");
569   //
570   //    Pythia
571   AliGenPythia* pythia = new AliGenPythia(-1); 
572   pythia->SetMomentumRange(0, 999999.);
573   pythia->SetThetaRange(0., 180.);
574   pythia->SetYRange(-12.,12.);
575   pythia->SetPtRange(0,1000.);
576   pythia->SetProcess(kPyMb);
577   pythia->SetEnergyCMS(energy);
578   //    Tune
579   //    109     D6T : Rick Field's CDF Tune D6T (NB: needs CTEQ6L pdfs externally)
580   pythia->SetTune(109); // F I X 
581   pythia->SetStrucFunc(kCTEQ6l);
582   //
583   return pythia;
584 }
585
586 AliGenerator* MbPythiaTunePerugia0()
587 {
588   comment = comment.Append(" pp: Pythia low-pt (Perugia0)");
589   //
590   //    Pythia
591   AliGenPythia* pythia = new AliGenPythia(-1); 
592   pythia->SetMomentumRange(0, 999999.);
593   pythia->SetThetaRange(0., 180.);
594   pythia->SetYRange(-12.,12.);
595   pythia->SetPtRange(0,1000.);
596   pythia->SetProcess(kPyMb);
597   pythia->SetEnergyCMS(energy);
598   //    Tune
599   //    320     Perugia 0
600   pythia->SetTune(320); 
601   pythia->UseNewMultipleInteractionsScenario();
602   pythia->SetCrossingAngle(0,0.000280);
603
604   //
605   return pythia;
606 }
607
608
609 AliGenerator* MbPythiaTuneATLAS()
610 {
611   comment = comment.Append(" pp: Pythia low-pt");
612   //
613   //    Pythia
614   AliGenPythia* pythia = new AliGenPythia(-1); 
615   pythia->SetMomentumRange(0, 999999.);
616   pythia->SetThetaRange(0., 180.);
617   pythia->SetYRange(-12.,12.);
618   pythia->SetPtRange(0,1000.);
619   pythia->SetProcess(kPyMb);
620   pythia->SetEnergyCMS(energy);
621   //    Tune
622   //    C   306 ATLAS-CSC: Arthur Moraes' (new) ATLAS tune (needs CTEQ6L externally)
623   pythia->SetTune(306);
624   pythia->SetStrucFunc(kCTEQ6l);
625   //
626   return pythia;
627 }
628
629 AliGenerator* MbPythiaTuneATLAS_Flat()
630 {
631   AliGenPythia* pythia = MbPythiaTuneATLAS();
632       
633   comment = comment.Append("; flat multiplicity distribution");
634       
635   // set high multiplicity trigger
636   // this weight achieves a flat multiplicity distribution
637   Double_t cont[] =
638     {0, 
639      0.234836, 0.103484, 0.00984802, 0.0199906, 0.0260018, 0.0208481, 0.0101477, 0.00146998, -0.00681513, -0.0114928,
640      -0.0113352, -0.0102012, -0.00895238, -0.00534961, -0.00261648, -0.000819048, 0.00130816, 0.00177978, 0.00373838, 0.00566255,
641      0.00628156, 0.00687458, 0.00668017, 0.00702917, 0.00810848, 0.00876167, 0.00832783, 0.00848518, 0.0107709, 0.0106849,
642      0.00933702, 0.00912525, 0.0106553, 0.0102785, 0.0101756, 0.010962, 0.00957103, 0.00970448, 0.0117133, 0.012271,
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.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113,
649      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
650      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
651      0};
652
653   Double_t err[] =
654     {0, 
655      0.00168216, 0.000743548, 0.00103125, 0.00108605, 0.00117101, 0.00124577, 0.00129119, 0.00128341, 0.00121421, 0.00112431,
656      0.00100588, 0.000895078, 0.000790314, 0.000711673, 0.000634547, 0.000589133, 0.000542763, 0.000509552, 0.000487375, 0.000468906, 
657      0.000460196, 0.000455806, 0.00044843, 0.000449317, 0.00045007, 0.000458016, 0.000461167, 0.000474742, 0.00050161, 0.00051637, 
658      0.000542336, 0.000558854, 0.000599169, 0.000611982, 0.000663855, 0.000696322, 0.000722825, 0.000771686, 0.000838023, 0.000908317, 
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.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003,
665      0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
666      0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
667      0};
668
669   TH1F *weight = new TH1F("newweight","newweight",120,-0.5,119.5);
670
671   weight->SetContent(cont);
672   weight->SetError(err);
673         
674   Int_t limit = weight->GetRandom();
675   pythia->SetTriggerChargedMultiplicity(limit, 1.4);
676       
677   comment = comment.Append(Form("; multiplicity threshold set to %d in |eta| < 1.4", limit));
678
679   return pythia;
680 }
681
682 AliGenerator* MbPhojet()
683 {
684   comment = comment.Append(" pp: Pythia low-pt");
685   //
686   //    DPMJET
687 #if defined(__CINT__)
688   gSystem->Load("libdpmjet");      // Parton density functions
689   gSystem->Load("libTDPMjet");      // Parton density functions
690 #endif
691   AliGenDPMjet* dpmjet = new AliGenDPMjet(-1); 
692   dpmjet->SetMomentumRange(0, 999999.);
693   dpmjet->SetThetaRange(0., 180.);
694   dpmjet->SetYRange(-12.,12.);
695   dpmjet->SetPtRange(0,1000.);
696   dpmjet->SetProcess(kDpmMb);
697   dpmjet->SetEnergyCMS(energy);
698   dpmjet->SetCrossingAngle(0,0.000280);
699   return dpmjet;
700 }
701
702 void ProcessEnvironmentVars()
703 {
704   // Run type
705   if (gSystem->Getenv("CONFIG_RUN_TYPE")) {
706     for (Int_t iRun = 0; iRun < kRunMax; iRun++) {
707       if (strcmp(gSystem->Getenv("CONFIG_RUN_TYPE"), pprRunName[iRun])==0) {
708         proc = (PDC06Proc_t)iRun;
709         cout<<"Run type set to "<<pprRunName[iRun]<<endl;
710       }
711     }
712   }
713
714   // Field
715   if (gSystem->Getenv("CONFIG_FIELD")) {
716     for (Int_t iField = 0; iField < kFieldMax; iField++) {
717       if (strcmp(gSystem->Getenv("CONFIG_FIELD"), pprField[iField])==0) {
718         mag = (Mag_t)iField;
719         cout<<"Field set to "<<pprField[iField]<<endl;
720       }
721     }
722   }
723
724   // Energy
725   if (gSystem->Getenv("CONFIG_ENERGY")) {
726     energy = atoi(gSystem->Getenv("CONFIG_ENERGY"));
727     cout<<"Energy set to "<<energy<<" GeV"<<endl;
728   }
729
730   // Random Number seed
731   if (gSystem->Getenv("CONFIG_SEED")) {
732     seed = atoi(gSystem->Getenv("CONFIG_SEED"));
733   }
734
735   // Run number
736   if (gSystem->Getenv("DC_RUN")) {
737     runNumber = atoi(gSystem->Getenv("DC_RUN"));
738   }
739
740   // Physics lists
741   if (gSystem->Getenv("CONFIG_PHYSICSLIST")) {
742     for (Int_t iList = 0; iList < kListMax; iList++) {
743       if (strcmp(gSystem->Getenv("CONFIG_PHYSICSLIST"), physicsListName[iList])==0){
744         physicslist = (PhysicsList_t)iList;
745         cout<<"Physics list set to "<< physicsListName[iList]<<endl;
746       }
747     }
748   }
749
750 }