PWGPP/TPC/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/PWGPP/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 "AliESDfriendTrack.h"
  47 #include "AliESDfriend.h"
  48 #include "AliLog.h"
  49 #include "AliMCEvent.h"
  50 #include "AliMCParticle.h"
  51 #include "AliHeader.h"
  52 #include "AliGenEventHeader.h"
  53 #include "AliStack.h"
  54 #include "AliMCInfoCuts.h"
  55 #include "AliRecInfoCuts.h"
  56 #include "AliTracker.h"
  57 #include "AliTreeDraw.h"
  58
  59 using namespace std;
  60
  61 ClassImp(AliPerformanceRes)
  62 Double_t          AliPerformanceRes::fgkMergeEntriesCut=5000000.; //5*10**6 tracks (small default to keep default memory foorprint low)
  63
  64 //_____________________________________________________________________________
  65 AliPerformanceRes::AliPerformanceRes(const Char_t* name, const Char_t* title, Int_t analysisMode, Bool_t hptGenerator):
  66   AliPerformanceObject(name,title),
  67   fResolHisto(0),
  68   fPullHisto(0),
  69
  70   // Cuts
  71   fCutsRC(0),
  72   fCutsMC(0),
  73
  74   // histogram folder
  75   fAnalysisFolder(0)
  76 {
  77   // named constructor
  78   //
  79   SetAnalysisMode(analysisMode);
  80   SetHptGenerator(hptGenerator);
  81
  82   Init();
  83 }
  84
  85 //_____________________________________________________________________________
  86 AliPerformanceRes::~AliPerformanceRes()
  87 {
  88   // destructor
  89
  90   if(fResolHisto) delete fResolHisto; fResolHisto=0;
  91   if(fPullHisto)  delete fPullHisto;  fPullHisto=0;
  92
  93   if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
  94 }
  95
  96 //_____________________________________________________________________________
  97 void AliPerformanceRes::Init(){
  98
  99   //
 100   // histogram bining
 101   //
 102
 103   // set pt bins
 104   Int_t nPtBins = 50;
 105   Double_t ptMin = 1.e-1, ptMax = 20.;
 106
 107   Double_t *binsPt = 0;
 108
 109   if (IsHptGenerator())  {
 110         ptMax = 100.;
 111   }
 112    binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
 113
 114   Double_t yMin = -0.02, yMax = 0.02;
 115   Double_t zMin = -12.0, zMax = 12.0;
 116   if(GetAnalysisMode() == 3) { // TrackRef coordinate system
 117     yMin = -100.; yMax = 100.;
 118     zMin = -100.; zMax = 100.;
 119   }
 120
 121   // res_y:res_z:res_phi,res_lambda:res_pt:y:z:eta:phi:pt
 122   Int_t binsResolHisto[10]={100,100,100,100,100,25,50,144,30,nPtBins};
 123   Double_t minResolHisto[10]={-1.,-1.,-0.03,-0.03,-0.2, yMin, zMin, 0., -1.5, ptMin};
 124   Double_t maxResolHisto[10]={ 1., 1., 0.03, 0.03, 0.2, yMax, zMax, 2.*TMath::Pi(), 1.5, ptMax};
 125
 126   fResolHisto = new THnSparseF("fResolHisto","res_y:res_z:res_phi:res_lambda:res_pt:y:z:phi:eta:pt",10,binsResolHisto,minResolHisto,maxResolHisto);
 127
 128   //fResolHisto->SetBinEdges(9,binsPt);
 129   fResolHisto->GetAxis(9)->Set(nPtBins,binsPt);
 130
 131   fResolHisto->GetAxis(0)->SetTitle("y-y_{mc} (cm)");
 132   fResolHisto->GetAxis(1)->SetTitle("z-z_{mc} (cm)");
 133   fResolHisto->GetAxis(2)->SetTitle("#phi-#phi_{mc} (rad)");
 134   fResolHisto->GetAxis(3)->SetTitle("#lambda-#lambda_{mc} (rad)");
 135   fResolHisto->GetAxis(4)->SetTitle("(p_{T}/p_{Tmc}-1)");
 136   fResolHisto->GetAxis(5)->SetTitle("y_{mc} (cm)");
 137   fResolHisto->GetAxis(6)->SetTitle("z_{mc} (cm)");
 138   fResolHisto->GetAxis(7)->SetTitle("#phi_{mc} (rad)");
 139   fResolHisto->GetAxis(8)->SetTitle("#eta_{mc}");
 140   fResolHisto->GetAxis(9)->SetTitle("p_{Tmc} (GeV/c)");
 141   fResolHisto->Sumw2();
 142
 143   ////pull_y:pull_z:pull_phi:pull_lambda:pull_1pt:y:z:eta:phi:pt
 144   //Int_t binsPullHisto[10]={100,100,100,100,100,50,50,30,144,nPtBins};
 145   //Double_t minPullHisto[10]={-5.,-5.,-5.,-5.,-5.,yMin, zMin,-1.5, 0., ptMin};
 146   //Double_t maxPullHisto[10]={ 5., 5., 5., 5., 5., yMax, zMax, 1.5, 2.*TMath::Pi(),ptMax};
 147   //fPullHisto = new THnSparseF("fPullHisto","pull_y:pull_z:pull_phi:pull_lambda:pull_1pt:y:z:eta:phi:pt",10,binsPullHisto,minPullHisto,maxPullHisto);
 148
 149   //pull_y:pull_z:pull_snp:pull_tgl:pull_1pt:y:z:snp:tgl:1pt
 150   //Int_t binsPullHisto[10]={100,100,100,100,100,50,50,50,50,nPtBins};
 151   //Double_t minPullHisto[10]={-5.,-5.,-5.,-5.,-5.,yMin, zMin,-1., -2.0, ptMin};
 152   //Double_t maxPullHisto[10]={ 5., 5., 5., 5., 5., yMax, zMax, 1., 2.0, ptMax};
 153   Int_t binsPullHisto[10]={100,100,100,100,100,50,50,50,50,20};
 154   Double_t minPullHisto[10]={-5.,-5.,-5.,-5.,-5.,yMin, zMin,-1., -2.0, 0.};
 155   Double_t maxPullHisto[10]={ 5., 5., 5., 5., 5., yMax, zMax, 1., 2.0, 10.};
 156   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);
 157
 158   /*
 159   if(!IsHptGenerator()) fPullHisto->SetBinEdges(9,bins1Pt);
 160   fPullHisto->GetAxis(0)->SetTitle("(y-y_{mc})/#sigma");
 161   fPullHisto->GetAxis(1)->SetTitle("(z-z_{mc})/#sigma");
 162   fPullHisto->GetAxis(2)->SetTitle("(#phi-#phi_{mc})/#sigma");
 163   fPullHisto->GetAxis(3)->SetTitle("(#lambda-#lambda_{mc})/#sigma");
 164   fPullHisto->GetAxis(4)->SetTitle("(p_{Tmc}/p_{T}-1)/#sigma");
 165   fPullHisto->GetAxis(5)->SetTitle("y_{mc} (cm)");
 166   fPullHisto->GetAxis(6)->SetTitle("z_{mc} (cm)");
 167   fPullHisto->GetAxis(7)->SetTitle("#eta_{mc}");
 168   fPullHisto->GetAxis(8)->SetTitle("#phi_{mc} (rad)");
 169   fPullHisto->GetAxis(9)->SetTitle("p_{Tmc} (GeV/c)");
 170   fPullHisto->Sumw2();
 171   */
 172
 173   fPullHisto->GetAxis(9)->Set(nPtBins,binsPt);
 174
 175   fPullHisto->GetAxis(0)->SetTitle("(y-y_{mc})/#sigma");
 176   fPullHisto->GetAxis(1)->SetTitle("(z-z_{mc})/#sigma");
 177   fPullHisto->GetAxis(2)->SetTitle("(sin#phi-sin#phi_{mc})/#sigma");
 178   fPullHisto->GetAxis(3)->SetTitle("(tan#lambda-tan#lambda_{mc})/#sigma");
 179   fPullHisto->GetAxis(4)->SetTitle("(p_{Tmc}/p_{T}-1)/#sigma");
 180   fPullHisto->GetAxis(5)->SetTitle("y_{mc} (cm)");
 181   fPullHisto->GetAxis(6)->SetTitle("z_{mc} (cm)");
 182   fPullHisto->GetAxis(7)->SetTitle("sin#phi_{mc}");
 183   fPullHisto->GetAxis(8)->SetTitle("tan#lambda_{mc}");
 184   fPullHisto->GetAxis(9)->SetTitle("1/p_{Tmc} (GeV/c)^{-1}");
 185   fPullHisto->Sumw2();
 186
 187   // Init cuts
 188   if(!fCutsMC)
 189     AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
 190   if(!fCutsRC)
 191     AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
 192
 193   // init folder
 194   fAnalysisFolder = CreateFolder("folderRes","Analysis Resolution Folder");
 195 }
 196
 197 //_____________________________________________________________________________
 198 void AliPerformanceRes::ProcessTPC(AliStack* const stack, AliESDtrack *const esdTrack, AliESDEvent* const esdEvent)
 199 {
 200   if(!esdEvent) return;
 201   if(!esdTrack) return;
 202
 203   if( IsUseTrackVertex() )
 204   {
 205     // Relate TPC inner params to prim. vertex
 206     const AliESDVertex *vtxESD = esdEvent->GetPrimaryVertexTracks();
 207     Double_t x[3]; esdTrack->GetXYZ(x);
 208     Double_t b[3]; AliTracker::GetBxByBz(x,b);
 209     Bool_t isOK = esdTrack->RelateToVertexTPCBxByBz(vtxESD, b, kVeryBig);
 210     if(!isOK) return;
 211
 212     /*
 213       // JMT -- recaluclate DCA for HLT if not present
 214       if ( dca[0] == 0. && dca[1] == 0. ) {
 215         track->GetDZ( vtxESD->GetX(), vtxESD->GetY(), vtxESD->GetZ(), esdEvent->GetMagneticField(), dca );
 216       }
 217     */
 218   }
 219
 220   // Fill TPC only resolution comparison information
 221   const AliExternalTrackParam* tmpTrack = esdTrack->GetTPCInnerParam();
 222   if(!tmpTrack) return;
 223
 224   AliExternalTrackParam track = *tmpTrack;
 225
 226   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
 227   esdTrack->GetImpactParametersTPC(dca,cov);
 228
 229   //
 230   // Fill rec vs MC information
 231   //
 232   if(!stack) return;
 233   Int_t label = esdTrack->GetTPCLabel(); //Use TPC-only label for TPC-only resolution analysis
 234   if (label <= 0) return;
 235   TParticle* particle = stack->Particle(label);
 236   if(!particle) return;
 237   if(!particle->GetPDG()) return;
 238   if(particle->GetPDG()->Charge()==0) return;
 239   //printf("charge %d \n",particle->GetPDG()->Charge());
 240
 241   // Only 5 charged particle species (e,mu,pi,K,p)
 242   if (fCutsMC->IsPdgParticle(TMath::Abs(particle->GetPdgCode())) == kFALSE) return;
 243
 244   // exclude electrons
 245   if (fCutsMC->GetEM()==TMath::Abs(particle->GetPdgCode())) return;
 246
 247   Float_t deltaPtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC;
 248   Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC;
 249
 250   Float_t mceta =  particle->Eta();
 251   Float_t mcphi =  particle->Phi();
 252   if(mcphi<0) mcphi += 2.*TMath::Pi();
 253   Float_t mcpt = particle->Pt();
 254   Float_t mcsnp = TMath::Sin(TMath::ATan2(particle->Py(),particle->Px()));
 255   Float_t mctgl = TMath::Tan(TMath::ATan2(particle->Pz(),particle->Pt()));
 256
 257   // nb. TPC clusters cut
 258   if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;
 259
 260   // select primaries
 261   if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ())
 262   {
 263     if(mcpt == 0) return;
 264         double Bz = esdEvent->GetMagneticField();
 265
 266         Double_t mclocal[4]; //Rotated x,y,px,py mc-coordinates - the MC data should be rotated since the track is propagated best along x
 267         Double_t c = TMath::Cos(track.GetAlpha());
 268         Double_t s = TMath::Sin(track.GetAlpha());
 269         Double_t x = particle->Vx();
 270         Double_t y = particle->Vy();
 271         mclocal[0] = x*c + y*s;
 272         mclocal[1] =-x*s + y*c;
 273         Double_t px = particle->Px();
 274         Double_t py = particle->Py();
 275         mclocal[2] = px*c + py*s;
 276         mclocal[3] =-px*s + py*c;
 277     Float_t mcsnplocal = TMath::Sin(TMath::ATan2(mclocal[3],mclocal[2]));
 278
 279
 280         track.AliExternalTrackParam::PropagateTo(mclocal[0],Bz);
 281
 282     deltaYTPC= track.GetY()-mclocal[1];
 283     deltaZTPC = track.GetZ()-particle->Vz();
 284     deltaLambdaTPC = TMath::ATan2(track.Pz(),track.Pt())-TMath::ATan2(particle->Pz(),particle->Pt());
 285         //See comments in ProcessInnerTPC for remarks on local and global momentum coordinates for deltaPhi / pullSnp calculation
 286     deltaPhiTPC = TMath::ATan2(track.Py(),track.Px())-TMath::ATan2(particle->Py(),particle->Px());
 287     //delta1PtTPC = (track.OneOverPt()-1./mcpt)*mcpt;
 288     deltaPtTPC = (track.Pt()-mcpt) / mcpt;
 289
 290     pullYTPC= deltaYTPC / TMath::Sqrt(track.GetSigmaY2());
 291     pullZTPC = deltaZTPC / TMath::Sqrt(track.GetSigmaZ2());
 292
 293     //Double_t sigma_lambda = 1./(1.+track.GetTgl()*track.GetTgl()) * TMath::Sqrt(track.GetSigmaTgl2());
 294     //Double_t sigma_phi = 1./TMath::Sqrt(1-track.GetSnp()*track.GetSnp()) * TMath::Sqrt(track.GetSigmaSnp2());
 295     pullPhiTPC = (track.GetSnp() - mcsnplocal) / TMath::Sqrt(track.GetSigmaSnp2());
 296     pullLambdaTPC = (track.GetTgl() - mctgl) / TMath::Sqrt(track.GetSigmaTgl2());
 297
 298     //pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(track.GetSigmaTgl2());
 299     //pullPhiTPC = deltaPhiTPC / TMath::Sqrt(track.GetSigmaSnp2());
 300     if (mcpt) pull1PtTPC = (track.OneOverPt()-1./mcpt) / TMath::Sqrt(track.GetSigma1Pt2());
 301     else pull1PtTPC = 0.;
 302
 303     Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,deltaPtTPC,particle->Vy(),particle->Vz(),mcphi,mceta,mcpt};
 304     fResolHisto->Fill(vResolHisto);
 305
 306     Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,particle->Vy(),particle->Vz(),mcsnp,mctgl,1./mcpt};
 307     fPullHisto->Fill(vPullHisto);
 308   }
 309 }
 310
 311 //_____________________________________________________________________________
 312 void AliPerformanceRes::ProcessTPCITS(AliStack* const stack, AliESDtrack *const esdTrack, AliESDEvent* const esdEvent)
 313 {
 314   // Fill resolution comparison information (TPC+ITS)
 315   if(!esdEvent) return;
 316   if(!esdTrack) return;
 317
 318   if( IsUseTrackVertex() )
 319   {
 320     // Relate TPC inner params to prim. vertex
 321     const AliESDVertex *vtxESD = esdEvent->GetPrimaryVertexTracks();
 322     Double_t x[3]; esdTrack->GetXYZ(x);
 323     Double_t b[3]; AliTracker::GetBxByBz(x,b);
 324     Bool_t isOK = esdTrack->RelateToVertexBxByBz(vtxESD, b, kVeryBig);
 325     if(!isOK) return;
 326
 327     /*
 328       // JMT -- recaluclate DCA for HLT if not present
 329       if ( dca[0] == 0. && dca[1] == 0. ) {
 330         track->GetDZ( vtxESD->GetX(), vtxESD->GetY(), vtxESD->GetZ(), esdEvent->GetMagneticField(), dca );
 331       }
 332     */
 333   }
 334
 335   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
 336   esdTrack->GetImpactParameters(dca,cov);
 337
 338   //
 339   // Fill rec vs MC information
 340   //
 341   if(!stack) return;
 342
 343   Int_t label = TMath::Abs(esdTrack->GetLabel()); //Use global label for combined resolution analysis
 344   TParticle* particle = stack->Particle(label);
 345   if(!particle) return;
 346   if(!particle->GetPDG()) return;
 347   if(particle->GetPDG()->Charge()==0) return;
 348   //printf("charge %d \n",particle->GetPDG()->Charge());
 349
 350
 351   // Only 5 charged particle species (e,mu,pi,K,p)
 352   if (fCutsMC->IsPdgParticle(TMath::Abs(particle->GetPdgCode())) == kFALSE) return;
 353
 354   // exclude electrons
 355   if (fCutsMC->GetEM()==TMath::Abs(particle->GetPdgCode())) return;
 356
 357   Float_t mceta =  particle->Eta();
 358   Float_t mcphi =  particle->Phi();
 359   if(mcphi<0) mcphi += 2.*TMath::Pi();
 360   Float_t mcpt = particle->Pt();
 361   Float_t mcsnp = TMath::Sin(TMath::ATan2(particle->Py(),particle->Px()));
 362   Float_t mctgl = TMath::Tan(TMath::ATan2(particle->Pz(),particle->Pt()));
 363
 364   if ((esdTrack->GetStatus()&AliESDtrack::kTPCrefit)==0) return; // TPC refit
 365   if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; // min. nb. TPC clusters
 366   if(esdTrack->GetITSclusters(0)<fCutsRC->GetMinNClustersITS()) return;  // min. nb. ITS clusters
 367
 368   Float_t deltaPtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC;
 369   Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC;
 370
 371   // select primaries
 372   if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ())
 373   {
 374     if(mcpt == 0) return;
 375
 376     deltaYTPC= esdTrack->GetY()-particle->Vy();
 377     deltaZTPC = esdTrack->GetZ()-particle->Vz();
 378     deltaLambdaTPC = TMath::ATan2(esdTrack->Pz(),esdTrack->Pt())-TMath::ATan2(particle->Pz(),particle->Pt());
 379     deltaPhiTPC = TMath::ATan2(esdTrack->Py(),esdTrack->Px())-TMath::ATan2(particle->Py(),particle->Px());
 380     //delta1PtTPC = (esdTrack->OneOverPt()-1./mcpt)*mcpt;
 381     deltaPtTPC = (esdTrack->Pt()-mcpt) / mcpt;
 382
 383     pullYTPC= (esdTrack->GetY()-particle->Vy()) / TMath::Sqrt(esdTrack->GetSigmaY2());
 384     pullZTPC = (esdTrack->GetZ()-particle->Vz()) / TMath::Sqrt(esdTrack->GetSigmaZ2());
 385
 386     //Double_t sigma_lambda = 1./(1.+esdTrack->GetTgl()*esdTrack->GetTgl()) * TMath::Sqrt(esdTrack->GetSigmaTgl2());
 387     //Double_t sigma_phi = 1./TMath::Sqrt(1-esdTrack->GetSnp()*esdTrack->GetSnp()) * TMath::Sqrt(esdTrack->GetSigmaSnp2());
 388     pullPhiTPC = (esdTrack->GetSnp() - mcsnp) / TMath::Sqrt(esdTrack->GetSigmaSnp2());
 389     pullLambdaTPC = (esdTrack->GetTgl() - mctgl) / TMath::Sqrt(esdTrack->GetSigmaTgl2());
 390
 391     //pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(esdTrack->GetSigmaTgl2());
 392     //pullPhiTPC = deltaPhiTPC / TMath::Sqrt(esdTrack->GetSigmaSnp2());
 393     if (mcpt) pull1PtTPC = (esdTrack->OneOverPt()-1./mcpt) / TMath::Sqrt(esdTrack->GetSigma1Pt2());
 394     else pull1PtTPC = 0.;
 395
 396     Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,deltaPtTPC,particle->Vy(),particle->Vz(),mcphi,mceta,mcpt};
 397     fResolHisto->Fill(vResolHisto);
 398
 399     Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,particle->Vy(),particle->Vz(),mcsnp,mctgl,1./mcpt};
 400     fPullHisto->Fill(vPullHisto);
 401
 402
 403     /*
 404     deltaYTPC= esdTrack->GetY()-particle->Vy();
 405     deltaZTPC = esdTrack->GetZ()-particle->Vz();
 406     deltaLambdaTPC = TMath::ATan2(esdTrack->Pz(),esdTrack->Pt())-TMath::ATan2(particle->Pz(),particle->Pt());
 407     deltaPhiTPC = TMath::ATan2(esdTrack->Py(),esdTrack->Px())-TMath::ATan2(particle->Py(),particle->Px());
 408     delta1PtTPC = (esdTrack->OneOverPt()-1./mcpt)*mcpt;
 409
 410     pullYTPC= (esdTrack->GetY()-particle->Vy()) / TMath::Sqrt(esdTrack->GetSigmaY2());
 411     pullZTPC = (esdTrack->GetZ()-particle->Vz()) / TMath::Sqrt(esdTrack->GetSigmaZ2());
 412
 413     Double_t sigma_lambda = 1./(1.+esdTrack->GetTgl()*esdTrack->GetTgl()) * TMath::Sqrt(esdTrack->GetSigmaTgl2());
 414     Double_t sigma_phi = 1./TMath::Sqrt(1-esdTrack->GetSnp()*esdTrack->GetSnp()) * TMath::Sqrt(esdTrack->GetSigmaSnp2());
 415     pullLambdaTPC = deltaLambdaTPC / sigma_lambda;
 416     pullPhiTPC = deltaPhiTPC / sigma_phi;
 417
 418     //pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(esdTrack->GetSigmaTgl2());
 419     //pullPhiTPC = deltaPhiTPC / TMath::Sqrt(esdTrack->GetSigmaSnp2());
 420     if (mcpt) pull1PtTPC = (esdTrack->OneOverPt()-1./mcpt) / TMath::Sqrt(esdTrack->GetSigma1Pt2());
 421     else pull1PtTPC = 0.;
 422
 423     Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,delta1PtTPC,particle->Vy(),particle->Vz(),mceta,mcphi,mcpt};
 424     fResolHisto->Fill(vResolHisto);
 425
 426     Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,particle->Vy(),particle->Vz(),mceta,mcphi,mcpt};
 427     fPullHisto->Fill(vPullHisto);
 428     */
 429   }
 430 }
 431
 432 //_____________________________________________________________________________
 433 void AliPerformanceRes::ProcessConstrained(AliStack* const stack, AliESDtrack *const esdTrack, AliESDEvent* const esdEvent)
 434 {
 435   // Fill resolution comparison information (constarained parameters)
 436   //
 437   if(!esdEvent) return;
 438   if(!esdTrack) return;
 439
 440   if( IsUseTrackVertex() )
 441   {
 442     // Relate TPC inner params to prim. vertex
 443     const AliESDVertex *vtxESD = esdEvent->GetPrimaryVertexTracks();
 444     Double_t x[3]; esdTrack->GetXYZ(x);
 445     Double_t b[3]; AliTracker::GetBxByBz(x,b);
 446     Bool_t isOK = esdTrack->RelateToVertexBxByBz(vtxESD, b, kVeryBig);
 447     if(!isOK) return;
 448
 449     /*
 450       // JMT -- recaluclate DCA for HLT if not present
 451       if ( dca[0] == 0. && dca[1] == 0. ) {
 452         track->GetDZ( vtxESD->GetX(), vtxESD->GetY(), vtxESD->GetZ(), esdEvent->GetMagneticField(), dca );
 453       }
 454     */
 455   }
 456
 457
 458   const AliExternalTrackParam * track = esdTrack->GetConstrainedParam();
 459   if(!track) return;
 460
 461   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
 462   esdTrack->GetImpactParameters(dca,cov);
 463
 464   //
 465   // Fill rec vs MC information
 466   //
 467   if(!stack) return;
 468
 469   Int_t label = TMath::Abs(esdTrack->GetLabel());
 470   TParticle* particle = stack->Particle(label);
 471   if(!particle) return;
 472   if(!particle->GetPDG()) return;
 473   if(particle->GetPDG()->Charge()==0) return;
 474   //printf("charge %d \n",particle->GetPDG()->Charge());
 475
 476   // Only 5 charged particle species (e,mu,pi,K,p)
 477   if (fCutsMC->IsPdgParticle(TMath::Abs(particle->GetPdgCode())) == kFALSE) return;
 478
 479   // exclude electrons
 480   if (fCutsMC->GetEM()==TMath::Abs(particle->GetPdgCode())) return;
 481
 482   Float_t mceta =  particle->Eta();
 483   Float_t mcphi =  particle->Phi();
 484   if(mcphi<0) mcphi += 2.*TMath::Pi();
 485   Float_t mcpt = particle->Pt();
 486   Float_t mcsnp = TMath::Sin(TMath::ATan2(particle->Py(),particle->Px()));
 487   Float_t mctgl = TMath::Tan(TMath::ATan2(particle->Pz(),particle->Pt()));
 488
 489   if ((esdTrack->GetStatus()&AliESDtrack::kTPCrefit)==0) return; // TPC refit
 490   if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; // min. nb. TPC clusters
 491
 492   Float_t deltaPtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC;
 493   Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC;
 494
 495   // select primaries
 496   if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ())
 497   {
 498
 499     if(mcpt == 0) return;
 500
 501     deltaYTPC= track->GetY()-particle->Vy();
 502     deltaZTPC = track->GetZ()-particle->Vz();
 503     deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(particle->Pz(),particle->Pt());
 504     deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(particle->Py(),particle->Px());
 505     //delta1PtTPC = (track->OneOverPt()-1./mcpt)*mcpt;
 506     deltaPtTPC = (track->Pt()-mcpt) / mcpt;
 507
 508     pullYTPC= (track->GetY()-particle->Vy()) / TMath::Sqrt(track->GetSigmaY2());
 509     pullZTPC = (track->GetZ()-particle->Vz()) / TMath::Sqrt(track->GetSigmaZ2());
 510
 511     //Double_t sigma_lambda = 1./(1.+track->GetTgl()*track->GetTgl()) * TMath::Sqrt(track->GetSigmaTgl2());
 512     //Double_t sigma_phi = 1./TMath::Sqrt(1-track->GetSnp()*track->GetSnp()) * TMath::Sqrt(track->GetSigmaSnp2());
 513     pullPhiTPC = (track->GetSnp() - mcsnp) / TMath::Sqrt(track->GetSigmaSnp2());
 514     pullLambdaTPC = (track->GetTgl() - mctgl) / TMath::Sqrt(track->GetSigmaTgl2());
 515
 516     //pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(track->GetSigmaTgl2());
 517     //pullPhiTPC = deltaPhiTPC / TMath::Sqrt(track->GetSigmaSnp2());
 518     if (mcpt) pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2());
 519     else pull1PtTPC = 0.;
 520
 521     Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,deltaPtTPC,particle->Vy(),particle->Vz(),mcphi,mceta,mcpt};
 522     fResolHisto->Fill(vResolHisto);
 523
 524     Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,particle->Vy(),particle->Vz(),mcsnp,mctgl,1./mcpt};
 525     fPullHisto->Fill(vPullHisto);
 526
 527     /*
 528
 529     deltaYTPC= track->GetY()-particle->Vy();
 530     deltaZTPC = track->GetZ()-particle->Vz();
 531     deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(particle->Pz(),particle->Pt());
 532     deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(particle->Py(),particle->Px());
 533     delta1PtTPC = (track->OneOverPt()-1./mcpt)*mcpt;
 534
 535     pullYTPC= (track->GetY()-particle->Vy()) / TMath::Sqrt(track->GetSigmaY2());
 536     pullZTPC = (track->GetZ()-particle->Vz()) / TMath::Sqrt(track->GetSigmaZ2());
 537
 538     Double_t sigma_lambda = 1./(1.+track->GetTgl()*track->GetTgl()) * TMath::Sqrt(track->GetSigmaTgl2());
 539     Double_t sigma_phi = 1./TMath::Sqrt(1-track->GetSnp()*track->GetSnp()) * TMath::Sqrt(track->GetSigmaSnp2());
 540     pullLambdaTPC = deltaLambdaTPC / sigma_lambda;
 541     pullPhiTPC = deltaPhiTPC / sigma_phi;
 542
 543     //pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(track->GetSigmaTgl2());
 544     //pullPhiTPC = deltaPhiTPC / TMath::Sqrt(track->GetSigmaSnp2());
 545
 546     if (mcpt) pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2());
 547     else pull1PtTPC = 0.;
 548
 549     Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,delta1PtTPC,particle->Vy(),particle->Vz(),mceta,mcphi,mcpt};
 550     fResolHisto->Fill(vResolHisto);
 551
 552     Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,particle->Vy(),particle->Vz(),mceta,mcphi,mcpt};
 553     fPullHisto->Fill(vPullHisto);
 554
 555     */
 556   }
 557 }
 558
 559 //_____________________________________________________________________________
 560 void AliPerformanceRes::ProcessInnerTPC(AliMCEvent *const mcEvent, AliESDtrack *const esdTrack, AliESDEvent* const esdEvent)
 561 {
 562   //
 563   // Fill resolution comparison information (inner params at TPC reference point)
 564   //
 565   if(!esdEvent) return;
 566   if(!esdTrack) return;
 567
 568   if( IsUseTrackVertex() )
 569   {
 570     // Relate TPC inner params to prim. vertex
 571     const AliESDVertex *vtxESD = esdEvent->GetPrimaryVertexTracks();
 572     Double_t x[3]; esdTrack->GetXYZ(x);
 573     Double_t b[3]; AliTracker::GetBxByBz(x,b);
 574     Bool_t isOK = esdTrack->RelateToVertexTPCBxByBz(vtxESD, b, kVeryBig);
 575     if(!isOK) return;
 576
 577     /*
 578       // JMT -- recaluclate DCA for HLT if not present
 579       if ( dca[0] == 0. && dca[1] == 0. ) {
 580         track->GetDZ( vtxESD->GetX(), vtxESD->GetY(), vtxESD->GetZ(), esdEvent->GetMagneticField(), dca );
 581       }
 582     */
 583   }
 584
 585   const AliExternalTrackParam * innerParam = esdTrack->GetInnerParam();
 586   if(!innerParam) return;
 587
 588   // create new AliExternalTrackParam
 589   AliExternalTrackParam *track = new AliExternalTrackParam(*innerParam);
 590   if(!track) return;
 591
 592   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
 593   esdTrack->GetImpactParametersTPC(dca,cov);
 594
 595   //
 596   // Fill rec vs MC information
 597   //
 598   if(!mcEvent) return;
 599
 600   Int_t label = esdTrack->GetTPCLabel(); //Use TPC-only label for TPC-only resolution analysis
 601   if (label <= 0) return;
 602   AliMCParticle *mcParticle = (AliMCParticle*) mcEvent->GetTrack(label);
 603   if(!mcParticle) return;
 604
 605   // get the first TPC track reference
 606   AliTrackReference *ref0 = GetFirstTPCTrackRef(mcParticle);
 607   if(!ref0) return;
 608
 609   // Only 5 charged particle species (e,mu,pi,K,p)
 610   TParticle *particle = mcParticle->Particle();
 611   if(!particle) return;
 612
 613   if (fCutsMC->IsPdgParticle(TMath::Abs(particle->GetPdgCode())) == kFALSE) return;
 614
 615   // exclude electrons
 616   if (fCutsMC->GetEM()==TMath::Abs(particle->GetPdgCode())) return;
 617
 618   Double_t mclocal[4]; //Rotated x,y,Px,Py mc-coordinates - the MC data should be rotated since the track is propagated best along x
 619   Double_t c = TMath::Cos(track->GetAlpha());
 620   Double_t s = TMath::Sin(track->GetAlpha());
 621   Double_t x = ref0->X();
 622   Double_t y = ref0->Y();
 623   mclocal[0] = x*c + y*s;
 624   mclocal[1] =-x*s + y*c;
 625   Double_t px = ref0->Px();
 626   Double_t py = ref0->Py();
 627   mclocal[2] = px*c + py*s;
 628   mclocal[3] =-px*s + py*c;
 629
 630   //Double_t xyz[3] = {ref0->X(),ref0->Y(),ref0->Z()};
 631   // propagate track to the radius of the first track reference within TPC
 632   //Double_t trRadius = TMath::Sqrt(xyz[1] * xyz[1] + xyz[0] * xyz[0]);
 633   Double_t field[3]; track->GetBxByBz(field);
 634   if (TGeoGlobalMagField::Instance()->GetField() == NULL) {
 635     Error("ProcessInnerTPC", "Magnetic Field not set");
 636   }
 637   Bool_t isOK = track->PropagateToBxByBz(mclocal[0],field);
 638   if(!isOK) {return;}
 639   //track->AliExternalTrackParam::PropagateTo(mclocal[0],esdEvent->GetMagneticField());  //Use different propagation since above methods returns zero B field and does not work
 640
 641   Float_t mceta =  -TMath::Log(TMath::Tan(0.5 * ref0->Theta()));
 642   Float_t mcphi =  ref0->Phi();
 643   if(mcphi<0) mcphi += 2.*TMath::Pi();
 644   Float_t mcpt = ref0->Pt();
 645   Float_t mcsnp = TMath::Sin(TMath::ATan2(ref0->Py(),ref0->Px()));
 646   Float_t mcsnplocal = TMath::Sin(TMath::ATan2(mclocal[3],mclocal[2]));
 647   Float_t mctgl = TMath::Tan(TMath::ATan2(ref0->Pz(),ref0->Pt()));
 648
 649   if ((esdTrack->GetStatus()&AliESDtrack::kTPCrefit)==0) return; // TPC refit
 650   if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; // min. nb. TPC clusters
 651
 652   Float_t deltaPtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC;
 653   Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC;
 654
 655   // select primaries
 656   Bool_t isPrimary;
 657   if ( IsUseTrackVertex() )
 658   {
 659           isPrimary = TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ();
 660   }
 661   else
 662   {
 663           //If Track vertex is not used the above check does not work, hence we use the MC reference track
 664           isPrimary = label < mcEvent->Stack()->GetNprimary();
 665   }
 666   if(isPrimary)
 667   {
 668     if(mcpt == 0) return;
 669
 670     deltaYTPC= track->GetY()-mclocal[1];
 671     deltaZTPC = track->GetZ()-ref0->Z();
 672     deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(ref0->Pz(),ref0->Pt());
 673         //track->Px() etc returns momentum in global coordinates, hence mc momentum must not be rotated.
 674     deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(ref0->Py(),ref0->Px());
 675     //delta1PtTPC = (track->OneOverPt()-1./mcpt)*mcpt;
 676     deltaPtTPC = (track->Pt()-mcpt) / mcpt;
 677
 678     pullYTPC= deltaYTPC / TMath::Sqrt(track->GetSigmaY2());
 679     pullZTPC = deltaZTPC / TMath::Sqrt(track->GetSigmaZ2());
 680
 681     //Double_t sigma_lambda = 1./(1.+track->GetTgl()*track->GetTgl()) * TMath::Sqrt(track->GetSigmaTgl2());
 682     //Double_t sigma_phi = 1./TMath::Sqrt(1-track->GetSnp()*track->GetSnp()) * TMath::Sqrt(track->GetSigmaSnp2());
 683         //track->GetSnp returns value in local coordinates, hence here, in contrast to deltaPhiTPC, the mc momentum must be rotated
 684     pullPhiTPC = (track->GetSnp() - mcsnplocal) / TMath::Sqrt(track->GetSigmaSnp2());
 685     pullLambdaTPC = (track->GetTgl() - mctgl) / TMath::Sqrt(track->GetSigmaTgl2());
 686
 687     //pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(track->GetSigmaTgl2());
 688     //pullPhiTPC = deltaPhiTPC / TMath::Sqrt(track->GetSigmaSnp2());
 689     if (mcpt) pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2());
 690     else pull1PtTPC = 0.;
 691
 692     Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,deltaPtTPC,ref0->Y(),ref0->Z(),mcphi,mceta,mcpt};
 693     fResolHisto->Fill(vResolHisto);
 694
 695     Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,ref0->Y(),ref0->Z(),mcsnp,mctgl,1./mcpt};
 696     fPullHisto->Fill(vPullHisto);
 697   }
 698
 699   if(track) delete track;
 700 }
 701
 702 //_____________________________________________________________________________
 703 void AliPerformanceRes::ProcessOuterTPC(AliMCEvent *const mcEvent, AliESDtrack *const esdTrack, AliESDfriendTrack *const friendTrack, AliESDEvent* const esdEvent)
 704 {
 705   //
 706   // Fill resolution comparison information (outer params at TPC reference point)
 707   //
 708   if(!friendTrack) return;
 709   if(!esdEvent) return;
 710   if(!esdTrack) return;
 711
 712   if( IsUseTrackVertex() )
 713   {
 714     // Relate TPC inner params to prim. vertex
 715     const AliESDVertex *vtxESD = esdEvent->GetPrimaryVertexTracks();
 716     Double_t x[3]; esdTrack->GetXYZ(x);
 717     Double_t b[3]; AliTracker::GetBxByBz(x,b);
 718     Bool_t isOK = esdTrack->RelateToVertexTPCBxByBz(vtxESD, b, kVeryBig);
 719     if(!isOK) return;
 720
 721     /*
 722       // JMT -- recaluclate DCA for HLT if not present
 723       if ( dca[0] == 0. && dca[1] == 0. ) {
 724         track->GetDZ( vtxESD->GetX(), vtxESD->GetY(), vtxESD->GetZ(), esdEvent->GetMagneticField(), dca );
 725       }
 726     */
 727   }
 728
 729   const AliExternalTrackParam * outerParam = friendTrack->GetTPCOut();
 730   if(!outerParam) return;
 731
 732   // create new AliExternalTrackParam
 733   AliExternalTrackParam *track = new AliExternalTrackParam(*outerParam);
 734   if(!track) return;
 735
 736   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
 737   esdTrack->GetImpactParametersTPC(dca,cov);
 738
 739   //
 740   // Fill rec vs MC information
 741   //
 742   if(!mcEvent) return;
 743
 744   Int_t label = esdTrack->GetTPCLabel(); //Use TPC-only label for TPC-only resolution analysis
 745   if (label <= 0) return;
 746   AliMCParticle *mcParticle = (AliMCParticle*) mcEvent->GetTrack(label);
 747   if(!mcParticle) return;
 748
 749   // get the last TPC track reference
 750   AliTrackReference *ref0 = GetLastTPCTrackRef(mcParticle);
 751   if(!ref0) return;
 752
 753   // Only 5 charged particle species (e,mu,pi,K,p)
 754   TParticle *particle = mcParticle->Particle();
 755   if(!particle) return;
 756
 757   if (fCutsMC->IsPdgParticle(TMath::Abs(particle->GetPdgCode())) == kFALSE) return;
 758
 759   // exclude electrons
 760   if (fCutsMC->GetEM()==TMath::Abs(particle->GetPdgCode())) return;
 761
 762   // calculate alpha angle
 763   Double_t xyz[3] = {ref0->X(),ref0->Y(),ref0->Z()};
 764   Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
 765   //if (alpha<0) alpha += TMath::TwoPi();
 766
 767   // rotate outer track to local coordinate system
 768   // and propagate to the radius of the last track reference of TPC
 769   Double_t trRadius = TMath::Sqrt(xyz[1] * xyz[1] + xyz[0] * xyz[0]);
 770   //Bool_t isOK = track->Propagate(alpha,trRadius,AliTracker::GetBz());
 771   Double_t field[3]; track->GetBxByBz(field);
 772   Bool_t isOK = track->PropagateBxByBz(alpha,trRadius,field);
 773   if(!isOK) return;
 774
 775   Float_t mceta =  -TMath::Log(TMath::Tan(0.5 * ref0->Theta()));
 776   Float_t mcphi =  ref0->Phi();
 777   if(mcphi<0) mcphi += 2.*TMath::Pi();
 778   Float_t mcpt = ref0->Pt();
 779   Float_t mcsnp = TMath::Sin(TMath::ATan2(ref0->Py(),ref0->Px()));
 780   Float_t mctgl = TMath::Tan(TMath::ATan2(ref0->Pz(),ref0->Pt()));
 781
 782   if ((esdTrack->GetStatus()&AliESDtrack::kTPCrefit)==0) return; // TPC refit
 783   if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; // min. nb. TPC clusters
 784
 785   Float_t deltaPtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC;
 786   Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC;
 787
 788   // select primaries
 789   if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ())
 790   {
 791     if(mcpt == 0) return;
 792
 793     deltaYTPC= track->GetY(); // already rotated
 794     deltaZTPC = track->GetZ()-ref0->Z();
 795     deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(ref0->Pz(),ref0->Pt());
 796     deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(ref0->Py(),ref0->Px());
 797     //delta1PtTPC = (track->OneOverPt()-1./mcpt)*mcpt;
 798     deltaPtTPC = (track->Pt()-mcpt) / mcpt;
 799
 800     pullYTPC= track->GetY() / TMath::Sqrt(track->GetSigmaY2());
 801     pullZTPC = (track->GetZ()-ref0->Z()) / TMath::Sqrt(track->GetSigmaZ2());
 802
 803     //Double_t sigma_lambda = 1./(1.+track->GetTgl()*track->GetTgl()) * TMath::Sqrt(track->GetSigmaTgl2());
 804     //Double_t sigma_phi = 1./TMath::Sqrt(1-track->GetSnp()*track->GetSnp()) * TMath::Sqrt(track->GetSigmaSnp2());
 805     pullPhiTPC = (track->GetSnp() - mcsnp) / TMath::Sqrt(track->GetSigmaSnp2());
 806     pullLambdaTPC = (track->GetTgl() - mctgl) / TMath::Sqrt(track->GetSigmaTgl2());
 807
 808     //pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(track->GetSigmaTgl2());
 809     //pullPhiTPC = deltaPhiTPC / TMath::Sqrt(track->GetSigmaSnp2());
 810     if (mcpt) pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2());
 811     else pull1PtTPC = 0.;
 812
 813     Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,deltaPtTPC,ref0->Y(),ref0->Z(),mcphi,mceta,mcpt};
 814     fResolHisto->Fill(vResolHisto);
 815
 816     Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,ref0->Y(),ref0->Z(),mcsnp,mctgl,1./mcpt};
 817     fPullHisto->Fill(vPullHisto);
 818   }
 819
 820   if(track) delete track;
 821 }
 822
 823 //_____________________________________________________________________________
 824 AliTrackReference * AliPerformanceRes::GetFirstTPCTrackRef(AliMCParticle *mcParticle)
 825 {
 826   // return first TPC track reference
 827   if(!mcParticle) return 0;
 828
 829   // find first track reference
 830   // check direction to select proper reference point for looping tracks
 831   Int_t nTrackRef = mcParticle->GetNumberOfTrackReferences();
 832   AliTrackReference *ref = 0;
 833   AliTrackReference *refIn = 0;
 834   for (Int_t iref = 0; iref < nTrackRef; iref++) {
 835     ref = mcParticle->GetTrackReference(iref);
 836     if(ref && (ref->DetectorId()==AliTrackReference::kTPC))
 837     {
 838       Float_t dir = ref->X()*ref->Px()+ref->Y()*ref->Py();
 839       if(dir < 0.) break;
 840
 841       refIn = ref;
 842       break;
 843     }
 844   }
 845
 846 return refIn;
 847 }
 848
 849 //_____________________________________________________________________________
 850 AliTrackReference * AliPerformanceRes::GetLastTPCTrackRef(AliMCParticle *mcParticle)
 851 {
 852   // return last TPC track reference
 853   if(!mcParticle) return 0;
 854
 855   // find last track reference
 856   // check direction to select proper reference point for looping tracks
 857   Int_t nTrackRef = mcParticle->GetNumberOfTrackReferences();
 858   AliTrackReference *ref = 0;
 859   AliTrackReference *refOut = 0;
 860   for (Int_t iref = 0; iref < nTrackRef; iref++) {
 861     ref = mcParticle->GetTrackReference(iref);
 862     if(ref && (ref->DetectorId()==AliTrackReference::kTPC)) {
 863       Float_t dir=ref->X()*ref->Px()+ref->Y()*ref->Py();
 864       if(dir< 0.0) break;
 865       refOut = ref;
 866     }
 867   }
 868
 869 return refOut;
 870 }
 871
 872 //_____________________________________________________________________________
 873 void AliPerformanceRes::Exec(AliMCEvent* const mcEvent, AliESDEvent *const esdEvent, AliESDfriend *const esdFriend, const Bool_t bUseMC, const Bool_t bUseESDfriend)
 874 {
 875   // Process comparison information
 876   if(!esdEvent)
 877   {
 878     Error("Exec","esdEvent not available");
 879     return;
 880   }
 881   AliHeader* header = 0;
 882   AliGenEventHeader* genHeader = 0;
 883   AliStack* stack = 0;
 884   TArrayF vtxMC(3);
 885
 886   if(bUseMC)
 887   {
 888     if(!mcEvent) {
 889       Error("Exec","mcEvent not available");
 890       return;
 891     }
 892     // get MC event header
 893     header = mcEvent->Header();
 894     if (!header) {
 895       Error("Exec","Header not available");
 896       return;
 897     }
 898     // MC particle stack
 899     stack = mcEvent->Stack();
 900     if (!stack) {
 901       Error("Exec","Stack not available");
 902       return;
 903     }
 904     // get MC vertex
 905     genHeader = header->GenEventHeader();
 906     if (!genHeader) {
 907       Error("Exec","Could not retrieve genHeader from Header");
 908       return;
 909     }
 910     genHeader->PrimaryVertex(vtxMC);
 911   }
 912   else {
 913     Error("Exec","MC information required!");
 914     return;
 915   }
 916
 917   // use ESD friends
 918   if(bUseESDfriend) {
 919     if(!esdFriend) {
 920       Error("Exec","esdFriend not available");
 921       return;
 922     }
 923   }
 924
 925   // get event vertex
 926   const AliESDVertex *vtxESD = NULL;
 927   if( IsUseTrackVertex() )
 928   {
 929     // track vertex
 930     vtxESD = esdEvent->GetPrimaryVertexTracks();
 931         if(vtxESD && (vtxESD->GetStatus()<=0)) return;
 932   }
 933   // Coverity - removed else branch as vtxESD is not further used in method
 934   //  else {
 935   //    // TPC track vertex
 936   //    vtxESD = esdEvent->GetPrimaryVertexTPC();
 937   //  }
 938
 939
 940
 941
 942   //  Process events
 943   for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
 944   {
 945     AliESDtrack *track = esdEvent->GetTrack(iTrack);
 946     if(!track) continue;
 947
 948     AliESDfriendTrack *friendTrack=0;
 949
 950
 951     Int_t label = TMath::Abs(track->GetLabel());
 952     if ( label > stack->GetNtrack() )
 953     {
 954       ULong_t status = track->GetStatus();
 955       printf ("Error : ESD MCLabel %d - StackSize %d - Status %lu \n",
 956                track->GetLabel(), stack->GetNtrack(), status );
 957       printf(" NCluster %d \n", track->GetTPCclusters(0) );
 958       /*
 959       if ((status&AliESDtrack::kTPCrefit)== 0 ) printf("   kTPCrefit \n");
 960       if ((status&AliESDtrack::kTPCin)== 0 )    printf("   kTPCin \n");
 961       if ((status&AliESDtrack::kTPCout)== 0 )   printf("   kTPCout \n");
 962       if ((status&AliESDtrack::kTRDrefit)== 0 ) printf("   kTRDrefit \n");
 963       if ((status&AliESDtrack::kTRDin)== 0 )    printf("   kTRDin \n");
 964       if ((status&AliESDtrack::kTRDout)== 0 )   printf("   kTRDout \n");
 965       if ((status&AliESDtrack::kITSrefit)== 0 ) printf("   kITSrefit \n");
 966       if ((status&AliESDtrack::kITSin)== 0 )    printf("   kITSin \n");
 967       if ((status&AliESDtrack::kITSout)== 0 )   printf("   kITSout \n");
 968       */
 969
 970       continue;
 971     }
 972
 973         if (label == 0) continue;               //Cannot distinguish between track or fake track
 974         if (track->GetLabel() < 0) continue; //Do not consider fake tracks
 975
 976     if(GetAnalysisMode() == 0) ProcessTPC(stack,track,esdEvent);
 977     else if(GetAnalysisMode() == 1) ProcessTPCITS(stack,track,esdEvent);
 978     else if(GetAnalysisMode() == 2) ProcessConstrained(stack,track,esdEvent);
 979     else if(GetAnalysisMode() == 3) ProcessInnerTPC(mcEvent,track,esdEvent);
 980     else if(GetAnalysisMode() == 4) {
 981
 982     if(bUseESDfriend) {
 983       friendTrack=esdFriend->GetTrack(iTrack);
 984       if(!friendTrack) continue;
 985     }
 986
 987         ProcessOuterTPC(mcEvent,track,friendTrack,esdEvent);
 988         }
 989     else {
 990       printf("ERROR: AnalysisMode %d \n",fAnalysisMode);
 991       return;
 992     }
 993   }
 994 }
 995
 996 //_____________________________________________________________________________
 997 TH1F* AliPerformanceRes::MakeResol(TH2F * his, Int_t integ, Bool_t type, Int_t cut){
 998   // Create resolution histograms
 999
1000   TH1F *hisr, *hism;
1001   if (!gPad) new TCanvas;
1002   hisr = AliTreeDraw::CreateResHistoII(his,&hism,integ,kTRUE,cut);
1003   if (type) return hism;
1004   else
1005     return hisr;
1006 }
1007
1008 //_____________________________________________________________________________
1009 void AliPerformanceRes::Analyse() {
1010   // Analyse comparison information and store output histograms
1011   // in the folder "folderRes"
1012   //
1013   TH1::AddDirectory(kFALSE);
1014   TH1F *h=0;
1015   TH2F *h2D=0;
1016   TObjArray *aFolderObj = new TObjArray;
1017   if(!aFolderObj) return;
1018
1019   // write results in the folder
1020   TCanvas * c = new TCanvas("Phi resol Tan","Phi resol Tan");
1021   c->cd();
1022
1023   char name[256];
1024   char title[256];
1025
1026   for(Int_t i=0; i<5; i++)
1027   {
1028     for(Int_t j=5; j<10; j++)
1029     {
1030       if(j!=8) fResolHisto->GetAxis(8)->SetRangeUser(-0.9,0.89); // eta window
1031       //if(j!=8) fResolHisto->GetAxis(8)->SetRangeUser(0.0,0.89); // eta window
1032       else fResolHisto->GetAxis(8)->SetRangeUser(-1.5,1.49);
1033       if(GetAnalysisMode() == 3) fResolHisto->GetAxis(5)->SetRangeUser(-80.,79.99); // y range
1034
1035       h2D = (TH2F*)fResolHisto->Projection(i,j);
1036
1037       h = AliPerformanceRes::MakeResol(h2D,1,0,100);
1038       snprintf(name,256,"h_res_%d_vs_%d",i,j);
1039       h->SetName(name);
1040
1041       h->GetXaxis()->SetTitle(fResolHisto->GetAxis(j)->GetTitle());
1042       snprintf(title,256,"%s %s",fResolHisto->GetAxis(i)->GetTitle(),"(resolution)");
1043       h->GetYaxis()->SetTitle(title);
1044       snprintf(title,256,"%s vs %s",title,fResolHisto->GetAxis(j)->GetTitle());
1045       h->SetTitle(title);
1046
1047       if(j==9) h->SetBit(TH1::kLogX);
1048       aFolderObj->Add(h);
1049
1050       h = AliPerformanceRes::MakeResol(h2D,1,1,100);
1051       //h = (TH1F*)arr->At(1);
1052       snprintf(name,256,"h_mean_res_%d_vs_%d",i,j);
1053       h->SetName(name);
1054
1055       h->GetXaxis()->SetTitle(fResolHisto->GetAxis(j)->GetTitle());
1056       snprintf(title,256,"%s %s",fResolHisto->GetAxis(i)->GetTitle(),"(mean)");
1057       h->GetYaxis()->SetTitle(title);
1058
1059       snprintf(title,256,"%s vs %s",title,fResolHisto->GetAxis(j)->GetTitle());
1060       h->SetTitle(title);
1061
1062       if(j==9) h->SetBit(TH1::kLogX);
1063       aFolderObj->Add(h);
1064
1065       fResolHisto->GetAxis(8)->SetRangeUser(-1.5,1.5);
1066       fResolHisto->GetAxis(9)->SetRangeUser(0.1,100.);
1067
1068       //
1069       if(j!=8) fPullHisto->GetAxis(8)->SetRangeUser(-0.9,0.89); // eta window
1070       //if(j!=8) fPullHisto->GetAxis(8)->SetRangeUser(0.0,0.89);    // eta window
1071       else  fPullHisto->GetAxis(8)->SetRangeUser(-1.5,1.49);      // eta window
1072       fPullHisto->GetAxis(9)->SetRangeUser(0.,9.99);            // 1./pt threshold
1073       if(GetAnalysisMode() == 3) fPullHisto->GetAxis(5)->SetRangeUser(-80.,79.99); // y range
1074
1075       h2D = (TH2F*)fPullHisto->Projection(i,j);
1076
1077       h = AliPerformanceRes::MakeResol(h2D,1,0,100);
1078       snprintf(name,256,"h_pull_%d_vs_%d",i,j);
1079       h->SetName(name);
1080
1081       h->GetXaxis()->SetTitle(fPullHisto->GetAxis(j)->GetTitle());
1082       snprintf(title,256,"%s %s",fPullHisto->GetAxis(i)->GetTitle(),"(resolution)");
1083       h->GetYaxis()->SetTitle(title);
1084       snprintf(title,256,"%s vs %s",title,fPullHisto->GetAxis(j)->GetTitle());
1085       h->SetTitle(title);
1086
1087       //if(j==9) h->SetBit(TH1::kLogX);
1088       aFolderObj->Add(h);
1089
1090       h = AliPerformanceRes::MakeResol(h2D,1,1,100);
1091       snprintf(name,256,"h_mean_pull_%d_vs_%d",i,j);
1092       h->SetName(name);
1093
1094       h->GetXaxis()->SetTitle(fPullHisto->GetAxis(j)->GetTitle());
1095       snprintf(title,256,"%s %s",fPullHisto->GetAxis(i)->GetTitle(),"(mean)");
1096       h->GetYaxis()->SetTitle(title);
1097       snprintf(title,256,"%s vs %s",title,fPullHisto->GetAxis(j)->GetTitle());
1098       h->SetTitle(title);
1099
1100       //if(j==9) h->SetBit(TH1::kLogX);
1101       aFolderObj->Add(h);
1102     }
1103   }
1104
1105   for (Int_t i = 0;i < fResolHisto->GetNdimensions();i++)
1106   {
1107           fResolHisto->GetAxis(i)->SetRange(1,0);                               //Reset Range
1108   }
1109   for (Int_t i = 0;i < fPullHisto->GetNdimensions();i++)
1110   {
1111           fPullHisto->GetAxis(i)->SetRange(1,0);                                //Reset Range
1112   }
1113
1114   // export objects to analysis folder
1115   fAnalysisFolder = ExportToFolder(aFolderObj);
1116
1117   // delete only TObjArray
1118   if(aFolderObj) delete aFolderObj;
1119 }
1120
1121 //_____________________________________________________________________________
1122 TFolder* AliPerformanceRes::ExportToFolder(TObjArray * array)
1123 {
1124   // recreate folder avery time and export objects to new one
1125   //
1126   AliPerformanceRes * comp=this;
1127   TFolder *folder = comp->GetAnalysisFolder();
1128
1129   TString name, title;
1130   TFolder *newFolder = 0;
1131   Int_t i = 0;
1132   Int_t size = array->GetSize();
1133
1134   if(folder) {
1135      // get name and title from old folder
1136      name = folder->GetName();
1137      title = folder->GetTitle();
1138
1139          // delete old one
1140      delete folder;
1141
1142          // create new one
1143      newFolder = CreateFolder(name.Data(),title.Data());
1144      newFolder->SetOwner();
1145
1146          // add objects to folder
1147      while(i < size) {
1148            newFolder->Add(array->At(i));
1149            i++;
1150          }
1151   }
1152
1153 return newFolder;
1154 }
1155
1156 //_____________________________________________________________________________
1157 Long64_t AliPerformanceRes::Merge(TCollection* const list)
1158 {
1159   // Merge list of objects (needed by PROOF)
1160
1161   if (!list)
1162   return 0;
1163
1164   if (list->IsEmpty())
1165   return 1;
1166
1167   TIterator* iter = list->MakeIterator();
1168   TObject* obj = 0;
1169
1170   // collection of generated histograms
1171   Int_t count=0;
1172   while((obj = iter->Next()) != 0)
1173   {
1174   AliPerformanceRes* entry = dynamic_cast<AliPerformanceRes*>(obj);
1175   if (entry == 0) continue;
1176   if (fResolHisto->GetEntries()<fgkMergeEntriesCut){
1177     fResolHisto->Add(entry->fResolHisto);
1178     fPullHisto->Add(entry->fPullHisto);
1179   }
1180
1181   count++;
1182   }
1183
1184 return count;
1185 }
1186
1187 //_____________________________________________________________________________
1188 TFolder* AliPerformanceRes::CreateFolder(TString name,TString title) {
1189 // create folder for analysed histograms
1190 //
1191 TFolder *folder = 0;
1192   folder = new TFolder(name.Data(),title.Data());
1193
1194   return folder;
1195 }