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