diff --git a/app/CMakeLists.txt b/app/CMakeLists.txt
index 1d488e3..1afec7d 100644
--- a/app/CMakeLists.txt
+++ b/app/CMakeLists.txt
@@ -1,217 +1,225 @@
# Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
################################################################################
# This file is part of NUISANCE.
#
# NUISANCE is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# NUISANCE is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with NUISANCE. If not, see .
################################################################################
set(TARGETS_TO_BUILD)
if(USE_MINIMIZER)
add_executable(nuismin nuismin.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuismin)
target_link_libraries(nuismin ${MODULETargets})
target_link_libraries(nuismin ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(nuismin ${ROOT_LIBS})
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(nuismin PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
add_executable(nuissplines nuissplines.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuissplines)
target_link_libraries(nuissplines ${MODULETargets})
target_link_libraries(nuissplines ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(nuissplines ${ROOT_LIBS})
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(nuissplines PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
endif()
include_directories(${RWENGINE_INCLUDE_DIRECTORIES})
include_directories(${CMAKE_SOURCE_DIR}/src/Routines)
include_directories(${CMAKE_SOURCE_DIR}/src/InputHandler)
include_directories(${CMAKE_SOURCE_DIR}/src/Genie)
include_directories(${CMAKE_SOURCE_DIR}/src/FitBase)
include_directories(${CMAKE_SOURCE_DIR}/src/Statistical)
include_directories(${CMAKE_SOURCE_DIR}/src/Utils)
include_directories(${CMAKE_SOURCE_DIR}/src/Config)
include_directories(${CMAKE_SOURCE_DIR}/src/Logger)
include_directories(${CMAKE_SOURCE_DIR}/src/Splines)
include_directories(${CMAKE_SOURCE_DIR}/src/Reweight)
include_directories(${CMAKE_SOURCE_DIR}/src/FCN)
include_directories(${CMAKE_SOURCE_DIR}/src/MCStudies)
include_directories(${CMAKE_SOURCE_DIR}/src/Smearceptance)
include_directories(${EXP_INCLUDE_DIRECTORIES})
if (USE_NuWro AND NOT NUWRO_BUILT_FROM_FILE)
add_executable(nuwro_nuisance nuwro_NUISANCE.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuwro_nuisance)
target_link_libraries(nuwro_nuisance ${MODULETargets})
target_link_libraries(nuwro_nuisance ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(nuwro_nuisance ${ROOT_LIBS})
include_directories(${CMAKE_SOURCE_DIR}/src/FitBase)
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(nuwro_nuisance PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
endif()
if (USE_NEUT)
add_executable(neut_nuisance neut_NUISANCE.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};neut_nuisance)
target_link_libraries(neut_nuisance ${MODULETargets})
target_link_libraries(neut_nuisance ${CMAKE_DEPENDLIB_FLAGS})
target_link_libraries(neut_nuisance ${ROOT_LIBS})
include_directories(${CMAKE_SOURCE_DIR}/src/FitBase)
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(neut_nuisance PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
endif()
if (BUILD_GEVGEN)
add_executable(gevgen_nuisance gEvGen_NUISANCE.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};gevgen_nuisance)
target_link_libraries(gevgen_nuisance ${MODULETargets})
target_link_libraries(gevgen_nuisance ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(gevgen_nuisance ${ROOT_LIBS})
include_directories(${CMAKE_SOURCE_DIR}/src/FitBase)
include_directories(${GENIE_INCLUDES}/Apps)
include_directories(${GENIE_INCLUDES}/FluxDrivers)
include_directories(${GENIE_INCLUDES}/EVGDrivers)
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(gevgen_nuisance PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
add_executable(gevgen_nuisance_mixed gEvGen_NUISANCE_MIXED.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};gevgen_nuisance_mixed)
target_link_libraries(gevgen_nuisance_mixed ${MODULETargets})
target_link_libraries(gevgen_nuisance_mixed ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(gevgen_nuisance_mixed ${ROOT_LIBS})
include_directories(${CMAKE_SOURCE_DIR}/src/FitBase)
include_directories(${GENIE_INCLUDES}/Apps)
include_directories(${GENIE_INCLUDES}/FluxDrivers)
include_directories(${GENIE_INCLUDES}/EVGDrivers)
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(gevgen_nuisance_mixed PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
endif()
if (USE_GiBUU)
add_executable(DumpGiBUUEvents DumpGiBUUEvents.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};DumpGiBUUEvents)
target_link_libraries(DumpGiBUUEvents ${MODULETargets})
target_link_libraries(DumpGiBUUEvents ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(DumpGiBUUEvents ${ROOT_LIBS})
include_directories(${CMAKE_SOURCE_DIR}/src/FitBase)
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(DumpGiBUUEvents PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
endif()
add_executable(nuiscomp nuiscomp.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuiscomp)
target_link_libraries(nuiscomp ${MODULETargets})
target_link_libraries(nuiscomp ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(nuiscomp ${ROOT_LIBS})
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(nuiscomp PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
add_executable(nuisflat nuisflat.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuisflat)
target_link_libraries(nuisflat ${MODULETargets})
target_link_libraries(nuisflat ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(nuisflat ${ROOT_LIBS})
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(nuisflat PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
add_executable(nuissmear nuissmear.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuissmear)
target_link_libraries(nuissmear ${MODULETargets})
target_link_libraries(nuissmear ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(nuissmear ${ROOT_LIBS})
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(nuissmear PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
add_executable(nuissyst nuissyst.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuissyst)
target_link_libraries(nuissyst ${MODULETargets})
target_link_libraries(nuissyst ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(nuissyst ${ROOT_LIBS})
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(nuissyst PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
add_executable(nuisbayes nuisbayes.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuisbayes)
target_link_libraries(nuisbayes ${MODULETargets})
target_link_libraries(nuisbayes ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(nuisbayes ${ROOT_LIBS})
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(nuisbayes PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
if(USE_GENIE)
add_executable(PrepareGENIE PrepareGENIE.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};PrepareGENIE)
target_link_libraries(PrepareGENIE ${MODULETargets})
target_link_libraries(PrepareGENIE ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(PrepareGENIE ${ROOT_LIBS})
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(PrepareGENIE PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
endif()
if(USE_NEUT)
add_executable(PrepareNEUT PrepareNEUT.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};PrepareNEUT)
target_link_libraries(PrepareNEUT ${MODULETargets})
target_link_libraries(PrepareNEUT ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(PrepareNEUT ${ROOT_LIBS})
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(PrepareNEUT PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
endif()
# PREPARE NUWRO
# Commented out for the time being until it is finished..
if(USE_NuWro)
add_executable(PrepareNuwro PrepareNuwroEvents.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};PrepareNuwro)
target_link_libraries(PrepareNuwro ${MODULETargets})
target_link_libraries(PrepareNuwro ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(PrepareNuwro ${ROOT_LIBS})
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(PrepareNuwro PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
endif()
+add_executable(nuisbac nuisbac.cxx)
+set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuisbac)
+target_link_libraries(nuisbac ${MODULETargets})
+target_link_libraries(nuisbac ${CMAKE_DEPENDLIB_FLAGS})
+if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
+ set_target_properties(nuisbac PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
+endif()
+
install(TARGETS ${TARGETS_TO_BUILD} DESTINATION bin)
#add_executable(DumpROOTClassesFromVector DumpROOTClassesFromVector.cxx)
# #Strip out -lNuWro_event1
# string(REPLACE "-lNuWro_event1" "" NWEVSTRIPPED_CDF ${CMAKE_DEPENDLIB_FLAGS})
# cmessage(DEBUG "Attempted to strip out nuwro library: \"${CMAKE_DEPENDLIB_FLAGS}\" -> \"${NWEVSTRIPPED_CDF}\"")
# add_executable(PrepareNEUT PrepareNEUT.cxx)
# target_link_libraries(DumpROOTClassesFromVector ${MODULETargets})
# target_link_libraries(DumpROOTClassesFromVector ${NWEVSTRIPPED_CDF})
# if(NOT CMAKE_LINK_FLAGS STREQUAL "")
# set_target_properties(DumpROOTClassesFromVector PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
# endif()
#install(TARGETS DumpROOTClassesFromVector DESTINATION bin)
diff --git a/app/nuisbac.cxx b/app/nuisbac.cxx
new file mode 100644
index 0000000..0b6df2f
--- /dev/null
+++ b/app/nuisbac.cxx
@@ -0,0 +1,9 @@
+#include "Initialiser.h"
+
+int main(){
+
+ for(size_t i = 0; i < 100; ++i){
+ nuisance_init();
+ }
+
+}
diff --git a/cmake/ROOTSetup.cmake b/cmake/ROOTSetup.cmake
index ebf9dee..37f6157 100644
--- a/cmake/ROOTSetup.cmake
+++ b/cmake/ROOTSetup.cmake
@@ -1,162 +1,167 @@
# Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
################################################################################
# This file is part of NUISANCE.
#
# NUISANCE is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# NUISANCE is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with NUISANCE. If not, see .
################################################################################
if ( NOT DEFINED ENV{ROOTSYS} )
cmessage (FATAL_ERROR "$ROOTSYS is not defined, please set up ROOT first.")
else()
cmessage(STATUS "Using ROOT installed at $ENV{ROOTSYS}")
set(CMAKE_ROOTSYS $ENV{ROOTSYS})
endif()
# Get cflags from ROOT
execute_process (COMMAND root-config
--cflags OUTPUT_VARIABLE ROOT_CXX_FLAGS_RAW OUTPUT_STRIP_TRAILING_WHITESPACE)
string(REPLACE " " ";" ROOT_CXX_FLAGS "${ROOT_CXX_FLAGS_RAW}")
# Get libdir from ROOT
execute_process (COMMAND root-config
--libdir OUTPUT_VARIABLE ROOT_LIBDIR OUTPUT_STRIP_TRAILING_WHITESPACE)
# Get version from ROOT
execute_process (COMMAND root-config
--version OUTPUT_VARIABLE ROOT_VERSION OUTPUT_STRIP_TRAILING_WHITESPACE)
# Get features from ROOT
execute_process (COMMAND root-config
--features OUTPUT_VARIABLE ROOT_FEATURES OUTPUT_STRIP_TRAILING_WHITESPACE)
LIST(APPEND EXTRA_LINK_DIRS ${ROOT_LIBDIR})
LIST(APPEND ROOT_LIBS
Core
Cint
RIO
XMLIO
Net
Hist
Graf
Graf3d
Gpad
Tree
Rint
Postscript
Matrix
Physics
MathCore
Thread
EG
Geom
GenVector)
if(USE_MINIMIZER)
if("${ROOT_FEATURES}" MATCHES "minuit2")
cmessage(STATUS "ROOT built with MINUIT2 support")
LIST(APPEND EXTRA_CXX_FLAGS -D__MINUIT2_ENABLED__)
else()
cmessage(FATAL_ERROR "ROOT built without MINUIT2 support but minimizer functionality requested. Either configure with -DUSE_MINIMIZER=FALSE or use a version of ROOT with MINUIT2 support.")
endif()
- string(REGEX MATCH "5.34/([0-9]+)" ROOTVERSMATCH ${ROOT_VERSION})
- if(NOT ROOTVERSMATCH OR ${CMAKE_MATCH_1} LESS "19")
- cmessage(FATAL_ERROR "ROOT Version: ${ROOT_VERSION} has out of date minimizer interface, but minimizer functionality requested. Please configure with -DUSE_MINIMIZER=FALSE or update to 5.34/19 or greater to enable minimization features.")
+ string(REGEX MATCH "6\..*+)" ROOTVERSIXMATCH ${ROOT_VERSION})
+ if(ROOTVERSIXMATCH)
+ cmessage(STATUS "Using ROOT6, We are essentiall flying blind here.")
+ else()
+ string(REGEX MATCH "5.34/([0-9]+)" ROOTVERSMATCH ${ROOT_VERSION})
+ if(NOT ROOTVERSMATCH OR ${CMAKE_MATCH_1} LESS "19")
+ cmessage(FATAL_ERROR "ROOT Version: ${ROOT_VERSION} has out of date minimizer interface, but minimizer functionality requested. Please configure with -DUSE_MINIMIZER=FALSE or update to 5.34/19 or greater to enable minimization features.")
+ endif()
endif()
endif()
if("${ROOT_FEATURES}" MATCHES "opengl")
cmessage(STATUS "ROOT built with OpenGL support")
LIST(APPEND ROOT_LIBS RGL)
endif()
if(DEFINED NEED_ROOTPYTHIA6 AND NEED_ROOTPYTHIA6)
LIST(APPEND ROOT_LIBS EGPythia6 Pythia6)
endif()
cmessage ( STATUS "[ROOT]: root-config --version: ${ROOT_VERSION} ")
cmessage ( STATUS "[ROOT]: root-config --cflags : ${ROOT_CXX_FLAGS} ")
cmessage ( STATUS "[ROOT]: root-config --libs : ${ROOT_LD_FLAGS} ")
LIST(APPEND EXTRA_CXX_FLAGS ${ROOT_CXX_FLAGS})
#Helper functions for building dictionaries
function(GenROOTDictionary OutputDictName Header LinkDef)
get_directory_property(incdirs INCLUDE_DIRECTORIES)
string(REPLACE ";" ";-I" LISTDIRINCLUDES "-I${incdirs}")
string(REPLACE " " ";" LISTCPPFLAGS "${EXTRA_CXX_FLAGS}")
#ROOT5 CINT cannot handle it.
list(REMOVE_ITEM LISTCPPFLAGS "-std=c++11")
message(STATUS "LISTCPPFLAGS: ${LISTCPPFLAGS}")
message(STATUS "LISTINCLUDES: ${LISTDIRINCLUDES}")
#Learn how to generate the Dict.cxx and Dict.hxx
add_custom_command(
OUTPUT "${OutputDictName}.cxx" "${OutputDictName}.h"
COMMAND rootcint
ARGS -f ${OutputDictName}.cxx -c
-p ${LISTDIRINCLUDES} ${LISTCPPFLAGS} ${Header} ${LinkDef}
DEPENDS ${Header};${LinkDef})
endfunction()
function(BuildROOTProject ProjectName InputFile CommaSeparatedClassesToDump LIBLINKMODE)
string(REPLACE "," ";" HeadersToDump ${CommaSeparatedClassesToDump})
set(OUTPUTFILES ${CMAKE_BINARY_DIR}/${ProjectName}/${ProjectName}ProjectSource.cxx
${CMAKE_BINARY_DIR}/${ProjectName}/${ProjectName}LinkDef.h
${CMAKE_BINARY_DIR}/${ProjectName}/${ProjectName}ProjectHeaders.h
${CMAKE_BINARY_DIR}/${ProjectName}/${ProjectName}ProjectInstances.h)
cmessage(STATUS "As part of ROOT project: ${ProjectName}")
foreach (header ${HeadersToDump})
LIST(APPEND OUTPUTFILES "${CMAKE_BINARY_DIR}/${ProjectName}/${header}.h")
cmessage(STATUS "Will generate: ${CMAKE_BINARY_DIR}/${ProjectName}/${header}.h")
endforeach()
add_custom_command(
OUTPUT ${OUTPUTFILES}
COMMAND ${CMAKE_BINARY_DIR}/src/Utils/DumpROOTClassesFromVector
ARGS ${InputFile}
${CMAKE_BINARY_DIR}/${ProjectName}
${CommaSeparatedClassesToDump}
VERBATIM
DEPENDS DumpROOTClassesFromVector)
add_custom_target(${ProjectName}_sources
DEPENDS ${OUTPUTFILES})
GenROOTDictionary(
${CMAKE_BINARY_DIR}/${ProjectName}/${ProjectName}ProjectDict
${CMAKE_BINARY_DIR}/${ProjectName}/${ProjectName}ProjectHeaders.h
${CMAKE_BINARY_DIR}/${ProjectName}/${ProjectName}LinkDef.h
)
add_custom_target(${ProjectName}ProjectDict
DEPENDS
${CMAKE_BINARY_DIR}/${ProjectName}/${ProjectName}ProjectDict.cxx
${CMAKE_BINARY_DIR}/${ProjectName}/${ProjectName}ProjectDict.h )
# add_dependencies(${ProjectName}ProjectDict ${ProjectName}_sources)
#ProjectSource.cxx includes ProjectDict.cxx, so no need to add to compilation.
set(ROAA_SOURCEFILES
${CMAKE_BINARY_DIR}/${ProjectName}/${ProjectName}ProjectSource.cxx)
add_library(${ProjectName} ${LIBLINKMODE} ${ROAA_SOURCEFILES})
add_dependencies(${ProjectName} ${ProjectName}ProjectDict)
endfunction()
diff --git a/parameters/dial_conversion.card b/parameters/dial_conversion.card
index b34a7d9..0227a07 100644
--- a/parameters/dial_conversion.card
+++ b/parameters/dial_conversion.card
@@ -1,23 +1,38 @@
# par Name Units Nominal FracErr ConvFunc
neut_parameter MaCCQE GeV 1.21*(1.0+x*0.16)
+neut_parameter MaNFFRES GeV 0.95*(1.0+x*0.157894737)
+neut_parameter CA5RES x 1.01*(1.0+x*0.247524752)
+neut_parameter BgSclRES x 1.30*(1.0+x*0.153846154)
+
+neut_parameter FrAbs_pi x 1.1*(1.0+x*0.5)
+neut_parameter FrInelLow_pi x 1*(1.0+x*0.5)
+neut_parameter FrInelHigh_pi x 1.8*(1.0+x*0.3)
+neut_parameter FrPiProd_pi x 1*(1.0+x*0.5)
+neut_parameter FrCExLow_pi x 1*(1.0+x*0.5)
+neut_parameter FrCExHigh_pi x 1.8*(1.0+x*0.3)
+
niwg_parameter NIWGMEC_Norm_C12 % 100.0*(1.0+x)
niwg_parameter VecFFCCQE MDLQE x
niwg_parameter CCQEFermiSurfMom MeV 217*(1.0+x*0.15)
-
t2k_parameter NIWG2014a_pF_C12 MeV 217*(1.0+x*0.073733)
t2k_parameter NIWG2014a_pF_O16 MeV 225*(1.0+x*0.071111)
t2k_parameter NIWGMEC_PDDWeight_C12 x 0.5*(1.0+x*1.0)
t2k_parameter NIWGMEC_PDDWeight_O16 x 0.5*(1.0+x*1.0)
t2k_parameter NXSec_MaCCQE GeV 1.21*(1.0+x*0.165289256)
t2k_parameter NXSec_MaNFFRES GeV 0.95*(1.0+x*0.157894737)
t2k_parameter NXSec_CA5RES x 1.01*(1.0+x*0.247524752)
t2k_parameter NXSec_BgSclCCRES x 1.30*(1.0+x*0.153846154)
t2k_parameter NIWG_Effective_rpaCCQE_A x 1.0*(1.0 + x*0.1)
t2k_parameter NIWG_Effective_rpaCCQE_B x 1.0*(1.0 + x*0.1)
t2k_parameter NIWG_Effective_rpaCCQE_C x 1.0*(1.0 + x*0.1)
t2k_parameter NIWG_Effective_rpaCCQE_D x 1.0*(1.0 + x*0.1)
t2k_parameter NIWG_Effective_rpaCCQE_U x 1.0*(1.0 + x*0.1)
+
+nuwro_parameter kNuwro_Ma_CCQE MeV 1200*(1.0+x*0.160)
+nuwro_parameter kNuwro_MaRES GeV 0.940*(1.0+x*0.1)
+nuwro_parameter kNuwro_CA5 x 1.19*(1.0+x*0.1)
+nuwro_parameter kNuwro_SPPBkgScale x 1.0*(1.0+x*0.26)
diff --git a/src/Electron/CLAS6-EG2_Accepter.cxx b/src/Electron/CLAS6-EG2_Accepter.cxx
index d2c4307..3b910ee 100644
--- a/src/Electron/CLAS6-EG2_Accepter.cxx
+++ b/src/Electron/CLAS6-EG2_Accepter.cxx
@@ -1,177 +1,252 @@
#include "ISmearcepter.h"
#include "TH3D.h"
#include "TRandom3.h"
#include
-#define DEBUG_CLASACCEPT
+// #define DEBUG_CLASACCEPT 1
+
+struct EffMap {
+ TH3D *Generated;
+ TH3D *Accepted;
+
+ void Build(TFile *inpF, std::string const &GenName,
+ std::string const &AccName) {
+ Generated = dynamic_cast(inpF->Get(GenName.c_str()));
+ Accepted = dynamic_cast(inpF->Get(AccName.c_str()));
+
+ if (!Generated) {
+ std::cout << "[ERROR]: Could not retrieve \"accepted\" histogram: \""
+ << AccName << "\" from file: \"" << inpF->GetName() << "\"."
+ << std::endl;
+ exit(1);
+ }
+ if (!Accepted) {
+ std::cout << "[ERROR]: Could not retrieve \"generated\" histogram: \""
+ << GenName << "\" from file: \"" << inpF->GetName() << "\"."
+ << std::endl;
+ exit(1);
+ }
+
+ Generated = static_cast(Generated->Clone());
+ Generated->SetDirectory(NULL);
+ Accepted = static_cast(Accepted->Clone());
+ Accepted->SetDirectory(NULL);
+ }
+
+ double GetAccRatio(double p_GeV, double costheta, double phi_deg) {
+ // For a bin in phase space defined by p, cost, phi:
+ // Find number of generated events
+ Int_t pbin = Generated->GetXaxis()->FindBin(p_GeV);
+
+ Int_t tbin = Generated->GetYaxis()->FindBin(costheta);
+ Int_t phibin = Generated->GetZaxis()->FindBin(phi_deg);
+
+ if (((pbin == 0) || (pbin == (Generated->GetXaxis()->GetNbins() + 1))) ||
+ ((tbin == 0) || (tbin == (Generated->GetYaxis()->GetNbins() + 1))) ||
+ ((phibin == 0) ||
+ (phibin == (Generated->GetZaxis()->GetNbins() + 1)))) {
+ return 0;
+ }
+
+ double num_gen = Generated->GetBinContent(pbin, tbin, phibin);
+ // Find number of accepted events
+ pbin = Accepted->GetXaxis()->FindBin(p_GeV);
+ tbin = Accepted->GetYaxis()->FindBin(costheta);
+ phibin = Accepted->GetZaxis()->FindBin(phi_deg);
+ double num_acc = Accepted->GetBinContent(pbin, tbin, phibin);
+ double acc_ratio = double(num_acc) / double(num_gen);
+
+ if (((pbin == 0) || (pbin == (Accepted->GetXaxis()->GetNbins() + 1))) ||
+ ((tbin == 0) || (tbin == (Accepted->GetYaxis()->GetNbins() + 1))) ||
+ ((phibin == 0) || (phibin == (Accepted->GetZaxis()->GetNbins() + 1)))) {
+ return 0;
+ }
+
+ if ((acc_ratio != 0 && !std::isnormal(acc_ratio)) || (acc_ratio > 1)) {
+ std::cout << "[BINS]: p " << Generated->GetXaxis()->GetBinLowEdge(1)
+ << " -- "
+ << Generated->GetXaxis()->GetBinUpEdge(
+ Generated->GetXaxis()->GetNbins())
+ << ", cost " << Generated->GetYaxis()->GetBinLowEdge(1)
+ << " -- "
+ << Generated->GetYaxis()->GetBinUpEdge(
+ Generated->GetYaxis()->GetNbins())
+ << ", phi " << Generated->GetZaxis()->GetBinLowEdge(1) << " -- "
+ << Generated->GetZaxis()->GetBinUpEdge(
+ Generated->GetZaxis()->GetNbins())
+ << ". " << std::endl
+ << "[ERROR]: Bad acceptance ratio: " << acc_ratio << " = "
+ << num_acc << " / " << num_gen << ". (" << p_GeV << ", "
+ << costheta << ", " << phi_deg << ")." << std::endl;
+ exit(1);
+ }
+ return acc_ratio;
+ }
+};
class CLASAccepter : public ISmearcepter {
TRandom3 rand;
// Maps a particle PDG to the relevant generated and accepted histograms from
// the input map.
- std::map > Acceptance;
+ std::map Acceptance;
public:
CLASAccepter() { ElementName = "CLASAccepter"; }
void SpecifcSetup(nuiskey &nk) {
rand.~TRandom3();
new (&rand) TRandom3();
InstanceName = nk.GetS("name");
std::string const &mapfile = nk.GetS("map");
if (!mapfile.length()) {
std::cout << "[ERROR]: No input file specified by \"map\" attribute."
<< std::endl;
exit(1);
}
TFile *f = new TFile(mapfile.c_str());
if (!f || !f->IsOpen()) {
std::cout << "[ERROR]: Could not open root file specified by \"map\" "
"attribute: \""
<< mapfile << "\"" << std::endl;
exit(1);
}
std::vector accepts = nk.GetListOfChildNodes("accept");
for (auto &acc : accepts) {
std::string const &genStr = acc.GetS("generated");
std::string const &accStr = acc.GetS("accepted");
if (!genStr.length() || !accStr.length()) {
std::cout << "[ERROR]: expected accept node to contain both "
"\"generated\" and \"accepted\" attributes."
<< std::endl;
exit(1);
}
std::string const &pdgs_s = acc.GetS("PDG");
std::vector pdgs_i = GeneralUtils::ParseToInt(pdgs_s, ",");
if (!pdgs_i.size()) {
std::cout
<< "[ERROR]: Could not find any applicable particle PDG codes."
<< std::endl;
exit(1);
}
- std::pair genacc;
-
- genacc.first = dynamic_cast(f->Get(accStr.c_str()));
- genacc.second = dynamic_cast(f->Get(genStr.c_str()));
-
- if (!genacc.first) {
- std::cout << "[ERROR]: Could not retrieve \"accepted\" histogram: \""
- << accStr << "\" from file: \"" << mapfile << "\"."
- << std::endl;
- exit(1);
- }
- if (!genacc.second) {
- std::cout << "[ERROR]: Could not retrieve \"accepted\" histogram: \""
- << genStr << "\" from file: \"" << mapfile << "\"."
- << std::endl;
- exit(1);
- }
+ EffMap ef;
- genacc.first = static_cast(genacc.first->Clone());
- genacc.first->SetDirectory(NULL);
- genacc.second = static_cast(genacc.second->Clone());
- genacc.second->SetDirectory(NULL);
+ ef.Build(f, genStr, accStr);
for (size_t pdg_it = 0; pdg_it < pdgs_i.size(); ++pdg_it) {
if (Acceptance.count(pdgs_i[pdg_it])) {
std::cout
<< "[WARN]: Acceptance map already contains acceptance for PDG: "
<< pdgs_i[pdg_it] << ". Overwriting..." << std::endl;
}
- Acceptance[pdgs_i[pdg_it]] = genacc;
+ Acceptance[pdgs_i[pdg_it]] = ef;
}
}
std::cout << "Loaded " << Acceptance.size()
<< " particle acceptance definitions." << std::endl;
f->Close();
delete f;
}
RecoInfo *Smearcept(FitEvent *fe) {
RecoInfo *ri = new RecoInfo();
for (size_t p_it = 0; p_it < fe->NParticles(); ++p_it) {
FitParticle *fp = fe->GetParticle(p_it);
int PDG = fp->PDG();
- double p = fp->P3().Mag();
+ double p = fp->P3().Mag() * 1E-3;
double cost = fp->P3().CosTheta();
- double phi = fp->P3().Phi();
+ double phi = fp->P3().Phi() * (180.0 / TMath::Pi()) + 150.0;
-#ifdef DEBUG_CLASACCEPT
+#if DEBUG_CLASACCEPT > 1
std::cout << std::endl;
std::cout << "[" << p_it << "]: " << PDG << ", " << fp->Status() << ", "
<< fp->E() << " -- KE:" << fp->KE() << " Mom: " << p
<< std::flush;
#endif
if (fp->Status() != kFinalState) {
-#ifdef DEBUG_CLASACCEPT
+#if DEBUG_CLASACCEPT > 1
std::cout << " -- Not final state." << std::flush;
#endif
continue;
}
if (!Acceptance.count(PDG)) {
-#ifdef DEBUG_CLASACCEPT
+#if DEBUG_CLASACCEPT > 1
std::cout << " -- Unknown acceptance." << std::flush;
#endif
continue;
}
- std::pair acc = Acceptance[PDG];
-
- // For a bin in phase space defined by p, cost, phi:
- // Find number of generated events
- Int_t pbin = acc.second->GetXaxis()->FindBin(p);
- Int_t tbin = acc.second->GetYaxis()->FindBin(cost);
- Int_t phibin = acc.second->GetZaxis()->FindBin(phi);
- double num_gen = acc.second->GetBinContent(pbin, tbin, phibin);
- // Find number of accepted events
- pbin = acc.first->GetXaxis()->FindBin(p);
- tbin = acc.first->GetYaxis()->FindBin(cost);
- phibin = acc.first->GetZaxis()->FindBin(phi);
- double num_acc = acc.first->GetBinContent(pbin, tbin, phibin);
- double acc_ratio = double(num_acc) / double(num_gen);
+ EffMap eff = Acceptance[PDG];
+
+ double acc_ratio = eff.GetAccRatio(p, cost, phi);
bool accepted = (rand.Uniform() < acc_ratio);
if (accepted) {
#ifdef DEBUG_CLASACCEPT
+ std::cout << "(" << p << ", " << cost << ", " << phi << ")."
+ << std::endl;
std::cout << " -- Reconstructed with probability: " << acc_ratio
<< std::flush;
#endif
ri->RecObjMom.push_back(fp->P3());
ri->RecObjClass.push_back(fp->PDG());
continue;
}
#ifdef DEBUG_CLASACCEPT
+ std::cout << "(" << p << ", " << cost << ", " << phi << ")." << std::endl;
std::cout << " -- Rejected with probability: " << acc_ratio << std::flush;
#endif
-#ifdef DEBUG_CLASACCEPT
- std::cout << std::endl;
-#endif
}
#ifdef DEBUG_CLASACCEPT
- std::cout << "Reconstructed " << ri->RecObjMom.size() << " particles. "
- << std::endl;
+ std::cout << std::endl;
+
+ if (ri->RecObjMom.size()) {
+ std::cout << "Reconstructed " << ri->RecObjMom.size() << " particles. "
+ << std::endl;
+ }
#endif
return ri;
}
+
+ double GetEfficiency(FitEvent *fe) {
+ double effweight = 1;
+ for (size_t p_it = 0; p_it < fe->NParticles(); ++p_it) {
+ FitParticle *fp = fe->GetParticle(p_it);
+
+ int PDG = fp->PDG();
+ double p = fp->P3().Mag() * 1E-3;
+ double cost = fp->P3().CosTheta();
+ double phi = fp->P3().Phi() * (180.0 / TMath::Pi()) + 150.0;
+ if (fp->Status() != kFinalState) {
+ continue;
+ }
+ if (!Acceptance.count(PDG)) {
+ continue;
+ }
+ EffMap eff = Acceptance[PDG];
+
+ effweight *= eff.GetAccRatio(p, cost, phi);
+ }
+ return effweight;
+ }
};
diff --git a/src/FitBase/CustomVariableBoxes.h b/src/FitBase/CustomVariableBoxes.h
index 8bb8b41..e900871 100644
--- a/src/FitBase/CustomVariableBoxes.h
+++ b/src/FitBase/CustomVariableBoxes.h
@@ -1,21 +1,28 @@
#ifndef CUSTOM_VARIABLES_BOX_H
#define CUSTOM_VARIABLES_BOX_H
#include "MeasurementVariableBox.h"
#include "MeasurementVariableBox1D.h"
#include "MeasurementVariableBox2D.h"
/*!
* \addtogroup FitBase
* @{
*/
/// Custom box used to also save Q2 for each event.
class Q2VariableBox1D : public MeasurementVariableBox1D {
-public:
- inline Q2VariableBox1D() { Reset(); };
- inline void Reset() { fQ2 = -999.9; }
- double fQ2;
+ public:
+ inline Q2VariableBox1D() { Reset(); };
+ inline MeasurementVariableBox* CloneSignalBox() {
+ Q2VariableBox1D* box = new Q2VariableBox1D();
+ box->fX = this->fX;
+ box->fSampleWeight = this->fSampleWeight;
+ box->fQ2 = this->fQ2;
+ return box;
+ };
+ inline void Reset() { fQ2 = -999.9; }
+ double fQ2;
};
#endif
diff --git a/src/InputHandler/GIBUUInputHandler.cxx b/src/InputHandler/GIBUUInputHandler.cxx
index c0e8724..f61559f 100644
--- a/src/InputHandler/GIBUUInputHandler.cxx
+++ b/src/InputHandler/GIBUUInputHandler.cxx
@@ -1,301 +1,302 @@
#ifdef __GiBUU_ENABLED__
#include "GIBUUInputHandler.h"
#include "InputUtils.h"
GIBUUGeneratorInfo::~GIBUUGeneratorInfo() { DeallocateParticleStack(); }
void GIBUUGeneratorInfo::AddBranchesToTree(TTree* tn) {
// tn->Branch("NEUTParticleN", fNEUTParticleN, "NEUTParticleN/I");
// tn->Branch("NEUTParticleStatusCode", fNEUTParticleStatusCode,
// "NEUTParticleStatusCode[NEUTParticleN]/I");
// tn->Branch("NEUTParticleAliveCode", fNEUTParticleAliveCode,
// "NEUTParticleAliveCode[NEUTParticleN]/I");
}
void GIBUUGeneratorInfo::SetBranchesFromTree(TTree* tn) {
// tn->SetBranchAddress("NEUTParticleN", &fNEUTParticleN );
// tn->SetBranchAddress("NEUTParticleStatusCode", &fNEUTParticleStatusCode );
// tn->SetBranchAddress("NEUTParticleAliveCode", &fNEUTParticleAliveCode );
}
void GIBUUGeneratorInfo::AllocateParticleStack(int stacksize) {
// fNEUTParticleN = 0;
// fNEUTParticleStatusCode = new int[stacksize];
// fNEUTParticleStatusCode = new int[stacksize];
}
void GIBUUGeneratorInfo::DeallocateParticleStack() {
// delete fNEUTParticleStatusCode;
// delete fNEUTParticleAliveCode;
}
void GIBUUGeneratorInfo::FillGeneratorInfo(GiBUUStdHepReader* nevent) {
Reset();
// for (int i = 0; i < nevent->Npart(); i++) {
// fNEUTParticleStatusCode[i] = nevent->PartInfo(i)->fStatus;
// fNEUTParticleAliveCode[i] = nevent->PartInfo(i)->fIsAlive;
// fNEUTParticleN++;
// }
}
void GIBUUGeneratorInfo::Reset() {
// for (int i = 0; i < fNEUTParticleN; i++) {
// fNEUTParticleStatusCode[i] = -1;
// fNEUTParticleAliveCode[i] = 9;
// }
// fNEUTParticleN = 0;
}
GIBUUInputHandler::GIBUUInputHandler(std::string const& handle,
std::string const& rawinputs) {
LOG(SAM) << "Creating GiBUUInputHandler : " << handle << std::endl;
// Run a joint input handling
fName = handle;
fEventType = kGiBUU;
fGIBUUTree = new TChain("giRooTracker");
// Loop over all inputs and grab flux, eventhist, and nevents
std::vector inputs = InputUtils::ParseInputFileList(rawinputs);
for (size_t inp_it = 0; inp_it < inputs.size(); ++inp_it) {
// Open File for histogram access
LOG(SAM) << "Opening event file " << inputs[inp_it] << std::endl;
TFile* inp_file = new TFile(inputs[inp_it].c_str(), "READ");
if ((!inp_file) || (!inp_file->IsOpen())) {
THROW("GiBUU file !IsOpen() at : '"
<< inputs[inp_it] << "'" << std::endl
<< "Check that your file paths are correct and the file exists!");
}
int NFluxes = bool(dynamic_cast(inp_file->Get("numu_flux"))) +
bool(dynamic_cast(inp_file->Get("numub_flux"))) +
bool(dynamic_cast(inp_file->Get("nue_flux"))) +
bool(dynamic_cast(inp_file->Get("nueb_flux"))) +
bool(dynamic_cast(inp_file->Get("e_flux")));
if (NFluxes != 1) {
THROW("Found " << NFluxes << " input fluxes in " << inputs[inp_it]
<< ". The NUISANCE GiBUU interface expects to be "
"passed multiple species vectors as separate "
"input files like: "
"\"GiBUU:(MINERVA_FHC_numu_evts.root,MINERVA_FHC_"
"numubar_evts.root,[...])\"");
}
// Get Flux/Event hist
TH1D* fluxhist = dynamic_cast(inp_file->Get("flux"));
TH1D* eventhist = dynamic_cast(inp_file->Get("evt"));
if (!fluxhist || !eventhist) {
ERROR(FTL, "Input File Contents: " << inputs[inp_it]);
inp_file->ls();
THROW(
"GiBUU FILE doesn't contain flux/xsec info. You may have to "
"regenerate your MC!");
}
// Get N Events
TTree* giRooTracker = dynamic_cast(inp_file->Get("giRooTracker"));
if (!giRooTracker) {
ERROR(FTL,
"giRooTracker Tree not located in NEUT file: " << inputs[inp_it]);
THROW("Check your inputs, they may need to be completely regenerated!");
throw;
}
int nevents = giRooTracker->GetEntries();
if (nevents <= 0) {
THROW("Trying to a TTree with "
<< nevents << " to TChain from : " << inputs[inp_it]);
}
// Register input to form flux/event rate hists
RegisterJointInput(inputs[inp_it], nevents, fluxhist, eventhist);
// Add To TChain
fGIBUUTree->AddFile(inputs[inp_it].c_str());
}
// Registor all our file inputs
SetupJointInputs();
// Create Fit Event
fNUISANCEEvent = new FitEvent();
fGiReader = new GiBUUStdHepReader();
fGiReader->SetBranchAddresses(fGIBUUTree);
fNUISANCEEvent->HardReset();
};
FitEvent* GIBUUInputHandler::GetNuisanceEvent(const UInt_t entry,
const bool lightweight) {
// Check out of bounds
if (entry >= (UInt_t)fNEvents) return NULL;
// Read Entry from TTree to fill NEUT Vect in BaseFitEvt;
fGIBUUTree->GetEntry(entry);
// Run NUISANCE Vector Filler
if (!lightweight) {
CalcNUISANCEKinematics();
}
#ifdef __PROB3PP_ENABLED__
else {
for (int i = 0; i < fGiReader->StdHepN; i++) {
int state = GetGIBUUParticleStatus(fGiReader->StdHepStatus[i],
fGiReader->StdHepPdg[i]);
if (state != kInitialState) {
continue;
}
if (std::count(PhysConst::pdg_neutrinos, PhysConst::pdg_neutrinos + 4,
fGiReader->StdHepPdg[i])) {
fNUISANCEEvent->probe_E = fGiReader->StdHepP4[i][3] * 1.E3;
fNUISANCEEvent->probe_pdg = fGiReader->StdHepPdg[i];
break;
}
}
}
#endif
fNUISANCEEvent->InputWeight *= GetInputWeight(entry);
+ fNUISANCEEvent->GiRead = fGiReader;
return fNUISANCEEvent;
}
int GetGIBUUParticleStatus(int status, int pdg) {
int state = kUndefinedState;
switch (status) {
case 0: // Incoming
case 11: // Struck nucleon
state = kInitialState;
break;
case 1: // Good Final State
state = kFinalState;
break;
default: // Other
break;
}
// Set Nuclear States Flag
if (pdg > 1000000) {
if (state == kInitialState)
state = kNuclearInitial;
else if (state == kFinalState)
state = kNuclearRemnant;
else
state = kUndefinedState;
}
return state;
}
void GIBUUInputHandler::CalcNUISANCEKinematics() {
// Reset all variables
fNUISANCEEvent->ResetEvent();
FitEvent* evt = fNUISANCEEvent;
evt->Mode = fGiReader->GiBUU2NeutCode;
evt->fEventNo = 0.0;
evt->fTotCrs = 0;
evt->fTargetA = 0.0; // Change to get these from nuclear remnant.
evt->fTargetZ = 0.0;
evt->fTargetH = 0;
evt->fBound = 0.0;
// Extra GiBUU Input Weight
evt->InputWeight = fGiReader->EvtWght;
// Check Stack N
int npart = fGiReader->StdHepN;
int kmax = evt->kMaxParticles;
if ((UInt_t)npart > (UInt_t)kmax) {
ERROR(WRN, "GiBUU has too many particles. Expanding Stack.");
fNUISANCEEvent->ExpandParticleStack(npart);
}
// Create Stack
evt->fNParticles = 0;
for (int i = 0; i < npart; i++) {
// State
int state = GetGIBUUParticleStatus(fGiReader->StdHepStatus[i],
fGiReader->StdHepPdg[i]);
int curpart = evt->fNParticles;
// Set State
evt->fParticleState[evt->fNParticles] = state;
// Mom
evt->fParticleMom[curpart][0] = fGiReader->StdHepP4[i][0] * 1.E3;
evt->fParticleMom[curpart][1] = fGiReader->StdHepP4[i][1] * 1.E3;
evt->fParticleMom[curpart][2] = fGiReader->StdHepP4[i][2] * 1.E3;
evt->fParticleMom[curpart][3] = fGiReader->StdHepP4[i][3] * 1.E3;
// PDG
evt->fParticlePDG[curpart] = fGiReader->StdHepPdg[i];
// Add to total particles
evt->fNParticles++;
}
// Run Initial, FSI, Final, Other ordering.
fNUISANCEEvent->OrderStack();
FitParticle* ISNeutralLepton =
fNUISANCEEvent->GetHMISParticle(PhysConst::pdg_neutrinos);
if (ISNeutralLepton) {
fNUISANCEEvent->probe_E = ISNeutralLepton->E();
fNUISANCEEvent->probe_pdg = ISNeutralLepton->PDG();
}
return;
}
void GIBUUInputHandler::Print() {}
void GIBUUInputHandler::SetupJointInputs() {
if (jointeventinputs.size() <= 1) {
jointinput = false;
} else if (jointeventinputs.size() > 1) {
jointinput = true;
jointindexswitch = 0;
}
fMaxEvents = FitPar::Config().GetParI("MAXEVENTS");
if (fMaxEvents != -1 and jointeventinputs.size() > 1) {
THROW("Can only handle joint inputs when config MAXEVENTS = -1!");
}
for (size_t i = 0; i < jointeventinputs.size(); i++) {
double scale = double(fNEvents) / fEventHist->Integral("width");
scale *= jointfluxinputs.at(i)->Integral("width");
jointindexscale.push_back(scale);
}
fEventHist->SetNameTitle((fName + "_EVT").c_str(), (fName + "_EVT").c_str());
fFluxHist->SetNameTitle((fName + "_FLUX").c_str(), (fName + "_FLUX").c_str());
// Setup Max Events
if (fMaxEvents > 1 && fMaxEvents < fNEvents) {
if (LOG_LEVEL(SAM)) {
std::cout << "\t\t|-> Read Max Entries : " << fMaxEvents << std::endl;
}
fNEvents = fMaxEvents;
}
// Print out Status
if (LOG_LEVEL(SAM)) {
std::cout << "\t\t|-> Total Entries : " << fNEvents << std::endl
<< "\t\t|-> Event Integral : "
<< fEventHist->Integral("width") * 1.E-38 << " events/nucleon"
<< std::endl
<< "\t\t|-> Flux Integral : " << fFluxHist->Integral("width")
<< " /cm2" << std::endl
<< "\t\t|-> Event/Flux : "
<< fEventHist->Integral("width") * 1.E-38 /
fFluxHist->Integral("width")
<< " cm2/nucleon" << std::endl;
}
}
#endif
diff --git a/src/InputHandler/NuWroInputHandler.cxx b/src/InputHandler/NuWroInputHandler.cxx
index f4c62ee..384a326 100644
--- a/src/InputHandler/NuWroInputHandler.cxx
+++ b/src/InputHandler/NuWroInputHandler.cxx
@@ -1,473 +1,479 @@
#ifdef __NUWRO_ENABLED__
#include "NuWroInputHandler.h"
#include "InputUtils.h"
NuWroGeneratorInfo::~NuWroGeneratorInfo() { delete fNuWroParticlePDGs; }
void NuWroGeneratorInfo::AddBranchesToTree(TTree* tn) {
tn->Branch("NuWroParticlePDGs", &fNuWroParticlePDGs, "NuWroParticlePDGs/I");
}
void NuWroGeneratorInfo::SetBranchesFromTree(TTree* tn) {
tn->SetBranchAddress("NuWroParticlePDGs", &fNuWroParticlePDGs);
}
void NuWroGeneratorInfo::AllocateParticleStack(int stacksize) {
fNuWroParticlePDGs = new int[stacksize];
}
void NuWroGeneratorInfo::DeallocateParticleStack() {
delete fNuWroParticlePDGs;
}
void NuWroGeneratorInfo::FillGeneratorInfo(event* e) { Reset(); }
void NuWroGeneratorInfo::Reset() {
for (int i = 0; i < kMaxParticles; i++) {
fNuWroParticlePDGs[i] = 0;
}
}
int event1_nof(event* e, int pdg) {
int c = 0;
for (size_t i = 0; i < e->out.size(); i++)
if (e->out[i].pdg == pdg) c++;
return c;
}
NuWroInputHandler::NuWroInputHandler(std::string const& handle,
std::string const& rawinputs) {
LOG(SAM) << "Creating NuWroInputHandler : " << handle << std::endl;
// Run a joint input handling
fName = handle;
fMaxEvents = FitPar::Config().GetParI("MAXEVENTS");
fSaveExtra = false; // FitPar::Config().GetParB("NuWroSaveExtra");
// Setup the TChain
fNuWroTree = new TChain("treeout");
// Loop over all inputs and grab flux, eventhist, and nevents
std::vector inputs = InputUtils::ParseInputFileList(rawinputs);
for (size_t inp_it = 0; inp_it < inputs.size(); ++inp_it) {
// Open File for histogram access
TFile* inp_file = new TFile(inputs[inp_it].c_str(), "READ");
if (!inp_file or inp_file->IsZombie()) {
ERR(FTL) << "nuwro File IsZombie() at " << inputs[inp_it] << std::endl;
throw;
}
// Get Flux/Event hist
TH1D* fluxhist = (TH1D*)inp_file->Get(
(PlotUtils::GetObjectWithName(inp_file, "FluxHist")).c_str());
TH1D* eventhist = (TH1D*)inp_file->Get(
(PlotUtils::GetObjectWithName(inp_file, "EvtHist")).c_str());
if (!fluxhist or !eventhist) {
ERR(FTL) << "nuwro FILE doesn't contain flux/xsec info" << std::endl;
if (FitPar::Config().GetParB("regennuwro")) {
ERR(FTL) << "Regen NuWro has not been added yet. Email the developers!"
<< std::endl;
// ProcessNuWroInputFlux(inputs[inp_it]);
throw;
} else {
ERR(FTL) << "If you would like NUISANCE to generate these for you "
<< "please set parameter regennuwro=1 and re-run."
<< std::endl;
throw;
}
}
// Get N Events
TTree* nuwrotree = (TTree*)inp_file->Get("treeout");
if (!nuwrotree) {
ERR(FTL) << "treeout not located in nuwro file! " << inputs[inp_it]
<< std::endl;
throw;
}
int nevents = nuwrotree->GetEntries();
// Register input to form flux/event rate hists
RegisterJointInput(inputs[inp_it], nevents, fluxhist, eventhist);
// Add to TChain
fNuWroTree->Add(inputs[inp_it].c_str());
}
// Registor all our file inputs
SetupJointInputs();
// Setup Events
fNuWroEvent = NULL;
fNuWroTree->SetBranchAddress("e", &fNuWroEvent);
fNuWroTree->GetEntry(0);
fNUISANCEEvent = new FitEvent();
fNUISANCEEvent->fType = kNUWRO;
fNUISANCEEvent->fNuwroEvent = fNuWroEvent;
fNUISANCEEvent->HardReset();
if (fSaveExtra) {
fNuWroInfo = new NuWroGeneratorInfo();
fNUISANCEEvent->AddGeneratorInfo(fNuWroInfo);
}
};
NuWroInputHandler::~NuWroInputHandler() {
if (fNuWroTree) delete fNuWroTree;
}
void NuWroInputHandler::CreateCache() {
// fNuWroTree->SetCacheEntryRange(0, fNEvents);
// fNuWroTree->AddBranchToCache("*", 1);
// fNuWroTree->SetCacheSize(fCacheSize);
}
void NuWroInputHandler::RemoveCache() {
// fNuWroTree->SetCacheEntryRange(0, fNEvents);
// fNuWroTree->AddBranchToCache("*", 0);
// fNuWroTree->SetCacheSize(0);
}
void NuWroInputHandler::ProcessNuWroInputFlux(const std::string file) {}
FitEvent* NuWroInputHandler::GetNuisanceEvent(const UInt_t entry,
const bool lightweight) {
// Catch too large entries
if (entry >= (UInt_t)fNEvents) return NULL;
// Read Entry from TTree to fill NEUT Vect in BaseFitEvt;
fNuWroTree->GetEntry(entry);
// Run NUISANCE Vector Filler
if (!lightweight) {
CalcNUISANCEKinematics();
}
#ifdef __PROB3PP_ENABLED__
for (size_t i = 0; i < fNUISANCEEvent->fNuwroEvent->in.size(); i++) {
if (std::count(PhysConst::pdg_neutrinos, PhysConst::pdg_neutrinos + 4,
fNUISANCEEvent->fNuwroEvent->in[i].pdg)) {
fNUISANCEEvent->probe_E = fNUISANCEEvent->fNuwroEvent->in[i].t;
fNUISANCEEvent->probe_pdg = fNUISANCEEvent->fNuwroEvent->in[i].pdg;
break;
}
}
#endif
// Setup Input scaling for joint inputs
fNUISANCEEvent->InputWeight = GetInputWeight(entry);
#ifdef __USE_NUWRO_SRW_EVENTS__
if (!rwEvs.size()) {
fNuwroParams = fNuWroEvent->par;
}
-
if (entry >= rwEvs.size()) {
rwEvs.push_back(BaseFitEvt());
+ rwEvs.back().fType = kNUWRO;
+ rwEvs.back().Mode = fNUISANCEEvent->Mode;
rwEvs.back().fNuwroSRWEvent = SRW::SRWEvent(*fNuWroEvent);
rwEvs.back().fNuwroEvent = NULL;
rwEvs.back().fNuwroParams = &fNuwroParams;
rwEvs.back().probe_E = rwEvs.back().fNuwroSRWEvent.NeutrinoEnergy;
rwEvs.back().probe_pdg = rwEvs.back().fNuwroSRWEvent.NeutrinoPDG;
}
+
+ fNUISANCEEvent->fNuwroSRWEvent = SRW::SRWEvent(*fNuWroEvent);
+ fNUISANCEEvent->fNuwroParams = &fNuwroParams;
+ fNUISANCEEvent->probe_E = fNUISANCEEvent->fNuwroSRWEvent.NeutrinoEnergy;
+ fNUISANCEEvent->probe_pdg = fNUISANCEEvent->fNuwroSRWEvent.NeutrinoPDG;
#endif
return fNUISANCEEvent;
}
int NuWroInputHandler::ConvertNuwroMode(event* e) {
Int_t proton_pdg, neutron_pdg, pion_pdg, pion_plus_pdg, pion_minus_pdg,
lambda_pdg, eta_pdg, kaon_pdg, kaon_plus_pdg;
proton_pdg = 2212;
eta_pdg = 221;
neutron_pdg = 2112;
pion_pdg = 111;
pion_plus_pdg = 211;
pion_minus_pdg = -211;
// O_16_pdg = 100069; // oznacznie z Neuta
lambda_pdg = 3122;
kaon_pdg = 311;
kaon_plus_pdg = 321;
if (e->flag.qel) // kwiazielastyczne oddziaływanie
{
if (e->flag.anty) // jeśli jest to oddziaływanie z antyneutrinem
{
if (e->flag.cc)
return -1;
else {
if (event1_nof(e, proton_pdg))
return -51;
else if (event1_nof(e, neutron_pdg))
return -52; // sprawdzam dodatkowo ?
}
} else // oddziaływanie z neutrinem
{
if (e->flag.cc)
return 1;
else {
if (event1_nof(e, proton_pdg))
return 51;
else if (event1_nof(e, neutron_pdg))
return 52;
}
}
}
if (e->flag.mec) {
if (e->flag.anty)
return -2;
else
return 2;
}
if (e->flag.res) // rezonansowa produkcja: pojedynczy pion, pojed.eta, kaon,
// multipiony
{
Int_t liczba_pionow, liczba_kaonow;
liczba_pionow = event1_nof(e, pion_pdg) + event1_nof(e, pion_plus_pdg) +
event1_nof(e, pion_minus_pdg);
liczba_kaonow = event1_nof(e, kaon_pdg) + event1_nof(e, kaon_pdg);
if (liczba_pionow > 1 || liczba_pionow == 0) // multipiony
{
if (e->flag.anty) {
if (e->flag.cc)
return -21;
else
return -41;
} else {
if (e->flag.cc)
return 21;
else
return 41;
}
}
if (liczba_pionow == 1) {
if (e->flag.anty) // jeśli jest to oddziaływanie z antyneutrinem
{
if (e->flag.cc) {
if (event1_nof(e, neutron_pdg) && event1_nof(e, pion_minus_pdg))
return -11;
if (event1_nof(e, neutron_pdg) && event1_nof(e, pion_pdg)) return -12;
if (event1_nof(e, proton_pdg) && event1_nof(e, pion_minus_pdg))
return -13;
} else {
if (event1_nof(e, proton_pdg)) {
if (event1_nof(e, pion_minus_pdg))
return -33;
else if (event1_nof(e, pion_pdg))
return -32;
} else if (event1_nof(e, neutron_pdg)) {
if (event1_nof(e, pion_plus_pdg))
return -34;
else if (event1_nof(e, pion_pdg))
return -31;
}
}
} else // oddziaływanie z neutrinem
{
if (e->flag.cc) {
if (event1_nof(e, proton_pdg) && event1_nof(e, pion_plus_pdg))
return 11;
if (event1_nof(e, proton_pdg) && event1_nof(e, pion_pdg)) return 12;
if (event1_nof(e, neutron_pdg) && event1_nof(e, pion_plus_pdg))
return 13;
} else {
if (event1_nof(e, proton_pdg)) {
if (event1_nof(e, pion_minus_pdg))
return 33;
else if (event1_nof(e, pion_pdg))
return 32;
} else if (event1_nof(e, neutron_pdg)) {
if (event1_nof(e, pion_plus_pdg))
return 34;
else if (event1_nof(e, pion_pdg))
return 31;
}
}
}
}
if (event1_nof(e, eta_pdg)) // produkcja rezonansowa ety
{
if (e->flag.anty) // jeśli jest to oddziaływanie z antyneutrinem
{
if (e->flag.cc)
return -22;
else {
if (event1_nof(e, neutron_pdg))
return -42;
else if (event1_nof(e, proton_pdg))
return -43; // sprawdzam dodatkowo ?
}
} else // oddziaływanie z neutrinem
{
if (e->flag.cc)
return 22;
else {
if (event1_nof(e, neutron_pdg))
return 42;
else if (event1_nof(e, proton_pdg))
return 43;
}
}
}
if (event1_nof(e, lambda_pdg) == 1 &&
liczba_kaonow == 1) // produkcja rezonansowa kaonu
{
if (e->flag.anty) // jeśli jest to oddziaływanie z antyneutrinem
{
if (e->flag.cc && event1_nof(e, kaon_pdg))
return -23;
else {
if (event1_nof(e, kaon_pdg))
return -44;
else if (event1_nof(e, kaon_plus_pdg))
return -45;
}
} else // oddziaływanie z neutrinem
{
if (e->flag.cc && event1_nof(e, kaon_plus_pdg))
return 23;
else {
if (event1_nof(e, kaon_pdg))
return 44;
else if (event1_nof(e, kaon_plus_pdg))
return 45;
}
}
}
}
if (e->flag.coh) // koherentne oddziaływanie tylko na O(16)
{
Int_t _target;
_target = e->par.nucleus_p + e->par.nucleus_n; // liczba masowa O(16)
if (_target == 16) {
if (e->flag.anty) // jeśli jest to oddziaływanie z antyneutrinem
{
if (e->flag.cc && event1_nof(e, pion_minus_pdg))
return -16;
else if (event1_nof(e, pion_pdg))
return -36;
} else // oddziaływanie z neutrinem
{
if (e->flag.cc && event1_nof(e, pion_plus_pdg))
return 16;
else if (event1_nof(e, pion_pdg))
return 36;
}
}
}
// gleboko nieelastyczne rozpraszanie
if (e->flag.dis) {
if (e->flag.anty) {
if (e->flag.cc)
return -26;
else
return -46;
} else {
if (e->flag.cc)
return 26;
else
return 46;
}
}
return 9999;
}
void NuWroInputHandler::CalcNUISANCEKinematics() {
// std::cout << "NuWro Event Address " << fNuWroEvent << std::endl;
// Reset all variables
fNUISANCEEvent->ResetEvent();
FitEvent* evt = fNUISANCEEvent;
// Sort Event Info
evt->Mode = ConvertNuwroMode(fNuWroEvent);
if (abs(evt->Mode) > 60) {
evt->Mode = 0;
}
evt->fEventNo = 0.0;
evt->fTotCrs = 0.0;
evt->fTargetA = fNuWroEvent->par.nucleus_p + fNuWroEvent->par.nucleus_n;
evt->fTargetZ = fNuWroEvent->par.nucleus_p;
evt->fTargetH = 0;
evt->fBound = (evt->fTargetA) == 1;
// Check Particle Stack
UInt_t npart_in = fNuWroEvent->in.size();
UInt_t npart_out = fNuWroEvent->out.size();
UInt_t npart_post = fNuWroEvent->post.size();
UInt_t npart = npart_in + npart_out + npart_post;
UInt_t kmax = evt->kMaxParticles;
if (npart > kmax) {
ERR(WRN) << "NUWRO has too many particles. Expanding stack." << std::endl;
fNUISANCEEvent->ExpandParticleStack(npart);
}
// Sort Particles
evt->fNParticles = 0;
std::vector::iterator p_iter;
// Initial State
for (p_iter = fNuWroEvent->in.begin(); p_iter != fNuWroEvent->in.end();
p_iter++) {
AddNuWroParticle(fNUISANCEEvent, (*p_iter), kInitialState);
}
// FSI State
// for (size_t i = 0; i < npart_in; i++ ) {
// AddNuWroParticle(fNUISANCEEvent, (*p_iter), kFSIState);
// }
// Final State
for (p_iter = fNuWroEvent->post.begin(); p_iter != fNuWroEvent->post.end();
p_iter++) {
AddNuWroParticle(fNUISANCEEvent, (*p_iter), kFinalState);
}
// Fill Generator Info
if (fSaveExtra) fNuWroInfo->FillGeneratorInfo(fNuWroEvent);
// Run Initial, FSI, Final, Other ordering.
fNUISANCEEvent->OrderStack();
FitParticle* ISNeutralLepton =
fNUISANCEEvent->GetHMISParticle(PhysConst::pdg_neutrinos);
if (ISNeutralLepton) {
fNUISANCEEvent->probe_E = ISNeutralLepton->E();
fNUISANCEEvent->probe_pdg = ISNeutralLepton->PDG();
}
return;
}
void NuWroInputHandler::AddNuWroParticle(FitEvent* evt, particle& p,
int state) {
// Add Mom
evt->fParticleMom[evt->fNParticles][0] = static_cast(p).x;
evt->fParticleMom[evt->fNParticles][1] = static_cast(p).y;
evt->fParticleMom[evt->fNParticles][2] = static_cast(p).z;
evt->fParticleMom[evt->fNParticles][3] = static_cast(p).t;
// Status/PDG
evt->fParticleState[evt->fNParticles] = state;
evt->fParticlePDG[evt->fNParticles] = p.pdg;
// Add to particle count
evt->fNParticles++;
}
void NuWroInputHandler::Print() {}
#endif
diff --git a/src/InputHandler/StdHepEvt.cxx b/src/InputHandler/StdHepEvt.cxx
index dfa93f8..02e821a 100644
--- a/src/InputHandler/StdHepEvt.cxx
+++ b/src/InputHandler/StdHepEvt.cxx
@@ -1,135 +1,141 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#include
#include
#include
#include "StdHepEvt.h"
// Include logging
#include "FitLogger.h"
StdHepReader::StdHepReader(){};
bool StdHepReader::SetBranchAddresses(TChain *chain) {
bool ok = true;
int SBAStatus = 0;
SBAStatus = chain->SetBranchAddress("StdHepN", &StdHepN);
ok = ok && (SBAStatus || SBAStatus == 5);
if (!(!SBAStatus || SBAStatus == 5)) {
ERR(WRN) << "Failed to set branch address for \"StdHepN\": " << SBAStatus
<< std::endl;
}
SBAStatus = chain->SetBranchAddress("StdHepPdg", StdHepPdg);
ok = ok && (SBAStatus || SBAStatus == 5);
if (!(!SBAStatus || SBAStatus == 5)) {
ERR(WRN) << "Failed to set branch address for \"StdHepPdg\": " << SBAStatus
<< std::endl;
}
SBAStatus = chain->SetBranchAddress("StdHepStatus", StdHepStatus);
ok = ok && (SBAStatus || SBAStatus == 5);
if (!(!SBAStatus || SBAStatus == 5)) {
ERR(WRN) << "Failed to set branch address for \"StdHepStatus\": "
<< SBAStatus << std::endl;
}
SBAStatus = chain->SetBranchAddress("StdHepP4", StdHepP4);
ok = ok && (SBAStatus || SBAStatus == 5);
if (!(!SBAStatus || SBAStatus == 5)) {
ERR(WRN) << "Failed to set branch address for \"StdHepP4\": " << SBAStatus
<< std::endl;
}
return ok;
}
bool GiBUUStdHepReader::SetBranchAddresses(TChain *chain) {
bool ok = true;
int SBAStatus = 0;
ok = ok && StdHepReader::SetBranchAddresses(chain);
SBAStatus = chain->SetBranchAddress("GiBUU2NeutCode", &GiBUU2NeutCode);
ok = ok && (SBAStatus || SBAStatus == 5);
if (!(!SBAStatus || SBAStatus == 5)) {
ERR(WRN) << "Failed to set branch address for \"GiBUU2NeutCode\": "
<< SBAStatus << std::endl;
}
+ SBAStatus = chain->SetBranchAddress("GiBUUReactionCode", &GiBUUReactionCode);
+ ok = ok && (SBAStatus || SBAStatus == 5);
+ if (!(!SBAStatus || SBAStatus == 5)) {
+ ERR(WRN) << "Failed to set branch address for \"GiBUUReactionCode\": "
+ << SBAStatus << std::endl;
+ }
SBAStatus = chain->SetBranchAddress("EvtWght", &EvtWght);
ok = ok && (SBAStatus || SBAStatus == 5);
if (!(!SBAStatus || SBAStatus == 5)) {
ERR(WRN) << "Failed to set branch address for \"EvtWght\": " << SBAStatus
<< std::endl;
}
return ok;
}
std::string NegSpacer(double const &num) { return (num >= 0) ? " " : ""; }
std::ostream &operator<<(std::ostream &os, TLorentzVector const &tlv) {
std::streamsize prec = os.precision();
std::ios_base::fmtflags flags = os.flags();
os.precision(2);
os.flags(std::ios::scientific);
os << "[" << NegSpacer(tlv[0]) << tlv[0] << "," << NegSpacer(tlv[1]) << tlv[1]
<< "," << NegSpacer(tlv[2]) << tlv[2] << "," << NegSpacer(tlv[3]) << tlv[3]
<< ":M(" << tlv.M() << ")]";
os.precision(prec);
os.flags(flags);
return os;
}
std::string WriteGiBUUEvent(GiBUUStdHepReader const &gi) {
std::stringstream ss("");
ss << "[INFO]: contained " << gi.StdHepN
<< ", Event Weight: " << std::setprecision(3) << gi.EvtWght
<< ", NeutConventionReactionCode: " << gi.GiBUU2NeutCode
<< "\n\t[Lep In](" << std::setw(3)
<< gi.StdHepPdg[0] << ") "
<< TLorentzVector(gi.StdHepP4[0][StdHepReader::kStdHepIdxPx],
gi.StdHepP4[0][StdHepReader::kStdHepIdxPy],
gi.StdHepP4[0][StdHepReader::kStdHepIdxPz],
gi.StdHepP4[0][StdHepReader::kStdHepIdxE])
<< std::endl;
ss << "\t[Target] : " << gi.StdHepPdg[1] << std::endl;
ss << "\t[Nuc In] : "
<< TLorentzVector(gi.StdHepP4[3][StdHepReader::kStdHepIdxPx],
gi.StdHepP4[3][StdHepReader::kStdHepIdxPy],
gi.StdHepP4[3][StdHepReader::kStdHepIdxPz],
gi.StdHepP4[3][StdHepReader::kStdHepIdxE])
<< " (" << std::setw(4) << gi.StdHepPdg[3] << ")" << std::endl;
for (Int_t stdHepInd = 4; stdHepInd < gi.StdHepN; ++stdHepInd) {
ss << "\t[" << std::setw(2) << (stdHepInd - (4)) << "](" << std::setw(5)
<< gi.StdHepPdg[stdHepInd] << ") "
<< TLorentzVector(gi.StdHepP4[stdHepInd][StdHepReader::kStdHepIdxPx],
gi.StdHepP4[stdHepInd][StdHepReader::kStdHepIdxPy],
gi.StdHepP4[stdHepInd][StdHepReader::kStdHepIdxPz],
gi.StdHepP4[stdHepInd][StdHepReader::kStdHepIdxE])
<< std::endl;
}
ss << "\t[Lep Out](" << std::setw(3)
<< gi.StdHepPdg[2] << ") "
<< TLorentzVector(gi.StdHepP4[2][StdHepReader::kStdHepIdxPx],
gi.StdHepP4[2][StdHepReader::kStdHepIdxPy],
gi.StdHepP4[2][StdHepReader::kStdHepIdxPz],
gi.StdHepP4[2][StdHepReader::kStdHepIdxE])
<< std::endl;
return ss.str();
}
diff --git a/src/InputHandler/StdHepEvt.h b/src/InputHandler/StdHepEvt.h
index 198279b..96f280c 100644
--- a/src/InputHandler/StdHepEvt.h
+++ b/src/InputHandler/StdHepEvt.h
@@ -1,111 +1,116 @@
#ifndef __STDHEPEVT_SEEN__
#define __STDHEPEVT_SEEN__
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#include "TChain.h"
#include "TLorentzVector.h"
struct StdHepReader {
public:
const static int kStdHepIdxPx = 0;
const static int kStdHepIdxPy = 1;
const static int kStdHepIdxPz = 2;
const static int kStdHepIdxE = 3;
const static int kStdHepNPmax = 100;
StdHepReader();
///\brief The number of StdHep particles in this event.
Int_t StdHepN;
///\brief The PDG codes of particles in this event.
///
/// This is determined from the GiBUU particle number by
/// GiBUUUtils::GiBUUToPDG.
///\warning This is not a one-to-one mapping, e.g. resonances are not uniquely
/// determined by the GiBUU scheme.
Int_t StdHepPdg[kStdHepNPmax]; //[StdHepN]
///\brief The StdHep Status of particles in this event.
///
/// Status Codes in use:
/// - -1: Initial state real particle.
/// - 1: Final state real particle.
Int_t StdHepStatus[kStdHepNPmax]; //[StdHepN]
///\brief Four momentum for particles in this event.
Double_t StdHepP4[kStdHepNPmax][4];
bool SetBranchAddresses(TChain*);
};
struct GiBUUStdHepReader : public StdHepReader {
GiBUUStdHepReader() : StdHepReader(){};
+ ///\brief GiBUU interaction code
+ ///
+ /// See https://gibuu.hepforge.org/trac/wiki/LesHouches for details
+ Int_t GiBUUReactionCode;
+
///\brief NEUT equivalent reaction code.
/// CC:
/// * 1 : QE
/// * 2 : 2p2h
/// * 10 : Single pion background (non-resonant)
/// * 11 : Delta++ ( -11 : Delta- for nubar)
/// * 12 : Delta+ (-12 : Delta0 for nubar)
/// * 21 : Multi pion production
/// * 26 : DIS
/// * 4 : Higher resonance, charge: -1
/// * 5 : Higher resonance, charge: 0
/// * 6 : Higher resonance, charge: +1
/// * 7 : Higher resonance, charge: +2
///
/// NC:
/// * 30 : Single pion background (non-resonant)
/// * 31 : Delta0
/// * 32 : Delta+
/// * 41 : Multi pion production
/// * 42 : 2p2h
/// * 46 : DIS
/// * 47 : Higher resonance, charge: -1
/// * 48 : Higher resonance, charge: 0
/// * 49 : Higher resonance, charge: +1
/// * 50 : Higher resonance, charge: +2
/// * 51 : NCEL proton-target
/// * 52 : NCEL neutron-target
///
Int_t GiBUU2NeutCode;
///\brief The total XSec weighting that should be applied to this event.
Double_t EvtWght;
///\brief Weighting which takes account of multiple input numu species.
///
/// Defined such that W_numu + W_numubar = 1
Double_t SpeciesWght_numu;
///\brief Weighting which takes account of multiple input nue species.
///
/// Defined such that W_nue + W_nuebar = 1
Double_t SpeciesWght_nue;
///\brief Weighting which takes account of multiple input neutrino species.
///
/// Defined such that \Sum_species W_species = 1
Double_t SpeciesWght;
bool SetBranchAddresses(TChain*);
};
std::string WriteGiBUUEvent(GiBUUStdHepReader const& gi);
#endif
diff --git a/src/Logger/Initialiser.cxx b/src/Logger/Initialiser.cxx
index 923805a..8cb09f2 100644
--- a/src/Logger/Initialiser.cxx
+++ b/src/Logger/Initialiser.cxx
@@ -1,98 +1,101 @@
#include "Initialiser.h"
-void RunNuisance(){
- std::cout << "Starting NUISANCE" << std::endl;
-}
+void RunNuisance() { std::cout << "Starting NUISANCE" << std::endl; }
struct LetterBackronym {
LetterBackronym(size_t n, std::string const &b, float p = 1.0,
std::string const &t = "") {
NUsed = n;
Backkie = b;
ProbAccept = p;
TagLine = t;
};
size_t NUsed;
float ProbAccept;
std::string Backkie;
std::string TagLine;
};
-__attribute__((constructor)) void nuisance_init(void) {
+__attribute__((constructor)) void constructor(void) { nuisance_init(); }
+
+void nuisance_init(void) {
std::vector > Letters;
for (size_t i = 0; i < 8; ++i) {
Letters.push_back(std::vector());
}
Letters[0].push_back(LetterBackronym(2, "Neutrino"));
Letters[0].push_back(LetterBackronym(3, "NUIsance", 0.2));
Letters[2].push_back(LetterBackronym(1, "Interaction"));
Letters[3].push_back(LetterBackronym(1, "Systematics"));
Letters[3].push_back(LetterBackronym(
- 1, "Synthesiser", 0.2, "Playing on the comparisons you want to see"));
+ 1, "Synthesiser", 0.2, "Playing on the comparisons you want to see"));
Letters[4].push_back(LetterBackronym(2, "ANalyser"));
Letters[4].push_back(LetterBackronym(1, "Aggregating", 0.5));
Letters[4].push_back(LetterBackronym(3, "from A-Neutrino sCattering", 1,
"You can always find a frame"));
Letters[5].push_back(
- LetterBackronym(1, "New", 1, "The freshest comparisons"));
+ LetterBackronym(1, "New", 1, "The freshest comparisons"));
Letters[6].push_back(LetterBackronym(1, "by Comparing"));
Letters[6].push_back(LetterBackronym(1, "Constraints from"));
Letters[7].push_back(LetterBackronym(1, "Experiments"));
std::vector TagLines;
TagLines.push_back("Fit and compare.");
std::stringstream back("");
TRandom3 tr;
tr.SetSeed();
for (size_t i = 0; i < 8;) {
LetterBackronym const &let = Letters[i][tr.Integer(Letters[i].size())];
if (tr.Uniform() > let.ProbAccept) {
continue;
}
back << let.Backkie << " ";
i += let.NUsed;
if (let.TagLine.length()) {
TagLines.push_back(let.TagLine);
}
}
std::string Name = "Nuisance";
std::string TagL = TagLines[tr.Integer(TagLines.size())];
std::vector > >
- OneBlob;
+ OneBlob;
OneBlob.push_back(
- std::make_pair("NUISANCE", std::make_pair("", "FiXing your Neutrinos")));
+ std::make_pair("NUISANCE", std::make_pair("", "FiXing your Neutrinos")));
if (tr.Uniform() < 0.01) {
std::pair > const &blob =
- OneBlob[tr.Integer(OneBlob.size())];
+ OneBlob[tr.Integer(OneBlob.size())];
Name = blob.first;
back.str("");
back << blob.second.first;
TagL = blob.second.second;
}
- std::cout << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
- "%%%%%%%%%%%%%%%"
- "%%"
- << std::endl
- << "%% Welcome to " << Name << ": \033[5m" << back.str()
- << "\033[0m-- " << TagL << std::endl
- << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
- "%%%%%%%%%%%%%%%"
- "%%"
- << std::endl;
-}
\ No newline at end of file
+ // std::cout <<
+ // "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
+ // "%%%%%%%%%%%%%%%"
+ // "%%"
+ // << std::endl
+ // << "%% Welcome to " << Name << ": \033[5m" << back.str()
+ // << "\033[0m-- " << TagL << std::endl
+ // <<
+ // "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
+ // "%%%%%%%%%%%%%%%"
+ // "%%"
+ // << std::endl;
+ std::cout << Name << ": " << back.str() << " -- " << TagL << std::endl;
+}
diff --git a/src/Logger/Initialiser.h b/src/Logger/Initialiser.h
index 104ead9..4386235 100644
--- a/src/Logger/Initialiser.h
+++ b/src/Logger/Initialiser.h
@@ -1,20 +1,21 @@
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include