debug message removed
[u/mrichter/AliRoot.git] / STEER / AliKFParticleTest.C
CommitLineData
cc7cb084 1
2//---------------------------------------------------------------------------------
3// The example of usage of AliKFParticle & AliKFVertex classes for V0 analysis
4// .
5// @author S.Gorbunov, I.Kisel
6// @version 1.0
7// @since 13.05.07
8//
9// The AliKFParticleTest macro contains a toy V0 finder for ESD tracks.
10// At the first step, event primary vertex is reconstructed.
11// At the second step, ideal PID hypothesis are assigned to all the particles.
12// At the third step, V0 candidates are constructed for each pair
13// of positive and negative particles.
14// V0 candidate considered as good when it passes Chi^2 cut and
15// it is placed >= 3 Sigma away from the primary vertex.
16// Invariant mass distribution for all good V0 candidates is plotted.
17//
18// -= Copyright &copy ALICE HLT Group =-
19//_________________________________________________________________________________
20
21
22
23#if !defined( __CINT__) || defined(__MAKECINT__)
24
25 //gSystem->AddIncludePath("-I\"$ALICE_ROOT/include\"");
26
27 #include <Riostream.h>
28 #include <TTree.h>
29 #include <TFile.h>
30 #include <TCanvas.h>
31 #include <TH1D.h>
32 #include <TLatex.h>
33 #include <TDatabasePDG.h>
34 #include <TParticle.h>
35
36 #include "AliRun.h"
37 #include "AliStack.h"
86ad5fcb 38 #include "AliESDEvent.h"
cc7cb084 39 #include "AliTracker.h"
40 #include "AliKFParticle.h"
41 #include "AliKFVertex.h"
42
43#endif
44
45
46Bool_t old_file = 1;
47TCanvas *canvas=0;
48TH1D *histoMass[3]={0,0,0};
49Int_t histoPDG [4]= {22,310,3122,421};
50char* histoName[4]= {"#gamma","K^{0}_{s}","#Lambda","D^{0}"};
51TLatex *histoText[4]={0,0,0,0};
52
53void DrawV0()
54{
55 //* Draw the invariant mass histogram
56
57 if( !canvas ) return;
58 histoMass[0]->Draw();
59 for( Int_t i=1; i<3; i++ ) histoMass[i]->Draw("same");
60 for( Int_t i=0; i<4; i++ ){
61 Double_t mass = TDatabasePDG::Instance()->GetParticle(histoPDG[i])->Mass();
62 Int_t bin = histoMass[2]->FindBin(mass) -5;
63 if( bin<0 ) bin =0;
64 Double_t max = 0;
65 for( Int_t j=bin; j<bin+10; j++ )
66 if( max<histoMass[2]->GetBinContent(j) ) max = histoMass[2]->GetBinContent(j);
67 histoText[i]->SetY(max);
68 histoText[i]->SetX(mass+.05);
69 if(max>0) histoText[i]->Draw();
70 }
71 if( canvas ) canvas->Update();
72}
73
74void StartV0()
75{
76 //* (Re)create histograms and fill them from temporary file
77
78 TDirectory *curr = gDirectory;
79 canvas = new TCanvas();
80 histoMass[0] = new TH1D("massAll","AliKFParticle test", 500,0,3);
81 histoMass[0]->SetXTitle("V^{0} invariant mass [GeV]");
82 histoMass[1] = new TH1D("massMulti","V^{+-} contributions", 500,0,3);
83 histoMass[2] = new TH1D("massV0","V^{0} signal", 500,0,3);
84 histoMass[0]->SetLineColor(8);
85 histoMass[0]->SetFillColor(8);
86 histoMass[1]->SetLineColor(kMagenta);
87 histoMass[1]->SetFillColor(kMagenta);
88 histoMass[2]->SetLineColor(kBlue);
89 histoMass[2]->SetFillColor(kBlue);
90 for( Int_t i=0; i<4; i++ ){
91 histoText[i] = new TLatex(0,0,histoName[i]);
92 histoText[i]->SetTextColor(kBlue);
93 }
94 TFile* file = new TFile("AliKFParticleTest.root", "READ" );
95 if( !old_file && file ){
96 file->cd();
97 for( Int_t i=0; i<3; i++ ){
98 TObject *old = file->FindObjectAny(histoMass[i]->GetName());
99 if( old && histoMass[i] ) histoMass[i]->Add( (TH1D*) old);
100 }
101 file->Close();
102 }
103 curr->cd();
104 DrawV0();
105}
106
107void EndV0()
108{
109 //* End of gAlice -> store all histos in temporary file, clean memory
110
111 TDirectory *curr = gDirectory;
112 TFile* file = new TFile("AliKFParticleTest.root", "RECREATE" );
113 old_file = 0;
114 file->cd();
115 for( Int_t i=0; i<3; i++ ) histoMass[i]->Write();
116 file->Close();
117 for( Int_t i=0; i<3; i++ ) delete histoMass[i];
118 for( Int_t i=0; i<4; i++ ) delete histoText[i];
119 delete canvas;
120 curr->cd();
121}
122
123class TESDTrackInfo
124{
125public:
126 TESDTrackInfo(){}
127 AliKFParticle fParticle; //* assigned KFParticle
128 Bool_t fPrimUsedFlag; //* flag shows that the particle was used for primary vertex fit
129 Bool_t fOK; //* is the track good enough
130 Int_t mcPDG; //* Monte Carlo PDG code
131 Int_t mcMotherID; //* Monte Carlo ID of its mother
132};
133
86ad5fcb 134void RunV0( AliESDEvent *event )
cc7cb084 135{
136 //* V0 finder
137
138 static Int_t iEvent = 0;
139 cout<<"Event "<<iEvent++<<endl;
140
141 if( !gAlice ) return;
33c3c91a 142 AliRunLoader *rl = AliRunLoader::Instance();
cc7cb084 143
144 AliStack *stack = rl->Stack();
145 if( !stack ) return;
146
147 Int_t nESDTracks=event->GetNumberOfTracks();
148 if( nESDTracks>1000 ) nESDTracks=1000;
149
150
151 TESDTrackInfo ESDTrackInfo[1000]; //* parallel to ESD tracks
152
153 //* Fill ESDTrackInfo array
154
155 for (Int_t iTr=0; iTr<nESDTracks; iTr++)
156 {
157 TESDTrackInfo &info = ESDTrackInfo[iTr];
158 info.fOK = 0;
159 info.fPrimUsedFlag = 0;
160 info.mcPDG = -1;
161 info.mcMotherID = -1;
162
163 //* track quality check
164
165 AliESDtrack *pTrack = event->GetTrack(iTr);
166 if( !pTrack ) continue;
167 if (pTrack->GetKinkIndex(0)>0) continue;
168 if ( !( pTrack->GetStatus()&AliESDtrack::kITSrefit ) ) continue;
169 Int_t indi[12];
170 if( pTrack->GetITSclusters(indi) <5 ) continue;
6c23ffed 171 Int_t PDG = ( pTrack->GetSigned1Pt() <0 ) ?321 :211;
cc7cb084 172
173 //* take MC PDG
174 {
175 Int_t mcID = TMath::Abs(pTrack->GetLabel());
176 TParticle * part = stack->Particle(TMath::Abs(mcID));
177 if( part ){
178 info.mcPDG = part->GetPdgCode();
179 PDG = info.mcPDG;
180 if( mcID>=0 ) info.mcMotherID = part->GetFirstMother();
181 }
182 }
183
184 //* Construct KFParticle for the track
185
186 info.fParticle = AliKFParticle( *pTrack, PDG );
187 info.fOK = 1;
188 }
189
190 //* Find event primary vertex
191
192 AliKFVertex primVtx;
193 {
194 const AliKFParticle * vSelected[1000]; //* Selected particles for vertex fit
195 Int_t vIndex[1000]; //* Indices of selected particles
196 Bool_t vFlag[1000]; //* Flags returned by the vertex finder
197
198 Int_t nSelected = 0;
199 for( Int_t i = 0; i<nESDTracks; i++){
200 if(ESDTrackInfo[i].fOK ){
201 vSelected[nSelected] = &(ESDTrackInfo[i].fParticle);
202 vIndex[nSelected] = i;
203 nSelected++;
204 }
205 }
206 primVtx.ConstructPrimaryVertex( vSelected, nSelected, vFlag, 3. );
207 for( Int_t i = 0; i<nSelected; i++){
208 if( vFlag[i] ) ESDTrackInfo[vIndex[i]].fPrimUsedFlag = 1;
209 }
210 if( primVtx.GetNDF() <1 ) return; //* Less then two tracks in primary vertex
211 }
212
213 //* V0 finder
214
215 for( Int_t iTr = 0; iTr<nESDTracks; iTr++ ){ //* first daughter
216
217 TESDTrackInfo &info = ESDTrackInfo[iTr];
218 if( !info.fOK ) continue;
219
220 for( Int_t jTr = iTr+1; jTr<nESDTracks; jTr++ ){ //* second daughter
221 TESDTrackInfo &jnfo = ESDTrackInfo[jTr];
222 if( !jnfo.fOK ) continue;
223
224 //* check for different charge
225
226 if( info.fParticle.GetQ() == jnfo.fParticle.GetQ() ) continue;
227
228 //* construct V0 mother
229
616ffc76 230 AliKFParticle V0( info.fParticle, jnfo.fParticle );
cc7cb084 231
232 //* check V0 Chi^2
233
234 if( V0.GetNDF()<1 ) continue;
235 if( TMath::Sqrt(TMath::Abs(V0.GetChi2()/V0.GetNDF())) >3. ) continue;
236
237 //* subtruct daughters from primary vertex
238
239 AliKFVertex primVtxCopy = primVtx;
240
241 if( info.fPrimUsedFlag ) primVtxCopy -= info.fParticle;
242 if( jnfo.fPrimUsedFlag ) primVtxCopy -= jnfo.fParticle;
243
244 //* Check V0 Chi^2 deviation from primary vertex
245
246 if( V0.GetDeviationFromVertex( primVtxCopy ) >3. ) continue;
247
248 //* Add V0 to primary vertex to improve the primary vertex resolution
249
250 primVtxCopy += V0;
251
252 //* Set production vertex for V0
253
254 V0.SetProductionVertex( primVtxCopy );
255
256 //* Check chi^2 for a case
257
258 if( TMath::Sqrt( TMath::Abs(V0.GetChi2()/V0.GetNDF()) >3. )) continue;
259
260 //* Get V0 decay length with estimated error
261
262 Double_t length, sigmaLength;
263 if( V0.GetDecayLength( length, sigmaLength ) ) continue;
264
265 //* Reject V0 if it decays too close to the primary vertex
266
267 if( length <3.*sigmaLength ) continue;
268
269 //* Get V0 invariant mass and plot it
270
271 Double_t mass, sigmaMass;
272 if( V0.GetMass( mass, sigmaMass ) ) continue;
273 histoMass[0]->Fill(mass);
274
275 //* Fill signal histograms using MC information
276
277 if( info.mcMotherID==jnfo.mcMotherID && info.mcMotherID>=0 ){
278 TParticle *mother = stack->Particle(info.mcMotherID);
279 if( mother && TMath::Abs(mother->GetPdgCode())!=21 ){
280 histoMass[1]->Fill(mass);
281 if( mother->GetNDaughters()==2 ){
282 histoMass[2]->Fill(mass);
283 }
284 cout<<"PDG V0,pi,pj, ndaughters, mc mass, reco mass = "<<mother->GetPdgCode()<<","<<info.mcPDG<<","<<jnfo.mcPDG<<", "
285 << mother->GetNDaughters()<<", "<<mother->GetMass()<<", "<<mass<<endl;
286 }
287 }
288 }
289 }
290 if( iEvent %100 == 0 || (iEvent<100 && iEvent %10==0)) DrawV0();
291}
292
293
294
616ffc76 295Int_t AliKFParticleTest(Int_t n1=0,Int_t n2=1000,char *dire="/d/alice10/sma/my_v4-05-Release/pp/"){
cc7cb084 296 //* Main macro
297
298 // LOOP OVER SERIES OF DIRECTORIES
299
300 for (Int_t ifi=n1; ifi<=n2; ifi++) {
301
302 char nstring[5], filename[100], esdfile[100];
303 sprintf(nstring,"%3.3d",ifi);
304 sprintf(filename,"%s%s/galice.root",dire,nstring);
305 sprintf(esdfile,"%s%s/AliESDs.root",dire,nstring);
306
307 cout <<" Opening "<<filename<<"\nand ESD "<<esdfile<<endl;
308
309 if (gAlice) {
33c3c91a 310 delete AliRunLoader::Instance();
cc7cb084 311 delete gAlice;
312 gAlice=0;
313 }
314
315 AliRunLoader *rl = AliRunLoader::Open(filename);
316 if ( !rl ) {
317 ::Error("AliKFParticleTest.C","Can not open session !");
318 continue;
319 }
320 if (rl->LoadgAlice()) {
321 ::Error("AliKFParticleTest.C","LoadgAlice returned error !");
322 delete rl;
323 continue;
324 }
325 if (rl->LoadHeader()) {
326 ::Error("AliKFParticleTest.C","LoadHeader returned error !");
327 delete rl;
328 continue;
329 }
330 rl->LoadKinematics();
331 AliTracker::SetFieldMap(gAlice->Field(),1);
332 AliKFParticle::SetField( AliTracker::GetBz() );
333
334 //---------------------------------------//
335 // //
336 // ESD file //
337 // //
338 //---------------------------------------//
339
340 // Open file with the ESD
341 TFile *ef=TFile::Open(esdfile);
342
343 //Check if the file could be opened
344 if (!ef || !ef->IsOpen()) {cerr<<"Error open AliESDs.root !\n"; continue ;}
345
346 //create event object
86ad5fcb 347 AliESD *event = new AliESDEvent;
cc7cb084 348
349 //Set pointer to the esd tree in the file
350 TTree* tree = (TTree*) ef->Get("esdTree");
351
352 //check if the tree exists
353 if (!tree) {cerr<<"no ESD tree found\n"; continue;};
354
355 //Set pointer to the esd object in the tree
86ad5fcb 356 event->ReadFromTree(event);
cc7cb084 357
358 //Number of events
359 Int_t nevents=tree->GetEntriesFast();
360 cout << "Number of events: " << nevents << endl;
361
362 StartV0();
363 for (Int_t iev=0; iev<nevents; iev++){
364 tree->GetEvent(iev);
365 rl->GetEvent(iev);
366 RunV0(event);
367 }
368 EndV0();
369 delete event;
370 ef->Close();
371 }
372 StartV0();
373 return 0;
374}