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