PWGPP/TPC/AliTPCcalibResidualPID.cxx

   1 /**************************************************************************
   2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: Yvonne Pachmayer <pachmay@physi.uni-heidelberg.de>             *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15 //
  16 // The task:
  17 // stores TPC PID quantities in a THnSparse
  18 // output can then be used for e.g. dEdx calibration
  19 //
  20 // Author:
  21 // Yvonne Pachmayer <pachmay@physi.uni-heidelberg.de>
  22 //
  23
  24
  25 #include "AliTPCcalibResidualPID.h"
  26 #include "TChain.h"
  27 #include "AliAnalysisManager.h"
  28 #include "AliESDEvent.h"
  29 #include "AliAODEvent.h"
  30 #include "AliMCEvent.h"
  31 #include "AliMCParticle.h"
  32 //#include "AliStack.h"
  33 #include "AliPID.h"
  34 #include "AliTPCcalibResidualPID.h"
  35 #include "AliESDtrack.h"
  36 #include "AliESDtrackCuts.h"
  37 #include "AliESDpid.h"
  38 #include "AliESDInputHandler.h"
  39 #include "AliESDv0KineCuts.h"
  40 #include "AliESDv0.h"
  41 #include "AliCentrality.h"
  42 #include "THnSparse.h"
  43 #include "TH2D.h"
  44 #include "TCanvas.h"
  45 #include "TGraphErrors.h"
  46 #include "TMultiGraph.h"
  47 #include "TAxis.h"
  48 #include "TFile.h"
  49 #include "TSpline.h"
  50 #include "TStyle.h"
  51 #include "AliTPCdEdxInfo.h"
  52 #include "TString.h"
  53 #include "TFitResult.h"
  54 #include "TH1F.h"
  55 #include "TLegend.h"
  56 #include "TVirtualFitter.h"
  57
  58 using namespace std;
  59
  60 ClassImp(AliTPCcalibResidualPID)
  61
  62 Double_t AliTPCcalibResidualPID::fgCutGeo = 1.;
  63 Double_t AliTPCcalibResidualPID::fgCutNcr = 0.85;
  64 Double_t AliTPCcalibResidualPID::fgCutNcl = 0.7;
  65
  66 //________________________________________________________________________
  67 AliTPCcalibResidualPID::AliTPCcalibResidualPID()
  68   : AliAnalysisTaskSE(), fESD(0), fMC(0), fOutputContainer(0), fESDtrackCuts(0), fESDtrackCutsV0(0), fESDpid(0),
  69     fUseTPCCutMIGeo(kFALSE),
  70     fUseMCinfo(kTRUE),
  71     fIsPbpOrpPb(kFALSE),
  72     fZvtxCutEvent(9999.0),
  73     fV0KineCuts(0x0),
  74     fNumTagsStored(0),
  75     fV0tags(0x0),
  76     fV0motherIndex(0x0),
  77     fV0motherPDG(0x0),
  78     fProduceAllPadTypes(0),
  79     fProduceGlobal(0),
  80     fProduceShortPads(0),
  81     fProduceMediumPads(0),
  82     fProduceLongPads(0),
  83     fProduceOroc(0),
  84     fHistPidQA(0),
  85     fHistPidQAshort(0),
  86     fHistPidQAmedium(0),
  87     fHistPidQAlong(0),
  88     fHistPidQAoroc(0),
  89     fProduceTPCSignalSparse(0),
  90     fCorrectdEdxEtaDependence(0),
  91     fCorrectdEdxMultiplicityDependence(0),
  92     fThnspTpc(0),
  93     fQAList(0x0),
  94     fhInvMassGamma(0x0),
  95     fhInvMassK0s(0x0),
  96     fhInvMassLambda(0x0),
  97     fhInvMassAntiLambda(0x0),
  98     fhArmenterosAll(0x0),
  99     fhArmenterosGamma(0x0),
 100     fhArmenterosK0s(0x0),
 101     fhArmenterosLambda(0x0),
 102     fhArmenterosAntiLambda(0x0),
 103     fHistSharedClusQAV0Pi(0x0),
 104     fHistSharedClusQAV0Pr(0x0),
 105     fHistSharedClusQAV0El(0x0)
 106 {
 107   // default Constructor
 108    /* fast compilation test
 109      gSystem->Load("libANALYSIS");
 110      gSystem->Load("libANALYSISalice");
 111      .L /lustre/alice/akalweit/train/trunk/util/statsQA/AliTPCcalibResidualPID.cxx++
 112    */
 113
 114 }
 115
 116
 117 //________________________________________________________________________
 118 AliTPCcalibResidualPID::AliTPCcalibResidualPID(const char *name)
 119   : AliAnalysisTaskSE(name), fESD(0), fMC(0), fOutputContainer(0), fESDtrackCuts(0), fESDtrackCutsV0(0), fESDpid(0),
 120     fUseTPCCutMIGeo(kFALSE),
 121     fUseMCinfo(kTRUE),
 122     fIsPbpOrpPb(kFALSE),
 123     fZvtxCutEvent(9999.0),
 124     fV0KineCuts(0x0),
 125     fNumTagsStored(0),
 126     fV0tags(0x0),
 127     fV0motherIndex(0x0),
 128     fV0motherPDG(0x0),
 129     fProduceAllPadTypes(0),
 130     fProduceGlobal(0),
 131     fProduceShortPads(0),
 132     fProduceMediumPads(0),
 133     fProduceLongPads(0),
 134     fProduceOroc(0),
 135     fHistPidQA(0),
 136     fHistPidQAshort(0),
 137     fHistPidQAmedium(0),
 138     fHistPidQAlong(0),
 139     fHistPidQAoroc(0),
 140     fProduceTPCSignalSparse(0),
 141     fCorrectdEdxEtaDependence(0),
 142     fCorrectdEdxMultiplicityDependence(0),
 143     fThnspTpc(0),
 144     fQAList(0x0),
 145     fhInvMassGamma(0x0),
 146     fhInvMassK0s(0x0),
 147     fhInvMassLambda(0x0),
 148     fhInvMassAntiLambda(0x0),
 149     fhArmenterosAll(0x0),
 150     fhArmenterosGamma(0x0),
 151     fhArmenterosK0s(0x0),
 152     fhArmenterosLambda(0x0),
 153     fhArmenterosAntiLambda(0x0),
 154     fHistSharedClusQAV0Pi(0x0),
 155     fHistSharedClusQAV0Pr(0x0),
 156     fHistSharedClusQAV0El(0x0)
 157 {
 158
 159   //fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE);
 160   //fESDtrackCutsV0 = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kFALSE);
 161
 162   // Constructor
 163   DefineInput(0, TChain::Class());
 164   DefineOutput(1, TObjArray::Class());
 165   DefineOutput(2, TObjArray::Class());
 166
 167 }
 168
 169
 170 //_________________________________________________
 171 AliTPCcalibResidualPID::~AliTPCcalibResidualPID()
 172 {
 173   delete fOutputContainer;
 174   fOutputContainer = 0;
 175
 176   delete fQAList;
 177   fQAList = 0;
 178
 179   delete fESDtrackCuts;
 180   fESDtrackCuts = 0;
 181
 182   delete fESDtrackCutsV0;
 183   fESDtrackCutsV0 = 0;
 184
 185   delete fV0KineCuts;
 186   fV0KineCuts = 0;
 187
 188   delete fV0tags;
 189   fV0tags = 0;
 190   fNumTagsStored = 0;
 191
 192   delete fV0motherIndex;
 193   fV0motherIndex = 0;
 194
 195   delete fV0motherPDG;
 196   fV0motherPDG = 0;
 197 }
 198
 199
 200 //________________________________________________________________________
 201 void AliTPCcalibResidualPID::UserCreateOutputObjects()
 202 {
 203
 204     // THnSparse binning
 205     const Int_t kNdim = 9;
 206     //v0id,  dEdx,  TPCsigele,  TPCsigpion,  TOFbit,  eta, TPCclus, centr, p
 207     Int_t bins[kNdim] =    {    4,    250,    200,    200,          8,    20,        50,   20,   100};
 208     Double_t xmin[kNdim] = {  0.5,     30,   -10.,   -10.,       -1.5,   -1.,       60.,   0.,   0.1};
 209     Double_t xmax[kNdim] = {  4.5,   500.,    10.,    10.,        6.5,    1.,      160.,  100,   4};
 210     fThnspTpc= new THnSparseF("tpcsignals", "TPC signal;v0id;tpc signal;tpcsigele;tpcsigpion;tofbit;eta;tpcclus;centr;p (GeV/c)", kNdim, bins, xmin, xmax);
 211     BinLogAxis(fThnspTpc, 8);
 212
 213     //
 214     // 0.ptot, 1.tpcSig, 2.particle ID, 3. assumed particle, 4. nSigmaTPC (4x), 5. nSigmaTOF (4x), 6. centrality
 215     // Concerning 2 (particle ID):
 216     // - in case of MC: Contains MC ID. Bin 1-4 (<=> Slot 0-3): el, pi, ka, pr
 217     // - 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
 218     //
 219
 220     Int_t    binsHistQA[7] = {135, 1980,    4,    5, 40, 10,   40};
 221     Double_t xminHistQA[7] = {0.1,   20, -0.5, -0.5, -10, -5,   0.};
 222     Double_t xmaxHistQA[7] = {50., 2000,  3.5,  4.5,  10,  5, 20000};
 223     fHistPidQA = new THnSparseF("fHistPidQA","PID QA",7,binsHistQA,xminHistQA,xmaxHistQA);
 224     BinLogAxis(fHistPidQA, 0);
 225
 226     //
 227     fHistPidQAshort  = new THnSparseF("fHistPidQAshort" ,"PID QA -- short pads",7,binsHistQA,xminHistQA,xmaxHistQA);
 228     fHistPidQAmedium = new THnSparseF("fHistPidQAmedium","PID QA -- med pads",7,binsHistQA,xminHistQA,xmaxHistQA);
 229     fHistPidQAlong   = new THnSparseF("fHistPidQAlong"  ,"PID QA -- long pads",7,binsHistQA,xminHistQA,xmaxHistQA);
 230     fHistPidQAoroc   = new THnSparseF("fHistPidQAoroc"  ,"PID QA -- oroc",7,binsHistQA,xminHistQA,xmaxHistQA);
 231     BinLogAxis(fHistPidQAshort, 0);
 232     BinLogAxis(fHistPidQAmedium, 0);
 233     BinLogAxis(fHistPidQAlong, 0);
 234     BinLogAxis(fHistPidQAoroc, 0);
 235     //
 236     fOutputContainer = new TObjArray(2);
 237     fOutputContainer->SetName(GetName());
 238     fOutputContainer->SetOwner();
 239
 240     if(fProduceTPCSignalSparse)fOutputContainer->Add(fThnspTpc);
 241     if(fProduceGlobal)fOutputContainer->Add(fHistPidQA);
 242     if(fProduceAllPadTypes || fProduceShortPads)fOutputContainer->Add(fHistPidQAshort);
 243     if(fProduceAllPadTypes || fProduceMediumPads)fOutputContainer->Add(fHistPidQAmedium);
 244     if(fProduceAllPadTypes || fProduceLongPads)fOutputContainer->Add(fHistPidQAlong);
 245     if(fProduceAllPadTypes || fProduceOroc)fOutputContainer->Add(fHistPidQAoroc);
 246
 247     // V0 Kine cuts
 248     fV0KineCuts = new AliESDv0KineCuts;
 249     fV0KineCuts->SetGammaCutChi2NDF(5.);
 250     // Only accept V0el with prod. radius within 45 cm -> PID will by systematically biased for larger values!
 251     Float_t gammaProdVertexRadiusCuts[2] = { 3.0, 45. };
 252     fV0KineCuts->SetGammaCutVertexR(&gammaProdVertexRadiusCuts[0]);
 253
 254
 255     fQAList = new TObjArray(4);
 256     fQAList->SetName(GetName());
 257     fQAList->SetOwner();
 258
 259     fhInvMassGamma      = new TH1F("fhInvMassGamma", "Invariant Mass of gammas; m_{ee} (GeV/#it{c}^{2}); Entries", 200, 0., 0.2);
 260     fhInvMassK0s        = new TH1F("fhInvMassK0s", "Invariant Mass of K0s; m_{#pi#pi} (GeV/#it{c}^{2}); Entries;", 200, 0.45, 0.55);
 261     fhInvMassLambda     = new TH1F("fhInvMassLambda", "Invariant Mass of lambdas; m_{p#pi^{-}} (GeV/#it{c}^{2}); Entries", 200, 1.05, 1.15);
 262     fhInvMassAntiLambda = new TH1F("fhInvMassAntiLambda", "Invariant Mass of anti-lambdas; m_{#pi^{+}#bar{p}} (GeV/#it{c}^{2}); Entries",
 263                                    200, 1.05, 1.15);
 264
 265     fQAList->Add(fhInvMassGamma);
 266     fQAList->Add(fhInvMassK0s);
 267     fQAList->Add(fhInvMassLambda);
 268     fQAList->Add(fhInvMassAntiLambda);
 269
 270
 271     fhArmenterosAll = new TH2F("fhArmenterosAll",
 272                                "Armenteros plot all V0s;#alpha = (#it{p}^{+}_{L}-#it{p}^{-}_{L})/(#it{p}^{+}_{L}+#it{p}^{-}_{L});#it{q}_{T} (GeV/#it{c})",
 273                                200, -1., 1., 200, 0., 0.4);
 274     fhArmenterosGamma = new TH2F("fhArmenterosGamma",
 275                                  "Armenteros plot Gamma;#alpha = (#it{p}^{+}_{L}-#it{p}^{-}_{L})/(#it{p}^{+}_{L}+#it{p}^{-}_{L});#it{q}_{T} (GeV/#it{c})",
 276                                  200, -1., 1., 200, 0., 0.4);
 277     fhArmenterosK0s = new TH2F("fhArmenterosK0s",
 278                                "Armenteros plot K0s;#alpha = (#it{p}^{+}_{L}-#it{p}^{-}_{L})/(#it{p}^{+}_{L}+#it{p}^{-}_{L});#it{q}_{T} (GeV/#it{c})",
 279                                200, -1., 1., 200, 0., 0.4);
 280     fhArmenterosLambda = new TH2F("fhArmenterosLambda",
 281                                  "Armenteros plot lambda;#alpha = (#it{p}^{+}_{L}-#it{p}^{-}_{L})/(#it{p}^{+}_{L}+#it{p}^{-}_{L});#it{q}_{T} (GeV/#it{c})",
 282                                  200, -1., 1., 200, 0., 0.4);
 283     fhArmenterosAntiLambda = new TH2F("fhArmenterosAntiLambda",
 284                                       "Armenteros plot anti-lambda;#alpha = (#it{p}^{+}_{L}-#it{p}^{-}_{L})/(#it{p}^{+}_{L}+#it{p}^{-}_{L});#it{q}_{T} (GeV/#it{c})",
 285                                       200, -1., 1., 200, 0., 0.4);
 286
 287     fQAList->Add(fhArmenterosAll);
 288     fQAList->Add(fhArmenterosGamma);
 289     fQAList->Add(fhArmenterosK0s);
 290     fQAList->Add(fhArmenterosLambda);
 291     fQAList->Add(fhArmenterosAntiLambda);
 292
 293
 294     const Int_t dimQASharedClusters = 4;
 295     const TString axisTitles[dimQASharedClusters] = { "#it{p}_{TPC} (GeV/#it{c})", "#it{#Delta'}", "#it{N}_{shared cl}", "Pad row"};
 296     Int_t    binsHistQASharedClusters[dimQASharedClusters] = { 100,  100, 160, 160};
 297     Double_t xminHistQASharedClusters[dimQASharedClusters] = { 0.3, 0.5,  0,   -1};
 298     Double_t xmaxHistQASharedClusters[dimQASharedClusters] = { 20,  1.5, 160, 159};
 299
 300     fHistSharedClusQAV0Pi = new THnSparseF("fHistSharedClusQAV0Pi" ,"PID QA shared clusters - V0 pi", dimQASharedClusters,
 301                                            binsHistQASharedClusters, xminHistQASharedClusters, xmaxHistQASharedClusters);
 302     BinLogAxis(fHistSharedClusQAV0Pi, 0);
 303     for (Int_t i = 0; i < dimQASharedClusters; i++)
 304       fHistSharedClusQAV0Pi->GetAxis(i)->SetTitle(axisTitles[i].Data());
 305     fQAList->Add(fHistSharedClusQAV0Pi);
 306
 307     fHistSharedClusQAV0Pr = new THnSparseF("fHistSharedClusQAV0Pr" ,"PID QA shared clusters - V0 pr", dimQASharedClusters,
 308                                            binsHistQASharedClusters, xminHistQASharedClusters, xmaxHistQASharedClusters);
 309     BinLogAxis(fHistSharedClusQAV0Pr, 0);
 310     for (Int_t i = 0; i < dimQASharedClusters; i++)
 311       fHistSharedClusQAV0Pi->GetAxis(i)->SetTitle(axisTitles[i].Data());
 312     fQAList->Add(fHistSharedClusQAV0Pr);
 313
 314     fHistSharedClusQAV0El = new THnSparseF("fHistSharedClusQAV0El" ,"PID QA shared clusters - V0 el", dimQASharedClusters,
 315                                            binsHistQASharedClusters, xminHistQASharedClusters, xmaxHistQASharedClusters);
 316     BinLogAxis(fHistSharedClusQAV0El, 0);
 317     for (Int_t i = 0; i < dimQASharedClusters; i++)
 318       fHistSharedClusQAV0Pi->GetAxis(i)->SetTitle(axisTitles[i].Data());
 319     fQAList->Add(fHistSharedClusQAV0El);
 320
 321     PostData(1,fOutputContainer);
 322     PostData(2,fQAList);
 323 }
 324
 325
 326 //_____________________________________________________________________________
 327 void AliTPCcalibResidualPID::UserExec(Option_t *)
 328 {
 329     //calls the Process function
 330     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 331
 332     if (!esdH) {
 333       printf("ERROR: Could not get ESDInputHandler \n");
 334     }
 335     else fESD = esdH->GetEvent();
 336
 337     // If MC, forward MCevent
 338     fMC = dynamic_cast<AliMCEvent*>(MCEvent());
 339     //
 340     Process(fESD, fMC);
 341     //
 342     PostData(1,fOutputContainer);
 343     PostData(2,fQAList);
 344 }
 345
 346
 347
 348 //________________________________________________________________________
 349 void AliTPCcalibResidualPID::Process(AliESDEvent *const esdEvent, AliMCEvent *const mcEvent)
 350 {
 351
 352   //called for each event
 353   if (!esdEvent) {
 354     Printf("ERROR: esdEvent not available");
 355     return;
 356   }
 357
 358   if(!((AliESDInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))) printf("ESD inputhandler not available \n");
 359
 360   if (!fESDpid) {
 361     fESDpid = ((AliESDInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetESDpid();
 362   }
 363
 364   if (!fESDpid || !fV0KineCuts)
 365     return;
 366
 367
 368   Float_t centralityFper=99;
 369
 370   AliCentrality *esdCentrality = esdEvent->GetCentrality();
 371   centralityFper = esdCentrality->GetCentralityPercentile("V0M");
 372
 373   if (!GetVertexIsOk(esdEvent))
 374     return;
 375
 376   const AliESDVertex* fESDvertex = esdEvent->GetPrimaryVertexTracks();
 377   if (!fESDvertex)
 378     return;
 379
 380   Int_t ncontr = fESDvertex->GetNContributors();
 381
 382   if (ncontr <= 0)
 383     return;
 384
 385   // Fill V0 arrays with V0s
 386   FillV0PIDlist(esdEvent);
 387
 388   // Array with flags wheter QA for this V0 was already done or not
 389   const Int_t numV0s = esdEvent->GetNumberOfV0s();
 390   Bool_t v0QAadded[numV0s];
 391   for (Int_t i = 0; i < numV0s; i++)
 392     v0QAadded[i] = kFALSE;
 393
 394   Int_t nTotTracks = esdEvent->GetNumberOfTracks();
 395   for (Int_t iTracks = 0; iTracks < nTotTracks; iTracks++){//begin track loop
 396     AliESDtrack *trackESD = esdEvent->GetTrack(iTracks);
 397     if(!trackESD) {
 398       Printf("ERROR: Could not receive track %d (esd loop)", iTracks);
 399       continue;
 400     }
 401     if((TMath::Abs(trackESD->Eta())) > 0.9)
 402       continue;
 403
 404      // Do not distinguish positively and negatively charged V0's
 405     Char_t v0tagAllCharges = TMath::Abs(GetV0tag(iTracks));
 406     if (v0tagAllCharges == -99) {
 407       AliError(Form("Problem with V0 tag list (requested V0 track for track %d from %d (list status %d))!", iTracks, esdEvent->GetNumberOfTracks(),
 408                     fV0tags != 0x0));
 409       continue;
 410     }
 411
 412     Bool_t isV0el = v0tagAllCharges == 14;
 413     Bool_t isV0pi = v0tagAllCharges == 15;
 414     Bool_t isV0pr = v0tagAllCharges == 16;
 415     Bool_t isV0 = isV0el || isV0pi || isV0pr;
 416
 417     if (mcEvent && fUseMCinfo) {
 418       // For MC, do not look for V0s, i.e. demand the non-V0 track cuts
 419       if (fESDtrackCuts && !fESDtrackCuts->AcceptTrack(trackESD)) continue;
 420
 421       if (fUseTPCCutMIGeo) {
 422         // If cut on geometry is active, don't cut on number of clusters, since such a cut is already included.
 423         if (!TPCCutMIGeo(trackESD, esdEvent))
 424           continue;
 425       }
 426       else {
 427         // If cut on geometry is not active, always cut on num clusters
 428         if (trackESD->GetTPCsignalN() < 60)
 429           continue;
 430       }
 431     }
 432     else {
 433       // For data, take V0 AND non-V0 with separate cuts
 434       //if (!isV0 && fESDtrackCuts && !fESDtrackCuts->AcceptTrack(trackESD)) continue;
 435       if (isV0) {
 436         if (fESDtrackCutsV0 && !fESDtrackCutsV0->AcceptTrack(trackESD))
 437           continue;
 438       }
 439       else {
 440         if (fESDtrackCuts && !fESDtrackCuts->AcceptTrack(trackESD))
 441           continue;
 442       }
 443
 444       if (fUseTPCCutMIGeo) {
 445         // If cut on geometry is active, don't cut on number of clusters, since such a cut is already included.
 446         if (!isV0) {
 447           if (!TPCCutMIGeo(trackESD, esdEvent))
 448             continue;
 449         }
 450         else {
 451           // One should not cut on geometry for V0's since they have different topology. (Loosely) cut on num of clusters instead.
 452           if (trackESD->GetTPCsignalN() < 60)
 453             continue;
 454         }
 455       }
 456       else {
 457         // If cut on geometry is not active, always cut on num clusters
 458         if (trackESD->GetTPCsignalN() < 60)
 459           continue;
 460       }
 461     }
 462
 463     const AliExternalTrackParam *paramIn = trackESD->GetInnerParam();
 464     Float_t precin=-1;
 465     if(paramIn) {
 466       precin=paramIn->GetP();
 467     } else {
 468       continue;
 469     }
 470     Int_t precdefault=CompareFloat(precin,-1);
 471     if(precdefault==1) continue; // momentum cut
 472     //
 473
 474     AliMCParticle* mcTrack = 0x0;
 475     Int_t particleID = -1;
 476     Int_t v0id = 0;
 477
 478     if (mcEvent && fUseMCinfo) {
 479       // MC - particle ID = MC ID
 480       Int_t label = trackESD->GetLabel();
 481
 482       if (label < 0)
 483         continue; // If MC is available, reject tracks with negative ESD label
 484
 485       mcTrack = dynamic_cast<AliMCParticle*>(mcEvent->GetTrack(TMath::Abs(label)));
 486       if (!mcTrack) {
 487         Printf("ERROR: Could not receive mcTrack with label %d for track %d", label, iTracks);
 488         continue;
 489       }
 490
 491       /*// Only accept MC primaries
 492       if (!mcEvent->Stack()->IsPhysicalPrimary(TMath::Abs(label))) {
 493         continue;
 494       }*/
 495
 496       Int_t pdgAbs = TMath::Abs(mcTrack->PdgCode());
 497
 498       if (pdgAbs == 11) { // electron
 499         particleID = 0;
 500       }
 501       else if (pdgAbs == 211) { // pion
 502         particleID = 1;
 503       }
 504       else if (pdgAbs == 321) { // kaon
 505         particleID = 2;
 506       }
 507       else if (pdgAbs == 2212) { // proton
 508         particleID = 3;
 509       }
 510       else
 511         continue; // Reject all other particles
 512     }
 513     else {
 514       // Data - particle ID = V0 ID (if V0)
 515       if (isV0pi) { // pion
 516         particleID = 2;
 517         v0id = 2;
 518       }
 519       else if (isV0el) { // electron
 520         particleID = 1;
 521         v0id = 1;
 522       }
 523       else if (isV0pr) { // proton
 524         particleID = 3;
 525         v0id = 3;
 526       }
 527       else { // No V0-ID available (species must be determined via TPC, TOF, ...)
 528         particleID = 0;
 529       }
 530     }
 531
 532
 533     Float_t tpcsignal=trackESD->GetTPCsignal();
 534     Int_t tofbit=-1;
 535     Int_t iTOFpid=0;
 536     Int_t ikTIME=0;
 537     Double_t tpcnsigmaele=-10;
 538     Double_t tpcnsigmapion=-10;
 539     //
 540     if ((trackESD->GetStatus() & AliESDtrack::kTOFpid)) iTOFpid = 1;
 541     if ((trackESD->GetStatus() & AliESDtrack::kTIME)) ikTIME = 1;
 542     Float_t time0 = fESDpid->GetTOFResponse().GetTimeZero();
 543     //
 544     if((iTOFpid==1) &&(ikTIME==1)){
 545       Double_t tofnsigmaele= fESDpid->NumberOfSigmasTOF(trackESD,AliPID::kElectron, time0);
 546       Double_t tofnsigmapion=fESDpid->NumberOfSigmasTOF(trackESD,AliPID::kPion, time0);
 547       if (TMath::Abs(tofnsigmapion)<3) tofbit = 1;
 548       if (TMath::Abs(tofnsigmaele)<3)  tofbit = 0;
 549       tpcnsigmaele=fESDpid->NumberOfSigmasTPC(trackESD,AliPID::kElectron);
 550       tpcnsigmapion=fESDpid->NumberOfSigmasTPC(trackESD,AliPID::kPion);
 551     }
 552     //
 553     Int_t tpcnclusN=trackESD->GetTPCsignalN();
 554     Double_t eta=trackESD->Eta();
 555
 556     //
 557     if(fProduceTPCSignalSparse){
 558       Double_t contentSignal[9];
 559       contentSignal[0]=v0id;
 560       contentSignal[1]=tpcsignal;
 561       contentSignal[2]=tpcnsigmaele;
 562       contentSignal[3]=tpcnsigmapion;
 563       contentSignal[4]=tofbit;
 564       contentSignal[5]=eta;
 565       contentSignal[6]=tpcnclusN;
 566       contentSignal[7]=centralityFper;
 567       contentSignal[8]=precin;
 568       //
 569       if (fThnspTpc->GetEntries() < 1e6) fThnspTpc->Fill(contentSignal);
 570     }//should it be created or not
 571
 572
 573
 574     //
 575     // 2nd THnSparse
 576     //
 577     Double_t tpcQA[5] = {fESDpid->NumberOfSigmasTPC(trackESD, AliPID::kElectron),
 578                          fESDpid->NumberOfSigmasTPC(trackESD, AliPID::kPion),
 579                          fESDpid->NumberOfSigmasTPC(trackESD, AliPID::kKaon),
 580                          fESDpid->NumberOfSigmasTPC(trackESD, AliPID::kProton),
 581                          0};
 582     Double_t tofQA[5] = {fESDpid->NumberOfSigmasTOF(trackESD, AliPID::kElectron, time0),
 583                          fESDpid->NumberOfSigmasTOF(trackESD, AliPID::kPion, time0),
 584                          fESDpid->NumberOfSigmasTOF(trackESD, AliPID::kKaon, time0),
 585                          fESDpid->NumberOfSigmasTOF(trackESD, AliPID::kProton, time0),
 586                          0 };
 587
 588     // id 5 is just again Kaons in restricted eta range
 589     tpcQA[4] = tpcQA[2];
 590     tofQA[4] = tofQA[2];
 591
 592     //
 593     // dE/dx in different pad regions
 594     //
 595     AliTPCdEdxInfo * infoTpcPid = trackESD->GetTPCdEdxInfo();
 596     Double32_t signal[4]; Char_t ncl[3]; Char_t nrows[3];
 597     if (infoTpcPid) {
 598       infoTpcPid->GetTPCSignalRegionInfo(signal, ncl, nrows);
 599     } else {
 600       for(Int_t iarr = 0; iarr < 3; iarr++) {
 601         signal[iarr] = 0;
 602         ncl[iarr] = 0;
 603         nrows[iarr] = 0;
 604       }
 605       signal[3] = 0;
 606     }
 607     //
 608     UInt_t status = trackESD->GetStatus();
 609     Bool_t hasTOFout  = status&AliESDtrack::kTOFout;
 610     Bool_t hasTOFtime = status&AliESDtrack::kTIME;
 611     Bool_t hasTOF     = kFALSE;
 612     if (hasTOFout && hasTOFtime) hasTOF = kTRUE;
 613     Float_t length = trackESD->GetIntegratedLength();
 614     // Length check only for primaries!
 615     if (length < 350 && !isV0) hasTOF = kFALSE;
 616     //
 617
 618     if (!hasTOF) {
 619       // Make sure that number of sigmas is large in this case, so that the track will be rejected if a TOF cut is applied
 620       for (Int_t i = 0; i < 5; i++) {
 621         tofQA[i] = 999;
 622       }
 623     }
 624
 625
 626     Double_t processedTPCsignal[5] = { tpcsignal, tpcsignal, tpcsignal, tpcsignal, tpcsignal };
 627
 628
 629     if (fCorrectdEdxEtaDependence && fCorrectdEdxMultiplicityDependence) {
 630       processedTPCsignal[0] = fESDpid->GetTPCResponse().GetEtaAndMultiplicityCorrectedTrackdEdx(trackESD, AliPID::kElectron,
 631                                                                                                 AliTPCPIDResponse::kdEdxDefault);
 632       processedTPCsignal[1] = fESDpid->GetTPCResponse().GetEtaAndMultiplicityCorrectedTrackdEdx(trackESD, AliPID::kPion,
 633                                                                                                 AliTPCPIDResponse::kdEdxDefault);
 634       processedTPCsignal[2] = fESDpid->GetTPCResponse().GetEtaAndMultiplicityCorrectedTrackdEdx(trackESD, AliPID::kKaon,
 635                                                                                                 AliTPCPIDResponse::kdEdxDefault);
 636       processedTPCsignal[3] = fESDpid->GetTPCResponse().GetEtaAndMultiplicityCorrectedTrackdEdx(trackESD, AliPID::kProton,
 637                                                                                                 AliTPCPIDResponse::kdEdxDefault);
 638     }
 639     else if (fCorrectdEdxEtaDependence) {
 640       processedTPCsignal[0] = fESDpid->GetTPCResponse().GetEtaCorrectedTrackdEdx(trackESD, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault);
 641       processedTPCsignal[1] = fESDpid->GetTPCResponse().GetEtaCorrectedTrackdEdx(trackESD, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault);
 642       processedTPCsignal[2] = fESDpid->GetTPCResponse().GetEtaCorrectedTrackdEdx(trackESD, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault);
 643       processedTPCsignal[3] = fESDpid->GetTPCResponse().GetEtaCorrectedTrackdEdx(trackESD, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault);
 644     }
 645     else if (fCorrectdEdxMultiplicityDependence) {
 646       processedTPCsignal[0] = fESDpid->GetTPCResponse().GetMultiplicityCorrectedTrackdEdx(trackESD, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault);
 647       processedTPCsignal[1] = fESDpid->GetTPCResponse().GetMultiplicityCorrectedTrackdEdx(trackESD, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault);
 648       processedTPCsignal[2] = fESDpid->GetTPCResponse().GetMultiplicityCorrectedTrackdEdx(trackESD, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault);
 649       processedTPCsignal[3] = fESDpid->GetTPCResponse().GetMultiplicityCorrectedTrackdEdx(trackESD, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault);
 650     }
 651
 652     // id 5 is just again Kaons in restricted eta range
 653     processedTPCsignal[4] = processedTPCsignal[2];
 654
 655     for(Int_t iPart = 0; iPart < 5; iPart++) {
 656       // 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)
 657       if (iPart == 4 && ((mcEvent && fUseMCinfo) || abs(trackESD->Eta()) > 0.2))
 658         continue;
 659
 660       Double_t vecHistQA[7] = {precin, processedTPCsignal[iPart], particleID, iPart, tpcQA[iPart], tofQA[iPart], nTotTracks};
 661       if (fProduceGlobal) fHistPidQA->Fill(vecHistQA);
 662       vecHistQA[1] = signal[0]; vecHistQA[4] = ncl[0];
 663       if ((fProduceAllPadTypes || fProduceShortPads)) fHistPidQAshort->Fill(vecHistQA);
 664       //
 665       vecHistQA[1] = signal[1]; vecHistQA[4] = ncl[1];
 666       if ((fProduceAllPadTypes || fProduceMediumPads)) fHistPidQAmedium->Fill(vecHistQA);
 667       //
 668       vecHistQA[1] = signal[2]; vecHistQA[4] = ncl[2];
 669       if ((fProduceAllPadTypes || fProduceLongPads)) fHistPidQAlong->Fill(vecHistQA);
 670       //
 671       vecHistQA[1] = signal[3]; vecHistQA[4] = nrows[1] + nrows[2];
 672       if ((fProduceAllPadTypes || fProduceOroc)) fHistPidQAoroc->Fill(vecHistQA);
 673     }
 674
 675     if (isV0) {
 676       // Fill QA hists for shared clusters dE/dx
 677       Double_t vecHistQA[4] = { precin, -1, trackESD->GetTPCSharedMap().CountBits(), -1 };
 678       THnSparse* currHist = fHistSharedClusQAV0Pi;
 679
 680       if (isV0pi) {
 681         Double_t expectedDeDx = fESDpid->GetTPCResponse().GetExpectedSignal(trackESD, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault,
 682                                                                             fCorrectdEdxEtaDependence, fCorrectdEdxMultiplicityDependence);
 683         if (expectedDeDx > 0 ) {
 684           vecHistQA[1] = tpcsignal / expectedDeDx;
 685           currHist = fHistSharedClusQAV0Pi;
 686         }
 687       }
 688       else if (isV0pr) {
 689         Double_t expectedDeDx = fESDpid->GetTPCResponse().GetExpectedSignal(trackESD, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault,
 690                                                                             fCorrectdEdxEtaDependence, fCorrectdEdxMultiplicityDependence);
 691         if (expectedDeDx > 0 ) {
 692           vecHistQA[1] = tpcsignal / expectedDeDx;
 693           currHist = fHistSharedClusQAV0Pr;
 694         }
 695       }
 696       else if (isV0el) {
 697         Double_t expectedDeDx = fESDpid->GetTPCResponse().GetExpectedSignal(trackESD, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault,
 698                                                                             fCorrectdEdxEtaDependence, fCorrectdEdxMultiplicityDependence);
 699         if (expectedDeDx > 0 ) {
 700           vecHistQA[1] = tpcsignal / expectedDeDx;
 701           currHist = fHistSharedClusQAV0El;
 702         }
 703       }
 704
 705       Int_t iRowInd = -1;
 706       // iRowInd == -1 for "all rows w/o multiple counting"
 707       currHist->Fill(vecHistQA);
 708
 709       // Fill hist for every pad row with shared cluster
 710       for (iRowInd = 0; iRowInd < 159; iRowInd++) {
 711         if (trackESD->GetTPCSharedMap().TestBitNumber(iRowInd)) {
 712           vecHistQA[3] = iRowInd;
 713           currHist->Fill(vecHistQA);
 714         }
 715       }
 716
 717
 718
 719       // Check, whether the QA for this V0 mother was already filled into the histograms (is the case if the other daughter has already
 720       // been processed). If not, set flag to kTRUE and fill the QA.
 721       const Int_t iV0 = GetV0motherIndex(iTracks);
 722       if (!v0QAadded[iV0]) {
 723         v0QAadded[iV0] = kTRUE;
 724
 725         AliESDv0* esdV0 = (AliESDv0*)esdEvent->GetV0(iV0);
 726
 727         if (!esdV0) {
 728           printf("Error: V0 tagged, but does not exist!\n");
 729           continue;
 730         }
 731
 732         Float_t armVar[2] = {0.0,0.0};
 733         fV0KineCuts->Armenteros(esdV0, armVar);
 734
 735         const Int_t motherPDG = GetV0motherPDG(iTracks);
 736
 737         // Mother and daughter match the requirements, otherwise it wouldn't be tagged as a V0. Now just fill the QA histos.
 738         fhArmenterosAll->Fill(armVar[0], armVar[1]);
 739         //if (esdV0->TestBit(BIT(14))) {
 740         if (TMath::Abs(motherPDG) == 22) {
 741           fhInvMassGamma->Fill(esdV0->GetEffMass(AliPID::kElectron, AliPID::kElectron));
 742           fhArmenterosGamma->Fill(armVar[0], armVar[1]);
 743         }
 744         else if (TMath::Abs(motherPDG) == 310) {
 745         //else if (esdV0->TestBit(BIT(15))) {
 746           fhInvMassK0s->Fill(esdV0->GetEffMass(AliPID::kPion, AliPID::kPion));
 747           fhArmenterosK0s->Fill(armVar[0], armVar[1]);
 748         }
 749         else if (motherPDG == 3122) {
 750         //else if (esdV0->TestBit(BIT(16))) {
 751           fhInvMassLambda->Fill(esdV0->GetEffMass(AliPID::kProton, AliPID::kPion));
 752           fhArmenterosLambda->Fill(armVar[0], armVar[1]);
 753         }
 754         else if (motherPDG == -3122) {
 755         //else if (esdV0->TestBit(BIT(17))) {
 756           fhInvMassAntiLambda->Fill(esdV0->GetEffMass(AliPID::kPion, AliPID::kProton));
 757           fhArmenterosAntiLambda->Fill(armVar[0], armVar[1]);
 758         }
 759       }
 760     }
 761   } //end track loop
 762
 763
 764   // Clear the V0 PID arrays
 765   ClearV0PIDlist();
 766 }
 767
 768
 769 //________________________________________________________________________
 770 Int_t  AliTPCcalibResidualPID::CompareFloat(Float_t f1, Float_t f2) const
 771 {
 772     //compares if the Float_t f1 is equal with f2 and returns 1 if true and 0 if false
 773     Float_t precision = 0.00001;
 774     if (((f1 - precision) < f2) &&
 775         ((f1 + precision) > f2))
 776     {
 777         return 1;
 778     }
 779     else
 780     {
 781         return 0;
 782     }
 783 }
 784
 785
 786 //________________________________________________________________________
 787 void AliTPCcalibResidualPID::Terminate(const Option_t *)
 788 {
 789
 790
 791 }
 792
 793
 794
 795 //________________________________________________________________________
 796 void AliTPCcalibResidualPID::BinLogAxis(const THnSparseF *h, Int_t axisNumber) {
 797   //
 798   // Method for the correct logarithmic binning of histograms
 799   //
 800   TAxis *axis = h->GetAxis(axisNumber);
 801   int bins = axis->GetNbins();
 802
 803   Double_t from = axis->GetXmin();
 804   Double_t to = axis->GetXmax();
 805   Double_t *newBins = new Double_t[bins + 1];
 806
 807   newBins[0] = from;
 808   Double_t factor = pow(to/from, 1./bins);
 809
 810   for (int i = 1; i <= bins; i++) {
 811    newBins[i] = factor * newBins[i-1];
 812   }
 813   axis->Set(bins, newBins);
 814   delete [] newBins;
 815
 816 }
 817
 818
 819 //________________________________________________________________________
 820 Double_t* AliTPCcalibResidualPID::ExtractResidualPID(THnSparseF * histPidQA, const Bool_t useV0s, const Char_t * outFile,
 821                                                      const Char_t * type, const Char_t * period, const Char_t * pass,
 822                                                      const Char_t * system, const Double_t * initialParameters,
 823                                                      const Char_t *dedxtype,
 824                                                      AliTPCcalibResidualPID::FitType fitType) {
 825   //
 826   // (1.) obtain residual graphs
 827   //
 828   TObjArray * arr = 0x0;
 829
 830   Bool_t isMC = kFALSE;
 831   if (strcmp(type, "MC") == 0) {
 832     arr = GetResidualGraphsMC(histPidQA, system);
 833     isMC = kTRUE;
 834   }
 835   else if (strcmp(type, "DATA") == 0) {
 836     arr = GetResidualGraphs(histPidQA, system, useV0s);
 837   }
 838   else {
 839     Printf("ERROR - ExtractResidualPID: Unknown type \"%s\" - must be \"MC\" or \"DATA\"!", type);
 840
 841     return 0x0;
 842   }
 843
 844   Bool_t isPPb = kFALSE;
 845   if (strcmp(system, "PPB") == 0 || strcmp(system, "PBP") == 0) {
 846     Printf("p-Pb/Pb-p detected - Applying special handling!");
 847     isPPb = kTRUE;
 848   }
 849   //
 850   // (2.) get old-style Bethe-Bloch parameters
 851   //
 852   TF1* parametrisation = FitBB(arr, isMC, isPPb, useV0s, initialParameters, fitType);
 853   //
 854   // (3.) obtain response functions to OADB
 855   //
 856   TObjArray* arrResponse = GetResponseFunctions(parametrisation, arr, type, period, pass, system, dedxtype);
 857   //
 858   // (4.) write results to file and exit
 859   //
 860   if (arrResponse) {
 861     TFile outputFile(outFile,"RECREATE");
 862     arrResponse->Write();
 863     outputFile.Close();
 864   }
 865   //
 866
 867   return parametrisation->GetParameters();
 868 }
 869
 870
 871 //________________________________________________________________________
 872 TObjArray * AliTPCcalibResidualPID::GetSeparation(THnSparseF * histPidQA, Int_t kParticle1, Int_t kParticle2) {
 873
 874   //get THnSparse
 875   THnSparse * hist = histPidQA;
 876
 877
 878   Float_t nSigmaPlus1 = 0, nSigmaMinus1 = 0; //sigma of first particle in the TPC
 879   Float_t nSigmaPlus2 = 0, nSigmaMinus2 = 0; //sigma of second particle in the TPC
 880
 881   if(kParticle1 == 0){ // electron
 882         nSigmaMinus1   =  -1.999;
 883     nSigmaPlus1   =  2.999;
 884   } else if(kParticle1 == 1){ //pion
 885         nSigmaMinus1   =  -1.999;
 886     nSigmaPlus1   =  2.999;
 887   } else if(kParticle1 == 2){ //kaons
 888         nSigmaMinus1   =  -2.999;
 889     nSigmaPlus1   =    2.999;
 890   } else if(kParticle1 == 3){ //protons
 891         nSigmaMinus1   =  -2.999;
 892     nSigmaPlus1   =    2.999;
 893   }
 894
 895   //
 896   // 1. select and fit first particle
 897   //
 898   hist->GetAxis(3)->SetRangeUser(kParticle1,kParticle1);
 899   hist->GetAxis(4)->SetRangeUser(nSigmaMinus1,nSigmaPlus1);
 900   hist->GetAxis(5)->SetRangeUser(-2.999,2.999); // 3sigma in TOF
 901   TH2D * histPart1 = (TH2D*) hist->Projection(1,0)->Clone("histPart1");
 902   histPart1->SetTitle(";p /(GeV/c) ; TPCSignal /(a.u.)");
 903   histPart1->RebinX(4);
 904
 905   TObjArray arr;
 906   histPart1->FitSlicesY(0,0,-1,10,"QNR",&arr);
 907   TH1D * PartMean1 = (TH1D *) arr.At(1)->Clone("PartMean1");
 908   TH1D * PartSigma1 = (TH1D *) arr.At(2)->Clone("PartSigma1");
 909   PartMean1->SetTitle(";p /(GeV/c) ; Mean (Gauss)");
 910   PartSigma1->SetTitle(";p /(GeV/c) ; Sigma (Gauss)");
 911
 912   histPart1->SetMarkerStyle(22);
 913   PartMean1->SetMarkerStyle(21);
 914   hist->GetAxis(4)->SetRange(0,-1); // RESET RANGES
 915   hist->GetAxis(5)->SetRange(0,-1); // RESET RANGES
 916
 917   if(kParticle2==0){ // electron
 918         nSigmaMinus2   =  -1.999;
 919     nSigmaPlus2   =    2.999;
 920   } else if(kParticle2 == 1){ //pion
 921         nSigmaMinus2   =  -1.999;
 922     nSigmaPlus2   =    2.999;
 923   } else if(kParticle2 == 2){ //kaons
 924         nSigmaMinus2   =  -2.999;
 925     nSigmaPlus2   =    2.999;
 926   } else if(kParticle2 == 3){ //protons
 927         nSigmaMinus2   =  -2.999;
 928     nSigmaPlus2   =    2.999;
 929   }
 930
 931   //
 932   // 2. select and fit second particle
 933   //
 934   hist->GetAxis(3)->SetRangeUser(kParticle2,kParticle2);
 935   hist->GetAxis(4)->SetRangeUser(nSigmaMinus2,nSigmaPlus2);
 936   hist->GetAxis(5)->SetRangeUser(-2.999,2.999); // 3sigma in TOF
 937   TH2D * histPart2 = (TH2D*)hist->Projection(1,0)->Clone("histPart2");
 938   histPart2->RebinX(4);
 939   histPart2->FitSlicesY(0,0,-1,10,"QNR",&arr);
 940   TH1D * PartMean2 =  (TH1D *) arr.At(1)->Clone("PartMean2");
 941   TH1D * PartSigma2 =  (TH1D *) arr.At(2)->Clone("PartSigma2");
 942   PartMean2->SetTitle(";p /(GeV/c) ; Mean (Gauss)");
 943   PartSigma2->SetTitle(";p /(GeV/c) ; Sigma (Gauss)");
 944   PartMean2->SetMarkerStyle(20);
 945   hist->GetAxis(4)->SetRange(0,-1); // RESET RANGES
 946   hist->GetAxis(5)->SetRange(0,-1); // RESET RANGES
 947
 948   //
 949   // 3. separation
 950   //
 951
 952         TH1F *fHistSeparation=(TH1F*)PartMean1->Clone("fHistSeparation"); //to get same binning
 953         fHistSeparation->SetMarkerStyle(22);
 954         const Int_t Nbins = PartMean1->GetNbinsX();
 955         fHistSeparation->SetTitle(";p /(GeV/c) ; Separation");
 956
 957         Float_t DeltaMean[Nbins] ;
 958         Float_t DeltaSigma[Nbins];
 959
 960         for(Int_t i=0 ; i<Nbins; i++){
 961
 962         DeltaMean[i] = TMath::Abs(PartMean1->GetBinContent(i) - PartMean2->GetBinContent(i));
 963         DeltaSigma[i] = TMath::Abs((PartSigma1->GetBinContent(i) + PartSigma2->GetBinContent(i)))/2.;
 964
 965         if(!(TMath::Abs(DeltaSigma[i])<0.000001))fHistSeparation->SetBinContent(i,DeltaMean[i]/DeltaSigma[i]);
 966         }//for(Int_t i=0 ; i<Nbins ; i++)
 967
 968         TObjArray *array = new TObjArray();
 969         array->Add(histPart1);
 970         array->Add(histPart2);
 971         array->Add(PartMean1);
 972         array->Add(PartSigma1);
 973         array->Add(PartMean2);
 974         array->Add(PartSigma2);
 975         array->Add(fHistSeparation);
 976         return array;
 977
 978 }
 979
 980
 981 //________________________________________________________________________
 982 TObjArray * AliTPCcalibResidualPID::GetResidualGraphs(THnSparseF * histPidQA, const Char_t * system, Bool_t useV0s) {
 983   //
 984   // Extracts the residual graphs from THnSparse created from data (NOT MC!)
 985   //
 986
 987   const TString momTitle = "#it{p}_{TPC} (GeV/#it{c})";
 988   const TString dEdxTitle = "d#it{E}/d#it{x} (arb. unit)";
 989
 990   Bool_t isPPb = kFALSE;
 991   if (strcmp(system, "PPB") == 0 || strcmp(system, "PBP") == 0) {
 992     Printf("p-Pb/Pb-p detected - Applying special handling!");
 993     isPPb = kTRUE;
 994   }
 995
 996   // the following parameters have to be extracted from the data itself
 997   //
 998
 999   Float_t nSigmaPionMinus   =  -3.999;
1000   Float_t nSigmaPionPlus    =   3.999;
1001   Float_t nSigmaKaonMinus   =  -4.199;
1002   Float_t nSigmaKaonPlus    =  4.7999;
1003   Float_t nSigmaElectronMinus = -1.999;
1004   Float_t nSigmaElectronPlus  =  4.999;
1005
1006   Int_t cutForFitting = 10;
1007   Double_t heightFractionForFittingRange = 0.1;
1008   //
1009   THnSparse * hist = histPidQA;
1010   //
1011   TCanvas * canvasQAtpc = new TCanvas("canvasQAtpcResGraph","Control canvas for residual graphs (TPC)",100,10,1380,800);
1012   canvasQAtpc->Divide(2,2);
1013
1014   TCanvas * canvasQAtof = new TCanvas("canvasQAtofResGraph","Control canvas for residual graphs (TOF)",100,10,1380,800);
1015   canvasQAtof->Divide(2,2);
1016
1017   TCanvas * canvasQAv0 = new TCanvas("canvasQAv0ResGraph","Control canvas for residual graphs (V0)",100,10,1380,800);
1018   canvasQAv0->Divide(2,2);
1019
1020   TCanvas * canvasQAv0plusTOF = new TCanvas("canvasQAv0plusTOFResGraph","Control canvas for residual graphs (V0+TOF)",100,10,1380,800);
1021   canvasQAv0plusTOF->Divide(2,2);
1022
1023
1024   TCanvas * canvasQAv0DeDxPurityEl = new TCanvas("canvasQAv0DeDxPurityEl","Control canvas for residual graphs (V0 el)",100,10,600,400);
1025   TCanvas * canvasQAv0DeDxPurityPi = new TCanvas("canvasQAv0DeDxPurityPi","Control canvas for residual graphs (V0 pi)",100,10,600,400);
1026   TCanvas * canvasQAv0DeDxPurityPr = new TCanvas("canvasQAv0DeDxPurityPr","Control canvas for residual graphs (V0 pr)",100,10,600,400);
1027
1028   canvasQAv0DeDxPurityEl->SetLogx();
1029   canvasQAv0DeDxPurityPi->SetLogx();
1030   canvasQAv0DeDxPurityPr->SetLogx();
1031
1032   canvasQAv0DeDxPurityEl->SetLogz();
1033   canvasQAv0DeDxPurityPi->SetLogz();
1034   canvasQAv0DeDxPurityPr->SetLogz();
1035
1036   canvasQAv0DeDxPurityEl->SetTopMargin(0.03);
1037   canvasQAv0DeDxPurityPi->SetTopMargin(0.03);
1038   canvasQAv0DeDxPurityPr->SetTopMargin(0.03);
1039
1040   for (Int_t i = 1; i <= 4; i++) {
1041     canvasQAtpc->GetPad(i)->SetGrid(1, 1);
1042     canvasQAtof->GetPad(i)->SetGrid(1, 1);
1043     canvasQAv0->GetPad(i)->SetGrid(1, 1);
1044     canvasQAv0plusTOF->GetPad(i)->SetGrid(1, 1);
1045
1046     canvasQAtpc->GetPad(i)->SetLogz();
1047     canvasQAtof->GetPad(i)->SetLogz();
1048     canvasQAv0->GetPad(i)->SetLogz();
1049     canvasQAv0plusTOF->GetPad(i)->SetLogz();
1050
1051     canvasQAtpc->GetPad(i)->SetLogx();
1052     canvasQAtof->GetPad(i)->SetLogx();
1053     canvasQAv0->GetPad(i)->SetLogx();
1054     canvasQAv0plusTOF->GetPad(i)->SetLogx();
1055   }
1056   //
1057   // 1. select and fit electrons
1058   //
1059   hist->GetAxis(2)->SetRangeUser(0,0); // Non-V0s
1060   hist->GetAxis(3)->SetRangeUser(0,0); // electrons
1061   hist->GetAxis(4)->SetRangeUser(nSigmaElectronMinus,nSigmaElectronPlus); // 3sigma in TPC
1062   hist->GetAxis(5)->SetRangeUser(-2.999,2.999); // 3sigma in TOF
1063   TH2D * histElectron = hist->Projection(1,0);
1064   histElectron->SetName("histElectron");
1065   histElectron->RebinX(4);
1066   histElectron->GetXaxis()->SetRangeUser(0.2,3.0);
1067   TObjArray arr;
1068   FitSlicesY(histElectron, heightFractionForFittingRange, cutForFitting, "QNR", &arr);
1069   TH1D * electronPoints = (TH1D *) arr.At(1)->Clone();
1070   //
1071   TGraphErrors * electronGraph = new TGraphErrors(electronPoints);
1072   electronGraph->SetName("electronGraph");
1073   for(Int_t ip = 0; ip < electronGraph->GetN(); ip ++) {
1074     Bool_t removePoint = electronGraph->GetY()[ip] < 10 ||
1075       electronGraph->GetEY()[ip]/electronGraph->GetY()[ip] > 0.05 ||
1076       electronGraph->GetX()[ip] > 2.0 ||
1077       electronGraph->GetX()[ip] < 0.5;
1078     if (removePoint) {
1079       electronGraph->RemovePoint(ip);
1080       ip--;
1081       continue;
1082     }
1083   }
1084   //
1085   canvasQAtof->cd(1);
1086   histElectron->SetTitle("electrons");
1087   histElectron->GetYaxis()->SetRangeUser(50, 120);
1088   histElectron->GetXaxis()->SetTitle(momTitle.Data());
1089   histElectron->GetYaxis()->SetTitle(dEdxTitle.Data());
1090   histElectron->GetXaxis()->SetMoreLogLabels(kTRUE);
1091   histElectron->GetXaxis()->SetNoExponent(kTRUE);
1092   histElectron->Draw("colz");
1093   electronPoints->SetMarkerStyle(24);
1094   electronPoints->Draw("same");
1095   hist->GetAxis(4)->SetRange(0,-1); // RESET RANGES
1096   hist->GetAxis(5)->SetRange(0,-1); // RESET RANGES
1097   //
1098   // 2. protons
1099   //
1100   hist->GetAxis(2)->SetRangeUser(0,0); // Non-V0s
1101   hist->GetAxis(3)->SetRangeUser(3,3); // protons
1102   //
1103   //hist->GetAxis(4)->SetRangeUser(nSigmaProtonMinus,nSigmaProtonPlus); // protons --> not reliable, use PATTERN RECOGNITION
1104   TH2D * histProtonTPC = hist->Projection(1,0);
1105   histProtonTPC->SetName("histProtonTPC");
1106   histProtonTPC->GetXaxis()->SetRangeUser(0.15,0.7);
1107   histProtonTPC->GetYaxis()->SetRangeUser(50, hist->GetAxis(1)->GetBinUpEdge(hist->GetAxis(1)->GetNbins()));
1108
1109   // PATTERN RECOGNITION
1110   //OLD TF1 betaSq("betaSq","50./TMath::Power(x,1.3)",0.1,10);
1111   // Add some additional Erf - cuts away kaons for pPb and PbPb and does not hurt in pp
1112   TF1 betaSq("betaSq","50./TMath::Power(x,1.3) + (1-TMath::Erf((x-0.2) / 0.075)) * 100",0.1,10);
1113   for(Int_t ix = 1; ix <= histProtonTPC->GetXaxis()->GetNbins(); ix++) {
1114      for(Int_t jy = 1; jy <= histProtonTPC->GetYaxis()->GetNbins(); jy++) {
1115        Float_t yPos = histProtonTPC->GetYaxis()->GetBinCenter(jy);
1116        Float_t xPos = histProtonTPC->GetXaxis()->GetBinCenter(ix);
1117        Int_t bin = histProtonTPC->GetBin(ix,jy);
1118        if (yPos < betaSq.Eval(xPos)) histProtonTPC->SetBinContent(bin,0);
1119      }
1120   }
1121   //
1122   //TODO Use likelihood option -> Inside fitting range ~no contamination, but little statistics at low momenta requires L
1123   FitSlicesY(histProtonTPC, heightFractionForFittingRange, cutForFitting, "QNRL", &arr);
1124   TH1D * protonPointsTPC = (TH1D *) arr.At(1)->Clone();
1125   //
1126   TGraphErrors * protonTpcGraph = new TGraphErrors(protonPointsTPC);
1127   protonTpcGraph->SetName("protonTpcGraph");
1128   for(Int_t ip = 0; ip < protonTpcGraph->GetN(); ip ++) {
1129     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
1130                         || 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
1131                         //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
1132     if (removePoint) {
1133       protonTpcGraph->RemovePoint(ip);
1134       ip--;
1135       continue;
1136     }
1137   }
1138   //
1139   hist->GetAxis(5)->SetRangeUser(-2.999,2.999); // 3sigma in TOF
1140   TH2D * histProtonTOF = hist->Projection(1,0);
1141   histProtonTOF->SetName("histProtonTOF");
1142   histProtonTOF->GetXaxis()->SetRangeUser(0.6,2.5);
1143
1144   // Same pattern recognition as before
1145   for(Int_t ix = 1; ix <= histProtonTOF->GetXaxis()->GetNbins(); ix++) {
1146     for(Int_t jy = 1; jy <= histProtonTOF->GetYaxis()->GetNbins(); jy++) {
1147       Float_t yPos = histProtonTOF->GetYaxis()->GetBinCenter(jy);
1148       Float_t xPos = histProtonTOF->GetXaxis()->GetBinCenter(ix);
1149       Int_t bin = histProtonTOF->GetBin(ix,jy);
1150       if (yPos < betaSq.Eval(xPos)) histProtonTOF->SetBinContent(bin,0);
1151     }
1152   }
1153
1154   FitSlicesY(histProtonTOF, heightFractionForFittingRange, cutForFitting, "QNR", &arr);
1155
1156   TH1D * protonPointsTOF =  (TH1D *)arr.At(1)->Clone();
1157   //
1158   TGraphErrors * protonTofGraph = new TGraphErrors(protonPointsTOF);
1159   protonTofGraph->SetName("protonTofGraph");
1160   for(Int_t ip = 0; ip < protonTofGraph->GetN(); ip ++) {
1161     // If V0s are to be used, TPC+TOF protons are only for reference and can be plotted to higher momenta
1162     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;
1163     if (removePoint) {
1164       protonTofGraph->RemovePoint(ip);
1165       ip--;
1166       continue;
1167     }
1168   }
1169   //
1170   canvasQAtpc->cd(2);
1171   histProtonTPC->SetTitle("protons");
1172   histProtonTPC->GetXaxis()->SetTitle(momTitle.Data());
1173   histProtonTPC->GetYaxis()->SetTitle(dEdxTitle.Data());
1174   histProtonTPC->GetXaxis()->SetMoreLogLabels(kTRUE);
1175   histProtonTPC->GetXaxis()->SetNoExponent(kTRUE);
1176   histProtonTPC->Draw("colz");
1177   betaSq.DrawCopy("same");
1178   protonPointsTPC->SetMarkerStyle(20);
1179   protonPointsTOF->SetMarkerStyle(24);
1180   protonPointsTOF->Draw("same");
1181   protonPointsTPC->Draw("same");
1182   //
1183   protonTpcGraph->SetMarkerStyle(26);
1184   protonTpcGraph->SetMarkerColor(kMagenta);
1185   protonTpcGraph->DrawClone("p");
1186   protonTofGraph->SetMarkerStyle(25);
1187   protonTofGraph->SetMarkerColor(kMagenta);
1188   protonTofGraph->DrawClone("p");
1189
1190   canvasQAtof->cd(2);
1191   histProtonTOF->SetTitle("protons");
1192   histProtonTOF->GetYaxis()->SetRangeUser(30, 250);
1193   histProtonTOF->GetXaxis()->SetTitle(momTitle.Data());
1194   histProtonTOF->GetYaxis()->SetTitle(dEdxTitle.Data());
1195   histProtonTOF->GetXaxis()->SetMoreLogLabels(kTRUE);
1196   histProtonTOF->GetXaxis()->SetNoExponent(kTRUE);
1197   histProtonTOF->Draw("colz");
1198   betaSq.DrawCopy("same");
1199   protonPointsTOF->Draw("same");
1200   protonPointsTPC->Draw("same");
1201   //
1202   protonTpcGraph->DrawClone("p");
1203   protonTofGraph->DrawClone("p");
1204
1205   hist->GetAxis(4)->SetRange(0,-1); // RESET RANGES
1206   hist->GetAxis(5)->SetRange(0,-1); // RESET RANGES
1207   //
1208   // 3. pions
1209   //
1210   hist->GetAxis(2)->SetRangeUser(0,0); // Non-V0s
1211   hist->GetAxis(3)->SetRangeUser(1,1); // pions
1212   //
1213   Double_t lowerPionThreshold = 0.3;
1214
1215   hist->GetAxis(4)->SetRangeUser(nSigmaPionMinus,nSigmaPionPlus); // pions
1216   TH2D * histPionTPC = hist->Projection(1,0);
1217   histPionTPC->SetName("histPionTPC");
1218   histPionTPC->GetXaxis()->SetRangeUser(0.15,0.7);
1219   FitSlicesY(histPionTPC, heightFractionForFittingRange, cutForFitting, "QNR", &arr);
1220
1221   // In case of really bad splines comment last 5 lines and use the following 6 lines:
1222   //TH2D * histPionTPC = hist->Projection(1,0);
1223   //histPionTPC->SetName("histPionTPC");
1224   //histPionTPC->GetYaxis()->SetRangeUser(30,100);
1225   //histPionTPC->GetXaxis()->SetRangeUser(0.15,0.7);
1226   //FitSlicesY(histPionTPC, heightFractionForFittingRange, cutForFitting, "QNR", &arr);
1227   //lowerPionThreshold = 0.15;
1228
1229
1230   TH1D * pionPointsTPC =  (TH1D *) arr.At(1)->Clone();
1231   //
1232   TGraphErrors * pionTpcGraph = new TGraphErrors(pionPointsTPC);
1233   pionTpcGraph->SetName("pionTpcGraph");
1234   for(Int_t ip = 0; ip < pionTpcGraph->GetN(); ip ++) {
1235     Bool_t removePoint = pionTpcGraph->GetY()[ip] < 10 || pionTpcGraph->GetEY()[ip]/pionTpcGraph->GetY()[ip] > 0.02 || pionTpcGraph->GetX()[ip] > 0.5
1236                          || pionTpcGraph->GetX()[ip] < lowerPionThreshold; // Exclude contamination by electrons
1237     if (removePoint) {
1238       pionTpcGraph->RemovePoint(ip);
1239       ip--;
1240       continue;
1241     }
1242   }
1243   //
1244   hist->GetAxis(5)->SetRangeUser(-2.999,2.999); // 3sigma in TOF
1245   TH2D * histPionTOF = hist->Projection(1,0);
1246   histPionTOF->SetName("histPionTOF");
1247   histPionTOF->GetXaxis()->SetRangeUser(0.5,1.1);
1248   FitSlicesY(histPionTOF, heightFractionForFittingRange, cutForFitting, "QNR", &arr);
1249   TH1D * pionPointsTOF = (TH1D *) arr.At(1)->Clone();
1250   TGraphErrors * pionTofGraph = new TGraphErrors(pionPointsTOF);
1251   pionTofGraph->SetName("pionTofGraph");
1252   for(Int_t ip = 0; ip < pionTofGraph->GetN(); ip ++) {
1253     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
1254     if (removePoint) {
1255       pionTofGraph->RemovePoint(ip);
1256       ip--;
1257       continue;
1258     }
1259   }
1260   canvasQAtpc->cd(3);
1261   histPionTPC->GetYaxis()->SetRangeUser(30, 90);
1262   histPionTPC->SetTitle("pions");
1263   histPionTPC->GetXaxis()->SetTitle(momTitle.Data());
1264   histPionTPC->GetYaxis()->SetTitle(dEdxTitle.Data());
1265   histPionTPC->GetXaxis()->SetMoreLogLabels(kTRUE);
1266   histPionTPC->GetXaxis()->SetNoExponent(kTRUE);
1267   histPionTPC->Draw("colz");
1268   pionPointsTPC->SetMarkerStyle(20);
1269   pionPointsTOF->SetMarkerStyle(24);
1270   pionPointsTOF->Draw("same");
1271   pionPointsTPC->Draw("same");
1272   //
1273   pionTpcGraph->SetMarkerStyle(26);
1274   pionTpcGraph->SetMarkerColor(kMagenta);
1275   pionTpcGraph->DrawClone("p");
1276   pionTofGraph->SetMarkerStyle(25);
1277   pionTofGraph->SetMarkerColor(kMagenta);
1278   pionTofGraph->DrawClone("p");
1279
1280   canvasQAtof->cd(3);
1281   histPionTOF->GetYaxis()->SetRangeUser(30, 80);
1282   histPionTOF->GetXaxis()->SetTitle(momTitle.Data());
1283   histPionTOF->GetYaxis()->SetTitle(dEdxTitle.Data());
1284   histPionTOF->GetXaxis()->SetMoreLogLabels(kTRUE);
1285   histPionTOF->GetXaxis()->SetNoExponent(kTRUE);
1286   histPionTOF->Draw("colz");
1287   histPionTOF->SetTitle("pions");
1288   pionPointsTOF->Draw("same");
1289   pionPointsTPC->Draw("same");
1290   //
1291   pionTpcGraph->DrawClone("p");
1292   pionTofGraph->DrawClone("p");
1293
1294
1295   hist->GetAxis(4)->SetRange(0,-1); // RESET RANGES
1296   hist->GetAxis(5)->SetRange(0,-1); // RESET RANGES
1297   //
1298   // 4. kaons
1299   //
1300   hist->GetAxis(2)->SetRangeUser(0,0); // Non-V0s
1301   hist->GetAxis(3)->SetRangeUser(2,2); // kaons
1302   //
1303   hist->GetAxis(4)->SetRangeUser(nSigmaKaonMinus,nSigmaKaonPlus); // kaons
1304   TH2D * histKaonTPC = hist->Projection(1,0);
1305   histKaonTPC->SetName("histKaonTPC");
1306   histKaonTPC->GetXaxis()->SetRangeUser(0.12,0.6999);
1307   FitSlicesY(histKaonTPC, heightFractionForFittingRange, cutForFitting, "QNR", &arr);
1308   TH1D * kaonPointsTPC = (TH1D*) arr.At(1)->Clone();
1309   //
1310   TGraphErrors * kaonTpcGraph = new TGraphErrors(kaonPointsTPC);
1311   kaonTpcGraph->SetName("kaonTpcGraph");
1312   for(Int_t ip = 0; ip < kaonTpcGraph->GetN(); ip ++) {
1313     Bool_t removePoint = kaonTpcGraph->GetY()[ip] < 10 || kaonTpcGraph->GetEY()[ip]/kaonTpcGraph->GetY()[ip] > 0.02
1314                          || kaonTpcGraph->GetX()[ip] < 0.12 || kaonTpcGraph->GetX()[ip] > 0.319; //TODO BEN was > 0.399
1315     if (removePoint) {
1316       kaonTpcGraph->RemovePoint(ip);
1317       ip--;
1318       continue;
1319     }
1320   }
1321   //
1322
1323   hist->GetAxis(5)->SetRangeUser(-2.999,2.999); // sigma in TOF
1324
1325   /*
1326   hist->GetAxis(5)->SetRangeUser(-1.999,1.999); // sigma in TOF
1327
1328   if (hist->GetAxis(3)->GetNbins() >= 5) {
1329     printf("Using restricted eta range for TPC+TOF Kaons!\n");
1330     hist->GetAxis(3)->SetRangeUser(4,4);
1331   }
1332   else {
1333     printf("WARNING: Data points for restricted eta range for TPC+TOF Kaons not available! Using full eta range (worse w.r.t. contamination)\n");
1334   }*/
1335
1336   TH2D * histKaonTOF = hist->Projection(1,0);
1337   histKaonTOF->SetName("histKaonTOF");
1338   histKaonTOF->GetXaxis()->SetRangeUser(0.3,1.5);
1339   FitSlicesY(histKaonTOF, heightFractionForFittingRange, cutForFitting, "QNR", &arr);
1340   TH1D * kaonPointsTOF = (TH1D*) arr.At(1)->Clone();
1341   //
1342   TGraphErrors * kaonTofGraph = new TGraphErrors(kaonPointsTOF);
1343   kaonTofGraph->SetName("kaonTofGraph");
1344   //
1345   for(Int_t ip = 0; ip < kaonTofGraph->GetN(); ip ++) {
1346     Bool_t removePoint = kaonTofGraph->GetY()[ip] < 10 || kaonTofGraph->GetEY()[ip]/kaonTofGraph->GetY()[ip] > 0.02 || kaonTofGraph->GetX()[ip] > 1.0
1347                          || kaonTofGraph->GetX()[ip] < 0.36; //TODO BEN NOW was 0.5
1348     if (removePoint) {
1349       kaonTofGraph->RemovePoint(ip);
1350       ip--;
1351       continue;
1352     }
1353   }
1354   //
1355   canvasQAtpc->cd(4);
1356   histKaonTPC->SetTitle("kaons");
1357   histKaonTPC->GetYaxis()->SetRangeUser(30, 500);
1358   histKaonTPC->GetXaxis()->SetTitle(momTitle.Data());
1359   histKaonTPC->GetYaxis()->SetTitle(dEdxTitle.Data());
1360   histKaonTPC->GetXaxis()->SetMoreLogLabels(kTRUE);
1361   histKaonTPC->GetXaxis()->SetNoExponent(kTRUE);
1362   histKaonTPC->Draw("colz");
1363   kaonPointsTPC->SetMarkerStyle(20);
1364   kaonPointsTOF->SetMarkerStyle(24);
1365   kaonPointsTOF->Draw("same");
1366   kaonPointsTPC->Draw("same");
1367   //
1368   kaonTpcGraph->SetMarkerStyle(26);
1369   kaonTpcGraph->SetMarkerColor(kMagenta);
1370   kaonTpcGraph->DrawClone("p");
1371   kaonTofGraph->SetMarkerStyle(25);
1372   kaonTofGraph->SetMarkerColor(kMagenta);
1373   kaonTofGraph->DrawClone("p");
1374
1375
1376   canvasQAtof->cd(4);
1377   histKaonTOF->GetYaxis()->SetRangeUser(30, 250);
1378   histKaonTOF->SetTitle("kaons");
1379   histKaonTOF->GetXaxis()->SetTitle(momTitle.Data());
1380   histKaonTOF->GetYaxis()->SetTitle(dEdxTitle.Data());
1381   histKaonTOF->GetXaxis()->SetMoreLogLabels(kTRUE);
1382   histKaonTOF->GetXaxis()->SetNoExponent(kTRUE);
1383   histKaonTOF->Draw("colz");
1384   kaonPointsTOF->Draw("same");
1385   kaonPointsTPC->Draw("same");
1386   kaonTpcGraph->DrawClone("p");
1387   kaonTofGraph->DrawClone("p");
1388
1389   hist->GetAxis(4)->SetRange(0,-1); // RESET RANGES
1390   hist->GetAxis(5)->SetRange(0,-1); // RESET RANGES
1391
1392
1393   //
1394   // 5. V0 electrons
1395   //
1396   hist->GetAxis(2)->SetRangeUser(1,1); // V0 electrons
1397   hist->GetAxis(3)->SetRangeUser(0,0); // electrons
1398
1399   //
1400   TH2D * histElectronV0 = hist->Projection(1,0);
1401   histElectronV0->SetName("histElectronV0");
1402   //histElectronV0->RebinX(2);
1403   histElectronV0->GetXaxis()->SetRangeUser(0.15,5.0);
1404   FitSlicesY(histElectronV0, heightFractionForFittingRange, cutForFitting, "QNR", &arr);
1405   TH1D * electronPointsV0 = (TH1D*) arr.At(1)->Clone();
1406   //
1407   TGraphErrors * electronV0Graph = new TGraphErrors(electronPointsV0);
1408   electronV0Graph->SetName("electronV0Graph");
1409   for(Int_t ip = 0; ip < electronV0Graph->GetN(); ip ++) {
1410     Bool_t removePoint = electronV0Graph->GetY()[ip] < 10 || electronV0Graph->GetEY()[ip]/electronV0Graph->GetY()[ip] > 0.02;
1411
1412     if (removePoint) {
1413       electronV0Graph->RemovePoint(ip);
1414       ip--;
1415       continue;
1416     }
1417   }
1418   //
1419   canvasQAv0->cd(1);
1420   histElectronV0->SetTitle("V0 electrons");
1421   histElectronV0->GetYaxis()->SetRangeUser(50, 120);
1422   histElectronV0->GetXaxis()->SetTitle(momTitle.Data());
1423   histElectronV0->GetYaxis()->SetTitle(dEdxTitle.Data());
1424   histElectronV0->GetXaxis()->SetMoreLogLabels(kTRUE);
1425   histElectronV0->GetXaxis()->SetNoExponent(kTRUE);
1426   histElectronV0->SetStats(0);
1427   histElectronV0->Draw("colz");
1428   electronPointsV0->SetMarkerStyle(34);
1429   electronPointsV0->Draw("same");
1430   electronGraph->SetMarkerStyle(25);
1431   electronGraph->SetMarkerColor(kMagenta);
1432   electronGraph->DrawClone("p");
1433
1434   canvasQAv0DeDxPurityEl->cd();
1435   gStyle->SetOptTitle(0);
1436   TH1D* hElNew = (TH1D*)histElectronV0->DrawClone("colz");
1437   hElNew->GetXaxis()->SetTitleOffset(1.0);
1438   hElNew->SetTitle("");
1439   electronPointsV0->DrawClone("same");
1440   gStyle->SetOptTitle(1);
1441
1442
1443
1444   canvasQAtof->cd(1);
1445   electronV0Graph->SetMarkerStyle(24);
1446   electronV0Graph->SetMarkerColor(kMagenta);
1447   electronV0Graph->DrawClone("p");
1448
1449   //
1450   // 6. V0 pions
1451   //
1452   hist->GetAxis(2)->SetRangeUser(2,2); // V0 pions
1453   hist->GetAxis(3)->SetRangeUser(1,1); // pions
1454   //
1455   TH2D * histPionV0 = hist->Projection(1,0);
1456   histPionV0->SetName("histPionV0");
1457   histPionV0->GetXaxis()->SetRangeUser(0.15,10.);
1458   FitSlicesY(histPionV0, heightFractionForFittingRange, cutForFitting, "QNR", &arr);
1459   TH1D * pionPointsV0 = (TH1D*) arr.At(1)->Clone();
1460   //
1461   TGraphErrors * pionV0Graph = new TGraphErrors(pionPointsV0);
1462   pionV0Graph->SetName("pionV0Graph");
1463   for(Int_t ip = 0; ip < pionV0Graph->GetN(); ip ++) {
1464     Bool_t removePoint = pionV0Graph->GetY()[ip] < 10 || pionV0Graph->GetEY()[ip]/pionV0Graph->GetY()[ip] > 0.02;
1465
1466     if (removePoint) {
1467       pionV0Graph->RemovePoint(ip);
1468       ip--;
1469       continue;
1470     }
1471   }
1472   //
1473   canvasQAv0->cd(3);
1474   histPionV0->SetTitle("V0 pions");
1475   histPionV0->GetYaxis()->SetRangeUser(30, 100);
1476   histPionV0->GetXaxis()->SetTitle(momTitle.Data());
1477   histPionV0->GetYaxis()->SetTitle(dEdxTitle.Data());
1478   histPionV0->GetXaxis()->SetMoreLogLabels(kTRUE);
1479   histPionV0->GetXaxis()->SetNoExponent(kTRUE);
1480   histPionV0->Draw("colz");
1481   pionPointsV0->SetMarkerStyle(34);
1482   pionPointsV0->Draw("same");
1483   pionTofGraph->DrawClone("p");
1484   pionTpcGraph->DrawClone("p");
1485
1486   canvasQAv0DeDxPurityPi->cd();
1487   gStyle->SetOptTitle(0);
1488   TH1D* hPiNew = (TH1D*)histPionV0->DrawClone("colz");
1489   hPiNew->GetXaxis()->SetTitleOffset(1.0);
1490   hPiNew->SetTitle("");
1491   pionPointsV0->DrawClone("same");
1492   gStyle->SetOptTitle(1);
1493
1494   canvasQAtof->cd(3);
1495   pionV0Graph->SetMarkerStyle(24);
1496   pionV0Graph->SetMarkerColor(kMagenta);
1497   pionV0Graph->DrawClone("p");
1498
1499   canvasQAtpc->cd(3);
1500   pionV0Graph->DrawClone("p");
1501
1502   //
1503   // 6. V0 protons
1504   //
1505   hist->GetAxis(2)->SetRangeUser(3,3); // V0 protons
1506   hist->GetAxis(3)->SetRangeUser(3,3); // protons
1507   //
1508   TH2D * histProtonV0 = hist->Projection(1,0);
1509   histProtonV0->SetName("histProtonV0");
1510   histProtonV0->GetXaxis()->SetRangeUser(0.15,10.);
1511
1512   // Remove contamination due to el and pions and low momenta because there the statistics
1513   // for the V0 protons goes down and becomes comparable with the contamination
1514   // -> This spoils the fit, if the contamination is not removed
1515
1516   // PATTERN RECOGNITION
1517   Int_t correctionUpToBin = histProtonV0->GetXaxis()->FindBin(0.6);
1518
1519   for(Int_t ix = 1; ix <= correctionUpToBin; ix++) {
1520      for(Int_t jy = 1; jy <= histProtonV0->GetYaxis()->GetNbins(); jy++) {
1521        Float_t yPos = histProtonV0->GetYaxis()->GetBinCenter(jy);
1522        Float_t xPos = histProtonV0->GetXaxis()->GetBinCenter(ix);
1523        Int_t bin = histProtonV0->GetBin(ix,jy);
1524        if (yPos < betaSq.Eval(xPos)) histProtonV0->SetBinContent(bin,0);
1525      }
1526   }
1527
1528   /*
1529   Int_t correctionUpToBin = histProtonV0->GetXaxis()->FindBin(0.4);
1530   Double_t threshold = 120.;
1531
1532   for(Int_t ix = 1; ix <= correctionUpToBin; ix++) {
1533      for(Int_t jy = 1; jy <= histProtonV0->GetYaxis()->GetNbins(); jy++) {
1534        Float_t yPos = histProtonV0->GetYaxis()->GetBinCenter(jy);
1535        Int_t bin = histProtonV0->GetBin(ix,jy);
1536        if (yPos < threshold) histProtonV0->SetBinContent(bin,0);
1537        else break;
1538      }
1539   }
1540   */
1541
1542   FitSlicesY(histProtonV0, heightFractionForFittingRange, cutForFitting, "QNR", &arr);
1543   TH1D * protonPointsV0 = (TH1D*) arr.At(1)->Clone();
1544   //
1545   TGraphErrors * protonV0Graph = new TGraphErrors(protonPointsV0);
1546   protonV0Graph->SetName("protonV0Graph");
1547   for(Int_t ip = 0; ip < protonV0Graph->GetN(); ip ++) {
1548     Bool_t removePoint = protonV0Graph->GetY()[ip] < 10 || protonV0Graph->GetEY()[ip]/protonV0Graph->GetY()[ip] > 0.02;
1549
1550     if (removePoint) {
1551       protonV0Graph->RemovePoint(ip);
1552       ip--;
1553       continue;
1554     }
1555   }
1556   //
1557   canvasQAv0->cd(2);
1558   histProtonV0->SetTitle("V0 protons");
1559   histProtonV0->GetXaxis()->SetTitle(momTitle.Data());
1560   histProtonV0->GetYaxis()->SetTitle(dEdxTitle.Data());
1561   histProtonV0->GetXaxis()->SetMoreLogLabels(kTRUE);
1562   histProtonV0->GetXaxis()->SetNoExponent(kTRUE);
1563   histProtonV0->Draw("colz");
1564   protonPointsV0->SetMarkerStyle(34);
1565   protonPointsV0->Draw("same");
1566   protonTofGraph->DrawClone("p");
1567   protonTpcGraph->DrawClone("p");
1568
1569   canvasQAv0DeDxPurityPr->cd();
1570   gStyle->SetOptTitle(0);
1571   TH1D* hPrNew = (TH1D*)histProtonV0->DrawClone("colz");
1572   hPrNew->GetXaxis()->SetTitleOffset(1.0);
1573   hPrNew->SetTitle("");
1574   protonPointsV0->DrawClone("same");
1575   gStyle->SetOptTitle(1);
1576
1577   canvasQAtof->cd(2);
1578   protonV0Graph->SetMarkerStyle(24);
1579   protonV0Graph->SetMarkerColor(kMagenta);
1580   protonV0Graph->DrawClone("p");
1581
1582   canvasQAtpc->cd(2);
1583   protonV0Graph->DrawClone("p");
1584
1585
1586   hist->GetAxis(2)->SetRange(0,-1); // RESET RANGES
1587   hist->GetAxis(3)->SetRange(0,-1); // RESET RANGES
1588
1589
1590   //
1591   // 5. V0 electrons + TOF
1592   //
1593   hist->GetAxis(2)->SetRangeUser(1,1); // V0 electrons
1594   hist->GetAxis(3)->SetRangeUser(0,0); // electrons
1595   hist->GetAxis(5)->SetRangeUser(-2.999,2.999); // 3sigma in TOF
1596   //
1597   TH2D * histElectronV0plusTOF = hist->Projection(1,0);
1598   histElectronV0plusTOF->SetName("histElectronV0plusTOF");
1599   histElectronV0plusTOF->RebinX(2);
1600   histElectronV0plusTOF->GetXaxis()->SetRangeUser(0.2,5.0);
1601   FitSlicesY(histElectronV0plusTOF, heightFractionForFittingRange, cutForFitting, "QNR", &arr);
1602   TH1D * electronPointsV0plusTOF = (TH1D*) arr.At(1)->Clone();
1603   //
1604   TGraphErrors * electronV0plusTOFGraph = new TGraphErrors(electronPointsV0plusTOF);
1605   electronV0plusTOFGraph->SetName("electronV0plusTOFGraph");
1606   for(Int_t ip = 0; ip < electronV0plusTOFGraph->GetN(); ip ++) {
1607     Bool_t removePoint = electronV0plusTOFGraph->GetY()[ip] < 10 || electronV0plusTOFGraph->GetEY()[ip]/electronV0plusTOFGraph->GetY()[ip] > 0.02;
1608
1609     if (removePoint) {
1610       electronV0plusTOFGraph->RemovePoint(ip);
1611       ip--;
1612       continue;
1613     }
1614   }
1615   //
1616   canvasQAv0plusTOF->cd(1);
1617   histElectronV0plusTOF->SetTitle("V0+TOF electrons");
1618   histElectronV0plusTOF->GetYaxis()->SetRangeUser(50, 120);
1619   histElectronV0plusTOF->GetXaxis()->SetTitle(momTitle.Data());
1620   histElectronV0plusTOF->GetYaxis()->SetTitle(dEdxTitle.Data());
1621   histElectronV0plusTOF->GetXaxis()->SetMoreLogLabels(kTRUE);
1622   histElectronV0plusTOF->GetXaxis()->SetNoExponent(kTRUE);
1623   histElectronV0plusTOF->Draw("colz");
1624   electronPointsV0plusTOF->SetMarkerStyle(29);
1625   electronPointsV0plusTOF->Draw("same");
1626   electronGraph->DrawClone("p");
1627   electronV0Graph->DrawClone("p");
1628
1629   canvasQAv0->cd(1);
1630   electronV0plusTOFGraph->SetMarkerStyle(30);
1631   electronV0plusTOFGraph->SetMarkerColor(kMagenta);
1632   electronV0plusTOFGraph->DrawClone("p");
1633
1634   canvasQAtof->cd(1);
1635   electronV0plusTOFGraph->DrawClone("p");
1636
1637   //
1638   // 6. V0 pions + TOF
1639   //
1640   hist->GetAxis(2)->SetRangeUser(2,2); // V0 pions
1641   hist->GetAxis(3)->SetRangeUser(1,1); // pions
1642   hist->GetAxis(5)->SetRangeUser(-2.999,2.999); // 3sigma in TOF
1643   //
1644   TH2D * histPionV0plusTOF = hist->Projection(1,0);
1645   histPionV0plusTOF->SetName("histPionV0plusTOF");
1646   histPionV0plusTOF->GetXaxis()->SetRangeUser(0.15,10.);
1647   FitSlicesY(histPionV0plusTOF, heightFractionForFittingRange, cutForFitting, "QNR", &arr);
1648   TH1D * pionPointsV0plusTOF = (TH1D*) arr.At(1)->Clone();
1649   //
1650   TGraphErrors * pionV0plusTOFGraph = new TGraphErrors(pionPointsV0plusTOF);
1651   pionV0plusTOFGraph->SetName("pionV0plusTOFGraph");
1652   for(Int_t ip = 0; ip < pionV0plusTOFGraph->GetN(); ip ++) {
1653     Bool_t removePoint = pionV0plusTOFGraph->GetY()[ip] < 10 || pionV0plusTOFGraph->GetEY()[ip]/pionV0plusTOFGraph->GetY()[ip] > 0.02;
1654
1655     if (removePoint) {
1656       pionV0plusTOFGraph->RemovePoint(ip);
1657       ip--;
1658       continue;
1659     }
1660   }
1661   //
1662   canvasQAv0plusTOF->cd(3);
1663   histPionV0plusTOF->SetTitle("V0+TOF pions");
1664   histPionV0plusTOF->GetYaxis()->SetRangeUser(30, 100);
1665   histPionV0plusTOF->GetXaxis()->SetTitle(momTitle.Data());
1666   histPionV0plusTOF->GetYaxis()->SetTitle(dEdxTitle.Data());
1667   histPionV0plusTOF->GetXaxis()->SetMoreLogLabels(kTRUE);
1668   histPionV0plusTOF->GetXaxis()->SetNoExponent(kTRUE);
1669   histPionV0plusTOF->Draw("colz");
1670   pionPointsV0plusTOF->SetMarkerStyle(29);
1671   pionPointsV0plusTOF->Draw("same");
1672   pionV0Graph->DrawClone("p");
1673   pionTofGraph->DrawClone("p");
1674   pionTpcGraph->DrawClone("p");
1675
1676   canvasQAv0->cd(3);
1677   pionV0plusTOFGraph->SetMarkerStyle(30);
1678   pionV0plusTOFGraph->SetMarkerColor(kMagenta);
1679   pionV0plusTOFGraph->DrawClone("p");
1680
1681   canvasQAtof->cd(3);
1682   pionV0plusTOFGraph->DrawClone("p");
1683
1684   canvasQAtpc->cd(3);
1685   pionV0plusTOFGraph->DrawClone("p");
1686
1687   //
1688   // 6. V0 protons + TOF
1689   //
1690   hist->GetAxis(2)->SetRangeUser(3,3); // V0 protons
1691   hist->GetAxis(3)->SetRangeUser(3,3); // protons
1692   hist->GetAxis(5)->SetRangeUser(-2.999,2.999); // 3sigma in TOF
1693   //
1694   TH2D * histProtonV0plusTOF = hist->Projection(1,0);
1695   histProtonV0plusTOF->SetName("histProtonV0plusTOF");
1696   histProtonV0plusTOF->GetXaxis()->SetRangeUser(0.15,10.);
1697
1698   // Remove contamination due to el and pions and low momenta because there the statistics
1699   // for the V0 protons goes down and becomes comparable with the contamination
1700   // -> This spoils the fit, if the contamination is not removed
1701
1702   // PATTERN RECOGNITION
1703   correctionUpToBin = histProtonV0plusTOF->GetXaxis()->FindBin(0.6);
1704
1705   for(Int_t ix = 1; ix <= correctionUpToBin; ix++) {
1706      for(Int_t jy = 1; jy <= histProtonV0plusTOF->GetYaxis()->GetNbins(); jy++) {
1707        Float_t yPos = histProtonV0plusTOF->GetYaxis()->GetBinCenter(jy);
1708        Float_t xPos = histProtonV0plusTOF->GetXaxis()->GetBinCenter(ix);
1709        Int_t bin = histProtonV0plusTOF->GetBin(ix,jy);
1710        if (yPos < betaSq.Eval(xPos)) histProtonV0plusTOF->SetBinContent(bin,0);
1711      }
1712   }
1713
1714   FitSlicesY(histProtonV0plusTOF, heightFractionForFittingRange, cutForFitting, "QNR", &arr);
1715   TH1D * protonPointsV0plusTOF = (TH1D*) arr.At(1)->Clone();
1716   //
1717   TGraphErrors * protonV0plusTOFGraph = new TGraphErrors(protonPointsV0plusTOF);
1718   protonV0plusTOFGraph->SetName("protonV0plusTOFGraph");
1719   for(Int_t ip = 0; ip < protonV0plusTOFGraph->GetN(); ip ++) {
1720     Bool_t removePoint = protonV0plusTOFGraph->GetY()[ip] < 10 || protonV0plusTOFGraph->GetEY()[ip]/protonV0plusTOFGraph->GetY()[ip] > 0.02;
1721
1722     if (removePoint) {
1723       protonV0plusTOFGraph->RemovePoint(ip);
1724       ip--;
1725       continue;
1726     }
1727   }
1728   //
1729   canvasQAv0plusTOF->cd(2);
1730   histProtonV0plusTOF->SetTitle("V0+TOF protons");
1731   histProtonV0plusTOF->GetXaxis()->SetTitle(momTitle.Data());
1732   histProtonV0plusTOF->GetYaxis()->SetTitle(dEdxTitle.Data());
1733   histProtonV0plusTOF->GetXaxis()->SetMoreLogLabels(kTRUE);
1734   histProtonV0plusTOF->GetXaxis()->SetNoExponent(kTRUE);
1735   histProtonV0plusTOF->Draw("colz");
1736   protonPointsV0plusTOF->SetMarkerStyle(29);
1737   protonPointsV0plusTOF->Draw("same");
1738   protonV0Graph->DrawClone("p");
1739   protonTofGraph->DrawClone("p");
1740   protonTpcGraph->DrawClone("p");
1741
1742   canvasQAv0->cd(2);
1743   protonV0plusTOFGraph->SetMarkerStyle(30);
1744   protonV0plusTOFGraph->SetMarkerColor(kMagenta);
1745   protonV0plusTOFGraph->DrawClone("p");
1746
1747   canvasQAtof->cd(2);
1748   protonV0plusTOFGraph->DrawClone("p");
1749
1750   canvasQAtpc->cd(2);
1751   protonV0plusTOFGraph->DrawClone("p");
1752
1753   hist->GetAxis(2)->SetRange(0,-1); // RESET RANGES
1754   hist->GetAxis(3)->SetRange(0,-1); // RESET RANGES
1755   hist->GetAxis(5)->SetRange(0,-1); // RESET RANGES
1756
1757
1758   // Clone (V0, V0+TOF) graphs: Remove some points, which are expected to deviate from the Bethe-Bloch curve, from the original graphs, so that
1759   // these graphs can be used for the BB fit. The original graph will keep all point in order to determine the correction and response
1760   // functions
1761   TGraphErrors * electronV0GraphForBBfit = new TGraphErrors(*electronV0Graph);
1762   electronV0GraphForBBfit->SetName("electronV0GraphForBBfit");
1763
1764   for(Int_t ip = 0; ip < electronV0GraphForBBfit->GetN(); ip ++) {
1765     Bool_t removePoint = electronV0GraphForBBfit->GetX()[ip] < (isPPb ? 1.2 : 0.4);// || electronV0GraphForBBfit->GetX()[ip] > 2.2;
1766
1767     if (removePoint) {
1768       electronV0GraphForBBfit->RemovePoint(ip);
1769       ip--;
1770       continue;
1771     }
1772   }
1773
1774   TGraphErrors * pionV0GraphForBBfit = new TGraphErrors(*pionV0Graph);
1775   pionV0GraphForBBfit->SetName("pionV0GraphForBBfit");
1776
1777   for(Int_t ip = 0; ip < pionV0GraphForBBfit->GetN(); ip ++) {
1778     Bool_t removePoint = pionV0GraphForBBfit->GetX()[ip] < (isPPb ? 1.0 : 0.76);// || pionV0GraphForBBfit->GetX()[ip] > 6.; //TODO NOW was < 0.4 before
1779     //Bool_t removePoint = pionV0GraphForBBfit->GetX()[ip] < 0.4 || pionV0GraphForBBfit->GetX()[ip] > 0.399;// In this case NOT used for BB fit
1780     if (removePoint) {
1781       pionV0GraphForBBfit->RemovePoint(ip);
1782       ip--;
1783       continue;
1784     }
1785   }
1786
1787   TGraphErrors * protonV0GraphForBBfit = new TGraphErrors(*protonV0Graph);
1788   protonV0GraphForBBfit->SetName("protonV0GraphForBBfit");
1789
1790   for(Int_t ip = 0; ip < protonV0GraphForBBfit->GetN(); ip ++) {
1791     Bool_t removePoint = protonV0GraphForBBfit->GetX()[ip] < 0.9 || protonV0GraphForBBfit->GetX()[ip] > 5.0; //TODO NOW was 2.5 before
1792     //Bool_t removePoint = protonV0GraphForBBfit->GetX()[ip] < 0.9 || protonV0GraphForBBfit->GetX()[ip] > 0.599;// In this case NOT used for BB fit
1793     if (removePoint) {
1794       protonV0GraphForBBfit->RemovePoint(ip);
1795       ip--;
1796       continue;
1797     }
1798   }
1799
1800
1801   // V0 plus TOF
1802   TGraphErrors * electronV0plusTOFGraphForBBfit = new TGraphErrors(*electronV0plusTOFGraph);
1803   electronV0plusTOFGraphForBBfit->SetName("electronV0plusTOFGraphForBBfit");
1804
1805   for(Int_t ip = 0; ip < electronV0plusTOFGraphForBBfit->GetN(); ip ++) {
1806     Bool_t removePoint = electronV0plusTOFGraphForBBfit->GetX()[ip] < (isPPb ? 1.2 : 0.6);//0.4 || electronV0plusTOFGraphForBBfit->GetX()[ip] > 1.3799;
1807     if (removePoint) {
1808       electronV0plusTOFGraphForBBfit->RemovePoint(ip);
1809       ip--;
1810       continue;
1811     }
1812   }
1813
1814   TGraphErrors * pionV0plusTOFGraphForBBfit = new TGraphErrors(*pionV0plusTOFGraph);
1815   pionV0plusTOFGraphForBBfit->SetName("pionV0plusTOFGraphForBBfit");
1816
1817   for(Int_t ip = 0; ip < pionV0plusTOFGraphForBBfit->GetN(); ip ++) {
1818     Bool_t removePoint = pionV0plusTOFGraphForBBfit->GetX()[ip] < 0.4;// || pionV0plusTOFGraphForBBfit->GetX()[ip] > 6.;
1819     if (removePoint) {
1820       pionV0plusTOFGraphForBBfit->RemovePoint(ip);
1821       ip--;
1822       continue;
1823     }
1824   }
1825
1826   TGraphErrors * protonV0plusTOFGraphForBBfit = new TGraphErrors(*protonV0plusTOFGraph);
1827   protonV0plusTOFGraphForBBfit->SetName("protonV0plusTOFGraphForBBfit");
1828
1829   for(Int_t ip = 0; ip < protonV0plusTOFGraphForBBfit->GetN(); ip ++) {
1830     Bool_t removePoint = protonV0plusTOFGraphForBBfit->GetX()[ip] < 0.9 || protonV0plusTOFGraphForBBfit->GetX()[ip] > 2.5;
1831     if (removePoint) {
1832       protonV0plusTOFGraphForBBfit->RemovePoint(ip);
1833       ip--;
1834       continue;
1835     }
1836   }
1837
1838
1839   //
1840   // rescale graphs to betaGamma
1841   //
1842   kaonTofGraph->SetMarkerStyle(21);
1843   kaonTpcGraph->SetMarkerStyle(22);
1844   for(Int_t ip = 0; ip < kaonTofGraph->GetN(); ip ++) {
1845     kaonTofGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kKaon);
1846     kaonTofGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kKaon);
1847   }
1848   for(Int_t ip = 0; ip < kaonTpcGraph->GetN(); ip ++) {
1849     kaonTpcGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kKaon);
1850     kaonTpcGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kKaon);
1851   }
1852   //
1853   electronGraph->SetMarkerStyle(22);
1854   electronV0Graph->SetMarkerStyle(34);
1855   electronV0GraphForBBfit->SetMarkerStyle(34);
1856   electronV0plusTOFGraph->SetMarkerStyle(30);
1857   electronV0plusTOFGraphForBBfit->SetMarkerStyle(30);
1858   electronGraph->SetMarkerColor(2);
1859   electronV0Graph->SetMarkerColor(2);
1860   electronV0GraphForBBfit->SetMarkerColor(2);
1861   electronV0plusTOFGraph->SetMarkerColor(2);
1862   electronV0plusTOFGraphForBBfit->SetMarkerColor(2);
1863   for(Int_t ip = 0; ip < electronGraph->GetN(); ip ++) {
1864     electronGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kElectron);
1865     electronGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kElectron);
1866   }
1867   for(Int_t ip = 0; ip < electronV0Graph->GetN(); ip ++) {
1868     electronV0Graph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kElectron);
1869     electronV0Graph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kElectron);
1870   }
1871   for(Int_t ip = 0; ip < electronV0GraphForBBfit->GetN(); ip ++) {
1872     electronV0GraphForBBfit->GetX()[ip] /= AliPID::ParticleMass(AliPID::kElectron);
1873     electronV0GraphForBBfit->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kElectron);
1874   }
1875   for(Int_t ip = 0; ip < electronV0plusTOFGraph->GetN(); ip ++) {
1876     electronV0plusTOFGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kElectron);
1877     electronV0plusTOFGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kElectron);
1878   }
1879   for(Int_t ip = 0; ip < electronV0plusTOFGraphForBBfit->GetN(); ip ++) {
1880     electronV0plusTOFGraphForBBfit->GetX()[ip] /= AliPID::ParticleMass(AliPID::kElectron);
1881     electronV0plusTOFGraphForBBfit->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kElectron);
1882   }
1883   //
1884   pionTofGraph->SetMarkerStyle(21);
1885   pionTpcGraph->SetMarkerStyle(22);
1886   pionV0Graph->SetMarkerStyle(34);
1887   pionV0GraphForBBfit->SetMarkerStyle(34);
1888   pionV0plusTOFGraph->SetMarkerStyle(30);
1889   pionV0plusTOFGraphForBBfit->SetMarkerStyle(30);
1890   pionTofGraph->SetMarkerColor(3);
1891   pionTpcGraph->SetMarkerColor(3);
1892   pionV0Graph->SetMarkerColor(3);
1893   pionV0GraphForBBfit->SetMarkerColor(3);
1894   pionV0plusTOFGraph->SetMarkerColor(3);
1895   pionV0plusTOFGraphForBBfit->SetMarkerColor(3);
1896   for(Int_t ip = 0; ip < pionTofGraph->GetN(); ip ++) {
1897     pionTofGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kPion);
1898     pionTofGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kPion);
1899   }
1900   for(Int_t ip = 0; ip < pionTpcGraph->GetN(); ip ++) {
1901     pionTpcGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kPion);
1902     pionTpcGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kPion);
1903   }
1904   for(Int_t ip = 0; ip < pionV0Graph->GetN(); ip ++) {
1905     pionV0Graph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kPion);
1906     pionV0Graph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kPion);
1907   }
1908   for(Int_t ip = 0; ip < pionV0GraphForBBfit->GetN(); ip ++) {
1909     pionV0GraphForBBfit->GetX()[ip] /= AliPID::ParticleMass(AliPID::kPion);
1910     pionV0GraphForBBfit->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kPion);
1911   }
1912   for(Int_t ip = 0; ip < pionV0plusTOFGraph->GetN(); ip ++) {
1913     pionV0plusTOFGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kPion);
1914     pionV0plusTOFGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kPion);
1915   }
1916   for(Int_t ip = 0; ip < pionV0plusTOFGraphForBBfit->GetN(); ip ++) {
1917     pionV0plusTOFGraphForBBfit->GetX()[ip] /= AliPID::ParticleMass(AliPID::kPion);
1918     pionV0plusTOFGraphForBBfit->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kPion);
1919   }
1920   //
1921   protonTofGraph->SetMarkerStyle(21);
1922   protonTpcGraph->SetMarkerStyle(22);
1923   protonV0Graph->SetMarkerStyle(34);
1924   protonV0GraphForBBfit->SetMarkerStyle(34);
1925   protonV0plusTOFGraph->SetMarkerStyle(30);
1926   protonV0plusTOFGraphForBBfit->SetMarkerStyle(30);
1927   protonTofGraph->SetMarkerColor(4);
1928   protonTpcGraph->SetMarkerColor(4);
1929   protonV0Graph->SetMarkerColor(4);
1930   protonV0GraphForBBfit->SetMarkerColor(4);
1931   protonV0plusTOFGraph->SetMarkerColor(4);
1932   protonV0plusTOFGraphForBBfit->SetMarkerColor(4);
1933   for(Int_t ip = 0; ip < protonTofGraph->GetN(); ip ++) {
1934     protonTofGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kProton);
1935     protonTofGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kProton);
1936   }
1937   for(Int_t ip = 0; ip < protonTpcGraph->GetN(); ip ++) {
1938     protonTpcGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kProton);
1939     protonTpcGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kProton);
1940   }
1941   for(Int_t ip = 0; ip < protonV0Graph->GetN(); ip ++) {
1942     protonV0Graph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kProton);
1943     protonV0Graph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kProton);
1944   }
1945   for(Int_t ip = 0; ip < protonV0GraphForBBfit->GetN(); ip ++) {
1946     protonV0GraphForBBfit->GetX()[ip] /= AliPID::ParticleMass(AliPID::kProton);
1947     protonV0GraphForBBfit->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kProton);
1948   }
1949   for(Int_t ip = 0; ip < protonV0plusTOFGraph->GetN(); ip ++) {
1950     protonV0plusTOFGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kProton);
1951     protonV0plusTOFGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kProton);
1952   }
1953   for(Int_t ip = 0; ip < protonV0plusTOFGraphForBBfit->GetN(); ip ++) {
1954     protonV0plusTOFGraphForBBfit->GetX()[ip] /= AliPID::ParticleMass(AliPID::kProton);
1955     protonV0plusTOFGraphForBBfit->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kProton);
1956   }
1957   //
1958   //
1959   TGraphErrors * graphAll = new TGraphErrors();
1960   TList * listColl = new TList();
1961   if (!useV0s) {
1962     //listColl->Add(protonTpcGraph);
1963     //listColl->Add(kaonTpcGraph);
1964     listColl->Add(pionTpcGraph);
1965     listColl->Add(electronGraph);
1966     listColl->Add(protonTofGraph);
1967     listColl->Add(kaonTofGraph);
1968     listColl->Add(pionTofGraph);
1969   }
1970   else {
1971     listColl->Add(pionV0GraphForBBfit);
1972     //listColl->Add(electronV0GraphForBBfit);
1973     listColl->Add(protonV0GraphForBBfit);
1974
1975     //listColl->Add(pionV0plusTOFGraphForBBfit);
1976     listColl->Add(electronV0plusTOFGraphForBBfit);
1977     //listColl->Add(protonV0plusTOFGraphForBBfit);
1978   }
1979   MergeGraphErrors(graphAll, listColl);
1980   graphAll->SetNameTitle("beamDataPoints","beamDataPoints");
1981
1982   // Create graph out of TPC only, V0, V0+TOF and TPC+TOF to determine correction functions
1983   TList * listCollPr = new TList();
1984   if (useV0s) {
1985     listCollPr->Add(protonTpcGraph);
1986     //listCollPr->Add(protonTofGraph);
1987
1988     TGraphErrors * protonV0GraphCut = new TGraphErrors(*protonV0Graph);
1989     protonV0GraphCut->SetName("protonV0GraphCut");
1990     for(Int_t ip = 0; ip < protonV0GraphCut->GetN(); ip ++) {
1991       Double_t mom = protonV0GraphCut->GetX()[ip] * AliPID::ParticleMass(AliPID::kProton);
1992       //Bool_t removePoint = mom >= 0.6 || mom < 0.35; // Will take TPC protons instead (statistics) for very low momenta
1993       Bool_t removePoint = mom < (isPPb ? 0.6 : 0.35); // Will take TPC protons instead (statistics) for very low momenta;
1994                                                        // for pPb low statistics for the moment, therefore, go a bit further with TPC protons
1995
1996
1997       if (removePoint) {
1998         protonV0GraphCut->RemovePoint(ip);
1999         ip--;
2000         continue;
2001       }
2002     }
2003     /*
2004     TGraphErrors * protonV0plusTOFGraphCut = new TGraphErrors(*protonV0plusTOFGraph);
2005     protonV0plusTOFGraphCut->SetName("protonV0plusTOFGraphCut");
2006     for(Int_t ip = 0; ip < protonV0plusTOFGraphCut->GetN(); ip ++) {
2007       Double_t mom = protonV0plusTOFGraphCut->GetX()[ip] * AliPID::ParticleMass(AliPID::kProton);
2008       Bool_t removePoint = mom < 0.6;
2009
2010       if (removePoint) {
2011         protonV0plusTOFGraphCut->RemovePoint(ip);
2012         ip--;
2013         continue;
2014       }
2015     }*/
2016
2017     listCollPr->Add(protonV0GraphCut);
2018     //listCollPr->Add(protonV0plusTOFGraphCut);
2019   }
2020   else {
2021     listCollPr->Add(protonTpcGraph);
2022     listCollPr->Add(protonTofGraph);
2023   }
2024   TGraphErrors * graphPrCombined = new TGraphErrors();
2025   MergeGraphErrors(graphPrCombined, listCollPr);
2026   graphPrCombined->SetMarkerStyle(22);
2027   graphPrCombined->SetMarkerColor(4);
2028   graphPrCombined->SetNameTitle("Graph_Protons_Combined", "Graph_Protons_Combined");
2029
2030   TList * listCollPi = new TList();
2031   if (useV0s) {
2032     //listCollPi->Add(pionTpcGraph);
2033     //listCollPi->Add(pionTofGraph);
2034
2035     TGraphErrors * pionV0GraphCut = new TGraphErrors(*pionV0Graph);
2036     pionV0GraphCut->SetName("pionV0GraphCut");
2037     for(Int_t ip = 0; ip < pionV0GraphCut->GetN(); ip ++) {
2038       //Double_t mom = pionV0GraphCut->GetX()[ip] * AliPID::ParticleMass(AliPID::kPion);
2039       Bool_t removePoint = kFALSE;//mom >= 0.4;
2040                            //|| mom < 0.2;//TODO NOW NEEDED?
2041       if (removePoint) {
2042         pionV0GraphCut->RemovePoint(ip);
2043         ip--;
2044         continue;
2045       }
2046     }
2047     /*
2048     TGraphErrors * pionV0plusTOFGraphCut = new TGraphErrors(*pionV0plusTOFGraph);
2049     pionV0plusTOFGraphCut->SetName("pionV0plusTOFGraphCut");
2050     for(Int_t ip = 0; ip < pionV0plusTOFGraphCut->GetN(); ip ++) {
2051       Double_t mom = pionV0plusTOFGraphCut->GetX()[ip] * AliPID::ParticleMass(AliPID::kPion);
2052       Bool_t removePoint = mom < 0.4;
2053
2054       if (removePoint) {
2055         pionV0plusTOFGraphCut->RemovePoint(ip);
2056         ip--;
2057         continue;
2058       }
2059     }*/
2060
2061     listCollPi->Add(pionV0GraphCut);
2062     //listCollPi->Add(pionV0plusTOFGraphCut);
2063   }
2064   else {
2065     listCollPi->Add(pionTpcGraph);
2066     listCollPi->Add(pionTofGraph);
2067   }
2068   TGraphErrors * graphPiCombined = new TGraphErrors();
2069   MergeGraphErrors(graphPiCombined, listCollPi);
2070   graphPiCombined->SetMarkerStyle(22);
2071   graphPiCombined->SetMarkerColor(3);
2072   graphPiCombined->SetNameTitle("Graph_Pions_Combined", "Graph_Pions_Combined");
2073
2074   TList * listCollKa = new TList();
2075   listCollKa->Add(kaonTpcGraph);
2076   listCollKa->Add(kaonTofGraph);
2077   TGraphErrors * graphKaCombined = new TGraphErrors();
2078   MergeGraphErrors(graphKaCombined, listCollKa);
2079   graphKaCombined->SetMarkerStyle(22);
2080   graphKaCombined->SetMarkerColor(5);
2081   graphKaCombined->SetNameTitle("Graph_Kaons_Combined", "Graph_Kaons_Combined");
2082
2083   TList * listCollEl = new TList();
2084   if (useV0s) {
2085     //listCollEl->Add(electronGraph);
2086
2087     TGraphErrors * electronV0GraphCut = new TGraphErrors(*electronV0Graph);
2088     electronV0GraphCut->SetName("electronV0GraphCut");
2089     /*
2090     for(Int_t ip = 0; ip < electronV0GraphCut->GetN(); ip ++) {
2091       Double_t mom = electronV0GraphCut->GetX()[ip] * AliPID::ParticleMass(AliPID::kElectron);
2092       Bool_t removePoint = mom > 0.3;
2093       if (removePoint) {
2094         electronV0GraphCut->RemovePoint(ip);
2095         ip--;
2096         continue;
2097       }
2098     }
2099
2100     TGraphErrors * electronV0plusTOFGraphCut = new TGraphErrors(*electronV0plusTOFGraph);
2101     electronV0plusTOFGraphCut->SetName("electronV0plusTOFGraphCut");
2102     for(Int_t ip = 0; ip < electronV0plusTOFGraphCut->GetN(); ip ++) {
2103       Double_t mom = electronV0plusTOFGraphCut->GetX()[ip] * AliPID::ParticleMass(AliPID::kElectron);
2104       Bool_t removePoint = mom <= 0.4;
2105
2106       if (removePoint) {
2107         electronV0plusTOFGraphCut->RemovePoint(ip);
2108         ip--;
2109         continue;
2110       }
2111     }
2112     */
2113     listCollEl->Add(electronV0GraphCut);
2114     //listCollEl->Add(electronV0plusTOFGraphCut);
2115   }
2116   else {
2117     listCollEl->Add(electronGraph);
2118   }
2119   TGraphErrors * graphElCombined = new TGraphErrors();
2120   MergeGraphErrors(graphElCombined, listCollEl);
2121   graphElCombined->SetMarkerStyle(22);
2122   graphElCombined->SetMarkerColor(2);
2123   graphElCombined->SetNameTitle("Graph_Electrons_Combined", "Graph_Electrons_Combined");
2124
2125
2126
2127   //
2128   //  return array with all graphs
2129   //
2130   TObjArray * arrGraphs = new TObjArray();
2131   arrGraphs->Add(graphAll);
2132   //
2133   arrGraphs->Add(electronGraph);
2134   arrGraphs->Add(electronV0Graph);
2135   arrGraphs->Add(electronV0plusTOFGraph);
2136   arrGraphs->Add(electronV0GraphForBBfit);
2137   arrGraphs->Add(pionTpcGraph);
2138   arrGraphs->Add(pionTofGraph);
2139   arrGraphs->Add(pionV0Graph);
2140   arrGraphs->Add(pionV0plusTOFGraph);
2141   arrGraphs->Add(pionV0GraphForBBfit);
2142   arrGraphs->Add(kaonTpcGraph);
2143   arrGraphs->Add(kaonTofGraph);
2144   arrGraphs->Add(protonTpcGraph);
2145   arrGraphs->Add(protonTofGraph);
2146   arrGraphs->Add(protonV0Graph);
2147   arrGraphs->Add(protonV0plusTOFGraph);
2148   arrGraphs->Add(protonV0GraphForBBfit);
2149   //
2150   arrGraphs->Add(graphElCombined);
2151   arrGraphs->Add(graphPiCombined);
2152   arrGraphs->Add(graphKaCombined);
2153   arrGraphs->Add(graphPrCombined);
2154   //
2155   //
2156
2157   canvasQAtpc->SaveAs("splines_QA_ResidualGraphsTPC.root");
2158   canvasQAtof->SaveAs("splines_QA_ResidualGraphsTOF.root");
2159   canvasQAv0->SaveAs("splines_QA_ResidualGraphsV0.root");
2160   canvasQAv0plusTOF->SaveAs("splines_QA_ResidualGraphsV0plusTOF.root");
2161
2162
2163   canvasQAv0DeDxPurityEl->SaveAs("V0_dEdx_purity_el.root");
2164   canvasQAv0DeDxPurityPi->SaveAs("V0_dEdx_purity_Pi.root");
2165   canvasQAv0DeDxPurityPr->SaveAs("V0_dEdx_purity_Pr.root");
2166
2167   return arrGraphs;
2168
2169 }
2170
2171
2172 //________________________________________________________________________
2173 TObjArray * AliTPCcalibResidualPID::GetResidualGraphsMC(THnSparseF * histPidQA, const Char_t * /*system*/) {
2174   //
2175   // Extracts the residual graphs from THnSparse created from MC (NOT DATA!)
2176   //
2177
2178   const TString momTitle = "#it{p}_{TPC} (GeV/#it{c})";
2179   const TString dEdxTitle = "d#it{E}/d#it{x} (arb. unit)";
2180
2181   Int_t cutForFitting = 10;
2182   Double_t heightFractionForFittingRange = 0.1;
2183
2184   //
2185   THnSparse * hist = histPidQA;
2186   //
2187   TCanvas * canvasQAmc = new TCanvas("canvasQAmcResGraph","Control canvas for residual graphs (MC)",100,10,1380,800);
2188   canvasQAmc->Divide(2,2);
2189
2190   for (Int_t i = 1; i <= 4; i++) {
2191     canvasQAmc->GetPad(i)->SetGrid(1, 1);
2192     canvasQAmc->GetPad(i)->SetLogz();
2193     canvasQAmc->GetPad(i)->SetLogx();
2194   }
2195   //
2196   // 1. select and fit electrons
2197   //
2198   // 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
2199   hist->GetAxis(2)->SetRangeUser(0,0); // electrons
2200   hist->GetAxis(3)->SetRangeUser(0,0); // electrons (used for nsigma and for the eta correction of the signal)
2201   TH2D * histElectron = hist->Projection(1,0);
2202   histElectron->SetName("histElectron");
2203   histElectron->RebinX(2);
2204   histElectron->GetXaxis()->SetRangeUser(0.15,8.0);
2205   TObjArray arr;
2206   FitSlicesY(histElectron, heightFractionForFittingRange, cutForFitting, "QNRL", &arr);
2207   TH1D * electronPoints = (TH1D *) arr.At(1)->Clone();
2208   //
2209   TGraphErrors * electronGraph = new TGraphErrors(electronPoints);
2210   electronGraph->SetName("electronGraphMC");
2211   for(Int_t ip = 0; ip < electronGraph->GetN(); ip ++) {
2212     Bool_t removePoint = electronGraph->GetY()[ip] < 10 || electronGraph->GetEY()[ip]/electronGraph->GetY()[ip] > 0.03
2213                          || electronGraph->GetX()[ip] < 0.15;// || electronGraph->GetX()[ip] > 3.1;
2214     if (removePoint) {
2215       electronGraph->RemovePoint(ip);
2216       ip--;
2217       continue;
2218     }
2219   }
2220   //
2221   canvasQAmc->cd(1);
2222   histElectron->SetTitle("electrons");
2223   histElectron->GetYaxis()->SetRangeUser(50, 120);
2224   histElectron->GetXaxis()->SetTitle(momTitle.Data());
2225   histElectron->GetYaxis()->SetTitle(dEdxTitle.Data());
2226   histElectron->GetXaxis()->SetMoreLogLabels(kTRUE);
2227   histElectron->GetXaxis()->SetNoExponent(kTRUE);
2228   histElectron->Draw("colz");
2229   electronPoints->SetMarkerStyle(24);
2230   electronPoints->Draw("same");
2231   //
2232   // 2. protons
2233   //
2234   hist->GetAxis(2)->SetRangeUser(3,3); // protons
2235   hist->GetAxis(3)->SetRangeUser(3,3); // protons (used for nsigma and for the eta correction of the signal)
2236   //
2237   TH2D * histProton = hist->Projection(1,0);
2238   histProton->SetName("histProton");
2239   histProton->GetXaxis()->SetRangeUser(0.15,8);
2240   histProton->GetYaxis()->SetRangeUser(10, hist->GetAxis(1)->GetBinUpEdge(hist->GetAxis(1)->GetNbins()));
2241   FitSlicesY(histProton, heightFractionForFittingRange, cutForFitting, "QNRL", &arr);
2242   TH1D * protonPoints = (TH1D *) arr.At(1)->Clone();
2243   //
2244   TGraphErrors * protonGraph = new TGraphErrors(protonPoints);
2245   protonGraph->SetName("protonGraphMC");
2246   for(Int_t ip = 0; ip < protonGraph->GetN(); ip ++) {
2247     Bool_t removePoint = protonGraph->GetY()[ip] < 10 || protonGraph->GetEY()[ip]/protonGraph->GetY()[ip] > 0.02 ||
2248                          protonGraph->GetX()[ip] < 0.15 || protonGraph->GetX()[ip] > 6;
2249     if (removePoint) {
2250       protonGraph->RemovePoint(ip);
2251       ip--;
2252       continue;
2253     }
2254   }
2255   //
2256   canvasQAmc->cd(2);
2257   histProton->SetTitle("protons");
2258   histProton->GetXaxis()->SetTitle(momTitle.Data());
2259   histProton->GetYaxis()->SetTitle(dEdxTitle.Data());
2260   histProton->GetXaxis()->SetMoreLogLabels(kTRUE);
2261   histProton->GetXaxis()->SetNoExponent(kTRUE);
2262   histProton->Draw("colz");
2263   protonPoints->SetMarkerStyle(20);
2264   protonPoints->Draw("same");
2265
2266   //
2267   // 3. pions
2268   //
2269   hist->GetAxis(2)->SetRangeUser(1,1); // pions
2270   hist->GetAxis(3)->SetRangeUser(1,1); // pions (used for nsigma and for the eta correction of the signal)
2271   //
2272   TH2D * histPion = hist->Projection(1,0);
2273   histPion->SetName("histPion");
2274   histPion->GetXaxis()->SetRangeUser(0.15,50);
2275   FitSlicesY(histPion, heightFractionForFittingRange, cutForFitting, "QNRL", &arr);
2276   TH1D * pionPoints =  (TH1D *) arr.At(1)->Clone();
2277   //
2278   TGraphErrors * pionGraph = new TGraphErrors(pionPoints);
2279   pionGraph->SetName("pionGraphMC");
2280   for(Int_t ip = 0; ip < pionGraph->GetN(); ip ++) {
2281     Bool_t removePoint = pionGraph->GetY()[ip] < 10 || pionGraph->GetEY()[ip]/pionGraph->GetY()[ip] > 0.005
2282                          || pionGraph->GetX()[ip] < 0.15 || pionGraph->GetX()[ip] > 30;
2283     if (removePoint) {
2284       pionGraph->RemovePoint(ip);
2285       ip--;
2286       continue;
2287     }
2288   }
2289   //
2290   canvasQAmc->cd(3);
2291   histPion->GetYaxis()->SetRangeUser(30, 90);
2292   histPion->GetXaxis()->SetTitle(momTitle.Data());
2293   histPion->GetYaxis()->SetTitle(dEdxTitle.Data());
2294   histPion->GetXaxis()->SetMoreLogLabels(kTRUE);
2295   histPion->GetXaxis()->SetNoExponent(kTRUE);
2296   histPion->Draw("colz");
2297   histPion->SetTitle("pions");
2298   pionPoints->SetMarkerStyle(20);
2299   pionPoints->Draw("same");
2300
2301   //
2302   // 4. kaons
2303   //
2304   hist->GetAxis(2)->SetRangeUser(2,2); // kaons
2305   hist->GetAxis(3)->SetRangeUser(2,2); // kaons (used for nsigma and for the eta correction of the signal)
2306   //
2307   TH2D * histKaon = hist->Projection(1,0);
2308   histKaon->SetName("histKaon");
2309   histKaon->GetXaxis()->SetRangeUser(0.15,8);
2310   FitSlicesY(histKaon, heightFractionForFittingRange, cutForFitting, "QNRL", &arr);
2311   TH1D * kaonPoints = (TH1D*) arr.At(1)->Clone();
2312   //
2313   TGraphErrors * kaonGraph = new TGraphErrors(kaonPoints);
2314   kaonGraph->SetName("kaonGraphMC");
2315   for(Int_t ip = 0; ip < kaonGraph->GetN(); ip ++) {
2316     Bool_t removePoint = kaonGraph->GetY()[ip] < 10 || kaonGraph->GetEY()[ip]/kaonGraph->GetY()[ip] > 0.02
2317                          || kaonGraph->GetX()[ip] < 0.15;
2318     if (removePoint) {
2319       kaonGraph->RemovePoint(ip);
2320       ip--;
2321       continue;
2322     }
2323   }
2324   //
2325   canvasQAmc->cd(4);
2326   histKaon->SetTitle("kaons");
2327   histKaon->GetYaxis()->SetRangeUser(30, 500);
2328   histKaon->GetXaxis()->SetTitle(momTitle.Data());
2329   histKaon->GetYaxis()->SetTitle(dEdxTitle.Data());
2330   histKaon->GetXaxis()->SetMoreLogLabels(kTRUE);
2331   histKaon->GetXaxis()->SetNoExponent(kTRUE);
2332   histKaon->Draw("colz");
2333   kaonPoints->SetMarkerStyle(20);
2334   kaonPoints->Draw("same");
2335
2336
2337
2338   hist->GetAxis(2)->SetRange(0,-1); // RESET RANGES
2339   hist->GetAxis(3)->SetRange(0,-1); // RESET RANGES
2340
2341
2342   // Clone graphs to have graphs with the same names as for the data case -> same subsequent functions can be used to determine
2343   // the correction functions and, finally, the response functions.
2344   // In addition: Remove some points, which are expected to deviate from the Bethe-Bloch curve, from the original graphs, so that
2345   // these graphs can be used for the BB fit
2346   TGraphErrors * graphPrCombined = new TGraphErrors(*protonGraph);
2347   graphPrCombined->SetMarkerStyle(22);
2348   graphPrCombined->SetMarkerColor(4);
2349   graphPrCombined->SetNameTitle("Graph_Protons_Combined", "Graph_Protons_Combined");
2350
2351   TGraphErrors * graphPiCombined = new TGraphErrors(*pionGraph);
2352   graphPiCombined->SetMarkerStyle(22);
2353   graphPiCombined->SetMarkerColor(3);
2354   graphPiCombined->SetNameTitle("Graph_Pions_Combined", "Graph_Pions_Combined");
2355
2356   TGraphErrors * graphKaCombined = new TGraphErrors(*kaonGraph);
2357   graphKaCombined->SetMarkerStyle(22);
2358   graphKaCombined->SetMarkerColor(5);
2359   graphKaCombined->SetNameTitle("Graph_Kaons_Combined", "Graph_Kaons_Combined");
2360
2361   TGraphErrors * graphElCombined = new TGraphErrors(*electronGraph);
2362   graphElCombined->SetMarkerStyle(22);
2363   graphElCombined->SetMarkerColor(2);
2364   graphElCombined->SetNameTitle("Graph_Electrons_Combined", "Graph_Electrons_Combined");
2365
2366   for(Int_t ip = 0; ip < electronGraph->GetN(); ip ++) {
2367     Bool_t removePoint = electronGraph->GetX()[ip] < 2.5;// || electronGraph->GetX()[ip] > 5.0;
2368     //Bool_t removePoint = electronGraph->GetX()[ip] < 1.2 || electronGraph->GetX()[ip] > 5.0;
2369     if (removePoint) {
2370       electronGraph->RemovePoint(ip);
2371       ip--;
2372       continue;
2373     }
2374   }
2375
2376   for(Int_t ip = 0; ip < protonGraph->GetN(); ip ++) {
2377     Bool_t removePoint = protonGraph->GetX()[ip] < 0.9
2378                          || protonGraph->GetX()[ip] > 2.5;//3.5; //TODO NEW in order not to get into conflict with pions/kaons
2379     if (removePoint) {
2380       protonGraph->RemovePoint(ip);
2381       ip--;
2382       continue;
2383     }
2384   }
2385
2386   for(Int_t ip = 0; ip < pionGraph->GetN(); ip ++) {
2387     Bool_t removePoint = pionGraph->GetX()[ip] < 2.1;//TODO APR18 0.8;
2388     if (removePoint) {
2389       pionGraph->RemovePoint(ip);
2390       ip--;
2391       continue;
2392     }
2393   }
2394
2395   for(Int_t ip = 0; ip < kaonGraph->GetN(); ip ++) {
2396     Bool_t removePoint = kaonGraph->GetX()[ip] < 1.25 || kaonGraph->GetX()[ip] > 4.0;// < 0.7;// || kaonGraph->GetX()[ip] > 4;
2397     if (removePoint) {
2398       kaonGraph->RemovePoint(ip);
2399       ip--;
2400       continue;
2401     }
2402   }
2403   //
2404   // rescale graphs to betaGamma
2405   //
2406   kaonGraph->SetMarkerStyle(22);
2407   for(Int_t ip = 0; ip < kaonGraph->GetN(); ip ++) {
2408     kaonGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kKaon);
2409     kaonGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kKaon);
2410   }
2411   for(Int_t ip = 0; ip < graphKaCombined->GetN(); ip ++) {
2412     graphKaCombined->GetX()[ip] /= AliPID::ParticleMass(AliPID::kKaon);
2413     graphKaCombined->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kKaon);
2414   }
2415   //
2416   electronGraph->SetMarkerStyle(22);
2417   electronGraph->SetMarkerColor(2);
2418   for(Int_t ip = 0; ip < electronGraph->GetN(); ip ++) {
2419     electronGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kElectron);
2420     electronGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kElectron);
2421   }
2422   for(Int_t ip = 0; ip < graphElCombined->GetN(); ip ++) {
2423     graphElCombined->GetX()[ip] /= AliPID::ParticleMass(AliPID::kElectron);
2424     graphElCombined->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kElectron);
2425   }
2426   //
2427   pionGraph->SetMarkerStyle(22);
2428   pionGraph->SetMarkerColor(3);
2429   for(Int_t ip = 0; ip < pionGraph->GetN(); ip ++) {
2430     pionGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kPion);
2431     pionGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kPion);
2432   }
2433   for(Int_t ip = 0; ip < graphPiCombined->GetN(); ip ++) {
2434     graphPiCombined->GetX()[ip] /= AliPID::ParticleMass(AliPID::kPion);
2435     graphPiCombined->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kPion);
2436   }
2437   //
2438   protonGraph->SetMarkerStyle(22);
2439   protonGraph->SetMarkerColor(4);
2440   for(Int_t ip = 0; ip < protonGraph->GetN(); ip ++) {
2441     protonGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kProton);
2442     protonGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kProton);
2443   }
2444   for(Int_t ip = 0; ip < graphPrCombined->GetN(); ip ++) {
2445     graphPrCombined->GetX()[ip] /= AliPID::ParticleMass(AliPID::kProton);
2446     graphPrCombined->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kProton);
2447   }
2448   //
2449   //
2450
2451   // Store graphs without further restriction
2452   TGraphErrors * graphPrAll = new TGraphErrors(*graphPrCombined);
2453   graphPrAll->SetNameTitle("Protons_MC_all", "Protons_MC_all");
2454
2455   TGraphErrors * graphPiAll = new TGraphErrors(*graphPiCombined);
2456   graphPiAll->SetNameTitle("Pions_MC_all", "Pions_MC_all");
2457
2458   TGraphErrors * graphKaAll = new TGraphErrors(*graphKaCombined);
2459   graphKaAll->SetNameTitle("Kaons_MC_all", "Kaons_MC_all");
2460
2461   TGraphErrors * graphElAll = new TGraphErrors(*graphElCombined);
2462   graphElAll->SetNameTitle("Electrons_MC_all", "Electrons_MC_all");
2463
2464   //
2465   //
2466   TGraphErrors * graphAll = new TGraphErrors();
2467   TList * listColl = new TList();
2468   listColl->Add(electronGraph);
2469   listColl->Add(protonGraph);
2470   listColl->Add(kaonGraph);
2471   listColl->Add(pionGraph);
2472   MergeGraphErrors(graphAll, listColl);
2473   graphAll->SetNameTitle("beamDataPoints","beamDataPoints");
2474
2475   //
2476   //  return array with all graphs
2477   //
2478   TObjArray * arrGraphs = new TObjArray();
2479   arrGraphs->Add(graphAll);
2480   //
2481   arrGraphs->Add(electronGraph);
2482   arrGraphs->Add(pionGraph);
2483   arrGraphs->Add(kaonGraph);
2484   arrGraphs->Add(protonGraph);
2485   //
2486   arrGraphs->Add(graphElAll);
2487   arrGraphs->Add(graphPiAll);
2488   arrGraphs->Add(graphKaAll);
2489   arrGraphs->Add(graphPrAll);
2490   //
2491   arrGraphs->Add(graphElCombined);
2492   arrGraphs->Add(graphPiCombined);
2493   arrGraphs->Add(graphKaCombined);
2494   arrGraphs->Add(graphPrCombined);
2495   //
2496   //
2497
2498   canvasQAmc->SaveAs("splines_QA_ResidualGraphsTPC.root");
2499
2500   return arrGraphs;
2501
2502 }
2503
2504
2505 //________________________________________________________________________
2506 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) {
2507   //
2508   //
2509
2510   Bool_t isMC = kFALSE;
2511   if (strcmp(type, "MC") == 0) {
2512     isMC = kTRUE;
2513   }
2514   else if (strcmp(type, "DATA") == 0) {
2515     isMC = kFALSE;
2516   }
2517   else {
2518     Printf("ERROR - GetResponseFunctions: Unknown type \"%s\" - must be \"MC\" or \"DATA\"!", type);
2519
2520     return 0x0;
2521   }
2522
2523   Bool_t isPbPb = kFALSE;
2524   Bool_t isPPb = kFALSE;
2525   if (strcmp(system, "PBPB") == 0) {
2526     Printf("PbPb detected - Applying special handling!");
2527     isPbPb = kTRUE;
2528   }
2529   else if (strcmp(system, "PPB") == 0 || strcmp(system, "PBP") == 0) {
2530     Printf("p-Pb/Pb-p detected - Applying special handling!");
2531     isPPb = kTRUE;
2532   }
2533
2534   TCanvas* canvasQA[4] = { 0x0, };
2535   for (Int_t i = 0; i < 4; i++) {
2536     canvasQA[i] = new TCanvas(Form("canvasQAres_%d", i), "Control canvas for residual polynomial",100,10,1380,800);
2537     canvasQA[i]->SetGrid(1, 1);
2538     canvasQA[i]->SetLogx();
2539   }
2540   /*
2541   TCanvas * canvasQA = new TCanvas("canvasQAres","Control canvas for residual polynomial",100,10,1380,800);
2542   canvasQA->Divide(2,2);
2543
2544   for (Int_t i = 1; i <= 4; i++) {
2545     canvasQA->GetPad(i)->SetGrid(1, 1);
2546   }
2547   */
2548   //
2549   // smallest needed (100MeV proton) bg = 0.1/0.938(=AliPID::ParticleMass(AliPID::kProton)) = 0.1,
2550   // largest needed 100 GeV electron, bg = 100/0.000511(=AliPID::ParticleMass(AliPID::kElectron)) = 200 000
2551   //
2552   Double_t from = 0.1;
2553   Double_t to = 200000;
2554
2555   TF1* funcBB = new TF1(*parametrisation);
2556   funcBB->SetLineColor(2);
2557   funcBB->SetLineWidth(1);
2558
2559   TString fitFuncString = "[0] + [1] / x + [2] / (x**2) + [3] / (x**3) + [4] / (x**4) + [5] / (x**5) + [6] * x";
2560   TString fitFuncString2 = "pol6";
2561   //
2562   // 1. extract proton corrections
2563   //
2564   TGraphErrors *  graphProtonTPC = (TGraphErrors *) inputGraphs->FindObject("Graph_Protons_Combined");
2565   TGraphErrors *  graphProtonTPCfinal = new TGraphErrors(*graphProtonTPC);
2566   graphProtonTPCfinal->SetName("graphProtonTPCfinal");
2567   for(Int_t ip = 0; ip < graphProtonTPC->GetN(); ip ++) {
2568     graphProtonTPC->GetY()[ip] -= funcBB->Eval(graphProtonTPC->GetX()[ip]);
2569     graphProtonTPC->GetY()[ip] /= funcBB->Eval(graphProtonTPC->GetX()[ip]);
2570     graphProtonTPC->GetEY()[ip] /= funcBB->Eval(graphProtonTPC->GetX()[ip]);;
2571   }
2572   canvasQA[0]->cd();
2573   //canvasQA->cd(1);
2574   //gPad->SetLogx();
2575   graphProtonTPC->GetXaxis()->SetTitle("#beta#gamma");
2576   graphProtonTPC->GetXaxis()->SetMoreLogLabels(kTRUE);
2577   graphProtonTPC->GetXaxis()->SetNoExponent(kTRUE);
2578   graphProtonTPC->GetYaxis()->SetLabelSize(0.04);
2579   graphProtonTPC->GetYaxis()->SetTitleOffset(0.85);
2580   graphProtonTPC->GetXaxis()->SetTitleOffset(1.0);
2581   graphProtonTPC->GetYaxis()->SetTitle("(data - fit) / fit");
2582   graphProtonTPC->Draw("ap");
2583
2584   if (graphProtonTPC->GetN() <= 0)  {
2585     Printf("ERROR - GetResponseFunctions: Proton graph has no data points!");
2586
2587     return 0x0;
2588   }
2589   graphProtonTPC->Sort(); // Sort points along x. Next, the very first point will be used to determine the starting point of the correction function
2590   TF1 * funcCorrProton = new TF1("funcCorrProton", fitFuncString2.Data(), TMath::Max(graphProtonTPC->GetX()[0], 0.15),1.0); // TODO BEN was 0.18 - 0.85
2591   graphProtonTPC->Fit(funcCorrProton, "QREX0M");
2592   //
2593   // 2. extract kaon corrections
2594   //
2595   TGraphErrors *  graphKaonTPC = (TGraphErrors *) inputGraphs->FindObject("Graph_Kaons_Combined");
2596   TGraphErrors *  graphKaonTPCfinal = new TGraphErrors(*graphKaonTPC);
2597   graphKaonTPCfinal->SetName("graphKaonTPCfinal");
2598   for(Int_t ip = 0; ip < graphKaonTPC->GetN(); ip ++) {
2599     graphKaonTPC->GetY()[ip] -= funcBB->Eval(graphKaonTPC->GetX()[ip]);
2600     graphKaonTPC->GetY()[ip] /= funcBB->Eval(graphKaonTPC->GetX()[ip]);
2601     graphKaonTPC->GetEY()[ip] /= funcBB->Eval(graphKaonTPC->GetX()[ip]);;
2602   }
2603   canvasQA[1]->cd();
2604   //canvasQA->cd(2);
2605   //gPad->SetLogx();
2606   graphKaonTPC->GetXaxis()->SetTitle("#beta#gamma");
2607   graphKaonTPC->GetYaxis()->SetTitle("(data - fit) / fit");
2608   graphKaonTPC->GetXaxis()->SetMoreLogLabels(kTRUE);
2609   graphKaonTPC->GetXaxis()->SetNoExponent(kTRUE);
2610   graphKaonTPC->GetYaxis()->SetLabelSize(0.04);
2611   graphKaonTPC->GetYaxis()->SetTitleOffset(0.85);
2612   graphKaonTPC->GetXaxis()->SetTitleOffset(1.0);
2613   graphKaonTPC->Draw("ap");
2614
2615   if (graphKaonTPC->GetN() <= 0)  {
2616     Printf("ERROR - GetResponseFunctions: Kaon graph has no data points!");
2617
2618     return 0x0;
2619   }
2620   graphKaonTPC->Sort(); // Sort points along x. Next, the very first point will be used to determine the starting point of the correction function
2621
2622   TF1 * funcCorrKaon = 0x0;
2623
2624   if (isMC) {
2625     funcCorrKaon = new TF1("funcCorrKaon", fitFuncString.Data(), TMath::Max(graphKaonTPC->GetX()[0], 0.25), isPPb ? 3.44 : 1.981);
2626     graphKaonTPC->Fit(funcCorrKaon, "QREX0M");
2627   }
2628   else {
2629     // In case of data there are sometimes problems to fit the shape with one function (could describe the overall deviation,
2630     // but will not cover all the details).
2631     // Nevertheless, this shape (including most of the "details") can be fitted with the following approach with two functions
2632     TF1 * funcCorrKaon1 = new TF1("funcCorrKaon1", fitFuncString.Data(),
2633                                   TMath::Max(graphKaonTPC->GetX()[0], 0.25), 1.981);
2634     graphKaonTPC->Fit(funcCorrKaon1, "QREX0M", "same", TMath::Max(graphKaonTPC->GetX()[0], 0.25), 1.0);
2635
2636     TF1 * funcCorrKaon2 = new TF1("funcCorrKaon2", fitFuncString2.Data(), TMath::Max(graphKaonTPC->GetX()[0], 0.25),  1.981);
2637     graphKaonTPC->Fit(funcCorrKaon2, "QREX0M", "same", (isMC ? 1.981 : 1.0), 1.981);
2638
2639     funcCorrKaon = new TF1("funcCorrKaon",
2640                            "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))",
2641                            TMath::Max(graphKaonTPC->GetX()[0], 0.25), 1.981);
2642
2643     for (Int_t i = funcCorrKaon1->GetNpar(), j = 0; i < funcCorrKaon1->GetNpar() + funcCorrKaon2->GetNpar(); i++, j++) {
2644       funcCorrKaon->SetParameter(j, funcCorrKaon2->GetParameter(j));
2645     }
2646     funcCorrKaon->SetLineColor(kRed);
2647     funcCorrKaon->GetHistogram()->DrawClone("csame");
2648     //funcCorrKaon->Draw("same");
2649   }
2650   /*TODO
2651   TF1 * funcCorrKaon = new TF1("funcCorrKaon", fitFuncString.Data(),//TODO BEN was fitFuncString2
2652                                TMath::Max(graphKaonTPC->GetX()[0], 0.25), (isMC ? 1.981 : 1.0)); //TODO BEN was 0.79 for data and 1.45 for MC
2653   graphKaonTPC->Fit(funcCorrKaon, "QREX0M");
2654   */
2655   //
2656   // 3. extract pion corrections
2657   //
2658   TGraphErrors *  graphPionTPC = (TGraphErrors *)  inputGraphs->FindObject("Graph_Pions_Combined");
2659   TGraphErrors *  graphPionTPCfinal = new TGraphErrors(*graphPionTPC);
2660   graphPionTPCfinal->SetName("graphPionTPCfinal");
2661   for(Int_t ip = 0; ip < graphPionTPC->GetN(); ip ++) {
2662     graphPionTPC->GetY()[ip] -= funcBB->Eval(graphPionTPC->GetX()[ip]);
2663     graphPionTPC->GetY()[ip] /= funcBB->Eval(graphPionTPC->GetX()[ip]);
2664     graphPionTPC->GetEY()[ip] /= funcBB->Eval(graphPionTPC->GetX()[ip]);
2665   }
2666   canvasQA[2]->cd();
2667   //canvasQA->cd(3);
2668   //gPad->SetLogx();
2669   graphPionTPC->GetXaxis()->SetTitle("#beta#gamma");
2670   graphPionTPC->GetYaxis()->SetTitle("(data - fit) / fit");
2671   graphPionTPC->GetXaxis()->SetMoreLogLabels(kTRUE);
2672   graphPionTPC->GetXaxis()->SetNoExponent(kTRUE);
2673   graphPionTPC->GetYaxis()->SetLabelSize(0.04);
2674   graphPionTPC->GetYaxis()->SetTitleOffset(0.85);
2675   graphPionTPC->GetXaxis()->SetTitleOffset(1.0);
2676   graphPionTPC->Draw("ap");
2677
2678   if (graphPionTPC->GetN() <= 0)  {
2679     Printf("ERROR - GetResponseFunctions: Pion graph has no data points!");
2680
2681     return 0x0;
2682   }
2683   graphPionTPC->Sort(); // Sort points along x. Next, the very first point will be used to determine the starting point of the correction function
2684   // In case of PbPb (data, not MC) only correct down to 0.2 GeV/c; otherwise: contamination due to electrons
2685   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));
2686   graphPionTPC->Fit(funcCorrPion, "QREX0M");
2687   //
2688   // 4. extract electron corrections
2689   //
2690   TGraphErrors *  graphElectronTPC = (TGraphErrors *)  inputGraphs->FindObject("Graph_Electrons_Combined");
2691   TGraphErrors *  graphElectronTPCfinal = new TGraphErrors(*graphElectronTPC);
2692   graphElectronTPCfinal->SetName("graphElectronTPCfinal");
2693   for(Int_t ip = 0; ip < graphElectronTPC->GetN(); ip ++) {
2694     graphElectronTPC->GetY()[ip] -= funcBB->Eval(graphElectronTPC->GetX()[ip]);
2695     graphElectronTPC->GetY()[ip] /= funcBB->Eval(graphElectronTPC->GetX()[ip]);
2696     graphElectronTPC->GetEY()[ip] /= funcBB->Eval(graphElectronTPC->GetX()[ip]);;
2697   }
2698   canvasQA[3]->cd();
2699   //canvasQA->cd(4);
2700   //gPad->SetLogx();
2701   graphElectronTPC->GetXaxis()->SetTitle("#beta#gamma");
2702   graphElectronTPC->GetYaxis()->SetTitle("(data - fit) / fit");
2703   graphElectronTPC->GetXaxis()->SetMoreLogLabels(kTRUE);
2704   graphElectronTPC->GetXaxis()->SetNoExponent(kTRUE);
2705   graphElectronTPC->GetYaxis()->SetLabelSize(0.04);
2706   graphElectronTPC->GetYaxis()->SetTitleOffset(0.85);
2707   graphElectronTPC->GetXaxis()->SetTitleOffset(1.0);
2708   graphElectronTPC->Draw("ap");
2709
2710   if (graphElectronTPC->GetN() <= 0)  {
2711     Printf("ERROR - GetResponseFunctions: Electron graph has no data points!");
2712
2713     return 0x0;
2714   }
2715   graphElectronTPC->Sort(); // Sort points along x. Next, the very first point will be used to determine the starting point of the correction function
2716   // In case of PbPb (data, not MC) only correct down to 0.2 GeV/c; otherwise: contamination due to pions
2717   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
2718   // NOTE: For data, the results are almost the same for fitFuncString and fitFuncString2. Maybe, fitFuncString2 is slightly better.
2719   graphElectronTPC->Fit(funcCorrElectron, "QREX0M");
2720   //
2721   // EXTRACT GRAPHS AND PUT THEM TO THE TOBJARRAY
2722   //
2723   const Int_t nBins = 500;
2724   Double_t xBetaGamma[nBins];
2725   Double_t yProton[nBins];
2726   Double_t yKaon[nBins];
2727   Double_t yPion[nBins];
2728   Double_t yElectron[nBins];
2729   Double_t yDefault[nBins];
2730   //
2731   //
2732   //
2733   xBetaGamma[0] = from;
2734   Double_t factor = pow(to/from, 1./nBins);
2735
2736   //
2737   for(Int_t kk = 0; kk < nBins; kk++) {
2738     if (kk > 0) xBetaGamma[kk] = factor * xBetaGamma[kk-1];
2739     yProton[kk] =  funcBB->Eval(xBetaGamma[kk])/50.;
2740     yPion[kk] =  funcBB->Eval(xBetaGamma[kk])/50.;
2741     yKaon[kk] =  funcBB->Eval(xBetaGamma[kk])/50.;
2742     yElectron[kk] =  funcBB->Eval(xBetaGamma[kk])/50.;
2743     yDefault[kk] = funcBB->Eval(xBetaGamma[kk])/50.;
2744
2745     // Added by Ben
2746     Double_t widthFactor = 0.020;
2747     Double_t smoothProton   = 0.5*(TMath::Erf((funcCorrProton->GetXmax()-xBetaGamma[kk])*AliPID::ParticleMass(AliPID::kProton)/widthFactor) + 1);
2748     Double_t smoothKaon     = 0.5*(TMath::Erf((funcCorrKaon->GetXmax()-xBetaGamma[kk])*AliPID::ParticleMass(AliPID::kKaon)/widthFactor) + 1);
2749     Double_t smoothPion     = 0.5*(TMath::Erf((funcCorrPion->GetXmax()-xBetaGamma[kk])*AliPID::ParticleMass(AliPID::kPion)/widthFactor) + 1);
2750     Double_t smoothElectron = 0.5*(TMath::Erf((funcCorrElectron->GetXmax()-xBetaGamma[kk])*AliPID::ParticleMass(AliPID::kElectron)/widthFactor) + 1);
2751
2752     if (xBetaGamma[kk] > funcCorrProton->GetXmax())
2753       yProton[kk] *= (1 + smoothProton*funcCorrProton->Eval(funcCorrProton->GetXmax()));
2754     // Correction is so large that one runs into trouble at low bg for Protons.
2755     // The deviation is smaller if one takes the lower bound of the correction function
2756     else if (xBetaGamma[kk] < funcCorrProton->GetXmin())
2757       yProton[kk] *= (1 + smoothProton*funcCorrProton->Eval(funcCorrProton->GetXmin()));
2758     else
2759       yProton[kk] *= (1 + smoothProton*funcCorrProton->Eval(xBetaGamma[kk]));
2760
2761     if (xBetaGamma[kk] > funcCorrKaon->GetXmax())
2762       yKaon[kk] *= (1 + smoothKaon*funcCorrKaon->Eval(funcCorrKaon->GetXmax()));
2763     // Correction is so large that one runs into trouble at low bg for Kaons.
2764     // The deviation is smaller if one takes the lower bound of the correction function
2765     else if (xBetaGamma[kk] < funcCorrKaon->GetXmin())
2766       yKaon[kk] *= (1 + smoothKaon*funcCorrKaon->Eval(funcCorrKaon->GetXmin()));
2767     else
2768       yKaon[kk] *= (1 + smoothKaon*funcCorrKaon->Eval(xBetaGamma[kk]));
2769
2770     if (xBetaGamma[kk] > funcCorrElectron->GetXmax())
2771       yElectron[kk] *= (1 + smoothElectron*funcCorrElectron->Eval(funcCorrElectron->GetXmax()));
2772     else if (xBetaGamma[kk] < funcCorrElectron->GetXmin())
2773       yElectron[kk] *= (1 + smoothElectron*funcCorrElectron->Eval(funcCorrElectron->GetXmin()));
2774     else
2775       yElectron[kk] *= (1 + smoothElectron*funcCorrElectron->Eval(xBetaGamma[kk]));
2776
2777     // Only true for LHC10d.pass?: Seems not to be needed because any (very small) deviations are most likely due to electron contamination(BEN)
2778     if (xBetaGamma[kk] > funcCorrPion->GetXmax())
2779       yPion[kk] *= (1 + smoothPion*funcCorrPion->Eval(funcCorrPion->GetXmax()));
2780     else if (xBetaGamma[kk] < funcCorrPion->GetXmin())
2781       yPion[kk] *= (1 + smoothPion*funcCorrPion->Eval(funcCorrPion->GetXmin()));
2782     else
2783       yPion[kk] *= (1 + smoothPion*funcCorrPion->Eval(xBetaGamma[kk]));
2784
2785
2786     /* Removed by Ben
2787     Double_t smoothProton = 0.5*(TMath::Erf((funcCorrProton->GetXmax()-xBetaGamma[kk])/0.002) + 1);
2788     Double_t smoothKaon   = 0.5*(TMath::Erf((funcCorrKaon->GetXmax()-xBetaGamma[kk])/0.002) + 1);
2789     Double_t smoothPion   = 0.5*(TMath::Erf((funcCorrPion->GetXmax()-xBetaGamma[kk])/0.002) + 1);
2790
2791     if (xBetaGamma[kk] < funcCorrProton->GetXmax() && xBetaGamma[kk] > funcCorrProton->GetXmin()) yProton[kk] *= (1 + smoothProton*funcCorrProton->Eval(xBetaGamma[kk]));
2792     if (xBetaGamma[kk] < funcCorrKaon->GetXmax() && xBetaGamma[kk] > funcCorrKaon->GetXmin()) yKaon[kk] *= (1 + smoothKaon*funcCorrKaon->Eval(xBetaGamma[kk]));
2793     if (xBetaGamma[kk] < funcCorrPion->GetXmax() && xBetaGamma[kk] > funcCorrPion->GetXmin()) yPion[kk] *= (1 + smoothPion*funcCorrPion->Eval(xBetaGamma[kk]));
2794     //if (xBetaGamma[kk] < funcCorrElectron->GetXmax()) yElectron[kk] *= (1 + funcCorrElectron->Eval(xBetaGamma[kk]));
2795     //
2796     if (xBetaGamma[kk] < funcCorrProton->GetXmin()) yProton[kk] *= (1 + funcCorrProton->Eval(funcCorrProton->GetXmin()));
2797     if (xBetaGamma[kk] < funcCorrKaon->GetXmin()) yKaon[kk] *= (1 + funcCorrKaon->Eval(funcCorrKaon->GetXmin()));
2798     if (xBetaGamma[kk] < funcCorrPion->GetXmin()) yPion[kk] *= (1 + funcCorrPion->Eval(funcCorrPion->GetXmin()));
2799     */
2800   }
2801   //
2802   TGraph * graphProton = new TGraph(nBins, xBetaGamma, yProton);
2803   TGraph * graphElectron = new TGraph(nBins, xBetaGamma, yElectron);
2804   TGraph * graphKaon = new TGraph(nBins, xBetaGamma, yKaon);
2805   TGraph * graphPion = new TGraph(nBins, xBetaGamma, yPion);
2806   TGraph * graphDefault = new TGraph(nBins, xBetaGamma, yDefault);
2807   //
2808   //
2809   TSpline3 * splineProton = new TSpline3("splineProton", graphProton);
2810   TSpline3 * splineElectron = new TSpline3("splineElectron", graphElectron);
2811   TSpline3 * splineKaon = new TSpline3("splineKaon", graphKaon);
2812   TSpline3 * splinePion = new TSpline3("splinePion", graphPion);
2813   TSpline3 * splineDefault = new TSpline3("splineDefault", graphDefault);
2814   //
2815   //
2816   splineProton->SetNameTitle(Form("TSPLINE3_%s_PROTON_%s_%s_%s_%sMEAN",type,period,pass,system,dedxtype),
2817                              Form("TSPLINE3_%s_PROTON_%s_%s_%s_%sMEAN",type,period,pass,system,dedxtype));
2818   splineElectron->SetNameTitle(Form("TSPLINE3_%s_ELECTRON_%s_%s_%s_%sMEAN",type,period,pass,system,dedxtype),
2819                                Form("TSPLINE3_%s_ELECTRON_%s_%s_%s_%sMEAN",type,period,pass,system,dedxtype));
2820   splineKaon->SetNameTitle(Form("TSPLINE3_%s_KAON_%s_%s_%s_%sMEAN",type,period,pass,system,dedxtype),
2821                            Form("TSPLINE3_%s_KAON_%s_%s_%s_%sMEAN",type,period,pass,system,dedxtype));
2822   splinePion->SetNameTitle(Form("TSPLINE3_%s_PION_%s_%s_%s_%sMEAN",type,period,pass,system,dedxtype),
2823                            Form("TSPLINE3_%s_PION_%s_%s_%s_%sMEAN",type,period,pass,system,dedxtype));
2824   splineDefault->SetNameTitle(Form("TSPLINE3_%s_ALL_%s_%s_%s_%sMEAN",type,period,pass,system,dedxtype),
2825                               Form("TSPLINE3_%s_ALL_%s_%s_%s_%sMEAN",type,period,pass,system,dedxtype));
2826   //
2827   TObjArray * arrResponse = new TObjArray();
2828   arrResponse->SetName("TPCpidResponseFunctions");
2829   arrResponse->AddLast(splineProton);
2830   arrResponse->AddLast(splineElectron);
2831   arrResponse->AddLast(splineKaon);
2832   arrResponse->AddLast(splinePion);
2833   arrResponse->AddLast(splineDefault);
2834   //
2835
2836   // Draw deviation from final results
2837   for(Int_t ip = 0; ip < graphProtonTPCfinal->GetN(); ip ++) {
2838     graphProtonTPCfinal->GetY()[ip] -= 50. * splineProton->Eval(graphProtonTPCfinal->GetX()[ip]);
2839     graphProtonTPCfinal->GetY()[ip] /= 50. * splineProton->Eval(graphProtonTPCfinal->GetX()[ip]);
2840     graphProtonTPCfinal->GetEY()[ip] /= 50. * splineProton->Eval(graphProtonTPCfinal->GetX()[ip]);;
2841   }
2842   canvasQA[0]->cd();
2843   //canvasQA->cd(1);
2844   graphProtonTPCfinal->SetMarkerStyle(26);
2845   graphProtonTPCfinal->Draw("psame");
2846
2847   for(Int_t ip = 0; ip < graphKaonTPCfinal->GetN(); ip ++) {
2848     graphKaonTPCfinal->GetY()[ip] -= 50. * splineKaon->Eval(graphKaonTPCfinal->GetX()[ip]);
2849     graphKaonTPCfinal->GetY()[ip] /= 50. * splineKaon->Eval(graphKaonTPCfinal->GetX()[ip]);
2850     graphKaonTPCfinal->GetEY()[ip] /= 50. * splineKaon->Eval(graphKaonTPCfinal->GetX()[ip]);;
2851   }
2852   canvasQA[1]->cd();
2853   //canvasQA->cd(2);
2854   graphKaonTPCfinal->SetMarkerStyle(26);
2855   graphKaonTPCfinal->Draw("psame");
2856
2857   for(Int_t ip = 0; ip < graphPionTPCfinal->GetN(); ip ++) {
2858     graphPionTPCfinal->GetY()[ip] -= 50. * splinePion->Eval(graphPionTPCfinal->GetX()[ip]);
2859     graphPionTPCfinal->GetY()[ip] /= 50. * splinePion->Eval(graphPionTPCfinal->GetX()[ip]);
2860     graphPionTPCfinal->GetEY()[ip] /= 50. * splinePion->Eval(graphPionTPCfinal->GetX()[ip]);;
2861   }
2862   canvasQA[2]->cd();
2863   //canvasQA->cd(3);
2864   graphPionTPCfinal->SetMarkerStyle(26);
2865   graphPionTPCfinal->Draw("psame");
2866
2867   for(Int_t ip = 0; ip < graphElectronTPCfinal->GetN(); ip ++) {
2868     graphElectronTPCfinal->GetY()[ip] -= 50. * splineElectron->Eval(graphElectronTPCfinal->GetX()[ip]);
2869     graphElectronTPCfinal->GetY()[ip] /= 50. * splineElectron->Eval(graphElectronTPCfinal->GetX()[ip]);
2870     graphElectronTPCfinal->GetEY()[ip] /= 50. * splineElectron->Eval(graphElectronTPCfinal->GetX()[ip]);;
2871   }
2872   canvasQA[3]->cd();
2873   //canvasQA->cd(4);
2874   graphElectronTPCfinal->SetMarkerStyle(26);
2875   graphElectronTPCfinal->Draw("psame");
2876
2877   TFile* fSave = TFile::Open("splines_QA_ResidualPolynomials.root", "RECREATE");
2878   fSave->cd();
2879   for (Int_t i = 0; i < 4; i++)
2880     canvasQA[i]->Write();
2881
2882   fSave->Close();
2883   //canvasQA->SaveAs("splines_QA_ResidualPolynomials.root");
2884
2885
2886   delete funcBB;
2887
2888   return arrResponse;
2889 }
2890
2891
2892 //________________________________________________________________________
2893 TF1* AliTPCcalibResidualPID::FitBB(TObjArray* inputGraphs, Bool_t isMC, Bool_t isPPb, const Bool_t useV0s,
2894                                    const Double_t * initialParameters, AliTPCcalibResidualPID::FitType fitType) {
2895   //
2896   // Fit Bethe-Bloch parametrisation to data points (inputGraphs)
2897   //
2898   const Int_t nPar = 6;
2899   TGraphErrors * graphAll = (TGraphErrors *) inputGraphs->FindObject("beamDataPoints");
2900   //
2901   Float_t from = 0.9; //TODO ADJUST -> Very important
2902   Float_t to = graphAll->GetXaxis()->GetXmax() * 2;
2903
2904
2905
2906   TF1* funcBB = 0x0;
2907
2908   Double_t parametersBBForward[nPar] = { 0, };
2909
2910   TVirtualFitter::SetMaxIterations(5e6);
2911
2912   if (fitType == AliTPCcalibResidualPID::kSaturatedLund) {
2913     printf("Fit function: Saturated Lund\n");
2914
2915     funcBB = new TF1("SaturatedLund", SaturatedLund, from, to, nPar);
2916     //Double_t parametersBB[nPar] = {34.0446, 8.42221, 4.16724, 1.29473, 80.6663, 0}; //Xianguos values
2917     //Double_t parametersBB[nPar] = {35.5, 8.7, 2.0, 1.09, 75.6, 0}; // No saturation
2918     Double_t parametersBB[nPar] = {61.0, 8.7, 1.86, 0.85, 113.4, -38}; // Yields reasonable results for data and MC ~ all periods
2919
2920     if (isPPb && !isMC) {
2921       parametersBB[0] = 51.6;
2922       parametersBB[1] = 9.7;
2923       parametersBB[2] = 1.62;
2924       parametersBB[3] = 0.99;
2925       parametersBB[4] = 104.4;
2926       parametersBB[5] = -27.0;
2927     }
2928     if (isMC) {
2929       parametersBB[0] = 41.4;
2930       parametersBB[1] = 8.6;
2931       parametersBB[2] = 2.2;
2932       parametersBB[3] = 0.92;
2933       parametersBB[4] = 90.4;
2934       parametersBB[5] = -20.0;
2935     }
2936     funcBB->SetParameters(parametersBB);
2937
2938     Double_t parameterErrorsBB[nPar] = {5., 0.5, 0.2, 0.05, 10, 10};
2939     funcBB->SetParErrors(parameterErrorsBB);
2940
2941     for (Int_t i = 0; i < nPar; i++)
2942       parametersBBForward[i] = parametersBB[i];
2943   }
2944   else if (fitType == AliTPCcalibResidualPID::kLund) {
2945     printf("Fit function: Lund\n");
2946     printf("****WARNING: Fit settings not tuned for this fit function, only for saturated Lund!\n");
2947
2948     funcBB = new TF1("SaturatedLund", SaturatedLund, from, to, nPar);
2949     //Double_t parametersBB[nPar] = {34.0446, 8.42221, 4.16724, 1.29473, 80.6663, 0}; //Xianguos values
2950     //Double_t parametersBB[nPar] = {35.5, 8.7, 2.0, 1.09, 75.6, 0}; // No saturation
2951     Double_t parametersBB[nPar] = {35.0, 7., 1.86, 1., 75., 0}; // Yields reasonable results for data and MC ~ all periods
2952
2953     funcBB->SetParameters(parametersBB);
2954
2955     Double_t parameterErrorsBB[nPar] = {5., 0.5, 0.2, 0.05, 10, 10};
2956     funcBB->SetParErrors(parameterErrorsBB);
2957
2958     funcBB->FixParameter(5, 0.0); // No saturation
2959
2960     for (Int_t i = 0; i < nPar; i++)
2961       parametersBBForward[i] = parametersBB[i];
2962   }
2963   else if (fitType == kAlephWithAdditionalParam) {
2964     printf("Fit function: Aleph with additional parameter\n");
2965     printf("****WARNING: Fit settings not tuned for this fit function, only for saturated Lund!\n");
2966
2967     // NOTE: This form is equivalent to the original form, but with parameter [2] redefined for better numerical stability.
2968     // The additional parameter [5] has been introduced later and is unity originally. It seems not to be needed and is, thus,
2969     // fixed to unity
2970     funcBB = new TF1("funcAleph",
2971                      "[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);
2972                      //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);
2973
2974     Double_t parametersBB[nPar] = {2.6,14.3,-15,2.2,2.7, 0.06};
2975     //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};
2976     funcBB->SetParameters(parametersBB);
2977
2978     for (Int_t i = 0; i < nPar; i++)
2979       parametersBBForward[i] = parametersBB[i];
2980   }
2981   else if (fitType == AliTPCcalibResidualPID::kAleph) {
2982     printf("Fit function: Aleph\n");
2983     printf("****WARNING: Fit settings not tuned for this fit function, only for saturated Lund!\n");
2984
2985     // NOTE: This form is equivalent to the original form, but with parameter [2] redefined for better numerical stability.
2986     // The additional parameter [5] has been introduced later and is unity originally. It seems not to be needed and is, thus,
2987     // fixed to unity
2988     funcBB = new TF1("funcAleph",
2989                      "[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);
2990                      //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);
2991     //TEST Double_t parametersBB[nPar] = {1.2, 26.0, -30.0, -2.15, 5.6, 1};
2992     Double_t parametersBB[nPar] = {1.25, 27.5, -29.0, 2.2, 5.2, 1};
2993
2994     // For [5] floating: Double_t parametersBB[nPar] = {2.42,15.2,-16,-2.24,2.8, 0.057};
2995     //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};
2996     //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};
2997     funcBB->SetParameters(parametersBB);
2998
2999     // Fix parameter 5 to original value of unity
3000     funcBB->FixParameter(5, 1);
3001
3002     for (Int_t i = 0; i < nPar; i++)
3003       parametersBBForward[i] = parametersBB[i];
3004   }
3005   else {
3006     printf("Error: fitType not supported!\n");
3007     return 0x0;
3008   }
3009
3010   funcBB->SetLineColor(2);
3011   funcBB->SetLineWidth(1);
3012
3013   // Override initial parameters, if user provides some
3014   if (initialParameters) {
3015     printf("Setting user initial parameters!\n");
3016     funcBB->SetParameters(initialParameters);
3017   }
3018
3019   //
3020   //
3021   //
3022   TGraphErrors * graphDelta = new TGraphErrors(*graphAll);
3023   graphDelta->SetName("graphDelta");
3024
3025   // In MC case: Remove errors from fit -> Some data points with extremely small errors otherwise completely dominate the fit,
3026   // but these points could still be wrong due to systematics (low p effects, e.g.)
3027   if (isMC) {
3028     for(Int_t ip = 0; ip < graphAll->GetN(); ip++) {
3029       graphAll->SetPointError(ip, 0, 0);
3030     }
3031   }
3032   else if (useV0s) {
3033     // Increase weight in plateau by reducing errors. Neccessary because errors are large there compared to other data points
3034     for(Int_t ip = 0; ip < graphAll->GetN(); ip++) {
3035       if (graphAll->GetX()[ip] >= 2500) {
3036         graphAll->SetPointError(ip, 0, graphAll->GetEY()[ip] / 6.0);
3037       }
3038       // Same for pions in the very rel. rise
3039       if (graphAll->GetX()[ip] >= 25 && graphAll->GetX()[ip] < 1000) {
3040         graphAll->SetPointError(ip, 0, graphAll->GetEY()[ip] / 6.0);
3041       }
3042       // Same for protons in the MIP region
3043       if (graphAll->GetX()[ip] >= 2 && graphAll->GetX()[ip] < 4) {
3044         graphAll->SetPointError(ip, 0, graphAll->GetEY()[ip] / 6.0);
3045       }
3046     }
3047   }
3048   graphAll->Fit(funcBB, "REX0M");
3049   funcBB->SetRange(from, to);
3050   funcBB->GetParameters(parametersBBForward);
3051   //
3052   //
3053   //
3054   for(Int_t ip = 0; ip < graphDelta->GetN(); ip++) {
3055     graphDelta->GetY()[ip] -= funcBB->Eval(graphDelta->GetX()[ip]);
3056     graphDelta->GetY()[ip] /= funcBB->Eval(graphDelta->GetX()[ip]);
3057     graphDelta->GetEY()[ip] /= funcBB->Eval(graphDelta->GetX()[ip]);
3058   }
3059   TCanvas * canvDelta_1 = new TCanvas("canvDelta_1","control histogram for Bethe-Bloch fit 1",100,10,1380,800);
3060   TCanvas * canvDelta_2 = new TCanvas("canvDelta_2","control histogram for Bethe-Bloch fit 2",100,10,1380,800);
3061   canvDelta_1->SetGrid(1, 1);
3062   canvDelta_1->SetLogx();
3063   canvDelta_2->SetGrid(1, 1);
3064   canvDelta_2->SetLogx();
3065   /*
3066   canvDelta->Divide(2, 1);
3067   canvDelta->GetPad(1)->SetGrid(1, 1);
3068   canvDelta->GetPad(1)->SetLogx();
3069   canvDelta->GetPad(2)->SetGrid(1, 1);
3070   canvDelta->GetPad(2)->SetLogx();
3071   */
3072
3073   canvDelta_1->cd();
3074   //canvDelta->cd(1);
3075   TH1F *hBBdummy=new TH1F("hBBdummy","BB fit;#beta#gamma;#LTdE/dx#GT (arb. unit)",100,.8,1e4);
3076   hBBdummy->SetMinimum(45);
3077   hBBdummy->SetMaximum(120);
3078   hBBdummy->GetXaxis()->SetTitleOffset(1.1);
3079   hBBdummy->SetStats(kFALSE);
3080   hBBdummy->Draw();
3081
3082   graphAll->SetTitle("BB multi-graph fit");
3083   graphAll->SetMarkerStyle(22);
3084   graphAll->SetMarkerColor(kMagenta);
3085   graphAll->Draw("p");
3086
3087   TLegend *leg=new TLegend(.7,.12,.89,isMC?.4:.6);
3088   leg->SetFillColor(10);
3089   leg->SetBorderSize(1);
3090
3091   //overlay individual graphs
3092   TGraph *gr=0x0;
3093   if (isMC) {
3094     gr=(TGraph*)inputGraphs->FindObject("Protons_MC_all");
3095     gr->SetMarkerColor(kRed);
3096     gr->SetMarkerStyle(24);
3097     gr->Draw("p");
3098     leg->AddEntry(gr,"p (MC)","p");
3099     gr=(TGraph*)inputGraphs->FindObject("Pions_MC_all");
3100     gr->SetMarkerColor(kBlue);
3101     gr->SetMarkerStyle(24);
3102     gr->Draw("p");
3103     leg->AddEntry(gr,"#pi (MC)","p");
3104     gr=(TGraph*)inputGraphs->FindObject("Kaons_MC_all");
3105     gr->SetMarkerColor(kGreen);
3106     gr->SetMarkerStyle(24);
3107     gr->Draw("p");
3108     leg->AddEntry(gr,"K (MC)","p");
3109     gr=(TGraph*)inputGraphs->FindObject("Electrons_MC_all");
3110     gr->SetMarkerColor(kBlack);
3111     gr->SetMarkerStyle(24);
3112     gr->Draw("p");
3113     leg->AddEntry(gr,"e (MC)","p");
3114   }
3115   else {
3116     gr=(TGraph*)inputGraphs->FindObject("protonTpcGraph");
3117     gr->SetMarkerColor(kRed);
3118     gr->SetMarkerStyle(26);
3119     gr->Draw("p");
3120     leg->AddEntry(gr,"p (TPC)","p");
3121     gr=(TGraph*)inputGraphs->FindObject("protonTofGraph");
3122     gr->SetMarkerColor(kRed);
3123     gr->SetMarkerStyle(25);
3124     gr->Draw("p");
3125     leg->AddEntry(gr,"p (TOF)","p");
3126     gr=(TGraph*)inputGraphs->FindObject("protonV0Graph");//ForBBfit");
3127     gr->SetMarkerColor(kRed);
3128     gr->SetMarkerStyle(24);
3129     gr->Draw("p");
3130     leg->AddEntry(gr,"p (V0)","p");
3131     gr=(TGraph*)inputGraphs->FindObject("protonV0plusTOFGraph");//ForBBfit");
3132     gr->SetMarkerColor(kRed);
3133     gr->SetMarkerStyle(30);
3134     gr->Draw("p");
3135     leg->AddEntry(gr,"p (V0+TOF)","p");
3136     gr=(TGraph*)inputGraphs->FindObject("pionTpcGraph");
3137     gr->SetMarkerColor(kBlue);
3138     gr->SetMarkerStyle(26);
3139     gr->Draw("p");
3140     leg->AddEntry(gr,"#pi (TPC)","p");
3141     gr=(TGraph*)inputGraphs->FindObject("pionTofGraph");
3142     gr->SetMarkerColor(kBlue);
3143     gr->SetMarkerStyle(25);
3144     gr->Draw("p");
3145     leg->AddEntry(gr,"#pi (TOF)","p");
3146     gr=(TGraph*)inputGraphs->FindObject("pionV0Graph");//ForBBfit");
3147     gr->SetMarkerColor(kBlue);
3148     gr->SetMarkerStyle(24);
3149     gr->Draw("p");
3150     leg->AddEntry(gr,"#pi (V0)","p");
3151     gr=(TGraph*)inputGraphs->FindObject("pionV0plusTOFGraph");//ForBBfit");
3152     gr->SetMarkerColor(kBlue);
3153     gr->SetMarkerStyle(30);
3154     gr->Draw("p");
3155     leg->AddEntry(gr,"#pi (V0+TOF)","p");
3156     gr=(TGraph*)inputGraphs->FindObject("kaonTpcGraph");
3157     gr->SetMarkerColor(kGreen);
3158     gr->SetMarkerStyle(26);
3159     gr->Draw("p");
3160     leg->AddEntry(gr,"K (TPC)","p");
3161     gr=(TGraph*)inputGraphs->FindObject("kaonTofGraph");
3162     gr->SetMarkerColor(kGreen);
3163     gr->SetMarkerStyle(25);
3164     gr->Draw("p");
3165     leg->AddEntry(gr,"K (TOF)","p");
3166     gr=(TGraph*)inputGraphs->FindObject("electronGraph");
3167     gr->SetMarkerColor(kBlack);
3168     gr->SetMarkerStyle(25);
3169     gr->Draw("p");
3170     leg->AddEntry(gr,"e","p");
3171     gr=(TGraph*)inputGraphs->FindObject("electronV0Graph");//ForBBfit");
3172     gr->SetMarkerColor(kBlack);
3173     gr->SetMarkerStyle(24);
3174     gr->Draw("p");
3175     leg->AddEntry(gr,"e (V0)","p");
3176     gr=(TGraph*)inputGraphs->FindObject("electronV0plusTOFGraph");//ForBBfit");
3177     gr->SetMarkerColor(kBlack);
3178     gr->SetMarkerStyle(30);
3179     gr->Draw("p");
3180     leg->AddEntry(gr,"e (V0+TOF)","p");
3181   }
3182
3183   graphAll->GetXaxis()->SetTitle("#beta#gamma");
3184   graphAll->GetYaxis()->SetTitle("<dE/dx> (arb. unit)");
3185   graphAll->Draw("p");
3186   funcBB->GetHistogram()->DrawClone("csame");
3187
3188   leg->AddEntry(graphAll,"used","p");
3189   leg->AddEntry(funcBB,"fit","l");
3190   leg->Draw("same");
3191
3192   canvDelta_2->cd();
3193   //canvDelta->cd(2);
3194   TH1F *hBBresdummy=new TH1F("hBBresdummy","residuals of BB fit;#beta#gamma;(data-fit)/fit",100,.8,1e4);
3195   hBBresdummy->SetMinimum(-0.04);
3196   hBBresdummy->SetMaximum(0.04);
3197   hBBresdummy->GetXaxis()->SetTitleOffset(1.1);
3198   hBBresdummy->GetYaxis()->SetTitleOffset(0.8);
3199   hBBresdummy->SetStats(kFALSE);
3200   hBBresdummy->Draw();
3201
3202   graphDelta->SetTitle("residuals of BB fit to multi-graph");
3203   graphDelta->GetXaxis()->SetTitle("#beta#gamma");
3204   graphDelta->GetYaxis()->SetTitle("(data - fit) / fit");
3205   graphDelta->SetMarkerStyle(22);
3206   graphDelta->SetMarkerColor(4);
3207   graphDelta->Draw("p");
3208
3209   TFile* fSave = TFile::Open("splines_QA_BetheBlochFit.root", "RECREATE");
3210   fSave->cd();
3211   canvDelta_1->Write();
3212   canvDelta_2->Write();
3213
3214   TString fitResults = "Fit results:\n";
3215   for (Int_t i = 0; i < nPar; i++) {
3216     fitResults.Append(Form("par%d:\t%f +- %f\n", i, funcBB->GetParameters()[i], funcBB->GetParErrors()[i]));
3217   }
3218
3219   TNamed* settings = new TNamed(fitResults.Data(), fitResults.Data());
3220   settings->Write();
3221   fSave->Close();
3222
3223
3224   printf("\n\n%s", fitResults.Data());
3225
3226   return funcBB;
3227 }
3228
3229
3230 //________________________________________________________________________
3231 Double_t AliTPCcalibResidualPID::Lund(Double_t* xx, Double_t* par)
3232 {
3233   // bg is beta*gamma
3234   const Double_t bg = xx[0];
3235
3236   const Double_t beta2 = bg*bg / (1.0 + bg*bg);
3237
3238   const Double_t a = par[0];
3239   const Double_t b = par[1];
3240   const Double_t c = par[2];
3241   const Double_t e = par[3];
3242   const Double_t f = par[4];
3243
3244   const Double_t d = TMath::Exp(c*(a - f) / b);
3245
3246   const Double_t powbg = TMath::Power(1.0 + bg, c);
3247
3248   const Double_t value = a / TMath::Power(beta2,e) +
3249     b/c * TMath::Log(powbg / (1.0 + d*powbg));
3250
3251   return value;
3252 }
3253
3254
3255 //________________________________________________________________________
3256 Double_t AliTPCcalibResidualPID::SaturatedLund(Double_t* xx, Double_t* par)
3257 {
3258   const Double_t qq = Lund(xx, par);
3259   return qq * TMath::Exp(par[5] / qq);
3260 }
3261
3262
3263 //________________________________________________________________________
3264 void AliTPCcalibResidualPID::FitSlicesY(TH2 *hist, Double_t heightFractionForRange, Int_t cutThreshold, TString fitOption, TObjArray *arr)
3265 {
3266   //heightPercentageForRange
3267   // custom slices fit
3268   //
3269
3270   if (!hist) return;
3271   if (!arr) return;
3272
3273   // If old style is to be used
3274   /*
3275   hist->FitSlicesY(0, 0, -1, cutThreshold, fitOption.Data(), &arr);
3276   return;
3277   */
3278
3279
3280   arr->Clear();
3281   arr->SetOwner();
3282   arr->Expand(4);
3283
3284   TAxis *axis=hist->GetXaxis();
3285   const TArrayD *bins = axis->GetXbins();
3286   TH1D** hList = new TH1D*[4];
3287
3288   for (Int_t i = 0; i < 4; i++) {
3289     delete gDirectory->FindObject(Form("%s_%d", hist->GetName(), i));
3290
3291     if (bins->fN == 0) {
3292       hList[i] = new TH1D(Form("%s_%d", hist->GetName(), i), i < 3 ? Form("Fitted value of par[%d]", i) : "Chi2/NDF",
3293                           hist->GetNbinsX(), axis->GetXmin(), axis->GetXmax());
3294     } else {
3295       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);
3296     }
3297
3298     (*arr)[i] = hList[i];
3299   }
3300
3301   for (Int_t ibin=axis->GetFirst(); ibin<=axis->GetLast(); ++ibin){
3302     TH1 *h=hist->ProjectionY("_temp",ibin,ibin);
3303     if (!h)
3304       continue;
3305
3306     if (h->GetEntries() < cutThreshold) {
3307       delete h;
3308       continue;
3309     }
3310
3311     // Average around maximum bin -> Might catch some outliers
3312     Int_t maxBin = h->GetMaximumBin();
3313     Double_t maxVal = h->GetBinContent(maxBin);
3314
3315     if (maxVal < 2) { // It could happen that all entries are in overflow/underflow; don't fit in this case
3316       delete h;
3317       continue;
3318     }
3319
3320     UChar_t usedBins = 1;
3321     if (maxBin > 1) {
3322       maxVal += h->GetBinContent(maxBin - 1);
3323       usedBins++;
3324     }
3325     if (maxBin < h->GetNbinsX()) {
3326       maxVal += h->GetBinContent(maxBin + 1);
3327       usedBins++;
3328     }
3329     maxVal /= usedBins;
3330
3331     Double_t thresholdFraction = heightFractionForRange * maxVal;
3332     Int_t lowThrBin = TMath::Max(1, h->FindFirstBinAbove(thresholdFraction));
3333     Int_t highThrBin = TMath::Min(h->GetNbinsX(), h->FindLastBinAbove(thresholdFraction));
3334
3335     Double_t lowThreshold = h->GetBinCenter(lowThrBin);
3336     Double_t highThreshold = h->GetBinCenter(highThrBin);
3337
3338     TFitResultPtr res = h->Fit("gaus", Form("%sS", fitOption.Data()), "", lowThreshold, highThreshold);
3339
3340     if ((Int_t)res == 0) {
3341       Int_t resBin = ibin;
3342       hList[0]->SetBinContent(resBin,res->GetParams()[0]);
3343       hList[0]->SetBinError(resBin,res->GetErrors()[0]);
3344       hList[1]->SetBinContent(resBin,res->GetParams()[1]);
3345       hList[1]->SetBinError(resBin,res->GetErrors()[1]);
3346       hList[2]->SetBinContent(resBin,res->GetParams()[2]);
3347       hList[2]->SetBinError(resBin,res->GetErrors()[2]);
3348       hList[3]->SetBinContent(resBin,res->Ndf()>0?res->Chi2()/res->Ndf():0);
3349     }
3350
3351     delete h;
3352   }
3353 }
3354
3355
3356 //______________________________________________________________________________
3357 Bool_t AliTPCcalibResidualPID::GetVertexIsOk(AliVEvent* event) const
3358 {
3359   AliAODEvent* aod = 0x0;
3360   AliESDEvent* esd = 0x0;
3361
3362   aod = dynamic_cast<AliAODEvent*>(event);
3363   if (!aod) {
3364     esd = dynamic_cast<AliESDEvent*>(event);
3365
3366     if (!esd) {
3367       AliError("Event seems to be neither AOD nor ESD!");
3368       return kFALSE;
3369     }
3370   }
3371
3372   if (fIsPbpOrpPb) {
3373     const AliVVertex* trkVtx = (aod ? dynamic_cast<const AliVVertex*>(aod->GetPrimaryVertex()) :
3374                                       dynamic_cast<const AliVVertex*>(esd->GetPrimaryVertex()));
3375
3376     if (!trkVtx || trkVtx->GetNContributors() <= 0)
3377       return kFALSE;
3378
3379     TString vtxTtl = trkVtx->GetTitle();
3380     if (!vtxTtl.Contains("VertexerTracks"))
3381       return kFALSE;
3382
3383     Float_t zvtx = trkVtx->GetZ();
3384     const AliVVertex* spdVtx = (aod ? dynamic_cast<const AliVVertex*>(aod->GetPrimaryVertexSPD()) :
3385                                       dynamic_cast<const AliVVertex*>(esd->GetPrimaryVertexSPD()));
3386     if (spdVtx->GetNContributors() <= 0)
3387       return kFALSE;
3388
3389     TString vtxTyp = spdVtx->GetTitle();
3390     Double_t cov[6] = {0};
3391     spdVtx->GetCovarianceMatrix(cov);
3392     Double_t zRes = TMath::Sqrt(cov[5]);
3393     if (vtxTyp.Contains("vertexer:Z") && (zRes > 0.25))
3394       return kFALSE;
3395
3396     if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ()) > 0.5)
3397       return kFALSE;
3398
3399     if (TMath::Abs(zvtx) > fZvtxCutEvent) //Default: 10 cm
3400       return kFALSE;
3401
3402     return kTRUE;
3403   }
3404
3405
3406   // pp and PbPb
3407   const AliVVertex* primaryVertex = (aod ? dynamic_cast<const AliVVertex*>(aod->GetPrimaryVertex()) :
3408                                            dynamic_cast<const AliVVertex*>(esd->GetPrimaryVertexTracks()));
3409
3410   if (!primaryVertex || primaryVertex->GetNContributors() <= 0)
3411     return kFALSE;
3412
3413   if (TMath::Abs(primaryVertex->GetZ()) >= fZvtxCutEvent) //Default: 10 cm
3414     return kFALSE;
3415
3416   return kTRUE;
3417 }
3418
3419
3420 //______________________________________________________________________________
3421 void AliTPCcalibResidualPID::FillV0PIDlist(AliESDEvent* event)
3422 {
3423   //
3424   // Fill the PID tag list
3425   //
3426
3427   // If no event forwarded as parameter (default), cast current input event.
3428   // Dynamic cast to ESD events (DO NOTHING for AOD events)
3429   if (!event)
3430     event = dynamic_cast<AliESDEvent *>(InputEvent());
3431
3432   // If this fails, do nothing
3433   if (!event)
3434     return;
3435
3436   if (!fV0KineCuts) {
3437     AliError("V0KineCuts not available!");
3438     return;
3439   }
3440
3441   TString beamType(event->GetBeamType());
3442
3443   if (beamType.CompareTo("Pb-Pb") == 0 || beamType.CompareTo("A-A") == 0) {
3444     fV0KineCuts->SetMode(AliESDv0KineCuts::kPurity, AliESDv0KineCuts::kPbPb);
3445   }
3446   else {
3447     fV0KineCuts->SetMode(AliESDv0KineCuts::kPurity, AliESDv0KineCuts::kPP);
3448   }
3449
3450   // V0 selection
3451   // set event
3452   fV0KineCuts->SetEvent(event);
3453
3454   const Int_t numTracks = event->GetNumberOfTracks();
3455   fV0tags = new Char_t[numTracks];
3456   for (Int_t i = 0; i < numTracks; i++)
3457     fV0tags[i] = 0;
3458
3459   fV0motherIndex = new Int_t[numTracks];
3460   for (Int_t i = 0; i < numTracks; i++)
3461     fV0motherIndex[i] = -1;
3462
3463   fV0motherPDG = new Int_t[numTracks];
3464   for (Int_t i = 0; i < numTracks; i++)
3465     fV0motherPDG[i] = 0;
3466
3467   fNumTagsStored = numTracks;
3468
3469
3470
3471
3472
3473   // loop over V0 particles
3474   for (Int_t iv0 = 0; iv0 < event->GetNumberOfV0s(); iv0++) {
3475     AliESDv0* v0 = (AliESDv0*)event->GetV0(iv0);
3476
3477     if (!v0)
3478       continue;
3479
3480     // Reject onFly V0's <-> Only take V0's from offline V0 finder
3481     if (v0->GetOnFlyStatus())
3482       continue;
3483
3484     // Get the particle selection
3485     Bool_t foundV0 = kFALSE;
3486     Int_t pdgV0 = 0, pdgP = 0, pdgN = 0;
3487
3488     foundV0 = fV0KineCuts->ProcessV0(v0, pdgV0, pdgP, pdgN);
3489     if (!foundV0)
3490       continue;
3491
3492     Int_t iTrackP = v0->GetPindex();  // positive track
3493     Int_t iTrackN = v0->GetNindex();  // negative track
3494
3495     /*
3496     AliESDtrack* trackP = event->GetTrack(iTrackP);
3497     AliESDtrack* trackN = event->GetTrack(iTrackN);
3498
3499     if (!trackP || !trackN)
3500       continue;
3501
3502
3503
3504     Float_t xy = 999, z = 999;
3505     trackP->GetImpactParameters(xy, z);
3506
3507     Bool_t reject = kFALSE;
3508     if (TMath::Abs(z) > 1)
3509       continue;
3510
3511     trackN->GetImpactParameters(xy, z);
3512     if (TMath::Abs(z) > 1)
3513       continue;
3514     */
3515
3516
3517     // Fill the Object arrays
3518     // positive particles
3519     if (pdgP == -11) {
3520       fV0tags[iTrackP] = 14;
3521     }
3522     else if (pdgP == 211) {
3523       fV0tags[iTrackP] = 15;
3524     }
3525     else if(pdgP == 2212) {
3526       fV0tags[iTrackP] = 16;
3527     }
3528
3529     fV0motherIndex[iTrackP] = iv0;
3530     fV0motherPDG[iTrackP] = pdgV0;
3531
3532     // negative particles
3533     if( pdgN == 11){
3534       fV0tags[iTrackN] = -14;
3535     }
3536     else if( pdgN == -211){
3537       fV0tags[iTrackN] = -15;
3538     }
3539     else if( pdgN == -2212){
3540       fV0tags[iTrackN] = -16;
3541     }
3542
3543     fV0motherIndex[iTrackN] = iv0;
3544     fV0motherPDG[iTrackN] = pdgV0;
3545
3546   }
3547 }
3548
3549
3550 //______________________________________________________________________________
3551 void AliTPCcalibResidualPID::ClearV0PIDlist()
3552 {
3553   //
3554   // Clear the PID tag list
3555   //
3556
3557   delete fV0tags;
3558   fV0tags = 0;
3559
3560   delete fV0motherIndex;
3561   fV0motherIndex = 0;
3562
3563   delete fV0motherPDG;
3564   fV0motherPDG = 0;
3565
3566   fNumTagsStored = 0;
3567 }
3568
3569
3570 //______________________________________________________________________________
3571 Char_t AliTPCcalibResidualPID::GetV0tag(Int_t trackIndex) const
3572 {
3573   //
3574   // Get the tag for the corresponding trackIndex. Returns -99 in case of invalid index/tag list.
3575   //
3576
3577   if (trackIndex < 0 || trackIndex >= fNumTagsStored || !fV0tags)
3578     return -99;
3579
3580   return fV0tags[trackIndex];
3581 }
3582
3583
3584 //______________________________________________________________________________
3585 Int_t AliTPCcalibResidualPID::GetV0motherIndex(Int_t trackIndex) const
3586 {
3587   //
3588   // Get the index of the V0 mother for the corresponding trackIndex. Returns -99 in case of invalid index/mother index list.
3589   //
3590
3591   if (trackIndex < 0 || trackIndex >= fNumTagsStored || !fV0motherIndex)
3592     return -99;
3593
3594   return fV0motherIndex[trackIndex];
3595 }
3596
3597
3598 //______________________________________________________________________________
3599 Int_t AliTPCcalibResidualPID::GetV0motherPDG(Int_t trackIndex) const
3600 {
3601   //
3602   // Get the PDG code of the V0 mother for the corresponding trackIndex. Returns 0 in case of invalid index/mother index list.
3603   //
3604
3605   if (trackIndex < 0 || trackIndex >= fNumTagsStored || !fV0motherPDG)
3606     return 0;
3607
3608   return fV0motherPDG[trackIndex];
3609 }
3610
3611
3612 //______________________________________________________________________________
3613 Int_t AliTPCcalibResidualPID::MergeGraphErrors(TGraphErrors* mergedGraph, TCollection* li)
3614 {
3615   // Adds all graphs with errors from the collection to this graph.
3616   // Returns the total number of poins in the result or -1 in case of an error.
3617
3618   // TODO "Normal" merge of latest root trunk will also take into account the error bars,
3619   // so this function can be replaced by the normal root function with the latest root version.
3620
3621   TIter next(li);
3622   while (TObject* o = next()) {
3623     TGraph *g = dynamic_cast<TGraph*>(o);
3624     if (!g) {
3625       Printf("ERROR: Cannot merge an object which doesn't inherit from TGraph found in the list");
3626       return -1;
3627     }
3628     int n0 = mergedGraph->GetN();
3629     int n1 = n0+g->GetN();
3630     mergedGraph->Set(n1);
3631     Double_t * x = g->GetX();
3632     Double_t * y = g->GetY();
3633     Double_t * ex = g->GetEX();
3634     Double_t * ey = g->GetEY();
3635     for (Int_t i = 0 ; i < g->GetN(); i++) {
3636       mergedGraph->SetPoint(n0+i, x[i], y[i]);
3637       Double_t exPoint = ex ? ex[i] : 0;
3638       Double_t eyPoint = ey ? ey[i] : 0;
3639       mergedGraph->SetPointError(n0+i, exPoint, eyPoint);
3640     }
3641   }
3642   return mergedGraph->GetN();
3643 }
3644
3645 //________________________________________________________________________
3646 Bool_t AliTPCcalibResidualPID::TPCCutMIGeo(const AliVTrack* track, const AliVEvent* evt, TTreeStream* streamer)
3647 {
3648   //
3649   // TPC Cut MIGeo
3650   //
3651
3652   if (!track || !evt)
3653     return kFALSE;
3654
3655   const Short_t sign = track->Charge();
3656   Double_t xyz[50];
3657   Double_t pxpypz[50];
3658   Double_t cv[100];
3659
3660   track->GetXYZ(xyz);
3661   track->GetPxPyPz(pxpypz);
3662
3663   AliExternalTrackParam* par = new AliExternalTrackParam(xyz, pxpypz, cv, sign);
3664   const AliESDtrack dummy;
3665
3666   const Double_t magField = evt->GetMagneticField();
3667   Double_t varGeom = dummy.GetLengthInActiveZone(par, 3, 236, magField, 0, 0);
3668   Double_t varNcr  = track->GetTPCClusterInfo(3, 1);
3669   Double_t varNcls = track->GetTPCsignalN();
3670
3671   const Double_t varEval = 130. - 5. * TMath::Abs(1. / track->Pt());
3672   Bool_t cutGeom   = varGeom > fgCutGeo * varEval;
3673   Bool_t cutNcr    = varNcr  > fgCutNcr * varEval;
3674   Bool_t cutNcls   = varNcls > fgCutNcl * varEval;
3675
3676   Bool_t kout = cutGeom && cutNcr && cutNcls;
3677
3678   if (streamer) {
3679     Double_t dedx = track->GetTPCsignal();
3680
3681     (*streamer)<<"tree"<<
3682       "param.="<< par<<
3683       "varGeom="<<varGeom<<
3684       "varNcr="<<varNcr<<
3685       "varNcls="<<varNcls<<
3686
3687       "cutGeom="<<cutGeom<<
3688       "cutNcr="<<cutNcr<<
3689       "cutNcls="<<cutNcls<<
3690
3691       "kout="<<kout<<
3692       "dedx="<<dedx<<
3693
3694       "\n";
3695   }
3696
3697   delete par;
3698
3699   return kout;
3700 }