]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGPP/EvTrkSelection/AliAnalysisTrackingUncertainties.cxx
Dump object functionality
[u/mrichter/AliRoot.git] / PWGPP / EvTrkSelection / AliAnalysisTrackingUncertainties.cxx
index eaba463bf1294bac50c0cd11211dbd9b85f4b296..be651d5d35fd18dc9878d272683643c6dffd0202 100644 (file)
@@ -64,7 +64,8 @@ AliAnalysisTrackingUncertainties::AliAnalysisTrackingUncertainties()
     fListHist(0),
     fESDtrackCuts(0),
     fMatchTr(),
-    fMatchChi()
+    fMatchChi(),
+    fExcludeMomFromChi2ITSTPC(0) 
 {
   // default Constructor
   /* fast compilation test
@@ -86,7 +87,8 @@ AliAnalysisTrackingUncertainties::AliAnalysisTrackingUncertainties(const char *n
     fListHist(0),
     fESDtrackCuts(0),
     fMatchTr(),
-    fMatchChi()
+    fMatchChi(),
+    fExcludeMomFromChi2ITSTPC(0) 
 {
   //
   // standard constructur which should be used
@@ -120,37 +122,75 @@ void AliAnalysisTrackingUncertainties::UserCreateOutputObjects()
   fListHist = new TList();
   fListHist->SetOwner(kTRUE);
   //
-  // basic QA and statistics histograms
+  // (1.) basic QA and statistics histograms
   //
   TH2F * histVertexSelection = new TH2F("histVertexSelection", "vertex selection; vertex z (cm); accepted/rejected", 100, -50., 50., 2, -0.5, 1.5);
   fListHist->Add(histVertexSelection);
   //
-  // track cut variation histograms
+  // (2.) track cut variation histograms
+  //
   InitializeTrackCutHistograms();
   //
+  // (3.) ITS -> TPC matching histograms
+  //
+  //  Int_t    binsMatch[kNumberOfAxes] = { 10,   50,    20,            18,  6};
+  //  Double_t minMatch[kNumberOfAxes]  = {  0,  0.1,    -1,             0, -0.5};
+  //  Double_t maxMatch[kNumberOfAxes]  = {200,   20,    +1, 2*TMath::Pi(),  5.5};
+  //
+  //  TString axisNameMatch[kNumberOfAxes]  = {"matchChi2","pT","eta","phi","pid"};
+  //  TString axisTitleMatch[kNumberOfAxes] = {"matchChi2","pT","eta","phi","pid"};
+  //
+  //  THnF * hBestMatch = new THnF("hBestMatch","ITS -> TPC matching ",kNumberOfAxes, binsMatch, minMatch, maxMatch);
+  //  for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
+  //    hBestMatch->GetAxis(iaxis)->SetName(axisNameMatch[iaxis]);
+  // hBestMatch->GetAxis(iaxis)->SetTitle(axisTitleMatch[iaxis]);
+  //  }
+  //  BinLogAxis(hBestMatch, 1);
+  //  fListHist->Add(hBestMatch);
+  //
   //
-  // matching histograms
   //
   const int nbPt=40;
   const double ptMax=5;
   //
   TH2F * hNMatch    = new TH2F("hNMatch","N Matches",nbPt,0,ptMax,kMaxMatch+1,-0.5,kMaxMatch+0.5);
-  TH2F * hBestMatch = new TH2F("hBestMatch","Best Match Chi2",nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2);
+  TH2F * hBestMatch = new TH2F("hBestMatch","Best Match Chi2",nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2); // OB
+  TH2F * hBestMatch_cuts = new TH2F("hBestMatch_cuts","Best Match Chi2",nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2);
   TH2F * hAllMatch  = new TH2F("hAllMatch","All Matches Chi2",nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2);
   TH2F * hAllMatchGlo  = new TH2F("hAllMatchGlo","All Matches Chi2",nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2);
+  TH2F * hPtCorr_ITSTPC = new TH2F("hPtCorr_ITSTPC","PtCorr",nbPt,0,ptMax,nbPt,0,ptMax);
+  TH2F * hdPtRel_ITSTPC = new TH2F("hdPtRel_ITSTPC","dPt/pt",nbPt,0,ptMax,2*nbPt+1,-0.4*ptMax,0.4*ptMax);
+  TH2F * hdInvPtRel_ITSTPC = new TH2F("hdInvPtRel_ITSTPC","pt*dPt^{-1}",nbPt,0,ptMax,2*nbPt+1,-0.4*ptMax,0.4*ptMax);
+  //
   fListHist->Add(hNMatch);
   fListHist->Add(hBestMatch);
+  fListHist->Add(hBestMatch_cuts);
   fListHist->Add(hAllMatch);
   fListHist->Add(hAllMatchGlo);
+  fListHist->Add(hPtCorr_ITSTPC);
+  fListHist->Add(hdPtRel_ITSTPC);
+  fListHist->Add(hdInvPtRel_ITSTPC);
+  //
   //
   TH2F * hNMatchBg    = new TH2F("hNMatchBg","N Matches",nbPt,0,ptMax,kMaxMatch+1,-0.5,kMaxMatch+0.5);
   TH2F * hBestMatchBg = new TH2F("hBestMatchBg","Best Match Chi2",nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2);
+  TH2F * hBestMatchBg_cuts = new TH2F("hBestMatchBg_cuts","Best Match Chi2",nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2); // OB
   TH2F * hAllMatchBg  = new TH2F("hAllMatchBg","All Matches Chi2",nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2);
   TH2F * hAllMatchGloBg  = new TH2F("hAllMatchGloBg","All Matches Chi2",nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2);
+  TH2F * hdPtRelBg_ITSTPC = new TH2F("hdPtRelBg_ITSTPC","dPt/pt",nbPt,0,ptMax,2*nbPt+1,-0.4*ptMax,0.4*ptMax);
+  TH2F * hdInvPtRelBg_ITSTPC = new TH2F("hdInvPtRelBg_ITSTPC","pt*dPt^{-1}",nbPt,0,ptMax,2*nbPt+1,-0.4*ptMax,0.4*ptMax);
+
+  //cout<<" here 0 : hdPtRelBg_ITSTPC "<<hdPtRelBg_ITSTPC<<" hdInvPtRelBg_ITSTPC "<<hdInvPtRelBg_ITSTPC<<endl;
+
   fListHist->Add(hNMatchBg);
   fListHist->Add(hBestMatchBg);
+  fListHist->Add(hBestMatchBg_cuts);
   fListHist->Add(hAllMatchBg);
   fListHist->Add(hAllMatchGloBg);
+   fListHist->Add(hdPtRelBg_ITSTPC);
+  fListHist->Add(hdInvPtRelBg_ITSTPC);
+  //add default track cuts in the output list
+  fListHist->Add(fESDtrackCuts);
   //
   // post data
   //
@@ -247,32 +287,78 @@ void  AliAnalysisTrackingUncertainties::ProcessItsTpcMatching(){
   //
   // initialize histograms
   //
-  TH2F * hNMatch      = (TH2F*) fListHist->FindObject("hNMatch");
-  TH2F * hBestMatch   = (TH2F*) fListHist->FindObject("hBestMatch");
-  TH2F * hAllMatch    = (TH2F*) fListHist->FindObject("hAllMatch");
-  TH2F * hAllMatchGlo = (TH2F*) fListHist->FindObject("hAllMatchGlo");  
+  TH2F * hNMatch         = (TH2F*) fListHist->FindObject("hNMatch");
+  TH2F * hBestMatch      = (TH2F*) fListHist->FindObject("hBestMatch");
+  TH2F * hBestMatch_cuts = (TH2F*) fListHist->FindObject("hBestMatch_cuts");
+  TH2F * hAllMatch       = (TH2F*) fListHist->FindObject("hAllMatch");
+  TH2F * hAllMatchGlo    = (TH2F*) fListHist->FindObject("hAllMatchGlo");  
+  TH2F * hPtCorr_ITSTPC  = (TH2F*) fListHist->FindObject("hPtCorr_ITSTPC");
+  TH2F * hdPtRel_ITSTPC  = (TH2F*) fListHist->FindObject("hdPtRel_ITSTPC");
+  TH2F * hdInvPtRel_ITSTPC = (TH2F*) fListHist->FindObject("hdInvPtRel_ITSTPC");
+
   //
-  TH2F * hNMatchBg      = (TH2F*) fListHist->FindObject("hNMatchBg");
-  TH2F * hBestMatchBg   = (TH2F*) fListHist->FindObject("hBestMatchBg");
-  TH2F * hAllMatchBg    = (TH2F*) fListHist->FindObject("hAllMatchBg");
-  TH2F * hAllMatchGloBg = (TH2F*) fListHist->FindObject("hAllMatchGloBg");    
+  TH2F * hNMatchBg          = (TH2F*) fListHist->FindObject("hNMatchBg");
+  TH2F * hBestMatchBg       = (TH2F*) fListHist->FindObject("hBestMatchBg");
+  TH2F * hBestMatchBg_cuts  = (TH2F*) fListHist->FindObject("hBestMatchBg_cuts");
+  TH2F * hAllMatchBg        = (TH2F*) fListHist->FindObject("hAllMatchBg");
+  TH2F * hAllMatchGloBg     = (TH2F*) fListHist->FindObject("hAllMatchGloBg");    
+  TH2F * hdPtRelBg_ITSTPC    = (TH2F*) fListHist->FindObject("hdPtRelBg_ITSTPC");
+  TH2F * hdInvPtRelBg_ITSTPC = (TH2F*) fListHist->FindObject("hdInvPtRelBg_ITSTPC");
+
+  //cout<<" here 1: hdPtRelBg_ITSTPC "<<hdPtRelBg_ITSTPC<<" hdInvPtRelBg_ITSTPC "<<hdInvPtRelBg_ITSTPC<<endl;
+
   //
   for (int it=0;it<ntr;it++) {
     AliESDtrack* trSA = fESD->GetTrack(it);
     if (!trSA->IsOn(AliESDtrack::kITSpureSA) || !trSA->IsOn(AliESDtrack::kITSrefit)) continue;
     double pt = trSA->Pt();
+
+    // OB - fiducial eta and pt cuts
+    Double_t etaSA = trSA->Eta();
+    // std::cout<<" etaSA "<<etaSA<<std::endl; // eta range up to +/- 1.4
+
+    Double_t ptSA  = trSA->Pt();
+
+    if(TMath::Abs(etaSA)>0.8) continue;
+
     //
     Int_t nmatch = 0;
-    for (int i=kMaxMatch;i--;) {fMatchChi[i]=0; fMatchTr[i]=0;} // reset array
-    for (int it1=0;it1<ntr;it1++) { 
+    for (int i=kMaxMatch;i--;) {fMatchChi[i]=0; fMatchTr[i]=0;}
+    for (int it1=0;it1<ntr;it1++){
+
+      //std::cout<<" here 0, it1 "<<it1<<" it "<<it<<std::endl;
+
       if (it1==it) continue;
+      
+      //std::cout<<" here 2, it1 "<<it1<<" it "<<it<<std::endl;
+
       AliESDtrack* trESD = fESD->GetTrack(it1);
       if (!trESD->IsOn(AliESDtrack::kTPCrefit)) continue;
-      Match(trSA,trESD, nmatch);
+      
+      //std::cout<<" call match: it1 "<<it1<<" it "<<it<<" nmatch "<<nmatch<<std::endl;
+      Match(trSA,trESD, nmatch, fExcludeMomFromChi2ITSTPC);
+      //std::cout<<" left match: it1 "<<it1<<" it "<<it<<" nmatch "<<nmatch<<std::endl;
     }
     //
+    
+    // std::cout<<" if "<<it<<" filling nmatch "<<nmatch<<" best chi2 "<<fMatchChi[0]<<std::endl;
+    
     hNMatch->Fill(pt,nmatch);
-    if (nmatch>0) hBestMatch->Fill(pt,fMatchChi[0]);
+    if (nmatch>0){
+      hBestMatch->Fill(pt,fMatchChi[0]);
+      hPtCorr_ITSTPC->Fill(pt,fMatchTr[0]->Pt()); 
+      hdPtRel_ITSTPC->Fill(pt,(pt-fMatchTr[0]->Pt())/pt); 
+      hdInvPtRel_ITSTPC->Fill(pt,pt*( 1/pt - (1/fMatchTr[0]->Pt()) )); 
+    }
+    
+    if (nmatch>0 && fESDtrackCuts){
+      
+      if(fESDtrackCuts->AcceptTrack(fMatchTr[0])){
+       hBestMatch_cuts->Fill(pt,fMatchChi[0]);
+      }
+    }
+    
+    //
     for (int imt=nmatch;imt--;) {
       hAllMatch->Fill(pt,fMatchChi[imt]);
       if (fMatchTr[imt]->IsOn(AliESDtrack::kITSrefit)) hAllMatchGlo->Fill(pt,fMatchChi[imt]);
@@ -284,11 +370,22 @@ void  AliAnalysisTrackingUncertainties::ProcessItsTpcMatching(){
       if (it1==it) continue;
       AliESDtrack* trESD = fESD->GetTrack(it1);
       if (!trESD->IsOn(AliESDtrack::kTPCrefit)) continue;
-      Match(trSA,trESD, nmatch, TMath::Pi());
+      Match(trSA,trESD, nmatch, fExcludeMomFromChi2ITSTPC, TMath::Pi());
     }
     //
     hNMatchBg->Fill(pt,nmatch);
-    if (nmatch>0) hBestMatchBg->Fill(pt,fMatchChi[0]);
+    if (nmatch>0){
+      hBestMatchBg->Fill(pt,fMatchChi[0]);
+      hdPtRelBg_ITSTPC->Fill(pt,(pt-fMatchTr[0]->Pt())/pt); 
+      hdInvPtRelBg_ITSTPC->Fill(pt,pt*( 1/pt - (1/fMatchTr[0]->Pt()) )); 
+    }
+
+    if (nmatch>0 && fESDtrackCuts){
+      if(fESDtrackCuts->AcceptTrack(fMatchTr[0])){
+       hBestMatchBg_cuts->Fill(pt,fMatchChi[0]);
+      }
+    }
+
     for (int imt=nmatch;imt--;) {
       hAllMatchBg->Fill(pt,fMatchChi[imt]);
       if (fMatchTr[imt]->IsOn(AliESDtrack::kITSrefit)) hAllMatchGloBg->Fill(pt,fMatchChi[imt]);
@@ -300,7 +397,8 @@ void  AliAnalysisTrackingUncertainties::ProcessItsTpcMatching(){
 }
 
 
-void AliAnalysisTrackingUncertainties::Match(const AliESDtrack* tr0, const AliESDtrack* tr1, Int_t &nmatch, Double_t rotate) {
+void AliAnalysisTrackingUncertainties::Match(const AliESDtrack* tr0, const AliESDtrack* tr1, Int_t& nmatch, 
+                                            Bool_t excludeMom, Double_t rotate) {
   //
   // check if two tracks are matching, possible rotation for combinatoric backgr.
   // 
@@ -320,7 +418,30 @@ void AliAnalysisTrackingUncertainties::Match(const AliESDtrack* tr0, const AliES
   if (!trtpc.Rotate(tr0->GetAlpha())) return;
   if (!trtpc.PropagateTo(tr0->GetX(),bField)) return;
   double chi2 = tr0->GetPredictedChi2(&trtpc);
+
+  //std::cout<<" in Match, nmatch "<<nmatch<<" par[4] before "<<trtpc.GetParameter()[4]<<" chi2 "<<chi2<<endl;
+
+  // OB chi2 excluding pt 
+  if(excludeMom){
+    ((double*)trtpc.GetParameter())[4] = tr0->GetParameter()[4]; // set ITS mom equal TPC mom
+    chi2 = tr0->GetPredictedChi2(&trtpc);
+
+    //std::cout<<" in Match, nmatch "<<nmatch<<" par[4] after "<<trtpc.GetParameter()[4]<<" tr0 mom "<<tr0->GetParameter()[4]
+    //      <<" chi2 "<<chi2<<std::endl;
+  }
+
+
   if (chi2>kMaxChi2) return;
+
+  // std::cout<<" found good match, tr1 "<<tr1<<" chi2 "<<chi2<<std::endl;
+  // std::cout<<" before: fMatchChi[0]  "<<fMatchChi[0]<<" [1] "<<fMatchChi[1]
+  //      <<" [2]  "<<fMatchChi[2]<<" [3] "<<fMatchChi[3]
+  //      <<" [4]  "<<fMatchChi[4]<<std::endl; 
+
+  // std::cout<<" before: fMatchTr[0]  "<<fMatchTr[0]<<" [1] "<<fMatchTr[1]
+  //      <<" [2]  "<<fMatchTr[2]<<" [3] "<<fMatchTr[3]
+  //      <<" [4]  "<<fMatchTr[4]<<std::endl; 
+
   //
   int ins;
   for (ins=0;ins<nmatch;ins++) if (chi2<fMatchChi[ins]) break;
@@ -347,10 +468,15 @@ void AliAnalysisTrackingUncertainties::ProcessTrackCutVariation() {
   //
   // initialize histograms
   //
-  THnF * histNcl     = (THnF *) fListHist->FindObject("histNcl");
-  THnF * histChi2Tpc = (THnF *) fListHist->FindObject("histChi2Tpc");
-  THnF * histDcaZ    = (THnF *) fListHist->FindObject("histDcaZ");
-  THnF * histSpd     = (THnF *) fListHist->FindObject("histSpd");
+  THnF * histNcl         = (THnF *) fListHist->FindObject("histNcl");
+  THnF * histChi2Tpc     = (THnF *) fListHist->FindObject("histChi2Tpc");
+  THnF * histDcaZ        = (THnF *) fListHist->FindObject("histDcaZ");
+  THnF * histSpd         = (THnF *) fListHist->FindObject("histSpd");
+  THnF * histNcr         = (THnF *) fListHist->FindObject("histNcr");
+  THnF * histCRoverFC    = (THnF *) fListHist->FindObject("histCRoverFC");
+  THnF * histChi2Its     = (THnF *) fListHist->FindObject("histChi2Its");
+  THnF * histTpcLength   = (THnF *) fListHist->FindObject("histTpcLength");
+  THnF * histTpcItsMatch = (THnF *) fListHist->FindObject("histTpcItsMatch");
   //
   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z for the vertex cut
   //
@@ -360,17 +486,41 @@ void AliAnalysisTrackingUncertainties::ProcessTrackCutVariation() {
     //
     // relevant variables
     //
-    Int_t nclsTPC = track->GetTPCncls();
-    Float_t pT  = track->Pt();
-    Float_t eta = track->Eta();
-    Float_t phi = track->Phi();
-    Float_t chi2TPC = track->GetTPCchi2();
-    Double_t pid = Double_t(GetPid(track));
+    //Double_t pid        = Double_t(GetPid(track));
+    //
+    Int_t nclsTPC       = track->GetTPCncls();
+    Float_t pT          = track->Pt();
+    Float_t eta         = track->Eta();
+    Float_t phi         = track->Phi();
+    Float_t chi2TPC     = track->GetTPCchi2();
+    Float_t ncrTPC      = track->GetTPCCrossedRows();
+    Int_t nclsTPCF      = track->GetTPCNclsF(); 
+    Float_t nCRoverFC   = track->GetTPCCrossedRows();
+    Double_t chi2ITS    = track->GetITSchi2();
+    Int_t nclsITS       = track->GetITSclusters(0);
+    Float_t tpcLength   = 0.;
+
+    if (track->GetInnerParam() && track->GetESDEvent()) {
+      tpcLength = track->GetLengthInActiveZone(1, 1.8, 220, track->GetESDEvent()->GetMagneticField());
+    }
+
     if (nclsTPC != 0) {
       chi2TPC /= nclsTPC; 
     } else {
       chi2TPC = 999.;
     }
+
+    if (nclsTPCF !=0) {
+      nCRoverFC /= nclsTPCF;
+    } else {
+      nCRoverFC = 999.;
+    }
+
+    if (nclsITS != 0){
+      chi2ITS /= nclsITS;
+    }else {
+      chi2ITS = 999.;
+    }
     //
     track->GetImpactParameters(dca, cov);
     //
@@ -379,8 +529,10 @@ void AliAnalysisTrackingUncertainties::ProcessTrackCutVariation() {
     Int_t minNclsTPC = fESDtrackCuts->GetMinNClusterTPC();
     fESDtrackCuts->SetMinNClustersTPC(0);
     if (fESDtrackCuts->AcceptTrack(track)) {
-      Double_t vecHistNcl[kNumberOfAxes] = {nclsTPC, pT, eta, phi, pid};
-      histNcl->Fill(vecHistNcl);
+      for(Int_t iPid = 0; iPid < 6; iPid++) {
+       Double_t vecHistNcl[kNumberOfAxes] = {nclsTPC, pT, eta, phi, iPid};
+       if (IsConsistentWithPid(iPid, track)) histNcl->Fill(vecHistNcl);
+      }
     }
     fESDtrackCuts->SetMinNClustersTPC(minNclsTPC);
     //
@@ -389,8 +541,10 @@ void AliAnalysisTrackingUncertainties::ProcessTrackCutVariation() {
     Float_t maxChi2 = fESDtrackCuts->GetMaxChi2PerClusterTPC();
     fESDtrackCuts->SetMaxChi2PerClusterTPC(999.);
     if (fESDtrackCuts->AcceptTrack(track)) {
-      Double_t vecHistChi2Tpc[kNumberOfAxes] = {chi2TPC, pT, eta, phi, pid};
-      histChi2Tpc->Fill(vecHistChi2Tpc);
+      for(Int_t iPid = 0; iPid < 6; iPid++) {
+       Double_t vecHistChi2Tpc[kNumberOfAxes] = {chi2TPC, pT, eta, phi, iPid};
+       if (IsConsistentWithPid(iPid, track)) histChi2Tpc->Fill(vecHistChi2Tpc);
+      }
     }
     fESDtrackCuts->SetMaxChi2PerClusterTPC(maxChi2);
     //
@@ -399,8 +553,10 @@ void AliAnalysisTrackingUncertainties::ProcessTrackCutVariation() {
     Float_t maxDcaZ = fESDtrackCuts->GetMaxDCAToVertexZ();
     fESDtrackCuts->SetMaxDCAToVertexZ(999.);
     if (fESDtrackCuts->AcceptTrack(track)) {
-      Double_t vecHistDcaZ[kNumberOfAxes] = {TMath::Abs(dca[1]), pT, eta, phi, pid};
-      histDcaZ->Fill(vecHistDcaZ);
+      for(Int_t iPid = 0; iPid < 6; iPid++) {
+       Double_t vecHistDcaZ[kNumberOfAxes] = {TMath::Abs(dca[1]), pT, eta, phi, iPid};
+       if (IsConsistentWithPid(iPid, track)) histDcaZ->Fill(vecHistDcaZ);
+      }
     }
     fESDtrackCuts->SetMaxDCAToVertexZ(maxDcaZ);
     //
@@ -410,12 +566,106 @@ void AliAnalysisTrackingUncertainties::ProcessTrackCutVariation() {
     if (fESDtrackCuts->AcceptTrack(track)) {
       Int_t hasPoint = 0;
       if (track->HasPointOnITSLayer(0) || track->HasPointOnITSLayer(1)) hasPoint = 1;
-      Double_t vecHistSpd[kNumberOfAxes] = {hasPoint, pT, eta, phi, pid};
-      histSpd->Fill(vecHistSpd);
+      for(Int_t iPid = 0; iPid < 6; iPid++) {
+       Double_t vecHistSpd[kNumberOfAxes] = {hasPoint, pT, eta, phi, iPid};
+       if (IsConsistentWithPid(iPid, track)) histSpd->Fill(vecHistSpd);
+      }
     }
     fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
-
-
+    //
+    // (5.) fill number of crossed rows histogram
+    //
+    //Int_t minNcrTPC = fESDtrackCuts->GetMinNCrossedRowsTPC(); //wait for getter in ESDtrackCuts
+    Int_t minNcrTPC = 0;  //for now use standard value from 2010 !!
+    fESDtrackCuts->SetMinNCrossedRowsTPC(0);
+    if (fESDtrackCuts->AcceptTrack(track)) {
+      for(Int_t iPid = 0; iPid < 6; iPid++) {
+       Double_t vecHistNcr[kNumberOfAxes] = {ncrTPC, pT, eta, phi, iPid};
+       if (IsConsistentWithPid(iPid, track)) histNcr->Fill(vecHistNcr);
+      }
+    }
+    fESDtrackCuts->SetMinNCrossedRowsTPC(minNcrTPC);
+    //
+    // (6.) fill crossed rows over findable clusters histogram
+    //
+    //Int_t minCRoverFC = fESDtrackCuts->GetMinRatioCrossedRowsOverFindableClustersTPC(); //wait for getter in ESDtrackCuts
+    Int_t minCRoverFC = 0.; //for now use standard value from 2010 !!
+    fESDtrackCuts->SetMinRatioCrossedRowsOverFindableClustersTPC(0.);
+    if (fESDtrackCuts->AcceptTrack(track)) {
+       for(Int_t iPid = 0; iPid < 6; iPid++) {
+        Double_t vecHistCRoverFC[kNumberOfAxes] = {nCRoverFC, pT, eta, phi, iPid};
+        if (IsConsistentWithPid(iPid, track)) histCRoverFC->Fill(vecHistCRoverFC);
+       }
+    }
+    fESDtrackCuts->SetMinRatioCrossedRowsOverFindableClustersTPC(minCRoverFC);
+    //
+    // (7.) fill chi2 ITS histogram
+    //
+    Float_t maxChi2ITS = fESDtrackCuts->GetMaxChi2PerClusterITS();
+    fESDtrackCuts->SetMaxChi2PerClusterITS(999.);
+    if (fESDtrackCuts->AcceptTrack(track)) {
+      for(Int_t iPid = 0; iPid < 6; iPid++) {
+       Double_t vecHistChi2ITS[kNumberOfAxes] = {chi2ITS, pT, eta, phi, iPid};
+       if (IsConsistentWithPid(iPid, track)) histChi2Its->Fill(vecHistChi2ITS);
+      }
+    }
+    fESDtrackCuts->SetMaxChi2PerClusterITS(maxChi2ITS);
+    //
+    // (8.) fill active length in TPC histogram
+    //
+    Int_t minTpcLength = fESDtrackCuts->GetMinLengthActiveVolumeTPC();
+    fESDtrackCuts->SetMinLengthActiveVolumeTPC(0);
+    if (fESDtrackCuts->AcceptTrack(track)) {
+      for(Int_t iPid = 0; iPid < 6; iPid++) {
+       Double_t vecHistTpcLength[kNumberOfAxes] = {tpcLength, pT, eta, phi, iPid};
+       if (IsConsistentWithPid(iPid, track)) histTpcLength->Fill(vecHistTpcLength);
+      }
+    }
+    fESDtrackCuts->SetMinLengthActiveVolumeTPC(minTpcLength);
+    //
+    // (9.) fill TPC->ITS matching efficiency histogram
+    //
+    Bool_t isMatched = kFALSE;
+    // remove all ITS requirements
+    //
+    // Leonardo and Emilia: 
+    //  -> if MC is available: fill it only for true primaries, 
+    //        --to be done for every cut?
+    //  -> Postprocessing: plot histogram with 1 divided by histogram with 0 as a function of pT/eta/phi
+    //  -> Do we want to remove the DCA cut?
+    Bool_t refit=fESDtrackCuts->GetRequireITSRefit();
+    Float_t chi2tpc= fESDtrackCuts->GetMaxChi2TPCConstrainedGlobal();
+    Float_t chi2its= fESDtrackCuts->GetMaxChi2PerClusterITS();
+    //TString str = fESDtrackCuts->GetMaxDCAToVertexXYPtDep();
+    
+    fESDtrackCuts->SetRequireITSRefit(kFALSE);
+    fESDtrackCuts->SetMaxChi2TPCConstrainedGlobal(99999.);
+    fESDtrackCuts->SetMaxChi2PerClusterITS(999999.);
+       //TString str = fESDtrackCuts->GetMaxDCAToVertexXYPtDep();
+    //fESDtrackCuts->SetMaxDCAToVertexXYPtDep();
+    fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff);
+    
+    if (fESDtrackCuts->AcceptTrack(track)) {
+      for(Int_t iPid = 0; iPid < 6; iPid++) {
+       Double_t vecHistTpcItsMatch[kNumberOfAxes] = {isMatched, pT, eta, phi, iPid};
+       if (IsConsistentWithPid(iPid, track)) histTpcItsMatch->Fill(vecHistTpcItsMatch); // fill with 1 here
+      }
+    }
+    //apply back the cuts
+    fESDtrackCuts->SetRequireITSRefit(refit);
+    fESDtrackCuts->SetMaxChi2TPCConstrainedGlobal(chi2tpc);
+    fESDtrackCuts->SetMaxChi2PerClusterITS(chi2its);
+    //fESDtrackCuts->SetMaxDCAToVertexXYPtDep(str.Data());
+    fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
+    //set is matched
+    isMatched=kTRUE;
+    if (fESDtrackCuts->AcceptTrack(track)) {
+      for(Int_t iPid = 0; iPid < 6; iPid++) {
+       Double_t vecHistTpcItsMatch[kNumberOfAxes] = {isMatched, pT, eta, phi, iPid};
+       if (IsConsistentWithPid(iPid, track)) histTpcItsMatch->Fill(vecHistTpcItsMatch); // fill with 0 here
+      }
+    }
+    
   } // end of track loop
 
 
@@ -511,9 +761,9 @@ AliAnalysisTrackingUncertainties::ESpecies_t AliAnalysisTrackingUncertainties::G
     if((isProton = IsProton(tr,useTPCTOF))) nspec++;
     if(nspec != 1) return kUndef;   // No decision or ambiguous decision;
     if(isElectron) return kSpecElectron;
-    if(isPion) return kSpecElectron;
-    if(isProton) return kSpecElectron;
-    if(isKaon) return kSpecElectron;
+    if(isPion) return kSpecPion;
+    if(isProton) return kSpecProton;
+    if(isKaon) return kSpecKaon;
     return kUndef;
 }
 
@@ -543,30 +793,54 @@ Bool_t AliAnalysisTrackingUncertainties::IsElectron(const AliESDtrack * const tr
 }
 
 //________________________________________________________________________
-Bool_t AliAnalysisTrackingUncertainties::IsPion(const AliESDtrack * const /*tr*/, Bool_t /*useTPCPTOF*/) const{
-    //
-    // Selectron of pion candidates
-    // @TODO: To be implemented
-    //
-    return kFALSE;
+Bool_t AliAnalysisTrackingUncertainties::IsConsistentWithPid(Int_t type, const AliESDtrack * const tr) {
+  //
+  // just check if the PID is consistent with a given hypothesis in order to 
+  // investigate effects which are only dependent on the energy loss.
+  //
+  if (type == kSpecPion)     return IsPion(tr);
+  if (type == kSpecKaon)     return IsKaon(tr);
+  if (type == kSpecProton)   return IsProton(tr);
+  if (type == kSpecElectron) return IsElectron(tr);
+  if (type == kAll)          return kTRUE;
+  return kFALSE;
+
 }
 
 //________________________________________________________________________
-Bool_t AliAnalysisTrackingUncertainties::IsKaon(const AliESDtrack * const /*tr*/, Bool_t /*useTPCPTOF*/) const {
-    //
-    // Selection of kaon candidates
-    // @TODO: To be implemented
-    //
-    return kFALSE;
+Bool_t AliAnalysisTrackingUncertainties::IsPion(const AliESDtrack * const tr, Bool_t /*useTPCPTOF*/) const{
+  //
+  // Selectron of pion candidates
+  // @TODO: To be implemented
+  //
+  Float_t nsigmaPionTPC = fESDpid->NumberOfSigmasTPC(tr, AliPID::kPion);
+  if (TMath::Abs(nsigmaPionTPC) < 3) return kTRUE;
+  return kFALSE;
+
 }
 
 //________________________________________________________________________
-Bool_t AliAnalysisTrackingUncertainties::IsProton(const AliESDtrack * const /*tr*/, Bool_t /*useTPCPTOF*/) const{
-    // 
-    // Selection of proton candidates
-    // @TODO: To be implemented
-    //
-    return kFALSE;
+Bool_t AliAnalysisTrackingUncertainties::IsKaon(const AliESDtrack * const tr, Bool_t /*useTPCPTOF*/) const {
+  //
+  // Selection of kaon candidates
+  // @TODO: To be implemented
+  //
+  Float_t nsigmaKaonTPC = fESDpid->NumberOfSigmasTPC(tr, AliPID::kKaon);
+  if (TMath::Abs(nsigmaKaonTPC) < 3) return kTRUE;
+  return kFALSE;
+
+}
+
+//________________________________________________________________________
+Bool_t AliAnalysisTrackingUncertainties::IsProton(const AliESDtrack * const tr, Bool_t /*useTPCPTOF*/) const{
+  // 
+  // Selection of proton candidates
+  // @TODO: To be implemented
+  //
+  Float_t nsigmaProtonTPC = fESDpid->NumberOfSigmasTPC(tr, AliPID::kProton);
+  if (TMath::Abs(nsigmaProtonTPC) < 3) return kTRUE;
+  return kFALSE;
+  
 }
 
 //________________________________________________________________________
@@ -576,7 +850,7 @@ void AliAnalysisTrackingUncertainties::InitializeTrackCutHistograms() {
   //
   //
   // (1.) number of clusters       
-  //                               0-ncl, 1-pt, 2-eta,         3-phi, 4-pid(0,unid,etc.)
+  //                               0-ncl, 1-pt, 2-eta,         3-phi, 4-pid(0-3 -> electron-proton, 4 -> undef, 5 -> all)
   Int_t    binsNcl[kNumberOfAxes] = { 40,   50,    20,            18,  6};
   Double_t minNcl[kNumberOfAxes]  = { 20,  0.1,    -1,             0, -0.5};
   Double_t maxNcl[kNumberOfAxes]  = {160,   20,    +1, 2*TMath::Pi(),  5.5};
@@ -594,7 +868,7 @@ void AliAnalysisTrackingUncertainties::InitializeTrackCutHistograms() {
   }
   //
   // (2.) chi2/cls-TPC            
-  //                                  0-chi2, 1-pt, 2-eta,         3-phi, 4-pid(0,unid,etc.)
+  //                                  0-chi2, 1-pt, 2-eta,         3-phi, 4-pid(0-3 -> electron-proton, 4 -> undef, 5 -> all)
   Int_t    binsChi2Tpc[kNumberOfAxes] = { 40,   50,    20,            18,  6};
   Double_t minChi2Tpc[kNumberOfAxes]  = {  0,  0.1,    -1,             0, -0.5};
   Double_t maxChi2Tpc[kNumberOfAxes]  = {  8,   20,    +1, 2*TMath::Pi(),  5.5};
@@ -602,7 +876,7 @@ void AliAnalysisTrackingUncertainties::InitializeTrackCutHistograms() {
   TString axisNameChi2Tpc[kNumberOfAxes]  = {"chi2tpc","pT","eta","phi","pid"};
   TString axisTitleChi2Tpc[kNumberOfAxes] = {"chi2tpc","pT","eta","phi","pid"};
   //
-  THnF * histChi2Tpc = new THnF("histChi2Tpc","number of clusters histogram",kNumberOfAxes, binsChi2Tpc, minChi2Tpc, maxChi2Tpc);
+  THnF * histChi2Tpc = new THnF("histChi2Tpc","chi2 per cls. in TPC",kNumberOfAxes, binsChi2Tpc, minChi2Tpc, maxChi2Tpc);
   BinLogAxis(histChi2Tpc, 1);
   fListHist->Add(histChi2Tpc);
   //
@@ -612,7 +886,7 @@ void AliAnalysisTrackingUncertainties::InitializeTrackCutHistograms() {
   }
   //
   // (3.) dca_z
-  //                               0-dcaZ, 1-pt, 2-eta,         3-phi, 4-pid(0,unid,etc.)
+  //                               0-dcaZ, 1-pt, 2-eta,         3-phi, 4-pid(0-3 -> electron-proton, 4 -> undef, 5 -> all)
   Int_t    binsDcaZ[kNumberOfAxes] = { 20,   50,    20,            18,  6};
   Double_t minDcaZ[kNumberOfAxes]  = {  0,  0.1,    -1,             0, -0.5};
   Double_t maxDcaZ[kNumberOfAxes]  = {  4,   20,    +1, 2*TMath::Pi(),  5.5};
@@ -620,7 +894,7 @@ void AliAnalysisTrackingUncertainties::InitializeTrackCutHistograms() {
   TString axisNameDcaZ[kNumberOfAxes]  = {"dcaZ","pT","eta","phi","pid"};
   TString axisTitleDcaZ[kNumberOfAxes] = {"dcaZ","pT","eta","phi","pid"};
   //
-  THnF * histDcaZ = new THnF("histDcaZ","number of clusters histogram",kNumberOfAxes, binsDcaZ, minDcaZ, maxDcaZ);
+  THnF * histDcaZ = new THnF("histDcaZ","dca_z to prim. vtx.",kNumberOfAxes, binsDcaZ, minDcaZ, maxDcaZ);
   BinLogAxis(histDcaZ, 1);
   fListHist->Add(histDcaZ);
   //
@@ -630,7 +904,7 @@ void AliAnalysisTrackingUncertainties::InitializeTrackCutHistograms() {
   }
   //
   // (4.) hit in SPD layer
-  //                              0-spdHit, 1-pt, 2-eta,         3-phi, 4-pid(0,unid,etc.)
+  //                              0-spdHit, 1-pt, 2-eta,         3-phi, 4-pid(0-3 -> electron-proton, 4 -> undef, 5 -> all)
   Int_t    binsSpd[kNumberOfAxes] = {    2,   50,    20,            18,  6};
   Double_t minSpd[kNumberOfAxes]  = { -0.5,  0.1,    -1,             0, -0.5};
   Double_t maxSpd[kNumberOfAxes]  = {  1.5,   20,    +1, 2*TMath::Pi(),  5.5};
@@ -638,7 +912,7 @@ void AliAnalysisTrackingUncertainties::InitializeTrackCutHistograms() {
   TString axisNameSpd[kNumberOfAxes]  = {"spdHit","pT","eta","phi","pid"};
   TString axisTitleSpd[kNumberOfAxes] = {"spdHit","pT","eta","phi","pid"};
   //
-  THnF * histSpd = new THnF("histSpd","number of clusters histogram",kNumberOfAxes, binsSpd, minSpd, maxSpd);
+  THnF * histSpd = new THnF("histSpd","hit in SPD layer or not",kNumberOfAxes, binsSpd, minSpd, maxSpd);
   BinLogAxis(histSpd, 1);
   fListHist->Add(histSpd);
   //
@@ -646,6 +920,96 @@ void AliAnalysisTrackingUncertainties::InitializeTrackCutHistograms() {
     histSpd->GetAxis(iaxis)->SetName(axisNameSpd[iaxis]);
     histSpd->GetAxis(iaxis)->SetTitle(axisTitleSpd[iaxis]);
   }
+  //
+  // (5.) number of crossed rows       
+  //                               0-ncr, 1-pt, 2-eta,         3-phi, 4-pid(0,unid,etc.)
+  Int_t    binsNcr[kNumberOfAxes] = { 40,   50,    20,            18,  6};
+  Double_t minNcr[kNumberOfAxes]  = { 20,  0.1,    -1,             0, -0.5};
+  Double_t maxNcr[kNumberOfAxes]  = {160,   20,    +1, 2*TMath::Pi(),  5.5};
+  //
+  TString axisNameNcr[kNumberOfAxes]  = {"Ncr","pT","eta","phi","pid"};
+  TString axisTitleNcr[kNumberOfAxes] = {"Ncr","pT","eta","phi","pid"};
+  //
+  THnF * histNcr = new THnF("histNcr","number of crossed rows TPC histogram",kNumberOfAxes, binsNcr, minNcr, maxNcr);
+  BinLogAxis(histNcr, 1);
+  fListHist->Add(histNcr);
+  //
+  for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
+    histNcr->GetAxis(iaxis)->SetName(axisNameNcr[iaxis]);
+    histNcr->GetAxis(iaxis)->SetTitle(axisTitleNcr[iaxis]);
+  }
+  //
+  // (6.) ratio crossed rows over findable clusters       
+  //                                0-CRoverFC, 1-pt, 2-eta,         3-phi, 4-pid(0,unid,etc.)
+  Int_t    binsCRoverFC[kNumberOfAxes] = {  26,   50,    20,            18,  6};
+  Double_t minCRoverFC[kNumberOfAxes]  = { 0.4,  0.1,    -1,             0, -0.5};
+  Double_t maxCRoverFC[kNumberOfAxes]  = { 1.8,   20,    +1, 2*TMath::Pi(),  5.5};
+  //
+  TString axisNameCRoverFC[kNumberOfAxes]  = {"CRoverFC","pT","eta","phi","pid"};
+  TString axisTitleCRoverFC[kNumberOfAxes] = {"CRoverFC","pT","eta","phi","pid"};
+  //
+  THnF * histCRoverFC = new THnF("histCRoverFC","number of crossed rows over findable clusters histogram",kNumberOfAxes, binsCRoverFC, minCRoverFC, maxCRoverFC);
+  BinLogAxis(histCRoverFC, 1);
+  fListHist->Add(histCRoverFC);
+  //
+  for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
+    histCRoverFC->GetAxis(iaxis)->SetName(axisNameCRoverFC[iaxis]);
+    histCRoverFC->GetAxis(iaxis)->SetTitle(axisTitleCRoverFC[iaxis]);
+  }
+  //
+  // (7.) max chi2 / ITS cluster       
+  //                               0-Chi2Its, 1-pt, 2-eta,         3-phi, 4-pid(0,unid,etc.)
+  Int_t    binsChi2Its[kNumberOfAxes] = { 25,   50,    20,            18,  6};
+  Double_t minChi2Its[kNumberOfAxes]  = {  0,  0.1,    -1,             0, -0.5};
+  Double_t maxChi2Its[kNumberOfAxes]  = { 50,   20,    +1, 2*TMath::Pi(),  5.5};
+  //
+  TString axisNameChi2Its[kNumberOfAxes]  = {"Chi2Its","pT","eta","phi","pid"};
+  TString axisTitleChi2Its[kNumberOfAxes] = {"Chi2Its","pT","eta","phi","pid"};
+  //
+  THnF * histChi2Its = new THnF("histChi2Its","number of crossed rows TPC histogram",kNumberOfAxes, binsChi2Its, minChi2Its, maxChi2Its);
+  BinLogAxis(histChi2Its, 1);
+  fListHist->Add(histChi2Its);
+  //
+  for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
+    histChi2Its->GetAxis(iaxis)->SetName(axisNameChi2Its[iaxis]);
+    histChi2Its->GetAxis(iaxis)->SetTitle(axisTitleChi2Its[iaxis]);
+  }
+  //
+  // (8.) tpc active volume length       
+  //                                0-TpcLength, 1-pt, 2-eta,         3-phi, 4-pid(0,unid,etc.)
+  Int_t    binsTpcLength[kNumberOfAxes] = {  40,   50,    20,            18,  6};
+  Double_t minTpcLength[kNumberOfAxes]  = {   0,  0.1,    -1,             0, -0.5};
+  Double_t maxTpcLength[kNumberOfAxes]  = { 170,   20,    +1, 2*TMath::Pi(),  5.5};
+  //
+  TString axisNameTpcLength[kNumberOfAxes]  = {"TpcLength","pT","eta","phi","pid"};
+  TString axisTitleTpcLength[kNumberOfAxes] = {"TpcLength","pT","eta","phi","pid"};
+  //
+  THnF * histTpcLength = new THnF("histTpcLength","number of crossed rows TPC histogram",kNumberOfAxes, binsTpcLength, minTpcLength, maxTpcLength);
+  BinLogAxis(histTpcLength, 1);
+  fListHist->Add(histTpcLength);
+  //
+  for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
+    histTpcLength->GetAxis(iaxis)->SetName(axisNameTpcLength[iaxis]);
+    histTpcLength->GetAxis(iaxis)->SetTitle(axisTitleTpcLength[iaxis]);
+  }
+  //
+  // (9.) match TPC->ITS
+  //                                  0-is matched, 1-pt, 2-eta,         3-phi, 4-pid(0-3 -> electron-proton, 4 -> undef, 5 -> all)
+  Int_t    binsTpcItsMatch[kNumberOfAxes] = {    2,   50,    20,            18,  6};
+  Double_t minTpcItsMatch[kNumberOfAxes]  = { -0.5,  0.1,    -1,             0, -0.5};
+  Double_t maxTpcItsMatch[kNumberOfAxes]  = {  1.5,   20,    +1, 2*TMath::Pi(),  5.5};
+  //
+  TString axisNameTpcItsMatch[kNumberOfAxes]  = {"isMatched","pT","eta","phi","pid"};
+  TString axisTitleTpcItsMatch[kNumberOfAxes] = {"isMatched","pT","eta","phi","pid"};
+  //
+  THnF * histTpcItsMatch = new THnF("histTpcItsMatch","TPC -> ITS matching",kNumberOfAxes, binsTpcItsMatch, minTpcItsMatch, maxTpcItsMatch);
+  BinLogAxis(histTpcItsMatch, 1);
+  fListHist->Add(histTpcItsMatch);
+  //
+  for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
+    histTpcItsMatch->GetAxis(iaxis)->SetName(axisNameTpcItsMatch[iaxis]);
+    histTpcItsMatch->GetAxis(iaxis)->SetTitle(axisTitleTpcItsMatch[iaxis]);
+  }