NumuExtrap.C
Go to the documentation of this file.
1 #include "CAFAna/Core/Spectrum.h"
3 
4 #include "CAFAna/Cuts/Cuts.h"
9 
11 #include "CAFAna/Vars/Vars.h"
15 #include "CAFAna/Vars/XsecTunes.h"
16 
20 
21 #include "CAFAna/Core/Loaders.h"
23 #include "TFile.h"
24 
25 #include <iostream>
26 #include <cmath>
27 using namespace ana;
28 
29 void NumuExtrap(bool isFHC )
30 {
31  // What do my arguments mean?
32  std::string beam = "rhc";
33  if (isFHC) beam = "fhc";
34 
35  // Set my axis.
37 
38  // Set the default cuts to the combined quantile cuts for the Full FOM
41 
43 
44  // FD Files are the tag, so swap by "beam"
45  std::string ldrFDNonSwap = "prod_caf_R19-11-18-prod5reco.f_fd_genie_N1810j0211a_nonswap_" +beam+"_nova_v08_full_v1";
46  std::string ldrFDFlxSwap = "prod_caf_R19-11-18-prod5reco.f_fd_genie_N1810j0211a_fluxswap_"+beam+"_nova_v08_full_v1";
47  std::string ldrFDTauSwap = "prod_caf_R19-11-18-prod5reco.f_fd_genie_N1810j0211a_tau_" +beam+"_nova_v08_full_v1";
48  // ND files have different tags, so need a switch
49  std::string ldrNDNonSwap = "prod_caf_R19-11-18-prod5reco.d_nd_genie_N1810j0211a_nonswap_rhc_nova_v08_full_v1";
50  std::string ldrNDData = "prod_caf_R19-11-18-prod5reco.g_nd_numi_rhc_full_v1_goodruns";
51  if ( isFHC ){
52  ldrNDNonSwap = "prod_caf_R19-11-18-prod5reco.d.h_nd_genie_N1810j0211a_nonswap_fhc_nova_v08_full_v1";
53  ldrNDData = "prod_caf_R19-11-18-prod5reco.d.f.h_nd_numi_fhc_period235910_v1";
54  }
55 
61 
63 
64  //std::string infile_quant = (std::string)std::getenv("NUMUDATA_DIR")+"/lib/ana2020/Quantiles/quantiles_"+beam+"_full_numu2020.root";
65  std::string infile_quant = "/cvmfs/nova.opensciencegrid.org/externals/numudata/v00.03/NULL/lib/ana2020/Quantiles/quantiles_"+beam+"_full_numu2020.root";
66  TFile* infile = new TFile(infile_quant.c_str());
67 
68  TH2* FDSpec2D = (TH2*)infile->FindObjectAny("FDSpec2D_LoosePTP");
69  std::vector<Cut> HadEFracQuantCuts = QuantileCutsFromTH2(FDSpec2D, axis_numu, kHadEFracAxis, 4);
70  infile->Close();
71  delete infile;
72 
73  std::vector<IPrediction*> pred_extrap;
74  std::vector<IPrediction*> pred_noextrap;
75 
76  const Var kMyWeight = kPPFXFluxCVWgt*kXSecCVWgt2020;
77 
78  // --- Make a combined quantile Extrap Generator and NoExtrap Generator
79  NumuExtrapGenerator *decomp = new NumuExtrapGenerator(axis_numu, cut_numufd, cut_numund, kNoShift, kMyWeight);
80  pred_extrap .push_back(decomp ->Generate(loaders).release());
81 
82  NoExtrapGenerator *noextrap = new NoExtrapGenerator (axis_numu, cut_numufd, kMyWeight);
83  pred_noextrap.push_back(noextrap->Generate(loaders).release());
84 
85  // --- Loop through my Quantile cuts.
86  for(int quant=0; quant<4; ++quant) {
87 
88  Cut ThisCutND = cut_numund && HadEFracQuantCuts[quant];
89  Cut ThisCutFD = cut_numufd && HadEFracQuantCuts[quant];
90 
91  NumuExtrapGenerator *decomp_q = new NumuExtrapGenerator(axis_numu, ThisCutFD, ThisCutND, kNoShift, kMyWeight);
92  pred_extrap .push_back(decomp_q ->Generate(loaders).release());
93 
94  NoExtrapGenerator *noextrap_q = new NoExtrapGenerator (axis_numu, ThisCutFD, kMyWeight);
95  pred_noextrap.push_back(noextrap_q->Generate(loaders).release());
96  }
97 
98  loaders.Go();
99 
100  TFile* outFile = new TFile(("fd_numu_preds_extrap_"+beam+".root").c_str(),"RECREATE");
101  outFile->cd();
102  for(int quant=0; quant<5; quant++){
103  pred_extrap [quant]->SaveTo(outFile, Form("extrap_q%d" , quant));
104  pred_noextrap[quant]->SaveTo(outFile, Form("noextrap_q%d", quant));
105  }
106  outFile->Close();
107 }
108 
Near Detector underground.
Definition: SREnums.h:10
Far Detector at Ash River.
Definition: SREnums.h:11
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const Cut cut_numund
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:23
const HistAxis kHadEFracAxis("E_{had.} / E_{#nu}", Binning::Simple(200, 0, 1), kHadEFrac)
HistAxis that implements Hadronic Energy fraction binning used by L.Vinton to derive Hadronic Energy ...
Definition: HistAxes.h:30
const Cut kNumu2020ND
Definition: NumuCuts2020.h:57
void NumuExtrap(bool isFHC)
Definition: NumuExtrap.C:29
Generates FD-only predictions (no extrapolation)
Generates extrapolated Numu predictions.
const Cut cut_numufd
void Go()
Call Go() on all the loaders.
Definition: Loaders.cxx:162
void SetSpillCut(const SpillCut &cut)
Definition: Loaders.cxx:121
TFile * outFile
Definition: PlotXSec.C:135
string infile
const Cut kNumu2020FD
Definition: NumuCuts2020.h:59
const HistAxis kNumuCCOptimisedAxis("Reconstructed Neutrino Energy (GeV)", kNumuCCEOptimisedBinning, kCCE)
HistAxis that implements optimised numuCCE from L. Vinton. See docdb 16332. This was close to &#39;custC&#39;...
Definition: HistAxes.h:24
static bool isFHC
const SystShifts kNoShift
Definition: SystShifts.cxx:22
std::unique_ptr< IPrediction > Generate(Loaders &loaders, const SystShifts &shiftMC=kNoShift) const override
const HistAxis axis_numu
std::vector< Loaders * > loaders
Definition: syst_header.h:386
std::vector< Cut > QuantileCutsFromTH2(TH2 *quantileHist, const HistAxis &independentAxis, const HistAxis &quantileAxis, const unsigned int &numQuantiles, const bool verbose)
: Do the same as the QuantileCuts function but taking in the TH2 instead of making it...
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
void SetLoaderPath(const std::string &path, caf::Det_t det, DataMC datamc, DataSource src=kBeam, SwappingConfig swap=kNonSwap)
Configure loader via wildcard path.
Definition: Loaders.cxx:25
const Var kXSecCVWgt2020
Definition: XsecTunes.h:105
enum BeamMode string