diff --git a/cmake/cacheVariables.cmake b/cmake/cacheVariables.cmake
index 9d51ba7..94d4019 100644
--- a/cmake/cacheVariables.cmake
+++ b/cmake/cacheVariables.cmake
@@ -1,211 +1,213 @@
# Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
################################################################################
# This file is part of NUISANCE.
#
# NUISANCE is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# NUISANCE is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with NUISANCE. If not, see .
################################################################################
function(CheckAndSetDefaultEnv VARNAME DEFAULT CACHETYPE DOCSTRING ENVNAME)
#cmessage(DEBUG "Trying to assign variable ${VARNAME} into the cache.")
if(NOT DEFINED ${VARNAME})
if(DEFINED ENV{${ENVNAME}} AND NOT $ENV{${ENVNAME}} STREQUAL "")
set(${VARNAME} $ENV{${ENVNAME}} CACHE ${CACHETYPE} ${DOCSTRING})
cmessage(DEBUG " Read ${VARNAME} from ENVVAR ${ENVNAME} as $ENV{${ENVNAME}}.")
else()
set(${VARNAME} ${DEFAULT} CACHE ${CACHETYPE} ${DOCSTRING})
endif()
else()
set(${VARNAME} ${${VARNAME}} CACHE ${CACHETYPE} ${DOCSTRING})
unset(${VARNAME})
endif()
cmessage(CACHE "--Set cache variable: \"${VARNAME}\" to \"${${VARNAME}}\", in cache ${CACHETYPE}.")
endfunction()
function(CheckAndSetDefaultCache VARNAME DEFAULT CACHETYPE DOCSTRING)
# cmessage(DEBUG "Trying to assign variable ${VARNAME} into the cache.")
if(NOT DEFINED ${VARNAME})
set(${VARNAME} ${DEFAULT} CACHE ${CACHETYPE} ${DOCSTRING})
else()
set(${VARNAME} ${${VARNAME}} CACHE ${CACHETYPE} ${DOCSTRING})
unset(${VARNAME})
endif()
cmessage(CACHE "--Set cache variable: \"${VARNAME}\" to \"${${VARNAME}}\", in cache ${CACHETYPE}.")
endfunction()
function(CheckAndSetDefault VARNAME DEFAULT)
# cmessage(DEBUG "Trying to assign variable ${VARNAME}.")
if(NOT DEFINED ${VARNAME})
set(${VARNAME} ${DEFAULT} PARENT_SCOPE)
set(${VARNAME} ${DEFAULT})
endif()
cmessage(CACHE "--Set variable: \"${VARNAME}\" to \"${${VARNAME}}\".")
endfunction()
CheckAndSetDefaultCache(VERBOSE TRUE BOOL "Whether to configure loudly.")
set (CMAKE_SKIP_BUILD_RPATH TRUE)
#Changes default install path to be a subdirectory of the build dir.
#Can set build dir at configure time with -DCMAKE_INSTALL_PREFIX=/install/path
if(CMAKE_INSTALL_PREFIX STREQUAL "" OR CMAKE_INSTALL_PREFIX STREQUAL
"/usr/local")
set(CMAKE_INSTALL_PREFIX "${CMAKE_BINARY_DIR}/${CMAKE_SYSTEM_NAME}")
elseif(NOT DEFINED CMAKE_INSTALL_PREFIX)
set(CMAKE_INSTALL_PREFIX "${CMAKE_BINARY_DIR}/${CMAKE_SYSTEM_NAME}")
endif()
if(CMAKE_BUILD_TYPE STREQUAL "")
set(CMAKE_BUILD_TYPE DEBUG)
elseif(NOT DEFINED CMAKE_BUILD_TYPE)
set(CMAKE_BUILD_TYPE DEBUG)
endif()
+CheckAndSetDefaultCache(EXTRA_SETUP_SCRIPT "" PATH "The path to an extra script to inject into the NUISANCE setup script. <>")
+
CheckAndSetDefaultCache(USE_MINIMIZER TRUE INTERNAL "Whether we are using the ROOT minimization libraries. ")
CheckAndSetDefaultCache(USE_HEPMC FALSE BOOL "Whether to enable HepMC input support. ")
CheckAndSetDefaultEnv(HEPMC "" PATH "Path to HepMC source tree root directory. Overrides environment variable \$HEPMC <>" HEPMC)
CheckAndSetDefaultCache(HEPMC_MOMUNIT "GEV" STRING "HepMC momentum units [MEV|GEV]. ")
CheckAndSetDefaultCache(HEPMC_LENUNIT "CM" STRING "HepMC momentum units [MM|CM]. ")
CheckAndSetDefaultCache(HEPMC_USED_EP FALSE INTERNAL "Whether we built HepMC or not. ")
CheckAndSetDefaultCache(USE_NEUT FALSE BOOL "Whether to enable NEUT (reweight) support. Requires external libraries. ")
CheckAndSetDefaultEnv(NEUT_ROOT "" PATH "Path to NEUT source tree root directory. Overrides environment variable \$NEUT_ROOT <>" NEUT_ROOT)
CheckAndSetDefaultEnv(CERN "" PATH "Path to CERNLIB source tree root directory that NEUT was built against. Overrides environment variable \$CERN <>" CERN)
CheckAndSetDefaultEnv(CERN_LEVEL "" STRING "CERNLIB Library version. Overrides environment variable \$CERN_LEVEL <>" CERN_LEVEL)
CheckAndSetDefaultCache(USE_NuWro FALSE BOOL "Whether to enable NuWro support. ")
CheckAndSetDefaultEnv(NUWRO "" PATH "Path to NuWro source tree root directory. Overrides environment variable \$NUWRO <>" NUWRO)
CheckAndSetDefaultEnv(NUWRO_INC "" PATH "Path to NuWro installed includes directory, needs to contain \"params_all.h\". Overrides environment variable \$NUWRO_INC <>" NUWRO_INC)
CheckAndSetDefaultCache(NUWRO_INPUT_FILE "" FILEPATH "Path to an input NuWro event vector, which can be used to build NuWro i/o libraries. <>")
CheckAndSetDefaultCache(NUWRO_BUILT_FROM_FILE FALSE INTERNAL "Whether the NuWro libraries were built by NUISANCE. ")
CheckAndSetDefaultCache(USE_NuWro_RW FALSE BOOL "Whether to try and build support for NuWro reweighting. ")
CheckAndSetDefaultCache(USE_NuWro_SRW_Event FALSE BOOL "Whether to use cut down NuWro reweight event format. Requires NuWro reweight. ")
CheckAndSetDefaultCache(USE_GENIE FALSE BOOL "Whether to enable GENIE (reweight) support. Requires external libraries. ")
CheckAndSetDefaultCache(GENIE_VERSION "AUTO" STRING "GENIE Version ")
CheckAndSetDefaultEnv(GENIE "" PATH "Path to GENIE source tree root directory. Overrides environment variable \$GENIE <>" GENIE)
CheckAndSetDefaultEnv(LHAPDF_LIB "" PATH "Path to pre-built LHAPDF libraries. Overrides environment variable \$LHAPDF_LIB. <>" LHAPDF_LIB)
CheckAndSetDefaultEnv(LHAPDF_INC "" PATH "Path to installed LHAPDF headers. Overrides environment variable \$LHAPDF_INC. <>" LHAPDF_INC)
CheckAndSetDefaultEnv(LHAPATH "" PATH "Path to LHA PDF inputs. Overrides environment variable \$LHAPATH. <>" LHAPATH)
CheckAndSetDefaultEnv(LIBXML2_LIB "" PATH "Path to pre-built LIBXML2 libraries. Overrides environment variable \$LIBXML2_LIB. <>" LIBXML2_LIB)
CheckAndSetDefaultEnv(LIBXML2_INC "" PATH "Path to installed LIBXML2 headers. Overrides environment variable \$LIBXML2_INC. <>" LIBXML2_INC)
CheckAndSetDefaultEnv(LOG4CPP_LIB "" PATH "Path to pre-built LOG4CPP libraries. Overrides environment variable \$LOG4CPP_LIB. <>" LOG4CPP_LIB)
CheckAndSetDefaultEnv(LOG4CPP_INC "" PATH "Path to installed LOG4CPP headers. Overrides environment variable \$LOG4CPP_INC. <>" LOG4CPP_INC)
CheckAndSetDefaultCache(BUILD_GEVGEN FALSE BOOL "Whether to build nuisance_gevgen app.")
CheckAndSetDefaultCache(USE_T2K FALSE BOOL "Whether to enable T2KReWeight support. Requires external libraries. ")
CheckAndSetDefaultEnv(T2KREWEIGHT "" PATH "Path to installed T2KREWEIGHTReWeight. Overrides environment variable \$T2KREWEIGHT. <>" T2KREWEIGHT)
CheckAndSetDefaultCache(USE_NIWG FALSE BOOL "Whether to enable (T2K) NIWG ReWeight support. Requires external libraries. ")
CheckAndSetDefaultEnv(NIWG_ROOT "" PATH "Path to installed NIWGReWeight. Overrides environment variable \$NIWG. <>" NIWG)
CheckAndSetDefaultCache(USE_MINERvA_RW FALSE BOOL "Whether to enable MINERvA ReWeight support. ")
CheckAndSetDefaultEnv(PYTHIA6 "" PATH "Path to directory containing libPythia6.so. Overrides environment variable \$PYTHIA6 <>" PYTHIA6)
CheckAndSetDefaultEnv(PYTHIA8 "" PATH "Path to directory containing libPythia8.so. Overrides environment variable \$PYTHIA8 <>" PYTHIA8)
CheckAndSetDefaultCache(USE_PYTHIA8 FALSE BOOL "Whether to enable PYTHIA8 event support. ")
CheckAndSetDefaultCache(USE_GiBUU TRUE BOOL "Whether to enable GiBUU event support. ")
CheckAndSetDefaultCache(BUILD_GiBUU FALSE BOOL "Whether to build supporting GiBUU event tools along with a patched version of GiBUU. ")
CheckAndSetDefaultCache(USE_NUANCE TRUE BOOL "Whether to enable NUANCE event support. ")
CheckAndSetDefaultCache(USE_PROB3PP FALSE BOOL "Whether to download and compile in Prob3++ support. ")
CheckAndSetDefaultCache(NO_EXTERNAL_UPDATE TRUE BOOL "Whether to perform the update target for external dependencies. ")
CheckAndSetDefaultCache(USE_GPERFTOOLS FALSE BOOL "Whether to compile in google performance tools. ")
CheckAndSetDefault(NEED_PYTHIA6 FALSE)
CheckAndSetDefault(NEED_PYTHIA8 FALSE)
CheckAndSetDefault(NEED_ROOTEVEGEN FALSE)
CheckAndSetDefault(NEED_ROOTPYTHIA6 FALSE)
CheckAndSetDefaultCache(USE_OMP FALSE BOOL "Whether to enable multicore features (there currently are none...). ")
CheckAndSetDefaultCache(USE_DYNSAMPLES FALSE BOOL "Whether to enable the dynamic sample loader. ")
CheckAndSetDefault(NO_EXPERIMENTS FALSE)
cmessage(STATUS "NO_EXPERIMENTS: ${NO_EXPERIMENTS}")
CheckAndSetDefaultCache(NO_ANL ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build ANL samples. <-DNO_EXPERIMENTS=FALSE>")
CheckAndSetDefaultCache(NO_ArgoNeuT ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build ArgoNeuT samples. <-DNO_EXPERIMENTS=FALSE>")
CheckAndSetDefaultCache(NO_BEBC ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build BEBC samples. <-DNO_EXPERIMENTS=FALSE>")
CheckAndSetDefaultCache(NO_BNL ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build BNL samples. <-DNO_EXPERIMENTS=FALSE>")
CheckAndSetDefaultCache(NO_FNAL ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build FNAL samples. <-DNO_EXPERIMENTS=FALSE>")
CheckAndSetDefaultCache(NO_GGM ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build GGM samples. <-DNO_EXPERIMENTS=FALSE>")
CheckAndSetDefaultCache(NO_K2K ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build K2K samples. <-DNO_EXPERIMENTS=FALSE>")
CheckAndSetDefaultCache(NO_MINERvA ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build MINERvA samples. <-DNO_EXPERIMENTS=FALSE>")
CheckAndSetDefaultCache(NO_MiniBooNE ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build MiniBooNE samples. <-DNO_EXPERIMENTS=FALSE>")
CheckAndSetDefaultCache(NO_T2K ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build T2K samples. <-DNO_EXPERIMENTS=FALSE>")
CheckAndSetDefaultCache(NO_SciBooNE ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build SciBooNE samples. <-DNO_EXPERIMENTS=FALSE>")
function(SAYVARS)
LIST(APPEND VARS
USE_HEPMC
HEPMC
HEPMC_MOMUNIT
HEPMC_LENUNIT
HEPMC_USED_EP
USE_NEUT
NEUT_ROOT
CERN
CERN_LEVEL
USE_NuWro
NUWRO
NUWRO_INC
NUWRO_INPUT_FILE
NUWRO_BUILT_FROM_FILE
USE_GENIE
GENIE
LHAPDF_LIB
LHAPDF_INC
LIBXML2_LIB
LIBXML2_INC
LOG4CPP_LIB
GENIE_LOG4CPP_INC
BUILD_GEVGEN
USE_T2K
USE_NIWG
USE_GiBUU
BUILD_GiBUU
USE_NUANCE
NO_EXTERNAL_UPDATE
USE_GPERFTOOLS
NO_ANL
NO_ArgoNeuT
NO_BEBC
NO_BNL
NO_FNAL
NO_GGM
NO_K2K
NO_MINERvA
NO_MiniBooNE
NO_T2K
NO_SciBooNE)
foreach(v ${VARS})
if(DEFINED ${v})
cmessage(DEBUG "VARIABLE: \"${v}\" = \"${${v}}\"")
endif()
endforeach(v)
endfunction()
diff --git a/cmake/setup.sh.in b/cmake/setup.sh.in
index cd1adcf..a090bf6 100644
--- a/cmake/setup.sh.in
+++ b/cmake/setup.sh.in
@@ -1,146 +1,155 @@
# 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 .
################################################################################
#!/bin/sh
### Adapted from https://unix.stackexchange.com/questions/4965/keep-duplicates-out-of-path-on-source
function add_to_PATH () {
for d; do
d=$(cd -- "$d" && { pwd -P || pwd; }) 2>/dev/null # canonicalize symbolic links
if [ -z "$d" ]; then continue; fi # skip nonexistent directory
if [ "$d" == "/usr/bin" ] || [ "$d" == "/usr/bin64" ] || [ "$d" == "/usr/local/bin" ] || [ "$d" == "/usr/local/bin64" ]; then
case ":$PATH:" in
*":$d:"*) :;;
*) export PATH=$PATH:$d;;
esac
else
case ":$PATH:" in
*":$d:"*) :;;
*) export PATH=$d:$PATH;;
esac
fi
done
}
function add_to_LD_LIBRARY_PATH () {
for d; do
d=$(cd -- "$d" && { pwd -P || pwd; }) 2>/dev/null # canonicalize symbolic links
if [ -z "$d" ]; then continue; fi # skip nonexistent directory
if [ "$d" == "/usr/lib" ] || [ "$d" == "/usr/lib64" ] || [ "$d" == "/usr/local/lib" ] || [ "$d" == "/usr/local/lib64" ]; then
case ":$LD_LIBRARY_PATH:" in
*":$d:"*) :;;
*) export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$d;;
esac
else
case ":$LD_LIBRARY_PATH:" in
*":$d:"*) :;;
*) export LD_LIBRARY_PATH=$d:$LD_LIBRARY_PATH;;
esac
fi
done
}
+if [ "@EXTRA_SETUP_SCRIPT@" ]; then
+ if [ ! -e @EXTRA_SETUP_SCRIPT@ ]; then
+ echo "[WARN]: Extra setup script \"@EXTRA_SETUP_SCRIPT@\" requested, but could not be found. Skipping..."
+ else
+ echo "[INFO]: Sourcing extra setup from \"@EXTRA_SETUP_SCRIPT@\"."
+ . @EXTRA_SETUP_SCRIPT@
+ fi
+fi
+
add_to_PATH "@CMAKE_INSTALL_PREFIX@/bin"
add_to_LD_LIBRARY_PATH "@CMAKE_INSTALL_PREFIX@/lib"
if [ ! "${ROOTSYS}" ]; then
echo "[INFO]: Sourcing ROOT from: @CMAKE_ROOTSYS@"
source "@CMAKE_ROOTSYS@/bin/thisroot.sh"
fi
if [ "@USE_NEUT@" != "FALSE" ]; then
echo "[INFO]: Adding NEUT library paths to the environment."
export NEUT_ROOT=@NEUT_ROOT@
export CERN=@CERN@
export CERN_LEVEL=@CERN_LEVEL@
add_to_LD_LIBRARY_PATH "${NEUT_ROOT}/lib/Linux_pc" "${NEUT_ROOT}/src/reweight"
fi
if [ "@USE_NuWro@" != "FALSE" ]; then
if [ "@NUWRO_BUILT_FROM_FILE@" == "FALSE" ]; then
echo "[INFO]: Adding NuWro library paths to the environment."
export NUWRO="@NUWRO@"
add_to_PATH "@NUWRO@/bin"
add_to_LD_LIBRARY_PATH "@NUWRO@/build/@CMAKE_SYSTEM_NAME@/lib"
if [ "@NUWRO_INC@" ]; then
export NUWRO_INC=@NUWRO_INC@
fi
else
echo "[INFO]: NuWro support included from input event file."
fi
fi
if [ "@NEED_PYTHIA6@" != "FALSE" ]; then
echo "[INFO]: Adding PYTHIA6 library paths to the environment."
export PYTHIA6="@PYTHIA6@"
add_to_LD_LIBRARY_PATH "@PYTHIA6@"
fi
if [ "@USE_GENIE@" != "FALSE" ]; then
echo "[INFO]: Adding GENIE paths to the environment."
export GENIE="@GENIE@"
export LHAPDF_LIB="@LHAPDF_LIB@"
export LHAPDF_INC="@LHAPDF_INC@"
export LIBXML2_LIB="@LIBXML2_LIB@"
export LIBXML2_INC="@LIBXML2_INC@"
export LOG4CPP_LIB="@LOG4CPP_LIB@"
export LOG4CPP_INC="@LOG4CPP_INC@"
if [ "@LHAPATH@" ]; then
export LHAPATH="@LHAPATH@"
fi
add_to_PATH "@GENIE@/bin"
add_to_LD_LIBRARY_PATH "@GENIE@/lib" "@LHAPDF_LIB@" "@LIBXML2_LIB@" "@LOG4CPP_LIB@"
fi
if [ "@USE_NIWG@" != "FALSE" ]; then
echo "[INFO]: Adding NIWG paths to the environment."
export NIWG=@NIWG_ROOT@
export NIWGREWEIGHT_INPUTS=@NIWG_ROOT@/inputs
add_to_LD_LIBRARY_PATH "@NIWG_ROOT@"
fi
if [ "@USE_T2K@" != "FALSE" ]; then
echo "[INFO]: Adding T2K paths to the environment."
export T2KREWEIGHT=@T2KREWEIGHT@
add_to_LD_LIBRARY_PATH "@T2KREWEIGHT@/lib"
fi
if [ "@BUILD_GiBUU@" != "FALSE" ]; then
echo "[INFO]: Sourcing GiBUU tools."
source @CMAKE_BINARY_DIR@/GiBUUTools/src/GiBUUTools-build/Linux/setup.sh
fi
export NUISANCE="@CMAKE_SOURCE_DIR@"
diff --git a/src/MCStudies/Smearceptance_Tester.cxx b/src/MCStudies/Smearceptance_Tester.cxx
index 5d11066..db06681 100644
--- a/src/MCStudies/Smearceptance_Tester.cxx
+++ b/src/MCStudies/Smearceptance_Tester.cxx
@@ -1,790 +1,790 @@
// 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 "Smearceptance_Tester.h"
#include "SmearceptanceUtils.h"
#include "Smearcepterton.h"
-#define DEBUG_SMEARTESTER
+//#define DEBUG_SMEARTESTER
//********************************************************************
/// @brief Class to perform smearceptance MC Studies on a custom measurement
Smearceptance_Tester::Smearceptance_Tester(nuiskey samplekey) {
//********************************************************************
samplekey.Print();
// Sample overview ---------------------------------------------------
std::string descrip =
"Simple measurement class for producing an event summary tree of smeared "
"events.\n";
if (Config::HasPar("NPOT")) {
samplekey.SetS("NPOT", Config::GetParS("NPOT"));
}
if (Config::HasPar("FluxIntegralOverride")) {
samplekey.SetS("FluxIntegralOverride",
Config::GetParS("FluxIntegralOverride"));
}
if (Config::HasPar("TargetVolume")) {
samplekey.SetS("TargetVolume", Config::GetParS("TargetVolume"));
}
if (Config::HasPar("TargetMaterialDensity")) {
samplekey.SetS("TargetMaterialDensity",
Config::GetParS("TargetMaterialDensity"));
}
OutputSummaryTree = true;
if (Config::HasPar("smear.OutputSummaryTree")) {
OutputSummaryTree = Config::GetParI("smear.OutputSummaryTree");
}
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetTitle("Smearceptance Studies");
fSettings.SetDescription(descrip);
fSettings.SetXTitle("XXX");
fSettings.SetYTitle("Number of events");
fSettings.SetEnuRange(0.0, 1E5);
fSettings.SetAllowedTypes("EVT/SHAPE/DIAG", "EVT/SHAPE/DIAG");
fSettings.DefineAllowedTargets("*");
fSettings.DefineAllowedSpecies("*");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor =
(GetEventHistogram()->Integral("width") * 1E-38 / (fNEvents + 0.)) /
TotalIntegratedFlux();
// Measurement Details
std::vector splitName = GeneralUtils::ParseToStr(fName, "_");
size_t firstUS = fName.find_first_of("_");
std::string smearceptorName = samplekey.GetS("smearceptor");
QLOG(SAM, "Using smearceptor: " << smearceptorName
<< " (parsed from: " << fName << ").");
fDataHist = new TH1D(("empty_data"), ("empty-data"), 1, 0, 1);
SetupDefaultHist();
fFullCovar = StatUtils::MakeDiagonalCovarMatrix(fDataHist);
covar = StatUtils::GetInvert(fFullCovar);
eventVariables = NULL;
QLOG(SAM, "Smearceptance Flux Scaling Factor = " << fScaleFactor);
if (fScaleFactor <= 0.0) {
ERROR(WRN, "SCALE FACTOR TOO LOW ");
sleep(20);
}
// Setup our TTrees
AddEventVariablesToTree();
smearceptor = &Smearcepterton::Get().GetSmearcepter(smearceptorName);
Int_t RecNBins = 20, TrueNBins = 20;
double RecBinL = 0xdeadbeef, TrueBinL = 0, RecBinH = 10, TrueBinH = 10;
if (Config::HasPar("smear.reconstructed.binning")) {
std::vector args = GeneralUtils::ParseToStr(
Config::GetParS("smear.reconstructed.binning"), ",");
RecNBins = GeneralUtils::StrToInt(args[0]);
RecBinL = GeneralUtils::StrToDbl(args[1]);
RecBinH = GeneralUtils::StrToDbl(args[2]);
TrueNBins = RecNBins;
TrueBinL = RecBinL;
TrueBinH = RecBinH;
}
if (Config::HasPar("smear.true.binning")) {
std::vector args =
GeneralUtils::ParseToStr(Config::GetParS("smear.true.binning"), ",");
TrueNBins = GeneralUtils::StrToInt(args[0]);
TrueBinL = GeneralUtils::StrToDbl(args[1]);
TrueBinH = GeneralUtils::StrToDbl(args[2]);
}
SVDTruncation = 0;
if (Config::HasPar("smear.true.binning")) {
SVDTruncation = Config::GetParI("smear.SVD.truncation");
QLOG(SAM, "Applying SVD truncation of: " << SVDTruncation)
}
ETrueDistrib = NULL;
ETrueDistrib_noweight = NULL;
ERecDistrib = NULL;
RecoSmear = NULL;
if (RecBinL != 0xdeadbeef) {
QLOG(SAM, "Using binning True: " << TrueNBins << ", [" << TrueBinL << " -- "
<< TrueBinH << "], Rec: " << RecNBins
<< ", [" << RecBinL << " -- " << RecBinH
<< "]");
ETrueDistrib = new TH1D("ELep_rate", ";True E_{#nu};Count", TrueNBins,
TrueBinL, TrueBinH);
ETrueDistrib_noweight =
new TH1D("ELep_rate_noweight", ";True E_{#nu};Count", TrueNBins,
TrueBinL, TrueBinH);
ERecDistrib = new TH1D("ELepRec_rate", ";Rec E_{#nu};Count", RecNBins,
RecBinL, RecBinH);
ETrueDistrib->Sumw2();
ERecDistrib->Sumw2();
RecoSmear =
new TH2D("ELepHadVis_Recon", ";True E_{#nu};Recon. E_{#nu}", RecNBins,
RecBinL, RecBinH, TrueNBins, TrueBinL, TrueBinH);
RecoSmear->Sumw2();
}
// Final setup ---------------------------------------------------
FinaliseMeasurement();
}
void Smearceptance_Tester::AddEventVariablesToTree() {
if (OutputSummaryTree) {
// Setup the TTree to save everything
if (!eventVariables) {
Config::Get().out->cd();
eventVariables =
new TTree((fName + "_VARS").c_str(), (fName + "_VARS").c_str());
}
LOG(SAM) << "Adding Event Variables" << std::endl;
eventVariables->Branch("Omega_true", &Omega_true, "Omega_true/F");
eventVariables->Branch("Q2_true", &Q2_true, "Q2_true/F");
eventVariables->Branch("Mode_true", &Mode_true, "Mode_true/I");
eventVariables->Branch("EISLep_true", &EISLep_true, "EISLep_true/F");
eventVariables->Branch("HMFS_mu_true", &HMFS_mu_true);
eventVariables->Branch("HMFS_pip_true", &HMFS_pip_true);
eventVariables->Branch("HMFS_pim_true", &HMFS_pim_true);
eventVariables->Branch("HMFS_cpi_true", &HMFS_cpi_true);
eventVariables->Branch("HMFS_p_true", &HMFS_p_true);
eventVariables->Branch("KEFSHad_cpip_true", &KEFSHad_cpip_true,
"KEFSHad_cpip_true/F");
eventVariables->Branch("KEFSHad_cpim_true", &KEFSHad_cpim_true,
"KEFSHad_cpim_true/F");
eventVariables->Branch("KEFSHad_cpi_true", &KEFSHad_cpi_true,
"KEFSHad_cpi_true/F");
eventVariables->Branch("TEFSHad_pi0_true", &TEFSHad_pi0_true,
"TEFSHad_pi0_true/F");
eventVariables->Branch("KEFSHad_p_true", &KEFSHad_p_true,
"KEFSHad_p_true/F");
eventVariables->Branch("KEFSHad_n_true", &KEFSHad_n_true,
"KEFSHad_n_true/F");
eventVariables->Branch("EFSHad_true", &EFSHad_true, "EFSHad_true/F");
eventVariables->Branch("EFSChargedEMHad_true", &EFSChargedEMHad_true,
"EFSChargedEMHad_true/F");
eventVariables->Branch("EFSLep_true", &EFSLep_true, "EFSLep_true/F");
eventVariables->Branch("EFSgamma_true", &EFSgamma_true, "EFSgamma_true/F");
eventVariables->Branch("PDGISLep_true", &PDGISLep_true, "PDGISLep_true/I");
eventVariables->Branch("PDGFSLep_true", &PDGFSLep_true, "PDGFSLep_true/I");
eventVariables->Branch("Nprotons_true", &Nprotons_true, "Nprotons_true/I");
eventVariables->Branch("Nneutrons_true", &Nneutrons_true,
"Nneutrons_true/I");
eventVariables->Branch("Ncpiplus_true", &Ncpiplus_true, "Ncpiplus_true/I");
eventVariables->Branch("Ncpiminus_true", &Ncpiminus_true,
"Ncpiminus_true/I");
eventVariables->Branch("Ncpi_true", &Ncpi_true, "Ncpi_true/I");
eventVariables->Branch("Npi0_true", &Npi0_true, "Npi0_true/I");
eventVariables->Branch("HMFS_mu_rec", &HMFS_mu_rec);
eventVariables->Branch("HMFS_pip_rec", &HMFS_pip_rec);
eventVariables->Branch("HMFS_pim_rec", &HMFS_pim_rec);
eventVariables->Branch("HMFS_cpi_rec", &HMFS_cpi_rec);
eventVariables->Branch("HMFS_p_rec", &HMFS_p_rec);
eventVariables->Branch("KEFSHad_cpip_rec", &KEFSHad_cpip_rec,
"KEFSHad_cpip_rec/F");
eventVariables->Branch("KEFSHad_cpim_rec", &KEFSHad_cpim_rec,
"KEFSHad_cpim_rec/F");
eventVariables->Branch("KEFSHad_cpi_rec", &KEFSHad_cpi_rec,
"KEFSHad_cpi_rec/F");
eventVariables->Branch("TEFSHad_pi0_rec", &TEFSHad_pi0_rec,
"TEFSHad_pi0_rec/F");
eventVariables->Branch("KEFSHad_p_rec", &KEFSHad_p_rec, "KEFSHad_p_rec/F");
eventVariables->Branch("KEFSHad_n_rec", &KEFSHad_n_rec, "KEFSHad_n_rec/F");
eventVariables->Branch("EFSHad_rec", &EFSHad_rec, "EFSHad_rec/F");
eventVariables->Branch("EFSLep_rec", &EFSLep_rec, "EFSLep_rec/F");
eventVariables->Branch("EFSVis_cpip", &EFSVis_cpip, "EFSVis_cpip/F");
eventVariables->Branch("EFSVis_cpim", &EFSVis_cpim, "EFSVis_cpim/F");
eventVariables->Branch("EFSVis_cpi", &EFSVis_cpi, "EFSVis_cpi/F");
eventVariables->Branch("EFSVis_pi0", &EFSVis_pi0, "EFSVis_pi0/F");
eventVariables->Branch("EFSVis_p", &EFSVis_p, "EFSVis_p/F");
eventVariables->Branch("EFSVis_n", &EFSVis_n, "EFSVis_n/F");
eventVariables->Branch("EFSVis_gamma", &EFSVis_gamma, "EFSVis_gamma/F");
eventVariables->Branch("EFSVis_other", &EFSVis_other, "EFSVis_other/F");
eventVariables->Branch("EFSVis", &EFSVis, "EFSVis/F");
eventVariables->Branch("FSCLep_seen", &FSCLep_seen, "FSCLep_seen/I");
eventVariables->Branch("Nprotons_seen", &Nprotons_seen, "Nprotons_seen/I");
eventVariables->Branch("Nneutrons_seen", &Nneutrons_seen,
"Nneutrons_seen/I");
eventVariables->Branch("Ncpip_seen", &Ncpip_seen, "Ncpip_seen/I");
eventVariables->Branch("Ncpim_seen", &Ncpim_seen, "Ncpim_seen/I");
eventVariables->Branch("Ncpi_seen", &Ncpi_seen, "Ncpi_seen/I");
eventVariables->Branch("Npi0_seen", &Npi0_seen, "Npi0_seen/I");
eventVariables->Branch("Nothers_seen", &Nothers_seen, "Nothers_seen/I");
eventVariables->Branch("EISLep_QE_rec", &EISLep_QE_rec, "EISLep_QE_rec/F");
eventVariables->Branch("EISLep_LepHad_rec", &EISLep_LepHad_rec,
"EISLep_LepHad_rec/F");
eventVariables->Branch("EISLep_LepHadVis_rec", &EISLep_LepHadVis_rec,
"EISLep_LepHadVis_rec/F");
eventVariables->Branch("Nprotons_contributed", &Nprotons_contributed,
"Nprotons_contributed/I");
eventVariables->Branch("Nneutrons_contributed", &Nneutrons_contributed,
"Nneutrons_contributed/I");
eventVariables->Branch("Ncpip_contributed", &Ncpip_contributed,
"Ncpip_contributed/I");
eventVariables->Branch("Ncpim_contributed", &Ncpim_contributed,
"Ncpim_contributed/I");
eventVariables->Branch("Ncpi_contributed", &Ncpi_contributed,
"Ncpi_contributed/I");
eventVariables->Branch("Npi0_contributed", &Npi0_contributed,
"Npi0_contributed/I");
eventVariables->Branch("Ngamma_contributed", &Ngamma_contributed,
"Ngamma_contributed/I");
eventVariables->Branch("Nothers_contibuted", &Nothers_contibuted,
"Nothers_contibuted/I");
eventVariables->Branch("Weight", &Weight, "Weight/F");
eventVariables->Branch("RWWeight", &RWWeight, "RWWeight/F");
eventVariables->Branch("InputWeight", &InputWeight, "InputWeight/F");
eventVariables->Branch("FluxWeight", &FluxWeight, "FluxWeight/F");
eventVariables->Branch("EffWeight", &EffWeight, "EffWeight/F");
xsecScaling = fScaleFactor;
eventVariables->Branch("xsecScaling", &xsecScaling, "xsecScaling/F");
eventVariables->Branch("flagCCINC_true", &flagCCINC_true,
"flagCCINC_true/O");
eventVariables->Branch("flagCC0Pi_true", &flagCC0Pi_true,
"flagCC0Pi_true/O");
eventVariables->Branch("flagCC1Pi_true", &flagCC1Pi_true,
"flagCC1Pi_true/O");
eventVariables->Branch("flagCCINC_rec", &flagCCINC_rec, "flagCCINC_rec/O");
eventVariables->Branch("flagCC0Pi_rec", &flagCC0Pi_rec, "flagCC0Pi_rec/O");
eventVariables->Branch("flagCC1Pi_rec", &flagCC1Pi_rec, "flagCC1Pi_rec/O");
}
PredEvtRateWeight = 1;
if (fEvtRateScaleFactor != 0xdeadbeef) {
if (OutputSummaryTree) {
eventVariables->Branch("PredEvtRateWeight", &PredEvtRateWeight,
"PredEvtRateWeight/F");
}
PredEvtRateWeight = fScaleFactor * fEvtRateScaleFactor;
}
}
template
int CountNPdgsSeen(RecoInfo ri, int const (&pdgs)[N]) {
int sum = 0;
for (size_t pdg_it = 0; pdg_it < N; ++pdg_it) {
sum +=
std::count(ri.RecObjClass.begin(), ri.RecObjClass.end(), pdgs[pdg_it]);
}
return sum;
}
template
int CountNNotPdgsSeen(RecoInfo ri, int const (&pdgs)[N]) {
int sum = 0;
for (size_t pdg_it = 0; pdg_it < N; ++pdg_it) {
sum +=
(std::count(ri.RecObjClass.begin(), ri.RecObjClass.end(), pdgs[pdg_it])
? 0
: 1);
}
return sum;
}
template
int CountNPdgsContributed(RecoInfo ri, int const (&pdgs)[N]) {
int sum = 0;
for (size_t pdg_it = 0; pdg_it < N; ++pdg_it) {
sum += std::count(ri.TrueContribPDGs.begin(), ri.TrueContribPDGs.end(),
pdgs[pdg_it]);
}
return sum;
}
template
int CountNNotPdgsContributed(RecoInfo ri, int const (&pdgs)[N]) {
int sum = 0;
for (size_t pdg_it = 0; pdg_it < N; ++pdg_it) {
sum += (std::count(ri.TrueContribPDGs.begin(), ri.TrueContribPDGs.end(),
pdgs[pdg_it])
? 0
: 1);
}
return sum;
}
TLorentzVector GetHMFSRecParticles(RecoInfo ri, int pdg) {
TLorentzVector mom(0, 0, 0, 0);
for (size_t p_it = 0; p_it < ri.RecObjMom.size(); ++p_it) {
if ((ri.RecObjClass[p_it] == pdg) &&
(mom.Mag() < ri.RecObjMom[p_it].Mag())) {
mom.SetXYZM(ri.RecObjMom[p_it].X(), ri.RecObjMom[p_it].Y(),
ri.RecObjMom[p_it].Z(),
PhysConst::GetMass(ri.RecObjClass[p_it]) * 1.0E3);
}
}
return mom;
}
template
double SumKE_RecoInfo(RecoInfo ri, int const (&pdgs)[N], double mass) {
double sum = 0;
for (size_t p_it = 0; p_it < ri.RecObjMom.size(); ++p_it) {
if (!std::count(pdgs, pdgs + N,
ri.RecObjClass[p_it])) { // If we don't care about this
// particle type.
continue;
}
sum += sqrt(ri.RecObjMom[p_it].Mag2() + mass * mass) - mass;
}
return sum;
}
template
double SumTE_RecoInfo(RecoInfo ri, int const (&pdgs)[N], double mass) {
double sum = 0;
for (size_t p_it = 0; p_it < ri.RecObjMom.size(); ++p_it) {
if (!std::count(pdgs, pdgs + N,
ri.RecObjClass[p_it])) { // If we don't care about this
// particle type.
continue;
}
sum += sqrt(ri.RecObjMom[p_it].Mag2() + mass * mass);
}
return sum;
}
template
double SumVisE_RecoInfo(RecoInfo ri, int const (&pdgs)[N]) {
double sum = 0;
for (size_t p_it = 0; p_it < ri.RecVisibleEnergy.size(); ++p_it) {
if (!std::count(pdgs, pdgs + N,
ri.TrueContribPDGs[p_it])) { // If we don't care about this
// particle type.
continue;
}
sum += ri.RecVisibleEnergy[p_it];
}
return sum;
}
template
double SumVisE_RecoInfo_NotPdgs(RecoInfo ri, int const (&pdgs)[N]) {
double sum = 0;
for (size_t p_it = 0; p_it < ri.RecVisibleEnergy.size(); ++p_it) {
if (std::count(pdgs, pdgs + N,
ri.TrueContribPDGs[p_it])) { // If we know about this
// particle type.
continue;
}
sum += ri.RecVisibleEnergy[p_it];
}
return sum;
}
//********************************************************************
void Smearceptance_Tester::FillEventVariables(FitEvent *event) {
//********************************************************************
static int const cpipPDG[] = {211};
static int const cpimPDG[] = {-211};
static int const pi0PDG[] = {111};
static int const ProtonPDG[] = {2212};
static int const NeutronPDG[] = {2112};
static int const GammaPDG[] = {22};
static int const CLeptonPDGs[] = {11, 13, 15};
static int const ExplicitPDGs[] = {211, -211, 111, 2212, 2112, 22,
11, 13, 15, 12, 14, 16};
RecoInfo *ri = smearceptor->Smearcept(event);
HMFS_mu_true = TLorentzVector(0, 0, 0, 0);
HMFS_mu_rec = TLorentzVector(0, 0, 0, 0);
FitParticle *fsMu = event->GetHMFSMuon();
if (fsMu) {
HMFS_mu_true = fsMu->P4();
HMFS_mu_rec = GetHMFSRecParticles(*ri, 13);
}
HMFS_pip_true = TLorentzVector(0, 0, 0, 0);
HMFS_pip_rec = TLorentzVector(0, 0, 0, 0);
FitParticle *fsPip = event->GetHMFSPiPlus();
if (fsPip) {
HMFS_pip_true = fsPip->P4();
HMFS_pip_rec = GetHMFSRecParticles(*ri, 211);
}
HMFS_pim_true = TLorentzVector(0, 0, 0, 0);
HMFS_pim_rec = TLorentzVector(0, 0, 0, 0);
FitParticle *fsPim = event->GetHMFSPiMinus();
if (fsPim) {
HMFS_pim_true = fsPim->P4();
HMFS_pim_rec = GetHMFSRecParticles(*ri, -211);
}
HMFS_cpi_true = TLorentzVector(0, 0, 0, 0);
HMFS_cpi_rec = TLorentzVector(0, 0, 0, 0);
if (fsPip || fsPim) {
if (!fsPip) {
HMFS_cpi_true = HMFS_pim_true;
HMFS_cpi_rec = HMFS_pim_rec;
} else if (!fsPim) {
HMFS_cpi_true = HMFS_pip_true;
HMFS_cpi_rec = HMFS_pip_rec;
} else {
HMFS_cpi_true =
(fsPip->p2() > fsPim->p2()) ? HMFS_pip_true : HMFS_pim_true;
HMFS_cpi_rec = (fsPip->p2() > fsPim->p2()) ? HMFS_pip_rec : HMFS_pim_rec;
}
}
HMFS_p_true = TLorentzVector(0, 0, 0, 0);
HMFS_p_rec = TLorentzVector(0, 0, 0, 0);
FitParticle *fsP = event->GetHMFSProton();
if (fsP) {
HMFS_p_true = fsP->P4();
HMFS_p_rec = GetHMFSRecParticles(*ri, 2212);
}
TLorentzVector FourMomentumTransfer =
(event->GetHMISAnyLeptons()->P4() - event->GetHMFSAnyLeptons()->P4());
Omega_true = FourMomentumTransfer.E();
Q2_true = -1 * FourMomentumTransfer.Mag2();
Mode_true = event->Mode;
EISLep_true = event->GetHMISAnyLeptons()->E();
KEFSHad_cpip_true = FitUtils::SumTE_PartVect(event->GetAllFSPiPlus());
KEFSHad_cpim_true = FitUtils::SumTE_PartVect(event->GetAllFSPiMinus());
KEFSHad_cpi_true = KEFSHad_cpip_true + KEFSHad_cpim_true;
TEFSHad_pi0_true = FitUtils::SumTE_PartVect(event->GetAllFSPiZero());
KEFSHad_p_true = FitUtils::SumKE_PartVect(event->GetAllFSProton());
KEFSHad_n_true = FitUtils::SumKE_PartVect(event->GetAllFSNeutron());
EFSHad_true =
KEFSHad_cpi_true + TEFSHad_pi0_true + KEFSHad_p_true + KEFSHad_n_true;
EFSChargedEMHad_true = KEFSHad_cpi_true + TEFSHad_pi0_true + KEFSHad_p_true;
EFSLep_true = event->GetHMFSAnyLeptons()->E();
EFSgamma_true = FitUtils::SumTE_PartVect(event->GetAllFSPhoton());
PDGISLep_true = event->GetHMISAnyLeptons()->PDG();
PDGFSLep_true = event->GetHMFSAnyLeptons()->PDG();
Nprotons_true = event->GetAllFSProton().size();
Nneutrons_true = event->GetAllFSNeutron().size();
Ncpiplus_true = event->GetAllFSPiPlus().size();
Ncpiminus_true = event->GetAllFSPiMinus().size();
Ncpi_true = Ncpiplus_true + Ncpiminus_true;
Npi0_true = event->GetAllFSPiZero().size();
KEFSHad_cpip_rec =
SumKE_RecoInfo(*ri, cpipPDG, PhysConst::mass_cpi * PhysConst::mass_MeV);
KEFSHad_cpim_rec =
SumKE_RecoInfo(*ri, cpimPDG, PhysConst::mass_cpi * PhysConst::mass_MeV);
KEFSHad_cpi_rec = KEFSHad_cpip_rec + KEFSHad_cpim_rec;
TEFSHad_pi0_rec =
SumTE_RecoInfo(*ri, pi0PDG, PhysConst::mass_pi0 * PhysConst::mass_MeV);
KEFSHad_p_rec = SumKE_RecoInfo(*ri, ProtonPDG,
PhysConst::mass_proton * PhysConst::mass_MeV);
KEFSHad_n_rec = SumKE_RecoInfo(*ri, NeutronPDG,
PhysConst::mass_neutron * PhysConst::mass_MeV);
EFSHad_rec =
KEFSHad_cpi_rec + TEFSHad_pi0_rec + KEFSHad_p_rec + KEFSHad_n_rec;
TLorentzVector FSLepMom_rec(0, 0, 0, 0);
if (event->GetHMFSAnyLeptons()) {
FSLepMom_rec = GetHMFSRecParticles(*ri, event->GetHMFSAnyLeptons()->PDG());
EFSLep_rec = FSLepMom_rec.E();
} else {
EFSLep_rec = 0;
}
EFSVis_cpip = SumVisE_RecoInfo(*ri, cpipPDG);
EFSVis_cpim = SumVisE_RecoInfo(*ri, cpimPDG);
EFSVis_cpi = EFSVis_cpip + EFSVis_cpim;
EFSVis_pi0 = SumVisE_RecoInfo(*ri, pi0PDG);
EFSVis_p = SumVisE_RecoInfo(*ri, ProtonPDG);
EFSVis_n = SumVisE_RecoInfo(*ri, NeutronPDG);
EFSVis_gamma = SumVisE_RecoInfo(*ri, GammaPDG);
EFSVis_other = SumVisE_RecoInfo_NotPdgs(*ri, ExplicitPDGs);
EFSVis = EFSVis_cpi + EFSVis_pi0 + EFSVis_p + EFSVis_n + EFSVis_gamma;
FSCLep_seen = CountNPdgsSeen(*ri, CLeptonPDGs);
Nprotons_seen = CountNPdgsSeen(*ri, ProtonPDG);
Nneutrons_seen = CountNPdgsSeen(*ri, NeutronPDG);
Ncpip_seen = CountNPdgsSeen(*ri, cpipPDG);
Ncpim_seen = CountNPdgsSeen(*ri, cpimPDG);
Ncpi_seen = Ncpip_seen + Ncpim_seen;
Npi0_seen = CountNPdgsSeen(*ri, pi0PDG);
Nothers_seen = CountNNotPdgsSeen(*ri, ExplicitPDGs);
if (FSCLep_seen && (FSLepMom_rec.Mag() > 1E-8)) {
EISLep_QE_rec =
FitUtils::EnuQErec(FSLepMom_rec.Mag() / 1000.0, FSLepMom_rec.CosTheta(),
34, PDGFSLep_true > 0) *
1000.0;
} else {
EISLep_QE_rec = 0;
}
EISLep_LepHad_rec = EFSLep_rec + EFSHad_rec;
EISLep_LepHadVis_rec = EFSLep_rec + EFSHad_rec + EFSVis;
Nprotons_contributed = CountNPdgsContributed(*ri, ProtonPDG);
Nneutrons_contributed = CountNPdgsContributed(*ri, NeutronPDG);
Ncpip_contributed = CountNPdgsContributed(*ri, cpipPDG);
Ncpim_contributed = CountNPdgsContributed(*ri, cpimPDG);
Ncpi_contributed = Ncpip_contributed + Ncpim_contributed;
Npi0_contributed = CountNPdgsContributed(*ri, pi0PDG);
Ngamma_contributed = CountNPdgsContributed(*ri, GammaPDG);
Nothers_contibuted = CountNNotPdgsContributed(*ri, ExplicitPDGs);
Weight = event->RWWeight * event->InputWeight;
RWWeight = event->RWWeight;
InputWeight = event->InputWeight;
FluxWeight = GetFluxHistogram()->GetBinContent(
GetFluxHistogram()->FindBin(EISLep_true)) /
GetFluxHistogram()->Integral();
EffWeight = ri->Weight;
flagCCINC_true = PDGFSLep_true & 1;
flagCC0Pi_true = (Ncpi_true + Npi0_true) == 0;
flagCC1Pi_true = (Ncpi_true + Npi0_true) == 1;
flagCCINC_rec = FSCLep_seen && PDGFSLep_true & 1;
flagCC0Pi_rec = ((Ncpi_seen + Npi0_seen) == 0) && flagCCINC_rec;
flagCC1Pi_rec = ((Ncpi_seen + Npi0_seen) == 1) && flagCCINC_rec;
if (OutputSummaryTree) {
// Fill the eventVariables Tree
eventVariables->Fill();
}
if (RecoSmear) {
RecoSmear->Fill(EISLep_true / 1000.0,
flagCCINC_rec ? EISLep_LepHadVis_rec / 1000.0 : -1, Weight);
ETrueDistrib_noweight->Fill(EISLep_true / 1000.0,
flagCCINC_true ? Weight : 0);
ETrueDistrib->Fill(EISLep_true / 1000.0,
flagCCINC_true ? Weight * PredEvtRateWeight : 0);
ERecDistrib->Fill(EISLep_LepHadVis_rec / 1000.0,
flagCCINC_rec ? Weight * PredEvtRateWeight : 0);
}
};
//********************************************************************
void Smearceptance_Tester::Write(std::string drawOpt) {
//********************************************************************
if (OutputSummaryTree) {
// First save the TTree
eventVariables->Write();
}
// Save Flux and Event Histograms too
GetInput()->GetFluxHistogram()->Write();
GetInput()->GetEventHistogram()->Write();
if (!RecoSmear) {
return;
}
TH2D *SmearMatrix_ev =
static_cast(RecoSmear->Clone("ELepHadVis_Smear_ev"));
for (Int_t trueAxis_it = 1;
trueAxis_it < RecoSmear->GetXaxis()->GetNbins() + 1; ++trueAxis_it) {
double NEISLep = ETrueDistrib_noweight->GetBinContent(trueAxis_it);
for (Int_t recoAxis_it = 1;
recoAxis_it < RecoSmear->GetYaxis()->GetNbins() + 1; ++recoAxis_it) {
if (NEISLep > std::numeric_limits::epsilon()) {
SmearMatrix_ev->SetBinContent(
trueAxis_it, recoAxis_it,
SmearMatrix_ev->GetBinContent(trueAxis_it, recoAxis_it) / NEISLep);
}
}
}
ETrueDistrib_noweight->Write();
ETrueDistrib->Write();
ERecDistrib->Write();
RecoSmear->Write();
SmearMatrix_ev->Write();
TH2D *ResponseMatrix_ev =
SmearceptanceUtils::SVDGetInverse(SmearMatrix_ev, SVDTruncation);
ResponseMatrix_ev->SetName("ResponseMatrix_ev");
ResponseMatrix_ev->Write();
#ifdef DEBUG_SMEARTESTER
TMatrixD SmearMatrix_ev_md = SmearceptanceUtils::GetMatrix(SmearMatrix_ev);
TH1D *SmearedEvt = static_cast(ERecDistrib->Clone());
SmearedEvt->SetNameTitle("SmearedEvt", ";Rec E_{#nu}; count");
SmearceptanceUtils::PushTH1ThroughMatrixWithErrors(
ETrueDistrib, SmearedEvt, SmearMatrix_ev_md, 5000, false);
SmearedEvt->Write();
SmearedEvt->Scale(1, "width");
SmearedEvt->SetName("SmearedEvt_bw");
SmearedEvt->Write();
#endif
#ifdef __PROB3PP_ENABLED__
FitWeight *fw = FitBase::GetRW();
if (fw->HasRWEngine(kOSCILLATION)) {
OscWeightEngine *oscWE =
dynamic_cast(fw->GetRWEngine(kOSCILLATION));
TGraph POsc;
POsc.Set(1E4 - 1);
double min = ETrueDistrib->GetXaxis()->GetBinLowEdge(1);
double step = (ETrueDistrib->GetXaxis()->GetBinUpEdge(
ETrueDistrib->GetXaxis()->GetNbins()) -
ETrueDistrib->GetXaxis()->GetBinLowEdge(1)) /
double(1E4);
for (size_t i = 1; i < 1E4; ++i) {
double enu = min + i * step;
double ow = oscWE->CalcWeight(enu, 14);
if (ow != ow) {
std::cout << "Bad osc weight for ENu: " << enu << std::endl;
}
POsc.SetPoint(i - 1, enu, ow);
}
POsc.Write("POsc", TObject::kOverwrite);
}
#endif
TMatrixD ResponseMatrix_evt_md =
SmearceptanceUtils::GetMatrix(ResponseMatrix_ev);
TH1D *Unfolded_enu_obs = static_cast(ETrueDistrib->Clone());
Unfolded_enu_obs->SetNameTitle("UnfoldedENu_evt", ";True E_{#nu};count");
SmearceptanceUtils::PushTH1ThroughMatrixWithErrors(
ERecDistrib, Unfolded_enu_obs, ResponseMatrix_evt_md, 5000, false);
Unfolded_enu_obs->Write();
Unfolded_enu_obs->Scale(1, "width");
Unfolded_enu_obs->SetName("UnfoldedENu_evt_bw");
Unfolded_enu_obs->Write();
ETrueDistrib->Scale(1, "width");
ETrueDistrib->SetName("ELep_rate_bw");
ETrueDistrib->Write();
ERecDistrib->Scale(1, "width");
ERecDistrib->SetName("ELepRec_rate_bw");
ERecDistrib->Write();
}
// -------------------------------------------------------------------
// Purely MC Plot
// Following functions are just overrides to handle this
// -------------------------------------------------------------------
//********************************************************************
/// Everything is classed as signal...
bool Smearceptance_Tester::isSignal(FitEvent *event) {
//********************************************************************
(void)event;
return true;
};
//********************************************************************
void Smearceptance_Tester::ScaleEvents() {
//********************************************************************
// Saving everything to a TTree so no scaling required
return;
}
//********************************************************************
void Smearceptance_Tester::ApplyNormScale(float norm) {
//********************************************************************
// Saving everything to a TTree so no scaling required
fCurrentNorm = norm;
return;
}
//********************************************************************
void Smearceptance_Tester::FillHistograms() {
//********************************************************************
// No Histograms need filling........
return;
}
//********************************************************************
void Smearceptance_Tester::ResetAll() {
//********************************************************************
if (OutputSummaryTree) {
eventVariables->Reset();
}
return;
}
//********************************************************************
float Smearceptance_Tester::GetChi2() {
//********************************************************************
// No Likelihood to test, purely MC
return 0.0;
}