PWGPP/EvTrkSelection/AliAnalysisTrackingUncertainties.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, proviyaded that the above copyright notice appears in all *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purapose. It is         *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //////////////////////////////////////////////////////////////////////////////
  17 //                                                                          //
  18 // Analysis task for the systematic study of the uncertainties related to   //
  19 // the tracking and ITS-TPC matching efficiency for different particle      //
  20 // species.                                                                 //
  21 //                                                                          //
  22 //////////////////////////////////////////////////////////////////////////////
  23
  24
  25 #include "Riostream.h"
  26 #include "TH1F.h"
  27 #include "TH2D.h"
  28 #include "TH3F.h"
  29 #include "THn.h"
  30 #include "TList.h"
  31 #include "TMath.h"
  32 #include "TParticlePDG.h"
  33 //
  34 #include "AliAnalysisTaskSE.h"
  35 #include "AliAnalysisManager.h"
  36 #include "AliPID.h"
  37 #include "AliESDtrackCuts.h"
  38 #include "AliESDVertex.h"
  39 #include "AliESDEvent.h"
  40 #include "AliESDInputHandler.h"
  41 #include "AliESDtrack.h"
  42 #include "AliESDpid.h"
  43 #include "AliESDUtils.h"
  44 #include "AliMCEventHandler.h"
  45 #include "AliMCEvent.h"
  46 #include "AliStack.h"
  47 #include "AliLog.h"
  48 //
  49 #include "AliAnalysisTrackingUncertainties.h"
  50
  51
  52 ClassImp(AliAnalysisTrackingUncertainties)
  53
  54 // histogram constants
  55 const Int_t kNumberOfAxes = 5;
  56
  57 //________________________________________________________________________
  58 AliAnalysisTrackingUncertainties::AliAnalysisTrackingUncertainties()
  59   : AliAnalysisTaskSE("TaskTestPA"),
  60     fESD(0),
  61     fESDpid(0),
  62     fUtils(0),
  63     fMCtrue(0),
  64     fListHist(0),
  65     fESDtrackCuts(0),
  66     fMatchTr(),
  67     fMatchChi(),
  68     fExcludeMomFromChi2ITSTPC(0)
  69 {
  70   // default Constructor
  71   /* fast compilation test
  72
  73      gSystem->Load("libANALYSIS");
  74      gSystem->Load("libANALYSISalice");
  75      .L AliAnalysisTrackingUncertainties.cxx++
  76    */
  77 }
  78
  79
  80 //________________________________________________________________________
  81 AliAnalysisTrackingUncertainties::AliAnalysisTrackingUncertainties(const char *name)
  82   : AliAnalysisTaskSE(name),
  83     fESD(0),
  84     fESDpid(0),
  85     fUtils(0),
  86     fMCtrue(0),
  87     fListHist(0),
  88     fESDtrackCuts(0),
  89     fMatchTr(),
  90     fMatchChi(),
  91     fExcludeMomFromChi2ITSTPC(0)
  92 {
  93   //
  94   // standard constructur which should be used
  95   //
  96   fMCtrue = kTRUE;
  97   //
  98   // create track cuts
  99   //
 100   fESDtrackCuts = new AliESDtrackCuts("AliESDtrackCuts","AliESDtrackCuts");
 101   fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
 102   fESDtrackCuts->SetEtaRange(-1., 1.);
 103   //
 104   // analysis utils if needed
 105   //
 106   fUtils = new AliAnalysisUtils();
 107   //
 108
 109   // Output slot #0 writes into a TList container
 110   DefineOutput(1, TList::Class());
 111
 112 }
 113
 114
 115 //________________________________________________________________________
 116 void AliAnalysisTrackingUncertainties::UserCreateOutputObjects()
 117 {
 118   //
 119   // Create histograms
 120   // Called once
 121   //
 122   fListHist = new TList();
 123   fListHist->SetOwner(kTRUE);
 124   //
 125   // (1.) basic QA and statistics histograms
 126   //
 127   TH2F * histVertexSelection = new TH2F("histVertexSelection", "vertex selection; vertex z (cm); accepted/rejected", 100, -50., 50., 2, -0.5, 1.5);
 128   fListHist->Add(histVertexSelection);
 129   //
 130   // (2.) track cut variation histograms
 131   //
 132   InitializeTrackCutHistograms();
 133   //
 134   // (3.) ITS -> TPC matching histograms
 135   //
 136   //  Int_t    binsMatch[kNumberOfAxes] = { 10,   50,    20,            18,  6};
 137   //  Double_t minMatch[kNumberOfAxes]  = {  0,  0.1,    -1,             0, -0.5};
 138   //  Double_t maxMatch[kNumberOfAxes]  = {200,   20,    +1, 2*TMath::Pi(),  5.5};
 139   //
 140   //  TString axisNameMatch[kNumberOfAxes]  = {"matchChi2","pT","eta","phi","pid"};
 141   //  TString axisTitleMatch[kNumberOfAxes] = {"matchChi2","pT","eta","phi","pid"};
 142   //
 143   //  THnF * hBestMatch = new THnF("hBestMatch","ITS -> TPC matching ",kNumberOfAxes, binsMatch, minMatch, maxMatch);
 144   //  for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
 145   //    hBestMatch->GetAxis(iaxis)->SetName(axisNameMatch[iaxis]);
 146   // hBestMatch->GetAxis(iaxis)->SetTitle(axisTitleMatch[iaxis]);
 147   //  }
 148   //  BinLogAxis(hBestMatch, 1);
 149   //  fListHist->Add(hBestMatch);
 150   //
 151   //
 152   //
 153   const int nbPt=40;
 154   const double ptMax=5;
 155   //
 156   TH2F * hNMatch    = new TH2F("hNMatch","N Matches",nbPt,0,ptMax,kMaxMatch+1,-0.5,kMaxMatch+0.5);
 157   TH2F * hBestMatch = new TH2F("hBestMatch","Best Match Chi2",nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2); // OB
 158   TH2F * hBestMatch_cuts = new TH2F("hBestMatch_cuts","Best Match Chi2",nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2);
 159   TH2F * hAllMatch  = new TH2F("hAllMatch","All Matches Chi2",nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2);
 160   TH2F * hAllMatchGlo  = new TH2F("hAllMatchGlo","All Matches Chi2",nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2);
 161   TH2F * hPtCorr_ITSTPC = new TH2F("hPtCorr_ITSTPC","PtCorr",nbPt,0,ptMax,nbPt,0,ptMax);
 162   TH2F * hdPtRel_ITSTPC = new TH2F("hdPtRel_ITSTPC","dPt/pt",nbPt,0,ptMax,2*nbPt+1,-0.4*ptMax,0.4*ptMax);
 163   TH2F * hdInvPtRel_ITSTPC = new TH2F("hdInvPtRel_ITSTPC","pt*dPt^{-1}",nbPt,0,ptMax,2*nbPt+1,-0.4*ptMax,0.4*ptMax);
 164   //
 165   fListHist->Add(hNMatch);
 166   fListHist->Add(hBestMatch);
 167   fListHist->Add(hBestMatch_cuts);
 168   fListHist->Add(hAllMatch);
 169   fListHist->Add(hAllMatchGlo);
 170   fListHist->Add(hPtCorr_ITSTPC);
 171   fListHist->Add(hdPtRel_ITSTPC);
 172   fListHist->Add(hdInvPtRel_ITSTPC);
 173   //
 174   //
 175   TH2F * hNMatchBg    = new TH2F("hNMatchBg","N Matches",nbPt,0,ptMax,kMaxMatch+1,-0.5,kMaxMatch+0.5);
 176   TH2F * hBestMatchBg = new TH2F("hBestMatchBg","Best Match Chi2",nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2);
 177   TH2F * hBestMatchBg_cuts = new TH2F("hBestMatchBg_cuts","Best Match Chi2",nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2); // OB
 178   TH2F * hAllMatchBg  = new TH2F("hAllMatchBg","All Matches Chi2",nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2);
 179   TH2F * hAllMatchGloBg  = new TH2F("hAllMatchGloBg","All Matches Chi2",nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2);
 180   TH2F * hdPtRelBg_ITSTPC = new TH2F("hdPtRelBg_ITSTPC","dPt/pt",nbPt,0,ptMax,2*nbPt+1,-0.4*ptMax,0.4*ptMax);
 181   TH2F * hdInvPtRelBg_ITSTPC = new TH2F("hdInvPtRelBg_ITSTPC","pt*dPt^{-1}",nbPt,0,ptMax,2*nbPt+1,-0.4*ptMax,0.4*ptMax);
 182
 183   //cout<<" here 0 : hdPtRelBg_ITSTPC "<<hdPtRelBg_ITSTPC<<" hdInvPtRelBg_ITSTPC "<<hdInvPtRelBg_ITSTPC<<endl;
 184
 185   fListHist->Add(hNMatchBg);
 186   fListHist->Add(hBestMatchBg);
 187   fListHist->Add(hBestMatchBg_cuts);
 188   fListHist->Add(hAllMatchBg);
 189   fListHist->Add(hAllMatchGloBg);
 190    fListHist->Add(hdPtRelBg_ITSTPC);
 191   fListHist->Add(hdInvPtRelBg_ITSTPC);
 192   //add default track cuts in the output list
 193   fListHist->Add(fESDtrackCuts);
 194   //
 195   // post data
 196   //
 197   PostData(1, fListHist);
 198
 199
 200 }
 201
 202
 203
 204 //________________________________________________________________________
 205 void AliAnalysisTrackingUncertainties::UserExec(Option_t *)
 206 {
 207   //
 208   // main event loop
 209   //
 210   fESD = dynamic_cast<AliESDEvent*>( InputEvent() );
 211   if (!fESD) {
 212     PostData(1, fListHist);
 213     return;
 214   }
 215   //
 216   if (!fESDpid) fESDpid = ((AliESDInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetESDpid();
 217   //
 218   // Check Monte Carlo information and other access first:
 219   //
 220   AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 221   if (!eventHandler) {
 222     fMCtrue = kFALSE;
 223   }
 224   //
 225   // extract generated particles information
 226   //
 227   AliMCEvent* mcEvent = 0x0;
 228   AliStack* stack = 0x0;
 229   if (eventHandler) mcEvent = eventHandler->MCEvent();
 230   if (!mcEvent) {
 231     if (fMCtrue) return;
 232   }
 233   if (fMCtrue) {
 234     //
 235     stack = mcEvent->Stack();
 236     if (!stack) return;
 237     //
 238     for(Int_t i = 0; i < stack->GetNtrack(); i++) {
 239       /* at the moment nothing is needed here
 240       TParticle * trackMC = stack->Particle(i);
 241       Double_t rap = trackMC->Eta();
 242       Double_t y = trackMC->Y();
 243       Double_t pT  = trackMC->Pt();
 244       */
 245       //
 246     }
 247
 248   }
 249   //
 250   if (!fESDtrackCuts) {
 251     PostData(1, fListHist);
 252     return;
 253   }
 254   //
 255   // monitor vertex position and selection
 256   //
 257   TH2F * histVertexSelection = (TH2F *) fListHist->FindObject("histVertexSelection");
 258   //
 259   Float_t vertexZ = 0.;
 260   if (IsVertexAccepted(fESD, vertexZ)) {
 261     histVertexSelection->Fill(vertexZ, 0);
 262   } else {
 263     histVertexSelection->Fill(vertexZ, 1);
 264     return;
 265   }
 266   //
 267   // fill track cut variation histograms
 268   //
 269   ProcessTrackCutVariation();
 270   //
 271   // fill ITS->TPC matching histograms
 272   //
 273   ProcessItsTpcMatching();
 274
 275   // Post output data
 276   PostData(1, fListHist);
 277
 278 }
 279
 280
 281 //________________________________________________________________________
 282 void  AliAnalysisTrackingUncertainties::ProcessItsTpcMatching(){
 283   //
 284   // check how many its-sa tracks get matched to TPC
 285   //
 286   int ntr = fESD->GetNumberOfTracks();
 287   //
 288   // initialize histograms
 289   //
 290   TH2F * hNMatch         = (TH2F*) fListHist->FindObject("hNMatch");
 291   TH2F * hBestMatch      = (TH2F*) fListHist->FindObject("hBestMatch");
 292   TH2F * hBestMatch_cuts = (TH2F*) fListHist->FindObject("hBestMatch_cuts");
 293   TH2F * hAllMatch       = (TH2F*) fListHist->FindObject("hAllMatch");
 294   TH2F * hAllMatchGlo    = (TH2F*) fListHist->FindObject("hAllMatchGlo");
 295   TH2F * hPtCorr_ITSTPC  = (TH2F*) fListHist->FindObject("hPtCorr_ITSTPC");
 296   TH2F * hdPtRel_ITSTPC  = (TH2F*) fListHist->FindObject("hdPtRel_ITSTPC");
 297   TH2F * hdInvPtRel_ITSTPC = (TH2F*) fListHist->FindObject("hdInvPtRel_ITSTPC");
 298
 299   //
 300   TH2F * hNMatchBg          = (TH2F*) fListHist->FindObject("hNMatchBg");
 301   TH2F * hBestMatchBg       = (TH2F*) fListHist->FindObject("hBestMatchBg");
 302   TH2F * hBestMatchBg_cuts  = (TH2F*) fListHist->FindObject("hBestMatchBg_cuts");
 303   TH2F * hAllMatchBg        = (TH2F*) fListHist->FindObject("hAllMatchBg");
 304   TH2F * hAllMatchGloBg     = (TH2F*) fListHist->FindObject("hAllMatchGloBg");
 305   TH2F * hdPtRelBg_ITSTPC    = (TH2F*) fListHist->FindObject("hdPtRelBg_ITSTPC");
 306   TH2F * hdInvPtRelBg_ITSTPC = (TH2F*) fListHist->FindObject("hdInvPtRelBg_ITSTPC");
 307
 308   //cout<<" here 1: hdPtRelBg_ITSTPC "<<hdPtRelBg_ITSTPC<<" hdInvPtRelBg_ITSTPC "<<hdInvPtRelBg_ITSTPC<<endl;
 309
 310   //
 311   for (int it=0;it<ntr;it++) {
 312     AliESDtrack* trSA = fESD->GetTrack(it);
 313     if (!trSA->IsOn(AliESDtrack::kITSpureSA) || !trSA->IsOn(AliESDtrack::kITSrefit)) continue;
 314     double pt = trSA->Pt();
 315
 316     // OB - fiducial eta and pt cuts
 317     Double_t etaSA = trSA->Eta();
 318     // std::cout<<" etaSA "<<etaSA<<std::endl; // eta range up to +/- 1.4
 319
 320     Double_t ptSA  = trSA->Pt();
 321
 322     if(TMath::Abs(etaSA)>0.8) continue;
 323
 324     //
 325     Int_t nmatch = 0;
 326     for (int i=kMaxMatch;i--;) {fMatchChi[i]=0; fMatchTr[i]=0;}
 327     for (int it1=0;it1<ntr;it1++){
 328
 329       //std::cout<<" here 0, it1 "<<it1<<" it "<<it<<std::endl;
 330
 331       if (it1==it) continue;
 332
 333       //std::cout<<" here 2, it1 "<<it1<<" it "<<it<<std::endl;
 334
 335       AliESDtrack* trESD = fESD->GetTrack(it1);
 336       if (!trESD->IsOn(AliESDtrack::kTPCrefit)) continue;
 337
 338       //std::cout<<" call match: it1 "<<it1<<" it "<<it<<" nmatch "<<nmatch<<std::endl;
 339       Match(trSA,trESD, nmatch, fExcludeMomFromChi2ITSTPC);
 340       //std::cout<<" left match: it1 "<<it1<<" it "<<it<<" nmatch "<<nmatch<<std::endl;
 341     }
 342     //
 343
 344     // std::cout<<" if "<<it<<" filling nmatch "<<nmatch<<" best chi2 "<<fMatchChi[0]<<std::endl;
 345
 346     hNMatch->Fill(pt,nmatch);
 347     if (nmatch>0){
 348       hBestMatch->Fill(pt,fMatchChi[0]);
 349       hPtCorr_ITSTPC->Fill(pt,fMatchTr[0]->Pt());
 350       hdPtRel_ITSTPC->Fill(pt,(pt-fMatchTr[0]->Pt())/pt);
 351       hdInvPtRel_ITSTPC->Fill(pt,pt*( 1/pt - (1/fMatchTr[0]->Pt()) ));
 352     }
 353
 354     if (nmatch>0 && fESDtrackCuts){
 355
 356       if(fESDtrackCuts->AcceptTrack(fMatchTr[0])){
 357         hBestMatch_cuts->Fill(pt,fMatchChi[0]);
 358       }
 359     }
 360
 361     //
 362     for (int imt=nmatch;imt--;) {
 363       hAllMatch->Fill(pt,fMatchChi[imt]);
 364       if (fMatchTr[imt]->IsOn(AliESDtrack::kITSrefit)) hAllMatchGlo->Fill(pt,fMatchChi[imt]);
 365     }
 366     //
 367     nmatch = 0;
 368     for (int i=kMaxMatch;i--;) {fMatchChi[i]=0; fMatchTr[i]=0;}
 369     for (int it1=0;it1<ntr;it1++) {
 370       if (it1==it) continue;
 371       AliESDtrack* trESD = fESD->GetTrack(it1);
 372       if (!trESD->IsOn(AliESDtrack::kTPCrefit)) continue;
 373       Match(trSA,trESD, nmatch, fExcludeMomFromChi2ITSTPC, TMath::Pi());
 374     }
 375     //
 376     hNMatchBg->Fill(pt,nmatch);
 377     if (nmatch>0){
 378       hBestMatchBg->Fill(pt,fMatchChi[0]);
 379       hdPtRelBg_ITSTPC->Fill(pt,(pt-fMatchTr[0]->Pt())/pt);
 380       hdInvPtRelBg_ITSTPC->Fill(pt,pt*( 1/pt - (1/fMatchTr[0]->Pt()) ));
 381     }
 382
 383     if (nmatch>0 && fESDtrackCuts){
 384       if(fESDtrackCuts->AcceptTrack(fMatchTr[0])){
 385         hBestMatchBg_cuts->Fill(pt,fMatchChi[0]);
 386       }
 387     }
 388
 389     for (int imt=nmatch;imt--;) {
 390       hAllMatchBg->Fill(pt,fMatchChi[imt]);
 391       if (fMatchTr[imt]->IsOn(AliESDtrack::kITSrefit)) hAllMatchGloBg->Fill(pt,fMatchChi[imt]);
 392     }
 393     //
 394   }
 395
 396
 397 }
 398
 399
 400 void AliAnalysisTrackingUncertainties::Match(const AliESDtrack* tr0, const AliESDtrack* tr1, Int_t& nmatch,
 401                                              Bool_t excludeMom, Double_t rotate) {
 402   //
 403   // check if two tracks are matching, possible rotation for combinatoric backgr.
 404   //
 405   Float_t bField = fESD->GetMagneticField();
 406   //
 407   const AliExternalTrackParam* trtpc0 = tr1->GetInnerParam();
 408   if (!trtpc0) return;
 409   AliExternalTrackParam trtpc(*trtpc0);
 410   //
 411   if (TMath::Abs(rotate)>1e-5) {
 412     const double *par = trtpc.GetParameter();
 413     const double *cov = trtpc.GetCovariance();
 414     double alp = trtpc.GetAlpha() + rotate;
 415     trtpc.Set(trtpc.GetX(),alp,par,cov);
 416   }
 417   //
 418   if (!trtpc.Rotate(tr0->GetAlpha())) return;
 419   if (!trtpc.PropagateTo(tr0->GetX(),bField)) return;
 420   double chi2 = tr0->GetPredictedChi2(&trtpc);
 421
 422   //std::cout<<" in Match, nmatch "<<nmatch<<" par[4] before "<<trtpc.GetParameter()[4]<<" chi2 "<<chi2<<endl;
 423
 424   // OB chi2 excluding pt
 425   if(excludeMom){
 426     ((double*)trtpc.GetParameter())[4] = tr0->GetParameter()[4]; // set ITS mom equal TPC mom
 427     chi2 = tr0->GetPredictedChi2(&trtpc);
 428
 429     //std::cout<<" in Match, nmatch "<<nmatch<<" par[4] after "<<trtpc.GetParameter()[4]<<" tr0 mom "<<tr0->GetParameter()[4]
 430     //       <<" chi2 "<<chi2<<std::endl;
 431   }
 432
 433
 434   if (chi2>kMaxChi2) return;
 435
 436   // std::cout<<" found good match, tr1 "<<tr1<<" chi2 "<<chi2<<std::endl;
 437   // std::cout<<" before: fMatchChi[0]  "<<fMatchChi[0]<<" [1] "<<fMatchChi[1]
 438   //       <<" [2]  "<<fMatchChi[2]<<" [3] "<<fMatchChi[3]
 439   //       <<" [4]  "<<fMatchChi[4]<<std::endl;
 440
 441   // std::cout<<" before: fMatchTr[0]  "<<fMatchTr[0]<<" [1] "<<fMatchTr[1]
 442   //       <<" [2]  "<<fMatchTr[2]<<" [3] "<<fMatchTr[3]
 443   //       <<" [4]  "<<fMatchTr[4]<<std::endl;
 444
 445   //
 446   int ins;
 447   for (ins=0;ins<nmatch;ins++) if (chi2<fMatchChi[ins]) break;
 448   if (ins>=kMaxMatch) return;
 449
 450   for (int imv=nmatch;imv>ins;imv--) {
 451     if (imv>=kMaxMatch) continue;
 452     fMatchTr[imv]  = fMatchTr[imv-1];
 453     fMatchChi[imv] = fMatchChi[imv-1];
 454   }
 455   fMatchTr[ins] = tr1;
 456   fMatchChi[ins] = chi2;
 457   nmatch++;
 458   if (nmatch>=kMaxMatch) nmatch = kMaxMatch;
 459   //
 460 }
 461
 462
 463 //________________________________________________________________________
 464 void AliAnalysisTrackingUncertainties::ProcessTrackCutVariation() {
 465   //
 466   // fill track cut variation histograms - undo cuts step-by-step and fill histograms
 467   //
 468   //
 469   // initialize histograms
 470   //
 471   THnF * histNcl         = (THnF *) fListHist->FindObject("histNcl");
 472   THnF * histChi2Tpc     = (THnF *) fListHist->FindObject("histChi2Tpc");
 473   THnF * histDcaZ        = (THnF *) fListHist->FindObject("histDcaZ");
 474   THnF * histSpd         = (THnF *) fListHist->FindObject("histSpd");
 475   THnF * histNcr         = (THnF *) fListHist->FindObject("histNcr");
 476   THnF * histCRoverFC    = (THnF *) fListHist->FindObject("histCRoverFC");
 477   THnF * histChi2Its     = (THnF *) fListHist->FindObject("histChi2Its");
 478   THnF * histTpcLength   = (THnF *) fListHist->FindObject("histTpcLength");
 479   THnF * histTpcItsMatch = (THnF *) fListHist->FindObject("histTpcItsMatch");
 480   //
 481   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z for the vertex cut
 482   //
 483   for (Int_t i=0;i<fESD->GetNumberOfTracks();++i) {
 484     //
 485     AliESDtrack *track =fESD->GetTrack(i);
 486     //
 487     // relevant variables
 488     //
 489     //Double_t pid        = Double_t(GetPid(track));
 490     //
 491     Int_t nclsTPC       = track->GetTPCncls();
 492     Float_t pT          = track->Pt();
 493     Float_t eta         = track->Eta();
 494     Float_t phi         = track->Phi();
 495     Float_t chi2TPC     = track->GetTPCchi2();
 496     Float_t ncrTPC      = track->GetTPCCrossedRows();
 497     Int_t nclsTPCF      = track->GetTPCNclsF();
 498     Float_t nCRoverFC   = track->GetTPCCrossedRows();
 499     Double_t chi2ITS    = track->GetITSchi2();
 500     Int_t nclsITS       = track->GetITSclusters(0);
 501     Float_t tpcLength   = 0.;
 502
 503     if (track->GetInnerParam() && track->GetESDEvent()) {
 504       tpcLength = track->GetLengthInActiveZone(1, 1.8, 220, track->GetESDEvent()->GetMagneticField());
 505     }
 506
 507     if (nclsTPC != 0) {
 508       chi2TPC /= nclsTPC;
 509     } else {
 510       chi2TPC = 999.;
 511     }
 512
 513     if (nclsTPCF !=0) {
 514       nCRoverFC /= nclsTPCF;
 515     } else {
 516       nCRoverFC = 999.;
 517     }
 518
 519     if (nclsITS != 0){
 520       chi2ITS /= nclsITS;
 521     }else {
 522       chi2ITS = 999.;
 523     }
 524     //
 525     track->GetImpactParameters(dca, cov);
 526     //
 527     // (1.) fill number of clusters histogram
 528     //
 529     Int_t minNclsTPC = fESDtrackCuts->GetMinNClusterTPC();
 530     fESDtrackCuts->SetMinNClustersTPC(0);
 531     if (fESDtrackCuts->AcceptTrack(track)) {
 532       for(Int_t iPid = 0; iPid < 6; iPid++) {
 533         Double_t vecHistNcl[kNumberOfAxes] = {static_cast<Double_t>(nclsTPC), pT, eta, phi, static_cast<Double_t>(iPid)};
 534         if (IsConsistentWithPid(iPid, track)) histNcl->Fill(vecHistNcl);
 535       }
 536     }
 537     fESDtrackCuts->SetMinNClustersTPC(minNclsTPC);
 538     //
 539     // (2.) fill chi2 TPC histogram
 540     //
 541     Float_t maxChi2 = fESDtrackCuts->GetMaxChi2PerClusterTPC();
 542     fESDtrackCuts->SetMaxChi2PerClusterTPC(999.);
 543     if (fESDtrackCuts->AcceptTrack(track)) {
 544       for(Int_t iPid = 0; iPid < 6; iPid++) {
 545         Double_t vecHistChi2Tpc[kNumberOfAxes] = {chi2TPC, pT, eta, phi, static_cast<Double_t>(iPid)};
 546         if (IsConsistentWithPid(iPid, track)) histChi2Tpc->Fill(vecHistChi2Tpc);
 547       }
 548     }
 549     fESDtrackCuts->SetMaxChi2PerClusterTPC(maxChi2);
 550     //
 551     // (3.) fill dca_z histogram
 552     //
 553     Float_t maxDcaZ = fESDtrackCuts->GetMaxDCAToVertexZ();
 554     fESDtrackCuts->SetMaxDCAToVertexZ(999.);
 555     if (fESDtrackCuts->AcceptTrack(track)) {
 556       for(Int_t iPid = 0; iPid < 6; iPid++) {
 557         Double_t vecHistDcaZ[kNumberOfAxes] = {TMath::Abs(dca[1]), pT, eta, phi, static_cast<Double_t>(iPid)};
 558         if (IsConsistentWithPid(iPid, track)) histDcaZ->Fill(vecHistDcaZ);
 559       }
 560     }
 561     fESDtrackCuts->SetMaxDCAToVertexZ(maxDcaZ);
 562     //
 563     // (4.) fill hit in SPD histogram
 564     //
 565     fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff);
 566     if (fESDtrackCuts->AcceptTrack(track)) {
 567       Int_t hasPoint = 0;
 568       if (track->HasPointOnITSLayer(0) || track->HasPointOnITSLayer(1)) hasPoint = 1;
 569       for(Int_t iPid = 0; iPid < 6; iPid++) {
 570         Double_t vecHistSpd[kNumberOfAxes] = {static_cast<Double_t>(hasPoint), pT, eta, phi, static_cast<Double_t>(iPid)};
 571         if (IsConsistentWithPid(iPid, track)) histSpd->Fill(vecHistSpd);
 572       }
 573     }
 574     fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
 575     //
 576     // (5.) fill number of crossed rows histogram
 577     //
 578     //Int_t minNcrTPC = fESDtrackCuts->GetMinNCrossedRowsTPC(); //wait for getter in ESDtrackCuts
 579     Int_t minNcrTPC = 0;  //for now use standard value from 2010 !!
 580     fESDtrackCuts->SetMinNCrossedRowsTPC(0);
 581     if (fESDtrackCuts->AcceptTrack(track)) {
 582       for(Int_t iPid = 0; iPid < 6; iPid++) {
 583         Double_t vecHistNcr[kNumberOfAxes] = {static_cast<Double_t>(ncrTPC), pT, eta, phi, static_cast<Double_t>(iPid)};
 584         if (IsConsistentWithPid(iPid, track)) histNcr->Fill(vecHistNcr);
 585       }
 586     }
 587     fESDtrackCuts->SetMinNCrossedRowsTPC(minNcrTPC);
 588     //
 589     // (6.) fill crossed rows over findable clusters histogram
 590     //
 591     //Int_t minCRoverFC = fESDtrackCuts->GetMinRatioCrossedRowsOverFindableClustersTPC(); //wait for getter in ESDtrackCuts
 592     Int_t minCRoverFC = 0.; //for now use standard value from 2010 !!
 593     fESDtrackCuts->SetMinRatioCrossedRowsOverFindableClustersTPC(0.);
 594     if (fESDtrackCuts->AcceptTrack(track)) {
 595        for(Int_t iPid = 0; iPid < 6; iPid++) {
 596          Double_t vecHistCRoverFC[kNumberOfAxes] = {static_cast<Double_t>(nCRoverFC), pT, eta, phi, static_cast<Double_t>(iPid)};
 597          if (IsConsistentWithPid(iPid, track)) histCRoverFC->Fill(vecHistCRoverFC);
 598        }
 599     }
 600     fESDtrackCuts->SetMinRatioCrossedRowsOverFindableClustersTPC(minCRoverFC);
 601     //
 602     // (7.) fill chi2 ITS histogram
 603     //
 604     Float_t maxChi2ITS = fESDtrackCuts->GetMaxChi2PerClusterITS();
 605     fESDtrackCuts->SetMaxChi2PerClusterITS(999.);
 606     if (fESDtrackCuts->AcceptTrack(track)) {
 607       for(Int_t iPid = 0; iPid < 6; iPid++) {
 608         Double_t vecHistChi2ITS[kNumberOfAxes] = {chi2ITS, pT, eta, phi, static_cast<Double_t>(iPid)};
 609         if (IsConsistentWithPid(iPid, track)) histChi2Its->Fill(vecHistChi2ITS);
 610       }
 611     }
 612     fESDtrackCuts->SetMaxChi2PerClusterITS(maxChi2ITS);
 613     //
 614     // (8.) fill active length in TPC histogram
 615     //
 616     Int_t minTpcLength = fESDtrackCuts->GetMinLengthActiveVolumeTPC();
 617     fESDtrackCuts->SetMinLengthActiveVolumeTPC(0);
 618     if (fESDtrackCuts->AcceptTrack(track)) {
 619       for(Int_t iPid = 0; iPid < 6; iPid++) {
 620         Double_t vecHistTpcLength[kNumberOfAxes] = {tpcLength, pT, eta, phi, static_cast<Double_t>(iPid)};
 621         if (IsConsistentWithPid(iPid, track)) histTpcLength->Fill(vecHistTpcLength);
 622       }
 623     }
 624     fESDtrackCuts->SetMinLengthActiveVolumeTPC(minTpcLength);
 625     //
 626     // (9.) fill TPC->ITS matching efficiency histogram
 627     //
 628     Bool_t isMatched = kFALSE;
 629     // remove all ITS requirements
 630     //
 631     // Leonardo and Emilia:
 632     //  -> if MC is available: fill it only for true primaries,
 633     //        --to be done for every cut?
 634     //  -> Postprocessing: plot histogram with 1 divided by histogram with 0 as a function of pT/eta/phi
 635     //  -> Do we want to remove the DCA cut?
 636     Bool_t refit=fESDtrackCuts->GetRequireITSRefit();
 637     Float_t chi2tpc= fESDtrackCuts->GetMaxChi2TPCConstrainedGlobal();
 638     Float_t chi2its= fESDtrackCuts->GetMaxChi2PerClusterITS();
 639     //TString str = fESDtrackCuts->GetMaxDCAToVertexXYPtDep();
 640
 641     fESDtrackCuts->SetRequireITSRefit(kFALSE);
 642     fESDtrackCuts->SetMaxChi2TPCConstrainedGlobal(99999.);
 643     fESDtrackCuts->SetMaxChi2PerClusterITS(999999.);
 644         //TString str = fESDtrackCuts->GetMaxDCAToVertexXYPtDep();
 645     //fESDtrackCuts->SetMaxDCAToVertexXYPtDep();
 646     fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff);
 647
 648     if (fESDtrackCuts->AcceptTrack(track)) {
 649       for(Int_t iPid = 0; iPid < 6; iPid++) {
 650         Double_t vecHistTpcItsMatch[kNumberOfAxes] = {static_cast<Double_t>(isMatched), pT, eta, phi, static_cast<Double_t>(iPid)};
 651         if (IsConsistentWithPid(iPid, track)) histTpcItsMatch->Fill(vecHistTpcItsMatch); // fill with 1 here
 652       }
 653     }
 654     //apply back the cuts
 655     fESDtrackCuts->SetRequireITSRefit(refit);
 656     fESDtrackCuts->SetMaxChi2TPCConstrainedGlobal(chi2tpc);
 657     fESDtrackCuts->SetMaxChi2PerClusterITS(chi2its);
 658     //fESDtrackCuts->SetMaxDCAToVertexXYPtDep(str.Data());
 659     fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
 660     //set is matched
 661     isMatched=kTRUE;
 662     if (fESDtrackCuts->AcceptTrack(track)) {
 663       for(Int_t iPid = 0; iPid < 6; iPid++) {
 664         Double_t vecHistTpcItsMatch[kNumberOfAxes] = {static_cast<Double_t>(isMatched), pT, eta, phi, static_cast<Double_t>(iPid)};
 665         if (IsConsistentWithPid(iPid, track)) histTpcItsMatch->Fill(vecHistTpcItsMatch); // fill with 0 here
 666       }
 667     }
 668
 669   } // end of track loop
 670
 671
 672 }
 673
 674
 675
 676 //________________________________________________________________________
 677 void AliAnalysisTrackingUncertainties::Terminate(Option_t *)
 678 {
 679   // Draw result to the screen
 680   // Called once at the end of the query
 681
 682
 683 }
 684
 685
 686 //________________________________________________________________________
 687 void AliAnalysisTrackingUncertainties::BinLogAxis(const THn *h, Int_t axisNumber) {
 688   //
 689   // Method for the correct logarithmic binning of histograms
 690   //
 691   TAxis *axis = h->GetAxis(axisNumber);
 692   int bins = axis->GetNbins();
 693
 694   Double_t from = axis->GetXmin();
 695   Double_t to = axis->GetXmax();
 696   Double_t *newBins = new Double_t[bins + 1];
 697
 698   newBins[0] = from;
 699   Double_t factor = pow(to/from, 1./bins);
 700
 701   for (int i = 1; i <= bins; i++) {
 702    newBins[i] = factor * newBins[i-1];
 703   }
 704   axis->Set(bins, newBins);
 705   delete [] newBins;
 706
 707 }
 708
 709
 710 //________________________________________________________________________
 711 Bool_t AliAnalysisTrackingUncertainties::IsVertexAccepted(AliESDEvent * esd, Float_t &vertexZ) {
 712   //
 713   // function to check if a proper vertex is reconstructed and write z-position in vertexZ
 714   //
 715   vertexZ = -999.;
 716   Bool_t vertexOkay = kFALSE;
 717   const AliESDVertex *vertex = esd->GetPrimaryVertexTracks();
 718   if (vertex->GetNContributors() < 1) {
 719     //
 720     vertex = esd->GetPrimaryVertexSPD();
 721     if (vertex->GetNContributors() < 1) {
 722       vertexOkay = kFALSE; }
 723     else {
 724       vertexOkay = kTRUE;
 725     }
 726     //
 727     TString vtxTyp = vertex->GetTitle();
 728     Double_t cov[6]={0};
 729     vertex->GetCovarianceMatrix(cov);
 730     Double_t zRes = TMath::Sqrt(cov[5]);
 731     if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) vertexOkay = kFALSE;
 732   }
 733   else {
 734     vertexOkay = kTRUE;
 735   }
 736
 737   vertexZ = vertex->GetZ();
 738   return vertexOkay;
 739
 740 }
 741
 742 //________________________________________________________________________
 743 AliAnalysisTrackingUncertainties::ESpecies_t AliAnalysisTrackingUncertainties::GetPid(const AliESDtrack * const tr, Bool_t useTPCTOF) const {
 744     //
 745     // Determine particle species for a given track
 746     // Two approaches can be used: As default the selection is done using TPC-only, in addition
 747     // the TOF usage is optional. In case of TPC-TOF, a valid TOF signal has to be provided for
 748     // the given track. The identification is delegated to helper function for each species.
 749     // Tracks which are selected as more than one species (ambiguous decision) are rejected.
 750     //
 751     // @Return: Particles species (kUndef in case no identification is possible)
 752     //
 753     if(!fESDpid) return kUndef;
 754     if(useTPCTOF && !(tr->GetStatus() & AliVTrack::kTOFpid)) return kUndef;
 755
 756     Bool_t isElectron(kFALSE), isPion(kFALSE), isKaon(kFALSE), isProton(kFALSE);
 757     Int_t nspec(0);
 758     if((isElectron = IsElectron(tr, useTPCTOF))) nspec++;
 759     if((isPion = IsPion(tr, useTPCTOF))) nspec++;
 760     if((isKaon = IsKaon(tr, useTPCTOF))) nspec++;
 761     if((isProton = IsProton(tr,useTPCTOF))) nspec++;
 762     if(nspec != 1) return kUndef;   // No decision or ambiguous decision;
 763     if(isElectron) return kSpecElectron;
 764     if(isPion) return kSpecPion;
 765     if(isProton) return kSpecProton;
 766     if(isKaon) return kSpecKaon;
 767     return kUndef;
 768 }
 769
 770 //________________________________________________________________________
 771 Bool_t AliAnalysisTrackingUncertainties::IsElectron(const AliESDtrack * const tr, Bool_t useTPCTOF) const {
 772     //
 773     // Selection of electron candidates using the upper half of the TPC sigma band, starting at
 774     // the mean ignoring its shift, and going up to 3 sigma above the mean. In case TOF information
 775     // is available, tracks which are incompatible with electrons within 3 sigma are rejected. If
 776     // no TOF information is used, the momentum regions where the kaon and the proton line cross
 777     // the electron line are cut out using a 3 sigma cut around the kaon or proton line.
 778     //
 779
 780     Float_t nsigmaElectronTPC = fESDpid->NumberOfSigmasTPC(tr, AliPID::kElectron);
 781     if(nsigmaElectronTPC < 0 || nsigmaElectronTPC  > 3) return kFALSE;
 782
 783     if(useTPCTOF){
 784         Float_t nsigmaElectronTOF = fESDpid->NumberOfSigmasTOF(tr, AliPID::kElectron);
 785         if(TMath::Abs(nsigmaElectronTOF) > 3) return kFALSE;
 786         else return kTRUE;
 787     } else {
 788         Float_t nsigmaKaonTPC = fESDpid->NumberOfSigmasTPC(tr, AliPID::kKaon),
 789                 nsigmaProtonTPC =fESDpid->NumberOfSigmasTPC(tr, AliPID::kProton);
 790         if(TMath::Abs(nsigmaKaonTPC < 3) || TMath::Abs(nsigmaProtonTPC < 3)) return kFALSE;
 791         else return kTRUE;
 792     }
 793 }
 794
 795 //________________________________________________________________________
 796 Bool_t AliAnalysisTrackingUncertainties::IsConsistentWithPid(Int_t type, const AliESDtrack * const tr) {
 797   //
 798   // just check if the PID is consistent with a given hypothesis in order to
 799   // investigate effects which are only dependent on the energy loss.
 800   //
 801   if (type == kSpecPion)     return IsPion(tr);
 802   if (type == kSpecKaon)     return IsKaon(tr);
 803   if (type == kSpecProton)   return IsProton(tr);
 804   if (type == kSpecElectron) return IsElectron(tr);
 805   if (type == kAll)          return kTRUE;
 806   return kFALSE;
 807
 808 }
 809
 810 //________________________________________________________________________
 811 Bool_t AliAnalysisTrackingUncertainties::IsPion(const AliESDtrack * const tr, Bool_t /*useTPCPTOF*/) const{
 812   //
 813   // Selectron of pion candidates
 814   // @TODO: To be implemented
 815   //
 816   Float_t nsigmaPionTPC = fESDpid->NumberOfSigmasTPC(tr, AliPID::kPion);
 817   if (TMath::Abs(nsigmaPionTPC) < 3) return kTRUE;
 818   return kFALSE;
 819
 820 }
 821
 822 //________________________________________________________________________
 823 Bool_t AliAnalysisTrackingUncertainties::IsKaon(const AliESDtrack * const tr, Bool_t /*useTPCPTOF*/) const {
 824   //
 825   // Selection of kaon candidates
 826   // @TODO: To be implemented
 827   //
 828   Float_t nsigmaKaonTPC = fESDpid->NumberOfSigmasTPC(tr, AliPID::kKaon);
 829   if (TMath::Abs(nsigmaKaonTPC) < 3) return kTRUE;
 830   return kFALSE;
 831
 832 }
 833
 834 //________________________________________________________________________
 835 Bool_t AliAnalysisTrackingUncertainties::IsProton(const AliESDtrack * const tr, Bool_t /*useTPCPTOF*/) const{
 836   //
 837   // Selection of proton candidates
 838   // @TODO: To be implemented
 839   //
 840   Float_t nsigmaProtonTPC = fESDpid->NumberOfSigmasTPC(tr, AliPID::kProton);
 841   if (TMath::Abs(nsigmaProtonTPC) < 3) return kTRUE;
 842   return kFALSE;
 843
 844 }
 845
 846 //________________________________________________________________________
 847 void AliAnalysisTrackingUncertainties::InitializeTrackCutHistograms() {
 848   //
 849   // create histograms for the track cut studies
 850   //
 851   //
 852   // (1.) number of clusters
 853   //                               0-ncl, 1-pt, 2-eta,         3-phi, 4-pid(0-3 -> electron-proton, 4 -> undef, 5 -> all)
 854   Int_t    binsNcl[kNumberOfAxes] = { 40,   50,    20,            18,  6};
 855   Double_t minNcl[kNumberOfAxes]  = { 20,  0.1,    -1,             0, -0.5};
 856   Double_t maxNcl[kNumberOfAxes]  = {160,   20,    +1, 2*TMath::Pi(),  5.5};
 857   //
 858   TString axisNameNcl[kNumberOfAxes]  = {"ncl","pT","eta","phi","pid"};
 859   TString axisTitleNcl[kNumberOfAxes] = {"ncl","pT","eta","phi","pid"};
 860   //
 861   THnF * histNcl = new THnF("histNcl","number of clusters histogram",kNumberOfAxes, binsNcl, minNcl, maxNcl);
 862   BinLogAxis(histNcl, 1);
 863   fListHist->Add(histNcl);
 864   //
 865   for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
 866     histNcl->GetAxis(iaxis)->SetName(axisNameNcl[iaxis]);
 867     histNcl->GetAxis(iaxis)->SetTitle(axisTitleNcl[iaxis]);
 868   }
 869   //
 870   // (2.) chi2/cls-TPC
 871   //                                  0-chi2, 1-pt, 2-eta,         3-phi, 4-pid(0-3 -> electron-proton, 4 -> undef, 5 -> all)
 872   Int_t    binsChi2Tpc[kNumberOfAxes] = { 40,   50,    20,            18,  6};
 873   Double_t minChi2Tpc[kNumberOfAxes]  = {  0,  0.1,    -1,             0, -0.5};
 874   Double_t maxChi2Tpc[kNumberOfAxes]  = {  8,   20,    +1, 2*TMath::Pi(),  5.5};
 875   //
 876   TString axisNameChi2Tpc[kNumberOfAxes]  = {"chi2tpc","pT","eta","phi","pid"};
 877   TString axisTitleChi2Tpc[kNumberOfAxes] = {"chi2tpc","pT","eta","phi","pid"};
 878   //
 879   THnF * histChi2Tpc = new THnF("histChi2Tpc","chi2 per cls. in TPC",kNumberOfAxes, binsChi2Tpc, minChi2Tpc, maxChi2Tpc);
 880   BinLogAxis(histChi2Tpc, 1);
 881   fListHist->Add(histChi2Tpc);
 882   //
 883   for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
 884     histChi2Tpc->GetAxis(iaxis)->SetName(axisNameChi2Tpc[iaxis]);
 885     histChi2Tpc->GetAxis(iaxis)->SetTitle(axisTitleChi2Tpc[iaxis]);
 886   }
 887   //
 888   // (3.) dca_z
 889   //                               0-dcaZ, 1-pt, 2-eta,         3-phi, 4-pid(0-3 -> electron-proton, 4 -> undef, 5 -> all)
 890   Int_t    binsDcaZ[kNumberOfAxes] = { 20,   50,    20,            18,  6};
 891   Double_t minDcaZ[kNumberOfAxes]  = {  0,  0.1,    -1,             0, -0.5};
 892   Double_t maxDcaZ[kNumberOfAxes]  = {  4,   20,    +1, 2*TMath::Pi(),  5.5};
 893   //
 894   TString axisNameDcaZ[kNumberOfAxes]  = {"dcaZ","pT","eta","phi","pid"};
 895   TString axisTitleDcaZ[kNumberOfAxes] = {"dcaZ","pT","eta","phi","pid"};
 896   //
 897   THnF * histDcaZ = new THnF("histDcaZ","dca_z to prim. vtx.",kNumberOfAxes, binsDcaZ, minDcaZ, maxDcaZ);
 898   BinLogAxis(histDcaZ, 1);
 899   fListHist->Add(histDcaZ);
 900   //
 901   for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
 902     histDcaZ->GetAxis(iaxis)->SetName(axisNameDcaZ[iaxis]);
 903     histDcaZ->GetAxis(iaxis)->SetTitle(axisTitleDcaZ[iaxis]);
 904   }
 905   //
 906   // (4.) hit in SPD layer
 907   //                              0-spdHit, 1-pt, 2-eta,         3-phi, 4-pid(0-3 -> electron-proton, 4 -> undef, 5 -> all)
 908   Int_t    binsSpd[kNumberOfAxes] = {    2,   50,    20,            18,  6};
 909   Double_t minSpd[kNumberOfAxes]  = { -0.5,  0.1,    -1,             0, -0.5};
 910   Double_t maxSpd[kNumberOfAxes]  = {  1.5,   20,    +1, 2*TMath::Pi(),  5.5};
 911   //
 912   TString axisNameSpd[kNumberOfAxes]  = {"spdHit","pT","eta","phi","pid"};
 913   TString axisTitleSpd[kNumberOfAxes] = {"spdHit","pT","eta","phi","pid"};
 914   //
 915   THnF * histSpd = new THnF("histSpd","hit in SPD layer or not",kNumberOfAxes, binsSpd, minSpd, maxSpd);
 916   BinLogAxis(histSpd, 1);
 917   fListHist->Add(histSpd);
 918   //
 919   for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
 920     histSpd->GetAxis(iaxis)->SetName(axisNameSpd[iaxis]);
 921     histSpd->GetAxis(iaxis)->SetTitle(axisTitleSpd[iaxis]);
 922   }
 923   //
 924   // (5.) number of crossed rows
 925   //                               0-ncr, 1-pt, 2-eta,         3-phi, 4-pid(0,unid,etc.)
 926   Int_t    binsNcr[kNumberOfAxes] = { 40,   50,    20,            18,  6};
 927   Double_t minNcr[kNumberOfAxes]  = { 20,  0.1,    -1,             0, -0.5};
 928   Double_t maxNcr[kNumberOfAxes]  = {160,   20,    +1, 2*TMath::Pi(),  5.5};
 929   //
 930   TString axisNameNcr[kNumberOfAxes]  = {"Ncr","pT","eta","phi","pid"};
 931   TString axisTitleNcr[kNumberOfAxes] = {"Ncr","pT","eta","phi","pid"};
 932   //
 933   THnF * histNcr = new THnF("histNcr","number of crossed rows TPC histogram",kNumberOfAxes, binsNcr, minNcr, maxNcr);
 934   BinLogAxis(histNcr, 1);
 935   fListHist->Add(histNcr);
 936   //
 937   for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
 938     histNcr->GetAxis(iaxis)->SetName(axisNameNcr[iaxis]);
 939     histNcr->GetAxis(iaxis)->SetTitle(axisTitleNcr[iaxis]);
 940   }
 941   //
 942   // (6.) ratio crossed rows over findable clusters
 943   //                                0-CRoverFC, 1-pt, 2-eta,         3-phi, 4-pid(0,unid,etc.)
 944   Int_t    binsCRoverFC[kNumberOfAxes] = {  26,   50,    20,            18,  6};
 945   Double_t minCRoverFC[kNumberOfAxes]  = { 0.4,  0.1,    -1,             0, -0.5};
 946   Double_t maxCRoverFC[kNumberOfAxes]  = { 1.8,   20,    +1, 2*TMath::Pi(),  5.5};
 947   //
 948   TString axisNameCRoverFC[kNumberOfAxes]  = {"CRoverFC","pT","eta","phi","pid"};
 949   TString axisTitleCRoverFC[kNumberOfAxes] = {"CRoverFC","pT","eta","phi","pid"};
 950   //
 951   THnF * histCRoverFC = new THnF("histCRoverFC","number of crossed rows over findable clusters histogram",kNumberOfAxes, binsCRoverFC, minCRoverFC, maxCRoverFC);
 952   BinLogAxis(histCRoverFC, 1);
 953   fListHist->Add(histCRoverFC);
 954   //
 955   for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
 956     histCRoverFC->GetAxis(iaxis)->SetName(axisNameCRoverFC[iaxis]);
 957     histCRoverFC->GetAxis(iaxis)->SetTitle(axisTitleCRoverFC[iaxis]);
 958   }
 959   //
 960   // (7.) max chi2 / ITS cluster
 961   //                               0-Chi2Its, 1-pt, 2-eta,         3-phi, 4-pid(0,unid,etc.)
 962   Int_t    binsChi2Its[kNumberOfAxes] = { 25,   50,    20,            18,  6};
 963   Double_t minChi2Its[kNumberOfAxes]  = {  0,  0.1,    -1,             0, -0.5};
 964   Double_t maxChi2Its[kNumberOfAxes]  = { 50,   20,    +1, 2*TMath::Pi(),  5.5};
 965   //
 966   TString axisNameChi2Its[kNumberOfAxes]  = {"Chi2Its","pT","eta","phi","pid"};
 967   TString axisTitleChi2Its[kNumberOfAxes] = {"Chi2Its","pT","eta","phi","pid"};
 968   //
 969   THnF * histChi2Its = new THnF("histChi2Its","number of crossed rows TPC histogram",kNumberOfAxes, binsChi2Its, minChi2Its, maxChi2Its);
 970   BinLogAxis(histChi2Its, 1);
 971   fListHist->Add(histChi2Its);
 972   //
 973   for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
 974     histChi2Its->GetAxis(iaxis)->SetName(axisNameChi2Its[iaxis]);
 975     histChi2Its->GetAxis(iaxis)->SetTitle(axisTitleChi2Its[iaxis]);
 976   }
 977   //
 978   // (8.) tpc active volume length
 979   //                                0-TpcLength, 1-pt, 2-eta,         3-phi, 4-pid(0,unid,etc.)
 980   Int_t    binsTpcLength[kNumberOfAxes] = {  40,   50,    20,            18,  6};
 981   Double_t minTpcLength[kNumberOfAxes]  = {   0,  0.1,    -1,             0, -0.5};
 982   Double_t maxTpcLength[kNumberOfAxes]  = { 170,   20,    +1, 2*TMath::Pi(),  5.5};
 983   //
 984   TString axisNameTpcLength[kNumberOfAxes]  = {"TpcLength","pT","eta","phi","pid"};
 985   TString axisTitleTpcLength[kNumberOfAxes] = {"TpcLength","pT","eta","phi","pid"};
 986   //
 987   THnF * histTpcLength = new THnF("histTpcLength","number of crossed rows TPC histogram",kNumberOfAxes, binsTpcLength, minTpcLength, maxTpcLength);
 988   BinLogAxis(histTpcLength, 1);
 989   fListHist->Add(histTpcLength);
 990   //
 991   for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
 992     histTpcLength->GetAxis(iaxis)->SetName(axisNameTpcLength[iaxis]);
 993     histTpcLength->GetAxis(iaxis)->SetTitle(axisTitleTpcLength[iaxis]);
 994   }
 995   //
 996   // (9.) match TPC->ITS
 997   //                                  0-is matched, 1-pt, 2-eta,         3-phi, 4-pid(0-3 -> electron-proton, 4 -> undef, 5 -> all)
 998   Int_t    binsTpcItsMatch[kNumberOfAxes] = {    2,   50,    20,            18,  6};
 999   Double_t minTpcItsMatch[kNumberOfAxes]  = { -0.5,  0.1,    -1,             0, -0.5};
1000   Double_t maxTpcItsMatch[kNumberOfAxes]  = {  1.5,   20,    +1, 2*TMath::Pi(),  5.5};
1001   //
1002   TString axisNameTpcItsMatch[kNumberOfAxes]  = {"isMatched","pT","eta","phi","pid"};
1003   TString axisTitleTpcItsMatch[kNumberOfAxes] = {"isMatched","pT","eta","phi","pid"};
1004   //
1005   THnF * histTpcItsMatch = new THnF("histTpcItsMatch","TPC -> ITS matching",kNumberOfAxes, binsTpcItsMatch, minTpcItsMatch, maxTpcItsMatch);
1006   BinLogAxis(histTpcItsMatch, 1);
1007   fListHist->Add(histTpcItsMatch);
1008   //
1009   for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
1010     histTpcItsMatch->GetAxis(iaxis)->SetName(axisNameTpcItsMatch[iaxis]);
1011     histTpcItsMatch->GetAxis(iaxis)->SetTitle(axisTitleTpcItsMatch[iaxis]);
1012   }
1013
1014
1015
1016
1017 }
1018