Updated version of the proof-enabled comparison macros: now available for ITS, TPC...
authorhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 22 Jul 2008 09:07:14 +0000 (09:07 +0000)
committerhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 22 Jul 2008 09:07:14 +0000 (09:07 +0000)
15 files changed:
PWG1/comparison/AF-v4-12.par [deleted file]
PWG1/comparison/AliITSComparisonTask.cxx [new file with mode: 0644]
PWG1/comparison/AliITSComparisonTask.h [new file with mode: 0644]
PWG1/comparison/AliMCComparisonTrack.cxx
PWG1/comparison/AliMCComparisonTrack.h
PWG1/comparison/AliTOFComparisonTask.cxx [new file with mode: 0644]
PWG1/comparison/AliTOFComparisonTask.h [new file with mode: 0644]
PWG1/comparison/AliTPCComparisonTask.cxx
PWG1/comparison/AliTPCComparisonTask.h
PWG1/comparison/AliTRDComparisonTask.cxx [new file with mode: 0644]
PWG1/comparison/AliTRDComparisonTask.h [new file with mode: 0644]
PWG1/comparison/runProofITSComparison.C [copied from PWG1/comparison/runProof.C with 54% similarity]
PWG1/comparison/runProofTOFComparison.C [copied from PWG1/comparison/runProof.C with 54% similarity]
PWG1/comparison/runProofTPCComparison.C [copied from PWG1/comparison/runProof.C with 66% similarity]
PWG1/comparison/runProofTRDComparison.C [moved from PWG1/comparison/runProof.C with 54% similarity]

diff --git a/PWG1/comparison/AF-v4-12.par b/PWG1/comparison/AF-v4-12.par
deleted file mode 100644 (file)
index 1cd13e9..0000000
Binary files a/PWG1/comparison/AF-v4-12.par and /dev/null differ
diff --git a/PWG1/comparison/AliITSComparisonTask.cxx b/PWG1/comparison/AliITSComparisonTask.cxx
new file mode 100644 (file)
index 0000000..2613cf1
--- /dev/null
@@ -0,0 +1,417 @@
+#include "TChain.h"
+#include "TTree.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TList.h"
+#include "TClonesArray.h"
+#include "TCanvas.h"
+#include "TStyle.h"
+#include "TLine.h"
+#include "TText.h"
+#include "TFile.h"
+#include "TBenchmark.h"
+
+#include "AliAnalysisTask.h"
+#include "AliAnalysisManager.h"
+
+#include "AliLog.h"
+#include "AliVEvent.h"
+#include "AliESDEvent.h"
+#include "AliMCEvent.h"
+#include "AliESDtrack.h"
+#include "AliTrackReference.h"
+#include "AliMCComparisonTrack.h"
+
+#include "AliITSComparisonTask.h"
+
+ClassImp(AliITSComparisonTask)
+
+AliITSComparisonTask::AliITSComparisonTask()
+  : AliAnalysisTaskSE("AliITSComaprisonTask"),
+    fListOfHistos(0),
+    hgood(0),
+    hfound(0),
+    hfake(0),
+    hp(0),
+    hl(0),
+    hpt(0),
+    htip(0),
+    he(0),
+    hep(0),
+    hgoodPhi(0),
+    hfoundPhi(0),
+    hlip(0)
+{
+  // Default constructor
+  AliInfo("Default constructor AliITSComparisonTask");
+  // 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());
+}
+
+
+
+AliITSComparisonTask::AliITSComparisonTask(const char* name)
+  : AliAnalysisTaskSE(name),
+    fListOfHistos(0),
+    hgood(0),
+    hfound(0),
+    hfake(0),
+    hp(0),
+    hl(0),
+    hpt(0),
+    htip(0),
+    he(0),
+    hep(0),
+    hgoodPhi(0),
+    hfoundPhi(0),
+    hlip(0)
+{
+  // Constructor
+  AliInfo("Constructor AliITSComparisonTask");
+  // 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());
+}
+
+
+
+void AliITSComparisonTask::UserCreateOutputObjects() 
+{
+  // Create histograms
+  // Called once
+  AliInfo("AliITSComparisonTask::UserCreateOutputObjects");
+  // Create output container
+  fListOfHistos = new TList();
+  
+  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.);
+  htip = new TH1F("htip", "Transverse impact parameter", 40, -800., 800.);
+  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, 0., 2.*TMath::Pi());
+  hfoundPhi = new TH1F("hfoundPhi", "Phi for found tracks", 90, 0., 2.*TMath::Pi());
+  hlip = new TH1F("hlip","Longitudinal impact parameter", 40, -800., 800.);
+  
+  
+  fListOfHistos->Add(hgood);
+  fListOfHistos->Add(hfound);
+  fListOfHistos->Add(hfake);
+  fListOfHistos->Add(hp);
+  fListOfHistos->Add(hl);
+  fListOfHistos->Add(hpt);
+  fListOfHistos->Add(htip);
+  fListOfHistos->Add(he);
+  fListOfHistos->Add(hep);
+  fListOfHistos->Add(hgoodPhi);
+  fListOfHistos->Add(hfoundPhi);
+  fListOfHistos->Add(hlip);
+}
+
+
+void AliITSComparisonTask::UserExec(Option_t *) 
+{
+  // Main loop
+  // Called for each event
+  
+    // MC information
+  AliMCEvent* mcEvent = MCEvent();
+  if (!mcEvent) {
+    Printf("ERROR: Could not retrieve MC event");
+    return;
+  }
+   
+  Int_t nt = 0;
+    
+  TClonesArray dummy("AliMCComparisonTrack",1000), *refs=&dummy;
+    
+  // 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) {
+      Printf("ERROR: Could not receive track %d (mc loop)", iTracks);
+      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;
+    if (TMath::Abs(vz) > 50.) continue; 
+    
+    // Loop over Track References
+    Bool_t LabelTPC = kFALSE;
+    AliTrackReference* trackRef = 0;
+    UInt_t ITSLayerMap = 0;
+    for (Int_t iTrackRef = 0; iTrackRef  < track->GetNumberOfTrackReferences(); iTrackRef++) {
+      trackRef = track->GetTrackReference(iTrackRef);
+      if(trackRef) {
+       Int_t detectorId = trackRef->DetectorId();
+       if (detectorId == AliTrackReference::kITS) {
+         Float_t Radius = trackRef->R();
+         if (Radius > 2.5 && Radius < 5.5)
+           ITSLayerMap |= 1;
+         else if (Radius > 5.5 && Radius < 8.5)
+           ITSLayerMap |= 1 << 1;
+         else if (Radius > 13.5 && Radius < 16.5)
+           ITSLayerMap |= 1 << 2;
+         else if (Radius > 22. && Radius < 26.)
+           ITSLayerMap |= 1 << 3;
+         else if (Radius > 36. && Radius < 41.)
+           ITSLayerMap |= 1 << 4;
+         else if (Radius > 42. && Radius < 46.)
+           ITSLayerMap |= 1 << 5;
+         else {
+           Printf("Wrong radius %f ", Radius);
+           return;
+         }
+
+       }
+       if (detectorId == AliTrackReference::kTPC) {        
+         LabelTPC = kTRUE;
+         break;
+       }
+      }      
+    }    // track references loop   
+    
+    // Skip tracks that passed not all ITS layers 
+    if (ITSLayerMap != 0x3F) continue;
+    
+    // "Good" tracks
+    if (LabelTPC) {
+      AliMCComparisonTrack* ref = new((*refs)[nt]) AliMCComparisonTrack();
+      ref->SetLabel(iTracks);
+      TParticle* particle = track->Particle();
+      Int_t pdg = particle->GetPdgCode();
+      ref->SetPDG(pdg); 
+      ref->SetPz(track->Pz());
+      Float_t Pt = track->Pt();
+      ref->SetPt(Pt);
+      hgood->Fill(Pt);
+      Float_t phig = track->Phi();
+      ref->SetPhi(phig);
+      hgoodPhi->Fill(phig);
+      nt++;  
+    }  
+  }    // track loop 
+  
+  
+  // ESD information  
+  AliVEvent* event = InputEvent();
+  if (!event) {
+    Printf("ERROR: Could not retrieve event");
+    return;
+  }
+    
+  AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
+  
+  Bool_t iFound;
+  Int_t nfound = 0;
+  Int_t nfake = 0;
+  Int_t nlost = 0;
+  
+  Int_t MCgoods = refs->GetEntriesFast();
+  
+  // Loop over all "good" MC tracks
+  for (Int_t k = 0; k < MCgoods; k++) {
+    AliMCComparisonTrack* ref = (AliMCComparisonTrack*)refs->UncheckedAt(k); 
+    Int_t MCLabel = ref->GetLabel();
+    Float_t ptg = ref->GetPt(); 
+    Float_t PhiG = ref->GetPhi();
+    iFound = kFALSE;
+    for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
+      AliESDtrack* track = esd->GetTrack(iTrack);  
+      if (!track) {
+        Printf("ERROR: Could not receive track %d", iTrack);
+        continue;
+      }
+      
+      //      if (! track->IsOn(AliESDtrack::kTPCrefit)) continue;
+      if (! track->IsOn(AliESDtrack::kITSrefit)) continue;
+      
+      Int_t Label =  track->GetLabel();
+      
+      if (MCLabel == TMath::Abs(Label)) {        
+        if (MCLabel == Label) {
+         nfound++;
+         hfound->Fill(ptg);
+         hfoundPhi->Fill(PhiG);
+       } 
+       else {
+         nfake++;
+         hfake->Fill(ptg);
+       }
+       iFound = kTRUE;
+
+       Float_t phi = track->Phi();
+       Double_t pt = track->Pt();
+       
+        Float_t lam = TMath::ATan2(track->Pz(),pt); 
+        Float_t pt_1 = 1/pt;
+       
+       Float_t phig = ref->GetPhi();
+        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 d,z; track->GetImpactParameters(d,z);
+        htip->Fill(10000.*d);
+        hlip->Fill(10000.*z);
+       
+       Float_t mom = pt/TMath::Cos(lam);
+        Float_t dedx = track->GetTPCsignal();
+        hep->Fill(mom, dedx, 1.);
+       
+       Int_t pdg = ref->GetPDG();
+        if (TMath::Abs(pdg )== 211) { // pions
+         if (mom > 0.4 && mom < 0.5) {
+            he->Fill(dedx, 1.);
+          }
+       }   
+          
+        break; 
+      }
+    }  
+    if (!iFound) {
+      nlost++;
+    } 
+  }
+     
+  Printf(" Results: " );
+  Printf(" Found %d Fake %d Lost %d ", nfound, nfake, nlost);
+  
+  refs->Clear();
+  
+  // Post output data.
+  PostData(1, fListOfHistos);
+}      
+
+//________________________________________________________________________
+void AliITSComparisonTask::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));  
+  hp = dynamic_cast<TH1F*>(fListOfHistos->At(3));
+  hl = dynamic_cast<TH1F*>(fListOfHistos->At(4));
+  hpt = dynamic_cast<TH1F*>(fListOfHistos->At(5));
+  htip = dynamic_cast<TH1F*>(fListOfHistos->At(6));
+  he = dynamic_cast<TH1F*>(fListOfHistos->At(7));
+  hep = dynamic_cast<TH2F*>(fListOfHistos->At(8));
+  hgoodPhi = dynamic_cast<TH1F*>(fListOfHistos->At(9));
+  hfoundPhi = dynamic_cast<TH1F*>(fListOfHistos->At(10));
+  hlip = dynamic_cast<TH1F*>(fListOfHistos->At(11));
+  
+  gStyle->SetOptStat(111110);
+  gStyle->SetOptFit(1);
+  
+  TCanvas* c1 = new TCanvas("c1", "", 0, 0, 700, 850);
+  
+  Int_t minc = 33;
+  
+  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* 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);
+  htip->SetFillColor(6); htip->SetXTitle("(micron)");
+  if (htip->GetEntries() < minc) htip->Draw(); else htip->Fit("gaus");
+  hlip->Draw("same"); 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., .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("AliITSComparison.root","RECREATE");
+  c1->Write();
+  c2->Write();
+  c3->Write();
+  fc.Close();
+}
diff --git a/PWG1/comparison/AliITSComparisonTask.h b/PWG1/comparison/AliITSComparisonTask.h
new file mode 100644 (file)
index 0000000..2a519c5
--- /dev/null
@@ -0,0 +1,48 @@
+#ifndef AliITSComparisonTask_cxx
+#define AliITSComparisonTask_cxx
+
+class TList;
+class TH1F;
+class TH2F;
+class TClonesArray;
+
+#ifdef __MAKECINT__
+#pragma link C++ class AliMCComparisonTrack+;
+#endif
+
+#include "AliAnalysisTaskSE.h"
+
+class AliITSComparisonTask: public AliAnalysisTaskSE
+{
+  public:
+    AliITSComparisonTask();
+    AliITSComparisonTask(const char* name);
+    virtual ~AliITSComparisonTask() {}
+    
+    virtual void UserCreateOutputObjects();
+    virtual void UserExec(Option_t *option);
+    virtual void Terminate(Option_t *);
+  
+  private:
+    TList* fListOfHistos;
+    
+    TH1F* hgood;
+    TH1F* hfound;
+    TH1F* hfake;
+    TH1F* hp;
+    TH1F* hl;
+    TH1F* hpt;
+    TH1F* htip;
+    TH1F* he;
+    TH2F* hep;
+    TH1F* hgoodPhi;
+    TH1F* hfoundPhi;
+    TH1F* hlip;
+    
+    AliITSComparisonTask(const AliITSComparisonTask&); // not implemented
+    AliITSComparisonTask& operator=(const AliITSComparisonTask&); // not implemented
+
+    ClassDef(AliITSComparisonTask, 1); // example of analysis 
+};
+
+#endif
index 05a4029..9dc9224 100644 (file)
@@ -1,5 +1,6 @@
 #include "AliMCComparisonTrack.h"
 
 AliMCComparisonTrack::AliMCComparisonTrack() :
-  MCLabel(0), MCPdg(0), TPCPx(0), TPCPy(0), TPCPz(0)
+  MCLabel(0), MCPdg(0), Pz(0), Pt(0), Phi(0), 
+  LocalX(0), LocalY(0), Z(0)
   {}
index 8fe575d..066a554 100644 (file)
@@ -14,29 +14,53 @@ class AliMCComparisonTrack: public TObject
       {MCLabel = Label;}
     void SetPDG(Int_t PDG) 
       {MCPdg = PDG;}
-    void SetTPCMomentum(Float_t Px, Float_t Py, Float_t Pz)
-      {TPCPx = Px; TPCPy = Py; TPCPz = Pz;}
+    void SetPz(Float_t PZ)
+      {Pz = PZ;}
+    void SetPt(Float_t Ptrans)
+      {Pt = Ptrans;}
+    void SetPhi(Float_t PhiAngle)
+      {Phi = PhiAngle;}
+    void SetLocalX(Float_t localX)
+      {LocalX =localX;}
+    void SetLocalY(Float_t localY)
+      {LocalY =localY;}  
+    void SetZ(Float_t Zcoor)
+      {Z = Zcoor;}
       
     Int_t GetLabel()
       {return MCLabel;}  
     Int_t GetPDG()
       {return MCPdg;}
-    Float_t GetTPCPx()
-      {return TPCPx;}
-    Float_t GetTPCPy()
-      {return TPCPy;}  
-    Float_t GetTPCPz()
-      {return TPCPz;}  
+    Float_t GetPz()
+      {return Pz;}  
+    Float_t GetPt()
+      {return Pt;}
+    Float_t GetPhi()
+      {return Phi;}
+    Float_t GetLocalX()
+      {return LocalX;}
+    Float_t GetLocalY()
+      {return LocalY;}      
+    Float_t GetZ()
+      {return Z;}
   
   private:
       // track label
     Int_t MCLabel;
-      // particle PDG
+      // PDG particle code
     Int_t MCPdg;
-      // momentum at the entrance of the TPC
-    Float_t TPCPx;
-    Float_t TPCPy;
-    Float_t TPCPz;
+      // z-component of momentum
+    Float_t Pz;
+      // component of momentum in the transverse plane
+    Float_t Pt; 
+      // phi angle of the particle direction
+    Float_t Phi;
+      // local x-coordinate
+    Float_t LocalX;
+      // local y-coordinate
+    Float_t LocalY;
+      // z-coordinate
+    Float_t Z;
     
     ClassDef(AliMCComparisonTrack, 1); 
 };
diff --git a/PWG1/comparison/AliTOFComparisonTask.cxx b/PWG1/comparison/AliTOFComparisonTask.cxx
new file mode 100644 (file)
index 0000000..078e211
--- /dev/null
@@ -0,0 +1,365 @@
+#include "TChain.h"
+#include "TTree.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TList.h"
+#include "TClonesArray.h"
+#include "TCanvas.h"
+#include "TStyle.h"
+#include "TLine.h"
+#include "TText.h"
+#include "TFile.h"
+#include "TBenchmark.h"
+
+#include "AliAnalysisTask.h"
+#include "AliAnalysisManager.h"
+
+#include "AliLog.h"
+#include "AliVEvent.h"
+#include "AliESDEvent.h"
+#include "AliMCEvent.h"
+#include "AliESDtrack.h"
+#include "AliTrackReference.h"
+#include "AliMCComparisonTrack.h"
+
+#include "AliTOFComparisonTask.h"
+
+ClassImp(AliTOFComparisonTask)
+
+AliTOFComparisonTask::AliTOFComparisonTask()
+  : AliAnalysisTaskSE("AliTOFComparisonTask"),
+    fListOfHistos(0),
+    hgood(0),
+    hfound(0),
+    hfake(0),
+    hgoodPhi(0),
+    hfoundPhi(0),
+    hgoodl(0),
+    hfakel(0),
+    hfoundl(0)
+{
+  // Default constructor
+  AliInfo("Defaulut Constructor AliTOFComparisonTask");
+  // 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());
+}
+
+AliTOFComparisonTask::AliTOFComparisonTask(const char* name)
+  : AliAnalysisTaskSE(name),
+    fListOfHistos(0),
+    hgood(0),
+    hfound(0),
+    hfake(0),
+    hgoodPhi(0),
+    hfoundPhi(0),
+    hgoodl(0),
+    hfakel(0),
+    hfoundl(0)
+{
+  // Constructor
+  AliInfo("Constructor AliTOFComparisonTask");
+  // 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());
+}
+
+
+
+void AliTOFComparisonTask::UserCreateOutputObjects() 
+{
+  // Create histograms
+  // Called once
+  AliInfo("AliTOFComparisonTask::UserCreateOutputObjects");
+  // Create output container
+  fListOfHistos = new TList();
+  
+  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);
+  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());
+  hgoodl = new TH1F("hgoodl", "Good tracks", 30, -1., 1.);
+  hfakel = new TH1F("hfakel", "Mismatched tracks", 30, -1., 1.);
+  hfoundl = new TH1F("hfoundl", "Matched tracks", 30, -1., 1.);
+  
+  fListOfHistos->Add(hgood);
+  fListOfHistos->Add(hfound);
+  fListOfHistos->Add(hfake);
+  fListOfHistos->Add(hgoodPhi);
+  fListOfHistos->Add(hfoundPhi);
+  fListOfHistos->Add(hgoodl);
+  fListOfHistos->Add(hfakel);
+  fListOfHistos->Add(hfoundl);
+}
+
+
+void AliTOFComparisonTask::UserExec(Option_t *) 
+{
+  // Main loop
+  // Called for each event
+
+  // MC information
+  AliMCEvent* mcEvent = MCEvent();
+  if (!mcEvent) {
+    Printf("ERROR: Could not retrieve MC event");
+    return;
+  }
+  Int_t nt = 0;
+    
+  TClonesArray dummy("AliMCComparisonTrack",1000), *refs=&dummy;
+    
+  // 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) {
+      Printf("ERROR: Could not receive track %d (mc loop)", iTracks);
+      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;
+    if (TMath::Abs(vz) > 50.) continue; 
+    
+    if (TMath::Abs(vx) > 0.1) continue;
+    if (TMath::Abs(vy) > 0.1) continue;
+       
+       
+    // Loop over Track References
+    Bool_t LabelTPC = kFALSE;
+    Bool_t LabelTOF = kFALSE;
+    AliTrackReference* trackRef = 0;
+    UInt_t ITSLayerMap = 0;
+    for (Int_t iTrackRef = 0; iTrackRef  < track->GetNumberOfTrackReferences(); iTrackRef++) {
+      trackRef = track->GetTrackReference(iTrackRef);
+      if(trackRef) {
+        Int_t detectorId = trackRef->DetectorId();
+       if (detectorId == AliTrackReference::kITS) {
+          Float_t Radius = trackRef->R();
+         if (Radius > 2.5 && Radius < 5.5)
+            ITSLayerMap |= 1;
+          else if (Radius > 5.5 && Radius < 8.5)
+            ITSLayerMap |= 1 << 1;
+          else if (Radius > 13.5 && Radius < 16.5)
+            ITSLayerMap |= 1 << 2;
+          else if (Radius > 22. && Radius < 26.)
+            ITSLayerMap |= 1 << 3;
+          else if (Radius > 36. && Radius < 41.)
+            ITSLayerMap |= 1 << 4;
+          else if (Radius > 42. && Radius < 46.)
+            ITSLayerMap |= 1 << 5;
+          else {
+            Printf("Wrong radius %f ", Radius);
+            return;
+          }
+
+        }
+       if (detectorId == AliTrackReference::kTPC) {        
+         LabelTPC = kTRUE;
+       }
+       if (detectorId == AliTrackReference::kTOF) {
+         LabelTOF = kTRUE;
+       }
+      }      
+    }    // track references loop   
+    
+    // Skip tracks that passed not all ITS layers 
+    if (ITSLayerMap != 0x3F) continue;
+    
+    // "Good" tracks
+    if (LabelTPC && LabelTOF) {
+      AliMCComparisonTrack* ref = new((*refs)[nt]) AliMCComparisonTrack();
+      ref->SetLabel(iTracks);
+      TParticle* particle = track->Particle();
+      Int_t pdg = particle->GetPdgCode();
+      ref->SetPDG(pdg); 
+      ref->SetPz(track->Pz());
+      Float_t Pt = track->Pt();
+      ref->SetPt(Pt);
+      hgood->Fill(Pt);
+      Float_t phig = track->Phi();
+      ref->SetPhi(phig);
+      hgoodPhi->Fill(phig);
+      Double_t tgl = track->Pz()/Pt;
+      hgoodl->Fill(tgl);
+      nt++;  
+    }      
+  }    // track loop 
+  
+  // ESD information  
+  AliVEvent* event = InputEvent();
+  if (!event) {
+    Printf("ERROR: Could not retrieve event");
+    return;
+  }
+    
+  AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
+  
+  Bool_t iFound;
+  Int_t nfound = 0;
+  Int_t nfake = 0;
+  Int_t nlost = 0;
+  
+  Int_t MCgoods = refs->GetEntriesFast();
+    
+  // Loop over all "good" MC tracks
+  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 = ref->GetPt(); 
+    Float_t Phig = ref->GetPhi();
+    iFound = kFALSE;
+    Int_t iTrack;
+    AliESDtrack* track = 0;
+    for (iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
+      track = esd->GetTrack(iTrack);  
+      if (!track) {
+        Printf("ERROR: Could not receive track %d", iTrack);
+        continue;
+      }
+      
+      Int_t Label =  track->GetLabel();
+      if (MCLabel != TMath::Abs(Label)) continue;
+  
+      if (track->GetTOFsignal() <= 0.) continue;
+      
+      Double_t tgl = ref->GetPz()/ref->GetPt();      
+      Int_t TOFLabel[3];  
+      track->GetTOFLabel(TOFLabel);            
+      if (TOFLabel[0] != MCLabel)
+        if (TOFLabel[1] != MCLabel)
+         if (TOFLabel[2] != MCLabel) {
+           hfake->Fill(ptg); 
+           hfakel->Fill(tgl);
+           nfake++;
+           break;
+         }
+      hfound->Fill(ptg); 
+      hfoundl->Fill(tgl);
+      hfoundPhi->Fill(Phig);
+      nfound++;          
+      break;      
+    } 
+    if(iTrack == esd->GetNumberOfTracks()) {
+      Printf("Not matched: MCLabel %d ", MCLabel);
+      nlost++;
+      if (track) {
+        Printf("  kITSout %d kTPCout %d kTRDout %d kTIME %d",
+              (track->GetStatus()&AliESDtrack::kITSout),
+              (track->GetStatus()&AliESDtrack::kTPCout),
+              (track->GetStatus()&AliESDtrack::kTRDout),
+              (track->GetStatus()&AliESDtrack::kTIME));
+      }
+      else
+        Printf(" No ESD track !");
+    }  
+  }    
+     
+  Printf(" Results: " );
+  Printf(" Good %d Found %d Fake %d Lost %d", nt, nfound, nfake, nlost);
+    
+  refs->Clear();
+
+  // Post output data.
+  PostData(1, fListOfHistos);
+
+}      
+
+
+void AliTOFComparisonTask::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));  
+  hgoodPhi = dynamic_cast<TH1F*>(fListOfHistos->At(3));
+  hfoundPhi = dynamic_cast<TH1F*>(fListOfHistos->At(4));
+  hgoodl = dynamic_cast<TH1F*>(fListOfHistos->At(5));
+  hfakel = dynamic_cast<TH1F*>(fListOfHistos->At(6));
+  hfoundl = dynamic_cast<TH1F*>(fListOfHistos->At(7));
+  
+  gStyle->SetOptStat(111110);
+  gStyle->SetOptFit(1); 
+  
+  TCanvas* c1 = new TCanvas("c1","",0,0,600,900);
+  
+  TPad* p1 = new TPad("p1", "", 0., 0.5, 1., 1.); p1->Draw();  
+  p1->cd(); p1->SetFillColor(42); p1->SetFrameFillColor(10);
+  
+  hfound->Sumw2(); hgood->Sumw2(); hfake->Sumw2();
+  TH1F* hgp = new TH1F("hgp", " ", 34, 0.2, 7.0);
+  TH1F* hfp = new TH1F("hfp", "Probability of mismatching", 34, 0.2, 7.0);
+  hgp->Divide(hfound, hgood, 1., 1., "b");
+  hfp->Divide(hfake, hgood, 1., 1., "b");
+  hgp->SetLineColor(4); hgp->SetLineWidth(2);
+  hgp->SetMaximum(1.4);
+  hgp->SetYTitle("Tracking efficiency");
+  hgp->SetXTitle("Pt (GeV/c)");
+  hfp->SetFillColor(1); hfp->SetFillStyle(3013); hfp->SetLineWidth(2);
+  hfp->SetLineColor(2);
+  
+  hgp->Draw();
+  hfp->Draw("histsame");
+  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");
+  c1->cd();
+  
+  TPad* p2 = new TPad("p2", "", 0., 0., 1., 0.5); p2->Draw();
+  p2->cd(); p2->SetFillColor(42); p2->SetFrameFillColor(10);  
+  
+  hfoundl->Sumw2(); hgoodl->Sumw2(); hfakel->Sumw2();
+  TH1F* hgl = new TH1F("hgl", "", 30, -1., 1.);
+  TH1F* hfl = new TH1F("hfl", "Probability of mismatching", 30, -1., 1.);
+  hgl->Divide(hfoundl, hgoodl, 1., 1., "b");
+  hfl->Divide(hfakel, hgoodl, 1., 1., "b");
+  hgl->SetLineColor(4); hgl->SetLineWidth(2);
+  hgl->SetMaximum(1.4);
+  hgl->SetYTitle("Tracking efficiency");
+  hgl->SetXTitle("Tan(lambda)");
+  hfl->SetFillColor(1); hfl->SetFillStyle(3013); hfl->SetLineWidth(2); 
+  hfl->SetLineColor(2);
+  
+  hgl->Draw();
+  hfl->Draw("histsame");
+  TLine *line3 = new TLine(-1,1.0,1,1.0); line3->SetLineStyle(4);
+  line3->Draw("same");
+  TLine *line4 = new TLine(-1,0.9,1,0.9); line4->SetLineStyle(4);
+  line4->Draw("same");
+  
+  c1->Update();
+  
+  TCanvas* c2 = new TCanvas("c2", "", 10, 10, 510, 510);
+  c2->SetFillColor(42); c2->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("AliTOFComparison.root","RECREATE");
+  c1->Write();
+  c2->Write();
+  fc.Close();  
+}
diff --git a/PWG1/comparison/AliTOFComparisonTask.h b/PWG1/comparison/AliTOFComparisonTask.h
new file mode 100644 (file)
index 0000000..215c2fb
--- /dev/null
@@ -0,0 +1,44 @@
+#ifndef AliTOFComparisonTask_cxx
+#define AliTOFComparisonTask_cxx
+
+class TList;
+class TH1F;
+class TH2F;
+class TClonesArray;
+
+#ifdef __MAKECINT__
+#pragma link C++ class AliMCComparisonTrack+;
+#endif
+
+#include "AliAnalysisTaskSE.h"
+
+class AliTOFComparisonTask: public AliAnalysisTaskSE
+{
+  public:
+    AliTOFComparisonTask();
+    AliTOFComparisonTask(const char* name);
+    virtual ~AliTOFComparisonTask() {}
+    
+    virtual void UserCreateOutputObjects();
+    virtual void UserExec(Option_t *option);
+    virtual void Terminate(Option_t *);
+  
+  private:
+    TList* fListOfHistos;
+    
+    TH1F* hgood;
+    TH1F* hfound;
+    TH1F* hfake;
+    TH1F* hgoodPhi;
+    TH1F* hfoundPhi;
+    TH1F* hgoodl;
+    TH1F* hfakel;
+    TH1F* hfoundl;
+    
+    AliTOFComparisonTask(const AliTOFComparisonTask&); // not implemented
+    AliTOFComparisonTask& operator=(const AliTOFComparisonTask&); // not implemented
+
+    ClassDef(AliTOFComparisonTask, 1); // example of analysis 
+};
+
+#endif
index f303ba6..b0cb745 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();
 
 }
index 8838d2c..a32c4b6 100644 (file)
@@ -15,7 +15,8 @@ class TClonesArray;
 class AliTPCComparisonTask: public AliAnalysisTaskSE
 {
   public:
-    AliTPCComparisonTask(const char* name = "AliTPCComparisonTask");
+    AliTPCComparisonTask();
+    AliTPCComparisonTask(const char* name);
     virtual ~AliTPCComparisonTask() {}
     
     virtual void UserCreateOutputObjects();
diff --git a/PWG1/comparison/AliTRDComparisonTask.cxx b/PWG1/comparison/AliTRDComparisonTask.cxx
new file mode 100644 (file)
index 0000000..d7d0dc7
--- /dev/null
@@ -0,0 +1,436 @@
+#include "TChain.h"
+#include "TTree.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TList.h"
+#include "TClonesArray.h"
+#include "TCanvas.h"
+#include "TStyle.h"
+#include "TLine.h"
+#include "TText.h"
+#include "TFile.h"
+#include "TBenchmark.h"
+
+#include "AliAnalysisTask.h"
+#include "AliAnalysisManager.h"
+
+#include "AliLog.h"
+#include "AliVEvent.h"
+#include "AliESDEvent.h"
+#include "AliMCEvent.h"
+#include "AliESDtrack.h"
+#include "AliTrackReference.h"
+#include "AliMCComparisonTrack.h"
+
+#include "AliTRDComparisonTask.h"
+
+ClassImp(AliTRDComparisonTask)
+
+AliTRDComparisonTask::AliTRDComparisonTask()
+  : AliAnalysisTaskSE("AliTRDComparisonTask"),
+    fListOfHistos(0),
+    hgood(0),
+    hfound(0),
+    hfake(0),
+    hp(0),
+    hl(0),
+    hpt(0),
+    hmpt(0),
+    he(0),
+    hep(0),
+    hgoodPhi(0),
+    hfoundPhi(0),
+    hz(0),
+    hc(0)
+{
+  // Default constructor
+  AliInfo("Default constructor AliTRDComparisonTask");
+  // 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());
+}
+
+AliTRDComparisonTask::AliTRDComparisonTask(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),
+    hz(0),
+    hc(0)
+{
+  // Constructor
+  AliInfo("Constructor AliTRDComparisonTask");
+  // 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());
+}
+
+
+
+void AliTRDComparisonTask::UserCreateOutputObjects() 
+{
+  // Create histograms
+  // Called once
+  AliInfo("AliTRDComparisonTask::UserCreateOutputObjects");
+  // Create output container
+  fListOfHistos = new TList();
+  
+  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", "Y and Z resolution", 30, -30., 30.);
+  he = new TH1F("he", "dE/dX for pions with 0.4<p<0.5 GeV/c", 50, 0., 1000.);
+  hep = new TH2F("hep", "dE/dX vs momentum", 50, 0., 2., 50, 0., 2000.);
+  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());
+  hz = new TH1F("hz", "Z resolution", 30, -30., 30.);
+  hc = new TH1F("hc", "Number of the assigned clusters", 25, 110., 135.);
+  
+  fListOfHistos->Add(hgood);
+  fListOfHistos->Add(hfound);
+  fListOfHistos->Add(hfake);
+  fListOfHistos->Add(hp);
+  fListOfHistos->Add(hl);
+  fListOfHistos->Add(hpt);
+  fListOfHistos->Add(hmpt);
+  fListOfHistos->Add(he);
+  fListOfHistos->Add(hep);
+  fListOfHistos->Add(hgoodPhi);
+  fListOfHistos->Add(hfoundPhi);
+  fListOfHistos->Add(hz);
+  fListOfHistos->Add(hc);  
+}
+
+
+void AliTRDComparisonTask::UserExec(Option_t *) 
+{
+
+  // Main loop
+  // Called for each event
+  
+  // MC information
+  AliMCEvent* mcEvent = MCEvent();
+  if (!mcEvent) {
+    Printf("ERROR: Could not retrieve MC event");
+    return;
+  }
+  
+  Int_t nt = 0;
+    
+  TClonesArray dummy("AliMCComparisonTrack",1000), *refs=&dummy;
+    
+  // 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) {
+      Printf("ERROR: Could not receive track %d (mc loop)", iTracks);
+      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;
+    if (TMath::Abs(vz) > 50.) continue; 
+    
+    // Check TRD information       
+    Int_t nRefs = track->GetNumberOfTrackReferences();
+    if (nRefs <= 1) continue;
+
+    AliTrackReference* trackRef0 = 0;
+    for (Int_t iTrackRef = 0; iTrackRef < nRefs; iTrackRef++) {
+      trackRef0 = track->GetTrackReference(iTrackRef);
+      if(trackRef0) {
+        Int_t detectorId = trackRef0->DetectorId();
+        if (detectorId == AliTrackReference::kTRD) {
+         break;
+        }
+       trackRef0 = 0;
+      }
+    }
+    if (!trackRef0) continue;
+
+    if (trackRef0->LocalX() > 300.) continue;
+    
+    AliTrackReference* trackRef = 0;
+    for (Int_t iTrackRef = nRefs - 1; iTrackRef >= 0; --iTrackRef) {
+      trackRef = track->GetTrackReference(iTrackRef);
+      if (!trackRef) continue;
+      if (trackRef->LocalX() > 363. &&
+         trackRef->DetectorId() == AliTrackReference::kTRD) break;       
+      trackRef = 0;           
+    }
+    if (!trackRef) continue;
+     
+    if (TMath::Abs(trackRef0->Alpha() - trackRef->Alpha()) > 1e-5) continue;
+    
+    // Check TPC information
+    Bool_t LabelTPC = kFALSE;
+    for (Int_t iTrackRef = 0; iTrackRef  < nRefs; iTrackRef++) {
+      trackRef = track->GetTrackReference(iTrackRef);
+      if(trackRef) {
+       Int_t detectorId = trackRef->DetectorId();
+       if (detectorId == AliTrackReference::kTPC) {        
+         LabelTPC = kTRUE;
+         break;
+       }
+      }      
+    }    // track references loop   
+    
+    // "Good" tracks
+    if (LabelTPC) {
+      AliMCComparisonTrack* ref = new((*refs)[nt]) AliMCComparisonTrack();
+      ref->SetLabel(iTracks);
+      TParticle* particle = track->Particle();
+      Int_t pdg = particle->GetPdgCode();
+      ref->SetPDG(pdg); 
+      ref->SetPz(trackRef->Pz());
+      Float_t Pt = trackRef->Pt();
+      ref->SetPt(Pt);
+      hgood->Fill(Pt);
+      Float_t phig = trackRef->Phi();
+      ref->SetPhi(phig);
+      hgoodPhi->Fill(phig);
+      ref->SetLocalX(trackRef->LocalX());
+      ref->SetLocalY(trackRef->LocalY());
+      ref->SetZ(trackRef->Z());
+      nt++;  
+    }  
+  }    // track loop 
+
+  
+  // ESD information  
+  AliVEvent* event = InputEvent();
+  if (!event) {
+    Printf("ERROR: Could not retrieve event");
+    return;
+  }
+    
+  AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
+  
+  Bool_t iFound;
+  Int_t nfound = 0;
+  Int_t nfake = 0;
+  Int_t nlost = 0;
+  
+  Int_t MCgoods = refs->GetEntriesFast();
+    
+  // Loop over all "good" MC tracks
+  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 = ref->GetPt(); 
+    Float_t PhiG = ref->GetPhi();
+    iFound = kFALSE;
+    for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
+      AliESDtrack* track = esd->GetTrack(iTrack);  
+      if (!track) {
+        Printf("ERROR: Could not receive track %d", iTrack);
+        continue;
+      }
+      
+      if (! track->IsOn(AliESDtrack::kTRDout)) continue;
+      const AliExternalTrackParam * outorig = track->GetOuterParam();
+      if (!outorig) continue;
+      
+      Int_t Label =  track->GetTRDLabel();      
+
+      if (MCLabel == TMath::Abs(Label)) {        
+        if (MCLabel == Label) {
+         nfound++;
+         hfound->Fill(ptg);
+         hfoundPhi->Fill(PhiG);
+       } 
+       else {
+         nfake++;
+         hfake->Fill(ptg);
+       }
+       iFound = kTRUE;
+       
+       AliExternalTrackParam out(*outorig);
+
+        Double_t bz = esd->GetMagneticField();
+        Double_t xg = ref->GetLocalX();
+        out.PropagateTo(xg,bz);
+       
+        Float_t phi = TMath::ASin(out.GetSnp()) + out.GetAlpha();
+        if (phi < 0) phi += 2*TMath::Pi();
+        if (phi >= 2*TMath::Pi()) phi -= 2*TMath::Pi();
+        Float_t phig = ref->GetPhi();
+        hp->Fill((phi - phig)*1000.);
+       
+        Float_t lam = TMath::ATan(out.GetTgl()); 
+        Float_t lamg = TMath::ATan2(ref->GetPz(), ptg);
+        hl->Fill((lam - lamg)*1000.);
+       
+        hc->Fill(track->GetTRDclusters(0));
+       
+        Float_t pt_1 = out.OneOverPt();
+        hpt->Fill((pt_1 - 1/ptg)/(1/ptg)*100.);
+
+        Float_t y = out.GetY();
+        Float_t yg = ref->GetLocalY();
+        hmpt->Fill(10*(y - yg));
+
+        Float_t z = out.GetZ();
+        Float_t zg = ref->GetZ();
+        hz->Fill(10.*(z - zg));
+
+        Float_t mom = 1./(pt_1*TMath::Cos(lam));
+        Float_t dedx = track->GetTRDsignal();
+       Printf (" dedx = %f ", dedx);
+        hep->Fill(mom, dedx, 1.);
+       
+        Int_t pdg = ref->GetPDG();
+        if (TMath::Abs(pdg)==211) //pions
+         if (mom>0.4 && mom<0.5) he->Fill(dedx,1.);
+          
+        break; 
+      }
+    }  
+    if (!iFound) {
+      nlost++;
+    } 
+  }
+     
+  Printf(" Results: " );
+  Printf(" Good %d Found %d Fake %d Lost %d ", nt, nfound, nfake, nlost);
+    
+  refs->Clear();
+
+  // Post output data.
+  PostData(1, fListOfHistos);  
+}      
+
+
+void AliTRDComparisonTask::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));  
+  hp = dynamic_cast<TH1F*>(fListOfHistos->At(3));
+  hl = dynamic_cast<TH1F*>(fListOfHistos->At(4));
+  hpt = dynamic_cast<TH1F*>(fListOfHistos->At(5));
+  hmpt = dynamic_cast<TH1F*>(fListOfHistos->At(6));
+  he = dynamic_cast<TH1F*>(fListOfHistos->At(7));
+  hep = dynamic_cast<TH2F*>(fListOfHistos->At(8));
+  hgoodPhi = dynamic_cast<TH1F*>(fListOfHistos->At(9));
+  hfoundPhi = dynamic_cast<TH1F*>(fListOfHistos->At(10));
+  hz = dynamic_cast<TH1F*>(fListOfHistos->At(11));
+  hc = dynamic_cast<TH1F*>(fListOfHistos->At(12));
+
+  gStyle->SetOptStat(111110);
+  gStyle->SetOptFit(1);
+  
+  TCanvas* c1 = new TCanvas("c1", "", 0, 0, 700, 850);
+
+  Int_t minc = 33;
+  
+  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* 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("(mm)");
+  if (hmpt->GetEntries() < minc) hmpt->Draw(); else hmpt->Fit("gaus");
+  hz->Draw("same"); 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->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., .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->SetMarkerStyle(8); hep->SetMarkerSize(0.4);
+  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, 0., 2.*TMath::Pi());
+  hgphi->Divide(hfoundPhi, hgoodPhi, 1., 1., "B"); 
+  hgphi->SetYTitle("Tracking efficiency");
+  hgphi->SetXTitle("Phi (rad)");
+  hgphi->Draw();
+
+  TFile fc("AliTRDComparison.root","RECREATE");
+  c1->Write();
+  c2->Write();
+  c3->Write();
+  fc.Close();  
+}
diff --git a/PWG1/comparison/AliTRDComparisonTask.h b/PWG1/comparison/AliTRDComparisonTask.h
new file mode 100644 (file)
index 0000000..7025263
--- /dev/null
@@ -0,0 +1,49 @@
+#ifndef AliTRDComparisonTask_cxx
+#define AliTRDComparisonTask_cxx
+
+class TList;
+class TH1F;
+class TH2F;
+class TClonesArray;
+
+#ifdef __MAKECINT__
+#pragma link C++ class AliMCComparisonTrack+;
+#endif
+
+#include "AliAnalysisTaskSE.h"
+
+class AliTRDComparisonTask: public AliAnalysisTaskSE
+{
+  public:
+    AliTRDComparisonTask();
+    AliTRDComparisonTask(const char* name);
+    virtual ~AliTRDComparisonTask() {}
+    
+    virtual void UserCreateOutputObjects();
+    virtual void UserExec(Option_t *option);
+    virtual void Terminate(Option_t *);
+  
+  private:
+    TList* fListOfHistos;
+    
+    TH1F* hgood;
+    TH1F* hfound;
+    TH1F* hfake;
+    TH1F* hp;
+    TH1F* hl;
+    TH1F* hpt;
+    TH1F* hmpt;
+    TH1F* he;
+    TH2F* hep;
+    TH1F* hgoodPhi;
+    TH1F* hfoundPhi;
+    TH1F* hz;
+    TH1F* hc;
+    
+    AliTRDComparisonTask(const AliTRDComparisonTask&); // not implemented
+    AliTRDComparisonTask& operator=(const AliTRDComparisonTask&); // not implemented
+
+    ClassDef(AliTRDComparisonTask, 1); // example of analysis 
+};
+
+#endif
similarity index 54%
copy from PWG1/comparison/runProof.C
copy to PWG1/comparison/runProofITSComparison.C
index a3192d9..757c4cf 100644 (file)
@@ -1,14 +1,14 @@
-void runProof(const char * dataset = "/PWG0/COMMON/run30000X_10TeV_0.5T", Long64_t nentries=1234567890, Long64_t firstentry=0)
+void runProofITSComparison(const char * dataset ="/PWG0/COMMON/run30000X_10TeV_0.5T",Long64_t nentries=1000, Long64_t firstentry=0)
 {
   // Connect to Proof
   TProof::Open("lxb6046");
 
   // Upload and enable packages: please use the correct version!
-  gProof->UploadPackage("AF-v4-12");
-  gProof->EnablePackage("AF-v4-12");
+  gProof->UploadPackage("AF-v4-14");
+  gProof->EnablePackage("AF-v4-14");
 
   // Create the analysis manager
-  AliAnalysisManager *mgr = new AliAnalysisManager("AliTPCComparison");
+  AliAnalysisManager *mgr = new AliAnalysisManager("AliITSComparison");
 
   AliVEventHandler* esdH = new AliESDInputHandler();
   mgr->SetInputEventHandler(esdH);
@@ -19,15 +19,18 @@ void runProof(const char * dataset = "/PWG0/COMMON/run30000X_10TeV_0.5T", Long64
 
   // Create task
   gProof->Load("AliMCComparisonTrack.cxx++g");
-  gProof->Load("AliTPCComparisonTask.cxx++g");
-  AliAnalysisTask *task = new AliTPCComparisonTask("AliTPCComparisonTask");
+  gProof->Load("AliITSComparisonTask.cxx++g");
+  AliAnalysisTask *task = new AliITSComparisonTask("AliITSComparisonTask");
 
   // Add task
   mgr->AddTask(task);
 
   // Create containers for input/output
-  AliAnalysisDataContainer *cinput = mgr->CreateContainer("cchain", TChain::Class(), AliAnalysisManager::kInputContainer);
-  AliAnalysisDataContainer *coutput = mgr->CreateContainer("coutput", TList::Class(), AliAnalysisManager::kOutputContainer, "Hist.root");
+  AliAnalysisDataContainer* cinput = 
+    mgr->CreateContainer("cchain", TChain::Class(), AliAnalysisManager::kInputContainer);
+  AliAnalysisDataContainer* coutput = 
+    mgr->CreateContainer("coutput", TList::Class(), 
+    AliAnalysisManager::kOutputContainer, "AliITSComparisonHist.root");
 
   // Connect input/output
   mgr->ConnectInput(task, 0, cinput);
@@ -44,3 +47,4 @@ void runProof(const char * dataset = "/PWG0/COMMON/run30000X_10TeV_0.5T", Long64
 
   mgr->StartAnalysis("proof",dataset,nentries,firstentry);
 }
+
similarity index 54%
copy from PWG1/comparison/runProof.C
copy to PWG1/comparison/runProofTOFComparison.C
index a3192d9..ac2f506 100644 (file)
@@ -1,14 +1,14 @@
-void runProof(const char * dataset = "/PWG0/COMMON/run30000X_10TeV_0.5T", Long64_t nentries=1234567890, Long64_t firstentry=0)
+void runProofTOFComparison(const char * dataset="/PWG0/COMMON/run30000X_10TeV_0.5T",Long64_t nentries=1000, Long64_t firstentry=0)
 {
   // Connect to Proof
   TProof::Open("lxb6046");
 
   // Upload and enable packages: please use the correct version!
-  gProof->UploadPackage("AF-v4-12");
-  gProof->EnablePackage("AF-v4-12");
+  gProof->UploadPackage("AF-v4-14");
+  gProof->EnablePackage("AF-v4-14");
 
   // Create the analysis manager
-  AliAnalysisManager *mgr = new AliAnalysisManager("AliTPCComparison");
+  AliAnalysisManager *mgr = new AliAnalysisManager("AliTOFComparison");
 
   AliVEventHandler* esdH = new AliESDInputHandler();
   mgr->SetInputEventHandler(esdH);
@@ -19,15 +19,18 @@ void runProof(const char * dataset = "/PWG0/COMMON/run30000X_10TeV_0.5T", Long64
 
   // Create task
   gProof->Load("AliMCComparisonTrack.cxx++g");
-  gProof->Load("AliTPCComparisonTask.cxx++g");
-  AliAnalysisTask *task = new AliTPCComparisonTask("AliTPCComparisonTask");
+  gProof->Load("AliTOFComparisonTask.cxx++g");
+  AliAnalysisTask *task = new AliTOFComparisonTask("AliTOFComparisonTask");
 
   // Add task
   mgr->AddTask(task);
 
   // Create containers for input/output
-  AliAnalysisDataContainer *cinput = mgr->CreateContainer("cchain", TChain::Class(), AliAnalysisManager::kInputContainer);
-  AliAnalysisDataContainer *coutput = mgr->CreateContainer("coutput", TList::Class(), AliAnalysisManager::kOutputContainer, "Hist.root");
+  AliAnalysisDataContainer* cinput = 
+    mgr->CreateContainer("cchain", TChain::Class(), AliAnalysisManager::kInputContainer);
+  AliAnalysisDataContainer* coutput = 
+    mgr->CreateContainer("coutput", TList::Class(), 
+    AliAnalysisManager::kOutputContainer, "AliTOFComparisonHist.root");
 
   // Connect input/output
   mgr->ConnectInput(task, 0, cinput);
@@ -44,3 +47,4 @@ void runProof(const char * dataset = "/PWG0/COMMON/run30000X_10TeV_0.5T", Long64
 
   mgr->StartAnalysis("proof",dataset,nentries,firstentry);
 }
+
similarity index 66%
copy from PWG1/comparison/runProof.C
copy to PWG1/comparison/runProofTPCComparison.C
index a3192d9..5ce0986 100644 (file)
@@ -1,11 +1,11 @@
-void runProof(const char * dataset = "/PWG0/COMMON/run30000X_10TeV_0.5T", Long64_t nentries=1234567890, Long64_t firstentry=0)
+void runProofTPCComparison(const char * dataset = "/PWG0/COMMON/run30000X_10TeV_0.5T",Long64_t nentries=1000, Long64_t firstentry=0)
 {
   // Connect to Proof
   TProof::Open("lxb6046");
 
   // Upload and enable packages: please use the correct version!
-  gProof->UploadPackage("AF-v4-12");
-  gProof->EnablePackage("AF-v4-12");
+  gProof->UploadPackage("AF-v4-14");
+  gProof->EnablePackage("AF-v4-14");
 
   // Create the analysis manager
   AliAnalysisManager *mgr = new AliAnalysisManager("AliTPCComparison");
@@ -26,8 +26,11 @@ void runProof(const char * dataset = "/PWG0/COMMON/run30000X_10TeV_0.5T", Long64
   mgr->AddTask(task);
 
   // Create containers for input/output
-  AliAnalysisDataContainer *cinput = mgr->CreateContainer("cchain", TChain::Class(), AliAnalysisManager::kInputContainer);
-  AliAnalysisDataContainer *coutput = mgr->CreateContainer("coutput", TList::Class(), AliAnalysisManager::kOutputContainer, "Hist.root");
+  AliAnalysisDataContainer *cinput = 
+    mgr->CreateContainer("cchain", TChain::Class(), AliAnalysisManager::kInputContainer);
+  AliAnalysisDataContainer *coutput = 
+    mgr->CreateContainer("coutput", TList::Class(), 
+    AliAnalysisManager::kOutputContainer, "AliTPCComparisonHist.root");
 
   // Connect input/output
   mgr->ConnectInput(task, 0, cinput);
@@ -44,3 +47,4 @@ void runProof(const char * dataset = "/PWG0/COMMON/run30000X_10TeV_0.5T", Long64
 
   mgr->StartAnalysis("proof",dataset,nentries,firstentry);
 }
+
similarity index 54%
rename from PWG1/comparison/runProof.C
rename to PWG1/comparison/runProofTRDComparison.C
index a3192d9..47f6f44 100644 (file)
@@ -1,14 +1,14 @@
-void runProof(const char * dataset = "/PWG0/COMMON/run30000X_10TeV_0.5T", Long64_t nentries=1234567890, Long64_t firstentry=0)
+void runProofTRDComparison(const char *dataset="/PWG0/COMMON/run30000X_10TeV_0.5T",Long64_t nentries=1000, Long64_t firstentry=0)
 {
   // Connect to Proof
   TProof::Open("lxb6046");
 
   // Upload and enable packages: please use the correct version!
-  gProof->UploadPackage("AF-v4-12");
-  gProof->EnablePackage("AF-v4-12");
+  gProof->UploadPackage("AF-v4-14");
+  gProof->EnablePackage("AF-v4-14");
 
   // Create the analysis manager
-  AliAnalysisManager *mgr = new AliAnalysisManager("AliTPCComparison");
+  AliAnalysisManager *mgr = new AliAnalysisManager("AliTRDComparison");
 
   AliVEventHandler* esdH = new AliESDInputHandler();
   mgr->SetInputEventHandler(esdH);
@@ -19,15 +19,18 @@ void runProof(const char * dataset = "/PWG0/COMMON/run30000X_10TeV_0.5T", Long64
 
   // Create task
   gProof->Load("AliMCComparisonTrack.cxx++g");
-  gProof->Load("AliTPCComparisonTask.cxx++g");
-  AliAnalysisTask *task = new AliTPCComparisonTask("AliTPCComparisonTask");
+  gProof->Load("AliTRDComparisonTask.cxx++g");
+  AliAnalysisTask *task = new AliTRDComparisonTask("AliTRDComparisonTask");
 
   // Add task
   mgr->AddTask(task);
 
   // Create containers for input/output
-  AliAnalysisDataContainer *cinput = mgr->CreateContainer("cchain", TChain::Class(), AliAnalysisManager::kInputContainer);
-  AliAnalysisDataContainer *coutput = mgr->CreateContainer("coutput", TList::Class(), AliAnalysisManager::kOutputContainer, "Hist.root");
+  AliAnalysisDataContainer* cinput = 
+    mgr->CreateContainer("cchain", TChain::Class(), AliAnalysisManager::kInputContainer);
+  AliAnalysisDataContainer* coutput = 
+    mgr->CreateContainer("coutput", TList::Class(), 
+    AliAnalysisManager::kOutputContainer, "AliTRDComparisonHist.root");
 
   // Connect input/output
   mgr->ConnectInput(task, 0, cinput);
@@ -44,3 +47,4 @@ void runProof(const char * dataset = "/PWG0/COMMON/run30000X_10TeV_0.5T", Long64
 
   mgr->StartAnalysis("proof",dataset,nentries,firstentry);
 }
+