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