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