- ATO-57, ATO-71 calibResidual: Added tree output to study local track density effect...

[u/mrichter/AliRoot.git] / PWGPP / TPC / AliTPCcalibResidualPID.cxx

/**************************************************************************
* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
*                                                                        *
* Author: Yvonne Pachmayer <pachmay@physi.uni-heidelberg.de>             *
* Contributors are mentioned in the code where appropriate.              *
*                                                                        *
* Permission to use, copy, modify and distribute this software and its   *
* documentation strictly for non-commercial purposes is hereby granted   *
* without fee, provided that the above copyright notice appears in all   *
* copies and that both the copyright notice and this permission notice   *
* appear in the supporting documentation. The authors make no claims     *
* about the suitability of this software for any purpose. It is          *
* provided "as is" without express or implied warranty.                  *
**************************************************************************/
//
// The task:
// stores TPC PID quantities in a THnSparse
// output can then be used for e.g. dEdx calibration
//
// Author:
// Yvonne Pachmayer <pachmay@physi.uni-heidelberg.de>
//


#include "AliTPCcalibResidualPID.h"
#include "TChain.h"
#include "AliAnalysisManager.h"
#include "AliESDEvent.h"
#include "AliAODEvent.h"
#include "AliMCEvent.h"
#include "AliMCParticle.h"
//#include "AliStack.h"
#include "AliPID.h"
#include "AliTPCcalibResidualPID.h"
#include "AliESDtrack.h"
#include "AliESDtrackCuts.h"
#include "AliPIDResponse.h"
#include "AliInputEventHandler.h"
#include "AliESDInputHandler.h"
#include "AliESDv0KineCuts.h"
#include "AliESDv0.h"
#include "AliCentrality.h"
#include "THnSparse.h"
#include "TH2D.h"
#include "TCanvas.h"
#include "TGraphErrors.h"
#include "TMultiGraph.h"
#include "TAxis.h"
#include "TFile.h"
#include "TSpline.h"
#include "TStyle.h"
#include "AliTPCdEdxInfo.h"
#include "TString.h"
#include "TFitResult.h"
#include "TH1F.h"
#include "TLegend.h"
#include "TVirtualFitter.h"
#include "TTree.h"

using namespace std;

ClassImp(AliTPCcalibResidualPID)

Double_t AliTPCcalibResidualPID::fgCutGeo = 1.;   
Double_t AliTPCcalibResidualPID::fgCutNcr = 0.85; 
Double_t AliTPCcalibResidualPID::fgCutNcl = 0.7;  

//________________________________________________________________________
AliTPCcalibResidualPID::AliTPCcalibResidualPID()
  : AliAnalysisTaskSE(), fESD(0), fMC(0), fOutputContainer(0), fESDtrackCuts(0), fESDtrackCutsV0(0), fPIDResponse(0),
    fNumEtaCorrReqErrorsIssued(0),
    fNumMultCorrReqErrorsIssued(0),
    fUseTPCCutMIGeo(kFALSE),
    fUseMCinfo(kTRUE),
    fIsPbpOrpPb(kFALSE),
    fZvtxCutEvent(9999.0),
    fV0KineCuts(0x0),
    fCutOnProdRadiusForV0el(kTRUE),
    fNumTagsStored(0),
    fV0tags(0x0),
    fV0motherIndex(0x0),
    fV0motherPDG(0x0),
    fProduceAllPadTypes(0),
    fProduceGlobal(0),
    fProduceShortPads(0),
    fProduceMediumPads(0),
    fProduceLongPads(0),
    fProduceOroc(0),
    fHistPidQA(0), 
    fHistPidQAshort(0),
    fHistPidQAmedium(0),
    fHistPidQAlong(0),
    fHistPidQAoroc(0),
    fProduceTPCSignalSparse(0),
    fCorrectdEdxEtaDependence(0),
    fCorrectdEdxMultiplicityDependence(0),
    fThnspTpc(0),
    fWriteAdditionalOutput(kFALSE),
    fQAList(0x0),
    fhInvMassGamma(0x0),
    fhInvMassK0s(0x0),
    fhInvMassLambda(0x0),
    fhInvMassAntiLambda(0x0),
    fhArmenterosAll(0x0),
    fhArmenterosGamma(0x0),
    fhArmenterosK0s(0x0),
    fhArmenterosLambda(0x0),
    fhArmenterosAntiLambda(0x0),
    fHistSharedClusQAV0Pi(0x0),
    fHistSharedClusQAV0Pr(0x0),
    fHistSharedClusQAV0El(0x0),
    fTreeV0El(0x0),
    fTreeV0Pi(0x0),
    fTreeV0Pr(0x0),
    fTree_dEdx_tr(0.),
    fTree_dEdx_nb(0.),
    fTree_dEdx_vs(0.),
    fTree_dEdxExpected_tr(0.),
    fTree_p_TPC_tr(0.),
    fTree_p_TPC_nb(0.),
    fTree_p_TPC_vs(0.),
    fTree_BtimesChargeOverPt_tr(0.),
    fTree_BtimesChargeOverPt_nb(0.),
    fTree_BtimesChargeOverPt_vs(0.),
    fTree_tanTheta_tr(0.),
    fTree_tanTheta_nb(0.),
    fTree_tanTheta_vs(0.),
    fTree_distance_nb(0.),
    fTree_distance_vs(0.)
{
  // default Constructor
   /* fast compilation test
     gSystem->Load("libANALYSIS");
     gSystem->Load("libANALYSISalice");
     .L /lustre/alice/akalweit/train/trunk/util/statsQA/AliTPCcalibResidualPID.cxx++
   */

}


//________________________________________________________________________
AliTPCcalibResidualPID::AliTPCcalibResidualPID(const char *name)
  : AliAnalysisTaskSE(name), fESD(0), fMC(0), fOutputContainer(0), fESDtrackCuts(0), fESDtrackCutsV0(0), fPIDResponse(0),
    fNumEtaCorrReqErrorsIssued(0),
    fNumMultCorrReqErrorsIssued(0),
    fUseTPCCutMIGeo(kFALSE),
    fUseMCinfo(kTRUE),
    fIsPbpOrpPb(kFALSE),
    fZvtxCutEvent(9999.0),
    fV0KineCuts(0x0),
    fCutOnProdRadiusForV0el(kTRUE),
    fNumTagsStored(0),
    fV0tags(0x0),
    fV0motherIndex(0x0),
    fV0motherPDG(0x0),
    fProduceAllPadTypes(0),
    fProduceGlobal(0),
    fProduceShortPads(0),
    fProduceMediumPads(0),
    fProduceLongPads(0),
    fProduceOroc(0),
    fHistPidQA(0),
    fHistPidQAshort(0),
    fHistPidQAmedium(0),
    fHistPidQAlong(0), 
    fHistPidQAoroc(0),
    fProduceTPCSignalSparse(0),
    fCorrectdEdxEtaDependence(0),
    fCorrectdEdxMultiplicityDependence(0),
    fThnspTpc(0),
    fWriteAdditionalOutput(kFALSE),
    fQAList(0x0),
    fhInvMassGamma(0x0),
    fhInvMassK0s(0x0),
    fhInvMassLambda(0x0),
    fhInvMassAntiLambda(0x0),
    fhArmenterosAll(0x0),
    fhArmenterosGamma(0x0),
    fhArmenterosK0s(0x0),
    fhArmenterosLambda(0x0),
    fhArmenterosAntiLambda(0x0),
    fHistSharedClusQAV0Pi(0x0),
    fHistSharedClusQAV0Pr(0x0),
    fHistSharedClusQAV0El(0x0),
    fTreeV0El(0x0),
    fTreeV0Pi(0x0),
    fTreeV0Pr(0x0),
    fTree_dEdx_tr(0.),
    fTree_dEdx_nb(0.),
    fTree_dEdx_vs(0.),
    fTree_dEdxExpected_tr(0.),
    fTree_p_TPC_tr(0.),
    fTree_p_TPC_nb(0.),
    fTree_p_TPC_vs(0.),
    fTree_BtimesChargeOverPt_tr(0.),
    fTree_BtimesChargeOverPt_nb(0.),
    fTree_BtimesChargeOverPt_vs(0.),
    fTree_tanTheta_tr(0.),
    fTree_tanTheta_nb(0.),
    fTree_tanTheta_vs(0.),
    fTree_distance_nb(0.),
    fTree_distance_vs(0.)
{

  //fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE);
  //fESDtrackCutsV0 = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kFALSE);

  // Constructor
  DefineInput(0, TChain::Class());
  DefineOutput(1, TObjArray::Class());
  DefineOutput(2, TObjArray::Class());
  DefineOutput(3, TTree::Class());
  DefineOutput(4, TTree::Class());
  DefineOutput(5, TTree::Class());

}


//_________________________________________________
AliTPCcalibResidualPID::~AliTPCcalibResidualPID()
{
  delete fOutputContainer;
  fOutputContainer = 0;
  
  delete fQAList;
  fQAList = 0;
  
  /*
  delete fTreeV0El;
  fTreeV0El = 0;
  
  delete fTreeV0Pi;
  fTreeV0Pi = 0;
  
  delete fTreeV0Pr;
  fTreeV0Pr = 0;
  */
  
  delete fESDtrackCuts;
  fESDtrackCuts = 0;
  
  delete fESDtrackCutsV0;
  fESDtrackCutsV0 = 0;
  
  delete fV0KineCuts;
  fV0KineCuts = 0;
  
  delete fV0tags;
  fV0tags = 0;
  fNumTagsStored = 0;
  
  delete fV0motherIndex;
  fV0motherIndex = 0;
  
  delete fV0motherPDG;
  fV0motherPDG = 0;
}


//________________________________________________________________________
void AliTPCcalibResidualPID::UserCreateOutputObjects()
{
  AliInputEventHandler* inputHandler = dynamic_cast<AliInputEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
  if (!inputHandler)
    printf("Inputhandler not available \n");
  else
    fPIDResponse = inputHandler->GetPIDResponse();
  
  // THnSparse binning
  const Int_t kNdim = 9;
  //v0id,  dEdx,  TPCsigele,  TPCsigpion,  TOFbit,  eta, TPCclus, centr, p
  Int_t bins[kNdim] =    {    4,    250,    200,    200,          8,    20,        50,   20,   100};
  Double_t xmin[kNdim] = {  0.5,     30,   -10.,   -10.,       -1.5,   -1.,       60.,   0.,   0.1};
  Double_t xmax[kNdim] = {  4.5,   500.,    10.,    10.,        6.5,    1.,      160.,  100,   4};
  fThnspTpc= new THnSparseF("tpcsignals", "TPC signal;v0id;tpc signal;tpcsigele;tpcsigpion;tofbit;eta;tpcclus;centr;p (GeV/c)", kNdim, bins, xmin, xmax);
  BinLogAxis(fThnspTpc, 8);

  //
  // 0.ptot, 1.tpcSig, 2.particle ID, 3. assumed particle, 4. nSigmaTPC (4x), 5. nSigmaTOF (4x), 6. centrality
  // Concerning 2 (particle ID): 
  // - in case of MC: Contains MC ID. Bin 1-4 (<=> Slot 0-3): el, pi, ka, pr
  // - in case of data: Contains V0 particles + bin with all others: Bin 1-4 (<=> Slot 0-3): All non-V0s, V0-el, V0-pi, V0-pr
  //
      
  Int_t    binsHistQA[7] = {135, 1980,    4,    5, 40, 10,   40};
  Double_t xminHistQA[7] = {0.1,   20, -0.5, -0.5, -10, -5,   0.};
  Double_t xmaxHistQA[7] = {50., 2000,  3.5,  4.5,  10,  5, 20000};
  fHistPidQA = new THnSparseF("fHistPidQA","PID QA",7,binsHistQA,xminHistQA,xmaxHistQA);
  BinLogAxis(fHistPidQA, 0);

  //
  fHistPidQAshort  = new THnSparseF("fHistPidQAshort" ,"PID QA -- short pads",7,binsHistQA,xminHistQA,xmaxHistQA);
  fHistPidQAmedium = new THnSparseF("fHistPidQAmedium","PID QA -- med pads",7,binsHistQA,xminHistQA,xmaxHistQA);
  fHistPidQAlong   = new THnSparseF("fHistPidQAlong"  ,"PID QA -- long pads",7,binsHistQA,xminHistQA,xmaxHistQA);
  fHistPidQAoroc   = new THnSparseF("fHistPidQAoroc"  ,"PID QA -- oroc",7,binsHistQA,xminHistQA,xmaxHistQA);
  BinLogAxis(fHistPidQAshort, 0);
  BinLogAxis(fHistPidQAmedium, 0);
  BinLogAxis(fHistPidQAlong, 0);
  BinLogAxis(fHistPidQAoroc, 0);
  //
  fOutputContainer = new TObjArray(2);
  fOutputContainer->SetName(GetName());
  fOutputContainer->SetOwner();

  if(fProduceTPCSignalSparse)fOutputContainer->Add(fThnspTpc);
  if(fProduceGlobal)fOutputContainer->Add(fHistPidQA);
  if(fProduceAllPadTypes || fProduceShortPads)fOutputContainer->Add(fHistPidQAshort);
  if(fProduceAllPadTypes || fProduceMediumPads)fOutputContainer->Add(fHistPidQAmedium);
  if(fProduceAllPadTypes || fProduceLongPads)fOutputContainer->Add(fHistPidQAlong);
  if(fProduceAllPadTypes || fProduceOroc)fOutputContainer->Add(fHistPidQAoroc);
  
  // V0 Kine cuts 
  fV0KineCuts = new AliESDv0KineCuts;
  fV0KineCuts->SetGammaCutChi2NDF(5.);
  
  if (fCutOnProdRadiusForV0el) {
    // Only accept V0el with prod. radius within 45 cm -> PID will by systematically biased for larger values!
    Float_t gammaProdVertexRadiusCuts[2] = { 3.0, 45. }; 
    fV0KineCuts->SetGammaCutVertexR(&gammaProdVertexRadiusCuts[0]);
  }
  
  
  if (fWriteAdditionalOutput) {
    fQAList = new TObjArray(4);
    fQAList->SetName(GetName());
    fQAList->SetOwner();
    
    fhInvMassGamma      = new TH1F("fhInvMassGamma", "Invariant Mass of gammas; m_{ee} (GeV/#it{c}^{2}); Entries", 200, 0., 0.2);
    fhInvMassK0s        = new TH1F("fhInvMassK0s", "Invariant Mass of K0s; m_{#pi#pi} (GeV/#it{c}^{2}); Entries;", 200, 0.45, 0.55);
    fhInvMassLambda     = new TH1F("fhInvMassLambda", "Invariant Mass of lambdas; m_{p#pi^{-}} (GeV/#it{c}^{2}); Entries", 200, 1.05, 1.15);
    fhInvMassAntiLambda = new TH1F("fhInvMassAntiLambda", "Invariant Mass of anti-lambdas; m_{#pi^{+}#bar{p}} (GeV/#it{c}^{2}); Entries",
                                  200, 1.05, 1.15);
    
    fQAList->Add(fhInvMassGamma);
    fQAList->Add(fhInvMassK0s);
    fQAList->Add(fhInvMassLambda);
    fQAList->Add(fhInvMassAntiLambda);
    
    
    fhArmenterosAll = new TH2F("fhArmenterosAll",
                              "Armenteros plot all V0s;#alpha = (#it{p}^{+}_{L}-#it{p}^{-}_{L})/(#it{p}^{+}_{L}+#it{p}^{-}_{L});#it{q}_{T} (GeV/#it{c})",
                              200, -1., 1., 200, 0., 0.4);
    fhArmenterosGamma = new TH2F("fhArmenterosGamma",
                                "Armenteros plot Gamma;#alpha = (#it{p}^{+}_{L}-#it{p}^{-}_{L})/(#it{p}^{+}_{L}+#it{p}^{-}_{L});#it{q}_{T} (GeV/#it{c})",
                                200, -1., 1., 200, 0., 0.4);
    fhArmenterosK0s = new TH2F("fhArmenterosK0s",
                              "Armenteros plot K0s;#alpha = (#it{p}^{+}_{L}-#it{p}^{-}_{L})/(#it{p}^{+}_{L}+#it{p}^{-}_{L});#it{q}_{T} (GeV/#it{c})",
                              200, -1., 1., 200, 0., 0.4);
    fhArmenterosLambda = new TH2F("fhArmenterosLambda",
                                "Armenteros plot lambda;#alpha = (#it{p}^{+}_{L}-#it{p}^{-}_{L})/(#it{p}^{+}_{L}+#it{p}^{-}_{L});#it{q}_{T} (GeV/#it{c})",
                                200, -1., 1., 200, 0., 0.4);
    fhArmenterosAntiLambda = new TH2F("fhArmenterosAntiLambda",
                                      "Armenteros plot anti-lambda;#alpha = (#it{p}^{+}_{L}-#it{p}^{-}_{L})/(#it{p}^{+}_{L}+#it{p}^{-}_{L});#it{q}_{T} (GeV/#it{c})",
                                      200, -1., 1., 200, 0., 0.4);

    fQAList->Add(fhArmenterosAll);
    fQAList->Add(fhArmenterosGamma);
    fQAList->Add(fhArmenterosK0s);
    fQAList->Add(fhArmenterosLambda);
    fQAList->Add(fhArmenterosAntiLambda);
    
    
    const Int_t dimQASharedClusters = 4;
    const TString axisTitles[dimQASharedClusters] = { "#it{p}_{TPC} (GeV/#it{c})", "#it{#Delta'}", "#it{N}_{shared cl}", "Pad row"};
    Int_t    binsHistQASharedClusters[dimQASharedClusters] = { 100,  100, 160, 160};
    Double_t xminHistQASharedClusters[dimQASharedClusters] = { 0.3, 0.5,  0,   -1};
    Double_t xmaxHistQASharedClusters[dimQASharedClusters] = { 20,  1.5, 160, 159};
    
    fHistSharedClusQAV0Pi = new THnSparseF("fHistSharedClusQAV0Pi" ,"PID QA shared clusters - V0 pi", dimQASharedClusters,
                                          binsHistQASharedClusters, xminHistQASharedClusters, xmaxHistQASharedClusters);
    BinLogAxis(fHistSharedClusQAV0Pi, 0);
    for (Int_t i = 0; i < dimQASharedClusters; i++)
      fHistSharedClusQAV0Pi->GetAxis(i)->SetTitle(axisTitles[i].Data());
    fQAList->Add(fHistSharedClusQAV0Pi);
    
    fHistSharedClusQAV0Pr = new THnSparseF("fHistSharedClusQAV0Pr" ,"PID QA shared clusters - V0 pr", dimQASharedClusters,
                                          binsHistQASharedClusters, xminHistQASharedClusters, xmaxHistQASharedClusters);
    BinLogAxis(fHistSharedClusQAV0Pr, 0);
    for (Int_t i = 0; i < dimQASharedClusters; i++)
      fHistSharedClusQAV0Pi->GetAxis(i)->SetTitle(axisTitles[i].Data());
    fQAList->Add(fHistSharedClusQAV0Pr);

    fHistSharedClusQAV0El = new THnSparseF("fHistSharedClusQAV0El" ,"PID QA shared clusters - V0 el", dimQASharedClusters,
                                          binsHistQASharedClusters, xminHistQASharedClusters, xmaxHistQASharedClusters);
    BinLogAxis(fHistSharedClusQAV0El, 0);
    for (Int_t i = 0; i < dimQASharedClusters; i++)
      fHistSharedClusQAV0Pi->GetAxis(i)->SetTitle(axisTitles[i].Data());
    fQAList->Add(fHistSharedClusQAV0El);
    
    OpenFile(3);
    fTreeV0El = new TTree("treeV0El", "V0 el together with info from closest neighbour");
    fTreeV0El->Branch("dEdx_tr", &fTree_dEdx_tr);
    fTreeV0El->Branch("dEdx_nb", &fTree_dEdx_nb);
    fTreeV0El->Branch("dEdx_vs", &fTree_dEdx_vs);
    fTreeV0El->Branch("dEdxExpected_tr", &fTree_dEdxExpected_tr);
    fTreeV0El->Branch("p_TPC_tr", &fTree_p_TPC_tr);
    fTreeV0El->Branch("p_TPC_nb", &fTree_p_TPC_nb);
    fTreeV0El->Branch("p_TPC_vs", &fTree_p_TPC_vs);
    fTreeV0El->Branch("BtimesChargeOverPt_tr", &fTree_BtimesChargeOverPt_tr);
    fTreeV0El->Branch("BtimesChargeOverPt_nb", &fTree_BtimesChargeOverPt_nb);
    fTreeV0El->Branch("BtimesChargeOverPt_vs", &fTree_BtimesChargeOverPt_vs);
    fTreeV0El->Branch("tanTheta_tr", &fTree_tanTheta_tr);
    fTreeV0El->Branch("tanTheta_nb", &fTree_tanTheta_nb);
    fTreeV0El->Branch("tanTheta_vs", &fTree_tanTheta_vs);
    fTreeV0El->Branch("distance_nb", &fTree_distance_nb);
    fTreeV0El->Branch("distance_vs", &fTree_distance_vs);
    
    OpenFile(4);
    fTreeV0Pi = new TTree("treeV0Pi", "V0 pi together with info from closest neighbour");
    fTreeV0Pi->Branch("dEdx_tr", &fTree_dEdx_tr);
    fTreeV0Pi->Branch("dEdx_nb", &fTree_dEdx_nb);
    fTreeV0Pi->Branch("dEdx_vs", &fTree_dEdx_vs);
    fTreeV0Pi->Branch("dEdxExpected_tr", &fTree_dEdxExpected_tr);
    fTreeV0Pi->Branch("p_TPC_tr", &fTree_p_TPC_tr);
    fTreeV0Pi->Branch("p_TPC_nb", &fTree_p_TPC_nb);
    fTreeV0Pi->Branch("p_TPC_vs", &fTree_p_TPC_vs);
    fTreeV0Pi->Branch("BtimesChargeOverPt_tr", &fTree_BtimesChargeOverPt_tr);
    fTreeV0Pi->Branch("BtimesChargeOverPt_nb", &fTree_BtimesChargeOverPt_nb);
    fTreeV0Pi->Branch("BtimesChargeOverPt_vs", &fTree_BtimesChargeOverPt_vs);
    fTreeV0Pi->Branch("tanTheta_tr", &fTree_tanTheta_tr);
    fTreeV0Pi->Branch("tanTheta_nb", &fTree_tanTheta_nb);
    fTreeV0Pi->Branch("tanTheta_vs", &fTree_tanTheta_vs);
    fTreeV0Pi->Branch("distance_nb", &fTree_distance_nb);
    fTreeV0Pi->Branch("distance_vs", &fTree_distance_vs);
    
    OpenFile(5);
    fTreeV0Pr = new TTree("treeV0Pr", "V0 pr together with info from closest neighbour");
    fTreeV0Pr->Branch("dEdx_tr", &fTree_dEdx_tr);
    fTreeV0Pr->Branch("dEdx_nb", &fTree_dEdx_nb);
    fTreeV0Pr->Branch("dEdx_vs", &fTree_dEdx_vs);
    fTreeV0Pr->Branch("dEdxExpected_tr", &fTree_dEdxExpected_tr);
    fTreeV0Pr->Branch("p_TPC_tr", &fTree_p_TPC_tr);
    fTreeV0Pr->Branch("p_TPC_nb", &fTree_p_TPC_nb);
    fTreeV0Pr->Branch("p_TPC_vs", &fTree_p_TPC_vs);
    fTreeV0Pr->Branch("BtimesChargeOverPt_tr", &fTree_BtimesChargeOverPt_tr);
    fTreeV0Pr->Branch("BtimesChargeOverPt_nb", &fTree_BtimesChargeOverPt_nb);
    fTreeV0Pr->Branch("BtimesChargeOverPt_vs", &fTree_BtimesChargeOverPt_vs);
    fTreeV0Pr->Branch("tanTheta_tr", &fTree_tanTheta_tr);
    fTreeV0Pr->Branch("tanTheta_nb", &fTree_tanTheta_nb);
    fTreeV0Pr->Branch("tanTheta_vs", &fTree_tanTheta_vs);
    fTreeV0Pr->Branch("distance_nb", &fTree_distance_nb);
    fTreeV0Pr->Branch("distance_vs", &fTree_distance_vs);
  }
  
  PostData(1,fOutputContainer);
  if (fWriteAdditionalOutput) {
    PostData(2,fQAList);
    PostData(3,fTreeV0El);
    PostData(4,fTreeV0Pi);
    PostData(5,fTreeV0Pr);
  }
}


//_____________________________________________________________________________
void AliTPCcalibResidualPID::UserExec(Option_t *)
{
    //calls the Process function
    AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());

    if (!esdH) {
      printf("ERROR: Could not get ESDInputHandler \n");
    }
    else fESD = esdH->GetEvent();
    
    // If MC, forward MCevent
    fMC = dynamic_cast<AliMCEvent*>(MCEvent());
    //
    Process(fESD, fMC);
    //
    PostData(1,fOutputContainer);
    if (fWriteAdditionalOutput) {
      PostData(2,fQAList);
      PostData(3,fTreeV0El);
      PostData(4,fTreeV0Pi);
      PostData(5,fTreeV0Pr);
    }
}



//________________________________________________________________________
void AliTPCcalibResidualPID::Process(AliESDEvent *const esdEvent, AliMCEvent *const mcEvent)
{
  
  //called for each event
  if (!esdEvent) {
    Printf("ERROR: esdEvent not available"); 
    return;
  }

  if (!fPIDResponse || !fV0KineCuts) {
    Printf("ERROR: No PIDresponse and/or v0KineCuts!");
    return;
  }

  Float_t centralityFper=99;

  AliCentrality *esdCentrality = esdEvent->GetCentrality();
  centralityFper = esdCentrality->GetCentralityPercentile("V0M");

  if (!GetVertexIsOk(esdEvent))
    return;
  
  const AliESDVertex* fESDvertex = esdEvent->GetPrimaryVertexTracks(); 
  if (!fESDvertex)
    return;
  
  Int_t ncontr = fESDvertex->GetNContributors();
  
  if (ncontr <= 0)
    return;

  // Fill V0 arrays with V0s
  FillV0PIDlist(esdEvent);
  
  // Array with flags wheter QA for this V0 was already done or not
  const Int_t numV0s = esdEvent->GetNumberOfV0s();
  Bool_t v0QAadded[numV0s];
  for (Int_t i = 0; i < numV0s; i++)
    v0QAadded[i] = kFALSE;
  
  Int_t nTotTracks = esdEvent->GetNumberOfTracks();
  const Int_t nTotESDTracks = esdEvent->GetNumberOfESDTracks();
  
  const Double_t magField = esdEvent->GetMagneticField();
  
  Bool_t etaCorrAvail = fPIDResponse->UseTPCEtaCorrection();
  Bool_t multCorrAvail = fPIDResponse->UseTPCMultiplicityCorrection();
  
  for (Int_t iTracks = 0; iTracks < nTotTracks; iTracks++){//begin track loop 
    AliESDtrack *trackESD = esdEvent->GetTrack(iTracks);
    if(!trackESD) {
      Printf("ERROR: Could not receive track %d (esd loop)", iTracks);
      continue;
    }
    if((TMath::Abs(trackESD->Eta())) > 0.9)
      continue;
    
     // Do not distinguish positively and negatively charged V0's
    Char_t v0tagAllCharges = TMath::Abs(GetV0tag(iTracks));
    if (v0tagAllCharges == -99) {
      AliError(Form("Problem with V0 tag list (requested V0 track for track %d from %d (list status %d))!", iTracks, esdEvent->GetNumberOfTracks(),
                    fV0tags != 0x0));
      continue;
    }
    
    Bool_t isV0el = v0tagAllCharges == 14;
    Bool_t isV0pi = v0tagAllCharges == 15;
    Bool_t isV0pr = v0tagAllCharges == 16;
    Bool_t isV0 = isV0el || isV0pi || isV0pr;
    
    if (mcEvent && fUseMCinfo) {
      // For MC, do not look for V0s, i.e. demand the non-V0 track cuts
      if (fESDtrackCuts && !fESDtrackCuts->AcceptTrack(trackESD)) continue;
      
      if (fUseTPCCutMIGeo) {
        // If cut on geometry is active, don't cut on number of clusters, since such a cut is already included.
        if (!TPCCutMIGeo(trackESD, esdEvent))
          continue;
      }
      else {
        // If cut on geometry is not active, always cut on num clusters
        if (trackESD->GetTPCsignalN() < 60)
          continue;
      }
    }
    else {
      // For data, take V0 AND non-V0 with separate cuts
      //if (!isV0 && fESDtrackCuts && !fESDtrackCuts->AcceptTrack(trackESD)) continue;
      if (isV0) {
        if (fESDtrackCutsV0 && !fESDtrackCutsV0->AcceptTrack(trackESD))
          continue;
      }
      else {
        if (fESDtrackCuts && !fESDtrackCuts->AcceptTrack(trackESD))
          continue;
      }
        
      if (fUseTPCCutMIGeo) {
        // If cut on geometry is active, don't cut on number of clusters, since such a cut is already included.
        if (!isV0) {
          if (!TPCCutMIGeo(trackESD, esdEvent))
            continue;
        }
        else {
          // One should not cut on geometry for V0's since they have different topology. (Loosely) cut on num of clusters instead.
          if (trackESD->GetTPCsignalN() < 60)
            continue;
        }
      }
      else {
        // If cut on geometry is not active, always cut on num clusters
        if (trackESD->GetTPCsignalN() < 60)
          continue;
      }
    }

    const AliExternalTrackParam *paramIn = trackESD->GetInnerParam();
    Float_t precin=-1;
    if(paramIn) {
      precin=paramIn->GetP();
    } else {
      printf("Inner param not found for track %d - skipping!\n", iTracks);
      continue;
    }
    Int_t precdefault=CompareFloat(precin,-1);
    if(precdefault==1) continue; // momentum cut
    //
    
    AliMCParticle* mcTrack = 0x0;
    Int_t particleID = -1;
    Int_t v0id = 0;
    
    if (mcEvent && fUseMCinfo) {
      // MC - particle ID = MC ID
      Int_t label = trackESD->GetLabel();
      
      if (label < 0)
        continue; // If MC is available, reject tracks with negative ESD label
        
      mcTrack = dynamic_cast<AliMCParticle*>(mcEvent->GetTrack(TMath::Abs(label)));
      if (!mcTrack) {
        Printf("ERROR: Could not receive mcTrack with label %d for track %d", label, iTracks);
        continue;
      }
      
      /*// Only accept MC primaries
      if (!mcEvent->Stack()->IsPhysicalPrimary(TMath::Abs(label))) {
        continue;
      }*/
      
      Int_t pdgAbs = TMath::Abs(mcTrack->PdgCode());
      
      if (pdgAbs == 11) { // electron
        particleID = 0;
      }
      else if (pdgAbs == 211) { // pion
        particleID = 1;
      }
      else if (pdgAbs == 321) { // kaon
        particleID = 2;
      }
      else if (pdgAbs == 2212) { // proton
        particleID = 3;
      }
      else
        continue; // Reject all other particles
    }
    else {
      // Data - particle ID = V0 ID (if V0)
      if (isV0pi) { // pion
        particleID = 2;
        v0id = 2;
      }
      else if (isV0el) { // electron
        particleID = 1;
        v0id = 1;
      }
      else if (isV0pr) { // proton
        particleID = 3;
        v0id = 3;
      }
      else { // No V0-ID available (species must be determined via TPC, TOF, ...)
        particleID = 0;
      }
    }


    Float_t tpcsignal=trackESD->GetTPCsignal();
    Int_t tofbit=-1;
    Int_t iTOFpid=0;
    Int_t ikTIME=0;
    Double_t tpcnsigmaele=-10;
    Double_t tpcnsigmapion=-10;
    //
    if ((trackESD->GetStatus() & AliESDtrack::kTOFpid)) iTOFpid = 1;
    if ((trackESD->GetStatus() & AliESDtrack::kTIME)) ikTIME = 1;
    Float_t time0 = fPIDResponse->GetTOFResponse().GetTimeZero();
    //
    if((iTOFpid==1) &&(ikTIME==1)){
      Double_t tofnsigmaele= fPIDResponse->NumberOfSigmasTOF(trackESD,AliPID::kElectron, time0);
      Double_t tofnsigmapion=fPIDResponse->NumberOfSigmasTOF(trackESD,AliPID::kPion, time0);
      if (TMath::Abs(tofnsigmapion)<3) tofbit = 1;
      if (TMath::Abs(tofnsigmaele)<3)  tofbit = 0;
      tpcnsigmaele=fPIDResponse->NumberOfSigmasTPC(trackESD,AliPID::kElectron);
      tpcnsigmapion=fPIDResponse->NumberOfSigmasTPC(trackESD,AliPID::kPion);
    }
    //
    Int_t tpcnclusN=trackESD->GetTPCsignalN();
    Double_t eta=trackESD->Eta();
    
    //
    if(fProduceTPCSignalSparse){
      Double_t contentSignal[9];
      contentSignal[0]=v0id;
      contentSignal[1]=tpcsignal;
      contentSignal[2]=tpcnsigmaele;
      contentSignal[3]=tpcnsigmapion;
      contentSignal[4]=tofbit;
      contentSignal[5]=eta;
      contentSignal[6]=tpcnclusN;
      contentSignal[7]=centralityFper;
      contentSignal[8]=precin;
      //
      if (fThnspTpc->GetEntries() < 1e6) fThnspTpc->Fill(contentSignal);
    }//should it be created or not



    //
    // 2nd THnSparse
    //
    Double_t tpcQA[5] = {fPIDResponse->NumberOfSigmasTPC(trackESD, AliPID::kElectron),
			 fPIDResponse->NumberOfSigmasTPC(trackESD, AliPID::kPion),
			 fPIDResponse->NumberOfSigmasTPC(trackESD, AliPID::kKaon),
			 fPIDResponse->NumberOfSigmasTPC(trackESD, AliPID::kProton),
			 0};
    Double_t tofQA[5] = {fPIDResponse->NumberOfSigmasTOF(trackESD, AliPID::kElectron, time0),
			 fPIDResponse->NumberOfSigmasTOF(trackESD, AliPID::kPion, time0),
			 fPIDResponse->NumberOfSigmasTOF(trackESD, AliPID::kKaon, time0),
			 fPIDResponse->NumberOfSigmasTOF(trackESD, AliPID::kProton, time0),
			 0 };
       
    // id 5 is just again Kaons in restricted eta range
    tpcQA[4] = tpcQA[2];
    tofQA[4] = tofQA[2];
    
    //
    // dE/dx in different pad regions
    //
    AliTPCdEdxInfo * infoTpcPid = trackESD->GetTPCdEdxInfo();
    Double32_t signal[4]; Char_t ncl[3]; Char_t nrows[3];
    if (infoTpcPid) {
      infoTpcPid->GetTPCSignalRegionInfo(signal, ncl, nrows);
    } else {
      for(Int_t iarr = 0; iarr < 3; iarr++) {
	signal[iarr] = 0;
	ncl[iarr] = 0;
	nrows[iarr] = 0;
      }
      signal[3] = 0;
    }
    //
    UInt_t status = trackESD->GetStatus();
    Bool_t hasTOFout  = status&AliESDtrack::kTOFout; 
    Bool_t hasTOFtime = status&AliESDtrack::kTIME;
    Bool_t hasTOF     = kFALSE;
    if (hasTOFout && hasTOFtime) hasTOF = kTRUE;
    Float_t length = trackESD->GetIntegratedLength();
    // Length check only for primaries!
    if (length < 350 && !isV0) hasTOF = kFALSE;
    //
    
    if (!hasTOF) {
      // Make sure that number of sigmas is large in this case, so that the track will be rejected if a TOF cut is applied
      for (Int_t i = 0; i < 5; i++) {
        tofQA[i] = 999;
      }
    }


    Double_t processedTPCsignal[5] = { tpcsignal, tpcsignal, tpcsignal, tpcsignal, tpcsignal };
  
    if (fCorrectdEdxEtaDependence && fNumEtaCorrReqErrorsIssued < 23 && !etaCorrAvail) {
      AliError("TPC eta correction requested, but was not activated in PID response (most likely not available)!");
      fNumEtaCorrReqErrorsIssued++;
      if (fNumEtaCorrReqErrorsIssued == 23)
        AliError("Ignoring this error from now on!");
    }
    
    if (fCorrectdEdxMultiplicityDependence && fNumMultCorrReqErrorsIssued < 23 && !multCorrAvail) {
      AliError("TPC multiplicity correction requested, but was not activated in PID response (most likely not available)!");
      fNumMultCorrReqErrorsIssued++;
      if (fNumMultCorrReqErrorsIssued == 23)
        AliError("Ignoring this error from now on!");
    }
    
    if (fCorrectdEdxEtaDependence && etaCorrAvail && fCorrectdEdxMultiplicityDependence && multCorrAvail) {
      processedTPCsignal[0] = fPIDResponse->GetTPCResponse().GetEtaAndMultiplicityCorrectedTrackdEdx(trackESD, AliPID::kElectron,
                                                                                                     AliTPCPIDResponse::kdEdxDefault);
      processedTPCsignal[1] = fPIDResponse->GetTPCResponse().GetEtaAndMultiplicityCorrectedTrackdEdx(trackESD, AliPID::kPion,
                                                                                                     AliTPCPIDResponse::kdEdxDefault);
      processedTPCsignal[2] = fPIDResponse->GetTPCResponse().GetEtaAndMultiplicityCorrectedTrackdEdx(trackESD, AliPID::kKaon,
                                                                                                     AliTPCPIDResponse::kdEdxDefault);
      processedTPCsignal[3] = fPIDResponse->GetTPCResponse().GetEtaAndMultiplicityCorrectedTrackdEdx(trackESD, AliPID::kProton,
                                                                                                     AliTPCPIDResponse::kdEdxDefault);
    }
    else if (fCorrectdEdxEtaDependence && etaCorrAvail) {
      processedTPCsignal[0] = fPIDResponse->GetTPCResponse().GetEtaCorrectedTrackdEdx(trackESD, AliPID::kElectron,
                                                                                      AliTPCPIDResponse::kdEdxDefault);
      processedTPCsignal[1] = fPIDResponse->GetTPCResponse().GetEtaCorrectedTrackdEdx(trackESD, AliPID::kPion,
                                                                                      AliTPCPIDResponse::kdEdxDefault);
      processedTPCsignal[2] = fPIDResponse->GetTPCResponse().GetEtaCorrectedTrackdEdx(trackESD, AliPID::kKaon,
                                                                                      AliTPCPIDResponse::kdEdxDefault);
      processedTPCsignal[3] = fPIDResponse->GetTPCResponse().GetEtaCorrectedTrackdEdx(trackESD, AliPID::kProton,
                                                                                      AliTPCPIDResponse::kdEdxDefault);
    }
    else if (fCorrectdEdxMultiplicityDependence && multCorrAvail) {
      processedTPCsignal[0] = fPIDResponse->GetTPCResponse().GetMultiplicityCorrectedTrackdEdx(trackESD, AliPID::kElectron,
                                                                                               AliTPCPIDResponse::kdEdxDefault);
      processedTPCsignal[1] = fPIDResponse->GetTPCResponse().GetMultiplicityCorrectedTrackdEdx(trackESD, AliPID::kPion,
                                                                                               AliTPCPIDResponse::kdEdxDefault);
      processedTPCsignal[2] = fPIDResponse->GetTPCResponse().GetMultiplicityCorrectedTrackdEdx(trackESD, AliPID::kKaon,
                                                                                               AliTPCPIDResponse::kdEdxDefault);
      processedTPCsignal[3] = fPIDResponse->GetTPCResponse().GetMultiplicityCorrectedTrackdEdx(trackESD, AliPID::kProton,
                                                                                               AliTPCPIDResponse::kdEdxDefault);
    }
    
    // id 5 is just again Kaons in restricted eta range
    processedTPCsignal[4] = processedTPCsignal[2];
    
    for(Int_t iPart = 0; iPart < 5; iPart++) {
      // Only accept "Kaons" within |eta| < 0.2 for index 4 in case of data (no contamination in case of MC, so this index is not used)
      if (iPart == 4 && ((mcEvent && fUseMCinfo) || abs(trackESD->Eta()) > 0.2))
        continue;
      
      Double_t vecHistQA[7] = {precin, processedTPCsignal[iPart], (Double_t)particleID, (Double_t)iPart, tpcQA[iPart], tofQA[iPart],
                               (Double_t)nTotESDTracks};
      if (fProduceGlobal) fHistPidQA->Fill(vecHistQA);
      vecHistQA[1] = signal[0]; vecHistQA[4] = ncl[0];
      if ((fProduceAllPadTypes || fProduceShortPads)) fHistPidQAshort->Fill(vecHistQA);
      //
      vecHistQA[1] = signal[1]; vecHistQA[4] = ncl[1];
      if ((fProduceAllPadTypes || fProduceMediumPads)) fHistPidQAmedium->Fill(vecHistQA);
      //
      vecHistQA[1] = signal[2]; vecHistQA[4] = ncl[2];
      if ((fProduceAllPadTypes || fProduceLongPads)) fHistPidQAlong->Fill(vecHistQA);
      //
      vecHistQA[1] = signal[3]; vecHistQA[4] = nrows[1] + nrows[2];
      if ((fProduceAllPadTypes || fProduceOroc)) fHistPidQAoroc->Fill(vecHistQA);      
    }
    
    if (fWriteAdditionalOutput && isV0) {
      
      // Find the closest neighbour and fill the information into the trees.
      fTree_p_TPC_tr = precin;
      fTree_dEdx_tr = tpcsignal;
      fTree_BtimesChargeOverPt_tr = TMath::Abs(magField) * trackESD->GetSigned1Pt();
      fTree_tanTheta_tr = trackESD->GetInnerParam()->GetTgl();
      
      Int_t closestNeighbourIndex = -1;
      const Double_t tpcInnerRadius = 85.;
      const Int_t parIndexGlobalPhi = 8;
      const Double_t z_tr = trackESD->GetInnerParam()->GetZ();
      const Double_t phi_tr = trackESD->GetInnerParam()->GetParameterAtRadius(tpcInnerRadius, magField, parIndexGlobalPhi);
      
      const Double_t distanceCut = 12.; // Cut on distance between track to closest neighbour
      
      Double_t z_nb = 9999;
      Double_t phi_nb = 9999;
      Double_t delta_z = 9999;
      Double_t delta_phi = 9999;
      
      fTree_distance_nb = -9999;
      
      for (Int_t jTracks = 0; jTracks < nTotTracks; jTracks++) {
        if (iTracks == jTracks)
          continue; // Exclude track we are looking at right now from neighbour list
        
        AliESDtrack *track_nb= esdEvent->GetTrack(jTracks);
        if(!track_nb)
          continue;
        
        const AliExternalTrackParam *paramIn_nb = track_nb->GetInnerParam();
        if (!paramIn_nb)
          continue;
        
        z_nb =  paramIn_nb->GetZ();
        delta_z = TMath::Abs(z_nb - z_tr);
        
        if (delta_z >= distanceCut)
          continue;
        
        phi_nb = paramIn_nb->GetParameterAtRadius(tpcInnerRadius, magField, parIndexGlobalPhi);
        delta_phi = TMath::Abs(phi_nb - phi_tr);
        
        const Double_t delta_r_phi = tpcInnerRadius * delta_phi;
        
        if (delta_r_phi >= distanceCut)
          continue;
        
        
        const Double_t distance = TMath::Sqrt(delta_z * delta_z + delta_r_phi * delta_r_phi);
        
        if (distance >= distanceCut)
          continue;
        
        if (closestNeighbourIndex < 0 || distance < fTree_distance_nb) {
          fTree_distance_nb = distance;
          closestNeighbourIndex = jTracks;
        }
      }
      
      if (closestNeighbourIndex >= 0 && closestNeighbourIndex < nTotTracks) {
        AliESDtrack *track_nb= esdEvent->GetTrack(closestNeighbourIndex);
        fTree_dEdx_nb = track_nb->GetTPCsignal();
        fTree_p_TPC_nb = track_nb->GetInnerParam()->GetP();
        fTree_BtimesChargeOverPt_nb = TMath::Abs(magField) * track_nb->GetSigned1Pt();
        fTree_tanTheta_nb = track_nb->GetInnerParam()->GetTgl();
      }
      else {
        // No neighbour in this range -> nevertheless store the track, but set neighbour values to invalid
        fTree_dEdx_nb = -999.;
        fTree_p_TPC_nb = -999.;
        fTree_BtimesChargeOverPt_nb = -999.;
        fTree_tanTheta_nb = -999.;
        fTree_distance_nb = -999.;
      }
      
      
      // Find the other V0 daughter and store its values
      const Int_t motherIndex = GetV0motherIndex(iTracks);
      AliESDv0* v0Mother = (AliESDv0*)esdEvent->GetV0(motherIndex);
      if (!v0Mother) {
        printf("Error: Track tagged as V0 daughter, but V0 mother not found!\n");
        continue;
      }
      
      const Int_t iTrackP = v0Mother->GetPindex();  // positive track
      const Int_t iTrackN = v0Mother->GetNindex();  // negative track_nb
      
      Int_t iTrackV0sister = -1;
      if (iTrackP == iTracks)
        iTrackV0sister = iTrackN;
      else if (iTrackN == iTracks)
        iTrackV0sister = iTrackP;
      else {
        printf("Error: V0 sister relations are wrong!\n");
        continue;
      }
      
      fTree_dEdx_vs = -999.;
      fTree_p_TPC_vs = -999.;
      fTree_BtimesChargeOverPt_vs = -999.;
      fTree_tanTheta_vs = -999.;
      fTree_distance_vs = -999.;
        
      AliESDtrack *track_vs = esdEvent->GetTrack(iTrackV0sister);
      if (track_vs) {
        const AliExternalTrackParam *paramIn_vs = track_vs->GetInnerParam();
        if (paramIn_vs) {
          fTree_dEdx_vs = track_vs->GetTPCsignal();
          fTree_p_TPC_vs = paramIn_vs->GetP();
          fTree_BtimesChargeOverPt_vs = TMath::Abs(magField) * track_vs->GetSigned1Pt();
          fTree_tanTheta_vs = paramIn_vs->GetTgl();
          
          // Calculate distance
          const Double_t z_vs =  paramIn_vs->GetZ();
          const Double_t delta_z_vs = TMath::Abs(z_vs - z_tr);
          
          const Double_t phi_vs = paramIn_vs->GetParameterAtRadius(tpcInnerRadius, magField, parIndexGlobalPhi);
          const Double_t delta_phi_vs = TMath::Abs(phi_vs - phi_tr);
          
          const Double_t delta_r_phi_vs = tpcInnerRadius * delta_phi_vs;
          
          fTree_distance_vs = TMath::Sqrt(delta_z_vs * delta_z_vs + delta_r_phi_vs * delta_r_phi_vs);
        }
      }
      
      if (isV0el) {
        fTree_dEdxExpected_tr = fPIDResponse->GetTPCResponse().GetExpectedSignal(trackESD, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault, 
                                                                                 fPIDResponse->UseTPCEtaCorrection(),
                                                                                 fPIDResponse->UseTPCMultiplicityCorrection());
        if (fTree_dEdxExpected_tr > 0)
          fTreeV0El->Fill();
      }
      else if (isV0pi) {
        fTree_dEdxExpected_tr = fPIDResponse->GetTPCResponse().GetExpectedSignal(trackESD, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault, 
                                                                                 fPIDResponse->UseTPCEtaCorrection(),
                                                                                 fPIDResponse->UseTPCMultiplicityCorrection());
        if (fTree_dEdxExpected_tr > 0)
          fTreeV0Pi->Fill();
      }
      else if (isV0pr) {
        fTree_dEdxExpected_tr = fPIDResponse->GetTPCResponse().GetExpectedSignal(trackESD, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault, 
                                                                                 fPIDResponse->UseTPCEtaCorrection(),
                                                                                 fPIDResponse->UseTPCMultiplicityCorrection());
        if (fTree_dEdxExpected_tr > 0)
          fTreeV0Pr->Fill();
      }
      
      // Fill QA hists for shared clusters dE/dx
      Double_t vecHistQA[4] = { precin, -1., (Double_t)trackESD->GetTPCSharedMap().CountBits(), -1. };
      THnSparse* currHist = fHistSharedClusQAV0Pi;
      
      if (isV0pi) {
        Double_t expectedDeDx = fPIDResponse->GetTPCResponse().GetExpectedSignal(trackESD, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault,
                                                                                 fPIDResponse->UseTPCEtaCorrection(),
                                                                                 fPIDResponse->UseTPCMultiplicityCorrection()); 
        if (expectedDeDx > 0 ) {
          vecHistQA[1] = tpcsignal / expectedDeDx;
          currHist = fHistSharedClusQAV0Pi;
        }
      }
      else if (isV0pr) {
        Double_t expectedDeDx = fPIDResponse->GetTPCResponse().GetExpectedSignal(trackESD, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault,
                                                                                 fPIDResponse->UseTPCEtaCorrection(),
                                                                                 fPIDResponse->UseTPCMultiplicityCorrection()); 
        if (expectedDeDx > 0 ) {
          vecHistQA[1] = tpcsignal / expectedDeDx;
          currHist = fHistSharedClusQAV0Pr;
        }
      }
      else if (isV0el) {
        Double_t expectedDeDx = fPIDResponse->GetTPCResponse().GetExpectedSignal(trackESD, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault,
                                                                                 fPIDResponse->UseTPCEtaCorrection(), 
                                                                                 fPIDResponse->UseTPCMultiplicityCorrection()); 
        if (expectedDeDx > 0 ) {
          vecHistQA[1] = tpcsignal / expectedDeDx;
          currHist = fHistSharedClusQAV0El;
        }
      }
      
      Int_t iRowInd = -1;
      // iRowInd == -1 for "all rows w/o multiple counting"
      currHist->Fill(vecHistQA);
      
      // Fill hist for every pad row with shared cluster
      for (iRowInd = 0; iRowInd < 159; iRowInd++) {
        if (trackESD->GetTPCSharedMap().TestBitNumber(iRowInd)) {
          vecHistQA[3] = iRowInd;
          currHist->Fill(vecHistQA);
        }
      }
      
      
      
      // Check, whether the QA for this V0 mother was already filled into the histograms (is the case if the other daughter has already
      // been processed). If not, set flag to kTRUE and fill the QA.
      const Int_t iV0 = GetV0motherIndex(iTracks);
      if (!v0QAadded[iV0]) {
        v0QAadded[iV0] = kTRUE;
      
        AliESDv0* esdV0 = (AliESDv0*)esdEvent->GetV0(iV0);
        
        if (!esdV0) {
          printf("Error: V0 tagged, but does not exist!\n");
        }
        else {
          Float_t armVar[2] = {0.0,0.0};
          fV0KineCuts->Armenteros(esdV0, armVar);
          
          const Int_t motherPDG = GetV0motherPDG(iTracks);
          
          // Mother and daughter match the requirements, otherwise it wouldn't be tagged as a V0. Now just fill the QA histos.
          fhArmenterosAll->Fill(armVar[0], armVar[1]);
          //if (esdV0->TestBit(BIT(14))) {
          if (TMath::Abs(motherPDG) == 22) {
            fhInvMassGamma->Fill(esdV0->GetEffMass(AliPID::kElectron, AliPID::kElectron));
            fhArmenterosGamma->Fill(armVar[0], armVar[1]);
          }
          else if (TMath::Abs(motherPDG) == 310) {
          //else if (esdV0->TestBit(BIT(15))) {
            fhInvMassK0s->Fill(esdV0->GetEffMass(AliPID::kPion, AliPID::kPion));
            fhArmenterosK0s->Fill(armVar[0], armVar[1]);
          }
          else if (motherPDG == 3122) {
          //else if (esdV0->TestBit(BIT(16))) {
            fhInvMassLambda->Fill(esdV0->GetEffMass(AliPID::kProton, AliPID::kPion));
            fhArmenterosLambda->Fill(armVar[0], armVar[1]);
          }
          else if (motherPDG == -3122) {
          //else if (esdV0->TestBit(BIT(17))) {
            fhInvMassAntiLambda->Fill(esdV0->GetEffMass(AliPID::kPion, AliPID::kProton));
            fhArmenterosAntiLambda->Fill(armVar[0], armVar[1]);
          }
        }
      }
    }
  } //end track loop 
  

  // Clear the V0 PID arrays
  ClearV0PIDlist();
}      


//________________________________________________________________________
Int_t  AliTPCcalibResidualPID::CompareFloat(Float_t f1, Float_t f2) const
{
    //compares if the Float_t f1 is equal with f2 and returns 1 if true and 0 if false
    Float_t precision = 0.00001;
    if (((f1 - precision) < f2) &&
	((f1 + precision) > f2))
    {
	return 1;
    }
    else
    {
	return 0;
    }
}


//________________________________________________________________________
void AliTPCcalibResidualPID::Terminate(const Option_t *)
{


}



//________________________________________________________________________
void AliTPCcalibResidualPID::BinLogAxis(const THnSparseF *h, Int_t axisNumber) {
  //
  // Method for the correct logarithmic binning of histograms
  //
  TAxis *axis = h->GetAxis(axisNumber);
  int bins = axis->GetNbins();

  Double_t from = axis->GetXmin();
  Double_t to = axis->GetXmax();
  Double_t *newBins = new Double_t[bins + 1];
   
  newBins[0] = from;
  Double_t factor = pow(to/from, 1./bins);
  
  for (int i = 1; i <= bins; i++) {
   newBins[i] = factor * newBins[i-1];
  }
  axis->Set(bins, newBins);
  delete [] newBins;
  
}


//________________________________________________________________________
Double_t* AliTPCcalibResidualPID::ExtractResidualPID(THnSparseF * histPidQA, const Bool_t useV0s, const Char_t * outFile,
                                                     const Char_t * type, const Char_t * period, const Char_t * pass,
                                                     const Char_t * system, const Double_t * initialParameters,
                                                     const Char_t *dedxtype,
                                                     AliTPCcalibResidualPID::FitType fitType) {
  //
  // (1.) obtain residual graphs
  //
  TObjArray * arr = 0x0;
  
  Bool_t isMC = kFALSE;
  if (strcmp(type, "MC") == 0) {
    arr = GetResidualGraphsMC(histPidQA, system);
    isMC = kTRUE;
  }
  else if (strcmp(type, "DATA") == 0) {
    arr = GetResidualGraphs(histPidQA, system, useV0s);
  }
  else {
    Printf("ERROR - ExtractResidualPID: Unknown type \"%s\" - must be \"MC\" or \"DATA\"!", type);
    
    return 0x0;
  }
  
  Bool_t isPPb = kFALSE;
  if (strcmp(system, "PPB") == 0 || strcmp(system, "PBP") == 0) {
    Printf("p-Pb/Pb-p detected - Applying special handling!");
    isPPb = kTRUE;
  }
  //
  // (2.) get old-style Bethe-Bloch parameters
  //
  TF1* parametrisation = FitBB(arr, isMC, isPPb, useV0s, initialParameters, fitType);
  //
  // (3.) obtain response functions to OADB
  //
  TObjArray* arrResponse = GetResponseFunctions(parametrisation, arr, type, period, pass, system, dedxtype);
  //
  // (4.) write results to file and exit
  //
  if (arrResponse) {
    TFile outputFile(outFile,"RECREATE");
    arrResponse->Write();
    outputFile.Close();
  }
  //
  
  return parametrisation->GetParameters();
}


//________________________________________________________________________
TObjArray * AliTPCcalibResidualPID::GetSeparation(THnSparseF * histPidQA, Int_t kParticle1, Int_t kParticle2) {

  //get THnSparse
  THnSparse * hist = histPidQA;


  Float_t nSigmaPlus1 = 0, nSigmaMinus1 = 0; //sigma of first particle in the TPC
  Float_t nSigmaPlus2 = 0, nSigmaMinus2 = 0; //sigma of second particle in the TPC
  
  if(kParticle1 == 0){ // electron
	nSigmaMinus1   =  -1.999;
    nSigmaPlus1   =  2.999;
  } else if(kParticle1 == 1){ //pion
	nSigmaMinus1   =  -1.999;
    nSigmaPlus1   =  2.999;
  } else if(kParticle1 == 2){ //kaons
	nSigmaMinus1   =  -2.999;
    nSigmaPlus1   =    2.999;
  } else if(kParticle1 == 3){ //protons
	nSigmaMinus1   =  -2.999;
    nSigmaPlus1   =    2.999;
  }
  
  //
  // 1. select and fit first particle
  //
  hist->GetAxis(3)->SetRangeUser(kParticle1,kParticle1);
  hist->GetAxis(4)->SetRangeUser(nSigmaMinus1,nSigmaPlus1); 
  hist->GetAxis(5)->SetRangeUser(-2.999,2.999); // 3sigma in TOF
  TH2D * histPart1 = (TH2D*) hist->Projection(1,0)->Clone("histPart1");
  histPart1->SetTitle(";p /(GeV/c) ; TPCSignal /(a.u.)");
  histPart1->RebinX(4);

  TObjArray arr;
  histPart1->FitSlicesY(0,0,-1,10,"QNR",&arr);
  TH1D * PartMean1 = (TH1D *) arr.At(1)->Clone("PartMean1");
  TH1D * PartSigma1 = (TH1D *) arr.At(2)->Clone("PartSigma1");
  PartMean1->SetTitle(";p /(GeV/c) ; Mean (Gauss)");
  PartSigma1->SetTitle(";p /(GeV/c) ; Sigma (Gauss)");

  histPart1->SetMarkerStyle(22);
  PartMean1->SetMarkerStyle(21);
  hist->GetAxis(4)->SetRange(0,-1); // RESET RANGES
  hist->GetAxis(5)->SetRange(0,-1); // RESET RANGES
  
  if(kParticle2==0){ // electron
	nSigmaMinus2   =  -1.999;
    nSigmaPlus2   =    2.999;
  } else if(kParticle2 == 1){ //pion
	nSigmaMinus2   =  -1.999;
    nSigmaPlus2   =    2.999;
  } else if(kParticle2 == 2){ //kaons
	nSigmaMinus2   =  -2.999;
    nSigmaPlus2   =    2.999;
  } else if(kParticle2 == 3){ //protons
	nSigmaMinus2   =  -2.999;
    nSigmaPlus2   =    2.999;
  }

  //
  // 2. select and fit second particle
  //
  hist->GetAxis(3)->SetRangeUser(kParticle2,kParticle2);
  hist->GetAxis(4)->SetRangeUser(nSigmaMinus2,nSigmaPlus2); 
  hist->GetAxis(5)->SetRangeUser(-2.999,2.999); // 3sigma in TOF
  TH2D * histPart2 = (TH2D*)hist->Projection(1,0)->Clone("histPart2");
  histPart2->RebinX(4);
  histPart2->FitSlicesY(0,0,-1,10,"QNR",&arr);
  TH1D * PartMean2 =  (TH1D *) arr.At(1)->Clone("PartMean2");
  TH1D * PartSigma2 =  (TH1D *) arr.At(2)->Clone("PartSigma2");
  PartMean2->SetTitle(";p /(GeV/c) ; Mean (Gauss)");
  PartSigma2->SetTitle(";p /(GeV/c) ; Sigma (Gauss)");
  PartMean2->SetMarkerStyle(20);
  hist->GetAxis(4)->SetRange(0,-1); // RESET RANGES
  hist->GetAxis(5)->SetRange(0,-1); // RESET RANGES

  //
  // 3. separation
  //

	TH1F *fHistSeparation=(TH1F*)PartMean1->Clone("fHistSeparation"); //to get same binning
        fHistSeparation->SetMarkerStyle(22);
	const Int_t Nbins = PartMean1->GetNbinsX();
	fHistSeparation->SetTitle(";p /(GeV/c) ; Separation");

	Float_t DeltaMean[Nbins] ;
	Float_t DeltaSigma[Nbins];

	for(Int_t i=0 ; i<Nbins; i++){
	
	DeltaMean[i] = TMath::Abs(PartMean1->GetBinContent(i) - PartMean2->GetBinContent(i));
	DeltaSigma[i] = TMath::Abs((PartSigma1->GetBinContent(i) + PartSigma2->GetBinContent(i)))/2.;

 	if(!(TMath::Abs(DeltaSigma[i])<0.000001))fHistSeparation->SetBinContent(i,DeltaMean[i]/DeltaSigma[i]);
	}//for(Int_t i=0 ; i<Nbins ; i++)

	TObjArray *array = new TObjArray();
	array->Add(histPart1);
	array->Add(histPart2);
	array->Add(PartMean1);
	array->Add(PartSigma1);
	array->Add(PartMean2);
	array->Add(PartSigma2);
 	array->Add(fHistSeparation);
	return array;

}


//________________________________________________________________________
TObjArray * AliTPCcalibResidualPID::GetResidualGraphs(THnSparseF * histPidQA, const Char_t * system, Bool_t useV0s) {
  //
  // Extracts the residual graphs from THnSparse created from data (NOT MC!)
  //
  
  const TString momTitle = "#it{p}_{TPC} (GeV/#it{c})";
  const TString dEdxTitle = "d#it{E}/d#it{x} (arb. unit)";
  
  Bool_t isPPb = kFALSE;
  if (strcmp(system, "PPB") == 0 || strcmp(system, "PBP") == 0) {
    Printf("p-Pb/Pb-p detected - Applying special handling!");
    isPPb = kTRUE;
  }
  
  // the following parameters have to be extracted from the data itself
  //
    
  Float_t nSigmaPionMinus   =  -3.999;
  Float_t nSigmaPionPlus    =   3.999;
  Float_t nSigmaKaonMinus   =  -4.199;
  Float_t nSigmaKaonPlus    =  4.7999;
  Float_t nSigmaElectronMinus = -1.999;
  Float_t nSigmaElectronPlus  =  4.999;
  
  Int_t cutForFitting = 10;
  Double_t heightFractionForFittingRange = 0.1;
  //
  THnSparse * hist = histPidQA;
  //
  TCanvas * canvasQAtpc = new TCanvas("canvasQAtpcResGraph","Control canvas for residual graphs (TPC)",100,10,1380,800);
  canvasQAtpc->Divide(2,2);
  
  TCanvas * canvasQAtof = new TCanvas("canvasQAtofResGraph","Control canvas for residual graphs (TOF)",100,10,1380,800);
  canvasQAtof->Divide(2,2);
  
  TCanvas * canvasQAv0 = new TCanvas("canvasQAv0ResGraph","Control canvas for residual graphs (V0)",100,10,1380,800);
  canvasQAv0->Divide(2,2);
  
  TCanvas * canvasQAv0plusTOF = new TCanvas("canvasQAv0plusTOFResGraph","Control canvas for residual graphs (V0+TOF)",100,10,1380,800);
  canvasQAv0plusTOF->Divide(2,2);
  
  
  TCanvas * canvasQAv0DeDxPurityEl = new TCanvas("canvasQAv0DeDxPurityEl","Control canvas for residual graphs (V0 el)",100,10,600,400);
  TCanvas * canvasQAv0DeDxPurityPi = new TCanvas("canvasQAv0DeDxPurityPi","Control canvas for residual graphs (V0 pi)",100,10,600,400);
  TCanvas * canvasQAv0DeDxPurityPr = new TCanvas("canvasQAv0DeDxPurityPr","Control canvas for residual graphs (V0 pr)",100,10,600,400);
  
  canvasQAv0DeDxPurityEl->SetLogx();
  canvasQAv0DeDxPurityPi->SetLogx();
  canvasQAv0DeDxPurityPr->SetLogx();
  
  canvasQAv0DeDxPurityEl->SetLogz();
  canvasQAv0DeDxPurityPi->SetLogz();
  canvasQAv0DeDxPurityPr->SetLogz();
  
  canvasQAv0DeDxPurityEl->SetTopMargin(0.03);
  canvasQAv0DeDxPurityPi->SetTopMargin(0.03);
  canvasQAv0DeDxPurityPr->SetTopMargin(0.03);
  
  for (Int_t i = 1; i <= 4; i++) {
    canvasQAtpc->GetPad(i)->SetGrid(1, 1);
    canvasQAtof->GetPad(i)->SetGrid(1, 1);
    canvasQAv0->GetPad(i)->SetGrid(1, 1);
    canvasQAv0plusTOF->GetPad(i)->SetGrid(1, 1);

    canvasQAtpc->GetPad(i)->SetLogz();
    canvasQAtof->GetPad(i)->SetLogz();
    canvasQAv0->GetPad(i)->SetLogz();
    canvasQAv0plusTOF->GetPad(i)->SetLogz();

    canvasQAtpc->GetPad(i)->SetLogx();
    canvasQAtof->GetPad(i)->SetLogx();
    canvasQAv0->GetPad(i)->SetLogx();
    canvasQAv0plusTOF->GetPad(i)->SetLogx();
  }
  //
  // 1. select and fit electrons
  //
  hist->GetAxis(2)->SetRangeUser(0,0); // Non-V0s
  hist->GetAxis(3)->SetRangeUser(0,0); // electrons
  hist->GetAxis(4)->SetRangeUser(nSigmaElectronMinus,nSigmaElectronPlus); // 3sigma in TPC 
  hist->GetAxis(5)->SetRangeUser(-2.999,2.999); // 3sigma in TOF
  TH2D * histElectron = hist->Projection(1,0);
  histElectron->SetName("histElectron");
  histElectron->RebinX(4);
  histElectron->GetXaxis()->SetRangeUser(0.2,3.0);
  TObjArray arr;
  FitSlicesY(histElectron, heightFractionForFittingRange, cutForFitting, "QNR", &arr);
  TH1D * electronPoints = (TH1D *) arr.At(1)->Clone();
  //
  TGraphErrors * electronGraph = new TGraphErrors(electronPoints);
  electronGraph->SetName("electronGraph");
  for(Int_t ip = 0; ip < electronGraph->GetN(); ip ++) {
    Bool_t removePoint = electronGraph->GetY()[ip] < 10 || 
      electronGraph->GetEY()[ip]/electronGraph->GetY()[ip] > 0.05 || 
      electronGraph->GetX()[ip] > 2.0 || 
      electronGraph->GetX()[ip] < 0.5; 
    if (removePoint) {
      electronGraph->RemovePoint(ip);
      ip--;
      continue;
    }
  }
  //
  canvasQAtof->cd(1);
  histElectron->SetTitle("electrons");
  histElectron->GetYaxis()->SetRangeUser(50, 120);
  histElectron->GetXaxis()->SetTitle(momTitle.Data());
  histElectron->GetYaxis()->SetTitle(dEdxTitle.Data());
  histElectron->GetXaxis()->SetMoreLogLabels(kTRUE);
  histElectron->GetXaxis()->SetNoExponent(kTRUE);
  histElectron->Draw("colz");
  electronPoints->SetMarkerStyle(24);
  electronPoints->Draw("same");
  hist->GetAxis(4)->SetRange(0,-1); // RESET RANGES
  hist->GetAxis(5)->SetRange(0,-1); // RESET RANGES
  //
  // 2. protons
  //
  hist->GetAxis(2)->SetRangeUser(0,0); // Non-V0s
  hist->GetAxis(3)->SetRangeUser(3,3); // protons
  //
  //hist->GetAxis(4)->SetRangeUser(nSigmaProtonMinus,nSigmaProtonPlus); // protons --> not reliable, use PATTERN RECOGNITION
  TH2D * histProtonTPC = hist->Projection(1,0);
  histProtonTPC->SetName("histProtonTPC");
  histProtonTPC->GetXaxis()->SetRangeUser(0.15,0.7);
  histProtonTPC->GetYaxis()->SetRangeUser(50, hist->GetAxis(1)->GetBinUpEdge(hist->GetAxis(1)->GetNbins()));
  
  // PATTERN RECOGNITION
  //OLD TF1 betaSq("betaSq","50./TMath::Power(x,1.3)",0.1,10); 
  // Add some additional Erf - cuts away kaons for pPb and PbPb and does not hurt in pp
  TF1 betaSq("betaSq","50./TMath::Power(x,1.3) + (1-TMath::Erf((x-0.2) / 0.075)) * 100",0.1,10);
  for(Int_t ix = 1; ix <= histProtonTPC->GetXaxis()->GetNbins(); ix++) {
     for(Int_t jy = 1; jy <= histProtonTPC->GetYaxis()->GetNbins(); jy++) {
       Float_t yPos = histProtonTPC->GetYaxis()->GetBinCenter(jy);
       Float_t xPos = histProtonTPC->GetXaxis()->GetBinCenter(ix);
       Int_t bin = histProtonTPC->GetBin(ix,jy);
       if (yPos < betaSq.Eval(xPos)) histProtonTPC->SetBinContent(bin,0);
     }
  }
  //
  //TODO Use likelihood option -> Inside fitting range ~no contamination, but little statistics at low momenta requires L
  FitSlicesY(histProtonTPC, heightFractionForFittingRange, cutForFitting, "QNRL", &arr);
  TH1D * protonPointsTPC = (TH1D *) arr.At(1)->Clone();
  //
  TGraphErrors * protonTpcGraph = new TGraphErrors(protonPointsTPC);
  protonTpcGraph->SetName("protonTpcGraph");
  for(Int_t ip = 0; ip < protonTpcGraph->GetN(); ip ++) {
    Bool_t removePoint = protonTpcGraph->GetY()[ip] < 10 || protonTpcGraph->GetEY()[ip]/protonTpcGraph->GetY()[ip] > 0.05 // TODO Larger tolerance, since this is only for the correction function, where 5% are still better than no data point
                        || protonTpcGraph->GetX()[ip] > (useV0s ? (isPPb ? 0.499 : 0.349) : 0.649); // If V0s are to be used, don't use TPC only protons to too high momenta; for pPb low statistics for the moment, therefore, go a bit further with TPC protons
                        //TODO NOW -> Changed range of TPC-signal -> Hopefully, sufficient statistics now for very low momenta!|| protonTpcGraph->GetX()[ip] < 0.19; // Added this cut because almost always bad statistics below 0.19
    if (removePoint) {
      protonTpcGraph->RemovePoint(ip);
      ip--;
      continue;
    }
  }
  //
  hist->GetAxis(5)->SetRangeUser(-2.999,2.999); // 3sigma in TOF
  TH2D * histProtonTOF = hist->Projection(1,0);
  histProtonTOF->SetName("histProtonTOF");
  histProtonTOF->GetXaxis()->SetRangeUser(0.6,2.5);
  
  // Same pattern recognition as before
  for(Int_t ix = 1; ix <= histProtonTOF->GetXaxis()->GetNbins(); ix++) {
    for(Int_t jy = 1; jy <= histProtonTOF->GetYaxis()->GetNbins(); jy++) {
      Float_t yPos = histProtonTOF->GetYaxis()->GetBinCenter(jy);
      Float_t xPos = histProtonTOF->GetXaxis()->GetBinCenter(ix);
      Int_t bin = histProtonTOF->GetBin(ix,jy);
      if (yPos < betaSq.Eval(xPos)) histProtonTOF->SetBinContent(bin,0);
    }
  }
  
  FitSlicesY(histProtonTOF, heightFractionForFittingRange, cutForFitting, "QNR", &arr);

  TH1D * protonPointsTOF =  (TH1D *)arr.At(1)->Clone();
  //
  TGraphErrors * protonTofGraph = new TGraphErrors(protonPointsTOF);
  protonTofGraph->SetName("protonTofGraph");
  for(Int_t ip = 0; ip < protonTofGraph->GetN(); ip ++) {
    // If V0s are to be used, TPC+TOF protons are only for reference and can be plotted to higher momenta
    Bool_t removePoint = protonTofGraph->GetY()[ip] < 10 || protonTofGraph->GetEY()[ip]/protonTofGraph->GetY()[ip] > 0.02 || protonTofGraph->GetX()[ip] > (useV0s ? 3 : 2) || protonTofGraph->GetX()[ip] < 0.65;
    if (removePoint) {
      protonTofGraph->RemovePoint(ip);
      ip--;
      continue;
    }
  }
  //
  canvasQAtpc->cd(2);
  histProtonTPC->SetTitle("protons");
  histProtonTPC->GetXaxis()->SetTitle(momTitle.Data());
  histProtonTPC->GetYaxis()->SetTitle(dEdxTitle.Data());
  histProtonTPC->GetXaxis()->SetMoreLogLabels(kTRUE);
  histProtonTPC->GetXaxis()->SetNoExponent(kTRUE);
  histProtonTPC->Draw("colz");
  betaSq.DrawCopy("same");
  protonPointsTPC->SetMarkerStyle(20);
  protonPointsTOF->SetMarkerStyle(24);
  protonPointsTOF->Draw("same");
  protonPointsTPC->Draw("same");
  //
  protonTpcGraph->SetMarkerStyle(26);
  protonTpcGraph->SetMarkerColor(kMagenta);
  protonTpcGraph->DrawClone("p");
  protonTofGraph->SetMarkerStyle(25);
  protonTofGraph->SetMarkerColor(kMagenta);
  protonTofGraph->DrawClone("p");
  
  canvasQAtof->cd(2);
  histProtonTOF->SetTitle("protons");
  histProtonTOF->GetYaxis()->SetRangeUser(30, 250);
  histProtonTOF->GetXaxis()->SetTitle(momTitle.Data());
  histProtonTOF->GetYaxis()->SetTitle(dEdxTitle.Data());
  histProtonTOF->GetXaxis()->SetMoreLogLabels(kTRUE);
  histProtonTOF->GetXaxis()->SetNoExponent(kTRUE);
  histProtonTOF->Draw("colz");
  betaSq.DrawCopy("same");
  protonPointsTOF->Draw("same");
  protonPointsTPC->Draw("same");
  //
  protonTpcGraph->DrawClone("p");
  protonTofGraph->DrawClone("p");
  
  hist->GetAxis(4)->SetRange(0,-1); // RESET RANGES
  hist->GetAxis(5)->SetRange(0,-1); // RESET RANGES
  //
  // 3. pions
  //
  hist->GetAxis(2)->SetRangeUser(0,0); // Non-V0s
  hist->GetAxis(3)->SetRangeUser(1,1); // pions
  //
  Double_t lowerPionThreshold = 0.3;

  hist->GetAxis(4)->SetRangeUser(nSigmaPionMinus,nSigmaPionPlus); // pions
  TH2D * histPionTPC = hist->Projection(1,0);
  histPionTPC->SetName("histPionTPC");
  histPionTPC->GetXaxis()->SetRangeUser(0.15,0.7);
  FitSlicesY(histPionTPC, heightFractionForFittingRange, cutForFitting, "QNR", &arr);

  // In case of really bad splines comment last 5 lines and use the following 6 lines:
  //TH2D * histPionTPC = hist->Projection(1,0);
  //histPionTPC->SetName("histPionTPC");
  //histPionTPC->GetYaxis()->SetRangeUser(30,100);
  //histPionTPC->GetXaxis()->SetRangeUser(0.15,0.7);
  //FitSlicesY(histPionTPC, heightFractionForFittingRange, cutForFitting, "QNR", &arr);
  //lowerPionThreshold = 0.15;


  TH1D * pionPointsTPC =  (TH1D *) arr.At(1)->Clone();
  //
  TGraphErrors * pionTpcGraph = new TGraphErrors(pionPointsTPC);
  pionTpcGraph->SetName("pionTpcGraph");
  for(Int_t ip = 0; ip < pionTpcGraph->GetN(); ip ++) {
    Bool_t removePoint = pionTpcGraph->GetY()[ip] < 10 || pionTpcGraph->GetEY()[ip]/pionTpcGraph->GetY()[ip] > 0.02 || pionTpcGraph->GetX()[ip] > 0.5
                         || pionTpcGraph->GetX()[ip] < lowerPionThreshold; // Exclude contamination by electrons
    if (removePoint) {
      pionTpcGraph->RemovePoint(ip);
      ip--;
      continue;
    }
  }
  //
  hist->GetAxis(5)->SetRangeUser(-2.999,2.999); // 3sigma in TOF
  TH2D * histPionTOF = hist->Projection(1,0);
  histPionTOF->SetName("histPionTOF");
  histPionTOF->GetXaxis()->SetRangeUser(0.5,1.1);
  FitSlicesY(histPionTOF, heightFractionForFittingRange, cutForFitting, "QNR", &arr);
  TH1D * pionPointsTOF = (TH1D *) arr.At(1)->Clone();
  TGraphErrors * pionTofGraph = new TGraphErrors(pionPointsTOF);
  pionTofGraph->SetName("pionTofGraph");
  for(Int_t ip = 0; ip < pionTofGraph->GetN(); ip ++) {
    Bool_t removePoint = pionTofGraph->GetY()[ip] < 10 || pionTofGraph->GetEY()[ip]/pionTofGraph->GetY()[ip] > 0.02 || pionTofGraph->GetX()[ip] > 1.1 || pionTofGraph->GetX()[ip] < 0.5; // Changed by Ben (was 0.35) to avoid problems with TOF efficiency/systematics
    if (removePoint) {
      pionTofGraph->RemovePoint(ip);
      ip--;
      continue;
    }
  }
  canvasQAtpc->cd(3);
  histPionTPC->GetYaxis()->SetRangeUser(30, 90);
  histPionTPC->SetTitle("pions");
  histPionTPC->GetXaxis()->SetTitle(momTitle.Data());
  histPionTPC->GetYaxis()->SetTitle(dEdxTitle.Data());
  histPionTPC->GetXaxis()->SetMoreLogLabels(kTRUE);
  histPionTPC->GetXaxis()->SetNoExponent(kTRUE);
  histPionTPC->Draw("colz");
  pionPointsTPC->SetMarkerStyle(20);
  pionPointsTOF->SetMarkerStyle(24);
  pionPointsTOF->Draw("same");
  pionPointsTPC->Draw("same");
  //
  pionTpcGraph->SetMarkerStyle(26);
  pionTpcGraph->SetMarkerColor(kMagenta);
  pionTpcGraph->DrawClone("p");
  pionTofGraph->SetMarkerStyle(25);
  pionTofGraph->SetMarkerColor(kMagenta);
  pionTofGraph->DrawClone("p");
  
  canvasQAtof->cd(3);
  histPionTOF->GetYaxis()->SetRangeUser(30, 80);
  histPionTOF->GetXaxis()->SetTitle(momTitle.Data());
  histPionTOF->GetYaxis()->SetTitle(dEdxTitle.Data());
  histPionTOF->GetXaxis()->SetMoreLogLabels(kTRUE);
  histPionTOF->GetXaxis()->SetNoExponent(kTRUE);
  histPionTOF->Draw("colz");
  histPionTOF->SetTitle("pions");
  pionPointsTOF->Draw("same");
  pionPointsTPC->Draw("same");
  //
  pionTpcGraph->DrawClone("p");
  pionTofGraph->DrawClone("p");

  
  hist->GetAxis(4)->SetRange(0,-1); // RESET RANGES
  hist->GetAxis(5)->SetRange(0,-1); // RESET RANGES
  //
  // 4. kaons
  //
  hist->GetAxis(2)->SetRangeUser(0,0); // Non-V0s
  hist->GetAxis(3)->SetRangeUser(2,2); // kaons
  //
  hist->GetAxis(4)->SetRangeUser(nSigmaKaonMinus,nSigmaKaonPlus); // kaons
  TH2D * histKaonTPC = hist->Projection(1,0);
  histKaonTPC->SetName("histKaonTPC");
  histKaonTPC->GetXaxis()->SetRangeUser(0.12,0.6999);
  FitSlicesY(histKaonTPC, heightFractionForFittingRange, cutForFitting, "QNR", &arr);
  TH1D * kaonPointsTPC = (TH1D*) arr.At(1)->Clone();
  //
  TGraphErrors * kaonTpcGraph = new TGraphErrors(kaonPointsTPC);
  kaonTpcGraph->SetName("kaonTpcGraph");
  for(Int_t ip = 0; ip < kaonTpcGraph->GetN(); ip ++) {
    Bool_t removePoint = kaonTpcGraph->GetY()[ip] < 10 || kaonTpcGraph->GetEY()[ip]/kaonTpcGraph->GetY()[ip] > 0.02
                         || kaonTpcGraph->GetX()[ip] < 0.12 || kaonTpcGraph->GetX()[ip] > 0.319; //TODO BEN was > 0.399
    if (removePoint) {
      kaonTpcGraph->RemovePoint(ip);
      ip--;
      continue;
    }
  }
  //
  
  hist->GetAxis(5)->SetRangeUser(-2.999,2.999); // sigma in TOF
  
  /*
  hist->GetAxis(5)->SetRangeUser(-1.999,1.999); // sigma in TOF
  
  if (hist->GetAxis(3)->GetNbins() >= 5) {
    printf("Using restricted eta range for TPC+TOF Kaons!\n");
    hist->GetAxis(3)->SetRangeUser(4,4);
  }
  else {
    printf("WARNING: Data points for restricted eta range for TPC+TOF Kaons not available! Using full eta range (worse w.r.t. contamination)\n");
  }*/
  
  TH2D * histKaonTOF = hist->Projection(1,0);
  histKaonTOF->SetName("histKaonTOF");
  histKaonTOF->GetXaxis()->SetRangeUser(0.3,1.5);
  FitSlicesY(histKaonTOF, heightFractionForFittingRange, cutForFitting, "QNR", &arr);
  TH1D * kaonPointsTOF = (TH1D*) arr.At(1)->Clone();
  //
  TGraphErrors * kaonTofGraph = new TGraphErrors(kaonPointsTOF);
  kaonTofGraph->SetName("kaonTofGraph");
  //
  for(Int_t ip = 0; ip < kaonTofGraph->GetN(); ip ++) {
    Bool_t removePoint = kaonTofGraph->GetY()[ip] < 10 || kaonTofGraph->GetEY()[ip]/kaonTofGraph->GetY()[ip] > 0.02 || kaonTofGraph->GetX()[ip] > 1.0
                         || kaonTofGraph->GetX()[ip] < 0.36; //TODO BEN NOW was 0.5
    if (removePoint) {
      kaonTofGraph->RemovePoint(ip);
      ip--;
      continue;
    }
  }
  //
  canvasQAtpc->cd(4);
  histKaonTPC->SetTitle("kaons");
  histKaonTPC->GetYaxis()->SetRangeUser(30, 500);
  histKaonTPC->GetXaxis()->SetTitle(momTitle.Data());
  histKaonTPC->GetYaxis()->SetTitle(dEdxTitle.Data());
  histKaonTPC->GetXaxis()->SetMoreLogLabels(kTRUE);
  histKaonTPC->GetXaxis()->SetNoExponent(kTRUE);
  histKaonTPC->Draw("colz");
  kaonPointsTPC->SetMarkerStyle(20);
  kaonPointsTOF->SetMarkerStyle(24);
  kaonPointsTOF->Draw("same");
  kaonPointsTPC->Draw("same");
  //
  kaonTpcGraph->SetMarkerStyle(26);
  kaonTpcGraph->SetMarkerColor(kMagenta);
  kaonTpcGraph->DrawClone("p");
  kaonTofGraph->SetMarkerStyle(25);
  kaonTofGraph->SetMarkerColor(kMagenta);
  kaonTofGraph->DrawClone("p");

  
  canvasQAtof->cd(4);
  histKaonTOF->GetYaxis()->SetRangeUser(30, 250);
  histKaonTOF->SetTitle("kaons");
  histKaonTOF->GetXaxis()->SetTitle(momTitle.Data());
  histKaonTOF->GetYaxis()->SetTitle(dEdxTitle.Data());
  histKaonTOF->GetXaxis()->SetMoreLogLabels(kTRUE);
  histKaonTOF->GetXaxis()->SetNoExponent(kTRUE);
  histKaonTOF->Draw("colz");
  kaonPointsTOF->Draw("same");
  kaonPointsTPC->Draw("same");
  kaonTpcGraph->DrawClone("p");
  kaonTofGraph->DrawClone("p");
  
  hist->GetAxis(4)->SetRange(0,-1); // RESET RANGES
  hist->GetAxis(5)->SetRange(0,-1); // RESET RANGES
  
  
  //
  // 5. V0 electrons
  //
  hist->GetAxis(2)->SetRangeUser(1,1); // V0 electrons
  hist->GetAxis(3)->SetRangeUser(0,0); // electrons
  
  //
  TH2D * histElectronV0 = hist->Projection(1,0);
  histElectronV0->SetName("histElectronV0");
  //histElectronV0->RebinX(2);
  histElectronV0->GetXaxis()->SetRangeUser(0.15,5.0);
  FitSlicesY(histElectronV0, heightFractionForFittingRange, cutForFitting, "QNR", &arr);
  TH1D * electronPointsV0 = (TH1D*) arr.At(1)->Clone();
  //
  TGraphErrors * electronV0Graph = new TGraphErrors(electronPointsV0);
  electronV0Graph->SetName("electronV0Graph");
  for(Int_t ip = 0; ip < electronV0Graph->GetN(); ip ++) {
    Bool_t removePoint = electronV0Graph->GetY()[ip] < 10 || electronV0Graph->GetEY()[ip]/electronV0Graph->GetY()[ip] > 0.02;

    if (removePoint) {
      electronV0Graph->RemovePoint(ip);
      ip--;
      continue;
    }
  }
  //
  canvasQAv0->cd(1);
  histElectronV0->SetTitle("V0 electrons");
  histElectronV0->GetYaxis()->SetRangeUser(50, 120);
  histElectronV0->GetXaxis()->SetTitle(momTitle.Data());
  histElectronV0->GetYaxis()->SetTitle(dEdxTitle.Data());
  histElectronV0->GetXaxis()->SetMoreLogLabels(kTRUE);
  histElectronV0->GetXaxis()->SetNoExponent(kTRUE);
  histElectronV0->SetStats(0);
  histElectronV0->Draw("colz");
  electronPointsV0->SetMarkerStyle(34);
  electronPointsV0->Draw("same");
  electronGraph->SetMarkerStyle(25);
  electronGraph->SetMarkerColor(kMagenta);
  electronGraph->DrawClone("p");
  
  canvasQAv0DeDxPurityEl->cd();
  gStyle->SetOptTitle(0);
  TH1D* hElNew = (TH1D*)histElectronV0->DrawClone("colz");
  hElNew->GetXaxis()->SetTitleOffset(1.0);
  hElNew->SetTitle("");
  electronPointsV0->DrawClone("same");
  gStyle->SetOptTitle(1);



  canvasQAtof->cd(1);
  electronV0Graph->SetMarkerStyle(24);
  electronV0Graph->SetMarkerColor(kMagenta);
  electronV0Graph->DrawClone("p");
  
  //
  // 6. V0 pions
  //
  hist->GetAxis(2)->SetRangeUser(2,2); // V0 pions
  hist->GetAxis(3)->SetRangeUser(1,1); // pions
  //
  TH2D * histPionV0 = hist->Projection(1,0);
  histPionV0->SetName("histPionV0");
  histPionV0->GetXaxis()->SetRangeUser(0.15,10.);
  FitSlicesY(histPionV0, heightFractionForFittingRange, cutForFitting, "QNR", &arr);
  TH1D * pionPointsV0 = (TH1D*) arr.At(1)->Clone();
  //
  TGraphErrors * pionV0Graph = new TGraphErrors(pionPointsV0);
  pionV0Graph->SetName("pionV0Graph");
  for(Int_t ip = 0; ip < pionV0Graph->GetN(); ip ++) {
    Bool_t removePoint = pionV0Graph->GetY()[ip] < 10 || pionV0Graph->GetEY()[ip]/pionV0Graph->GetY()[ip] > 0.02;

    if (removePoint) {
      pionV0Graph->RemovePoint(ip);
      ip--;
      continue;
    }
  }
  //
  canvasQAv0->cd(3);
  histPionV0->SetTitle("V0 pions");
  histPionV0->GetYaxis()->SetRangeUser(30, 100);
  histPionV0->GetXaxis()->SetTitle(momTitle.Data());
  histPionV0->GetYaxis()->SetTitle(dEdxTitle.Data());
  histPionV0->GetXaxis()->SetMoreLogLabels(kTRUE);
  histPionV0->GetXaxis()->SetNoExponent(kTRUE);
  histPionV0->Draw("colz");
  pionPointsV0->SetMarkerStyle(34);
  pionPointsV0->Draw("same");
  pionTofGraph->DrawClone("p");
  pionTpcGraph->DrawClone("p");
  
  canvasQAv0DeDxPurityPi->cd();
  gStyle->SetOptTitle(0);
  TH1D* hPiNew = (TH1D*)histPionV0->DrawClone("colz");
  hPiNew->GetXaxis()->SetTitleOffset(1.0);
  hPiNew->SetTitle("");
  pionPointsV0->DrawClone("same");
  gStyle->SetOptTitle(1);
  
  canvasQAtof->cd(3);
  pionV0Graph->SetMarkerStyle(24);
  pionV0Graph->SetMarkerColor(kMagenta);
  pionV0Graph->DrawClone("p");

  canvasQAtpc->cd(3);
  pionV0Graph->DrawClone("p");
  
  //
  // 6. V0 protons
  //
  hist->GetAxis(2)->SetRangeUser(3,3); // V0 protons
  hist->GetAxis(3)->SetRangeUser(3,3); // protons
  //
  TH2D * histProtonV0 = hist->Projection(1,0);
  histProtonV0->SetName("histProtonV0");
  histProtonV0->GetXaxis()->SetRangeUser(0.15,10.);
  
  // Remove contamination due to el and pions and low momenta because there the statistics
  // for the V0 protons goes down and becomes comparable with the contamination
  // -> This spoils the fit, if the contamination is not removed
  
  // PATTERN RECOGNITION
  Int_t correctionUpToBin = histProtonV0->GetXaxis()->FindBin(0.6);
  
  for(Int_t ix = 1; ix <= correctionUpToBin; ix++) {
     for(Int_t jy = 1; jy <= histProtonV0->GetYaxis()->GetNbins(); jy++) {
       Float_t yPos = histProtonV0->GetYaxis()->GetBinCenter(jy);
       Float_t xPos = histProtonV0->GetXaxis()->GetBinCenter(ix);
       Int_t bin = histProtonV0->GetBin(ix,jy);
       if (yPos < betaSq.Eval(xPos)) histProtonV0->SetBinContent(bin,0);
     }
  }
  
  /*
  Int_t correctionUpToBin = histProtonV0->GetXaxis()->FindBin(0.4);
  Double_t threshold = 120.;
  
  for(Int_t ix = 1; ix <= correctionUpToBin; ix++) {
     for(Int_t jy = 1; jy <= histProtonV0->GetYaxis()->GetNbins(); jy++) {
       Float_t yPos = histProtonV0->GetYaxis()->GetBinCenter(jy);
       Int_t bin = histProtonV0->GetBin(ix,jy);
       if (yPos < threshold) histProtonV0->SetBinContent(bin,0);
       else break;
     }
  }
  */
  
  FitSlicesY(histProtonV0, heightFractionForFittingRange, cutForFitting, "QNR", &arr);
  TH1D * protonPointsV0 = (TH1D*) arr.At(1)->Clone();
  //
  TGraphErrors * protonV0Graph = new TGraphErrors(protonPointsV0);
  protonV0Graph->SetName("protonV0Graph");
  for(Int_t ip = 0; ip < protonV0Graph->GetN(); ip ++) {
    Bool_t removePoint = protonV0Graph->GetY()[ip] < 10 || protonV0Graph->GetEY()[ip]/protonV0Graph->GetY()[ip] > 0.02;
                         
    if (removePoint) {
      protonV0Graph->RemovePoint(ip);
      ip--;
      continue;
    }
  }
  //
  canvasQAv0->cd(2);
  histProtonV0->SetTitle("V0 protons");
  histProtonV0->GetXaxis()->SetTitle(momTitle.Data());
  histProtonV0->GetYaxis()->SetTitle(dEdxTitle.Data());
  histProtonV0->GetXaxis()->SetMoreLogLabels(kTRUE);
  histProtonV0->GetXaxis()->SetNoExponent(kTRUE);
  histProtonV0->Draw("colz");
  protonPointsV0->SetMarkerStyle(34);
  protonPointsV0->Draw("same");
  protonTofGraph->DrawClone("p");
  protonTpcGraph->DrawClone("p");
  
  canvasQAv0DeDxPurityPr->cd();
  gStyle->SetOptTitle(0);
  TH1D* hPrNew = (TH1D*)histProtonV0->DrawClone("colz");
  hPrNew->GetXaxis()->SetTitleOffset(1.0);
  hPrNew->SetTitle("");
  protonPointsV0->DrawClone("same");
  gStyle->SetOptTitle(1);

  canvasQAtof->cd(2);
  protonV0Graph->SetMarkerStyle(24);
  protonV0Graph->SetMarkerColor(kMagenta);
  protonV0Graph->DrawClone("p");

  canvasQAtpc->cd(2);
  protonV0Graph->DrawClone("p");
  
  
  hist->GetAxis(2)->SetRange(0,-1); // RESET RANGES
  hist->GetAxis(3)->SetRange(0,-1); // RESET RANGES
  
  
  //
  // 5. V0 electrons + TOF
  //
  hist->GetAxis(2)->SetRangeUser(1,1); // V0 electrons
  hist->GetAxis(3)->SetRangeUser(0,0); // electrons
  hist->GetAxis(5)->SetRangeUser(-2.999,2.999); // 3sigma in TOF
  //
  TH2D * histElectronV0plusTOF = hist->Projection(1,0);
  histElectronV0plusTOF->SetName("histElectronV0plusTOF");
  histElectronV0plusTOF->RebinX(2);
  histElectronV0plusTOF->GetXaxis()->SetRangeUser(0.2,5.0);
  FitSlicesY(histElectronV0plusTOF, heightFractionForFittingRange, cutForFitting, "QNR", &arr);
  TH1D * electronPointsV0plusTOF = (TH1D*) arr.At(1)->Clone();
  //
  TGraphErrors * electronV0plusTOFGraph = new TGraphErrors(electronPointsV0plusTOF);
  electronV0plusTOFGraph->SetName("electronV0plusTOFGraph");
  for(Int_t ip = 0; ip < electronV0plusTOFGraph->GetN(); ip ++) {
    Bool_t removePoint = electronV0plusTOFGraph->GetY()[ip] < 10 || electronV0plusTOFGraph->GetEY()[ip]/electronV0plusTOFGraph->GetY()[ip] > 0.02;

    if (removePoint) {
      electronV0plusTOFGraph->RemovePoint(ip);
      ip--;
      continue;
    }
  }
  //
  canvasQAv0plusTOF->cd(1);
  histElectronV0plusTOF->SetTitle("V0+TOF electrons");
  histElectronV0plusTOF->GetYaxis()->SetRangeUser(50, 120);
  histElectronV0plusTOF->GetXaxis()->SetTitle(momTitle.Data());
  histElectronV0plusTOF->GetYaxis()->SetTitle(dEdxTitle.Data());
  histElectronV0plusTOF->GetXaxis()->SetMoreLogLabels(kTRUE);
  histElectronV0plusTOF->GetXaxis()->SetNoExponent(kTRUE);
  histElectronV0plusTOF->Draw("colz");
  electronPointsV0plusTOF->SetMarkerStyle(29);
  electronPointsV0plusTOF->Draw("same");
  electronGraph->DrawClone("p");
  electronV0Graph->DrawClone("p");

  canvasQAv0->cd(1);
  electronV0plusTOFGraph->SetMarkerStyle(30);
  electronV0plusTOFGraph->SetMarkerColor(kMagenta);
  electronV0plusTOFGraph->DrawClone("p");

  canvasQAtof->cd(1);
  electronV0plusTOFGraph->DrawClone("p");
  
  //
  // 6. V0 pions + TOF
  //
  hist->GetAxis(2)->SetRangeUser(2,2); // V0 pions
  hist->GetAxis(3)->SetRangeUser(1,1); // pions
  hist->GetAxis(5)->SetRangeUser(-2.999,2.999); // 3sigma in TOF
  //
  TH2D * histPionV0plusTOF = hist->Projection(1,0);
  histPionV0plusTOF->SetName("histPionV0plusTOF");
  histPionV0plusTOF->GetXaxis()->SetRangeUser(0.15,10.);
  FitSlicesY(histPionV0plusTOF, heightFractionForFittingRange, cutForFitting, "QNR", &arr);
  TH1D * pionPointsV0plusTOF = (TH1D*) arr.At(1)->Clone();
  //
  TGraphErrors * pionV0plusTOFGraph = new TGraphErrors(pionPointsV0plusTOF);
  pionV0plusTOFGraph->SetName("pionV0plusTOFGraph");
  for(Int_t ip = 0; ip < pionV0plusTOFGraph->GetN(); ip ++) {
    Bool_t removePoint = pionV0plusTOFGraph->GetY()[ip] < 10 || pionV0plusTOFGraph->GetEY()[ip]/pionV0plusTOFGraph->GetY()[ip] > 0.02;

    if (removePoint) {
      pionV0plusTOFGraph->RemovePoint(ip);
      ip--;
      continue;
    }
  }
  //
  canvasQAv0plusTOF->cd(3);
  histPionV0plusTOF->SetTitle("V0+TOF pions");
  histPionV0plusTOF->GetYaxis()->SetRangeUser(30, 100);
  histPionV0plusTOF->GetXaxis()->SetTitle(momTitle.Data());
  histPionV0plusTOF->GetYaxis()->SetTitle(dEdxTitle.Data());
  histPionV0plusTOF->GetXaxis()->SetMoreLogLabels(kTRUE);
  histPionV0plusTOF->GetXaxis()->SetNoExponent(kTRUE);
  histPionV0plusTOF->Draw("colz");
  pionPointsV0plusTOF->SetMarkerStyle(29);
  pionPointsV0plusTOF->Draw("same");
  pionV0Graph->DrawClone("p");
  pionTofGraph->DrawClone("p");
  pionTpcGraph->DrawClone("p");
  
  canvasQAv0->cd(3);
  pionV0plusTOFGraph->SetMarkerStyle(30);
  pionV0plusTOFGraph->SetMarkerColor(kMagenta);
  pionV0plusTOFGraph->DrawClone("p");
  
  canvasQAtof->cd(3);
  pionV0plusTOFGraph->DrawClone("p");
  
  canvasQAtpc->cd(3);
  pionV0plusTOFGraph->DrawClone("p");
  
  //
  // 6. V0 protons + TOF
  //
  hist->GetAxis(2)->SetRangeUser(3,3); // V0 protons
  hist->GetAxis(3)->SetRangeUser(3,3); // protons
  hist->GetAxis(5)->SetRangeUser(-2.999,2.999); // 3sigma in TOF
  //
  TH2D * histProtonV0plusTOF = hist->Projection(1,0);
  histProtonV0plusTOF->SetName("histProtonV0plusTOF");
  histProtonV0plusTOF->GetXaxis()->SetRangeUser(0.15,10.);
  
  // Remove contamination due to el and pions and low momenta because there the statistics
  // for the V0 protons goes down and becomes comparable with the contamination
  // -> This spoils the fit, if the contamination is not removed
  
  // PATTERN RECOGNITION
  correctionUpToBin = histProtonV0plusTOF->GetXaxis()->FindBin(0.6);
  
  for(Int_t ix = 1; ix <= correctionUpToBin; ix++) {
     for(Int_t jy = 1; jy <= histProtonV0plusTOF->GetYaxis()->GetNbins(); jy++) {
       Float_t yPos = histProtonV0plusTOF->GetYaxis()->GetBinCenter(jy);
       Float_t xPos = histProtonV0plusTOF->GetXaxis()->GetBinCenter(ix);
       Int_t bin = histProtonV0plusTOF->GetBin(ix,jy);
       if (yPos < betaSq.Eval(xPos)) histProtonV0plusTOF->SetBinContent(bin,0);
     }
  }
  
  FitSlicesY(histProtonV0plusTOF, heightFractionForFittingRange, cutForFitting, "QNR", &arr);
  TH1D * protonPointsV0plusTOF = (TH1D*) arr.At(1)->Clone();
  //
  TGraphErrors * protonV0plusTOFGraph = new TGraphErrors(protonPointsV0plusTOF);
  protonV0plusTOFGraph->SetName("protonV0plusTOFGraph");
  for(Int_t ip = 0; ip < protonV0plusTOFGraph->GetN(); ip ++) {
    Bool_t removePoint = protonV0plusTOFGraph->GetY()[ip] < 10 || protonV0plusTOFGraph->GetEY()[ip]/protonV0plusTOFGraph->GetY()[ip] > 0.02;
                         
    if (removePoint) {
      protonV0plusTOFGraph->RemovePoint(ip);
      ip--;
      continue;
    }
  }
  //
  canvasQAv0plusTOF->cd(2);
  histProtonV0plusTOF->SetTitle("V0+TOF protons");
  histProtonV0plusTOF->GetXaxis()->SetTitle(momTitle.Data());
  histProtonV0plusTOF->GetYaxis()->SetTitle(dEdxTitle.Data());
  histProtonV0plusTOF->GetXaxis()->SetMoreLogLabels(kTRUE);
  histProtonV0plusTOF->GetXaxis()->SetNoExponent(kTRUE);
  histProtonV0plusTOF->Draw("colz");
  protonPointsV0plusTOF->SetMarkerStyle(29);
  protonPointsV0plusTOF->Draw("same");
  protonV0Graph->DrawClone("p");
  protonTofGraph->DrawClone("p");
  protonTpcGraph->DrawClone("p");

  canvasQAv0->cd(2);
  protonV0plusTOFGraph->SetMarkerStyle(30);
  protonV0plusTOFGraph->SetMarkerColor(kMagenta);
  protonV0plusTOFGraph->DrawClone("p");
  
  canvasQAtof->cd(2);
  protonV0plusTOFGraph->DrawClone("p");

  canvasQAtpc->cd(2);
  protonV0plusTOFGraph->DrawClone("p");
  
  hist->GetAxis(2)->SetRange(0,-1); // RESET RANGES
  hist->GetAxis(3)->SetRange(0,-1); // RESET RANGES
  hist->GetAxis(5)->SetRange(0,-1); // RESET RANGES
  
  
  // Clone (V0, V0+TOF) graphs: Remove some points, which are expected to deviate from the Bethe-Bloch curve, from the original graphs, so that
  // these graphs can be used for the BB fit. The original graph will keep all point in order to determine the correction and response
  // functions
  TGraphErrors * electronV0GraphForBBfit = new TGraphErrors(*electronV0Graph);
  electronV0GraphForBBfit->SetName("electronV0GraphForBBfit");
  
  for(Int_t ip = 0; ip < electronV0GraphForBBfit->GetN(); ip ++) {
    Bool_t removePoint = electronV0GraphForBBfit->GetX()[ip] < (isPPb ? 1.2 : 0.4);// || electronV0GraphForBBfit->GetX()[ip] > 2.2;
    
    if (removePoint) {
      electronV0GraphForBBfit->RemovePoint(ip);
      ip--;
      continue;
    }
  }
  
  TGraphErrors * pionV0GraphForBBfit = new TGraphErrors(*pionV0Graph);
  pionV0GraphForBBfit->SetName("pionV0GraphForBBfit");
  
  for(Int_t ip = 0; ip < pionV0GraphForBBfit->GetN(); ip ++) {
    Bool_t removePoint = pionV0GraphForBBfit->GetX()[ip] < (isPPb ? 1.0 : 0.76);// || pionV0GraphForBBfit->GetX()[ip] > 6.; //TODO NOW was < 0.4 before
    //Bool_t removePoint = pionV0GraphForBBfit->GetX()[ip] < 0.4 || pionV0GraphForBBfit->GetX()[ip] > 0.399;// In this case NOT used for BB fit
    if (removePoint) {
      pionV0GraphForBBfit->RemovePoint(ip);
      ip--;
      continue;
    }
  }
  
  TGraphErrors * protonV0GraphForBBfit = new TGraphErrors(*protonV0Graph);
  protonV0GraphForBBfit->SetName("protonV0GraphForBBfit");
  
  for(Int_t ip = 0; ip < protonV0GraphForBBfit->GetN(); ip ++) {
    Bool_t removePoint = protonV0GraphForBBfit->GetX()[ip] < 0.9 || protonV0GraphForBBfit->GetX()[ip] > 5.0; //TODO NOW was 2.5 before
    //Bool_t removePoint = protonV0GraphForBBfit->GetX()[ip] < 0.9 || protonV0GraphForBBfit->GetX()[ip] > 0.599;// In this case NOT used for BB fit
    if (removePoint) {
      protonV0GraphForBBfit->RemovePoint(ip);
      ip--;
      continue;
    }
  }
  
  
  // V0 plus TOF
  TGraphErrors * electronV0plusTOFGraphForBBfit = new TGraphErrors(*electronV0plusTOFGraph);
  electronV0plusTOFGraphForBBfit->SetName("electronV0plusTOFGraphForBBfit");
  
  for(Int_t ip = 0; ip < electronV0plusTOFGraphForBBfit->GetN(); ip ++) {
    Bool_t removePoint = electronV0plusTOFGraphForBBfit->GetX()[ip] < (isPPb ? 1.2 : 0.6);//0.4 || electronV0plusTOFGraphForBBfit->GetX()[ip] > 1.3799;
    if (removePoint) {
      electronV0plusTOFGraphForBBfit->RemovePoint(ip);
      ip--;
      continue;
    }
  }
  
  TGraphErrors * pionV0plusTOFGraphForBBfit = new TGraphErrors(*pionV0plusTOFGraph);
  pionV0plusTOFGraphForBBfit->SetName("pionV0plusTOFGraphForBBfit");
  
  for(Int_t ip = 0; ip < pionV0plusTOFGraphForBBfit->GetN(); ip ++) {
    Bool_t removePoint = pionV0plusTOFGraphForBBfit->GetX()[ip] < 0.4;// || pionV0plusTOFGraphForBBfit->GetX()[ip] > 6.;
    if (removePoint) {
      pionV0plusTOFGraphForBBfit->RemovePoint(ip);
      ip--;
      continue;
    }
  }
  
  TGraphErrors * protonV0plusTOFGraphForBBfit = new TGraphErrors(*protonV0plusTOFGraph);
  protonV0plusTOFGraphForBBfit->SetName("protonV0plusTOFGraphForBBfit");
  
  for(Int_t ip = 0; ip < protonV0plusTOFGraphForBBfit->GetN(); ip ++) {
    Bool_t removePoint = protonV0plusTOFGraphForBBfit->GetX()[ip] < 0.9 || protonV0plusTOFGraphForBBfit->GetX()[ip] > 2.5;
    if (removePoint) {
      protonV0plusTOFGraphForBBfit->RemovePoint(ip);
      ip--;
      continue;
    }
  }
  
  
  //
  // rescale graphs to betaGamma
  //
  kaonTofGraph->SetMarkerStyle(21);
  kaonTpcGraph->SetMarkerStyle(22);
  for(Int_t ip = 0; ip < kaonTofGraph->GetN(); ip ++) {
    kaonTofGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kKaon);
    kaonTofGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kKaon);
  }
  for(Int_t ip = 0; ip < kaonTpcGraph->GetN(); ip ++) {
    kaonTpcGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kKaon);
    kaonTpcGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kKaon);
  }
  //
  electronGraph->SetMarkerStyle(22);
  electronV0Graph->SetMarkerStyle(34);
  electronV0GraphForBBfit->SetMarkerStyle(34);
  electronV0plusTOFGraph->SetMarkerStyle(30);
  electronV0plusTOFGraphForBBfit->SetMarkerStyle(30);
  electronGraph->SetMarkerColor(2);
  electronV0Graph->SetMarkerColor(2);
  electronV0GraphForBBfit->SetMarkerColor(2);
  electronV0plusTOFGraph->SetMarkerColor(2);
  electronV0plusTOFGraphForBBfit->SetMarkerColor(2);
  for(Int_t ip = 0; ip < electronGraph->GetN(); ip ++) {
    electronGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kElectron);
    electronGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kElectron);
  }
  for(Int_t ip = 0; ip < electronV0Graph->GetN(); ip ++) {
    electronV0Graph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kElectron);
    electronV0Graph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kElectron);
  }
  for(Int_t ip = 0; ip < electronV0GraphForBBfit->GetN(); ip ++) {
    electronV0GraphForBBfit->GetX()[ip] /= AliPID::ParticleMass(AliPID::kElectron);
    electronV0GraphForBBfit->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kElectron);
  }
  for(Int_t ip = 0; ip < electronV0plusTOFGraph->GetN(); ip ++) {
    electronV0plusTOFGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kElectron);
    electronV0plusTOFGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kElectron);
  }
  for(Int_t ip = 0; ip < electronV0plusTOFGraphForBBfit->GetN(); ip ++) {
    electronV0plusTOFGraphForBBfit->GetX()[ip] /= AliPID::ParticleMass(AliPID::kElectron);
    electronV0plusTOFGraphForBBfit->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kElectron);
  }
  //
  pionTofGraph->SetMarkerStyle(21);
  pionTpcGraph->SetMarkerStyle(22);
  pionV0Graph->SetMarkerStyle(34);
  pionV0GraphForBBfit->SetMarkerStyle(34);
  pionV0plusTOFGraph->SetMarkerStyle(30);
  pionV0plusTOFGraphForBBfit->SetMarkerStyle(30);
  pionTofGraph->SetMarkerColor(3);
  pionTpcGraph->SetMarkerColor(3);
  pionV0Graph->SetMarkerColor(3);
  pionV0GraphForBBfit->SetMarkerColor(3);
  pionV0plusTOFGraph->SetMarkerColor(3);
  pionV0plusTOFGraphForBBfit->SetMarkerColor(3);
  for(Int_t ip = 0; ip < pionTofGraph->GetN(); ip ++) {
    pionTofGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kPion);
    pionTofGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kPion);
  }
  for(Int_t ip = 0; ip < pionTpcGraph->GetN(); ip ++) {
    pionTpcGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kPion);
    pionTpcGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kPion);
  }
  for(Int_t ip = 0; ip < pionV0Graph->GetN(); ip ++) {
    pionV0Graph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kPion);
    pionV0Graph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kPion);
  }
  for(Int_t ip = 0; ip < pionV0GraphForBBfit->GetN(); ip ++) {
    pionV0GraphForBBfit->GetX()[ip] /= AliPID::ParticleMass(AliPID::kPion);
    pionV0GraphForBBfit->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kPion);
  }
  for(Int_t ip = 0; ip < pionV0plusTOFGraph->GetN(); ip ++) {
    pionV0plusTOFGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kPion);
    pionV0plusTOFGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kPion);
  }
  for(Int_t ip = 0; ip < pionV0plusTOFGraphForBBfit->GetN(); ip ++) {
    pionV0plusTOFGraphForBBfit->GetX()[ip] /= AliPID::ParticleMass(AliPID::kPion);
    pionV0plusTOFGraphForBBfit->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kPion);
  }
  //
  protonTofGraph->SetMarkerStyle(21);
  protonTpcGraph->SetMarkerStyle(22);
  protonV0Graph->SetMarkerStyle(34);
  protonV0GraphForBBfit->SetMarkerStyle(34);
  protonV0plusTOFGraph->SetMarkerStyle(30);
  protonV0plusTOFGraphForBBfit->SetMarkerStyle(30);
  protonTofGraph->SetMarkerColor(4);
  protonTpcGraph->SetMarkerColor(4);
  protonV0Graph->SetMarkerColor(4);
  protonV0GraphForBBfit->SetMarkerColor(4);
  protonV0plusTOFGraph->SetMarkerColor(4);
  protonV0plusTOFGraphForBBfit->SetMarkerColor(4);
  for(Int_t ip = 0; ip < protonTofGraph->GetN(); ip ++) {
    protonTofGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kProton);
    protonTofGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kProton);
  }
  for(Int_t ip = 0; ip < protonTpcGraph->GetN(); ip ++) {
    protonTpcGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kProton);
    protonTpcGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kProton);
  }
  for(Int_t ip = 0; ip < protonV0Graph->GetN(); ip ++) {
    protonV0Graph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kProton);
    protonV0Graph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kProton);
  }
  for(Int_t ip = 0; ip < protonV0GraphForBBfit->GetN(); ip ++) {
    protonV0GraphForBBfit->GetX()[ip] /= AliPID::ParticleMass(AliPID::kProton);
    protonV0GraphForBBfit->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kProton);
  }
  for(Int_t ip = 0; ip < protonV0plusTOFGraph->GetN(); ip ++) {
    protonV0plusTOFGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kProton);
    protonV0plusTOFGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kProton);
  }
  for(Int_t ip = 0; ip < protonV0plusTOFGraphForBBfit->GetN(); ip ++) {
    protonV0plusTOFGraphForBBfit->GetX()[ip] /= AliPID::ParticleMass(AliPID::kProton);
    protonV0plusTOFGraphForBBfit->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kProton);
  }
  //
  //
  TGraphErrors * graphAll = new TGraphErrors();
  TList * listColl = new TList();
  if (!useV0s) {
    //listColl->Add(protonTpcGraph);
    //listColl->Add(kaonTpcGraph);  
    listColl->Add(pionTpcGraph);
    listColl->Add(electronGraph);
    listColl->Add(protonTofGraph); 
    listColl->Add(kaonTofGraph);
    listColl->Add(pionTofGraph);
  }
  else {
    listColl->Add(pionV0GraphForBBfit);
    //listColl->Add(electronV0GraphForBBfit);
    listColl->Add(protonV0GraphForBBfit);
    
    //listColl->Add(pionV0plusTOFGraphForBBfit);
    listColl->Add(electronV0plusTOFGraphForBBfit);
    //listColl->Add(protonV0plusTOFGraphForBBfit);
  }
  MergeGraphErrors(graphAll, listColl);
  graphAll->SetNameTitle("beamDataPoints","beamDataPoints");
  
  // Create graph out of TPC only, V0, V0+TOF and TPC+TOF to determine correction functions
  TList * listCollPr = new TList();
  if (useV0s) {
    listCollPr->Add(protonTpcGraph);
    //listCollPr->Add(protonTofGraph);
    
    TGraphErrors * protonV0GraphCut = new TGraphErrors(*protonV0Graph);
    protonV0GraphCut->SetName("protonV0GraphCut");
    for(Int_t ip = 0; ip < protonV0GraphCut->GetN(); ip ++) {
      Double_t mom = protonV0GraphCut->GetX()[ip] * AliPID::ParticleMass(AliPID::kProton);
      //Bool_t removePoint = mom >= 0.6 || mom < 0.35; // Will take TPC protons instead (statistics) for very low momenta
      Bool_t removePoint = mom < (isPPb ? 0.6 : 0.35); // Will take TPC protons instead (statistics) for very low momenta;
                                                       // for pPb low statistics for the moment, therefore, go a bit further with TPC protons
      
      
      if (removePoint) {
        protonV0GraphCut->RemovePoint(ip);
        ip--;
        continue;
      }
    }
    /*
    TGraphErrors * protonV0plusTOFGraphCut = new TGraphErrors(*protonV0plusTOFGraph);
    protonV0plusTOFGraphCut->SetName("protonV0plusTOFGraphCut");
    for(Int_t ip = 0; ip < protonV0plusTOFGraphCut->GetN(); ip ++) {
      Double_t mom = protonV0plusTOFGraphCut->GetX()[ip] * AliPID::ParticleMass(AliPID::kProton);
      Bool_t removePoint = mom < 0.6;
      
      if (removePoint) {
        protonV0plusTOFGraphCut->RemovePoint(ip);
        ip--;
        continue;
      }
    }*/
  
    listCollPr->Add(protonV0GraphCut);
    //listCollPr->Add(protonV0plusTOFGraphCut);
  }
  else {
    listCollPr->Add(protonTpcGraph);
    listCollPr->Add(protonTofGraph);
  }
  TGraphErrors * graphPrCombined = new TGraphErrors();
  MergeGraphErrors(graphPrCombined, listCollPr);
  graphPrCombined->SetMarkerStyle(22);
  graphPrCombined->SetMarkerColor(4);
  graphPrCombined->SetNameTitle("Graph_Protons_Combined", "Graph_Protons_Combined");
  
  TList * listCollPi = new TList();
  if (useV0s) {
    //listCollPi->Add(pionTpcGraph);
    //listCollPi->Add(pionTofGraph);
    
    TGraphErrors * pionV0GraphCut = new TGraphErrors(*pionV0Graph);
    pionV0GraphCut->SetName("pionV0GraphCut");
    for(Int_t ip = 0; ip < pionV0GraphCut->GetN(); ip ++) {
      //Double_t mom = pionV0GraphCut->GetX()[ip] * AliPID::ParticleMass(AliPID::kPion);
      Bool_t removePoint = kFALSE;//mom >= 0.4;
                           //|| mom < 0.2;//TODO NOW NEEDED?
      if (removePoint) {
        pionV0GraphCut->RemovePoint(ip);
        ip--;
        continue;
      }
    }
    /*
    TGraphErrors * pionV0plusTOFGraphCut = new TGraphErrors(*pionV0plusTOFGraph);
    pionV0plusTOFGraphCut->SetName("pionV0plusTOFGraphCut");
    for(Int_t ip = 0; ip < pionV0plusTOFGraphCut->GetN(); ip ++) {
      Double_t mom = pionV0plusTOFGraphCut->GetX()[ip] * AliPID::ParticleMass(AliPID::kPion);
      Bool_t removePoint = mom < 0.4;
      
      if (removePoint) {
        pionV0plusTOFGraphCut->RemovePoint(ip);
        ip--;
        continue;
      }
    }*/
  
    listCollPi->Add(pionV0GraphCut);
    //listCollPi->Add(pionV0plusTOFGraphCut);
  }
  else {
    listCollPi->Add(pionTpcGraph);
    listCollPi->Add(pionTofGraph);
  }
  TGraphErrors * graphPiCombined = new TGraphErrors();
  MergeGraphErrors(graphPiCombined, listCollPi);
  graphPiCombined->SetMarkerStyle(22);
  graphPiCombined->SetMarkerColor(3);
  graphPiCombined->SetNameTitle("Graph_Pions_Combined", "Graph_Pions_Combined");
  
  TList * listCollKa = new TList();
  listCollKa->Add(kaonTpcGraph);
  listCollKa->Add(kaonTofGraph);
  TGraphErrors * graphKaCombined = new TGraphErrors();
  MergeGraphErrors(graphKaCombined, listCollKa);
  graphKaCombined->SetMarkerStyle(22);
  graphKaCombined->SetMarkerColor(5);
  graphKaCombined->SetNameTitle("Graph_Kaons_Combined", "Graph_Kaons_Combined");
  
  TList * listCollEl = new TList();
  if (useV0s) {
    //listCollEl->Add(electronGraph);
    
    TGraphErrors * electronV0GraphCut = new TGraphErrors(*electronV0Graph);
    electronV0GraphCut->SetName("electronV0GraphCut");
    /*
    for(Int_t ip = 0; ip < electronV0GraphCut->GetN(); ip ++) {
      Double_t mom = electronV0GraphCut->GetX()[ip] * AliPID::ParticleMass(AliPID::kElectron);
      Bool_t removePoint = mom > 0.3;
      if (removePoint) {
        electronV0GraphCut->RemovePoint(ip);
        ip--;
        continue;
      }
    }
    
    TGraphErrors * electronV0plusTOFGraphCut = new TGraphErrors(*electronV0plusTOFGraph);
    electronV0plusTOFGraphCut->SetName("electronV0plusTOFGraphCut");
    for(Int_t ip = 0; ip < electronV0plusTOFGraphCut->GetN(); ip ++) {
      Double_t mom = electronV0plusTOFGraphCut->GetX()[ip] * AliPID::ParticleMass(AliPID::kElectron);
      Bool_t removePoint = mom <= 0.4;
      
      if (removePoint) {
        electronV0plusTOFGraphCut->RemovePoint(ip);
        ip--;
        continue;
      }
    }
    */
    listCollEl->Add(electronV0GraphCut);
    //listCollEl->Add(electronV0plusTOFGraphCut);
  }
  else {
    listCollEl->Add(electronGraph);
  }
  TGraphErrors * graphElCombined = new TGraphErrors();
  MergeGraphErrors(graphElCombined, listCollEl);
  graphElCombined->SetMarkerStyle(22);
  graphElCombined->SetMarkerColor(2);
  graphElCombined->SetNameTitle("Graph_Electrons_Combined", "Graph_Electrons_Combined");
  
  
  
  //
  //  return array with all graphs
  //
  TObjArray * arrGraphs = new TObjArray();
  arrGraphs->Add(graphAll);
  //  
  arrGraphs->Add(electronGraph);
  arrGraphs->Add(electronV0Graph);
  arrGraphs->Add(electronV0plusTOFGraph);
  arrGraphs->Add(electronV0GraphForBBfit);
  arrGraphs->Add(pionTpcGraph);
  arrGraphs->Add(pionTofGraph);
  arrGraphs->Add(pionV0Graph);
  arrGraphs->Add(pionV0plusTOFGraph);
  arrGraphs->Add(pionV0GraphForBBfit);
  arrGraphs->Add(kaonTpcGraph);
  arrGraphs->Add(kaonTofGraph);
  arrGraphs->Add(protonTpcGraph);
  arrGraphs->Add(protonTofGraph);
  arrGraphs->Add(protonV0Graph);
  arrGraphs->Add(protonV0plusTOFGraph);
  arrGraphs->Add(protonV0GraphForBBfit);
  //
  arrGraphs->Add(graphElCombined);
  arrGraphs->Add(graphPiCombined);
  arrGraphs->Add(graphKaCombined);
  arrGraphs->Add(graphPrCombined);
  //
  //
  
  canvasQAtpc->SaveAs("splines_QA_ResidualGraphsTPC.root");
  canvasQAtof->SaveAs("splines_QA_ResidualGraphsTOF.root");
  canvasQAv0->SaveAs("splines_QA_ResidualGraphsV0.root");
  canvasQAv0plusTOF->SaveAs("splines_QA_ResidualGraphsV0plusTOF.root");
  
  
  canvasQAv0DeDxPurityEl->SaveAs("V0_dEdx_purity_el.root");
  canvasQAv0DeDxPurityPi->SaveAs("V0_dEdx_purity_Pi.root");
  canvasQAv0DeDxPurityPr->SaveAs("V0_dEdx_purity_Pr.root");

  return arrGraphs;

}


//________________________________________________________________________
TObjArray * AliTPCcalibResidualPID::GetResidualGraphsMC(THnSparseF * histPidQA, const Char_t * /*system*/) {
  //
  // Extracts the residual graphs from THnSparse created from MC (NOT DATA!)
  //
  
  const TString momTitle = "#it{p}_{TPC} (GeV/#it{c})";
  const TString dEdxTitle = "d#it{E}/d#it{x} (arb. unit)";
  
  Int_t cutForFitting = 10;
  Double_t heightFractionForFittingRange = 0.1; 
  
  //
  THnSparse * hist = histPidQA;
  //
  TCanvas * canvasQAmc = new TCanvas("canvasQAmcResGraph","Control canvas for residual graphs (MC)",100,10,1380,800);
  canvasQAmc->Divide(2,2);
  
  for (Int_t i = 1; i <= 4; i++) {
    canvasQAmc->GetPad(i)->SetGrid(1, 1);
    canvasQAmc->GetPad(i)->SetLogz();
    canvasQAmc->GetPad(i)->SetLogx();
  }
  //
  // 1. select and fit electrons
  //
  // In the following, axis 3 will also be limited in range, although nsigma itself is not used. This is to avoid multiple counting of entries
  hist->GetAxis(2)->SetRangeUser(0,0); // electrons
  hist->GetAxis(3)->SetRangeUser(0,0); // electrons (used for nsigma and for the eta correction of the signal)
  TH2D * histElectron = hist->Projection(1,0);
  histElectron->SetName("histElectron");
  histElectron->RebinX(2);
  histElectron->GetXaxis()->SetRangeUser(0.15,8.0);
  TObjArray arr;
  FitSlicesY(histElectron, heightFractionForFittingRange, cutForFitting, "QNRL", &arr);
  TH1D * electronPoints = (TH1D *) arr.At(1)->Clone();
  //
  TGraphErrors * electronGraph = new TGraphErrors(electronPoints);
  electronGraph->SetName("electronGraphMC");
  for(Int_t ip = 0; ip < electronGraph->GetN(); ip ++) {
    Bool_t removePoint = electronGraph->GetY()[ip] < 10 || electronGraph->GetEY()[ip]/electronGraph->GetY()[ip] > 0.03
                         || electronGraph->GetX()[ip] < 0.15;// || electronGraph->GetX()[ip] > 3.1;
    if (removePoint) {
      electronGraph->RemovePoint(ip);
      ip--;
      continue;
    }
  }
  //
  canvasQAmc->cd(1);
  histElectron->SetTitle("electrons");
  histElectron->GetYaxis()->SetRangeUser(50, 120);
  histElectron->GetXaxis()->SetTitle(momTitle.Data());
  histElectron->GetYaxis()->SetTitle(dEdxTitle.Data());
  histElectron->GetXaxis()->SetMoreLogLabels(kTRUE);
  histElectron->GetXaxis()->SetNoExponent(kTRUE);
  histElectron->Draw("colz");
  electronPoints->SetMarkerStyle(24);
  electronPoints->Draw("same");
  //
  // 2. protons
  //
  hist->GetAxis(2)->SetRangeUser(3,3); // protons
  hist->GetAxis(3)->SetRangeUser(3,3); // protons (used for nsigma and for the eta correction of the signal)
  //
  TH2D * histProton = hist->Projection(1,0);
  histProton->SetName("histProton");
  histProton->GetXaxis()->SetRangeUser(0.15,8);
  histProton->GetYaxis()->SetRangeUser(10, hist->GetAxis(1)->GetBinUpEdge(hist->GetAxis(1)->GetNbins()));
  FitSlicesY(histProton, heightFractionForFittingRange, cutForFitting, "QNRL", &arr);
  TH1D * protonPoints = (TH1D *) arr.At(1)->Clone();
  //
  TGraphErrors * protonGraph = new TGraphErrors(protonPoints);
  protonGraph->SetName("protonGraphMC");
  for(Int_t ip = 0; ip < protonGraph->GetN(); ip ++) {
    Bool_t removePoint = protonGraph->GetY()[ip] < 10 || protonGraph->GetEY()[ip]/protonGraph->GetY()[ip] > 0.02 || 
                         protonGraph->GetX()[ip] < 0.15 || protonGraph->GetX()[ip] > 6;
    if (removePoint) {
      protonGraph->RemovePoint(ip);
      ip--;
      continue;
    }
  }
  //
  canvasQAmc->cd(2);
  histProton->SetTitle("protons");
  histProton->GetXaxis()->SetTitle(momTitle.Data());
  histProton->GetYaxis()->SetTitle(dEdxTitle.Data());
  histProton->GetXaxis()->SetMoreLogLabels(kTRUE);
  histProton->GetXaxis()->SetNoExponent(kTRUE);
  histProton->Draw("colz");
  protonPoints->SetMarkerStyle(20);
  protonPoints->Draw("same");
  
  //
  // 3. pions
  //
  hist->GetAxis(2)->SetRangeUser(1,1); // pions
  hist->GetAxis(3)->SetRangeUser(1,1); // pions (used for nsigma and for the eta correction of the signal)
  //
  TH2D * histPion = hist->Projection(1,0);
  histPion->SetName("histPion");
  histPion->GetXaxis()->SetRangeUser(0.15,50);
  FitSlicesY(histPion, heightFractionForFittingRange, cutForFitting, "QNRL", &arr);
  TH1D * pionPoints =  (TH1D *) arr.At(1)->Clone();
  //
  TGraphErrors * pionGraph = new TGraphErrors(pionPoints);
  pionGraph->SetName("pionGraphMC");
  for(Int_t ip = 0; ip < pionGraph->GetN(); ip ++) {
    Bool_t removePoint = pionGraph->GetY()[ip] < 10 || pionGraph->GetEY()[ip]/pionGraph->GetY()[ip] > 0.005
                         || pionGraph->GetX()[ip] < 0.15 || pionGraph->GetX()[ip] > 30;
    if (removePoint) {
      pionGraph->RemovePoint(ip);
      ip--;
      continue;
    }
  }
  //
  canvasQAmc->cd(3);
  histPion->GetYaxis()->SetRangeUser(30, 90);
  histPion->GetXaxis()->SetTitle(momTitle.Data());
  histPion->GetYaxis()->SetTitle(dEdxTitle.Data());
  histPion->GetXaxis()->SetMoreLogLabels(kTRUE);
  histPion->GetXaxis()->SetNoExponent(kTRUE);
  histPion->Draw("colz");
  histPion->SetTitle("pions");
  pionPoints->SetMarkerStyle(20);
  pionPoints->Draw("same");
  
  //
  // 4. kaons
  //
  hist->GetAxis(2)->SetRangeUser(2,2); // kaons
  hist->GetAxis(3)->SetRangeUser(2,2); // kaons (used for nsigma and for the eta correction of the signal)
  //
  TH2D * histKaon = hist->Projection(1,0);
  histKaon->SetName("histKaon");
  histKaon->GetXaxis()->SetRangeUser(0.15,8);
  FitSlicesY(histKaon, heightFractionForFittingRange, cutForFitting, "QNRL", &arr);
  TH1D * kaonPoints = (TH1D*) arr.At(1)->Clone();
  //
  TGraphErrors * kaonGraph = new TGraphErrors(kaonPoints);
  kaonGraph->SetName("kaonGraphMC");
  for(Int_t ip = 0; ip < kaonGraph->GetN(); ip ++) {
    Bool_t removePoint = kaonGraph->GetY()[ip] < 10 || kaonGraph->GetEY()[ip]/kaonGraph->GetY()[ip] > 0.02
                         || kaonGraph->GetX()[ip] < 0.15;
    if (removePoint) {
      kaonGraph->RemovePoint(ip);
      ip--;
      continue;
    }
  }
  //
  canvasQAmc->cd(4);
  histKaon->SetTitle("kaons");
  histKaon->GetYaxis()->SetRangeUser(30, 500);
  histKaon->GetXaxis()->SetTitle(momTitle.Data());
  histKaon->GetYaxis()->SetTitle(dEdxTitle.Data());
  histKaon->GetXaxis()->SetMoreLogLabels(kTRUE);
  histKaon->GetXaxis()->SetNoExponent(kTRUE);
  histKaon->Draw("colz");
  kaonPoints->SetMarkerStyle(20);
  kaonPoints->Draw("same");
  
  
  
  hist->GetAxis(2)->SetRange(0,-1); // RESET RANGES
  hist->GetAxis(3)->SetRange(0,-1); // RESET RANGES
  
  
  // Clone graphs to have graphs with the same names as for the data case -> same subsequent functions can be used to determine
  // the correction functions and, finally, the response functions.
  // In addition: Remove some points, which are expected to deviate from the Bethe-Bloch curve, from the original graphs, so that
  // these graphs can be used for the BB fit
  TGraphErrors * graphPrCombined = new TGraphErrors(*protonGraph);
  graphPrCombined->SetMarkerStyle(22);
  graphPrCombined->SetMarkerColor(4);
  graphPrCombined->SetNameTitle("Graph_Protons_Combined", "Graph_Protons_Combined");
  
  TGraphErrors * graphPiCombined = new TGraphErrors(*pionGraph);
  graphPiCombined->SetMarkerStyle(22);
  graphPiCombined->SetMarkerColor(3);
  graphPiCombined->SetNameTitle("Graph_Pions_Combined", "Graph_Pions_Combined");
  
  TGraphErrors * graphKaCombined = new TGraphErrors(*kaonGraph);
  graphKaCombined->SetMarkerStyle(22);
  graphKaCombined->SetMarkerColor(5);
  graphKaCombined->SetNameTitle("Graph_Kaons_Combined", "Graph_Kaons_Combined");
  
  TGraphErrors * graphElCombined = new TGraphErrors(*electronGraph);
  graphElCombined->SetMarkerStyle(22);
  graphElCombined->SetMarkerColor(2);
  graphElCombined->SetNameTitle("Graph_Electrons_Combined", "Graph_Electrons_Combined");
  
  for(Int_t ip = 0; ip < electronGraph->GetN(); ip ++) {
    Bool_t removePoint = electronGraph->GetX()[ip] < 2.5;// || electronGraph->GetX()[ip] > 5.0;
    //Bool_t removePoint = electronGraph->GetX()[ip] < 1.2 || electronGraph->GetX()[ip] > 5.0;
    if (removePoint) {
      electronGraph->RemovePoint(ip);
      ip--;
      continue;
    }
  }
  
  for(Int_t ip = 0; ip < protonGraph->GetN(); ip ++) {
    Bool_t removePoint = protonGraph->GetX()[ip] < 0.9
                         || protonGraph->GetX()[ip] > 2.5;//3.5; //TODO NEW in order not to get into conflict with pions/kaons
    if (removePoint) {
      protonGraph->RemovePoint(ip);
      ip--;
      continue;
    }
  }
  
  for(Int_t ip = 0; ip < pionGraph->GetN(); ip ++) {
    Bool_t removePoint = pionGraph->GetX()[ip] < 2.1;//TODO APR18 0.8;
    if (removePoint) {
      pionGraph->RemovePoint(ip);
      ip--;
      continue;
    }
  }
  
  for(Int_t ip = 0; ip < kaonGraph->GetN(); ip ++) {
    Bool_t removePoint = kaonGraph->GetX()[ip] < 1.25 || kaonGraph->GetX()[ip] > 4.0;// < 0.7;// || kaonGraph->GetX()[ip] > 4;
    if (removePoint) {
      kaonGraph->RemovePoint(ip);
      ip--;
      continue;
    }
  }
  //
  // rescale graphs to betaGamma
  //
  kaonGraph->SetMarkerStyle(22);
  for(Int_t ip = 0; ip < kaonGraph->GetN(); ip ++) {
    kaonGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kKaon);
    kaonGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kKaon);
  }
  for(Int_t ip = 0; ip < graphKaCombined->GetN(); ip ++) {
    graphKaCombined->GetX()[ip] /= AliPID::ParticleMass(AliPID::kKaon);
    graphKaCombined->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kKaon);
  }
  //
  electronGraph->SetMarkerStyle(22);
  electronGraph->SetMarkerColor(2);
  for(Int_t ip = 0; ip < electronGraph->GetN(); ip ++) {
    electronGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kElectron);
    electronGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kElectron);
  }
  for(Int_t ip = 0; ip < graphElCombined->GetN(); ip ++) {
    graphElCombined->GetX()[ip] /= AliPID::ParticleMass(AliPID::kElectron);
    graphElCombined->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kElectron);
  }
  //
  pionGraph->SetMarkerStyle(22);
  pionGraph->SetMarkerColor(3);
  for(Int_t ip = 0; ip < pionGraph->GetN(); ip ++) {
    pionGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kPion);
    pionGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kPion);
  }
  for(Int_t ip = 0; ip < graphPiCombined->GetN(); ip ++) {
    graphPiCombined->GetX()[ip] /= AliPID::ParticleMass(AliPID::kPion);
    graphPiCombined->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kPion);
  }
  //
  protonGraph->SetMarkerStyle(22);
  protonGraph->SetMarkerColor(4);
  for(Int_t ip = 0; ip < protonGraph->GetN(); ip ++) {
    protonGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kProton);
    protonGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kProton);
  }
  for(Int_t ip = 0; ip < graphPrCombined->GetN(); ip ++) {
    graphPrCombined->GetX()[ip] /= AliPID::ParticleMass(AliPID::kProton);
    graphPrCombined->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kProton);
  }
  //
  //
  
  // Store graphs without further restriction
  TGraphErrors * graphPrAll = new TGraphErrors(*graphPrCombined);
  graphPrAll->SetNameTitle("Protons_MC_all", "Protons_MC_all");
  
  TGraphErrors * graphPiAll = new TGraphErrors(*graphPiCombined);
  graphPiAll->SetNameTitle("Pions_MC_all", "Pions_MC_all");
  
  TGraphErrors * graphKaAll = new TGraphErrors(*graphKaCombined);
  graphKaAll->SetNameTitle("Kaons_MC_all", "Kaons_MC_all");
  
  TGraphErrors * graphElAll = new TGraphErrors(*graphElCombined);
  graphElAll->SetNameTitle("Electrons_MC_all", "Electrons_MC_all");
  
  //
  //
  TGraphErrors * graphAll = new TGraphErrors();
  TList * listColl = new TList();
  listColl->Add(electronGraph);
  listColl->Add(protonGraph); 
  listColl->Add(kaonGraph);
  listColl->Add(pionGraph);
  MergeGraphErrors(graphAll, listColl);
  graphAll->SetNameTitle("beamDataPoints","beamDataPoints");
  
  //
  //  return array with all graphs
  //
  TObjArray * arrGraphs = new TObjArray();
  arrGraphs->Add(graphAll);
  //  
  arrGraphs->Add(electronGraph);
  arrGraphs->Add(pionGraph);
  arrGraphs->Add(kaonGraph);
  arrGraphs->Add(protonGraph);
  //
  arrGraphs->Add(graphElAll);
  arrGraphs->Add(graphPiAll);
  arrGraphs->Add(graphKaAll);
  arrGraphs->Add(graphPrAll);  
  //
  arrGraphs->Add(graphElCombined);
  arrGraphs->Add(graphPiCombined);
  arrGraphs->Add(graphKaCombined);
  arrGraphs->Add(graphPrCombined);
  //
  //
  
  canvasQAmc->SaveAs("splines_QA_ResidualGraphsTPC.root");

  return arrGraphs;

}


//________________________________________________________________________
TObjArray * AliTPCcalibResidualPID::GetResponseFunctions(TF1* parametrisation, TObjArray* inputGraphs, const Char_t * type, const Char_t * period, const Char_t * pass, const Char_t * system, const Char_t* dedxtype) {
  //
  //
  
  Bool_t isMC = kFALSE;
  if (strcmp(type, "MC") == 0) {
    isMC = kTRUE;
  }
  else if (strcmp(type, "DATA") == 0) {
    isMC = kFALSE;
  }
  else {
    Printf("ERROR - GetResponseFunctions: Unknown type \"%s\" - must be \"MC\" or \"DATA\"!", type);
    
    return 0x0;
  }
  
  Bool_t isPbPb = kFALSE;
  Bool_t isPPb = kFALSE;
  if (strcmp(system, "PBPB") == 0) {
    Printf("PbPb detected - Applying special handling!");
    isPbPb = kTRUE;
  }
  else if (strcmp(system, "PPB") == 0 || strcmp(system, "PBP") == 0) {
    Printf("p-Pb/Pb-p detected - Applying special handling!");
    isPPb = kTRUE;
  }
  
  TCanvas* canvasQA[4] = { 0x0, };
  for (Int_t i = 0; i < 4; i++) {
    canvasQA[i] = new TCanvas(Form("canvasQAres_%d", i), "Control canvas for residual polynomial",100,10,1380,800);
    canvasQA[i]->SetGrid(1, 1);
    canvasQA[i]->SetLogx();
  }
  /*
  TCanvas * canvasQA = new TCanvas("canvasQAres","Control canvas for residual polynomial",100,10,1380,800);
  canvasQA->Divide(2,2);
  
  for (Int_t i = 1; i <= 4; i++) {
    canvasQA->GetPad(i)->SetGrid(1, 1);
  }
  */
  //
  // smallest needed (100MeV proton) bg = 0.1/0.938(=AliPID::ParticleMass(AliPID::kProton)) = 0.1, 
  // largest needed 100 GeV electron, bg = 100/0.000511(=AliPID::ParticleMass(AliPID::kElectron)) = 200 000
  //
  Double_t from = 0.1;
  Double_t to = 200000;
  
  TF1* funcBB = new TF1(*parametrisation);
  funcBB->SetLineColor(2);
  funcBB->SetLineWidth(1);
  
  TString fitFuncString = "[0] + [1] / x + [2] / (x**2) + [3] / (x**3) + [4] / (x**4) + [5] / (x**5) + [6] * x";
  TString fitFuncString2 = "pol6";
  //
  // 1. extract proton corrections
  //
  TGraphErrors *  graphProtonTPC = (TGraphErrors *) inputGraphs->FindObject("Graph_Protons_Combined");
  TGraphErrors *  graphProtonTPCfinal = new TGraphErrors(*graphProtonTPC);
  graphProtonTPCfinal->SetName("graphProtonTPCfinal");
  for(Int_t ip = 0; ip < graphProtonTPC->GetN(); ip ++) {
    graphProtonTPC->GetY()[ip] -= funcBB->Eval(graphProtonTPC->GetX()[ip]);
    graphProtonTPC->GetY()[ip] /= funcBB->Eval(graphProtonTPC->GetX()[ip]);
    graphProtonTPC->GetEY()[ip] /= funcBB->Eval(graphProtonTPC->GetX()[ip]);;
  }
  canvasQA[0]->cd();
  //canvasQA->cd(1);
  //gPad->SetLogx();
  graphProtonTPC->GetXaxis()->SetTitle("#beta#gamma");
  graphProtonTPC->GetXaxis()->SetMoreLogLabels(kTRUE);
  graphProtonTPC->GetXaxis()->SetNoExponent(kTRUE);
  graphProtonTPC->GetYaxis()->SetLabelSize(0.04);
  graphProtonTPC->GetYaxis()->SetTitleOffset(0.85);
  graphProtonTPC->GetXaxis()->SetTitleOffset(1.0);
  graphProtonTPC->GetYaxis()->SetTitle("(data - fit) / fit");
  graphProtonTPC->Draw("ap");
  
  if (graphProtonTPC->GetN() <= 0)  {
    Printf("ERROR - GetResponseFunctions: Proton graph has no data points!");
    
    return 0x0;
  }
  graphProtonTPC->Sort(); // Sort points along x. Next, the very first point will be used to determine the starting point of the correction function
  TF1 * funcCorrProton = new TF1("funcCorrProton", fitFuncString2.Data(), TMath::Max(graphProtonTPC->GetX()[0], 0.15),1.0); // TODO BEN was 0.18 - 0.85
  graphProtonTPC->Fit(funcCorrProton, "QREX0M");
  //
  // 2. extract kaon corrections
  //
  TGraphErrors *  graphKaonTPC = (TGraphErrors *) inputGraphs->FindObject("Graph_Kaons_Combined");
  TGraphErrors *  graphKaonTPCfinal = new TGraphErrors(*graphKaonTPC);
  graphKaonTPCfinal->SetName("graphKaonTPCfinal");
  for(Int_t ip = 0; ip < graphKaonTPC->GetN(); ip ++) {
    graphKaonTPC->GetY()[ip] -= funcBB->Eval(graphKaonTPC->GetX()[ip]);
    graphKaonTPC->GetY()[ip] /= funcBB->Eval(graphKaonTPC->GetX()[ip]);
    graphKaonTPC->GetEY()[ip] /= funcBB->Eval(graphKaonTPC->GetX()[ip]);;
  }
  canvasQA[1]->cd();
  //canvasQA->cd(2);
  //gPad->SetLogx();
  graphKaonTPC->GetXaxis()->SetTitle("#beta#gamma");
  graphKaonTPC->GetYaxis()->SetTitle("(data - fit) / fit");
  graphKaonTPC->GetXaxis()->SetMoreLogLabels(kTRUE);
  graphKaonTPC->GetXaxis()->SetNoExponent(kTRUE);
  graphKaonTPC->GetYaxis()->SetLabelSize(0.04);
  graphKaonTPC->GetYaxis()->SetTitleOffset(0.85);
  graphKaonTPC->GetXaxis()->SetTitleOffset(1.0);
  graphKaonTPC->Draw("ap");
  
  if (graphKaonTPC->GetN() <= 0)  {
    Printf("ERROR - GetResponseFunctions: Kaon graph has no data points!");
    
    return 0x0;
  }
  graphKaonTPC->Sort(); // Sort points along x. Next, the very first point will be used to determine the starting point of the correction function
  
  TF1 * funcCorrKaon = 0x0;
  
  if (isMC) {
    funcCorrKaon = new TF1("funcCorrKaon", fitFuncString.Data(), TMath::Max(graphKaonTPC->GetX()[0], 0.25), isPPb ? 3.44 : 1.981);
    graphKaonTPC->Fit(funcCorrKaon, "QREX0M");
  }
  else {
    // In case of data there are sometimes problems to fit the shape with one function (could describe the overall deviation,
    // but will not cover all the details).
    // Nevertheless, this shape (including most of the "details") can be fitted with the following approach with two functions
    TF1 * funcCorrKaon1 = new TF1("funcCorrKaon1", fitFuncString.Data(),
                                  TMath::Max(graphKaonTPC->GetX()[0], 0.25), 1.981); 
    graphKaonTPC->Fit(funcCorrKaon1, "QREX0M", "same", TMath::Max(graphKaonTPC->GetX()[0], 0.25), 1.0);
    
    TF1 * funcCorrKaon2 = new TF1("funcCorrKaon2", fitFuncString2.Data(), TMath::Max(graphKaonTPC->GetX()[0], 0.25),  1.981);
    graphKaonTPC->Fit(funcCorrKaon2, "QREX0M", "same", (isMC ? 1.981 : 1.0), 1.981);
    
    funcCorrKaon = new TF1("funcCorrKaon",
                           "funcCorrKaon1 * 0.5*(1.+TMath::Erf((1 - x) / 0.1)) + ([0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x+[5]*x*x*x*x*x+[6]*x*x*x*x*x*x) * 0.5*(1.+TMath::Erf((x - 1) / 0.1))",
                           TMath::Max(graphKaonTPC->GetX()[0], 0.25), 1.981);
  
    for (Int_t i = funcCorrKaon1->GetNpar(), j = 0; i < funcCorrKaon1->GetNpar() + funcCorrKaon2->GetNpar(); i++, j++) {
      funcCorrKaon->SetParameter(j, funcCorrKaon2->GetParameter(j));
    }
    funcCorrKaon->SetLineColor(kRed);
    funcCorrKaon->GetHistogram()->DrawClone("csame");
    //funcCorrKaon->Draw("same");
  }
  /*TODO
  TF1 * funcCorrKaon = new TF1("funcCorrKaon", fitFuncString.Data(),//TODO BEN was fitFuncString2
                               TMath::Max(graphKaonTPC->GetX()[0], 0.25), (isMC ? 1.981 : 1.0)); //TODO BEN was 0.79 for data and 1.45 for MC
  graphKaonTPC->Fit(funcCorrKaon, "QREX0M");
  */
  //
  // 3. extract pion corrections
  //
  TGraphErrors *  graphPionTPC = (TGraphErrors *)  inputGraphs->FindObject("Graph_Pions_Combined");
  TGraphErrors *  graphPionTPCfinal = new TGraphErrors(*graphPionTPC);
  graphPionTPCfinal->SetName("graphPionTPCfinal");
  for(Int_t ip = 0; ip < graphPionTPC->GetN(); ip ++) {
    graphPionTPC->GetY()[ip] -= funcBB->Eval(graphPionTPC->GetX()[ip]);
    graphPionTPC->GetY()[ip] /= funcBB->Eval(graphPionTPC->GetX()[ip]);
    graphPionTPC->GetEY()[ip] /= funcBB->Eval(graphPionTPC->GetX()[ip]);
  }
  canvasQA[2]->cd();
  //canvasQA->cd(3);
  //gPad->SetLogx();
  graphPionTPC->GetXaxis()->SetTitle("#beta#gamma");
  graphPionTPC->GetYaxis()->SetTitle("(data - fit) / fit");
  graphPionTPC->GetXaxis()->SetMoreLogLabels(kTRUE);
  graphPionTPC->GetXaxis()->SetNoExponent(kTRUE);
  graphPionTPC->GetYaxis()->SetLabelSize(0.04);
  graphPionTPC->GetYaxis()->SetTitleOffset(0.85);
  graphPionTPC->GetXaxis()->SetTitleOffset(1.0);
  graphPionTPC->Draw("ap");
  
  if (graphPionTPC->GetN() <= 0)  {
    Printf("ERROR - GetResponseFunctions: Pion graph has no data points!");
    
    return 0x0;
  }
  graphPionTPC->Sort(); // Sort points along x. Next, the very first point will be used to determine the starting point of the correction function
  // In case of PbPb (data, not MC) only correct down to 0.2 GeV/c; otherwise: contamination due to electrons
  TF1 * funcCorrPion = new TF1("funcCorrPion", fitFuncString.Data(), ((isPbPb && !isMC) ? (0.2 / AliPID::ParticleMass(AliPID::kPion)) : graphPionTPC->GetX()[0]), isMC ? (isPPb ? 12.8 : (isPbPb ? 12.8 : 7.1)) : (isPPb ? 7.68 : 7.1)); 
  graphPionTPC->Fit(funcCorrPion, "QREX0M");
  //
  // 4. extract electron corrections
  //
  TGraphErrors *  graphElectronTPC = (TGraphErrors *)  inputGraphs->FindObject("Graph_Electrons_Combined");
  TGraphErrors *  graphElectronTPCfinal = new TGraphErrors(*graphElectronTPC);
  graphElectronTPCfinal->SetName("graphElectronTPCfinal");
  for(Int_t ip = 0; ip < graphElectronTPC->GetN(); ip ++) {
    graphElectronTPC->GetY()[ip] -= funcBB->Eval(graphElectronTPC->GetX()[ip]);
    graphElectronTPC->GetY()[ip] /= funcBB->Eval(graphElectronTPC->GetX()[ip]);
    graphElectronTPC->GetEY()[ip] /= funcBB->Eval(graphElectronTPC->GetX()[ip]);;
  }
  canvasQA[3]->cd();
  //canvasQA->cd(4);
  //gPad->SetLogx();
  graphElectronTPC->GetXaxis()->SetTitle("#beta#gamma");
  graphElectronTPC->GetYaxis()->SetTitle("(data - fit) / fit");
  graphElectronTPC->GetXaxis()->SetMoreLogLabels(kTRUE);
  graphElectronTPC->GetXaxis()->SetNoExponent(kTRUE);
  graphElectronTPC->GetYaxis()->SetLabelSize(0.04);
  graphElectronTPC->GetYaxis()->SetTitleOffset(0.85);
  graphElectronTPC->GetXaxis()->SetTitleOffset(1.0);
  graphElectronTPC->Draw("ap");
  
  if (graphElectronTPC->GetN() <= 0)  {
    Printf("ERROR - GetResponseFunctions: Electron graph has no data points!");
    
    return 0x0;
  }
  graphElectronTPC->Sort(); // Sort points along x. Next, the very first point will be used to determine the starting point of the correction function
  // In case of PbPb (data, not MC) only correct down to 0.2 GeV/c; otherwise: contamination due to pions
  TF1 * funcCorrElectron = new TF1("funcCorrElectron", fitFuncString.Data(), (!isMC && isPbPb ? (0.2 / AliPID::ParticleMass(AliPID::kElectron)) :graphElectronTPC->GetX()[0]), (isMC ? 3565 : (isPPb ? 2900 : 1920/*970*/)));// TODO was 1800 for pp data
  // NOTE: For data, the results are almost the same for fitFuncString and fitFuncString2. Maybe, fitFuncString2 is slightly better.
  graphElectronTPC->Fit(funcCorrElectron, "QREX0M");
  //
  // EXTRACT GRAPHS AND PUT THEM TO THE TOBJARRAY
  //
  const Int_t nBins = 500;
  Double_t xBetaGamma[nBins];
  Double_t yProton[nBins];
  Double_t yKaon[nBins];
  Double_t yPion[nBins];
  Double_t yElectron[nBins];
  Double_t yDefault[nBins];
  //
  // 
  //
  xBetaGamma[0] = from;
  Double_t factor = pow(to/from, 1./nBins);
  
  //
  for(Int_t kk = 0; kk < nBins; kk++) {
    if (kk > 0) xBetaGamma[kk] = factor * xBetaGamma[kk-1];
    yProton[kk] =  funcBB->Eval(xBetaGamma[kk])/50.;
    yPion[kk] =  funcBB->Eval(xBetaGamma[kk])/50.;
    yKaon[kk] =  funcBB->Eval(xBetaGamma[kk])/50.;
    yElectron[kk] =  funcBB->Eval(xBetaGamma[kk])/50.;
    yDefault[kk] = funcBB->Eval(xBetaGamma[kk])/50.;

    // Added by Ben
    Double_t widthFactor = 0.020;
    Double_t smoothProton   = 0.5*(TMath::Erf((funcCorrProton->GetXmax()-xBetaGamma[kk])*AliPID::ParticleMass(AliPID::kProton)/widthFactor) + 1);
    Double_t smoothKaon     = 0.5*(TMath::Erf((funcCorrKaon->GetXmax()-xBetaGamma[kk])*AliPID::ParticleMass(AliPID::kKaon)/widthFactor) + 1);
    Double_t smoothPion     = 0.5*(TMath::Erf((funcCorrPion->GetXmax()-xBetaGamma[kk])*AliPID::ParticleMass(AliPID::kPion)/widthFactor) + 1);
    Double_t smoothElectron = 0.5*(TMath::Erf((funcCorrElectron->GetXmax()-xBetaGamma[kk])*AliPID::ParticleMass(AliPID::kElectron)/widthFactor) + 1);

    if (xBetaGamma[kk] > funcCorrProton->GetXmax())
      yProton[kk] *= (1 + smoothProton*funcCorrProton->Eval(funcCorrProton->GetXmax()));
    // Correction is so large that one runs into trouble at low bg for Protons.
    // The deviation is smaller if one takes the lower bound of the correction function
    else if (xBetaGamma[kk] < funcCorrProton->GetXmin())
      yProton[kk] *= (1 + smoothProton*funcCorrProton->Eval(funcCorrProton->GetXmin()));
    else
      yProton[kk] *= (1 + smoothProton*funcCorrProton->Eval(xBetaGamma[kk]));
    
    if (xBetaGamma[kk] > funcCorrKaon->GetXmax())
      yKaon[kk] *= (1 + smoothKaon*funcCorrKaon->Eval(funcCorrKaon->GetXmax()));
    // Correction is so large that one runs into trouble at low bg for Kaons.
    // The deviation is smaller if one takes the lower bound of the correction function
    else if (xBetaGamma[kk] < funcCorrKaon->GetXmin())
      yKaon[kk] *= (1 + smoothKaon*funcCorrKaon->Eval(funcCorrKaon->GetXmin()));
    else
      yKaon[kk] *= (1 + smoothKaon*funcCorrKaon->Eval(xBetaGamma[kk]));

    if (xBetaGamma[kk] > funcCorrElectron->GetXmax())
      yElectron[kk] *= (1 + smoothElectron*funcCorrElectron->Eval(funcCorrElectron->GetXmax()));
    else if (xBetaGamma[kk] < funcCorrElectron->GetXmin())
      yElectron[kk] *= (1 + smoothElectron*funcCorrElectron->Eval(funcCorrElectron->GetXmin()));
    else
      yElectron[kk] *= (1 + smoothElectron*funcCorrElectron->Eval(xBetaGamma[kk]));
    
    // Only true for LHC10d.pass?: Seems not to be needed because any (very small) deviations are most likely due to electron contamination(BEN)
    if (xBetaGamma[kk] > funcCorrPion->GetXmax())
      yPion[kk] *= (1 + smoothPion*funcCorrPion->Eval(funcCorrPion->GetXmax()));
    else if (xBetaGamma[kk] < funcCorrPion->GetXmin())
      yPion[kk] *= (1 + smoothPion*funcCorrPion->Eval(funcCorrPion->GetXmin()));
    else
      yPion[kk] *= (1 + smoothPion*funcCorrPion->Eval(xBetaGamma[kk]));

    
    /* Removed by Ben
    Double_t smoothProton = 0.5*(TMath::Erf((funcCorrProton->GetXmax()-xBetaGamma[kk])/0.002) + 1);
    Double_t smoothKaon   = 0.5*(TMath::Erf((funcCorrKaon->GetXmax()-xBetaGamma[kk])/0.002) + 1);
    Double_t smoothPion   = 0.5*(TMath::Erf((funcCorrPion->GetXmax()-xBetaGamma[kk])/0.002) + 1);

    if (xBetaGamma[kk] < funcCorrProton->GetXmax() && xBetaGamma[kk] > funcCorrProton->GetXmin()) yProton[kk] *= (1 + smoothProton*funcCorrProton->Eval(xBetaGamma[kk])); 
    if (xBetaGamma[kk] < funcCorrKaon->GetXmax() && xBetaGamma[kk] > funcCorrKaon->GetXmin()) yKaon[kk] *= (1 + smoothKaon*funcCorrKaon->Eval(xBetaGamma[kk])); 
    if (xBetaGamma[kk] < funcCorrPion->GetXmax() && xBetaGamma[kk] > funcCorrPion->GetXmin()) yPion[kk] *= (1 + smoothPion*funcCorrPion->Eval(xBetaGamma[kk])); 
    //if (xBetaGamma[kk] < funcCorrElectron->GetXmax()) yElectron[kk] *= (1 + funcCorrElectron->Eval(xBetaGamma[kk])); 
    //
    if (xBetaGamma[kk] < funcCorrProton->GetXmin()) yProton[kk] *= (1 + funcCorrProton->Eval(funcCorrProton->GetXmin())); 
    if (xBetaGamma[kk] < funcCorrKaon->GetXmin()) yKaon[kk] *= (1 + funcCorrKaon->Eval(funcCorrKaon->GetXmin())); 
    if (xBetaGamma[kk] < funcCorrPion->GetXmin()) yPion[kk] *= (1 + funcCorrPion->Eval(funcCorrPion->GetXmin())); 
    */
  }
  //
  TGraph * graphProton = new TGraph(nBins, xBetaGamma, yProton);
  TGraph * graphElectron = new TGraph(nBins, xBetaGamma, yElectron);
  TGraph * graphKaon = new TGraph(nBins, xBetaGamma, yKaon);
  TGraph * graphPion = new TGraph(nBins, xBetaGamma, yPion);
  TGraph * graphDefault = new TGraph(nBins, xBetaGamma, yDefault);
  //
  //
  TSpline3 * splineProton = new TSpline3("splineProton", graphProton);
  TSpline3 * splineElectron = new TSpline3("splineElectron", graphElectron);
  TSpline3 * splineKaon = new TSpline3("splineKaon", graphKaon);
  TSpline3 * splinePion = new TSpline3("splinePion", graphPion);
  TSpline3 * splineDefault = new TSpline3("splineDefault", graphDefault);
  //
  //
  splineProton->SetNameTitle(Form("TSPLINE3_%s_PROTON_%s_%s_%s_%sMEAN",type,period,pass,system,dedxtype),
                             Form("TSPLINE3_%s_PROTON_%s_%s_%s_%sMEAN",type,period,pass,system,dedxtype));
  splineElectron->SetNameTitle(Form("TSPLINE3_%s_ELECTRON_%s_%s_%s_%sMEAN",type,period,pass,system,dedxtype),
                               Form("TSPLINE3_%s_ELECTRON_%s_%s_%s_%sMEAN",type,period,pass,system,dedxtype));
  splineKaon->SetNameTitle(Form("TSPLINE3_%s_KAON_%s_%s_%s_%sMEAN",type,period,pass,system,dedxtype),
                           Form("TSPLINE3_%s_KAON_%s_%s_%s_%sMEAN",type,period,pass,system,dedxtype));
  splinePion->SetNameTitle(Form("TSPLINE3_%s_PION_%s_%s_%s_%sMEAN",type,period,pass,system,dedxtype),
                           Form("TSPLINE3_%s_PION_%s_%s_%s_%sMEAN",type,period,pass,system,dedxtype));
  splineDefault->SetNameTitle(Form("TSPLINE3_%s_ALL_%s_%s_%s_%sMEAN",type,period,pass,system,dedxtype),
                              Form("TSPLINE3_%s_ALL_%s_%s_%s_%sMEAN",type,period,pass,system,dedxtype));
  //
  TObjArray * arrResponse = new TObjArray();
  arrResponse->SetName("TPCpidResponseFunctions");
  arrResponse->AddLast(splineProton);
  arrResponse->AddLast(splineElectron);
  arrResponse->AddLast(splineKaon);
  arrResponse->AddLast(splinePion);
  arrResponse->AddLast(splineDefault);
  //
  
  // Draw deviation from final results
  for(Int_t ip = 0; ip < graphProtonTPCfinal->GetN(); ip ++) {
    graphProtonTPCfinal->GetY()[ip] -= 50. * splineProton->Eval(graphProtonTPCfinal->GetX()[ip]);
    graphProtonTPCfinal->GetY()[ip] /= 50. * splineProton->Eval(graphProtonTPCfinal->GetX()[ip]);
    graphProtonTPCfinal->GetEY()[ip] /= 50. * splineProton->Eval(graphProtonTPCfinal->GetX()[ip]);;
  }
  canvasQA[0]->cd();
  //canvasQA->cd(1);
  graphProtonTPCfinal->SetMarkerStyle(26);
  graphProtonTPCfinal->Draw("psame");
  
  for(Int_t ip = 0; ip < graphKaonTPCfinal->GetN(); ip ++) {
    graphKaonTPCfinal->GetY()[ip] -= 50. * splineKaon->Eval(graphKaonTPCfinal->GetX()[ip]);
    graphKaonTPCfinal->GetY()[ip] /= 50. * splineKaon->Eval(graphKaonTPCfinal->GetX()[ip]);
    graphKaonTPCfinal->GetEY()[ip] /= 50. * splineKaon->Eval(graphKaonTPCfinal->GetX()[ip]);;
  }
  canvasQA[1]->cd();
  //canvasQA->cd(2);
  graphKaonTPCfinal->SetMarkerStyle(26);
  graphKaonTPCfinal->Draw("psame");
  
  for(Int_t ip = 0; ip < graphPionTPCfinal->GetN(); ip ++) {
    graphPionTPCfinal->GetY()[ip] -= 50. * splinePion->Eval(graphPionTPCfinal->GetX()[ip]);
    graphPionTPCfinal->GetY()[ip] /= 50. * splinePion->Eval(graphPionTPCfinal->GetX()[ip]);
    graphPionTPCfinal->GetEY()[ip] /= 50. * splinePion->Eval(graphPionTPCfinal->GetX()[ip]);;
  }
  canvasQA[2]->cd();
  //canvasQA->cd(3);
  graphPionTPCfinal->SetMarkerStyle(26);
  graphPionTPCfinal->Draw("psame");
  
  for(Int_t ip = 0; ip < graphElectronTPCfinal->GetN(); ip ++) {
    graphElectronTPCfinal->GetY()[ip] -= 50. * splineElectron->Eval(graphElectronTPCfinal->GetX()[ip]);
    graphElectronTPCfinal->GetY()[ip] /= 50. * splineElectron->Eval(graphElectronTPCfinal->GetX()[ip]);
    graphElectronTPCfinal->GetEY()[ip] /= 50. * splineElectron->Eval(graphElectronTPCfinal->GetX()[ip]);;
  }
  canvasQA[3]->cd();
  //canvasQA->cd(4);
  graphElectronTPCfinal->SetMarkerStyle(26);
  graphElectronTPCfinal->Draw("psame");
  
  TFile* fSave = TFile::Open("splines_QA_ResidualPolynomials.root", "RECREATE");
  fSave->cd();
  for (Int_t i = 0; i < 4; i++)
    canvasQA[i]->Write();  
  
  fSave->Close();
  //canvasQA->SaveAs("splines_QA_ResidualPolynomials.root");
  
  
  delete funcBB;
  
  return arrResponse;
}


//________________________________________________________________________
TF1* AliTPCcalibResidualPID::FitBB(TObjArray* inputGraphs, Bool_t isMC, Bool_t isPPb, const Bool_t useV0s, 
                                   const Double_t * initialParameters, AliTPCcalibResidualPID::FitType fitType) {
  //
  // Fit Bethe-Bloch parametrisation to data points (inputGraphs)
  //
  const Int_t nPar = 6;
  TGraphErrors * graphAll = (TGraphErrors *) inputGraphs->FindObject("beamDataPoints");
  //
  Float_t from = 0.9; //TODO ADJUST -> Very important
  Float_t to = graphAll->GetXaxis()->GetXmax() * 2;
  
  
  
  TF1* funcBB = 0x0;
  
  Double_t parametersBBForward[nPar] = { 0, };
  
  TVirtualFitter::SetMaxIterations(5e6);
  
  if (fitType == AliTPCcalibResidualPID::kSaturatedLund) {
    printf("Fit function: Saturated Lund\n");
    
    funcBB = new TF1("SaturatedLund", SaturatedLund, from, to, nPar);
    //Double_t parametersBB[nPar] = {34.0446, 8.42221, 4.16724, 1.29473, 80.6663, 0}; //Xianguos values
    //Double_t parametersBB[nPar] = {35.5, 8.7, 2.0, 1.09, 75.6, 0}; // No saturation
    Double_t parametersBB[nPar] = {61.0, 8.7, 1.86, 0.85, 113.4, -38}; // Yields reasonable results for data and MC ~ all periods
    
    if (isPPb && !isMC) {
      parametersBB[0] = 51.6;
      parametersBB[1] = 9.7;
      parametersBB[2] = 1.62;
      parametersBB[3] = 0.99;
      parametersBB[4] = 104.4;
      parametersBB[5] = -27.0;
    }
    if (isMC) {
      parametersBB[0] = 41.4;
      parametersBB[1] = 8.6;
      parametersBB[2] = 2.2;
      parametersBB[3] = 0.92;
      parametersBB[4] = 90.4;
      parametersBB[5] = -20.0;
    }
    funcBB->SetParameters(parametersBB);
    
    Double_t parameterErrorsBB[nPar] = {5., 0.5, 0.2, 0.05, 10, 10};
    funcBB->SetParErrors(parameterErrorsBB);
    
    for (Int_t i = 0; i < nPar; i++)
      parametersBBForward[i] = parametersBB[i];
  }
  else if (fitType == AliTPCcalibResidualPID::kLund) {
    printf("Fit function: Lund\n");
    printf("****WARNING: Fit settings not tuned for this fit function, only for saturated Lund!\n");

    funcBB = new TF1("SaturatedLund", SaturatedLund, from, to, nPar);
    //Double_t parametersBB[nPar] = {34.0446, 8.42221, 4.16724, 1.29473, 80.6663, 0}; //Xianguos values
    //Double_t parametersBB[nPar] = {35.5, 8.7, 2.0, 1.09, 75.6, 0}; // No saturation
    Double_t parametersBB[nPar] = {35.0, 7., 1.86, 1., 75., 0}; // Yields reasonable results for data and MC ~ all periods
    
    funcBB->SetParameters(parametersBB);
    
    Double_t parameterErrorsBB[nPar] = {5., 0.5, 0.2, 0.05, 10, 10};
    funcBB->SetParErrors(parameterErrorsBB);
    
    funcBB->FixParameter(5, 0.0); // No saturation
    
    for (Int_t i = 0; i < nPar; i++)
      parametersBBForward[i] = parametersBB[i];
  }
  else if (fitType == kAlephWithAdditionalParam) {
    printf("Fit function: Aleph with additional parameter\n");
    printf("****WARNING: Fit settings not tuned for this fit function, only for saturated Lund!\n");
    
    // NOTE: This form is equivalent to the original form, but with parameter [2] redefined for better numerical stability.
    // The additional parameter [5] has been introduced later and is unity originally. It seems not to be needed and is, thus,
    // fixed to unity
    funcBB = new TF1("funcAleph",
                     "[0]/TMath::Power(TMath::Sqrt(1. + x*x)/x , [3])*([1]-[2]-[5]*TMath::Power(TMath::Sqrt(1. + x*x)/x , [3])-TMath::Log(1. + TMath::Power(x, -[4])*TMath::Exp(-[2])))", from, to);
                     //OLD"[0]*([1]*TMath::Power(TMath::Sqrt(1 + x*x)/x , [3]) - [5] - TMath::Power(TMath::Sqrt(1 + x*x)/x , [3])*TMath::Log(TMath::Exp([2]) + 1/TMath::Power(x, [4])))", from, to); 
    
    Double_t parametersBB[nPar] = {2.6,14.3,-15,2.2,2.7, 0.06};
    //OLD with different sign for [3]: Double_t parametersBB[nPar] = {0.0762*50.3,10.632,TMath::Log(1.34e-05),1.863,1.948, 1};
    funcBB->SetParameters(parametersBB);
    
    for (Int_t i = 0; i < nPar; i++)
      parametersBBForward[i] = parametersBB[i];
  }
  else if (fitType == AliTPCcalibResidualPID::kAleph) {
    printf("Fit function: Aleph\n");
    printf("****WARNING: Fit settings not tuned for this fit function, only for saturated Lund!\n");
    
    // NOTE: This form is equivalent to the original form, but with parameter [2] redefined for better numerical stability.
    // The additional parameter [5] has been introduced later and is unity originally. It seems not to be needed and is, thus,
    // fixed to unity
    funcBB = new TF1("funcAleph",
                     "[0]/TMath::Power(x/TMath::Sqrt(1. + x*x) , [3])*([1]-[2]-[5]*TMath::Power(x/TMath::Sqrt(1. + x*x) , [3])-TMath::Log(1. + TMath::Power(x, -[4])*TMath::Exp(-[2])))", from, to);
                     //OLD"[0]*([1]*TMath::Power(TMath::Sqrt(1 + x*x)/x , [3]) - [5] - TMath::Power(TMath::Sqrt(1 + x*x)/x , [3])*TMath::Log(TMath::Exp([2]) + 1/TMath::Power(x, [4])))", from, to); 
    //TEST Double_t parametersBB[nPar] = {1.2, 26.0, -30.0, -2.15, 5.6, 1};
    Double_t parametersBB[nPar] = {1.25, 27.5, -29.0, 2.2, 5.2, 1};
    
    // For [5] floating: Double_t parametersBB[nPar] = {2.42,15.2,-16,-2.24,2.8, 0.057};
    //OLD with different sign for [3] and [5] not fix: Double_t parametersBB[nPar] = {2.6,14.3,-15,2.2,2.7, 0.06};
    //OLD with different sign for [3]: Double_t parametersBB[nPar] = {0.0762*50.3,10.632,TMath::Log(1.34e-05),1.863,1.948, 1};
    funcBB->SetParameters(parametersBB);
    
    // Fix parameter 5 to original value of unity
    funcBB->FixParameter(5, 1); 
    
    for (Int_t i = 0; i < nPar; i++)
      parametersBBForward[i] = parametersBB[i];
  }
  else {
    printf("Error: fitType not supported!\n");
    return 0x0;
  }
  
  funcBB->SetLineColor(2);
  funcBB->SetLineWidth(1);
  
  // Override initial parameters, if user provides some
  if (initialParameters) {
    printf("Setting user initial parameters!\n");
    funcBB->SetParameters(initialParameters);
  }
  
  //
  //
  //
  TGraphErrors * graphDelta = new TGraphErrors(*graphAll);
  graphDelta->SetName("graphDelta");
  
  // In MC case: Remove errors from fit -> Some data points with extremely small errors otherwise completely dominate the fit,
  // but these points could still be wrong due to systematics (low p effects, e.g.)
  if (isMC) {
    for(Int_t ip = 0; ip < graphAll->GetN(); ip++) {
      graphAll->SetPointError(ip, 0, 0);
    }
  }
  else if (useV0s) {
    // Increase weight in plateau by reducing errors. Neccessary because errors are large there compared to other data points
    for(Int_t ip = 0; ip < graphAll->GetN(); ip++) {
      if (graphAll->GetX()[ip] >= 2500) { 
        graphAll->SetPointError(ip, 0, graphAll->GetEY()[ip] / 6.0); 
      }
      // Same for pions in the very rel. rise
      if (graphAll->GetX()[ip] >= 25 && graphAll->GetX()[ip] < 1000) {
        graphAll->SetPointError(ip, 0, graphAll->GetEY()[ip] / 6.0); 
      }
      // Same for protons in the MIP region
      if (graphAll->GetX()[ip] >= 2 && graphAll->GetX()[ip] < 4) {
        graphAll->SetPointError(ip, 0, graphAll->GetEY()[ip] / 6.0); 
      }
    }
  }
  graphAll->Fit(funcBB, "REX0M");  
  funcBB->SetRange(from, to);
  funcBB->GetParameters(parametersBBForward);
  //
  //
  //
  for(Int_t ip = 0; ip < graphDelta->GetN(); ip++) {
    graphDelta->GetY()[ip] -= funcBB->Eval(graphDelta->GetX()[ip]);
    graphDelta->GetY()[ip] /= funcBB->Eval(graphDelta->GetX()[ip]);
    graphDelta->GetEY()[ip] /= funcBB->Eval(graphDelta->GetX()[ip]);
  }
  TCanvas * canvDelta_1 = new TCanvas("canvDelta_1","control histogram for Bethe-Bloch fit 1",100,10,1380,800);
  TCanvas * canvDelta_2 = new TCanvas("canvDelta_2","control histogram for Bethe-Bloch fit 2",100,10,1380,800);
  canvDelta_1->SetGrid(1, 1);
  canvDelta_1->SetLogx();
  canvDelta_2->SetGrid(1, 1);
  canvDelta_2->SetLogx();
  /*
  canvDelta->Divide(2, 1);
  canvDelta->GetPad(1)->SetGrid(1, 1);
  canvDelta->GetPad(1)->SetLogx();
  canvDelta->GetPad(2)->SetGrid(1, 1);
  canvDelta->GetPad(2)->SetLogx();
  */

  canvDelta_1->cd();
  //canvDelta->cd(1);
  TH1F *hBBdummy=new TH1F("hBBdummy","BB fit;#beta#gamma;#LTdE/dx#GT (arb. unit)",100,.8,1e4);
  hBBdummy->SetMinimum(45);
  hBBdummy->SetMaximum(120);
  hBBdummy->GetXaxis()->SetTitleOffset(1.1);
  hBBdummy->SetStats(kFALSE);
  hBBdummy->Draw();

  graphAll->SetTitle("BB multi-graph fit");
  graphAll->SetMarkerStyle(22);
  graphAll->SetMarkerColor(kMagenta);
  graphAll->Draw("p");

  TLegend *leg=new TLegend(.7,.12,.89,isMC?.4:.6);
  leg->SetFillColor(10);
  leg->SetBorderSize(1);

  //overlay individual graphs
  TGraph *gr=0x0;
  if (isMC) {
    gr=(TGraph*)inputGraphs->FindObject("Protons_MC_all");
    gr->SetMarkerColor(kRed);
    gr->SetMarkerStyle(24);
    gr->Draw("p");
    leg->AddEntry(gr,"p (MC)","p");
    gr=(TGraph*)inputGraphs->FindObject("Pions_MC_all");
    gr->SetMarkerColor(kBlue);
    gr->SetMarkerStyle(24);
    gr->Draw("p");
    leg->AddEntry(gr,"#pi (MC)","p");
    gr=(TGraph*)inputGraphs->FindObject("Kaons_MC_all");
    gr->SetMarkerColor(kGreen);
    gr->SetMarkerStyle(24);
    gr->Draw("p");
    leg->AddEntry(gr,"K (MC)","p");
    gr=(TGraph*)inputGraphs->FindObject("Electrons_MC_all");
    gr->SetMarkerColor(kBlack);
    gr->SetMarkerStyle(24);
    gr->Draw("p");
    leg->AddEntry(gr,"e (MC)","p");
  }
  else {
    gr=(TGraph*)inputGraphs->FindObject("protonTpcGraph");
    gr->SetMarkerColor(kRed);
    gr->SetMarkerStyle(26);
    gr->Draw("p");
    leg->AddEntry(gr,"p (TPC)","p");
    gr=(TGraph*)inputGraphs->FindObject("protonTofGraph");
    gr->SetMarkerColor(kRed);
    gr->SetMarkerStyle(25);
    gr->Draw("p");
    leg->AddEntry(gr,"p (TOF)","p");
    gr=(TGraph*)inputGraphs->FindObject("protonV0Graph");//ForBBfit");
    gr->SetMarkerColor(kRed);
    gr->SetMarkerStyle(24);
    gr->Draw("p");
    leg->AddEntry(gr,"p (V0)","p");
    gr=(TGraph*)inputGraphs->FindObject("protonV0plusTOFGraph");//ForBBfit");
    gr->SetMarkerColor(kRed);
    gr->SetMarkerStyle(30);
    gr->Draw("p");
    leg->AddEntry(gr,"p (V0+TOF)","p");
    gr=(TGraph*)inputGraphs->FindObject("pionTpcGraph");
    gr->SetMarkerColor(kBlue);
    gr->SetMarkerStyle(26);
    gr->Draw("p");
    leg->AddEntry(gr,"#pi (TPC)","p");
    gr=(TGraph*)inputGraphs->FindObject("pionTofGraph");
    gr->SetMarkerColor(kBlue);
    gr->SetMarkerStyle(25);
    gr->Draw("p");
    leg->AddEntry(gr,"#pi (TOF)","p");
    gr=(TGraph*)inputGraphs->FindObject("pionV0Graph");//ForBBfit");
    gr->SetMarkerColor(kBlue);
    gr->SetMarkerStyle(24);
    gr->Draw("p");
    leg->AddEntry(gr,"#pi (V0)","p");
    gr=(TGraph*)inputGraphs->FindObject("pionV0plusTOFGraph");//ForBBfit");
    gr->SetMarkerColor(kBlue);
    gr->SetMarkerStyle(30);
    gr->Draw("p");
    leg->AddEntry(gr,"#pi (V0+TOF)","p");
    gr=(TGraph*)inputGraphs->FindObject("kaonTpcGraph");
    gr->SetMarkerColor(kGreen);
    gr->SetMarkerStyle(26);
    gr->Draw("p");
    leg->AddEntry(gr,"K (TPC)","p");
    gr=(TGraph*)inputGraphs->FindObject("kaonTofGraph");
    gr->SetMarkerColor(kGreen);
    gr->SetMarkerStyle(25);
    gr->Draw("p");
    leg->AddEntry(gr,"K (TOF)","p");
    gr=(TGraph*)inputGraphs->FindObject("electronGraph");
    gr->SetMarkerColor(kBlack);
    gr->SetMarkerStyle(25);
    gr->Draw("p");
    leg->AddEntry(gr,"e","p");
    gr=(TGraph*)inputGraphs->FindObject("electronV0Graph");//ForBBfit");
    gr->SetMarkerColor(kBlack);
    gr->SetMarkerStyle(24);
    gr->Draw("p");
    leg->AddEntry(gr,"e (V0)","p");
    gr=(TGraph*)inputGraphs->FindObject("electronV0plusTOFGraph");//ForBBfit");
    gr->SetMarkerColor(kBlack);
    gr->SetMarkerStyle(30);
    gr->Draw("p");
    leg->AddEntry(gr,"e (V0+TOF)","p");
  }
  
  graphAll->GetXaxis()->SetTitle("#beta#gamma");
  graphAll->GetYaxis()->SetTitle("<dE/dx> (arb. unit)");
  graphAll->Draw("p");
  funcBB->GetHistogram()->DrawClone("csame");
 
  leg->AddEntry(graphAll,"used","p");
  leg->AddEntry(funcBB,"fit","l");
  leg->Draw("same");

  canvDelta_2->cd();
  //canvDelta->cd(2);
  TH1F *hBBresdummy=new TH1F("hBBresdummy","residuals of BB fit;#beta#gamma;(data-fit)/fit",100,.8,1e4);
  hBBresdummy->SetMinimum(-0.04);
  hBBresdummy->SetMaximum(0.04);
  hBBresdummy->GetXaxis()->SetTitleOffset(1.1);
  hBBresdummy->GetYaxis()->SetTitleOffset(0.8);
  hBBresdummy->SetStats(kFALSE);
  hBBresdummy->Draw();

  graphDelta->SetTitle("residuals of BB fit to multi-graph");
  graphDelta->GetXaxis()->SetTitle("#beta#gamma");
  graphDelta->GetYaxis()->SetTitle("(data - fit) / fit");
  graphDelta->SetMarkerStyle(22);
  graphDelta->SetMarkerColor(4);
  graphDelta->Draw("p");
  
  TFile* fSave = TFile::Open("splines_QA_BetheBlochFit.root", "RECREATE");
  fSave->cd();
  canvDelta_1->Write();
  canvDelta_2->Write();
  
  TString fitResults = "Fit results:\n";
  for (Int_t i = 0; i < nPar; i++) {
    fitResults.Append(Form("par%d:\t%f +- %f\n", i, funcBB->GetParameters()[i], funcBB->GetParErrors()[i]));
  }
  
  TNamed* settings = new TNamed(fitResults.Data(), fitResults.Data());
  settings->Write();
  fSave->Close();
  
  
  printf("\n\n%s", fitResults.Data());
  
  return funcBB;
}


//________________________________________________________________________
Double_t AliTPCcalibResidualPID::Lund(Double_t* xx, Double_t* par)
{
  // bg is beta*gamma
  const Double_t bg = xx[0];

  const Double_t beta2 = bg*bg / (1.0 + bg*bg);
  
  const Double_t a = par[0];
  const Double_t b = par[1];
  const Double_t c = par[2];
  const Double_t e = par[3];
  const Double_t f = par[4];
  
  const Double_t d = TMath::Exp(c*(a - f) / b);
 
  const Double_t powbg = TMath::Power(1.0 + bg, c);
 
  const Double_t value = a / TMath::Power(beta2,e) +
    b/c * TMath::Log(powbg / (1.0 + d*powbg));
    
  return value;
}


//________________________________________________________________________
Double_t AliTPCcalibResidualPID::SaturatedLund(Double_t* xx, Double_t* par)
{
  const Double_t qq = Lund(xx, par);
  return qq * TMath::Exp(par[5] / qq);
}


//________________________________________________________________________
void AliTPCcalibResidualPID::FitSlicesY(TH2 *hist, Double_t heightFractionForRange, Int_t cutThreshold, TString fitOption, TObjArray *arr)
{
  //heightPercentageForRange
  // custom slices fit
  //

  if (!hist) return;
  if (!arr) return;

  // If old style is to be used
  /*
  hist->FitSlicesY(0, 0, -1, cutThreshold, fitOption.Data(), &arr);
  return;
  */
  
  
  arr->Clear();
  arr->SetOwner();
  arr->Expand(4);
  
  TAxis *axis=hist->GetXaxis();
  const TArrayD *bins = axis->GetXbins();
  TH1D** hList = new TH1D*[4];
  
  for (Int_t i = 0; i < 4; i++) {
    delete gDirectory->FindObject(Form("%s_%d", hist->GetName(), i));
    
    if (bins->fN == 0) {
      hList[i] = new TH1D(Form("%s_%d", hist->GetName(), i), i < 3 ? Form("Fitted value of par[%d]", i) : "Chi2/NDF",
                          hist->GetNbinsX(), axis->GetXmin(), axis->GetXmax());
    } else {
      hList[i] = new TH1D(Form("%s_%d", hist->GetName(), i), i < 3 ? Form("Fitted value of par[%d]", i) : "Chi2/NDF", hist->GetNbinsX(), bins->fArray);
    }
    
    (*arr)[i] = hList[i];
  }
  
  for (Int_t ibin=axis->GetFirst(); ibin<=axis->GetLast(); ++ibin){
    TH1 *h=hist->ProjectionY("_temp",ibin,ibin);
    if (!h)
      continue;
    
    if (h->GetEntries() < cutThreshold) {
      delete h;
      continue;
    }
    
    // Average around maximum bin -> Might catch some outliers
    Int_t maxBin = h->GetMaximumBin();
    Double_t maxVal = h->GetBinContent(maxBin);
    
    if (maxVal < 2) { // It could happen that all entries are in overflow/underflow; don't fit in this case
      delete h;
      continue;
    }
    
    UChar_t usedBins = 1;
    if (maxBin > 1) {
      maxVal += h->GetBinContent(maxBin - 1);
      usedBins++;
    }
    if (maxBin < h->GetNbinsX()) {
      maxVal += h->GetBinContent(maxBin + 1);
      usedBins++;
    }
    maxVal /= usedBins;
    
    Double_t thresholdFraction = heightFractionForRange * maxVal; 
    Int_t lowThrBin = TMath::Max(1, h->FindFirstBinAbove(thresholdFraction));
    Int_t highThrBin = TMath::Min(h->GetNbinsX(), h->FindLastBinAbove(thresholdFraction));
    
    Double_t lowThreshold = h->GetBinCenter(lowThrBin);
    Double_t highThreshold = h->GetBinCenter(highThrBin);

    TFitResultPtr res = h->Fit("gaus", Form("%sS", fitOption.Data()), "", lowThreshold, highThreshold);
    
    if ((Int_t)res == 0) {
      Int_t resBin = ibin;
      hList[0]->SetBinContent(resBin,res->GetParams()[0]);
      hList[0]->SetBinError(resBin,res->GetErrors()[0]);
      hList[1]->SetBinContent(resBin,res->GetParams()[1]);
      hList[1]->SetBinError(resBin,res->GetErrors()[1]);
      hList[2]->SetBinContent(resBin,res->GetParams()[2]);
      hList[2]->SetBinError(resBin,res->GetErrors()[2]);
      hList[3]->SetBinContent(resBin,res->Ndf()>0?res->Chi2()/res->Ndf():0);
    }
    
    delete h;
  }
}


//______________________________________________________________________________
Bool_t AliTPCcalibResidualPID::GetVertexIsOk(AliVEvent* event) const
{
  AliAODEvent* aod = 0x0;
  AliESDEvent* esd = 0x0;
  
  aod = dynamic_cast<AliAODEvent*>(event);
  if (!aod) {
    esd = dynamic_cast<AliESDEvent*>(event);
    
    if (!esd) {
      AliError("Event seems to be neither AOD nor ESD!");
      return kFALSE;
    }
  }
    
  if (fIsPbpOrpPb) {
    const AliVVertex* trkVtx = (aod ? dynamic_cast<const AliVVertex*>(aod->GetPrimaryVertex()) :
                                      dynamic_cast<const AliVVertex*>(esd->GetPrimaryVertex()));
      
    if (!trkVtx || trkVtx->GetNContributors() <= 0)
      return kFALSE;
      
    TString vtxTtl = trkVtx->GetTitle();
    if (!vtxTtl.Contains("VertexerTracks"))
      return kFALSE;
      
    Float_t zvtx = trkVtx->GetZ();
    const AliVVertex* spdVtx = (aod ? dynamic_cast<const AliVVertex*>(aod->GetPrimaryVertexSPD()) :
                                      dynamic_cast<const AliVVertex*>(esd->GetPrimaryVertexSPD()));
    if (spdVtx->GetNContributors() <= 0)
      return kFALSE;
      
    TString vtxTyp = spdVtx->GetTitle();
    Double_t cov[6] = {0};
    spdVtx->GetCovarianceMatrix(cov);
    Double_t zRes = TMath::Sqrt(cov[5]);
    if (vtxTyp.Contains("vertexer:Z") && (zRes > 0.25))
      return kFALSE;
      
    if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ()) > 0.5)
      return kFALSE;

    if (TMath::Abs(zvtx) > fZvtxCutEvent) //Default: 10 cm
      return kFALSE;
      
    return kTRUE;
  }
    
  
  // pp and PbPb
  const AliVVertex* primaryVertex = (aod ? dynamic_cast<const AliVVertex*>(aod->GetPrimaryVertex()) :
                                           dynamic_cast<const AliVVertex*>(esd->GetPrimaryVertexTracks()));
    
  if (!primaryVertex || primaryVertex->GetNContributors() <= 0)
    return kFALSE;
      
  if (TMath::Abs(primaryVertex->GetZ()) >= fZvtxCutEvent) //Default: 10 cm
    return kFALSE;
  
  return kTRUE;
}


//______________________________________________________________________________
void AliTPCcalibResidualPID::FillV0PIDlist(AliESDEvent* event)
{
  //
  // Fill the PID tag list
  //
  
  // If no event forwarded as parameter (default), cast current input event.
  // Dynamic cast to ESD events (DO NOTHING for AOD events)
  if (!event)
    event = dynamic_cast<AliESDEvent *>(InputEvent());
  
  // If this fails, do nothing
  if (!event)
    return;
  
  if (!fV0KineCuts) {
    AliError("V0KineCuts not available!");
    return;
  }
  
  TString beamType(event->GetBeamType());
  
  if (beamType.CompareTo("Pb-Pb") == 0 || beamType.CompareTo("A-A") == 0) {
    fV0KineCuts->SetMode(AliESDv0KineCuts::kPurity, AliESDv0KineCuts::kPbPb);
  }
  else {
    fV0KineCuts->SetMode(AliESDv0KineCuts::kPurity, AliESDv0KineCuts::kPP); 
  }

  // V0 selection
  // set event
  fV0KineCuts->SetEvent(event);

  const Int_t numTracks = event->GetNumberOfTracks();
  fV0tags = new Char_t[numTracks];
  for (Int_t i = 0; i < numTracks; i++)
    fV0tags[i] = 0;
  
  fV0motherIndex = new Int_t[numTracks];
  for (Int_t i = 0; i < numTracks; i++)
    fV0motherIndex[i] = -1;
  
  fV0motherPDG = new Int_t[numTracks];
  for (Int_t i = 0; i < numTracks; i++)
    fV0motherPDG[i] = 0;
  
  fNumTagsStored = numTracks;
  
  
  
  
  
  // loop over V0 particles
  for (Int_t iv0 = 0; iv0 < event->GetNumberOfV0s(); iv0++) {
    AliESDv0* v0 = (AliESDv0*)event->GetV0(iv0);
 
    if (!v0)
      continue;
    
    // Reject onFly V0's <-> Only take V0's from offline V0 finder
    if (v0->GetOnFlyStatus())
      continue; 
  
    // Get the particle selection 
    Bool_t foundV0 = kFALSE;
    Int_t pdgV0 = 0, pdgP = 0, pdgN = 0;
    
    foundV0 = fV0KineCuts->ProcessV0(v0, pdgV0, pdgP, pdgN);
    if (!foundV0)
      continue;
    
    Int_t iTrackP = v0->GetPindex();  // positive track
    Int_t iTrackN = v0->GetNindex();  // negative track

    /*
    AliESDtrack* trackP = event->GetTrack(iTrackP);
    AliESDtrack* trackN = event->GetTrack(iTrackN);
    
    if (!trackP || !trackN)
      continue;
    
    
    
    Float_t xy = 999, z = 999;
    trackP->GetImpactParameters(xy, z);
    
    Bool_t reject = kFALSE;
    if (TMath::Abs(z) > 1)
      continue;
    
    trackN->GetImpactParameters(xy, z);
    if (TMath::Abs(z) > 1)
      continue;
    */
    
    
    // Fill the Object arrays
    // positive particles
    if (pdgP == -11) {
      fV0tags[iTrackP] = 14;
    }
    else if (pdgP == 211) {
      fV0tags[iTrackP] = 15;
    }
    else if(pdgP == 2212) {
      fV0tags[iTrackP] = 16;
    }
    
    fV0motherIndex[iTrackP] = iv0;
    fV0motherPDG[iTrackP] = pdgV0;

    // negative particles
    if( pdgN == 11){
      fV0tags[iTrackN] = -14;
    }
    else if( pdgN == -211){
      fV0tags[iTrackN] = -15;
    }
    else if( pdgN == -2212){
      fV0tags[iTrackN] = -16;
    }
    
    fV0motherIndex[iTrackN] = iv0;
    fV0motherPDG[iTrackN] = pdgV0;

  }
}


//______________________________________________________________________________
void AliTPCcalibResidualPID::ClearV0PIDlist()
{
  //
  // Clear the PID tag list
  //

  delete fV0tags;
  fV0tags = 0;
  
  delete fV0motherIndex;
  fV0motherIndex = 0;
  
  delete fV0motherPDG;
  fV0motherPDG = 0;
  
  fNumTagsStored = 0;
}


//______________________________________________________________________________
Char_t AliTPCcalibResidualPID::GetV0tag(Int_t trackIndex) const
{
  //
  // Get the tag for the corresponding trackIndex. Returns -99 in case of invalid index/tag list.
  //
  
  if (trackIndex < 0 || trackIndex >= fNumTagsStored || !fV0tags)
    return -99;
  
  return fV0tags[trackIndex];
}


//______________________________________________________________________________
Int_t AliTPCcalibResidualPID::GetV0motherIndex(Int_t trackIndex) const
{
  //
  // Get the index of the V0 mother for the corresponding trackIndex. Returns -99 in case of invalid index/mother index list.
  //
  
  if (trackIndex < 0 || trackIndex >= fNumTagsStored || !fV0motherIndex)
    return -99;
  
  return fV0motherIndex[trackIndex];
}


//______________________________________________________________________________
Int_t AliTPCcalibResidualPID::GetV0motherPDG(Int_t trackIndex) const
{
  //
  // Get the PDG code of the V0 mother for the corresponding trackIndex. Returns 0 in case of invalid index/mother index list.
  //
  
  if (trackIndex < 0 || trackIndex >= fNumTagsStored || !fV0motherPDG)
    return 0;
  
  return fV0motherPDG[trackIndex];
}


//______________________________________________________________________________
Int_t AliTPCcalibResidualPID::MergeGraphErrors(TGraphErrors* mergedGraph, TCollection* li) 
{
  // Adds all graphs with errors from the collection to this graph.
  // Returns the total number of poins in the result or -1 in case of an error.
  
  // TODO "Normal" merge of latest root trunk will also take into account the error bars,
  // so this function can be replaced by the normal root function with the latest root version.
  
  TIter next(li);
  while (TObject* o = next()) {
    TGraph *g = dynamic_cast<TGraph*>(o);
    if (!g) {
      Printf("ERROR: Cannot merge an object which doesn't inherit from TGraph found in the list");
      return -1;
    }
    int n0 = mergedGraph->GetN();
    int n1 = n0+g->GetN();
    mergedGraph->Set(n1);
    Double_t * x = g->GetX();
    Double_t * y = g->GetY();
    Double_t * ex = g->GetEX();
    Double_t * ey = g->GetEY();
    for (Int_t i = 0 ; i < g->GetN(); i++) {
      mergedGraph->SetPoint(n0+i, x[i], y[i]);
      Double_t exPoint = ex ? ex[i] : 0;
      Double_t eyPoint = ey ? ey[i] : 0;
      mergedGraph->SetPointError(n0+i, exPoint, eyPoint);
    }
  }
  return mergedGraph->GetN();
}

//________________________________________________________________________
Bool_t AliTPCcalibResidualPID::TPCCutMIGeo(const AliVTrack* track, const AliVEvent* evt, TTreeStream* streamer)
{
  //
  // TPC Cut MIGeo
  //

  if (!track || !evt)
    return kFALSE;
  
  const Short_t sign = track->Charge();
  Double_t xyz[50];
  Double_t pxpypz[50];
  Double_t cv[100];

  track->GetXYZ(xyz);
  track->GetPxPyPz(pxpypz);

  AliExternalTrackParam* par = new AliExternalTrackParam(xyz, pxpypz, cv, sign);
  const AliESDtrack dummy;

  const Double_t magField = evt->GetMagneticField();
  Double_t varGeom = dummy.GetLengthInActiveZone(par, 3, 236, magField, 0, 0);
  Double_t varNcr  = track->GetTPCClusterInfo(3, 1);
  Double_t varNcls = track->GetTPCsignalN();

  const Double_t varEval = 130. - 5. * TMath::Abs(1. / track->Pt());
  Bool_t cutGeom   = varGeom > fgCutGeo * varEval;
  Bool_t cutNcr    = varNcr  > fgCutNcr * varEval;
  Bool_t cutNcls   = varNcls > fgCutNcl * varEval;

  Bool_t kout = cutGeom && cutNcr && cutNcls;

  if (streamer) {
    Double_t dedx = track->GetTPCsignal();

    (*streamer)<<"tree"<<
      "param.="<< par<<
      "varGeom="<<varGeom<<
      "varNcr="<<varNcr<<
      "varNcls="<<varNcls<<
      
      "cutGeom="<<cutGeom<<
      "cutNcr="<<cutNcr<<
      "cutNcls="<<cutNcls<<
      
      "kout="<<kout<<
      "dedx="<<dedx<<
      
      "\n";
  }
  
  delete par;
  
  return kout;
}