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 CovarianceMatrixFit fitting
5 //
6 
7 #include "TH2.h"
8 #include "TFile.h"
9 #include "TKey.h"
10 #include <algorithm>
11 
15 
20 
21 #include "OscLib/OscCalcPMNS_NSI.h"
22 #include "OscLib/OscCalcSterile.h"
23 #include "Utilities/func/MathUtil.h"
24 
25 
26 namespace cmf {
27 
29  std::string const& shortName)
30  : fFileName(fileName)
31  , fShortName(shortName)
32  {
33  }
34 
35  //-----------------------------------------------------------------------------
37  std::string const& shortName)
38  {
41  }
42 
43  //-----------------------------------------------------------------------------
45  {
46  // Someone already called us
47  if(!fBeamSystHists.empty()) return;
48 
49  //open file
50  TFile fin (fFileName.c_str(),"read");
51  if(fin.IsZombie()){
52  LOG_WARNING("ShifterAndWeighter")
53  << "Warning: couldn't open beam systematic file "
54  << fFileName;
55  return;
56  }
57 
58  //Check that desired systematic is in file
59  // Loop through the detector folder, over all histograms
60  TIter iterHist(gDirectory->GetListOfKeys());
61  TKey* keyHist;
62  int foundHisto = 0;
63  bool isSeparatedByFlavor = false;
64 
66 
67  while((keyHist = (TKey*)iterHist())) {
68  TString histName = keyHist->GetName(); // Get a histogram name
69 
70  LOG_DEBUG("ShifterAndWeighter")
71  << "histName is "
72  << histName;
73 
74  if(histName.Contains(fShortName)){ //is it the desired syst
75  ++foundHisto;
76 
77  key = (histName.Contains("FHC") ) ? "FHC" : "RHC";
78  key += (histName.Contains("ND") ) ? "ND" : "FD";
79 
80  // order matters in these if/else if statements because
81  // anumu contains numu, anue contains nue but not the
82  // other way around
83  if (histName.Contains("anumu")) key += "anumu";
84  else if(histName.Contains("numu" )) key += "numu";
85  else if(histName.Contains("anue" )) key += "anue";
86  else if(histName.Contains("nue" )) key += "nue";
87 
88  key += (histName.Contains("minus")) ? "minus" : "plus";
89 
90  //store relevant histograms
91  fBeamSystHists[key] = (TH1D*) fin.Get(histName)->Clone();
92  // disassociate it from the file it came from
93  fBeamSystHists[key]->SetDirectory(nullptr);
94  }
95  if(histName.Contains("numu") || histName.Contains("nue")){
96  isSeparatedByFlavor = true;
97  }
98  }
99 
100  if(!foundHisto){
101  throw cet::exception("BeamSystFileError")
102  << "Beam systematic "
103  << fShortName
104  << " is not in this file "
105  << fFileName;
106  }
107 
108  if(!isSeparatedByFlavor){
109  throw cet::exception("BeamSystFileError")
110  << "Beam systematic "
111  << fShortName
112  << " doesn't have required flavor information";
113  }
114  fin.Close();
115  }
116 
117  //----------------------------------------------------------------------------
118  std::map<float, float> BeamSyst::CalcBeamSystWeights(cmf::DetType_t const& curDet,
119  cmf::BeamType_t const& beamType,
120  double energy,
121  int pdgOrig)
122  {
123  // Need different histogram per horn current
124  std::string key = (beamType == cmf::kFHC) ? "FHC" : "RHC";
125 
126  // Need different histogram per detector
127  key += (curDet == cmf::kNEARDET) ? "ND" : "FD";
128 
129  //Find the neutrino flavor
130  switch(pdgOrig){
131  case 14: key += "numu"; break;
132  case 12: key += "nue"; break;
133  case -14: key += "anumu"; break;
134  case -12: key += "anue"; break;
135  case 13:
136  //Cosmic muon
137  return { {-1.0, 1.0}, {0.0, 1.0}, {1.0, 1.0} };
138  default:
139  LOG_WARNING("BeamSystFileError")
140  << "Wrong nu flavor; "
141  << pdgOrig
142  << " ignore";
143 
144  return { {-1.0, 1.0}, {0.0, 1.0}, {1.0, 1.0}};
145  }
146 
147  //Find the right histogram for the detector and sign
148  TH1D* hMinus = fBeamSystHists[key + "minus"];
149  TH1D* hPlus = fBeamSystHists[key + "plus" ];
150 
151  double scaleMinus = 1.;
152 
153  if(hMinus){
154  if(energy > hMinus->GetXaxis()->GetXmin() &&
155  energy < hMinus->GetXaxis()->GetXmax())
156  scaleMinus = hMinus->Interpolate(energy);
157  }
158 
159  double scalePlus = 1.;
160 
161  if(hPlus){
162  if(energy > hPlus->GetXaxis()->GetXmin() &&
163  energy < hPlus->GetXaxis()->GetXmax())
164  scalePlus = hPlus->Interpolate(energy);
165  }
166 
167  return { {-1.0, scaleMinus}, {0.0, 1.0}, {1.0, scalePlus} };
168  } // std::map<double, double> BeamSystWeightFn::CalcWeights(const cmf::Event&)
169 
170  //-------------------------------------------------------------------------
171  void CalibSyst::LoadCalibSystRatios(std::string const& CalibFilename)
172  {
173  // Someone already called us
174  if(!fCalibSystRatios.empty()) return;
175 
176  LOG_DEBUG("ShifterAndWeighter")
177  << "LoadCalibSystRatios";
178 
179  std::string fileCalib;
180 
181  LOG_DEBUG("ShifterAndWeighter")
182  << "Looking for calib syst file "
183  << CalibFilename;
184 
185  cet::search_path sp("FW_SEARCH_PATH");
186  if(!sp.find_file(CalibFilename, fileCalib)){
187  throw cet::exception("CalibSystematicFileError")
188  << "Cannot find calib systematics file "
189  << CalibFilename
190  << " bail ungracefully\n\n"
191  << __FILE__ << ":" << __LINE__;
192  }
193 
194  LOG_DEBUG("ShifterAndWeighter")
195  << "Found "
196  << CalibFilename
197  << " as "
198  << fileCalib;
199 
200  TFile fin(fileCalib.c_str(), "read");
201  if(fin.IsZombie()){
202  throw cet::exception("Weighters")
203  << "Warning: couldn't open "
204  << fileCalib
205  << ". Crashing";
206  }
207 
208  LOG_DEBUG("ShifterAndWeighter")
209  << "Opened file, apparently";
210 
211  // make some string collections with the combinations we need
212  // for the names of the histograms
213  std::vector<std::string> detectors({"Near",
214  "Far"});
215 
216  std::vector<std::string> selections({"NuESel_LowPID",
217  "NuESel_HighPID",
218  "NuESel_Peripheral",
219  "NuMuSelQ1",
220  "NuMuSelQ2",
221  "NuMuSelQ3",
222  "NuMuSelQ4",
223  "NCSel"});
224 
225  std::vector<std::string> systematics;
226 
227 #ifdef PRODUCTION5
228  LOG_VERBATIM("p5") << __LINE__;
229  systematics.push_back("calibration");
230 #else
231  LOG_VERBATIM("p5") << __LINE__;
232  systematics.push_back("lightmodel-lightup-calibdown");
233  systematics.push_back("lightmodel-lightdown-calibup");
234  systematics.push_back("ckv-proton-shift-down");
235  systematics.push_back("calib-shift-xyview-pos-offset");
236  systematics.push_back("calib-shift-xyview-neg-offset");
237  systematics.push_back("calib-shift-func");
238 #endif
239  std::vector<std::string> beamTypes({"FHC",
240  "RHC"});
241 
242  std::string ratioName;
243  TH1D* ratioHist;
244  for(auto const& det : detectors){
245  for(auto const& sel : selections){
246 
247  if(det == "Near" &&
248  sel == "NuESel_Peripheral") continue;
249 
250  for(auto const& beam : beamTypes){
251  for(auto const& syst : systematics){
252 
253 #ifdef PRODUCTION5
254  ratioName = syst;
255  ratioName += "Total";
256  ratioName += beam;
257  ratioName += det;
258  ratioName += sel;
259  ratioName += "Plus1Sigma";
260 
261  ratioHist = dynamic_cast<TH1D*>(fin.Get((ratioName).c_str()));
262 #else
263  ratioName = det;
264  ratioName += beam;
265  ratioName += sel;
266  ratioName += syst;
267 
268  // get the ratio from the TFile
269  ratioHist = dynamic_cast<TH1D*>(fin.Get((ratioName + "Ratio").c_str()));
270 #endif
271 
272  if(ratioHist == nullptr){
273  LOG_WARNING("CalibSyst")
274  << "failed to find "
275  << ratioName
276  << " in "
277  << fileCalib
278  << " try the next systematic ";
279 
280  continue;
281  } // end if the requested ratio is not there
282 
283  // loop over the histogram and put the ratios wrt nominal in
284  // the map
285  for(int i = 1; i < ratioHist->GetNbinsX() + 1; ++i){
286  fCalibSystRatios[ratioName].push_back(ratioHist->GetBinContent(i));
287 
288  // if(sel == selections.front() ||
289  // sel == selections.back() )
290  // LOG_VERBATIM("CalibSystRatios")
291  // << ratioName
292  // << " "
293  // << i
294  // << " "
295  // << ratioHist->GetBinCenter(i)
296  // << " "
297  // << fCalibSystRatios[ratioName].back();
298  }
299 
300  } // end loop over systematic names
301  } // end loop over beam types
302  } // end loop over selections
303  } // end loop over detectors
304 
305  LOG_DEBUG("CalibSyst") << "Calib histograms loaded!";
306  }
307 
308  //----------------------------------------------------------------------
309  float CalibSyst::CalibSystRatio(std::string const& systName,
310  cmf::MetaData const& currentMD,
311  double energy)
312  {
313  // use the metadata to figure out which ratios we need
314  std::string ratioName;
315 
316 #ifdef PRODUCTION5
317  ratioName = systName;
318  ratioName += "Total";
319  ratioName += cmf::cBeamType_Strings[currentMD.BeamType()];
320  ratioName += currentMD.DetectorString();
321  ratioName += cmf::cSelectionType_Strings[currentMD.selectionType];
322  ratioName += "Plus1Sigma";
323 #else
324  ratioName += currentMD.DetectorString();
325  ratioName += cmf::cBeamType_Strings[currentMD.BeamType()];
326  ratioName += cmf::cSelectionType_Strings[currentMD.selectionType];
327  ratioName += systName;
328 #endif
329 
330  // find the energy bin for this energy
331  int bin = cmf::CovarianceBinUtility::Instance()->EnergyToBin(energy, currentMD);
332 
333  if(fCalibSystRatios.count(ratioName) < 1){
334  LOG_WARNING("CalibSystRatio")
335  << "requested calibration ratio for "
336  << ratioName
337  << " not found, returning 1";
338 
339  return 1;
340  }
341 
342  return fCalibSystRatios[ratioName][bin];
343  }
344 
345  //-------------------------------------------------------------------------
347  {
348  // Someone already called us
349  if(!fGeniePCAHists.empty()) return;
350 
351  const std::string kGeniePCASystFile = "CAFAna/data/xs/genie_small_pc_shifts_fn_2018.root";
352 
353  std::string fileGeniePCA;
354 
355  cet::search_path sp("FW_SEARCH_PATH");
356  if ( !sp.find_file(kGeniePCASystFile, fileGeniePCA) ) {
357  throw cet::exception("GeniePCASystematicFileError")
358  << "Cannot find Genie PCA systematics file "
359  << fileGeniePCA
360  << " bail ungracefully\n\n"
361  << __FILE__ << ":" << __LINE__;
362  }
363 
364  TFile fin(fileGeniePCA.c_str(), "read");
365  if (fin.IsZombie()){
366  throw cet::exception("Weighters")
367  << "Warning: couldn't open "
368  << fileGeniePCA
369  << ". Crashing!";
370  }
371 
372  //Check that desired systematic is in file
373  // Loop through the detector folder, over all histograms
374  TIter iterDir(gDirectory->GetListOfKeys());
375  TKey* keyDir;
376  int foundHist = 0;
377 
379 
380  while((keyDir = (TKey*)iterDir())) {
381  TString histName = keyDir->GetName();
382  if(histName.Contains(shortName)) {
383  foundHist++;
384 
385  key = (histName.Contains("FHC") ) ? "FHC" : "RHC";
386  key += (histName.Contains("ND") ) ? "ND" : "FD";
387  key += (histName.Contains("Numu") ) ? "NumuSel" : "NueSel";
388 
389  if (histName.Contains("SigQE" )) key += "SigQE";
390  else if(histName.Contains("SigNonQE")) key += "SigNonQE";
391  else if(histName.Contains("Other" )) key += "BkgFlav";
392  else if(histName.Contains("NC" )) key += "NC";
393 
394  key += (histName.Contains("minus")) ? "minus" : "plus";
395 
396  // grab relevant histogram
397  TH1D *hist = (TH1D*)fin.Get(histName)->Clone();
398 
399  if(hist) {
400  // store relevant histograms
401  fGeniePCAHists[key] = hist;
402 
403  // disassociate from file
404  fGeniePCAHists[key]->SetDirectory(nullptr);
405  }
406  else {
407  throw cet::exception("Weighters")
408  << "Genie PCA histogram "
409  << histName
410  << " does not exist! Something has gone terribly wrong!";
411  } // end if statement checking hist exists
412  } // end hist name check
413  } // end while loop
414 
415  LOG_DEBUG("GeniePCASysts") << "Genie PCA histograms loaded!";
416  }
417 
418  //----------------------------------------------------------------------
419  std::map<float, float> GeniePCASyst::CalcGeniePCASystWeights(cmf::MetaData const& currentMD,
420  int mode,
421  double energy)
422  {
424 
425  key = (currentMD.BeamType() == cmf::kFHC) ? "FHC" : "RHC";
426  key += (currentMD.detector == cmf::kNEARDET) ? "ND" : "FD";
427  if ( currentMD.IsNuMuSelected() ) key += "NumuSel";
428  else if( currentMD.IsNuESelected() ) key += "NueSel";
429  else if( currentMD.IsNCSelected() ) key += "NCSel";
430 
431  // if it's any of the interactions we're not interested in then just return 1
432  cmf::InteractionType_t interaction = currentMD.interactionType;
433 
434  if (currentMD.IsNuMuSelected() &&
435  interaction == cmf::kNuMuCC &&
436  mode == simb::kQE) key += "SigQE";
437  else if(currentMD.IsNuESelected() &&
438  interaction == cmf::kNuECC &&
439  mode == simb::kQE) key += "SigQE";
440  else if(currentMD.IsNuMuSelected() &&
441  interaction == cmf::kNuMuCC &&
442  mode != simb::kQE) key += "SigNonQE";
443  else if(currentMD.IsNuESelected() &&
444  interaction == cmf::kNuECC &&
445  mode != simb::kQE) key += "SigNonQE";
446  else if(interaction == cmf::kNC) key += "NC";
447  else key += "BkgFlav";
448 
449  // for some reason background histograms don't exist for numu so ignore these
450  if( fGeniePCAHists.count(key + "minus") < 1 ||
451  fGeniePCAHists.count(key + "plus") < 1)
452  return { {-1.0, 1.0}, {0.0, 1.0}, {1.0, 1.0} };
453 
454  // find the right histogram for the detector and sign
455  TH1D* hMinus = fGeniePCAHists[key + "minus"];
456  TH1D* hPlus = fGeniePCAHists[key + "plus" ];
457 
458  double scaleMinus = 1.;
459 
460  if(hMinus){
461  if(energy > hMinus->GetXaxis()->GetXmin() &&
462  energy < hMinus->GetXaxis()->GetXmax())
463  scaleMinus = hMinus->Interpolate(energy);
464  }
465 
466  double scalePlus = 1.;
467 
468  if(hPlus){
469  if(energy > hPlus->GetXaxis()->GetXmin() &&
470  energy < hPlus->GetXaxis()->GetXmax())
471  scalePlus = hPlus->Interpolate(energy);
472  }
473 
474  return { {-1.0, scaleMinus}, {0.0, 1.0}, {1.0, scalePlus} };
475  }
476 
477  //=======================================================================//
478  //------- GENERIC/GENERAL PURPOSE CODE (INITIALISATION ETC.) ------------//
479  //=======================================================================//
480 
481 
483 
484  //----------------------------------------------------------------------------
486  {
487  if(!gShifterAndWeighter) gShifterAndWeighter = new ShifterAndWeighter();
488 
489  return gShifterAndWeighter;
490  }
491 
492  //----------------------------------------------------------------------------
494  : fOscCalc (nullptr)
495  , fCalcType(cmf::kPMNS)
496  {
497  // create graphs for the muon catcher positive and negative energy shifts,
498  // then use those to create functions that can be used during the evaluation
499  // of the shifts
500  std::vector<double> initialEnergy({0.00872008, 0.0640129, 0.119306, 0.174598, 0.229891,
501  0.285184, 0.340477, 0.616941, 0.672234, 0.727526,
502  0.782819, 0.838112, 0.893405, 0.948697, 1.00399,
503  1.05928, 1.11458, 1.16987, 1.36119, 1.3964,
504  1.43162, 1.46684, 1.50206, 1.53728, 1.57249,
505  1.60771, 1.64293, 1.67815, 1.71337, 1.74858,
506  1.7838, 1.81902, 1.85424, 1.88946, 1.92467,
507  1.95989, 1.99511, 2.03033, 2.06555, 2.10076,
508  2.20642, 2.24164, 2.27685, 2.31207, 2.34729,
509  2.38251, 2.41773, 2.45294, 2.48816, 2.52338,
510  2.5586, 2.59381, 2.62903, 2.66425, 2.69947,
511  2.73469, 2.7699, 2.80512, 2.84034, 2.87556,
512  2.91078, 2.94599, 2.98121, 3.01643, 3.05165,
513  3.08687, 3.12208, 3.1573, 3.19252, 3.22774,
514  3.26296, 3.29817});
515  std::vector<double> posShift({7.6579, 1.90696, 1.48663, 1.33252, 1.25254,
516  1.20358, 1.17052, 1.09411, 1.08636, 1.0798,
517  1.07416, 1.06927, 1.06498, 1.0612, 1.05783,
518  1.05481, 1.05209, 1.04963, 1.02717, 1.02648,
519  1.02583, 1.02521, 1.02462, 1.02405, 1.02352,
520  1.023, 1.02251, 1.02204, 1.02158, 1.02115,
521  1.02073, 1.02033, 1.01994, 1.01957, 1.01921,
522  1.01887, 1.01853, 1.01821, 1.0179, 1.0176,
523  1.01676, 1.0165, 1.01624, 1.01599, 1.01575,
524  1.01552, 1.01529, 1.01508, 1.01486, 1.01465,
525  1.01445, 1.01426, 1.01407, 1.01388, 1.0137,
526  1.01352, 1.01335, 1.01318, 1.01302, 1.01286,
527  1.0127, 1.01255, 1.0124, 1.01226, 1.01212,
528  1.01198, 1.01184, 1.01171, 1.01158, 1.01146,
529  1.01133, 1.01121});
530  std::vector<double> negShift({0, 0, 0.513372, 0.66748, 0.747457,
531  0.796421, 0.829482, 0.905895, 0.913635, 0.920199,
532  0.925835, 0.930728, 0.935016, 0.938803, 0.942173,
533  0.945192, 0.947911, 0.950373, 0.972833, 0.973519,
534  0.97417, 0.97479, 0.975381, 0.975945, 0.976484,
535  0.976999, 0.977492, 0.977964, 0.978417, 0.978852,
536  0.97927, 0.979671, 0.980057, 0.980429, 0.980787,
537  0.981132, 0.981465, 0.981787, 0.982097, 0.982397,
538  0.98324, 0.983504, 0.983759, 0.984006, 0.984246,
539  0.984479, 0.984705, 0.984925, 0.985138, 0.985346,
540  0.985547, 0.985743, 0.985934, 0.98612, 0.986301,
541  0.986478, 0.98665, 0.986817, 0.986981, 0.98714,
542  0.987296, 0.987448, 0.987596, 0.987741, 0.987882,
543  0.988021, 0.988156, 0.988288, 0.988417, 0.988543,
544  0.988667, 0.988788});
545 
546  // Make TGraphs for the positive and negative shifts
547  fNDMuonCatcherPlusGr = std::make_unique<TGraph>(initialEnergy.size(), &initialEnergy[0], &posShift[0]);
548  fNDMuonCatcherMinusGr = std::make_unique<TGraph>(initialEnergy.size(), &initialEnergy[0], &negShift[0]);
549 
550  // Use the TGraphs to make the stored TF1s
551  fNDMuonCatcherEShiftPlus = std::make_unique<TF1>("posMuonEnergyShift",[&](double*x, double *p){ return fNDMuonCatcherPlusGr ->Eval(x[0]); }, 0., 3.5, 1);
552  fNDMuonCatcherEShiftMinus = std::make_unique<TF1>("negMuonEnergyShift",[&](double*x, double *p){ return fNDMuonCatcherMinusGr->Eval(x[0]); }, 0., 3.5, 1);
553 
554  }
555 
556  //----------------------------------------------------------------------------
558  cmf::ParameterSpaceLoc const& detLoc,
559  std::string const& calibFileName,
560  std::vector<WeighterInfo> & weighters,
561  std::map<std::string, uint8_t> & shifters)
562  {
563  // loop over the systematic parameters from the input point and evaluate the
564  // name associated with each one, then add the appropriate function to the
565  // list.
566  weighters.clear();
567  shifters .clear();
568 
569  std::string systName;
570  WeighterInfo weighter;
571  for(auto const& itr : detLoc){
572  if(itr.second.IsOscPar()) continue;
573  systName = itr.first;
574 
575  if(systName.find("Wgt") != std::string::npos){
576  weighter.wgtName = systName;
577 
578  // check all the systematic parameter names and add the necessary functions
579  // for weighting if there is a function for that name.
580  if (systName == "AbsNormWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::AbsNormWeight, this, systName, cmf::kNuMuSelection);
581  else if(systName == "BirksNearSystWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::BirksNormWeight, this, systName, cmf::kNEARDET);
582  else if(systName == "BirksFarSystWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::BirksNormWeight, this, systName, cmf::kFARDET );
583  else if(systName == "CosmicMuNormWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::CosmicMuNormWeight, this);
584  else if(systName == "CosmicOverlayNormWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::CosmicOverlayNormWeight, this);
585  else if(systName == "CCQEPauliSupViaKFWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
586  else if(systName == "COHCCScale2018Wgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::COHScale2018Weight, this, systName, cmf::kNuMuCC);
587  else if(systName == "COHNCScale2018Wgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::COHScale2018Weight, this, systName, cmf::kNC);
588  else if(systName == "DISvnCC1piWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::DISvnCC1piWeight, this);
589  else if(systName == "FHCNormWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::HornNormWeight, this, systName, cmf::kFHC );
590  else if(systName == "FormZoneWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
591  else if(systName == "FrInel_piWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
592  else if(systName == "FrCEx_NWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
593  else if(systName == "FrElas_NWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
594  else if(systName == "FrAbs_NWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
595  else if(systName == "MaCCQEReducedWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::MaCCQEReducedWeight, this);
596  else if(systName == "MaCCQE_shapeWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
597  else if(systName == "MaCCRESWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
598  else if(systName == "MaNCRESWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
599  else if(systName == "MaNCELWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
600  else if(systName == "MECShapeNuWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::MECShapeWeight, this, false);
601  else if(systName == "MECShapeAntiNuWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::MECShapeWeight, this, true );
602  else if(systName == "MECInitNPFracNuWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::MECInitStateNPFracWeight, this, false);
603  else if(systName == "MECInitNPFracAntiNuWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::MECInitStateNPFracWeight, this, true );
604  else if(systName == "MECEnuShapeNuWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::MECEnuShapeWeight, this, false);
605  else if(systName == "MECEnuShapeAntiNuWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::MECEnuShapeWeight, this, true );
606  else if(systName == "MECNormWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::MECNormWeight, this);
607  else if(systName == "MFP_NWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
608  else if(systName == "MvCCRESWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
609  else if(systName == "MvNCRESWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
610  else if(systName == "NCNormWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::NCNormWeight, this);
611  else if(systName == "NoiseFarSystWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::NoiseFarNormWeight, this);
612  else if(systName == "NormCCQEWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
613  else if(systName == "NueAbsNormWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::AbsNormWeight, this, systName, cmf::kNuESelection );
614  else if(systName == "NueAcceptSignalKin2018FHCWgt") weighter.wgtFn = std::bind(&ShifterAndWeighter::NueAcceptSignalKin2018Weight, this, cmf::kFHC);
615  else if(systName == "NueAcceptSignalKin2018RHCWgt") weighter.wgtFn = std::bind(&ShifterAndWeighter::NueAcceptSignalKin2018Weight, this, cmf::kRHC);
616  else if(systName == "NueAcceptBkg2018FHCWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::NueAcceptBkg2018Weight, this, cmf::kFHC);
617  else if(systName == "NueAcceptBkg2018RHCWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::NueAcceptBkg2018Weight, this, cmf::kRHC);
618  else if(systName == "NueRelNormWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::RelNormWeight, this, systName, cmf::kNuESelection );
619  else if(systName == "NueFDRockMuonWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::NueFDRockMuonWeight, this);
620  else if(systName == "OtherGENIEWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
621  else if(systName == "RadCorrNueWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::RadCorrNueWeight, this);
622  else if(systName == "RadCorrNueBarWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::RadCorrNueBarWeight, this);
623  else if(systName == "RelNormWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::RelNormWeight, this, systName, cmf::kNuMuSelection);
624  else if(systName == "RHCNormWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::HornNormWeight, this, systName, cmf::kRHC );
625  else if(systName == "RockNormOnFidWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::RockNormOnFidWeight, this);
626  else if(systName == "RPACCQEWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::RPACCQEWeight, this);
627  else if(systName == "RPARESWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::RPARESWeight, this);
628  else if(systName == "RPACCQEshapeEnhWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::RPACCQEshapeEnhWeight, this);
629  else if(systName == "RPACCQEshapeSuppWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::RPACCQEshapeSuppWeight, this);
630  else if(systName == "SecondClassCurrWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::SecondClassCurrWeight, this);
631  else if(systName == "SumSmallXSecNumu2017Wgt") weighter.wgtFn = std::bind(&ShifterAndWeighter::SumSmallXSec2017Weight, this, systName, cmf::kNuMuSelection);
632  else if(systName == "SumSmallXSecNue2017Wgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::SumSmallXSec2017Weight, this, systName, cmf::kNuESelection);
633  else if(systName == "XsecCVWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::XSecCVPPFXWeight, this);
634  else if(systName == "TauScaleSystWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::TauScaleSystWeight, this);
635  else if(systName.find("BeamSystWgt") != std::string::npos){
636  auto detId = (systName.find("Far") != std::string::npos) ? cmf::kFARDET : cmf::kNEARDET;
637  weighter.wgtFn = std::bind(&ShifterAndWeighter::BeamSystWeight, this, systName, detId);
638  if(systName.find("2017") != std::string::npos) this->LoadBeamSystHists2017();
639  }
640  else if(systName.find("Calibration") != std::string::npos ||
641  systName.find("lightmodel") != std::string::npos ||
642  systName.find("ckv") != std::string::npos ){
643  weighter.wgtFn = std::bind(&ShifterAndWeighter::CalibSystWeight, this, systName);
644  this->LoadCalibrationSystRatios(calibFileName);
645  }
646  else if(systName.find("FluxPCA") != std::string::npos){
647  auto num = std::stol(systName.substr(7,1).c_str(), nullptr, 0);
648  weighter.wgtFn = std::bind(&ShifterAndWeighter::FluxPCAWeight, this, num);
649  this->LoadFluxBeamSystHists();
650  }
651 #ifdef PRODUCTION5
652  // 2020
653  else if (systName == "ZExpAxialFFSyst2020_eV1Wgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
654  else if (systName == "ZExpAxialFFSyst2020_eV2Wgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
655  else if (systName == "ZExpAxialFFSyst2020_eV3Wgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
656  else if (systName == "ZExpAxialFFSyst2020_eV4Wgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
657  else if (systName == "ZNormCCQE" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
658  else if (systName == "RESLowQ2SuppSyst2020Wgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
659  else if (systName == "DISvnCC1piWgt_2020Wgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
660  else if (systName == "hNFSISyst2020_MFPWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
661  else if (systName == "hNFSISyst2020_EV1Wgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
662  else if (systName == "MECEnuShapeSyst2020NuWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
663  else if (systName == "MECEnuShapeSyst2020AntiNuWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
664  else if (systName == "MECShapeSyst2020NuWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
665  else if (systName == "MECShapeSyst2020AntiNuWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
666  else if (systName == "MECInitStateNPFracSyst2020NuWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
667  else if (systName == "MECInitStateNPFracSyst2020AntiNuWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
668  else if (systName == "RadCorrNueWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
669  else if (systName == "RadCorrNuebarWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
670  else if (systName == "2ndClassCurrsWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
671  else if (systName == "GENIEPCA0Wgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
672  else if (systName == "GENIEPCA1Wgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
673  else if (systName == "GENIEPCA2Wgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
674  else if (systName == "GENIEPCA3Wgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
675  else if (systName == "GENIEPCA4Wgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
676  else if (systName == "GENIEPCA5Wgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
677  else if (systName == "GENIEPCA6Wgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
678  else if (systName == "GENIEPCA7Wgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
679  else if (systName == "GENIEPCA8Wgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
680  else if (systName == "GENIEPCA9Wgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
681  else if (systName == "GENIEPCA10Wgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
682  else if (systName == "GENIEPCA11Wgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
683  else if (systName == "UnCorrNDMuEScaleSyst2020Wgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
684  else if (systName == "UnCorrMuCatMuESyst2020Wgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
685  else if (systName == "PileupMuESyst2020Wgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
686  else if (systName == "CorrMuEScaleSyst2020Wgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
687  else if (systName == "NeutronVisEScalePrimariesSyst2018Wgt") weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
688  else if (systName == "MichelTaggingSyst2020Wgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
689  else if (systName == "Ana2020NormFHCWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
690  else if (systName == "Ana2020NormRHCWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
691  else if (systName == "TauScaleSystWgt" ) weighter.wgtFn = std::bind(&ShifterAndWeighter::GENIESystWeight, this, systName);
692 #else
693  else if(systName.find("GeniePCA") != std::string::npos){
694  auto num = std::stol(systName.substr(8,1).c_str(), nullptr, 0);
695  weighter.wgtFn = std::bind(&ShifterAndWeighter::GeniePCASystWeight, this, num);
696  this->LoadGeniePCASystHists();
697  }
698 #endif
699  else{
700  LOG_WARNING("ShifterAndWeighter")
701  << "Unable to determine weighting function for "
702  << systName;
703  continue;
704  }
705  LOG_DEBUG("ShifterAndWeighter")
706  << "Adding weighter "
707  << weighter.wgtName
708  << " to detector "
710 
711  weighters.push_back(weighter);
712 
713  } // the parameter was for a weight
714  else if(systName.find("Shifter") != std::string::npos){
715 
716  LOG_DEBUG("ShifterAndWeighter")
717  << "Adding shifter for "
718  << systName
719  << " to detector "
720  << det;
721 
722  this->FillShifterInfo(systName,
723  det);
724  } // end if this is a shifter
725  } // end loop over systematic parameters
726 
727  }
728 
729  //----------------------------------------------------------------------------
731  fhicl::ParameterSet const& parset)
732  {
733  // this method sets which systematics are to be used in calculating weights
734  // by creating a set of std::functions corresponding to the weighting functions
735  // that are to be turned on
736  // the code is a bit tedious, but it gets the job done
737 
738  auto calibFileName = parset.get<std::string>("CalibSystFile", "CovarianceMatrixFit/data/calibration/CalibrationSystematics2019.root");
739  auto calcString = parset.get<std::string>("OscCalculatorType", "PMNS");
740  fTrueEOscillationWeight = parset.get<bool >("TrueEOscCalc", true);
741 
742  cmf::CalcType_t calcType = cmf::kPMNS;
743  if (calcString == "PMNS_NSI" ) calcType = cmf::kNSI;
744  else if(calcString == "PMNS_Sterile") calcType = cmf::kSterile;
745 
746  // setup the oscillation calculator to use
747  this->OscillationCalculatorToUse(calcType);
748 
749  // go ahead and set the current values
750  this->SetCurrentVals(loc);
751 
753  loc.NDLoc,
754  calibFileName,
758  loc.FDLoc,
759  calibFileName,
762 
763 // for(auto const& itr : fShiftersToUseND){
764 // LOG_VERBATIM("ShifterAndWeighter")
765 // << "ND shifter param: "
766 // << itr.first;
767 // for(auto const& si : itr.second)
768 // LOG_VERBATIM("ShifterAndWeighter")
769 // << " associated shifter "
770 // << si.fShifterName
771 // << " "
772 // << si.fDetector;
773 // }
774 // for(auto const& itr : fShiftersToUseFD){
775 // LOG_VERBATIM("ShifterAndWeighter")
776 // << "FD shifter param: "
777 // << itr.first;
778 // for(auto const& si : itr.second)
779 // LOG_VERBATIM("ShifterAndWeighter")
780 // << " associated shifter "
781 // << si.fShifterName
782 // << " "
783 // << si.fDetector;
784 // }
785 
786  LOG_DEBUG("ShifterAndWeighter")
787  << "Done initialising shifters and weighters";
788 
789  }
790 
791  //----------------------------------------------------------------------------
793  cmf::DetType_t const& det)
794  {
795  uint8_t param;
796  if (name.find("HadE") != std::string::npos) param = cmf::kHad_RecoE;
797  else if(name.find("LepE") != std::string::npos) param = cmf::kLep_RecoE;
798  else if(name.find("NuE" ) != std::string::npos) param = cmf::kNu_RecoE;
799 
800  if (det == cmf::kNEARDET) fShiftersToUseND.emplace(name, param);
801  else if(det == cmf::kFARDET ) fShiftersToUseFD.emplace(name, param);
802  }
803 
804  //----------------------------------------------------------------------------
806  {
807  LOG_VERBATIM("ShifterAndWeighter")
808  << "Current shifter/weighter point:\n"
809  << fCurrentLoc;
810  }
811 
812  //----------------------------------------------------------------------------
813  // Use this version if there are any parameters with different values in the ND
814  // and FD
816  {
817  fCurrentLoc.NDLoc = curLoc.NDLoc;
818  fCurrentLoc.FDLoc = curLoc.FDLoc;
819  //this->ReportCurrentVals();
820 
821  this->SetOscillationVals();
822  }
823 
824  //----------------------------------------------------------------------------
825  // Only use this version if you are sure the ND and FD parameter values are the
826  // same when the same parameter is present and changing during a fit in both
827  // detectors
828  void ShifterAndWeighter::SetCurrentVals(std::vector<std::string> const& pars,
829  const double* vals)
830  {
831 
832  for(size_t pv = 0; pv < pars.size(); ++pv){
833  fCurrentLoc.FDLoc.find(pars[pv])->second.SetValue(vals[pv]);
834  fCurrentLoc.NDLoc.find(pars[pv])->second.SetValue(vals[pv]);
835  }
836  //this->ReportCurrentVals();
837 
838  this->SetOscillationVals();
839  }
840 
841  //----------------------------------------------------------------------------
843  {
844  // get the oscillation parameters from the current point, we haven't set
845  // the detector yet, but the only thing that changes is the
846  // baseline which we change with the SetCurrentEvent method
847  for(auto const& itr : fCurrentLoc.FDLoc){
848 
850  if (itr.first == "L" ) fOscCalc->SetL (itr.second.Value());
851  else if(itr.first == "Rho" ) fOscCalc->SetRho (itr.second.Value());
852  else if(itr.first == "Dmsq21") fOscCalc->SetDmsq21(itr.second.Value());
853  else if(itr.first == "Dmsq32") fOscCalc->SetDmsq32(itr.second.Value());
854  else if(itr.first == "Th12" ) fOscCalc->SetTh12 (itr.second.Value());
855  else if(itr.first == "Th13" ) fOscCalc->SetTh13 (itr.second.Value());
856  else if(itr.first == "Th23" ) fOscCalc->SetTh23 (itr.second.Value());
857  else if(itr.first == "dCP" ) fOscCalc->SetdCP (itr.second.Value());
858  if(fCalcType == cmf::kNSI){
859  //Eps = NSI epsilon parameters
860  if (itr.first == "Eps_ee" ) dynamic_cast<osc::OscCalcPMNS_NSI*>(fOscCalc)->SetEps_ee (itr.second.Value());
861  else if(itr.first == "Eps_emu" ) dynamic_cast<osc::OscCalcPMNS_NSI*>(fOscCalc)->SetEps_emu (itr.second.Value());
862  else if(itr.first == "Eps_etau" ) dynamic_cast<osc::OscCalcPMNS_NSI*>(fOscCalc)->SetEps_etau (itr.second.Value());
863  else if(itr.first == "Eps_mumu" ) dynamic_cast<osc::OscCalcPMNS_NSI*>(fOscCalc)->SetEps_mumu (itr.second.Value());
864  else if(itr.first == "Eps_mutau" ) dynamic_cast<osc::OscCalcPMNS_NSI*>(fOscCalc)->SetEps_mutau (itr.second.Value());
865  else if(itr.first == "Eps_tautau" ) dynamic_cast<osc::OscCalcPMNS_NSI*>(fOscCalc)->SetEps_tautau (itr.second.Value());
866  else if(itr.first == "Delta_emu" ) dynamic_cast<osc::OscCalcPMNS_NSI*>(fOscCalc)->SetDelta_emu (itr.second.Value());
867  else if(itr.first == "Delta_etau" ) dynamic_cast<osc::OscCalcPMNS_NSI*>(fOscCalc)->SetDelta_etau (itr.second.Value());
868  else if(itr.first == "Delta_mutau") dynamic_cast<osc::OscCalcPMNS_NSI*>(fOscCalc)->SetDelta_mutau(itr.second.Value());
869  } // end if NSI oscillations
870  } // end if PMNS or NSI oscillations
871  else if(fCalcType == cmf::kSterile){
872  if (itr.first == "Th12" ) dynamic_cast<osc::OscCalcSterile*>(fOscCalc)->SetAngle(1, 2, itr.second.Value());
873  else if(itr.first == "Th23" ) dynamic_cast<osc::OscCalcSterile*>(fOscCalc)->SetAngle(2, 3, itr.second.Value());
874  else if(itr.first == "Th24" ) dynamic_cast<osc::OscCalcSterile*>(fOscCalc)->SetAngle(2, 4, itr.second.Value());
875  else if(itr.first == "Th34" ) dynamic_cast<osc::OscCalcSterile*>(fOscCalc)->SetAngle(3, 4, itr.second.Value());
876  else if(itr.first == "dCP" ) dynamic_cast<osc::OscCalcSterile*>(fOscCalc)->SetDelta(1, 3, itr.second.Value());
877  else if(itr.first == "Dmsq21") dynamic_cast<osc::OscCalcSterile*>(fOscCalc)->SetDm (2, itr.second.Value());
878  else if(itr.first == "Dmsq32") dynamic_cast<osc::OscCalcSterile*>(fOscCalc)->SetDm (3, itr.second.Value());
879  else if(itr.first == "Dmsq41") dynamic_cast<osc::OscCalcSterile*>(fOscCalc)->SetDm (4, itr.second.Value());
880  else if(itr.first == "L" ) dynamic_cast<osc::OscCalcSterile*>(fOscCalc)->SetL ( itr.second.Value());
881  else if(itr.first == "Rho" ) dynamic_cast<osc::OscCalcSterile*>(fOscCalc)->SetRho ( itr.second.Value());
882  } // end if sterile oscillations
883  } // end loop over FD location
884  }
885 
886  //----------------------------------------------------------------------------
888  cmf::MetaData const& md)
889  {
890  fCurrentEvent = curEvent.get();
891  fCurrentMD = md;
892 
893  LOG_DEBUG("SetCurrentEvent")
894  << " Metadata is "
895  << md.ToString();
896 
897  // check the validity of the mcVals before trying to get weights, if it is
898  // a nullptr, we have a problem.
899  if( !fCurrentEvent->MCValsFilled() && md.isMC)
900  throw cet::exception("ShifterAndWeighter")
901  << "Event does not have a valid MCVarVals unique_ptr, bail";
902 
903  // the data values better always be valid
904  if( !fCurrentEvent->DataValsFilled() )
905  throw cet::exception("ShifterAndWeighter")
906  << "Event does not have a valid DataVarVals unique_ptr, bail";
907 
908  // reset the oscillation baseline based on the detector
909  fOscCalc->SetL(this->CurrentParameterValue("L"));
910 
911  }
912 
913  //----------------------------------------------------------------------------
915  {
916  return (fCurrentMD.detector == cmf::kNEARDET) ?
918  }
919 
920  //----------------------------------------------------------------------------
921  float ShifterAndWeighter::ShiftForParameter(std::map<std::string, uint8_t> const& shifterMap,
923  uint8_t const& param)
924  {
925 
926  float totalShift = 1.;
927  float shiftScale;
928  float curShift;
929 
930  for(auto const& itr : shifterMap){
931 
932  if(itr.second != param) continue;
933 
934  auto const& sigma = loc.find(itr.first)->second.Sigma();
935 
936  curShift = this->CurrentParameterValue(itr.first);
937 
938  if(itr.first.find("MuonCatcherEShifter") != std::string::npos &&
940  // The plan is to determine how much of the energy is contained in the
941  // ND muon catcher and then scale that amount by the product of the
942  // number of sigma and the shift in the scale due to a 10.5 cm shift in
943  // the muon catcher pathlength The shift in pathlength is accounted as a
944  // shift in muon catcher energy from range.
945  // find the sigma for the muon catcher from the TF1s first, then multiply that
946  // by the fraction of the total energy that is in the muon catcher
947  //
948  // For the near detector if we are using relative energy shifts
949  // we want to shift things down by 1/2 sigma
951  shiftScale *= curShift;
952  shiftScale *= (curShift > 0.) ? fNDMuonCatcherEShiftPlus->Eval(fCurrentEvent->DataVals().val_at(cmf::kLep_RecoE, fCurrentMD))
954  shiftScale += 1.;
955  }
956  else if(itr.first.find("RelRecoLepEShifter") != std::string::npos){
957  // if this shifter is the relative leptonic one, set the scale
958  // to be 0.5 in the FD and -0.5 in the ND, ie the FD is shifted up by
959  // 1/2 sigma and the ND is shifted down by 1/2 sigma
960  shiftScale = sigma * curShift;
961  shiftScale *= (fCurrentMD.detector == cmf::kNEARDET) ? -0.5 : 0.5;
962  shiftScale += 1.;
963  }
964  else
965  shiftScale = 1. + sigma * curShift;
966 
967  // LOG_VERBATIM("ShifterAndWeighterTotalShift")
968  // << " associated shifter "
969  // << infoItr.fShifterName
970  // << " "
971  // << infoItr.fDetector
972  // << " "
973  // << fCurrentMD.detector
974  // << " "
975  // << sigma
976  // << " "
977  // << this->CurrentParameterValue(infoItr.fShifterName)
978  // << " "
979  // << totalShift;
980  totalShift *= shiftScale;
981  }
982 
983  return totalShift;
984  }
985 
986  //----------------------------------------------------------------------------
988  cmf::MetaData const& md,
989  cmf::WeightType_t const wgtType)
990  {
991  // get ready to use the passed in event
992  this->SetCurrentEvent(curEvent, md);
993 
994  // set the weight to be the fake weight value, which is 1 for
995  // data.
996  double wgt = curEvent->DataVals().val_at(cmf::kFake_Weight, md);
997 
998  // if the mcVals is null or the metadata claims it isn't MC,
999  // return the fake weight for the event. If the event is data
1000  // then the value is always 1, if it is fake data, then the
1001  // value is the proper weighting
1002  if(!curEvent->MCValsFilled() || !md.isMC)
1003  return wgt;
1004 
1005  // if(md.selectionType == cmf::kNuESelectionPeripheral)
1006  // LOG_VERBATIM("ShifterAndWeighter")
1007  // << "weight from systematics is "
1008  // << wgt;
1009 
1010  // check if we only want weights from the systematic parameters or only
1011  // from the oscillation parameters
1012  if(wgtType == cmf::kSystWgtOnly){
1013  wgt *= this->TotalWeightFromWeighters();
1014  }
1015  else if(wgtType == cmf::kOscWgtOnly){
1016  wgt *= this->OscillationWeight();
1017  }
1018  else
1019  wgt *= this->OscillationWeight() * this->TotalWeightFromWeighters();
1020 
1021  // if(md.selectionType == cmf::kNuESelectionPeripheral)
1022  // LOG_VERBATIM("ShifterAndWeighter")
1023  // << "weight from oscillations is "
1024  // << this->OscillationWeight();
1025 
1026  LOG_DEBUG("ShifterAndWeighter")
1027  << "Oscillation weight is "
1028  << wgt
1029  << " energy: "
1031  << " pdg: "
1033  << " orig pdg: "
1035  << " detector: "
1036  << fCurrentMD.detector;
1037 
1038  // if(md.selectionType == cmf::kNuESelectionPeripheral)
1039  // LOG_VERBATIM("ShifterAndWeighter")
1040  // << "total weight is "
1041  // << wgt;
1042 
1043  return wgt;
1044  }
1045 
1046  //----------------------------------------------------------------------------
1048  {
1049  double wgt = 1.;
1050 
1052 
1053  //std::string ss;
1054  for(auto itr : weighters){
1055  // ss += "\t" + itr.wgtName +
1056  // " " + std::to_string(itr.wgtFn()) +
1057  // " " + std::to_string(wgt) +
1058  // " " + std::to_string(wgt * itr.wgtFn()) +
1059  // "\n";
1060 
1061  wgt *= itr.wgtFn();
1062  }
1063  // LOG_VERBATIM("SAWWeighterInfo")
1064  // << ss;
1065 
1066  return wgt;
1067  }
1068 
1069  //----------------------------------------------------------------------------
1071  {
1072  // nothing to change
1073  if(fOscCalc != nullptr && fCalcType == calcType) return;
1074  else{
1075  LOG_WARNING("ShifterAndWeighter")
1076  << "resetting oscillation calculator to type "
1077  << calcType
1078  << " from "
1079  << fCalcType
1080  << " are you sure you meant to change types in this job?";
1081 
1082  delete fOscCalc;
1083  fOscCalc = nullptr;
1084  }
1085 
1086  fCalcType = calcType;
1087 
1089  else if(fCalcType == cmf::kNSI ) fOscCalc = new osc::OscCalcPMNS_NSI();
1090  else if(fCalcType == cmf::kSterile){
1091  fOscCalc = new osc::OscCalcSterile();
1092  dynamic_cast<osc::OscCalcSterile*>(fOscCalc)->SetNFlavors(4);
1093  }
1094  }
1095 
1096  //----------------------------------------------------------------------------
1098  {
1099  // don't bother doing anything if we are looking at the near detector
1100  // and standard oscillations, it is just a waste of time.
1101  if(fCalcType == cmf::kPMNS &&
1102  fCurrentMD.detector == cmf::kNEARDET) return 1.;
1103 
1105 
1106  return this->OscillationWeightBinCenter();
1107  }
1108  //----------------------------------------------------------------------------
1110  {
1111  double E = fCurrentEvent->MCVals().val_at(cmf::kTrueE);
1113  int pdgOrig = fCurrentEvent->MCVals().val_at(cmf::kTruePDGOrig);
1114 
1115  // only attempt to weight neutrinos for oscillations
1116  if(std::abs(pdg) != 12 &&
1117  std::abs(pdg) != 14 &&
1118  std::abs(pdg) != 16) return 1.;
1119 
1120  return fOscCalc->P(pdgOrig, pdg, E);
1121 
1122  }
1123 
1124  //----------------------------------------------------------------------------
1126  {
1127  if(fTrueEnergyBins.empty()){
1128  const int numTrueEnergyBins = 100;
1129 
1130  // How many edges to generate. Allow room for 0-Emin bin
1131  const double N = numTrueEnergyBins-1;
1132  const double A = N * 0.5; // 500 MeV: there's really no events below there
1133 
1134  fTrueEnergyBins.insert(0);
1135  for(int i = 1; i < N - 1; ++i) fTrueEnergyBins.insert(A/i);
1136  fTrueEnergyBins.insert(120.);
1137  }
1138 
1139  double E = fCurrentEvent->MCVals().val_at(cmf::kTrueE);
1141  int pdgOrig = fCurrentEvent->MCVals().val_at(cmf::kTruePDGOrig);
1142 
1143  // get the relevant bin for the current neutrino, then get the
1144  // energy at the center of that bin
1145  double eUpper = *(fTrueEnergyBins.rbegin());
1146  double eLower = *(fTrueEnergyBins.begin());
1147 
1148  // find the bin edge that goes after the value of E
1149  auto itr = fTrueEnergyBins.upper_bound(E);
1150  if(itr != fTrueEnergyBins.end()) eUpper = *itr;
1151  else --itr;
1152 
1153  // now find the edge before by just decrementing the iterator
1154  --itr;
1155  if(itr != fTrueEnergyBins.end()) eLower = *itr;
1156 
1157  // get the midpoint between the upper and lower values,
1158  // ie take their mean
1159  E = 0.5 * (eUpper + eLower);
1160 
1161  LOG_DEBUG("ShifterAndWeighter")
1162  << "true energy: "
1164  << " central bin energy "
1165  << E;
1166 
1167  // only attempt to weight neutrinos for oscillations
1168  if(std::abs(pdg) != 12 &&
1169  std::abs(pdg) != 14 &&
1170  std::abs(pdg) != 16) return 1.;
1171 
1172  return fOscCalc->P(pdgOrig, pdg, E);
1173 
1174  }
1175 
1176  //----------------------------------------------------------------------------
1178  std::map<float, float> const& wgtsBySigma)
1179  {
1180  auto mapEnd = wgtsBySigma.end();
1181 
1182  // no way this is a valid number and is obviously stupid
1183  // not using a std::numeric_limits value to avoid rounding to 0
1184  double badVal = -99999999999999.;
1185  double minus2Val = (wgtsBySigma.find(-2.) != mapEnd) ? static_cast<double>(wgtsBySigma.find(-2.)->second) : badVal;
1186  double minus1Val = (wgtsBySigma.find(-1.) != mapEnd) ? static_cast<double>(wgtsBySigma.find(-1.)->second) : badVal;
1187  // double zeroVal = 1.;
1188  double zeroVal = (wgtsBySigma.find( 0.) != mapEnd) ? static_cast<double>(wgtsBySigma.find( 0.)->second) : 1.;
1189  double plus1Val = (wgtsBySigma.find( 1.) != mapEnd) ? static_cast<double>(wgtsBySigma.find( 1.)->second) : badVal;
1190  double plus2Val = (wgtsBySigma.find( 2.) != mapEnd) ? static_cast<double>(wgtsBySigma.find( 2.)->second) : badVal;
1191  double slope = 0.;
1192  double intercept = 0.;
1193 
1194  //Debugging RPA weights
1195  LOG_DEBUG("ShifterAndWeighter")
1196  << "\tCalcLinearInterpWgt: "
1197  << " -2:" << minus2Val
1198  << ", -1:" << minus1Val
1199  << ", 0:" << zeroVal
1200  << ", 1:" << plus1Val
1201  << ", 2:" << plus2Val;
1202 
1203  /*
1204  if(wgtsBySigma.find(-1.) != mapEnd)
1205  LOG_VERBATIM("ShifterAndWeighter")
1206  << "\t\t 0 val from "
1207  << wgtsBySigma.find(0.)->second;
1208  */
1209 
1210  // The following code assumes that at least the -1 sigma, +1 sigma
1211  // values are in the map. If they are not, we have a problem
1212  if(minus1Val == badVal ||
1213  plus1Val == badVal){
1214  LOG_WARNING("ShifterAndWeighterCalcLinearInterpWgt")
1215  << "No values associated with either -1 sigma or +1 sigma "
1216  << "so we cannot calculate a weight, return 1.";
1217  return 1.;
1218  }
1219 
1220  if(sigma <= -1){
1221  // It is possible there is no -2 sigma value from the map,
1222  // in that case just interpolate from 0 to -1 sigma
1223  if(minus2Val != badVal){
1224  slope = (minus1Val - minus2Val);
1225  intercept = minus2Val - slope * (-2.);
1226  }
1227  else{
1228  slope = (zeroVal - minus1Val);
1229  intercept = minus1Val - slope * (-1.);
1230  }
1231  }
1232  else if(sigma > -1. && sigma <= 0.){
1233  slope = (zeroVal - minus1Val);
1234  intercept = minus1Val - slope * (-1.);
1235  }
1236  else if(sigma > 0. && sigma <= 1.){
1237  slope = (plus1Val - zeroVal);
1238  intercept = plus1Val - slope * 1.;
1239  }
1240  else if(sigma > 1.){
1241  // It is possible there is no +2 sigma value from the map,
1242  // in that case just interpolate from 0 to +1 sigma
1243  if(plus2Val != badVal){
1244  slope = (plus2Val - plus1Val);
1245  intercept = plus2Val - slope * 2.;
1246  }
1247  else{
1248  slope = (plus1Val - zeroVal);
1249  intercept = plus1Val - slope * 1.;
1250  }
1251  }
1252 
1253  double wgt = sigma * slope + intercept;
1254 
1255  // Don't let weights be negative
1256  if(wgt < 0) wgt = 0.;
1257 
1258  return wgt;
1259  }
1260 
1261  //--------------------------------------------------------------------------
1263  {
1264  // Someone already called us
1265  if(!fBeamSyst2017.fBeamSystHists.empty()) return;
1266 
1268 
1269  // TODO: make this a tunable parameter
1270  std::string shortName("enuMF_FHC_totErr");
1271 
1272  // constructor decides if initialized value is a path or an environment variable
1273  cet::search_path sp("FW_SEARCH_PATH");
1274  if ( !sp.find_file(std::string("CAFAna/data/beam/TABeamSyst.root"), fileName) ) {
1275  throw cet::exception("BeamSystFileError")
1276  << "cannot find Beam Systematics file CAFAna/data/beam/TABeamSyst.root "
1277  << fileName
1278  << " bail ungracefully\n"
1279  << "(RJN Guess: Did you run setup for this release?)\n"
1280  << __FILE__ << ":" << __LINE__;
1281  }
1282 
1283  fBeamSyst2017.SetFileAndHistNames(fileName, shortName);
1285  }
1286 
1287  //----------------------------------------------------------------------------
1289  {
1290  // Someone already called us
1291  if(!fVecFluxPrincipals.empty()) return;
1292 
1293  cet::search_path sp("FW_SEARCH_PATH");
1294 
1295  // the central value smooth PPFX weights appear to be no longer used
1296 // std::string fileNameCV;
1297 // std::string shortNameCV("ppfx_cv");
1298 //
1299 // // constructor decides if initialized value is a path or an environment variable
1300 // if ( !sp.find_file(std::string("CAFAna/data/beam/ppfx_smooth_multiverse_weights_2018.root"), fileNameCV) ) {
1301 // throw cet::exception("BeamSystFileError")
1302 // << "cannot find Beam Systematics file CAFAna/data/beam/ppfx_smooth_multiverse_weights_2018.root "
1303 // << fileNameCV
1304 // << " bail ungracefully\n"
1305 // << __FILE__ << ":" << __LINE__;
1306 // }
1307 //
1308 // fPPFXSmoothCVWgt.SetFileAndHistNames(fileNameCV,shortNameCV);
1309 // fPPFXSmoothCVWgt.LoadBeamSystHists();
1310 
1311  std::string fileNamePCA;
1312  std::string shortNamePCA;
1313  if ( !sp.find_file(std::string("CAFAna/data/beam/ppfx_hadp_beam_pc_shifts_fn_2018.root"), fileNamePCA) ) {
1314  throw cet::exception("BeamSystFileError")
1315  << "cannot find Beam Systematics file CAFAna/data/beam/ppfx_hadp_beam_pc_shifts_fn_2018.root "
1316  << fileNamePCA
1317  << " bail ungracefully\n"
1318  << __FILE__ << ":" << __LINE__;
1319  }
1320 
1321  for(int index = 0; index < 5; ++index){
1322  ///For now 5 is just an arbitrary number, should be configurable
1323  std::string PCIdxStr = TString::Format("%02d", index).Data();
1324  shortNamePCA="ppfx_hadp_beam_pc"+PCIdxStr;
1325 
1326  fVecFluxPrincipals.emplace_back(BeamSyst(fileNamePCA,shortNamePCA));
1327  fVecFluxPrincipals[index].LoadBeamSystHists();
1328 
1329  LOG_DEBUG("ShifterAndWeighter")
1330  << "setting up BeamSyst object for "
1331  << fileNamePCA
1332  << " "
1333  << shortNamePCA
1334  << " "
1335  << fVecFluxPrincipals.size();
1336  }
1337 
1338  //Can also add a similar loop for the univeses
1339  //will leave this commented out for now
1340  // std::string shortNameUniv;
1341  // for(int index=0;index<100;index++) {
1342  // std::string UniIdxStr = TString::Format("%02d", index).Data();
1343  // shortNameUniv="ppfx_univ"+UniIdxStr;
1344  // fVecPPFXUniverses.push_back(BeamSyst(fileNameCV,shortNameUniv));
1345  // }
1346 
1347 
1348  }
1349 
1350  //------------------------------------------------------------------------
1352  {
1353  fCalibSyst.LoadCalibSystRatios(CalibFilename);
1354  }
1355 
1356  //------------------------------------------------------------------------
1358  // someone already called us
1359  if(!fVecGeniePCASysts.empty()) return;
1360 
1361  for(int index = 0; index < 5; ++index) {
1362  // 5 chosen arbitrarily, should probably make configurable at some point
1363  std::string PCIdxStr = TString::Format("%02d", index).Data();
1364  std::string shortName = "genie_small_pc"+PCIdxStr;
1365  fVecGeniePCASysts.emplace_back(GeniePCASyst());
1366 
1367  fVecGeniePCASysts[index].LoadGeniePCASystHists(shortName);
1368  }
1369 
1370  }
1371 
1372  //----------------------------------------------------------------------------
1374  cmf::SelectionType_t const& sel)
1375  {
1376  // check for the correct selection for the current event with this function
1377  // sel is one of the two generic selection types, kNuMuSelection or kNuESelection
1378  if(fCurrentMD.IsNuMuSelected() &&
1379  sel == cmf::kNuMuSelection &&
1380  wgtType.find("NuMu") == std::string::npos) return 1.;
1381  else if(fCurrentMD.IsNuESelected() &&
1382  sel == cmf::kNuESelection &&
1383  wgtType.find("Nue") == std::string::npos) return 1.;
1384  else if(fCurrentMD.IsNCSelected() &&
1385  sel == cmf::kNCSelection &&
1386  wgtType.find("NC") == std::string::npos) return 1.;
1387 
1388  auto const& sigma = cmf::ParameterUtility::Instance()->ParameterInfo(wgtType).Sigma();
1389 
1390  LOG_DEBUG("AbsNormWeight")
1391  << " Sigma Value"
1392  << sigma;
1393 
1394  return 1. + sigma * this->CurrentParameterValue(wgtType);
1395 
1396  }
1397 
1398  //----------------------------------------------------------------------------
1400  cmf::DetType_t const& det)
1401  {
1402  if(det != fCurrentMD.detector) return 1.;
1403 
1404  std::map<float, float> wgtBySigma;
1405 
1406  if(sysType.find("2017") != std::string::npos) {
1407  wgtBySigma = fBeamSyst2017.CalcBeamSystWeights(det,
1408  fCurrentMD.BeamType(),
1411  }
1412 
1413  return this->CalcLinearInterpWgt(this->CurrentParameterValue(sysType),
1414  wgtBySigma);
1415  }
1416 
1417  //----------------------------------------------------------------------------
1419  cmf::DetType_t const& det)
1420  {
1421  if(det != fCurrentMD.detector) return 1.;
1422 
1423  auto const& sigma = cmf::ParameterUtility::Instance()->ParameterInfo(wgtType).Sigma();
1424 
1425  return 1. + sigma * this->CurrentParameterValue(wgtType);
1426  }
1427 
1428  //-------------------------------------------------------------------------
1430  {
1431  // The calibration shape weights are from the calib-shift-func files
1432 
1433  // The relative and absolute calibration weights are based off of the
1434  // calib-shift-xyview-neg(pos)-offset data sets
1435  // For the absolute calibration both ND and FD are shifted
1436  // according to the same underlying distribution, ie calib-shift-xyview-pos or -neg
1437  // For the relative calibration the ND is shifted according to one, but the FD is
1438  // shifted according to the other. The other wrinkle is that the relative
1439  // calibration only shifts by 1/2 sigma in each detector
1440 
1441  // need to butcher these wgtType string to get the weights by sigma
1442  double sigShift = this->CurrentParameterValue(wgtType);
1443 
1444  // first off, figure out if we are a relative or absolution calibration shift
1445  auto systDataSetName = wgtType.substr(0, wgtType.size() - 3);
1446 
1447  //calib-shift-xyview-neg-offsetWgt calib-shift-func
1448  if(systDataSetName == "RelativeCalibration"){
1449  sigShift *= 0.5;
1451  systDataSetName = (sigShift < 0.) ? "calib-shift-xyview-neg-offset" : "calib-shift-xyview-pos-offset";
1452  else
1453  systDataSetName = (sigShift < 0.) ? "calib-shift-xyview-pos-offset" : "calib-shift-xyview-neg-offset";
1454  }
1455  else if(systDataSetName == "AbsoluteCalibration"){
1456  if(sigShift < 0)
1457  systDataSetName = "calib-shift-xyview-neg-offset";
1458  else
1459  systDataSetName = "calib-shift-xyview-pos-offset";
1460  }
1461  else if(systDataSetName == "CalibrationShape"){
1462  systDataSetName = "calib-shift-func";
1463  }
1464  else if (systDataSetName == "Calibration2020")
1465  systDataSetName = "calibration";
1466 
1467  // now we can get the amount of a 1 sigma shift for the given
1468  // reconstructed energy
1469  auto oneSigmaWgt = fCalibSyst.CalibSystRatio(systDataSetName,
1470  fCurrentMD,
1472 
1473  // the calibration systematics are based on the ratio of 1sigma shifted histograms to the nominal
1474  // A value of oneSigmaWgt > 1 means (1 - oneSigmaWgt) < 0 and 1 - (1 - oneSigmaWgt) will give the
1475  // correct weighting of just oneSigmaWgt if sigShift = 1. If oneSigmaWgt < 1 then the returned
1476  // expression still does the correct thing.
1477 
1478  return 1. - (1. - oneSigmaWgt) * sigShift;
1479  }
1480 
1481  //-----------------------------------------------------------------------------
1483  cmf::InteractionType_t const& intType)
1484  {
1485  // if event isn't coherent return 1
1489 
1490  // check that we have the right interaction type for this weight
1491  if((wgtType.find("NC") == std::string::npos && fCurrentMD.interactionType == cmf::kNC) ||
1492  (wgtType.find("CC") != std::string::npos && fCurrentMD.interactionType == cmf::kNC) ) return 1.;
1493 
1494  auto const& sigma = cmf::ParameterUtility::Instance()->ParameterInfo(wgtType).Sigma();
1495 
1496  auto parVal = this->CurrentParameterValue(wgtType);
1497 
1498  double weight = 1. + sigma * parVal;
1499 
1500  LOG_DEBUG("COHCCScale")
1501  << wgtType
1502  << weight;
1503 
1504  if(weight >= 0.) return weight;
1505  else return 0.;
1506  }
1507 
1508  //----------------------------------------------------------------------------
1510  {
1511  // cosmic muon normalization is only for the far detector
1512  if(fCurrentMD.detector == cmf::kNEARDET) return 1.;
1513 
1514  auto const& sigma = cmf::ParameterUtility::Instance()->ParameterInfo("CosmicMuNormWgt").Sigma();
1516  return 1. + sigma * this->CurrentParameterValue("CosmicMuNormWgt");
1517  }
1518 
1519  return 1.0;
1520  }
1521 
1522  //----------------------------------------------------------------------------
1524  {
1527  auto const& sigma = cmf::ParameterUtility::Instance()->ParameterInfo("CosmicOverlayNormWgt").Sigma();
1528 
1529  return 1. + sigma * this->CurrentParameterValue("CosmicOverlayNormWgt");
1530  }
1531 
1532  return 1.0;
1533  }
1534 
1535  //----------------------------------------------------------------------------
1537  {
1538  auto parVal = this->CurrentParameterValue("DISvnCC1piWgt");
1539  if(parVal > 0) return 1. + parVal * (1 - fCurrentEvent->MCVals().val_at(cmf::kDISvnCC1pi_Weight));
1540 
1541  return 1.;
1542  }
1543 
1544  //----------------------------------------------------------------------------
1546  {
1547  if(index > fVecFluxPrincipals.size() ) return 1;
1548 
1549  /*if(fCurrentMD.BeamType() == cmf::kUnknownBeam &&
1550  fCurrentEvent->eventId) {
1551  if(fCurrentEvent->eventId->run)
1552  LOG_VERBATIM("SJB") << "UNKNOWN BEAM!"
1553  << "Current run: "
1554  << fCurrentEvent->eventId->run;
1555  else
1556  LOG_VERBATIM("SJB") << "Can't find run information! What!";
1557 
1558  }*/
1559 
1560  auto wgtBySigma = fVecFluxPrincipals[index].CalcBeamSystWeights(fCurrentMD.detector,
1561  fCurrentMD.BeamType(),
1562  fCurrentEvent->MCVals().val_at(cmf::kTrueE), //True Energy
1564 
1565 
1566 // for(auto const& itr : wgtBySigma)
1567 // LOG_DEBUG("ShifterAndWeighter")
1568 // << "flux weights for PCA "
1569 // << index
1570 // << " "
1571 // << fCurrentEvent->MCVals().val_at(cmf::kTrueE)
1572 // << " "
1573 // << itr.first
1574 // << " : "
1575 // << itr.second;
1576 
1577  std::string systname;
1578  std::string PCIdxStr = TString::Format("%01d", (int)index).Data();
1579  systname="FluxPCA"+PCIdxStr+"Wgt";
1580 
1581  LOG_DEBUG("ShifterAndWeighter")
1582  << "flux PCA weight for "
1583  << systname
1584  << " and "
1585  << this->CurrentParameterValue(systname)
1586  << " sigma is "
1587  << this->CalcLinearInterpWgt(this->CurrentParameterValue(systname),
1588  wgtBySigma);
1589  return this->CalcLinearInterpWgt(this->CurrentParameterValue(systname),
1590  wgtBySigma);
1591  }
1592 
1593  //-------------------------------------------------------------------------
1595  {
1598 
1599  auto wgtBySigma = fVecGeniePCASysts[index].CalcGeniePCASystWeights(fCurrentMD, mode, energy);
1600 
1601  std::string systName;
1602  std::string PCIdxStr = TString::Format("%01d", (int)index).Data();
1603  systName = "GeniePCA" + PCIdxStr + "Wgt";
1604 
1605  return this->CalcLinearInterpWgt(this->CurrentParameterValue(systName),
1606  wgtBySigma);
1607  }
1608 
1609  //----------------------------------------------------------------------------
1611  {
1612  // return 1 if this is a cosmic muon
1613  if(fCurrentMD.interactionType == cmf::kCosmicMuon) return 1.;
1614 
1615  // if this weight has NC or CC in the name and this interaction type
1616  // is not that, return 1
1617  if((systName.find("CC") != std::string::npos && fCurrentMD.interactionType == cmf::kNC) ||
1618  (systName.find("NC") == std::string::npos && fCurrentMD.interactionType == cmf::kNC) ) return 1.;
1619 
1620  // may need to come up with default for cosmic muons, however I think the fact that
1621  // the default weights for every GENIE variable is a collection of 1's should
1622  // mean this is OK
1623  std::map<float, float> wgtBySigma = fCurrentEvent->MCVals().SystVarWgts(cmf::VarNameToKey(systName.substr(0, systName.size() - 3)));
1624 
1625  // for(auto itr : wgtBySigma)
1626  //LOG_VERBATIM("GENIESystWeight")
1627  //<< systName
1628  //<< " weights "
1629  //<< itr.first
1630  //<< " "
1631  //<< itr.second;
1632 
1633  return this->CalcLinearInterpWgt(this->CurrentParameterValue(systName),
1634  wgtBySigma);
1635  }
1636 
1637 // //----------------------------------------------------------------------------
1638 // double ShifterAndWeighter::HornCurrentWeight()
1639 // {
1640 // double wgt = 1;
1641 //
1642 // if(fHornCurrentHists.empty()) {
1643 // std::string fileName;
1644 // cet::search_path sp("FW_SEARCH_PATH");
1645 // if ( !sp.find_file(std::string("CovarianceMatrixFit/data/horncurrent/190kAHornCurrentHists.root"), fileName) ) {
1646 // throw cet::exception("HornCurrentHistError")
1647 // << "Cannot find horn current histograms.";
1648 // }
1649 //
1650 // TFile fin(fileName.c_str(), "READ");
1651 // if(fin.IsZombie()) {
1652 // throw cet::exception("HornCurrentHistError")
1653 // << "Couldn't open horn current file " << fileName;
1654 // }
1655 //
1656 // TString histNameND = "rat190ND";
1657 // TString histNameFD = "rat190FD";
1658 // TString histNameNDbar = "rat190NDan";
1659 // TString histNameFDbar = "rat190FDan";
1660 //
1661 // TH1D *hHornCurrentND = (TH1D*)fin.Get(histNameND)->Clone();
1662 // TH1D *hHornCurrentFD = (TH1D*)fin.Get(histNameFD)->Clone();
1663 // TH1D *hHornCurrentNDbar = (TH1D*)fin.Get(histNameNDbar)->Clone();
1664 // TH1D *hHornCurrentFDbar = (TH1D*)fin.Get(histNameFDbar)->Clone();
1665 //
1666 // hHornCurrentND ->SetDirectory(nullptr);
1667 // hHornCurrentFD ->SetDirectory(nullptr);
1668 // hHornCurrentNDbar->SetDirectory(nullptr);
1669 // hHornCurrentFDbar->SetDirectory(nullptr);
1670 //
1671 // fHornCurrentHists.push_back(hHornCurrentND );
1672 // fHornCurrentHists.push_back(hHornCurrentFD );
1673 // fHornCurrentHists.push_back(hHornCurrentNDbar);
1674 // fHornCurrentHists.push_back(hHornCurrentFDbar);
1675 // }
1676 //
1677 // double E = fCurrentEvent->MCVals().val_at(cmf::kTrueE);
1678 // int beamType = fCurrentMD.BeamType();
1679 // int pid = fCurrentEvent->MCVals().val_at(cmf::kTruePDG);
1680 //
1681 // if(fCurrentMD.interactionType == cmf::kCosmicMuon) return wgt;
1682 // if(E >= 20) return wgt;
1683 // if( (beamType == cmf::kFHC && pid < 0 && E >= 8) ||
1684 // (beamType == cmf::kRHC && pid > 0 && E >= 8) ) return wgt;
1685 // if(beamType != cmf::kFHC && beamType != cmf::kRHC) return wgt;
1686 //
1687 // if(fCurrentMD.detector == cmf::kNEARDET) {
1688 // if( (beamType == cmf::kFHC && pid > 0) || (beamType == cmf::kRHC && pid < 0) )
1689 // wgt = fHornCurrentHists.at(0)->Interpolate(E);
1690 // if( (beamType == cmf::kFHC && pid < 0) || (beamType == cmf::kRHC && pid > 0) )
1691 // wgt = fHornCurrentHists.at(2)->Interpolate(E);
1692 // }
1693 // if(fCurrentMD.detector == cmf::kFARDET) {
1694 // if( (beamType == cmf::kFHC && pid > 0) || (beamType == cmf::kRHC && pid < 0) )
1695 // wgt = fHornCurrentHists.at(1)->Interpolate(E);
1696 // if( (beamType == cmf::kFHC && pid < 0) || (beamType == cmf::kRHC && pid > 0) )
1697 // wgt = fHornCurrentHists.at(3)->Interpolate(E);
1698 // }
1699 //
1700 // LOG_DEBUG("HornCurrentWeight") << "Horn current weight is: " << wgt;
1701 //
1702 // return wgt;
1703 // }
1704 
1705  //----------------------------------------------------------------------------
1706  // Function to apply a normalization weight only to a given horn current
1707  // There is no check on the selection type for the current event, but that
1708  // could be added in the future
1710  cmf::BeamType_t const& beam)
1711  {
1712  // check the beam type is correct for the current event
1713  if(beam != fCurrentMD.BeamType() ) return 1.;
1714 
1715  auto const& sigma = cmf::ParameterUtility::Instance()->ParameterInfo(wgtType).Sigma();
1716 
1717  LOG_DEBUG("HornNormWeight")
1718  << " Sigma Value "
1719  << sigma
1720  << " beam type "
1721  << beam;
1722 
1723  return 1. + sigma * this->CurrentParameterValue(wgtType);
1724  }
1725 
1726  //----------------------------------------------------------------------------
1728  {
1729  // it doesn't make sense to ask for a CC weight if the event is an NC
1730  // or cosmic muon
1733 
1734  auto parVal = this->CurrentParameterValue("MaCCQEReducedWgt");
1735  std::map<float, float> wgtBySigma = fCurrentEvent->MCVals().SystVarWgts(cmf::VarNameToKey("MaCCQE"));
1736 
1737  // mysterious 2018 correction factor (M_A^QE CV shift from 0.99 to 1.04)
1738  const double correctionInSigma = (1.04 - 0.99) / 0.25;
1739  double cv_weight = 1. + correctionInSigma * (wgtBySigma.at(1) - 1);
1740 
1741  double genie_sigma = (parVal > 0) ? parVal * (0.05 / 0.25) : parVal * (0.05 / 0.15);
1742 
1743  if( cv_weight > 0) {
1744  for(auto wgt : wgtBySigma) wgt.second /= cv_weight;
1745 
1746  double reduced_error = 0.05;
1747  double cv_ma = 1.04;
1748  double genie_ma = 0.99;
1749 
1750  double shifted_ma = cv_ma * ( 1 + parVal * reduced_error );
1751  double genie_frac_shift = ( shifted_ma - genie_ma ) / genie_ma;
1752  double genie_error = ( genie_frac_shift > 0 ) ? 0.25 : 0.15;
1753  genie_sigma = genie_frac_shift / genie_error;
1754  }
1755 
1756  return this->CalcLinearInterpWgt(genie_sigma, wgtBySigma);
1757  }
1758 
1759  //----------------------------------------------------------------------------
1761  {
1765  ((!anti && fCurrentEvent->MCVals().val_at(cmf::kTruePDG) > 0) ||
1766  (anti && fCurrentEvent->MCVals().val_at(cmf::kTruePDG) < 0)) );
1767  }
1768 
1769  //----------------------------------------------------------------------------
1771  {
1772  if( this->UseMECWeight(anti) ){
1773 
1774  std::string sign;
1775  if(!anti) sign = "Nu";
1776  else sign = "AntiNu";
1777 
1779  auto parVal = this->CurrentParameterValue("MECShape"+sign+"Wgt");
1780  double wgt = 1.;
1781 
1782  LOG_DEBUG("ShifterAndWeighter")
1783  << "Found MEC event with "
1784  << " \n-- EmpiricalMEC_Weight (nomWgt) : " << fCurrentEvent->MCVals().val_at(cmf::kEmpiricalMEC_Weight)
1785  << " \n-- EmpiricalMECtoGENIEQE_Weight : " << fCurrentEvent->MCVals().val_at(cmf::kEmpiricalMECtoGENIEQE_Weight)
1786  << " \n-- EmpiricalMECToGENIERES_Weight: " << fCurrentEvent->MCVals().val_at(cmf::kEmpiricalMECtoGENIERES_Weight)
1787  << " \n-- parVal for MECShape" << sign << "Wgt: " << parVal
1788  << " \n-- PDG : " << fCurrentEvent->MCVals().val_at(cmf::kTruePDG);
1789 
1790 
1791  auto wgtVar = (parVal < 0) ?
1794 
1795  if(nomWgt > 0)
1796  wgt = 1 + std::abs(parVal) * (wgtVar / nomWgt - 1);
1797 
1798  LOG_DEBUG("ShifterAndWeighter")
1799  << "wgt: " << wgt << " parVal: " << std::abs(parVal) << " wgtVar: " << wgtVar << " nomWgt: " << nomWgt;
1800 
1801  // don't allow negative weights
1802  if(wgt < 0)
1803  wgt = 0;
1804 
1805  // don't allow crazy weights
1806  if(wgt > 10)
1807  wgt = 10;
1808 
1809 
1810  return wgt;
1811  }
1812 
1813  return 1.0;
1814  }
1815 
1816  //----------------------------------------------------------------------------
1818  {
1819  if( this->UseMECWeight(anti) ){
1820 
1821  LOG_DEBUG("ShifterAndWeighter")
1822  << "Found MEC event with "
1823  << " MEC weight " << fCurrentEvent->MCVals().val_at(cmf::kEmpiricalMEC_Weight)
1826  << " PDG " << fCurrentEvent->MCVals().val_at(cmf::kTruePDG);
1827 
1828  std::string sign;
1829  if(!anti) sign = "Nu";
1830  else sign = "AntiNu";
1831 
1832  auto parVal = this->CurrentParameterValue("MECEnuShape"+sign+"Wgt");
1833  double nuE = fCurrentEvent->MCVals().val_at(cmf::kTrueE);
1834  if(nuE < 0) return 1.0;
1835 
1836  // this function is the 1 sigma bound around central value weight 1
1837  static const TF1 rwgtFn("MECRwgtFn", "1/(2.5*x+1)");
1838  double wgt = 1 + parVal * rwgtFn.Eval(nuE);
1839 
1840  // no negative weights
1841  if(wgt < 0)
1842  wgt = 0;
1843 
1844  LOG_DEBUG("ShifterAndWeighter")
1845  << "\tMECEnuShapeWgt " << wgt
1846  << " (Enu = " << nuE
1847  << ", rwgtFn.Eval(nuE) = " << rwgtFn.Eval(nuE)
1848  << ", parval = " << parVal << ")";
1849  return wgt;
1850  }
1851 
1852  return 1.0;
1853  }
1854 
1855  //----------------------------------------------------------------------------
1857  {
1858  if( this->UseMECWeight(anti) ){
1859  std::string sign;
1860  if(!anti) sign = "Nu";
1861  else sign = "AntiNu";
1862 
1863  auto parVal = this->CurrentParameterValue("MECInitNPFrac"+sign+"Wgt");
1864  const double nominalNPfrac = 0.8; // see DocDB 18741
1865  double newNPfrac = nominalNPfrac + (0.1*parVal);
1866  double wgt;
1867 
1868  // only physical between 0 and 1
1869  if(newNPfrac > 1)
1870  newNPfrac = 1;
1871  else if(newNPfrac < 0)
1872  newNPfrac = 0;
1873 
1874  if(fCurrentEvent->MCVals().val_at(cmf::kTrueHitNuc) == 2000000201)
1875  wgt = newNPfrac / nominalNPfrac;
1876  else
1877  wgt = (1 - newNPfrac) / (1 - nominalNPfrac);
1878 
1879  LOG_DEBUG("ShifterAndWeighter")
1880  << "\tMECInitStateNPFracWgt " << wgt
1881  << " (nominalNPfrac = " << nominalNPfrac
1882  << ", newNPfrac = " << newNPfrac
1883  << ", parval = " << parVal << ")";
1884  return wgt;
1885  }
1886 
1887  return 1.0;
1888  }
1889 
1890  //----------------------------------------------------------------------------
1892  {
1893  auto const& sigma = cmf::ParameterUtility::Instance()->ParameterInfo("MECNormWgt").Sigma();
1894 
1896  return 1. + sigma * this->CurrentParameterValue("MECNormWgt");
1897  }
1898 
1899  return 1.0;
1900  }
1901 
1902  //----------------------------------------------------------------------------
1904  {
1905  auto const& sigma = cmf::ParameterUtility::Instance()->ParameterInfo("NCNormWgt").Sigma();
1907  return 1. + sigma * this->CurrentParameterValue("NCNormWgt");
1908 
1909  return 1.0;
1910  }
1911 
1912  //----------------------------------------------------------------------------
1914  {
1915  // this one only applies to the far detector
1916  if(fCurrentMD.detector == cmf::kNEARDET) return 1.;
1917 
1918  auto const& sigma = cmf::ParameterUtility::Instance()->ParameterInfo("NoiseFarSystWgt").Sigma();
1919 
1920  return 1. + sigma * this->CurrentParameterValue("NoiseFarSystWgt");
1921  }
1922 
1923  //----------------------------------------------------------------------------
1924  // For CAFAna version of this systematic, see CAFAna/Systs/NueAcceptSysts.cxx
1926  {
1927  //Only apply weight to FD, nue selected events
1929  !fCurrentMD.IsNuESelected() ) return 1.;
1930 
1931  //If this is the FHC systematic, only apply to FHC, and vice versa
1932  if(beamType != fCurrentMD.BeamType())
1933  return 1.;
1934 
1935  // LOG_VERBATIM("ShifterAndWeighter")
1936  // << "NueAcceptSignalKin2018Weight time!";
1937 
1938  //TODO: Check that it contains a true neutrino (CAF version is if(sr->mc.nnu return; )
1939 
1940  //Check flavors and that it's CC
1944  return 1;
1945 
1946  // Grab the relevant histogram (if you haven't already): code stolen from LoadSABeamSystHists function,
1947  // I hope that's still a valid way to get histograms
1948  if((beamType == cmf::kFHC && !fNueAcceptSignalKin2018FHCWgtHist) ||
1949  (beamType == cmf::kRHC && !fNueAcceptSignalKin2018RHCWgtHist)){
1951  std::string shortName("Cos_FD_KinematicsCorrection_" + cmf::cBeamType_Strings[beamType]);
1952 
1953  // constructor decides if initialized value is a path or an environment variable
1954  cet::search_path sp("FW_SEARCH_PATH");
1955  if ( !sp.find_file(std::string("CAFAna/nue/Ana2018/AcceptSysts/FD_KinematicsCorrection_" + cmf::cBeamType_Strings[beamType] + ".root"), fileName) ) {
1956  throw cet::exception("NueExtrapSystFileError")
1957  << "cannot find Nue Extrap Systematics file CAFAna/nue/Ana2018/AcceptSysts/FD_KinematicsCorrection_"
1958  << cmf::cBeamType_Strings[beamType]
1959  << ".root "
1960  << fileName
1961  << " bail ungracefully\n"
1962  << "(RJN Guess: Did you run setup for this release?)\n"
1963  << __FILE__ << ":" << __LINE__;
1964  }
1965 
1966  TFile fin(fileName.c_str(), "read");
1967  if(fin.IsZombie()) {
1968  LOG_WARNING("ShifterAndWeighter")
1969  << "Warning: couldn't open nue extrapolation systematic file "
1970  << fileName;
1971  }
1972 
1973  if(beamType == cmf::kFHC){
1974  fNueAcceptSignalKin2018FHCWgtHist = (TH1*)fin.Get(shortName.c_str())->Clone();
1975  fNueAcceptSignalKin2018FHCWgtHist->SetDirectory(nullptr);
1976  }
1977  else{
1978  fNueAcceptSignalKin2018RHCWgtHist = (TH1*)fin.Get(shortName.c_str())->Clone();
1979  fNueAcceptSignalKin2018RHCWgtHist->SetDirectory(nullptr);
1980  }
1981 
1982  fin.Close();
1983  }
1984 
1985  // use the selection to find an offset to the energy to get value from the
1986  // histogram
1987  double selEOffset;
1989  selEOffset = 21.;
1991  selEOffset = 2.;
1993  selEOffset = 20.;
1994  else{
1995  selEOffset = -5;
1996  LOG_DEBUG("ShifterAndWeighter")
1997  << "Could not find selection type: " << fCurrentMD.selectionType;
1998  }
1999 
2000  // There are 6 energy bins for each selection type, each energy bin is 0.5 GeV wide, range
2001  // is 1 - 4 GeV in the FD
2003  LOG_DEBUG("ShifterAndWeighter")
2004  << "selBin = "
2005  << selEOffset
2006  << ", anaBin = "
2007  << anaBin;
2008 
2009  double FDSignalWeight = 1;
2010  if(beamType == cmf::kFHC && fNueAcceptSignalKin2018FHCWgtHist->GetBinContent(anaBin) != 0)
2011  FDSignalWeight = fNueAcceptSignalKin2018FHCWgtHist->GetBinContent(anaBin);
2012  else if(beamType == cmf::kRHC && fNueAcceptSignalKin2018RHCWgtHist->GetBinContent(anaBin) != 0)
2013  FDSignalWeight = fNueAcceptSignalKin2018RHCWgtHist->GetBinContent(anaBin);
2014 
2015  //Find number of sigma to shift by
2016  double parVal = 0;
2017  if(beamType == cmf::kFHC)
2018  parVal = this->CurrentParameterValue("NueAcceptSignalKin2018FHCWgt");
2019  else
2020  parVal = this->CurrentParameterValue("NueAcceptSignalKin2018RHCWgt");
2021 
2022  LOG_DEBUG("ShifterAndWeighter")
2023  << "return weight "
2024  << 1 + (FDSignalWeight - 1) * parVal;
2025 
2026  return 1 + (FDSignalWeight - 1) * parVal;
2027 
2028  }
2029 
2030  //----------------------------------------------------------------------------
2031  // For CAFAna version of this systematic, see CAFAna/Systs/NueAcceptSysts.cxx
2033  {
2034 
2035  //Only apply weight to FD, nue selected events
2039 
2040  //If this is the FHC systematic, only apply to FHC, and vice versa
2041  if(beamType != fCurrentMD.BeamType())
2042  return 1.;
2043 
2044  double FDBkgCorr = (beamType == cmf::kFHC) ? 0.015 : 0.041;
2046  FDBkgCorr = (beamType == cmf::kFHC) ? 0.012 : 0.023;
2047 
2048  //TODO: Check that it contains a true neutrino (CAF version is if(sr->mc.nnu return; )
2049 
2050  //Check if neutrino, but not numu->nue signal, otherwise left unaltered
2056  return 1;
2057 
2058  //Find number of sigma to shift by
2059  double parVal = 0;
2060  if(beamType == cmf::kFHC)
2061  parVal = this->CurrentParameterValue("NueAcceptBkg2018FHCWgt");
2062  else
2063  parVal = this->CurrentParameterValue("NueAcceptBkg2018RHCWgt");
2064 
2065  return 1 + FDBkgCorr * parVal;
2066 
2067  }
2068 
2069  //----------------------------------------------------------------------------
2070  // This weight accounts for possible rock events that sneak into the nue
2071  // peripheral sample
2073  {
2076  fCurrentMD.fileType != cmf::kRockFluxSwap )) return 1.;
2077 
2078  double parVal = this->CurrentParameterValue("NueFDRockMuonWgt");
2079 
2080  // the maximum effect is a 100% uncertainty - we can't go below
2081  // a weight of 0, but we could go above a weight of 2
2082  if(parVal < -1.) parVal = -1.;
2083 
2084  return 1. + parVal;
2085  }
2086 
2087  //----------------------------------------------------------------------------
2089  {
2090  double parVal = this->CurrentParameterValue("RadCorrNueWgt");
2091 
2094  return 1.0 + 0.02 * parVal;
2095  else return 1.0;
2096  }
2097 
2098  //----------------------------------------------------------------------------
2100  {
2101  double parVal = this->CurrentParameterValue("RadCorrNueBarWgt");
2102 
2105  return 1.0 + 0.02 * parVal;
2106  else return 1.0;
2107  }
2108 
2109  //----------------------------------------------------------------------------
2111  cmf::SelectionType_t const& sel)
2112  {
2113  // relative normalization is only for the far detector
2114  if(fCurrentMD.detector == cmf::kNEARDET) return 1.;
2115 
2116  // sel is one of the generic selection types, kNCSelection, kNuMuSelection or kNuESelection
2117  bool wantNue = (sel == cmf::kNuESelection ) && (wgtType.find("Nue") != std::string::npos);
2118  bool wantNumu = (sel == cmf::kNuMuSelection) && (wgtType.find("NuMu") != std::string::npos);
2119  bool wantNC = (sel == cmf::kNCSelection) && (wgtType.find("NC") != std::string::npos);
2120 
2121  // check for the correct selection for the current event with this function
2122  // return 1 if we have a numu selected event but want the nue relative normalization
2123  // or if we have a nue selected event but want the numu relative normalization
2124  if (fCurrentMD.IsNuMuSelected() && !wantNumu) return 1.;
2125  else if(fCurrentMD.IsNuESelected() && !wantNue ) return 1.;
2126  else if(fCurrentMD.IsNCSelected() && !wantNC ) return 1.;
2127 
2128  auto const& sigma = cmf::ParameterUtility::Instance()->ParameterInfo(wgtType).Sigma();
2129 
2130  return 1. + sigma * this->CurrentParameterValue(wgtType);
2131  }
2132 
2133  //----------------------------------------------------------------------------
2135  {
2136  auto parVal = this->CurrentParameterValue("RPARESWgt");
2137  if(parVal > 0) return 1. + parVal * (1 - fCurrentEvent->MCVals().val_at(cmf::kRPARES_Weight));
2138 
2139  return 1.;
2140  }
2141 
2142  //----------------------------------------------------------------------------
2144  {
2147 
2148  std::map<float, float> wgtBySigma = fCurrentEvent->MCVals().SystVarWgts(cmf::VarNameToKey("RPACCQEshapeEnh"));
2149 
2150 
2151  auto initWgt = this->CalcLinearInterpWgt(this->CurrentParameterValue("RPACCQEshapeEnhWgt"),
2152  wgtBySigma);
2153  // if(initWgt < 0.01){
2154  // std::string wgtStr("Check initial RPACCQEshapeEnhMap ");
2155  // for(auto const& itr : wgtBySigma) wgtStr += std::to_string(itr.first) + " " + std::to_string(itr.second) + " ";
2156 
2157  // LOG_VERBATIM("SAWRPACCQEshapeEnhWeight")
2158  // << wgtStr;
2159 
2160  // wgtStr = "SAWRPACCQEshapeEnhWeight " + std::to_string(this->CurrentParameterValue("RPACCQEshapeEnhWgt"));
2161  // wgtStr += " " + std::to_string(initWgt) + " ";
2162  // for(auto const& itr : wgtBySigma) wgtStr += std::to_string(itr.first) + " " + std::to_string(itr.second) + " ";
2163 
2164  // LOG_VERBATIM("SAWRPACCQEshapeEnhWeight")
2165  // << wgtStr;
2166  // }
2167 
2168  // problem: for the SA, it was safe to assume that a 0sigma shift meant a weight of 1
2169  // BUT in A2017, the central value for RPACCQESupp/EnhWeights is the value of RPACCQEWeight
2170  return initWgt * fCurrentEvent->MCVals().val_at(cmf::kRPACCQE_Weight);
2171  }
2172 
2173  //----------------------------------------------------------------------------
2175  {
2176  auto parVal = this->CurrentParameterValue("RPACCQEWgt");
2177  if(parVal > 0) return 1. + parVal * (1 - fCurrentEvent->MCVals().val_at(cmf::kRPACCQE_Weight));
2178 
2179  return 1.;
2180  }
2181 
2182  //----------------------------------------------------------------------------
2184  {
2187 
2188  std::map<float, float> wgtBySigma = fCurrentEvent->MCVals().SystVarWgts(cmf::VarNameToKey("RPACCQEshapeSupp"));
2189 
2190  double wgt = this->CalcLinearInterpWgt(this->CurrentParameterValue("RPACCQEshapeSuppWgt"),
2191  wgtBySigma);
2192 
2193  //See comment in RPACCQEshapeEnhWeight()
2195 
2196  // LOG_VERBATIM("ShifterAndWeighter")
2197  // << "RPACCQEshapeSuppWeight = " << wgt
2198  // << ", from (" << wgtBySigma.at(-1.)
2199  // << ", " << wgtBySigma.at( 0.)
2200  // << ", " << wgtBySigma.at( 1.)
2201  // << ")";
2202  return wgt;
2203  }
2204 
2205  //----------------------------------------------------------------------------
2207  {
2208  // rock muon normalization is only for the near detector
2209  if(fCurrentMD.detector == cmf::kFARDET) return 1.;
2210 
2211  auto const& sigma = cmf::ParameterUtility::Instance()->ParameterInfo("RockNormOnFidWgt").Sigma();
2212 
2213  return 1. + sigma * this->CurrentParameterValue("RockNormOnFidWgt");
2214  }
2215 
2216  //-----------------------------------------------------------------------------
2218  {
2219  auto const& sigma = cmf::ParameterUtility::Instance()->ParameterInfo("TauScaleSystWgt").Sigma();
2220 
2223 
2224  double weight = 1. + sigma * this->CurrentParameterValue("TauScaleSystWgt");
2225 
2226  return std::max(0., weight);
2227  }
2228 
2229  return 1.0;
2230  }
2231 
2232  //==============================================================================//
2233  //----------------- SYSTEMATICS USED IN THE NUMU ANALYSIS ----------------------//
2234  //==============================================================================//
2235 
2236  //==============================================================================//
2237  //----------------- SYSTEMATICS USED IN THE NUE ANALYSIS ----------------------//
2238  //==============================================================================//
2239 
2240  //----------------------------------------------------------------------------
2242  {
2243  double parVal = this->CurrentParameterValue("SecondClassCurrWgt");
2244 
2247  return 1.0 + 0.02 * parVal;
2248  else if(fCurrentEvent->MCVals().val_at(cmf::kTrueCCNC) == 1. * simb::kCC &&
2250  return 1.0 - 0.02 * parVal;
2251  else return 1.0;
2252  }
2253 
2254  //----------------------------------------------------------------------------
2256  cmf::SelectionType_t const& selType)
2257  {
2258  //This should ALWAYS be correlated with:
2259  // RPACCQEshapeSuppWeight
2260  // MECInitStateNPFrac
2261  // MECEnuShape
2262 
2263  double tempwgt = 1.;
2264  double parVal = this->CurrentParameterValue(wgtType);
2265 
2266  //All the GENIE style wgts
2267  tempwgt *= this->GENIESystWeight(wgtType);
2268  LOG_DEBUG("ShifterAndWeighter")
2269  << "SumSmallXSec: with GENIE wgts "
2270  << tempwgt;
2271 
2272  //DIS
2273  tempwgt *= 1. + parVal * (1 - fCurrentEvent->MCVals().val_at(cmf::kOtherDIS_Weight));
2274 
2275  // if(fCurrentEvent->MCVals().val_at(cmf::kOtherDIS_Weight) != 1.)
2276  // LOG_VERBATIM("ShifterAndWeighter")
2277  // << "\twith OtherDIS wgts "
2278  // << tempwgt;
2279 
2280  if(selType == cmf::kNuMuSelection){
2281  //RadCorrNue and RadCorrNueBar
2282  //There should be a NuE version of this weight WITHOUT these systematics.
2283  //For joint analyses, USE THE NUE VERSION OF THE SUMSMALLXSEC WEIGHT
2286  tempwgt *= 1 + 0.02 * parVal;
2287  //SimpleSecondClassCurrentsSecSyst
2288  if(fCurrentEvent->MCVals().val_at(cmf::kTruePDG) == 12) tempwgt *= 1 + 0.02 * parVal;
2289  if(fCurrentEvent->MCVals().val_at(cmf::kTruePDG) == -12) tempwgt *= 1 - 0.02 * parVal;
2290  LOG_DEBUG("ShifterAndWeighter")
2291  << "\twith RadCorr wgts "
2292  << tempwgt;
2293  }
2294  } // end if numu selection
2295 
2296  return tempwgt;
2297  }
2298 
2299  //----------------------------------------------------------------------------
2301  {
2303  }
2304 
2305 
2306 
2307 } // cmf
T max(const caf::Proxy< T > &a, T b)
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
TString fin
Definition: Style.C:24
std::string fFileName
cmf::Parameter const & ParameterInfo(std::string const &parName)
Optimized version of OscCalcPMNS.
void OscillationCalculatorToUse(cmf::CalcType_t calcType)
const XML_Char * name
Definition: expat.h:151
MCVarVals const & MCVals() const
Definition: Event.h:80
virtual void SetL(double L)=0
fileName
Definition: plotROC.py:78
enum cmf::interaction_type InteractionType_t
BeamSyst()=default
std::unique_ptr< Event > EventPtr
Definition: Event.h:91
cmf::Location fCurrentLoc
the current point in parameter space
double CalcLinearInterpWgt(double sigma, std::map< float, float > const &wgtsBySigma)
_OscCalcPMNSOpt< double > OscCalcPMNSOpt
float ShiftForParameter(std::map< std::string, uint8_t > const &shifterMap, cmf::ParameterSpaceLoc const &loc, uint8_t const &param)
bool MCValsFilled() const
Definition: Event.h:73
const std::string cSelectionType_Strings[12]
Definition: Constants.h:73
cmf::DetType_t detector
Definition: Structs.h:114
virtual void SetDmsq21(const T &dmsq21)=0
const Var weight
std::set< double > fTrueEnergyBins
the true energy binning used by CAF
void InitShiftsAndWeightsForDetector(cmf::DetType_t const &det, cmf::ParameterSpaceLoc const &detLoc, std::string const &calibFileName, std::vector< WeighterInfo > &weighters, std::map< std::string, uint8_t > &shifters)
double GeniePCASystWeight(size_t index)
const char * p
Definition: xmltok.h:285
std::map< float, float > CalcGeniePCASystWeights(cmf::MetaData const &currentMD, int mode, double energy)
std::map< float, float > CalcBeamSystWeights(cmf::DetType_t const &curDet, cmf::BeamType_t const &beamType, double energy, int pdgOrig)
string fFileName
Definition: nue_pid_effs.C:43
virtual void SetTh13(const T &th13)=0
Adapt the PMNS_Sterile calculator to standard interface.
std::map< float, float > SystVarWgts(uint8_t const &varkey) const
Definition: VarVals.cxx:53
ParameterSpaceLoc FDLoc
Definition: Parameter.h:31
void SetCurrentEvent(cmf::EventPtr const &curEvent, cmf::MetaData const &md)
enum cmf::det_type DetType_t
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
void abs(TH1 *hist)
ParameterSpaceLoc NDLoc
Definition: Parameter.h:30
void FillShifterInfo(std::string const &name, cmf::DetType_t const &det)
enum cmf::sel_type SelectionType_t
double COHScale2018Weight(std::string const &wgtType, cmf::InteractionType_t const &intType)
double MECInitStateNPFracWeight(bool anti)
std::string find_file(std::string const &filename) const
virtual void SetDmsq32(const T &dmsq32)=0
float abs(float number)
Definition: d0nt_math.hpp:39
std::map< std::string, uint8_t > fShiftersToUseFD
std::map< std::string, uint8_t > fShiftersToUseND
double MECEnuShapeWeight(bool anti)
std::vector< WeighterInfo > fWeightersToUseND
std::string DetectorString() const
Definition: Structs.h:149
int EnergyToBin(double const &energy, cmf::MetaData const &md)
void LoadCalibSystRatios(std::string const &CalibFilename)
float CalibSystRatio(std::string const &systName, cmf::MetaData const &currentMD, double energy)
static ParameterUtility * Instance()
std::unique_ptr< TF1 > fNDMuonCatcherEShiftMinus
negative shift in the muon catcher energy
float const & Sigma() const
Definition: Parameter.h:62
static uint8_t VarNameToKey(std::string const &name)
Definition: StaticFuncs.h:394
double FluxPCAWeight(size_t index)
cmf::SelectionType_t selectionType
Definition: Structs.h:116
double NueAcceptSignalKin2018Weight(cmf::BeamType_t const &beamType)
double BirksNormWeight(std::string const &wgtType, cmf::DetType_t const &det)
enum cmf::weight_type WeightType_t
Float_t E
Definition: plot.C:20
DataVarVals const & DataVals() const
Definition: Event.h:79
T get(std::string const &key) const
Definition: ParameterSet.h:231
double energy
Definition: plottest35.C:25
static ShifterAndWeighter * gShifterAndWeighter
double AbsNormWeight(std::string const &wgtType, cmf::SelectionType_t const &sel)
double RelNormWeight(std::string const &wgtType, cmf::SelectionType_t const &sel)
bool IsNuMuSelected() const
Definition: Structs.cxx:219
double const CurrentParameterValue(std::string const &parName)
void SetRho(const float &r)
Definition: HoughPeak.h:94
static ShifterAndWeighter * Instance()
float val_at(uint8_t const &varkey, cmf::MetaData const &md) const
Definition: VarVals.cxx:173
bool DataValsFilled() const
Definition: Event.h:72
float bin[41]
Definition: plottest35.C:14
const ana::Var wgt
double BeamSystWeight(std::string const &wgtType, cmf::DetType_t const &det)
double sigma(TH1F *hist, double percentile)
std::map< std::string, TH1D * > fBeamSystHists
#define LOG_WARNING(category)
void SetFileAndHistNames(std::string const &filename, std::string const &shortName)
std::string ToString() const
Definition: Structs.cxx:133
virtual T P(int flavBefore, int flavAfter, double E)=0
E in GeV; flavors as PDG codes (so, neg==>antinu)
Module to combine a set of results into a single file currently only does one data product type at a ...
Definition: Event.cxx:24
cmf::Event * fCurrentEvent
the current event to weight/shift
double SumSmallXSec2017Weight(std::string const &wgtType, cmf::SelectionType_t const &selType)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
virtual void SetRho(double rho)=0
cmf::FileType_t fileType
Definition: Structs.h:115
static const double A
Definition: Units.h:82
void InitShiftsAndWeightsToUse(cmf::Location const &loc, fhicl::ParameterSet const &parset)
charged current coherent pion
Definition: MCNeutrino.h:139
double HornNormWeight(std::string const &wgtType, cmf::BeamType_t const &beam)
std::string pars("Th23Dmsq32")
cmf::CalcType_t fCalcType
type of oscillation calculator
cmf::MetaData fCurrentMD
the metadata for the current event
float val_at(uint8_t const &varkey) const
Access a Var by name.
Definition: VarVals.cxx:18
int num
Definition: f2_nu.C:119
virtual void SetTh23(const T &th23)=0
enum cmf::beam_type BeamType_t
void LoadGeniePCASystHists(std::string const &shortName)
std::vector< GeniePCASyst > fVecGeniePCASysts
This is each of the Genie PCA weights.
std::unique_ptr< TF1 > fNDMuonCatcherEShiftPlus
positive shift in the muon catcher energy
std::vector< WeighterInfo > fWeightersToUseFD
bool fTrueEOscillationWeight
flag to calculate the oscillation weight at the true E_nu vs bin center
string shortName
THUMBNAIL BLOCK: We need to make a thumbnail for each.
bool systematics
Definition: fnexvscaf.C:31
std::vector< BeamSyst > fVecFluxPrincipals
This is each of the Flux PCA weights.
std::map< std::string, Parameter > ParameterSpaceLoc
Definition: Parameter.h:15
std::string fShortName
float FractionalShift(uint8_t const &param)
std::unique_ptr< TGraph > fNDMuonCatcherMinusGr
Graph for the negative shift in the muon catcher.
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
enum cmf::calc_type CalcType_t
cmf::BeamType_t BeamType() const
Definition: Structs.cxx:178
#define LOG_VERBATIM(category)
double MECShapeWeight(bool anti)
cmf::InteractionType_t interactionType
Definition: Structs.h:117
const std::string cDetType_Strings[5]
Definition: Constants.h:509
const std::string cBeamType_Strings[4]
Definition: Constants.h:29
double Weight(cmf::EventPtr const &curEvent, cmf::MetaData const &md, cmf::WeightType_t wgtType=kAllWgt)
void SetCurrentVals(cmf::Location const &loc)
void LoadCalibrationSystRatios(const std::string &CalibFilename)
std::unique_ptr< TGraph > fNDMuonCatcherPlusGr
Graph for the positive shift in the muon catcher.
osc::IOscCalcAdjustable * fOscCalc
the oscillation calculator
double CalibSystWeight(std::string const &wgtType)
double NueAcceptBkg2018Weight(cmf::BeamType_t const &beamType)
def sign(x)
Definition: canMan.py:197
virtual void SetTh12(const T &th12)=0
virtual void SetdCP(const T &dCP)=0
static CovarianceBinUtility * Instance()
double GENIESystWeight(std::string const &systName)
bool IsNuESelected() const
Definition: Structs.cxx:225
bool IsNCSelected() const
Definition: Structs.cxx:231