diff --git a/Apps/CMakeLists.txt b/Apps/CMakeLists.txt index 6cf68401..ce26e4fa 100644 --- a/Apps/CMakeLists.txt +++ b/Apps/CMakeLists.txt @@ -8,9 +8,9 @@ foreach(app LikelihoodScan Variations ) - + add_executable(${app} ${app}.cpp) - target_link_libraries(${app} SamplesDUNE MaCh3DUNE::All Minuit) + target_link_libraries(${app} SamplesDUNE MaCh3DUNE::All ROOT::Minuit) add_dependencies(MaCh3DUNEApps ${app}) diff --git a/Apps/SigmaVariation.cpp b/Apps/SigmaVariation.cpp index 45650c9c..693c3c93 100644 --- a/Apps/SigmaVariation.cpp +++ b/Apps/SigmaVariation.cpp @@ -12,84 +12,36 @@ #include #include -#include "Samples/MaCh3DUNEFactory.h" -#include "Samples/StructsDUNE.h" #include "Fitters/MaCh3Factory.h" +#include "Samples/MaCh3DUNEFactory.h" int main(int argc, char * argv[]) { - auto FitManager = MaCh3ManagerFactory(argc, argv); - //############################################################################################################################### - - //DB Sigma variations in units of each parameters Sigma - std::vector sigmaVariations = {-3, -1, 0, 1, 3}; + auto FitManager = MaCh3ManagerFactory(argc, argv); + auto OutputFileName = FitManager->raw()["General"]["OutputFile"].as(); - //############################################################################################################################### + //#################################################################################### //Create sample handler + parameter_handler objects auto [param_handler, samples] = MaCh3DuneFactory(FitManager); - //############################################################################################################################### - //Perform reweight and print total integral - - MACH3LOG_INFO("======================================================="); - for(SampleHandlerBase* handler: samples){ - handler->Reweight(); - for (int iSample=0;iSampleGetNSamples();iSample++) { - MACH3LOG_INFO("Event rate for {} : {:<5.2f}", handler->GetSampleTitle(iSample), handler->GetMCHist(iSample)->Integral()); - } - } - MACH3LOG_INFO("======================================================="); - - 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()); + //########################################################################################################### + //MCMC - int nPars = param_handler->GetNumParams(); - for (int iPar=0;iParGetParFancyName(iPar); - double VarInit = param_handler->GetParPreFit(iPar); - double VarSigma = param_handler->GetDiagonalError(iPar); + auto MaCh3Fitter = std::make_unique(FitManager.get()); - MACH3LOG_INFO("\tParameter : {:<30} - Variations around value : {:<10.7f} , in units of 1 Sigma : {:<10.7f}",ParName,VarInit,VarSigma); + //Add systematic objects + MaCh3Fitter->AddSystObj(param_handler.get()); - File->cd(); - File->mkdir(ParName.c_str()); - File->cd(ParName.c_str()); - - for (size_t iSigVar=0;iSigVarGetLowerBound(iPar)) VarVal = param_handler->GetLowerBound(iPar); - if (VarVal > param_handler->GetUpperBound(iPar)) VarVal = param_handler->GetUpperBound(iPar); - - MACH3LOG_INFO("\t\tVariation {:<5.3f} - Parameter Value : {:<10.7f}", - sigmaVariations[iSigVar], VarVal); - param_handler->SetParProp(iPar, VarVal); - - for (auto handler : samples) { - for (int iSample = 0; iSample < handler->GetNSamples(); iSample++) { - std::string SampleName = handler->GetSampleTitle(iSample); - - File->cd(ParName.c_str()); - if (iSigVar == 0) { - File->mkdir((ParName + "/" + SampleName).c_str()); - } - File->cd((ParName + "/" + SampleName).c_str()); - - handler->Reweight(); - auto Hist = handler->GetMCHist(iSample); - MACH3LOG_INFO("\t\t\tSample : {:<30} - Integral : {:<10}", SampleName, - Hist->Integral()); - - Hist->Write(Form("Variation_%i", (int)iSigVar)); - } - } - } + //Add samples + for(auto Sample : samples){ + MaCh3Fitter->AddSampleHandler(Sample); + } - param_handler->SetParProp(iPar, VarInit); + MaCh3Fitter->RunSigmaVar(); - MACH3LOG_INFO("======================================================="); - } + //Writing the memory usage at the end to eventually spot some nasty leak + MACH3LOG_WARN("\033[0;31mCurrent Total RAM usage is {:.2f} GB\033[0m", M3::Utils::getValue("VmRSS") / 1048576.0); + MACH3LOG_WARN("\033[0;31mOut of Total available RAM {:.2f} GB\033[0m", M3::Utils::getValue("MemTotal") / 1048576.0); return 0; } diff --git a/CMakeLists.txt b/CMakeLists.txt index b8710a2e..e6b82e26 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -39,24 +39,8 @@ file( ) include(${CMAKE_CURRENT_BINARY_DIR}/cmake/CPM.cmake) -CPMFindPackage( - NAME CMakeModules - GIT_TAG stable - GITHUB_REPOSITORY NuHepMC/CMakeModules - DOWNLOAD_ONLY -) -include(${CMakeModules_SOURCE_DIR}/NuHepMCModules.cmake) - #=================== ROOT -include(ROOT) - -if(NOT TARGET ROOT::ROOT) - cmessage(FATAL_ERROR "MaCh3 Expected dependency target: ROOT::ROOT") -endif() - -if(DEFINED ROOT_CXX_STANDARD AND ROOT_CXX_STANDARD GREATER CMAKE_CXX_STANDARD) - set(CMAKE_CXX_STANDARD ${ROOT_CXX_STANDARD}) -endif() +find_package(ROOT 6 REQUIRED) #=================== DUNEAnaObj if(NOT DEFINED DUNE_ANAOBJ_BRANCH) @@ -79,22 +63,44 @@ if(MaCh3_DUNE_USE_SRProxy) endif() if(NOT TARGET duneanaobj::all) - cmessage(FATAL_ERROR "MaCh3 DUNE Expected dependency target: duneanaobj::all") + message(FATAL_ERROR "MaCh3 DUNE Expected dependency target: duneanaobj::all") endif() else() include(fetchduneanaobj) fetchduneanaobj(${DUNE_ANAOBJ_BRANCH}) endif() +include(FetchContent) +FetchContent_Declare( + duneanafluxtools + GIT_REPOSITORY https://github.com/DUNE/duneanafluxtools.git + GIT_TAG v0.9.0 +) +# we don't need to build the caf interface +SET(duneanafluxtools_USE_DUNEANAOBJ OFF) +FetchContent_MakeAvailable(duneanafluxtools) + +CPMAddPackage( + NAME duneanasel + GITHUB_REPOSITORY DUNE/duneanasel + GIT_TAG main + DOWNLOAD_ONLY YES +) +add_library(duneanasel_all INTERFACE) +target_include_directories(duneanasel_all INTERFACE ${duneanasel_SOURCE_DIR}) +add_library(duneanasel::all ALIAS duneanasel_all) + ################################## MaCh3 ###################################### #If MaCh3 was sourced find it, otherwise use CPM SET(MaCh3_FOUND FALSE) -find_package(MaCh3 2.6.0 EXACT QUIET) +# don't find mach3 core from a previous build in this directory or sad times occur +LIST(APPEND CMAKE_IGNORE_PATH ${CMAKE_CURRENT_BINARY_DIR}) +find_package(MaCh3 2.6.0 CONFIG 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") + message(STATUS "Didn't find MaCh3, attempting to use built in MaCh3") if(NOT DEFINED MaCh3_GPU_ENABLED) set(MaCh3_GPU_ENABLED ON) @@ -143,7 +149,7 @@ else() endif() if(NOT TARGET MaCh3::All) - cmessage(FATAL_ERROR "MaCh3 DUNE Expected dependency target: MaCh3::All, did you source setup.MaCh3.sh") + message(FATAL_ERROR "MaCh3 DUNE Expected dependency target: MaCh3::All, did you source setup.MaCh3.sh") endif() ############################ C++ Compiler #################################### @@ -153,7 +159,7 @@ if(NOT DEFINED MACH3_CXX_STANDARD) else() SET(CMAKE_CXX_STANDARD ${MACH3_CXX_STANDARD}) endif() -cmessage(STATUS "CMAKE CXX Standard: ${CMAKE_CXX_STANDARD}") +message(STATUS "CMAKE CXX Standard: ${CMAKE_CXX_STANDARD}") ############################ FLAGS #################################### add_library(DUNEMaCh3Warnings INTERFACE) @@ -185,9 +191,9 @@ target_compile_options(DUNEMaCh3Warnings INTERFACE LIST(APPEND ALL_FEATURES ) -cmessage(STATUS "MaCh3DUNE Features: ") +message(STATUS "MaCh3DUNE Features: ") foreach(f ${ALL_FEATURES}) - cmessage(STATUS " ${f}: ${MaCh3DUNE_${f}_ENABLED}") + message(STATUS " ${f}: ${MaCh3DUNE_${f}_ENABLED}") endforeach() #KS: Options to print dependency graph @@ -217,13 +223,13 @@ set_target_properties(MaCh3DUNECompilerOptions PROPERTIES EXPORT_NAME CompilerOp target_include_directories(MaCh3DUNECompilerOptions INTERFACE $ - $ $) install(TARGETS MaCh3DUNECompilerOptions EXPORT mach3dune-targets LIBRARY DESTINATION lib/) +add_subdirectory(Systematics) add_subdirectory(Splines) add_subdirectory(Samples) add_subdirectory(Apps) @@ -237,41 +243,41 @@ if(NOT DEFINED MaCh3_PREFIX OR MaCh3_PREFIX STREQUAL "") # Try to detect if MaCh3 was fetched by CPM if(EXISTS "${CMAKE_BINARY_DIR}/_deps/mach3-src") set(MaCh3_PREFIX "${CMAKE_BINARY_DIR}/_deps/mach3-src") - cmessage(STATUS "Detected CPM-fetched MaCh3: ${MaCh3_PREFIX}") + message(STATUS "Detected CPM-fetched MaCh3: ${MaCh3_PREFIX}") # Try to get from MaCh3 package config elseif(DEFINED MaCh3_DIR) get_filename_component(MaCh3_PREFIX "${MaCh3_DIR}" DIRECTORY) get_filename_component(MaCh3_PREFIX "${MaCh3_PREFIX}" DIRECTORY) - cmessage(STATUS "Detected installed MaCh3 from MaCh3_DIR: ${MaCh3_PREFIX}") + message(STATUS "Detected installed MaCh3 from MaCh3_DIR: ${MaCh3_PREFIX}") # Try to get from target properties elseif(TARGET MaCh3::All) get_target_property(MACH3_INCLUDE_DIRS MaCh3::All INTERFACE_INCLUDE_DIRECTORIES) if(MACH3_INCLUDE_DIRS) list(GET MACH3_INCLUDE_DIRS 0 FIRST_INCLUDE_DIR) get_filename_component(MaCh3_PREFIX "${FIRST_INCLUDE_DIR}" DIRECTORY) - cmessage(STATUS "Detected MaCh3 from target include dirs: ${MaCh3_PREFIX}") + message(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=") + message(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}") + message(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}") + message(STATUS "Set MaCh3_PREFIX: ${MaCh3_PREFIX}") endif() if(PYTHON_ENABLED) - cmessage(STATUS "PYTHON IS ENABLED!!") + message(STATUS "PYTHON IS ENABLED!!") set(CMAKE_INSTALL_RPATH "$ORIGIN/../lib") - cmessage(STATUS "CMAKE_INSTALL_RPATH is ${CMAKE_INSTALL_RPATH}") + message(STATUS "CMAKE_INSTALL_RPATH is ${CMAKE_INSTALL_RPATH}") set(PYMACH3_PATH ${CMAKE_BINARY_DIR}/python/pyMaCh3) add_subdirectory(python) endif() @@ -291,6 +297,7 @@ file(MAKE_DIRECTORY ${CMAKE_BINARY_DIR}/output) install(DIRECTORY Inputs DESTINATION ${CMAKE_BINARY_DIR}) install(DIRECTORY Configs DESTINATION ${CMAKE_BINARY_DIR}) install(DIRECTORY Utils DESTINATION ${CMAKE_BINARY_DIR}) +install(DIRECTORY plotting DESTINATION ${CMAKE_BINARY_DIR}) install(EXPORT mach3dune-targets FILE MaCh3DUNETargets.cmake diff --git a/Configs/CovObjs/BeamParameters2026.yaml b/Configs/CovObjs/BeamParameters2026.yaml new file mode 100644 index 00000000..28a68cb3 --- /dev/null +++ b/Configs/CovObjs/BeamParameters2026.yaml @@ -0,0 +1,994 @@ +--- +Systematics: + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: TargetUpstreamDegredation + ParameterName: TargetUpstreamDegredation + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: TargetTiltTransverseY + ParameterName: TargetTiltTransverseY + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: TargetTiltTransverseX + ParameterName: TargetTiltTransverseX + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: TargetLength + ParameterName: TargetLength + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: TargetDisplaceTransverseY + ParameterName: TargetDisplaceTransverseY + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: TargetDisplaceTransverseX + ParameterName: TargetDisplaceTransverseX + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: TargetDensity + ParameterName: TargetDensity + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: ProtonBeamTransverseY + ParameterName: ProtonBeamTransverseY + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: ProtonBeamTransverseX + ParameterName: ProtonBeamTransverseX + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: ProtonBeamRadius + ParameterName: ProtonBeamRadius + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: ProtonBeamAngleY + ParameterName: ProtonBeamAngleY + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: ProtonBeamAngleX + ParameterName: ProtonBeamAngleX + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: HornWaterLayerThickness + ParameterName: HornWaterLayerThickness + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: HornCurrent + ParameterName: HornCurrent + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: HornCTiltTransverseY + ParameterName: HornCTiltTransverseY + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: HornCTiltTransverseX + ParameterName: HornCTiltTransverseX + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: HornCEllipticityXInducedBField + ParameterName: HornCEllipticityXInducedBField + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: HornCEccentricityXInducedBField + ParameterName: HornCEccentricityXInducedBField + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: HornCDisplaceLongitudinalZ + ParameterName: HornCDisplaceLongitudinalZ + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: HornBTiltTransverseY + ParameterName: HornBTiltTransverseY + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: HornBTiltTransverseX + ParameterName: HornBTiltTransverseX + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: HornBEllipticityXInducedBField + ParameterName: HornBEllipticityXInducedBField + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: HornBDisplaceLongitudinalZ + ParameterName: HornBDisplaceLongitudinalZ + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: HornATiltTransverseY + ParameterName: HornATiltTransverseY + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: HornATiltTransverseX + ParameterName: HornATiltTransverseX + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: HornAEllipticityXInducedBField + ParameterName: HornAEllipticityXInducedBField + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: HornAEccentricityXInducedBField + ParameterName: HornAEccentricityXInducedBField + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: HornADisplaceLongitudinalZ + ParameterName: HornADisplaceLongitudinalZ + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: DecayPipeTiltY + ParameterName: DecayPipeTiltY + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: DecayPipeTiltX + ParameterName: DecayPipeTiltX + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: DecayPipeRadius + ParameterName: DecayPipeRadius + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: DecayPipeLength + ParameterName: DecayPipeLength + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: DecayPipeGeoBField + ParameterName: DecayPipeGeoBField + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: DecayPipeEllipticalCrossSectionYB + ParameterName: DecayPipeEllipticalCrossSectionYB + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: DecayPipeEllipticalCrossSectionXA + ParameterName: DecayPipeEllipticalCrossSectionXA + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: DecayPipeDisplaceTransverseY + ParameterName: DecayPipeDisplaceTransverseY + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: DecayPipeDisplaceTransverseX + ParameterName: DecayPipeDisplaceTransverseX + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: DecayPipe3SegmentBowingY + ParameterName: DecayPipe3SegmentBowingY + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: DecayPipe3SegmentBowingX + ParameterName: DecayPipe3SegmentBowingX + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: HornCDisplaceTransverseY + ParameterName: HornCDisplaceTransverseY + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: HornBDisplaceTransverseY + ParameterName: HornBDisplaceTransverseY + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: HornADisplaceTransverseY + ParameterName: HornADisplaceTransverseY + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: HornCDisplaceTransverseX + ParameterName: HornCDisplaceTransverseX + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: HornBDisplaceTransverseX + ParameterName: HornBDisplaceTransverseX + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: HornADisplaceTransverseX + ParameterName: HornADisplaceTransverseX + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + + + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: Flux_HadProd_Param_0 + ParameterName: Flux_HadProd_Param_0 + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: Flux_HadProd_Param_1 + ParameterName: Flux_HadProd_Param_1 + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: Flux_HadProd_Param_2 + ParameterName: Flux_HadProd_Param_2 + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: Flux_HadProd_Param_3 + ParameterName: Flux_HadProd_Param_3 + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: Flux_HadProd_Param_4 + ParameterName: Flux_HadProd_Param_4 + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: Flux_HadProd_Param_5 + ParameterName: Flux_HadProd_Param_5 + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: Flux_HadProd_Param_6 + ParameterName: Flux_HadProd_Param_6 + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: Flux_HadProd_Param_7 + ParameterName: Flux_HadProd_Param_7 + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: Flux_HadProd_Param_8 + ParameterName: Flux_HadProd_Param_8 + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: Flux_HadProd_Param_9 + ParameterName: Flux_HadProd_Param_9 + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: Flux_HadProd_Param_10 + ParameterName: Flux_HadProd_Param_10 + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: Flux_HadProd_Param_11 + ParameterName: Flux_HadProd_Param_11 + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: Flux_HadProd_Param_12 + ParameterName: Flux_HadProd_Param_12 + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: Flux_HadProd_Param_13 + ParameterName: Flux_HadProd_Param_13 + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: Flux_HadProd_Param_14 + ParameterName: Flux_HadProd_Param_14 + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: Flux_HadProd_Param_15 + ParameterName: Flux_HadProd_Param_15 + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: Flux_HadProd_Param_16 + ParameterName: Flux_HadProd_Param_16 + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: Flux_HadProd_Param_17 + ParameterName: Flux_HadProd_Param_17 + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: Flux_HadProd_Param_18 + ParameterName: Flux_HadProd_Param_18 + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: Flux_HadProd_Param_19 + ParameterName: Flux_HadProd_Param_19 + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional + - Systematic: + SampleNames: ["*"] + Error: 1 + FlatPrior: false + Names: + FancyName: Flux_HadProd_Param_20 + ParameterName: Flux_HadProd_Param_20 + ParameterBounds: [ -3, 3 ] + ParameterGroup: Flux + ParameterValues: + Generated: 0.0 + PreFitValue: 0.0 + StepScale: + MCMC: 0.02 + Type: Functional diff --git a/Configs/EventRates_BeamFDStandardRecord.yaml b/Configs/EventRates_BeamFDStandardRecord.yaml new file mode 100644 index 00000000..e5035587 --- /dev/null +++ b/Configs/EventRates_BeamFDStandardRecord.yaml @@ -0,0 +1,63 @@ +--- +General: + OutputFile: "DuneEventRates.root" + DUNESamples: + - "Configs/Samples/SampleHandler_BeamFDStandardRecord.yaml" + # Nu-FIT + OscillationParameters: [0.310, 0.582, 0.0224, 7.39E-5, 2.525E-3, -2.498, 1284.9, 2.848] #CUDAProb3 + # OscillationParameters: [0.310, 0.582, 0.0224, 7.39E-5, 2.525E-3, -2.498, 1284.9, 2.848, 0.5] #NuFAST + + # T2K-like best-fit + #OscillationParameters: [0.307, 0.528, 0.0218, 7.53e-5, 2.509e-3, -1.601, 1284.9, 2.848] #CUDAProb3 + #OscillationParameters: [0.307, 0.528, 0.0218, 7.53e-5, 2.509e-3, -1.601, 1284.9, 2.848, 0.5] #NuFAST + + OscillatorConfigName: "Configs/OscillatorObj.yaml" + + Systematics: + XsecCovFile: + - "Configs/CovObjs/OscCov_PDG2021_v2.yaml" + - "Configs/CovObjs/BeamParameters2026.yaml" + XsecCovName: "xsec_cov" + XsecStepScale: 0.05 + + + Fitter: + FitTestLikelihood: false + MCMC: + NSteps: 2000 + AutoSave: 10000 + Output: + FileName: "TestEventRates.root" + OUTPUTNAME: "TestLLH.root" + ProcessMCMC: No + Seed: 0 + Debug: No + +SigmaVar: + # whether to separate by mach3 mode + PlotByMode: false + # whether to seaparate by osc channel + PlotByChannel: false + # list of parameters to skip, can include glob patterns e.g: MyParamPrefix_* + # SkipVector: [sin2th_12, sin2th_23, sin2th_13, delm2_12, delm2_23, delta_cp] + # parameter values + SigmaArray: [-1,-0.5,0,0.5,1] + +# LLH Scan settings +LLHScan: + # don't make LLH for certain systematic groups + # LLHScanSkipVector: ["b_"] + # Save individual sample contributions to the LLH + LLHScanBySample: false + # Number of points we run LLH scans for each systematic + LLHScanPoints: 10 + # Number of points for 2D scan in each direction + #2DLLHScanPoints: 20 + # Add physical ranges for special params whose priors/errors may not be well well-defined + ScanRanges: + #delta_cp: [-3.14,3.14] + +# Tests +Tests: + SkipChecks: False + TestResultsFile: tests/test_results/EventRatesV2TDR.txt diff --git a/Configs/Samples/SampleHandler_BeamFDStandardRecord.yaml b/Configs/Samples/SampleHandler_BeamFDStandardRecord.yaml new file mode 100644 index 00000000..a7c2c4fa --- /dev/null +++ b/Configs/Samples/SampleHandler_BeamFDStandardRecord.yaml @@ -0,0 +1,76 @@ +--- +SampleHandlerName: "BeamFDStandardRecord" + +MaCh3ModeConfig: "Configs/CovObjs/MaCh3Modes.yaml" +BinningFile: "Configs/Samples/BeamParameterBinning.yaml" + +InputFiles: + PrepSplineFile: false + MCEvents: + - Tag: numode_nonswap + FileLocation: "Inputs/DUNE_CAF_files/fardet-vd_mc_list-merge_physics_cafmaker_dunevd10kt_1x8x6_3view_30deg_runreco-nuenergy_geov3_root-tuple_merged_skip000000_lim000060_final_20241024T024241.root" + # If set to greater than 0, then events are skipped during reading with + # probability DownsampleFraction + # DownsampleFraction: 0.5 + # uncomment this if you just want to run with all the files + # - tag: "fhc_nonswap" + # file_location: "Inputs/DUNE_CAF_files/fardet-vd_mc_list-merge_physics_cafmaker_dunevd10kt_1x8x6_3view_30deg_runreco-nuenergy_geov3_root-tuple_merged_skip*.root"} + +NuOsc: + NuOscConfigFile: "Configs/NuOsc/CUDAProb3Linear.yaml" + EqualBinningPerOscChannel: false + +Samples: [ + BeamFD_numode_CCNumuLike, +] + +POT: + numode_nonswap: 1.36E23 + +BeamFD_numode_CCNumuLike: + SampleName: BeamFD_numode_CCNumuLike + SampleTitle: BeamFD_numode_CCNumuLike + + MCEventsTags: [ numode_nonswap, ] + + SelectionCuts: + # enum Sample { kRejected = 0, kNuMuCCLike, kNuECCLike, kNCLike }; + - KinematicStr: RecoSample + Bounds: [0.5, 1.5] + + Binning: + VarStr: RecoNeutrinoEnergy + Uniform: true + BinEdges: { linspace: { nb: 40, low: 0, up: 8 } } + + InputFiles: + mtupleprefix: "" + mtuplesuffix: "" + splineprefix: "" + splinesuffix: "" + + OscChannels: + - Name: "RHC_numu_x_numu" + LatexName: "RHC_numu_x_numu" + mtuplefile: [""] + splinefile: "" + nutype: 14 + oscnutype: 14 + - Name: "RHC_nue_x_nue" + LatexName: "RHC_nue_x_nue" + mtuplefile: [""] + splinefile: "" + nutype: 12 + oscnutype: 12 + - Name: "RHC_numubar_x_numubar" + LatexName: "RHC_numubar_x_numubar" + mtuplefile: [""] + splinefile: "" + nutype: -14 + oscnutype: -14 + - Name: "RHC_nuebar_x_nuebar" + LatexName: "RHC_nuebar_x_nuebar" + mtuplefile: [""] + splinefile: "" + nutype: -12 + oscnutype: -12 diff --git a/Samples/BeamFDStandardRecord/CMakeLists.txt b/Samples/BeamFDStandardRecord/CMakeLists.txt new file mode 100644 index 00000000..08fd4031 --- /dev/null +++ b/Samples/BeamFDStandardRecord/CMakeLists.txt @@ -0,0 +1,31 @@ +set(SOURCE_FILES + ReadEvents.cpp + SampleHandlerBeamFDStandardRecord.cpp +) +set(HEADERS + EventInfo.h + Projections.h + ReadEvents.h + SampleHandlerBeamFDStandardRecord.h +) + +add_library(SampleHandlerBeamFDStandardRecord SHARED + ${SOURCE_FILES} +) + +set_target_properties(SampleHandlerBeamFDStandardRecord PROPERTIES + PUBLIC_HEADER "${HEADERS}" + EXPORT_NAME SampleHandlerBeamFDStandardRecord) + +target_include_directories(SampleHandlerBeamFDStandardRecord PUBLIC + $ + $) + +install(TARGETS SampleHandlerBeamFDStandardRecord + EXPORT mach3dune-targets + LIBRARY DESTINATION lib/ + PUBLIC_HEADER DESTINATION include/Samples/BeamFDStandardRecord) + +target_link_libraries(SampleHandlerBeamFDStandardRecord PUBLIC SplinesDUNE SystematicsDUNE MaCh3::All MaCh3DUNECompilerOptions) +target_link_libraries(SampleHandlerBeamFDStandardRecord PRIVATE duneanaobj::StandardRecord duneanasel::all DUNEMaCh3Warnings ROOT::TreePlayer) + diff --git a/Samples/BeamFDStandardRecord/EventInfo.h b/Samples/BeamFDStandardRecord/EventInfo.h new file mode 100644 index 00000000..df371b8d --- /dev/null +++ b/Samples/BeamFDStandardRecord/EventInfo.h @@ -0,0 +1,57 @@ +#pragma once + +#include +#include +#include + +namespace dune::beamfd { + +struct CAFEventInfo { + struct Truth { + + int generator_mode; + double tgt_a; + double mach3_mode; + double is_cc; + + struct Neutrino { + int pdg, pdg_unosc; + double e; + } nu; + + } truth; + + struct Reconstructed { + + enum ESample { kRejected = 0, kNuMuCCLike, kNuECCLike, kNCLike }; + double sample; + + double e_nu; + + std::array vtx_pos_cm; + + } reco; +}; + +struct EventInfo : public CAFEventInfo { + + size_t tag_id; + bool is_numode; + int sample; + + struct SystInfo { + struct Flux { + std::vector focussing_weights; + std::vector hadprod_weights; + + double total_weight; + } flux; + + } syst; + + struct Weights { + double pot; + } weights; +}; + +} // namespace dune::beamfd diff --git a/Samples/BeamFDStandardRecord/Projections.h b/Samples/BeamFDStandardRecord/Projections.h new file mode 100644 index 00000000..79b6a984 --- /dev/null +++ b/Samples/BeamFDStandardRecord/Projections.h @@ -0,0 +1,68 @@ +#pragma once + +#include "Samples/BeamFDStandardRecord/EventInfo.h" + +#include "Manager/MaCh3Exception.h" +#include "Manager/MaCh3Logger.h" + +#include +#include + +namespace dune::beamfd { + +// below is ugly, but lets us define it only once and get the enum and both +// maps https://en.wikipedia.org/wiki/X_macro +#define LIST_OF_VARIABLES \ + X(TrueNeutrinoEnergy) \ + X(OscillationChannel) \ + X(TargetNucleus) \ + X(IsCC) \ + X(Mode) \ + X(RecoSample) \ + X(RecoNeutrinoEnergy) + +#define X(a) k##a, + +/// @brief Enum to identify kinematics +enum KinematicTypes { LIST_OF_VARIABLES }; + +#undef X +#define X(a) {#a, k##a}, +const std::unordered_map KinematicParametersDUNE = { + LIST_OF_VARIABLES}; + +#undef X +#define X(a) {k##a, #a}, +const std::unordered_map ReversedKinematicParametersDUNE = { + LIST_OF_VARIABLES}; + +#undef X +#undef LIST_OF_VARIABLES + +inline const double *ResolveKinematicEventMember(KinematicTypes KinPar, + EventInfo const &ev) { + switch (KinPar) { + case kTrueNeutrinoEnergy: + return &ev.truth.nu.e; + case kOscillationChannel: + return &ev.truth.mach3_mode; + case kTargetNucleus: + return &ev.truth.tgt_a; + case kIsCC: + return &ev.truth.is_cc; + case kMode: + return &ev.truth.mach3_mode; + + + case kRecoSample: + return &ev.reco.sample; + case kRecoNeutrinoEnergy: + return &ev.reco.e_nu; + + default: + MACH3LOG_ERROR("Did not recognise Kinematic Parameter type..."); + throw MaCh3Exception(__FILE__, __LINE__); + } +} + +} // namespace dune::beamfd diff --git a/Samples/BeamFDStandardRecord/ReadEvents.cpp b/Samples/BeamFDStandardRecord/ReadEvents.cpp new file mode 100644 index 00000000..a5ada620 --- /dev/null +++ b/Samples/BeamFDStandardRecord/ReadEvents.cpp @@ -0,0 +1,85 @@ +#include "Samples/BeamFDStandardRecord/ReadEvents.h" + +#include "Manager/Manager.h" + +_MaCh3_Safe_Include_Start_ //{ +#include "duneanaobj/StandardRecord/StandardRecord.h" +_MaCh3_Safe_Include_End_ //} + +#include "duneanasel/common/TruthTools.h" +#include "duneanasel/fd/beam/Observables.h" +#include "duneanasel/fd/beam/Selections.h" + +#include "TTreeReader.h" +#include "TTreeReaderValue.h" + +#include +#include + +namespace dune::beamfd { + + float GetPOT(TTree & tree) { + double pot = 0; + TTreeReader metardr(&tree); + TTreeReaderValue entry_pot(metardr, "pot"); + while (metardr.Next()) { +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wnull-dereference" + pot += *entry_pot; +#pragma GCC diagnostic pop + } + return float(pot); + } + + std::vector ReadEvents(TTree & tree, float downsamplefraction) { + + std::random_device r; + std::default_random_engine e1(r()); + std::uniform_real_distribution<> dis(0, 1); + + // Reco Variables + TTreeReader caf_reader(&tree); + + TTreeReaderValue sr(caf_reader, "rec"); + + std::vector events; + + size_t sr_it = 0; + while (caf_reader.Next()) { + + if (dis(e1) < downsamplefraction) { + continue; + } + + if (!sr->common.ixn.pandora.size()) { // empty event + continue; + } + + auto const &ixn = sr->common.ixn.pandora[0]; + EventInfo ev; + + auto sample = sel::beam::FD1::numode::ApplySelection(ixn); + ev.reco.sample = sample; + ev.reco.e_nu = proj::beam::FD1::ENuReco(ixn, sample); + auto const &det_pos_cm = proj::beam::det_pos_cm(ixn); + ev.reco.vtx_pos_cm = {det_pos_cm.x, det_pos_cm.y, det_pos_cm.z}; + + auto const &tixn = sr->mc.nu.front(); + + ev.truth.generator_mode = tixn.mode; + ev.truth.is_cc = tixn.iscc; + ev.truth.tgt_a = 40; // should do properly + + ev.truth.nu.pdg = tixn.pdg; + ev.truth.nu.pdg_unosc = tixn.pdgorig; + ev.truth.nu.e = tixn.E; + + if (ev.reco.sample > 0) { + events.push_back(ev); + } + + sr_it++; + } + return events; + } +} // namespace dune::beamfd diff --git a/Samples/BeamFDStandardRecord/ReadEvents.h b/Samples/BeamFDStandardRecord/ReadEvents.h new file mode 100644 index 00000000..dd0e86a8 --- /dev/null +++ b/Samples/BeamFDStandardRecord/ReadEvents.h @@ -0,0 +1,12 @@ +#pragma once + +#include "Samples/BeamFDStandardRecord/EventInfo.h" + +#include "TTree.h" + +namespace dune::beamfd { + +float GetPOT(TTree &); +std::vector ReadEvents(TTree &, float downsamplefraction = 0); + +} // namespace dune::beamfd diff --git a/Samples/BeamFDStandardRecord/SampleHandlerBeamFDStandardRecord.cpp b/Samples/BeamFDStandardRecord/SampleHandlerBeamFDStandardRecord.cpp new file mode 100644 index 00000000..21d3c448 --- /dev/null +++ b/Samples/BeamFDStandardRecord/SampleHandlerBeamFDStandardRecord.cpp @@ -0,0 +1,291 @@ +#include "Samples/BeamFDStandardRecord/SampleHandlerBeamFDStandardRecord.h" + +#include "Samples/BeamFDStandardRecord/ReadEvents.h" + +#include "Systematics/Beam.h" + +#include "Splines/BinnedSplineHandlerDUNE.h" + +#include + +#include + +namespace dune::beamfd { + +SampleHandlerBeamFDStandardRecord::SampleHandlerBeamFDStandardRecord( + std::string mc_version_, ParameterHandlerGeneric *ParHandler_, + const std::shared_ptr &Oscillator_) + : SampleHandlerBase(mc_version_, ParHandler_, Oscillator_) { + + KinematicParameters = &KinematicParametersDUNE; + ReversedKinematicParameters = &ReversedKinematicParametersDUNE; + + Initialise(); +} + +void SampleHandlerBeamFDStandardRecord::Init() {} + +void SampleHandlerBeamFDStandardRecord::SetupSplines() { + if (!ParHandler) { + return; + } + + ///@todo move all of the spline setup into core + int num_splines = 0; + for (int iSubSample = 0; iSubSample < int(SampleDetails.size()); + iSubSample++) { + num_splines += ParHandler->GetNumParamsFromSampleName( + GetSampleTitle(iSubSample), kSpline); + } + + if (num_splines > 0) { + MACH3LOG_INFO( + "Found {} splines for this sample so I will create a spline object", + num_splines); + SplineHandler = std::unique_ptr( + new BinnedSplineHandlerDUNE(ParHandler, Modes.get())); + InitialiseSplineObject(); + } else { + MACH3LOG_INFO( + "Found {} splines for this sample so I will not load or " + "evaluate splines", + ParHandler->GetNumParamsFromSampleName(SampleHandlerName, kSpline)); + SplineHandler = nullptr; + } +} + +void SampleHandlerBeamFDStandardRecord::RegisterFunctionalParameters() { + if (!ParHandler) { + return; + } + + if (ParHandler->GetNumParFromGroup("Flux")) { + RegisterIndividualFunctionalParameter( + DUNEMCEvents, syst::GetFluxFocussingParamNames(), + [](std::vector const &par_vals, EventInfo &ev) { + for (size_t i = 0; i < par_vals.size(); ++i) { + ev.syst.flux.total_weight *= + 1 + (par_vals[i] * (ev.syst.flux.focussing_weights[i] - 1)); + } + }); + + RegisterIndividualFunctionalParameter( + DUNEMCEvents, syst::GetFluxHadProdParamNames(), + [](std::vector const &par_vals, EventInfo &ev) { + for (size_t i = 0; i < par_vals.size(); ++i) { + ev.syst.flux.total_weight *= + 1 + (par_vals[i] * (ev.syst.flux.hadprod_weights[i] - 1)); + } + }); + } +} + +void SampleHandlerBeamFDStandardRecord::ResetShifts(int iEvent) { + auto &ev = DUNEMCEvents[iEvent]; + + // flux weights + ev.syst.flux.total_weight = 1.0; +} + +void SampleHandlerBeamFDStandardRecord::AddAdditionalWeightPointers() { + for (size_t i = 0; i < DUNEMCEvents.size(); ++i) { + MCEvents[i].total_weight_pointers.push_back(&(DUNEMCEvents[i].weights.pot)); + MCEvents[i].total_weight_pointers.push_back( + &(DUNEMCEvents[i].syst.flux.total_weight)); + } +} + +int SampleHandlerBeamFDStandardRecord::SetupExperimentMC() { + + MACH3LOG_INFO( + "-------------------------------------------------------------------"); + + bool do_flux_systematics = + ParHandler && ParHandler->GetNumParFromGroup("Flux"); + + std::vector> RecoSampleRanges; + + auto sample_name = Get>( + SampleManager->raw()["Samples"], __FILE__, __LINE__); + + for (int i = 0; i < GetNSamples(); i++) { + auto first_cut = Get( + SampleManager->raw()[sample_name[i]]["SelectionCuts"][0], __FILE__, + __LINE__); + auto kinstr = + Get(first_cut["KinematicStr"], __FILE__, __LINE__); + + if (kinstr == ReversedKinematicParametersDUNE.at(kRecoSample)) { + RecoSampleRanges.push_back( + Get>(first_cut["Bounds"], __FILE__, __LINE__)); + } else { + MACH3LOG_ERROR("Expected to only find a single Selection cut cutting on " + "{}. But found a cut on {}", + ReversedKinematicParametersDUNE.at(kRecoSample), kinstr); + throw MaCh3Exception(__FILE__, __LINE__); + } + } + + std::map>> + input_mc_event_descriptors; + + for (auto const &file_descriptor : + SampleManager->raw()["InputFiles"]["MCEvents"]) { + + auto tag = Get(file_descriptor["Tag"], __FILE__, __LINE__); + auto file_location = + Get(file_descriptor["FileLocation"], __FILE__, __LINE__); + auto downsamplefraction = GetFromManager( + file_descriptor["DownsampleFraction"], 0.0, __FILE__, __LINE__); + + input_mc_event_descriptors[tag].push_back( + std::make_pair(file_location, downsamplefraction)); + + MACH3LOG_INFO("-- Found input event descriptor: Tag: {}, FileLocation: {}, " + "DownsampleFraction: {}", + tag, file_location, downsamplefraction); + } + + for (auto const &[tag, input_files] : input_mc_event_descriptors) { + + float tag_pot = + Get(SampleManager->raw()["POT"][tag], __FILE__, __LINE__); + size_t tag_id = mc_tags.size(); + mc_tags.push_back(MCTag{tag, tag_pot}); + + bool is_numode = tag.find("numode") != std::string::npos; + + float tag_input_pot = 0; + + for (auto const &[filename, downsamplefraction] : input_files) { + if (filename.empty()) { + MACH3LOG_INFO("-- -- Skipping empty filename entry"); + continue; + } + TChain MetaChain("cafmaker/meta"); + MACH3LOG_INFO("-- -- Adding file descriptor to Meta TChain: {}", + filename); + if (!MetaChain.Add(filename.c_str(), -1)) { + MACH3LOG_ERROR("Could not add file {} to TChain, please check the file " + "exists and is readable", + filename); + throw MaCh3Exception(__FILE__, __LINE__); + } + + float file_descriptor_pot = GetPOT(MetaChain); + tag_input_pot += file_descriptor_pot * (1 - downsamplefraction); + MACH3LOG_INFO( + "-- -- Read {:.3G} input POT (with downsample weight of: {:.2f})", + file_descriptor_pot, (1 - downsamplefraction)); + } + + MACH3LOG_INFO("-- Read {:.3G} total POT for tag: {}, which has analysis " + "POT of {:.3G}", + tag_input_pot, tag, tag_pot); + + for (auto const &[filename, downsamplefraction] : input_files) { + if (filename.empty()) { + MACH3LOG_INFO("-- -- Skipping empty filename entry"); + continue; + } + + TChain CAFChain("cafmaker/cafTree"); + MACH3LOG_INFO("-- -- Adding file descriptor to cafTree TChain: {}", + filename); + if (!CAFChain.Add(filename.c_str(), -1)) { + MACH3LOG_ERROR("Could not add file {} to TChain, please check the file " + "exists and is readable", + filename); + throw MaCh3Exception(__FILE__, __LINE__); + } + + auto sample_evs = ReadEvents(CAFChain, downsamplefraction); + MACH3LOG_INFO( + "-- -- Read: {}/{} events (with downsample fraction of: {:.2f})", + sample_evs.size(), CAFChain.GetEntries(), downsamplefraction); + + // fix up any analysis specific information + for (auto &ev : sample_evs) { + + ev.tag_id = tag_id; + ev.is_numode = is_numode; + + ev.sample = -1; + for (size_t i = 0; i < RecoSampleRanges.size(); ++i) { + if ((ev.reco.sample > RecoSampleRanges[i][0]) && + (ev.reco.sample < RecoSampleRanges[i][1])) { + ev.sample = int(i); + } + } + + if (ev.sample < 0) { + continue; + } + + ev.truth.mach3_mode = + Modes->GetModeFromGenerator(std::abs(ev.truth.generator_mode)); + if (!ev.truth.is_cc) { + // Account for no ability to distinguish CC/NC + ev.truth.mach3_mode += 14; + } + if (ev.truth.mach3_mode > 15) { + // Account for no NCSingleKaon + ev.truth.mach3_mode -= 1; + } + + ev.weights.pot = tag_pot / tag_input_pot; + + ev.syst.flux.total_weight = 1; + if (do_flux_systematics) { + std::tie(ev.syst.flux.focussing_weights, + ev.syst.flux.hadprod_weights) = + syst::GetFluxVariationWeights(ev.truth.nu.pdg_unosc, + ev.truth.nu.e, true, is_numode); + } + + DUNEMCEvents.emplace_back(std::move(ev)); + } + } + } + + return int(DUNEMCEvents.size()); +} + +void SampleHandlerBeamFDStandardRecord::SetupMC() { + + size_t iEvent = 0; + for (auto const &ev : DUNEMCEvents) { + MCEvents[iEvent].isNC = !ev.truth.is_cc; + + MCEvents[iEvent].enu_true = ev.truth.nu.e; + MCEvents[iEvent].nupdg = ev.truth.nu.pdg; + MCEvents[iEvent].nupdgUnosc = ev.truth.nu.pdg_unosc; + + MCEvents[iEvent].NominalSample = ev.sample; + + iEvent++; + } +} + +void SampleHandlerBeamFDStandardRecord::InititialiseData() { + // Reweight MC to match + Reweight(); + // set asimov data + for (int iSample = 0; iSample < GetNSamples(); iSample++) { + AddData(iSample, GetMCArray(iSample)); + } +} + +const double *SampleHandlerBeamFDStandardRecord::GetPointerToKinematicParameter( + int KinematicVariable, int iEvent) const { + KinematicTypes KinPar = static_cast(KinematicVariable); + return ResolveKinematicEventMember(KinPar, DUNEMCEvents[iEvent]); +} + +double SampleHandlerBeamFDStandardRecord::ReturnKinematicParameter( + int KinematicVariable, int iEvent) const { + KinematicTypes KinPar = static_cast(KinematicVariable); + return *GetPointerToKinematicParameter(KinPar, iEvent); +} + +} // namespace dune::beamfd diff --git a/Samples/BeamFDStandardRecord/SampleHandlerBeamFDStandardRecord.h b/Samples/BeamFDStandardRecord/SampleHandlerBeamFDStandardRecord.h new file mode 100644 index 00000000..4424be86 --- /dev/null +++ b/Samples/BeamFDStandardRecord/SampleHandlerBeamFDStandardRecord.h @@ -0,0 +1,54 @@ +#pragma once + +#include "Samples/BeamFDStandardRecord/EventInfo.h" +#include "Samples/BeamFDStandardRecord/Projections.h" +#include "Samples/SampleHandlerBase.h" + +namespace dune::beamfd { + +class SampleHandlerBeamFDStandardRecord : public SampleHandlerBase { + +public: + SampleHandlerBeamFDStandardRecord( + std::string mc_version, ParameterHandlerGeneric *xsec_cov, + const std::shared_ptr &Oscillator); + + ~SampleHandlerBeamFDStandardRecord() {} + + std::vector DUNEMCEvents; + + struct MCTag { + std::string name; + float POT; + }; + std::vector mc_tags; + + int GetTagIndex(std::string const &tag){ + int index = 0; + for(auto const &t : mc_tags){ + if(t.name == tag){ + return index; + } + index++; + } + return -1; + } + + void CleanMemoryBeforeFit() {} + +protected: + void Init() override; + int SetupExperimentMC() override; + void SetupMC() override; + void AddAdditionalWeightPointers() override; + void SetupSplines() override; + void RegisterFunctionalParameters() override; + void ResetShifts(int iEvent) override; + void InititialiseData(); + + double ReturnKinematicParameter(int KinematicVariable, int iEvent) const; + + const double *GetPointerToKinematicParameter(int KinematicVariable, + int iEvent) const; +}; +} // namespace dune::beamfd diff --git a/Samples/CMakeLists.txt b/Samples/CMakeLists.txt index fd1160d5..1d9bc0b0 100644 --- a/Samples/CMakeLists.txt +++ b/Samples/CMakeLists.txt @@ -1,3 +1,5 @@ +add_subdirectory(BeamFDStandardRecord) + set(HEADERS ${CMAKE_CURRENT_SOURCE_DIR}/SampleHandlerAtm.h ${CMAKE_CURRENT_SOURCE_DIR}/SampleHandlerBeamFD.h @@ -21,7 +23,7 @@ set_target_properties(SamplesDUNE PROPERTIES PUBLIC_HEADER "${HEADERS}" EXPORT_NAME SamplesDUNE) -target_link_libraries(SamplesDUNE PUBLIC SplinesDUNE MaCh3::All MaCh3DUNECompilerOptions) +target_link_libraries(SamplesDUNE PUBLIC SplinesDUNE SampleHandlerBeamFDStandardRecord 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 84d354d9..cf1a2adc 100644 --- a/Samples/MaCh3DUNEFactory.cpp +++ b/Samples/MaCh3DUNEFactory.cpp @@ -2,6 +2,7 @@ // DUNE Handlers #include "Samples/SampleHandlerBeamFD.h" +#include "Samples/BeamFDStandardRecord/SampleHandlerBeamFDStandardRecord.h" #include "Samples/SampleHandlerBeamND.h" #include "Samples/SampleHandlerBeamNDGAr.h" #include "Samples/SampleHandlerAtm.h" @@ -12,6 +13,8 @@ SampleHandlerBase* GetMaCh3DuneInstance(std::string SampleType, std::string Samp SampleHandlerBase *Sample; if (SampleType == "BeamFD") { Sample = new SampleHandlerBeamFD(SampleConfig, param_handler.get(), BeamOscillator_); + } else if (SampleType == "BeamFDStandardRecord") { + Sample = new dune::beamfd::SampleHandlerBeamFDStandardRecord(SampleConfig, param_handler.get(), BeamOscillator_); } else if (SampleType == "BeamND") { if (beamNDCov.NDCov_FHC == nullptr || beamNDCov.NDCov_RHC == nullptr || beamNDCov.NDCov_all == nullptr) { diff --git a/Systematics/Beam.cpp b/Systematics/Beam.cpp new file mode 100644 index 00000000..4c58c56d --- /dev/null +++ b/Systematics/Beam.cpp @@ -0,0 +1,52 @@ +#include "Systematics/Beam.h" + +#include "duneanafluxtools/FluxWeighter.h" + +namespace dune::syst { +std::pair, std::vector> +GetFluxVariationWeights(int nu_pdg, double enu_true_GeV, bool is_FD, + bool is_numode, double off_axis_pos_m) { + + auto const &fw = FluxWeighter::Get(); + auto nucfg = fw.GetNuConfig(nu_pdg, !is_FD, is_numode); + + auto flux_focussing_systbin = + fw.GetFocussingBin(nu_pdg, enu_true_GeV, off_axis_pos_m, nucfg); + + auto flux_hadprod_systbin = + fw.GetFocussingBin(nu_pdg, enu_true_GeV, off_axis_pos_m, nucfg); + + std::pair, std::vector> weights; + for (size_t i = 0; i < fw.GetNFocussingParams(); i++) { + weights.first.push_back( + float(fw.GetFluxFocussingWeight(i, 1, flux_focussing_systbin, nucfg))); + } + for (size_t i = 0; i < fw.GetNHadProdPCAComponents(); i++) { + weights.second.push_back( + float(fw.GetFluxHadProdWeight(i, 1, flux_hadprod_systbin, nucfg))); + } + return weights; +} + +std::vector GetFluxFocussingParamNames() { + auto const &fw = FluxWeighter::Get(); + size_t nfocus_par = fw.GetNFocussingParams(); + std::vector focussing_par_names; + for (size_t focussing_par = 0; focussing_par < nfocus_par; focussing_par++) { + focussing_par_names.push_back(fw.GetFocussingParamName(focussing_par)); + } + return focussing_par_names; +} + +std::vector GetFluxHadProdParamNames() { + auto const &fw = FluxWeighter::Get(); + size_t Nhadprod_par = fw.GetNHadProdPCAComponents(); + std::vector hadprod_par_names; + for (size_t hadprod_par = 0; hadprod_par < Nhadprod_par; hadprod_par++) { + hadprod_par_names.push_back("Flux_HadProd_Param_" + + std::to_string(hadprod_par)); + } + return hadprod_par_names; +} + +} // namespace dune::syst diff --git a/Systematics/Beam.h b/Systematics/Beam.h new file mode 100644 index 00000000..fe9704e3 --- /dev/null +++ b/Systematics/Beam.h @@ -0,0 +1,16 @@ +#pragma once + +#include +#include +#include + +namespace dune::syst { + +std::vector GetFluxFocussingParamNames(); +std::vector GetFluxHadProdParamNames(); + +std::pair, std::vector> +GetFluxVariationWeights(int nu_pdg, double enu_true_GeV, bool is_FD, + bool is_numode, double off_axis_pos_m = 0); + +} // namespace dune::syst diff --git a/Systematics/CMakeLists.txt b/Systematics/CMakeLists.txt new file mode 100644 index 00000000..d3a6dfe9 --- /dev/null +++ b/Systematics/CMakeLists.txt @@ -0,0 +1,33 @@ +set(HEADERS + ${CMAKE_CURRENT_SOURCE_DIR}/Beam.h +) +set(SOURCE_FILES + Beam.cpp +) + +add_library(SystematicsDUNE SHARED + ${SOURCE_FILES} +) + +set_target_properties(SystematicsDUNE PROPERTIES + PUBLIC_HEADER "${HEADERS}" + EXPORT_NAME SystematicsDUNE) + +target_link_libraries(SystematicsDUNE PUBLIC MaCh3DUNECompilerOptions) +target_link_libraries(SystematicsDUNE PRIVATE duneanafluxtools::all DUNEMaCh3Warnings) + +target_include_directories(SystematicsDUNE PUBLIC + $ + $) + +install(TARGETS SystematicsDUNE + EXPORT mach3dune-targets + LIBRARY DESTINATION lib/ + PUBLIC_HEADER DESTINATION include/Samples) + +add_library(MACH3DUNE::SystematicsDUNE ALIAS SystematicsDUNE) + +add_executable(DumpBeamSystematicRatios DumpBeamSystematicRatios.cpp) +target_link_libraries(DumpBeamSystematicRatios SystematicsDUNE ROOT::Hist) + +install(TARGETS DumpBeamSystematicRatios DESTINATION bin) diff --git a/Systematics/DumpBeamSystematicRatios.cpp b/Systematics/DumpBeamSystematicRatios.cpp new file mode 100644 index 00000000..c61fa1ae --- /dev/null +++ b/Systematics/DumpBeamSystematicRatios.cpp @@ -0,0 +1,40 @@ +#include "TFile.h" +#include "TH1.h" + +#include "Systematics/Beam.h" + +int main() { + + TFile out("BeamSystematicRatios.root", "RECREATE"); + + for (auto const &spec : std::vector>{ + {"numu", 14}, {"nue", 12}, {"numubar", -14}, {"nuebar", -12}}) { + std::vector foc, hp; + + for (auto const pn : dune::syst::GetFluxFocussingParamNames()) { + foc.push_back( + new TH1D(Form("%s_numode_fd_%s", spec.first.c_str(), pn.c_str()), "", + 1000, 0, 10)); + } + + for (auto const pn : dune::syst::GetFluxHadProdParamNames()) { + hp.push_back( + new TH1D(Form("%s_numode_fd_%s", spec.first.c_str(), pn.c_str()), "", + 1000, 0, 10)); + } + + for (int i = 0; i < 1000; ++i) { + double enu = (i + 0.5) * 10.0 / 1000.0; + auto ws = + dune::syst::GetFluxVariationWeights(spec.second, enu, true, true); + for (int j = 0; j < ws.first.size(); ++j) { + foc[j]->SetBinContent(i + 1, ws.first[j]); + } + for (int j = 0; j < ws.second.size(); ++j) { + hp[j]->SetBinContent(i + 1, ws.second[j]); + } + } + } + out.Write(); + out.Close(); +} diff --git a/cmake/Modules/fetchduneanaobj.cmake b/cmake/Modules/fetchduneanaobj.cmake index 4650aa65..894750e8 100644 --- a/cmake/Modules/fetchduneanaobj.cmake +++ b/cmake/Modules/fetchduneanaobj.cmake @@ -11,6 +11,7 @@ macro(fetchduneanaobj DUNE_ANAOBJ_BRANCH) # Check if population has already been performed FetchContent_GetProperties(duneanaobj) if(NOT duneanaobj_POPULATED) + cmake_policy(SET CMP0169 OLD) # Fetch the content using previously declared details FetchContent_Populate(duneanaobj) @@ -63,5 +64,6 @@ macro(fetchduneanaobj DUNE_ANAOBJ_BRANCH) ) add_library(duneanaobj::StandardRecord ALIAS duneanaobj_StandardRecord) + add_library(duneanaobj::all ALIAS duneanaobj_StandardRecord) endmacro() diff --git a/cmake/Templates/setup.MaCh3DUNE.sh.in b/cmake/Templates/setup.MaCh3DUNE.sh.in index 1be82e68..ac7b3ff8 100644 --- a/cmake/Templates/setup.MaCh3DUNE.sh.in +++ b/cmake/Templates/setup.MaCh3DUNE.sh.in @@ -72,3 +72,7 @@ if [ -f "${MACH3DUNE_ROOT}/bin/setup.MaCh3.sh" ]; then echo "Sourcing MaCh3 Core" source "${MACH3DUNE_ROOT}/bin/setup.MaCh3.sh" fi + +if [ -z "${duneanafluxtools_ROOT}" ]; then + export duneanafluxtools_ROOT=${MACH3DUNE_ROOT} +fi diff --git a/plotting/PlottingConfig.yaml b/plotting/PlottingConfig.yaml new file mode 100644 index 00000000..91eacd20 --- /dev/null +++ b/plotting/PlottingConfig.yaml @@ -0,0 +1,90 @@ +ManagerOptions: + ## The translation config file to be used by the InputManager class to read files from each fitter + ## leave blank to use the default ({MACH3_ROOT}/plotting/universalTranslator.yaml) + translationConfig: "" + + ## The style definition config file to be used by the StyleManager class, sets things like colour palettes, fancy names for parameters etc. + ## leave blank to use the default ({MACH3_ROOT}/plotting/StyleConfig.yaml) + styleConfig: "" + +######################################### +PlotLLH: + defaultOutputName: "LLHScan" + + colorPalette: "OkabeItoReduced" + + ratioPlotSplit: 0.3 + yTitleOffset: 1.25 + sampleLabelThreshold: 0.0275 + lineWidth: 3 + totalOnSplitPlots: false + sameAxis: true + parameterTags: [] + sampleTags: [] + +######################################### +GetPostfitParamPlots: + defaultOutputName: "PostFitParamPlots" + + ## Some stylistic choices + colorPalette: "BlackBlueGreen" + prefitHistStyle: "redHatchedError" + plotLineWidth: 2 + refLineWidth: 2 + + detParYRange_low : 0.25 + detParYRange_high: 1.6 + oscParYRange_low : -3.0 + oscParYRange_high: 3.0 + + ## groups of flux parameters to be plotted together + fluxGroups: [ + ] + + ## define the groups of flux parameters, it goes groupName: [, , [, ]] + EXAMPLE_FULX_GROUP: ["COOL GROUP NAME", [0.7, 1.3], "EXAMPLE_BINNING", [0, 10]] + + ## the energy binnings to be used for flux plots + FluxBinning: + EXAMPLE_BINNING : [0.3, 0.4, 0.5, 0.6, 0.7, 1.0, 1.5, 2.5, 3.5, 5.0, 7.0, 30.0] + + + ## groups of parameters to plot together + paramGroups: [ + "Cross_Section", + "Detector_Systematic", + "Osc_Mixing_Angles", + "Osc_Mass_Splitting" + ] + + ## define the groups of parameters, should go : [, [, ], ["ParamName1", "ParamName2" ...] + Cross_Section: ["Cross Section", [-1.5, 2.0], + [ + "Norm_Param_0", + "BinnedSplineParam1", + "EResLep", + "EResTot" + ] + ] + + Detector_Systematic: ["Detector Systematics", [-1.5, 2.0], + [ + "Norm_Param_1", + "BinnedSplineParam2", + ] + ] + + Osc_Mixing_Angles: ["Oscillations - Mixing Angles", [0.0, 1.0], + [ + "sin2th_12", + "sin2th_23", + "sin2th_13", + ] + ] + + Osc_Mass_Splitting: ["Oscillations - Mass Splitting", [-0.02, 0.02], + [ + "delm2_12", + "delm2_23", + ] + ] diff --git a/plotting/StyleConfig.yaml b/plotting/StyleConfig.yaml new file mode 100644 index 00000000..87d0fb89 --- /dev/null +++ b/plotting/StyleConfig.yaml @@ -0,0 +1,115 @@ +######################################### +## can specify colour palettes that can be used here +## these are some nice ones but you can easily add more +ColorPallettes: + + OkabeItoReduced: [ + ## define the colour palette taken from https://www.color-hex.com/color-palette/49436 + + # the number of colours to include in the palette + [255.0], + + # now the RGB values for the palette + [ 0.0 , 0.25 , 0.5 , 0.75 ], ## stops + [ 0.83203125 , 0.796875 , 0.0 , 0.9375 ], ## Reds + [ 0.3671875 , 0.47265625 , 0.4453125 , 0.890625 ], ## greens + [ 0.0 , 0.65234375 , 0.6953125 , 0.2578125 ], ## blues + ] + + OkabeIto: [ + ## full Okabe-Ito colourblind friendly colour pallette defined at https://jfly.uni-koeln.de/color/ + [255.0], + [ 0.0 , 0.15 , 0.3 , 0.45 , 0.57 , 0.71 , 0.86 , 1.0 ], ## stops + [ 0.0 , 0.9 , 0.35 , 0.0 , 0.95 , 0.0 , 0.80 , 0.80 ], ## Reds + [ 0.0 , 0.6 , 0.7 , 0.6 , 0.90 , 0.45 , 0.40 , 0.60 ], ## greens + [ 0.0 , 0.0 , 0.9 , 0.5 , 0.25 , 0.70 , 0.00 , 0.70 ], ## blues + ] + + OkabeItoDiscrete: [ + ## full Okabe-Ito colourblind friendly colour pallette defined at https://jfly.uni-koeln.de/color/ + [8.0], + [ 0.0 , 0.15 , 0.3 , 0.45 , 0.57 , 0.71 , 0.86 , 1.0 ], ## stops + [ 0.0 , 0.9 , 0.35 , 0.0 , 0.95 , 0.0 , 0.80 , 0.80 ], ## Reds + [ 0.0 , 0.6 , 0.7 , 0.6 , 0.90 , 0.45 , 0.40 , 0.60 ], ## greens + [ 0.0 , 0.0 , 0.9 , 0.5 , 0.25 , 0.70 , 0.00 , 0.70 ], ## blues + ] + + BlackBlueGreen: [ + [4.0], + [ 0.0 , 0.33 , 0.66 , 1.0 ], ## stops + [ 0.0 , 0.0 , 0.0 , 0.0 ], ## Reds + [ 0.0 , 0.0 , 1.0 , 0.0 ], ## greens + [ 0.0 , 1.0 , 0.0 , 0.0 ], ## blues + ] + + RedWhiteBlue: [ + ## blue to white to red gradient, usually used for cov matrices and sigma vatiations + [255.0], + [ 0.00, 0.25, 0.50, 0.75, 1.00 ], + [ 0.00, 0.25, 1.00, 1.00, 0.50 ], + [ 0.00, 0.25, 1.00, 0.25, 0.00 ], + [ 0.50, 1.00, 1.00, 0.25, 0.00 ], + ] + + +######################################### +TH1Styles: +## define custom styles for TH1 type hists + redHatchedError: + MarkerColor: 632 ## marker colour (EM: there must be a way of getting root colours from strings, that would be a nice feature here) + MarkerStyle: 7 ## marker style + FillColor: 632 ## fill colour + FillStyle: 3003 ## fill style + LineColor: 632 ## line Colour + LineStyle: 1 ## line style + + +######################################### +## define prettier names for parameters and samples to use in plots +## can include TLatex for extra fanciness +PrettyNames: + ## first nice names for parameters + parameters: + Norm_Param_0: "#pi FSI" + Norm_Param_1: "det syst 0" + Norm_Param_2: "#Phi_{0}" + + BinnedSplineParam4: "XSec 0" + BinnedSplineParam5: "det syst 1" + + EResLep: "Energy Resolution Lepton" + EResTot: "Energy Resolution Total" + + sin2th_12: "sin^{2}#theta_{12}" + sin2th_23: "sin^{2}#theta_{23}" + sin2th_13: "sin^{2}#theta_{13}" + delm2_12: "#Delta m^{2}_{12}" + delm2_23: "#Delta m^{2}_{23}" + delta_cp: "#delta_{cp}" + + baseline: "L" + density: "#rho" + + ## now same for samples + samples: + SAMPLE_1: "sample 1" + SAMPLE_2: "sample 2" + SAMPLE_3: "sample 3" + + Tutorial_Beam: "FD Beam" + Tutorial_ATM: "FD ATM" + ND_CC0pi: "ND CC0#pi" + ND_CC1pip: "ND CC1#pi^{+}" + ND_CCPi0: "ND CC#pi^{0}" + ND_NC0pi: "ND NC0#pi" + ND_NC1pi0: "ND NC1#pi^{0}" + + ## now same for kinematics + kinematics: + RecoNeutrinoEnergy: "E_{#nu}^{reco}" + TrueQ2: "Q^{2}_{true}" + +######################################### +BinWidthScaleFactor: + RecoNeutrinoEnergy: 10 + TrueQ2: 100 diff --git a/plotting/universalTranslator.yaml b/plotting/universalTranslator.yaml new file mode 100644 index 00000000..8b3ca415 --- /dev/null +++ b/plotting/universalTranslator.yaml @@ -0,0 +1,237 @@ + +## Specify what each fitter looks like so that the plotting code knows where to look for different outputs +FitterSpec: + + ## first we declare what fitters the plotting code should be aware of + fitters: ["FITTER_1"] + + ## now for each one specified we give some details + FITTER_1: + + ################################################################### + ##### First we specify where to look for log likelihood scans ##### + ################################################################### + + ## Specify what type of objects are used to save the LLH scans + ## EM: could maybe just get this from the root file + ## on the other hand, might be useful if one day we move away from root, + ## then could specify different file types and just specify how to read them in the input manager + LLHObjectType: "TH1D" + + ## tell the code where to find likelihood scans for sample, penalty and total likelihood + sample_LLH: + ## we define a string that tells the code what path to look in to find likelihood scans + ## the token {PARAMETER} will be expanded to the fitter specific parameter name of a parameter + ## when the code is looking for scans for that parameter. + ## for example if there were parameters {p1, p2, p3} which in FITTER_1 are called {param1, param2, param3} + ## then the path specified below would tell the code to look for objects in the root file called: + ## - likelihood_scans/Sample_LLH/param1_sam + ## - likelihood_scans/Sample_LLH/param2_sam + ## - likelihood_scans/Sample_LLH/param3_sam + location: ["Sample_LLH/{PARAMETER}_sam"] + + penalty_LLH: + ## Additionally we can split the string two pars by using the ':' token + ## The first part should be a path to the directory within the root file + ## in this case that path is 'likelihood_scans/Sample_LLH' + ## The second part is a string to be matched to the end of the objects in that directory + ## this can be useful if there are a nuber of prefixes attached to the start of different parameters + ## e.g. re-using our example above, if the root file contained objects called + ## - likelihood_scans/Penalty_LLH/BLAAA_param1_pen + ## - likelihood_scans/Penalty_LLH/BLEEE_param2_pen + ## - likelihood_scans/Penalty_LLH/BLOOO_param3_pen + ## then the following path would correctly find the scans + location: ["xsec_LLH:{PARAMETER}_xsec",] + + total_LLH: + ## Fitters may split parameters across multiple different directories + ## + ## e.g + ## - likelihood_scans/TOTAL_LLH_DIR1/BLAAA_param1_pen + ## - likelihood_scans/TOTAL_LLH_DIR2/BLEEE_param2_pen + ## - likelihood_scans/TOTAL_LLH_DIR3/BLOOO_param3_pen + ## + ## to deal with this we can specify multiple locations to look for LLH scans as below + ## the code will search all of them for each parameter and will shout if it finds more than 1 matching pattern + location: [ + "Total_LLH:{PARAMETER}_full", + "Total_LLH_DIR_2:{PARAMETER}_full", + ] + + ## We also like to be able to beak our likelihood down by sample + ## here we tell the code where to find the per-sample likelihood scans + bySample_LLH: + ## we specify the locations in the same way as above, the only difference is that now there is an additional + ## token, {SAMPLE}, which will be expanded to the fitter specific sample name when the code + ## is looking for LLH scans for a particular sample + location: [ + "{SAMPLE}_LLH:{PARAMETER}{SAMPLE}", + ] + + + + ###################################################################### + ##### Now we specify where to look for processed post-fit errors ##### + ###################################################################### + + ## We often have multiple different types of errors that we would like to be able to plot + ## here we specify all the ones that the fitter might have + AvailablePostFitErrorTypes: ["Gauss", "HPD", "Prior"] + + ## We also specify one to use as a default, that will be assumed within the plotting code + ## if one is not specified + defaultPostFitErrorType: "HPD" + + ## Now we tell the code where to find the values for each of the error types + ## **** WARNING *** currently it is assumed that the values and errors are stored in a TH1D + ## with the bin labels specifying the name of the parameter. This will likely not be the + ## case for every fitter encountered. Should probably add the possibility of specifying + ### other formats. + postFitErrorTypes: + Gauss: + ## EM: specify possible locations and then scan through all of them when reading the input + ## code will look through the TH1D at each location and check it for the parameter it is looking for + location: ["param_spline_gaus", + "param_xsec_gaus" + ] + + HPD: + location: ["param_spline_HPD", + "param_xsec_HPD" + ] + + Prior: + location: ["param_spline_prefit", + "param_xsec_prefit" + ] + + 1dPosteriors: + location: ["Post_1d_hists:{PARAMETER}", + ] + + MCMCsteps: + location: [ + "posteriors", + ] + +## Now we define what samples are available for this particular experiment +Samples: + ## Here we define the master list of all the available samples + Samples: ["Tutorial_Beam", + "Tutorial_ATM", + "ND_CC0pi", + "ND_CC1pip", + "ND_CCPi0", + "ND_NC0pi", + "ND_NC1pi0", + ] + + SAMPLE_1: + FITTER_1: + + ## We can specify "tags" to apply to a sample + ## This allows you to later on in your plotting scripts get, for example, all the samples which, for example, + ## are related to some particular sub-detector by adding a tag for that detector + ## Or to samples which correspond to a particular tagged particle + tags: ["subDet1", "particleTag"] + + SAMPLE_2: + FITTER_1: + + SAMPLE_3: + Fitter_1: + + +## Now we define the parameters that are modelled in the experiment +Parameters: + + ## First set out the master list of parameters + Parameters: [ + "Norm_Param_0", + "Norm_Param_1", + "Norm_Param_2", + + "BinnedSplineParam1", + "BinnedSplineParam2", + "BinnedSplineParam3", + "BinnedSplineParam4", + "BinnedSplineParam5", + + "EResLep", + "EResTot", + + "sin2th_12", + "sin2th_23", + "sin2th_13", + "delm2_12", + "delm2_23", + "delta_cp", + + "baseline", + "density" + ] + + ## Now we tell the input manager how it should translate the parameter names to each of the fitters + ## Let's assume again that FITTER_1 uses the same parameter names defined above and so needs no + ## further specification, but that the names in FITTER_2 follow a different pattern, and also use + ## different names for the XSEC parameters when dealing with LLH scans and post fit errors + + Norm_Param_0: + ## We can also specify "tags" to apply to a parameter + ## This allows you to later on in your plotting scripts get, for example, all parameters which have the tag "pion" + ## and plot them together, or all parameters with both the tags "pion" and "FSI" + ## Let's pretend that this is a xsec parameter relating to pion FSI + tags: ["norm", "xsec", "FSI", "pion"] + + Norm_Param_1: + ## lets pretend this one is a detector systematic parameter + tags: ["systematic", "norm", "det-syst"] + + Norm_Param_2: + ## lets pretend this one is a flux parameter + tags: ["systematic", "norm", "flux"] + + BinnedSplineParam1: + tags: ["systematic", "spline", "xsec"] + + BinnedSplineParam2: + tags: ["systematic", "spline", "det-syst"] + + BinnedSplineParam3: + tags: ["systematic", "spline", "det-syst"] + + BinnedSplineParam4: + tags: ["systematic", "spline", "det-syst"] + + BinnedSplineParam5: + tags: ["systematic", "spline", "det-syst"] + + EResLep: + tags: ["systematic", "func", "xsec"] + + EResTot: + tags: ["systematic", "func", "xsec"] + + sin2th_12: + tags: ["oscillation", "angle"] + + sin2th_23: + tags: ["oscillation", "angle"] + + sin2th_13: + tags: ["oscillation", "angle"] + + delm2_12: + tags: ["oscillation", "mass-diff"] + + delm2_23: + tags: ["oscillation", "mass-diff"] + + delta_cp: + tags: ["oscillation"] + + baseline: + tags: ["oscillation", "experiment"] + + density: + tags: ["oscillation", "experiment"]