From 6d1c79ca3b0dac358da13a40626ec133f10fa05e Mon Sep 17 00:00:00 2001 From: marian Date: Wed, 2 Apr 2008 21:07:45 +0000 Subject: [PATCH] Add new type oh histograms and analysis All resolution plots also in 1/pt binning (Jacek Otwinowski, Marian Ivanov) --- PWG1/AliComparisonDCA.cxx | 560 ++++++------ PWG1/AliComparisonDEdx.cxx | 573 ++++++------ PWG1/AliComparisonEff.cxx | 1682 ++++++++++++++++++------------------ PWG1/AliComparisonRes.cxx | 909 +++++++++++-------- PWG1/AliComparisonRes.h | 178 ++-- PWG1/AliComparisonTask.cxx | 447 +++++----- PWG1/AliMCInfoCuts.h | 178 ++-- 7 files changed, 2372 insertions(+), 2155 deletions(-) diff --git a/PWG1/AliComparisonDCA.cxx b/PWG1/AliComparisonDCA.cxx index d5f0ca83750..fe5cec011a1 100644 --- a/PWG1/AliComparisonDCA.cxx +++ b/PWG1/AliComparisonDCA.cxx @@ -1,279 +1,281 @@ -//------------------------------------------------------------------------------ -// 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 - -#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()GetMaxR() && TMath::Abs(infoMC->GetParticle().Vz())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()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(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(); -} +//------------------------------------------------------------------------------ +// 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 + +#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",40,-2,2, 10,0.3,3, 100,-1,1); + fD0TanSPtB1->SetXTitle("tan(#theta)"); + fD0TanSPtB1->SetYTitle("#sqrt{p_{t}(GeV/c)}"); + fD0TanSPtB1->SetZTitle("DCA_{xy}"); + + fD1TanSPtB1 = new TH3F("DCAzTanSPtB1","DCAzTanSPt",40,-2,2, 10,0.3,3, 100,-1,1); + fD1TanSPtB1->SetXTitle("tan(#theta)"); + fD1TanSPtB1->SetYTitle("#sqrt(p_{t}(GeV/c))"); + fD1TanSPtB1->SetZTitle("DCA_{z}"); + + fD0TanSPtL1 = new TH3F("DCAyTanSPtL1","DCAyTanSPt",40,-2,2, 10,0.3,3, 100,-1,1); + fD0TanSPtL1->SetXTitle("tan(#theta)"); + fD0TanSPtL1->SetYTitle("#sqrt{p_{t}(GeV/c)}"); + fD0TanSPtL1->SetZTitle("DCA_{xy}"); + + fD1TanSPtL1 = new TH3F("DCAzTanSPtL1","DCAzTanSPt",40,-2,2, 10,0.3,3, 100, -1,1); + fD1TanSPtL1->SetXTitle("tan(#theta)"); + fD1TanSPtL1->SetYTitle("#sqrt{p_{t}(GeV/c)}"); + fD1TanSPtL1->SetZTitle("DCA_{z}"); + + fD0TanSPtInTPC = new TH3F("DCAyTanSPtInTPC","DCAyTanSPt",40,-2,2, 10,0.3,3, 100,-1,1); + fD0TanSPtInTPC->SetXTitle("tan(#theta)"); + fD0TanSPtInTPC->SetYTitle("#sqrt{p_{t}(GeV/c)}"); + fD0TanSPtInTPC->SetZTitle("DCA_{xy}"); + + fD1TanSPtInTPC = new TH3F("DCAzTanSPtInTPC","DCAzTanSPt",40,-2,2, 10,0.3,3, 100, -1,1); + fD1TanSPtInTPC->SetXTitle("tan(#theta)"); + fD1TanSPtInTPC->SetYTitle("#sqrt{p_{t}(GeV/c)}"); + 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 + AliExternalTrackParam *track = 0; + Double_t kRadius = 3.0; // beam pipe radius + Double_t kMaxStep = 5.0; // max step + Double_t field = AliTracker::GetBz(); // nominal Bz field [kG] + 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 + + Float_t mcpt = infoMC->GetParticle().Pt(); + Float_t tantheta = TMath::Tan(infoMC->GetParticle().Theta()-TMath::Pi()*0.5); + Float_t spt = TMath::Sqrt(mcpt); + + // distance to Prim. vertex + const Double_t* dv = infoMC->GetVDist(); + + Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])GetMaxR() && TMath::Abs(dv[2])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()GetMinNClustersTPC()) return; + if (!infoRC->GetESDtrack()->GetConstrainedParam()) return; + + // 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] ); + + // calculate track parameters at vertex + if (infoRC->GetESDtrack()->GetTPCInnerParam()) + { + if ((track = new AliExternalTrackParam(*infoRC->GetESDtrack()->GetTPCInnerParam())) != 0 ) + { + Bool_t bStatus = AliTracker::PropagateTrackTo(track,kRadius,infoMC->GetMass(),kMaxStep,kTRUE); + Bool_t bDCAStatus = track->PropagateToDCA(fVertex,field,kMaxD,dca,cov); + + if(bStatus && bDCAStatus) { + 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(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("#sqrt{p_{t}(GeV/c)}"); + gr->GetYaxis()->SetTitle("#sigmaDCA (cm)"); + gr->GetHistogram()->Write("DCAResolSPTTan"); + + gPad->Clear(); + gr0->Draw("al*"); + + fp->Close(); +} diff --git a/PWG1/AliComparisonDEdx.cxx b/PWG1/AliComparisonDEdx.cxx index 50e8d295604..a60aae7a204 100644 --- a/PWG1/AliComparisonDEdx.cxx +++ b/PWG1/AliComparisonDEdx.cxx @@ -1,285 +1,288 @@ -//------------------------------------------------------------------------------ -// 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 - -#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()GetMaxR() && TMath::Abs(infoMC->GetParticle().Vz())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()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()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(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(""); - 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(""); - 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(); -} +//------------------------------------------------------------------------------ +// 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 + +#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); + Float_t mprim = infoMC->GetPrim(); + + // distance to Prim. vertex + const Double_t* dv = infoMC->GetVDist(); + + Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])GetMaxR() && TMath::Abs(dv[2])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()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()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(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(""); + 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(""); + 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/AliComparisonEff.cxx b/PWG1/AliComparisonEff.cxx index 70e1f30412b..46278ffb55c 100644 --- a/PWG1/AliComparisonEff.cxx +++ b/PWG1/AliComparisonEff.cxx @@ -1,841 +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 - -#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()GetMaxR() && TMath::Abs(infoMC->GetParticle().Vz())GetMaxVz(); - Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])GetMaxR() && TMath::Abs(dv[2])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)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 " <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)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(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(); -} +//------------------------------------------------------------------------------ +// 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 + +#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 = AliTracker::GetBz(); // nominal Bz field [kG] + 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 + + // distance to Prim. vertex + const Double_t* dv = infoMC->GetVDist(); + + Float_t mcpt = infoMC->GetParticle().Pt(); + Float_t tantheta = TMath::Tan(infoMC->GetParticle().Theta()-TMath::Pi()*0.5); + Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])GetMaxR() && TMath::Abs(dv[2])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 + // Pdg convension is different for hadrons and leptons + // (e.g. K+/K- = 321/-321; e+/e- = -11/11 ) + Double_t pid = -1; + if( TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetEM() ) pid = 0; + if( TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetMuM() ) 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)GetMaxAbsTanTheta()) + { + if (infoRC->GetESDtrack()->GetTPCInnerParam()) + { + if ((track = new AliExternalTrackParam(*infoRC->GetESDtrack()->GetTPCInnerParam())) != 0 ) + { + + Bool_t bStatus = AliTracker::PropagateTrackTo(track,kRadius,infoMC->GetMass(),kMaxStep,kTRUE); + Bool_t bDCAStatus = track->PropagateToDCA(fVertex,field,kMaxD,dca,cov); + + if(bStatus && bDCAStatus) { + // + 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 " <Fill(mcpt); + 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)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(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/AliComparisonRes.cxx b/PWG1/AliComparisonRes.cxx index 70593fe8a93..8e715135bdb 100644 --- a/PWG1/AliComparisonRes.cxx +++ b/PWG1/AliComparisonRes.cxx @@ -1,358 +1,551 @@ -//------------------------------------------------------------------------------ -// 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 - -#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()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(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()GetMaxR() && TMath::Abs(infoMC->GetParticle().Vz())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()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(); -} - - - +//------------------------------------------------------------------------------ +// 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 + +#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 "AliESDVertex.h" +#include "AliRecInfoCuts.h" +#include "AliMCInfoCuts.h" +#include "AliLog.h" +#include "AliTracker.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 + + // + // Parametrisation histograms + // + fC1Pt2Resol1PtTPC(0), + fC1Pt2Resol1PtTPCITS(0), + fCYResol1PtTPC(0), + fCYResol1PtTPCITS(0), + fCZResol1PtTPC(0), + fCZResol1PtTPCITS(0), + fCPhiResol1PtTPC(0), + fCPhiResol1PtTPCITS(0), + fCThetaResol1PtTPC(0), + fCThetaResol1PtTPCITS(0), + + // vertex + 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); +} + +//_____________________________________________________________________________ +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; + + // Parametrisation histograms + if(fC1Pt2Resol1PtTPC) delete fC1Pt2Resol1PtTPC; fC1Pt2Resol1PtTPC=0; + if(fC1Pt2Resol1PtTPCITS) delete fC1Pt2Resol1PtTPCITS; fC1Pt2Resol1PtTPCITS=0; + if(fCYResol1PtTPC) delete fCYResol1PtTPC; fCYResol1PtTPC=0; + if(fCYResol1PtTPCITS) delete fCYResol1PtTPCITS; fCYResol1PtTPCITS=0; + if(fCZResol1PtTPC) delete fCZResol1PtTPC; fCZResol1PtTPC=0; + if(fCZResol1PtTPCITS) delete fCZResol1PtTPCITS; fCZResol1PtTPCITS=0; + if(fCPhiResol1PtTPC) delete fCPhiResol1PtTPC; fCPhiResol1PtTPC=0; + if(fCPhiResol1PtTPCITS) delete fCPhiResol1PtTPCITS; fCPhiResol1PtTPCITS=0; + if(fCThetaResol1PtTPC) delete fCThetaResol1PtTPC; fCThetaResol1PtTPC=0; + if(fCThetaResol1PtTPCITS) delete fCThetaResol1PtTPCITS; fCThetaResol1PtTPCITS=0; + + if(fVertex) delete fVertex; fVertex=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 pull",10, 0.1,3,200,-6,6); + fPtPullLPT->SetXTitle("p_{t}"); + fPtPullLPT->SetYTitle("#Deltap_{t}/#Sigma"); + + fPtPullHPT = new TH2F("Pt_pull_hpt","pt pull",10,2,100,200,-6,6); + fPtPullHPT->SetXTitle("p_{t}"); + fPtPullHPT->SetYTitle("#Deltap_{t}/#Sigma"); + + // Parametrisation histograms + fC1Pt2Resol1PtTPC = new TH2F("fC1Pt2Resol1PtTPC","(1/mcpt-1/pt)/(1+1/mcpt)^2 vs 1/pt)",100,0,10,200,-0.020,0.020); + fC1Pt2Resol1PtTPC->SetXTitle("1/mcp_{t}"); + fC1Pt2Resol1PtTPC->SetYTitle("(1/mcp_{t}-1/p_{t})/(1+1/mcp_{t})^2)"); + + fC1Pt2Resol1PtTPCITS = new TH2F("fC1Pt2Resol1PtTPCITS","(1/mcpt-1/pt)/(1+1/mcpt)^2 vs 1/pt)",100,0,10,200,-0.020,0.020); + fC1Pt2Resol1PtTPCITS->SetXTitle("1/mcp_{t}"); + fC1Pt2Resol1PtTPCITS->SetYTitle("(1/mcp_{t}-1/p_{t})/(1+1/mcp_{t})^2)"); + + fCYResol1PtTPC = new TH2F("fCYResol1PtTPC","fCYResol1PtTPC",100, 0,10,200,-0.010,0.010); + fCYResol1PtTPC->SetXTitle("1/mcpt"); + fCYResol1PtTPC->SetYTitle("#DeltaY"); + + fCYResol1PtTPCITS = new TH2F("fCYResol1PtTPCITS","fCYResol1PtTPCITS",100, 0,10,200,-0.010,0.010); + fCYResol1PtTPCITS->SetXTitle("1/mcpt"); + fCYResol1PtTPCITS->SetYTitle("#DeltaY"); + + fCZResol1PtTPC = new TH2F("fCZResol1PtTPC","fCZResol1PtTPC",50, 0,10,200,-0.020,0.020); + fCZResol1PtTPC->SetXTitle("1/mcpt"); + fCZResol1PtTPC->SetYTitle("#DeltaZ"); + + fCZResol1PtTPCITS = new TH2F("fCZResol1PtTPCITS","fCZResol1PtTPCITS",50, 0,10,200,-0.020,0.020); + fCZResol1PtTPCITS->SetXTitle("1/mcpt"); + fCZResol1PtTPCITS->SetYTitle("#DeltaZ"); + + fCPhiResol1PtTPC = new TH2F("fCPhiResol1PtTPC","fCPhiResol1PtTPC",50, 0,10,200,-0.005,0.005); + fCPhiResol1PtTPC->SetXTitle("1/mcpt"); + fCPhiResol1PtTPC->SetYTitle("#Delta#phi"); + + fCPhiResol1PtTPCITS = new TH2F("fCPhiResol1PtTPCITS","fCPhiResol1PtTPCITS",50, 0,10,200,-0.005,0.005); + fCPhiResol1PtTPCITS->SetXTitle("1/mcpt"); + fCPhiResol1PtTPCITS->SetYTitle("#Delta#phi"); + + fCThetaResol1PtTPC = new TH2F("fCThetaResol1PtTPC","fCThetaResol1PtTPC",50, 0,10,200,-0.005,0.005); + fCThetaResol1PtTPC->SetXTitle("1/mcpt"); + fCThetaResol1PtTPC->SetYTitle("#Delta#theta"); + + fCThetaResol1PtTPCITS = new TH2F("fCThetaResol1PtTPCITS","fCThetaResol1PtTPCITS",50, 0,10,200,-0.005,0.005); + fCThetaResol1PtTPCITS->SetXTitle("1/mcpt"); + fCThetaResol1PtTPCITS->SetYTitle("#Delta#theta"); +} + +//_____________________________________________________________________________ +void AliComparisonRes::InitCuts() +{ + // 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(); + + // distance to Prim. vertex + const Double_t* dv = infoMC->GetVDist(); + + Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])GetMaxR() && TMath::Abs(dv[2])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()GetMinNClustersTPC()) return; + + // 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] ); + + 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); +} + +//_____________________________________________________________________________ +void AliComparisonRes::ProcessConstrained(AliMCInfo* infoMC, AliESDRecInfo *infoRC) +{ + // Fill resolution comparison information (constarained parameters) + // + AliExternalTrackParam *track = 0; + Double_t kRadius = 3.0; // beam pipe radius + Double_t kMaxStep = 5.0; // max step + Double_t field = AliTracker::GetBz(); // nominal Bz field [kG] + 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 + + Float_t mcpt = infoMC->GetParticle().Pt(); + Float_t tantheta = TMath::Tan(infoMC->GetParticle().Theta()-TMath::Pi()*0.5); + + // distance to Prim. vertex + const Double_t* dv = infoMC->GetVDist(); + Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])GetMaxR() && TMath::Abs(dv[2])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()GetMinNClustersTPC()) return; + if (!infoRC->GetESDtrack()->GetConstrainedParam()) return; + + // 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] ); + + // 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()); + + + //Float_t delta1Pt2 = (1/mcpt-infoRC->GetESDtrack()->OneOverPt())/TMath::Power(1+1/mcpt,2)); + Float_t delta1Pt2 = (1/mcpt-cparam->OneOverPt())/TMath::Power(1+1/mcpt,2); + + Float_t deltaY1Pt = (infoMC->GetParticle().Vy()-cparam->GetY()) / (0.2+1/mcpt); + Float_t deltaZ1Pt = (infoMC->GetParticle().Vz()-cparam->GetZ()) / (0.2+1/mcpt); + Float_t deltaPhi1Pt = deltaPhi / (0.1+1/mcpt); + Float_t deltaTheta1Pt = deltaTan / (0.1+1/mcpt); + + // calculate track parameters at vertex + if (infoRC->GetESDtrack()->GetTPCInnerParam()) + { + if ((track = new AliExternalTrackParam(*infoRC->GetESDtrack()->GetTPCInnerParam())) != 0 ) + { + Bool_t bStatus = AliTracker::PropagateTrackTo(track,kRadius,infoMC->GetMass(),kMaxStep,kTRUE); + Bool_t bDCAStatus = track->PropagateToDCA(fVertex,field,kMaxD,dca,cov); + + // Fill parametrisation histograms (only TPC track) + if(bStatus && bDCAStatus) + { + fC1Pt2Resol1PtTPC->Fill(1/mcpt,delta1Pt2); + fCYResol1PtTPC->Fill(1/mcpt,deltaY1Pt); + fCZResol1PtTPC->Fill(1/mcpt,deltaZ1Pt); + fCPhiResol1PtTPC->Fill(1/mcpt,deltaPhi1Pt); + fCThetaResol1PtTPC->Fill(1/mcpt,deltaTheta1Pt); + } + delete track; + } + } + + // TPC and ITS (nb. of clusters >2) in the system + if(infoRC->GetStatus(1)==3 && infoRC->GetESDtrack()->GetITSclusters(clusterITS)>2) + { + fC1Pt2Resol1PtTPCITS->Fill(1/mcpt,delta1Pt2); + fCYResol1PtTPCITS->Fill(1/mcpt,deltaY1Pt); + fCZResol1PtTPCITS->Fill(1/mcpt,deltaZ1Pt); + fCPhiResol1PtTPCITS->Fill(1/mcpt,deltaPhi1Pt); + fCThetaResol1PtTPCITS->Fill(1/mcpt,deltaTheta1Pt); + } + + // 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"); + // + hiss = comp->MakeResol(comp->fC1Pt2Resol1PtTPC,1,0); + hiss->SetXTitle("1/mcp_{t}"); + hiss->SetYTitle("1/mcp_{t}-1/p_{t}/(1+1/p_{t})^2"); + hiss->Draw(); + hiss->Write("Pt2Resol1PtTPC"); + fC1Pt2Resol1PtTPC->Write(); + + hiss = comp->MakeResol(comp->fC1Pt2Resol1PtTPCITS,1,0); + hiss->SetXTitle("1/mcp_{t}"); + hiss->SetYTitle("1/mcp_{t}-1/p_{t}/(1+1/p_{t})^2"); + hiss->Draw(); + hiss->Write("Pt2Resol1PtTPCITS"); + fC1Pt2Resol1PtTPCITS->Write(); + // + hiss = comp->MakeResol(comp->fCYResol1PtTPC,1,0); + hiss->SetXTitle("1/mcp_{t}"); + hiss->SetYTitle("(mcy-y)/(0.2+1/mcp_{t})"); + hiss->Draw(); + hiss->Write("CYResol1PtTPC"); + fCYResol1PtTPC->Write(); + + hiss = comp->MakeResol(comp->fCYResol1PtTPCITS,1,0); + hiss->SetXTitle("1/mcp_{t}"); + hiss->SetYTitle("(mcy-y)/(0.2+1/mcp_{t})"); + hiss->Draw(); + hiss->Write("CYResol1PtTPCITS"); + fCYResol1PtTPCITS->Write(); + // + hiss = comp->MakeResol(comp->fCZResol1PtTPC,1,0); + hiss->SetXTitle("1/mcp_{t}"); + hiss->SetYTitle("(mcz-z)/(0.2+1/mcp_{t})"); + hiss->Draw(); + hiss->Write("CZResol1PtTPC"); + fCZResol1PtTPC->Write(); + + hiss = comp->MakeResol(comp->fCZResol1PtTPCITS,1,0); + hiss->SetXTitle("1/mcp_{t}"); + hiss->SetYTitle("(mcz-z)/(0.2+1/mcp_{t})"); + hiss->Draw(); + hiss->Write("CZResol1PtTPCITS"); + fCZResol1PtTPCITS->Write(); + // + hiss = comp->MakeResol(comp->fCPhiResol1PtTPC,1,0); + hiss->SetXTitle("1/mcp_{t}"); + hiss->SetYTitle("(mc#phi-#phi)/(0.1+1/mcp_{t})"); + hiss->Draw(); + hiss->Write("CPhiResol1PtTPC"); + fCPhiResol1PtTPC->Write(); + + hiss = comp->MakeResol(comp->fCPhiResol1PtTPCITS,1,0); + hiss->SetXTitle("1/mcp_{t}"); + hiss->SetYTitle("(mc#phi-#phi)/(0.1+1/mcp_{t})"); + hiss->Draw(); + hiss->Write("CPhiResol1PtTPCITS"); + fCPhiResol1PtTPCITS->Write(); + // + hiss = comp->MakeResol(comp->fCThetaResol1PtTPC,1,0); + hiss->SetXTitle("1/mcp_{t}"); + hiss->SetYTitle("(mc#theta-#theta)/(0.1+1/mcp_{t})"); + hiss->Draw(); + hiss->Write("CThetaResol1PtTPC"); + fCThetaResol1PtTPC->Write(); + + hiss = comp->MakeResol(comp->fCThetaResol1PtTPCITS,1,0); + hiss->SetXTitle("1/mcp_{t}"); + hiss->SetYTitle("(mc#theta-#theta)/(0.1+1/mcp_{t})"); + hiss->Draw(); + hiss->Write("CThetaResol1PtTPCITS"); + fCThetaResol1PtTPCITS->Write(); + + fp->Close(); +} + +//_____________________________________________________________________________ +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(obj); + if (entry == 0) continue; + + 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); + + // Histograms for 1/pt parameterisation + fC1Pt2Resol1PtTPC->Add(entry->fC1Pt2Resol1PtTPC); + fCYResol1PtTPC->Add(entry->fCYResol1PtTPC); + fCZResol1PtTPC->Add(entry->fCZResol1PtTPC); + fCPhiResol1PtTPC->Add(entry->fCPhiResol1PtTPC); + fCThetaResol1PtTPC->Add(entry->fCThetaResol1PtTPC); + + fC1Pt2Resol1PtTPCITS->Add(entry->fC1Pt2Resol1PtTPCITS); + fCYResol1PtTPCITS->Add(entry->fCYResol1PtTPCITS); + fCZResol1PtTPCITS->Add(entry->fCZResol1PtTPCITS); + fCPhiResol1PtTPCITS->Add(entry->fCPhiResol1PtTPCITS); + fCThetaResol1PtTPCITS->Add(entry->fCThetaResol1PtTPCITS); + + count++; + } + +return count; +} diff --git a/PWG1/AliComparisonRes.h b/PWG1/AliComparisonRes.h index df1e3fcb96e..3c3881be35b 100644 --- a/PWG1/AliComparisonRes.h +++ b/PWG1/AliComparisonRes.h @@ -1,80 +1,98 @@ -#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 +#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; +class AliESDVertex; + +#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 + + // + // Histograms for track resolution parameterisation + // + TH2F* fC1Pt2Resol1PtTPC; //-> (1/mcpt-1/pt)/(1+1/pt)^2 vs 1/pt (TPC) + TH2F* fC1Pt2Resol1PtTPCITS; //-> (1/mcpt-1/pt)/(1+1/pt)^2 vs 1/pt (TPC+ITS) + TH2F* fCYResol1PtTPC; //-> (mcy-y)/(0.2+1/pt) vs 1/pt (TPC) + TH2F* fCYResol1PtTPCITS; //-> (mcy-y)/(0.2+1/pt) vs 1/pt (TPC + ITS) + TH2F* fCZResol1PtTPC; //-> (mcz-z)/(0.2+1/pt) vs 1/pt (TPC) + TH2F* fCZResol1PtTPCITS; //-> (mcz-z)/(0.2+1/pt) vs 1/pt (TPC+ITS) + TH2F* fCPhiResol1PtTPC; //-> (mcphi-phi)/(0.1+1/pt) vs 1/pt (TPC) + TH2F* fCPhiResol1PtTPCITS; //-> (mcphi-phi)/(0.1+1/pt) vs 1/pt (TPC+ITS) + TH2F* fCThetaResol1PtTPC; //-> (mctheta-theta)/(0.1+1/pt) vs 1/pt (TPC) + TH2F* fCThetaResol1PtTPCITS; //-> (mctheta-theta)/(0.1+1/pt) vs 1/pt (TPC+ITS) + + AliESDVertex *fVertex; //! + + // 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 index def4f22e3f9..1d329f34695 100644 --- a/PWG1/AliComparisonTask.cxx +++ b/PWG1/AliComparisonTask.cxx @@ -1,223 +1,224 @@ -//------------------------------------------------------------------------------ -// 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 (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 (GetOutputData(0)); - if (!fOutput) { - Printf("ERROR: fOutput not available"); - return; - } - - fCompRes = dynamic_cast (fOutput->FindObject("AliComparisonRes")); - if (!fCompRes) { - Printf("WARNING: AliComparisonRes not available"); - return; - } - - fCompEff = dynamic_cast (fOutput->FindObject("AliComparisonEff")); - if (!fCompEff) { - Printf("WARNING: AliComparisonEff not available"); - return; - } - - fCompDEdx = dynamic_cast (fOutput->FindObject("AliComparisonDEdx")); - if (!fCompDEdx) { - Printf("WARNING: AliComparisonDEdx not available"); - return; - } - - fCompDCA = dynamic_cast (fOutput->FindObject("AliComparisonDCA")); - if (!fCompDCA) { - Printf("WARNING: AliComparisonDCA not available"); - return; - } - - fOutput->Write(); - out->Close(); -} +//------------------------------------------------------------------------------ +// 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 (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 (GetOutputData(0)); + if (!fOutput) { + Printf("ERROR: fOutput not available"); + return; + } + + fCompRes = dynamic_cast (fOutput->FindObject("AliComparisonRes")); + if (!fCompRes) { + Printf("WARNING: AliComparisonRes not available"); + return; + } + + fCompEff = dynamic_cast (fOutput->FindObject("AliComparisonEff")); + if (!fCompEff) { + Printf("WARNING: AliComparisonEff not available"); + return; + } + + fCompDEdx = dynamic_cast (fOutput->FindObject("AliComparisonDEdx")); + if (!fCompDEdx) { + Printf("WARNING: AliComparisonDEdx not available"); + return; + } + + fCompDCA = dynamic_cast (fOutput->FindObject("AliComparisonDCA")); + if (!fCompDCA) { + Printf("WARNING: AliComparisonDCA not available"); + return; + } + + fOutput->Write(); + out->Close(); +} diff --git a/PWG1/AliMCInfoCuts.h b/PWG1/AliMCInfoCuts.h index 55eb9bcd39c..4806166cff0 100644 --- a/PWG1/AliMCInfoCuts.h +++ b/PWG1/AliMCInfoCuts.h @@ -1,89 +1,89 @@ -#ifndef ALIMCINFOCUTS_H -#define ALIMCINFOCUTS_H - -//------------------------------------------------------------------------------ -// Class to keep selection cuts for MC tracks. -// -// Author: J.Otwinowski 04/02/2008 -//------------------------------------------------------------------------------ - -#include "AliAnalysisCuts.h" - -class TArrayI; - -class AliMCInfoCuts : public AliAnalysisCuts -{ -public: - AliMCInfoCuts(const Char_t* name ="AliMCInfoCuts", const Char_t *title =""); - virtual ~AliMCInfoCuts(); - - // setters - void SetMinRowsWithDigits(const Int_t min=0) {fMinRowsWithDigits = min;} - void SetMaxR(const Float_t max=1e99) {fMaxR = max;} - void SetMaxVz(const Float_t max=1e99) {fMaxVz = max;} - void SetRangeTPCSignal(const Float_t min=0, const Float_t max=1e99) {fMinTPCSignal = min; fMaxTPCSignal = max;} - - // getters - Int_t GetMinRowsWithDigits() const {return fMinRowsWithDigits;} - Float_t GetMaxR() const {return fMaxR;} - Float_t GetMaxVz() const {return fMaxVz;} - Float_t GetMinTPCSignal() const {return fMinTPCSignal;} - Float_t GetMaxTPCSignal() const {return fMaxTPCSignal;} - - Float_t GetEP() const {return ep;} - Float_t GetEM() const {return em;} - Float_t GetMuP() const {return mup;} - Float_t GetMuM() const {return mum;} - Float_t GetPiP() const {return pip;} - Float_t GetPiM() const {return pim;} - Float_t GetKP() const {return kp;} - Float_t GetKM() const {return km;} - Float_t GetProt() const {return prot;} - Float_t GetProtBar() const {return protbar;} - - // cuts init function - void Init(); - - // check MC tracks - Bool_t IsSelected(TObject *) {return kTRUE;} - - // add particle to array - void AddPdgParticle(Int_t idx=-1, Int_t pdgcode=0) const; - - // check particle in array - Bool_t IsPdgParticle(Int_t pdgcode=0) const; - - // Merge output objects (needed by PROOF) - virtual Long64_t Merge(TCollection* list); - -private: - Int_t fMinRowsWithDigits; // min. number of TPC raws with digits - Float_t fMaxR; // max. R distance from MC vertex - Float_t fMaxVz; // max. Z distance from MC vertex - Float_t fMinTPCSignal; // min. TPC Signal calculated from Bethe Bloch formula - Float_t fMaxTPCSignal; // max. TPC Signal calculated from Bethe Bloch formula - - TArrayI* aTrackParticles; // array of tracked particles - - // PDG tracked particles (later added to aTrackParticles) - enum enumData { - kNParticles = 10, // number of particles below - ep = 11, - em = -11, - mup = 13, - mum = -13, - pip = 211, - pim = -211, - kp = 321, - km = -321, - prot = 2212, - protbar = -2212 - }; - - AliMCInfoCuts(const AliMCInfoCuts&); // not implemented - AliMCInfoCuts& operator=(const AliMCInfoCuts&); // not implemented - - ClassDef(AliMCInfoCuts, 1) -}; - -#endif // ALIMCINFOCUTS_H +#ifndef ALIMCINFOCUTS_H +#define ALIMCINFOCUTS_H + +//------------------------------------------------------------------------------ +// Class to keep selection cuts for MC tracks. +// +// Author: J.Otwinowski 04/02/2008 +//------------------------------------------------------------------------------ + +#include "AliAnalysisCuts.h" + +class TArrayI; + +class AliMCInfoCuts : public AliAnalysisCuts +{ +public: + AliMCInfoCuts(const Char_t* name ="AliMCInfoCuts", const Char_t *title =""); + virtual ~AliMCInfoCuts(); + + // setters + void SetMinRowsWithDigits(const Int_t min=0) {fMinRowsWithDigits = min;} + void SetMaxR(const Float_t max=1e99) {fMaxR = max;} + void SetMaxVz(const Float_t max=1e99) {fMaxVz = max;} + void SetRangeTPCSignal(const Float_t min=0, const Float_t max=1e99) {fMinTPCSignal = min; fMaxTPCSignal = max;} + + // getters + Int_t GetMinRowsWithDigits() const {return fMinRowsWithDigits;} + Float_t GetMaxR() const {return fMaxR;} + Float_t GetMaxVz() const {return fMaxVz;} + Float_t GetMinTPCSignal() const {return fMinTPCSignal;} + Float_t GetMaxTPCSignal() const {return fMaxTPCSignal;} + + Float_t GetEP() const {return ep;} + Float_t GetEM() const {return em;} + Float_t GetMuP() const {return mup;} + Float_t GetMuM() const {return mum;} + Float_t GetPiP() const {return pip;} + Float_t GetPiM() const {return pim;} + Float_t GetKP() const {return kp;} + Float_t GetKM() const {return km;} + Float_t GetProt() const {return prot;} + Float_t GetProtBar() const {return protbar;} + + // cuts init function + void Init(); + + // check MC tracks + Bool_t IsSelected(TObject *) {return kTRUE;} + + // add particle to array + void AddPdgParticle(Int_t idx=-1, Int_t pdgcode=0) const; + + // check particle in array + Bool_t IsPdgParticle(Int_t pdgcode=0) const; + + // Merge output objects (needed by PROOF) + virtual Long64_t Merge(TCollection* list); + +private: + Int_t fMinRowsWithDigits; // min. number of TPC raws with digits + Float_t fMaxR; // max. R distance from MC vertex + Float_t fMaxVz; // max. Z distance from MC vertex + Float_t fMinTPCSignal; // min. TPC Signal calculated from Bethe Bloch formula + Float_t fMaxTPCSignal; // max. TPC Signal calculated from Bethe Bloch formula + + TArrayI* aTrackParticles; // array of tracked particles + + // PDG tracked particles (later added to aTrackParticles) + enum enumData { + kNParticles = 10, // number of particles below + ep = -11, + em = 11, + mup = -13, + mum = 13, + pip = 211, + pim = -211, + kp = 321, + km = -321, + prot = 2212, + protbar = -2212 + }; + + AliMCInfoCuts(const AliMCInfoCuts&); // not implemented + AliMCInfoCuts& operator=(const AliMCInfoCuts&); // not implemented + + ClassDef(AliMCInfoCuts, 1) +}; + +#endif // ALIMCINFOCUTS_H -- 2.43.5