ITS/UPGRADE/macros/efficiency.C

   1
   2 #if !defined(__CINT__) || defined(__MAKECINT__)
   3 #include <TH1F.h>
   4 #include <TH1I.h>
   5 #include <TH2D.h>
   6 #include <TFile.h>
   7 #include <Riostream.h>
   8 #include <TParticle.h>
   9 #include <TCanvas.h>
  10 #include <TLegend.h>
  11 #include <TSystem.h>
  12 #include <TROOT.h>
  13 #include <TClonesArray.h>
  14 #include <TTree.h>
  15 #include <TArrayF.h>
  16 #include <TStyle.h>
  17 #include "AliRun.h"
  18 #include "AliRunLoader.h"
  19 #include "AliStack.h"
  20 #include "AliESDEvent.h"
  21 #include "AliESDtrack.h"
  22 #include "AliTrackReference.h"
  23 #include "AliITSsegmentationUpgrade.h"
  24 #endif
  25
  26 Int_t minNpoints=3;
  27 Bool_t IsTrackable(TClonesArray *trackRef);
  28 TArrayF radii;
  29
  30 void efficiency(char *title="",Bool_t isPrimary=kTRUE){
  31
  32   gSystem->Load("libITSUpgradeBase");
  33   gSystem->Load("libITSUpgradeRec");
  34   gSystem->Load("libITSUpgradeSim");
  35
  36   // setting up the configuration to be considered
  37   AliITSsegmentationUpgrade *seg = new AliITSsegmentationUpgrade();
  38   if(!seg){
  39     printf("no segmentation info available... Exiting");
  40     return;
  41   }
  42   radii.Set(seg->GetNLayers());
  43   for(Int_t l=0; l<seg->GetNLayers(); l++){
  44     radii.AddAt(seg->GetRadius(l),l);
  45   }
  46   delete seg;
  47
  48   TString titlet(title);
  49   titlet.ReplaceAll(" ","");
  50
  51   gROOT->SetStyle("Plain");
  52   gStyle->SetPalette(1);
  53   gStyle->SetOptStat(0);
  54   gStyle->SetOptTitle(0);
  55
  56
  57   //booking plots
  58   // pt
  59   Int_t nbins= 64;
  60   Double_t ptmax = 1.6;
  61
  62   TH1F *hPtRec = new TH1F("hPtRec"," p_{T} distribution of reconstructed tracks",nbins,0,ptmax);
  63   hPtRec->SetXTitle("p_{T}");
  64   hPtRec->SetYTitle(Form("entries / %3.0f MeV/c",1000*(ptmax/nbins)));
  65
  66   TH1F *hPtRecMc = new TH1F("hPtRecMc"," MC p_{T} distribution of reconstructed tracks",nbins,0,ptmax);
  67   hPtRecMc->SetXTitle("p_{T} GeV/c");
  68   hPtRecMc->SetYTitle(Form("entries / %3.0f MeV/c",1000*(ptmax/nbins)));
  69
  70   TH1F *hPtRecGood = new TH1F("hPtRecGood"," P_{t} distribution of reconstructed tracks [positive labels]",nbins,0,ptmax);
  71   hPtRecGood->SetXTitle("p_{T} GeV/c");
  72   hPtRecGood->SetYTitle(Form("entries / %3.0f MeV/c",1000*(ptmax/nbins)));
  73
  74   TH1F *hPtMc  = new TH1F("hPtMc", " p_{T} distribution of MC tracks",nbins,0,ptmax);
  75   hPtMc->SetXTitle("p_{T} GeV/c");
  76   hPtMc->SetYTitle(Form("entries / %3.0f MeV/c",1000*(ptmax/nbins)));
  77
  78   // Eta
  79   Int_t nbinsEta = 24;
  80   Double_t eta[2]={-1.2,1.2};
  81
  82   TH1F *hEtaRec = new TH1F("hEtaRec"," #eta distribution of reconstructed tracks",nbinsEta,eta[0],eta[1]);
  83   hEtaRec->SetXTitle("#eta");
  84   hEtaRec->SetYTitle(Form("entries / %3.3f ",((eta[1]-eta[0])/nbinsEta)));
  85
  86   TH1F *hEtaRecMc = new TH1F("hEtaRecMc"," #eta Mc distribution of reconstructed tracks",nbinsEta,eta[0],eta[1]);
  87   hEtaRecMc->SetXTitle("#eta");
  88   hEtaRecMc->SetYTitle(Form("entries / %3.3f ",((eta[1]-eta[0])/nbinsEta)));
  89
  90   TH1F *hEtaRecGood = new TH1F("hEtaRecGood"," #eta distribution of Good reconstructed tracks",nbinsEta,eta[0],eta[1]);
  91   hEtaRecGood->SetXTitle("#eta");
  92   hEtaRecGood->SetYTitle(Form("entries / %3.3f ",((eta[1]-eta[0])/nbinsEta)));
  93
  94   TH1F *hEtaMc  = new TH1F("hEtaMc", " #eta distribution of MC tracks",nbinsEta,eta[0],eta[1]);
  95   hEtaMc->SetXTitle("#eta");
  96   hEtaMc->SetYTitle(Form("entries / %3.3f ",((eta[1]-eta[0])/nbinsEta)));
  97
  98
  99   Double_t ptlimit=0.0;
 100
 101   Double_t cont=0;
 102   Double_t ntrack=0;
 103
 104   // Init simulation
 105   gAlice=NULL;
 106   AliRunLoader* runLoader = AliRunLoader::Open("galice.root");
 107   if (!runLoader) {
 108     Error("OpengAlice", "no files, please check the path");
 109     return ;
 110   }
 111   runLoader->LoadgAlice();
 112   gAlice = runLoader->GetAliRun();
 113   runLoader->LoadHeader();
 114   runLoader->LoadKinematics();
 115   runLoader->LoadTrackRefs();
 116
 117   //Trackref
 118   TTree *trackRefTree = 0x0;
 119   TClonesArray *trackRef = new TClonesArray("AliTrackReference",1000);
 120
 121   // ESD
 122   TFile* esdFile = TFile::Open("AliESDs.root");
 123   if (!esdFile || !esdFile->IsOpen()) {
 124     Error("CheckESD", "opening ESD file %s failed", "AliESDs.root");
 125     return ;
 126   }
 127   AliESDEvent * esd = new AliESDEvent;
 128   TTree* tree = (TTree*) esdFile->Get("esdTree");
 129   if (!tree) {
 130     Error("CheckESD", "no ESD tree found");
 131     return ;
 132   }
 133   esd->ReadFromTree(tree);
 134
 135   printf("setting the low pt value to %f \n",ptlimit);
 136   // Loop over events
 137
 138   for(Int_t i=0; i<runLoader->GetNumberOfEvents(); i++){
 139     //for(Int_t i=0; i< 1; i++){
 140     cout << "Event "<< i << endl;
 141     runLoader->GetEvent(i);
 142
 143     AliStack *stack = runLoader->Stack();
 144     trackRefTree=runLoader->TreeTR();
 145     TBranch *br = trackRefTree->GetBranch("TrackReferences");
 146     if(!br) {
 147       printf("no TR branch available , exiting \n");
 148       return;
 149     }
 150     br->SetAddress(&trackRef);
 151
 152     Int_t nParticles = stack->GetNtrack();
 153     trackRefTree=runLoader->TreeTR();
 154     trackRefTree->SetBranchAddress("TrackReferences",&trackRef);
 155     for(Int_t p=0; p<nParticles; p++){
 156       if(stack->Particle(p)->Pt()<ptlimit) continue;
 157       if(isPrimary && !stack->IsPhysicalPrimary(p)) continue;
 158       trackRefTree->GetEntry(stack->TreeKEntry(p));
 159
 160       if(IsTrackable(trackRef)) {
 161         cont++;
 162         hPtMc->Fill((stack->Particle(p))->Pt());
 163         hEtaMc->Fill((stack->Particle(p))->Eta());
 164       }
 165     }//loop on kine particle
 166     cout << " # trackable "<< cont;
 167
 168     // ESD
 169     tree->GetEvent(i);
 170     for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
 171       AliESDtrack* track = esd->GetTrack(iTrack);
 172       if(track->Pt()<ptlimit) continue;
 173       if(isPrimary && !stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) continue;
 174       if(track->GetNcls(0)<minNpoints) continue; // selection compatible with IsSelected
 175       ntrack++;
 176       hPtRec->Fill(track->Pt());
 177       hEtaRec->Fill(track->Eta());
 178       if(track->GetLabel()>=0) {
 179         hPtRecGood->Fill(stack->Particle(track->GetLabel())->Pt());
 180         hEtaRecGood->Fill(stack->Particle(track->GetLabel())->Eta());
 181       }
 182       hPtRecMc->Fill(stack->Particle(TMath::Abs(track->GetLabel()))->Pt());
 183       hEtaRecMc->Fill(stack->Particle(TMath::Abs(track->GetLabel()))->Eta());
 184     }// loop on tracks
 185     cout << " tracks " << ntrack << endl;
 186
 187   }//loop events
 188
 189   cout << "Global Efficiency  "<< ntrack/cont << endl;
 190
 191   TCanvas *c = new TCanvas();
 192   c->cd();
 193   c->cd()->SetLogy();
 194   c->cd()->SetGridy();
 195   c->cd()->SetGridx();
 196   c->cd()->SetObjectStat(0);
 197   hPtMc->SetMarkerStyle(20);
 198   hPtMc->Draw("err");
 199   hPtRec->SetLineColor(kRed);
 200   hPtRec->SetMarkerStyle(20);
 201   hPtRec->SetMarkerColor(kRed);
 202   hPtRec->Sumw2();
 203   hPtRec->Draw("sameerr");
 204   hPtRecMc->SetLineColor(kBlue);
 205   hPtRecMc->SetMarkerStyle(20);
 206   hPtRecMc->SetMarkerColor(kBlue);
 207   hPtRecMc->Sumw2();
 208   hPtRecMc->Draw("sameerr");
 209   hPtRecGood->SetLineColor(kGreen);
 210   hPtRecGood->SetMarkerStyle(20);
 211   hPtRecGood->SetMarkerColor(kGreen);
 212   hPtRecGood->Sumw2();
 213   hPtRecGood->Draw("sameerr");
 214
 215   TLegend *leg = new TLegend(0.5,0.7,0.99,0.95);
 216   leg->SetHeader(Form("%s",title));
 217   leg->AddEntry(hPtRecGood,"Good tracks (positive labels)","p");
 218   leg->AddEntry(hPtRec,"Reconstructed tracks (Rec p_{T})","p");
 219   leg->AddEntry(hPtRecMc,"Reconstructed tracks ","p");
 220   leg->AddEntry(hPtMc,"MC trackable particles","p");
 221   leg->Draw();
 222   c->SaveAs(Form("ptdist_%s.png",titlet.Data()));
 223
 224   TCanvas *cc = new TCanvas();
 225   cc->cd();
 226   cc->cd()->SetGridy();
 227   cc->cd()->SetGridx();
 228
 229   Double_t max = 1.2;
 230
 231   TH1F *effrec=new TH1F("hEffRec","tracking efficiency of reconstructed tracks",nbins,0,ptmax);
 232   effrec->Divide(hPtRecMc,hPtMc,1,1,"b");
 233   effrec->SetXTitle("p_{T} GeV/c");
 234   effrec->SetLineColor(kBlue);
 235   effrec->SetMarkerStyle(20);
 236   effrec->SetMarkerColor(kBlue);
 237   effrec->SetMaximum(max);
 238   effrec->SetMinimum(0);
 239   effrec->DrawCopy("err");
 240
 241
 242   TH1F *effrecgood=new TH1F("hEffRecGood","tracking efficiency of reconstructed tracks [positive labels]",nbins,0,ptmax);
 243   effrecgood->SetXTitle("P_{T}");
 244   effrecgood->Divide(hPtRecGood,hPtMc,1,1,"b");
 245   effrecgood->SetLineColor(kGreen);
 246   effrecgood->SetMarkerStyle(20);
 247   effrecgood->SetMarkerColor(kGreen);
 248   effrecgood->SetMaximum(max);
 249   effrecgood->SetMinimum(0);
 250   effrecgood->DrawCopy("sameerr");
 251
 252   TH1F *purity=new TH1F("purity","track purity [positive labels / all tracks]",nbins,0,ptmax);
 253   purity->Divide(hPtRecGood,hPtRecMc,1,1,"b");
 254   purity->SetLineColor(kGray);
 255   purity->SetMarkerStyle(28);
 256   purity->SetMarkerColor(kOrange);
 257   purity->SetMaximum(max);
 258   purity->SetMinimum(0);
 259   purity->DrawCopy("sameerr");
 260
 261   TH1F *fakes=new TH1F("fakes","fake ratio [negative labels / reconstructable]",nbins,0,ptmax);
 262
 263   fakes->Add(hPtRecMc,hPtRecGood,1,-1);
 264   fakes->Divide(hPtMc);
 265   fakes->SetLineColor(kViolet);
 266   fakes->SetLineStyle(3);
 267   fakes->SetLineWidth(3);
 268   fakes->SetMarkerColor(kViolet);
 269   fakes->SetMaximum(max);
 270   fakes->SetMinimum(0);
 271   fakes->DrawCopy("samehist");
 272
 273   TLegend *legEff = new TLegend(0.1,0.83,0.65,0.99);
 274   legEff->SetHeader(Form("%s",title));
 275   legEff->AddEntry(effrec,"Overall Efficiency ","p");
 276   legEff->AddEntry(effrecgood,"Efficiency good tracks (positive labels)","p");
 277   legEff->AddEntry(purity,"Purity (positive labels / All tracks)","p");
 278   legEff->AddEntry(fakes,"Fakes","l");
 279   legEff->Draw();
 280
 281   cc->SaveAs(Form("efficiency_%s.png",titlet.Data()));
 282
 283
 284   TCanvas *cEta = new TCanvas();
 285   cEta->cd();
 286   cEta->cd()->SetGridy();
 287   cEta->cd()->SetGridx();
 288   cEta->cd()->SetObjectStat(0);
 289   hEtaMc->SetMarkerStyle(20);
 290   hEtaMc->Draw("err");
 291   hEtaMc->SetMinimum(0);
 292   hEtaRec->SetLineColor(kRed);
 293   hEtaRec->SetMarkerStyle(20);
 294   hEtaRec->SetMarkerColor(kRed);
 295   hEtaRec->Sumw2();
 296   hEtaRec->Draw("sameerr");
 297   hEtaRecMc->SetLineColor(kBlue);
 298   hEtaRecMc->SetMarkerStyle(20);
 299   hEtaRecMc->SetMarkerColor(kBlue);
 300   hEtaRecMc->Sumw2();
 301   hEtaRecMc->Draw("sameerr");
 302   hEtaRecGood->SetLineColor(kGreen);
 303   hEtaRecGood->SetMarkerStyle(20);
 304   hEtaRecGood->SetMarkerColor(kGreen);
 305   hEtaRecGood->Sumw2();
 306   hEtaRecGood->Draw("sameerr");
 307
 308   TLegend *legEta = new TLegend(0.2,0.3,0.59,0.55);
 309   legEta->SetHeader(Form("%s",title));
 310   legEta->AddEntry(hEtaRecGood,"Good tracks (positive labels)","p");
 311   legEta->AddEntry(hEtaRec,"Reconstructed tracks (Rec p_{T})","p");
 312   legEta->AddEntry(hEtaRecMc,"Reconstructed tracks ","p");
 313   legEta->AddEntry(hEtaMc,"MC trackable particles","p");
 314   legEta->Draw();
 315   cEta->SaveAs(Form("etadist_%s.png",title));
 316
 317   TCanvas *ccEta = new TCanvas();
 318   ccEta->cd();
 319   ccEta->cd()->SetGridy();
 320   ccEta->cd()->SetGridx();
 321
 322   TH1F *effrecEta=new TH1F("hEffRecEta","#eta tracking efficiency of reconstructed tracks",nbinsEta,eta[0],eta[1]);
 323   effrecEta->Divide(hEtaRecMc,hEtaMc,1,1,"b");
 324   effrecEta->SetXTitle("#eta");
 325   effrecEta->SetLineColor(kBlue);
 326   effrecEta->SetMarkerStyle(20);
 327   effrecEta->SetMarkerColor(kBlue);
 328   effrecEta->SetMaximum(max);
 329   effrecEta->SetMinimum(0);
 330   effrecEta->DrawCopy("err");
 331
 332
 333   TH1F *effrecgoodEta=new TH1F("hEffRecGoodEta","#eta tracking efficiency of reconstructed tracks [positive labels]",nbinsEta,eta[0],eta[1]);
 334   effrecgoodEta->SetXTitle("eta");
 335   effrecgoodEta->Divide(hEtaRecGood,hEtaMc,1,1,"b");
 336   effrecgoodEta->SetLineColor(kGreen);
 337   effrecgoodEta->SetMarkerStyle(20);
 338   effrecgoodEta->SetMarkerColor(kGreen);
 339   effrecgoodEta->SetMaximum(max);
 340   effrecgoodEta->SetMinimum(0);
 341   effrecgoodEta->DrawCopy("sameerr");
 342
 343   TH1F *purityEta=new TH1F("purityEta","#eta track purity [positive labels / all tracks]",nbinsEta,eta[0],eta[1]);
 344   purityEta->Divide(hEtaRecGood,hEtaRecMc,1,1,"b");
 345   purityEta->SetLineColor(kGray);
 346   purityEta->SetMarkerStyle(28);
 347   purityEta->SetMarkerColor(kOrange);
 348   purityEta->SetMaximum(max);
 349   purityEta->SetMinimum(0);
 350   purityEta->DrawCopy("sameerr");
 351
 352   TH1F *fakesEta=new TH1F("fakesEta","#eta fake ratio [negative labels / reconstructable]",nbinsEta,eta[0],eta[1]);
 353   fakesEta->Add(hEtaRecMc,hEtaRecGood,1,-1);
 354   fakesEta->Divide(hEtaMc);
 355   fakesEta->SetLineColor(kViolet);
 356   fakesEta->SetLineStyle(3);
 357   fakesEta->SetLineWidth(3);
 358   fakesEta->SetMarkerColor(kViolet);
 359   fakesEta->SetMaximum(max);
 360   fakesEta->SetMinimum(0);
 361   fakesEta->DrawCopy("samehist");
 362
 363   TLegend *legEffEta = new TLegend(0.3,0.3,0.75,0.55);
 364   legEffEta->SetHeader(Form("%s",title));
 365   legEffEta->AddEntry(effrecEta,"Overall Efficiency ","p");
 366   legEffEta->AddEntry(effrecgoodEta,"Efficiency good tracks (positive labels)","p");
 367   legEffEta->AddEntry(purityEta,"Purity (positive labels / All tracks)","p");
 368   legEffEta->AddEntry(fakesEta,"Fakes","l");
 369   legEffEta->Draw();
 370
 371   ccEta->SaveAs(Form("efficiencyEta_%s.png",titlet.Data()));
 372
 373
 374 }// main
 375
 376 //___________________________________________________
 377 Bool_t IsTrackable(TClonesArray *trackRef){
 378
 379   Bool_t isOk=kFALSE;
 380
 381   Int_t nTrackRef =0;
 382   TArrayF isInLayer(radii.GetSize());
 383   for(Int_t l=0; l<radii.GetSize(); l++ ) isInLayer.AddAt(0,l);
 384   for(Int_t nT =0; nT<trackRef->GetEntries(); nT++){
 385     AliTrackReference *trR = (AliTrackReference*)trackRef->At(nT);
 386     if(!trR) continue;
 387     if(trR->DetectorId()!=0)continue;
 388     Double_t rPart = TMath::Sqrt(trR->X()*trR->X()+trR->Y()*trR->Y());
 389     for(Int_t iLay=0;iLay<radii.GetSize();iLay++){
 390       if(TMath::Abs(rPart-radii.At(iLay))<0.01) isInLayer.AddAt(1,iLay);
 391     }
 392   }
 393
 394   for(Int_t iLayer=0;iLayer<radii.GetSize();iLayer++){
 395     if(isInLayer.At(iLayer)) nTrackRef++;
 396   }
 397
 398   if(nTrackRef>=minNpoints) isOk=kTRUE;
 399   return isOk;
 400 }
 401