d053c91812e2ecd78d1acdc9d76491b91cf08a03
[u/mrichter/AliRoot.git] / test / genkine / gen / fastGen.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
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
26 Float_t EtaToTheta(Float_t arg);
27 void GetFinalDecayProducts(Int_t ind, AliStack & stack , TArrayI & ar);
28
29 void fastGen(Int_t nev = 1, char* filename = "galice.root")
30 {
31   AliPDG::AddParticlesToPdgDataBase();
32   TDatabasePDG::Instance();
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();
75   py->SetMDME(737,1,0); //forbid D*+->D+ + pi0
76   py->SetMDME(738,1,0);//forbid D*+->D+ + gamma
77
78   // Forbid all D0 decays except D0->K- pi+
79   for(Int_t d=747; d<=762; d++){ 
80     py->SetMDME(d,1,0);
81   }
82   // decay 763 is D0->K- pi+
83   for(Int_t d=764; d<=807; d++){ 
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;
106     Int_t ndstar = 0;
107    
108    
109
110     //-------------------------------------------------------------------------------------
111
112     while(!ndstar) {
113       // Selection of events with D*
114       stack->Reset();
115       stack->ConnectTree(rl->TreeK());
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();
179   //  Write file
180   rl->WriteHeader("OVERWRITE");
181   gener->Write();
182   rl->Write();
183 }
184
185
186
187 Float_t EtaToTheta(Float_t arg){
188   return (180./TMath::Pi())*2.*atan(exp(-arg));
189 }
190
191
192 void 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 }