Macro to compite pion, kaon and proton spectra with nsigma cuts (Leonardo)
authorfprino <prino@to.infn.it>
Mon, 11 Aug 2014 14:57:02 +0000 (16:57 +0200)
committerfprino <prino@to.infn.it>
Mon, 11 Aug 2014 14:57:02 +0000 (16:57 +0200)
PWGLF/SPECTRA/PiKaPr/ITSsa/MakeCorrectedSpectraITSsaNsigma.C [new file with mode: 0644]

diff --git a/PWGLF/SPECTRA/PiKaPr/ITSsa/MakeCorrectedSpectraITSsaNsigma.C b/PWGLF/SPECTRA/PiKaPr/ITSsa/MakeCorrectedSpectraITSsaNsigma.C
new file mode 100644 (file)
index 0000000..5f3eaff
--- /dev/null
@@ -0,0 +1,897 @@
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TF1.h>
+#include <TPad.h>
+#include <TGraphErrors.h>
+#include <TROOT.h>
+#include <TFile.h>
+#include <TTree.h>
+#include <TGrid.h>
+#include <TGridResult.h>
+#include <TMath.h>
+#include <TCanvas.h>
+#include <TStyle.h>
+#include <TLatex.h>
+#include <TImage.h>
+#include <TPaveText.h>
+#include <TLine.h>
+#endif
+
+void RatioPlot(TH1F **num,TH1F **den,TString t1,TString t2,TString opt,Double_t Low[],Double_t Up[]);
+void SetDrawAtt(Int_t markerstyle,Int_t markercolor,Int_t markersize,Int_t linecolor,Int_t linewidth,TH1 *h1);
+void ResetOutOfRange(TH1F *histo, Int_t ipart, Double_t lowRange[3], Double_t upRange[3]);
+
+void MakeCorrectedSpectraITSsaNsigma( TString period="LHC10d",
+                                     TString MCname="LHC10f6a",
+                                     Float_t DCAcut=7,
+                                     Int_t lowmult=-1,
+                                     Int_t upmult=-1
+                                     )    
+{
+
+  // parameters
+
+  Bool_t MakeMCPlots=1;
+  Bool_t MakeSystErr=0;
+  Bool_t NormalizeData=1;
+  Bool_t NormalizeToINEL=0;
+  Bool_t SaveOnlyAsymm=0;
+  Bool_t CorrectOnlyForEfficiency=0;
+  Bool_t MarekNorm=0;
+  Bool_t DCACorrection=1;
+  Bool_t GFCorrection=1;
+
+  //MC Histos
+  TH1F *fHistPrimMCposBefEvSel[3];
+  TH1F *fHistPrimMCnegBefEvSel[3];
+  TH1F *fHistPrimMCpos[3];
+  TH1F *fHistPrimMCneg[3];
+  TH1F *hHistPosNSigma[3];
+  TH1F *hHistNegNSigma[3];
+  TH1F *hHistPosNSigmaMean[3];
+  TH1F *hHistNegNSigmaMean[3];
+  TH1F *hHistPosNSigmaPrim[3];
+  TH1F *hHistNegNSigmaPrim[3];
+  TH1F *hHistPosNSigmaPrimMean[3];
+  TH1F *hHistNegNSigmaPrimMean[3];
+  TH1F *hHistPosNSigmaPrimMC[3];
+  TH1F *hHistNegNSigmaPrimMC[3];
+  TH1F *hHistPosNSigmaPrimMCMean[3];
+  TH1F *hHistNegNSigmaPrimMCMean[3];
+  TH1F *hHistPosNSigmaMC[3];
+  TH1F *hHistNegNSigmaMC[3];
+  TH1F *hHistPosNSigmaMCMean[3];
+  TH1F *hHistNegNSigmaMCMean[3];
+  TH1F *hCorrFactPos[3];
+  TH1F *hCorrFactNeg[3];
+  TH1F *hCorrFactPosMean[3];
+  TH1F *hCorrFactNegMean[3];
+  TH1F *hEffPos[3];
+  TH1F *hEffNeg[3];
+  TH1F *hEffPosMean[3];
+  TH1F *hEffNegMean[3];
+  
+  //DATA Histos
+  TH1F *hHistPosNSigmaDATA[3];
+  TH1F *hHistNegNSigmaDATA[3];
+  TH1F *hHistPosNSigmaMeanDATA[3];
+  TH1F *hHistNegNSigmaMeanDATA[3];
+  TH1F *hITSsaPos[3];
+  TH1F *hITSsaNeg[3];
+  TH1F *hITSsaPosMean[3];
+  TH1F *hITSsaNegMean[3];
+  TH1F *hITSsaRawPos[3];
+  TH1F *hITSsaRawNeg[3];
+  TH1F *hITSsaRawPosMean[3];
+  TH1F *hITSsaRawNegMean[3];
+  
+  //  Float_t nmin=1.5;
+  Double_t Low[3]={0.10,0.2,0.3};
+  Double_t Up[3]={0.4,0.5,0.7};
+  Double_t lowRange[3]={0.10,0.2,0.3};
+  Double_t upRange[3]={0.6,0.6,0.6};
+
+  
+
+  const Int_t nbinspt=22;
+  Double_t xbins[nbinspt+1]={0.08,0.10,0.12,0.14,0.16,0.18,0.20,0.25,0.30,0.35,0.40,0.45,0.50,0.55,0.60,0.65,0.70,0.75,0.80,0.85,0.90,0.95,1.0};
+    
+  TString fname_MC=Form("./%s/AnalysisResults.root",MCname.Data());
+  cout<<fname_MC.Data()<<endl;
+
+  TString foutMC=Form("outSpectraMC_%s_%s_%.0fDCA.root",period.Data(),MCname.Data(),DCAcut);
+  TFile *outMC=new TFile(foutMC.Data(),"recreate");
+  outMC->Close();
+  delete outMC;
+
+  TString foutDATA;
+  if(CorrectOnlyForEfficiency){
+    foutDATA=Form("outSpectraData_%s_%s_%.0fDCA.root",period.Data(),MCname.Data(),DCAcut);
+  }else{
+    foutDATA=Form("outCorrSpectraData_%s_%s_%.0fDCA.root",period.Data(),MCname.Data(),DCAcut);
+  }
+
+  TString fnameDATA=Form("%s/AnalysisResults.root",period.Data());
+  cout<<fnameDATA.Data()<<endl;
+  TFile *outDATA=new TFile(foutDATA.Data(),"recreate");
+  outDATA->Close();
+  delete outDATA;
+  
+
+  TString lnameMC=Form("clistITSsaMult%ito%i",lowmult,upmult);
+  Printf("\n------------------ READING MC  %s -------------------\n",lnameMC.Data());
+  TFile *finMC = new TFile(fname_MC.Data());
+  TDirectory *dirFileMC=(TDirectory*)finMC->Get("PWG2SpectraITSsa");
+  TList * cOutputMC = (TList*)dirFileMC->Get(lnameMC.Data());
+  //cOutputMC->Print();
+  
+  TH1F *fHistNEventsMC = (TH1F*)cOutputMC->FindObject("fHistNEvents");
+  Float_t normMC=(Float_t)1/fHistNEventsMC->GetBinContent(fHistNEventsMC->FindBin(1.));
+  printf("- NormalizationMC (NEvents): %.0f \n",fHistNEventsMC->GetBinContent(fHistNEventsMC->FindBin(1.))); 
+  
+  for(Int_t ipart=0;ipart<=2;ipart++){
+    fHistPrimMCposBefEvSel[ipart]=(TH1F*)cOutputMC->FindObject(Form("fHistPrimMCposBefEvSel%i",ipart));//primaries generated (before event selection)
+    fHistPrimMCnegBefEvSel[ipart]=(TH1F*)cOutputMC->FindObject(Form("fHistPrimMCnegBefEvSel%i",ipart));
+    fHistPrimMCpos[ipart]=(TH1F*)cOutputMC->FindObject(Form("fHistPrimMCpos%i",ipart)); //primaries generated (after event selection)
+    fHistPrimMCneg[ipart]=(TH1F*)cOutputMC->FindObject(Form("fHistPrimMCneg%i",ipart));
+    hHistPosNSigma[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistPosNSigma%i",ipart));//NSigma recontructed, no MCtruth
+    hHistNegNSigma[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistNegNSigma%i",ipart));
+    hHistPosNSigmaPrimMC[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistPosNSigmaPrimMC%i",ipart));//NSigma recontructed, MCtruth
+    hHistNegNSigmaPrimMC[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistNegNSigmaPrimMC%i",ipart));
+    hHistPosNSigmaMC[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistPosNSigmaMC%i",ipart));//NSigma recontructed, MCtruth
+    hHistNegNSigmaMC[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistNegNSigmaMC%i",ipart));
+    hHistPosNSigmaPrim[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistPosNSigmaPrim%i",ipart));//NSigma recontructed, truth
+    hHistNegNSigmaPrim[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistNegNSigmaPrim%i",ipart));
+    hHistPosNSigmaMean[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistPosNSigmaMean%i",ipart));//NSigmaMean recontructed, no MCtruth
+    hHistNegNSigmaMean[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistNegNSigmaMean%i",ipart));
+    hHistPosNSigmaPrimMCMean[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistPosNSigmaPrimMCMean%i",ipart));//NSigmaMean recontructed, no MCtruth
+    hHistNegNSigmaPrimMCMean[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistNegNSigmaPrimMCMean%i",ipart));
+    hHistPosNSigmaMCMean[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistPosNSigmaMCMean%i",ipart));//NSigmaMean recontructed, no MCtruth
+    hHistNegNSigmaMCMean[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistNegNSigmaMCMean%i",ipart));
+    hHistPosNSigmaPrimMean[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistPosNSigmaPrimMean%i",ipart));//NSigmaMean recontructed, no MCtruth
+    hHistNegNSigmaPrimMean[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistNegNSigmaPrimMean%i",ipart));
+  }
+    
+  for(Int_t ipart=0;ipart<=2;ipart++){
+    hCorrFactPos[ipart]=new TH1F(Form("hCorrFactPos%i",ipart),Form("hCorrFactPos%i",ipart),nbinspt,xbins);
+    hCorrFactNeg[ipart]=new TH1F(Form("hCorrFactNeg%i",ipart),Form("hCorrFactNeg%i",ipart),nbinspt,xbins);
+    hCorrFactPosMean[ipart]=new TH1F(Form("hCorrFactPosMean%i",ipart),Form("hCorrFactPosMean%i",ipart),nbinspt,xbins);
+    hCorrFactNegMean[ipart]=new TH1F(Form("hCorrFactNegMean%i",ipart),Form("hCorrFactNegMean%i",ipart),nbinspt,xbins);
+    hEffPos[ipart]=new TH1F(Form("hEffPos%i",ipart),Form("hEffPos%i",ipart),nbinspt,xbins);
+    hEffNeg[ipart]=new TH1F(Form("hEffNeg%i",ipart),Form("hEffNeg%i",ipart),nbinspt,xbins);
+    hEffPosMean[ipart]=new TH1F(Form("hEffPosMean%i",ipart),Form("hEffPosMean%i",ipart),nbinspt,xbins);
+    hEffNegMean[ipart]=new TH1F(Form("hEffNegMean%i",ipart),Form("hEffNegMean%i",ipart),nbinspt,xbins);
+  }
+    
+  for(Int_t ipart=0;ipart<=2;ipart++){     
+    hHistPosNSigma[ipart]->Sumw2();
+    hHistNegNSigma[ipart]->Sumw2();
+    hHistPosNSigmaMean[ipart]->Sumw2();
+    hHistNegNSigmaMean[ipart]->Sumw2();
+    hHistPosNSigmaPrimMC[ipart]->Sumw2();
+    hHistNegNSigmaPrimMC[ipart]->Sumw2();
+    hHistPosNSigmaMC[ipart]->Sumw2();
+    hHistNegNSigmaMC[ipart]->Sumw2();
+    hHistPosNSigmaPrimMCMean[ipart]->Sumw2();
+    hHistNegNSigmaPrimMCMean[ipart]->Sumw2();
+    hHistPosNSigmaPrim[ipart]->Sumw2();
+    hHistNegNSigmaPrim[ipart]->Sumw2();
+    hHistPosNSigmaPrimMean[ipart]->Sumw2();
+    hHistNegNSigmaPrimMean[ipart]->Sumw2();
+    hHistPosNSigmaMCMean[ipart]->Sumw2();
+    hHistNegNSigmaMCMean[ipart]->Sumw2();
+    
+    hCorrFactPos[ipart]->Divide(hHistPosNSigma[ipart],fHistPrimMCposBefEvSel[ipart],1,1,"B");
+    hCorrFactNeg[ipart]->Divide(hHistNegNSigma[ipart],fHistPrimMCnegBefEvSel[ipart],1,1,"B");
+    hCorrFactPosMean[ipart]->Divide(hHistPosNSigmaMean[ipart],fHistPrimMCposBefEvSel[ipart],1,1,"B");
+    hCorrFactNegMean[ipart]->Divide(hHistNegNSigmaMean[ipart],fHistPrimMCnegBefEvSel[ipart],1,1,"B");
+    
+    // hCorrFactPos[ipart]->Divide(hHistPosNSigmaPrim[ipart],fHistPrimMCposBefEvSel[ipart],1,1,"B");
+    // hCorrFactNeg[ipart]->Divide(hHistNegNSigmaPrim[ipart],fHistPrimMCnegBefEvSel[ipart],1,1,"B");
+    // hCorrFactPosMean[ipart]->Divide(hHistPosNSigmaPrimMean[ipart],fHistPrimMCposBefEvSel[ipart],1,1,"B");
+    // hCorrFactNegMean[ipart]->Divide(hHistNegNSigmaPrimMean[ipart],fHistPrimMCnegBefEvSel[ipart],1,1,"B");
+    
+    //hCorrFactPos[ipart]->Divide(hHistPosNSigmaPrimMC[ipart],fHistPrimMCposBefEvSel[ipart],1,1,"B");
+    //hCorrFactNeg[ipart]->Divide(hHistNegNSigmaPrimMC[ipart],fHistPrimMCnegBefEvSel[ipart],1,1,"B");
+    //hCorrFactPosMean[ipart]->Divide(hHistPosNSigmaPrimMCMean[ipart],fHistPrimMCposBefEvSel[ipart],1,1,"B");
+    //hCorrFactNegMean[ipart]->Divide(hHistNegNSigmaPrimMCMean[ipart],fHistPrimMCnegBefEvSel[ipart],1,1,"B");
+    
+    hCorrFactPos[ipart]->SetTitle(Form("Correction Factor NSigma Pos%i",ipart));
+    hCorrFactNeg[ipart]->SetTitle(Form("Correction Factor NSigma Neg%i",ipart));
+    hCorrFactPosMean[ipart]->SetTitle(Form("Corr Fact NSigma Pos%i",ipart));
+    hCorrFactNegMean[ipart]->SetTitle(Form("Corr Fact NSigma Neg%i",ipart));
+    hCorrFactPos[ipart]->SetName(Form("hCorrFactPos%i",ipart));
+    hCorrFactNeg[ipart]->SetName(Form("hCorrFactNeg%i",ipart));
+    hCorrFactPosMean[ipart]->SetName(Form("hCorrFactPosMean%i",ipart));
+    hCorrFactNegMean[ipart]->SetName(Form("hCorrFactNegMean%i",ipart));
+    SetDrawAtt(21+ipart,ipart+4,1,ipart+4,1,hCorrFactPos[ipart]);
+    SetDrawAtt(25+ipart,ipart+4,1,ipart+4,1,hCorrFactNeg[ipart]);
+    SetDrawAtt(21+ipart,ipart+1,1,ipart+1,1,hCorrFactPosMean[ipart]);
+    SetDrawAtt(25+ipart,ipart+1,1,ipart+1,1,hCorrFactNegMean[ipart]); 
+    ResetOutOfRange(hCorrFactPos[ipart],ipart,lowRange,upRange);
+    ResetOutOfRange(hCorrFactNeg[ipart],ipart,lowRange,upRange);
+    ResetOutOfRange(hCorrFactPosMean[ipart],ipart,lowRange,upRange);
+    ResetOutOfRange(hCorrFactNegMean[ipart],ipart,lowRange,upRange);
+    
+    hHistPosNSigmaPrimMC[ipart]->Sumw2();
+    hHistNegNSigmaPrimMC[ipart]->Sumw2();
+    hHistPosNSigmaPrimMCMean[ipart]->Sumw2();
+    hHistNegNSigmaPrimMCMean[ipart]->Sumw2();
+    hEffPos[ipart]->Divide(hHistPosNSigmaPrimMC[ipart],fHistPrimMCpos[ipart],1,1,"B");
+    hEffNeg[ipart]->Divide(hHistNegNSigmaPrimMC[ipart],fHistPrimMCneg[ipart],1,1,"B");
+    hEffPosMean[ipart]->Divide(hHistPosNSigmaPrimMCMean[ipart],fHistPrimMCpos[ipart],1,1,"B");
+    hEffNegMean[ipart]->Divide(hHistNegNSigmaPrimMCMean[ipart],fHistPrimMCneg[ipart],1,1,"B");
+    hEffPos[ipart]->SetTitle(Form("Efficiency NSigma Pos%i",ipart));
+    hEffNeg[ipart]->SetTitle(Form("Efficiency NSigma Neg%i",ipart));
+    hEffPosMean[ipart]->SetTitle(Form("Efficiency NSigma Asymm Pos%i",ipart));
+    hEffNegMean[ipart]->SetTitle(Form("Efficiency NSigma Asymm Neg%i",ipart));
+    hEffPos[ipart]->SetName(Form("hEffPos%i",ipart));
+    hEffNeg[ipart]->SetName(Form("hEffNeg%i",ipart));
+    hEffPosMean[ipart]->SetName(Form("hEffPosMean%i",ipart));
+    hEffNegMean[ipart]->SetName(Form("hEffNegMean%i",ipart));
+    
+    SetDrawAtt(21+ipart,ipart+4,1,ipart+4,1,hEffPos[ipart]);
+    SetDrawAtt(25+ipart,ipart+4,1,ipart+4,1,hEffNeg[ipart]);
+    SetDrawAtt(21+ipart,ipart+1,1,ipart+1,1,hEffPosMean[ipart]);
+    SetDrawAtt(25+ipart,ipart+1,1,ipart+1,1,hEffNegMean[ipart]);
+    ResetOutOfRange(hEffPos[ipart],ipart,lowRange,upRange);
+    ResetOutOfRange(hEffNeg[ipart],ipart,lowRange,upRange);
+    ResetOutOfRange(hEffPosMean[ipart],ipart,lowRange,upRange);
+    ResetOutOfRange(hEffNegMean[ipart],ipart,lowRange,upRange);
+      
+      
+    hCorrFactPos[ipart]->SetMinimum(0.);
+    hCorrFactNeg[ipart]->SetMinimum(0.);
+    hCorrFactPosMean[ipart]->SetMinimum(0.);
+    hCorrFactNegMean[ipart]->SetMinimum(0.);
+    hEffPos[ipart]->SetMinimum(0.);
+    hEffNeg[ipart]->SetMinimum(0.);
+    hEffPosMean[ipart]->SetMinimum(0.);
+    hEffNegMean[ipart]->SetMinimum(0.);
+    hCorrFactPos[ipart]->SetMaximum(1.);
+    hCorrFactNeg[ipart]->SetMaximum(1.);
+    hCorrFactPosMean[ipart]->SetMaximum(1.);
+    hCorrFactNegMean[ipart]->SetMaximum(1.);
+    hEffPos[ipart]->SetMaximum(1.);
+    hEffNeg[ipart]->SetMaximum(1.);
+    hEffPosMean[ipart]->SetMaximum(1.);
+    hEffNegMean[ipart]->SetMaximum(1.);
+  }
+    
+  if(MakeMCPlots){
+    TCanvas *cCorr=new TCanvas("cCorr","cCorr");
+    cCorr->SetGridy();
+    for(Int_t ipart=0;ipart<=2;ipart++){
+      if(ipart==0)hCorrFactPos[ipart]->Draw("P");
+      else hCorrFactPos[ipart]->Draw("PSAME");
+      hCorrFactNeg[ipart]->Draw("PSAME");
+    }
+    cCorr->BuildLegend();
+    
+    TCanvas *cCorrMean=new TCanvas("cCorrMean","cCorrMean");
+    cCorrMean->SetGridy();
+    for(Int_t ipart=0;ipart<=2;ipart++){
+      if(ipart==0){
+       TH1F * hempty=(TH1F*)hCorrFactPosMean[ipart]->Clone();
+       hempty->SetXTitle("P_{T} (GeV/c)");
+       hempty->SetYTitle("Acc x #epsilon x Contamination");
+       hempty->Reset("all");
+       hempty->Draw();
+       
+      }
+      hCorrFactPosMean[ipart]->Draw("PSAME");
+      hCorrFactNegMean[ipart]->Draw("PSAME");
+    }
+    cCorrMean->BuildLegend();
+    
+    TCanvas *cCorrcfPos=new TCanvas("cCorrcfPos","cCorrcfPos");
+    cCorrcfPos->SetGridy();
+    for(Int_t ipart=0;ipart<=2;ipart++){
+      if(ipart==0)hCorrFactPos[ipart]->Draw("P");
+      else hCorrFactPos[ipart]->Draw("PSAME");
+      hCorrFactPosMean[ipart]->Draw("PSAME");
+    }
+    cCorrcfPos->BuildLegend();
+    
+    TCanvas *cCorrcfNeg=new TCanvas("cCorrcfNeg","cCorrcfNeg");
+    cCorrcfNeg->SetGridy();
+    for(Int_t ipart=0;ipart<=2;ipart++){
+      if(ipart==0)hCorrFactNeg[ipart]->Draw("P");
+      else hCorrFactNeg[ipart]->Draw("PSAME");
+      hCorrFactNegMean[ipart]->Draw("PSAME");
+    }
+    cCorrcfNeg->BuildLegend();
+    
+    TCanvas *cEff=new TCanvas("cEff","cEff");
+    cEff->SetGridy();
+    for(Int_t ipart=0;ipart<=2;ipart++){
+      if(ipart==0)hEffPos[ipart]->Draw("P");
+      else hEffPos[ipart]->Draw("PSAME");
+      hEffNeg[ipart]->Draw("PSAME");
+    }
+    cEff->BuildLegend();
+    
+    TCanvas *cEffMean=new TCanvas("cEffMean","cEffMean");
+    cEffMean->SetGridy();
+    for(Int_t ipart=0;ipart<=2;ipart++){
+      if(ipart==0)hEffPosMean[ipart]->Draw("P");
+      else hEffPosMean[ipart]->Draw("PSAME");
+      hEffNegMean[ipart]->Draw("PSAME");
+    }
+    cEffMean->BuildLegend();
+    
+    RatioPlot(hEffNegMean,hEffPosMean,"Neg","Pos","Efficiency",Low,Up);
+    
+    
+    TCanvas *cEffcfPos=new TCanvas("cEffcfPos","cEffcfPos");
+    cEffcfPos->SetGridy();
+    for(Int_t ipart=0;ipart<=2;ipart++){
+      if(ipart==0)hEffPos[ipart]->Draw("P");
+      else hEffPos[ipart]->Draw("PSAME");
+      hEffPosMean[ipart]->Draw("PSAME");
+    }
+    cEffcfPos->BuildLegend();
+    
+    TCanvas *cEffcfNeg=new TCanvas("cEffcfNeg","cEffcfNeg");
+    cEffcfNeg->SetGridy();
+    for(Int_t ipart=0;ipart<=2;ipart++){
+      if(ipart==0)hEffNeg[ipart]->Draw("P");
+      else hEffNeg[ipart]->Draw("PSAME");
+      hEffNegMean[ipart]->Draw("PSAME");
+    }
+    cEffcfNeg->BuildLegend();
+  }
+    
+    
+    
+    
+  for(Int_t ipart=0;ipart<=2;ipart++) {
+    hHistPosNSigma[ipart]->Scale(normMC,"width");
+    hHistNegNSigma[ipart]->Scale(normMC,"width");
+    hHistPosNSigmaMean[ipart]->Scale(normMC,"width");
+    hHistNegNSigmaMean[ipart]->Scale(normMC,"width");
+  }
+    
+
+  //////////////////////////////Prim/All Sec/All SecStr/All
+  TH1F *hPrimPos[3];
+  TH1F *hPrimNeg[3];
+  TH1F *hSecPos[3];
+  TH1F *hSecNeg[3];
+  TH1F *hSecStrPos[3];
+  TH1F *hSecStrNeg[3];
+  TH1F *hAllPos[3];
+  TH1F *hAllNeg[3];
+  TH1F *hAllSecPos[3];
+  TH1F *hAllSecNeg[3];
+    
+  for(Int_t ipart=0;ipart<=2;ipart++){
+    hPrimPos[ipart]=(TH1F*)cOutputMC->FindObject(Form("fHistPrimMCposReco%i",ipart));
+    hPrimNeg[ipart]=(TH1F*)cOutputMC->FindObject(Form("fHistPrimMCnegReco%i",ipart));
+    hSecPos[ipart]=(TH1F*)cOutputMC->FindObject(Form("fHistSecMatMCposReco%i",ipart));
+    hSecNeg[ipart]=(TH1F*)cOutputMC->FindObject(Form("fHistSecMatMCnegReco%i",ipart));
+    hSecStrPos[ipart]=(TH1F*)cOutputMC->FindObject(Form("fHistSecStrMCposReco%i",ipart));
+    hSecStrNeg[ipart]=(TH1F*)cOutputMC->FindObject(Form("fHistSecStrMCnegReco%i",ipart));
+    
+    hAllSecPos[ipart]=(TH1F*)hSecPos[ipart]->Clone(Form("hAllSecPos%i",ipart));
+    hAllSecPos[ipart]->Add(hSecStrPos[ipart]);
+    hAllSecNeg[ipart]=(TH1F*)hSecNeg[ipart]->Clone(Form("hAllSecNeg%i",ipart));
+    hAllSecNeg[ipart]->Add(hSecStrNeg[ipart]);
+    hAllPos[ipart]=(TH1F*)hPrimPos[ipart]->Clone(Form("hAllPos%i",ipart));
+    hAllPos[ipart]->Add(hSecPos[ipart]);
+    hAllPos[ipart]->Add(hSecStrPos[ipart]);
+    hAllNeg[ipart]=(TH1F*)hPrimNeg[ipart]->Clone(Form("hAllNeg%i",ipart));
+    hAllNeg[ipart]->Add(hSecNeg[ipart]);
+    hAllNeg[ipart]->Add(hSecStrNeg[ipart]);
+    
+    hPrimPos[ipart]->Divide(hAllPos[ipart]);
+    hSecPos[ipart]->Divide(hAllPos[ipart]);
+    hSecStrPos[ipart]->Divide(hAllPos[ipart]);
+    hPrimNeg[ipart]->Divide(hAllNeg[ipart]);
+    hSecNeg[ipart]->Divide(hAllNeg[ipart]);
+    hSecStrNeg[ipart]->Divide(hAllNeg[ipart]);
+    
+    hAllSecPos[ipart]->Divide(hAllPos[ipart]);
+    hAllSecNeg[ipart]->Divide(hAllNeg[ipart]);
+    for(Int_t binnum=1;binnum<=hAllSecNeg[ipart]->GetNbinsX();binnum++){
+      hAllSecPos[ipart]->SetBinContent(binnum,1-hAllSecPos[ipart]->GetBinContent(binnum));
+      hAllSecNeg[ipart]->SetBinContent(binnum,1-hAllSecNeg[ipart]->GetBinContent(binnum));
+    }
+  }
+  
+  /////////////////////////////////////////////
+  TFile *outMC2=new TFile(foutMC.Data(),"update");
+  TList *lMC=new TList();
+  lMC->SetOwner();
+  lMC->SetName(Form("MC_Mult%ito%i",lowmult,upmult));
+  for(Int_t ipart=0;ipart<=2;ipart++){
+    lMC->Add(hCorrFactPos[ipart]);
+    lMC->Add(hCorrFactNeg[ipart]);
+    lMC->Add(hCorrFactPosMean[ipart]);
+    lMC->Add(hCorrFactNegMean[ipart]);
+    lMC->Add(hEffPos[ipart]);
+    lMC->Add(hEffNeg[ipart]);
+    lMC->Add(hEffPosMean[ipart]);
+    lMC->Add(hEffNegMean[ipart]);
+    lMC->Add(hHistPosNSigma[ipart]);
+    lMC->Add(hHistNegNSigma[ipart]);
+    lMC->Add(hHistPosNSigmaMean[ipart]);
+    lMC->Add(hHistNegNSigmaMean[ipart]);
+    lMC->Add(hPrimPos[ipart]);
+    lMC->Add(hPrimNeg[ipart]);
+    lMC->Add(hSecPos[ipart]);
+    lMC->Add(hSecNeg[ipart]);
+    lMC->Add(hSecStrPos[ipart]);
+    lMC->Add(hSecStrNeg[ipart]);
+    lMC->Add(hAllSecPos[ipart]);
+    lMC->Add(hAllSecNeg[ipart]);
+  }
+  lMC->Write(Form("MC_Mult%ito%i",lowmult,upmult),1);
+  outMC2->Close();
+  delete outMC2;
+  
+    
+  TString lnameDATA=Form("clistITSsaMult%ito%i",lowmult,upmult);
+  Printf("\n\n------------------ READING DATA  %s -------------------\n",lnameDATA.Data());
+  TFile *finDATA = new TFile(fnameDATA.Data());
+  TDirectory *dirFileDATA=(TDirectory*)finDATA->Get("PWG2SpectraITSsa");
+  TList *cOutputDATA = (TList*)dirFileDATA->Get(lnameDATA.Data());
+  
+  
+  for(Int_t ipart=0;ipart<=2;ipart++){
+    hHistPosNSigmaDATA[ipart]=(TH1F*)cOutputDATA->FindObject(Form("hHistPosNSigma%i",ipart));
+    hHistNegNSigmaDATA[ipart]=(TH1F*)cOutputDATA->FindObject(Form("hHistNegNSigma%i",ipart));
+    hHistPosNSigmaMeanDATA[ipart]=(TH1F*)cOutputDATA->FindObject(Form("hHistPosNSigmaMean%i",ipart));
+    hHistNegNSigmaMeanDATA[ipart]=(TH1F*)cOutputDATA->FindObject(Form("hHistNegNSigmaMean%i",ipart)); 
+    
+    hITSsaPos[ipart]=(TH1F*)hHistPosNSigmaDATA[ipart]->Clone(Form("hITSsaPos%i",ipart));
+    hITSsaNeg[ipart]=(TH1F*)hHistNegNSigmaDATA[ipart]->Clone(Form("hITSsaNeg%i",ipart));
+    hITSsaPosMean[ipart]=(TH1F*)hHistPosNSigmaMeanDATA[ipart]->Clone(Form("hITSsaPosAsymm%i",ipart));
+    hITSsaNegMean[ipart]=(TH1F*)hHistNegNSigmaMeanDATA[ipart]->Clone(Form("hITSsaNegAsymm%i",ipart));
+    hITSsaRawPos[ipart]=(TH1F*)hHistPosNSigmaDATA[ipart]->Clone(Form("hITSsaRawPos%i",ipart));
+    hITSsaRawNeg[ipart]=(TH1F*)hHistNegNSigmaDATA[ipart]->Clone(Form("hITSsaRawNeg%i",ipart));
+    hITSsaRawPosMean[ipart]=(TH1F*)hHistPosNSigmaMeanDATA[ipart]->Clone(Form("hITSsaRawPosAsymm%i",ipart));
+    hITSsaRawNegMean[ipart]=(TH1F*)hHistNegNSigmaMeanDATA[ipart]->Clone(Form("hITSsaRawNegAsymm%i",ipart));
+    
+    hITSsaPos[ipart]->Sumw2();
+    hITSsaNeg[ipart]->Sumw2();
+    hITSsaPosMean[ipart]->Sumw2();
+    hITSsaNegMean[ipart]->Sumw2();
+    hITSsaRawPos[ipart]->Sumw2();
+    hITSsaRawNeg[ipart]->Sumw2();
+    hITSsaRawPosMean[ipart]->Sumw2();
+    hITSsaRawNegMean[ipart]->Sumw2();
+    
+    hITSsaPos[ipart]->Divide(hCorrFactPos[ipart]);
+    hITSsaNeg[ipart]->Divide(hCorrFactNeg[ipart]);
+    hITSsaPosMean[ipart]->Divide(hCorrFactPosMean[ipart]);
+    hITSsaNegMean[ipart]->Divide(hCorrFactNegMean[ipart]);
+    
+    SetDrawAtt(21+ipart,ipart+4,1,ipart+4,1,hITSsaPos[ipart]);
+    SetDrawAtt(25+ipart,ipart+4,1,ipart+4,1,hITSsaNeg[ipart]);
+    SetDrawAtt(21+ipart,ipart+1,1,ipart+1,1,hITSsaPosMean[ipart]);
+    SetDrawAtt(25+ipart,ipart+1,1,ipart+1,1,hITSsaNegMean[ipart]); 
+    SetDrawAtt(21+ipart,ipart+4,1,ipart+4,1,hITSsaRawPos[ipart]);
+    SetDrawAtt(25+ipart,ipart+4,1,ipart+4,1,hITSsaRawNeg[ipart]);
+    SetDrawAtt(21+ipart,ipart+1,1,ipart+1,1,hITSsaRawPosMean[ipart]);
+    SetDrawAtt(25+ipart,ipart+1,1,ipart+1,1,hITSsaRawNegMean[ipart]); 
+  }
+  
+  if(!CorrectOnlyForEfficiency){
+    
+    
+    if(DCACorrection){
+      printf("- Correction from DCA\n");
+      
+      // Correction factor based on fit to DCA distr for P and Pibar
+      TString fnameDCAPosP=Form("./DCACorr%s_%s_%.0fDCA_PosP_TFraction.root",period.Data(),MCname.Data(),DCAcut);
+      TFile *fDCAPosP= new TFile(fnameDCAPosP.Data());
+      TH1F * hDCAPosP= (TH1F*)fDCAPosP->Get("hPrimAllDATAMC");
+      //TH1F * hDCAPosP= (TH1F*)fDCAPosP->Get("hPrimAllDATA");
+      TString fnameDCANegP=Form("./DCACorr%s_%s_%.0fDCA_NegP_TFraction.root",period.Data(),MCname.Data(),DCAcut);
+      TFile *fDCANegP= new TFile(fnameDCANegP.Data());
+      TH1F * hDCANegP= (TH1F*)fDCANegP->Get("hPrimAllDATAMC");
+      //TH1F * hDCANegP= (TH1F*)fDCANegP->Get("hPrimAllDATA");
+      
+      for(Int_t pbin=0;pbin<=hDCAPosP->GetNbinsX();pbin++)hDCAPosP->SetBinError(pbin,0);
+      for(Int_t pbin=0;pbin<=hDCANegP->GetNbinsX();pbin++)hDCANegP->SetBinError(pbin,0);
+      
+      hITSsaPos[2]->Multiply(hDCAPosP);
+      hITSsaNeg[2]->Multiply(hDCANegP);
+      hITSsaPosMean[2]->Multiply(hDCAPosP);
+      hITSsaNegMean[2]->Multiply(hDCANegP);
+      
+      hDCAPosP->Draw("");
+      hDCANegP->Draw("same"); 
+
+    }
+    else{
+      Printf("NO DCA correction for protons");
+    }
+    
+    if(GFCorrection){
+      printf("- Geant/FLUKA correction\n"); 
+      //GeantFluka correction for Pi and K
+      TString fnameGeanFlukaPi="$ALICE_ROOT/PWGLF/SPECTRA/PiKaPr/ITSsa/RootFilesGeantFlukaCorrection/correctionForCrossSection.211.root";
+      TFile *fGeanFlukaPi= new TFile(fnameGeanFlukaPi.Data());
+      TH1F *hGeantFlukaPiPos=(TH1F*)fGeanFlukaPi->Get("gHistCorrectionForCrossSectionParticles");
+      TH1F *hGeantFlukaPiNeg=(TH1F*)fGeanFlukaPi->Get("gHistCorrectionForCrossSectionAntiParticles");
+      for(Int_t binPi=0;binPi<=hITSsaPos[0]->GetNbinsX();binPi++){
+       Float_t FlukaCorrPiPos=hGeantFlukaPiPos->GetBinContent(hGeantFlukaPiPos->FindBin(hITSsaPos[0]->GetBinCenter(binPi)));
+       Float_t FlukaCorrPiNeg=hGeantFlukaPiNeg->GetBinContent(hGeantFlukaPiNeg->FindBin(hITSsaNeg[0]->GetBinCenter(binPi)));
+       hITSsaPos[0]->SetBinContent(binPi,hITSsaPos[0]->GetBinContent(binPi)*FlukaCorrPiPos);
+       hITSsaPos[0]->SetBinError(binPi,hITSsaPos[0]->GetBinError(binPi)*FlukaCorrPiPos);
+       hITSsaNeg[0]->SetBinContent(binPi,hITSsaNeg[0]->GetBinContent(binPi)*FlukaCorrPiNeg);
+       hITSsaNeg[0]->SetBinError(binPi,hITSsaNeg[0]->GetBinError(binPi)*FlukaCorrPiNeg);
+       hITSsaPosMean[0]->SetBinContent(binPi,hITSsaPosMean[0]->GetBinContent(binPi)*FlukaCorrPiPos);
+       hITSsaPosMean[0]->SetBinError(binPi,hITSsaPosMean[0]->GetBinError(binPi)*FlukaCorrPiPos);
+       hITSsaNegMean[0]->SetBinContent(binPi,hITSsaNegMean[0]->GetBinContent(binPi)*FlukaCorrPiNeg);
+       hITSsaNegMean[0]->SetBinError(binPi,hITSsaNegMean[0]->GetBinError(binPi)*FlukaCorrPiNeg);
+      }
+      
+      TString fnameGeanFlukaK="$ALICE_ROOT/PWGLF/SPECTRA/PiKaPr/ITSsa/RootFilesGeantFlukaCorrection/correctionForCrossSection.321.root";
+      TFile *fGeanFlukaK= new TFile(fnameGeanFlukaK.Data());
+      TH1F *hGeantFlukaKPos=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionParticles");
+      TH1F *hGeantFlukaKNeg=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionAntiParticles");
+      for(Int_t binK=0;binK<=hITSsaPos[1]->GetNbinsX();binK++){
+       Float_t FlukaCorrKPos=hGeantFlukaKPos->GetBinContent(hGeantFlukaKPos->FindBin(hITSsaPos[1]->GetBinCenter(binK)));
+       Float_t FlukaCorrKNeg=hGeantFlukaKNeg->GetBinContent(hGeantFlukaKNeg->FindBin(hITSsaNeg[1]->GetBinCenter(binK)));
+       hITSsaPos[1]->SetBinContent(binK,hITSsaPos[1]->GetBinContent(binK)*FlukaCorrKPos);
+       hITSsaNeg[1]->SetBinContent(binK,hITSsaNeg[1]->GetBinContent(binK)*FlukaCorrKNeg);
+       hITSsaPos[1]->SetBinError(binK,hITSsaPos[1]->GetBinError(binK)*FlukaCorrKPos);
+       hITSsaNeg[1]->SetBinError(binK,hITSsaNeg[1]->GetBinError(binK)*FlukaCorrKNeg);
+       hITSsaPosMean[1]->SetBinContent(binK,hITSsaPosMean[1]->GetBinContent(binK)*FlukaCorrKPos);
+       hITSsaNegMean[1]->SetBinContent(binK,hITSsaNegMean[1]->GetBinContent(binK)*FlukaCorrKNeg);
+       hITSsaPosMean[1]->SetBinError(binK,hITSsaPosMean[1]->GetBinError(binK)*FlukaCorrKPos);
+       hITSsaNegMean[1]->SetBinError(binK,hITSsaNegMean[1]->GetBinError(binK)*FlukaCorrKNeg);
+      }
+      
+      //Geant Fluka for P
+      // ITS specific file for protons/antiprotons
+      Int_t kNCharge=2;
+      Int_t kPos=0;
+      Int_t kNeg=1;
+      TFile* fITS = new TFile ("$ALICE_ROOT/PWGLF/SPECTRA/PiKaPr/ITSsa/RootFilesGeantFlukaCorrection/correctionForCrossSectionITS_20100719.root");
+      // TH2D * hCorrFlukaITS[kNCharge];
+      TH2D * hCorrFlukaITS[2];
+      hCorrFlukaITS[kPos] = (TH2D*)fITS->Get("gHistCorrectionForCrossSectionProtons");
+      hCorrFlukaITS[kNeg] = (TH2D*)fITS->Get("gHistCorrectionForCrossSectionAntiProtons");
+      for(Int_t icharge = 0; icharge < kNCharge; icharge++){
+       Int_t nbins = hITSsaPos[2]->GetNbinsX();
+       Int_t nbinsy=hCorrFlukaITS[icharge]->GetNbinsY();
+       //for(Int_t binP=0;binP<=hITSsaPos[2]->GetNbinsX();binP++){
+       for(Int_t ibin = 0; ibin < nbins; ibin++){
+         Float_t pt = hITSsaPos[2]->GetBinCenter(ibin);
+         Float_t minPtCorrection = hCorrFlukaITS[icharge]->GetYaxis()->GetBinLowEdge(1);
+         Float_t maxPtCorrection = hCorrFlukaITS[icharge]->GetYaxis()->GetBinLowEdge(nbinsy+1);
+         if (pt < minPtCorrection) pt = minPtCorrection+0.0001;
+         if (pt > maxPtCorrection) pt = maxPtCorrection;
+         Float_t correction = hCorrFlukaITS[icharge]->GetBinContent(1,hCorrFlukaITS[icharge]->GetYaxis()->FindBin(pt));
+         cout<<correction<<"     charge "<<icharge<<endl;
+         if(icharge==0){
+           if (correction != 0) {// If the bin is empty this is a  0
+             hITSsaPos[2]->SetBinContent(ibin,hITSsaPos[2]->GetBinContent(ibin)*correction);
+             hITSsaPos[2]->SetBinError(ibin,hITSsaPos[2]->GetBinError  (ibin)*correction);
+             hITSsaPosMean[2]->SetBinContent(ibin,hITSsaPosMean[2]->GetBinContent(ibin)*correction);
+             hITSsaPosMean[2]->SetBinError(ibin,hITSsaPosMean[2]->GetBinError  (ibin)*correction);
+             
+           }else if (hITSsaPos[2]->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
+             cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for protons secondaries, ITS, " << endl;
+             cout << " Bin content: " << hITSsaPos[2]->GetBinContent(ibin)  << endl;
+           }
+         }
+         if(icharge==1){
+           if (correction != 0) {// If the bin is empty this is a  0
+             hITSsaNeg[2]->SetBinContent(ibin,hITSsaNeg[2]->GetBinContent(ibin)*correction);
+             hITSsaNeg[2]->SetBinError(ibin,hITSsaNeg[2]->GetBinError  (ibin)*correction);
+             hITSsaNegMean[2]->SetBinContent(ibin,hITSsaNegMean[2]->GetBinContent(ibin)*correction);
+             hITSsaNegMean[2]->SetBinError(ibin,hITSsaNegMean[2]->GetBinError  (ibin)*correction);
+           }else if (hITSsaNeg[2]->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
+             cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for Antiprotons secondaries, ITS, " << endl;
+             cout << " Bin content: " << hITSsaNeg[2]->GetBinContent(ibin)  << endl;
+           }
+         }
+         
+       }       
+      }
+    }else{
+      Printf("SPECTRA ARE NOT CORRECTED FOR GEANT/FLUKA");
+      
+    }
+    
+  }else{
+    Printf("SPECTRA ARE NOT CORRECTED FOR DCA FIT + GEANT/FLUKA");
+  }
+
+  Float_t norm;
+  TH1F *fHistNEvents =0x0;
+  if(MarekNorm){
+    Printf("Normalization to the number of events with good vertex");
+    TFile *finMAREK = new TFile(fnameDATA.Data());
+    TDirectory *dirFileMAREK=(TDirectory*)finMAREK->Get("PWG2SpectraITSTPC");
+    if(!dirFileMAREK)Printf("!dirFileMAREK");
+    TString lnameMAREK=Form("outputlow%dup%dHI0",lowmult,upmult);
+    TList *cOutputMAREK = (TList*)dirFileMAREK->Get(lnameMAREK.Data());
+    if(!cOutputMAREK)Printf("!cOutputMAREK");
+    TH1F *StatsHistMAREK=(TH1F*)cOutputMAREK->FindObject("StatsHist");
+    norm=(Float_t)1/StatsHistMAREK->GetBinContent(3);
+  }
+  else{
+    Printf("Normalization to the number of events after phys sel");
+    fHistNEvents = (TH1F*)cOutputDATA->FindObject("fHistNEvents");
+    norm=(Float_t)1/fHistNEvents->GetBinContent(fHistNEvents->FindBin(1.));
+  }
+  if(NormalizeData)printf("- Normalization (NEvents): %.0f \n",1/norm); 
+  else printf("- Normalization (NEvents): %.0f !!!!!!!SPECTRA ARE NOT NORMALIZED!!!!!!\n",fHistNEvents->GetBinContent(fHistNEvents->FindBin(1.))); 
+  for(Int_t ipart=0;ipart<=2;ipart++){
+    if(NormalizeData){
+      hITSsaPos[ipart]->Scale(norm,"width");
+      hITSsaNeg[ipart]->Scale(norm,"width");
+      hITSsaPosMean[ipart]->Scale(norm,"width");
+      hITSsaNegMean[ipart]->Scale(norm,"width");
+      hITSsaRawPos[ipart]->Scale(norm,"width");
+      hITSsaRawNeg[ipart]->Scale(norm,"width");
+      hITSsaRawPosMean[ipart]->Scale(norm,"width");
+      hITSsaRawNegMean[ipart]->Scale(norm,"width");
+    }
+    else{
+      hITSsaPos[ipart]->Scale(1,"width");
+      hITSsaNeg[ipart]->Scale(1,"width");
+      hITSsaPosMean[ipart]->Scale(1,"width");
+      hITSsaNegMean[ipart]->Scale(1,"width");
+      hITSsaRawPos[ipart]->Scale(1,"width");
+      hITSsaRawNeg[ipart]->Scale(1,"width");
+      hITSsaRawPosMean[ipart]->Scale(1,"width");
+      hITSsaRawNegMean[ipart]->Scale(1,"width");
+    }
+    if(NormalizeToINEL){
+      Float_t INEL=0.86;
+      printf("Normalization to INEL : %.2f \n",INEL);
+      hITSsaPos[ipart]->Scale(INEL,"");
+      hITSsaNeg[ipart]->Scale(INEL,"");
+      hITSsaPosMean[ipart]->Scale(INEL,"");
+      hITSsaNegMean[ipart]->Scale(INEL,"");
+      hITSsaRawPos[ipart]->Scale(INEL,"");
+      hITSsaRawNeg[ipart]->Scale(INEL,"");
+      hITSsaRawPosMean[ipart]->Scale(INEL,"");
+      hITSsaRawNegMean[ipart]->Scale(INEL,"");
+    }
+    else printf("Normalization to INEL ignored \n");
+    
+    hITSsaPos[ipart]->SetMinimum(0.00001);
+    hITSsaNeg[ipart]->SetMinimum(0.00001);
+    hITSsaPosMean[ipart]->SetMinimum(0.00001);
+    hITSsaNegMean[ipart]->SetMinimum(0.00001);
+    hITSsaRawPos[ipart]->SetMinimum(0.00001);
+    hITSsaRawNeg[ipart]->SetMinimum(0.00001);
+    hITSsaRawPosMean[ipart]->SetMinimum(0.00001);
+    hITSsaRawNegMean[ipart]->SetMinimum(0.00001);
+    hITSsaPos[ipart]->GetXaxis()->SetTitle("Pt [GeV/c]");
+    hITSsaNeg[ipart]->GetXaxis()->SetTitle("Pt [GeV/c]");
+    hITSsaPosMean[ipart]->GetXaxis()->SetTitle("Pt [GeV/c]");
+    hITSsaNegMean[ipart]->GetXaxis()->SetTitle("Pt [GeV/c]");
+    hITSsaRawPos[ipart]->GetXaxis()->SetTitle("Pt [GeV/c]");
+    hITSsaRawNeg[ipart]->GetXaxis()->SetTitle("Pt [GeV/c]");
+    hITSsaRawPosMean[ipart]->GetXaxis()->SetTitle("Pt [GeV/c]");
+    hITSsaRawNegMean[ipart]->GetXaxis()->SetTitle("Pt [GeV/c]");
+    hITSsaPos[ipart]->SetTitle(Form("ITSsa NSigma Symm Pos%i",ipart));
+    hITSsaNeg[ipart]->SetTitle(Form("ITSsa NSigma Symm Neg%i",ipart));
+    hITSsaPosMean[ipart]->SetTitle(Form("ITSsa NSigma Asymm Pos%i",ipart));
+    hITSsaNegMean[ipart]->SetTitle(Form("ITSsa NSigma Asymm Neg%i",ipart));
+    hITSsaRawPos[ipart]->SetTitle(Form("ITSsaRaw NSigma Symm Pos%i",ipart));
+    hITSsaRawNeg[ipart]->SetTitle(Form("ITSsaRaw NSigma Symm Neg%i",ipart));
+    hITSsaRawPosMean[ipart]->SetTitle(Form("ITSsaRaw NSigma Asymm Pos%i",ipart));
+    hITSsaRawNegMean[ipart]->SetTitle(Form("ITSsaRaw NSigma Asymm Neg%i",ipart));
+      
+
+    ResetOutOfRange(hITSsaPos[ipart],ipart,lowRange,upRange);
+    ResetOutOfRange(hITSsaNeg[ipart],ipart,lowRange,upRange);
+    ResetOutOfRange(hITSsaPosMean[ipart],ipart,lowRange,upRange);
+    ResetOutOfRange(hITSsaNegMean[ipart],ipart,lowRange,upRange);
+
+  }
+  TH1F *hSystematicsITSsaPos[3];
+  TH1F *hSystematicsITSsaNeg[3];
+  TH1F *hSystematicsITSsaPosMean[3];
+  TH1F *hSystematicsITSsaNegMean[3];
+  
+  if(MakeSystErr){
+    printf("- making Systematic Error\n");
+    //preliminary systematic error
+    TFile *fsys = new TFile("RootFilesSystError/ITSsa-systematics-20101014.root");
+    //fsys->Print("all");
+    //cout<<fsys->GetSize()<<endl;
+    TH1F* hSystPos[3];
+    TH1F* hSystNeg[3];
+    hSystPos[0]=(TH1F*)fsys->Get("hSystTotPosPion");
+    hSystNeg[0]=(TH1F*)fsys->Get("hSystTotNegPion");
+    hSystPos[1]=(TH1F*)fsys->Get("hSystTotPosKaon");
+    hSystNeg[1]=(TH1F*)fsys->Get("hSystTotNegKaon");
+    hSystPos[2]=(TH1F*)fsys->Get("hSystTotPosProton");
+    hSystNeg[2]=(TH1F*)fsys->Get("hSystTotNegProton");
+
+    for(Int_t ipart=0;ipart<=2;ipart++){
+      hSystematicsITSsaPos[ipart]=(TH1F*)hITSsaPos[ipart]->Clone(Form("hSystematicsITSsaPos%i",ipart));  
+      hSystematicsITSsaNeg[ipart]=(TH1F*)hITSsaNeg[ipart]->Clone(Form("hSystematicsITSsaNeg%i",ipart));  
+      hSystematicsITSsaPosMean[ipart]=(TH1F*)hITSsaPosMean[ipart]->Clone(Form("hSystematicsITSsaPosAsymm%i",ipart));  
+      hSystematicsITSsaNegMean[ipart]=(TH1F*)hITSsaNegMean[ipart]->Clone(Form("hSystematicsITSsaNegAsymm%i",ipart));  
+    
+      for(Int_t ibin=0;ibin<=hSystematicsITSsaPos[ipart]->GetNbinsX();ibin++){
+       Float_t syserr=hSystPos[ipart]->GetBinContent(hSystPos[ipart]->FindBin(hSystematicsITSsaPos[ipart]->GetBinCenter(ibin)));
+       hSystematicsITSsaPos[ipart]->SetBinError(ibin,syserr*hSystematicsITSsaPos[ipart]->GetBinContent(ibin));
+       hSystematicsITSsaPosMean[ipart]->SetBinError(ibin,syserr*hSystematicsITSsaPosMean[ipart]->GetBinContent(ibin));
+      }
+      for(Int_t ibin=0;ibin<=hSystematicsITSsaNeg[ipart]->GetNbinsX();ibin++){
+       Float_t syserr=hSystNeg[ipart]->GetBinContent(hSystNeg[ipart]->FindBin(hSystematicsITSsaNeg[ipart]->GetBinCenter(ibin)));
+       hSystematicsITSsaNeg[ipart]->SetBinError(ibin,syserr*hSystematicsITSsaNeg[ipart]->GetBinContent(ibin));
+       hSystematicsITSsaNegMean[ipart]->SetBinError(ibin,syserr*hSystematicsITSsaNegMean[ipart]->GetBinContent(ibin));
+      }
+    }
+  }
+    
+  TCanvas *cSpectraSymm=new TCanvas("cSpectraSymm","cSpectraSymm");
+  cSpectraSymm->SetGridy();
+  for(Int_t ipart=0;ipart<=2;ipart++){
+    if(ipart==0)hITSsaPos[ipart]->DrawCopy("P");
+    else hITSsaPos[ipart]->DrawCopy("PSAME");
+    hITSsaNeg[ipart]->DrawCopy("PSAME");
+  }
+  cSpectraSymm->BuildLegend();
+  
+  TCanvas *cSpectraAsymm=new TCanvas("cSpectraAsymm","cSpectraAsymm");
+  cSpectraAsymm->SetGridy();
+  for(Int_t ipart=0;ipart<=2;ipart++){
+    if(ipart==0)hITSsaPosMean[ipart]->DrawCopy("P");
+    else hITSsaPosMean[ipart]->DrawCopy("PSAME");
+    hITSsaNegMean[ipart]->DrawCopy("PSAME");
+  }
+  cSpectraAsymm->BuildLegend();
+  
+
+  
+  ///////////////////////////  SCALING SPECTRA FOR TRACKING EFFICIENCY   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  Float_t ScalingFactor=1.015;
+  Printf("SCALING FACTOR: %f",ScalingFactor);
+  for(Int_t ipart=0;ipart<=2;ipart++){
+    hITSsaPosMean[ipart]->Scale(ScalingFactor);
+    Printf("Scaling %s by %f",hITSsaPosMean[ipart]->GetName(),ScalingFactor);
+    hITSsaNegMean[ipart]->Scale(ScalingFactor);
+    Printf("Scaling %s by %f",hITSsaNegMean[ipart]->GetName(),ScalingFactor);
+    hITSsaPos[ipart]->Scale(ScalingFactor);
+    Printf("Scaling %s by %f",hITSsaPos[ipart]->GetName(),ScalingFactor);
+    hITSsaNeg[ipart]->Scale(ScalingFactor);
+    Printf("Scaling %s by %f",hITSsaNeg[ipart]->GetName(),ScalingFactor);
+  }
+  
+  TFile *outDATA2=new TFile(foutDATA.Data(),"update");
+  TList *lDATA=new TList();
+  lDATA->SetOwner();
+  lDATA->SetName(Form("DATA_Mult%ito%i",lowmult,upmult));
+  if(SaveOnlyAsymm){
+    for(Int_t ipart=0;ipart<=2;ipart++)  lDATA->Add(hITSsaPosMean[ipart]);
+    for(Int_t ipart=0;ipart<=2;ipart++)  lDATA->Add(hITSsaNegMean[ipart]);
+  }else{
+    for(Int_t ipart=0;ipart<=2;ipart++){
+      lDATA->Add(fHistPrimMCposBefEvSel[ipart]);
+      lDATA->Add(fHistPrimMCnegBefEvSel[ipart]);
+      lDATA->Add(hHistPosNSigmaPrimMean[ipart]);
+      lDATA->Add(hHistNegNSigmaPrimMean[ipart]);
+      lDATA->Add(hHistPosNSigmaPrimMCMean[ipart]);
+      lDATA->Add(hHistNegNSigmaPrimMCMean[ipart]);
+      lDATA->Add(hITSsaPos[ipart]);
+      lDATA->Add(hITSsaNeg[ipart]);
+      lDATA->Add(hITSsaPosMean[ipart]);
+      lDATA->Add(hITSsaNegMean[ipart]);
+      lDATA->Add(hITSsaRawPos[ipart]);
+      lDATA->Add(hITSsaRawNeg[ipart]);
+      lDATA->Add(hITSsaRawPosMean[ipart]);
+      lDATA->Add(hITSsaRawNegMean[ipart]);
+    }
+    if(MakeSystErr){
+      for(Int_t ipart=0;ipart<=2;ipart++){
+       lDATA->Add(hSystematicsITSsaPos[ipart]);
+       lDATA->Add(hSystematicsITSsaNeg[ipart]);
+       lDATA->Add( hSystematicsITSsaPosMean[ipart]);
+       lDATA->Add( hSystematicsITSsaNegMean[ipart]);
+      }
+    }
+  }
+  lDATA->Write(Form("DATA_Mult%ito%i",lowmult,upmult),1);
+  outDATA2->Close();
+  delete outDATA2;
+  
+  
+  RatioPlot(hITSsaNeg,hITSsaPos,"Neg","Pos","Symm NSigma",Low,Up);
+  RatioPlot(hITSsaNegMean,hITSsaPosMean,"Neg","Pos","Assymm NSigma",Low,Up);
+  RatioPlot(hITSsaPos,hITSsaPosMean,"Symm","Asymm","Positive",Low,Up);
+  RatioPlot(hITSsaNeg,hITSsaNegMean,"Symm","Asymm","Negative",Low,Up);
+  
+  
+} //end main
+
+
+void ResetOutOfRange(TH1F *histo, Int_t ipart, Double_t lowRange[3], Double_t upRange[3]){
+  // set to -1 the bin contents out of the selcted pt range
+  for(Int_t ibin=0;ibin<=histo->GetNbinsX();ibin++){
+    if(histo->GetBinCenter(ibin)<lowRange[ipart] || histo->GetBinCenter(ibin)>upRange[ipart]) histo->SetBinContent(ibin,-1);
+  }
+}
+
+void RatioPlot(TH1F **num,TH1F **den,TString t1,TString t2,TString opt,Double_t Low[],Double_t Up[])
+{
+  gROOT->SetStyle("Plain");
+  gStyle->SetOptTitle(0);
+  TH1F *Num[3];
+  TString title=Form("%s/%s - %s",t1.Data(),t2.Data(),opt.Data());
+  TCanvas *c=new TCanvas(title.Data(),title.Data());
+  c->Divide(3,2);
+  for(Int_t ipart=0;ipart<=2;ipart++){
+    num[ipart]->GetXaxis()->SetRangeUser(Low[ipart],Up[ipart]-0.00001);
+    num[ipart]->SetMarkerColor(2);
+    num[ipart]->SetLineColor(2);
+    num[ipart]->SetStats(kFALSE);
+    den[ipart]->GetXaxis()->SetRangeUser(Low[ipart],Up[ipart]-0.00001);
+    den[ipart]->SetMarkerColor(1);
+    den[ipart]->SetLineColor(1);
+    den[ipart]->SetStats(kFALSE);
+    c->cd(ipart+1);
+    gPad->SetGridy();
+    //gPad->SetLogy();
+    num[ipart]->Draw("");
+    den[ipart]->Draw("same");
+    TPaveText *tpave=new TPaveText(0.1,0.2,0.5,0.29,"brNDC");
+    tpave->SetBorderSize(0);
+    tpave->SetFillStyle(0);
+    tpave->SetFillColor(0);
+    tpave->SetTextColor(2);
+    tpave->SetTextSize(.04);
+    TText *txt1=tpave->AddText(t1.Data());
+    txt1->SetTextFont(62);
+    txt1->SetTextColor(2);
+    TText *txt2=tpave->AddText(t2.Data());
+    txt2->SetTextFont(62);
+    txt2->SetTextColor(1);
+    tpave->Draw();
+    Num[ipart]=(TH1F*)num[ipart]->Clone(title.Data());
+    Num[ipart]->SetTitle(title.Data());
+    Num[ipart]->Divide(den[ipart]);
+    Num[ipart]->GetXaxis()->SetRangeUser(Low[ipart],Up[ipart]-0.00001);
+    Num[ipart]->SetMarkerColor(4);
+    Num[ipart]->SetLineColor(4);
+    Num[ipart]->SetMaximum(1.2);
+    Num[ipart]->SetMinimum(.8);
+    c->cd(ipart+4);
+    gPad->SetGridy();
+    Num[ipart]->SetStats(kFALSE);
+    Num[ipart]->Draw();
+    TPaveText *tpave_1=new TPaveText(0.3,0.7,0.7,0.99,"brNDC");
+    tpave_1->SetBorderSize(0);
+    tpave_1->SetFillStyle(0);
+    tpave_1->SetFillColor(0);
+    tpave_1->SetTextColor(2);
+    tpave_1->SetTextSize(.04);
+    TText *txt1_1=tpave_1->AddText(title.Data());
+    txt1_1->SetTextFont(62);
+    tpave_1->Draw();
+    // TLine line = TLine(Low[ipart],1,Up[ipart],1);
+    // line.SetLineColor(2);
+    // line.SetLineStyle(2);
+    // line.SetLineWidth(2);
+    // line.DrawClone("same");
+  }
+}
+
+
+
+
+ void SetDrawAtt(Int_t markerstyle,Int_t markercolor,Int_t markersize,Int_t linecolor,Int_t linewidth,TH1 *h1)
+{ 
+  h1->SetMarkerStyle(markerstyle);
+  h1->SetMarkerColor(markercolor);
+  h1->SetMarkerSize(markersize);
+  h1->SetLineColor(linecolor);
+  h1->SetLineWidth(linewidth);
+}
+
+