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