cafe_FD_predictions.C
Go to the documentation of this file.
1 // This macro creates unextrapolated predictions at the far detector for the
2 // nue selections (low, high, peripheral) and numu selections (Q1-Q4). Both
3 // unoscillated and oscillated predictions are created and written into a ROOT file.
4 //
5 // The predictions from CAF are loaded up from SAM datasets to avoid having weights
6 // from a particular GENIE version baked into them.
7 
8 #include <string>
9 #include <vector>
10 #include <set>
11 
12 #include "TFile.h"
13 #include "TH2.h"
14 
15 #include "CAFAna/Analysis/Calcs.h"
19 #include "CAFAna/Vars/FitVars.h"
20 #include "CAFAna/numu/Analysis2018/includes.h"
23 #include "CAFAna/Vars/NueEnergy2018.h"
24 #include "NuXAna/Vars/NusVars.h"
27 
28 #include "Utilities/rootlogon.C"
29 #include "OscLib/IOscCalc.h"
30 
31 using namespace ana;
32 
33 // convenience enumeration for each selection type
34 // The CAFAna prediction files put all three nue selections
35 // into a single histogram
36 typedef enum SelType{
37  lNue = 0,
42 } SelType_t;
43 
44 std::vector<std::string> gSelTypeNames({"NuE", "NuMuQ1", "NuMuQ2", "NuMuQ3", "NuMuQ4"});
45 
46 //-----------------------------------------------------------------------------
47 // convenience struct to hold information about the selection, horn current,
48 // and prediction file name.
50 
51  PredictionInfo(std::string const& hornCur,
52  SelType_t const& selType)
53  : curName (hornCur)
54  , selection(selType)
55  {
56  fhc = (curName.find("FHC") != std::string::npos) ? true : false;
57  std::string recolet("d");
58  std::string lowerCur("fhc");
59  std::string lowerSel("numu");
60 
61  if(fhc){
63  }
64  else{
65  recolet = "e";
66  lowerCur = "rhc";
68  }
69 
70  // define cuts and vars we will need
72  Var energy = kCCE;
74 
75  if(selection == lNue){
76  bins = kNue2018Binning;
77  cuts = {kNue2018FDAllSamples && kHasNeutrino && kInBeamSpill};
78  energy = kNue2018AnaBin;
79  lowerSel = "nue";
80  }
81  else{
82  std::string quantFile = "/cvmfs/nova.opensciencegrid.org/externals/numudata/v00.00/NULL/ana2018/Quantiles/";
83  quantFile += (fhc) ? "quantiles__fhc_full__numu2018.root" : "quantiles__rhc_full__numu2018.root";
84  TFile qtf(quantFile.c_str(), "READ");
85  TH2F *FDSpec2D = dynamic_cast<TH2F*>(qtf.Get("FDSpec2D"));
86  std::vector<Cut> numuQuantCuts = QuantileCutsFromTH2(FDSpec2D, kNumuCCOptimisedAxis, kHadEFracAxis, 4);
87 
88  if(selection == lNumuQ1){
89  cuts = {kNumuCutFD2018 && kHasNeutrino && kInBeamSpill && numuQuantCuts[0]};
90  }
91  else if(selection == lNumuQ2){
92  cuts = {kNumuCutFD2018 && kHasNeutrino && kInBeamSpill && numuQuantCuts[1]};
93  }
94  else if(selection == lNumuQ3){
95  cuts = {kNumuCutFD2018 && kHasNeutrino && kInBeamSpill && numuQuantCuts[2]};
96  }
97  else if(selection == lNumuQ4){
98  cuts = {kNumuCutFD2018 && kHasNeutrino && kInBeamSpill && numuQuantCuts[3]};
99  }
100  }
101 
102  // create Spectrum loaders from the concat datasets
103  // TODO: no cosmic backgrounds or rock backgrounds can be loaded in this way
104  std::string base = "prod_sumdecaf_R17-11-14-prod4reco." + recolet + "_fd_genie_";
105  SpectrumLoader nonSwapLoader (base + "nonswap_" + lowerCur + "_nova_v08_full_v1_" + lowerSel + "2018");
106  SpectrumLoader fluxSwapLoader(base + "fluxswap_" + lowerCur + "_nova_v08_full_v1_" + lowerSel + "2018");
107  SpectrumLoader tauSwapLoader (base + "tau_" + lowerCur + "_nova_v08_full_v1_" + lowerSel + "2018");
108 
109  // now make the no extrapolation prediction thing
110  cafPred = new PredictionNoExtrap(nonSwapLoader,
111  fluxSwapLoader,
112  tauSwapLoader,
113  "E_{reco} (GeV)",
114  bins,
115  energy,
116  cuts,
117  kNoShift,
119 
120  // apparently we make the loaders do stuff after we give them to the prediction...
121  nonSwapLoader .Go();
122  fluxSwapLoader.Go();
123  tauSwapLoader .Go();
124  }
125 
126  //---------------------------------------------------------------
128  {
129  std::string totName = (gSelTypeNames[selection] + curName);
130  totName += (osc) ? "Oscillated" : "Unoscillated";
131  return totName;
132  }
133 
134  //---------------------------------------------------------------
137  {
138  TH1D* hist = nullptr;
139  if(osc){
140  hist = cafPred->Predict(calc).ToTH1(pot);
141  }
142  else
143  hist = cafPred->PredictUnoscillated().ToTH1(pot);
144 
145  // set these histograms free from the input file by calling
146  // SetDirectory(0)
147  hist->SetDirectory(0);
148 
149  // the bool indicates if the histograms are oscillated or not
150  hist->SetName((HistName(osc)).c_str());
151 
152  return hist;
153  }
154 
155  //---------------------------------------------------------------
156 
157  bool fhc; ///< true for FHC, false for RHC
158  double pot; ///< total POT for horn current
159  SelType_t selection; ///< numu + quantile or nue
160  std::string curName; ///< string horn current value
161  PredictionNoExtrap* cafPred; ///< CAFAna prediction object
162 };
163 
164 //-----------------------------------------------------------------------------
165 // Set the oscillation parameters. Perhaps we want to set these on the command
166 // line at some point
168 {
169  calc->SetdCP(0.17);
170  calc->SetTh13(0.145426);
171  calc->SetTh23(0.865744);
172  calc->SetDmsq32(2.51e-3);
173 }
174 
175 //-----------------------------------------------------------------------------
176 // main function to create the prediction histograms
177 // bool fhc: true if we want FHC events
178 // int sel: the selection type according to SelType_t defined above
179 void cafe_FD_predictions(bool fhc,
180  int sel)
181 {
182  std::string curStr = (fhc) ? "FHC" : "RHC";
183 
184  // prediction information for the desired selection type
185  // and horn current.
186  PredictionInfo pi(curStr, (SelType_t)sel);
187 
188  // Set up the oscillation calculator
191 
192  // open the output file.
193  std::string outFileName("caf_fd_predictions_");
194  outFileName += (fhc) ? "fhc_" : "rhc_";
195  outFileName += gSelTypeNames[sel];
196  outFileName += ".root";
197  TFile* outfile = new TFile(outFileName.c_str(), "RECREATE");
198 
199  // get and write prediction histograms for the oscillated and
200  // unoscillated cases
201 
202  auto oscHist = pi.SpectrumFromPrediction(true, calc);
203  oscHist->SetDirectory(outfile);
204  oscHist->Write();
205 
206  auto noOscHist = pi.SpectrumFromPrediction(false, calc);
207  noOscHist->SetDirectory(outfile);
208  noOscHist->Write();
209 
210  // write and close the output file.
211  outfile->Write();
212  outfile->Close();
213 }
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
PredictionNoExtrap * cafPred
CAFAna prediction object.
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
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
virtual void SetTh13(const T &th13)=0
osc::OscCalcDumb calc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
const XML_Char int const XML_Char int const XML_Char * base
Definition: expat.h:331
virtual void SetDmsq32(const T &dmsq32)=0
const Var kNue2018AnaBin([](const caf::SRProxy *sr){int selBin=kNue2018SelectionBin(sr);float nuE=kNueEnergy2018(sr);int nuEBin=nuE/0.5;assert(nuEBin<=8 &&"An event with nuE > 4.5 should never happen");int anaBin=9 *selBin+nuEBin;return anaBin;})
Use this Analysis Binning for Ana2018, official Binning.
Definition: NueCuts2018.h:171
std::string HistName(bool osc)
enum SelType SelType_t
std::vector< std::string > gSelTypeNames({"NuE","NuMuQ1","NuMuQ2","NuMuQ3","NuMuQ4"})
const Var kCCE
Definition: NumuVars.h:21
double energy
Definition: plottest35.C:25
#define pot
double pot
total POT for horn current
virtual void Go() override
Load all the registered spectra.
const Binning kNue2018Binning
Definition: NueCuts2018.cxx:93
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
Oscillation probability calculators.
Definition: Calcs.h:5
void SetOscillationParameters(osc::IOscCalcAdjustable *calc)
const SystShifts kNoShift
Definition: SystShifts.cxx:22
const Cut kHasNeutrino([](const caf::SRProxy *sr){return(sr->mc.nnu!=0);})
Check if MC slice has neutrino information (useful for in-and-out tests)
Definition: TruthCuts.h:61
const double kAna2019RHCPOT
Definition: Exposures.h:224
const Binning bins
Definition: NumuCC_CPiBin.h:8
void cafe_FD_predictions(bool fhc, int sel)
SelType_t selection
numu + quantile or nue
std::string curName
string horn current value
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const Var kXSecCVWgt2018
Definition: XsecTunes.h:48
virtual void SetTh23(const T &th23)=0
TH1D * SpectrumFromPrediction(bool osc, osc::IOscCalcAdjustable *calc)
const Binning kNumuCCEOptimisedBinning
Optimised binning for numuCCE from L. Vinton. See docdb 16332. This was close to &#39;custC&#39; in that talk...
Definition: Binnings.cxx:28
const Cut kNue2018FDAllSamples
Definition: NueCuts2018.h:77
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 Cut kInBeamSpill([](const caf::SRProxy *sr){if(sr->spill.run > util::kLastBadTimingRun) return kInBeamSpill_main(sr);else return kInBeamSpill_main(sr)||kInBeamSpill_shifted(sr);}, [](const caf::SRSpillProxy *spill){if(spill->run > util::kLastBadTimingRun) return kInBeamSpill_main.Livetime(spill);else return kInBeamSpill_main.Livetime(spill)+kInBeamSpill_shifted.Livetime(spill);}, [](const caf::SRSpillProxy *spill) -> double{return spill->spillpot;})
Does the event fall inside the window we call the beam spill?
Definition: TimingCuts.h:8
const double kAna2019FHCPOT
Definition: Exposures.h:223
const Cut kNumuCutFD2018
Definition: NumuCuts2018.h:39
Prediction that just uses FD MC, with no extrapolation.
Float_t e
Definition: plot.C:35
bool fhc
true for FHC, false for RHC
PredictionInfo(std::string const &hornCur, SelType_t const &selType)
FILE * outfile
Definition: dump_event.C:13
virtual void SetdCP(const T &dCP)=0
enum BeamMode string