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