TAmpt/macros/fastGenAmpt.C

   1 // Example: generation of kinematics tree with selected properties.
   2 // Below we select events containing the decays D* -> D0 pi, D0 -> K- pi+
   3 // inside the barrel part of the ALICE detector (45 < theta < 135)
   4
   5 #if !defined(__CINT__) || defined(__MAKECINT__)
   6 #include <Riostream.h>
   7 #include <TH1F.h>
   8 #include <TStopwatch.h>
   9 #include <TDatime.h>
  10 #include <TRandom.h>
  11 #include <TDatabasePDG.h>
  12 #include <TParticle.h>
  13 #include <TArrayI.h>
  14 #include <TPDGCode.h>
  15 #include "AliGenerator.h"
  16 #include "AliPDG.h"
  17 #include "AliRunLoader.h"
  18 #include "AliDecayer.h"
  19 #include "AliDecayerPythia.h"
  20 #include "AliRun.h"
  21 #include "AliStack.h"
  22 #include "AliHeader.h"
  23 #include "AliGenAmpt.h"
  24 #endif
  25
  26 void fastGenAmpt(Int_t nev = 10, const char* filename = "galice.root")
  27 {
  28   AliPDG::AddParticlesToPdgDataBase();
  29   TDatabasePDG::Instance();
  30
  31   // Run loader
  32   AliRunLoader* rl = AliRunLoader::Open(filename,"FASTRUN","recreate");
  33   rl->SetCompressionLevel(2);
  34   rl->SetNumberOfEventsPerFile(nev);
  35   rl->LoadKinematics("RECREATE");
  36   rl->MakeTree("E");
  37   rl->SetNumberOfEventsPerRun(nev);
  38   gAlice->SetRunLoader(rl);
  39
  40   // Create stack
  41   rl->MakeStack();
  42   AliStack* stack = rl->Stack();
  43
  44   // Header
  45   AliHeader* header = rl->GetHeader();
  46
  47   AliGenAmpt *genHi = new AliGenAmpt(-1);
  48 //=============================================================================
  49 // THE DECAYER
  50   AliDecayer *decayer = new AliDecayerPythia();
  51   cout << "*****************************************" << endl;
  52   genHi->SetForceDecay( kHadronicD );
  53   genHi->SetDecayer( decayer );
  54 //=============================================================================
  55   genHi->SetEnergyCMS(2760);
  56   genHi->SetReferenceFrame("CMS");
  57   genHi->SetProjectile("A", 208, 82);
  58   genHi->SetTarget    ("A", 208, 82);
  59   genHi->SetPtHardMin (3);
  60   //genHi->SetImpactParameterRange(9.,9.5);
  61   genHi->SetImpactParameterRange(0.,20.0);
  62   genHi->SetJetQuenching(1); // enable jet quenching
  63   genHi->SetShadowing(1);    // enable shadowing
  64   genHi->SetDecaysOff(1);    // neutral pion and heavy particle decays switched off
  65   genHi->SetSpectators(1);   // track spectators
  66   genHi->Init();
  67
  68   rl->CdGAFile();
  69
  70   TStopwatch timer;
  71   timer.Start();
  72   for (Int_t iev = 0; iev < nev; iev++) {
  73     cout <<"Event number "<< iev << endl;
  74     //  Initialize event
  75     header->Reset(0,iev);
  76     rl->SetEventNumber(iev);
  77     stack->Reset();
  78     rl->MakeTree("K");
  79     Int_t nprim = 0;
  80     Int_t ntrial = 0;
  81     Int_t ndstar = 0;
  82     stack->Reset();
  83     stack->ConnectTree(rl->TreeK());
  84     genHi->Generate();
  85     ntrial++;
  86     nprim = stack->GetNprimary();
  87     cout << "Number of particles " << nprim << endl;
  88     cout << "Number of trials " << ntrial << endl;
  89     header->SetNprimary(stack->GetNprimary());
  90     header->SetNtrack(stack->GetNtrack());
  91     stack->FinishEvent();
  92     header->SetStack(stack);
  93     rl->TreeE()->Fill();
  94     rl->WriteKinematics("OVERWRITE");
  95   } // event loop
  96   timer.Stop();
  97   timer.Print();
  98   genHi->FinishRun();
  99   rl->WriteHeader("OVERWRITE");
 100   genHi->Write();
 101   rl->Write();
 102 }