ShifterAndWeighter.cxx
Go to the documentation of this file.
1 
2 //
3 // ShifterAndWeighter.cxx
4 // Implementation of singleton to calculate weights and keep track of shifts
5 // for CovarianceMatrixFit fitting
6 //
7 
8 #include "TH2.h"
9 #include "TFile.h"
10 #include "TKey.h"
11 #include <algorithm>
12 
16 
22 
23 #include "OscLib/OscCalcPMNS_NSI.h"
24 #include "OscLib/OscCalcSterile.h"
26 #include "Utilities/func/MathUtil.h"
27 
28 
29 namespace cmf {
30 
32 
33  //----------------------------------------------------------------------------
35  {
36  if(!gShifterAndWeighter) gShifterAndWeighter = new ShifterAndWeighter();
37 
38  return gShifterAndWeighter;
39  }
40 
41  //----------------------------------------------------------------------------
43  : fOscCalc (nullptr)
44  , fCalcType(cmf::kPMNS)
45  {}
46 
47  //---------------------------------------------------------------------------
49 
50  TH1D* ratioHist = nullptr;
51  std::string ratioName = "";
52  auto const& dbs = cmf::SelectionUtility::Instance()->SelectionsToUse();
53  for(auto const& itr : dbs){
54  for(auto const& sel : itr.Selections()){
55 
56  std::string strippedName = systName;
57  if ((systName.substr(systName.length()-2, systName.length()) == "ND") |
58  (systName.substr(systName.length()-2, systName.length()) == "FD"))
59  strippedName = systName.substr(0, systName.length()-2);
60 
61  // TODO: get the biggest excursion of +/- 1 sigma
62  ratioName = strippedName
63  + "Total"
64  + cmf::cBeamType_Strings[itr.BeamType()]
65  + cmf::cDetType_Strings[itr.Detector()]
67  + "Plus1Sigma";
68 
69  ratioHist = (TH1D*)file->Get(ratioName.c_str());
70 
71  if (ratioHist == nullptr)
72  throw cet::exception("HistLookupError")
73  << "Histogram: "
74  << ratioName
75  << " not found in file "
76  << file->GetName()
77  << ". Bailing";
78 
79  for(int i = 1; i < ratioHist->GetNbinsX() + 1; ++i){
80  fileSystVec[ratioName].push_back(ratioHist->GetBinContent(i));
81 
82  }
83  }
84  }
85  }
86 
87  //---------------------------------------------------------------------------
89  std::string filePath;
90  cet::search_path sp("FW_SEARCH_PATH");
91  if(!sp.find_file(calibFileName, filePath)){
92  throw cet::exception("FileLookupError")
93  << "Calibration systematics file: "
94  << calibFileName
95  << " cannot be found in FW_SEARCH_PATH. Bailing.";
96  }
97 
98  calibSystFile = new TFile(filePath.c_str(), "READ");
99 
100  if (calibSystFile->IsZombie()){
101  throw cet::exception("FileOpenError")
102  << "Cannot open file: "
103  << filePath
104  << ". Bailing";
105  }
106 
107  }
108 
109  //---------------------------------------------------------------------------
111 
112  std::string filePath;
113  cet::search_path sp("FW_SEARCH_PATH");
114  if(!sp.find_file(eshiftFileName, filePath)){
115  throw cet::exception("FileLookupError")
116  << "EShift systematics file: "
117  << eshiftFileName
118  << " cannot be found in FW_SEARCH_PATH. Bailing.";
119  }
120 
121  eshiftSystFile = new TFile(filePath.c_str(), "READ");
122 
123  if (eshiftSystFile->IsZombie()){
124  throw cet::exception("FileOpenError")
125  << "Cannot open file: "
126  << filePath
127  << ". Bailing";
128  }
129 
130  }
131 
132 
133  //----------------------------------------------------------------------------
135  cmf::ParameterSpaceLoc const& detLoc,
136  std::vector<WeighterInfo> & weighters,
137  std::map<std::string, uint8_t> & shifters)
138  {
139  // loop over the systematic parameters from the input point and evaluate the
140  // name associated with each one, then add the appropriate function to the
141  // list.
142  weighters.clear();
143  shifters .clear();
144 
145  std::string systName;
146  WeighterInfo weighter;
148  for(auto const& itr : detLoc){
149  if(itr.second.IsOscPar()) continue;
150  systName = itr.first;
151  parType = itr.second.ParType();
152 
153  if(systName.find("Shifter") != std::string::npos){
154 
155  LOG_DEBUG("ShifterAndWeighter")
156  << "Adding shifter for "
157  << systName
158  << " to detector "
159  << det;
160 
161  this->FillShifterInfo(systName, det);
162  } // end if this is a shifter
163  else{
164  weighter.wgtName = systName;
165  // check all the systematic parameter names and add the necessary functions
166  // for weighting if there is a function for that name.
167  if(parType == cmf::kCalibSystPar || parType == cmf::kEShiftSystPar){
168  if (eshiftSystFile==nullptr && parType == cmf::kEShiftSystPar) {
169  this->InitEShiftSystFile();
170  this->LoadRatiosFromFile(eshiftSystFile, systName);
171  }
172  if (calibSystFile==nullptr && parType == cmf::kCalibSystPar){
173  this->InitCalibSystFile();
174  this->LoadRatiosFromFile(calibSystFile, systName);
175  }
176  weighter.wgtFn = std::bind(&ShifterAndWeighter::FileSystWeight, this, systName);
177  }
178  else if(parType == cmf::kCMFNormSystPar)
179  weighter.wgtFn = std::bind(&ShifterAndWeighter::NormSystWeight, this, systName);
180  else if(parType == cmf::kCMFSystPar)
181  weighter.wgtFn = std::bind(&ShifterAndWeighter::CMFSystVarWeight, this, systName);
182  else{
183  LOG_WARNING("ShifterAndWeighter")
184  << "Unable to determine weighting function for "
185  << systName;
186  continue;
187  }
188  LOG_DEBUG("ShifterAndWeighter")
189  << "Adding weighter "
190  << weighter.wgtName
191  << " to detector "
193  << "\nthe current weight type is "
194  << fWeightType;
195 
196  weighters.push_back(weighter);
197 
198  } // the parameter was for a weight
199  } // end loop over systematic parameters
200 
201  }
202 
203  //----------------------------------------------------------------------------
205  fhicl::ParameterSet const& parset)
206  {
207  // this method sets which systematics are to be used in calculating weights
208  // by creating a set of std::functions corresponding to the weighting functions
209  // that are to be turned on
210  // the code is a bit tedious, but it gets the job done
211 
212  calibFileName = parset.get<std::string>("CalibSystFile" , "CovarianceMatrixFit/data/calibration/CalibrationSystematics2020.root");
213  eshiftFileName = parset.get<std::string>("EShiftSystFile" , "CovarianceMatrixFit/data/eshift/EShiftSystematics2020.root");
214  auto calcString = parset.get<std::string>("OscCalculatorType", "PMNS");
215  fTrueEOscillationWeight = parset.get<bool >("TrueEOscCalc" , true);
216 
218  auto weightString = parset.get<std::string>("WeightType", "OscWeightsOnly");
219  if (weightString == "SystAndOscWeights") fWeightType = cmf::kSystAndOscWgt;
220  else if(weightString == "SystWeightsOnly" ) fWeightType = cmf::kSystWgtOnly;
221 
222  cmf::CalcType_t calcType = cmf::kPMNS;
223  if (calcString == "PMNS_NSI" ) calcType = cmf::kNSI;
224  else if(calcString == "PMNS_Sterile") calcType = cmf::kSterile;
225 
226  // setup the oscillation calculator to use
227  this->OscillationCalculatorToUse(calcType);
228 
229  // go ahead and set the current values
230  this->SetCurrentVals(loc);
231 
233  loc.NDLocation(),
237  loc.FDLocation(),
240 
241  /*
242  for(auto const& itr : fShiftersToUseND){
243  LOG_VERBATIM("ShifterAndWeighter")
244  << "ND shifter param: "
245  << itr.first;
246  for(auto const& si : itr.second)
247  LOG_VERBATIM("ShifterAndWeighter")
248  << " associated shifter "
249  << si.fShifterName
250  << " "
251  << si.fDetector;
252  }
253  for(auto const& itr : fShiftersToUseFD){
254  LOG_VERBATIM("ShifterAndWeighter")
255  << "FD shifter param: "
256  << itr.first;
257  for(auto const& si : itr.second)
258  LOG_VERBATIM("ShifterAndWeighter")
259  << " associated shifter "
260  << si.fShifterName
261  << " "
262  << si.fDetector;
263  }
264  */
265 
266  LOG_DEBUG("ShifterAndWeighter")
267  << "Done initialising shifters and weighters";
268 
269  }
270 
271  //----------------------------------------------------------------------------
273  cmf::DetType_t const& det)
274  {
275  uint8_t param;
276  if (name.find("HadE") != std::string::npos) param = cmf::kHad_RecoE;
277  else if(name.find("LepE") != std::string::npos) param = cmf::kLep_RecoE;
278  else if(name.find("NuE" ) != std::string::npos) param = cmf::kNu_RecoE;
279 
280  if (det == cmf::kNEARDET) fShiftersToUseND.emplace(name, param);
281  else if(det == cmf::kFARDET ) fShiftersToUseFD.emplace(name, param);
282  }
283 
284  //----------------------------------------------------------------------------
286  {
287  LOG_VERBATIM("ShifterAndWeighter")
288  << "Current shifter/weighter point:\n"
289  << fCurrentLoc;
290  }
291 
292  //----------------------------------------------------------------------------
293  // Use this version if there are any parameters with different values in the ND
294  // and FD
296  {
298  curLoc.FDLocation());
299  //this->ReportCurrentVals();
300 
301  this->SetOscillationVals();
302  }
303 
304  //----------------------------------------------------------------------------
305  // Only use this version if you are sure the ND and FD parameter values are the
306  // same when the same parameter is present and changing during a fit in both
307  // detectors
308  void ShifterAndWeighter::SetCurrentVals(std::vector<std::string> const& pars,
309  const double* vals)
310  {
311 
312  for(size_t pv = 0; pv < pars.size(); ++pv){
313  fCurrentLoc.SetParameterValue(pars[pv], vals[pv], cmf::kFARDET);
314  fCurrentLoc.SetParameterValue(pars[pv], vals[pv], cmf::kNEARDET);
315  }
316  this->ReportCurrentVals();
317 
318  this->SetOscillationVals();
319  }
320 
321  //----------------------------------------------------------------------------
323  {
324  // get the oscillation parameters from the current point, we haven't set
325  // the detector yet, but the only thing that changes is the
326  // baseline which we change with the SetCurrentEvent method
327  for(auto const& itr : fCurrentLoc.FDLocation()){
328 
330  if (itr.first == "L" ) fOscCalc->SetL (itr.second.Value());
331  else if(itr.first == "Rho" ) fOscCalc->SetRho (itr.second.Value());
332  else if(itr.first == "Dmsq21") fOscCalc->SetDmsq21(itr.second.Value());
333  else if(itr.first == "Dmsq32") fOscCalc->SetDmsq32(itr.second.Value());
334  else if(itr.first == "Th12" ) fOscCalc->SetTh12 (itr.second.Value());
335  else if(itr.first == "Th13" ) fOscCalc->SetTh13 (itr.second.Value());
336  else if(itr.first == "Th23" ) fOscCalc->SetTh23 (itr.second.Value());
337  else if(itr.first == "dCP" ) fOscCalc->SetdCP (itr.second.Value());
338  if(fCalcType == cmf::kNSI){
339  //Eps = NSI epsilon parameters
340  if (itr.first == "Eps_ee" ) dynamic_cast<osc::OscCalcPMNS_NSI*>(fOscCalc)->SetEps_ee (itr.second.Value());
341  else if(itr.first == "Eps_emu" ) dynamic_cast<osc::OscCalcPMNS_NSI*>(fOscCalc)->SetEps_emu (itr.second.Value());
342  else if(itr.first == "Eps_etau" ) dynamic_cast<osc::OscCalcPMNS_NSI*>(fOscCalc)->SetEps_etau (itr.second.Value());
343  else if(itr.first == "Eps_mumu" ) dynamic_cast<osc::OscCalcPMNS_NSI*>(fOscCalc)->SetEps_mumu (itr.second.Value());
344  else if(itr.first == "Eps_mutau" ) dynamic_cast<osc::OscCalcPMNS_NSI*>(fOscCalc)->SetEps_mutau (itr.second.Value());
345  else if(itr.first == "Eps_tautau" ) dynamic_cast<osc::OscCalcPMNS_NSI*>(fOscCalc)->SetEps_tautau (itr.second.Value());
346  else if(itr.first == "Delta_emu" ) dynamic_cast<osc::OscCalcPMNS_NSI*>(fOscCalc)->SetDelta_emu (itr.second.Value());
347  else if(itr.first == "Delta_etau" ) dynamic_cast<osc::OscCalcPMNS_NSI*>(fOscCalc)->SetDelta_etau (itr.second.Value());
348  else if(itr.first == "Delta_mutau") dynamic_cast<osc::OscCalcPMNS_NSI*>(fOscCalc)->SetDelta_mutau(itr.second.Value());
349  } // end if NSI oscillations
350  } // end if PMNS or NSI oscillations
351  else if(fCalcType == cmf::kSterile){
352  if (itr.first == "Th12" ) dynamic_cast<osc::OscCalcSterile*>(fOscCalc)->SetAngle(1, 2, itr.second.Value());
353  else if(itr.first == "Th23" ) dynamic_cast<osc::OscCalcSterile*>(fOscCalc)->SetAngle(2, 3, itr.second.Value());
354  else if(itr.first == "Th24" ) dynamic_cast<osc::OscCalcSterile*>(fOscCalc)->SetAngle(2, 4, itr.second.Value());
355  else if(itr.first == "Th34" ) dynamic_cast<osc::OscCalcSterile*>(fOscCalc)->SetAngle(3, 4, itr.second.Value());
356  else if(itr.first == "dCP" ) dynamic_cast<osc::OscCalcSterile*>(fOscCalc)->SetDelta(1, 3, itr.second.Value());
357  else if(itr.first == "d24" ) dynamic_cast<osc::OscCalcSterile*>(fOscCalc)->SetDelta(2, 4, itr.second.Value());
358  else if(itr.first == "Dmsq21") dynamic_cast<osc::OscCalcSterile*>(fOscCalc)->SetDm (2, itr.second.Value());
359  else if(itr.first == "Dmsq32") dynamic_cast<osc::OscCalcSterile*>(fOscCalc)->SetDm (3, itr.second.Value());
360  else if(itr.first == "Dmsq41") dynamic_cast<osc::OscCalcSterile*>(fOscCalc)->SetDm (4, itr.second.Value());
361  else if(itr.first == "L" ) dynamic_cast<osc::OscCalcSterile*>(fOscCalc)->SetL ( itr.second.Value());
362  else if(itr.first == "Rho" ) dynamic_cast<osc::OscCalcSterile*>(fOscCalc)->SetRho ( itr.second.Value());
363  } // end if sterile oscillations
364  } // end loop over FD location
365  }
366 
367  //----------------------------------------------------------------------------
369  cmf::MetaData const& md)
370  {
371  fCurrentEvent = curEvent.get();
372  fCurrentMD = md;
373 
374  LOG_DEBUG("SetCurrentEvent")
375  << " Metadata is "
376  << md.ToString();
377 
378  // for(auto const& itr : fCurrentEvent->MCVals().SystematicShifts()){
379  // LOG_VERBATIM("EventListManipulator")
380  // << "ShifterAndWeighter::SetCurrentEvent() systematic: "
381  // << itr;
382  // }
383 
384  // check the validity of the mcVals before trying to get weights, if it is
385  // a nullptr, we have a problem.
386  if( !fCurrentEvent->MCValsFilled() && md.isMC)
387  throw cet::exception("ShifterAndWeighter")
388  << "Event does not have a valid MCVarVals unique_ptr, bail";
389 
390  // the data values better always be valid
391  if( !fCurrentEvent->DataValsFilled() )
392  throw cet::exception("ShifterAndWeighter")
393  << "Event does not have a valid DataVarVals unique_ptr, bail";
394 
395  // reset the oscillation baseline based on the detector
396  fOscCalc->SetL(this->CurrentParameterValue("L"));
397 
398  }
399 
400  //----------------------------------------------------------------------------
402  {
403  return (fCurrentMD.detector == cmf::kNEARDET) ?
405  }
406 
407  //----------------------------------------------------------------------------
408  float ShifterAndWeighter::ShiftForParameter(std::map<std::string, uint8_t> const& shifterMap,
410  uint8_t const& param)
411  {
412 
413  float totalShift = 1.;
414  float shiftScale;
415  float curShift;
416 
417  for(auto const& itr : shifterMap){
418 
419  if(itr.second != param) continue;
420 
421  auto const& sigma = loc.find(itr.first)->second.Sigma();
422 
423  curShift = this->CurrentParameterValue(itr.first);
424 
425  shiftScale = 1. + sigma * curShift;
426 
427  // LOG_VERBATIM("ShifterAndWeighterTotalShift")
428  // << " associated shifter "
429  // << infoItr.fShifterName
430  // << " "
431  // << infoItr.fDetector
432  // << " "
433  // << fCurrentMD.detector
434  // << " "
435  // << sigma
436  // << " "
437  // << this->CurrentParameterValue(infoItr.fShifterName)
438  // << " "
439  // << totalShift;
440  totalShift *= shiftScale;
441  }
442 
443  return totalShift;
444  }
445 
446  //----------------------------------------------------------------------------
448  cmf::MetaData const& md)
449  {
450  // get ready to use the passed in event
451  this->SetCurrentEvent(curEvent, md);
452 
453  // set the weight to be the fake weight value, which is 1 for
454  // data.
455  double wgt = curEvent->DataVals().val_at(cmf::kFake_Weight, md);
456 
457 
458  // if the mcVals is null or the metadata claims it isn't MC,
459  // return the fake weight for the event. If the event is data
460  // then the value is always 1, if it is fake data, then the
461  // value is the proper weighting
462  if(!curEvent->MCValsFilled() ||
463  !md.isMC ||
464  md.IsCosmicMuon() )
465  return wgt;
466 
467  //if(md.selectionType == cmf::kNuESelectionPeripheral)
468  // LOG_VERBATIM("ShifterAndWeighter")
469  // << "weight from systematics is "
470  // << wgt;
471 
472  // check if we only want weights from the systematic parameters or only
473  // from the oscillation parameters
475  wgt *= this->TotalWeightFromWeighters();
476  }
477  else if(fWeightType == cmf::kOscWgtOnly){
478  wgt *= this->OscillationWeight();
479  }
480  else
481  wgt *= this->OscillationWeight() * this->TotalWeightFromWeighters();
482 
483  //if(md.selectionType == cmf::kNuESelectionPeripheral)
484  // LOG_VERBATIM("ShifterAndWeighter")
485  // << "weight from oscillations is "
486  // << this->OscillationWeight();
487 
488  LOG_DEBUG("ShifterAndWeighter")
489  << "Oscillation weight is "
490  << wgt
491  << " energy: "
493  << " pdg: "
495  << " orig pdg: "
497  << " detector: "
498  << fCurrentMD.detector;
499 
500  //if(md.selectionType == cmf::kNuESelectionPeripheral)
501  // LOG_VERBATIM("ShifterAndWeighter")
502  // << "total weight is "
503  // << wgt;
504 
505  return wgt;
506  }
507 
508  //----------------------------------------------------------------------------
510  {
511  double wgt = 1.;
512 
514 
515  //std::string ss;
516  for(auto itr : weighters){
517  // ss += "\t" + itr.wgtName +
518  // " " + std::to_string(itr.wgtFn()) +
519  // " " + std::to_string(wgt) +
520  // " " + std::to_string(wgt * itr.wgtFn()) +
521  // "\n";
522 
523  wgt *= itr.wgtFn();
524  }
525  // LOG_VERBATIM("SAWWeighterInfo")
526  // << ss;
527 
528  return wgt;
529  }
530 
531  //----------------------------------------------------------------------------
533  {
534  // nothing to change
535  if(fOscCalc != nullptr && fCalcType == calcType) return;
536  else{
537  LOG_WARNING("ShifterAndWeighter")
538  << "resetting oscillation calculator to type "
539  << calcType
540  << " from "
541  << fCalcType
542  << " are you sure you meant to change types in this job?";
543 
544  delete fOscCalc;
545  fOscCalc = nullptr;
546  }
547 
548  fCalcType = calcType;
549 
551  else if(fCalcType == cmf::kNSI ) fOscCalc = new osc::OscCalcPMNS_NSI();
552  else if(fCalcType == cmf::kSterile){
554  dynamic_cast<osc::OscCalcSterile*>(fOscCalc)->SetNFlavors(4);
555  }
556  }
557 
558  //----------------------------------------------------------------------------
560  {
561  // don't bother doing anything if we are looking at the near detector
562  // and standard oscillations, it is just a waste of time.
563  if(fCalcType == cmf::kPMNS &&
564  fCurrentMD.detector == cmf::kNEARDET) return 1.;
565 
567 
568  return this->OscillationWeightBinCenter();
569  }
570  //----------------------------------------------------------------------------
572  {
573  double E = fCurrentEvent->MCVals().val_at(cmf::kTrueE);
575  int pdgOrig = fCurrentEvent->MCVals().val_at(cmf::kTruePDGOrig);
576 
577  // only attempt to weight neutrinos for oscillations
578  if(std::abs(pdg) != 12 &&
579  std::abs(pdg) != 14 &&
580  std::abs(pdg) != 16) return 1.;
581 
582  return fOscCalc->P(pdgOrig, pdg, E);
583 
584  }
585 
586  //----------------------------------------------------------------------------
588  {
589  if(fTrueEnergyBins.empty()){
590  const int numTrueEnergyBins = 100;
591 
592  // How many edges to generate. Allow room for 0-Emin bin
593  const double N = numTrueEnergyBins-1;
594  const double A = N * 0.5; // 500 MeV: there's really no events below there
595 
596  fTrueEnergyBins.insert(0);
597  for(int i = 1; i < N - 1; ++i) fTrueEnergyBins.insert(A/i);
598  fTrueEnergyBins.insert(120.);
599  }
600 
601  double E = fCurrentEvent->MCVals().val_at(cmf::kTrueE);
603  int pdgOrig = fCurrentEvent->MCVals().val_at(cmf::kTruePDGOrig);
604 
605  // get the relevant bin for the current neutrino, then get the
606  // energy at the center of that bin
607  double eUpper = *(fTrueEnergyBins.rbegin());
608  double eLower = *(fTrueEnergyBins.begin());
609 
610  // find the bin edge that goes after the value of E
611  auto itr = fTrueEnergyBins.upper_bound(E);
612  if(itr != fTrueEnergyBins.end()) eUpper = *itr;
613  else --itr;
614 
615  // now find the edge before by just decrementing the iterator
616  --itr;
617  if(itr != fTrueEnergyBins.end()) eLower = *itr;
618 
619  // get the midpoint between the upper and lower values,
620  // ie take their mean
621  E = 0.5 * (eUpper + eLower);
622 
623  LOG_DEBUG("ShifterAndWeighter")
624  << "true energy: "
626  << " central bin energy "
627  << E;
628 
629  // only attempt to weight neutrinos for oscillations
630  if(std::abs(pdg) != 12 &&
631  std::abs(pdg) != 14 &&
632  std::abs(pdg) != 16) return 1.;
633 
634  return fOscCalc->P(pdgOrig, pdg, E);
635 
636  }
637 
638  //----------------------------------------------------------------------------
640  std::map<float, float> const& wgtsBySigma)
641  {
642  auto mapEnd = wgtsBySigma.end();
643 
644  // no way this is a valid number and is obviously stupid
645  // not using a std::numeric_limits value to avoid rounding to 0
646  double badVal = -99999999999999.;
647  double minus2Val = (wgtsBySigma.find(-2.) != mapEnd) ? static_cast<double>(wgtsBySigma.find(-2.)->second) : badVal;
648  double minus1Val = (wgtsBySigma.find(-1.) != mapEnd) ? static_cast<double>(wgtsBySigma.find(-1.)->second) : badVal;
649  // double zeroVal = 1.;
650  double zeroVal = (wgtsBySigma.find( 0.) != mapEnd) ? static_cast<double>(wgtsBySigma.find( 0.)->second) : 1.;
651  double plus1Val = (wgtsBySigma.find( 1.) != mapEnd) ? static_cast<double>(wgtsBySigma.find( 1.)->second) : badVal;
652  double plus2Val = (wgtsBySigma.find( 2.) != mapEnd) ? static_cast<double>(wgtsBySigma.find( 2.)->second) : badVal;
653  double slope = 0.;
654  double intercept = 0.;
655 
656  //Debugging RPA weights
657  LOG_DEBUG("ShifterAndWeighter")
658  << "\tCalcLinearInterpWgt: "
659  << " -2:" << minus2Val
660  << ", -1:" << minus1Val
661  << ", 0:" << zeroVal
662  << ", 1:" << plus1Val
663  << ", 2:" << plus2Val;
664 
665  /*
666  if(wgtsBySigma.find(-1.) != mapEnd)
667  LOG_VERBATIM("ShifterAndWeighter")
668  << "\t\t 0 val from "
669  << wgtsBySigma.find(0.)->second;
670  */
671 
672  // The following code assumes that at least the -1 sigma, +1 sigma
673  // values are in the map. If they are not, we have a problem
674  if(minus1Val == badVal ||
675  plus1Val == badVal){
676  LOG_WARNING("ShifterAndWeighterCalcLinearInterpWgt")
677  << "No values associated with either -1 sigma or +1 sigma "
678  << "so we cannot calculate a weight, return 1.";
679  return 1.;
680  }
681 
682  if(sigma <= -1){
683  // It is possible there is no -2 sigma value from the map,
684  // in that case just interpolate from 0 to -1 sigma
685  if(minus2Val != badVal){
686  slope = (minus1Val - minus2Val);
687  intercept = minus2Val - slope * (-2.);
688  }
689  else{
690  slope = (zeroVal - minus1Val);
691  intercept = minus1Val - slope * (-1.);
692  }
693  }
694  else if(sigma > -1. && sigma <= 0.){
695  slope = (zeroVal - minus1Val);
696  intercept = minus1Val - slope * (-1.);
697  }
698  else if(sigma > 0. && sigma <= 1.){
699  slope = (plus1Val - zeroVal);
700  intercept = plus1Val - slope * 1.;
701  }
702  else if(sigma > 1.){
703  // It is possible there is no +2 sigma value from the map,
704  // in that case just interpolate from 0 to +1 sigma
705  if(plus2Val != badVal){
706  slope = (plus2Val - plus1Val);
707  intercept = plus2Val - slope * 2.;
708  }
709  else{
710  slope = (plus1Val - zeroVal);
711  intercept = plus1Val - slope * 1.;
712  }
713  }
714 
715  double wgt = sigma * slope + intercept;
716 
717  // Don't let weights be negative
718  if(wgt < 0) wgt = 0.;
719 
720  return wgt;
721  }
722 
723  //-------------------------------------------------------------------------
725  // we have a separate handling for normalisation systematics because
726  // there's a bunch of numbers it's just not worth storing in the event
727  // list
728  double oneSigmaWgt = 1.0;
729 
730  double normPOT = 0.0055;
731  double normNDMass = 0.003;
732  double normFDMass = 0.003;
733  double normML = 0.008;
734  double normNDNumuContext = fCurrentMD.BeamType() == cmf::kFHC ? 0.0021 : 0.0048;
735  double normNDNueContext = fCurrentMD.BeamType() == cmf::kFHC ? 0.0041 : 0.0053;
736  double normNDNCContext = fCurrentMD.BeamType() == cmf::kFHC ? 0.01 : 0.00;
737 
738  bool isND = fCurrentMD.detector == cmf::kNEARDET ? true : false;
739  bool isFD = fCurrentMD.detector == cmf::kFARDET ? true : false;
740  bool isFHC = fCurrentMD.BeamType() == cmf::kFHC ? true : false;
741  bool isRHC = fCurrentMD.BeamType() == cmf::kRHC ? true : false;
742  bool isNuMu = false;
743  bool isNuE = false;
744  bool isNC = false;
745 
751  isNuMu = true;
755  isNuE = true;
757  isNC = true;
758 
759  if (systName == "NormalizationPOT") oneSigmaWgt += normPOT;
760  else if (systName == "NormalizationNDMass" && isND) oneSigmaWgt += normNDMass;
761  else if (systName == "NormalizationFDMass" && isFD) oneSigmaWgt += normFDMass;
762  else if (systName == "NormalizationML" && isFD) oneSigmaWgt += normML;
763  else if (systName == "NormalizationFHCNuMuContext" && isND && isNuMu && isFHC) oneSigmaWgt += normNDNumuContext;
764  else if (systName == "NormalizationRHCNuMuContext" && isND && isNuMu && isRHC) oneSigmaWgt += normNDNumuContext;
765  else if (systName == "NormalizationFHCNuEContext" && isND && isNuE && isFHC) oneSigmaWgt += normNDNueContext;
766  else if (systName == "NormalizationRHCNuEContext" && isND && isNuE && isRHC) oneSigmaWgt += normNDNueContext;
767  else if (systName == "NormalizationFHCNCContext" && isND && isNC && isFHC) oneSigmaWgt += normNDNCContext;
768  else if (systName == "NormalizationRHCNCContext" && isND && isNC && isRHC) oneSigmaWgt += normNDNCContext;
769 
770  float sigShift = this->CurrentParameterValue(systName);
771  return 1. - (1. - oneSigmaWgt) * sigShift;
772 
773  }
774 
775  //-------------------------------------------------------------------------
777  {
778  bool isWrongDet = false;
779  if (systName.find("ND") != std::string::npos && fCurrentMD.detector == cmf::kFARDET) isWrongDet = true;
780  if (systName.find("FD") != std::string::npos && fCurrentMD.detector == cmf::kNEARDET) isWrongDet = true;
781 
782  std::string strippedName = systName;
783  if ((systName.substr(systName.length()-2, systName.length()) == "ND") |
784  (systName.substr(systName.length()-2, systName.length()) == "FD"))
785  strippedName = systName.substr(0, systName.length()-2);
786 
787  std::string ratioName = strippedName
788  + "Total"
792  + "Plus1Sigma";
793 
796 
797  float oneSigmaWgt = isWrongDet ? 1 : fileSystVec[ratioName][bin];
798  float sigShift = this->CurrentParameterValue(systName);
799 
800  return 1. - (1. - oneSigmaWgt) * sigShift;
801  }
802 
803  //----------------------------------------------------------------------------
805  {
806  // return 1 if this is a cosmic muon
808 
809  // if this weight has NC or CC in the name and this interaction type
810  // is not that, return 1
811  if((systName.find("CC") != std::string::npos && fCurrentMD.interactionType == cmf::kNC) ||
812  (systName.find("NC") == std::string::npos && fCurrentMD.interactionType == cmf::kNC) ) return 1.;
813 
814  // may need to come up with default for cosmic muons, however I think the fact that
815  // the default weights for every GENIE variable is a collection of 1's should
816  // mean this is OK
817 
818  std::map<float, float> wgtBySigma = fCurrentEvent->MCVals().SystVarWgts(cmf::VarNameToKey(systName));
819 
820  //for(auto itr : wgtBySigma)
821  // LOG_VERBATIM("CMFSystVarWeight")
822  // << systName
823  // << " weights "
824  // << itr.first
825  // << " "
826  // << itr.second;
827 
828  return this->CalcLinearInterpWgt(this->CurrentParameterValue(systName),
829  wgtBySigma);
830  }
831 
832  //----------------------------------------------------------------------------
834  {
836  }
837 
838 } // cmf
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
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
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)
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:93
cmf::DetType_t detector
Definition: Structs.h:114
virtual void SetDmsq21(const T &dmsq21)=0
cmf::WeightType_t fWeightType
which type of weights to use - oscillation, systematic or both
std::set< double > fTrueEnergyBins
the true energy binning used by CAF
static SelectionUtility * Instance()
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:47
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
double Weight(cmf::EventPtr const &curEvent, cmf::MetaData const &md)
void FillShifterInfo(std::string const &name, cmf::DetType_t const &det)
double NormSystWeight(std::string const &systName)
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
enum cmf::par_type ParType_t
std::map< std::string, uint8_t > fShiftersToUseFD
std::map< std::string, uint8_t > fShiftersToUseND
std::vector< WeighterInfo > fWeightersToUseND
double CMFSystVarWeight(std::string const &systName)
int EnergyToBin(double const &energy, cmf::MetaData const &md)
static uint8_t VarNameToKey(std::string const &name)
Definition: StaticFuncs.h:381
cmf::SelectionType_t selectionType
Definition: Structs.h:116
cmf::ParameterSpaceLoc const & NDLocation() const
Definition: Parameter.h:160
Float_t E
Definition: plot.C:20
DataVarVals const & DataVals() const
Definition: Event.h:79
void InitShiftsAndWeightsForDetector(cmf::DetType_t const &det, cmf::ParameterSpaceLoc const &detLoc, std::vector< WeighterInfo > &weighters, std::map< std::string, uint8_t > &shifters)
T get(std::string const &key) const
Definition: ParameterSet.h:231
cmf::ParameterSpaceLoc const & FDLocation() const
Definition: Parameter.h:161
double energy
Definition: plottest35.C:25
static ShifterAndWeighter * gShifterAndWeighter
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:161
bool DataValsFilled() const
Definition: Event.h:72
A re-optimized version of OscCalcPMNSOpt.
float bin[41]
Definition: plottest35.C:14
const ana::Var wgt
static bool isFHC
double sigma(TH1F *hist, double percentile)
#define LOG_WARNING(category)
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
std::map< std::string, std::vector< float > > fileSystVec
virtual void SetRho(double rho)=0
static const double A
Definition: Units.h:82
void InitShiftsAndWeightsToUse(cmf::Location const &loc, fhicl::ParameterSet const &parset)
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
virtual void SetTh23(const T &th23)=0
std::vector< WeighterInfo > fWeightersToUseFD
bool fTrueEOscillationWeight
flag to calculate the oscillation weight at the true E_nu vs bin center
TFile * file
Definition: cellShifts.C:17
const DetBeamSelSet & SelectionsToUse()
std::map< std::string, Parameter > ParameterSpaceLoc
Definition: Parameter.h:19
float FractionalShift(uint8_t const &param)
enum cmf::calc_type CalcType_t
void LoadRatiosFromFile(TFile *file, std::string systName)
cmf::BeamType_t BeamType() const
Definition: Structs.cxx:178
#define LOG_VERBATIM(category)
cmf::InteractionType_t interactionType
Definition: Structs.h:117
const std::string cDetType_Strings[5]
Definition: Constants.h:29
const std::string cBeamType_Strings[4]
Definition: Constants.h:51
double FileSystWeight(std::string const &systName)
void SetCurrentVals(cmf::Location const &loc)
bool IsCosmicMuon() const
Definition: Structs.h:139
osc::IOscCalcAdjustable * fOscCalc
the oscillation calculator
void SetParameterValue(std::string const &parName, double const &value, cmf::DetType_t const &det)
Definition: Parameter.cxx:236
virtual void SetTh12(const T &th12)=0
virtual void SetdCP(const T &dCP)=0
static CovarianceBinUtility * Instance()
enum BeamMode string