MFT/Config.C

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