ShifterAndWeighter.cxx
Go to the documentation of this file.
1 //
2 // ShifterAndWeighter.cxx
3 // Implementation of singleton to calculate weights and keep track of shifts
4 // for FNEX fitting
5 //
6 
7 #include "TH2.h"
8 #include "TFile.h"
9 #include "TKey.h"
10 #include "TCanvas.h"
11 #include "TF1.h"
12 #include <algorithm>
13 
15 
17 
18 #include "Utilities/func/MathUtil.h"
19 
22 #include "FNEX/core/VarVals.h"
23 
27 
28 namespace fnex {
29 
30  //=======================================================================//
31  //------- GENERIC/GENERAL PURPOSE CODE (INITIALISATION ETC.) ------------//
32  //=======================================================================//
33 
34 
36 
37  //----------------------------------------------------------------------------
39  {
40  if(!gShifterAndWeighter) gShifterAndWeighter = new ShifterAndWeighter();
41 
42  return gShifterAndWeighter;
43  }
44 
45  //----------------------------------------------------------------------------
47  {
48  return;
49  }
50 
51  //----------------------------------------------------------------------------
53  {
54  return;
55  }
56 
57  //----------------------------------------------------------------------------
59  std::string CalibFilename)
60  {
61  // this method sets which systematics are to be used in calculating weights
62  // by creating a set of std::functions corresponding to the weighting functions
63  // that are to be turned on
64  // the code is a bit tedious, but it gets the job done
65 
66  // go ahead and set the current values
67  this->SetCurrentVals(curPoint);
68 
69  // loop over the systematic parameters from the input point and evaluate the
70  // name associated with each one, then add the appropriate function to the
71  // list.
72  fWeightersToUseND.clear();
73  fWeightersToUseFD.clear();
74  fShiftersToUseND .clear();
75  fShiftersToUseFD .clear();
76 
77  // Reset the Nue signal extrapolation weights
78  fSigExtrapWeights.clear();
79 
80  std::string systName;
82  WeighterInfo weighter;
83  for(auto const& itr : fCurrentPoint.SystematicPulls()){
84  systName = itr.first.first;
85  det = itr.first.second;
86 
87  if(systName.find("Wgt") != std::string::npos){
88  weighter.wgtName = systName;
89 
90  // check all the systematic parameter names and add the necessary functions
91  // for weighting if there is a function for that name.
92 
93  if(systName.compare("TauScaleSystWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::TauScaleSystWeight, this);
94  else if(systName.compare("NCNormWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::NCNormWeight, this);
95  else if(systName.compare("MECNormWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::MECNormWeight, this);
96  else if(systName.compare("CosmicMuNormWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::CosmicMuNormWeight, this);
97  else if(systName.compare("NueAbsNormWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::AbsNormWeight, this, "NueAbsNormWgt" , fnex::kNuESelection );
98  else if(systName.compare("AbsNormWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::AbsNormWeight, this, "AbsNormWgt" , fnex::kNuMuSelection);
99  else if(systName.compare("FHCNormWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::HornNormWeight, this, "FHCNormWgt" , fnex::kFHC );
100  else if(systName.compare("RHCNormWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::HornNormWeight, this, "RHCNormWgt" , fnex::kRHC );
101  else if(systName.compare("NueRelNormWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::RelNormWeight, this, "NueRelNormWgt" , fnex::kNuESelection );
102  else if(systName.compare("RelNormWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::RelNormWeight, this, "RelNormWgt" , fnex::kNuMuSelection);
103  else if(systName.compare("RockNormOnFidWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::RockNormOnFidWeight, this);
104  else if(systName.compare("CalibNearSystWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::CalibNormWeight, this, "CalibNearSystWgt", novadaq::cnv::kNEARDET);
105  else if(systName.compare("CalibFarSystWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::CalibNormWeight, this, "CalibFarSystWgt" , novadaq::cnv::kFARDET );
106  else if(systName.compare("BirksNearSystWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::BirksNormWeight, this, "BirksNearSystWgt", novadaq::cnv::kNEARDET);
107  else if(systName.compare("BirksFarSystWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::BirksNormWeight, this, "BirksFarSystWgt" , novadaq::cnv::kFARDET );
108  else if(systName.compare("NoiseFarSystWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::NoiseFarNormWeight, this);
109  else if(systName.compare("XsecCVWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::XSecCVPPFXWeight, this);
110  else if(systName.compare("RPACCQEWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::RPACCQEWeight, this);
111  else if(systName.compare("RPARESWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::RPARESWeight, this);
112  else if(systName.compare("RPACCQEshapeEnhWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::RPACCQEshapeEnhWeight, this);
113  else if(systName.compare("RPACCQEshapeSuppWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::RPACCQEshapeSuppWeight, this);
114  else if(systName.compare("DISvnCC1piWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::DISvnCC1piWeight, this);
115  else if(systName.compare("MECShapeNuWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::MECShapeWeight, this, false);
116  else if(systName.compare("MECShapeAntiNuWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::MECShapeWeight, this, true );
117  else if(systName.compare("MECInitNPFracNuWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::MECInitStateNPFracWeight, this, false);
118  else if(systName.compare("MECInitNPFracAntiNuWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::MECInitStateNPFracWeight, this, true );
119  else if(systName.compare("MECEnuShapeNuWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::MECEnuShapeWeight, this, false);
120  else if(systName.compare("MECEnuShapeAntiNuWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::MECEnuShapeWeight, this, true );
121  else if(systName.compare("MaCCQEReducedWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::MaCCQEReducedWeight, this);
122  else if(systName.compare("CosmicOverlayNormWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::CosmicOverlayNormWeight, this);
123  else if(systName.compare("NueBirksCSystCVNWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::NueBirksCWeight, this);
124  else if(systName.compare("NueShiftXSystCVNWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::NueShiftXWeight, this);
125  else if(systName.compare("NueShiftYSystCVNWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::NueShiftYWeight, this);
126  else if(systName.compare("NueXYAbsSystCVNWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::NueShiftXYAbsWeight, this);
127  else if(systName.compare("NueXYRelSystCVNWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::NueShiftXYRelWeight, this);
128  else if(systName.compare("RadCorrNueWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::RadCorrNueWeight, this);
129  else if(systName.compare("RadCorrNueBarWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::RadCorrNueBarWeight, this);
130  else if(systName.compare("SecondClassCurrWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::SecondClassCurrWeight, this);
131  else if(systName.compare("NueExtrapBkg2017Wgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::NueExtrapBkg2017Weight, this);
132  else if(systName.compare("NueExtrapSig2017Wgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::NueExtrapSig2017Weight, this);
133  else if(systName.compare("NearBeamSystWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::BeamSystWeight, this, "NearBeamSystWgt" , novadaq::cnv::kNEARDET);
134  else if(systName.compare("FarBeamSystWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::BeamSystWeight, this, "FarBeamSystWgt" , novadaq::cnv::kFARDET );
135  else if(systName.compare("NearBeamSystWgt2017" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::BeamSystWeight, this, "NearBeamSystWgt2017" , novadaq::cnv::kNEARDET);
136  else if(systName.compare("FarBeamSystWgt2017" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::BeamSystWeight, this, "FarBeamSystWgt2017" , novadaq::cnv::kFARDET );
137  else if(systName.compare("PPFXSmoothCVWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::PPFXSmoothCVWeight, this);
138  else if(systName.compare("FluxPCA0Wgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::FluxPCAWeight, this, 0);
139  else if(systName.compare("FluxPCA1Wgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::FluxPCAWeight, this, 1);
140  else if(systName.compare("FluxPCA2Wgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::FluxPCAWeight, this, 2);
141  else if(systName.compare("FluxPCA3Wgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::FluxPCAWeight, this, 3);
142  else if(systName.compare("FluxPCA4Wgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::FluxPCAWeight, this, 4);
143  else if(systName.compare("GeniePCA0Wgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GeniePCASystWeight, this, 0);
144  else if(systName.compare("GeniePCA1Wgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GeniePCASystWeight, this, 1);
145  else if(systName.compare("GeniePCA2Wgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GeniePCASystWeight, this, 2);
146  else if(systName.compare("GeniePCA3Wgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GeniePCASystWeight, this, 3);
147  else if(systName.compare("GeniePCA4Wgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GeniePCASystWeight, this, 4);
148  else if(systName.compare("NearCalibShapeWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::CalibSystWeight, this, "NearCalibShapeWgt", novadaq::cnv::kNEARDET);
149  else if(systName.compare("FarCalibShapeWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::CalibSystWeight, this, "FarCalibShapeWgt", novadaq::cnv::kFARDET );
150  else if(systName.compare("NearLightLevelWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::CalibSystWeight, this, "NearLightLevelWgt", novadaq::cnv::kNEARDET);
151  else if(systName.compare("FarLightLevelWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::CalibSystWeight, this, "FarLightLevelWgt", novadaq::cnv::kFARDET );
152  else if(systName.compare("NearCherenkovWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::CalibSystWeight, this, "NearCherenkovWgt", novadaq::cnv::kNEARDET);
153  else if(systName.compare("FarCherenkovWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::CalibSystWeight, this, "FarCherenkovWgt", novadaq::cnv::kFARDET );
154  else if(systName.compare("LightLevelWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::CalibSystWeight, this, "LightLevelWgt", novadaq::cnv::kFARDET );
155  else if(systName.compare("AbsCalibrationWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::CalibSystWeight, this, "AbsCalibrationWgt", novadaq::cnv::kFARDET );
156  else if(systName.compare("RelCalibrationWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::CalibSystWeight, this, "RelCalibrationWgt", novadaq::cnv::kFARDET );
157  else if(systName.compare("CalibShapeWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::CalibSystWeight, this, "CalibShapeWgt", novadaq::cnv::kFARDET );
158  else if(systName.compare("CherenkovWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::CalibSystWeight, this, "CherenkovWgt", novadaq::cnv::kFARDET );
159  else if(systName.compare("MaCCQE_shapeWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::MaCCQEShapeWeight, this);
160  else if(systName.compare("NormCCQEWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::NormCCQEWeight, this);
161  else if(systName.compare("MaCCRESWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::MaCCRESWeight, this);
162  else if(systName.compare("MvCCRESWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::MvCCRESWeight, this);
163  else if(systName.compare("MaNCRESWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::MaNCRESWeight, this);
164  else if(systName.compare("MvNCRESWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::MvNCRESWeight, this);
165  else if(systName.compare("MaNCELWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::MaNCELWeight, this);
166  else if(systName.compare("MFP_NWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::MFP_NWeight, this);
167  else if(systName.compare("CCQEPauliSupViaKFWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::CCQEPauliSupViaKFWeight, this);
168  else if(systName.compare("OtherGENIEWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::OtherGENIEWeight, this);
169  else if(systName.compare("SumSmallXSecNumu2017Wgt") == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::SumSmallXSecNumu2017Weight, this);
170  else if(systName.compare("SumSmallXSecNue2017Wgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::SumSmallXSecNue2017Weight, this);
171  else if(systName.compare("NueAcceptSignalKin2018FHCWgt") == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::NueAcceptSignalKin2018Weight, this, true);
172  else if(systName.compare("NueAcceptSignalKin2018RHCWgt") == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::NueAcceptSignalKin2018Weight, this, false);
173  else if(systName.compare("NueAcceptBkg2018FHCWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::NueAcceptBkg2018Weight, this, true);
174  else if(systName.compare("NueAcceptBkg2018RHCWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::NueAcceptBkg2018Weight, this, false);
175  else if(systName.compare("FormZoneWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::FormZoneWeight, this);
176  else if(systName.compare("FrInel_piWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::FrInelpiWeight, this);
177  else if(systName.compare("FrCEx_NWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::FrCExNWeight, this);
178  else if(systName.compare("FrElas_NWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::FrElasNWeight, this);
179  else if(systName.compare("FrAbs_NWgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::FrAbsNWeight, this);
180  else if(systName.compare("COHCCScale2018Wgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::COHCCScale2018Weight, this);
181  else if(systName.compare("COHNCScale2018Wgt" ) == 0 ) weighter.wgtFn = std::bind(&ShifterAndWeighter::COHNCScale2018Weight, this);
182  else{
183  LOG_WARNING("ShifterAndWeighter")
184  << "Unable to determine weighting function for "
185  << systName;
186  continue;
187  }
188 
189  LOG_DEBUG("ShifterAndWeighter")
190  << "Adding weighter "
191  << weighter.wgtName
192  << " to detector "
193  << itr.first.second;
194 
195  if (det == novadaq::cnv::kNEARDET) fWeightersToUseND.push_back(weighter);
196  else if(det == novadaq::cnv::kFARDET ) fWeightersToUseFD.push_back(weighter);
197 
198  // Decide which histograms we need
199  if(systName.find("BeamSystWgt" ) != std::string::npos ) {
200  if(systName.find("2017") != std::string::npos) this->LoadBeamSystHists2017();
201  else this->LoadSABeamSystHists();
202  }
203  else if(systName.find("FluxPCA" ) != std::string::npos )
204  this->LoadFluxBeamSystHists();
205  else if(systName.find("GeniePCA") != std::string::npos )
206  this->LoadGeniePCASystHists();
207  else if(systName.find("NueBirks") != std::string::npos ||
208  systName.find("NueShift") != std::string::npos )
209  this->LoadNueSystHists();
210  else if(systName.find("CalibShape") != std::string::npos ||
211  systName.find("LightLevel") != std::string::npos ||
212  systName.find("Cherenkov" ) != std::string::npos ||
213  systName.find("Calibration" ) != std::string::npos)
214  this->LoadCalibrationSystHists(CalibFilename);
215 
216  } // the parameter was for a weight
217  else if(systName.find("Shifter") != std::string::npos){
218 
219  LOG_DEBUG("ShifterAndWeighter")
220  << "Adding shifter for "
221  << systName
222  << " to detector "
223  << det;
224 
225  if(systName.compare("RecoHadEShifter") == 0){
226  this->FillShifterInfo(systName,
228  det);
229  }
230  else if(systName.compare("RecoHadEShifter2017") == 0){
231  this->FillShifterInfo(systName,
233  det);
234  }
235  else if(systName.compare("RecoLepEShifter") == 0){
236  this->FillShifterInfo(systName,
238  det);
239  }
240  else if(systName.compare("RecoLepEShifter2017") == 0){
241  this->FillShifterInfo(systName,
243  det);
244  }
245  else if(systName.compare("RelRecoHadEShifter") == 0){
246  this->FillShifterInfo(systName,
248  det);
249  }
250  else if(systName.compare("RelRecoHadEShifter2017") == 0){
251  this->FillShifterInfo(systName,
253  det);
254  }
255  else if(systName.compare("RelRecoLepEShifter") == 0){
256  this->FillShifterInfo(systName,
258  det);
259  }
260  else if(systName.compare("RelRecoLepEShifter2017") == 0){
261  this->FillShifterInfo(systName,
263  det);
264  }
265  else if(systName.compare("RecoNuECalibFarShifter") == 0){
266  this->FillShifterInfo(systName,
268  det);
269  this->FillShifterInfo(systName,
271  det);
272  }
273  else if(systName.compare("RecoNuEBirksFarShifter") == 0){
274  this->FillShifterInfo(systName,
276  det);
277  this->FillShifterInfo(systName,
279  det);
280  }
281  else if(systName.compare("RecoNuECalibNearShifter") == 0){
282  this->FillShifterInfo(systName,
284  det);
285  this->FillShifterInfo(systName,
287  det);
288  }
289  else if(systName.compare("RecoNuEBirksNearShifter") == 0){
290  this->FillShifterInfo(systName,
292  det);
293  this->FillShifterInfo(systName,
295  det);
296  }
297  else if(systName.compare("RecoNuENoiseFarShifter") == 0){
298  this->FillShifterInfo(systName,
300  det);
301  this->FillShifterInfo(systName,
303  det);
304  }
305  else if(systName.compare("RecoNuEShifter") == 0){
306  this->FillShifterInfo(systName,
308  det);
309  this->FillShifterInfo(systName,
311  det);
312  }
313 
314  } // end if this is a shifter
315  } // end loop over systematic parameters
316 
317 // for(auto const& itr : fShiftersToUseND){
318 // LOG_VERBATIM("ShifterAndWeighter")
319 // << "ND shifter param: "
320 // << itr.first;
321 // for(auto const& si : itr.second)
322 // LOG_VERBATIM("ShifterAndWeighter")
323 // << " associated shifter "
324 // << si.fShifterName
325 // << " "
326 // << si.fDetector;
327 // }
328 // for(auto const& itr : fShiftersToUseFD){
329 // LOG_VERBATIM("ShifterAndWeighter")
330 // << "FD shifter param: "
331 // << itr.first;
332 // for(auto const& si : itr.second)
333 // LOG_VERBATIM("ShifterAndWeighter")
334 // << " associated shifter "
335 // << si.fShifterName
336 // << " "
337 // << si.fDetector;
338 // }
339 
340  LOG_DEBUG("ShifterAndWeighter")
341  << "Done initialising shifters and weighters";
342 
343  return;
344  }
345 
346  //----------------------------------------------------------------------------
348  unsigned char const& param,
349  novadaq::cnv::DetId const& det)
350  {
351 
352  if(det == novadaq::cnv::kNEARDET) fShiftersToUseND[param].emplace(det, name);
353  else if(det == novadaq::cnv::kFARDET ) fShiftersToUseFD[param].emplace(det, name);
354 
355  return;
356  }
357 
358 
359  //----------------------------------------------------------------------------
361  {
362  LOG_DEBUG("ShifterAndWeighter")
363  << "Previous shifter/weighter point was\n"
364  << fCurrentPoint;
365 
366  fCurrentPoint = curPoint;
367 
368  LOG_DEBUG("ShifterAndWeighter")
369  << "Now shifter/weighter point is\n"
370  << fCurrentPoint;
371 
372  // get the oscillation parameters from the current point, we haven't set
373  // the detector yet, but the only thing that changes is the
374  // baseline
376  LOG_DEBUG("ShifterAndWeighter")
377  << itr.first
378  << " "
379  << itr.second;
380 
381  if(itr.first.compare("L" ) == 0) fOscCalc.SetL (itr.second);
382  else if(itr.first.compare("Rho" ) == 0) fOscCalc.SetRho (itr.second);
383  else if(itr.first.compare("Dmsq21" ) == 0) fOscCalc.SetDmsq21 (itr.second);
384  else if(itr.first.compare("Dmsq32" ) == 0) fOscCalc.SetDmsq32 (itr.second);
385  else if(itr.first.compare("Th12" ) == 0) fOscCalc.SetTh12 (itr.second);
386  else if(itr.first.compare("Th13" ) == 0) fOscCalc.SetTh13 (itr.second);
387  else if(itr.first.compare("Th23" ) == 0) fOscCalc.SetTh23 (itr.second);
388  else if(itr.first.compare("dCP" ) == 0) fOscCalc.SetdCP (itr.second);
389 // else if(itr.first.compare("Eps_ee" ) == 0) fOscCalc.SetEps_ee (itr.second); //Eps = NSI epsilon parameters
390 // else if(itr.first.compare("Eps_emu" ) == 0) fOscCalc.SetEps_emu (itr.second);
391 // else if(itr.first.compare("Eps_etau" ) == 0) fOscCalc.SetEps_etau (itr.second);
392 // else if(itr.first.compare("Eps_mumu" ) == 0) fOscCalc.SetEps_mumu (itr.second);
393 // else if(itr.first.compare("Eps_mutau" ) == 0) fOscCalc.SetEps_mutau (itr.second);
394 // else if(itr.first.compare("Eps_tautau" ) == 0) fOscCalc.SetEps_tautau (itr.second);
395 // else if(itr.first.compare("Delta_emu" ) == 0) fOscCalc.SetDelta_emu (itr.second);
396 // else if(itr.first.compare("Delta_etau" ) == 0) fOscCalc.SetDelta_etau (itr.second);
397 // else if(itr.first.compare("Delta_mutau" ) == 0) fOscCalc.SetDelta_mutau (itr.second);
398  }
399 
400  return;
401  }
402 
403  //----------------------------------------------------------------------------
405  fnex::MetaData const& md)
406  {
407  fCurrentEvent = curEvent.get();
408  fCurrentMD = md;
409 
410  LOG_DEBUG("SetCurrentEvent")
411  << " Metadata is "
412  << md.ToString();
413 
414  // check the validity of the mcVals before trying to get weights, if it is
415  // a nullptr, we have a problem.
416  if(fCurrentEvent->mcVals.get() == nullptr && md.isMC)
417  throw cet::exception("ShifterAndWeighter")
418  << "Event does not have a valid MCVarVals unique_ptr, bail";
419 
420  // the data values better always be valid
421  if(fCurrentEvent->dataVals.get() == nullptr)
422  throw cet::exception("ShifterAndWeighter")
423  << "Event does not have a valid DataVarVals unique_ptr, bail";
424 
425  // reset the oscillation baseline based on the detector
427 
428  return;
429  }
430 
431  //----------------------------------------------------------------------------
432  double ShifterAndWeighter::ShiftAmount(unsigned char const& param)
433  {
434  // loop over the shifters for the current detector and get the product of the
435  // shift amount
436  double totalShift = 1.;
437 
439  totalShift = this->TotalShift(fShiftersToUseND, param);
441  totalShift = this->TotalShift(fShiftersToUseFD, param);
442 
443  return totalShift;
444  }
445 
446  //----------------------------------------------------------------------------
447  double ShifterAndWeighter::TotalShift(std::map<unsigned char, std::set<ShifterInfo> > const& shiftersToUse,
448  unsigned char const& param)
449  {
450 
451  // find the set of shifters for this parameter
452  auto itr = shiftersToUse.find(param);
453  if(itr == shiftersToUse.end()){
454  LOG_DEBUG("ShifterAndWeighterTotalShift")
455  << "could not find shifter set for "
456  << fnex::KeyToVarName(param)
457  << " "
458  << param
459  << " return total shift of 1.";
460  return 1.;
461  }
462 
463  std::set<ShifterInfo> shifterInfos = itr->second;
464 
465  // The SystematicSigmas should trigger an exception if the shifter is not
466  // available, but then again, we are iterating over the shifters from the
467  // input point, so it better be. Same is true for getting the parameter
468  // value
469  double totalShift = 1.;
470 
471  for(auto const& itr : shifterInfos){
472  auto const& sigma = fnex::SystematicSigmas::Instance()->Sigma(itr.fShifterName);
473 
474  // we want to treat the relative energy shifters in an anti-correlated fashion
475  if(itr.fShifterName == "RelRecoLepEShifter2017" || itr.fShifterName == "RelRecoHadEShifter2017") {
477  // TODO: hardcoded muon catcher sigma because I can't think how best to do this
478  // properly right now, fix this later!!
479  double mcFrac = fCurrentEvent->dataVals->val_at(fnex::kLep_RecoE_MCFrac, fCurrentMD);
480  totalShift *= 1. - ((1 - mcFrac) * (sigma/2) + mcFrac * 0.0075) * fCurrentPoint.ParameterValue(std::make_pair(itr.fShifterName, fCurrentMD.detector));
481  }
483  totalShift *= 1. + (sigma/2) * fCurrentPoint.ParameterValue(std::make_pair(itr.fShifterName, fCurrentMD.detector));
484  }
486  itr.fShifterName == "RecoLepEShifter2017") {
487  double mcFrac = fCurrentEvent->dataVals->val_at(fnex::kLep_RecoE_MCFrac, fCurrentMD);
488  // TODO: hardcoded muon catcher sigma because I can't think how best to do this
489  // properly right now, fix this later!!
490  totalShift *= 1. + ((1 - mcFrac) * sigma + mcFrac * 0.0069) * fCurrentPoint.ParameterValue(std::make_pair(itr.fShifterName, fCurrentMD.detector));
491  }
492  else
493  totalShift *= 1. + sigma * fCurrentPoint.ParameterValue(std::make_pair(itr.fShifterName, fCurrentMD.detector));
494 
495  // LOG_VERBATIM("ShifterAndWeighterTotalShift")
496  // << " associated shifter "
497  // << itr.fShifterName
498  // << " "
499  // << itr.fDetector
500  // << " "
501  // << fCurrentMD.detector
502  // << " "
503  // << sigma
504  // << " "
505  // << fCurrentPoint.ParameterValue(std::make_pair(itr.fShifterName, fCurrentMD.detector))
506  // << " "
507  // << totalShift;
508 
509  }
510 
511  return totalShift;
512  }
513 
514  //----------------------------------------------------------------------------
516  fnex::MetaData const& md,
517  fnex::WeightType_t const wgtType)
518  {
519  double wgt = 1.;
520 
521  // if the mcVals is null or the metadata claims it isn't MC,
522  // return the fake weight for the event. If the event is data
523  // then the value is always 1, if it is fake data, then the
524  // value is the proper weighting
525  if(curEvent->mcVals.get() == nullptr || !md.isMC)
526  return curEvent->dataVals->val_at(fnex::kFake_Weight, md);
527 
528  // get ready to use the passed in event
529  this->SetCurrentEvent(curEvent, md);
530 
531  // calculate the weights for different systematic uncertainties, starting
532  // with the normalization ones and moving from there
533  if(wgtType != fnex::kOscWgtOnly) wgt *= this->TotalWeightFromWeighters();
534 
535  // if(md.selectionType == fnex::kNuESelectionPeripheral)
536  // LOG_VERBATIM("ShifterAndWeighter")
537  // << "weight from systematics is "
538  // << wgt;
539 
540  // if we only want weights from systematics, we are done
541  if(wgtType == fnex::kSystWgtOnly) return wgt;
542 
543  // calculate the oscillation parameters weight
544  if(wgtType != fnex::kSystWgtOnly) wgt *= this->OscillationWeight();
545 
546  // if(md.selectionType == fnex::kNuESelectionPeripheral)
547  // LOG_VERBATIM("ShifterAndWeighter")
548  // << "weight from oscillations is "
549  // << this->OscillationWeight();
550 
551  LOG_DEBUG("ShifterAndWeighter")
552  << "Oscillation weight is "
553  << wgt
554  << " energy: "
555  << fCurrentEvent->mcVals->val_at(fnex::kTrueE)
556  << " pdg: "
557  << fCurrentEvent->mcVals->val_at(fnex::kTruePDG)
558  << " orig pdg: "
560  << " detector: "
561  << fCurrentMD.detector;
562 
563  // apply any weights that may be due to extrapolation techniques
564  if(wgtType == fnex::kAllWgt) wgt *= this->ExtrapolationWeight();
565 
566  // if(md.selectionType == fnex::kNuESelectionPeripheral)
567  // LOG_VERBATIM("ShifterAndWeighter")
568  // << "total weight is "
569  // << wgt;
570 
571  return wgt;
572  }
573 
574  //----------------------------------------------------------------------------
576  {
577  double wgt = 1.;
578  double tmpWgt = 1.;
579 
581 
582  for(auto itr : weighters){
583  tmpWgt = itr.wgtFn();
584  wgt *= tmpWgt;
585  }
586 
587  return wgt;
588  }
589 
590  //----------------------------------------------------------------------------
592  {
593  double wgt = 1;
594 
595  if(fHornCurrentHists.empty()) {
597  cet::search_path sp("FW_SEARCH_PATH");
598  if ( !sp.find_file(std::string("FNEX/data/horncurrent/190kAHornCurrentHists.root"), fileName) ) {
599  throw cet::exception("HornCurrentHistError")
600  << "Cannot find horn current histograms.";
601  }
602 
603  TFile fin(fileName.c_str(), "READ");
604  if(fin.IsZombie()) {
605  throw cet::exception("HornCurrentHistError")
606  << "Couldn't open horn current file " << fileName;
607  }
608 
609  TString histNameND = "rat190ND";
610  TString histNameFD = "rat190FD";
611  TString histNameNDbar = "rat190NDan";
612  TString histNameFDbar = "rat190FDan";
613 
614  TH1D *hHornCurrentND = (TH1D*)fin.Get(histNameND)->Clone();
615  TH1D *hHornCurrentFD = (TH1D*)fin.Get(histNameFD)->Clone();
616  TH1D *hHornCurrentNDbar = (TH1D*)fin.Get(histNameNDbar)->Clone();
617  TH1D *hHornCurrentFDbar = (TH1D*)fin.Get(histNameFDbar)->Clone();
618 
619  hHornCurrentND->SetDirectory(0);
620  hHornCurrentFD->SetDirectory(0);
621  hHornCurrentNDbar->SetDirectory(0);
622  hHornCurrentFDbar->SetDirectory(0);
623 
624  fHornCurrentHists.push_back(hHornCurrentND );
625  fHornCurrentHists.push_back(hHornCurrentFD );
626  fHornCurrentHists.push_back(hHornCurrentNDbar);
627  fHornCurrentHists.push_back(hHornCurrentFDbar);
628  }
629 
630  double E = fCurrentEvent->mcVals->val_at(fnex::kTrueE);
631  int beamType = fCurrentMD.BeamType();
632  int pid = fCurrentEvent->mcVals->val_at(fnex::kTruePDG);
633 
635  if(E >= 20) return wgt;
636  if( (beamType == fnex::kFHC && pid < 0 && E >= 8) ||
637  (beamType == fnex::kRHC && pid > 0 && E >= 8) ) return wgt;
638  if(beamType != fnex::kFHC && beamType != fnex::kRHC) return wgt;
639 
641  if( (beamType == fnex::kFHC && pid > 0) || (beamType == fnex::kRHC && pid < 0) )
642  wgt = fHornCurrentHists.at(0)->Interpolate(E);
643  if( (beamType == fnex::kFHC && pid < 0) || (beamType == fnex::kRHC && pid > 0) )
644  wgt = fHornCurrentHists.at(2)->Interpolate(E);
645  }
647  if( (beamType == fnex::kFHC && pid > 0) || (beamType == fnex::kRHC && pid < 0) )
648  wgt = fHornCurrentHists.at(1)->Interpolate(E);
649  if( (beamType == fnex::kFHC && pid < 0) || (beamType == fnex::kRHC && pid > 0) )
650  wgt = fHornCurrentHists.at(3)->Interpolate(E);
651  }
652 
653  LOG_DEBUG("HornCurrentWeight") << "Horn current weight is: " << wgt;
654 
655  return wgt;
656  }
657 
658  //----------------------------------------------------------------------------
660  {
661  double wgt = 1;
662 
663  // need to see if it is the nue FD signal and only apply the signal weight
664  // to the correct metadata, ie a flux swap file and nue selected events
668  wgt *= this->NueSignalExtrapWeight();
669 
670  //need to see if it is the numu need to see if it is the
671  //to the correct metadata, ie a beam files and numu selected
672  //events in the FD
676  wgt *= this->NuMuSignalExtrapWeight();
677 
678  return wgt;
679  }
680 
681  //----------------------------------------------------------------------------
683  {
684  double E = fCurrentEvent->mcVals->val_at(fnex::kTrueE);
685  int pdg = fCurrentEvent->mcVals->val_at(fnex::kTruePDG);
686  int pdgOrig = fCurrentEvent->mcVals->val_at(fnex::kTruePDGOrig);
687 
688  // only attempt to weight neutrinos for oscillations
689  if(std::abs(pdg) != 12 &&
690  std::abs(pdg) != 14 &&
691  std::abs(pdg) != 16) return 1.;
692 
693 // Let's comment following part of the code because it is interferring in FNEX fit jobs
694 // if(fCurrentPoint.fUseBinnedOscillations) {
695 // if(!hCAFTrueEBins) {
696 //
697 // LOG_VERBATIM("OscillationWeight")
698 // << "We are using binned oscillations so generating true energy histogram for oscillation weights.";
699 //
700 // // emulate CAF true energy binning
701 // const int kNTrueEBins = 100;
702 // Double_t edges[kNTrueEBins+1];
703 //
704 // const double kEmin = 0.5; // not really events below 500 MeV
705 //
706 // const double kNEdges = kNTrueEBins - 1;
707 // const double kA = kNEdges * kEmin;
708 //
709 // edges[0] = 0;
710 // for(int i = 1; i <= kNEdges; ++i)
711 // edges[kNTrueEBins-i] = kA/i;
712 //
713 // edges[kNTrueEBins] = 120;
714 //
715 // hCAFTrueEBins = new TH1D("hBinning", "hBinning", kNTrueEBins, edges);
716 // }
717 //
718 // // grab the center of the CAF true E histogram to give us the true energy
719 // E = hCAFTrueEBins->GetBinCenter(hCAFTrueEBins->FindBin(E));
720 // }
721 
722  return fOscCalc.P(pdgOrig, pdg, E);
723 
724  }
725 
726  //============================================================================//
727  //-------------------- WEIGHT INTERPOLATION CALCULATIONS ---------------------//
728  //============================================================================//
729 
730  //----------------------------------------------------------------------------
732  std::map<float, float> const& wgtsBySigma)
733  {
734  double wgt = 1.;
735 
736  auto mapEnd = wgtsBySigma.end();
737  double badVal = std::numeric_limits<double>::min();
738  double minus2Val = (wgtsBySigma.find(-2.) != mapEnd) ? static_cast<double>(wgtsBySigma.find(-2.)->second) : badVal;
739  double minus1Val = (wgtsBySigma.find(-1.) != mapEnd) ? static_cast<double>(wgtsBySigma.find(-1.)->second) : badVal;
740  // double zeroVal = 1.;
741  double zeroVal = (wgtsBySigma.find( 0.) != mapEnd) ? static_cast<double>(wgtsBySigma.find( 0.)->second) : 1.;
742  double plus1Val = (wgtsBySigma.find( 1.) != mapEnd) ? static_cast<double>(wgtsBySigma.find( 1.)->second) : badVal;
743  double plus2Val = (wgtsBySigma.find( 2.) != mapEnd) ? static_cast<double>(wgtsBySigma.find( 2.)->second) : badVal;
744  double slope = 0.;
745  double intercept = 0.;
746 
747  //Debugging RPA weights
748  LOG_DEBUG("ShifterAndWeighter")
749  << "\tCalcLinearInterpWgt: "
750  << " -2:" << minus2Val
751  << ", -1:" << minus1Val
752  << ", 0:" << zeroVal
753  << ", 1:" << plus1Val
754  << ", 2:" << plus2Val;
755  if (wgtsBySigma.find(-1.) != mapEnd)
756  LOG_DEBUG("ShifterAndWeighter")
757  << "\t\t 0 val from " << wgtsBySigma.find( 0.)->second;
758 
759  // The following code assumes that at least the -1 sigma, +1 sigma
760  // values are in the map. If they are not, we have a problem
761  if(minus1Val == badVal ||
762  plus1Val == badVal){
763  LOG_WARNING("ShifterAndWeighterCalcLinearInterpWgt")
764  << "No values associated with either -1 sigma or +1 sigma "
765  << "so we cannot calculate a weight, return 1.";
766  return 1.;
767  }
768 
769  if(sigma <= -1){
770  // It is possible there is no -2 sigma value from the map,
771  // in that case just interpolate from 0 to -1 sigma
772  if(minus2Val != badVal){
773  slope = (minus1Val - minus2Val);
774  intercept = minus2Val - slope * (-2.);
775  }
776  else{
777  slope = (zeroVal - minus1Val);
778  intercept = minus1Val - slope * (-1.);
779  }
780  }
781  else if(sigma > -1. && sigma <= 0.){
782  slope = (zeroVal - minus1Val);
783  intercept = minus1Val - slope * (-1.);
784  }
785  else if(sigma > 0. && sigma <= 1.){
786  slope = (plus1Val - zeroVal);
787  intercept = plus1Val - slope * 1.;
788  }
789  else if(sigma > 1.){
790  // It is possible there is no +2 sigma value from the map,
791  // in that case just interpolate from 0 to +1 sigma
792  if(plus2Val != badVal){
793  slope = (plus2Val - plus1Val);
794  intercept = plus2Val - slope * 2.;
795  }
796  else{
797  slope = (plus1Val - zeroVal);
798  intercept = plus1Val - slope * 1.;
799  }
800  }
801 
802  wgt = sigma * slope + intercept;
803 
804  // Don't let weights be negative
805  if(wgt < 0) wgt = 0.;
806 
807  return wgt;
808  }
809 
810  //----------------------------------------------------------------------------
811  /// The interpolation done here is a cubic spline interpolation with "natural" boundary conditions.
812  /// (See https://en.wikipedia.org/wiki/Spline_interpolation).
813  /// This algorithm is due to J. Douglas Faires, Richard L. Burden in "Numerical Methods", Fourth Edition
814  /// (see section 3.5 and pseudocode therein).
815  /// This implementation assumes that you will provide a symmetrical number of integral sigmas
816  /// (except you don't have to include weight(0 sigma) = 1.0; that's assumed):
817  /// i.e., {-2, -1, 1, 2}.
818  /// Such an assumption eliminates the need to store the knots themselves, saving a lot of memory.
819  std::vector<std::array<float, 4>> ShifterAndWeighter::CalcInterpFitCoeffs(std::map<float, float> const& wgtsBySigma)
820  {
821  // -1 is to match notation in cited algorithm pseudocode text,
822  // making N = number of spline functions.
823 
824  const auto N = wgtsBySigma.size() - 1;
825 
826  //std::cout << " for sigma = " << this->fSigma << ":";
827  //if(this->fSigma != 0)
828  //{
829  // std::cout << " input weights: ";
830  // for (const auto & wgtPair : wgtsBySigma)
831  // std::cout << "(" << wgtPair.first << "," << wgtPair.second << ")\t";
832  // std::cout << std::endl;
833  //}
834 
835  // there are N=(# points - 1) sets of 4 coefficients.
836  // they are calculated as Legendre polynomials (expanded around the knots)
837  // and at the end we'll convert them to polynoimals expanded around x=0
838  // so that we don't have to save the knots themselves.
839  std::vector<std::array<float, 4>> legendreCoeffs(N);
840 
841  // the first (constant) term in the Legendre expansion
842  // is always the function value. we need to look both ways
843  // in this term so just fill it first.
844  auto it_prevKnot = wgtsBySigma.begin();
845  auto it_thisKnot = wgtsBySigma.begin();
846  auto it_nextKnot = wgtsBySigma.begin();
847  it_nextKnot++;
848 
849  // loop number 1 calculates the working values ...
850  std::vector<float> h(N), z(N); // working values
851  std::vector<float> mu(N);
852  z[0] = 0;
853  mu[0] = 0;
854  legendreCoeffs[0][0] = it_thisKnot->second; // i.e., a_0 = f(x_0)
855  legendreCoeffs.back()[2] = 0; // i.e., c_n = 0
856  for(unsigned int knotNum=0; it_nextKnot != wgtsBySigma.end(); ++it_thisKnot, ++it_nextKnot, ++knotNum){
857  //std::cout << " knot number " << knotNum << ":" << std::endl;
858  //std::cout << " determining intermediate values for knot #" << knotNum << std::endl;
859  if ( std::abs(it_nextKnot->first - it_thisKnot->first - 1) > 0.0001 )
860  throw cet::exception("InterpSystWgtFn::CalcInterpFitCoeffs()")
861  << "provided sigmas are required to be integers and separated by exactly 1 sigma";
862 
863  legendreCoeffs[knotNum+1][0] = it_nextKnot->second; // i.e., a_(i+1) = f(x_(i+1))
864  h[knotNum] = it_nextKnot->first - it_thisKnot->first; // i.e., h_i = x_(i+1) - x_i
865 
866  if (knotNum == 0) continue;
867 
868  //std::cout << " a_(i-1) = " << legendreCoeffs[knotNum-1][0] << std::endl;
869  //std::cout << " a_(i+1) = " << legendreCoeffs[knotNum+1][0] << std::endl;
870  //std::cout << " h_(i-1) = " << h[knotNum-1] << std::endl;
871 
872  float alpha = 3 * ( ( ( legendreCoeffs[knotNum+1][0] - legendreCoeffs[knotNum][0]) / h[knotNum] )
873  - ( ( legendreCoeffs[knotNum][0] - legendreCoeffs[knotNum-1][0]) / h[knotNum-1] ) );
874  //std::cout << " alpha = " << alpha << std::endl;
875 
876  // it_prevKnot is ok here since we don't get to this point in the loop the first time through.
877  // similarly, this 'l' is always for i >= 1
878  float l = 2 * (it_nextKnot->first - it_prevKnot->first) - h[knotNum-1] * mu[knotNum-1];
879  mu[knotNum] = h[knotNum] / l;
880  z[knotNum] = (alpha - h[knotNum-1]*z[knotNum-1]) / l;
881 
882  //std::cout << " l = " << l << std::endl;
883  //std::cout << " mu = " << mu[knotNum] << std::endl;
884  //std::cout << " z = " << z[knotNum] << std::endl;
885 
886  it_prevKnot = it_thisKnot;
887  }
888 
889  // ... loop #2 runs 'backwards' and fills the Legendre polynomial coefficients...
890  // (note weird loop condition because unsigned ints will underflow when pushed 'below' zero,
891  // but signed int will give comparison warnings)
892  for (unsigned int knotNum = N - 1; knotNum <= N - 1 && knotNum > 0; --knotNum){
893  // std::cout << " determining Legendre coefficients for knot #" << knotNum << std::endl;
894  auto & thisKnotCoeffs = legendreCoeffs[knotNum];
895  auto a_j = thisKnotCoeffs[0];
896  auto a_jPlus1 = knotNum+1 == N ? wgtsBySigma.rbegin()->second : legendreCoeffs[knotNum+1][0];
897  auto c_jPlus1 = knotNum+1 == N ? 0 : legendreCoeffs[knotNum+1][2];
898  auto h_j = h[knotNum];
899  thisKnotCoeffs[2] = z[knotNum] - mu[knotNum] * c_jPlus1;
900  thisKnotCoeffs[1] = (a_jPlus1 - a_j)/h_j - h_j*(c_jPlus1 + 2*thisKnotCoeffs[2])/3;
901  thisKnotCoeffs[3] = (c_jPlus1 - thisKnotCoeffs[2]) / (3*h_j);
902  }
903 
904  return legendreCoeffs;
905  }
906 
907  //----------------------------------------------------------------------------
908  // be sure that the interpolation function here matches the decomposition in CalcInterpFitCoeffs()...
910  std::map<float, float> const& wgtsBySigma)
911  {
912  auto coeffs = this->CalcInterpFitCoeffs(wgtsBySigma);
913  //LOG_VERBATIM("Weighters") << " Sigma = " << sigma;
914 
915  // There are N sets of coefficients for Legendre polynomials.
916  // We insisted in CalcInterpFitCoeffs() that they be separated
917  // by exactly 1 sigma, so the have x-values
918  // -N/2 sigma, -N/2+1 sigma, ..., 0 sigma, ..., N/2-1 sigma, N/2 sigma
919  // We need to determine which knot is the correct one for the requested number of sigma,
920  // which enables us to use the correct spline segment for interpolation.
921  unsigned int knotNum;
922  double maxKnotXval = coeffs.size()/2.;
923 
924  // if outside the range, use the end knot
925  // (the spline extrapolates linearly from the slope at the end knots.
926  // look up 'natural' boundary conditions for cubic splines
927  // if this scares you.)
928  if (sigma > maxKnotXval - 1) // -1 because the correct knot is the lower bound (see next comment)
929  knotNum = coeffs.size()-1;
930  else if (sigma < -maxKnotXval)
931  knotNum = 0;
932  // otherwise, find the correct knot.
933  // the right set of coefficients
934  // is the one that corresponds to the knot
935  // whose x-val is immediately below the requested sigma.
936  else
937  {
938  // x-vals are in units of sigma.
939  // note that this is really
940  // floor(sigma - (-maxKnotXval))
941  knotNum = std::floor(sigma + maxKnotXval);
942  }
943 
944  const auto & cs = coeffs[knotNum];
945  // the coefficients are for Legendre polynomials,
946  // which are of the form sum_{n} (a_n * (x-x_n)^n)
947  // since I insisted the sigmas be separated by 1 here,
948  // the x_n just counts from the lowest value up to
949  // the selected knot number by integers.
950  double xprime = sigma - (-maxKnotXval + knotNum);
951  double wgt = cs[0] + cs[1] * xprime + cs[2] * util::sqr(xprime) + cs[3] * util::cube(xprime);
952  /*if(sigma != 0 && sigma != 1){
953  std::cout << " for sigma = " << sigma << ":";
954  std::cout << " chose knot #" << knotNum+1 << " of " << coeffs.size() << std::endl;
955  std::cout << " coeffs are: ";
956  for (const auto c : cs)
957  std::cout << c << " ";
958  std::cout << std::endl;
959  std::cout << " x-x_i = " << xprime << std::endl;
960  std::cout << " interpolated weight = " << wgt << std::endl;
961  std::cout << "wgt = " << wgt << std::endl;
962  }*/
963 
964  return wgt;
965  }
966 
967  //=============================================================================//
968  //----------------------------- GENERAL WEIGHTS -------------------------------//
969  //=============================================================================//
970 
971  //----------------------------------------------------------------------------
973  {
975  }
976 
977  //----------------------------------------------------------------------------
979  {
980 
981  //File Service to save ND ratio plots
983  art::TFileDirectory ratios_dir = fTFS->mkdir("ND_Prediction_ratios");
984 
985  TAxis reweightAxis(480, 0, 30);
986  TAxis * pAxis = &reweightAxis;
987 
988  // We need to determine the ratio as a function of true energy for each
989  // epoch in the ND with the numu selection
990  std::unordered_map<int, double> dataPOT;
991  std::unordered_map<int, double> mcPOT;
992  std::unordered_map<int, double> dataEvents;
993  std::unordered_map<int, double> mcEvents;
994 
995  // figure out the number of bins for the reco energy histograms
996  // and then increase that by a factor of 4
997  auto axis = pAxis;
998  int bins = axis->GetNbins(); //4. * axis->GetNbins();
999  double low = axis->GetXmin();
1000  double high = axis->GetXmax();
1001 
1002  // at some point we will want to keep track of the data and MC epoch by epoch
1003  // but for now, we will just lump them all together
1004  //std::unordered_map<int, TH1F> dataReco;
1005  //std::unordered_map<int, TH2F> mcTrueVsReco;
1006  // make the histograms for this epoch
1007  TH1F dataReco("dataReco",
1008  ";Energy (GeV);Events",
1009  bins, low, high);
1010  TH2F mcTrueVsReco("mcTrueVsReco",
1011  ";True Energy (GeV);Reco Energy (GeV)",
1012  bins, low, high,
1013  bins, low, high);
1014 
1015  double potWeight = 1.;
1016  double wgt = 1.;
1017 
1018  // get the data POT for each epoch, only care about the ND
1019  // numu selection here.
1020  for(auto const& itr : *pEventListMap){
1021  auto md = itr.first;
1022 
1023  if(md.detector == novadaq::cnv::kNEARDET &&
1024  md.IsNuMuSelected()){
1025 
1026  if(!md.isMC){
1027  dataPOT.emplace(md.epoch, itr.second.ListSpillSummary().goodPOT);
1028  dataEvents.emplace(md.epoch, itr.second.size());
1029  }
1030  else{
1031  if(mcPOT.count(md.epoch) < 1){
1032  mcPOT .emplace(md.epoch, 0.0);
1033  mcEvents.emplace(md.epoch, 0.0);
1034  }
1035  mcPOT[md.epoch] += itr.second.ListSpillSummary().goodPOT;
1036  mcEvents[md.epoch] += itr.second.size();
1037  } // end if MC
1038  } // end if the ND and nu mu selection
1039  } // end loop over the event list map
1040 
1041  /*
1042  //Debugging output
1043  for(auto e : dataPOT){
1044 
1045  LOG_VERBATIM("CorrSpec_PropDecomp")
1046  << e.first
1047  << " **** DATA events / POT : "
1048  << dataEvents[e.first] << " / " << dataPOT[e.first]
1049  << " = " << dataEvents[e.first]/dataPOT[e.first];
1050  LOG_VERBATIM("CorrSpec_PropDecomp")
1051  << e.first
1052  << " **** MC events / MC : "
1053  << mcEvents[e.first] << " / " << mcPOT[e.first]
1054  << " = " << mcEvents[e.first]/mcPOT[e.first];
1055 
1056 
1057  }*/
1058 
1059  std::unique_ptr<fnex::DataVarVals> dataVals;
1060  std::unique_ptr<fnex::MCVarVals> mcVals;
1061  for(auto const& itr : *pEventListMap){
1062  auto md = itr.first;
1063 
1064  // We only want numu selected events from the ND to make the
1065  // prediction for the FD signal
1066  if(md.detector != novadaq::cnv::kNEARDET ||
1067  !md.IsNuMuSelected()) continue;
1068 
1069  auto dataPOTItr = dataPOT.find(md.epoch);
1070  if(dataPOTItr == dataPOT.end()){
1071  LOG_VERBATIM("ShifterAndWeighter")
1072  << md.ToString()
1073  << " ***** cannot find a corresponding number of data POT for this epoch"
1074  << " will not use this MC list";
1075  continue;
1076  }
1077 
1078  potWeight = dataPOTItr->second/itr.second.ListSpillSummary().goodPOT;
1079 
1080  // fill the histograms
1081  for(auto ev : itr.second ){
1082 
1083  // make a copy of the event variables so that any shifting we do is
1084  // not applied again (ie twice) in the Spectrum
1085  dataVals.reset(new fnex::DataVarVals(*(ev->dataVals)));
1086  if(ev->mcVals)
1087  mcVals.reset(new fnex::MCVarVals(*(ev->mcVals)));
1088 
1089  // TODO: handle the shifts for the event here
1090  // shifts of the reco neutrino energy are handled by just calling
1091  // for the fnex::kNu_RecoE variable of the dataVals because
1092  // that method in turn requests the amount of energy shift from this
1093  // singleton, so just make sure the stored input point is current before
1094  // calling this method
1095 
1096  if(!md.isMC){
1097 
1098  wgt = dataVals->val_at(fnex::kFake_Weight, md);
1099  // Sweet mother of mercy, never remove the following line
1100  // KM Sep 14 2016
1101  dataReco.Fill(dataVals->val_at(fnex::kNu_RecoE, md), wgt);
1102  }
1103  else{
1104  if(mcVals.get() == nullptr)
1105  throw cet::exception("AddNearTrueEnergyRatioWeighter")
1106  << "MC event does not have a valid MCVarVals unique_ptr, bail";
1107 
1108  this->SetCurrentEvent(ev, md);
1109 
1110  wgt = this->TotalWeightFromWeighters();
1111 
1112  mcTrueVsReco.Fill(mcVals->val_at(fnex::kTrueE),
1113  dataVals->val_at(fnex::kNu_RecoE, md),
1114  wgt * potWeight);
1115  }//end if mc events
1116  } // end loop over events in this list
1117  } // end loop over event lists
1118 
1119  dataVals.reset(nullptr);
1120  mcVals .reset(nullptr);
1121 
1122  // Now find the ratio and add the weights
1123  // The weights should be stored in the fSigExtrapWeights set
1124 
1125  // The weight is a ratio which is a function of the true energy.
1126  // The numerator for the weight is found by taking the sum over all
1127  // reconstructed energies contributing to the desired true energy.
1128  // The sum done over the product of the ND data to MC for each reco
1129  // energy bin with the number of events in the desired set of reco and
1130  // true energy bins.
1131  // The denominator is the number of MC events in the given true energy bin
1132 
1133  double data_mc_ratio = 0.;
1134  std::vector<double> numuFDPredict(bins, 0.);
1135 
1136  int bin = 0;
1137  double signalPredWgt = 0.;
1138 
1139  auto hreco = mcTrueVsReco.ProjectionY();
1140  auto htrue = mcTrueVsReco.ProjectionX();
1141 
1142  if(fSaveRatios){
1143  // Save the ND Data as a function of reco. neutrino energy
1144  TH1F datareco(dataReco);
1145  datareco.SetName("NDData_Reco");
1146  ratios_dir.make<TH1F>(dataReco);
1147 
1148  // Save the ND true vs reco before the ND Data/MC reco reweighting
1149  TH2F mcTrueVsReco_before_rwt(mcTrueVsReco);
1150  mcTrueVsReco_before_rwt.SetName("NDmcTrueVsReco");
1151  ratios_dir.make<TH2F>(mcTrueVsReco_before_rwt);
1152 
1153  // Save the ND MC as a function of reco. neutrino energy
1154  htrue->SetName("NDMC_True");
1155  ratios_dir.make<TH1D>(*htrue);
1156 
1157  // Save the ND MC as a function of true neutrino energy
1158  hreco->SetName("NDMC_Reco");
1159  ratios_dir.make<TH1D>(*hreco);
1160 
1161  // Save the ND Data/MC as a function of reco. neutrino energy
1162  TH1D* nd_data_over_mc = (TH1D*)datareco.Clone("nd_data_over_mc");
1163  nd_data_over_mc ->SetName("ND_Data_over_MC_Reco");
1164  nd_data_over_mc ->Divide(hreco);
1165  ratios_dir.make<TH1D>(*nd_data_over_mc);
1166  }
1167  // Weight the MC reco values by the Data/MC ratio for each reco bin
1168  // and sum the contributions from all reco bins into each truth bin
1169  for(int tbin = 1; tbin < bins + 1; ++tbin){
1170 
1171  // loop over the the data and mc histograms to get the
1172  // data/mc ratios and weight the 2D histogram with the ratio
1173  // for all true values
1174  for(int rbin = 1; rbin < bins + 1; ++rbin){
1175  bin = mcTrueVsReco.GetBin(tbin, rbin);
1176 
1177  data_mc_ratio = (hreco->GetBinContent(rbin) > 0) ? dataReco.GetBinContent(rbin) / hreco->GetBinContent(rbin) : 0.;
1178 
1179  numuFDPredict[tbin - 1] += data_mc_ratio * mcTrueVsReco.GetBinContent(bin);
1180 
1181  /*if(tbin == 10)
1182  LOG_VERBATIM("CorrSpec_PropDecomp")
1183  << low + (rbin * (high - low)/bins)
1184  << " Data/MC = "
1185  << data_mc_ratio
1186  << " "
1187  << numuFDPredict[tbin - 1];*/
1188 
1189 
1190  }//end of loop over reco energy bins
1191 
1192  // now the prediction in a given true energy bin is the
1193  // ratio of the weighted sum we just found to the raw MC expectation
1194  signalPredWgt = (htrue->GetBinContent(tbin) > 0.) ? numuFDPredict[tbin - 1] / htrue->GetBinContent(tbin) : 0.;
1195 
1196  LOG_DEBUG("CorrSpec_PropDecomp")
1197  << low + (tbin * (high - low)/bins)
1198  << " numuFDPred = "
1199  << numuFDPredict[tbin - 1]
1200  << " true MC = "
1201  << htrue->GetBinContent(tbin)
1202  << " signal pred wgt = "
1203  << signalPredWgt;
1204 
1205  fSigExtrapWeights.emplace(low + (tbin * (high - low)/bins), signalPredWgt);
1206  }//end loop over true energy bins
1207 
1208  //Save histograms used to make nue predictions
1209  if(fSaveRatios){
1210  //Make a histograms with nue weights and save it
1211  TH1F NDSignalWeightsTrueE("NDSignalWeightsTrueE",
1212  ";True Energy (GeV); ND^{Prediction}/ND^{MC}",
1213  bins, low, high);
1214 
1215  //loop over fSigExtrapWeights and fill the histogram
1216  for(auto it : fSigExtrapWeights) NDSignalWeightsTrueE.Fill(it.eTrueLow, it.weight);
1217  ratios_dir.make<TH1F>(NDSignalWeightsTrueE);
1218 
1219  }//only make and fill histograms once
1220  fSaveRatios = false;
1221 
1222  return;
1223  }
1224 
1225  //==========================================================================//
1226  //------------------ SYSTEMATICS USED IN BOTH ANALYSES ---------------------//
1227  //==========================================================================//
1228 
1229  //---------------------------------------------------------------------------
1231  fnex::SelectionType_t const& sel)
1232  {
1233  if (sel == fnex::kNuMuSelection &&
1234  !fCurrentMD.IsNuMuSelected() ) return 1.;
1235  else if(sel == fnex::kNuESelection &&
1236  !fCurrentMD.IsNuESelected() ) return 1.;
1237 
1238  auto const& sigma = fnex::SystematicSigmas::Instance()->Sigma(wgtType);
1239 
1240  LOG_DEBUG("NormBothWeight")
1241  << "Sigma value: "
1242  << sigma;
1243 
1245  }
1246 
1247  //--------------------------------------------------------------------------
1249  {
1250  // Someone already called us
1251  if(!fBeamSyst2017.fBeamSystHists.empty()) return;
1252 
1254 
1255  // TODO: make this a tunable parameter
1256  std::string shortName("enuMF_FHC_totErr");
1257 
1258  // constructor decides if initialized value is a path or an environment variable
1259  cet::search_path sp("FW_SEARCH_PATH");
1260  if ( !sp.find_file(std::string("CAFAna/data/beam/TABeamSyst.root"), fileName) ) {
1261  throw cet::exception("BeamSystFileError")
1262  << "cannot find Beam Systematics file CAFAna/data/beam/TABeamSyst.root "
1263  << fileName
1264  << " bail ungracefully\n"
1265  << "(RJN Guess: Did you run setup for this release?)\n"
1266  << __FILE__ << ":" << __LINE__;
1267  }
1268 
1269  fBeamSyst2017.setFileAndHistNames(fileName, shortName);
1271  }
1272 
1273  //----------------------------------------------------------------------------
1275  {
1276  // Someone already called us
1277  if(!fVecFluxPrincipals.empty()) return;
1278 
1279  std::string fileNameCV;
1280 
1281  std::string shortNameCV("ppfx_cv");
1282 
1283  // constructor decides if initialized value is a path or an environment variable
1284  cet::search_path sp("FW_SEARCH_PATH");
1285  if ( !sp.find_file(std::string("CAFAna/data/beam/ppfx_smooth_multiverse_weights_2018.root"), fileNameCV) ) {
1286  throw cet::exception("BeamSystFileError")
1287  << "cannot find Beam Systematics file CAFAna/data/beam/ppfx_smooth_multiverse_weights_2018.root "
1288  << fileNameCV
1289  << " bail ungracefully\n"
1290  << "(RJN Guess: Did you run setup for this release?)\n"
1291  << __FILE__ << ":" << __LINE__;
1292  }
1293 
1294  fPPFXSmoothCVWgt.setFileAndHistNames(fileNameCV,shortNameCV);
1296 
1297 
1298  std::string fileNamePCA;
1299  std::string shortNamePCA;
1300  if ( !sp.find_file(std::string("CAFAna/data/beam/ppfx_hadp_beam_pc_shifts_fn_2018.root"), fileNamePCA) ) {
1301  throw cet::exception("BeamSystFileError")
1302  << "cannot find Beam Systematics file CAFAna/data/beam/ppfx_hadp_beam_pc_shifts_fn_2018.root "
1303  << fileNamePCA
1304  << " bail ungracefully\n"
1305  << "(RJN Guess: Did you run setup for this release?)\n"
1306  << __FILE__ << ":" << __LINE__;
1307  }
1308 
1309  for(int index=0;index<5;index++) {
1310  ///For now 5 is just an arbitrary number, should be configurable
1311  std::string PCIdxStr = TString::Format("%02d", index).Data();
1312  shortNamePCA="ppfx_hadp_beam_pc"+PCIdxStr;
1313  fVecFluxPrincipals.push_back(BeamSyst(fileNamePCA,shortNamePCA));
1314  fVecFluxPrincipals[index].LoadBeamSystHists();
1315  }
1316 
1317  //Can also add a similar loop for the univeses
1318  //will leave this commented out for now
1319  // std::string shortNameUniv;
1320  // for(int index=0;index<100;index++) {
1321  // std::string UniIdxStr = TString::Format("%02d", index).Data();
1322  // shortNameUniv="ppfx_univ"+UniIdxStr;
1323  // fVecPPFXUniverses.push_back(BeamSyst(fileNameCV,shortNameUniv));
1324  // }
1325 
1326 
1327  }
1328 
1329  //--------------------------------------------------------------------------
1331  {
1332  //Default constructor, may need to zero some things
1333  }
1334 
1335  //-------------------------------------------------------------------------
1337  {
1338  // Someone already called us
1339  if(!fCalibSystHists.empty()) return;
1340 
1341  LOG_DEBUG("ShifterAndWeighter")
1342  << "LoadCalibSystHists";
1343 
1344  //const std::string kCalibSystFile = "FNEX/data/calibration/CalibrationSystematics2017.root";
1345 
1346  //loop over FHC and RHC
1347  for (int curr = 0; curr < kNumCurrents; ++curr){
1348 
1349  std::string fileCalib;
1350  std::string tempfileCalib;
1351  // if(curr == kFHC)
1352  tempfileCalib = CalibFilename;
1353 
1354  LOG_DEBUG("ShifterAndWeighter")
1355  << "curr = " << curr
1356  << ". Looking for file "
1357  << tempfileCalib;
1358 
1359  cet::search_path sp("FW_SEARCH_PATH");
1360  if(!sp.find_file(tempfileCalib, fileCalib)){
1361  throw cet::exception("CalibSystematicFileError")
1362  << "Cannot find calib systematics file "
1363  << tempfileCalib
1364  << " bail ungracefully\n\n"
1365  << __FILE__ << ":" << __LINE__;
1366  }
1367 
1368  LOG_DEBUG("ShifterAndWeighter")
1369  << "Found "
1370  << tempfileCalib
1371  << " as "
1372  << fileCalib;
1373 
1374  TFile fin(fileCalib.c_str(), "read");
1375  if(fin.IsZombie()){
1376  throw cet::exception("Weighters")
1377  << "Warning: couldn't open "
1378  << fileCalib
1379  << ". Crashing";
1380  }
1381 
1382  LOG_DEBUG("ShifterAndWeighter")
1383  << "Opened file, apparently";
1384 
1385  //Check that desired systematic is in file
1386  // Loop through the detector folder, over all histograms
1387  TIter iterDir(gDirectory->GetListOfKeys());
1388  TKey *keyDir;
1389 
1390  std::vector <std::string> systNames{"CalibShape",
1391  "LightLevel",
1392  "Cherenkov",
1393  "AbsCalibration",
1394  "RelCalibration"};
1395 
1396  for(auto systName : systNames){
1397  LOG_DEBUG("ShifterAndWeighter")
1398  << "systName = " << systName
1399  << ", kNumCurrents = " << kNumCurrents
1400  << ", kNumDets = " << kNumDets
1401  << ", kNumQuantiles = " << kNumQuantiles
1402  << ", kNumSigns = " << kNumSigns
1403  << ", kNumInteractions = " << kNumInteractions;
1404  fCalibSystHists[systName].resize(kNumCurrents);
1405  fCalibSystHists[systName][curr].resize(kNumDets);
1406  for(size_t d = 0; d < fCalibSystHists[systName][curr].size(); ++d){
1407  fCalibSystHists[systName][curr][d].resize(kNumQuantiles);
1408  for(size_t q = 0; q < fCalibSystHists[systName][curr][d].size(); ++q){
1409  fCalibSystHists[systName][curr][d][q].resize(kNumSigns);
1410  for(size_t pm = 0; pm < fCalibSystHists[systName][curr][d][q].size(); ++pm){
1411  fCalibSystHists[systName][curr][d][q][pm].resize(kNumInteractions);
1412  }
1413  }
1414  }
1415  }
1416 
1417  while((keyDir = (TKey *) iterDir())){
1418  TString dirName = keyDir->GetName(); // Get directory name
1419  LOG_DEBUG("ShifterAndWeighter")
1420  << "Looking in directory " << dirName;
1421  // for(auto systName : systNames) {
1422  // if(dirName.Contains(systName)){ //is it the desired syst
1423  TDirectory *dir = fin.GetDirectory(dirName);
1424  TIter iterSystDir(dir->GetListOfKeys());
1425  TKey *keySystDir;
1426  while((keySystDir = (TKey *) iterSystDir())){
1427  TString systDirName = keySystDir->GetName();
1428  LOG_DEBUG("ShifterAndWeighter")
1429  << "\tLooking in directory " << systDirName;
1430  TDirectory *systDir = dir->GetDirectory(systDirName);
1431  for(auto systName : systNames){
1432  if(systDirName.Contains(systName)){
1433  TIter iterHist(systDir->GetListOfKeys());
1434  TKey *keyHist;
1435 
1436  while((keyHist = (TKey *) iterHist())){
1437  TString histName = keyHist->GetName();
1438  LOG_DEBUG("ShifterAndWeighter")
1439  << "histName is " << histName;
1440 
1441  int det = 0;
1442  if(histName.Contains("ND")) det = kND;
1443  if(histName.Contains("FD")) det = kFD;
1444 
1445  int sign = 0;
1446  if(histName.Contains("minus")) sign = kMinus;
1447  if(histName.Contains("plus")) sign = kPlus;
1448 
1449  int quant = 0;
1450  if(histName.Contains("quant1")) quant = kQ1;
1451  if(histName.Contains("quant2")) quant = kQ2;
1452  if(histName.Contains("quant3")) quant = kQ3;
1453  if(histName.Contains("quant4")) quant = kQ4;
1454 
1455  int interaction = 0;
1456  if(histName.Contains("numuCC") ) interaction = kNumu;
1457  if(histName.Contains("nueCC") ) interaction = kNue;
1458  if(histName.Contains("anumuCC")) interaction = kAntiNumu;
1459  if(histName.Contains("anueCC") ) interaction = kAntiNue;
1460  if(histName.Contains("NC") ) interaction = kNC;
1461 
1462  LOG_DEBUG("ShifterAndWeighter")
1463  << "Get histogram " << histName;
1464  // grab the histogram
1465  TH1D *hist = (TH1D *) dir->Get(histName)->Clone();
1466 
1467  if(hist){
1468  // disassociate histogram from file
1469  hist->SetDirectory(0);
1470  //store relevant histograms
1471  LOG_DEBUG("ShifterAndWeighter")
1472  << "Store as fCalibSystHists "
1473  << systName
1474  << ", "
1475  << curr
1476  << ", "
1477  << det
1478  << ", "
1479  << quant
1480  << ", "
1481  << sign
1482  << ", "
1483  << interaction;
1484 
1485  fCalibSystHists[systName][curr][det][quant][sign][interaction] = hist;
1486  LOG_DEBUG("ShifterAndWeighter")
1487  << "Stored";
1488  }
1489  else{
1490  throw cet::exception("Weighters")
1491  << "Calibration systematic histogram "
1492  << histName
1493  << " doesn't seem to exist. Oh no.";
1494  } // end if statement checking histogram exists
1495  } // end loop through histograms in dir
1496  } // end check of hist name for syst
1497  } // end loop thru syst names
1498  } // end while loop through directories
1499  }//end loop over horn currents
1500  }
1501  LOG_VERBATIM("CalibSyst") << "Calib histograms loaded!";
1502  }
1503  //-------------------------------------------------------------------------
1505  {
1506  // Someone already called us
1507  if(!fCalibSystHists.empty()) return;
1508 
1509  LOG_VERBATIM("ShifterAndWeighter")
1510  << "Load2018CalibSystHists";
1511 
1512  std::string fileCalib;
1513  // std::string tempfileCalib;
1514  // tempfileCalib = CalibFilename;
1515 
1516  LOG_DEBUG("ShifterAndWeighter")
1517  << "Looking for calib syst file "
1518  << CalibFilename;
1519 
1520  cet::search_path sp("FW_SEARCH_PATH");
1521  if(!sp.find_file(CalibFilename, fileCalib)){
1522  throw cet::exception("CalibSystematicFileError")
1523  << "Cannot find calib systematics file "
1524  << CalibFilename
1525  << " bail ungracefully\n\n"
1526  << __FILE__ << ":" << __LINE__;
1527  }
1528 
1529  LOG_DEBUG("ShifterAndWeighter")
1530  << "Found "
1531  << CalibFilename
1532  << " as "
1533  << fileCalib;
1534 
1535  TFile fin(fileCalib.c_str(), "read");
1536  if(fin.IsZombie()){
1537  throw cet::exception("Weighters")
1538  << "Warning: couldn't open "
1539  << fileCalib
1540  << ". Crashing";
1541  }
1542 
1543  LOG_DEBUG("ShifterAndWeighter")
1544  << "Opened file, apparently";
1545 
1546  //Check that desired systematic is in file
1547  // Loop through the detector folder, over all histograms
1548  TIter iterDir(gDirectory->GetListOfKeys());
1549  TKey *keyDir;
1550  TKey *keySystDir;
1551 
1552  std::vector <std::string> systNames{"CalibShape",
1553  "LightLevel",
1554  "Cherenkov",
1555  "AbsCalibration",
1556  "RelCalibration"};
1557 
1558  for(auto systName : systNames){
1559  LOG_DEBUG("ShifterAndWeighter")
1560  << "systName = " << systName
1561  << ", kNumCurrents = " << kNumCurrents
1562  << ", kNumSels = " << kNumSels
1563  << ", kNumQuantiles = " << kNumQuantiles
1564  << ", kNumSigns = " << kNumSigns
1565  << ", kNumInteractions = " << kNumInteractions;
1566  //resize to fhc/rhc
1567  fCalibSystHists[systName].resize(kNumCurrents);
1568  //resize to nue/numu selection
1569  for(size_t c = 0; c < fCalibSystHists[systName].size(); ++c){
1570  fCalibSystHists[systName][c].resize(kNumSels);
1571  //resize to quantile OR low/mid/high/peri cvn
1572  for(size_t s = 0; s < fCalibSystHists[systName][c].size(); ++s){
1573  fCalibSystHists[systName][c][s].resize(kNumQuantiles);
1574  //resize to positive/negative sigma shift
1575  for(size_t q = 0; q < fCalibSystHists[systName][c][s].size(); ++q){
1576  fCalibSystHists[systName][c][s][q].resize(kNumSigns);
1577  //resize to number of different true interactions
1578  for(size_t pm = 0; pm < fCalibSystHists[systName][c][s][q].size(); ++pm){
1579  fCalibSystHists[systName][c][s][q][pm].resize(kNumInteractions);
1580  }
1581  }
1582  }
1583  }
1584  }
1585 
1586  //will loop over numuFHC, numuRHC, nueFHC, nueRHC
1587  while((keyDir = (TKey *) iterDir())){
1588  TString dirName = keyDir->GetName(); // Get directory name
1589  LOG_DEBUG("ShifterAndWeighter")
1590  << "Looking in directory " << dirName;
1591  TDirectory *dir = fin.GetDirectory(dirName);
1592  TIter iterSystDir(dir->GetListOfKeys());
1593 
1594  while((keySystDir = (TKey *) iterSystDir())){
1595  TString systDirName = keySystDir->GetName(); // Get directory name
1596  LOG_DEBUG("ShifterAndWeighter")
1597  << "\tin subdirectory "
1598  << systDirName;
1599 
1600  for(auto systName : systNames){
1601  LOG_DEBUG("ShifterAndWeighter")
1602  << "\tlooking for syst "
1603  << systName;
1604  if(systDirName.Contains(systName)){ //is it the desired syst
1605  LOG_DEBUG("ShifterAndWeighter")
1606  << "\tFound " << systName;
1607  TDirectory *systDir = dir->GetDirectory(systDirName);
1608  LOG_DEBUG("ShifterAndWeighter")
1609  << "\tGot directory " << systDirName;
1610  TIter iterHist(systDir->GetListOfKeys());
1611  TKey *keyHist;
1612 
1613  while((keyHist = (TKey *) iterHist())){
1614  TString histName = keyHist->GetName();
1615  LOG_DEBUG("ShifterAndWeighter")
1616  << "histName is "
1617  << histName;
1618 
1619  // int det = 0;
1620  // if(histName.Contains("ND")) det = kND;
1621  // if(histName.Contains("FD")) det = kFD;
1622 
1623  int curr = 0;
1624  if(dirName.Contains("FHC")) curr = kFHC;
1625  if(dirName.Contains("RHC")) curr = kRHC;
1626 
1627  int sel = 0;
1628  if(dirName.Contains("Nue")) sel = kNueSel;
1629  if(dirName.Contains("Numu")) sel = kNumuSel;
1630 
1631  int quant = 0;
1632  //numu quantiles
1633  if(histName.Contains("quant1")) quant = kQ1;
1634  if(histName.Contains("quant2")) quant = kQ2;
1635  if(histName.Contains("quant3")) quant = kQ3;
1636  if(histName.Contains("quant4")) quant = kQ4;
1637  //nue cvn bins
1638  if(histName.Contains("Low")) quant = kQ1;
1639  if(histName.Contains("Mid")) quant = kQ2;
1640  if(histName.Contains("High")) quant = kQ3;
1641  if(histName.Contains("Peripheral")) quant = kQ4;
1642 
1643  int sign = 0;
1644  if(histName.Contains("minus")) sign = kMinus;
1645  if(histName.Contains("plus")) sign = kPlus;
1646 
1647  int interaction = 0;
1648  if(histName.Contains("numuCC")) interaction = kNumu;
1649  if(histName.Contains("nueCC")) interaction = kNue;
1650  if(histName.Contains("anumuCC")) interaction = kAntiNumu;
1651  if(histName.Contains("anueCC")) interaction = kAntiNue;
1652  if(histName.Contains("NC")) interaction = kNC;
1653 
1654  LOG_DEBUG("ShifterAndWeighter")
1655  << "Get histogram " << histName;
1656  // grab the histogram
1657  TH1D *hist = (TH1D *) systDir->Get(histName)->Clone();
1658 
1659  if(hist){
1660  // disassociate histogram from file
1661  hist->SetDirectory(0);
1662  //store relevant histograms
1663  LOG_DEBUG("ShifterAndWeighter")
1664  << "Store as fCalibSystHists "
1665  << systName
1666  << ", "
1667  << curr
1668  << ", "
1669  << sel
1670  << ", "
1671  << quant
1672  << ", "
1673  << sign
1674  << ", "
1675  << interaction;
1676 
1677  fCalibSystHists[systName][curr][sel][quant][sign][interaction] = hist;
1678  LOG_DEBUG("ShifterAndWeighter")
1679  << "Stored";
1680  }
1681  else{
1682  throw cet::exception("Weighters")
1683  << "Calibration systematic histogram "
1684  << histName
1685  << " doesn't seem to exist. Oh no.";
1686  } // end if statement checking histogram exists
1687  } // end loop through histograms in dir
1688  } // end check of hist name for syst
1689  } // end loop thru syst names
1690  } // end while loop through syst directories
1691  }//end while loop through syst directories
1692 
1693  LOG_DEBUG("CalibSyst") << "Calib histograms loaded!";
1694  }
1695 
1696  //----------------------------------------------------------------------
1697  std::map<float, float> CalibSyst::CalcCalibSystWeights(std::string const& systName,
1698  fnex::MetaData currentMD,
1699  double trueE)
1700  {
1701  //FHC or RHC?
1702  int curr = currentMD.BeamType();
1703  //if it's an unknown beam type, default to using FHC
1704  if(curr > kRHC) curr = kFHC;
1705  LOG_DEBUG("ShifterAndWeighter")
1706  <<"Starting to CalcCalibSystWeights for "
1707  << systName
1708  << ", "
1709  << curr;
1710 
1711  // Need different histogram per detector, don't bother doing
1712  // anything if we have a bogus detector
1713  // int det = currentMD.detector;
1714  // if(det != novadaq::cnv::kNEARDET && det != novadaq::cnv::kFARDET) {
1715  // LOG_DEBUG("CalibSystFileError")
1716  // << "Wrong detector; ignore";
1717  // return { {-1.0, 1.0}, {0.0, 1.0}, {1.0, 1.0} };
1718  // }
1719  // else if(det == novadaq::cnv::kNEARDET) det = kND;
1720  // else if(det == novadaq::cnv::kFARDET ) det = kFD;
1721  // LOG_DEBUG("ShifterAndWeighter")
1722  // << "Detector " << det;
1723 
1724  // Need different histograms based on selection
1725  int sel = 0;
1726  int subsel = -999.;//quantile or cvn
1727  if(currentMD.IsNuMuSelected()){
1728  sel = kNumuSel;
1729  // find which quantile this event is in
1730  if (currentMD.selectionType == fnex::kNuMuSelectionQ1) subsel = kQ1;
1731  else if(currentMD.selectionType == fnex::kNuMuSelectionQ2) subsel = kQ2;
1732  else if(currentMD.selectionType == fnex::kNuMuSelectionQ3) subsel = kQ3;
1733  else if(currentMD.selectionType == fnex::kNuMuSelectionQ4) subsel = kQ4;
1734  else {
1735  LOG_DEBUG("CalibSystFileError")
1736  << "Selection type "
1738  << " is not a quantile; ignore";
1739  return {{-1.0, 1.0},
1740  {0.0, 1.0},
1741  {1.0, 1.0}};
1742  }
1743  }
1744  else if(currentMD.IsNuESelected()){
1745  sel = kNueSel;
1746  // find which cvn bin this event is in
1747  if (currentMD.selectionType == fnex::kNuESelectionLowPID) subsel = kQ1;
1748  else if(currentMD.selectionType == fnex::kNuESelectionMidPID) subsel = kQ2;
1749  else if(currentMD.selectionType == fnex::kNuESelectionHighPID) subsel = kQ3;
1750  else if(currentMD.selectionType == fnex::kNuESelectionPeripheral) subsel = kQ4;
1751  else {
1752  LOG_DEBUG("CalibSystFileError")
1753  << "Selection type "
1755  << " is not a CVN bin; ignore";
1756  return {{-1.0, 1.0},
1757  {0.0, 1.0},
1758  {1.0, 1.0}};
1759  }
1760 
1761  }
1762 
1763  // if it's any of the interactions we're not interested in then just return 1
1764  int interaction = currentMD.interactionType;
1765  if(interaction == fnex::kCosmicMuon ||
1766  interaction == fnex::kNuTauCC ||
1767  interaction == fnex::kNuTauBarCC )
1768  return { {-1.0, 1.0}, {0.0, 1.0}, {1.0, 1.0} };
1769  else if(interaction == fnex::kNuMuCC ) interaction = kNumu;
1770  else if(interaction == fnex::kNuMuBarCC) interaction = kAntiNumu;
1771  else if(interaction == fnex::kNuECC ) interaction = kNue;
1772  else if(interaction == fnex::kNuEBarCC ) interaction = kAntiNue;
1773  else if(interaction == fnex::kNC ) interaction = kNC;
1774  else {
1775  LOG_DEBUG("CalibSystFileError")
1776  << "Not expecting interaction type"
1777  << interaction
1778  << ", so we're ignoring it.";
1779  return { {-1.0, 1.0}, {0.0, 1.0}, {1.0, 1.0} };
1780  }
1781 
1782  // there are no NC histograms for CalibrationShape with the
1783  // nue selection so if we have that combination, just return
1784  // no change
1785  if(sel == kNueSel &&
1786  interaction == kNC &&
1787  systName == "CalibShape")
1788  return { {-1.0, 1.0}, {0.0, 1.0}, {1.0, 1.0} };
1789 
1790  // find the right histogram for the detector and sign
1791  double scale[2] = {1, 1};
1792  TH1D* hPlusMinus[2] = {0, 0};
1793 
1794  int sign = kPlus;
1795  while(sign >= kMinus){
1796 
1797  // cherenkov and calibshape only have histograms for +1 sigma
1798  if(sign == kMinus && (systName == "Cherenkov" || systName == "CalibShape")) {
1799  scale[sign] = 1.0;
1800  break;
1801  }
1802 
1803  LOG_DEBUG("ShifterAndWeighter")
1804  << "Looking for fCalibSystHists[" << systName
1805  << "][" << curr
1806  << "][" << sel
1807  << "][" << subsel
1808  << "][" << sign
1809  << "][" << interaction
1810  << "]";
1811  hPlusMinus[sign] = fCalibSystHists.find(systName)->second[curr][sel][subsel][sign][interaction];
1812 
1813  if(hPlusMinus[sign] == 0){
1814  LOG_WARNING("CalibSystFileError")
1815  << ": Can't find desired histogram for sign: "
1816  << sign
1817  << " systematic: "
1818  << systName
1819  << "[" << curr
1820  << "][" << sel
1821  << "][" << subsel
1822  << "][" << sign
1823  << "][" << interaction
1824  << "], ignore";
1825  // return { {-1.0, 1.0}, {0.0, 1.0}, {1.0, 1.0} };
1826  scale[sign] = 1.0;
1827  }
1828 
1829  //Only use weights if the energy is in the range of the histogram
1830  else if (trueE > hPlusMinus[sign]->GetXaxis()->GetXmin() &&
1831  trueE < hPlusMinus[sign]->GetXaxis()->GetXmax() ){
1832  scale[sign] = hPlusMinus[sign]->Interpolate(trueE);
1833 
1834  // assume a scale of 0 really means that there just
1835  // isn't a change to be applied
1836  if(scale[sign] == 0.) scale[sign] = 1.;
1837  }
1838  --sign;
1839  }
1840  LOG_DEBUG("ShifterAndWeighter")
1841  << "Neg scale: " << scale[kMinus]
1842  << "\t Pos scale: " << scale[kPlus];
1843 
1844  // Only want light level systematic to be one-sided but don't want to totally
1845  // restructure code since minus > plus for it
1846  if(systName == "LightLevel") {
1847  scale[kPlus] = scale[kMinus];
1848  scale[kMinus] = 1.0;
1849  }
1850 
1851  return { {-1.0, scale[kMinus]}, {0.0, 1.0}, {1.0, scale[kPlus]} };
1852  }
1853 
1854  //-------------------------------------------------------------------------
1856  novadaq::cnv::DetId const& det)
1857  {
1858  if(fCurrentMD.detector != det) return 1.0;
1859 
1860  // need to butcher these strings to get the weights by sigma
1861  std::string systName = wgtType.substr(0, wgtType.size() - 3);
1862  // if (det == novadaq::cnv::kNEARDET) systName.erase(0, 4);
1863  // else if(det == novadaq::cnv::kFARDET ) systName.erase(0, 3);
1864 
1865  double trueE = fCurrentEvent->mcVals->val_at(fnex::kTrueE);
1866  auto wgtBySigma = fCalibSystHists.CalcCalibSystWeights(systName, fCurrentMD, trueE);
1867 
1868  double wgt = this->CalcLinearInterpWgt(fCurrentPoint.ParameterValue(std::make_pair(wgtType, det)),
1869  wgtBySigma);
1870 
1871  // if(fCurrentMD.selectionType == fnex::kNuESelectionPeripheral){
1872  // for(auto const& itr : wgtBySigma){
1873  // LOG_VERBATIM("ShifterAndWeighter")
1874  // << systName
1875  // << " sigma "
1876  // << itr.first
1877  // << " scale "
1878  // << itr.second
1879  // << " total wgt "
1880  // << wgt;
1881  // }
1882  // }
1883 
1884  return wgt;
1885  }
1886 
1887  //------------------------------------------------------------------------
1889  fCalibSystHists.Load2018CalibSystHists(CalibFilename);
1890  }
1891 
1892  //--------------------------------------------------------------------------
1894  {
1895  //Default constructor, may need to zero some things
1896  }
1897 
1898  //-------------------------------------------------------------------------
1900  {
1901  // Someone already called us
1902  if(!fGeniePCAHists.empty()) return;
1903 
1904  const std::string kGeniePCASystFile = "CAFAna/data/xs/genie_small_pc_shifts_fn_2018.root";
1905 
1906  std::string fileGeniePCA;
1907 
1908  cet::search_path sp("FW_SEARCH_PATH");
1909  if ( !sp.find_file(kGeniePCASystFile, fileGeniePCA) ) {
1910  throw cet::exception("GeniePCASystematicFileError")
1911  << "Cannot find Genie PCA systematics file "
1912  << fileGeniePCA
1913  << " bail ungracefully\n\n"
1914  << __FILE__ << ":" << __LINE__;
1915  }
1916 
1917  TFile fin(fileGeniePCA.c_str(), "read");
1918  if (fin.IsZombie()){
1919  throw cet::exception("Weighters")
1920  << "Warning: couldn't open "
1921  << fileGeniePCA
1922  << ". Crashing!";
1923  }
1924 
1925  // resize histogram vectors appropriately
1926  // (too many nested for loops? i have no idea what you're talking about)
1927  fGeniePCAHists.resize(kNumFlux);
1928  for(size_t f = 0; f < fGeniePCAHists.size(); ++f) {
1929  fGeniePCAHists[f].resize(kNumDets);
1930  for(size_t d = 0; d < fGeniePCAHists[f].size(); ++d) {
1931  fGeniePCAHists[f][d].resize(kNumSels);
1932  for(size_t s = 0; s < fGeniePCAHists[f][d].size(); ++s) {
1933  fGeniePCAHists[f][d][s].resize(kNumInteractions);
1934  for(size_t i = 0; i < fGeniePCAHists[f][d][s].size(); ++i) {
1935  fGeniePCAHists[f][d][s][i].resize(kNumSigns);
1936  }
1937  }
1938  }
1939  }
1940 
1941  //Check that desired systematic is in file
1942  // Loop through the detector folder, over all histograms
1943  TIter iterDir(gDirectory->GetListOfKeys());
1944  TKey* keyDir;
1945  int foundHist = 0;
1946 
1947  while((keyDir = (TKey*)iterDir())) {
1948  TString histName = keyDir->GetName();
1949  if(histName.Contains(shortName)) {
1950  foundHist++;
1951 
1952  int flux = 0;
1953  if(histName.Contains("FHC")) flux = kFHC;
1954  if(histName.Contains("RHC")) flux = kRHC;
1955 
1956  int det = 0;
1957  if(histName.Contains("ND")) det = kND;
1958  if(histName.Contains("FD")) det = kFD;
1959 
1960  int sel = 0;
1961  if(histName.Contains("Numu")) sel = kNumuSel;
1962  if(histName.Contains("Nue" )) sel = kNueSel;
1963 
1964  int flav = 0;
1965  if(histName.Contains("SigQE" )) flav = kSigQE;
1966  if(histName.Contains("SigNonQE")) flav = kSigNonQE;
1967  if(histName.Contains("Other" )) flav = kBkgFlav;
1968  if(histName.Contains("NC" )) flav = kNC;
1969 
1970  int sign = 0;
1971  if(histName.Contains("minus")) sign = kMinus;
1972  if(histName.Contains("plus" )) sign = kPlus;
1973 
1974  // grab relevant histogram
1975  TH1D *hist = (TH1D*)fin.Get(histName)->Clone();
1976 
1977  if(hist) {
1978  // store relevant histograms
1979  fGeniePCAHists[flux][det][sel][flav][sign] = hist;
1980 
1981  // disassociate from file
1982  fGeniePCAHists[flux][det][sel][flav][sign]->SetDirectory(0);
1983  } else {
1984  throw cet::exception("Weighters")
1985  << "Genie PCA histogram "
1986  << histName
1987  << " does not exist! Something has gone terribly wrong!";
1988  } // end if statement checking hist exists
1989  } // end hist name check
1990  } // end while loop
1991 
1992  LOG_DEBUG("GeniePCASysts") << "Genie PCA histograms loaded!";
1993  }
1994 
1995  //----------------------------------------------------------------------
1996  std::map<float, float> GeniePCASyst::CalcGeniePCASystWeights(fnex::MetaData currentMD, int mode, double energy)
1997  {
1998  // Need different histogram per detector, don't bother doing
1999  // anything if we have a bogus detector
2000  int det = currentMD.detector;
2001  if(det != novadaq::cnv::kNEARDET && det != novadaq::cnv::kFARDET) {
2002  LOG_VERBATIM("GeniePCASystFileError")
2003  << det
2004  << " is invalid detector; ignore";
2005  return { {-1.0, 1.0}, {0.0, 1.0}, {1.0, 1.0} };
2006  }
2007  else if(det == novadaq::cnv::kNEARDET) det = kND;
2008  else if(det == novadaq::cnv::kFARDET ) det = kFD;
2009 
2010  // different histograms for beam type
2011  int beamType = currentMD.BeamType();
2012  if(beamType != fnex::kFHC && beamType != fnex::kRHC) {
2013  LOG_VERBATIM("GeniePCASystFileError")
2014  << beamType
2015  << " not a valid beam type! IGNORE.";
2016  return { {-1.0, 1.0}, {0.0, 1.0}, {1.0, 1.0} };
2017  }
2018  else if(beamType == fnex::kFHC) beamType = kFHC;
2019  else if(beamType == fnex::kRHC) beamType = kRHC;
2020 
2021  // different histograms for different selections
2022  int sel = currentMD.selectionType;
2023  if(!currentMD.IsNuMuSelected() && !currentMD.IsNuESelected()) {
2024  LOG_VERBATIM("GeniePCASystFileError")
2025  << sel
2026  << " not a valid selection type! IGNOOORE.";
2027  return { {-1.0, 1.0}, {0.0, 1.0}, {1.0, 1.0} };
2028  }
2029  else if(currentMD.IsNuMuSelected()) sel = kNumuSel;
2030  else if(currentMD.IsNuESelected ()) sel = kNueSel;
2031 
2032  // if it's any of the interactions we're not interested in then just return 1
2033  int interaction = currentMD.interactionType;
2034  if (currentMD.IsNuMuSelected() && interaction == fnex::kNuMuCC &&
2035  mode == simb::kQE)
2036  interaction = kSigQE;
2037  else if(currentMD.IsNuESelected() && interaction == fnex::kNuECC &&
2038  mode == simb::kQE)
2039  interaction = kSigQE;
2040  else if(currentMD.IsNuMuSelected() && interaction == fnex::kNuMuCC &&
2041  mode != simb::kQE)
2042  interaction = kSigNonQE;
2043  else if(currentMD.IsNuESelected() && interaction == fnex::kNuECC &&
2044  mode != simb::kQE)
2045  interaction = kSigNonQE;
2046  else if(interaction == fnex::kNC)
2047  interaction = kNC;
2048  else interaction = kBkgFlav;
2049 
2050  // for some reason background histograms don't exist for numu so ignore these
2051  if(interaction == kBkgFlav && sel == kNumuSel)
2052  return { {-1.0, 1.0}, {0.0, 1.0}, {1.0, 1.0} };
2053 
2054  // find the right histogram for the detector and sign
2055  double scale[2] = {1, 1};
2056  TH1D* hPlusMinus[2] = {0};
2057 
2058  // and now deal with the number of sigma
2059  int numSigns = kNumSigns;
2060 
2061  for(int sign = kPlus; sign < numSigns; ++sign){
2062  hPlusMinus[sign] = fGeniePCAHists[beamType][det][sel][interaction][sign];
2063 
2064  if(hPlusMinus[sign] == 0){
2065  LOG_WARNING("GeniePCASystFileError")
2066  << ": Can't find desired histogram; ignore";
2067  return { {-1.0, 1.0}, {0.0, 1.0}, {1.0, 1.0} };
2068  }
2069 
2070  // only use weights if the energy is in the range of the histogram
2071  else if (energy > hPlusMinus[sign]->GetXaxis()->GetXmin() &&
2072  energy < hPlusMinus[sign]->GetXaxis()->GetXmax() ){
2073  scale[sign] = hPlusMinus[sign]->Interpolate(energy);
2074 
2075  }
2076 
2077  }
2078  return { {-1.0, scale[kMinus]}, {0.0, 1.0}, {1.0, scale[kPlus]} };
2079  }
2080  //-------------------------------------------------------------------------
2082  {
2084  int mode = fCurrentEvent->mcVals->val_at(fnex::kTrueIntMode);
2085 
2086  auto wgtBySigma = fVecGeniePCASysts[index].CalcGeniePCASystWeights(fCurrentMD, mode, energy);
2087 
2088  std::string systName;
2089  std::string PCIdxStr = TString::Format("%01d", (int)index).Data();
2090  systName = "GeniePCA" + PCIdxStr + "Wgt";
2091 
2093  wgtBySigma);
2094  }
2095 
2096  //------------------------------------------------------------------------
2098  // someone already called us
2099  if(!fVecGeniePCASysts.empty()) return;
2100 
2101  for(int index = 0; index < 5; ++index) {
2102  // 5 chosen arbitrarily, should probably make configurable at some point
2103  std::string PCIdxStr = TString::Format("%02d", index).Data();
2104  std::string shortName = "genie_small_pc"+PCIdxStr;
2105  fVecGeniePCASysts.push_back(GeniePCASyst());
2106 
2107  fVecGeniePCASysts[index].LoadGeniePCASystHists(shortName);
2108  }
2109 
2110  }
2111 
2112  //-------------------------------------------------------------------------
2114  {
2115  //Default constructor, may need to zero some things
2116 
2117  }
2118 
2120  {
2122  fShortName=shortName;
2123  }
2124 
2125 
2127  {
2129  fShortName=shortName;
2130  }
2131 
2132 
2134  {
2135  // Someone already called us
2136  if(!fBeamSystHists.empty()) return;
2137 
2138  fBeamSystHists.resize(kNumFlux);
2139  for(size_t h = 0; h < fBeamSystHists.size(); ++h){
2140  fBeamSystHists[h].resize(kNumDets);
2141  for(size_t d = 0; d < fBeamSystHists[h].size(); ++d){
2142  fBeamSystHists[h][d].resize(kNumSigns);
2143  for(size_t pm = 0; pm < fBeamSystHists[h][d].size(); ++pm){
2144  fBeamSystHists[h][d][pm].resize(kNumFlavors);
2145  }
2146  }
2147  }
2148 
2149  //open file
2150  TFile fin (fFileName.c_str(),"read");
2151  if(fin.IsZombie()){
2152  LOG_WARNING("ShifterAndWeighter")
2153  << "Warning: couldn't open beam systematic file "
2154  << fFileName;
2155  return;
2156  }
2157 
2158  //Check that desired systematic is in file
2159  // Loop through the detector folder, over all histograms
2160  TIter iterHist(gDirectory->GetListOfKeys());
2161  TKey* keyHist;
2162  int foundHisto = 0;
2163  bool isSeparatedByFlavor = false;
2164 
2165  while((keyHist = (TKey*)iterHist())) {
2166  TString histName = keyHist->GetName(); // Get a histogram name
2167  if(histName.Contains(fShortName)){ //is it the desired syst
2168  ++foundHisto;
2169 
2170  int horn = 0;
2171  if(histName.Contains("FHC")) horn = kFHC;
2172  if(histName.Contains("RHC")) horn = kRHC;
2173 
2174  int det = 0;
2175  if(histName.Contains("ND")) det = kND;
2176  if(histName.Contains("FD")) det = kFD;
2177 
2178  int sign = 0;
2179  if(histName.Contains("minus")) sign = kMinus;
2180  if(histName.Contains("plus" )) sign = kPlus;
2181 
2182  int flav = 0;
2183  if(histName.Contains("numu" )) flav = kNumu;
2184  if(histName.Contains("nue" )) flav = kNue;
2185  if(histName.Contains("anumu")) flav = kAntiNumu;
2186  if(histName.Contains("anue" )) flav = kAntiNue;
2187 
2188  //store relevant histograms
2189  fBeamSystHists[horn][det][sign][flav] = (TH1D*) fin.Get(histName)->Clone();
2190  // disassociate it from the file it came from
2191  fBeamSystHists[horn][det][sign][flav]->SetDirectory(0);
2192  }
2193  if(histName.Contains("numu")){
2194  isSeparatedByFlavor=true;
2195  }
2196  }
2197 
2198  if (foundHisto == 0) {
2199  throw cet::exception("BeamSystFileError")
2200  << "Beam systematic "
2201  << fShortName
2202  << " is not in this file";
2203  }
2204 
2205  if (!isSeparatedByFlavor) {
2206  throw cet::exception("BeamSystFileError")
2207  << "Beam systematic "
2208  << fShortName
2209  << " doesn't have required flavor information";
2210  }
2211  fin.Close();
2212  }
2213 
2214  //----------------------------------------------------------------------------
2215  std::map<float, float> BeamSyst::CalcBeamSystWeights(novadaq::cnv::DetId const &curDet, fnex::BeamType_t const &beamType, double energy, int pdgOrig)
2216  {
2217  // Need different histogram per detector, don't bother doing
2218  // anything if we have a bogus detector
2219  int det = kNumDets;
2220  if(curDet == novadaq::cnv::kNEARDET)
2221  det = kND;
2222  else if(curDet == novadaq::cnv::kFARDET)
2223  det = kFD;
2224  else{
2225  LOG_VERBATIM("BeamSystFileError")
2226  << "Wrong detector; ignore";
2227  return { {-1.0, 1.0}, {0.0, 1.0}, {1.0, 1.0} };
2228  }
2229 
2230 
2231  //Find the neutrino flavor
2232  int flav = 0;
2233  switch(pdgOrig){
2234  case 14: flav = kNumu; break;
2235  case 12: flav = kNue; break;
2236  case -14: flav = kAntiNumu; break;
2237  case -12: flav = kAntiNue; break;
2238  case 13:
2239  //Cosmic muon
2240  return { {-1.0, 1.0}, {0.0, 1.0}, {1.0, 1.0} };
2241  default:
2242  LOG_WARNING("BeamSystFileError")
2243  << "Wrong nu flavor; "
2244  << pdgOrig
2245  << " ignore";
2246 
2247  return { {-1.0, 1.0}, {0.0, 1.0}, {1.0, 1.0}};
2248  }
2249 
2250  // Need different histogram per horn current
2251  int horn = kNumFlux;
2252  if(beamType == fnex::kFHC)
2253  horn = kFHC;
2254  else if(beamType == fnex::kRHC)
2255  horn = kRHC;
2256  else{
2257  LOG_VERBATIM("BeamSystFileError")
2258  << "Wrong horn type: "
2259  << beamType
2260  << "; ignore";
2261 
2262  return { {-1.0, 1.0}, {0.0, 1.0}, {1.0, 1.0} };
2263  }
2264 
2265  //Find the right histogram for the detector and sign
2266  double scale[2] = {1, 1};
2267  TH1D* hPlusMinus[2] = {0};
2268  for(int sign = kMinus; sign < kNumSigns; ++sign){
2269  hPlusMinus[sign] = fBeamSystHists[horn][det][sign][flav];
2270 
2271  if(hPlusMinus[sign] == 0){
2272  LOG_WARNING("BeamSystFileError")
2273  << ": Can't find desired histogram; ignore";
2274  return { {-1.0, 1.0}, {0.0, 1.0}, {1.0, 1.0} };
2275  }
2276 
2277  //Only use weights if the energy is in the range of the histogram
2278  else if (energy > hPlusMinus[sign]->GetXaxis()->GetXmin() &&
2279  energy < hPlusMinus[sign]->GetXaxis()->GetXmax() ){
2280  scale[sign] = hPlusMinus[sign]->Interpolate(energy);
2281  }
2282  }
2283 
2284  return { {-1.0, scale[kMinus]}, {0.0, 1.0}, {1.0, scale[kPlus]} };
2285  } // std::map<double, double> BeamSystWeightFn::CalcWeights(const fnex::Event&)
2286 
2287  //----------------------------------------------------------------------------
2289  novadaq::cnv::DetId const& det)
2290  {
2291  if(det != fCurrentMD.detector) return 1.;
2292 
2293  std::map<float, float> wgtBySigma;
2294 
2295  if(sysType.find("2017") != std::string::npos) {
2296  wgtBySigma = fBeamSyst2017.CalcBeamSystWeights(det,
2297  fCurrentMD.BeamType(),
2298  fCurrentEvent->mcVals->val_at(fnex::kTrueE),
2300  } else {
2301  wgtBySigma = fSABeamSyst.CalcBeamSystWeights(det,
2302  fCurrentMD.BeamType(),
2303  fCurrentEvent->mcVals->val_at(fnex::kTrueE), //True Energy
2304  fCurrentEvent->mcVals->val_at(fnex::kTruePDGOrig)); //Flavour
2305  }
2306 
2307 // return this->CalcInterpolatedWgt(fCurrentPoint.ParameterValue(std::make_pair(sysType, fCurrentMD.detector)),
2308 // wgtBySigma);
2310  wgtBySigma);
2311  }
2312 
2313  //----------------------------------------------------------------------------
2315  {
2316  if(index > fVecFluxPrincipals.size() ) return 1;
2317 
2318  /*if(fCurrentMD.BeamType() == fnex::kUnknownBeam &&
2319  fCurrentEvent->eventId) {
2320  if(fCurrentEvent->eventId->run)
2321  LOG_VERBATIM("SJB") << "UNKNOWN BEAM!"
2322  << "Current run: "
2323  << fCurrentEvent->eventId->run;
2324  else
2325  LOG_VERBATIM("SJB") << "Can't find run information! What!";
2326 
2327  }*/
2328 
2329  auto wgtBySigma = fVecFluxPrincipals[index].CalcBeamSystWeights(fCurrentMD.detector,
2330  fCurrentMD.BeamType(),
2331  fCurrentEvent->mcVals->val_at(fnex::kTrueE), //True Energy
2332  fCurrentEvent->mcVals->val_at(fnex::kTruePDGOrig)); //Flavour
2333 
2334  std::string systname;
2335  std::string PCIdxStr = TString::Format("%01d", (int)index).Data();
2336  systname="FluxPCA"+PCIdxStr+"Wgt";
2337 
2339  wgtBySigma);
2340  }
2341 
2342  //----------------------------------------------------------------------------
2344  {
2345  // may need to come up with default for cosmic muons, however I think the fact that
2346  // the default weights for every GENIE variable is a collection of 1's should
2347  // mean this is OK
2348  std::map<float, float> wgtBySigma = fCurrentEvent->mcVals->GENIEWgts(systName.substr(0, systName.size() - 3));
2349 
2350  // for(auto itr : wgtBySigma)
2351  //LOG_VERBATIM("GENIESystWeight")
2352  //<< systName
2353  //<< " weights "
2354  //<< itr.first
2355  //<< " "
2356  //<< itr.second;
2357 
2358 // return this->CalcInterpolatedWgt(fCurrentPoint.ParameterValue(std::make_pair(systName, fCurrentMD.detector)),
2359 // wgtBySigma);
2360 
2362  wgtBySigma);
2363 
2364  }
2365 
2366  //----------------------------------------------------------------------------
2368  {
2369  if((fCurrentEvent->mcVals->val_at(fnex::kEmpiricalMEC_Weight) != 1 ||
2372  ((!anti && fCurrentEvent->mcVals->val_at(fnex::kTruePDG) > 0) ||
2373  (anti && fCurrentEvent->mcVals->val_at(fnex::kTruePDG) < 0)) )
2374  {
2375 
2376  std::string sign;
2377  if(!anti) sign = "Nu";
2378  else sign = "AntiNu";
2379 
2380  auto nomWgt = fCurrentEvent->mcVals->val_at(fnex::kEmpiricalMEC_Weight);
2381  auto parVal = fCurrentPoint.ParameterValue(std::make_pair("MECShape"+sign+"Wgt", fCurrentMD.detector));
2382  double wgt = 1.;
2383 
2384  auto wgtVar = (parVal < 0) ?
2387 
2388  if(nomWgt > 0)
2389  wgt = 1 + std::abs(parVal) * (wgtVar / nomWgt - 1);
2390 
2391  // don't allow negative weights
2392  if(wgt < 0)
2393  wgt = 0;
2394 
2395  // don't allow crazy weights
2396  if(wgt > 10)
2397  wgt = 10;
2398 
2399  return wgt;
2400  }
2401 
2402  return 1.0;
2403  }
2404 
2405  //----------------------------------------------------------------------------
2407  {
2408  auto parVal = fCurrentPoint.ParameterValue(std::make_pair("RPARESWgt", fCurrentMD.detector));
2409  if(parVal > 0) return 1. + parVal * (1 - fCurrentEvent->mcVals->val_at(fnex::kRPARES_Weight));
2410 
2411  return 1.;
2412  }
2413 
2414  //----------------------------------------------------------------------------
2416  {
2417  auto parVal = fCurrentPoint.ParameterValue(std::make_pair("DISvnCC1piWgt", fCurrentMD.detector));
2418  if(parVal > 0) return 1. + parVal * (1 - fCurrentEvent->mcVals->val_at(fnex::kDISvnCC1pi_Weight));
2419 
2420  return 1.;
2421  }
2422 
2423  //----------------------------------------------------------------------------
2425  {
2428 
2429  // problem: for the SA, it was safe to assume that a 0sigma shift meant a weight of 1
2430  // BUT in A2017, the central value for RPACCQESupp/EnhWeights is the value of RPACCQEWeight
2431  std::map<float, float> wgtBySigma = fCurrentEvent->mcVals->GENIEWgts("RPACCQEshapeEnh");
2432  wgtBySigma.at(0.) = fCurrentEvent->mcVals->val_at(fnex::kRPACCQE_Weight);
2433 
2435  wgtBySigma) / wgtBySigma.find(0.)->second;
2436  }
2437 
2438  //----------------------------------------------------------------------------
2440  {
2441  // if(std::abs(fCurrentEvent->mcVals->val_at(fnex::kTrueIntType)) == 1. * simb::kMEC &&
2442  if((fCurrentEvent->mcVals->val_at(fnex::kEmpiricalMEC_Weight) != 1 ||
2445  ((!anti && fCurrentEvent->mcVals->val_at(fnex::kTruePDG) > 0) ||
2446  (anti && fCurrentEvent->mcVals->val_at(fnex::kTruePDG) < 0)) )
2447  {
2448 
2449  LOG_DEBUG("ShifterAndWeighter")
2450  << "Found MEC event with "
2451  << " MEC weight " << fCurrentEvent->mcVals->val_at(fnex::kEmpiricalMEC_Weight)
2452  << " MECtoGENIEQE " << fCurrentEvent->mcVals->val_at(fnex::kEmpiricalMECtoGENIEQE_Weight)
2453  << " MECtoGENIERES " << fCurrentEvent->mcVals->val_at(fnex::kEmpiricalMECtoGENIERES_Weight)
2454  << " PDG " << fCurrentEvent->mcVals->val_at(fnex::kTruePDG);
2455 
2456  std::string sign;
2457  if(!anti) sign = "Nu";
2458  else sign = "AntiNu";
2459 
2460  auto parVal = fCurrentPoint.ParameterValue(std::make_pair("MECEnuShape"+sign+"Wgt", fCurrentMD.detector));
2461  double nuE = fCurrentEvent->mcVals->val_at(fnex::kTrueE);
2462  if(nuE < 0) return 1.0;
2463 
2464  // this function is the 1 sigma bound around central value weight 1
2465  static const TF1 rwgtFn("MECRwgtFn", "1/(2.5*x+1)");
2466  double wgt = 1 + parVal * rwgtFn.Eval(nuE);
2467 
2468  // no negative weights
2469  if(wgt < 0)
2470  wgt = 0;
2471 
2472  LOG_DEBUG("ShifterAndWeighter")
2473  << "\tMECEnuShapeWgt " << wgt
2474  << " (Enu = " << nuE
2475  << ", rwgtFn.Eval(nuE) = " << rwgtFn.Eval(nuE)
2476  << ", parval = " << parVal << ")";
2477  return wgt;
2478  }
2479 
2480  return 1.0;
2481  }
2482 
2483  //----------------------------------------------------------------------------
2485  {
2486  // it doesn't make sense to ask for a CC weight if the event is an NC
2487  // or cosmic muon
2490 
2491  auto parVal = fCurrentPoint.ParameterValue(std::make_pair("MaCCQEReducedWgt", fCurrentMD.detector));
2492  std::map<float, float> wgtBySigma = fCurrentEvent->mcVals->GENIEWgts("MaCCQE");
2493 
2494  // mysterious 2018 correction factor (M_A^QE CV shift from 0.99 to 1.04)
2495  const double correctionInSigma = (1.04 - 0.99) / 0.25;
2496  double cv_weight = 1. + correctionInSigma * (wgtBySigma.at(1) - 1);
2497 
2498  double genie_sigma = (parVal > 0) ? parVal * (0.05 / 0.25) : parVal * (0.05 / 0.15);
2499 
2500  if( cv_weight > 0) {
2501  for(auto wgt : wgtBySigma) wgt.second /= cv_weight;
2502 
2503  double reduced_error = 0.05;
2504  double cv_ma = 1.04;
2505  double genie_ma = 0.99;
2506 
2507  double shifted_ma = cv_ma * ( 1 + parVal * reduced_error );
2508  double genie_frac_shift = ( shifted_ma - genie_ma ) / genie_ma;
2509  double genie_error = ( genie_frac_shift > 0 ) ? 0.25 : 0.15;
2510  genie_sigma = genie_frac_shift / genie_error;
2511  }
2512 
2513  return this->CalcLinearInterpWgt(genie_sigma, wgtBySigma);
2514  }
2515 
2516  //----------------------------------------------------------------------------
2518  {
2519  // it doesn't make sense to ask for a CC weight if the event is an NC
2520  // or cosmic muon
2523 
2524  return this->GENIESystWeight("MaCCRESWgt");
2525  }
2526 
2527  //----------------------------------------------------------------------------
2529  {
2530  // it doesn't make sense to ask for a CC weight if the event is an NC
2531  // or cosmic muon
2534 
2535  return this->GENIESystWeight("MvCCRESWgt");
2536  }
2537 
2538  //----------------------------------------------------------------------------
2540  {
2541  // it doesn't make sense to ask for a NC weight if the event isn't NC
2542  if(fCurrentMD.interactionType != fnex::kNC) return 1.;
2543 
2544  return this->GENIESystWeight("MaNCRESWgt");
2545  }
2546 
2547  //-----------------------------------------------------------------------------
2549  {
2550  // if event isn't coherent or is NC, return 1
2553  fCurrentMD.interactionType == fnex::kNC) return 1.;
2554 
2555  auto const& sigma = fnex::SystematicSigmas::Instance()->Sigma("COHCCScale2018Wgt");
2556 
2557  auto parVal = fCurrentPoint.ParameterValue(std::make_pair("COHCCScale2018Wgt", fCurrentMD.detector));
2558 
2559  double weight = 1. + sigma * parVal;
2560 
2561  LOG_DEBUG("COHCCScale") << "COHCCScale2018Weight is "
2562  << weight;
2563 
2564  if(weight >= 0.) return weight;
2565  else return 0.;
2566  }
2567 
2568  //----------------------------------------------------------------------------
2570  {
2571  // if event isn't coherent or isn't NC, return 1
2574  fCurrentMD.interactionType != fnex::kNC) return 1.;
2575 
2576  auto const& sigma = fnex::SystematicSigmas::Instance()->Sigma("COHNCScale2018Wgt");
2577 
2578  auto parVal = fCurrentPoint.ParameterValue(std::make_pair("COHNCScale2018Wgt", fCurrentMD.detector));
2579 
2580  double weight = 1. + sigma * parVal;
2581 
2582  LOG_DEBUG("COHNCScale") << "COHNCScale2018Weight is "
2583  << weight;
2584 
2585  if(weight >= 0.) return weight;
2586  else return 0.;
2587  }
2588 
2589  //-----------------------------------------------------------------------------
2591  {
2592  auto const& sigma = fnex::SystematicSigmas::Instance()->Sigma("TauScaleSystWgt");
2593 
2594  if(std::abs(fCurrentEvent->mcVals->val_at(fnex::kTruePDG)) == 16. &&
2595  fCurrentEvent->mcVals->val_at(fnex::kTrueCCNC) == 1. * simb::kCC ){
2596 
2597  double weight = 1. + sigma * fCurrentPoint.ParameterValue(std::make_pair("TauScaleSystWgt", fCurrentMD.detector));
2598 
2599  return std::max(0., weight);
2600  }
2601 
2602  return 1.0;
2603  }
2604 
2605  //==============================================================================//
2606  //----------------- SYSTEMATICS USED IN THE NUMU ANALYSIS ----------------------//
2607  //==============================================================================//
2608 
2609  //----------------------------------------------------------------------------
2611  {
2612  //This should ALWAYS be correlated with:
2613  // RPACCQEshapeSuppWeight
2614  // MECInitStateNPFrac
2615  // MECEnuShape
2616 
2617  double tempwgt = 1.;
2618  double parVal = fCurrentPoint.ParameterValue(std::make_pair("SumSmallXSecNumu2017Wgt", fCurrentMD.detector));
2619 
2620  //All the GENIE style wgts
2621  tempwgt *= this->GENIESystWeight("SumSmallXSecNumu2017Wgt");
2622  LOG_DEBUG("ShifterAndWeighter")
2623  << "SumSmallXSec: with GENIE wgts "
2624  << tempwgt;
2625 
2626  //DIS
2627  tempwgt *= 1. + parVal * (1 - fCurrentEvent->mcVals->val_at(fnex::kOtherDIS_Weight));
2628  if(fCurrentEvent->mcVals->val_at(fnex::kOtherDIS_Weight) != 1.)
2629  LOG_DEBUG("ShifterAndWeighter")
2630  << "\twith OtherDIS wgts "
2631  << tempwgt;
2632 
2633  //RadCorrNue and RadCorrNueBar
2634  //There should be a NuE version of this weight WITHOUT these systematics.
2635  //For joint analyses, USE THE NUE VERSION OF THE SUMSMALLXSEC WEIGHT
2636  if(fCurrentEvent->mcVals->val_at(fnex::kTrueCCNC) == 1. * simb::kCC &&
2637  (abs(fCurrentEvent->mcVals->val_at(fnex::kTruePDG))) == 12){
2638  tempwgt *= 1 + 0.02 * parVal;
2639  //SimpleSecondClassCurrentsSecSyst
2640  if(fCurrentEvent->mcVals->val_at(fnex::kTruePDG) == 12) tempwgt *= 1 + 0.02 * parVal;
2641  if(fCurrentEvent->mcVals->val_at(fnex::kTruePDG) == -12) tempwgt *= 1 - 0.02 * parVal;
2642  LOG_DEBUG("ShifterAndWeighter")
2643  << "\twith RadCorr wgts "
2644  << tempwgt;
2645  }
2646 
2647  return tempwgt;
2648  }
2649 
2650 
2651  //----------------------------------------------------------------------------
2653  {
2654  // if(std::abs(fCurrentEvent->mcVals->val_at(fnex::kTrueIntType)) == 1. * simb::kMEC &&
2655  if((fCurrentEvent->mcVals->val_at(fnex::kEmpiricalMEC_Weight) != 1 ||
2658  ((!anti && fCurrentEvent->mcVals->val_at(fnex::kTruePDG) > 0) ||
2659  (anti && fCurrentEvent->mcVals->val_at(fnex::kTruePDG) < 0)) )
2660  {
2661  std::string sign;
2662  if(!anti) sign = "Nu";
2663  else sign = "AntiNu";
2664 
2665  auto parVal = fCurrentPoint.ParameterValue(std::make_pair("MECInitNPFrac"+sign+"Wgt", fCurrentMD.detector));
2666  const double nominalNPfrac = 0.8; // see DocDB 18741
2667  double newNPfrac = nominalNPfrac + (0.1*parVal);
2668  double wgt = 1;
2669 
2670  // only physical between 0 and 1
2671  if(newNPfrac > 1)
2672  newNPfrac = 1;
2673  else if(newNPfrac < 0)
2674  newNPfrac = 0;
2675 
2676  if(fCurrentEvent->mcVals->val_at(fnex::kTrueHitNuc) == 2000000201)
2677  wgt = newNPfrac / nominalNPfrac;
2678  else
2679  wgt = (1 - newNPfrac) / (1 - nominalNPfrac);
2680 
2681  LOG_DEBUG("ShifterAndWeighter")
2682  << "\tMECInitStateNPFracWgt " << wgt
2683  << " (nominalNPfrac = " << nominalNPfrac
2684  << ", newNPfrac = " << newNPfrac
2685  << ", parval = " << parVal << ")";
2686  return wgt;
2687  }
2688 
2689  return 1.0;
2690  }
2691 
2692  //==============================================================================//
2693  //----------------- SYSTEMATICS USED IN THE NUE ANALYSIS ----------------------//
2694  //==============================================================================//
2695 
2696  //----------------------------------------------------------------------------
2698  {
2699  double parVal = fCurrentPoint.ParameterValue(std::make_pair("RadCorrNueWgt", fCurrentMD.detector));
2700 
2701  if(fCurrentEvent->mcVals->val_at(fnex::kTrueCCNC) == 1. * simb::kCC &&
2702  fCurrentEvent->mcVals->val_at(fnex::kTruePDG ) == 12 )
2703  return 1.0 + 0.02 * parVal;
2704  else return 1.0;
2705  }
2706 
2707  //----------------------------------------------------------------------------
2709  {
2710  double parVal = fCurrentPoint.ParameterValue(std::make_pair("RadCorrNueBarWgt", fCurrentMD.detector));
2711 
2712  if(fCurrentEvent->mcVals->val_at(fnex::kTrueCCNC) == 1. * simb::kCC &&
2713  fCurrentEvent->mcVals->val_at(fnex::kTruePDG ) == -12 )
2714  return 1.0 + 0.02 * parVal;
2715  else return 1.0;
2716  }
2717 
2718  //----------------------------------------------------------------------------
2720  {
2721  double parVal = fCurrentPoint.ParameterValue(std::make_pair("SecondClassCurrWgt", fCurrentMD.detector));
2722 
2723  if(fCurrentEvent->mcVals->val_at(fnex::kTrueCCNC) == 1. * simb::kCC &&
2724  fCurrentEvent->mcVals->val_at(fnex::kTruePDG ) == 12 )
2725  return 1.0 + 0.02 * parVal;
2726  else if(fCurrentEvent->mcVals->val_at(fnex::kTrueCCNC) == 1. * simb::kCC &&
2727  fCurrentEvent->mcVals->val_at(fnex::kTruePDG ) == -12 )
2728  return 1.0 - 0.02 * parVal;
2729  else return 1.0;
2730  }
2731 
2733  {
2734  // Apply weights to the ND as ND makes the FD signal extrapolation weights
2735  if(fCurrentMD.detector == novadaq::cnv::kFARDET) return 1.0;
2736 
2737  double parVal = fCurrentPoint.ParameterValue(std::make_pair("NueExtrapSig2017Wgt", fCurrentMD.detector));
2738 
2739  // Make sure it's a muon neutrino
2740  if(abs(fCurrentEvent->mcVals->val_at(fnex::kTruePDG)) != 14) return 1.0;
2741 
2742  // Check whether it's affecting nue background
2743 // if(fCurrentEvent->mcVals->val_at(fnex::kNuE_CVN) >= 0.5) return 1.0;
2744  if(fCurrentEvent->dataVals->val_at(fnex::kNuE_CVN, fCurrentMD) >= 0.5) return 1.0;
2745 
2746  // Grab the relevant histogram
2747  if(!fNueExtrapWgtHist) {
2749  std::string shortName("trueQ2_weight_AllSamples");
2750 
2751  // constructor decides if initialized value is a path or an environment variable
2752  cet::search_path sp("FW_SEARCH_PATH");
2753  if ( !sp.find_file(std::string("CAFAna/nue/NDtoFD_kin_combo_weights.root"), fileName) ) {
2754  throw cet::exception("NueExtrapSystFileError")
2755  << "cannot find Nue Extrap Systematics file CAFAna/nue/NDtoFD_kin_combo_weights.root "
2756  << fileName
2757  << " bail ungracefully\n"
2758  << "(RJN Guess: Did you run setup for this release?)\n"
2759  << __FILE__ << ":" << __LINE__;
2760  }
2761 
2762  TFile fin(fileName.c_str(), "read");
2763  if(fin.IsZombie()) {
2764  LOG_WARNING("ShifterAndWeighter")
2765  << "Warning: couldn't open nue extrapolation systematic file "
2766  << fileName;
2767  }
2768 
2769  fNueExtrapWgtHist = (TH1*)fin.Get(shortName.c_str())->Clone();
2770 
2771  // disassociate from file
2772  fNueExtrapWgtHist->SetDirectory(0);
2773  fin.Close();
2774  }
2775 
2776  // Get reconstructed Q^2
2777  double kinVar = fCurrentEvent->dataVals->val_at(fnex::kRecoQ2, fCurrentMD);
2778  if(kinVar == std::numeric_limits<float>::min())
2779  LOG_WARNING("ShifterAndWeighter")
2780  << "Something has gone wrong, kRecoQ2 should have a proper value.";
2781  double kinWeight = 1.0;
2782 
2783  if(kinVar >= fNueExtrapWgtHist->GetXaxis()->GetXmin() &&
2784  kinVar <= fNueExtrapWgtHist->GetXaxis()->GetXmax() )
2785  kinWeight = fNueExtrapWgtHist->GetBinContent(fNueExtrapWgtHist->GetXaxis()->FindBin(kinVar));
2786  else if(kinVar > fNueExtrapWgtHist->GetXaxis()->GetXmax())
2787  kinWeight = fNueExtrapWgtHist->GetBinContent(fNueExtrapWgtHist->GetNbinsX());
2788  else if(kinVar < fNueExtrapWgtHist->GetXaxis()->GetXmin())
2789  kinWeight = fNueExtrapWgtHist->GetBinContent(1);
2790 
2791  double weight = std::max(0., kinWeight - 1);
2792 
2793  return 1 + weight * parVal;
2794  }
2795 
2796  //---------------------------------------------------------------------------
2798  {
2799  // This weight only applies to the FD
2800  if(fCurrentMD.detector == novadaq::cnv::kNEARDET) return 1.0;
2801 
2802  double parVal = fCurrentPoint.ParameterValue(std::make_pair("NueExtrapBkg2017Wgt", fCurrentMD.detector));
2803 
2804  const double coreWgt = 0.012;
2805  const double periWgt = 0.013;
2806 
2807  // Make sure it's not numu->nue signal
2808  if(abs(fCurrentEvent->mcVals->val_at(fnex::kTruePDG )) == 12 &&
2809  abs(fCurrentEvent->mcVals->val_at(fnex::kTruePDGOrig)) == 14 &&
2810  fCurrentEvent ->mcVals->val_at(fnex::kTrueCCNC ) == 1. * simb::kCC )
2811  return 1.0;
2812 
2813  // Check whether it's the peripheral or core sample
2817  return 1.0 + coreWgt * parVal;
2819  return 1.0 + periWgt * parVal;
2820  else
2821  return 1.0;
2822  }
2823 
2824  //----------------------------------------------------------------------------
2826  {
2827  // This should ALWAYS be correlated with:
2828  // RPACCQEshapeSuppWeight
2829  // MECInitStateNPFrac
2830  // MECEnuShape
2831 
2832  double tempwgt = 1.;
2833  double parVal = fCurrentPoint.ParameterValue(std::make_pair("SumSmallXSecNue2017Wgt", fCurrentMD.detector));
2834 
2835  // All the GENIE style wgts (these are the same as for numu so just grab the numu object)
2836  tempwgt *= this->GENIESystWeight("SumSmallXSecNumu2017Wgt");
2837  LOG_DEBUG("ShifterAndWeighter")
2838  << "SumSmallXSec: with GENIE wgts "
2839  << tempwgt;
2840 
2841  // DIS
2842  tempwgt *= 1. + parVal * (1 - fCurrentEvent->mcVals->val_at(fnex::kOtherDIS_Weight));
2843  if(fCurrentEvent->mcVals->val_at(fnex::kOtherDIS_Weight) != 1.)
2844  LOG_DEBUG("ShifterAndWeighter")
2845  << "\twith OtherDIS wgts "
2846  << tempwgt;
2847 
2848  return tempwgt;
2849  }
2850 
2851  //============================================================================//
2852  //------------------------ OLD/DEPRECATED SYSTEMATICS ------------------------//
2853  //---------------- (OR SYSTEMATICS INCLUDED IN SUMMED SYSTEMATIC) ------------//
2854  //============================================================================//
2855 
2856  //----------------------------------------------------------------------------
2858  {
2859  auto const& sigma = fnex::SystematicSigmas::Instance()->Sigma("NCNormWgt");
2860  if(fCurrentEvent->mcVals->val_at(fnex::kTrueCCNC) == 1. * simb::kNC)
2861  return 1. + sigma * fCurrentPoint.ParameterValue(std::make_pair("NCNormWgt", fCurrentMD.detector));
2862 
2863  return 1.0;
2864  }
2865 
2866  //----------------------------------------------------------------------------
2868  {
2869  auto const& sigma = fnex::SystematicSigmas::Instance()->Sigma("MECNormWgt");
2870  if(std::abs(fCurrentEvent->mcVals->val_at(fnex::kTrueIntType)) == 1000.){
2871  return 1. + sigma * fCurrentPoint.ParameterValue(std::make_pair("MECNormWgt", fCurrentMD.detector));
2872  }
2873 
2874  return 1.0;
2875  }
2876 
2877  //----------------------------------------------------------------------------
2879  {
2880  // cosmic muon normalization is only for the far detector
2881  if(fCurrentMD.detector == novadaq::cnv::kNEARDET) return 1.;
2882 
2883  auto const& sigma = fnex::SystematicSigmas::Instance()->Sigma("CosmicMuNormWgt");
2885  return 1. + sigma * fCurrentPoint.ParameterValue(std::make_pair("CosmicMuNormWgt", fCurrentMD.detector));
2886  }
2887 
2888  return 1.0;
2889  }
2890 
2891  //----------------------------------------------------------------------------
2893  fnex::SelectionType_t const& sel)
2894  {
2895  // check for the correct selection for the current event with this function
2896  // sel is one of the two generic selection types, kNuMuSelection or kNuESelection
2897  if(fCurrentMD.IsNuMuSelected() &&
2898  sel == fnex::kNuMuSelection &&
2899  wgtType.find("Nue") != std::string::npos) return 1.;
2900  else if(fCurrentMD.IsNuESelected() &&
2901  sel == fnex::kNuESelection &&
2902  wgtType.find("Nue") == std::string::npos) return 1.;
2903 
2904  auto const& sigma = fnex::SystematicSigmas::Instance()->Sigma(wgtType);
2905 
2906  LOG_DEBUG("AbsNormWeight")
2907  << " Sigma Value"
2908  << sigma;
2909 
2911 
2912  }
2913 
2914  //----------------------------------------------------------------------------
2916  fnex::SelectionType_t const& sel)
2917  {
2918  // relative normalization is only for the far detector
2919  if(fCurrentMD.detector == novadaq::cnv::kNEARDET) return 1.;
2920 
2921  // sel is one of the two generic selection types, kNuMuSelection or kNuESelection
2922  bool wantNue = (sel == fnex::kNuESelection ) && (wgtType.find("Nue") != std::string::npos);
2923  bool wantNumu = (sel == fnex::kNuMuSelection) && (wgtType.find("Nue") == std::string::npos);
2924 
2925  // check for the correct selection for the current event with this function
2926  // return 1 if we have a numu selected event but want the nue relative normalization
2927  // or if we have a nue selected event but want the numu relative normalization
2928  if(fCurrentMD.IsNuMuSelected() && wantNue) return 1.;
2929  else if(fCurrentMD.IsNuESelected() && wantNumu) return 1.;
2930 
2931  auto const& sigma = fnex::SystematicSigmas::Instance()->Sigma(wgtType);
2932 
2934  }
2935 
2936  //----------------------------------------------------------------------------
2937  // Function to apply a normalization weight only to a given horn current
2938  // There is no check on the selection type for the current event, but that
2939  // could be added in the future
2941  fnex::BeamType_t const& beam)
2942  {
2943  // check the beam type is correct for the current event
2944  if(beam != fCurrentMD.BeamType() ) return 1.;
2945 
2946  auto const& sigma = fnex::SystematicSigmas::Instance()->Sigma(wgtType);
2947 
2948  LOG_DEBUG("HornNormWeight")
2949  << " Sigma Value "
2950  << sigma
2951  << " beam type "
2952  << beam;
2953 
2955  }
2956 
2957  //----------------------------------------------------------------------------
2959  {
2960  // rock muon normalization is only for the near detector
2961  if(fCurrentMD.detector == novadaq::cnv::kFARDET) return 1.;
2962 
2963  auto const& sigma = fnex::SystematicSigmas::Instance()->Sigma("RockNormOnFidWgt");
2964 
2965  return 1. + sigma * fCurrentPoint.ParameterValue(std::make_pair("RockNormOnFidWgt", fCurrentMD.detector));
2966  }
2967 
2968  //----------------------------------------------------------------------------
2969  // This one is not in SystematicParameters.fcl, so comment it out for now.
2970 // double ShifterAndWeighter::NearIntensityNormWeight(novadaq::cnv::DetId const& det)
2971 // {
2972 // auto const& sigma = fnex::SystematicSigmas::Instance()->Sigma("IntensityNormWgt");
2973 //
2974 // return 1. + sigma * fCurrentPoint.ParameterValue(std::make_pair("IntensityNormWgt", det));
2975 // }
2976 
2977  //----------------------------------------------------------------------------
2979  novadaq::cnv::DetId const& det)
2980  {
2981  if(det != fCurrentMD.detector) return 1.;
2982 
2983  auto const& sigma = fnex::SystematicSigmas::Instance()->Sigma(wgtType);
2984 
2986  }
2987 
2988  //----------------------------------------------------------------------------
2990  novadaq::cnv::DetId const& det)
2991  {
2992  if(det != fCurrentMD.detector) return 1.;
2993 
2994  auto const& sigma = fnex::SystematicSigmas::Instance()->Sigma(wgtType);
2995 
2997  }
2998 
2999  //----------------------------------------------------------------------------
3001  {
3002  // this one only applies to the far detector
3003  if(fCurrentMD.detector == novadaq::cnv::kNEARDET) return 1.;
3004 
3005  auto const& sigma = fnex::SystematicSigmas::Instance()->Sigma("NoiseFarSystWgt");
3006 
3007  return 1. + sigma * fCurrentPoint.ParameterValue(std::make_pair("NoiseFarSystWgt", fCurrentMD.detector));
3008  }
3009 
3010  //----------------------------------------------------------------------------
3012  {
3013  auto parVal = fCurrentPoint.ParameterValue(std::make_pair("RPACCQEWgt", fCurrentMD.detector));
3014  if(parVal > 0) return 1. + parVal * (1 - fCurrentEvent->mcVals->val_at(fnex::kRPACCQE_Weight));
3015 
3016  return 1.;
3017  }
3018 
3019  //----------------------------------------------------------------------------
3021  {
3024 
3025  //See comment in RPACCQEshapeEnhWeight()
3026  std::map<float, float> wgtBySigma = fCurrentEvent->mcVals->GENIEWgts("RPACCQEshapeSupp");
3027  wgtBySigma.at(0.) = fCurrentEvent->mcVals->val_at(fnex::kRPACCQE_Weight);
3028 
3029  double wgt = this->CalcLinearInterpWgt(fCurrentPoint.ParameterValue(std::make_pair("RPACCQEshapeSuppWgt", fCurrentMD.detector)),
3030  wgtBySigma) / wgtBySigma.find(0.)->second;
3031  // LOG_VERBATIM("ShifterAndWeighter")
3032  // << "RPACCQEshapeSuppWeight = " << wgt
3033  // << ", from (" << wgtBySigma.at(-1.)
3034  // << ", " << wgtBySigma.at( 0.)
3035  // << ", " << wgtBySigma.at( 1.)
3036  // << ")";
3037  return wgt;
3038  }
3039 
3040  //----------------------------------------------------------------------------
3042  {
3045  auto const& sigma = fnex::SystematicSigmas::Instance()->Sigma("CosmicOverlayNormWgt");
3046 
3047  return 1. + sigma * fCurrentPoint.ParameterValue(std::make_pair("CosmicOverlayNormWgt", fCurrentMD.detector));
3048  }
3049 
3050  return 1.0;
3051  }
3052 
3053  //----------------------------------------------------------------------------
3055  {
3056 
3057  SigExtrapWgt nsw(fCurrentEvent->mcVals->val_at(fnex::kTrueE), 0.);
3058 
3059  // find the location of the true energy in the set
3060  auto loc = fSigExtrapWeights.lower_bound(nsw);
3061 
3062  if(loc == fSigExtrapWeights.end()){
3063  LOG_DEBUG("NueSignalExtrapWeightFn")
3064  << nsw.eTrueLow
3065  << ": energy is not found in weight set, returning 0";
3066  return nsw.weight;
3067  }
3068 
3069  return loc->weight;
3070  }
3071 
3072  //----------------------------------------------------------------------------
3074  {
3075 
3076  SigExtrapWgt nsw(fCurrentEvent->mcVals->val_at(fnex::kTrueE), 1.);
3077 
3078  // find the location of the true energy in the set
3079  auto loc = fSigExtrapWeights.lower_bound(nsw);
3080 
3081  if(loc == fSigExtrapWeights.end()){
3082  LOG_DEBUG("NuMuSignalExtrapWeightFn")
3083  << nsw.eTrueLow
3084  << ": energy is not found in weight set, returning 1";
3085  return nsw.weight;
3086  }
3087 
3088  return loc->weight;
3089  }
3090 
3091  //----------------------------------------------------------------------------
3093  {
3094  // Someone already called us
3095  if(!fNueSystHists.empty()) return;
3096 
3097  // TODO: Ideally, the path and file names for these files would be a tunable
3098  // parameter and not hardcoded
3099  const std::string kNueSASystDir = "CAFAna/nue/SecondAna/Systs";
3100  std::vector<std::string> fileNames;
3101  // We currently aren't using the Birks B systematic
3102  //fileNames.push_back(kNueSASystDir + "/nueSA_syst_BirksB.root");
3103  fileNames.push_back(kNueSASystDir + "/nueSA_syst_BirksC.root");
3104  fileNames.push_back(kNueSASystDir + "/nueSA_syst_ShiftX.root");
3105  fileNames.push_back(kNueSASystDir + "/nueSA_syst_ShiftY.root");
3106  fileNames.push_back(kNueSASystDir + "/nueSA_syst_ShiftXYAbs.root");
3107  fileNames.push_back(kNueSASystDir + "/nueSA_syst_ShiftXYRel.root");
3108 
3109  std::vector<std::string> systNames{"NueBirksCSystCVNWgt",
3110  "NueShiftXSystCVNWgt",
3111  "NueShiftYSystCVNWgt",
3112  "NueXYAbsSystCVNWgt",
3113  "NueXYRelSystCVNWgt"};
3114 
3115  // for now just use the CVN pid histograms because that is all NOvA currently
3116  // uses
3117  std::string pid("CVN");
3119  for(size_t f = 0; f < fileNames.size(); ++f){
3120  // std::cout << "opening file: " << fFileName.c_str() << "\n";
3121 
3122  cet::search_path sp("FW_SEARCH_PATH");
3123  if ( !sp.find_file(std::string(fileNames[f]), histFile) ) {
3124  throw cet::exception("NueSystematicFileError")
3125  << "cannot find nue systematics file "
3126  << fileNames[f]
3127  << " bail ungracefully\n\n"
3128  << __FILE__ << ":" << __LINE__;
3129  }
3130 
3131  TFile fin(histFile.c_str(), "read");
3132  if (fin.IsZombie()){
3133  throw cet::exception("Weighters")
3134  << "Warning: couldn't open "
3135  << histFile
3136  << ". Crashing";
3137  }
3138 
3139  std::vector<std::string> channels = {"numutonumucc", "numutonuecc", "numutonutaucc",
3140  "nuetonumucc", "nuetonuecc", "nuetonutaucc", "nc"};
3141  const std::vector<int> sigmas = {-3, -2, -1, 0, +1, +2, +3};
3142  const std::vector<std::string> sigmaStrs = {"minus3", "minus2", "minus1", "zero",
3143  "plus1", "plus2", "plus3"};
3144 
3145  for (size_t i = 0; i < channels.size(); ++i){
3146  std::vector< std::pair<int, TH1D*> > curHists;
3147  for (size_t j = 0; j < sigmas.size(); ++j){
3148  // Find histogram for given OscChannel at j sigma
3149  std::string hName = TString::Format("fdratio_%s_%s_%s",
3150  pid.c_str(),
3151  channels[i].c_str(),
3152  sigmaStrs[j].c_str()).Data();
3153  TH1D* h = (TH1D*)fin.Get(hName.c_str());
3154  // std::cout << hName.c_str() << "\t" << h << "\n";
3155  if(!h && (j == 3 || j == 4)){
3156  // Surely it's a bad thing if we can't find nom, +-1 sigma hists
3157  throw cet::exception("FNEXWeighters")
3158  << "Error: can't find necessary "
3159  << hName
3160  << " histogram in file "
3161  << histFile
3162  << ". Crashing";
3163  }
3164  if(!h) continue; // Missing other shifts is fine
3165  h->SetDirectory(0);
3166  curHists.emplace_back(sigmas[j], h);
3167  } // end loop over sigmas
3168 
3169  fNueSystHists[systNames[f]].push_back(curHists);
3170  }// end loop over channels
3171  } // end loop over file names
3172 
3173  return;
3174  }
3175 
3176  //----------------------------------------------------------------------
3178  {
3179  int pdg = fCurrentEvent->mcVals->val_at(fnex::kTruePDG );
3180  int pdgOrig = fCurrentEvent->mcVals->val_at(fnex::kTruePDGOrig);
3181  int isCC = fCurrentEvent->mcVals->val_at(fnex::kTrueCCNC );
3182 
3183  if (!isCC) return NueSystsSAW::OscChannel::kNC;
3184  if (std::abs(pdgOrig) == 12){
3185  if (std::abs(pdg) == 12) return NueSystsSAW::OscChannel::kNuEToNuECC;
3186  if (std::abs(pdg) == 14) return NueSystsSAW::OscChannel::kNuEToNuMuCC;
3187  if (std::abs(pdg) == 16) return NueSystsSAW::OscChannel::kNuEToNuTauCC;
3188  }
3189  if (std::abs(pdgOrig) == 14){
3190  if (std::abs(pdg) == 12) return NueSystsSAW::OscChannel::kNuMuToNuECC;
3191  if (std::abs(pdg) == 14) return NueSystsSAW::OscChannel::kNuMuToNuMuCC;
3192  if (std::abs(pdg) == 16) return NueSystsSAW::OscChannel::kNuMuToNuTauCC;
3193  }
3194 
3195  throw cet::exception("NueSystematicFromHistogram")
3196  << "Unknown Oscillation Channel";
3197 
3199  }
3200 
3201  //----------------------------------------------------------------------
3202  std::map<float, float> ShifterAndWeighter::CalcNueSystWeights(std::string const& systName)
3203  {
3204  //Before we go anything else we have to check if it is a cosmic muon because the cosmic muon is missing
3205  //some of the variables that we rely upon and causes exceptions to be thrown
3206  int pdgOrig = fCurrentEvent->mcVals->val_at(fnex::kTruePDGOrig);
3207 
3208  if(pdgOrig == 13) {
3209  //Cosmic muon
3210  return { {-1.0, 1.0}, {0.0, 1.0}, {1.0, 1.0} };
3211  }
3212 
3213  // make a fake MetaData for this function, it should be MC with the nue selection
3216  md.isMC = true;
3217 
3218  //Now we need to calculate the weight for this event
3219  //So we need to know what the OscChannel is
3220  //And then calcuate the weight from the correct energy+pid bin
3221  // TODO: at some point we may need to be more specific about the selection type
3223  double nueenergy = fCurrentEvent->dataVals->val_at(fnex::kNu_RecoE, md);
3224  double pid = fCurrentEvent->dataVals->val_at(fnex::kNuE_CVN, md);
3225 
3226  auto chanHist = fNueSystHists.find(systName)->second;
3227  const int bin = chanHist[chan][0].second->FindBin(pid, nueenergy);
3228  int sigma = -1;
3229  int negIdx = 0;
3230  if (sigma < chanHist[chan].front().first) {
3231  negIdx = 0;
3232  }
3233  else if (sigma > chanHist[chan].back().first) {
3234  negIdx = chanHist[chan].size()-1;
3235 
3236  }
3237  else{
3238  for (int i = 0; i < (int)chanHist[chan].size()-1; i++){
3239  if (sigma >= chanHist[chan][i].first){
3240  negIdx = i;
3241  break;
3242  }
3243  }
3244  }
3245  // std::cout << chan << "\t" << pid << "\t" << nueenergy << "\t" << bin << "\t" << sigma << "\t" << negIdx
3246  // << "\t" << fHists[chan].front().first << "\t" << fHists[chan].back().first << "\n";
3247 
3248  sigma=+1;
3249  int posIdx=0;
3250  if (sigma < chanHist[chan].front().first) {
3251  // std::cout << "a) posIdx=0 for sigma= " << sigma << "\t" << fHists[chan].front().first << "\n";
3252  posIdx = 0;
3253  }
3254  else if (sigma > chanHist[chan].back().first) {
3255  // std::cout << "b) posIdx=" << posIdx << " for sigma= " << sigma << "\t" << fHists[chan].back().first << "\t" << fHists[chan].size() << "\n";
3256  posIdx = chanHist[chan].size()-1; //RJN change... Why was this -2??
3257  }
3258  else{
3259  for (int i = (int)chanHist[chan].size()-1; i>=0; i--){
3260  if (sigma >= chanHist[chan][i].first){
3261  posIdx = i;
3262  // std::cout << "c) posIdx=" << posIdx << " for sigma= " << sigma << "\t" << fHists[chan][i].first << "\t" << i << "\t" << fHists[chan].size() << "\n";
3263 
3264  break;
3265  }
3266  }
3267  }
3268 
3269  int negSigma = chanHist[chan][negIdx].first;
3270  int posSigma = chanHist[chan][posIdx].first;
3271 
3272  float negValue = chanHist[chan][negIdx].second->GetBinContent(bin);
3273  float posValue = chanHist[chan][posIdx].second->GetBinContent(bin);
3274  // std::cout << fHists[chan][negIdx].second->GetName() << "\t" << fHists[chan][posIdx].second->GetName() << "\n\t"
3275  // << chan << "\t" << pid << "\t" << nueenergy << "\t" << bin
3276  // << "\t" << negValue << "\t" << posValue << "\t" << negSigma << "\t" << posSigma << "\n";
3277 
3278  if(negSigma == -1 && posSigma == 1) {
3279  return { {-1.0, 1. + negValue}, {0.0, 1.0}, {1.0, 1. + posValue} };
3280  }
3281  else if(negSigma == 0 && posSigma == 1) {
3282  return { {-1.0, 1. - posValue}, {0.0, 1.0}, {1.0, 1. + posValue} };
3283  }
3284  return { {-1.0, 1.}, {0.0, 1.0}, {1.0 , 1.}};
3285  }
3286 
3287  //----------------------------------------------------------------------------
3289  {
3290  // do nothing if in the near detector or the selection is not nue
3294 
3295  auto wgtBySigma = this->CalcNueSystWeights(systName);
3296 
3297 // return this->CalcInterpolatedWgt(fCurrentPoint.ParameterValue(std::make_pair(systName, fCurrentMD.detector)),
3298 // wgtBySigma);
3300  wgtBySigma);
3301  }
3302 
3303  //----------------------------------------------------------------------------
3304  // We aren't currently using this one.
3306  {
3307  return 1.;
3308  }
3309 
3310  //----------------------------------------------------------------------------
3312  {
3313  return this->NueHistSystWeight("NueBirksCSystCVNWgt");
3314  }
3315 
3316  //----------------------------------------------------------------------------
3318  {
3319  return this->NueHistSystWeight("NueShiftXSystCVNWgt");
3320  }
3321 
3322  //----------------------------------------------------------------------------
3324  {
3325  return this->NueHistSystWeight("NueShiftYSystCVNWgt");
3326  }
3327 
3328  //----------------------------------------------------------------------------
3330  {
3331  return this->NueHistSystWeight("NueXYAbsSystCVNWgt");
3332  }
3333 
3334  //----------------------------------------------------------------------------
3336  {
3337  return this->NueHistSystWeight("NueXYRelSystCVNWgt");
3338  }
3339 
3340  //----------------------------------------------------------------------------
3342  {
3343  // Someone already called us
3344  if(!fSABeamSyst.fBeamSystHists.empty()) return;
3345 
3347 
3348  // TODO: make this a tunable parameter
3349  std::string shortName("TransportPlusNA49_NOvA");
3350 
3351  // constructor decides if initialized value is a path or an environment variable
3352  cet::search_path sp("FW_SEARCH_PATH");
3353  if ( !sp.find_file(std::string("CAFAna/data/beam/beamSyst_histos.root"), fileName) ) {
3354  throw cet::exception("BeamSystFileError")
3355  << "cannot find Beam Systematics file CAFAna/data/beam/beamSyst_histos.root "
3356  << fileName
3357  << " bail ungracefully\n"
3358  << "(RJN Guess: Did you run setup for this release?)\n"
3359  << __FILE__ << ":" << __LINE__;
3360  }
3361 
3362  fSABeamSyst.setFileAndHistNames(fileName, shortName);
3364  }
3365 
3366  //----------------------------------------------------------------------------
3367  // deprecated since this is included in the combination XSec/PPFX CV weight
3369  {
3370 
3372  fCurrentMD.BeamType(),
3373  fCurrentEvent->mcVals->val_at(fnex::kTrueE), //True Energy
3374  fCurrentEvent->mcVals->val_at(fnex::kTruePDGOrig)); //Flavour
3375 
3377  wgtBySigma);
3378  }
3379 
3380  //----------------------------------------------------------------------------
3382  {
3383  // it doesn't make sense to ask for a CC weight if the event is an NC
3384  // or cosmic muon
3387 
3388  return this->GENIESystWeight("MaCCQE_shapeWgt");
3389  }
3390 
3391  //----------------------------------------------------------------------------
3393  {
3394  // it doesn't make sense to ask for a CC weight if the event is an NC
3395  // or cosmic muon
3398 
3399  return this->GENIESystWeight("NormCCQEWgt");
3400  }
3401 
3402  //----------------------------------------------------------------------------
3404  {
3405  // it doesn't make sense to ask for a NC weight if the event isn't NC
3406  if(fCurrentMD.interactionType != fnex::kNC) return 1.;
3407 
3408  return this->GENIESystWeight("MvNCRESWgt");
3409  }
3410 
3411  //----------------------------------------------------------------------------
3413  {
3414  // it doesn't make sense to ask for a NC weight if the event isn't NC
3415  if(fCurrentMD.interactionType != fnex::kNC) return 1.;
3416 
3417  return this->GENIESystWeight("MaNCELWgt");
3418  }
3419 
3420  //----------------------------------------------------------------------------
3422  {
3423  return this->GENIESystWeight("MFP_NWgt");
3424  }
3425 
3426  //----------------------------------------------------------------------------
3428  {
3429  // it doesn't make sense to ask for a CC weight if the event is an NC
3430  // or cosmic muon
3433 
3434  return this->GENIESystWeight("CCQEPauliSupViaKFWgt");
3435  }
3436 
3437  //----------------------------------------------------------------------------
3439  {
3440  return this->GENIESystWeight("OtherGENIEWgt");
3441  }
3442 
3443  //----------------------------------------------------------------------------
3444  // For CAFAna version of this systematic, see CAFAna/Systs/NueAcceptSysts.cxx
3446  {
3447  //Only apply weight to FD, nue selected events
3451 
3452  //If this is the FHC systematic, only apply to FHC, and vice versa
3453  if((isFHC && fCurrentMD.BeamType() != fnex::kFHC) ||
3454  (!isFHC && fCurrentMD.BeamType() == fnex::kFHC))
3455  return 1;
3456 
3457  // LOG_VERBATIM("ShifterAndWeighter")
3458  // << "NueAcceptSignalKin2018Weight time!";
3459 
3460  //TODO: Check that it contains a true neutrino (CAF version is if(sr->mc.nnu == 0) return; )
3461 
3462  //Check flavours and that it's CC
3463  if(abs(fCurrentEvent->mcVals->val_at(fnex::kTruePDG)) != 12 ||
3464  abs(fCurrentEvent->mcVals->val_at(fnex::kTruePDGOrig)) != 14 ||
3466  return 1;
3467 
3468  // Grab the relevant histogram (if you haven't already): code stolen from LoadSABeamSystHists function, I hope that's still a valid way to get histograms
3469  if((isFHC && !fNueAcceptSignalKin2018FHCWgtHist) ||
3470  (!isFHC && !fNueAcceptSignalKin2018RHCWgtHist))
3471  {
3473  std::string shortName("Cos_FD_KinematicsCorrection_FHC");
3474  if(!isFHC) shortName = "Cos_FD_KinematicsCorrection_RHC";
3475 
3476  // constructor decides if initialized value is a path or an environment variable
3477  cet::search_path sp("FW_SEARCH_PATH");
3478  if ( isFHC && !sp.find_file(std::string("CAFAna/nue/Ana2018/AcceptSysts/FD_KinematicsCorrection_FHC.root"), fileName) ) {
3479  throw cet::exception("NueExtrapSystFileError")
3480  << "cannot find Nue Extrap Systematics file CAFAna/nue/Ana2018/AcceptSysts/FD_KinematicsCorrection_FHC.root "
3481  << fileName
3482  << " bail ungracefully\n"
3483  << "(RJN Guess: Did you run setup for this release?)\n"
3484  << __FILE__ << ":" << __LINE__;
3485  }
3486  else if (!isFHC && !sp.find_file(std::string("CAFAna/nue/Ana2018/AcceptSysts/FD_KinematicsCorrection_RHC.root"), fileName) ) {
3487  throw cet::exception("NueExtrapSystFileError")
3488  << "cannot find Nue Extrap Systematics file CAFAna/nue/Ana2018/AcceptSysts/FD_KinematicsCorrection_RHC.root "
3489  << fileName
3490  << " bail ungracefully\n"
3491  << "(RJN Guess: Did you run setup for this release?)\n"
3492  << __FILE__ << ":" << __LINE__;
3493  }
3494 
3495 
3496  TFile fin(fileName.c_str(), "read");
3497  if(fin.IsZombie()) {
3498  LOG_WARNING("ShifterAndWeighter")
3499  << "Warning: couldn't open nue extrapolation systematic file "
3500  << fileName;
3501  }
3502 
3503  if( isFHC) {
3504  fNueAcceptSignalKin2018FHCWgtHist = (TH1*)fin.Get(shortName.c_str())->Clone();
3505  fNueAcceptSignalKin2018FHCWgtHist->SetDirectory(0);
3506  } else {
3507  fNueAcceptSignalKin2018RHCWgtHist = (TH1*)fin.Get(shortName.c_str())->Clone();
3508  fNueAcceptSignalKin2018RHCWgtHist->SetDirectory(0);
3509  }
3510 
3511  fin.Close();
3512  }
3513 
3514  //find selection to find which bin to get from the histogram
3515  double selBin;
3517  selBin = 2.5;
3519  selBin = 0.5;
3520  // else if (fCurrentMD.selectionType == fnex::kNuESelectionMidPID)
3521  // selBin = 1.5;
3523  selBin = 1.5;
3524  else
3525  {
3526  selBin = -5;
3527  LOG_DEBUG("ShifterAndWeighter")
3528  << "Could not find selection type: " << fCurrentMD.selectionType;
3529  }
3530  int nuEBin = fCurrentEvent->dataVals->val_at(fnex::kNu_RecoE, fCurrentMD) / 0.5;
3531  int anaBin = 9*selBin + nuEBin - 3;//bins always seem to be 3 off. I don't know why, but I'm just adding a -3 to fix it
3532  LOG_DEBUG("ShifterAndWeighter")
3533  << "selBin = " << selBin
3534  << ", nuEBin = " << nuEBin
3535  << ", anaBin = " << anaBin;
3536 
3537  double FDSignalWeight = 1;
3538  if(isFHC && fNueAcceptSignalKin2018FHCWgtHist->GetBinContent(anaBin) != 0)
3539  FDSignalWeight = fNueAcceptSignalKin2018FHCWgtHist->GetBinContent(anaBin);
3540  else if(!isFHC && fNueAcceptSignalKin2018RHCWgtHist->GetBinContent(anaBin) != 0)
3541  FDSignalWeight = fNueAcceptSignalKin2018RHCWgtHist->GetBinContent(anaBin);
3542 
3543  //Find number of sigma to shift by
3544  double parVal = 0;
3545  if(isFHC)
3546  parVal = fCurrentPoint.ParameterValue(std::make_pair("NueAcceptSignalKin2018FHCWgt", fCurrentMD.detector));
3547  else
3548  parVal = fCurrentPoint.ParameterValue(std::make_pair("NueAcceptSignalKin2018RHCWgt", fCurrentMD.detector));
3549 
3550  LOG_DEBUG("ShifterAndWeighter")
3551  << "return weight " << 1+(FDSignalWeight-1)*parVal;
3552  return 1+(FDSignalWeight-1)*parVal;
3553 
3554  }
3555 
3556  //----------------------------------------------------------------------------
3557  // For CAFAna version of this systematic, see CAFAna/Systs/NueAcceptSysts.cxx
3559  {
3560 
3561  //Only apply weight to FD, nue selected events
3565 
3566  //If this is the FHC systematic, only apply to FHC, and vice versa
3567  if((isFHC && fCurrentMD.BeamType() != fnex::kFHC) ||
3568  (!isFHC && fCurrentMD.BeamType() == fnex::kFHC))
3569  return 1;
3570 
3571  double core = .015;
3572  double peri = .012;
3573  if(!isFHC) {
3574  core = 0.041;
3575  peri = 0.023;
3576  }
3577 
3578  //TODO: Check that it contains a true neutrino (CAF version is if(sr->mc.nnu == 0) return; )
3579 
3580  //Check if neutrino, but not numu->nue signal, otherwise left unaltered
3581  if((abs(fCurrentEvent->mcVals->val_at(fnex::kTruePDG)) == 12 &&
3582  abs(fCurrentEvent->mcVals->val_at(fnex::kTruePDGOrig)) == 14 &&
3584  (abs(fCurrentEvent->mcVals->val_at(fnex::kTruePDG)) == 16 &&
3586  return 1;
3587 
3588  double FDBkgCorr = core;
3590  FDBkgCorr = peri;
3591 
3592  //Find number of sigma to shift by
3593  double parVal = 0;
3594  if(isFHC)
3595  parVal = fCurrentPoint.ParameterValue(std::make_pair("NueAcceptBkg2018FHCWgt", fCurrentMD.detector));
3596  else
3597  parVal = fCurrentPoint.ParameterValue(std::make_pair("NueAcceptBkg2018RHCWgt", fCurrentMD.detector));
3598 
3599  return 1 + FDBkgCorr*parVal;
3600 
3601  }
3602 
3603  //----------------------------------------------------------------------------
3605  {
3606  return this->GENIESystWeight("FormZoneWgt");
3607  }
3608 
3609  //----------------------------------------------------------------------------
3611  {
3612  return this->GENIESystWeight("FrInel_piWgt");
3613  }
3614 
3615  //----------------------------------------------------------------------------
3617  {
3618  return this->GENIESystWeight("FrCEx_NWgt");
3619  }
3620 
3621  //----------------------------------------------------------------------------
3623  {
3624  return this->GENIESystWeight("FrElas_NWgt");
3625  }
3626 
3627  //----------------------------------------------------------------------------
3629  {
3630  return this->GENIESystWeight("FrAbs_NWgt");
3631  }
3632 
3633 
3634 } // fnex
void SetCurrentVals(fnex::InputPoint const &curPoint)
double FluxPCAWeight(size_t index)
T max(const caf::Proxy< T > &a, T b)
double CalibNormWeight(std::string const &wgtType, novadaq::cnv::DetId const &det)
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
TString fin
Definition: Style.C:24
double NormBothWeight(std::string const &wgtType, fnex::SelectionType_t const &sel)
const XML_Char * name
Definition: expat.h:151
std::vector< WeighterInfo > fWeightersToUseND
static const unsigned char kNu_RecoE
Definition: VarVals.h:36
std::map< float, float > CalcCalibSystWeights(std::string const &systName, fnex::MetaData currentMD, double trueE)
static const unsigned char kNuE_CVN
Definition: VarVals.h:32
void setFileAndHistNames(std::string filename, std::string shortName)
enum fnex::weight_type WeightType_t
fileName
Definition: plotROC.py:78
bool IsNuESelected() const
Definition: Structs.cxx:265
void SetTh13(const T &th13) override
double const ParameterValue(ParameterDet const &parDet) const
Definition: InputPoint.cxx:296
bool IsNuMuSelected() const
Definition: Structs.cxx:256
set< int >::iterator it
static const unsigned char kEmpiricalMEC_Weight
Definition: VarVals.h:60
void SetL(double L) override
std::set< SigExtrapWgt > fSigExtrapWeights
The histgram with weights as a function of true energy.
double NueAcceptBkg2018Weight(bool isFHC)
static const unsigned char kOtherDIS_Weight
Definition: VarVals.h:63
const Var weight
static const unsigned char kDISvnCC1pi_Weight
Definition: VarVals.h:59
BeamSyst fPPFXSmoothCVWgt
This is the PPFX Central Value.
std::map< float, float > CalcGeniePCASystWeights(fnex::MetaData currentMD, int mode, double energy)
Base container for the MC related Vars that constitute an event.
Definition: VarVals.h:930
static const unsigned char kEmpiricalMECtoGENIERES_Weight
Definition: VarVals.h:62
void InitShiftsAndWeightsToUse(fnex::InputPoint const &curPoint, std::string CalibFilename="FNEX/data/calibration/CalibrationSystematics2018.root")
TFileDirectory mkdir(std::string const &dir, std::string const &descr="")
void LoadCalibSystHists(std::string CalibFilename)
NueSystsSAW::OscChannel GetOscChannel()
static const unsigned char kRecoQ2
Definition: VarVals.h:38
std::vector< std::array< float, 4 > > CalcInterpFitCoeffs(std::map< float, float > const &wgtsBySigma)
static const unsigned char kXSecCVPPFX_Weight
Definition: VarVals.h:56
string fFileName
Definition: nue_pid_effs.C:43
osc::OscCalcPMNSOpt fOscCalc
the oscillation calculator
static const unsigned char kTrueHitNuc
Definition: VarVals.h:52
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
void abs(TH1 *hist)
std::map< fnex::MetaData, fnex::EventList > EventListMap
Definition: Event.h:186
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
static const unsigned char kFake_Weight
Definition: VarVals.h:39
Create a list of fnex::Events to be used in fits.
void LoadCalibrationSystHists(std::string CalibFilename)
T cube(T x)
More efficient cube function than pow(x,3)
Definition: MathUtil.h:26
double GENIESystWeight(std::string const &systName)
std::map< unsigned char, std::set< ShifterInfo > > fShiftersToUseND
static ShifterAndWeighter * Instance()
std::string find_file(std::string const &filename) const
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
Loaders::FluxType flux
float abs(float number)
Definition: d0nt_math.hpp:39
static const unsigned char kRPARES_Weight
Definition: VarVals.h:58
double TotalShift(std::map< unsigned char, std::set< ShifterInfo > > const &shiftersToUse, unsigned char const &param)
std::vector< BeamSyst > fVecFluxPrincipals
This is each of the Flux PCA weights.
static const unsigned char kTruePDGOrig
Definition: VarVals.h:46
static const unsigned char kTrueCCNC
Definition: VarVals.h:44
const XML_Char * s
Definition: expat.h:262
double GeniePCASystWeight(size_t index)
Double_t scale
Definition: plot.C:25
std::string const ToString() const
Definition: Structs.cxx:114
Far Detector at Ash River, MN.
std::map< ToFCounter, std::vector< unsigned int > > channels
double MECEnuShapeWeight(bool anti)
double NueHistSystWeight(std::string const &systName)
bool fSaveRatios
save the ND weights used to make the signal prediction during extrapolation
enum fnex::sel_type SelectionType_t
Float_t E
Definition: plot.C:20
fnex::MetaData fCurrentMD
the metadata for the current event
std::vector< std::vector< std::vector< std::vector< TH1D * > > > > fBeamSystHists
ParameterMap const ParameterSpaceLocation(novadaq::cnv::DetId const &det) const
Definition: InputPoint.cxx:330
static const unsigned char kLep_RecoE_MCFrac
Definition: VarVals.h:37
double Weight(fnex::EventPtr const &curEvent, fnex::MetaData const &md, fnex::WeightType_t const wgtType=kAllWgt)
std::vector< TH1D * > fHornCurrentHists
correl_xv GetXaxis() -> SetDecimals()
double energy
Definition: plottest35.C:25
std::map< std::string, std::vector< std::vector< std::pair< int, TH1D * > > > > fNueSystHists
Near Detector in the NuMI cavern.
void Load2018CalibSystHists(std::string CalibFilename)
Float_t d
Definition: plot.C:236
const int kNumQuantiles
Definition: getTimePeak.C:42
T P(int flavBefore, int flavAfter, double E) override
E in GeV; flavors as PDG codes (so, neg==>antinu)
std::vector< GeniePCASyst > fVecGeniePCASysts
This is each of the Genie PCA weights.
void FillShifterInfo(std::string const &name, unsigned char const &param, novadaq::cnv::DetId const &det)
const double j
Definition: BetheBloch.cxx:29
const Binning bins
static const unsigned char kTrueIntType
Definition: VarVals.h:45
void SetTh23(const T &th23) override
novadaq::cnv::DetId detector
Definition: Structs.h:50
double ShiftAmount(unsigned char const &param)
std::map< float, float > CalcBeamSystWeights(novadaq::cnv::DetId const &curDet, fnex::BeamType_t const &beamType, double energy, int pdgOrig)
const int kNumSels
Definition: vars.h:43
float bin[41]
Definition: plottest35.C:14
static ShifterAndWeighter * gShifterAndWeighter
const ana::Var wgt
z
Definition: test.py:28
fnex::FileType_t fileType
Definition: Structs.h:51
void SetDmsq21(const T &dmsq21) override
static bool isFHC
double sigma(TH1F *hist, double percentile)
#define LOG_WARNING(category)
double RelNormWeight(std::string const &wgtType, fnex::SelectionType_t const &sel)
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
void SetDmsq32(const T &dmsq32) override
fnex::SelectionType_t selectionType
Definition: Structs.h:52
double AbsNormWeight(std::string const &wgtType, fnex::SelectionType_t const &sel)
fnex::InputPoint fCurrentPoint
the current point in parameter space
std::map< unsigned char, std::set< ShifterInfo > > fShiftersToUseFD
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
charged current coherent pion
Definition: MCNeutrino.h:139
void SetCurrentEvent(fnex::EventPtr const &curEvent, fnex::MetaData const &md)
double CalcLinearInterpWgt(double sigma, std::map< float, float > const &wgtsBySigma)
static const unsigned char kEmpiricalMECtoGENIEQE_Weight
Definition: VarVals.h:61
double CalcInterpolatedWgt(double sigma, std::map< float, float > const &wgtsBySigma)
fnex::InteractionType_t interactionType
Definition: Structs.h:53
TDirectory * dir
Definition: macro.C:5
static const unsigned char kRPACCQE_Weight
Definition: VarVals.h:57
double HornNormWeight(std::string const &wgtType, fnex::BeamType_t const &beam)
void CalcNearTrueEnergyRatioWeights(const fnex::EventListMap *pEventListMap)
std::shared_ptr< Event > EventPtr
Definition: Event.h:50
std::map< std::string, std::string > systNames
Definition: PlotUnfolding.C:32
double BirksNormWeight(std::string const &wgtType, novadaq::cnv::DetId const &det)
std::vector< WeighterInfo > fWeightersToUseFD
void SetdCP(const T &dCP) override
fvar< T > floor(const fvar< T > &x)
Definition: floor.hpp:11
void SetTh12(const T &th12) override
static const unsigned char kTrueE
Definition: VarVals.h:42
fnex::Event * fCurrentEvent
the current event to weight/shift
string shortName
THUMBNAIL BLOCK: We need to make a thumbnail for each.
enum fnex::beam_type BeamType_t
void LoadGeniePCASystHists(std::string shortName)
static const unsigned char kHad_RecoE
Definition: VarVals.h:34
double CalibSystWeight(std::string const &wgtType, novadaq::cnv::DetId const &det)
double BeamSystWeight(std::string const &wgtType, novadaq::cnv::DetId const &det)
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
std::string KeyToVarName(unsigned char const &key)
Definition: VarVals.h:445
static const unsigned char kTrueIntMode
Definition: VarVals.h:53
std::unique_ptr< MCVarVals > mcVals
Definition: Event.h:36
#define LOG_VERBATIM(category)
fnex::BeamType_t const BeamType() const
Definition: Structs.cxx:165
double NueAcceptSignalKin2018Weight(bool isFHC)
const std::string cSelectionType_Strings[11]
Definition: Constants.h:101
double MECInitStateNPFracWeight(bool anti)
double const Sigma(std::string const &name) const
static const unsigned char kLep_RecoE
Definition: VarVals.h:35
Base container for the MC related Vars that constitute an event.
Definition: VarVals.h:874
std::map< float, float > CalcNueSystWeights(std::string const &systName)
def sign(x)
Definition: canMan.py:197
void SetRho(double rho) override
std::unique_ptr< DataVarVals > dataVals
Definition: Event.h:35
ParameterMap const SystematicPulls(novadaq::cnv::DetId const &det) const
Definition: InputPoint.cxx:364
static SystematicSigmas * Instance()
static const unsigned char kTruePDG
Definition: VarVals.h:43