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