[u/mrichter/AliRoot.git] / macros / Config_PythiaHeavyFlavours.C

//*******************************************************************//
// Configuration file for charm / beauty generation with PYTHIA      //
//                                                                   //
// The parameters have been tuned in order to reproduce the inclusive//
// heavy quark pt distribution given by the NLO pQCD calculation by  //
// Mangano, Nason and Ridolfi.                                       //
//                                                                   //
// For details and for the NORMALIZATION of the yields see:          //
//   N.Carrer and A.Dainese,                                         //
//   "Charm and beauty production at the LHC",                       //
//   ALICE-INT-2003-019, [arXiv:hep-ph/0311225];                     //
//   PPR Chapter 6.6, CERN/LHCC 2005-030 (2005).                     //
//*******************************************************************//
#if !defined(__CINT__) || defined(__MAKECINT__)
#include <Riostream.h>
#include <TRandom.h>
#include <TDatime.h>
#include <TSystem.h>
#include <TVirtualMC.h>
#include <TGeant3TGeo.h>
#include "STEER/AliRunLoader.h"
#include "STEER/AliRun.h"
#include "STEER/AliConfig.h"
#include "PYTHIA6/AliDecayerPythia.h"
#include "PYTHIA6/AliGenPythia.h"
#include "STEER/AliMagF.h"
#include "STRUCT/AliBODY.h"
#include "STRUCT/AliMAG.h"
#include "STRUCT/AliABSOv3.h"
#include "STRUCT/AliDIPOv3.h"
#include "STRUCT/AliHALLv3.h"
#include "STRUCT/AliFRAMEv2.h"
#include "STRUCT/AliSHILv3.h"
#include "STRUCT/AliPIPEv3.h"
#include "ITS/AliITSv11Hybrid.h"
#include "TPC/AliTPCv2.h"
#include "TOF/AliTOFv6T0.h"
#include "HMPID/AliHMPIDv3.h"
#include "ZDC/AliZDCv3.h"
#include "TRD/AliTRDv1.h"
#include "FMD/AliFMDv0.h"
#include "MUON/AliMUONv1.h"
#include "PHOS/AliPHOSv1.h"
#include "PMD/AliPMDv1.h"
#include "T0/AliT0v1.h"
#include "ACORDE/AliACORDEv1.h"
#endif

//--- Heavy Flavour Production ---
enum ProcessHvFl_t 
{
  kCharmPbPb5500,  kCharmpPb8800,  kCharmpp14000,  kCharmpp14000wmi,
  kCharmSemiElpp14000wmi,
  kD0PbPb5500,     kD0pPb8800,     kD0pp14000,
  kDPlusPbPb5500,  kDPluspPb8800,  kDPluspp14000,
  kDPlusStrangePbPb5500, kDPlusStrangepPb8800, kDPlusStrangepp14000,
  kBeautyPbPb5500, kBeautypPb8800, kBeautypp14000, kBeautypp14000wmi,
  kBeautySemiElpp14000wmi
};
//--- Decay Mode ---
enum DecayHvFl_t 
{
  kNature,  kHadr, kSemiEl, kSemiMu
};
//--- Rapidity Cut ---
enum YCut_t
{
  kFull, kBarrel, kMuonArm
};
//--- Trigger config ---
enum TrigConf_t
{
    kDefaultPPTrig, kDefaultPbPbTrig
};
const char * TrigConfName[] = {
    "p-p","Pb-Pb"
};

//--- Functions ---
class AliGenPythia;
AliGenPythia *PythiaHVQ(ProcessHvFl_t proc);
void          LoadPythia();


// This part for configuration
static ProcessHvFl_t procHvFl = kCharmpp14000wmi;
static DecayHvFl_t   decHvFl  = kNature; 
static YCut_t        ycut     = kFull;
static AliMagF::BMap_t mag    = AliMagF::k5kG; 
static TrigConf_t    trig     = kDefaultPbPbTrig; // default PbPb trigger configuration
// nEvts = -1  : you get 1 QQbar pair and all the fragmentation and 
//               decay chain
// nEvts = N>0 : you get N charm / beauty Hadrons 
Int_t nEvts = -1; 
// stars = kTRUE : all heavy resonances and their decay stored
//       = kFALSE: only final heavy hadrons and their decays stored
Bool_t stars = kTRUE;

// To be used only with kCharmppMNRwmi and kBeautyppMNRwmi
// To get a "reasonable" agreement with MNR results, events have to be 
// generated with the minimum ptHard set to 2.76 GeV.
// To get a "perfect" agreement with MNR results, events have to be 
// generated in four ptHard bins with the following relative 
// normalizations:
//  CHARM
// 2.76-3 GeV: 25%
//    3-4 GeV: 40%
//    4-8 GeV: 29%
//     >8 GeV:  6%
//  BEAUTY
// 2.76-4 GeV:  5% 
//    4-6 GeV: 31%
//    6-8 GeV: 28%
//     >8 GeV: 36%
Float_t ptHardMin =  2.76;
Float_t ptHardMax = -1.;


// Comment line
static TString comment;

void Config()
{
 
  //========================//
  // Set Random Number seed //
  //========================//
  TDatime dt;
  UInt_t curtime=dt.Get();
  UInt_t procid=gSystem->GetPid();
  UInt_t seed=curtime-procid;

  //  gRandom->SetSeed(seed);
  //  cerr<<"Seed for random number generation= "<<seed<<endl; 
  gRandom->SetSeed(12345);
  
  // Load Pythia libraries
  LoadPythia();
  // libraries required by geant321
#if defined(__CINT__)
  gSystem->Load("libgeant321");
#endif

  new TGeant3TGeo("C++ Interface to Geant3");

  if(!AliCDBManager::Instance()->IsDefaultStorageSet()){
    AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
    AliCDBManager::Instance()->SetRun(0);
  }
  
  //=======================================================================
  //  Create the output file

   
  AliRunLoader* rl=0x0;

  cout<<"Config.C: Creating Run Loader ..."<<endl;
  rl = AliRunLoader::Open("galice.root",
			  AliConfig::GetDefaultEventFolderName(),
			  "recreate");
  if (rl == 0x0)
    {
      gAlice->Fatal("Config.C","Can not instatiate the Run Loader");
      return;
    }
  rl->SetCompressionLevel(2);
  rl->SetNumberOfEventsPerFile(1000);
  gAlice->SetRunLoader(rl);
  // gAlice->SetGeometryFromFile("geometry.root");
  // gAlice->SetGeometryFromCDB();

  // Set the trigger configuration
  AliSimulation::Instance()->SetTriggerConfig(TrigConfName[trig]);
  cout<<"Trigger configuration is set to  "<<TrigConfName[trig]<<endl;

  //
  //=======================================================================
  // ************* STEERING parameters FOR ALICE SIMULATION **************
  // --- Specify event type to be tracked through the ALICE setup
  // --- All positions are in cm, angles in degrees, and P and E in GeV


    gMC->SetProcess("DCAY",1);
    gMC->SetProcess("PAIR",1);
    gMC->SetProcess("COMP",1);
    gMC->SetProcess("PHOT",1);
    gMC->SetProcess("PFIS",0);
    gMC->SetProcess("DRAY",0);
    gMC->SetProcess("ANNI",1);
    gMC->SetProcess("BREM",1);
    gMC->SetProcess("MUNU",1);
    gMC->SetProcess("CKOV",1);
    gMC->SetProcess("HADR",1);
    gMC->SetProcess("LOSS",2);
    gMC->SetProcess("MULS",1);
    gMC->SetProcess("RAYL",1);

    Float_t cut = 1.e-3;        // 1MeV cut by default
    Float_t tofmax = 1.e10;

    gMC->SetCut("CUTGAM", cut);
    gMC->SetCut("CUTELE", cut);
    gMC->SetCut("CUTNEU", cut);
    gMC->SetCut("CUTHAD", cut);
    gMC->SetCut("CUTMUO", cut);
    gMC->SetCut("BCUTE",  cut); 
    gMC->SetCut("BCUTM",  cut); 
    gMC->SetCut("DCUTE",  cut); 
    gMC->SetCut("DCUTM",  cut); 
    gMC->SetCut("PPCUTM", cut);
    gMC->SetCut("TOFMAX", tofmax); 


  // Set External decayer //
  //======================//
  TVirtualMCDecayer* decayer = new AliDecayerPythia();
  // DECAYS
  //
  switch(decHvFl) {
  case kNature:
    decayer->SetForceDecay(kAll);
    break;
  case kHadr:
    decayer->SetForceDecay(kHadronicD);
    break;
  case kSemiEl:
    decayer->SetForceDecay(kSemiElectronic);
    break;
  case kSemiMu:
    decayer->SetForceDecay(kSemiMuonic);
    break;
  }
  decayer->Init();
  gMC->SetExternalDecayer(decayer);

  //=========================//
  // Generator Configuration //
  //=========================//
  AliGenPythia *pythia = PythiaHVQ(procHvFl);
  // FeedDown option
  pythia->SetFeedDownHigherFamily(kFALSE);
  // Stack filling option
  if(!stars) pythia->SetStackFillOpt(AliGenPythia::kParentSelection);
  // Set Count mode
  if(nEvts>0) pythia->SetCountMode(AliGenPythia::kCountParents);
  // DECAYS
  //
  switch(decHvFl) {
  case kNature:
    pythia->SetForceDecay(kAll);
    break;
  case kHadr:
    pythia->SetForceDecay(kHadronicD);
    break;
  case kSemiEl:
    pythia->SetForceDecay(kSemiElectronic);
    break;
  case kSemiMu:
    pythia->SetForceDecay(kSemiMuonic);
    break;
  }
  // GEOM & KINE CUTS
  //
  pythia->SetMomentumRange(0,99999999);
  pythia->SetPhiRange(0., 360.);
  pythia->SetThetaRange(0,180);
  switch(ycut) {
  case kFull:
    pythia->SetYRange(-999,999);
    break;
  case kBarrel:
    pythia->SetYRange(-2,2);
    break;
  case kMuonArm:
    pythia->SetYRange(1,6);
    break;
  }

  // PRIMARY VERTEX
  //
  pythia->SetOrigin(0, 0, 0);    // vertex position
  pythia->SetSigma(0, 0, 5.3);   // Sigma in (X,Y,Z) (cm) on IP position
  pythia->SetCutVertexZ(1.);     // Truncate at 1 sigma
  pythia->SetVertexSmear(kPerEvent); 


  pythia->SetTrackingFlag(0);
  // Specify GEANT tracking limits (Rmax, Zmax)
  //gAlice->GetMCApp()->TrackingLimits(90.,1.0e10);


  pythia->Init();

  // FIELD
  //    
  if (mag == AliMagF::k2kG) {
    comment = comment.Append(" | L3 field 0.2 T");
  } else if (mag == AliMagF::k5kG) {
    comment = comment.Append(" | L3 field 0.5 T");
  }
  printf("\n \n Comment: %s \n \n", comment.Data());
    
  AliMagF* field = new AliMagF("Maps","Maps",  1., 1.,  mag);
  field->SetL3ConstField(0); //Using const. field in the barrel
  TGeoGlobalMagField::Instance()->SetField(field);

  rl->CdGAFile();

  // By default all ALICE is switched off
  Int_t iABSO=0;
  Int_t iACORDE=0;
  Int_t iDIPO=0;
  Int_t iFMD=0;
  Int_t iFRAME=0;
  Int_t iHALL=0;
  Int_t iITS=0;
  Int_t iMAG=0;
  Int_t iMUON=0;
  Int_t iPHOS=0;
  Int_t iPIPE=0;
  Int_t iPMD=0;
  Int_t iHMPID=0;
  Int_t iSHIL=0;
  Int_t iT0=0;
  Int_t iTOF=0;
  Int_t iTPC=0;
  Int_t iTRD=0;
  Int_t iZDC=0;

  //=================== Alice BODY parameters =============================
  AliBODY *BODY = new AliBODY("BODY","Alice envelop");

  if(iMAG) {
    //=================== MAG parameters ============================
    // --- Start with Magnet since detector layouts may be depending ---
    // --- on the selected Magnet dimensions ---
    AliMAG *MAG  = new AliMAG("MAG","Magnet");
  }


  if(iABSO) {
    //=================== ABSO parameters ============================
    AliABSO *ABSO  = new AliABSOv3("ABSO","Muon Absorber");
  }

  if(iDIPO) {
    //=================== DIPO parameters ============================

    AliDIPO *DIPO  = new AliDIPOv3("DIPO","Dipole version 3");
  }

  if(iHALL) {
    //=================== HALL parameters ============================

    AliHALL *HALL  = new AliHALLv3("HALL","Alice Hall");
  }


  if(iFRAME) {
    //=================== FRAME parameters ============================

    AliFRAME *FRAME  = new AliFRAMEv2("FRAME","Space Frame");
    FRAME->SetHoles(1);
  }

  if(iSHIL) {
    //=================== SHIL parameters ============================

    AliSHIL *SHIL  = new AliSHILv3("SHIL","Shielding");
  }


  if(iPIPE) {
    //=================== PIPE parameters ============================

    AliPIPE *PIPE  = new AliPIPEv3("PIPE","Beam Pipe");
  }


  if(iITS) {
    //=================== ITS parameters ============================

    AliITS *ITS  = new AliITSv11Hybrid("ITS","ITS v11Hybrid");
  }
  

  if(iTPC) {
    //============================ TPC parameters ===================
    AliTPC *TPC  = new AliTPCv2("TPC","Default");
  }


  if(iTOF) {
    //=================== TOF parameters ============================
    AliTOF *TOF  = new AliTOFv6T0("TOF","normal TOF");
  }

  if(iHMPID) {
    //=================== HMPID parameters ===========================
    AliHMPID *HMPID  = new AliHMPIDv3("HMPID","normal HMPID");    

  }


  if(iZDC) {
    //=================== ZDC parameters ============================

    AliZDC *ZDC  = new AliZDCv3("ZDC","normal ZDC");
  }

  if(iACORDE) {
    //=================== ACORDE parameters ============================

    AliACORDE *ACORDE  = new AliACORDEv1("ACORDE","normal ACORDE");
  }

  if(iTRD) {
    //=================== TRD parameters ============================
  
    AliTRD *TRD  = new AliTRDv1("TRD","TRD slow simulator");
  }

  if(iFMD) {
    //=================== FMD parameters ============================

    AliFMD *FMD  = new AliFMDv0("FMD","normal FMD");
  }

  if(iMUON) {
    //=================== MUON parameters ===========================
    AliMUON *MUON  = new AliMUONv1("MUON","default");
  }
 
  //=================== PHOS parameters ===========================

  if(iPHOS) {
    AliPHOS *PHOS  = new AliPHOSv1("PHOS","GPS2");
  }


  //=================== ACORDE parameters ===========================

  if(iACORDE) {
    AliACORDE *ACORDE  = new AliACORDEv1("ACORDE","Normal ACORDEGPS2");
  }


  if(iPMD) {
    //=================== PMD parameters ============================

    AliPMD *PMD  = new AliPMDv1("PMD","normal PMD");
    PMD->SetPAR(1., 1., 0.8, 0.02);
    PMD->SetIN(6., 18., -580., 27., 27.);
    PMD->SetGEO(0.0, 0.2, 4.);
    PMD->SetPadSize(0.8, 1.0, 1.0, 1.5);

  }

  if(iT0) {
    //=================== T0 parameters ============================
    AliT0 *T0  = new AliT0v1("T0","T0 Detector");
  }

         
}
//
//           PYTHIA
//
AliGenPythia *PythiaHVQ(ProcessHvFl_t proc) {

  AliGenPythia * gener = 0x0;

  switch(proc) {
  case kCharmPbPb5500:
    comment = comment.Append(" Charm in Pb-Pb at 5.5 TeV");
    gener = new AliGenPythia(nEvts);
    gener->SetProcess(kPyCharmPbPbMNR);
    gener->SetStrucFunc(kCTEQ4L);
    gener->SetPtHard(2.1,-1.0);
    gener->SetEnergyCMS(5500.);
    gener->SetNuclei(208,208);
    break;
  case kCharmpPb8800:
    comment = comment.Append(" Charm in p-Pb at 8.8 TeV");
    gener = new AliGenPythia(nEvts);
    gener->SetProcess(kPyCharmpPbMNR);
    gener->SetStrucFunc(kCTEQ4L);
    gener->SetPtHard(2.1,-1.0);
    gener->SetEnergyCMS(8800.);
    gener->SetProjectile("P",1,1);
    gener->SetTarget("Pb",208,82);
    break;
  case kCharmpp14000:
    comment = comment.Append(" Charm in pp at 14 TeV");
    gener = new AliGenPythia(nEvts);
    gener->SetProcess(kPyCharmppMNR);
    gener->SetStrucFunc(kCTEQ4L);
    gener->SetPtHard(2.1,-1.0);
    gener->SetEnergyCMS(14000.);
    break;
  case kCharmpp14000wmi:
    comment = comment.Append(" Charm in pp at 14 TeV with mult. interactions");
    gener = new AliGenPythia(-1);
    gener->SetProcess(kPyCharmppMNRwmi);
    gener->SetStrucFunc(kCTEQ5L);
    gener->SetPtHard(ptHardMin,ptHardMax);
    gener->SetEnergyCMS(14000.);
    break;
  case kCharmSemiElpp14000wmi:
    comment = comment.Append(" Charm in pp at 14 TeV with mult. interactions");
    gener = new AliGenPythia(-1);
    gener->SetProcess(kPyCharmppMNRwmi);
    gener->SetStrucFunc(kCTEQ5L);
    gener->SetPtHard(ptHardMin,ptHardMax);
    gener->SetEnergyCMS(14000.);
    gener->SetCutOnChild(1);
    gener->SetPdgCodeParticleforAcceptanceCut(11);
    gener->SetChildYRange(-1,1);
    gener->SetChildPtRange(1,10000.);
    break;
  case kD0PbPb5500:
    comment = comment.Append(" D0 in Pb-Pb at 5.5 TeV");
    gener = new AliGenPythia(nEvts);
    gener->SetProcess(kPyD0PbPbMNR);
    gener->SetStrucFunc(kCTEQ4L);
    gener->SetPtHard(2.1,-1.0);
    gener->SetEnergyCMS(5500.);
    gener->SetNuclei(208,208);
    break;
  case kD0pPb8800:
    comment = comment.Append(" D0 in p-Pb at 8.8 TeV");
    gener = new AliGenPythia(nEvts);
    gener->SetProcess(kPyD0pPbMNR);
    gener->SetStrucFunc(kCTEQ4L);
    gener->SetPtHard(2.1,-1.0);
    gener->SetEnergyCMS(8800.);
    gener->SetProjectile("P",1,1);
    gener->SetTarget("Pb",208,82);
    break;
  case kD0pp14000:
    comment = comment.Append(" D0 in pp at 14 TeV");
    gener = new AliGenPythia(nEvts);
    gener->SetProcess(kPyD0ppMNR);
    gener->SetStrucFunc(kCTEQ4L);
    gener->SetPtHard(2.1,-1.0);
    gener->SetEnergyCMS(14000.);
    break;
  case kDPlusPbPb5500:
    comment = comment.Append(" DPlus in Pb-Pb at 5.5 TeV");
    gener = new AliGenPythia(nEvts);
    gener->SetProcess(kPyDPlusPbPbMNR);
    gener->SetStrucFunc(kCTEQ4L);
    gener->SetPtHard(2.1,-1.0);
    gener->SetEnergyCMS(5500.);
    gener->SetNuclei(208,208);
    break;
  case kDPluspPb8800:
    comment = comment.Append(" DPlus in p-Pb at 8.8 TeV");
    gener = new AliGenPythia(nEvts);
    gener->SetProcess(kPyDPluspPbMNR);
    gener->SetStrucFunc(kCTEQ4L);
    gener->SetPtHard(2.1,-1.0);
    gener->SetEnergyCMS(8800.);
    gener->SetProjectile("P",1,1);
    gener->SetTarget("Pb",208,82);
    break;
  case kDPluspp14000:
    comment = comment.Append(" DPlus in pp at 14 TeV");
    gener = new AliGenPythia(nEvts);
    gener->SetProcess(kPyDPlusppMNR);
    gener->SetStrucFunc(kCTEQ4L);
    gener->SetPtHard(2.1,-1.0);
    gener->SetEnergyCMS(14000.);
    break;
 case kDPlusStrangePbPb5500:
    comment = comment.Append(" DPlusStrange in Pb-Pb at 5.5 TeV");
    gener = new AliGenPythia(nEvts);
    gener->SetProcess(kPyDPlusStrangePbPbMNR);
    gener->SetStrucFunc(kCTEQ4L);
    gener->SetPtHard(2.1,-1.0);
    gener->SetEnergyCMS(5500.);
    gener->SetNuclei(208,208);
    break;
  case kDPlusStrangepPb8800:
    comment = comment.Append(" DPlusStrange in p-Pb at 8.8 TeV");
    gener = new AliGenPythia(nEvts);
    gener->SetProcess(kPyDPlusStrangepPbMNR);
    gener->SetStrucFunc(kCTEQ4L);
    gener->SetPtHard(2.1,-1.0);
    gener->SetEnergyCMS(8800.);
    gener->SetProjectile("P",1,1);
    gener->SetTarget("Pb",208,82);
    break;
  case kDPlusStrangepp14000:
    comment = comment.Append(" DPlusStrange in pp at 14 TeV");
    gener = new AliGenPythia(nEvts);
    gener->SetProcess(kPyDPlusStrangeppMNR);
    gener->SetStrucFunc(kCTEQ4L);
    gener->SetPtHard(2.1,-1.0);
    gener->SetEnergyCMS(14000.);
    break;
  case kBeautyPbPb5500:
    comment = comment.Append(" Beauty in Pb-Pb at 5.5 TeV");
    gener = new AliGenPythia(nEvts);
    gener->SetProcess(kPyBeautyPbPbMNR);
    gener->SetStrucFunc(kCTEQ4L);
    gener->SetPtHard(2.75,-1.0);
    gener->SetEnergyCMS(5500.);
    gener->SetNuclei(208,208);
    break;
  case kBeautypPb8800:
    comment = comment.Append(" Beauty in p-Pb at 8.8 TeV");
    gener = new AliGenPythia(nEvts);
    gener->SetProcess(kPyBeautypPbMNR);
    gener->SetStrucFunc(kCTEQ4L);
    gener->SetPtHard(2.75,-1.0);
    gener->SetEnergyCMS(8800.);
    gener->SetProjectile("P",1,1);
    gener->SetTarget("Pb",208,82);
    break;
  case kBeautypp14000:
    comment = comment.Append(" Beauty in pp at 14 TeV");
    gener = new AliGenPythia(nEvts);
    gener->SetProcess(kPyBeautyppMNR);
    gener->SetStrucFunc(kCTEQ4L);
    gener->SetPtHard(2.75,-1.0);
    gener->SetEnergyCMS(14000.);
    break;
  case kBeautypp14000wmi:
    comment = comment.Append(" Beauty in pp at 14 TeV with mult. interactions");
    gener = new AliGenPythia(-1);
    gener->SetProcess(kPyBeautyppMNRwmi);
    gener->SetStrucFunc(kCTEQ5L);
    gener->SetPtHard(ptHardMin,ptHardMax);
    gener->SetEnergyCMS(14000.);
    break;
  case kBeautySemiElpp14000wmi:
    comment = comment.Append(" Beauty in pp at 14 TeV with mult. interactions");
    gener = new AliGenPythia(-1);
    gener->SetProcess(kPyBeautyppMNRwmi);
    gener->SetStrucFunc(kCTEQ5L);
    gener->SetPtHard(ptHardMin,ptHardMax);
    gener->SetEnergyCMS(14000.);
    gener->SetCutOnChild(1);
    gener->SetPdgCodeParticleforAcceptanceCut(11);
    gener->SetChildYRange(-1,1);
    gener->SetChildPtRange(1,10000.);
    break;
  }

  return gener;
}


void LoadPythia()
{
    // Load Pythia related libraries
    gSystem->Load("liblhapdf.so");      // Parton density functions
    gSystem->Load("libEGPythia6.so");   // TGenerator interface
    gSystem->Load("libpythia6.so");     // Pythia
    gSystem->Load("libAliPythia6.so");  // ALICE specific implementations
}
Commit	Line	Data
6ccd5a72	1	//*******************************************************************//
	2	// Configuration file for charm / beauty generation with PYTHIA //
	3	// //
	4	// The parameters have been tuned in order to reproduce the inclusive//
	5	// heavy quark pt distribution given by the NLO pQCD calculation by //
	6	// Mangano, Nason and Ridolfi. //
	7	// //
	8	// For details and for the NORMALIZATION of the yields see: //
72fc043d	9	// N.Carrer and A.Dainese, //
72fc043d	10	// "Charm and beauty production at the LHC", //
09977179	11	// ALICE-INT-2003-019, [arXiv:hep-ph/0311225]; //
09977179	12	// PPR Chapter 6.6, CERN/LHCC 2005-030 (2005). //
6ccd5a72	13	//*******************************************************************//
de78d77f	14	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
	15	#include <Riostream.h>
	16	#include <TRandom.h>
	17	#include <TDatime.h>
	18	#include <TSystem.h>
	19	#include <TVirtualMC.h>
dcec124d	20	#include <TGeant3TGeo.h>
de78d77f	21	#include "STEER/AliRunLoader.h"
	22	#include "STEER/AliRun.h"
	23	#include "STEER/AliConfig.h"
	24	#include "PYTHIA6/AliDecayerPythia.h"
	25	#include "PYTHIA6/AliGenPythia.h"
f7a1cc68	26	#include "STEER/AliMagF.h"
de78d77f	27	#include "STRUCT/AliBODY.h"
de78d77f	28	#include "STRUCT/AliMAG.h"
96fb3353	29	#include "STRUCT/AliABSOv3.h"
	30	#include "STRUCT/AliDIPOv3.h"
	31	#include "STRUCT/AliHALLv3.h"
de78d77f	32	#include "STRUCT/AliFRAMEv2.h"
96fb3353	33	#include "STRUCT/AliSHILv3.h"
96fb3353	34	#include "STRUCT/AliPIPEv3.h"
99887042	35	#include "ITS/AliITSv11Hybrid.h"
de78d77f	36	#include "TPC/AliTPCv2.h"
0d7ef405	37	#include "TOF/AliTOFv6T0.h"
99887042	38	#include "HMPID/AliHMPIDv3.h"
9d07603a	39	#include "ZDC/AliZDCv3.h"
de78d77f	40	#include "TRD/AliTRDv1.h"
	41	#include "FMD/AliFMDv0.h"
	42	#include "MUON/AliMUONv1.h"
	43	#include "PHOS/AliPHOSv1.h"
	44	#include "PMD/AliPMDv1.h"
ababa197	45	#include "T0/AliT0v1.h"
b384f8a4	46	#include "ACORDE/AliACORDEv1.h"
de78d77f	47	#endif
6ccd5a72	48
	49	//--- Heavy Flavour Production ---
	50	enum ProcessHvFl_t
	51	{
09977179	52	kCharmPbPb5500, kCharmpPb8800, kCharmpp14000, kCharmpp14000wmi,
7e35f2ed	53	kCharmSemiElpp14000wmi,
6ccd5a72	54	kD0PbPb5500, kD0pPb8800, kD0pp14000,
09977179	55	kDPlusPbPb5500, kDPluspPb8800, kDPluspp14000,
84388d23	56	kDPlusStrangePbPb5500, kDPlusStrangepPb8800, kDPlusStrangepp14000,
7e35f2ed	57	kBeautyPbPb5500, kBeautypPb8800, kBeautypp14000, kBeautypp14000wmi,
7e35f2ed	58	kBeautySemiElpp14000wmi
6ccd5a72	59	};
	60	//--- Decay Mode ---
	61	enum DecayHvFl_t
	62	{
	63	kNature, kHadr, kSemiEl, kSemiMu
	64	};
	65	//--- Rapidity Cut ---
	66	enum YCut_t
	67	{
	68	kFull, kBarrel, kMuonArm
	69	};
dcec124d	70	//--- Trigger config ---
	71	enum TrigConf_t
	72	{
	73	kDefaultPPTrig, kDefaultPbPbTrig
	74	};
	75	const char * TrigConfName[] = {
	76	"p-p","Pb-Pb"
	77	};
6b612c2d	78
de78d77f	79	//--- Functions ---
1faeb480	80	class AliGenPythia;
de78d77f	81	AliGenPythia *PythiaHVQ(ProcessHvFl_t proc);
ee621efa	82	void LoadPythia();
de78d77f	83
6ccd5a72	84
6ccd5a72	85	// This part for configuration
7e35f2ed	86	static ProcessHvFl_t procHvFl = kCharmpp14000wmi;
6ccd5a72	87	static DecayHvFl_t decHvFl = kNature;
6ccd5a72	88	static YCut_t ycut = kFull;
f7a1cc68	89	static AliMagF::BMap_t mag = AliMagF::k5kG;
dcec124d	90	static TrigConf_t trig = kDefaultPbPbTrig; // default PbPb trigger configuration
6ccd5a72	91	// nEvts = -1 : you get 1 QQbar pair and all the fragmentation and
	92	// decay chain
	93	// nEvts = N>0 : you get N charm / beauty Hadrons
	94	Int_t nEvts = -1;
	95	// stars = kTRUE : all heavy resonances and their decay stored
	96	// = kFALSE: only final heavy hadrons and their decays stored
	97	Bool_t stars = kTRUE;
	98
09977179	99	// To be used only with kCharmppMNRwmi and kBeautyppMNRwmi
	100	// To get a "reasonable" agreement with MNR results, events have to be
	101	// generated with the minimum ptHard set to 2.76 GeV.
	102	// To get a "perfect" agreement with MNR results, events have to be
	103	// generated in four ptHard bins with the following relative
	104	// normalizations:
	105	// CHARM
	106	// 2.76-3 GeV: 25%
	107	// 3-4 GeV: 40%
	108	// 4-8 GeV: 29%
	109	// >8 GeV: 6%
	110	// BEAUTY
	111	// 2.76-4 GeV: 5%
	112	// 4-6 GeV: 31%
	113	// 6-8 GeV: 28%
	114	// >8 GeV: 36%
	115	Float_t ptHardMin = 2.76;
	116	Float_t ptHardMax = -1.;
	117
6ccd5a72	118
	119	// Comment line
	120	static TString comment;
	121
	122	void Config()
	123	{
	124
	125	//========================//
	126	// Set Random Number seed //
	127	//========================//
	128	TDatime dt;
	129	UInt_t curtime=dt.Get();
	130	UInt_t procid=gSystem->GetPid();
	131	UInt_t seed=curtime-procid;
	132
de78d77f	133	// gRandom->SetSeed(seed);
72fc043d	134	// cerr<<"Seed for random number generation= "<<seed<<endl;
6ccd5a72	135	gRandom->SetSeed(12345);
72fc043d	136
ee621efa	137	// Load Pythia libraries
ee621efa	138	LoadPythia();
6ccd5a72	139	// libraries required by geant321
de78d77f	140	#if defined(__CINT__)
6ccd5a72	141	gSystem->Load("libgeant321");
de78d77f	142	#endif
6ccd5a72	143
b0470150	144	new TGeant3TGeo("C++ Interface to Geant3");
6ccd5a72	145
458041af	146	if(!AliCDBManager::Instance()->IsDefaultStorageSet()){
162637e4	147	AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
458041af	148	AliCDBManager::Instance()->SetRun(0);
	149	}
	150
6ccd5a72	151	//=======================================================================
6ccd5a72	152	// Create the output file
de78d77f	153
	154
	155	AliRunLoader* rl=0x0;
	156
	157	cout<<"Config.C: Creating Run Loader ..."<<endl;
	158	rl = AliRunLoader::Open("galice.root",
f5a857b2	159	AliConfig::GetDefaultEventFolderName(),
de78d77f	160	"recreate");
	161	if (rl == 0x0)
	162	{
	163	gAlice->Fatal("Config.C","Can not instatiate the Run Loader");
	164	return;
	165	}
	166	rl->SetCompressionLevel(2);
7e35f2ed	167	rl->SetNumberOfEventsPerFile(1000);
de78d77f	168	gAlice->SetRunLoader(rl);
458041af	169	// gAlice->SetGeometryFromFile("geometry.root");
458041af	170	// gAlice->SetGeometryFromCDB();
de78d77f	171
dcec124d	172	// Set the trigger configuration
28da60d3	173	AliSimulation::Instance()->SetTriggerConfig(TrigConfName[trig]);
dcec124d	174	cout<<"Trigger configuration is set to "<<TrigConfName[trig]<<endl;
dcec124d	175
6ccd5a72	176	//
	177	//=======================================================================
	178	// *********** STEERING parameters FOR ALICE SIMULATION ************
	179	// --- Specify event type to be tracked through the ALICE setup
	180	// --- All positions are in cm, angles in degrees, and P and E in GeV
6a3df6c9	181
	182
	183	gMC->SetProcess("DCAY",1);
	184	gMC->SetProcess("PAIR",1);
	185	gMC->SetProcess("COMP",1);
	186	gMC->SetProcess("PHOT",1);
	187	gMC->SetProcess("PFIS",0);
	188	gMC->SetProcess("DRAY",0);
	189	gMC->SetProcess("ANNI",1);
	190	gMC->SetProcess("BREM",1);
	191	gMC->SetProcess("MUNU",1);
	192	gMC->SetProcess("CKOV",1);
	193	gMC->SetProcess("HADR",1);
	194	gMC->SetProcess("LOSS",2);
	195	gMC->SetProcess("MULS",1);
	196	gMC->SetProcess("RAYL",1);
	197
	198	Float_t cut = 1.e-3; // 1MeV cut by default
	199	Float_t tofmax = 1.e10;
	200
	201	gMC->SetCut("CUTGAM", cut);
	202	gMC->SetCut("CUTELE", cut);
	203	gMC->SetCut("CUTNEU", cut);
	204	gMC->SetCut("CUTHAD", cut);
	205	gMC->SetCut("CUTMUO", cut);
	206	gMC->SetCut("BCUTE", cut);
	207	gMC->SetCut("BCUTM", cut);
	208	gMC->SetCut("DCUTE", cut);
	209	gMC->SetCut("DCUTM", cut);
	210	gMC->SetCut("PPCUTM", cut);
	211	gMC->SetCut("TOFMAX", tofmax);
	212
	213
	214
	215
6ccd5a72	216	// Set External decayer //
	217	//======================//
	218	TVirtualMCDecayer* decayer = new AliDecayerPythia();
	219	// DECAYS
	220	//
	221	switch(decHvFl) {
	222	case kNature:
	223	decayer->SetForceDecay(kAll);
	224	break;
	225	case kHadr:
	226	decayer->SetForceDecay(kHadronicD);
	227	break;
	228	case kSemiEl:
	229	decayer->SetForceDecay(kSemiElectronic);
	230	break;
	231	case kSemiMu:
	232	decayer->SetForceDecay(kSemiMuonic);
	233	break;
	234	}
	235	decayer->Init();
	236	gMC->SetExternalDecayer(decayer);
	237
	238	//=========================//
	239	// Generator Configuration //
	240	//=========================//
	241	AliGenPythia *pythia = PythiaHVQ(procHvFl);
	242	// FeedDown option
	243	pythia->SetFeedDownHigherFamily(kFALSE);
	244	// Stack filling option
	245	if(!stars) pythia->SetStackFillOpt(AliGenPythia::kParentSelection);
	246	// Set Count mode
	247	if(nEvts>0) pythia->SetCountMode(AliGenPythia::kCountParents);
	248	// DECAYS
	249	//
	250	switch(decHvFl) {
	251	case kNature:
	252	pythia->SetForceDecay(kAll);
	253	break;
	254	case kHadr:
	255	pythia->SetForceDecay(kHadronicD);
	256	break;
	257	case kSemiEl:
	258	pythia->SetForceDecay(kSemiElectronic);
	259	break;
	260	case kSemiMu:
	261	pythia->SetForceDecay(kSemiMuonic);
	262	break;
	263	}
	264	// GEOM & KINE CUTS
	265	//
	266	pythia->SetMomentumRange(0,99999999);
0ac28de4	267	pythia->SetPhiRange(0., 360.);
6ccd5a72	268	pythia->SetThetaRange(0,180);
	269	switch(ycut) {
	270	case kFull:
	271	pythia->SetYRange(-999,999);
	272	break;
	273	case kBarrel:
	274	pythia->SetYRange(-2,2);
	275	break;
	276	case kMuonArm:
	277	pythia->SetYRange(1,6);
	278	break;
	279	}
	280
	281	// PRIMARY VERTEX
	282	//
	283	pythia->SetOrigin(0, 0, 0); // vertex position
	284	pythia->SetSigma(0, 0, 5.3); // Sigma in (X,Y,Z) (cm) on IP position
	285	pythia->SetCutVertexZ(1.); // Truncate at 1 sigma
	286	pythia->SetVertexSmear(kPerEvent);
	287
	288
	289
	290	pythia->SetTrackingFlag(0);
	291	// Specify GEANT tracking limits (Rmax, Zmax)
e196c4dc	292	//gAlice->GetMCApp()->TrackingLimits(90.,1.0e10);
6ccd5a72	293
	294
	295	pythia->Init();
	296
	297	// FIELD
	298	//
f7a1cc68	299	if (mag == AliMagF::k2kG) {
6ccd5a72	300	comment = comment.Append(" \| L3 field 0.2 T");
f7a1cc68	301	} else if (mag == AliMagF::k5kG) {
6ccd5a72	302	comment = comment.Append(" \| L3 field 0.5 T");
6ccd5a72	303	}
de78d77f	304	printf("\n \n Comment: %s \n \n", comment.Data());
6ccd5a72	305
4642ac4b	306	AliMagF* field = new AliMagF("Maps","Maps", 1., 1., mag);
e196c4dc	307	field->SetL3ConstField(0); //Using const. field in the barrel
f7a1cc68	308	TGeoGlobalMagField::Instance()->SetField(field);
6ccd5a72	309
f7a1cc68	310	rl->CdGAFile();
6ccd5a72	311
	312	// By default all ALICE is switched off
	313	Int_t iABSO=0;
b384f8a4	314	Int_t iACORDE=0;
6ccd5a72	315	Int_t iDIPO=0;
	316	Int_t iFMD=0;
	317	Int_t iFRAME=0;
	318	Int_t iHALL=0;
	319	Int_t iITS=0;
	320	Int_t iMAG=0;
	321	Int_t iMUON=0;
	322	Int_t iPHOS=0;
	323	Int_t iPIPE=0;
	324	Int_t iPMD=0;
de1857bb	325	Int_t iHMPID=0;
6ccd5a72	326	Int_t iSHIL=0;
ababa197	327	Int_t iT0=0;
6ccd5a72	328	Int_t iTOF=0;
	329	Int_t iTPC=0;
	330	Int_t iTRD=0;
	331	Int_t iZDC=0;
	332
	333	//=================== Alice BODY parameters =============================
	334	AliBODY *BODY = new AliBODY("BODY","Alice envelop");
	335
	336	if(iMAG) {
	337	//=================== MAG parameters ============================
	338	// --- Start with Magnet since detector layouts may be depending ---
	339	// --- on the selected Magnet dimensions ---
	340	AliMAG *MAG = new AliMAG("MAG","Magnet");
	341	}
	342
	343
	344	if(iABSO) {
	345	//=================== ABSO parameters ============================
96fb3353	346	AliABSO *ABSO = new AliABSOv3("ABSO","Muon Absorber");
6ccd5a72	347	}
	348
	349	if(iDIPO) {
	350	//=================== DIPO parameters ============================
	351
96fb3353	352	AliDIPO *DIPO = new AliDIPOv3("DIPO","Dipole version 3");
6ccd5a72	353	}
	354
	355	if(iHALL) {
	356	//=================== HALL parameters ============================
	357
96fb3353	358	AliHALL *HALL = new AliHALLv3("HALL","Alice Hall");
6ccd5a72	359	}
	360
	361
	362	if(iFRAME) {
	363	//=================== FRAME parameters ============================
	364
	365	AliFRAME *FRAME = new AliFRAMEv2("FRAME","Space Frame");
bab433b3	366	FRAME->SetHoles(1);
6ccd5a72	367	}
	368
	369	if(iSHIL) {
	370	//=================== SHIL parameters ============================
	371
96fb3353	372	AliSHIL *SHIL = new AliSHILv3("SHIL","Shielding");
6ccd5a72	373	}
	374
	375
	376	if(iPIPE) {
	377	//=================== PIPE parameters ============================
	378
96fb3353	379	AliPIPE *PIPE = new AliPIPEv3("PIPE","Beam Pipe");
6ccd5a72	380	}
	381
	382
	383	if(iITS) {
e89fa83b	384	//=================== ITS parameters ============================
6ccd5a72	385
99887042	386	AliITS *ITS = new AliITSv11Hybrid("ITS","ITS v11Hybrid");
6ccd5a72	387	}
	388
	389
	390	if(iTPC) {
26ca4b4a	391	//============================ TPC parameters ===================
6ccd5a72	392	AliTPC *TPC = new AliTPCv2("TPC","Default");
6ccd5a72	393	}
	394
	395
	396	if(iTOF) {
	397	//=================== TOF parameters ============================
0d7ef405	398	AliTOF *TOF = new AliTOFv6T0("TOF","normal TOF");
6ccd5a72	399	}
6ccd5a72	400
de1857bb	401	if(iHMPID) {
de1857bb	402	//=================== HMPID parameters ===========================
99887042	403	AliHMPID *HMPID = new AliHMPIDv3("HMPID","normal HMPID");
6ccd5a72	404
	405	}
	406
	407
	408	if(iZDC) {
	409	//=================== ZDC parameters ============================
	410
9d07603a	411	AliZDC *ZDC = new AliZDCv3("ZDC","normal ZDC");
6ccd5a72	412	}
6ccd5a72	413
b384f8a4	414	if(iACORDE) {
b384f8a4	415	//=================== ACORDE parameters ============================
6ccd5a72	416
b384f8a4	417	AliACORDE *ACORDE = new AliACORDEv1("ACORDE","normal ACORDE");
6ccd5a72	418	}
	419
	420	if(iTRD) {
	421	//=================== TRD parameters ============================
	422
	423	AliTRD *TRD = new AliTRDv1("TRD","TRD slow simulator");
6ccd5a72	424	}
	425
	426	if(iFMD) {
	427	//=================== FMD parameters ============================
	428
	429	AliFMD *FMD = new AliFMDv0("FMD","normal FMD");
	430	}
	431
	432	if(iMUON) {
	433	//=================== MUON parameters ===========================
	434	AliMUON *MUON = new AliMUONv1("MUON","default");
	435	}
	436
	437	//=================== PHOS parameters ===========================
	438
	439	if(iPHOS) {
	440	AliPHOS *PHOS = new AliPHOSv1("PHOS","GPS2");
	441	}
	442
	443
b384f8a4	444	//=================== ACORDE parameters ===========================
6ccd5a72	445
b384f8a4	446	if(iACORDE) {
b384f8a4	447	AliACORDE *ACORDE = new AliACORDEv1("ACORDE","Normal ACORDEGPS2");
6ccd5a72	448	}
	449
	450
	451	if(iPMD) {
	452	//=================== PMD parameters ============================
	453
	454	AliPMD *PMD = new AliPMDv1("PMD","normal PMD");
	455	PMD->SetPAR(1., 1., 0.8, 0.02);
	456	PMD->SetIN(6., 18., -580., 27., 27.);
	457	PMD->SetGEO(0.0, 0.2, 4.);
	458	PMD->SetPadSize(0.8, 1.0, 1.0, 1.5);
	459
	460	}
	461
ababa197	462	if(iT0) {
	463	//=================== T0 parameters ============================
	464	AliT0 *T0 = new AliT0v1("T0","T0 Detector");
6ccd5a72	465	}
	466
	467
	468	}
	469	//
	470	// PYTHIA
	471	//
	472	AliGenPythia *PythiaHVQ(ProcessHvFl_t proc) {
	473
de78d77f	474	AliGenPythia * gener = 0x0;
de78d77f	475
6ccd5a72	476	switch(proc) {
	477	case kCharmPbPb5500:
	478	comment = comment.Append(" Charm in Pb-Pb at 5.5 TeV");
de78d77f	479	gener = new AliGenPythia(nEvts);
6ccd5a72	480	gener->SetProcess(kPyCharmPbPbMNR);
	481	gener->SetStrucFunc(kCTEQ4L);
	482	gener->SetPtHard(2.1,-1.0);
	483	gener->SetEnergyCMS(5500.);
	484	gener->SetNuclei(208,208);
	485	break;
	486	case kCharmpPb8800:
	487	comment = comment.Append(" Charm in p-Pb at 8.8 TeV");
de78d77f	488	gener = new AliGenPythia(nEvts);
6ccd5a72	489	gener->SetProcess(kPyCharmpPbMNR);
	490	gener->SetStrucFunc(kCTEQ4L);
	491	gener->SetPtHard(2.1,-1.0);
	492	gener->SetEnergyCMS(8800.);
	493	gener->SetProjectile("P",1,1);
	494	gener->SetTarget("Pb",208,82);
	495	break;
	496	case kCharmpp14000:
	497	comment = comment.Append(" Charm in pp at 14 TeV");
de78d77f	498	gener = new AliGenPythia(nEvts);
6ccd5a72	499	gener->SetProcess(kPyCharmppMNR);
	500	gener->SetStrucFunc(kCTEQ4L);
	501	gener->SetPtHard(2.1,-1.0);
	502	gener->SetEnergyCMS(14000.);
	503	break;
09977179	504	case kCharmpp14000wmi:
	505	comment = comment.Append(" Charm in pp at 14 TeV with mult. interactions");
	506	gener = new AliGenPythia(-1);
	507	gener->SetProcess(kPyCharmppMNRwmi);
	508	gener->SetStrucFunc(kCTEQ5L);
	509	gener->SetPtHard(ptHardMin,ptHardMax);
	510	gener->SetEnergyCMS(14000.);
	511	break;
7e35f2ed	512	case kCharmSemiElpp14000wmi:
	513	comment = comment.Append(" Charm in pp at 14 TeV with mult. interactions");
	514	gener = new AliGenPythia(-1);
	515	gener->SetProcess(kPyCharmppMNRwmi);
	516	gener->SetStrucFunc(kCTEQ5L);
	517	gener->SetPtHard(ptHardMin,ptHardMax);
	518	gener->SetEnergyCMS(14000.);
	519	gener->SetCutOnChild(1);
	520	gener->SetPdgCodeParticleforAcceptanceCut(11);
	521	gener->SetChildYRange(-1,1);
	522	gener->SetChildPtRange(1,10000.);
	523	break;
6ccd5a72	524	case kD0PbPb5500:
6ccd5a72	525	comment = comment.Append(" D0 in Pb-Pb at 5.5 TeV");
de78d77f	526	gener = new AliGenPythia(nEvts);
6ccd5a72	527	gener->SetProcess(kPyD0PbPbMNR);
	528	gener->SetStrucFunc(kCTEQ4L);
	529	gener->SetPtHard(2.1,-1.0);
	530	gener->SetEnergyCMS(5500.);
	531	gener->SetNuclei(208,208);
	532	break;
	533	case kD0pPb8800:
	534	comment = comment.Append(" D0 in p-Pb at 8.8 TeV");
de78d77f	535	gener = new AliGenPythia(nEvts);
6ccd5a72	536	gener->SetProcess(kPyD0pPbMNR);
	537	gener->SetStrucFunc(kCTEQ4L);
	538	gener->SetPtHard(2.1,-1.0);
	539	gener->SetEnergyCMS(8800.);
	540	gener->SetProjectile("P",1,1);
	541	gener->SetTarget("Pb",208,82);
	542	break;
	543	case kD0pp14000:
	544	comment = comment.Append(" D0 in pp at 14 TeV");
de78d77f	545	gener = new AliGenPythia(nEvts);
6ccd5a72	546	gener->SetProcess(kPyD0ppMNR);
6ccd5a72	547	gener->SetStrucFunc(kCTEQ4L);
72fc043d	548	gener->SetPtHard(2.1,-1.0);
	549	gener->SetEnergyCMS(14000.);
	550	break;
	551	case kDPlusPbPb5500:
	552	comment = comment.Append(" DPlus in Pb-Pb at 5.5 TeV");
	553	gener = new AliGenPythia(nEvts);
	554	gener->SetProcess(kPyDPlusPbPbMNR);
	555	gener->SetStrucFunc(kCTEQ4L);
	556	gener->SetPtHard(2.1,-1.0);
	557	gener->SetEnergyCMS(5500.);
	558	gener->SetNuclei(208,208);
	559	break;
	560	case kDPluspPb8800:
	561	comment = comment.Append(" DPlus in p-Pb at 8.8 TeV");
	562	gener = new AliGenPythia(nEvts);
	563	gener->SetProcess(kPyDPluspPbMNR);
	564	gener->SetStrucFunc(kCTEQ4L);
	565	gener->SetPtHard(2.1,-1.0);
	566	gener->SetEnergyCMS(8800.);
	567	gener->SetProjectile("P",1,1);
	568	gener->SetTarget("Pb",208,82);
	569	break;
	570	case kDPluspp14000:
	571	comment = comment.Append(" DPlus in pp at 14 TeV");
	572	gener = new AliGenPythia(nEvts);
	573	gener->SetProcess(kPyDPlusppMNR);
	574	gener->SetStrucFunc(kCTEQ4L);
6ccd5a72	575	gener->SetPtHard(2.1,-1.0);
	576	gener->SetEnergyCMS(14000.);
	577	break;
84388d23	578	case kDPlusStrangePbPb5500:
	579	comment = comment.Append(" DPlusStrange in Pb-Pb at 5.5 TeV");
	580	gener = new AliGenPythia(nEvts);
	581	gener->SetProcess(kPyDPlusStrangePbPbMNR);
	582	gener->SetStrucFunc(kCTEQ4L);
	583	gener->SetPtHard(2.1,-1.0);
	584	gener->SetEnergyCMS(5500.);
	585	gener->SetNuclei(208,208);
	586	break;
	587	case kDPlusStrangepPb8800:
	588	comment = comment.Append(" DPlusStrange in p-Pb at 8.8 TeV");
	589	gener = new AliGenPythia(nEvts);
	590	gener->SetProcess(kPyDPlusStrangepPbMNR);
	591	gener->SetStrucFunc(kCTEQ4L);
	592	gener->SetPtHard(2.1,-1.0);
	593	gener->SetEnergyCMS(8800.);
	594	gener->SetProjectile("P",1,1);
	595	gener->SetTarget("Pb",208,82);
	596	break;
	597	case kDPlusStrangepp14000:
	598	comment = comment.Append(" DPlusStrange in pp at 14 TeV");
	599	gener = new AliGenPythia(nEvts);
	600	gener->SetProcess(kPyDPlusStrangeppMNR);
	601	gener->SetStrucFunc(kCTEQ4L);
	602	gener->SetPtHard(2.1,-1.0);
	603	gener->SetEnergyCMS(14000.);
	604	break;
6ccd5a72	605	case kBeautyPbPb5500:
6ccd5a72	606	comment = comment.Append(" Beauty in Pb-Pb at 5.5 TeV");
7005f584	607	gener = new AliGenPythia(nEvts);
6ccd5a72	608	gener->SetProcess(kPyBeautyPbPbMNR);
	609	gener->SetStrucFunc(kCTEQ4L);
	610	gener->SetPtHard(2.75,-1.0);
	611	gener->SetEnergyCMS(5500.);
	612	gener->SetNuclei(208,208);
	613	break;
	614	case kBeautypPb8800:
	615	comment = comment.Append(" Beauty in p-Pb at 8.8 TeV");
de78d77f	616	gener = new AliGenPythia(nEvts);
6ccd5a72	617	gener->SetProcess(kPyBeautypPbMNR);
	618	gener->SetStrucFunc(kCTEQ4L);
	619	gener->SetPtHard(2.75,-1.0);
	620	gener->SetEnergyCMS(8800.);
	621	gener->SetProjectile("P",1,1);
	622	gener->SetTarget("Pb",208,82);
	623	break;
	624	case kBeautypp14000:
	625	comment = comment.Append(" Beauty in pp at 14 TeV");
de78d77f	626	gener = new AliGenPythia(nEvts);
6ccd5a72	627	gener->SetProcess(kPyBeautyppMNR);
	628	gener->SetStrucFunc(kCTEQ4L);
	629	gener->SetPtHard(2.75,-1.0);
	630	gener->SetEnergyCMS(14000.);
	631	break;
09977179	632	case kBeautypp14000wmi:
	633	comment = comment.Append(" Beauty in pp at 14 TeV with mult. interactions");
	634	gener = new AliGenPythia(-1);
	635	gener->SetProcess(kPyBeautyppMNRwmi);
	636	gener->SetStrucFunc(kCTEQ5L);
	637	gener->SetPtHard(ptHardMin,ptHardMax);
	638	gener->SetEnergyCMS(14000.);
	639	break;
7e35f2ed	640	case kBeautySemiElpp14000wmi:
	641	comment = comment.Append(" Beauty in pp at 14 TeV with mult. interactions");
	642	gener = new AliGenPythia(-1);
	643	gener->SetProcess(kPyBeautyppMNRwmi);
	644	gener->SetStrucFunc(kCTEQ5L);
	645	gener->SetPtHard(ptHardMin,ptHardMax);
	646	gener->SetEnergyCMS(14000.);
	647	gener->SetCutOnChild(1);
	648	gener->SetPdgCodeParticleforAcceptanceCut(11);
	649	gener->SetChildYRange(-1,1);
	650	gener->SetChildPtRange(1,10000.);
	651	break;
6ccd5a72	652	}
	653
	654	return gener;
	655	}
	656
	657
	658
ee621efa	659	void LoadPythia()
	660	{
	661	// Load Pythia related libraries
	662	gSystem->Load("liblhapdf.so"); // Parton density functions
	663	gSystem->Load("libEGPythia6.so"); // TGenerator interface
	664	gSystem->Load("libpythia6.so"); // Pythia
	665	gSystem->Load("libAliPythia6.so"); // ALICE specific implementations
	666	}
	667
6ccd5a72	668
6ccd5a72	669