make_nueFDprediction_kinematics_RHC_REW.C

Go to the documentation of this file.

// Analysis
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Analysis/Prod5Loaders.h"
#include "CAFAna/Analysis/Style.h"

// Core
#include "CAFAna/Core/Binning.h"
#include "CAFAna/Core/Cut.h"
#include "CAFAna/Core/HistAxis.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/SpectrumLoaderBase.h"

// Cuts
#include "CAFAna/Cuts/Cuts.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "3FlavorAna/Cuts/NumuCuts2020.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "3FlavorAna/Cuts/NueCuts2020.h"

// Decomp
#include "CAFAna/Decomp/ProportionalDecomp.h"
#include "3FlavorAna/Decomp/BENDecomp.h"
#include "3FlavorAna/Decomp/MichelDecomp.h"
#include "3FlavorAna/Decomp/NumuDecomp.h"

// Extrap
#include "CAFAna/Extrap/ModularExtrap.h"

// Predict
#include "CAFAna/Prediction/PredictionExtrap.h"
#include "3FlavorAna/Prediction/PredictionGenerator3Flavor.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "3FlavorAna/Prediction/PredictionExtendToPeripheral.h"

// Systs
#include "CAFAna/Systs/Systs.h"

// Vars
#include "CAFAna/Vars/XsecTunes.h"
#include "CAFAna/Vars/TruthVars.h"
#include "3FlavorAna/Vars/HistAxes.h"
#include "3FlavorAna/Vars/NueVars.h"
#include "CAFAna/Vars/Vars.h"
#include "3FlavorAna/Vars/NueVarsExtra.h"
#include "CAFAna/Vars/PPFXWeights.h"

#include "OscLib/OscCalcDumb.h"
#include "StandardRecord/Proxy/SRProxy.h"

#include <iostream>
#include <iomanip>

using namespace ana;

osc::OscCalcDumb calc;

//sample : ALL/CORE/PERI and the corresponding Weight file name. 
void make_nueFDprediction_kinematics_RHC_REW(const std::string sample, const std::string wfile)

{  
  const std::string outfilename = ("FDprediction_kinematics_RHC_REW_"+sample+".root").c_str();

  //For grid
  //  std::string wfile_name = ("/pnfs/nova/scratch/users/zvallari/AcceptSyst/"+wfile).c_str();
  //  TFile *weight_file = TFile::Open(pnfs2xrootd(wfile_name).c_str() );

  //For non-grid
  std::string wfile_name = ("/nova/ana/users/zvallari/AcceptSyst2020/"+wfile).c_str();
  TFile *weight_file = TFile::Open(wfile_name.c_str()); 

  // construction of trueQ2 weight
  TH1 *h;
  gDirectory->GetObject(("trueQ2_weight_"+sample).c_str(),h);
    
  std::cout << h->GetNbinsX() << endl;
  int nbins = h->GetNbinsX();
  double bins[nbins], binslow[nbins], binshigh[nbins];
  for(int i = 1; i <= nbins; ++i){
    bins[i-1] = h->GetBinContent(i);
    binslow[i-1] = h->GetBinLowEdge(i);
    binshigh[i-1] = h->GetBinLowEdge(i) + h->GetBinWidth(i);
  }
  
  const Var kTrueQ2Weight([nbins, &binslow, &binshigh, &bins](const caf::SRProxy *sr){
      if(sr->hdr.det == caf::kFARDET) return 1.0;
      if(!kNumu2020ND(sr)) {
        std::cout << "THERE IS AN EVENT WHICH DOES NOT PASS NUMU ND CUT IN NUE SIGNAL EXTRAP!!!!" <<  std::endl;
        exit(0);
        return 1.0;
      }
      for(int i = 0; i < nbins; ++i){
        if((kRecoQ2(sr) > binslow[i]) && (kRecoQ2(sr) < binshigh[i])){
          return bins[i];         
        }
      }
      return 1.0;
    }
    );
  
  //construction of PtP weight
  TH1 *hptp;
  gDirectory->GetObject(("PtP_weight_"+sample).c_str(),hptp);
  
  std::cout << hptp->GetNbinsX() << endl;
  int nbins_ptp = hptp->GetNbinsX();
  double bins_ptp[nbins_ptp], binslow_ptp[nbins_ptp], binshigh_ptp[nbins_ptp]; 
  for(int i = 1; i <= nbins_ptp; ++i){
    bins_ptp[i-1] = hptp->GetBinContent(i);
    binslow_ptp[i-1] = hptp->GetBinLowEdge(i);
    binshigh_ptp[i-1] = hptp->GetBinLowEdge(i) + hptp->GetBinWidth(i);
  }
  
  const Var kPtPWeight([nbins_ptp, &binslow_ptp, &binshigh_ptp, &bins_ptp](const caf::SRProxy *sr){
      if(sr->hdr.det == caf::kFARDET) return 1.0;
      for(int i = 0; i < nbins_ptp; ++i){
        if((kPtP(sr) > binslow_ptp[i]) && (kPtP(sr) < binshigh_ptp[i])){return bins_ptp[i];}
      }
      return 1.0;
    }
    );
  
  //construction of Cos weight
  TH1 *hcos;
  gDirectory->GetObject(("CosNumi_weight_"+sample).c_str(),hcos);
  
  std::cout << hcos->GetNbinsX() << endl;
  int nbins_cos = hptp->GetNbinsX();
  double bins_cos[nbins_cos], binslow_cos[nbins_cos], binshigh_cos[nbins_cos]; 
  for(int i = 1; i <= nbins_cos; ++i){
    bins_cos[i-1] = hcos->GetBinContent(i);
    binslow_cos[i-1] = hcos->GetBinLowEdge(i);
        binshigh_cos[i-1] = hcos->GetBinLowEdge(i) + hcos->GetBinWidth(i);
    }
  
  const Var kCosWeight([nbins_cos, &binslow_cos, &binshigh_cos, &bins_cos](const caf::SRProxy *sr){
      if(sr->hdr.det == caf::kFARDET) return 1.0;
      for(int i = 0; i < nbins_cos; ++i){
        if((kCosNumi(sr) > binslow_cos[i]) && (kCosNumi(sr) < binshigh_cos[i])){return bins_cos[i];}
      }
      return 1.0;
    }
    );
  
  // close weights file
  weight_file->Close();

  Prod5NomLoaders loaders(ana::kNueConcat, ana::Loaders::kRHC);
  loaders.SetSpillCut(kStandardSpillCuts);

  struct GenDef{
    const IPredictionGenerator *gen;
    const TString cutname;
    const TString varname;
  };

  std::vector <GenDef> gens;
  std::vector <IPrediction * > predictions;

  gens.push_back({new NoExtrapGenerator(kNue2020Axis, kNue2020FDAllSamples, kXSecCVWgt2020*kPPFXFluxCVWgt),
        sample.c_str(), "nueAxis_NoExtrap"}); 

      
  gens.push_back({new NueSignalExtrapGenerator( kNue2020Axis,
                                                kNumuCCOptimisedAxis2020,
                                                kNue2020FDAllSamples,
                                                kNumu2020ND,
                                                kNoShift,
                                                kXSecCVWgt2020*kPPFXFluxCVWgt*kTrueQ2Weight),
        sample.c_str(), "nueAxis_NueSignalExtrap_Q2Weight"});

  gens.push_back({new NueSignalExtrapGenerator( kNue2020Axis,
                                                kNumuCCOptimisedAxis2020,
                                                kNue2020FDAllSamples,
                                                kNumu2020ND,
                                                kNoShift,
                                                kXSecCVWgt2020*kPPFXFluxCVWgt*kCosWeight),
        sample.c_str(), "nueAxis_NueSignalExtrap_CosWeight"});

  gens.push_back({new NueSignalExtrapGenerator( kNue2020Axis,
                                                kNumuCCOptimisedAxis2020,
                                                kNue2020FDAllSamples,
                                                kNumu2020ND,
                                                kNoShift,
                                                kXSecCVWgt2020*kPPFXFluxCVWgt*kPtPWeight),
        sample.c_str(), "nueAxis_NueSignalExtrap_PtPWeight"});
  

  for(auto & gen:gens){
      predictions.push_back(gen.gen->Generate(loaders).release());
  }

  loaders.Go();

  TFile* file = new TFile(outfilename.c_str(),"RECREATE");
  for (int i = 0 ; i < (int)gens.size(); ++i){
    auto dir = file->GetDirectory(gens[i].cutname);
    if(!dir) dir = file->mkdir(gens[i].cutname);
    predictions[i]->SaveTo(dir, gens[i].varname);
  }

  file->Close();
  
}