]> git.uio.no Git - u/mrichter/AliRoot.git/blame - HLT/TPCLib/macros/AliHLTTPCTrackerEvaluation.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
update
[u/mrichter/AliRoot.git] / HLT / TPCLib / macros / AliHLTTPCTrackerEvaluation.C
CommitLineData
924ba08d 1// $Id$
2/*
3 * Modified version of $ALICE_ROOT/TPC/AliTPCComparison.C for evaluating the
4 * performance of the CA tracker.
5 *
6 * Usage:
7 * <pre>
8 * aliroot AliHLTTPCTrackerEvaluation.C
9 * </pre>
10 *
11 * documentation to be filled in when the macro is finalized
12 *
13 *
14 * @ingroup alihlt_tpc
15 */
16#if defined(__CINT__) && !defined(__MAKECINT__)
17
18//Int_t AliHLTTPCTrackerEvaluation(const Char_t *dir=".", const char* input="AliHLTESDs.root", Float_t ptLowerCut=0.1, Float_t ptUpperCut=10., Bool_t bDedxAndClus=kFALSE)
19{
20 gSystem->AddIncludePath("-I$ROOTSYS/include -I$ALICE_ROOT/include -I$ALICE_ROOT -I$ALICE_ROOT/TPC -I$ALICE/geant3/TGeant3");
21 gROOT->LoadMacro("AliHLTTPCTrackerEvaluation.C++");
22
23 AliHLTTPCTrackerEvaluation();
24 // the attempt to pass the arguments with a self-called macro hasn't worked so far
25 //AliHLTTPCTrackerEvaluation(dir, input, ptLowerCut, ptUpperCut, bDedxAndClus);
26
27}
28#else
29
30#include <AliLog.h>
31#include <TMath.h>
32#include <TError.h>
33#include <Riostream.h>
34#include <TH1F.h>
35#include <TH2F.h>
36#include <TF1.h>
37#include <TTree.h>
38#include <TParticle.h>
39#include <TCanvas.h>
40#include <TLegend.h>
41#include <TLine.h>
42#include <TText.h>
43#include <TBenchmark.h>
44#include <TStyle.h>
45#include <TFile.h>
46#include <TROOT.h>
47
48#include "AliStack.h"
49#include "AliHeader.h"
50#include "AliTrackReference.h"
51#include "AliRunLoader.h"
52#include "AliRun.h"
53#include "AliESDEvent.h"
54#include "AliESDtrack.h"
55#include "AliTPCtrack.h"
56
57#include "AliSimDigits.h"
58#include "AliTPC.h"
59#include "AliTPCParamSR.h"
60#include "AliTPCClustersArray.h"
61#include "AliTPCClustersRow.h"
62#include "AliTPCcluster.h"
63#include "AliTPCLoader.h"
64#include "TParticlePDG.h"
65#include "TDatabasePDG.h"
66
67#include "AliCDBManager.h"
68#include "AliTPCcalibDB.h"
bda5a962 69#include "TProfile.h"
924ba08d 70
71//_______________________________________________________________________________________
72
73
74Int_t GoodTPCTracks(const Char_t *dir = ".");
75
76extern AliRun *gAlice;
77extern TBenchmark *gBenchmark;
78extern TROOT *gROOT;
79
80static Int_t allgood = 0;
81static Int_t allselected = 0;
82static Int_t allfound = 0;
83
84
85
86Int_t AliHLTTPCTrackerEvaluation(const Char_t *dir=".", const char* input="AliHLTESDs.root", Float_t ptLowerCut=0.1, Float_t ptUpperCut=10., Bool_t bDedxAndClus=kFALSE){
87
88 if(!input){
89 cerr << "please specify an input file" << endl;
90 return 1;
91 }
92
93 gBenchmark->Start("AliHLTTPCTrackerEvaluation");
94
95 ::Info("AliHLTTPCTrackerEvaluation.C","Calculating reconstruction efficiency...");
96
97 gSystem->SetIncludePath("-I$ROOTSYS/include -I$ALICE_ROOT/include -I$ALICE_ROOT -I$ALICE/geant3/TGeant3");
bda5a962 98
99 gStyle->SetTitleOffset(1.2,"X");
100 gStyle->SetTitleOffset(1.1,"Y");
101
924ba08d 102
103 // --------------- Efficiency histograms ------------------------
104
105 TH1F *hClus = new TH1F("hClus","# clusters on reconstructed tracks",300,0,300);
106
107 TH1F *hGood = new TH1F("hGood", "P_{t} distribution of MC tracks (selected for efficiency calculation)", 30, 0, 6);
108 TH1F *hFake = new TH1F("hFake", "P_{t} distribution of fake tracks", 30, 0, 6);
109 TH1F *hFound = new TH1F("hFound", "P_{t} distribution of reconstructed tracks", 30, 0, 6);
110 TH1F *hFoundMult = new TH1F("hFoundMult","P_{t} distribution of reconstructed tracks with clones", 30, 0, 6);
111 TH1F *hClone = new TH1F("hClone", "P_{t} distribution of clone tracks", 30, 0, 6);
112
113 TH1F *hEff = new TH1F("hEff", "Reconstruction efficiency based on selected MC tracks", 30, 0, 6);
114 TH1F *hFakeEff = new TH1F("hFakeEff", "Efficiency for fake tracks", 30, 0, 6);
115 TH1F *hCloneEff = new TH1F("hCloneEff","Efficiency for clone tracks", 30, 0, 6);
116
117 hEff->SetLineWidth(2);
118 hEff->SetMaximum(1.4);
119 hEff->SetYTitle("efficiency");
120 hEff->GetYaxis()->CenterTitle();
121 hEff->SetXTitle("P_{t} (GeV/c)");
122
123 hFakeEff->SetLineColor(kRed);
124 hFakeEff->SetLineWidth(2);
125
126 hCloneEff->SetLineColor(kBlue);
127 hCloneEff->SetLineWidth(2);
128
129 TH2F *hDedx = new TH2F("hDedx","dE/dx vs. momentum", 50,0.,2.,50,0.,400.);
130 hDedx->SetMarkerStyle(8);
131 hDedx->SetMarkerSize(0.4);
132 hDedx->SetXTitle("P (Gev/c)");
133 hDedx->SetYTitle("dE/dx (a.u.)");
134
135 TH1F *hnhit_ref = new TH1F("hnhit_ref","# clusters on reco tracks, used in fit performance",300,0,300);
136 TH1F *he = new TH1F("he","dE/dx for pions with 0.4<p<0.5 GeV/c",50,0.,100.);
137
138 // --------------- Resolution histograms ---------------------------
139
140 TH2F *hPt = new TH2F("hPt", "", 30, 0, 6, 50, -3, 3);
141 TH2F *hPhi = new TH2F("hPhi", "", 30, 0, 6, 50, -16, 16);
142 TH2F *hLambda = new TH2F("hLambda", "", 30, 0, 6, 50, -10, 10);
143 TH2F *hY = new TH2F("hY", "", 30, 0, 6, 30, -5, 5);
bda5a962 144 TH2F *hZ = new TH2F("hZ", "", 30, 0, 250, 30, -5, 5); // fix the x axis, if this is z dependent
924ba08d 145
bda5a962 146 TProfile *hResPt = new TProfile("hResPt", "P_{t} resolution (%)", 30, 0, 6);
147 TProfile *hResPhi = new TProfile("hResPhi", "#phi resolution (mrad)", 30, 0, 6);
148 TProfile *hResLambda = new TProfile("hResLambda", "#lambda resolution (mrad)",30, 0, 6);
149 TProfile *hResY = new TProfile("hResY", "Y resolution (mm)", 30, 0, 6);
150 TProfile *hResZ = new TProfile("hResZ", "Z resolution (mm)", 30, 0, 250); // select appropriate x axis, KKK
924ba08d 151
152 hResPt->SetXTitle("P_{t} (GeV/c)");
bda5a962 153 hResPt->SetYTitle("#sigma_{(P_{t}-P_{t_{MC}})/P_{t_{MC}}} (%)");
154 hResPt->GetYaxis()->CenterTitle(true);
924ba08d 155 hResPt->GetYaxis()->CenterTitle(true);
156
157 hResPhi->SetXTitle("P_{t} (GeV/c)");
bda5a962 158 hResPhi->SetYTitle("#sigma_{(#phi_{rec}-#phi_{MC})} (mrad)");
924ba08d 159 hResPhi->GetYaxis()->CenterTitle(true);
160
161 hResLambda->SetXTitle("P_{t} (GeV/c)");
bda5a962 162 hResLambda->SetYTitle("#sigma_{(#lambda_{rec}-#lambda_{MC})} (mrad)"); // KKK perhaps we should add the definition on the histo pad
924ba08d 163 hResLambda->GetYaxis()->CenterTitle(true);
164
165 hResY->SetXTitle("P_{t} (GeV/c)");
bda5a962 166 hResY->SetYTitle("#sigma_{(Y_{rec}-Y_{MC})} (mm)");
924ba08d 167 hResY->GetYaxis()->CenterTitle(true);
168
bda5a962 169 hResZ->SetXTitle("Z (mm)");
170 hResZ->SetYTitle("#sigma_{(Z_{rec}-Z_{MC})} (mm)");
171 hResZ->GetYaxis()->CenterTitle(true);
924ba08d 172
173 // KKK I leave the Z coordinate up to you. You can make the hResZ Z dependent.
174
bda5a962 175 /*
176
177 TH1D *hProjPt[15], *hProjPhi[15], *hProjLambda[15], *hProjY[15], *hProjZ[15]; // I am not sure we need these histos to be pointers
924ba08d 178 Char_t name[15]; Char_t title[100];
bda5a962 179
180 for(Int_t i=0; i<15; i++){
924ba08d 181 sprintf(name,"hProjPt%i",i+1);
182 sprintf(title,"(Pt_{MC}-Pt_{Rec})/Pt_{MC} @ Pt#in[%f, %f] GeV/c", i-0.5, i+0.5);
183 hProjPt[i] = new TH1D(name, title, 50, -3, -3);
184
185 sprintf(name,"hProjPhi%i",i+1);
bda5a962 186 sprintf(title,"(#phi_{MC}-#phi_{Rec}) @ #phi#in[%f, %f] GeV/c", i-0.5, i+0.5);
924ba08d 187 hProjPhi[i] = new TH1D(name, title, 50, -16, -16);
188
189 sprintf(name,"hProjLambda%i",i+1);
bda5a962 190 sprintf(title,"(#lambda_{MC}-#lambda_{Rec}) @ #lambda#in[%f, %f] GeV/c", i-0.5, i+0.5);
924ba08d 191 hProjLambda[i] = new TH1D(name, title, 50, -16, -16);
192
193 sprintf(name,"hProjY%i",i+1);
bda5a962 194 sprintf(title,"(Y_{MC}-Y_{Rec}) @ Y#in[%f, %f] GeV/c", i-0.5, i+0.5);
924ba08d 195 hProjY[i] = new TH1D(name, title, 50, -16, -16);
196
197 // here you add the hProjZ[] histograms KKK
198 }
bda5a962 199 */
924ba08d 200
201
202 // --------------- Pull variable histograms ------------------------
203
204 TH1F *hpullPhi = new TH1F("hpullPhi", "SinPhi pull", 30,-10.,10.);
205 TH1F *hpullY = new TH1F("hpullY", "Y pull", 30,-10.,10.);
206 TH1F *hpullZ = new TH1F("hpullZ", "Z pull", 30,-10.,10.);
207 TH1F *hpullDzds = new TH1F("hpullDzds","Dzds pull", 30,-10.,10.);
208 TH1F *hpullK = new TH1F("hpullK", "Kappa pull", 30,-10.,10.);
209
210//---------------------------------------------------------------------------------------
211
212 Char_t fname[100];
213 sprintf(fname,"%s/GoodTPCTracks.root",dir);
214
215 TFile *refFile = TFile::Open(fname,"old");
216 if(!refFile || !refFile->IsOpen()){
217 ::Info("AliHLTTPCTrackerEvaluation.C","Marking good tracks (will take a while)...");
218 if(GoodTPCTracks(dir)){
219 ::Error("AliHLTTPCTrackerEvaluation.C","Cannot generate the reference file!");
220 return 1;
221 }
222 }
223
224 refFile = TFile::Open(fname,"old");
225 if(!refFile || !refFile->IsOpen()){
226 ::Error("AliHLTTPCTrackerEvaluation.C","Cannot open the reference file!");
227 return 2;
228 }
229
230 //------------ access the contents of the tree created by the function GoodTracks() -------------
231
232 TTree *tpcTree = (TTree*)refFile->Get("tpcTree");
233 if(!tpcTree){
234 ::Error("AliHLTTPCTrackerEvaluation.C","Cannot get the reference tree!");
235 return 3;
236 }
237
238 TBranch *branch = tpcTree->GetBranch("TPC");
239 if(!branch){
240 ::Error("AliHLTTPCTrackerEvaluation.C","Cannot get the TPC branch!");
241 return 4;
242 }
243
244 TClonesArray dummy("AliTrackReference",1000), *refs = &dummy;
245 branch->SetAddress(&refs);
246
247 //------------ access the HLT output --------------------
248
249 sprintf(fname,"%s/%s", dir, input);
250 TFile *ef = TFile::Open(fname);
251
252 if( !ef || !ef->IsOpen() ){
253 ::Error("AliHLTTPCTrackerEvaluation.C","Cannot open [%s]!",fname);
254 return 5;
255 }
256
257 AliESDEvent *event = new AliESDEvent();
258 TTree *esdTree = NULL;
259 TString tsinput = input;
260
261 if(tsinput.CompareTo("AliESDs.root")==0){
262 esdTree = (TTree*)ef->Get("HLTesdTree");
263 if(!esdTree){
264 ::Error("AliHLTTPCTrackerEvaluation.C", "no HLTESD tree found");
265 return 6;
266 }
267 event->ReadFromTree(esdTree);
268 } else if(tsinput.CompareTo("AliHLTESDs.root")==0){
269 esdTree = (TTree*)ef->Get("esdTree");
270 if(!esdTree){
271 ::Error("AliHLTTPCTrackerEvaluation.C", "no ESD tree found");
272 return 7;
273 }
274 event->ReadFromTree(esdTree);
275 } else return 8;
276
277
278 //---------- Loop over events -------------------
279
280 Int_t iEvent = 0;
281 while(esdTree->GetEvent(iEvent)){
282
283 Int_t nEntries = event->GetNumberOfTracks();
bda5a962 284 cout << "********* Processing event number: " << iEvent << ", " << nEntries << " reconstructed tracks *******\n" << endl;
924ba08d 285
286 allfound += nEntries;
287
288 if(tpcTree->GetEvent(iEvent++)==0){
289 cerr << "No reconstructable tracks !\n";
290 continue;
291 }
292
293 Int_t ngood = refs->GetEntriesFast(); // access the track references
bda5a962 294 cout << ngood << " good MC tracks" << endl;
924ba08d 295 allgood += ngood;
296
297 const Int_t MAX = 15000;
298
299 Int_t track_notfound[MAX], itrack_notfound = 0;
300 Int_t track_fake[MAX], itrack_fake = 0;
301 Int_t track_multifound[MAX], track_multifound_n[MAX], itrack_multifound = 0;
302
303 for(Int_t i=0; i<nEntries; i++){
304
305 hClus->Fill(event->GetTrack(i)->GetTPCNcls()); // filled but not displayed
bda5a962 306 cout << "TPC label for track " << i << " = " << event->GetTrack(i)->GetTPCLabel() << ", # clusters " << event->GetTrack(i)->GetTPCNcls() << endl;
924ba08d 307 }
308
309 Int_t i;
310 for(Int_t k=0; k<ngood; k++){ // read the MC information
311
312 AliTrackReference *ref = (AliTrackReference*)refs->UncheckedAt(k);
313
314 Int_t label = ref->Label();
315 Int_t tlabel = -1;
316 Float_t ptMC = TMath::Sqrt(ref->Px()*ref->Px() + ref->Py()*ref->Py());
317 Float_t pMC = TMath::Sqrt(ref->Px()*ref->Px() + ref->Py()*ref->Py()+ref->Pz()*ref->Pz());
318
319 if (ptMC<1e-33 || ptMC<ptLowerCut || ptMC>ptUpperCut) continue; // first condition is for tracks not crossing padrow 0
bda5a962 320
924ba08d 321 allselected++;
322
323 hGood->Fill(ptMC);
324
325 for(i=0; i<nEntries; i++){
326 Int_t ttlabel = event->GetTrack(i)->GetTPCLabel();
327 if( label==TMath::Abs(ttlabel) && ttlabel<0 ){
328 track_fake[itrack_fake++] = label;
329 hFake->Fill(ptMC);
330 }
331 }
332
333 AliESDtrack *esdtrack = 0;
334 for(i=0; i<nEntries; i++){
bda5a962 335 esdtrack = event->GetTrack(i);
336 tlabel = esdtrack->GetTPCLabel();
337 if(label==tlabel) break;
924ba08d 338 }
339
340 if(i==nEntries){
bda5a962 341 track_notfound[itrack_notfound++] = label;
342 cout << "MC track " << label << " not reconstructed" << endl;
343 continue;
924ba08d 344 }
345
346 Int_t micount = 0;
347 Int_t mi;
348 AliESDtrack *mitrack;
349
350 for(mi=0; mi<nEntries; mi++){
351 mitrack = event->GetTrack(mi);
352 if(label==mitrack->GetTPCLabel()) micount++;
353 }
354
355 if(micount>1){
356 track_multifound[itrack_multifound] = label;
357 track_multifound_n[itrack_multifound] = micount;
358 itrack_multifound++;
359 hClone->Fill(ptMC,micount-1);
360 }
361
362 //if((mitrack->GetStatus()&AliESDtrack::kTPCrefit)==0) continue;
363 //if((mitrack->GetStatus()&AliESDtrack::kTPCin)==0) continue;
bda5a962 364 cout << "MC track " << label << " reconstructed "<<micount<<" times" << endl;
365
924ba08d 366 if(label==tlabel){
367 hFound->Fill(ptMC);
368 //if( micount<1 ) cout<<"ERROR!!!!"<<endl;
369 hFoundMult->Fill(ptMC,micount);
370 }
371
372 AliExternalTrackParam tpctrack =*(esdtrack->GetInnerParam());
373
374 if( TMath::Abs(tpctrack.GetSigned1Pt()) <1.e-8 ) continue;
375
376 Double_t mcxyz[3] = { ref->X(), ref->Y(), ref->Z() };
377 Double_t mclocal[3] = { ref->X(), ref->Y(), ref->Z() };
378
379 if(1){
380 Double_t c = TMath::Cos(tpctrack.GetAlpha());
381 Double_t s = TMath::Sin(tpctrack.GetAlpha());
382 Double_t x = mclocal[0];
383 mclocal[0] = x*c + mclocal[1]*s;
384 mclocal[1] =-x*s + mclocal[1]*c;
385 }
386
387 if(0){
388 Double_t local[3] = { tpctrack.GetX(),tpctrack.GetY(),tpctrack.GetZ() };
389// cout << "label: " << label << endl;
390// cout << "Z diff = " << local[2] - mclocal[2] << endl;
391// cout << "orig rec: " << local[0]<<" "<<local[1]<<" "<<local[2] <<" "<<TMath::Sqrt(local[0]*local[0]+local[1]*local[1])<<endl;
392// cout << "rotated mc: " << mclocal[0]<<" "<<mclocal[1]<<" "<<mclocal[2] <<" "<<TMath::Sqrt(mclocal[0]*mclocal[0]+mclocal[1]*mclocal[1])<<endl;
393// cout << "orig mc: " << mcxyz[0]<<" "<<mcxyz[1]<<" "<<mcxyz[2]<<endl;
394// cout << "alpha = " << tpctrack.GetAlpha()/TMath::Pi()*180.<<" "<<ref->Alpha()/TMath::Pi()*180.<<endl;
395 }
396
397 //cout << "X = " << mclocal[0] << " / " << tpctrack.GetX() << endl;
398 tpctrack.AliExternalTrackParam::PropagateTo(mclocal[0],5.);
399
400 Double_t xyz[3], pxpypz[3];
401
402 tpctrack.GetXYZ(xyz); // GetInnerXYZ
403 tpctrack.GetPxPyPz(pxpypz); // GetInnerPxPyPz
404
405 Double_t local[3] = { tpctrack.GetX(),tpctrack.GetY(),tpctrack.GetZ() };
406 Float_t phiRec = TMath::ATan2(pxpypz[1],pxpypz[0]);
407
408 if(phiRec <= TMath::Pi()) phiRec += 2*TMath::Pi();
409 if(phiRec >= TMath::Pi()) phiRec -= 2*TMath::Pi();
410
411 Double_t ptRec = TMath::Sqrt(pxpypz[0]*pxpypz[0]+pxpypz[1]*pxpypz[1]);
412 Double_t pRec = TMath::Sqrt(pxpypz[0]*pxpypz[0]+pxpypz[1]*pxpypz[1]+pxpypz[2]*pxpypz[2]);
413
414 Float_t lambdaRec = TMath::ATan2(pxpypz[2],ptRec);
415 Float_t pt_1 = 1./ptRec;
416 //Float_t pts_1 = tpctrack->GetSigned1Pt();
417 //Float_t ptMCs_1 = ref->/ptMC;
418
419
420
421 Int_t pdg = (Int_t)ref->GetLength(); // particle PDG
422 const Double_t kCLight = 0.000299792458;
423 //Double_t Bz = 5.;
424 Double_t q = 0.;
425
426 if( TMath::Abs(pdg)<10000 ){
427 TParticlePDG *ppdg = TDatabasePDG::Instance()->GetParticle(pdg);
428 if( ppdg ) q = ppdg->Charge()/3.;
429 }
430
431 Double_t kappaMC = q/ptMC; // /Bz/kCLight;
432 hnhit_ref->Fill(esdtrack->GetTPCNcls());
433
434 hPt->Fill(ptMC, (ptRec - ptMC)/(ptMC)*100.);
435
436 Float_t phiMC = TMath::ATan2(ref->Py(),ref->Px());
437 hPhi->Fill(ptMC, (phiRec - phiMC)/phiMC*1000.);
438
439 Float_t lambdaMC = TMath::ATan2(ref->Pz(),ptMC);
440 hLambda->Fill(ptMC, (lambdaRec - lambdaMC)/lambdaMC*1000.);
441
442 hY->Fill(ptMC, (local[1]-mclocal[1])*10.);
bda5a962 443 hZ->Fill( fabs(mclocal[2]), (local[2]-mclocal[2] ) *10.); // Please check hY and hZ! KKK
924ba08d 444
445 if( tpctrack.GetSigmaY2()>0 && finite(tpctrack.GetSigmaY2()) ) hpullY ->Fill( (local[1]-mclocal[1])/TMath::Sqrt(TMath::Abs(tpctrack.GetSigmaY2())) );
446 if( tpctrack.GetSigmaZ2()>0 && finite(tpctrack.GetSigmaZ2()) ) hpullZ ->Fill( (local[2]-mclocal[2])/TMath::Sqrt(TMath::Abs(tpctrack.GetSigmaZ2())) );
447 if( tpctrack.GetSigmaSnp2()>0 && finite(tpctrack.GetSigmaSnp2()) ) hpullPhi ->Fill( (phiRec-phiMC)/TMath::Sqrt(TMath::Abs(tpctrack.GetSigmaSnp2())) );
448 if( tpctrack.GetSigmaTgl2()>0 && finite(tpctrack.GetSigmaTgl2()) ) hpullDzds->Fill( (lambdaRec-lambdaMC)/TMath::Sqrt(TMath::Abs(tpctrack.GetSigmaTgl2())) );
449
450 if( tpctrack.GetSigma1Pt2()>0 && finite(tpctrack.GetSigma1Pt2()) )
451 if( q!=0 ) hpullK->Fill((tpctrack.GetSigned1Pt()-kappaMC)/TMath::Sqrt(TMath::Abs(tpctrack.GetSigma1Pt2())));
452
453 Float_t dedx = esdtrack->GetTPCsignal();
454 hDedx->Fill(pRec,dedx,1.);
455
456 if(TMath::Abs(pdg)==211) // pions
457 if(pRec>0.4 && pRec<0.5){
458 he->Fill(dedx,1.);
459 }
460 } // end of loop over "good" tracks
461
462 /*
463 cout<<"\nList of Not found tracks :\n";
464 for ( i = 0; i< itrack_notfound; i++){
465 cout<<track_notfound[i]<<"\t";
466 if ((i%5)==4) cout<<"\n";
467 }
468 cout<<"\nList of fake tracks :\n";
469 for ( i = 0; i< itrack_fake; i++){
470 cout<<track_fake[i]<<"\t";
471 if ((i%5)==4) cout<<"\n";
472 }
473 cout<<"\nList of multiple found tracks :\n";
474 for ( i=0; i<itrack_multifound; i++) {
475 cout<<"id. "<<track_multifound[i]
476 <<" found - "<<track_multifound_n[i]<<"times\n";
477 }
478
479 cout<<"Number of found tracks ="<<nEntries<<endl;
480 cout<<"Number of \"good\" tracks ="<<ngood<<endl;
481 */
482
483 refs->Clear();
484 } // end of the loop over events
485
486 delete event;
487 delete esdTree;
488 ef->Close();
489
490 delete tpcTree;
491 refFile->Close();
492
493 Stat_t nGoodMC = hGood->GetEntries();
494 Stat_t nFakes = hFake->GetEntries();
495 Stat_t nClones = hClone->GetEntries();
496 Stat_t nRec = hFound->GetEntries();
497 Stat_t nRecMult = hFoundMult->GetEntries();
498
499
500 if(nGoodMC!=0) ::Info("\n\nAliHLTTPCTrackerEvaluation","Integral efficiency in (%) is about : %f\n", nRec/nGoodMC*100.);
501 if(nRecMult+nFakes!=0) ::Info("AliHLTTPCTrackerEvaluation","Integral fake rate in (%) is about : %f\n", nFakes/(nRecMult+nFakes)*100.);
502 if(nRecMult+nFakes!=0) ::Info("AliHLTTPCTrackerEvaluation","Integral clone rate in (%) is about : %f\n", nClones/(nRecMult+nFakes)*100.);
503
504 ::Info("AliHLTTPCTrackerEvaluation", "Total number of selected MC tracks : %d\n", allgood);
505 ::Info("AliHLTTPCTrackerEvaluation", "Total number of selected MC tracks with momentum cuts : %d\n", allselected);
506 ::Info("AliHLTTPCTrackerEvaluation", "Total number of reconstructed tracks : %d\n", allfound);
507
508
509 // ---------- Plotting the histograms ----------------------------------------------------------------------------
510
511 hFound ->Sumw2();
512 hFoundMult->Sumw2();
513 hGood ->Sumw2();
514 hFake ->Sumw2();
515 hClone ->Sumw2();
516
517 hEff ->Divide(hFound,hGood, 1,1.,"b");
518 hFakeEff ->Divide(hFake, hGood, 1,1.,"b");
519 hCloneEff->Divide(hClone,hFoundMult,1,1.,"b");
520 // KKK are you sure this is the correct definition? When plotted, it looks like we have 70% clones in the lower momenta...
521 // perhaps fakes and clones should be displayed in another histogram??
522
523 gROOT->SetStyle("Plain");
bda5a962 524 gStyle->SetOptStat(""); //gStyle->SetOptStat(111110);
924ba08d 525 gStyle->SetOptFit(0); //gStyle->SetOptFit(1);
526 gStyle->SetPalette(1);
bda5a962 527
528
924ba08d 529 //--------------- canvas 0 for EFFICIENCY -------------------
530
531 TCanvas *c0 = new TCanvas("c0","reconstruction efficiency",500,450);
532
533 c0->cd();
534 hEff ->Draw();
535 hCloneEff->Draw("histsame");
536 hFakeEff ->Draw("histsame");
537
538 TLegend *leg = new TLegend(0.3,0.5,0.85,0.8);
539 leg->SetFillColor(10);
540 leg->SetLineColor(10);
541 leg->AddEntry(hEff, "reconstruction eff.", "l");
542 leg->AddEntry(hCloneEff, "clone contribution", "l");
543 leg->AddEntry(hFakeEff, "fake contribution", "l");
544 leg->Draw("same");
545
546 TLine *line1 = new TLine(0,1,6,1);
547 line1->SetLineStyle(4);
548 line1->Draw("same");
549
550 Float_t fakeData = nFakes/(nRecMult+nFakes)*100.;
551 Float_t cloneData = nClones/(nRecMult+nFakes)*100.;
552 Char_t fakeBuf[100];
553 Char_t cloneBuf[100];
554
555 sprintf(fakeBuf,"%d", (Int_t)fakeData);
556 TString fakeStr = fakeBuf;
557 fakeStr.Append(" %");
558
559 sprintf(cloneBuf,"%d", (Int_t)cloneData);
560 TString cloneStr = cloneBuf;
561 cloneStr.Append(" %");
562
563 TText *text = new TText(1.4,0.1,fakeStr.Data());
564 text->SetTextSize(0.05);
565 text->SetTextColor(kRed);
566 text->Draw();
567 text = new TText(0.6,0.3,cloneStr.Data());
568 text->SetTextSize(0.05);
569 text->SetTextColor(kBlue);
570 text->Draw();
571
572
573 //--------------- canvas 1 for RESOLUTIONS ----------------
574
575
bda5a962 576 gStyle->SetOptStat(0);
924ba08d 577 TCanvas *c1 = new TCanvas("c1","resolutions",1100,900); // please add the Y and Z projections
bda5a962 578
924ba08d 579 TF1 *fGauss = new TF1("gauss","gaus");
bda5a962 580 for(Int_t i=0; i<15; i++){
581 Float_t pMin = i*6/15.;// -.5;
582 Float_t pMax = (i+1)*6/15.;// +.5;
583 Int_t binx1 = hPt->GetXaxis()->FindBin(pMin);
584 Int_t binx2 = hPt->GetXaxis()->FindBin(pMax);
585 Float_t zMin = i*250./15.;
586 Float_t zMax = (i+1)*250./15.;
587 Int_t binz1 = hZ->GetXaxis()->FindBin(zMin);
588 Int_t binz2 = hZ->GetXaxis()->FindBin(zMax);
589 //if(i>0) binx1++;
924ba08d 590
bda5a962 591 TH1D *p;
592 p = (hPt ->ProjectionY("", binx1, binx2));
593 //cout<<i<<" "<<pMin<<" "<<pMax<<" "<<binx1<<" "<<binx2<<": "<<p->GetEntries()<<endl;
594 if(p->GetEntries()>50){
924ba08d 595 fGauss->SetParameter(2,0);
bda5a962 596 p->Fit("gauss","Q");
597 hResPt->Fill(pMin, fGauss->GetParameter(2));fGauss->GetParError(2);
598 }
599
600 //KKK I am resetting only the sigma and its error, perhaps we should do this with all the fit parameters ???
601 // SG I don't know;
602
603 p = (hPhi ->ProjectionY("", binx1, binx2));
604 if(p->GetEntries()>50){
605 fGauss->SetParameter(2,0);
606 p->Fit("gauss","Q");
607 hResPhi->Fill(pMin, fGauss->GetParameter(2)); fGauss->GetParError(2);
608 }
609
610 p = (hLambda->ProjectionY("", binx1, binx2));
611 if(p->GetEntries()>50){
612 fGauss->SetParameter(2,0);
613 p->Fit("gauss","Q");
614 hResLambda->Fill(pMin, fGauss->GetParameter(2)); fGauss->GetParError(2);
615 }
616
617 p = (hY ->ProjectionY("", binx1, binx2));
618 if(p->GetEntries()>50){
619 fGauss->SetParameter(2,0);
620 p->Fit("gauss","Q");
621 hResY->Fill(pMin, fGauss->GetParameter(2)); fGauss->GetParError(2);
622 }
924ba08d 623
bda5a962 624 p = (hZ ->ProjectionY("", binz1, binz2));
625 //cout<<i<<" "<<zMin<<" "<<zMax<<" "<<binz1<<" "<<binz2<<endl;
626 //cout<<p->GetEntries()<<endl;
627 if(p->GetEntries()>50){
628 fGauss->SetParameter(2,0);
629 p->Fit("gauss","Q");
630 hResZ->Fill(zMin, fGauss->GetParameter(2)); fGauss->GetParError(2);
631 }
924ba08d 632 }
924ba08d 633
bda5a962 634 c1->Clear(); c1->Divide(3,2);
924ba08d 635
bda5a962 636 TVirtualPad *p = c1->cd(1);
637 p->SetGridy();
638 hResPt->SetMarkerStyle(20);
639 hResPt->Draw("P"); // KKK I haven't filled the errors for the sigma, perhaps a TGraph would make them easier to plot
640
641 p = c1->cd(2);
642 p->SetGridy();
643 hResPhi->SetMarkerStyle(20);
644 hResPhi->Draw("P");
645
646 p=c1->cd(3);
647 p->SetGridy();
648 hResLambda->SetMarkerStyle(20);
649 hResLambda->Draw("P");
650
651 p=c1->cd(4);
652 p->SetGridy();
653 hResY->SetMarkerStyle(20);
654 hResY->Draw("P");
655
656 p=c1->cd(5);
657 p->SetGridy();
658 hResZ->SetMarkerStyle(20);
659 hResZ->Draw("P");
660 c1->Update();
924ba08d 661// pad 6 stays empty
662
663
664 //----------------- optional canvases 2 and 3 for dE/dx and clusters ----------------
665
666 if(bDedxAndClus){
667 TCanvas *c2 = new TCanvas("c2","dE/dx",500,450);
668 c2->cd();
669 hDedx->Draw();
670
671 TCanvas *c3 = new TCanvas("c","clusters",500,450);
672 c3->cd();
673 hClus->Draw();
674 hnhit_ref->SetLineColor(kRed);
675 hnhit_ref->Draw("same");
676 }
677
678
679 //----------------- canvas 3 for PULL VARIABLES --------------------
680
681 TCanvas *c4 = new TCanvas("c4","pull variables",1100,900);
bda5a962 682
683 gStyle->SetOptFit(2);
684 gStyle->SetOptStat("e");
685
924ba08d 686 c4->Divide(3,2);
687
688 c4->cd(1);
689 hpullY->Draw();
690 hpullY->Fit("gaus","Q");
691
692 c4->cd(2);
693 hpullZ->Draw();
694 hpullZ->Fit("gaus","Q");
695
696 c4->cd(3);
697 hpullPhi->Draw();
698 hpullPhi->Fit("gaus","Q");
699
700 c4->cd(4);
701 hpullDzds->Draw();
702 hpullDzds->Fit("gaus","Q");
703
704 c4->cd(5);
705 hpullK->Draw();
706 hpullK->Fit("gaus","Q");
707
708
709
710//-------------------- OUTPUT FILE ----------------------------------------
711
712 //TFile fc("CAtrackerEvaluation.root","RECREATE");
713 //c1->Write();
714 //c2->Write();
715 //c3->Write();
716 //fc.Close(); // KKK commented out temporarily
717
718 gBenchmark->Stop("AliHLTTPCTrackerEvaluation");
719 gBenchmark->Show("AliHLTTPCTrackerEvaluation");
720
721 return 0;
722}
723
724
725//============================================================================
726
727Int_t GoodTPCTracks(const Char_t *dir){
728
729 Char_t fname[100];
730 sprintf(fname,"%s/galice.root",dir);
731
732 AliRunLoader *runLoader = AliRunLoader::Open(fname,"COMPARISON");
733 if (!runLoader) {
734 ::Error("GoodTPCTracks","Cannot start session!");
735 return 1;
736 }
737
738 // load MC information
739
740 runLoader->LoadgAlice();
741 runLoader->LoadHeader();
742 runLoader->LoadKinematics();
743 runLoader->LoadTrackRefs();
744
745 AliTPCLoader *tpcLoader = (AliTPCLoader*)runLoader->GetLoader("TPCLoader");
746 if (tpcLoader == 0x0) {
747 ::Error("GoodTPCTracks","Cannot find TPCLoader!");
748 delete runLoader;
749 return 2;
750 }
751
752 AliTPC *TPC = (AliTPC*)runLoader->GetAliRun()->GetDetector("TPC");
753 Int_t tpcVersion = TPC->IsVersion();
754
755 //cout << "TPC version "<< tpcVersion << " has been found!\n";
756
757 if (tpcVersion==1) tpcLoader->LoadRecPoints();
758 else if (tpcVersion==2) tpcLoader->LoadDigits();
759 else {
760 ::Error("GoodTPCTracks","Wrong TPC version!");
761 delete runLoader;
762 return 3;
763 }
764
765 AliCDBManager *manager = AliCDBManager::Instance();
766 manager->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
767 manager->SetRun(0);
768
769 // loading the TPC parameters
770 AliTPCParamSR *digp = (AliTPCParamSR*)(AliTPCcalibDB::Instance()->GetParameters());
771 if (!digp) {
772 ::Error("AliHLTTPCTrackerEvaluation.C","TPC parameters have not been found!");
773 delete runLoader;
774 return 4;
775 }
776
777 Int_t nrow_up = digp->GetNRowUp(); //get the number of pad rows in up sector
778 Int_t nrows = digp->GetNRowLow()+nrow_up; //get the number of pad rows in low sector and add
779 Int_t zeroSup = digp->GetZeroSup();
780 Int_t gap = Int_t(0.125*nrows), shift = Int_t(0.5*gap);
781 Int_t good_number = Int_t(0.4*nrows); // will be used for selecting tracks with hits in 40% of the rows
782
783 Int_t nEvents = runLoader->GetNumberOfEvents();
784 ::Info("GoodTPCTracks","Number of events : %d\n", nEvents);
785
786 sprintf(fname,"%s/GoodTPCTracks.root",dir);
787 TFile *file = TFile::Open(fname,"recreate");
788
789 TClonesArray dummy("AliTrackReference",1000), *refs = &dummy;
790 TTree tpcTree("tpcTree","Tree with information about the reconstructable TPC tracks");
791 tpcTree.Branch("TPC",&refs);
792
793
794
795 //-------------- Loop over generated events ------------------------------
796
797 for(Int_t iEvent=0; iEvent<nEvents; iEvent++){
798
799 Int_t nt = 0;
800 refs->Clear();
801
802 Int_t i;
803
804 runLoader->GetEvent(iEvent);
805 file->cd();
806
807 Int_t nParticles = runLoader->GetHeader()->GetNtrack();
808 //cout << "Event " << iEvent << ", Number of particles: " << nParticles << endl;
809
810 // ------- Set the selection criteria for the MC sample that will be used for the efficiency -------
811
812 Int_t *good = new Int_t[nParticles]; for(i=0; i<nParticles; i++) good[i] = 0;
813
814 switch (tpcVersion){
815 case 1: // not used any more
816 {
817 AliTPCClustersArray *pca = new AliTPCClustersArray, &ca=*pca;
818 ca.Setup(digp);
819 ca.SetClusterType("AliTPCcluster");
820 ca.ConnectTree(tpcLoader->TreeR());
821 Int_t nrows=Int_t(ca.GetTree()->GetEntries());
822 for (Int_t n=0; n<nrows; n++) {
823 AliSegmentID *s=ca.LoadEntry(n);
824 Int_t sec,row;
825 digp->AdjustSectorRow(s->GetID(),sec,row);
826 AliTPCClustersRow &clrow = *ca.GetRow(sec,row);
827 Int_t ncl=clrow.GetArray()->GetEntriesFast();
828 while (ncl--) {
829 AliTPCcluster *c=(AliTPCcluster*)clrow[ncl];
830 Int_t label=c->GetLabel(0);
831 if (label<0) continue; //noise cluster
832 label=TMath::Abs(label);
833
834 if (sec>=digp->GetNInnerSector())
835 if (row==nrow_up-1) good[label]|=0x4000;
836 if (sec>=digp->GetNInnerSector())
837 if (row==nrow_up-1-gap) good[label]|=0x1000;
838
839 if (sec>=digp->GetNInnerSector())
840 if (row==nrow_up-1-shift) good[label]|=0x2000;
841 if (sec>=digp->GetNInnerSector())
842 if (row==nrow_up-1-gap-shift) good[label]|=0x800;
843
844 good[label]++;
845 }
846 ca.ClearRow(sec,row);
847 }
848 delete pca;
849 }
850 break;
851
852 case 2:
853 {
854 TTree *TD = tpcLoader->TreeD(); // get the digits tree
855
856 AliSimDigits da, *digits = &da;
857 TD->GetBranch("Segment")->SetAddress(&digits);
858
859 Int_t *count = new Int_t[nParticles];
860 Int_t i;
861 for (i=0; i<nParticles; i++) count[i] = 0;
862
863 Int_t sectors_by_rows = (Int_t)TD->GetEntries();
864
865 for(i=0; i<sectors_by_rows; i++){
866
867 if (!TD->GetEvent(i)) continue;
868 Int_t sec,row;
869
870 digp->AdjustSectorRow(digits->GetID(),sec,row);
871 if(digits->First()){
872 do {
873 Int_t it = digits->CurrentRow();
874 Int_t ip = digits->CurrentColumn();
875 Short_t dig = digits->GetDigit(it,ip);
876 Int_t idx0 = digits->GetTrackID(it,ip,0);
877 Int_t idx1 = digits->GetTrackID(it,ip,1);
878 Int_t idx2 = digits->GetTrackID(it,ip,2);
879 if(idx0>=0 && dig>=zeroSup && idx0<nParticles) count[idx0]+=1;
880 if(idx1>=0 && dig>=zeroSup && idx1<nParticles) count[idx1]+=1;
881 if(idx2>=0 && dig>=zeroSup && idx2<nParticles) count[idx2]+=1;
882 } while (digits->Next());
883 }
884
885 for(Int_t j=0; j<nParticles; j++){
886
887 if(count[j]>1){
888 if(sec>=digp->GetNInnerSector())
889 if(row==nrow_up-1) good[j]|=0x4000;
890 if(sec>=digp->GetNInnerSector())
891 if(row==nrow_up-1-gap) good[j]|=0x1000;
892
893 if(sec>=digp->GetNInnerSector())
894 if(row==nrow_up-1-shift) good[j]|=0x2000;
895 if(sec>=digp->GetNInnerSector())
896 if(row==nrow_up-1-gap-shift) good[j]|=0x800;
897
898 good[j]++;
899 }
900 count[j] = 0;
901 } // end of loop over particles
902 } // end of loop over sectors_by_rows
903 delete [] count;
904 } // end of case 2
905 break;
906 } // end of switch
907
908
909
910
911 // ---------------- Select sample of MC particles that will be used for forming the efficiency ------------------
912
913 AliStack *stack = runLoader->Stack();
914
915 for(i=0; i<nParticles; i++){ // loop over stack
916
917 if ((good[i]&0x5000) != 0x5000)//SG!!!
918 if ((good[i]&0x2800) != 0x2800) continue;
919 if ((good[i]&0x7FF ) < good_number) continue;
920 // N TPC rows with digits >= 64 => common checks of cluster finder & tracker
921 // There are digits in rows (159&&139) || (149&&129)
922
923
924 TParticle *part = (TParticle*)stack->Particle(i);
925 if(part == 0x0){
926 cerr << "Cannot get particle "<< i << endl;
927 continue;
928 }
929
930 if(part->Pt()<0.100) continue; // skip particles with pt below 0.1 GeV/c
931 if(TMath::Abs(part->Pz()/part->Pt())>0.999) continue; // reject tracks outside the TPC acceptance, below 45 degrees
932
933 Double_t vx = part->Vx();
934 Double_t vy = part->Vy();
935 Double_t vz = part->Vz();
936 //if (TMath::Sqrt(vx*vx+vy*vy)>3.5) continue; // vertex cuts
937 //if (TMath::Abs(vz) > 50.) continue;
938
939 AliTrackReference *ref = new((*refs)[nt])AliTrackReference();
940
941 ref->SetLabel(i);
942 Int_t pdg = part->GetPdgCode();
943 ref->SetLength(pdg); //This will the particle's PDG !
944 ref->SetMomentum(0.,0.,0.);
945 ref->SetPosition(0.,0.,0.);
946
947 nt++;
948 } // end of loop over stack
949
950 //----------- Check whether there is also information at the entrance of the TPC
951
952 TTree *TR = runLoader->TreeTR();
953 TBranch *branch = TR->GetBranch("TrackReferences");
954 if(branch==0){
955 ::Error("GoodTPCTracks","No track references!");
956 delete runLoader;
957 return 5;
958 }
959
960 TClonesArray tpcdummy("AliTrackReference",1000), *tpcRefs = &tpcdummy;
961 branch->SetAddress(&tpcRefs);
962
963 for(Int_t r=0; r<(Int_t)TR->GetEntries(); r++){
964
965 //cerr<<r<<' '<<(Int_t)TR->GetEntries()<<'\r';
966 TR->GetEvent(r);
967
968 if(!tpcRefs->GetEntriesFast()) continue;
969
970 AliTrackReference *tpcRef = 0x0;
971 for(Int_t iref=0; iref<tpcRefs->GetEntriesFast(); ++iref){
972
973 tpcRef = (AliTrackReference*)tpcRefs->UncheckedAt(iref);
974 if(tpcRef->DetectorId() == AliTrackReference::kTPC) break;
975 tpcRef = 0x0;
976 }
977
978 if(!tpcRef) continue;
979
980 Int_t j;
981 AliTrackReference *ref = 0;
982
983 for(j=0; j<nt; j++){
984
985 ref = (AliTrackReference*)refs->UncheckedAt(j);
986 if(ref->Label()==tpcRef->Label()) break;
987 ref = 0;
988 }
989 if(ref==0) continue;
990
991 ref->SetMomentum(tpcRef->Px(),tpcRef->Py(),tpcRef->Pz());
992 ref->SetPosition(tpcRef->X(), tpcRef->Y(), tpcRef->Z());
993
994 tpcRefs->Clear();
995 }
996
997 tpcTree.Fill();
998
999 delete [] good;
1000 } //---------- end of the loop over generated events
1001
1002 tpcTree.Write();
1003 file->Close();
1004
1005 delete runLoader;
1006 return 0;
1007}
1008#endif