NusSysts.cxx
Go to the documentation of this file.
2 
3 #include "CAFAna/Systs/Systs.h"
4 #include "CAFAna/Vars/Vars.h"
5 
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 NusSystFromHist kNusOscParamSysts (kNusAna01SystFAlt, "EX", "OscParam", "Oscillation Parameters");
19  const NusSystFromHist kNusBeamSysts (kNusAna01SystFile, "EX", "Beam", "All Beam");
20  const NusSystFromHist kNusBirksSyst (kNusAna01SystFile, "EX", "Birks", "Birks C");
21  const NusSystFromHist kNusCalibFlatSyst (kNusAna01SystFile, "EX", "CalFlat", "Flat Miscalibration");
22  const NusSystFromHist kNusCalibRelSyst (kNusAna01SystFile, "EX", "CalRel", "Relative Detector Miscalibration");
23  const NusSystFromHist kNusCalibSlopeXSyst(kNusAna01SystFile, "EX", "CalSlopeX", "Sloped Miscalibration, X");
24  const NusSystFromHist kNusCalibSlopeYSyst(kNusAna01SystFile, "EX", "CalSlopeY", "Sloped Miscalibration, Y");
25  const NusSystFromHist kNusGENIESmallSysts(kNusAna01SystFile, "EX", "GENIESm", "Summed small GENIE Systs");
26  const NusSystFromHist kNusNDRockSyst (kNusAna01SystFile, "EX", "NDNDRock", "ND Rock");
27  const NusSystFromHist kNusNueCCSyst (kNusAna01SystFile, "EX", "NueCC", "#nu_{e} CC Background");
28  const NusSystFromHist kNusNumuCCSyst (kNusAna01SystFile, "EX", "NumuCC", "#nu_{#mu} CC Background");
29 
30  const NusSystFromHist kNusNDBeamSysts (kNusAna01SystFile, "ND", "NDBeam", "All Beam");
31  const NusSystFromHist kNusNDBirksSyst (kNusAna01SystFile, "ND", "NDBirks", "Birks C");
32  const NusSystFromHist kNusNDCalibFlatSyst (kNusAna01SystFile, "ND", "NDCalFlat", "Flat Miscalibration");
33  const NusSystFromHist kNusNDCalibSlopeXSyst(kNusAna01SystFile, "ND", "NDCalSlopeX", "Sloped Miscalibration, X");
34  const NusSystFromHist kNusNDCalibSlopeYSyst(kNusAna01SystFile, "ND", "NDCalSlopeY", "Sloped Miscalibration, Y");
35  const NusSystFromHist kNusNDGENIESmallSysts(kNusAna01SystFile, "ND", "NDGENIESm", "Summed small GENIE Systs");
36  const NusSystFromHist kNusNDNDRockSyst (kNusAna01SystFile, "ND", "NDRock", "ND Rock");
37  const NusSystFromHist kNusNDNueCCSyst (kNusAna01SystFile, "ND", "NDNueCC", "#nu_{e} CC Background");
38  const NusSystFromHist kNusNDNumuCCSyst (kNusAna01SystFile, "ND", "NDNumuCC", "#nu_{#mu} CC Background");
39 
40  // Systematics applied as flat shifts
41  const NusFlatSyst kNusMCStatsSyst("mcstat", "MC Stats", 2.0, 4.8);
42  const NusFlatSyst kNusNDContSyst ("ndcont", "ND Containment", 1.0, 0.6);
43  const NusFlatSyst kNusNormSyst ("normNus","Normalization", 4.9, 4.9);
44 
45  //----------------------------------------------------------------------
47  const std::string &det,
48  const std::string &shortname,
49  const std::string &latexname)
50  : ISyst(shortname, latexname),
51  fFileName(fname),
52  fDet(det)
53  {
54  }
55 
56  //----------------------------------------------------------------------
58  {
59  }
60 
61  //----------------------------------------------------------------------
63  {
64  // Someone already called us
65  if(!fHists.empty()) return;
66 
67  TFile fin(fFileName.c_str(), "read");
68  if(fin.IsZombie()) {
69  std::cout << "Warning: couldn't open " << fFileName
70  << ". Crashing" << std::endl;
71  abort();
72  }
73 
74  std::vector<std::string> channels = {"NC", "BG"};
75 
76  const std::vector<int> sigmas = {-1, 0, +1};
77  const std::vector<std::string> sigstr = {"-1", "0", "+1"};
78 
79  for(int i_chan = 0; i_chan < (int)channels.size(); ++i_chan) {
80  std::vector< std::pair<int,TH1D*> > curHists;
81 
82  for(int i_sig = 0; i_sig < (int)sigmas.size(); ++i_sig) {
83  // Find histogram for given NusChannel at given sigma
84  std::string hName = TString::Format("h%s_%s_%s_%s",
85  channels[i_chan].c_str(),
86  fDet.c_str(),
87  ShortName().c_str(),
88  sigstr[i_sig].c_str()).Data();
89  TH1D* h = (TH1D*)fin.Get(hName.c_str());
90  if(!h) {
91  // Surely it's a bad thing if we can't find nom, +-1 sigma hists
92  std::cout << "Error: can't find necessary " << hName
93  << " histogram in file " << fFileName
94  << ". Crashing" << std::endl;
95  abort();
96  }
97 
98  h->SetDirectory(0);
99  curHists.emplace_back(sigmas[i_sig], h);
100  }
101 
102  fHists.push_back(curHists);
103  }
104  }
105 
106  //----------------------------------------------------------------------
108  {
109  if(sr->mc.nnu == 0) { return NusChannel::kNC; }
110  if(!sr->mc.nu[0].iscc) { return NusChannel::kNC; }
111  if( sr->mc.nu[0].iscc) { return NusChannel::kBG; }
112  assert(0 && "Unknown Oscillation Channel");
113  }
114 
115  //----------------------------------------------------------------------
117  double sigma,
118  double calE) const
119  {
120  LoadHists();
121 
122  const int bin = fHists[chan][0].second->FindBin(calE);
123 
124  int LowIdx = 0;
125  if (sigma < fHists[chan].front().first)
126  LowIdx = 0;
127  else if (sigma >= fHists[chan].back().first)
128  LowIdx = fHists[chan].size()-2;
129  else{
130  for (int i = 0; i < (int)fHists[chan].size()-1; i++){
131  if (sigma >= fHists[chan][i].first){
132  LowIdx = i;
133  break;
134  }
135  }
136  }
137 
138  // Why would we have templates differing by more than 1 sigma?
139  // fracpart below assumes this
140  assert(fHists[chan][LowIdx+1].first - fHists[chan][LowIdx].first == 1);
141 
142  const double fracpart = sigma - fHists[chan][LowIdx].first;
143  const double ret = fracpart* fHists[chan][LowIdx+1].second->GetBinContent(bin) +
144  (1-fracpart)*fHists[chan][LowIdx] .second->GetBinContent(bin);
145 
146  return std::max(0., ret); // Keep the LL from blowing up
147  }
148 
149  //----------------------------------------------------------------------
151  caf::SRProxy* sr, double& weight) const
152  {
153  if(sr->hdr.det == caf::kFARDET &&
154  fDet.compare("FD") != 0 &&
155  fDet.compare("EX") != 0) { return; }
156  if(sr->hdr.det == caf::kNEARDET &&
157  fDet.compare("ND") != 0) { return; }
158 
159  // Find the energy estimate for this slice
160  double calE = kCaloE(sr);
161 
162  // Find the signal or background oscillation channel
164 
165  weight *= WeightFor(chan, sigma, calE);
166  }
167 
168  //----------------------------------------------------------------------
170  const std::string &latexname,
171  const double& ncWei,
172  const double& bgWei)
173  : ISyst(shortname, latexname),
174  fNCWei(ncWei/100.),
175  fBGWei(bgWei/100.)
176  {
177  }
178 
179  //----------------------------------------------------------------------
181  {
182  }
183 
184  //----------------------------------------------------------------------
186  caf::SRProxy* sr, double& weight) const
187  {
188  // Takes in output from extrap - only shift FD spectra
189  if(sr->hdr.det != caf::kFARDET) return;
190  if(sr->mc.nnu == 0) return;
191  if(sigma == 0) return;
192 
193  weight *= ((sr->mc.nu[0].iscc) ?
194  1. + fBGWei*sigma :
195  1. + fNCWei*sigma);
196 
197  return;
198  }
199 
200  //----------------------------------------------------------------------
201  /// Get a vector of all the nus group systs
202  std::vector<const ISyst*> getAllNusSysts()
203  {
204  std::vector<const ISyst*> systs;
205 
206  // systs.push_back(&kNusOscParamSysts);
207  systs.push_back(&kNusBeamSysts);
208  systs.push_back(&kNusBirksSyst);
209  systs.push_back(&kNusCalibFlatSyst);
210  systs.push_back(&kNusCalibRelSyst);
211  systs.push_back(&kNusCalibSlopeXSyst);
212  systs.push_back(&kNusCalibSlopeYSyst);
213  systs.push_back(&kNusGENIESmallSysts);
214  systs.push_back(&kNusMCStatsSyst);
215  systs.push_back(&kNusNDContSyst);
216  systs.push_back(&kNusNDRockSyst);
217  systs.push_back(&kNusNormSyst);
218  systs.push_back(&kNusNueCCSyst);
219  systs.push_back(&kNusNumuCCSyst);
220 
221  return systs;
222  }
223 }
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
double fNCWei
Definition: NusSysts.h:56
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
std::string fFileName
Definition: NusSysts.h:31
NusFlatSyst(const std::string &shortname, const std::string &latexname, const double &ncWei, const double &bgWei)
Definition: NusSysts.cxx:169
const NusSystFromHist kNusNDBirksSyst(kNusAna01SystFile,"ND","NDBirks","Birks C")
Definition: NusSysts.h:81
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
const NusFlatSyst kNusNDContSyst("ndcont","ND Containment", 1.0, 0.6)
Definition: NusSysts.h:92
Proxy for caf::StandardRecord.
Definition: SRProxy.h:1993
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:570
string fFileName
Definition: nue_pid_effs.C:43
double fBGWei
Definition: NusSysts.h:57
caf::Proxy< short int > nnu
Definition: SRProxy.h:569
const NusSystFromHist kNusNDCalibSlopeYSyst(kNusAna01SystFile,"ND","NDCalSlopeY","Sloped Miscalibration, Y")
Definition: NusSysts.h:84
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
const NusSystFromHist kNusCalibFlatSyst(kNusAna01SystFile,"EX","CalFlat","Flat Miscalibration")
Definition: NusSysts.h:71
const NusSystFromHist kNusNDNDRockSyst(kNusAna01SystFile,"ND","NDRock","ND Rock")
Definition: NusSysts.h:86
const NusSystFromHist kNusBeamSysts(kNusAna01SystFile,"EX","Beam","All Beam")
Definition: NusSysts.h:69
const NusSystFromHist kNusNDCalibSlopeXSyst(kNusAna01SystFile,"ND","NDCalSlopeX","Sloped Miscalibration, X")
Definition: NusSysts.h:83
std::vector< std::vector< std::pair< int, TH1D * > > > fHists
Definition: NusSysts.h:40
std::vector< const ISyst * > getAllNusSysts()
Get a vector of all the nus group systs.
Definition: NusSysts.cxx:202
std::map< ToFCounter, std::vector< unsigned int > > channels
const Var kCaloE
Summed calorimetric energy of all hits in slice, uncorrected.
Definition: Vars.cxx:52
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Definition: NusSysts.cxx:185
void LoadHists() const
Definition: NusSysts.cxx:62
const NusSystFromHist kNusBirksSyst(kNusAna01SystFile,"EX","Birks","Birks C")
Definition: NusSysts.h:70
const NusFlatSyst kNusMCStatsSyst("mcstat","MC Stats", 2.0, 4.8)
Definition: NusSysts.h:91
float bin[41]
Definition: plottest35.C:14
const NusSystFromHist kNusGENIESmallSysts(kNusAna01SystFile,"EX","GENIESm","Summed small GENIE Systs")
Definition: NusSysts.h:75
const std::string kNusAna01SystFile
Definition: NusSysts.h:62
double sigma(TH1F *hist, double percentile)
const NusSystFromHist kNusNueCCSyst(kNusAna01SystFile,"EX","NueCC","#nu_{e} CC Background")
Definition: NusSysts.h:77
OStream cout
Definition: OStream.cxx:6
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2005
const NusSystFromHist kNusNDGENIESmallSysts(kNusAna01SystFile,"ND","NDGENIESm","Summed small GENIE Systs")
Definition: NusSysts.h:85
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Definition: NusSysts.cxx:150
NusSystFromHist(const std::string &fname, const std::string &det, const std::string &shortname, const std::string &latexname)
Definition: NusSysts.cxx:46
NusChannel
Definition: NusSysts.h:13
const NusSystFromHist kNusCalibRelSyst(kNusAna01SystFile,"EX","CalRel","Relative Detector Miscalibration")
Definition: NusSysts.h:72
const NusFlatSyst kNusNormSyst("normNus","Normalization", 4.9, 4.9)
Definition: NusSysts.h:93
double WeightFor(NusChannel chan, double sigma, double nueenergy) const
Definition: NusSysts.cxx:116
assert(nhit_max >=nhit_nbins)
const NusSystFromHist kNusNDRockSyst(kNusAna01SystFile,"EX","NDNDRock","ND Rock")
Definition: NusSysts.h:76
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
std::string fDet
Definition: NusSysts.h:32
const NusSystFromHist kNusNDNumuCCSyst(kNusAna01SystFile,"ND","NDNumuCC","#nu_{#mu} CC Background")
Definition: NusSysts.h:88
const NusSystFromHist kNusCalibSlopeXSyst(kNusAna01SystFile,"EX","CalSlopeX","Sloped Miscalibration, X")
Definition: NusSysts.h:73
const NusSystFromHist kNusCalibSlopeYSyst(kNusAna01SystFile,"EX","CalSlopeY","Sloped Miscalibration, Y")
Definition: NusSysts.h:74
const NusSystFromHist kNusNDBeamSysts(kNusAna01SystFile,"ND","NDBeam","All Beam")
Definition: NusSysts.h:80
static constexpr Double_t sr
Definition: Munits.h:164
NusChannel GetNusChannel(caf::SRProxy *sr) const
Definition: NusSysts.cxx:107
const NusSystFromHist kNusNDCalibFlatSyst(kNusAna01SystFile,"ND","NDCalFlat","Flat Miscalibration")
Definition: NusSysts.h:82
caf::Proxy< caf::Det_t > det
Definition: SRProxy.h:227
const NusSystFromHist kNusNumuCCSyst(kNusAna01SystFile,"EX","NumuCC","#nu_{#mu} CC Background")
Definition: NusSysts.h:78
const NusSystFromHist kNusNDNueCCSyst(kNusAna01SystFile,"ND","NDNueCC","#nu_{e} CC Background")
Definition: NusSysts.h:87