diff --git a/src/InputHandler/GENIEInputHandler.cxx b/src/InputHandler/GENIEInputHandler.cxx
index aa9a02a..3ed0af4 100644
--- a/src/InputHandler/GENIEInputHandler.cxx
+++ b/src/InputHandler/GENIEInputHandler.cxx
@@ -1,543 +1,545 @@
// Copyright 2016 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 .
*******************************************************************************/
#ifdef __GENIE_ENABLED__
#include "GENIEInputHandler.h"
#include "InputUtils.h"
#ifdef __DUNERWT_ENABLED__
#include "systematicstools/utility/ParameterAndProviderConfigurationUtility.hh"
#include "fhiclcpp/make_ParameterSet.h"
#endif
GENIEGeneratorInfo::~GENIEGeneratorInfo() { DeallocateParticleStack(); }
void GENIEGeneratorInfo::AddBranchesToTree(TTree *tn) {
tn->Branch("GenieParticlePDGs", &fGenieParticlePDGs, "GenieParticlePDGs/I");
}
void GENIEGeneratorInfo::SetBranchesFromTree(TTree *tn) {
tn->SetBranchAddress("GenieParticlePDGs", &fGenieParticlePDGs);
}
void GENIEGeneratorInfo::AllocateParticleStack(int stacksize) {
fGenieParticlePDGs = new int[stacksize];
}
void GENIEGeneratorInfo::DeallocateParticleStack() {
delete fGenieParticlePDGs;
}
void GENIEGeneratorInfo::FillGeneratorInfo(NtpMCEventRecord *ntpl) {
Reset();
// Check for GENIE Event
if (!ntpl)
return;
if (!ntpl->event)
return;
// Cast Event Record
GHepRecord *ghep = static_cast(ntpl->event);
if (!ghep)
return;
// Fill Particle Stack
GHepParticle *p = 0;
TObjArrayIter iter(ghep);
// Loop over all particles
int i = 0;
while ((p = (dynamic_cast((iter).Next())))) {
if (!p)
continue;
// Get PDG
fGenieParticlePDGs[i] = p->Pdg();
i++;
}
}
void GENIEGeneratorInfo::Reset() {
for (int i = 0; i < kMaxParticles; i++) {
fGenieParticlePDGs[i] = 0;
}
}
GENIEInputHandler::GENIEInputHandler(std::string const &handle,
std::string const &rawinputs) {
LOG(SAM) << "Creating GENIEInputHandler : " << handle << std::endl;
// Run a joint input handling
fName = handle;
// Setup the TChain
fGENIETree = new TChain("gtree");
fSaveExtra = FitPar::Config().GetParB("SaveExtraGenie");
fCacheSize = FitPar::Config().GetParI("CacheSize");
fMaxEvents = FitPar::Config().GetParI("MAXEVENTS");
// Loop over all inputs and grab flux, eventhist, and nevents
std::vector inputs = InputUtils::ParseInputFileList(rawinputs);
for (size_t inp_it = 0; inp_it < inputs.size(); ++inp_it) {
// Open File for histogram access
TFile *inp_file = new TFile(
InputUtils::ExpandInputDirectories(inputs[inp_it]).c_str(), "READ");
if (!inp_file or inp_file->IsZombie()) {
THROW("GENIE File IsZombie() at : '"
<< inputs[inp_it] << "'" << std::endl
<< "Check that your file paths are correct and the file exists!"
<< std::endl
<< "$ ls -lh " << inputs[inp_it]);
}
// Get Flux/Event hist
TH1D *fluxhist = (TH1D *)inp_file->Get("nuisance_flux");
TH1D *eventhist = (TH1D *)inp_file->Get("nuisance_events");
if (!fluxhist or !eventhist) {
ERROR(FTL, "Input File Contents: " << inputs[inp_it]);
inp_file->ls();
THROW("GENIE FILE doesn't contain flux/xsec info."
<< std::endl
<< "Try running the app PrepareGENIE first on :" << inputs[inp_it]
<< std::endl
<< "$ PrepareGENIE -h");
}
// Get N Events
TTree *genietree = (TTree *)inp_file->Get("gtree");
if (!genietree) {
ERROR(FTL, "gtree not located in GENIE file: " << inputs[inp_it]);
THROW("Check your inputs, they may need to be completely regenerated!");
throw;
}
int nevents = genietree->GetEntries();
if (nevents <= 0) {
THROW("Trying to a TTree with "
<< nevents << " to TChain from : " << inputs[inp_it]);
}
// Register input to form flux/event rate hists
RegisterJointInput(inputs[inp_it], nevents, fluxhist, eventhist);
// Add To TChain
fGENIETree->AddFile(inputs[inp_it].c_str());
}
// Registor all our file inputs
SetupJointInputs();
// Assign to tree
fEventType = kGENIE;
fGenieNtpl = NULL;
fGENIETree->SetBranchAddress("gmcrec", &fGenieNtpl);
// Libraries should be seen but not heard...
StopTalking();
fGENIETree->GetEntry(0);
StartTalking();
#ifndef __DUNERWT_ENABLED__
// Create Fit Event
fNUISANCEEvent = new FitEvent();
fNUISANCEEvent->SetGenieEvent(fGenieNtpl);
if (fSaveExtra) {
fGenieInfo = new GENIEGeneratorInfo();
fNUISANCEEvent->AddGeneratorInfo(fGenieInfo);
}
fNUISANCEEvent->HardReset();
#else
std::vector HandlerOpts = Config::QueryKeys("GENIEInputHandler");
fUseCache = HandlerOpts.size() && HandlerOpts.front().Has("UseCache") &&
HandlerOpts.front().GetB("UseCache");
DUNERwtCachedResponseReader = nullptr;
HaveCachedResponseReader = false;
if (fUseCache && (inputs.size() == 1)) {
std::vector DuneRwtCacheParams =
Config::QueryKeys("DUNERwtResponseCache");
for (nuiskey &key : DuneRwtCacheParams) {
if (key.Has("Input") && (key.GetS("Input") == inputs.front()) &&
key.Has("CacheFile") && key.Has("ParameterFHiCL")) {
fhicl::ParameterSet ps =
fhicl::make_ParameterSet(key.GetS("ParameterFHiCL"));
fhicl::ParameterSet syst_providers = ps.get(
"generated_systematic_provider_configuration");
systtools::param_header_map_t configuredParameterHeaders =
systtools::BuildParameterHeaders(syst_providers);
DUNERwtCachedResponseReader =
std::make_unique>(
key.GetS("CacheFile"), "resp_tree",
configuredParameterHeaders.size());
HaveCachedResponseReader = true;
break;
}
}
} else {
fNUISANCEEvent = new FitEvent();
fNUISANCEEvent->SetGenieEvent(fGenieNtpl);
if (fSaveExtra) {
fGenieInfo = new GENIEGeneratorInfo();
fNUISANCEEvent->AddGeneratorInfo(fGenieInfo);
}
fNUISANCEEvent->HardReset();
}
#endif
};
GENIEInputHandler::~GENIEInputHandler() {
// if (fGenieGHep) delete fGenieGHep;
// if (fGenieNtpl) delete fGenieNtpl;
// if (fGENIETree) delete fGENIETree;
// if (fGenieInfo) delete fGenieInfo;
}
void GENIEInputHandler::CreateCache() {
if (fCacheSize > 0) {
// fGENIETree->SetCacheEntryRange(0, fNEvents);
fGENIETree->AddBranchToCache("*", 1);
fGENIETree->SetCacheSize(fCacheSize);
}
}
void GENIEInputHandler::RemoveCache() {
// fGENIETree->SetCacheEntryRange(0, fNEvents);
fGENIETree->AddBranchToCache("*", 0);
fGENIETree->SetCacheSize(0);
}
FitEvent *GENIEInputHandler::GetNuisanceEvent(const UInt_t entry,
const bool lightweight) {
if (entry >= (UInt_t)fNEvents)
return NULL;
+#ifdef __DUNERWT_ENABLED__
// Reduce memory pressure from the cache by clearing out the last entry each
// time.
if (entry && rwEvs[entry - 1].NParticles()) {
rwEvs[entry - 1].DeallocateParticleStack();
}
+#endif
// Read Entry from TTree to fill NEUT Vect in BaseFitEvt;
fGENIETree->GetEntry(entry);
#ifdef __DUNERWT_ENABLED__
if (entry >= rwEvs.size()) {
rwEvs.push_back(FitEvent());
if (HaveCachedResponseReader) {
rwEvs.back().DUNERwtPolyResponses =
DUNERwtCachedResponseReader->GetEventResponse(entry);
rwEvs.back().HasDUNERwtPolyResponses = true;
}
}
rwEvs[entry].SetGenieEvent(fGenieNtpl);
fNUISANCEEvent = &rwEvs[entry];
#endif
// Run NUISANCE Vector Filler
if (!lightweight) {
CalcNUISANCEKinematics();
}
#ifdef __PROB3PP_ENABLED__
else {
// Check for GENIE Event
if (!fGenieNtpl)
return NULL;
if (!fGenieNtpl->event)
return NULL;
// Cast Event Record
fGenieGHep = fGenieNtpl->event;
if (!fGenieGHep)
return NULL;
TObjArrayIter iter(fGenieGHep);
genie::GHepParticle *p;
while ((p = (dynamic_cast((iter).Next())))) {
if (!p) {
continue;
}
// Get Status
int state = GetGENIEParticleStatus(p, fNUISANCEEvent->Mode);
if (state != genie::kIStInitialState) {
continue;
}
fNUISANCEEvent->probe_E = p->E() * 1.E3;
fNUISANCEEvent->probe_pdg = p->Pdg();
break;
}
}
#endif
// Setup Input scaling for joint inputs
fNUISANCEEvent->InputWeight = GetInputWeight(entry);
return fNUISANCEEvent;
}
int GENIEInputHandler::GetGENIEParticleStatus(genie::GHepParticle *p,
int mode) {
/*
kIStUndefined = -1,
kIStInitialState = 0, / generator-level initial state /
kIStStableFinalState = 1, / generator-level final state:
particles to be tracked by detector-level MC /
kIStIntermediateState = 2,
kIStDecayedState = 3,
kIStCorrelatedNucleon = 10,
kIStNucleonTarget = 11,
kIStDISPreFragmHadronicState = 12,
kIStPreDecayResonantState = 13,
kIStHadronInTheNucleus = 14, / hadrons inside the nucleus:
marked for hadron transport modules to act on /
kIStFinalStateNuclearRemnant = 15, / low energy nuclear fragments
entering the record collectively as a 'hadronic blob' pseudo-particle /
kIStNucleonClusterTarget = 16, // for composite nucleons before
phase space decay
*/
int state = kUndefinedState;
switch (p->Status()) {
case genie::kIStNucleonTarget:
case genie::kIStInitialState:
case genie::kIStCorrelatedNucleon:
case genie::kIStNucleonClusterTarget:
state = kInitialState;
break;
case genie::kIStStableFinalState:
state = kFinalState;
break;
case genie::kIStHadronInTheNucleus:
if (abs(mode) == 2)
state = kInitialState;
else
state = kFSIState;
break;
case genie::kIStPreDecayResonantState:
case genie::kIStDISPreFragmHadronicState:
case genie::kIStIntermediateState:
state = kFSIState;
break;
case genie::kIStFinalStateNuclearRemnant:
case genie::kIStUndefined:
case genie::kIStDecayedState:
default:
break;
}
// Flag to remove nuclear part in genie
if (p->Pdg() > 1000000) {
if (state == kInitialState)
state = kNuclearInitial;
else if (state == kFinalState)
state = kNuclearRemnant;
}
return state;
}
#endif
#ifdef __GENIE_ENABLED__
int GENIEInputHandler::ConvertGENIEReactionCode(GHepRecord *gheprec) {
// Electron Scattering
if (gheprec->Summary()->ProcInfo().IsEM()) {
if (gheprec->Summary()->InitState().ProbePdg() == 11) {
if (gheprec->Summary()->ProcInfo().IsQuasiElastic())
return 1;
else if (gheprec->Summary()->ProcInfo().IsMEC())
return 2;
else if (gheprec->Summary()->ProcInfo().IsResonant())
return 13;
else if (gheprec->Summary()->ProcInfo().IsDeepInelastic())
return 26;
else {
ERROR(WRN,
"Unknown GENIE Electron Scattering Mode!"
<< std::endl
<< "ScatteringTypeId = "
<< gheprec->Summary()->ProcInfo().ScatteringTypeId() << " "
<< "InteractionTypeId = "
<< gheprec->Summary()->ProcInfo().InteractionTypeId()
<< std::endl
<< genie::ScatteringType::AsString(
gheprec->Summary()->ProcInfo().ScatteringTypeId())
<< " "
<< genie::InteractionType::AsString(
gheprec->Summary()->ProcInfo().InteractionTypeId())
<< " " << gheprec->Summary()->ProcInfo().IsMEC());
return 0;
}
}
// Weak CC
} else if (gheprec->Summary()->ProcInfo().IsWeakCC()) {
// CC MEC
if (gheprec->Summary()->ProcInfo().IsMEC()) {
if (pdg::IsNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return 2;
else if (pdg::IsAntiNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return -2;
// CC OTHER
} else {
return utils::ghep::NeutReactionCode(gheprec);
}
// Weak NC
} else if (gheprec->Summary()->ProcInfo().IsWeakNC()) {
// NC MEC
if (gheprec->Summary()->ProcInfo().IsMEC()) {
if (pdg::IsNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return 32;
else if (pdg::IsAntiNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return -32;
// NC OTHER
} else {
return utils::ghep::NeutReactionCode(gheprec);
}
}
return 0;
}
void GENIEInputHandler::CalcNUISANCEKinematics() {
// Reset all variables
fNUISANCEEvent->ResetEvent();
// Check for GENIE Event
if (!fGenieNtpl)
return;
if (!fGenieNtpl->event)
return;
// Cast Event Record
fGenieGHep = fGenieNtpl->event;
if (!fGenieGHep)
return;
// Convert GENIE Reaction Code
fNUISANCEEvent->Mode = ConvertGENIEReactionCode(fGenieGHep);
// Set Event Info
fNUISANCEEvent->fEventNo = 0.0;
fNUISANCEEvent->fTotCrs = fGenieGHep->XSec();
fNUISANCEEvent->fTargetA = 0.0;
fNUISANCEEvent->fTargetZ = 0.0;
fNUISANCEEvent->fTargetH = 0;
fNUISANCEEvent->fBound = 0.0;
fNUISANCEEvent->InputWeight =
1.0; //(1E+38 / genie::units::cm2) * fGenieGHep->XSec();
// Get N Particle Stack
unsigned int npart = fGenieGHep->GetEntries();
unsigned int kmax = fNUISANCEEvent->kMaxParticles;
if (npart > kmax) {
fNUISANCEEvent->ExpandParticleStack(npart);
}
// Fill Particle Stack
GHepParticle *p = 0;
TObjArrayIter iter(fGenieGHep);
fNUISANCEEvent->fNParticles = 0;
// Loop over all particles
while ((p = (dynamic_cast((iter).Next())))) {
if (!p)
continue;
// Get Status
int state = GetGENIEParticleStatus(p, fNUISANCEEvent->Mode);
// Remove Undefined
if (kRemoveUndefParticles && state == kUndefinedState)
continue;
// Remove FSI
if (kRemoveFSIParticles && state == kFSIState)
continue;
if (kRemoveNuclearParticles &&
(state == kNuclearInitial || state == kNuclearRemnant))
continue;
// Fill Vectors
int curpart = fNUISANCEEvent->fNParticles;
fNUISANCEEvent->fParticleState[curpart] = state;
// Mom
fNUISANCEEvent->fParticleMom[curpart][0] = p->Px() * 1.E3;
fNUISANCEEvent->fParticleMom[curpart][1] = p->Py() * 1.E3;
fNUISANCEEvent->fParticleMom[curpart][2] = p->Pz() * 1.E3;
fNUISANCEEvent->fParticleMom[curpart][3] = p->E() * 1.E3;
// PDG
fNUISANCEEvent->fParticlePDG[curpart] = p->Pdg();
// Add to N particle count
fNUISANCEEvent->fNParticles++;
// Extra Check incase GENIE fails.
if ((UInt_t)fNUISANCEEvent->fNParticles == kmax) {
ERR(WRN) << "Number of GENIE Particles exceeds maximum (Max: " << kmax
<< ", GHEP: " << fGenieGHep->GetEntries()
<< ", Added: " << fNUISANCEEvent->fNParticles << ")!"
<< std::endl;
ERR(WRN) << "Extend kMax, or run without including FSI particles!"
<< std::endl;
break;
}
}
// Fill Extra Stack
if (fSaveExtra)
fGenieInfo->FillGeneratorInfo(fGenieNtpl);
// Run Initial, FSI, Final, Other ordering.
fNUISANCEEvent->OrderStack();
FitParticle *ISNeutralLepton =
fNUISANCEEvent->GetHMISParticle(PhysConst::pdg_neutrinos);
if (ISNeutralLepton) {
fNUISANCEEvent->probe_E = ISNeutralLepton->E();
fNUISANCEEvent->probe_pdg = ISNeutralLepton->PDG();
}
return;
}
void GENIEInputHandler::Print() {}
#endif
diff --git a/src/Reweight/DUNERwtWeightEngine.cxx b/src/Reweight/DUNERwtWeightEngine.cxx
index 00f5d3d..a660446 100644
--- a/src/Reweight/DUNERwtWeightEngine.cxx
+++ b/src/Reweight/DUNERwtWeightEngine.cxx
@@ -1,164 +1,165 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
-
+#ifdef __DUNERWT_ENABLED__
/*******************************************************************************
* 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 .
*******************************************************************************/
#include "DUNERwtWeightEngine.h"
#include
#include
DUNERwtWeightEngine::DUNERwtWeightEngine() { Config(); }
void DUNERwtWeightEngine::Config() {
std::vector DuneRwtParam = Config::QueryKeys("DUNERwt");
if (DuneRwtParam.size() < 1) {
ERROR(WRN, "Oscillation parameters specified but no OscParam element "
"configuring the experimental characteristics found.\nExpect at "
"least . Pausing for "
"10...");
sleep(10);
return;
}
std::string fhicl_name = DuneRwtParam.front().GetS("ConfigFHiCL");
DUNErwt.LoadConfiguration(fhicl_name);
}
int DUNERwtWeightEngine::ConvDial(std::string name) {
if (!DUNErwt.HaveHeader(name)) {
THROW("DUNERwtWeightEngine passed dial: "
<< name << " that it does not understand.");
}
return DUNErwt.GetHeaderId(name);
}
void DUNERwtWeightEngine::IncludeDial(std::string name, double startval) {
EnabledParams.push_back({systtools::paramId_t(ConvDial(name)), startval});
}
void DUNERwtWeightEngine::SetDialValue(int nuisenum, double val) {
systtools::paramId_t DuneRwtEnum = (nuisenum % 1000);
systtools::ParamValue &pval =
GetParamElementFromContainer(EnabledParams, DuneRwtEnum);
fHasChanged = (pval.val - val) > std::numeric_limits::epsilon();
pval.val = val;
}
void DUNERwtWeightEngine::SetDialValue(std::string name, double val) {
if (!IsDialIncluded(name)) {
THROW("DUNERwtWeightEngine passed dial: " << name
<< " that is not enabled.");
}
systtools::ParamValue &pval =
GetParamElementFromContainer(EnabledParams, ConvDial(name));
fHasChanged = (pval.val - val) > std::numeric_limits::epsilon();
pval.val = val;
}
bool DUNERwtWeightEngine::IsDialIncluded(std::string name) {
return IsDialIncluded(ConvDial(name));
}
bool DUNERwtWeightEngine::IsDialIncluded(int nuisenum) {
systtools::paramId_t DuneRwtEnum = (nuisenum % 1000);
return systtools::ContainterHasParam(EnabledParams, DuneRwtEnum);
}
double DUNERwtWeightEngine::GetDialValue(std::string name) {
if (!IsDialIncluded(name)) {
THROW("DUNERwtWeightEngine passed dial: " << name
<< " that is not enabled.");
}
systtools::ParamValue &pval =
GetParamElementFromContainer(EnabledParams, ConvDial(name));
return pval.val;
}
double DUNERwtWeightEngine::GetDialValue(int nuisenum) {
if (!IsDialIncluded(nuisenum)) {
THROW("DUNERwtWeightEngine passed dial: " << nuisenum
<< " that is not enabled.");
}
systtools::paramId_t DuneRwtEnum = (nuisenum % 1000);
systtools::ParamValue &pval =
GetParamElementFromContainer(EnabledParams, DuneRwtEnum);
return pval.val;
}
void DUNERwtWeightEngine::Reconfigure(bool silent) { fHasChanged = false; };
bool DUNERwtWeightEngine::NeedsEventReWeight() {
if (fHasChanged) {
return true;
}
return false;
}
double DUNERwtWeightEngine::CalcWeight(BaseFitEvt *evt) {
double weight = 1;
if (evt->HasDUNERwtPolyResponses) {
for (size_t i = 0; i < EnabledParams.size(); ++i) {
if (DUNErwt.IsSplineParam(EnabledParams[i].pid)) {
if (!ContainterHasParam(evt->DUNERwtPolyResponses,
EnabledParams[i].pid)) {
continue;
}
weight *= GetParamElementFromContainer(evt->DUNERwtPolyResponses,
EnabledParams[i].pid)
.resp.eval(EnabledParams[i].val);
} else {
if (!evt->HasDUNERwtResponses) {
evt->DUNERwtResponses =
DUNErwt.GetEventResponses(*evt->genie_event->event);
evt->HasDUNERwtResponses = true;
}
weight *= DUNErwt.GetDiscreteResponse(EnabledParams[i].pid,
int(EnabledParams[i].val),
evt->DUNERwtResponses);
}
}
} else {
if (!evt->HasDUNERwtResponses) {
evt->DUNERwtResponses =
DUNErwt.GetEventResponses(*evt->genie_event->event);
evt->HasDUNERwtResponses = true;
}
for (size_t i = 0; i < EnabledParams.size(); ++i) {
if (DUNErwt.IsSplineParam(EnabledParams[i].pid)) {
weight *= DUNErwt.GetParameterResponse(
EnabledParams[i].pid, EnabledParams[i].val, evt->DUNERwtResponses);
} else {
weight *= DUNErwt.GetDiscreteResponse(EnabledParams[i].pid,
int(EnabledParams[i].val),
evt->DUNERwtResponses);
}
}
}
return weight;
}
void DUNERwtWeightEngine::Print() {
std::cout << "DUNERwtWeightEngine: " << std::endl;
}
+#endif
diff --git a/src/Reweight/DUNERwtWeightEngine.h b/src/Reweight/DUNERwtWeightEngine.h
index 908f1d3..5c73f5a 100644
--- a/src/Reweight/DUNERwtWeightEngine.h
+++ b/src/Reweight/DUNERwtWeightEngine.h
@@ -1,64 +1,66 @@
#ifndef DUNERWTWEIGHTENGINE_SEEN
#define DUNERWTWEIGHTENGINE_SEEN
+#ifdef __DUNERWT_ENABLED__
// Copyright 2016 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 .
*******************************************************************************/
#include "WeightEngineBase.h"
#include "systematicstools/interface/types.hh"
#include "nusystematics/artless/response_helper.hh"
#include
class DUNERwtWeightEngine : public WeightEngineBase {
public:
DUNERwtWeightEngine();
nusyst::response_helper DUNErwt;
systtools::param_value_list_t EnabledParams;
void Config();
int ConvDial(std::string name);
// 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);
bool NeedsEventReWeight();
double CalcWeight(BaseFitEvt* evt);
void Print();
};
#endif
+#endif
diff --git a/src/Reweight/FitWeight.cxx b/src/Reweight/FitWeight.cxx
index 436f333..cf33305 100644
--- a/src/Reweight/FitWeight.cxx
+++ b/src/Reweight/FitWeight.cxx
@@ -1,293 +1,300 @@
#include "FitWeight.h"
+#ifdef __DUNERWT_ENABLED__
#include "DUNERwtWeightEngine.h"
+#endif
#include "GENIEWeightEngine.h"
#include "LikelihoodWeightEngine.h"
#include "ModeNormEngine.h"
#include "NEUTWeightEngine.h"
#include "NIWGWeightEngine.h"
#include "NUISANCEWeightEngine.h"
#include "NuWroWeightEngine.h"
#include "OscWeightEngine.h"
#include "SampleNormEngine.h"
#include "SplineWeightEngine.h"
#include "T2KWeightEngine.h"
void FitWeight::AddRWEngine(int type) {
switch (type) {
case kNEUT:
fAllRW[type] = new NEUTWeightEngine("neutrw");
break;
case kNUWRO:
fAllRW[type] = new NuWroWeightEngine("nuwrorw");
break;
case kGENIE:
fAllRW[type] = new GENIEWeightEngine("genierw");
break;
case kNORM:
fAllRW[type] = new SampleNormEngine("normrw");
break;
case kLIKEWEIGHT:
fAllRW[type] = new LikelihoodWeightEngine("likerw");
break;
case kT2K:
fAllRW[type] = new T2KWeightEngine("t2krw");
break;
case kCUSTOM:
fAllRW[type] = new NUISANCEWeightEngine("nuisrw");
break;
case kSPLINEPARAMETER:
fAllRW[type] = new SplineWeightEngine("splinerw");
break;
case kNIWG:
fAllRW[type] = new NIWGWeightEngine("niwgrw");
break;
case kOSCILLATION:
fAllRW[type] = new OscWeightEngine();
break;
case kMODENORM:
fAllRW[type] = new ModeNormEngine();
break;
+#ifdef __DUNERWT_ENABLED__
case kDUNERwt:
fAllRW[type] = new DUNERwtWeightEngine();
break;
+#endif
default:
THROW("CANNOT ADD RW Engine for unknown dial type: " << type);
break;
}
}
WeightEngineBase *FitWeight::GetRWEngine(int type) {
if (HasRWEngine(type)) {
return fAllRW[type];
}
THROW("CANNOT get RW Engine for dial type: " << type);
}
bool FitWeight::HasRWEngine(int type) {
switch (type) {
case kNEUT:
case kNUWRO:
case kGENIE:
case kNORM:
case kLIKEWEIGHT:
case kT2K:
case kCUSTOM:
case kSPLINEPARAMETER:
case kNIWG:
case kOSCILLATION:
- case kDUNERwt: {
+#ifdef __DUNERWT_ENABLED__
+ case kDUNERwt:
+#endif
+ {
return fAllRW.count(type);
}
default: { THROW("CANNOT get RW Engine for dial type: " << type); }
}
}
void FitWeight::IncludeDial(std::string name, std::string type, double val) {
// Should register the dial here.
int typeenum = Reweight::ConvDialType(type);
IncludeDial(name, typeenum, val);
}
void FitWeight::IncludeDial(std::string name, int dialtype, double val) {
// Setup RW Engine Pointer
if (fAllRW.find(dialtype) == fAllRW.end()) {
AddRWEngine(dialtype);
}
// Get the dial enum
int nuisenum;
if (dialtype != kDUNERwt) {
nuisenum = Reweight::ConvDial(name, dialtype);
} else {
#ifdef __DUNERWT_ENABLED__
DUNERwtWeightEngine *drw =
static_cast(fAllRW[kDUNERwt]);
- nuisenum = drw->ConvDial(name) + dialtype*1000;
+ nuisenum = drw->ConvDial(name) + dialtype * 1000;
#else
- THROW("[ERROR]: Not built with DUNEReWeight support but ")
+ THROW("[ERROR]: Not built with DUNEReWeight support.");
#endif
}
if (nuisenum == -1) {
THROW("NUISENUM == " << nuisenum << " for " << name);
}
WeightEngineBase *rw = fAllRW[dialtype];
// Include the dial
rw->IncludeDial(name, val);
// Set Dial Value
if (val != -9999.9) {
rw->SetDialValue(name, val);
}
// Sort Maps
fAllEnums[name] = nuisenum;
fAllValues[nuisenum] = val;
// Sort Lists
fNameList.push_back(name);
fEnumList.push_back(nuisenum);
fValueList.push_back(val);
}
void FitWeight::Reconfigure(bool silent) {
// Reconfigure all added RW engines
for (std::map::iterator iter = fAllRW.begin();
iter != fAllRW.end(); iter++) {
(*iter).second->Reconfigure(silent);
}
}
void FitWeight::SetDialValue(std::string name, double val) {
// Add extra check, if name not found look for one with name in it.
int nuisenum = fAllEnums[name];
SetDialValue(nuisenum, val);
}
// Allow for name aswell using GlobalList to determine sample name.
void FitWeight::SetDialValue(int nuisenum, double val) {
// Conv dial type
int dialtype = Reweight::GetDialType(nuisenum);
if (fAllRW.find(dialtype) == fAllRW.end()) {
THROW("Cannot find RW Engine for dialtype = "
<< dialtype << ", " << Reweight::RemoveDialType(nuisenum));
}
// Get RW Engine for this dial
fAllRW[dialtype]->SetDialValue(nuisenum, val);
fAllValues[nuisenum] = val;
// Update ValueList
for (size_t i = 0; i < fEnumList.size(); i++) {
if (fEnumList[i] == nuisenum) {
fValueList[i] = val;
}
}
}
void FitWeight::SetAllDials(const double *x, int n) {
for (size_t i = 0; i < (UInt_t)n; i++) {
int rwenum = fEnumList[i];
SetDialValue(rwenum, x[i]);
}
Reconfigure();
}
double FitWeight::GetDialValue(std::string name) {
// Add extra check, if name not found look for one with name in it.
int nuisenum = fAllEnums[name];
return GetDialValue(nuisenum);
}
double FitWeight::GetDialValue(int nuisenum) { return fAllValues[nuisenum]; }
int FitWeight::GetDialPos(std::string name) {
int rwenum = fAllEnums[name];
return GetDialPos(rwenum);
}
int FitWeight::GetDialPos(int nuisenum) {
for (size_t i = 0; i < fEnumList.size(); i++) {
if (fEnumList[i] == nuisenum) {
return i;
}
}
ERR(FTL) << "No Dial Found! " << std::endl;
throw;
return -1;
}
bool FitWeight::DialIncluded(std::string name) {
return (fAllEnums.find(name) != fAllEnums.end());
}
bool FitWeight::DialIncluded(int rwenum) {
return (fAllValues.find(rwenum) != fAllValues.end());
}
double FitWeight::CalcWeight(BaseFitEvt *evt) {
double rwweight = 1.0;
for (std::map::iterator iter = fAllRW.begin();
iter != fAllRW.end(); iter++) {
double w = (*iter).second->CalcWeight(evt);
rwweight *= w;
}
return rwweight;
}
void FitWeight::UpdateWeightEngine(const double *x) {
size_t count = 0;
for (std::vector::iterator iter = fEnumList.begin();
iter != fEnumList.end(); iter++) {
SetDialValue((*iter), x[count]);
count++;
}
}
void FitWeight::GetAllDials(double *x, int n) {
for (int i = 0; i < n; i++) {
x[i] = GetDialValue(fEnumList[i]);
}
}
// bool FitWeight::NeedsEventReWeight(const double* x) {
// bool haschange = false;
// size_t count = 0;
// // Compare old to new and decide if RW needed.
// for (std::vector::iterator iter = fEnumList.begin();
// iter != fEnumList.end(); iter++) {
// int nuisenum = (*iter);
// int type = (nuisenum / 1000) - (nuisenum % 1000);
// // Compare old to new
// double oldval = GetDialValue(nuisenum);
// double newval = x[count];
// if (oldval != newval) {
// if (fAllRW[type]->NeedsEventReWeight()) {
// haschange = true;
// }
// }
// count++;
// }
// return haschange;
// }
double FitWeight::GetSampleNorm(std::string name) {
if (name.empty())
return 1.0;
// Find norm dial
if (fAllEnums.find(name + "_norm") != fAllEnums.end()) {
if (fAllValues.find(fAllEnums[name + "_norm"]) != fAllValues.end()) {
return fAllValues[fAllEnums[name + "_norm"]];
} else {
return 1.0;
}
} else {
return 1.0;
}
}
void FitWeight::Print() {
LOG(REC) << "Fit Weight State: " << std::endl;
for (size_t i = 0; i < fNameList.size(); i++) {
LOG(REC) << " -> Par " << i << ". " << fNameList[i] << " " << fValueList[i]
<< std::endl;
}
}
diff --git a/src/Reweight/GENIEWeightEngine.cxx b/src/Reweight/GENIEWeightEngine.cxx
index 4586ae4..b6dbd5a 100644
--- a/src/Reweight/GENIEWeightEngine.cxx
+++ b/src/Reweight/GENIEWeightEngine.cxx
@@ -1,244 +1,260 @@
#include "GENIEWeightEngine.h"
+#ifdef __GENIE_EMP_MECRW_ENABLED
+#include "ReWeight/GReWeightXSecEmpiricalMEC.h"
+#endif
+
GENIEWeightEngine::GENIEWeightEngine(std::string name) {
#ifdef __GENIE_ENABLED__
- // Setup the NEUT Reweight engien
- fCalcName = name;
- LOG(FIT) << "Setting up GENIE RW : " << fCalcName << std::endl;
-
- // Create RW Engine suppressing cout
- StopTalking();
- fGenieRW = new genie::rew::GReWeight();
-
- // Get List of Vetos (Just for debugging)
- std::string rw_engine_list = FitPar::Config().GetParS("FitWeight_fGenieRW_veto");
- bool xsec_ncel = rw_engine_list.find("xsec_ncel") == std::string::npos;
- bool xsec_ccqe = rw_engine_list.find("xsec_ccqe") == std::string::npos;
- bool xsec_coh = rw_engine_list.find("xsec_coh") == std::string::npos;
- bool xsec_nnres = rw_engine_list.find("xsec_nonresbkg") == std::string::npos;
- bool xsec_nudis = rw_engine_list.find("nuclear_dis") == std::string::npos;
- bool xsec_resdec = rw_engine_list.find("hadro_res_decay") == std::string::npos;
- bool xsec_fzone = rw_engine_list.find("hadro_intranuke") == std::string::npos;
- bool xsec_intra = rw_engine_list.find("hadro_fzone") == std::string::npos;
- bool xsec_agky = rw_engine_list.find("hadro_agky") == std::string::npos;
- bool xsec_qevec = rw_engine_list.find("xsec_ccqe_vec") == std::string::npos;
- bool xsec_dis = rw_engine_list.find("xsec_dis") == std::string::npos;
- bool xsec_nc = rw_engine_list.find("xsec_nc") == std::string::npos;
- bool xsec_ccres = rw_engine_list.find("xsec_ccres") == std::string::npos;
- bool xsec_ncres = rw_engine_list.find("xsec_ncres") == std::string::npos;
- bool xsec_nucqe = rw_engine_list.find("nuclear_qe") == std::string::npos;
- bool xsec_qeaxial = rw_engine_list.find("xsec_ccqe_axial") == std::string::npos;
-
- // Now actually add the RW Calcs
- if (xsec_ncel)
- fGenieRW->AdoptWghtCalc("xsec_ncel", new genie::rew::GReWeightNuXSecNCEL);
- if (xsec_ccqe){
- fGenieRW->AdoptWghtCalc("xsec_ccqe", new genie::rew::GReWeightNuXSecCCQE);
- // (dynamic_cast (fGenieRW->WghtCalc("xsec_ccqe")))
- // ->SetXSecModel( FitPar::Config().GetParS("GENIEXSecModelCCQE") );
- }
- if (xsec_coh){
- fGenieRW->AdoptWghtCalc("xsec_coh", new genie::rew::GReWeightNuXSecCOH());
- // (dynamic_cast (fGenieRW->WghtCalc("xsec_coh")))
- // ->SetXSecModel( FitPar::Config().GetParS("GENIEXSecModelCOH") );
- }
-
- if (xsec_nnres)
- fGenieRW->AdoptWghtCalc("xsec_nonresbkg",
- new genie::rew::GReWeightNonResonanceBkg);
- if (xsec_nudis)
- fGenieRW->AdoptWghtCalc("nuclear_dis", new genie::rew::GReWeightDISNuclMod);
- if (xsec_resdec)
- fGenieRW->AdoptWghtCalc("hadro_res_decay",
- new genie::rew::GReWeightResonanceDecay);
- if (xsec_fzone)
- fGenieRW->AdoptWghtCalc("hadro_fzone", new genie::rew::GReWeightFZone);
- if (xsec_intra)
- fGenieRW->AdoptWghtCalc("hadro_intranuke", new genie::rew::GReWeightINuke);
- if (xsec_agky)
- fGenieRW->AdoptWghtCalc("hadro_agky", new genie::rew::GReWeightAGKY);
- if (xsec_qevec)
- fGenieRW->AdoptWghtCalc("xsec_ccqe_vec",
- new genie::rew::GReWeightNuXSecCCQEvec);
+ // Setup the NEUT Reweight engien
+ fCalcName = name;
+ LOG(FIT) << "Setting up GENIE RW : " << fCalcName << std::endl;
+
+ // Create RW Engine suppressing cout
+ StopTalking();
+ fGenieRW = new genie::rew::GReWeight();
+
+ // Get List of Vetos (Just for debugging)
+ std::string rw_engine_list =
+ FitPar::Config().GetParS("FitWeight_fGenieRW_veto");
+ bool xsec_ncel = rw_engine_list.find("xsec_ncel") == std::string::npos;
+ bool xsec_ccqe = rw_engine_list.find("xsec_ccqe") == std::string::npos;
+ bool xsec_coh = rw_engine_list.find("xsec_coh") == std::string::npos;
+ bool xsec_nnres = rw_engine_list.find("xsec_nonresbkg") == std::string::npos;
+ bool xsec_nudis = rw_engine_list.find("nuclear_dis") == std::string::npos;
+ bool xsec_resdec =
+ rw_engine_list.find("hadro_res_decay") == std::string::npos;
+ bool xsec_fzone = rw_engine_list.find("hadro_intranuke") == std::string::npos;
+ bool xsec_intra = rw_engine_list.find("hadro_fzone") == std::string::npos;
+ bool xsec_agky = rw_engine_list.find("hadro_agky") == std::string::npos;
+ bool xsec_qevec = rw_engine_list.find("xsec_ccqe_vec") == std::string::npos;
+ bool xsec_dis = rw_engine_list.find("xsec_dis") == std::string::npos;
+ bool xsec_nc = rw_engine_list.find("xsec_nc") == std::string::npos;
+ bool xsec_ccres = rw_engine_list.find("xsec_ccres") == std::string::npos;
+ bool xsec_ncres = rw_engine_list.find("xsec_ncres") == std::string::npos;
+ bool xsec_nucqe = rw_engine_list.find("nuclear_qe") == std::string::npos;
+ bool xsec_qeaxial =
+ rw_engine_list.find("xsec_ccqe_axial") == std::string::npos;
+#ifdef __GENIE_EMP_MECRW_ENABLED
+ bool xsec_empMEC = rw_engine_list.find("xsec_empMEC") == std::string::npos;
+#endif
+
+ // Now actually add the RW Calcs
+ if (xsec_ncel)
+ fGenieRW->AdoptWghtCalc("xsec_ncel", new genie::rew::GReWeightNuXSecNCEL);
+ if (xsec_ccqe) {
+ fGenieRW->AdoptWghtCalc("xsec_ccqe", new genie::rew::GReWeightNuXSecCCQE);
+ // (dynamic_cast (fGenieRW->WghtCalc("xsec_ccqe")))
+ // ->SetXSecModel( FitPar::Config().GetParS("GENIEXSecModelCCQE") );
+ }
+#ifdef __GENIE_EMP_MECRW_ENABLED
+ if (xsec_empMEC) {
+ fGenieRW->AdoptWghtCalc("xsec_empMEC",
+ new genie::rew::GReWeightXSecEmpiricalMEC);
+ }
+#endif
+ if (xsec_coh) {
+ fGenieRW->AdoptWghtCalc("xsec_coh", new genie::rew::GReWeightNuXSecCOH());
+ // (dynamic_cast (fGenieRW->WghtCalc("xsec_coh")))
+ // ->SetXSecModel( FitPar::Config().GetParS("GENIEXSecModelCOH") );
+ }
+
+ if (xsec_nnres)
+ fGenieRW->AdoptWghtCalc("xsec_nonresbkg",
+ new genie::rew::GReWeightNonResonanceBkg);
+ if (xsec_nudis)
+ fGenieRW->AdoptWghtCalc("nuclear_dis", new genie::rew::GReWeightDISNuclMod);
+ if (xsec_resdec)
+ fGenieRW->AdoptWghtCalc("hadro_res_decay",
+ new genie::rew::GReWeightResonanceDecay);
+ if (xsec_fzone)
+ fGenieRW->AdoptWghtCalc("hadro_fzone", new genie::rew::GReWeightFZone);
+ if (xsec_intra)
+ fGenieRW->AdoptWghtCalc("hadro_intranuke", new genie::rew::GReWeightINuke);
+ if (xsec_agky)
+ fGenieRW->AdoptWghtCalc("hadro_agky", new genie::rew::GReWeightAGKY);
+ if (xsec_qevec)
+ fGenieRW->AdoptWghtCalc("xsec_ccqe_vec",
+ new genie::rew::GReWeightNuXSecCCQEvec);
#if __GENIE_VERSION__ >= 212
- if (xsec_qeaxial)
- fGenieRW->AdoptWghtCalc("xsec_ccqe_axial",
- new genie::rew::GReWeightNuXSecCCQEaxial);
+ if (xsec_qeaxial)
+ fGenieRW->AdoptWghtCalc("xsec_ccqe_axial",
+ new genie::rew::GReWeightNuXSecCCQEaxial);
#endif
- if (xsec_dis)
- fGenieRW->AdoptWghtCalc("xsec_dis", new genie::rew::GReWeightNuXSecDIS);
- if (xsec_nc)
- fGenieRW->AdoptWghtCalc("xsec_nc", new genie::rew::GReWeightNuXSecNC);
- if (xsec_ccres){
+ if (xsec_dis)
+ fGenieRW->AdoptWghtCalc("xsec_dis", new genie::rew::GReWeightNuXSecDIS);
+ if (xsec_nc)
+ fGenieRW->AdoptWghtCalc("xsec_nc", new genie::rew::GReWeightNuXSecNC);
+ if (xsec_ccres) {
#if __GENIE_VERSION__ < 213
- fGenieRW->AdoptWghtCalc("xsec_ccres", new genie::rew::GReWeightNuXSecCCRES);
+ fGenieRW->AdoptWghtCalc("xsec_ccres", new genie::rew::GReWeightNuXSecCCRES);
#else
- fGenieRW->AdoptWghtCalc("xsec_ccres", new genie::rew::GReWeightNuXSecCCRES( FitPar::Config().GetParS("GENIEXSecModelCCRES"), "Default"));
+ fGenieRW->AdoptWghtCalc(
+ "xsec_ccres",
+ new genie::rew::GReWeightNuXSecCCRES(
+ FitPar::Config().GetParS("GENIEXSecModelCCRES"), "Default"));
#endif
- }
-
- if (xsec_ncres)
- fGenieRW->AdoptWghtCalc("xsec_ncres", new genie::rew::GReWeightNuXSecNCRES);
- if (xsec_nucqe)
- fGenieRW->AdoptWghtCalc("nuclear_qe", new genie::rew::GReWeightFGM);
-
- GReWeightNuXSecCCQE * rwccqe =
- dynamic_cast (fGenieRW->WghtCalc("xsec_ccqe"));
- rwccqe->SetMode(GReWeightNuXSecCCQE::kModeMa);
-
- // Default to include shape and normalization changes for CCRES (can be changed downstream if desired)
- GReWeightNuXSecCCRES * rwccres =
- dynamic_cast (fGenieRW->WghtCalc("xsec_ccres"));
- std::string marestype = FitPar::Config().GetParS("GENIEWeightEngine_CCRESMode");
- if (!marestype.compare("kModeNormAndMaMvShape")){ rwccres->SetMode(GReWeightNuXSecCCRES::kModeNormAndMaMvShape); }
- else if (!marestype.compare("kModeMaMv")){ rwccres->SetMode(GReWeightNuXSecCCRES::kModeMaMv); }
- else {
- THROW("Unkown MARES Mode in GENIE Weight Engine : " << marestype );
- }
-
- // Default to include shape and normalization changes for NCRES (can be changed downstream if desired)
- GReWeightNuXSecNCRES * rwncres =
- dynamic_cast (fGenieRW->WghtCalc("xsec_ncres"));
- rwncres->SetMode(GReWeightNuXSecNCRES::kModeMaMv);
-
- // Default to include shape and normalization changes for DIS (can be changed downstream if desired)
- GReWeightNuXSecDIS * rwdis =
- dynamic_cast (fGenieRW->WghtCalc("xsec_dis"));
- rwdis->SetMode(GReWeightNuXSecDIS::kModeABCV12u);
-
- // Set Abs Twk Config
- fIsAbsTwk = (FitPar::Config().GetParB("setabstwk"));
-
- // allow cout again
- StartTalking();
+ }
+
+ if (xsec_ncres)
+ fGenieRW->AdoptWghtCalc("xsec_ncres", new genie::rew::GReWeightNuXSecNCRES);
+ if (xsec_nucqe)
+ fGenieRW->AdoptWghtCalc("nuclear_qe", new genie::rew::GReWeightFGM);
+
+ GReWeightNuXSecCCQE *rwccqe =
+ dynamic_cast(fGenieRW->WghtCalc("xsec_ccqe"));
+ rwccqe->SetMode(GReWeightNuXSecCCQE::kModeMa);
+
+ // Default to include shape and normalization changes for CCRES (can be
+ // changed downstream if desired)
+ GReWeightNuXSecCCRES *rwccres =
+ dynamic_cast(fGenieRW->WghtCalc("xsec_ccres"));
+ std::string marestype =
+ FitPar::Config().GetParS("GENIEWeightEngine_CCRESMode");
+ if (!marestype.compare("kModeNormAndMaMvShape")) {
+ rwccres->SetMode(GReWeightNuXSecCCRES::kModeNormAndMaMvShape);
+ } else if (!marestype.compare("kModeMaMv")) {
+ rwccres->SetMode(GReWeightNuXSecCCRES::kModeMaMv);
+ } else {
+ THROW("Unkown MARES Mode in GENIE Weight Engine : " << marestype);
+ }
+
+ // Default to include shape and normalization changes for NCRES (can be
+ // changed downstream if desired)
+ GReWeightNuXSecNCRES *rwncres =
+ dynamic_cast(fGenieRW->WghtCalc("xsec_ncres"));
+ rwncres->SetMode(GReWeightNuXSecNCRES::kModeMaMv);
+
+ // Default to include shape and normalization changes for DIS (can be changed
+ // downstream if desired)
+ GReWeightNuXSecDIS *rwdis =
+ dynamic_cast(fGenieRW->WghtCalc("xsec_dis"));
+ rwdis->SetMode(GReWeightNuXSecDIS::kModeABCV12u);
+
+ // Set Abs Twk Config
+ fIsAbsTwk = (FitPar::Config().GetParB("setabstwk"));
+
+ // allow cout again
+ StartTalking();
#else
- ERR(FTL) << "GENIE RW NOT ENABLED" << std::endl;
+ ERR(FTL) << "GENIE RW NOT ENABLED" << std::endl;
#endif
};
-
void GENIEWeightEngine::IncludeDial(std::string name, double startval) {
#ifdef __GENIE_ENABLED__
- // Get First enum
- int nuisenum = Reweight::ConvDial(name, kGENIE);
-
- // Setup Maps
- fEnumIndex[nuisenum];// = std::vector(0);
- fNameIndex[name]; // = std::vector(0);
+ // Get First enum
+ int nuisenum = Reweight::ConvDial(name, kGENIE);
- // Split by commas
- std::vector allnames = GeneralUtils::ParseToStr(name, ",");
- for (uint i = 0; i < allnames.size(); i++) {
- std::string singlename = allnames[i];
+ // Setup Maps
+ fEnumIndex[nuisenum]; // = std::vector(0);
+ fNameIndex[name]; // = std::vector(0);
- // Get RW
- genie::rew::GSyst_t rwsyst = GSyst::FromString(singlename);
+ // Split by commas
+ std::vector allnames = GeneralUtils::ParseToStr(name, ",");
+ for (uint i = 0; i < allnames.size(); i++) {
+ std::string singlename = allnames[i];
- // Fill Maps
- int index = fValues.size();
- fValues.push_back(0.0);
- fGENIESysts.push_back(rwsyst);
+ // Get RW
+ genie::rew::GSyst_t rwsyst = GSyst::FromString(singlename);
- // Initialize dial
- std::cout << "Registering " << singlename << " from " << name << std::endl;
- fGenieRW->Systematics().Init( fGENIESysts[index] );
+ // Fill Maps
+ int index = fValues.size();
+ fValues.push_back(0.0);
+ fGENIESysts.push_back(rwsyst);
- // If Absolute
- if (fIsAbsTwk) {
- GSystUncertainty::Instance()->SetUncertainty( rwsyst, 1.0, 1.0 );
- }
+ // Initialize dial
+ std::cout << "Registering " << singlename << " from " << name << std::endl;
+ fGenieRW->Systematics().Init(fGENIESysts[index]);
- // Setup index
- fEnumIndex[nuisenum].push_back(index);
- fNameIndex[name].push_back(index);
+ // If Absolute
+ if (fIsAbsTwk) {
+ GSystUncertainty::Instance()->SetUncertainty(rwsyst, 1.0, 1.0);
+ }
- }
+ // Setup index
+ fEnumIndex[nuisenum].push_back(index);
+ fNameIndex[name].push_back(index);
+ }
- // Set Value if given
- if (startval != -999.9) {
- SetDialValue(nuisenum, startval);
- }
+ // Set Value if given
+ if (startval != -999.9) {
+ SetDialValue(nuisenum, startval);
+ }
#endif
};
-
void GENIEWeightEngine::SetDialValue(int nuisenum, double val) {
#ifdef __GENIE_ENABLED__
- std::vector indices = fEnumIndex[nuisenum];
- for (uint i = 0; i < indices.size(); i++) {
- fValues[indices[i]] = val;
- fGenieRW->Systematics().Set( fGENIESysts[indices[i]], val);
- }
+ std::vector indices = fEnumIndex[nuisenum];
+ for (uint i = 0; i < indices.size(); i++) {
+ fValues[indices[i]] = val;
+ fGenieRW->Systematics().Set(fGENIESysts[indices[i]], val);
+ }
#endif
}
void GENIEWeightEngine::SetDialValue(std::string name, double val) {
#ifdef __GENIE_ENABLED__
- std::vector indices = fNameIndex[name];
- for (uint i = 0; i < indices.size(); i++) {
- fValues[indices[i]] = val;
- fGenieRW->Systematics().Set(fGENIESysts[indices[i]], val);
- }
+ std::vector indices = fNameIndex[name];
+ for (uint i = 0; i < indices.size(); i++) {
+ fValues[indices[i]] = val;
+ fGenieRW->Systematics().Set(fGENIESysts[indices[i]], val);
+ }
#endif
}
void GENIEWeightEngine::Reconfigure(bool silent) {
#ifdef __GENIE_ENABLED__
- // Hush now...
- if (silent) StopTalking();
-
- // Reconf
- fGenieRW->Reconfigure();
- fGenieRW->Print();
-
- // Shout again
- if (silent) StartTalking();
+ // Hush now...
+ if (silent)
+ StopTalking();
+
+ // Reconf
+ fGenieRW->Reconfigure();
+ fGenieRW->Print();
+
+ // Shout again
+ if (silent)
+ StartTalking();
#endif
}
-
-double GENIEWeightEngine::CalcWeight(BaseFitEvt* evt) {
- double rw_weight = 1.0;
+double GENIEWeightEngine::CalcWeight(BaseFitEvt *evt) {
+ double rw_weight = 1.0;
#ifdef __GENIE_ENABLED__
- // Skip Non GENIE
- if (evt->fType != kGENIE) return 1.0;
-
- // Make nom weight
- if (!evt) {
- THROW("evt not found : " << evt);
- }
-
- if (!(evt->genie_event)) {
- THROW("evt->genie_event not found!" << evt->genie_event);
- }
-
- if (!(evt->genie_event->event)) {
- THROW("evt->genie_event->event GHepRecord not found!" << (evt->genie_event->event));
- }
-
- if (!fGenieRW) {
- THROW("GENIE RW Not Found!" << fGenieRW);
- }
-
- rw_weight = fGenieRW->CalcWeight(*(evt->genie_event->event));
- // std::cout << "Returning GENIE Weight for electron scattering = " << rw_weight << std::endl;
+ // Skip Non GENIE
+ if (evt->fType != kGENIE)
+ return 1.0;
+
+ // Make nom weight
+ if (!evt) {
+ THROW("evt not found : " << evt);
+ }
+
+ if (!(evt->genie_event)) {
+ THROW("evt->genie_event not found!" << evt->genie_event);
+ }
+
+ if (!(evt->genie_event->event)) {
+ THROW("evt->genie_event->event GHepRecord not found!"
+ << (evt->genie_event->event));
+ }
+
+ if (!fGenieRW) {
+ THROW("GENIE RW Not Found!" << fGenieRW);
+ }
+
+ rw_weight = fGenieRW->CalcWeight(*(evt->genie_event->event));
+ // std::cout << "Returning GENIE Weight for electron scattering = " <<
+ // rw_weight << std::endl;
#endif
- // Return rw_weight
- return rw_weight;
+ // Return rw_weight
+ return rw_weight;
}
-
-
-
-
-
-
-
-
-
-