diff --git a/src/InputHandler/FitParticle.h b/src/InputHandler/FitParticle.h index 5776d5e..34c6362 100644 --- a/src/InputHandler/FitParticle.h +++ b/src/InputHandler/FitParticle.h @@ -1,109 +1,111 @@ // Copyright 2016-2021 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #ifndef FITPARTICLE_H_SEEN #define FITPARTICLE_H_SEEN /*! * \addtogroup InputHandler * @{ */ #include "TLorentzVector.h" +#include + /// Partle state flags for its position in the event enum particle_state{ kUndefinedState = 5, kInitialState = 0, kFSIState = 1, kFinalState = 2, kNuclearInitial = 3, kNuclearRemnant = 4, }; /// Condensed FitParticle class which acts a common format between the generators class FitParticle { public: /// Create particle of given pdg from momentum variables and state FitParticle(double x, double y, double z, double t, int pdg, Int_t state); /// Create empty particle (zero momentum) FitParticle(){}; ~FitParticle(){}; /// Used to change values after creation void SetValues(double x, double y, double z, double t, int pdg, Int_t state); /// Return Status Code according to particle_state enum inline int Status (void) const { return fStatus; }; /// Get Particle PDF inline int PDG (void) const { return fPID; }; /// Check if particle makes it to final state inline bool IsFinalState (void) const { return (fStatus == kFinalState); }; /// Was particle involved in intermediate state inline bool IsFSIState (void) const { return (fStatus == kFSIState); }; /// Was particle incoming inline bool IsInitialState (void) const { return (fStatus == kInitialState); }; /// Get Mass inline double M (void){ return fP.Mag(); }; /// Get Kinetic Energy inline double KE (void){ return fP.E() - fP.Mag(); }; /// Get Total Energy inline double E (void){ return fP.E(); }; /// Get 4 Momentum inline TLorentzVector P4(void) {return fP;}; /// Get 3 Momentum inline TVector3 P3(void) {return fP.Vect();}; /// Get 3 momentum magnitude inline double p(void) { return fP.Vect().Mag(); }; /// Get 3 momentum magnitude squared inline double p2(void) { return fP.Vect().Mag2(); }; /// Data Members TLorentzVector fP; ///< Particle 4 Momentum int fPID; ///< Particle PDG Code int fIsAlive; ///< Whether the particle is alive at the end of the event (Yes 1, No 0, Other? -1) int fNEUTStatusCode; ///< Particle Status (Incoming 1, FSI 2, Outgoing 0, Other 3) double fMass; ///< Particle Mass int fStatus; ///< State corresponding to particle_state enum bool fIsPrimary; ///< Primary target }; inline std::ostream& operator<<(std::ostream& os, FitParticle const& p){ return os << " Particle[pdgc:" << p.fPID << ", stat:"<. *******************************************************************************/ #include "nusystematicsWeightEngine.h" #include #include nusystematicsWeightEngine::nusystematicsWeightEngine() { fUseCV = false; Config(); } void nusystematicsWeightEngine::Config() { std::vector DuneRwtParam = Config::QueryKeys("DUNERwt"); if (DuneRwtParam.size() < 1) { NUIS_ABORT( "Instantiaged nusystematicsWeightEngine but without specifying a " "DUNERwt element that leads the way to the configuration."); } std::string fhicl_name = DuneRwtParam.front().GetS("ConfigFHiCL"); DUNErwt.LoadConfiguration(fhicl_name); } systtools::paramId_t const kNuSystCVResponse = 999; int nusystematicsWeightEngine::ConvDial(std::string name) { if (name == "NuSystCVResponse") { return kNuSystCVResponse; } if (!DUNErwt.HaveHeader(name)) { NUIS_ABORT("nusystematicsWeightEngine passed dial: " << name << " that it does not understand."); } + NUIS_LOG(FIT, "Added NuSyst param, " << name << " with ID: " << DUNErwt.GetHeaderId(name)); return DUNErwt.GetHeaderId(name); } void nusystematicsWeightEngine::IncludeDial(std::string name, double startval) { systtools::paramId_t DuneRwtEnum(ConvDial(name)); if (DuneRwtEnum == kNuSystCVResponse) { fUseCV = true; } EnabledParams.push_back({DuneRwtEnum, startval}); } void nusystematicsWeightEngine::SetDialValue(int nuisenum, double val) { systtools::paramId_t DuneRwtEnum = (nuisenum % NUIS_DIAL_OFFSET); if (DuneRwtEnum == kNuSystCVResponse) { return; } systtools::ParamValue &pval = GetParamElementFromContainer(EnabledParams, DuneRwtEnum); fHasChanged = (pval.val - val) > std::numeric_limits::epsilon(); pval.val = val; } void nusystematicsWeightEngine::SetDialValue(std::string name, double val) { if (!IsDialIncluded(name)) { NUIS_ABORT("nusystematicsWeightEngine passed dial: " << name << " that is not enabled."); } systtools::paramId_t DuneRwtEnum(ConvDial(name)); if (DuneRwtEnum == kNuSystCVResponse) { return; } systtools::ParamValue &pval = GetParamElementFromContainer(EnabledParams, DuneRwtEnum); fHasChanged = (pval.val - val) > std::numeric_limits::epsilon(); pval.val = val; } bool nusystematicsWeightEngine::IsDialIncluded(std::string name) { return IsDialIncluded(ConvDial(name)); } bool nusystematicsWeightEngine::IsDialIncluded(int nuisenum) { systtools::paramId_t DuneRwtEnum = (nuisenum % NUIS_DIAL_OFFSET); if (DuneRwtEnum == kNuSystCVResponse) { return fUseCV; } return systtools::ContainterHasParam(EnabledParams, DuneRwtEnum); } double nusystematicsWeightEngine::GetDialValue(std::string name) { if (!IsDialIncluded(name)) { NUIS_ABORT("nusystematicsWeightEngine passed dial: " << name << " that is not enabled."); } systtools::ParamValue &pval = GetParamElementFromContainer(EnabledParams, ConvDial(name)); return pval.val; } double nusystematicsWeightEngine::GetDialValue(int nuisenum) { if (!IsDialIncluded(nuisenum)) { NUIS_ABORT("nusystematicsWeightEngine passed dial: " << nuisenum << " that is not enabled."); } systtools::paramId_t DuneRwtEnum = (nuisenum % NUIS_DIAL_OFFSET); if (DuneRwtEnum == kNuSystCVResponse) { return 1; } systtools::ParamValue &pval = GetParamElementFromContainer(EnabledParams, DuneRwtEnum); return pval.val; } void nusystematicsWeightEngine::Reconfigure(bool silent) { fHasChanged = false; }; bool nusystematicsWeightEngine::NeedsEventReWeight() { if (fHasChanged) { return true; } return false; } double nusystematicsWeightEngine::CalcWeight(BaseFitEvt *evt) { systtools::event_unit_response_w_cv_t responses = DUNErwt.GetEventVariationAndCVResponse(*evt->genie_event->event); double weight = 1; for (auto const &resp : responses) { if (!DUNErwt.IsWeightResponse(resp.pid)) { continue; } if (fUseCV) { weight *= resp.CV_response; } else { // This is very inefficient for fitting, as it recalculates the // spline every time. - systtools::ParamValue const &pval = - GetParamElementFromContainer(EnabledParams, resp.pid); + + //this is a completely backwards way of doing this loop, but the whole thing is broken anyway. + size_t index = GetParamContainerIndex(EnabledParams, resp.pid); + if(index != systtools::kParamUnhandled){ weight *= (resp.CV_response * - DUNErwt.GetParameterResponse(resp.pid, pval.val, + DUNErwt.GetParameterResponse(resp.pid, EnabledParams[index].val, systtools::event_unit_response_t{ {resp.pid, resp.responses}})); + } } } return weight; } void nusystematicsWeightEngine::Print() { std::cout << "nusystematicsWeightEngine: " << std::endl; }