[u/mrichter/AliRoot.git] / MFT / Config.C

#if !defined(__CINT__) || defined(__MAKECINT__)
#include <Riostream.h>
#include <TRandom.h>
#include <TDatime.h>
#include <TSystem.h>
#include <TVirtualMC.h>
#include <TGeant3TGeo.h>
#include "AliRunLoader.h"
#include "AliRun.h"
#include "AliConfig.h"
#include "AliDecayerPythia.h"
#include "AliGenPythia.h"
#include "AliGenDPMjet.h"
#include "AliMagFCheb.h"
#include "AliBODY.h"
#include "AliMAG.h"
#include "AliABSOv3.h"
#include "AliDIPOv3.h"
#include "AliHALLv3.h"
#include "AliFRAMEv2.h"
#include "AliSHILv3.h"
#include "AliPIPEv3.h"
#include "AliPIPEv4.h"
#include "AliITSv11.h"
#include "AliTPCv2.h"
#include "AliTOFv6T0.h"
#include "AliHMPIDv3.h"
#include "AliZDCv3.h"
#include "AliTRDv1.h"
#include "AliTRDgeometry.h"
#include "AliFMDv1.h"
#include "AliMUONv1.h"
#include "AliPHOSv1.h"
#include "AliPHOSSimParam.h"
#include "AliPMDv1.h"
#include "AliT0v1.h"
#include "AliEMCALv2.h"
#include "AliACORDEv1.h"
#include "AliVZEROv7.h"
#include "AliMFT.h"
#endif

// Git push Test: R. Tieulent

enum PDCProc_t {kGenBox,
		kGenMuonLMR,
		kGenParamJpsi,
		kGenPionKaon,
		kGenCorrHF,
                kPythia6,
		kPythiaPerugia0, 
		kPythiaPerugia0Jpsi2mu, 
		kPythiaPerugia0BtoJpsi2mu, 
		kHijing, 
		kHijing2500,
		kHijing2500Cocktail};

const Char_t* pprRunName[] = {"kGenBox",
			      "kGenMuonLMR",
			      "kGenParamJpsi",
			      "kGenPionKaon",
			      "kGenCorrHF",
			      "kPythia6",
			      "kPythiaPerugia0", 
			      "kPythiaPerugia0Jpsi2mu", 
			      "kPythiaPerugia0BtoJpsi2mu", 
			      "kHijing", 
			      "kHijing2500", 
			      "kHijing2500Cocktail"};

enum Mag_t { kNoField, k5kG, kFieldMax };

const Char_t* pprField[] = { "kNoField", "k5kG", "kFieldMax" };

void LoadLibs();

// ----------------------- Generator, field, beam energy,... ------------------------------------------------------------
static PDCProc_t     proc     = kGenBox;
static PDCProc_t     signal   = kGenBox;    // only in case kHijing2500Cocktail is the proc
static Mag_t         mag      = k5kG;
static Float_t       energy   = 5500.; // energy in CMS
static Float_t       bMin     = 0.;
static Float_t       bMax =   = 5.; // 0-5 fm corresponds to around 0-10% (see https://twiki.cern.ch/twiki/bin/viewauth/ALICE/CentStudies#Tables_with_centrality_bins_for)
static Double_t      JpsiPol  = 0; // Jpsi polarisation
static Bool_t        JpsiHarderPt = kFALSE; // Jpsi harder pt spectrum (8.8 TeV)
// ----------------------------------------------------------------------------------------------------------------------

static TString comment;

//====================================================================================================================================================

void Config() {

  //  AliLog::SetClassDebugLevel("AliMFT", 1);

  LoadLibs();

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

  // Create the output file

  AliRunLoader* rl=0x0;

  printf("Config.C: Creating Run Loader ...");
  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);

  // ************* 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); 
  
  TVirtualMCDecayer *decayer = new AliDecayerPythia();
  decayer->SetForceDecay(kAll);
  decayer->Init();
  gMC->SetExternalDecayer(decayer);
  
  // Generator
  AliGenerator* gener = 0x0;
  if (proc == kPythia6)                        gener = MbPythia(); 
  else if (proc == kPythiaPerugia0)            gener = MbPythiaTunePerugia0();
  else if (proc == kHijing)                    gener = Hijing();	
  else if (proc == kHijing2500)                gener = Hijing2500();	
  else if (proc == kHijing2500Cocktail)        gener = Hijing2500Cocktail();
  else if (proc == kGenBox)                    gener = GenBox();
  else if (proc == kGenMuonLMR)                gener = GenMuonLMR();
  else if (proc == kGenParamJpsi)              gener = GenParamJpsi();
  else if (proc == kGenCorrHF)                 gener = GenCorrHF();
  else if (proc == kGenPionKaon)               gener = GenParamPionKaon();
  else if (proc == kPythiaPerugia0BtoJpsi2mu)  gener = MbPythiaTunePerugia0BtoJpsi2mu();

  // Size of the interaction diamond
  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();

  printf("\n \n Comment: %s \n \n", comment.Data());

  //  TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", -1., -1., AliMagF::k5kG, AliMagF::kBeamTypeAA, 2750.));
  TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", -1., -1., AliMagF::k5kG, AliMagF::kBeamTypepp, 7000.));

  rl->CdGAFile();
  
  // Detector Setup

  Int_t iABSO  = 1;
  Int_t iDIPO  = 1;
  Int_t iHALL  = 1;
  Int_t iMUON  = 1;
  Int_t iPIPE  = 1;
  Int_t iSHIL  = 1;
  Int_t iT0    = 0;
  Int_t iVZERO = 1;
  Int_t iMFT   = 1;
  Int_t iACORDE= 0;
  Int_t iEMCAL = 0;
  Int_t iFMD   = 0;
  Int_t iFRAME = 0;
  Int_t iITS   = 0;
  Int_t iMAG   = 0;
  Int_t iPHOS  = 0;
  Int_t iPMD   = 0;
  Int_t iHMPID = 0;
  Int_t iTOF   = 0;
  Int_t iTPC   = 0;
  Int_t iTRD   = 0;
  Int_t iZDC   = 0;
  

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

  if (iMAG)       AliMAG    *MAG    = new AliMAG("MAG", "Magnet");
  if (iABSO)      AliABSO   *ABSO   = new AliABSOv3("ABSO", "Muon Absorber");
  if (iDIPO)      AliDIPO   *DIPO   = new AliDIPOv3("DIPO", "Dipole version 3");
  if (iHALL)      AliHALL   *HALL   = new AliHALLv3("HALL", "Alice Hall");
  if (iSHIL)      AliSHIL   *SHIL   = new AliSHILv3("SHIL", "Shielding Version 3");
  if (iITS)       AliITS    *ITS    = new AliITSv11("ITS","ITS v11");
  if (iTPC)       AliTPC    *TPC    = new AliTPCv2("TPC", "Default");
  if (iTOF)       AliTOF    *TOF    = new AliTOFv6T0("TOF", "normal TOF");
  if (iHMPID)     AliHMPID  *HMPID  = new AliHMPIDv3("HMPID", "normal HMPID");
  if (iFMD)       AliFMD    *FMD    = new AliFMDv1("FMD", "normal FMD");
  if (iPHOS)      AliPHOS   *PHOS   = new AliPHOSv1("PHOS", "noCPV_Modules123");
  if (iPMD)       AliPMD    *PMD    = new AliPMDv1("PMD", "normal PMD");
  if (iT0)        AliT0     *T0     = new AliT0v1("T0", "T0 Detector");
  if (iEMCAL)     AliEMCAL  *EMCAL  = new AliEMCALv2("EMCAL", "EMCAL_FIRSTYEARV1");
  if (iACORDE)    AliACORDE *ACORDE = new AliACORDEv1("ACORDE", "normal ACORDE");
  if (iVZERO)     AliVZERO  *VZERO  = new AliVZEROv7("VZERO", "normal VZERO");
  if (iFRAME) {
    AliFRAMEv2 *FRAME = new AliFRAMEv2("FRAME", "Space Frame");
    FRAME->SetHoles(1);
  }
  if (iPIPE) {
    //    AliPIPE *PIPE = new AliPIPEv3("PIPE", "Beam Pipe");
    AliPIPE *PIPE = new AliPIPEv4("PIPE", "Beam Pipe");
  }
  if (iZDC) {
    AliZDC *ZDC = new AliZDCv3("ZDC", "normal ZDC");
    ZDC->SetSpectatorsTrack();	
    ZDC->SetLumiLength(0.);
  }
  if (iTRD) {
    AliTRD *TRD = new AliTRDv1("TRD", "TRD slow simulator");
  }
  if (iMUON) {
    AliMUON *MUON = new AliMUONv1("MUON", "default");
    MUON->SetTriggerEffCells(1); // not needed if raw masks 
    Char_t* digitstore="AliMUONDigitStoreV2S";    
    MUON->SetDigitStoreClassName(digitstore);
  }
  if (iMFT) {
    AliMFT *MFT = new AliMFT("MFT", "normal MFT");
  }

  TIter next(gAlice->Modules());
  AliModule *detector;
  printf("gAlice->Modules:\n");
  while((detector = (AliModule*)next())) printf("%s\n",detector->GetName());

}

//====================================================================================================================================================

AliGenerator* GenBox() {

  AliGenBox *gener = new AliGenBox(2);
  gener->SetMomentumRange(10.0, 10.1);
  gener->SetPhiRange(0., 360.);         
  gener->SetThetaRange(171.0,178.0);
  gener->SetPart(kMuonMinus);           // Muons
  gener->SetOrigin(0., 0., 0.);         // vertex position
  
  return gener;
  
}

//====================================================================================================================================================

AliGenerator* GenMuonLMR() {
  
  AliGenMUONLMR *gener = new AliGenMUONLMR();
  gener->SetMomentumRange(0,999);
  gener->SetPtRange(0,100.);
  gener->SetYRange(-4.5, -2.0);
  gener->SetChildThetaRange(171.0,178.0);
  gener->SetChildMomentumRange(4.0, 999.);
  gener->SetOrigin(0.0, 0.0, 0.0);             // vertex position
  gener->SetSigma(0.0, 0.0, 0.0);              // vertex position smearing
  enum {kEta2Body, kEtaDalitz, kRho2Body, kOmega2Body, kOmegaDalitz, kPhi2Body, kEtaPrimeDalitz, kPionLMR, kKaonLMR}; 
  gener->GenerateSingleProcess(kOmega2Body, 10);
  gener->SetCutOnChild(1);

  return gener;

}

//====================================================================================================================================================

AliGenerator* GenParamJpsi() {

  AliGenParam *gener = new AliGenParam(5, AliGenMUONlib::kJpsi);
  gener->SetMomentumRange(0,999);
  gener->SetPtRange(0,100.);
  gener->SetYRange(-4.0, -2.5);
  gener->SetPhiRange(0., 360.);
  gener->SetChildThetaRange(171.0,177.0);
  gener->SetChildMomentumRange(5.0, 999.);
  gener->SetOrigin(0.0, 0.0, 0.0);          // vertex position
  gener->SetSigma(0.0, 0.0, 0.0);           // Sigma in (X,Y,Z) (cm) on IP position
  gener->SetForceDecay(kDiMuon);
  gener->SetTrackingFlag(1);
  gener->SetCutOnChild(1);

  return gener;

}

//====================================================================================================================================================

AliGenerator* GenParamPionKaon() {
  
  AliGenParamPionsKaons *gener = new AliGenParamPionsKaons(100,"$ALICE_ROOT/MFT/PionKaonKinematics.root");
  gener->SetPtRange(0, 5.);
  gener->SetPhiRange(0., 360.);
  gener->SetYRange(-10., 0.);
  gener->SetOrigin(0.0, 0.0, 0.0);  // vertex position
  gener->SetSigma(0.0, 0.0, 0.0);   // vertex position smearing
  //  gener->SetCutOnChild(1);

  return gener;

}

//====================================================================================================================================================

AliGenerator* GenCorrHF() {
  
  AliGenCorrHF *gener = new AliGenCorrHF(1, 4, 6);  // for charm, 1 pair per event
  // AliGenCorrHF *gener = new AliGenCorrHF(1, 5, 6);  // for beauty, 1 pair per event
  
  gener->SetMomentumRange(0,9999);
  gener->SetCutOnChild(1);          // 1/0 means cuts on children enable/disable
  gener->SetChildThetaRange(171.0,178.0);
  gener->SetChildMomentumRange(4.0, 999.);
  gener->SetOrigin(0,0,0);          //vertex position    
  gener->SetForceDecay(kSemiMuonic);
  gener->SetTrackingFlag(1);
  gener->Init();
  
  return gener;
  
}

//====================================================================================================================================================

AliGenerator* MbPythia() {
  
  comment = comment.Append(" pp: Pythia low-pt");
  
  //    Pythia
  AliGenPythia* pythia = new AliGenPythia(-1); 
  pythia->SetMomentumRange(0, 999999.);
  //  pythia->SetThetaRange(0., 180.);
  //  pythia->SetChildYRange(-12.,0.);
  //  pythia->SetPtRange(0,1000.);
  //  pythia->SetCutOnChild(1);
  pythia->SetProcess(kPyMb);
  pythia->SetEnergyCMS(energy);
  pythia->SetForceDecay(kSemiMuonic);
  
  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* MbPythiaTunePerugia0Jpsi2mu() {
  
  comment = comment.Append("Jpsi forced to dimuons");
  AliGenParam *jpsi=0x0;
  if (JpsiHarderPt) jpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp 8.8", "Jpsi");  // 8.8 TeV
  else jpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp 7", "Jpsi");  // 7 TeV
  jpsi->SetPtRange(0.,999.);
  jpsi->SetYRange(-1.0, 1.0);
  jpsi->SetPhiRange(0.,360.);
  jpsi->SetOrigin(0,0,0);
  jpsi->SetForceDecay(kDiMuon);
  return jpsi;

}

//====================================================================================================================================================

AliGenerator* MbPythiaTunePerugia0BtoJpsi2mu() {

  comment = comment.Append(" pp: Pythia (Perugia0) BtoJpsi (1 bbbar per event, 1 b-hadron, 1 J/psi");
  
  //    Pythia
  AliGenPythia* pythia = new AliGenPythia(-1);
  pythia->SetMomentumRange(0, 999999.);
  pythia->SetYRange(-4.5, -2.0);
  pythia->SetPtRange(0,1000.);
  pythia->SetProcess(kPyBeautyppMNRwmi);
  pythia->SetEnergyCMS(energy);
  //    Tune
  //    320     Perugia 0
  pythia->SetTune(320);
  pythia->UseNewMultipleInteractionsScenario();
  //
  //    decays
  pythia->SetCutOnChild(1);
  pythia->SetPdgCodeParticleforAcceptanceCut(443);
  pythia->SetChildYRange(-4.5, -2.0);
  pythia->SetChildPtRange(0,10000.);
  //
  //    decays
  pythia->SetForceDecay(kBJpsiDiMuon);
  
  return pythia;

}

//====================================================================================================================================================

AliGenerator* HijingParam() {

  AliGenHIJINGpara *gener = new AliGenHIJINGpara(2000);
  gener->SetMomentumRange(0,999);
  gener->SetPhiRange(0,360);
  gener->SetThetaRange(171, 179);
  gener->SetOrigin(0,0,0);        // vertex position
  gener->SetSigma(0,0,0);         // Sigma in (X,Y,Z) (cm) on IP position
  gGener = gener;
  
  return gener;

}

//====================================================================================================================================================

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->SetShadowing(1);
  // Don't track spectators
  gener->SetSpectators(0);
  // kinematic selection
  gener->SetSelectAll(0);
  return gener;
}

//====================================================================================================================================================

AliGenerator* Hijing2500() {
  
  AliGenHijing *gener = (AliGenHijing*) Hijing();
  gener->SetJetQuenching(0);	
  gener->SetPtHardMin (4.5);
  return gener;
  
}

//====================================================================================================================================================

AliGenerator* Hijing2500Cocktail() {
  
  comment = comment.Append(" PbPb: Hjing2500 at 5.5 + muon signals");

  AliGenCocktail *cocktail = new AliGenCocktail();
  cocktail->SetProjectile("A", 208, 82);
  cocktail->SetTarget    ("A", 208, 82);
  cocktail->SetEnergyCMS(energy);

  // 1 Hijing event: provides underlying event and collision geometry 
  AliGenHijing *hijing = Hijing2500();
  cocktail->AddGenerator(hijing,"hijing",1);

  // generator for the custom signal
  AliGenerator* signalGen = 0x0;      
  switch (signal) {
  case 0:
    signalGen = GenBox();
    break;
  case 1:
    signalGen = GenMuonLMR(); 
    break;
  case 2:
    signalGen = GenCorrHF();
    break;
  default:
    signalGen = GenBox();
    break;
  }
  cocktail->AddGenerator(signalGen, "signal", 1);

  cocktail->SetTrackingFlag(1);

  return cocktail;

}

//====================================================================================================================================================

AliGenerator* QEDGeneratorPbPb() {
  
  AliGenEpEmv1 *generQED = new AliGenEpEmv1();
  generQED->SetEnergyCMS(energy);
  generQED->SetProjectile("A", 208, 82);
  generQED->SetTarget ("A", 208, 82);
  generQED->SetYRange(-8.5, 0.);
  generQED->SetPtRange(0.4e-3, 1.0); // Set pt limits (GeV) for e+-
  generQED->SetDebug(0);            // Set debug level (0 = silent)

}

//====================================================================================================================================================

void LoadLibs() {

#if defined(__CINT__)

  gSystem->Load("liblhapdf");      // Parton density functions
  gSystem->Load("libEGPythia6");   // TGenerator interface
  if (proc == kPythia6) {
    gSystem->Load("libpythia6");        // Pythia 6.2
    gSystem->Load("libAliPythia6");     // ALICE specific implementations
  } 
  else {
    gSystem->Load("libpythia6.4.21");   // Pythia 6.4
    gSystem->Load("libAliPythia6");     // ALICE specific implementations	
  }
  
  if (proc == kHijing || proc == kHijing2500 || proc == kHijing2500Cocktail) {
    gSystem->Load("libhijing");	
    gSystem->Load("libTHijing");
  } 
  
  gSystem->Load("libgeant321");
  
#endif

}

//====================================================================================================================================================
Commit	Line	Data
c95698be	1	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
	2	#include <Riostream.h>
	3	#include <TRandom.h>
	4	#include <TDatime.h>
	5	#include <TSystem.h>
	6	#include <TVirtualMC.h>
	7	#include <TGeant3TGeo.h>
	8	#include "AliRunLoader.h"
	9	#include "AliRun.h"
	10	#include "AliConfig.h"
	11	#include "AliDecayerPythia.h"
	12	#include "AliGenPythia.h"
	13	#include "AliGenDPMjet.h"
	14	#include "AliMagFCheb.h"
	15	#include "AliBODY.h"
	16	#include "AliMAG.h"
	17	#include "AliABSOv3.h"
	18	#include "AliDIPOv3.h"
	19	#include "AliHALLv3.h"
	20	#include "AliFRAMEv2.h"
	21	#include "AliSHILv3.h"
	22	#include "AliPIPEv3.h"
	23	#include "AliPIPEv4.h"
	24	#include "AliITSv11.h"
	25	#include "AliTPCv2.h"
	26	#include "AliTOFv6T0.h"
	27	#include "AliHMPIDv3.h"
	28	#include "AliZDCv3.h"
	29	#include "AliTRDv1.h"
	30	#include "AliTRDgeometry.h"
	31	#include "AliFMDv1.h"
	32	#include "AliMUONv1.h"
	33	#include "AliPHOSv1.h"
	34	#include "AliPHOSSimParam.h"
	35	#include "AliPMDv1.h"
	36	#include "AliT0v1.h"
	37	#include "AliEMCALv2.h"
	38	#include "AliACORDEv1.h"
	39	#include "AliVZEROv7.h"
	40	#include "AliMFT.h"
	41	#endif
	42
98acd653	43	// Git push Test: R. Tieulent
98acd653	44
c95698be	45	enum PDCProc_t {kGenBox,
c95698be	46	kGenMuonLMR,
b5ab1ac4	47	kGenParamJpsi,
6e8e4c4a	48	kGenPionKaon,
c95698be	49	kGenCorrHF,
	50	kPythia6,
	51	kPythiaPerugia0,
a2b7dc2a	52	kPythiaPerugia0Jpsi2mu,
a2b7dc2a	53	kPythiaPerugia0BtoJpsi2mu,
c95698be	54	kHijing,
	55	kHijing2500,
	56	kHijing2500Cocktail};
	57
	58	const Char_t* pprRunName[] = {"kGenBox",
	59	"kGenMuonLMR",
b5ab1ac4	60	"kGenParamJpsi",
6e8e4c4a	61	"kGenPionKaon",
c95698be	62	"kGenCorrHF",
	63	"kPythia6",
	64	"kPythiaPerugia0",
a2b7dc2a	65	"kPythiaPerugia0Jpsi2mu",
a2b7dc2a	66	"kPythiaPerugia0BtoJpsi2mu",
c95698be	67	"kHijing",
	68	"kHijing2500",
	69	"kHijing2500Cocktail"};
	70
	71	enum Mag_t { kNoField, k5kG, kFieldMax };
	72
	73	const Char_t* pprField[] = { "kNoField", "k5kG", "kFieldMax" };
	74
	75	void LoadLibs();
	76
	77	// ----------------------- Generator, field, beam energy,... ------------------------------------------------------------
3a11a1cf	78	static PDCProc_t proc = kGenBox;
c95698be	79	static PDCProc_t signal = kGenBox; // only in case kHijing2500Cocktail is the proc
c95698be	80	static Mag_t mag = k5kG;
a2b7dc2a	81	static Float_t energy = 5500.; // energy in CMS
c95698be	82	static Float_t bMin = 0.;
	83	static Float_t bMax = = 5.; // 0-5 fm corresponds to around 0-10% (see https://twiki.cern.ch/twiki/bin/viewauth/ALICE/CentStudies#Tables_with_centrality_bins_for)
	84	static Double_t JpsiPol = 0; // Jpsi polarisation
	85	static Bool_t JpsiHarderPt = kFALSE; // Jpsi harder pt spectrum (8.8 TeV)
	86	// ----------------------------------------------------------------------------------------------------------------------
	87
	88	static TString comment;
	89
	90	//====================================================================================================================================================
	91
	92	void Config() {
	93
	94	// AliLog::SetClassDebugLevel("AliMFT", 1);
	95
	96	LoadLibs();
	97
	98	new TGeant3TGeo("C++ Interface to Geant3");
	99
	100	// Create the output file
	101
	102	AliRunLoader* rl=0x0;
	103
cabfe25b	104	printf("Config.C: Creating Run Loader ...");
c95698be	105	rl = AliRunLoader::Open("galice.root", AliConfig::GetDefaultEventFolderName(), "recreate");
	106	if (rl == 0x0) {
	107	gAlice->Fatal("Config.C","Can not instatiate the Run Loader");
	108	return;
	109	}
	110	rl->SetCompressionLevel(2);
	111	rl->SetNumberOfEventsPerFile(1000);
	112	gAlice->SetRunLoader(rl);
	113
	114	// *********** STEERING parameters FOR ALICE SIMULATION ************
	115	// --- Specify event type to be tracked through the ALICE setup
	116	// --- All positions are in cm, angles in degrees, and P and E in GeV
	117
	118	gMC->SetProcess("DCAY",1);
	119	gMC->SetProcess("PAIR",1);
	120	gMC->SetProcess("COMP",1);
	121	gMC->SetProcess("PHOT",1);
	122	gMC->SetProcess("PFIS",0);
	123	gMC->SetProcess("DRAY",0);
	124	gMC->SetProcess("ANNI",1);
	125	gMC->SetProcess("BREM",1);
	126	gMC->SetProcess("MUNU",1);
	127	gMC->SetProcess("CKOV",1);
	128	gMC->SetProcess("HADR",1);
	129	gMC->SetProcess("LOSS",2);
	130	gMC->SetProcess("MULS",1);
	131	gMC->SetProcess("RAYL",1);
	132
	133	Float_t cut = 1.e-3; // 1MeV cut by default
	134	Float_t tofmax = 1.e10;
	135
	136	gMC->SetCut("CUTGAM", cut);
	137	gMC->SetCut("CUTELE", cut);
	138	gMC->SetCut("CUTNEU", cut);
	139	gMC->SetCut("CUTHAD", cut);
	140	gMC->SetCut("CUTMUO", cut);
	141	gMC->SetCut("BCUTE", cut);
	142	gMC->SetCut("BCUTM", cut);
	143	gMC->SetCut("DCUTE", cut);
	144	gMC->SetCut("DCUTM", cut);
	145	gMC->SetCut("PPCUTM", cut);
	146	gMC->SetCut("TOFMAX", tofmax);
	147
	148	TVirtualMCDecayer *decayer = new AliDecayerPythia();
	149	decayer->SetForceDecay(kAll);
	150	decayer->Init();
	151	gMC->SetExternalDecayer(decayer);
	152
	153	// Generator
	154	AliGenerator* gener = 0x0;
a2b7dc2a	155	if (proc == kPythia6) gener = MbPythia();
	156	else if (proc == kPythiaPerugia0) gener = MbPythiaTunePerugia0();
	157	else if (proc == kHijing) gener = Hijing();
	158	else if (proc == kHijing2500) gener = Hijing2500();
	159	else if (proc == kHijing2500Cocktail) gener = Hijing2500Cocktail();
	160	else if (proc == kGenBox) gener = GenBox();
	161	else if (proc == kGenMuonLMR) gener = GenMuonLMR();
	162	else if (proc == kGenParamJpsi) gener = GenParamJpsi();
	163	else if (proc == kGenCorrHF) gener = GenCorrHF();
	164	else if (proc == kGenPionKaon) gener = GenParamPionKaon();
	165	else if (proc == kPythiaPerugia0BtoJpsi2mu) gener = MbPythiaTunePerugia0BtoJpsi2mu();
c95698be	166
	167	// Size of the interaction diamond
	168	Float_t sigmaz = 5.4 / TMath::Sqrt(2.); // [cm]
	169	Float_t betast = 3.5; // beta* [m]
	170	Float_t eps = 3.75e-6; // emittance [m]
	171	Float_t gamma = energy / 2.0 / 0.938272; // relativistic gamma [1]
	172	Float_t sigmaxy = TMath::Sqrt(eps * betast / gamma) / TMath::Sqrt(2.) * 100.; // [cm]
	173
	174	printf("\n \n Diamond size x-y: %10.3e z: %10.3e\n \n", sigmaxy, sigmaz);
	175
	176	// gener->SetSigma(sigmaxy, sigmaxy, sigmaz); // Sigma in (X,Y,Z) (cm) on IP position
	177	// gener->SetVertexSmear(kPerEvent);
	178	gener->Init();
	179
	180	printf("\n \n Comment: %s \n \n", comment.Data());
	181
	182	// TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", -1., -1., AliMagF::k5kG, AliMagF::kBeamTypeAA, 2750.));
	183	TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", -1., -1., AliMagF::k5kG, AliMagF::kBeamTypepp, 7000.));
	184
	185	rl->CdGAFile();
	186
	187	// Detector Setup
	188
	189	Int_t iABSO = 1;
	190	Int_t iDIPO = 1;
	191	Int_t iHALL = 1;
	192	Int_t iMUON = 1;
	193	Int_t iPIPE = 1;
	194	Int_t iSHIL = 1;
	195	Int_t iT0 = 0;
	196	Int_t iVZERO = 1;
	197	Int_t iMFT = 1;
	198	Int_t iACORDE= 0;
	199	Int_t iEMCAL = 0;
	200	Int_t iFMD = 0;
	201	Int_t iFRAME = 0;
	202	Int_t iITS = 0;
	203	Int_t iMAG = 0;
	204	Int_t iPHOS = 0;
	205	Int_t iPMD = 0;
	206	Int_t iHMPID = 0;
	207	Int_t iTOF = 0;
	208	Int_t iTPC = 0;
	209	Int_t iTRD = 0;
	210	Int_t iZDC = 0;
	211
	212
	213	AliBODY *BODY = new AliBODY("BODY", "Alice envelop");
	214
	215	if (iMAG) AliMAG *MAG = new AliMAG("MAG", "Magnet");
	216	if (iABSO) AliABSO *ABSO = new AliABSOv3("ABSO", "Muon Absorber");
	217	if (iDIPO) AliDIPO *DIPO = new AliDIPOv3("DIPO", "Dipole version 3");
	218	if (iHALL) AliHALL *HALL = new AliHALLv3("HALL", "Alice Hall");
	219	if (iSHIL) AliSHIL *SHIL = new AliSHILv3("SHIL", "Shielding Version 3");
	220	if (iITS) AliITS *ITS = new AliITSv11("ITS","ITS v11");
	221	if (iTPC) AliTPC *TPC = new AliTPCv2("TPC", "Default");
	222	if (iTOF) AliTOF *TOF = new AliTOFv6T0("TOF", "normal TOF");
	223	if (iHMPID) AliHMPID *HMPID = new AliHMPIDv3("HMPID", "normal HMPID");
	224	if (iFMD) AliFMD *FMD = new AliFMDv1("FMD", "normal FMD");
	225	if (iPHOS) AliPHOS *PHOS = new AliPHOSv1("PHOS", "noCPV_Modules123");
	226	if (iPMD) AliPMD *PMD = new AliPMDv1("PMD", "normal PMD");
	227	if (iT0) AliT0 *T0 = new AliT0v1("T0", "T0 Detector");
	228	if (iEMCAL) AliEMCAL *EMCAL = new AliEMCALv2("EMCAL", "EMCAL_FIRSTYEARV1");
	229	if (iACORDE) AliACORDE *ACORDE = new AliACORDEv1("ACORDE", "normal ACORDE");
230	if (iVZERO) AliVZERO *VZERO = new AliVZEROv7("VZERO", "normal VZERO");
231	if (iFRAME) {
232	AliFRAMEv2 *FRAME = new AliFRAMEv2("FRAME", "Space Frame");
233	FRAME->SetHoles(1);
234	}
235	if (iPIPE) {
236	// AliPIPE *PIPE = new AliPIPEv3("PIPE", "Beam Pipe");
11ae1092	237	AliPIPE *PIPE = new AliPIPEv4("PIPE", "Beam Pipe");
c95698be	238	}
	239	if (iZDC) {
	240	AliZDC *ZDC = new AliZDCv3("ZDC", "normal ZDC");
	241	ZDC->SetSpectatorsTrack();
	242	ZDC->SetLumiLength(0.);
	243	}
	244	if (iTRD) {
	245	AliTRD *TRD = new AliTRDv1("TRD", "TRD slow simulator");
	246	}
	247	if (iMUON) {
	248	AliMUON *MUON = new AliMUONv1("MUON", "default");
	249	MUON->SetTriggerEffCells(1); // not needed if raw masks
	250	Char_t* digitstore="AliMUONDigitStoreV2S";
	251	MUON->SetDigitStoreClassName(digitstore);
	252	}
	253	if (iMFT) {
	254	AliMFT *MFT = new AliMFT("MFT", "normal MFT");
c95698be	255	}
	256
	257	TIter next(gAlice->Modules());
	258	AliModule *detector;
	259	printf("gAlice->Modules:\n");
	260	while((detector = (AliModule*)next())) printf("%s\n",detector->GetName());
	261
	262	}
	263
	264	//====================================================================================================================================================
	265
	266	AliGenerator* GenBox() {
	267
	268	AliGenBox *gener = new AliGenBox(2);
	269	gener->SetMomentumRange(10.0, 10.1);
	270	gener->SetPhiRange(0., 360.);
	271	gener->SetThetaRange(171.0,178.0);
	272	gener->SetPart(kMuonMinus); // Muons
	273	gener->SetOrigin(0., 0., 0.); // vertex position
	274
	275	return gener;
	276
	277	}
	278
	279	//====================================================================================================================================================
	280
	281	AliGenerator* GenMuonLMR() {
	282
	283	AliGenMUONLMR *gener = new AliGenMUONLMR();
	284	gener->SetMomentumRange(0,999);
	285	gener->SetPtRange(0,100.);
a2b7dc2a	286	gener->SetYRange(-4.5, -2.0);
	287	gener->SetChildThetaRange(171.0,178.0);
	288	gener->SetChildMomentumRange(4.0, 999.);
c95698be	289	gener->SetOrigin(0.0, 0.0, 0.0); // vertex position
c95698be	290	gener->SetSigma(0.0, 0.0, 0.0); // vertex position smearing
c95698be	291	enum {kEta2Body, kEtaDalitz, kRho2Body, kOmega2Body, kOmegaDalitz, kPhi2Body, kEtaPrimeDalitz, kPionLMR, kKaonLMR};
a2b7dc2a	292	gener->GenerateSingleProcess(kOmega2Body, 10);
b5ab1ac4	293	gener->SetCutOnChild(1);
	294
	295	return gener;
	296
	297	}
	298
	299	//====================================================================================================================================================
	300
	301	AliGenerator* GenParamJpsi() {
	302
	303	AliGenParam *gener = new AliGenParam(5, AliGenMUONlib::kJpsi);
	304	gener->SetMomentumRange(0,999);
	305	gener->SetPtRange(0,100.);
	306	gener->SetYRange(-4.0, -2.5);
	307	gener->SetPhiRange(0., 360.);
ada051d2	308	gener->SetChildThetaRange(171.0,177.0);
ada051d2	309	gener->SetChildMomentumRange(5.0, 999.);
b5ab1ac4	310	gener->SetOrigin(0.0, 0.0, 0.0); // vertex position
b5ab1ac4	311	gener->SetSigma(0.0, 0.0, 0.0); // Sigma in (X,Y,Z) (cm) on IP position
b5ab1ac4	312	gener->SetForceDecay(kDiMuon);
b5ab1ac4	313	gener->SetTrackingFlag(1);
c95698be	314	gener->SetCutOnChild(1);
	315
	316	return gener;
	317
	318	}
	319
	320	//====================================================================================================================================================
	321
6e8e4c4a	322	AliGenerator* GenParamPionKaon() {
	323
	324	AliGenParamPionsKaons *gener = new AliGenParamPionsKaons(100,"$ALICE_ROOT/MFT/PionKaonKinematics.root");
	325	gener->SetPtRange(0, 5.);
	326	gener->SetPhiRange(0., 360.);
	327	gener->SetYRange(-10., 0.);
	328	gener->SetOrigin(0.0, 0.0, 0.0); // vertex position
	329	gener->SetSigma(0.0, 0.0, 0.0); // vertex position smearing
	330	// gener->SetCutOnChild(1);
	331
	332	return gener;
	333
	334	}
	335
	336	//====================================================================================================================================================
	337
c95698be	338	AliGenerator* GenCorrHF() {
c95698be	339
a2b7dc2a	340	AliGenCorrHF *gener = new AliGenCorrHF(1, 4, 6); // for charm, 1 pair per event
a2b7dc2a	341	// AliGenCorrHF *gener = new AliGenCorrHF(1, 5, 6); // for beauty, 1 pair per event
c95698be	342
	343	gener->SetMomentumRange(0,9999);
	344	gener->SetCutOnChild(1); // 1/0 means cuts on children enable/disable
a2b7dc2a	345	gener->SetChildThetaRange(171.0,178.0);
a2b7dc2a	346	gener->SetChildMomentumRange(4.0, 999.);
c95698be	347	gener->SetOrigin(0,0,0); //vertex position
	348	gener->SetForceDecay(kSemiMuonic);
	349	gener->SetTrackingFlag(1);
	350	gener->Init();
	351
	352	return gener;
	353
	354	}
	355
	356	//====================================================================================================================================================
	357
	358	AliGenerator* MbPythia() {
	359
	360	comment = comment.Append(" pp: Pythia low-pt");
	361
	362	// Pythia
	363	AliGenPythia* pythia = new AliGenPythia(-1);
	364	pythia->SetMomentumRange(0, 999999.);
	365	// pythia->SetThetaRange(0., 180.);
	366	// pythia->SetChildYRange(-12.,0.);
	367	// pythia->SetPtRange(0,1000.);
	368	// pythia->SetCutOnChild(1);
	369	pythia->SetProcess(kPyMb);
	370	pythia->SetEnergyCMS(energy);
	371	pythia->SetForceDecay(kSemiMuonic);
	372
	373	return pythia;
	374	}
	375
	376	//====================================================================================================================================================
	377
	378	AliGenerator* MbPythiaTunePerugia0() {
	379
	380	comment = comment.Append(" pp: Pythia low-pt (Perugia0)");
	381
	382	// Pythia
	383	AliGenPythia* pythia = new AliGenPythia(-1);
	384	pythia->SetMomentumRange(0, 999999.);
	385	pythia->SetThetaRange(0., 180.);
	386	pythia->SetYRange(-12.,12.);
	387	pythia->SetPtRange(0,1000.);
	388	pythia->SetProcess(kPyMb);
	389	pythia->SetEnergyCMS(energy);
	390	// Tune
	391	// 320 Perugia 0
	392	pythia->SetTune(320);
	393	pythia->UseNewMultipleInteractionsScenario();
	394
	395	return pythia;
	396	}
	397
	398	//====================================================================================================================================================
	399
a2b7dc2a	400	AliGenerator* MbPythiaTunePerugia0Jpsi2mu() {
c95698be	401
a2b7dc2a	402	comment = comment.Append("Jpsi forced to dimuons");
c95698be	403	AliGenParam *jpsi=0x0;
a2b7dc2a	404	if (JpsiHarderPt) jpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp 8.8", "Jpsi"); // 8.8 TeV
c95698be	405	else jpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp 7", "Jpsi"); // 7 TeV
	406	jpsi->SetPtRange(0.,999.);
	407	jpsi->SetYRange(-1.0, 1.0);
	408	jpsi->SetPhiRange(0.,360.);
a2b7dc2a	409	jpsi->SetOrigin(0,0,0);
a2b7dc2a	410	jpsi->SetForceDecay(kDiMuon);
c95698be	411	return jpsi;
	412
	413	}
	414
	415	//====================================================================================================================================================
	416
a2b7dc2a	417	AliGenerator* MbPythiaTunePerugia0BtoJpsi2mu() {
c95698be	418
a2b7dc2a	419	comment = comment.Append(" pp: Pythia (Perugia0) BtoJpsi (1 bbbar per event, 1 b-hadron, 1 J/psi");
c95698be	420
	421	// Pythia
	422	AliGenPythia* pythia = new AliGenPythia(-1);
	423	pythia->SetMomentumRange(0, 999999.);
a2b7dc2a	424	pythia->SetYRange(-4.5, -2.0);
c95698be	425	pythia->SetPtRange(0,1000.);
	426	pythia->SetProcess(kPyBeautyppMNRwmi);
	427	pythia->SetEnergyCMS(energy);
	428	// Tune
	429	// 320 Perugia 0
	430	pythia->SetTune(320);
	431	pythia->UseNewMultipleInteractionsScenario();
	432	//
	433	// decays
	434	pythia->SetCutOnChild(1);
	435	pythia->SetPdgCodeParticleforAcceptanceCut(443);
a2b7dc2a	436	pythia->SetChildYRange(-4.5, -2.0);
c95698be	437	pythia->SetChildPtRange(0,10000.);
	438	//
	439	// decays
a2b7dc2a	440	pythia->SetForceDecay(kBJpsiDiMuon);
c95698be	441
c95698be	442	return pythia;
a2b7dc2a	443
c95698be	444	}
	445
	446	//====================================================================================================================================================
	447
	448	AliGenerator* HijingParam() {
	449
	450	AliGenHIJINGpara *gener = new AliGenHIJINGpara(2000);
	451	gener->SetMomentumRange(0,999);
	452	gener->SetPhiRange(0,360);
	453	gener->SetThetaRange(171, 179);
	454	gener->SetOrigin(0,0,0); // vertex position
	455	gener->SetSigma(0,0,0); // Sigma in (X,Y,Z) (cm) on IP position
	456	gGener = gener;
	457
	458	return gener;
	459
	460	}
	461
	462	//====================================================================================================================================================
	463
	464	AliGenerator* Hijing() {
	465
	466	AliGenHijing *gener = new AliGenHijing(-1);
	467	// centre of mass energy
	468	gener->SetEnergyCMS(energy);
	469	gener->SetImpactParameterRange(bMin, bMax);
	470	// reference frame
	471	gener->SetReferenceFrame("CMS");
	472	// projectile
	473	gener->SetProjectile("A", 208, 82);
	474	gener->SetTarget ("A", 208, 82);
	475	// tell hijing to keep the full parent child chain
	476	gener->KeepFullEvent();
	477	// enable jet quenching
	478	gener->SetJetQuenching(1);
	479	// enable shadowing
	480	gener->SetShadowing(1);
	481	// Don't track spectators
	482	gener->SetSpectators(0);
	483	// kinematic selection
	484	gener->SetSelectAll(0);
	485	return gener;
	486	}
	487
	488	//====================================================================================================================================================
	489
	490	AliGenerator* Hijing2500() {
	491
	492	AliGenHijing gener = (AliGenHijing) Hijing();
	493	gener->SetJetQuenching(0);
	494	gener->SetPtHardMin (4.5);
	495	return gener;
	496
	497	}
	498
	499	//====================================================================================================================================================
	500
	501	AliGenerator* Hijing2500Cocktail() {
	502
	503	comment = comment.Append(" PbPb: Hjing2500 at 5.5 + muon signals");
	504
	505	AliGenCocktail *cocktail = new AliGenCocktail();
	506	cocktail->SetProjectile("A", 208, 82);
	507	cocktail->SetTarget ("A", 208, 82);
508	cocktail->SetEnergyCMS(energy);
509
510	// 1 Hijing event: provides underlying event and collision geometry
511	AliGenHijing *hijing = Hijing2500();
512	cocktail->AddGenerator(hijing,"hijing",1);
513
514	// generator for the custom signal
515	AliGenerator* signalGen = 0x0;
516	switch (signal) {
517	case 0:
518	signalGen = GenBox();
519	break;
520	case 1:
521	signalGen = GenMuonLMR();
522	break;
523	case 2:
524	signalGen = GenCorrHF();
525	break;
526	default:
527	signalGen = GenBox();
528	break;
529	}
530	cocktail->AddGenerator(signalGen, "signal", 1);
531
532	cocktail->SetTrackingFlag(1);
533
534	return cocktail;
535
536	}
537
538	//====================================================================================================================================================
539
7cf09ddf AU	540	AliGenerator* QEDGeneratorPbPb() {
	541
	542	AliGenEpEmv1 *generQED = new AliGenEpEmv1();
	543	generQED->SetEnergyCMS(energy);
	544	generQED->SetProjectile("A", 208, 82);
	545	generQED->SetTarget ("A", 208, 82);
	546	generQED->SetYRange(-8.5, 0.);
	547	generQED->SetPtRange(0.4e-3, 1.0); // Set pt limits (GeV) for e+-
	548	generQED->SetDebug(0); // Set debug level (0 = silent)
	549
	550	}
	551
	552	//====================================================================================================================================================
	553
c95698be	554	void LoadLibs() {
	555
	556	#if defined(__CINT__)
	557
	558	gSystem->Load("liblhapdf"); // Parton density functions
	559	gSystem->Load("libEGPythia6"); // TGenerator interface
	560	if (proc == kPythia6) {
	561	gSystem->Load("libpythia6"); // Pythia 6.2
	562	gSystem->Load("libAliPythia6"); // ALICE specific implementations
	563	}
	564	else {
	565	gSystem->Load("libpythia6.4.21"); // Pythia 6.4
	566	gSystem->Load("libAliPythia6"); // ALICE specific implementations
	567	}
	568
	569	if (proc == kHijing \|\| proc == kHijing2500 \|\| proc == kHijing2500Cocktail) {
	570	gSystem->Load("libhijing");
	571	gSystem->Load("libTHijing");
	572	}
	573
	574	gSystem->Load("libgeant321");
	575
	576	#endif
	577
	578	}
	579
	580	//====================================================================================================================================================
	581