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 
13 #include "art_root_io/TFileDirectory.h"
16 
22 
23 #include "OscLib/OscCalcPMNS_NSI.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  {
137  std::string systName;
138  WeighterInfo weighter;
140 
141  weighter.detector = det;
142  for(auto const& itr : detLoc.Parameters()){
143  if(itr.IsOscPar()) continue;
144  systName = itr.Name();
145  parType = itr.ParType();
146 
147  if(systName.find("Shifter") != std::string::npos){
148 
149  MF_LOG_DEBUG("ShifterAndWeighter")
150  << "Adding shifter for "
151  << systName
152  << " to detector "
153  << det;
154 
155  this->FillShifterInfo(systName, det);
156  } // end if this is a shifter
157  else{
158  weighter.wgtName = systName;
159  // check all the systematic parameter names and add the necessary functions
160  // for weighting if there is a function for that name.
161  if(parType == cmf::kCalibSystPar || parType == cmf::kEShiftSystPar){
162  if (eshiftSystFile==nullptr && parType == cmf::kEShiftSystPar) {
163  this->InitEShiftSystFile();
164  this->LoadRatiosFromFile(eshiftSystFile, systName);
165  }
166  if (calibSystFile==nullptr && parType == cmf::kCalibSystPar){
167  this->InitCalibSystFile();
168  this->LoadRatiosFromFile(calibSystFile, systName);
169  }
170  weighter.wgtFn = std::bind(&ShifterAndWeighter::FileSystWeight, this, systName);
171  }
172  else if(parType == cmf::kCMFNormSystPar)
173  weighter.wgtFn = std::bind(&ShifterAndWeighter::NormSystWeight, this, systName);
174  else if(parType == cmf::kCMFSystPar)
175  weighter.wgtFn = std::bind(&ShifterAndWeighter::CMFSystVarWeight, this, systName);
176  else if (parType == cmf::kCMFWeightPar)
177  weighter.wgtFn = std::bind(&ShifterAndWeighter::XSecCVPPFXWeight, this, systName);
178  else{
179  MF_LOG_WARNING("ShifterAndWeighter")
180  << "Unable to determine weighting function for "
181  << systName;
182  continue;
183  }
184  MF_LOG_DEBUG("ShifterAndWeighter")
185  << "Adding weighter "
186  << weighter.wgtName
187  << " to detector "
189  << "\nthe current weight type is "
190  << fWeightType;
191 
192  fWeightersToUse.emplace_back(weighter);
193 
194  } // the parameter was for a weight
195  } // end loop over systematic parameters
196 
197  }
198 
199  //----------------------------------------------------------------------------
201  fhicl::ParameterSet const& parset)
202  {
203  // this method sets which systematics are to be used in calculating weights
204  // by creating a set of std::functions corresponding to the weighting functions
205  // that are to be turned on
206  // the code is a bit tedious, but it gets the job done
207 
208  // clear shiftts & weights to use first so we don't double count anything
209  fWeightersToUse.clear();
210  fShiftersToUse.clear();
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  else if(weightString == "NoWeights" ) fWeightType = cmf::kNoWgt;
222 
223  cmf::CalcType_t calcType = cmf::kPMNS;
224  if (calcString == "PMNS_NSI" ) calcType = cmf::kNSI;
225  else if(calcString == "PMNS_Sterile") calcType = cmf::kSterile;
226 
227  // setup the oscillation calculator to use
228  this->OscillationCalculatorToUse(calcType);
229 
230  // go ahead and set the current values
231  this->SetCurrentVals(loc);
232 
233  // TODO make this more elegant
234  // Always set the nova ND and FD, check if we want the MINOS
235  // detectors too
236  for(auto const& detItr : cmf::SelectionUtility::Instance()->DetectorsToUse()){
237  this->InitShiftsAndWeightsForDetector(detItr,
238  loc.DetectorLocation(detItr));
239  }
240 
241  /*
242  for(auto const& itr : fShiftersToUse){
243  MF_LOG_VERBATIM("ShifterAndWeighter")
244  << cmf::cDetType_Strings[itr.detector]
245  << " shifter param: "
246  << itr.shiftName;
247  }
248  */
249 
250  MF_LOG_DEBUG("ShifterAndWeighter")
251  << "Done initialising shifters and weighters";
252 
253  }
254 
255  //----------------------------------------------------------------------------
257  cmf::DetType_t const& det)
258  {
259  uint8_t param;
260  if (name.find("HadE") != std::string::npos) param = cmf::kHad_RecoE;
261  else if(name.find("LepE") != std::string::npos) param = cmf::kLep_RecoE;
262  else if(name.find("NuE" ) != std::string::npos) param = cmf::kCMF_RecoE;
263 
264  fShiftersToUse.emplace_back(det, name, param);
265  }
266 
267  //----------------------------------------------------------------------------
269  {
270  MF_LOG_VERBATIM("ShifterAndWeighter")
271  << "Current shifter/weighter point:\n"
272  << fCurrentLoc;
273  }
274 
275  //----------------------------------------------------------------------------
276  // Use this version if there are any parameters with different values in the ND
277  // and FD
279  {
280  MF_LOG_DEBUG("SAW")
281  << "setting current location in SAW to"
282  << curLoc;
283 
285  //this->ReportCurrentVals();
286 
287  this->SetOscillationVals();
288  }
289 
290  //----------------------------------------------------------------------------
291  // Only use this version if you are sure the ND and FD parameter values are the
292  // same when the same parameter is present and changing during a fit in both
293  // detectors
294  void ShifterAndWeighter::SetCurrentVals(std::vector<std::string> const& pars,
295  const double* vals)
296  {
297 
298  for(auto const& det : cmf::SelectionUtility::Instance()->DetectorsToUse()){
299  for(size_t pv = 0; pv < pars.size(); ++pv){
300  fCurrentLoc.SetParameterValue(pars[pv], vals[pv], det);
301  }
302  }
303  //this->ReportCurrentVals();
304 
305  this->SetOscillationVals();
306  }
307 
308  //----------------------------------------------------------------------------
310  {
311  // get the oscillation parameters from the current point, we haven't set
312  // the detector yet, but the only thing that changes is the
313  // baseline which we change with the SetCurrentEvent method
314  for(auto const& itr : fCurrentLoc.OscillationParameters()){
315 
317  if ((cmf::OscParm_t)itr.Key() == cmf::kL ) fOscCalc->SetL (itr.Value());
318  else if((cmf::OscParm_t)itr.Key() == cmf::kRho ) fOscCalc->SetRho (itr.Value());
319  else if((cmf::OscParm_t)itr.Key() == cmf::kDmsq21) fOscCalc->SetDmsq21(itr.Value());
320  else if((cmf::OscParm_t)itr.Key() == cmf::kDmsq32) fOscCalc->SetDmsq32(itr.Value());
321  else if((cmf::OscParm_t)itr.Key() == cmf::kTh12 ) fOscCalc->SetTh12 (itr.Value());
322  else if((cmf::OscParm_t)itr.Key() == cmf::kTh13 ) fOscCalc->SetTh13 (itr.Value());
323  else if((cmf::OscParm_t)itr.Key() == cmf::kTh23 ) fOscCalc->SetTh23 (itr.Value());
324  else if((cmf::OscParm_t)itr.Key() == cmf::kdCP ) fOscCalc->SetdCP (itr.Value());
325  if(fCalcType == cmf::kNSI){
326  //Eps = NSI epsilon parameters
327  if ((cmf::OscParm_t)itr.Key() == cmf::kEps_ee ) dynamic_cast<osc::OscCalcPMNS_NSI*>(fOscCalc)->SetEps_ee (itr.Value());
328  else if((cmf::OscParm_t)itr.Key() == cmf::kEps_emu ) dynamic_cast<osc::OscCalcPMNS_NSI*>(fOscCalc)->SetEps_emu (itr.Value());
329  else if((cmf::OscParm_t)itr.Key() == cmf::kEps_etau ) dynamic_cast<osc::OscCalcPMNS_NSI*>(fOscCalc)->SetEps_etau (itr.Value());
330  else if((cmf::OscParm_t)itr.Key() == cmf::kEps_mumu ) dynamic_cast<osc::OscCalcPMNS_NSI*>(fOscCalc)->SetEps_mumu (itr.Value());
331  else if((cmf::OscParm_t)itr.Key() == cmf::kEps_mutau ) dynamic_cast<osc::OscCalcPMNS_NSI*>(fOscCalc)->SetEps_mutau (itr.Value());
332  else if((cmf::OscParm_t)itr.Key() == cmf::kEps_tautau ) dynamic_cast<osc::OscCalcPMNS_NSI*>(fOscCalc)->SetEps_tautau (itr.Value());
333  else if((cmf::OscParm_t)itr.Key() == cmf::kDelta_emu ) dynamic_cast<osc::OscCalcPMNS_NSI*>(fOscCalc)->SetDelta_emu (itr.Value());
334  else if((cmf::OscParm_t)itr.Key() == cmf::kDelta_etau ) dynamic_cast<osc::OscCalcPMNS_NSI*>(fOscCalc)->SetDelta_etau (itr.Value());
335  else if((cmf::OscParm_t)itr.Key() == cmf::kDelta_mutau) dynamic_cast<osc::OscCalcPMNS_NSI*>(fOscCalc)->SetDelta_mutau(itr.Value());
336  } // end if NSI oscillations
337  } // end if PMNS or NSI oscillations
338  else if(fCalcType == cmf::kSterile){
339  if ((cmf::OscParm_t)itr.Key() == cmf::kRho ) dynamic_cast<osc::OscCalcSterileEigen*>(fOscCalc)->SetRho ( itr.Value());
340  else if((cmf::OscParm_t)itr.Key() == cmf::kL ) dynamic_cast<osc::OscCalcSterileEigen*>(fOscCalc)->SetL ( itr.Value());
341  else if((cmf::OscParm_t)itr.Key() == cmf::kTh12 ) dynamic_cast<osc::OscCalcSterileEigen*>(fOscCalc)->SetAngle(1, 2, itr.Value());
342  else if((cmf::OscParm_t)itr.Key() == cmf::kTh13 ) dynamic_cast<osc::OscCalcSterileEigen*>(fOscCalc)->SetAngle(1, 3, itr.Value());
343  else if((cmf::OscParm_t)itr.Key() == cmf::kTh23 ) dynamic_cast<osc::OscCalcSterileEigen*>(fOscCalc)->SetAngle(2, 3, itr.Value());
344  else if((cmf::OscParm_t)itr.Key() == cmf::kTh24 ) dynamic_cast<osc::OscCalcSterileEigen*>(fOscCalc)->SetAngle(2, 4, itr.Value());
345  else if((cmf::OscParm_t)itr.Key() == cmf::kTh34 ) dynamic_cast<osc::OscCalcSterileEigen*>(fOscCalc)->SetAngle(3, 4, itr.Value());
346  else if((cmf::OscParm_t)itr.Key() == cmf::kdCP ) dynamic_cast<osc::OscCalcSterileEigen*>(fOscCalc)->SetDelta(1, 3, itr.Value());
347  else if((cmf::OscParm_t)itr.Key() == cmf::kd24 ) dynamic_cast<osc::OscCalcSterileEigen*>(fOscCalc)->SetDelta(2, 4, itr.Value());
348  else if((cmf::OscParm_t)itr.Key() == cmf::kDmsq41) dynamic_cast<osc::OscCalcSterileEigen*>(fOscCalc)->SetDm (4, itr.Value());
349  else if((cmf::OscParm_t)itr.Key() == cmf::kDmsq21) dynamic_cast<osc::OscCalcSterileEigen*>(fOscCalc)->SetDm (2, itr.Value());
350  else if((cmf::OscParm_t)itr.Key() == cmf::kDmsq32){
351  dynamic_cast<osc::OscCalcSterileEigen*>(fOscCalc)->SetDm(3, itr.Value() + fCurrentLoc.AsOscillationParameterMap().find(cmf::kDmsq21)->second);
352  }
353  } // end if sterile oscillations
354  } // end loop over detector location
355  }
356 
357  //----------------------------------------------------------------------------
359  cmf::MetaData const& md)
360  {
361  fCurrentEvent = curEvent.get();
362  fCurrentMD = md;
363 
364  MF_LOG_DEBUG("SetCurrentEvent")
365  << " Metadata is "
366  << md.ToString();
367 
368  // for(auto const& itr : fCurrentEvent->MCVals().SystematicShifts()){
369  // MF_LOG_VERBATIM("EventListManipulator")
370  // << "ShifterAndWeighter::SetCurrentEvent() systematic: "
371  // << itr;
372  // }
373 
374  // check the validity of the mcVals before trying to get weights, if it is
375  // a nullptr, we have a problem.
376  if( !fCurrentEvent->MCValsFilled() && md.isMC)
377  throw cet::exception("ShifterAndWeighter")
378  << "Event does not have a valid MCVarVals unique_ptr, bail";
379 
380  // the data values better always be valid
381  if( !fCurrentEvent->DataValsFilled() )
382  throw cet::exception("ShifterAndWeighter")
383  << "Event does not have a valid DataVarVals unique_ptr, bail";
384 
385  // reset the oscillation baseline based on the detector
386  fOscCalc->SetL(this->CurrentParameterValue("L"));
387 
388  }
389 
390  //----------------------------------------------------------------------------
392  {
393  for(auto const& itr : fCurrentLoc.DetectorLocations()){
394  if(fCurrentMD.detector == itr.Detector())
395  return this->ShiftForParameter(itr, param);
396  }
397 
398  // if we got here, we don't have the desired detector in the locations, a problem
399  throw cet::exception("ShifterAndWeighter")
400  << "no detector location found for "
402 
403  // make the compiler happy
404  return this->ShiftForParameter(fCurrentLoc.DetectorLocations().front(), param);
405  }
406 
407  //----------------------------------------------------------------------------
409  uint8_t const& param)
410  {
411 
412  float totalShift = 1.;
413  float shiftScale;
414  float curShift;
415 
416  for(auto const& itr : fShiftersToUse){
417 
418  if(itr.key != param ||
419  itr.detector != loc.Detector() ) continue;
420 
421  auto const& sigma = loc.Parameter(itr.shiftName).Sigma();
422 
423  curShift = this->CurrentParameterValue(itr.shiftName);
424 
425  shiftScale = 1. + sigma * curShift;
426 
427  // MF_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  if(fWeightType == cmf::kNoWgt) return 1.;
454 
455  // set the weight to be the fake weight value, which is 1 for
456  // data.
457  double wgt = curEvent->DataVals().FakeWeight();
458 
459 
460  // if the mcVals is null or the metadata claims it isn't MC,
461  // return the fake weight for the event. If the event is data
462  // then the value is always 1, if it is fake data, then the
463  // value is the proper weighting
464  if(!curEvent->MCValsFilled() ||
465  !md.isMC ||
466  md.IsCosmicMuon() )
467  return wgt;
468 
469  //if(md.selectionType == cmf::kNuESelectionPeripheral)
470  // MF_LOG_VERBATIM("ShifterAndWeighter")
471  // << "weight from systematics is "
472  // << wgt;
473 
474  // check if we only want weights from the systematic parameters or only
475  // from the oscillation parameters
477  wgt *= this->TotalWeightFromWeighters();
478  }
479  else if(fWeightType == cmf::kOscWgtOnly){
480  wgt *= this->OscillationWeight();
481  }
482  else
483  wgt *= this->OscillationWeight() * this->TotalWeightFromWeighters();
484 
485  // if(md.selectionType == cmf::kNuESelectionPeripheral)
486  // MF_LOG_VERBATIM("ShifterAndWeighter")
487  // << "weight from oscillations is "
488  // << this->OscillationWeight();
489 
490  MF_LOG_DEBUG("ShifterAndWeighter")
491  << "Oscillation weight is "
492  << wgt
493  << " energy: "
495  << " pdg: "
497  << " orig pdg: "
499  << " detector: "
500  << fCurrentMD.detector;
501 
502  // if(md.selectionType == cmf::kNuESelectionPeripheral)
503  // MF_LOG_VERBATIM("ShifterAndWeighter")
504  // << "total weight is "
505  // << wgt;
506 
507  return wgt;
508  }
509 
510  //----------------------------------------------------------------------------
512  {
513  double wgt = 1.;
514 
515  std::string ss;
516  for(auto const& detWgt : fWeightersToUse){
517 
518  if(detWgt.detector != fCurrentMD.detector) continue;
519  ss += "\t" + detWgt.wgtName +
520  " " + detWgt.detector +
521  " " + std::to_string(detWgt.wgtFn()) +
522  " " + std::to_string(wgt) +
523  " " + std::to_string(wgt * detWgt.wgtFn()) +
524  "";
525  wgt *= detWgt.wgtFn();
526  }
527  // MF_LOG_VERBATIM("SAWWeighterInfo")
528  // << "detWgt:"
529  // << ss;
530 
531  return wgt;
532  }
533 
534  //----------------------------------------------------------------------------
536  {
537  // nothing to change
538  if(fOscCalc != nullptr && fCalcType == calcType) return;
539  else{
540  MF_LOG_WARNING("ShifterAndWeighter")
541  << "resetting oscillation calculator to type "
542  << calcType
543  << " from "
544  << fCalcType
545  << " are you sure you meant to change types in this job?";
546 
547  delete fOscCalc;
548  fOscCalc = nullptr;
549  }
550 
551  fCalcType = calcType;
552 
554  else if(fCalcType == cmf::kNSI ) fOscCalc = new osc::OscCalcPMNS_NSI();
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 &&
565  fCurrentMD.detector == cmf::kMINOSNEAR)) return 1.;
566 
568 
569  return this->OscillationWeightBinCenter();
570  }
571  //----------------------------------------------------------------------------
573  {
574  double E = fCurrentEvent->MCVals().val_at(cmf::kTrueE);
576  int pdgOrig = fCurrentEvent->MCVals().val_at(cmf::kTruePDGOrig);
577 
578  // only attempt to weight neutrinos for oscillations
579  if(std::abs(pdg) != 12 &&
580  std::abs(pdg) != 14 &&
581  std::abs(pdg) != 16) return 1.;
582 
583  return fOscCalc->P(pdgOrig, pdg, E);
584 
585  }
586 
587  //----------------------------------------------------------------------------
589  {
590  if(fTrueEnergyBins.empty()){
591  const int numTrueEnergyBins = 100;
592 
593  // How many edges to generate. Allow room for 0-Emin bin
594  const double N = numTrueEnergyBins-1;
595  const double A = N * 0.5; // 500 MeV: there's really no events below there
596 
597  fTrueEnergyBins.insert(0);
598  for(int i = 1; i < N - 1; ++i) fTrueEnergyBins.insert(A/i);
599  fTrueEnergyBins.insert(120.);
600  }
601 
602  double E = fCurrentEvent->MCVals().val_at(cmf::kTrueE);
604  int pdgOrig = fCurrentEvent->MCVals().val_at(cmf::kTruePDGOrig);
605 
606  // get the relevant bin for the current neutrino, then get the
607  // energy at the center of that bin
608  double eUpper = *(fTrueEnergyBins.rbegin());
609  double eLower = *(fTrueEnergyBins.begin());
610 
611  // find the bin edge that goes after the value of E
612  auto itr = fTrueEnergyBins.upper_bound(E);
613  if(itr != fTrueEnergyBins.end()) eUpper = *itr;
614  else --itr;
615 
616  // now find the edge before by just decrementing the iterator
617  --itr;
618  if(itr != fTrueEnergyBins.end()) eLower = *itr;
619 
620  // get the midpoint between the upper and lower values,
621  // ie take their mean
622  E = 0.5 * (eUpper + eLower);
623 
624  MF_LOG_DEBUG("ShifterAndWeighter")
625  << "true energy: "
627  << " central bin energy "
628  << E;
629 
630  // only attempt to weight neutrinos for oscillations
631  if(std::abs(pdg) != 12 &&
632  std::abs(pdg) != 14 &&
633  std::abs(pdg) != 16) return 1.;
634 
635  return fOscCalc->P(pdgOrig, pdg, E);
636 
637  }
638 
639  //----------------------------------------------------------------------------
641  std::map<float, float> const& wgtsBySigma)
642  {
643  auto mapEnd = wgtsBySigma.end();
644 
645  // no way this is a valid number and is obviously stupid
646  // not using a std::numeric_limits value to avoid rounding to 0
647  double badVal = -99999999999999.;
648  double minus2Val = (wgtsBySigma.find(-2.) != mapEnd) ? static_cast<double>(wgtsBySigma.find(-2.)->second) : badVal;
649  double minus1Val = (wgtsBySigma.find(-1.) != mapEnd) ? static_cast<double>(wgtsBySigma.find(-1.)->second) : badVal;
650  // double zeroVal = 1.;
651  double zeroVal = (wgtsBySigma.find( 0.) != mapEnd) ? static_cast<double>(wgtsBySigma.find( 0.)->second) : 1.;
652  double plus1Val = (wgtsBySigma.find( 1.) != mapEnd) ? static_cast<double>(wgtsBySigma.find( 1.)->second) : badVal;
653  double plus2Val = (wgtsBySigma.find( 2.) != mapEnd) ? static_cast<double>(wgtsBySigma.find( 2.)->second) : badVal;
654  double slope = 0.;
655  double intercept = 0.;
656 
657  //Debugging RPA weights
658  MF_LOG_DEBUG("ShifterAndWeighter")
659  << "\tCalcLinearInterpWgt: "
660  << " -2:" << minus2Val
661  << ", -1:" << minus1Val
662  << ", 0:" << zeroVal
663  << ", 1:" << plus1Val
664  << ", 2:" << plus2Val;
665 
666  /*
667  if(wgtsBySigma.find(-1.) != mapEnd)
668  MF_LOG_VERBATIM("ShifterAndWeighter")
669  << "\t\t 0 val from "
670  << wgtsBySigma.find(0.)->second;
671  */
672 
673  // The following code assumes that at least the -1 sigma, +1 sigma
674  // values are in the map. If they are not, we have a problem
675  if(minus1Val == badVal ||
676  plus1Val == badVal){
677  MF_LOG_WARNING("ShifterAndWeighterCalcLinearInterpWgt")
678  << "No values associated with either -1 sigma or +1 sigma "
679  << "so we cannot calculate a weight, return 1.";
680  return 1.;
681  }
682 
683  if(sigma <= -1){
684  // It is possible there is no -2 sigma value from the map,
685  // in that case just interpolate from 0 to -1 sigma
686  if(minus2Val != badVal){
687  slope = (minus1Val - minus2Val);
688  intercept = minus2Val - slope * (-2.);
689  }
690  else{
691  slope = (zeroVal - minus1Val);
692  intercept = minus1Val - slope * (-1.);
693  }
694  }
695  else if(sigma > -1. && sigma <= 0.){
696  slope = (zeroVal - minus1Val);
697  intercept = minus1Val - slope * (-1.);
698  }
699  else if(sigma > 0. && sigma <= 1.){
700  slope = (plus1Val - zeroVal);
701  intercept = plus1Val - slope * 1.;
702  }
703  else if(sigma > 1.){
704  // It is possible there is no +2 sigma value from the map,
705  // in that case just interpolate from 0 to +1 sigma
706  if(plus2Val != badVal){
707  slope = (plus2Val - plus1Val);
708  intercept = plus2Val - slope * 2.;
709  }
710  else{
711  slope = (plus1Val - zeroVal);
712  intercept = plus1Val - slope * 1.;
713  }
714  }
715 
716  double wgt = sigma * slope + intercept;
717 
718  // Don't let weights be negative
719  if(wgt < 0) wgt = 0.;
720 
721  return wgt;
722  }
723 
724  //-------------------------------------------------------------------------
726  // we have a separate handling for normalisation systematics because
727  // there's a bunch of numbers it's just not worth storing in the event
728  // list
729  double oneSigmaWgt = 1.0;
730 
731  double normPOT = 0.0055;
732  double normNDMass = 0.003;
733  double normFDMass = 0.003;
734  double normML = 0.008;
735  double normNDNumuContext = fCurrentMD.BeamType() == cmf::kFHC ? 0.0021 : 0.0048;
736  double normNDNueContext = fCurrentMD.BeamType() == cmf::kFHC ? 0.0041 : 0.0053;
737  double normNDNCContext = fCurrentMD.BeamType() == cmf::kFHC ? 0.01 : 0.00;
738 
739  // TODO need to account for MINOS detectors perhaps
740  bool isND = fCurrentMD.detector == cmf::kNEARDET ? true : false;
741  bool isFD = fCurrentMD.detector == cmf::kFARDET ? true : false;
742  bool isFHC = fCurrentMD.BeamType() == cmf::kFHC ? true : false;
743  bool isRHC = fCurrentMD.BeamType() == cmf::kRHC ? true : false;
744  bool isNuMu = false;
745  bool isNuE = false;
746  bool isNC = false;
747 
753  isNuMu = true;
757  isNuE = true;
759  isNC = true;
760 
761  if (systName == "NormalizationPOT") oneSigmaWgt += normPOT;
762  else if (systName == "NormalizationNDMass" && isND) oneSigmaWgt += normNDMass;
763  else if (systName == "NormalizationFDMass" && isFD) oneSigmaWgt += normFDMass;
764  else if (systName == "NormalizationML" && isFD) oneSigmaWgt += normML;
765  else if (systName == "NormalizationFHCNuMuContext" && isND && isNuMu && isFHC) oneSigmaWgt += normNDNumuContext;
766  else if (systName == "NormalizationRHCNuMuContext" && isND && isNuMu && isRHC) oneSigmaWgt += normNDNumuContext;
767  else if (systName == "NormalizationFHCNuEContext" && isND && isNuE && isFHC) oneSigmaWgt += normNDNueContext;
768  else if (systName == "NormalizationRHCNuEContext" && isND && isNuE && isRHC) oneSigmaWgt += normNDNueContext;
769  else if (systName == "NormalizationFHCNCContext" && isND && isNC && isFHC) oneSigmaWgt += normNDNCContext;
770  else if (systName == "NormalizationRHCNCContext" && isND && isNC && isRHC) oneSigmaWgt += normNDNCContext;
771 
772  float sigShift = this->CurrentParameterValue(systName);
773  return 1. - (1. - oneSigmaWgt) * sigShift;
774 
775  }
776 
777  //-------------------------------------------------------------------------
779  {
780  // TODO may need to account for MINOS detectors
781  bool isWrongDet = false;
782  if (systName.find("ND") != std::string::npos && fCurrentMD.detector == cmf::kFARDET) isWrongDet = true;
783  if (systName.find("FD") != std::string::npos && fCurrentMD.detector == cmf::kNEARDET) isWrongDet = true;
784 
785  std::string strippedName = systName;
786  if ((systName.substr(systName.length()-2, systName.length()) == "ND") |
787  (systName.substr(systName.length()-2, systName.length()) == "FD"))
788  strippedName = systName.substr(0, systName.length()-2);
789 
790  std::string ratioName = strippedName
791  + "Total"
795  + "Plus1Sigma";
796 
798  int logicalBin = fCurrentEvent->DataVals().LogicalBin();
800 
801  float oneSigmaWgt = isWrongDet ? 1 : fileSystVec[ratioName][bin];
802  float sigShift = this->CurrentParameterValue(systName);
803 
804  return 1. - (1. - oneSigmaWgt) * sigShift;
805  }
806 
807  //----------------------------------------------------------------------------
809  {
810  // return 1 if this is a cosmic muon
812 
813  // if this weight has NC or CC in the name and this interaction type
814  // is not that, return 1
815  if((systName.find("CC") != std::string::npos && fCurrentMD.interactionType == cmf::kNC) ||
816  (systName.find("NC") == std::string::npos && fCurrentMD.interactionType == cmf::kNC) ) return 1.;
817 
818  // may need to come up with default for cosmic muons, however I think the fact that
819  // the default weights for every GENIE variable is a collection of 1's should
820  // mean this is OK
821 
822  std::map<float, float> wgtBySigma = fCurrentEvent->MCVals().SystVarWgts(cmf::VarNameToKey(systName));
823 
824  // for(auto itr : wgtBySigma)
825  //MF_LOG_VERBATIM("GENIESystWeight")
826  //<< systName
827  //<< " weights "
828  //<< itr.first
829  //<< " "
830  //<< itr.second;
831 
832  return this->CalcLinearInterpWgt(this->CurrentParameterValue(systName),
833  wgtBySigma);
834  }
835 
836  //----------------------------------------------------------------------------
838  {
839  double thisWeight = 1.0;
840  if (systName.find("XsecCV") != std::string::npos) thisWeight *= fCurrentEvent->MCVals().val_at(cmf::kXSecCV2020Wgt);
841  if (systName.find("PPFXFluxCV") != std::string::npos) thisWeight *= fCurrentEvent->MCVals().val_at(cmf::kPPFXFluxCVWgt);
842 
843  return thisWeight;
844  }
845 
846 } // cmf
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
enum cmf::osc_params OscParm_t
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)
std::vector< cmf::Parameter > const & Parameters() const
Definition: Parameter.h:43
DetBeamSelSet const & SelectionsToUse()
bool MCValsFilled() const
Definition: Event.h:73
int LogicalBinToSelectionBin(int logicalBin, cmf::MetaData md)
double XSecCVPPFXWeight(std::string const &systName)
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
const std::vector< std::string > cDetType_Strings({"UnknownDet","NearDet","FarDet","MINOSNear","MINOSFar","AllDetectors"})
Float_t ss
Definition: plot.C:24
static SelectionUtility * Instance()
virtual void SetTh13(const T &th13)=0
cmf::DetType_t const Detector() const
Definition: Parameter.h:44
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)
cmf::OscillationParameterMap AsOscillationParameterMap(cmf::DetType_t const &det=cmf::kFARDET) const
Definition: Parameter.cxx:321
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
cmf::Parameter const & Parameter(std::string const &parName) const
Definition: Parameter.cxx:30
enum cmf::par_type ParType_t
Eigen-based 3+1 sterile oscillation calculator.
std::vector< ShifterInfo > fShiftersToUse
double CMFSystVarWeight(std::string const &systName)
float const & Sigma() const
Definition: Parameter.h:138
static uint8_t VarNameToKey(std::string const &name)
Definition: StaticFuncs.h:395
cmf::SelectionType_t selectionType
Definition: Structs.h:116
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
std::vector< WeighterInfo > fWeightersToUse
double const CurrentParameterValue(std::string const &parName)
const std::vector< std::string > cSelectionType_Strings({"NuESel_AllPID","NuESel_LowPID","NuESel_MidPID","NuESel_HighPID","NuESel_Peripheral","NuMuSel","NuMuSelQ1","NuMuSelQ2","NuMuSelQ3","NuMuSelQ4","NCSel","UnknownSel"})
const std::vector< std::string > cBeamType_Strings({"FHC","RHC","0HC","UnknownBeam"})
float LogicalBin() const
Definition: VarVals.h:294
static ShifterAndWeighter * Instance()
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)
std::string ToString() const
Definition: Structs.cxx:135
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
std::vector< ParameterSpaceLoc > const & DetectorLocations() const
Definition: Parameter.h:193
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
#define MF_LOG_VERBATIM(category)
float val_at(uint8_t const &varkey) const
Access a Var by name.
Definition: VarVals.cxx:18
virtual void SetTh23(const T &th23)=0
#define MF_LOG_DEBUG(id)
float ShiftForParameter(cmf::ParameterSpaceLoc const &loc, uint8_t const &param)
cmf::ParameterSpaceLoc const & DetectorLocation(cmf::DetType_t const &det) const
Definition: Parameter.cxx:339
bool fTrueEOscillationWeight
flag to calculate the oscillation weight at the true E_nu vs bin center
TFile * file
Definition: cellShifts.C:17
std::vector< Parameter > OscillationParameters(cmf::DetType_t const &det=cmf::kFARDET) const
Definition: Parameter.cxx:357
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:184
#define MF_LOG_WARNING(category)
cmf::InteractionType_t interactionType
Definition: Structs.h:117
float RecoEnergy() const
Definition: VarVals.h:293
double FileSystWeight(std::string const &systName)
void SetCurrentVals(cmf::Location const &loc)
void InitShiftsAndWeightsForDetector(cmf::DetType_t const &det, cmf::ParameterSpaceLoc const &detLoc)
std::string to_string(ModuleType const mt)
Definition: ModuleType.h:34
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:302
virtual void SetTh12(const T &th12)=0
virtual void SetdCP(const T &dCP)=0
static CovarianceBinUtility * Instance()
enum BeamMode string