diff --git a/CMakeLists.txt b/CMakeLists.txt index 39354c6d..7ea3014f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -268,6 +268,7 @@ install(TARGETS MaCh3DUNECompilerOptions LIBRARY DESTINATION lib/) add_subdirectory(Splines) +add_subdirectory(Systematics) add_subdirectory(Samples) add_subdirectory(Apps) diff --git a/Systematics/CMakeLists.txt b/Systematics/CMakeLists.txt new file mode 100644 index 00000000..0615ddc2 --- /dev/null +++ b/Systematics/CMakeLists.txt @@ -0,0 +1 @@ +add_subdirectory(Flux) \ No newline at end of file diff --git a/Systematics/Flux/CMakeLists.txt b/Systematics/Flux/CMakeLists.txt new file mode 100644 index 00000000..98ec37e1 --- /dev/null +++ b/Systematics/Flux/CMakeLists.txt @@ -0,0 +1,15 @@ +add_library(OffAxisFluxUncertaintyHelper SHARED OffAxisFluxUncertaintyHelper.cxx) +target_link_libraries(OffAxisFluxUncertaintyHelper PUBLIC ROOT::Hist) + +set_target_properties(OffAxisFluxUncertaintyHelper PROPERTIES + PUBLIC_HEADER "OffAxisFluxUncertaintyHelper.h" + EXPORT_NAME OffAxisFluxUncertaintyHelper) + +target_include_directories(OffAxisFluxUncertaintyHelper PUBLIC + $ + $) + +install(TARGETS OffAxisFluxUncertaintyHelper + EXPORT mach3dune-targets + LIBRARY DESTINATION lib/ + PUBLIC_HEADER DESTINATION include/Systematics/Flux) \ No newline at end of file diff --git a/Systematics/Flux/FluxParameters_FD_and_PRISM.yaml b/Systematics/Flux/FluxParameters_FD_and_PRISM.yaml new file mode 100644 index 00000000..9b539ef5 --- /dev/null +++ b/Systematics/Flux/FluxParameters_FD_and_PRISM.yaml @@ -0,0 +1,1171 @@ +Systematics: +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: TargetUpstreamDegredation + ParameterName: TargetUpstreamDegredation + ParameterBounds: + - 0 + - 1.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: TargetTiltTransverseY + ParameterName: TargetTiltTransverseY + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: TargetTiltTransverseX + ParameterName: TargetTiltTransverseX + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: TargetLength + ParameterName: TargetLength + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: TargetDisplaceTransverseY + ParameterName: TargetDisplaceTransverseY + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: TargetDisplaceTransverseX + ParameterName: TargetDisplaceTransverseX + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: TargetDensity + ParameterName: TargetDensity + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: ProtonBeamTransverseY + ParameterName: ProtonBeamTransverseY + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: ProtonBeamTransverseX + ParameterName: ProtonBeamTransverseX + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: ProtonBeamRadius + ParameterName: ProtonBeamRadius + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: ProtonBeamAngleY + ParameterName: ProtonBeamAngleY + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: ProtonBeamAngleX + ParameterName: ProtonBeamAngleX + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HornWaterLayerThickness + ParameterName: HornWaterLayerThickness + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HornCurrent + ParameterName: HornCurrent + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HornCTiltTransverseY + ParameterName: HornCTiltTransverseY + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HornCTiltTransverseX + ParameterName: HornCTiltTransverseX + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HornCEllipticityXInducedBField + ParameterName: HornCEllipticityXInducedBField + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HornCEccentricityXInducedBField + ParameterName: HornCEccentricityXInducedBField + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HornCDisplaceLongitudinalZ + ParameterName: HornCDisplaceLongitudinalZ + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HornBTiltTransverseY + ParameterName: HornBTiltTransverseY + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HornBTiltTransverseX + ParameterName: HornBTiltTransverseX + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HornBEllipticityXInducedBField + ParameterName: HornBEllipticityXInducedBField + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HornBDisplaceLongitudinalZ + ParameterName: HornBDisplaceLongitudinalZ + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HornATiltTransverseY + ParameterName: HornATiltTransverseY + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HornATiltTransverseX + ParameterName: HornATiltTransverseX + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HornAEllipticityXInducedBField + ParameterName: HornAEllipticityXInducedBField + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HornAEccentricityXInducedBField + ParameterName: HornAEccentricityXInducedBField + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HornADisplaceLongitudinalZ + ParameterName: HornADisplaceLongitudinalZ + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: DecayPipeTiltY + ParameterName: DecayPipeTiltY + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: DecayPipeTiltX + ParameterName: DecayPipeTiltX + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: DecayPipeRadius + ParameterName: DecayPipeRadius + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: DecayPipeLength + ParameterName: DecayPipeLength + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: DecayPipeGeoBField + ParameterName: DecayPipeGeoBField + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: DecayPipeEllipticalCrossSectionYB + ParameterName: DecayPipeEllipticalCrossSectionYB + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: DecayPipeEllipticalCrossSectionXA + ParameterName: DecayPipeEllipticalCrossSectionXA + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: DecayPipeDisplaceTransverseY + ParameterName: DecayPipeDisplaceTransverseY + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: DecayPipeDisplaceTransverseX + ParameterName: DecayPipeDisplaceTransverseX + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: DecayPipe3SegmentBowingY + ParameterName: DecayPipe3SegmentBowingY + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: DecayPipe3SegmentBowingX + ParameterName: DecayPipe3SegmentBowingX + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HornCDisplaceTransverseY + ParameterName: HornCDisplaceTransverseY + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HornBDisplaceTransverseY + ParameterName: HornBDisplaceTransverseY + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HornADisplaceTransverseY + ParameterName: HornADisplaceTransverseY + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HornCDisplaceTransverseX + ParameterName: HornCDisplaceTransverseX + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HornBDisplaceTransverseX + ParameterName: HornBDisplaceTransverseX + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HornADisplaceTransverseX + ParameterName: HornADisplaceTransverseX + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HadronProduction_pca_19 + ParameterName: HadronProduction_pca_19 + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HadronProduction_pca_18 + ParameterName: HadronProduction_pca_18 + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HadronProduction_pca_17 + ParameterName: HadronProduction_pca_17 + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HadronProduction_pca_16 + ParameterName: HadronProduction_pca_16 + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HadronProduction_pca_15 + ParameterName: HadronProduction_pca_15 + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HadronProduction_pca_14 + ParameterName: HadronProduction_pca_14 + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HadronProduction_pca_13 + ParameterName: HadronProduction_pca_13 + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HadronProduction_pca_12 + ParameterName: HadronProduction_pca_12 + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HadronProduction_pca_11 + ParameterName: HadronProduction_pca_11 + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HadronProduction_pca_10 + ParameterName: HadronProduction_pca_10 + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HadronProduction_pca_9 + ParameterName: HadronProduction_pca_9 + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HadronProduction_pca_8 + ParameterName: HadronProduction_pca_8 + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HadronProduction_pca_7 + ParameterName: HadronProduction_pca_7 + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HadronProduction_pca_6 + ParameterName: HadronProduction_pca_6 + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HadronProduction_pca_5 + ParameterName: HadronProduction_pca_5 + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HadronProduction_pca_4 + ParameterName: HadronProduction_pca_4 + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HadronProduction_pca_3 + ParameterName: HadronProduction_pca_3 + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HadronProduction_pca_2 + ParameterName: HadronProduction_pca_2 + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HadronProduction_pca_1 + ParameterName: HadronProduction_pca_1 + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional +- Systematic: + Error: 1.0 + FlatPrior: false + Names: + FancyName: HadronProduction_pca_0 + ParameterName: HadronProduction_pca_0 + ParameterBounds: + - -3.0 + - 3.0 + ParameterGroup: Flux + ParameterValues: + Generated: 0 + PreFitValue: 0 + SampleNames: + - '*' + StepScale: + MCMC: 1 + Type: Functional diff --git a/Systematics/Flux/FluxVariationsValidation.pdf b/Systematics/Flux/FluxVariationsValidation.pdf new file mode 100644 index 00000000..e807370f Binary files /dev/null and b/Systematics/Flux/FluxVariationsValidation.pdf differ diff --git a/Systematics/Flux/OffAxisFluxUncertaintyHelper.cxx b/Systematics/Flux/OffAxisFluxUncertaintyHelper.cxx new file mode 100644 index 00000000..dac7e4b4 --- /dev/null +++ b/Systematics/Flux/OffAxisFluxUncertaintyHelper.cxx @@ -0,0 +1,563 @@ +#include "Systematics/Flux/OffAxisFluxUncertaintyHelper.h" + +#include "TAxis.h" +#include "TFile.h" +#include "TH1.h" + +#include +#include +#include +#include + +static OffAxisFluxUncertaintyHelper *globalFluxHelper = nullptr; + +OffAxisFluxUncertaintyHelper const &OffAxisFluxUncertaintyHelper::Get() { + if (!globalFluxHelper) { + globalFluxHelper = new OffAxisFluxUncertaintyHelper(); + globalFluxHelper->Initialize( + std::string(std::getenv("MACH3")) + + "/Systematics/Flux/flux_variations_FD_and_PRISM_2023.root"); + } + return *globalFluxHelper; +} + +bool check_axis(TAxis *ref, TAxis *other) { + if (ref->GetNbins() != other->GetNbins()) { + std::cout << "-- Mismatched TAxis: Reference NBins = " << ref->GetNbins() + << ", Other NBins = " << other->GetNbins() << std::endl; + return false; + } + + for (int bi = 0; bi < ref->GetNbins(); ++bi) { + if (std::fabs(ref->GetBinLowEdge(bi + 1) - other->GetBinLowEdge(bi + 1)) > + 1E-8) { + std::cout << "-- Mismatched TAxis: Low edge of bin " << bi + << ", Reference low edge: " << ref->GetBinLowEdge(bi + 1) + << ", other low edge: " << other->GetBinLowEdge(bi + 1) + << std::endl; + return false; + } + if (std::fabs(ref->GetBinUpEdge(bi + 1) - other->GetBinUpEdge(bi + 1)) > + 1E-8) { + std::cout << "-- Mismatched TAxis: Up edge of bin " << bi + << ", Reference up edge: " << ref->GetBinUpEdge(bi + 1) + << ", other up edge: " << other->GetBinUpEdge(bi + 1) + << std::endl; + return false; + } + } + return true; +} + +void OffAxisFluxUncertaintyHelper::Initialize(std::string const &filename, + bool verbose) { + + std::string ND_detector_tag = "ND"; + std::string ND_SpecHCRun_detector_tag = "_specrun_"; + std::string FD_detector_tag = "FD"; + std::string nu_mode_beam_tag = "nu"; + std::string nubar_mode_beam_tag = "nubar"; + std::string numu_species_tag = "numu"; + std::string nue_species_tag = "nue"; + std::string numubar_species_tag = "numubar"; + std::string nuebar_species_tag = "nuebar"; + + static std::string const location_tags[] = {ND_detector_tag, FD_detector_tag}; + static std::string const beam_mode_tags[] = {nu_mode_beam_tag, + nubar_mode_beam_tag}; + static std::string const species_tags[] = {numu_species_tag, nue_species_tag, + numubar_species_tag, + nuebar_species_tag}; + static std::string const spec_run_tags[] = {"_", ND_SpecHCRun_detector_tag}; + + if (verbose) { + std::cout << "Reading inputs from " << filename << std::endl; + } + + TFile *inpF = new TFile(filename.c_str()); + if (!inpF || !inpF->IsOpen()) { + std::cout << "[ERROR]: Couldn't open input file: " << filename << std::endl; + exit(1); + } + + { + focussing.OffAxisTAxes.clear(); + focussing.NDuncerts.clear(); + focussing.FDuncerts.clear(); + focussing.UncertLabels.clear(); + + std::string input_dir = "FluxParameters/Focussing"; + TDirectory *d = inpF->GetDirectory(input_dir.c_str()); + if (!d) { + std::cout << "[ERROR]: Couldn't open directory : " << input_dir + << " in input file: " << filename << std::endl; + exit(1); + } + + int param_id = 0; + + for (TObject *key : *d->GetListOfKeys()) { + + std::string syst(key->GetName()); + + if (verbose) { + std::cout << " Found parameter directory: " << input_dir << "/" << syst + << std::endl; + } + + int nucfg = kND_numu_numode; // = 0 + + std::stringstream input_dir_i(""); + input_dir_i << input_dir << (input_dir.size() ? "/" : "") << syst << "/"; + TDirectory *param_d = inpF->GetDirectory(input_dir_i.str().c_str()); + + std::unique_ptr oa_axis = + std::unique_ptr((TAxis *)param_d->Get("OffAxisTAxis")); + if (param_id && + !check_axis(focussing.OffAxisTAxes.front().get(), oa_axis.get())) { + exit(1); + } + focussing.OffAxisTAxes.emplace_back(std::move(oa_axis)); + + focussing.NDuncerts.emplace_back(); + focussing.FDuncerts.emplace_back(); + + std::vector>> param_NDuncerts; + std::vector> param_FDuncerts; + + focussing.UncertLabels.push_back(syst); + + for (size_t lt_it = 0; lt_it < 2; ++lt_it) { + std::string const &location_tag = location_tags[lt_it]; + for (size_t bm_it = 0; bm_it < 2; ++bm_it) { + std::string const &beam_mode_tag = beam_mode_tags[bm_it]; + for (size_t sr_it = 0; sr_it < 2; ++sr_it) { + std::string const &spec_run_tag = spec_run_tags[sr_it]; + for (size_t sp_it = 0; sp_it < 4; ++sp_it) { + std::string const &species_tag = species_tags[sp_it]; + + std::string hname = location_tag + "_" + beam_mode_tag + + spec_run_tag + species_tag; + + if (verbose) { + std::cout << " Reading nu configuration: " << hname + << std::endl; + } + + focussing.NDuncerts.back().emplace_back(); + focussing.FDuncerts.back().emplace_back(nullptr); + + if (nucfg < kFD_numu_numode) { + // Is ND (Any horn current for now) + // ND dir name is the same as the hist name + std::string nd_dir = input_dir_i.str() + hname; + std::vector> AllOffAxisShifts = + GetNDOffAxisShifts(inpF, nd_dir, hname); + + // check that for a given off axis slice, the energy binning is + // the same for every parameter + if (param_id) { + if (focussing.NDuncerts.front().at(nucfg).size() != + AllOffAxisShifts.size()) { + std::cout << "[ERROR]: Found differing number of off axis " + "slices for two flux focussing parameters: " + << focussing.UncertLabels.front() << " and " + << syst << std::endl; + } + for (size_t oab_i = 0; oab_i < AllOffAxisShifts.size(); + ++oab_i) { + if (!check_axis(focussing.NDuncerts.front() + .at(nucfg) + .at(oab_i) + ->GetXaxis(), + AllOffAxisShifts.at(oab_i)->GetXaxis())) { + std::cout << "[ERROR]: Found differing energy binning " + "for off axis bin " + << oab_i + << "slices for two flux focussing parameters: " + << focussing.UncertLabels.front() << " and " + << syst << std::endl; + exit(1); + } + } + } + + if (verbose) { + std::cout << " Found " << AllOffAxisShifts.size() + << " off axis bins!" << std::endl; + } + + focussing.NDuncerts.back().at(nucfg) = + std::move(AllOffAxisShifts); + nucfg += 1; + } else if (spec_run_tag == "_") { // Is FD and not 280kA run + std::unique_ptr h_in = + std::unique_ptr((TH1 *)param_d->Get(hname.c_str())); + if (!h_in) { + std::cout << "[WARN] Cannot find " << hname << std::endl; + continue; + } + + h_in->SetDirectory(nullptr); + + if (verbose) { + std::cout << " Found an FD variation." << std::endl; + } + + focussing.FDuncerts.back().at(nucfg) = std::move(h_in); + nucfg += 1; + } + } + } + } + } + param_id += 1; + } + focussing.NParams = param_id; + } + + { // a horrific stop gap until we have inputs with matching binning + + std::string ND_SpecHCRun_detector_tag = "_280kA_"; + + static std::string const spec_run_tags[] = {"_", ND_SpecHCRun_detector_tag}; + + hadprod.OffAxisTAxes.clear(); + hadprod.NDuncerts.clear(); + hadprod.FDuncerts.clear(); + + std::string input_dir = "FluxParameters/HadronProduction"; + TDirectory *d = inpF->GetDirectory(input_dir.c_str()); + + if (!d) { + std::cout << "[ERROR]: Couldn't open directory : " << input_dir + << " in input file: " << filename << std::endl; + exit(1); + } + + int param_id = 0; + + for (size_t i = 0; i < 20; ++i) { + + std::string syst = "HadronProduction_pca_" + std::to_string(param_id); + + if (verbose) { + std::cout << " Found parameter directory: " << input_dir << "/" << syst + << std::endl; + } + + int nucfg = kND_numu_numode; // = 0 + + std::stringstream input_dir_i(""); + input_dir_i << input_dir << (input_dir.size() ? "/" : "") << syst << "/"; + TDirectory *param_d = inpF->GetDirectory(input_dir_i.str().c_str()); + + hadprod.NDuncerts.emplace_back(); + hadprod.FDuncerts.emplace_back(); + hadprod.OffAxisTAxes.emplace_back(); + + std::vector>> param_NDuncerts; + std::vector> param_FDuncerts; + + for (size_t lt_it = 0; lt_it < 2; ++lt_it) { + std::string const &location_tag = location_tags[lt_it]; + for (size_t bm_it = 0; bm_it < 2; ++bm_it) { + std::string const &beam_mode_tag = beam_mode_tags[bm_it]; + for (size_t sr_it = 0; sr_it < 2; ++sr_it) { + std::string const &spec_run_tag = spec_run_tags[sr_it]; + for (size_t sp_it = 0; sp_it < 4; ++sp_it) { + std::string const &species_tag = species_tags[sp_it]; + + std::string hname = location_tag + "_" + beam_mode_tag + + spec_run_tag + species_tag; + + if (verbose) { + std::cout << " Reading nu configuration: " << hname + << std::endl; + } + + hadprod.NDuncerts.back().emplace_back(); + hadprod.FDuncerts.back().emplace_back(nullptr); + + if (nucfg < kFD_numu_numode) { + // Is ND (Any horn current for now) + // ND dir name is the same as the hist name + std::string nd_dir = input_dir_i.str() + hname; + std::vector> AllOffAxisShifts = + GetNDOffAxisShifts(inpF, nd_dir, hname); + + std::unique_ptr oa_axis = std::unique_ptr( + (TAxis *)inpF->GetDirectory(nd_dir.c_str()) + ->Get("OffAxisTAxis")); + if (param_id && + !check_axis(hadprod.OffAxisTAxes.front()[nucfg].get(), + oa_axis.get())) { + exit(1); + } + hadprod.OffAxisTAxes.back().emplace_back(std::move(oa_axis)); + + if (param_id && + !check_axis( + hadprod.NDuncerts.front().at(nucfg).front()->GetXaxis(), + AllOffAxisShifts.front()->GetXaxis())) { + exit(1); + } + + if (verbose) { + std::cout << " Found " << AllOffAxisShifts.size() + << " off axis bins!" << std::endl; + } + + hadprod.NDuncerts.back().at(nucfg) = + std::move(AllOffAxisShifts); + nucfg += 1; + } else if (spec_run_tag == "_") { // Is FD and not 280kA run + std::unique_ptr h_in = + std::unique_ptr((TH1 *)param_d->Get(hname.c_str())); + if (!h_in) { + std::cout << "[WARN] Cannot find " << hname << std::endl; + continue; + } + + h_in->SetDirectory(nullptr); + + if (verbose) { + std::cout << " Found an FD variation." << std::endl; + } + + hadprod.FDuncerts.back().at(nucfg) = std::move(h_in); + nucfg += 1; + } + } + } + } + } + param_id += 1; + } + hadprod.NPCAComponents = param_id; + } +} + +int OffAxisFluxUncertaintyHelper::GetNuConfig(int nu_pdg, bool IsND, + bool IsNuMode, + bool isSpecHCRun) const { + + int nucfg = kUnhandled; + + switch (nu_pdg) { + case 14: { + if (IsND) { + nucfg = IsNuMode + ? (isSpecHCRun ? kND_SpecHCRun_numu_numode : kND_numu_numode) + : (isSpecHCRun ? kND_SpecHCRun_numu_nubarmode + : kND_numu_nubarmode); + } else { + nucfg = IsNuMode ? kFD_numu_numode : kFD_numu_nubarmode; + } + break; + } + case -14: { + if (IsND) { + nucfg = IsNuMode ? (isSpecHCRun ? kND_SpecHCRun_numubar_numode + : kND_numubar_numode) + : (isSpecHCRun ? kND_SpecHCRun_numubar_nubarmode + : kND_numubar_nubarmode); + } else { + nucfg = IsNuMode ? kFD_numubar_numode : kFD_numubar_nubarmode; + } + break; + } + case 12: { + if (IsND) { + nucfg = + IsNuMode + ? (isSpecHCRun ? kND_SpecHCRun_nue_numode : kND_nue_numode) + : (isSpecHCRun ? kND_SpecHCRun_nue_nubarmode : kND_nue_nubarmode); + } else { + nucfg = IsNuMode ? kFD_nue_numode : kFD_nue_nubarmode; + } + break; + } + case -12: { + if (IsND) { + nucfg = IsNuMode ? (isSpecHCRun ? kND_SpecHCRun_nuebar_numode + : kND_nuebar_numode) + : (isSpecHCRun ? kND_SpecHCRun_nuebar_nubarmode + : kND_nuebar_nubarmode); + } else { + nucfg = IsNuMode ? kFD_nuebar_numode : kFD_nuebar_nubarmode; + } + break; + } + } + + return nucfg; +} + +std::vector> +OffAxisFluxUncertaintyHelper::GetNDOffAxisShifts(TFile *f, std::string nd_dir, + std::string hname) const { + std::vector> OffAxisUncerts; + + TDirectory *d = f->GetDirectory(nd_dir.c_str()); + if (!d) { + std::cout << "[ERROR]: Failed to open directory: " << nd_dir << std::endl; + exit(1); + } + int n_offaxis = d->GetListOfKeys()->GetSize(); + + for (int oa = 0; oa < n_offaxis; oa++) { + std::string hname_oa = hname + "_" + std::to_string(oa); + std::unique_ptr h_oa = + std::unique_ptr((TH1 *)d->Get(hname_oa.c_str())); + if (!h_oa) { + break; + } + h_oa->SetDirectory(nullptr); + OffAxisUncerts.emplace_back(std::move(h_oa)); + } + + return OffAxisUncerts; +} + +int OffAxisFluxUncertaintyHelper::GetFocussingBin(int nu_config, double enu_GeV, + double off_axis_pos_m) const { + + if (nu_config < kFD_numu_numode) { + // Sign flip in off axis position + int bin_oa = focussing.OffAxisTAxes.front()->FindFixBin(off_axis_pos_m); + if ((bin_oa == 0) || + (bin_oa == (focussing.OffAxisTAxes.front()->GetNbins() + 1))) { + // flow bin + return kInvalidBin; + } + bin_oa--; + auto &ehist = focussing.NDuncerts.front().at(nu_config).at(bin_oa); + int bin_e = ehist->FindFixBin(enu_GeV); + if ((bin_e == 0) || (bin_e == (ehist->GetXaxis()->GetNbins() + 1))) { + // flow bin + return kInvalidBin; + } + return bin_oa * 1000 + bin_e; + + } else { + if (focussing.FDuncerts.front().size() < 1) { + std::cout << "[WARN] no param_id" << std::endl; + } + if (!focussing.FDuncerts.front().at(nu_config)) { + std::cout << "[WARN] no nu_config = " << nu_config << std::endl; + } + auto &ehist = focussing.FDuncerts.front().at(nu_config); + int bin_e = ehist->FindFixBin(enu_GeV); + if ((bin_e == 0) || (bin_e == (ehist->GetXaxis()->GetNbins() + 1))) { + // flow bin + return kInvalidBin; + } + return bin_e; + } +} + +int OffAxisFluxUncertaintyHelper::GetHadProdBin(int nu_config, double enu_GeV, + double off_axis_pos_m) const { + if (nu_config < kFD_numu_numode) { + int bin_oa = + hadprod.OffAxisTAxes.front()[nu_config]->FindFixBin(off_axis_pos_m); + if ((bin_oa == 0) || + (bin_oa == (hadprod.OffAxisTAxes.front()[nu_config]->GetNbins() + 1))) { + // flow bin + return kInvalidBin; + } + bin_oa--; + auto &ehist = hadprod.NDuncerts.front().at(nu_config).at(bin_oa); + int bin_e = ehist->FindFixBin(enu_GeV); + if ((bin_e == 0) || (bin_e == (ehist->GetXaxis()->GetNbins() + 1))) { + // flow bin + return kInvalidBin; + } + return bin_oa * 1000 + bin_e; + + } else { + if (hadprod.FDuncerts.front().size() < 1) { + std::cout << "[WARN] no param_id" << std::endl; + } + if (!hadprod.FDuncerts.front().at(nu_config)) { + std::cout << "[WARN] no nu_config = " << nu_config << std::endl; + } + auto &ehist = hadprod.FDuncerts.front().at(nu_config); + int bin_e = ehist->FindFixBin(enu_GeV); + if ((bin_e == 0) || (bin_e == (ehist->GetXaxis()->GetNbins() + 1))) { + // flow bin + return kInvalidBin; + } + return bin_e; + } +} + +double OffAxisFluxUncertaintyHelper::GetFluxFocussingWeight(size_t param_id, + double param_val, + int nucfg, + int bin) const { + if (nucfg == kUnhandled) { + return 1; + } + + if (bin == kInvalidBin) { + return 1; + } + + if (nucfg < kFD_numu_numode) { + int bin_oa = bin / 1000; + int bin_e = bin % 1000; + return 1 + param_val * focussing.NDuncerts.at(param_id) + .at(nucfg) + .at(bin_oa) + ->GetBinContent(bin_e); + } else { + return 1 + + param_val * + focussing.FDuncerts.at(param_id).at(nucfg)->GetBinContent(bin); + } +} + +double OffAxisFluxUncertaintyHelper::GetFluxHadProdWeight(size_t param_id, + double param_val, + int nucfg, + int bin) const { + if (nucfg == kUnhandled) { + return 1; + } + + if (bin == kInvalidBin) { + return 1; + } + + if (nucfg < kFD_numu_numode) { + int bin_oa = bin / 1000; + int bin_e = bin % 1000; + return 1 + param_val * hadprod.NDuncerts.at(param_id) + .at(nucfg) + .at(bin_oa) + ->GetBinContent(bin_e); + } else { + return 1 + param_val * + hadprod.FDuncerts.at(param_id).at(nucfg)->GetBinContent(bin); + } +} + +std::vector +OffAxisFluxUncertaintyHelper::GetFluxFocussingOffAxisBinning() { + std::vector rtn{focussing.OffAxisTAxes[0]->GetBinLowEdge(1)}; + for (int i = 0; i < focussing.OffAxisTAxes[0]->GetNbins(); ++i) { + rtn.push_back(focussing.OffAxisTAxes[0]->GetBinUpEdge(i + 1)); + } + return rtn; +} +std::vector +OffAxisFluxUncertaintyHelper::GetFluxHadProdOffAxisBinning(int nu_config) { + std::vector rtn{ + hadprod.OffAxisTAxes[0][nu_config]->GetBinLowEdge(1)}; + for (int i = 0; i < hadprod.OffAxisTAxes[0][nu_config]->GetNbins(); ++i) { + rtn.push_back(hadprod.OffAxisTAxes[0][nu_config]->GetBinUpEdge(i + 1)); + } + return rtn; +} diff --git a/Systematics/Flux/OffAxisFluxUncertaintyHelper.h b/Systematics/Flux/OffAxisFluxUncertaintyHelper.h new file mode 100644 index 00000000..a37da501 --- /dev/null +++ b/Systematics/Flux/OffAxisFluxUncertaintyHelper.h @@ -0,0 +1,110 @@ +#pragma once + +// forward declaration +class TFile; +class TH1; +class TAxis; + +#include +#include +#include +#include +#include + +// adapted from CAFAna class from: +// https://github.com/DUNE/lblpwgtools/commit/195f456add17d2b9af67d6c9fc05cb035b843fe0 +class OffAxisFluxUncertaintyHelper { +public: + ~OffAxisFluxUncertaintyHelper() {} + + // nu_config options + static int const kND_numu_numode = 0; + static int const kND_nue_numode = 1; + static int const kND_numubar_numode = 2; + static int const kND_nuebar_numode = 3; + + static int const kND_SpecHCRun_numu_numode = 4; + static int const kND_SpecHCRun_nue_numode = 5; + static int const kND_SpecHCRun_numubar_numode = 6; + static int const kND_SpecHCRun_nuebar_numode = 7; + + static int const kND_numu_nubarmode = 8; + static int const kND_nue_nubarmode = 9; + static int const kND_numubar_nubarmode = 10; + static int const kND_nuebar_nubarmode = 11; + + static int const kND_SpecHCRun_numu_nubarmode = 12; + static int const kND_SpecHCRun_nue_nubarmode = 13; + static int const kND_SpecHCRun_numubar_nubarmode = 14; + static int const kND_SpecHCRun_nuebar_nubarmode = 15; + + static int const kFD_numu_numode = 16; + static int const kFD_nue_numode = 17; + static int const kFD_numubar_numode = 18; + static int const kFD_nuebar_numode = 19; + + static int const kFD_numu_nubarmode = 20; + static int const kFD_nue_nubarmode = 21; + static int const kFD_numubar_nubarmode = 22; + static int const kFD_nuebar_nubarmode = 23; + + static int const kUnhandled = 24; + + static int const kInvalidBin = std::numeric_limits::max(); + + void Initialize(std::string const &filename, bool verbose = false); + + static OffAxisFluxUncertaintyHelper const &Get(); + size_t GetNFocussingParams() const { return focussing.NDuncerts.size(); } + std::string GetFocussingParamName(size_t i) const { + return focussing.UncertLabels.at(i); + } + + size_t GetNHadProdPCAComponents() const { return hadprod.NDuncerts.size(); } + + int GetNuConfig(int nu_pdg, bool IsND, bool IsNuMode, + bool isSpecHCRun = false) const; + + std::vector> + GetNDOffAxisShifts(TFile *f, std::string nd_dir, std::string hname) const; + + int GetFocussingBin(int nu_config, double enu_GeV, + double off_axis_pos_m) const; + + int GetHadProdBin(int nu_config, double enu_GeV, double off_axis_pos_m) const; + + double GetFluxFocussingWeight(size_t param_id, double param_val, + int nu_config, int bin) const; + + double GetFluxHadProdWeight(size_t param_id, double param_val, int nu_config, + int bin) const; + + std::vector GetFluxFocussingOffAxisBinning(); + std::vector GetFluxHadProdOffAxisBinning(int nu_config); + + struct { + size_t NParams; + + // param + std::vector> OffAxisTAxes; + // param | nucfg | OA bin | E bin + std::vector>>> NDuncerts; + // param | nucfg | E bin + std::vector>> FDuncerts; + + std::vector UncertLabels; + + } focussing; + + struct { + size_t NPCAComponents; + + // param | nucfg + std::vector>> OffAxisTAxes; + // param | nucfg | OA bin | E bin + std::vector>>> NDuncerts; + // param | nucfg | E bin + std::vector>> FDuncerts; + + } hadprod; +}; diff --git a/Systematics/Flux/README.md b/Systematics/Flux/README.md new file mode 100644 index 00000000..8e48099f --- /dev/null +++ b/Systematics/Flux/README.md @@ -0,0 +1,114 @@ +# Flux Systematic Handling + +## Input Provenance + +The flux ratios used here were taken from CAFAna/prism in October 2025. The +inputs we use here were synthesised by combining two sets of histograms, an +older set ([`flux_shifts_OffAxis.root`](https://github.com/DUNE/lblpwgtools/blob/b9fd2dc0acca52a6c20a246cb3ceea313c0a4529/CAFAna/Systs/flux_shifts_OffAxis.root)) +containing the flux ratios from throws of the hadron production model (via +the package PPFX) and an updated set +([`flux_shifts_OffAxis2023.root`](https://github.com/DUNE/lblpwgtools/blob/b9fd2dc0acca52a6c20a246cb3ceea313c0a4529/CAFAna/Systs/flux_shifts_OffAxis2023.root)) +containing re-evaluated focussing, alignment and target condition systematics. +The auxilliary script [`combine.C`](combine.C) contains the logic for creating +our single input file. Note that the hadron production systematics were stored +in a non-standard ROOT TH2 format and require extra code to extract, +`combine.C` converts them to standard root objects in an equivalent format to +the focussing, alignment, and target condition inputs. + +## Usage in Samples + +The main interface to the systematic parameter responses is through the +[`OffAxisFluxUncertaintyHelper`](OffAxisFluxUncertaintyHelper.h) class. + +*At initialisation time*, when events are read in, samples should retrieve two +systematic bin identifiers for each event, one for the focussing, alignment, +and target condition parameters and one for the hadron production parameters. +A 'neutrino configuration' integer is used to identify the beam mode, neutrino +species, and detector and should also be stored with the event to avoid +re-calculation at step time. An example snippet for doing so is shown below: + +```c++ +//for each event, i + duneobj->flux_syst_nu_config.push_back( + flux_helper->GetNuConfig(duneobj->nupdgUnosc[i], isND, isFHC)); + + // xdir in detsim and flux syst inputs are opposite. + double syst_xpos_m = -(_det_x + _vtx_x)/100.0; + duneobj->flux_focussing_syst_bin.push_back(flux_helper->GetFocussingBin( + duneobj->flux_syst_nu_config.back(), enu_GeV, syst_xpos_m)); + duneobj->flux_hadprod_syst_bin.push_back(flux_helper->GetHadProdBin( + duneobj->flux_syst_nu_config.back(), enu_GeV, syst_xpos_m)); +``` + +The exact names of various variables may change in your particular sample, +but hopefully, the usage is clear. + +*At step time*, when weights should be retrieved, the previously determined +neutrino configuration and systematic bin can be proferred to the interface, +along with parameter identifiers and values, in exchange for response weights +to be used in predicting varied event rates. An example snippet is shown +below: + +```c++ +for (int i = 0; i < int(flux_helper->GetNFocussingParams()); i++) { + //... + double w = flux_helper->GetFluxFocussingWeight( + i, *par, dunemcSamples[iSample].flux_syst_nu_config[iEvent], + dunemcSamples[iSample].flux_focussing_syst_bin[iEvent]); + //... +} +//... +for (int i = 0; i < int(flux_helper->GetNHadProdPCAComponents()); i++) { + //... + double w = flux_helper->GetFluxHadProdWeight( + i, *par, dunemcSamples[iSample].flux_syst_nu_config[iEvent], + dunemcSamples[iSample].flux_hadprod_syst_bin[iEvent]); + //... +} +``` + +The response weight from each parameter represents independent variations and +so can safetly multiplied together to produce a total event weight for some +multi-parameter variation. + +## Systematic Parameter Configuration + +A `yaml` file containing the parameter configurations can be found in +[FluxParameters_FD_and_PRISM.yaml](FluxParameters_FD_and_PRISM.yaml). A python +script [`emit_flux_yaml.py`](emit_flux_yaml.py) is provided to emit this +configuration given the input file of varied flux ratios +([flux_variations_FD_and_PRISM_2023.root](flux_variations_FD_and_PRISM_2023.root)). + +## Validations + +The [`OffAxisFluxUncertaintyHelper`](OffAxisFluxUncertaintyHelper.h) class has +been validated in MaCh3_DUNE by comparing the histograms directly from the +input file with ratios built via direct calling of the interface code (outside +of MaCh3) with the results of running [SigmaVariation](SigmaVariation.cpp) +inside of MaCh3. At the time of committing, these validations passed for ND on +axis and one off axis position. A script is provided to emit validation plots w +hen proffered an output file from a SigmaVariation run and the output of +running the interface testing macro [`test_read_flux_systs.C`](test_read_flux_systs.C). + +A pdf containing validation plots can be found here: +[FluxVariationsValidation](FluxVariationsValidation.pdf). + +## Auxilliary Scripts + +* [`combine.C`](combine.C) + + A script to combine two sets of input files from CAFAna into a single set + of inputs used by this library. Requires + [TH2Jagged](https://github.com/luketpickering/TH2Jagged) to be built as a + subdirectory of the directory that it is run in. +* [`test_read_flux_systs.C`](test_read_flux_systs.C) + + A script to test that the interface code can read the input histograms. + Builds the interface code standalone with cling, and then loops over all the + parameters writing out some binning and parameter name information to prove + that it's managed to read the inputs. Also outputs a root file containing some + flux ratios retrieved via the interface to test that the right inputs are + being read for the right physics quantities. Used by `plot_flux_sigvar.py`. +* [`emit_flux_yaml.py`](emit_flux_yaml.py) + + A script to read the file containing the input histograms and dump out a + systematic parameter `yaml` file. +* [`plot_flux_sigvar.py`](plot_flux_sigvar.py) + + Validation plotting script, described in [Validations](#validations). diff --git a/Systematics/Flux/combine.C b/Systematics/Flux/combine.C new file mode 100644 index 00000000..a24f22d2 --- /dev/null +++ b/Systematics/Flux/combine.C @@ -0,0 +1,86 @@ +#include "TH2Jagged/build/Linux/include/TH2Jagged.h" + +#pragma cling load("TH2Jagged/build/Linux/lib/libTH2Jagged.so") + +void TH2JToDir(TH2Jagged *thj, TDirectory *dir) { + dir->WriteTObject(thj->fUniformAxis.Clone("OffAxisTAxis"), "OffAxisTAxis"); + for (size_t ui = 0; ui < thj->fUniformAxis.GetNbins(); ++ui) { + auto thf = thj->NonUniformSlice(ui); + dir->WriteTObject( + thf, (std::string(thj->GetName()) + "_" + std::to_string(ui)).c_str()); + delete thf; + } +} + +// adapted from https://root.cern/doc/v632/copyFiles_8C.html +void CopyDir(TDirectory *source, TDirectory *dest, bool pca_limit) { + + // loop on all entries of this directory + TKey *key; + // Loop in reverse order to make sure that the order of cycles is + // preserved. + TIter nextkey(source->GetListOfKeys(), kIterBackward); + while ((key = (TKey *)nextkey())) { + + const char *classname = key->GetClassName(); + TClass *cl = gROOT->GetClass(classname); + + if (!cl) + continue; + + if (cl->InheritsFrom(TDirectory::Class())) { + + if (pca_limit) { + if (std::string(key->GetName()).substr(0, 9) == "param_pca") { + int pca_num = std::stoi(std::string(key->GetName()).substr(10)); + if (pca_num >= 20) { + continue; + } + } else { + continue; + } + } + + std::string dirname = key->GetName(); + if (dirname.substr(0, 6) == "param_") { + dirname = dirname.substr(6); + } + if (dirname.substr(0, 3) == "pca") { + dirname = "HadronProduction_" + dirname; + } + + std::cout << "Copying directory: " << source->GetName() << "/" << dirname + << " to " << dest->GetName() << std::endl; + CopyDir(source->GetDirectory(key->GetName()), + dest->mkdir(dirname.c_str()), pca_limit); + } else { + if (std::string(key->GetName()) == "param_names") { + continue; + } + std::cout << " writing object: " << key->GetName() << " of type " + << cl->GetName() << " to " << dest->GetName() << std::endl; + if (std::string(cl->GetName()) == "TH2Jagged") { + TH2JToDir(static_cast *>(key->ReadObj()), + dest->mkdir(key->GetName())); + } else { + dest->WriteTObject(key->ReadObj(), key->GetName()); + } + } + } +} + +void combine() { + TFile fcomb("flux_variations_FD_and_PRISM_2023.root", "Recreate"); + TFile ffoc("flux_shifts_OffAxis2023.root", "READ"); + + TFile fhp("flux_shifts_OffAxis.root", "READ"); + + CopyDir(ffoc.GetDirectory("FluxParameters"), + fcomb.mkdir("FluxParameters")->mkdir("Focussing"), false); + + CopyDir(fhp.GetDirectory("FluxParameters"), + fcomb.GetDirectory("FluxParameters")->mkdir("HadronProduction"), + true); + + fcomb.Write(); +} \ No newline at end of file diff --git a/Systematics/Flux/emit_flux_yaml.py b/Systematics/Flux/emit_flux_yaml.py new file mode 100644 index 00000000..1482edce --- /dev/null +++ b/Systematics/Flux/emit_flux_yaml.py @@ -0,0 +1,52 @@ +# --- +# Systematics: +# - Systematic: +# SampleNames: ["ND_*", "FD_*"] +# Error: 1.0 +# FlatPrior: false +# Names: +# FancyName: MAQE +# ParameterName: MAQE +# ParameterBounds: +# - -4.0 +# - 4.0 +# ParameterGroup: Xsec +# ParameterValues: +# Generated: 0.0 +# PreFitValue: 0.0 +# SplineInformation: +# Mode: +# - 0 +# SplineName: maqe +# StepScale: +# MCMC: 0.002 +# Type: Spline + +import yaml + +import ROOT + +Systematics = [] + +fin = ROOT.TFile.Open("flux_variations_FD_and_PRISM_2023.root") + +dfp = fin.GetDirectory("FluxParameters") +for pgroups in dfp.GetListOfKeys(): + for p in dfp.GetDirectory(pgroups.GetName()).GetListOfKeys(): + Systematic = {} + Systematic["SampleNames"] = ["*"] + Systematic["Error"] = 1.0 + Systematic["FlatPrior"] = False + Systematic["Names"] = { "FancyName": p.GetName(), "ParameterName": p.GetName() } + Systematic["ParameterBounds"] = [-3.0, 3.0] + Systematic["ParameterGroup"] = "Flux" + Systematic["ParameterValues"] = { "Generated": 0, "PreFitValue": 0 } + Systematic["StepScale"] = { "MCMC": 1 } + Systematic["Type"] = "Functional" + + if p.GetName() == "TargetUpstreamDegredation": + Systematic["ParameterBounds"] = [0, 1.0] + + Systematics.append({"Systematic": Systematic}) + +print(yaml.dump({"Systematics": Systematics})) \ No newline at end of file diff --git a/Systematics/Flux/flux_variations_FD_and_PRISM_2023.root b/Systematics/Flux/flux_variations_FD_and_PRISM_2023.root new file mode 100644 index 00000000..e037e05b Binary files /dev/null and b/Systematics/Flux/flux_variations_FD_and_PRISM_2023.root differ diff --git a/Systematics/Flux/plot_flux_sigvar.py b/Systematics/Flux/plot_flux_sigvar.py new file mode 100644 index 00000000..4bf8bd54 --- /dev/null +++ b/Systematics/Flux/plot_flux_sigvar.py @@ -0,0 +1,117 @@ +import uproot +import sys + +import matplotlib.pyplot as plt +from matplotlib.backends.backend_pdf import PdfPages + +import numpy as np + +Systematics = [] + +with uproot.open("flux_variations_FD_and_PRISM_2023.root") as inpfile: + with uproot.open("flux_1sigma_weights.root") as valifile: + with uproot.open(sys.argv[1]) as varfile: + + for pgk in inpfile["FluxParameters"].keys(recursive=False): + pg = pgk.split(";")[0] + for pk in inpfile["FluxParameters"][pg].keys(recursive=False): + p = pk.split(";")[0] + pdir = inpfile["FluxParameters"][pg][p] + offaxis_axis = pdir["OffAxisTAxis"] \ + if (pg == "Focussing") else \ + pdir["ND_nu_numu"]["OffAxisTAxis"] + + num_hists = np.count_nonzero([ x for x in pdir["ND_nu_numu"].keys(recursive=False) if (x.find("ND_nu_numu") == 0) ]) + hists = [ pdir["ND_nu_numu"][f"ND_nu_numu_{x}"] for x in range(num_hists) ] + Systematics.append({ "name": p, + "offaxis_axis": offaxis_axis, + "hists": hists}) + + + Systematics.reverse() + pp = PdfPages('FluxVariationsValidation.pdf') + for syst in Systematics: + print(syst["name"]) + + vardir_0m = varfile[syst["name"]]["ND_FHC_CCnumu_0m"] + vardir_12m = varfile[syst["name"]]["ND_FHC_CCnumu_12m"] + + fig, axs = plt.subplots(ncols=3, nrows=2, figsize=(18, 12)) + + fig.text(0.5,0.95,"Varied Parameter: " + syst["name"], horizontalalignment="center", size="xx-large") + + labels = [r"$-3\sigma$",r"$-1\sigma$",r"Nominal",r"$+1\sigma$",r"$+3\sigma$"] + chweel = ["firebrick", "tomato", "black", "deepskyblue", "mediumblue"] + lswheel = ["solid", "dashed", "dotted", "dashdot"] + for i in range(5): + bin_heights, bin_edges = vardir_0m[f"Variation_{i}"].to_numpy() + axs[0,0].stairs(bin_heights, bin_edges, label=labels[i], color=chweel[i]) + + if i != 2: + bin_heights_ref, bin_edges_ref = vardir_0m["Variation_2"].to_numpy() + axs[1,0].stairs(bin_heights/bin_heights_ref, bin_edges, color=chweel[i]) + + bin_heights, bin_edges = vardir_12m[f"Variation_{i}"].to_numpy() + axs[0,1].stairs(bin_heights, bin_edges, label=labels[i], color=chweel[i]) + + if i != 2: + bin_heights_ref, bin_edges_ref = vardir_12m["Variation_2"].to_numpy() + axs[1,1].stairs(bin_heights/bin_heights_ref, bin_edges, color=chweel[i]) + + for i,oad in enumerate([0, 12]): + off_axis_edges = syst["offaxis_axis"].edges() + oa_bin = np.argmax(off_axis_edges >= oad) + bin_heights, bin_edges = syst["hists"][oa_bin].to_numpy() + enu_lim = np.argmax(bin_edges > 5) + axs[0,2].stairs(bin_heights[:enu_lim-1], bin_edges[:enu_lim], color=chweel[3], \ + linestyle=lswheel[i], label=f"${off_axis_edges[oa_bin]} -- {off_axis_edges[oa_bin+1]}$ m") + + bin_heights, bin_edges = valifile[syst["name"]][f"ND_Numu_{oad}m"].to_numpy() + axs[1,2].stairs(bin_heights, bin_edges, color=chweel[3], \ + linestyle=lswheel[i], label=f"{oad} m") + + axs[1,0].set_ylim([0.8,1.2]) + + axs[0,2].set_ylim([-0.2, 0.2]) + axs[1,1].set_ylim([0.8,1.2]) + axs[1,2].set_ylim([0.8,1.2]) + + axs[0,0].set_xlim([0,5]) + axs[1,0].set_xlim([0,5]) + + axs[0,1].set_xlim([0,5]) + axs[1,1].set_xlim([0,5]) + + axs[0,2].set_xlim([0,5]) + axs[1,2].set_xlim([0,5]) + + fig.text(0.2,0.89,"MaCh3 SigVar On Axis",horizontalalignment="left", size="large") + fig.text(0.2,0.47,"MaCh3 SigVar On Axis Ratio",horizontalalignment="left", size="large") + + fig.text(0.4,0.89,"MaCh3 SigVar 12 m Off Axis",horizontalalignment="left", size="large") + fig.text(0.4,0.47,"MaCh3 SigVar 12 m Off Axis Ratio",horizontalalignment="left", size="large") + + fig.text(0.65,0.89,"Histograms from inputs",horizontalalignment="left", size="large") + fig.text(0.65,0.47,"Weights directly from Weight Interface",horizontalalignment="left", size="large") + + axs[0,0].legend() + axs[0,2].legend() + axs[1,2].legend() + + axs[0,0].set_ylabel("Event Rate", size="x-large") + axs[1,0].set_ylabel("Varied/Nominal Event Rate", size="x-large") + axs[0,1].set_ylabel("Event Rate", size="x-large") + axs[1,1].set_ylabel("Varied/Nominal Event Rate", size="x-large") + + axs[0,2].set_ylabel("Input Ratio", size="x-large") + axs[1,2].set_ylabel("Systematic Weight", size="x-large") + + axs[0,0].set_xlabel(r"$E_\nu$", size="x-large") + axs[1,0].set_xlabel(r"$E_\nu$", size="x-large") + axs[0,1].set_xlabel(r"$E_\nu$", size="x-large") + axs[1,1].set_xlabel(r"$E_\nu$", size="x-large") + axs[0,2].set_xlabel(r"$E_\nu$", size="x-large") + axs[1,2].set_xlabel(r"$E_\nu$", size="x-large") + pp.savefig(fig) + plt.close(fig) + pp.close() \ No newline at end of file diff --git a/Systematics/Flux/test_read_flux_systs.C b/Systematics/Flux/test_read_flux_systs.C new file mode 100644 index 00000000..8ebb2577 --- /dev/null +++ b/Systematics/Flux/test_read_flux_systs.C @@ -0,0 +1,76 @@ +{ + gROOT->ProcessLine(".include ../../"); + gROOT->ProcessLine(".L OffAxisFluxUncertaintyHelper.cxx+"); + + gROOT->ProcessLine(R"( + OffAxisFluxUncertaintyHelper flux_helper; + flux_helper.Initialize("flux_variations_FD_and_PRISM_2023.root", true); + + + for(size_t i = 0; i < flux_helper.GetNFocussingParams(); ++i){ + std::cout << i << ": " << flux_helper.GetFocussingParamName(i) + << std::endl; + + std::cout << " Binning: [ "; + for(auto be : flux_helper.GetFluxFocussingOffAxisBinning()){ + std::cout << be << ", "; + } + std::cout << "]" << std::endl; + } + + std::cout << "Have " << flux_helper.GetNHadProdPCAComponents() + << " hadron production PCA components." << std::endl; + + for(int nucfg = 0; nucfg < OffAxisFluxUncertaintyHelper::kFD_numu_numode; ++nucfg){ + std::cout << "HadProd Off Axis Binning for nucfg: " << nucfg << ", [ "; + for(auto be : flux_helper.GetFluxHadProdOffAxisBinning(nucfg)){ + std::cout << be << ", "; + } + std::cout << "]" << std::endl; + } + + TFile fout("flux_1sigma_weights.root","RECREATE"); + + for(size_t i = 0; i < flux_helper.GetNFocussingParams(); ++i){ + auto dir = fout.mkdir(flux_helper.GetFocussingParamName(i).c_str()); + + for(double oa_m : std::vector{0.0, 12.0}){ + TH1D hist(("ND_Numu_" + std::to_string(int(oa_m)) + "m").c_str(), + (flux_helper.GetFocussingParamName(i) +";Enu;Weight").c_str(), + 1000,0,5); + + int nucfg = 0; + for(int bi = 0; bi < 1000; ++bi){ + int syst_bi = flux_helper.GetFocussingBin(nucfg, + hist.GetXaxis()->GetBinCenter(bi+1), oa_m); + hist.SetBinContent(bi+1, flux_helper.GetFluxFocussingWeight(i, 1, + nucfg, syst_bi)); + } + dir->WriteTObject(&hist, hist.GetName()); + } + } + + for(size_t i = 0; i < flux_helper.GetNHadProdPCAComponents(); ++i){ + auto dir = fout.mkdir(("HadronProduction_pca_" + std::to_string(i)).c_str()); + + for(double oa_m : std::vector{0.0, 12.0}){ + TH1D hist(("ND_Numu_" + std::to_string(int(oa_m)) + "m").c_str(), + ("HadronProduction_pca_" + std::to_string(i)).c_str(), + 1000,0,5); + + int nucfg = 0; + for(int bi = 0; bi < 1000; ++bi){ + int syst_bi = flux_helper.GetHadProdBin(nucfg, + hist.GetXaxis()->GetBinCenter(bi+1), oa_m); + hist.SetBinContent(bi+1, flux_helper.GetFluxHadProdWeight(i, 1, + nucfg, syst_bi)); + } + dir->WriteTObject(&hist, hist.GetName()); + } + } + + fout.Write(); + fout.Close(); + + )"); +} \ No newline at end of file