Updated version of the proof-enabled comparison macros: now available for ITS, TPC...
[u/mrichter/AliRoot.git] / PWG1 / comparison / AliTPCComparisonTask.cxx
index f303ba6be95050a793786d43f61443fd7d969eb0..b0cb745bc1092c039a4d68e4c1960dd971b7f637 100644 (file)
@@ -18,8 +18,7 @@
 #include "AliVEvent.h"
 #include "AliESDEvent.h"
 #include "AliMCEvent.h"
-#include "AliAODEvent.h"
-
+#include "AliESDtrack.h"
 #include "AliTrackReference.h"
 #include "AliMCComparisonTrack.h"
 
 
 ClassImp(AliTPCComparisonTask)
 
+
+AliTPCComparisonTask::AliTPCComparisonTask()
+  : AliAnalysisTaskSE("AliTPCComparisonTask"),
+    fListOfHistos(0),
+    hgood(0),
+    hfound(0),
+    hfake(0),
+    hp(0),
+    hl(0),
+    hpt(0),
+    hmpt(0),
+    he(0),
+    hep(0),
+    hgoodPhi(0),
+    hfoundPhi(0)
+{
+  // Default constructor
+  AliInfo("Default constructor AliTPCComparisonTask");
+  // Define input and output slots here
+  // Input slot #0 works with a TChain
+  DefineInput(0, TChain::Class());
+  // Output slot #1 TList
+  DefineOutput(1, TList::Class());
+}
+
+
 AliTPCComparisonTask::AliTPCComparisonTask(const char* name)
-  : AliAnalysisTaskSE(name), fListOfHistos(0), hgood(0), hfound(0),
-  hfake(0), hp(0), hl(0), hpt(0), hmpt(0), he(0), hep(0), hgoodPhi(0),
-  hfoundPhi(0)
+  : AliAnalysisTaskSE(name),
+    fListOfHistos(0),
+    hgood(0),
+    hfound(0),
+    hfake(0),
+    hp(0),
+    hl(0),
+    hpt(0),
+    hmpt(0),
+    he(0),
+    hep(0),
+    hgoodPhi(0),
+    hfoundPhi(0)
 {
-    // Constructor
+  // Constructor
   AliInfo("Constructor AliTPCComparisonTask");
-    // Define input and output slots here
-    // Input slot #0 works with a TChain
+  // Define input and output slots here
+  // Input slot #0 works with a TChain
   DefineInput(0, TChain::Class());
-    // Output slot #1 TList
+  // Output slot #1 TList
   DefineOutput(1, TList::Class());
 }
 
 
-
 void AliTPCComparisonTask::UserCreateOutputObjects() 
 {
-    // Create histograms
-    // Called once
+  // Create histograms
+  // Called once
   AliInfo("AliTPCComparisonTask::UserCreateOutputObjects");
-    // Create output container
+  // Create output container
   fListOfHistos = new TList();
   
-  hgood = new TH1F("hgood", "Good tracks", 30, 0.2, 6.1);
-  hfound= new TH1F("hfound", "Found tracks", 30, 0.2, 6.1);
-  hfake = new TH1F("hfake",  "Fake tracks", 30, 0.2, 6.1);
+  hgood = new TH1F("hgood", "Pt for good tracks", 34, 0.2, 7.0);
+  hfound = new TH1F("hfound", "Pt for found tracks", 34, 0.2, 7.0);
+  hfake = new TH1F("hfake",  "Pt for fake tracks", 34, 0.2, 7.0);
   hp = new TH1F("hp", "PHI resolution", 50, -20., 20.);
   hl = new TH1F("hl", "LAMBDA resolution", 50, -20., 20.);
   hpt = new TH1F("hpt", "Relative Pt resolution", 30, -10., 10.);
   hmpt = new TH1F("hmpt", "Relative Pt resolution (pt>4GeV/c)", 30, -60., 60.);
   he = new TH1F("he", "dE/dX for pions with 0.4<p<0.5 GeV/c", 50, 0., 100.);
   hep = new TH2F("hep", "dE/dX vs momentum", 50, 0., 2., 50, 0., 400.);
-  hgoodPhi = new TH1F("hgoodPhi", "Phi for Good tracks", 90, -TMath::Pi(), TMath::Pi());
-  hfoundPhi = new TH1F("hfoundPhi", "Phi for Found tracks", 90, -TMath::Pi(), TMath::Pi());
+  hgoodPhi = new TH1F("hgoodPhi", "Phi for good tracks", 90, 0., 2.*TMath::Pi());
+  hfoundPhi = new TH1F("hfoundPhi", "Phi for found tracks", 90, 0., 2.*TMath::Pi());
   
   fListOfHistos->Add(hgood);
   fListOfHistos->Add(hfound);
@@ -79,80 +113,76 @@ void AliTPCComparisonTask::UserCreateOutputObjects()
 
 void AliTPCComparisonTask::UserExec(Option_t *) 
 {
-    // Main loop
-    // Called for each event
-
-    // MC information
+  // Main loop
+  // Called for each event
+  
+  // MC information
   AliMCEvent* mcEvent = MCEvent();
   if (!mcEvent) {
-     Printf("ERROR: Could not retrieve MC event");
-     return;
+    Printf("ERROR: Could not retrieve MC event");
+    return;
   }
-
-  Printf("MC particles: %d", mcEvent->GetNumberOfTracks());
   
   Int_t nt = 0;
-    
+  
   TClonesArray dummy("AliMCComparisonTrack",1000), *refs=&dummy;
-    
-  for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) 
-  {
+  
+  // Loop over all MC tracks and select "good" tracks
+  for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {
     AliMCParticle* track = mcEvent->GetTrack(iTracks);
-    if (!track) 
-    {
+    if (!track) {
       Printf("ERROR: Could not receive track %d (mc loop)", iTracks);
       continue;
     }
     
-      // Track selection
-    if (track->Pt()<0.100) continue;    
-    if (TMath::Abs(track->Pz()/track->Pt())>0.999) continue;
+    // Track selection
+    if (track->Pt() < 0.2) continue; 
+    if (track->Pt() > 7.0) continue;   
+    if (TMath::Abs(track->Pz()/track->Pt()) > 0.999) continue;
     
-    Double_t vx=track->Xv(),vy=track->Yv(),vz=track->Zv();
-    if (TMath::Sqrt(vx*vx+vy*vy)>3.5) continue;
+    Double_t vx = track->Xv(),vy = track->Yv(),vz = track->Zv();
+    if (TMath::Sqrt(vx*vx+vy*vy) > 3.5) continue;
     if (TMath::Abs(vz) > 50.) continue; 
     
+    // Loop over Track References
     Bool_t LabelTPC = kFALSE;
     AliTrackReference* trackRef = 0;
-    for (Int_t iTrackRef = 0; iTrackRef  < track->GetNumberOfTrackReferences(); iTrackRef++)
-    {
-       trackRef = track->GetTrackReference(iTrackRef);
-       if(trackRef)
-       {
-          Int_t detectorId = trackRef->DetectorId(); 
-         if (detectorId == AliTrackReference::kTPC)
-         {         
-           LabelTPC = kTRUE;
-           break;
-         }
-       }      
+    for (Int_t iTrackRef = 0; iTrackRef  < track->GetNumberOfTrackReferences(); iTrackRef++) {
+      trackRef = track->GetTrackReference(iTrackRef);
+      if(trackRef) {
+       Int_t detectorId = trackRef->DetectorId(); 
+       if (detectorId == AliTrackReference::kTPC) {        
+         LabelTPC = kTRUE;
+         break;
+       }
+      }      
     }
-    if (LabelTPC)
-    {
-      AliMCComparisonTrack *ref=new((*refs)[nt]) AliMCComparisonTrack();
+    // "Good" tracks
+    if (LabelTPC) {
+      AliMCComparisonTrack* ref = new((*refs)[nt]) AliMCComparisonTrack();
       ref->SetLabel(iTracks);
-      TParticle* particle= track->Particle();
+      TParticle* particle = track->Particle();
       Int_t pdg = particle->GetPdgCode();
       ref->SetPDG(pdg); 
-      ref->SetTPCMomentum(trackRef->Px(),trackRef->Py(),trackRef->Pz());
-      Float_t Pt = TMath::Sqrt(trackRef->Px()*trackRef->Px() + 
-        trackRef->Py()*trackRef->Py());
+      ref->SetPz(trackRef->Pz());
+      Float_t Pt = trackRef->Pt();
+      ref->SetPt(Pt);
       hgood->Fill(Pt);
-      Float_t phig=TMath::ATan2(trackRef->Py(),trackRef->Px());
+      Float_t phig = trackRef->Phi();
+      ref->SetPhi(phig);
       hgoodPhi->Fill(phig);
       nt++;  
     }    
   } //track loop 
-
   
-    // ESD information  
-  AliVEvent *event = InputEvent();
+  
+  // ESD information  
+  AliVEvent* event = InputEvent();
   if (!event) {
-     Printf("ERROR: Could not retrieve event");
-     return;
+    Printf("ERROR: Could not retrieve event");
+    return;
   }
-  
-  
+    
   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
   
   Bool_t iFound;
@@ -162,36 +192,31 @@ void AliTPCComparisonTask::UserExec(Option_t *)
   
   Int_t MCgoods = refs->GetEntriesFast();
   
-  for (Int_t k = 0; k < MCgoods; k++) 
-  {
-    AliMCComparisonTrack *ref=(AliMCComparisonTrack*)refs->UncheckedAt(k); 
+  for (Int_t k = 0; k < MCgoods; k++) {
+    AliMCComparisonTrack* ref = (AliMCComparisonTrack*)refs->UncheckedAt(k);
+    if (!ref) continue;
     Int_t MCLabel = ref->GetLabel();
-    Float_t ptg = 
-      TMath::Sqrt(ref->GetTPCPx()*ref->GetTPCPx() + 
-      ref->GetTPCPy()*ref->GetTPCPy());
-    Float_t PhiG=TMath::ATan2(ref->GetTPCPy(),ref->GetTPCPx());
+    Float_t ptg = ref->GetPt();
+    Float_t Phig = ref->GetPhi();
 
     iFound = kFALSE;
-    for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) 
-    {
+    for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
       AliESDtrack* track = esd->GetTrack(iTrack);  
-      if (!track) 
-      {
+      if (!track) {
         Printf("ERROR: Could not receive track %d", iTrack);
         continue;
       }
-      Int_t TPCLabel =  track->GetTPCLabel();
       
-      if (MCLabel == TMath::Abs(TPCLabel))
-      {          
-        if (MCLabel == TPCLabel)
-       {
+      if (! track->IsOn(AliESDtrack::kTPCrefit)) continue;     
+      
+      Int_t TPCLabel =  track->GetTPCLabel();      
+      if (MCLabel == TMath::Abs(TPCLabel)) {     
+        if (MCLabel == TPCLabel) {
          nfound++;
          hfound->Fill(ptg);
-         hfoundPhi->Fill(PhiG);
+         hfoundPhi->Fill(Phig);
        } 
-       else
-       {
+       else {
          nfake++;
          hfake->Fill(ptg);
        }
@@ -199,69 +224,62 @@ void AliTPCComparisonTask::UserExec(Option_t *)
        
       
         Double_t pxpypz[3]; track->GetInnerPxPyPz(pxpypz);
-        Float_t phi=TMath::ATan2(pxpypz[1],pxpypz[0]);
-        if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
-        if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
-        Double_t pt=TMath::Sqrt(pxpypz[0]*pxpypz[0]+pxpypz[1]*pxpypz[1]);
-        Float_t lam=TMath::ATan2(pxpypz[2],pt); 
-        Float_t pt_1=1/pt;
-
-        Int_t pdg=(Int_t)ref->GetPDG();  
-        if (TMath::Abs(pdg)==11 && ptg>4.) 
-       {
-           //high momentum electrons
+        Float_t phi = TMath::ATan2(pxpypz[1],pxpypz[0]);
+       if (phi < 0) phi += 2*TMath::Pi();
+        Double_t pt = TMath::Sqrt(pxpypz[0]*pxpypz[0]+pxpypz[1]*pxpypz[1]);
+        Float_t lam = TMath::ATan2(pxpypz[2],pt); 
+        Float_t pt_1 = 1/pt;
+
+        Int_t pdg = ref->GetPDG();  
+        if (TMath::Abs(pdg) == 11 && ptg>4.) {
+         //high momentum electrons
           hmpt->Fill((pt_1 - 1/ptg)/(1/ptg)*100.);
         } 
-       else 
-       {
-          Float_t phig=TMath::ATan2(ref->GetTPCPy(),ref->GetTPCPx());
-          hp->Fill((phi - phig)*1000.);          
-         Float_t lamg=TMath::ATan2(ref->GetTPCPz(),ptg);
+       else {
+          hp->Fill((phi - Phig)*1000.);          
+         Float_t lamg = TMath::ATan2(ref->GetPz(),ptg);
           hl->Fill((lam - lamg)*1000.);
           hpt->Fill((pt_1 - 1/ptg)/(1/ptg)*100.);
-        }      
+        }
        
-       Float_t mom=pt/TMath::Cos(lam);
-        Float_t dedx=track->GetTPCsignal();
-        hep->Fill(mom,dedx,1.);
-        if (TMath::Abs(pdg)==211) //pions
-          if (mom>0.4 && mom<0.5) 
-          {
-             he->Fill(dedx,1.);
-           }
-          
+       Float_t mom = pt/TMath::Cos(lam);
+        Float_t dedx = track->GetTPCsignal();
+        hep->Fill(mom, dedx, 1.);
+        if (TMath::Abs(pdg) == 211) { // pions
+         if (mom > 0.4 && mom < 0.5) {
+            he->Fill(dedx,1.);
+          }
+       }   
           
         break; 
       }
     }  
-    if (!iFound)
-    {
+    if (!iFound) {
       nlost++;
     } 
   }
      
   Printf(" Results: " );
   Printf(" Found %d Fake %d Lost %d ", nfound, nfake, nlost);
-    
+  
   refs->Clear();
-
+  
   // Post output data.
   PostData(1, fListOfHistos);
 }      
 
-//________________________________________________________________________
+
 void AliTPCComparisonTask::Terminate(Option_t *) 
 {
   // Draw result to the screen
   // Called once at the end of the query 
-
+  
   fListOfHistos = dynamic_cast<TList*>(GetOutputData(1));
   if (!fListOfHistos) {
     Printf("ERROR: fListOfHistos not available");
     return;
   }
 
-
   hgood = dynamic_cast<TH1F*>(fListOfHistos->At(0));
   hfound = dynamic_cast<TH1F*>(fListOfHistos->At(1));
   hfake = dynamic_cast<TH1F*>(fListOfHistos->At(2));  
@@ -274,92 +292,93 @@ void AliTPCComparisonTask::Terminate(Option_t *)
   hgoodPhi = dynamic_cast<TH1F*>(fListOfHistos->At(9));
   hfoundPhi = dynamic_cast<TH1F*>(fListOfHistos->At(10));
 
-
   gStyle->SetOptStat(111110);
   gStyle->SetOptFit(1);
    
-  TCanvas* c1=new TCanvas("c1","",0,0,700,850);
+  TCanvas* c1 = new TCanvas("c1", "", 0, 0, 700, 850);
 
-  Int_t minc=33;
+  Int_t minc = 33;
   
-   TPad *p1=new TPad("p1","",0,0.3,.5,.6); p1->Draw();
-   p1->cd(); p1->SetFillColor(42); p1->SetFrameFillColor(10);
-   hp->SetFillColor(4);  hp->SetXTitle("(mrad)");
-   if (hp->GetEntries()<minc) hp->Draw(); else hp->Fit("gaus"); c1->cd();
-
-   TPad *p2=new TPad("p2","",0.5,.3,1,.6); p2->Draw();
-   p2->cd(); p2->SetFillColor(42); p2->SetFrameFillColor(10);
-   hl->SetXTitle("(mrad)");
-   if (hl->GetEntries()<minc) hl->Draw(); else hl->Fit("gaus"); c1->cd();
-
-   TPad *p3=new TPad("p3","",0,0,0.5,0.3); p3->Draw();
-   p3->cd(); p3->SetFillColor(42); p3->SetFrameFillColor(10);
-   hpt->SetXTitle("(%)");
-   if (hpt->GetEntries()<minc) hpt->Draw(); else hpt->Fit("gaus"); c1->cd();
-
-   TPad *p4=new TPad("p4","",0.5,0,1,0.3); p4->Draw();
-   p4->cd(); p4->SetFillColor(42); p4->SetFrameFillColor(10);
-   hmpt->SetXTitle("(%)");
-   if (hmpt->GetEntries()<minc) hmpt->Draw(); else hmpt->Fit("gaus"); c1->cd(); 
-
+  TPad* p1 = new TPad("p1", "", 0., 0.3, 0.5, 0.6); p1->Draw();
+  p1->cd(); p1->SetFillColor(42); p1->SetFrameFillColor(10);
+  hp->SetFillColor(4); hp->SetXTitle("(mrad)");
+  if (hp->GetEntries() < minc) hp->Draw(); else hp->Fit("gaus"); c1->cd();
+
+  TPad* p2 = new TPad("p2", "", 0.5, 0.3, 1, 0.6); p2->Draw();
+  p2->cd(); p2->SetFillColor(42); p2->SetFrameFillColor(10);
+  hl->SetFillColor(4); hl->SetXTitle("(mrad)");
+  if (hl->GetEntries() < minc) hl->Draw(); else hl->Fit("gaus"); c1->cd();
   
-   TPad *p5=new TPad("p5","",0,0.6,1,1); p5->Draw(); p5->cd();
-   p5->SetFillColor(41); p5->SetFrameFillColor(10);
-   hfound->Sumw2(); hgood->Sumw2(); hfake->Sumw2();
-   TH1F* hg = new TH1F("hg", "Efficiency for good tracks", 30, 0.2, 6.1);
-   TH1F* hf = new TH1F("hf", "Efficiency for fake tracks", 30, 0.2, 6.1);
-   hg->Divide(hfound,hgood,1.,1.,"B");
-   hf->Divide(hfake,hgood,1.,1.,"B");
-   hg->SetMaximum(1.4);
-   hg->SetYTitle("Tracking efficiency");
-   hg->SetXTitle("Pt (GeV/c)");
-   hg->Draw();
-
-   TLine *line1 = new TLine(0.1,1.0,6.1,1.0); line1->SetLineStyle(4);
-   line1->Draw("same");
-   TLine *line2 = new TLine(0.1,0.9,6.1,0.9); line2->SetLineStyle(4);
-   line2->Draw("same");
-
-   hf->SetFillColor(1);
-   hf->SetFillStyle(3013);
-   hf->SetLineColor(2);
-   hf->SetLineWidth(2);
-   hf->Draw("histsame");
-   TText *text = new TText(0.461176,0.248448,"Fake tracks");
-   text->SetTextSize(0.05);
-   text->Draw();
-   text = new TText(0.453919,1.11408,"Good tracks");
-   text->SetTextSize(0.05);
-   text->Draw();
-     
-   TCanvas *c2=new TCanvas("c2","",320,32,530,590);
-
-   TPad *p6=new TPad("p6","",0.,0.,1.,.5); p6->Draw();
-   p6->cd(); p6->SetFillColor(42); p6->SetFrameFillColor(10);
-   he->SetFillColor(2); he->SetFillStyle(3005);
-   he->SetXTitle("Arbitrary Units");
-   if (he->GetEntries()<minc) he->Draw(); else he->Fit("gaus"); c2->cd();
-
-   TPad *p7=new TPad("p7","",0.,0.5,1.,1.); p7->Draw();
-   p7->cd(); p7->SetFillColor(42); p7->SetFrameFillColor(10);
-   hep->SetFillColor(2); hep->SetFillStyle(3005);
-   hep->SetXTitle("p (GeV/c)"); hep->SetYTitle("dE/dX (Arb. Units)");
-   hep->Draw(); c1->cd();
-   
-   TCanvas* c3=new TCanvas("c3","",10,10,510,510);
-   c3->SetFillColor(42); c3->SetFrameFillColor(10);
-   hfoundPhi->Sumw2(); hgoodPhi->Sumw2();
-   TH1F* hgphi = new TH1F("hgphi", "Efficiency for good tracks (Phi)", 
-                         90, -TMath::Pi(), TMath::Pi());
-   hgphi->Divide(hfoundPhi,hgoodPhi,1.,1.,"B"); 
-   hgphi->SetYTitle("Tracking efficiency");
-   hgphi->SetXTitle("Phi (rad)");
-   hgphi->Draw();
-
-   TFile fc("AliTPCComparison.root","RECREATE");
-   c1->Write();
-   c2->Write();
-   c3->Write();
-   fc.Close();
+  TPad* p3 = new TPad("p3", "", 0., 0., 0.5, 0.3); p3->Draw();
+  p3->cd(); p3->SetFillColor(42); p3->SetFrameFillColor(10);
+  hpt->SetFillColor(2); hpt->SetXTitle("(%)");
+  if (hpt->GetEntries() < minc) hpt->Draw(); else hpt->Fit("gaus"); c1->cd();
+  
+  TPad* p4 = new TPad("p4", "", 0.5, 0., 1., 0.3); p4->Draw();
+  p4->cd(); p4->SetFillColor(42); p4->SetFrameFillColor(10);
+  hmpt->SetFillColor(6); hmpt->SetXTitle("(%)");
+  if (hmpt->GetEntries() < minc) hmpt->Draw(); else hmpt->Fit("gaus"); c1->cd(); 
+  
+  TPad* p5 = new TPad("p5", "", 0., 0.6, 1., 1.); p5->Draw(); p5->cd();
+  p5->SetFillColor(41); p5->SetFrameFillColor(10);
+  hfound->Sumw2(); hgood->Sumw2(); hfake->Sumw2();
+  TH1F* hg = new TH1F("hg", "Efficiency for good tracks", 34, 0.2, 7.0);
+  TH1F* hf = new TH1F("hf", "Efficiency for fake tracks", 34, 0.2, 7.0);
+  hg->Divide(hfound, hgood, 1., 1., "B");
+  hf->Divide(hfake, hgood, 1., 1., "B");
+  hg->SetLineColor(4); hg->SetLineWidth(2);
+  hg->SetMaximum(1.4);
+  hg->SetYTitle("Tracking efficiency");
+  hg->SetXTitle("Pt (GeV/c)");
+  hg->Draw();
+  
+  TLine* line1 = new TLine(0.2, 1.0, 7.0, 1.0); line1->SetLineStyle(4);
+  line1->Draw("same");
+  TLine* line2 = new TLine(0.2, 0.9, 7.0, 0.9); line2->SetLineStyle(4);
+  line2->Draw("same");
+  
+  hf->SetFillColor(1);
+  hf->SetFillStyle(3013);
+  hf->SetLineColor(2);
+  hf->SetLineWidth(2);
+  hf->Draw("histsame");
+  TText* text = new TText(0.461176, 0.248448, "Fake tracks");
+  text->SetTextSize(0.05);
+  text->Draw();
+  text = new TText(0.453919, 1.11408, "Good tracks");
+  text->SetTextSize(0.05);
+  text->Draw();
+  
+  TCanvas* c2 = new TCanvas("c2", "", 320, 32, 530, 590);
+  
+  TPad* p6 = new TPad("p6", "", 0., 0., 1., 0.5); p6->Draw();
+  p6->cd(); p6->SetFillColor(42); p6->SetFrameFillColor(10);
+  he->SetFillColor(2); he->SetFillStyle(3005);
+  he->SetXTitle("Arbitrary Units");
+  if (he->GetEntries() < minc) he->Draw(); else he->Fit("gaus"); c2->cd();
+  
+  TPad* p7 = new TPad("p7", "", 0., 0.5, 1., 1.); p7->Draw();
+  p7->cd(); p7->SetFillColor(42); p7->SetFrameFillColor(10);
+  hep->SetFillColor(2); hep->SetFillStyle(3005);
+  hep->SetMarkerStyle(8);
+  hep->SetMarkerSize(0.4);
+  hep->SetXTitle("p (GeV/c)"); hep->SetYTitle("dE/dX (Arb. Units)");
+  hep->Draw(); c1->cd();
+  
+  TCanvas* c3 = new TCanvas("c3", "", 10, 10, 510, 510);
+  c3->SetFillColor(42); c3->SetFrameFillColor(10);
+  hfoundPhi->Sumw2(); hgoodPhi->Sumw2();
+  TH1F* hgphi = new TH1F("hgphi", "Efficiency for good tracks (Phi)", 
+                        90, 0., 2.*TMath::Pi());
+  hgphi->Divide(hfoundPhi, hgoodPhi, 1., 1., "B"); 
+  hgphi->SetYTitle("Tracking efficiency");
+  hgphi->SetXTitle("Phi (rad)");
+  hgphi->Draw();
+  
+  TFile fc("AliTPCComparison.root", "RECREATE");
+  c1->Write();
+  c2->Write();
+  c3->Write();
+  fc.Close();
 
 }