C++ TH1D::Print方法代码示例

void drawHistos(TString rfile) {
  
  TH1::AddDirectory(0);

  TFile f(rfile);
  f.ls();
  TH1D *hpp = (TH1D*)f.Get("hpp");
  TH1D *hpn = (TH1D*)f.Get("hpn");
  TH1D *hnp = (TH1D*)f.Get("hnp");
  TH1D *hnn = (TH1D*)f.Get("hnn");

  hpp->Print();
  hnn->Print();
  hpn->Print();
  hnp->Print();
 
  TCanvas *cv = new TCanvas("cv", "Delta t", 500, 700);
  cv->Divide(1,2);
  
  cv->cd(1);
  drawOverlay(hpp, hnn);
  cv->cd(2);
  drawOverlay(hpn, hnp);

  f.Close();
}

void compareOldandNew()
{
  TCanvas * c1 = new TCanvas("c1","c1",800,600);
  TLatex lat;  

  int centBounds[6] = {60,80,100,120,140,200};
  for(int i = 0 ; i<5; i++){
  const char * histName = "pPbPbp_FF";
  TFile * fold = TFile::Open("FragmentationFunctionsUE2_nominal.root","read");
  TFile * fnew = TFile::Open("FragmentationFunctionsUE2.root","read");
  TH1D * num = (TH1D*)fnew->Get(Form("%s_%d_%d",histName,centBounds[i],centBounds[i+1]));
  TH1D * den = (TH1D*)fold->Get(Form("%s_%d_%d",histName,centBounds[i],centBounds[i+1]));
  //den->Add(num);
  num->Divide(den);
  //num->Scale(2);
  num->GetYaxis()->SetTitle("(no L2 residual)/(Nominal)");
  num->GetYaxis()->SetRangeUser(0.7,1.3);
  num->GetXaxis()->SetTitle("p_{T}");
  num->GetXaxis()->SetRangeUser(0.5,200);
  num->Print("All");
  num->Draw();
  
  lat.DrawLatex(1,1.2,Form("%d < p_{T}^{jet} < %d",centBounds[i],centBounds[i+1]));

  c1->SetLogx();
  c1->SaveAs(Form("diffPlots/noL2Residual_%d_%d.png",centBounds[i],centBounds[i+1]));
  c1->SaveAs(Form("diffPlots/noL2Residual_%d_%d.pdf",centBounds[i],centBounds[i+1]));
  fnew->Close();
  fold->Close();
  }  

  /*const char * histName = "PbPbTrackSpectrum_0_5";
  TFile * fnew = TFile::Open("Spectra_Jun9_noChi2Cut.root","read");
  TFile * fold = TFile::Open("Spectra_Jun9_withChi2Cut.root","read");
  TH1D * num = (TH1D*)fnew->Get(histName);
  TH1D * den = (TH1D*)fold->Get(histName);
  //den->Add(num);
  num->Divide(den);
  //num->Scale(2);
  num->GetYaxis()->SetTitle("(No Chi2)/(With Chi2)");
  num->GetYaxis()->SetRangeUser(0.7,1.5);
  num->GetXaxis()->SetRangeUser(0.7,350);
  num->Print("All");
  num->Draw();
  c1->SetLogx();
  c1->SaveAs("plots/comparisonPlots/Chi2CutTest_PbPb_0_5.png");
  c1->SaveAs("plots/comparisonPlots/Chi2CutTest_PbPb_0_5.pdf");
  */

  //c1->SaveAs("plots/comparisonPlots/ppChargeFraction2.C");
}

void MakeClosurePlot(TH1D* hObs, TH1D* hPred, TString name, bool logy=true) {
  set_style(hObs,"data_obs");
  set_style(hPred,"pred");

  // Setup canvas and pads
  TCanvas * c1 = new TCanvas("c1", "c1", 700, 700);
  c1->SetFillStyle(4000);
  TPad * pad1 = new TPad("pad1", "top pad" , 0.0, 0.3, 1.0, 1.0);
  pad1->SetFillStyle(4000);
  pad1->SetBottomMargin(0.0);
  pad1->Draw();
  TPad * pad2 = new TPad("pad2", "bottom pad", 0.0, 0.0, 1.0, 0.3);
  pad2->SetFillStyle(4000);
  pad2->SetTopMargin(0.0);
  pad2->SetBottomMargin(0.35);
  pad2->Draw();
  pad1->cd();
  pad1->SetLogy(logy);

  TH1D * staterr = (TH1D *) hPred->Clone("staterr");
  //  staterr->Sumw2();
  //staterr->SetFillColor(kRed);
  staterr->SetFillColor(kGray+3);
  staterr->SetMarkerSize(0);
  staterr->SetFillStyle(3013);

  //cout << "Compute ratio hist..." << endl;
  TH1D * ratio = (TH1D *) hPred->Clone("ratio");
  //ratio->Sumw2();
  ratio->SetMarkerSize(0.8);
  ratio->SetMarkerColor(1);
  //ratio->SetMarkerSize(0.5);
  ratio->Divide(hObs, hPred, 1., 1., "B");
  ratio->Print("all");
  TH1D * ratiostaterr = (TH1D *) hObs->Clone("ratiostaterr");
  //  ratiostaterr->Sumw2();
  ratiostaterr->SetStats(0);
  ratiostaterr->SetTitle(hPred->GetTitle());
  ratiostaterr->GetYaxis()->SetTitle("Obs/Pred.");
  ratiostaterr->SetMaximum(2.2);
  ratiostaterr->SetMinimum(0);
  ratiostaterr->SetMarkerSize(0);
  //ratiostaterr->SetFillColor(kRed);
  ratiostaterr->SetFillColor(kGray+3);
  ratiostaterr->SetFillStyle(3013);
  ratiostaterr->GetXaxis()->SetLabelSize(0.2);
  ratiostaterr->GetXaxis()->SetLabelOffset(0.03);
  ratiostaterr->GetXaxis()->SetTitleSize(0.14);
  ratiostaterr->GetXaxis()->SetTitleOffset(1.10);
  ratiostaterr->GetYaxis()->SetLabelSize(0.10);
  ratiostaterr->GetYaxis()->SetTitleSize(0.12);
  ratiostaterr->GetYaxis()->SetTitleOffset(0.6);
  ratiostaterr->GetYaxis()->SetNdivisions(505);
  TLine* ratiounity = new TLine(hPred->GetBinLowEdge(1),1,hPred->GetBinLowEdge(hPred->GetNbinsX()+1),1);
  ratiounity->SetLineStyle(2);
  for (Int_t i = 0; i < hPred->GetNbinsX()+2; i++) {
    // if (hObs->GetBinContent(i)>0.) ratio->SetBinError(i, hObs->GetBinError(i)/hObs->GetBinContent(i)); // just the fractional uncertainty on the observation
    ratiostaterr->SetBinContent(i, 1.0);
    if (hPred->GetBinContent(i) > 1e-6) { //< not empty
      double binerror = hPred->GetBinError(i) / hPred->GetBinContent(i);
      ratiostaterr->SetBinError(i, binerror);
    } else {
      ratiostaterr->SetBinError(i, 999.);
    }
  }
  TH1D * ratiosysterr = (TH1D *) ratiostaterr->Clone("ratiosysterr");
  //  ratiosysterr->Sumw2();
  ratiosysterr->SetMarkerSize(0);
  ratiosysterr->SetFillColor(kYellow-4);
  //ratiosysterr->SetFillStyle(3002);
  ratiosysterr->SetFillStyle(1001);
  for (Int_t i = 0; i < hPred->GetNbinsX()+2; i++) {
    if (hPred->GetBinContent(i) > 1e-6) { //< not empty
      double binerror2 = (pow(hPred->GetBinError(i), 2));
      double binerror = sqrt(binerror2);
      ratiosysterr->SetBinError(i, binerror / hPred->GetBinContent(i));
    }
  }


  double max = hObs->GetMaximum();
  if (hPred->GetMaximum() > max) max = hPred->GetMaximum();

  if (logy) {
    hPred->SetMaximum(max*10);
    hPred->SetMinimum(0.09);
  } else {
    hPred->SetMaximum(max*1.5);
    hPred->SetMinimum(0.);
  }

  hPred->GetYaxis()->SetTitle("Events / bin");
  hPred->Draw("hist");
  hObs->Draw("e1,same");
  staterr->Draw("e2 same");
  hPred->GetXaxis()->SetTitleSize(0.035);
  hPred->GetXaxis()->SetLabelSize(0.035);
  /* hObs->GetYaxis()->SetTitleSize(0.035); */
  /* hObs->GetYaxis()->SetLabelSize(0.035); */

//.........这里部分代码省略.........

//.........这里部分代码省略.........

  GenHist->SetTitle("");
  GenHist->SetMarkerStyle(20);
  GenHist->SetMarkerSize(1.2);
  GenHist->SetLineColor(1);
  GenHist->GetXaxis()->SetTitle(xtitlename);
  GenHist->GetYaxis()->SetTitle(ytitlename);
  GenHist->Scale(lumi/lumi_ref);
  EstHist->Scale(lumi/lumi_ref);
  TH1D * GenHist_Normalize = static_cast<TH1D*>(GenHist->Clone("GenHist_Normalize"));
  GenHist_Normalize->SetMaximum(ymax_top);
  GenHist_Normalize->SetMinimum(ymin_top);
  ex1->Draw();
  //GenHist_Normalize->GetListOfFunctions()->Add(ex1);
  GenHist_Normalize->DrawCopy("e");

  EstHist->SetFillStyle(3144);
  EstHist->SetFillColor(kRed-10);
  EstHist->SetMarkerStyle(20);
  EstHist->SetMarkerSize(0.0001);
  TH1D * EstHist_Normalize = static_cast<TH1D*>(EstHist->Clone("EstHist_Normalize"));
  ex2->Draw();
  //EstHist_Normalize->GetListOfFunctions()->Add(ex2);
  EstHist_Normalize->DrawCopy("e2same");
  //EstHist_Normalize->DrawCopy("esame");

  TH1D *EstHist_Normalize_Clone = (TH1D*)EstHist_Normalize->Clone();
  for(int i=1; i<208; i++) {
    EstHist_Normalize_Clone->SetBinError(i,0);
  }
  EstHist_Normalize_Clone->SetFillColor(kWhite);
  EstHist_Normalize_Clone->Draw("esame");

  GenHist->Print("all");
  EstHist->Print("all");
  
  //
  // Re-draw to have "expectation" on top of "prediction"
  ex1->Draw();
  GenHist_Normalize->DrawCopy("esame");
  //

  TString line = "";
  sprintf(tempname,"%8.1f",lumi);
  line+=tempname;
  line+=" fb^{-1} (13 TeV)";
  
  int iPeriod = 0;    // 1=7TeV, 2=8TeV, 3=7+8TeV, 7=7+8+13TeV, 0=free form (uses lumi_sqrtS)
  int iPos=0;
    
  writeExtraText = true;  
  if(doDataVsMC) extraText   = "        Preliminary";
  else extraText   = "        Simulation";
  //float extraTextFont = 52;  // default is helvetica-italics

  // text sizes and text offsets with respect to the top frame
  // in unit of the top margin size
  //lumiTextSize     = 0.5;
  //float lumiTextOffset   = 0.2;
  //cmsTextSize      = 0.65;
  //float cmsTextOffset    = 0.1;  // only used in outOfFrame version
  
  //relPosX    = 0.045;
  //relPosY    = 0.035;
  //relExtraDY = 1.2;
  

Scale_ByLumi(double newLum=3, int imode=1){   // Scale histogram contents and uncertainty to xxx fb^{-1}

  char tempname[200];
  TH1D * temphist;

  const int sampleSize=11;
  string sample[sampleSize]={"TTbar",
			     "WJet_100_200","WJet_200_400","WJet_400_600","WJet_600_inf",
			     "t_top","tW_top","t_antitop","tW_antitop"};

  /*
  TTbar:            806.1 (pb), 
  Wjets HT=100-200: 1.23 * 1817.0 (pb),
  Wjets HT=200-400: 1.23 * 471.6  (pb),
  Wjets HT=400-600: 1.23 * 55.61  (pb),
  Wjets HT=600-inf: 1.23 * 18.81  (pb),
  Single top T_s:       2   (pb)
  Single top Tbar_s:    1   (pb)
  Single top T_t:     103.4 (pb)
  Single top Tbar_t:   61.6 (pb)
  Single top T_tW:     35   (pb)
  Single top Tbar_tW:  35   (pb)
  */
  
  double sample_xsec[sampleSize]={806.1, 
				  1.23*1817.0, 1.23*471.6, 1.23*55.61, 1.23*18.81,
                                  2.,103.4,35.,1.,61.3,35.};  // in pb
  double sample_nevents[sampleSize]={25348009,
				     5257344,4936077,4640594,4581841,
				     500000,3941400,986100,
				     250000,1999800,971800};  // total number of events in ntuples
  double sample_lumi[sampleSize];

  //
  // Looping over different samples
  //
  for(int filei=0;filei<sampleSize;filei++){

    //#######
    // Prediction 
    //#######

    cout << " \n Working on sample: " << sample[filei] << endl;
    sprintf(tempname,"TauHad2/HadTauEstimation_%s_org.root",sample[filei].c_str());
    TFile * file_in  = new TFile(tempname,"READ");
    if (imode==1)      sprintf(tempname,"TauHad2/HadTauEstimation_%s_LumiScaledTo%.1ffbinv.root",sample[filei].c_str(),newLum);
    else if (imode==0) sprintf(tempname,"TauHad2/HadTauEstimation_%s_%.1ffbinv_OrigRelStatError.root",sample[filei].c_str(),newLum);
    else { printf("This mode is not supported\n"); break; }
    TFile * file_out = new TFile(tempname,"UPDATE");

    //#######
    // QCD_Low
    //#######
   
    temphist = (TH1D*) file_in->Get("QCD_Low")->Clone();

    cout << " QCD_Low: Before fix: \n ############### \n ";
    temphist->GetXaxis()->SetRange(1,1);
    if (filei==0) temphist->Print("all");
    else temphist->Print("range");

    //sprintf(tempname,"allEvents/PreSel/MHT_PreSel_allEvents");
    //double nevents_from_root_file = (* (TH1D*)file_in->Get(tempname)).GetEntries();
    sprintf(tempname,"cutflow_preselection");
    double nevents_from_root_file = (* (TH1D*)file_in->Get(tempname)).GetBinContent(1);
  
    double originalSampleLum = nevents_from_root_file / sample_xsec[filei] / 1000.; // original sample luminosity in terms of fb^{-1}
                                                                                   // 1/1000. to convert pb^{-1} to fb^{-1}
    printf("Sample name = %10s, sample lumi = %8.3f, total events = %8.0f, events from allEvent folder = %8.0f\n",
	   sample[filei].c_str(),originalSampleLum,sample_nevents[filei],nevents_from_root_file
	   );

    for (int ibin=0; ibin<temphist->GetNbinsX()+2; ibin++){ // scan including underflow and overflow bins
      double bin=temphist->GetBinContent(ibin);
      temphist->SetBinContent(ibin, bin*newLum/originalSampleLum );
      double binEr=temphist->GetBinError(ibin);
      if (imode==1)      temphist->SetBinError(ibin,   binEr*sqrt(newLum/originalSampleLum) );
      else if (imode==0) temphist->SetBinError(ibin,   binEr*newLum/originalSampleLum );
    }

    cout << " QCD_Low: After fix: \n ############### \n ";
    if (filei==0) temphist->Print("all");
    else temphist->Print("range");
    temphist->GetXaxis()->SetRange();

    temphist->Write();

    //#######
    // QCD_Up
    //#######

    temphist = (TH1D*) file_in->Get("QCD_Up")->Clone();

    cout << " QCD_Up: Before fix: \n ############### \n ";
    temphist->GetXaxis()->SetRange(1,1);
    temphist->Print("range");

    double originalSampleLum = sample_nevents[filei] / sample_xsec[filei] / 1000.; // original sample luminosity in terms of fb^{-1}

    for (int ibin=0; ibin<temphist->GetNbinsX()+2; ibin++){ // scan including underflow and overflow bins
//.........这里部分代码省略.........

// ============================================================================
/// finalize the algorithm
// ============================================================================
StatusCode Aida2Root::finalize()
{

  always() << "Get the native ROOT representation of histograms!" << endmsg ;

  {  // loop over all 1D-histograms
    for ( List::const_iterator ipath = m_1Ds.begin() ;
          m_1Ds.end() != ipath ; ++ipath )
    {
      /// retrieve the historam by full path:
      AIDA::IHistogram1D* aida = 0 ;
      StatusCode sc = histoSvc()->retrieveObject( *ipath , aida ) ;
      if ( sc.isFailure() || 0 == aida )
      { return Error ( "Unable to retrieve 1D-histogram '" + (*ipath) + "'"  ) ; }
      /// convert it to ROOT
      TH1D* root = Gaudi::Utils::Aida2ROOT::aida2root ( aida ) ;
      if ( 0 == root )
      { return Error ( "Unable to convert to ROOT the 1D-histogram '"+(*ipath)+"'") ; }
      /// use the native printout from ROOT
      info() << "The native ROOT printout for 1D-histogram '" << (*ipath) << "':" << endmsg ;
      root->Print() ;

      info () << " |  Compare       | AIDA/HistoStats |     ROOT/TH1    |      Delta      | "  << endmsg ;
      const std::string format = " | %1$-14.14s | %2$ 15.8g | %3$- 15.8g | %4$= 15.8g | "  ;
      info () << print
        ( Gaudi::Utils::HistoStats::mean        ( aida ) ,
          root->GetMean      ()    , "'mean'"        , format  ) << endmsg ;
      info () << print
        ( Gaudi::Utils::HistoStats::meanErr     ( aida ) ,
          root->GetMeanError ()    , "'meanErr'"     , format  ) << endmsg ;
      info () << print
        ( Gaudi::Utils::HistoStats::rms         ( aida ) ,
          root->GetRMS       ()    , "'rms'"         , format  ) << endmsg ;
      info () << print
        ( Gaudi::Utils::HistoStats::rmsErr      ( aida ) ,
          root->GetRMSError  ()    , "'rmsErr'"      , format  ) << endmsg ;
      info () << print
        ( Gaudi::Utils::HistoStats::skewness    ( aida ) ,
          root->GetSkewness ()     , "'skewness'"    , format  ) << endmsg ;
      info () << print
        ( Gaudi::Utils::HistoStats::skewnessErr ( aida ) ,
          root->GetSkewness ( 11 ) , "'skewnessErr'" , format  ) << endmsg ;
      info () << print
        ( Gaudi::Utils::HistoStats::kurtosis  ( aida ) ,
          root->GetKurtosis ()     , "'kurtosis'"    , format  ) << endmsg ;
      info () << print
        ( Gaudi::Utils::HistoStats::kurtosisErr ( aida ) ,
          root->GetKurtosis ( 11 ) , "'kurtosisErr'" , format  ) << endmsg ;
    }
  }

  { // loop over all 2D-histograms
    for ( List::const_iterator ipath = m_2Ds.begin() ;
          m_2Ds.end() != ipath ; ++ipath )
    {
      /// retrieve the historam by full path:
      AIDA::IHistogram2D* aida = 0 ;
      StatusCode sc = histoSvc()->retrieveObject( *ipath , aida ) ;
      if ( sc.isFailure() || 0 == aida )
      { return Error ( "Unable to retrieve 2D-histogram '" + (*ipath) + "'"  ) ; }
      /// convert it to ROOT
      TH2D* root = Gaudi::Utils::Aida2ROOT::aida2root ( aida ) ;
      if ( 0 == root )
      { return Error ( "Unable to convert to ROOT the 2D-histogram '"+(*ipath)+"'") ; }
      /// use the native printout from ROOT
      info() << "The native ROOT printout for 2D-histogram '" << (*ipath) << "':" << endmsg ;
      root->Print() ;
    }
  }

  { // loop over all 3D-histograms
    for ( List::const_iterator ipath = m_3Ds.begin() ;
          m_3Ds.end() != ipath ; ++ipath )
    {
      /// retrieve the historam by full path:
      AIDA::IHistogram3D* aida = 0 ;
      StatusCode sc = histoSvc()->retrieveObject( *ipath , aida ) ;
      if ( sc.isFailure() || 0 == aida )
      { return Error ( "Unable to retrieve 3D-histogram '" + (*ipath) + "'"  ) ; }
      /// convert it to ROOT
	TH3D* root = Gaudi::Utils::Aida2ROOT::aida2root ( aida ) ;
      if ( 0 == root )
      { return Error ( "Unable to convert to ROOT the 3D-histogram '"+(*ipath)+"'") ; }
      /// use the native printout from ROOT
      info() << "The native ROOT printout for 3D-histogram '" << (*ipath) << "':" << endmsg ;
      root->Print() ;
    }
  }


  { // loop over all 1D-profiles
    for ( List::const_iterator ipath = m_1Ps.begin() ;
          m_1Ps.end() != ipath ; ++ipath )
    {
      /// retrieve the historam by full path:
      AIDA::IProfile1D* aida = 0 ;
      StatusCode sc = histoSvc()->retrieveObject( *ipath , aida ) ;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  leg->AddEntry(s.ppUsedByTrigger[3],"Jet80 trigger","f");
  leg->AddEntry((TObject*)0,"|#eta|<1","");
  leg->Draw("same");
  
  
 
  c1->SaveAs("plots/png/ppTrack_FullSpectrum.png");
  c1->SaveAs("plots/pdf/ppTrack_FullSpectrum.pdf");
  c1->SetLogy(0);
  for(int i = 0; i<s.nTriggers; i++) s.ppUsedByTrigger[i]->Divide(s.pp);
  s.ppUsedByTrigger[0]->GetYaxis()->SetRangeUser(0,2);
  s.ppUsedByTrigger[0]->GetYaxis()->SetTitle("Relative Contribution to Bin");
  s.ppUsedByTrigger[0]->Draw("HIST");
  for(int i = 1; i<s.nTriggers; i++) s.ppUsedByTrigger[i]->Draw("HIST same");
  s.ppUsedByTrigger[0]->Draw("sameaxis");
  leg->Draw("same");
  c1->SaveAs("plots/png/ppTrack_FullSpectrum_relativeContribution.png");
  c1->SaveAs("plots/pdf/ppTrack_FullSpectrum_relativeContribution.pdf");
  c1->SetLogy();


  TH1D * jtVsTrk = (TH1D*)s.pp->Clone("jtVsTrkRatio");
  
  //xt plot
  TFile * outF = TFile::Open("Spectra.root","UPDATE");
  TH1D * pp_perMBTrigger_xt = (TH1D*)s.pp_perMBTrigger->Clone("pp_perMBTrigger_xt");
  pp_perMBTrigger_xt->Scale(TMath::Power(5020,4.9));
  TH1D * pp_xt = new TH1D("ppTrackSpectrum_xt",";x_{T};#sqrt{s}^{4.9}E#frac{d^{3}#sigma}{d^{3}p} (mb/GeV^{2})",s.ntrkBins,s.xt_xtrkbins);
  for(int i = 1; i< pp_xt->GetSize()-1; i++){
    pp_xt->SetBinContent(i,pp_perMBTrigger_xt->GetBinContent(i));
    pp_xt->SetBinError(i,pp_perMBTrigger_xt->GetBinError(i));
  }
  pp_xt->Draw();
  pp_xt->Print("All");
  pp_xt->Write();
  c1->SaveAs("plots/png/ppTrack_xt_FullSpectrum.png");
  c1->SaveAs("plots/pdf/ppTrack_xt_FullSpectrum.pdf");
  outF->Close();


  //RpPb bin shift correction
  
  //interpolation points
  float ppintyvals[37+1] = {4.1202/0.5946, 2.8116/0.6211, 1.9778/0.6552, 1.4330/0.6984, 1.0405/0.7219, 0.7719/0.7515, 0.5816/0.7809, 0.3893/0.825, 0.23381/0.866, 0.14395/0.901, 0.09103/0.925, 0.05937/0.965, 0.03906/0.984, 0.014787/1.023, 0.003806/1.052, 0.001181/1.056, 0.0004290/1.048, 0.0001787/1.054, 0.00008152/1.031, 0.00002216/1.023, 0.000004653/1.036, 0.000001402/1.054, 0.0000003180/1.072, 0.00000006850/1.142, 0.00000001971/1.189, 0.000000006013/1.259, 0.000000001910/1.308, 0.0000000007181/1.342, 0.0000000002083/1.382,0.00000000005311/1.407,0.00000000001639/1.363,0.000000000005354/1.381,0.000000000001709/1.316,0,0 };//first 2 points were 4.1202/0.5946*0.999, 2.8116/0.6211*1.003
  //float ppintyvals_pPb[s.ntrkBins+1] = {4.1202*0.999, 2.8116*1.003, 1.9778*1.002, 1.4330*0.994 , 1.0405*1.002, 0.7719*1.002, 0.5816*1.002, 0.3893*1.001, 0.23381*1.000, 0.14395*1.004, 0.09103*1.005, 0.05937*0.997, 0.03906*0.998, 0.014787*1.062, 0.003806*1.042, 0.001181*1.033, 0.0004290*1.024, 0.0001787*1.018, 0.00008152*1.015, 0.00002216* 1.109, 0.000004653*1.061, 0.000001402*1.040, 0.0000003180*1.111, 0.00000006850*1.065, 0.00000001971*1.044, 0.000000006013*1.051, 0.000000001910*1.033, 0.0000000007181*1.024, 0.0000000002083*1.078,0.00000000005311*1.05,0.00000000001639*1.034386,0.000000000005354*1.047316,0.000000000001709*1.032760,0,0 };
  float ppintyvals_pPb[37+1] = { 4.1202, 2.8116,1.9778, 1.4330 , 1.0405, 0.7719, 0.5816, 0.3893, 0.23381, 0.14395, 0.09103, 0.05937, 0.03906, 0.014787, 0.003806, 0.001181, 0.0004290, 0.0001787, 0.00008152, 0.00002216, 0.000004653, 0.000001402, 0.0000003180, 0.00000006850, 0.00000001971, 0.000000006013, 0.000000001910, 0.0000000007181, 0.0000000002083,0.00000000005311,0.00000000001639,0.000000000005354,0.000000000001709,0,0 };//first 2 points were 4.1202, 2.8116
  float ppintyvals_RpPb[37+1] = { 0.5946, 0.6211, 0.6552, 0.6984,0.7219,0.7515, 0.7809, 0.825, 0.866, 0.901, 0.925, 0.965, 0.984, 1.023, 1.052,1.056, 1.048, 1.054, 1.031, 1.023, 1.036, 1.054, 1.072, 1.142, 1.189, 1.259, 1.308, 1.342,1.382,1.407,1.363,1.381,1.316,0,0 };//first 2 points were 0.5946, 0.6211
  float ppintyerrRpPb[37+1] = {0.0027/0.5946,0.0028/0.6211, 0.0030/0.6552,0.0032/0.6984,0.0033/0.7219,0.0034/0.7515,0.0036/0.7809,0.003/0.825,0.003/0.866,0.003/0.901,0.003/0.925,0.003/0.965,0.003/0.984,0.002/1.023,0.002/1.052,0.002/1.056,0.002/1.048,0.003/1.054,0.003/1.031,0.003/1.023,0.005/1.036,0.008/1.054,0.001/1.072,0.002/1.142,0.001/1.189,0.002/1.259,0.002/1.308,0.004/1.342,0.004/1.382,0.008/1.407,0.013/1.363,0.019/1.381,0.030/1.316,0,0 };//first 2 points were 0.0027/0.5946,0.0028/0.6211
  float ppintyerr_pPb[37+1] = {0.0188/4.1202,0.0128/2.8116, 0.0090/1.9778,0.0065/1.4330,0.0047/1.0405,0.0035/0.7719,0.0026/0.5816,0.0013/0.3893,0.00076/0.23381,0.00047/0.14395,0.00029/0.09103,0.00019/0.05937,0.00013/0.03906,0.000026/0.014787,0.000007/0.003806,0.000002/0.001181,0.0000009/0.0004290,0.0000004/0.0001787,0.00000025/0.00008152,0.00000025/0.00002216,0.00000006/0.000004653,0.000000023/0.000001402,0.000000011/0.0000003180,0.0000000004/0.00000006850,0.00000000011/0.00000001971,0.00000000002/0.000000006013,0.000000000009/0.000000001910,0.000000000004/0.0000000007181,0.0000000000020/0.0000000002083,0.0000000000007/0.00000000005311,0.00000000000031/0.00000000001639,0.00000000000016/0.000000000005354,0.000000000000073/0.000000000001709,0,0 };//first 2 points were 0.0188/4.1202,0.0128/2.8116
  TH1D * RpPb = new TH1D("pp_interpolation",";pt;E#frac{d^{3}N}{dPtdyd#phi}",s.ntrkBins,s.xtrkbins);
  TH1D * RpPbSpectrum = new TH1D("RpPb_interp",";pt;RpPb",s.ntrkBins,s.xtrkbins);
  TH1D * pPbSpectrum = new TH1D("pPb_data",";pt;E#frac{d^{3}N}{dPtdyd#phi}",s.ntrkBins,s.xtrkbins);
  for(int i=1; i<s.ntrkBins+1; i++) RpPb->SetBinContent(i,ppintyvals[i-1]);
  for(int i=1; i<s.ntrkBins+1; i++) RpPb->SetBinError(i,TMath::Power(TMath::Abs(TMath::Power(ppintyerrRpPb[i-1],2)-TMath::Power(ppintyerr_pPb[i-1],2)),0.5)*ppintyvals[i-1]);
  for(int i=1; i<s.ntrkBins+1; i++) RpPbSpectrum->SetBinContent(i,ppintyvals_RpPb[i-1]);
  for(int i=1; i<s.ntrkBins+1; i++) RpPbSpectrum->SetBinError(i,ppintyvals_RpPb[i-1]*ppintyerrRpPb[i-1]);
  for(int i=1; i<s.ntrkBins+1; i++) pPbSpectrum->SetBinContent(i,ppintyvals_pPb[i-1]);
  for(int i=1; i<s.ntrkBins+1; i++) pPbSpectrum->SetBinError(i,ppintyvals_pPb[i-1]*ppintyerr_pPb[i-1]);
  
  //float ppScale = RpPb->GetBinContent(17)/s.pp->GetBinContent(17);
  //s.pp->Scale(ppScale); 

  TFile * binShift = TFile::Open("BinningAndResolutionCorrection_TrackTrigger.root","read");
  TH1D *hBinningAndResol_extendedPt = (TH1D*)binShift->Get("hPt_pseudo2_copy1");
  hBinningAndResol_extendedPt->SetDirectory(0);
  binShift->Close();