Adding Analisys tasks for performance study
authormarian <marian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 14 Mar 2008 15:54:41 +0000 (15:54 +0000)
committermarian <marian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 14 Mar 2008 15:54:41 +0000 (15:54 +0000)
(Jacek Otwinowski)

12 files changed:
PWG1/AliComparisonDCA.cxx [new file with mode: 0644]
PWG1/AliComparisonDCA.h [new file with mode: 0644]
PWG1/AliComparisonDEdx.cxx [new file with mode: 0644]
PWG1/AliComparisonDEdx.h [new file with mode: 0644]
PWG1/AliComparisonEff.cxx [new file with mode: 0644]
PWG1/AliComparisonEff.h [new file with mode: 0644]
PWG1/AliComparisonRes.cxx [new file with mode: 0644]
PWG1/AliComparisonRes.h [new file with mode: 0644]
PWG1/AliComparisonTask.cxx [new file with mode: 0644]
PWG1/AliComparisonTask.h [new file with mode: 0644]
PWG1/PWG1LinkDef.h
PWG1/libPWG1.pkg

diff --git a/PWG1/AliComparisonDCA.cxx b/PWG1/AliComparisonDCA.cxx
new file mode 100644 (file)
index 0000000..d5f0ca8
--- /dev/null
@@ -0,0 +1,279 @@
+//------------------------------------------------------------------------------
+// Implementation of AliComparisonDCA class. It keeps information from 
+// comparison of reconstructed and MC particle tracks. In addtion, 
+// it keeps selection cuts used during comparison. The comparison 
+// information is stored in the ROOT histograms. Analysis of these 
+// histograms can be done by using Analyse() class function. The result of 
+// the analysis (histograms) are stored in the output picture_dca.root file.
+//  
+// Author: J.Otwinowski 04/02/2008 
+//------------------------------------------------------------------------------
+
+/*
+  // after running analysis, read the file, and get component
+  gSystem->Load("libPWG1.so");
+  TFile f("Output.root");
+  AliComparisonDCA * comp = (AliComparisonDCA*)f.Get("AliComparisonDCA");
+
+  // analyse comparison data (output stored in pictures_dca.root)
+  comp->Analyse();
+
+  // TPC track length parameterisation
+  TF1 fl("fl","((min(250./(abs(x+0.000001)),250)-90))",0,2);  // length function
+  TF1 fl2("fl2","[0]/((min(250./(abs(x+0.000001)),250)-90))^[1]",0,2);
+  fl2.SetParameter(1,1);
+  fl2.SetParameter(0,1);
+*/
+
+#include <iostream>
+
+#include "TFile.h"
+#include "TCint.h"
+#include "TH3F.h"
+#include "TH2F.h"
+#include "TF1.h"
+#include "TProfile.h"
+#include "TProfile2D.h"
+#include "TGraph2D.h"
+#include "TCanvas.h"
+#include "TGraph.h"
+// 
+#include "AliTracker.h"   
+#include "AliESDEvent.h"   
+#include "AliESD.h"
+#include "AliESDfriend.h"
+#include "AliESDfriendTrack.h"
+#include "AliRecInfoCuts.h" 
+#include "AliMCInfoCuts.h" 
+#include "AliLog.h" 
+#include "AliESDVertex.h" 
+//
+#include "AliMathBase.h"
+#include "AliTreeDraw.h" 
+
+#include "AliMCInfo.h" 
+#include "AliESDRecInfo.h" 
+#include "AliComparisonDCA.h" 
+
+using namespace std;
+
+ClassImp(AliComparisonDCA)
+
+//_____________________________________________________________________________
+AliComparisonDCA::AliComparisonDCA():
+  TNamed("AliComparisonDCA","AliComparisonDCA"),
+
+  // DCA histograms
+  fD0TanSPtB1(0),
+  fD1TanSPtB1(0),
+  fD0TanSPtL1(0),
+  fD1TanSPtL1(0),
+  fD0TanSPtInTPC(0),
+  fD1TanSPtInTPC(0),
+  fVertex(0),
+
+  // Cuts 
+  fCutsRC(0), 
+  fCutsMC(0)  
+{
+  InitHisto();
+  InitCuts();
+
+  // vertex (0,0,0)
+  fVertex = new AliESDVertex();
+  fVertex->SetXv(0.0);
+  fVertex->SetYv(0.0);
+  fVertex->SetZv(0.0);
+}
+
+//_____________________________________________________________________________
+AliComparisonDCA::~AliComparisonDCA()
+{
+  //
+  if(fD0TanSPtB1) delete fD0TanSPtB1; fD0TanSPtB1=0;
+  if(fD1TanSPtB1) delete fD1TanSPtB1; fD1TanSPtB1=0;
+  if(fD0TanSPtL1) delete fD0TanSPtL1; fD0TanSPtL1=0;
+  if(fD1TanSPtL1) delete fD1TanSPtL1; fD1TanSPtL1=0;
+  if(fD0TanSPtInTPC) delete fD0TanSPtInTPC; fD0TanSPtInTPC=0;
+  if(fD1TanSPtInTPC) delete fD1TanSPtInTPC; fD1TanSPtInTPC=0;
+  if(fVertex) delete fVertex; fVertex=0;
+}
+
+//_____________________________________________________________________________
+void AliComparisonDCA::InitHisto()
+{
+  // DCA histograms
+  fD0TanSPtB1 = new TH3F("DCAyTanSPtB1","DCAyTanSPt",20,1,2, 10,0.3,2, 100,-4,4);
+  fD0TanSPtB1->SetXTitle("tan(#theta)");
+  fD0TanSPtB1->SetYTitle("sqrt(p_{t})");
+  fD0TanSPtB1->SetZTitle("DCA_{xy}");
+
+  fD1TanSPtB1 = new TH3F("DCAzTanSPtB1","DCAzTanSPt",20,1,2, 10,0.3,2, 100,-4,4);
+  fD1TanSPtB1->SetXTitle("tan(#theta)");
+  fD1TanSPtB1->SetYTitle("sqrt(p_{t})");
+  fD1TanSPtB1->SetZTitle("DCA_{z}");
+
+  fD0TanSPtL1 = new TH3F("DCAyTanSPtL1","DCAyTanSPt",20,0,1, 10,0.3,2, 100,-0.1,0.1);
+  fD0TanSPtL1->SetXTitle("tan(#theta)");
+  fD0TanSPtL1->SetYTitle("sqrt(p_{t})");
+  fD0TanSPtL1->SetZTitle("DCA_{xy}");
+
+  fD1TanSPtL1 = new TH3F("DCAzTanSPtL1","DCAzTanSPt",20,0,1, 10,0.3,2, 100, -0.1,0.1);
+  fD1TanSPtL1->SetXTitle("tan(#theta)");
+  fD1TanSPtL1->SetYTitle("sqrt(p_{t})");
+  fD1TanSPtL1->SetZTitle("DCA_{z}");
+
+  fD0TanSPtInTPC = new TH3F("DCAyTanSPtInTPC","DCAyTanSPt",100,0,100, 10,0.3,2, 100,-0.1,0.1);
+  fD0TanSPtInTPC->SetXTitle("tan(#theta)");
+  fD0TanSPtInTPC->SetYTitle("sqrt(p_{t})");
+  fD0TanSPtInTPC->SetZTitle("DCA_{xy}");
+
+  fD1TanSPtInTPC = new TH3F("DCAzTanSPtInTPC","DCAzTanSPt",100,0,100, 10,0.3,2, 100, -0.1,0.1);
+  fD1TanSPtInTPC->SetXTitle("tan(#theta)");
+  fD1TanSPtInTPC->SetYTitle("sqrt(p_{t})");
+  fD1TanSPtInTPC->SetZTitle("DCA_{z}");
+}
+
+//_____________________________________________________________________________
+void AliComparisonDCA::InitCuts()
+{
+  if(!fCutsMC) 
+    AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
+  if(!fCutsRC) 
+    AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
+}
+
+//_____________________________________________________________________________
+void AliComparisonDCA::Process(AliMCInfo* infoMC, AliESDRecInfo *infoRC)
+{
+  // Fill DCA comparison information
+   
+  Float_t mcpt = infoMC->GetParticle().Pt();
+  Float_t tantheta = TMath::Tan(infoMC->GetParticle().Theta()-TMath::Pi()*0.5);
+  Bool_t isPrim = infoMC->GetParticle().R()<fCutsMC->GetMaxR() && TMath::Abs(infoMC->GetParticle().Vz())<fCutsMC->GetMaxVz() ;
+  Float_t spt = TMath::Sqrt(mcpt);
+
+  AliExternalTrackParam *track = 0;
+  Double_t kRadius    = 3.0;      // beam pipe radius
+  Double_t kMaxStep   = 5.0;      // max step
+  Double_t field      = 0.4;      // mag. field
+  Double_t kMaxD      = 123456.0; // max distance
+
+  Int_t clusterITS[200];
+  Double_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
+  const Double_t* dv;
+
+  // Check selection cuts
+  if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return; 
+  if (!isPrim) return;
+  if (infoRC->GetStatus(1)!=3) return;
+  if (!infoRC->GetESDtrack()) return;  
+  if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;
+  if (!infoRC->GetESDtrack()->GetConstrainedParam()) return;
+
+  // calculate and set prim. vertex
+  dv = infoMC->GetVDist(); // distance to Prim. vertex
+
+  fVertex->SetXv( infoMC->GetParticle().Vx() - dv[0] );
+  fVertex->SetYv( infoMC->GetParticle().Vy() - dv[1] );
+  fVertex->SetZv( infoMC->GetParticle().Vz() - dv[2] );
+
+  // calculate track parameters at vertex
+
+  if (infoRC->GetESDtrack()->GetTPCInnerParam())
+  {
+    if ((track = new AliExternalTrackParam(*infoRC->GetESDtrack()->GetTPCInnerParam())) != 0 )
+    {
+      AliTracker::PropagateTrackTo(track,kRadius,infoMC->GetMass(),kMaxStep,kTRUE);
+      track->PropagateToDCA(fVertex,field,kMaxD,dca,cov);
+
+      fD0TanSPtInTPC->Fill(tantheta,spt,dca[0]);
+      fD1TanSPtInTPC->Fill(tantheta,spt,dca[1]);
+
+         delete track;
+    }
+  }
+
+ if(infoRC->GetESDtrack()->GetITSclusters(clusterITS)==0){
+    fD0TanSPtB1->Fill(tantheta,spt,dca[0]);
+    fD1TanSPtB1->Fill(tantheta,spt,dca[1]);
+  }
+    fD0TanSPtL1->Fill(tantheta,spt,dca[0]);
+    fD1TanSPtL1->Fill(tantheta,spt,dca[1]);
+}
+
+//_____________________________________________________________________________
+Long64_t AliComparisonDCA::Merge(TCollection* list) 
+{
+  // Merge list of objects (needed by PROOF)
+
+  if (!list)
+  return 0;
+
+  if (list->IsEmpty())
+  return 1;
+
+  TIterator* iter = list->MakeIterator();
+  TObject* obj = 0;
+
+  // collection of generated histograms
+  Int_t count=0;
+  while((obj = iter->Next()) != 0) 
+  {
+    AliComparisonDCA* entry = dynamic_cast<AliComparisonDCA*>(obj);
+    if (entry == 0) { 
+      Error("Add","Attempt to add object of class: %s to a %s",
+         obj->ClassName(),this->ClassName());
+      return -1;
+    }
+
+    fD0TanSPtB1->Add(entry->fD0TanSPtB1);
+    fD1TanSPtB1->Add(entry->fD1TanSPtB1);
+    fD0TanSPtL1->Add(entry->fD0TanSPtL1);
+    fD1TanSPtL1->Add(entry->fD1TanSPtL1);
+    fD0TanSPtInTPC->Add(entry->fD0TanSPtInTPC);
+    fD1TanSPtInTPC->Add(entry->fD1TanSPtInTPC);
+
+    count++;
+  }
+
+return count;
+}
+
+//_____________________________________________________________________________
+void AliComparisonDCA::Exec(AliMCInfo* infoMC, AliESDRecInfo *infoRC){
+  // Process comparison information
+  Process(infoMC,infoRC);
+}
+
+//_____________________________________________________________________________
+void AliComparisonDCA::Analyse()
+{
+  // Analyse output histograms
+  
+  AliComparisonDCA * comp=this;
+  TGraph2D * gr=0;
+  TGraph * gr0=0;
+
+  TFile *fp = new TFile("pictures_dca.root","recreate");
+  fp->cd();
+
+  TCanvas * c = new TCanvas("DCA","DCA resloution");
+  c->cd();
+
+  // DCA resolution
+  gr0 = AliMathBase::MakeStat1D(comp->fD0TanSPtB1,2,5);
+  gr0->GetXaxis()->SetTitle("Tan(#theta)");
+  gr0->GetYaxis()->SetTitle("#sigmaDCA (cm)");
+  gr0->Write("DCAResolTan");
+  //
+  gr = AliMathBase::MakeStat2D(comp->fD0TanSPtB1,4,2,5); 
+  gr->GetXaxis()->SetTitle("Tan(#theta)");
+  gr->GetYaxis()->SetTitle("#sigmaDCA (cm)");
+  gr->GetHistogram()->Write("DCAResolSPTTan");
+
+  gPad->Clear();
+  gr0->Draw("al*");
+
+  fp->Close();
+}
diff --git a/PWG1/AliComparisonDCA.h b/PWG1/AliComparisonDCA.h
new file mode 100644 (file)
index 0000000..eaf8d9b
--- /dev/null
@@ -0,0 +1,72 @@
+#ifndef ALICOMPARISONDCA_H
+#define ALICOMPARISONDCA_H
+
+//------------------------------------------------------------------------------
+// Class to keep information from comparison of 
+// reconstructed and MC particle tracks (DCA - Distance of Closest Approach 
+// to the vertex).   
+// 
+// Author: J.Otwinowski 04/02/2008 
+//------------------------------------------------------------------------------
+
+class TFile;
+class AliMCInfo;
+class AliESDRecInfo;
+class AliESDEvent; 
+class AliESD;
+class AliESDfriend;
+class AliRecInfoCuts;
+class AliMCInfoCuts;
+class TH1I;
+class TH3F;
+class TH3;
+class TProfile;
+class TProfile2D;
+class AliESDVertex;
+
+#include "TNamed.h"
+
+class AliComparisonDCA : public TNamed {
+public :
+  AliComparisonDCA(); 
+  ~AliComparisonDCA();
+  void      InitHisto();
+  void      InitCuts();
+  void      Exec(AliMCInfo* infoMC, AliESDRecInfo *infoRC);
+  void      Process(AliMCInfo* infoMC, AliESDRecInfo *infoRC);
+  
+  // Selection cuts
+  void SetAliRecInfoCuts(AliRecInfoCuts* cuts=0) {fCutsRC = cuts;}
+  void SetAliMCInfoCuts(AliMCInfoCuts* cuts=0) {fCutsMC = cuts;}  
+
+  AliRecInfoCuts*  GetAliRecInfoCuts() const {return fCutsRC;}
+  AliMCInfoCuts*   GetAliMCInfoCuts()  const {return fCutsMC;}
+
+  // Merge output objects (needed by PROOF) 
+  virtual Long64_t Merge(TCollection* list);
+
+  // Analyse output histograms
+  void Analyse();
+
+private:
+  // DCA resolution
+  TH3F  *fD0TanSPtB1;     //-> distance to vertex y (no ITS clusters) 
+  TH3F  *fD1TanSPtB1;     //-> distance to vertex z (no ITS clusters) 
+  TH3F  *fD0TanSPtL1;     //-> distance to vertex y  
+  TH3F  *fD1TanSPtL1;     //-> distance to vertex z 
+  TH3F  *fD0TanSPtInTPC;  //-> distance to vertex y (Inner TPC track parameters) 
+  TH3F  *fD1TanSPtInTPC;  //-> distance to vertex z (Inner TPC track parameters)
+
+  AliESDVertex *fVertex;  //! 
+
+  // Global cuts objects
+  AliRecInfoCuts*  fCutsRC; // selection cuts for reconstructed tracks
+  AliMCInfoCuts*  fCutsMC;  // selection cuts for MC tracks
+
+  AliComparisonDCA(const AliComparisonDCA&); // not implemented
+  AliComparisonDCA& operator=(const AliComparisonDCA&); // not implemented
+
+  ClassDef(AliComparisonDCA,1);
+};
+
+#endif
diff --git a/PWG1/AliComparisonDEdx.cxx b/PWG1/AliComparisonDEdx.cxx
new file mode 100644 (file)
index 0000000..50e8d29
--- /dev/null
@@ -0,0 +1,285 @@
+//------------------------------------------------------------------------------
+// Implementation of AliComparisonDEdx class. It keeps information from 
+// comparison of reconstructed and MC particle tracks. In addtion, 
+// it keeps selection cuts used during comparison. The comparison 
+// information is stored in the ROOT histograms. Analysis of these 
+// histograms can be done by using Analyse() class function. The result of 
+// the analysis (histograms) are stored in the output picture_dedx.root file.
+//  
+// Author: J.Otwinowski 04/02/2008 
+//------------------------------------------------------------------------------
+
+
+/*
+  //after running analysis, read the file, and get component
+  gSystem->Load("libPWG1.so");
+  TFile f("Output.root");
+  AliComparisonDEdx * comp = (AliComparisonDEdx*)f.Get("AliComparisonDEdx");
+
+  // analyse comparison data (output stored in pictures_dedx.root)
+  comp->Analyse();
+
+  // TPC track length parameterisation
+  TF1 fl("fl","((min(250./(abs(x+0.000001)),250)-90))",0,2);  // length function
+  TF1 fl2("fl2","[0]/((min(250./(abs(x+0.000001)),250)-90))^[1]",0,2);
+  fl2.SetParameter(1,1);
+  fl2.SetParameter(0,1);
+*/
+
+#include <iostream>
+
+#include "TFile.h"
+#include "TCint.h"
+#include "TH3F.h"
+#include "TH2F.h"
+#include "TF1.h"
+#include "TProfile.h"
+#include "TProfile2D.h"
+#include "TGraph2D.h"
+#include "TCanvas.h"
+#include "TGraph.h"
+//
+#include "AliESDEvent.h"
+#include "AliESD.h"
+#include "AliESDfriend.h"
+#include "AliESDfriendTrack.h"
+#include "AliRecInfoCuts.h" 
+#include "AliMCInfoCuts.h" 
+#include "AliLog.h" 
+//
+#include "AliMathBase.h"
+#include "AliTreeDraw.h"
+//#include "TStatToolkit.h"
+
+#include "AliMCInfo.h" 
+#include "AliESDRecInfo.h" 
+#include "AliComparisonDEdx.h" 
+
+using namespace std;
+
+ClassImp(AliComparisonDEdx)
+
+//_____________________________________________________________________________
+AliComparisonDEdx::AliComparisonDEdx():
+  TNamed("AliComparisonDEdx","AliComparisonDEdx"),
+
+  // dEdx 
+  fTPCSignalNormTan(0), 
+  fTPCSignalNormSPhi(0),
+  fTPCSignalNormTPhi(0), 
+  //
+  fTPCSignalNormTanSPhi(0),
+  fTPCSignalNormTanTPhi(0),
+  fTPCSignalNormTanSPt(0), 
+  
+  // Cuts 
+  fCutsRC(0), 
+  fCutsMC(0),
+  fMCPtMin(0),
+  fMCAbsTanThetaMax(0),
+  fMCPdgCode(0)
+{
+  InitHisto();
+  InitCuts();
+}
+
+//_____________________________________________________________________________
+AliComparisonDEdx::~AliComparisonDEdx(){
+   
+  if(fTPCSignalNormTan)  delete fTPCSignalNormTan; fTPCSignalNormTan=0; 
+  if(fTPCSignalNormSPhi) delete fTPCSignalNormSPhi; fTPCSignalNormSPhi=0;
+  if(fTPCSignalNormTPhi) delete fTPCSignalNormTPhi; fTPCSignalNormTPhi=0;
+  //
+  if(fTPCSignalNormTanSPhi) delete fTPCSignalNormTanSPhi; fTPCSignalNormTanSPhi=0;
+  if(fTPCSignalNormTanTPhi) delete fTPCSignalNormTanTPhi; fTPCSignalNormTanTPhi=0;
+  if(fTPCSignalNormTanSPt)  delete fTPCSignalNormTanSPt; fTPCSignalNormTanSPt=0;
+}
+
+//_____________________________________________________________________________
+void AliComparisonDEdx::InitHisto()
+{
+  // Init histograms
+  
+  // TPC dEdx
+  fTPCSignalNormTan = new TH2F("CdEdxTan","CdEdxTan",50, -2,2,  40,30,70); 
+  fTPCSignalNormTan->SetXTitle("tan(#theta)");
+  fTPCSignalNormTan->SetYTitle("rec. dE/dx / calc. dE/dx");
+
+  fTPCSignalNormSPhi   = new TH2F("CdEdxSPhi","CdEdxSPhi",10,0.0,1,40,30,70);
+  fTPCSignalNormSPhi->SetXTitle("sin(#phi)");
+  fTPCSignalNormSPhi->SetYTitle("rec. dE/dx / calc. dE/dx");
+
+  fTPCSignalNormTPhi   = new TH2F("CdEdxTPhi","CdEdxTPhi",10,0.0,2,40,30,70); 
+  fTPCSignalNormTPhi->SetXTitle("tan(#phi)");
+  fTPCSignalNormTPhi->SetYTitle("rec. dE/dx / calc. dE/dx");
+
+  fTPCSignalNormTanSPhi= new TH3F("CdEdxTanSPhi","CdEdxTanSPhi",20, -2,2, 10,0.0 ,1,  40,30,70);
+  fTPCSignalNormTanSPhi->SetXTitle("tan(#theta)");
+  fTPCSignalNormTanSPhi->SetYTitle("sin(#phi)");
+  fTPCSignalNormTanSPhi->SetZTitle("rec. dE/dx / calc. dE/dx");
+
+  fTPCSignalNormTanTPhi= new TH3F("CdEdxTanTPhi","CdEdxTanTPhi",20, -2,2, 10,0.0 ,1,  40,30,70);
+  fTPCSignalNormTanTPhi->SetXTitle("tan(#theta)");
+  fTPCSignalNormTanTPhi->SetYTitle("tan(#phi)");
+  fTPCSignalNormTanTPhi->SetZTitle("rec. dE/dx / calc. dE/dx");
+
+  fTPCSignalNormTanSPt= new TH3F("CdEdxTanSPt","CdEdxTanSPt",20, -2,2, 10,0.3 ,3, 40,30,70); 
+  fTPCSignalNormTanSPt->SetXTitle("tan(#theta)");
+  fTPCSignalNormTanSPt->SetYTitle("#sqrt(p_{t})");
+  fTPCSignalNormTanSPt->SetZTitle("rec. dE/dx / calc. dE/dx");
+}
+
+void AliComparisonDEdx::InitCuts()
+{
+  // Init cuts
+  if(!fCutsMC) 
+    AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
+  if(!fCutsRC) 
+    AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
+}
+
+//_____________________________________________________________________________
+void AliComparisonDEdx::Process(AliMCInfo* infoMC, AliESDRecInfo *infoRC){
+
+  // Fill dE/dx  comparison information
+  
+  Float_t mcpt = infoMC->GetParticle().Pt();
+  Float_t tantheta = TMath::Tan(infoMC->GetParticle().Theta()-TMath::Pi()*0.5);
+  Bool_t isPrim = infoMC->GetParticle().R()<fCutsMC->GetMaxR() && TMath::Abs(infoMC->GetParticle().Vz())<fCutsMC->GetMaxVz() ;
+
+  Float_t mprim = infoMC->GetPrim();
+  
+  // Check selection cuts 
+  if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return; 
+  if (!isPrim) return;
+  if (infoRC->GetStatus(1)!=3) return;
+  if (!infoRC->GetESDtrack()) return;  
+  if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;
+  if (!infoRC->GetESDtrack()->GetConstrainedParam()) return;
+  //if (mprim>1.4) return;
+  //if (mprim<0.5) return;
+  if (mprim > fCutsMC->GetMaxTPCSignal()) return;
+  if (mprim < fCutsMC->GetMinTPCSignal()) return;
+  if (infoRC->GetESDtrack()->GetTPCsignalN()<fCutsRC->GetMinTPCsignalN()) return;
+  //
+  Float_t ratio = infoRC->GetESDtrack()->GetTPCsignal()/infoMC->GetPrim();
+  Float_t sphi =  infoRC->GetESDtrack()->GetInnerParam()->GetSnp();
+  Float_t tphi =  sphi/TMath::Sqrt(1-sphi*sphi);
+
+  if (TMath::Abs(infoMC->GetParticle().GetPdgCode()) != GetMCPdgCode()) return;
+  //if (mcpt>0.5) {
+  if (mcpt > GetMCPtMin()) {
+    fTPCSignalNormTan->Fill(tantheta,ratio);    // only subset
+  }
+
+  //if (TMath::Abs(tantheta)<0.5){
+  if (TMath::Abs(tantheta) < GetMCAbsTanThetaMax()){
+    fTPCSignalNormSPhi->Fill(sphi,ratio);       // only subset
+    fTPCSignalNormTPhi->Fill(tphi,ratio);       // only subset
+  }
+  fTPCSignalNormTanSPhi->Fill(tantheta,sphi,ratio);    
+  fTPCSignalNormTanTPhi->Fill(tantheta,tphi,ratio);    
+  fTPCSignalNormTanSPt->Fill(tantheta,TMath::Sqrt(mcpt),ratio);
+}
+
+//_____________________________________________________________________________
+Long64_t AliComparisonDEdx::Merge(TCollection* list) 
+{
+  // Merge list of objects (needed by PROOF)
+
+  if (!list)
+  return 0;
+
+  if (list->IsEmpty())
+  return 1;
+
+  TIterator* iter = list->MakeIterator();
+  TObject* obj = 0;
+
+  // collection of generated histograms
+  Int_t count=0;
+  while((obj = iter->Next()) != 0) 
+  {
+    AliComparisonDEdx* entry = dynamic_cast<AliComparisonDEdx*>(obj);
+    if (entry == 0) { 
+      Error("Add","Attempt to add object of class: %s to a %s",
+         obj->ClassName(),this->ClassName());
+      return -1;
+  }
+
+  fTPCSignalNormTan->Add(entry->fTPCSignalNormTan);
+  fTPCSignalNormSPhi->Add(entry->fTPCSignalNormSPhi);
+  fTPCSignalNormTPhi->Add(entry->fTPCSignalNormTPhi);
+  //
+  fTPCSignalNormTanSPhi->Add(entry->fTPCSignalNormTanSPhi);
+  fTPCSignalNormTanTPhi->Add(entry->fTPCSignalNormTanTPhi);
+  fTPCSignalNormTanSPt->Add(entry->fTPCSignalNormTanSPt);
+
+  count++;
+  }
+
+return count;
+}
+
+//_____________________________________________________________________________
+void AliComparisonDEdx::Exec(AliMCInfo* infoMC, AliESDRecInfo *infoRC)
+{
+  // Process comparison information
+  Process(infoMC,infoRC);
+}
+
+//_____________________________________________________________________________
+TH1F* AliComparisonDEdx::MakeResol(TH2F * his, Int_t integ, Bool_t type)
+{
+  // Make resolution histograms
+  TH1F *hisr, *hism;
+  if (!gPad) new TCanvas;
+  hisr = AliTreeDraw::CreateResHistoI(his,&hism,integ);
+  if (type) return hism;
+  else 
+    return hisr;
+}
+
+//_____________________________________________________________________________
+void AliComparisonDEdx::Analyse()
+{
+  // Analyse output histograms
+  
+  AliComparisonDEdx * comp=this;
+  TH1F *hiss=0;
+  TGraph2D * gr=0;
+
+  TFile *fp = new TFile("pictures_dedx.root","recreate");
+  fp->cd();
+
+  TCanvas * c = new TCanvas("TPCdedx","TPC dedx");
+  c->cd();
+
+  hiss = comp->MakeResol(comp->fTPCSignalNormTan,4,0);
+  hiss->SetXTitle("Tan(#theta)");
+  hiss->SetYTitle("#sigma_{dEdx}");
+  hiss->Draw();
+  hiss->Write("TPCdEdxResolTan");
+  //
+  hiss = comp->MakeResol(comp->fTPCSignalNormTan,4,1); 
+  hiss->SetXTitle("Tan(#theta)");
+  hiss->SetYTitle("<dEdx>");
+  hiss->Draw(); 
+  hiss->Write("TPCdEdxMeanTan");
+  //
+  gr = AliMathBase::MakeStat2D(comp->fTPCSignalNormTanSPt,3,1,4);
+  gr->GetXaxis()->SetTitle("Tan(#theta)");
+  gr->GetYaxis()->SetTitle("#sqrt{p_{t}(GeV)}");
+  gr->GetZaxis()->SetTitle("<dEdx>");
+  gr->Draw("colz"); 
+  gr->GetHistogram()->Write("TPCdEdxMeanTanPt_1");
+  //
+  gr = AliMathBase::MakeStat2D(comp->fTPCSignalNormTanSPt,3,1,5);
+  gr->GetXaxis()->SetTitle("Tan(#theta)");
+  gr->GetYaxis()->SetTitle("#sqrt{p_{t}(GeV)}");
+  gr->GetZaxis()->SetTitle("#sigma_{dEdx}");
+  gr->Draw("colz"); 
+  gr->GetHistogram()->Write("TPCdEdxMeanTanPt_2");
+
+  fp->Close();
+}
diff --git a/PWG1/AliComparisonDEdx.h b/PWG1/AliComparisonDEdx.h
new file mode 100644 (file)
index 0000000..6b6a345
--- /dev/null
@@ -0,0 +1,85 @@
+#ifndef ALICOMPARISONDEdx_H
+#define ALICOMPARISONDEdx_H
+
+//------------------------------------------------------------------------------
+// Class to keep information from comparison of 
+// reconstructed and MC particle tracks (TPC dE/dx).   
+// 
+// Author: J.Otwinowski 04/02/2008 
+//------------------------------------------------------------------------------
+
+class TFile;
+class AliMCInfo;
+class AliESDRecInfo;
+class AliESDEvent; 
+class AliESD;
+class AliESDfriend;
+class AliRecInfoCuts;
+class AliMCInfoCuts;
+class TH1I;
+class TH3F;
+class TH3;
+class TProfile;
+class TProfile2D;
+class TGraph2D;
+class TGraph; 
+
+#include "TNamed.h"
+
+class AliComparisonDEdx : public TNamed {
+public :
+  AliComparisonDEdx(); 
+  ~AliComparisonDEdx();
+  void      InitHisto();
+  void      InitCuts();
+  void      Exec(AliMCInfo* infoMC, AliESDRecInfo *infoRC);
+  void      Process(AliMCInfo* infoMC, AliESDRecInfo *infoRC);
+
+  // Selection cuts
+  void SetAliRecInfoCuts(AliRecInfoCuts* cuts=0) {fCutsRC = cuts;}
+  void SetAliMCInfoCuts(AliMCInfoCuts* cuts=0)   {fCutsMC = cuts;} 
+
+  void SetMCPtMin(const Float_t cuts=0) {fMCPtMin = cuts;} 
+  void SetMCAbsTanThetaMax(const Float_t cuts=1e99) {fMCAbsTanThetaMax = cuts;} 
+  void SetMCPdgCode(const Int_t cuts=0) {fMCPdgCode = cuts;} 
+
+  AliRecInfoCuts*  GetAliRecInfoCuts() const {return fCutsRC;}      
+  AliMCInfoCuts*   GetAliMCInfoCuts()  const {return fCutsMC;}     
+  Float_t GetMCPtMin() const {return fMCPtMin;}
+  Float_t GetMCAbsTanThetaMax() const {return fMCAbsTanThetaMax;}
+  Int_t GetMCPdgCode() const {return fMCPdgCode;} 
+
+  static TH1F*     MakeResol(TH2F * his, Int_t integ, Bool_t type); 
+
+  // Merge output objects (needed by PROOF) 
+  virtual Long64_t Merge(TCollection* list);
+
+  // Analyse output histograms
+  void Analyse();
+
+private:
+
+  // TPC dE/dx 
+  TH2F* fTPCSignalNormTan;    //-> TPC signal normalized to the calculated MC signal 
+  TH2F* fTPCSignalNormSPhi;   //-> TPC signal normalized to the calculated MC signal
+  TH2F* fTPCSignalNormTPhi;   //-> TPC signal normalized to the calculated MC signal
+  //
+  TH3F* fTPCSignalNormTanSPhi;   //-> TPC signal normalized to the calculated MC signal
+  TH3F* fTPCSignalNormTanTPhi;   //-> TPC signal normalized to the calculated MC signal
+  TH3F* fTPCSignalNormTanSPt;    //-> TPC signal normalized to the calculated MC signal
+  
+  // Selection cuts
+  AliRecInfoCuts*  fCutsRC; // selection cuts for reconstructed tracks
+  AliMCInfoCuts*   fCutsMC; // selection cuts for MC tracks
+
+  Float_t fMCPtMin;               // min. MC pt cut
+  Float_t fMCAbsTanThetaMax;      // max. MC abs[tan(theta)] cut
+  Int_t fMCPdgCode;               // selected particle with Pdg code
+
+  AliComparisonDEdx(const AliComparisonDEdx&); // not implemented
+  AliComparisonDEdx& operator=(const AliComparisonDEdx&); // not implemented
+
+  ClassDef(AliComparisonDEdx,1);
+};
+
+#endif
diff --git a/PWG1/AliComparisonEff.cxx b/PWG1/AliComparisonEff.cxx
new file mode 100644 (file)
index 0000000..70e1f30
--- /dev/null
@@ -0,0 +1,841 @@
+//------------------------------------------------------------------------------
+// Implementation of AliComparisonEff class. It keeps information from 
+// comparison of reconstructed and MC particle tracks. In addtion, 
+// it keeps selection cuts used during comparison. The comparison 
+// information is stored in the ROOT histograms. Analysis of these 
+// histograms can be done by using Analyse() class function. The result of 
+// the analysis (histograms) are stored in the output picture_eff.root file.
+//  
+// Author: J.Otwinowski 04/02/2008 
+//------------------------------------------------------------------------------
+
+/*
+  // after running analysis, read the file, and get component
+  gSystem->Load("libPWG1.so");
+  TFile f("Output.root");
+  AliComparisonEff * comp = (AliComparisonEff*)f.Get("AliComparisonEff");
+
+
+  // analyse comparison data (output stored in pictures_eff.root)
+  comp->Analyse();
+  // paramtetrisation of the TPC track length (for information only)
+  TF1 fl("fl","((min(250./(abs(x+0.000001)),250)-90))",0,2);  // length function
+  TF1 fl2("fl2","[0]/((min(250./(abs(x+0.000001)),250)-90))^[1]",0,2);
+  fl2.SetParameter(1,1);
+  fl2.SetParameter(0,1);
+*/
+
+
+#include <iostream>
+
+#include "TFile.h"
+#include "TCint.h"
+#include "TH3F.h"
+#include "TH2F.h"
+#include "TF1.h"
+#include "TProfile.h"
+#include "TProfile2D.h"
+#include "TGraph2D.h"
+#include "TCanvas.h"
+#include "TGraph.h"
+// 
+#include "AliESDEvent.h"  
+#include "AliESD.h"
+#include "AliESDfriend.h"
+#include "AliESDfriendTrack.h"
+#include "AliRecInfoCuts.h" 
+#include "AliMCInfoCuts.h" 
+#include "AliLog.h" 
+//
+#include "AliMathBase.h"
+#include "AliTreeDraw.h" 
+#include "AliMagFMaps.h" 
+#include "AliESDVertex.h" 
+#include "AliExternalTrackParam.h" 
+#include "AliTracker.h" 
+
+#include "AliMCInfo.h" 
+#include "AliESDRecInfo.h" 
+#include "AliComparisonEff.h" 
+
+using namespace std;
+
+
+ClassImp(AliComparisonEff)
+
+//_____________________________________________________________________________
+AliComparisonEff::AliComparisonEff():
+  TNamed("AliComparisonEff","AliComparisonEff"),
+
+  // histograms
+  fMCPt(0),
+  fMCRecPt(0),
+  fMCRecPrimPt(0),
+  fMCRecSecPt(0),
+
+  fEffTPCPt(0),
+  fEffTPCPtMC(0),
+  fEffTPCPtF(0),
+  //
+  fEffTPCPt_P(0),
+  fEffTPCPtMC_P(0),
+  fEffTPCPtF_P(0),
+  //
+  fEffTPCPt_Pi(0),
+  fEffTPCPtMC_Pi(0),
+  fEffTPCPtF_Pi(0),
+  //
+  fEffTPCPt_K(0),
+  fEffTPCPtMC_K(0),
+  fEffTPCPtF_K(0),
+  fEffTPCTan(0),
+  fEffTPCTanMC(0),
+  fEffTPCTanF(0),
+  //
+  fEffTPCPtTan(0),
+  fEffTPCPtTanMC(0),
+  fEffTPCPtTanF(0),
+
+  // Cuts 
+  fCutsRC(0), 
+  fCutsMC(0),
+
+  fVertex(0)
+{
+  // init vertex
+  fVertex = new AliESDVertex();
+  fVertex->SetXv(0.0); fVertex->SetYv(0.0); fVertex->SetZv(0.0); 
+
+  for(Int_t i=0; i<4; ++i)
+  {
+    fTPCPtDCASigmaIdeal[i]=0;
+    fTPCPtDCASigmaFull[i]=0;
+    fTPCPtDCASigmaDay0[i]=0;
+
+    fTPCPtDCAXY[i]=0;
+    fTPCPtDCAZ[i]=0;
+
+       fTPCPtDCASigmaIdealPid[i]=0;
+       fTPCPtDCASigmaFullPid[i]=0;
+       fTPCPtDCASigmaDay0Pid[i]=0;
+
+       fTPCPtDCAXYPid[i]=0;   
+       fTPCPtDCAZPid[i]=0; 
+  }
+  InitHisto();
+  InitCuts();
+}
+
+//_____________________________________________________________________________
+AliComparisonEff::~AliComparisonEff(){
+
+  // 
+  if(fMCPt)  delete  fEffTPCPt; fEffTPCPt=0;
+  if(fMCRecPt)  delete  fMCRecPt; fMCRecPt=0;
+  if(fMCRecPrimPt)  delete  fMCRecPrimPt; fMCRecPrimPt=0;
+  if(fMCRecSecPt)  delete  fMCRecSecPt; fMCRecSecPt=0;
+
+  // 
+  if(fEffTPCPt)   delete  fEffTPCPt;   fEffTPCPt=0;
+  if(fEffTPCPtMC) delete  fEffTPCPtMC; fEffTPCPtMC=0;
+  if(fEffTPCPtF)  delete  fEffTPCPtF;  fEffTPCPtF=0;
+
+  // 
+  if(fEffTPCPt_P)   delete  fEffTPCPt_P;   fEffTPCPt_P=0;
+  if(fEffTPCPtMC_P) delete  fEffTPCPtMC_P; fEffTPCPtMC_P=0;
+  if(fEffTPCPtF_P)  delete  fEffTPCPtF_P;  fEffTPCPtF_P=0;
+
+  // 
+  if(fEffTPCPt_Pi)   delete  fEffTPCPt_Pi;   fEffTPCPt_Pi=0;
+  if(fEffTPCPtMC_Pi) delete  fEffTPCPtMC_Pi; fEffTPCPtMC_Pi=0;
+  if(fEffTPCPtF_Pi)  delete  fEffTPCPtF_Pi;  fEffTPCPtF_Pi=0;
+
+  // 
+  if(fEffTPCPt_K)   delete  fEffTPCPt_K;   fEffTPCPt_K=0;
+  if(fEffTPCPtMC_K) delete  fEffTPCPtMC_K; fEffTPCPtMC_K=0;
+  if(fEffTPCPtF_K)  delete  fEffTPCPtF_K;  fEffTPCPtF_K=0;
+
+  // 
+  if(fEffTPCTan)   delete  fEffTPCTan;   fEffTPCTan=0;
+  if(fEffTPCTanMC) delete  fEffTPCTanMC; fEffTPCTanMC=0;
+  if(fEffTPCTanF)  delete  fEffTPCTanF;  fEffTPCTanF=0;
+
+  //
+  if(fEffTPCPtTan)   delete  fEffTPCPtTan;   fEffTPCPtTan=0;
+  if(fEffTPCPtTanMC) delete  fEffTPCPtTanMC; fEffTPCPtTanMC=0;
+  if(fEffTPCPtTanF)  delete  fEffTPCPtTanF;  fEffTPCPtTanF=0;
+
+  for(Int_t i=0; i<4; ++i)
+  {
+    if(fTPCPtDCASigmaIdeal[i]) delete  fTPCPtDCASigmaIdeal[i];  fTPCPtDCASigmaIdeal[i]=0;
+    if(fTPCPtDCASigmaFull[i]) delete  fTPCPtDCASigmaFull[i];  fTPCPtDCASigmaFull[i]=0;
+    if(fTPCPtDCASigmaDay0[i]) delete  fTPCPtDCASigmaDay0[i];  fTPCPtDCASigmaDay0[i]=0;
+
+    if(fTPCPtDCAXY[i]) delete  fTPCPtDCAXY[i];  fTPCPtDCAXY[i]=0;
+    if(fTPCPtDCAZ[i]) delete  fTPCPtDCAZ[i];  fTPCPtDCAZ[i]=0;
+
+       if(fTPCPtDCASigmaIdealPid[i]) delete  fTPCPtDCASigmaIdealPid[i];  fTPCPtDCASigmaIdealPid[i]=0;
+       if(fTPCPtDCASigmaFullPid[i]) delete  fTPCPtDCASigmaFullPid[i];  fTPCPtDCASigmaFullPid[i]=0;
+       if(fTPCPtDCASigmaDay0Pid[i]) delete  fTPCPtDCASigmaDay0Pid[i];  fTPCPtDCASigmaDay0Pid[i]=0;
+
+       if(fTPCPtDCAXYPid[i]) delete  fTPCPtDCAXYPid[i]; fTPCPtDCAXYPid[i]=0;
+       if(fTPCPtDCAZPid[i]) delete   fTPCPtDCAZPid[i];  fTPCPtDCAZPid[i]=0;
+  }
+}
+
+//_____________________________________________________________________________
+void AliComparisonEff::InitHisto(){
+
+  // Init histograms
+  //
+  fMCPt = new TH1F("fMCPt","fMCPt",50,0.1,3);           
+  fMCPt->SetXTitle("p_{t}");
+  fMCPt->SetYTitle("yield");
+
+  fMCRecPt = new TH1F("fMCRecPt","fMCRecPt",50,0.1,3);           
+  fMCRecPt->SetXTitle("p_{t}");
+  fMCRecPt->SetYTitle("yield");
+
+  fMCRecPrimPt = new TH1F("fMCRecPrimPt","fMCRecPrimPt",50,0.1,3);           
+  fMCRecPrimPt->SetXTitle("p_{t}");
+  fMCRecPrimPt->SetYTitle("yield");
+
+  fMCRecSecPt = new TH1F("fMCRecSecPt","fMCRecSecPt",50,0.1,3);           
+  fMCRecSecPt->SetXTitle("p_{t}");
+  fMCRecSecPt->SetYTitle("yield");
+
+  // Efficiency as function of pt
+  fEffTPCPt = new TProfile("Eff_pt","Eff_Pt",50,0.1,3);           
+  fEffTPCPt->SetXTitle("p_{t}");
+  fEffTPCPt->SetYTitle("TPC Efficiency");
+
+  fEffTPCPtMC = new TProfile("MC_Eff_pt","MC_Eff_Pt",50,0.1,3);   
+  fEffTPCPtMC->SetXTitle("p_{t}");
+  fEffTPCPtMC->SetYTitle("MC TPC Efficiency");
+
+  fEffTPCPtF = new TProfile("F_Eff_pt","F_Eff_Pt",50,0.1,3);     
+  fEffTPCPtF->SetXTitle("p_{t}");
+  fEffTPCPtF->SetYTitle("TPC Findable Efficiency");
+
+  // Efficiency as function of pt protons
+  fEffTPCPt_P = new TProfile("Eff_pt_P","Eff_Pt_P",50,0.1,3);           
+  fEffTPCPt_P->SetXTitle("p_{t}");
+  fEffTPCPt_P->SetYTitle("TPC Efficiency");
+
+  fEffTPCPtMC_P = new TProfile("MC_Eff_pt_P","MC_Eff_Pt_P",50,0.1,3);   
+  fEffTPCPtMC_P->SetXTitle("p_{t}");
+  fEffTPCPtMC_P->SetYTitle("MC TPC Efficiency");
+
+  fEffTPCPtF_P = new TProfile("F_Eff_pt_P","F_Eff_Pt_P",50,0.1,3);     
+  fEffTPCPtF_P->SetXTitle("p_{t}");
+  fEffTPCPtF_P->SetYTitle("TPC Findable Efficiency");
+
+  // Efficiency as function of pt pions
+  fEffTPCPt_Pi = new TProfile("Eff_pt_Pi","Eff_Pit_Pi",50,0.1,3);           
+  fEffTPCPt_Pi->SetXTitle("p_{t}");
+  fEffTPCPt_Pi->SetYTitle("TPC Efficiency");
+
+  fEffTPCPtMC_Pi = new TProfile("MC_Eff_pt_Pi","MC_Eff_Pit_Pi",50,0.1,3);   
+  fEffTPCPtMC_Pi->SetXTitle("p_{t}");
+  fEffTPCPtMC_Pi->SetYTitle("MC TPC Efficiency");
+
+  fEffTPCPtF_Pi = new TProfile("F_Eff_pt_Pi","F_Eff_Pit_Pi",50,0.1,3);     
+  fEffTPCPtF_Pi->SetXTitle("p_{t}");
+  fEffTPCPtF_Pi->SetYTitle("TPC Findable Efficiency");
+
+  // Efficiency as function of pt kaons
+  fEffTPCPt_K = new TProfile("Eff_pt_K","Eff_Kt_K",50,0.1,3);           
+  fEffTPCPt_K->SetXTitle("p_{t}");
+  fEffTPCPt_K->SetYTitle("TPC Efficiency");
+
+  fEffTPCPtMC_K = new TProfile("MC_Eff_pt_K","MC_Eff_Kt_K",50,0.1,3);   
+  fEffTPCPtMC_K->SetXTitle("p_{t}");
+  fEffTPCPtMC_K->SetYTitle("MC TPC Efficiency");
+
+  fEffTPCPtF_K = new TProfile("F_Eff_pt_K","F_Eff_Kt_K",50,0.1,3);     
+  fEffTPCPtF_K->SetXTitle("p_{t}");
+  fEffTPCPtF_K->SetYTitle("TPC Findable Efficiency");
+
+  // Efficiency as function of tan(theta) 
+  fEffTPCTan = new TProfile("Eff_tan","Eff_tan",50,-2.5,2.5);         
+  fEffTPCTan->SetXTitle("tan(#theta)");
+  fEffTPCTan->SetYTitle("TPC Efficiency");
+
+  fEffTPCTanMC = new TProfile("MC_Eff_tan","MC_Eff_tan",50,-2.5,2.5); 
+  fEffTPCTanMC->SetXTitle("tan(#theta)");
+  fEffTPCTanMC->SetYTitle("MC TPC Efficiency");
+
+  fEffTPCTanF = new TProfile("F_Eff_tan","F_Eff_tan",50,-2.5,2.5);   
+  fEffTPCPtF->SetXTitle("tan(#theta)");
+  fEffTPCPtF->SetYTitle("TPC Findable Efficiency");
+
+  // Efficiency as function of pt and tan(theta) 
+  fEffTPCPtTan = new TProfile2D("Eff_pt_tan","Eff_Pt_tan",10,0.1,3,20,-2.,2.);
+  fEffTPCPtTan->SetXTitle("tan(#theta)");
+  fEffTPCPtTan->SetYTitle("p_{t}");
+
+  fEffTPCPtTanMC = new TProfile2D("MC_Eff_pt_tan_MC","MC Eff Pt",10,0.1,3,20,-2.,2.);
+  fEffTPCPtTanMC->SetXTitle("tan(#theta)");
+  fEffTPCPtTanMC->SetYTitle("p_{t}");
+
+  fEffTPCPtTanF = new TProfile2D("MC_Eff_pt_tan_F","MC Eff Pt",10,0.1,3,20,-2.,2.);
+  fEffTPCPtTanF->SetXTitle("tan(#theta)");
+  fEffTPCPtTanF->SetYTitle("p_{t}");
+
+  char name[256];
+  for(Int_t i=0; i<4; ++i)
+  {
+    sprintf(name, "fTPCPtDCASigmaIdeal_%d",i);
+    fTPCPtDCASigmaIdeal[i] = new TH2F(name,name,50,0.1,3,100,0,100);
+
+    sprintf(name, "fTPCPtDCASigmaFull_%d",i);
+    fTPCPtDCASigmaFull[i] = new TH2F(name,name,50,0.1,3,100,0,100);
+
+    sprintf(name, "fTPCPtDCASigmaDay0_%d",i);
+    fTPCPtDCASigmaDay0[i] = new TH2F(name,name,50,0.1,3,100,0,100);
+
+    sprintf(name, "fTPCPtDCAXY_%d",i);
+    fTPCPtDCAXY[i]= new TH2F(name,name,50,0.1,3,100,0,100);
+    sprintf(name, "fTPCPtDCAZ_%d",i);
+    fTPCPtDCAZ[i]= new TH2F(name,name,50,0.1,3,100,0,100);
+
+    sprintf(name, "fTPCPtDCASigmaIdealPid_%d",i);
+       fTPCPtDCASigmaIdealPid[i] = new TH3F(name,name,50,0.1,3,100,0,100,5,0,5);
+
+    sprintf(name, "fTPCPtDCASigmaFullPid_%d",i);
+       fTPCPtDCASigmaFullPid[i]= new TH3F(name,name,50,0.1,3,100,0,100,5,0,5);
+
+    sprintf(name, "fTPCPtDCASigmaDay0Pid_%d",i);
+       fTPCPtDCASigmaDay0Pid[i]= new TH3F(name,name,50,0.1,3,100,0,100,5,0,5);
+
+    sprintf(name, "fTPCPtDCAXYPid_%d",i);
+       fTPCPtDCAXYPid[i]= new TH3F(name,name,50,0.1,3,100,0,100,5,0,5);
+
+    sprintf(name, "fTPCPtDCAZPid_%d",i);
+       fTPCPtDCAZPid[i]= new TH3F(name,name,50,0.1,3,100,0,100,5,0,5);
+  }
+}
+
+//_____________________________________________________________________________
+void AliComparisonEff::InitCuts()
+{
+
+  if(!fCutsMC) 
+    AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
+  if(!fCutsRC) 
+    AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
+}
+
+//_____________________________________________________________________________
+void AliComparisonEff::Process(AliMCInfo* infoMC, AliESDRecInfo *infoRC)
+{
+  // Fill efficiency comparison information
+  
+  AliExternalTrackParam *track = 0;
+  Double_t kRadius    = 3.0;      // beam pipe radius
+  Double_t kMaxStep   = 5.0;      // max step
+  Double_t field      = 0.4;      // mag. field 
+  Double_t kMaxD      = 123456.0; // max distance
+
+  // systematics
+  const Double_t kSigma2Full_xy  = 0.25; // ExB full systematics  [cm]
+  const Double_t kSigma2Full_z  =  5.0;  // [cm] 
+
+  const Double_t kSigma2Day0_xy  = 0.02; //  ExB  [cm]
+  const Double_t kSigma2Day0_z  =  0.1; //  [cm] goofie  
+
+  //  
+  Double_t     DCASigmaIdeal=0;
+  Double_t  DCASigmaFull=0;
+  Double_t  DCASigmaDay0=0;
+
+  Double_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
+  const Double_t* dv;
+
+  // distance to Prim. vertex 
+  dv = infoMC->GetVDist(); 
+
+  Float_t mcpt = infoMC->GetParticle().Pt();
+  Float_t tantheta = TMath::Tan(infoMC->GetParticle().Theta()-TMath::Pi()*0.5);
+  //Bool_t isPrim = infoMC->GetParticle().R()<fCutsMC->GetMaxR() && TMath::Abs(infoMC->GetParticle().Vz())<fCutsMC->GetMaxVz();
+  Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])<fCutsMC->GetMaxR() && TMath::Abs(dv[2])<fCutsMC->GetMaxVz();
+  // calculate and set prim. vertex
+  fVertex->SetXv( infoMC->GetParticle().Vx() - dv[0] );
+  fVertex->SetYv( infoMC->GetParticle().Vy() - dv[1] );
+  fVertex->SetZv( infoMC->GetParticle().Vz() - dv[2] );
+  
+  // Check selection cuts
+  if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return; 
+
+  // transform Pdg to Pid
+  Double_t pid = -1;
+  if( TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetEP() ) pid = 0; 
+  if( TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetMuP() ) pid = 1; 
+  if( TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetPiP() ) pid = 2; 
+  if( TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetKP() ) pid = 3; 
+  if( TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetProt() ) pid = 4; 
+
+  //cout << "dv[0] " << dv[0] << " dv[1] " << dv[1]  <<  " dv[2] " << dv[2] << endl; 
+  //cout << "v[0] " << fVertex->GetXv()  << " v[1] " << fVertex->GetYv()  <<  " v[2] " << fVertex->GetZv()<< endl; 
+  if (TMath::Abs(tantheta)<fCutsRC->GetMaxAbsTanTheta())
+  {
+  
+    // MC pt
+    fMCPt->Fill(mcpt); 
+
+
+    if (infoRC->GetESDtrack()->GetTPCInnerParam()) 
+    {
+      if ((track = new AliExternalTrackParam(*infoRC->GetESDtrack()->GetTPCInnerParam())) != 0 )
+      {
+        AliTracker::PropagateTrackTo(track,kRadius,infoMC->GetMass(),kMaxStep,kTRUE);
+        track->PropagateToDCA(fVertex,field,kMaxD,dca,cov);
+
+               //
+               cov[2] = track->GetCovariance()[2];
+
+           // Eff = infoRC->GetStatus(1)==3 && isPrim / isPrim;
+        // Pt vs ( dca[0]^2/cov[0]^2 + dca[1]^2/cov[2]^2 ) 
+        // Pt vs ( dca[0]^2/(cov[0]^2 + kSigma2Full_xy)  + dca[1]^2/( cov[2]^2 + kSigma2Full_z ) 
+        // Pt vs ( dca[0]^2/(cov[0]^2 + kSigma2_xy)  + dca[1]^2/( cov[2]^2 + kSigma2_z ) 
+
+         if(cov[0]>0.0 && cov[2]>0.0)
+                {
+                        DCASigmaIdeal = TMath::Power(dca[0],2)/cov[0] 
+                                                                         + TMath::Power(dca[1],2)/cov[2]; 
+
+                        DCASigmaFull = TMath::Power(dca[0],2)/(cov[0]+kSigma2Full_xy) 
+                                                                        + TMath::Power(dca[1],2)/(cov[2]+kSigma2Full_z); 
+
+                        DCASigmaDay0 = TMath::Power(dca[0],2)/(cov[0]+kSigma2Day0_xy) 
+                                                                        + TMath::Power(dca[1],2)/(cov[2]+kSigma2Day0_z); 
+
+              //cout << "dca[0] " << dca[0]  << " dca[1] " << dca[1]  << endl; 
+              //cout << "cov[0] " << cov[0]  << " cov[2] " << cov[2]  << endl; 
+              //cout << "DCASigmaIdeal " << DCASigmaIdeal  << " DCASigmaFull " << DCASigmaFull  << " DCASigmaDay0 "  <<DCASigmaDay0 <<  endl; 
+            //cout << " -------------------------------------------------------- "<<  endl; 
+                }
+
+         if(infoRC->GetStatus(1)==3) fMCRecPt->Fill(mcpt); 
+         if(infoRC->GetStatus(1)==3 && isPrim) fMCRecPrimPt->Fill(mcpt); 
+         if(infoRC->GetStatus(1)==3 && !isPrim) fMCRecSecPt->Fill(mcpt); 
+
+         if(isPrim)
+                {
+           fTPCPtDCASigmaIdeal[0]->Fill(mcpt,DCASigmaIdeal);
+           fTPCPtDCASigmaFull[0]->Fill(mcpt,DCASigmaFull);
+           fTPCPtDCASigmaDay0[0]->Fill(mcpt,DCASigmaDay0);
+
+           fTPCPtDCAXY[0]->Fill(mcpt,dca[0]);
+           fTPCPtDCAZ[0]->Fill(mcpt,dca[1]);
+
+           fTPCPtDCASigmaIdealPid[0]->Fill(mcpt,DCASigmaIdeal,pid);
+           fTPCPtDCASigmaFullPid[0]->Fill(mcpt,DCASigmaFull,pid);
+           fTPCPtDCASigmaDay0Pid[0]->Fill(mcpt,DCASigmaDay0,pid);
+
+           fTPCPtDCAXYPid[0]->Fill(mcpt,dca[0],pid);
+           fTPCPtDCAZPid[0]->Fill(mcpt,dca[1],pid);
+                          
+                  if(infoRC->GetStatus(1)==3)
+                  {
+             fTPCPtDCASigmaIdeal[1]->Fill(mcpt,DCASigmaIdeal);
+             fTPCPtDCASigmaFull[1]->Fill(mcpt,DCASigmaFull);
+             fTPCPtDCASigmaDay0[1]->Fill(mcpt,DCASigmaDay0);
+
+             fTPCPtDCAXY[1]->Fill(mcpt,dca[0]);
+             fTPCPtDCAZ[1]->Fill(mcpt,dca[1]);
+
+             fTPCPtDCASigmaIdealPid[1]->Fill(mcpt,DCASigmaIdeal,pid);
+             fTPCPtDCASigmaFullPid[1]->Fill(mcpt,DCASigmaFull,pid);
+             fTPCPtDCASigmaDay0Pid[1]->Fill(mcpt,DCASigmaDay0,pid);
+
+             fTPCPtDCAXYPid[1]->Fill(mcpt,dca[0],pid);
+             fTPCPtDCAZPid[1]->Fill(mcpt,dca[1],pid);
+           }
+                }
+
+        // Cont = infoRC->GetStatus(1)==3 && !isPrim / infoRC->GetStatus(1)==3   
+        // Pt vs ( dz[0]^2/cov[0]^2 + dz[1]^2/cov[2]^2 ) 
+        // Pt vs ( dz[0]^2/(cov[0]^2 + kSigma2Full_xy)  + dz[1]^2/( cov[2]^2 + kSigma2Full_z ) 
+        // Pt vs ( dz[0]^2/(cov[0]^2 + kSigma2_xy)  + dz[1]^2/( cov[2]^2 + kSigma2_z ) 
+
+                if(infoRC->GetStatus(1)==3) 
+                {
+                  fTPCPtDCASigmaIdeal[2]->Fill(mcpt,DCASigmaIdeal);
+           fTPCPtDCASigmaFull[2]->Fill(mcpt,DCASigmaFull);
+           fTPCPtDCASigmaDay0[2]->Fill(mcpt,DCASigmaDay0);
+
+           fTPCPtDCAXY[2]->Fill(mcpt,dca[0]);
+           fTPCPtDCAZ[2]->Fill(mcpt,dca[1]);
+
+           fTPCPtDCASigmaIdealPid[2]->Fill(mcpt,DCASigmaIdeal,pid);
+           fTPCPtDCASigmaFullPid[2]->Fill(mcpt,DCASigmaFull,pid);
+           fTPCPtDCASigmaDay0Pid[2]->Fill(mcpt,DCASigmaDay0,pid);
+
+           fTPCPtDCAXYPid[2]->Fill(mcpt,dca[0],pid);
+           fTPCPtDCAZPid[2]->Fill(mcpt,dca[1],pid);
+                  if(isPrim==0)
+                  {
+             fTPCPtDCASigmaIdeal[3]->Fill(mcpt,DCASigmaIdeal);
+             fTPCPtDCASigmaFull[3]->Fill(mcpt,DCASigmaFull);
+             fTPCPtDCASigmaDay0[3]->Fill(mcpt,DCASigmaDay0);
+
+             fTPCPtDCAXY[3]->Fill(mcpt,dca[0]);
+             fTPCPtDCAZ[3]->Fill(mcpt,dca[1]);
+
+             fTPCPtDCASigmaIdealPid[3]->Fill(mcpt,DCASigmaIdeal,pid);
+             fTPCPtDCASigmaFullPid[3]->Fill(mcpt,DCASigmaFull,pid);
+             fTPCPtDCASigmaDay0Pid[3]->Fill(mcpt,DCASigmaDay0,pid);
+
+             fTPCPtDCAXYPid[3]->Fill(mcpt,dca[0],pid);
+             fTPCPtDCAZPid[3]->Fill(mcpt,dca[1],pid);
+           }
+            }
+          delete track;
+       }
+     } 
+        else 
+        {
+       if(isPrim)
+          {
+            fTPCPtDCASigmaIdeal[0]->Fill(mcpt,0.0);
+         fTPCPtDCASigmaFull[0]->Fill(mcpt,0.0);
+         fTPCPtDCASigmaDay0[0]->Fill(mcpt,0.0);
+
+         fTPCPtDCAXY[0]->Fill(mcpt,0.0);
+         fTPCPtDCAZ[0]->Fill(mcpt,0.0);
+
+         fTPCPtDCASigmaIdealPid[0]->Fill(mcpt,0.0,pid);
+         fTPCPtDCASigmaFullPid[0]->Fill(mcpt,0.0,pid);
+         fTPCPtDCASigmaDay0Pid[0]->Fill(mcpt,0.0,pid);
+
+         fTPCPtDCAXYPid[0]->Fill(mcpt,0.0,pid);
+         fTPCPtDCAZPid[0]->Fill(mcpt,0.0,pid);
+          }
+     }
+  }
+
+  // only primary particles
+  if (!isPrim) return;
+
+  // pt
+  if (TMath::Abs(tantheta)<fCutsRC->GetMaxAbsTanTheta()){
+
+    fEffTPCPt->Fill(mcpt, infoRC->GetStatus(1)==3);
+    fEffTPCPtMC->Fill(mcpt, infoMC->GetRowsWithDigits()>fCutsMC->GetMinRowsWithDigits());
+    if (infoMC->GetRowsWithDigits()>fCutsMC->GetMinRowsWithDigits()){
+      fEffTPCPtF->Fill(mcpt, infoRC->GetStatus(1)==3);
+    }
+
+    // protons
+    if(TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetProt()) { 
+          fEffTPCPt_P->Fill(mcpt, infoRC->GetStatus(1)==3);
+       fEffTPCPtMC_P->Fill(mcpt, infoMC->GetRowsWithDigits()>fCutsMC->GetMinRowsWithDigits());
+
+       if(infoMC->GetRowsWithDigits()>fCutsMC->GetMinRowsWithDigits()) {
+         fEffTPCPtF_P->Fill(mcpt, infoRC->GetStatus(1)==3);
+       }
+       }
+
+    // pions
+    if(TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetPiP()) {
+         fEffTPCPt_Pi->Fill(mcpt, infoRC->GetStatus(1)==3);
+      fEffTPCPtMC_Pi->Fill(mcpt, infoMC->GetRowsWithDigits()>fCutsMC->GetMinRowsWithDigits());
+
+       if(infoMC->GetRowsWithDigits()>fCutsMC->GetMinRowsWithDigits()) {
+         fEffTPCPtF_Pi->Fill(mcpt, infoRC->GetStatus(1)==3);
+       }
+       }
+
+       // kaons
+    if(TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetKP()) {
+         fEffTPCPt_K->Fill(mcpt, infoRC->GetStatus(1)==3);
+      fEffTPCPtMC_K->Fill(mcpt, infoMC->GetRowsWithDigits()>fCutsMC->GetMinRowsWithDigits());
+
+       if(infoMC->GetRowsWithDigits()>fCutsMC->GetMinRowsWithDigits()) {
+         fEffTPCPtF_K->Fill(mcpt, infoRC->GetStatus(1)==3);
+       }
+       }
+  }
+
+  // theta
+  if (TMath::Abs(mcpt)>fCutsRC->GetPtMin()){
+    fEffTPCTan->Fill(tantheta, infoRC->GetStatus(1)==3);
+    fEffTPCTanMC->Fill(tantheta, infoMC->GetRowsWithDigits()>fCutsMC->GetMinRowsWithDigits());
+    if (infoMC->GetRowsWithDigits()>fCutsMC->GetMinRowsWithDigits()){
+      fEffTPCTanF->Fill(tantheta, infoRC->GetStatus(1)==3);
+    }
+  }
+
+  // pt-theta
+  fEffTPCPtTan->Fill(mcpt,tantheta,infoRC->GetStatus(1)==3);
+  fEffTPCPtTanMC->Fill(mcpt,tantheta,infoMC->GetRowsWithDigits()>fCutsMC->GetMinRowsWithDigits()); 
+  if (infoMC->GetRowsWithDigits()>fCutsMC->GetMinRowsWithDigits()){
+    fEffTPCPtTanF->Fill(mcpt,tantheta,infoRC->GetStatus(1)==3); 
+  }
+}
+
+//_____________________________________________________________________________
+void AliComparisonEff::Exec(AliMCInfo* infoMC, AliESDRecInfo *infoRC)
+{
+  // Process comparison information
+  Process(infoMC,infoRC);
+}
+
+//_____________________________________________________________________________
+Long64_t AliComparisonEff::Merge(TCollection* list) 
+{
+  // Merge list of objects (needed by PROOF)
+
+  if (!list)
+  return 0;
+
+  if (list->IsEmpty())
+  return 1;
+
+  TIterator* iter = list->MakeIterator();
+  TObject* obj = 0;
+
+  // collection of generated histograms
+
+  Int_t count=0;
+  while((obj = iter->Next()) != 0) 
+  {
+    AliComparisonEff* entry = dynamic_cast<AliComparisonEff*>(obj);
+    if (entry == 0) 
+      continue; 
+  
+    fMCPt->Add(entry->fMCPt);
+    fMCRecPt->Add(entry->fMCRecPt);
+    fMCRecPrimPt->Add(entry->fMCRecPrimPt);
+    fMCRecSecPt->Add(entry->fMCRecSecPt);
+
+    fEffTPCPt->Add(entry->fEffTPCPt);
+       fEffTPCPtMC->Add(entry->fEffTPCPtMC);
+       fEffTPCPtF->Add(entry->fEffTPCPtF);
+
+    fEffTPCPt_P->Add(entry->fEffTPCPt_P);
+       fEffTPCPtMC_P->Add(entry->fEffTPCPtMC_P);
+       fEffTPCPtF_P->Add(entry->fEffTPCPtF_P);
+
+    fEffTPCPt_Pi->Add(entry->fEffTPCPt_Pi);
+       fEffTPCPtMC_Pi->Add(entry->fEffTPCPtMC_Pi);
+       fEffTPCPtF_Pi->Add(entry->fEffTPCPtF_Pi);
+
+    fEffTPCPt_K->Add(entry->fEffTPCPt_K);
+       fEffTPCPtMC_K->Add(entry->fEffTPCPtMC_K);
+       fEffTPCPtF_K->Add(entry->fEffTPCPtF_K);
+
+       fEffTPCTan->Add(entry->fEffTPCTan);
+       fEffTPCTanMC->Add(entry->fEffTPCTanMC);
+       fEffTPCTanF->Add(entry->fEffTPCTanF);
+         
+       fEffTPCPtTan->Add(entry->fEffTPCPtTan);
+       fEffTPCPtTanMC->Add(entry->fEffTPCPtTanMC);
+       fEffTPCPtTanF->Add(entry->fEffTPCPtTanF);
+    
+    for(Int_t i=0; i<4; ++i)
+    {
+      fTPCPtDCASigmaIdeal[i]->Add(entry->fTPCPtDCASigmaIdeal[i]);
+      fTPCPtDCASigmaFull[i]->Add(entry->fTPCPtDCASigmaFull[i]);
+      fTPCPtDCASigmaDay0[i]->Add(entry->fTPCPtDCASigmaDay0[i]);
+
+      fTPCPtDCAXY[i]->Add(entry->fTPCPtDCAXY[i]);
+      fTPCPtDCAZ[i]->Add(entry->fTPCPtDCAXY[i]);
+
+      fTPCPtDCASigmaIdealPid[i]->Add(entry->fTPCPtDCASigmaIdealPid[i]);
+      fTPCPtDCASigmaFullPid[i]->Add(entry->fTPCPtDCASigmaFullPid[i]);
+      fTPCPtDCASigmaDay0Pid[i]->Add(entry->fTPCPtDCASigmaDay0Pid[i]);
+
+      fTPCPtDCAXYPid[i]->Add(entry->fTPCPtDCAXYPid[i]);
+      fTPCPtDCAZPid[i]->Add(entry->fTPCPtDCAXYPid[i]);
+    }
+
+  count++;
+  }
+
+return count;
+}
+//_____________________________________________________________________________
+void AliComparisonEff::Analyse() 
+{
+  // Analyse output histograms
+  
+  AliComparisonEff * comp=this;
+
+  // calculate efficiency and contamination (4 sigma) 
+
+  TH1 *h_sigmaidealpid[20];
+  TH1 *h_sigmafullpid[20];
+  TH1 *h_sigmaday0pid[20];
+
+  //TH1 *h_sigmaday0pidclone[20];
+
+  TH1 *h_sigmaidealpidtot[4];
+  TH1 *h_sigmafullpidtot[4];
+  TH1 *h_sigmaday0pidtot[4];
+
+  //TH1 *h_sigmaday0pidtotclone[4];
+
+  char name[256];
+  char name1[256];
+  Int_t idx;
+
+  for(Int_t i=0; i<4; ++i)
+  {
+     //total
+     comp->fTPCPtDCASigmaIdealPid[i]->GetYaxis()->SetRange(1,4);
+     sprintf(name,"h_sigmaidealpidtot_%d",i);
+     h_sigmaidealpidtot[i] = comp->fTPCPtDCASigmaIdealPid[i]->Project3D();
+     h_sigmaidealpidtot[i]->SetName(name);
+
+     comp->fTPCPtDCASigmaFullPid[i]->GetYaxis()->SetRange(1,4);
+     sprintf(name,"h_sigmafullpidtot_%d",i);
+     h_sigmafullpidtot[i] = comp->fTPCPtDCASigmaFullPid[i]->Project3D();
+     h_sigmafullpidtot[i]->SetName(name);
+
+     comp->fTPCPtDCASigmaDay0Pid[i]->GetYaxis()->SetRange(1,4);
+     sprintf(name,"h_sigmaday0pidtot_%d",i);
+     h_sigmaday0pidtot[i] = comp->fTPCPtDCASigmaDay0Pid[i]->Project3D();
+     h_sigmaday0pidtot[i]->SetName(name);
+
+     // pid wise
+     for(Int_t j=0; j<5; ++j)
+        {
+       idx = i*5 + j;
+
+       comp->fTPCPtDCASigmaIdealPid[i]->GetYaxis()->SetRange(1,4);
+       comp->fTPCPtDCASigmaIdealPid[i]->GetZaxis()->SetRange(j+1,j+1);
+
+       sprintf(name,"h_sigmaidealpid_%d",idx);
+       h_sigmaidealpid[idx] = comp->fTPCPtDCASigmaIdealPid[i]->Project3D();
+       h_sigmaidealpid[idx]->SetName(name);
+          
+
+       comp->fTPCPtDCASigmaFullPid[i]->GetYaxis()->SetRange(1,4);
+       comp->fTPCPtDCASigmaFullPid[i]->GetZaxis()->SetRange(j+1,j+1);
+
+       sprintf(name,"h_sigmafullpid_%d",idx);
+       h_sigmafullpid[idx] = comp->fTPCPtDCASigmaFullPid[i]->Project3D();
+       h_sigmafullpid[idx]->SetName(name);
+       
+       comp->fTPCPtDCASigmaDay0Pid[i]->GetYaxis()->SetRange(1,4);
+       comp->fTPCPtDCASigmaDay0Pid[i]->GetZaxis()->SetRange(j+1,j+1);
+
+       sprintf(name,"h_sigmaday0pid_%d",idx);
+       h_sigmaday0pid[idx] = comp->fTPCPtDCASigmaDay0Pid[i]->Project3D();
+       h_sigmaday0pid[idx]->SetName(name);
+
+       } 
+  }
+
+  // calculate efficiency and contamination (all pids)
+  h_sigmaidealpidtot[0]->Sumw2();
+  h_sigmaidealpidtot[1]->Divide(h_sigmaidealpidtot[0]);
+  h_sigmaidealpidtot[2]->Sumw2();
+  h_sigmaidealpidtot[3]->Divide(h_sigmaidealpidtot[2]);
+
+  h_sigmafullpidtot[0]->Sumw2();
+  h_sigmafullpidtot[1]->Divide(h_sigmafullpidtot[0]);
+  h_sigmafullpidtot[2]->Sumw2();
+  h_sigmafullpidtot[3]->Divide(h_sigmafullpidtot[2]);
+
+  h_sigmaday0pidtot[0]->Sumw2();
+  h_sigmaday0pidtot[1]->Divide(h_sigmaday0pidtot[0]);
+  h_sigmaday0pidtot[2]->Sumw2();
+  h_sigmaday0pidtot[3]->Divide(h_sigmaday0pidtot[2]);
+
+  // calculate efficiency pid wise
+  for(Int_t idx = 0; idx<5; idx++)
+  {
+    h_sigmaidealpid[idx]->Sumw2();
+    h_sigmaidealpid[idx+5]->Divide(h_sigmaidealpid[idx]);
+
+    h_sigmafullpid[idx]->Sumw2();
+    h_sigmafullpid[idx+5]->Divide(h_sigmafullpid[idx]);
+
+    h_sigmaday0pid[idx]->Sumw2();
+    h_sigmaday0pid[idx+5]->Divide(h_sigmaday0pid[idx]);
+  }
+
+  // calculate cont. pid wise
+  for(Int_t idx = 0; idx<5; idx++)
+  {
+    h_sigmaidealpid[idx+15]->Divide(h_sigmaidealpidtot[2]);
+    h_sigmafullpid[idx+15]->Divide(h_sigmafullpidtot[2]);
+    h_sigmaday0pid[idx+15]->Divide(h_sigmaday0pidtot[2]);
+  }
+
+  // write results
+  TFile *fp = new TFile("pictures_eff.root","recreate");
+  fp->cd();
+
+  TCanvas * c = new TCanvas("Efficiency","Track efficiency");
+  c->cd();
+
+  fMCPt->Write();
+  fMCRecPt->Write();
+  fMCRecPrimPt->Write();
+  fMCRecSecPt->Write();
+
+  for(int i = 0; i<4;i++)  
+  {
+   comp->fTPCPtDCASigmaIdealPid[i]->GetYaxis()->SetRange();
+   comp->fTPCPtDCASigmaIdealPid[i]->GetZaxis()->SetRange();
+   comp->fTPCPtDCASigmaFullPid[i]->GetYaxis()->SetRange();
+   comp->fTPCPtDCASigmaFullPid[i]->GetZaxis()->SetRange();
+   comp->fTPCPtDCASigmaDay0Pid[i]->GetYaxis()->SetRange();
+   comp->fTPCPtDCASigmaDay0Pid[i]->GetZaxis()->SetRange();
+
+    comp->fTPCPtDCASigmaIdealPid[i]->Write();
+    comp->fTPCPtDCASigmaFullPid[i]->Write();
+    comp->fTPCPtDCASigmaDay0Pid[i]->Write();
+  }
+  //
+  comp->fEffTPCTanF->SetXTitle("Tan(#theta)");
+  comp->fEffTPCTanF->SetYTitle("eff_{findable}");
+  comp->fEffTPCTanF->Write("EffTanFindable");
+  //
+  comp->fEffTPCTan->SetXTitle("Tan(#theta)");
+  comp->fEffTPCTan->SetYTitle("eff_{all}");
+  comp->fEffTPCTan->Write("EffTanAll");
+
+  h_sigmaidealpidtot[1]->Write("Eff_SigmaIdeal");
+  h_sigmaidealpidtot[3]->Write("Cont_SigmaIdeal");
+
+  h_sigmafullpidtot[1]->Write("Eff_SigmaFull");
+  h_sigmafullpidtot[3]->Write("Cont_SigmaFull");
+
+  h_sigmaday0pidtot[1]->Write("Eff_SigmaDay0");
+  h_sigmaday0pidtot[3]->Write("Cont_SigmaDay0");
+
+  for(Int_t idx = 0; idx<5; idx++)
+  {
+    sprintf(name,"Eff_SigmaIdeal_%d",idx);
+    sprintf(name1,"Cont_SigmaIdeal_%d",idx);
+
+    h_sigmaidealpid[idx+5]->Write(name);
+    h_sigmaidealpid[idx+15]->Write(name1);
+
+    sprintf(name,"Eff_SigmaFull_%d",idx);
+    sprintf(name1,"Cont_SigmaFull_%d",idx);
+
+    h_sigmafullpid[idx+5]->Write(name);
+    h_sigmafullpid[idx+15]->Write(name1);
+
+    sprintf(name,"Eff_SigmaDay0_%d",idx);
+    sprintf(name1,"Cont_SigmaDay0_%d",idx);
+
+    h_sigmaday0pid[idx+5]->Write(name);
+    h_sigmaday0pid[idx+15]->Write(name1);
+  }
+
+  //
+  fp->Close();
+}
diff --git a/PWG1/AliComparisonEff.h b/PWG1/AliComparisonEff.h
new file mode 100644 (file)
index 0000000..6f7ccbb
--- /dev/null
@@ -0,0 +1,120 @@
+#ifndef ALICOMPARISONEFF_H
+#define ALICOMPARISONEFF_H
+
+//------------------------------------------------------------------------------
+// Class to keep information from comparison of 
+// reconstructed and MC particle tracks (TPC efficiency).   
+// 
+// Author: J.Otwinowski 04/02/2008 
+//------------------------------------------------------------------------------
+
+class TFile;
+class AliMCInfo;
+class AliESDRecInfo;
+class AliESDEvent; 
+class AliESD;
+class AliESDfriend;
+class AliRecInfoCuts;
+class AliMCInfoCuts;
+class TH1I;
+class TH3F;
+class TH3;
+class TProfile;
+class TProfile2D;
+class TGraph2D;
+class TGraph; 
+class TGeoManager; 
+class TStatToolkit; 
+class AliMagFMaps;
+class AliESDVertex;
+
+#include "TNamed.h"
+
+class AliComparisonEff : public TNamed {
+public :
+  AliComparisonEff(); 
+  ~AliComparisonEff();
+  void InitHisto();
+  void InitCuts();
+  void Exec(AliMCInfo* infoMC, AliESDRecInfo *infoRC);
+  void Process(AliMCInfo* infoMC, AliESDRecInfo *infoRC);
+  
+  // Selection cuts
+  void SetAliRecInfoCuts(AliRecInfoCuts* cuts=0) {fCutsRC = cuts;}
+  void SetAliMCInfoCuts(AliMCInfoCuts* cuts=0) {fCutsMC = cuts;} 
+  
+  AliRecInfoCuts*  GetAliRecInfoCuts() const {return fCutsRC;} 
+  AliMCInfoCuts*   GetAliMCInfoCuts()  const {return fCutsMC;}
+
+  // Merge output objects (needed by PROOF) 
+  virtual Long64_t Merge(TCollection* list);
+
+  // Analyse output histograms 
+  void Analyse();
+
+private:
+
+  // Control histograms
+  TH1F *fMCPt;
+  TH1F *fMCRecPt;
+  TH1F *fMCRecPrimPt;
+  TH1F *fMCRecSecPt;
+  
+  TProfile* fEffTPCPt;      //->TPC efficiency as function of Pt (tan+-1)
+  TProfile* fEffTPCPtMC;    //->MC -TPC efficiency as function of Pt (tan+-1)
+  TProfile* fEffTPCPtF;     //->efficiency for findable tracks
+
+  TProfile* fEffTPCPt_P;    //->TPC efficiency as function of Pt (tan+-1) - Protons
+  TProfile* fEffTPCPtMC_P;  //->MC -TPC efficiency as function of Pt (tan+-1) - Protons
+  TProfile* fEffTPCPtF_P;   //->efficiency for findable tracks - Protons
+
+  TProfile* fEffTPCPt_Pi;   //->TPC efficiency as function of Pt (tan+-1) - Pions
+  TProfile* fEffTPCPtMC_Pi; //->MC -TPC efficiency as function of Pt (tan+-1) - Pions
+  TProfile* fEffTPCPtF_Pi;  //->efficiency for findable tracks - Pions
+
+  TProfile* fEffTPCPt_K;    //->TPC efficiency as function of Pt (tan+-1) - Kaons
+  TProfile* fEffTPCPtMC_K;  //->MC -TPC efficiency as function of Pt (tan+-1) - Kaons
+  TProfile* fEffTPCPtF_K;   //->efficiency for findable tracks - Kaons
+
+  //
+  TProfile* fEffTPCTan;     //->TPC efficiency as function of Tan (pt>0.15
+  TProfile* fEffTPCTanMC;   //->MC -TPC efficiency as function of Tan (pt>0.15)
+  TProfile* fEffTPCTanF;    //->efficiency for findable tracks Tan (pt>0.15)
+  //
+  TProfile2D* fEffTPCPtTan;    //->TPC efficiency as function of Pt and tan
+  TProfile2D* fEffTPCPtTanMC;  //->MC -TPC efficiency as function of Pt and tan
+  TProfile2D* fEffTPCPtTanF;   //->TPC efficiency as function of Pt and tan
+
+  // idx - 0 (isPrim), idx - 1 (isPrim && infoRC->GetStatus(1)==3)
+  // idx - 2 (infoRC->GetStatus(1)==3),  idx - 3 (infoRC->GetStatus(1)==3 && !isPrim )
+  //
+  
+  TH2F* fTPCPtDCASigmaIdeal[4]; //->TPC efficiency vs Pt vs DCA/Sigma (tan+-1)
+  TH2F* fTPCPtDCASigmaFull[4];  //->TPC efficiency vs Pt vs DCA/Sigma (tan+-1, full systematics)
+  TH2F* fTPCPtDCASigmaDay0[4];  //->TPC efficiency vs Pt vs DCA/Sigma (tan+-1, goofie systematics)
+
+  TH2F* fTPCPtDCAXY[4];     //->TPC efficiency as Pt vs DCA_XY (tan+-1)
+  TH2F* fTPCPtDCAZ[4];      //->TPC efficiency as Pt vs DCA_Z (tan+-1)
+
+  // Pid = 0 - electrons,  1 - muons, 2 - kaons, 3 - pions, 4 - protons   
+  TH3F* fTPCPtDCASigmaIdealPid[4]; //->TPC efficiency vs Pt vs DCA/Sigma (tan+-1)
+  TH3F* fTPCPtDCASigmaFullPid[4];  //->TPC efficiency vs Pt vs DCA/Sigma (tan+-1, full systematics)
+  TH3F* fTPCPtDCASigmaDay0Pid[4];  //->TPC efficiency vs Pt vs DCA/Sigma (tan+-1, goofie systematics)
+
+  TH3F* fTPCPtDCAXYPid[4];     //->TPC efficiency as Pt vs DCA_XY (tan+-1)
+  TH3F* fTPCPtDCAZPid[4];      //->TPC efficiency as Pt vs DCA_Z (tan+-1)
+
+  // Global cuts objects
+  AliRecInfoCuts* fCutsRC;     // selection cuts for reconstructed tracks
+  AliMCInfoCuts*  fCutsMC;     // selection cuts for MC tracks
+
+  // Magnet (needed for DCA calculations) 
+  AliESDVertex* fVertex;  //! 
+  
+  AliComparisonEff(const AliComparisonEff&); // not implemented
+  AliComparisonEff& operator=(const AliComparisonEff&); // not implemented
+
+  ClassDef(AliComparisonEff,1);
+};
+
+#endif
diff --git a/PWG1/AliComparisonRes.cxx b/PWG1/AliComparisonRes.cxx
new file mode 100644 (file)
index 0000000..70593fe
--- /dev/null
@@ -0,0 +1,358 @@
+//------------------------------------------------------------------------------
+// Implementation of AliComparisonRes class. It keeps information from 
+// comparison of reconstructed and MC particle tracks. In addtion, 
+// it keeps selection cuts used during comparison. The comparison 
+// information is stored in the ROOT histograms. Analysis of these 
+// histograms can be done by using Analyse() class function. The result of 
+// the analysis (histograms) are stored in the output picture_res.root file.
+//  
+// Author: J.Otwinowski 04/02/2008 
+//------------------------------------------------------------------------------
+
+/*
+  //after running analysis, read the file, and get component
+  gSystem->Load("libPWG1.so");
+  TFile f("Output.root");
+  AliComparisonRes * comp = (AliComparisonRes*)f.Get("AliComparisonRes");
+
+  // analyse comparison data (output stored in pictures_res.root)
+  comp->Analyse();
+  
+  // paramtetrisation of the TPC track length (for information only) 
+  TF1 fl("fl","((min(250./(abs(x+0.000001)),250)-90))",0,2);  // TPC track length function
+  TF1 fl2("fl2","[0]/((min(250./(abs(x+0.000001)),250)-90))^[1]",0,2);
+  fl2.SetParameter(1,1);
+  fl2.SetParameter(0,1);
+*/
+
+#include <iostream>
+
+#include "TFile.h"
+#include "TCint.h"
+#include "TH3F.h"
+#include "TH2F.h"
+#include "TF1.h"
+#include "TProfile.h"
+#include "TProfile2D.h"
+#include "TGraph2D.h"
+#include "TCanvas.h"
+#include "TGraph.h"
+
+#include "AliESDEvent.h"   
+#include "AliESD.h"
+#include "AliESDfriend.h"
+#include "AliESDfriendTrack.h"
+#include "AliRecInfoCuts.h" 
+#include "AliMCInfoCuts.h" 
+#include "AliLog.h" 
+
+#include "AliMathBase.h"
+#include "AliTreeDraw.h" 
+
+#include "AliMCInfo.h" 
+#include "AliESDRecInfo.h" 
+#include "AliComparisonRes.h" 
+
+using namespace std;
+
+ClassImp(AliComparisonRes)
+
+//_____________________________________________________________________________
+AliComparisonRes::AliComparisonRes():
+  TNamed("AliComparisonRes","AliComparisonRes"),
+
+  // Resolution 
+  fPtResolLPT(0),        // pt resolution - low pt
+  fPtResolHPT(0),        // pt resolution - high pt 
+  fPtPullLPT(0),         // pt resolution - low pt
+  fPtPullHPT(0),         // pt resolution - high pt 
+  //
+  // Resolution constrained param
+  //
+  fCPhiResolTan(0),  // angular resolution -  constrained
+  fCTanResolTan(0),  // angular resolution -  constrained
+  fCPtResolTan(0),   // pt resolution      -  constrained
+  fCPhiPullTan(0),   // angular resolution -  constrained
+  fCTanPullTan(0),   // angular resolution -  constrained
+  fCPtPullTan(0),    // pt resolution      -  constrained
+  // Cuts 
+  fCutsRC(0),  
+  fCutsMC(0)  
+{
+  InitHisto();
+  InitCuts();
+}
+
+//_____________________________________________________________________________
+AliComparisonRes::~AliComparisonRes(){
+  
+  // Resolution histograms
+  if(fPtResolLPT) delete  fPtResolLPT; fPtResolLPT=0;     
+  if(fPtResolHPT) delete  fPtResolHPT; fPtResolHPT=0;    
+  if(fPtPullLPT)  delete  fPtPullLPT;  fPtPullLPT=0;    
+  if(fPtPullHPT)  delete  fPtPullHPT;  fPtPullHPT=0;   
+
+  // Resolution histograms (constrained param)
+  if(fCPhiResolTan) delete fCPhiResolTan; fCPhiResolTan=0;
+  if(fCTanResolTan) delete fCTanResolTan; fCTanResolTan=0;
+  if(fCPtResolTan)  delete fCPtResolTan;  fCPtResolTan=0; 
+  if(fCPhiPullTan)  delete fCPhiPullTan;  fCPhiPullTan=0;
+  if(fCTanPullTan)  delete fCTanPullTan;  fCTanPullTan=0;
+  if(fCPtPullTan)   delete fCPtPullTan;   fCPtPullTan=0;
+}
+
+//_____________________________________________________________________________
+void AliComparisonRes::InitHisto(){
+
+  // Init histograms
+  fCPhiResolTan = new TH2F("CPhiResolTan","CPhiResolTan",50, -2,2,200,-0.025,0.025);   
+  fCPhiResolTan->SetXTitle("tan(#theta)");
+  fCPhiResolTan->SetYTitle("#Delta#phi");
+
+  fCTanResolTan = new TH2F("CTanResolTan","CTanResolTan",50, -2,2,200,-0.025,0.025);
+  fCTanResolTan->SetXTitle("tan(#theta)");
+  fCTanResolTan->SetYTitle("#Delta#theta");
+
+  fCPtResolTan=new TH2F("CPtResol","CPtResol",50, -2,2,200,-0.2,0.2);    
+  fCPtResolTan->SetXTitle("Tan(#theta)");
+  fCPtResolTan->SetYTitle("#Deltap_{t}/p_{t}");
+
+  fCPhiPullTan = new TH2F("CPhiPullTan","CPhiPullTan",50, -2,2,200,-5,5);   
+  fCPhiPullTan->SetXTitle("Tan(#theta)");
+  fCPhiPullTan->SetYTitle("#Delta#phi/#Sigma");
+
+  fCTanPullTan = new TH2F("CTanPullTan","CTanPullTan",50, -2,2,200,-5,5);
+  fCTanPullTan->SetXTitle("Tan(#theta)");
+  fCTanPullTan->SetYTitle("#Delta#theta/#Sigma");
+
+  fCPtPullTan=new TH2F("CPtPull","CPtPull",50, -2,2,200,-5,5);    
+  fCPtPullTan->SetXTitle("Tan(#theta)");
+  fCPtPullTan->SetYTitle("(1/mcp_{t}-1/p_{t})/#Sigma");
+
+  fPtResolLPT = new TH2F("Pt_resol_lpt","pt resol",10, 0.1,3,200,-0.2,0.2);
+  fPtResolLPT->SetXTitle("p_{t}");
+  fPtResolLPT->SetYTitle("#Deltap_{t}/p_{t}");
+
+  fPtResolHPT = new TH2F("Pt_resol_hpt","pt resol",10, 2,100,200,-0.3,0.3);  
+  fPtResolHPT->SetXTitle("p_{t}");
+  fPtResolHPT->SetYTitle("#Deltap_{t}/p_{t}");
+
+  fPtPullLPT = new TH2F("Pt_pull_lpt","pt pool",10, 0.1,3,200,-6,6);
+  fPtPullLPT->SetXTitle("p_{t}");
+  fPtPullLPT->SetYTitle("#Deltap_{t}/#Sigma");
+
+  fPtPullHPT = new TH2F("Pt_pull_hpt","pt pool",10, 2,100,200,-6,6);  
+  fPtPullHPT->SetXTitle("p_{t}");
+  fPtPullHPT->SetYTitle("#Deltap_{t}/#Sigma");
+}
+
+//_____________________________________________________________________________
+void AliComparisonRes::InitCuts()
+{
+
+ /*
+    fCutsRC = new AliRecInfoCuts();
+  if(fCutsRC) {
+    fCutsRC->SetPtRange(0.15,200.0);
+    fCutsRC->SetMaxAbsTanTheta(1.0);
+    fCutsRC->SetMinNClustersTPC(10);
+    fCutsRC->SetMinTPCsignalN(50);
+  } else {
+    AliDebug(AliLog::kError, "ERROR: Cannot create AliRecInfoCuts object");
+    return;
+  }
+
+  fCutsMC = new AliMCInfoCuts();
+  if(fCutsMC) {
+    fCutsMC->SetMinRowsWithDigits(50);
+    fCutsMC->SetMaxR(0.1);
+    fCutsMC->SetMaxVz(10.);
+    fCutsMC->SetRangeTPCSignal(0.5,1.4);
+  } else {
+  AliDebug(AliLog::kError, "ERROR: Cannot AliMCInfoCuts object");
+  return;
+  }
+ */
+  // Init cuts 
+  if(!fCutsMC) 
+    AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
+  if(!fCutsRC) 
+    AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
+}
+
+//_____________________________________________________________________________
+void AliComparisonRes::Process(AliMCInfo* infoMC, AliESDRecInfo *infoRC)
+{
+  // Fill resolution comparison information 
+
+  Float_t mcpt = infoMC->GetParticle().Pt();
+  Bool_t isPrim = infoMC->GetParticle().R() < fCutsMC->GetMaxR() && 
+                  TMath::Abs(infoMC->GetParticle().Vz()) < fCutsMC->GetMaxVz() ;
+
+  // Check selection cuts
+  if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return; 
+  if (!isPrim) return;
+  if (infoRC->GetStatus(1)==0) return;
+  if (!infoRC->GetESDtrack()) return;  
+  if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;
+  
+  Float_t deltaPt= (mcpt-infoRC->GetESDtrack()->Pt())/mcpt;  
+  Float_t pullPt= (1/mcpt-infoRC->GetESDtrack()->OneOverPt())/TMath::Sqrt(infoRC->GetESDtrack()->GetSigma1Pt2());  
+
+  // Fill histograms
+  fPtResolLPT->Fill(mcpt,deltaPt);
+  fPtResolHPT->Fill(mcpt,deltaPt);
+  fPtPullLPT->Fill(mcpt,pullPt);
+  fPtPullHPT->Fill(mcpt,pullPt);  
+}
+
+Long64_t AliComparisonRes::Merge(TCollection* list) 
+{
+  // Merge list of objects (needed by PROOF)
+
+  if (!list)
+  return 0;
+
+  if (list->IsEmpty())
+  return 1;
+
+  TIterator* iter = list->MakeIterator();
+  TObject* obj = 0;
+
+  // collection of generated histograms
+  Int_t count=0;
+  while((obj = iter->Next()) != 0) 
+  {
+  AliComparisonRes* entry = dynamic_cast<AliComparisonRes*>(obj);
+  if (entry == 0) { 
+    Error("Add","Attempt to add object of class: %s to a %s",
+       obj->ClassName(),this->ClassName());
+    return -1;
+  }
+
+  fPtResolLPT->Add(entry->fPtResolLPT);
+  fPtResolHPT->Add(entry->fPtResolHPT);
+  fPtPullLPT->Add(entry->fPtPullLPT);
+  fPtPullHPT->Add(entry->fPtPullHPT);
+
+  // Resolution histograms (constrained param)
+  fCPhiResolTan->Add(entry->fCPhiResolTan);
+  fCTanResolTan->Add(entry->fCTanResolTan);
+  fCPtResolTan->Add(entry->fCPtResolTan);
+  fCPhiPullTan->Add(entry->fCPhiPullTan);
+  fCTanPullTan->Add(entry->fCTanPullTan);
+  fCPtPullTan->Add(entry->fCPtPullTan);
+
+  count++;
+  }
+
+return count;
+}
+
+
+
+
+
+//_____________________________________________________________________________
+void AliComparisonRes::ProcessConstrained(AliMCInfo* infoMC, AliESDRecInfo *infoRC)
+{
+  // Fill resolution comparison information (constarained parameters) 
+  
+  Float_t mcpt = infoMC->GetParticle().Pt();
+  Float_t tantheta = TMath::Tan(infoMC->GetParticle().Theta()-TMath::Pi()*0.5);
+  Bool_t isPrim = infoMC->GetParticle().R()<fCutsMC->GetMaxR() && TMath::Abs(infoMC->GetParticle().Vz())<fCutsMC->GetMaxVz();
+  
+  // Check selection cuts
+  if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return; 
+  if (!isPrim) return;
+  if (infoRC->GetStatus(1)!=3) return;
+  if (!infoRC->GetESDtrack()) return;  
+  if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;
+  if (!infoRC->GetESDtrack()->GetConstrainedParam()) return;
+  
+  // constrained parameters resolution
+  const AliExternalTrackParam * cparam = infoRC->GetESDtrack()->GetConstrainedParam();
+  Float_t deltaCPt= (mcpt-cparam->Pt())/mcpt;  
+  Float_t pullCPt= (1/mcpt-cparam->OneOverPt())/
+    TMath::Sqrt(cparam->GetSigma1Pt2());          
+  Float_t deltaPhi = TMath::ATan2(cparam->Py(),cparam->Px())-
+    TMath::ATan2(infoMC->GetParticle().Py(),infoMC->GetParticle().Px());
+  Float_t pullPhi = deltaPhi/TMath::Sqrt(cparam->GetSigmaSnp2()); 
+
+  Float_t deltaTan = TMath::ATan2(cparam->Pz(),cparam->Pt())-TMath::ATan2(infoMC->GetParticle().Pz(),infoMC->GetParticle().Pt());
+  Float_t pullTan = deltaPhi/TMath::Sqrt(cparam->GetSigmaSnp2()); 
+
+  // Fill histograms
+  fCPtResolTan->Fill(tantheta,deltaCPt);
+  fCPtPullTan->Fill(tantheta,pullCPt);
+  fCPhiResolTan->Fill(tantheta,deltaPhi);
+  fCPhiPullTan->Fill(tantheta,pullPhi);
+  fCTanResolTan->Fill(tantheta,deltaTan);
+  fCTanPullTan->Fill(tantheta,pullTan);
+}
+//_____________________________________________________________________________
+void AliComparisonRes::Exec(AliMCInfo* infoMC, AliESDRecInfo *infoRC){
+  
+  // Process comparison information 
+  Process(infoMC,infoRC);
+  ProcessConstrained(infoMC,infoRC);
+
+}
+
+//_____________________________________________________________________________
+TH1F* AliComparisonRes::MakeResol(TH2F * his, Int_t integ, Bool_t type){
+  // Create resolution histograms
+  TH1F *hisr, *hism;
+  if (!gPad) new TCanvas;
+  hisr = AliTreeDraw::CreateResHistoI(his,&hism,integ);
+  if (type) return hism;
+  else 
+    return hisr;
+}
+
+//_____________________________________________________________________________
+void AliComparisonRes::Analyse(){
+  // Analyse comparison information and store output histograms 
+  // in the "pictures_res.root" file 
+  
+  AliComparisonRes * comp=this;
+  TH1F *hiss=0;
+
+  TFile *fp = new TFile("pictures_res.root","recreate");
+  fp->cd();
+
+  TCanvas * c = new TCanvas("Phi resol Tan","Phi resol Tan");
+  c->cd();
+  //
+  hiss = comp->MakeResol(comp->fCPtResolTan,1,0);
+  hiss->SetXTitle("Tan(#theta)");
+  hiss->SetYTitle("#sigmap_{t}/p_{t}");
+  hiss->Draw(); 
+  hiss->Write("CptResolTan");
+  //
+  hiss = comp->MakeResol(comp->fCPhiResolTan,1,0);
+  hiss->SetXTitle("Tan(#theta)");
+  hiss->SetYTitle("#sigma#phi (rad)");
+  hiss->Draw();
+  hiss->Write("PhiResolTan");
+  //
+  hiss = comp->MakeResol(comp->fCTanResolTan,1,0);
+  hiss->SetXTitle("Tan(#theta)");
+  hiss->SetYTitle("#sigma#theta (rad)");
+  hiss->Draw();
+  hiss->Write("ThetaResolTan");
+  //
+  hiss = comp->MakeResol(comp->fCPtPullTan,1,0);
+  hiss->SetXTitle("Tan(#theta)");
+  hiss->SetYTitle("1/mcp_{t}-1/p_{t}/#Sigma(1/p_{t})");
+  hiss->Draw();
+  hiss->Write("CptPullTan");
+
+  fp->Close();
+}
+
+
+
diff --git a/PWG1/AliComparisonRes.h b/PWG1/AliComparisonRes.h
new file mode 100644 (file)
index 0000000..df1e3fc
--- /dev/null
@@ -0,0 +1,80 @@
+#ifndef ALICOMPARISONRES_H
+#define ALICOMPARISONRES_H
+
+//------------------------------------------------------------------------------
+// Class to keep information from comparison of 
+// reconstructed and MC particle tracks (TPC resolution).   
+// 
+// Author: J.Otwinowski 04/02/2008 
+//------------------------------------------------------------------------------
+
+class TFile;
+class AliMCInfo;
+class AliESDRecInfo;
+class AliESDEvent; 
+class AliESD;
+class AliESDfriend;
+class AliMCInfoCuts;
+class AliRecInfoCuts;
+class TH1I;
+class TH3F;
+class TH3;
+class TProfile;
+class TProfile2D;
+
+#include "TNamed.h"
+
+class AliComparisonRes : public TNamed {
+public :
+  AliComparisonRes(); 
+  virtual ~AliComparisonRes();
+  void      InitHisto();
+  void      InitCuts();
+  void      Exec(AliMCInfo* infoMC, AliESDRecInfo *infoRC);
+  void      ProcessConstrained(AliMCInfo* infoMC, AliESDRecInfo *infoRC);
+  void      Process(AliMCInfo* infoMC, AliESDRecInfo *infoRC);
+
+  // Selection cuts
+  void SetAliRecInfoCuts(AliRecInfoCuts* cuts=0) {fCutsRC = cuts;}   
+  void SetAliMCInfoCuts(AliMCInfoCuts* cuts=0) {fCutsMC = cuts;}  
+   
+  AliRecInfoCuts*  GetAliRecInfoCuts() const {return fCutsRC;}  
+  AliMCInfoCuts*   GetAliMCInfoCuts()  const {return fCutsMC;}  
+
+  // Merge output objects (needed by PROOF) 
+  virtual Long64_t Merge(TCollection* list);
+
+  // Analyse output histograms
+  void Analyse();
+  static TH1F*       MakeResol(TH2F * his, Int_t integ, Bool_t type); 
+
+
+private:
+  //
+  // Control histograms
+  //
+  TH2F* fPtResolLPT;        //-> pt resolution - low pt
+  TH2F* fPtResolHPT;        //-> pt resolution - high pt 
+  TH2F* fPtPullLPT;         //-> pt resolution - low pt
+  TH2F* fPtPullHPT;         //-> pt resolution - high pt 
+  //
+  // Resolution constrained param
+  //
+  TH2F   *fCPhiResolTan;   //-> angular resolution -  constrained
+  TH2F   *fCTanResolTan;   //-> angular resolution -  constrained
+  TH2F   *fCPtResolTan;    //-> pt resolution      -  constrained
+  TH2F   *fCPhiPullTan;    //-> angular resolution -  constrained
+  TH2F   *fCTanPullTan;    //-> angular resolution -  constrained
+  TH2F   *fCPtPullTan;     //-> pt resolution      -  constrained
+
+  // Global cuts objects
+  AliRecInfoCuts*  fCutsRC;      // selection cuts for reconstructed tracks
+  AliMCInfoCuts*  fCutsMC;       // selection cuts for MC tracks
+
+  AliComparisonRes(const AliComparisonRes&); // not implemented
+  AliComparisonRes& operator=(const AliComparisonRes&); // not implemented
+
+  ClassDef(AliComparisonRes,1);
+};
+
+#endif
diff --git a/PWG1/AliComparisonTask.cxx b/PWG1/AliComparisonTask.cxx
new file mode 100644 (file)
index 0000000..def4f22
--- /dev/null
@@ -0,0 +1,223 @@
+//------------------------------------------------------------------------------
+// Implementation of the AliComparisonTask class. It compares properties of the 
+// reconstructed and MC particle tracks under several conditions. 
+// As the input it requires the TTree with AliRecInfo and AliMCInfo branches. 
+// The comparison output histograms are stored 
+// in the comparison objects: AliComparisonRes, AliComparisonEff, 
+// AliComparisonDEdx and AliComparisonDCA. Each of these objects also contains 
+// selection cuts which were used during filling the histograms.
+// 
+// Author: J.Otwinowski 04/02/2008 
+//------------------------------------------------------------------------------
+
+#include "iostream"
+
+#include "TChain.h"
+#include "TTree.h"
+#include "TH1F.h"
+#include "TCanvas.h"
+#include "TList.h"
+#include "TFile.h"
+
+#include "AliAnalysisTask.h"
+#include "AliAnalysisManager.h"
+#include "AliESDEvent.h"
+#include "AliESDInputHandler.h"
+#include "AliESDVertex.h"
+#include "AliMagFMaps.h"
+#include "AliTracker.h"
+#include "AliGeomManager.h"
+
+#include "AliMCInfo.h"
+#include "AliESDRecInfo.h"
+#include "AliMCInfoCuts.h"
+#include "AliRecInfoCuts.h"
+#include "AliComparisonRes.h"
+#include "AliComparisonEff.h"
+#include "AliComparisonDEdx.h"
+#include "AliComparisonDCA.h"
+#include "AliComparisonTask.h"
+
+using namespace std;
+
+ClassImp(AliComparisonTask)
+
+Int_t AliComparisonTask::evtNumber = 0;
+
+//_____________________________________________________________________________
+AliComparisonTask::AliComparisonTask(const char *name) 
+  : AliAnalysisTask(name, "")
+  , fTree(0)
+  , fInfoMC(0)
+  , fInfoRC(0)
+  , fCompRes(0)
+  , fCompEff(0)
+  , fCompDEdx(0)
+  , fCompDCA(0)
+  , fOutput(0)
+  , fMagField(0)
+  , fMagFMap(0)
+  , fGeom(0)
+{
+  // Constructor
+
+  // Define input and output slots here
+  DefineInput(0, TChain::Class());
+  DefineOutput(0, TList::Class());
+
+  // set default mag. field
+  SetMagField();
+  
+  // set default geometry
+  SetGeometry();
+}
+
+//_____________________________________________________________________________
+AliComparisonTask::~AliComparisonTask()
+{
+  if(fOutput)   delete fOutput;  fOutput =0; 
+  if(fMagFMap)  delete fMagFMap;  fMagFMap =0; 
+}
+
+//_____________________________________________________________________________
+void AliComparisonTask::ConnectInputData(Option_t *) 
+{
+  // Connect input data 
+  // Called once
+
+  fTree = dynamic_cast<TTree*> (GetInputData(0));
+  if (!fTree) {
+    Printf("ERROR: Could not read chain from input slot 0");
+  } else {
+    fTree->SetBranchStatus("*",1);
+  }
+
+  if(fTree->GetBranch("MC") &&  fTree->GetBranch("RC")) {
+    fTree->GetBranch("MC")->SetAddress(&fInfoMC);
+    fTree->GetBranch("RC")->SetAddress(&fInfoRC);
+  } else {
+      Printf("ERROR: Could not get MC and RC branches");
+  }
+  
+  // set mag. field map 
+  fMagFMap = new AliMagFMaps("Maps","Maps", 2, 1., 10., fMagField);
+  AliTracker::SetFieldMap(fMagFMap,kFALSE);
+
+  // set geommetry
+  AliGeomManager::LoadGeometry(fGeom);
+}
+
+//_____________________________________________________________________________
+void AliComparisonTask::CreateOutputObjects()
+{
+  // Create histograms
+  // Called once
+  fOutput = new TList;
+  fOutput->SetOwner();
+
+  if(fCompRes) fOutput->Add(fCompRes);
+  else 
+     Printf("WARNING: AliComparisonRes is not added to the output");
+
+  if(fCompEff) fOutput->Add(fCompEff);
+  else 
+    Printf("WARNING: AliComparisonEff is not added to the output");
+
+  if(fCompDEdx) fOutput->Add(fCompDEdx);
+  else 
+    Printf("WARNING: AliComparisonDEdx is not added to the output");
+
+  if(fCompDCA) fOutput->Add(fCompDCA);
+  else 
+     Printf("WARNING: AliComparisonDCA is not added to the output");
+}
+
+//_____________________________________________________________________________
+Bool_t AliComparisonTask::ReadEntry(Int_t evt) 
+{
+  Long64_t centry = fTree->LoadTree(evt);
+  if(centry < 0) return kFALSE;
+
+  if(fTree->GetBranch("MC") &&  fTree->GetBranch("RC")) {
+    fTree->GetBranch("MC")->SetAddress(&fInfoMC);
+    fTree->GetBranch("RC")->SetAddress(&fInfoRC);
+  } else {
+      Printf("ERROR: Could not get MC and RC branches");
+         return kFALSE;
+  }
+  fTree->GetEntry(evt);
+
+return kTRUE;
+}
+//_____________________________________________________________________________
+void AliComparisonTask::Exec(Option_t *) 
+{
+  // Main loop
+  // Called for each event
+
+  if (!fInfoMC && !fInfoRC) {
+    Printf("ERROR: fInfoMC && fInfoRC not available");
+    return;
+  }
+
+  // Process comparison
+  Bool_t status = ReadEntry(evtNumber);
+  if(status == kTRUE) 
+  {
+     if(fCompRes)  fCompRes->Exec(fInfoMC,fInfoRC);
+     if(fCompEff)  fCompEff->Exec(fInfoMC,fInfoRC);
+     if(fCompDEdx) fCompDEdx->Exec(fInfoMC,fInfoRC);
+     if(fCompDCA)  fCompDCA->Exec(fInfoMC,fInfoRC);
+  }
+
+  if( !( evtNumber % 10000) ) { 
+    cout << evtNumber << endl;
+  }
+  evtNumber++;
+
+  // Post output data.
+  PostData(0, fOutput);
+}
+
+//_____________________________________________________________________________
+void AliComparisonTask::Terminate(Option_t *) 
+{
+  // Called once at the end of the event loop
+  cout << "Terminate " << endl;
+
+  TFile *out = new TFile("Output.root","RECREATE");
+  out->cd();
+
+  fOutput = dynamic_cast<TList*> (GetOutputData(0));
+  if (!fOutput) {
+    Printf("ERROR: fOutput not available");
+    return;
+  }
+
+  fCompRes = dynamic_cast<AliComparisonRes*> (fOutput->FindObject("AliComparisonRes"));
+  if (!fCompRes) {
+    Printf("WARNING: AliComparisonRes not available");
+    return;
+  }
+
+  fCompEff = dynamic_cast<AliComparisonEff*> (fOutput->FindObject("AliComparisonEff"));
+  if (!fCompEff) {
+    Printf("WARNING: AliComparisonEff not available");
+    return;
+  }
+   
+  fCompDEdx = dynamic_cast<AliComparisonDEdx*> (fOutput->FindObject("AliComparisonDEdx"));
+  if (!fCompDEdx) {
+    Printf("WARNING: AliComparisonDEdx not available");
+    return;
+  }
+
+  fCompDCA = dynamic_cast<AliComparisonDCA*> (fOutput->FindObject("AliComparisonDCA"));
+  if (!fCompDCA) {
+    Printf("WARNING: AliComparisonDCA not available");
+    return;
+  }
+
+  fOutput->Write();
+  out->Close();
+}
diff --git a/PWG1/AliComparisonTask.h b/PWG1/AliComparisonTask.h
new file mode 100644 (file)
index 0000000..bf5537a
--- /dev/null
@@ -0,0 +1,62 @@
+#ifndef ALICOMPARISONRESTASK_H
+#define ALICOMPARISONRESTASK_H
+
+//------------------------------------------------------------------------------
+// Class to compare properties of reconstructed and MC particle tracks. 
+// 
+// Author: J.Otwinowski 04/02/2008 
+//------------------------------------------------------------------------------
+
+class AliComparisonRes;
+class AliComparisonEff;
+class AliComparisonDEdx;
+class AliComparisonDCA;
+class AliMagFMaps;
+class TList;
+
+#include "AliAnalysisTask.h"
+
+class AliComparisonTask : public AliAnalysisTask {
+ public:
+  AliComparisonTask(const char *name = "AliComparisonTask");
+  virtual ~AliComparisonTask();
+  
+  virtual void   ConnectInputData(Option_t *);
+  virtual void   CreateOutputObjects();
+  virtual void   Exec(Option_t *option);
+  virtual void   Terminate(Option_t *);
+
+  // Read TTree entry (event by event)
+  Bool_t  ReadEntry(Int_t evt);
+
+  // Set comparison objects
+  void SetAliComparisonRes(AliComparisonRes* comp) {fCompRes = comp;}
+  void SetAliComparisonEff(AliComparisonEff* comp) {fCompEff = comp;}
+  void SetAliComparisonDEdx(AliComparisonDEdx* comp) {fCompDEdx = comp;}
+  void SetAliComparisonDCA(AliComparisonDCA* comp) {fCompDCA = comp;}
+
+  void SetMagField(Int_t mag = 2) {fMagField = mag;}
+  void SetGeometry(char* geom = "/d/alice12/jacek/sim/v4-10-Release/pp/0/geometry.root")  {fGeom = geom;}
+
+ private:
+  TTree* fTree;                   //! input tree
+  AliMCInfo *fInfoMC;             //! AliMCInfo object
+  AliESDRecInfo *fInfoRC;         //! AliESDRecInfo object
+  AliComparisonRes* fCompRes;     // TPC resolution comparison object
+  AliComparisonEff* fCompEff;     // TPC efficiency comparison object
+  AliComparisonDEdx* fCompDEdx;   // TPC DEdx comparison object
+  AliComparisonDCA* fCompDCA;     // TPC DCA comparison object
+
+  TList* fOutput;                 //! list send on output slot 0
+  static Int_t evtNumber;         //! event number
+  Int_t  fMagField;               //! mag. field (0 - 0.2 T, 1 - 0.4 T, 2 - 0.5 T) 
+  AliMagFMaps *fMagFMap;          //! mag. field map 
+  const char *fGeom;              //! ROOT file with detector geometry
+
+  AliComparisonTask(const AliComparisonTask&); // not implemented
+  AliComparisonTask& operator=(const AliComparisonTask&); // not implemented
+  
+  ClassDef(AliComparisonTask, 1); // example of analysis
+};
+
+#endif
index dc30eda8666fae3a63d623f61494ee28be48c331..49ab539b9fbba1396427007d6cc30d6624b66622 100644 (file)
 //
 #pragma link C++ class AliRecInfoMaker+;
 #pragma link C++ class AliComparisonDraw+;
-#pragma link C++ class AliGenInfoTask+;
+
+#pragma link C++ class AliRecInfoCuts+;
+#pragma link C++ class AliMCInfoCuts+;
+#pragma link C++ class AliComparisonTask+;
+#pragma link C++ class AliComparisonRes+;
+#pragma link C++ class AliComparisonEff+;
+#pragma link C++ class AliComparisonDEdx+;
+#pragma link C++ class AliComparisonDCA+;
+
 
 #endif
index d92cd494b5324ddfc60252606dc77c07cc09b050..9ece35b00f7c47be1a560f7b40c49c2c4c16cc61 100644 (file)
@@ -9,10 +9,17 @@ SRCS:=        AliTreeDraw.cxx \
        AliESDRecKinkInfo.cxx \
        AliRecInfoMaker.cxx \
        AliComparisonDraw.cxx \
-       AliGenInfoTask.cxx
+       AliRecInfoCuts.cxx \
+       AliMCInfoCuts.cxx \
+       AliComparisonTask.cxx \
+       AliComparisonRes.cxx \
+       AliComparisonEff.cxx \
+       AliComparisonDEdx.cxx \
+       AliComparisonDCA.cxx
+
  
 HDRS:= $(SRCS:.cxx=.h)
 
 DHDR:= PWG1LinkDef.h
 
-EINCLUDE:= STEER TPC ITS TRD ANALYSIS
+EINCLUDE:= STEER TPC ITS TRD PWG0