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