/***************************************************************************
 *  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;
    
}