3 * Modified version of $ALICE_ROOT/TPC/AliTPCComparison.C for evaluating the
4 * performance of the CA tracker.
8 * aliroot AliHLTTPCTrackerEvaluation.C
11 * documentation to be filled in when the macro is finalized
16 #if defined(__CINT__) && !defined(__MAKECINT__)
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)
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++");
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);
33 #include <Riostream.h>
38 #include <TParticle.h>
43 #include <TBenchmark.h>
49 #include "AliHeader.h"
50 #include "AliTrackReference.h"
51 #include "AliRunLoader.h"
53 #include "AliESDEvent.h"
54 #include "AliESDtrack.h"
55 #include "AliTPCtrack.h"
57 #include "AliSimDigits.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"
67 #include "AliCDBManager.h"
68 #include "AliTPCcalibDB.h"
70 //_______________________________________________________________________________________
73 Int_t GoodTPCTracks(const Char_t *dir = ".");
75 extern AliRun *gAlice;
76 extern TBenchmark *gBenchmark;
79 static Int_t allgood = 0;
80 static Int_t allselected = 0;
81 static Int_t allfound = 0;
85 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){
88 cerr << "please specify an input file" << endl;
92 gBenchmark->Start("AliHLTTPCTrackerEvaluation");
94 ::Info("AliHLTTPCTrackerEvaluation.C","Calculating reconstruction efficiency...");
96 gSystem->SetIncludePath("-I$ROOTSYS/include -I$ALICE_ROOT/include -I$ALICE_ROOT -I$ALICE/geant3/TGeant3");
98 // --------------- Efficiency histograms ------------------------
100 TH1F *hClus = new TH1F("hClus","# clusters on reconstructed tracks",300,0,300);
102 TH1F *hGood = new TH1F("hGood", "P_{t} distribution of MC tracks (selected for efficiency calculation)", 30, 0, 6);
103 TH1F *hFake = new TH1F("hFake", "P_{t} distribution of fake tracks", 30, 0, 6);
104 TH1F *hFound = new TH1F("hFound", "P_{t} distribution of reconstructed tracks", 30, 0, 6);
105 TH1F *hFoundMult = new TH1F("hFoundMult","P_{t} distribution of reconstructed tracks with clones", 30, 0, 6);
106 TH1F *hClone = new TH1F("hClone", "P_{t} distribution of clone tracks", 30, 0, 6);
108 TH1F *hEff = new TH1F("hEff", "Reconstruction efficiency based on selected MC tracks", 30, 0, 6);
109 TH1F *hFakeEff = new TH1F("hFakeEff", "Efficiency for fake tracks", 30, 0, 6);
110 TH1F *hCloneEff = new TH1F("hCloneEff","Efficiency for clone tracks", 30, 0, 6);
112 hEff->SetLineWidth(2);
113 hEff->SetMaximum(1.4);
114 hEff->SetYTitle("efficiency");
115 hEff->GetYaxis()->CenterTitle();
116 hEff->SetXTitle("P_{t} (GeV/c)");
118 hFakeEff->SetLineColor(kRed);
119 hFakeEff->SetLineWidth(2);
121 hCloneEff->SetLineColor(kBlue);
122 hCloneEff->SetLineWidth(2);
124 TH2F *hDedx = new TH2F("hDedx","dE/dx vs. momentum", 50,0.,2.,50,0.,400.);
125 hDedx->SetMarkerStyle(8);
126 hDedx->SetMarkerSize(0.4);
127 hDedx->SetXTitle("P (Gev/c)");
128 hDedx->SetYTitle("dE/dx (a.u.)");
130 TH1F *hnhit_ref = new TH1F("hnhit_ref","# clusters on reco tracks, used in fit performance",300,0,300);
131 TH1F *he = new TH1F("he","dE/dx for pions with 0.4<p<0.5 GeV/c",50,0.,100.);
133 // --------------- Resolution histograms ---------------------------
135 TH2F *hPt = new TH2F("hPt", "", 30, 0, 6, 50, -3, 3);
136 TH2F *hPhi = new TH2F("hPhi", "", 30, 0, 6, 50, -16, 16);
137 TH2F *hLambda = new TH2F("hLambda", "", 30, 0, 6, 50, -10, 10);
138 TH2F *hY = new TH2F("hY", "", 30, 0, 6, 30, -5, 5);
139 TH2F *hZ = new TH2F("hZ", "", 30, 0, 6, 30, -5, 5); // fix the x axis, if this is z dependent
141 TH1F *hResPt = new TH1F("hResPt", "P_{t} resolution", 30, 0, 6);
142 TH1F *hResPhi = new TH1F("hResPhi", "#phi resolution", 30, 0, 6);
143 TH1F *hResLambda = new TH1F("hResLambda", "#lambda resolution",30, 0, 6);
144 TH1F *hResY = new TH1F("hResY", "Y resolution", 30, 0, 6);
145 //TH1F *hResZ = new TH1F("hResZ", "Z resolution", 30, 0, 6); // select appropriate x axis, KKK
147 hResPt->SetXTitle("P_{t} (GeV/c)");
148 hResPt->SetYTitle("#sigma_{(P_{t}-P_{MC})/P_{MC}} (%)");
149 hResPt->GetYaxis()->CenterTitle(true);
151 hResPhi->SetXTitle("P_{t} (GeV/c)");
152 hResPhi->SetYTitle("#sigma_{(#phi_{rec}-#phi_{MC})/#phi_{MC}} (%)");
153 hResPhi->GetYaxis()->CenterTitle(true);
155 hResLambda->SetXTitle("P_{t} (GeV/c)");
156 hResLambda->SetYTitle("#sigma_{(#lambda_{rec}-#lambda_{MC})/#lambda_{MC}} (%)"); // KKK perhaps we should add the definition on the histo pad
157 hResLambda->GetYaxis()->CenterTitle(true);
159 hResY->SetXTitle("P_{t} (GeV/c)");
160 hResY->SetYTitle("#sigma_{(Y_{rec}-Y_{MC})/Y_{MC}} (%)");
161 hResY->GetYaxis()->CenterTitle(true);
163 // hResZ->SetXTitle("Z (mm)");
164 // hResZ->SetYTitle("#sigma_{(Z_{rec}-Z_{MC})/Z_{MC}} (%)");
165 // hResZ->GetYaxis()->CenterTitle(true);
167 // KKK I leave the Z coordinate up to you. You can make the hResZ Z dependent.
170 TH1D *hProjPt[6], *hProjPhi[6], *hProjLambda[6], *hProjY[6], *hProjZ[6]; // I am not sure we need these histos to be pointers
171 Char_t name[15]; Char_t title[100];
173 for(Int_t i=0; i<6; i++){
174 sprintf(name,"hProjPt%i",i+1);
175 sprintf(title,"(Pt_{MC}-Pt_{Rec})/Pt_{MC} @ Pt#in[%f, %f] GeV/c", i-0.5, i+0.5);
176 hProjPt[i] = new TH1D(name, title, 50, -3, -3);
178 sprintf(name,"hProjPhi%i",i+1);
179 sprintf(title,"(#phi_{MC}-#phi_{Rec})/#phi_{MC} @ #phi#in[%f, %f] GeV/c", i-0.5, i+0.5);
180 hProjPhi[i] = new TH1D(name, title, 50, -16, -16);
182 sprintf(name,"hProjLambda%i",i+1);
183 sprintf(title,"(#lambda_{MC}-#lambda_{Rec})/#lambda_{MC} @ #lambda#in[%f, %f] GeV/c", i-0.5, i+0.5);
184 hProjLambda[i] = new TH1D(name, title, 50, -16, -16);
186 sprintf(name,"hProjY%i",i+1);
187 sprintf(title,"(Y_{MC}-Y_{Rec})/Y_{MC} @ Y#in[%f, %f] GeV/c", i-0.5, i+0.5);
188 hProjY[i] = new TH1D(name, title, 50, -16, -16);
190 // here you add the hProjZ[] histograms KKK
195 // --------------- Pull variable histograms ------------------------
197 TH1F *hpullPhi = new TH1F("hpullPhi", "SinPhi pull", 30,-10.,10.);
198 TH1F *hpullY = new TH1F("hpullY", "Y pull", 30,-10.,10.);
199 TH1F *hpullZ = new TH1F("hpullZ", "Z pull", 30,-10.,10.);
200 TH1F *hpullDzds = new TH1F("hpullDzds","Dzds pull", 30,-10.,10.);
201 TH1F *hpullK = new TH1F("hpullK", "Kappa pull", 30,-10.,10.);
203 //---------------------------------------------------------------------------------------
206 sprintf(fname,"%s/GoodTPCTracks.root",dir);
208 TFile *refFile = TFile::Open(fname,"old");
209 if(!refFile || !refFile->IsOpen()){
210 ::Info("AliHLTTPCTrackerEvaluation.C","Marking good tracks (will take a while)...");
211 if(GoodTPCTracks(dir)){
212 ::Error("AliHLTTPCTrackerEvaluation.C","Cannot generate the reference file!");
217 refFile = TFile::Open(fname,"old");
218 if(!refFile || !refFile->IsOpen()){
219 ::Error("AliHLTTPCTrackerEvaluation.C","Cannot open the reference file!");
223 //------------ access the contents of the tree created by the function GoodTracks() -------------
225 TTree *tpcTree = (TTree*)refFile->Get("tpcTree");
227 ::Error("AliHLTTPCTrackerEvaluation.C","Cannot get the reference tree!");
231 TBranch *branch = tpcTree->GetBranch("TPC");
233 ::Error("AliHLTTPCTrackerEvaluation.C","Cannot get the TPC branch!");
237 TClonesArray dummy("AliTrackReference",1000), *refs = &dummy;
238 branch->SetAddress(&refs);
240 //------------ access the HLT output --------------------
242 sprintf(fname,"%s/%s", dir, input);
243 TFile *ef = TFile::Open(fname);
245 if( !ef || !ef->IsOpen() ){
246 ::Error("AliHLTTPCTrackerEvaluation.C","Cannot open [%s]!",fname);
250 AliESDEvent *event = new AliESDEvent();
251 TTree *esdTree = NULL;
252 TString tsinput = input;
254 if(tsinput.CompareTo("AliESDs.root")==0){
255 esdTree = (TTree*)ef->Get("HLTesdTree");
257 ::Error("AliHLTTPCTrackerEvaluation.C", "no HLTESD tree found");
260 event->ReadFromTree(esdTree);
261 } else if(tsinput.CompareTo("AliHLTESDs.root")==0){
262 esdTree = (TTree*)ef->Get("esdTree");
264 ::Error("AliHLTTPCTrackerEvaluation.C", "no ESD tree found");
267 event->ReadFromTree(esdTree);
271 //---------- Loop over events -------------------
274 while(esdTree->GetEvent(iEvent)){
276 Int_t nEntries = event->GetNumberOfTracks();
277 //cout << "********* Processing event number: " << iEvent << ", " << nEntries << " reconstructed tracks *******\n" << endl;
279 allfound += nEntries;
281 if(tpcTree->GetEvent(iEvent++)==0){
282 cerr << "No reconstructable tracks !\n";
286 Int_t ngood = refs->GetEntriesFast(); // access the track references
287 //cout << ngood << " good MC tracks" << endl;
290 const Int_t MAX = 15000;
292 Int_t track_notfound[MAX], itrack_notfound = 0;
293 Int_t track_fake[MAX], itrack_fake = 0;
294 Int_t track_multifound[MAX], track_multifound_n[MAX], itrack_multifound = 0;
296 for(Int_t i=0; i<nEntries; i++){
298 hClus->Fill(event->GetTrack(i)->GetTPCNcls()); // filled but not displayed
299 //cout << "TPC label for track " << i << " = " << event->GetTrack(i)->GetTPCLabel() << ", # clusters " << event->GetTrack(i)->GetTPCNcls() << endl;
303 for(Int_t k=0; k<ngood; k++){ // read the MC information
305 AliTrackReference *ref = (AliTrackReference*)refs->UncheckedAt(k);
307 Int_t label = ref->Label();
309 Float_t ptMC = TMath::Sqrt(ref->Px()*ref->Px() + ref->Py()*ref->Py());
310 Float_t pMC = TMath::Sqrt(ref->Px()*ref->Px() + ref->Py()*ref->Py()+ref->Pz()*ref->Pz());
312 if (ptMC<1e-33 || ptMC<ptLowerCut || ptMC>ptUpperCut) continue; // first condition is for tracks not crossing padrow 0
317 for(i=0; i<nEntries; i++){
318 Int_t ttlabel = event->GetTrack(i)->GetTPCLabel();
319 if( label==TMath::Abs(ttlabel) && ttlabel<0 ){
320 track_fake[itrack_fake++] = label;
325 AliESDtrack *esdtrack = 0;
326 for(i=0; i<nEntries; i++){
327 esdtrack = event->GetTrack(i);
328 tlabel = esdtrack->GetTPCLabel();
329 if(label==tlabel) break;
333 track_notfound[itrack_notfound++] = label;
334 //cout << "MC track " << label << " not reconstructed" << endl;
340 AliESDtrack *mitrack;
342 for(mi=0; mi<nEntries; mi++){
343 mitrack = event->GetTrack(mi);
344 if(label==mitrack->GetTPCLabel()) micount++;
348 track_multifound[itrack_multifound] = label;
349 track_multifound_n[itrack_multifound] = micount;
351 hClone->Fill(ptMC,micount-1);
354 //if((mitrack->GetStatus()&AliESDtrack::kTPCrefit)==0) continue;
355 //if((mitrack->GetStatus()&AliESDtrack::kTPCin)==0) continue;
359 //if( micount<1 ) cout<<"ERROR!!!!"<<endl;
360 hFoundMult->Fill(ptMC,micount);
363 AliExternalTrackParam tpctrack =*(esdtrack->GetInnerParam());
365 if( TMath::Abs(tpctrack.GetSigned1Pt()) <1.e-8 ) continue;
367 Double_t mcxyz[3] = { ref->X(), ref->Y(), ref->Z() };
368 Double_t mclocal[3] = { ref->X(), ref->Y(), ref->Z() };
371 Double_t c = TMath::Cos(tpctrack.GetAlpha());
372 Double_t s = TMath::Sin(tpctrack.GetAlpha());
373 Double_t x = mclocal[0];
374 mclocal[0] = x*c + mclocal[1]*s;
375 mclocal[1] =-x*s + mclocal[1]*c;
379 Double_t local[3] = { tpctrack.GetX(),tpctrack.GetY(),tpctrack.GetZ() };
380 // cout << "label: " << label << endl;
381 // cout << "Z diff = " << local[2] - mclocal[2] << endl;
382 // cout << "orig rec: " << local[0]<<" "<<local[1]<<" "<<local[2] <<" "<<TMath::Sqrt(local[0]*local[0]+local[1]*local[1])<<endl;
383 // cout << "rotated mc: " << mclocal[0]<<" "<<mclocal[1]<<" "<<mclocal[2] <<" "<<TMath::Sqrt(mclocal[0]*mclocal[0]+mclocal[1]*mclocal[1])<<endl;
384 // cout << "orig mc: " << mcxyz[0]<<" "<<mcxyz[1]<<" "<<mcxyz[2]<<endl;
385 // cout << "alpha = " << tpctrack.GetAlpha()/TMath::Pi()*180.<<" "<<ref->Alpha()/TMath::Pi()*180.<<endl;
388 //cout << "X = " << mclocal[0] << " / " << tpctrack.GetX() << endl;
389 tpctrack.AliExternalTrackParam::PropagateTo(mclocal[0],5.);
391 Double_t xyz[3], pxpypz[3];
393 tpctrack.GetXYZ(xyz); // GetInnerXYZ
394 tpctrack.GetPxPyPz(pxpypz); // GetInnerPxPyPz
396 Double_t local[3] = { tpctrack.GetX(),tpctrack.GetY(),tpctrack.GetZ() };
397 Float_t phiRec = TMath::ATan2(pxpypz[1],pxpypz[0]);
399 if(phiRec <= TMath::Pi()) phiRec += 2*TMath::Pi();
400 if(phiRec >= TMath::Pi()) phiRec -= 2*TMath::Pi();
402 Double_t ptRec = TMath::Sqrt(pxpypz[0]*pxpypz[0]+pxpypz[1]*pxpypz[1]);
403 Double_t pRec = TMath::Sqrt(pxpypz[0]*pxpypz[0]+pxpypz[1]*pxpypz[1]+pxpypz[2]*pxpypz[2]);
405 Float_t lambdaRec = TMath::ATan2(pxpypz[2],ptRec);
406 Float_t pt_1 = 1./ptRec;
407 //Float_t pts_1 = tpctrack->GetSigned1Pt();
408 //Float_t ptMCs_1 = ref->/ptMC;
412 Int_t pdg = (Int_t)ref->GetLength(); // particle PDG
413 const Double_t kCLight = 0.000299792458;
417 if( TMath::Abs(pdg)<10000 ){
418 TParticlePDG *ppdg = TDatabasePDG::Instance()->GetParticle(pdg);
419 if( ppdg ) q = ppdg->Charge()/3.;
422 Double_t kappaMC = q/ptMC; // /Bz/kCLight;
423 hnhit_ref->Fill(esdtrack->GetTPCNcls());
425 hPt->Fill(ptMC, (ptRec - ptMC)/(ptMC)*100.);
427 Float_t phiMC = TMath::ATan2(ref->Py(),ref->Px());
428 hPhi->Fill(ptMC, (phiRec - phiMC)/phiMC*1000.);
430 Float_t lambdaMC = TMath::ATan2(ref->Pz(),ptMC);
431 hLambda->Fill(ptMC, (lambdaRec - lambdaMC)/lambdaMC*1000.);
433 hY->Fill(ptMC, (local[1]-mclocal[1])*10.);
434 /*if(mclocal[2]>0 )*/ //hZ->Fill( (local[2]-mclocal[2] ) *10.); // Please check hY and hZ! KKK
436 if( tpctrack.GetSigmaY2()>0 && finite(tpctrack.GetSigmaY2()) ) hpullY ->Fill( (local[1]-mclocal[1])/TMath::Sqrt(TMath::Abs(tpctrack.GetSigmaY2())) );
437 if( tpctrack.GetSigmaZ2()>0 && finite(tpctrack.GetSigmaZ2()) ) hpullZ ->Fill( (local[2]-mclocal[2])/TMath::Sqrt(TMath::Abs(tpctrack.GetSigmaZ2())) );
438 if( tpctrack.GetSigmaSnp2()>0 && finite(tpctrack.GetSigmaSnp2()) ) hpullPhi ->Fill( (phiRec-phiMC)/TMath::Sqrt(TMath::Abs(tpctrack.GetSigmaSnp2())) );
439 if( tpctrack.GetSigmaTgl2()>0 && finite(tpctrack.GetSigmaTgl2()) ) hpullDzds->Fill( (lambdaRec-lambdaMC)/TMath::Sqrt(TMath::Abs(tpctrack.GetSigmaTgl2())) );
441 if( tpctrack.GetSigma1Pt2()>0 && finite(tpctrack.GetSigma1Pt2()) )
442 if( q!=0 ) hpullK->Fill((tpctrack.GetSigned1Pt()-kappaMC)/TMath::Sqrt(TMath::Abs(tpctrack.GetSigma1Pt2())));
444 Float_t dedx = esdtrack->GetTPCsignal();
445 hDedx->Fill(pRec,dedx,1.);
447 if(TMath::Abs(pdg)==211) // pions
448 if(pRec>0.4 && pRec<0.5){
451 } // end of loop over "good" tracks
454 cout<<"\nList of Not found tracks :\n";
455 for ( i = 0; i< itrack_notfound; i++){
456 cout<<track_notfound[i]<<"\t";
457 if ((i%5)==4) cout<<"\n";
459 cout<<"\nList of fake tracks :\n";
460 for ( i = 0; i< itrack_fake; i++){
461 cout<<track_fake[i]<<"\t";
462 if ((i%5)==4) cout<<"\n";
464 cout<<"\nList of multiple found tracks :\n";
465 for ( i=0; i<itrack_multifound; i++) {
466 cout<<"id. "<<track_multifound[i]
467 <<" found - "<<track_multifound_n[i]<<"times\n";
470 cout<<"Number of found tracks ="<<nEntries<<endl;
471 cout<<"Number of \"good\" tracks ="<<ngood<<endl;
475 } // end of the loop over events
484 Stat_t nGoodMC = hGood->GetEntries();
485 Stat_t nFakes = hFake->GetEntries();
486 Stat_t nClones = hClone->GetEntries();
487 Stat_t nRec = hFound->GetEntries();
488 Stat_t nRecMult = hFoundMult->GetEntries();
491 if(nGoodMC!=0) ::Info("\n\nAliHLTTPCTrackerEvaluation","Integral efficiency in (%) is about : %f\n", nRec/nGoodMC*100.);
492 if(nRecMult+nFakes!=0) ::Info("AliHLTTPCTrackerEvaluation","Integral fake rate in (%) is about : %f\n", nFakes/(nRecMult+nFakes)*100.);
493 if(nRecMult+nFakes!=0) ::Info("AliHLTTPCTrackerEvaluation","Integral clone rate in (%) is about : %f\n", nClones/(nRecMult+nFakes)*100.);
495 ::Info("AliHLTTPCTrackerEvaluation", "Total number of selected MC tracks : %d\n", allgood);
496 ::Info("AliHLTTPCTrackerEvaluation", "Total number of selected MC tracks with momentum cuts : %d\n", allselected);
497 ::Info("AliHLTTPCTrackerEvaluation", "Total number of reconstructed tracks : %d\n", allfound);
500 // ---------- Plotting the histograms ----------------------------------------------------------------------------
508 hEff ->Divide(hFound,hGood, 1,1.,"b");
509 hFakeEff ->Divide(hFake, hGood, 1,1.,"b");
510 hCloneEff->Divide(hClone,hFoundMult,1,1.,"b");
511 // KKK are you sure this is the correct definition? When plotted, it looks like we have 70% clones in the lower momenta...
512 // perhaps fakes and clones should be displayed in another histogram??
514 gROOT->SetStyle("Plain");
515 gStyle->SetOptStat(0); //gStyle->SetOptStat(111110);
516 gStyle->SetOptFit(0); //gStyle->SetOptFit(1);
517 gStyle->SetPalette(1);
519 //--------------- canvas 0 for EFFICIENCY -------------------
521 TCanvas *c0 = new TCanvas("c0","reconstruction efficiency",500,450);
525 hCloneEff->Draw("histsame");
526 hFakeEff ->Draw("histsame");
528 TLegend *leg = new TLegend(0.3,0.5,0.85,0.8);
529 leg->SetFillColor(10);
530 leg->SetLineColor(10);
531 leg->AddEntry(hEff, "reconstruction eff.", "l");
532 leg->AddEntry(hCloneEff, "clone contribution", "l");
533 leg->AddEntry(hFakeEff, "fake contribution", "l");
536 TLine *line1 = new TLine(0,1,6,1);
537 line1->SetLineStyle(4);
540 Float_t fakeData = nFakes/(nRecMult+nFakes)*100.;
541 Float_t cloneData = nClones/(nRecMult+nFakes)*100.;
543 Char_t cloneBuf[100];
545 sprintf(fakeBuf,"%d", (Int_t)fakeData);
546 TString fakeStr = fakeBuf;
547 fakeStr.Append(" %");
549 sprintf(cloneBuf,"%d", (Int_t)cloneData);
550 TString cloneStr = cloneBuf;
551 cloneStr.Append(" %");
553 TText *text = new TText(1.4,0.1,fakeStr.Data());
554 text->SetTextSize(0.05);
555 text->SetTextColor(kRed);
557 text = new TText(0.6,0.3,cloneStr.Data());
558 text->SetTextSize(0.05);
559 text->SetTextColor(kBlue);
563 //--------------- canvas 1 for RESOLUTIONS ----------------
566 TCanvas *c1 = new TCanvas("c1","resolutions",1100,900); // please add the Y and Z projections
568 TF1 *fGauss = new TF1("gauss","gaus");
569 for(Int_t i=0; i<6; i++){
570 Float_t pMin = i - 0.5;
571 Float_t pMax = i + 0.5;
572 Int_t binx1 = hPt->GetXaxis()->FindBin(pMin);
573 Int_t binx2 = hPt->GetXaxis()->FindBin(pMax);
576 *hProjPt[i] = *(hPt ->ProjectionY("", binx1, binx2));
577 *hProjPhi[i] = *(hPhi ->ProjectionY("", binx1, binx2));
578 *hProjLambda[i] = *(hLambda->ProjectionY("", binx1, binx2));
579 *hProjY[i] = *(hY ->ProjectionY("", binx1, binx2));
580 // *hProjZ[i] = *(hZ ->ProjectionY("", binx1, binx2));
582 if(hProjPt[i]->GetEntries()>30) hProjPt[i]->Fit("gauss","Q");
583 hResPt->Fill(i+0.5, fGauss->GetParameter(2));fGauss->GetParError(2);
584 fGauss->SetParameter(2,0); fGauss->SetParError(2,0);
585 //KKK I am resetting only the sigma and its error, perhaps we should do this with all the fit parameters ???
587 if(hProjPhi[i]->GetEntries()>30) hProjPhi[i]->Fit("gauss","Q");
588 hResPhi->Fill(i+0.5, fGauss->GetParameter(2)); fGauss->GetParError(2);
589 fGauss->SetParameter(2,0); fGauss->SetParError(2,0);
591 if(hProjLambda[i]->GetEntries()>30) hProjLambda[i]->Fit("gauss","Q");
592 hResLambda->Fill(i+0.5, fGauss->GetParameter(2)); fGauss->GetParError(2);
593 fGauss->SetParameter(2,0);
595 if(hProjY[i]->GetEntries()>30) hProjY[i]->Fit("gauss","Q");
596 hResY->Fill(i+0.5, fGauss->GetParameter(2)); fGauss->GetParError(2);
597 fGauss->SetParameter(2,0); fGauss->SetParError(2,0);
599 // if(hProjZ[i]->GetEntries()>30) hProjZ[i]->Fit("gauss","Q");
600 // hResZ->Fill(i+0.5, fGauss->GetParameter(2)); fGaus->GetParError(2);
601 // fGauss->SetParameter(2,0); KKK
605 c1->Clear(); c1->Divide(3,2);
608 hResPt->Draw(); // KKK I haven't filled the errors for the sigma, perhaps a TGraph would make them easier to plot
620 // hResZ->Draw(); // KKK
625 //----------------- optional canvases 2 and 3 for dE/dx and clusters ----------------
628 TCanvas *c2 = new TCanvas("c2","dE/dx",500,450);
632 TCanvas *c3 = new TCanvas("c","clusters",500,450);
635 hnhit_ref->SetLineColor(kRed);
636 hnhit_ref->Draw("same");
640 //----------------- canvas 3 for PULL VARIABLES --------------------
642 TCanvas *c4 = new TCanvas("c4","pull variables",1100,900);
647 hpullY->Fit("gaus","Q");
651 hpullZ->Fit("gaus","Q");
655 hpullPhi->Fit("gaus","Q");
659 hpullDzds->Fit("gaus","Q");
663 hpullK->Fit("gaus","Q");
667 //-------------------- OUTPUT FILE ----------------------------------------
669 //TFile fc("CAtrackerEvaluation.root","RECREATE");
673 //fc.Close(); // KKK commented out temporarily
675 gBenchmark->Stop("AliHLTTPCTrackerEvaluation");
676 gBenchmark->Show("AliHLTTPCTrackerEvaluation");
682 //============================================================================
684 Int_t GoodTPCTracks(const Char_t *dir){
687 sprintf(fname,"%s/galice.root",dir);
689 AliRunLoader *runLoader = AliRunLoader::Open(fname,"COMPARISON");
691 ::Error("GoodTPCTracks","Cannot start session!");
695 // load MC information
697 runLoader->LoadgAlice();
698 runLoader->LoadHeader();
699 runLoader->LoadKinematics();
700 runLoader->LoadTrackRefs();
702 AliTPCLoader *tpcLoader = (AliTPCLoader*)runLoader->GetLoader("TPCLoader");
703 if (tpcLoader == 0x0) {
704 ::Error("GoodTPCTracks","Cannot find TPCLoader!");
709 AliTPC *TPC = (AliTPC*)runLoader->GetAliRun()->GetDetector("TPC");
710 Int_t tpcVersion = TPC->IsVersion();
712 //cout << "TPC version "<< tpcVersion << " has been found!\n";
714 if (tpcVersion==1) tpcLoader->LoadRecPoints();
715 else if (tpcVersion==2) tpcLoader->LoadDigits();
717 ::Error("GoodTPCTracks","Wrong TPC version!");
722 AliCDBManager *manager = AliCDBManager::Instance();
723 manager->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
726 // loading the TPC parameters
727 AliTPCParamSR *digp = (AliTPCParamSR*)(AliTPCcalibDB::Instance()->GetParameters());
729 ::Error("AliHLTTPCTrackerEvaluation.C","TPC parameters have not been found!");
734 Int_t nrow_up = digp->GetNRowUp(); //get the number of pad rows in up sector
735 Int_t nrows = digp->GetNRowLow()+nrow_up; //get the number of pad rows in low sector and add
736 Int_t zeroSup = digp->GetZeroSup();
737 Int_t gap = Int_t(0.125*nrows), shift = Int_t(0.5*gap);
738 Int_t good_number = Int_t(0.4*nrows); // will be used for selecting tracks with hits in 40% of the rows
740 Int_t nEvents = runLoader->GetNumberOfEvents();
741 ::Info("GoodTPCTracks","Number of events : %d\n", nEvents);
743 sprintf(fname,"%s/GoodTPCTracks.root",dir);
744 TFile *file = TFile::Open(fname,"recreate");
746 TClonesArray dummy("AliTrackReference",1000), *refs = &dummy;
747 TTree tpcTree("tpcTree","Tree with information about the reconstructable TPC tracks");
748 tpcTree.Branch("TPC",&refs);
752 //-------------- Loop over generated events ------------------------------
754 for(Int_t iEvent=0; iEvent<nEvents; iEvent++){
761 runLoader->GetEvent(iEvent);
764 Int_t nParticles = runLoader->GetHeader()->GetNtrack();
765 //cout << "Event " << iEvent << ", Number of particles: " << nParticles << endl;
767 // ------- Set the selection criteria for the MC sample that will be used for the efficiency -------
769 Int_t *good = new Int_t[nParticles]; for(i=0; i<nParticles; i++) good[i] = 0;
772 case 1: // not used any more
774 AliTPCClustersArray *pca = new AliTPCClustersArray, &ca=*pca;
776 ca.SetClusterType("AliTPCcluster");
777 ca.ConnectTree(tpcLoader->TreeR());
778 Int_t nrows=Int_t(ca.GetTree()->GetEntries());
779 for (Int_t n=0; n<nrows; n++) {
780 AliSegmentID *s=ca.LoadEntry(n);
782 digp->AdjustSectorRow(s->GetID(),sec,row);
783 AliTPCClustersRow &clrow = *ca.GetRow(sec,row);
784 Int_t ncl=clrow.GetArray()->GetEntriesFast();
786 AliTPCcluster *c=(AliTPCcluster*)clrow[ncl];
787 Int_t label=c->GetLabel(0);
788 if (label<0) continue; //noise cluster
789 label=TMath::Abs(label);
791 if (sec>=digp->GetNInnerSector())
792 if (row==nrow_up-1) good[label]|=0x4000;
793 if (sec>=digp->GetNInnerSector())
794 if (row==nrow_up-1-gap) good[label]|=0x1000;
796 if (sec>=digp->GetNInnerSector())
797 if (row==nrow_up-1-shift) good[label]|=0x2000;
798 if (sec>=digp->GetNInnerSector())
799 if (row==nrow_up-1-gap-shift) good[label]|=0x800;
803 ca.ClearRow(sec,row);
811 TTree *TD = tpcLoader->TreeD(); // get the digits tree
813 AliSimDigits da, *digits = &da;
814 TD->GetBranch("Segment")->SetAddress(&digits);
816 Int_t *count = new Int_t[nParticles];
818 for (i=0; i<nParticles; i++) count[i] = 0;
820 Int_t sectors_by_rows = (Int_t)TD->GetEntries();
822 for(i=0; i<sectors_by_rows; i++){
824 if (!TD->GetEvent(i)) continue;
827 digp->AdjustSectorRow(digits->GetID(),sec,row);
830 Int_t it = digits->CurrentRow();
831 Int_t ip = digits->CurrentColumn();
832 Short_t dig = digits->GetDigit(it,ip);
833 Int_t idx0 = digits->GetTrackID(it,ip,0);
834 Int_t idx1 = digits->GetTrackID(it,ip,1);
835 Int_t idx2 = digits->GetTrackID(it,ip,2);
836 if(idx0>=0 && dig>=zeroSup && idx0<nParticles) count[idx0]+=1;
837 if(idx1>=0 && dig>=zeroSup && idx1<nParticles) count[idx1]+=1;
838 if(idx2>=0 && dig>=zeroSup && idx2<nParticles) count[idx2]+=1;
839 } while (digits->Next());
842 for(Int_t j=0; j<nParticles; j++){
845 if(sec>=digp->GetNInnerSector())
846 if(row==nrow_up-1) good[j]|=0x4000;
847 if(sec>=digp->GetNInnerSector())
848 if(row==nrow_up-1-gap) good[j]|=0x1000;
850 if(sec>=digp->GetNInnerSector())
851 if(row==nrow_up-1-shift) good[j]|=0x2000;
852 if(sec>=digp->GetNInnerSector())
853 if(row==nrow_up-1-gap-shift) good[j]|=0x800;
858 } // end of loop over particles
859 } // end of loop over sectors_by_rows
868 // ---------------- Select sample of MC particles that will be used for forming the efficiency ------------------
870 AliStack *stack = runLoader->Stack();
872 for(i=0; i<nParticles; i++){ // loop over stack
874 if ((good[i]&0x5000) != 0x5000)//SG!!!
875 if ((good[i]&0x2800) != 0x2800) continue;
876 if ((good[i]&0x7FF ) < good_number) continue;
877 // N TPC rows with digits >= 64 => common checks of cluster finder & tracker
878 // There are digits in rows (159&&139) || (149&&129)
881 TParticle *part = (TParticle*)stack->Particle(i);
883 cerr << "Cannot get particle "<< i << endl;
887 if(part->Pt()<0.100) continue; // skip particles with pt below 0.1 GeV/c
888 if(TMath::Abs(part->Pz()/part->Pt())>0.999) continue; // reject tracks outside the TPC acceptance, below 45 degrees
890 Double_t vx = part->Vx();
891 Double_t vy = part->Vy();
892 Double_t vz = part->Vz();
893 //if (TMath::Sqrt(vx*vx+vy*vy)>3.5) continue; // vertex cuts
894 //if (TMath::Abs(vz) > 50.) continue;
896 AliTrackReference *ref = new((*refs)[nt])AliTrackReference();
899 Int_t pdg = part->GetPdgCode();
900 ref->SetLength(pdg); //This will the particle's PDG !
901 ref->SetMomentum(0.,0.,0.);
902 ref->SetPosition(0.,0.,0.);
905 } // end of loop over stack
907 //----------- Check whether there is also information at the entrance of the TPC
909 TTree *TR = runLoader->TreeTR();
910 TBranch *branch = TR->GetBranch("TrackReferences");
912 ::Error("GoodTPCTracks","No track references!");
917 TClonesArray tpcdummy("AliTrackReference",1000), *tpcRefs = &tpcdummy;
918 branch->SetAddress(&tpcRefs);
920 for(Int_t r=0; r<(Int_t)TR->GetEntries(); r++){
922 //cerr<<r<<' '<<(Int_t)TR->GetEntries()<<'\r';
925 if(!tpcRefs->GetEntriesFast()) continue;
927 AliTrackReference *tpcRef = 0x0;
928 for(Int_t iref=0; iref<tpcRefs->GetEntriesFast(); ++iref){
930 tpcRef = (AliTrackReference*)tpcRefs->UncheckedAt(iref);
931 if(tpcRef->DetectorId() == AliTrackReference::kTPC) break;
935 if(!tpcRef) continue;
938 AliTrackReference *ref = 0;
942 ref = (AliTrackReference*)refs->UncheckedAt(j);
943 if(ref->Label()==tpcRef->Label()) break;
948 ref->SetMomentum(tpcRef->Px(),tpcRef->Py(),tpcRef->Pz());
949 ref->SetPosition(tpcRef->X(), tpcRef->Y(), tpcRef->Z());
957 } //---------- end of the loop over generated events