PWG1/AliPerformanceRes.cxx

   1 //------------------------------------------------------------------------------
   2 // Implementation of AliPerformanceRes 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/graphs) are stored in the folder which is
   8 // a data member of AliPerformanceRes.
   9 //
  10 // Author: J.Otwinowski 04/02/2008
  11 //------------------------------------------------------------------------------
  12
  13 /*
  14
  15   // after running comparison task, read the file, and get component
  16   gROOT->LoadMacro("$ALICE_ROOT/PWG1/Macros/LoadMyLibs.C");
  17   LoadMyLibs();
  18
  19   TFile f("Output.root");
  20   AliPerformanceRes * compObj = (AliPerformanceRes*)coutput->FindObject("AliPerformanceRes");
  21
  22   // analyse comparison data
  23   compObj->Analyse();
  24
  25   // the output histograms/graphs will be stored in the folder "folderRes"
  26   compObj->GetAnalysisFolder()->ls("*");
  27
  28   // user can save whole comparison object (or only folder with anlysed histograms)
  29   // in the seperate output file (e.g.)
  30   TFile fout("Analysed_Res.root","recreate");
  31   compObj->Write(); // compObj->GetAnalysisFolder()->Write();
  32   fout.Close();
  33
  34 */
  35
  36 #include "TCanvas.h"
  37 #include "TH1.h"
  38 #include "TH2.h"
  39 #include "TAxis.h"
  40 #include "TF1.h"
  41
  42 #include "AliPerformanceRes.h"
  43 #include "AliESDEvent.h"
  44 #include "AliESDVertex.h"
  45 #include "AliESDtrack.h"
  46 #include "AliLog.h"
  47 #include "AliMCEvent.h"
  48 #include "AliMCParticle.h"
  49 #include "AliHeader.h"
  50 #include "AliGenEventHeader.h"
  51 #include "AliStack.h"
  52 #include "AliMCInfoCuts.h"
  53 #include "AliRecInfoCuts.h"
  54 #include "AliTracker.h"
  55 #include "AliTreeDraw.h"
  56
  57 using namespace std;
  58
  59 ClassImp(AliPerformanceRes)
  60
  61 //_____________________________________________________________________________
  62 AliPerformanceRes::AliPerformanceRes():
  63   AliPerformanceObject("AliPerformanceRes"),
  64   fResolHisto(0),
  65   fPullHisto(0),
  66
  67   // Cuts
  68   fCutsRC(0),
  69   fCutsMC(0),
  70
  71   // histogram folder
  72   fAnalysisFolder(0)
  73 {
  74   Init();
  75 }
  76
  77 //_____________________________________________________________________________
  78 AliPerformanceRes::AliPerformanceRes(Char_t* name="AliPerformanceRes", Char_t* title="AliPerformanceRes",Int_t analysisMode=0,Bool_t hptGenerator=kFALSE):
  79   AliPerformanceObject(name,title),
  80   fResolHisto(0),
  81   fPullHisto(0),
  82
  83   // Cuts
  84   fCutsRC(0),
  85   fCutsMC(0),
  86
  87   // histogram folder
  88   fAnalysisFolder(0)
  89 {
  90   // named constructor
  91   //
  92   SetAnalysisMode(analysisMode);
  93   SetHptGenerator(hptGenerator);
  94
  95   Init();
  96 }
  97
  98 //_____________________________________________________________________________
  99 AliPerformanceRes::~AliPerformanceRes()
 100 {
 101   // destructor
 102
 103   if(fResolHisto) delete fResolHisto; fResolHisto=0;
 104   if(fPullHisto)  delete fPullHisto;  fPullHisto=0;
 105
 106   if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
 107 }
 108
 109 //_____________________________________________________________________________
 110 void AliPerformanceRes::Init(){
 111
 112   //
 113   // histogram bining
 114   //
 115
 116   // set pt bins
 117   Int_t nPtBins = 50;
 118   Double_t ptMin = 1.e-2, ptMax = 10.;
 119
 120   Double_t *binsPt = 0;
 121   if (IsHptGenerator())  {
 122     nPtBins = 100; ptMax = 100.;
 123     binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
 124   } else {
 125     binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
 126   }
 127
 128   /*
 129   Int_t nPtBins = 31;
 130   Double_t binsPt[32] = {0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.7,0.8,0.9,1.0,1.2,1.4,1.6,1.8,2.0,2.25,2.5,2.75,3.,3.5,4.,5.,6.,8.,10.};
 131   Double_t ptMin = 0., ptMax = 10.;
 132
 133   if(IsHptGenerator() == kTRUE) {
 134     nPtBins = 100;
 135     ptMin = 0.; ptMax = 100.;
 136   }
 137   */
 138
 139   Double_t yMin = -0.02, yMax = 0.02;
 140   Double_t zMin = -12.0, zMax = 12.0;
 141   if(GetAnalysisMode() == 3) { // TrackRef coordinate system
 142     yMin = -100.; yMax = 100.;
 143     zMin = -100.; zMax = 100.;
 144   }
 145
 146   // res_y:res_z:res_phi,res_lambda:res_pt:y:z:eta:phi:pt
 147   Int_t binsResolHisto[10]={100,100,100,100,100,25,50,144,30,nPtBins};
 148   Double_t minResolHisto[10]={-1.,-1.,-0.03,-0.03,-0.2, yMin, zMin, 0., -1.5, ptMin};
 149   Double_t maxResolHisto[10]={ 1., 1., 0.03, 0.03, 0.2, yMax, zMax, 2.*TMath::Pi(), 1.5, ptMax};
 150
 151   fResolHisto = new THnSparseF("fResolHisto","res_y:res_z:res_phi:res_lambda:res_pt:y:z:phi:eta:pt",10,binsResolHisto,minResolHisto,maxResolHisto);
 152   fResolHisto->SetBinEdges(9,binsPt);
 153
 154   fResolHisto->GetAxis(0)->SetTitle("y-y_{mc} (cm)");
 155   fResolHisto->GetAxis(1)->SetTitle("z-z_{mc} (cm)");
 156   fResolHisto->GetAxis(2)->SetTitle("#phi-#phi_{mc} (rad)");
 157   fResolHisto->GetAxis(3)->SetTitle("#lambda-#lambda_{mc} (rad)");
 158   fResolHisto->GetAxis(4)->SetTitle("(p_{T}/p_{Tmc}-1)");
 159   fResolHisto->GetAxis(5)->SetTitle("y_{mc} (cm)");
 160   fResolHisto->GetAxis(6)->SetTitle("z_{mc} (cm)");
 161   fResolHisto->GetAxis(7)->SetTitle("#phi_{mc} (rad)");
 162   fResolHisto->GetAxis(8)->SetTitle("#eta_{mc}");
 163   fResolHisto->GetAxis(9)->SetTitle("p_{Tmc} (GeV/c)");
 164   fResolHisto->Sumw2();
 165
 166   ////pull_y:pull_z:pull_phi:pull_lambda:pull_1pt:y:z:eta:phi:pt
 167   //Int_t binsPullHisto[10]={100,100,100,100,100,50,50,30,144,nPtBins};
 168   //Double_t minPullHisto[10]={-5.,-5.,-5.,-5.,-5.,yMin, zMin,-1.5, 0., ptMin};
 169   //Double_t maxPullHisto[10]={ 5., 5., 5., 5., 5., yMax, zMax, 1.5, 2.*TMath::Pi(),ptMax};
 170   //fPullHisto = new THnSparseF("fPullHisto","pull_y:pull_z:pull_phi:pull_lambda:pull_1pt:y:z:eta:phi:pt",10,binsPullHisto,minPullHisto,maxPullHisto);
 171
 172   //pull_y:pull_z:pull_snp:pull_tgl:pull_1pt:y:z:snp:tgl:1pt
 173   Int_t binsPullHisto[10]={100,100,100,100,100,50,50,50,50,nPtBins};
 174   Double_t minPullHisto[10]={-5.,-5.,-5.,-5.,-5.,yMin, zMin,-1., -2.0, ptMin};
 175   Double_t maxPullHisto[10]={ 5., 5., 5., 5., 5., yMax, zMax, 1., 2.0, ptMax};
 176   fPullHisto = new THnSparseF("fPullHisto","pull_y:pull_z:pull_y:pull_z:pull_snp:pull_tgl:pull_1pt:y:z:snp:tgl:1pt",10,binsPullHisto,minPullHisto,maxPullHisto);
 177
 178   /*
 179   if(!IsHptGenerator()) fPullHisto->SetBinEdges(9,bins1Pt);
 180   fPullHisto->GetAxis(0)->SetTitle("(y-y_{mc})/#sigma");
 181   fPullHisto->GetAxis(1)->SetTitle("(z-z_{mc})/#sigma");
 182   fPullHisto->GetAxis(2)->SetTitle("(#phi-#phi_{mc})/#sigma");
 183   fPullHisto->GetAxis(3)->SetTitle("(#lambda-#lambda_{mc})/#sigma");
 184   fPullHisto->GetAxis(4)->SetTitle("(p_{Tmc}/p_{T}-1)/#sigma");
 185   fPullHisto->GetAxis(5)->SetTitle("y_{mc} (cm)");
 186   fPullHisto->GetAxis(6)->SetTitle("z_{mc} (cm)");
 187   fPullHisto->GetAxis(7)->SetTitle("#eta_{mc}");
 188   fPullHisto->GetAxis(8)->SetTitle("#phi_{mc} (rad)");
 189   fPullHisto->GetAxis(9)->SetTitle("p_{Tmc} (GeV/c)");
 190   fPullHisto->Sumw2();
 191   */
 192
 193   fPullHisto->GetAxis(0)->SetTitle("(y-y_{mc})/#sigma");
 194   fPullHisto->GetAxis(1)->SetTitle("(z-z_{mc})/#sigma");
 195   fPullHisto->GetAxis(2)->SetTitle("(sin#phi-sin#phi_{mc})/#sigma");
 196   fPullHisto->GetAxis(3)->SetTitle("(tan#lambda-tan#lambda_{mc})/#sigma");
 197   fPullHisto->GetAxis(4)->SetTitle("(p_{Tmc}/p_{T}-1)/#sigma");
 198   fPullHisto->GetAxis(5)->SetTitle("y_{mc} (cm)");
 199   fPullHisto->GetAxis(6)->SetTitle("z_{mc} (cm)");
 200   fPullHisto->GetAxis(7)->SetTitle("sin#phi_{mc}");
 201   fPullHisto->GetAxis(8)->SetTitle("tan#lambda_{mc}");
 202   fPullHisto->GetAxis(9)->SetTitle("1/p_{Tmc} (GeV/c)^{-1}");
 203   fPullHisto->Sumw2();
 204
 205   // Init cuts
 206   if(!fCutsMC)
 207     AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
 208   if(!fCutsRC)
 209     AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
 210
 211   // init folder
 212   fAnalysisFolder = CreateFolder("folderRes","Analysis Resolution Folder");
 213 }
 214
 215 //_____________________________________________________________________________
 216 void AliPerformanceRes::ProcessTPC(AliStack* const stack, AliESDtrack *const esdTrack)
 217 {
 218   if(!esdTrack) return;
 219
 220   // Fill TPC only resolution comparison information
 221   const AliExternalTrackParam *track = esdTrack->GetTPCInnerParam();
 222   if(!track) return;
 223
 224   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
 225   esdTrack->GetImpactParametersTPC(dca,cov);
 226
 227   //
 228   // Fill rec vs MC information
 229   //
 230   if(!stack) return;
 231   Int_t label = TMath::Abs(esdTrack->GetLabel());
 232   if(label == 0) return;
 233
 234   TParticle* particle = stack->Particle(label);
 235   if(!particle) return;
 236   if(!particle->GetPDG()) return;
 237   if(particle->GetPDG()->Charge()==0) return;
 238   //printf("charge %d \n",particle->GetPDG()->Charge());
 239
 240   Float_t deltaPtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC;
 241   Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC;
 242
 243   Float_t mceta =  particle->Eta();
 244   Float_t mcphi =  particle->Phi();
 245   if(mcphi<0) mcphi += 2.*TMath::Pi();
 246   Float_t mcpt = particle->Pt();
 247   Float_t mcsnp = TMath::Sin(TMath::ATan2(particle->Py(),particle->Px()));
 248   Float_t mctgl = TMath::Tan(TMath::ATan2(particle->Pz(),particle->Pt()));
 249
 250   // nb. TPC clusters cut
 251   if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;
 252
 253   // select primaries
 254   if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ())
 255   {
 256     if(mcpt == 0) return;
 257
 258     deltaYTPC= track->GetY()-particle->Vy();
 259     deltaZTPC = track->GetZ()-particle->Vz();
 260     deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(particle->Pz(),particle->Pt());
 261     deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(particle->Py(),particle->Px());
 262     //delta1PtTPC = (track->OneOverPt()-1./mcpt)*mcpt;
 263     deltaPtTPC = (track->Pt()-mcpt) / mcpt;
 264
 265     pullYTPC= (track->GetY()-particle->Vy()) / TMath::Sqrt(track->GetSigmaY2());
 266     pullZTPC = (track->GetZ()-particle->Vz()) / TMath::Sqrt(track->GetSigmaZ2());
 267
 268     //Double_t sigma_lambda = 1./(1.+track->GetTgl()*track->GetTgl()) * TMath::Sqrt(track->GetSigmaTgl2());
 269     //Double_t sigma_phi = 1./TMath::Sqrt(1-track->GetSnp()*track->GetSnp()) * TMath::Sqrt(track->GetSigmaSnp2());
 270     pullPhiTPC = (track->GetSnp() - mcsnp) / TMath::Sqrt(track->GetSigmaSnp2());
 271     pullLambdaTPC = (track->GetTgl() - mctgl) / TMath::Sqrt(track->GetSigmaTgl2());
 272
 273     //pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(track->GetSigmaTgl2());
 274     //pullPhiTPC = deltaPhiTPC / TMath::Sqrt(track->GetSigmaSnp2());
 275     if (mcpt) pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2());
 276     else pull1PtTPC = 0.;
 277
 278     Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,deltaPtTPC,particle->Vy(),particle->Vz(),mcphi,mceta,mcpt};
 279     fResolHisto->Fill(vResolHisto);
 280
 281     Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,particle->Vy(),particle->Vz(),mcsnp,mctgl,1./mcpt};
 282     fPullHisto->Fill(vPullHisto);
 283   }
 284 }
 285
 286 //_____________________________________________________________________________
 287 void AliPerformanceRes::ProcessTPCITS(AliStack* const stack, AliESDtrack *const esdTrack)
 288 {
 289   // Fill resolution comparison information (TPC+ITS)
 290   if(!esdTrack) return;
 291
 292   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
 293   esdTrack->GetImpactParameters(dca,cov);
 294
 295   //
 296   // Fill rec vs MC information
 297   //
 298   if(!stack) return;
 299
 300   Int_t label = TMath::Abs(esdTrack->GetLabel());
 301   if(label == 0) return;
 302
 303   TParticle* particle = stack->Particle(label);
 304   if(!particle) return;
 305   if(particle->GetPDG()->Charge()==0) return;
 306   //printf("charge %d \n",particle->GetPDG()->Charge());
 307
 308   Float_t mceta =  particle->Eta();
 309   Float_t mcphi =  particle->Phi();
 310   if(mcphi<0) mcphi += 2.*TMath::Pi();
 311   Float_t mcpt = particle->Pt();
 312   Float_t mcsnp = TMath::Sin(TMath::ATan2(particle->Py(),particle->Px()));
 313   Float_t mctgl = TMath::Tan(TMath::ATan2(particle->Pz(),particle->Pt()));
 314
 315   if ((esdTrack->GetStatus()&AliESDtrack::kTPCrefit)==0) return; // TPC refit
 316   if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; // min. nb. TPC clusters
 317   Int_t clusterITS[200];
 318   if(esdTrack->GetITSclusters(clusterITS)<fCutsRC->GetMinNClustersITS()) return;  // min. nb. ITS clusters
 319
 320   Float_t deltaPtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC;
 321   Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC;
 322
 323   // select primaries
 324   if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ())
 325   {
 326     if(mcpt == 0) return;
 327
 328     deltaYTPC= esdTrack->GetY()-particle->Vy();
 329     deltaZTPC = esdTrack->GetZ()-particle->Vz();
 330     deltaLambdaTPC = TMath::ATan2(esdTrack->Pz(),esdTrack->Pt())-TMath::ATan2(particle->Pz(),particle->Pt());
 331     deltaPhiTPC = TMath::ATan2(esdTrack->Py(),esdTrack->Px())-TMath::ATan2(particle->Py(),particle->Px());
 332     //delta1PtTPC = (esdTrack->OneOverPt()-1./mcpt)*mcpt;
 333     deltaPtTPC = (esdTrack->Pt()-mcpt) / mcpt;
 334
 335     pullYTPC= (esdTrack->GetY()-particle->Vy()) / TMath::Sqrt(esdTrack->GetSigmaY2());
 336     pullZTPC = (esdTrack->GetZ()-particle->Vz()) / TMath::Sqrt(esdTrack->GetSigmaZ2());
 337
 338     //Double_t sigma_lambda = 1./(1.+esdTrack->GetTgl()*esdTrack->GetTgl()) * TMath::Sqrt(esdTrack->GetSigmaTgl2());
 339     //Double_t sigma_phi = 1./TMath::Sqrt(1-esdTrack->GetSnp()*esdTrack->GetSnp()) * TMath::Sqrt(esdTrack->GetSigmaSnp2());
 340     pullPhiTPC = (esdTrack->GetSnp() - mcsnp) / TMath::Sqrt(esdTrack->GetSigmaSnp2());
 341     pullLambdaTPC = (esdTrack->GetTgl() - mctgl) / TMath::Sqrt(esdTrack->GetSigmaTgl2());
 342
 343     //pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(esdTrack->GetSigmaTgl2());
 344     //pullPhiTPC = deltaPhiTPC / TMath::Sqrt(esdTrack->GetSigmaSnp2());
 345     if (mcpt) pull1PtTPC = (esdTrack->OneOverPt()-1./mcpt) / TMath::Sqrt(esdTrack->GetSigma1Pt2());
 346     else pull1PtTPC = 0.;
 347
 348     Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,deltaPtTPC,particle->Vy(),particle->Vz(),mcphi,mceta,mcpt};
 349     fResolHisto->Fill(vResolHisto);
 350
 351     Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,particle->Vy(),particle->Vz(),mcsnp,mctgl,1./mcpt};
 352     fPullHisto->Fill(vPullHisto);
 353
 354
 355     /*
 356     deltaYTPC= esdTrack->GetY()-particle->Vy();
 357     deltaZTPC = esdTrack->GetZ()-particle->Vz();
 358     deltaLambdaTPC = TMath::ATan2(esdTrack->Pz(),esdTrack->Pt())-TMath::ATan2(particle->Pz(),particle->Pt());
 359     deltaPhiTPC = TMath::ATan2(esdTrack->Py(),esdTrack->Px())-TMath::ATan2(particle->Py(),particle->Px());
 360     delta1PtTPC = (esdTrack->OneOverPt()-1./mcpt)*mcpt;
 361
 362     pullYTPC= (esdTrack->GetY()-particle->Vy()) / TMath::Sqrt(esdTrack->GetSigmaY2());
 363     pullZTPC = (esdTrack->GetZ()-particle->Vz()) / TMath::Sqrt(esdTrack->GetSigmaZ2());
 364
 365     Double_t sigma_lambda = 1./(1.+esdTrack->GetTgl()*esdTrack->GetTgl()) * TMath::Sqrt(esdTrack->GetSigmaTgl2());
 366     Double_t sigma_phi = 1./TMath::Sqrt(1-esdTrack->GetSnp()*esdTrack->GetSnp()) * TMath::Sqrt(esdTrack->GetSigmaSnp2());
 367     pullLambdaTPC = deltaLambdaTPC / sigma_lambda;
 368     pullPhiTPC = deltaPhiTPC / sigma_phi;
 369
 370     //pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(esdTrack->GetSigmaTgl2());
 371     //pullPhiTPC = deltaPhiTPC / TMath::Sqrt(esdTrack->GetSigmaSnp2());
 372     if (mcpt) pull1PtTPC = (esdTrack->OneOverPt()-1./mcpt) / TMath::Sqrt(esdTrack->GetSigma1Pt2());
 373     else pull1PtTPC = 0.;
 374
 375     Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,delta1PtTPC,particle->Vy(),particle->Vz(),mceta,mcphi,mcpt};
 376     fResolHisto->Fill(vResolHisto);
 377
 378     Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,particle->Vy(),particle->Vz(),mceta,mcphi,mcpt};
 379     fPullHisto->Fill(vPullHisto);
 380     */
 381   }
 382 }
 383
 384 //_____________________________________________________________________________
 385 void AliPerformanceRes::ProcessConstrained(AliStack* const stack, AliESDtrack *const esdTrack)
 386 {
 387   // Fill resolution comparison information (constarained parameters)
 388   //
 389   if(!esdTrack) return;
 390
 391   const AliExternalTrackParam * track = esdTrack->GetConstrainedParam();
 392   if(!track) return;
 393
 394   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
 395   esdTrack->GetImpactParameters(dca,cov);
 396
 397   //
 398   // Fill rec vs MC information
 399   //
 400   if(!stack) return;
 401
 402   Int_t label = TMath::Abs(esdTrack->GetLabel());
 403   if(label == 0) return;
 404
 405   TParticle* particle = stack->Particle(label);
 406   if(!particle) return;
 407   if(particle->GetPDG()->Charge()==0) return;
 408   //printf("charge %d \n",particle->GetPDG()->Charge());
 409
 410   Float_t mceta =  particle->Eta();
 411   Float_t mcphi =  particle->Phi();
 412   if(mcphi<0) mcphi += 2.*TMath::Pi();
 413   Float_t mcpt = particle->Pt();
 414   Float_t mcsnp = TMath::Sin(TMath::ATan2(particle->Py(),particle->Px()));
 415   Float_t mctgl = TMath::Tan(TMath::ATan2(particle->Pz(),particle->Pt()));
 416
 417   if ((esdTrack->GetStatus()&AliESDtrack::kTPCrefit)==0) return; // TPC refit
 418   if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; // min. nb. TPC clusters
 419
 420   Float_t deltaPtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC;
 421   Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC;
 422
 423   // select primaries
 424   if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ())
 425   {
 426
 427     if(mcpt == 0) return;
 428
 429     deltaYTPC= track->GetY()-particle->Vy();
 430     deltaZTPC = track->GetZ()-particle->Vz();
 431     deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(particle->Pz(),particle->Pt());
 432     deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(particle->Py(),particle->Px());
 433     //delta1PtTPC = (track->OneOverPt()-1./mcpt)*mcpt;
 434     deltaPtTPC = (track->Pt()-mcpt) / mcpt;
 435
 436     pullYTPC= (track->GetY()-particle->Vy()) / TMath::Sqrt(track->GetSigmaY2());
 437     pullZTPC = (track->GetZ()-particle->Vz()) / TMath::Sqrt(track->GetSigmaZ2());
 438
 439     //Double_t sigma_lambda = 1./(1.+track->GetTgl()*track->GetTgl()) * TMath::Sqrt(track->GetSigmaTgl2());
 440     //Double_t sigma_phi = 1./TMath::Sqrt(1-track->GetSnp()*track->GetSnp()) * TMath::Sqrt(track->GetSigmaSnp2());
 441     pullPhiTPC = (track->GetSnp() - mcsnp) / TMath::Sqrt(track->GetSigmaSnp2());
 442     pullLambdaTPC = (track->GetTgl() - mctgl) / TMath::Sqrt(track->GetSigmaTgl2());
 443
 444     //pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(track->GetSigmaTgl2());
 445     //pullPhiTPC = deltaPhiTPC / TMath::Sqrt(track->GetSigmaSnp2());
 446     if (mcpt) pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2());
 447     else pull1PtTPC = 0.;
 448
 449     Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,deltaPtTPC,particle->Vy(),particle->Vz(),mcphi,mceta,mcpt};
 450     fResolHisto->Fill(vResolHisto);
 451
 452     Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,particle->Vy(),particle->Vz(),mcsnp,mctgl,1./mcpt};
 453     fPullHisto->Fill(vPullHisto);
 454
 455     /*
 456
 457     deltaYTPC= track->GetY()-particle->Vy();
 458     deltaZTPC = track->GetZ()-particle->Vz();
 459     deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(particle->Pz(),particle->Pt());
 460     deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(particle->Py(),particle->Px());
 461     delta1PtTPC = (track->OneOverPt()-1./mcpt)*mcpt;
 462
 463     pullYTPC= (track->GetY()-particle->Vy()) / TMath::Sqrt(track->GetSigmaY2());
 464     pullZTPC = (track->GetZ()-particle->Vz()) / TMath::Sqrt(track->GetSigmaZ2());
 465
 466     Double_t sigma_lambda = 1./(1.+track->GetTgl()*track->GetTgl()) * TMath::Sqrt(track->GetSigmaTgl2());
 467     Double_t sigma_phi = 1./TMath::Sqrt(1-track->GetSnp()*track->GetSnp()) * TMath::Sqrt(track->GetSigmaSnp2());
 468     pullLambdaTPC = deltaLambdaTPC / sigma_lambda;
 469     pullPhiTPC = deltaPhiTPC / sigma_phi;
 470
 471     //pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(track->GetSigmaTgl2());
 472     //pullPhiTPC = deltaPhiTPC / TMath::Sqrt(track->GetSigmaSnp2());
 473
 474     if (mcpt) pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2());
 475     else pull1PtTPC = 0.;
 476
 477     Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,delta1PtTPC,particle->Vy(),particle->Vz(),mceta,mcphi,mcpt};
 478     fResolHisto->Fill(vResolHisto);
 479
 480     Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,particle->Vy(),particle->Vz(),mceta,mcphi,mcpt};
 481     fPullHisto->Fill(vPullHisto);
 482
 483     */
 484   }
 485 }
 486
 487 //_____________________________________________________________________________
 488 void AliPerformanceRes::ProcessInnerTPC(AliMCEvent *const mcEvent, AliESDtrack *const esdTrack)
 489 {
 490   //
 491   // Fill resolution comparison information (inner params at TPC reference point)
 492   //
 493   if(!esdTrack) return;
 494
 495   const AliExternalTrackParam * innerParam = esdTrack->GetInnerParam();
 496   if(!innerParam) return;
 497
 498   // create new AliExternalTrackParam
 499   AliExternalTrackParam *track = new AliExternalTrackParam(*innerParam);
 500   if(!track) return;
 501
 502   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
 503   esdTrack->GetImpactParameters(dca,cov);
 504
 505   //
 506   // Fill rec vs MC information
 507   //
 508   if(!mcEvent) return;
 509
 510   Int_t label = TMath::Abs(esdTrack->GetLabel());
 511   if(label == 0) return;
 512   AliMCParticle *mcParticle = mcEvent->GetTrack(label);
 513   if(!mcParticle) return;
 514
 515   // get the first TPC track reference
 516   AliTrackReference *ref0 = GetFirstTPCTrackRef(mcParticle);
 517   if(!ref0) return;
 518
 519   // calculate alpha angle
 520   Double_t xyz[3] = {ref0->X(),ref0->Y(),ref0->Z()};
 521   Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
 522   //if (alpha<0) alpha += TMath::TwoPi();
 523
 524   // rotate inner track to local coordinate system
 525   // and propagate to the radius of the first track referenco of TPC
 526   Double_t trRadius = TMath::Sqrt(xyz[1] * xyz[1] + xyz[0] * xyz[0]);
 527   Bool_t isOK = track->Propagate(alpha,trRadius,AliTracker::GetBz());
 528   if(!isOK) return;
 529
 530   Float_t mceta =  -TMath::Log(TMath::Tan(0.5 * ref0->Theta()));
 531   Float_t mcphi =  ref0->Phi();
 532   if(mcphi<0) mcphi += 2.*TMath::Pi();
 533   Float_t mcpt = ref0->Pt();
 534   Float_t mcsnp = TMath::Sin(TMath::ATan2(ref0->Py(),ref0->Px()));
 535   Float_t mctgl = TMath::Tan(TMath::ATan2(ref0->Pz(),ref0->Pt()));
 536
 537   if ((esdTrack->GetStatus()&AliESDtrack::kTPCrefit)==0) return; // TPC refit
 538   if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; // min. nb. TPC clusters
 539
 540   Float_t deltaPtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC;
 541   Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC;
 542
 543   // select primaries
 544   if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ())
 545   {
 546     if(mcpt == 0) return;
 547
 548     deltaYTPC= track->GetY(); // already rotated
 549     deltaZTPC = track->GetZ()-ref0->Z();
 550     deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(ref0->Pz(),ref0->Pt());
 551     deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(ref0->Py(),ref0->Px());
 552     //delta1PtTPC = (track->OneOverPt()-1./mcpt)*mcpt;
 553     deltaPtTPC = (track->Pt()-mcpt) / mcpt;
 554
 555     pullYTPC= track->GetY() / TMath::Sqrt(track->GetSigmaY2());
 556     pullZTPC = (track->GetZ()-ref0->Z()) / TMath::Sqrt(track->GetSigmaZ2());
 557
 558     //Double_t sigma_lambda = 1./(1.+track->GetTgl()*track->GetTgl()) * TMath::Sqrt(track->GetSigmaTgl2());
 559     //Double_t sigma_phi = 1./TMath::Sqrt(1-track->GetSnp()*track->GetSnp()) * TMath::Sqrt(track->GetSigmaSnp2());
 560     pullPhiTPC = (track->GetSnp() - mcsnp) / TMath::Sqrt(track->GetSigmaSnp2());
 561     pullLambdaTPC = (track->GetTgl() - mctgl) / TMath::Sqrt(track->GetSigmaTgl2());
 562
 563     //pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(track->GetSigmaTgl2());
 564     //pullPhiTPC = deltaPhiTPC / TMath::Sqrt(track->GetSigmaSnp2());
 565     if (mcpt) pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2());
 566     else pull1PtTPC = 0.;
 567
 568     Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,deltaPtTPC,ref0->Y(),ref0->Z(),mcphi,mceta,mcpt};
 569     fResolHisto->Fill(vResolHisto);
 570
 571     Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,ref0->Y(),ref0->Z(),mcsnp,mctgl,1./mcpt};
 572     fPullHisto->Fill(vPullHisto);
 573
 574
 575     /*
 576
 577     //deltaYTPC= track->GetY()-ref0->Y();
 578     deltaYTPC= track->GetY(); // it corresponds to deltaY after alpha rotation
 579     deltaZTPC = track->GetZ()-ref0->Z();
 580     deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(ref0->Pz(),ref0->Pt());
 581     deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(ref0->Py(),ref0->Px());
 582     if(mcpt) delta1PtTPC = (track->OneOverPt()-1./mcpt)*mcpt;
 583     else delta1PtTPC = 0.;
 584
 585     //pullYTPC= (track->GetY()-ref0->Y()) / TMath::Sqrt(track->GetSigmaY2());
 586     pullYTPC= track->GetY() / TMath::Sqrt(track->GetSigmaY2()); //it corresponds to deltaY after alpha rotation
 587     pullZTPC = (track->GetZ()-ref0->Z()) / TMath::Sqrt(track->GetSigmaZ2());
 588
 589     //pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(track->GetSigmaTgl2());
 590     //pullPhiTPC = deltaPhiTPC / TMath::Sqrt(track->GetSigmaSnp2());
 591
 592     Double_t sigma_lambda = 1./(1.+track->GetTgl()*track->GetTgl()) * TMath::Sqrt(track->GetSigmaTgl2());
 593     Double_t sigma_phi = 1./TMath::Sqrt(1-track->GetSnp()*track->GetSnp()) * TMath::Sqrt(track->GetSigmaSnp2());
 594     pullLambdaTPC = deltaLambdaTPC / sigma_lambda;
 595     pullPhiTPC = deltaPhiTPC / sigma_phi;
 596
 597     if(mcpt) pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2());
 598     else pull1PtTPC = 0.;
 599
 600     Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,delta1PtTPC,ref0->Y(),ref0->Z(),mceta,mcphi,mcpt};
 601     fResolHisto->Fill(vResolHisto);
 602
 603     Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,ref0->Y(),ref0->Z(),mceta,mcphi,mcpt};
 604     fPullHisto->Fill(vPullHisto);
 605
 606     */
 607   }
 608
 609   if(track) delete track;
 610 }
 611
 612 //_____________________________________________________________________________
 613 AliTrackReference * AliPerformanceRes::GetFirstTPCTrackRef(AliMCParticle *mcParticle)
 614 {
 615   // return first TPC track reference
 616   if(!mcParticle) return 0;
 617
 618   Int_t nTrackRef = mcParticle->GetNumberOfTrackReferences();
 619   AliTrackReference *ref0 = 0;
 620   for (Int_t iref = 0; iref < nTrackRef; iref++) {
 621     ref0 = mcParticle->GetTrackReference(iref);
 622     if(ref0 && (ref0->DetectorId()==AliTrackReference::kTPC))
 623       break;
 624   }
 625
 626 return ref0;
 627 }
 628
 629 //_____________________________________________________________________________
 630 void AliPerformanceRes::Exec(AliMCEvent* const mcEvent, AliESDEvent* const esdEvent, const Bool_t bUseMC)
 631 {
 632   // Process comparison information
 633   //
 634   if(!esdEvent)
 635   {
 636       AliDebug(AliLog::kError, "esdEvent not available");
 637       return;
 638   }
 639   AliHeader* header = 0;
 640   AliGenEventHeader* genHeader = 0;
 641   AliStack* stack = 0;
 642   TArrayF vtxMC(3);
 643
 644   if(bUseMC)
 645   {
 646     if(!mcEvent) {
 647       AliDebug(AliLog::kError, "mcEvent not available");
 648       return;
 649     }
 650
 651     // get MC event header
 652     header = mcEvent->Header();
 653     if (!header) {
 654       AliDebug(AliLog::kError, "Header not available");
 655       return;
 656     }
 657     // MC particle stack
 658     stack = mcEvent->Stack();
 659     if (!stack) {
 660       AliDebug(AliLog::kError, "Stack not available");
 661       return;
 662     }
 663
 664     // get MC vertex
 665     genHeader = header->GenEventHeader();
 666     if (!genHeader) {
 667       AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
 668       return;
 669     }
 670     genHeader->PrimaryVertex(vtxMC);
 671
 672   } // end bUseMC
 673
 674   //  Process events
 675   for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
 676   {
 677     AliESDtrack*track = esdEvent->GetTrack(iTrack);
 678     if(!track) continue;
 679
 680     if(GetAnalysisMode() == 0) ProcessTPC(stack,track);
 681     else if(GetAnalysisMode() == 1) ProcessTPCITS(stack,track);
 682     else if(GetAnalysisMode() == 2) ProcessConstrained(stack,track);
 683     else if(GetAnalysisMode() == 3) ProcessInnerTPC(mcEvent,track);
 684     else {
 685       printf("ERROR: AnalysisMode %d \n",fAnalysisMode);
 686       return;
 687     }
 688   }
 689 }
 690
 691 //_____________________________________________________________________________
 692 TH1F* AliPerformanceRes::MakeResol(TH2F * his, Int_t integ, Bool_t type, Int_t cut){
 693   // Create resolution histograms
 694
 695   TH1F *hisr, *hism;
 696   if (!gPad) new TCanvas;
 697   hisr = AliTreeDraw::CreateResHistoII(his,&hism,integ,kTRUE,cut);
 698   if (type) return hism;
 699   else
 700     return hisr;
 701 }
 702
 703 //_____________________________________________________________________________
 704 void AliPerformanceRes::Analyse() {
 705   // Analyse comparison information and store output histograms
 706   // in the folder "folderRes"
 707   //
 708   TH1::AddDirectory(kFALSE);
 709   TH1F *h=0;
 710   TH2F *h2D=0;
 711   TObjArray *aFolderObj = new TObjArray;
 712
 713   // write results in the folder
 714   TCanvas * c = new TCanvas("Phi resol Tan","Phi resol Tan");
 715   c->cd();
 716
 717   char name[256];
 718   char title[256];
 719
 720   for(Int_t i=0; i<5; i++)
 721   {
 722     for(Int_t j=5; j<10; j++)
 723     {
 724       if(j!=8) fResolHisto->GetAxis(8)->SetRangeUser(-0.9,0.9); // eta window
 725       fResolHisto->GetAxis(9)->SetRangeUser(0.16,10.); // pt threshold
 726       if(GetAnalysisMode() == 3) fResolHisto->GetAxis(5)->SetRangeUser(-80.,80.); // y range
 727
 728       h2D = (TH2F*)fResolHisto->Projection(i,j);
 729
 730       h = AliPerformanceRes::MakeResol(h2D,1,0,100);
 731       sprintf(name,"h_res_%d_vs_%d",i,j);
 732       h->SetName(name);
 733
 734       h->GetXaxis()->SetTitle(fResolHisto->GetAxis(j)->GetTitle());
 735       sprintf(title,"%s %s",fResolHisto->GetAxis(i)->GetTitle(),"(resolution)");
 736       h->GetYaxis()->SetTitle(title);
 737       sprintf(title,"%s vs %s",title,fResolHisto->GetAxis(j)->GetTitle());
 738       h->SetTitle(title);
 739
 740       if(j==9) h->SetBit(TH1::kLogX);
 741       aFolderObj->Add(h);
 742
 743       h = AliPerformanceRes::MakeResol(h2D,1,1,100);
 744       //h = (TH1F*)arr->At(1);
 745       sprintf(name,"h_mean_res_%d_vs_%d",i,j);
 746       h->SetName(name);
 747
 748       h->GetXaxis()->SetTitle(fResolHisto->GetAxis(j)->GetTitle());
 749       sprintf(title,"%s %s",fResolHisto->GetAxis(i)->GetTitle(),"(mean)");
 750       h->GetYaxis()->SetTitle(title);
 751
 752       sprintf(title,"%s vs %s",title,fResolHisto->GetAxis(j)->GetTitle());
 753       h->SetTitle(title);
 754
 755       if(j==9) h->SetBit(TH1::kLogX);
 756       aFolderObj->Add(h);
 757
 758       fResolHisto->GetAxis(8)->SetRangeUser(-1.5,1.5);
 759       fResolHisto->GetAxis(9)->SetRangeUser(0.0,10.);
 760
 761       //
 762       if(j!=8) fPullHisto->GetAxis(8)->SetRangeUser(-1.0,1.0); // theta window
 763       else  fPullHisto->GetAxis(8)->SetRangeUser(-1.5,1.5); // theta window
 764       fPullHisto->GetAxis(9)->SetRangeUser(0.0,10.);  // pt threshold
 765       if(GetAnalysisMode() == 3) fPullHisto->GetAxis(5)->SetRangeUser(-80.,80.); // y range
 766
 767       h2D = (TH2F*)fPullHisto->Projection(i,j);
 768
 769       h = AliPerformanceRes::MakeResol(h2D,1,0,100);
 770       sprintf(name,"h_pull_%d_vs_%d",i,j);
 771       h->SetName(name);
 772
 773       h->GetXaxis()->SetTitle(fPullHisto->GetAxis(j)->GetTitle());
 774       sprintf(title,"%s %s",fPullHisto->GetAxis(i)->GetTitle(),"(resolution)");
 775       h->GetYaxis()->SetTitle(title);
 776       sprintf(title,"%s vs %s",title,fPullHisto->GetAxis(j)->GetTitle());
 777       h->SetTitle(title);
 778
 779       //if(j==9) h->SetBit(TH1::kLogX);
 780       aFolderObj->Add(h);
 781
 782       h = AliPerformanceRes::MakeResol(h2D,1,1,100);
 783       sprintf(name,"h_mean_pull_%d_vs_%d",i,j);
 784       h->SetName(name);
 785
 786       h->GetXaxis()->SetTitle(fPullHisto->GetAxis(j)->GetTitle());
 787       sprintf(title,"%s %s",fPullHisto->GetAxis(i)->GetTitle(),"(mean)");
 788       h->GetYaxis()->SetTitle(title);
 789       sprintf(title,"%s vs %s",title,fPullHisto->GetAxis(j)->GetTitle());
 790       h->SetTitle(title);
 791
 792       //if(j==9) h->SetBit(TH1::kLogX);
 793       aFolderObj->Add(h);
 794     }
 795   }
 796
 797   // export objects to analysis folder
 798   fAnalysisFolder = ExportToFolder(aFolderObj);
 799
 800   // delete only TObjArray
 801   if(aFolderObj) delete aFolderObj;
 802 }
 803
 804 //_____________________________________________________________________________
 805 TFolder* AliPerformanceRes::ExportToFolder(TObjArray * array)
 806 {
 807   // recreate folder avery time and export objects to new one
 808   //
 809   AliPerformanceRes * comp=this;
 810   TFolder *folder = comp->GetAnalysisFolder();
 811
 812   TString name, title;
 813   TFolder *newFolder = 0;
 814   Int_t i = 0;
 815   Int_t size = array->GetSize();
 816
 817   if(folder) {
 818      // get name and title from old folder
 819      name = folder->GetName();
 820      title = folder->GetTitle();
 821
 822          // delete old one
 823      delete folder;
 824
 825          // create new one
 826      newFolder = CreateFolder(name.Data(),title.Data());
 827      newFolder->SetOwner();
 828
 829          // add objects to folder
 830      while(i < size) {
 831            newFolder->Add(array->At(i));
 832            i++;
 833          }
 834   }
 835
 836 return newFolder;
 837 }
 838
 839 //_____________________________________________________________________________
 840 Long64_t AliPerformanceRes::Merge(TCollection* const list)
 841 {
 842   // Merge list of objects (needed by PROOF)
 843
 844   if (!list)
 845   return 0;
 846
 847   if (list->IsEmpty())
 848   return 1;
 849
 850   TIterator* iter = list->MakeIterator();
 851   TObject* obj = 0;
 852
 853   // collection of generated histograms
 854   Int_t count=0;
 855   while((obj = iter->Next()) != 0)
 856   {
 857   AliPerformanceRes* entry = dynamic_cast<AliPerformanceRes*>(obj);
 858   if (entry == 0) continue;
 859
 860   fResolHisto->Add(entry->fResolHisto);
 861   fPullHisto->Add(entry->fPullHisto);
 862
 863   count++;
 864   }
 865
 866 return count;
 867 }
 868
 869 //_____________________________________________________________________________
 870 TFolder* AliPerformanceRes::CreateFolder(TString name,TString title) {
 871 // create folder for analysed histograms
 872 //
 873 TFolder *folder = 0;
 874   folder = new TFolder(name.Data(),title.Data());
 875
 876   return folder;
 877 }