make_pi0_xcheck.C

Go to the documentation of this file.

#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 "CAFAna/Core/Cut.h"
#include "CAFAna/Cuts/NueCutsFirstAna.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 "CAFAna/Vars/PPFXWeights.h"

#include "CAFAna/nue/pi0xcheck/make_pi0_xcheck.h"

#include "StandardRecord/Proxy/SRProxy.h"

using namespace ana;

void make_pi0_xcheck(bool isFHC = 1)
{
  const std::string tag = isFHC ? "FHC" : "RHC";

  std::string fData = isFHC ? "prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_numi_fhc_full_v1_addShortSimpleCVN_goodruns" :
    "prod_caf_R17-11-14-prod4reco.neutron-respin.b_nd_numi_rhc_full_v1_goodruns" ;
  std::string fMC   = isFHC ? "prod_caf_R17-11-14-prod4reco.neutron-respin.c_nd_genie_nonswap_fhc_nova_v08_full_v1" :
    "prod_caf_R17-11-14-prod4reco.neutron-respin.c_nd_genie_nonswap_rhc_nova_v08_full_v1" ;

  SpectrumLoaderBase* loader = new SpectrumLoader(fData);
  SpectrumLoaderBase* loaderMC = new SpectrumLoader(fMC);
  
  std::cout << "MC:   " << fMC << std::endl;
  std::cout << "Data: " << fData << std::endl;
  
  loader->SetSpillCut(kStandardSpillCuts);
  loaderMC->SetSpillCut(kStandardDQCuts);

  const Var weight = kXSecCVWgt2018*kPPFXFluxCVWgt;
  
  Spectrum* dataSpects[kNSels][kNPlots];
  Spectrum* mcSpects[kNSels][kNPlots];
  Spectrum* bkgSpects[kNSels][kNPlots];
  
  Spectrum* dataSpects2[kNSels][kN2D];
  Spectrum* mcSpects2[kNSels][kN2D];
  Spectrum* bkgSpects2[kNSels][kN2D];

  Spectrum* dataSpectsMulti[kNSels][kNMulti];
  Spectrum* mcSpectsMulti[kNSels][kNMulti];
  Spectrum* bkgSpectsMulti[kNSels][kNMulti];

  for(unsigned int i = 0; i < kNSels; ++i){
    for(unsigned int j = 0; j < kNPlots; ++j){
      dataSpects[i][j] = new Spectrum(plots[j].label, plots[j].bins, *loader, plots[j].var, sels[i].cut);
      mcSpects[i][j]   = new Spectrum(plots[j].label, plots[j].bins, *loaderMC, plots[j].var, sels[i].cut, 
                                      kNoShift, weight);
      bkgSpects[i][j]  = new Spectrum(plots[j].label, plots[j].bins, *loaderMC, plots[j].var, sels[i].cut && !kTruePi0, 
                                      kNoShift, weight);
    }
  }

  for(unsigned int i = 0; i < kNSels; ++i){
    for(unsigned int j = 0; j < kNMulti; ++j){
      dataSpectsMulti[i][j] = new Spectrum(plotsMulti[j].label, plotsMulti[j].bins, *loader, plotsMulti[j].var, sels[i].cut);
      mcSpectsMulti[i][j]   = new Spectrum(plotsMulti[j].label, plotsMulti[j].bins, *loaderMC, plotsMulti[j].var, sels[i].cut, 
                                           kNoShift, weight);
      bkgSpectsMulti[i][j]  = new Spectrum(plotsMulti[j].label, plotsMulti[j].bins, *loaderMC, plotsMulti[j].var, sels[i].cut && !kTruePi0, 
                                           kNoShift, weight);
    }
  }

  for(unsigned int i = 0; i < kNSels; ++i){
    for(unsigned int j = 0; j < kN2D; ++j){
      dataSpects2[i][j] = new Spectrum(plots2d[j].label1, plots2d[j].label2, *loader, plots2d[j].bins1, plots2d[j].var1, 
                                       plots2d[j].bins2, plots2d[j].var2, sels[i].cut);
      mcSpects2[i][j] = new Spectrum(plots2d[j].label1, plots2d[j].label2, *loaderMC, plots2d[j].bins1, plots2d[j].var1, 
                                     plots2d[j].bins2, plots2d[j].var2, sels[i].cut, 
                                     kNoShift, weight);
      bkgSpects2[i][j] = new Spectrum(plots2d[j].label1, plots2d[j].label2, *loaderMC, plots2d[j].bins1, plots2d[j].var1, 
                                      plots2d[j].bins2, plots2d[j].var2, sels[i].cut && !kTruePi0, 
                                      kNoShift, weight);
    }
  }

  loader->Go();
  loaderMC->Go();

  TFile fout(("make_pi0_xcheck_output_"+tag+".root").c_str(), "RECREATE");

  for(unsigned int i=0;i<kNSels;++i){
    for(unsigned int j=0;j<kNPlots;++j){
      std::string mysuffix=plots[j].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]->SaveTo(&fout, d.c_str());
      mcSpects[i][j]->SaveTo(&fout, m.c_str());
      bkgSpects[i][j]->SaveTo(&fout, b.c_str());
    }
  }

  for(unsigned int i=0;i<kNSels;++i){
    for(unsigned int j=0;j<kN2D;++j){
      std::string mysuffix=plots2d[j].suffix + "_" + sels[i].suffix;
      
      std::string d = tag+"_data_"+mysuffix;
      std::string m = tag+"_mc_"+mysuffix;
      std::string b = tag+"_bg_"+mysuffix;
      
      dataSpects2[i][j]->SaveTo(&fout, d.c_str());
      mcSpects2[i][j]->SaveTo(&fout, m.c_str());
      bkgSpects2[i][j]->SaveTo(&fout, b.c_str());
    }
  }

  for(unsigned int i=0;i<kNSels;++i){
    for(unsigned int j=0;j<kNMulti;++j){
      std::string mysuffix=plotsMulti[j].suffix + "_" + sels[i].suffix;
      
      std::string d = tag+"_data_"+mysuffix;
      std::string m = tag+"_mc_"+mysuffix;
      std::string b = tag+"_bg_"+mysuffix;
      
      dataSpectsMulti[i][j]->SaveTo(&fout, d.c_str());
      mcSpectsMulti[i][j]->SaveTo(&fout, m.c_str());
      bkgSpectsMulti[i][j]->SaveTo(&fout, b.c_str());
    }
  }
}