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