Updates for PWGLF-QA
authorlramona <ramona.lea@cern.ch>
Wed, 5 Feb 2014 16:22:15 +0000 (17:22 +0100)
committerlramona <ramona.lea@cern.ch>
Wed, 5 Feb 2014 16:32:47 +0000 (17:32 +0100)
PWGLF/QATasks/AliAnalysisTaskQAHighPtDeDx.cxx
PWGLF/QATasks/macros/AddTaskQAMultistrange.C [new file with mode: 0644]
PWGLF/QATasks/macros/AddTaskQAPhi.C [new file with mode: 0644]
PWGLF/QATasks/post/PostProcessQAHighPtDeDx.C [new file with mode: 0644]
PWGLF/QATasks/post/PostProcessQAMultistrange.C [new file with mode: 0644]
PWGLF/QATasks/post/PostProcessQAPhi.C [new file with mode: 0644]

index 270a696..bcb568e 100644 (file)
@@ -136,6 +136,32 @@ AliAnalysisTaskQAHighPtDeDx::AliAnalysisTaskQAHighPtDeDx():
 
 {
   //default constructor
+  for(Int_t i=0;i<9;++i){
+    
+    hMIPVsNch[i]=0;//TH2D, MIP vs Nch for different eta intervals
+    pMIPVsNch[i]=0;//TProfile, MIP vs Nch for different eta intervals
+    hMIPVsPhi[i]=0;//TH2D, MIP vs phi for different eta intervals
+    pMIPVsPhi[i]=0;//TProfile, MIP vs phi for different eta intervals
+    hPlateauVsPhi[i]=0;//TH2D, dE/dx vs Phi, electrons 0.4<p<0.6 GeV/c
+    pPlateauVsPhi[i]=0;//TProfile, dE/dx vs Phi, electrons 0.4<p<0.6 GeV/c
+    histPiV0[i]=0;//TH2D, dE/dx vs p, pi id by V0s
+    histpPiV0[i]=0;//TH1D, pi id by V0s
+    histPV0[i]=0;// TH2D, dE/dx vs p, p id by V0s
+    histpPV0[i]=0;// TH1D, p id by V0s
+    histAllCh[i]=0;//TH2D, dE/dx vs p for all charged particles
+    histPiTof[i]=0;//TH2D, dE/dx vs p for a "clean" sample of pions, beta>1
+    histpPiTof[i]=0;//TH1D, for a "clean" sample of pions, beta>1
+    histEV0[i]=0;
+    
+    for(Int_t pid=0;pid<7;++pid){
+      hMcIn[pid][i]=0;
+      hMcOut[pid][i]=0;
+    }
+    
+  }
+  
+  
+
 }
 
 
diff --git a/PWGLF/QATasks/macros/AddTaskQAMultistrange.C b/PWGLF/QATasks/macros/AddTaskQAMultistrange.C
new file mode 100644 (file)
index 0000000..dcbadbc
--- /dev/null
@@ -0,0 +1,77 @@
+AliAnalysisTaskQAMultistrange *AddTaskQAMultistrange( Int_t    minnTPCcls             = 70,
+                                                      Float_t  centrlowlim            = 0.,
+                                                      Float_t  centruplim             = 90.,
+                                                      TString  centrest               = "V0M",
+                                                      Bool_t   kusecleaning           = kTRUE, 
+                                                      Float_t  vtxlim                 = 10.,
+                                                      TString  collidingSystem        = "PbPb",
+                                                      Bool_t   SDDSelection           = kFALSE,
+                                                      Bool_t   withSDD                = kFALSE,
+                                                      Float_t  minptondaughtertracks  = 1.,
+                                                      Float_t  etacutondaughtertracks = 9999999.) {
+
+   // Creates, configures and attaches to the train a cascades check task.
+   // Get the pointer to the existing analysis manager via the static access method.
+   //==============================================================================
+   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+   if (!mgr) {
+      ::Error("AddTaskQAMultistrange", "No analysis manager to connect to.");
+      return NULL;
+   }   
+
+   // Check the analysis type using the event handlers connected to the analysis manager.
+   //==============================================================================
+   if (!mgr->GetInputEventHandler()) {
+      ::Error("AddTaskQAMultistrange", "This task requires an input event handler");
+      return NULL;
+   }   
+   TString type = mgr->GetInputEventHandler()->GetDataType(); // can be "ESD" or "AOD"
+
+   // Create and configure the task
+   AliAnalysisTaskQAMultistrange *taskcheckcascade = new AliAnalysisTaskQAMultistrange("TaskCheckCascade");
+
+   taskcheckcascade->SetAnalysisType               (type);
+   taskcheckcascade->SetCollidingSystem            (collidingSystem);
+   taskcheckcascade->SetSDDselection               (SDDSelection);
+   taskcheckcascade->SetQualityCutZprimVtxPos      (kTRUE);             // selects vertices in +-10cm
+   taskcheckcascade->SetQualityCutNoTPConlyPrimVtx (kTRUE);             // retains only events with tracking + SPD vertex
+   taskcheckcascade->SetQualityCutTPCrefit         (kTRUE);             // requires TPC refit flag to be true to select a track
+   taskcheckcascade->SetQualityCutnTPCcls          (kTRUE);             // rejects tracks that have less than n clusters in the TPC
+   taskcheckcascade->SetQualityCutMinnTPCcls       (minnTPCcls);        // minimum number of TPC clusters to accept daughter tracks
+   taskcheckcascade->SetQualityCutPileup           (kFALSE);
+   taskcheckcascade->SetwithSDD                    (withSDD);
+   taskcheckcascade->SetCentralityLowLim           (centrlowlim);       // setting centrality selection vriables
+   taskcheckcascade->SetCentralityUpLim            (centruplim);
+   taskcheckcascade->SetCentralityEst              (centrest);
+   taskcheckcascade->SetUseCleaning                (kusecleaning);
+   taskcheckcascade->SetVertexRange                (vtxlim);
+   taskcheckcascade->SetMinptCutOnDaughterTracks   (minptondaughtertracks);  
+   taskcheckcascade->SetEtaCutOnDaughterTracks     (etacutondaughtertracks);
+   taskcheckcascade->SelectCollisionCandidates();
+
+   mgr->AddTask(taskcheckcascade);
+
+   // Create ONLY the output containers for the data produced by the task.
+   // Get and connect other common input/output containers via the manager as below
+   //==============================================================================
+
+   // User file name (if need be)
+   
+   TString outputFileName = AliAnalysisManager::GetCommonFileName();
+   
+   outputFileName += ":PWGLFStrangeness.outputCheckCascade";
+   
+   Printf("AddTaskCheckCascade - Set OutputFileName : \n %s\n", outputFileName.Data() );
+
+   AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("cfcontCuts",
+                                                             AliCFContainer::Class(),
+                                                             AliAnalysisManager::kOutputContainer,
+                                                             outputFileName );
+
+   
+   mgr->ConnectInput( taskcheckcascade, 0, mgr->GetCommonInputContainer());
+   mgr->ConnectOutput(taskcheckcascade, 1, coutput1);
+   
+   return taskcheckcascade;
+}   
+
diff --git a/PWGLF/QATasks/macros/AddTaskQAPhi.C b/PWGLF/QATasks/macros/AddTaskQAPhi.C
new file mode 100644 (file)
index 0000000..790d210
--- /dev/null
@@ -0,0 +1,5 @@
+AliAnalysisTask *AddTaskQAPhi(const char *dataType = "Rsn_pp_ESD_MINI")
+{
+   gROOT->LoadMacro(gSystem->ExpandPathName("$ALICE_ROOT/PWGLF/RESONANCES/macros/lego_train/AddRsnTaskTrain.C"));
+   return AddRsnTaskTrain(dataType,"Phi","PhiNsigma:KTPCnsig30","","RsnTrainSettingsExtra.C","10.0,-1,0,-1,1,0,10,1,-1,-1,0,1");
+}
diff --git a/PWGLF/QATasks/post/PostProcessQAHighPtDeDx.C b/PWGLF/QATasks/post/PostProcessQAHighPtDeDx.C
new file mode 100644 (file)
index 0000000..f5171b4
--- /dev/null
@@ -0,0 +1,3524 @@
+/*
+
+  Use AliRoot because of AliXRDPROOFtoolkit:
+
+
+
+  gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/Base")
+  gSystem->AddIncludePath("-I$ALICE_ROOT/PWGLF/SPECTRA/IdentifiedHighPt/macros")
+  gSystem->AddIncludePath("-I$ALICE_ROOT/PWGLF/SPECTRA/IdentifiedHighPt/grid")
+  gSystem->AddIncludePath("-I$ALICE_ROOT/PWGLF/SPECTRA/IdentifiedHighPt/lib")
+  gROOT->SetMacroPath(".:$ALICE_ROOT/PWGLF/SPECTRA/IdentifiedHighPt/macros:$ALICE_ROOT/PWGLF/SPECTRA/IdentifiedHighPt/grid:$ALICE_ROOT/PWGLF/SPECTRA/IdentifiedHighPt/lib/")
+  .L my_functions.C+
+  .L my_tools.C+
+  .L PostProcessQAHighPtDeDx.C+
+
+   PlotQA("FileRoot/AnalysisResults.root")
+   MakeFitsExternalData("FileRoot/AnalysisResults.root", "HistosForBB")
+
+  MakeFitsV0s("FileRoot/AnalysisResults.root", "HistosForBB/PrimaryElectrons.root", "HistosForBB",0)
+  MakeFitsV0s("FileRoot/AnalysisResults.root", "HistosForBB/PrimaryElectrons.root", "HistosForBB",1)
+  MakeFitsV0s("FileRoot/AnalysisResults.root", "HistosForBB/PrimaryElectrons.root", "HistosForBB",2)
+  MakeFitsV0s("FileRoot/AnalysisResults.root", "HistosForBB/PrimaryElectrons.root", "HistosForBB",3)
+
+
+  PlotParametrizations("HistosForBB")
+
+
+
+
+  FitDeDxVsP("FileRoot/AnalysisResults.root", 3.0, 10.0, 0, 6, 13, kTRUE,  0,  2, 0,1, 0, "HistosForBB/hres_0_5_02.root",27);
+  FitDeDxVsP("FileRoot/AnalysisResults.root", 3.0, 10.0, 0, 6, 13, kTRUE,  2,  4, 0,1, 0, "HistosForBB/hres_0_5_24.root",27);
+  FitDeDxVsP("FileRoot/AnalysisResults.root", 3.0, 10.0, 0, 6, 13, kTRUE,  4,  6, 0,1, 0, "HistosForBB/hres_0_5_46.root",27);
+  FitDeDxVsP("FileRoot/AnalysisResults.root", 3.0, 10.0, 0, 6, 13, kTRUE,  6,  8, 0,1, 0, "HistosForBB/hres_0_5_68.root",27);
+
+
+  MakeNSigmaPlot("FileRoot/AnalysisResults.root","fitparameters/MB/02_dataPbPb.root",2,50, 0,  "02_dataPbPb.root");
+  MakeNSigmaPlot("FileRoot/AnalysisResults.root","fitparameters/MB/24_dataPbPb.root",2,50, 1,  "24_dataPbPb.root");
+  MakeNSigmaPlot("FileRoot/AnalysisResults.root","fitparameters/MB/46_dataPbPb.root",2,50, 2,  "46_dataPbPb.root");
+  MakeNSigmaPlot("FileRoot/AnalysisResults.root","fitparameters/MB/68_dataPbPb.root",2,50, 3,  "68_dataPbPb.root");
+
+
+  PlotNSigma("nsigma_results")
+
+
+
+  //add particle fractions vs p
+
+
+  ExtractUncPartFractvsP("FileRoot/AnalysisResults.root", 2, 50,0,0, "results/eta02","fractions");
+  ExtractUncPartFractvsP("FileRoot/AnalysisResults.root", 2, 50,0,1, "results/eta24","fractions");
+  ExtractUncPartFractvsP("FileRoot/AnalysisResults.root", 2, 50,0,2, "results/eta46","fractions");
+  ExtractUncPartFractvsP("FileRoot/AnalysisResults.root", 2, 50,0,3, "results/eta68","fractions");
+
+
+*/
+#include <TLegend.h>
+#include <TFile.h>
+#include <TList.h>
+#include <TTree.h>
+#include <TH2.h>
+#include <TF1.h>
+#include <TF2.h>
+#include <TH1.h>
+#include <TCanvas.h>
+#include <TProfile.h>
+#include <TClonesArray.h>
+#include <TObjString.h>
+#include <TString.h>
+#include <TROOT.h>
+#include <TStyle.h>
+#include <TLatex.h>
+#include <TGraphErrors.h>
+#include <TSpectrum.h>
+#include <TFitResultPtr.h>
+#include <TFitResult.h>
+
+#include "my_tools.C"
+#include "my_functions.C"
+
+#include <iostream>
+#include <fstream>
+#include <string>
+
+using namespace std;
+
+
+
+
+
+//const Char_t* ending[4] = {"02", "24", "46", "68"};
+
+const Char_t * legEtaCut[4]=
+  {"|#eta_{lab}|<0.2",  "0.2<|#eta_{lab}|<0.4",   "0.4<|#eta_{lab}|<0.6",   "0.6<|#eta_{lab}|<0.8"};
+
+
+const Char_t* endingCent[1] = 
+  {"MB"};
+const Char_t* CentLatex[1] =
+  {"MB"};
+
+  const Char_t *Pid[4]={"Pion", "Kaon", "Proton", "Electron"};
+//const Char_t *Pid[3]={"Pion","Kaon","Proton"};
+const Char_t *PidLatex[3]={"#pi^{+} + #pi^{-}","K^{+} + K^{-}","p + #bar{p}"};
+
+
+const Color_t PidColor[3]={2,kGreen,4};
+
+TF1* piFunc = 0;
+TF1* kFunc  = 0;
+TF1* pFunc = 0;
+TF1* eFunc = 0;
+TF1* sigmaFunc = 0;
+
+const Int_t nPtBins = 68;
+
+Double_t xBins[nPtBins+1] = {0. ,  0.05, 0.1,  0.15, 0.2,  0.25, 0.3,  0.35, 0.4,  0.45,
+                            0.5,  0.55, 0.6,  0.65, 0.7,  0.75, 0.8,  0.85, 0.9,  0.95,
+                            1.0,  1.1 , 1.2,  1.3 , 1.4,  1.5 , 1.6,  1.7 , 1.8,  1.9 ,
+                            2.0,  2.2 , 2.4,  2.6 , 2.8,  3.0 , 3.2,  3.4 , 3.6,  3.8 ,
+                            4.0,  4.5 , 5.0,  5.5 , 6.0,  6.5 , 7.0,  8.0 , 9.0,  10.0,
+                            11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 18.0, 20.0, 22.0, 24.0,
+                            26.0, 28.0, 30.0, 32.0, 34.0, 36.0, 40.0, 45.0, 50.0 };
+
+
+
+
+const Int_t nPtBinsV0s = 25;
+Double_t ptBinsV0s[nPtBinsV0s+1] = { 0.0 , 0.1 , 0.2 , 0.3 , 0.4 , 0.5 , 0.6 , 0.7 , 0.8 , 0.9 , 1.0 ,
+                                    1.2 , 1.4 , 1.6 , 1.8 , 2.0 , 2.5 , 3.0 , 3.5 , 4.0 , 5.0 , 7.0 ,
+                                    9.0 , 12.0, 15.0, 20.0 };
+
+
+TH2D * AddTwoSameBinningTH2D(TH2D *hPos, TH2D *hNeg, const Char_t *nameHist);
+
+
+void PlotNSigma(const Char_t *path){
+
+  gStyle->SetCanvasColor(10);
+  gStyle->SetFrameFillColor(10);
+  gStyle->SetOptStat(0);
+  gStyle->SetOptTitle(0);
+  gStyle->SetOptFit(0);
+
+
+
+  TFile *fin[2]={0,0};
+  
+  fin[0] = FindFileFresh(Form("%s/02_dataPbPb.root",path));
+  if(!fin[0])
+    return;
+
+  fin[1] = FindFileFresh(Form("%s/68_dataPbPb.root",path));
+  if(!fin[1])
+    return;
+
+  TH1D *hNPiP[2]={0,0}; 
+  TH1D *hNPiK[2]={0,0}; 
+  TH1D *hNPK[2]={0,0}; 
+  Int_t colorestmp[2]={2,4};
+  for(Int_t i=0; i<2; ++i){
+
+    hNPiP[i]=(TH1D *)fin[i]->Get("hPionProtonSeparation");
+    hNPiP[i]->SetLineColor(colorestmp[i]);
+    hNPiP[i]->SetLineWidth(2);
+    hNPiK[i]=(TH1D *)fin[i]->Get("hPionKaonSeparation");
+    hNPiK[i]->SetLineColor(colorestmp[i]);
+    hNPiK[i]->SetLineWidth(2);
+    hNPK[i]=(TH1D *)fin[i]->Get("hKaonProtonSeparation");
+    hNPK[i]->SetLineColor(colorestmp[i]);
+    hNPK[i]->SetLineWidth(2);
+
+  }
+
+
+  TH1D *hframe[3]={0,0,0};
+  for(Int_t i=0; i<3; ++i){
+
+    hframe[i] = new TH1D(Form("hframe%d",i),";#it{p} (GeV/#it{c});",10,3,15);
+    hframe[i]->GetYaxis()->SetRangeUser(-0.01,7.01);
+    hframe[i]->GetXaxis()->SetTitleSize(0.05);
+    hframe[i]->GetYaxis()->SetTitleSize(0.05);
+  }
+
+
+  TCanvas* vC1 = new TCanvas("vC1","vC1",200,10,1050,400);
+  TPad *pad[3]={0,0,0};
+  for(Int_t i=0; i<3; ++i){
+
+    pad[i]=new TPad(Form("pad%d",i),"pad1",0.01+0.33*i,0.01,0.33+0.33*i,0.99,0);
+    pad[i]->SetLeftMargin(0.12);
+    pad[i]->SetRightMargin(0.01);
+    pad[i]->SetTopMargin(0.01);
+    pad[i]->SetBottomMargin(0.1);
+    pad[i]->SetBorderSize(0);
+    pad[i]->SetBorderMode(0);
+
+  }
+
+  
+  
+  for(Int_t i =0 ; i<3; ++i){
+    
+    
+
+
+    vC1->cd();
+    pad[i]->Draw();
+    
+    
+    pad[i]->cd();
+    
+    
+    pad[i]->SetTickx(1);
+    pad[i]->SetTicky(1);
+    
+    if(i==0)
+      hframe[i]->GetYaxis()->SetTitle("Separation #pi/p, N_{#sigma}");
+    if(i==1)
+      hframe[i]->GetYaxis()->SetTitle("Separation #pi/K, N_{#sigma}");
+    if(i==2)
+      hframe[i]->GetYaxis()->SetTitle("Separation p/K, N_{#sigma}");
+
+
+    hframe[i]->Draw();
+
+    if(i==0){
+      for(Int_t j =0 ; j<2; ++j){
+
+       hNPiP[j]->Draw("samel");
+       
+      }
+      TF1 *fpip1=new TF1("fpip1","5.0+pol0",5,15);
+      fpip1->SetLineColor(1);
+      fpip1->SetLineWidth(2);
+      fpip1->SetLineStyle(2);
+      fpip1->Draw("same");
+
+      TF1 *fpip2=new TF1("fpip2","3.5+pol0",5,15);
+      fpip2->SetLineColor(1);
+      fpip2->SetLineWidth(2);
+      fpip2->SetLineStyle(4);
+      fpip2->Draw("same");
+
+      TLatex* latex0 = new TLatex();
+      latex0->SetNDC();
+      latex0->SetTextAlign(22);
+      latex0->SetTextFont(42);
+      latex0->SetTextSize(0.05);
+      latex0->DrawLatex(0.5,0.77,"pp @ 2.76 TeV, 0.6<|#eta|<0.8 (final)");
+
+      latex0->DrawLatex(0.5,0.57,"Pb-Pb, 0-5%, 0.6<|#eta|<0.8 (final)");
+
+    }
+
+    if(i==1){
+      for(Int_t j =0 ; j<2; ++j){
+       hNPiK[j]->Draw("samel");
+       
+      }
+
+
+      TF1 *fpik1=new TF1("fpik1","3.5+pol0",5,15);
+      fpik1->SetLineColor(1);
+      fpik1->SetLineWidth(2);
+      fpik1->SetLineStyle(2);
+      fpik1->Draw("same");
+
+      TF1 *fpik2=new TF1("fpik2","2.2+pol0",5,15);
+      fpik2->SetLineColor(1);
+      fpik2->SetLineWidth(2);
+      fpik2->SetLineStyle(4);
+      fpik2->Draw("same");
+
+      TLegend* legend = new TLegend(0.51, 0.65, 0.85, 0.85);    
+      legend->SetBorderSize(0);
+      legend->SetFillColor(0);
+      legend->SetTextFont(42);
+      legend->SetTextSize(0.05);
+      legend->SetLineColor(kWhite);
+      legend->SetLineStyle(3);
+      legend->SetShadowColor(kWhite);
+      legend->SetFillColor(kWhite);
+      legend->SetFillStyle(0);
+
+      legend->AddEntry(hNPiK[0], "|#eta|<0.2", "L");
+      legend->AddEntry(hNPiK[1], "0.6<|#eta|<0.8", "L");
+      legend->Draw();
+
+
+    }
+
+    if(i==2){
+      for(Int_t j =0 ; j<2; ++j){
+       hNPK[j]->Draw("samel");
+       
+      }
+
+      TF1 *fpk1=new TF1("fpk1","1.8+pol0",5,15);
+      fpk1->SetLineColor(1);
+      fpk1->SetLineWidth(2);
+      fpk1->SetLineStyle(2);
+      fpk1->Draw("same");
+
+      TF1 *fpk2=new TF1("fpk2","1.2+pol0",5,15);
+      fpk2->SetLineColor(1);
+      fpk2->SetLineWidth(2);
+      fpk2->SetLineStyle(4);
+      fpk2->Draw("same");
+
+    }
+
+  }
+
+
+    
+  vC1->SaveAs("NSigma.png");
+  vC1->SaveAs("NSigma.eps");
+    
+
+
+
+
+
+
+}
+
+
+
+//____________________________________________________________________________
+void MakeNSigmaPlot(const Char_t* inFile, const Char_t* fitFileName,
+                   Double_t ptStart, Double_t ptStop, Int_t i_eta,
+                   const Char_t* outFileName=0)
+{
+
+  const Char_t* ending[4] = {"02", "24", "46", "68"};
+  gStyle->SetOptStat(0);
+  
+  if(fitFileName) {
+    
+    TFile* fitFile = FindFileFresh(fitFileName);
+    if(!fitFile)
+      return;
+    DeDxFitInfo* fitPar = (DeDxFitInfo*)fitFile->Get("fitInfo");
+    fitPar->Print();
+    
+    fixMIP      = fitPar->MIP;
+    fixPlateau  = fitPar->plateau;
+
+    Double_t dedxPar[6]  = {0, 0, 0, 0, 0, 0};
+    Double_t sigmaPar[8] = {0, 0, 0, 0, 0, 0, 0, 0};
+    
+    dedxPar[0] = fitPar->optionDeDx;
+    for(Int_t i = 0; i < fitPar->nDeDxPar; i++) {
+      dedxPar[i+1] = fitPar->parDeDx[i];
+    }
+
+    sigmaPar[0] = fitPar->optionSigma;
+    for(Int_t i = 0; i < fitPar->nSigmaPar; i++) {
+      sigmaPar[i+1] = fitPar->parSigma[i];
+    }
+
+    piFunc = new TF1("piFunc", FitFunc, 0, 100, fitPar->nDeDxPar+1);
+    piFunc->SetParameters(dedxPar);
+
+    kFunc = new TF1("kFunc", FitFunc, 0, 100, fitPar->nDeDxPar+1);
+    kFunc->SetParameters(dedxPar);
+    kFunc->SetParameter(0, kFunc->GetParameter(0)+10);
+
+    pFunc = new TF1("pFunc", FitFunc, 0, 100, fitPar->nDeDxPar+1);
+    pFunc->SetParameters(dedxPar);
+    pFunc->SetParameter(0, pFunc->GetParameter(0)+20);
+
+    eFunc = new TF1("eFunc", FitFunc, 0, 100, fitPar->nDeDxPar+1);
+    eFunc->SetParameters(dedxPar);
+    eFunc->SetParameter(0, eFunc->GetParameter(0)+30);
+
+    sigmaFunc = new TF1("sigmaFunc", SigmaFunc, 0, 100, fitPar->nSigmaPar+1); 
+    sigmaFunc->SetParameters(sigmaPar);
+  }
+  kFunc->SetLineColor(kGreen+2);
+  pFunc->SetLineColor(4);
+  eFunc->SetLineColor(kMagenta);
+
+  /*
+  TString baseName(gSystem->BaseName(calibFileName));
+  baseName.ReplaceAll(".root", "");
+  
+  TFile* calibFile = FindFileFresh(calibFileName);
+  if(!calibFile)
+    return;
+  AliHighPtDeDxCalib* calib = (AliHighPtDeDxCalib*)GetObject(calibFile, filter, usePhiCut, run, etaAbs, etaLow, etaHigh);
+  if(calib) {    
+    fixMIP      = calib->GetHistDeDx(kTRUE, 0)->GetMean();
+    fixPlateau  = calib->GetHistDeDx(kFALSE, 0)->GetMean();
+    cout << "Plateau: " << fixPlateau << endl;
+  } else {
+    calib = (AliHighPtDeDxCalib*)calibFile->Get("v0phicut");
+    fixMIP = 50.0;
+    fixPlateau  = 83.461;
+  }
+  calib->Print();
+
+  hDeDxVsP = calib->GetHistDeDxVsP(0);
+  */
+
+  TFile* dedxFile = FindFileFresh(inFile);
+  if(!dedxFile)
+    return;
+
+  TList * list = (TList *)dedxFile->Get("outputdedx");
+  TH2D *hDeDxVsPPlus = 0;
+  hDeDxVsPPlus = (TH2D *)list->FindObject(Form("histAllCh%s",ending[i_eta]));
+  hDeDxVsPPlus->Sumw2();
+
+  TH2D *hDeDxVsPMinus = 0;
+  hDeDxVsPMinus = (TH2D *)list->FindObject(Form("histAllCh%s",ending[i_eta]));
+  hDeDxVsPMinus->Sumw2();
+  
+  gSystem->Exec("rm *.root");
+  TFile *fplus = new TFile("hist2Dplus.root","recreate");
+  fplus->cd();
+  hDeDxVsPPlus->SetName(Form("histAllCh%s", ending[i_eta]));
+  hDeDxVsPPlus->Write();
+  fplus->Close();
+  
+  TFile *fminus = new TFile("hist2Dminus.root","recreate");
+  fminus->cd();
+  hDeDxVsPMinus->SetName(Form("histAllCh%s", ending[i_eta]));
+  hDeDxVsPMinus->Write();
+  fminus->Close();
+  
+  gSystem->Exec("hadd hist2Ddedx.root hist2Dplus.root hist2Dminus.root");
+  
+  TFile *ffinal = TFile::Open("hist2Ddedx.root"); 
+  hDeDxVsP = (TH2D *)ffinal->Get(Form("histAllCh%s", ending[i_eta]));
+
+
+
+
+   //in this case, sigma is not the relative resolution
+
+  // Root is a bit stupid with finidng bins so we have to add and subtract a
+  // little to be sure we get the right bin as we typically put edges as
+  // limits
+  const Int_t binStart = hDeDxVsP->GetXaxis()->FindBin(ptStart+0.01);
+  ptStart = hDeDxVsP->GetXaxis()->GetBinLowEdge(binStart);
+  const Int_t binStop  = hDeDxVsP->GetXaxis()->FindBin(ptStop-0.01);
+  ptStop  = hDeDxVsP->GetXaxis()->GetBinUpEdge(binStop);
+  //  const Int_t nBins    = binStop - binStart + 1;
+
+  cout << "Doing fits from pTlow = " << ptStart << " (bin: " << binStart
+       << ") to pThigh = " << ptStop << " (bin: " << binStop << ")" << endl;
+
+
+  TH1D* hPionRatio =(TH1D*)hDeDxVsP->ProjectionX("hPionRatio", 1, 1);
+  hPionRatio->Reset();
+
+  TH1D* hPionKaonSeparation = (TH1D*)hPionRatio->Clone("hPionKaonSeparation");
+  hPionKaonSeparation->SetTitle(Form("#pi/K separation, %s; #it{p} (GeV/#it{c}); N_{#sigma}",legEtaCut[i_eta]));
+  TH1D* hPionProtonSeparation = (TH1D*)hPionRatio->Clone("hPionProtonSeparation");
+  hPionProtonSeparation->SetTitle(Form("#pi/p separation, %s; #it{p} (GeV/#it{c}); N_{#sigma}",legEtaCut[i_eta]));
+  TH1D* hKaonProtonSeparation = (TH1D*)hPionRatio->Clone("hKaonProtonSeparation");
+  hKaonProtonSeparation->SetTitle(Form("K/p separation, %s; #it{p} (GeV/#it{c}); N_{#sigma}",legEtaCut[i_eta]));
+
+
+
+
+  for(Int_t bin = binStart; bin <= binStop; bin++){
+    
+    Double_t pt=hPionRatio->GetBinCenter(bin);
+
+    Double_t dEdx_pi = piFunc->Eval(pt);
+    Double_t Sigma_pi = sigmaFunc->Eval(piFunc->Eval(pt));
+
+    Double_t dEdx_k = kFunc->Eval(pt);
+    Double_t Sigma_k = sigmaFunc->Eval(kFunc->Eval(pt));
+
+    Double_t dEdx_p = pFunc->Eval(pt);
+    Double_t Sigma_p = sigmaFunc->Eval(pFunc->Eval(pt));
+
+    Double_t N1 = (dEdx_pi - dEdx_k)/( (Sigma_pi + Sigma_k)/2.0 );
+    Double_t N2 = (dEdx_pi - dEdx_p)/( (Sigma_pi + Sigma_p)/2.0 );
+    Double_t N3 = (dEdx_k - dEdx_p)/( (Sigma_k + Sigma_p)/2.0 );
+
+    hPionKaonSeparation->SetBinContent(bin,N1);
+    hPionProtonSeparation->SetBinContent(bin,N2);
+    hKaonProtonSeparation->SetBinContent(bin,N3);
+
+  }
+
+
+
+  if(outFileName) {
+
+    CreateDir("nsigma_results");
+    TFile* fileOut = new TFile(Form("nsigma_results/%s", outFileName), "RECREATE");
+
+    hPionKaonSeparation->Write();
+    hPionProtonSeparation->Write();
+    hKaonProtonSeparation->Write();
+
+
+
+    fileOut->Close();
+  }
+
+
+
+}
+
+
+
+//________________________________________________________________________
+void PlotParametrizations(const Char_t *path){
+
+
+  gStyle->SetCanvasColor(10);
+  gStyle->SetFrameFillColor(10);
+  gStyle->SetOptStat(0);
+  gStyle->SetOptTitle(0);
+  gStyle->SetOptFit(0);
+
+  TFile *fin[2]={0,0};
+
+
+  fin[0] = FindFileFresh(Form("%s/hres_0_5_02.root",path));
+  if(!fin[0])
+    return;
+
+  fin[1] = FindFileFresh(Form("%s/hres_0_5_68.root",path));
+  if(!fin[1])
+    return;
+
+
+  TH1D *hframeRes = new TH1D("hframeRes","; #LT dE/dx #GT; #sigma/#LT dE/dx #GT",50,40,79);
+  hframeRes->GetYaxis()->SetRangeUser(0,0.105);
+  hframeRes->GetYaxis()->SetTitleSize(0.05);
+  hframeRes->GetXaxis()->SetTitleSize(0.05);
+
+  TH1D *hframeBB = new TH1D("hframeBB","; #beta#gamma;#LT dE/dx #GT",50,2,59);
+  hframeBB->GetYaxis()->SetRangeUser(40,81);
+  hframeBB->GetYaxis()->SetTitleSize(0.05);
+  hframeBB->GetXaxis()->SetTitleSize(0.05);
+
+  TGraphErrors *gRes[2] = {0,0}; 
+
+  gRes[0] = (TGraphErrors *)fin[0]->Get("res_allpions");
+  gRes[1] = (TGraphErrors *)fin[1]->Get("res_allpions");
+
+
+  TGraphErrors *gBB[2] = {0,0}; 
+
+  gBB[0] = (TGraphErrors *)fin[0]->Get("gBB");
+  gBB[1] = (TGraphErrors *)fin[1]->Get("gBB");
+
+  Int_t colorestmp[2]={2,4};
+
+
+
+  TCanvas* vC1 = new TCanvas("vC1","vC1",200,10,700,400);
+  TPad * pad1=new TPad("pad1","pad1",0.01,0.01,0.5,0.99,0);
+  TPad * pad2=new TPad("pad2","pad2",0.5,0.01,0.99,0.99,0);
+  
+  pad1->SetLeftMargin(0.12);
+  pad1->SetRightMargin(0.01);
+  pad1->SetTopMargin(0.01);
+  pad1->SetBottomMargin(0.1);
+  pad1->SetBorderSize(0);
+  pad1->SetBorderMode(0);
+  
+  pad2->SetLeftMargin(0.12);
+  pad2->SetRightMargin(0.01);
+  pad2->SetBottomMargin(0.1);
+  pad2->SetTopMargin(0);
+  pad2->SetBorderSize(0);
+  vC1->cd();
+  pad1->Draw();
+  
+  
+  pad1->cd();
+
+
+  pad1->SetTickx(1);
+  pad1->SetTicky(1);
+
+
+
+
+
+
+  hframeRes->Draw();
+
+  for(Int_t i =0 ; i<2; ++i){
+    gRes[i]->SetMarkerStyle(20);
+    gRes[i]->SetLineWidth(2);
+    gRes[i]->SetLineColor(colorestmp[i]);
+    gRes[i]->SetMarkerColor(colorestmp[i]);
+    gRes[i]->Draw("samep");
+  }
+
+
+  TF1 *f0005 = new TF1("f0005","pol2",50,80);
+  //  f0005->SetParameter(0,1.20000000000000000e+01);
+  f0005->SetParameter(0,6.08803411659195118e-02);
+  f0005->SetParameter(1,2.34760597152924448e-04);
+  f0005->SetParameter(2,-2.81841683363273653e-06);
+
+  f0005->SetLineColor(1);
+  f0005->SetLineWidth(2);
+  f0005->SetLineStyle(2);
+  f0005->Draw("samel");
+
+
+  TF1 *fpp = new TF1("fpp","pol2",50,80);
+  fpp->SetParameter(0,1.06773399595441187e-01);
+  fpp->SetParameter(1,-1.55087701882609280e-03);
+  fpp->SetParameter(2,1.01579672586848698e-05);
+
+  fpp->SetLineColor(1);
+  fpp->SetLineWidth(2);
+  fpp->SetLineStyle(3);
+  fpp->Draw("samel");
+
+
+  TLegend* legendRes = new TLegend(0.21, 0.15, 0.65, 0.35);    
+  legendRes->SetBorderSize(0);
+  legendRes->SetFillColor(0);
+  legendRes->SetTextFont(42);
+  legendRes->SetTextSize(0.05);
+  legendRes->SetLineColor(kWhite);
+  legendRes->SetLineStyle(3);
+  legendRes->SetShadowColor(kWhite);
+  legendRes->SetFillColor(kWhite);
+  legendRes->SetFillStyle(0);
+  legendRes->SetHeader("Final, 0.6<|#eta|<0.8");
+  legendRes->AddEntry(f0005, "Pb-Pb 0-5%", "L");
+  legendRes->AddEntry(fpp, "pp @ 2.76 TeV", "L");
+  legendRes->Draw();
+
+  vC1->cd();
+  pad2->Draw();
+  
+  
+  pad2->cd();
+
+
+  pad2->SetTickx(1);
+  pad2->SetTicky(1);
+
+
+
+
+
+
+  hframeBB->Draw();
+
+  for(Int_t i =0 ; i<2; ++i){
+    gBB[i]->SetMarkerStyle(20);
+    gBB[i]->SetLineWidth(2);
+    gBB[i]->SetLineColor(colorestmp[i]);
+    gBB[i]->SetMarkerColor(colorestmp[i]);
+    gBB[i]->Draw("samep");
+  }
+
+  TLegend* legend = new TLegend(0.51, 0.8, 0.85, 0.95);    
+  legend->SetBorderSize(0);
+  legend->SetFillColor(0);
+  legend->SetTextFont(42);
+  legend->SetTextSize(0.05);
+  legend->SetLineColor(kWhite);
+  legend->SetLineStyle(3);
+  legend->SetShadowColor(kWhite);
+  legend->SetFillColor(kWhite);
+  legend->SetFillStyle(0);
+  
+  legend->AddEntry(gBB[0], "|#eta|<0.2", "P");
+  legend->AddEntry(gBB[1], "0.6<|#eta|<0.8", "P");
+  legend->Draw();
+
+
+
+
+  vC1->SaveAs("Parametrizations.png");
+  vC1->SaveAs("Parametrizations.eps");
+
+
+
+}
+//____________________________________________________________________________
+void ExtractUncPartFractvsP(const Char_t * inFile,  Double_t ptStart, Double_t ptStop,
+                           Int_t i_cent=1,
+                           Int_t i_eta=1,
+                           const Char_t* dirName="debugfits", const Char_t* outFileName=0)
+
+
+{
+
+
+  const Char_t* ending[4] = {"02", "24", "46", "68"};
+  const Double_t etaHigh[4]={2,4,6,8};
+  gStyle->SetOptStat(0);
+  
+  
+  
+  TFile* fitFile = FindFileFresh(Form("fitparameters/MB/%s_dataPbPb.root",ending[i_eta]));
+  if(!fitFile)
+    return;
+  
+  DeDxFitInfo* fitPar = (DeDxFitInfo*)fitFile->Get("fitInfo");
+  fitPar->Print();
+  
+  fixMIP      = fitPar->MIP;
+  
+
+  //return;
+
+  Double_t dedxPar[6]  = {0, 0, 0, 0, 0, 0};
+  Double_t sigmaPar[8] = {0, 0, 0, 0, 0, 0, 0, 0};
+  
+  dedxPar[0] = fitPar->optionDeDx;
+  for(Int_t i = 0; i < fitPar->nDeDxPar; i++) {
+    dedxPar[i+1] = fitPar->parDeDx[i];
+    cout<<"idedx_par"<<i+1<<"="<<dedxPar[i+1]<<endl;
+
+  }
+  fixPlateau  = dedxPar[5];
+  cout<<"fixPlateau="<<fixPlateau<<"  fixMIP="<<fixMIP<<endl;
+  //return;
+  
+  sigmaPar[0] = fitPar->optionSigma;
+  for(Int_t i = 0; i < fitPar->nSigmaPar; i++) {
+    sigmaPar[i+1] = fitPar->parSigma[i];
+  }
+  
+  piFunc = new TF1("piFunc", FitFunc, 0, 100, fitPar->nDeDxPar+1);
+  piFunc->SetParameters(dedxPar);
+  
+  
+  kFunc = new TF1("kFunc", FitFunc, 0, 100, fitPar->nDeDxPar+1);
+  kFunc->SetParameters(dedxPar);
+  kFunc->SetParameter(0, kFunc->GetParameter(0)+10);
+  
+  pFunc = new TF1("pFunc", FitFunc, 0, 100, fitPar->nDeDxPar+1);
+  pFunc->SetParameters(dedxPar);
+  pFunc->SetParameter(0, pFunc->GetParameter(0)+20);
+  
+  eFunc = new TF1("eFunc", FitFunc, 0, 100, fitPar->nDeDxPar+1);
+  eFunc->SetParameters(dedxPar);
+  eFunc->SetParameter(0, eFunc->GetParameter(0)+30);
+  
+  sigmaFunc = new TF1("sigmaFunc", SigmaFunc, 0, 100, fitPar->nSigmaPar+1); 
+  sigmaFunc->SetParameters(sigmaPar);
+  
+  
+  kFunc->SetLineColor(kGreen+2);
+  pFunc->SetLineColor(4);
+  eFunc->SetLineColor(kMagenta);
+  /*
+    TFile* calibFile = FindFileFresh(inFile);
+    if(!calibFile)
+    return;
+    hDeDxVsP = (TH2D *)calibFile->Get(Form("hDeDxVsP%s_%s", ending[i_eta], endingCent[i_cent]));
+    hDeDxVsP->Sumw2();
+  */
+
+  TFile* dedxFile = FindFileFresh(inFile);
+  if(!dedxFile)
+    return;
+
+  TList * list = (TList *)dedxFile->Get("outputdedx");
+  TH2D *hDeDxVsPPlus = 0;
+  hDeDxVsPPlus = (TH2D *)list->FindObject(Form("histAllCh%s",ending[i_eta]));
+  hDeDxVsPPlus->Sumw2();
+
+  TH2D *hDeDxVsPMinus = 0;
+  hDeDxVsPMinus = (TH2D *)list->FindObject(Form("histAllCh%s",ending[i_eta]));
+  hDeDxVsPMinus->Sumw2();
+  
+  gSystem->Exec("rm *.root");
+  TFile *fplus = new TFile("hist2Dplus.root","recreate");
+  fplus->cd();
+  hDeDxVsPPlus->SetName(Form("histAllCh%s", ending[i_eta]));
+  hDeDxVsPPlus->Write();
+  fplus->Close();
+  
+  TFile *fminus = new TFile("hist2Dminus.root","recreate");
+  fminus->cd();
+  hDeDxVsPMinus->SetName(Form("histAllCh%s", ending[i_eta]));
+  hDeDxVsPMinus->Write();
+  fminus->Close();
+  
+  gSystem->Exec("hadd hist2Ddedx.root hist2Dplus.root hist2Dminus.root");
+  
+  TFile *ffinal = TFile::Open("hist2Ddedx.root"); 
+  hDeDxVsP = (TH2D *)ffinal->Get(Form("histAllCh%s", ending[i_eta]));
+  
+
+
+
+
+
+
+  CreateDir(Form("fit_results/%s_%s",dirName,endingCent[i_cent]));
+  
+  TCanvas* cDeDxVsPLogX = new TCanvas("cDeDxVsPLogX", "; #it{p} (GeV/#it{c}) ;d#it{E}/d#it{x}", 500, 400);
+  cDeDxVsPLogX->Clear();
+  cDeDxVsPLogX->cd();
+  cDeDxVsPLogX->SetLogz();
+  
+  TH2D *hDeDxVsP2=(TH2D *)hDeDxVsP->Clone("h2DP");
+  hDeDxVsP2->SetTitle("; #it{p} (GeV/#it{c}) ;d#it{E}/d#it{x}");
+  hDeDxVsP2->GetXaxis()->SetRangeUser(2,19.5);
+  hDeDxVsP2->Draw("COLZ");
+  
+  piFunc->SetLineColor(1);
+  piFunc->Draw("samel");
+  
+  eFunc->SetLineColor(1);
+  eFunc->Draw("samel");
+  
+  kFunc->SetLineColor(1);
+  kFunc->Draw("samel");
+  
+  pFunc->SetLineColor(1);
+  pFunc->Draw("samel");
+  
+  TLatex* latex0 = new TLatex();
+  latex0->SetNDC();
+  latex0->SetTextAlign(22);
+  latex0->SetTextSize(0.05);
+  latex0->DrawLatex(0.7,0.94+0.035,Form("p-Pb, %s",legEtaCut[i_eta]));
+  latex0->DrawLatex(0.71,0.89+0.035,Form("#sqrt{s_{NN}}=5.02 TeV, %s",CentLatex[i_cent]));
+
+  
+  
+  cDeDxVsPLogX->SaveAs(Form("fit_results/%s_%s/%s.png",dirName,endingCent[i_cent],cDeDxVsPLogX->GetName()));
+  cDeDxVsPLogX->SaveAs(Form("fit_results/%s_%s/%s.eps",dirName,endingCent[i_cent],cDeDxVsPLogX->GetName()));
+  //in this case, sigma is not the relative resolution
+  
+  // Root is a bit stupid with finidng bins so we have to add and subtract a
+  // little to be sure we get the right bin as we typically put edges as
+  // limits
+  const Int_t binStart = hDeDxVsP->GetXaxis()->FindBin(ptStart+0.01);
+  ptStart = hDeDxVsP->GetXaxis()->GetBinLowEdge(binStart);
+  const Int_t binStop  = hDeDxVsP->GetXaxis()->FindBin(ptStop-0.01);
+  ptStop  = hDeDxVsP->GetXaxis()->GetBinUpEdge(binStop);
+  //  const Int_t nBins    = binStop - binStart + 1;
+  
+  cout << "Doing fits from pTlow = " << ptStart << " (bin: " << binStart
+       << ") to pThigh = " << ptStop << " (bin: " << binStop << ")" << endl;
+  
+  
+  //cross check
+  TCanvas* cFits = new TCanvas("cFits", "Fit comparison to data", 1200, 800);
+  cFits->Clear();
+  cFits->Divide(7, 4);
+  
+  
+  TF1* pion = new TF1("pion", "gausn", 40, 100);
+  
+  TF1* kaon = new TF1("kaon", "gausn", 40, 100);
+  
+  TF1* proton = new TF1("proton", "gausn", 40, 100);
+  //proton->SetLineWidth(2);
+  //proton->SetLineColor(kBlue);
+  TF1* electron = new TF1("electron", "gausn", 40, 100);
+  //electron->SetLineWidth(2);
+  //electron->SetLineColor(kMagenta);
+  //  TF1* total = new TF1("total", "gausn(0)+gausn(3)+gausn(6)+gausn(9)", -30, 30);
+  TF1* total = 0;
+  total = new TF1("total", "([0]+[12])*gausn(1)+[4]*gausn(5)+[8]*gausn(9)+[12]*gausn(13)", 40, 100);
+  
+  total->SetLineColor(kBlack);
+  total->SetLineWidth(2);
+  total->SetLineStyle(2);
+  
+  TLegend* legend = new TLegend(0.11, 0.6, 0.35, 0.85);    
+  legend->SetBorderSize(0);
+  legend->SetFillColor(0);
+  legend->AddEntry(total, "4-Gauss fit", "L");
+  legend->AddEntry(pion, "#pi", "L");
+  legend->AddEntry(kaon, "K", "L");
+  legend->AddEntry(proton, "p", "L");
+  legend->AddEntry(electron, "e", "L");
+  
+  TCanvas* cSingleFit = new TCanvas("cSingleFit", "single fit");
+  cSingleFit->SetLeftMargin(0.11);
+  cSingleFit->SetRightMargin(0.03);
+  cSingleFit->SetTopMargin(0.01);
+  cSingleFit->SetBottomMargin(0.1);
+  cSingleFit->SetBorderSize(0);
+  cSingleFit->SetBorderMode(0);
+  
+  
+  TH1D* hPionRatio =(TH1D*)hDeDxVsP->ProjectionX("hPionRatio", 1, 1);
+  hPionRatio->SetTitle("particle fractions; p [GeV/c]; particle fraction");
+  hPionRatio->Reset();
+  TH1D* hKaonRatio   = (TH1D*)hPionRatio->Clone("hKaonRatio");
+  TH1D* hProtonRatio = (TH1D*)hPionRatio->Clone("hProtonRatio");
+  TH1D* hElectronRatio = (TH1D*)hPionRatio->Clone("hElectronRatio");
+  
+    
+    
+  TH1D* hPionYield =(TH1D*)hDeDxVsP->ProjectionX("hPionYield", 1, 1);
+  hPionYield->SetTitle("particle fractions; p [GeV/c]; particle fraction");
+  hPionYield->Reset();
+  TH1D* hKaonYield   = (TH1D*)hPionYield->Clone("hKaonYield");
+  TH1D* hProtonYield = (TH1D*)hPionYield->Clone("hProtonYield");
+  TH1D* hElectronYield = (TH1D*)hPionYield->Clone("hElectronYield");
+  
+  
+  TF1 *fmip=new TF1("fmip","pol0",0,1);
+  fmip->SetParameter(0,fixMIP);
+  
+  
+  
+  TF1* fElectronFraction = 0;
+  if(etaHigh[i_eta]==8||etaHigh[i_eta]==-6)
+    fElectronFraction = new TF1("fElectronFraction", "[0]+(x<9.0)*[1]*(x-9.0)", 0, ptStop);
+  if(etaHigh[i_eta]==6||etaHigh[i_eta]==-4)
+    fElectronFraction = new TF1("fElectronFraction", "[0]+(x<8.0)*[1]*(x-8.0)", 0, ptStop);
+  if(etaHigh[i_eta]==4||etaHigh[i_eta]==-2)
+    fElectronFraction = new TF1("fElectronFraction", "[0]+(x<7.5)*[1]*(x-7.5)", 0, ptStop);
+  if(etaHigh[i_eta]==2||etaHigh[i_eta]==0)
+    fElectronFraction = new TF1("fElectronFraction", "[0]+(x<7.0)*[1]*(x-7.0)", 0, ptStop);
+  
+  //TF1* fElectronFraction = new TF1("fElectronFraction", "[0]+(x<7.0)*[1]*(x-7.0)", 0, ptStop);
+  fElectronFraction->SetParameters(0.01, 0.0);
+  
+  
+  TH1D *dEdxpoints[binStop-binStart];
+  TF1 *fPion[binStop-binStart];
+  TF1 *fKaon[binStop-binStart];
+  TF1 *fProton[binStop-binStart];
+  TF1 *fElectron[binStop-binStart];
+  TF1 *fTotal[binStop-binStart];
+  
+  
+  for(Int_t bin = binStart; bin <= binStop; bin++){
+      
+    Double_t pt=hPionRatio->GetBinCenter(bin);
+    
+    cout << "Making projection for bin: " << bin << endl;
+    
+    const Int_t j = bin-binStart;
+    
+    TH1D* hDeltaPiVsPtProj =(TH1D*)hDeDxVsP->ProjectionY(Form("hDeltaPiVsPtProj%d", bin), bin, bin);
+    //    hDeltaPiVsPtProj->GetXaxis()->SetRangeUser(-25, 20);
+    hDeltaPiVsPtProj->GetXaxis()->SetRangeUser(40, 90);
+    hDeltaPiVsPtProj->SetTitle(Form("%.1f<p<%.1f GeV/c", 
+                                   hDeDxVsP->GetXaxis()->GetBinLowEdge(bin),
+                                   hDeDxVsP->GetXaxis()->GetBinUpEdge(bin)));
+    
+    dEdxpoints[j]=0;
+    dEdxpoints[j]=hDeltaPiVsPtProj;
+    fPion[j]=0;
+    fPion[j]=new TF1(Form("fPiongauss_%d",bin), "gausn", 40, 100);
+    
+    fKaon[j]=0;
+    fKaon[j]=new TF1(Form("fKaongauss_%d",bin), "gausn", 40, 100);
+    
+    fProton[j]=0;
+    fProton[j]=new TF1(Form("fProtongauss_%d",bin), "gausn", 40, 100);
+    
+    fElectron[j]=0;
+    fElectron[j]=new TF1(Form("fElectrongauss_%d",bin), "gausn", 40, 100);
+    
+    fTotal[j]=0;
+    fTotal[j]=new TF1(Form("fTotalgauss_%d",bin), "([0]+[12])*gausn(1)+[4]*gausn(5)+[8]*gausn(9)+[12]*gausn(13)", 40, 100);
+    
+    
+    Double_t all =  hDeltaPiVsPtProj->GetEntries();
+    
+    
+    const Int_t nPar = 16;
+    Double_t gausParams[nPar] = { 
+      0.59,
+      all,
+      piFunc->Eval(pt), 
+      //fpi->Eval(pt),
+      sigmaFunc->Eval(piFunc->Eval(pt)),
+      
+      0.3,
+      all,
+      kFunc->Eval(pt), 
+      //fp->Eval(pt),
+      sigmaFunc->Eval(kFunc->Eval(pt)),
+      
+      0.1,
+      all,
+      pFunc->Eval(pt), 
+      //fp->Eval(pt),
+      sigmaFunc->Eval(pFunc->Eval(pt)),
+      
+      0.01,
+      all,
+      eFunc->Eval(pt),
+      //fpi->Eval(pt),
+      sigmaFunc->Eval(eFunc->Eval(pt)),
+    };
+    
+   
+    for(Int_t ipar=0;ipar<nPar;++ipar)
+      cout<<gausParams[ipar]<<endl;
+    
+    
+    
+    cFits->cd();
+    cFits->cd(j + 1);
+    
+    total->SetParameters(gausParams);
+    
+    
+    for(Int_t i = 0; i < nPar; i++) {
+      if(i!=0 && i!=4 && i!=8 && i!=12)
+       total->FixParameter(i, gausParams[i]);
+      else
+       continue;
+    }
+    
+    total->SetParLimits(7, gausParams[7]-0.05*gausParams[7],gausParams[7]+0.05*gausParams[7]);
+    total->SetParLimits(8, 0.005,0.4);
+    total->SetParLimits(4, 0.005,0.6);
+    
+    if(bin==48) {
+      if(etaHigh[i_eta]==8||etaHigh[i_eta]==-6)
+       hElectronRatio->Fit(fElectronFraction, "N", "", 4.0, 10.0);
+      if(etaHigh[i_eta]==6||etaHigh[i_eta]==-4)
+       hElectronRatio->Fit(fElectronFraction, "N", "", 3.66, 10.0);
+      if(etaHigh[i_eta]==4||etaHigh[i_eta]==-2)
+       hElectronRatio->Fit(fElectronFraction, "N", "", 3.33, 10.0);
+      if(etaHigh[i_eta]==2||etaHigh[i_eta]==0)
+       hElectronRatio->Fit(fElectronFraction, "N", "", 3.0, 10.0);
+      fElectronFraction->SetRange(0, ptStop);
+    }
+    
+    if(bin>=48) {
+      total->FixParameter(12, fElectronFraction->Eval(hDeDxVsP->GetXaxis()->GetBinCenter(bin)));
+    }
+     
+      
+    hDeltaPiVsPtProj->Fit(total, "0L");
+    
+    
+    hDeltaPiVsPtProj->Draw();
+    total->DrawCopy("same");    
+    
+    Double_t parametersOut[nPar];
+    total->GetParameters(parametersOut);
+    const Double_t* parameterErrorsOut = total->GetParErrors();
+    
+    for(Int_t i = 0; i < nPar; i++) 
+      cout << parametersOut[i] << ", ";
+    cout << endl;
+    fTotal[j]->SetParameters(parametersOut);
+    
+    
+
+    pion->SetParameters(&parametersOut[1]);
+    pion->SetParameter(0,all*(parametersOut[0]+parametersOut[12]));
+    
+    fPion[j]->SetParameters(&parametersOut[1]);;
+    fPion[j]->SetParameter(0,all*(parametersOut[0]+parametersOut[12]));
+    
+    
+    hPionRatio->SetBinContent(bin, parametersOut[0]);
+    hPionRatio->SetBinError(bin, parameterErrorsOut[0]);
+    hPionYield->SetBinContent(bin, parametersOut[0]*all);
+    hPionYield->SetBinError(bin, parameterErrorsOut[0]*all);
+      
+    kaon->SetParameters(&parametersOut[5]);
+    kaon->SetParameter(0,all*parametersOut[4]);
+    
+    fKaon[j]->SetParameters(&parametersOut[5]);
+    fKaon[j]->SetParameter(0,all*parametersOut[4]);
+    
+    
+    hKaonRatio->SetBinContent(bin, parametersOut[4]);
+    hKaonRatio->SetBinError(bin, parameterErrorsOut[4]);
+    hKaonYield->SetBinContent(bin, parametersOut[4]*all);
+    hKaonYield->SetBinError(bin, parameterErrorsOut[4]*all);
+    
+    proton->SetParameters(&parametersOut[9]);
+    proton->SetParameter(0,all*parametersOut[8]);
+    //proton->DrawCopy("same");
+    
+    fProton[j]->SetParameters(&parametersOut[9]);
+    fProton[j]->SetParameter(0,all*parametersOut[8]);
+    
+    
+    hProtonRatio->SetBinContent(bin, parametersOut[8]);
+    hProtonRatio->SetBinError(bin, parameterErrorsOut[8]);
+    hProtonYield->SetBinContent(bin, parametersOut[8]*all);
+    hProtonYield->SetBinError(bin, parameterErrorsOut[8]*all);
+    
+    electron->SetParameters(&parametersOut[13]);
+    electron->SetParameter(0,all*parametersOut[12]);
+    
+    fElectron[j]->SetParameters(&parametersOut[13]);
+    fElectron[j]->SetParameter(0,all*parametersOut[12]);
+      
+      
+    //electron->DrawCopy("same");
+    hElectronRatio->SetBinContent(bin, parametersOut[12]);
+    hElectronRatio->SetBinError(bin, parameterErrorsOut[12]);
+    hElectronYield->SetBinContent(bin, parametersOut[12]*all);
+    hElectronYield->SetBinError(bin, parameterErrorsOut[12]*all);
+    
+    //DrawALICELogo(kFALSE, 0.71, 0.76, 0.82, 0.98);
+    cSingleFit->cd();
+    cSingleFit->Clear();
+    cSingleFit->SetTickx(1);
+    cSingleFit->SetTicky(1);
+    cSingleFit->SetLogy(0);
+    
+    //    cSingleFit->SetLogy();
+    hDeltaPiVsPtProj->GetXaxis()->SetRangeUser(40,90);
+    hDeltaPiVsPtProj->SetMarkerStyle(25);
+    hDeltaPiVsPtProj->SetMarkerColor(1);
+    hDeltaPiVsPtProj->SetLineColor(1);
+    hDeltaPiVsPtProj->SetLineWidth(2);
+    hDeltaPiVsPtProj->GetYaxis()->SetLabelFont(42);
+    hDeltaPiVsPtProj->GetXaxis()->SetLabelFont(42);
+    hDeltaPiVsPtProj->GetZaxis()->SetLabelFont(42);
+    hDeltaPiVsPtProj->GetYaxis()->SetTitleFont(42);
+    hDeltaPiVsPtProj->GetXaxis()->SetTitleFont(42);
+    hDeltaPiVsPtProj->GetZaxis()->SetTitleFont(42);
+    hDeltaPiVsPtProj->GetYaxis()->SetTitle("Entries");
+    hDeltaPiVsPtProj->GetXaxis()->SetTitle("d#it{E}/d#it{x} in TPC (arb. units)");
+    hDeltaPiVsPtProj->GetXaxis()->SetTitleSize(0.05);
+    hDeltaPiVsPtProj->GetXaxis()->SetTitleOffset(0.9);
+    hDeltaPiVsPtProj->GetYaxis()->SetTitleSize(0.05);
+    hDeltaPiVsPtProj->GetYaxis()->SetTitleOffset(0.9);
+    hDeltaPiVsPtProj->Draw();
+    
+    total->SetLineColor(kGray+1);
+    total->SetLineWidth(3);
+    total->SetLineStyle(1);
+    total->DrawCopy("same");
+    
+    pion->GetXaxis()->SetRangeUser(40,90);
+    pion->SetLineColor(2);
+    pion->SetLineWidth(2);
+    //pion->SetLineStyle(2);
+    pion->DrawCopy("same");
+    
+    kaon->SetLineColor(kGreen+2);
+    kaon->SetLineWidth(2);
+    //kaon->SetLineStyle(2);
+    kaon->DrawCopy("same");
+    
+    
+    proton->SetLineColor(4);
+    proton->SetLineWidth(2);
+    //proton->SetLineStyle(2);
+    proton->DrawCopy("same");
+    
+    
+    TLatex* latex = new TLatex();
+    //  latex->SetTextFont(132);
+    latex->SetNDC();
+    latex->SetTextAlign(22);
+    
+    latex->SetTextSize(0.05);
+    latex->SetTextFont(42);
+    latex->SetTextColor(kRed+2);
+    //latex->DrawLatex(0.8,0.3,"pp");
+    //latex->DrawLatex(0.81,0.25,"#sqrt{s}=2.76 TeV");
+    latex->DrawLatex(0.8,0.25+0.035,Form("p-Pb, %s",CentLatex[i_cent]));
+    latex->DrawLatex(0.81,0.2+0.035,"#sqrt{s_{NN}}=5.02 TeV");
+    
+    latex->SetTextColor(1);
+    latex->SetTextSize(0.04);
+    //latex->DrawLatex(0.81,0.2,"25/07/2012");
+    //latex->DrawLatex(0.80,0.3+0.036,"31/10/2012");
+    
+    latex->SetTextSize(0.05);
+    latex->DrawLatex(0.8,0.71,legEtaCut[i_eta]);
+    latex->DrawLatex(0.8,0.64,Form("%.1f<#it{p}<%.1f GeV/#it{c}", hDeDxVsP->GetXaxis()->GetBinLowEdge(bin), hDeDxVsP->GetXaxis()->GetBinUpEdge(bin)));
+    
+    
+    
+    //cSingleFit->SaveAs(Form("%s/ptspectrum_bin%d.gif", dirName, bin));
+    cSingleFit->SaveAs(Form("fit_results/%s_%s/ptspectrum_bin%d.png", dirName, endingCent[i_cent],bin));
+    cSingleFit->SaveAs(Form("fit_results/%s_%s/ptspectrum_bin%d.pdf", dirName, endingCent[i_cent],bin));
+    cSingleFit->SaveAs(Form("fit_results/%s_%s/ptspectrum_bin%d.eps", dirName, endingCent[i_cent],bin));
+    //    legend->Draw();
+    
+    
+    
+      
+  }
+
+    
+  TCanvas* cRatio = new TCanvas("cRatio", "ratios/all vs p; pt", 600, 400);
+  
+  cRatio->Clear();
+  cRatio->SetGridy();
+  hElectronRatio->GetYaxis()->SetRangeUser(-0.05, 0.1);
+  hElectronRatio->DrawCopy("P E");
+  fElectronFraction->DrawCopy("SAME");
+  
+  gROOT->ProcessLine(".x drawText.C");
+  cRatio->SaveAs(Form("fit_results/%s_%s/electron_ratio.gif", dirName,endingCent[i_cent]));
+  cRatio->SaveAs(Form("fit_results/%s_%s/electron_ratio.pdf", dirName,endingCent[i_cent]));
+  
+  cRatio->Clear();
+  cRatio->SetGridy(kFALSE);
+  hPionRatio->SetMarkerStyle(20);
+  hPionRatio->SetMarkerColor(2);
+  hPionRatio->SetLineColor(2);
+  hPionRatio->GetXaxis()->SetRangeUser(0.0, ptStop-0.001);
+  hPionRatio->GetYaxis()->SetRangeUser(0.0, 1.0);
+  hPionRatio->DrawCopy("P E");
+  hPionYield->SetMarkerStyle(20);
+  hPionYield->SetMarkerColor(2);
+  hPionYield->SetLineColor(2);
+  hPionYield->GetXaxis()->SetRangeUser(0.0, ptStop-0.001);
+  
+  hKaonRatio->SetMarkerStyle(20);
+  hKaonRatio->SetMarkerColor(kGreen);
+  hKaonRatio->SetLineColor(kGreen);
+  hKaonRatio->GetXaxis()->SetRangeUser(0.0, ptStop-0.001);
+  hKaonYield->SetMarkerStyle(20);
+  hKaonYield->SetMarkerColor(kGreen);
+  hKaonYield->SetLineColor(kGreen);
+  hKaonYield->GetXaxis()->SetRangeUser(0.0, ptStop-0.001);
+  hKaonRatio->DrawCopy("SAME P E");
+  
+  hProtonRatio->SetMarkerStyle(20);
+  hProtonRatio->SetMarkerColor(4);
+  hProtonRatio->SetLineColor(4);
+  hProtonRatio->GetXaxis()->SetRangeUser(0.0, ptStop-0.001);
+  hProtonYield->SetMarkerStyle(20);
+  hProtonYield->SetMarkerColor(4);
+  hProtonYield->SetLineColor(4);
+  hProtonYield->GetXaxis()->SetRangeUser(0.0, ptStop-0.001);
+  hProtonRatio->DrawCopy("SAME P E");
+  
+  hElectronRatio->SetMarkerStyle(20);
+  hElectronRatio->SetMarkerColor(6);
+  hElectronRatio->SetLineColor(6);
+  hElectronRatio->GetXaxis()->SetRangeUser(0.0, ptStop-0.001);
+  hElectronRatio->DrawCopy("SAME P E");
+  hElectronYield->SetMarkerStyle(20);
+  hElectronYield->SetMarkerColor(6);
+  hElectronYield->SetLineColor(6);
+  hElectronYield->GetXaxis()->SetRangeUser(0.0, ptStop-0.001);
+  
+  TH1D *hallRec=0;
+  
+  hallRec=(TH1D *)hPionRatio->Clone("hRecAll");
+  hallRec->Add(hKaonRatio);
+  hallRec->Add(hProtonRatio);
+  hallRec->SetLineColor(kYellow-3);
+  hallRec->SetMarkerColor(kYellow-3);
+  hallRec->DrawCopy("SAME P E");
+  
+  
+  gROOT->ProcessLine(".x drawText.C");
+  cRatio->SaveAs(Form("fit_results/%s_%s/particle_ratios.gif", dirName,endingCent[i_cent]));
+  cRatio->SaveAs(Form("fit_results/%s_%s/particle_ratios.pdf", dirName,endingCent[i_cent]));
+  
+  
+
+  if(outFileName) {
+    
+    CreateDir(Form("fit_results/fit_yields_results_%s",endingCent[i_cent]));
+    TFile* fileOut = new TFile(Form("fit_results/fit_yields_results_%s/%s%s.root", endingCent[i_cent], outFileName, ending[i_eta]), "RECREATE");
+    
+    
+    hPionRatio->Write();
+    hKaonRatio->Write();
+    hProtonRatio->Write();
+    hElectronRatio->Write();
+    
+    hPionYield->Write();
+    hKaonYield->Write();
+    hProtonYield->Write();
+    hElectronYield->Write();
+    
+    
+    
+    fElectronFraction->Write();
+    
+    for(Int_t bin = binStart; bin <= binStop; bin++){
+      
+      Int_t j = bin-binStart;
+      dEdxpoints[j]->Write();
+      fPion[j]->Write();
+      fKaon[j]->Write();
+      fProton[j]->Write();
+      fElectron[j]->Write();
+      fTotal[j]->Write();
+      
+      
+    }
+    fmip->Write();
+    fileOut->Close();
+  }
+  
+  
+  
+  
+}
+
+
+
+
+
+
+//____________________________________________________________________________
+void FitDeDxVsP(const Char_t * dedxFileName,
+               Double_t pStart, Double_t pStop,  Int_t i_cent,
+               Int_t optionDeDx, Int_t optionSigma,
+               Bool_t performFit = kFALSE,
+               Int_t etaLow=0, Int_t etaHigh=8, 
+               Bool_t fixAllParBB=kFALSE,
+               Bool_t fixAllParSigma=kFALSE,
+               Bool_t fixKtoZero=kFALSE,
+               const Char_t* initialFitFileName=0,
+               Int_t fixParametersDedx=-1,
+               Int_t fixParametersSigma=-1)
+{
+  //gStyle->SetOptStat(0);
+
+
+
+  TFile* dedxFile = FindFileFresh(dedxFileName);
+  if(!dedxFile)
+    return;
+
+  TList * list = (TList *)dedxFile->Get("outputdedx");
+  TH2D *hDeDxVsPPlus = 0;
+  hDeDxVsPPlus = (TH2D *)list->FindObject(Form("histAllCh%d%d", etaLow, etaHigh));
+  hDeDxVsPPlus->Sumw2();
+
+  TH2D *hDeDxVsPMinus = 0;
+  hDeDxVsPMinus = (TH2D *)list->FindObject(Form("histAllCh%d%d", etaHigh, etaLow));
+  hDeDxVsPMinus->Sumw2();
+
+  gSystem->Exec("rm *.root");
+  TFile *fplus = new TFile("hist2Dplus.root","recreate");
+  fplus->cd();
+  hDeDxVsPPlus->SetName(Form("histAllCh%d%d", etaLow, etaHigh));
+  hDeDxVsPPlus->Write();
+  fplus->Close();
+
+  TFile *fminus = new TFile("hist2Dminus.root","recreate");
+  fminus->cd();
+  hDeDxVsPMinus->SetName(Form("histAllCh%d%d", etaLow, etaHigh));
+  hDeDxVsPMinus->Write();
+  fminus->Close();
+
+  gSystem->Exec("hadd hist2Ddedx.root hist2Dplus.root hist2Dminus.root");
+
+  TFile *ffinal = TFile::Open("hist2Ddedx.root"); 
+  hDeDxVsP = (TH2D *)ffinal->Get(Form("histAllCh%d%d", etaLow, etaHigh));
+
+
+
+  CreateDir("old");
+  gSystem->Exec(Form("mv results/calibdedx_%d%d/* old/",etaLow, etaHigh));
+
+
+  TCanvas* cDeDxVsP = new TCanvas("cDeDxVsP", "dE/dx vs p", 400, 300);
+  cDeDxVsP->Clear();
+  cDeDxVsP->cd();
+  cDeDxVsP->SetLogz();
+  hDeDxVsP->SetTitle(0);
+  hDeDxVsP->Draw("COLZ");
+
+  TCanvas* cDeDxVsPLogX = new TCanvas("cDeDxVsPLogX", "dE/dx vs p", 400, 300);
+  cDeDxVsPLogX->Clear();
+  cDeDxVsPLogX->cd();
+  cDeDxVsPLogX->SetLogz();
+  cDeDxVsPLogX->SetLogx();
+  hDeDxVsP->Draw("COLZ");
+
+  // Root is a bit stupid with finidng bins so we have to add and subtract a
+  // little to be sure we get the right bin as we typically put edges as
+  // limits
+  const Int_t binStart = hDeDxVsP->GetXaxis()->FindBin(pStart+0.01);
+  pStart = hDeDxVsP->GetXaxis()->GetBinLowEdge(binStart);
+  const Int_t binStop  = hDeDxVsP->GetXaxis()->FindBin(pStop-0.01);
+  pStop = hDeDxVsP->GetXaxis()->GetBinUpEdge(binStop);
+  const Int_t nBins    = binStop - binStart + 1;
+
+  cout << "Doing 2d fit from pTlow = " << pStart << " (bin: " << binStart
+       << ") to pThigh = " << pStop << " (bin: " << binStop << ")" << endl;
+  
+  // Double_t binSize = (histdEdxvsP->GetXaxis()->GetXmax() - histdEdxvsP->GetXaxis()->GetXmin())/ histdEdxvsP->GetXaxis()->GetNbins();
+  
+  Double_t parDeDx[6]  = {0, 0, 0, 0, 0, 0};
+  Double_t parSigma[8] = {0, 0, 0, 0, 0, 0, 0, 0};
+  
+  const Int_t nLocalParams  = 3; // pi, K, p yields
+  Int_t nDeDxPar      = 0;
+  Int_t nSigmaPar     = 0;
+
+  switch(optionDeDx) {
+    
+  case 1:
+    nDeDxPar = 2;
+    parDeDx[0] = 39.7;
+    parDeDx[1] =  6.3;
+    break;
+  case 2:
+    nDeDxPar = 1;
+    parDeDx[0] =  7.3;
+    break;
+  case 3:
+    nDeDxPar = 2;
+    parDeDx[0] =  6.85097;
+    parDeDx[1] =  29.5933;
+    break;
+  case 4:
+    nDeDxPar = 3;
+    parDeDx[0] =  31.2435;
+    parDeDx[1] =  9.73037;
+    parDeDx[2] =  1.95832;
+    break;
+  case 5:
+    nDeDxPar = 4;
+    parDeDx[0] =  31.384;
+    parDeDx[1] =  9.5253;
+    parDeDx[2] =  7.3542;
+    parDeDx[3] =  1.4;
+    break;
+  case 6:
+    nDeDxPar = 5;
+    parDeDx[0] =  31.384;
+    parDeDx[1] =  9.5253;
+    parDeDx[2] =  7.3542;
+    parDeDx[3] =  1.5;
+    parDeDx[4] =  81;
+    break;
+  case 7:
+    nDeDxPar = 3;
+    parDeDx[0] =  31.3;
+    parDeDx[1] =  9.5;
+    parDeDx[2] =  1.4;
+    break;
+  default:
+
+    cout << "optionDeDx does not support option: " << optionDeDx << endl;
+    return;
+    break;
+  }
+
+  switch(optionSigma) {
+    
+  case 1:
+    nSigmaPar = 1;
+    parSigma[0] = 3.0;
+    break;
+  case 2:
+    nSigmaPar = 1;
+    parSigma[0] = 0.0655688;
+    break;
+  case 3:
+    nSigmaPar = 2;
+    parSigma[0] = 0.06;
+    parSigma[1] = -0.001;
+  case 4:
+    nSigmaPar = 2;
+    parSigma[0] = 0.06;
+    parSigma[1] = 1.0;
+    break;
+  case 5:
+    nSigmaPar = 2;
+    parSigma[0] = 0.06;
+    parSigma[1] = 1.0;
+    break;
+  case 6:
+    nSigmaPar = 2;
+    parSigma[0] = 0.06;
+    parSigma[1] = 0.0;
+    break;
+  case 7:
+    nSigmaPar = 3;
+    parSigma[0] = 0.06;
+    parSigma[1] = 0.0;
+    parSigma[2] = 2.0;
+    break;
+  case 8:
+    nSigmaPar = 3;
+    parSigma[0] = 0.06;
+    parSigma[1] = 0.0;
+    parSigma[2] = 2.0;
+    break;
+  case 9://works for 0-5 %
+    nSigmaPar = 3;
+    parSigma[0] = 1.88409e-01;
+    parSigma[1] = -3.84073e-03;
+    parSigma[2] = 3.03110e-05;
+    break;
+  case 10://works for 0-5 %
+    nSigmaPar = 6;
+    parSigma[0]=7.96380e-03;
+    parSigma[1]=8.09869e-05;
+    parSigma[2]=6.29707e-06;
+    parSigma[3]=-3.27295e-08;
+    parSigma[4]=-1.20200e+03;
+    parSigma[5]=-3.97089e+01;
+    break;
+  case 11://works for 0-5 %
+    nSigmaPar = 6;
+    parSigma[0]=7.96380e-03;
+    parSigma[1]=8.09869e-05;
+    parSigma[2]=6.29707e-06;
+    parSigma[3]=-3.27295e-08;
+    parSigma[4]=-1.20200e+03;
+    parSigma[5]=-3.97089e+01;
+    break;
+  case 12://works for 0-5 %
+    nSigmaPar = 3;
+    parSigma[0]=7.96380e-03;
+    parSigma[1]=8.09869e-05;
+    parSigma[2]=6.29707e-06;
+    break;
+  case 13://works for 0-5 %
+    nSigmaPar = 3;
+    parSigma[0]=7.96380e-03;
+    parSigma[1]=8.09869e-05;
+    parSigma[2]=6.29707e-06;
+    break;
+  case 14:
+    nSigmaPar = 2;
+    parSigma[0]=7.96380e-03;
+    parSigma[1]=8.09869e-05;
+    break;
+  case 15:
+    nSigmaPar = 3;
+    parSigma[0]=7.96380e-03;
+    parSigma[1]=8.09869e-05;
+    break;
+
+
+
+
+  default:
+
+    cout << "optionSigma does not support option: " << optionSigma << endl;
+    return;
+    break;
+  }
+
+  if(initialFitFileName) {
+    
+  
+
+    TFile* fresBB = FindFileFresh(initialFitFileName);
+    if(!fresBB)
+      return;
+    TF1 *fres=(TF1 *)fresBB->Get("sigmaRes");
+    TF1 *fBB=(TF1 *)fresBB->Get("dedxFunc");
+
+    Int_t nparres=fres->GetNpar();
+    Int_t nparBB=fBB->GetNpar();
+    for(Int_t ires=0;ires<nparres;++ires){
+      parSigma[ires] = fres->GetParameter(ires+1);
+      //cout<<fres->GetParameter(ires+1)<<endl;
+    }
+    for(Int_t iBB=0;iBB<nparBB;++iBB){
+      parDeDx[iBB] = fBB->GetParameter(iBB+1);
+      cout<<parDeDx[iBB]<<endl;
+
+    }
+    //Esto lo puse apenas
+    fixPlateau  = parDeDx[4];
+    fixMIP = fBB->Eval(3.5);
+    cout<<"Hello!!"<<"   fixPlateau="<<fixPlateau<<"  fixMIP="<<fixMIP<<endl;
+    
+    //fixMIP      = fitPar->MIP;
+    //return;
+  }
+
+
+  const Int_t nGlobalParams = 2  // binStart, binStop, 
+    + 2 + nDeDxPar               // optionDeDx, nDeDxPar, dedxpar0 ....
+    + 2 + nSigmaPar;             // nSigmaPar, sigmapar0 .....
+  
+  const Int_t nParams = nBins*nLocalParams + nGlobalParams;
+  
+  
+  TF2* fitAll = new TF2("fitAll", fitf3G, pStart, pStop, 30, 90, nParams);
+  Double_t parametersIn[nParams]; 
+  
+  parametersIn[0] = binStart;
+  fitAll->SetParName(0,"binStart");
+  fitAll->FixParameter(0, parametersIn[0]);
+
+  parametersIn[1] = binStop;
+  fitAll->SetParName(1,"binStop");
+  fitAll->FixParameter(1, parametersIn[1]);
+
+  // dE/dx parameters
+  parametersIn[2] = nDeDxPar;
+  fitAll->SetParName(2,"nDeDxPar");
+  fitAll->FixParameter(2, parametersIn[2]);
+
+  parametersIn[3] = optionDeDx;
+  fitAll->SetParName(3,"optionDeDx");
+  fitAll->FixParameter(3, parametersIn[3]);
+
+  for(Int_t i = 0; i < nDeDxPar; i++) {
+
+    parametersIn[4+i] = parDeDx[i];
+    fitAll->SetParName(4+i,Form("dEdXpar%d", i));
+    Int_t bit = TMath::Nint(TMath::Power(2, i));
+    if(fixParametersDedx>=0)
+      if (fixParametersDedx==0 || fixParametersDedx&bit) {
+       
+       cout << "Fixing dE/dx parameter " << i << endl;
+       fitAll->FixParameter(4+i, parametersIn[4+i]);
+      }
+    
+    if(fixAllParBB) {
+
+      fitAll->FixParameter(4+i, parametersIn[4+i]);
+      }
+  }
+
+  // sigma parameters
+  parametersIn[4+nDeDxPar] = nSigmaPar;
+  fitAll->SetParName(4+nDeDxPar,"nSigmaPar");
+  fitAll->FixParameter(4+nDeDxPar, parametersIn[4+nDeDxPar]);
+
+  parametersIn[5+nDeDxPar] = optionSigma;
+  fitAll->SetParName(5+nDeDxPar,"optionSigma");
+  fitAll->FixParameter(5+nDeDxPar, parametersIn[5+nDeDxPar]);
+
+  for(Int_t i = 0; i < nSigmaPar; i++) {
+
+    parametersIn[6+nDeDxPar+i] = parSigma[i];
+    fitAll->SetParName(6+nDeDxPar+i,Form("sigmapar%d", i));
+    Int_t bit = TMath::Nint(TMath::Power(2, i));
+    if(fixParametersSigma>=0)
+      if (fixParametersSigma==0 || fixParametersSigma&bit) {
+       
+       cout << "Fixing sigma parameter " << i << endl;
+       fitAll->FixParameter(6+nDeDxPar+i, parametersIn[6+nDeDxPar+i]);
+      }
+    
+    if(fixAllParSigma) {
+      
+      fitAll->FixParameter(6+nDeDxPar+i, parametersIn[6+nDeDxPar+i]);
+      }
+  }
+  
+  // Initial yields
+
+  for(Int_t bin = binStart; bin <= binStop; bin++) {
+    
+    TH1D* hDeDxVsPProj =(TH1D*)hDeDxVsP->ProjectionY("hDeDxVsPProj", bin, bin);
+    
+    const Int_t offset = nGlobalParams + (bin - binStart)*nLocalParams;
+    const Double_t all = hDeDxVsPProj->Integral();
+
+    fitAll->SetParName(offset + 0, Form("piYield%d", bin));
+    fitAll->SetParName(offset + 1, Form("kYield%d", bin));
+    fitAll->SetParName(offset + 2, Form("pYield%d", bin));
+    fitAll->SetParLimits(offset + 0, 0, 10*all);
+    fitAll->SetParLimits(offset + 2, 0, 10*all);
+
+    if(fixKtoZero) {
+      parametersIn[offset + 0] = (all)*0.5;
+      parametersIn[offset + 1] = (all)*0.0;
+      fitAll->FixParameter(offset + 1, 0.0);
+      parametersIn[offset + 2] = (all)*0.5;
+    } else {
+      parametersIn[offset + 0] = (all)*0.6;
+      parametersIn[offset + 1] = (all)*0.3;
+      fitAll->SetParLimits(offset + 1, 0, 10*all);
+      parametersIn[offset + 2] = (all)*0.1;
+    }    
+    // fitAll->SetParLimits(offset + 0, 0., histdEdxvsPpy->GetEntries());
+    // fitAll->SetParLimits(offset + 1, 0., histdEdxvsPpy->GetEntries());
+    // fitAll->SetParLimits(offset + 2, 0., histdEdxvsPpy->GetEntries());    
+  }
+  
+  fitAll->SetParameters(parametersIn);
+  fitAll->Print();
+
+  Bool_t converged = kFALSE;
+
+  if(performFit) {
+    for(Int_t i = 0; i < 4; i++) {
+      TFitResultPtr fitResultPtr =  hDeDxVsP->Fit(fitAll, "0S", "", pStart+0.01, pStop-0.01);
+      if(!fitResultPtr->Status()) {
+       //      if(!fitResultPtr->Status() && !fitResultPtr->CovMatrixStatus()) {
+       
+       converged = kTRUE;
+       break;
+      }
+    }
+  }
+  // else we just draw how the results looks with the input parameters
+
+  Double_t parametersOut[nParams];
+  fitAll->GetParameters(parametersOut);
+  const Double_t* parameterErrorsOut = fitAll->GetParErrors();
+
+  // overlay the fitfunction
+  
+
+  TF1* fDeDxPi = new TF1("fDeDxPi", FitFunc, 0, 50, nDeDxPar+1); // +1 because of option! 
+  fDeDxPi->SetParameters(&parametersOut[3]);
+  fDeDxPi->SetLineWidth(2);
+  cDeDxVsP->cd();
+  fDeDxPi->Draw("SAME");
+  cDeDxVsPLogX->cd();
+  fDeDxPi->Draw("SAME");
+
+  TF1* fDeDxK = new TF1("fDeDxK", FitFunc, 0, 50, nDeDxPar+1); 
+  fDeDxK->SetParameters(&parametersOut[3]);
+  fDeDxK->SetParameter(0, fDeDxK->GetParameter(0)+10);
+  fDeDxK->SetLineWidth(2);
+  cDeDxVsP->cd();
+  fDeDxK->Draw("SAME");
+  cDeDxVsPLogX->cd();
+  fDeDxK->Draw("SAME");
+
+  TF1* fDeDxP = new TF1("fDeDxP", FitFunc, 0, 50, nDeDxPar+1); 
+  fDeDxP->SetParameters(&parametersOut[3]);
+  fDeDxP->SetParameter(0, fDeDxP->GetParameter(0)+20);
+  fDeDxP->SetLineWidth(2);
+  cDeDxVsP->cd();
+  fDeDxP->Draw("SAME");
+  cDeDxVsPLogX->cd();
+  fDeDxP->Draw("SAME");
+
+  TF1* fDeDxE = new TF1("fDeDxE", FitFunc, 0, 50, nDeDxPar+1); 
+  fDeDxE->SetParameters(&parametersOut[3]);
+  fDeDxE->SetParameter(0, fDeDxE->GetParameter(0)+30);
+  fDeDxE->SetLineWidth(2);
+  cDeDxVsP->cd();
+  fDeDxE->Draw("SAME");
+  cDeDxVsPLogX->cd();
+  fDeDxE->Draw("SAME");
+
+  TF1* fSigma = new TF1("fSigma", SigmaFunc, 0, 50, nSigmaPar+1); 
+  fSigma->SetParameters(&parametersOut[5 + nDeDxPar]);
+
+  //  fitAll->Draw("same cont3"); 
+
+  CreateDir(Form("results/%s/calibdedx_%d%d",endingCent[i_cent],etaLow,etaHigh));
+
+  //CreateDir(Form("results/calibdedx_%d%d",etaLow,etaHigh));
+
+  cDeDxVsP->cd();
+  cDeDxVsP->Modified();
+  cDeDxVsP->Update();
+  gROOT->ProcessLine(".x drawText.C");
+  cDeDxVsP->SaveAs(Form("results/%s/calibdedx_%d%d/dedx_vs_p.gif",endingCent[i_cent],etaLow,etaHigh));
+  cDeDxVsP->SaveAs(Form("results/%s/calibdedx_%d%d/dedx_vs_p.pdf",endingCent[i_cent],etaLow,etaHigh));
+
+  cDeDxVsPLogX->cd();
+  cDeDxVsPLogX->Modified();
+  cDeDxVsPLogX->Update();
+  gROOT->ProcessLine(".x drawText.C");
+  cDeDxVsPLogX->SaveAs(Form("results/%s/calibdedx_%d%d/dedx_vs_p_logx.gif",endingCent[i_cent],etaLow,etaHigh));
+  cDeDxVsPLogX->SaveAs(Form("results/%s/calibdedx_%d%d/dedx_vs_p_logx.pdf",endingCent[i_cent],etaLow,etaHigh));
+
+  //cross check
+  TCanvas* cFits = new TCanvas("cFits", "Fit comparison to data", 1200, 800);
+
+  cFits->Clear();
+  cFits->Divide(7, 5);
+
+  TF1* pion = new TF1("pion", "gausn", 30, 90);
+  pion->SetLineWidth(2);
+  pion->SetLineColor(kRed);
+  TF1* kaon = new TF1("kaon", "gausn", 30, 90);
+  kaon->SetLineWidth(2);
+  kaon->SetLineColor(kGreen);
+  TF1* proton = new TF1("proton", "gausn", 30, 90);
+  proton->SetLineWidth(2);
+  proton->SetLineColor(kBlue);
+  TF1* total = new TF1("total", "gausn(0)+gausn(3)+gausn(6)", 30, 90);
+  total->SetLineColor(kBlack);
+  total->SetLineWidth(2);
+  total->SetLineStyle(2);
+
+  TLegend* legend = new TLegend(0.11, 0.6, 0.35, 0.85);    
+  legend->SetBorderSize(0);
+  legend->SetFillColor(0);
+  legend->AddEntry(total, "3-Gauss fit", "L");
+  legend->AddEntry(pion, "#pi", "L");
+  legend->AddEntry(kaon, "K", "L");
+  legend->AddEntry(proton, "p", "L");
+
+
+  TCanvas* cSingleFit = new TCanvas("cSingleFit", "single fit");
+  cSingleFit->SetLogy(1);
+
+  TH1D* hPionRatio =(TH1D*)hDeDxVsP->ProjectionX("hPionRatio", 1, 1);
+  hPionRatio->Reset();
+  hPionRatio->GetXaxis()->SetRangeUser(pStart+0.001, pStop-0.001);
+  hPionRatio->GetYaxis()->SetRangeUser(0.0, 1.0);
+  hPionRatio->SetTitle("particle fractions; p [GeV/c]; particle fraction");
+  TH1D* hKaonRatio   = (TH1D*)hPionRatio->Clone("hKaonRatio");
+  TH1D* hProtonRatio = (TH1D*)hPionRatio->Clone("hProtonRatio");
+  
+  for(Int_t bin = binStart; bin <= binStop; bin++){
+    
+    cout << "Making projection for bin: " << bin << endl;
+    
+    const Int_t j = bin-binStart;
+    cFits->cd();
+    cFits->cd(j + 1);
+    
+    TH1D* hDeDxVsPProj =(TH1D*)hDeDxVsP->ProjectionY(Form("hDeDxVsPProj%d", bin), bin, bin);
+    //    hDeDxVsPProj->GetXaxis()->SetRangeUser(40, 90);
+    hDeDxVsPProj->SetTitle(Form("%.1f<p<%.1f GeV/c", 
+                               hDeDxVsP->GetXaxis()->GetBinLowEdge(bin),
+                               hDeDxVsP->GetXaxis()->GetBinUpEdge(bin)));
+    hDeDxVsPProj->Draw();
+    
+    const Int_t offset = nGlobalParams + j*nLocalParams; 
+    const Double_t p = hDeDxVsP->GetXaxis()->GetBinCenter(bin);
+    const Double_t pKeff = p*piMass/kMass; // corresponding p of a pion with same dE/dx
+    const Double_t pPeff = p*piMass/pMass; // corresponding p of a pion with same dE/dx
+    const Double_t meanDeDxPi = fDeDxPi->Eval(p);
+    const Double_t meanDeDxK  = fDeDxPi->Eval(pKeff);
+    const Double_t meanDeDxP  = fDeDxPi->Eval(pPeff);
+    Double_t gausParams[9] = { 
+      parametersOut[offset + 0], 
+      meanDeDxPi, 
+      fSigma->Eval(meanDeDxPi) ,
+      parametersOut[offset + 1], 
+      meanDeDxK, 
+      fSigma->Eval(meanDeDxK) ,
+      parametersOut[offset + 2], 
+      meanDeDxP, 
+      fSigma->Eval(meanDeDxP) ,
+    };
+
+    for(Int_t i = 0; i < 9; i++) 
+      cout << gausParams[i] << ", ";
+
+    cout << endl;
+    
+    pion->SetParameters(&gausParams[0]);
+    pion->DrawCopy("same");
+    Double_t all =  hDeDxVsPProj->Integral();
+    hPionRatio->SetBinContent(bin, parametersOut[offset + 0]/all);
+    hPionRatio->SetBinError(bin, parameterErrorsOut[offset + 0]/all);
+
+    kaon->SetParameters(&gausParams[3]);
+    kaon->DrawCopy("same");
+    hKaonRatio->SetBinContent(bin, parametersOut[offset + 1]/all);
+    hKaonRatio->SetBinError(bin, parameterErrorsOut[offset + 1]/all);
+    
+    proton->SetParameters(&gausParams[6]);
+    proton->DrawCopy("same");
+    hProtonRatio->SetBinContent(bin, parametersOut[offset + 2]/all);
+    hProtonRatio->SetBinError(bin, parameterErrorsOut[offset + 2]/all);
+    
+    total->SetParameters(gausParams);
+    total->DrawCopy("same");
+
+    cSingleFit->cd();
+    cSingleFit->Clear();
+    //    cSingleFit->SetLogy();
+    hDeDxVsPProj->Draw();
+    pion->DrawCopy("same");
+    kaon->DrawCopy("same");
+    proton->DrawCopy("same");
+    total->DrawCopy("same");
+    
+    gROOT->ProcessLine(".x drawText.C(2)");
+
+
+    cSingleFit->SaveAs(Form("results/%s/calibdedx_%d%d/ptspectrum_bin%d.gif",endingCent[i_cent], etaLow,etaHigh, bin));
+    cSingleFit->SaveAs(Form("results/%s/calibdedx_%d%d/ptspectrum_bin%d.pdf",endingCent[i_cent],etaLow, etaHigh, bin));
+    //    legend->Draw();
+
+
+  }
+
+  TCanvas* cRatio = new TCanvas("cRatio", "ratios/all vs p", 600, 400);
+  cRatio->Clear();
+  hPionRatio->SetMaximum(0.8);
+  hPionRatio->SetMarkerStyle(20);
+  hPionRatio->SetMarkerColor(2);
+  hPionRatio->Draw("P E");
+
+  hKaonRatio->SetMarkerStyle(20);
+  hKaonRatio->SetMarkerColor(3);
+  hKaonRatio->Draw("SAME P E");
+
+  hProtonRatio->SetMarkerStyle(20);
+  hProtonRatio->SetMarkerColor(4);
+  hProtonRatio->Draw("SAME P E");
+  gROOT->ProcessLine(".x drawText.C(2)");
+  cRatio->SaveAs(Form("results/%s/calibdedx_%d%d/particle_ratios.gif",endingCent[i_cent],etaLow,etaHigh));
+  cRatio->SaveAs(Form("results/%s/calibdedx_%d%d/particle_ratios.pdf",endingCent[i_cent],etaLow,etaHigh));
+
+
+
+  //
+  // Store the <dE/dx> parameters in a file that we can get them back to use for the Delta-pi!
+  //
+  DeDxFitInfo* fitInfo = new DeDxFitInfo();
+  fitInfo->MIP     = fixMIP;
+  //fitInfo->plateau = 80; 
+  fitInfo->optionDeDx = optionDeDx; 
+  fitInfo->nDeDxPar = nDeDxPar; 
+  for(Int_t i = 0; i < nDeDxPar; i++) {
+    fitInfo->parDeDx[i] = fDeDxPi->GetParameter(i+1); // 1st parameter is option
+  }
+  //fitInfo->plateau = fDeDxPi->GetParameter(5);//lo puse 16/09/13
+  fitInfo->optionSigma = optionSigma; 
+  fitInfo->nSigmaPar = nSigmaPar; 
+  for(Int_t i = 0; i < nSigmaPar; i++) {
+    cout<<"my sugma param="<<fSigma->GetParameter(i+1)<<endl;
+    fitInfo->parSigma[i] = fSigma->GetParameter(i+1); // 1st parameter is option 
+  }
+
+
+  fitInfo->Print();
+
+
+  if(converged) {
+
+    cout << "Fit converged and error matrix was ok" << endl;
+  } else {
+
+    cout << "WARNING!!!!!!!!!!!!!!! Fit did not converge, or error matrix was  not ok!" << endl;
+  }
+
+
+  CreateDir(Form("fitparameters/%s",endingCent[i_cent]));
+  //cout<<"flag 0"<<"    gSystem->BaseName(calibFileName))="<< gSystem->BaseName(calibFileName)<<endl;
+  TFile* outFile = new TFile(Form("fitparameters/%s/%d%d_dataPbPb.root", endingCent[i_cent], etaLow, etaHigh), "RECREATE");
+  outFile->cd();
+  fitInfo->Write("fitInfo");
+  cout<<"flag 1"<<endl;
+
+  outFile->Close();
+
+  if(converged) {
+
+    cout << "Fit converged and error matrix was ok" << endl;
+  } else {
+
+    cout << "WARNING!!!!!!!!!!!!!!! Fit did not converge, or error matrix was  not ok!" << endl;
+  }
+}
+
+
+
+
+
+void MakeFitsV0s(const Char_t * fileInputName, const Char_t * fileE, const Char_t *outDir, Int_t index_eta){
+  /*
+  gStyle->SetPalette(1);
+  gStyle->SetCanvasColor(10);
+  gStyle->SetFrameFillColor(10);
+  gStyle->SetOptStat(0);
+  gStyle->SetOptTitle(0);
+  gStyle->SetOptFit(0);
+  */
+  TFile *f12=TFile::Open(fileE);
+  TF1 *fPlateau = (TF1 *)f12->Get(Form("fGauss%d",index_eta));
+  /*
+  TH1D * hEV0s=(TH1D *)f12->Get("hdedx_vs_p_e");
+  TGraphErrors *gResE=(TGraphErrors    *)f12->Get("gRes_Electrons");
+  */
+
+  const Char_t * endNamePos[4]={"02","24","46","68"};
+  const Char_t * endNameNeg[4]={"20","42","64","86"};
+
+
+
+  TFile *f1=TFile::Open(fileInputName);
+  TList * list = (TList *)f1->Get("outputdedx");
+
+
+  //positive eta
+  TH2D *hPi2DPos=(TH2D *)list->FindObject(Form("histPiTof%s",endNamePos[index_eta]));
+  TH2D *hPi2DV0sPos=(TH2D *)list->FindObject(Form("histPiV0%s",endNamePos[index_eta]));
+  TH2D *hP2DPos=(TH2D *)list->FindObject(Form("histPV0%s",endNamePos[index_eta]));
+  TH1D *hPi2DpPos=(TH1D *)list->FindObject(Form("histPiTof1D%s",endNamePos[index_eta]));
+  TH1D *hPi2DV0spPos=(TH1D *)list->FindObject(Form("histPiV01D%s",endNamePos[index_eta]));
+  TH1D *hP2DpPos=(TH1D *)list->FindObject(Form("histPV01D%s",endNamePos[index_eta]));
+
+
+  gSystem->Exec("rm *.root");
+
+  TFile *fpos = new TFile("histtmppos.root","recreate");
+  fpos->cd();
+  hPi2DPos->SetName(Form("histPiTof%s",endNamePos[index_eta]));
+  hPi2DV0sPos->SetName(Form("histPiV0%s",endNamePos[index_eta]));
+  hP2DPos->SetName(Form("histPV0%s",endNamePos[index_eta]));
+  hPi2DpPos->SetName(Form("histPiTof1D%s",endNamePos[index_eta]));
+  hPi2DV0spPos->SetName(Form("histPiV01D%s",endNamePos[index_eta]));
+  hP2DpPos->SetName(Form("histPV01D%s",endNamePos[index_eta]));
+  hPi2DPos->Write();
+  hPi2DV0sPos->Write();
+  hP2DPos->Write();
+  hPi2DpPos->Write();
+  hPi2DV0spPos->Write();
+  hP2DpPos->Write();
+  fpos->Close();
+
+
+  //negative eta
+  TH2D *hPi2DNeg=(TH2D *)list->FindObject(Form("histPiTof%s",endNameNeg[index_eta]));
+  TH2D *hPi2DV0sNeg=(TH2D *)list->FindObject(Form("histPiV0%s",endNameNeg[index_eta]));
+  TH2D *hP2DNeg=(TH2D *)list->FindObject(Form("histPV0%s",endNameNeg[index_eta]));
+  TH1D *hPi2DpNeg=(TH1D *)list->FindObject(Form("histPiTof1D%s",endNameNeg[index_eta]));
+  TH1D *hPi2DV0spNeg=(TH1D *)list->FindObject(Form("histPiV01D%s",endNameNeg[index_eta]));
+  TH1D *hP2DpNeg=(TH1D *)list->FindObject(Form("histPV01D%s",endNameNeg[index_eta]));
+
+
+
+  TFile *fneg = new TFile("histtmpneg.root","recreate");
+  fneg->cd();
+  hPi2DNeg->SetName(Form("histPiTof%s",endNamePos[index_eta]));
+  hPi2DV0sNeg->SetName(Form("histPiV0%s",endNamePos[index_eta]));
+  hP2DNeg->SetName(Form("histPV0%s",endNamePos[index_eta]));
+  hPi2DpNeg->SetName(Form("histPiTof1D%s",endNamePos[index_eta]));
+  hPi2DV0spNeg->SetName(Form("histPiV01D%s",endNamePos[index_eta]));
+  hP2DpNeg->SetName(Form("histPV01D%s",endNamePos[index_eta]));
+  hPi2DNeg->Write();
+  hPi2DV0sNeg->Write();
+  hP2DNeg->Write();
+  hPi2DpNeg->Write();
+  hPi2DV0spNeg->Write();
+  hP2DpNeg->Write();
+  fneg->Close();
+
+  //pos+neg eta
+  TH2D *hPi2D=0;
+  TH2D *hPi2DV0s=0;
+  TH2D *hP2D=0;
+  TH1D *hPi2Dp=0;
+  TH1D *hPi2DV0sp=0;
+  TH1D *hP2Dp=0;
+
+
+  gSystem->Exec(Form("hadd tmp%d.root histtmppos.root histtmpneg.root",index_eta));
+
+  TFile *ffinalv0 = TFile::Open(Form("tmp%d.root",index_eta));
+  hPi2D=(TH2D *)ffinalv0->Get(Form("histPiTof%s",endNamePos[index_eta]));
+  hPi2DV0s=(TH2D *)ffinalv0->Get(Form("histPiV0%s",endNamePos[index_eta]));
+  hP2D=(TH2D *)ffinalv0->Get(Form("histPV0%s",endNamePos[index_eta]));
+  hPi2Dp=(TH1D *)ffinalv0->Get(Form("histPiTof1D%s",endNamePos[index_eta]));
+  hPi2DV0sp=(TH1D *)ffinalv0->Get(Form("histPiV01D%s",endNamePos[index_eta]));
+  hP2Dp=(TH1D *)ffinalv0->Get(Form("histPV01D%s",endNamePos[index_eta]));
+
+
+
+  TList *l1=new TList();
+  l1->SetOwner();
+
+
+  TGraphErrors* graphRes = new TGraphErrors();
+  Int_t nERes = 0;
+
+  TGraphErrors* graphBB = new TGraphErrors();
+  Int_t nBB = 0;
+
+  TGraphErrors* graphBBpitof = new TGraphErrors();
+  Int_t nBBpitof = 0;
+  TGraphErrors* graphBBpiv0s = new TGraphErrors();
+  Int_t nBBpiv0s = 0;
+  TGraphErrors* graphBBpv0s = new TGraphErrors();
+  Int_t nBBpv0s = 0;
+
+
+
+  TF1* fFit1D = new TF1("fFit1D", "gausn(0)", 40, 90);
+
+
+  //TF1* pionsv0s = new TF1("pionsv0s", "gausn", 56, 80);
+  TF1* pionsv0s = new TF1("pionsv0s", "gausn", 40, 80);
+  pionsv0s->SetLineColor(kRed);
+  pionsv0s->SetLineWidth(2);
+  pionsv0s->SetLineStyle(2);
+
+  //TF1* protonsv0s = new TF1("protonsv0s", "gausn", 40, 55);
+  TF1* protonsv0s = new TF1("protonsv0s", "gausn", 40, 80);
+  protonsv0s->SetLineColor(kBlue);
+  protonsv0s->SetLineWidth(2);
+  protonsv0s->SetLineStyle(2);
+
+
+  TF1* total = new TF1("total", "gausn(0)+gausn(3)", 40, 80);
+  total->SetLineColor(1);
+  total->SetLineWidth(2);
+  total->SetLineStyle(2);
+
+
+
+  Int_t nbins=hPi2D->GetXaxis()->GetNbins();
+
+
+
+  for(Int_t bin=1;bin<=nbins;++bin){
+    
+
+    TH1D *hproy_pion=0;
+    TH1D *hproy_proton=0;
+
+    Double_t p=hPi2D->GetXaxis()->GetBinCenter(bin);
+    Double_t lowp = hPi2D->GetXaxis()->GetBinLowEdge(bin);
+    Double_t widthp = hPi2D->GetXaxis()->GetBinWidth(bin);
+    hPi2Dp->GetXaxis()->SetRangeUser(lowp,lowp+widthp);
+    Double_t meanp= hPi2Dp->GetMean();
+    Double_t meanpE= hPi2Dp->GetMeanError();
+
+
+    cout<<"p="<<p<<endl;
+
+    if(p>2.0&&p<9.0){
+      
+      cout<<bin<<endl;
+
+      hproy_pion=(TH1D *)hPi2D->ProjectionY(Form("hproy_pionTof_%d",bin),bin,bin);
+
+      cout<<bin<<endl;
+
+      hproy_pion->SetTitle(Form("pionsv0s: %.1f < p < %.1f GeV/c; dE/dx (GeV/c); Entries",
+                               hPi2D->GetXaxis()->GetBinLowEdge(bin),
+                               hPi2D->GetXaxis()->GetBinUpEdge(bin)));
+      cout<<bin<<endl;
+      hproy_proton=(TH1D *)hP2D->ProjectionY(Form("hproy_proton_%d",bin),bin,bin);
+
+      cout<<"nbinsproj="<<hproy_proton->GetNbinsX()<<endl;
+      //hproy_pion->Draw();
+
+
+      TSpectrum *s = new TSpectrum(2);
+      Int_t nfound = s->Search(hproy_pion,2,"",0.15);
+      cout<<"!!!!!!!!!!!!!!         nfound="<<nfound<<endl;
+      Int_t npeaks = 0;
+      Float_t *xpeaks = s->GetPositionX();
+
+
+      for (Int_t p1=0;p1<nfound;p1++) {
+       //Float_t xp = xpeaks[p1];
+       npeaks++;
+      }
+
+      if(npeaks==0)
+       continue;
+      if(npeaks==1){
+
+       fFit1D->SetParameters(hproy_pion->GetEntries(), xpeaks[0], hproy_pion->GetRMS());
+       hproy_pion->Fit(fFit1D, "Q");
+       hproy_pion->Fit(fFit1D, "Q", "", xpeaks[0]-hproy_pion->GetRMS(),
+                       xpeaks[0]+hproy_pion->GetRMS());
+       
+       hproy_pion->Fit(fFit1D, "R", "", fFit1D->GetParameter(1)-1.0*fFit1D->GetParameter(2),
+                       fFit1D->GetParameter(1)+2*fFit1D->GetParameter(2));
+
+
+
+
+         graphBB->SetPoint(nBB,meanp/0.14, fFit1D->GetParameter(1));
+         graphBB->SetPointError(nBB,meanpE, fFit1D->GetParError(1));
+         nBB++;
+         
+         graphBBpitof->SetPoint(nBBpitof,meanp/0.14, fFit1D->GetParameter(1));
+         graphBBpitof->SetPointError(nBBpitof,meanpE, fFit1D->GetParError(1));
+         nBBpitof++;
+         
+         if(p>3&&p<9){
+           graphRes->SetPoint(nERes,fFit1D->GetParameter(1),fFit1D->GetParameter(2)/fFit1D->GetParameter(1));
+           graphRes->SetPointError(nERes,fFit1D->GetParError(1),fFit1D->GetParError(2)/fFit1D->GetParameter(1));
+           nERes++;
+         }
+         
+
+
+
+      }
+      else if(npeaks==2){
+       Int_t bin1 = hproy_pion->GetXaxis()->FindBin(xpeaks[0]);
+       Int_t bin2 = hproy_pion->GetXaxis()->FindBin(xpeaks[1]);
+
+       Int_t bin1C = hproy_pion->GetBinContent(bin1);
+       Int_t bin2C = hproy_pion->GetBinContent(bin2);
+
+       cout<<"bin1="<<bin1<<"  content="<<bin1C<<endl;
+       cout<<"bin2="<<bin2<<"  content="<<bin2C<<endl;
+
+       
+       if(bin==24){
+         
+         hproy_pion->Fit(fFit1D, "R", "", 65, 85);
+         graphBB->SetPoint(nBB,meanp/0.14, fFit1D->GetParameter(1));
+         graphBB->SetPointError(nBB,meanpE, fFit1D->GetParError(1));
+         nBB++;
+         
+         graphBBpitof->SetPoint(nBBpitof,meanp/0.14, fFit1D->GetParameter(1));
+         graphBBpitof->SetPointError(nBBpitof,meanpE, fFit1D->GetParError(1));
+         nBBpitof++;
+           
+         
+         
+       }
+       
+
+       else if(bin1C>bin2C){//pion signal larger
+         
+         
+         fFit1D->SetParameters(hproy_pion->GetEntries(), xpeaks[0], hproy_pion->GetRMS());
+         hproy_pion->Fit(fFit1D, "Q");
+         hproy_pion->Fit(fFit1D, "Q", "", xpeaks[0]-hproy_pion->GetRMS(),
+                         xpeaks[0]+hproy_pion->GetRMS());
+         
+         
+         
+         
+         hproy_pion->Fit(fFit1D, "R", "", fFit1D->GetParameter(1)-1.0*fFit1D->GetParameter(2),
+                         fFit1D->GetParameter(1)+2.5*fFit1D->GetParameter(2));
+
+         
+         
+         
+         cout<<"hello"<<endl;
+         //protonsv0s->SetParameters(hproy_proton->GetEntries(), xpeaks[1], hproy_proton->GetRMS());
+         
+         
+         hproy_proton->Fit(protonsv0s, "Q");
+
+         
+         hproy_proton->Fit(protonsv0s, "Q", "", xpeaks[1]-hproy_proton->GetRMS(),
+                           xpeaks[1]+hproy_proton->GetRMS());
+         
+         
+         hproy_proton->Fit(protonsv0s, "R", "", xpeaks[1]-2*protonsv0s->GetParameter(2),
+                         xpeaks[1]+2*protonsv0s->GetParameter(2));
+         
+         protonsv0s->SetParameters(hproy_proton->GetEntries(), xpeaks[1], protonsv0s->GetParameter(2));
+         hproy_proton->Fit(protonsv0s, "R", "", xpeaks[1]-2.0*protonsv0s->GetParameter(2),
+                           xpeaks[1]+2.0*protonsv0s->GetParameter(2));
+         
+         total->SetParameter(1,protonsv0s->GetParameter(1));
+         total->SetParameter(2,protonsv0s->GetParameter(2));
+         total->SetParameter(4,fFit1D->GetParameter(1));
+         total->SetParameter(5,fFit1D->GetParameter(2));
+         hproy_pion->Fit(total,"R");
+
+         Float_t reducedCh2=total->GetChisquare()/total->GetNDF();
+         if(reducedCh2>1.0){
+           hproy_pion->Fit(fFit1D, "R", "", fFit1D->GetParameter(1)-2*fFit1D->GetParameter(2),
+                           fFit1D->GetParameter(1)+2*fFit1D->GetParameter(2));
+           //hproy_pion->Fit(fFit1D, "R", "", fFit1D->GetParameter(1)-0.5*fFit1D->GetParameter(2),
+           hproy_pion->Fit(fFit1D, "R", "", fFit1D->GetParameter(1)-1.0*fFit1D->GetParameter(2),
+                           fFit1D->GetParameter(1)+2*fFit1D->GetParameter(2));
+
+
+           graphBB->SetPoint(nBB,meanp/0.14, fFit1D->GetParameter(1));
+           graphBB->SetPointError(nBB,meanpE, fFit1D->GetParError(1));
+           nBB++;
+           
+
+           graphBBpitof->SetPoint(nBBpitof,meanp/0.14, fFit1D->GetParameter(1));
+           graphBBpitof->SetPointError(nBBpitof,meanpE, fFit1D->GetParError(1));
+           nBBpitof++;
+           
+           if(p>3.0&&p<9){
+             graphRes->SetPoint(nERes,fFit1D->GetParameter(1),fFit1D->GetParameter(2)/fFit1D->GetParameter(1));
+             graphRes->SetPointError(nERes,fFit1D->GetParError(1),fFit1D->GetParError(2)/fFit1D->GetParameter(1));
+             nERes++;
+           }
+           
+           
+
+         }else{
+           
+           graphBB->SetPoint(nBB,meanp/0.14, total->GetParameter(4));
+           graphBB->SetPointError(nBB,meanpE, total->GetParError(4));
+           nBB++;
+           
+           graphBBpitof->SetPoint(nBBpitof,meanp/0.14, total->GetParameter(4));
+           graphBBpitof->SetPointError(nBBpitof,meanpE, total->GetParError(4));
+           nBBpitof++;
+
+           if(p>3.0&&p<9){
+             graphRes->SetPoint(nERes,total->GetParameter(4),total->GetParameter(5)/total->GetParameter(4));
+             
+             graphRes->SetPointError(nERes,total->GetParError(4),total->GetParError(5)/total->GetParameter(4));
+             nERes++;
+           }
+           
+           
+           
+           
+         }
+         
+         
+         
+         
+       }else{
+         
+         fFit1D->SetParameters(hproy_pion->GetEntries(), xpeaks[1], hproy_pion->GetRMS());
+         hproy_pion->Fit(fFit1D, "Q");
+         hproy_pion->Fit(fFit1D, "Q", "", xpeaks[1]-hproy_pion->GetRMS(),
+                         xpeaks[0]+hproy_pion->GetRMS());
+         
+         hproy_pion->Fit(fFit1D, "R", "", fFit1D->GetParameter(1)-2*fFit1D->GetParameter(2),
+                         fFit1D->GetParameter(1)+2*fFit1D->GetParameter(2));
+         
+
+         protonsv0s->SetParameters(hproy_proton->GetEntries(), xpeaks[0], hproy_proton->GetRMS());
+         hproy_proton->Fit(protonsv0s, "Q");
+         hproy_proton->Fit(protonsv0s, "Q", "", xpeaks[0]-hproy_proton->GetRMS(),
+                         xpeaks[0]+hproy_proton->GetRMS());
+         
+         hproy_proton->Fit(protonsv0s, "R", "", xpeaks[0]-2*protonsv0s->GetParameter(2),
+                         xpeaks[0]+2*protonsv0s->GetParameter(2));
+
+         protonsv0s->SetParameters(hproy_proton->GetEntries(), xpeaks[0], protonsv0s->GetParameter(2));
+         hproy_proton->Fit(protonsv0s, "R", "", xpeaks[0]-2*protonsv0s->GetParameter(2),
+                         xpeaks[0]+2*protonsv0s->GetParameter(2));
+
+
+         //graphBB->SetPoint(nBB,p/0.14, xpeaks[1]);
+         //graphBB->SetPointError(nBB,0, fFit1D->GetParError(1));
+         //nBB++;
+         total->SetParameter(1,protonsv0s->GetParameter(1));
+         total->SetParameter(2,protonsv0s->GetParameter(2));
+         total->SetParameter(4,fFit1D->GetParameter(1));
+         total->SetParameter(5,fFit1D->GetParameter(2));
+         hproy_pion->Fit(total,"R");
+
+         //graphBB->SetPoint(nBB,p/0.14, xpeaks[1]);
+         //graphBB->SetPointError(nBB,0, total->GetParError(4));
+         graphBB->SetPoint(nBB,meanp/0.14, fFit1D->GetParameter(1));
+         graphBB->SetPointError(nBB,meanpE, fFit1D->GetParError(1));
+         nBB++;
+
+         graphBBpitof->SetPoint(nBBpitof,meanp/0.14, fFit1D->GetParameter(1));
+         graphBBpitof->SetPointError(nBBpitof,meanpE, fFit1D->GetParError(1));
+         nBBpitof++;
+         
+
+         if(p>3.0&&p<9){
+           graphRes->SetPoint(nERes,fFit1D->GetParameter(1),fFit1D->GetParameter(2)/fFit1D->GetParameter(1));
+           graphRes->SetPointError(nERes,fFit1D->GetParError(1),fFit1D->GetParError(2)/fFit1D->GetParameter(1));
+           nERes++;
+         }
+
+
+
+       }
+
+
+
+    
+
+
+
+
+
+
+
+      }
+      
+      l1->Add(hproy_pion);
+      
+      
+    }
+    else
+      continue;
+    
+  }
+
+
+
+  graphRes->SetTitle(";#LT dE/dx  #GT;#sigma/#LT dE/dx  #GT");
+  graphRes->SetName("gRes_PionsTof");
+  
+  graphBBpitof->SetTitle(";#beta*#gamma; #LT dE/dx  #GT");
+  graphBBpitof->SetName("BB_piTOF");
+  graphBBpitof->SetMarkerColor(kGreen+2);
+  graphBBpitof->SetLineColor(kGreen+2);
+  graphBBpitof->SetMarkerStyle(21);
+
+  l1->Add(graphRes);
+  l1->Add(graphBBpitof);
+
+  
+
+  //pi by V0s
+
+  TGraphErrors* graphResV0s = new TGraphErrors();
+  nERes = 0;
+
+  nbins=hPi2DV0s->GetXaxis()->GetNbins();
+
+
+
+  for(Int_t bin=1;bin<=nbins;++bin){
+    
+  
+
+    TH1D *hproy_pionV0s=0;
+    TH1D *hproy_proton=0;
+    Double_t p=hPi2DV0s->GetXaxis()->GetBinCenter(bin);
+
+    Double_t lowp = hPi2DV0s->GetXaxis()->GetBinLowEdge(bin);
+    Double_t widthp = hPi2DV0s->GetXaxis()->GetBinWidth(bin);
+    hPi2DV0sp->GetXaxis()->SetRangeUser(lowp,lowp+widthp);
+    Double_t meanp= hPi2DV0sp->GetMean();
+    Double_t meanpE= hPi2DV0sp->GetMeanError();
+
+
+
+    if(p>2.0&&p<9){
+
+      hproy_pionV0s=(TH1D *)hPi2DV0s->ProjectionY(Form("hproy_pionV0s_%d",bin),bin,bin);
+      hproy_pionV0s->SetTitle(Form("pionsv0s: %.1f < p < %.1f GeV/c; dE/dx (GeV/c); Entries",
+                               hPi2DV0s->GetXaxis()->GetBinLowEdge(bin),
+                               hPi2DV0s->GetXaxis()->GetBinUpEdge(bin)));
+
+     fFit1D->SetParameters(hproy_pionV0s->GetEntries(), hproy_pionV0s->GetMean(), hproy_pionV0s->GetRMS());
+      hproy_pionV0s->Fit(fFit1D, "Q");
+      hproy_pionV0s->Fit(fFit1D, "Q", "", hproy_pionV0s->GetMean()-hproy_pionV0s->GetRMS(),
+                        hproy_pionV0s->GetMean()+hproy_pionV0s->GetRMS());
+      
+      
+      //      hproy_pionV0s->Fit(fFit1D, "Q", "", fFit1D->GetParameter(1)-1.5*fFit1D->GetParameter(2),
+      //             fFit1D->GetParameter(1)+3.0*fFit1D->GetParameter(2));
+
+      hproy_pionV0s->Fit(fFit1D, "R", "", fFit1D->GetParameter(1)-1.0*fFit1D->GetParameter(2),
+                        fFit1D->GetParameter(1)+2.0*fFit1D->GetParameter(2));
+      
+      
+      hproy_proton=(TH1D *)hP2D->ProjectionY(Form("hproy_proton_%d",bin),bin,bin);
+      hproy_proton->SetTitle(Form("protonsv0s: %.1f < p < %.1f GeV/c; dE/dx; Entries",
+                                 hPi2D->GetXaxis()->GetBinLowEdge(bin),
+                                 hPi2D->GetXaxis()->GetBinUpEdge(bin)));
+      
+      //protonsv0s->SetParameters(hproy_proton->GetEntries(), hproy_proton->GetMean(), hproy_proton->GetRMS());
+      hproy_proton->Fit(protonsv0s,"R");
+      //hproy_proton->Fit(protonsv0s,"R", "",protonsv0s->GetParameter(1)-3.0*protonsv0s->GetParameter(2),protonsv0s->GetParameter(1)+1.0*protonsv0s->GetParameter(2));
+      hproy_proton->Fit(protonsv0s,"R","",protonsv0s->GetParameter(1)-2.0*protonsv0s->GetParameter(2),protonsv0s->GetParameter(1)+0.8*protonsv0s->GetParameter(2));
+
+      total->SetParameter(0,0.1*hproy_pionV0s->GetEntries());
+      total->SetParameter(1,protonsv0s->GetParameter(1));
+      total->SetParameter(2,protonsv0s->GetParameter(2));
+      total->SetParameter(3,0.9*hproy_pionV0s->GetEntries());
+      total->SetParameter(4,fFit1D->GetParameter(1));
+      total->SetParameter(5,fFit1D->GetParameter(2));
+      
+      hproy_pionV0s->Fit(total,"0R");
+      
+      cout<<"p="<<p<<"   yield0="<<total->GetParameter(0)<<"  yield1="<<total->GetParameter(3)<<endl;
+      cout<<"p="<<p<<"   mean0="<<total->GetParameter(1)<<"  mean1="<<total->GetParameter(4)<<endl;
+      cout<<"p="<<p<<"   sigma0="<<total->GetParameter(1)<<"  sigma1="<<total->GetParameter(5)<<endl;
+      
+      total->SetParameter(1,total->GetParameter(1));
+      total->SetParameter(2,total->GetParameter(2));
+      total->SetParameter(4,total->GetParameter(4));
+      total->SetParameter(5,total->GetParameter(5));
+      
+      hproy_pionV0s->Fit(total,"R");
+
+      
+      if(p>4.0){
+
+
+       //return;
+
+       graphResV0s->SetPoint(nERes,total->GetParameter(4),total->GetParameter(5)/total->GetParameter(4));
+       graphResV0s->SetPointError(nERes,total->GetParError(4),total->GetParError(5)/total->GetParameter(4));
+       nERes++;
+       
+       //graphBB->SetPoint(nBB,p/0.14, total->GetParameter(4));
+       //graphBB->SetPointError(nBB,0, total->GetParError(4));
+       //nBB++;
+       
+      
+       graphBBpiv0s->SetPoint(nBBpiv0s,meanp/0.14, total->GetParameter(4));
+       graphBBpiv0s->SetPointError(nBBpiv0s,meanpE,total->GetParError(4));
+       nBBpiv0s++;
+       
+       l1->Add(hproy_pionV0s);
+       
+      }
+      
+      
+      
+    }
+    
+  }
+  
+  
+  graphResV0s->SetTitle(";#LT dE/dx  #GT;#sigma/#LT dE/dx  #GT");
+  graphResV0s->SetName("gRes_PionsV0sArmenteros");
+  l1->Add(graphResV0s);
+
+  graphBBpiv0s->SetTitle(";#beta*#gamma; #LT dE/dx  #GT");
+  graphBBpiv0s->SetName("BB_piV0s");
+  graphBBpiv0s->SetMarkerColor(2);
+  graphBBpiv0s->SetLineColor(2);
+  graphBBpiv0s->SetMarkerStyle(24);
+  l1->Add(graphBBpiv0s);
+  
+  //protons
+  
+  
+
+  
+  TGraphErrors* graphPRes = new TGraphErrors();
+  nERes = 0;
+  TGraphErrors* graphdEdxvsP = new TGraphErrors();
+  TGraphErrors* graphSigmavsP = new TGraphErrors();
+  graphdEdxvsP->SetName("graphdEdxvsP_p");
+  graphSigmavsP->SetName("graphSigmavsP_p");
+
+  nbins=hPi2D->GetXaxis()->GetNbins();
+  for(Int_t bin=1;bin<=nbins;++bin){
+    
+    TH1D *hproy_pion=0;
+    TH1D *hproy_proton=0;
+
+    
+    Double_t p=hPi2D->GetXaxis()->GetBinCenter(bin);
+    Double_t lowp = hP2D->GetXaxis()->GetBinLowEdge(bin);
+    Double_t widthp = hP2D->GetXaxis()->GetBinWidth(bin);
+    hP2Dp->GetXaxis()->SetRangeUser(lowp,lowp+widthp);
+    Double_t meanp= hP2Dp->GetMean();
+    Double_t meanpE= hP2Dp->GetMeanError();
+
+
+    if(p<2.0||p>5.0)
+      continue;
+
+
+    if(p>2){
+
+      cout<<"p="<<hPi2D->GetXaxis()->GetBinCenter(bin)<<endl;
+
+      hproy_pion=(TH1D *)hPi2D->ProjectionY(Form("hproy_pion_%d",bin),bin,bin);
+      hproy_pion->SetTitle(Form("pionsv0s: %.1f < p < %.1f GeV/c; dE/dx (GeV/c); Entries",
+                               hPi2D->GetXaxis()->GetBinLowEdge(bin),
+                               hPi2D->GetXaxis()->GetBinUpEdge(bin)));
+      
+      hproy_pion->Fit(pionsv0s,"R");
+      
+      hproy_proton=(TH1D *)hP2D->ProjectionY(Form("hproy_proton_%d",bin),bin,bin);
+      hproy_proton->SetTitle(Form("protonsv0s: %.1f < p < %.1f GeV/c; dE/dx (GeV/c); Entries",
+                                 hPi2D->GetXaxis()->GetBinLowEdge(bin),
+                                 hPi2D->GetXaxis()->GetBinUpEdge(bin)));
+      
+
+
+      hproy_proton->Fit(protonsv0s,"R");
+
+      hproy_proton->Fit(protonsv0s,"R","",protonsv0s->GetParameter(1)-2.0*protonsv0s->GetParameter(2),protonsv0s->GetParameter(1)+0.8*protonsv0s->GetParameter(2));
+
+
+      Float_t reducedCh2=protonsv0s->GetChisquare()/protonsv0s->GetNDF();
+      //if(reducedCh2<3)
+      if(reducedCh2>3){
+       total->SetParameter(1,protonsv0s->GetParameter(1));
+       total->SetParameter(2,protonsv0s->GetParameter(2));
+       total->SetParameter(4,pionsv0s->GetParameter(1));
+       total->SetParameter(5,pionsv0s->GetParameter(2));
+       
+       hproy_proton->Fit(total,"R");
+      }
+      //if(total->GetParameter(1)<=0)continue; 
+
+      if(p>2&&p<7){
+       if(reducedCh2<5){
+
+
+
+         if(p>3.0&&p<5){
+           graphPRes->SetPoint(nERes,protonsv0s->GetParameter(1),protonsv0s->GetParameter(2)/protonsv0s->GetParameter(1));
+           graphPRes->SetPointError(nERes,protonsv0s->GetParError(1),protonsv0s->GetParError(2)/protonsv0s->GetParameter(1));
+           nERes++;
+         }
+         graphBB->SetPoint(nBB,meanp/0.943, protonsv0s->GetParameter(1));
+         graphBB->SetPointError(nBB,meanpE,protonsv0s->GetParError(1));
+
+         graphBBpv0s->SetPoint(nBBpv0s,meanp/0.943, protonsv0s->GetParameter(1));
+         graphBBpv0s->SetPointError(nBBpv0s,meanpE,protonsv0s->GetParError(1));
+         nBBpv0s++;
+
+         graphdEdxvsP->SetPoint(nBB,p, protonsv0s->GetParameter(1));
+         graphdEdxvsP->SetPointError(nBB,0.001, protonsv0s->GetParError(1));
+         graphSigmavsP->SetPoint(nBB,p, protonsv0s->GetParameter(2));
+         graphSigmavsP->SetPointError(nBB,0.0001, protonsv0s->GetParError(2));
+         nBB++;
+
+       }
+       else{
+
+         if(p>3.0&&p<5){
+           graphPRes->SetPoint(nERes,total->GetParameter(1),total->GetParameter(2)/total->GetParameter(1));
+           graphPRes->SetPointError(nERes,total->GetParError(1),total->GetParError(2)/total->GetParameter(1));
+           nERes++;
+         }
+
+         //graphBB->SetPoint(nBB,p/0.943, total->GetParameter(1));
+         //graphBB->SetPointError(nBB,0,total->GetParError(1));
+
+         graphBBpv0s->SetPoint(nBBpv0s,meanp/0.943, protonsv0s->GetParameter(1));
+         graphBBpv0s->SetPointError(nBBpv0s,meanpE,protonsv0s->GetParError(1));
+         nBBpv0s++;
+
+         
+         graphdEdxvsP->SetPoint(nBB,p, total->GetParameter(1));
+         graphdEdxvsP->SetPointError(nBB,0.001, total->GetParError(1));
+         graphSigmavsP->SetPoint(nBB,p, total->GetParameter(2));
+         graphSigmavsP->SetPointError(nBB,0.0001, total->GetParError(2));
+
+         graphBB->SetPoint(nBB,meanp/0.943, protonsv0s->GetParameter(1));
+         graphBB->SetPointError(nBB,meanpE,protonsv0s->GetParError(1));
+
+         nBB++;
+       }
+
+       
+      }
+      /*else{
+       
+       graphPRes->SetPoint(nERes,protonsv0s->GetParameter(1),protonsv0s->GetParameter(2)/protonsv0s->GetParameter(1));
+       graphPRes->SetPointError(nERes,protonsv0s->GetParError(1),protonsv0s->GetParError(2)/protonsv0s->GetParameter(1));
+       
+       
+       nERes++;
+       
+       }*/
+
+      l1->Add(hproy_proton);
+
+
+    }
+    
+  }
+
+  graphPRes->SetTitle(";#LT dE/dx  #GT;#sigma/#LT dE/dx  #GT");
+  graphPRes->SetName("gRes_Protons");
+
+  l1->Add(graphPRes);
+
+  graphBBpv0s->SetTitle(";#beta*#gamma; #LT dE/dx  #GT");
+  graphBBpv0s->SetName("BB_pv0s");
+  graphBBpv0s->SetMarkerColor(4);
+  graphBBpv0s->SetLineColor(4);
+  graphBBpv0s->SetMarkerStyle(24);
+  l1->Add(graphBBpv0s);
+  
+  
+
+  graphBB->SetTitle(";#beta*#gamma; #LT dE/dx  #GT");
+  graphBB->SetName("gBB");
+
+  Double_t dedxPar[10]  = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+
+  dedxPar[0] =  0;
+  dedxPar[1] =  31.384;
+  dedxPar[2] =  9.5253;
+  dedxPar[3] =  7.3542;
+  dedxPar[4] =  1.5;
+  
+  TF1 *dedxFunc = new TF1("dedxFunc", FitFunc, 2.5, 60, 6);
+  dedxFunc->SetParameters(dedxPar);
+  // Fit betagamma (option+40)
+  dedxFunc->SetParameter(0, 6+40);
+  dedxFunc->FixParameter(0, 6+40);
+  //parameters for MB
+  dedxFunc->SetParameter(1, 33.7538);
+  dedxFunc->SetParameter(2, 8.05926);
+  dedxFunc->SetParameter(3, 1.69954);
+  dedxFunc->SetParameter(4, 1.37793);
+  dedxFunc->SetParameter(5, 78.0465);
+
+  cout << "!!!!!!!!!!!!!!!!!                    Setting new plateau=" << fPlateau->GetParameter(1) <<  endl;
+  dedxFunc->SetParameter(5, fPlateau->GetParameter(1));
+  dedxFunc->FixParameter(5, fPlateau->GetParameter(1));
+
+  graphBB->Fit(dedxFunc, "0R", "SAME");
+  graphBB->Fit(dedxFunc, "0R", "SAME");
+  graphBB->Fit(dedxFunc, "0R", "SAME");
+
+  graphBB->Draw();
+
+  for(Int_t i=0;i<dedxFunc->GetNpar();++i)
+    dedxFunc->SetParameter(i,dedxFunc->GetParameter(i));
+  
+  l1->Add(dedxFunc);
+  l1->Add(graphBB);
+
+
+
+  //Fit to get resolution
+  TGraphErrors *g1=0;
+  g1=new TGraphErrors();
+  g1->SetName("res_allpions");
+
+  Int_t npionts=0;
+
+  for(Int_t ipoint = 0; ipoint < graphPRes->GetN(); ipoint++){
+    g1->SetPoint(npionts,graphPRes->GetX()[ipoint],graphPRes->GetY()[ipoint]);
+    g1->SetPointError(npionts,graphPRes->GetEX()[ipoint],graphPRes->GetEY()[ipoint]);
+    npionts++;
+  }
+  
+  //  pions TOF
+  
+  for(Int_t ipoint = 0; ipoint < graphRes->GetN(); ipoint++){
+    g1->SetPoint(npionts,graphRes->GetX()[ipoint],graphRes->GetY()[ipoint]);
+    g1->SetPointError(npionts,graphRes->GetEX()[ipoint],graphRes->GetEY()[ipoint]);
+    npionts++;
+    }
+  //electrons
+  g1->SetPoint(npionts,fPlateau->GetParameter(1),fPlateau->GetParameter(2)/fPlateau->GetParameter(1));
+  g1->SetPointError(npionts,fPlateau->GetParError(1),fPlateau->GetParError(2)/fPlateau->GetParameter(1));
+
+  Float_t sigmaPar[7];
+  sigmaPar[0]=12;
+  sigmaPar[1]=7.96380e-03;
+  sigmaPar[2]=8.09869e-05;
+  sigmaPar[3]=6.29707e-06;
+  sigmaPar[4]=-3.27295e-08;
+  sigmaPar[5]=-1.20200e+03;
+  sigmaPar[6]=-3.97089e+01;
+
+
+  TF1 *sigmarrrrr=new TF1("sigmaRes",SigmaFunc,47.5, 85,4);
+  sigmarrrrr->FixParameter(0,sigmaPar[0]);
+  g1->Fit(sigmarrrrr, "0R", "SAME");
+  g1->Fit(sigmarrrrr, "0R", "SAME");
+  g1->Fit(sigmarrrrr, "0R", "SAME");
+
+  TF1 *fdedxvspP=new TF1("fdedxvspP","pol0",3.0,5.5);
+  graphdEdxvsP->Fit(fdedxvspP, "R", "SAME");
+  TF1 *fsigmavspP=new TF1("fsigmavspP","pol0",3.0,5.5);
+  graphSigmavsP->Fit(fsigmavspP, "R", "SAME");
+
+
+
+  for(Int_t i=0;i<6;++i)
+    sigmarrrrr->SetParameter(i,sigmarrrrr->GetParameter(i));
+
+  l1->Add(graphdEdxvsP);
+  l1->Add(graphSigmavsP);
+
+  l1->Add(fdedxvspP);
+  l1->Add(fsigmavspP);
+
+  l1->Add(sigmarrrrr);
+  l1->Add(g1);
+
+  TFile *fout=new TFile(Form("%s/hres_0_5_%s.root",outDir,endNamePos[index_eta]),"recreate");
+  fout->cd();
+  l1->Write();
+  fout->Close();
+
+
+
+
+}
+
+
+
+void MakeFitsExternalData(const Char_t * inFile, const Char_t * outDir){
+
+
+  TFile *fIn = TFile::Open(inFile);  
+  TList * list = (TList *)fIn->Get("outputdedx");
+  //plateau, electrons 0.4<p<0.6 GeV/c
+  TH2D *hEVsEta = (TH2D *)list->FindObject("hPlateauVsEta");
+  hEVsEta->Sumw2();
+  //four th1, projections
+  TH1D *hdEdxE[4]={0,0,0,0};// |eta|
+  TH1D *hdEdxEpos[4]={0,0,0,0};// eta>0
+  TF1 *fGaussE[4]={0,0,0,0};
+  Int_t index_max = 1;
+  for(Int_t i=0; i<4; ++i){
+
+    Int_t index_min_negeta = i+index_max;
+    Int_t index_max_negeta = i+index_max+1;
+    Int_t index_min_poseta = 17-(i+index_max+1);
+    Int_t index_max_poseta = 17-(i+index_max);
+
+    hdEdxE[i] = (TH1D *)hEVsEta->ProjectionY(Form("hdEdxNegEta%d",i),index_min_negeta,index_max_negeta);
+    hdEdxEpos[i] = (TH1D *)hEVsEta->ProjectionY(Form("hdEdxPosEta%d",i),index_min_poseta,index_max_poseta);
+    hdEdxE[i]->Add(hdEdxEpos[i]);
+
+    //Fitting the signal
+    fGaussE[i]=new TF1(Form("fGauss%d",i),"gaus");
+    hdEdxE[i]->Fit(fGaussE[i],"0","");
+    hdEdxE[i]->Fit(fGaussE[i],"","", fGaussE[i]->GetParameter(1)-1.5*fGaussE[i]->GetParameter(2), fGaussE[i]->GetParameter(1)+1.5*fGaussE[i]->GetParameter(2));
+
+    index_max++;
+
+  }
+
+  CreateDir(outDir);
+
+  TFile *fout = 0;
+
+  if(outDir){
+
+    fout = new TFile(Form("%s/PrimaryElectrons.root",outDir),"recreate");
+    fout->cd();
+    for(Int_t i=0; i<4; ++i){
+      hdEdxE[i]->Write();
+      fGaussE[i]->Write();
+    }
+    fout->Close();
+
+
+  }
+
+
+
+}
+//__________________________________________
+void PlotQA(const Char_t * inFile){
+
+  //gStyle->SetPalette(1);
+  gStyle->SetCanvasColor(10);
+  gStyle->SetFrameFillColor(10);
+  gStyle->SetOptStat(0);
+  gStyle->SetOptTitle(0);
+  gStyle->SetOptFit(0);
+  
+  const Char_t* ending[8] = { "86", "64", "42", "20", "02", "24", "46", "68" };
+  const Char_t* LatexEta[8] = { "-0.8<#eta<-0.6", "-0.6<#eta<-0.4", "-0.4<#eta<-0.2", "-0.2<#eta<0", 
+                               "0<#eta<0.2", "0.2<#eta<0.4", "0.4<#eta<0.6", "0.6<#eta<0.8" };
+  TFile *fIn = TFile::Open(inFile);
+  if(!fIn)
+    return;
+
+  TList * list = (TList *)fIn->Get("outputdedx");
+
+
+  TH2D *hMIP[8];
+  TH2D *hPlateau[8];
+  TProfile *pMIP[8];
+  TProfile *pPlateau[8];
+
+  TH2D *hMIPVsNch[8];
+  TProfile *pMIPVsNch[8];
+
+
+  //Secondary pions
+  TH2D *hPiV0s2D[8];
+  TH1D *hPion3[8];
+  TH1D *hPion2GeV = new TH1D("hPion2GeV","hPion2GeV",8,-0.8,0.8);
+
+
+  for(Int_t index_eta = 0; index_eta < 8; ++index_eta){
+    hMIP[index_eta] = (TH2D *)list->FindObject(Form("hMIPVsPhi%s",ending[index_eta]));
+    hPlateau[index_eta] = (TH2D *)list->FindObject(Form("hPlateauVsPhi%s",ending[index_eta]));
+    pMIP[index_eta] = (TProfile *)list->FindObject(Form("pMIPVsPhi%s",ending[index_eta]));
+    pPlateau[index_eta] = (TProfile *)list->FindObject(Form("pPlateauVsPhi%s",ending[index_eta]));
+    hPiV0s2D[index_eta] = (TH2D *)list->FindObject(Form("histPiV0%s",ending[index_eta]));
+    hPion3[index_eta] = (TH1D *)hPiV0s2D[index_eta]->ProjectionY(Form("Pion3%s",ending[index_eta]),16,16);
+
+    hMIPVsNch[index_eta] = (TH2D *)list->FindObject(Form("hMIPVsNch%s",ending[index_eta]));
+    pMIPVsNch[index_eta] = (TProfile *)list->FindObject(Form("pMIPVsNch%s",ending[index_eta]));
+
+
+
+    TF1 *fgaus = 0;
+    fgaus = new TF1("fgaus","gausn",40,80);
+    hPion3[index_eta]->Fit(fgaus,"0","",40,80);
+
+    hPion2GeV->SetBinContent(index_eta+1,fgaus->GetParameter(1));
+    hPion2GeV->SetBinError(index_eta+1,fgaus->GetParError(1));
+  }
+
+
+
+  TH2D *hMIPVsEta = (TH2D *)list->FindObject("hMIPVsEta");
+  TProfile *pMIPVsEta = (TProfile *)list->FindObject("pMIPVsEta");
+  TProfile *pMIPVsEtaV0s = (TProfile *)list->FindObject("pMIPVsEtaV0s");
+  TProfile *pPlateauVsEta = (TProfile *)list->FindObject("pPlateauVsEta");
+  TH2D *hPlateauVsEta = (TH2D *)list->FindObject("hPlateauVsEta");
+
+
+
+  const Int_t nPadX = 4;
+  const Int_t nPadY = 2;
+  Float_t noMargin    = 0.000;
+  Float_t topMargin   = 0.01;
+  Float_t botMargin   = 0.1;
+  Float_t leftMargin  = 0.05;
+  Float_t rightMargin = 0.01;
+  Float_t width       = (1-rightMargin-leftMargin)/nPadX;
+  Float_t height      = (1-botMargin-topMargin)/nPadY;
+  Float_t shift = 0.05;
+    
+  
+  TCanvas* cMIP = new TCanvas("cMIP2", "Raa Pions", 900, 500);
+  TPad* padMIP[nPadX*nPadY];
+  //Float_t shift = 0.1;
+  //Float_t shift = 0.05;
+  const char* yTitleMIP = "d#it{E}/d#it{x}_{MIP}";
+  const char* xTitleMIP = "#phi (rad)";
+  cMIP->cd();
+  
+  
+  TLatex* latex = new TLatex();
+  latex->SetNDC();
+  latex->SetTextAlign(22);
+  latex->SetTextAngle(0);
+  latex->SetTextFont(62);
+  
+  latex->DrawLatex(0.53,0.04,xTitleMIP);
+  latex->SetTextAngle(90);
+  latex->SetTextSize(0.055);
+  latex->DrawLatex(0.025,0.5,yTitleMIP);
+  latex->SetTextSize(0.05);
+
+  
+  for (Int_t i = 1; i < 9; i++) {
+    
+    Int_t iX = (i-1)%nPadX;
+    Int_t iY = (i-1)/nPadX;
+    Float_t x1 = leftMargin + iX*width;
+    if(iX==2)
+      x1-=0.01;
+    else if(iX==1)
+      x1+=0.01;
+    Float_t x2 = leftMargin + (iX + 1)*width;
+     if(iX==0)
+       x2+=0.01;
+     else if(iX==1)
+       x2-=0.01;
+     Float_t y1 = 1 - topMargin - (iY +1)*height;
+     Float_t y2 = 1 - topMargin - iY*height;
+     padMIP[i-1] = new TPad(Form("padRaa%d",i),"", x1, y1, x2, y2, 0, 0, 0);
+     Float_t mBot = noMargin; 
+     Float_t mTop = noMargin; 
+     Float_t mLeft = noMargin; 
+     Float_t mRight = noMargin; 
+     if(iY==0)       mTop   = shift;
+     if(iY==nPadY-1) mBot   = shift;
+     if(iX==0)       mLeft  = 0.08;
+     if(iX==nPadX-1) mRight = 0.08;
+     padMIP[i-1]->SetBottomMargin(mBot);
+     padMIP[i-1]->SetTopMargin(mTop);
+     padMIP[i-1]->SetLeftMargin(mLeft);
+     padMIP[i-1]->SetRightMargin(mRight);
+     
+     cMIP->cd(); 
+     
+     padMIP[i-1]->Draw();
+     
+  }
+
+  latex->SetTextAngle(0);
+  latex->SetTextSize(0.085);
+  
+  Int_t index_dec = 7;
+
+  for(Int_t index_eta = 0; index_eta < 8; ++index_eta){
+    padMIP[index_eta]->cd();
+    padMIP[index_eta]->SetTickx(1);
+    padMIP[index_eta]->SetTicky(1);
+    padMIP[index_eta]->SetLogy(0);
+    padMIP[index_eta]->SetLogx(0);
+    hMIP[index_eta]->GetXaxis()->SetLabelSize(0.06);
+    hMIP[index_eta]->GetYaxis()->SetLabelSize(0.06);
+    hMIP[index_eta]->GetXaxis()->SetTitleSize(0);
+    hMIP[index_eta]->GetYaxis()->SetTitleSize(0);
+
+    if(index_eta<4){
+      hMIP[index_eta]->Draw("colz"); 
+      pMIP[index_eta]->Draw("samep"); 
+      latex->DrawLatex(0.5,0.1,LatexEta[index_eta]);
+    }
+    else{
+      hMIP[index_dec]->Draw("colz");
+      pMIP[index_dec]->Draw("samep"); 
+      latex->DrawLatex(0.5,0.1,LatexEta[index_dec]);
+      index_dec--;
+    }
+
+  }
+
+
+
+  cMIP->SaveAs("MIPvsPhi.png");
+
+  //Plateau vs phi, different eta intervals 
+  TCanvas* cPlateau = new TCanvas("cPlateau2", "Raa Pions", 900, 500);
+  TPad* padPlateau[nPadX*nPadY];
+  const char* yTitlePlateau = "d#it{E}/d#it{x}_{Plateau}";
+  const char* xTitlePlateau = "#phi (rad)";
+  cPlateau->cd();
+  
+  
+  //TLatex* latex = new TLatex();
+  latex->SetNDC();
+  latex->SetTextAlign(22);
+  latex->SetTextAngle(0);
+  latex->SetTextFont(62);
+  
+  latex->DrawLatex(0.53,0.04,xTitlePlateau);
+  latex->SetTextAngle(90);
+  latex->SetTextSize(0.055);
+  latex->DrawLatex(0.025,0.5,yTitlePlateau);
+  latex->SetTextSize(0.05);
+
+  
+  for (Int_t i = 1; i < 9; i++) {
+    
+    Int_t iX = (i-1)%nPadX;
+    Int_t iY = (i-1)/nPadX;
+    Float_t x1 = leftMargin + iX*width;
+    if(iX==2)
+      x1-=0.01;
+    else if(iX==1)
+      x1+=0.01;
+    Float_t x2 = leftMargin + (iX + 1)*width;
+     if(iX==0)
+       x2+=0.01;
+     else if(iX==1)
+       x2-=0.01;
+     Float_t y1 = 1 - topMargin - (iY +1)*height;
+     Float_t y2 = 1 - topMargin - iY*height;
+     padPlateau[i-1] = new TPad(Form("padPlateau%d",i),"", x1, y1, x2, y2, 0, 0, 0);
+     Float_t mBot = noMargin; 
+     Float_t mTop = noMargin; 
+     Float_t mLeft = noMargin; 
+     Float_t mRight = noMargin; 
+     if(iY==0)       mTop   = shift;
+     if(iY==nPadY-1) mBot   = shift;
+     if(iX==0)       mLeft  = 0.08;
+     if(iX==nPadX-1) mRight = 0.08;
+     padPlateau[i-1]->SetBottomMargin(mBot);
+     padPlateau[i-1]->SetTopMargin(mTop);
+     padPlateau[i-1]->SetLeftMargin(mLeft);
+     padPlateau[i-1]->SetRightMargin(mRight);
+     
+     cPlateau->cd(); 
+     
+     padPlateau[i-1]->Draw();
+     
+  }
+
+  latex->SetTextAngle(0);
+  latex->SetTextSize(0.085);
+  
+  index_dec = 7;
+
+  for(Int_t index_eta = 0; index_eta < 8; ++index_eta){
+    padPlateau[index_eta]->cd();
+    padPlateau[index_eta]->SetTickx(1);
+    padPlateau[index_eta]->SetTicky(1);
+    padPlateau[index_eta]->SetLogy(0);
+    padPlateau[index_eta]->SetLogx(0);
+    hPlateau[index_eta]->GetXaxis()->SetLabelSize(0.06);
+
+    hPlateau[index_eta]->GetYaxis()->SetLabelSize(0.06);
+    hPlateau[index_eta]->GetXaxis()->SetTitleSize(0);
+    hPlateau[index_eta]->GetYaxis()->SetTitleSize(0);
+
+    if(index_eta<4){
+      hPlateau[index_eta]->Draw("colz"); 
+      pPlateau[index_eta]->Draw("samep"); 
+      latex->DrawLatex(0.5,0.1,LatexEta[index_eta]);
+    }
+    else{
+      hPlateau[index_dec]->Draw("colz");
+      pPlateau[index_dec]->Draw("samep"); 
+      latex->DrawLatex(0.5,0.1,LatexEta[index_dec]);
+      index_dec--;
+    }
+
+  }
+
+
+
+  cPlateau->SaveAs("PlateauVsPhi.png");
+
+
+  //MIP vs Nch
+  TCanvas* cMIPVsNch = new TCanvas("cMIPVsNch2", "Raa Pions", 900, 500);
+  TPad* padMIPVsNch[nPadX*nPadY];
+  const char* yTitleMIPVsNch = "d#it{E}/d#it{x}_{MIP}";
+  const char* xTitleMIPVsNch = "TPC-track multiplicity";
+  cMIPVsNch->cd();
+  
+  
+  latex->SetNDC();
+  latex->SetTextAlign(22);
+  latex->SetTextAngle(0);
+  latex->SetTextFont(62);
+  
+  latex->DrawLatex(0.53,0.04,xTitleMIPVsNch);
+  latex->SetTextAngle(90);
+  latex->SetTextSize(0.055);
+  latex->DrawLatex(0.025,0.5,yTitleMIPVsNch);
+  latex->SetTextSize(0.05);
+
+  
+  for (Int_t i = 1; i < 9; i++) {
+    
+    Int_t iX = (i-1)%nPadX;
+    Int_t iY = (i-1)/nPadX;
+    Float_t x1 = leftMargin + iX*width;
+    if(iX==2)
+      x1-=0.01;
+    else if(iX==1)
+      x1+=0.01;
+    Float_t x2 = leftMargin + (iX + 1)*width;
+     if(iX==0)
+       x2+=0.01;
+     else if(iX==1)
+       x2-=0.01;
+     Float_t y1 = 1 - topMargin - (iY +1)*height;
+     Float_t y2 = 1 - topMargin - iY*height;
+     padMIPVsNch[i-1] = new TPad(Form("padMIPVsNch%d",i),"", x1, y1, x2, y2, 0, 0, 0);
+     Float_t mBot = noMargin; 
+     Float_t mTop = noMargin; 
+     Float_t mLeft = noMargin; 
+     Float_t mRight = noMargin; 
+     if(iY==0)       mTop   = shift;
+     if(iY==nPadY-1) mBot   = shift;
+     if(iX==0)       mLeft  = 0.08;
+     if(iX==nPadX-1) mRight = 0.08;
+     padMIPVsNch[i-1]->SetBottomMargin(mBot);
+     padMIPVsNch[i-1]->SetTopMargin(mTop);
+     padMIPVsNch[i-1]->SetLeftMargin(mLeft);
+     padMIPVsNch[i-1]->SetRightMargin(mRight);
+     
+     if(i>=5&&i<9)
+       padMIPVsNch[i-1]->SetBottomMargin(mBot+0.05);
+
+     cMIPVsNch->cd(); 
+     
+     padMIPVsNch[i-1]->Draw();
+     
+  }
+
+  latex->SetTextAngle(0);
+  latex->SetTextSize(0.085);
+  
+  index_dec = 7;
+
+  for(Int_t index_eta = 0; index_eta < 8; ++index_eta){
+    padMIPVsNch[index_eta]->cd();
+
+    padMIPVsNch[index_eta]->SetTickx(1);
+    padMIPVsNch[index_eta]->SetTicky(1);
+    padMIPVsNch[index_eta]->SetLogy(0);
+    padMIPVsNch[index_eta]->SetLogx(1);
+    hMIPVsNch[index_eta]->GetXaxis()->SetLabelSize(0.08);
+
+    hMIPVsNch[index_eta]->GetXaxis()->SetRangeUser(10,2000);
+    hMIPVsNch[index_eta]->GetYaxis()->SetRangeUser(43.5,54.5);
+    hMIPVsNch[index_eta]->GetYaxis()->SetLabelSize(0.06);
+    hMIPVsNch[index_eta]->GetXaxis()->SetTitleSize(0);
+    hMIPVsNch[index_eta]->GetYaxis()->SetTitleSize(0);
+
+    if(index_eta<4){
+      hMIPVsNch[index_eta]->Draw("colz"); 
+      pMIPVsNch[index_eta]->Draw("samep"); 
+      latex->DrawLatex(0.5,0.1,LatexEta[index_eta]);
+    }
+    else{
+      hMIPVsNch[index_dec]->Draw("colz");
+      pMIPVsNch[index_dec]->Draw("samep"); 
+      latex->DrawLatex(0.5,0.2,LatexEta[index_dec]);
+      index_dec--;
+    }
+
+  }
+
+
+
+
+
+
+
+
+  cMIPVsNch->SaveAs("MIPvsNch.png");
+
+
+
+
+
+
+
+
+
+  //eta dependence, scaling
+  const Int_t nPadXS = 3;
+  const Int_t nPadYS = 1;
+  Float_t noMarginS    = 0.0;
+  Float_t topMarginS   = 0.01;
+  Float_t botMarginS   = 0.1;
+  Float_t leftMarginS  = 0.05;
+  Float_t rightMarginS = 0.01;
+  Float_t widthS       = (1-rightMarginS-leftMarginS)/nPadXS;
+  Float_t heightS      = (1-botMarginS-topMarginS)/nPadYS;
+  Float_t shiftS = 0.05;
+    
+  
+  TCanvas* cMIPEta = new TCanvas("cMIPEta", "MIP eta", 900, 500);
+  TPad* padMIPEta[nPadXS*nPadYS];
+  //Float_t shift = 0.1;
+  //Float_t shift = 0.05;
+  const char* yTitleMIPS = "d#it{E}/d#it{x}";
+  const char* xTitleMIPS = "#eta";
+  cMIPEta->cd();
+  
+  
+  //TLatexlatex = new TLatex();
+  latex->SetNDC();
+  latex->SetTextAlign(22);
+  latex->SetTextAngle(0);
+  latex->SetTextFont(42);
+  
+  latex->DrawLatex(0.53,0.04,xTitleMIPS);
+  latex->SetTextAngle(90);
+  latex->SetTextSize(0.055);
+  latex->DrawLatex(0.025,0.5,yTitleMIPS);
+  latex->SetTextSize(0.05);
+
+  
+  for (Int_t i = 1; i < 4; i++) {
+    
+    Int_t iX = (i-1)%nPadXS;
+    Int_t iY = (i-1)/nPadXS;
+    Float_t x1 = leftMarginS + iX*widthS;
+    if(iX==2)
+      x1-=0.01;
+    else if(iX==1)
+      x1+=0.01;
+    Float_t x2 = leftMarginS + (iX + 1)*widthS;
+     if(iX==0)
+       x2+=0.01;
+     else if(iX==1)
+       x2-=0.01;
+     Float_t y1 = 1 - topMarginS - (iY +1)*heightS;
+     Float_t y2 = 1 - topMarginS - iY*heightS;
+     padMIPEta[i-1] = new TPad(Form("padMIPEta%d",i),"", x1, y1, x2, y2, 0, 0, 0);
+     Float_t mBot = noMarginS; 
+     Float_t mTop = noMarginS; 
+     Float_t mLeft = noMarginS; 
+     Float_t mRight = noMarginS; 
+     if(iY==0)       mTop   = shiftS;
+     if(iY==nPadYS-1) mBot   = shiftS;
+     if(iX==0)       mLeft  = 0.08;
+     if(iX==nPadXS-1) mRight = 0.08;
+     padMIPEta[i-1]->SetBottomMargin(mBot);
+     padMIPEta[i-1]->SetTopMargin(mTop);
+     padMIPEta[i-1]->SetLeftMargin(mLeft);
+     padMIPEta[i-1]->SetRightMargin(mRight);
+     
+     cMIPEta->cd(); 
+     
+     padMIPEta[i-1]->Draw();
+     
+  }
+
+  latex->SetTextAngle(0);
+  latex->SetTextSize(0.085);
+  
+
+  TF1 *fetapos = new TF1("fetapos","pol3",0.0,0.8);
+  TF1 *fetaneg = new TF1("fetaneg","pol3",-0.8,0.0);
+
+  pMIPVsEta->Fit(fetaneg,"0","",-0.8,0.0);
+  pMIPVsEta->Fit(fetapos,"0","",0.0,0.8);
+  TLegend *leg1=new TLegend(0.15,0.5,0.4,0.7);
+
+  fetapos->SetLineColor(1);
+  fetapos->SetLineStyle(1);
+  fetapos->SetLineWidth(2);
+
+  fetaneg->SetLineColor(1);
+  fetaneg->SetLineStyle(1);
+  fetaneg->SetLineWidth(2);
+
+  TH1D *hframe20 = new TH1D("hframe20","hframe20",8,-0.8,0.8);
+
+  TH1D *hMIPPiPrim = new TH1D("hMIPPiPrim","",16,-0.8,0.8);
+  for(Int_t bin =1; bin <= hMIPPiPrim->GetNbinsX(); ++bin){
+    Double_t eta = pMIPVsEta->GetBinCenter(bin);
+    Double_t dedx = pMIPVsEta->GetBinContent(bin);
+    Double_t e_dedx = pMIPVsEta->GetBinError(bin);
+    if(eta<0){
+      dedx *= 50/fetaneg->Eval(eta);
+      e_dedx *= 50/fetaneg->Eval(eta);
+    }else{
+      dedx *= 50/fetapos->Eval(eta);
+      e_dedx *= 50/fetapos->Eval(eta);
+    }
+    hMIPPiPrim->SetBinContent(bin,dedx);
+    hMIPPiPrim->SetBinError(bin,e_dedx);
+
+  }
+
+  TH1D *hPlateauPiPrim = new TH1D("hPlateauPiPrim","",16,-0.8,0.8);
+  for(Int_t bin =1; bin <= hPlateauPiPrim->GetNbinsX(); ++bin){
+    Double_t eta = pPlateauVsEta->GetBinCenter(bin);
+    Double_t dedx = pPlateauVsEta->GetBinContent(bin);
+    Double_t e_dedx = pPlateauVsEta->GetBinError(bin);
+    if(eta<0){
+      dedx *= 50/fetaneg->Eval(eta);
+      e_dedx *= 50/fetaneg->Eval(eta);
+    }else{
+      dedx *= 50/fetapos->Eval(eta);
+      e_dedx *= 50/fetapos->Eval(eta);
+    }
+    hPlateauPiPrim->SetBinContent(bin,dedx);
+    hPlateauPiPrim->SetBinError(bin,e_dedx);
+
+  }
+
+
+  TH1D *h2GeVPiSec = new TH1D("h2GeVPiSec","",8,-0.8,0.8);
+  for(Int_t bin =1; bin <= h2GeVPiSec->GetNbinsX(); ++bin){
+    Double_t eta = hPion2GeV->GetBinCenter(bin);
+    Double_t dedx = hPion2GeV->GetBinContent(bin);
+    Double_t e_dedx = hPion2GeV->GetBinError(bin);
+    if(eta<0){
+      dedx *= 50/fetaneg->Eval(eta);
+      e_dedx *= 50/fetaneg->Eval(eta);
+    }else{
+      dedx *= 50/fetapos->Eval(eta);
+      e_dedx *= 50/fetapos->Eval(eta);
+    }
+    h2GeVPiSec->SetBinContent(bin,dedx);
+    h2GeVPiSec->SetBinError(bin,e_dedx);
+
+  }
+
+
+
+
+
+
+
+
+
+  
+  for(Int_t index_eta = 0; index_eta < 3; ++index_eta){
+    padMIPEta[index_eta]->cd();
+    padMIPEta[index_eta]->SetTickx(1);
+    padMIPEta[index_eta]->SetTicky(1);
+    padMIPEta[index_eta]->SetLogy(0);
+    padMIPEta[index_eta]->SetLogx(0);
+    if(index_eta==0){
+      hframe20->GetYaxis()->SetRangeUser(30,95);
+      hframe20->GetXaxis()->SetLabelSize(0.06);
+      hframe20->GetXaxis()->SetLabelOffset(0.004);
+      hframe20->GetYaxis()->SetLabelSize(0.06);
+      hframe20->GetXaxis()->SetTitleSize(0);
+      hframe20->GetYaxis()->SetTitleSize(0);
+      hframe20->Draw();
+      hMIPVsEta->Draw("colz same");
+      hPlateauVsEta->Draw("colz same");
+
+      latex->SetTextSize(0.05);
+      latex->DrawLatex(0.55,0.45,"MIP Pions, 0.4<#it{p}<0.6 GeV/c");
+      latex->DrawLatex(0.55,0.85,"Electrons, 0.4<#it{p}<0.6 GeV/c");
+
+      }
+
+    if(index_eta==1){
+      pMIPVsEtaV0s->GetXaxis()->SetLabelSize(0.06);
+      pMIPVsEtaV0s->GetYaxis()->SetLabelSize(0.06);
+      pMIPVsEtaV0s->GetXaxis()->SetLabelOffset(0.004);
+      pMIPVsEtaV0s->GetXaxis()->SetTitleSize(0);
+      pMIPVsEtaV0s->GetYaxis()->SetTitleSize(0);
+      pMIPVsEtaV0s->GetYaxis()->SetRangeUser(30,95);
+      pMIPVsEtaV0s->SetMarkerStyle(29);
+      pMIPVsEtaV0s->SetMarkerColor(4);
+      pMIPVsEtaV0s->SetLineColor(4);
+      pMIPVsEtaV0s->Draw();
+
+      pMIPVsEta->SetMarkerStyle(25);
+      pMIPVsEta->SetMarkerColor(1);
+      pMIPVsEta->SetLineColor(1);
+      pMIPVsEta->Draw("samep");
+      pPlateauVsEta->SetMarkerStyle(21);
+      pPlateauVsEta->SetMarkerColor(kBlue);
+      pPlateauVsEta->SetLineColor(4);
+      pPlateauVsEta->Draw("samep");
+
+      //fetapos->Draw("same");
+      //fetaneg->Draw("same");
+
+
+      pMIPVsEta->Draw("samep");
+      hPion2GeV->SetMarkerStyle(20);
+      hPion2GeV->SetMarkerColor(4);
+      hPion2GeV->SetLineColor(4);
+      hPion2GeV->Draw("samep");
+      pPlateauVsEta->Draw("samep");
+
+    
+      leg1->SetTextFont(42);
+      leg1->SetTextSize(0.045);
+      leg1->SetLineColor(kWhite);
+      leg1->SetLineStyle(3);
+      leg1->SetShadowColor(kWhite);
+      leg1->SetFillColor(kWhite);
+      leg1->SetFillStyle(0);
+      leg1->AddEntry(pMIPVsEta,"Primary pions at MIP","P");
+      leg1->AddEntry(pMIPVsEtaV0s,"Secondary pions at MIP","P");
+      leg1->AddEntry(hPion2GeV,"Secondary pions at #beta#gamma ~ 14","P");
+      leg1->AddEntry(pPlateauVsEta,"Electrons,  at #beta#gamma ~ 800","P");
+      leg1->Draw();
+
+
+
+    }
+
+
+
+    if(index_eta==2){
+      hMIPPiPrim->GetXaxis()->SetLabelSize(0.06);
+      hMIPPiPrim->GetYaxis()->SetLabelSize(0.06);
+      hMIPPiPrim->GetXaxis()->SetLabelOffset(0.004);
+      hMIPPiPrim->GetXaxis()->SetTitleSize(0);
+      hMIPPiPrim->GetYaxis()->SetTitleSize(0);
+      hMIPPiPrim->GetYaxis()->SetRangeUser(30,95);
+      hMIPPiPrim->SetMarkerStyle(29);
+      hMIPPiPrim->SetMarkerColor(4);
+      hMIPPiPrim->SetLineColor(4);
+      hMIPPiPrim->Draw();
+
+
+      hPlateauPiPrim->SetMarkerStyle(21);
+      hPlateauPiPrim->SetLineColor(4);
+      hPlateauPiPrim->SetMarkerColor(4);
+      hPlateauPiPrim->Draw("samep");
+      h2GeVPiSec->SetMarkerStyle(20);
+      h2GeVPiSec->SetLineColor(4);
+      h2GeVPiSec->SetMarkerColor(4);
+      h2GeVPiSec->Draw("samep");
+
+      latex->DrawLatex(0.45,0.85,"#LT d#it{E}/d#it{x} #GT #times 50 / #LT d#it{E}/d#it{x} #GT_{primary #pi at MIP}");
+
+    }
+
+
+
+
+
+
+
+
+
+
+
+  }
+
+
+  cMIPEta->SaveAs("MIPvsEta.png");
+
+
+  
+  
+}
+TH2D * AddTwoSameBinningTH2D(TH2D *hPos, TH2D *hNeg, const Char_t *nameHist){
+  
+  TH2D *hout = (TH2D *)hPos->Clone(nameHist);
+  hout->Reset();
+
+
+  for(Int_t binx = 1; binx <= hPos->GetNbinsX(); ++binx){
+    
+    for(Int_t biny = 1; biny <= hPos->GetNbinsY(); ++biny){
+      
+      Double_t yield1 =  hPos->GetBinContent(binx,biny);
+      Double_t yield2 =  hNeg->GetBinContent(binx,biny);
+      Double_t e_yield1 =  hPos->GetBinError(binx,biny);
+      Double_t e_yield2 =  hNeg->GetBinError(binx,biny);
+
+      Double_t yield = yield1 + yield2;
+      Double_t e_yield = TMath::Sqrt( e_yield1*e_yield1 + e_yield2*e_yield2 );
+
+      hout->SetBinContent(binx,biny,yield);
+      hout->SetBinError(binx,biny,e_yield);
+
+    }
+    
+  }
+
+  hout->SetEntries( hPos->GetEntries()+hNeg->GetEntries() );
+
+  return hout;
+
+}
diff --git a/PWGLF/QATasks/post/PostProcessQAMultistrange.C b/PWGLF/QATasks/post/PostProcessQAMultistrange.C
new file mode 100644 (file)
index 0000000..1678251
--- /dev/null
@@ -0,0 +1,488 @@
+//////////////////////////////////////////////////
+//
+//  This macro was written by Domenico Colella (domenico.colella@cern.ch).
+//  12 November 2013
+//
+//   ------ Arguments
+//   --  icasType          =  0) Xi- 1) Xi+ 2) Omega- 3) Omega+
+//   --  collidingsystem   =  0) PbPb  1) pp
+//   --  fileDir           =  "Input file directory"
+//   --  filein            =  "Input file name"
+//
+//   ------ QATask output content
+//   The output produced by the QATask is a CFContainer with 4 steps and 21 variables.
+//   The meaning of each variable within the container are listed here:
+//   --  0   = Max DCA Cascade Daughters                 pp: 2.0     PbPb: 0.3
+//   --  1   = Min DCA Bach To PV                        pp: 0.01    PbPb: 0.03
+//   --  2   = Min Cascade Cosine Of PA                  pp: 0.98    PbPb: 0.999
+//   --  3   = Min Cascade Radius Fid. Vol.              pp: 0.2     PbPb: 0.9
+//   --  4   = Window Invariant Mass Lambda              pp: 0.008   PbPb: 0.0008
+//   --  5   = Max DCA V0 Daughters                      pp: 1.5     PbPb: 1.0
+//   --  6   = Min V0 Cosine Of PA To PV                 pp: pT dep. PbPb: 0.98
+//   --  7   = Min V0 Radius Fid. Vol.                   pp: 0.2     PbPb: 0.9
+//   --  8   = Min DCA V0 To PV                          pp: 0.01    PbPb: 0.05
+//   --  9   = Min DCA Pos To PV                         pp: 0.05    PbPb: 0.1
+//   --  10  = Min DCA Neg To PV                         pp: 0.05    PbPb: 0.1
+//   --  11  = Invariant Mass distribution for Xi
+//   --  12  = Invariant Mass distribution for Omega
+//   --  13  = Transverse Momentum distribution
+//   --  14  = Rapidity distribution for Xi
+//   --  15  = Rapidity distribution for Omega
+//   --  16  = Proper length distribution for the cascade
+//   --  17  = Proper length distribution for the V0
+//   --  18  = Min V0 Cosine Of PA To Xi Vertex         pp: pT dep. PbPb: pT dep.
+//   --  19  = Centrality
+//   --  20  = ESD track multiplicity
+//   The last two variables are empty in the case proton-proton collisions.
+//
+//   ------ Present Macro Check
+//   With this macro are produced the plots of the distributions for the topological
+//   variables used during the reconstruction of the cascades and defined in the two
+//   classes AliCascadeVertexer.cxx and AliV0vertexer.cxx contained in /STEER/ESD/ 
+//   folder in Aliroot.
+//
+//   -- First Canvas: DCA cascade daughters, Bachelor IP to PV, Cascade cosine of PA
+//                    Cascade radius of fiducial volume, Invariant mass Lambda,
+//                    DCA V0 daughters.
+//   -- Second Canvas: V0 cosine of PA to PV, Min V0 Radius Fid. Vol., Min DCA V0 To PV
+//                     Min DCA Pos To PV, Min DCA Neg To PV, V0 cosine of PA to XiV
+//
+//   In this plots, in correspondence to the min/max cut value adopted in the reconstruction
+//   a line is drawn, to check if there is same problem in the cuts definition.
+//
+//   Also, other specific distribution for the selected cascades are ploted as: the
+//   invariant mass distribution, the transverse momentum distribution, the rapidity
+//   distribution, proper length distribution for the cascade and the v0.
+//
+//   -- Third Canvas: InvMass, Transverse momentum, Cascade proper length, V0 proper length
+//
+//   In the end, only for thr PbPb collision the centrality distribution and the
+//   track multiplicity distribution are sored.
+//
+//   -- Fourth Canvas: Centrality, track multiplicity
+//
+//
+//////////////////////////////////////////////////////
+
+
+
+
+class AliCFContainer;
+
+void PostProcessQAMultistrange(Int_t   icasType        = 0,                             // 0) Xi- 1) Xi+ 2) Omega- 3) Omega+
+                               Int_t   collidingsystem = 0,                             // 0) PbPb  1) pp
+                               Char_t *fileDir         = "./",                          // Input file directory
+                               Char_t *filein          = "AnalysisResults_AOD.root"     // Input file name
+                              ) {
+
+
+      /////////////
+      gStyle->SetOptStat(1110);
+      gStyle->SetOptStat(kFALSE);
+      gStyle->SetOptTitle(kFALSE);
+      gStyle->SetFrameLineWidth(2.5);
+      gStyle->SetCanvasColor(0);
+      gStyle->SetPadColor(0);
+      gStyle->SetHistLineWidth(2.5);
+      gStyle->SetLabelSize(0.05, "x");
+      gStyle->SetLabelSize(0.05, "y");
+      gStyle->SetTitleSize(0.05, "x");
+      gStyle->SetTitleSize(0.05, "y");
+      gStyle->SetTitleOffset(1.1, "x");
+      gStyle->SetPadBottomMargin(0.14);
+      gSystem->Load("libANALYSIS.so");
+      gSystem->Load("libANALYSISalice.so");
+      gSystem->Load("libCORRFW.so");
+
+
+     TFile *f1 = new TFile(Form("%s/%s",fileDir,filein));
+     AliCFContainer *cf = (AliCFContainer*) (f1->Get("PWGLFStrangeness.outputCheckCascade/fCFContCascadeCuts"));
+   
+
+     //DEEFINE TEXT
+     TLatex* t1 = new TLatex(0.6,0.55,"#color[3]{OK!!}");
+     t1->SetTextSize(0.1);
+     t1->SetNDC();
+     TLatex* t2 = new TLatex(0.6,0.55,"#color[2]{NOT OK!!}");
+     t2->SetTextSize(0.1);
+     t2->SetNDC();
+     t2->SetTextColor(2);
+     TLatex* tcasc;
+     if      (icasType == 0) tcasc = new TLatex(0.8,0.7,"#color[1]{#Xi^{-}}");
+     else if (icasType == 1) tcasc = new TLatex(0.8,0.7,"#color[1]{#Xi^{+}}");
+     else if (icasType == 2) tcasc = new TLatex(0.8,0.7,"#color[1]{#Omega^{-}}");
+     else if (icasType == 3) tcasc = new TLatex(0.8,0.7,"#color[1]{#Omega^{+}}");
+     tcasc->SetTextSize(0.2);
+     tcasc->SetNDC();
+     tcasc->SetTextColor(2);
+     TLatex* tpdgmass;
+     if      (icasType == 0) tpdgmass = new TLatex(0.55,0.7,"#color[1]{PDG mass: 1.321 GeV/c^{2}}");
+     else if (icasType == 1) tpdgmass = new TLatex(0.55,0.7,"#color[1]{PDG mass: 1.321 GeV/c^{2}}");
+     else if (icasType == 2) tpdgmass = new TLatex(0.55,0.7,"#color[1]{PDG mass: 1.672 GeV/c^{2}}");
+     else if (icasType == 3) tpdgmass = new TLatex(0.55,0.7,"#color[1]{PDG mass: 1.672 GeV/c^{2}}");
+     tpdgmass->SetTextSize(0.07);
+     tpdgmass->SetNDC();
+     tpdgmass->SetTextColor(2);
+     //DEFINE 1st CANVAS AND DRAW PLOTS
+     TCanvas *c1 = new TCanvas("c1","",1200,800);
+     c1->Divide(2,3); 
+       //Pad 1: DCA cascade daughters
+       c1->cd(1);
+       gPad->SetLogy();
+       TH1D *hvar0 = cf->ShowProjection(0,icasType);
+       hvar0->Draw("histo");
+       Double_t x0;
+       if      (collidingsystem == 0) x0 = 0.3;
+       else if (collidingsystem == 1) x0 = 2.0;
+       TLine *line0 = new TLine(x0,0.,x0,hvar0->GetBinContent(hvar0->GetMaximumBin()));
+       line0->SetLineColor(kRed);
+       line0->SetLineStyle(9);
+       line0->SetLineWidth(2.0);
+       line0->Draw("same");
+          Bool_t check_0 = checkOverTheLimit(hvar0, x0);
+          if (check_0) { cout<<"The cut is OK!!"<<endl; t1->Draw(); }
+          else         { cout<<"The cut is NOT OK!!"<<endl; t2->Draw(); }
+       tcasc->Draw();
+       //Pad 2: Bachelor IP to PV
+       c1->cd(2);
+       gPad->SetLogy();
+       TH1D *hvar1 = cf->ShowProjection(1,icasType);
+       hvar1->GetXaxis()->SetRangeUser(0.,0.24);
+       hvar1->Draw("histo");
+       Double_t x1;
+       if      (collidingsystem == 0) x1 = 0.03;
+       else if (collidingsystem == 1) x1 = 0.01;
+       TLine *line1 = new TLine(x1,0.,x1,hvar1->GetBinContent(hvar1->GetMaximumBin()));
+       line1->SetLineColor(kRed);
+       line1->SetLineStyle(9);
+       line1->SetLineWidth(2.0);
+       line1->Draw("same");
+          Bool_t check_1 = checkUnderTheLimit(hvar1, x1);
+          if (check_1) { cout<<"The cut is OK!!"<<endl; t1->Draw(); }
+          else         { cout<<"The cut is NOT OK!!"<<endl; t2->Draw(); }
+       //Pad 3: Cascade cosine of Pointing Angle
+       c1->cd(3);
+       gPad->SetLogy();
+       TH1D *hvar2 = cf->ShowProjection(2,icasType);
+       Double_t max2 = hvar2->GetBinContent(hvar2->GetMaximumBin());
+       hvar2->GetYaxis()->SetRangeUser(0.01,max2*1.5);
+       hvar2->Draw("histo");
+       Double_t x2;
+       if      (collidingsystem == 0) x2 = 0.999;
+       else if (collidingsystem == 1) x2 = 0.98;
+       TLine *line2 = new TLine(x2,0.,x2,hvar2->GetBinContent(hvar2->GetMaximumBin()));
+       line2->SetLineColor(kRed);
+       line2->SetLineStyle(9);
+       line2->SetLineWidth(2.0);
+       line2->Draw("same");
+       line1->Draw("same");
+          Bool_t check_2 = checkUnderTheLimit(hvar2, x2);
+          if (check_2) { cout<<"The cut is OK!!"<<endl; t1->Draw(); }
+          else         { cout<<"The cut is NOT OK!!"<<endl; t2->Draw(); }
+       //Pad 4: Cascade radius of fiducial volume
+       c1->cd(4);
+       gPad->SetLogy();
+       TH1D *hvar3 = cf->ShowProjection(3,icasType);
+       hvar3->GetXaxis()->SetRangeUser(0.,3.8);
+       hvar3->Draw("histo");
+       Double_t x3;
+       if      (collidingsystem == 0) x3 = 0.9;
+       else if (collidingsystem == 1) x3 = 0.2;
+       TLine *line3 = new TLine(x3,0.,x3,hvar3->GetBinContent(hvar3->GetMaximumBin()));
+       line3->SetLineColor(kRed);
+       line3->SetLineStyle(9);
+       line3->SetLineWidth(2.0);
+       line3->Draw("same");
+          Bool_t check_3 = checkUnderTheLimit(hvar3, x3);
+          if (check_3) { cout<<"The cut is OK!!"<<endl; t1->Draw(); }
+          else         { cout<<"The cut is NOT OK!!"<<endl; t2->Draw(); }
+       //Pad 5: Invariant mass Lambda
+       c1->cd(5);
+       TH1D *hvar4 = cf->ShowProjection(4,icasType);
+       hvar4->Draw("histo");
+       Double_t x41 = 1.116 + 0.008;
+       TLine *line41 = new TLine(x41,0.,x41,hvar4->GetBinContent(hvar4->GetMaximumBin()));
+       line41->SetLineColor(kRed);
+       line41->SetLineStyle(9);
+       line41->SetLineWidth(2.0);
+       line41->Draw("same");
+       Double_t x42 = 1.115 - 0.008;
+       TLine *line42 = new TLine(x42,0.,x42,hvar4->GetBinContent(hvar4->GetMaximumBin()));
+       line42->SetLineColor(kRed);
+       line42->SetLineStyle(9);
+       line42->SetLineWidth(2.0);
+       line42->Draw("same");
+          Bool_t check_4_1 = checkUnderTheLimit(hvar3, x3);
+          Bool_t check_4_2 = checkOverTheLimit(hvar0, x0);
+          if (check_4_1 && check_4_2) { cout<<"The cut is OK!!"<<endl; t1->Draw(); }
+          else                        { cout<<"The cut is NOT OK!!"<<endl; t2->Draw(); }
+       //Pad 6: DCA V0 daughters
+       c1->cd(6);
+       gPad->SetLogy();
+       TH1D *hvar5 = cf->ShowProjection(5,icasType);
+       hvar5->Draw("histo");
+       Double_t x5;
+       if      (collidingsystem == 0) x5 = 1.0;
+       else if (collidingsystem == 1) x5 = 1.5;
+       TLine *line5 = new TLine(x5,0.,x5,hvar5->GetBinContent(hvar5->GetMaximumBin()));
+       line5->SetLineColor(kRed);
+       line5->SetLineStyle(9);
+       line5->SetLineWidth(2.0);
+       line5->Draw("same");
+          Bool_t check_5 = checkOverTheLimit(hvar5, x5);
+          if (check_5) { cout<<"The cut is OK!!"<<endl; t1->Draw(); }
+          else         { cout<<"The cut is NOT OK!!"<<endl; t2->Draw(); }
+     c1->SaveAs("fig_lf_Multistrange.pdf(");
+    
+
+     //DEFINE 2st CANVAS AND DRAW PLOTS
+     TCanvas *c2 = new TCanvas("c2","",1200,800);
+     c2->Divide(2,3);
+       //Pad 1: V0 cosine of Pointing Angle to PV
+       c2->cd(1);
+       gPad->SetLogy();
+       TH1D *hvar6 = cf->ShowProjection(6,icasType);
+       Double_t max6 = hvar6->GetBinContent(hvar6->GetMaximumBin());
+       hvar6->GetYaxis()->SetRangeUser(0.01,max6*1.5);
+       hvar6->Draw("histo");
+       //Pad 2: Min V0 Radius Fid. Vol.  
+       c2->cd(2);
+       gPad->SetLogy();
+       TH1D *hvar7 = cf->ShowProjection(7,icasType);
+       hvar7->GetXaxis()->SetRangeUser(0.,3.0);
+       hvar7->Draw("histo");
+       Double_t x7;
+       if      (collidingsystem == 0) x7 = 0.9;
+       else if (collidingsystem == 1) x7 = 0.2;
+       TLine *line7 = new TLine(x7,0.,x7,hvar7->GetBinContent(hvar7->GetMaximumBin()));
+       line7->SetLineColor(kRed);
+       line7->SetLineStyle(9);
+       line7->SetLineWidth(2.0);
+       line7->Draw("same");
+          Bool_t check_7 = checkUnderTheLimit(hvar7, x7);
+          if (check_7) { cout<<"The cut is OK!!"<<endl; t1->Draw(); }
+          else         { cout<<"The cut is NOT OK!!"<<endl; t2->Draw(); }
+       //Pad3: Min DCA V0 To PV
+       c2->cd(3);
+       gPad->SetLogy();
+       TH1D *hvar8 = cf->ShowProjection(8,icasType);
+       hvar8->GetXaxis()->SetRangeUser(0.,0.3);
+       hvar8->Draw("histo");
+       Double_t x8;
+       if      (collidingsystem == 0) x8 = 0.05;
+       else if (collidingsystem == 1) x8 = 0.01;
+       TLine *line8 = new TLine(x8,0.,x8,hvar8->GetBinContent(hvar8->GetMaximumBin()));
+       line8->SetLineColor(kRed);
+       line8->SetLineStyle(9);
+       line8->SetLineWidth(2.0);
+       line8->Draw("same");
+          Bool_t check_8 = checkUnderTheLimit(hvar8, x8);
+          if (check_8) { cout<<"The cut is OK!!"<<endl; t1->Draw(); }
+          else         { cout<<"The cut is NOT OK!!"<<endl; t2->Draw(); }
+       //Pad 4: Min DCA Pos To PV
+       c2->cd(4);
+       gPad->SetLogy();
+       TH1D *hvar9 = cf->ShowProjection(9,icasType);
+       hvar9->GetXaxis()->SetRangeUser(0.,0.2);
+       hvar9->Draw("histo");
+       Double_t x9;
+       if      (collidingsystem == 0) x9 = 0.1;
+       else if (collidingsystem == 1) x9 = 0.05;
+       TLine *line9 = new TLine(x9,0.,x9,hvar9->GetBinContent(hvar9->GetMaximumBin()));
+       line9->SetLineColor(kRed);
+       line9->SetLineStyle(9);
+       line9->SetLineWidth(2.0);
+       line9->Draw("same");
+          Bool_t check_9 = checkUnderTheLimit(hvar9, x9);
+          if (check_9) { cout<<"The cut is OK!!"<<endl; t1->Draw(); }
+          else         { cout<<"The cut is NOT OK!!"<<endl; t2->Draw(); }
+       //Pad 5: Min DCA Neg To PV
+       c2->cd(5);
+       gPad->SetLogy();
+       TH1D *hvar10 = cf->ShowProjection(10,icasType);
+       hvar10->GetXaxis()->SetRangeUser(0.,0.2);
+       hvar10->Draw("histo");
+       Double_t x10;
+       if      (collidingsystem == 0) x10 = 0.1;
+       else if (collidingsystem == 1) x10 = 0.05;
+       TLine *line10 = new TLine(x10,0.,x10,hvar10->GetBinContent(hvar10->GetMaximumBin()));
+       line10->SetLineColor(kRed);
+       line10->SetLineStyle(9);
+       line10->SetLineWidth(2.0);
+       line10->Draw("same");
+          Bool_t check_10 = checkUnderTheLimit(hvar10, x10);
+          if (check_10) { cout<<"The cut is OK!!"<<endl; t1->Draw(); }
+          else         { cout<<"The cut is NOT OK!!"<<endl; t2->Draw(); }
+       //Pad 6: V0 cosine of Pointing Angle to XiV
+       c2->cd(6);
+       gPad->SetLogy();
+       TH1D *hvar20 = cf->ShowProjection(18,icasType);
+       Double_t max20 = hvar20->GetBinContent(hvar20->GetMaximumBin());
+       hvar20->GetYaxis()->SetRangeUser(0.01,max20*1.5);
+       hvar20->Draw("histo");
+     c2->SaveAs("fig_lf_Multistrange.pdf");
+
+     //DEFINE 3st CANVAS AND DRAW PLOTS
+     TCanvas *c3 = new TCanvas("c3","",1200,800);
+     c3->Divide(2,3);
+       //Pad 1: InvMass
+       c3->cd(1);
+       TH1D *hvar12 = cf->ShowProjection(11+icasType/2,icasType);
+       hvar12->Draw("histo");
+       tpdgmass->Draw(); 
+       TLine *linemass;
+       if      (icasType == 0) linemass = new TLine(1.32171,0.,1.32171,0.5*hvar12->GetBinContent(hvar12->GetMaximumBin()));
+       else if (icasType == 1) linemass = new TLine(1.32171,0.,1.32171,0.5*hvar12->GetBinContent(hvar12->GetMaximumBin()));
+       else if (icasType == 2) linemass = new TLine(1.67245,0.,1.67245,0.5*hvar12->GetBinContent(hvar12->GetMaximumBin()));
+       else if (icasType == 3) linemass = new TLine(1.67245,0.,1.67245,0.5*hvar12->GetBinContent(hvar12->GetMaximumBin()));
+       linemass->SetLineColor(kRed);
+       linemass->SetLineStyle(1);
+       linemass->SetLineWidth(2.0);
+       linemass->Draw("same");
+       //Pad 2: Transverse momentum
+       c3->cd(2);
+       TH1D *hvar13 = cf->ShowProjection(13,icasType);
+       hvar13->Draw("histo");
+       //Pad 3: Y
+       c3->cd(3);
+       TH1D *hvar14 = cf->ShowProjection(14+icasType/2,icasType);
+       hvar14->Draw("histo");
+       //Pad 4: Cascade proper length
+       c3->cd(4);
+       TH1D *hvar18;
+       hvar18 = cf->ShowProjection(16,icasType);
+       hvar18->GetXaxis()->SetRangeUser(0.,90.);
+       hvar18->Draw("histo");
+       //Pad 5: V0 proper length 
+       c3->cd(5);
+       TH1D *hvar19;
+       hvar19 = cf->ShowProjection(17,icasType);
+       hvar19->GetXaxis()->SetRangeUser(0.,90.);
+       hvar19->Draw("histo");
+      //Pad 6
+      // empty 
+     if      (collidingsystem == 1) c3->SaveAs("fig_lf_Multistrange.pdf)");
+     else if (collidingsystem == 0) c3->SaveAs("fig_lf_Multistrange.pdf");
+
+    
+     //DEFINE 4st CANVAS AND DRAW PLOTS
+    TCanvas *c4 = new TCanvas("c4","",600,400);
+    c4->Divide(2,1);
+      //Pad1: invariant mass fit
+      c4->cd(1);
+      TH1D *hvar18 = cf->ShowProjection(11+icasType/2,icasType);
+      hvar18->Draw("histo");
+       // - SOME PARAMETER VALUE
+       Bool_t kfitgauss = kFALSE;
+       Bool_t kfitleft  = kFALSE;
+       Bool_t kfitright = kFALSE;
+       Int_t  ptbinNarrowY = 0;
+       if (icasType < 2) ptbinNarrowY = 10;   // 6;
+       else              ptbinNarrowY =  3;   // 2;
+       // - SOME DEFINITIONS
+       Float_t lowlimmass;
+       Float_t uplimmass;
+       Float_t lowgausslim;
+       Float_t upgausslim;
+       if (icasType==0||icasType==1) {
+           lowlimmass=1.30;
+           uplimmass=1.34;
+           lowgausslim=1.312;
+           upgausslim=1.332;
+       } else {
+           lowlimmass=1.645;
+           uplimmass=1.70;
+           lowgausslim=1.668;
+           upgausslim=1.678;
+       }
+       TF1*  fitinvmass = new TF1("fitinvmass","gaus(0)+pol2(3)",lowlimmass,uplimmass);
+       fitinvmass->SetParName(0, "cnstntG");
+       fitinvmass->SetParName(1, "meanG");
+       fitinvmass->SetParName(2, "sigmaG");
+       fitinvmass->SetParLimits(0,0.,500000.);
+       if (icasType==0||icasType==1) {
+           fitinvmass->SetParameter(1, 1.32171);
+           fitinvmass->SetParLimits(1, 1.31,1.33);
+           fitinvmass->SetParLimits(2,0.001,0.005);
+       } else {
+           fitinvmass->SetParameter(1, 1.67245);
+           fitinvmass->SetParLimits(1, 1.664,1.68);
+           fitinvmass->SetParLimits(2,0.0008,0.006);
+       }
+       hvar18->Fit("fitinvmass","rimeN");
+       fitinvmass->SetLineColor(kRed);
+       fitinvmass->Draw("same");
+       Float_t meanGauss   = fitinvmass->GetParameter(1);
+       Float_t sigmaGauss  = fitinvmass->GetParameter(2);
+       cout<<"Mean: "<<meanGauss<<endl;
+       cout<<"Sigma: "<<sigmaGauss<<endl;
+     //Pad2: Text
+      c4->cd(2);
+       Float_t refwidth = 0.002;
+      TPaveText *pave1 = new TPaveText(0.05,0.3,0.95,0.5);
+      pave1->SetFillColor(0);
+      pave1->SetTextSize(0.04);
+      pave1->SetTextAlign(12);
+      if (icasType < 2) pave1->AddText("PDG mass: 1.32171 GeV/c^{2}");
+      else              pave1->AddText("PDG mass: 1.67245 GeV/c^{2}");
+      pave1->AddText(Form("#color[1]{Mass form Fit: %.5f #pm %.5f GeV/c^{2}}",meanGauss,sigmaGauss));
+      if (sigmaGauss > refwidth - 0.0003 && sigmaGauss < refwidth + 0.0003) pave1->AddText("#color[3]{OK!! The width is compatible with standard.}");
+      else                                                                  pave1->AddText("#color[2]{NOT OK!! Problem.}");
+      pave1->Draw();
+      cout<<"   "<<refwidth - 0.0003<<"<"<<sigmaGauss<<"<"<<refwidth + 0.0003<<endl;
+    
+     //DEFINE 5st CANVAS AND DRAW PLOTS
+     if (collidingsystem == 0) {
+         TCanvas *c5 = new TCanvas("c5","",1200,270);
+         c5->Divide(2,1);
+            //Pad 1: centrality
+            c5->cd(1);
+            TH1D *hvar16 = cf->ShowProjection(19,icasType);
+            hvar16->Draw("histo");
+            //Pad 2: track multiplicity
+            c5->cd(2);
+            TH1D *hvar17 = cf->ShowProjection(20,icasType);
+            hvar17->Draw("histo");
+         c5->SaveAs("fig_lf_Multistrange.pdf)");
+     }
+    
+
+}
+
+
+
+
+Bool_t checkUnderTheLimit(TH1D *lHist, Double_t limit) {
+
+      Int_t binlimit = lHist->FindBin(limit);
+      Bool_t checkOk = kTRUE;
+      for (Int_t i = 1; i < binlimit; i++) {
+           Int_t content = 0;
+           content = lHist->GetBinContent(i);
+           if (content != 0) checkOk = kFALSE;
+           //cout<<"Content bin "<<i<<": "<<content<<endl;
+      }
+      return checkOk;
+
+}
+
+
+Bool_t checkOverTheLimit(TH1D *lHist, Double_t limit) {
+
+      Int_t binlimit = lHist->FindBin(limit);
+      Int_t lastbin  = lHist->GetNbinsX();
+      Bool_t checkOk = kTRUE;
+      for (Int_t i = binlimit; i < lastbin+1; i++) {
+           Int_t content = 0;
+           content = lHist->GetBinContent(i);
+           if (content != 0) checkOk = kFALSE;
+           //cout<<"Content bin "<<i<<": "<<content<<endl;
+      }
+      return checkOk;
+
+}
+
+
+
diff --git a/PWGLF/QATasks/post/PostProcessQAPhi.C b/PWGLF/QATasks/post/PostProcessQAPhi.C
new file mode 100644 (file)
index 0000000..c018965
--- /dev/null
@@ -0,0 +1,748 @@
+void PostProcessQAPhi(char* system="pp276",char* name_fin="AnalysisResults.root",char* name_fout="fig_lf_Phi",char* name_list="RsnHistMini_Phi_PhiNsigma_KTPCnsig30",char* name_hist_base="RsnMini_phi.PhiNsigma_KTPCnsig30"){
+  //This macro was written by Anders Knospe (anders.knospe@cern.ch).
+  //5 November 2013
+
+  //FOR DETAILED INSTRUCTIONS, SEE THE END OF THIS DOCUMENT.
+
+  //Arguments:
+  //system: string giving the collision system ("pp276","pp7", or "PbPb276"), used to figure out what the expected yield should be.  If the system you want is unavailable, you will need to calculate your own expected yield.
+  //name_fin: name of input file
+  //name_fout: base name of output files (without suffix)
+  //name_list: name of the list in fin that contains histograms, may be different depending on the cuts applied to the kaon daughters
+  //name_hist_base: base name of the THnSparse histograms, may be different depending on the cuts applied to the kaon daughters
+
+  //This Macro does the following:
+  //1.) plots the standard histograms showing the number of events that pass selection of criteria and the number of events in each multiplicity or centrality bin
+  //2.) plots invariant-mass histograms (unlike-charge, like-charge, background-subtracted)
+  //3.) fits the background-subtracted invariant-mass distribution.  In case there is a problem with the "default" fit, it does multiple fits using different residual backgrounds and fit regions.  All fits are plotted and all fit results are saved.
+  //4.) prints out relevant information: pT range, phi yield per event, phi mass, and phi width (taken from the default fit)
+
+  gStyle->SetOptStat(0);
+
+  //----- get input histograms -----
+
+  TFile* fin=TFile::Open(name_fin);//open input file
+  if(!fin) return;
+
+  TList* l=(TList*) fin->Get(name_list);
+  if(!l){cerr<<"Error in QAphi(): missing input list "<<name_list<<" in file "<<name_fin<<".  Stopping."<<endl; return;}
+
+  TH1D *hevent,*hmult,*hu,*hm,*hp,*hl,*hs,*hs_plot;
+  THnSparse *su,*sm,*sp;
+
+  hevent=(TH1D*) l->FindObject("hEventStat");//standard histogram in resonance package, shows number of events after various selections
+  if(!hevent){cerr<<"Error in QAphi(): missing input histogram hEventStat in file "<<name_fin<<".  Stopping."<<endl; return;}
+
+  hmult=(TH1D*) l->FindObject("hAEventsVsMulti");//standard histogram in resonance package, shows number of events in each multiplicity (for pp) or centrality (for Pb-Pb) bin
+  if(!hmult){cerr<<"Error in QAphi(): missing input histogram hAEventsVsMulti in file "<<name_fin<<".  Stopping."<<endl; return;}
+  double nevt=hmult->Integral();
+
+  su=(THnSparse*) l->FindObject(Form("%s_Unlike",name_hist_base));//invariant-mass histogram, unlike-charge
+  if(!su){cerr<<"Error in QAphi(): missing input histogram "<<name_hist_base<<"_Unlike in file "<<name_fin<<".  Stopping."<<endl; return;}
+
+  sm=(THnSparse*) l->FindObject(Form("%s_LikeMM",name_hist_base));//invariant-mass histogram, like-charge (K-K-)
+  if(!sm){cerr<<"Error in QAphi(): missing input histogram "<<name_hist_base<<"_LikeMM in file "<<name_fin<<".  Stopping."<<endl; return;}
+
+  sp=(THnSparse*) l->FindObject(Form("%s_LikePP",name_hist_base));//invariant-mass histogram, like-charge (K+K+)
+  if(!sp){cerr<<"Error in QAphi(): missing input histogram "<<name_hist_base<<"_LikePP in file "<<name_fin<<".  Stopping."<<endl; return;}
+
+  TFile* fout=new TFile(Form("%s.root",name_fout),"RECREATE","HistoFile");//open output file
+
+  bool ptOK=SetPtRange(su);//Was the requested pT range set correctly?
+  SetPtRange(sm);
+  SetPtRange(sp);
+
+  double dy;
+  bool yOK=SetRapidityRange(su,dy);//Is the expected rapidity range (|y|<0.5) being used?  Fill value of dy.
+  SetRapidityRange(sm,dy);
+  SetRapidityRange(sp,dy);
+
+  //----- plot the event and multiplicity/centrality histograms -----
+
+  TCanvas* c=new TCanvas("c","",10,10,1500,500);
+  c->SetFillColor(0);
+
+  TPad* p1=new TPad("p1","",0.,0.,1./3,1.);
+  p1->SetFillColor(0);
+
+  TPad* p2=new TPad("p2","",1./3,0.,2./3,1.);
+  p2->SetFillColor(0);
+
+  TPad* p3=new TPad("p3","",2./3,0.,1.,1.);
+  p3->SetFillColor(0);
+  p3->SetLogy();
+
+  p1->cd();
+  hevent->SetLineColor(1);
+  hevent->Draw();
+
+  p2->cd();
+  hmult->SetLineColor(1);
+  int j,k,xmax=0;
+  for(j=1;j<=hmult->GetNbinsX();j++) if(hmult->GetBinContent(j)>0.) xmax=j;
+  xmax=(int) (1.1*xmax);
+  if(xmax>hmult->GetNbinsX()) xmax=hmult->GetNbinsX();
+  hmult->GetXaxis()->SetRange(1,xmax);
+  hmult->Draw();
+
+  p3->cd();
+  hmult->Draw();//same as p2, but log scale on y axis
+
+  c->cd();
+  p1->Draw();
+  p2->Draw();
+  p3->Draw();
+
+  c->SaveAs(Form("%s.eps(",name_fout));
+
+  delete p1;
+  delete p2;
+  delete p3;
+  delete c;
+
+  fout->cd();
+  hevent->Write();
+  hmult->Write();
+
+  //----- get invariant-mass histograms -----
+
+  //Project the THnSparse histograms onto their x axes.
+  hu=(TH1D*) su->Projection(0,"e");
+  hu->SetName("mass_unlike");
+  hm=(TH1D*) sm->Projection(0,"e");
+  hm->SetName("mass_likeMM");
+  hp=(TH1D*) sp->Projection(0,"e");
+  hp->SetName("mass_likePP");
+
+  double A,UA,B,UB;
+
+  hl=(TH1D*) hm->Clone("mass_like");//total like-charge histogram
+  for(j=1;j<=hl->GetNbinsX();j++){
+    A=hm->GetBinContent(j);
+    B=hp->GetBinContent(j);
+    hl->SetBinContent(j,2*sqrt(A*B));//note: 2sqrt(n{K-K-}*n{K+K+})
+    hl->SetBinError(j,sqrt(A+B));
+  }
+  hl->SetTitle("Like-Charge Combinatorial Background");
+
+  hs=(TH1D*) hu->Clone("mass_signal");
+  for(j=1;j<=hs->GetNbinsX();j++){
+    A=hu->GetBinContent(j); UA=hu->GetBinError(j);
+    B=hl->GetBinContent(j); UB=hl->GetBinError(j);
+    hs->SetBinContent(j,A-B);//subtract the combinatorial background
+    hs->SetBinError(j,sqrt(UA*UA+UB*UB));
+  }
+
+  double Imin=1.01,Imax=1.03,dx=hs->GetXaxis()->GetBinWidth(1);
+
+  hb=(TH1D*) hs->Clone("mass_background");//temporary histogram with the peak removed; used to fit the residual background
+  for(j=1;j<=hb->GetNbinsX();j++){
+    A=hb->GetXaxis()->GetBinCenter(j);
+    if(A>Imin && A<Imax){hb->SetBinContent(j,0); hb->SetBinError(j,0);}
+  }
+
+  hs_plot=(TH1D*) hs->Clone("mass_signal_plot");//temporary histogram for plotting
+
+  //----- plot the invariant-mass histograms -----
+
+  c=new TCanvas("c","",10,10,1000,500);
+  c->SetFillColor(0);
+
+  p1=new TPad("p11","",0.,0.,0.5,1.);
+  p1->SetFillColor(0);
+  p1->SetRightMargin(0.05);
+
+  p2=new TPad("p21","",0.5,0.,1.,1.);
+  p2->SetFillColor(0);
+  p2->SetRightMargin(0.05);
+
+  p1->cd();
+  hu->SetLineColor(1); hu->SetMarkerColor(1);
+  hl->SetLineColor(2); hl->SetMarkerColor(2);
+  hu->SetTitle("Invariant Mass (KK)     ");
+  hu->SetXTitle("Invariant Mass (KK) [GeV/#it{c}^{2}]");
+  hu->Draw();
+  hl->Draw("same");
+
+  TLatex* t1=new TLatex(0.7,0.96,"Unlike-Charge");
+  t1->SetTextSize(0.04);
+  t1->SetNDC();
+  t1->Draw();
+
+  TLatex* t2=new TLatex(0.7,0.91,"Like-Charge");
+  t2->SetTextSize(0.04);
+  t2->SetNDC();
+  t2->SetTextColor(2);
+  t2->Draw();
+
+  p2->cd();
+  hs->SetLineColor(1); hs->SetMarkerColor(1);
+  hs->SetTitle("Combinatorial Background Subtracted");
+  hs->SetXTitle("Invariant Mass (KK) [GeV/#it{c}^{2}]");
+  hs->Draw();
+
+  c->cd();
+  p1->Draw();
+  p2->Draw();
+
+  c->SaveAs(Form("%s.eps",name_fout));
+
+  delete p1;
+  delete p2;
+  delete c;
+
+  fout->cd();
+  hu->Write();
+  hl->Write();
+  hs->Write();
+
+  //----- fit the signal histogram -----
+
+  c=new TCanvas("c","",10,10,500,500);
+  c->SetFillColor(0);
+
+  //fit peak and extract yield, mass, and width
+  //The peak is fit using a polynomial for the residual background, plus a Voigtian peak.  A Voigtian peak is the convolution of a (non-relativistic) Breit-Wigner peak and a Gaussian (to describe detector effects).  The resolution (Gaussian sigma) of the Voigtian peak will be fixed.
+  //Ideally, we would only need to fit once.  However, the fits are not always stable so we try a variety of fits and record all results.  A similar procedure was followed in the analysis of phi mesons in Pb+Pb collisions (2010 data).
+
+  int jBF;//index controlling the order of the residual background polynomial
+  int nBF=3;//number of different orders used for the residual background polynomial; currently 3 (linear, quadratic, cubic)
+  int jFR;//index controlling the fit region
+  int nFR=4;//number of different fit regions
+  double urb;
+
+  //This TH2D is just an array of numbers to hold the different fit results and associated quantities.  The x-axis gives the name of the fit (order of residual background and number of the fit region), the y-axis gives the name of the quantity stored.
+  TH2D* r=new TH2D("fit_results","",nBF*nFR,0,nBF*nFR, 16,0,16);
+  r->GetYaxis()->SetBinLabel(1,"peak integral");
+  r->GetYaxis()->SetBinLabel(2,"mass");
+  r->GetYaxis()->SetBinLabel(3,"width");
+  r->GetYaxis()->SetBinLabel(4,"resolution (fixed)");
+  r->GetYaxis()->SetBinLabel(5,"p0");
+  r->GetYaxis()->SetBinLabel(6,"p1");
+  r->GetYaxis()->SetBinLabel(7,"p2");
+  r->GetYaxis()->SetBinLabel(8,"p3");
+  r->GetYaxis()->SetBinLabel(9,"chi2");
+  r->GetYaxis()->SetBinLabel(10,"NDF");
+  r->GetYaxis()->SetBinLabel(11,"yield (fit)");
+  r->GetYaxis()->SetBinLabel(12,"yield/event (fit)");
+  r->GetYaxis()->SetBinLabel(13,"peak correction factor");
+  r->GetYaxis()->SetBinLabel(14,"yield (bin counting)");
+  r->GetYaxis()->SetBinLabel(15,"yield/event (bin conunting)");
+  r->GetYaxis()->SetBinLabel(16,"yield/event (bin conunting + tails)");
+
+  TF1 *g[3][4],*gp[3][4],*gb[3][4];
+  TFitResultPtr fr;
+  int status;
+  double pc;
+
+  cerr<<"Please wait a moment while the peak fits are performed."<<endl;
+
+  for(jBF=0;jBF<nBF;jBF++) for(jFR=0;jFR<nFR;jFR++){//loops: change the order of the residual background polynomial and the fit region
+      if(!jFR){A=0.995; B=1.06;}//fr0
+      else if(jFR==1){A=1.; B=1.06;}//fr1
+      else if(jFR==2){A=0.995; B=1.07;}//fr2
+      else if(jFR==3){A=1.; B=1.07;}//fr3
+
+      if(!jBF){//pol1
+       g[jBF][jFR]=new TF1(Form("fit_pol1_fr%i",jFR),"[0]*TMath::Voigt(x-[1],[3],[2])+[4]+[5]*x",A,B);
+       gb[jBF][jFR]=new TF1(Form("fit_back_pol1_fr%i",jFR),"pol1",A,B);
+       gp[jBF][jFR]=new TF1(Form("fit_peak_pol1_fr%i",jFR),"[0]*TMath::Voigt(x-[1],[3],[2])",A,B);
+      }else if(jBF==1){//pol2
+       g[jBF][jFR]=new TF1(Form("fit_pol2_fr%i",jFR),"[0]*TMath::Voigt(x-[1],[3],[2])+[4]+[5]*x+[6]*x*x",A,B);
+       gb[jBF][jFR]=new TF1(Form("fit_back_pol2_fr%i",jFR),"pol2",A,B);
+       gp[jBF][jFR]=new TF1(Form("fit_peak_pol2_fr%i",jFR),"[0]*TMath::Voigt(x-[1],[3],[2])",A,B);
+      }else if(jBF==2){//pol3
+       g[jBF][jFR]=new TF1(Form("fit_pol3_fr%i",jFR),"[0]*TMath::Voigt(x-[1],[3],[2])+[4]+[5]*x+[6]*x*x+[7]*x*x*x",A,B);
+       gb[jBF][jFR]=new TF1(Form("fit_back_pol3_fr%i",jFR),"pol3",A,B);
+       gp[jBF][jFR]=new TF1(Form("fit_peak_pol3_fr%i",jFR),"[0]*TMath::Voigt(x-[1],[3],[2])",A,B);
+      }
+
+      for(j=0;j<100;j++){//initial fit of residual background
+       fr=hb->Fit(gb[jBF][jFR],"RSQ");
+       status=fr->Status();
+       if(!status) break;
+      }
+      urb=gb[jBF][jFR]->IntegralError(Imin,Imax)/dx;//estimated uncertainty of residual background integral
+
+      j=hs->GetXaxis()->FindBin(Imin+1e-5); k=hs->GetXaxis()->FindBin(Imax-1e-5);
+      g[jBF][jFR]->SetParameter(0,hs->Integral(j,k)*dx);//initial value of peak integral
+      g[jBF][jFR]->SetParameter(1,1.019455); g[jBF][jFR]->FixParameter(1,1.019455);//fix mass to vacuum value
+      g[jBF][jFR]->SetParameter(2,0.00426); g[jBF][jFR]->FixParameter(2,0.00426);//fix width to vacuum value
+      g[jBF][jFR]->SetParameter(3,0.0011); g[jBF][jFR]->FixParameter(3,0.0011);//fix resolution to 1.1 MeV/c^2 (will not be changed)
+      for(j=0;j<=jBF+1;j++){//fix residual background to the initial fit
+       A=gb[jBF][jFR]->GetParameter(j);
+       g[jBF][jFR]->SetParameter(j+4,A);
+       g[jBF][jFR]->SetParError(j+4,gb[jBF][jFR]->GetParError(j));
+       g[jBF][jFR]->FixParameter(j+4,A);
+      }
+
+      cerr<<"Fitting pol"<<jBF+1<<"_fr"<<jFR<<endl;
+
+      for(j=0;j<100;j++){//fit with only the peak integral free
+       fr=hs->Fit(g[jBF][jFR],"RSQ");
+       status=fr->Status();
+       if(!status) break;
+      }
+
+      g[jBF][jFR]->ReleaseParameter(2);//release width parameter
+      for(j=0;j<100;j++){
+       fr=hs->Fit(g[jBF][jFR],"RSQ");
+       status=fr->Status();
+       if(!status) break;
+      }
+
+      g[jBF][jFR]->ReleaseParameter(1);//release mass parameter
+      for(j=0;j<100;j++){
+       fr=hs->Fit(g[jBF][jFR],"RSQ");
+       status=fr->Status();
+       if(!status) break;
+      }
+
+      g[jBF][jFR]->ReleaseParameter(4);//release residual background constant parameter
+      for(j=0;j<100;j++){
+       fr=hs->Fit(g[jBF][jFR],"RSQ");
+       status=fr->Status();
+       if(!status) break;
+      }
+
+      for(j=1;j<=gb[jBF][jFR]->GetNpar();j++) g[jBF][jFR]->ReleaseParameter(4+j);//release other residual background parameters
+      for(j=0;j<100;j++){
+       fr=hs->Fit(g[jBF][jFR],"RSQ");
+       status=fr->Status();
+       if(!status) break;
+      }
+
+      for(j=0;j<100;j++){//redo fit using "I" option, typically does not affect the result very much, but best to do it anyway.  Comment this out if you are debugging and want to save time.
+       fr=hs->Fit(g[jBF][jFR],"RSQI");
+       status=fr->Status();
+       if(!status) break;
+      }
+
+      fout->cd();
+      g[jBF][jFR]->Write();//save the combined fit
+
+      for(j=0;j<gb[jBF][jFR]->GetNpar();j++) gb[jBF][jFR]->SetParameter(j,g[jBF][jFR]->GetParameter(j+4));//put the background parameters into the function gb
+      for(j=0;j<4;j++) gp[jBF][jFR]->SetParameter(j,g[jBF][jFR]->GetParameter(j));//put the peak parameters into the funciton gp
+
+      j=nFR*jBF+jFR+1;//bin for the x-axis of results histogram
+      r->GetXaxis()->SetBinLabel(j,Form("pol%i_fr%i",jBF+1,jFR));
+      for(k=0;k<g[jBF][jFR]->GetNpar();k++){//store the values of the peak parameters
+       r->SetBinContent(j,k+1,g[jBF][jFR]->GetParameter(k));
+       r->SetBinError(j,k+1,g[jBF][jFR]->GetParError(k));
+      }
+      r->SetBinContent(j,9,g[jBF][jFR]->GetChisquare());//store chi^2
+      r->SetBinContent(j,10,g[jBF][jFR]->GetNDF());//store number of degrees of freedom
+
+      A=r->GetBinContent(j,1)-gp[jBF][jFR]->Integral(0.,2*0.493667);//subtract integral of peak below kinematic cutoff
+      UA=A*r->GetBinError(j,1)/r->GetBinContent(j,1);
+      r->SetBinContent(j,11,A/dx); r->SetBinError(j,11,UA/dx);//peak yield extracted from the fit
+      r->SetBinContent(j,12,A/dx/nevt/dy); r->SetBinError(j,12,UA/dx/nevt/dy);//peak yield per event extracted from the fit, corrected for dy
+
+      B=gp[jBF][jFR]->Integral(Imin,Imax);
+      pc=B/A;//peak correction factor: the yield in the interval (Imin,Imax) divided by the total integral above the kinematic cutoff, used to correct the yield from bin counting to account for the yield in the tails
+      r->SetBinContent(j,13,pc);//store peak correction factor
+
+      A=UA=0;
+      for(k=hs->GetXaxis()->FindBin(Imin+1e-5);k<=hs->GetXaxis()->FindBin(Imax-1e-5);k++){//yield from bin counting
+       A+=hs->GetBinContent(k);
+       UA+=pow(hs->GetBinError(k),2);
+      }
+      A-=gb[jBF][jFR]->Integral(Imin,Imax)/dx;//subtract residual background integral
+      UA+=urb*urb;
+
+      r->SetBinContent(j,14,A); r->SetBinError(j,14,sqrt(UA));//yield from bin counting
+      r->SetBinContent(j,15,A/nevt/dy); r->SetBinError(j,15,sqrt(UA)/nevt/dy);//yield per event from bin counting, corrected for dy
+      r->SetBinContent(j,16,A/nevt/pc/dy); r->SetBinError(j,16,sqrt(UA)/nevt/pc/dy);//yield per event from bin counting, corrected for dy, corrected to account for yield in tails
+    }
+
+  fout->cd();
+  r->Write();//save the results histogram in case it is needed later
+
+  delete c;
+
+  //----- plot the fits -----
+
+  c=new TCanvas("c","",10,10,1000,1000);
+  c->SetFillColor(0);
+
+  p1=new TPad("p1_2","",0,0.5,0.5,1);
+  p1->SetFillColor(0);
+  p1->SetRightMargin(0.05);
+
+  p2=new TPad("p2_2","",0.5,0.5,1,1);
+  p2->SetFillColor(0);
+  p2->SetRightMargin(0.05);
+
+  p3=new TPad("p3_2","",0,0,0.5,0.5);
+  p3->SetFillColor(0);
+  p3->SetRightMargin(0.05);
+
+  TPad* p4=new TPad("p4_2","",0.5,0,1,0.5);
+  p4->SetFillColor(0);
+
+  hs_plot->GetXaxis()->SetRangeUser(0.995,1.07-1e-5);
+  hs_plot->SetLineColor(1); hs_plot->SetMarkerColor(1);
+  hs_plot->SetTitle(hs->GetTitle());
+  hs_plot->SetXTitle(hs->GetXaxis()->GetTitle());
+
+  for(jBF=0;jBF<nBF;jBF++) for(jFR=0;jFR<nFR;jFR++){
+      //assign styles and colors to the functions
+      g[jBF][jFR]->SetNpx(300);
+      if(!jFR) g[jBF][jFR]->SetLineColor(2);
+      else if(jFR==1) g[jBF][jFR]->SetLineColor(TColor::GetColor("#ffaa00"));
+      else if(jFR==2) g[jBF][jFR]->SetLineColor(TColor::GetColor("#238e23"));
+      else if(jFR==3) g[jBF][jFR]->SetLineColor(4);
+      gb[jBF][jFR]->SetLineColor(g[jBF][jFR]->GetLineColor());
+      gb[jBF][jFR]->SetLineStyle(2);
+    }
+
+  p1->cd();
+  hs_plot->Draw();
+  for(jFR=0;jFR<nFR;jFR++) gb[0][jFR]->Draw("same");
+  for(jFR=0;jFR<nFR;jFR++) g[0][jFR]->Draw("same");
+  t1=new TLatex(0.94,0.87,"linear RB (pol1)"); t1->SetTextAlign(32);
+  t1->SetNDC();
+  t1->Draw();
+
+  p2->cd();
+  hs_plot->Draw();
+  for(jFR=0;jFR<nFR;jFR++) gb[1][jFR]->Draw("same");
+  for(jFR=0;jFR<nFR;jFR++) g[1][jFR]->Draw("same");
+  t2=new TLatex(0.94,0.87,"quadratic RB (pol2)"); t2->SetTextAlign(32);
+  t2->SetNDC();
+  t2->Draw();
+
+  p3->cd();
+  hs_plot->Draw();
+  for(jFR=0;jFR<nFR;jFR++) gb[2][jFR]->Draw("same");
+  for(jFR=0;jFR<nFR;jFR++) g[2][jFR]->Draw("same");
+  TLatex* t3=new TLatex(0.94,0.87,"cubic RB (pol3)"); t3->SetTextAlign(32);
+  t3->SetNDC();
+  t3->Draw();
+
+  //additional information in pad p4
+
+  p4->cd();
+  TLine* dummy1=new TLine(0,0,1,0);
+  dummy1->SetLineColor(1); dummy1->SetLineWidth(2);
+  TLine* dummy2=new TLine(0,0,1,0);
+  dummy2->SetLineColor(1); dummy2->SetLineWidth(2); dummy2->SetLineStyle(gb[0][0]->GetLineStyle());
+  TLegend* L=new TLegend(0.3,0.69,0.7,0.99);
+  L->SetFillColor(0);
+  L->AddEntry(g[0][0],"Fit Region 0 (0.995,1.06)","l");
+  L->AddEntry(g[0][1],"Fit Region 1 (1,1.06)","l");
+  L->AddEntry(g[0][2],"Fit Region 2 (0.995,1.07)","l");
+  L->AddEntry(g[0][3],"Fit Region 3 (1,1.07)","l");
+  L->AddEntry(dummy1,"Combined Fit","l");
+  L->AddEntry(dummy2,"Residual Background Fit","l");
+  L->Draw();
+
+  //----- check results -----
+
+  double ex_yield=0;
+  if(!strcmp(system,"pp276")) ex_yield=1.416e-3;
+  else if(!strcmp(system,"pp7")) ex_yield=1.706e-3;
+  else if(!strcmp(system,"PbPb276")) ex_yield=0.245;
+  else{
+    cerr<<"Warning in QAphi(): Unknown collision system "<<system<<".  Expected yield not known."<<endl;
+    ex_yield=1e10;
+  }
+  double yield=r->GetBinContent(8,16);
+  int status_yield;
+  if(yield/ex_yield<1.5 && yield/ex_yield>0.5) status_yield=0;//OK
+  else if(yield/ex_yield<5 && yield/ex_yield>0.2) status_yield=1;//problem
+  else status_yield=2;//big problem
+
+  double ex_mass=1.019455;
+  double mass=r->GetBinContent(8,2);
+  int status_mass;
+  if(fabs(mass-ex_mass)<1) status_mass=0;//OK
+  else if(fabs(mass-ex_mass)<3) status_mass=1;//problem
+  else status_mass=2;//big problem
+
+  double ex_width=0.00426;
+  double width=r->GetBinContent(8,3);
+  int status_width;
+  if(fabs(width-ex_width)<1) status_width=0;//OK
+  else if(fabs(width-ex_width)<2) status_width=1;//problem
+  else status_width=2;//big problem;
+
+  double tx=0.005,ty=0.64,dty=0.06;
+
+  TString s;
+
+  j=su->GetAxis(1)->GetFirst(); k=su->GetAxis(1)->GetLast();
+  s.Form("%1.2f < #it{p}_{T} < %1.2f GeV/#it{c}, ",su->GetAxis(1)->GetBinLowEdge(j),su->GetAxis(1)->GetBinLowEdge(k+1));
+  j=su->GetAxis(2)->GetFirst(); k=su->GetAxis(2)->GetLast();
+  s.Append(Form("%1.2f < #it{y} < %1.2f",su->GetAxis(2)->GetBinLowEdge(j),su->GetAxis(2)->GetBinLowEdge(k+1)));
+  if(!ptOK || !yOK) s.Append(" [PROBLEM]");
+  TLatex* a1=new TLatex(tx,ty,s.Data());
+  SetText(a1,ptOK);
+  a1->Draw();
+
+  TLatex* a2=new TLatex(tx,ty-1*dty,Form("#phi Yield/event = %1.5e #pm %1.5e (%s)",r->GetBinContent(8,16),r->GetBinError(8,16),r->GetXaxis()->GetBinLabel(8)));
+  SetText(a2,status_yield);
+  a2->Draw();
+
+  s.Form("Expected Yield = %1.5e",ex_yield);
+  if(status_yield) s.Append(" [PROBLEM]");
+  TLatex* a3=new TLatex(tx+0.1,ty-2*dty,s.Data());
+  SetText(a3,status_yield);
+  a3->Draw();
+
+  TLatex* a4=new TLatex(tx,ty-3*dty,Form("#phi Mass = %1.5f #pm %1.5f GeV/#it{c}^{2}",r->GetBinContent(8,2),r->GetBinError(8,2)));
+  SetText(a4,status_mass);
+  a4->Draw();
+
+  s.Form("PDG Mass = 1.019455");
+  if(status_mass) s.Append(" [PROBLEM]");
+  TLatex* a5=new TLatex(tx+0.1,ty-4*dty,s.Data());
+  SetText(a5,status_mass);
+  a5->Draw();
+
+  TLatex* a6=new TLatex(tx,ty-5*dty,Form("#phi Width = %1.5f #pm %1.5f MeV/#it{c}^{2}",r->GetBinContent(8,3)*1000.,r->GetBinError(8,3)*1000.));
+  SetText(a6,status_width);
+  a6->Draw();
+
+  s.Form("PDG Width = 4.26");
+  if(status_width) s.Append(" [PROBLEM]");
+  TLatex* a7=new TLatex(tx+0.1,ty-6*dty,s.Data());
+  SetText(a7,status_width);
+  a7->Draw();
+
+  if(!ptOK || !yOK || status_yield || status_mass || status_width){
+    status=2;
+    s.Form("THERE ARE PROBLEMS!!!");
+  }else{
+    status=0;
+    s.Form("Things look OK.");
+  }
+
+  TLatex* a8=new TLatex(tx,ty-7*dty,s.Data());
+  SetText(a8,status); a8->SetTextSize(0.05);
+  a8->Draw();
+
+  if(status) s.Form("PLEASE CHECK THE ERROR MESSAGES.");
+  else s.Form("");
+
+  TLatex* a9=new TLatex(tx,ty-8*dty,s.Data());
+  SetText(a9,status); a9->SetTextSize(0.05);
+  a9->Draw();
+
+  c->cd();
+  p1->Draw();
+  p2->Draw();
+  p3->Draw();
+  p4->Draw();
+
+  c->SaveAs(Form("%s.eps)",name_fout));
+
+  //----- print out results -----
+
+  j=su->GetAxis(1)->GetFirst(); k=su->GetAxis(1)->GetLast();
+  printf("%1.2f < pT < %1.2f GeV/c, ",su->GetAxis(1)->GetBinLowEdge(j),su->GetAxis(1)->GetBinLowEdge(k+1));
+  j=su->GetAxis(2)->GetFirst(); k=su->GetAxis(2)->GetLast();
+  printf("%1.2f < y < %1.2f\n",su->GetAxis(2)->GetBinLowEdge(j),su->GetAxis(2)->GetBinLowEdge(k+1));
+  printf("phi Yield/event = %1.5e +/- %1.5e (%s)\n",r->GetBinContent(8,16),r->GetBinError(8,16),r->GetXaxis()->GetBinLabel(8));
+  if(status_yield) printf("*** PROBLEM: phi yield too far from expected value (%1.5e).\n",ex_yield);
+  printf("phi Mass = %1.5f +/- %1.5f GeV/c^2 (PDG Value = 1.019455)\n",r->GetBinContent(8,2),r->GetBinError(8,2));
+  if(status_mass) printf("*** PROBLEM: phi mass too far from expected value.\n");
+  printf("phi Width = %1.5f +/- %1.5f MeV/c^2 (PDG Value = 4.26)\n",r->GetBinContent(8,3)*1000.,r->GetBinError(8,3)*1000.);
+  if(status_width) printf("*** PROBLEM: phi width too far from expected value.\n");
+  if(status) printf("\n*** One or more parameters does not have the expected value.  Please check the error messages and read the instructions at the bottom of the macro file.\n");
+
+  //----- finish -----
+
+  cerr<<"\nMacro has finished successfully."<<endl;
+
+  fin->Close();
+  fout->Close();
+
+  return;
+}
+
+
+bool SetPtRange(THnSparse* h){
+  TAxis* a=h->GetAxis(1);
+  double min=0.5,max=1.5;
+
+  int b1=a->FindBin(1.00001*min);
+  int b2=a->FindBin(0.99999*max);
+  a->SetRange(b1,b2);
+  b1=a->GetFirst();
+  b2=a->GetLast();
+
+  //min and/or max are not bin boundaries
+  if(fabs(a->GetBinLowEdge(b1)-min)>1e-5*min || fabs(a->GetBinLowEdge(b2+1)-max)>1e-5*max){
+    cerr<<"Error in QAphi(): pT range cannot be set to ("<<min<<","<<max<<").  The yield may therefore be different from the expected value."<<endl;
+    return false;
+  }
+
+  return true;
+}
+
+bool SetRapidityRange(THnSparse* h,double& dy){
+  TAxis* a=h->GetAxis(2);
+  double min=-0.5,max=0.5;
+
+  int b1=a->FindBin(min+1e-5);
+  int b2=a->FindBin(max-1e-5);
+  a->SetRange(b1,b2);
+  b1=a->GetFirst();
+  b2=a->GetLast();
+  dy=a->GetBinLowEdge(b2+1)-a->GetBinLowEdge(b1);
+
+  //min and/or max are not bin boundaries
+  if(fabs(a->GetBinLowEdge(b1)-min)>1e-5 || fabs(a->GetBinLowEdge(b2+1)-max)>1e-5){
+    cerr<<"Error in QAphi(): rapidity range cannot be set to |y|<0.5.  Although this macro does correct for the rapidity range dy, using a different rapidity range may cause the yield to be different from the expected value."<<endl;
+    return false;
+  }
+
+  return true;
+}
+
+
+void SetText(TLatex* t,int flag){
+  t->SetTextSize(0.04);
+  if(!flag) t->SetTextColor(TColor::GetColor("#238e23"));
+  else if(flag==1) t->SetTextColor(TColor::GetColor("#cc5500"));
+  else t->SetTextColor(2);
+  return;
+}
+
+
+void SetText(TLatex* t,bool flag){
+  if(flag) SetText(t,0);
+  else SetText(t,2);
+  return;
+}
+
+
+/*-----------------------------------------------
+
+INSTRUCTIONS:
+
+*Introduction:
+
+This macro is designed to read the output of the analysis task to check that the yield per event, mass, and width of the phi meson in the production are within acceptable ranges.  The analysis task constructs invariant-mass distributions of charged-kaon pairs in the vicinity of the mass of the phi meson (1.019455 GeV/c^2).  The branching ratio of the phi->K-K+ decay is 0.489.  In this macro, a combinatorial background is constructed using the like-charge distributions (K-K- and K+K+) and subtracted from the unlike-charge (K-K+) distribution.  This leaves a phi peak sitting on top of a residual background.  The background-subtracted distribution is fit with a Voigtian peak (a convolution of a non-relativistic Breit-Wigner peak with a Gaussian to account for detector effects) plus a polynomial to describe the residual background.  Because these fits sometimes fail, a total of 12 fits are performed.  There are 3 choices of residual background polynomial (linear, quadratic, or cubic) and 4 different fit regions (defined in the macro).  For more information on the fitting procedure, see refs. [2] and [3].  The yield, mass, and width of the phi meson extracted from the default fit (quadratic residual background, fit region 3) are compared to the expected values.
+
+*Output:
+
+This macro produces an output PDF file and an output ROOT file.
+
+**Output PDF:
+
+The output PDF consists of three canvases, organized into panels as follows:
+
+Canvas 1: [[1a][1b][1c]]
+Canvas 2: [[2a][2b]]
+Canvas 3: [[3a][3b]]
+          [[3c][3d]]
+
+The contents of the panels is:
+1a.) histogram hEventStat (see description below)
+1b.) histogram hAEventsVsMulti (see description below)
+1c.) histogram hAEventsVsMulti, logarithmic y-axis
+2a.) unlike-charge invariant-mass distribution and like-charge combinatorial background
+2b.) background-subtracted invariant-mass distribution
+3a.) background-subtracted invariant-mass distribution with fits.  The fits shown here assume a linear residual background.  A fit is shown for each of the four different fit regions.  Fit region 3 (blue) is plotted on top.  The dashed lines are the residual background.
+3b.) same as 3a, but with a quadratic residual background.  Note that the blue fit in this panel is the default fit, from which the yield, mass, and width are extracted.
+3c.) same as 3a, but with a cubic residual background
+3d.) a legend describing the fit functions and a printout of information.  The printed information should be green.  If it is not, there is a problem.
+
+**Output ROOT file:
+
+The output ROOT file contains the following:
+1.) histogram hEventStat (see description below)
+2.) histogram hAEventsVsMulti (see description below)
+3.) histogram mass_unlike: the unlike-charge invariant-mass distribution
+4.) histogram mass_like: the like-charge combinatorial background
+5.) histogram mass_signal: background-subtracted invariant-mass distribution (mass_unlike - mass_like)
+6.) functions fit_pol*_fr*: the fit functions for mass_signal.  pol1: linear residual background, pol2: quadratic residual background, pol3: cubic residual background.  fr* indicates the fit region.
+7.) histogram fit_results:  contains the results of the different fits (see below)
+
+***Histogram fit_results:
+
+This is an array of numbers (with uncertainties) stored in a TH2D.  The x-axis is the name of the fit (e.g., "pol2_fr3") and the y-axis gives the name of the quantity stored.  The quantities stored (numbered along the y-axis) are:
+
+1.) peak integral: integral of the phi peak in the fit function (parameter [0])
+2.) peak mass: parameter [1]
+3.) peak width: parameter [2]
+4.) peak resolution: parameter [3]: currently fixed to 1.1 MeV/c^2, but stored anyway
+5-8.) coefficients of the residual background polynomial (parameters [4] through [7])
+9.) fit chi^2
+10.) fit number of degrees of freedom
+11.) yield extracted from the fit: parameter [0] minus the yield below the kinematic cutoff and corrected for the width of the invariant-mass bins
+12.) yield/event (fit): the value above divided by the number of events and the width of the rapidity range (dy)
+13.) peak correction factor: the fraction of the peak that lies inside the range (1.01,1.03) GeV/c^2
+14.) yield (bin counting): the integral of the mass_signal histogram for (1.01,1.03) minus the integral of the residual background over the same interval
+15.) yield/event (bin counting) the value above divided by the number of events and the width of the rapidity range (dy)
+16.) yield/event (bin counting + tails) the value above further corrected by dividing by the peak correction factor (now accounting for the yield in the tails)
+
+The default fit is "pol2_fr3".  Quantity 2 is reported as the mass, quantity 3 is reported as the width, and quantity 16 is reported as the yield.  So, to find the default mass (which is printed out) do fit_results->GetBinContent(8,2).  The default width is fit_results->GetBinContent(8,3) and the default yield is fit_results->GetBinContent(8,16).
+
+*More Information:
+
+In the analysis task, the kaons should be selected using track selection cuts and a 3-sigma cut (about the kaon mean) on TPC dE/dx.  See the analysis task for the exact cuts.
+
+The output of the analysis task is expected to contain:
+1.) the histogram hEventStat, which gives the number of events that pass different selection criteria
+2.) the histogram hAEventsVsMulti, which gives the number of events in each multiplicity (for pp) or centrality (for Pb-Pb) bin
+3.) three THnSparse histograms with three dimensions and the suffixes "Unlike", "LikePP", and "LikeMM"
+
+For each THnSparse, axis 0 is the KK invariant mass, axis 1 is the transverse momentum, and axis 2 is the rapidity.  This macro attempts to select the range 0.5<pT<1.5 GeV/c; if this is not possible an error message will be generated, as the expected yields were calculated assuming that pT range.  This macro uses the full rapidity range on axis 2.  The expected yields were calculated assuming a rapidity range of |y|<0.5.  The macro corrects for the width of the rapidity range (dy), but using a different rapidity range may cause the yield to deviate from the expected value.
+
+*Expected Values
+
+**Expected Yields
+
+If the yield is less than 50% or more than 150% of the expected value, it is flagged as an orange problem.  If the yield deviates from the expected value by more than a factor of 5, it is flagged as a red problem.
+
+The expected yields differ depending on the collision system.
+
+pp Collisions at 7 TeV (system="pp7"): The expected yield of 1.706e-3 is computed based on the published (ref. [1]) phi spectrum.  The yield per event for 0.5<pT<1.5 GeV/c is multiplied by the efficiency and the branching ratio.
+
+pp Collisions at 2.76 TeV (system="pp276"): The expected yield of 1.416e-3 is computed by extrapolating from the pp 7 TeV yield, assuming that the phi yield scales as s^0.1 (energy dependence taken from ref. [2]).
+
+Pb-Pb Collisions at 2.76 TeV (system="PbPb276"): The expected yield of 0.245 is computed based on the nearly published (ref. [2]) phi spectrum.
+
+**Expected Mass
+
+The expected phi mass is the PDG value: 1.019455 GeV/c^2.  A deviation from this value of more than 1 MeV/c^2 is flagged as an orange problem.  A deviation from this value of more than 3 MeV/c^2 is flagged as a red problem.
+
+**Expected Width
+
+The expected phi width is the PDG value: 4.26 MeV/c^2.  A deviation from this value of more than 1 MeV/c^2 is flagged as an orange problem.  A deviation from this value of more than 3 MeV/c^2 is flagged as a red problem.
+
+*Troubleshooting
+
+If a parameter does not have the expected value, it will be flagged.  In the output PDF file, a problematic parameter will be colored orange or red, with red indicating a larger deviation from the expected value.
+
+If the yield, mass, or width deviates from the expected value, you should check to see if the default fit is bad.  Look at the blue fit in the upper right plot of the third canvas in the PDF file (panel 3b).  This is the default fit.  Does it describe the data?  Does the residual background (dashed line) appear to be reasonable.  If the default fit looks bad, look through the other fits and see if you can find a fit that looks OK.  You can read out the yield, mass, and width from the fit_results histogram.  Compare these values to the expected values.  It might also make sense to compare the results of all fits to the expected value.  Are the fit parameters always outside the acceptable range?
+
+If the phi yield deviates from the expected value, the following should be noted.  The yield depends on
+1.) collision system and energy
+2.) the triggers used
+3.) multiplicity or centrality range
+4.) cuts used to select the decay daughters
+5.) pT and rapidity range for phi mesons
+
+The expected yields were calculated for
+1.) pp and Pb-Pb collisions at 2.76 TeV, pp collisions at 7 TeV
+2.) minimum-bias triggers
+3.) all multiplicities for pp, centrality 0-90% for Pb-Pb
+4.) the track selection cuts used in the original version of the analysis task
+5.) 0.5<pT<1.5 GeV/c and |y|<0.5 for phi mesons
+
+Changes to any of these criteria can cause the yield to change.  If the criteria you are using do not exactly match the criteria for the expected yield, you will need to correct for those differences.  The macro corrects for the rapidity range, so SMALL changes to the rapidity window will be accounted for to some extent.
+
+If the width deviates from the expected value, it is possible that the resolution is not 1.1 MeV/c^2 (the value to which it is currently fixed).  Studying this issue is beyond the scope of these instructions, but you should be aware of the possibility.
+
+*References
+
+[1]: B. Abelev et al. (ALICE Collaboration), Eur. Phys. J. C 72, 2183 (2012)
+[2]: Paper in preparation: "K*(892)^0 and phi(1020) resonances in Pb-Pb collisions at 2.76 TeV", by the ALICE Collaboration, intended for publication in Phys. Rev. C (2014)
+[3]: Analysis Note: A. G. Knospe, "Yield of phi mesons at low pT in Pb-Pb collisions at 2.76 TeV (2010 data)", ALICE-ANA-2012-300, https://aliceinfo.cern.ch/Notes/node/42
+  -----------------------------------------------*/
+