diff --git a/src/Reweight/NOvARwgtEngine.cxx b/src/Reweight/NOvARwgtEngine.cxx index 1277db9..6002f79 100644 --- a/src/Reweight/NOvARwgtEngine.cxx +++ b/src/Reweight/NOvARwgtEngine.cxx @@ -1,232 +1,316 @@ #include "NOvARwgtEngine.h" #include "NOvARwgt/interfaces/GenieInterface.h" #include "NOvARwgt/rwgt/tunes/Tunes2017.h" #include "NOvARwgt/rwgt/tunes/Tunes2018.h" -#include "NOvARwgt/rwgt/tunes/Tunes2019.h" +#include "NOvARwgt/rwgt/tunes/Tunes2020.h" #include "NOvARwgt/rwgt/tunes/TunesSA.h" #include "NOvARwgt/rwgt/EventRecord.h" -#include "NOvARwgt/rwgt/Tune.h" #include "NOvARwgt/rwgt/ISystKnob.h" +#include "NOvARwgt/rwgt/Tune.h" + +// GENIEv3 +#include "Framework/Utils/AppInit.h" +#include "Framework/Utils/RunOpt.h" static size_t const kCVTune2017 = 100; static size_t const kCVTune2018 = 200; -static size_t const kCVTune2019 = 300; +static size_t const kCVTune2020 = 300; static size_t const kCVTuneSA = 400; static size_t const kNoSuchKnob = std::numeric_limits::max(); +novarwgt::Tune const *TuneFactory(size_t e) { + switch (e) { + case kCVTune2017: { + return &novarwgt::kCVTune2017; + } + case kCVTune2018: { + return &novarwgt::kCVTune2018; + } + case kCVTune2020: { + return &novarwgt::kCVTune2020; + } + case kCVTuneSA: { + return &novarwgt::kCVTuneSA; + } + default: { + return NULL; + } + } +} + +void NOvARwgtEngine::InitializeKnobs() { + + fTunes.push_back(TuneFactory(kCVTune2017)); + fTuneEnums[kCVTune2017] = 0; + fTunes.push_back(TuneFactory(kCVTune2018)); + fTuneEnums[kCVTune2018] = 1; + fTunes.push_back(TuneFactory(kCVTune2020)); + fTuneEnums[kCVTune2020] = 2; + fTunes.push_back(TuneFactory(kCVTuneSA)); + fTuneEnums[kCVTuneSA] = 3; + fTuneInUse = {false, false, false, false}; + fTuneValues = {0, 0, 0, 0}; + + size_t ctr = 0; + + NUIS_LOG(FIT, "NOvARwgt kCVTune2017 sub-knobs:"); + size_t tune_ctr = 0; + for (auto &k : novarwgt::kCVTune2017.SystKnobs()) { + NUIS_LOG(FIT, "\t" << k.first); + fKnobs.push_back(k.second); + fKnobTunes.push_back(&novarwgt::kCVTune2017); + fKnobTuneidx.push_back(0); + fKnobEnums[kCVTune2017 + 1 + tune_ctr++] = ctr++; + fKnobInUse.push_back(false); + fKnobValues.push_back(0); + } + + tune_ctr = 0; + NUIS_LOG(FIT, "NOvARwgt kCVTune2018 sub-knobs:"); + for (auto &k : novarwgt::kCVTune2018.SystKnobs()) { + NUIS_LOG(FIT, "\t" << k.first); + fKnobs.push_back(k.second); + fKnobTunes.push_back(&novarwgt::kCVTune2018); + fKnobTuneidx.push_back(1); + fKnobEnums[kCVTune2018 + 1 + tune_ctr++] = ctr++; + fKnobInUse.push_back(false); + fKnobValues.push_back(0); + } + + tune_ctr = 0; + NUIS_LOG(FIT, "NOvARwgt kCVTune2020 sub-knobs:"); + for (auto &k : novarwgt::kCVTune2020.SystKnobs()) { + NUIS_LOG(FIT, "\t" << k.first); + fKnobs.push_back(k.second); + fKnobTunes.push_back(&novarwgt::kCVTune2020); + fKnobTuneidx.push_back(2); + fKnobEnums[kCVTune2020 + 1 + tune_ctr++] = ctr++; + fKnobInUse.push_back(false); + fKnobValues.push_back(0); + } + + tune_ctr = 0; + NUIS_LOG(FIT, "NOvARwgt kCVTuneSA sub-knobs:"); + for (auto &k : novarwgt::kCVTuneSA.SystKnobs()) { + NUIS_LOG(FIT, "\t" << k.first); + fKnobs.push_back(k.second); + fKnobTunes.push_back(&novarwgt::kCVTuneSA); + fKnobTuneidx.push_back(3); + fKnobEnums[kCVTuneSA + 1 + tune_ctr++] = ctr++; + fKnobInUse.push_back(false); + fKnobValues.push_back(0); + } +} + +void NOvARwgtEngine::InitializeGENIE() { + genie::RunOpt *grunopt = genie::RunOpt::Instance(); + grunopt->EnableBareXSecPreCalc(true); + grunopt->SetEventGeneratorList(Config::GetParS("GENIEEventGeneratorList")); + if (!Config::HasPar("GENIETune")) { + NUIS_ABORT( + "GENIE tune was not specified, this is required when reweighting GENIE " + "V3+ events. Add a config parameter like: to your nuisance card."); + } + grunopt->SetTuneName(Config::GetParS("GENIETune")); + grunopt->BuildTune(); + std::string genv = + std::string(getenv("GENIE")) + "/config/Messenger_laconic.xml"; + genie::utils::app_init::MesgThresholds(genv); +} + size_t NOvARwgtEngine::GetWeightGeneratorIndex(std::string const &strname) { - int upos = strname.find_first_of("_"); + size_t upos = strname.find_first_of("_"); if (strname.find("CVTune2017") == 0) { if (upos == std::string::npos) { return kCVTune2017; } std::string knobname = strname.substr(upos + 1); - if (novarwgt::kCVTune2017.SystKnobs().contains(knobname)) { - auto loc = std::find(fKnobs.begin(), fKnobs.end(), - novarwgt::kCVTune2017.SystKnobs()[knobname]); - if (loc == fKnobs.end()) { - size_t rtn = kCVTune2017 + 1 + fKnobs.size(); - fKnobs.push_back(novarwgt::kCVTune2017.SystKnobs()[knobname]); - return rtn; - } else { - return kCVTune2017 + 1 + std::distance(fKnobs.begin(), loc); - } + if (novarwgt::kCVTune2017.SystKnobs().count(knobname)) { + auto loc = novarwgt::kCVTune2017.SystKnobs().find(knobname); + + return kCVTune2017 + 1 + + std::distance(novarwgt::kCVTune2017.SystKnobs().begin(), loc); } + } else if (strname.find("CVTune2018") == 0) { if (upos == std::string::npos) { return kCVTune2018; } std::string knobname = strname.substr(upos + 1); - if (novarwgt::kCVTune2018.SystKnobs().contains(knobname)) { - auto loc = std::find(fKnobs.begin(), fKnobs.end(), - novarwgt::kCVTune2018.SystKnobs()[knobname]); - if (loc == fKnobs.end()) { - size_t rtn = kCVTune2018 + 1 + fKnobs.size(); - fKnobs.push_back(novarwgt::kCVTune2018.SystKnobs()[knobname]); - return rtn; - } else { - return kCVTune2018 + 1 + std::distance(fKnobs.begin(), loc); - } - } - } else if (strname.find("CVTune2019") == 0) { + if (novarwgt::kCVTune2018.SystKnobs().count(knobname)) { + auto loc = novarwgt::kCVTune2018.SystKnobs().find(knobname); + + return kCVTune2018 + 1 + + std::distance(novarwgt::kCVTune2018.SystKnobs().begin(), loc); + } + } else if (strname.find("CVTune2020") == 0) { if (upos == std::string::npos) { - return kCVTune2019; + return kCVTune2020; } std::string knobname = strname.substr(upos + 1); - if (novarwgt::kCVTune2019.SystKnobs().contains(knobname)) { - auto loc = std::find(fKnobs.begin(), fKnobs.end(), - novarwgt::kCVTune2019.SystKnobs()[knobname]); - if (loc == fKnobs.end()) { - size_t rtn = kCVTune2019 + 1 + fKnobs.size(); - fKnobs.push_back(novarwgt::kCVTune2019.SystKnobs()[knobname]); - return rtn; - } else { - return kCVTune2019 + 1 + std::distance(fKnobs.begin(), loc); - } + if (novarwgt::kCVTune2020.SystKnobs().count(knobname)) { + auto loc = novarwgt::kCVTune2020.SystKnobs().find(knobname); + + return kCVTune2020 + 1 + + std::distance(novarwgt::kCVTune2020.SystKnobs().begin(), loc); } } else if (strname.find("CVTuneSA") == 0) { if (upos == std::string::npos) { return kCVTuneSA; } std::string knobname = strname.substr(upos + 1); - if (novarwgt::kCVTuneSA.SystKnobs().contains(knobname)) { - auto loc = std::find(fKnobs.begin(), fKnobs.end(), - novarwgt::kCVTuneSA.SystKnobs()[knobname]); - if (loc == fKnobs.end()) { - size_t rtn = kCVTuneSA + 1 + fKnobs.size(); - fKnobs.push_back(novarwgt::kCVTuneSA.SystKnobs()[knobname]); - return rtn; - } else { - return kCVTuneSA + 1 + std::distance(fKnobs.begin(), loc); - } + if (novarwgt::kCVTuneSA.SystKnobs().count(knobname)) { + auto loc = novarwgt::kCVTuneSA.SystKnobs().find(knobname); + + return kCVTuneSA + 1 + + std::distance(novarwgt::kCVTuneSA.SystKnobs().begin(), loc); } } return kNoSuchKnob; } -novarwgt::Tune const *TuneFactory(size_t e) { - switch (e) { - case kTune2017: { - return &novarwgt::kCVTune2017; - } - case kTune2018: { - return &novarwgt::kCVTune2018; - } - case kTune2019: { - return &novarwgt::kCVTune2019; - } - case kTuneSA: { - return &novarwgt::kCVTuneSA; - } - default: { return NULL; } - } -} - void NOvARwgtEngine::IncludeDial(std::string name, double startval) { - size_t we_e = GetWeightGeneratorIndex(name); - if (we_e == kNoSuchKnob) { + size_t we_indx = GetWeightGeneratorIndex(name); + if (we_indx == kNoSuchKnob) { NUIS_ABORT("[ERROR]: Invalid NOvARwgt Engine name: " << name); } - bool IsTune = !(we_e % 100); + bool IsTune = !(we_indx % 100); if (IsTune) { - if (fTuneEnums.find(we_e) != fTuneEnums.end()) { - NUIS_ABORT("[ERROR]: NOvARwgt Tune name: " << name - << " already included."); - } - fTuneEnums[we_e] = fTunes.size(); - fTunes.push_back(TuneFactory(we_e)); - fTuneValues.push_back(startval); + auto tune_idx = fTuneEnums[we_indx]; + fTuneValues[tune_idx] = startval; + fTuneInUse[tune_idx] = true; } else { - // The GetWeightGeneratorIndex caches the Knob so we don't have to. - fKnobEnums[we_e] = fKnobValues.size(); - fKnobValues.push_back(startval); + auto knob_idx = fKnobEnums[we_indx]; + fKnobValues[knob_idx] = startval; + fKnobInUse[knob_idx] = true; } }; void NOvARwgtEngine::SetDialValue(int nuisenum, double val) { size_t we_indx = (nuisenum % 1000); - bool IsTune = !(we_e % 100); + bool IsTune = !(we_indx % 100); if (IsTune) { - if (!fTuneEnums.count(we_indx)) { + auto tune_idx = fTuneEnums[we_indx]; + if (!fTuneInUse[tune_idx]) { NUIS_ABORT("[ERROR]: SetDialValue for NOvARwgt dial: " << we_indx << " but that tune hasn't been included yet."); } - size_t engine_index = fTuneEnums[we_indx]; - fTuneValues[engine_index] = val; + fTuneValues[tune_idx] = val; } else { - if (!fKnobEnums.count(we_indx)) { + auto knob_idx = fKnobEnums[we_indx]; + if (!fKnobInUse[knob_idx]) { NUIS_ABORT("[ERROR]: SetDialValue for NOvARwgt dial: " << we_indx << " but that Knob hasn't been included yet."); } - size_t engine_index = fKnobEnums[we_indx]; - fKnobValues[engine_index] = val; + fKnobValues[knob_idx] = val; } } void NOvARwgtEngine::SetDialValue(std::string name, double val) { SetDialValue(GetWeightGeneratorIndex(name), val); } bool NOvARwgtEngine::IsDialIncluded(std::string name) { return IsDialIncluded(GetWeightGeneratorIndex(name)); } bool NOvARwgtEngine::IsDialIncluded(int nuisenum) { size_t we_indx = (nuisenum % 1000); - bool IsTune = !(we_e % 100); + bool IsTune = !(we_indx % 100); if (IsTune) { - return fTuneEnums.count(we_indx); + auto tune_idx = fTuneEnums[we_indx]; + return fTuneInUse[tune_idx]; } else { - return fKnobEnums.count(we_indx); + auto knob_idx = fKnobEnums[we_indx]; + return fKnobInUse[knob_idx]; } } double NOvARwgtEngine::GetDialValue(std::string name) { - return GetDialValueGetWeightGeneratorIndex(name)); + return GetDialValue(GetWeightGeneratorIndex(name)); } double NOvARwgtEngine::GetDialValue(int nuisenum) { size_t we_indx = (nuisenum % 1000); if (we_indx == kNoSuchKnob) { - NUIS_ABORT("[ERROR]: Invalid NOvARwgt Engine name: " << name); + NUIS_ABORT("[ERROR]: Invalid NOvARwgt Engine enum: " << nuisenum); } - bool IsTune = !(we_e % 100); + bool IsTune = !(we_indx % 100); if (IsTune) { - if (!fTuneEnums.count(we_indx)) { - NUIS_ABORT("[ERROR]: SetDialValue for NOvARwgt dial: " + auto tune_idx = fTuneEnums[we_indx]; + if (!fTuneInUse[tune_idx]) { + NUIS_ABORT("[ERROR]: GetDialValue for NOvARwgt dial: " << we_indx << " but that tune hasn't been included yet."); } - return fTuneValues[fTuneEnums[we_indx]]; + return fTuneValues[tune_idx]; } else { - - if (!fKnobEnums.count(we_indx)) { - NUIS_ABORT("[ERROR]: SetDialValue for NOvARwgt dial: " - << we_indx << " but that engine hasn't been included yet."); + auto knob_idx = fKnobEnums[we_indx]; + if (!fKnobInUse[knob_idx]) { + NUIS_ABORT("[ERROR]: GetDialValue for NOvARwgt dial: " + << we_indx << " but that Knob hasn't been included yet."); } - return fKnobValues[fKnobEnums[we_indx]]; + return fKnobValues[knob_idx]; } } double NOvARwgtEngine::CalcWeight(BaseFitEvt *evt) { double rw_weight = 1.0; // Make nom weight if (!evt) { NUIS_ABORT("evt not found : " << evt); } // Skip Non GENIE if (evt->fType != kGENIE) return 1.0; if (!(evt->genie_event)) { NUIS_ABORT("evt->genie_event not found!" << evt->genie_event); } if (!(evt->genie_event->event)) { NUIS_ABORT("evt->genie_event->event GHepRecord not found!" << (evt->genie_event->event)); } novarwgt::EventRecord rcd = novarwgt::ConvertGenieEvent(evt->genie_event->event); - for (size_t w_it = 0; w_it < fKnobs.size(); ++w_it) { - rw_weight *= fKnobs[w_it]->CalcWeight(fKnobValues[w_it], rcd); + for (size_t k_it = 0; k_it < fKnobs.size(); ++k_it) { + if (!fKnobInUse[k_it]) { + continue; + } + + double wght = fKnobTunes[k_it]->EventSystKnobWeight( + fKnobs[k_it]->Name(), fKnobValues[k_it], rcd, {}, false); + + // have to divide out the CV weight for this, ugly hack because the last + // parameter doesn't do what I want + if (fTuneInUse[fKnobTuneidx[k_it]] && fTuneValues[fKnobTuneidx[k_it]]) { + wght /= fKnobTunes[k_it]->EventSystKnobWeight(fKnobs[k_it]->Name(), 0, + rcd, {}, false); + } + + rw_weight *= wght; } - for (size_t w_it = 0; w_it < fTunes.size(); ++w_it) { - if (fTuneValues[w_it] == - 0) { // if a dial is set to 0, don't include its weight + for (size_t k_it = 0; k_it < fTunes.size(); ++k_it) { + if (!fTuneInUse[k_it]) { + continue; + } + if (!fTuneValues[k_it]) { continue; } - rw_weight *= fTunes[w_it]->EventWeight(rcd); + double wght = fTunes[k_it]->EventWeight(rcd); + rw_weight *= wght; } return rw_weight; } diff --git a/src/Reweight/NOvARwgtEngine.h b/src/Reweight/NOvARwgtEngine.h index 33d96d1..9f3dbeb 100644 --- a/src/Reweight/NOvARwgtEngine.h +++ b/src/Reweight/NOvARwgtEngine.h @@ -1,41 +1,50 @@ #pragma once #include "WeightEngineBase.h" namespace novarwgt { class ISystKnob; class Tune; -} +} // namespace novarwgt class NOvARwgtEngine : public WeightEngineBase { public: - NOvARwgtEngine(){}; + NOvARwgtEngine() { + InitializeKnobs(); + InitializeGENIE(); + }; virtual ~NOvARwgtEngine(){}; + void InitializeKnobs(); + void InitializeGENIE(); static size_t GetWeightGeneratorIndex(std::string const &); // Functions requiring Override void IncludeDial(std::string name, double startval); void SetDialValue(int nuisenum, double val); void SetDialValue(std::string name, double val); bool IsDialIncluded(std::string name); bool IsDialIncluded(int nuisenum); double GetDialValue(std::string name); double GetDialValue(int nuisenum); void Reconfigure(bool silent){}; double CalcWeight(BaseFitEvt *evt); bool NeedsEventReWeight() { return true; } std::map fTuneEnums; std::vector fTunes; std::vector fTuneValues; + std::vector fTuneInUse; std::map fKnobEnums; std::vector fKnobs; + std::vector fKnobTunes; + std::vector fKnobTuneidx; std::vector fKnobValues; + std::vector fKnobInUse; };