PWG1/AliComparisonRes.cxx

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