-//------------------------------------------------------------------------------\r
-// Implementation of AliComparisonRes class. It keeps information from \r
-// comparison of reconstructed and MC particle tracks. In addtion, \r
-// it keeps selection cuts used during comparison. The comparison \r
-// information is stored in the ROOT histograms. Analysis of these \r
-// histograms can be done by using Analyse() class function. The result of \r
-// the analysis (histograms) are stored in the output picture_res.root file.\r
-// \r
-// Author: J.Otwinowski 04/02/2008 \r
-//------------------------------------------------------------------------------\r
-\r
-/*\r
- //after running analysis, read the file, and get component\r
- gSystem->Load("libPWG1.so");\r
- TFile f("Output.root");\r
- AliComparisonRes * comp = (AliComparisonRes*)f.Get("AliComparisonRes");\r
-\r
- // analyse comparison data (output stored in pictures_res.root)\r
- comp->Analyse();\r
- \r
- // paramtetrisation of the TPC track length (for information only) \r
- TF1 fl("fl","((min(250./(abs(x+0.000001)),250)-90))",0,2); // TPC track length function\r
- TF1 fl2("fl2","[0]/((min(250./(abs(x+0.000001)),250)-90))^[1]",0,2);\r
- fl2.SetParameter(1,1);\r
- fl2.SetParameter(0,1);\r
-*/\r
-\r
-#include <iostream>\r
-\r
-#include "TFile.h"\r
-#include "TCint.h"\r
-#include "TH3F.h"\r
-#include "TH2F.h"\r
-#include "TF1.h"\r
-#include "TProfile.h"\r
-#include "TProfile2D.h"\r
-#include "TGraph2D.h"\r
-#include "TCanvas.h"\r
-#include "TGraph.h"\r
-\r
-#include "AliESDEvent.h" \r
-#include "AliESD.h"\r
-#include "AliESDfriend.h"\r
-#include "AliESDfriendTrack.h"\r
-#include "AliESDVertex.h"\r
-#include "AliRecInfoCuts.h" \r
-#include "AliMCInfoCuts.h" \r
-#include "AliLog.h" \r
-#include "AliTracker.h" \r
-\r
-#include "AliMathBase.h"\r
-#include "AliTreeDraw.h" \r
-\r
-#include "AliMCInfo.h" \r
-#include "AliESDRecInfo.h" \r
-#include "AliComparisonRes.h" \r
-\r
-using namespace std;\r
-\r
-ClassImp(AliComparisonRes)\r
-\r
-//_____________________________________________________________________________\r
-AliComparisonRes::AliComparisonRes():\r
- TNamed("AliComparisonRes","AliComparisonRes"),\r
-\r
- // Resolution \r
- fPtResolLPT(0), // pt resolution - low pt\r
- fPtResolHPT(0), // pt resolution - high pt \r
- fPtPullLPT(0), // pt resolution - low pt\r
- fPtPullHPT(0), // pt resolution - high pt \r
- //\r
- // Resolution constrained param\r
- //\r
- fCPhiResolTan(0), // angular resolution - constrained\r
- fCTanResolTan(0), // angular resolution - constrained\r
- fCPtResolTan(0), // pt resolution - constrained\r
- fCPhiPullTan(0), // angular resolution - constrained\r
- fCTanPullTan(0), // angular resolution - constrained\r
- fCPtPullTan(0), // pt resolution - constrained\r
-\r
- //\r
- // Parametrisation histograms\r
- //\r
-\r
- f1Pt2Resol1PtTPC(0),\r
- f1Pt2Resol1PtTPCITS(0),\r
- fYResol1PtTPC(0),\r
- fYResol1PtTPCITS(0),\r
- fZResol1PtTPC(0),\r
- fZResol1PtTPCITS(0),\r
- fPhiResol1PtTPC(0),\r
- fPhiResol1PtTPCITS(0),\r
- fThetaResol1PtTPC(0),\r
- fThetaResol1PtTPCITS(0),\r
-\r
- // constrained\r
- fC1Pt2Resol1PtTPC(0),\r
- fC1Pt2Resol1PtTPCITS(0),\r
- fCYResol1PtTPC(0),\r
- fCYResol1PtTPCITS(0),\r
- fCZResol1PtTPC(0),\r
- fCZResol1PtTPCITS(0),\r
- fCPhiResol1PtTPC(0),\r
- fCPhiResol1PtTPCITS(0),\r
- fCThetaResol1PtTPC(0),\r
- fCThetaResol1PtTPCITS(0),\r
-\r
- // vertex\r
- fVertex(0),\r
- \r
- // Cuts \r
- fCutsRC(0), \r
- fCutsMC(0) \r
-{\r
- InitHisto();\r
- InitCuts();\r
- \r
- // vertex (0,0,0)\r
- fVertex = new AliESDVertex();\r
- fVertex->SetXv(0.0);\r
- fVertex->SetYv(0.0);\r
- fVertex->SetZv(0.0);\r
-}\r
-\r
-//_____________________________________________________________________________\r
-AliComparisonRes::~AliComparisonRes(){\r
- \r
- // Resolution histograms\r
- if(fPtResolLPT) delete fPtResolLPT; fPtResolLPT=0; \r
- if(fPtResolHPT) delete fPtResolHPT; fPtResolHPT=0; \r
- if(fPtPullLPT) delete fPtPullLPT; fPtPullLPT=0; \r
- if(fPtPullHPT) delete fPtPullHPT; fPtPullHPT=0; \r
-\r
- // Resolution histograms (constrained param)\r
- if(fCPhiResolTan) delete fCPhiResolTan; fCPhiResolTan=0;\r
- if(fCTanResolTan) delete fCTanResolTan; fCTanResolTan=0;\r
- if(fCPtResolTan) delete fCPtResolTan; fCPtResolTan=0; \r
- if(fCPhiPullTan) delete fCPhiPullTan; fCPhiPullTan=0;\r
- if(fCTanPullTan) delete fCTanPullTan; fCTanPullTan=0;\r
- if(fCPtPullTan) delete fCPtPullTan; fCPtPullTan=0;\r
-\r
- // Parametrisation histograms\r
- // \r
- if(f1Pt2Resol1PtTPC) delete f1Pt2Resol1PtTPC; f1Pt2Resol1PtTPC=0;\r
- if(f1Pt2Resol1PtTPCITS) delete f1Pt2Resol1PtTPCITS; f1Pt2Resol1PtTPCITS=0;\r
- if(fYResol1PtTPC) delete fYResol1PtTPC; fYResol1PtTPC=0;\r
- if(fYResol1PtTPCITS) delete fYResol1PtTPCITS; fYResol1PtTPCITS=0;\r
- if(fZResol1PtTPC) delete fZResol1PtTPC; fZResol1PtTPC=0;\r
- if(fZResol1PtTPCITS) delete fZResol1PtTPCITS; fZResol1PtTPCITS=0;\r
- if(fPhiResol1PtTPC) delete fPhiResol1PtTPC; fPhiResol1PtTPC=0;\r
- if(fPhiResol1PtTPCITS) delete fPhiResol1PtTPCITS; fPhiResol1PtTPCITS=0;\r
- if(fThetaResol1PtTPC) delete fThetaResol1PtTPC; fThetaResol1PtTPC=0;\r
- if(fThetaResol1PtTPCITS) delete fThetaResol1PtTPCITS; fThetaResol1PtTPCITS=0;\r
-\r
- // constrained\r
- if(fC1Pt2Resol1PtTPC) delete fC1Pt2Resol1PtTPC; fC1Pt2Resol1PtTPC=0;\r
- if(fC1Pt2Resol1PtTPCITS) delete fC1Pt2Resol1PtTPCITS; fC1Pt2Resol1PtTPCITS=0;\r
- if(fCYResol1PtTPC) delete fCYResol1PtTPC; fCYResol1PtTPC=0;\r
- if(fCYResol1PtTPCITS) delete fCYResol1PtTPCITS; fCYResol1PtTPCITS=0;\r
- if(fCZResol1PtTPC) delete fCZResol1PtTPC; fCZResol1PtTPC=0;\r
- if(fCZResol1PtTPCITS) delete fCZResol1PtTPCITS; fCZResol1PtTPCITS=0;\r
- if(fCPhiResol1PtTPC) delete fCPhiResol1PtTPC; fCPhiResol1PtTPC=0;\r
- if(fCPhiResol1PtTPCITS) delete fCPhiResol1PtTPCITS; fCPhiResol1PtTPCITS=0;\r
- if(fCThetaResol1PtTPC) delete fCThetaResol1PtTPC; fCThetaResol1PtTPC=0;\r
- if(fCThetaResol1PtTPCITS) delete fCThetaResol1PtTPCITS; fCThetaResol1PtTPCITS=0;\r
-\r
- if(fVertex) delete fVertex; fVertex=0;\r
-\r
-}\r
-\r
-//_____________________________________________________________________________\r
-void AliComparisonRes::InitHisto(){\r
-\r
- // Init histograms\r
- fCPhiResolTan = new TH2F("CPhiResolTan","CPhiResolTan",50, -2,2,200,-0.025,0.025); \r
- fCPhiResolTan->SetXTitle("tan(#theta)");\r
- fCPhiResolTan->SetYTitle("#Delta#phi");\r
-\r
- fCTanResolTan = new TH2F("CTanResolTan","CTanResolTan",50, -2,2,200,-0.025,0.025);\r
- fCTanResolTan->SetXTitle("tan(#theta)");\r
- fCTanResolTan->SetYTitle("#Delta#theta");\r
-\r
- fCPtResolTan=new TH2F("CPtResol","CPtResol",50, -2,2,200,-0.2,0.2); \r
- fCPtResolTan->SetXTitle("Tan(#theta)");\r
- fCPtResolTan->SetYTitle("#Deltap_{t}/p_{t}");\r
-\r
- fCPhiPullTan = new TH2F("CPhiPullTan","CPhiPullTan",50, -2,2,200,-5,5); \r
- fCPhiPullTan->SetXTitle("Tan(#theta)");\r
- fCPhiPullTan->SetYTitle("#Delta#phi/#Sigma");\r
-\r
- fCTanPullTan = new TH2F("CTanPullTan","CTanPullTan",50, -2,2,200,-5,5);\r
- fCTanPullTan->SetXTitle("Tan(#theta)");\r
- fCTanPullTan->SetYTitle("#Delta#theta/#Sigma");\r
-\r
- fCPtPullTan=new TH2F("CPtPull","CPtPull",50, -2,2,200,-5,5); \r
- fCPtPullTan->SetXTitle("Tan(#theta)");\r
- fCPtPullTan->SetYTitle("(1/mcp_{t}-1/p_{t})/#Sigma");\r
-\r
- fPtResolLPT = new TH2F("Pt_resol_lpt","pt resol",10, 0.1,3,200,-0.2,0.2);\r
- fPtResolLPT->SetXTitle("p_{t}");\r
- fPtResolLPT->SetYTitle("#Deltap_{t}/p_{t}");\r
-\r
- fPtResolHPT = new TH2F("Pt_resol_hpt","pt resol",10, 2,100,200,-0.3,0.3); \r
- fPtResolHPT->SetXTitle("p_{t}");\r
- fPtResolHPT->SetYTitle("#Deltap_{t}/p_{t}");\r
-\r
- fPtPullLPT = new TH2F("Pt_pull_lpt","pt pull",10, 0.1,3,200,-6,6);\r
- fPtPullLPT->SetXTitle("p_{t}");\r
- fPtPullLPT->SetYTitle("#Deltap_{t}/#Sigma");\r
-\r
- fPtPullHPT = new TH2F("Pt_pull_hpt","pt pull",10,2,100,200,-6,6); \r
- fPtPullHPT->SetXTitle("p_{t}");\r
- fPtPullHPT->SetYTitle("#Deltap_{t}/#Sigma");\r
-\r
- //\r
- // Parametrisation histograms\r
- // \r
-\r
- f1Pt2Resol1PtTPC = new TH2F("f1Pt2Resol1PtTPC","(1/mcpt-1/pt)/(1+1/mcpt)^2 vs 1/pt)",100,0,10,200,-0.020,0.020); \r
- f1Pt2Resol1PtTPC->SetXTitle("1/mcp_{t}");\r
- f1Pt2Resol1PtTPC->SetYTitle("(1/mcp_{t}-1/p_{t})/(1+1/mcp_{t})^2)");\r
-\r
- f1Pt2Resol1PtTPCITS = new TH2F("f1Pt2Resol1PtTPCITS","(1/mcpt-1/pt)/(1+1/mcpt)^2 vs 1/pt)",100,0,10,200,-0.020,0.020); \r
- f1Pt2Resol1PtTPCITS->SetXTitle("1/mcp_{t}");\r
- f1Pt2Resol1PtTPCITS->SetYTitle("(1/mcp_{t}-1/p_{t})/(1+1/mcp_{t})^2)");\r
-\r
- fYResol1PtTPC = new TH2F("fYResol1PtTPC","fYResol1PtTPC",100, 0,10,200,-0.010,0.010); \r
- fYResol1PtTPC->SetXTitle("1/mcpt");\r
- fYResol1PtTPC->SetYTitle("#DeltaY");\r
-\r
- fYResol1PtTPCITS = new TH2F("fYResol1PtTPCITS","fYResol1PtTPCITS",100, 0,10,200,-0.010,0.010); \r
- fYResol1PtTPCITS->SetXTitle("1/mcpt");\r
- fYResol1PtTPCITS->SetYTitle("#DeltaY");\r
-\r
- fZResol1PtTPC = new TH2F("fZResol1PtTPC","fZResol1PtTPC",50, 0,10,200,-0.020,0.020); \r
- fZResol1PtTPC->SetXTitle("1/mcpt");\r
- fZResol1PtTPC->SetYTitle("#DeltaZ");\r
-\r
- fZResol1PtTPCITS = new TH2F("fZResol1PtTPCITS","fZResol1PtTPCITS",50, 0,10,200,-0.020,0.020); \r
- fZResol1PtTPCITS->SetXTitle("1/mcpt");\r
- fZResol1PtTPCITS->SetYTitle("#DeltaZ");\r
-\r
- fPhiResol1PtTPC = new TH2F("fPhiResol1PtTPC","fPhiResol1PtTPC",50, 0,10,200,-0.005,0.005); \r
- fPhiResol1PtTPC->SetXTitle("1/mcpt");\r
- fPhiResol1PtTPC->SetYTitle("#Delta#phi");\r
-\r
- fPhiResol1PtTPCITS = new TH2F("fPhiResol1PtTPCITS","fPhiResol1PtTPCITS",50, 0,10,200,-0.005,0.005); \r
- fPhiResol1PtTPCITS->SetXTitle("1/mcpt");\r
- fPhiResol1PtTPCITS->SetYTitle("#Delta#phi");\r
-\r
- fThetaResol1PtTPC = new TH2F("fThetaResol1PtTPC","fThetaResol1PtTPC",50, 0,10,200,-0.005,0.005); \r
- fThetaResol1PtTPC->SetXTitle("1/mcpt");\r
- fThetaResol1PtTPC->SetYTitle("#Delta#theta");\r
-\r
- fThetaResol1PtTPCITS = new TH2F("fThetaResol1PtTPCITS","fThetaResol1PtTPCITS",50, 0,10,200,-0.005,0.005); \r
- fThetaResol1PtTPCITS->SetXTitle("1/mcpt");\r
- fThetaResol1PtTPCITS->SetYTitle("#Delta#theta");\r
- \r
- // constrained\r
- fC1Pt2Resol1PtTPC = new TH2F("fC1Pt2Resol1PtTPC","(1/mcpt-1/pt)/(1+1/mcpt)^2 vs 1/pt)",100,0,10,200,-0.020,0.020); \r
- fC1Pt2Resol1PtTPC->SetXTitle("1/mcp_{t}");\r
- fC1Pt2Resol1PtTPC->SetYTitle("(1/mcp_{t}-1/p_{t})/(1+1/mcp_{t})^2)");\r
-\r
- fC1Pt2Resol1PtTPCITS = new TH2F("fC1Pt2Resol1PtTPCITS","(1/mcpt-1/pt)/(1+1/mcpt)^2 vs 1/pt)",100,0,10,200,-0.020,0.020); \r
- fC1Pt2Resol1PtTPCITS->SetXTitle("1/mcp_{t}");\r
- fC1Pt2Resol1PtTPCITS->SetYTitle("(1/mcp_{t}-1/p_{t})/(1+1/mcp_{t})^2)");\r
-\r
- fCYResol1PtTPC = new TH2F("fCYResol1PtTPC","fCYResol1PtTPC",100, 0,10,200,-0.010,0.010); \r
- fCYResol1PtTPC->SetXTitle("1/mcpt");\r
- fCYResol1PtTPC->SetYTitle("#DeltaY");\r
-\r
- fCYResol1PtTPCITS = new TH2F("fCYResol1PtTPCITS","fCYResol1PtTPCITS",100, 0,10,200,-0.010,0.010); \r
- fCYResol1PtTPCITS->SetXTitle("1/mcpt");\r
- fCYResol1PtTPCITS->SetYTitle("#DeltaY");\r
-\r
- fCZResol1PtTPC = new TH2F("fCZResol1PtTPC","fCZResol1PtTPC",50, 0,10,200,-0.020,0.020); \r
- fCZResol1PtTPC->SetXTitle("1/mcpt");\r
- fCZResol1PtTPC->SetYTitle("#DeltaZ");\r
-\r
- fCZResol1PtTPCITS = new TH2F("fCZResol1PtTPCITS","fCZResol1PtTPCITS",50, 0,10,200,-0.020,0.020); \r
- fCZResol1PtTPCITS->SetXTitle("1/mcpt");\r
- fCZResol1PtTPCITS->SetYTitle("#DeltaZ");\r
-\r
- fCPhiResol1PtTPC = new TH2F("fCPhiResol1PtTPC","fCPhiResol1PtTPC",50, 0,10,200,-0.005,0.005); \r
- fCPhiResol1PtTPC->SetXTitle("1/mcpt");\r
- fCPhiResol1PtTPC->SetYTitle("#Delta#phi");\r
-\r
- fCPhiResol1PtTPCITS = new TH2F("fCPhiResol1PtTPCITS","fCPhiResol1PtTPCITS",50, 0,10,200,-0.005,0.005); \r
- fCPhiResol1PtTPCITS->SetXTitle("1/mcpt");\r
- fCPhiResol1PtTPCITS->SetYTitle("#Delta#phi");\r
-\r
- fCThetaResol1PtTPC = new TH2F("fCThetaResol1PtTPC","fCThetaResol1PtTPC",50, 0,10,200,-0.005,0.005); \r
- fCThetaResol1PtTPC->SetXTitle("1/mcpt");\r
- fCThetaResol1PtTPC->SetYTitle("#Delta#theta");\r
-\r
- fCThetaResol1PtTPCITS = new TH2F("fCThetaResol1PtTPCITS","fCThetaResol1PtTPCITS",50, 0,10,200,-0.005,0.005); \r
- fCThetaResol1PtTPCITS->SetXTitle("1/mcpt");\r
- fCThetaResol1PtTPCITS->SetYTitle("#Delta#theta");\r
-}\r
-\r
-//_____________________________________________________________________________\r
-void AliComparisonRes::InitCuts()\r
-{\r
- // Init cuts \r
- if(!fCutsMC) \r
- AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");\r
- if(!fCutsRC) \r
- AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");\r
-}\r
-\r
-//_____________________________________________________________________________\r
-void AliComparisonRes::Process(AliMCInfo* infoMC, AliESDRecInfo *infoRC)\r
-{\r
- // Fill resolution comparison information \r
- AliExternalTrackParam *track = 0;\r
- Double_t kRadius = 3.0; // beam pipe radius\r
- Double_t kMaxStep = 5.0; // max step\r
- Double_t field = AliTracker::GetBz(); // nominal Bz field [kG]\r
- Double_t kMaxD = 123456.0; // max distance\r
-\r
- Int_t clusterITS[200];\r
- Double_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z\r
-\r
- Float_t mcpt = infoMC->GetParticle().Pt();\r
-\r
- // distance to Prim. vertex \r
- const Double_t* dv = infoMC->GetVDist(); \r
- Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])<fCutsMC->GetMaxR() && TMath::Abs(dv[2])<fCutsMC->GetMaxVz();\r
-\r
- // Check selection cuts\r
- if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return; \r
- if (!isPrim) return;\r
- //if (infoRC->GetStatus(1)==0) return;\r
- if (infoRC->GetStatus(1)!=3) return; // TPC refit\r
- if (!infoRC->GetESDtrack()) return; \r
- if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;\r
-\r
- // calculate and set prim. vertex\r
- fVertex->SetXv( infoMC->GetParticle().Vx() - dv[0] );\r
- fVertex->SetYv( infoMC->GetParticle().Vy() - dv[1] );\r
- fVertex->SetZv( infoMC->GetParticle().Vz() - dv[2] );\r
-\r
- Float_t deltaPt= (mcpt-infoRC->GetESDtrack()->Pt())/mcpt; \r
- Float_t pullPt= (1/mcpt-infoRC->GetESDtrack()->OneOverPt())/TMath::Sqrt(infoRC->GetESDtrack()->GetSigma1Pt2()); \r
- Float_t deltaPhi = TMath::ATan2(infoRC->GetESDtrack()->Py(),infoRC->GetESDtrack()->Px())-\r
- TMath::ATan2(infoMC->GetParticle().Py(),infoMC->GetParticle().Px());\r
-\r
- Float_t deltaTan = TMath::ATan2(infoRC->GetESDtrack()->Pz(),infoRC->GetESDtrack()->Pt())-\r
- TMath::ATan2(infoMC->GetParticle().Pz(),infoMC->GetParticle().Pt());\r
-\r
- Float_t delta1Pt2 = (1/mcpt-infoRC->GetESDtrack()->OneOverPt())/TMath::Power(1+1/mcpt,2); \r
-\r
- Float_t deltaY1Pt = (infoMC->GetParticle().Vy()-infoRC->GetESDtrack()->GetY()) / (0.2+1/mcpt);\r
- Float_t deltaZ1Pt = (infoMC->GetParticle().Vz()-infoRC->GetESDtrack()->GetZ()) / (0.2+1/mcpt);\r
- Float_t deltaPhi1Pt = deltaPhi / (0.1+1/mcpt);\r
- Float_t deltaTheta1Pt = deltaTan / (0.1+1/mcpt);\r
-\r
- // calculate track parameters at vertex\r
- if (infoRC->GetESDtrack()->GetTPCInnerParam())\r
- {\r
- if ((track = new AliExternalTrackParam(*infoRC->GetESDtrack()->GetTPCInnerParam())) != 0 )\r
- {\r
- Bool_t bStatus = AliTracker::PropagateTrackTo(track,kRadius,infoMC->GetMass(),kMaxStep,kTRUE);\r
- Bool_t bDCAStatus = track->PropagateToDCA(fVertex,field,kMaxD,dca,cov);\r
-\r
- // Fill parametrisation histograms (only TPC track)\r
- if(bStatus && bDCAStatus) \r
- {\r
- f1Pt2Resol1PtTPC->Fill(1/mcpt,delta1Pt2);\r
- fYResol1PtTPC->Fill(1/mcpt,deltaY1Pt);\r
- fZResol1PtTPC->Fill(1/mcpt,deltaZ1Pt);\r
- fPhiResol1PtTPC->Fill(1/mcpt,deltaPhi1Pt);\r
- fThetaResol1PtTPC->Fill(1/mcpt,deltaTheta1Pt);\r
- }\r
- delete track;\r
- }\r
- }\r
-\r
- // TPC and ITS (nb. of clusters >2) in the system\r
- if(infoRC->GetESDtrack()->GetITSclusters(clusterITS)>2) \r
- {\r
- f1Pt2Resol1PtTPCITS->Fill(1/mcpt,delta1Pt2);\r
- fYResol1PtTPCITS->Fill(1/mcpt,deltaY1Pt);\r
- fZResol1PtTPCITS->Fill(1/mcpt,deltaZ1Pt);\r
- fPhiResol1PtTPCITS->Fill(1/mcpt,deltaPhi1Pt);\r
- fThetaResol1PtTPCITS->Fill(1/mcpt,deltaTheta1Pt);\r
- }\r
-\r
- // Fill histograms\r
- fPtResolLPT->Fill(mcpt,deltaPt);\r
- fPtResolHPT->Fill(mcpt,deltaPt);\r
- fPtPullLPT->Fill(mcpt,pullPt);\r
- fPtPullHPT->Fill(mcpt,pullPt); \r
-}\r
-\r
-//_____________________________________________________________________________\r
-void AliComparisonRes::ProcessConstrained(AliMCInfo* infoMC, AliESDRecInfo *infoRC)\r
-{\r
- // Fill resolution comparison information (constarained parameters) \r
- //\r
- AliExternalTrackParam *track = 0;\r
- Double_t kRadius = 3.0; // beam pipe radius\r
- Double_t kMaxStep = 5.0; // max step\r
- Double_t field = AliTracker::GetBz(); // nominal Bz field [kG]\r
- Double_t kMaxD = 123456.0; // max distance\r
-\r
- Int_t clusterITS[200];\r
- Double_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z\r
- \r
- Float_t mcpt = infoMC->GetParticle().Pt();\r
- Float_t tantheta = TMath::Tan(infoMC->GetParticle().Theta()-TMath::Pi()*0.5);\r
-\r
- // distance to Prim. vertex \r
- const Double_t* dv = infoMC->GetVDist(); \r
- Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])<fCutsMC->GetMaxR() && TMath::Abs(dv[2])<fCutsMC->GetMaxVz();\r
- \r
- // Check selection cuts\r
- if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return; \r
- if (!isPrim) return;\r
- if (infoRC->GetStatus(1)!=3) return;\r
- if (!infoRC->GetESDtrack()) return; \r
- if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;\r
- if (!infoRC->GetESDtrack()->GetConstrainedParam()) return;\r
-\r
-// calculate and set prim. vertex\r
- fVertex->SetXv( infoMC->GetParticle().Vx() - dv[0] );\r
- fVertex->SetYv( infoMC->GetParticle().Vy() - dv[1] );\r
- fVertex->SetZv( infoMC->GetParticle().Vz() - dv[2] );\r
-\r
- // constrained parameters resolution\r
- const AliExternalTrackParam * cparam = infoRC->GetESDtrack()->GetConstrainedParam();\r
- Float_t deltaPt= (mcpt-cparam->Pt())/mcpt; \r
- Float_t pullPt= (1/mcpt-cparam->OneOverPt())/TMath::Sqrt(cparam->GetSigma1Pt2()); \r
- Float_t deltaPhi = TMath::ATan2(cparam->Py(),cparam->Px())-\r
- TMath::ATan2(infoMC->GetParticle().Py(),infoMC->GetParticle().Px());\r
- Float_t pullPhi = deltaPhi/TMath::Sqrt(cparam->GetSigmaSnp2()); \r
- Float_t deltaTan = TMath::ATan2(cparam->Pz(),cparam->Pt())-TMath::ATan2(infoMC->GetParticle().Pz(),infoMC->GetParticle().Pt());\r
- Float_t pullTan = deltaPhi/TMath::Sqrt(cparam->GetSigmaSnp2()); \r
-\r
-\r
- Float_t delta1Pt2 = (1/mcpt-cparam->OneOverPt())/TMath::Power(1+1/mcpt,2); \r
-\r
- Float_t deltaY1Pt = (infoMC->GetParticle().Vy()-cparam->GetY()) / (0.2+1/mcpt);\r
- Float_t deltaZ1Pt = (infoMC->GetParticle().Vz()-cparam->GetZ()) / (0.2+1/mcpt);\r
- Float_t deltaPhi1Pt = deltaPhi / (0.1+1/mcpt);\r
- Float_t deltaTheta1Pt = deltaTan / (0.1+1/mcpt);\r
-\r
- // calculate track parameters at vertex\r
- if (infoRC->GetESDtrack()->GetTPCInnerParam())\r
- {\r
- if ((track = new AliExternalTrackParam(*infoRC->GetESDtrack()->GetTPCInnerParam())) != 0 )\r
- {\r
- Bool_t bStatus = AliTracker::PropagateTrackTo(track,kRadius,infoMC->GetMass(),kMaxStep,kTRUE);\r
- Bool_t bDCAStatus = track->PropagateToDCA(fVertex,field,kMaxD,dca,cov);\r
-\r
- // Fill parametrisation histograms (only TPC track)\r
- if(bStatus && bDCAStatus) \r
- {\r
- fC1Pt2Resol1PtTPC->Fill(1/mcpt,delta1Pt2);\r
- fCYResol1PtTPC->Fill(1/mcpt,deltaY1Pt);\r
- fCZResol1PtTPC->Fill(1/mcpt,deltaZ1Pt);\r
- fCPhiResol1PtTPC->Fill(1/mcpt,deltaPhi1Pt);\r
- fCThetaResol1PtTPC->Fill(1/mcpt,deltaTheta1Pt);\r
- }\r
- delete track;\r
- }\r
- }\r
-\r
- // TPC and ITS (nb. of clusters >2) in the system\r
- if(infoRC->GetESDtrack()->GetITSclusters(clusterITS)>2) \r
- {\r
- fC1Pt2Resol1PtTPCITS->Fill(1/mcpt,delta1Pt2);\r
- fCYResol1PtTPCITS->Fill(1/mcpt,deltaY1Pt);\r
- fCZResol1PtTPCITS->Fill(1/mcpt,deltaZ1Pt);\r
- fCPhiResol1PtTPCITS->Fill(1/mcpt,deltaPhi1Pt);\r
- fCThetaResol1PtTPCITS->Fill(1/mcpt,deltaTheta1Pt);\r
- }\r
-\r
- // Fill histograms\r
- fCPtResolTan->Fill(tantheta,deltaPt);\r
- fCPtPullTan->Fill(tantheta,pullPt);\r
- fCPhiResolTan->Fill(tantheta,deltaPhi);\r
- fCPhiPullTan->Fill(tantheta,pullPhi);\r
- fCTanResolTan->Fill(tantheta,deltaTan);\r
- fCTanPullTan->Fill(tantheta,pullTan);\r
-}\r
- \r
-//_____________________________________________________________________________\r
-void AliComparisonRes::Exec(AliMCInfo* infoMC, AliESDRecInfo *infoRC){\r
- \r
- // Process comparison information \r
- Process(infoMC,infoRC);\r
- ProcessConstrained(infoMC,infoRC);\r
-}\r
-\r
-//_____________________________________________________________________________\r
-TH1F* AliComparisonRes::MakeResol(TH2F * his, Int_t integ, Bool_t type){\r
- // Create resolution histograms\r
- \r
- TH1F *hisr, *hism;\r
- if (!gPad) new TCanvas;\r
- hisr = AliTreeDraw::CreateResHistoI(his,&hism,integ);\r
- if (type) return hism;\r
- else \r
- return hisr;\r
-}\r
-\r
-//_____________________________________________________________________________\r
-void AliComparisonRes::Analyse(){\r
- // Analyse comparison information and store output histograms \r
- // in the "pictures_res.root" file \r
- \r
- AliComparisonRes * comp=this;\r
- TH1F *hiss=0;\r
-\r
- TFile *fp = new TFile("pictures_res.root","recreate");\r
- fp->cd();\r
-\r
- TCanvas * c = new TCanvas("Phi resol Tan","Phi resol Tan");\r
- c->cd();\r
- //\r
- hiss = comp->MakeResol(comp->fCPtResolTan,1,0);\r
- hiss->SetXTitle("Tan(#theta)");\r
- hiss->SetYTitle("#sigmap_{t}/p_{t}");\r
- hiss->Draw(); \r
- hiss->Write("CptResolTan");\r
- //\r
- hiss = comp->MakeResol(comp->fCPhiResolTan,1,0);\r
- hiss->SetXTitle("Tan(#theta)");\r
- hiss->SetYTitle("#sigma#phi (rad)");\r
- hiss->Draw();\r
- hiss->Write("PhiResolTan");\r
- //\r
- hiss = comp->MakeResol(comp->fCTanResolTan,1,0);\r
- hiss->SetXTitle("Tan(#theta)");\r
- hiss->SetYTitle("#sigma#theta (rad)");\r
- hiss->Draw();\r
- hiss->Write("ThetaResolTan");\r
- //\r
- hiss = comp->MakeResol(comp->fCPtPullTan,1,0);\r
- hiss->SetXTitle("Tan(#theta)");\r
- hiss->SetYTitle("1/mcp_{t}-1/p_{t}/#Sigma(1/p_{t})");\r
- hiss->Draw();\r
- hiss->Write("CptPullTan");\r
- //\r
- hiss = comp->MakeResol(comp->fC1Pt2Resol1PtTPC,1,0);\r
- hiss->SetXTitle("1/mcp_{t}");\r
- hiss->SetYTitle("1/mcp_{t}-1/p_{t}/(1+1/p_{t})^2");\r
- hiss->Draw();\r
- hiss->Write("C1Pt2Resol1PtTPC");\r
- fC1Pt2Resol1PtTPC->Write();\r
-\r
- hiss = comp->MakeResol(comp->fC1Pt2Resol1PtTPCITS,1,0);\r
- hiss->SetXTitle("1/mcp_{t}");\r
- hiss->SetYTitle("1/mcp_{t}-1/p_{t}/(1+1/p_{t})^2");\r
- hiss->Draw();\r
- hiss->Write("C1Pt2Resol1PtTPCITS");\r
- fC1Pt2Resol1PtTPCITS->Write();\r
- //\r
- hiss = comp->MakeResol(comp->fCYResol1PtTPC,1,0);\r
- hiss->SetXTitle("1/mcp_{t}");\r
- hiss->SetYTitle("(mcy-y)/(0.2+1/mcp_{t})");\r
- hiss->Draw();\r
- hiss->Write("CYResol1PtTPC");\r
- fCYResol1PtTPC->Write();\r
-\r
- hiss = comp->MakeResol(comp->fCYResol1PtTPCITS,1,0);\r
- hiss->SetXTitle("1/mcp_{t}");\r
- hiss->SetYTitle("(mcy-y)/(0.2+1/mcp_{t})");\r
- hiss->Draw();\r
- hiss->Write("CYResol1PtTPCITS");\r
- fCYResol1PtTPCITS->Write();\r
- //\r
- hiss = comp->MakeResol(comp->fCZResol1PtTPC,1,0);\r
- hiss->SetXTitle("1/mcp_{t}");\r
- hiss->SetYTitle("(mcz-z)/(0.2+1/mcp_{t})");\r
- hiss->Draw();\r
- hiss->Write("CZResol1PtTPC");\r
- fCZResol1PtTPC->Write();\r
-\r
- hiss = comp->MakeResol(comp->fCZResol1PtTPCITS,1,0);\r
- hiss->SetXTitle("1/mcp_{t}");\r
- hiss->SetYTitle("(mcz-z)/(0.2+1/mcp_{t})");\r
- hiss->Draw();\r
- hiss->Write("CZResol1PtTPCITS");\r
- fCZResol1PtTPCITS->Write();\r
- //\r
- hiss = comp->MakeResol(comp->fCPhiResol1PtTPC,1,0);\r
- hiss->SetXTitle("1/mcp_{t}");\r
- hiss->SetYTitle("(mc#phi-#phi)/(0.1+1/mcp_{t})");\r
- hiss->Draw();\r
- hiss->Write("CPhiResol1PtTPC");\r
- fCPhiResol1PtTPC->Write();\r
-\r
- hiss = comp->MakeResol(comp->fCPhiResol1PtTPCITS,1,0);\r
- hiss->SetXTitle("1/mcp_{t}");\r
- hiss->SetYTitle("(mc#phi-#phi)/(0.1+1/mcp_{t})");\r
- hiss->Draw();\r
- hiss->Write("CPhiResol1PtTPCITS");\r
- fCPhiResol1PtTPCITS->Write();\r
- //\r
- hiss = comp->MakeResol(comp->fCThetaResol1PtTPC,1,0);\r
- hiss->SetXTitle("1/mcp_{t}");\r
- hiss->SetYTitle("(mc#theta-#theta)/(0.1+1/mcp_{t})");\r
- hiss->Draw();\r
- hiss->Write("CThetaResol1PtTPC");\r
- fCThetaResol1PtTPC->Write();\r
-\r
- hiss = comp->MakeResol(comp->fCThetaResol1PtTPCITS,1,0);\r
- hiss->SetXTitle("1/mcp_{t}");\r
- hiss->SetYTitle("(mc#theta-#theta)/(0.1+1/mcp_{t})");\r
- hiss->Draw();\r
- hiss->Write("CThetaResol1PtTPCITS");\r
- fCThetaResol1PtTPCITS->Write();\r
-\r
- //\r
- hiss = comp->MakeResol(comp->f1Pt2Resol1PtTPC,1,0);\r
- hiss->SetXTitle("1/mcp_{t}");\r
- hiss->SetYTitle("1/mcp_{t}-1/p_{t}/(1+1/p_{t})^2");\r
- hiss->Draw();\r
- hiss->Write("OnePt2Resol1PtTPC");\r
- f1Pt2Resol1PtTPC->Write();\r
-\r
- hiss = comp->MakeResol(comp->f1Pt2Resol1PtTPCITS,1,0);\r
- hiss->SetXTitle("1/mcp_{t}");\r
- hiss->SetYTitle("1/mcp_{t}-1/p_{t}/(1+1/p_{t})^2");\r
- hiss->Draw();\r
- hiss->Write("OnePt2Resol1PtTPCITS");\r
- f1Pt2Resol1PtTPCITS->Write();\r
- //\r
- hiss = comp->MakeResol(comp->fYResol1PtTPC,1,0);\r
- hiss->SetXTitle("1/mcp_{t}");\r
- hiss->SetYTitle("(mcy-y)/(0.2+1/mcp_{t})");\r
- hiss->Draw();\r
- hiss->Write("YResol1PtTPC");\r
- fYResol1PtTPC->Write();\r
-\r
- hiss = comp->MakeResol(comp->fYResol1PtTPCITS,1,0);\r
- hiss->SetXTitle("1/mcp_{t}");\r
- hiss->SetYTitle("(mcy-y)/(0.2+1/mcp_{t})");\r
- hiss->Draw();\r
- hiss->Write("YResol1PtTPCITS");\r
- fYResol1PtTPCITS->Write();\r
- //\r
- hiss = comp->MakeResol(comp->fZResol1PtTPC,1,0);\r
- hiss->SetXTitle("1/mcp_{t}");\r
- hiss->SetYTitle("(mcz-z)/(0.2+1/mcp_{t})");\r
- hiss->Draw();\r
- hiss->Write("ZResol1PtTPC");\r
- fZResol1PtTPC->Write();\r
-\r
- hiss = comp->MakeResol(comp->fZResol1PtTPCITS,1,0);\r
- hiss->SetXTitle("1/mcp_{t}");\r
- hiss->SetYTitle("(mcz-z)/(0.2+1/mcp_{t})");\r
- hiss->Draw();\r
- hiss->Write("ZResol1PtTPCITS");\r
- fZResol1PtTPCITS->Write();\r
- //\r
- hiss = comp->MakeResol(comp->fPhiResol1PtTPC,1,0);\r
- hiss->SetXTitle("1/mcp_{t}");\r
- hiss->SetYTitle("(mc#phi-#phi)/(0.1+1/mcp_{t})");\r
- hiss->Draw();\r
- hiss->Write("PhiResol1PtTPC");\r
- fPhiResol1PtTPC->Write();\r
-\r
- hiss = comp->MakeResol(comp->fPhiResol1PtTPCITS,1,0);\r
- hiss->SetXTitle("1/mcp_{t}");\r
- hiss->SetYTitle("(mc#phi-#phi)/(0.1+1/mcp_{t})");\r
- hiss->Draw();\r
- hiss->Write("PhiResol1PtTPCITS");\r
- fPhiResol1PtTPCITS->Write();\r
- //\r
- hiss = comp->MakeResol(comp->fThetaResol1PtTPC,1,0);\r
- hiss->SetXTitle("1/mcp_{t}");\r
- hiss->SetYTitle("(mc#theta-#theta)/(0.1+1/mcp_{t})");\r
- hiss->Draw();\r
- hiss->Write("ThetaResol1PtTPC");\r
- fThetaResol1PtTPC->Write();\r
-\r
- hiss = comp->MakeResol(comp->fThetaResol1PtTPCITS,1,0);\r
- hiss->SetXTitle("1/mcp_{t}");\r
- hiss->SetYTitle("(mc#theta-#theta)/(0.1+1/mcp_{t})");\r
- hiss->Draw();\r
- hiss->Write("ThetaResol1PtTPCITS");\r
- fThetaResol1PtTPCITS->Write();\r
-\r
- fp->Close();\r
-}\r
-\r
-//_____________________________________________________________________________\r
-Long64_t AliComparisonRes::Merge(TCollection* list) \r
-{\r
- // Merge list of objects (needed by PROOF)\r
-\r
- if (!list)\r
- return 0;\r
-\r
- if (list->IsEmpty())\r
- return 1;\r
-\r
- TIterator* iter = list->MakeIterator();\r
- TObject* obj = 0;\r
-\r
- // collection of generated histograms\r
- Int_t count=0;\r
- while((obj = iter->Next()) != 0) \r
- {\r
- AliComparisonRes* entry = dynamic_cast<AliComparisonRes*>(obj);\r
- if (entry == 0) continue; \r
-\r
- fPtResolLPT->Add(entry->fPtResolLPT);\r
- fPtResolHPT->Add(entry->fPtResolHPT);\r
- fPtPullLPT->Add(entry->fPtPullLPT);\r
- fPtPullHPT->Add(entry->fPtPullHPT);\r
-\r
- // Histograms for 1/pt parameterisation\r
- f1Pt2Resol1PtTPC->Add(entry->f1Pt2Resol1PtTPC);\r
- fYResol1PtTPC->Add(entry->fYResol1PtTPC);\r
- fZResol1PtTPC->Add(entry->fZResol1PtTPC);\r
- fPhiResol1PtTPC->Add(entry->fPhiResol1PtTPC);\r
- fThetaResol1PtTPC->Add(entry->fThetaResol1PtTPC);\r
-\r
- f1Pt2Resol1PtTPCITS->Add(entry->f1Pt2Resol1PtTPCITS);\r
- fYResol1PtTPCITS->Add(entry->fYResol1PtTPCITS);\r
- fZResol1PtTPCITS->Add(entry->fZResol1PtTPCITS);\r
- fPhiResol1PtTPCITS->Add(entry->fPhiResol1PtTPCITS);\r
- fThetaResol1PtTPCITS->Add(entry->fThetaResol1PtTPCITS);\r
-\r
- // Resolution histograms (constrained param)\r
- fCPhiResolTan->Add(entry->fCPhiResolTan);\r
- fCTanResolTan->Add(entry->fCTanResolTan);\r
- fCPtResolTan->Add(entry->fCPtResolTan);\r
- fCPhiPullTan->Add(entry->fCPhiPullTan);\r
- fCTanPullTan->Add(entry->fCTanPullTan);\r
- fCPtPullTan->Add(entry->fCPtPullTan);\r
-\r
- // Histograms for 1/pt parameterisation (constrained)\r
- fC1Pt2Resol1PtTPC->Add(entry->fC1Pt2Resol1PtTPC);\r
- fCYResol1PtTPC->Add(entry->fCYResol1PtTPC);\r
- fCZResol1PtTPC->Add(entry->fCZResol1PtTPC);\r
- fCPhiResol1PtTPC->Add(entry->fCPhiResol1PtTPC);\r
- fCThetaResol1PtTPC->Add(entry->fCThetaResol1PtTPC);\r
-\r
- fC1Pt2Resol1PtTPCITS->Add(entry->fC1Pt2Resol1PtTPCITS);\r
- fCYResol1PtTPCITS->Add(entry->fCYResol1PtTPCITS);\r
- fCZResol1PtTPCITS->Add(entry->fCZResol1PtTPCITS);\r
- fCPhiResol1PtTPCITS->Add(entry->fCPhiResol1PtTPCITS);\r
- fCThetaResol1PtTPCITS->Add(entry->fCThetaResol1PtTPCITS);\r
-\r
- count++;\r
- }\r
-\r
-return count;\r
-}\r
+//------------------------------------------------------------------------------
+// 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/graphs) are stored in the folder which is
+// a data member of AliComparisonRes.
+//
+// Author: J.Otwinowski 04/02/2008
+//------------------------------------------------------------------------------
+
+/*
+
+ // after running comparison task, read the file, and get component
+ gROOT->LoadMacro("$ALICE_ROOT/PWG1/Macros/LoadMyLibs.C");
+ LoadMyLibs();
+
+ TFile f("Output.root");
+ AliComparisonRes * compObj = (AliComparisonRes*)f.Get("AliComparisonRes");
+
+ // analyse comparison data
+ compObj->Analyse();
+
+ // the output histograms/graphs will be stored in the folder "folderRes"
+ compObj->GetAnalysisFolder()->ls("*");
+
+ // user can save whole comparison object (or only folder with anlysed histograms)
+ // in the seperate output file (e.g.)
+ TFile fout("Analysed_Res.root","recreate");
+ compObj->Write(); // compObj->GetAnalysisFolder()->Write();
+ fout.Close();
+
+*/
+
+#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 "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():
+ AliComparisonObject("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
+ fPhiResolTan(0), //-> angular resolution
+ fTanResolTan(0), //-> angular resolution
+ fPhiPullTan(0), //-> angular resolution
+ fTanPullTan(0), //-> angular resolution
+
+ //
+ // 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
+ //
+
+ f1Pt2ResolS1PtTPC(0),
+ f1Pt2ResolS1PtTPCITS(0),
+ fYResolS1PtTPC(0),
+ fYResolS1PtTPCITS(0),
+ fZResolS1PtTPC(0),
+ fZResolS1PtTPCITS(0),
+ fPhiResolS1PtTPC(0),
+ fPhiResolS1PtTPCITS(0),
+ fThetaResolS1PtTPC(0),
+ fThetaResolS1PtTPCITS(0),
+
+ // constrained
+ fC1Pt2ResolS1PtTPC(0),
+ fC1Pt2ResolS1PtTPCITS(0),
+ fCYResolS1PtTPC(0),
+ fCYResolS1PtTPCITS(0),
+ fCZResolS1PtTPC(0),
+ fCZResolS1PtTPCITS(0),
+ fCPhiResolS1PtTPC(0),
+ fCPhiResolS1PtTPCITS(0),
+ fCThetaResolS1PtTPC(0),
+ fCThetaResolS1PtTPCITS(0),
+
+ // vertex
+ fVertex(0),
+
+ // Cuts
+ fCutsRC(0),
+ fCutsMC(0),
+
+ // histogram folder
+ fAnalysisFolder(0)
+{
+ Init();
+
+ // 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;
+
+ if(fPhiResolTan) delete fPhiResolTan; fPhiResolTan=0;
+ if(fTanResolTan) delete fTanResolTan; fTanResolTan=0;
+ if(fPhiPullTan) delete fPhiPullTan; fPhiPullTan=0;
+ if(fTanPullTan) delete fTanPullTan; fTanPullTan=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(f1Pt2ResolS1PtTPC) delete f1Pt2ResolS1PtTPC; f1Pt2ResolS1PtTPC=0;
+ if(f1Pt2ResolS1PtTPCITS) delete f1Pt2ResolS1PtTPCITS; f1Pt2ResolS1PtTPCITS=0;
+ if(fYResolS1PtTPC) delete fYResolS1PtTPC; fYResolS1PtTPC=0;
+ if(fYResolS1PtTPCITS) delete fYResolS1PtTPCITS; fYResolS1PtTPCITS=0;
+ if(fZResolS1PtTPC) delete fZResolS1PtTPC; fZResolS1PtTPC=0;
+ if(fZResolS1PtTPCITS) delete fZResolS1PtTPCITS; fZResolS1PtTPCITS=0;
+ if(fPhiResolS1PtTPC) delete fPhiResolS1PtTPC; fPhiResolS1PtTPC=0;
+ if(fPhiResolS1PtTPCITS) delete fPhiResolS1PtTPCITS; fPhiResolS1PtTPCITS=0;
+ if(fThetaResolS1PtTPC) delete fThetaResolS1PtTPC; fThetaResolS1PtTPC=0;
+ if(fThetaResolS1PtTPCITS) delete fThetaResolS1PtTPCITS; fThetaResolS1PtTPCITS=0;
+
+ // constrained
+ if(fC1Pt2ResolS1PtTPC) delete fC1Pt2ResolS1PtTPC; fC1Pt2ResolS1PtTPC=0;
+ if(fC1Pt2ResolS1PtTPCITS) delete fC1Pt2ResolS1PtTPCITS; fC1Pt2ResolS1PtTPCITS=0;
+ if(fCYResolS1PtTPC) delete fCYResolS1PtTPC; fCYResolS1PtTPC=0;
+ if(fCYResolS1PtTPCITS) delete fCYResolS1PtTPCITS; fCYResolS1PtTPCITS=0;
+ if(fCZResolS1PtTPC) delete fCZResolS1PtTPC; fCZResolS1PtTPC=0;
+ if(fCZResolS1PtTPCITS) delete fCZResolS1PtTPCITS; fCZResolS1PtTPCITS=0;
+ if(fCPhiResolS1PtTPC) delete fCPhiResolS1PtTPC; fCPhiResolS1PtTPC=0;
+ if(fCPhiResolS1PtTPCITS) delete fCPhiResolS1PtTPCITS; fCPhiResolS1PtTPCITS=0;
+ if(fCThetaResolS1PtTPC) delete fCThetaResolS1PtTPC; fCThetaResolS1PtTPC=0;
+ if(fCThetaResolS1PtTPCITS) delete fCThetaResolS1PtTPCITS; fCThetaResolS1PtTPCITS=0;
+
+ if(fVertex) delete fVertex; fVertex=0;
+
+ if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
+}
+
+//_____________________________________________________________________________
+void AliComparisonRes::Init(){
+
+ // 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");
+
+ fPhiResolTan = new TH2F("PhiResolTan","PhiResolTan",50, -2,2,200,-0.025,0.025);
+ fPhiResolTan->SetXTitle("tan(#theta)");
+ fPhiResolTan->SetYTitle("#Delta#phi");
+
+ fTanResolTan = new TH2F("TanResolTan","TanResolTan",50, -2,2,200,-0.025,0.025);
+ fTanResolTan->SetXTitle("tan(#theta)");
+ fTanResolTan->SetYTitle("#Delta#theta");
+
+ fPhiPullTan = new TH2F("PhiPullTan","PhiPullTan",50, -2,2,200,-5,5);
+ fPhiPullTan->SetXTitle("Tan(#theta)");
+ fPhiPullTan->SetYTitle("#Delta#phi/#Sigma");
+
+ fTanPullTan = new TH2F("TanPullTan","TanPullTan",50, -2,2,200,-5,5);
+ fTanPullTan->SetXTitle("Tan(#theta)");
+ fTanPullTan->SetYTitle("#Delta#theta/#Sigma");
+
+ //
+ // Parametrisation histograms
+ //
+
+ f1Pt2ResolS1PtTPC = new TH2F("f1Pt2ResolS1PtTPC","(1/mcpt-1/pt)/(1+1/mcpt)^2 vs sqrt(1/pt))",100,0,3,200,-0.010,0.010);
+ f1Pt2ResolS1PtTPC->SetXTitle("#sqrt{1/mcp_{t}}");
+ f1Pt2ResolS1PtTPC->SetYTitle("(1/mcp_{t}-1/p_{t})/(1+1/mcp_{t})^2)");
+
+ f1Pt2ResolS1PtTPCITS = new TH2F("f1Pt2ResolS1PtTPCITS","(1/mcpt-1/pt)/(1+1/mcpt)^2 vs sqrt(1/pt))",100,0,3,200,-0.010,0.010);
+ f1Pt2ResolS1PtTPCITS->SetXTitle("#sqrt{1/mcp_{t}}");
+ f1Pt2ResolS1PtTPCITS->SetYTitle("(1/mcp_{t}-1/p_{t})/(1+1/mcp_{t})^2)");
+
+ fYResolS1PtTPC = new TH2F("fYResolS1PtTPC","fYResolS1PtTPC",100, 0,3,200,-1.0,1.0);
+ fYResolS1PtTPC->SetXTitle("#sqrt{1/mcp_{t}}");
+ fYResolS1PtTPC->SetYTitle("#DeltaY");
+
+ fYResolS1PtTPCITS = new TH2F("fYResolS1PtTPCITS","fYResolS1PtTPCITS",100, 0,3,200,-0.05,0.05);
+ fYResolS1PtTPCITS->SetXTitle("#sqrt{1/mcp_{t}}");
+ fYResolS1PtTPCITS->SetYTitle("#DeltaY");
+
+ fZResolS1PtTPC = new TH2F("fZResolS1PtTPC","fZResolS1PtTPC",100, 0,3,200,-1.0,1.0);
+ fZResolS1PtTPC->SetXTitle("#sqrt{1/mcp_{t}}");
+ fZResolS1PtTPC->SetYTitle("#DeltaZ");
+
+ fZResolS1PtTPCITS = new TH2F("fZResolS1PtTPCITS","fZResolS1PtTPCITS",100, 0,3,200,-0.05,0.05);
+ fZResolS1PtTPCITS->SetXTitle("#sqrt{1/mcp_{t}}");
+ fZResolS1PtTPCITS->SetYTitle("#DeltaZ");
+
+ fPhiResolS1PtTPC = new TH2F("fPhiResolS1PtTPC","fPhiResolS1PtTPC",100, 0,3,200,-0.025,0.025);
+ fPhiResolS1PtTPC->SetXTitle("#sqrt{1/mcp_{t}}");
+ fPhiResolS1PtTPC->SetYTitle("#Delta#phi");
+
+ fPhiResolS1PtTPCITS = new TH2F("fPhiResolS1PtTPCITS","fPhiResolS1PtTPCITS",100, 0,3,200,-0.01,0.01);
+ fPhiResolS1PtTPCITS->SetXTitle("#sqrt{1/mcp_{t}}");
+ fPhiResolS1PtTPCITS->SetYTitle("#Delta#phi");
+
+ fThetaResolS1PtTPC = new TH2F("fThetaResolS1PtTPC","fThetaResolS1PtTPC",100, 0,3,200,-0.025,0.025);
+ fThetaResolS1PtTPC->SetXTitle("#sqrt{1/mcp_{t}}");
+ fThetaResolS1PtTPC->SetYTitle("#Delta#theta");
+
+ fThetaResolS1PtTPCITS = new TH2F("fThetaResolS1PtTPCITS","fThetaResolS1PtTPCITS",100, 0,3,200,-0.01,0.01);
+ fThetaResolS1PtTPCITS->SetXTitle("#sqrt{1/mcp_{t}}");
+ fThetaResolS1PtTPCITS->SetYTitle("#Delta#theta");
+
+ // constrained
+ fC1Pt2ResolS1PtTPC = new TH2F("fC1Pt2ResolS1PtTPC","(1/mcpt-1/pt)/(1+1/mcpt)^2 vs 1/pt)",100,0,3,200,-0.010,0.010);
+ fC1Pt2ResolS1PtTPC->SetXTitle("#sqrt{1/mcp_{t}}");
+ fC1Pt2ResolS1PtTPC->SetYTitle("(1/mcp_{t}-1/p_{t})/(1+1/mcp_{t})^2)");
+
+ fC1Pt2ResolS1PtTPCITS = new TH2F("fC1Pt2ResolS1PtTPCITS","(1/mcpt-1/pt)/(1+1/mcpt)^2 vs 1/pt)",100,0,3,200,-0.010,0.010);
+ fC1Pt2ResolS1PtTPCITS->SetXTitle("#sqrt{1/mcp_{t}}");
+ fC1Pt2ResolS1PtTPCITS->SetYTitle("(1/mcp_{t}-1/p_{t})/(1+1/mcp_{t})^2)");
+
+ fCYResolS1PtTPC = new TH2F("fCYResolS1PtTPC","fCYResolS1PtTPC",100, 0,3,200,-1.0,1.0);
+ fCYResolS1PtTPC->SetXTitle("#sqrt{1/mcp_{t}}");
+ fCYResolS1PtTPC->SetYTitle("#DeltaY");
+
+ fCYResolS1PtTPCITS = new TH2F("fCYResolS1PtTPCITS","fCYResolS1PtTPCITS",100, 0,3,200,-0.01,0.01);
+ fCYResolS1PtTPCITS->SetXTitle("#sqrt{1/mcp_{t}}");
+ fCYResolS1PtTPCITS->SetYTitle("#DeltaY");
+
+ fCZResolS1PtTPC = new TH2F("fCZResolS1PtTPC","fCZResolS1PtTPC",100, 0,3,200,-1.0,1.0);
+ fCZResolS1PtTPC->SetXTitle("#sqrt{1/mcp_{t}}");
+ fCZResolS1PtTPC->SetYTitle("#DeltaZ");
+
+ fCZResolS1PtTPCITS = new TH2F("fCZResolS1PtTPCITS","fCZResolS1PtTPCITS",100, 0,3,200,-0.025,0.025);
+ fCZResolS1PtTPCITS->SetXTitle("#sqrt{1/mcp_{t}}");
+ fCZResolS1PtTPCITS->SetYTitle("#DeltaZ");
+
+ fCPhiResolS1PtTPC = new TH2F("fCPhiResolS1PtTPC","fCPhiResolS1PtTPC",100, 0,3,200,-0.025,0.025);
+ fCPhiResolS1PtTPC->SetXTitle("#sqrt{1/mcp_{t}}");
+ fCPhiResolS1PtTPC->SetYTitle("#Delta#phi");
+
+ fCPhiResolS1PtTPCITS = new TH2F("fCPhiResolS1PtTPCITS","fCPhiResolS1PtTPCITS",100, 0,3,200,-0.003,0.003);
+ fCPhiResolS1PtTPCITS->SetXTitle("#sqrt{1/mcp_{t}}");
+ fCPhiResolS1PtTPCITS->SetYTitle("#Delta#phi");
+
+ fCThetaResolS1PtTPC = new TH2F("fCThetaResolS1PtTPC","fCThetaResolS1PtTPC",100, 0,3,200,-0.025,0.025);
+ fCThetaResolS1PtTPC->SetXTitle("#sqrt{1/mcp_{t}}");
+ fCThetaResolS1PtTPC->SetYTitle("#Delta#theta");
+
+ fCThetaResolS1PtTPCITS = new TH2F("fCThetaResolS1PtTPCITS","fCThetaResolS1PtTPCITS",100, 0,3,200,-0.005,0.005);
+ fCThetaResolS1PtTPCITS->SetXTitle("#sqrt{1/mcp_{t}}");
+ fCThetaResolS1PtTPCITS->SetYTitle("#Delta#theta");
+
+ // Init cuts
+ if(!fCutsMC)
+ AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
+ if(!fCutsRC)
+ AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
+
+ // init folder
+ fAnalysisFolder = CreateFolder("folderRes","Analysis Resolution Folder");
+}
+
+//_____________________________________________________________________________
+void AliComparisonRes::Process(AliMCInfo* infoMC, AliESDRecInfo *infoRC)
+{
+ // Fill resolution comparison information
+ AliExternalTrackParam *track = 0;
+ 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 deltaPt, pullPt, deltaPhi,pullPhi, deltaTan, pullTan, delta1Pt2, deltaY1Pt, deltaZ1Pt, deltaPhi1Pt, deltaTheta1Pt;
+ Float_t deltaPtTPC, pullPtTPC, deltaPhiTPC, deltaTanTPC, delta1Pt2TPC, deltaY1PtTPC, deltaZ1PtTPC, deltaPhi1PtTPC, deltaTheta1PtTPC;
+
+ Float_t mcpt = infoMC->GetParticle().Pt();
+ Float_t s1mcpt = TMath::Sqrt(1./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])<fCutsMC->GetMaxR() && TMath::Abs(dv[2])<fCutsMC->GetMaxVz();
+
+ // Check selection cuts
+ if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return;
+ if (!isPrim) return;
+ if (infoRC->GetStatus(1)!=3) return; // TPC refit
+ if (!infoRC->GetESDtrack()) return;
+ if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->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] );
+
+ deltaPt= (mcpt-infoRC->GetESDtrack()->Pt())/mcpt;
+ pullPt= (1/mcpt-infoRC->GetESDtrack()->OneOverPt())/TMath::Sqrt(infoRC->GetESDtrack()->GetSigma1Pt2());
+ deltaPhi = TMath::ATan2(infoRC->GetESDtrack()->Py(),infoRC->GetESDtrack()->Px())-
+ TMath::ATan2(infoMC->GetParticle().Py(),infoMC->GetParticle().Px());
+ pullPhi = deltaPhi/TMath::Sqrt(infoRC->GetESDtrack()->GetSigmaSnp2());
+
+ deltaTan = TMath::ATan2(infoRC->GetESDtrack()->Pz(),infoRC->GetESDtrack()->Pt())-
+ TMath::ATan2(infoMC->GetParticle().Pz(),infoMC->GetParticle().Pt());
+ pullTan = deltaPhi/TMath::Sqrt(infoRC->GetESDtrack()->GetSigmaSnp2());
+
+ delta1Pt2 = (1/mcpt-infoRC->GetESDtrack()->OneOverPt())/TMath::Power(1+1/mcpt,2);
+ deltaY1Pt = (infoMC->GetParticle().Vy()-infoRC->GetESDtrack()->GetY()) / (0.2+1/mcpt);
+ deltaZ1Pt = (infoMC->GetParticle().Vz()-infoRC->GetESDtrack()->GetZ()) / (0.2+1/mcpt);
+ deltaPhi1Pt = deltaPhi / (0.1+1/mcpt);
+ deltaTheta1Pt = deltaTan / (0.1+1/mcpt);
+
+ // calculate track parameters at vertex
+ const AliExternalTrackParam *innerTPC = 0;
+ if ((innerTPC = infoRC->GetESDtrack()->GetTPCInnerParam()) != 0)
+ {
+ if ((track = new AliExternalTrackParam(*infoRC->GetESDtrack()->GetTPCInnerParam())) != 0 )
+ {
+ Bool_t bDCAStatus = track->PropagateToDCA(fVertex,field,kMaxD,dca,cov);
+
+ // Fill parametrisation histograms (only TPC track)
+ if(bDCAStatus)
+ {
+ deltaPtTPC= (mcpt-innerTPC->Pt())/mcpt;
+ pullPtTPC= (1/mcpt-innerTPC->OneOverPt())/TMath::Sqrt(innerTPC->GetSigma1Pt2());
+ deltaPhiTPC = TMath::ATan2(innerTPC->Py(),innerTPC->Px())-
+ TMath::ATan2(infoMC->GetParticle().Py(),infoMC->GetParticle().Px());
+
+ deltaTanTPC = TMath::ATan2(innerTPC->Pz(),innerTPC->Pt())-
+ TMath::ATan2(infoMC->GetParticle().Pz(),infoMC->GetParticle().Pt());
+
+ delta1Pt2TPC = (1/mcpt-innerTPC->OneOverPt())/TMath::Power(1+1/mcpt,2);
+ deltaY1PtTPC= (infoMC->GetParticle().Vy()-innerTPC->GetY()) / (0.2+1/mcpt);
+ deltaZ1PtTPC = (infoMC->GetParticle().Vz()-innerTPC->GetZ()) / (0.2+1/mcpt);
+ deltaPhi1PtTPC = deltaPhiTPC / (0.1+1/mcpt);
+ deltaTheta1PtTPC = deltaTanTPC / (0.1+1/mcpt);
+
+ f1Pt2ResolS1PtTPC->Fill(s1mcpt,delta1Pt2TPC);
+ fYResolS1PtTPC->Fill(s1mcpt,deltaY1PtTPC);
+ fZResolS1PtTPC->Fill(s1mcpt,deltaZ1PtTPC);
+ fPhiResolS1PtTPC->Fill(s1mcpt,deltaPhi1PtTPC);
+ fThetaResolS1PtTPC->Fill(s1mcpt,deltaTheta1PtTPC);
+ }
+ delete track;
+ }
+ }
+
+ // TPC and ITS (nb. of clusters >2) in the system
+ if(infoRC->GetESDtrack()->GetITSclusters(clusterITS)>2)
+ {
+ f1Pt2ResolS1PtTPCITS->Fill(s1mcpt,delta1Pt2);
+ fYResolS1PtTPCITS->Fill(s1mcpt,deltaY1Pt);
+ fZResolS1PtTPCITS->Fill(s1mcpt,deltaZ1Pt);
+ fPhiResolS1PtTPCITS->Fill(s1mcpt,deltaPhi1Pt);
+ fThetaResolS1PtTPCITS->Fill(s1mcpt,deltaTheta1Pt);
+ }
+
+ // Fill histograms
+ fPtResolLPT->Fill(mcpt,deltaPt);
+ fPtResolHPT->Fill(mcpt,deltaPt);
+ fPtPullLPT->Fill(mcpt,pullPt);
+ fPtPullHPT->Fill(mcpt,pullPt);
+
+ fPhiResolTan->Fill(tantheta,deltaPhi);
+ fPhiPullTan->Fill(tantheta,pullPhi);
+ fTanResolTan->Fill(tantheta,deltaTan);
+ fTanPullTan->Fill(tantheta,pullTan);
+
+}
+
+//_____________________________________________________________________________
+void AliComparisonRes::ProcessConstrained(AliMCInfo* infoMC, AliESDRecInfo *infoRC)
+{
+ // Fill resolution comparison information (constarained parameters)
+ //
+ AliExternalTrackParam *track = 0;
+ 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 deltaPt, pullPt, deltaPhi, pullPhi, deltaTan, pullTan, delta1Pt2, deltaY1Pt, deltaZ1Pt, deltaPhi1Pt, deltaTheta1Pt;
+ Float_t deltaPtTPC, pullPtTPC, deltaPhiTPC, deltaTanTPC, delta1Pt2TPC, deltaY1PtTPC, deltaZ1PtTPC, deltaPhi1PtTPC, deltaTheta1PtTPC;
+
+ Float_t mcpt = infoMC->GetParticle().Pt();
+ Float_t s1mcpt = TMath::Sqrt(1./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])<fCutsMC->GetMaxR() && TMath::Abs(dv[2])<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;
+
+// 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();
+ deltaPt= (mcpt-cparam->Pt())/mcpt;
+ pullPt= (1/mcpt-cparam->OneOverPt())/TMath::Sqrt(cparam->GetSigma1Pt2());
+ deltaPhi = TMath::ATan2(cparam->Py(),cparam->Px())-
+ TMath::ATan2(infoMC->GetParticle().Py(),infoMC->GetParticle().Px());
+ pullPhi = deltaPhi/TMath::Sqrt(cparam->GetSigmaSnp2());
+ deltaTan = TMath::ATan2(cparam->Pz(),cparam->Pt())-TMath::ATan2(infoMC->GetParticle().Pz(),infoMC->GetParticle().Pt());
+ pullTan = deltaPhi/TMath::Sqrt(cparam->GetSigmaSnp2());
+
+
+ delta1Pt2 = (1/mcpt-cparam->OneOverPt())/TMath::Power(1+1/mcpt,2);
+
+ deltaY1Pt = (infoMC->GetParticle().Vy()-cparam->GetY()) / (0.2+1/mcpt);
+ deltaZ1Pt = (infoMC->GetParticle().Vz()-cparam->GetZ()) / (0.2+1/mcpt);
+ deltaPhi1Pt = deltaPhi / (0.1+1/mcpt);
+ deltaTheta1Pt = deltaTan / (0.1+1/mcpt);
+
+ // calculate track parameters at vertex
+ const AliExternalTrackParam *innerTPC = 0;
+ if ((innerTPC = infoRC->GetESDtrack()->GetTPCInnerParam()) != 0)
+ {
+ if ((track = new AliExternalTrackParam(*infoRC->GetESDtrack()->GetTPCInnerParam())) != 0 )
+ {
+ Bool_t bDCAStatus = track->PropagateToDCA(fVertex,field,kMaxD,dca,cov);
+
+ // Fill parametrisation histograms (only TPC track)
+ if(bDCAStatus)
+ {
+ deltaPtTPC= (mcpt-innerTPC->Pt())/mcpt;
+ pullPtTPC= (1/mcpt-innerTPC->OneOverPt())/TMath::Sqrt(innerTPC->GetSigma1Pt2());
+ deltaPhiTPC = TMath::ATan2(innerTPC->Py(),innerTPC->Px())-
+ TMath::ATan2(infoMC->GetParticle().Py(),infoMC->GetParticle().Px());
+
+ deltaTanTPC = TMath::ATan2(innerTPC->Pz(),innerTPC->Pt())-
+ TMath::ATan2(infoMC->GetParticle().Pz(),infoMC->GetParticle().Pt());
+
+ delta1Pt2TPC = (1/mcpt-innerTPC->OneOverPt())/TMath::Power(1+1/mcpt,2);
+ deltaY1PtTPC= (infoMC->GetParticle().Vy()-innerTPC->GetY()) / (0.2+1/mcpt);
+ deltaZ1PtTPC = (infoMC->GetParticle().Vz()-innerTPC->GetZ()) / (0.2+1/mcpt);
+ deltaPhi1PtTPC = deltaPhiTPC / (0.1+1/mcpt);
+ deltaTheta1PtTPC = deltaTanTPC / (0.1+1/mcpt);
+
+ fC1Pt2ResolS1PtTPC->Fill(s1mcpt,delta1Pt2TPC);
+ fCYResolS1PtTPC->Fill(s1mcpt,deltaY1PtTPC);
+ fCZResolS1PtTPC->Fill(s1mcpt,deltaZ1PtTPC);
+ fCPhiResolS1PtTPC->Fill(s1mcpt,deltaPhi1PtTPC);
+ fCThetaResolS1PtTPC->Fill(s1mcpt,deltaTheta1PtTPC);
+ }
+ delete track;
+ }
+ }
+
+ // TPC and ITS (nb. of clusters >2) in the system
+ if(infoRC->GetESDtrack()->GetITSclusters(clusterITS)>2)
+ {
+ fC1Pt2ResolS1PtTPCITS->Fill(s1mcpt,delta1Pt2);
+ fCYResolS1PtTPCITS->Fill(s1mcpt,deltaY1Pt);
+ fCZResolS1PtTPCITS->Fill(s1mcpt,deltaZ1Pt);
+ fCPhiResolS1PtTPCITS->Fill(s1mcpt,deltaPhi1Pt);
+ fCThetaResolS1PtTPCITS->Fill(s1mcpt,deltaTheta1Pt);
+ }
+
+ // Fill histograms
+ fCPtResolTan->Fill(tantheta,deltaPt);
+ fCPtPullTan->Fill(tantheta,pullPt);
+ 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 folder "folderRes"
+ //
+
+ TH1::AddDirectory(kFALSE);
+
+ AliComparisonRes * comp=this;
+ TH1F *hiss=0;
+ TObjArray *aFolderObj = new TObjArray;
+
+ // write results in the folder
+
+ 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->SetName("CptResolTan");
+
+ aFolderObj->Add(hiss);
+
+ //
+ hiss = comp->MakeResol(comp->fPtResolLPT,1,0);
+ hiss->SetXTitle("p_{t} (GeV)");
+ hiss->SetYTitle("#sigmap_{t}/p_{t}");
+ hiss->Draw();
+ hiss->SetName("PtResolPt");
+
+ aFolderObj->Add(hiss);
+
+ //
+ hiss = comp->MakeResol(comp->fPtResolLPT,1,1);
+ hiss->SetXTitle("p_{t} (GeV)");
+ hiss->SetYTitle("mean #Deltap_{t}/p_{t}");
+ hiss->Draw();
+ hiss->SetName("PtMeanPt");
+
+ aFolderObj->Add(hiss);
+
+ //
+ hiss = comp->MakeResol(comp->fPhiResolTan,1,0);
+ hiss->SetXTitle("Tan(#theta)");
+ hiss->SetYTitle("#sigma#phi (rad)");
+ hiss->Draw();
+ hiss->SetName("PhiResolTan");
+
+ aFolderObj->Add(hiss);
+ //
+ hiss = comp->MakeResol(comp->fTanResolTan,1,0);
+ hiss->SetXTitle("Tan(#theta)");
+ hiss->SetYTitle("#sigma#theta (rad)");
+ hiss->Draw();
+ hiss->SetName("ThetaResolTan");
+
+ aFolderObj->Add(hiss);
+
+ //
+ hiss = comp->MakeResol(comp->fPhiResolTan,1,1);
+ hiss->SetXTitle("Tan(#theta)");
+ hiss->SetYTitle("mean #Delta#phi (rad)");
+ hiss->Draw();
+ hiss->SetName("PhiMeanTan");
+
+ aFolderObj->Add(hiss);
+ //
+ hiss = comp->MakeResol(comp->fTanResolTan,1,1);
+ hiss->SetXTitle("Tan(#theta)");
+ hiss->SetYTitle("mean #Delta#theta (rad)");
+ hiss->Draw();
+ hiss->SetName("ThetaMeanTan");
+
+ aFolderObj->Add(hiss);
+
+ //
+ hiss = comp->MakeResol(comp->fCPhiResolTan,1,0);
+ hiss->SetXTitle("Tan(#theta)");
+ hiss->SetYTitle("#sigma#phi (rad)");
+ hiss->Draw();
+ hiss->SetName("CPhiResolTan");
+
+ aFolderObj->Add(hiss);
+ //
+ hiss = comp->MakeResol(comp->fCTanResolTan,1,0);
+ hiss->SetXTitle("Tan(#theta)");
+ hiss->SetYTitle("#sigma#theta (rad)");
+ hiss->Draw();
+ hiss->SetName("CThetaResolTan");
+
+ aFolderObj->Add(hiss);
+ //
+ 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->SetName("CptPullTan");
+
+ aFolderObj->Add(hiss);
+ //
+ hiss = comp->MakeResol(comp->fC1Pt2ResolS1PtTPC,1,0);
+ hiss->SetXTitle("#sqrt(1/mcp_{t})");
+ hiss->SetYTitle("1/mcp_{t}-1/p_{t}/(1+1/p_{t})^2");
+ hiss->Draw();
+ hiss->SetName("C1Pt2ResolS1PtTPC");
+
+ aFolderObj->Add(hiss);
+
+ hiss = comp->MakeResol(comp->fC1Pt2ResolS1PtTPCITS,1,0);
+ hiss->SetXTitle("#sqrt(1/mcp_{t})");
+ hiss->SetYTitle("1/mcp_{t}-1/p_{t}/(1+1/p_{t})^2");
+ hiss->Draw();
+ hiss->SetName("C1Pt2ResolS1PtTPCITS");
+
+ aFolderObj->Add(hiss);
+ //
+ hiss = comp->MakeResol(comp->fCYResolS1PtTPC,1,0);
+ hiss->SetXTitle("#sqrt(1/mcp_{t})");
+ hiss->SetYTitle("(mcy-y)/(0.2+1/mcp_{t})");
+ hiss->Draw();
+ hiss->SetName("CYResolS1PtTPC");
+
+ aFolderObj->Add(hiss);
+
+ hiss = comp->MakeResol(comp->fCYResolS1PtTPCITS,1,0);
+ hiss->SetXTitle("#sqrt(1/mcp_{t})");
+ hiss->SetYTitle("(mcy-y)/(0.2+1/mcp_{t})");
+ hiss->Draw();
+ hiss->SetName("CYResolS1PtTPCITS");
+
+ aFolderObj->Add(hiss);
+ //
+ hiss = comp->MakeResol(comp->fCZResolS1PtTPC,1,0);
+ hiss->SetXTitle("#sqrt(1/mcp_{t})");
+ hiss->SetYTitle("(mcz-z)/(0.2+1/mcp_{t})");
+ hiss->Draw();
+ hiss->SetName("CZResolS1PtTPC");
+
+ aFolderObj->Add(hiss);
+
+ hiss = comp->MakeResol(comp->fCZResolS1PtTPCITS,1,0);
+ hiss->SetXTitle("#sqrt(1/mcp_{t})");
+ hiss->SetYTitle("(mcz-z)/(0.2+1/mcp_{t})");
+ hiss->Draw();
+ hiss->SetName("CZResolS1PtTPCITS");
+
+ aFolderObj->Add(hiss);
+ //
+ hiss = comp->MakeResol(comp->fCPhiResolS1PtTPC,1,0);
+ hiss->SetXTitle("#sqrt(1/mcp_{t})");
+ hiss->SetYTitle("(mc#phi-#phi)/(0.1+1/mcp_{t})");
+ hiss->Draw();
+ hiss->SetName("CPhiResolS1PtTPC");
+
+ aFolderObj->Add(hiss);
+
+ hiss = comp->MakeResol(comp->fCPhiResolS1PtTPCITS,1,0);
+ hiss->SetXTitle("#sqrt(1/mcp_{t})");
+ hiss->SetYTitle("(mc#phi-#phi)/(0.1+1/mcp_{t})");
+ hiss->Draw();
+ hiss->SetName("CPhiResolS1PtTPCITS");
+
+ aFolderObj->Add(hiss);
+ //
+ hiss = comp->MakeResol(comp->fCThetaResolS1PtTPC,1,0);
+ hiss->SetXTitle("#sqrt(1/mcp_{t})");
+ hiss->SetYTitle("(mc#theta-#theta)/(0.1+1/mcp_{t})");
+ hiss->Draw();
+ hiss->SetName("CThetaResolS1PtTPC");
+
+ aFolderObj->Add(hiss);
+
+ hiss = comp->MakeResol(comp->fCThetaResolS1PtTPCITS,1,0);
+ hiss->SetXTitle("#sqrt(1/mcp_{t})");
+ hiss->SetYTitle("(mc#theta-#theta)/(0.1+1/mcp_{t})");
+ hiss->Draw();
+ hiss->SetName("CThetaResolS1PtTPCITS");
+
+ aFolderObj->Add(hiss);
+
+ //
+ hiss = comp->MakeResol(comp->f1Pt2ResolS1PtTPC,1,0);
+ hiss->SetXTitle("#sqrt(1/mcp_{t})");
+ hiss->SetYTitle("1/mcp_{t}-1/p_{t}/(1+1/p_{t})^2");
+ hiss->Draw();
+ hiss->SetName("OnePt2ResolS1PtTPC");
+
+ aFolderObj->Add(hiss);
+
+ hiss = comp->MakeResol(comp->f1Pt2ResolS1PtTPCITS,1,0);
+ hiss->SetXTitle("#sqrt(1/mcp_{t})");
+ hiss->SetYTitle("1/mcp_{t}-1/p_{t}/(1+1/p_{t})^2");
+ hiss->Draw();
+ hiss->SetName("OnePt2ResolS1PtTPCITS");
+
+ aFolderObj->Add(hiss);
+ //
+ hiss = comp->MakeResol(comp->fYResolS1PtTPC,1,0);
+ hiss->SetXTitle("#sqrt(1/mcp_{t})");
+ hiss->SetYTitle("(mcy-y)/(0.2+1/mcp_{t})");
+ hiss->Draw();
+ hiss->SetName("YResolS1PtTPC");
+
+ aFolderObj->Add(hiss);
+
+ hiss = comp->MakeResol(comp->fYResolS1PtTPCITS,1,0);
+ hiss->SetXTitle("#sqrt(1/mcp_{t})");
+ hiss->SetYTitle("(mcy-y)/(0.2+1/mcp_{t})");
+ hiss->Draw();
+ hiss->SetName("YResolS1PtTPCITS");
+
+ aFolderObj->Add(hiss);
+ //
+ hiss = comp->MakeResol(comp->fZResolS1PtTPC,1,0);
+ hiss->SetXTitle("#sqrt(1/mcp_{t})");
+ hiss->SetYTitle("(mcz-z)/(0.2+1/mcp_{t})");
+ hiss->Draw();
+ hiss->SetName("ZResolS1PtTPC");
+
+ aFolderObj->Add(hiss);
+
+ hiss = comp->MakeResol(comp->fZResolS1PtTPCITS,1,0);
+ hiss->SetXTitle("#sqrt(1/mcp_{t})");
+ hiss->SetYTitle("(mcz-z)/(0.2+1/mcp_{t})");
+ hiss->Draw();
+ hiss->SetName("ZResolS1PtTPCITS");
+
+ aFolderObj->Add(hiss);
+ //
+ hiss = comp->MakeResol(comp->fPhiResolS1PtTPC,1,0);
+ hiss->SetXTitle("#sqrt(1/mcp_{t})");
+ hiss->SetYTitle("(mc#phi-#phi)/(0.1+1/mcp_{t})");
+ hiss->Draw();
+ hiss->SetName("PhiResolS1PtTPC");
+
+ aFolderObj->Add(hiss);
+
+ hiss = comp->MakeResol(comp->fPhiResolS1PtTPCITS,1,0);
+ hiss->SetXTitle("#sqrt(1/mcp_{t})");
+ hiss->SetYTitle("(mc#phi-#phi)/(0.1+1/mcp_{t})");
+ hiss->Draw();
+ hiss->SetName("PhiResolS1PtTPCITS");
+
+ aFolderObj->Add(hiss);
+ //
+ hiss = comp->MakeResol(comp->fThetaResolS1PtTPC,1,0);
+ hiss->SetXTitle("#sqrt(1/mcp_{t})");
+ hiss->SetYTitle("(mc#theta-#theta)/(0.1+1/mcp_{t})");
+ hiss->Draw();
+ hiss->SetName("ThetaResolS1PtTPC");
+
+ aFolderObj->Add(hiss);
+
+ hiss = comp->MakeResol(comp->fThetaResolS1PtTPCITS,1,0);
+ hiss->SetXTitle("#sqrt(1/mcp_{t})");
+ hiss->SetYTitle("(mc#theta-#theta)/(0.1+1/mcp_{t})");
+ hiss->Draw();
+ hiss->SetName("ThetaResolS1PtTPCITS");
+
+ aFolderObj->Add(hiss);
+
+ // export objects to analysis folder
+ fAnalysisFolder = ExportToFolder(aFolderObj);
+
+ // delete only TObjArray
+ if(aFolderObj) delete aFolderObj;
+}
+
+//_____________________________________________________________________________
+TFolder* AliComparisonRes::ExportToFolder(TObjArray * array)
+{
+ // recreate folder avery time and export objects to new one
+ //
+ AliComparisonRes * comp=this;
+ TFolder *folder = comp->GetAnalysisFolder();
+
+ TString name, title;
+ TFolder *newFolder = 0;
+ Int_t i = 0;
+ Int_t size = array->GetSize();
+
+ if(folder) {
+ // get name and title from old folder
+ name = folder->GetName();
+ title = folder->GetTitle();
+
+ // delete old one
+ delete folder;
+
+ // create new one
+ newFolder = CreateFolder(name.Data(),title.Data());
+ newFolder->SetOwner();
+
+ // add objects to folder
+ while(i < size) {
+ newFolder->Add(array->At(i));
+ i++;
+ }
+ }
+
+return newFolder;
+}
+
+//_____________________________________________________________________________
+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) continue;
+
+ fPtResolLPT->Add(entry->fPtResolLPT);
+ fPtResolHPT->Add(entry->fPtResolHPT);
+ fPtPullLPT->Add(entry->fPtPullLPT);
+ fPtPullHPT->Add(entry->fPtPullHPT);
+ fPhiResolTan->Add(entry->fPhiResolTan);
+ fTanResolTan->Add(entry->fTanResolTan);
+ fPhiPullTan->Add(entry->fPhiPullTan);
+ fTanPullTan->Add(entry->fTanPullTan);
+
+ // Histograms for 1/pt parameterisation
+ f1Pt2ResolS1PtTPC->Add(entry->f1Pt2ResolS1PtTPC);
+ fYResolS1PtTPC->Add(entry->fYResolS1PtTPC);
+ fZResolS1PtTPC->Add(entry->fZResolS1PtTPC);
+ fPhiResolS1PtTPC->Add(entry->fPhiResolS1PtTPC);
+ fThetaResolS1PtTPC->Add(entry->fThetaResolS1PtTPC);
+
+ f1Pt2ResolS1PtTPCITS->Add(entry->f1Pt2ResolS1PtTPCITS);
+ fYResolS1PtTPCITS->Add(entry->fYResolS1PtTPCITS);
+ fZResolS1PtTPCITS->Add(entry->fZResolS1PtTPCITS);
+ fPhiResolS1PtTPCITS->Add(entry->fPhiResolS1PtTPCITS);
+ fThetaResolS1PtTPCITS->Add(entry->fThetaResolS1PtTPCITS);
+
+ // 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 (constrained)
+ fC1Pt2ResolS1PtTPC->Add(entry->fC1Pt2ResolS1PtTPC);
+ fCYResolS1PtTPC->Add(entry->fCYResolS1PtTPC);
+ fCZResolS1PtTPC->Add(entry->fCZResolS1PtTPC);
+ fCPhiResolS1PtTPC->Add(entry->fCPhiResolS1PtTPC);
+ fCThetaResolS1PtTPC->Add(entry->fCThetaResolS1PtTPC);
+
+ fC1Pt2ResolS1PtTPCITS->Add(entry->fC1Pt2ResolS1PtTPCITS);
+ fCYResolS1PtTPCITS->Add(entry->fCYResolS1PtTPCITS);
+ fCZResolS1PtTPCITS->Add(entry->fCZResolS1PtTPCITS);
+ fCPhiResolS1PtTPCITS->Add(entry->fCPhiResolS1PtTPCITS);
+ fCThetaResolS1PtTPCITS->Add(entry->fCThetaResolS1PtTPCITS);
+
+ count++;
+ }
+
+return count;
+}
+
+//_____________________________________________________________________________
+TFolder* AliComparisonRes::CreateFolder(TString name,TString title) {
+// create folder for analysed histograms
+//
+TFolder *folder = 0;
+ folder = new TFolder(name.Data(),title.Data());
+
+ return folder;
+}