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)
5 #if !defined(__CINT__) || defined(__MAKECINT__)
8 #include <TStopwatch.h>
11 #include <TDatabasePDG.h>
12 #include <TParticle.h>
15 #include "AliGenerator.h"
17 #include "AliRunLoader.h"
20 #include "AliHeader.h"
21 #include "PYTHIA6/AliGenPythia.h"
22 #include "PYTHIA6/AliPythia.h"
25 Float_t EtaToTheta(Float_t arg);
26 void GetFinalDecayProducts(Int_t ind, AliStack & stack , TArrayI & ar);
28 void fastGen(Int_t nev = 1, char* filename = "galice.root")
30 AliPDG::AddParticlesToPdgDataBase();
31 TDatabasePDG* pPdg = TDatabasePDG::Instance();
36 AliRunLoader* rl = AliRunLoader::Open("galice.root","FASTRUN","recreate");
38 rl->SetCompressionLevel(2);
39 rl->SetNumberOfEventsPerFile(nev);
40 rl->LoadKinematics("RECREATE");
42 gAlice->SetRunLoader(rl);
46 AliStack* stack = rl->Stack();
49 AliHeader* header = rl->GetHeader();
51 // Create and Initialize Generator
53 // Example of charm generation taken from Config_PythiaHeavyFlavours.C
54 AliGenPythia *gener = new AliGenPythia(-1);
55 gener->SetEnergyCMS(14000.);
56 gener->SetMomentumRange(0,999999);
57 gener->SetPhiRange(0., 360.);
58 gener->SetThetaRange(0.,180.);
59 // gener->SetProcess(kPyCharmppMNR); // Correct Pt distribution, wrong mult
60 gener->SetProcess(kPyMb); // Correct multiplicity, wrong Pt
61 gener->SetStrucFunc(kCTEQ4L);
62 gener->SetPtHard(2.1,-1.0);
63 gener->SetFeedDownHigherFamily(kFALSE);
64 gener->SetStack(stack);
70 // Forbid some decays. Do it after gener->Init(0, because
71 // the initialization of the generator includes reading of the decay table.
73 AliPythia * py= AliPythia::Instance();
74 py->SetMDME(731,1,0); //forbid D*+->D+ + pi0
75 py->SetMDME(732,1,0);//forbid D*+->D+ + gamma
77 // Forbid all D0 decays except D0->K- pi+
78 for(Int_t d=741; d<=756; d++){
81 // decay 757 is D0->K- pi+
82 for(Int_t d=758; d<=801; d++){
92 for (Int_t iev = 0; iev < nev; iev++) {
94 cout <<"Event number "<< iev << endl;
98 rl->SetEventNumber(iev);
105 Int_t minmult = 1000;
110 //-------------------------------------------------------------------------------------
113 // Selection of events with multiplicity
114 // bigger than "minmult"
118 nprim = stack->GetNprimary();
120 for(Int_t ipart =0; ipart < nprim; ipart++){
121 TParticle * part = stack->Particle(ipart);
124 if (TMath::Abs(part->GetPdgCode())== 413) {
128 GetFinalDecayProducts(ipart,*stack,daughtersId);
130 Bool_t kineOK = kTRUE;
132 Double_t thetaMin = TMath::Pi()/4;
133 Double_t thetaMax = 3*TMath::Pi()/4;
135 for (Int_t id=1; id<=daughtersId[0]; id++) {
136 TParticle * daughter = stack->Particle(daughtersId[id]);
142 Double_t theta = daughter->Theta();
143 if (theta<thetaMin || theta>thetaMax) {
149 if (!kineOK) continue;
159 cout << "Number of particles " << nprim << endl;
160 cout << "Number of trials " << ntrial << endl;
163 header->SetNprimary(stack->GetNprimary());
164 header->SetNtrack(stack->GetNtrack());
167 stack->FinishEvent();
168 header->SetStack(stack);
170 rl->WriteKinematics("OVERWRITE");
182 rl->WriteHeader("OVERWRITE");
189 Float_t EtaToTheta(Float_t arg){
190 return (180./TMath::Pi())*2.*atan(exp(-arg));
194 void GetFinalDecayProducts(Int_t ind, AliStack & stack , TArrayI & ar){
196 // Recursive algorithm to get the final decay products of a particle
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
203 if (ind<0 || ind>stack.GetNtrack()) {
204 cerr << "Invalid index of the particle " << ind << endl;
207 if (ar.GetSize()==0) {
212 TParticle * part = stack.Particle(ind);
214 Int_t iFirstDaughter = part->GetFirstDaughter();
215 if( iFirstDaughter<0) {
216 // This particle is a final decay product, add its index to the array
218 if (ar.GetSize() <= ar[0]) ar.Set(ar.GetSize()+10); // resize if needed
223 Int_t iLastDaughter = part->GetLastDaughter();
225 for (Int_t id=iFirstDaughter; id<=iLastDaughter;id++) {
226 // Now search for final decay products of the daughters
227 GetFinalDecayProducts(id,stack,ar);