Nus17Systs.cxx
Go to the documentation of this file.
2 
3 #include "CAFAna/Systs/Systs.h"
4 #include "CAFAna/Vars/Vars.h"
5 #include "NuXAna/Vars/NusVars.h"
7 #include "Utilities/func/MathUtil.h"
8 
9 #include "TFile.h"
10 #include "TH1.h"
11 
12 #include <cassert>
13 #include <cmath>
14 
15 
16 namespace ana
17 {
18  //const Nus17SystFromHist kNus17OscParamSysts (kNusAna01SystFAlt, "EX", "OscParam", "Oscillation Parameters");
19 /* const Nus17SystFromHist kNus17BeamSysts (kNusAna17SystFile, "EX", "Beam17", "All Beam");
20  const Nus17SystFromHist kNus17CherenkovSyst (kNusAna17SystFile, "EX", "Cherenkov17", "Cherenkov");
21  const Nus17SystFromHist kNus17Cherenkov_mSyst (kNusAna17SystFile, "EX", "Cherenkov_m17", "Light level method");
22  const Nus17SystFromHist kNus17CalibFlatSyst (kNusAna17SystFile, "EX", "CalFlat17", "Flat Miscalibration");
23  const Nus17SystFromHist kNus17CalibRelFDSyst (kNusAna17SystFile, "EX", "CalRelFD17", "Relative Detector FD Miscalibration");
24  const Nus17SystFromHist kNus17CalibShapeSyst (kNusAna17SystFile, "EX", "CalShape17", "Calibration shape");
25  const Nus17SystFromHist kNus17GENIESmallSysts (kNusAna17SystFile, "EX", "GENIESm17", "Summed small GENIE Systs");
26  const Nus17SystFromHist kNus17LightLSyst (kNusAna17SystFile, "EX", "CalLight17", "Lightlevel syst");
27  const Nus17SystFromHist kNus17NueCCSyst (kNusAna17SystFile, "EX", "NueCC17", "#nu_{e} CC Background");
28  const Nus17SystFromHist kNus17NumuCCSyst (kNusAna17SystFile, "EX", "NumuCC17", "#nu_{#mu} CC Background");
29  const Nus17SystFromHist kNus17NCSyst (kNusAna17SystFile, "EX", "NC17", "NC Background");
30  const Nus17SystFromHist kNus17AccSyst (kNusAna17SystFile, "EX", "NDAcc17", "NC Background"); //Caveat: Latex name is "acceptance syst". To match with tat of ISyst object created, the wrong Latex name is retained. But it is not troublesome. Will update soon.
31  const Nus17SystFromHist kNus17CalibRelNDSyst (kNusAna17SystFile, "EX", "CalRelND17", "Relative Detector NDMiscalibration");
32 */
33 
34 /* const Nus17SystFromHist kNus17NDBeamSysts (kNusAna17SystFileND, "ND", "NDBeam17", "All Beam");
35  const Nus17SystFromHist kNus17NDCherenkovSyst (kNusAna17SystFileND, "ND", "NDCherenkov17", "Cherenkov");
36  const Nus17SystFromHist kNus17NDCalibFlatSyst (kNusAna17SystFileND, "ND", "NDCalFlat17", "Flat Miscalibration");
37  const Nus17SystFromHist kNus17NDCalibShapeSyst (kNusAna17SystFileND, "ND", "NDCalShape17", "Calibration shape");
38  const Nus17SystFromHist kNus17NDGENIESmallSysts (kNusAna17SystFileND, "ND", "NDGENIESm17", "Summed small GENIE Systs");
39  const Nus17SystFromHist kNus17NDCalibRelSyst (kNusAna17SystFileND, "ND", "NDCalRel17", "Relative systematics");
40  const Nus17SystFromHist kNus17NDNueCCSyst (kNusAna17SystFileND, "ND", "NDNueCC17", "#nu_{e} CC Background");
41  const Nus17SystFromHist kNus17NDNumuCCSyst (kNusAna17SystFileND, "ND", "NDNumuCC17", "#nu_{#mu} CC Background");
42  const Nus17SystFromHist kNus17NDNCSyst (kNusAna17SystFileND, "ND", "NDNC17", "NC Background");
43  const Nus17SystFromHist kNus17NDCherenkov_mSyst (kNusAna17SystFileND, "ND", "NDCherenkov_m17", "Light level method");
44  const Nus17SystFromHist kNus17NDLightLSyst (kNusAna17SystFileND, "ND", "NDCalLight17", "Lightlevel syst");
45 */
46 
47  // const Nus17SystFromHist kNus17FDBeamSysts (kNusAna17SystFileND, "FD", "NDBeam17", "All Beam");
48  // const Nus17SystFromHist kNus17FDCherenkovSyst (kNusAna17SystFileND, "FD", "NDCherenkov17", "Cherenkov");
49  // const Nus17SystFromHist kNus17FDCalibFlatSyst (kNusAna17SystFileND, "FD", "NDCalFlat17", "Flat Miscalibration");
50  // const Nus17SystFromHist kNus17FDCalibShapeSyst (kNusAna17SystFileND, "FD", "NDCalShape17", "Calibration shape");
51  // const Nus17SystFromHist kNus17FDGENIESmallSysts (kNusAna17SystFileND, "FD", "NDGENIESm17", "Summed small GENIE Systs");
52  // const Nus17SystFromHist kNus17FDCalibRelSyst (kNusAna17SystFileND, "FD", "NDCalRel17", "Relative systematics");
53  // const Nus17SystFromHist kNus17FDCherenkov_mSyst (kNusAna17SystFileND, "FD", "NDCherenkov_m17", "Light level method");
54  // const Nus17SystFromHist kNus17FDLightLSyst (kNusAna17SystFileND, "FD", "NDCalLight17", "Lightlevel syst");
55 
56 
57 
58 
59  // Systematics applied as flat shifts
60 // const NusFlatSyst kNusMCStatsSyst("mcstat", "MC Stats", 2.0, 4.8);
61 // const NusFlatSyst kNusNDContSyst ("ndcont", "ND Containment", 1.0, 0.6);
62  const Nus17FlatSyst kNus17NormSyst ("normNus17","Normalization", 4.9, 4.9);
63  const Nus17FlatSyst kNus17OscParamSysts ("oscparm", "Oscillation_Parameter",0.7,10.7);
64  const Nus17FlatSyst kNusFDPPFXSyst ("ppfx_fd","Fluxweight_fd",0.96,2.8);
65 
66  //----------------------------------------------------------------------
68  const std::string &det,
69  const std::string &shortname,
70  const std::string &latexname)
71  : ISyst(shortname, latexname),
72  fFileName(fname),
73  fDet(det)
74  {
75  }
76 
77  //----------------------------------------------------------------------
79  {
80  }
81 
82  //----------------------------------------------------------------------
84  {
85  // Someone already called us
86  if(!fHists.empty()) return;
87 
88  TFile fin(fFileName.c_str(), "read");
89  if(fin.IsZombie()) {
90  std::cout << "Warning: couldn't open " << fFileName
91  << ". Crashing" << std::endl;
92  abort();
93  }
94 
95  std::vector<std::string> channels = {"NC", "BG"};
96 
97  const std::vector<int> sigmas = {-1, 0, +1};
98  const std::vector<std::string> sigstr = {"-1", "0", "+1"};
99 
100  for(int i_chan = 0; i_chan < (int)channels.size(); ++i_chan) {
101  std::vector< std::pair<int,TH1D*> > curHists;
102 
103  for(int i_sig = 0; i_sig < (int)sigmas.size(); ++i_sig) {
104  // Find histogram for given NusChannel at given sigma
105  std::string hName = TString::Format("h%s_%s_%s_%s",
106  channels[i_chan].c_str(),
107  fDet.c_str(),
108  ShortName().c_str(),
109  sigstr[i_sig].c_str()).Data();
110  TH1D* h = (TH1D*)fin.Get(hName.c_str());
111  if(!h) {
112  // Surely it's a bad thing if we can't find nom, +-1 sigma hists
113  std::cout << "Error: can't find necessary " << hName
114  << " histogram in file " << fFileName
115  << ". Crashing" << std::endl;
116  abort();
117  }
118 
119  h->SetDirectory(0);
120  curHists.emplace_back(sigmas[i_sig], h);
121  }
122 
123  fHists.push_back(curHists);
124  }
125  }
126 
127  //----------------------------------------------------------------------
129  {
130  if(sr->mc.nnu == 0) { return NusChannel::kNC; }
131  if(!sr->mc.nu[0].iscc) { return NusChannel::kNC; }
132  if( sr->mc.nu[0].iscc) { return NusChannel::kBG; }
133  assert(0 && "Unknown Oscillation Channel");
134  }
135 
136  //----------------------------------------------------------------------
138  double sigma,
139  double calE) const
140  {
141  LoadHists();
142 
143  const int bin = fHists[chan][0].second->FindBin(calE);
144 
145  int LowIdx = 0;
146  if (sigma < fHists[chan].front().first)
147  LowIdx = 0;
148  else if (sigma >= fHists[chan].back().first)
149  LowIdx = fHists[chan].size()-2;
150  else{
151  for (int i = 0; i < (int)fHists[chan].size()-1; i++){
152  if (sigma >= fHists[chan][i].first){
153  LowIdx = i;
154  break;
155  }
156  }
157  }
158 
159  // Why would we have templates differing by more than 1 sigma?
160  // fracpart below assumes this
161  assert(fHists[chan][LowIdx+1].first - fHists[chan][LowIdx].first == 1);
162 
163  const double fracpart = sigma - fHists[chan][LowIdx].first;
164  const double ret = fracpart* fHists[chan][LowIdx+1].second->GetBinContent(bin) +
165  (1-fracpart)*fHists[chan][LowIdx] .second->GetBinContent(bin);
166 
167  return std::max(0., ret); // Keep the LL from blowing up
168  }
169 
170  //----------------------------------------------------------------------
172  caf::SRProxy* sr, double& weight) const
173  {
174  if(sr->hdr.det == caf::kFARDET &&
175  fDet.compare("FD") != 0 &&
176  fDet.compare("EX") != 0) { return; }
177  if(sr->hdr.det == caf::kNEARDET &&
178  fDet.compare("ND") != 0) { return; }
179 
180  // Find the energy estimate for this slice
181  double calE = kNus17Energy(sr);
182 
183 
184  // Find the signal or background oscillation channel
186 
187  weight *= WeightFor(chan, sigma, calE);
188  }
189 
190  //----------------------------------------------------------------------
192  const std::string &latexname,
193  const double& ncWei,
194  const double& bgWei)
195  : ISyst(shortname, latexname),
196  fNCWei(ncWei/100.),
197  fBGWei(bgWei/100.)
198  {
199  }
200 
201  //----------------------------------------------------------------------
203  {
204  }
205 
206  //----------------------------------------------------------------------
208  caf::SRProxy* sr, double& weight) const
209  {
210  // Takes in output from extrap - only shift FD spectra
211  if(sr->hdr.det != caf::kFARDET) return;
212  if(sr->mc.nnu == 0) return;
213  if(sigma == 0) return;
214 
215  weight *= ((sr->mc.nu[0].iscc) ?
216  1. + fBGWei*sigma :
217  1. + fNCWei*sigma);
218 
219  return;
220  }
221 
222  //----------------------------------------------------------------------
223  /// Get a vector of all the nus group systs
224 /* std::vector<const ISyst*> getAllNusFDSysts()
225  {
226  std::vector<const ISyst*> systs;
227 
228  systs.push_back(&kNus17OscParamSysts);
229  systs.push_back(&kNusFDPPFXSyst);
230  systs.push_back(&kNus17BeamSysts);
231  systs.push_back(&kNus17CherenkovSyst);
232  systs.push_back(&kNus17CalibFlatSyst);
233  systs.push_back(&kNus17CalibRelNDSyst);
234  systs.push_back(&kNus17CalibRelFDSyst);
235  systs.push_back(&kNus17CalibShapeSyst);
236  systs.push_back(&kNus17Cherenkov_mSyst);
237  systs.push_back(&kNus17GENIESmallSysts);
238  systs.push_back(&kNus17LightLSyst);
239  systs.push_back(&kNus17NumuCCSyst);
240  systs.push_back(&kNus17NCSyst);
241  systs.push_back(&kNus17AccSyst);
242  return systs;
243  }*/
244 //..............................................................................................
245 /* std::vector<const ISyst*> getAllNusNDSysts()
246  {
247  std::vector<const ISyst*> systs;
248 
249  systs.push_back(&kNus17OscParamSysts);
250  systs.push_back(&kNus17NDBeamSysts);
251  systs.push_back(&kNus17NDCherenkovSyst);
252  systs.push_back(&kNus17NDCalibFlatSyst);
253  systs.push_back(&kNus17NDCalibRelSyst);
254  systs.push_back(&kNus17NDCalibShapeSyst);
255  systs.push_back(&kNus17NDCherenkov_mSyst);
256  systs.push_back(&kNus17NDGENIESmallSysts);
257  systs.push_back(&kNus17NDLightLSyst);
258  systs.push_back(&kNus17NDNueCCSyst);
259  systs.push_back(&kNus17NDNumuCCSyst);
260  systs.push_back(&kNus17NDNCSyst);
261  systs.push_back(&kNus17NormSyst);
262  return systs;
263  }*/
264  //................................................................................................................
265 // std::vector<const ISyst*> getAllNusFDRSysts()
266 // {
267 // std::vector<const ISyst*> systs;
268 
269 // systs.push_back(&kNus17OscParamSysts);
270 // systs.push_back(&kNus17FDBeamSysts);
271 // systs.push_back(&kNus17FDCherenkovSyst);
272 // systs.push_back(&kNus17FDCalibFlatSyst);
273 // systs.push_back(&kNus17FDCalibRelSyst);
274 // systs.push_back(&kNus17FDCalibShapeSyst);
275 // systs.push_back(&kNus17FDCherenkov_mSyst);
276 // systs.push_back(&kNus17FDGENIESmallSysts);
277 // systs.push_back(&kNus17FDLightLSyst);
278 // systs.push_back(&kNus17NormSyst);
279 // return systs;
280 // }
281 
282 
283 }
Near Detector underground.
Definition: SREnums.h:10
T max(const caf::Proxy< T > &a, T b)
TString fin
Definition: Style.C:24
Far Detector at Ash River.
Definition: SREnums.h:11
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
virtual const std::string & ShortName() const final
The name printed out to the screen.
Definition: ISyst.h:27
const Var weight
caf::Proxy< caf::SRHeader > hdr
Definition: SRProxy.h:2004
Proxy for caf::StandardRecord.
Definition: SRProxy.h:1993
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:570
Nus17FlatSyst(const std::string &shortname, const std::string &latexname, const double &ncWei, const double &bgWei)
Definition: Nus17Systs.cxx:191
string fFileName
Definition: nue_pid_effs.C:43
const Nus17FlatSyst kNus17OscParamSysts("oscparm","Oscillation_Parameter", 0.7, 10.7)
caf::Proxy< short int > nnu
Definition: SRProxy.h:569
const Nus17FlatSyst kNus17NormSyst("normNus17","Normalization", 4.9, 4.9)
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
NusChannel GetNusChannel(caf::SRProxy *sr) const
Definition: Nus17Systs.cxx:128
const Var kNus17Energy([](const caf::SRProxy *sr){double cale=sr->slc.calE;double recoE=0.;if(sr->hdr.det==caf::kFARDET) recoE=FDscaleCalE17 *cale;if(sr->hdr.det==caf::kNEARDET) recoE=NDscaleCalE17 *cale;return recoE;})
Definition: NusVars.h:62
std::map< ToFCounter, std::vector< unsigned int > > channels
void LoadHists() const
Definition: Nus17Systs.cxx:83
double WeightFor(NusChannel chan, double sigma, double nueenergy) const
Definition: Nus17Systs.cxx:137
float bin[41]
Definition: plottest35.C:14
double sigma(TH1F *hist, double percentile)
OStream cout
Definition: OStream.cxx:6
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2005
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
NusChannel
Definition: NusSysts.h:13
std::vector< std::vector< std::pair< int, TH1D * > > > fHists
Definition: Nus17Systs.h:38
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Definition: Nus17Systs.cxx:171
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Definition: Nus17Systs.cxx:207
assert(nhit_max >=nhit_nbins)
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
static constexpr Double_t sr
Definition: Munits.h:164
const Nus17FlatSyst kNusFDPPFXSyst("ppfx_fd","Fluxweight_fd", 0.96, 2.8)
Nus17SystFromHist(const std::string &fname, const std::string &det, const std::string &shortname, const std::string &latexname)
Definition: Nus17Systs.cxx:67
caf::Proxy< caf::Det_t > det
Definition: SRProxy.h:227
std::string fFileName
Definition: Nus17Systs.h:29