TDPMjet/ConfigDPMjet.C

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