Corrected library names and paths to macros
[u/mrichter/AliRoot.git] / test / genkine / gen / fastGen.C
CommitLineData
3928b038 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
257e1b2b 5// To be able to compile, please add -I<ALIROOT_INSTALL_PATH>/include
6// to your .rootrc
3928b038 7#if !defined(__CINT__) || defined(__MAKECINT__)
8#include <Riostream.h>
9#include <TH1F.h>
10#include <TStopwatch.h>
11#include <TDatime.h>
12#include <TRandom.h>
13#include <TDatabasePDG.h>
14#include <TParticle.h>
15#include <TArrayI.h>
eb83ad15 16#include <TTree.h>
3928b038 17
18#include "AliGenerator.h"
19#include "AliPDG.h"
20#include "AliRunLoader.h"
21#include "AliRun.h"
22#include "AliStack.h"
23#include "AliHeader.h"
257e1b2b 24#include "AliGenPythia.h"
25#include "AliPythia.h"
3928b038 26#endif
27
28Float_t EtaToTheta(Float_t arg);
29void GetFinalDecayProducts(Int_t ind, AliStack & stack , TArrayI & ar);
30
5c4f93c0 31void fastGen(Int_t nev = 1, const char* filename = "galice.root")
3928b038 32{
33 AliPDG::AddParticlesToPdgDataBase();
87feb3d7 34 TDatabasePDG::Instance();
3928b038 35
36
37
38 // Run loader
39 AliRunLoader* rl = AliRunLoader::Open("galice.root","FASTRUN","recreate");
40
41 rl->SetCompressionLevel(2);
42 rl->SetNumberOfEventsPerFile(nev);
43 rl->LoadKinematics("RECREATE");
44 rl->MakeTree("E");
45 gAlice->SetRunLoader(rl);
46
47 // Create stack
48 rl->MakeStack();
49 AliStack* stack = rl->Stack();
50
51 // Header
52 AliHeader* header = rl->GetHeader();
53
54 // Create and Initialize Generator
55
56 // Example of charm generation taken from Config_PythiaHeavyFlavours.C
57 AliGenPythia *gener = new AliGenPythia(-1);
58 gener->SetEnergyCMS(14000.);
59 gener->SetMomentumRange(0,999999);
60 gener->SetPhiRange(0., 360.);
61 gener->SetThetaRange(0.,180.);
62 // gener->SetProcess(kPyCharmppMNR); // Correct Pt distribution, wrong mult
63 gener->SetProcess(kPyMb); // Correct multiplicity, wrong Pt
64 gener->SetStrucFunc(kCTEQ4L);
65 gener->SetPtHard(2.1,-1.0);
66 gener->SetFeedDownHigherFamily(kFALSE);
67 gener->SetStack(stack);
68 gener->Init();
69
70 // Go to galice.root
71 rl->CdGAFile();
72
73 // Forbid some decays. Do it after gener->Init(0, because
74 // the initialization of the generator includes reading of the decay table.
75
76 AliPythia * py= AliPythia::Instance();
f1047fdc 77 py->SetMDME(737,1,0); //forbid D*+->D+ + pi0
78 py->SetMDME(738,1,0);//forbid D*+->D+ + gamma
3928b038 79
80 // Forbid all D0 decays except D0->K- pi+
f1047fdc 81 for(Int_t d=747; d<=762; d++){
3928b038 82 py->SetMDME(d,1,0);
83 }
f1047fdc 84 // decay 763 is D0->K- pi+
85 for(Int_t d=764; d<=807; d++){
3928b038 86 py->SetMDME(d,1,0);
87 }
88
89 //
90 // Event Loop
91 //
92
93 TStopwatch timer;
94 timer.Start();
95 for (Int_t iev = 0; iev < nev; iev++) {
96
97 cout <<"Event number "<< iev << endl;
98
99 // Initialize event
100 header->Reset(0,iev);
101 rl->SetEventNumber(iev);
102 stack->Reset();
103 rl->MakeTree("K");
104
105 // Generate event
106 Int_t nprim = 0;
107 Int_t ntrial = 0;
3928b038 108 Int_t ndstar = 0;
109
110
111
112 //-------------------------------------------------------------------------------------
113
114 while(!ndstar) {
87feb3d7 115 // Selection of events with D*
3928b038 116 stack->Reset();
87feb3d7 117 stack->ConnectTree(rl->TreeK());
3928b038 118 gener->Generate();
119 ntrial++;
120 nprim = stack->GetNprimary();
121
122 for(Int_t ipart =0; ipart < nprim; ipart++){
123 TParticle * part = stack->Particle(ipart);
124 if(part) {
125
126 if (TMath::Abs(part->GetPdgCode())== 413) {
127
128 TArrayI daughtersId;
129
130 GetFinalDecayProducts(ipart,*stack,daughtersId);
131
132 Bool_t kineOK = kTRUE;
133
134 Double_t thetaMin = TMath::Pi()/4;
135 Double_t thetaMax = 3*TMath::Pi()/4;
136
137 for (Int_t id=1; id<=daughtersId[0]; id++) {
138 TParticle * daughter = stack->Particle(daughtersId[id]);
139 if (!daughter) {
140 kineOK = kFALSE;
141 break;
142 }
143
144 Double_t theta = daughter->Theta();
145 if (theta<thetaMin || theta>thetaMax) {
146 kineOK = kFALSE;
147 break;
148 }
149 }
150
151 if (!kineOK) continue;
152
153 part->Print();
154 ndstar++;
155
156 }
157 }
158 }
159 }
160
161 cout << "Number of particles " << nprim << endl;
162 cout << "Number of trials " << ntrial << endl;
163
164 // Finish event
165 header->SetNprimary(stack->GetNprimary());
166 header->SetNtrack(stack->GetNtrack());
167
168 // I/O
169 stack->FinishEvent();
170 header->SetStack(stack);
171 rl->TreeE()->Fill();
172 rl->WriteKinematics("OVERWRITE");
173
174 } // event loop
175 timer.Stop();
176 timer.Print();
177
178 // Termination
179 // Generator
180 gener->FinishRun();
3928b038 181 // Write file
182 rl->WriteHeader("OVERWRITE");
183 gener->Write();
184 rl->Write();
185}
186
187
188
189Float_t EtaToTheta(Float_t arg){
190 return (180./TMath::Pi())*2.*atan(exp(-arg));
191}
192
193
194void GetFinalDecayProducts(Int_t ind, AliStack & stack , TArrayI & ar){
195
196 // Recursive algorithm to get the final decay products of a particle
197 //
198 // ind is the index of the particle in the AliStack
199 // stack is the particle stack from the generator
200 // ar contains the indexes of the final decay products
201 // ar[0] is the number of final decay products
202
203 if (ind<0 || ind>stack.GetNtrack()) {
204 cerr << "Invalid index of the particle " << ind << endl;
205 return;
206 }
207 if (ar.GetSize()==0) {
208 ar.Set(10);
209 ar[0] = 0;
210 }
211
212 TParticle * part = stack.Particle(ind);
213
214 Int_t iFirstDaughter = part->GetFirstDaughter();
215 if( iFirstDaughter<0) {
216 // This particle is a final decay product, add its index to the array
217 ar[0]++;
218 if (ar.GetSize() <= ar[0]) ar.Set(ar.GetSize()+10); // resize if needed
219 ar[ar[0]] = ind;
220 return;
221 }
222
223 Int_t iLastDaughter = part->GetLastDaughter();
224
225 for (Int_t id=iFirstDaughter; id<=iLastDaughter;id++) {
226 // Now search for final decay products of the daughters
227 GetFinalDecayProducts(id,stack,ar);
228 }
229}