MakeCovarSim.C
Go to the documentation of this file.
1 /// MakeCovarSim.C by J. Hewes <jhewes15@fnal.gov>
2 
5 
11 
12 #include "CAFAna/Cuts/SpillCuts.h"
15 #include "NuXAna/Cuts/NusCuts18.h"
17 
18 #include "NuXAna/Vars/NusVars.h"
20 #include "CAFAna/Vars/HistAxes.h"
23 #include "NuXAna/Vars/HistAxes.h"
24 
27 
28 #include "CAFAna/Analysis/Calcs.h"
31 #include "CAFAna/Core/Sample.h"
32 
34 
35 #include "OscLib/OscCalcSterile.h"
36 
37 using namespace ana;
38 
39 void MakeCovarSim(TString opt)
40 {
41  DontAddDirectory guard;
42 
43  // Set detector and analysis type
44 
45  bool nd;
46  covmx::Detector c_det;
47  if (CheckOption(opt, "neardet")) {
48  nd = true;
50  }
51  else if (CheckOption(opt, "fardet")) {
52  nd = false;
54  }
55  else assert(false and "Option string must include detector: \"neardet\" or \"fardet\".");
56  std::string detector = nd ? "neardet" : "fardet";
57 
59  covmx::Analysis c_ana;
60  if (CheckOption(opt, "nue")) {
61  ana = "nue";
63  }
64  else if (CheckOption(opt, "numu")) {
65  ana = "numu";
67  }
68  else if (CheckOption(opt, "nus")) {
69  ana = "nus";
70  c_ana = covmx::Analysis::kNC;
71  }
72  else assert(false and "Option string must include analysis: \"nue\", \"numu\" or \"nus\".");
73 
74  bool rhc;
75  covmx::Polarity c_pol;
76  if (CheckOption(opt, "fhc")) {
77  rhc = false;
78  c_pol = covmx::Polarity::kFHC;
79  }
80  else if (CheckOption(opt, "rhc")) {
81  rhc = true;
82  c_pol = covmx::Polarity::kRHC;
83  }
84  else assert(false and "Option string must include beam polarity: \"fhc\" or \"rhc\".");
85  std::string polarity = rhc ? "rhc" : "fhc";
86 
87  covmx::Sample sample(c_ana, c_pol, c_det);
88 
89  bool concat;
90  if (CheckOption(opt, "fullcaf")) concat = false;
91  else if (CheckOption(opt, "concat")) concat = true;
92  else assert(false and "Option string must include file type: \"concat\" or \"fullcaf\".");
93 
94  std::string syst_type;
95  if (CheckOption(opt, "nosysts")) syst_type = "nosysts";
96  else if (CheckOption(opt, "filesysts")) syst_type = "systs";
97  else if (CheckOption(opt, "geniesysts")) syst_type = "genie";
98  else assert(false and "Option string must include syst type: \"nosysts\", \"filesysts\" or \"geniesysts\".");
99 
100  std::string optstring = std::string(opt.Data());
101 
102  size_t quantile = 999;
103  if (ana == "numu") {
104  size_t qpos = optstring.find("quantile");
105  if (qpos != std::string::npos) {
106  if (optstring.substr(qpos+8, 1) != "=")
107  assert(false and "The quantile option must take the form \"quantile=xxx\".");
108  size_t start = qpos+9;
109  size_t end = optstring.find("_", qpos);
110  if (end != std::string::npos)
111  quantile = std::stod(optstring.substr(start, end-start));
112  else
113  quantile = std::stod(optstring.substr(start));
114  std::cout << "Using quantile value " << quantile << std::endl;
115  }
116  }
117 
118  size_t usesyst = 999;
119  size_t spos = optstring.find("usesyst");
120  if (spos != std::string::npos) {
121  if (optstring.substr(spos+7, 1) != "=")
122  assert(false and "The usesyst option must take the form \"usesyst=xxx\".");
123  size_t start = spos+8;
124  size_t end = optstring.find("_", spos);
125  if (end != std::string::npos)
126  usesyst = std::stod(optstring.substr(start, end-start));
127  else
128  usesyst = std::stod(optstring.substr(start));
129  std::cout << "Using systematic in position " << usesyst << std::endl;
130  }
131  // else assert(false && "You must provide a \"usesyst\" option!");
132 
133  // Important variables
134  const Var kReweight = kXSecCVWgt2018 * kPPFXFluxCVWgt;
136 
137  // Set up systematics
138  std::vector<const ISyst*> all_systs, systs;
139  if (syst_type == "systs")
140  all_systs = GetSystsFromFile(sample);
141  if (syst_type == "genie")
142  all_systs = getAllXsecSysts_2018();
143 
144  if (all_systs.size() > 0) {
145 
146  if (usesyst == 999) {
147  std::cout << "No \"usesyst\" option specified, so running for all systs simultaneously!" << std::endl;
148  systs = all_systs;
149  }
150 
151  else if (usesyst >= all_systs.size()) {
152  std::ostringstream err;
153  err << "You asked for the systematic in position "
154  << usesyst << ", but there are only "
155  << all_systs.size() << " systematics in the vector provided.";
156  assert(false and err.str().c_str());
157  }
158 
159  else {
160  const ISyst* syst = all_systs[usesyst];
161  std::cout << "Running for systematic " << syst->ShortName() << std::endl;
162  systs = { syst };
163  }
164 
165  } // if allsysts not empty
166 
167  // Set CAF type
168  enum ECAFType caf_type;
169  if (!concat) caf_type = ECAFType::kFullCAF;
170  else if (ana == "nue") caf_type = ECAFType::kNueConcat;
171  else if (ana == "numu") caf_type = ECAFType::kNumuConcat;
172  else if (ana == "nus") caf_type = ECAFType::kNusConcat;
173  else assert(false and "Couldn't find valid caf type.");
174 
175  // List systematics
176  // std::cout << std::endl << "Enabled systematics:" << std::endl;
177  // for (auto syst : systs)
178  // std::cout << " " << syst->ShortName() << std::endl;
179  // std::cout << std::endl;
180 
181  // Set axes and syst type
182  const HistAxis* axis = nullptr;
183  const Cut* cut = nullptr;
184 
185  // Nue axis and cuts
186  if (ana == "nue") {
187  axis = new const HistAxis(kNue2018AxisMergedPeripheral);
188  if (nd) cut = new const Cut(kNue2018NDCVNSsb);
189  else cut = new const Cut(kNue2018FDAllSamples);
190  }
191 
192  // Numu axis and cuts
193  else if (ana == "numu") {
194  axis = new const HistAxis(kNumuCCOptimisedAxis);
195  if (quantile > 3)
196  assert(false and "For a numu prediction you must define a quantile between 0 and 3!");
197  std::string quant_dir = "/pnfs/nova/persistent/users/jhewes15/quantiles";
198  std::string infile_quant = quant_dir + "/quantiles__" + polarity + "_full__numu2018.root";
199  std::cout << "\n\n --- get quantile cuts from file ---" << std::endl;
200  TFile* infile = TFile::Open( pnfs2xrootd(infile_quant).c_str() );
201  TH2* FDSpec2D = (TH2*)infile->FindObjectAny("FDSpec2D");
202  std::vector<Cut> HadEFracQuantCuts = QuantileCutsFromTH2(FDSpec2D, kNumuCCOptimisedAxis, kHadEFracAxis, 4);
203  infile->Close();
204  delete infile;
205  if (nd) cut = new const Cut(kNumuCutND2018 && HadEFracQuantCuts[quantile]);
206  else cut = new const Cut(kNumuCutFD2018 && HadEFracQuantCuts[quantile]);
207  }
208 
209  // Neutral current axis and cuts
210 
211  else if (ana == "nus") {
212  if (nd) axis = new const HistAxis(kNCNDAxisE);
213  else axis = new const HistAxis(kNCFDAxisE);
214  if (nd) cut = new const Cut(kNus18ND);
215  else cut = new const Cut(kNus18FD);
216  }
217 
218  // Define the loader
219  Prod4NomLoaders loaders(caf_type,
221  if (nd) loaders.SetND(true);
223 
224  // Create the prediction
225  IPredictionGenerator* gen = nullptr;
226  if (nd) gen = new NDPredictionGenerator(*axis, *cut, kNoShift, kReweight);
227  else gen = new FDPredictionGenerator(*axis, *cut, kNoShift, kReweight);
228  IPrediction* pred = gen->Generate(loaders).release();
229  PredictionInterp* pred_interp = nullptr;
230  if (systs.size() > 0)
231  pred_interp = new PredictionInterp(systs, calc, *gen, loaders);
232  loaders.Go();
233 
234  // Initialise output file
235  std::string systname;
236  if (systs.size() == 1)
237  systname = systs[0]->ShortName();
238  else if (systs.size() == 0)
239  systname = "nosysts";
240  else
241  systname = syst_type;
242 
243  std::ostringstream f_out_name;
244  f_out_name << "pred_" << ana << "_" << polarity << "_" << detector;
245  if (ana == "numu") f_out_name << "_q" << quantile;
246  f_out_name << "_" << systname << ".root";
247  TFile* f_out = TFile::Open(f_out_name.str().c_str(), "recreate");
248 
249  // Save prediction to file
250  std::ostringstream dirname;
251  dirname << "pred_" << ana << "_" << polarity << "_" << detector;
252  if (ana == "numu") dirname << "_q" << quantile;
253  dirname << "_" << systname;
254 
255  pred->SaveTo(f_out, dirname.str().c_str());
256 
257  if (pred_interp) {
258  dirname.str("");
259  dirname << "pred_interp_" << ana << "_" << polarity << "_" << detector;
260  if (ana == "numu") dirname << "_q" << quantile;
261  dirname << "_" << systname;
262  pred_interp->SaveTo(f_out, dirname.str().c_str());
263  }
264 
265 } // function MakeCovarSim
266 
void MakeCovarSim(TString opt)
Definition: MakeCovarSim.C:39
Implements systematic errors by interpolation between shifted templates.
_HistAxis< Var > HistAxis
Definition: HistAxis.h:103
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::vector< const ISyst * > getAllXsecSysts_2018()
Get master XSec syst list for 2018 analyses.
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
double stod(const std::string &val)
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
virtual const std::string & ShortName() const final
The name printed out to the screen.
Definition: ISyst.h:27
const Cut kNumuCutND2018
Definition: NumuCuts2018.h:41
const Cut kNue2018NDCVNSsb
Definition: NueCuts2018.h:153
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
Adapt the PMNS_Sterile calculator to standard interface.
virtual void SaveTo(TDirectory *dir, const std::string &name) const
std::string pnfs2xrootd(std::string loc, bool unauth)
Definition: UtilsExt.cxx:237
osc::OscCalcDumb calc
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
const Cut kNus18FD
Definition: NusCuts18.h:100
void Go()
Call Go() on all the loaders.
Definition: Loaders.cxx:162
void SetSpillCut(const SpillCut &cut)
Definition: Loaders.cxx:121
virtual void SaveTo(TDirectory *dir, const std::string &name) const override
osc::OscCalcSterile * DefaultSterileCalc(int nflavors)
Create a sterile calculator with default assumptions for all parameters.
Definition: Calcs.cxx:97
const HistAxis kNCNDAxisE("Energy Deposited in Scintillator (GeV)", kNCNDBinning, kNus18Energy)
NC covariance hist axes.
Definition: HistAxes.h:20
string infile
_Cut< caf::SRProxy > Cut
Representation of a cut (selection) to be applied to a caf::StandardRecord object.
Definition: Cut.h:96
const HistAxis kNue2018AxisMergedPeripheral("NuE Energy / Analysis Bin", kNue2018BinningMergedPeripheral, kNue2018AnaBinMergedPeripheral)
Definition: NueCuts2018.h:181
const Cut kNus18ND
Full Nus18 ND analysis selection.
Definition: NusCuts18.h:137
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
bool CheckOption(TString opts, TString opt)
unsigned int quantile(std::vector< std::vector< double >> v, float E, float hadEfrac)
Definition: Toy_analyses.C:480
const SystShifts kNoShift
Definition: SystShifts.cxx:22
OStream cout
Definition: OStream.cxx:6
const Cut cut
Definition: exporter_fd.C:30
const HistAxis kNCFDAxisE("Energy Deposited in Scintillator (GeV)", kNCFDBinning, kNus18Energy)
Definition: HistAxes.h:21
virtual std::unique_ptr< IPrediction > Generate(Loaders &loaders, const SystShifts &shiftMC=kNoShift) const =0
const Var kXSecCVWgt2018
Definition: XsecTunes.h:49
std::vector< Loaders * > loaders
Definition: syst_header.h:386
void SetND(bool nd)
Definition: Loaders.h:65
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...
assert(nhit_max >=nhit_nbins)
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
const Cut kNumuCutFD2018
Definition: NumuCuts2018.h:39
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
Given loaders and an MC shift, Generate() generates an IPrediction.
For nominal spectra and reweighting systs (xsec/flux)
Definition: Prod4Loaders.h:96
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
Generates Near Detector predictions.
std::vector< const ISyst * > GetSystsFromFile(covmx::Sample sample)
Get systematics for a given sample.
std::vector< std::string > all_systs
ECAFType
Definition: Loaders.h:19
enum BeamMode string