plot_pi0_xcheck.C File Reference

#include <sstream>
#include <iomanip>
#include <fstream>
#include "TCanvas.h"
#include "TH1.h"
#include "TF1.h"
#include "TFile.h"
#include "TLegend.h"
#include "TLegendEntry.h"
#include "TArrayD.h"
#include "TAxis.h"
#include "TLatex.h"
#include "TStyle.h"
#include "TLine.h"
#include "TProfile.h"
#include "CAFAna/Core/Cut.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Cuts/NueCutsFirstAna.h"
#include "3FlavorAna/Plotting/NuePlotStyle.h"
#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Analysis/Style.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Vars/Vars.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "StandardRecord/Proxy/SRProxy.h"
#include "/nova/app/users/grohmc/pi_S18-02-09/CAFAna/nue/pi0xcheck/make_pi0_xcheck.h"

Go to the source code of this file.

Functions
void	Preliminary ()
void	legend (TH1 *data, TH1 *mc, TH1 *bkgd, double x1=0.58, double y1=0.64, double x2=0.86, double y2=0.84, bool isData=1)
void	HCtag (TString tag, double pot, double size=0.8)
void	tex (std::string str)
void	FitValues (TF1 *data, TF1 *mc, TH1 *hdata, TH1 *hmc, double x1=0.50, double y1=0.30, double x2=0.85, double y2=0.55)
std::string	ToString (double value)
TCanvas *	RatioPlot (std::vector< TH1 * > topHistos, std::vector< TString > topOption, std::vector< TH1 * > bottomHistos, std::vector< TString > bottomOption)
void	plot_pi0_xcheck (const std::string inFile, const bool isFHC)

Function Documentation

void FitValues	(	TF1 *	data,
		TF1 *	mc,
		TH1 *	hdata,
		TH1 *	hmc,
		double	x1 = 0.50,
		double	y1 = 0.30,
		double	x2 = 0.85,
		double	y2 = 0.55
	)

Definition at line 93 of file plot_pi0_xcheck.C.

References Munits::mm, Munits::ms, x1, submit_syst::x2, y1, and submit_syst::y2.

Referenced by plot_pi0_xcheck().

                                                                                                                     {

  TLegend *leg_dm = new TLegend(x1,y1,x2,y2);
    
  char dmean[50], dsigma[50], mmean[50], msigma[50];
  sprintf(dmean, "Data #mu: #bf{%2.1f} #pm %2.1f MeV", data->GetParameter(1), data->GetParError(1));
  sprintf(dsigma,"Data #sigma:   %2.1f #pm %2.1f MeV", data->GetParameter(2), data->GetParError(2));
  sprintf(mmean, "MC #mu:   #bf{%2.1f} #pm %2.1f MeV", mc->GetParameter(1), mc->GetParError(1));
  sprintf(msigma,"MC #sigma:     %2.1f #pm %2.1f MeV", mc->GetParameter(2), mc->GetParError(2));
  
  TLegendEntry* dm = leg_dm->AddEntry((TObject*)0, dmean, "");
  dm->SetTextColor(hdata->GetLineColor());
  TLegendEntry* ds = leg_dm->AddEntry((TObject*)0, dsigma, "");
  ds->SetTextColor(hdata->GetLineColor());
  TLegendEntry* mm = leg_dm->AddEntry((TObject*)0, mmean, "");
  mm->SetTextColor(hmc->GetLineColor());
  TLegendEntry* ms = leg_dm->AddEntry((TObject*)0, msigma, "");
  ms->SetTextColor(hmc->GetLineColor());

  leg_dm->SetBorderSize(0); //no border for legend
  leg_dm->SetFillColor(0);  //fill colour is white
  //leg_dm->SetTextSize(2/40.);
  leg_dm->Draw();
  
}

void HCtag	(	TString	tag,
		double	pot,
		double	size = 0.8
	)

Definition at line 64 of file plot_pi0_xcheck.C.

References confusionMatrixTree::out, pot, gen_hdf5record::size, string, tex(), and art::to_string().

Referenced by plot_pi0_xcheck().

{
  std::stringstream out;
  out << std::setprecision(3);
  out << pot;
  std::string potstr = out.str();
  std::string coef = out.str().substr(0,4);
  std::string power = out.str().substr(6,2);
  potstr = coef+"x10^{"+power+"}";
  //potstr = "3.10x10^{20}";
  
  TLatex* tex = new TLatex(.125, .85, "#scale["+std::to_string(size)+"]{"+potstr+" POT-equiv "+tag+"}");
  tex->SetNDC();
  tex->SetTextSize(0.043);
  tex->SetTextFont(42);
  tex->SetTextAlign(12);
  tex->Draw();
}

void legend	(	TH1 *	data,
		TH1 *	mc,
		TH1 *	bkgd,
		double	x1 = 0.58,
		double	y1 = 0.64,
		double	x2 = 0.86,
		double	y2 = 0.84,
		bool	isData = 1
	)

Definition at line 48 of file plot_pi0_xcheck.C.

References MECModelEnuComparisons::leg, x1, submit_syst::x2, y1, and submit_syst::y2.

Referenced by plot_pi0_xcheck().

                                                                                                                    {
  TLegend *leg = new TLegend(x1,y1,x2,y2);

  if(isData){
    TLegendEntry* dl = leg->AddEntry(data,"Near Det Data","lpe");
    dl->SetTextColor(data->GetLineColor());
  }
  TLegendEntry* ml = leg->AddEntry(mc,"MC #pi^{0} Signal","l");
  ml->SetTextColor(mc->GetLineColor());
  TLegendEntry* bl = leg->AddEntry(bkgd,"MC Background","f");
  bl->SetTextColor(bkgd->GetLineColor());
  leg->SetBorderSize(0); //no border for legend
  leg->SetFillColor(0);  //fill colour is white
  leg->Draw();
}

void plot_pi0_xcheck	(	const std::string	inFile,
		const bool	isFHC
	)

Definition at line 208 of file plot_pi0_xcheck.C.

References b, plot_validation_datamc::c, demo5::c1, demo5::c2, chisquared::c3, ana::CenterTitles(), om::cout, d, MakeMiniprodValidationCuts::f, ana::FillWithDimColor(), FitValues(), submit_syst::fout, make_syst_table_plots::h, HCtag(), MECModelEnuComparisons::i, inFile, calib::j, kN2D, kNMulti, kNPlots, kNSels, ana::kTotalBackgroundColor, ana::kTotalMCColor, legend(), m, extractScale::mean, plots, plots2d, plotsMulti, POT, Preliminary(), ana::RatioPlot(), submit_cafana::rc, runNovaSAM::release, ana::sels, string, getGoodRuns4SAM::tag, tex(), ToString(), ana::Spectrum::ToTH1(), and ana::Spectrum::ToTH2().

{
  const std::string tag = isFHC ? "FHC" : "RHC";

  Spectrum* dataSpects[kNSels][kNPlots+kNMulti];
  Spectrum* mcSpects[kNSels][kNPlots+kNMulti];
  Spectrum* bkgSpects[kNSels][kNPlots+kNMulti];

  TFile fout(("plot_pi0_xcheck_output_"+tag+".root").c_str(), "RECREATE");
  
  double integrals[kNSels][kNPlots+kNMulti][3];
  double fwhmlo[kNSels];
  double fwhmhi[kNSels];

  for(unsigned int j=0;j<kNPlots+kNMulti;++j){
    for(unsigned int i=0;i<kNSels;++i){
      std::string mysuffix;
      if( j < kNPlots) mysuffix = plots[j].suffix + "_" + sels[i].suffix;
      else mysuffix = plotsMulti[j-kNPlots].suffix + "_" + sels[i].suffix;
      
      std::string d = tag+"_data_"+mysuffix;
      std::string m = tag+"_mc_"+mysuffix;
      std::string b = tag+"_bg_"+mysuffix;
      
      dataSpects[i][j] = LoadFromFile<Spectrum>(inFile, d).release();
      mcSpects[i][j] = LoadFromFile<Spectrum>(inFile, m).release();
      bkgSpects[i][j] = LoadFromFile<Spectrum>(inFile, b).release();
      
      TCanvas *c = new TCanvas(mysuffix.c_str(),mysuffix.c_str());
      TH1* hMC = mcSpects[i][j]->ToTH1(dataSpects[i][j]->POT());
      hMC->SetLineColor(kTotalMCColor);
      
      TH1* hData = dataSpects[i][j]->ToTH1(dataSpects[i][j]->POT());
      hData->SetMarkerStyle(kFullCircle);
      
      double mean, bin1, bin2, fwhm;
      mean  = hData->GetBinCenter(hData->GetMaximumBin());
      bin1  = hData->FindFirstBinAbove(hData->GetMaximum()/2);
      bin2  = hData->FindLastBinAbove(hData->GetMaximum()/2);
      fwhm  = hData->GetBinCenter(bin2) - hData->GetBinCenter(bin1);
      if(j==1){
        fwhmlo[i]=mean-fwhm/2;
        fwhmhi[i]=mean+fwhm/2;
      }

      TF1 *fData = new TF1("fData","gaus", mean-fwhm/2 , mean+fwhm/2);
      fData->SetLineColor(kBlack);
      if(j==1) hData->Fit("fData","R0","");
      
      mean  = hMC->GetBinCenter(hMC->GetMaximumBin());
      bin1  = hMC->FindFirstBinAbove(hMC->GetMaximum()/2);
      bin2  = hMC->FindLastBinAbove(hMC->GetMaximum()/2);
      fwhm  = hMC->GetBinCenter(bin2) - hMC->GetBinCenter(bin1);
      
      TF1 *fMC = new TF1("fMC","gaus", mean-fwhm/2 , mean+fwhm/2);
      fMC->SetLineColor(kTotalMCColor);
      if(j==1) hMC->Fit("fMC","R0+","");

      TH1* h = bkgSpects[i][j]->ToTH1(dataSpects[i][j]->POT());
      h->SetLineColor(kTotalBackgroundColor);
      FillWithDimColor(h);
      
      hMC->GetYaxis()->SetRangeUser(0,1.2*hMC->GetMaximum());
      if(j==1) hMC->GetYaxis()->SetTitle(("Events / "+ToString(hMC->GetBinWidth(1))+" MeV").c_str());
      else hMC->GetYaxis()->SetTitle("Events");

      hMC->Draw("hist");
      hData->Draw("ep same");
      h->Draw("hist same");
      
      CenterTitles(hMC);
      CenterTitles(hData);
      CenterTitles(h);
      
      if(j==1) FitValues(fData, fMC, hData, hMC);
      legend(hData,hMC,h);
      Preliminary();
      HCtag(tag,dataSpects[i][j]->POT());
      
      double mData = fData->GetParameter(1);
      double meData = fData->GetParError(1);
      double sData = fData->GetParameter(2);
      double seData = fData->GetParError(2);
      
      std::cout << "Data Mean:  " << mData << " +- " << meData << "\n";
      std::cout << "Data Sigma: " << sData << " +- " << seData << "\n";
      
      integrals[i][j][0] = hData->Integral();
      integrals[i][j][1] = hMC->Integral();
      integrals[i][j][2] = h->Integral();
      
      double mMC = fMC->GetParameter(1);
      double meMC = fMC->GetParError(1);
      double sMC = fMC->GetParameter(2);
      double seMC = fMC->GetParError(2);
      
      std::cout << "MC Mean:  " << mMC << " +- " << meMC << "\n";
      std::cout << "MC Sigma: " << sMC << " +- " << seMC << "\n";
      
      gPad->RedrawAxis();
      
      fout.cd();
      c->Write();
      c->Print(("plots/data_mc_"+mysuffix+"_"+tag+".root").c_str());
      c->Print(("plots/data_mc_"+mysuffix+"_"+tag+".pdf").c_str());

      //if(j == 1){
        TH1 *hRat = (TH1*)hData->Clone("hRat");
        hRat->Divide(hMC);
        TCanvas *rc = RatioPlot({hMC,hData,h},{"hist","ep","hist"},{hRat},{"ep"});
        
        if(j == 1) FitValues(fData, fMC, hData, hMC,0.48,0.40,0.85,0.55);
        legend(hData,hMC,h,0.56,0.64,0.86,0.79);
        HCtag(tag,dataSpects[i][j]->POT(),0.5);
        
        fout.cd();
        rc->Write();
        rc->Print(("plots/data_mc_rat_"+mysuffix+"_"+tag+".root").c_str());
        rc->Print(("plots/data_mc_rat_"+mysuffix+"_"+tag+".pdf").c_str());
        //}
    }
  }

  gStyle->SetPalette(51);

  Spectrum* dataSpect2d[kNSels][kN2D];
  Spectrum* mcSpect2d[kNSels][kN2D];
  Spectrum* bkgSpect2d[kNSels][kN2D];

  for(int j = 0;j < kN2D; ++j){
    for(int i = 0;i < kNSels; ++i){
      std::string mysuffix =  plots2d[j].suffix + "_" + sels[i].suffix;
      dataSpect2d[i][j] = LoadFromFile<Spectrum>(inFile, (tag+"_data_"+plots2d[j].suffix+"_"+sels[i].suffix).c_str()).release();
      mcSpect2d[i][j] = LoadFromFile<Spectrum>(inFile, (tag+"_mc_"+plots2d[j].suffix+"_"+sels[i].suffix).c_str()).release();
      bkgSpect2d[i][j] = LoadFromFile<Spectrum>(inFile, (tag+"_bg_"+plots2d[j].suffix+"_"+sels[i].suffix).c_str()).release();

      TH2* hData = dataSpect2d[i][j]->ToTH2(dataSpect2d[i][j]->POT());
      TH2* hMC = mcSpect2d[i][j]->ToTH2(dataSpect2d[i][j]->POT());
      TH2* h = bkgSpect2d[i][j]->ToTH2(dataSpect2d[i][j]->POT());

      if(j>2){
        TCanvas *c = new TCanvas(mysuffix.c_str(),mysuffix.c_str());
        TH1D *pData = new TH1D("pData","pData",plots2d[j].bins1.NBins(),plots2d[j].bins1.Min(),plots2d[j].bins1.Max());
        TH1D *pMC   = new TH1D("pMC","pMC",plots2d[j].bins1.NBins(),plots2d[j].bins1.Min(),plots2d[j].bins1.Max());

        for(int k=1;k<hData->GetNbinsX();++k){
          TH1D *hProj = hData->ProjectionY("test",k,k);

          double mean, bin1, bin2, fwhm;
          mean  = hProj->GetBinCenter(hProj->GetMaximumBin());
          bin1  = hProj->FindFirstBinAbove(hProj->GetMaximum()/2);
          bin2  = hProj->FindLastBinAbove(hProj->GetMaximum()/2);
          fwhm  = hProj->GetBinCenter(bin2) - hProj->GetBinCenter(bin1);          

          TF1 *f = new TF1("f","gaus", mean-fwhm/2 , mean+fwhm/2);
          if(hProj->GetEntries() != 0){
            hProj->Fit("f","RQ0","");
            pData->SetBinContent(k,f->GetParameter(1));
            pData->SetBinError(k,f->GetParameter(2));
          }
        }
        for(int k=1;k<hData->GetNbinsX();++k){
          TH1D *hProj = hMC->ProjectionY("test",k,k);

          double mean, bin1, bin2, fwhm;
          mean  = hProj->GetBinCenter(hProj->GetMaximumBin());
          bin1  = hProj->FindFirstBinAbove(hProj->GetMaximum()/2);
          bin2  = hProj->FindLastBinAbove(hProj->GetMaximum()/2);
          fwhm  = hProj->GetBinCenter(bin2) - hProj->GetBinCenter(bin1);          

          TF1 *f = new TF1("f","gaus", mean-fwhm/2 , mean+fwhm/2);
          if(hProj->GetEntries() != 0){
            hProj->Fit("f","RQ0","");
            pMC->SetBinContent(k,f->GetParameter(1));
            pMC->SetBinError(k,f->GetParameter(2));
          }
        }
        
        pData->SetMarkerStyle(kFullCircle);     
        pMC->SetLineColor(kTotalMCColor);

        pMC->SetTitle("");
        pMC->GetYaxis()->SetRangeUser(0,500);
        pMC->GetXaxis()->SetTitle(plots2d[j].label1.c_str());
        pMC->GetYaxis()->SetTitle(plots2d[j].label2.c_str());
        
        pMC->Draw("ep");
        pData->Draw("ep same");

        c->Update();

        TLine *lone = new TLine();
        lone->SetLineStyle(3);
        lone->SetLineColor(kGray+3);
        lone->DrawLine(c->GetUxmin(),fwhmlo[i],c->GetUxmax(),fwhmlo[i]);
        lone->DrawLine(c->GetUxmin(),fwhmhi[i],c->GetUxmax(),fwhmhi[i]);

        fout.cd();
        c->Write();
        c->Print(("plots/data_mc_"+mysuffix+"_"+tag+"_prof.root").c_str());
        c->Print(("plots/data_mc_"+mysuffix+"_"+tag+"_prof.pdf").c_str());
      }
      TCanvas *c1 = new TCanvas((mysuffix+"1").c_str(),(mysuffix+"1").c_str());
      hData->Draw("COLZ");
      tex("Data");
      gPad->SetLogz();
      fout.cd();
      c1->Write();
      c1->Print(("plots/data_mc_"+mysuffix+"_"+tag+"_data.root").c_str());
      c1->Print(("plots/data_mc_"+mysuffix+"_"+tag+"_data.pdf").c_str());
      
      TCanvas *c2 = new TCanvas((mysuffix+"2").c_str(),(mysuffix+"2").c_str());
      hMC->Draw("COLZ");
      tex("MC");
      gPad->SetLogz();
      fout.cd();
      c2->Write();
      c2->Print(("plots/data_mc_"+mysuffix+"_"+tag+"_mc.root").c_str());
      c2->Print(("plots/data_mc_"+mysuffix+"_"+tag+"_mc.pdf").c_str());
      
      TCanvas *c3 = new TCanvas((mysuffix+"3").c_str(),(mysuffix+"3").c_str());
      h->Draw("COLZ");
      tex("Background");
      gPad->SetLogz();
      fout.cd();
      c3->Write();
      c3->Print(("plots/data_mc_"+mysuffix+"_"+tag+"_bkg.root").c_str());
      c3->Print(("plots/data_mc_"+mysuffix+"_"+tag+"_bkg.pdf").c_str());
      
    }
  }

  ofstream txtout;
  txtout.open("integrals.csv");
  
  txtout<<"Selection,Plot,Data,MC,Bkg,Signal,\n";
  for(unsigned int s=0;s<kNPlots+kNMulti;++s){
    for(unsigned int k=0;k<kNSels;++k){
      if( s < kNPlots ) txtout<<sels[k].suffix+","+plots[s].suffix+",";
      else  txtout<<sels[k].suffix+","+plotsMulti[s-kNPlots].suffix+",";
      for(unsigned int i=0;i<3;++i){
        txtout<<integrals[k][s][i]<<",";
        if(i == 2){
          txtout<<integrals[k][s][1]-integrals[k][s][2]<<",";
        }
      }
      txtout<<"\n";
    }
  }
  txtout.close();

}

void Preliminary

(

)

Definition at line 38 of file plot_pi0_xcheck.C.

References prelim.

Referenced by plot_pi0_xcheck().

{
  TLatex* prelim = new TLatex(.9, .95, "NOvA Preliminary");
  prelim->SetTextColor(kAzure);
  prelim->SetNDC();
  prelim->SetTextSize(2/30.);
  prelim->SetTextAlign(32);
  prelim->Draw();
}

TCanvas* RatioPlot	(	std::vector< TH1 * >	topHistos,
		std::vector< TString >	topOption,
		std::vector< TH1 * >	bottomHistos,
		std::vector< TString >	bottomOption
	)

Definition at line 127 of file plot_pi0_xcheck.C.

References plot_validation_datamc::c, plot_validation_datamc::Clone(), make_syst_table_plots::h, h1, h3, plot_validation_datamc::p1, and plot_validation_datamc::p2.

{
  gStyle->SetTitleStyle(0);
  TCanvas *c = new TCanvas("c", "canvas", 600, 700);
  //define top and bottom pads
  TPad* p1 = new TPad("", "", 0, 0, 1, 1);
  TPad* p2 = new TPad("", "", 0, 0, 1, 1);
  p1->SetBottomMargin(.3);
  p2->SetTopMargin(.7);

  for(auto p:{p1,p2}){
    p->SetFillStyle(0);
    p->Draw();
  }

  // format histograms
  auto h1 = (TH1*) topHistos[0]->Clone();
  auto h3 = (TH1*) bottomHistos[0]->Clone();

  for(auto & h:{h1,h3}){
    h->SetStats(0);

    h->GetYaxis()->SetTitleSize(26);
    h->GetYaxis()->SetTitleFont(43);
    h->GetYaxis()->SetTitleOffset(1.3);
    h->GetYaxis()->SetLabelFont(43);

    h->GetXaxis()->SetTitleOffset(1.2);
    h->GetXaxis()->SetTitleSize(28);
    h->GetXaxis()->SetTitleFont(43);
    h->GetXaxis()->SetLabelFont(43);
  }

  h1->GetYaxis()->SetLabelSize(18);
  
  h1->GetXaxis()->SetLabelSize(0);
  h1->GetXaxis()->SetTitleSize(0);

  h3->SetTitle("");
  h3->GetYaxis()->SetTitle("Data/MC");
  h3->GetYaxis()->SetRangeUser(0.5,1.5);
  h3->GetYaxis()->SetLabelSize(15);
  h3->GetYaxis()->SetDecimals();

  h3->GetXaxis()->SetTitle(h1->GetXaxis()->GetTitle());
  h3->GetXaxis()->SetLabelSize(18);
  h3->GetXaxis()->SetLabelOffset(0.005);

  p1->cd();

  h1->Draw(topOption[0]);

  for(unsigned int ii = 1; ii < topHistos.size(); ++ii) {
    topHistos[ii]->Draw("same "+topOption[ii]);
  }
  h1->Draw("same "+topOption[0]);

  p1->RedrawAxis(); //white histograms cover axes

  p2->cd();
  h3->Draw(bottomOption[0]);
  for(unsigned int ii = 1; ii < bottomHistos.size(); ++ii){
    bottomHistos[ii]->Draw("same " + bottomOption[ii]);
  }
  h3->Draw("same " + bottomOption[0]);

  TLine *lone = new TLine();
  lone->SetLineStyle(3);
  lone->SetLineColor(kGray+3);
  p2->Update();
  lone->DrawLine(p2->GetUxmin(),1.0,p2->GetUxmax(),1.0);
  lone->DrawLine(p2->GetUxmin(),1.2,p2->GetUxmax(),1.2);
  lone->DrawLine(p2->GetUxmin(),0.8,p2->GetUxmax(),0.8);
  p2->RedrawAxis();

  return c;
}

void tex

(

std::string

str

)

Definition at line 83 of file plot_pi0_xcheck.C.

Referenced by HCtag(), and plot_pi0_xcheck().

                       {
  TLatex* tex = new TLatex(.1, .925, str.c_str());

  tex->SetNDC();
  tex->SetTextSize(0.043);
  tex->SetTextFont(42);
  tex->SetTextAlign(12);
  tex->Draw();
}

std::string ToString

(

double

value

)

Definition at line 119 of file plot_pi0_xcheck.C.

References confusionMatrixTree::out.

Referenced by cmf::MetaDataLite::Period(), cmf::MetaData::Period(), and plot_pi0_xcheck().

                                 {
  std::stringstream out;
  out << std::setprecision(3);
  out << value;
  return out.str();
}