d5a381f35dd3a1b841a44b0b12544d6e89592c4d
[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   geant4->ProcessGeantCommand("/mcPhysics/selectOpProcess Scintillation");
430   geant4->ProcessGeantCommand("/mcPhysics/setOpProcessActivation false");
431   geant4->ProcessGeantCommand("/mcVerbose/composedPhysicsList 2");  
432   geant4->ProcessGeantCommand("/mcTracking/skipNeutrino true");
433   // geant4->ProcessGeantCommand("/mcDet/setMaxStepInLowDensityMaterials 1 cm");
434
435
436   //
437   //=======================================================================
438   // ************* STEERING parameters FOR ALICE SIMULATION **************
439   // --- Specify event type to be tracked through the ALICE setup
440   // --- All positions are in cm, angles in degrees, and P and E in GeV
441
442
443   gMC->SetProcess("DCAY",1);
444   gMC->SetProcess("PAIR",1);
445   gMC->SetProcess("COMP",1);
446   gMC->SetProcess("PHOT",1);
447   gMC->SetProcess("PFIS",0);
448   gMC->SetProcess("DRAY",0);
449   gMC->SetProcess("ANNI",1);
450   gMC->SetProcess("BREM",1);
451   gMC->SetProcess("MUNU",1);
452   gMC->SetProcess("CKOV",1);
453   gMC->SetProcess("HADR",1);
454   gMC->SetProcess("LOSS",2);
455   gMC->SetProcess("MULS",1);
456   gMC->SetProcess("RAYL",1);
457
458   Float_t cut = 1.e-3;        // 1MeV cut by default
459   Float_t tofmax = 1.e10;
460
461   gMC->SetCut("CUTGAM", cut);
462   gMC->SetCut("CUTELE", cut);
463   gMC->SetCut("CUTNEU", cut);
464   gMC->SetCut("CUTHAD", cut);
465   gMC->SetCut("CUTMUO", cut);
466   gMC->SetCut("BCUTE",  cut); 
467   gMC->SetCut("BCUTM",  cut); 
468   gMC->SetCut("DCUTE",  cut); 
469   gMC->SetCut("DCUTM",  cut); 
470   gMC->SetCut("PPCUTM", cut);
471   gMC->SetCut("TOFMAX", tofmax); 
472
473
474
475
476   //======================//
477   // Set External decayer //
478   //======================//
479   TVirtualMCDecayer* decayer = new AliDecayerPythia();
480   decayer->SetForceDecay(kAll);
481   decayer->Init();
482   gMC->SetExternalDecayer(decayer);
483
484   //=========================//
485   // Generator Configuration //
486   //=========================//
487   AliGenerator* gener = 0x0;
488   
489   if (proc == kPythia6) {
490     gener = MbPythia();
491   } else if (proc == kPythia6D6T) {
492     gener = MbPythiaTuneD6T();
493   } else if (proc == kPythia6ATLAS) {
494     gener = MbPythiaTuneATLAS();
495   } else if (proc == kPythiaPerugia0) {
496     gener = MbPythiaTunePerugia0();
497   } else if (proc == kPythia6ATLAS_Flat) {
498     gener = MbPythiaTuneATLAS_Flat();
499   } else if (proc == kPhojet) {
500     gener = MbPhojet();
501   }
502   
503   
504   //
505   //
506   // Size of the interaction diamond
507   // Longitudinal
508   Float_t sigmaz  = 5.4 / TMath::Sqrt(2.); // [cm]
509   if (energy == 900)
510     //sigmaz  = 10.5 / TMath::Sqrt(2.); // [cm]
511     //sigmaz = 3.7;
512     if (energy == 7000)
513       sigmaz  = 6.3 / TMath::Sqrt(2.); // [cm]
514   
515   //
516   // Transverse
517
518   // beta*
519   Float_t betast                  = 10.0;  // beta* [m]
520   if (runNumber >= 117048) betast =  2.0;
521   if (runNumber >  122375) betast =  3.5;  // starting with fill 1179
522   //    
523   //
524   Float_t eps     = 5.0e-6;            // emittance [m]
525   Float_t gamma   = energy / 2.0 / 0.938272;  // relativistic gamma [1]
526   Float_t sigmaxy = TMath::Sqrt(eps * betast / gamma) / TMath::Sqrt(2.) * 100.;  // [cm]
527   printf("\n \n Diamond size x-y: %10.3e z: %10.3e\n \n", sigmaxy, sigmaz);
528     
529   gener->SetSigma(sigmaxy, sigmaxy, sigmaz);      // Sigma in (X,Y,Z) (cm) on IP position
530   gener->SetVertexSmear(kPerEvent);
531   gener->Init();
532
533   printf("\n \n Comment: %s \n \n", comment.Data());
534
535
536 }
537 //
538 //           PYTHIA
539 //
540
541 AliGenerator* MbPythia()
542 {
543   comment = comment.Append(" pp: Pythia low-pt");
544   //
545   //    Pythia
546   AliGenPythia* pythia = new AliGenPythia(-1); 
547   pythia->SetMomentumRange(0, 999999.);
548   pythia->SetThetaRange(0., 180.);
549   pythia->SetYRange(-12.,12.);
550   pythia->SetPtRange(0,1000.);
551   pythia->SetProcess(kPyMb);
552   pythia->SetEnergyCMS(energy);
553       
554   return pythia;
555 }
556
557 AliGenerator* MbPythiaTuneD6T()
558 {
559   comment = comment.Append(" pp: Pythia low-pt");
560   //
561   //    Pythia
562   AliGenPythia* pythia = new AliGenPythia(-1); 
563   pythia->SetMomentumRange(0, 999999.);
564   pythia->SetThetaRange(0., 180.);
565   pythia->SetYRange(-12.,12.);
566   pythia->SetPtRange(0,1000.);
567   pythia->SetProcess(kPyMb);
568   pythia->SetEnergyCMS(energy);
569   //    Tune
570   //    109     D6T : Rick Field's CDF Tune D6T (NB: needs CTEQ6L pdfs externally)
571   pythia->SetTune(109); // F I X 
572   pythia->SetStrucFunc(kCTEQ6l);
573   //
574   return pythia;
575 }
576
577 AliGenerator* MbPythiaTunePerugia0()
578 {
579   comment = comment.Append(" pp: Pythia low-pt (Perugia0)");
580   //
581   //    Pythia
582   AliGenPythia* pythia = new AliGenPythia(-1); 
583   pythia->SetMomentumRange(0, 999999.);
584   pythia->SetThetaRange(0., 180.);
585   pythia->SetYRange(-12.,12.);
586   pythia->SetPtRange(0,1000.);
587   pythia->SetProcess(kPyMb);
588   pythia->SetEnergyCMS(energy);
589   //    Tune
590   //    320     Perugia 0
591   pythia->SetTune(320); 
592   pythia->UseNewMultipleInteractionsScenario();
593   pythia->SetCrossingAngle(0,0.000280);
594
595   //
596   return pythia;
597 }
598
599
600 AliGenerator* MbPythiaTuneATLAS()
601 {
602   comment = comment.Append(" pp: Pythia low-pt");
603   //
604   //    Pythia
605   AliGenPythia* pythia = new AliGenPythia(-1); 
606   pythia->SetMomentumRange(0, 999999.);
607   pythia->SetThetaRange(0., 180.);
608   pythia->SetYRange(-12.,12.);
609   pythia->SetPtRange(0,1000.);
610   pythia->SetProcess(kPyMb);
611   pythia->SetEnergyCMS(energy);
612   //    Tune
613   //    C   306 ATLAS-CSC: Arthur Moraes' (new) ATLAS tune (needs CTEQ6L externally)
614   pythia->SetTune(306);
615   pythia->SetStrucFunc(kCTEQ6l);
616   //
617   return pythia;
618 }
619
620 AliGenerator* MbPythiaTuneATLAS_Flat()
621 {
622   AliGenPythia* pythia = MbPythiaTuneATLAS();
623       
624   comment = comment.Append("; flat multiplicity distribution");
625       
626   // set high multiplicity trigger
627   // this weight achieves a flat multiplicity distribution
628   Double_t cont[] =
629     {0, 
630      0.234836, 0.103484, 0.00984802, 0.0199906, 0.0260018, 0.0208481, 0.0101477, 0.00146998, -0.00681513, -0.0114928,
631      -0.0113352, -0.0102012, -0.00895238, -0.00534961, -0.00261648, -0.000819048, 0.00130816, 0.00177978, 0.00373838, 0.00566255,
632      0.00628156, 0.00687458, 0.00668017, 0.00702917, 0.00810848, 0.00876167, 0.00832783, 0.00848518, 0.0107709, 0.0106849,
633      0.00933702, 0.00912525, 0.0106553, 0.0102785, 0.0101756, 0.010962, 0.00957103, 0.00970448, 0.0117133, 0.012271,
634      0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113,
635      0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113,
636      0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113,
637      0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113,
638      0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113,
639      0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113, 0.0113,
640      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
641      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
642      0};
643
644   Double_t err[] =
645     {0, 
646      0.00168216, 0.000743548, 0.00103125, 0.00108605, 0.00117101, 0.00124577, 0.00129119, 0.00128341, 0.00121421, 0.00112431,
647      0.00100588, 0.000895078, 0.000790314, 0.000711673, 0.000634547, 0.000589133, 0.000542763, 0.000509552, 0.000487375, 0.000468906, 
648      0.000460196, 0.000455806, 0.00044843, 0.000449317, 0.00045007, 0.000458016, 0.000461167, 0.000474742, 0.00050161, 0.00051637, 
649      0.000542336, 0.000558854, 0.000599169, 0.000611982, 0.000663855, 0.000696322, 0.000722825, 0.000771686, 0.000838023, 0.000908317, 
650      0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003,
651      0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003,
652      0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003,
653      0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003,
654      0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003,
655      0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003,
656      0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
657      0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
658      0};
659
660   TH1F *weight = new TH1F("newweight","newweight",120,-0.5,119.5);
661
662   weight->SetContent(cont);
663   weight->SetError(err);
664         
665   Int_t limit = weight->GetRandom();
666   pythia->SetTriggerChargedMultiplicity(limit, 1.4);
667       
668   comment = comment.Append(Form("; multiplicity threshold set to %d in |eta| < 1.4", limit));
669
670   return pythia;
671 }
672
673 AliGenerator* MbPhojet()
674 {
675   comment = comment.Append(" pp: Pythia low-pt");
676   //
677   //    DPMJET
678 #if defined(__CINT__)
679   gSystem->Load("libdpmjet");      // Parton density functions
680   gSystem->Load("libTDPMjet");      // Parton density functions
681 #endif
682   AliGenDPMjet* dpmjet = new AliGenDPMjet(-1); 
683   dpmjet->SetMomentumRange(0, 999999.);
684   dpmjet->SetThetaRange(0., 180.);
685   dpmjet->SetYRange(-12.,12.);
686   dpmjet->SetPtRange(0,1000.);
687   dpmjet->SetProcess(kDpmMb);
688   dpmjet->SetEnergyCMS(energy);
689   dpmjet->SetCrossingAngle(0,0.000280);
690   return dpmjet;
691 }
692
693 void ProcessEnvironmentVars()
694 {
695   // Run type
696   if (gSystem->Getenv("CONFIG_RUN_TYPE")) {
697     for (Int_t iRun = 0; iRun < kRunMax; iRun++) {
698       if (strcmp(gSystem->Getenv("CONFIG_RUN_TYPE"), pprRunName[iRun])==0) {
699         proc = (PDC06Proc_t)iRun;
700         cout<<"Run type set to "<<pprRunName[iRun]<<endl;
701       }
702     }
703   }
704
705   // Field
706   if (gSystem->Getenv("CONFIG_FIELD")) {
707     for (Int_t iField = 0; iField < kFieldMax; iField++) {
708       if (strcmp(gSystem->Getenv("CONFIG_FIELD"), pprField[iField])==0) {
709         mag = (Mag_t)iField;
710         cout<<"Field set to "<<pprField[iField]<<endl;
711       }
712     }
713   }
714
715   // Energy
716   if (gSystem->Getenv("CONFIG_ENERGY")) {
717     energy = atoi(gSystem->Getenv("CONFIG_ENERGY"));
718     cout<<"Energy set to "<<energy<<" GeV"<<endl;
719   }
720
721   // Random Number seed
722   if (gSystem->Getenv("CONFIG_SEED")) {
723     seed = atoi(gSystem->Getenv("CONFIG_SEED"));
724   }
725
726   // Run number
727   if (gSystem->Getenv("DC_RUN")) {
728     runNumber = atoi(gSystem->Getenv("DC_RUN"));
729   }
730
731   // Physics lists
732   if (gSystem->Getenv("CONFIG_PHYSICSLIST")) {
733     for (Int_t iList = 0; iList < kListMax; iList++) {
734       if (strcmp(gSystem->Getenv("CONFIG_PHYSICSLIST"), physicsListName[iList])==0){
735         physicslist = (PhysicsList_t)iList;
736         cout<<"Physics list set to "<< physicsListName[iList]<<endl;
737       }
738     }
739   }
740
741 }