From 04ccf26024f45350f17b722a63542019eb8fdb71 Mon Sep 17 00:00:00 2001 From: Luke Pickering Date: Mon, 10 Aug 2026 11:20:07 +0100 Subject: [PATCH 1/2] Updates to Core 2.6.0: - updates CPM version to quieten warning due to core version update - disables MaCh3 WERROR to allow compilation ignoring may-be-used-uninitialized in core 2.6.0 - duneanaobj dependent build for non-SRProxy separated into new cmake macro in cmake/Modules/fetchduneanaobj.cmake - Fixes various uses of GetFromManager that now require __FILE__/__LINE__ arguments - Moves header inclusion to implementation file in MaCh3DUNEFactory - Updates SampleHandlerBeamFD to use 2.6.0 functional parameters --- Apps/EventRates.cpp | 8 +- Apps/Fit.cpp | 18 +- Apps/LikelihoodScan.cpp | 10 +- Apps/SigmaVariation.cpp | 10 +- CMakeLists.txt | 70 ++---- Samples/CMakeLists.txt | 4 +- Samples/MaCh3DUNEFactory.cpp | 24 +- Samples/MaCh3DUNEFactory.h | 6 - Samples/SampleHandlerBeamFD.cpp | 346 ++++++++++++--------------- Samples/SampleHandlerBeamFD.h | 58 +---- Samples/SampleHandlerBeamND.cpp | 30 +-- Splines/SplineHandlerFactoryDUNE.cpp | 2 +- cmake/Modules/fetchduneanaobj.cmake | 67 ++++++ 13 files changed, 305 insertions(+), 348 deletions(-) create mode 100644 cmake/Modules/fetchduneanaobj.cmake diff --git a/Apps/EventRates.cpp b/Apps/EventRates.cpp index 79abe345..7dee893b 100644 --- a/Apps/EventRates.cpp +++ b/Apps/EventRates.cpp @@ -67,8 +67,8 @@ int main(int argc, char * argv[]) { } } - std::string OutFileName = GetFromManager(FitManager->raw()["General"]["OutputFile"], "EventRatesOutput.root"); - Write1DHistogramsToFile(OutFileName, DUNEHists); + std::string OutFileName = GetFromManager(FitManager->raw()["General"]["OutputFile"], "EventRatesOutput.root", __FILE__, __LINE__); + Write1DHistogramsToFile(OutFileName, DUNEHists); Write1DHistogramsToPdf(OutFileName, DUNEHists); //############################################################################################################################### @@ -77,7 +77,7 @@ int main(int argc, char * argv[]) { MACH3LOG_INFO("========================================================================"); MACH3LOG_INFO("========================================================================"); MACH3LOG_INFO("Oscillation Mode Breakdown:"); - + for(auto handler : samples) { for (int iSample = 0; iSample < handler->GetNSamples(); iSample++) { MACH3LOG_INFO("======================"); @@ -90,7 +90,7 @@ int main(int argc, char * argv[]) { SelecChannel.LowerBound = iOscChan; SelecChannel.UpperBound = iOscChan+1; SelectionVec.push_back(SelecChannel); - + auto Hist = handler->Get1DVarHist(iSample, handler->GetKinVarName(iSample, 0),SelectionVec); MACH3LOG_INFO("{:<20} : {:<20} : {:<20.2f}",handler->GetSampleTitle(iSample),handler->GetFlavourName(iSample, iOscChan),Hist->Integral()); } diff --git a/Apps/Fit.cpp b/Apps/Fit.cpp index 3e6f21a9..f200e864 100644 --- a/Apps/Fit.cpp +++ b/Apps/Fit.cpp @@ -33,11 +33,11 @@ int main(int argc, char * argv[]) { for (auto handler : samples) { for (unsigned iSample = 0; iSample < handler->GetNSamples(); ++iSample) { - + std::string name = handler->GetSampleTitle(iSample); sample_names.push_back(name); TString NameTString = TString(name.c_str()); - + handler->Reweight(); PredictionHistograms.push_back(static_cast(handler->GetMCHist(iSample)->Clone(NameTString+"_DataHist"))); @@ -46,9 +46,9 @@ int main(int argc, char * argv[]) { } else if (handler->GetNDim(iSample) == 2){ handler->AddData(iSample, static_cast(PredictionHistograms.back())); } - + else { - MACH3LOG_ERROR("Unsupported number of dimensions > 2 - Quitting"); + MACH3LOG_ERROR("Unsupported number of dimensions > 2 - Quitting"); throw MaCh3Exception(__FILE__ , __LINE__ ); } @@ -57,16 +57,16 @@ int main(int argc, char * argv[]) { MACH3LOG_INFO("--------------"); } } - + //########################################################################################################### //MCMC auto MaCh3Fitter = MaCh3FitterFactory(FitManager.get()); - bool StartFromPreviousChain = GetFromManager(FitManager->raw()["General"]["StartFromPos"], false); + bool StartFromPreviousChain = GetFromManager(FitManager->raw()["General"]["StartFromPos"], false, __FILE__, __LINE__); //Start chain from random position unless continuing a chain if(!StartFromPreviousChain){ - if (!GetFromManager(FitManager->raw()["General"]["StatOnly"], false)) { + if (!GetFromManager(FitManager->raw()["General"]["StatOnly"], false, __FILE__, __LINE__)) { param_handler->ThrowParameters(); } } @@ -79,13 +79,13 @@ int main(int argc, char * argv[]) { MACH3LOG_INFO("MCMC getting starting position from: {}",PreviousChainPath); MaCh3Fitter->StartFromPreviousFit(PreviousChainPath); } - + //Add samples for(auto Sample : samples){ MaCh3Fitter->AddSampleHandler(Sample); } - + //Run fit MaCh3Fitter->RunMCMC(); diff --git a/Apps/LikelihoodScan.cpp b/Apps/LikelihoodScan.cpp index 0ea14fbe..f32176fd 100644 --- a/Apps/LikelihoodScan.cpp +++ b/Apps/LikelihoodScan.cpp @@ -20,8 +20,8 @@ int main(int argc, char * argv[]) { auto FitManager = MaCh3ManagerFactory(argc, argv); // 1D scan on by default, and 2D off - const bool do_1d_llhscan = GetFromManager(FitManager->raw()["General"]["1DLLHScan"], true); - const bool do_2d_llhscan = GetFromManager(FitManager->raw()["General"]["2DLLHScan"], false); + const bool do_1d_llhscan = GetFromManager(FitManager->raw()["General"]["1DLLHScan"], true, __FILE__, __LINE__); + const bool do_2d_llhscan = GetFromManager(FitManager->raw()["General"]["2DLLHScan"], false, __FILE__, __LINE__); if (!do_1d_llhscan && !do_2d_llhscan) { MACH3LOG_ERROR("Neither 1D or 2D llhscan enabled"); @@ -52,15 +52,15 @@ int main(int argc, char * argv[]) { auto MaCh3Fitter = MaCh3FitterFactory(FitManager.get()); //############################################################################################################################### - //Lets benefit from the core code utilities - + //Lets benefit from the core code utilities + //Add samples to FitterBase for(auto Sample : samples){ MaCh3Fitter->AddSampleHandler(Sample); } MaCh3Fitter->AddSystObj(param_handler.get()); - + if (do_1d_llhscan) { MaCh3Fitter->RunLLHScan(); } diff --git a/Apps/SigmaVariation.cpp b/Apps/SigmaVariation.cpp index 98ee70c7..45650c9c 100644 --- a/Apps/SigmaVariation.cpp +++ b/Apps/SigmaVariation.cpp @@ -42,21 +42,21 @@ int main(int argc, char * argv[]) { std::string OutputFileName = FitManager->raw()["General"]["OutputFile"].as(); TFile* File = TFile::Open(OutputFileName.c_str(),"RECREATE"); - + MACH3LOG_INFO("Starting Variations for covarianceBase Object: {}",param_handler->GetName()); - + int nPars = param_handler->GetNumParams(); for (int iPar=0;iParGetParFancyName(iPar); - double VarInit = param_handler->GetParInit(iPar); + double VarInit = param_handler->GetParPreFit(iPar); double VarSigma = param_handler->GetDiagonalError(iPar); - + MACH3LOG_INFO("\tParameter : {:<30} - Variations around value : {:<10.7f} , in units of 1 Sigma : {:<10.7f}",ParName,VarInit,VarSigma); File->cd(); File->mkdir(ParName.c_str()); File->cd(ParName.c_str()); - + for (size_t iSigVar=0;iSigVarGetLowerBound(iPar)) VarVal = param_handler->GetLowerBound(iPar); diff --git a/CMakeLists.txt b/CMakeLists.txt index 0e468897..b8710a2e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -34,7 +34,7 @@ LIST(APPEND CMAKE_MODULE_PATH ${CMAKE_CURRENT_LIST_DIR}/cmake/Modules) # download CPM.cmake file( DOWNLOAD - https://github.com/cpm-cmake/CPM.cmake/releases/download/v0.40.2/CPM.cmake + https://github.com/cpm-cmake/CPM.cmake/releases/download/v0.43.1/CPM.cmake ${CMAKE_CURRENT_BINARY_DIR}/cmake/CPM.cmake ) include(${CMAKE_CURRENT_BINARY_DIR}/cmake/CPM.cmake) @@ -68,7 +68,7 @@ option(MaCh3_DUNE_USE_SRProxy "Whether to build proxy classes for Standard Recor if(MaCh3_DUNE_USE_SRProxy) find_package(duneanaobj) - + if(NOT duneanaobj_FOUND) CPMFindPackage( NAME duneanaobj @@ -77,78 +77,44 @@ if(MaCh3_DUNE_USE_SRProxy) VERSION 4.0.0 ) endif() - + if(NOT TARGET duneanaobj::all) cmessage(FATAL_ERROR "MaCh3 DUNE Expected dependency target: duneanaobj::all") endif() else() - CPMAddPackage( - NAME duneanaobj - GIT_TAG ${DUNE_ANAOBJ_BRANCH} - GITHUB_REPOSITORY DUNE/duneanaobj - DOWNLOAD_ONLY YES - ) - include_directories(${duneanaobj_SOURCE_DIR}) - - ROOT_GENERATE_DICTIONARY(StandardRecordDict ${duneanaobj_SOURCE_DIR}/duneanaobj/StandardRecord/StandardRecord.h - LINKDEF ${duneanaobj_SOURCE_DIR}/duneanaobj/StandardRecord/classes_def.xml) - - file(GLOB SR_IMPL_FILES "${duneanaobj_SOURCE_DIR}/duneanaobj/StandardRecord/*.cxx") - LIST(APPEND SR_IMPL_FILES ${CMAKE_CURRENT_BINARY_DIR}/StandardRecordDict.cxx) - - add_library(duneanaobj_StandardRecord SHARED ${SR_IMPL_FILES}) - target_link_libraries(duneanaobj_StandardRecord PUBLIC ROOT::MathCore) - - target_include_directories(duneanaobj_StandardRecord PUBLIC - $ - $) - target_include_directories(duneanaobj_StandardRecord PRIVATE - $ #root puts this in the dictionary - ) - - set_target_properties(duneanaobj_StandardRecord PROPERTIES EXPORT_NAME all) - install(TARGETS duneanaobj_StandardRecord EXPORT mach3dune-targets DESTINATION lib) - - file(GLOB SR_HEADER_FILES "${duneanaobj_SOURCE_DIR}/duneanaobj/StandardRecord/*.h") - - install(FILES ${SR_HEADER_FILES} DESTINATION include/duneanaobj/StandardRecord) - install(FILES - ${CMAKE_CURRENT_BINARY_DIR}/libStandardRecordDict_rdict.pcm - ${CMAKE_CURRENT_BINARY_DIR}/libStandardRecordDict.rootmap - DESTINATION lib) - - add_library(duneanaobj::all ALIAS duneanaobj_StandardRecord) - + include(fetchduneanaobj) + fetchduneanaobj(${DUNE_ANAOBJ_BRANCH}) endif() ################################## MaCh3 ###################################### #If MaCh3 was sourced find it, otherwise use CPM SET(MaCh3_FOUND FALSE) -find_package(MaCh3 2.5.1 EXACT QUIET) -set(MaCh3_CORE_BRANCH v2.5.1 CACHE STRING "Specify the MaCh3 core branch to use") +find_package(MaCh3 2.6.0 EXACT QUIET) +set(MaCh3_CORE_BRANCH v2.6.0 CACHE STRING "Specify the MaCh3 core branch to use") message(STATUS "Using MaCh3_CORE_BRANCH: ${MaCh3_CORE_BRANCH}") if(NOT MaCh3_FOUND) cmessage(STATUS "Didn't find MaCh3, attempting to use built in MaCh3") - + if(NOT DEFINED MaCh3_GPU_ENABLED) set(MaCh3_GPU_ENABLED ON) endif() - + if(NOT DEFINED MaCh3_DEBUG_ENABLED) set(MaCh3_DEBUG_ENABLED FALSE) endif() - + if(NOT DEFINED MaCh3_MULTITHREAD_ENABLED) set(MaCh3_MULTITHREAD_ENABLED TRUE) endif() - + # Options list construction set(MaCh3_OPTIONS "MaCh3_GPU_ENABLED ${MaCh3_GPU_ENABLED}" "MaCh3_DEBUG_ENABLED ${MaCh3_DEBUG_ENABLED}" "MaCh3_MULTITHREAD_ENABLED ${MaCh3_MULTITHREAD_ENABLED}" "MaCh3_CORE_BRANCH ${MaCh3_CORE_BRANCH}" + "MaCh3_WERROR_ENABLED OFF" #required due to small uninitialized issue in MaCh3Core 2.6.0, revert when fixed ) # Add LOG_LEVEL if defined @@ -159,7 +125,7 @@ if(NOT MaCh3_FOUND) list(APPEND MaCh3_OPTIONS "NuFastLinear_ENABLED TRUE") list(APPEND MaCh3_OPTIONS "CUDAProb3_ENABLED TRUE") list(APPEND MaCh3_OPTIONS "CUDAProb3Linear_ENABLED TRUE") - + CPMAddPackage( NAME MaCh3 GIT_TAG ${MaCh3_CORE_BRANCH} @@ -171,7 +137,7 @@ else() ##KS: This ensure that all executables that are in core will be moved FILE(GLOB MaCh3Exe $ENV{MaCh3_ROOT}/Diagnostics/*) FILE(COPY ${MaCh3Exe} DESTINATION ${CMAKE_BINARY_DIR}/Diagnostics/) - + FILE(GLOB MaCh3Exe $ENV{MaCh3_ROOT}/plotting/*) FILE(COPY ${MaCh3Exe} DESTINATION ${CMAKE_BINARY_DIR}/plotting/) endif() @@ -286,16 +252,16 @@ if(NOT DEFINED MaCh3_PREFIX OR MaCh3_PREFIX STREQUAL "") cmessage(STATUS "Detected MaCh3 from target include dirs: ${MaCh3_PREFIX}") endif() endif() - + # Final validation if(NOT DEFINED MaCh3_PREFIX OR MaCh3_PREFIX STREQUAL "") cmessage(FATAL_ERROR "Could not determine MaCh3_PREFIX. Please set it manually with -DMaCh3_PREFIX=") endif() - + if(NOT EXISTS "${MaCh3_PREFIX}") cmessage(FATAL_ERROR "MaCh3_PREFIX does not exist: ${MaCh3_PREFIX}") endif() - + # Cache it for use in subdirectories set(MaCh3_PREFIX "${MaCh3_PREFIX}" CACHE PATH "MaCh3 installation or source directory" FORCE) cmessage(STATUS "Set MaCh3_PREFIX: ${MaCh3_PREFIX}") @@ -350,4 +316,4 @@ install(FILES ${CMAKE_BINARY_DIR}/MaCh3DUNEConfigVersion.cmake DESTINATION ${CMAKE_INSTALL_PREFIX}/lib/cmake/MaCh3DUNE) -include(mach3dune-config) \ No newline at end of file +include(mach3dune-config) diff --git a/Samples/CMakeLists.txt b/Samples/CMakeLists.txt index 7a06c74e..fd1160d5 100644 --- a/Samples/CMakeLists.txt +++ b/Samples/CMakeLists.txt @@ -21,8 +21,8 @@ set_target_properties(SamplesDUNE PROPERTIES PUBLIC_HEADER "${HEADERS}" EXPORT_NAME SamplesDUNE) -target_link_libraries(SamplesDUNE PUBLIC SplinesDUNE duneanaobj::all MaCh3::All MaCh3DUNECompilerOptions) -target_link_libraries(SamplesDUNE PRIVATE DUNEMaCh3Warnings) +target_link_libraries(SamplesDUNE PUBLIC SplinesDUNE MaCh3::All MaCh3DUNECompilerOptions) +target_link_libraries(SamplesDUNE PRIVATE duneanaobj::StandardRecord DUNEMaCh3Warnings) target_include_directories(SamplesDUNE PUBLIC $ diff --git a/Samples/MaCh3DUNEFactory.cpp b/Samples/MaCh3DUNEFactory.cpp index 91232d63..84d354d9 100644 --- a/Samples/MaCh3DUNEFactory.cpp +++ b/Samples/MaCh3DUNEFactory.cpp @@ -1,5 +1,11 @@ #include "Samples/MaCh3DUNEFactory.h" +// DUNE Handlers +#include "Samples/SampleHandlerBeamFD.h" +#include "Samples/SampleHandlerBeamND.h" +#include "Samples/SampleHandlerBeamNDGAr.h" +#include "Samples/SampleHandlerAtm.h" + // ############################################################### SampleHandlerBase* GetMaCh3DuneInstance(std::string SampleType, std::string SampleConfig, std::unique_ptr& param_handler, const std::shared_ptr& BeamOscillator_, const std::shared_ptr& AtmOscillator_, BeamNDCov beamNDCov) { // ############################################################### @@ -7,12 +13,12 @@ SampleHandlerBase* GetMaCh3DuneInstance(std::string SampleType, std::string Samp if (SampleType == "BeamFD") { Sample = new SampleHandlerBeamFD(SampleConfig, param_handler.get(), BeamOscillator_); } else if (SampleType == "BeamND") { - + if (beamNDCov.NDCov_FHC == nullptr || beamNDCov.NDCov_RHC == nullptr || beamNDCov.NDCov_all == nullptr) { MACH3LOG_ERROR("NDCov objects are not defined"); throw MaCh3Exception(__FILE__, __LINE__); } - Sample = new SampleHandlerBeamND(SampleConfig, param_handler.get(), beamNDCov); + Sample = new SampleHandlerBeamND(SampleConfig, param_handler.get(), beamNDCov); } else if (SampleType == "Atm") { Sample = new SampleHandlerAtm(SampleConfig, param_handler.get(), AtmOscillator_); } else if (SampleType == "BeamNDGAr") { @@ -22,7 +28,7 @@ SampleHandlerBase* GetMaCh3DuneInstance(std::string SampleType, std::string Samp MACH3LOG_ERROR("Invalid SampleType: {} defined in {}", SampleType, SampleConfig); throw MaCh3Exception(__FILE__, __LINE__); } - + return Sample; } @@ -58,7 +64,7 @@ BeamNDCov SetupBeamNDCov(std::unique_ptr &FitManager) } std::string NDCovMatrixFile = FitManager->raw()["General"]["Systematics"]["NDCovFile"].as(); - bool useCombinedNDCov = GetFromManager(FitManager->raw()["General"]["Systematics"]["UseCombinedNDCov"], true); + bool useCombinedNDCov = GetFromManager(FitManager->raw()["General"]["Systematics"]["UseCombinedNDCov"], true, __FILE__, __LINE__); auto NDCovFile = M3::Open(NDCovMatrixFile, "READ", __FILE__, __LINE__); @@ -98,7 +104,7 @@ std::vector MaCh3DuneSampleFactory(std::unique_ptr throw MaCh3Exception(__FILE__, __LINE__); } - + // ========================================================== // Setup oscillation handlers auto AtmOscHandler = SetupOscillationHandler(FitManager, param_handler, "ATM", "ATM"); @@ -122,7 +128,7 @@ std::vector MaCh3DuneSampleFactory(std::unique_ptr std::string SampleType = tempSampleManager->raw()["SampleHandlerName"].as(); auto sample = GetMaCh3DuneInstance(SampleType, DUNESampleConfigs[Sample_i], param_handler, BeamOscHandler, AtmOscHandler, beamNDCov); - + #if DEBUG_DUNE_WEIGHTS==1 // Pure for debugging, lets us set which weights we don't want via the manager sample->setWeightSwitchOffVector(FitManager->getWeightSwitchOffVector()); @@ -137,11 +143,11 @@ std::vector MaCh3DuneSampleFactory(std::unique_ptr std::pair, std::vector> MaCh3DuneFactory(std::unique_ptr &FitManager) { /// Generates a MaCh3 DUNE instance auto param_handler = MaCh3CovarianceFactory(FitManager.get(), "Xsec"); - + if (CheckNodeExists(FitManager->raw(), "General", "OscillationParameters")){ auto oscpars = Get>(FitManager->raw()["General"]["OscillationParameters"], __FILE__, __LINE__); param_handler->SetGroupOnlyParameters("Osc", oscpars); - } + } auto samples = MaCh3DuneSampleFactory(FitManager, param_handler); @@ -150,4 +156,4 @@ std::pair, std::vector& Oscillator_) : SampleHandlerBase(mc_version_, ParHandler_, Oscillator_) { KinematicParameters = &KinematicParametersDUNE; ReversedKinematicParameters = &ReversedKinematicParametersDUNE; - + Initialise(); downsamplingStep = 1; @@ -29,14 +29,14 @@ void SampleHandlerBeamFD::Init() { MACH3LOG_INFO("- iselike: {}", beamFDSampleDetails[i].iselike); } - downsamplingStep = GetFromManager(SampleManager->raw()["DownsamplingStep"], 1); + downsamplingStep = GetFromManager(SampleManager->raw()["DownsamplingStep"], 1, __FILE__ , __LINE__); if (downsamplingStep == 0) { throw MaCh3Exception(__FILE__, __LINE__, "Downsampling step cannot be zero. Please set it to a positive integer in the Beam FD sample config file." ); - } + } MACH3LOG_INFO("Beam FD downsampling step: {}", downsamplingStep); - + MACH3LOG_INFO("-------------------------------------------------------------------"); } @@ -46,7 +46,7 @@ void SampleHandlerBeamFD::InititialiseData() // ************************************************ // Reweight MC to match Reweight(); - // set asimov data + // set asimov data for (int iSample = 0; iSample < GetNSamples(); iSample++) { AddData(iSample, GetMCArray(iSample)); @@ -66,147 +66,147 @@ void SampleHandlerBeamFD::SetupSplines() { MACH3LOG_INFO("Found {} splines for this sample so I will not load or evaluate splines", ParHandler->GetNumParamsFromSampleName(SampleHandlerName, kSpline)); SplineHandler = nullptr; } - + return; } // === HH: Functional parameters === -void SampleHandlerBeamFD::TotalEScale(const double * par, std::size_t iEvent) { +void TotalEScale(double const & par_val, dunemc_beamfd &ev) { // Total energy scale uncertainties for anything but CC Numu, see: // https://github.com/DUNE/lblpwgtools/blob/3d475f50a998fbfa6266df9a0c4eb3056c0cdfe5/CAFAna/Systs/EnergySysts.h#L39 - dunemcSamples[iEvent].rw_erec_shifted += (*par) * dunemcSamples[iEvent].rw_erec_had; + ev.rw_erec_shifted += par_val * ev.rw_erec_had; } -void SampleHandlerBeamFD::TotalEScaleNotCCNumu(const double * par, std::size_t iEvent) { +void TotalEScaleNotCCNumu(double const & par_val, dunemc_beamfd &ev) { // A special case for Not (CC Numu), where we also scale Erec by lepton energy // Since we reconstruct muon energy in a different way, see: // https://github.com/DUNE/lblpwgtools/blob/3d475f50a998fbfa6266df9a0c4eb3056c0cdfe5/CAFAna/Systs/EnergySysts.h#L39 - dunemcSamples[iEvent].rw_erec_shifted += (*par) * dunemcSamples[iEvent].rw_erec_lep; + ev.rw_erec_shifted += par_val * ev.rw_erec_lep; } -void SampleHandlerBeamFD::TotalEScaleSqrt(const double * par, std::size_t iEvent) { - dunemcSamples[iEvent].rw_erec_shifted += (*par) * dunemcSamples[iEvent].rw_erec_had * dunemcSamples[iEvent].rw_erec_had_sqrt; +void TotalEScaleSqrt(double const & par_val, dunemc_beamfd &ev) { + ev.rw_erec_shifted += par_val * ev.rw_erec_had * ev.rw_erec_had_sqrt; } -void SampleHandlerBeamFD::TotalEScaleSqrtNotCCNumu(const double * par, std::size_t iEvent) { +void TotalEScaleSqrtNotCCNumu(double const & par_val, dunemc_beamfd &ev) { // See comments in TotalEScaleNotCCNumu - dunemcSamples[iEvent].rw_erec_shifted += (*par) * dunemcSamples[iEvent].rw_erec_lep * dunemcSamples[iEvent].rw_erec_lep_sqrt; + ev.rw_erec_shifted += par_val * ev.rw_erec_lep * ev.rw_erec_lep_sqrt; } -void SampleHandlerBeamFD::TotalEScaleInvSqrt(const double * par, std::size_t iEvent) { +void TotalEScaleInvSqrt(double const & par_val, dunemc_beamfd &ev) { // Erec/sqrt(Erec) = sqrt(Erec) - dunemcSamples[iEvent].rw_erec_shifted += (*par) * dunemcSamples[iEvent].rw_erec_had_sqrt; + ev.rw_erec_shifted += par_val * ev.rw_erec_had_sqrt; } -void SampleHandlerBeamFD::TotalEScaleInvSqrtNotCCNumu(const double * par, std::size_t iEvent) { +void TotalEScaleInvSqrtNotCCNumu(double const & par_val, dunemc_beamfd &ev) { // See comments in TotalEScaleNotCCNumu - dunemcSamples[iEvent].rw_erec_shifted += (*par) * dunemcSamples[iEvent].rw_erec_lep_sqrt; + ev.rw_erec_shifted += par_val * ev.rw_erec_lep_sqrt; } -void SampleHandlerBeamFD::HadEScale(const double * par, std::size_t iEvent) { - dunemcSamples[iEvent].rw_erec_shifted += (*par) * dunemcSamples[iEvent].rw_sum_ehad; +void HadEScale(double const & par_val, dunemc_beamfd &ev) { + ev.rw_erec_shifted += par_val * ev.rw_sum_ehad; } -void SampleHandlerBeamFD::HadEScaleSqrt(const double * par, std::size_t iEvent) { - dunemcSamples[iEvent].rw_erec_shifted += (*par) * dunemcSamples[iEvent].rw_sum_ehad * dunemcSamples[iEvent].rw_sum_ehad_sqrt; +void HadEScaleSqrt(double const & par_val, dunemc_beamfd &ev) { + ev.rw_erec_shifted += par_val * ev.rw_sum_ehad * ev.rw_sum_ehad_sqrt; } -void SampleHandlerBeamFD::HadEScaleInvSqrt(const double * par, std::size_t iEvent) { - dunemcSamples[iEvent].rw_erec_shifted += (*par) * dunemcSamples[iEvent].rw_sum_ehad_sqrt; +void HadEScaleInvSqrt(double const & par_val, dunemc_beamfd &ev) { + ev.rw_erec_shifted += par_val * ev.rw_sum_ehad_sqrt; } -void SampleHandlerBeamFD::MuEScale(const double * par, std::size_t iEvent) { +void MuEScale(double const & par_val, dunemc_beamfd &ev) { // HH TODO: Functionally this is the same as TotalEScaleNotCCNumu, not sure if this function is even needed - TotalEScaleNotCCNumu(par, iEvent); + TotalEScaleNotCCNumu(par_val, ev); } -void SampleHandlerBeamFD::MuEScaleSqrt(const double * par, std::size_t iEvent) { +void MuEScaleSqrt(double const & par_val, dunemc_beamfd &ev) { // See comments in MuEScale - TotalEScaleSqrtNotCCNumu(par, iEvent); + TotalEScaleSqrtNotCCNumu(par_val, ev); } -void SampleHandlerBeamFD::MuEScaleInvSqrt(const double * par, std::size_t iEvent) { +void MuEScaleInvSqrt(double const & par_val, dunemc_beamfd &ev) { // See comments in MuEScale - TotalEScaleInvSqrtNotCCNumu(par, iEvent); + TotalEScaleInvSqrtNotCCNumu(par_val, ev); } -void SampleHandlerBeamFD::NEScale(const double * par, std::size_t iEvent) { - dunemcSamples[iEvent].rw_erec_shifted += (*par) * dunemcSamples[iEvent].rw_eRecoN; +void NEScale(double const & par_val, dunemc_beamfd &ev) { + ev.rw_erec_shifted += par_val * ev.rw_eRecoN; } -void SampleHandlerBeamFD::NEScaleSqrt(const double * par, std::size_t iEvent) { - dunemcSamples[iEvent].rw_erec_shifted += (*par) * dunemcSamples[iEvent].rw_eRecoN * dunemcSamples[iEvent].rw_eRecoN_sqrt; +void NEScaleSqrt(double const & par_val, dunemc_beamfd &ev) { + ev.rw_erec_shifted += par_val * ev.rw_eRecoN * ev.rw_eRecoN_sqrt; } -void SampleHandlerBeamFD::NEScaleInvSqrt(const double * par, std::size_t iEvent) { - dunemcSamples[iEvent].rw_erec_shifted += (*par) * dunemcSamples[iEvent].rw_eRecoN_sqrt; +void NEScaleInvSqrt(double const & par_val, dunemc_beamfd &ev) { + ev.rw_erec_shifted += par_val * ev.rw_eRecoN_sqrt; } -void SampleHandlerBeamFD::EMEScale(const double * par, std::size_t iEvent) { - dunemcSamples[iEvent].rw_erec_shifted += (*par) * dunemcSamples[iEvent].rw_eRecoPi0; +void EMEScale(double const & par_val, dunemc_beamfd &ev) { + ev.rw_erec_shifted += par_val * ev.rw_eRecoPi0; } -void SampleHandlerBeamFD::EMEScaleCCNue(const double * par, std::size_t iEvent) { +void EMEScaleCCNue(double const & par_val, dunemc_beamfd &ev) { // Again this is the same as TotalEScaleNotCCNumu, not sure if this function is needed - TotalEScaleNotCCNumu(par, iEvent); + TotalEScaleNotCCNumu(par_val, ev); } -void SampleHandlerBeamFD::EMEScaleSqrt(const double * par, std::size_t iEvent) { - dunemcSamples[iEvent].rw_erec_shifted += (*par) * dunemcSamples[iEvent].rw_eRecoPi0 * dunemcSamples[iEvent].rw_eRecoPi0_sqrt; +void EMEScaleSqrt(double const & par_val, dunemc_beamfd &ev) { + ev.rw_erec_shifted += par_val * ev.rw_eRecoPi0 * ev.rw_eRecoPi0_sqrt; } -void SampleHandlerBeamFD::EMEScaleSqrtCCNue(const double * par, std::size_t iEvent) { +void EMEScaleSqrtCCNue(double const & par_val, dunemc_beamfd &ev) { // See comments in EMEScaleCCNue - TotalEScaleSqrtNotCCNumu(par, iEvent); + TotalEScaleSqrtNotCCNumu(par_val, ev); } -void SampleHandlerBeamFD::EMEScaleInvSqrt(const double * par, std::size_t iEvent) { - dunemcSamples[iEvent].rw_erec_shifted += (*par) * dunemcSamples[iEvent].rw_eRecoPi0_sqrt; +void EMEScaleInvSqrt(double const & par_val, dunemc_beamfd &ev) { + ev.rw_erec_shifted += par_val * ev.rw_eRecoPi0_sqrt; } -void SampleHandlerBeamFD::EMEScaleInvSqrtCCNue(const double * par, std::size_t iEvent) { +void EMEScaleInvSqrtCCNue(double const & par_val, dunemc_beamfd &ev) { // See comments in EMEScaleCCNue - TotalEScaleInvSqrtNotCCNumu(par, iEvent); + TotalEScaleInvSqrtNotCCNumu(par_val, ev); } -void SampleHandlerBeamFD::HadRes(const double * par, std::size_t iEvent) { +void HadRes(double const & par_val, dunemc_beamfd &ev) { // True sum - reco sum - dunemcSamples[iEvent].rw_erec_shifted += (*par) * (dunemcSamples[iEvent].rw_eP - + dunemcSamples[iEvent].rw_ePip - + dunemcSamples[iEvent].rw_ePim - - dunemcSamples[iEvent].rw_sum_ehad); + ev.rw_erec_shifted += par_val * (ev.rw_eP + + ev.rw_ePip + + ev.rw_ePim + - ev.rw_sum_ehad); } -void SampleHandlerBeamFD::MuRes(const double * par, std::size_t iEvent) { +void MuRes(double const & par_val, dunemc_beamfd &ev) { // True muon energy - reco muon energy - dunemcSamples[iEvent].rw_erec_shifted += (*par) * (dunemcSamples[iEvent].rw_LepE - dunemcSamples[iEvent].rw_erec_lep); + ev.rw_erec_shifted += par_val * (ev.rw_LepE - ev.rw_erec_lep); } -void SampleHandlerBeamFD::NRes(const double * par, std::size_t iEvent) { +void NRes(double const & par_val, dunemc_beamfd &ev) { // True neutron energy - reco neutron energy - dunemcSamples[iEvent].rw_erec_shifted += (*par) * (dunemcSamples[iEvent].rw_eN - dunemcSamples[iEvent].rw_eRecoN); + ev.rw_erec_shifted += par_val * (ev.rw_eN - ev.rw_eRecoN); } -void SampleHandlerBeamFD::EMRes(const double * par, std::size_t iEvent) { +void EMRes(double const & par_val, dunemc_beamfd &ev) { // True pi0 energy - reco pi0 energy - dunemcSamples[iEvent].rw_erec_shifted += (*par) * (dunemcSamples[iEvent].rw_ePi0 - dunemcSamples[iEvent].rw_eRecoPi0); + ev.rw_erec_shifted += par_val * (ev.rw_ePi0 - ev.rw_eRecoPi0); } -void SampleHandlerBeamFD::EMResCCNue(const double * par, std::size_t iEvent) { +void EMResCCNue(double const & par_val, dunemc_beamfd &ev) { // This is the same as MuRes, again not sure if this function is needed - MuRes(par, iEvent); + MuRes(par_val, ev); } -void SampleHandlerBeamFD::RecoCVNNumu(const double * par, std::size_t iEvent) { +void RecoCVNNumu(double const & par_val, dunemc_beamfd &ev) { // CVN numu uncertainty - dunemcSamples[iEvent].rw_cvnnumu_shifted += (*par); + ev.rw_cvnnumu_shifted += par_val; } -void SampleHandlerBeamFD::RecoCVNNue(const double * par, std::size_t iEvent) { +void RecoCVNNue(double const & par_val, dunemc_beamfd &ev) { // CVN nue uncertainty - dunemcSamples[iEvent].rw_cvnnue_shifted += (*par); + ev.rw_cvnnue_shifted += par_val; } void SampleHandlerBeamFD::RegisterFunctionalParameters() { @@ -214,119 +214,83 @@ void SampleHandlerBeamFD::RegisterFunctionalParameters() { // This function manually populates the map of functional parameters // Maps the name of the functional parameter to the pointer of the function - // This is the part where we manually enter things - // A lambda function has to be used so we can refer to a non-static member function - RegisterIndividualFunctionalParameter("TotalEScaleFD", - kTotalEScale, - [this](const double * par, std::size_t iEvent) { this->TotalEScale(par, iEvent); }); + RegisterIndividualFunctionalParameter(dunemcSamples, "TotalEScaleFD", + TotalEScale); - RegisterIndividualFunctionalParameter("TotalEScaleNotCCNumuFD", - kTotalEScaleNotCCNumu, - [this](const double * par, std::size_t iEvent) { this->TotalEScaleNotCCNumu(par, iEvent); }); + RegisterIndividualFunctionalParameter(dunemcSamples, "TotalEScaleNotCCNumuFD", + TotalEScaleNotCCNumu); - RegisterIndividualFunctionalParameter("TotalEScaleSqrtFD", - kTotalEScaleSqrt, - [this](const double * par, std::size_t iEvent) { this->TotalEScaleSqrt(par, iEvent); }); + RegisterIndividualFunctionalParameter(dunemcSamples, "TotalEScaleSqrtFD", + TotalEScaleSqrt); - RegisterIndividualFunctionalParameter("TotalEScaleSqrtNotCCNumuFD", - kTotalEScaleSqrtNotCCNumu, - [this](const double * par, std::size_t iEvent) { this->TotalEScaleSqrtNotCCNumu(par, iEvent); }); + RegisterIndividualFunctionalParameter( + dunemcSamples, "TotalEScaleSqrtNotCCNumuFD", TotalEScaleSqrtNotCCNumu); - RegisterIndividualFunctionalParameter("TotalEScaleInvSqrtFD", - kTotalEScaleInvSqrt, - [this](const double * par, std::size_t iEvent) { this->TotalEScaleInvSqrt(par, iEvent); }); + RegisterIndividualFunctionalParameter(dunemcSamples, "TotalEScaleInvSqrtFD", + TotalEScaleInvSqrt); - RegisterIndividualFunctionalParameter("TotalEScaleInvSqrtNotCCNumuFD", - kTotalEScaleInvSqrtNotCCNumu, - [this](const double * par, std::size_t iEvent) { this->TotalEScaleInvSqrtNotCCNumu(par, iEvent); }); + RegisterIndividualFunctionalParameter(dunemcSamples, + "TotalEScaleInvSqrtNotCCNumuFD", + TotalEScaleInvSqrtNotCCNumu); - RegisterIndividualFunctionalParameter("HadEScaleFD", - kHadEScale, - [this](const double * par, std::size_t iEvent) { this->HadEScale(par, iEvent); }); + RegisterIndividualFunctionalParameter(dunemcSamples, "HadEScaleFD", + HadEScale); - RegisterIndividualFunctionalParameter("HadEScaleSqrtFD", - kHadEScaleSqrt, - [this](const double * par, std::size_t iEvent) { this->HadEScaleSqrt(par, iEvent); }); + RegisterIndividualFunctionalParameter(dunemcSamples, "HadEScaleSqrtFD", + HadEScaleSqrt); - RegisterIndividualFunctionalParameter("HadEScaleInvSqrtFD", - kHadEScaleInvSqrt, - [this](const double * par, std::size_t iEvent) { this->HadEScaleInvSqrt(par, iEvent); }); + RegisterIndividualFunctionalParameter(dunemcSamples, "HadEScaleInvSqrtFD", + HadEScaleInvSqrt); - RegisterIndividualFunctionalParameter("MuEScaleFD", - kMuEScale, - [this](const double * par, std::size_t iEvent) { this->MuEScale(par, iEvent); }); + RegisterIndividualFunctionalParameter(dunemcSamples, "MuEScaleFD", MuEScale); - RegisterIndividualFunctionalParameter("MuEScaleSqrtFD", - kMuEScaleSqrt, - [this](const double * par, std::size_t iEvent) { this->MuEScaleSqrt(par, iEvent); }); + RegisterIndividualFunctionalParameter(dunemcSamples, "MuEScaleSqrtFD", + MuEScaleSqrt); - RegisterIndividualFunctionalParameter("MuEScaleInvSqrtFD", - kMuEScaleInvSqrt, - [this](const double * par, std::size_t iEvent) { this->MuEScaleInvSqrt(par, iEvent); }); + RegisterIndividualFunctionalParameter(dunemcSamples, "MuEScaleInvSqrtFD", + MuEScaleInvSqrt); - RegisterIndividualFunctionalParameter("NEScaleFD", - kNEScale, - [this](const double * par, std::size_t iEvent) { this->NEScale(par, iEvent); }); + RegisterIndividualFunctionalParameter(dunemcSamples, "NEScaleFD", NEScale); - RegisterIndividualFunctionalParameter("NEScaleSqrtFD", - kNEScaleSqrt, - [this](const double * par, std::size_t iEvent) { this->NEScaleSqrt(par, iEvent); }); + RegisterIndividualFunctionalParameter(dunemcSamples, "NEScaleSqrtFD", + NEScaleSqrt); - RegisterIndividualFunctionalParameter("NEScaleInvSqrtFD", - kNEScaleInvSqrt, - [this](const double * par, std::size_t iEvent) { this->NEScaleInvSqrt(par, iEvent); }); + RegisterIndividualFunctionalParameter(dunemcSamples, "NEScaleInvSqrtFD", + NEScaleInvSqrt); - RegisterIndividualFunctionalParameter("EMEScaleFD", - kEMEScale, - [this](const double * par, std::size_t iEvent) { this->EMEScale(par, iEvent); }); + RegisterIndividualFunctionalParameter(dunemcSamples, "EMEScaleFD", EMEScale); - RegisterIndividualFunctionalParameter("EMEScaleCCNueFD", - kEMEScaleCCNue, - [this](const double * par, std::size_t iEvent) { this->EMEScaleCCNue(par, iEvent); }); + RegisterIndividualFunctionalParameter(dunemcSamples, "EMEScaleCCNueFD", + EMEScaleCCNue); - RegisterIndividualFunctionalParameter("EMEScaleSqrtFD", - kEMEScaleSqrt, - [this](const double * par, std::size_t iEvent) { this->EMEScaleSqrt(par, iEvent); }); + RegisterIndividualFunctionalParameter(dunemcSamples, "EMEScaleSqrtFD", + EMEScaleSqrt); - RegisterIndividualFunctionalParameter("EMEScaleSqrtCCNueFD", - kEMEScaleSqrtCCNue, - [this](const double * par, std::size_t iEvent) { this->EMEScaleSqrtCCNue(par, iEvent); }); + RegisterIndividualFunctionalParameter(dunemcSamples, "EMEScaleSqrtCCNueFD", + EMEScaleSqrtCCNue); - RegisterIndividualFunctionalParameter("EMEScaleInvSqrtFD", - kEMEScaleInvSqrt, - [this](const double * par, std::size_t iEvent) { this->EMEScaleInvSqrt(par, iEvent); }); + RegisterIndividualFunctionalParameter(dunemcSamples, "EMEScaleInvSqrtFD", + EMEScaleInvSqrt); - RegisterIndividualFunctionalParameter("EMEScaleInvSqrtCCNueFD", - kEMEScaleInvSqrtCCNue, - [this](const double * par, std::size_t iEvent) { this->EMEScaleInvSqrtCCNue(par, iEvent); }); + RegisterIndividualFunctionalParameter(dunemcSamples, "EMEScaleInvSqrtCCNueFD", + EMEScaleInvSqrtCCNue); - RegisterIndividualFunctionalParameter("HadResFD", - kHadRes, - [this](const double * par, std::size_t iEvent) { this->HadRes(par, iEvent); }); + RegisterIndividualFunctionalParameter(dunemcSamples, "HadResFD", HadRes); - RegisterIndividualFunctionalParameter("MuResFD", - kMuRes, - [this](const double * par, std::size_t iEvent) { this->MuRes(par, iEvent); }); + RegisterIndividualFunctionalParameter(dunemcSamples, "MuResFD", MuRes); - RegisterIndividualFunctionalParameter("NResFD", - kNRes, - [this](const double * par, std::size_t iEvent) { this->NRes(par, iEvent); }); + RegisterIndividualFunctionalParameter(dunemcSamples, "NResFD", NRes); - RegisterIndividualFunctionalParameter("EMResFD", - kEMRes, - [this](const double * par, std::size_t iEvent) { this->EMRes(par, iEvent); }); + RegisterIndividualFunctionalParameter(dunemcSamples, "EMResFD", EMRes); - RegisterIndividualFunctionalParameter("EMResCCNueFD", - kEMResCCNue, - [this](const double * par, std::size_t iEvent) { this->EMResCCNue(par, iEvent); }); + RegisterIndividualFunctionalParameter(dunemcSamples, "EMResCCNueFD", + EMResCCNue); - RegisterIndividualFunctionalParameter("RecoCVNNumuFD", - kRecoCVNNumu, - [this](const double * par, std::size_t iEvent) { this->RecoCVNNumu(par, iEvent); }); + RegisterIndividualFunctionalParameter(dunemcSamples, "RecoCVNNumuFD", + RecoCVNNumu); - RegisterIndividualFunctionalParameter("RecoCVNNueFD", - kRecoCVNNue, - [this](const double * par, std::size_t iEvent) { this->RecoCVNNue(par, iEvent); }); + RegisterIndividualFunctionalParameter(dunemcSamples, "RecoCVNNueFD", + RecoCVNNue); MACH3LOG_INFO("Finished registering functional parameters"); } @@ -352,11 +316,11 @@ void SampleHandlerBeamFD::AddAdditionalWeightPointers() { int SampleHandlerBeamFD::SetupExperimentMC() { // dunemc_base *duneobj = &(dunemcSamples[iSample]); - + MACH3LOG_INFO("-------------------------------------------------------------------"); TChain* _data = new TChain("caf"); // Maps the file index within the TChain (GetTreeNumber()) to its sample index and - // per-file norm values. + // per-file norm values. std::vector fileIndexToSample; std::vector> fileIndexToNorm; // [norm_s, pot_s] for (size_t iSample=0; iSampleGet("norm"); @@ -384,8 +348,8 @@ int SampleHandlerBeamFD::SetupExperimentMC() { } } } - - + + if(_data){ MACH3LOG_INFO("Number of entries in TChain: {}", _data->GetEntries()); } @@ -415,7 +379,7 @@ int SampleHandlerBeamFD::SetupExperimentMC() { double _vtx_z; //Truth Variables - int _mode; + int _mode; double _ev; double _LepE; double _eP; @@ -427,7 +391,7 @@ int SampleHandlerBeamFD::SetupExperimentMC() { int _isCC; int _nuPDGunosc; int _nuPDG; - + _data->SetBranchStatus("*", 0); _data->SetBranchStatus("Ev", 1); _data->SetBranchAddress("Ev", &_ev); @@ -487,7 +451,7 @@ int SampleHandlerBeamFD::SetupExperimentMC() { _data->SetBranchStatus("vtx_y", 1); _data->SetBranchAddress("vtx_y", &_vtx_y); _data->SetBranchStatus("vtx_z", 1); - _data->SetBranchAddress("vtx_z", &_vtx_z); + _data->SetBranchAddress("vtx_z", &_vtx_z); size_t nEntries = static_cast(_data->GetEntries()); size_t nDownsampledEntries = nEntries / downsamplingStep; @@ -508,7 +472,7 @@ int SampleHandlerBeamFD::SetupExperimentMC() { negative_counts["rw_eRecoN"] = 0; negative_counts["rw_eRecoPi0"] = 0; negative_counts["rw_sum_ehad"] = 0; - + //FILL DUNE STRUCT for (unsigned int i = 0; i < nDownsampledEntries; ++i) { // Loop through tree _data->GetEntry(i * downsamplingStep); @@ -527,33 +491,33 @@ int SampleHandlerBeamFD::SetupExperimentMC() { // POT stuff dunemcSamples[i].norm_s = fileIndexToNorm[static_cast(_data->GetTreeNumber())][0]; // Norm in sample dunemcSamples[i].pot_s = fileIndexToNorm[static_cast(_data->GetTreeNumber())][1] * downsamplingStep; // POT in sample - + dunemcSamples[i].rw_cvnnumu = (_cvnnumu); dunemcSamples[i].rw_cvnnue = (_cvnnue); - dunemcSamples[i].rw_cvnnumu_shifted = (_cvnnumu); + dunemcSamples[i].rw_cvnnumu_shifted = (_cvnnumu); dunemcSamples[i].rw_cvnnue_shifted = (_cvnnue); if (iselike_temp) { dunemcSamples[i].rw_erec = (_erec_nue); - dunemcSamples[i].rw_erec_shifted = (_erec_nue); + dunemcSamples[i].rw_erec_shifted = (_erec_nue); dunemcSamples[i].rw_erec_had = (_erec_had_nue); dunemcSamples[i].rw_erec_lep = (_erec_lep_nue); } else { - dunemcSamples[i].rw_erec = (_erec); - dunemcSamples[i].rw_erec_shifted = (_erec); - dunemcSamples[i].rw_erec_had = (_erec_had); - dunemcSamples[i].rw_erec_lep = (_erec_lep); + dunemcSamples[i].rw_erec = (_erec); + dunemcSamples[i].rw_erec_shifted = (_erec); + dunemcSamples[i].rw_erec_had = (_erec_had); + dunemcSamples[i].rw_erec_lep = (_erec_lep); } - - dunemcSamples[i].rw_eRecoP = (_eRecoP); - dunemcSamples[i].rw_eRecoPip = (_eRecoPip); - dunemcSamples[i].rw_eRecoPim = (_eRecoPim); - dunemcSamples[i].rw_eRecoPi0 = (_eRecoPi0); - dunemcSamples[i].rw_eRecoN = (_eRecoN); - dunemcSamples[i].rw_LepE = (_LepE); - dunemcSamples[i].rw_eP = (_eP); - dunemcSamples[i].rw_ePip = (_ePip); - dunemcSamples[i].rw_ePim = (_ePim); - dunemcSamples[i].rw_ePi0 = (_ePi0); + + dunemcSamples[i].rw_eRecoP = (_eRecoP); + dunemcSamples[i].rw_eRecoPip = (_eRecoPip); + dunemcSamples[i].rw_eRecoPim = (_eRecoPim); + dunemcSamples[i].rw_eRecoPi0 = (_eRecoPi0); + dunemcSamples[i].rw_eRecoN = (_eRecoN); + dunemcSamples[i].rw_LepE = (_LepE); + dunemcSamples[i].rw_eP = (_eP); + dunemcSamples[i].rw_ePip = (_ePip); + dunemcSamples[i].rw_ePim = (_ePim); + dunemcSamples[i].rw_ePi0 = (_ePi0); dunemcSamples[i].rw_eN = (_eN); // HH: Add checks to make sure the energies are not negative @@ -572,7 +536,7 @@ int SampleHandlerBeamFD::SetupExperimentMC() { if (dunemcSamples[i].rw_eRecoPi0 < 0) { dunemcSamples[i].rw_eRecoPi0 = 0; negative_counts["rw_eRecoPi0"]++; - } + } dunemcSamples[i].rw_erec_had_sqrt = sqrt(dunemcSamples[i].rw_erec_had); dunemcSamples[i].rw_erec_lep_sqrt = sqrt(dunemcSamples[i].rw_erec_lep); @@ -597,15 +561,15 @@ int SampleHandlerBeamFD::SetupExperimentMC() { dunemcSamples[i].rw_trueccnumu = static_cast(dunemcSamples[i].rw_isCC==1 && abs(dunemcSamples[i].nupdg)==14); dunemcSamples[i].rw_trueccnue = static_cast(dunemcSamples[i].rw_isCC==1 && abs(dunemcSamples[i].nupdg)==12); - + //Assume everything is on Argon40 for now.... dunemcSamples[i].Target = kTarget_Ar; - + int M3Mode = Modes->GetModeFromGenerator(std::abs(_mode)); if (!_isCC) M3Mode += 14; //Account for no ability to distinguish CC/NC if (M3Mode > 15) M3Mode -= 1; //Account for no NCSingleKaon dunemcSamples[i].mode = M3Mode; - + dunemcSamples[i].flux_w = 1.0; } @@ -615,7 +579,7 @@ int SampleHandlerBeamFD::SetupExperimentMC() { MACH3LOG_WARN("Found {} negative values for {}.", pair.second, pair.first); } } - + delete _data; return static_cast(nDownsampledEntries); } @@ -623,7 +587,7 @@ int SampleHandlerBeamFD::SetupExperimentMC() { const double* SampleHandlerBeamFD::GetPointerToKinematicParameter(const int KinPar, const int iEvent) const { switch(KinPar){ case kTrueNeutrinoEnergy: - return &(dunemcSamples[iEvent].enu_true); + return &(dunemcSamples[iEvent].enu_true); case kRecoNeutrinoEnergy: return &(dunemcSamples[iEvent].rw_erec_shifted); break; @@ -643,16 +607,16 @@ const double* SampleHandlerBeamFD::GetPointerToKinematicParameter(const int KinP return &(dunemcSamples[iEvent].OscChannelIndex); case kIsFHC: return &(beamFDSampleDetails[MCEvents[iEvent].NominalSample].isFHC); - case kTrueCCnue: + case kTrueCCnue: return &(dunemcSamples[iEvent].rw_trueccnue); - case kTrueCCnumu: + case kTrueCCnumu: return &(dunemcSamples[iEvent].rw_trueccnumu); case kTargetNucleus: return &(dunemcSamples[iEvent].Target); default: MACH3LOG_ERROR("Did not recognise Kinematic Parameter type {}...", KinPar); throw MaCh3Exception(__FILE__, __LINE__); - } + } } @@ -663,8 +627,8 @@ double SampleHandlerBeamFD::ReturnKinematicParameter(const int KinematicVariable void SampleHandlerBeamFD::SetupMC() { // dunemc_base *duneobj = &(dunemcSamples[iSample]); - // FarDetectorCoreInfo *fdobj = &(MCEvents[iSample]); - + // FarDetectorCoreInfo *fdobj = &(MCEvents[iSample]); + for (unsigned int iEvent = 0; iEvent < GetNEvents(); ++iEvent) { MCEvents[iEvent].enu_true = dunemcSamples[iEvent].enu_true; MCEvents[iEvent].isNC = !(dunemcSamples[iEvent].rw_isCC); @@ -672,6 +636,6 @@ void SampleHandlerBeamFD::SetupMC() { MCEvents[iEvent].nupdgUnosc = dunemcSamples[iEvent].nupdgUnosc; MCEvents[iEvent].NominalSample = dunemcSamples[iEvent].SampleIndex; } - + } - + diff --git a/Samples/SampleHandlerBeamFD.h b/Samples/SampleHandlerBeamFD.h index a50d731e..1bbfbc79 100644 --- a/Samples/SampleHandlerBeamFD.h +++ b/Samples/SampleHandlerBeamFD.h @@ -56,15 +56,15 @@ class SampleHandlerBeamFD : virtual public SampleHandlerBase /// @brief Sets up pointers weights for each event (oscillation/xsec/etc.) void AddAdditionalWeightPointers(); void SetupSplines(); - + // === HH: Functional parameters === - enum FuncParEnum {kTotalEScale, kTotalEScaleNotCCNumu, - kTotalEScaleSqrt, kTotalEScaleSqrtNotCCNumu, + enum FuncParEnum {kTotalEScale, kTotalEScaleNotCCNumu, + kTotalEScaleSqrt, kTotalEScaleSqrtNotCCNumu, kTotalEScaleInvSqrt, kTotalEScaleInvSqrtNotCCNumu, kHadEScale, kHadEScaleSqrt, kHadEScaleInvSqrt, kMuEScale, kMuEScaleSqrt, kMuEScaleInvSqrt, kNEScale, kNEScaleSqrt, kNEScaleInvSqrt, - kEMEScale, kEMEScaleCCNue, + kEMEScale, kEMEScaleCCNue, kEMEScaleSqrt, kEMEScaleSqrtCCNue, kEMEScaleInvSqrt, kEMEScaleInvSqrtCCNue, kHadRes, kMuRes, kNRes, kEMRes, kEMResCCNue, @@ -73,57 +73,17 @@ class SampleHandlerBeamFD : virtual public SampleHandlerBase void RegisterFunctionalParameters() override; void ResetShifts(int iEvent) override; - // Global energy scale systematics - void TotalEScale(const double * par, std::size_t iEvent); - void TotalEScaleNotCCNumu(const double * par, std::size_t iEvent); - void TotalEScaleSqrt(const double * par, std::size_t iEvent); - void TotalEScaleSqrtNotCCNumu(const double * par, std::size_t iEvent); - void TotalEScaleInvSqrt(const double * par, std::size_t iEvent); - void TotalEScaleInvSqrtNotCCNumu(const double * par, std::size_t iEvent); - - // Particle specific energy uncertainties - // Charged hadron - void HadEScale(const double * par, std::size_t iEvent); - void HadEScaleSqrt(const double * par, std::size_t iEvent); - void HadEScaleInvSqrt(const double * par, std::size_t iEvent); - // Muons - void MuEScale(const double * par, std::size_t iEvent); - void MuEScaleSqrt(const double * par, std::size_t iEvent); - void MuEScaleInvSqrt(const double * par, std::size_t iEvent); - // Neutrons - void NEScale(const double * par, std::size_t iEvent); - void NEScaleSqrt(const double * par, std::size_t iEvent); - void NEScaleInvSqrt(const double * par, std::size_t iEvent); - // Electromagnetic showers - void EMEScale(const double * par, std::size_t iEvent); - void EMEScaleCCNue(const double * par, std::size_t iEvent); - void EMEScaleSqrt(const double * par, std::size_t iEvent); - void EMEScaleSqrtCCNue(const double * par, std::size_t iEvent); - void EMEScaleInvSqrt(const double * par, std::size_t iEvent); - void EMEScaleInvSqrtCCNue(const double * par, std::size_t iEvent); - - // Resolution uncertainties - void HadRes(const double * par, std::size_t iEvent); - void MuRes(const double * par, std::size_t iEvent); - void NRes(const double * par, std::size_t iEvent); - void EMRes(const double * par, std::size_t iEvent); - void EMResCCNue(const double * par, std::size_t iEvent); - - //Reconstruction (CVN) uncertainties - void RecoCVNNumu(const double * par, std::size_t iEvent); - void RecoCVNNue(const double * par, std::size_t iEvent); - /// @brief Returns pointer to kinemtatic parameter for event in Structs DUNE /// @param KinematicVariable Kinematic parameter ID as int /// @param iEvent Event ID - /// @return Value of kinematic parameter corresponding for a given event + /// @return Value of kinematic parameter corresponding for a given event double ReturnKinematicParameter (const int KinematicVariable, const int iEvent) const override; /// @brief Returns pointer to kinemtatic parameter for event in Structs DUNE /// @param KinematicVariable Kinematic parameter as double (gets cast -> int) /// @param iEvent Event ID /// @return Pointer to KinPar for a given event - const double* GetPointerToKinematicParameter(const int KinematicVariable, const int iEvent) const override; + const double* GetPointerToKinematicParameter(const int KinematicVariable, const int iEvent) const override; //DB functions which could be initialised to do something which is non-trivial /// @brief NOT IMPLEMENTED: Dunder method to calculate xsec weights @@ -131,8 +91,8 @@ class SampleHandlerBeamFD : virtual public SampleHandlerBase double CalcXsecWeightFunc(int iEvent) {(void)iEvent; return 1.;} // dunemc - /// DUNE MC sampels - std::vector dunemcSamples; + /// DUNE MC samples + std::vector dunemcSamples; std::vector beamFDSampleDetails; const std::unordered_map KinematicParametersDUNE = { @@ -174,7 +134,7 @@ class SampleHandlerBeamFD : virtual public SampleHandlerBase /// @brief Cleanup memory void CleanMemoryBeforeFit() override {}; }; - + #endif diff --git a/Samples/SampleHandlerBeamND.cpp b/Samples/SampleHandlerBeamND.cpp index 9f1ae71f..58fa490c 100644 --- a/Samples/SampleHandlerBeamND.cpp +++ b/Samples/SampleHandlerBeamND.cpp @@ -7,12 +7,12 @@ SampleHandlerBeamND::SampleHandlerBeamND(std::string mc_version_, ParameterHandl MACH3LOG_ERROR("You've passed me a nullptr to a ND covarince matrix... "); throw MaCh3Exception(__FILE__, __LINE__); } - + beamNDCov = beamNDCov_; KinematicParameters = &KinematicParametersDUNE; ReversedKinematicParameters = &ReversedKinematicParametersDUNE; - + Initialise(); downsamplingStep = 1; @@ -23,7 +23,7 @@ SampleHandlerBeamND::~SampleHandlerBeamND() { void SampleHandlerBeamND::Init() { beamNDSampleDetails.resize(GetNSamples()); - + auto EnabledSamples = Get>(SampleManager->raw()["Samples"], __FILE__ , __LINE__); for (int i = 0; i < GetNSamples(); i++){ @@ -32,7 +32,7 @@ void SampleHandlerBeamND::Init() { beamNDSampleDetails[i].iselike = SampleManager->raw()[TempTitle]["DUNESampleBools"]["iselike"].as(); beamNDSampleDetails[i].pot = SampleManager->raw()[TempTitle]["POT"].as(); - if (beamNDSampleDetails[i].isFHC) { + if (beamNDSampleDetails[i].isFHC) { beamNDSampleDetails[i].norm_s = (1e21/1.905e21); } else { beamNDSampleDetails[i].norm_s = (1e21/1.5e21); @@ -44,12 +44,12 @@ void SampleHandlerBeamND::Init() { MACH3LOG_INFO("- iselike: {}", beamNDSampleDetails[i].iselike); } - downsamplingStep = GetFromManager(SampleManager->raw()["DownsamplingStep"], 1); + downsamplingStep = GetFromManager(SampleManager->raw()["DownsamplingStep"], 1, __FILE__ , __LINE__); if (downsamplingStep == 0) { throw MaCh3Exception(__FILE__, __LINE__, "Downsampling step cannot be zero. Please set it to a positive integer in the Beam ND sample config file." ); - } + } MACH3LOG_INFO("Beam ND downsampling step: {}", downsamplingStep); MACH3LOG_INFO("-------------------------------------------------------------------"); @@ -95,7 +95,7 @@ void SampleHandlerBeamND::AddAdditionalWeightPointers() { int SampleHandlerBeamND::SetupExperimentMC() { // dunemc_base *duneobj = &(dunendmcSamples[iSample]); - + MACH3LOG_INFO("-------------------------------------------------------------------"); TChain* _data = new TChain("caf"); @@ -183,7 +183,7 @@ int SampleHandlerBeamND::SetupExperimentMC() { // POT stuff dunendmcSamples[i].norm_s = beamNDSampleDetails[sample_index].norm_s; dunendmcSamples[i].pot_s = beamNDSampleDetails[sample_index].pot_s * downsamplingStep; - + dunendmcSamples[i].rw_erec = _erec; dunendmcSamples[i].rw_erec_shifted = _erec; dunendmcSamples[i].rw_erec_lep = _erec_lep; @@ -194,7 +194,7 @@ int SampleHandlerBeamND::SetupExperimentMC() { dunendmcSamples[i].rw_isCC = _isCC; dunendmcSamples[i].rw_reco_q = _reco_q; dunendmcSamples[i].rw_berpaacvwgt = _BeRPA_cvwgt; - + //Assume everything is on Argon for now.... dunendmcSamples[i].Target = kTarget_Ar; @@ -202,19 +202,19 @@ int SampleHandlerBeamND::SetupExperimentMC() { if (!_isCC) M3Mode += 14; //Account for no ability to distinguish CC/NC if (M3Mode > 15) M3Mode -= 1; //Account for no NCSingleKaon dunendmcSamples[i].mode = M3Mode; - + dunendmcSamples[i].flux_w = 1.0; } - + //_sampleFile->Close(); _data->Reset(); delete _data; return static_cast(nDownsampledEntries); - + } -const double* SampleHandlerBeamND::GetPointerToKinematicParameter(const int KinPar, const int iEvent) const{ +const double* SampleHandlerBeamND::GetPointerToKinematicParameter(const int KinPar, const int iEvent) const{ switch(KinPar){ case kTrueNeutrinoEnergy: return &(dunendmcSamples[iEvent].enu_true); @@ -235,7 +235,7 @@ const double* SampleHandlerBeamND::GetPointerToKinematicParameter(const int KinP default: MACH3LOG_ERROR("Did not recognise Kinematic Parameter type..."); throw MaCh3Exception(__FILE__, __LINE__); - } + } } @@ -248,7 +248,7 @@ double SampleHandlerBeamND::ReturnKinematicParameter(const int KinematicVariable void SampleHandlerBeamND::SetupMC() { // dunemc_base *duneobj = &(dunendmcSamples[iSample]); // FarDetectorCoreInfo *fdobj = &(MCEvents[iSample]); - + for (unsigned int iEvent = 0; iEvent < GetNEvents(); ++iEvent) { MCEvents[iEvent].enu_true = dunendmcSamples[iEvent].enu_true; MCEvents[iEvent].isNC = !(dunendmcSamples[iEvent].rw_isCC); diff --git a/Splines/SplineHandlerFactoryDUNE.cpp b/Splines/SplineHandlerFactoryDUNE.cpp index 012965b2..3870bb56 100644 --- a/Splines/SplineHandlerFactoryDUNE.cpp +++ b/Splines/SplineHandlerFactoryDUNE.cpp @@ -78,7 +78,7 @@ SplineHandlerFactoryDUNE::SplineHandlerFactoryDUNE(ParameterHandlerGeneric* xsec for(uint i = 0; i < samplePars.size(); ++i) { splineWeightPtrs[i] = xsec_params->RetPointer(samplePars[i].index); } - monolith->setSplinePointers(splineWeightPtrs); + monolith->SetSplinePointers(splineWeightPtrs); fSplineHandler = std::move(monolith); break; diff --git a/cmake/Modules/fetchduneanaobj.cmake b/cmake/Modules/fetchduneanaobj.cmake new file mode 100644 index 00000000..4650aa65 --- /dev/null +++ b/cmake/Modules/fetchduneanaobj.cmake @@ -0,0 +1,67 @@ +macro(fetchduneanaobj DUNE_ANAOBJ_BRANCH) + + include(FetchContent) + FetchContent_Declare( + duneanaobj + GIT_REPOSITORY https://github.com/DUNE/duneanaobj.git + GIT_TAG ${DUNE_ANAOBJ_BRANCH} + ) + + + # Check if population has already been performed + FetchContent_GetProperties(duneanaobj) + if(NOT duneanaobj_POPULATED) + # Fetch the content using previously declared details + FetchContent_Populate(duneanaobj) + + endif() + + include_directories(${duneanaobj_SOURCE_DIR}) + + ROOT_GENERATE_DICTIONARY(StandardRecordDict + ${duneanaobj_SOURCE_DIR}/duneanaobj/StandardRecord/classes.h + LINKDEF + ${duneanaobj_SOURCE_DIR}/duneanaobj/StandardRecord/classes_def.xml) + + add_library(libStandardRecordDict SHARED ${CMAKE_CURRENT_BINARY_DIR}/StandardRecordDict.cxx) + target_link_libraries(libStandardRecordDict PUBLIC ROOT::RIO) + set_target_properties(libStandardRecordDict PROPERTIES OUTPUT_NAME "StandardRecordDict") + + target_include_directories(libStandardRecordDict PUBLIC + $ + $) + target_include_directories(libStandardRecordDict PRIVATE + $ #root puts this in the dictionary + ) + + file(GLOB SR_IMPL_FILES "${duneanaobj_SOURCE_DIR}/duneanaobj/StandardRecord/*.cxx") + add_library(duneanaobj_StandardRecord SHARED ${SR_IMPL_FILES}) + target_link_libraries(duneanaobj_StandardRecord PUBLIC libStandardRecordDict ROOT::MathCore ROOT::Physics) + + target_include_directories(duneanaobj_StandardRecord PUBLIC + $ + $) + + set_target_properties(duneanaobj_StandardRecord PROPERTIES EXPORT_NAME StandardRecord) + + file(GLOB SR_HEADER_FILES "${duneanaobj_SOURCE_DIR}/duneanaobj/StandardRecord/*.h") + + install(FILES ${SR_HEADER_FILES} DESTINATION include/duneanaobj/StandardRecord) + install(FILES + ${CMAKE_CURRENT_BINARY_DIR}/libStandardRecordDict_rdict.pcm + ${CMAKE_CURRENT_BINARY_DIR}/libStandardRecordDict.rootmap + DESTINATION lib) + + install(TARGETS duneanaobj_StandardRecord libStandardRecordDict + EXPORT duneanaobj-targets + LIBRARY DESTINATION lib) + + install(EXPORT duneanaobj-targets + FILE duneanaobjTargets.cmake + NAMESPACE duneanaobj:: + DESTINATION ${CMAKE_INSTALL_PREFIX}/ + ) + + add_library(duneanaobj::StandardRecord ALIAS duneanaobj_StandardRecord) + +endmacro() From f6295ad374c9736eadd1f88415bdd35546579f56 Mon Sep 17 00:00:00 2001 From: Luke Pickering Date: Mon, 10 Aug 2026 11:32:17 +0100 Subject: [PATCH 2/2] Update SampleHandlerBeamFD.h --- Samples/SampleHandlerBeamFD.h | 13 ------------- 1 file changed, 13 deletions(-) diff --git a/Samples/SampleHandlerBeamFD.h b/Samples/SampleHandlerBeamFD.h index 1bbfbc79..e11aaa96 100644 --- a/Samples/SampleHandlerBeamFD.h +++ b/Samples/SampleHandlerBeamFD.h @@ -57,19 +57,6 @@ class SampleHandlerBeamFD : virtual public SampleHandlerBase void AddAdditionalWeightPointers(); void SetupSplines(); - // === HH: Functional parameters === - enum FuncParEnum {kTotalEScale, kTotalEScaleNotCCNumu, - kTotalEScaleSqrt, kTotalEScaleSqrtNotCCNumu, - kTotalEScaleInvSqrt, kTotalEScaleInvSqrtNotCCNumu, - kHadEScale, kHadEScaleSqrt, kHadEScaleInvSqrt, - kMuEScale, kMuEScaleSqrt, kMuEScaleInvSqrt, - kNEScale, kNEScaleSqrt, kNEScaleInvSqrt, - kEMEScale, kEMEScaleCCNue, - kEMEScaleSqrt, kEMEScaleSqrtCCNue, - kEMEScaleInvSqrt, kEMEScaleInvSqrtCCNue, - kHadRes, kMuRes, kNRes, kEMRes, kEMResCCNue, - kRecoCVNNumu, kRecoCVNNue - }; void RegisterFunctionalParameters() override; void ResetShifts(int iEvent) override;