]> git.uio.no Git - u/mrichter/AliRoot.git/blob - test/vmctest/production/Config.C
Add process command which kills neutrinos in order to reduce simulation time
[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, kPythia6D6T_Flat, kPythia6ATLAS, kPythia6ATLAS_Flat, kPythia6_Perugia0, kPhojet, kRunMax
51 };
52
53 const char * pprRunName[] = {
54     "kPythia6", "kPythia6D6T", "kPythia6D6T_Flat", "kPythia6ATLAS", "kPythia6ATLAS_Flat", "kPythia6_Perugia0", "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_BERT_EMV_OPTICAL, CHIPS_OPTICAL, QGSP_BERT_CHIPS_OPTICAL, kListMax
69 };
70
71 const char * physicsListName[] = {
72   "QGSP_BERT_EMV",         "CHIPS",         "QGSP_BERT_CHIPS", 
73   "QGSP_BERT_EMV_OPTICAL", "CHIPS_OPTICAL", "QGSP_BERT_CHIPS_OPTICAL"
74 };
75
76 //--- Functions ---
77 class AliGenPythia;
78 AliGenerator *MbPythia();
79 AliGenerator *MbPythiaTuneD6T();
80 AliGenerator *MbPhojet();
81 void ProcessEnvironmentVars();
82
83 // Geterator, field, beam energy
84 static PDC06Proc_t   proc         = kPhojet;
85 static Mag_t         mag          = k5kG;
86 static Float_t       energy       = 10000; // energy in CMS
87 static PhysicsList_t physicslist  = QGSP_BERT_EMV;
88 //========================//
89 // Set Random Number seed //
90 //========================//
91 TDatime dt;
92 static UInt_t seed    = dt.Get();
93
94 // Comment line
95 static TString comment;
96
97 void Config()
98 {
99     
100
101   // Get settings from environment variables
102   ProcessEnvironmentVars();
103
104   gRandom->SetSeed(seed);
105   cerr<<"Seed for random number generation= "<<seed<<endl; 
106
107   // Libraries required by geant321
108 #if defined(__CINT__)
109   gSystem->Load("liblhapdf");       // Parton density functions
110   gSystem->Load("libEGPythia6");    // TGenerator interface
111
112   if (proc == kPythia6 || proc == kPhojet) {
113       gSystem->Load("libpythia6");        // Pythia 6.2 
114   } else {
115       gSystem->Load("libpythia6.4.21");   // Pythia 6.4
116   } 
117
118   gSystem->Load("libAliPythia6");   // ALICE specific implementations
119   //gSystem->Load("libgeant321");
120 #endif
121
122   //  new TGeant3TGeo("C++ Interface to Geant3");
123
124   //=======================================================================
125   //  Create the output file
126
127    
128   AliRunLoader* rl=0x0;
129
130   cout<<"Config.C: Creating Run Loader ..."<<endl;
131   rl = AliRunLoader::Open("galice.root",
132                           AliConfig::GetDefaultEventFolderName(),
133                           "recreate");
134   if (rl == 0x0)
135     {
136       gAlice->Fatal("Config.C","Can not instatiate the Run Loader");
137       return;
138     }
139   rl->SetCompressionLevel(2);
140   rl->SetNumberOfEventsPerFile(1000);
141   gAlice->SetRunLoader(rl);
142   // gAlice->SetGeometryFromFile("geometry.root");
143   // gAlice->SetGeometryFromCDB();
144   
145   // Set the trigger configuration: proton-proton
146   gAlice->SetTriggerDescriptor("p-p");
147   //  AliSimulation::Instance()->SetTriggerConfig("p-p");
148
149  
150   //
151   // FIELD
152   //
153   AliMagF* field = 0x0;
154   if (mag == kNoField) {
155     comment = comment.Append(" | L3 field 0.0 T");
156     field = new AliMagF("Maps","Maps", 0., 0., AliMagF::k5kGUniform,AliMagF::kBeamTypepp, energy/2.0);
157   } else if (mag == k5kG) {
158     comment = comment.Append(" | L3 field 0.5 T");
159     field = new AliMagF("Maps","Maps", -1., -1., AliMagF::k5kG, AliMagF::kBeamTypepp, energy/2.0);
160   }
161
162   printf("\n \n Comment: %s \n \n", comment.Data());
163
164   TGeoGlobalMagField::Instance()->SetField(field);
165     
166   rl->CdGAFile();
167   
168   Int_t iABSO  = 1;
169   Int_t iACORDE= 0;
170   Int_t iDIPO  = 1;
171   Int_t iEMCAL = 1;
172   Int_t iFMD   = 1;
173   Int_t iFRAME = 1;
174   Int_t iHALL  = 1;
175   Int_t iITS   = 1;
176   Int_t iMAG   = 1;
177   Int_t iMUON  = 1;
178   Int_t iPHOS  = 1;
179   Int_t iPIPE  = 1;
180   Int_t iPMD   = 1;
181   Int_t iHMPID = 1;
182   Int_t iSHIL  = 1;
183   Int_t iT0    = 1;
184   Int_t iTOF   = 1;
185   Int_t iTPC   = 1;
186   Int_t iTRD   = 1;
187   Int_t iVZERO = 1;
188   Int_t iZDC   = 1;
189   
190
191     //=================== Alice BODY parameters =============================
192     AliBODY *BODY = new AliBODY("BODY", "Alice envelop");
193
194
195     if (iMAG)
196     {
197         //=================== MAG parameters ============================
198         // --- Start with Magnet since detector layouts may be depending ---
199         // --- on the selected Magnet dimensions ---
200         AliMAG *MAG = new AliMAG("MAG", "Magnet");
201     }
202
203
204     if (iABSO)
205     {
206         //=================== ABSO parameters ============================
207         AliABSO *ABSO = new AliABSOv3("ABSO", "Muon Absorber");
208     }
209
210     if (iDIPO)
211     {
212         //=================== DIPO parameters ============================
213
214         AliDIPO *DIPO = new AliDIPOv3("DIPO", "Dipole version 3");
215     }
216
217     if (iHALL)
218     {
219         //=================== HALL parameters ============================
220
221         AliHALL *HALL = new AliHALLv3("HALL", "Alice Hall");
222     }
223
224
225     if (iFRAME)
226     {
227         //=================== FRAME parameters ============================
228
229         AliFRAMEv2 *FRAME = new AliFRAMEv2("FRAME", "Space Frame");
230         FRAME->SetHoles(1);
231     }
232
233     if (iSHIL)
234     {
235         //=================== SHIL parameters ============================
236
237         AliSHIL *SHIL = new AliSHILv3("SHIL", "Shielding Version 3");
238     }
239
240
241     if (iPIPE)
242     {
243         //=================== PIPE parameters ============================
244
245         AliPIPE *PIPE = new AliPIPEv3("PIPE", "Beam Pipe");
246     }
247  
248     if (iITS)
249     {
250         //=================== ITS parameters ============================
251
252         AliITS *ITS  = new AliITSv11Hybrid("ITS","ITS v11Hybrid");
253     }
254
255     if (iTPC)
256     {
257       //============================ TPC parameters =====================
258
259         AliTPC *TPC = new AliTPCv2("TPC", "Default");
260         TPC->SetPrimaryIonisation(); // not used with Geant3
261     }
262
263
264     if (iTOF) {
265         //=================== TOF parameters ============================
266
267         AliTOF *TOF = new AliTOFv6T0("TOF", "normal TOF");
268     }
269
270
271     if (iHMPID)
272     {
273         //=================== HMPID parameters ===========================
274
275         AliHMPID *HMPID = new AliHMPIDv3("HMPID", "normal HMPID");
276
277     }
278
279
280     if (iZDC)
281     {
282         //=================== ZDC parameters ============================
283
284         AliZDC *ZDC = new AliZDCv3("ZDC", "normal ZDC");
285     }
286
287     if (iTRD)
288     {
289         //=================== TRD parameters ============================
290
291         AliTRD *TRD = new AliTRDv1("TRD", "TRD slow simulator");
292         AliTRDgeometry *geoTRD = TRD->GetGeometry();
293         // Partial geometry: modules at 0,1,7,8,9,16,17
294         // starting at 3h in positive direction
295         geoTRD->SetSMstatus(2,0);
296         geoTRD->SetSMstatus(3,0);
297         geoTRD->SetSMstatus(4,0);
298         geoTRD->SetSMstatus(5,0);
299         geoTRD->SetSMstatus(6,0);
300         geoTRD->SetSMstatus(11,0);
301         geoTRD->SetSMstatus(12,0);
302         geoTRD->SetSMstatus(13,0);
303         geoTRD->SetSMstatus(14,0);
304         geoTRD->SetSMstatus(15,0);
305         geoTRD->SetSMstatus(16,0);
306     }
307
308     if (iFMD)
309     {
310         //=================== FMD parameters ============================
311
312         AliFMD *FMD = new AliFMDv1("FMD", "normal FMD");
313    }
314
315     if (iMUON)
316     {
317         //=================== MUON parameters ===========================
318         // New MUONv1 version (geometry defined via builders)
319
320         AliMUON *MUON = new AliMUONv1("MUON", "default");
321     }
322
323     if (iPHOS)
324     {
325         //=================== PHOS parameters ===========================
326
327         AliPHOS *PHOS = new AliPHOSv1("PHOS", "noCPV");
328         //Set simulation parameters different from the default ones.
329         AliPHOSSimParam* simEmc = AliPHOSSimParam::GetInstance() ;
330   
331         // APD noise of warm (+20C) PHOS:
332         // a2 = a1*(Y1/Y2)*(M1/M2), where a1 = 0.012 is APD noise at -25C,
333         // Y1 = 4.3 photo-electrons/MeV, Y2 = 1.7 p.e/MeV - light yields at -25C and +20C,
334         // M1 = 50, M2 = 50 - APD gain factors chosen for t1 = -25C and t2 = +20C,
335         // Y = MeanLightYield*APDEfficiency.
336
337         Float_t apdNoise = 0.012*2.5; 
338         simEmc->SetAPDNoise(apdNoise);
339
340         //Raw Light Yield at +20C
341         simEmc->SetMeanLightYield(18800);
342
343         //ADC channel width at +18C.
344         simEmc->SetADCchannelW(0.0125);
345     }
346
347
348     if (iPMD)
349     {
350         //=================== PMD parameters ============================
351
352         AliPMD *PMD = new AliPMDv1("PMD", "normal PMD");
353     }
354
355     if (iT0)
356     {
357         //=================== T0 parameters ============================
358         AliT0 *T0 = new AliT0v1("T0", "T0 Detector");
359     }
360
361     if (iEMCAL)
362     {
363         //=================== EMCAL parameters ============================
364
365         AliEMCAL *EMCAL = new AliEMCALv2("EMCAL", "EMCAL_FIRSTYEAR");
366     }
367
368      if (iACORDE)
369     {
370         //=================== ACORDE parameters ============================
371
372         AliACORDE *ACORDE = new AliACORDEv1("ACORDE", "normal ACORDE");
373     }
374
375      if (iVZERO)
376     {
377         //=================== ACORDE parameters ============================
378
379         AliVZERO *VZERO = new AliVZEROv7("VZERO", "normal VZERO");
380     }
381
382
383
384   // Load Geant4 + Geant4 VMC libraries
385   //
386   if (gClassTable->GetID("TGeant4") == -1) {
387     TString g4libsMacro = "$G4INSTALL/macro/g4libs.C";
388     //TString g4libsMacro = "$ALICE/geant4_vmc/examples/macro/g4libs.C";
389     // Load Geant4 libraries 
390     if (!gInterpreter->IsLoaded(g4libsMacro.Data())) {
391       gROOT->LoadMacro(g4libsMacro.Data());
392       gInterpreter->ProcessLine("g4libs()");
393     }
394   }    
395
396
397   // Create Geant4 VMC
398   //  
399   TGeant4 *geant4 = 0;
400   if ( ! gMC ) {
401     TG4RunConfiguration* runConfiguration=0x0;
402     for (Int_t iList = 0; iList < kListMax; iList++) {
403       if(iList<kListMax/2){
404         if(physicslist == iList){
405           runConfiguration = 
406             new TG4RunConfiguration("geomRoot", 
407                                     physicsListName[iList], 
408                                     "specialCuts+stackPopper+stepLimiter",
409                                     true);
410         }
411       }
412       else if(iList>=kListMax/2){//add "optical" PL to HadronPhysicsList
413         if(physicslist == iList){
414           runConfiguration = 
415             new TG4RunConfiguration("geomRoot", 
416                                     Form("%s+optical",physicsListName[iList-kListMax/2]), 
417                                     "specialCuts+stackPopper+stepLimiter",
418                                     true);
419         }
420       }
421     }
422     geant4 = new TGeant4("TGeant4", "The Geant4 Monte Carlo", runConfiguration);
423     cout << "Geant4 has been created." << endl;
424   } 
425   else {
426     cout << "Monte Carlo has been already created." << endl;
427   }  
428   
429   
430   
431   // Customization of Geant4 VMC
432   //
433   geant4->ProcessGeantCommand("/mcVerbose/all 1");  
434   geant4->ProcessGeantCommand("/mcVerbose/geometryManager 1");  
435   geant4->ProcessGeantCommand("/mcVerbose/opGeometryManager 1");  
436   geant4->ProcessGeantCommand("/mcTracking/loopVerbose 0");     
437   geant4->ProcessGeantCommand("/mcPhysics/rangeCuts 0.01 mm"); 
438   geant4->ProcessGeantCommand("/mcPhysics/selectOpProcess Scintillation");
439   geant4->ProcessGeantCommand("/mcPhysics/setOpProcessActivation false");
440   geant4->ProcessGeantCommand("/mcVerbose/composedPhysicsList 2");  
441   geant4->ProcessGeantCommand("/mcTracking/skipNeutrino true");
442
443
444
445
446
447  //
448   //=======================================================================
449   // ************* STEERING parameters FOR ALICE SIMULATION **************
450   // --- Specify event type to be tracked through the ALICE setup
451   // --- All positions are in cm, angles in degrees, and P and E in GeV
452
453
454     gMC->SetProcess("DCAY",1);
455     gMC->SetProcess("PAIR",1);
456     gMC->SetProcess("COMP",1);
457     gMC->SetProcess("PHOT",1);
458     gMC->SetProcess("PFIS",0);
459     gMC->SetProcess("DRAY",0);
460     gMC->SetProcess("ANNI",1);
461     gMC->SetProcess("BREM",1);
462     gMC->SetProcess("MUNU",1);
463     gMC->SetProcess("CKOV",1);
464     gMC->SetProcess("HADR",1);
465     gMC->SetProcess("LOSS",2);
466     gMC->SetProcess("MULS",1);
467     gMC->SetProcess("RAYL",1);
468
469     Float_t cut = 1.e-3;        // 1MeV cut by default
470     Float_t tofmax = 1.e10;
471
472     gMC->SetCut("CUTGAM", cut);
473     gMC->SetCut("CUTELE", cut);
474     gMC->SetCut("CUTNEU", cut);
475     gMC->SetCut("CUTHAD", cut);
476     gMC->SetCut("CUTMUO", cut);
477     gMC->SetCut("BCUTE",  cut); 
478     gMC->SetCut("BCUTM",  cut); 
479     gMC->SetCut("DCUTE",  cut); 
480     gMC->SetCut("DCUTM",  cut); 
481     gMC->SetCut("PPCUTM", cut);
482     gMC->SetCut("TOFMAX", tofmax); 
483
484
485
486
487   //======================//
488   // Set External decayer //
489   //======================//
490   TVirtualMCDecayer* decayer = new AliDecayerPythia();
491   decayer->SetForceDecay(kAll);
492   decayer->Init();
493   gMC->SetExternalDecayer(decayer);
494
495   //=========================//
496   // Generator Configuration //
497   //=========================//
498   AliGenerator* gener = 0x0;
499   
500   if (proc == kPythia6) {
501       gener = MbPythia();
502   } else if (proc == kPythia6D6T) {
503       gener = MbPythiaTuneD6T();
504   } else if (proc == kPythia6D6T_Flat) {
505       gener = MbPythiaTuneD6T_Flat();
506   } else if (proc == kPythia6ATLAS) {
507       gener = MbPythiaTuneATLAS();
508   } else if (proc == kPythia6ATLAS_Flat) {
509       gener = MbPythiaTuneATLAS_Flat();
510   } else if (proc == kPhojet) {
511       gener = MbPhojet();
512   } else if (proc == kPythia6_Perugia0) {
513       gener = MbPythiaTunePerugia0();
514   }
515
516   //
517   // PRIMARY VERTEX
518   //
519    Double_t xv = 0;
520    Double_t yv = 0;
521    Double_t zv = 0;
522
523    Double_t vtxPos[3];
524    Double_t vtxErr[3];
525    for (Int_t i = 0; i < 3; i++) {
526        vtxPos[i] = 0.;
527        vtxErr[i] = 0.;
528    }
529
530    if(AliCDBManager::Instance()->IsDefaultStorageSet()){
531       AliCDBEntry* entry = AliCDBManager::Instance()->Get("GRP/Calib/MeanVertexSPD");
532       AliESDVertex* vertex = dynamic_cast<AliESDVertex*> (entry->GetObject());
533       vertex->GetXYZ(vtxPos);
534       vertex->GetSigmaXYZ(vtxErr);
535     }
536
537     printf("Vertex position from OCDB entry: x = %13.3f, y = %13.3f, z = %13.3f (sigma = %13.3f)\n",
538     vtxPos[0], vtxPos[1], vtxPos[2], vtxErr[2]);
539
540      gener->SetOrigin(vtxPos[0], vtxPos[1], vtxPos[2]);   // vertex position
541   //
542   //
543   // Size of the interaction diamond
544   // Longitudinal
545      sigmaz = vtxErr[2];
546 //
547   // Transverse
548   Float_t betast  = 10;                 // beta* [m]
549   Float_t eps     = 2.5e-6;            // emittance [m]
550   Float_t gamma   = energy / 2.0 / 0.938272;  // relativistic gamma [1]
551   Float_t sigmaxy = TMath::Sqrt(eps * betast / gamma) / TMath::Sqrt(2.) * 100.;  // [cm]
552  
553   printf("\n \n Diamond size x-y: %10.3e z: %10.3e\n \n", sigmaxy, sigmaz);
554
555   gener->SetSigma(sigmaxy, sigmaxy, sigmaz);      // Sigma in (X,Y,Z) (cm) on IP position
556   gener->SetVertexSmear(kPerEvent);
557   gener->Init();
558
559
560
561
562
563
564
565
566
567
568
569
570
571 }
572 //
573 //           PYTHIA
574 //
575
576 AliGenerator* MbPythia()
577 {
578       comment = comment.Append(" pp: Pythia low-pt");
579 //
580 //    Pythia
581       AliGenPythia* pythia = new AliGenPythia(-1); 
582       pythia->SetMomentumRange(0, 999999.);
583       pythia->SetThetaRange(0., 180.);
584       pythia->SetYRange(-12.,12.);
585       pythia->SetPtRange(0,1000.);
586       pythia->SetProcess(kPyMb);
587       pythia->SetEnergyCMS(energy);
588       
589       return pythia;
590 }
591
592 AliGenerator* MbPythiaTuneD6T()
593 {
594       comment = comment.Append(" pp: Pythia low-pt");
595 //
596 //    Pythia
597       AliGenPythia* pythia = new AliGenPythia(-1); 
598       pythia->SetMomentumRange(0, 999999.);
599       pythia->SetThetaRange(0., 180.);
600       pythia->SetYRange(-12.,12.);
601       pythia->SetPtRange(0,1000.);
602       pythia->SetProcess(kPyMb);
603       pythia->SetEnergyCMS(energy);
604 //    Tune
605 //    109     D6T : Rick Field's CDF Tune D6T (NB: needs CTEQ6L pdfs externally)
606       pythia->SetTune(109); // F I X 
607       pythia->SetStrucFunc(kCTEQ6l);
608 //
609       return pythia;
610 }
611
612 AliGenerator* MbPythiaTunePerugia0()
613 {
614       comment = comment.Append(" pp: Pythia Minimum Bias Tune Perugia0");
615 //
616 //    Pythia
617       AliGenPythia* pythia = new AliGenPythia(-1); 
618       pythia->SetMomentumRange(0, 999999.);
619       pythia->SetThetaRange(0., 180.);
620       pythia->SetYRange(-12.,12.);
621       pythia->SetPtRange(0, 1000.);
622       pythia->SetProcess(kPyMb);
623       pythia->SetEnergyCMS(energy);
624 //    Tune
625       pythia->SetTune(320);
626 //
627       return pythia;
628 }
629
630 AliGenerator* MbPythiaTuneATLAS()
631 {
632       comment = comment.Append(" pp: Pythia low-pt");
633 //
634 //    Pythia
635       AliGenPythia* pythia = new AliGenPythia(-1); 
636       pythia->SetMomentumRange(0, 999999.);
637       pythia->SetThetaRange(0., 180.);
638       pythia->SetYRange(-12.,12.);
639       pythia->SetPtRange(0,1000.);
640       pythia->SetProcess(kPyMb);
641       pythia->SetEnergyCMS(energy);
642 //    Tune
643 //    C   306 ATLAS-CSC: Arthur Moraes' (new) ATLAS tune (needs CTEQ6L externally)
644       pythia->SetTune(306);
645       pythia->SetStrucFunc(kCTEQ6l);
646 //
647       return pythia;
648 }
649
650 AliGenerator* MbPythiaTuneATLAS_Flat()
651 {
652       AliGenPythia* pythia = MbPythiaTuneATLAS();
653       
654       comment = comment.Append("; flat multiplicity distribution");
655       
656       // set high multiplicity trigger
657       // this weight achieves a flat multiplicity distribution
658       TH1 *weight = new TH1D("weight","weight",201,-0.5,200.5);
659       weight->SetBinContent(1,5.49443);
660       weight->SetBinContent(2,8.770816);
661       weight->SetBinContent(6,0.4568624);
662       weight->SetBinContent(7,0.2919915);
663       weight->SetBinContent(8,0.6674189);
664       weight->SetBinContent(9,0.364737);
665       weight->SetBinContent(10,0.8818444);
666       weight->SetBinContent(11,0.531885);
667       weight->SetBinContent(12,1.035197);
668       weight->SetBinContent(13,0.9394057);
669       weight->SetBinContent(14,0.9643193);
670       weight->SetBinContent(15,0.94543);
671       weight->SetBinContent(16,0.9426507);
672       weight->SetBinContent(17,0.9423649);
673       weight->SetBinContent(18,0.789456);
674       weight->SetBinContent(19,1.149026);
675       weight->SetBinContent(20,1.100491);
676       weight->SetBinContent(21,0.6350525);
677       weight->SetBinContent(22,1.351941);
678       weight->SetBinContent(23,0.03233504);
679       weight->SetBinContent(24,0.9574557);
680       weight->SetBinContent(25,0.868133);
681       weight->SetBinContent(26,1.030998);
682       weight->SetBinContent(27,1.08897);
683       weight->SetBinContent(28,1.251382);
684       weight->SetBinContent(29,0.1391099);
685       weight->SetBinContent(30,1.192876);
686       weight->SetBinContent(31,0.448944);
687       weight->SetBinContent(32,1);
688       weight->SetBinContent(33,1);
689       weight->SetBinContent(34,1);
690       weight->SetBinContent(35,1);
691       weight->SetBinContent(36,0.9999997);
692       weight->SetBinContent(37,0.9999997);
693       weight->SetBinContent(38,0.9999996);
694       weight->SetBinContent(39,0.9999996);
695       weight->SetBinContent(40,0.9999995);
696       weight->SetBinContent(41,0.9999993);
697       weight->SetBinContent(42,1);
698       weight->SetBinContent(43,1);
699       weight->SetBinContent(44,1);
700       weight->SetBinContent(45,1);
701       weight->SetBinContent(46,1);
702       weight->SetBinContent(47,0.9999999);
703       weight->SetBinContent(48,0.9999998);
704       weight->SetBinContent(49,0.9999998);
705       weight->SetBinContent(50,0.9999999);
706       weight->SetBinContent(51,0.9999999);
707       weight->SetBinContent(52,0.9999999);
708       weight->SetBinContent(53,0.9999999);
709       weight->SetBinContent(54,0.9999998);
710       weight->SetBinContent(55,0.9999998);
711       weight->SetBinContent(56,0.9999998);
712       weight->SetBinContent(57,0.9999997);
713       weight->SetBinContent(58,0.9999996);
714       weight->SetBinContent(59,0.9999995);
715       weight->SetBinContent(60,1);
716       weight->SetBinContent(61,1);
717       weight->SetBinContent(62,1);
718       weight->SetBinContent(63,1);
719       weight->SetBinContent(64,1);
720       weight->SetBinContent(65,0.9999999);
721       weight->SetBinContent(66,0.9999998);
722       weight->SetBinContent(67,0.9999998);
723       weight->SetBinContent(68,0.9999999);
724       weight->SetBinContent(69,1);
725       weight->SetBinContent(70,1);
726       weight->SetBinContent(71,0.9999997);
727       weight->SetBinContent(72,0.9999995);
728       weight->SetBinContent(73,0.9999994);
729       weight->SetBinContent(74,1);
730       weight->SetBinContent(75,1);
731       weight->SetBinContent(76,1);
732       weight->SetBinContent(77,1);
733       weight->SetBinContent(78,0.9999999);
734       weight->SetBinContent(79,1);
735       weight->SetBinContent(80,1);
736       weight->SetEntries(526);
737         
738       Int_t limit = weight->GetRandom();
739       pythia->SetTriggerChargedMultiplicity(limit, 1.4);
740       
741       comment = comment.Append(Form("; multiplicity threshold set to %d in |eta| < 1.4", limit));
742
743       return pythia;
744 }
745 AliGenerator* MbPythiaTuneD6T_Flat()
746 {
747       AliGenPythia* pythia = MbPythiaTuneD6T();
748       
749       comment = comment.Append("; flat multiplicity distribution");
750       
751       // set high multiplicity trigger
752       // this weight achieves a flat multiplicity distribution
753       TH1 *weight = new TH1D("weight","weight",201,-0.5,200.5);
754       weight->SetBinContent(1,5.49443);
755       weight->SetBinContent(2,8.770816);
756       weight->SetBinContent(6,0.4568624);
757       weight->SetBinContent(7,0.2919915);
758       weight->SetBinContent(8,0.6674189);
759       weight->SetBinContent(9,0.364737);
760       weight->SetBinContent(10,0.8818444);
761       weight->SetBinContent(11,0.531885);
762       weight->SetBinContent(12,1.035197);
763       weight->SetBinContent(13,0.9394057);
764       weight->SetBinContent(14,0.9643193);
765       weight->SetBinContent(15,0.94543);
766       weight->SetBinContent(16,0.9426507);
767       weight->SetBinContent(17,0.9423649);
768       weight->SetBinContent(18,0.789456);
769       weight->SetBinContent(19,1.149026);
770       weight->SetBinContent(20,1.100491);
771       weight->SetBinContent(21,0.6350525);
772       weight->SetBinContent(22,1.351941);
773       weight->SetBinContent(23,0.03233504);
774       weight->SetBinContent(24,0.9574557);
775       weight->SetBinContent(25,0.868133);
776       weight->SetBinContent(26,1.030998);
777       weight->SetBinContent(27,1.08897);
778       weight->SetBinContent(28,1.251382);
779       weight->SetBinContent(29,0.1391099);
780       weight->SetBinContent(30,1.192876);
781       weight->SetBinContent(31,0.448944);
782       weight->SetBinContent(32,1);
783       weight->SetBinContent(33,1);
784       weight->SetBinContent(34,1);
785       weight->SetBinContent(35,1);
786       weight->SetBinContent(36,0.9999997);
787       weight->SetBinContent(37,0.9999997);
788       weight->SetBinContent(38,0.9999996);
789       weight->SetBinContent(39,0.9999996);
790       weight->SetBinContent(40,0.9999995);
791       weight->SetBinContent(41,0.9999993);
792       weight->SetBinContent(42,1);
793       weight->SetBinContent(43,1);
794       weight->SetBinContent(44,1);
795       weight->SetBinContent(45,1);
796       weight->SetBinContent(46,1);
797       weight->SetBinContent(47,0.9999999);
798       weight->SetBinContent(48,0.9999998);
799       weight->SetBinContent(49,0.9999998);
800       weight->SetBinContent(50,0.9999999);
801       weight->SetBinContent(51,0.9999999);
802       weight->SetBinContent(52,0.9999999);
803       weight->SetBinContent(53,0.9999999);
804       weight->SetBinContent(54,0.9999998);
805       weight->SetBinContent(55,0.9999998);
806       weight->SetBinContent(56,0.9999998);
807       weight->SetBinContent(57,0.9999997);
808       weight->SetBinContent(58,0.9999996);
809       weight->SetBinContent(59,0.9999995);
810       weight->SetBinContent(60,1);
811       weight->SetBinContent(61,1);
812       weight->SetBinContent(62,1);
813       weight->SetBinContent(63,1);
814       weight->SetBinContent(64,1);
815       weight->SetBinContent(65,0.9999999);
816       weight->SetBinContent(66,0.9999998);
817       weight->SetBinContent(67,0.9999998);
818       weight->SetBinContent(68,0.9999999);
819       weight->SetBinContent(69,1);
820       weight->SetBinContent(70,1);
821       weight->SetBinContent(71,0.9999997);
822       weight->SetBinContent(72,0.9999995);
823       weight->SetBinContent(73,0.9999994);
824       weight->SetBinContent(74,1);
825       weight->SetBinContent(75,1);
826       weight->SetBinContent(76,1);
827       weight->SetBinContent(77,1);
828       weight->SetBinContent(78,0.9999999);
829       weight->SetBinContent(79,1);
830       weight->SetBinContent(80,1);
831       weight->SetEntries(526);
832         
833       Int_t limit = weight->GetRandom();
834       pythia->SetTriggerChargedMultiplicity(limit, 1.4);
835       
836       comment = comment.Append(Form("; multiplicity threshold set to %d in |eta| < 1.4", limit));
837
838       return pythia;
839 }
840
841 AliGenerator* MbPhojet()
842 {
843       comment = comment.Append(" pp: Pythia low-pt");
844 //
845 //    DPMJET
846 #if defined(__CINT__)
847   gSystem->Load("libdpmjet");      // Parton density functions
848   gSystem->Load("libTDPMjet");      // Parton density functions
849 #endif
850       AliGenDPMjet* dpmjet = new AliGenDPMjet(-1); 
851       dpmjet->SetMomentumRange(0, 999999.);
852       dpmjet->SetThetaRange(0., 180.);
853       dpmjet->SetYRange(-12.,12.);
854       dpmjet->SetPtRange(0,1000.);
855       dpmjet->SetProcess(kDpmMb);
856       dpmjet->SetEnergyCMS(energy);
857
858       return dpmjet;
859 }
860
861 void ProcessEnvironmentVars()
862 {
863     // Run type
864     if (gSystem->Getenv("CONFIG_RUN_TYPE")) {
865       for (Int_t iRun = 0; iRun < kRunMax; iRun++) {
866         if (strcmp(gSystem->Getenv("CONFIG_RUN_TYPE"), pprRunName[iRun])==0) {
867           proc = (PDC06Proc_t)iRun;
868           cout<<"Run type set to "<<pprRunName[iRun]<<endl;
869         }
870       }
871     }
872
873     // Field
874     if (gSystem->Getenv("CONFIG_FIELD")) {
875       for (Int_t iField = 0; iField < kFieldMax; iField++) {
876         if (strcmp(gSystem->Getenv("CONFIG_FIELD"), pprField[iField])==0) {
877           mag = (Mag_t)iField;
878           cout<<"Field set to "<<pprField[iField]<<endl;
879         }
880       }
881     }
882
883     // Energy
884     if (gSystem->Getenv("CONFIG_ENERGY")) {
885       energy = atoi(gSystem->Getenv("CONFIG_ENERGY"));
886       cout<<"Energy set to "<<energy<<" GeV"<<endl;
887     }
888
889     // Random Number seed
890     if (gSystem->Getenv("CONFIG_SEED")) {
891       seed = atoi(gSystem->Getenv("CONFIG_SEED"));
892     }
893
894     // Physics lists
895     if (gSystem->Getenv("CONFIG_PHYSICSLIST")) {
896       for (Int_t iList = 0; iList < kListMax; iList++) {
897         if (strcmp(gSystem->Getenv("CONFIG_PHYSICSLIST"), physicsListName[iList])==0) {
898           physicslist = (PhysicsList_t)iList;
899           cout<<"Physics list set to "<< physicsListName[iList]<<endl;
900         }
901       }
902     }
903
904 }