macro to run generation only with choice of generators (S.Jena)
authorshahoian <ruben.shahoyan@cern.ch>
Fri, 21 Nov 2014 16:28:53 +0000 (17:28 +0100)
committershahoian <ruben.shahoyan@cern.ch>
Fri, 21 Nov 2014 17:36:49 +0000 (18:36 +0100)
test/genkine/gen/fastMcProduction.C [new file with mode: 0644]

diff --git a/test/genkine/gen/fastMcProduction.C b/test/genkine/gen/fastMcProduction.C
new file mode 100644 (file)
index 0000000..fc40955
--- /dev/null
@@ -0,0 +1,987 @@
+/***************************************************************************
+ *  Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                         *
+ *  Author: ALICE OFFLICE.                                                 *
+ *  Contributors are mentioned in the code where appropriate.              *
+ *                                                                         *
+ *  Permission to use, copy, modify and distribute this software and its   *
+ *  documentation strictly for non-commercial purposes is hereby granted   *
+ *  without fee, provided that the above copyright notice appears in all   *
+ *  copies and that both the copyright notice and this permission notice   *
+ *  appear in the supporting documentation. The authors make no claims     *
+ *  about the suitability of this software for any purpose. It is          *
+ *  provided "as is" without express or implied warranty.                  *
+ *                                                                         *
+ ***************************************************************************
+     $Satyajit Jena || alien:sjena Sun Apr 21 14:05:19 CEST 2013$
+                                                                          
+       Sterring macro for fast production of MC events - Followed from 
+       the original macro of fastGenAmpt.C     
+     
+     Implemented Generators: (Version 1.0: Sun Apr 21 14:05:19 CEST 2013)
+     ----------------------------------------------------------------
+      kPythia6,            kPythia8,               kPythia6D6T,       
+      kPythiaPerugia0,     kPythia6ATLAS,          
+      kPythiaJets,         
+      kPhojet,                
+      kDPMjet,             kDPMjet_pA, 
+      kAmptDefault,        kAmptStringMelting,          kAmptStringMeltingNoART,
+      kAmptpA,                kAmptFlow,
+      kAmptReducedFlow,
+
+     FIXME: 
+     kPythia6ATLAS_Flat, 
+     kHijing,             kHijing2000,            kHijing_pA,             
+
+                                                                
+ ***************************************************************************/
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <Riostream.h>
+#include <TStopwatch.h>
+#include <TDatime.h>
+#include <TRandom.h>
+#include <TFile.h>
+#include <TTree.h>
+#include <TMath.h>
+#include <TParticle.h>
+#include <TPDGCode.h>
+#include <TDatabasePDG.h>
+#include <TRandom3.h>
+#include <TChain.h>
+#include "AliRun.h"
+#include "AliRunLoader.h"
+#include "AliHeader.h"
+#include "AliStack.h"
+#include "AliPDG.h"
+#include "AliGenAmpt.h"
+#include "TAmpt.h"
+#include "PYTHIA6/AliDecayerPythia.h"
+#include <TVirtualMC.h>
+#include <TF1.h>
+#include "STEER/AliConfig.h"
+#include "STEER/AliGenerator.h"
+#include "STEER/AliLog.h"
+#include "EVGEN/AliGenHIJINGpara.h"
+#include "THijing/AliGenHijing.h"
+#include "EVGEN/AliGenCocktail.h"
+#include "EVGEN/AliGenSlowNucleons.h"
+#include "EVGEN/AliSlowNucleonModelExp.h"
+#include "EVGEN/AliGenParam.h"
+#include "EVGEN/AliGenMUONlib.h"
+#include "EVGEN/AliGenSTRANGElib.h"
+#include "EVGEN/AliGenMUONCocktail.h"
+#include "EVGEN/AliGenCocktail.h"
+#include "EVGEN/AliGenGeVSim.h"
+#include "EVGEN/AliGeVSimParticle.h"
+#include "PYTHIA6/AliGenPythia.h"
+#endif
+
+//___________________________________________________________________________
+enum PDC06Proc_t { 
+  kPythia6,
+  kPythia8, 
+  kPythia6D6T, 
+  kPythiaPerugia0, 
+  kPythia6ATLAS, 
+  kPythia6ATLAS_Flat, 
+  kPythiaJets,
+  kPhojet, 
+  kHijing, 
+  kHijing2000, 
+  kHijing_pA,
+  kDPMjet, 
+  kDPMjet_pA,
+  kAmptDefault,
+  kAmpt, 
+  kAmptpA,
+  kAmptFlowStringMelting,
+  kAmptStringMeltingNoART,
+  kAmptFlow,
+  kAmptReducedFlow,
+  kRunMax
+};
+
+//___________________________________________________________________________
+const char * pprRunName[] = {
+  "kPythia6", 
+  "kPythia8",
+  "kPythia6D6T", 
+  "kPythiaPerugia0", 
+  "kPythia6ATLAS", 
+  "kPythia6ATLAS_Flat", 
+  "kPythiaJets", 
+  "kPhojet",  
+  "kHijing", 
+  "kHijing2000", 
+  "kHijing_pA",
+  "kDPMjet", 
+  "kDPMjet_pA", 
+  "kAmpt",
+  "kAmptpA", 
+  "kAmptFlow",
+  "kAmptReducedFlow"  
+};
+
+enum PprTrigConf_t {kDefaultPPTrig, kDefaultPbPbTrig };
+const char * pprTrigConfName[] = {"p-p","Pb-Pb"};
+
+//___________________________________________________________________________
+void ProcessEnvironmentVars();
+class AliGenPythia;
+AliGenerator *MbPythia();
+AliGenerator *Pythia8();
+AliGenerator *MbPythiaTuneD6T();
+AliGenerator *MbPythiaTunePerugia0();
+AliGenerator *MbPythiaTuneATLAS();
+AliGenerator *MbPythiaTuneATLAS_Flat();
+AliGenerator *PythiaJets();
+AliGenerator *MbPhojet();
+AliGenerator *Hijing();
+AliGenerator *Hijing2000();
+AliGenerator *Hijing_pA(Bool_t kSlowN);
+AliGenerator *DPMjet();
+AliGenerator *DPMjet_pA(Bool_t fragments);
+AliGenerator *Ampt();
+AliGenerator *AmptpA();
+AliGenerator* AmptFlow();
+AliGenerator *AmptReducedFlow();
+
+
+//_________________________________________________________________________
+//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+// Geterator, field, beam energy
+
+//static Double_t    pBeamEnergy = 4000.0;  // Used during pA runs
+//static Double_t  energy  = 2.*pBeamEnergy*TMath::Sqrt(82./208.); //energy in CMS 
+
+static Double_t      energy    = 7000.;
+static PDC06Proc_t   proc      = kPythia6;
+static Float_t       bMin      = 0.;
+static Float_t       bMax      = 100.;
+static PprTrigConf_t strig     = kDefaultPPTrig; // default pp trigger configuration
+
+
+static Double_t  JpsiPol      = 0; // Jpsi polarisation
+static Bool_t    JpsiHarderPt = kFALSE; // Jpsi harder pt spectrum (8.8 TeV)
+static TString comment;
+//static PprTrigConf_t strig    = kDefaultPbPbTrig; // default pp trigger configuration
+TDatime dt; 
+static UInt_t seed    = dt.Get();
+
+//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+//_________________________________________________________________________
+void fastMcProduction(Int_t nev = 300) {
+
+  ProcessEnvironmentVars();
+
+  gRandom->SetSeed(seed);
+  cerr<<" Seed for random number generation = "<<seed<<endl; 
+  
+#if defined(__CINT__)
+  gSystem->Load("liblhapdf");  
+  gSystem->Load("libEGPythia6"); 
+  
+  if (proc == kPythia6 || proc == kPhojet || proc == kDPMjet || proc==kDPMjet_pA) {
+    gSystem->Load("libpythia6");        // Pythia 6.2
+    gSystem->Load("libAliPythia6");     // ALICE specific implementations
+  }
+    
+  if (proc == kHijing || proc == kHijing2000 || proc == kHijing_pA ) {
+    gSystem->Load("libhijing");        
+    gSystem->Load("libTHijing");
+  } 
+  
+  else if ( proc == kDPMjet || proc== kDPMjet_pA ) {
+    gSystem->Load("libdpmjet"); 
+    gSystem->Load("libTDPMjet");
+  } 
+  
+  else if (proc == kAmptDefault || kAmptFlowStringMelting || proc ==  kAmptStringMeltingNoART || proc == kAmptpA || proc == kAmptReducedFlow) {
+    gSystem->Load("libampt");  
+    gSystem->Load("libTAmpt");
+    gSystem->Load("libpythia6.so");
+    gSystem->Load("libAliPythia6.so");
+  } 
+
+  if (proc == kPythia8) {
+    gSystem->Load("libpythia8.so");
+    gSystem->Load("libAliPythia8.so");
+    gSystem->Setenv("PYTHIA8DATA", gSystem->ExpandPathName("$ALICE_ROOT/PYTHIA8/pythia8145/xmldoc"));
+    gSystem->Setenv("LHAPDF",      gSystem->ExpandPathName("$ALICE_ROOT/LHAPDF"));
+    gSystem->Setenv("LHAPATH",     gSystem->ExpandPathName("$ALICE_ROOT/LHAPDF/PDFsets"));
+  }
+#endif
+
+
+ AliGenerator* gener = 0x0;
+ cout<<"Run type set to ------------- "<<pprRunName[proc]<<"   " << proc << "    " << kDPMjet_pA<< endl;
+
+ if (proc == kPythia6) {
+   Printf("<<<<<<<<<<<<<<<<<<< Processing Mb. PYTHIA >>>>>>>>>>>>>>>>>>>>"); 
+   gener = MbPythia();
+ } 
+ else if (proc == kPythia8) {
+   Printf("<<<<<<<<<<<<<<<<<<< Processing Mb. Pythia 8 >>>>>>>>>>>>>>>>>>>>"); 
+   gener = Pythia8();
+ }  
+ else if (proc == kPythia6D6T) {
+   Printf("<<<<<<<<<<<<<<<<<<< Processing Mb. PYTHIA D6T >>>>>>>>>>>>>>>>>>>>"); 
+   gener = MbPythiaTuneD6T();
+ } 
+
+ else if (proc == kPythiaPerugia0) {
+   Printf("<<<<<<<<<<<<<<<<<<< Processing Mb. PYTHIA Perugia0 >>>>>>>>>>>>>>>>>>>>"); 
+   gener = MbPythiaTunePerugia0();
+ } 
+ else if (proc == kPythia6ATLAS) {
+   Printf("<<<<<<<<<<<<<<<<<<< Processing Mb. PYTHIA ATLAS>>>>>>>>>>>>>>>>>>>>"); 
+   gener = MbPythiaTuneATLAS();
+ } 
+ else if (proc == kPythia6ATLAS_Flat) {
+   Printf("<<<<<<<<<<<<<<<<<<< Processing Mb. PYTHIA ATLAS_FLAT >>>>>>>>>>>>>>>>>>>>"); 
+   gener = MbPythiaTuneATLAS_Flat();
+ } 
+ else if (proc == kPythiaJets ) {
+   Printf("<<<<<<<<<<<<<<<<<<< Processing Pythia Jets >>>>>>>>>>>>>>>>>>>>"); 
+   gener = PythiaJets();
+ } 
+ else if (proc == kPhojet) {
+   Printf("<<<<<<<<<<<<<<<<<<< Processing Mb. PHOJET >>>>>>>>>>>>>>>>>>>>"); 
+   gener = MbPhojet();
+ } 
+
+ else if (proc == kHijing) {
+   Printf("<<<<<<<<<<<<<<<<<<< Processing Mb. HIJING >>>>>>>>>>>>>>>>>>>>"); 
+   gener = Hijing();   
+ } 
+ else if (proc == kHijing2000) {
+   Printf("<<<<<<<<<<<<<<<<<<< Processing Mb. HIJING 2000 >>>>>>>>>>>>>>>>>>>>"); 
+   gener = Hijing2000();       
+ }
+ else if (proc ==kHijing_pA) {
+   Printf("<<<<<<<<<<<<<<<<<<< Processing Mb. pA Hijing >>>>>>>>>>>>>>>>>>>>"); 
+   gener = Hijing_pA(kTRUE);
+ } 
+ else if (proc == kDPMjet) {
+   Printf("<<<<<<<<<<<<<<<<<<< Processing  DMPJet  >>>>>>>>>>>>>>>>>>>>");
+   gener = DPMjet();   
+ } 
+
+ else if (proc == kDPMjet_pA) {
+   Printf("<<<<<<<<<<<<<<<<<<< Processing  DMPJet pA >>>>>>>>>>>>>>>>>>>>");
+   gener = DPMjet_pA(kFALSE);  
+ } 
+ else if (proc == kAmptDefault) {
+   Printf("<<<<<<<<<<<<<<<<<<< Processing  AMPT Default >>>>>>>>>>>>>>>>>>>>");
+   gener = AmptDefault();
+ } 
+
+ else if (proc == kAmptStringMelting) {
+     Printf("<<<<<<<<<<<<<<<<<<< Processing  AMPT With Flow  >>>>>>>>>>>>>>>>>>>>");
+     gener = AmptStringMelting();
+ }
+
+ else if (proc == kAmptStringMeltingNoART) {
+     Printf("<<<<<<<<<<<<<<<<<<< Processing  AMPT With Flow  >>>>>>>>>>>>>>>>>>>>");
+     gener = AmptStringMeltingNoART();
+ }
+
+ else if (proc == kAmptpA) {
+   Printf("<<<<<<<<<<<<<<<<<<< Processing  AMPT pA  >>>>>>>>>>>>>>>>>>>>");
+   gener = AmptpA();
+ } 
+
+ else if (proc == kAmptReducedFlow) {
+   // Specific Fastgen
+   Printf("<<<<<<<<<<<<<<<<<<< Processing  AMPT With Reduced Flow >>>>>>>>>>>>>>>>>>>>");
+   gener = AmptReducedFlow();
+ } 
+
+ else {
+   cout << "ERROR : Wrong Procss Selcted !!!" << endl;
+   return;
+ }
+
+
+ AliPDG::AddParticlesToPdgDataBase();
+ TDatabasePDG::Instance();
+ const char* filename = "galice.root";
+ AliRunLoader* rl = AliRunLoader::Open(filename,"FASTRUN","recreate");
+ rl->SetCompressionLevel(2);
+ rl->SetNumberOfEventsPerFile(nev);
+ rl->LoadKinematics("RECREATE");
+ rl->MakeTree("E");
+ gAlice->SetRunLoader(rl);
+ rl->MakeStack();
+ AliStack* stack = rl->Stack();
+ AliHeader* header = rl->GetHeader();
+ /*
+   Float_t sigmaz  = 5.4 / TMath::Sqrt(2.); // [cm]
+   Float_t betast  = 3.5;                      // beta* [m]
+   Float_t eps     = 3.75e-6;                   // emittance [m]
+   Float_t gamma   = energy / 2.0 / 0.938272;  // relativistic gamma [1]
+   Float_t sigmaxy = TMath::Sqrt(eps * betast / gamma) / TMath::Sqrt(2.) * 100.;  // [cm]
+   printf("\n \n Diamond size x-y: %10.3e z: %10.3e\n \n", sigmaxy, sigmaz);
+   gener->SetSigma(sigmaxy, sigmaxy, sigmaz);      // Sigma in (X,Y,Z) (cm) on IP position
+   gener->SetVertexSmear(kPerEvent);
+ */
+
+ gener->Init();
+ gener->SetStack(stack);
+ rl->CdGAFile();
+ TStopwatch timer;
+ timer.Start();
+ for (Int_t iev = 0; iev < nev; iev++) {
+   cout <<"============================================= Event number "<< iev << endl;
+   //  Initialize event
+   header->Reset(0,iev);
+   rl->SetEventNumber(iev);
+   stack->Reset();
+   rl->MakeTree("K");
+   Int_t nprim = 0;
+   Int_t ntrial = 0;
+   //  Int_t ndstar = 0;
+   stack->Reset();
+   stack->ConnectTree(rl->TreeK());
+   gener->Generate();
+   ntrial++;
+   nprim = stack->GetNprimary();
+   cout << "Number of particles " << nprim << endl;
+   cout << "Number of trials " << ntrial << endl;
+   header->SetNprimary(stack->GetNprimary());
+   header->SetNtrack(stack->GetNtrack());  
+   stack->FinishEvent();
+    header->SetStack(stack);
+    rl->TreeE()->Fill();
+    rl->WriteKinematics("OVERWRITE");
+ } // event loop
+  timer.Stop();
+  timer.Print();
+  gener->FinishRun();
+  rl->WriteHeader("OVERWRITE");
+  gener->Write();
+  rl->Write();
+}
+
+//___________________________________________________//
+void ProcessEnvironmentVars() {
+    // Run type
+    if (gSystem->Getenv("CONFIG_RUN_TYPE")) {
+      for (Int_t iRun = 0; iRun < kRunMax; iRun++) {
+        if (strcmp(gSystem->Getenv("CONFIG_RUN_TYPE"), pprRunName[iRun]) == 0) {
+          proc = (PDC06Proc_t)iRun;
+          cout<<"Run type set to "<<pprRunName[iRun]<<endl;
+        }
+      }
+    }
+
+    
+    // Energy
+    if (gSystem->Getenv("CONFIG_ENERGY")) {
+      energy = atoi(gSystem->Getenv("CONFIG_ENERGY"));
+      cout<<"Energy set to "<<energy<<" GeV"<<endl;
+    }
+
+    // Random Number seed
+    if (gSystem->Getenv("CONFIG_SEED")) {
+      seed = atoi(gSystem->Getenv("CONFIG_SEED"));
+    }
+
+    // Impact param
+    if (gSystem->Getenv("CONFIG_BMIN")) {
+      bMin = atof(gSystem->Getenv("CONFIG_BMIN"));
+    }
+
+    if (gSystem->Getenv("CONFIG_BMAX")) {
+      bMax = atof(gSystem->Getenv("CONFIG_BMAX"));
+    }
+    cout<<"Impact parameter in ["<<bMin<<","<<bMax<<"]"<<endl;
+}
+
+
+
+//______________________________________________________________________
+AliGenerator* MbPythia() // Mb Pythia
+{
+      comment = comment.Append(" pp: Pythia low-pt");
+      AliGenPythia* pythia = new AliGenPythia(-1); 
+      /* pythia->SetMomentumRange(0, 999999.);
+      pythia->SetThetaRange(0., 180.);
+      pythia->SetYRange(-12.,12.);
+      pythia->SetPtRange(0,1000.);*/
+      pythia->SetProcess(kPyMb);
+      pythia->SetEnergyCMS(energy);
+      
+      return pythia;
+}
+
+
+//______________________________________________________________________
+AliGenerator* Pythia8()
+{
+  AliGenPythiaPlus* gener = new AliGenPythiaPlus(AliPythia8::Instance());
+  gener->SetProcess(kPyMbDefault);
+  gener->SetEnergyCMS(energy);
+  gener->SetEventListRange(-1, 2);
+  return gener;
+}
+
+
+
+//______________________________________________________________________
+AliGenerator* MbPythiaTuneD6T()
+{
+      comment = comment.Append(" pp: Pythia low-pt");
+      AliGenPythia* pythia = new AliGenPythia(-1); 
+      pythia->SetMomentumRange(0, 999999.);
+      pythia->SetThetaRange(0., 180.);
+      pythia->SetYRange(-12.,12.);
+      pythia->SetPtRange(0,1000.);
+      pythia->SetProcess(kPyMb);
+      pythia->SetEnergyCMS(energy);
+//    Tune
+//    109     D6T : Rick Field's CDF Tune D6T (NB: needs CTEQ6L pdfs externally)
+      pythia->SetTune(109); // F I X 
+      pythia->SetStrucFunc(kCTEQ6l);
+//
+      return pythia;
+}
+
+//______________________________________________________________________
+AliGenerator* MbPythiaTunePerugia0()
+{
+      comment = comment.Append(" pp: Pythia low-pt (Perugia0)");
+//
+//    Pythia
+      AliGenPythia* pythia = new AliGenPythia(-1); 
+      /* pythia->SetMomentumRange(0, 999999.);
+      pythia->SetThetaRange(0., 180.);
+      pythia->SetYRange(-12.,12.);
+      pythia->SetPtRange(0,1000.);*/
+      pythia->SetProcess(kPyMb);
+      pythia->SetEnergyCMS(energy);
+//    Tune
+//    320     Perugia 0
+      pythia->SetTune(320); 
+      pythia->UseNewMultipleInteractionsScenario();
+//
+      return pythia;
+}
+
+//______________________________________________________________________
+AliGenerator* MbPythiaTuneATLAS()
+{
+      comment = comment.Append(" pp: Pythia low-pt");
+//
+//    Pythia
+      AliGenPythia* pythia = new AliGenPythia(-1); 
+      /*   pythia->SetMomentumRange(0, 999999.);
+      pythia->SetThetaRange(0., 180.);
+      pythia->SetYRange(-12.,12.);
+      pythia->SetPtRange(0,1000.);*/
+      pythia->SetProcess(kPyMb);
+      pythia->SetEnergyCMS(energy);
+//    Tune
+//    C   306 ATLAS-CSC: Arthur Moraes' (new) ATLAS tune (needs CTEQ6L externally)
+      pythia->SetTune(306);
+      pythia->SetStrucFunc(kCTEQ6l);
+//
+      return pythia;
+}
+
+
+//______________________________________________________________________
+AliGenerator* MbPythiaTuneATLAS_Flat()
+{
+  AliGenPythia* pythia = MbPythiaTuneATLAS();
+  
+  comment = comment.Append("; flat multiplicity distribution");
+  
+  // set high multiplicity trigger
+  // this weight achieves a flat multiplicity distribution
+  TH1 *weight = new TH1D("weight","weight",201,-0.5,200.5);
+  weight->SetBinContent(1,5.49443);
+  weight->SetBinContent(2,8.770816);
+  weight->SetBinContent(6,0.4568624);
+  weight->SetBinContent(7,0.2919915);
+  weight->SetBinContent(8,0.6674189);
+  weight->SetBinContent(9,0.364737);
+  weight->SetBinContent(10,0.8818444);
+  weight->SetBinContent(11,0.531885);
+  weight->SetBinContent(12,1.035197);
+  weight->SetBinContent(13,0.9394057);
+  weight->SetBinContent(14,0.9643193);
+  weight->SetBinContent(15,0.94543);
+  weight->SetBinContent(16,0.9426507);
+  weight->SetBinContent(17,0.9423649);
+  weight->SetBinContent(18,0.789456);
+  weight->SetBinContent(19,1.149026);
+  weight->SetBinContent(20,1.100491);
+  weight->SetBinContent(21,0.6350525);
+  weight->SetBinContent(22,1.351941);
+  weight->SetBinContent(23,0.03233504);
+  weight->SetBinContent(24,0.9574557);
+  weight->SetBinContent(25,0.868133);
+  weight->SetBinContent(26,1.030998);
+  weight->SetBinContent(27,1.08897);
+  weight->SetBinContent(28,1.251382);
+  weight->SetBinContent(29,0.1391099);
+  weight->SetBinContent(30,1.192876);
+  weight->SetBinContent(31,0.448944);
+  weight->SetBinContent(32,1);
+  weight->SetBinContent(33,1);
+  weight->SetBinContent(34,1);
+  weight->SetBinContent(35,1);
+  weight->SetBinContent(36,0.9999997);
+  weight->SetBinContent(37,0.9999997);
+  weight->SetBinContent(38,0.9999996);
+  weight->SetBinContent(39,0.9999996);
+  weight->SetBinContent(40,0.9999995);
+  weight->SetBinContent(41,0.9999993);
+  weight->SetBinContent(42,1);
+  weight->SetBinContent(43,1);
+  weight->SetBinContent(44,1);
+  weight->SetBinContent(45,1);
+  weight->SetBinContent(46,1);
+  weight->SetBinContent(47,0.9999999);
+  weight->SetBinContent(48,0.9999998);
+  weight->SetBinContent(49,0.9999998);
+  weight->SetBinContent(50,0.9999999);
+  weight->SetBinContent(51,0.9999999);
+  weight->SetBinContent(52,0.9999999);
+  weight->SetBinContent(53,0.9999999);
+  weight->SetBinContent(54,0.9999998);
+  weight->SetBinContent(55,0.9999998);
+  weight->SetBinContent(56,0.9999998);
+  weight->SetBinContent(57,0.9999997);
+  weight->SetBinContent(58,0.9999996);
+  weight->SetBinContent(59,0.9999995);
+  weight->SetBinContent(60,1);
+  weight->SetBinContent(61,1);
+  weight->SetBinContent(62,1);
+  weight->SetBinContent(63,1);
+  weight->SetBinContent(64,1);
+  weight->SetBinContent(65,0.9999999);
+  weight->SetBinContent(66,0.9999998);
+  weight->SetBinContent(67,0.9999998);
+  weight->SetBinContent(68,0.9999999);
+  weight->SetBinContent(69,1);
+  weight->SetBinContent(70,1);
+  weight->SetBinContent(71,0.9999997);
+  weight->SetBinContent(72,0.9999995);
+  weight->SetBinContent(73,0.9999994);
+  weight->SetBinContent(74,1);
+  weight->SetBinContent(75,1);
+  weight->SetBinContent(76,1);
+  weight->SetBinContent(77,1);
+  weight->SetBinContent(78,0.9999999);
+  weight->SetBinContent(79,1);
+  weight->SetBinContent(80,1);
+  weight->SetEntries(526);
+  Int_t limit = weight->GetRandom();
+  pythia->SetTriggerChargedMultiplicity(limit, 1.4);
+  comment = comment.Append(Form("; multiplicity threshold set to %d in |eta| < 1.4", limit));
+  return pythia;
+}
+
+//______________________________________________________________________
+AliGenerator* PythiaJets()
+{
+      comment = comment.Append(" pp: Pythia low-pt");
+//
+//    Pythia
+      AliGenPythia* pythia = new AliGenPythia(-1); 
+      /*   pythia->SetMomentumRange(0, 999999.);
+      pythia->SetThetaRange(0., 180.);
+      pythia->SetYRange(-12.,12.);
+      pythia->SetPtRange(0,1000.);*/
+      pythia->SetProcess(kPyJets);
+      pythia->SetEnergyCMS(energy);
+      pythia->SetStrucFunc(kCTEQ6l);
+      //  pythia->SetPtHard(50, 1000);
+//
+      return pythia;
+}
+
+
+
+//___________________________________________________________________
+AliGenerator* MbPhojet()
+{
+  comment = comment.Append(" pp: Pythia low-pt");
+#if defined(__CINT__)
+  gSystem->Load("libdpmjet");      // Parton density functions
+  gSystem->Load("libTDPMjet");      // Parton density functions
+#endif
+  AliGenDPMjet* dpmjet = new AliGenDPMjet(-1); 
+
+  /*  dpmjet->SetMomentumRange(0, 999999.);
+  dpmjet->SetThetaRange(0., 180.);
+  dpmjet->SetYRange(-12.,12.);
+  dpmjet->SetPtRange(0,1000.);*/
+
+  dpmjet->SetProcess(kDpmMb);
+  dpmjet->SetEnergyCMS(energy);
+  return dpmjet;
+}
+
+//__________________________________________________________________
+AliGenerator* Hijing()
+{
+    AliGenHijing *gener = new AliGenHijing(-1);
+    // centre of mass energy 
+    gener->SetEnergyCMS(energy);
+    gener->SetImpactParameterRange(bMin, bMax);        
+    // reference frame
+    gener->SetReferenceFrame("CMS");
+    // projectile
+    gener->SetProjectile("A", 208, 82);
+    gener->SetTarget    ("A", 208, 82);
+    // tell hijing to keep the full parent child chain
+    gener->KeepFullEvent();
+    // enable jet quenching
+    gener->SetJetQuenching(1);
+    // enable shadowing
+    gener->SetDecaysOff(0);
+    gener->SetShadowing(1);
+    // Don't track spectators
+    gener->SetSpectators(0);
+    // kinematic selection
+    gener->SetSelectAll(0);
+    return gener;
+}
+
+
+//__________________________________________________________________
+AliGenerator* Hijing2000()
+{
+  AliGenHijing *gener = (AliGenHijing*) Hijing();
+  gener->SetJetQuenching(0);   
+  gener->SetPtHardMin (2.3);
+  return gener;
+}
+
+
+
+//_____________________________________________________________
+AliGenerator* Hijing_pA(Bool_t kSlowN) {
+  AliGenCocktail *gener = 0x0;
+  if (kSlowN) {
+    gener  = new AliGenCocktail();
+    gener->SetProjectile("A", 208, 82);
+    gener->SetTarget    ("P",   1,  1);
+    gener->SetEnergyCMS(energy);
+  }
+  
+  AliGenHijing   *hijing = new AliGenHijing(-1);
+  // centre of mass energy 
+  hijing->SetEnergyCMS(energy);
+  // impact parameter range
+  hijing->SetImpactParameterRange(0., 100.);
+  // reference frame
+  hijing->SetReferenceFrame("CMS");
+  hijing->SetBoostLHC(1);
+  // projectile
+  hijing->SetProjectile("A", 208, 82);
+  hijing->SetTarget    ("P", 1, 1);
+  // tell hijing to keep the full parent child chain
+  hijing->KeepFullEvent();
+  // enable jet quenching
+  hijing->SetJetQuenching(0);
+  // enable shadowing
+  hijing->SetShadowing(1);
+  // kinematic selection
+  hijing->SetSelectAll(0);
+  
+  if (!kSlowN) {
+    // DO track spectators
+    hijing->SetSpectators(1);
+    return hijing;
+  }
+  else {
+    // Cocktail with slow nucleons generator
+    // DO NOT track spectators
+    hijing->SetSpectators(0);
+    AliGenSlowNucleons* gray   = new AliGenSlowNucleons(1);
+    AliCollisionGeometry* coll = hijing->CollisionGeometry();
+    AliSlowNucleonModelExp* model = new AliSlowNucleonModelExp();
+    //  Not yet in the release...
+    //   model->SetSaturation(kTRUE);
+    gray->SetSlowNucleonModel(model);
+    gray->SetTarget(208,82);
+    gray->SetThetaDist(1);
+    gray->SetProtonDirection(1);
+    //   gray->SetDebug(1);
+    gray->SetNominalCmsEnergy(2.*pBeamEnergy);
+    gray->NeedsCollisionGeometry();
+    gray->SetCollisionGeometry(coll);
+    
+    gener->AddGenerator(hijing, "Hijing pPb", 1);
+    gener->AddGenerator(gray, "Gray Particles", 1);
+    
+    return gener;
+  }
+}
+
+
+
+//__________________________________________________________________
+AliGenerator* DPMjet()
+{
+  AliGenDPMjet* dpmjet = new AliGenDPMjet(-1); 
+  dpmjet->SetEnergyCMS(energy);
+  dpmjet->SetProjectile("A", 208, 82);
+  dpmjet->SetTarget    ("A", 208, 82);
+  dpmjet->SetImpactParameterRange(bMin, bMax);
+  dpmjet->SetPi0Decay(0);
+  return dpmjet;
+}
+
+
+//____________________________________________
+AliGenerator* DPMjet_pA(Bool_t fragments)
+{
+  AliGenDPMjet *gener = new AliGenDPMjet(-1);
+  //  A-p
+  gener->SetProjectile("A", 208, 82);
+  gener->SetTarget("P", 1, 1);
+  //
+  gener->SetEnergyCMS(energy);
+  gener->SetProjectileBeamEnergy(82.*pBeamEnergy/208.);
+
+  //
+  gener->SetProcess(kDpmMb);
+  gener->SetImpactParameterRange(0., 100.);
+  // DO NOT BOOST !
+  //  gener->SetBoostLHC(1);
+  //
+  gener->SetFragmentProd(fragments);
+  
+  return gener;
+
+}
+
+//__________________________________________________________________
+AliGenerator* AmptDefault()
+{
+  AliGenAmpt *genHi = new AliGenAmpt(-1);
+  //=========================================================================
+  // THE DECAYER
+  AliDecayer *decayer = new AliDecayerPythia();
+  cout << "*****************************************" << endl;
+  genHi->SetForceDecay( kHadronicD );
+  genHi->SetDecayer( decayer );
+  //=========================================================================
+  genHi->SetEnergyCMS(energy);
+  genHi->SetReferenceFrame("CMS");
+  genHi->SetProjectile("A",208,82);
+  genHi->SetTarget("A",208,82);
+    
+  genHi->SetIsoft(1); //1: defaul - 4: string melting
+  genHi->SetStringFrag(0.5,0.9); //Lund string frgmentation parameters
+  genHi->SetPtHardMin (3);
+  //genHi->SetImpactParameterRange(9.,9.5);
+  genHi->SetImpactParameterRange(bMin,bMax);
+    
+  // Xmu = 3.2 fm^-1 and as = 0.33 ==> sigma_{partonic} = 1.5mb
+  // Ntmax = 150
+  // v_2{2} = 0.102105 +/- 0.000266894
+  // v_2{4} = 0.0829477 +/- 0.00106158
+    
+  genHi->SetNtMax(150);        // NTMAX: number of timesteps (D=150)
+  genHi->SetXmu(3.2264);        // parton screening mass in fm^(-1) (D=3.2264d0)
+    
+  genHi->SetJetQuenching(0);  // enable jet quenching
+  genHi->SetShadowing(1);     // enable shadowing
+  genHi->SetDecaysOff(1);     // neutral pion and heavy particle decays switched off
+  genHi->SetSpectators(0);    // track spectators
+  //Boost into LHC lab frame
+  genHi->SetBoostLHC(1);
+  //  genHi->Init();
+  genHi->SetRandomReactionPlane(kTRUE);
+  return genHi;
+}
+
+//__________________________________________________________________
+AliGenerator* AmptStringMelting()
+{
+  AliGenAmpt *genHi = new AliGenAmpt(-1);
+  //=========================================================================
+  // THE DECAYER
+  AliDecayer *decayer = new AliDecayerPythia();
+  cout << "*****************************************" << endl;
+  genHi->SetForceDecay( kHadronicD );
+  genHi->SetDecayer( decayer );
+  //=========================================================================
+  genHi->SetEnergyCMS(energy);
+  genHi->SetReferenceFrame("CMS");
+  genHi->SetProjectile("A",208,82);
+  genHi->SetTarget("A",208,82);
+    
+  genHi->SetIsoft(4); //1: defaul - 4: string melting
+  genHi->SetStringFrag(0.5,0.9); //Lund string frgmentation parameters
+  genHi->SetPtHardMin (3);
+  //genHi->SetImpactParameterRange(9.,9.5);
+  genHi->SetImpactParameterRange(bMin,bMax);
+
+  // Xmu = 3.2 fm^-1 and as = 0.33 ==> sigma_{partonic} = 1.5mb
+  // Ntmax = 150
+  // v_2{2} = 0.102105 +/- 0.000266894
+  // v_2{4} = 0.0829477 +/- 0.00106158
+
+  genHi->SetNtMax(150);        // NTMAX: number of timesteps (D=150)
+  genHi->SetXmu(3.2264);        // parton screening mass in fm^(-1) (D=3.2264d0)
+
+  genHi->SetJetQuenching(0);  // enable jet quenching
+  genHi->SetShadowing(1);     // enable shadowing
+  genHi->SetDecaysOff(1);     // neutral pion and heavy particle decays switched off
+  genHi->SetSpectators(0);    // track spectators
+  //Boost into LHC lab frame
+  genHi->SetBoostLHC(1);
+//  genHi->Init();
+  genHi->SetRandomReactionPlane(kTRUE);
+  return genHi;
+
+}
+
+//__________________________________________________________________
+AliGenerator* AmptStringMeltingNoART()
+{
+    AliGenAmpt *genHi = new AliGenAmpt(-1);
+    //=========================================================================
+    // THE DECAYER
+    AliDecayer *decayer = new AliDecayerPythia();
+    cout << "*****************************************" << endl;
+    genHi->SetForceDecay( kHadronicD );
+    genHi->SetDecayer( decayer );
+    //=========================================================================
+    genHi->SetEnergyCMS(energy);
+    genHi->SetReferenceFrame("CMS");
+    genHi->SetProjectile("A",208,82);
+    genHi->SetTarget("A",208,82);
+    
+    genHi->SetIsoft(4); //1: defaul - 4: string melting
+    genHi->SetStringFrag(0.5,0.9); //Lund string frgmentation parameters
+    genHi->SetPtHardMin (3);
+    //genHi->SetImpactParameterRange(9.,9.5);
+    genHi->SetImpactParameterRange(bMin,bMax);
+    
+    // Xmu = 3.2 fm^-1 and as = 0.33 ==> sigma_{partonic} = 1.5mb
+    // Ntmax = 150
+    // v_2{2} = 0.102105 +/- 0.000266894
+    // v_2{4} = 0.0829477 +/- 0.00106158
+    
+    genHi->SetNtMax(3);        // NTMAX: number of timesteps (D=150)
+    genHi->SetXmu(3.2264);        // parton screening mass in fm^(-1) (D=3.2264d0)
+    
+    genHi->SetJetQuenching(0);  // enable jet quenching
+    genHi->SetShadowing(1);     // enable shadowing
+    genHi->SetDecaysOff(1);     // neutral pion and heavy particle decays switched off
+    genHi->SetSpectators(0);    // track spectators
+    //Boost into LHC lab frame
+    genHi->SetBoostLHC(1);
+    //  genHi->Init();
+    genHi->SetRandomReactionPlane(kTRUE);
+    return genHi;
+    
+}
+
+
+//__________________________________________________________________
+AliGenerator* AmptReducedFlow()
+{
+  AliGenAmpt *genHi = new AliGenAmpt(-1);
+  //=========================================================================
+  // THE DECAYER
+  AliDecayer *decayer = new AliDecayerPythia();
+  cout << "*****************************************" << endl;
+  genHi->SetForceDecay( kHadronicD );
+  genHi->SetDecayer( decayer );
+  //=========================================================================
+  genHi->SetEnergyCMS(energy);
+  genHi->SetReferenceFrame("CMS");
+  genHi->SetProjectile("A",208,82);
+  genHi->SetTarget("A",208,82);
+    
+  genHi->SetIsoft(4); //1: defaul - 4: string melting
+  genHi->SetStringFrag(0.5,0.9); //Lund string frgmentation parameters
+  genHi->SetPtHardMin (3);
+  //genHi->SetImpactParameterRange(9.,9.5);
+  genHi->SetImpactParameterRange(bMin,bMax);
+
+  // Xmu = 12.4 fm^-1 and as = 0.33 ==> sigma_{partonic} = 0.1mb
+  // Ntmax = 20
+  // flow estimates from Q-cumulants
+  // (POI, without weights)
+  // v_2{2} = 0.0549735 +/- 0.000270249
+  // v_2{4} = 0.0421905 +/- 0.00189449
+
+  genHi->SetNtMax(20);        // NTMAX: number of timesteps (D=150)
+  genHi->SetXmu(12.4);        // parton screening mass in fm^(-1) (D=3.2264d0)
+
+  genHi->SetJetQuenching(0);  // enable jet quenching
+  genHi->SetShadowing(1);     // enable shadowing
+  genHi->SetDecaysOff(1);     // neutral pion and heavy particle decays switched off
+  genHi->SetSpectators(0);    // track spectators
+  //Boost into LHC lab frame
+  genHi->SetBoostLHC(1);
+ // genHi->Init();
+  genHi->SetRandomReactionPlane(kTRUE);
+  return genHi;
+
+}
+
+//__________________________________________________________________
+AliGenerator* AmptpA()
+{
+    AliGenAmpt *genHi = new AliGenAmpt(-1);
+    //=========================================================================
+    // THE DECAYER
+    AliDecayer *decayer = new AliDecayerPythia();
+    cout << "*****************************************" << endl;
+    genHi->SetForceDecay( kHadronicD );
+    genHi->SetDecayer( decayer );
+    //=========================================================================
+    genHi->SetEnergyCMS(energy);
+    genHi->SetReferenceFrame("CMS");
+    genHi->SetProjectile("A", 208, 82);
+    genHi->SetTarget    ("P", 1, 1);
+    genHi->SetIsoft(4); //1: defaul - 4: string melting
+    genHi->SetStringFrag(0.5,0.9); //Lund string frgmentation parameters
+    genHi->SetPtHardMin (3);
+    //genHi->SetImpactParameterRange(9.,9.5);
+    genHi->SetImpactParameterRange(bMin,bMax);
+    genHi->SetNtMax(1500); //NTMAX: number of timesteps (D=150)
+    genHi->SetXmu(3.2264); //parton screening mass in fm^(-1) (D=3.2264d0)
+    genHi->SetJetQuenching(0); // enable jet quenching
+    genHi->SetShadowing(1);    // enable shadowing
+    genHi->SetDecaysOff(1);    // neutral pion and heavy particle decays switched off
+    genHi->SetSpectators(0);   // track spectators
+    //Boost into LHC lab frame
+    genHi->SetBoostLHC(1);
+    //  genHi->Init();
+    genHi->SetRandomReactionPlane(kTRUE);
+    return genHi;
+    
+}