mrbrem_get_reweighted_spectra.C

Go to the documentation of this file.

#include "CAFAna/Core/Binning.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/Var.h"
#include "CAFAna/Core/MultiVar.h"

//#include "CAFAna/Cuts/BPFCuts.h"
//#include "CAFAna/Vars/BPFVars.h"

#include "CAFAna/Cuts/Cuts.h"
#include "3FlavorAna/Cuts/NueCuts2020.h"
#include "CAFAna/Cuts/TimingCuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "CAFAna/Cuts/TruthCuts.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/PPFXWeights.h"
#include "CAFAna/Vars/XsecTunes.h"

#include "CAFAna/Vars/Vars.h"
#include "3FlavorAna/Vars/NueVars.h"
#include "3FlavorAna/Vars/MRVars.h"
#include "3FlavorAna/Vars/NueVarsExtra.h"

#include "StandardRecord/Proxy/SRProxy.h"

#include <cassert>

#include "TFile.h"
#include "TCanvas.h"
#include "TH2.h"
#include "TH1.h"
#include "TStyle.h"
#include "TLegend.h"
#include "TLatex.h"
#include "TArrayD.h"
#include "TTree.h"
#include "TLeaf.h"
#include "stdio.h"
#include <sstream>
#include <string>
#include <iostream>
#include <fstream>
#include <string>


using namespace ana;
using namespace std;


void mrbrem_get_reweighted_spectra(std::string mode = "fhc")
{

  std::string filename = "${CONDOR_DIR_INPUT}/MRbrem_extract_weights_"+mode+".root"; // ${CONDOR_DIR_INPUT}
  struct HistDef
  {
     std::string name;
     HistAxis axis;
  };

  const int kNumVars = 1;
  const int kNumSels = 8;

  const HistDef defs[kNumVars] = {
        {"shwangle", {"Cos(#theta_{Z})", Binning::Simple(60,0,1.0), kShwAngle}},
  };

  struct CutDef
  {
     std::string name;
     Cut cut;
  };

  const CutDef sels[kNumSels] = {
      {"core_presel", kNue2020CorePresel_MRBrem},
      {"periph_presel", kNue2020PeripheralPresel_MRBrem},
      {"core_full", kNue2020FD_MRBrem},
      {"periph_full", kNue2020FDPeripheral_MRBrem},
      {"periph_preselbdt", kNue2020PeripheralPresel_MRBrem && kCosPIDPeriBDT > 0.57},
      {"periph_preselcvn", kNue2020PeripheralPresel_MRBrem && kCVNe_looseptp >= 0.97},
      {"periph_preselen", kNue2020PeripheralPresel_MRBrem && kShwE > 1},
      {"periph_fullen", kNue2020FDPeripheral_MRBrem && kShwE > 1},
  };

  Spectrum* spec_mrbrem_data[kNumSels][kNumVars];
  Spectrum* spec_mrbrem_cry[kNumSels][kNumVars];
  Spectrum* spec_genie_fd[kNumSels][kNumVars];

  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
      const char* selName = sels[selIdx].name.c_str();
      for (int varIdx=0; varIdx<kNumVars;varIdx++){
           const char* varName = defs[varIdx].name.c_str();
           spec_mrbrem_data[selIdx][varIdx] = LoadFromFile<Spectrum>(filename, TString::Format("%s/mrbrem_data_%s", selName, varName).Data()).release();
           spec_mrbrem_cry[selIdx][varIdx] = LoadFromFile<Spectrum>(filename, TString::Format("%s/mrbrem_mc_%s", selName, varName).Data()).release();
           spec_genie_fd[selIdx][varIdx] = LoadFromFile<Spectrum>(filename, TString::Format("%s/genie_fd_%s", selName, varName).Data()).release();
      }
  }

  double livetime = spec_mrbrem_data[0][0]->Livetime();
  double pot = spec_genie_fd[0][0]->POT();

  TH1D* hspec_mrbrem_data[kNumSels][kNumVars];
  TH1D* hspec_mrbrem_cry[kNumSels][kNumVars];
  TH1D* hspec_genie_fd[kNumSels][kNumVars];
  
  for (int selIdx = 0; selIdx < kNumSels; ++selIdx){
      for (int varIdx=0; varIdx<kNumVars;varIdx++){
           hspec_mrbrem_data[selIdx][varIdx] = spec_mrbrem_data[selIdx][varIdx]->ToTH1(livetime,kLivetime); 
           hspec_mrbrem_cry[selIdx][varIdx] = spec_mrbrem_cry[selIdx][varIdx]->ToTH1(livetime,kLivetime);
           hspec_genie_fd[selIdx][varIdx] = spec_genie_fd[selIdx][varIdx]->ToTH1(pot);
      }
  }

  TH1D* hgenie_core_presel = (TH1D*)hspec_genie_fd[0][0]->Clone();
  TH1D* hgenie_peri_presel = (TH1D*)hspec_genie_fd[1][0]->Clone();
  TH1D* hgenie_peri_preselbdt = (TH1D*)hspec_genie_fd[4][0]->Clone();
  TH1D* hgenie_peri_preselcvn = (TH1D*)hspec_genie_fd[5][0]->Clone();
  TH1D* hgenie_peri_preselen = (TH1D*)hspec_genie_fd[6][0]->Clone();
  TH1D* hgenie_core_presel_data = (TH1D*)hspec_genie_fd[0][0]->Clone();
  TH1D* hgenie_peri_presel_data = (TH1D*)hspec_genie_fd[1][0]->Clone();
  TH1D* hgenie_peri_preselbdt_data = (TH1D*)hspec_genie_fd[4][0]->Clone();
  TH1D* hgenie_peri_preselcvn_data = (TH1D*)hspec_genie_fd[5][0]->Clone();
  TH1D* hgenie_peri_preselen_data = (TH1D*)hspec_genie_fd[6][0]->Clone();

  hgenie_core_presel->Scale(hspec_mrbrem_cry[0][0]->Integral()/hspec_genie_fd[0][0]->Integral());
  hgenie_core_presel->Divide(hspec_mrbrem_cry[0][0]);

  hgenie_peri_presel->Scale(hspec_mrbrem_cry[1][0]->Integral()/hspec_genie_fd[1][0]->Integral());
  hgenie_peri_presel->Divide(hspec_mrbrem_cry[1][0]);

  hgenie_peri_preselbdt->Scale(hspec_mrbrem_cry[4][0]->Integral()/hspec_genie_fd[4][0]->Integral());
  hgenie_peri_preselbdt->Divide(hspec_mrbrem_cry[4][0]);

  hgenie_peri_preselcvn->Scale(hspec_mrbrem_cry[5][0]->Integral()/hspec_genie_fd[5][0]->Integral());
  hgenie_peri_preselcvn->Divide(hspec_mrbrem_cry[5][0]);

  hgenie_peri_preselen->Scale(hspec_mrbrem_cry[6][0]->Integral()/hspec_genie_fd[6][0]->Integral());
  hgenie_peri_preselen->Divide(hspec_mrbrem_cry[6][0]);

  hgenie_core_presel_data->Scale(hspec_mrbrem_data[0][0]->Integral()/hspec_genie_fd[0][0]->Integral());
  hgenie_core_presel_data->Divide(hspec_mrbrem_data[0][0]);

  hgenie_peri_presel_data->Scale(hspec_mrbrem_data[1][0]->Integral()/hspec_genie_fd[1][0]->Integral());
  hgenie_peri_presel_data->Divide(hspec_mrbrem_data[1][0]);

  hgenie_peri_preselbdt_data->Scale(hspec_mrbrem_data[4][0]->Integral()/hspec_genie_fd[4][0]->Integral());
  hgenie_peri_preselbdt_data->Divide(hspec_mrbrem_data[4][0]);

  hgenie_peri_preselcvn_data->Scale(hspec_mrbrem_data[5][0]->Integral()/hspec_genie_fd[5][0]->Integral());
  hgenie_peri_preselcvn_data->Divide(hspec_mrbrem_data[5][0]);

  hgenie_peri_preselen_data->Scale(hspec_mrbrem_data[6][0]->Integral()/hspec_genie_fd[6][0]->Integral());
  hgenie_peri_preselen_data->Divide(hspec_mrbrem_data[6][0]);

  std::string filenameCosmicsData = "prod_mrbremcaf_R19-11-18-prod5reco.i_fd_cosmic_fhc_full_v1_goodruns_underP9";
  std::string filenameCRYMC = "prod_mrbremcaf_R19-11-18-prod5reco.m_fd_cry_fhc_full_v3";
  std::string filenameGENIEFluxswap = "prod_caf_R19-11-18-prod5reco.f_fd_genie_N1810j0211a_fluxswap_fhc_nova_v08_full_v1";

  if (mode == "rhc"){
      filenameCosmicsData = "prod_mrbremcaf_R19-11-18-prod5reco.i_fd_cosmic_rhc_full_v1_goodruns";
      filenameCRYMC = "prod_mrbremcaf_R19-11-18-prod5reco.m_fd_cry_rhc_full_v3";
      filenameGENIEFluxswap = "prod_caf_R19-11-18-prod5reco.f_fd_genie_N1810j0211a_fluxswap_rhc_nova_v08_full_v1";
  }

  //CorePresel and PeripheralPresel weights
  const Var kShwAngleCorePresel1Dweights( [hgenie_core_presel](const caf::SRProxy* sr)
                                          {
                                            double myweight = 1.0;
                                            double myShwAngle = kShwAngle(sr);
                                            int bin = hgenie_core_presel->FindBin(myShwAngle);
                                            myweight = hgenie_core_presel->GetBinContent(bin);
                                            if(myweight>0) return myweight;
                                            else return 1.0;
                                          });


  const Var kShwAnglePeripheralPresel1Dweights( [hgenie_peri_presel](const caf::SRProxy* sr)
                                                {
                                                  double myweight = 1.0;
                                                  double myShwAngle = kShwAngle(sr);
                                                  int bin = hgenie_peri_presel->FindBin(myShwAngle);
                                                  myweight = hgenie_peri_presel->GetBinContent(bin);
                                                  if(myweight>0) return myweight;
                                                  else return 1.0;
                                                });

  const Var kShwAnglePeripheralPreselBDT1Dweights( [hgenie_peri_preselbdt](const caf::SRProxy* sr)
                                                   {
                                                     double myweight = 1.0;
                                                     double myShwAngle = kShwAngle(sr);
                                                     int bin = hgenie_peri_preselbdt->FindBin(myShwAngle);
                                                     myweight = hgenie_peri_preselbdt->GetBinContent(bin);
                                                     if(myweight>0) return myweight;
                                                     else return 1.0;
                                                   });

  const Var kShwAnglePeripheralPreselCVN1Dweights( [hgenie_peri_preselcvn](const caf::SRProxy* sr)
                                                   {
                                                     double myweight = 1.0;
                                                     double myShwAngle = kShwAngle(sr);
                                                     int bin = hgenie_peri_preselcvn->FindBin(myShwAngle);
                                                     myweight = hgenie_peri_preselcvn->GetBinContent(bin);
                                                     if(myweight>0) return myweight;
                                                     else return 1.0;
                                                   });

  const Var kShwAnglePeripheralPreselEn1Dweights( [hgenie_peri_preselen](const caf::SRProxy* sr)
                                                   {
                                                     double myweight = 1.0;
                                                     double myShwAngle = kShwAngle(sr);
                                                     int bin = hgenie_peri_preselen->FindBin(myShwAngle);
                                                     myweight = hgenie_peri_preselen->GetBinContent(bin);
                                                     if(myweight>0) return myweight;
                                                     else return 1.0;
                                                   });

  const Var kShwAngleCorePreselData1Dweights( [hgenie_core_presel_data](const caf::SRProxy* sr)
                                            {
                                              double myweight = 1.0;
                                              double myShwAngle = kShwAngle(sr);
                                              int bin = hgenie_core_presel_data->FindBin(myShwAngle);
                                              myweight = hgenie_core_presel_data->GetBinContent(bin);
                                              if(myweight>0) return myweight;
                                              else return 1.0;
                                            });


  const Var kShwAnglePeripheralPreselData1Dweights( [hgenie_peri_presel_data](const caf::SRProxy* sr)
                                                  {
                                                    double myweight = 1.0;
                                                    double myShwAngle = kShwAngle(sr);
                                                    int bin = hgenie_peri_presel_data->FindBin(myShwAngle);
                                                    myweight = hgenie_peri_presel_data->GetBinContent(bin);
                                                    if(myweight>0) return myweight;
                                                    else return 1.0;
                                                  });

  const Var kShwAnglePeripheralPreselBDTData1Dweights( [hgenie_peri_preselbdt_data](const caf::SRProxy* sr)
                                                     {
                                                       double myweight = 1.0;
                                                       double myShwAngle = kShwAngle(sr);
                                                       int bin = hgenie_peri_preselbdt_data->FindBin(myShwAngle);
                                                       myweight = hgenie_peri_preselbdt_data->GetBinContent(bin);
                                                       if(myweight>0) return myweight;
                                                       else return 1.0;
                                                     });

  const Var kShwAnglePeripheralPreselCVNData1Dweights( [hgenie_peri_preselcvn_data](const caf::SRProxy* sr)
                                                     {
                                                       double myweight = 1.0;
                                                       double myShwAngle = kShwAngle(sr);
                                                       int bin = hgenie_peri_preselcvn_data->FindBin(myShwAngle);
                                                       myweight = hgenie_peri_preselcvn_data->GetBinContent(bin);
                                                       if(myweight>0) return myweight;
                                                       else return 1.0;
                                                     });

  const Var kShwAnglePeripheralPreselEnData1Dweights( [hgenie_peri_preselen_data](const caf::SRProxy* sr)
                                                     {
                                                       double myweight = 1.0;
                                                       double myShwAngle = kShwAngle(sr);
                                                       int bin = hgenie_peri_preselen_data->FindBin(myShwAngle);
                                                       myweight = hgenie_peri_preselen_data->GetBinContent(bin);
                                                       if(myweight>0) return myweight;
                                                       else return 1.0;
                                                     });

  SpectrumLoader loaderCosmicsData(filenameCosmicsData, kCosmic);
  SpectrumLoader loaderCRYMC(filenameCRYMC, kCosmic);
  SpectrumLoader loaderGENIEFluxswap(filenameGENIEFluxswap);
  loaderGENIEFluxswap.SetSpillCut(kStandardSpillCuts);
  loaderCosmicsData.SetSpillCut(kStandardSpillCuts_FD_Prod4MotivatedDQ);

  const int kNumVars2 = 14;
  const HistDef defs_2round[kNumVars2] = {
        {"shwangle", {"Cos(#theta_{Z})", Binning::Simple(60,0,1.0), kShwAngle}},
        {"cvne", {"CVNe", Binning::Simple(120, 0.0, 1.0), kCVNe_looseptp}},
        {"corebdt", {"Core BDT", Binning::Simple(120, 0.0, 1.0), kCosPIDCoreBDT}},
        {"peribdt", {"Peripheral BDT", Binning::Simple(120, 0.0, 1.0), kCosPIDPeriBDT}},
        {"shwe", {"Shower Energy (GeV)", Binning::Simple(120, 0.0, 5.0), kShwE}},
        {"shwidth", {"Shower Width (cm)", Binning::Simple(100, 0.0, 10.0), kShwWidth}},
        {"shwlen", {"Shower Length (cm)", Binning::Simple(100, 0.0, 1000.0), kShwLen}},
        {"shwstx", {"Shower Start position X (cm)", Binning::Simple(64, -800, 800), kShwStartX}},
        {"shwsty", {"Shower Start position Y (cm)", Binning::Simple(64, -800, 800), kShwStartY}},
        {"shwstz", {"Shower Start position Z (cm)", Binning::Simple(104,0,5962), kShwStartZ}},
        {"nhit", {"NHit", Binning::Simple(120,0.0,150), kNHit}},
        {"ptp", {"pTp", Binning::Simple(120,0.0,1.0), kPtP}},
        {"asym", {"SparsenessAsymm", Binning::Simple(60,-1.0,1.0), kSparsenessAsymm}},
        {"dist", {"DistNearestWall (cm)", Binning::Simple(8,0,800), kDistNearestWall}},
  };

  struct CutDef2round
  {
     std::string name;
     Cut cut;
     Var weight;
     Var weightdata;
  };

  const int kNumSels2 = 10;

  const CutDef2round sels2[kNumSels2] = {
      {"core_presel", kNue2020CorePresel_MRBrem, kShwAngleCorePresel1Dweights, kShwAngleCorePreselData1Dweights},
      {"periph_presel", kNue2020PeripheralPresel_MRBrem, kShwAnglePeripheralPresel1Dweights, kShwAnglePeripheralPreselData1Dweights},
      {"core_full", kNue2020FD_MRBrem, kShwAngleCorePresel1Dweights, kShwAngleCorePreselData1Dweights },
      {"periph_full", kNue2020FDPeripheral_MRBrem, kShwAnglePeripheralPresel1Dweights, kShwAnglePeripheralPreselData1Dweights},
      {"periph_preselbdt", kNue2020PeripheralPresel_MRBrem && kCosPIDPeriBDT > 0.57, kShwAnglePeripheralPreselBDT1Dweights, kShwAnglePeripheralPreselBDTData1Dweights},
      {"periph_fullbdt", kNue2020FDPeripheral_MRBrem, kShwAnglePeripheralPreselBDT1Dweights, kShwAnglePeripheralPreselBDTData1Dweights},
      {"periph_preselcvn", kNue2020PeripheralPresel_MRBrem && kCVNe_looseptp >= 0.97, kShwAnglePeripheralPreselCVN1Dweights, kShwAnglePeripheralPreselCVNData1Dweights},
      {"periph_fullcvn", kNue2020FDPeripheral_MRBrem, kShwAnglePeripheralPreselCVN1Dweights, kShwAnglePeripheralPreselCVNData1Dweights},
      {"periph_preselen", kNue2020PeripheralPresel_MRBrem && kShwE > 1, kShwAnglePeripheralPreselEn1Dweights, kShwAnglePeripheralPreselEnData1Dweights},
      {"periph_fullen", kNue2020FDPeripheral_MRBrem && kShwE > 1, kShwAnglePeripheralPreselEn1Dweights, kShwAnglePeripheralPreselEnData1Dweights},
  };

  Spectrum* spec_mrbrem_data_wo_weight[kNumSels2][kNumVars2];
  Spectrum* spec_mrbrem_cry_wo_weight[kNumSels2][kNumVars2];
  Spectrum* spec_genie_fd_wo_weight[kNumSels2][kNumVars2];

  Spectrum* spec_mrbrem_data_w_weight[kNumSels2][kNumVars2];
  Spectrum* spec_mrbrem_cry_w_weight[kNumSels2][kNumVars2];

  for(int selIdx = 0; selIdx < kNumSels2; ++selIdx){
      for (int varIdx=0; varIdx<kNumVars2;varIdx++){
           const HistAxis& axis = defs_2round[varIdx].axis;
           spec_mrbrem_data_w_weight[selIdx][varIdx] = new Spectrum(loaderCosmicsData, axis, sels2[selIdx].cut && kInCosmicTimingWindow_FD_MR, kNoShift, sels2[selIdx].weightdata);
           spec_mrbrem_cry_w_weight[selIdx][varIdx] = new Spectrum(loaderCRYMC, axis, sels2[selIdx].cut && kInCosmicTimingWindow_FD_MR, kNoShift, sels2[selIdx].weight);
           spec_mrbrem_data_wo_weight[selIdx][varIdx] = new Spectrum(loaderCosmicsData, axis, sels2[selIdx].cut && kInCosmicTimingWindow_FD_MR);
           spec_mrbrem_cry_wo_weight[selIdx][varIdx] = new Spectrum(loaderCRYMC, axis, sels2[selIdx].cut && kInCosmicTimingWindow_FD_MR);
           spec_genie_fd_wo_weight[selIdx][varIdx] = new Spectrum(loaderGENIEFluxswap, axis, sels2[selIdx].cut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2020);
      }
  }

  loaderCosmicsData.Go();
  loaderCRYMC.Go();
  loaderGENIEFluxswap.Go();

  TFile *file = new TFile(("MRbrem_spectra_with_weights_"+mode+".root").c_str(), "recreate");
  for(int selIdx = 0; selIdx < kNumSels2; ++selIdx){
      TDirectory* d = file->mkdir(sels2[selIdx].name.c_str());
      for (int varIdx=0; varIdx<kNumVars2;varIdx++){
           const char* name = defs_2round[varIdx].name.c_str();
           spec_mrbrem_data_wo_weight[selIdx][varIdx]->SaveTo(d, TString::Format("mrbrem_data_%s", name));
           spec_mrbrem_cry_wo_weight[selIdx][varIdx]->SaveTo(d, TString::Format("mrbrem_mc_%s", name));
           spec_genie_fd_wo_weight[selIdx][varIdx]->SaveTo(d, TString::Format("genie_fd_%s", name));
           spec_mrbrem_data_w_weight[selIdx][varIdx]->SaveTo(d, TString::Format("mrbrem_data_w_weight_%s", name));
           spec_mrbrem_cry_w_weight[selIdx][varIdx]->SaveTo(d, TString::Format("mrbrem_mc_w_weight_%s", name));
      }
  }

  file->Close();
  return;

}