PWG1/AliComparisonRes.cxx

   1 //------------------------------------------------------------------------------
   2 // Implementation of AliComparisonRes class. It keeps information from
   3 // comparison of reconstructed and MC particle tracks. In addtion,
   4 // it keeps selection cuts used during comparison. The comparison
   5 // information is stored in the ROOT histograms. Analysis of these
   6 // histograms can be done by using Analyse() class function. The result of
   7 // the analysis (histograms) are stored in the output picture_res.root file.
   8 //
   9 // Author: J.Otwinowski 04/02/2008
  10 //------------------------------------------------------------------------------
  11
  12 /*
  13   //after running analysis, read the file, and get component
  14   gSystem->Load("libPWG1.so");
  15   TFile f("Output.root");
  16   AliComparisonRes * comp = (AliComparisonRes*)f.Get("AliComparisonRes");
  17
  18   // analyse comparison data (output stored in pictures_res.root)
  19   comp->Analyse();
  20
  21   // paramtetrisation of the TPC track length (for information only)
  22   TF1 fl("fl","((min(250./(abs(x+0.000001)),250)-90))",0,2);  // TPC track length function
  23   TF1 fl2("fl2","[0]/((min(250./(abs(x+0.000001)),250)-90))^[1]",0,2);
  24   fl2.SetParameter(1,1);
  25   fl2.SetParameter(0,1);
  26 */
  27
  28 #include <iostream>
  29
  30 #include "TFile.h"
  31 #include "TCint.h"
  32 #include "TH3F.h"
  33 #include "TH2F.h"
  34 #include "TF1.h"
  35 #include "TProfile.h"
  36 #include "TProfile2D.h"
  37 #include "TGraph2D.h"
  38 #include "TCanvas.h"
  39 #include "TGraph.h"
  40
  41 #include "AliESDEvent.h"
  42 #include "AliESD.h"
  43 #include "AliESDfriend.h"
  44 #include "AliESDfriendTrack.h"
  45 #include "AliESDVertex.h"
  46 #include "AliRecInfoCuts.h"
  47 #include "AliMCInfoCuts.h"
  48 #include "AliLog.h"
  49 #include "AliTracker.h"
  50
  51 #include "AliMathBase.h"
  52 #include "AliTreeDraw.h"
  53
  54 #include "AliMCInfo.h"
  55 #include "AliESDRecInfo.h"
  56 #include "AliComparisonRes.h"
  57
  58 using namespace std;
  59
  60 ClassImp(AliComparisonRes)
  61
  62 //_____________________________________________________________________________
  63 AliComparisonRes::AliComparisonRes():
  64   TNamed("AliComparisonRes","AliComparisonRes"),
  65
  66   // Resolution
  67   fPtResolLPT(0),        // pt resolution - low pt
  68   fPtResolHPT(0),        // pt resolution - high pt
  69   fPtPullLPT(0),         // pt resolution - low pt
  70   fPtPullHPT(0),         // pt resolution - high pt
  71   //
  72   // Resolution constrained param
  73   //
  74   fCPhiResolTan(0),  // angular resolution -  constrained
  75   fCTanResolTan(0),  // angular resolution -  constrained
  76   fCPtResolTan(0),   // pt resolution      -  constrained
  77   fCPhiPullTan(0),   // angular resolution -  constrained
  78   fCTanPullTan(0),   // angular resolution -  constrained
  79   fCPtPullTan(0),    // pt resolution      -  constrained
  80
  81   //
  82   // Parametrisation histograms
  83   //
  84
  85   f1Pt2Resol1PtTPC(0),
  86   f1Pt2Resol1PtTPCITS(0),
  87   fYResol1PtTPC(0),
  88   fYResol1PtTPCITS(0),
  89   fZResol1PtTPC(0),
  90   fZResol1PtTPCITS(0),
  91   fPhiResol1PtTPC(0),
  92   fPhiResol1PtTPCITS(0),
  93   fThetaResol1PtTPC(0),
  94   fThetaResol1PtTPCITS(0),
  95
  96   // constrained
  97   fC1Pt2Resol1PtTPC(0),
  98   fC1Pt2Resol1PtTPCITS(0),
  99   fCYResol1PtTPC(0),
 100   fCYResol1PtTPCITS(0),
 101   fCZResol1PtTPC(0),
 102   fCZResol1PtTPCITS(0),
 103   fCPhiResol1PtTPC(0),
 104   fCPhiResol1PtTPCITS(0),
 105   fCThetaResol1PtTPC(0),
 106   fCThetaResol1PtTPCITS(0),
 107
 108   // vertex
 109   fVertex(0),
 110
 111   // Cuts
 112   fCutsRC(0),
 113   fCutsMC(0)
 114 {
 115   InitHisto();
 116   InitCuts();
 117
 118   // vertex (0,0,0)
 119   fVertex = new AliESDVertex();
 120   fVertex->SetXv(0.0);
 121   fVertex->SetYv(0.0);
 122   fVertex->SetZv(0.0);
 123 }
 124
 125 //_____________________________________________________________________________
 126 AliComparisonRes::~AliComparisonRes(){
 127
 128   // Resolution histograms
 129   if(fPtResolLPT) delete  fPtResolLPT; fPtResolLPT=0;
 130   if(fPtResolHPT) delete  fPtResolHPT; fPtResolHPT=0;
 131   if(fPtPullLPT)  delete  fPtPullLPT;  fPtPullLPT=0;
 132   if(fPtPullHPT)  delete  fPtPullHPT;  fPtPullHPT=0;
 133
 134   // Resolution histograms (constrained param)
 135   if(fCPhiResolTan) delete fCPhiResolTan; fCPhiResolTan=0;
 136   if(fCTanResolTan) delete fCTanResolTan; fCTanResolTan=0;
 137   if(fCPtResolTan)  delete fCPtResolTan;  fCPtResolTan=0;
 138   if(fCPhiPullTan)  delete fCPhiPullTan;  fCPhiPullTan=0;
 139   if(fCTanPullTan)  delete fCTanPullTan;  fCTanPullTan=0;
 140   if(fCPtPullTan)   delete fCPtPullTan;   fCPtPullTan=0;
 141
 142   // Parametrisation histograms
 143   //
 144   if(f1Pt2Resol1PtTPC) delete f1Pt2Resol1PtTPC; f1Pt2Resol1PtTPC=0;
 145   if(f1Pt2Resol1PtTPCITS) delete f1Pt2Resol1PtTPCITS; f1Pt2Resol1PtTPCITS=0;
 146   if(fYResol1PtTPC) delete fYResol1PtTPC; fYResol1PtTPC=0;
 147   if(fYResol1PtTPCITS) delete fYResol1PtTPCITS; fYResol1PtTPCITS=0;
 148   if(fZResol1PtTPC) delete fZResol1PtTPC; fZResol1PtTPC=0;
 149   if(fZResol1PtTPCITS) delete fZResol1PtTPCITS; fZResol1PtTPCITS=0;
 150   if(fPhiResol1PtTPC) delete fPhiResol1PtTPC; fPhiResol1PtTPC=0;
 151   if(fPhiResol1PtTPCITS) delete fPhiResol1PtTPCITS; fPhiResol1PtTPCITS=0;
 152   if(fThetaResol1PtTPC) delete fThetaResol1PtTPC; fThetaResol1PtTPC=0;
 153   if(fThetaResol1PtTPCITS) delete fThetaResol1PtTPCITS; fThetaResol1PtTPCITS=0;
 154
 155   // constrained
 156   if(fC1Pt2Resol1PtTPC) delete fC1Pt2Resol1PtTPC; fC1Pt2Resol1PtTPC=0;
 157   if(fC1Pt2Resol1PtTPCITS) delete fC1Pt2Resol1PtTPCITS; fC1Pt2Resol1PtTPCITS=0;
 158   if(fCYResol1PtTPC) delete fCYResol1PtTPC; fCYResol1PtTPC=0;
 159   if(fCYResol1PtTPCITS) delete fCYResol1PtTPCITS; fCYResol1PtTPCITS=0;
 160   if(fCZResol1PtTPC) delete fCZResol1PtTPC; fCZResol1PtTPC=0;
 161   if(fCZResol1PtTPCITS) delete fCZResol1PtTPCITS; fCZResol1PtTPCITS=0;
 162   if(fCPhiResol1PtTPC) delete fCPhiResol1PtTPC; fCPhiResol1PtTPC=0;
 163   if(fCPhiResol1PtTPCITS) delete fCPhiResol1PtTPCITS; fCPhiResol1PtTPCITS=0;
 164   if(fCThetaResol1PtTPC) delete fCThetaResol1PtTPC; fCThetaResol1PtTPC=0;
 165   if(fCThetaResol1PtTPCITS) delete fCThetaResol1PtTPCITS; fCThetaResol1PtTPCITS=0;
 166
 167   if(fVertex) delete fVertex; fVertex=0;
 168
 169 }
 170
 171 //_____________________________________________________________________________
 172 void AliComparisonRes::InitHisto(){
 173
 174   // Init histograms
 175   fCPhiResolTan = new TH2F("CPhiResolTan","CPhiResolTan",50, -2,2,200,-0.025,0.025);
 176   fCPhiResolTan->SetXTitle("tan(#theta)");
 177   fCPhiResolTan->SetYTitle("#Delta#phi");
 178
 179   fCTanResolTan = new TH2F("CTanResolTan","CTanResolTan",50, -2,2,200,-0.025,0.025);
 180   fCTanResolTan->SetXTitle("tan(#theta)");
 181   fCTanResolTan->SetYTitle("#Delta#theta");
 182
 183   fCPtResolTan=new TH2F("CPtResol","CPtResol",50, -2,2,200,-0.2,0.2);
 184   fCPtResolTan->SetXTitle("Tan(#theta)");
 185   fCPtResolTan->SetYTitle("#Deltap_{t}/p_{t}");
 186
 187   fCPhiPullTan = new TH2F("CPhiPullTan","CPhiPullTan",50, -2,2,200,-5,5);
 188   fCPhiPullTan->SetXTitle("Tan(#theta)");
 189   fCPhiPullTan->SetYTitle("#Delta#phi/#Sigma");
 190
 191   fCTanPullTan = new TH2F("CTanPullTan","CTanPullTan",50, -2,2,200,-5,5);
 192   fCTanPullTan->SetXTitle("Tan(#theta)");
 193   fCTanPullTan->SetYTitle("#Delta#theta/#Sigma");
 194
 195   fCPtPullTan=new TH2F("CPtPull","CPtPull",50, -2,2,200,-5,5);
 196   fCPtPullTan->SetXTitle("Tan(#theta)");
 197   fCPtPullTan->SetYTitle("(1/mcp_{t}-1/p_{t})/#Sigma");
 198
 199   fPtResolLPT = new TH2F("Pt_resol_lpt","pt resol",10, 0.1,3,200,-0.2,0.2);
 200   fPtResolLPT->SetXTitle("p_{t}");
 201   fPtResolLPT->SetYTitle("#Deltap_{t}/p_{t}");
 202
 203   fPtResolHPT = new TH2F("Pt_resol_hpt","pt resol",10, 2,100,200,-0.3,0.3);
 204   fPtResolHPT->SetXTitle("p_{t}");
 205   fPtResolHPT->SetYTitle("#Deltap_{t}/p_{t}");
 206
 207   fPtPullLPT = new TH2F("Pt_pull_lpt","pt pull",10, 0.1,3,200,-6,6);
 208   fPtPullLPT->SetXTitle("p_{t}");
 209   fPtPullLPT->SetYTitle("#Deltap_{t}/#Sigma");
 210
 211   fPtPullHPT = new TH2F("Pt_pull_hpt","pt pull",10,2,100,200,-6,6);
 212   fPtPullHPT->SetXTitle("p_{t}");
 213   fPtPullHPT->SetYTitle("#Deltap_{t}/#Sigma");
 214
 215   //
 216   // Parametrisation histograms
 217   //
 218
 219   f1Pt2Resol1PtTPC = new TH2F("f1Pt2Resol1PtTPC","(1/mcpt-1/pt)/(1+1/mcpt)^2 vs 1/pt)",100,0,10,200,-0.010,0.010);
 220   f1Pt2Resol1PtTPC->SetXTitle("1/mcp_{t}");
 221   f1Pt2Resol1PtTPC->SetYTitle("(1/mcp_{t}-1/p_{t})/(1+1/mcp_{t})^2)");
 222
 223   f1Pt2Resol1PtTPCITS = new TH2F("f1Pt2Resol1PtTPCITS","(1/mcpt-1/pt)/(1+1/mcpt)^2 vs 1/pt)",100,0,10,200,-0.010,0.010);
 224   f1Pt2Resol1PtTPCITS->SetXTitle("1/mcp_{t}");
 225   f1Pt2Resol1PtTPCITS->SetYTitle("(1/mcp_{t}-1/p_{t})/(1+1/mcp_{t})^2)");
 226
 227   fYResol1PtTPC = new TH2F("fYResol1PtTPC","fYResol1PtTPC",100, 0,10,200,-1.0,1.0);
 228   fYResol1PtTPC->SetXTitle("1/mcpt");
 229   fYResol1PtTPC->SetYTitle("#DeltaY");
 230
 231   fYResol1PtTPCITS = new TH2F("fYResol1PtTPCITS","fYResol1PtTPCITS",100, 0,10,200,-0.05,0.05);
 232   fYResol1PtTPCITS->SetXTitle("1/mcpt");
 233   fYResol1PtTPCITS->SetYTitle("#DeltaY");
 234
 235   fZResol1PtTPC = new TH2F("fZResol1PtTPC","fZResol1PtTPC",100, 0,10,200,-1.0,1.0);
 236   fZResol1PtTPC->SetXTitle("1/mcpt");
 237   fZResol1PtTPC->SetYTitle("#DeltaZ");
 238
 239   fZResol1PtTPCITS = new TH2F("fZResol1PtTPCITS","fZResol1PtTPCITS",100, 0,10,200,-0.05,0.05);
 240   fZResol1PtTPCITS->SetXTitle("1/mcpt");
 241   fZResol1PtTPCITS->SetYTitle("#DeltaZ");
 242
 243   fPhiResol1PtTPC = new TH2F("fPhiResol1PtTPC","fPhiResol1PtTPC",100, 0,10,200,-0.025,0.025);
 244   fPhiResol1PtTPC->SetXTitle("1/mcpt");
 245   fPhiResol1PtTPC->SetYTitle("#Delta#phi");
 246
 247   fPhiResol1PtTPCITS = new TH2F("fPhiResol1PtTPCITS","fPhiResol1PtTPCITS",100, 0,10,200,-0.01,0.01);
 248   fPhiResol1PtTPCITS->SetXTitle("1/mcpt");
 249   fPhiResol1PtTPCITS->SetYTitle("#Delta#phi");
 250
 251   fThetaResol1PtTPC = new TH2F("fThetaResol1PtTPC","fThetaResol1PtTPC",100, 0,10,200,-0.025,0.025);
 252   fThetaResol1PtTPC->SetXTitle("1/mcpt");
 253   fThetaResol1PtTPC->SetYTitle("#Delta#theta");
 254
 255   fThetaResol1PtTPCITS = new TH2F("fThetaResol1PtTPCITS","fThetaResol1PtTPCITS",100, 0,10,200,-0.01,0.01);
 256   fThetaResol1PtTPCITS->SetXTitle("1/mcpt");
 257   fThetaResol1PtTPCITS->SetYTitle("#Delta#theta");
 258
 259   // constrained
 260   fC1Pt2Resol1PtTPC = new TH2F("fC1Pt2Resol1PtTPC","(1/mcpt-1/pt)/(1+1/mcpt)^2 vs 1/pt)",100,0,10,200,-0.010,0.010);
 261   fC1Pt2Resol1PtTPC->SetXTitle("1/mcp_{t}");
 262   fC1Pt2Resol1PtTPC->SetYTitle("(1/mcp_{t}-1/p_{t})/(1+1/mcp_{t})^2)");
 263
 264   fC1Pt2Resol1PtTPCITS = new TH2F("fC1Pt2Resol1PtTPCITS","(1/mcpt-1/pt)/(1+1/mcpt)^2 vs 1/pt)",100,0,10,200,-0.010,0.010);
 265   fC1Pt2Resol1PtTPCITS->SetXTitle("1/mcp_{t}");
 266   fC1Pt2Resol1PtTPCITS->SetYTitle("(1/mcp_{t}-1/p_{t})/(1+1/mcp_{t})^2)");
 267
 268   fCYResol1PtTPC = new TH2F("fCYResol1PtTPC","fCYResol1PtTPC",100, 0,10,200,-1.0,1.0);
 269   fCYResol1PtTPC->SetXTitle("1/mcpt");
 270   fCYResol1PtTPC->SetYTitle("#DeltaY");
 271
 272   fCYResol1PtTPCITS = new TH2F("fCYResol1PtTPCITS","fCYResol1PtTPCITS",100, 0,10,200,-0.05,0.05);
 273   fCYResol1PtTPCITS->SetXTitle("1/mcpt");
 274   fCYResol1PtTPCITS->SetYTitle("#DeltaY");
 275
 276   fCZResol1PtTPC = new TH2F("fCZResol1PtTPC","fCZResol1PtTPC",100, 0,10,200,-1.0,1.0);
 277   fCZResol1PtTPC->SetXTitle("1/mcpt");
 278   fCZResol1PtTPC->SetYTitle("#DeltaZ");
 279
 280   fCZResol1PtTPCITS = new TH2F("fCZResol1PtTPCITS","fCZResol1PtTPCITS",100, 0,10,200,-0.05,0.05);
 281   fCZResol1PtTPCITS->SetXTitle("1/mcpt");
 282   fCZResol1PtTPCITS->SetYTitle("#DeltaZ");
 283
 284   fCPhiResol1PtTPC = new TH2F("fCPhiResol1PtTPC","fCPhiResol1PtTPC",100, 0,10,200,-0.025,0.025);
 285   fCPhiResol1PtTPC->SetXTitle("1/mcpt");
 286   fCPhiResol1PtTPC->SetYTitle("#Delta#phi");
 287
 288   fCPhiResol1PtTPCITS = new TH2F("fCPhiResol1PtTPCITS","fCPhiResol1PtTPCITS",100, 0,10,200,-0.01,0.01);
 289   fCPhiResol1PtTPCITS->SetXTitle("1/mcpt");
 290   fCPhiResol1PtTPCITS->SetYTitle("#Delta#phi");
 291
 292   fCThetaResol1PtTPC = new TH2F("fCThetaResol1PtTPC","fCThetaResol1PtTPC",100, 0,10,200,-0.025,0.025);
 293   fCThetaResol1PtTPC->SetXTitle("1/mcpt");
 294   fCThetaResol1PtTPC->SetYTitle("#Delta#theta");
 295
 296   fCThetaResol1PtTPCITS = new TH2F("fCThetaResol1PtTPCITS","fCThetaResol1PtTPCITS",100, 0,10,200,-0.01,0.01);
 297   fCThetaResol1PtTPCITS->SetXTitle("1/mcpt");
 298   fCThetaResol1PtTPCITS->SetYTitle("#Delta#theta");
 299 }
 300
 301 //_____________________________________________________________________________
 302 void AliComparisonRes::InitCuts()
 303 {
 304   // Init cuts
 305   if(!fCutsMC)
 306     AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
 307   if(!fCutsRC)
 308     AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
 309 }
 310
 311 //_____________________________________________________________________________
 312 void AliComparisonRes::Process(AliMCInfo* infoMC, AliESDRecInfo *infoRC)
 313 {
 314   // Fill resolution comparison information
 315   AliExternalTrackParam *track = 0;
 316   Double_t kRadius    = 3.0;      // beam pipe radius
 317   Double_t kMaxStep   = 5.0;      // max step
 318   Double_t field      = AliTracker::GetBz(); // nominal Bz field [kG]
 319   Double_t kMaxD      = 123456.0; // max distance
 320
 321   Int_t clusterITS[200];
 322   Double_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
 323
 324   Float_t deltaPt, pullPt, deltaPhi, deltaTan, delta1Pt2, deltaY1Pt, deltaZ1Pt, deltaPhi1Pt, deltaTheta1Pt;
 325   Float_t deltaPtTPC, pullPtTPC, deltaPhiTPC, deltaTanTPC, delta1Pt2TPC, deltaY1PtTPC, deltaZ1PtTPC, deltaPhi1PtTPC, deltaTheta1PtTPC;
 326
 327   Float_t mcpt = infoMC->GetParticle().Pt();
 328
 329   // distance to Prim. vertex
 330   const Double_t* dv = infoMC->GetVDist();
 331   Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])<fCutsMC->GetMaxR() && TMath::Abs(dv[2])<fCutsMC->GetMaxVz();
 332
 333   // Check selection cuts
 334   if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return;
 335   if (!isPrim) return;
 336   //if (infoRC->GetStatus(1)==0) return;
 337   if (infoRC->GetStatus(1)!=3) return; // TPC refit
 338   if (!infoRC->GetESDtrack()) return;
 339   if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;
 340
 341   // calculate and set prim. vertex
 342   fVertex->SetXv( infoMC->GetParticle().Vx() - dv[0] );
 343   fVertex->SetYv( infoMC->GetParticle().Vy() - dv[1] );
 344   fVertex->SetZv( infoMC->GetParticle().Vz() - dv[2] );
 345
 346   deltaPt= (mcpt-infoRC->GetESDtrack()->Pt())/mcpt;
 347   pullPt= (1/mcpt-infoRC->GetESDtrack()->OneOverPt())/TMath::Sqrt(infoRC->GetESDtrack()->GetSigma1Pt2());
 348   deltaPhi = TMath::ATan2(infoRC->GetESDtrack()->Py(),infoRC->GetESDtrack()->Px())-
 349                      TMath::ATan2(infoMC->GetParticle().Py(),infoMC->GetParticle().Px());
 350
 351   deltaTan = TMath::ATan2(infoRC->GetESDtrack()->Pz(),infoRC->GetESDtrack()->Pt())-
 352                      TMath::ATan2(infoMC->GetParticle().Pz(),infoMC->GetParticle().Pt());
 353
 354   delta1Pt2 = (1/mcpt-infoRC->GetESDtrack()->OneOverPt())/TMath::Power(1+1/mcpt,2);
 355   deltaY1Pt = (infoMC->GetParticle().Vy()-infoRC->GetESDtrack()->GetY()) / (0.2+1/mcpt);
 356   deltaZ1Pt = (infoMC->GetParticle().Vz()-infoRC->GetESDtrack()->GetZ()) / (0.2+1/mcpt);
 357   deltaPhi1Pt = deltaPhi   / (0.1+1/mcpt);
 358   deltaTheta1Pt = deltaTan / (0.1+1/mcpt);
 359
 360   //track parameters at the first measured point (TPC)
 361   //Double_t param[5],x,alpha; // [0]-Y [cm],[1]-Z [cm],[2]-sin(phi),[3]-tan(theta),[4]-1/pt [1/GeV]
 362   //infoRC->GetESDtrack()->GetInnerExternalParameters(alpha,x,param);
 363   //const AliExternalTrackParam *innerTPC =  infoRC->GetESDtrack()->GetInnerParam();
 364   //const AliExternalTrackParam *innerTPC =  infoRC->GetESDtrack()->GetTPCInnerParam();
 365
 366
 367   // calculate track parameters at vertex
 368   const AliExternalTrackParam *innerTPC =  0;
 369   if ((innerTPC = infoRC->GetESDtrack()->GetTPCInnerParam()) != 0)
 370   {
 371     if ((track = new AliExternalTrackParam(*infoRC->GetESDtrack()->GetTPCInnerParam())) != 0 )
 372     {
 373       Bool_t bStatus = AliTracker::PropagateTrackTo(track,kRadius,infoMC->GetMass(),kMaxStep,kTRUE);
 374       Bool_t bDCAStatus = track->PropagateToDCA(fVertex,field,kMaxD,dca,cov);
 375
 376       // Fill parametrisation histograms (only TPC track)
 377       if(bStatus && bDCAStatus)
 378           {
 379                         deltaPtTPC= (mcpt-innerTPC->Pt())/mcpt;
 380                         pullPtTPC= (1/mcpt-innerTPC->OneOverPt())/TMath::Sqrt(innerTPC->GetSigma1Pt2());
 381                         deltaPhiTPC = TMath::ATan2(innerTPC->Py(),innerTPC->Px())-
 382                                                                 TMath::ATan2(infoMC->GetParticle().Py(),infoMC->GetParticle().Px());
 383
 384                         deltaTanTPC = TMath::ATan2(innerTPC->Pz(),innerTPC->Pt())-
 385                                                                 TMath::ATan2(infoMC->GetParticle().Pz(),infoMC->GetParticle().Pt());
 386
 387                         delta1Pt2TPC = (1/mcpt-innerTPC->OneOverPt())/TMath::Power(1+1/mcpt,2);
 388                         deltaY1PtTPC= (infoMC->GetParticle().Vy()-innerTPC->GetY()) / (0.2+1/mcpt);
 389                         deltaZ1PtTPC = (infoMC->GetParticle().Vz()-innerTPC->GetZ()) / (0.2+1/mcpt);
 390                         deltaPhi1PtTPC = deltaPhiTPC   / (0.1+1/mcpt);
 391                         deltaTheta1PtTPC = deltaTanTPC / (0.1+1/mcpt);
 392
 393                         f1Pt2Resol1PtTPC->Fill(1/mcpt,delta1Pt2TPC);
 394                         fYResol1PtTPC->Fill(1/mcpt,deltaY1PtTPC);
 395                         fZResol1PtTPC->Fill(1/mcpt,deltaZ1PtTPC);
 396                         fPhiResol1PtTPC->Fill(1/mcpt,deltaPhi1PtTPC);
 397                         fThetaResol1PtTPC->Fill(1/mcpt,deltaTheta1PtTPC);
 398           }
 399           delete track;
 400     }
 401   }
 402
 403   // TPC and ITS (nb. of clusters >2) in the system
 404   if(infoRC->GetESDtrack()->GetITSclusters(clusterITS)>2)
 405   {
 406       f1Pt2Resol1PtTPCITS->Fill(1/mcpt,delta1Pt2);
 407       fYResol1PtTPCITS->Fill(1/mcpt,deltaY1Pt);
 408       fZResol1PtTPCITS->Fill(1/mcpt,deltaZ1Pt);
 409       fPhiResol1PtTPCITS->Fill(1/mcpt,deltaPhi1Pt);
 410       fThetaResol1PtTPCITS->Fill(1/mcpt,deltaTheta1Pt);
 411   }
 412
 413   // Fill histograms
 414   fPtResolLPT->Fill(mcpt,deltaPt);
 415   fPtResolHPT->Fill(mcpt,deltaPt);
 416   fPtPullLPT->Fill(mcpt,pullPt);
 417   fPtPullHPT->Fill(mcpt,pullPt);
 418 }
 419
 420 //_____________________________________________________________________________
 421 void AliComparisonRes::ProcessConstrained(AliMCInfo* infoMC, AliESDRecInfo *infoRC)
 422 {
 423   // Fill resolution comparison information (constarained parameters)
 424   //
 425   AliExternalTrackParam *track = 0;
 426   Double_t kRadius    = 3.0;      // beam pipe radius
 427   Double_t kMaxStep   = 5.0;      // max step
 428   Double_t field      = AliTracker::GetBz(); // nominal Bz field [kG]
 429   Double_t kMaxD      = 123456.0; // max distance
 430
 431   Int_t clusterITS[200];
 432   Double_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
 433
 434   Float_t deltaPt, pullPt, deltaPhi, pullPhi, deltaTan, pullTan, delta1Pt2, deltaY1Pt, deltaZ1Pt, deltaPhi1Pt, deltaTheta1Pt;
 435   Float_t deltaPtTPC, pullPtTPC, deltaPhiTPC, deltaTanTPC, delta1Pt2TPC, deltaY1PtTPC, deltaZ1PtTPC, deltaPhi1PtTPC, deltaTheta1PtTPC;
 436
 437   Float_t mcpt = infoMC->GetParticle().Pt();
 438   Float_t tantheta = TMath::Tan(infoMC->GetParticle().Theta()-TMath::Pi()*0.5);
 439
 440   // distance to Prim. vertex
 441   const Double_t* dv = infoMC->GetVDist();
 442   Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])<fCutsMC->GetMaxR() && TMath::Abs(dv[2])<fCutsMC->GetMaxVz();
 443
 444   // Check selection cuts
 445   if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return;
 446   if (!isPrim) return;
 447   if (infoRC->GetStatus(1)!=3) return;
 448   if (!infoRC->GetESDtrack()) return;
 449   if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;
 450   if (!infoRC->GetESDtrack()->GetConstrainedParam()) return;
 451
 452 // calculate and set prim. vertex
 453   fVertex->SetXv( infoMC->GetParticle().Vx() - dv[0] );
 454   fVertex->SetYv( infoMC->GetParticle().Vy() - dv[1] );
 455   fVertex->SetZv( infoMC->GetParticle().Vz() - dv[2] );
 456
 457   // constrained parameters resolution
 458   const AliExternalTrackParam * cparam = infoRC->GetESDtrack()->GetConstrainedParam();
 459   deltaPt= (mcpt-cparam->Pt())/mcpt;
 460   pullPt= (1/mcpt-cparam->OneOverPt())/TMath::Sqrt(cparam->GetSigma1Pt2());
 461   deltaPhi = TMath::ATan2(cparam->Py(),cparam->Px())-
 462                      TMath::ATan2(infoMC->GetParticle().Py(),infoMC->GetParticle().Px());
 463   pullPhi = deltaPhi/TMath::Sqrt(cparam->GetSigmaSnp2());
 464   deltaTan = TMath::ATan2(cparam->Pz(),cparam->Pt())-TMath::ATan2(infoMC->GetParticle().Pz(),infoMC->GetParticle().Pt());
 465   pullTan = deltaPhi/TMath::Sqrt(cparam->GetSigmaSnp2());
 466
 467
 468   delta1Pt2 = (1/mcpt-cparam->OneOverPt())/TMath::Power(1+1/mcpt,2);
 469
 470   deltaY1Pt = (infoMC->GetParticle().Vy()-cparam->GetY()) / (0.2+1/mcpt);
 471   deltaZ1Pt = (infoMC->GetParticle().Vz()-cparam->GetZ()) / (0.2+1/mcpt);
 472   deltaPhi1Pt = deltaPhi   / (0.1+1/mcpt);
 473   deltaTheta1Pt = deltaTan / (0.1+1/mcpt);
 474
 475   // calculate track parameters at vertex
 476   const AliExternalTrackParam *innerTPC =  0;
 477   if ((innerTPC = infoRC->GetESDtrack()->GetTPCInnerParam()) != 0)
 478   {
 479     if ((track = new AliExternalTrackParam(*infoRC->GetESDtrack()->GetTPCInnerParam())) != 0 )
 480     {
 481       Bool_t bStatus = AliTracker::PropagateTrackTo(track,kRadius,infoMC->GetMass(),kMaxStep,kTRUE);
 482       Bool_t bDCAStatus = track->PropagateToDCA(fVertex,field,kMaxD,dca,cov);
 483
 484       // Fill parametrisation histograms (only TPC track)
 485       if(bStatus && bDCAStatus)
 486           {
 487                   deltaPtTPC= (mcpt-innerTPC->Pt())/mcpt;
 488                   pullPtTPC= (1/mcpt-innerTPC->OneOverPt())/TMath::Sqrt(innerTPC->GetSigma1Pt2());
 489                   deltaPhiTPC = TMath::ATan2(innerTPC->Py(),innerTPC->Px())-
 490                                                                 TMath::ATan2(infoMC->GetParticle().Py(),infoMC->GetParticle().Px());
 491
 492                   deltaTanTPC = TMath::ATan2(innerTPC->Pz(),innerTPC->Pt())-
 493                                                                 TMath::ATan2(infoMC->GetParticle().Pz(),infoMC->GetParticle().Pt());
 494
 495                   delta1Pt2TPC = (1/mcpt-innerTPC->OneOverPt())/TMath::Power(1+1/mcpt,2);
 496                   deltaY1PtTPC= (infoMC->GetParticle().Vy()-innerTPC->GetY()) / (0.2+1/mcpt);
 497                   deltaZ1PtTPC = (infoMC->GetParticle().Vz()-innerTPC->GetZ()) / (0.2+1/mcpt);
 498                   deltaPhi1PtTPC = deltaPhiTPC   / (0.1+1/mcpt);
 499                   deltaTheta1PtTPC = deltaTanTPC / (0.1+1/mcpt);
 500
 501           fC1Pt2Resol1PtTPC->Fill(1/mcpt,delta1Pt2TPC);
 502           fCYResol1PtTPC->Fill(1/mcpt,deltaY1PtTPC);
 503           fCZResol1PtTPC->Fill(1/mcpt,deltaZ1PtTPC);
 504           fCPhiResol1PtTPC->Fill(1/mcpt,deltaPhi1PtTPC);
 505           fCThetaResol1PtTPC->Fill(1/mcpt,deltaTheta1PtTPC);
 506           }
 507           delete track;
 508     }
 509   }
 510
 511  // TPC and ITS (nb. of clusters >2) in the system
 512   if(infoRC->GetESDtrack()->GetITSclusters(clusterITS)>2)
 513   {
 514       fC1Pt2Resol1PtTPCITS->Fill(1/mcpt,delta1Pt2);
 515       fCYResol1PtTPCITS->Fill(1/mcpt,deltaY1Pt);
 516       fCZResol1PtTPCITS->Fill(1/mcpt,deltaZ1Pt);
 517       fCPhiResol1PtTPCITS->Fill(1/mcpt,deltaPhi1Pt);
 518       fCThetaResol1PtTPCITS->Fill(1/mcpt,deltaTheta1Pt);
 519   }
 520
 521   // Fill histograms
 522   fCPtResolTan->Fill(tantheta,deltaPt);
 523   fCPtPullTan->Fill(tantheta,pullPt);
 524   fCPhiResolTan->Fill(tantheta,deltaPhi);
 525   fCPhiPullTan->Fill(tantheta,pullPhi);
 526   fCTanResolTan->Fill(tantheta,deltaTan);
 527   fCTanPullTan->Fill(tantheta,pullTan);
 528 }
 529
 530 //_____________________________________________________________________________
 531 void AliComparisonRes::Exec(AliMCInfo* infoMC, AliESDRecInfo *infoRC){
 532
 533   // Process comparison information
 534   Process(infoMC,infoRC);
 535   ProcessConstrained(infoMC,infoRC);
 536 }
 537
 538 //_____________________________________________________________________________
 539 TH1F* AliComparisonRes::MakeResol(TH2F * his, Int_t integ, Bool_t type){
 540   // Create resolution histograms
 541
 542   TH1F *hisr, *hism;
 543   if (!gPad) new TCanvas;
 544   hisr = AliTreeDraw::CreateResHistoI(his,&hism,integ);
 545   if (type) return hism;
 546   else
 547     return hisr;
 548 }
 549
 550 //_____________________________________________________________________________
 551 void AliComparisonRes::Analyse(){
 552   // Analyse comparison information and store output histograms
 553   // in the "pictures_res.root" file
 554
 555   AliComparisonRes * comp=this;
 556   TH1F *hiss=0;
 557
 558   TFile *fp = new TFile("pictures_res.root","recreate");
 559   fp->cd();
 560
 561   TCanvas * c = new TCanvas("Phi resol Tan","Phi resol Tan");
 562   c->cd();
 563   //
 564   hiss = comp->MakeResol(comp->fCPtResolTan,1,0);
 565   hiss->SetXTitle("Tan(#theta)");
 566   hiss->SetYTitle("#sigmap_{t}/p_{t}");
 567   hiss->Draw();
 568   hiss->Write("CptResolTan");
 569   //
 570   hiss = comp->MakeResol(comp->fCPhiResolTan,1,0);
 571   hiss->SetXTitle("Tan(#theta)");
 572   hiss->SetYTitle("#sigma#phi (rad)");
 573   hiss->Draw();
 574   hiss->Write("PhiResolTan");
 575   //
 576   hiss = comp->MakeResol(comp->fCTanResolTan,1,0);
 577   hiss->SetXTitle("Tan(#theta)");
 578   hiss->SetYTitle("#sigma#theta (rad)");
 579   hiss->Draw();
 580   hiss->Write("ThetaResolTan");
 581   //
 582   hiss = comp->MakeResol(comp->fCPtPullTan,1,0);
 583   hiss->SetXTitle("Tan(#theta)");
 584   hiss->SetYTitle("1/mcp_{t}-1/p_{t}/#Sigma(1/p_{t})");
 585   hiss->Draw();
 586   hiss->Write("CptPullTan");
 587   //
 588   hiss = comp->MakeResol(comp->fC1Pt2Resol1PtTPC,1,0);
 589   hiss->SetXTitle("1/mcp_{t}");
 590   hiss->SetYTitle("1/mcp_{t}-1/p_{t}/(1+1/p_{t})^2");
 591   hiss->Draw();
 592   hiss->Write("C1Pt2Resol1PtTPC");
 593   fC1Pt2Resol1PtTPC->Write();
 594
 595   hiss = comp->MakeResol(comp->fC1Pt2Resol1PtTPCITS,1,0);
 596   hiss->SetXTitle("1/mcp_{t}");
 597   hiss->SetYTitle("1/mcp_{t}-1/p_{t}/(1+1/p_{t})^2");
 598   hiss->Draw();
 599   hiss->Write("C1Pt2Resol1PtTPCITS");
 600   fC1Pt2Resol1PtTPCITS->Write();
 601   //
 602   hiss = comp->MakeResol(comp->fCYResol1PtTPC,1,0);
 603   hiss->SetXTitle("1/mcp_{t}");
 604   hiss->SetYTitle("(mcy-y)/(0.2+1/mcp_{t})");
 605   hiss->Draw();
 606   hiss->Write("CYResol1PtTPC");
 607   fCYResol1PtTPC->Write();
 608
 609   hiss = comp->MakeResol(comp->fCYResol1PtTPCITS,1,0);
 610   hiss->SetXTitle("1/mcp_{t}");
 611   hiss->SetYTitle("(mcy-y)/(0.2+1/mcp_{t})");
 612   hiss->Draw();
 613   hiss->Write("CYResol1PtTPCITS");
 614   fCYResol1PtTPCITS->Write();
 615   //
 616   hiss = comp->MakeResol(comp->fCZResol1PtTPC,1,0);
 617   hiss->SetXTitle("1/mcp_{t}");
 618   hiss->SetYTitle("(mcz-z)/(0.2+1/mcp_{t})");
 619   hiss->Draw();
 620   hiss->Write("CZResol1PtTPC");
 621   fCZResol1PtTPC->Write();
 622
 623   hiss = comp->MakeResol(comp->fCZResol1PtTPCITS,1,0);
 624   hiss->SetXTitle("1/mcp_{t}");
 625   hiss->SetYTitle("(mcz-z)/(0.2+1/mcp_{t})");
 626   hiss->Draw();
 627   hiss->Write("CZResol1PtTPCITS");
 628   fCZResol1PtTPCITS->Write();
 629   //
 630   hiss = comp->MakeResol(comp->fCPhiResol1PtTPC,1,0);
 631   hiss->SetXTitle("1/mcp_{t}");
 632   hiss->SetYTitle("(mc#phi-#phi)/(0.1+1/mcp_{t})");
 633   hiss->Draw();
 634   hiss->Write("CPhiResol1PtTPC");
 635   fCPhiResol1PtTPC->Write();
 636
 637   hiss = comp->MakeResol(comp->fCPhiResol1PtTPCITS,1,0);
 638   hiss->SetXTitle("1/mcp_{t}");
 639   hiss->SetYTitle("(mc#phi-#phi)/(0.1+1/mcp_{t})");
 640   hiss->Draw();
 641   hiss->Write("CPhiResol1PtTPCITS");
 642   fCPhiResol1PtTPCITS->Write();
 643   //
 644   hiss = comp->MakeResol(comp->fCThetaResol1PtTPC,1,0);
 645   hiss->SetXTitle("1/mcp_{t}");
 646   hiss->SetYTitle("(mc#theta-#theta)/(0.1+1/mcp_{t})");
 647   hiss->Draw();
 648   hiss->Write("CThetaResol1PtTPC");
 649   fCThetaResol1PtTPC->Write();
 650
 651   hiss = comp->MakeResol(comp->fCThetaResol1PtTPCITS,1,0);
 652   hiss->SetXTitle("1/mcp_{t}");
 653   hiss->SetYTitle("(mc#theta-#theta)/(0.1+1/mcp_{t})");
 654   hiss->Draw();
 655   hiss->Write("CThetaResol1PtTPCITS");
 656   fCThetaResol1PtTPCITS->Write();
 657
 658   //
 659   hiss = comp->MakeResol(comp->f1Pt2Resol1PtTPC,1,0);
 660   hiss->SetXTitle("1/mcp_{t}");
 661   hiss->SetYTitle("1/mcp_{t}-1/p_{t}/(1+1/p_{t})^2");
 662   hiss->Draw();
 663   hiss->Write("OnePt2Resol1PtTPC");
 664   f1Pt2Resol1PtTPC->Write();
 665
 666   hiss = comp->MakeResol(comp->f1Pt2Resol1PtTPCITS,1,0);
 667   hiss->SetXTitle("1/mcp_{t}");
 668   hiss->SetYTitle("1/mcp_{t}-1/p_{t}/(1+1/p_{t})^2");
 669   hiss->Draw();
 670   hiss->Write("OnePt2Resol1PtTPCITS");
 671   f1Pt2Resol1PtTPCITS->Write();
 672   //
 673   hiss = comp->MakeResol(comp->fYResol1PtTPC,1,0);
 674   hiss->SetXTitle("1/mcp_{t}");
 675   hiss->SetYTitle("(mcy-y)/(0.2+1/mcp_{t})");
 676   hiss->Draw();
 677   hiss->Write("YResol1PtTPC");
 678   fYResol1PtTPC->Write();
 679
 680   hiss = comp->MakeResol(comp->fYResol1PtTPCITS,1,0);
 681   hiss->SetXTitle("1/mcp_{t}");
 682   hiss->SetYTitle("(mcy-y)/(0.2+1/mcp_{t})");
 683   hiss->Draw();
 684   hiss->Write("YResol1PtTPCITS");
 685   fYResol1PtTPCITS->Write();
 686   //
 687   hiss = comp->MakeResol(comp->fZResol1PtTPC,1,0);
 688   hiss->SetXTitle("1/mcp_{t}");
 689   hiss->SetYTitle("(mcz-z)/(0.2+1/mcp_{t})");
 690   hiss->Draw();
 691   hiss->Write("ZResol1PtTPC");
 692   fZResol1PtTPC->Write();
 693
 694   hiss = comp->MakeResol(comp->fZResol1PtTPCITS,1,0);
 695   hiss->SetXTitle("1/mcp_{t}");
 696   hiss->SetYTitle("(mcz-z)/(0.2+1/mcp_{t})");
 697   hiss->Draw();
 698   hiss->Write("ZResol1PtTPCITS");
 699   fZResol1PtTPCITS->Write();
 700   //
 701   hiss = comp->MakeResol(comp->fPhiResol1PtTPC,1,0);
 702   hiss->SetXTitle("1/mcp_{t}");
 703   hiss->SetYTitle("(mc#phi-#phi)/(0.1+1/mcp_{t})");
 704   hiss->Draw();
 705   hiss->Write("PhiResol1PtTPC");
 706   fPhiResol1PtTPC->Write();
 707
 708   hiss = comp->MakeResol(comp->fPhiResol1PtTPCITS,1,0);
 709   hiss->SetXTitle("1/mcp_{t}");
 710   hiss->SetYTitle("(mc#phi-#phi)/(0.1+1/mcp_{t})");
 711   hiss->Draw();
 712   hiss->Write("PhiResol1PtTPCITS");
 713   fPhiResol1PtTPCITS->Write();
 714   //
 715   hiss = comp->MakeResol(comp->fThetaResol1PtTPC,1,0);
 716   hiss->SetXTitle("1/mcp_{t}");
 717   hiss->SetYTitle("(mc#theta-#theta)/(0.1+1/mcp_{t})");
 718   hiss->Draw();
 719   hiss->Write("ThetaResol1PtTPC");
 720   fThetaResol1PtTPC->Write();
 721
 722   hiss = comp->MakeResol(comp->fThetaResol1PtTPCITS,1,0);
 723   hiss->SetXTitle("1/mcp_{t}");
 724   hiss->SetYTitle("(mc#theta-#theta)/(0.1+1/mcp_{t})");
 725   hiss->Draw();
 726   hiss->Write("ThetaResol1PtTPCITS");
 727   fThetaResol1PtTPCITS->Write();
 728
 729   fp->Close();
 730 }
 731
 732 //_____________________________________________________________________________
 733 Long64_t AliComparisonRes::Merge(TCollection* list)
 734 {
 735   // Merge list of objects (needed by PROOF)
 736
 737   if (!list)
 738   return 0;
 739
 740   if (list->IsEmpty())
 741   return 1;
 742
 743   TIterator* iter = list->MakeIterator();
 744   TObject* obj = 0;
 745
 746   // collection of generated histograms
 747   Int_t count=0;
 748   while((obj = iter->Next()) != 0)
 749   {
 750   AliComparisonRes* entry = dynamic_cast<AliComparisonRes*>(obj);
 751   if (entry == 0) continue;
 752
 753   fPtResolLPT->Add(entry->fPtResolLPT);
 754   fPtResolHPT->Add(entry->fPtResolHPT);
 755   fPtPullLPT->Add(entry->fPtPullLPT);
 756   fPtPullHPT->Add(entry->fPtPullHPT);
 757
 758   // Histograms for 1/pt parameterisation
 759   f1Pt2Resol1PtTPC->Add(entry->f1Pt2Resol1PtTPC);
 760   fYResol1PtTPC->Add(entry->fYResol1PtTPC);
 761   fZResol1PtTPC->Add(entry->fZResol1PtTPC);
 762   fPhiResol1PtTPC->Add(entry->fPhiResol1PtTPC);
 763   fThetaResol1PtTPC->Add(entry->fThetaResol1PtTPC);
 764
 765   f1Pt2Resol1PtTPCITS->Add(entry->f1Pt2Resol1PtTPCITS);
 766   fYResol1PtTPCITS->Add(entry->fYResol1PtTPCITS);
 767   fZResol1PtTPCITS->Add(entry->fZResol1PtTPCITS);
 768   fPhiResol1PtTPCITS->Add(entry->fPhiResol1PtTPCITS);
 769   fThetaResol1PtTPCITS->Add(entry->fThetaResol1PtTPCITS);
 770
 771   // Resolution histograms (constrained param)
 772   fCPhiResolTan->Add(entry->fCPhiResolTan);
 773   fCTanResolTan->Add(entry->fCTanResolTan);
 774   fCPtResolTan->Add(entry->fCPtResolTan);
 775   fCPhiPullTan->Add(entry->fCPhiPullTan);
 776   fCTanPullTan->Add(entry->fCTanPullTan);
 777   fCPtPullTan->Add(entry->fCPtPullTan);
 778
 779   //  Histograms for 1/pt parameterisation (constrained)
 780   fC1Pt2Resol1PtTPC->Add(entry->fC1Pt2Resol1PtTPC);
 781   fCYResol1PtTPC->Add(entry->fCYResol1PtTPC);
 782   fCZResol1PtTPC->Add(entry->fCZResol1PtTPC);
 783   fCPhiResol1PtTPC->Add(entry->fCPhiResol1PtTPC);
 784   fCThetaResol1PtTPC->Add(entry->fCThetaResol1PtTPC);
 785
 786   fC1Pt2Resol1PtTPCITS->Add(entry->fC1Pt2Resol1PtTPCITS);
 787   fCYResol1PtTPCITS->Add(entry->fCYResol1PtTPCITS);
 788   fCZResol1PtTPCITS->Add(entry->fCZResol1PtTPCITS);
 789   fCPhiResol1PtTPCITS->Add(entry->fCPhiResol1PtTPCITS);
 790   fCThetaResol1PtTPCITS->Add(entry->fCThetaResol1PtTPCITS);
 791
 792   count++;
 793   }
 794
 795 return count;
 796 }