[u/mrichter/AliRoot.git] / TRD / Macros / ConfigPID.C

//
// Configuration for the Physics Data Challenge 2006 modified to generate
// e,mu,pi,K,p for PID studies of the TRD.
// Per event, AliGenBox is used to generate 100 particles per species
// (50 positive and 50 negative ones). Barrel detectors only.
// s.masciocchi@gsi.de
//

// One can use the configuration macro in compiled mode by
// root [0] gSystem->Load("libgeant321");
// root [0] gSystem->SetIncludePath("-I$ROOTSYS/include -I$ALICE_ROOT/include\
//                   -I$ALICE_ROOT -I$ALICE/geant3/TGeant3");
// root [0] .x grun.C(1,"Config_PDC06.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 "EVGEN/AliGenCocktail.h"
#include "EVGEN/AliGenParam.h"
#include "EVGEN/AliGenMUONlib.h"
#include "STEER/AliRunLoader.h"
#include "STEER/AliRun.h"
#include "STEER/AliConfig.h"
#include "PYTHIA6/AliDecayerPythia.h"
#include "PYTHIA6/AliGenPythia.h"
#include "STEER/AliMagFMaps.h"
#include "STRUCT/AliBODY.h"
#include "STRUCT/AliMAG.h"
#include "STRUCT/AliABSOv0.h"
#include "STRUCT/AliDIPOv2.h"
#include "STRUCT/AliHALL.h"
#include "STRUCT/AliFRAMEv2.h"
#include "STRUCT/AliSHILv2.h"
#include "STRUCT/AliPIPEv0.h"
#include "ITS/AliITSgeom.h"
#include "ITS/AliITSvPPRasymmFMD.h"
#include "TPC/AliTPCv2.h"
#include "TOF/AliTOFv5T0.h"
#include "RICH/AliRICHv1.h"
#include "ZDC/AliZDCv2.h"
#include "TRD/AliTRDv1.h"
#include "FMD/AliFMDv1.h"
#include "MUON/AliMUONv1.h"
#include "PHOS/AliPHOSv1.h"
#include "PMD/AliPMDv1.h"
#include "START/AliSTARTv1.h"
#include "EMCAL/AliEMCALv2.h"
#include "CRT/AliCRTv0.h"
#include "VZERO/AliVZEROv6.h"
#endif


enum PDC06Proc_t 
{
//--- Heavy Flavour Production ---
  kCharmPbPb5500,  kCharmpPb8800,  kCharmpp14000,  kCharmpp14000wmi,
  kD0PbPb5500,     kD0pPb8800,     kD0pp14000,
  kDPlusPbPb5500,  kDPluspPb8800,  kDPluspp14000,
  kBeautyPbPb5500, kBeautypPb8800, kBeautypp14000, kBeautypp14000wmi, 
// -- Pythia Mb
  kPyMbNoHvq, kPyOmegaPlus, kPyOmegaMinus, kRunMax
}; 

const char * pprRunName[] = {
  "kCharmPbPb5500",  "kCharmpPb8800",  "kCharmpp14000",  "kCharmpp14000wmi",
  "kD0PbPb5500",     "kD0pPb8800",     "kD0pp14000",
  "kDPlusPbPb5500",  "kDPluspPb8800",  "kDPluspp14000",
  "kBeautyPbPb5500", "kBeautypPb8800", "kBeautypp14000", "kBeautypp14000wmi", 
  "kPyMbNoHvq", "kPyOmegaPlus", "kPyOmegaMinus"
};


//--- Decay Mode ---
enum DecayHvFl_t 
{
  kNature,  kHadr, kSemiEl, kSemiMu
};
//--- Rapidity Cut ---
enum YCut_t
{
  kFull, kBarrel, kMuonArm
};
//--- Magnetic Field ---
enum Mag_t
{
    k2kG, k4kG, k5kG
};
//--- Functions ---
AliGenPythia *PythiaHVQ(PDC06Proc_t proc);
AliGenerator *MbCocktail();
AliGenerator *PyMbTriggered(Int_t pdg);
void ProcessEnvironmentVars();
Float_t EtaToTheta(Float_t arg){
  return (180./TMath::Pi())*2.*atan(exp(-arg));
}

// This part for configuration
static DecayHvFl_t   decHvFl   = kNature; 
static YCut_t        ycut      = kBarrel;
static Mag_t         mag       = k5kG; 
//========================//
// Set Random Number seed //
//========================//
TDatime dt;
static UInt_t seed    = dt.Get();

// 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()
{
 

  // Get settings from environment variables
  ProcessEnvironmentVars();

  gRandom->SetSeed(seed);
  cerr<<"Seed for random number generation= "<<seed<<endl; 

  // libraries required by geant321
#if defined(__CINT__)
  gSystem->Load("libgeant321");
#endif

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

  //=======================================================================
  //  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);

  //
  //=======================================================================
  // ************* 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 //
  //=========================//
  
  AliGenCocktail *gener  = new AliGenCocktail();
  Double_t momen         = 0.6;
  Int_t Npart            = 50;
  // Set pseudorapidity range from -1 to 1.
  Float_t thmin          = EtaToTheta(1.);   // theta min. <---> eta max
  Float_t thmax          = EtaToTheta(-1.);  // theta max. <---> eta min
  
  // electron generator
  AliGenBox *eminus= new AliGenBox(Npart);
  eminus->SetPart(11);
  eminus->SetMomentumRange(momen,momen);
  eminus->SetPhiRange(0,360);
  eminus->SetThetaRange(thmin,thmax);
  eminus->SetOrigin(0., 0., 0.);    // primary vertex position
  eminus->SetSigma(0.,0.,0.);       // sigma (x,y,z) in cm, on position

  // positron generator
  AliGenBox *eplus= new AliGenBox(Npart);
  eplus->SetPart(-11);
  eplus->SetMomentumRange(momen,momen);
  eplus->SetPhiRange(0,360);
  eplus->SetThetaRange(thmin,thmax);
  eplus->SetOrigin(0., 0., 0.);    // primary vertex position
  eplus->SetSigma(0.,0.,0.);       // sigma (x,y,z) in cm, on position

  // mu- generator
  AliGenBox *muminus= new AliGenBox(Npart);
  muminus->SetPart(13);
  muminus->SetMomentumRange(momen,momen);
  muminus->SetPhiRange(0,360);
  muminus->SetThetaRange(thmin,thmax);
  muminus->SetOrigin(0., 0., 0.);    // primary vertex position
  muminus->SetSigma(0.,0.,0.);       // sigma (x,y,z) in cm, on position

  // mu+ generator
  AliGenBox *muplus= new AliGenBox(Npart);
  muplus->SetPart(-13);
  muplus->SetMomentumRange(momen,momen);
  muplus->SetPhiRange(0,360);
  muplus->SetThetaRange(thmin,thmax);
  muplus->SetOrigin(0., 0., 0.);    // primary vertex position
  muplus->SetSigma(0.,0.,0.);       // sigma (x,y,z) in cm, on position

  // pi- generator
  AliGenBox *piminus= new AliGenBox(Npart);
  piminus->SetPart(-211);
  piminus->SetMomentumRange(momen,momen);
  piminus->SetPhiRange(0,360);
  piminus->SetThetaRange(thmin,thmax);
  piminus->SetOrigin(0., 0., 0.);    // primary vertex position
  piminus->SetSigma(0.,0.,0.);       // sigma (x,y,z) in cm, on position

  // pi+ generator
  AliGenBox *piplus= new AliGenBox(Npart);
  piplus->SetPart(211);
  piplus->SetMomentumRange(momen,momen);
  piplus->SetPhiRange(0,360);
  piplus->SetThetaRange(thmin,thmax);
  piplus->SetOrigin(0., 0., 0.);    // primary vertex position
  piplus->SetSigma(0.,0.,0.);       // sigma (x,y,z) in cm, on position

  // K- generator
  AliGenBox *kminus= new AliGenBox(Npart);
  kminus->SetPart(-321);
  kminus->SetMomentumRange(momen,momen);
  kminus->SetPhiRange(0,360);
  kminus->SetThetaRange(thmin,thmax);
  kminus->SetOrigin(0., 0., 0.);    // primary vertex position
  kminus->SetSigma(0.,0.,0.);       // sigma (x,y,z) in cm, on position

  // K+ generator
  AliGenBox *kplus= new AliGenBox(Npart);
  kplus->SetPart(321);
  kplus->SetMomentumRange(momen,momen);
  kplus->SetPhiRange(0,360);
  kplus->SetThetaRange(thmin,thmax);
  kplus->SetOrigin(0., 0., 0.);    // primary vertex position
  kplus->SetSigma(0.,0.,0.);       // sigma (x,y,z) in cm, on position

  // p generator
  AliGenBox *proton= new AliGenBox(Npart);
  proton->SetPart(2212);
  proton->SetMomentumRange(momen,momen);
  proton->SetPhiRange(0,360);
  proton->SetThetaRange(thmin,thmax);
  proton->SetOrigin(0., 0., 0.);    // primary vertex position
  proton->SetSigma(0.,0.,0.);       // sigma (x,y,z) in cm, on position

  // anti-p generator
  AliGenBox *aproton= new AliGenBox(Npart);
  aproton->SetPart(-2212);
  aproton->SetMomentumRange(momen,momen);
  aproton->SetPhiRange(0,360);
  aproton->SetThetaRange(thmin,thmax);
  aproton->SetOrigin(0., 0., 0.);    // primary vertex position
  aproton->SetSigma(0.,0.,0.);       // sigma (x,y,z) in cm, on position


  gener->AddGenerator(eminus,"e-",1);
  gener->AddGenerator(eplus,"e+",1);
  gener->AddGenerator(muminus,"mu-",1);
  gener->AddGenerator(muplus,"mu+",1);
  gener->AddGenerator(piminus,"pi-",1);
  gener->AddGenerator(piplus,"pi+",1);
  gener->AddGenerator(kminus,"K-",1);
  gener->AddGenerator(kplus,"K+",1);
  gener->AddGenerator(proton,"p",1);
  gener->AddGenerator(aproton,"anti-p",1);

  gener->Init();

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

  Int_t iABSO  = 1;
  Int_t iCRT   = 0;
  Int_t iDIPO  = 0;
  Int_t iEMCAL = 0;
  Int_t iFMD   = 0;
  Int_t iFRAME = 1;
  Int_t iHALL  = 1;
  Int_t iITS   = 1;
  Int_t iMAG   = 1;
  Int_t iMUON  = 0;
  Int_t iPHOS  = 1;
  Int_t iPIPE  = 1;
  Int_t iPMD   = 0;
  Int_t iRICH  = 1;
  Int_t iSHIL  = 0;
  Int_t iSTART = 0;
  Int_t iTOF   = 1;
  Int_t iTPC   = 1;
  Int_t iTRD   = 1;
  Int_t iVZERO = 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 AliABSOv0("ABSO", "Muon Absorber");
    }

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

        AliDIPO *DIPO = new AliDIPOv2("DIPO", "Dipole version 2");
    }

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

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


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

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

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

        AliSHIL *SHIL = new AliSHILv2("SHIL", "Shielding Version 2");
    }


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

        AliPIPE *PIPE = new AliPIPEv0("PIPE", "Beam Pipe");
    }
 
    if(iITS) {

    //=================== ITS parameters ============================
    //
    // As the innermost detector in ALICE, the Inner Tracking System "impacts" on
    // almost all other detectors. This involves the fact that the ITS geometry
    // still has several options to be followed in parallel in order to determine
    // the best set-up which minimizes the induced background. All the geometries
    // available to date are described in the following. Read carefully the comments
    // and use the default version (the only one uncommented) unless you are making
    // comparisons and you know what you are doing. In this case just uncomment the
    // ITS geometry you want to use and run Aliroot.
    //
    // Detailed geometries:         
    //
    //
    //AliITS *ITS  = new AliITSv5symm("ITS","Updated ITS TDR detailed version with symmetric services");
    //
    //AliITS *ITS  = new AliITSv5asymm("ITS","Updates ITS TDR detailed version with asymmetric services");
    //
	AliITSvPPRasymmFMD *ITS  = new AliITSvPPRasymmFMD("ITS","New ITS PPR detailed version with asymmetric services");
	ITS->SetMinorVersion(2);  // don't touch this parameter if you're not an ITS developer
	ITS->SetReadDet(kTRUE);	  // don't touch this parameter if you're not an ITS developer
    //    ITS->SetWriteDet("$ALICE_ROOT/ITS/ITSgeometry_vPPRasymm2.det");  // don't touch this parameter if you're not an ITS developer
	ITS->SetThicknessDet1(200.);   // detector thickness on layer 1 must be in the range [100,300]
	ITS->SetThicknessDet2(200.);   // detector thickness on layer 2 must be in the range [100,300]
	ITS->SetThicknessChip1(200.);  // chip thickness on layer 1 must be in the range [150,300]
	ITS->SetThicknessChip2(200.);  // chip thickness on layer 2 must be in the range [150,300]
	ITS->SetRails(0);	     // 1 --> rails in ; 0 --> rails out
	ITS->SetCoolingFluid(1);   // 1 --> water ; 0 --> freon

    // Coarse geometries (warning: no hits are produced with these coarse geometries and they unuseful 
    // for reconstruction !):
    //                                                     
    //
    //AliITSvPPRcoarseasymm *ITS  = new AliITSvPPRcoarseasymm("ITS","New ITS PPR coarse version with asymmetric services");
    //ITS->SetRails(0);                // 1 --> rails in ; 0 --> rails out
    //ITS->SetSupportMaterial(0);      // 0 --> Copper ; 1 --> Aluminum ; 2 --> Carbon
    //
    //AliITS *ITS  = new AliITSvPPRcoarsesymm("ITS","New ITS PPR coarse version with symmetric services");
    //ITS->SetRails(0);                // 1 --> rails in ; 0 --> rails out
    //ITS->SetSupportMaterial(0);      // 0 --> Copper ; 1 --> Aluminum ; 2 --> Carbon
    //                      
    //
    //
    // Geant3 <-> EUCLID conversion
    // ============================
    //
    // SetEUCLID is a flag to output (=1) or not to output (=0) both geometry and
    // media to two ASCII files (called by default ITSgeometry.euc and
    // ITSgeometry.tme) in a format understandable to the CAD system EUCLID.
    // The default (=0) means that you dont want to use this facility.
    //
	ITS->SetEUCLID(0);  
    }

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


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


    if (iRICH)
    {
        //=================== RICH parameters ===========================
        AliRICH *RICH = new AliRICHv1("RICH", "normal RICH");

    }


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

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

    if (iTRD)
    {
        //=================== TRD parameters ============================

        AliTRD *TRD = new AliTRDv1("TRD", "TRD slow simulator");
        AliTRDgeometry *geoTRD = TRD->GetGeometry();
 
    }

    if (iFMD)
    {
        //=================== FMD parameters ============================
	AliFMD *FMD = new AliFMDv1("FMD", "normal FMD");
   }

    if (iMUON)
    {
        //=================== MUON parameters ===========================
        // New MUONv1 version (geometry defined via builders)
        AliMUON *MUON = new AliMUONv1("MUON", "default");
    }
    //=================== PHOS parameters ===========================

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


    if (iPMD)
    {
        //=================== PMD parameters ============================
        AliPMD *PMD = new AliPMDv1("PMD", "normal PMD");
    }

    if (iSTART)
    {
        //=================== START parameters ============================
        AliSTART *START = new AliSTARTv1("START", "START Detector");
    }

    if (iEMCAL)
    {
        //=================== EMCAL parameters ============================
        AliEMCAL *EMCAL = new AliEMCALv2("EMCAL", "SHISH");
    }

     if (iCRT)
    {
        //=================== CRT parameters ============================
        AliCRT *CRT = new AliCRTv0("CRT", "normal ACORDE");
    }

     if (iVZERO)
    {
        //=================== CRT parameters ============================
        AliVZERO *VZERO = new AliVZEROv6("VZERO", "normal VZERO");
    }
}


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

    // Random Number seed
    if (gSystem->Getenv("CONFIG_SEED")) {
      seed = atoi(gSystem->Getenv("CONFIG_SEED"));
    }
}
Commit	Line	Data
ad05a23c	1	//
	2	// Configuration for the Physics Data Challenge 2006 modified to generate
	3	// e,mu,pi,K,p for PID studies of the TRD.
	4	// Per event, AliGenBox is used to generate 100 particles per species
	5	// (50 positive and 50 negative ones). Barrel detectors only.
	6	// s.masciocchi@gsi.de
	7	//
	8
defc9040	9	// One can use the configuration macro in compiled mode by
	10	// root [0] gSystem->Load("libgeant321");
	11	// root [0] gSystem->SetIncludePath("-I$ROOTSYS/include -I$ALICE_ROOT/include\
	12	// -I$ALICE_ROOT -I$ALICE/geant3/TGeant3");
ad05a23c	13	// root [0] .x grun.C(1,"Config_PDC06.C++")
defc9040	14
	15	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
	16	#include <Riostream.h>
	17	#include <TRandom.h>
ad05a23c	18	#include <TDatime.h>
defc9040	19	#include <TSystem.h>
	20	#include <TVirtualMC.h>
	21	#include <TGeant3TGeo.h>
ad05a23c	22	#include "EVGEN/AliGenCocktail.h"
	23	#include "EVGEN/AliGenParam.h"
	24	#include "EVGEN/AliGenMUONlib.h"
defc9040	25	#include "STEER/AliRunLoader.h"
	26	#include "STEER/AliRun.h"
	27	#include "STEER/AliConfig.h"
defc9040	28	#include "PYTHIA6/AliDecayerPythia.h"
defc9040	29	#include "PYTHIA6/AliGenPythia.h"
	30	#include "STEER/AliMagFMaps.h"
	31	#include "STRUCT/AliBODY.h"
	32	#include "STRUCT/AliMAG.h"
	33	#include "STRUCT/AliABSOv0.h"
	34	#include "STRUCT/AliDIPOv2.h"
	35	#include "STRUCT/AliHALL.h"
	36	#include "STRUCT/AliFRAMEv2.h"
	37	#include "STRUCT/AliSHILv2.h"
	38	#include "STRUCT/AliPIPEv0.h"
ad05a23c	39	#include "ITS/AliITSgeom.h"
defc9040	40	#include "ITS/AliITSvPPRasymmFMD.h"
defc9040	41	#include "TPC/AliTPCv2.h"
ad05a23c	42	#include "TOF/AliTOFv5T0.h"
defc9040	43	#include "RICH/AliRICHv1.h"
	44	#include "ZDC/AliZDCv2.h"
	45	#include "TRD/AliTRDv1.h"
	46	#include "FMD/AliFMDv1.h"
	47	#include "MUON/AliMUONv1.h"
defc9040	48	#include "PHOS/AliPHOSv1.h"
	49	#include "PMD/AliPMDv1.h"
	50	#include "START/AliSTARTv1.h"
ad05a23c	51	#include "EMCAL/AliEMCALv2.h"
defc9040	52	#include "CRT/AliCRTv0.h"
ad05a23c	53	#include "VZERO/AliVZEROv6.h"
defc9040	54	#endif
defc9040	55
ad05a23c	56
ad05a23c	57	enum PDC06Proc_t
defc9040	58	{
ad05a23c	59	//--- Heavy Flavour Production ---
	60	kCharmPbPb5500, kCharmpPb8800, kCharmpp14000, kCharmpp14000wmi,
	61	kD0PbPb5500, kD0pPb8800, kD0pp14000,
	62	kDPlusPbPb5500, kDPluspPb8800, kDPluspp14000,
	63	kBeautyPbPb5500, kBeautypPb8800, kBeautypp14000, kBeautypp14000wmi,
	64	// -- Pythia Mb
	65	kPyMbNoHvq, kPyOmegaPlus, kPyOmegaMinus, kRunMax
	66	};
	67
	68	const char * pprRunName[] = {
	69	"kCharmPbPb5500", "kCharmpPb8800", "kCharmpp14000", "kCharmpp14000wmi",
	70	"kD0PbPb5500", "kD0pPb8800", "kD0pp14000",
	71	"kDPlusPbPb5500", "kDPluspPb8800", "kDPluspp14000",
	72	"kBeautyPbPb5500", "kBeautypPb8800", "kBeautypp14000", "kBeautypp14000wmi",
	73	"kPyMbNoHvq", "kPyOmegaPlus", "kPyOmegaMinus"
defc9040	74	};
defc9040	75
defc9040	76
ad05a23c	77	//--- Decay Mode ---
ad05a23c	78	enum DecayHvFl_t
defc9040	79	{
ad05a23c	80	kNature, kHadr, kSemiEl, kSemiMu
defc9040	81	};
ad05a23c	82	//--- Rapidity Cut ---
ad05a23c	83	enum YCut_t
defc9040	84	{
ad05a23c	85	kFull, kBarrel, kMuonArm
defc9040	86	};
ad05a23c	87	//--- Magnetic Field ---
ad05a23c	88	enum Mag_t
defc9040	89	{
	90	k2kG, k4kG, k5kG
	91	};
ad05a23c	92	//--- Functions ---
	93	AliGenPythia *PythiaHVQ(PDC06Proc_t proc);
	94	AliGenerator *MbCocktail();
	95	AliGenerator *PyMbTriggered(Int_t pdg);
defc9040	96	void ProcessEnvironmentVars();
ad05a23c	97	Float_t EtaToTheta(Float_t arg){
	98	return (180./TMath::Pi())2.atan(exp(-arg));
	99	}
	100
	101	// This part for configuration
	102	static DecayHvFl_t decHvFl = kNature;
	103	static YCut_t ycut = kBarrel;
	104	static Mag_t mag = k5kG;
	105	//========================//
	106	// Set Random Number seed //
	107	//========================//
	108	TDatime dt;
	109	static UInt_t seed = dt.Get();
	110
	111	// nEvts = -1 : you get 1 QQbar pair and all the fragmentation and
	112	// decay chain
	113	// nEvts = N>0 : you get N charm / beauty Hadrons
	114	Int_t nEvts = -1;
	115	// stars = kTRUE : all heavy resonances and their decay stored
	116	// = kFALSE: only final heavy hadrons and their decays stored
	117	Bool_t stars = kTRUE;
	118
	119	// To be used only with kCharmppMNRwmi and kBeautyppMNRwmi
	120	// To get a "reasonable" agreement with MNR results, events have to be
	121	// generated with the minimum ptHard set to 2.76 GeV.
	122	// To get a "perfect" agreement with MNR results, events have to be
	123	// generated in four ptHard bins with the following relative
	124	// normalizations:
	125	// CHARM
	126	// 2.76-3 GeV: 25%
	127	// 3-4 GeV: 40%
	128	// 4-8 GeV: 29%
	129	// >8 GeV: 6%
	130	// BEAUTY
	131	// 2.76-4 GeV: 5%
	132	// 4-6 GeV: 31%
	133	// 6-8 GeV: 28%
	134	// >8 GeV: 36%
	135	Float_t ptHardMin = 2.76;
	136	Float_t ptHardMax = -1.;
	137
	138
	139	// Comment line
	140	static TString comment;
defc9040	141
	142	void Config()
	143	{
ad05a23c	144
defc9040	145
ad05a23c	146	// Get settings from environment variables
ad05a23c	147	ProcessEnvironmentVars();
defc9040	148
ad05a23c	149	gRandom->SetSeed(seed);
ad05a23c	150	cerr<<"Seed for random number generation= "<<seed<<endl;
defc9040	151
ad05a23c	152	// libraries required by geant321
defc9040	153	#if defined(__CINT__)
ad05a23c	154	gSystem->Load("libgeant321");
defc9040	155	#endif
defc9040	156
ad05a23c	157	new TGeant3TGeo("C++ Interface to Geant3");
defc9040	158
ad05a23c	159	//=======================================================================
ad05a23c	160	// Create the output file
defc9040	161
ad05a23c	162
ad05a23c	163	AliRunLoader* rl=0x0;
defc9040	164
ad05a23c	165	cout<<"Config.C: Creating Run Loader ..."<<endl;
	166	rl = AliRunLoader::Open("galice.root",
	167	AliConfig::GetDefaultEventFolderName(),
	168	"recreate");
	169	if (rl == 0x0)
	170	{
	171	gAlice->Fatal("Config.C","Can not instatiate the Run Loader");
	172	return;
	173	}
	174	rl->SetCompressionLevel(2);
	175	rl->SetNumberOfEventsPerFile(1000);
	176	gAlice->SetRunLoader(rl);
defc9040	177
ad05a23c	178	//
	179	//=======================================================================
	180	// *********** STEERING parameters FOR ALICE SIMULATION ************
	181	// --- Specify event type to be tracked through the ALICE setup
	182	// --- All positions are in cm, angles in degrees, and P and E in GeV
defc9040	183
defc9040	184
	185	gMC->SetProcess("DCAY",1);
	186	gMC->SetProcess("PAIR",1);
	187	gMC->SetProcess("COMP",1);
	188	gMC->SetProcess("PHOT",1);
	189	gMC->SetProcess("PFIS",0);
	190	gMC->SetProcess("DRAY",0);
	191	gMC->SetProcess("ANNI",1);
	192	gMC->SetProcess("BREM",1);
	193	gMC->SetProcess("MUNU",1);
	194	gMC->SetProcess("CKOV",1);
	195	gMC->SetProcess("HADR",1);
	196	gMC->SetProcess("LOSS",2);
	197	gMC->SetProcess("MULS",1);
	198	gMC->SetProcess("RAYL",1);
	199
	200	Float_t cut = 1.e-3; // 1MeV cut by default
	201	Float_t tofmax = 1.e10;
	202
	203	gMC->SetCut("CUTGAM", cut);
	204	gMC->SetCut("CUTELE", cut);
	205	gMC->SetCut("CUTNEU", cut);
	206	gMC->SetCut("CUTHAD", cut);
	207	gMC->SetCut("CUTMUO", cut);
	208	gMC->SetCut("BCUTE", cut);
	209	gMC->SetCut("BCUTM", cut);
	210	gMC->SetCut("DCUTE", cut);
	211	gMC->SetCut("DCUTM", cut);
	212	gMC->SetCut("PPCUTM", cut);
	213	gMC->SetCut("TOFMAX", tofmax);
	214
defc9040	215
defc9040	216
ad05a23c	217
	218	// Set External decayer //
	219	//======================//
	220	TVirtualMCDecayer* decayer = new AliDecayerPythia();
	221	// DECAYS
	222	//
	223	switch(decHvFl) {
	224	case kNature:
	225	decayer->SetForceDecay(kAll);
	226	break;
	227	case kHadr:
	228	decayer->SetForceDecay(kHadronicD);
	229	break;
	230	case kSemiEl:
	231	decayer->SetForceDecay(kSemiElectronic);
	232	break;
	233	case kSemiMu:
	234	decayer->SetForceDecay(kSemiMuonic);
	235	break;
	236	}
	237	decayer->Init();
	238	gMC->SetExternalDecayer(decayer);
	239
	240	//=========================//
	241	// Generator Configuration //
	242	//=========================//
	243
	244	AliGenCocktail *gener = new AliGenCocktail();
	245	Double_t momen = 0.6;
	246	Int_t Npart = 50;
	247	// Set pseudorapidity range from -1 to 1.
	248	Float_t thmin = EtaToTheta(1.); // theta min. <---> eta max
	249	Float_t thmax = EtaToTheta(-1.); // theta max. <---> eta min
	250
	251	// electron generator
	252	AliGenBox *eminus= new AliGenBox(Npart);
	253	eminus->SetPart(11);
	254	eminus->SetMomentumRange(momen,momen);
	255	eminus->SetPhiRange(0,360);
	256	eminus->SetThetaRange(thmin,thmax);
	257	eminus->SetOrigin(0., 0., 0.); // primary vertex position
	258	eminus->SetSigma(0.,0.,0.); // sigma (x,y,z) in cm, on position
	259
	260	// positron generator
	261	AliGenBox *eplus= new AliGenBox(Npart);
	262	eplus->SetPart(-11);
	263	eplus->SetMomentumRange(momen,momen);
	264	eplus->SetPhiRange(0,360);
	265	eplus->SetThetaRange(thmin,thmax);
	266	eplus->SetOrigin(0., 0., 0.); // primary vertex position
	267	eplus->SetSigma(0.,0.,0.); // sigma (x,y,z) in cm, on position
	268
	269	// mu- generator
	270	AliGenBox *muminus= new AliGenBox(Npart);
	271	muminus->SetPart(13);
	272	muminus->SetMomentumRange(momen,momen);
	273	muminus->SetPhiRange(0,360);
	274	muminus->SetThetaRange(thmin,thmax);
	275	muminus->SetOrigin(0., 0., 0.); // primary vertex position
	276	muminus->SetSigma(0.,0.,0.); // sigma (x,y,z) in cm, on position
	277
	278	// mu+ generator
	279	AliGenBox *muplus= new AliGenBox(Npart);
	280	muplus->SetPart(-13);
281	muplus->SetMomentumRange(momen,momen);
282	muplus->SetPhiRange(0,360);
283	muplus->SetThetaRange(thmin,thmax);
284	muplus->SetOrigin(0., 0., 0.); // primary vertex position
285	muplus->SetSigma(0.,0.,0.); // sigma (x,y,z) in cm, on position
286
287	// pi- generator
288	AliGenBox *piminus= new AliGenBox(Npart);
289	piminus->SetPart(-211);
290	piminus->SetMomentumRange(momen,momen);
291	piminus->SetPhiRange(0,360);
292	piminus->SetThetaRange(thmin,thmax);
293	piminus->SetOrigin(0., 0., 0.); // primary vertex position
294	piminus->SetSigma(0.,0.,0.); // sigma (x,y,z) in cm, on position
295
296	// pi+ generator
297	AliGenBox *piplus= new AliGenBox(Npart);
298	piplus->SetPart(211);
299	piplus->SetMomentumRange(momen,momen);
300	piplus->SetPhiRange(0,360);
301	piplus->SetThetaRange(thmin,thmax);
302	piplus->SetOrigin(0., 0., 0.); // primary vertex position
303	piplus->SetSigma(0.,0.,0.); // sigma (x,y,z) in cm, on position
304
305	// K- generator
306	AliGenBox *kminus= new AliGenBox(Npart);
307	kminus->SetPart(-321);
308	kminus->SetMomentumRange(momen,momen);
309	kminus->SetPhiRange(0,360);
310	kminus->SetThetaRange(thmin,thmax);
311	kminus->SetOrigin(0., 0., 0.); // primary vertex position
312	kminus->SetSigma(0.,0.,0.); // sigma (x,y,z) in cm, on position
313
314	// K+ generator
315	AliGenBox *kplus= new AliGenBox(Npart);
316	kplus->SetPart(321);
317	kplus->SetMomentumRange(momen,momen);
318	kplus->SetPhiRange(0,360);
319	kplus->SetThetaRange(thmin,thmax);
320	kplus->SetOrigin(0., 0., 0.); // primary vertex position
321	kplus->SetSigma(0.,0.,0.); // sigma (x,y,z) in cm, on position
322
323	// p generator
324	AliGenBox *proton= new AliGenBox(Npart);
325	proton->SetPart(2212);
326	proton->SetMomentumRange(momen,momen);
327	proton->SetPhiRange(0,360);
328	proton->SetThetaRange(thmin,thmax);
329	proton->SetOrigin(0., 0., 0.); // primary vertex position
330	proton->SetSigma(0.,0.,0.); // sigma (x,y,z) in cm, on position
331
332	// anti-p generator
333	AliGenBox *aproton= new AliGenBox(Npart);
334	aproton->SetPart(-2212);
335	aproton->SetMomentumRange(momen,momen);
336	aproton->SetPhiRange(0,360);
337	aproton->SetThetaRange(thmin,thmax);
338	aproton->SetOrigin(0., 0., 0.); // primary vertex position
339	aproton->SetSigma(0.,0.,0.); // sigma (x,y,z) in cm, on position
340
341
342
343	gener->AddGenerator(eminus,"e-",1);
344	gener->AddGenerator(eplus,"e+",1);
345	gener->AddGenerator(muminus,"mu-",1);
346	gener->AddGenerator(muplus,"mu+",1);
347	gener->AddGenerator(piminus,"pi-",1);
348	gener->AddGenerator(piplus,"pi+",1);
349	gener->AddGenerator(kminus,"K-",1);
350	gener->AddGenerator(kplus,"K+",1);
351	gener->AddGenerator(proton,"p",1);
352	gener->AddGenerator(aproton,"anti-p",1);
353
354	gener->Init();
355
356	// FIELD
357	//
358	if (mag == k2kG) {
359	comment = comment.Append(" \| L3 field 0.2 T");
360	} else if (mag == k4kG) {
361	comment = comment.Append(" \| L3 field 0.4 T");
362	} else if (mag == k5kG) {
363	comment = comment.Append(" \| L3 field 0.5 T");
364	}
365	printf("\n \n Comment: %s \n \n", comment.Data());
defc9040	366
ad05a23c	367	AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 2, 1., 10., mag);
	368	field->SetL3ConstField(0); //Using const. field in the barrel
	369	rl->CdGAFile();
	370	gAlice->SetField(field);
	371
	372	Int_t iABSO = 1;
	373	Int_t iCRT = 0;
	374	Int_t iDIPO = 0;
	375	Int_t iEMCAL = 0;
	376	Int_t iFMD = 0;
	377	Int_t iFRAME = 1;
	378	Int_t iHALL = 1;
	379	Int_t iITS = 1;
	380	Int_t iMAG = 1;
	381	Int_t iMUON = 0;
	382	Int_t iPHOS = 1;
	383	Int_t iPIPE = 1;
	384	Int_t iPMD = 0;
	385	Int_t iRICH = 1;
	386	Int_t iSHIL = 0;
	387	Int_t iSTART = 0;
	388	Int_t iTOF = 1;
	389	Int_t iTPC = 1;
	390	Int_t iTRD = 1;
	391	Int_t iVZERO = 0;
	392	Int_t iZDC = 0;
	393
defc9040	394
	395	//=================== Alice BODY parameters =============================
	396	AliBODY *BODY = new AliBODY("BODY", "Alice envelop");
	397
	398
	399	if (iMAG)
	400	{
	401	//=================== MAG parameters ============================
	402	// --- Start with Magnet since detector layouts may be depending ---
	403	// --- on the selected Magnet dimensions ---
	404	AliMAG *MAG = new AliMAG("MAG", "Magnet");
	405	}
	406
	407
	408	if (iABSO)
	409	{
	410	//=================== ABSO parameters ============================
	411	AliABSO *ABSO = new AliABSOv0("ABSO", "Muon Absorber");
	412	}
	413
	414	if (iDIPO)
	415	{
	416	//=================== DIPO parameters ============================
	417
	418	AliDIPO *DIPO = new AliDIPOv2("DIPO", "Dipole version 2");
	419	}
	420
	421	if (iHALL)
	422	{
	423	//=================== HALL parameters ============================
	424
	425	AliHALL *HALL = new AliHALL("HALL", "Alice Hall");
	426	}
	427
	428
	429	if (iFRAME)
	430	{
	431	//=================== FRAME parameters ============================
	432
	433	AliFRAMEv2 *FRAME = new AliFRAMEv2("FRAME", "Space Frame");
ad05a23c	434	FRAME->SetHoles(0);
defc9040	435	}
	436
	437	if (iSHIL)
	438	{
	439	//=================== SHIL parameters ============================
	440
	441	AliSHIL *SHIL = new AliSHILv2("SHIL", "Shielding Version 2");
	442	}
	443
	444
	445	if (iPIPE)
	446	{
	447	//=================== PIPE parameters ============================
	448
	449	AliPIPE *PIPE = new AliPIPEv0("PIPE", "Beam Pipe");
	450	}
	451
	452	if(iITS) {
	453
	454	//=================== ITS parameters ============================
	455	//
	456	// As the innermost detector in ALICE, the Inner Tracking System "impacts" on
	457	// almost all other detectors. This involves the fact that the ITS geometry
	458	// still has several options to be followed in parallel in order to determine
	459	// the best set-up which minimizes the induced background. All the geometries
	460	// available to date are described in the following. Read carefully the comments
	461	// and use the default version (the only one uncommented) unless you are making
	462	// comparisons and you know what you are doing. In this case just uncomment the
	463	// ITS geometry you want to use and run Aliroot.
	464	//
	465	// Detailed geometries:
	466	//
	467	//
	468	//AliITS *ITS = new AliITSv5symm("ITS","Updated ITS TDR detailed version with symmetric services");
	469	//
	470	//AliITS *ITS = new AliITSv5asymm("ITS","Updates ITS TDR detailed version with asymmetric services");
	471	//
	472	AliITSvPPRasymmFMD *ITS = new AliITSvPPRasymmFMD("ITS","New ITS PPR detailed version with asymmetric services");
	473	ITS->SetMinorVersion(2); // don't touch this parameter if you're not an ITS developer
	474	ITS->SetReadDet(kTRUE); // don't touch this parameter if you're not an ITS developer
	475	// ITS->SetWriteDet("$ALICE_ROOT/ITS/ITSgeometry_vPPRasymm2.det"); // don't touch this parameter if you're not an ITS developer
	476	ITS->SetThicknessDet1(200.); // detector thickness on layer 1 must be in the range [100,300]
	477	ITS->SetThicknessDet2(200.); // detector thickness on layer 2 must be in the range [100,300]
	478	ITS->SetThicknessChip1(200.); // chip thickness on layer 1 must be in the range [150,300]
	479	ITS->SetThicknessChip2(200.); // chip thickness on layer 2 must be in the range [150,300]
	480	ITS->SetRails(0); // 1 --> rails in ; 0 --> rails out
	481	ITS->SetCoolingFluid(1); // 1 --> water ; 0 --> freon
	482
	483	// Coarse geometries (warning: no hits are produced with these coarse geometries and they unuseful
	484	// for reconstruction !):
	485	//
	486	//
	487	//AliITSvPPRcoarseasymm *ITS = new AliITSvPPRcoarseasymm("ITS","New ITS PPR coarse version with asymmetric services");
	488	//ITS->SetRails(0); // 1 --> rails in ; 0 --> rails out
	489	//ITS->SetSupportMaterial(0); // 0 --> Copper ; 1 --> Aluminum ; 2 --> Carbon
	490	//
	491	//AliITS *ITS = new AliITSvPPRcoarsesymm("ITS","New ITS PPR coarse version with symmetric services");
	492	//ITS->SetRails(0); // 1 --> rails in ; 0 --> rails out
	493	//ITS->SetSupportMaterial(0); // 0 --> Copper ; 1 --> Aluminum ; 2 --> Carbon
	494	//
	495	//
	496	//
	497	// Geant3 <-> EUCLID conversion
	498	// ============================
499	//
500	// SetEUCLID is a flag to output (=1) or not to output (=0) both geometry and
501	// media to two ASCII files (called by default ITSgeometry.euc and
502	// ITSgeometry.tme) in a format understandable to the CAD system EUCLID.
503	// The default (=0) means that you dont want to use this facility.
504	//
505	ITS->SetEUCLID(0);
506	}
507
508	if (iTPC)
509	{
ad05a23c	510	//============================ TPC parameters =====================
defc9040	511	AliTPC *TPC = new AliTPCv2("TPC", "Default");
defc9040	512	}
	513
	514
	515	if (iTOF) {
	516	//=================== TOF parameters ============================
ad05a23c	517	AliTOF *TOF = new AliTOFv5T0("TOF", "normal TOF");
defc9040	518	}
	519
	520
	521	if (iRICH)
	522	{
	523	//=================== RICH parameters ===========================
	524	AliRICH *RICH = new AliRICHv1("RICH", "normal RICH");
	525
	526	}
	527
	528
	529	if (iZDC)
	530	{
	531	//=================== ZDC parameters ============================
	532
	533	AliZDC *ZDC = new AliZDCv2("ZDC", "normal ZDC");
	534	}
	535
	536	if (iTRD)
	537	{
	538	//=================== TRD parameters ============================
	539
	540	AliTRD *TRD = new AliTRDv1("TRD", "TRD slow simulator");
ad05a23c	541	AliTRDgeometry *geoTRD = TRD->GetGeometry();
ad05a23c	542
defc9040	543	}
	544
	545	if (iFMD)
	546	{
	547	//=================== FMD parameters ============================
	548	AliFMD *FMD = new AliFMDv1("FMD", "normal FMD");
	549	}
	550
	551	if (iMUON)
	552	{
	553	//=================== MUON parameters ===========================
	554	// New MUONv1 version (geometry defined via builders)
	555	AliMUON *MUON = new AliMUONv1("MUON", "default");
defc9040	556	}
	557	//=================== PHOS parameters ===========================
	558
	559	if (iPHOS)
	560	{
	561	AliPHOS *PHOS = new AliPHOSv1("PHOS", "IHEP");
	562	}
	563
	564
	565	if (iPMD)
	566	{
	567	//=================== PMD parameters ============================
	568	AliPMD *PMD = new AliPMDv1("PMD", "normal PMD");
	569	}
	570
	571	if (iSTART)
	572	{
	573	//=================== START parameters ============================
	574	AliSTART *START = new AliSTARTv1("START", "START Detector");
	575	}
	576
	577	if (iEMCAL)
	578	{
	579	//=================== EMCAL parameters ============================
ad05a23c	580	AliEMCAL *EMCAL = new AliEMCALv2("EMCAL", "SHISH");
defc9040	581	}
	582
	583	if (iCRT)
	584	{
	585	//=================== CRT parameters ============================
	586	AliCRT *CRT = new AliCRTv0("CRT", "normal ACORDE");
	587	}
	588
	589	if (iVZERO)
	590	{
	591	//=================== CRT parameters ============================
ad05a23c	592	AliVZERO *VZERO = new AliVZEROv6("VZERO", "normal VZERO");
defc9040	593	}
defc9040	594	}
defc9040	595
defc9040	596
	597	void ProcessEnvironmentVars()
	598	{
	599	// Run type
	600	if (gSystem->Getenv("CONFIG_RUN_TYPE")) {
	601	for (Int_t iRun = 0; iRun < kRunMax; iRun++) {
	602	if (strcmp(gSystem->Getenv("CONFIG_RUN_TYPE"), pprRunName[iRun])==0) {
ad05a23c	603	proc = (PDC06Proc_t)iRun;
defc9040	604	cout<<"Run type set to "<<pprRunName[iRun]<<endl;
	605	}
	606	}
	607	}
	608
	609	// Random Number seed
	610	if (gSystem->Getenv("CONFIG_SEED")) {
ad05a23c	611	seed = atoi(gSystem->Getenv("CONFIG_SEED"));
defc9040	612	}
defc9040	613	}
ad05a23c	614
	615
	616