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;
}