macros/Config_gener.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 #include <Riostream.h>
   3 #include <TRandom.h>
   4 #include <TSystem.h>
   5 #include <TVirtualMC.h>
   6 #include <TGeant3.h>
   7 #include "STEER/AliRunLoader.h"
   8 #include "STEER/AliRun.h"
   9 #include "STEER/AliConfig.h"
  10 #include "PYTHIA6/AliDecayerPythia.h"
  11 #include "EVGEN/AliGenCocktail.h"
  12 #include "EVGEN/AliGenFixed.h"
  13 #include "EVGEN/AliGenBox.h"
  14 #include "EVGEN/AliGenScan.h"
  15 #include "EVGEN/AliGenHIJINGpara.h"
  16 #include "THijing/AliGenHijing.h"
  17 #include "PYTHIA6/AliGenPythia.h"
  18 #include "THerwig/AliGenHerwig.h"
  19 #include "EVGEN/AliGenParam.h"
  20 #include "EVGEN/AliGenMUONlib.h"
  21 #include "EVGEN/AliGenPHOSlib.h"
  22 #include "EVGEN/AliGenGSIlib.h"
  23 #include "EVGEN/AliGenFLUKAsource.h"
  24 #include "EVGEN/AliGenExtFile.h"
  25 #include "EVGEN/AliGenHalo.h"
  26 #include "EVGEN/AliGenReaderTreeK.h"
  27 #include "STRUCT/AliBODY.h"
  28 #include "STRUCT/AliMAG.h"
  29 #endif
  30
  31 enum gentype_t {hijing, hijingParam, gun, box, pythia, herwig, mcatnlo,
  32                 param1, param2, param3, param4,
  33                 cocktail, fluka, halo, ntuple, scan};
  34
  35 gentype_t gentype = mcatnlo;
  36
  37 Int_t ntracks=1;
  38
  39 void Config()
  40 {
  41
  42   // Set Random Number seed
  43   gRandom->SetSeed(12345); //Set 0 to use the current time
  44   cout<<"Seed for random number generation= "<<gRandom->GetSeed()<<endl;
  45
  46
  47   // libraries required by geant321
  48 #if defined(__CINT__)
  49   gSystem->Load("libgeant321");
  50 #endif
  51
  52   new TGeant3TGeo("C++ Interface to Geant3");
  53
  54 //=======================================================================
  55 //  Create the output file
  56
  57   AliRunLoader* rl=0x0;
  58
  59   cout<<"Config.C: Creating Run Loader ..."<<endl;
  60   rl = AliRunLoader::Open("galice.root",
  61                           AliConfig::GetDefaultEventFolderName(),
  62                           "recreate");
  63   if (rl == 0x0)
  64     {
  65       gAlice->Fatal("Config.C","Can not instatiate the Run Loader");
  66       return;
  67     }
  68   rl->SetCompressionLevel(2);
  69   rl->SetNumberOfEventsPerFile(3);
  70   gAlice->SetRunLoader(rl);
  71
  72 //
  73 // Set External decayer
  74   AliDecayer* decayer = new AliDecayerPythia();
  75   decayer->SetForceDecay(kAll);
  76   decayer->Init();
  77   gMC->SetExternalDecayer(decayer);
  78
  79
  80 //
  81 //=======================================================================
  82 // ************* STEERING parameters FOR ALICE SIMULATION **************
  83 // --- Specify event type to be tracked through the ALICE setup
  84 // --- All positions are in cm, angles in degrees, and P and E in GeV
  85
  86
  87   gMC->SetProcess("DCAY",1);
  88   gMC->SetProcess("PAIR",1);
  89   gMC->SetProcess("COMP",1);
  90   gMC->SetProcess("PHOT",1);
  91   gMC->SetProcess("PFIS",0);
  92   gMC->SetProcess("DRAY",0);
  93   gMC->SetProcess("ANNI",1);
  94   gMC->SetProcess("BREM",1);
  95   gMC->SetProcess("MUNU",1);
  96   gMC->SetProcess("CKOV",1);
  97   gMC->SetProcess("HADR",1);
  98   gMC->SetProcess("LOSS",2);
  99   gMC->SetProcess("MULS",1);
 100   gMC->SetProcess("RAYL",1);
 101
 102   Float_t cut = 1.e-3;        // 1MeV cut by default
 103   Float_t tofmax = 1.e10;
 104
 105   gMC->SetCut("CUTGAM", cut);
 106   gMC->SetCut("CUTELE", cut);
 107   gMC->SetCut("CUTNEU", cut);
 108   gMC->SetCut("CUTHAD", cut);
 109   gMC->SetCut("CUTMUO", cut);
 110   gMC->SetCut("BCUTE",  cut);
 111   gMC->SetCut("BCUTM",  cut);
 112   gMC->SetCut("DCUTE",  cut);
 113   gMC->SetCut("DCUTM",  cut);
 114   gMC->SetCut("PPCUTM", cut);
 115   gMC->SetCut("TOFMAX", tofmax);
 116
 117
 118   AliGenerator * gGener = 0x0;
 119   switch(gentype)
 120     {
 121     case gun:
 122 //*********************************************
 123 // Example for Fixed Particle Gun
 124 //*********************************************
 125       {
 126         AliGenFixed *gener = new AliGenFixed(ntracks);
 127         gener->SetMomentum(50);
 128         gener->SetPhi(180.);
 129         gener->SetTheta(5.);
 130         gener->SetOrigin(0,0,0);        //vertex position
 131         gener->SetPart(kNeutron);
 132         gGener = gener;
 133       }
 134       break;
 135     case box:
 136 //*********************************************
 137 // Example for Moving Particle Gun            *
 138 //*********************************************
 139       {
 140         AliGenBox *gener = new AliGenBox(ntracks);
 141         gener->SetMomentumRange(3,4);
 142         gener->SetPhiRange(0,360);
 143         gener->SetThetaRange(90, 180. );
 144         gener->SetOrigin(0,0,0);
 145         //vertex position
 146         gener->SetSigma(0,0,0);         //Sigma in (X,Y,Z) (cm) on IP position
 147         gener->SetPart(kMuonPlus);
 148         gGener = gener;
 149       }
 150       break;
 151     case scan:
 152 //*********************************************
 153 // Scanning on a grid                         *
 154 //*********************************************
 155       {
 156         AliGenScan *gener = new AliGenScan(-1);
 157         gener->SetMomentumRange(4,4);
 158         gener->SetPhiRange(0,360);
 159         gener->SetThetaRange(9,9);
 160         //vertex position
 161         gener->SetSigma(6,6,0);         //Sigma in (X,Y,Z) (cm) on IP position
 162         gener->SetPart(kMuonPlus);
 163         gener->SetRange(20, -100, 100, 20, -100, 100, 1, 500, 500);
 164         gGener = gener;
 165       }
 166       break;
 167
 168     case hijingParam:
 169       {
 170         AliGenHIJINGpara *gener = new AliGenHIJINGpara(ntracks);
 171         gener->SetMomentumRange(0,999);
 172         gener->SetPhiRange(0,360);
 173         gener->SetThetaRange(2,10);
 174         gener->SetOrigin(0,0,0);        //vertex position
 175         gener->SetSigma(0,0,0);         //Sigma in (X,Y,Z) (cm) on IP position
 176         gGener = gener;
 177       }
 178       break;
 179     case hijing:
 180       {
 181         AliGenHijing *gener = new AliGenHijing(-1);
 182 // centre of mass energy
 183         gener->SetEnergyCMS(5500);
 184 // reference frame
 185         gener->SetReferenceFrame("CMS     ");
 186 // projectile
 187         gener->SetProjectile("A       ", 208, 82);
 188         gener->SetTarget    ("A       ", 208, 82);
 189 // impact parameter range
 190         gener->SetImpactParameterRange(0, 3.);
 191 // evaluate cross section before run
 192         gener->SetEvaluate(0);
 193 // tell hijing to keep the full parent child chain
 194         gener->KeepFullEvent();
 195 // enable jet quenching
 196         gener->SetJetQuenching(1);
 197 // enable shadowing
 198         gener->SetShadowing(1);
 199 // neutral pion and heavy particle decays switched off
 200         gener->SetDecaysOff(1);
 201 // trigger
 202         gener->SetTrigger(0);
 203 // kinematic selection
 204         gener->SetSelectAll(0);
 205 // momentum range
 206         gener->SetMomentumRange(0,999);
 207 // phi range
 208         gener->SetPhiRange(0.,360.);
 209 // theta range
 210         gener->SetThetaRange(0,180.);
 211 // select flavor (0: no, 4: charm+beauty, 5:beauty)
 212         gener->SetFlavor(0);
 213 //
 214         gener->SetOrigin(0., 0.0 ,0);
 215         gener->SetSigma(0,0,5.3);
 216         gener->SetVertexSmear(kPerEvent);
 217 // no tracking
 218         gener->SetTrackingFlag(0);
 219         gGener = gener;
 220       }
 221       break;
 222
 223     case pythia:
 224 //********************************************
 225 // Example for Charm  Production with Pythia *
 226 //********************************************
 227       {
 228         AliGenPythia *gener = new AliGenPythia(-1);
 229 //   final state kinematic cuts
 230         gener->SetMomentumRange(0,999);
 231         gener->SetPhiRange(0. ,360.);
 232         gener->SetThetaRange(0., 180.);
 233         gener->SetYRange(-10,10);
 234         gener->SetPtRange(0,100);
 235 //   vertex position and smearing
 236         gener->SetOrigin(0,0,0);       // vertex position
 237         gener->SetVertexSmear(kPerEvent);
 238         gener->SetSigma(0,0,5.6);      // Sigma in (X,Y,Z) (cm) on IP position
 239 //   Structure function. See the list in EVGEN/AliStructFuncType.h
 240         gener->SetStrucFunc(kGRVHO);
 241 // Select corection for nuclear structure functions
 242 //     gener->SetNuclei(208,208);
 243 //
 244 //   Process type. See the list in PYTHIA6/AliPythia.h
 245         gener->SetProcess(kPyBeauty);
 246 //
 247 //   Pt transfer of the hard scattering
 248         gener->SetPtHard(0.,5.);
 249 //   Decay type (semielectronic, semimuonic, nodecay)
 250         gener->SetForceDecay(kSemiElectronic);
 251 //   Centre of mass energy
 252         gener->SetEnergyCMS(5500.);
 253 //   No Tracking
 254         gener->SetTrackingFlag(0);
 255         gGener = gener;
 256       }
 257       break;
 258
 259     case herwig:
 260 //********************************************
 261 // Example for Charm  Production with Pythia *
 262 //********************************************
 263       {
 264         AliGenHerwig *gener = new AliGenHerwig(-1);
 265 //   final state kinematic cuts
 266         gener->SetMomentumRange(0,7000);
 267         gener->SetPhiRange(0. ,360.);
 268         gener->SetThetaRange(0., 180.);
 269         gener->SetYRange(-10,10);
 270         gener->SetPtRange(0,7000);
 271 //   vertex position and smearing
 272         gener->SetOrigin(0,0,0);       // vertex position
 273         gener->SetVertexSmear(kPerEvent);
 274         gener->SetSigma(0,0,5.6);      // Sigma in (X,Y,Z) (cm) on IP position
 275 //   Beam momenta
 276         gener->SetBeamMomenta(7000,7000);
 277 //   Beams
 278         gener->SetProjectile("P");
 279         gener->SetTarget("P");
 280 //   Structure function
 281         gener->SetStrucFunc(kGRVHO);
 282 //   Hard scatering
 283         gener->SetPtHardMin(200);
 284         gener->SetPtRMS(20);
 285 //   Min bias
 286         gener->SetProcess(8000);
 287 //   No Tracking
 288         gener->SetTrackingFlag(0);
 289         gGener = gener;
 290       }
 291       break;
 292
 293     case mcatnlo:
 294     {
 295       AliGenHerwig *gener = new AliGenHerwig(-1);
 296       gener->SetMomentumRange(0,7000);
 297       gener->SetPhiRange(0. ,360.);
 298       gener->SetThetaRange(0., 180.);
 299       gener->SetYRange(-10,10);
 300       gener->SetPtRange(0,7000);
 301       gener->SetOrigin(0,0,0);       // vertex position
 302       gener->SetVertexSmear(kPerEvent);
 303       gener->SetSigma(0,0,5.6);      // Sigma in (X,Y,Z) (cm) on IP position
 304       gener->SetBeamMomenta(7000,7000);
 305       gener->SetProjectile("P");
 306       gener->SetTarget("P");
 307       gener->SetStrucFunc(kCTEQ5M);
 308       gener->SetProcess(-1705);
 309       gener->SetHardProcessFile("sb.events");
 310       gener->SetEventListRange(0,1);
 311       gener->SetTrackingFlag(0);
 312       gGener = gener;
 313     }
 314     break;
 315
 316     case param1:
 317 //*******************************************************
 318 // Example for J/psi  Production from  Parameterisation
 319 // using default library (AliMUONlib)
 320 //*******************************************************
 321       {
 322         AliGenParam *gener =
 323           new AliGenParam(ntracks, AliGenMUONlib::kUpsilon);
 324         gener->SetMomentumRange(0,999);
 325         gener->SetPtRange(0,999);
 326         gener->SetPhiRange(0. , 360.);
 327         gener->SetYRange(2.5,4);
 328         gener->SetCutOnChild(1);
 329         gener->SetChildThetaRange(2,9);
 330         gener->SetOrigin(0,0,0);        //vertex position
 331         gener->SetSigma(0,0,5.3);       //Sigma in (X,Y,Z) (cm) on IP position
 332         gener->SetForceDecay(kDiMuon);
 333         gener->SetTrackingFlag(0);
 334         gGener = gener;
 335       }
 336       break;
 337
 338     case param2:
 339 //*******************************************************
 340 // Example for Omega  Production from  Parameterisation
 341 // specifying library.
 342 //*******************************************************
 343       {
 344         AliGenParam *gener = new AliGenParam(1000,new AliGenPHOSlib(),
 345                                              AliGenPHOSlib::kOmega);
 346         gener->SetWeighting(kNonAnalog);
 347         gener->SetForceDecay(kNoDecay);
 348         gener->SetPtRange(0,100);
 349         gener->SetThetaRange(45,135);
 350         gener->SetTrackingFlag(0);
 351         gGener = gener;
 352       }
 353       break;
 354
 355     case param3:
 356 //*******************************************************
 357 // Example for Upsilon  Production from  Parameterisation
 358 // specifying library.
 359 // GSI style
 360 //*******************************************************
 361       {
 362         AliGenParam *gener = new AliGenParam(1000,new AliGenGSIlib(),
 363                                              AliGenGSIlib::kUpsilon, "MUON");
 364         gener->SetMomentumRange(0,999);
 365         gener->SetPtRange(0,999);
 366         gener->SetPhiRange(0., 360.);
 367         gener->SetYRange(2.5,4);
 368         gener->SetCutOnChild(1);
 369         gener->SetChildThetaRange(2,9);
 370         gener->SetOrigin(0,0,0);        //vertex position
 371         gener->SetSigma(0,0,5.3);       //Sigma in (X,Y,Z) (cm) on IP position
 372         gener->SetForceDecay(kDiMuon);
 373         gener->SetTrackingFlag(0);
 374         gGener = gener;
 375       }
 376       break;
 377
 378     case param4:
 379 //*******************************************************
 380 // Example for Omega  Production from  Parameterisation
 381 // specifying library.
 382 // The alternative way.
 383 //*******************************************************
 384       {
 385         AliGenLib* Lib=new AliGenPHOSlib();
 386         Int_t iOmega = AliGenPHOSlib::kOmega;
 387         AliGenParam *gener = new AliGenParam(50, iOmega,
 388                                              Lib->GetPt(iOmega, ""),
 389                                              Lib->GetY (iOmega, ""),
 390                                              Lib->GetIp(iOmega, ""));
 391         gener->SetPtRange(0,999);
 392         gener->SetWeighting(kNonAnalog);
 393         gener->SetForceDecay(kNoDecay);
 394         gener->SetTrackingFlag(0);
 395         gGener = gener;
 396       }
 397       break;
 398
 399     case fluka:
 400 //*******************************************************
 401 // Example for a FLUKA Boundary Source                  *
 402 //*******************************************************
 403       {
 404         AliGenFLUKAsource *gener = new AliGenFLUKAsource(-1);
 405         gener->SetFileName("$(ALICE_ROOT)/data/all32.root");
 406         gener->SetPartFlag(9);
 407         gener->SetAgeMax(1.e-5);
 408 //  31.7 events
 409         gener->SetFraction(0.0315);
 410 //     gener->SetFraction(0.75*0.0315);
 411         rl->CdGAFile();
 412 //     gener->SetPartFlag(10);
 413         gener->SetMomentumRange(0,999);
 414         gener->SetPhiRange(0.,360.);
 415         gener->SetThetaRange(0., 180.);
 416         gener->SetAgeMax(1.e-5);
 417
 418 //  31.7 events
 419 //     gener->SetFraction(0.0315);
 420         gGener = gener;
 421       }
 422       break;
 423
 424     case ntuple:
 425 //*******************************************************
 426 // Example for reading from a external file                  *
 427 //*******************************************************
 428       {
 429         AliGenExtFile *gener = new AliGenExtFile(-1);
 430         gener->SetVertexSmear(kPerEvent);
 431         gener->SetTrackingFlag(1);
 432
 433         AliGenReaderTreeK * reader = new AliGenReaderTreeK();
 434         reader->SetFileName("$(ALICE_ROOT)/data/dtujet93.root");
 435         gener->SetReader(reader);
 436         gGener = gener;
 437       }
 438       break;
 439
 440     case halo:
 441 //*******************************************************
 442 // Example for Tunnel Halo Source                       *
 443 //*******************************************************
 444       {
 445         AliGenHalo *gener = new AliGenHalo(ntracks);
 446         gener->SetFileName("/h1/morsch/marsip/marsip5.mu");
 447         gGener = gener;
 448       }
 449       break;
 450
 451     case cocktail:
 452 //*******************************************************
 453 // Example for a Cocktail                               *
 454 //*******************************************************
 455       {
 456         AliGenCocktail *gener = new AliGenCocktail();
 457
 458         gener->SetPhiRange(0,360);
 459         gener->SetYRange(2.5,4);
 460         gener->SetThetaRange(2,9);
 461         gener->SetPtRange(0,10);
 462         gener->SetOrigin(0,0,0);        //vertex position
 463         gener->SetSigma(0,0,0);         //Sigma in (X,Y,Z) (cm) on IP position
 464         gener->SetMomentumRange(0,999);
 465
 466         AliGenParam *jpsi = new AliGenParam(1,AliGenMUONlib::kJpsi);
 467         jpsi->SetForceDecay(kDiMuon);
 468         jpsi->SetCutOnChild(1);
 469
 470
 471         AliGenFLUKAsource *bg = new AliGenFLUKAsource(-1);
 472         bg->AddFile("$(ALICE_ROOT)/data/all32.root");
 473         rl->CdGAFile();
 474         bg->SetPartFlag(9);
 475         bg->SetAgeMax(1.e-5);
 476 //  31.7 events
 477 //     gener->SetFraction(0.0315);
 478         bg->SetFraction(0.01*0.0315);
 479
 480         gener->AddGenerator(jpsi,"J/Psi", 1);
 481         gener->AddGenerator(bg,"Background",1);
 482
 483         gGener = gener;
 484       }
 485       break;
 486     }
 487
 488 // Activate this line if you want the vertex smearing to happen
 489 // track by track
 490 //
 491 // gener->SetVertexSmear(kPerTrack);
 492
 493   gGener->Init();
 494
 495   gAlice->SetField(-999,2);    //Specify maximum magnetic field in Tesla (neg. ==> default field)
 496
 497   Int_t iMAG=1;
 498   rl->CdGAFile();
 499
 500 //=================== Alice BODY parameters =============================
 501   AliBODY *BODY = new AliBODY("BODY","Alice envelop");
 502
 503
 504   if(iMAG) {
 505 //=================== MAG parameters ============================
 506 // --- Start with Magnet since detector layouts may be depending ---
 507 // --- on the selected Magnet dimensions ---
 508     AliMAG *MAG  = new AliMAG("MAG","Magnet");
 509   }
 510 }