diff --git a/cmake/MINERvASetup.cmake b/cmake/MINERvASetup.cmake new file mode 100644 index 0000000..0ccfcd9 --- /dev/null +++ b/cmake/MINERvASetup.cmake @@ -0,0 +1,19 @@ +# 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 . +################################################################################ +LIST(APPEND EXTRA_CXX_FLAGS -D__MINERVA_RW_ENABLED__) diff --git a/cmake/ReweightEnginesSetup.cmake b/cmake/ReweightEnginesSetup.cmake index 44c80ff..87f47d4 100644 --- a/cmake/ReweightEnginesSetup.cmake +++ b/cmake/ReweightEnginesSetup.cmake @@ -1,82 +1,87 @@ # 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 . ################################################################################ ################################## NEUT ###################################### if(USE_NEUT) include(${CMAKE_SOURCE_DIR}/cmake/NEUTSetup.cmake) cmessage(STATUS "Using NEUT Reweight engine.") set(USE_NEUT TRUE CACHE BOOL "Whether to enable NEUT (reweight) support. Requires external libraries. " FORCE) endif() ################################# NuWro ###################################### if(USE_NuWro) include(${CMAKE_SOURCE_DIR}/cmake/NuWroSetup.cmake) cmessage(STATUS "Using NuWro Reweight engine.") set(USE_NuWro TRUE CACHE BOOL "Whether to enable NuWro support. " FORCE) endif() ################################## GENIE ##################################### if(USE_GENIE) include(${CMAKE_SOURCE_DIR}/cmake/GENIESetup.cmake) cmessage(STATUS "Using GENIE Reweight engine.") set(USE_GENIE TRUE CACHE BOOL "Whether to enable GENIE (reweight) support. Requires external libraries. " FORCE) endif() ################################## NIWG ###################################### if(USE_NIWG) include(${CMAKE_SOURCE_DIR}/cmake/NIWGSetup.cmake) cmessage(STATUS "Using NIWG Reweight engine.") set(USE_NIWG TRUE CACHE BOOL "Whether to enable (T2K) NIWG ReWeight support. Requires external libraries. " FORCE) endif() ################################## T2K ###################################### if(USE_T2K) include(${CMAKE_SOURCE_DIR}/cmake/T2KSetup.cmake) cmessage(STATUS "Using T2K Reweight engine.") set(USE_T2K TRUE CACHE BOOL "Whether to enable T2KReWeight support. Requires external libraries. " FORCE) endif() - +################################## MINERvA ###################################### +if(USE_MINERvA_RW) + include(${CMAKE_SOURCE_DIR}/cmake/MINERvASetup.cmake) + cmessage(STATUS "Using MINERvA Reweight engine.") + set(USE_MINERvA_RW TRUE CACHE BOOL "Whether to enable MINERvA ReWeight support. " FORCE) +endif() ################################################################################ ################################ Prob3++ #################################### include(${CMAKE_SOURCE_DIR}/cmake/Prob3++Setup.cmake) ################################################################################ cmessage(STATUS "Reweight engine include directories: ${RWENGINE_INCLUDE_DIRECTORIES}") if(NEED_ROOTEVEGEN) cmessage(STATUS "Require ROOT eve generation libraries") LIST(REVERSE ROOT_LIBS) LIST(APPEND ROOT_LIBS Gui Ged Geom TreePlayer EG Eve) LIST(REVERSE ROOT_LIBS) endif() if(NEED_ROOTPYTHIA6) cmessage(STATUS "Require ROOT Pythia6 libraries") LIST(APPEND ROOT_LIBS EGPythia6 Pythia6) endif() LIST(APPEND EXTRA_LIBS ${ROOT_LIBS}) diff --git a/cmake/cacheVariables.cmake b/cmake/cacheVariables.cmake index 5e88e9c..1c75a96 100644 --- a/cmake/cacheVariables.cmake +++ b/cmake/cacheVariables.cmake @@ -1,206 +1,208 @@ # 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 . ################################################################################ function(CheckAndSetDefaultEnv VARNAME DEFAULT CACHETYPE DOCSTRING ENVNAME) #cmessage(DEBUG "Trying to assign variable ${VARNAME} into the cache.") if(NOT DEFINED ${VARNAME}) if(DEFINED ENV{${ENVNAME}} AND NOT $ENV{${ENVNAME}} STREQUAL "") set(${VARNAME} $ENV{${ENVNAME}} CACHE ${CACHETYPE} ${DOCSTRING}) cmessage(DEBUG " Read ${VARNAME} from ENVVAR ${ENVNAME} as $ENV{${ENVNAME}}.") else() set(${VARNAME} ${DEFAULT} CACHE ${CACHETYPE} ${DOCSTRING}) endif() else() set(${VARNAME} ${${VARNAME}} CACHE ${CACHETYPE} ${DOCSTRING}) unset(${VARNAME}) endif() cmessage(CACHE "--Set cache variable: \"${VARNAME}\" to \"${${VARNAME}}\", in cache ${CACHETYPE}.") endfunction() function(CheckAndSetDefaultCache VARNAME DEFAULT CACHETYPE DOCSTRING) # cmessage(DEBUG "Trying to assign variable ${VARNAME} into the cache.") if(NOT DEFINED ${VARNAME}) set(${VARNAME} ${DEFAULT} CACHE ${CACHETYPE} ${DOCSTRING}) else() set(${VARNAME} ${${VARNAME}} CACHE ${CACHETYPE} ${DOCSTRING}) unset(${VARNAME}) endif() cmessage(CACHE "--Set cache variable: \"${VARNAME}\" to \"${${VARNAME}}\", in cache ${CACHETYPE}.") endfunction() function(CheckAndSetDefault VARNAME DEFAULT) # cmessage(DEBUG "Trying to assign variable ${VARNAME}.") if(NOT DEFINED ${VARNAME}) set(${VARNAME} ${DEFAULT} PARENT_SCOPE) set(${VARNAME} ${DEFAULT}) endif() cmessage(CACHE "--Set variable: \"${VARNAME}\" to \"${${VARNAME}}\".") endfunction() CheckAndSetDefaultCache(VERBOSE TRUE BOOL "Whether to configure loudly.") set (CMAKE_SKIP_BUILD_RPATH TRUE) #Changes default install path to be a subdirectory of the build dir. #Can set build dir at configure time with -DCMAKE_INSTALL_PREFIX=/install/path if(CMAKE_INSTALL_PREFIX STREQUAL "" OR CMAKE_INSTALL_PREFIX STREQUAL "/usr/local") set(CMAKE_INSTALL_PREFIX "${CMAKE_BINARY_DIR}/${CMAKE_SYSTEM_NAME}") elseif(NOT DEFINED CMAKE_INSTALL_PREFIX) set(CMAKE_INSTALL_PREFIX "${CMAKE_BINARY_DIR}/${CMAKE_SYSTEM_NAME}") endif() if(CMAKE_BUILD_TYPE STREQUAL "") set(CMAKE_BUILD_TYPE DEBUG) elseif(NOT DEFINED CMAKE_BUILD_TYPE) set(CMAKE_BUILD_TYPE DEBUG) endif() CheckAndSetDefaultCache(USE_MINIMIZER TRUE INTERNAL "Whether we are using the ROOT minimization libraries. ") CheckAndSetDefaultCache(USE_HEPMC FALSE BOOL "Whether to enable HepMC input support. ") CheckAndSetDefaultEnv(HEPMC "" PATH "Path to HepMC source tree root directory. Overrides environment variable \$HEPMC <>" HEPMC) CheckAndSetDefaultCache(HEPMC_MOMUNIT "GEV" STRING "HepMC momentum units [MEV|GEV]. ") CheckAndSetDefaultCache(HEPMC_LENUNIT "CM" STRING "HepMC momentum units [MM|CM]. ") CheckAndSetDefaultCache(HEPMC_USED_EP FALSE INTERNAL "Whether we built HepMC or not. ") CheckAndSetDefaultCache(USE_NEUT FALSE BOOL "Whether to enable NEUT (reweight) support. Requires external libraries. ") CheckAndSetDefaultEnv(NEUT_ROOT "" PATH "Path to NEUT source tree root directory. Overrides environment variable \$NEUT_ROOT <>" NEUT_ROOT) CheckAndSetDefaultEnv(CERN "" PATH "Path to CERNLIB source tree root directory that NEUT was built against. Overrides environment variable \$CERN <>" CERN) CheckAndSetDefaultEnv(CERN_LEVEL "" STRING "CERNLIB Library version. Overrides environment variable \$CERN_LEVEL <>" CERN_LEVEL) CheckAndSetDefaultCache(USE_NuWro FALSE BOOL "Whether to enable NuWro support. ") CheckAndSetDefaultEnv(NUWRO "" PATH "Path to NuWro source tree root directory. Overrides environment variable \$NUWRO <>" NUWRO) CheckAndSetDefaultEnv(NUWRO_INC "" PATH "Path to NuWro installed includes directory, needs to contain \"params_all.h\". Overrides environment variable \$NUWRO_INC <>" NUWRO_INC) CheckAndSetDefaultCache(NUWRO_INPUT_FILE "" FILEPATH "Path to an input NuWro event vector, which can be used to build NuWro i/o libraries. <>") CheckAndSetDefaultCache(NUWRO_BUILT_FROM_FILE FALSE INTERNAL "Whether the NuWro libraries were built by NUISANCE. ") CheckAndSetDefaultCache(USE_NuWro_RW FALSE BOOL "Whether to try and build support for NuWro reweighting. ") CheckAndSetDefaultCache(USE_NuWro_SRW_Event FALSE BOOL "Whether to use cut down NuWro reweight event format. Requires NuWro reweight. ") CheckAndSetDefaultCache(USE_GENIE FALSE BOOL "Whether to enable GENIE (reweight) support. Requires external libraries. ") CheckAndSetDefaultEnv(GENIE "" PATH "Path to GENIE source tree root directory. Overrides environment variable \$GENIE <>" GENIE) CheckAndSetDefaultEnv(LHAPDF_LIB "" PATH "Path to pre-built LHAPDF libraries. Overrides environment variable \$LHAPDF_LIB. <>" LHAPDF_LIB) CheckAndSetDefaultEnv(LHAPDF_INC "" PATH "Path to installed LHAPDF headers. Overrides environment variable \$LHAPDF_INC. <>" LHAPDF_INC) CheckAndSetDefaultEnv(LHAPATH "" PATH "Path to LHA PDF inputs. Overrides environment variable \$LHAPATH. <>" LHAPATH) CheckAndSetDefaultEnv(LIBXML2_LIB "" PATH "Path to pre-built LIBXML2 libraries. Overrides environment variable \$LIBXML2_LIB. <>" LIBXML2_LIB) CheckAndSetDefaultEnv(LIBXML2_INC "" PATH "Path to installed LIBXML2 headers. Overrides environment variable \$LIBXML2_INC. <>" LIBXML2_INC) CheckAndSetDefaultEnv(LOG4CPP_LIB "" PATH "Path to pre-built LOG4CPP libraries. Overrides environment variable \$LOG4CPP_LIB. <>" LOG4CPP_LIB) CheckAndSetDefaultEnv(LOG4CPP_INC "" PATH "Path to installed LOG4CPP headers. Overrides environment variable \$LOG4CPP_INC. <>" LOG4CPP_INC) CheckAndSetDefaultCache(BUILD_GEVGEN FALSE BOOL "Whether to build nuisance_gevgen app.") CheckAndSetDefaultCache(USE_T2K FALSE BOOL "Whether to enable T2KReWeight support. Requires external libraries. ") CheckAndSetDefaultEnv(T2KREWEIGHT "" PATH "Path to installed T2KREWEIGHTReWeight. Overrides environment variable \$T2KREWEIGHT. <>" T2KREWEIGHT) CheckAndSetDefaultCache(USE_NIWG FALSE BOOL "Whether to enable (T2K) NIWG ReWeight support. Requires external libraries. ") CheckAndSetDefaultEnv(NIWG_ROOT "" PATH "Path to installed NIWGReWeight. Overrides environment variable \$NIWG. <>" NIWG) +CheckAndSetDefaultCache(USE_MINERvA_RW FALSE BOOL "Whether to enable MINERvA ReWeight support. ") + CheckAndSetDefaultEnv(PYTHIA6 "" PATH "Path to directory containing libPythia6.so. Overrides environment variable \$PYTHIA6 <>" PYTHIA6) CheckAndSetDefaultEnv(PYTHIA8 "" PATH "Path to directory containing libPythia8.so. Overrides environment variable \$PYTHIA8 <>" PYTHIA8) CheckAndSetDefaultCache(USE_PYTHIA8 FALSE BOOL "Whether to enable PYTHIA8 event support. ") CheckAndSetDefaultCache(USE_GiBUU TRUE BOOL "Whether to enable GiBUU event support. ") CheckAndSetDefaultCache(BUILD_GiBUU FALSE BOOL "Whether to build supporting GiBUU event tools along with a patched version of GiBUU. ") CheckAndSetDefaultCache(USE_NUANCE TRUE BOOL "Whether to enable NUANCE event support. ") CheckAndSetDefaultCache(USE_PROB3PP FALSE BOOL "Whether to download and compile in Prob3++ support. ") CheckAndSetDefaultCache(NO_EXTERNAL_UPDATE TRUE BOOL "Whether to perform the update target for external dependencies. ") CheckAndSetDefaultCache(USE_GPERFTOOLS FALSE BOOL "Whether to compile in google performance tools. ") CheckAndSetDefault(NEED_PYTHIA6 FALSE) CheckAndSetDefault(NEED_PYTHIA8 FALSE) CheckAndSetDefault(NEED_ROOTEVEGEN FALSE) CheckAndSetDefault(NEED_ROOTPYTHIA6 FALSE) CheckAndSetDefaultCache(USE_OMP FALSE BOOL "Whether to enable multicore features (there currently are none...). ") CheckAndSetDefault(NO_EXPERIMENTS FALSE) cmessage(STATUS "NO_EXPERIMENTS: ${NO_EXPERIMENTS}") CheckAndSetDefaultCache(NO_ANL ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build ANL samples. <-DNO_EXPERIMENTS=FALSE>") CheckAndSetDefaultCache(NO_ArgoNeuT ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build ArgoNeuT samples. <-DNO_EXPERIMENTS=FALSE>") CheckAndSetDefaultCache(NO_BEBC ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build BEBC samples. <-DNO_EXPERIMENTS=FALSE>") CheckAndSetDefaultCache(NO_BNL ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build BNL samples. <-DNO_EXPERIMENTS=FALSE>") CheckAndSetDefaultCache(NO_FNAL ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build FNAL samples. <-DNO_EXPERIMENTS=FALSE>") CheckAndSetDefaultCache(NO_GGM ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build GGM samples. <-DNO_EXPERIMENTS=FALSE>") CheckAndSetDefaultCache(NO_K2K ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build K2K samples. <-DNO_EXPERIMENTS=FALSE>") CheckAndSetDefaultCache(NO_MINERvA ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build MINERvA samples. <-DNO_EXPERIMENTS=FALSE>") CheckAndSetDefaultCache(NO_MiniBooNE ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build MiniBooNE samples. <-DNO_EXPERIMENTS=FALSE>") CheckAndSetDefaultCache(NO_T2K ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build T2K samples. <-DNO_EXPERIMENTS=FALSE>") CheckAndSetDefaultCache(NO_SciBooNE ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build SciBooNE samples. <-DNO_EXPERIMENTS=FALSE>") function(SAYVARS) LIST(APPEND VARS USE_HEPMC HEPMC HEPMC_MOMUNIT HEPMC_LENUNIT HEPMC_USED_EP USE_NEUT NEUT_ROOT CERN CERN_LEVEL USE_NuWro NUWRO NUWRO_INC NUWRO_INPUT_FILE NUWRO_BUILT_FROM_FILE USE_GENIE GENIE LHAPDF_LIB LHAPDF_INC LIBXML2_LIB LIBXML2_INC LOG4CPP_LIB GENIE_LOG4CPP_INC BUILD_GEVGEN USE_T2K USE_NIWG USE_GiBUU BUILD_GiBUU USE_NUANCE NO_EXTERNAL_UPDATE USE_GPERFTOOLS NO_ANL NO_ArgoNeuT NO_BEBC NO_BNL NO_FNAL NO_GGM NO_K2K NO_MINERvA NO_MiniBooNE NO_T2K NO_SciBooNE) foreach(v ${VARS}) if(DEFINED ${v}) cmessage(DEBUG "VARIABLE: \"${v}\" = \"${${v}}\"") endif() endforeach(v) endfunction() diff --git a/src/Reweight/CMakeLists.txt b/src/Reweight/CMakeLists.txt index 1f9270c..685133c 100644 --- a/src/Reweight/CMakeLists.txt +++ b/src/Reweight/CMakeLists.txt @@ -1,79 +1,81 @@ # 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 . ################################################################################ set(IMPLFILES GlobalDialList.cxx FitWeight.cxx WeightEngineBase.cxx NEUTWeightEngine.cxx NuWroWeightEngine.cxx GENIEWeightEngine.cxx WeightUtils.cxx SampleNormEngine.cxx LikelihoodWeightEngine.cxx SplineWeightEngine.cxx NUISANCESyst.cxx T2KWeightEngine.cxx NUISANCEWeightEngine.cxx NUISANCEWeightCalcs.cxx NIWGWeightEngine.cxx OscWeightEngine.cxx +MINERvAWeightCalcs.cxx ) set(HEADERFILES GlobalDialList.h FitWeight.h WeightEngineBase.h NEUTWeightEngine.h NuWroWeightEngine.h GENIEWeightEngine.h WeightUtils.h SampleNormEngine.h LikelihoodWeightEngine.h SplineWeightEngine.h NUISANCESyst.h T2KWeightEngine.h NUISANCEWeightEngine.h NUISANCEWeightCalcs.h NIWGWeightEngine.h OscWeightEngine.h +MINERvAWeightCalcs.h ) set(LIBNAME Reweight) if(CMAKE_BUILD_TYPE MATCHES DEBUG) add_library(${LIBNAME} STATIC ${IMPLFILES}) else(CMAKE_BUILD_TYPE MATCHES RELEASE) add_library(${LIBNAME} SHARED ${IMPLFILES}) endif() include_directories(${MINIMUM_INCLUDE_DIRECTORIES}) set_target_properties(${LIBNAME} PROPERTIES VERSION "${ExtFit_VERSION_MAJOR}.${ExtFit_VERSION_MINOR}.${ExtFit_VERSION_REVISION}") #set_target_properties(${LIBNAME} PROPERTIES LINK_FLAGS ${ROOT_LD_FLAGS}) if(DEFINED PROJECTWIDE_EXTRA_DEPENDENCIES) add_dependencies(${LIBNAME} ${PROJECTWIDE_EXTRA_DEPENDENCIES}) endif() install(TARGETS ${LIBNAME} DESTINATION lib) #Can uncomment this to install the headers... but is it really neccessary? install(FILES ${HEADERFILES} DESTINATION include) set(MODULETargets ${MODULETargets} ${LIBNAME} PARENT_SCOPE) diff --git a/src/Reweight/MINERvAWeightCalcs.cxx b/src/Reweight/MINERvAWeightCalcs.cxx index 5fefde9..ce70e83 100644 --- a/src/Reweight/MINERvAWeightCalcs.cxx +++ b/src/Reweight/MINERvAWeightCalcs.cxx @@ -1,345 +1,403 @@ #ifdef __MINERVA_RW_ENABLED__ #ifdef __GENIE_ENABLED__ #include "MINERvAWeightCalcs.h" #include "BaseFitEvt.h" namespace nuisance { namespace reweight { +//******************************************************* MINERvAReWeight_MEC::MINERvAReWeight_MEC() { + //******************************************************* fTwk_NormCCMEC = 0.0; fDef_NormCCMEC = 1.0; fCur_NormCCMEC = fDef_NormCCMEC; } -virtual ~MINERvAReWeight_MEC::MINERvAReWeight_MEC(){}; +MINERvAReWeight_MEC::~MINERvAReWeight_MEC(){}; double MINERvAReWeight_MEC::CalcWeight(BaseFitEvt* evt) { // Check GENIE if (evt->fType != kGENIE) return 1.0; // Extract the GENIE Record GHepRecord* ghep = static_cast(evt->genie_event->event); const Interaction* interaction = ghep->Summary(); const InitialState& init_state = interaction->InitState(); const ProcessInfo& proc_info = interaction->ProcInfo(); const Target& tgt = init_state.Tgt(); // If the event is not MEC this Calc doesn't handle it if (!proc_info.IsMEC()) return 1.0; // WEIGHT CALCULATIONS ------------- double w = 1.0; // CCMEC Dial - if (!proc_info.IsCC()) w *= fCur_NormCCMEC; + if (!proc_info.IsWeakCC()) w *= fCur_NormCCMEC; // Return Combined Weight return w; } void MINERvAReWeight_MEC::SetDialValue(std::string name, double val) { SetDialValue(Reweight::ConvDial(name, kCUSTOM), val); } void MINERvAReWeight_MEC::SetDialValue(int rwenum, double val) { // Check Handled int curenum = rwenum % 1000; if (!IsHandled(curenum)) return; // Set Values - if (curenum == kMINERvARW_NormCCMEC) fTwk_NormCCMEC = val; + if (curenum == kMINERvARW_NormCCMEC) { + fTwk_NormCCMEC = val; + fCur_NormCCMEC = fDef_NormCCMEC + fTwk_NormCCMEC; + } // Define Tweaked fTweaked = ((fTwk_NormCCMEC != 0.0)); } bool MINERvAReWeight_MEC::IsHandled(int rwenum) { int curenum = rwenum % 1000; switch (curenum) { case kMINERvARW_NormCCMEC: return true; default: return false; } } + +//******************************************************* +MINERvAReWeight_RES::MINERvAReWeight_RES() { + //******************************************************* + fTwk_NormCCRES = 0.0; + fDef_NormCCRES = 1.0; + fCur_NormCCRES = fDef_NormCCRES; +} + +MINERvAReWeight_RES::~MINERvAReWeight_RES(){}; + +double MINERvAReWeight_RES::CalcWeight(BaseFitEvt* evt) { + + // std::cout << "Caculating RES" << std::endl; + // Check GENIE + if (evt->fType != kGENIE) return 1.0; + + // Extract the GENIE Record + GHepRecord* ghep = static_cast(evt->genie_event->event); + const Interaction* interaction = ghep->Summary(); + const InitialState& init_state = interaction->InitState(); + const ProcessInfo& proc_info = interaction->ProcInfo(); + const Target& tgt = init_state.Tgt(); + + // If the event is not RES this Calc doesn't handle it + if (!proc_info.IsResonant()) return 1.0; + + // WEIGHT CALCULATIONS ------------- + double w = 1.0; + + // CCRES Dial + if (proc_info.IsWeakCC()) w *= fCur_NormCCRES; + + // Return Combined Weight + return w; +} + +void MINERvAReWeight_RES::SetDialValue(std::string name, double val) { + SetDialValue(Reweight::ConvDial(name, kCUSTOM), val); +} + +void MINERvAReWeight_RES::SetDialValue(int rwenum, double val) { + // Check Handled + int curenum = rwenum % 1000; + if (!IsHandled(curenum)) return; + + // Set Values + if (curenum == kMINERvARW_NormCCRES) { + fTwk_NormCCRES = val; + fCur_NormCCRES = fDef_NormCCRES + fTwk_NormCCRES; + } + + // Define Tweaked + fTweaked = ((fTwk_NormCCRES != 0.0)); +} + +bool MINERvAReWeight_RES::IsHandled(int rwenum) { + int curenum = rwenum % 1000; + + switch (curenum) { + case kMINERvARW_NormCCRES: + return true; + default: + return false; + } +} + +//******************************************************* +RikRPA::RikRPA() { + //******************************************************* + + // - Syst : kMINERvA_RikRPA_ApplyRPA + // - Type : Binary + // - Limits : 0.0 (false) -> 1.0 (true) + // - Default : 0.0 + fApplyDial_RPACorrection = false; + + // - Syst : kMINERvA_RikRPA_LowQ2 + // - Type : Absolute + // - Limits : 1.0 -> 1.0 + // - Default : 0.0 + // - Frac Error : 100% + fDefDial_RPALowQ2 = 0.0; + fCurDial_RPALowQ2 = fDefDial_RPALowQ2; + fErrDial_RPALowQ2 = 0.0; + + // - Syst : kMINERvA_RikRPA_HighQ2 + // - Type : Absolute + // - Limits : 1.0 -> 1.0 + // - Default : 0.0 + // - Frac Error : 100% + fDefDial_RPAHighQ2 = 0.0; + fCurDial_RPAHighQ2 = fDefDial_RPAHighQ2; + fErrDial_RPAHighQ2 = 1.0; + + fEventWeights = new double[5]; + + for (size_t i = 0; i < kMaxCalculators; i++) { + // fRPACalculators[i] = NULL; + } + + fTweaked = false; +} + +RikRPA::~RikRPA() { + delete fEventWeights; + + for (size_t i = 0; i < kMaxCalculators; i++) { + // if (fRPACalculators[i]) delete fRPACalculators[i]; + // fRPACalculators[i] = NULL; + } +} + +double RikRPA::CalcWeight(BaseFitEvt* evt) { + // LOG(FIT) << "Calculating RikRPA" << std::endl; + // Return 1.0 if not tweaked + if (!fTweaked) return 1.0; + + double w = 1.0; + + // Extract the GENIE Record + GHepRecord* ghep = static_cast(evt->genie_event->event); + const Interaction* interaction = ghep->Summary(); + const InitialState& init_state = interaction->InitState(); + const ProcessInfo& proc_info = interaction->ProcInfo(); + // const Kinematics & kine = interaction->Kine(); + // const XclsTag & xcls = interaction->ExclTag(); + const Target& tgt = init_state.Tgt(); + + // If not QE return 1.0 + // LOG(FIT) << "RikRPA : Event QE = " << proc_info.IsQuasiElastic() << + // std::endl; + if (!tgt.IsNucleus()) return 1.0; + if (!proc_info.IsQuasiElastic()) return 1.0; + + // Extract Beam and Target PDG + GHepParticle* neutrino = ghep->Probe(); + int bpdg = neutrino->Pdg(); + + GHepParticle* target = ghep->Particle(1); + assert(target); + int tpdg = target->Pdg(); + + // Find the enum we need + int calcenum = GetRPACalcEnum(bpdg, tpdg); + if (calcenum == -1) return 1.0; + + // Check we have the RPA Calc setup for this enum + // if not, set it up at that point + // if (!fRPACalculators[calcenum]) SetupRPACalculator(calcenum); + // weightRPA* rpacalc = fRPACalculators[calcenum]; + // if (!rpacalc) { + // THROW("Failed to grab the RPA Calculator : " << calcenum); + // } + + // Extract Q0-Q3 + GHepParticle* fsl = ghep->FinalStatePrimaryLepton(); + const TLorentzVector& k1 = *(neutrino->P4()); + const TLorentzVector& k2 = *(fsl->P4()); + double q0 = fabs((k1 - k2).E()); + double q3 = fabs((k1 - k2).Vect().Mag()); + + // Now use q0-q3 and RPA Calculator to fill fWeights + // LOG(FIT) << "Getting Weights = " << q0 << " " << q3 << std::endl; + // rpacalc->getWeight(q0, q3, fEventWeights); + + // Apply Interpolation (for the time being simple linear) + + // Syst Application : kMINERvA_RikRPA_ApplyRPA + if (fApplyDial_RPACorrection) { + w *= fEventWeights[0]; // CV + } + + LOG(FIT) << " fCurDial_RPALowQ2 = " << fCurDial_RPALowQ2 + << " fCurDial_RPAHighQ2 = " << fCurDial_RPAHighQ2 << " Weights " + << fEventWeights[0] << " " << fEventWeights[1] << " " + << fEventWeights[2] << " " << fEventWeights[3] << " " + << fEventWeights[4] << std::endl; + + // Syst Application : kMINERvA_RikRPA_LowQ2 + if (fabs(fCurDial_RPALowQ2) > 0.0) { + double interpw = fEventWeights[0]; + + if (fCurDial_RPALowQ2 > 0.0) { + interpw = fEventWeights[0] - (fEventWeights[0] - fEventWeights[1]) * + fCurDial_RPALowQ2; // WLow+ } else if + } else if (fCurDial_RPALowQ2 < 0.0) { + interpw = fEventWeights[0] - (fEventWeights[2] - fEventWeights[0]) * + fCurDial_RPALowQ2; // WLow- + } + w *= interpw / fEventWeights[0]; // Div by CV again + } + + // Syst Application : kMINERvA_RikRPA_HighQ2 + if (fabs(fCurDial_RPAHighQ2) > 0.0) { + double interpw = fEventWeights[0]; + if (fCurDial_RPAHighQ2 > 0.0) { + interpw = fEventWeights[0] - (fEventWeights[0] - fEventWeights[3]) * + fCurDial_RPAHighQ2; // WHigh+ } else + if (fCurDial_RPAHighQ2 < 0.0) { + interpw = fEventWeights[0] - (fEventWeights[4] - fEventWeights[0]) * + fCurDial_RPAHighQ2; // WHigh- + } + w *= interpw / fEventWeights[0]; // Div by CV again + } + } + + // LOG(FIT) << "RPA Weight = " << w << std::endl; + return w; } // namespace reweight -} // namespace nuisance -#endif -#endif +void RikRPA::SetDialValue(std::string name, double val) { + SetDialValue(Reweight::ConvDial(name, kCUSTOM), val); +} + +void RikRPA::SetDialValue(int rwenum, double val) { + int curenum = rwenum % 1000; + // Check Handled + if (!IsHandled(curenum)) return; + if (curenum == kMINERvARW_RikRPA_ApplyRPA) + fApplyDial_RPACorrection = (val > 0.5); + if (curenum == kMINERvARW_RikRPA_LowQ2) fCurDial_RPALowQ2 = val; + if (curenum == kMINERvARW_RikRPA_HighQ2) fCurDial_RPAHighQ2 = val; + + // Assign flag to say stuff has changed + fTweaked = (fApplyDial_RPACorrection || + fabs(fCurDial_RPAHighQ2 - fDefDial_RPAHighQ2) > 0.0 || + fabs(fCurDial_RPALowQ2 - fDefDial_RPALowQ2) > 0.0); +} +bool RikRPA::IsHandled(int rwenum) { + int curenum = rwenum % 1000; + switch (curenum) { + case kMINERvARW_RikRPA_ApplyRPA: + return true; + case kMINERvARW_RikRPA_LowQ2: + return true; + case kMINERvARW_RikRPA_HighQ2: + return true; + default: + return false; + } +} +void RikRPA::SetupRPACalculator(int calcenum) { + std::string rwdir = FitPar::GetDataBase() + "/reweight/MINERvA/RikRPA/"; + std::string fidir = ""; + switch (calcenum) { + case kNuMuC12: + fidir = "outNievesRPAratio-nu12C-20GeV-20170202.root"; + break; + + case kNuMuO16: + fidir = "outNievesRPAratio-nu16O-20GeV-20170202.root"; + break; + + case kNuMuAr40: + fidir = "outNievesRPAratio-nu40Ar-20GeV-20170202.root"; + break; + + case kNuMuCa40: + fidir = "outNievesRPAratio-nu40Ca-20GeV-20170202.root"; + break; + + case kNuMuFe56: + fidir = "outNievesRPAratio-nu56Fe-20GeV-20170202.root"; + break; + + case kNuMuBarC12: + fidir = "outNievesRPAratio-anu12C-20GeV-20170202.root"; + break; + + case kNuMuBarO16: + fidir = "outNievesRPAratio-anu16O-20GeV-20170202.root"; + break; + + case kNuMuBarAr40: + fidir = "outNievesRPAratio-anu40Ar-20GeV-20170202.root"; + break; + + case kNuMuBarCa40: + fidir = "outNievesRPAratio-anu40Ca-20GeV-20170202.root"; + break; + + case kNuMuBarFe56: + fidir = "outNievesRPAratio-anu56Fe-20GeV-20170202.root"; + break; + } + LOG(FIT) << "Loading RPA CALC : " << fidir << std::endl; + TDirectory* olddir = gDirectory; -// RikRPA::RikRPA() { + // fRPACalculators[calcenum] = new weightRPA(rwdir + "/" + fidir); + olddir->cd(); + return; +} -// // - Syst : kMINERvA_RikRPA_ApplyRPA -// // - Type : Binary -// // - Limits : 0.0 (false) -> 1.0 (true) -// // - Default : 0.0 -// fApplyDial_RPACorrection = false; +int RikRPA::GetRPACalcEnum(int bpdg, int tpdg) { + if (bpdg == 14 && tpdg == 1000060120) + return kNuMuC12; + else if (bpdg == 14 && tpdg == 1000080160) + return kNuMuO16; + else if (bpdg == 14 && tpdg == 1000180400) + return kNuMuAr40; + else if (bpdg == 14 && tpdg == 1000200400) + return kNuMuCa40; + else if (bpdg == 14 && tpdg == 1000280560) + return kNuMuFe56; + else if (bpdg == -14 && tpdg == 1000060120) + return kNuMuBarC12; + else if (bpdg == -14 && tpdg == 1000080160) + return kNuMuBarO16; + else if (bpdg == -14 && tpdg == 1000180400) + return kNuMuBarAr40; + else if (bpdg == -14 && tpdg == 1000200400) + return kNuMuBarCa40; + else if (bpdg == -14 && tpdg == 1000280560) + return kNuMuBarFe56; + else { + ERROR(WRN, "Unknown beam and target combination for RPA Calcs! " + << bpdg << " " << tpdg); + } -// // - Syst : kMINERvA_RikRPA_LowQ2 -// // - Type : Absolute -// // - Limits : 1.0 -> 1.0 -// // - Default : 0.0 -// // - Frac Error : 100% -// fDefDial_RPALowQ2 = 0.0; -// fCurDial_RPALowQ2 = fDefDial_RPALowQ2; -// fErrDial_RPALowQ2 = 0.0; + return -1; +} -// // - Syst : kMINERvA_RikRPA_HighQ2 -// // - Type : Absolute -// // - Limits : 1.0 -> 1.0 -// // - Default : 0.0 -// // - Frac Error : 100% -// fDefDial_RPAHighQ2 = 0.0; -// fCurDial_RPAHighQ2 = fDefDial_RPAHighQ2; -// fErrDial_RPAHighQ2 = 1.0; - -// fEventWeights = new double[5]; - -// for (size_t i = 0; i < kMaxCalculators; i++) { -// fRPACalculators[i] = NULL; -// } - -// fTweaked = false; -// } - -// RikRPA::~RikRPA() { - -// delete fEventWeights; - -// for (size_t i = 0; i < kMaxCalculators; i++) { -// if (fRPACalculators[i]) delete fRPACalculators[i]; -// fRPACalculators[i] = NULL; -// } - -// } - -// double RikRPA::CalcWeight(BaseFitEvt* evt) { -// // LOG(FIT) << "Calculating RikRPA" << std::endl; -// // Return 1.0 if not tweaked -// if (!fTweaked) return 1.0; - -// double w = 1.0; - -// // Extract the GENIE Record -// GHepRecord* ghep = static_cast(evt->genie_event->event); -// const Interaction* interaction = ghep->Summary(); -// const InitialState& init_state = interaction->InitState(); -// const ProcessInfo& proc_info = interaction->ProcInfo(); -// // const Kinematics & kine = interaction->Kine(); -// // const XclsTag & xcls = interaction->ExclTag(); -// const Target& tgt = init_state.Tgt(); - -// // If not QE return 1.0 -// // LOG(FIT) << "RikRPA : Event QE = " << proc_info.IsQuasiElastic() << -// std::endl; -// if (!tgt.IsNucleus()) return 1.0; -// if (!proc_info.IsQuasiElastic()) return 1.0; - -// // Extract Beam and Target PDG -// GHepParticle* neutrino = ghep->Probe(); -// int bpdg = neutrino->Pdg(); - -// GHepParticle* target = ghep->Particle(1); -// assert(target); -// int tpdg = target->Pdg(); - -// // Find the enum we need -// int calcenum = GetRPACalcEnum(bpdg, tpdg); -// if (calcenum == -1) return 1.0; - -// // Check we have the RPA Calc setup for this enum -// // if not, set it up at that point -// if (!fRPACalculators[calcenum]) SetupRPACalculator(calcenum); -// weightRPA* rpacalc = fRPACalculators[calcenum]; -// if (!rpacalc) { -// THROW("Failed to grab the RPA Calculator : " << calcenum); -// } - -// // Extract Q0-Q3 -// GHepParticle* fsl = ghep->FinalStatePrimaryLepton(); -// const TLorentzVector& k1 = *(neutrino->P4()); -// const TLorentzVector& k2 = *(fsl->P4()); -// double q0 = fabs((k1 - k2).E()); -// double q3 = fabs((k1 - k2).Vect().Mag()); - -// // Now use q0-q3 and RPA Calculator to fill fWeights -// // LOG(FIT) << "Getting Weights = " << q0 << " " << q3 << std::endl; -// rpacalc->getWeight(q0, q3, fEventWeights); - -// // Apply Interpolation (for the time being simple linear) - -// // Syst Application : kMINERvA_RikRPA_ApplyRPA -// if (fApplyDial_RPACorrection) { -// w *= fEventWeights[0]; // CV -// } - -// LOG(FIT) << " fCurDial_RPALowQ2 = " << fCurDial_RPALowQ2 -// << " fCurDial_RPAHighQ2 = " << fCurDial_RPAHighQ2 << " Weights " -// << fEventWeights[0] << " " << fEventWeights[1] << " " -// << fEventWeights[2] << " " << fEventWeights[3] << " " -// << fEventWeights[4] << std::endl; - -// // Syst Application : kMINERvA_RikRPA_LowQ2 -// if (fabs(fCurDial_RPALowQ2) > 0.0) { -// double interpw = fEventWeights[0]; - -// if (fCurDial_RPALowQ2 > 0.0) { -// interpw = fEventWeights[0] - -// (fEventWeights[0] - fEventWeights[1]) * -// fCurDial_RPALowQ2; // WLow+ } else if -// (fCurDial_RPALowQ2 < 0.0) { -// interpw = fEventWeights[0] - (fEventWeights[2] - fEventWeights[0]) * -// fCurDial_RPALowQ2; // WLow- -// } -// w *= interpw / fEventWeights[0]; // Div by CV again -// } - -// // Syst Application : kMINERvA_RikRPA_HighQ2 -// if (fabs(fCurDial_RPAHighQ2) > 0.0) { -// double interpw = fEventWeights[0]; -// if (fCurDial_RPAHighQ2 > 0.0) { -// interpw = fEventWeights[0] - -// (fEventWeights[0] - fEventWeights[3]) * -// fCurDial_RPAHighQ2; // WHigh+ } else -// if -// (fCurDial_RPAHighQ2 < 0.0) { -// interpw = fEventWeights[0] - (fEventWeights[4] - fEventWeights[0]) -// * -// fCurDial_RPAHighQ2; // WHigh- -// } -// w *= interpw / fEventWeights[0]; // Div by CV again -// } - -// // LOG(FIT) << "RPA Weight = " << w << std::endl; - -// return w; -// } - -// void RikRPA::SetDialValue(std::string name, double val) { -// SetDialValue(Reweight::ConvDial(name, kCUSTOM), val); -// } - -// void RikRPA::SetDialValue(int rwenum, double val) { -// int curenum = rwenum % 1000; - -// // Check Handled -// if (!IsHandled(curenum)) return; -// if (curenum == kMINERvA_RikRPA_ApplyRPA) -// fApplyDial_RPACorrection = (val > 0.5); -// if (curenum == kMINERvA_RikRPA_LowQ2) fCurDial_RPALowQ2 = val; -// if (curenum == kMINERvA_RikRPA_HighQ2) fCurDial_RPAHighQ2 = val; - -// // Assign flag to say stuff has changed -// fTweaked = (fApplyDial_RPACorrection || -// fabs(fCurDial_RPAHighQ2 - fDefDial_RPAHighQ2) > 0.0 || -// fabs(fCurDial_RPALowQ2 - fDefDial_RPALowQ2) > 0.0); -// } - -// bool RikRPA::IsHandled(int rwenum) { -// int curenum = rwenum % 1000; -// switch (curenum) { -// case kMINERvA_RikRPA_ApplyRPA: -// return true; -// case kMINERvA_RikRPA_LowQ2: -// return true; -// case kMINERvA_RikRPA_HighQ2: -// return true; -// default: -// return false; -// } -// } - -// void RikRPA::SetupRPACalculator(int calcenum) { -// std::string rwdir = FitPar::GetDataBase() + -// "/reweight/MINERvA/RikRPA/"; std::string fidir = ""; switch (calcenum) -// { -// case kNuMuC12: -// fidir = "outNievesRPAratio-nu12C-20GeV-20170202.root"; -// break; - -// case kNuMuO16: -// fidir = "outNievesRPAratio-nu16O-20GeV-20170202.root"; -// break; - -// case kNuMuAr40: -// fidir = "outNievesRPAratio-nu40Ar-20GeV-20170202.root"; -// break; - -// case kNuMuCa40: -// fidir = "outNievesRPAratio-nu40Ca-20GeV-20170202.root"; -// break; - -// case kNuMuFe56: -// fidir = "outNievesRPAratio-nu56Fe-20GeV-20170202.root"; -// break; - -// case kNuMuBarC12: -// fidir = "outNievesRPAratio-anu12C-20GeV-20170202.root"; -// break; - -// case kNuMuBarO16: -// fidir = "outNievesRPAratio-anu16O-20GeV-20170202.root"; -// break; - -// case kNuMuBarAr40: -// fidir = "outNievesRPAratio-anu40Ar-20GeV-20170202.root"; -// break; - -// case kNuMuBarCa40: -// fidir = "outNievesRPAratio-anu40Ca-20GeV-20170202.root"; -// break; - -// case kNuMuBarFe56: -// fidir = "outNievesRPAratio-anu56Fe-20GeV-20170202.root"; -// break; -// } - -// LOG(FIT) << "Loading RPA CALC : " << fidir << std::endl; -// TDirectory* olddir = gDirectory; - -// fRPACalculators[calcenum] = new weightRPA(rwdir + "/" + fidir); -// olddir->cd(); -// return; -// } - -// int RikRPA::GetRPACalcEnum(int bpdg, int tpdg) { -// if (bpdg == 14 && tpdg == 1000060120) -// return kNuMuC12; -// else if (bpdg == 14 && tpdg == 1000080160) -// return kNuMuO16; -// else if (bpdg == 14 && tpdg == 1000180400) -// return kNuMuAr40; -// else if (bpdg == 14 && tpdg == 1000200400) -// return kNuMuCa40; -// else if (bpdg == 14 && tpdg == 1000280560) -// return kNuMuFe56; -// else if (bpdg == -14 && tpdg == 1000060120) -// return kNuMuBarC12; -// else if (bpdg == -14 && tpdg == 1000080160) -// return kNuMuBarO16; -// else if (bpdg == -14 && tpdg == 1000180400) -// return kNuMuBarAr40; -// else if (bpdg == -14 && tpdg == 1000200400) -// return kNuMuBarCa40; -// else if (bpdg == -14 && tpdg == 1000280560) -// return kNuMuBarFe56; -// else { -// ERROR(WRN, "Unknown beam and target combination for RPA Calcs! " -// << bpdg << " " << tpdg); -// } - -// return -1; -// } -// } // namespace reweight -// } // namespace nuisance - -// #endif // __GENIE_ENABLED__ -// #endif //__MINERVA_RW_ENABLED__ +} // namespace reweight +} // namespace nuisance +#endif +#endif diff --git a/src/Reweight/MINERvAWeightCalcs.h b/src/Reweight/MINERvAWeightCalcs.h index ec059dc..8bcba23 100644 --- a/src/Reweight/MINERvAWeightCalcs.h +++ b/src/Reweight/MINERvAWeightCalcs.h @@ -1,96 +1,117 @@ #ifndef MINERVA_WEIGHT_CALCS #define MINERVA_WEIGHT_CALCS #include -// #ifdef __MINERVA_RW_ENABLED__ -// #ifdef __GENIE_ENABLED__ +#ifdef __MINERVA_RW_ENABLED__ +#ifdef __GENIE_ENABLED__ #include "Conventions/Units.h" #include "EVGCore/EventRecord.h" #include "FitEvent.h" #include "FitParameters.h" #include "GHEP/GHepParticle.h" #include "GHEP/GHepRecord.h" #include "GHEP/GHepUtils.h" #include "GeneralUtils.h" #include "NUISANCESyst.h" #include "NUISANCEWeightCalcs.h" #include "Ntuple/NtpMCEventRecord.h" #include "PDG/PDGUtils.h" #include "WeightUtils.h" // #include "weightRPA.h" using namespace genie; class BaseFitEvt; namespace nuisance { namespace reweight { // MEC Dials class MINERvAReWeight_MEC : public NUISANCEWeightCalc { public: MINERvAReWeight_MEC(); virtual ~MINERvAReWeight_MEC(); double CalcWeight(BaseFitEvt* evt); void SetDialValue(std::string name, double val); void SetDialValue(int rwenum, double val); bool IsHandled(int rwenum); double fTwk_NormCCMEC; double fCur_NormCCMEC; double fDef_NormCCMEC; -} - -// /// RPA Weight Calculator that applies RPA systematics -// /// to GENIE events. GENIE EVENTS ONLY! -// class RikRPA : public NUISANCEWeightCalc { -// public: -// RikRPA(); -// ~RikRPA(); - -// double CalcWeight(BaseFitEvt* evt); -// void SetDialValue(std::string name, double val); -// void SetDialValue(int rwenum, double val); -// bool IsHandled(int rwenum); - -// void SetupRPACalculator(int calcenum); -// int GetRPACalcEnum(int bpdg, int tpdg); - -// bool fApplyDial_RPACorrection; - -// double fTwkDial_RPALowQ2; -// double fDefDial_RPALowQ2; -// double fCurDial_RPALowQ2; -// double fErrDial_RPALowQ2; - -// double fTwkDial_RPAHighQ2; -// double fDefDial_RPAHighQ2; -// double fCurDial_RPAHighQ2; -// double fErrDial_RPAHighQ2; - -// double* fEventWeights; -// bool fTweaked; - -// const static int kMaxCalculators = 10; -// enum rpacalcenums { -// kNuMuC12, -// kNuMuO16, -// kNuMuAr40, -// kNuMuCa40, -// kNuMuFe56, -// kNuMuBarC12, -// kNuMuBarO16, -// kNuMuBarAr40, -// kNuMuBarCa40, -// kNuMuBarFe56 -// }; -// weightRPA* fRPACalculators[kMaxCalculators]; -// }; + bool fTweaked; + +}; + +// RES Dials +class MINERvAReWeight_RES : public NUISANCEWeightCalc { + public: + MINERvAReWeight_RES(); + virtual ~MINERvAReWeight_RES(); + + double CalcWeight(BaseFitEvt* evt); + void SetDialValue(std::string name, double val); + void SetDialValue(int rwenum, double val); + bool IsHandled(int rwenum); + + double fTwk_NormCCRES; + double fCur_NormCCRES; + double fDef_NormCCRES; + bool fTweaked; + + }; + +/// RPA Weight Calculator that applies RPA systematics +/// to GENIE events. GENIE EVENTS ONLY! +class RikRPA : public NUISANCEWeightCalc { +public: + RikRPA(); + ~RikRPA(); + + double CalcWeight(BaseFitEvt* evt); + void SetDialValue(std::string name, double val); + void SetDialValue(int rwenum, double val); + bool IsHandled(int rwenum); + + void SetupRPACalculator(int calcenum); + int GetRPACalcEnum(int bpdg, int tpdg); + + bool fApplyDial_RPACorrection; + + double fTwkDial_RPALowQ2; + double fDefDial_RPALowQ2; + double fCurDial_RPALowQ2; + double fErrDial_RPALowQ2; + + double fTwkDial_RPAHighQ2; + double fDefDial_RPAHighQ2; + double fCurDial_RPAHighQ2; + double fErrDial_RPAHighQ2; + + double* fEventWeights; + bool fTweaked; + + const static int kMaxCalculators = 10; + enum rpacalcenums { + kNuMuC12, + kNuMuO16, + kNuMuAr40, + kNuMuCa40, + kNuMuFe56, + kNuMuBarC12, + kNuMuBarO16, + kNuMuBarAr40, + kNuMuBarCa40, + kNuMuBarFe56 + }; + // weightRPA* fRPACalculators[kMaxCalculators]; +}; + }; // namespace reweight }; // namespace nuisance -// #endif // __GENIE_ENABLED__ -// #endif //__MINERVA_RW_ENABLED__ +#endif // __GENIE_ENABLED__ +#endif //__MINERVA_RW_ENABLED__ #endif diff --git a/src/Reweight/NUISANCESyst.cxx b/src/Reweight/NUISANCESyst.cxx index 44366e1..40aef5e 100644 --- a/src/Reweight/NUISANCESyst.cxx +++ b/src/Reweight/NUISANCESyst.cxx @@ -1,52 +1,60 @@ #include "NUISANCESyst.h" int Reweight::ConvertNUISANCEDial(std::string type) { for (int i = kUnkownNUISANCEDial + 1; i < kNUISANCEDial_LAST; i++) { if (!type.compare(ConvNUISANCEDial(i).c_str())) { return i; } } return kUnkownNUISANCEDial; }; std::string Reweight::ConvNUISANCEDial(int type) { switch (type) { case kGaussianCorr_CCQE_norm: { return "GaussianCorr_CCQE_norm"; } case kGaussianCorr_CCQE_tilt: { return "GaussianCorr_CCQE_tilt"; } case kGaussianCorr_CCQE_Pq0: { return "GaussianCorr_CCQE_Pq0"; } case kGaussianCorr_CCQE_Wq0: { return "GaussianCorr_CCQE_Wq0"; } case kGaussianCorr_CCQE_Pq3: { return "GaussianCorr_CCQE_Pq3"; } case kGaussianCorr_CCQE_Wq3: { return "GaussianCorr_CCQE_Wq3"; } case kGaussianCorr_2p2h_norm: { return "GaussianCorr_2p2h_norm"; } case kGaussianCorr_2p2h_tilt: { return "GaussianCorr_2p2h_tilt"; } case kGaussianCorr_2p2h_Pq0: { return "GaussianCorr_2p2h_Pq0"; } case kGaussianCorr_2p2h_Wq0: { return "GaussianCorr_2p2h_Wq0"; } case kGaussianCorr_2p2h_Pq3: { return "GaussianCorr_2p2h_Pq3"; } case kGaussianCorr_2p2h_Wq3: { return "GaussianCorr_2p2h_Wq3"; } case kGaussianCorr_2p2h_PPandNN_norm: { return "GaussianCorr_2p2h_PPandNN_norm"; } case kGaussianCorr_2p2h_PPandNN_tilt: { return "GaussianCorr_2p2h_PPandNN_tilt"; } case kGaussianCorr_2p2h_PPandNN_Pq0: { return "GaussianCorr_2p2h_PPandNN_Pq0"; } case kGaussianCorr_2p2h_PPandNN_Wq0: { return "GaussianCorr_2p2h_PPandNN_Wq0"; } case kGaussianCorr_2p2h_PPandNN_Pq3: { return "GaussianCorr_2p2h_PPandNN_Pq3"; } case kGaussianCorr_2p2h_PPandNN_Wq3: { return "GaussianCorr_2p2h_PPandNN_Wq3"; } case kGaussianCorr_2p2h_NP_norm: { return "GaussianCorr_2p2h_NP_norm"; } case kGaussianCorr_2p2h_NP_tilt: { return "GaussianCorr_2p2h_NP_tilt"; } case kGaussianCorr_2p2h_NP_Pq0: { return "GaussianCorr_2p2h_NP_Pq0"; } case kGaussianCorr_2p2h_NP_Wq0: { return "GaussianCorr_2p2h_NP_Wq0"; } case kGaussianCorr_2p2h_NP_Pq3: { return "GaussianCorr_2p2h_NP_Pq3"; } case kGaussianCorr_2p2h_NP_Wq3: { return "GaussianCorr_2p2h_NP_Wq3"; } case kGaussianCorr_CC1pi_norm: { return "GaussianCorr_CC1pi_norm"; } case kGaussianCorr_CC1pi_tilt: { return "GaussianCorr_CC1pi_tilt"; } case kGaussianCorr_CC1pi_Pq0: { return "GaussianCorr_CC1pi_Pq0"; } case kGaussianCorr_CC1pi_Wq0: { return "GaussianCorr_CC1pi_Wq0"; } case kGaussianCorr_CC1pi_Pq3: { return "GaussianCorr_CC1pi_Pq3"; } case kGaussianCorr_CC1pi_Wq3: { return "GaussianCorr_CC1pi_Wq3"; } case kGaussianCorr_AllowSuppression: { return "GaussianCorr_AllowSuppression"; } + case kModeNorm_NormRES: { return "NormRES"; } + case kMINERvARW_NormCCMEC: { return "MINERvARW_NormCCMEC"; } + case kMINERvARW_NormCCRES: { return "MINERvARW_NormCCRES"; } + + case kMINERvARW_RikRPA_ApplyRPA: { return "MINERvARW_RikRPA_ApplyRPA"; } + case kMINERvARW_RikRPA_LowQ2: { return "MINERvARW_RikRPA_LowQ2"; } + case kMINERvARW_RikRPA_HighQ2: { return "MINERvARW_RikRPA_HighQ2"; } + default: return "unknown_nuisance_dial"; } }; diff --git a/src/Reweight/NUISANCESyst.h b/src/Reweight/NUISANCESyst.h index 296c956..a3cfb74 100644 --- a/src/Reweight/NUISANCESyst.h +++ b/src/Reweight/NUISANCESyst.h @@ -1,50 +1,57 @@ #ifndef NUISANCESyst_H #define NUISANCESyst_H #include "GeneralUtils.h" namespace Reweight { enum NUISANCESyst { kUnkownNUISANCEDial = 0, kGaussianCorr_CCQE_norm, kGaussianCorr_CCQE_tilt, kGaussianCorr_CCQE_Pq0, kGaussianCorr_CCQE_Wq0, kGaussianCorr_CCQE_Pq3, kGaussianCorr_CCQE_Wq3, kGaussianCorr_2p2h_norm, kGaussianCorr_2p2h_tilt, kGaussianCorr_2p2h_Pq0, kGaussianCorr_2p2h_Wq0, kGaussianCorr_2p2h_Pq3, kGaussianCorr_2p2h_Wq3, kGaussianCorr_2p2h_PPandNN_norm, kGaussianCorr_2p2h_PPandNN_tilt, kGaussianCorr_2p2h_PPandNN_Pq0, kGaussianCorr_2p2h_PPandNN_Wq0, kGaussianCorr_2p2h_PPandNN_Pq3, kGaussianCorr_2p2h_PPandNN_Wq3, kGaussianCorr_2p2h_NP_norm, kGaussianCorr_2p2h_NP_tilt, kGaussianCorr_2p2h_NP_Pq0, kGaussianCorr_2p2h_NP_Wq0, kGaussianCorr_2p2h_NP_Pq3, kGaussianCorr_2p2h_NP_Wq3, kGaussianCorr_CC1pi_norm, kGaussianCorr_CC1pi_tilt, kGaussianCorr_CC1pi_Pq0, kGaussianCorr_CC1pi_Wq0, kGaussianCorr_CC1pi_Pq3, kGaussianCorr_CC1pi_Wq3, kGaussianCorr_AllowSuppression, kModeNorm_NormRES, + kMINERvARW_NormCCMEC, + kMINERvARW_NormCCRES, + + kMINERvARW_RikRPA_ApplyRPA, + kMINERvARW_RikRPA_LowQ2, + kMINERvARW_RikRPA_HighQ2, + kNUISANCEDial_LAST }; int ConvertNUISANCEDial(std::string type); std::string ConvNUISANCEDial(int type); }; #endif diff --git a/src/Reweight/NUISANCEWeightEngine.cxx b/src/Reweight/NUISANCEWeightEngine.cxx index e7344ca..116385c 100644 --- a/src/Reweight/NUISANCEWeightEngine.cxx +++ b/src/Reweight/NUISANCEWeightEngine.cxx @@ -1,129 +1,110 @@ #include "NUISANCEWeightEngine.h" #include "NUISANCEWeightCalcs.h" -NUISANCEWeightEngine::NUISANCEWeightEngine(std::string name) { +#ifdef __MINERVA_RW_ENABLED__ +#ifdef __GENIE_ENABLED__ +#include "MINERvAWeightCalcs.h" +#endif +#endif + +NUISANCEWeightEngine::NUISANCEWeightEngine(std::string name) { // Setup the NUISANCE Reweight engine fCalcName = name; LOG(FIT) << "Setting up NUISANCE Custom RW : " << fCalcName << std::endl; // Load in all Weight Calculations - fWeightCalculators.push_back( new GaussianModeCorr() ); - fWeightCalculators.push_back( new ModeNormCalc() ); + fWeightCalculators.push_back(new GaussianModeCorr()); + fWeightCalculators.push_back(new ModeNormCalc()); + +#ifdef __MINERVA_RW_ENABLED__ +#ifdef __GENIE_ENABLED__ + fWeightCalculators.push_back(new nuisance::reweight::MINERvAReWeight_MEC()); + fWeightCalculators.push_back(new nuisance::reweight::MINERvAReWeight_RES()); + fWeightCalculators.push_back(new nuisance::reweight::RikRPA()); +#endif +#endif + // Set Abs Twk Config fIsAbsTwk = true; - }; - void NUISANCEWeightEngine::IncludeDial(std::string name, double startval) { - // Get First enum int nuisenum = Reweight::ConvDial(name, kCUSTOM); // Setup Maps - fEnumIndex[nuisenum];// = std::vector(0); - fNameIndex[name]; // = std::vector(0); + fEnumIndex[nuisenum]; // = std::vector(0); + fNameIndex[name]; // = std::vector(0); // Split by commas std::vector allnames = GeneralUtils::ParseToStr(name, ","); for (uint i = 0; i < allnames.size(); i++) { std::string singlename = allnames[i]; // Get RW - int singleenum = Reweight::ConvDial(singlename, kCUSTOM); + int singleenum = Reweight::ConvDial(singlename, kCUSTOM); // Fill Maps int index = fValues.size(); fValues.push_back(0.0); fNUISANCEEnums.push_back(singleenum); // Initialize dial std::cout << "Registering " << singlename << " from " << name << std::endl; // Setup index fEnumIndex[nuisenum].push_back(index); fNameIndex[name].push_back(index); } // Set Value if given if (startval != -999.9) { SetDialValue(nuisenum, startval); } }; - void NUISANCEWeightEngine::SetDialValue(int nuisenum, double val) { std::vector indices = fEnumIndex[nuisenum]; for (uint i = 0; i < indices.size(); i++) { fValues[indices[i]] = val; } } void NUISANCEWeightEngine::SetDialValue(std::string name, double val) { std::vector indices = fNameIndex[name]; for (uint i = 0; i < indices.size(); i++) { fValues[indices[i]] = val; } } void NUISANCEWeightEngine::Reconfigure(bool silent) { - for (size_t i = 0; i < fNUISANCEEnums.size(); i++) { - - for (std::vector::iterator calciter = fWeightCalculators.begin(); + for (std::vector::iterator calciter = + fWeightCalculators.begin(); calciter != fWeightCalculators.end(); calciter++) { - - NUISANCEWeightCalc* nuiscalc = static_cast(*calciter); + NUISANCEWeightCalc* nuiscalc = + static_cast(*calciter); if (nuiscalc->IsHandled(fNUISANCEEnums[i])) { nuiscalc->SetDialValue(fNUISANCEEnums[i], fValues[i]); } } } } - - -double NUISANCEWeightEngine::CalcWeight(BaseFitEvt * evt) { +double NUISANCEWeightEngine::CalcWeight(BaseFitEvt* evt) { double rw_weight = 1.0; // Cast as usable class - for (std::vector::iterator iter = fWeightCalculators.begin(); + for (std::vector::iterator iter = + fWeightCalculators.begin(); iter != fWeightCalculators.end(); iter++) { NUISANCEWeightCalc* nuiscalc = static_cast(*iter); rw_weight *= nuiscalc->CalcWeight(evt); } // Return rw_weight return rw_weight; } - - - - - - - - - - - - - - - - - - - - - - - - - - - -