Corrected library names and paths to macros
[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 // To be able to compile, please add -I<ALIROOT_INSTALL_PATH>/include
6 // to your .rootrc
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>
16 #include <TTree.h>
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"
24 #include "AliGenPythia.h"
25 #include "AliPythia.h"
26 #endif
27
28 Float_t EtaToTheta(Float_t arg);
29 void GetFinalDecayProducts(Int_t ind, AliStack & stack , TArrayI & ar);
30
31 void fastGen(Int_t nev = 1, const char* filename = "galice.root")
32 {
33   AliPDG::AddParticlesToPdgDataBase();
34   TDatabasePDG::Instance();
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();
77   py->SetMDME(737,1,0); //forbid D*+->D+ + pi0
78   py->SetMDME(738,1,0);//forbid D*+->D+ + gamma
79
80   // Forbid all D0 decays except D0->K- pi+
81   for(Int_t d=747; d<=762; d++){ 
82     py->SetMDME(d,1,0);
83   }
84   // decay 763 is D0->K- pi+
85   for(Int_t d=764; d<=807; d++){ 
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;
108     Int_t ndstar = 0;
109    
110    
111
112     //-------------------------------------------------------------------------------------
113
114     while(!ndstar) {
115       // Selection of events with D*
116       stack->Reset();
117       stack->ConnectTree(rl->TreeK());
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();
181   //  Write file
182   rl->WriteHeader("OVERWRITE");
183   gener->Write();
184   rl->Write();
185 }
186
187
188
189 Float_t EtaToTheta(Float_t arg){
190   return (180./TMath::Pi())*2.*atan(exp(-arg));
191 }
192
193
194 void 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 }