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 <TTree.h>
  15 #include <TPDGCode.h>
  16 #include "AliGenerator.h"
  17 #include "AliPDG.h"
  18 #include "AliRunLoader.h"
  19 #include "AliDecayer.h"
  20 #include "AliDecayerPythia.h"
  21 #include "AliRun.h"
  22 #include "AliStack.h"
  23 #include "AliHeader.h"
  24 #include "AliGenAmpt.h"
  25 #endif
  26
  27 void fastGenAmpt(Int_t nev = 100, const char* filename = "galice.root")
  28 {
  29   AliPDG::AddParticlesToPdgDataBase();
  30   TDatabasePDG::Instance();
  31
  32   // Run loader
  33   AliRunLoader* rl = AliRunLoader::Open(filename,"FASTRUN","recreate");
  34   rl->SetCompressionLevel(2);
  35   rl->SetNumberOfEventsPerFile(nev);
  36   rl->LoadKinematics("RECREATE");
  37   rl->MakeTree("E");
  38   rl->SetNumberOfEventsPerRun(nev);
  39   gAlice->SetRunLoader(rl);
  40
  41   // Create stack
  42   rl->MakeStack();
  43   AliStack* stack = rl->Stack();
  44
  45   // Header
  46   AliHeader* header = rl->GetHeader();
  47
  48   AliGenAmpt *genHi = new AliGenAmpt(-1);
  49 //=============================================================================
  50 // THE DECAYER
  51   AliDecayer *decayer = new AliDecayerPythia();
  52   cout << "*****************************************" << endl;
  53   genHi->SetForceDecay( kHadronicD );
  54   genHi->SetDecayer( decayer );
  55 //=============================================================================
  56   genHi->SetEnergyCMS(2760);
  57   genHi->SetReferenceFrame("CMS");
  58   genHi->SetProjectile("A", 208, 82);
  59   genHi->SetTarget    ("A", 208, 82);
  60   genHi->SetIsoft(4); //1: defaul - 4: string melting
  61   genHi->SetPtHardMin (3);
  62   //genHi->SetImpactParameterRange(9.,9.5);
  63   genHi->SetImpactParameterRange(0.,20.0);
  64   genHi->SetNtMax(150); //NTMAX: number of timesteps (D=150)
  65   genHi->SetXmu(3.2264); //parton screening mass in fm^(-1) (D=3.2264d0)
  66   genHi->SetJetQuenching(0); // enable jet quenching
  67   genHi->SetShadowing(1);    // enable shadowing
  68   genHi->SetDecaysOff(1);    // neutral pion and heavy particle decays switched off
  69   genHi->SetSpectators(0);   // track spectators
  70   genHi->Init();
  71
  72   rl->CdGAFile();
  73
  74   TStopwatch timer;
  75   timer.Start();
  76   for (Int_t iev = 0; iev < nev; iev++) {
  77     cout <<"Event number "<< iev << endl;
  78     //  Initialize event
  79     header->Reset(0,iev);
  80     rl->SetEventNumber(iev);
  81     stack->Reset();
  82     rl->MakeTree("K");
  83     Int_t nprim = 0;
  84     Int_t ntrial = 0;
  85     Int_t ndstar = 0;
  86     stack->Reset();
  87     stack->ConnectTree(rl->TreeK());
  88     genHi->Generate();
  89     ntrial++;
  90     nprim = stack->GetNprimary();
  91     cout << "Number of particles " << nprim << endl;
  92     cout << "Number of trials " << ntrial << endl;
  93     header->SetNprimary(stack->GetNprimary());
  94     header->SetNtrack(stack->GetNtrack());
  95     stack->FinishEvent();
  96     header->SetStack(stack);
  97     rl->TreeE()->Fill();
  98     rl->WriteKinematics("OVERWRITE");
  99   } // event loop
 100   timer.Stop();
 101   timer.Print();
 102   genHi->FinishRun();
 103   rl->WriteHeader("OVERWRITE");
 104   genHi->Write();
 105   rl->Write();
 106 }