Adding TAmpt (Constantin)
[u/mrichter/AliRoot.git] / 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 "AliRun.h"
19 #include "AliStack.h"
20 #include "AliHeader.h"
21 #include "AliGenAmpt.h"
22 #endif
23
24 void fastGenAmpt(Int_t nev = 1, char* filename = "galice.root")
25 {
26   AliPDG::AddParticlesToPdgDataBase();
27   TDatabasePDG::Instance();
28
29   // Run loader
30   AliRunLoader* rl = AliRunLoader::Open(filename,"FASTRUN","recreate");
31   
32   rl->SetCompressionLevel(2);
33   rl->SetNumberOfEventsPerFile(nev);
34   rl->LoadKinematics("RECREATE");
35   rl->MakeTree("E");
36   rl->SetNumberOfEventsPerRun(nev);
37   gAlice->SetRunLoader(rl);
38   
39   // Create stack
40   rl->MakeStack();
41   AliStack* stack = rl->Stack();
42   
43   // Header
44   AliHeader* header = rl->GetHeader();
45   
46   AliGenAmpt *genHi = new AliGenAmpt(-1);
47   genHi->SetEnergyCMS(2760);
48   genHi->SetReferenceFrame("CMS");
49   genHi->SetProjectile("A", 208, 82);
50   genHi->SetTarget    ("A", 208, 82);
51   genHi->SetPtHardMin (3);
52   genHi->SetImpactParameterRange(9.,9.5);
53   genHi->SetJetQuenching(1); // enable jet quenching
54   genHi->SetShadowing(1);    // enable shadowing
55   genHi->SetDecaysOff(1);    // neutral pion and heavy particle decays switched off
56   genHi->SetSpectators(1);   // track spectators 
57   genHi->Init();
58
59   rl->CdGAFile();
60
61   TStopwatch timer;
62   timer.Start();
63   for (Int_t iev = 0; iev < nev; iev++) {
64     cout <<"Event number "<< iev << endl;
65     //  Initialize event
66     header->Reset(0,iev);
67     rl->SetEventNumber(iev);
68     stack->Reset();
69     rl->MakeTree("K");
70     Int_t nprim = 0;
71     Int_t ntrial = 0;
72     Int_t ndstar = 0;
73     stack->Reset();
74     stack->ConnectTree(rl->TreeK());
75     genHi->Generate();
76     ntrial++;
77     nprim = stack->GetNprimary();
78     cout << "Number of particles " << nprim << endl;
79     cout << "Number of trials " << ntrial << endl;
80     header->SetNprimary(stack->GetNprimary());
81     header->SetNtrack(stack->GetNtrack());  
82     stack->FinishEvent();
83     header->SetStack(stack);
84     rl->TreeE()->Fill();
85     rl->WriteKinematics("OVERWRITE");
86   } // event loop
87   timer.Stop();
88   timer.Print();
89   genHi->FinishRun();
90   rl->WriteHeader("OVERWRITE");
91   genHi->Write();
92   rl->Write();
93 }