test/vmctest/ppbench/gen.C

   1 // $Id$
   2 //
   3 // Generation of kinematics tree which can be then used as an extrenal
   4 // event generator.
   5 // According to: $ALICE_ROOT/test/genkine/gen/fastGen.C
   6 //
   7 // By I. Hrivnacova, IPN Orsay
   8
   9 #if !defined(__CINT__) || defined(__MAKECINT__)
  10 #include <Riostream.h>
  11 #include <TH1F.h>
  12 #include <TStopwatch.h>
  13 #include <TDatime.h>
  14 #include <TRandom.h>
  15 #include <TDatabasePDG.h>
  16 #include <TParticle.h>
  17 #include <TArrayI.h>
  18
  19 #include "AliGenerator.h"
  20 #include "AliPDG.h"
  21 #include "AliRunLoader.h"
  22 #include "AliRun.h"
  23 #include "AliStack.h"
  24 #include "AliHeader.h"
  25 #include "PYTHIA6/AliGenPythia.h"
  26 #include "PYTHIA6/AliPythia.h"
  27 #endif
  28
  29 void gen(Int_t nev = 1, char* filename = "galice.root")
  30 {
  31   // Load libraries
  32   // gSystem->SetIncludePath("-I$ROOTSYS/include -I$ALICE_ROOT/include -I$ALICE_ROOT");
  33   gSystem->Load("liblhapdf.so");      // Parton density functions
  34   gSystem->Load("libEGPythia6.so");   // TGenerator interface
  35   gSystem->Load("libpythia6.so");     // Pythia
  36   gSystem->Load("libAliPythia6.so");  // ALICE specific implementations
  37
  38   AliPDG::AddParticlesToPdgDataBase();
  39   TDatabasePDG::Instance();
  40
  41   // Run loader
  42   AliRunLoader* rl = AliRunLoader::Open("galice.root","FASTRUN","recreate");
  43
  44   rl->SetCompressionLevel(2);
  45   rl->SetNumberOfEventsPerFile(nev);
  46   rl->LoadKinematics("RECREATE");
  47   rl->MakeTree("E");
  48   gAlice->SetRunLoader(rl);
  49
  50   //  Create stack
  51   rl->MakeStack();
  52   AliStack* stack = rl->Stack();
  53
  54   //  Header
  55   AliHeader* header = rl->GetHeader();
  56
  57   //  Create and Initialize Generator
  58   gROOT->LoadMacro("$ALICE_ROOT/test/vmctest/ppbench/genPPbenchConfig.C");
  59   AliGenerator* gener = genPPbenchConfig();
  60
  61   // Go to galice.root
  62   rl->CdGAFile();
  63
  64   // Forbid some decays. Do it after gener->Init(0, because
  65   // the initialization of the generator includes reading of the decay table.
  66   // ...
  67
  68   //
  69   // Event Loop
  70   //
  71
  72   TStopwatch timer;
  73   timer.Start();
  74   for (Int_t iev = 0; iev < nev; iev++) {
  75
  76     cout <<"Event number "<< iev << endl;
  77
  78     // Initialize event
  79     header->Reset(0,iev);
  80     rl->SetEventNumber(iev);
  81     stack->Reset();
  82     rl->MakeTree("K");
  83
  84     // Generate event
  85     stack->Reset();
  86     stack->ConnectTree(rl->TreeK());
  87     gener->Generate();
  88     cout << "Number of particles " << stack->GetNprimary() << endl;
  89
  90     // Finish event
  91     header->SetNprimary(stack->GetNprimary());
  92     header->SetNtrack(stack->GetNtrack());
  93
  94     // I/O
  95     stack->FinishEvent();
  96     header->SetStack(stack);
  97     rl->TreeE()->Fill();
  98     rl->WriteKinematics("OVERWRITE");
  99
 100   } // event loop
 101   timer.Stop();
 102   timer.Print();
 103
 104   //                         Termination
 105   //  Generator
 106   gener->FinishRun();
 107   //  Write file
 108   rl->WriteHeader("OVERWRITE");
 109   gener->Write();
 110   rl->Write();
 111 }