diff --git a/CMakeLists.txt b/CMakeLists.txt
index 590b423..0154ce4 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,228 +1,232 @@
# 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 .
################################################################################
-cmake_minimum_required (VERSION 2.6 FATAL_ERROR)
+cmake_minimum_required (VERSION 2.8 FATAL_ERROR)
+
+#Use the compilers found in the path
+find_program(CMAKE_C_COMPILER NAMES $ENV{CC} gcc PATHS ENV PATH NO_DEFAULT_PATH)
+find_program(CMAKE_CXX_COMPILER NAMES $ENV{CXX} g++ PATHS ENV PATH NO_DEFAULT_PATH)
project(NUISANCE)
include(ExternalProject)
enable_language(Fortran)
set (NUISANCE_VERSION_MAJOR 2)
set (NUISANCE_VERSION_MINOR 7)
set (NUISANCE_VERSION_REVISION 0)
set (NUISANCE_VERSION_STRING "v${NUISANCE_VERSION_MAJOR}r${NUISANCE_VERSION_MINOR}")
if(${NUISANCE_VERSION_REVISION} STRGREATER "0")
set (NUISANCE_VERSION_STRING "${NUISANCE_VERSION_STRING}p${NUISANCE_VERSION_REVISION}")
endif()
#Set this to TRUE to enable build debugging messages
set(BUILD_DEBUG_MSGS TRUE)
include(${CMAKE_SOURCE_DIR}/cmake/cmessage.cmake)
include(${CMAKE_SOURCE_DIR}/cmake/cacheVariables.cmake)
cmessage(STATUS "CMAKE_INSTALL_PREFIX: \"${CMAKE_INSTALL_PREFIX}\"")
cmessage(STATUS "CMAKE_BUILD_TYPE: \"${CMAKE_BUILD_TYPE}\"")
################################################################################
# Check Dependencies
################################################################################
################################## ROOT ######################################
include(${CMAKE_SOURCE_DIR}/cmake/ROOTSetup.cmake)
################################# HEPMC ######################################
include(${CMAKE_SOURCE_DIR}/cmake/HepMC.cmake)
############################ Reweight Engines ################################
include(${CMAKE_SOURCE_DIR}/cmake/ReweightEnginesSetup.cmake)
############################ Other Generators ################################
include(${CMAKE_SOURCE_DIR}/cmake/GiBUUSetup.cmake)
if(USE_NUANCE)
LIST(APPEND EXTRA_CXX_FLAGS -D__NUANCE_ENABLED__)
endif()
################################# Pythia6/8 ####################################
include(${CMAKE_SOURCE_DIR}/cmake/pythia6Setup.cmake)
include(${CMAKE_SOURCE_DIR}/cmake/pythia8Setup.cmake)
################################# gperftools ###################################
include(${CMAKE_SOURCE_DIR}/cmake/gperfSetup.cmake)
if(NOT NOTEST)
enable_testing()
endif()
SET(GENERATOR_SUPPORT)
foreach(gen NEUT;NuWro;GENIE;GiBUU;NUANCE)
if(USE_${gen})
SET(GENERATOR_SUPPORT "${GENERATOR_SUPPORT}${gen} ")
endif()
endforeach(gen)
cmessage(STATUS "Generator Input Support: ${GENERATOR_SUPPORT}")
set(MINCODE
Routines
FCN)
set(CORE
MCStudies
Genie
FitBase
Config
Logger
InputHandler
Splines
Reweight
Utils
Statistical
#Devel
Smearceptance
)
LIST(APPEND ALLEXPERIMENTS
ANL
ArgoNeuT
BEBC
BNL
Electron
FNAL
GGM
K2K
MINERvA
MiniBooNE
SciBooNE
T2K)
foreach(exp ${ALLEXPERIMENTS})
if(NOT NO_${exp})
LIST(APPEND EXPERIMENTS_TO_BUILD ${exp})
else()
LIST(REVERSE EXTRA_CXX_FLAGS)
LIST(APPEND EXTRA_CXX_FLAGS -D__NO_${exp}__)
LIST(REVERSE EXTRA_CXX_FLAGS)
endif()
endforeach()
################################## COMPILER ####################################
include(${CMAKE_SOURCE_DIR}/cmake/c++CompilerSetup.cmake)
################################### doxygen ###################################
include(${CMAKE_SOURCE_DIR}/cmake/docsSetup.cmake)
################################################################################
set(MINIMUM_INCLUDE_DIRECTORIES)
LIST(APPEND MINIMUM_INCLUDE_DIRECTORIES
${RWENGINE_INCLUDE_DIRECTORIES}
${CMAKE_SOURCE_DIR}/src/FitBase
${CMAKE_SOURCE_DIR}/src/Reweight
${CMAKE_SOURCE_DIR}/src/InputHandler
${CMAKE_SOURCE_DIR}/src/Config
${CMAKE_SOURCE_DIR}/src/Logger
${CMAKE_SOURCE_DIR}/src/Statistical
${CMAKE_SOURCE_DIR}/src/Splines
${CMAKE_SOURCE_DIR}/src/Utils
${CMAKE_SOURCE_DIR}/src/Genie)
cmessage(DEBUG "Base include directories: ${MINIMUM_INCLUDE_DIRECTORIES}")
set(EXP_INCLUDE_DIRECTORIES)
foreach(edir ${EXPERIMENTS_TO_BUILD})
LIST(APPEND EXP_INCLUDE_DIRECTORIES ${CMAKE_SOURCE_DIR}/src/${edir})
endforeach()
cmessage(DEBUG "Included experiments: ${EXP_INCLUDE_DIRECTORIES}")
foreach(mdir ${MINCODE})
cmessage (DEBUG "Configuring directory: src/${mdir}")
add_subdirectory(src/${mdir})
endforeach()
foreach(edir ${EXPERIMENTS_TO_BUILD})
cmessage (DEBUG "Configuring directory: src/${edir}")
add_subdirectory(src/${edir})
endforeach()
foreach(cdir ${CORE})
cmessage (DEBUG "Configuring directory: src/${cdir}")
add_subdirectory(src/${cdir})
endforeach()
cmessage(DEBUG "Module targets: ${MODULETargets}")
add_subdirectory(app)
add_subdirectory(src/Tests)
configure_file(cmake/setup.sh.in
"${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/setup.sh" @ONLY)
install(FILES
"${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/setup.sh" DESTINATION
${CMAKE_INSTALL_PREFIX})
configure_file(cmake/MakeBinaryBlob.in
"${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/MakeBinaryBlob" @ONLY)
install(PROGRAMS
"${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/MakeBinaryBlob" DESTINATION
bin)
if(USE_DYNSAMPLES)
SET(ALL_INCLUDES ${MINIMUM_INCLUDE_DIRECTORIES})
LIST(APPEND ALL_INCLUDES ${CMAKE_SOURCE_DIR}/src/Smearceptance)
LIST(APPEND ALL_INCLUDES ${EXP_INCLUDE_DIRECTORIES})
string(REPLACE ";" " -I" ALL_INCLUDES_STR "${ALL_INCLUDES}")
cmessage(DEBUG ${CMAKE_DEPENDLIB_FLAGS})
string(REPLACE "-levent " "" CMAKE_DEPENDLIB_FLAGS_NEW ${CMAKE_DEPENDLIB_FLAGS})
set(CMAKE_DEPENDLIB_FLAGS ${CMAKE_DEPENDLIB_FLAGS_NEW})
cmessage(DEBUG ${CMAKE_DEPENDLIB_FLAGS})
string(REPLACE ";" " -l" ALL_MODULETARGETS_STR "${MODULETargets}")
configure_file(cmake/BuildDynamicSample.in
"${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/BuildDynamicSample" @ONLY)
install(PROGRAMS
"${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/BuildDynamicSample" DESTINATION
bin)
configure_file(cmake/BuildDynamicSmearcepter.in
"${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/BuildDynamicSmearcepter" @ONLY)
install(PROGRAMS
"${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/BuildDynamicSmearcepter" DESTINATION
bin)
endif()
install(PROGRAMS
"${PROJECT_SOURCE_DIR}/scripts/nuiscardgen" DESTINATION
bin)
install(PROGRAMS
"${PROJECT_SOURCE_DIR}/scripts/nuissamples" DESTINATION
bin)
diff --git a/cmake/ROOTSetup.cmake b/cmake/ROOTSetup.cmake
index 37f6157..caa6e8e 100644
--- a/cmake/ROOTSetup.cmake
+++ b/cmake/ROOTSetup.cmake
@@ -1,167 +1,169 @@
# 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 "6\..*+)" ROOTVERSIXMATCH ${ROOT_VERSION})
+ string(REGEX MATCH "6.*" ROOTVERSIXMATCH ${ROOT_VERSION})
if(ROOTVERSIXMATCH)
- cmessage(STATUS "Using ROOT6, We are essentiall flying blind here.")
+ cmessage(STATUS "Using ROOT6, We are essentially flying blind here.")
+ LIST(REMOVE_ITEM ROOT_LIBS Cint)
+ LIST(APPEND EXTRA_CXX_FLAGS -DROOT6_USE_FIT_FITTER_INTERFACE)
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/src/Electron/ElectronScattering_DurhamData.cxx b/src/Electron/ElectronScattering_DurhamData.cxx
index c171df5..7e7b616 100644
--- a/src/Electron/ElectronScattering_DurhamData.cxx
+++ b/src/Electron/ElectronScattering_DurhamData.cxx
@@ -1,486 +1,486 @@
// 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 "ElectronScattering_DurhamData.h"
//********************************************************************
ElectronScattering_DurhamData::ElectronScattering_DurhamData(
nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip =
"Electron Scattering Durham Data sample. \n"
"Target: Multiple \n"
"Flux: Energy should match data being handled \n"
"Signal: Any event with an electron in the final state \n";
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.DefineAllowedSpecies("electron");
fSettings.SetTitle("Electron");
fSettings.SetAllowedTypes("FIX/DIAG", "FIX,FREE,SHAPE/DIAG/NORM/MASK");
fSettings.SetXTitle("q0");
fSettings.SetYTitle("#sigma");
fIsNoWidth = true;
FinaliseSampleSettings();
// Plot Setup -------------------------------------------------------
SetDataFromName(fSettings.GetS("originalname"));
SetCovarFromDiagonal();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
// fScaleFactor = ((GetEventHistogram()->Integral("width") * 1E-38 / (fNEvents
// + 0.)) / TotalIntegratedFlux());
EnuMin = fZLowLim;
EnuMax = fZHighLim;
double sigscale = GetEventHistogram()->Integral() * 1E-38 / double(fNEvents) /
TotalIntegratedFlux();
// double dangle = 2 * M_PI * fabs((1. - cos(fYLowLim * M_PI / 180.)) - (1. -
// cos(fYHighLim * M_PI / 180.)));
// fScaleFactor = sigscale / dangle / fZCenter;
fScaleFactor = sigscale;
std::cout << "Event Integral = " << GetEventHistogram()->Integral()
<< std::endl;
std::cout << "Flux Integral = " << TotalIntegratedFlux() << std::endl;
std::cout << "FNEvents = " << fNEvents << std::endl;
std::cout << "Z Limits = " << fZLowLim << " " << fZHighLim << std::endl;
std::cout << "sigscale = " << sigscale << std::endl;
std::cout << "fZCenter = " << fZCenter << std::endl;
std::cout << "ScaleFactor = " << fScaleFactor << std::endl;
// Finish up
FinaliseMeasurement();
};
//********************************************************************
void ElectronScattering_DurhamData::SetDataFromName(std::string name) {
//********************************************************************
// Data Should be given in the format
// Electron_Z_A_Energy_Theta_Source
std::vector splitstring = GeneralUtils::ParseToStr(name, "_");
std::string zstring = splitstring[1];
std::string astring = splitstring[2];
std::string estring = splitstring[3];
std::string tstring = splitstring[4];
std::string sstring = splitstring[5];
fYCenter = GeneralUtils::StrToDbl(tstring);
fZCenter = GeneralUtils::StrToDbl(estring);
// Create effective E and Theta bin Edges
std::vector thetabinedges;
std::vector ebinedges;
int nthetabins = FitPar::Config().GetParI("Electron_NThetaBins");
int nebins = FitPar::Config().GetParI("Electron_NEnergyBins");
double thetawidth = FitPar::Config().GetParD("Electron_ThetaWidth");
double ewidth = FitPar::Config().GetParD("Electron_EnergyWidth");
for (int i = -nthetabins; i <= nthetabins; i++) {
thetabinedges.push_back(fYCenter + thetawidth * (double(i)));
}
for (int i = -nebins; i <= nebins; i++) {
double newval = fZCenter + ewidth * (double(i));
if (newval < 0.0) newval = 0.0;
if (newval < GetEventHistogram()->GetXaxis()->GetXmin())
newval = GetEventHistogram()->GetXaxis()->GetXmin();
if (newval > GetEventHistogram()->GetXaxis()->GetXmax())
newval = GetEventHistogram()->GetXaxis()->GetXmax();
if (std::find(ebinedges.begin(), ebinedges.end(), newval) !=
ebinedges.end())
continue;
ebinedges.push_back(newval);
}
// Determine target
std::string target = "";
if (!zstring.compare("6") && !astring.compare("12"))
target = "12C.dat";
else if (!zstring.compare("8") && !astring.compare("16"))
target = "16O.dat";
else {
ERR(FTL) << "Target not supported in electron scattering module!"
<< std::endl;
throw;
}
// Fill Data Points
std::string line;
std::ifstream mask((FitPar::GetDataBase() + "/Electron/" + target).c_str(),
- ifstream::in);
+ std::ifstream::in);
if (!mask.good()) {
ERR(FTL) << "Failed to open e-scattering database file: "
<< (FitPar::GetDataBase() + "/Electron/" + target) << std::endl;
throw;
}
int i = 0;
std::vector pointx;
std::vector errorx;
std::vector pointy;
std::vector errory;
double scalef = 1.E-38 * 1.E5;
while (std::getline(mask >> std::ws, line, '\n')) {
// std::cout << "Line = " << line << std::endl;
if (line.empty()) continue;
std::vector lineentries = GeneralUtils::ParseToStr(line, " ");
// std::cout << "Checking : " << line << std::endl;
if (zstring.compare(lineentries[0])) continue;
if (astring.compare(lineentries[1])) continue;
if (estring.compare(lineentries[2])) continue;
if (tstring.compare(lineentries[3])) continue;
if (sstring.compare(lineentries[7])) continue;
// std::cout << "Registering data point : " << line << std::endl;
// std::cout << "Adding Graph Point : " <<
// GeneralUtils::StrToDbl(lineentries[4]) << " " <<
// GeneralUtils::StrToDbl(lineentries[5]) << std::endl;
// Loop through x and y points and find a place to insert
if (pointx.empty()) {
pointx.push_back(GeneralUtils::StrToDbl(lineentries[4]));
errorx.push_back(0.0);
pointy.push_back(GeneralUtils::StrToDbl(lineentries[5]) * scalef);
errory.push_back(GeneralUtils::StrToDbl(lineentries[6]) * scalef);
} else {
for (size_t j = 0; j < pointx.size(); j++) {
if (GeneralUtils::StrToDbl(lineentries[4]) < pointx[j] && j == 0) {
// std::cout << "Inserting at start point iterator " << std::endl;
pointx.insert(pointx.begin() + j,
GeneralUtils::StrToDbl(lineentries[4]));
errorx.insert(errorx.begin() + j, 0.0);
pointy.insert(pointy.begin() + j,
GeneralUtils::StrToDbl(lineentries[5]) * scalef);
errory.insert(errory.begin() + j,
GeneralUtils::StrToDbl(lineentries[6]) * scalef);
break;
} else if (GeneralUtils::StrToDbl(lineentries[4]) > pointx[j] &&
j == pointx.size() - 1) {
// std::cout << "Pushing back data point " << std::endl;
pointx.push_back(GeneralUtils::StrToDbl(lineentries[4]));
errorx.push_back(0.0);
pointy.push_back(GeneralUtils::StrToDbl(lineentries[5]) * scalef);
errory.push_back(GeneralUtils::StrToDbl(lineentries[6]) * scalef);
break;
} else if (GeneralUtils::StrToDbl(lineentries[4]) > pointx[j - 1] &&
GeneralUtils::StrToDbl(lineentries[4]) < pointx[j]) {
// std::cout << "Inserting at point iterator = " << j << std::endl;
pointx.insert(pointx.begin() + j,
GeneralUtils::StrToDbl(lineentries[4]));
errorx.insert(errorx.begin() + j, 0.0);
pointy.insert(pointy.begin() + j,
GeneralUtils::StrToDbl(lineentries[5]) * scalef);
errory.insert(errory.begin() + j,
GeneralUtils::StrToDbl(lineentries[6]) * scalef);
break;
}
}
}
// pointx.push_back(GeneralUtils::StrToDbl(lineentries[4]));
// errorx.push_back(0.0);
// pointy.push_back(GeneralUtils::StrToDbl(lineentries[5]));
// errory.push_back(GeneralUtils::StrToDbl(lineentries[6]));
i++;
}
if (!pointx.size() || (pointx.size() != errorx.size()) || !pointy.size() ||
(pointy.size() != errory.size())) {
ERR(FTL) << "Failed to find dataset: " << name << "{"
<< "Z: " << zstring << ", A: " << astring << ", E: " << estring
<< ", CTheta: " << tstring << ", PubID: " << sstring
<< " } in file: "
<< (FitPar::GetDataBase() + "/Electron/" + target) << std::endl;
throw;
}
// for (uint i = 0; i < pointx.size(); i++) {
// std::cout << "Q0 Point " << i << " = " << pointx[i] << std::endl;
// }
fDataGraph = new TGraphErrors(pointx.size(), &pointx[0], &pointy[0],
&errorx[0], &errory[0]);
fDataGraph->SetNameTitle((fName + "_data_GRAPH").c_str(),
(fName + "_data_GRAPH").c_str());
// Now form an effective data and mc histogram
std::vector q0binedges;
const double* x = fDataGraph->GetX();
// Loop over graph and get mid way point between each data point.
for (int i = 0; i < fDataGraph->GetN(); i++) {
// std::cout << "X Point = " << x[i] << std::endl;
if (i == 0) {
// First point set lower width as half distance to point above
q0binedges.push_back(x[0] - ((x[1] - x[0]) / 2.0));
} else if (i == fDataGraph->GetN() - 1) {
// Last point set upper width as half distance to point above.
q0binedges.push_back(x[i] - ((x[i] - x[i - 1]) / 2.0));
q0binedges.push_back(x[i] + ((x[i] - x[i - 1]) / 2.0));
} else {
// Set half distance to point below
q0binedges.push_back(x[i] - ((x[i] - x[i - 1]) / 2.0));
}
}
// Bubble Sort
// for (uint i = 0; i < q0binedges.size(); i++) {
// std::cout << "Q0 Edge " << i << " = " << q0binedges[i] << std::endl;
// }
// for (uint i = 0; i < ebinedges.size(); i++) {
// std::cout << "e Edge " << i << " = " << ebinedges[i] << std::endl;
// }
// for (uint i = 0; i < thetabinedges.size(); i++) {
// std::cout << "theta Edge " << i << " = " << thetabinedges[i] <<
// std::endl;
// }
// Form the data hist, mchist, etc
fDataHist = new TH1D((fName + "_data").c_str(), (fName + "_data").c_str(),
q0binedges.size() - 1, &q0binedges[0]);
fMCHist = (TH1D*)fDataHist->Clone("MC");
const double* y = fDataGraph->GetY();
const double* ey = fDataGraph->GetEY();
for (int i = 0; i < fDataGraph->GetN(); i++) {
// std::cout << "Setting Data Bin " << i + 1 << " to " << y[i] << " +- " <<
// ey[i] << std::endl;
fDataHist->SetBinContent(i + 1, y[i]);
fDataHist->SetBinError(i + 1, ey[i]);
fMCHist->SetBinContent(i + 1, 0.0);
fMCHist->SetBinError(i + 1, 0.0);
}
fMCScan_Q0vsThetavsE =
new TH3D((fName + "_MC_q0vsthetavse").c_str(), "MC_q0vsthetavse",
q0binedges.size() - 1, &q0binedges[0], thetabinedges.size() - 1,
&thetabinedges[0], ebinedges.size() - 1, &ebinedges[0]);
fMCScan_Q0vsThetavsE->Reset();
fMCScan_Q0vsTheta = new TH2D(
(fName + "_MC_q0vstheta").c_str(), "MC_q0vstheta", q0binedges.size() - 1,
&q0binedges[0], thetabinedges.size() - 1, &thetabinedges[0]);
fMCScan_Q0vsTheta->Reset();
fMCScan_Q0vsE =
new TH2D((fName + "_MC_q0vse").c_str(), "MC_q0vse", q0binedges.size() - 1,
&q0binedges[0], ebinedges.size() - 1, &ebinedges[0]);
fMCScan_Q0vsE->Reset();
fXLowLim = fMCScan_Q0vsThetavsE->GetXaxis()->GetBinLowEdge(1);
fXHighLim = fMCScan_Q0vsThetavsE->GetXaxis()->GetBinLowEdge(
fMCScan_Q0vsThetavsE->GetNbinsX() + 2);
fYLowLim = fMCScan_Q0vsThetavsE->GetYaxis()->GetBinLowEdge(1);
fYHighLim = fMCScan_Q0vsThetavsE->GetYaxis()->GetBinLowEdge(
fMCScan_Q0vsThetavsE->GetNbinsY() + 2);
fZLowLim = fMCScan_Q0vsThetavsE->GetZaxis()->GetBinLowEdge(1);
fZHighLim = fMCScan_Q0vsThetavsE->GetZaxis()->GetBinLowEdge(
fMCScan_Q0vsThetavsE->GetNbinsZ() + 2);
std::cout << "Sample " << name << "initialised: "
<< "{" << fXLowLim << "--" << fXHighLim << ", " << fYLowLim << "--"
<< fYHighLim << ", " << fZLowLim << "--" << fZHighLim << "}"
<< std::endl;
}
//********************************************************************
void ElectronScattering_DurhamData::FillEventVariables(FitEvent* event) {
//********************************************************************
if (event->NumFSParticle(11) == 0) return;
FitParticle* ein = event->PartInfo(0);
FitParticle* eout = event->GetHMFSParticle(11);
double q0 = fabs(ein->fP.E() - eout->fP.E()) / 1000.0;
double E = ein->fP.E() / 1000.0;
double theta = ein->fP.Vect().Angle(eout->fP.Vect()) * 180. / M_PI;
fXVar = q0;
fYVar = theta;
fZVar = E;
return;
};
//********************************************************************
bool ElectronScattering_DurhamData::isSignal(FitEvent* event) {
//********************************************************************
if (event->NumFSParticle(11) == 0) {
std::cout << "Ev Cut due to no FS electron." << std::endl;
return false;
}
// std::cout << "fXVar = " << fXVar << " " << fXLowLim << " " << fXHighLim <<
// std::endl;
// std::cout << "fYVar = " << fYVar << " " << fYLowLim << " " << fYHighLim <<
// std::endl;
// std::cout << "fZVar = " << fZVar << " " << fZLowLim << " " << fZHighLim <<
// std::endl;
// std::cout << "iWeight: " << event->InputWeight << std::endl;
if (fXVar < fXLowLim or fXVar > fXHighLim) {
// std::cout << "Ev Cut due to X lims: " << fXLowLim << " -- " << fXHighLim
// << " !<> " << fXVar << std::endl;
return false;
}
if (fYVar < fYLowLim or fYVar > fYHighLim) {
// std::cout << "Ev Cut due to Y lims: " << fYLowLim << " -- " << fYHighLim
// << " !<> " << fXVar << std::endl;
return false;
}
if (fZVar < fZLowLim or fZVar > fZHighLim) {
// std::cout << "Ev Cut due to Z lims: " << fZLowLim << " -- " << fZHighLim
// << " !<> " << fXVar << std::endl;
return false;
}
return true;
};
//********************************************************************
void ElectronScattering_DurhamData::FillHistograms() {
//********************************************************************
Measurement1D::FillHistograms();
if (Signal) {
fMCScan_Q0vsThetavsE->Fill(fXVar, fYVar, fZVar);
fMCScan_Q0vsTheta->Fill(fXVar, fYVar);
fMCScan_Q0vsE->Fill(fXVar, fZVar);
}
}
// Weight = 1.0;
// if (Signal) {
// fMCHist->Fill(fXVar, Weight);
// fMCFine->Fill(fXVar, Weight);
// fMCStat->Fill(fXVar, 1.0);
// if (fMCHist_Modes) fMCHist_Modes->Fill(Mode, fXVar, Weight);
// }
// }
void ElectronScattering_DurhamData::ResetAll() {
Measurement1D::ResetAll();
fMCScan_Q0vsThetavsE->Reset();
fMCScan_Q0vsTheta->Reset();
fMCScan_Q0vsE->Reset();
}
void ElectronScattering_DurhamData::ApplyNormScale(double norm) {
Measurement1D::ApplyNormScale(norm);
fMCScan_Q0vsThetavsE->Scale(1.0 / norm);
fMCScan_Q0vsTheta->Scale(1.0 / norm);
fMCScan_Q0vsE->Scale(1.0 / norm);
}
//********************************************************************
void ElectronScattering_DurhamData::ScaleEvents() {
//********************************************************************
Measurement1D::ScaleEvents();
/*
fMCScan_Q0vsThetavsE->Scale(fScaleFactor, "width");
// Project into fMCScan_Q0vsTheta
for (int x = 0; x < fMCScan_Q0vsThetavsE->GetNbinsX(); x++) {
for (int y = 0; y < fMCScan_Q0vsThetavsE->GetNbinsY(); y++) {
double total = 0.;
for (int z = 0; z < fMCScan_Q0vsThetavsE->GetNbinsZ(); z++) {
double zwidth = fMCScan_Q0vsThetavsE->GetZaxis()->GetBinWidth(z + 1);
total += fMCScan_Q0vsThetavsE->GetBinContent(x + 1, y + 1, z + 1) *
zwidth;
}
fMCScan_Q0vsTheta->SetBinContent(x + 1, y + 1, total);
}
}
// Project into fMCScan_Q0vsE
for (int x = 0; x < fMCScan_Q0vsThetavsE->GetNbinsX(); x++) {
for (int z = 0; z < fMCScan_Q0vsThetavsE->GetNbinsZ(); z++) {
double total = 0.;
for (int y = 0; y < fMCScan_Q0vsThetavsE->GetNbinsY(); y++) {
double ywidth = fMCScan_Q0vsThetavsE->GetYaxis()->GetBinWidth(y + 1);
total += fMCScan_Q0vsThetavsE->GetBinContent(x + 1, y + 1, z + 1) *
ywidth;
}
fMCScan_Q0vsE->SetBinContent(x + 1, z + 1, total);
}
}
// Project fMCScan_Q0vsTheta into MC Hist
for (int x = 0; x < fMCScan_Q0vsTheta->GetNbinsX(); x++) {
double total = 0.;
for (int y = 0; y < fMCScan_Q0vsTheta->GetNbinsY(); y++) {
double ywidth = fMCScan_Q0vsTheta->GetYaxis()->GetBinWidth(y + 1);
total += fMCScan_Q0vsTheta->GetBinContent(x + 1, y + 1);
}
double xwidth = fMCScan_Q0vsTheta->GetXaxis()->GetBinWidth(x + 1);
fMCHist->SetBinContent(x + 1, total * xwidth);
}
fMCHist->Scale(fDataHist->Integral() / fMCHist->Integral());
*/
}
//********************************************************************
int ElectronScattering_DurhamData::GetNDOF() {
//********************************************************************
return fDataGraph->GetN();
}
void ElectronScattering_DurhamData::Write(std::string drawOpts) {
Measurement1D::Write(drawOpts);
fMCScan_Q0vsThetavsE->Write();
fMCScan_Q0vsTheta->Write();
fMCScan_Q0vsE->Write();
fDataGraph->Write();
}
double ElectronScattering_DurhamData::GetLikelihood() { return 0.0; }
void ElectronScattering_DurhamData::SetFitOptions(std::string opt) { return; }
TH1D* ElectronScattering_DurhamData::GetMCHistogram(void) { return fMCHist; }
TH1D* ElectronScattering_DurhamData::GetDataHistogram(void) {
return fDataHist;
}
diff --git a/src/FitBase/JointMeas1D.cxx b/src/FitBase/JointMeas1D.cxx
index 6700a10..6c3e699 100644
--- a/src/FitBase/JointMeas1D.cxx
+++ b/src/FitBase/JointMeas1D.cxx
@@ -1,2303 +1,2303 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This ile 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 "JointMeas1D.h"
//********************************************************************
JointMeas1D::JointMeas1D(void) {
//********************************************************************
// XSec Scalings
fScaleFactor = -1.0;
fCurrentNorm = 1.0;
// Histograms
fDataHist = NULL;
fDataTrue = NULL;
fMCHist = NULL;
fMCFine = NULL;
fMCWeighted = NULL;
fMaskHist = NULL;
// Covar
covar = NULL;
fFullCovar = NULL;
fCovar = NULL;
fInvert = NULL;
fDecomp = NULL;
// Fake Data
fFakeDataInput = "";
fFakeDataFile = NULL;
// Options
fDefaultTypes = "FIX/FULL/CHI2";
fAllowedTypes =
"FIX,FREE,SHAPE/FULL,DIAG/CHI2/NORM/ENUCORR/Q2CORR/ENU1D/MASK";
fIsFix = false;
fIsShape = false;
fIsFree = false;
fIsDiag = false;
fIsFull = false;
fAddNormPen = false;
fIsMask = false;
fIsChi2SVD = false;
fIsRawEvents = false;
fIsDifXSec = false;
fIsEnu1D = false;
// Inputs
fInput = NULL;
fRW = NULL;
// Extra Histograms
fMCHist_Modes = NULL;
for (std::vector::const_iterator iter = fSubChain.begin();
iter != fSubChain.end(); iter++) {
MeasurementBase* exp = *iter;
if (exp) delete exp;
}
fSubChain.clear();
// Flags for Joint Measurements
fIsRatio = false;
fIsSummed = false;
fSaveSubMeas = false;
fIsJoint = true;
}
//********************************************************************
void JointMeas1D::SetupDefaultHist() {
//********************************************************************
// Setup fMCHist
fMCHist = (TH1D*)fDataHist->Clone();
fMCHist->SetNameTitle((fName + "_MC").c_str(),
(fName + "_MC" + fPlotTitles).c_str());
// Setup fMCFine
Int_t nBins = fMCHist->GetNbinsX();
fMCFine = new TH1D(
(fName + "_MC_FINE").c_str(), (fName + "_MC_FINE" + fPlotTitles).c_str(),
nBins * 6, fMCHist->GetBinLowEdge(1), fMCHist->GetBinLowEdge(nBins + 1));
fMCStat = (TH1D*)fMCHist->Clone();
fMCStat->Reset();
fMCHist->Reset();
fMCFine->Reset();
// Setup the NEUT Mode Array
PlotUtils::CreateNeutModeArray((TH1D*)fMCHist, (TH1**)fMCHist_PDG);
PlotUtils::ResetNeutModeArray((TH1**)fMCHist_PDG);
// Setup bin masks using sample name
if (fIsMask) {
std::string maskloc = FitPar::Config().GetParDIR(fName + ".mask");
if (maskloc.empty()) {
maskloc = FitPar::GetDataBase() + "/masks/" + fName + ".mask";
}
SetBinMask(maskloc);
}
fMCHist_Modes = new TrueModeStack( (fName + "_MODES").c_str(), ("True Channels"), fMCHist);
SetAutoProcessTH1(fMCHist_Modes);
return;
}
//********************************************************************
JointMeas1D::~JointMeas1D(void) {
//********************************************************************
if (fDataHist) delete fDataHist;
if (fDataTrue) delete fDataTrue;
if (fMCHist) delete fMCHist;
if (fMCFine) delete fMCFine;
if (fMCWeighted) delete fMCWeighted;
if (fMaskHist) delete fMaskHist;
if (covar) delete covar;
if (fFullCovar) delete fFullCovar;
if (fCovar) delete fCovar;
if (fInvert) delete fInvert;
if (fDecomp) delete fDecomp;
for (std::vector::const_iterator iter = fSubChain.begin();
iter != fSubChain.end(); iter++) {
MeasurementBase* exp = *iter;
if (exp) delete exp;
}
fSubChain.clear();
}
//********************************************************************
SampleSettings JointMeas1D::LoadSampleSettings(nuiskey samplekey){
//********************************************************************
SampleSettings s = MeasurementBase::LoadSampleSettings(samplekey);
// Parse Inputs
fSubInFiles.clear();
std::vector entries = GeneralUtils::ParseToStr(s.GetS("input"), ";");
if (entries.size() < 2) {
ERR(FTL) << "Joint measurement expected to recieve at least two semi-colon "
"separated input files, but recieved: \""
<< s.GetS("input") << "\"" << std::endl;
throw;
}
std::vector first_file_descriptor =
GeneralUtils::ParseToStr(entries.front(), ":");
if (first_file_descriptor.size() != 2) {
ERR(FTL) << "Found Joint measurement where the input file had no type: \""
<< s.GetS("input") << "\", expected \"INPUTTYPE:File.root;File2.root\"."
<< std::endl;
throw;
}
std::string inpType = first_file_descriptor[0];
for (std::vector::iterator iter = entries.begin();
iter != entries.end(); iter++) {
if (GeneralUtils::ParseToStr(*iter, ":").size() != 2) {
std::stringstream ss("");
ss << inpType << ":" << (*iter);
fSubInFiles.push_back(ss.str());
} else {
fSubInFiles.push_back(*iter);
}
}
return s;
}
//********************************************************************
void JointMeas1D::FinaliseSampleSettings() {
//********************************************************************
// Setup naming + renaming
fName = fSettings.GetName();
fSettings.SetS("originalname", fName);
if (fSettings.Has("rename")) {
fName = fSettings.GetS("rename");
fSettings.SetS("name", fName);
}
// Setup all other options
LOG(SAM) << "Finalising Sample Settings: " << fName << std::endl;
if ((fSettings.GetS("originalname").find("Evt") != std::string::npos)) {
fIsRawEvents = true;
LOG(SAM) << "Found event rate measurement but using poisson likelihoods."
<< std::endl;
}
if (fSettings.GetS("originalname").find("XSec_1DEnu") != std::string::npos) {
fIsEnu1D = true;
LOG(SAM) << "::" << fName << "::" << std::endl;
LOG(SAM) << "Found XSec Enu measurement, applying flux integrated scaling, "
<< "not flux averaged!" << std::endl;
}
if (fIsEnu1D && fIsRawEvents) {
LOG(SAM) << "Found 1D Enu XSec distribution AND fIsRawEvents, is this "
"really correct?!"
<< std::endl;
LOG(SAM) << "Check experiment constructor for " << fName
<< " and correct this!" << std::endl;
LOG(SAM) << "I live in " << __FILE__ << ":" << __LINE__ << std::endl;
exit(-1);
}
// Parse Inputs
fSubInFiles.clear();
std::vector entries = GeneralUtils::ParseToStr(fSettings.GetS("input"), ";");
if (entries.size() < 2) {
ERR(FTL) << "Joint measurement expected to recieve at least two semi-colon "
"separated input files, but recieved: \""
<< fSettings.GetS("input") << "\"" << std::endl;
throw;
}
std::vector first_file_descriptor =
GeneralUtils::ParseToStr(entries.front(), ":");
if (first_file_descriptor.size() != 2) {
ERR(FTL) << "Found Joint measurement where the input file had no type: \""
<< fSettings.GetS("input") << "\", expected \"INPUTTYPE:File.root;File2.root\"."
<< std::endl;
throw;
}
std::string inpType = first_file_descriptor[0];
for (std::vector::iterator iter = entries.begin();
iter != entries.end(); iter++) {
if (GeneralUtils::ParseToStr(*iter, ":").size() != 2) {
std::stringstream ss("");
ss << inpType << ":" << (*iter);
fSubInFiles.push_back(ss.str());
} else {
fSubInFiles.push_back(*iter);
}
}
// Setup options
SetFitOptions(fDefaultTypes); // defaults
SetFitOptions(fSettings.GetS("type")); // user specified
EnuMin = GeneralUtils::StrToDbl(fSettings.GetS("enu_min"));
EnuMax = GeneralUtils::StrToDbl(fSettings.GetS("enu_max"));
if (fAddNormPen) {
if (fNormError <= 0.0) {
ERR(WRN) << "Norm error for class " << fName << " is 0.0!" << std::endl;
ERR(WRN) << "If you want to use it please add fNormError=VAL" << std::endl;
throw;
}
}
if (!fRW) fRW = FitBase::GetRW();
LOG(SAM) << "Finalised Sample Settings" << std::endl;
}
//********************************************************************
void JointMeas1D::SetDataFromTextFile(std::string datafile) {
//********************************************************************
LOG(SAM) << "Reading data from text file: " << datafile << std::endl;
fDataHist = PlotUtils::GetTH1DFromFile(datafile,
fSettings.GetName() + "_data",
fSettings.GetFullTitles());
}
//********************************************************************
void JointMeas1D::SetDataFromRootFile(std::string datafile,
std::string histname) {
//********************************************************************
LOG(SAM) << "Reading data from root file: " << datafile << ";" << histname << std::endl;
fDataHist = PlotUtils::GetTH1DFromRootFile(datafile, histname);
fDataHist->SetNameTitle((fSettings.GetName() + "_data").c_str(),
(fSettings.GetFullTitles()).c_str());
return;
};
//********************************************************************
void JointMeas1D::SetPoissonErrors() {
//********************************************************************
if (!fDataHist) {
ERR(FTL) << "Need a data hist to setup possion errors! " << std::endl;
ERR(FTL) << "Setup Data First!" << std::endl;
throw;
}
for (int i = 0; i < fDataHist->GetNbinsX() + 1; i++) {
fDataHist->SetBinError(i + 1, sqrt(fDataHist->GetBinContent(i + 1)));
}
}
//********************************************************************
void JointMeas1D::SetCovarFromDiagonal(TH1D* data) {
//********************************************************************
if (!data and fDataHist) {
data = fDataHist;
}
if (data) {
LOG(SAM) << "Setting diagonal covariance for: " << data->GetName() << std::endl;
fFullCovar = StatUtils::MakeDiagonalCovarMatrix(data);
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
} else {
ERR(FTL) << "No data input provided to set diagonal covar from!" << std::endl;
}
if (!fIsDiag) {
ERR(FTL) << "SetCovarMatrixFromDiag called for measurement "
<< "that is not set as diagonal." << std::endl;
throw;
}
}
//********************************************************************
void JointMeas1D::SetCovarFromTextFile(std::string covfile, int dim) {
//********************************************************************
LOG(SAM) << "Reading covariance from text file: " << covfile << std::endl;
fFullCovar = StatUtils::GetCovarFromTextFile(covfile, dim);
-
+
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
//********************************************************************
void JointMeas1D::SetCovarFromMultipleTextFiles(std::string covfiles, int dim) {
//********************************************************************
if (dim == -1) {
dim = fDataHist->GetNbinsX();
}
std::vector covList = GeneralUtils::ParseToStr(covfiles, ";");
fFullCovar = new TMatrixDSym(dim);
for (uint i = 0; i < covList.size(); ++i){
LOG(SAM) << "Reading covariance from text file: " << covList[i] << std::endl;
TMatrixDSym* temp_cov = StatUtils::GetCovarFromTextFile(covList[i], dim);
(*fFullCovar) += (*temp_cov);
delete temp_cov;
}
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
//********************************************************************
void JointMeas1D::SetCovarFromRootFile(std::string covfile, std::string histname) {
//********************************************************************
LOG(SAM) << "Reading covariance from text file: " << covfile << ";" << histname << std::endl;
fFullCovar = StatUtils::GetCovarFromRootFile(covfile, histname);
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
//********************************************************************
void JointMeas1D::SetCovarInvertFromTextFile(std::string covfile, int dim) {
//********************************************************************
LOG(SAM) << "Reading inverted covariance from text file: " << covfile << std::endl;
covar = StatUtils::GetCovarFromTextFile(covfile, dim);
fFullCovar = StatUtils::GetInvert(covar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
//********************************************************************
void JointMeas1D::SetCovarInvertFromRootFile(std::string covfile, std::string histname) {
//********************************************************************
LOG(SAM) << "Reading inverted covariance from text file: " << covfile << ";" << histname << std::endl;
covar = StatUtils::GetCovarFromRootFile(covfile, histname);
fFullCovar = StatUtils::GetInvert(covar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
//********************************************************************
void JointMeas1D::SetCorrelationFromTextFile(std::string covfile, int dim) {
//********************************************************************
if (dim == -1) dim = fDataHist->GetNbinsX();
LOG(SAM) << "Reading data correlations from text file: " << covfile << ";" << dim << std::endl;
TMatrixDSym* correlation = StatUtils::GetCovarFromTextFile(covfile, dim);
if (!fDataHist) {
ERR(FTL) << "Trying to set correlations from text file but there is no data to build it from. \n"
<< "In constructor make sure data is set before SetCorrelationFromTextFile is called. \n" << std::endl;
throw;
}
// Fill covar from data errors and correlations
fFullCovar = new TMatrixDSym(dim);
for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
for (int j = 0; j < fDataHist->GetNbinsX(); j++) {
(*fFullCovar)(i, j) = (*correlation)(i, j) * fDataHist->GetBinError(i + 1) * fDataHist->GetBinError(j + 1) * 1.E76;
}
}
// Fill other covars.
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
delete correlation;
}
//********************************************************************
void JointMeas1D::SetCorrelationFromMultipleTextFiles(std::string corrfiles, int dim) {
//********************************************************************
if (dim == -1) {
dim = fDataHist->GetNbinsX();
}
std::vector corrList = GeneralUtils::ParseToStr(corrfiles, ";");
fFullCovar = new TMatrixDSym(dim);
for (uint i = 0; i < corrList.size(); ++i){
LOG(SAM) << "Reading covariance from text file: " << corrList[i] << std::endl;
TMatrixDSym* temp_cov = StatUtils::GetCovarFromTextFile(corrList[i], dim);
for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
for (int j = 0; j < fDataHist->GetNbinsX(); j++) {
// Note that there is a factor of 1E76 here. It is very silly indeed.
// However, if you remove it, you also need to fix the factors of 1E38 added to the chi2 calculations!
(*temp_cov)(i, j) = (*temp_cov)(i, j) * fDataHist->GetBinError(i + 1) * fDataHist->GetBinError(j + 1) * 1E76;
}
}
(*fFullCovar) += (*temp_cov);
delete temp_cov;
}
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
//********************************************************************
void JointMeas1D::SetCorrelationFromRootFile(std::string covfile, std::string histname) {
//********************************************************************
LOG(SAM) << "Reading data correlations from text file: " << covfile << ";" << histname << std::endl;
TMatrixDSym* correlation = StatUtils::GetCovarFromRootFile(covfile, histname);
if (!fDataHist) {
ERR(FTL) << "Trying to set correlations from text file but there is no data to build it from. \n"
<< "In constructor make sure data is set before SetCorrelationFromTextFile is called. \n" << std::endl;
throw;
}
// Fill covar from data errors and correlations
fFullCovar = new TMatrixDSym(fDataHist->GetNbinsX());
for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
for (int j = 0; j < fDataHist->GetNbinsX(); j++) {
(*fFullCovar)(i, j) = (*correlation)(i, j) * fDataHist->GetBinError(i + 1) * fDataHist->GetBinError(j + 1) * 1.E76;
}
}
// Fill other covars.
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
delete correlation;
}
void JointMeas1D::SetShapeCovar(){
// Return if this is missing any pre-requisites
if (!fFullCovar) return;
if (!fDataHist) return;
// Also return if it's bloody stupid under the circumstances
if (fIsDiag) return;
fShapeCovar = StatUtils::ExtractShapeOnlyCovar(fFullCovar, fDataHist);
return;
}
//********************************************************************
void JointMeas1D::SetCholDecompFromTextFile(std::string covfile, int dim) {
//********************************************************************
LOG(SAM) << "Reading cholesky from text file: " << covfile << std::endl;
TMatrixD* temp = StatUtils::GetMatrixFromTextFile(covfile, dim, dim);
TMatrixD* trans = (TMatrixD*)temp->Clone();
trans->T();
(*trans) *= (*temp);
fFullCovar = new TMatrixDSym(dim, trans->GetMatrixArray(), "");
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
delete temp;
delete trans;
}
//********************************************************************
void JointMeas1D::SetCholDecompFromRootFile(std::string covfile, std::string histname) {
//********************************************************************
LOG(SAM) << "Reading cholesky decomp from root file: " << covfile << ";" << histname << std::endl;
TMatrixD* temp = StatUtils::GetMatrixFromRootFile(covfile, histname);
TMatrixD* trans = (TMatrixD*)temp->Clone();
trans->T();
(*trans) *= (*temp);
fFullCovar = new TMatrixDSym(temp->GetNrows(), trans->GetMatrixArray(), "");
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
delete temp;
delete trans;
}
//********************************************************************
void JointMeas1D::ScaleData(double scale) {
//********************************************************************
fDataHist->Scale(scale);
}
//********************************************************************
void JointMeas1D::ScaleCovar(double scale) {
//********************************************************************
(*fFullCovar) *= scale;
(*covar) *= 1.0 / scale;
(*fDecomp) *= sqrt(scale);
}
//********************************************************************
void JointMeas1D::SetBinMask(std::string maskfile) {
//********************************************************************
if (!fIsMask) return;
LOG(SAM) << "Reading bin mask from file: " << maskfile << std::endl;
// Create a mask histogram with dim of data
int nbins = fDataHist->GetNbinsX();
fMaskHist =
new TH1I((fSettings.GetName() + "_BINMASK").c_str(),
(fSettings.GetName() + "_BINMASK; Bin; Mask?").c_str(), nbins, 0, nbins);
std::string line;
- std::ifstream mask(maskfile.c_str(), ifstream::in);
+ std::ifstream mask(maskfile.c_str(), std::ifstream::in);
if (!mask.is_open()) {
LOG(FTL) << " Cannot find mask file." << std::endl;
throw;
}
while (std::getline(mask >> std::ws, line, '\n')) {
std::vector entries = GeneralUtils::ParseToInt(line, " ");
// Skip lines with poorly formatted lines
if (entries.size() < 2) {
LOG(WRN) << "JointMeas1D::SetBinMask(), couldn't parse line: " << line
<< std::endl;
continue;
}
// The first index should be the bin number, the second should be the mask
// value.
int val = 0;
if (entries[1] > 0) val = 1;
fMaskHist->SetBinContent(entries[0], val);
}
// Apply masking by setting masked data bins to zero
PlotUtils::MaskBins(fDataHist, fMaskHist);
return;
}
//********************************************************************
void JointMeas1D::FinaliseMeasurement() {
//********************************************************************
LOG(SAM) << "Finalising Measurement: " << fName << std::endl;
// Make sure data is setup
if (!fDataHist) {
ERR(FTL) << "No data has been setup inside " << fName << " constructor!" << std::endl;
throw;
}
// Make sure covariances are setup
if (!fFullCovar) {
fIsDiag = true;
SetCovarFromDiagonal(fDataHist);
}
if (!covar) {
covar = StatUtils::GetInvert(fFullCovar);
}
if (!fDecomp) {
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
// Push the diagonals of fFullCovar onto the data histogram
// Comment out until scaling is used consistently...
StatUtils::SetDataErrorFromCov(fDataHist, fFullCovar, 1E-38);
// Setup fMCHist from data
fMCHist = (TH1D*)fDataHist->Clone();
fMCHist->SetNameTitle((fSettings.GetName() + "_MC").c_str(),
(fSettings.GetFullTitles()).c_str());
fMCHist->Reset();
// Setup fMCFine
fMCFine = new TH1D("mcfine", "mcfine", fDataHist->GetNbinsX(),
fMCHist->GetBinLowEdge(1),
fMCHist->GetBinLowEdge(fDataHist->GetNbinsX() + 1));
fMCFine->SetNameTitle((fSettings.GetName() + "_MC_FINE").c_str(),
(fSettings.GetFullTitles()).c_str());
fMCFine->Reset();
// Setup MC Stat
fMCStat = (TH1D*)fMCHist->Clone();
fMCStat->Reset();
// Search drawopts for possible types to include by default
std::string drawopts = FitPar::Config().GetParS("drawopts");
if (drawopts.find("MODES") != std::string::npos) {
fMCHist_Modes = new TrueModeStack( (fSettings.GetName() + "_MODES").c_str(),
("True Channels"), fMCHist);
SetAutoProcessTH1(fMCHist_Modes);
}
// Setup bin masks using sample name
if (fIsMask) {
std::string curname = fName;
std::string origname = fSettings.GetS("originalname");
// Check rename.mask
std::string maskloc = FitPar::Config().GetParDIR(curname + ".mask");
// Check origname.mask
if (maskloc.empty()) maskloc = FitPar::Config().GetParDIR(origname + ".mask");
// Check database
if (maskloc.empty()) {
maskloc = FitPar::GetDataBase() + "/masks/" + origname + ".mask";
}
// Setup Bin Mask
SetBinMask(maskloc);
}
/*
if (fScaleFactor < 0) {
ERR(FTL) << "I found a negative fScaleFactor in " << __FILE__ << ":" << __LINE__ << std::endl;
ERR(FTL) << "fScaleFactor = " << fScaleFactor << std::endl;
ERR(FTL) << "EXITING" << std::endl;
throw;
}
*/
// Create and fill Weighted Histogram
if (!fMCWeighted) {
fMCWeighted = (TH1D*)fMCHist->Clone();
fMCWeighted->SetNameTitle((fName + "_MCWGHTS").c_str(),
(fName + "_MCWGHTS" + fPlotTitles).c_str());
fMCWeighted->GetYaxis()->SetTitle("Weighted Events");
}
}
//********************************************************************
void JointMeas1D::SetFitOptions(std::string opt) {
//********************************************************************
// Do nothing if default given
if (opt == "DEFAULT") return;
// CHECK Conflicting Fit Options
std::vector fit_option_allow =
GeneralUtils::ParseToStr(fAllowedTypes, "/");
for (UInt_t i = 0; i < fit_option_allow.size(); i++) {
std::vector fit_option_section =
GeneralUtils::ParseToStr(fit_option_allow.at(i), ",");
bool found_option = false;
for (UInt_t j = 0; j < fit_option_section.size(); j++) {
std::string av_opt = fit_option_section.at(j);
if (!found_option and opt.find(av_opt) != std::string::npos) {
found_option = true;
} else if (found_option and opt.find(av_opt) != std::string::npos) {
ERR(FTL) << "ERROR: Conflicting fit options provided: "
<< opt << std::endl
<< "Conflicting group = " << fit_option_section.at(i) << std::endl
<< "You should only supply one of these options in card file." << std::endl;
throw;
}
}
}
// Check all options are allowed
std::vector fit_options_input =
GeneralUtils::ParseToStr(opt, "/");
for (UInt_t i = 0; i < fit_options_input.size(); i++) {
if (fAllowedTypes.find(fit_options_input.at(i)) == std::string::npos) {
ERR(FTL) << "ERROR: Fit Option '" << fit_options_input.at(i)
<< "' Provided is not allowed for this measurement."
<< std::endl;
ERR(FTL) << "Fit Options should be provided as a '/' seperated list "
"(e.g. FREE/DIAG/NORM)"
<< std::endl;
ERR(FTL) << "Available options for " << fName << " are '" << fAllowedTypes
<< "'" << std::endl;
throw;
}
}
// Set TYPE
fFitType = opt;
// FIX,SHAPE,FREE
if (opt.find("FIX") != std::string::npos) {
fIsFree = fIsShape = false;
fIsFix = true;
} else if (opt.find("SHAPE") != std::string::npos) {
fIsFree = fIsFix = false;
fIsShape = true;
} else if (opt.find("FREE") != std::string::npos) {
fIsFix = fIsShape = false;
fIsFree = true;
}
// DIAG,FULL (or default to full)
if (opt.find("DIAG") != std::string::npos) {
fIsDiag = true;
fIsFull = false;
} else if (opt.find("FULL") != std::string::npos) {
fIsDiag = false;
fIsFull = true;
}
// CHI2/LL (OTHERS?)
if (opt.find("LOG") != std::string::npos) {
fIsChi2 = false;
ERR(FTL) << "No other LIKELIHOODS properly supported!" << std::endl;
ERR(FTL) << "Try to use a chi2!" << std::endl;
throw;
} else {
fIsChi2 = true;
}
// EXTRAS
if (opt.find("RAW") != std::string::npos) fIsRawEvents = true;
if (opt.find("DIF") != std::string::npos) fIsDifXSec = true;
if (opt.find("ENU1D") != std::string::npos) fIsEnu1D = true;
if (opt.find("NORM") != std::string::npos) fAddNormPen = true;
if (opt.find("MASK") != std::string::npos) fIsMask = true;
return;
};
//********************************************************************
void JointMeas1D::SetSmearingMatrix(std::string smearfile, int truedim,
int recodim) {
//********************************************************************
// The smearing matrix describes the migration from true bins (rows) to reco
// bins (columns)
// Counter over the true bins!
int row = 0;
std::string line;
- std::ifstream smear(smearfile.c_str(), ifstream::in);
+ std::ifstream smear(smearfile.c_str(), std::ifstream::in);
// Note that the smearing matrix may be rectangular.
fSmearMatrix = new TMatrixD(truedim, recodim);
if (smear.is_open())
LOG(SAM) << "Reading smearing matrix from file: " << smearfile << std::endl;
else
ERR(FTL) << "Smearing matrix provided is incorrect: " << smearfile
<< std::endl;
while (std::getline(smear >> std::ws, line, '\n')) {
int column = 0;
std::vector entries = GeneralUtils::ParseToDbl(line, " ");
for (std::vector::iterator iter = entries.begin();
iter != entries.end(); iter++) {
(*fSmearMatrix)(row, column) =
(*iter) / 100.; // Convert to fraction from
// percentage (this may not be
// general enough)
column++;
}
row++;
}
return;
}
//********************************************************************
void JointMeas1D::ApplySmearingMatrix() {
//********************************************************************
if (!fSmearMatrix) {
ERR(WRN) << fName
<< ": attempted to apply smearing matrix, but none was set"
<< std::endl;
return;
}
TH1D* unsmeared = (TH1D*)fMCHist->Clone();
TH1D* smeared = (TH1D*)fMCHist->Clone();
smeared->Reset();
// Loop over reconstructed bins
// true = row; reco = column
for (int rbin = 0; rbin < fSmearMatrix->GetNcols(); ++rbin) {
// Sum up the constributions from all true bins
double rBinVal = 0;
// Loop over true bins
for (int tbin = 0; tbin < fSmearMatrix->GetNrows(); ++tbin) {
rBinVal +=
(*fSmearMatrix)(tbin, rbin) * unsmeared->GetBinContent(tbin + 1);
}
smeared->SetBinContent(rbin + 1, rBinVal);
}
fMCHist = (TH1D*)smeared->Clone();
return;
}
/*
Reconfigure LOOP
*/
//********************************************************************
void JointMeas1D::ResetAll() {
//********************************************************************
fMCHist->Reset();
fMCFine->Reset();
fMCStat->Reset();
return;
};
//********************************************************************
void JointMeas1D::FillHistograms() {
//********************************************************************
if (Signal) {
fMCHist->Fill(fXVar, Weight);
fMCFine->Fill(fXVar, Weight);
fMCStat->Fill(fXVar, 1.0);
if (fMCHist_Modes) fMCHist_Modes->Fill(Mode, fXVar, Weight);
}
return;
};
//********************************************************************
void JointMeas1D::ScaleEvents() {
//********************************************************************
LOG(FIT) << "Scaling JointMeas1D" << std::endl;
// Fill MCWeighted;
for (int i = 0; i < fMCHist->GetNbinsX(); i++) {
fMCWeighted->SetBinContent(i + 1, fMCHist->GetBinContent(i + 1));
fMCWeighted->SetBinError(i + 1, fMCHist->GetBinError(i + 1));
}
// Setup Stat ratios for MC and MC Fine
double* statratio = new double[fMCHist->GetNbinsX()];
for (int i = 0; i < fMCHist->GetNbinsX(); i++) {
if (fMCHist->GetBinContent(i + 1) != 0) {
statratio[i] = fMCHist->GetBinError(i + 1) / fMCHist->GetBinContent(i + 1);
} else {
statratio[i] = 0.0;
}
}
double* statratiofine = new double[fMCFine->GetNbinsX()];
for (int i = 0; i < fMCFine->GetNbinsX(); i++) {
if (fMCFine->GetBinContent(i + 1) != 0) {
statratiofine[i] = fMCFine->GetBinError(i + 1) / fMCFine->GetBinContent(i + 1);
} else {
statratiofine[i] = 0.0;
}
}
// Scaling for raw event rates
if (fIsRawEvents) {
double datamcratio = fDataHist->Integral() / fMCHist->Integral();
fMCHist->Scale(datamcratio);
fMCFine->Scale(datamcratio);
if (fMCHist_Modes) fMCHist_Modes->Scale(datamcratio);
// Scaling for XSec as function of Enu
} else if (fIsEnu1D) {
PlotUtils::FluxUnfoldedScaling(fMCHist, GetFluxHistogram(),
GetEventHistogram(), fScaleFactor,
fNEvents);
PlotUtils::FluxUnfoldedScaling(fMCFine, GetFluxHistogram(),
GetEventHistogram(), fScaleFactor,
fNEvents);
// if (fMCHist_Modes) {
// PlotUtils::FluxUnfoldedScaling(fMCHist_Modes, GetFluxHistogram(),
// GetEventHistogram(), fScaleFactor,
// fNEvents);
// }
// Any other differential scaling
} else {
fMCHist->Scale(fScaleFactor, "width");
fMCFine->Scale(fScaleFactor, "width");
if (fMCHist_Modes) fMCHist_Modes->Scale(fScaleFactor, "width");
}
// Proper error scaling - ROOT Freaks out with xsec weights sometimes
for (int i = 0; i < fMCStat->GetNbinsX(); i++) {
fMCHist->SetBinError(i + 1, fMCHist->GetBinContent(i + 1) * statratio[i]);
}
for (int i = 0; i < fMCFine->GetNbinsX(); i++) {
fMCFine->SetBinError(i + 1, fMCFine->GetBinContent(i + 1) * statratiofine[i]);
}
// Clean up
delete statratio;
delete statratiofine;
return;
};
//********************************************************************
void JointMeas1D::ApplyNormScale(double norm) {
//********************************************************************
fCurrentNorm = norm;
fMCHist->Scale(1.0 / norm);
fMCFine->Scale(1.0 / norm);
return;
};
/*
Statistic Functions - Outsources to StatUtils
*/
//********************************************************************
int JointMeas1D::GetNDOF() {
//********************************************************************
int ndof = fDataHist->GetNbinsX();
if (fMaskHist) ndof -= fMaskHist->Integral();
return ndof;
}
//********************************************************************
double JointMeas1D::GetLikelihood() {
//********************************************************************
// If this is for a ratio, there is no data histogram to compare to!
if (fNoData || !fDataHist) return 0.;
// Apply Masking to MC if Required.
if (fIsMask and fMaskHist) {
PlotUtils::MaskBins(fMCHist, fMaskHist);
}
// Sort Shape Scaling
double scaleF = 0.0;
if (fIsShape) {
if (fMCHist->Integral(1, fMCHist->GetNbinsX(), "width")) {
scaleF = fDataHist->Integral(1, fDataHist->GetNbinsX(), "width") /
fMCHist->Integral(1, fMCHist->GetNbinsX(), "width");
fMCHist->Scale(scaleF);
fMCFine->Scale(scaleF);
}
}
// Likelihood Calculation
double stat = 0.;
if (fIsChi2) {
if (fIsRawEvents) {
stat = StatUtils::GetChi2FromEventRate(fDataHist, fMCHist, fMaskHist);
} else if (fIsDiag) {
stat = StatUtils::GetChi2FromDiag(fDataHist, fMCHist, fMaskHist);
} else if (!fIsDiag and !fIsRawEvents) {
stat = StatUtils::GetChi2FromCov(fDataHist, fMCHist, covar, fMaskHist);
}
}
// Sort Penalty Terms
if (fAddNormPen) {
double penalty =
(1. - fCurrentNorm) * (1. - fCurrentNorm) / (fNormError * fNormError);
stat += penalty;
}
// Return to normal scaling
if (fIsShape and !FitPar::Config().GetParB("saveshapescaling")) {
fMCHist->Scale(1. / scaleF);
fMCFine->Scale(1. / scaleF);
}
fLikelihood = stat;
return stat;
}
/*
Fake Data Functions
*/
//********************************************************************
void JointMeas1D::SetFakeDataValues(std::string fakeOption) {
//********************************************************************
// Setup original/datatrue
TH1D* tempdata = (TH1D*) fDataHist->Clone();
if (!fIsFakeData) {
fIsFakeData = true;
// Make a copy of the original data histogram.
if (!fDataOrig) fDataOrig = (TH1D*)fDataHist->Clone((fName + "_data_original").c_str());
} else {
ResetFakeData();
}
// Setup Inputs
fFakeDataInput = fakeOption;
LOG(SAM) << "Setting fake data from : " << fFakeDataInput << std::endl;
// From MC
if (fFakeDataInput.compare("MC") == 0) {
fDataHist = (TH1D*)fMCHist->Clone((fName + "_MC").c_str());
// Fake File
} else {
if (!fFakeDataFile) fFakeDataFile = new TFile(fFakeDataInput.c_str(), "READ");
fDataHist = (TH1D*)fFakeDataFile->Get((fName + "_MC").c_str());
}
// Setup Data Hist
fDataHist->SetNameTitle((fName + "_FAKE").c_str(),
(fName + fPlotTitles).c_str());
// Replace Data True
if (fDataTrue) delete fDataTrue;
fDataTrue = (TH1D*)fDataHist->Clone();
fDataTrue->SetNameTitle((fName + "_FAKE_TRUE").c_str(),
(fName + fPlotTitles).c_str());
// Make a new covariance for fake data hist.
int nbins = fDataHist->GetNbinsX();
double alpha_i = 0.0;
double alpha_j = 0.0;
for (int i = 0; i < nbins; i++) {
for (int j = 0; j < nbins; j++) {
alpha_i = fDataHist->GetBinContent(i + 1) / tempdata->GetBinContent(i + 1);
alpha_j = fDataHist->GetBinContent(j + 1) / tempdata->GetBinContent(j + 1);
(*fFullCovar)(i, j) = alpha_i * alpha_j * (*fFullCovar)(i, j);
}
}
// Setup Covariances
if (covar) delete covar;
covar = StatUtils::GetInvert(fFullCovar);
if (fDecomp) delete fDecomp;
fDecomp = StatUtils::GetInvert(fFullCovar);
delete tempdata;
return;
};
//********************************************************************
void JointMeas1D::ResetFakeData() {
//********************************************************************
if (fIsFakeData) {
if (fDataHist) delete fDataHist;
fDataHist = (TH1D*)fDataTrue->Clone((fSettings.GetName() + "_FKDAT").c_str());
}
}
//********************************************************************
void JointMeas1D::ResetData() {
//********************************************************************
if (fIsFakeData) {
if (fDataHist) delete fDataHist;
fDataHist = (TH1D*)fDataOrig->Clone((fSettings.GetName() + "_data").c_str());
}
fIsFakeData = false;
}
//********************************************************************
void JointMeas1D::ThrowCovariance() {
//********************************************************************
// Take a fDecomposition and use it to throw the current dataset.
// Requires fDataTrue also be set incase used repeatedly.
if (fDataHist) delete fDataHist;
fDataHist = StatUtils::ThrowHistogram(fDataTrue, fFullCovar);
return;
};
//********************************************************************
void JointMeas1D::ThrowDataToy(){
//********************************************************************
if (!fDataTrue) fDataTrue = (TH1D*) fDataHist->Clone();
if (fMCHist) delete fMCHist;
fMCHist = StatUtils::ThrowHistogram(fDataTrue, fFullCovar);
}
/*
Access Functions
*/
//********************************************************************
TH1D* JointMeas1D::GetMCHistogram() {
//********************************************************************
if (!fMCHist) return fMCHist;
std::ostringstream chi2;
chi2 << std::setprecision(5) << this->GetLikelihood();
int linecolor = kRed;
int linestyle = 1;
int linewidth = 1;
int fillcolor = 0;
int fillstyle = 1001;
if (fSettings.Has("linecolor")) linecolor = fSettings.GetI("linecolor");
if (fSettings.Has("linestyle")) linestyle = fSettings.GetI("linestyle");
if (fSettings.Has("linewidth")) linewidth = fSettings.GetI("linewidth");
if (fSettings.Has("fillcolor")) fillcolor = fSettings.GetI("fillcolor");
if (fSettings.Has("fillstyle")) fillstyle = fSettings.GetI("fillstyle");
fMCHist->SetTitle(chi2.str().c_str());
fMCHist->SetLineColor(linecolor);
fMCHist->SetLineStyle(linestyle);
fMCHist->SetLineWidth(linewidth);
fMCHist->SetFillColor(fillcolor);
fMCHist->SetFillStyle(fillstyle);
return fMCHist;
};
//********************************************************************
TH1D* JointMeas1D::GetDataHistogram() {
//********************************************************************
if (!fDataHist) return fDataHist;
int datacolor = kBlack;
int datastyle = 1;
int datawidth = 1;
if (fSettings.Has("datacolor")) datacolor = fSettings.GetI("datacolor");
if (fSettings.Has("datastyle")) datastyle = fSettings.GetI("datastyle");
if (fSettings.Has("datawidth")) datawidth = fSettings.GetI("datawidth");
fDataHist->SetLineColor(datacolor);
fDataHist->SetLineWidth(datawidth);
fDataHist->SetMarkerStyle(datastyle);
return fDataHist;
};
/*
Write Functions
*/
// Save all the histograms at once
//********************************************************************
void JointMeas1D::Write(std::string drawOpt) {
//********************************************************************
// Get Draw Options
drawOpt = FitPar::Config().GetParS("drawopts");
// Write Settigns
if (drawOpt.find("SETTINGS") != std::string::npos){
fSettings.Set("#chi^{2}",fLikelihood);
fSettings.Set("NDOF", this->GetNDOF() );
fSettings.Set("#chi^{2}/NDOF", fLikelihood / this->GetNDOF() );
fSettings.Write();
}
// Write Data/MC
GetDataHistogram()->Write();
GetMCHistogram()->Write();
// Write Fine Histogram
if (drawOpt.find("FINE") != std::string::npos)
GetFineList().at(0)->Write();
// Write Weighted Histogram
if (drawOpt.find("WEIGHTS") != std::string::npos && fMCWeighted)
fMCWeighted->Write();
// Save Flux/Evt if no event manager
if (!FitPar::Config().GetParB("EventManager")) {
if (drawOpt.find("FLUX") != std::string::npos && GetFluxHistogram())
GetFluxHistogram()->Write();
if (drawOpt.find("EVT") != std::string::npos && GetEventHistogram())
GetEventHistogram()->Write();
if (drawOpt.find("XSEC") != std::string::npos && GetEventHistogram())
GetEventHistogram()->Write();
}
// Write Mask
if (fIsMask && (drawOpt.find("MASK") != std::string::npos)) {
fMaskHist->Write();
}
// Write Covariances
if (drawOpt.find("COV") != std::string::npos && fFullCovar) {
PlotUtils::GetFullCovarPlot(fFullCovar, fSettings.GetName());
}
if (drawOpt.find("INVCOV") != std::string::npos && covar) {
PlotUtils::GetInvCovarPlot(covar, fSettings.GetName());
}
if (drawOpt.find("DECOMP") != std::string::npos && fDecomp) {
PlotUtils::GetDecompCovarPlot(fDecomp, fSettings.GetName());
}
// // Likelihood residual plots
// if (drawOpt.find("RESIDUAL") != std::string::npos) {
// WriteResidualPlots();
// }
// Ratio and Shape Plots
if (drawOpt.find("RATIO") != std::string::npos) {
WriteRatioPlot();
}
if (drawOpt.find("SHAPE") != std::string::npos) {
WriteShapePlot();
if (drawOpt.find("RATIO") != std::string::npos)
WriteShapeRatioPlot();
}
// // RATIO
// if (drawOpt.find("CANVMC") != std::string::npos) {
// TCanvas* c1 = WriteMCCanvas(fDataHist, fMCHist);
// c1->Write();
// delete c1;
// }
// // PDG
// if (drawOpt.find("CANVPDG") != std::string::npos && fMCHist_Modes) {
// TCanvas* c2 = WritePDGCanvas(fDataHist, fMCHist, fMCHist_Modes);
// c2->Write();
// delete c2;
// }
// Write Extra Histograms
AutoWriteExtraTH1();
WriteExtraHistograms();
if (fSaveSubMeas) {
for (std::vector::const_iterator expIter =
fSubChain.begin();
expIter != fSubChain.end(); expIter++) {
MeasurementBase* exp = *expIter;
exp->Write(drawOpt);
}
}
// Returning
LOG(SAM) << "Written Histograms: " << fName << std::endl;
return;
}
//********************************************************************
void JointMeas1D::WriteRatioPlot() {
//********************************************************************
// Setup mc data ratios
TH1D* dataRatio = (TH1D*)fDataHist->Clone((fName + "_data_RATIO").c_str());
TH1D* mcRatio = (TH1D*)fMCHist->Clone((fName + "_MC_RATIO").c_str());
// Extra MC Data Ratios
for (int i = 0; i < mcRatio->GetNbinsX(); i++) {
dataRatio->SetBinContent(i + 1, fDataHist->GetBinContent(i + 1) / fMCHist->GetBinContent(i + 1));
dataRatio->SetBinError(i + 1, fDataHist->GetBinError(i + 1) / fMCHist->GetBinContent(i + 1));
mcRatio->SetBinContent(i + 1, fMCHist->GetBinContent(i + 1) / fMCHist->GetBinContent(i + 1));
mcRatio->SetBinError(i + 1, fMCHist->GetBinError(i + 1) / fMCHist->GetBinContent(i + 1));
}
// Write ratios
mcRatio->Write();
dataRatio->Write();
delete mcRatio;
delete dataRatio;
}
//********************************************************************
void JointMeas1D::WriteShapePlot() {
//********************************************************************
TH1D* mcShape = (TH1D*)fMCHist->Clone((fName + "_MC_SHAPE").c_str());
double shapeScale = 1.0;
if (fIsRawEvents) {
shapeScale = fDataHist->Integral() / fMCHist->Integral();
} else {
shapeScale = fDataHist->Integral("width") / fMCHist->Integral("width");
}
mcShape->Scale(shapeScale);
std::stringstream ss;
ss << shapeScale;
mcShape->SetTitle(ss.str().c_str());
mcShape->SetLineWidth(3);
mcShape->SetLineStyle(7);
mcShape->Write();
delete mcShape;
}
//********************************************************************
void JointMeas1D::WriteShapeRatioPlot() {
//********************************************************************
// Get a mcshape histogram
TH1D* mcShape = (TH1D*)fMCHist->Clone((fName + "_MC_SHAPE").c_str());
double shapeScale = 1.0;
if (fIsRawEvents) {
shapeScale = fDataHist->Integral() / fMCHist->Integral();
} else {
shapeScale = fDataHist->Integral("width") / fMCHist->Integral("width");
}
mcShape->Scale(shapeScale);
// Create shape ratio histograms
TH1D* mcShapeRatio = (TH1D*)mcShape->Clone((fName + "_MC_SHAPE_RATIO").c_str());
TH1D* dataShapeRatio = (TH1D*)fDataHist->Clone((fName + "_data_SHAPE_RATIO").c_str());
// Divide the histograms
mcShapeRatio->Divide(mcShape);
dataShapeRatio->Divide(mcShape);
// Colour the shape ratio plots
mcShapeRatio->SetLineWidth(3);
mcShapeRatio->SetLineStyle(7);
mcShapeRatio->Write();
dataShapeRatio->Write();
delete mcShapeRatio;
delete dataShapeRatio;
}
// 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 "JointMeas1D.h"
/*
Constructor/Deconstuctor
*/
/*
Setup Functions
*/
//********************************************************************
void JointMeas1D::SetupMeasurement(std::string input, std::string type,
FitWeight* rw, std::string fkdt) {
//********************************************************************
// For joint samples, input files are given as a semi-colon seperated list.
// Parse this list and save it for later, and set up the types etc.
if (FitPar::Config().GetParB("EventManager")) {
ERR(FTL) << "Event Manager does not yet work with JointMeas1D Samples" << std::endl;
ERR(FTL) << "If you want good predictions for " << fName << " then run with it turned off! (-q EventManager=0)" << std::endl;
}
fSubInFiles.clear();
std::vector entries = GeneralUtils::ParseToStr(input, ";");
if (entries.size() < 2) {
ERR(FTL) << "Joint measurement expected to recieve at least two semi-colon "
"separated input files, but recieved: \""
<< input << "\"" << std::endl;
throw;
}
std::vector first_file_descriptor =
GeneralUtils::ParseToStr(entries.front(), ":");
if (first_file_descriptor.size() != 2) {
ERR(FTL) << "Found Joint measurement where the input file had no type: \""
<< input << "\", expected \"INPUTTYPE:File.root;File2.root\"."
<< std::endl;
throw;
}
std::string inpType = first_file_descriptor[0];
for (std::vector::iterator iter = entries.begin();
iter != entries.end(); iter++) {
if (GeneralUtils::ParseToStr(*iter, ":").size() != 2) {
std::stringstream ss("");
ss << inpType << ":" << (*iter);
fSubInFiles.push_back(ss.str());
} else {
fSubInFiles.push_back(*iter);
}
}
// Set Engine and Fake Data
fRW = rw;
fFakeDataInput = fkdt;
// Set Fit Options
SetFitOptions(type);
return;
}
/*
XSec Functions
*/
//********************************************************************
double JointMeas1D::TotalIntegratedFlux(std::string intOpt, double low,
double high) {
//********************************************************************
double totalflux = 0.0;
// Destroy the job for sub samples
for (std::vector::const_iterator expIter =
fSubChain.begin();
expIter != fSubChain.end(); expIter++) {
MeasurementBase* exp = *expIter;
double expflux = exp->TotalIntegratedFlux(intOpt, low, high);
// Fill flux options
if (fIsRatio) {
totalflux = expflux;
break;
}
if (fIsSummed) {
totalflux += expflux;
}
}
return totalflux;
}
/*
Reconfigure Functions
*/
//********************************************************************
void JointMeas1D::Reconfigure() {
//********************************************************************
for (std::vector::const_iterator expIter =
fSubChain.begin();
expIter != fSubChain.end(); expIter++) {
MeasurementBase* exp = *expIter;
exp->Reconfigure();
}
ConvertEventRates();
return;
}
//********************************************************************
void JointMeas1D::ConvertEventRates() {
//********************************************************************
// Apply Event Scaling
for (std::vector::const_iterator expIter =
fSubChain.begin();
expIter != fSubChain.end(); expIter++) {
MeasurementBase* exp = static_cast(*expIter);
exp->ScaleEvents();
}
// Joint function called by top level class
MakePlots();
// Do Final Normalisation
ApplyNormScale(fRW->GetSampleNorm(this->fName));
}
//********************************************************************
void JointMeas1D::MakePlots() {
//********************************************************************
// Reset the 1D histograms but not the subClasses
ResetAll();
// If Summed
if (fIsSummed) {
for (std::vector::const_iterator expIter =
fSubChain.begin();
expIter != fSubChain.end(); expIter++) {
MeasurementBase* exp = static_cast(*expIter);
this->fMCHist->Add(exp->GetMCList().at(0));
this->fMCFine->Add(exp->GetFineList().at(0));
}
return;
}
// If Ratio
if (fIsRatio) {
int sample = 0;
for (std::vector::const_iterator expIter =
fSubChain.begin();
expIter != fSubChain.end(); expIter++) {
MeasurementBase* exp = *expIter;
if (sample == 0) {
this->fMCHist->Add(exp->GetMCList().at(0));
this->fMCFine->Add(exp->GetFineList().at(0));
} else if (sample == 1) {
this->fMCHist->Divide(exp->GetMCList().at(0));
this->fMCFine->Divide(exp->GetFineList().at(0));
} else {
break;
}
sample++;
}
return;
}
return;
}
/*
Access Functions
*/
//********************************************************************
std::vector JointMeas1D::GetMCList() {
//********************************************************************
// Make Default Vector
std::vector tempVect;
tempVect.push_back(this->fMCHist);
// Return vector from all sub samples
for (std::vector::const_iterator expIter =
fSubChain.begin();
expIter != fSubChain.end(); expIter++) {
MeasurementBase* exp = *expIter;
std::vector subTempVect = exp->GetMCList();
for (UInt_t i = 0; i < subTempVect.size(); i++) {
tempVect.push_back(subTempVect.at(i));
}
}
return tempVect;
}
//********************************************************************
std::vector JointMeas1D::GetDataList() {
//********************************************************************
// Make Default Vector
std::vector tempVect;
tempVect.push_back(this->fDataHist);
// Return vector from all sub samples
for (std::vector::const_iterator expIter =
fSubChain.begin();
expIter != fSubChain.end(); expIter++) {
MeasurementBase* exp = *expIter;
std::vector subTempVect = exp->GetDataList();
for (UInt_t i = 0; i < subTempVect.size(); i++) {
tempVect.push_back(subTempVect.at(i));
}
}
return tempVect;
}
//********************************************************************
std::vector JointMeas1D::GetFineList() {
//********************************************************************
// Make Default Vector
std::vector tempVect;
tempVect.push_back(this->fMCFine);
// Return vector from all sub samples
for (std::vector::const_iterator expIter =
fSubChain.begin();
expIter != fSubChain.end(); expIter++) {
MeasurementBase* exp = *expIter;
std::vector subTempVect = exp->GetFineList();
for (UInt_t i = 0; i < subTempVect.size(); i++) {
tempVect.push_back(subTempVect.at(i));
}
}
return tempVect;
}
//********************************************************************
std::vector JointMeas1D::GetMaskList() {
//********************************************************************
// Make Default Vector
std::vector tempVect;
tempVect.push_back(this->fMaskHist);
// Return vector from all sub samples
for (std::vector::const_iterator expIter =
fSubChain.begin();
expIter != fSubChain.end(); expIter++) {
MeasurementBase* exp = *expIter;
std::vector subTempVect = exp->GetMaskList();
for (UInt_t i = 0; i < subTempVect.size(); i++) {
tempVect.push_back(subTempVect.at(i));
}
}
return tempVect;
}
//********************************************************************
std::vector JointMeas1D::GetFluxList() {
//********************************************************************
// Make Default Vector
std::vector tempVect;
tempVect.push_back(MeasurementBase::GetFluxHistogram());
// Return vector from all sub samples
for (std::vector::const_iterator expIter =
fSubChain.begin();
expIter != fSubChain.end(); expIter++) {
MeasurementBase* exp = *expIter;
std::vector subTempVect = exp->GetFluxList();
for (UInt_t i = 0; i < subTempVect.size(); i++) {
tempVect.push_back(subTempVect.at(i));
}
}
return tempVect;
}
//********************************************************************
std::vector JointMeas1D::GetEventRateList() {
//********************************************************************
// Make Default Vector
std::vector tempVect;
tempVect.push_back(MeasurementBase::GetEventHistogram());
// Return vector from all sub samples
for (std::vector::const_iterator expIter =
fSubChain.begin();
expIter != fSubChain.end(); expIter++) {
MeasurementBase* exp = *expIter;
std::vector subTempVect = exp->GetEventRateList();
for (UInt_t i = 0; i < subTempVect.size(); i++) {
tempVect.push_back(subTempVect.at(i));
}
}
return tempVect;
}
//********************************************************************
std::vector JointMeas1D::GetXSecList() {
//********************************************************************
// Make Default Vector
std::vector tempVect;
tempVect.push_back(MeasurementBase::GetXSecHistogram());
// Return vector from all sub samples
for (std::vector::const_iterator expIter =
fSubChain.begin();
expIter != fSubChain.end(); expIter++) {
MeasurementBase* exp = *expIter;
std::vector subTempVect = exp->GetXSecList();
for (UInt_t i = 0; i < subTempVect.size(); i++) {
tempVect.push_back(subTempVect.at(i));
}
}
return tempVect;
}
//********************************************************************
TH1D* JointMeas1D::GetCombinedFlux() {
//********************************************************************
TH1D* newflux = NULL;
int sample = 0;
for (std::vector::const_iterator expIter =
fSubChain.begin();
expIter != fSubChain.end(); expIter++) {
MeasurementBase* exp = *expIter;
// Get flux from experiment
std::vector fluxVect = exp->GetFluxList();
// Setup newflux
if (sample == 0) {
newflux = (TH1D*)fluxVect.at(0);
newflux->Reset();
}
// Add all fluxes
for (UInt_t i = 0; i < fluxVect.size(); i++) {
newflux->Add((TH1D*)fluxVect.at(i));
sample++;
}
}
if (!newflux){
ERR(FTL) << "No combined flux setup in JointMeas1D" << std::endl;
}
return newflux;
}
//********************************************************************
TH1D* JointMeas1D::GetCombinedEventRate() {
//********************************************************************
TH1D* newflux = NULL;
int sample = 0;
for (std::vector::const_iterator expIter =
fSubChain.begin();
expIter != fSubChain.end(); expIter++) {
MeasurementBase* exp = *expIter;
// Get flux from experiment
std::vector fluxVect = exp->GetFluxList();
// Setup newflux
if (sample == 0) {
newflux = (TH1D*)fluxVect.at(0);
newflux->Reset();
}
// Add all fluxes
for (UInt_t i = 0; i < fluxVect.size(); i++) {
newflux->Add(fluxVect.at(i));
sample++;
}
}
if (!newflux){
ERR(FTL) << "No combined event rate setup in JointMeas1D" << std::endl;
}
return newflux;
}
//********************************************************************
std::vector JointMeas1D::GetSubSamples() {
//********************************************************************
std::vector exps;
for (std::vector::const_iterator expIter =
fSubChain.begin();
expIter != fSubChain.end(); expIter++) {
exps.push_back(*expIter);
}
return exps;
}
//// CRAP TO BE REMOVED
//********************************************************************
void JointMeas1D::SetDataValues(std::string dataFile) {
//********************************************************************
// Override this function if the input file isn't in a suitable format
LOG(SAM) << "Reading data from: " << dataFile.c_str() << std::endl;
fDataHist =
PlotUtils::GetTH1DFromFile(dataFile, (fName + "_data"), fPlotTitles);
fDataTrue = (TH1D*)fDataHist->Clone();
// Number of data points is number of bins
fNDataPointsX = fDataHist->GetXaxis()->GetNbins();
return;
};
//********************************************************************
void JointMeas1D::SetDataFromDatabase(std::string inhistfile,
std::string histname) {
//********************************************************************
LOG(SAM) << "Filling histogram from " << inhistfile << "->" << histname
<< std::endl;
fDataHist = PlotUtils::GetTH1DFromRootFile(
(GeneralUtils::GetTopLevelDir() + "/data/" + inhistfile), histname);
fDataHist->SetNameTitle((fName + "_data").c_str(), (fName + "_data").c_str());
return;
};
//********************************************************************
void JointMeas1D::SetDataFromFile(std::string inhistfile,
std::string histname) {
//********************************************************************
LOG(SAM) << "Filling histogram from " << inhistfile << "->" << histname
<< std::endl;
fDataHist = PlotUtils::GetTH1DFromRootFile((inhistfile), histname);
fDataHist->SetNameTitle((fName + "_data").c_str(), (fName + "_data").c_str());
return;
};
//********************************************************************
void JointMeas1D::SetCovarMatrix(std::string covarFile) {
//********************************************************************
// Covariance function, only really used when reading in the MB Covariances.
TFile* tempFile = new TFile(covarFile.c_str(), "READ");
TH2D* covarPlot = new TH2D();
// TH2D* decmpPlot = new TH2D();
TH2D* covarInvPlot = new TH2D();
TH2D* fFullCovarPlot = new TH2D();
std::string covName = "";
std::string covOption = FitPar::Config().GetParS("thrown_covariance");
if (fIsShape || fIsFree) covName = "shp_";
if (fIsDiag)
covName += "diag";
else
covName += "full";
covarPlot = (TH2D*)tempFile->Get((covName + "cov").c_str());
covarInvPlot = (TH2D*)tempFile->Get((covName + "covinv").c_str());
if (!covOption.compare("SUB"))
fFullCovarPlot = (TH2D*)tempFile->Get((covName + "cov").c_str());
else if (!covOption.compare("FULL"))
fFullCovarPlot = (TH2D*)tempFile->Get("fullcov");
else
ERR(WRN) << "Incorrect thrown_covariance option in parameters."
<< std::endl;
int dim = int(fDataHist->GetNbinsX()); //-this->masked->Integral());
int covdim = int(fDataHist->GetNbinsX());
this->covar = new TMatrixDSym(dim);
fFullCovar = new TMatrixDSym(dim);
fDecomp = new TMatrixDSym(dim);
int row, column = 0;
row = 0;
column = 0;
for (Int_t i = 0; i < covdim; i++) {
// if (this->masked->GetBinContent(i+1) > 0) continue;
for (Int_t j = 0; j < covdim; j++) {
// if (this->masked->GetBinContent(j+1) > 0) continue;
(*this->covar)(row, column) = covarPlot->GetBinContent(i + 1, j + 1);
(*fFullCovar)(row, column) = fFullCovarPlot->GetBinContent(i + 1, j + 1);
column++;
}
column = 0;
row++;
}
// Set bin errors on data
if (!fIsDiag) {
StatUtils::SetDataErrorFromCov(fDataHist, fFullCovar);
}
// Get Deteriminant and inverse matrix
// fCovDet = this->covar->Determinant();
TDecompSVD LU = TDecompSVD(*this->covar);
this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), "");
return;
};
//********************************************************************
// Sets the covariance matrix from a provided file in a text format
// scale is a multiplicative pre-factor to apply in the case where the
// covariance is given in some unit (e.g. 1E-38)
void JointMeas1D::SetCovarMatrixFromText(std::string covarFile, int dim,
double scale) {
//********************************************************************
// Make a counter to track the line number
int row = 0;
std::string line;
- std::ifstream covarread(covarFile.c_str(), ifstream::in);
+ std::ifstream covarread(covarFile.c_str(), std::ifstream::in);
this->covar = new TMatrixDSym(dim);
fFullCovar = new TMatrixDSym(dim);
if (covarread.is_open())
LOG(SAM) << "Reading covariance matrix from file: " << covarFile
<< std::endl;
else
ERR(FTL) << "Covariance matrix provided is incorrect: " << covarFile
<< std::endl;
// Loop over the lines in the file
while (std::getline(covarread >> std::ws, line, '\n')) {
int column = 0;
// Loop over entries and insert them into matrix
std::vector entries = GeneralUtils::ParseToDbl(line, " ");
if (entries.size() <= 1) {
ERR(WRN) << "SetCovarMatrixFromText -> Covariance matrix only has <= 1 "
"entries on this line: "
<< row << std::endl;
}
for (std::vector::iterator iter = entries.begin();
iter != entries.end(); iter++) {
(*covar)(row, column) = *iter;
(*fFullCovar)(row, column) = *iter;
column++;
}
row++;
}
covarread.close();
// Scale the actualy covariance matrix by some multiplicative factor
(*fFullCovar) *= scale;
// Robust matrix inversion method
TDecompSVD LU = TDecompSVD(*this->covar);
// THIS IS ACTUALLY THE INVERSE COVARIANCE MATRIXA AAAAARGH
delete this->covar;
this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), "");
// Now need to multiply by the scaling factor
// If the covariance
(*this->covar) *= 1. / (scale);
return;
};
//********************************************************************
void JointMeas1D::SetCovarMatrixFromCorrText(std::string corrFile, int dim) {
//********************************************************************
// Make a counter to track the line number
int row = 0;
std::string line;
- std::ifstream corr(corrFile.c_str(), ifstream::in);
+ std::ifstream corr(corrFile.c_str(), std::ifstream::in);
this->covar = new TMatrixDSym(dim);
this->fFullCovar = new TMatrixDSym(dim);
if (corr.is_open())
LOG(SAM) << "Reading and converting correlation matrix from file: "
<< corrFile << std::endl;
else {
ERR(FTL) << "Correlation matrix provided is incorrect: " << corrFile
<< std::endl;
exit(-1);
}
while (std::getline(corr >> std::ws, line, '\n')) {
int column = 0;
// Loop over entries and insert them into matrix
// Multiply by the errors to get the covariance, rather than the correlation
// matrix
std::vector entries = GeneralUtils::ParseToDbl(line, " ");
for (std::vector::iterator iter = entries.begin();
iter != entries.end(); iter++) {
double val = (*iter) * this->fDataHist->GetBinError(row + 1) * 1E38 *
this->fDataHist->GetBinError(column + 1) * 1E38;
if (val == 0) {
ERR(FTL) << "Found a zero value in the covariance matrix, assuming "
"this is an error!"
<< std::endl;
exit(-1);
}
(*this->covar)(row, column) = val;
(*this->fFullCovar)(row, column) = val;
column++;
}
row++;
}
// Robust matrix inversion method
TDecompSVD LU = TDecompSVD(*this->covar);
delete this->covar;
this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), "");
return;
};
//********************************************************************
// FullUnits refers to if we have "real" unscaled units in the covariance matrix, e.g. 1E-76.
// If this is the case we need to scale it so that the chi2 contribution is correct
// NUISANCE internally assumes the covariance matrix has units of 1E76
void JointMeas1D::SetCovarFromDataFile(std::string covarFile,
std::string covName, bool FullUnits) {
//********************************************************************
LOG(SAM) << "Getting covariance from " << covarFile << "->" << covName
<< std::endl;
TFile* tempFile = new TFile(covarFile.c_str(), "READ");
TH2D* covPlot = (TH2D*)tempFile->Get(covName.c_str());
covPlot->SetDirectory(0);
// Scale the covariance matrix if it comes in normal units
if (FullUnits) {
covPlot->Scale(1.E76);
}
int dim = covPlot->GetNbinsX();
fFullCovar = new TMatrixDSym(dim);
for (int i = 0; i < dim; i++) {
for (int j = 0; j < dim; j++) {
(*fFullCovar)(i, j) = covPlot->GetBinContent(i + 1, j + 1);
}
}
this->covar = (TMatrixDSym*)fFullCovar->Clone();
fDecomp = (TMatrixDSym*)fFullCovar->Clone();
TDecompSVD LU = TDecompSVD(*this->covar);
this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), "");
TDecompChol LUChol = TDecompChol(*fDecomp);
LUChol.Decompose();
fDecomp = new TMatrixDSym(dim, LU.GetU().GetMatrixArray(), "");
return;
};
// std::vector JointMeas1D::GetMCList(void){
// std::vector temp;
// return temp;
// }
// std::vector JointMeas1D::GetDataList(void){
// std::vector temp;
// return temp;
// }
// std::vector JointMeas1D::GetMaskList(void){
// std::vector temp;
// return temp;
// }
// std::vector JointMeas1D::GetFineList(void){
// std::vector temp;
// return temp;
// }
diff --git a/src/FitBase/Measurement1D.cxx b/src/FitBase/Measurement1D.cxx
index cb2cdd3..5f88e1d 100644
--- a/src/FitBase/Measurement1D.cxx
+++ b/src/FitBase/Measurement1D.cxx
@@ -1,1903 +1,1903 @@
// Copyright 2016 L. Pickering, P. Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This ile 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 "Measurement1D.h"
//********************************************************************
Measurement1D::Measurement1D(void) {
//********************************************************************
// XSec Scalings
fScaleFactor = -1.0;
fCurrentNorm = 1.0;
// Histograms
fDataHist = NULL;
fDataTrue = NULL;
fMCHist = NULL;
fMCFine = NULL;
fMCWeighted = NULL;
fMaskHist = NULL;
// Covar
covar = NULL;
fFullCovar = NULL;
fShapeCovar = NULL;
fCovar = NULL;
fInvert = NULL;
fDecomp = NULL;
// Fake Data
fFakeDataInput = "";
fFakeDataFile = NULL;
// Options
fDefaultTypes = "FIX/FULL/CHI2";
fAllowedTypes =
"FIX,FREE,SHAPE/FULL,DIAG/CHI2/NORM/ENUCORR/Q2CORR/ENU1D/MASK/NOWIDTH";
fIsFix = false;
fIsShape = false;
fIsFree = false;
fIsDiag = false;
fIsFull = false;
fAddNormPen = false;
fIsMask = false;
fIsChi2SVD = false;
fIsRawEvents = false;
fIsNoWidth = false;
fIsDifXSec = false;
fIsEnu1D = false;
// Inputs
fInput = NULL;
fRW = NULL;
// Extra Histograms
fMCHist_Modes = NULL;
}
//********************************************************************
Measurement1D::~Measurement1D(void) {
//********************************************************************
if (fDataHist) delete fDataHist;
if (fDataTrue) delete fDataTrue;
if (fMCHist) delete fMCHist;
if (fMCFine) delete fMCFine;
if (fMCWeighted) delete fMCWeighted;
if (fMaskHist) delete fMaskHist;
if (covar) delete covar;
if (fFullCovar) delete fFullCovar;
if (fShapeCovar) delete fShapeCovar;
if (fCovar) delete fCovar;
if (fInvert) delete fInvert;
if (fDecomp) delete fDecomp;
}
//********************************************************************
void Measurement1D::FinaliseSampleSettings() {
//********************************************************************
MeasurementBase::FinaliseSampleSettings();
// Setup naming + renaming
fName = fSettings.GetName();
fSettings.SetS("originalname", fName);
if (fSettings.Has("rename")) {
fName = fSettings.GetS("rename");
fSettings.SetS("name", fName);
}
// Setup all other options
LOG(SAM) << "Finalising Sample Settings: " << fName << std::endl;
if ((fSettings.GetS("originalname").find("Evt") != std::string::npos)) {
fIsRawEvents = true;
LOG(SAM) << "Found event rate measurement but using poisson likelihoods."
<< std::endl;
}
if (fSettings.GetS("originalname").find("XSec_1DEnu") != std::string::npos) {
fIsEnu1D = true;
LOG(SAM) << "::" << fName << "::" << std::endl;
LOG(SAM) << "Found XSec Enu measurement, applying flux integrated scaling, "
<< "not flux averaged!" << std::endl;
}
if (fIsEnu1D && fIsRawEvents) {
LOG(SAM) << "Found 1D Enu XSec distribution AND fIsRawEvents, is this "
"really correct?!"
<< std::endl;
LOG(SAM) << "Check experiment constructor for " << fName
<< " and correct this!" << std::endl;
LOG(SAM) << "I live in " << __FILE__ << ":" << __LINE__ << std::endl;
exit(-1);
}
if (!fRW) fRW = FitBase::GetRW();
if (!fInput and !fIsJoint) SetupInputs(fSettings.GetS("input"));
// Setup options
SetFitOptions(fDefaultTypes); // defaults
SetFitOptions(fSettings.GetS("type")); // user specified
EnuMin = GeneralUtils::StrToDbl(fSettings.GetS("enu_min"));
EnuMax = GeneralUtils::StrToDbl(fSettings.GetS("enu_max"));
if (fAddNormPen) {
if (fNormError <= 0.0) {
ERR(WRN) << "Norm error for class " << fName << " is 0.0!" << std::endl;
ERR(WRN) << "If you want to use it please add fNormError=VAL" << std::endl;
throw;
}
}
}
//********************************************************************
void Measurement1D::CreateDataHistogram(int dimx, double* binx) {
//********************************************************************
if (fDataHist) delete fDataHist;
fDataHist = new TH1D( (fSettings.GetName() + "_data").c_str(), (fSettings.GetFullTitles()).c_str(),
dimx, binx) ;
}
//********************************************************************
void Measurement1D::SetDataFromTextFile(std::string datafile) {
//********************************************************************
LOG(SAM) << "Reading data from text file: " << datafile << std::endl;
fDataHist = PlotUtils::GetTH1DFromFile(datafile,
fSettings.GetName() + "_data",
fSettings.GetFullTitles());
}
//********************************************************************
void Measurement1D::SetDataFromRootFile(std::string datafile,
std::string histname) {
//********************************************************************
LOG(SAM) << "Reading data from root file: " << datafile << ";" << histname << std::endl;
fDataHist = PlotUtils::GetTH1DFromRootFile(datafile, histname);
fDataHist->SetNameTitle((fSettings.GetName() + "_data").c_str(),
(fSettings.GetFullTitles()).c_str());
return;
};
//********************************************************************
void Measurement1D::SetEmptyData(){
//********************************************************************
fDataHist = new TH1D("EMPTY_DATA","EMPTY_DATA",1,0.0,1.0);
}
//********************************************************************
void Measurement1D::SetPoissonErrors() {
//********************************************************************
if (!fDataHist) {
ERR(FTL) << "Need a data hist to setup possion errors! " << std::endl;
ERR(FTL) << "Setup Data First!" << std::endl;
throw;
}
for (int i = 0; i < fDataHist->GetNbinsX() + 1; i++) {
fDataHist->SetBinError(i + 1, sqrt(fDataHist->GetBinContent(i + 1)));
}
}
//********************************************************************
void Measurement1D::SetCovarFromDiagonal(TH1D* data) {
//********************************************************************
if (!data and fDataHist) {
data = fDataHist;
}
if (data) {
LOG(SAM) << "Setting diagonal covariance for: " << data->GetName() << std::endl;
fFullCovar = StatUtils::MakeDiagonalCovarMatrix(data);
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
} else {
ERR(FTL) << "No data input provided to set diagonal covar from!" << std::endl;
}
// if (!fIsDiag) {
// ERR(FTL) << "SetCovarMatrixFromDiag called for measurement "
// << "that is not set as diagonal." << std::endl;
// throw;
// }
}
//********************************************************************
void Measurement1D::SetCovarFromTextFile(std::string covfile, int dim) {
//********************************************************************
if (dim == -1) {
dim = fDataHist->GetNbinsX();
}
LOG(SAM) << "Reading covariance from text file: " << covfile << std::endl;
fFullCovar = StatUtils::GetCovarFromTextFile(covfile, dim);
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
//********************************************************************
void Measurement1D::SetCovarFromMultipleTextFiles(std::string covfiles, int dim) {
//********************************************************************
if (dim == -1) {
dim = fDataHist->GetNbinsX();
}
std::vector covList = GeneralUtils::ParseToStr(covfiles, ";");
fFullCovar = new TMatrixDSym(dim);
for (uint i = 0; i < covList.size(); ++i){
LOG(SAM) << "Reading covariance from text file: " << covList[i] << std::endl;
TMatrixDSym* temp_cov = StatUtils::GetCovarFromTextFile(covList[i], dim);
(*fFullCovar) += (*temp_cov);
delete temp_cov;
}
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
//********************************************************************
void Measurement1D::SetCovarFromRootFile(std::string covfile, std::string histname) {
//********************************************************************
LOG(SAM) << "Reading covariance from text file: " << covfile << ";" << histname << std::endl;
fFullCovar = StatUtils::GetCovarFromRootFile(covfile, histname);
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
//********************************************************************
void Measurement1D::SetCovarInvertFromTextFile(std::string covfile, int dim) {
//********************************************************************
if (dim == -1) {
dim = fDataHist->GetNbinsX();
}
LOG(SAM) << "Reading inverted covariance from text file: " << covfile << std::endl;
covar = StatUtils::GetCovarFromTextFile(covfile, dim);
fFullCovar = StatUtils::GetInvert(covar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
//********************************************************************
void Measurement1D::SetCovarInvertFromRootFile(std::string covfile, std::string histname) {
//********************************************************************
LOG(SAM) << "Reading inverted covariance from text file: " << covfile << ";" << histname << std::endl;
covar = StatUtils::GetCovarFromRootFile(covfile, histname);
fFullCovar = StatUtils::GetInvert(covar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
//********************************************************************
void Measurement1D::SetCorrelationFromTextFile(std::string covfile, int dim) {
//********************************************************************
if (dim == -1) dim = fDataHist->GetNbinsX();
LOG(SAM) << "Reading data correlations from text file: " << covfile << ";" << dim << std::endl;
TMatrixDSym* correlation = StatUtils::GetCovarFromTextFile(covfile, dim);
if (!fDataHist) {
ERR(FTL) << "Trying to set correlations from text file but there is no data to build it from. \n"
<< "In constructor make sure data is set before SetCorrelationFromTextFile is called. \n" << std::endl;
throw;
}
// Fill covar from data errors and correlations
fFullCovar = new TMatrixDSym(dim);
for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
for (int j = 0; j < fDataHist->GetNbinsX(); j++) {
(*fFullCovar)(i, j) = (*correlation)(i, j) * fDataHist->GetBinError(i + 1) * fDataHist->GetBinError(j + 1) * 1.E76;
}
}
// Fill other covars.
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
delete correlation;
}
//********************************************************************
void Measurement1D::SetCorrelationFromMultipleTextFiles(std::string corrfiles, int dim) {
//********************************************************************
if (dim == -1) {
dim = fDataHist->GetNbinsX();
}
std::vector corrList = GeneralUtils::ParseToStr(corrfiles, ";");
fFullCovar = new TMatrixDSym(dim);
for (uint i = 0; i < corrList.size(); ++i){
LOG(SAM) << "Reading covariance from text file: " << corrList[i] << std::endl;
TMatrixDSym* temp_cov = StatUtils::GetCovarFromTextFile(corrList[i], dim);
for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
for (int j = 0; j < fDataHist->GetNbinsX(); j++) {
(*temp_cov)(i, j) = (*temp_cov)(i, j) * fDataHist->GetBinError(i + 1) * fDataHist->GetBinError(j + 1) * 1.E76;
}
}
(*fFullCovar) += (*temp_cov);
delete temp_cov;
}
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
//********************************************************************
void Measurement1D::SetCorrelationFromRootFile(std::string covfile, std::string histname) {
//********************************************************************
LOG(SAM) << "Reading data correlations from text file: " << covfile << ";" << histname << std::endl;
TMatrixDSym* correlation = StatUtils::GetCovarFromRootFile(covfile, histname);
if (!fDataHist) {
ERR(FTL) << "Trying to set correlations from text file but there is no data to build it from. \n"
<< "In constructor make sure data is set before SetCorrelationFromTextFile is called. \n" << std::endl;
throw;
}
// Fill covar from data errors and correlations
fFullCovar = new TMatrixDSym(fDataHist->GetNbinsX());
for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
for (int j = 0; j < fDataHist->GetNbinsX(); j++) {
(*fFullCovar)(i, j) = (*correlation)(i, j) * fDataHist->GetBinError(i + 1) * fDataHist->GetBinError(j + 1) * 1.E76;
}
}
// Fill other covars.
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
delete correlation;
}
//********************************************************************
void Measurement1D::SetCholDecompFromTextFile(std::string covfile, int dim) {
//********************************************************************
if (dim == -1) {
dim = fDataHist->GetNbinsX();
}
LOG(SAM) << "Reading cholesky from text file: " << covfile << std::endl;
TMatrixD* temp = StatUtils::GetMatrixFromTextFile(covfile, dim, dim);
TMatrixD* trans = (TMatrixD*)temp->Clone();
trans->T();
(*trans) *= (*temp);
fFullCovar = new TMatrixDSym(dim, trans->GetMatrixArray(), "");
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
delete temp;
delete trans;
}
//********************************************************************
void Measurement1D::SetCholDecompFromRootFile(std::string covfile, std::string histname) {
//********************************************************************
LOG(SAM) << "Reading cholesky decomp from root file: " << covfile << ";" << histname << std::endl;
TMatrixD* temp = StatUtils::GetMatrixFromRootFile(covfile, histname);
TMatrixD* trans = (TMatrixD*)temp->Clone();
trans->T();
(*trans) *= (*temp);
fFullCovar = new TMatrixDSym(temp->GetNrows(), trans->GetMatrixArray(), "");
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
delete temp;
delete trans;
}
void Measurement1D::SetShapeCovar(){
// Return if this is missing any pre-requisites
if (!fFullCovar) return;
if (!fDataHist) return;
// Also return if it's bloody stupid under the circumstances
if (fIsDiag) return;
fShapeCovar = StatUtils::ExtractShapeOnlyCovar(fFullCovar, fDataHist);
return;
}
//********************************************************************
void Measurement1D::ScaleData(double scale) {
//********************************************************************
fDataHist->Scale(scale);
}
//********************************************************************
void Measurement1D::ScaleDataErrors(double scale) {
//********************************************************************
for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
fDataHist->SetBinError(i + 1, fDataHist->GetBinError(i + 1) * scale);
}
}
//********************************************************************
void Measurement1D::ScaleCovar(double scale) {
//********************************************************************
(*fFullCovar) *= scale;
(*covar) *= 1.0 / scale;
(*fDecomp) *= sqrt(scale);
}
//********************************************************************
void Measurement1D::SetBinMask(std::string maskfile) {
//********************************************************************
if (!fIsMask) return;
LOG(SAM) << "Reading bin mask from file: " << maskfile << std::endl;
// Create a mask histogram with dim of data
int nbins = fDataHist->GetNbinsX();
fMaskHist =
new TH1I((fSettings.GetName() + "_BINMASK").c_str(),
(fSettings.GetName() + "_BINMASK; Bin; Mask?").c_str(), nbins, 0, nbins);
std::string line;
- std::ifstream mask(maskfile.c_str(), ifstream::in);
+ std::ifstream mask(maskfile.c_str(), std::ifstream::in);
if (!mask.is_open()) {
LOG(FTL) << " Cannot find mask file." << std::endl;
throw;
}
while (std::getline(mask >> std::ws, line, '\n')) {
std::vector entries = GeneralUtils::ParseToInt(line, " ");
// Skip lines with poorly formatted lines
if (entries.size() < 2) {
LOG(WRN) << "Measurement1D::SetBinMask(), couldn't parse line: " << line
<< std::endl;
continue;
}
// The first index should be the bin number, the second should be the mask
// value.
int val = 0;
if (entries[1] > 0) val = 1;
fMaskHist->SetBinContent(entries[0], val);
}
// Apply masking by setting masked data bins to zero
PlotUtils::MaskBins(fDataHist, fMaskHist);
return;
}
//********************************************************************
void Measurement1D::FinaliseMeasurement() {
//********************************************************************
LOG(SAM) << "Finalising Measurement: " << fName << std::endl;
if (fSettings.GetB("onlymc")){
if (fDataHist) delete fDataHist;
fDataHist = new TH1D("empty_data","empty_data",1,0.0,1.0);
}
// Make sure data is setup
if (!fDataHist) {
ERR(FTL) << "No data has been setup inside " << fName << " constructor!" << std::endl;
throw;
}
// Make sure covariances are setup
if (!fFullCovar) {
fIsDiag = true;
SetCovarFromDiagonal(fDataHist);
}
if (!covar) {
covar = StatUtils::GetInvert(fFullCovar);
}
if (!fDecomp) {
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
// Push the diagonals of fFullCovar onto the data histogram
// Comment this out until the covariance/data scaling is consistent!
StatUtils::SetDataErrorFromCov(fDataHist, fFullCovar, 1E-38);
// If shape only, set covar and fDecomp using the shape-only matrix (if set)
if (fIsShape && fShapeCovar and FitPar::Config().GetParB("UseShapeCovar")){
if (covar) delete covar;
covar = StatUtils::GetInvert(fShapeCovar);
if (fDecomp) delete fDecomp;
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
// Setup fMCHist from data
fMCHist = (TH1D*)fDataHist->Clone();
fMCHist->SetNameTitle((fSettings.GetName() + "_MC").c_str(),
(fSettings.GetFullTitles()).c_str());
fMCHist->Reset();
// Setup fMCFine
fMCFine = new TH1D("mcfine", "mcfine", fDataHist->GetNbinsX() * 8,
fMCHist->GetBinLowEdge(1),
fMCHist->GetBinLowEdge(fDataHist->GetNbinsX() + 1));
fMCFine->SetNameTitle((fSettings.GetName() + "_MC_FINE").c_str(),
(fSettings.GetFullTitles()).c_str());
fMCFine->Reset();
// Setup MC Stat
fMCStat = (TH1D*)fMCHist->Clone();
fMCStat->Reset();
// Search drawopts for possible types to include by default
std::string drawopts = FitPar::Config().GetParS("drawopts");
if (drawopts.find("MODES") != std::string::npos) {
fMCHist_Modes = new TrueModeStack( (fSettings.GetName() + "_MODES").c_str(),
("True Channels"), fMCHist);
SetAutoProcessTH1(fMCHist_Modes, kCMD_Reset, kCMD_Norm, kCMD_Write);
}
// Setup bin masks using sample name
if (fIsMask) {
std::string curname = fName;
std::string origname = fSettings.GetS("originalname");
// Check rename.mask
std::string maskloc = FitPar::Config().GetParDIR(curname + ".mask");
// Check origname.mask
if (maskloc.empty()) maskloc = FitPar::Config().GetParDIR(origname + ".mask");
// Check database
if (maskloc.empty()) {
maskloc = FitPar::GetDataBase() + "/masks/" + origname + ".mask";
}
// Setup Bin Mask
SetBinMask(maskloc);
}
if (fScaleFactor < 0) {
ERR(FTL) << "I found a negative fScaleFactor in " << __FILE__ << ":" << __LINE__ << std::endl;
ERR(FTL) << "fScaleFactor = " << fScaleFactor << std::endl;
ERR(FTL) << "EXITING" << std::endl;
throw;
}
// Create and fill Weighted Histogram
if (!fMCWeighted) {
fMCWeighted = (TH1D*)fMCHist->Clone();
fMCWeighted->SetNameTitle((fName + "_MCWGHTS").c_str(),
(fName + "_MCWGHTS" + fPlotTitles).c_str());
fMCWeighted->GetYaxis()->SetTitle("Weighted Events");
}
}
//********************************************************************
void Measurement1D::SetFitOptions(std::string opt) {
//********************************************************************
// Do nothing if default given
if (opt == "DEFAULT") return;
// CHECK Conflicting Fit Options
std::vector fit_option_allow =
GeneralUtils::ParseToStr(fAllowedTypes, "/");
for (UInt_t i = 0; i < fit_option_allow.size(); i++) {
std::vector fit_option_section =
GeneralUtils::ParseToStr(fit_option_allow.at(i), ",");
bool found_option = false;
for (UInt_t j = 0; j < fit_option_section.size(); j++) {
std::string av_opt = fit_option_section.at(j);
if (!found_option and opt.find(av_opt) != std::string::npos) {
found_option = true;
} else if (found_option and opt.find(av_opt) != std::string::npos) {
ERR(FTL) << "ERROR: Conflicting fit options provided: "
<< opt << std::endl
<< "Conflicting group = " << fit_option_section.at(i) << std::endl
<< "You should only supply one of these options in card file." << std::endl;
throw;
}
}
}
// Check all options are allowed
std::vector fit_options_input =
GeneralUtils::ParseToStr(opt, "/");
for (UInt_t i = 0; i < fit_options_input.size(); i++) {
if (fAllowedTypes.find(fit_options_input.at(i)) == std::string::npos) {
ERR(FTL) << "ERROR: Fit Option '" << fit_options_input.at(i)
<< "' Provided is not allowed for this measurement."
<< std::endl;
ERR(FTL) << "Fit Options should be provided as a '/' seperated list "
"(e.g. FREE/DIAG/NORM)"
<< std::endl;
ERR(FTL) << "Available options for " << fName << " are '" << fAllowedTypes
<< "'" << std::endl;
throw;
}
}
// Set TYPE
fFitType = opt;
// FIX,SHAPE,FREE
if (opt.find("FIX") != std::string::npos) {
fIsFree = fIsShape = false;
fIsFix = true;
} else if (opt.find("SHAPE") != std::string::npos) {
fIsFree = fIsFix = false;
fIsShape = true;
} else if (opt.find("FREE") != std::string::npos) {
fIsFix = fIsShape = false;
fIsFree = true;
}
// DIAG,FULL (or default to full)
if (opt.find("DIAG") != std::string::npos) {
fIsDiag = true;
fIsFull = false;
} else if (opt.find("FULL") != std::string::npos) {
fIsDiag = false;
fIsFull = true;
}
// CHI2/LL (OTHERS?)
if (opt.find("LOG") != std::string::npos) {
fIsChi2 = false;
ERR(FTL) << "No other LIKELIHOODS properly supported!" << std::endl;
ERR(FTL) << "Try to use a chi2!" << std::endl;
throw;
} else {
fIsChi2 = true;
}
// EXTRAS
if (opt.find("RAW") != std::string::npos) fIsRawEvents = true;
if (opt.find("NOWIDTH") != std::string::npos) fIsNoWidth = true;
if (opt.find("DIF") != std::string::npos) fIsDifXSec = true;
if (opt.find("ENU1D") != std::string::npos) fIsEnu1D = true;
if (opt.find("NORM") != std::string::npos) fAddNormPen = true;
if (opt.find("MASK") != std::string::npos) fIsMask = true;
return;
};
//********************************************************************
void Measurement1D::SetSmearingMatrix(std::string smearfile, int truedim,
int recodim) {
//********************************************************************
// The smearing matrix describes the migration from true bins (rows) to reco
// bins (columns)
// Counter over the true bins!
int row = 0;
std::string line;
- std::ifstream smear(smearfile.c_str(), ifstream::in);
+ std::ifstream smear(smearfile.c_str(), std::ifstream::in);
// Note that the smearing matrix may be rectangular.
fSmearMatrix = new TMatrixD(truedim, recodim);
if (smear.is_open())
LOG(SAM) << "Reading smearing matrix from file: " << smearfile << std::endl;
else
ERR(FTL) << "Smearing matrix provided is incorrect: " << smearfile
<< std::endl;
while (std::getline(smear >> std::ws, line, '\n')) {
int column = 0;
std::vector entries = GeneralUtils::ParseToDbl(line, " ");
for (std::vector::iterator iter = entries.begin();
iter != entries.end(); iter++) {
(*fSmearMatrix)(row, column) =
(*iter) / 100.; // Convert to fraction from
// percentage (this may not be
// general enough)
column++;
}
row++;
}
return;
}
//********************************************************************
void Measurement1D::ApplySmearingMatrix() {
//********************************************************************
if (!fSmearMatrix) {
ERR(WRN) << fName
<< ": attempted to apply smearing matrix, but none was set"
<< std::endl;
return;
}
TH1D* unsmeared = (TH1D*)fMCHist->Clone();
TH1D* smeared = (TH1D*)fMCHist->Clone();
smeared->Reset();
// Loop over reconstructed bins
// true = row; reco = column
for (int rbin = 0; rbin < fSmearMatrix->GetNcols(); ++rbin) {
// Sum up the constributions from all true bins
double rBinVal = 0;
// Loop over true bins
for (int tbin = 0; tbin < fSmearMatrix->GetNrows(); ++tbin) {
rBinVal +=
(*fSmearMatrix)(tbin, rbin) * unsmeared->GetBinContent(tbin + 1);
}
smeared->SetBinContent(rbin + 1, rBinVal);
}
fMCHist = (TH1D*)smeared->Clone();
return;
}
/*
Reconfigure LOOP
*/
//********************************************************************
void Measurement1D::ResetAll() {
//********************************************************************
fMCHist->Reset();
fMCFine->Reset();
fMCStat->Reset();
return;
};
//********************************************************************
void Measurement1D::FillHistograms() {
//********************************************************************
if (Signal) {
fMCHist->Fill(fXVar, Weight);
fMCFine->Fill(fXVar, Weight);
fMCStat->Fill(fXVar, 1.0);
if (fMCHist_Modes) fMCHist_Modes->Fill(Mode, fXVar, Weight);
}
return;
};
//********************************************************************
void Measurement1D::ScaleEvents() {
//********************************************************************
// Fill MCWeighted;
// for (int i = 0; i < fMCHist->GetNbinsX(); i++) {
// fMCWeighted->SetBinContent(i + 1, fMCHist->GetBinContent(i + 1));
// fMCWeighted->SetBinError(i + 1, fMCHist->GetBinError(i + 1));
// }
// Setup Stat ratios for MC and MC Fine
double* statratio = new double[fMCHist->GetNbinsX()];
for (int i = 0; i < fMCHist->GetNbinsX(); i++) {
if (fMCHist->GetBinContent(i + 1) != 0) {
statratio[i] = fMCHist->GetBinError(i + 1) / fMCHist->GetBinContent(i + 1);
} else {
statratio[i] = 0.0;
}
}
double* statratiofine = new double[fMCFine->GetNbinsX()];
for (int i = 0; i < fMCFine->GetNbinsX(); i++) {
if (fMCFine->GetBinContent(i + 1) != 0) {
statratiofine[i] = fMCFine->GetBinError(i + 1) / fMCFine->GetBinContent(i + 1);
} else {
statratiofine[i] = 0.0;
}
}
// Scaling for raw event rates
if (fIsRawEvents) {
double datamcratio = fDataHist->Integral() / fMCHist->Integral();
fMCHist->Scale(datamcratio);
fMCFine->Scale(datamcratio);
if (fMCHist_Modes) fMCHist_Modes->Scale(datamcratio);
// Scaling for XSec as function of Enu
} else if (fIsEnu1D) {
PlotUtils::FluxUnfoldedScaling(fMCHist, GetFluxHistogram(),
GetEventHistogram(), fScaleFactor,
fNEvents);
PlotUtils::FluxUnfoldedScaling(fMCFine, GetFluxHistogram(),
GetEventHistogram(), fScaleFactor,
fNEvents);
// if (fMCHist_Modes) {
// PlotUtils::FluxUnfoldedScaling(fMCHist_Modes, GetFluxHistogram(),
// GetEventHistogram(), fScaleFactor,
// fNEvents);
// }
} else if (fIsNoWidth) {
fMCHist->Scale(fScaleFactor);
fMCFine->Scale(fScaleFactor);
if (fMCHist_Modes) fMCHist_Modes->Scale(fScaleFactor);
// Any other differential scaling
} else {
fMCHist->Scale(fScaleFactor, "width");
fMCFine->Scale(fScaleFactor, "width");
if (fMCHist_Modes) fMCHist_Modes->Scale(fScaleFactor, "width");
}
// Proper error scaling - ROOT Freaks out with xsec weights sometimes
for (int i = 0; i < fMCStat->GetNbinsX(); i++) {
fMCHist->SetBinError(i + 1, fMCHist->GetBinContent(i + 1) * statratio[i]);
}
for (int i = 0; i < fMCFine->GetNbinsX(); i++) {
fMCFine->SetBinError(i + 1, fMCFine->GetBinContent(i + 1) * statratiofine[i]);
}
// Clean up
delete statratio;
delete statratiofine;
return;
};
//********************************************************************
void Measurement1D::ApplyNormScale(double norm) {
//********************************************************************
fCurrentNorm = norm;
fMCHist->Scale(1.0 / norm);
fMCFine->Scale(1.0 / norm);
return;
};
/*
Statistic Functions - Outsources to StatUtils
*/
//********************************************************************
int Measurement1D::GetNDOF() {
//********************************************************************
int ndof = fDataHist->GetNbinsX();
if (fMaskHist and fIsMask) ndof -= fMaskHist->Integral();
return ndof;
}
//********************************************************************
double Measurement1D::GetLikelihood() {
//********************************************************************
// If this is for a ratio, there is no data histogram to compare to!
if (fNoData || !fDataHist) return 0.;
// Apply Masking to MC if Required.
if (fIsMask and fMaskHist) {
PlotUtils::MaskBins(fMCHist, fMaskHist);
}
// Sort Shape Scaling
double scaleF = 0.0;
// TODO Include !fIsRawEvents
if (fIsShape) {
if (fMCHist->Integral(1, fMCHist->GetNbinsX(), "width")) {
scaleF = fDataHist->Integral(1, fDataHist->GetNbinsX(), "width") /
fMCHist->Integral(1, fMCHist->GetNbinsX(), "width");
fMCHist->Scale(scaleF);
fMCFine->Scale(scaleF);
}
}
// Likelihood Calculation
double stat = 0.;
if (fIsChi2) {
if (fIsRawEvents) {
stat = StatUtils::GetChi2FromEventRate(fDataHist, fMCHist, fMaskHist);
} else if (fIsDiag) {
stat = StatUtils::GetChi2FromDiag(fDataHist, fMCHist, fMaskHist);
} else if (!fIsDiag and !fIsRawEvents) {
stat = StatUtils::GetChi2FromCov(fDataHist, fMCHist, covar, fMaskHist);
}
}
// Sort Penalty Terms
if (fAddNormPen) {
double penalty =
(1. - fCurrentNorm) * (1. - fCurrentNorm) / (fNormError * fNormError);
stat += penalty;
}
// Return to normal scaling
if (fIsShape) { // and !FitPar::Config().GetParB("saveshapescaling")) {
fMCHist->Scale(1. / scaleF);
fMCFine->Scale(1. / scaleF);
}
fLikelihood = stat;
return stat;
}
/*
Fake Data Functions
*/
//********************************************************************
void Measurement1D::SetFakeDataValues(std::string fakeOption) {
//********************************************************************
// Setup original/datatrue
TH1D* tempdata = (TH1D*) fDataHist->Clone();
if (!fIsFakeData) {
fIsFakeData = true;
// Make a copy of the original data histogram.
if (!fDataOrig) fDataOrig = (TH1D*)fDataHist->Clone((fName + "_data_original").c_str());
} else {
ResetFakeData();
}
// Setup Inputs
fFakeDataInput = fakeOption;
LOG(SAM) << "Setting fake data from : " << fFakeDataInput << std::endl;
// From MC
if (fFakeDataInput.compare("MC") == 0) {
fDataHist = (TH1D*)fMCHist->Clone((fName + "_MC").c_str());
// Fake File
} else {
if (!fFakeDataFile) fFakeDataFile = new TFile(fFakeDataInput.c_str(), "READ");
fDataHist = (TH1D*)fFakeDataFile->Get((fName + "_MC").c_str());
}
// Setup Data Hist
fDataHist->SetNameTitle((fName + "_FAKE").c_str(),
(fName + fPlotTitles).c_str());
// Replace Data True
if (fDataTrue) delete fDataTrue;
fDataTrue = (TH1D*)fDataHist->Clone();
fDataTrue->SetNameTitle((fName + "_FAKE_TRUE").c_str(),
(fName + fPlotTitles).c_str());
// Make a new covariance for fake data hist.
int nbins = fDataHist->GetNbinsX();
double alpha_i = 0.0;
double alpha_j = 0.0;
for (int i = 0; i < nbins; i++) {
for (int j = 0; j < nbins; j++) {
alpha_i = fDataHist->GetBinContent(i + 1) / tempdata->GetBinContent(i + 1);
alpha_j = fDataHist->GetBinContent(j + 1) / tempdata->GetBinContent(j + 1);
(*fFullCovar)(i, j) = alpha_i * alpha_j * (*fFullCovar)(i, j);
}
}
// Setup Covariances
if (covar) delete covar;
covar = StatUtils::GetInvert(fFullCovar);
if (fDecomp) delete fDecomp;
fDecomp = StatUtils::GetInvert(fFullCovar);
delete tempdata;
return;
};
//********************************************************************
void Measurement1D::ResetFakeData() {
//********************************************************************
if (fIsFakeData) {
if (fDataHist) delete fDataHist;
fDataHist = (TH1D*)fDataTrue->Clone((fSettings.GetName() + "_FKDAT").c_str());
}
}
//********************************************************************
void Measurement1D::ResetData() {
//********************************************************************
if (fIsFakeData) {
if (fDataHist) delete fDataHist;
fDataHist = (TH1D*)fDataOrig->Clone((fSettings.GetName() + "_data").c_str());
}
fIsFakeData = false;
}
//********************************************************************
void Measurement1D::ThrowCovariance() {
//********************************************************************
// Take a fDecomposition and use it to throw the current dataset.
// Requires fDataTrue also be set incase used repeatedly.
if (!fDataTrue) fDataTrue = (TH1D*) fDataHist->Clone();
if (fDataHist) delete fDataHist;
fDataHist = StatUtils::ThrowHistogram(fDataTrue, fFullCovar);
return;
};
//********************************************************************
void Measurement1D::ThrowDataToy(){
//********************************************************************
if (!fDataTrue) fDataTrue = (TH1D*) fDataHist->Clone();
if (fMCHist) delete fMCHist;
fMCHist = StatUtils::ThrowHistogram(fDataTrue, fFullCovar);
}
/*
Access Functions
*/
//********************************************************************
TH1D* Measurement1D::GetMCHistogram() {
//********************************************************************
if (!fMCHist) return fMCHist;
std::ostringstream chi2;
chi2 << std::setprecision(5) << this->GetLikelihood();
int linecolor = kRed;
int linestyle = 1;
int linewidth = 1;
int fillcolor = 0;
int fillstyle = 1001;
// if (fSettings.Has("linecolor")) linecolor = fSettings.GetI("linecolor");
// if (fSettings.Has("linestyle")) linestyle = fSettings.GetI("linestyle");
// if (fSettings.Has("linewidth")) linewidth = fSettings.GetI("linewidth");
// if (fSettings.Has("fillcolor")) fillcolor = fSettings.GetI("fillcolor");
// if (fSettings.Has("fillstyle")) fillstyle = fSettings.GetI("fillstyle");
fMCHist->SetTitle(chi2.str().c_str());
fMCHist->SetLineColor(linecolor);
fMCHist->SetLineStyle(linestyle);
fMCHist->SetLineWidth(linewidth);
fMCHist->SetFillColor(fillcolor);
fMCHist->SetFillStyle(fillstyle);
return fMCHist;
};
//********************************************************************
TH1D* Measurement1D::GetDataHistogram() {
//********************************************************************
if (!fDataHist) return fDataHist;
int datacolor = kBlack;
int datastyle = 1;
int datawidth = 1;
// if (fSettings.Has("datacolor")) datacolor = fSettings.GetI("datacolor");
// if (fSettings.Has("datastyle")) datastyle = fSettings.GetI("datastyle");
// if (fSettings.Has("datawidth")) datawidth = fSettings.GetI("datawidth");
fDataHist->SetLineColor(datacolor);
fDataHist->SetLineWidth(datawidth);
fDataHist->SetMarkerStyle(datastyle);
return fDataHist;
};
/*
Write Functions
*/
// Save all the histograms at once
//********************************************************************
void Measurement1D::Write(std::string drawOpt) {
//********************************************************************
// Get Draw Options
drawOpt = FitPar::Config().GetParS("drawopts");
// Write Settigns
if (drawOpt.find("SETTINGS") != std::string::npos){
fSettings.Set("#chi^{2}",fLikelihood);
fSettings.Set("NDOF", this->GetNDOF() );
fSettings.Set("#chi^{2}/NDOF", fLikelihood / this->GetNDOF() );
fSettings.Write();
}
// Write Data/MC
GetDataList().at(0)->Write();
GetMCList().at(0)->Write();
if((fEvtRateScaleFactor != 0xdeadbeef) && GetMCList().at(0)){
TH1D * PredictedEvtRate = static_cast(GetMCList().at(0)->Clone());
PredictedEvtRate->Scale(fEvtRateScaleFactor);
PredictedEvtRate->GetYaxis()->SetTitle("Predicted event rate");
PredictedEvtRate->Write();
}
// Write Fine Histogram
if (drawOpt.find("FINE") != std::string::npos)
GetFineList().at(0)->Write();
// Write Weighted Histogram
if (drawOpt.find("WEIGHTS") != std::string::npos && fMCWeighted)
fMCWeighted->Write();
// Save Flux/Evt if no event manager
if (!FitPar::Config().GetParB("EventManager")) {
if (drawOpt.find("FLUX") != std::string::npos && GetFluxHistogram())
GetFluxHistogram()->Write();
if (drawOpt.find("EVT") != std::string::npos && GetEventHistogram())
GetEventHistogram()->Write();
if (drawOpt.find("XSEC") != std::string::npos && GetEventHistogram())
GetXSecHistogram()->Write();
}
// Write Mask
if (fIsMask && (drawOpt.find("MASK") != std::string::npos)) {
fMaskHist->Write();
}
// Write Covariances
if (drawOpt.find("COV") != std::string::npos && fFullCovar) {
PlotUtils::GetFullCovarPlot(fFullCovar, fSettings.GetName());
}
if (drawOpt.find("INVCOV") != std::string::npos && covar) {
PlotUtils::GetInvCovarPlot(covar, fSettings.GetName());
}
if (drawOpt.find("DECOMP") != std::string::npos && fDecomp) {
PlotUtils::GetDecompCovarPlot(fDecomp, fSettings.GetName());
}
// // Likelihood residual plots
// if (drawOpt.find("RESIDUAL") != std::string::npos) {
// WriteResidualPlots();
// }
// Ratio and Shape Plots
if (drawOpt.find("RATIO") != std::string::npos) {
WriteRatioPlot();
}
if (drawOpt.find("SHAPE") != std::string::npos) {
WriteShapePlot();
if (drawOpt.find("RATIO") != std::string::npos)
WriteShapeRatioPlot();
}
// // RATIO
// if (drawOpt.find("CANVMC") != std::string::npos) {
// TCanvas* c1 = WriteMCCanvas(fDataHist, fMCHist);
// c1->Write();
// delete c1;
// }
// // PDG
// if (drawOpt.find("CANVPDG") != std::string::npos && fMCHist_Modes) {
// TCanvas* c2 = WritePDGCanvas(fDataHist, fMCHist, fMCHist_Modes);
// c2->Write();
// delete c2;
// }
// Write Extra Histograms
AutoWriteExtraTH1();
WriteExtraHistograms();
// Returning
LOG(SAM) << "Written Histograms: " << fName << std::endl;
return;
}
//********************************************************************
void Measurement1D::WriteRatioPlot() {
//********************************************************************
// Setup mc data ratios
TH1D* dataRatio = (TH1D*)fDataHist->Clone((fName + "_data_RATIO").c_str());
TH1D* mcRatio = (TH1D*)fMCHist->Clone((fName + "_MC_RATIO").c_str());
// Extra MC Data Ratios
for (int i = 0; i < mcRatio->GetNbinsX(); i++) {
dataRatio->SetBinContent(i + 1, fDataHist->GetBinContent(i + 1) / fMCHist->GetBinContent(i + 1));
dataRatio->SetBinError(i + 1, fDataHist->GetBinError(i + 1) / fMCHist->GetBinContent(i + 1));
mcRatio->SetBinContent(i + 1, fMCHist->GetBinContent(i + 1) / fMCHist->GetBinContent(i + 1));
mcRatio->SetBinError(i + 1, fMCHist->GetBinError(i + 1) / fMCHist->GetBinContent(i + 1));
}
// Write ratios
mcRatio->Write();
dataRatio->Write();
delete mcRatio;
delete dataRatio;
}
//********************************************************************
void Measurement1D::WriteShapePlot() {
//********************************************************************
TH1D* mcShape = (TH1D*)fMCHist->Clone((fName + "_MC_SHAPE").c_str());
TH1D* dataShape = (TH1D*)fDataHist->Clone((fName + "_data_SHAPE").c_str());
if (fShapeCovar) StatUtils::SetDataErrorFromCov(dataShape, fShapeCovar, 1E-38);
double shapeScale = 1.0;
if (fIsRawEvents) {
shapeScale = fDataHist->Integral() / fMCHist->Integral();
} else {
shapeScale = fDataHist->Integral("width") / fMCHist->Integral("width");
}
mcShape->Scale(shapeScale);
std::stringstream ss;
ss << shapeScale;
mcShape->SetTitle(ss.str().c_str());
mcShape->SetLineWidth(3);
mcShape->SetLineStyle(7);
mcShape->Write();
dataShape->Write();
delete mcShape;
}
//********************************************************************
void Measurement1D::WriteShapeRatioPlot() {
//********************************************************************
// Get a mcshape histogram
TH1D* mcShape = (TH1D*)fMCHist->Clone((fName + "_MC_SHAPE").c_str());
double shapeScale = 1.0;
if (fIsRawEvents) {
shapeScale = fDataHist->Integral() / fMCHist->Integral();
} else {
shapeScale = fDataHist->Integral("width") / fMCHist->Integral("width");
}
mcShape->Scale(shapeScale);
// Create shape ratio histograms
TH1D* mcShapeRatio = (TH1D*)mcShape->Clone((fName + "_MC_SHAPE_RATIO").c_str());
TH1D* dataShapeRatio = (TH1D*)fDataHist->Clone((fName + "_data_SHAPE_RATIO").c_str());
// Divide the histograms
mcShapeRatio->Divide(mcShape);
dataShapeRatio->Divide(mcShape);
// Colour the shape ratio plots
mcShapeRatio->SetLineWidth(3);
mcShapeRatio->SetLineStyle(7);
mcShapeRatio->Write();
dataShapeRatio->Write();
delete mcShapeRatio;
delete dataShapeRatio;
}
//// CRAP TO BE REMOVED
//********************************************************************
void Measurement1D::SetupMeasurement(std::string inputfile, std::string type,
FitWeight * rw, std::string fkdt) {
//********************************************************************
nuiskey samplekey = Config::CreateKey("sample");
samplekey.Set("name", fName);
samplekey.Set("type",type);
samplekey.Set("input",inputfile);
fSettings = LoadSampleSettings(samplekey);
// Reset everything to NULL
// Init();
// Check if name contains Evt, indicating that it is a raw number of events
// measurements and should thus be treated as once
fIsRawEvents = false;
if ((fName.find("Evt") != std::string::npos) && fIsRawEvents == false) {
fIsRawEvents = true;
LOG(SAM) << "Found event rate measurement but fIsRawEvents == false!"
<< std::endl;
LOG(SAM) << "Overriding this and setting fIsRawEvents == true!"
<< std::endl;
}
fIsEnu1D = false;
if (fName.find("XSec_1DEnu") != std::string::npos) {
fIsEnu1D = true;
LOG(SAM) << "::" << fName << "::" << std::endl;
LOG(SAM) << "Found XSec Enu measurement, applying flux integrated scaling, "
"not flux averaged!"
<< std::endl;
}
if (fIsEnu1D && fIsRawEvents) {
LOG(SAM) << "Found 1D Enu XSec distribution AND fIsRawEvents, is this "
"really correct?!"
<< std::endl;
LOG(SAM) << "Check experiment constructor for " << fName
<< " and correct this!" << std::endl;
LOG(SAM) << "I live in " << __FILE__ << ":" << __LINE__ << std::endl;
exit(-1);
}
fRW = rw;
if (!fInput and !fIsJoint) SetupInputs(inputfile);
// Set Default Options
SetFitOptions(fDefaultTypes);
// Set Passed Options
SetFitOptions(type);
// Still adding support for flat flux inputs
// // Set Enu Flux Scaling
// if (isFlatFluxFolding) this->Input()->ApplyFluxFolding(
// this->defaultFluxHist );
// FinaliseMeasurement();
}
//********************************************************************
void Measurement1D::SetupDefaultHist() {
//********************************************************************
// Setup fMCHist
fMCHist = (TH1D*)fDataHist->Clone();
fMCHist->SetNameTitle((fName + "_MC").c_str(),
(fName + "_MC" + fPlotTitles).c_str());
// Setup fMCFine
Int_t nBins = fMCHist->GetNbinsX();
fMCFine = new TH1D(
(fName + "_MC_FINE").c_str(), (fName + "_MC_FINE" + fPlotTitles).c_str(),
nBins * 6, fMCHist->GetBinLowEdge(1), fMCHist->GetBinLowEdge(nBins + 1));
fMCStat = (TH1D*)fMCHist->Clone();
fMCStat->Reset();
fMCHist->Reset();
fMCFine->Reset();
// Setup the NEUT Mode Array
PlotUtils::CreateNeutModeArray((TH1D*)fMCHist, (TH1**)fMCHist_PDG);
PlotUtils::ResetNeutModeArray((TH1**)fMCHist_PDG);
// Setup bin masks using sample name
if (fIsMask) {
std::string maskloc = FitPar::Config().GetParDIR(fName + ".mask");
if (maskloc.empty()) {
maskloc = FitPar::GetDataBase() + "/masks/" + fName + ".mask";
}
SetBinMask(maskloc);
}
fMCHist_Modes = new TrueModeStack( (fName + "_MODES").c_str(), ("True Channels"), fMCHist);
SetAutoProcessTH1(fMCHist_Modes, kCMD_Reset, kCMD_Norm, kCMD_Write);
return;
}
//********************************************************************
void Measurement1D::SetDataValues(std::string dataFile) {
//********************************************************************
// Override this function if the input file isn't in a suitable format
LOG(SAM) << "Reading data from: " << dataFile.c_str() << std::endl;
fDataHist =
PlotUtils::GetTH1DFromFile(dataFile, (fName + "_data"), fPlotTitles);
fDataTrue = (TH1D*)fDataHist->Clone();
// Number of data points is number of bins
fNDataPointsX = fDataHist->GetXaxis()->GetNbins();
return;
};
//********************************************************************
void Measurement1D::SetDataFromDatabase(std::string inhistfile,
std::string histname) {
//********************************************************************
LOG(SAM) << "Filling histogram from " << inhistfile << "->" << histname
<< std::endl;
fDataHist = PlotUtils::GetTH1DFromRootFile(
(GeneralUtils::GetTopLevelDir() + "/data/" + inhistfile), histname);
fDataHist->SetNameTitle((fName + "_data").c_str(), (fName + "_data").c_str());
return;
};
//********************************************************************
void Measurement1D::SetDataFromFile(std::string inhistfile,
std::string histname) {
//********************************************************************
LOG(SAM) << "Filling histogram from " << inhistfile << "->" << histname
<< std::endl;
fDataHist = PlotUtils::GetTH1DFromRootFile((inhistfile), histname);
fDataHist->SetNameTitle((fName + "_data").c_str(), (fName + "_data").c_str());
return;
};
//********************************************************************
void Measurement1D::SetCovarMatrix(std::string covarFile) {
//********************************************************************
// Covariance function, only really used when reading in the MB Covariances.
TFile* tempFile = new TFile(covarFile.c_str(), "READ");
TH2D* covarPlot = new TH2D();
// TH2D* decmpPlot = new TH2D();
TH2D* covarInvPlot = new TH2D();
TH2D* fFullCovarPlot = new TH2D();
std::string covName = "";
std::string covOption = FitPar::Config().GetParS("thrown_covariance");
if (fIsShape || fIsFree) covName = "shp_";
if (fIsDiag)
covName += "diag";
else
covName += "full";
covarPlot = (TH2D*)tempFile->Get((covName + "cov").c_str());
covarInvPlot = (TH2D*)tempFile->Get((covName + "covinv").c_str());
if (!covOption.compare("SUB"))
fFullCovarPlot = (TH2D*)tempFile->Get((covName + "cov").c_str());
else if (!covOption.compare("FULL"))
fFullCovarPlot = (TH2D*)tempFile->Get("fullcov");
else
ERR(WRN) << "Incorrect thrown_covariance option in parameters."
<< std::endl;
int dim = int(fDataHist->GetNbinsX()); //-this->masked->Integral());
int covdim = int(fDataHist->GetNbinsX());
this->covar = new TMatrixDSym(dim);
fFullCovar = new TMatrixDSym(dim);
fDecomp = new TMatrixDSym(dim);
int row, column = 0;
row = 0;
column = 0;
for (Int_t i = 0; i < covdim; i++) {
// if (this->masked->GetBinContent(i+1) > 0) continue;
for (Int_t j = 0; j < covdim; j++) {
// if (this->masked->GetBinContent(j+1) > 0) continue;
(*this->covar)(row, column) = covarPlot->GetBinContent(i + 1, j + 1);
(*fFullCovar)(row, column) = fFullCovarPlot->GetBinContent(i + 1, j + 1);
column++;
}
column = 0;
row++;
}
// Set bin errors on data
if (!fIsDiag) {
StatUtils::SetDataErrorFromCov(fDataHist, fFullCovar);
}
// Get Deteriminant and inverse matrix
// fCovDet = this->covar->Determinant();
TDecompSVD LU = TDecompSVD(*this->covar);
this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), "");
return;
};
//********************************************************************
// Sets the covariance matrix from a provided file in a text format
// scale is a multiplicative pre-factor to apply in the case where the
// covariance is given in some unit (e.g. 1E-38)
void Measurement1D::SetCovarMatrixFromText(std::string covarFile, int dim,
double scale) {
//********************************************************************
// Make a counter to track the line number
int row = 0;
std::string line;
- std::ifstream covarread(covarFile.c_str(), ifstream::in);
+ std::ifstream covarread(covarFile.c_str(), std::ifstream::in);
this->covar = new TMatrixDSym(dim);
fFullCovar = new TMatrixDSym(dim);
if (covarread.is_open())
LOG(SAM) << "Reading covariance matrix from file: " << covarFile
<< std::endl;
else
ERR(FTL) << "Covariance matrix provided is incorrect: " << covarFile
<< std::endl;
// Loop over the lines in the file
while (std::getline(covarread >> std::ws, line, '\n')) {
int column = 0;
// Loop over entries and insert them into matrix
std::vector entries = GeneralUtils::ParseToDbl(line, " ");
if (entries.size() <= 1) {
ERR(WRN) << "SetCovarMatrixFromText -> Covariance matrix only has <= 1 "
"entries on this line: "
<< row << std::endl;
}
for (std::vector::iterator iter = entries.begin();
iter != entries.end(); iter++) {
(*covar)(row, column) = *iter;
(*fFullCovar)(row, column) = *iter;
column++;
}
row++;
}
covarread.close();
// Scale the actualy covariance matrix by some multiplicative factor
(*fFullCovar) *= scale;
// Robust matrix inversion method
TDecompSVD LU = TDecompSVD(*this->covar);
// THIS IS ACTUALLY THE INVERSE COVARIANCE MATRIXA AAAAARGH
delete this->covar;
this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), "");
// Now need to multiply by the scaling factor
// If the covariance
(*this->covar) *= 1. / (scale);
return;
};
//********************************************************************
void Measurement1D::SetCovarMatrixFromCorrText(std::string corrFile, int dim) {
//********************************************************************
// Make a counter to track the line number
int row = 0;
std::string line;
- std::ifstream corr(corrFile.c_str(), ifstream::in);
+ std::ifstream corr(corrFile.c_str(), std::ifstream::in);
this->covar = new TMatrixDSym(dim);
this->fFullCovar = new TMatrixDSym(dim);
if (corr.is_open())
LOG(SAM) << "Reading and converting correlation matrix from file: "
<< corrFile << std::endl;
else {
ERR(FTL) << "Correlation matrix provided is incorrect: " << corrFile
<< std::endl;
exit(-1);
}
while (std::getline(corr >> std::ws, line, '\n')) {
int column = 0;
// Loop over entries and insert them into matrix
// Multiply by the errors to get the covariance, rather than the correlation
// matrix
std::vector entries = GeneralUtils::ParseToDbl(line, " ");
for (std::vector::iterator iter = entries.begin();
iter != entries.end(); iter++) {
double val = (*iter) * this->fDataHist->GetBinError(row + 1) * 1E38 *
this->fDataHist->GetBinError(column + 1) * 1E38;
if (val == 0) {
ERR(FTL) << "Found a zero value in the covariance matrix, assuming "
"this is an error!"
<< std::endl;
exit(-1);
}
(*this->covar)(row, column) = val;
(*this->fFullCovar)(row, column) = val;
column++;
}
row++;
}
// Robust matrix inversion method
TDecompSVD LU = TDecompSVD(*this->covar);
delete this->covar;
this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), "");
return;
};
//********************************************************************
// FullUnits refers to if we have "real" unscaled units in the covariance matrix, e.g. 1E-76.
// If this is the case we need to scale it so that the chi2 contribution is correct
// NUISANCE internally assumes the covariance matrix has units of 1E76
void Measurement1D::SetCovarFromDataFile(std::string covarFile,
std::string covName, bool FullUnits) {
//********************************************************************
LOG(SAM) << "Getting covariance from " << covarFile << "->" << covName
<< std::endl;
TFile* tempFile = new TFile(covarFile.c_str(), "READ");
TH2D* covPlot = (TH2D*)tempFile->Get(covName.c_str());
covPlot->SetDirectory(0);
// Scale the covariance matrix if it comes in normal units
if (FullUnits) {
covPlot->Scale(1.E76);
}
int dim = covPlot->GetNbinsX();
fFullCovar = new TMatrixDSym(dim);
for (int i = 0; i < dim; i++) {
for (int j = 0; j < dim; j++) {
(*fFullCovar)(i, j) = covPlot->GetBinContent(i + 1, j + 1);
}
}
this->covar = (TMatrixDSym*)fFullCovar->Clone();
fDecomp = (TMatrixDSym*)fFullCovar->Clone();
TDecompSVD LU = TDecompSVD(*this->covar);
this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), "");
TDecompChol LUChol = TDecompChol(*fDecomp);
LUChol.Decompose();
fDecomp = new TMatrixDSym(dim, LU.GetU().GetMatrixArray(), "");
return;
};
// //********************************************************************
// void Measurement1D::SetBinMask(std::string maskFile) {
// //********************************************************************
// // Create a mask histogram.
// int nbins = fDataHist->GetNbinsX();
// fMaskHist =
// new TH1I((fName + "_fMaskHist").c_str(),
// (fName + "_fMaskHist; Bin; Mask?").c_str(), nbins, 0, nbins);
// std::string line;
-// std::ifstream mask(maskFile.c_str(), ifstream::in);
+// std::ifstream mask(maskFile.c_str(), std::ifstream::in);
// if (mask.is_open())
// LOG(SAM) << "Reading bin mask from file: " << maskFile << std::endl;
// else
// LOG(FTL) << " Cannot find mask file." << std::endl;
// while (std::getline(mask >> std::ws, line, '\n')) {
// std::vector entries = GeneralUtils::ParseToInt(line, " ");
// // Skip lines with poorly formatted lines
// if (entries.size() < 2) {
// LOG(WRN) << "Measurement1D::SetBinMask(), couldn't parse line: " << line
// << std::endl;
// continue;
// }
// // The first index should be the bin number, the second should be the mask
// // value.
// fMaskHist->SetBinContent(entries[0], entries[1]);
// }
// // Set masked data bins to zero
// PlotUtils::MaskBins(fDataHist, fMaskHist);
// return;
// }
// //********************************************************************
// void Measurement1D::GetBinContents(std::vector& cont,
// std::vector& err) {
// //********************************************************************
// // Return a vector of the main bin contents
// for (int i = 0; i < fMCHist->GetNbinsX(); i++) {
// cont.push_back(fMCHist->GetBinContent(i + 1));
// err.push_back(fMCHist->GetBinError(i + 1));
// }
// return;
// };
/*
XSec Functions
*/
// //********************************************************************
// void Measurement1D::SetFluxHistogram(std::string fluxFile, int minE, int
// maxE,
// double fluxNorm) {
// //********************************************************************
// // Note this expects the flux bins to be given in terms of MeV
// LOG(SAM) << "Reading flux from file: " << fluxFile << std::endl;
// TGraph f(fluxFile.c_str(), "%lg %lg");
// fFluxHist =
// new TH1D((fName + "_flux").c_str(), (fName + "; E_{#nu} (GeV)").c_str(),
// f.GetN() - 1, minE, maxE);
// Double_t* yVal = f.GetY();
// for (int i = 0; i < fFluxHist->GetNbinsX(); ++i)
// fFluxHist->SetBinContent(i + 1, yVal[i] * fluxNorm);
// };
// //********************************************************************
// double Measurement1D::TotalIntegratedFlux(std::string intOpt, double low,
// double high) {
// //********************************************************************
// if (fInput->GetType() == kGiBUU) {
// return 1.0;
// }
// // The default case of low = -9999.9 and high = -9999.9
// if (low == -9999.9) low = this->EnuMin;
// if (high == -9999.9) high = this->EnuMax;
// int minBin = fFluxHist->GetXaxis()->FindBin(low);
// int maxBin = fFluxHist->GetXaxis()->FindBin(high);
// // Get integral over custom range
// double integral = fFluxHist->Integral(minBin, maxBin + 1, intOpt.c_str());
// return integral;
// };
diff --git a/src/FitBase/Measurement2D.cxx b/src/FitBase/Measurement2D.cxx
index 87f361f..77424bc 100644
--- a/src/FitBase/Measurement2D.cxx
+++ b/src/FitBase/Measurement2D.cxx
@@ -1,1961 +1,1961 @@
// 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 "Measurement2D.h"
#include "TDecompChol.h"
//********************************************************************
Measurement2D::Measurement2D(void) {
//********************************************************************
covar = NULL;
fDecomp = NULL;
fFullCovar = NULL;
fMCHist = NULL;
fMCFine = NULL;
fDataHist = NULL;
fMCHist_X = NULL;
fMCHist_Y = NULL;
fDataHist_X = NULL;
fDataHist_Y = NULL;
fMaskHist = NULL;
fMapHist = NULL;
fDataOrig = NULL;
fDataTrue = NULL;
fMCWeighted = NULL;
fDefaultTypes = "FIX/FULL/CHI2";
fAllowedTypes =
"FIX,FREE,SHAPE/FULL,DIAG/CHI2/NORM/ENUCORR/Q2CORR/ENU1D/FITPROJX/"
"FITPROJY";
fIsFix = false;
fIsShape = false;
fIsFree = false;
fIsDiag = false;
fIsFull = false;
fAddNormPen = false;
fIsMask = false;
fIsChi2SVD = false;
fIsRawEvents = false;
fIsDifXSec = false;
fIsEnu = false;
// XSec Scalings
fScaleFactor = -1.0;
fCurrentNorm = 1.0;
// Histograms
fDataHist = NULL;
fDataTrue = NULL;
fMCHist = NULL;
fMCFine = NULL;
fMCWeighted = NULL;
fMaskHist = NULL;
// Covar
covar = NULL;
fFullCovar = NULL;
fCovar = NULL;
fInvert = NULL;
fDecomp = NULL;
// Fake Data
fFakeDataInput = "";
fFakeDataFile = NULL;
// Options
fDefaultTypes = "FIX/FULL/CHI2";
fAllowedTypes =
"FIX,FREE,SHAPE/FULL,DIAG/CHI2/NORM/ENUCORR/Q2CORR/ENU1D/MASK";
fIsFix = false;
fIsShape = false;
fIsFree = false;
fIsDiag = false;
fIsFull = false;
fAddNormPen = false;
fIsMask = false;
fIsChi2SVD = false;
fIsRawEvents = false;
fIsDifXSec = false;
fIsEnu1D = false;
// Inputs
fInput = NULL;
fRW = NULL;
// Extra Histograms
fMCHist_Modes = NULL;
}
//********************************************************************
Measurement2D::~Measurement2D(void) {
//********************************************************************
if (fDataHist) delete fDataHist;
if (fDataTrue) delete fDataTrue;
if (fMCHist) delete fMCHist;
if (fMCFine) delete fMCFine;
if (fMCWeighted) delete fMCWeighted;
if (fMaskHist) delete fMaskHist;
if (covar) delete covar;
if (fFullCovar) delete fFullCovar;
if (fCovar) delete fCovar;
if (fInvert) delete fInvert;
if (fDecomp) delete fDecomp;
}
//********************************************************************
void Measurement2D::FinaliseSampleSettings() {
//********************************************************************
MeasurementBase::FinaliseSampleSettings();
// Setup naming + renaming
fName = fSettings.GetName();
fSettings.SetS("originalname", fName);
if (fSettings.Has("rename")) {
fName = fSettings.GetS("rename");
fSettings.SetS("name", fName);
}
// Setup all other options
LOG(SAM) << "Finalising Sample Settings: " << fName << std::endl;
if ((fSettings.GetS("originalname").find("Evt") != std::string::npos)) {
fIsRawEvents = true;
LOG(SAM) << "Found event rate measurement but using poisson likelihoods."
<< std::endl;
}
if (fSettings.GetS("originalname").find("XSec_1DEnu") != std::string::npos) {
fIsEnu1D = true;
LOG(SAM) << "::" << fName << "::" << std::endl;
LOG(SAM) << "Found XSec Enu measurement, applying flux integrated scaling, "
<< "not flux averaged!" << std::endl;
}
if (fIsEnu1D && fIsRawEvents) {
LOG(SAM) << "Found 1D Enu XSec distribution AND fIsRawEvents, is this "
"really correct?!"
<< std::endl;
LOG(SAM) << "Check experiment constructor for " << fName
<< " and correct this!" << std::endl;
LOG(SAM) << "I live in " << __FILE__ << ":" << __LINE__ << std::endl;
exit(-1);
}
if (!fRW) fRW = FitBase::GetRW();
if (!fInput) SetupInputs(fSettings.GetS("input"));
// Setup options
SetFitOptions(fDefaultTypes); // defaults
SetFitOptions(fSettings.GetS("type")); // user specified
EnuMin = GeneralUtils::StrToDbl(fSettings.GetS("enu_min"));
EnuMax = GeneralUtils::StrToDbl(fSettings.GetS("enu_max"));
if (fAddNormPen) {
fNormError = fSettings.GetNormError();
if (fNormError <= 0.0) {
ERR(WRN) << "Norm error for class " << fName << " is 0.0!" << std::endl;
ERR(WRN) << "If you want to use it please add fNormError=VAL" << std::endl;
throw;
}
}
}
void Measurement2D::CreateDataHistogram(int dimx, double* binx, int dimy, double* biny) {
if (fDataHist) delete fDataHist;
LOG(SAM) << "Creating Data Histogram dim : " << dimx << " " << dimy << std::endl;
fDataHist = new TH2D( (fSettings.GetName() + "_data").c_str(), (fSettings.GetFullTitles()).c_str(),
dimx - 1, binx, dimy - 1, biny );
}
void Measurement2D::SetDataFromTextFile(std::string datfile) {
// fDataHist = PlotUtils::GetTH2DFromTextFile(datfile,"");
}
void Measurement2D::SetDataFromRootFile(std::string datfile, std::string histname) {
fDataHist = PlotUtils::GetTH2DFromRootFile(datfile, histname);
}
void Measurement2D::SetDataValuesFromTextFile(std::string datfile, TH2D* hist) {
LOG(SAM) << "Setting data values from text file" << std::endl;
if (!hist) hist = fDataHist;
// Read TH2D From textfile
TH2D* valhist = (TH2D*) hist->Clone();
valhist->Reset();
PlotUtils::Set2DHistFromText(datfile, valhist, 1.0, true);
LOG(SAM) << " -> Filling values from read hist." << std::endl;
for (int i = 0; i < valhist->GetNbinsX(); i++) {
for (int j = 0; j < valhist->GetNbinsY(); j++) {
hist->SetBinContent(i + 1, j + 1, valhist->GetBinContent(i + 1, j + 1));
}
}
LOG(SAM) << " --> Done" << std::endl;
}
void Measurement2D::SetDataErrorsFromTextFile(std::string datfile, TH2D* hist) {
LOG(SAM) << "Setting data errors from text file" << std::endl;
if (!hist) hist = fDataHist;
// Read TH2D From textfile
TH2D* valhist = (TH2D*) hist->Clone();
valhist->Reset();
PlotUtils::Set2DHistFromText(datfile, valhist, 1.0);
// Fill Errors
LOG(SAM) << " -> Filling errors from read hist." << std::endl;
for (int i = 0; i < valhist->GetNbinsX(); i++) {
for (int j = 0; j < valhist->GetNbinsY(); j++) {
hist->SetBinError(i + 1, j + 1, valhist->GetBinContent(i + 1, j + 1));
}
}
LOG(SAM) << " --> Done" << std::endl;
}
void Measurement2D::SetMapValuesFromText(std::string dataFile) {
TH2D* hist = fDataHist;
std::vector edgex;
std::vector edgey;
for (int i = 0; i <= hist->GetNbinsX(); i++) edgex.push_back(hist->GetXaxis()->GetBinLowEdge(i + 1));
for (int i = 0; i <= hist->GetNbinsY(); i++) edgey.push_back(hist->GetYaxis()->GetBinLowEdge(i + 1));
fMapHist = new TH2I((fName + "_map").c_str(), (fName + fPlotTitles).c_str(),
edgex.size() - 1, &edgex[0], edgey.size() - 1, &edgey[0]);
LOG(SAM) << "Reading map from: " << dataFile << std::endl;
PlotUtils::Set2DHistFromText(dataFile, fMapHist, 1.0);
}
//********************************************************************
void Measurement2D::SetPoissonErrors() {
//********************************************************************
if (!fDataHist) {
ERR(FTL) << "Need a data hist to setup possion errors! " << std::endl;
ERR(FTL) << "Setup Data First!" << std::endl;
throw;
}
for (int i = 0; i < fDataHist->GetNbinsX() + 1; i++) {
fDataHist->SetBinError(i + 1, sqrt(fDataHist->GetBinContent(i + 1)));
}
}
//********************************************************************
void Measurement2D::SetCovarFromDiagonal(TH2D* data) {
//********************************************************************
if (!data and fDataHist) {
data = fDataHist;
}
if (data) {
LOG(SAM) << "Setting diagonal covariance for: " << data->GetName() << std::endl;
fFullCovar = StatUtils::MakeDiagonalCovarMatrix(data);
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
} else {
ERR(FTL) << "No data input provided to set diagonal covar from!" << std::endl;
}
// if (!fIsDiag) {
// ERR(FTL) << "SetCovarMatrixFromDiag called for measurement "
// << "that is not set as diagonal." << std::endl;
// throw;
// }
}
//********************************************************************
void Measurement2D::SetCovarFromTextFile(std::string covfile, int dim) {
//********************************************************************
if (dim == -1) {
dim = this->GetNDOF();
}
LOG(SAM) << "Reading covariance from text file: " << covfile << " " << dim << std::endl;
fFullCovar = StatUtils::GetCovarFromTextFile(covfile, dim);
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
//********************************************************************
void Measurement2D::SetCovarFromRootFile(std::string covfile, std::string histname) {
//********************************************************************
LOG(SAM) << "Reading covariance from text file: " << covfile << ";" << histname << std::endl;
fFullCovar = StatUtils::GetCovarFromRootFile(covfile, histname);
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
//********************************************************************
void Measurement2D::SetCovarInvertFromTextFile(std::string covfile, int dim) {
//********************************************************************
if (dim == -1) {
dim = this->GetNDOF();
}
LOG(SAM) << "Reading inverted covariance from text file: " << covfile << std::endl;
covar = StatUtils::GetCovarFromTextFile(covfile, dim);
fFullCovar = StatUtils::GetInvert(covar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
//********************************************************************
void Measurement2D::SetCovarInvertFromRootFile(std::string covfile, std::string histname) {
//********************************************************************
LOG(SAM) << "Reading inverted covariance from text file: " << covfile << ";" << histname << std::endl;
covar = StatUtils::GetCovarFromRootFile(covfile, histname);
fFullCovar = StatUtils::GetInvert(covar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
//********************************************************************
void Measurement2D::SetCorrelationFromTextFile(std::string covfile, int dim) {
//********************************************************************
if (dim == -1) dim = this->GetNDOF();
LOG(SAM) << "Reading data correlations from text file: " << covfile << ";" << dim << std::endl;
TMatrixDSym* correlation = StatUtils::GetCovarFromTextFile(covfile, dim);
if (!fDataHist) {
ERR(FTL) << "Trying to set correlations from text file but there is no data to build it from. \n"
<< "In constructor make sure data is set before SetCorrelationFromTextFile is called. \n" << std::endl;
throw;
}
// Fill covar from data errors and correlations
fFullCovar = new TMatrixDSym(dim);
for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
for (int j = 0; j < fDataHist->GetNbinsX(); j++) {
(*fFullCovar)(i, j) = (*correlation)(i, j) * fDataHist->GetBinError(i + 1) * fDataHist->GetBinError(j + 1) * 1.E76;
}
}
// Fill other covars.
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
delete correlation;
}
//********************************************************************
void Measurement2D::SetCorrelationFromRootFile(std::string covfile, std::string histname) {
//********************************************************************
LOG(SAM) << "Reading data correlations from text file: " << covfile << ";" << histname << std::endl;
TMatrixDSym* correlation = StatUtils::GetCovarFromRootFile(covfile, histname);
if (!fDataHist) {
ERR(FTL) << "Trying to set correlations from text file but there is no data to build it from. \n"
<< "In constructor make sure data is set before SetCorrelationFromTextFile is called. \n" << std::endl;
throw;
}
// Fill covar from data errors and correlations
fFullCovar = new TMatrixDSym(fDataHist->GetNbinsX());
for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
for (int j = 0; j < fDataHist->GetNbinsX(); j++) {
(*fFullCovar)(i, j) = (*correlation)(i, j) * fDataHist->GetBinError(i + 1) * fDataHist->GetBinError(j + 1) * 1.E76;
}
}
// Fill other covars.
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
delete correlation;
}
//********************************************************************
void Measurement2D::SetCholDecompFromTextFile(std::string covfile, int dim) {
//********************************************************************
if (dim == -1) {
dim = this->GetNDOF();
}
LOG(SAM) << "Reading cholesky from text file: " << covfile << " " << dim << std::endl;
TMatrixD* temp = StatUtils::GetMatrixFromTextFile(covfile, dim, dim);
TMatrixD* trans = (TMatrixD*)temp->Clone();
trans->T();
(*trans) *= (*temp);
fFullCovar = new TMatrixDSym(dim, trans->GetMatrixArray(), "");
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
delete temp;
delete trans;
}
//********************************************************************
void Measurement2D::SetCholDecompFromRootFile(std::string covfile, std::string histname) {
//********************************************************************
LOG(SAM) << "Reading cholesky decomp from root file: " << covfile << ";" << histname << std::endl;
TMatrixD* temp = StatUtils::GetMatrixFromRootFile(covfile, histname);
TMatrixD* trans = (TMatrixD*)temp->Clone();
trans->T();
(*trans) *= (*temp);
fFullCovar = new TMatrixDSym(temp->GetNrows(), trans->GetMatrixArray(), "");
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
delete temp;
delete trans;
}
//********************************************************************
void Measurement2D::ScaleData(double scale) {
//********************************************************************
fDataHist->Scale(scale);
}
//********************************************************************
void Measurement2D::ScaleDataErrors(double scale) {
//********************************************************************
for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
for (int j = 0; j < fDataHist->GetNbinsY(); j++) {
fDataHist->SetBinError(i + 1, j + 1, fDataHist->GetBinError(i + 1, j + 1) * scale);
}
}
}
//********************************************************************
void Measurement2D::ScaleCovar(double scale) {
//********************************************************************
(*fFullCovar) *= scale;
(*covar) *= 1.0 / scale;
(*fDecomp) *= sqrt(scale);
}
//********************************************************************
void Measurement2D::SetBinMask(std::string maskfile) {
//********************************************************************
if (!fIsMask) return;
LOG(SAM) << "Reading bin mask from file: " << maskfile << std::endl;
// Create a mask histogram with dim of data
int nbinsx = fDataHist->GetNbinsX();
int nbinxy = fDataHist->GetNbinsY();
fMaskHist =
new TH2I((fSettings.GetName() + "_BINMASK").c_str(),
(fSettings.GetName() + "_BINMASK; Bin; Mask?").c_str(), nbinsx, 0, nbinsx, nbinxy, 0, nbinxy);
std::string line;
- std::ifstream mask(maskfile.c_str(), ifstream::in);
+ std::ifstream mask(maskfile.c_str(), std::ifstream::in);
if (!mask.is_open()) {
LOG(FTL) << " Cannot find mask file." << std::endl;
throw;
}
while (std::getline(mask >> std::ws, line, '\n')) {
std::vector entries = GeneralUtils::ParseToInt(line, " ");
// Skip lines with poorly formatted lines
if (entries.size() < 2) {
LOG(WRN) << "Measurement2D::SetBinMask(), couldn't parse line: " << line
<< std::endl;
continue;
}
// The first index should be the bin number, the second should be the mask
// value.
int val = 0;
if (entries[2] > 0) val = 1;
fMaskHist->SetBinContent(entries[0], entries[1], val);
}
// Apply masking by setting masked data bins to zero
PlotUtils::MaskBins(fDataHist, fMaskHist);
return;
}
//********************************************************************
void Measurement2D::FinaliseMeasurement() {
//********************************************************************
LOG(SAM) << "Finalising Measurement: " << fName << std::endl;
if (fSettings.GetB("onlymc")) {
if (fDataHist) delete fDataHist;
fDataHist = new TH2D("empty_data", "empty_data", 1, 0.0, 1.0,1,0.0,1.0);
}
// Make sure data is setup
if (!fDataHist) {
ERR(FTL) << "No data has been setup inside " << fName << " constructor!" << std::endl;
throw;
}
// Make sure covariances are setup
if (!fFullCovar) {
fIsDiag = true;
SetCovarFromDiagonal(fDataHist);
}
if (!covar) {
covar = StatUtils::GetInvert(fFullCovar);
}
if (!fDecomp) {
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
// Setup fMCHist from data
fMCHist = (TH2D*)fDataHist->Clone();
fMCHist->SetNameTitle((fSettings.GetName() + "_MC").c_str(),
(fSettings.GetFullTitles()).c_str());
fMCHist->Reset();
// Setup fMCFine
fMCFine = new TH2D("mcfine", "mcfine", fDataHist->GetNbinsX() * 6,
fMCHist->GetXaxis()->GetBinLowEdge(1),
fMCHist->GetXaxis()->GetBinLowEdge(fDataHist->GetNbinsX() + 1),
fDataHist->GetNbinsY() * 6,
fMCHist->GetYaxis()->GetBinLowEdge(1),
fMCHist->GetYaxis()->GetBinLowEdge(fDataHist->GetNbinsY() + 1));
fMCFine->SetNameTitle((fSettings.GetName() + "_MC_FINE").c_str(),
(fSettings.GetFullTitles()).c_str());
fMCFine->Reset();
// Setup MC Stat
fMCStat = (TH2D*)fMCHist->Clone();
fMCStat->Reset();
// Search drawopts for possible types to include by default
std::string drawopts = FitPar::Config().GetParS("drawopts");
if (drawopts.find("MODES") != std::string::npos) {
fMCHist_Modes = new TrueModeStack( (fSettings.GetName() + "_MODES").c_str(),
("True Channels"), fMCHist);
SetAutoProcessTH1(fMCHist_Modes);
}
// Setup bin masks using sample name
if (fIsMask) {
std::string curname = fName;
std::string origname = fSettings.GetS("originalname");
// Check rename.mask
std::string maskloc = FitPar::Config().GetParDIR(curname + ".mask");
// Check origname.mask
if (maskloc.empty()) maskloc = FitPar::Config().GetParDIR(origname + ".mask");
// Check database
if (maskloc.empty()) {
maskloc = FitPar::GetDataBase() + "/masks/" + origname + ".mask";
}
// Setup Bin Mask
SetBinMask(maskloc);
}
if (fScaleFactor < 0) {
ERR(FTL) << "I found a negative fScaleFactor in " << __FILE__ << ":" << __LINE__ << std::endl;
ERR(FTL) << "fScaleFactor = " << fScaleFactor << std::endl;
ERR(FTL) << "EXITING" << std::endl;
throw;
}
// Create and fill Weighted Histogram
if (!fMCWeighted) {
fMCWeighted = (TH2D*)fMCHist->Clone();
fMCWeighted->SetNameTitle((fName + "_MCWGHTS").c_str(),
(fName + "_MCWGHTS" + fPlotTitles).c_str());
fMCWeighted->GetYaxis()->SetTitle("Weighted Events");
}
}
//********************************************************************
void Measurement2D::SetFitOptions(std::string opt) {
//********************************************************************
// Do nothing if default given
if (opt == "DEFAULT") return;
// CHECK Conflicting Fit Options
std::vector fit_option_allow =
GeneralUtils::ParseToStr(fAllowedTypes, "/");
for (UInt_t i = 0; i < fit_option_allow.size(); i++) {
std::vector fit_option_section =
GeneralUtils::ParseToStr(fit_option_allow.at(i), ",");
bool found_option = false;
for (UInt_t j = 0; j < fit_option_section.size(); j++) {
std::string av_opt = fit_option_section.at(j);
if (!found_option and opt.find(av_opt) != std::string::npos) {
found_option = true;
} else if (found_option and opt.find(av_opt) != std::string::npos) {
ERR(FTL) << "ERROR: Conflicting fit options provided: "
<< opt << std::endl
<< "Conflicting group = " << fit_option_section.at(i) << std::endl
<< "You should only supply one of these options in card file." << std::endl;
throw;
}
}
}
// Check all options are allowed
std::vector fit_options_input =
GeneralUtils::ParseToStr(opt, "/");
for (UInt_t i = 0; i < fit_options_input.size(); i++) {
if (fAllowedTypes.find(fit_options_input.at(i)) == std::string::npos) {
ERR(FTL) << "ERROR: Fit Option '" << fit_options_input.at(i)
<< "' Provided is not allowed for this measurement."
<< std::endl;
ERR(FTL) << "Fit Options should be provided as a '/' seperated list "
"(e.g. FREE/DIAG/NORM)"
<< std::endl;
ERR(FTL) << "Available options for " << fName << " are '" << fAllowedTypes
<< "'" << std::endl;
throw;
}
}
// Set TYPE
fFitType = opt;
// FIX,SHAPE,FREE
if (opt.find("FIX") != std::string::npos) {
fIsFree = fIsShape = false;
fIsFix = true;
} else if (opt.find("SHAPE") != std::string::npos) {
fIsFree = fIsFix = false;
fIsShape = true;
} else if (opt.find("FREE") != std::string::npos) {
fIsFix = fIsShape = false;
fIsFree = true;
}
// DIAG,FULL (or default to full)
if (opt.find("DIAG") != std::string::npos) {
fIsDiag = true;
fIsFull = false;
} else if (opt.find("FULL") != std::string::npos) {
fIsDiag = false;
fIsFull = true;
}
// CHI2/LL (OTHERS?)
if (opt.find("LOG") != std::string::npos) {
fIsChi2 = false;
ERR(FTL) << "No other LIKELIHOODS properly supported!" << std::endl;
ERR(FTL) << "Try to use a chi2!" << std::endl;
throw;
} else {
fIsChi2 = true;
}
// EXTRAS
if (opt.find("RAW") != std::string::npos) fIsRawEvents = true;
if (opt.find("DIF") != std::string::npos) fIsDifXSec = true;
if (opt.find("ENU1D") != std::string::npos) fIsEnu1D = true;
if (opt.find("NORM") != std::string::npos) fAddNormPen = true;
if (opt.find("MASK") != std::string::npos) fIsMask = true;
// Set TYPE
fFitType = opt;
// FIX,SHAPE,FREE
if (opt.find("FIX") != std::string::npos) {
fIsFree = fIsShape = false;
fIsFix = true;
} else if (opt.find("SHAPE") != std::string::npos) {
fIsFree = fIsFix = false;
fIsShape = true;
} else if (opt.find("FREE") != std::string::npos) {
fIsFix = fIsShape = false;
fIsFree = true;
}
// DIAG,FULL (or default to full)
if (opt.find("DIAG") != std::string::npos) {
fIsDiag = true;
fIsFull = false;
} else if (opt.find("FULL") != std::string::npos) {
fIsDiag = false;
fIsFull = true;
}
// CHI2/LL (OTHERS?)
if (opt.find("LOG") != std::string::npos)
fIsChi2 = false;
else
fIsChi2 = true;
// EXTRAS
if (opt.find("RAW") != std::string::npos) fIsRawEvents = true;
if (opt.find("DIF") != std::string::npos) fIsDifXSec = true;
if (opt.find("ENU1D") != std::string::npos) fIsEnu = true;
if (opt.find("NORM") != std::string::npos) fAddNormPen = true;
if (opt.find("MASK") != std::string::npos) fIsMask = true;
fIsProjFitX = (opt.find("FITPROJX") != std::string::npos);
fIsProjFitY = (opt.find("FITPROJY") != std::string::npos);
return;
};
/*
Reconfigure LOOP
*/
//********************************************************************
void Measurement2D::ResetAll() {
//********************************************************************
fMCHist->Reset();
fMCFine->Reset();
fMCStat->Reset();
return;
};
//********************************************************************
void Measurement2D::FillHistograms() {
//********************************************************************
if (Signal) {
fMCHist->Fill(fXVar, fYVar, Weight);
fMCFine->Fill(fXVar, fYVar, Weight);
fMCStat->Fill(fXVar, fYVar, 1.0);
if (fMCHist_Modes) fMCHist_Modes->Fill(Mode, fXVar, fYVar, Weight);
}
return;
};
//********************************************************************
void Measurement2D::ScaleEvents() {
//********************************************************************
// Fill MCWeighted;
// for (int i = 0; i < fMCHist->GetNbinsX(); i++) {
// fMCWeighted->SetBinContent(i + 1, fMCHist->GetBinContent(i + 1));
// fMCWeighted->SetBinError(i + 1, fMCHist->GetBinError(i + 1));
// }
// Setup Stat ratios for MC and MC Fine
double* statratio = new double[fMCHist->GetNbinsX()];
for (int i = 0; i < fMCHist->GetNbinsX(); i++) {
if (fMCHist->GetBinContent(i + 1) != 0) {
statratio[i] = fMCHist->GetBinError(i + 1) / fMCHist->GetBinContent(i + 1);
} else {
statratio[i] = 0.0;
}
}
double* statratiofine = new double[fMCFine->GetNbinsX()];
for (int i = 0; i < fMCFine->GetNbinsX(); i++) {
if (fMCFine->GetBinContent(i + 1) != 0) {
statratiofine[i] = fMCFine->GetBinError(i + 1) / fMCFine->GetBinContent(i + 1);
} else {
statratiofine[i] = 0.0;
}
}
// Scaling for raw event rates
if (fIsRawEvents) {
double datamcratio = fDataHist->Integral() / fMCHist->Integral();
fMCHist->Scale(datamcratio);
fMCFine->Scale(datamcratio);
if (fMCHist_Modes) fMCHist_Modes->Scale(datamcratio);
// Scaling for XSec as function of Enu
} else if (fIsEnu1D) {
PlotUtils::FluxUnfoldedScaling(fMCHist, GetFluxHistogram(),
GetEventHistogram(), fScaleFactor);
PlotUtils::FluxUnfoldedScaling(fMCFine, GetFluxHistogram(),
GetEventHistogram(), fScaleFactor);
// if (fMCHist_Modes) {
// PlotUtils::FluxUnfoldedScaling(fMCHist_Modes, GetFluxHistogram(),
// GetEventHistogram(), fScaleFactor,
// fNEvents);
// }
// Any other differential scaling
} else {
fMCHist->Scale(fScaleFactor, "width");
fMCFine->Scale(fScaleFactor, "width");
// if (fMCHist_Modes) fMCHist_Modes->Scale(fScaleFactor, "width");
}
// Proper error scaling - ROOT Freaks out with xsec weights sometimes
for (int i = 0; i < fMCStat->GetNbinsX(); i++) {
fMCHist->SetBinError(i + 1, fMCHist->GetBinContent(i + 1) * statratio[i]);
}
for (int i = 0; i < fMCFine->GetNbinsX(); i++) {
fMCFine->SetBinError(i + 1, fMCFine->GetBinContent(i + 1) * statratiofine[i]);
}
// Clean up
delete statratio;
delete statratiofine;
return;
};
//********************************************************************
void Measurement2D::ApplyNormScale(double norm) {
//********************************************************************
fCurrentNorm = norm;
fMCHist->Scale(1.0 / norm);
fMCFine->Scale(1.0 / norm);
return;
};
/*
Statistic Functions - Outsources to StatUtils
*/
//********************************************************************
int Measurement2D::GetNDOF() {
//********************************************************************
// Just incase it has gone...
if (!fDataHist) return -1;
int nDOF = 0;
// If datahist has no errors make sure we don't include those bins as they are
// not data points
for (int xBin = 0; xBin < fDataHist->GetNbinsX() + 1; ++xBin) {
for (int yBin = 0; yBin < fDataHist->GetNbinsY() + 1; ++yBin) {
if (fDataHist->GetBinError(xBin, yBin) != 0)
++nDOF;
}
}
// Account for possible bin masking
int nMasked = 0;
if (fMaskHist and fIsMask)
if (fMaskHist->Integral() > 0)
for (int xBin = 0; xBin < fMaskHist->GetNbinsX() + 1; ++xBin)
for (int yBin = 0; yBin < fMaskHist->GetNbinsY() + 1; ++yBin)
if (fMaskHist->GetBinContent(xBin, yBin) > 0.5) ++nMasked;
// Take away those masked DOF
if (fIsMask) {
nDOF -= nMasked;
}
return nDOF;
}
//********************************************************************
double Measurement2D::GetLikelihood() {
//********************************************************************
// If this is for a ratio, there is no data histogram to compare to!
if (fNoData || !fDataHist) return 0.;
// Fix weird masking bug
if (!fIsMask) {
if (fMaskHist) {
fMaskHist = NULL;
}
} else {
if (fMaskHist) {
PlotUtils::MaskBins(fMCHist, fMaskHist);
}
}
// if (fIsProjFitX or fIsProjFitY) return GetProjectedChi2();
// Scale up the results to match each other (Not using width might be
// inconsistent with Meas1D)
double scaleF = fDataHist->Integral() / fMCHist->Integral();
if (fIsShape) {
fMCHist->Scale(scaleF);
fMCFine->Scale(scaleF);
//PlotUtils::ScaleNeutModeArray((TH1**)fMCHist_PDG, scaleF);
}
if (!fMapHist) {
fMapHist = StatUtils::GenerateMap(fDataHist);
}
// Get the chi2 from either covar or diagonals
double chi2 = 0.0;
if (fIsChi2) {
if (fIsDiag) {
chi2 =
StatUtils::GetChi2FromDiag(fDataHist, fMCHist, fMapHist, fMaskHist);
} else {
chi2 = StatUtils::GetChi2FromCov(fDataHist, fMCHist, covar, fMapHist,
fMaskHist);
}
}
// Add a normal penalty term
if (fAddNormPen) {
chi2 +=
(1 - (fCurrentNorm)) * (1 - (fCurrentNorm)) / (fNormError * fNormError);
LOG(REC) << "Norm penalty = "
<< (1 - (fCurrentNorm)) * (1 - (fCurrentNorm)) /
(fNormError * fNormError)
<< std::endl;
}
// Adjust the shape back to where it was.
if (fIsShape and !FitPar::Config().GetParB("saveshapescaling")) {
fMCHist->Scale(1. / scaleF);
fMCFine->Scale(1. / scaleF);
}
fLikelihood = chi2;
return chi2;
}
/*
Fake Data Functions
*/
//********************************************************************
void Measurement2D::SetFakeDataValues(std::string fakeOption) {
//********************************************************************
// Setup original/datatrue
TH2D* tempdata = (TH2D*) fDataHist->Clone();
if (!fIsFakeData) {
fIsFakeData = true;
// Make a copy of the original data histogram.
if (!fDataOrig) fDataOrig = (TH2D*)fDataHist->Clone((fName + "_data_original").c_str());
} else {
ResetFakeData();
}
// Setup Inputs
fFakeDataInput = fakeOption;
LOG(SAM) << "Setting fake data from : " << fFakeDataInput << std::endl;
// From MC
if (fFakeDataInput.compare("MC") == 0) {
fDataHist = (TH2D*)fMCHist->Clone((fName + "_MC").c_str());
// Fake File
} else {
if (!fFakeDataFile) fFakeDataFile = new TFile(fFakeDataInput.c_str(), "READ");
fDataHist = (TH2D*)fFakeDataFile->Get((fName + "_MC").c_str());
}
// Setup Data Hist
fDataHist->SetNameTitle((fName + "_FAKE").c_str(),
(fName + fPlotTitles).c_str());
// Replace Data True
if (fDataTrue) delete fDataTrue;
fDataTrue = (TH2D*)fDataHist->Clone();
fDataTrue->SetNameTitle((fName + "_FAKE_TRUE").c_str(),
(fName + fPlotTitles).c_str());
// Make a new covariance for fake data hist.
int nbins = fDataHist->GetNbinsX() * fDataHist->GetNbinsY();
double alpha_i = 0.0;
double alpha_j = 0.0;
for (int i = 0; i < nbins; i++) {
for (int j = 0; j < nbins; j++) {
if (tempdata->GetBinContent(i + 1) && tempdata->GetBinContent(j + 1)) {
alpha_i = fDataHist->GetBinContent(i + 1) / tempdata->GetBinContent(i + 1);
alpha_j = fDataHist->GetBinContent(j + 1) / tempdata->GetBinContent(j + 1);
} else {
alpha_i = 0.0;
alpha_j = 0.0;
}
(*fFullCovar)(i, j) = alpha_i * alpha_j * (*fFullCovar)(i, j);
}
}
// Setup Covariances
if (covar) delete covar;
covar = StatUtils::GetInvert(fFullCovar);
if (fDecomp) delete fDecomp;
fDecomp = StatUtils::GetInvert(fFullCovar);
delete tempdata;
return;
};
//********************************************************************
void Measurement2D::ResetFakeData() {
//********************************************************************
if (fIsFakeData) {
if (fDataHist) delete fDataHist;
fDataHist = (TH2D*)fDataTrue->Clone((fSettings.GetName() + "_FKDAT").c_str());
}
}
//********************************************************************
void Measurement2D::ResetData() {
//********************************************************************
if (fIsFakeData) {
if (fDataHist) delete fDataHist;
fDataHist = (TH2D*)fDataOrig->Clone((fSettings.GetName() + "_data").c_str());
}
fIsFakeData = false;
}
//********************************************************************
void Measurement2D::ThrowCovariance() {
//********************************************************************
// Take a fDecomposition and use it to throw the current dataset.
// Requires fDataTrue also be set incase used repeatedly.
if (fDataHist) delete fDataHist;
fDataHist = StatUtils::ThrowHistogram(fDataTrue, fFullCovar);
return;
};
//********************************************************************
void Measurement2D::ThrowDataToy() {
//********************************************************************
if (!fDataTrue) fDataTrue = (TH2D*) fDataHist->Clone();
if (fMCHist) delete fMCHist;
fMCHist = StatUtils::ThrowHistogram(fDataTrue, fFullCovar);
}
/*
Access Functions
*/
//********************************************************************
TH2D* Measurement2D::GetMCHistogram() {
//********************************************************************
if (!fMCHist) return fMCHist;
std::ostringstream chi2;
chi2 << std::setprecision(5) << this->GetLikelihood();
int linecolor = kRed;
int linestyle = 1;
int linewidth = 1;
int fillcolor = 0;
int fillstyle = 1001;
if (fSettings.Has("linecolor")) linecolor = fSettings.GetI("linecolor");
if (fSettings.Has("linestyle")) linestyle = fSettings.GetI("linestyle");
if (fSettings.Has("linewidth")) linewidth = fSettings.GetI("linewidth");
if (fSettings.Has("fillcolor")) fillcolor = fSettings.GetI("fillcolor");
if (fSettings.Has("fillstyle")) fillstyle = fSettings.GetI("fillstyle");
fMCHist->SetTitle(chi2.str().c_str());
fMCHist->SetLineColor(linecolor);
fMCHist->SetLineStyle(linestyle);
fMCHist->SetLineWidth(linewidth);
fMCHist->SetFillColor(fillcolor);
fMCHist->SetFillStyle(fillstyle);
return fMCHist;
};
//********************************************************************
TH2D* Measurement2D::GetDataHistogram() {
//********************************************************************
if (!fDataHist) return fDataHist;
int datacolor = kBlack;
int datastyle = 1;
int datawidth = 1;
if (fSettings.Has("datacolor")) datacolor = fSettings.GetI("datacolor");
if (fSettings.Has("datastyle")) datastyle = fSettings.GetI("datastyle");
if (fSettings.Has("datawidth")) datawidth = fSettings.GetI("datawidth");
fDataHist->SetLineColor(datacolor);
fDataHist->SetLineWidth(datawidth);
fDataHist->SetMarkerStyle(datastyle);
return fDataHist;
};
/*
Write Functions
*/
// Save all the histograms at once
//********************************************************************
void Measurement2D::Write(std::string drawOpt) {
//********************************************************************
// Get Draw Options
drawOpt = FitPar::Config().GetParS("drawopts");
// Write Settigns
if (drawOpt.find("SETTINGS") != std::string::npos) {
fSettings.Set("#chi^{2}", fLikelihood);
fSettings.Set("NDOF", this->GetNDOF() );
fSettings.Set("#chi^{2}/NDOF", fLikelihood / this->GetNDOF() );
fSettings.Write();
}
// Write Data/MC
GetDataList().at(0)->Write();
GetMCList().at(0)->Write();
// Write Fine Histogram
if (drawOpt.find("FINE") != std::string::npos)
GetFineList().at(0)->Write();
// Write Weighted Histogram
if (drawOpt.find("WEIGHTS") != std::string::npos && fMCWeighted)
fMCWeighted->Write();
// Save Flux/Evt if no event manager
if (!FitPar::Config().GetParB("EventManager")) {
if (drawOpt.find("FLUX") != std::string::npos && GetFluxHistogram())
GetFluxHistogram()->Write();
if (drawOpt.find("EVT") != std::string::npos && GetEventHistogram())
GetEventHistogram()->Write();
if (drawOpt.find("XSEC") != std::string::npos && GetEventHistogram())
GetEventHistogram()->Write();
}
// Write Mask
if (fIsMask && (drawOpt.find("MASK") != std::string::npos)) {
fMaskHist->Write();
}
// Write Covariances
if (drawOpt.find("COV") != std::string::npos && fFullCovar) {
PlotUtils::GetFullCovarPlot(fFullCovar, fSettings.GetName());
}
if (drawOpt.find("INVCOV") != std::string::npos && covar) {
PlotUtils::GetInvCovarPlot(covar, fSettings.GetName());
}
if (drawOpt.find("DECOMP") != std::string::npos && fDecomp) {
PlotUtils::GetDecompCovarPlot(fDecomp, fSettings.GetName());
}
// // Likelihood residual plots
// if (drawOpt.find("RESIDUAL") != std::string::npos) {
// WriteResidualPlots();
// }
// // RATIO
// if (drawOpt.find("CANVMC") != std::string::npos) {
// TCanvas* c1 = WriteMCCanvas(fDataHist, fMCHist);
// c1->Write();
// delete c1;
// }
// // PDG
// if (drawOpt.find("CANVPDG") != std::string::npos && fMCHist_Modes) {
// TCanvas* c2 = WritePDGCanvas(fDataHist, fMCHist, fMCHist_Modes);
// c2->Write();
// delete c2;
// }
// Write Extra Histograms
AutoWriteExtraTH1();
WriteExtraHistograms();
/// 2D VERSION
// If null pointer return
if (!fMCHist and !fDataHist) {
LOG(SAM) << fName << "Incomplete histogram set!" << std::endl;
return;
}
// Config::Get().out->cd();
// Get Draw Options
drawOpt = FitPar::Config().GetParS("drawopts");
bool drawData = (drawOpt.find("DATA") != std::string::npos);
bool drawNormal = (drawOpt.find("MC") != std::string::npos);
bool drawEvents = (drawOpt.find("EVT") != std::string::npos);
bool drawXSec = (drawOpt.find("XSEC") != std::string::npos);
bool drawFine = (drawOpt.find("FINE") != std::string::npos);
bool drawRatio = (drawOpt.find("RATIO") != std::string::npos);
// bool drawModes = (drawOpt.find("MODES") != std::string::npos);
bool drawShape = (drawOpt.find("SHAPE") != std::string::npos);
bool residual = (drawOpt.find("RESIDUAL") != std::string::npos);
bool drawMatrix = (drawOpt.find("MATRIX") != std::string::npos);
bool drawFlux = (drawOpt.find("FLUX") != std::string::npos);
bool drawMask = (drawOpt.find("MASK") != std::string::npos);
bool drawMap = (drawOpt.find("MAP") != std::string::npos);
bool drawProj = (drawOpt.find("PROJ") != std::string::npos);
// bool drawCanvPDG = (drawOpt.find("CANVPDG") != std::string::npos);
bool drawCov = (drawOpt.find("COV") != std::string::npos);
bool drawSliceCanvYMC = (drawOpt.find("CANVYMC") != std::string::npos);
bool drawWeighted = (drawOpt.find("WGHT") != std::string::npos);
if (FitPar::Config().GetParB("EventManager")) {
drawFlux = false;
drawXSec = false;
drawEvents = false;
}
if (fMaskHist) fMaskHist->Write();
// Save standard plots
if (drawData) this->GetDataList().at(0)->Write();
if (drawNormal) this->GetMCList().at(0)->Write();
if (drawCov) {
TH2D(*fFullCovar).Write((fName + "_COV").c_str());
}
if (drawOpt.find("INVCOV") != std::string::npos) {
TH2D(*covar).Write((fName + "_INVCOV").c_str());
}
// Generate a simple map
if (!fMapHist) fMapHist = StatUtils::GenerateMap(fDataHist);
// Convert to 1D Lists
TH1D* data_1D = StatUtils::MapToTH1D(fDataHist, fMapHist);
TH1D* mc_1D = StatUtils::MapToTH1D(fMCHist, fMapHist);
TH1I* mask_1D = StatUtils::MapToMask(fMaskHist, fMapHist);
data_1D->Write();
mc_1D->Write();
if (mask_1D) {
mask_1D->Write();
TMatrixDSym* calc_cov =
StatUtils::ApplyInvertedMatrixMasking(covar, mask_1D);
TH1D* calc_data = StatUtils::ApplyHistogramMasking(data_1D, mask_1D);
TH1D* calc_mc = StatUtils::ApplyHistogramMasking(mc_1D, mask_1D);
TH2D* bin_cov = new TH2D(*calc_cov);
bin_cov->Write();
calc_data->Write();
calc_mc->Write();
delete mask_1D;
delete calc_cov;
delete calc_data;
delete calc_mc;
delete bin_cov;
}
delete data_1D;
delete mc_1D;
// Save only mc and data if splines
if (fEventType == 4 or fEventType == 3) {
return;
}
// Draw Extra plots
if (drawFine) this->GetFineList().at(0)->Write();
if (drawFlux and GetFluxHistogram()) {
GetFluxHistogram()->Write();
}
if (drawEvents and GetEventHistogram()) {
GetEventHistogram()->Write();
}
if (fIsMask and drawMask) {
fMaskHist->Write((fName + "_MSK").c_str()); //< save mask
}
if (drawMap) fMapHist->Write((fName + "_MAP").c_str()); //< save map
// // Save neut stack
// if (drawModes) {
// THStack combo_fMCHist_PDG = PlotUtils::GetNeutModeStack(
// (fName + "_MC_PDG").c_str(), (TH1**)fMCHist_PDG, 0);
// combo_fMCHist_PDG.Write();
// }
// Save Matrix plots
if (drawMatrix and fFullCovar and covar and fDecomp) {
TH2D cov = TH2D((*fFullCovar));
cov.SetNameTitle((fName + "_cov").c_str(),
(fName + "_cov;Bins; Bins;").c_str());
cov.Write();
TH2D covinv = TH2D((*this->covar));
covinv.SetNameTitle((fName + "_covinv").c_str(),
(fName + "_cov;Bins; Bins;").c_str());
covinv.Write();
TH2D covdec = TH2D((*fDecomp));
covdec.SetNameTitle((fName + "_covdec").c_str(),
(fName + "_cov;Bins; Bins;").c_str());
covdec.Write();
}
// Save ratio plots if required
if (drawRatio) {
// Needed for error bars
for (int i = 0; i < fMCHist->GetNbinsX() * fMCHist->GetNbinsY(); i++)
fMCHist->SetBinError(i + 1, 0.0);
fDataHist->GetSumw2();
fMCHist->GetSumw2();
// Create Ratio Histograms
TH2D* dataRatio = (TH2D*)fDataHist->Clone((fName + "_data_RATIO").c_str());
TH2D* mcRatio = (TH2D*)fMCHist->Clone((fName + "_MC_RATIO").c_str());
mcRatio->Divide(fMCHist);
dataRatio->Divide(fMCHist);
// Cancel bin errors on MC
for (int i = 0; i < mcRatio->GetNbinsX() * mcRatio->GetNbinsY(); i++) {
mcRatio->SetBinError(
i + 1, fMCHist->GetBinError(i + 1) / fMCHist->GetBinContent(i + 1));
}
mcRatio->SetMinimum(0);
mcRatio->SetMaximum(2);
dataRatio->SetMinimum(0);
dataRatio->SetMaximum(2);
mcRatio->Write();
dataRatio->Write();
delete mcRatio;
delete dataRatio;
}
// Save Shape Plots if required
if (drawShape) {
// Create Shape Histogram
TH2D* mcShape = (TH2D*)fMCHist->Clone((fName + "_MC_SHAPE").c_str());
double shapeScale = 1.0;
if (fIsRawEvents) {
shapeScale = fDataHist->Integral() / fMCHist->Integral();
} else {
shapeScale = fDataHist->Integral("width") / fMCHist->Integral("width");
}
mcShape->Scale(shapeScale);
mcShape->SetLineWidth(3);
mcShape->SetLineStyle(7); // dashes
mcShape->Write();
// Save shape ratios
if (drawRatio) {
// Needed for error bars
mcShape->GetSumw2();
// Create shape ratio histograms
TH2D* mcShapeRatio =
(TH2D*)mcShape->Clone((fName + "_MC_SHAPE_RATIO").c_str());
TH2D* dataShapeRatio =
(TH2D*)fDataHist->Clone((fName + "_data_SHAPE_RATIO").c_str());
// Divide the histograms
mcShapeRatio->Divide(mcShape);
dataShapeRatio->Divide(mcShape);
// Colour the shape ratio plots
mcShapeRatio->SetLineWidth(3);
mcShapeRatio->SetLineStyle(7); // dashes
mcShapeRatio->Write();
dataShapeRatio->Write();
delete mcShapeRatio;
delete dataShapeRatio;
}
delete mcShape;
}
// Save residual calculations of what contributed to the chi2 values.
if (residual) {
}
if (fIsProjFitX or fIsProjFitY or drawProj) {
// If not already made, make the projections
if (!fMCHist_X) {
PlotUtils::MatchEmptyBins(fDataHist, fMCHist);
fMCHist_X = PlotUtils::GetProjectionX(fMCHist, fMaskHist);
fMCHist_Y = PlotUtils::GetProjectionY(fMCHist, fMaskHist);
fDataHist_X = PlotUtils::GetProjectionX(fDataHist, fMaskHist);
fDataHist_Y = PlotUtils::GetProjectionY(fDataHist, fMaskHist);
double chi2X = StatUtils::GetChi2FromDiag(fDataHist_X, fMCHist_X);
double chi2Y = StatUtils::GetChi2FromDiag(fDataHist_Y, fMCHist_Y);
fMCHist_X->SetTitle(Form("%f", chi2X));
fMCHist_Y->SetTitle(Form("%f", chi2Y));
}
// Save the histograms
fDataHist_X->Write();
fMCHist_X->Write();
fDataHist_Y->Write();
fMCHist_Y->Write();
}
if (drawSliceCanvYMC or true) {
TCanvas* c1 = new TCanvas((fName + "_MC_CANV_Y").c_str(),
(fName + "_MC_CANV_Y").c_str(), 800, 600);
c1->Divide(int(sqrt(fDataHist->GetNbinsY() + 1)),
int(sqrt(fDataHist->GetNbinsY() + 1)));
TH2D* mcShape = (TH2D*)fMCHist->Clone((fName + "_MC_SHAPE").c_str());
double shapeScale =
fDataHist->Integral("width") / fMCHist->Integral("width");
mcShape->Scale(shapeScale);
mcShape->SetLineStyle(7);
c1->cd(1);
TLegend* leg = new TLegend(0.6, 0.6, 0.9, 0.9);
leg->AddEntry(fDataHist, (fName + " Data").c_str(), "ep");
leg->AddEntry(fMCHist, (fName + " MC").c_str(), "l");
leg->AddEntry(mcShape, (fName + " Shape").c_str(), "l");
leg->Draw("SAME");
/*
// Make Y slices
for (int i = 0; i < fDataHist->GetNbinY(); i++){
c1->cd(i+2);
TH1D* fDataHist_SliceY = PlotUtils::GetSliceY(fDataHist, i);
fDataHist_SliceY->Draw("E1");
TH1D* fMCHist_SliceY = PlotUtils::GetSliceY(fMCHist, i);
fMCHist_SliceY->Draw("SAME HIST C");
TH1D* mcShape_SliceY = PlotUtils::GetSliceY(mcShape, i);
mcShape_SliceY->Draw("SAME HIST C");
}
*/
c1->Write();
}
if (drawWeighted) {
fMCWeighted->Write();
}
// Returning
LOG(SAM) << "Written Histograms: " << fName << std::endl;
return;
// Returning
LOG(SAM) << "Written Histograms: " << fName << std::endl;
return;
}
/*
Setup Functions
*/
//********************************************************************
void Measurement2D::SetupMeasurement(std::string inputfile, std::string type,
FitWeight* rw, std::string fkdt) {
//********************************************************************
// Check if name contains Evt, indicating that it is a raw number of events
// measurements and should thus be treated as once
fIsRawEvents = false;
if ((fName.find("Evt") != std::string::npos) && fIsRawEvents == false) {
fIsRawEvents = true;
LOG(SAM) << "Found event rate measurement but fIsRawEvents == false!"
<< std::endl;
LOG(SAM) << "Overriding this and setting fIsRawEvents == true!"
<< std::endl;
}
fIsEnu = false;
if ((fName.find("XSec") != std::string::npos) &&
(fName.find("Enu") != std::string::npos)) {
fIsEnu = true;
LOG(SAM) << "::" << fName << "::" << std::endl;
LOG(SAM) << "Found XSec Enu measurement, applying flux integrated scaling, "
"not flux averaged!"
<< std::endl;
if (FitPar::Config().GetParB("EventManager")) {
ERR(FTL) << "Enu Measurements do not yet work with the Event Manager!"
<< std::endl;
ERR(FTL) << "If you want decent flux unfolded results please run in "
"series mode (-q EventManager=0)"
<< std::endl;
sleep(2);
}
}
if (fIsEnu && fIsRawEvents) {
LOG(SAM) << "Found 1D Enu XSec distribution AND fIsRawEvents, is this "
"really correct?!"
<< std::endl;
LOG(SAM) << "Check experiment constructor for " << fName
<< " and correct this!" << std::endl;
LOG(SAM) << "I live in " << __FILE__ << ":" << __LINE__ << std::endl;
exit(-1);
}
// Reset everything to NULL
fRW = rw;
// Setting up 2D Inputs
this->SetupInputs(inputfile);
// Set Default Options
SetFitOptions(fDefaultTypes);
// Set Passed Options
SetFitOptions(type);
}
//********************************************************************
void Measurement2D::SetupDefaultHist() {
//********************************************************************
// Setup fMCHist
fMCHist = (TH2D*)fDataHist->Clone();
fMCHist->SetNameTitle((fName + "_MC").c_str(),
(fName + "_MC" + fPlotTitles).c_str());
// Setup fMCFine
Int_t nBinsX = fMCHist->GetNbinsX();
Int_t nBinsY = fMCHist->GetNbinsY();
fMCFine = new TH2D((fName + "_MC_FINE").c_str(),
(fName + "_MC_FINE" + fPlotTitles).c_str(), nBinsX * 3,
fMCHist->GetXaxis()->GetBinLowEdge(1),
fMCHist->GetXaxis()->GetBinLowEdge(nBinsX + 1), nBinsY * 3,
fMCHist->GetYaxis()->GetBinLowEdge(1),
fMCHist->GetYaxis()->GetBinLowEdge(nBinsY + 1));
// Setup MC Stat
fMCStat = (TH2D*)fMCHist->Clone();
fMCStat->Reset();
// Setup the NEUT Mode Array
//PlotUtils::CreateNeutModeArray(fMCHist, (TH1**)fMCHist_PDG);
// Setup bin masks using sample name
if (fIsMask) {
std::string maskloc = FitPar::Config().GetParDIR(fName + ".mask");
if (maskloc.empty()) {
maskloc = FitPar::GetDataBase() + "/masks/" + fName + ".mask";
}
SetBinMask(maskloc);
}
return;
}
//********************************************************************
void Measurement2D::SetDataValues(std::string dataFile, std::string TH2Dname) {
//********************************************************************
if (dataFile.find(".root") == std::string::npos) {
ERR(FTL) << "Error! " << dataFile << " is not a .root file" << std::endl;
ERR(FTL) << "Currently only .root file reading is supported (MiniBooNE "
"CC1pi+ 2D), but implementing .txt should be dirt easy"
<< std::endl;
ERR(FTL) << "See me at " << __FILE__ << ":" << __LINE__ << std::endl;
exit(-1);
} else {
TFile* inFile = new TFile(dataFile.c_str(), "READ");
fDataHist = (TH2D*)(inFile->Get(TH2Dname.c_str())->Clone());
fDataHist->SetDirectory(0);
fDataHist->SetNameTitle((fName + "_data").c_str(),
(fName + "_MC" + fPlotTitles).c_str());
delete inFile;
}
return;
}
//********************************************************************
void Measurement2D::SetDataValues(std::string dataFile, double dataNorm,
std::string errorFile, double errorNorm) {
//********************************************************************
// Make a counter to track the line number
int yBin = 0;
std::string line;
- std::ifstream data(dataFile.c_str(), ifstream::in);
+ std::ifstream data(dataFile.c_str(), std::ifstream::in);
fDataHist = new TH2D((fName + "_data").c_str(), (fName + fPlotTitles).c_str(),
fNDataPointsX - 1, fXBins, fNDataPointsY - 1, fYBins);
if (data.is_open())
LOG(SAM) << "Reading data from: " << dataFile.c_str() << std::endl;
while (std::getline(data >> std::ws, line, '\n')) {
int xBin = 0;
// Loop over entries and insert them into the histogram
std::vector entries = GeneralUtils::ParseToDbl(line, " ");
for (std::vector::iterator iter = entries.begin();
iter != entries.end(); iter++) {
fDataHist->SetBinContent(xBin + 1, yBin + 1, (*iter) * dataNorm);
xBin++;
}
yBin++;
}
yBin = 0;
- std::ifstream error(errorFile.c_str(), ifstream::in);
+ std::ifstream error(errorFile.c_str(), std::ifstream::in);
if (error.is_open())
LOG(SAM) << "Reading errors from: " << errorFile.c_str() << std::endl;
while (std::getline(error >> std::ws, line, '\n')) {
int xBin = 0;
// Loop over entries and insert them into the histogram
std::vector entries = GeneralUtils::ParseToDbl(line, " ");
for (std::vector::iterator iter = entries.begin();
iter != entries.end(); iter++) {
fDataHist->SetBinError(xBin + 1, yBin + 1, (*iter) * errorNorm);
xBin++;
}
yBin++;
}
return;
};
//********************************************************************
void Measurement2D::SetDataValuesFromText(std::string dataFile,
double dataNorm) {
//********************************************************************
fDataHist = new TH2D((fName + "_data").c_str(), (fName + fPlotTitles).c_str(),
fNDataPointsX - 1, fXBins, fNDataPointsY - 1, fYBins);
LOG(SAM) << "Reading data from: " << dataFile << std::endl;
PlotUtils::Set2DHistFromText(dataFile, fDataHist, dataNorm, true);
return;
};
//********************************************************************
void Measurement2D::SetCovarMatrix(std::string covarFile) {
//********************************************************************
// Used to read a covariance matrix from a root file
TFile* tempFile = new TFile(covarFile.c_str(), "READ");
// Make plots that we want
TH2D* covarPlot = new TH2D();
// TH2D* decmpPlot = new TH2D();
TH2D* covarInvPlot = new TH2D();
TH2D* fFullCovarPlot = new TH2D();
// Get covariance options for fake data studies
std::string covName = "";
std::string covOption = FitPar::Config().GetParS("throw_covariance");
// Which matrix to get?
if (fIsShape || fIsFree) covName = "shp_";
if (fIsDiag)
covName += "diag";
else
covName += "full";
covarPlot = (TH2D*)tempFile->Get((covName + "cov").c_str());
covarInvPlot = (TH2D*)tempFile->Get((covName + "covinv").c_str());
// Throw either the sub matrix or the full matrix
if (!covOption.compare("SUB"))
fFullCovarPlot = (TH2D*)tempFile->Get((covName + "cov").c_str());
else if (!covOption.compare("FULL"))
fFullCovarPlot = (TH2D*)tempFile->Get("fullcov");
else
ERR(WRN) << " Incorrect thrown_covariance option in parameters."
<< std::endl;
// Bin masking?
int dim = int(fDataHist->GetNbinsX()); //-this->masked->Integral());
int covdim = int(fDataHist->GetNbinsX());
// Make new covars
this->covar = new TMatrixDSym(dim);
fFullCovar = new TMatrixDSym(dim);
fDecomp = new TMatrixDSym(dim);
// Full covariance values
int row, column = 0;
row = 0;
column = 0;
for (Int_t i = 0; i < covdim; i++) {
// masking can be dodgy
// if (this->masked->GetBinContent(i+1) > 0) continue;
for (Int_t j = 0; j < covdim; j++) {
// if (this->masked->GetBinContent(j+1) > 0) continue;
(*this->covar)(row, column) = covarPlot->GetBinContent(i + 1, j + 1);
(*fFullCovar)(row, column) = fFullCovarPlot->GetBinContent(i + 1, j + 1);
column++;
}
column = 0;
row++;
}
// Set bin errors on data
if (!fIsDiag) {
for (Int_t i = 0; i < fDataHist->GetNbinsX(); i++) {
fDataHist->SetBinError(
i + 1, sqrt((covarPlot->GetBinContent(i + 1, i + 1))) * 1E-38);
}
}
TDecompSVD LU = TDecompSVD(*this->covar);
this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), "");
tempFile->Close();
delete tempFile;
return;
};
//********************************************************************
void Measurement2D::SetCovarMatrixFromText(std::string covarFile, int dim) {
//********************************************************************
// Make a counter to track the line number
int row = 0;
std::string line;
- std::ifstream covar(covarFile.c_str(), ifstream::in);
+ std::ifstream covar(covarFile.c_str(), std::ifstream::in);
this->covar = new TMatrixDSym(dim);
fFullCovar = new TMatrixDSym(dim);
if (covar.is_open())
LOG(SAM) << "Reading covariance matrix from file: " << covarFile
<< std::endl;
while (std::getline(covar >> std::ws, line, '\n')) {
int column = 0;
// Loop over entries and insert them into matrix
// Multiply by the errors to get the covariance, rather than the correlation
// matrix
std::vector entries = GeneralUtils::ParseToDbl(line, " ");
for (std::vector::iterator iter = entries.begin();
iter != entries.end(); iter++) {
double val = (*iter) * fDataHist->GetBinError(row + 1) * 1E38 *
fDataHist->GetBinError(column + 1) * 1E38;
(*this->covar)(row, column) = val;
(*fFullCovar)(row, column) = val;
column++;
}
row++;
}
// Robust matrix inversion method
TDecompSVD LU = TDecompSVD(*this->covar);
this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), "");
return;
};
//********************************************************************
void Measurement2D::SetCovarMatrixFromChol(std::string covarFile, int dim) {
//********************************************************************
// Make a counter to track the line number
int row = 0;
std::string line;
- std::ifstream covarread(covarFile.c_str(), ifstream::in);
+ std::ifstream covarread(covarFile.c_str(), std::ifstream::in);
TMatrixD* newcov = new TMatrixD(dim, dim);
if (covarread.is_open())
LOG(SAM) << "Reading covariance matrix from file: " << covarFile
<< std::endl;
while (std::getline(covarread >> std::ws, line, '\n')) {
int column = 0;
// Loop over entries and insert them into matrix
// Multiply by the errors to get the covariance, rather than the correlation
// matrix
std::vector entries = GeneralUtils::ParseToDbl(line, " ");
for (std::vector::iterator iter = entries.begin();
iter != entries.end(); iter++) {
(*newcov)(row, column) = *iter;
column++;
}
row++;
}
covarread.close();
// Form full covariance
TMatrixD* trans = (TMatrixD*)(newcov)->Clone();
trans->T();
(*trans) *= (*newcov);
fFullCovar = new TMatrixDSym(dim, trans->GetMatrixArray(), "");
delete newcov;
delete trans;
// Robust matrix inversion method
TDecompChol LU = TDecompChol(*this->fFullCovar);
this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), "");
return;
};
// //********************************************************************
// void Measurement2D::SetMapValuesFromText(std::string dataFile) {
// //********************************************************************
// fMapHist = new TH2I((fName + "_map").c_str(), (fName + fPlotTitles).c_str(),
// fNDataPointsX - 1, fXBins, fNDataPointsY - 1, fYBins);
// LOG(SAM) << "Reading map from: " << dataFile << std::endl;
// PlotUtils::Set2DHistFromText(dataFile, fMapHist, 1.0);
// return;
// };
diff --git a/src/Reweight/WeightUtils.cxx b/src/Reweight/WeightUtils.cxx
index a17eb22..7204c64 100644
--- a/src/Reweight/WeightUtils.cxx
+++ b/src/Reweight/WeightUtils.cxx
@@ -1,615 +1,615 @@
#include "WeightUtils.h"
#include "FitLogger.h"
#ifdef __T2KREW_ENABLED__
#include "T2KGenieReWeight.h"
#include "T2KNIWGReWeight.h"
#include "T2KNIWGUtils.h"
#include "T2KNeutReWeight.h"
#include "T2KNeutUtils.h"
#include "T2KReWeight.h"
using namespace t2krew;
#endif
#ifdef __NIWG_ENABLED__
#include "NIWGReWeight.h"
#include "NIWGReWeight1piAngle.h"
#include "NIWGReWeight2010a.h"
#include "NIWGReWeight2012a.h"
#include "NIWGReWeight2014a.h"
#include "NIWGReWeightDeltaMass.h"
#include "NIWGReWeightEffectiveRPA.h"
#include "NIWGReWeightHadronMultSwitch.h"
#include "NIWGReWeightMEC.h"
#include "NIWGReWeightPiMult.h"
#include "NIWGReWeightProtonFSIbug.h"
#include "NIWGReWeightRPA.h"
#include "NIWGReWeightSpectralFunc.h"
#include "NIWGReWeightSplineEnu.h"
#include "NIWGSyst.h"
#include "NIWGSystUncertainty.h"
#endif
#ifdef __NEUT_ENABLED__
#include "NReWeight.h"
#include "NReWeightCasc.h"
#include "NReWeightNuXSecCCQE.h"
#include "NReWeightNuXSecCCRES.h"
#include "NReWeightNuXSecCOH.h"
#include "NReWeightNuXSecDIS.h"
#include "NReWeightNuXSecNC.h"
#include "NReWeightNuXSecNCEL.h"
#include "NReWeightNuXSecNCRES.h"
#include "NReWeightNuXSecRES.h"
#include "NReWeightNuclPiless.h"
#include "NSyst.h"
#include "NSystUncertainty.h"
#include "neutpart.h"
#include "neutvect.h"
#endif
#ifdef __NUWRO_ENABLED__
#include "event1.h"
#endif
#ifdef __NUWRO_REWEIGHT_ENABLED__
#include "NuwroReWeight.h"
#include "NuwroReWeight_FlagNorm.h"
#include "NuwroReWeight_QEL.h"
#include "NuwroReWeight_SPP.h"
#include "NuwroSyst.h"
#include "NuwroSystUncertainty.h"
#endif
#ifdef __GENIE_ENABLED__
#include "EVGCore/EventRecord.h"
#include "EVGCore/EventRecord.h"
#include "GHEP/GHepRecord.h"
#include "GSyst.h"
#include "GSystUncertainty.h"
#include "Ntuple/NtpMCEventRecord.h"
#include "ReWeight/GReWeight.h"
#include "ReWeight/GReWeightAGKY.h"
#include "ReWeight/GReWeightDISNuclMod.h"
#include "ReWeight/GReWeightFGM.h"
#include "ReWeight/GReWeightFZone.h"
#include "ReWeight/GReWeightINuke.h"
#include "ReWeight/GReWeightNonResonanceBkg.h"
#include "ReWeight/GReWeightNuXSecCCQE.h"
#include "ReWeight/GReWeightNuXSecCCQEvec.h"
#include "ReWeight/GReWeightNuXSecCCRES.h"
#include "ReWeight/GReWeightNuXSecCOH.h"
#include "ReWeight/GReWeightNuXSecDIS.h"
#include "ReWeight/GReWeightNuXSecNC.h"
#include "ReWeight/GReWeightNuXSecNCEL.h"
#include "ReWeight/GReWeightNuXSecNCRES.h"
#include "ReWeight/GReWeightResonanceDecay.h"
using namespace genie;
using namespace genie::rew;
#endif
#include "GlobalDialList.h"
#include "ModeNormEngine.h"
#include "NUISANCESyst.h"
#include "OscWeightEngine.h"
//********************************************************************
TF1 FitBase::GetRWConvFunction(std::string const &type,
std::string const &name) {
//********************************************************************
std::string dialfunc = "x";
std::string parType = type;
double low = -10000.0;
double high = 10000.0;
if (parType.find("parameter") == std::string::npos) parType += "_parameter";
std::string line;
- ifstream card(
+ std::ifstream card(
(GeneralUtils::GetTopLevelDir() + "/parameters/dial_conversion.card")
.c_str(),
- ifstream::in);
+ std::ifstream::in);
while (std::getline(card >> std::ws, line, '\n')) {
std::vector inputlist = GeneralUtils::ParseToStr(line, " ");
// Check the line length
if (inputlist.size() < 4) continue;
// Check whether this is a comment
if (inputlist[0].c_str()[0] == '#') continue;
// Check whether this is the correct parameter type
if (inputlist[0].compare(parType) != 0) continue;
// Check the parameter name
if (inputlist[1].compare(name) != 0) continue;
// inputlist[2] should be the units... ignore for now
dialfunc = inputlist[3];
// High and low are optional, check whether they exist
if (inputlist.size() > 4) low = GeneralUtils::StrToDbl(inputlist[4]);
if (inputlist.size() > 5) high = GeneralUtils::StrToDbl(inputlist[5]);
}
TF1 convfunc = TF1((name + "_convfunc").c_str(), dialfunc.c_str(), low, high);
return convfunc;
}
//********************************************************************
std::string FitBase::GetRWUnits(std::string const &type,
std::string const &name) {
//********************************************************************
std::string unit = "sig.";
std::string parType = type;
if (parType.find("parameter") == std::string::npos) {
parType += "_parameter";
}
std::string line;
std::ifstream card(
(GeneralUtils::GetTopLevelDir() + "/parameters/dial_conversion.card")
.c_str(),
- ifstream::in);
+ std::ifstream::in);
while (std::getline(card >> std::ws, line, '\n')) {
std::vector inputlist = GeneralUtils::ParseToStr(line, " ");
// Check the line length
if (inputlist.size() < 3) continue;
// Check whether this is a comment
if (inputlist[0].c_str()[0] == '#') continue;
// Check whether this is the correct parameter type
if (inputlist[0].compare(parType) != 0) continue;
// Check the parameter name
if (inputlist[1].compare(name) != 0) continue;
unit = inputlist[2];
break;
}
return unit;
}
//********************************************************************
double FitBase::RWAbsToSigma(std::string const &type, std::string const &name,
double val) {
//********************************************************************
TF1 f1 = GetRWConvFunction(type, name);
double conv_val = f1.GetX(val);
if (fabs(conv_val) < 1E-10) conv_val = 0.0;
std::cout << "AbsToSigma(" << name << ") = " << val << " -> " << conv_val
<< std::endl;
return conv_val;
}
//********************************************************************
double FitBase::RWSigmaToAbs(std::string const &type, std::string const &name,
double val) {
//********************************************************************
TF1 f1 = GetRWConvFunction(type, name);
double conv_val = f1.Eval(val);
return conv_val;
}
//********************************************************************
double FitBase::RWFracToSigma(std::string const &type, std::string const &name,
double val) {
//********************************************************************
TF1 f1 = GetRWConvFunction(type, name);
double conv_val = f1.GetX((val * f1.Eval(0.0)));
if (fabs(conv_val) < 1E-10) conv_val = 0.0;
return conv_val;
}
//********************************************************************
double FitBase::RWSigmaToFrac(std::string const &type, std::string const &name,
double val) {
//********************************************************************
TF1 f1 = GetRWConvFunction(type, name);
double conv_val = f1.Eval(val) / f1.Eval(0.0);
return conv_val;
}
int FitBase::ConvDialType(std::string const &type) {
if (!type.compare("neut_parameter"))
return kNEUT;
else if (!type.compare("niwg_parameter"))
return kNIWG;
else if (!type.compare("nuwro_parameter"))
return kNUWRO;
else if (!type.compare("t2k_parameter"))
return kT2K;
else if (!type.compare("genie_parameter"))
return kGENIE;
else if (!type.compare("custom_parameter"))
return kCUSTOM;
else if (!type.compare("norm_parameter"))
return kNORM;
else if (!type.compare("likeweight_parameter"))
return kLIKEWEIGHT;
else if (!type.compare("spline_parameter"))
return kSPLINEPARAMETER;
else if (!type.compare("osc_parameter"))
return kOSCILLATION;
else if (!type.compare("modenorm_parameter"))
return kMODENORM;
else
return kUNKNOWN;
}
std::string FitBase::ConvDialType(int type) {
switch (type) {
case kNEUT: {
return "neut_parameter";
}
case kNIWG: {
return "niwg_parameter";
}
case kNUWRO: {
return "nuwro_parameter";
}
case kT2K: {
return "t2k_parameter";
}
case kGENIE: {
return "genie_parameter";
}
case kNORM: {
return "norm_parameter";
}
case kCUSTOM: {
return "custom_parameter";
}
case kLIKEWEIGHT: {
return "likeweight_parameter";
}
case kSPLINEPARAMETER: {
return "spline_parameter";
}
case kOSCILLATION: {
return "osc_parameter";
}
case kMODENORM: {
return "modenorm_parameter";
}
default:
return "unknown_parameter";
}
}
int FitBase::GetDialEnum(std::string const &type, std::string const &name) {
return FitBase::GetDialEnum(FitBase::ConvDialType(type), name);
}
int FitBase::GetDialEnum(int type, std::string const &name) {
int offset = type * 1000;
int this_enum = Reweight::kNoDialFound; // Not Found
std::cout << "Getting dial enum " << type << " " << name << std::endl;
// Select Types
switch (type) {
// NEUT DIAL TYPE
case kNEUT: {
#ifdef __NEUT_ENABLED__
int neut_enum = (int)neut::rew::NSyst::FromString(name);
if (neut_enum != 0) {
this_enum = neut_enum + offset;
}
#else
this_enum = Reweight::kNoTypeFound; // Not enabled
#endif
break;
}
// NIWG DIAL TYPE
case kNIWG: {
#ifdef __NIWG_ENABLED__
int niwg_enum = (int)niwg::rew::NIWGSyst::FromString(name);
if (niwg_enum != 0) {
this_enum = niwg_enum + offset;
}
#else
this_enum = Reweight::kNoTypeFound;
#endif
break;
}
// NUWRO DIAL TYPE
case kNUWRO: {
#ifdef __NUWRO_REWEIGHT_ENABLED__
int nuwro_enum = (int)nuwro::rew::NuwroSyst::FromString(name);
if (nuwro_enum > 0) {
this_enum = nuwro_enum + offset;
}
#else
this_enum = Reweight::kNoTypeFound;
#endif
}
// GENIE DIAL TYPE
case kGENIE: {
#ifdef __GENIE_ENABLED__
int genie_enum = (int)genie::rew::GSyst::FromString(name);
if (genie_enum > 0) {
this_enum = genie_enum + offset;
}
#else
this_enum = Reweight::kNoTypeFound;
#endif
break;
}
case kCUSTOM: {
int custom_enum = 0; // PLACEHOLDER
this_enum = custom_enum + offset;
break;
}
// T2K DIAL TYPE
case kT2K: {
#ifdef __T2KREW_ENABLED__
int t2k_enum = (int)t2krew::T2KSyst::FromString(name);
if (t2k_enum > 0) {
this_enum = t2k_enum + offset;
}
#else
this_enum = Reweight::kNoTypeFound;
#endif
break;
}
case kNORM: {
if (gNormEnums.find(name) == gNormEnums.end()) {
gNormEnums[name] = gNormEnums.size() + 1 + offset;
}
this_enum = gNormEnums[name];
break;
}
case kLIKEWEIGHT: {
if (gLikeWeightEnums.find(name) == gLikeWeightEnums.end()) {
gLikeWeightEnums[name] = gLikeWeightEnums.size() + 1 + offset;
}
this_enum = gLikeWeightEnums[name];
break;
}
case kSPLINEPARAMETER: {
if (gSplineParameterEnums.find(name) == gSplineParameterEnums.end()) {
gSplineParameterEnums[name] = gSplineParameterEnums.size() + 1 + offset;
}
this_enum = gSplineParameterEnums[name];
break;
}
case kOSCILLATION: {
#ifdef __PROB3PP_ENABLED__
int oscEnum = OscWeightEngine::SystEnumFromString(name);
if (oscEnum != 0) {
this_enum = oscEnum + offset;
}
#else
this_enum = Reweight::kNoTypeFound; // Not enabled
#endif
}
case kMODENORM: {
size_t us_pos = name.find_first_of('_');
std::string numstr = name.substr(us_pos + 1);
int mode_num = std::atoi(numstr.c_str());
LOG(FTL) << "Getting mode num " << mode_num << std::endl;
if (!mode_num) {
ERR(FTL) << "Attempting to parse dial name: \"" << name
<< "\" as a mode norm dial but failed." << std::endl;
throw;
}
this_enum = 60 + mode_num + offset;
break;
}
}
// If Not Enabled
if (this_enum == Reweight::kNoTypeFound) {
ERR(FTL) << "RW Engine not supported for " << FitBase::ConvDialType(type)
<< std::endl;
ERR(FTL) << "Check dial " << name << std::endl;
}
// If Not Found
if (this_enum == Reweight::kNoDialFound) {
ERR(FTL) << "Dial " << name << " not found." << std::endl;
}
return this_enum;
}
int Reweight::ConvDialType(std::string const &type) {
return FitBase::ConvDialType(type);
}
std::string Reweight::ConvDialType(int type) {
return FitBase::ConvDialType(type);
}
int Reweight::GetDialType(int type) {
int t = (type / 1000);
return t > kMODENORM ? Reweight::kNoDialFound : t;
}
int Reweight::RemoveDialType(int type) { return (type % 1000); }
int Reweight::NEUTEnumFromName(std::string const &name) {
#ifdef __NEUT_ENABLED__
int neutenum = (int)neut::rew::NSyst::FromString(name);
return (neutenum > 0) ? neutenum : Reweight::kNoDialFound;
#else
return Reweight::kGeneratorNotBuilt;
#endif
}
int Reweight::NIWGEnumFromName(std::string const &name) {
#ifdef __NIWG_ENABLED__
int niwgenum = (int)niwg::rew::NIWGSyst::FromString(name);
return (niwgenum != 0) ? niwgenum : Reweight::kNoDialFound;
#else
return Reweight::kGeneratorNotBuilt;
#endif
}
int Reweight::NUWROEnumFromName(std::string const &name) {
#ifdef __NUWRO_REWEIGHT_ENABLED__
int nuwroenum = (int)nuwro::rew::NuwroSyst::FromString(name);
return (nuwroenum > 0) ? nuwroenum : Reweight::kNoDialFound;
#else
return Reweight::kGeneratorNotBuilt;
#endif
}
int Reweight::GENIEEnumFromName(std::string const &name) {
#ifdef __GENIE_ENABLED__
int genieenum = (int)genie::rew::GSyst::FromString(name);
return (genieenum > 0) ? genieenum : Reweight::kNoDialFound;
#else
return Reweight::kGeneratorNotBuilt;
#endif
}
int Reweight::T2KEnumFromName(std::string const &name) {
#ifdef __T2KREW_ENABLED__
int t2kenum = (int)t2krew::T2KSyst::FromString(name);
return (t2kenum > 0) ? t2kenum : Reweight::kNoDialFound;
#else
return Reweight::kGeneratorNotBuilt;
#endif
}
int Reweight::OscillationEnumFromName(std::string const &name) {
#ifdef __PROB3PP_ENABLED__
int oscEnum = OscWeightEngine::SystEnumFromString(name);
return (oscEnum > 0) ? oscEnum : Reweight::kNoDialFound;
#else
return Reweight::kGeneratorNotBuilt;
#endif
}
int Reweight::NUISANCEEnumFromName(std::string const &name, int type) {
int nuisenum = Reweight::DialList().EnumFromNameAndType(name, type);
return nuisenum;
}
int Reweight::CustomEnumFromName(std::string const &name) {
int custenum = Reweight::ConvertNUISANCEDial(name);
return custenum;
}
int Reweight::ConvDial(std::string const &name, std::string const &type,
bool exceptions) {
return Reweight::ConvDial(name, Reweight::ConvDialType(type), exceptions);
}
int Reweight::ConvDial(std::string const &fullname, int type, bool exceptions) {
std::string name =
GeneralUtils::ParseToStr(fullname, ",")[0]; // Only use first dial given
// Produce offset seperating each type.
int offset = type * 1000;
int genenum = Reweight::kNoDialFound;
switch (type) {
case kNEUT:
genenum = NEUTEnumFromName(name);
break;
case kNIWG:
genenum = NIWGEnumFromName(name);
break;
case kNUWRO:
genenum = NUWROEnumFromName(name);
break;
case kGENIE:
genenum = GENIEEnumFromName(name);
break;
case kT2K:
genenum = T2KEnumFromName(name);
break;
case kCUSTOM:
genenum = CustomEnumFromName(name);
break;
case kNORM:
case kLIKEWEIGHT:
case kSPLINEPARAMETER:
case kNEWSPLINE:
genenum = NUISANCEEnumFromName(name, type);
break;
case kOSCILLATION:
genenum = OscillationEnumFromName(name);
break;
case kMODENORM:
genenum = ModeNormEngine::SystEnumFromString(name);
break;
default:
genenum = Reweight::kNoTypeFound;
break;
}
// Throw if required.
if (exceptions) {
// If Not Enabled
if (genenum == Reweight::kGeneratorNotBuilt) {
ERR(FTL) << "RW Engine not supported for " << FitBase::ConvDialType(type)
<< std::endl;
ERR(FTL) << "Check dial " << name << std::endl;
throw;
}
// If no type enabled
if (genenum == Reweight::kNoTypeFound) {
ERR(FTL) << "Type mismatch inside ConvDialEnum" << std::endl;
throw;
}
// If Not Found
if (genenum == Reweight::kNoDialFound) {
ERR(FTL) << "Dial " << name << " not found." << std::endl;
throw;
}
}
// Add offset if no issue
int nuisenum = genenum;
if ((genenum != Reweight::kGeneratorNotBuilt) &&
(genenum != Reweight::kNoTypeFound) &&
(genenum != Reweight::kNoDialFound)) {
nuisenum += offset;
}
// Now register dial
Reweight::DialList().RegisterDialEnum(name, type, nuisenum);
return nuisenum;
}
std::string Reweight::ConvDial(int nuisenum) {
for (size_t i = 0; i < Reweight::DialList().fAllDialEnums.size(); i++) {
if (Reweight::DialList().fAllDialEnums[i] == nuisenum) {
return Reweight::DialList().fAllDialNames[i];
}
}
LOG(FIT) << "Cannot find dial with enum = " << nuisenum << std::endl;
return "";
}
diff --git a/src/Routines/MinimizerRoutines.h b/src/Routines/MinimizerRoutines.h
index 0733601..86b2d1d 100755
--- a/src/Routines/MinimizerRoutines.h
+++ b/src/Routines/MinimizerRoutines.h
@@ -1,270 +1,276 @@
// 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 .
*******************************************************************************/
#ifndef MINIMIZER_ROUTINES_H
#define MINIMIZER_ROUTINES_H
/*!
* \addtogroup Minimizer
* @{
*/
#include "TH1.h"
#include "TF1.h"
#include "TMatrixD.h"
#include "TVectorD.h"
+
+#ifdef ROOT6_USE_FIT_FITTER_INTERFACE
+#include "Fit/Fitter.h"
+#else
#include "Minuit2/FCNBase.h"
#include "TFitterMinuit.h"
+#endif
+
#include "TSystem.h"
#include "TFile.h"
#include "TProfile.h"
#include
#include
#include
#include
#include
#include "FitEvent.h"
#include "JointFCN.h"
#include "MinimizerFCN.h"
#include "Math/Minimizer.h"
#include "Math/Factory.h"
#include "Math/Functor.h"
#include "FitLogger.h"
#include "ParserUtils.h"
enum minstate {
kErrorStatus = -1,
kGoodStatus,
kFitError,
kNoChange,
kFitFinished,
kFitUnfinished,
kStateChange,
};
//*************************************
//! Collects all possible fit routines into a single class to avoid repeated code
class MinimizerRoutines{
//*************************************
public:
/*
Constructor/Destructor
*/
//! Constructor reads in arguments given at the command line for the fit here.
MinimizerRoutines(int argc, char* argv[]);
-
+
//! Default destructor
~MinimizerRoutines();
//! Reset everything to default/NULL
void Init();
-
+
/*
Input Functions
*/
//! Splits the arguments ready for initial setup
void ParseArgs(int argc, char* argv[]);
//! Sorts out configuration and verbosity right at the very start.
//! Calls readCard to set everything else up.
void InitialSetup();
//! Loops through each line of the card file and passes it to other read functions
void ReadCard(std::string cardfile);
//! Check for parameter string in the line and assign the correct type.
//! Fills maps for each of the parameters
int ReadParameters(std::string parstring);
//! Reads in fake parameters and assigns them (Requires the parameter to be included as a normal parameter as well)
int ReadFakeDataPars(std::string parstring);
//! Read in the samples so we can set up the free normalisation dials if required
int ReadSamples(std::string sampleString);
void SetupMinimizerFromXML();
/*
Setup Functions
*/
//! Setup the configuration given the arguments passed at the commandline and card file
void SetupConfig();
//! Setups up our custom RW engine with all the parameters passed in the card file
void SetupRWEngine();
//! Setups up the jointFCN.
void SetupFCN();
//! Sets up the minimizerObj for ROOT. there are cases where this is called repeatedly, e.g. If you are using a brute force scan before using Migrad.
void SetupFitter(std::string routine);
//! Set the current data histograms in each sample to the fake data.
void SetFakeData();
//! Setup the covariances with the correct dimensions. At the start this is either uncorrelated or merged given all the input covariances.
//! At the end of the fit this produces the blank covariances which can then be filled by the minimizerObj with best fit covariances.
void SetupCovariance();
/*
Fitting Functions
*/
//! Main function to actually start iterating over the different required fit routines
void Run();
//! Given a new map change the values that the RW engine is currently set to
void UpdateRWEngine(std::map& updateVals);
//! Given a single routine (see tutorial for options) run that fit routine now.
int RunFitRoutine(std::string routine);
//! Get the current state of minimizerObj and fill it into currentVals and currentNorms
void GetMinimizerState();
//! Print current value
void PrintState();
-
+
//! Performs a fit routine where the input.maxevents is set to a much lower value to try and move closer to the best fit minimum.
void LowStatRoutine(std::string routine);
//! Perform a chi2 scan in 1D around the current point
void Create1DScans();
//! Perform a chi2 scan in 2D around the current point
void Chi2Scan2D();
//! Currently a placeholder NEEDS UPDATING
void CreateContours();
//! If any currentVals are close to the limits set them to the limit and fix them
int FixAtLimit();
//! Throw the current covariance of dial values we have, and fill the thrownVals and thrownNorms maps.
//! If uniformly is true parameters will be thrown uniformly between their upper and lower limits.
void ThrowCovariance(bool uniformly);
//! Given the covariance we currently have generate error bands by throwing the covariance.
//! The FitPar config "error_uniform" defines whether to throw using the covariance or uniformly.
//! The FitPar config "error_throws" defines how many throws are needed.
//! Currently only supports TH1D plots.
void GenerateErrorBands();
/*
Write Functions
*/
//! Write plots and TTrees listing the minimizerObj result of the fit to file
void SaveMinimizerState();
//! Save the sample plots for current MC
//! dir if not empty forces plots to be saved in a subdirectory of outputfile
void SaveCurrentState(std::string subdir="");
//! Save starting predictions into a seperate folder
void SaveNominal();
//! Save predictions before the fit is ran into a seperate folder
void SavePrefit();
void SaveResults();
/*
MISC Functions
*/
//! Get previous fit status from a file
Int_t GetStatus();
/// Makes a histogram of likelihoods when throwing the data according to its statistics
void ThrowDataToys();
protected:
//! Our Custom ReWeight Object
FitWeight* rw;
std::string fOutputFile;
std::string fInputFile;
TFile* fInputRootFile;
TFile* fOutputRootFile;
//! Flag for whether the fit should be continued if an output file is already found.
bool fitContinue;
//! Minimizer Object for handling roots different minimizer methods
ROOT::Math::Minimizer* fMinimizer;
JointFCN* fSampleFCN;
MinimizerFCN* fMinimizerFCN;
ROOT::Math::Functor* fCallFunctor;
int nfreepars;
std::string fCardFile;
std::string fStrategy;
std::vector fRoutines;
std::string fAllowedRoutines;
-
+
std::string fFakeDataInput;
// Input Dial Vals
//! Vector of dial names
std::vector fParams;
std::map fStateVals;
std::map fStartVals;
std::map fCurVals;
std::map fErrorVals;
std::map fMinVals;
std::map fMaxVals;
std::map fStepVals;
std::map fTypeVals;
std::map fFixVals;
std::map fStartFixVals;
//! Vector of fake parameter names
std::map fFakeVals;
//! Map of thrown parameter names and values (After ThrowCovariance)
std::map fThrownVals;
TH2D* fCorrel;
TH2D* fDecomp;
TH2D* fCovar;
-
+
TH2D* fCorFree;
TH2D* fDecFree;
TH2D* fCovFree;
nuiskey fCompKey;
};
/*! @} */
#endif
diff --git a/src/Splines/SplineWriter.h b/src/Splines/SplineWriter.h
index ffc0655..87d7cc6 100644
--- a/src/Splines/SplineWriter.h
+++ b/src/Splines/SplineWriter.h
@@ -1,92 +1,98 @@
#ifndef SPLINEWRITER_H
#define SPLINEWRITER_H
#include "FitWeight.h"
#include "Spline.h"
#include "SplineUtils.h"
#ifdef __MINUIT2_ENABLED__
+
+#ifdef ROOT6_USE_FIT_FITTER_INTERFACE
+#include "Fit/Fitter.h"
+#else
#include "TFitterMinuit.h"
#endif
+#endif
+
class SplineFCN {
public:
SplineFCN(Spline* spl, std::vector > v, std::vector w) { fSpl = spl; fVal = v; fWeight = w; };
~SplineFCN() {};
double operator()(const double* x) const;
double DoEval(const double *x) const;
void SaveAs(std::string name, const float* fx);
void UpdateWeights(std::vector& w);
void SetCorrelated(bool state = true);
bool uncorrelated;
std::vector< std::vector > fVal;
std::vector< double > fWeight;
Spline* fSpl;
};
class SplineWriter : public SplineReader {
public:
SplineWriter(FitWeight* fw) {
fRW = fw;
fDrawSplines = FitPar::Config().GetParB("drawsplines");
};
~SplineWriter() {};
void SetupSplineSet();
void Write(std::string name);
void AddCoefficientsToTree(TTree* tree);
void FitSplinesForEvent(TCanvas* fitcanvas = NULL, bool saveplot = false);
void AddWeightsToTree(TTree* tr);
void ReadWeightsFromTree(TTree* tr);
void FitSplinesForEvent(double* weightvals, float* coeff);
void GetWeightsForEvent(FitEvent* event, double* weights);
void GetWeightsForEvent(FitEvent* event);
void ReconfigureSet(int iset);
double GetWeightForThisSet(FitEvent* event, int iset=-1);
void SetWeights(double* weights);
inline int GetNWeights(){return fParVect.size();};
inline int GetNPars(){ return fNCoEff;};
int fNCoEff;
// double* fCoEffStorer;
float* fCoEffStorer;
std::vector< std::vector > fParVect;
std::vector< int > fSetIndex;
double* fWeightList;
std::vector< std::vector > fValList;
int fCurrentSet;
FitWeight* fRW;
bool fDrawSplines;
std::vector fAllDrawnHists;
std::vector fAllDrawnGraphs;
#ifdef __MINUIT2_ENABLED__
std::map fSplineFCNs;
std::map fSplineFunctors;
std::map fSplineMinimizers;
#endif
// Spline* gSpline;
// Available Fitting Functions
void FitCoeff(Spline* spl, std::vector< std::vector >& v, std::vector& w, float* coeff, bool draw);
void FitCoeff1DGraph(Spline* spl, int n, double* x, double* y, float* coeff, bool draw);
void GetCoeff1DTSpline3(Spline* spl, int n, double* x, double* y, float* coeff, bool draw);
// void FitCoeff2DGraph(Spline* spl, std::vector< std::vector >& v, std::vector& w, float* coeff, bool draw);
void FitCoeffNDGraph(Spline* spl, std::vector< std::vector >& v, std::vector& w, float* coeff, bool draw);
void FitCoeff2DGraph(Spline* spl, int n, double* x, double* y, double* w, float* coeff, bool draw);
//double Func2DWrapper(double* x, double* p);
};
#endif
diff --git a/src/Utils/PlotUtils.cxx b/src/Utils/PlotUtils.cxx
index ec742e3..0f8bbcb 100644
--- a/src/Utils/PlotUtils.cxx
+++ b/src/Utils/PlotUtils.cxx
@@ -1,1001 +1,1001 @@
// 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 "PlotUtils.h"
#include "FitEvent.h"
#include "StatUtils.h"
// MOVE TO MB UTILS!
// This function is intended to be modified to enforce a consistent masking for
// all models.
TH2D* PlotUtils::SetMaskHist(std::string type, TH2D* data) {
TH2D* fMaskHist = (TH2D*)data->Clone("fMaskHist");
for (int xBin = 0; xBin < fMaskHist->GetNbinsX(); ++xBin) {
for (int yBin = 0; yBin < fMaskHist->GetNbinsY(); ++yBin) {
if (data->GetBinContent(xBin + 1, yBin + 1) == 0)
fMaskHist->SetBinContent(xBin + 1, yBin + 1, 0);
else
fMaskHist->SetBinContent(xBin + 1, yBin + 1, 0.5);
if (!type.compare("MB_numu_2D")) {
if (yBin == 19 && xBin < 8)
fMaskHist->SetBinContent(xBin + 1, yBin + 1, 1.0);
} else {
if (yBin == 19 && xBin < 11)
fMaskHist->SetBinContent(xBin + 1, yBin + 1, 1.0);
}
if (yBin == 18 && xBin < 3)
fMaskHist->SetBinContent(xBin + 1, yBin + 1, 1.0);
if (yBin == 17 && xBin < 2)
fMaskHist->SetBinContent(xBin + 1, yBin + 1, 1.0);
if (yBin == 16 && xBin < 1)
fMaskHist->SetBinContent(xBin + 1, yBin + 1, 1.0);
}
}
return fMaskHist;
};
// MOVE TO GENERAL UTILS?
bool PlotUtils::CheckObjectWithName(TFile* inFile, std::string substring) {
TIter nextkey(inFile->GetListOfKeys());
TKey* key;
while ((key = (TKey*)nextkey())) {
std::string test(key->GetName());
if (test.find(substring) != std::string::npos) return true;
}
return false;
};
// MOVE TO GENERAL UTILS?
std::string PlotUtils::GetObjectWithName(TFile* inFile, std::string substring) {
TIter nextkey(inFile->GetListOfKeys());
TKey* key;
std::string output = "";
while ((key = (TKey*)nextkey())) {
std::string test(key->GetName());
if (test.find(substring) != std::string::npos) output = test;
}
return output;
};
void PlotUtils::CreateNeutModeArray(TH1* hist, TH1* neutarray[]) {
for (int i = 0; i < 60; i++) {
neutarray[i] = (TH1*)hist->Clone(Form("%s_NMODE_%i", hist->GetName(), i));
}
return;
};
void PlotUtils::DeleteNeutModeArray(TH1* neutarray[]) {
for (int i = 0; i < 60; i++) {
delete neutarray[i];
}
return;
};
void PlotUtils::FillNeutModeArray(TH1D* hist[], int mode, double xval,
double weight) {
if (abs(mode) > 60) return;
hist[abs(mode)]->Fill(xval, weight);
return;
};
void PlotUtils::FillNeutModeArray(TH2D* hist[], int mode, double xval,
double yval, double weight) {
if (abs(mode) > 60) return;
hist[abs(mode)]->Fill(xval, yval, weight);
return;
};
THStack PlotUtils::GetNeutModeStack(std::string title, TH1* ModeStack[],
int option) {
(void)option;
THStack allmodes = THStack(title.c_str(), title.c_str());
for (int i = 0; i < 60; i++) {
allmodes.Add(ModeStack[i]);
}
// Credit to Clarence for copying all this out.
// CC
ModeStack[1]->SetTitle("CCQE");
ModeStack[1]->SetFillColor(kBlue);
// ModeStack[1]->SetFillStyle(3444);
ModeStack[1]->SetLineColor(kBlue);
ModeStack[2]->SetTitle("2p/2h Nieves");
ModeStack[2]->SetFillColor(kRed);
// ModeStack[2]->SetFillStyle(3344);
ModeStack[2]->SetLineColor(kRed);
// ModeStack[11]->SetTitle("#it{#nu + p #rightarrow l^{-} + p + #pi^{+}}");
ModeStack[11]->SetTitle("CC1#pi^{+} on p");
ModeStack[11]->SetFillColor(kGreen);
// ModeStack[11]->SetFillStyle(3004);
ModeStack[11]->SetLineColor(kGreen);
// ModeStack[12]->SetTitle("#it{#nu + n #rightarrow l^{-} + p + #pi^{0}}");
ModeStack[12]->SetTitle("CC1#pi^{0} on n");
ModeStack[12]->SetFillColor(kGreen + 3);
// ModeStack[12]->SetFillStyle(3005);
ModeStack[12]->SetLineColor(kGreen);
// ModeStack[13]->SetTitle("#it{#nu + n #rightarrow l^{-} + n + #pi^{+}}");
ModeStack[13]->SetTitle("CC1#pi^{+} on n");
ModeStack[13]->SetFillColor(kGreen - 2);
// ModeStack[13]->SetFillStyle(3004);
ModeStack[13]->SetLineColor(kGreen);
ModeStack[16]->SetTitle("CC coherent");
ModeStack[16]->SetFillColor(kBlue);
// ModeStack[16]->SetFillStyle(3644);
ModeStack[16]->SetLineColor(kBlue);
// ModeStack[17]->SetTitle("#it{#nu + n #rightarrow l^{-} + p + #gamma}");
ModeStack[17]->SetTitle("CC1#gamma");
ModeStack[17]->SetFillColor(kMagenta);
// ModeStack[17]->SetFillStyle(3001);
ModeStack[17]->SetLineColor(kMagenta);
ModeStack[21]->SetTitle("Multi #pi (1.3 < W < 2.0)");
ModeStack[21]->SetFillColor(kYellow);
// ModeStack[21]->SetFillStyle(3005);
ModeStack[21]->SetLineColor(kYellow);
// ModeStack[22]->SetTitle("#it{#nu + n #rightarrow l^{-} + p + #eta^{0}}");
ModeStack[22]->SetTitle("CC1#eta^{0} on n");
ModeStack[22]->SetFillColor(kYellow - 2);
// ModeStack[22]->SetFillStyle(3013);
ModeStack[22]->SetLineColor(kYellow - 2);
// ModeStack[23]->SetTitle("#it{#nu + n #rightarrow l^{-} + #Lambda +
// K^{+}}");
ModeStack[23]->SetTitle("CC1#Labda1K^{+}");
ModeStack[23]->SetFillColor(kYellow - 6);
// ModeStack[23]->SetFillStyle(3013);
ModeStack[23]->SetLineColor(kYellow - 6);
ModeStack[26]->SetTitle("DIS (W > 2.0)");
ModeStack[26]->SetFillColor(kRed);
// ModeStack[26]->SetFillStyle(3006);
ModeStack[26]->SetLineColor(kRed);
// NC
// ModeStack[31]->SetTitle("#it{#nu + n #rightarrow #nu + n + #pi^{0}}");
ModeStack[31]->SetTitle("NC1#pi^{0} on n");
ModeStack[31]->SetFillColor(kBlue);
// ModeStack[31]->SetFillStyle(3004);
ModeStack[31]->SetLineColor(kBlue);
// ModeStack[32]->SetTitle("#it{#nu + p #rightarrow #nu + p + #pi^{0}}");
ModeStack[32]->SetTitle("NC1#pi^{0} on p");
ModeStack[32]->SetFillColor(kBlue + 3);
// ModeStack[32]->SetFillStyle(3004);
ModeStack[32]->SetLineColor(kBlue + 3);
// ModeStack[33]->SetTitle("#it{#nu + n #rightarrow #nu + p + #pi^{-}}");
ModeStack[33]->SetTitle("NC1#pi^{-} on n");
ModeStack[33]->SetFillColor(kBlue - 2);
// ModeStack[33]->SetFillStyle(3005);
ModeStack[33]->SetLineColor(kBlue - 2);
// ModeStack[34]->SetTitle("#it{#nu + p #rightarrow #nu + n + #pi^{+}}");
ModeStack[34]->SetTitle("NC1#pi^{+} on p");
ModeStack[34]->SetFillColor(kBlue - 8);
// ModeStack[34]->SetFillStyle(3005);
ModeStack[34]->SetLineColor(kBlue - 8);
ModeStack[36]->SetTitle("NC Coherent");
ModeStack[36]->SetFillColor(kBlue + 8);
// ModeStack[36]->SetFillStyle(3644);
ModeStack[36]->SetLineColor(kBlue + 8);
// ModeStack[38]->SetTitle("#it{#nu + n #rightarrow #nu + n + #gamma}");
ModeStack[38]->SetTitle("NC1#gamma on n");
ModeStack[38]->SetFillColor(kMagenta);
// ModeStack[38]->SetFillStyle(3001);
ModeStack[38]->SetLineColor(kMagenta);
// ModeStack[39]->SetTitle("#it{#nu + p #rightarrow #nu + p + #gamma}");
ModeStack[39]->SetTitle("NC1#gamma on p");
ModeStack[39]->SetFillColor(kMagenta - 10);
// ModeStack[39]->SetFillStyle(3001);
ModeStack[39]->SetLineColor(kMagenta - 10);
ModeStack[41]->SetTitle("Multi #pi (1.3 < W < 2.0)");
ModeStack[41]->SetFillColor(kBlue - 10);
// ModeStack[41]->SetFillStyle(3005);
ModeStack[41]->SetLineColor(kBlue - 10);
// ModeStack[42]->SetTitle("#it{#nu + n #rightarrow #nu + n + #eta^{0}}");
ModeStack[42]->SetTitle("NC1#eta^{0} on n");
ModeStack[42]->SetFillColor(kYellow - 2);
// ModeStack[42]->SetFillStyle(3013);
ModeStack[42]->SetLineColor(kYellow - 2);
// ModeStack[43]->SetTitle("#it{#nu + p #rightarrow #nu + p + #eta^{0}}");
ModeStack[43]->SetTitle("NC1#eta^{0} on p");
ModeStack[43]->SetFillColor(kYellow - 4);
// ModeStack[43]->SetFillStyle(3013);
ModeStack[43]->SetLineColor(kYellow - 4);
// ModeStack[44]->SetTitle("#it{#nu + n #rightarrow #nu + #Lambda + K^{0}}");
ModeStack[44]->SetTitle("NC1#Lambda1K^{0} on n");
ModeStack[44]->SetFillColor(kYellow - 6);
// ModeStack[44]->SetFillStyle(3014);
ModeStack[44]->SetLineColor(kYellow - 6);
// ModeStack[45]->SetTitle("#it{#nu + p #rightarrow #nu + #Lambda + K^{+}}");
ModeStack[45]->SetTitle("NC1#Lambda1K^{+}");
ModeStack[45]->SetFillColor(kYellow - 10);
// ModeStack[45]->SetFillStyle(3014);
ModeStack[45]->SetLineColor(kYellow - 10);
ModeStack[46]->SetTitle("DIS (W > 2.0)");
ModeStack[46]->SetFillColor(kRed);
// ModeStack[46]->SetFillStyle(3006);
ModeStack[46]->SetLineColor(kRed);
// ModeStack[51]->SetTitle("#it{#nu + p #rightarrow #nu + p}");
ModeStack[51]->SetTitle("NC on p");
ModeStack[51]->SetFillColor(kBlack);
// ModeStack[51]->SetFillStyle(3444);
ModeStack[51]->SetLineColor(kBlack);
// ModeStack[52]->SetTitle("#it{#nu + n #rightarrow #nu + n}");
ModeStack[52]->SetTitle("NC on n");
ModeStack[52]->SetFillColor(kGray);
// ModeStack[52]->SetFillStyle(3444);
ModeStack[52]->SetLineColor(kGray);
return allmodes;
};
TLegend PlotUtils::GenerateStackLegend(THStack stack, int xlow, int ylow,
int xhigh, int yhigh) {
TLegend leg = TLegend(xlow, ylow, xhigh, yhigh);
TObjArray* histarray = stack.GetStack();
int nhist = histarray->GetEntries();
for (int i = 0; i < nhist; i++) {
TH1* hist = (TH1*)(histarray->At(i));
leg.AddEntry((hist), ((TH1*)histarray->At(i))->GetTitle(), "fl");
}
leg.SetName(Form("%s_LEG", stack.GetName()));
return leg;
};
void PlotUtils::ScaleNeutModeArray(TH1* hist[], double factor,
std::string option) {
for (int i = 0; i < 60; i++) {
if (hist[i]) hist[i]->Scale(factor, option.c_str());
}
return;
};
void PlotUtils::ResetNeutModeArray(TH1* hist[]) {
for (int i = 0; i < 60; i++) {
if (hist[i]) hist[i]->Reset();
}
return;
};
//********************************************************************
// This assumes the Enu axis is the x axis, as is the case for MiniBooNE 2D
// distributions
void PlotUtils::FluxUnfoldedScaling(TH2D* fMCHist, TH1D* fhist, TH1D* ehist,
double scalefactor) {
//********************************************************************
// Make clones to avoid changing stuff
TH1D* eventhist = (TH1D*)ehist->Clone();
TH1D* fFluxHist = (TH1D*)fhist->Clone();
// Undo width integral in SF
fMCHist->Scale(scalefactor /
eventhist->Integral(1, eventhist->GetNbinsX() + 1, "width"));
// Standardise The Flux
eventhist->Scale(1.0 / fFluxHist->Integral());
fFluxHist->Scale(1.0 / fFluxHist->Integral());
// Do interpolation for 2D plots?
// fFluxHist = PlotUtils::InterpolateFineHistogram(fFluxHist,100,"width");
// eventhist = PlotUtils::InterpolateFineHistogram(eventhist,100,"width");
// eventhist->Scale(1.0/fFluxHist->Integral());
// fFluxHist->Scale(1.0/fFluxHist->Integral());
// Scale fMCHist by eventhist integral
fMCHist->Scale(eventhist->Integral(1, eventhist->GetNbinsX() + 1));
// Now Get a flux PDF assuming X axis is Enu
TH1D* pdfflux = (TH1D*)fMCHist->ProjectionX()->Clone();
// pdfflux->Write( (std::string(fMCHist->GetName()) + "_PROJX").c_str());
pdfflux->Reset();
// Awful MiniBooNE Check for the time being
bool ismb =
std::string(fMCHist->GetName()).find("MiniBooNE") != std::string::npos;
for (int i = 0; i < pdfflux->GetNbinsX(); i++) {
double Ml = pdfflux->GetXaxis()->GetBinLowEdge(i + 1);
double Mh = pdfflux->GetXaxis()->GetBinLowEdge(i + 2);
// double Mc = pdfflux->GetXaxis()->GetBinCenter(i+1);
// double Mw = pdfflux->GetBinWidth(i+1);
double fluxint = 0.0;
// Scaling to match flux for MB
if (ismb) {
Ml /= 1.E3;
Mh /= 1.E3;
// Mc /= 1.E3;
// Mw /= 1.E3;
}
for (int j = 0; j < fFluxHist->GetNbinsX(); j++) {
// double Fc = fFluxHist->GetXaxis()->GetBinCenter(j+1);
double Fl = fFluxHist->GetXaxis()->GetBinLowEdge(j + 1);
double Fh = fFluxHist->GetXaxis()->GetBinLowEdge(j + 2);
double Fe = fFluxHist->GetBinContent(j + 1);
double Fw = fFluxHist->GetXaxis()->GetBinWidth(j + 1);
if (Fl >= Ml and Fh <= Mh) {
fluxint += Fe;
} else if (Fl < Ml and Fl < Mh and Fh > Ml and Fh < Mh) {
fluxint += Fe * (Fh - Ml) / Fw;
} else if (Fh > Mh and Fl < Mh and Fh > Ml and Fl > Ml) {
fluxint += Fe * (Mh - Fl) / Fw;
} else if (Ml >= Fl and Mh <= Fh) {
fluxint += Fe * (Mh - Ml) / Fw;
} else {
continue;
}
}
pdfflux->SetBinContent(i + 1, fluxint);
}
for (int i = 0; i < fMCHist->GetNbinsX(); i++) {
for (int j = 0; j < fMCHist->GetNbinsY(); j++) {
if (pdfflux->GetBinContent(i + 1) == 0.0) continue;
double binWidth = fMCHist->GetYaxis()->GetBinLowEdge(j + 2) -
fMCHist->GetYaxis()->GetBinLowEdge(j + 1);
fMCHist->SetBinContent(i + 1, j + 1,
fMCHist->GetBinContent(i + 1, j + 1) /
pdfflux->GetBinContent(i + 1) / binWidth);
fMCHist->SetBinError(i + 1, j + 1, fMCHist->GetBinError(i + 1, j + 1) /
pdfflux->GetBinContent(i + 1) /
binWidth);
}
}
delete eventhist;
delete fFluxHist;
return;
};
TH1D* PlotUtils::InterpolateFineHistogram(TH1D* hist, int res,
std::string opt) {
int nbins = hist->GetNbinsX();
double elow = hist->GetXaxis()->GetBinLowEdge(1);
double ehigh = hist->GetXaxis()->GetBinLowEdge(nbins + 1);
bool width = true; // opt.find("width") != std::string::npos;
TH1D* fine = new TH1D("fine", "fine", nbins * res, elow, ehigh);
TGraph* temp = new TGraph();
for (int i = 0; i < nbins; i++) {
double E = hist->GetXaxis()->GetBinCenter(i + 1);
double C = hist->GetBinContent(i + 1);
double W = hist->GetXaxis()->GetBinWidth(i + 1);
if (!width) W = 1.0;
if (W != 0.0) temp->SetPoint(temp->GetN(), E, C / W);
}
for (int i = 0; i < fine->GetNbinsX(); i++) {
double E = fine->GetXaxis()->GetBinCenter(i + 1);
double W = fine->GetBinWidth(i + 1);
if (!width) W = 1.0;
fine->SetBinContent(i + 1, temp->Eval(E, 0, "S") * W);
}
fine->Scale(hist->Integral(1, hist->GetNbinsX() + 1) /
fine->Integral(1, fine->GetNbinsX() + 1));
std::cout << "Interpolation Difference = "
<< fine->Integral(1, fine->GetNbinsX() + 1) << "/"
<< hist->Integral(1, hist->GetNbinsX() + 1) << std::endl;
return fine;
}
//********************************************************************
// This interpolates the flux by a TGraph instead of requiring the flux and MC
// flux to have the same binning
void PlotUtils::FluxUnfoldedScaling(TH1D* mcHist, TH1D* fhist, TH1D* ehist,
double scalefactor, int nevents) {
//********************************************************************
TH1D* eventhist = (TH1D*)ehist->Clone();
TH1D* fFluxHist = (TH1D*)fhist->Clone();
if (FitPar::Config().GetParB("save_flux_debug")) {
std::string name = std::string(mcHist->GetName());
mcHist->Write((name + "_UNF_MC").c_str());
fFluxHist->Write((name + "_UNF_FLUX").c_str());
eventhist->Write((name + "_UNF_EVT").c_str());
TH1D* scalehist = new TH1D("scalehist", "scalehist", 1, 0.0, 1.0);
scalehist->SetBinContent(1, scalefactor);
scalehist->SetBinContent(2, nevents);
scalehist->Write((name + "_UNF_SCALE").c_str());
}
// Undo width integral in SF
mcHist->Scale(scalefactor /
eventhist->Integral(1, eventhist->GetNbinsX() + 1, "width"));
// Standardise The Flux
eventhist->Scale(1.0 / fFluxHist->Integral());
fFluxHist->Scale(1.0 / fFluxHist->Integral());
// Scale mcHist by eventhist integral
mcHist->Scale(eventhist->Integral(1, eventhist->GetNbinsX() + 1));
// Now Get a flux PDF
TH1D* pdfflux = (TH1D*)mcHist->Clone();
pdfflux->Reset();
for (int i = 0; i < mcHist->GetNbinsX(); i++) {
double Ml = mcHist->GetXaxis()->GetBinLowEdge(i + 1);
double Mh = mcHist->GetXaxis()->GetBinLowEdge(i + 2);
// double Mc = mcHist->GetXaxis()->GetBinCenter(i+1);
// double Me = mcHist->GetBinContent(i+1);
// double Mw = mcHist->GetBinWidth(i+1);
double fluxint = 0.0;
for (int j = 0; j < fFluxHist->GetNbinsX(); j++) {
// double Fc = fFluxHist->GetXaxis()->GetBinCenter(j+1);
double Fl = fFluxHist->GetXaxis()->GetBinLowEdge(j + 1);
double Fh = fFluxHist->GetXaxis()->GetBinLowEdge(j + 2);
double Fe = fFluxHist->GetBinContent(j + 1);
double Fw = fFluxHist->GetXaxis()->GetBinWidth(j + 1);
if (Fl >= Ml and Fh <= Mh) {
fluxint += Fe;
} else if (Fl < Ml and Fl < Mh and Fh > Ml and Fh < Mh) {
fluxint += Fe * (Fh - Ml) / Fw;
} else if (Fh > Mh and Fl < Mh and Fh > Ml and Fl > Ml) {
fluxint += Fe * (Mh - Fl) / Fw;
} else if (Ml >= Fl and Mh <= Fh) {
fluxint += Fe * (Mh - Ml) / Fw;
} else {
continue;
}
}
pdfflux->SetBinContent(i + 1, fluxint);
}
// Scale MC hist by pdfflux
for (int i = 0; i < mcHist->GetNbinsX(); i++) {
if (pdfflux->GetBinContent(i + 1) == 0.0) continue;
mcHist->SetBinContent(
i + 1, mcHist->GetBinContent(i + 1) / pdfflux->GetBinContent(i + 1));
mcHist->SetBinError(
i + 1, mcHist->GetBinError(i + 1) / pdfflux->GetBinContent(i + 1));
}
delete eventhist;
delete fFluxHist;
return;
};
// MOVE TO GENERAL UTILS
//********************************************************************
void PlotUtils::Set2DHistFromText(std::string dataFile, TH2* hist, double norm,
bool skipbins) {
//********************************************************************
std::string line;
- std::ifstream data(dataFile.c_str(), ifstream::in);
+ std::ifstream data(dataFile.c_str(), std::ifstream::in);
int yBin = 0;
while (std::getline(data >> std::ws, line, '\n')) {
std::vector entries = GeneralUtils::ParseToDbl(line, " ");
// Loop over entries and insert them into the histogram
for (uint xBin = 0; xBin < entries.size(); xBin++) {
if (!skipbins or entries[xBin] != -1.0)
hist->SetBinContent(xBin + 1, yBin + 1, entries[xBin] * norm);
}
yBin++;
}
return;
}
// MOVE TO GENERAL UTILS
TH1D* PlotUtils::GetTH1DFromFile(std::string dataFile, std::string title,
std::string fPlotTitles,
std::string alt_name) {
TH1D* tempPlot;
// If format is a root file
if (dataFile.find(".root") != std::string::npos) {
TFile* temp_infile = new TFile(dataFile.c_str(), "READ");
tempPlot = (TH1D*)temp_infile->Get(title.c_str());
tempPlot->SetDirectory(0);
temp_infile->Close();
delete temp_infile;
// Else its a space seperated txt file
} else {
// Make a TGraph Errors
TGraphErrors* gr = new TGraphErrors(dataFile.c_str(), "%lg %lg %lg");
if (gr->IsZombie()) {
exit(-1);
}
double* bins = gr->GetX();
double* values = gr->GetY();
double* errors = gr->GetEY();
int npoints = gr->GetN();
// Fill the histogram from it
tempPlot = new TH1D(title.c_str(), title.c_str(), npoints - 1, bins);
for (int i = 0; i < npoints; ++i) {
tempPlot->SetBinContent(i + 1, values[i]);
// If only two columns are present in the input file, use the sqrt(values)
// as the error
// equivalent to assuming that the error is statistical
if (!errors[i])
tempPlot->SetBinError(i + 1, sqrt(values[i]));
else
tempPlot->SetBinError(i + 1, errors[i]);
}
delete gr;
}
// Allow alternate naming for root files
if (!alt_name.empty()) {
tempPlot->SetNameTitle(alt_name.c_str(), alt_name.c_str());
}
// Allow alternate axis titles
if (!fPlotTitles.empty()) {
tempPlot->SetNameTitle(
tempPlot->GetName(),
(std::string(tempPlot->GetTitle()) + fPlotTitles).c_str());
}
return tempPlot;
};
TH1D* PlotUtils::GetRatioPlot(TH1D* hist1, TH1D* hist2) {
// make copy of first hist
TH1D* new_hist = (TH1D*)hist1->Clone();
// Do bins and errors ourselves as scales can go awkward
for (int i = 0; i < new_hist->GetNbinsX(); i++) {
if (hist2->GetBinContent(i + 1) == 0.0) {
new_hist->SetBinContent(i + 1, 0.0);
}
new_hist->SetBinContent(
i + 1, hist1->GetBinContent(i + 1) / hist2->GetBinContent(i + 1));
new_hist->SetBinError(
i + 1, hist1->GetBinError(i + 1) / hist2->GetBinContent(i + 1));
}
return new_hist;
};
TH1D* PlotUtils::GetRenormalisedPlot(TH1D* hist1, TH1D* hist2) {
// make copy of first hist
TH1D* new_hist = (TH1D*)hist1->Clone();
if (hist1->Integral("width") == 0 or hist2->Integral("width") == 0) {
new_hist->Reset();
return new_hist;
}
Double_t scaleF = hist2->Integral("width") / hist1->Integral("width");
new_hist->Scale(scaleF);
return new_hist;
};
TH1D* PlotUtils::GetShapePlot(TH1D* hist1) {
// make copy of first hist
TH1D* new_hist = (TH1D*)hist1->Clone();
if (hist1->Integral("width") == 0) {
new_hist->Reset();
return new_hist;
}
Double_t scaleF1 = 1.0 / hist1->Integral("width");
new_hist->Scale(scaleF1);
return new_hist;
};
TH1D* PlotUtils::GetShapeRatio(TH1D* hist1, TH1D* hist2) {
TH1D* new_hist1 = GetShapePlot(hist1);
TH1D* new_hist2 = GetShapePlot(hist2);
// Do bins and errors ourselves as scales can go awkward
for (int i = 0; i < new_hist1->GetNbinsX(); i++) {
if (hist2->GetBinContent(i + 1) == 0) {
new_hist1->SetBinContent(i + 1, 0.0);
}
new_hist1->SetBinContent(i + 1, new_hist1->GetBinContent(i + 1) /
new_hist2->GetBinContent(i + 1));
new_hist1->SetBinError(
i + 1, new_hist1->GetBinError(i + 1) / new_hist2->GetBinContent(i + 1));
}
delete new_hist2;
return new_hist1;
};
TH2D* PlotUtils::GetCovarPlot(TMatrixDSym* cov, std::string name,
std::string title) {
TH2D* CovarPlot;
if (cov)
CovarPlot = new TH2D((*cov));
else
CovarPlot = new TH2D(name.c_str(), title.c_str(), 1, 0, 1, 1, 0, 1);
CovarPlot->SetName(name.c_str());
CovarPlot->SetTitle(title.c_str());
return CovarPlot;
}
TH2D* PlotUtils::GetFullCovarPlot(TMatrixDSym* cov, std::string name) {
return PlotUtils::GetCovarPlot(
cov, name + "_COV", name + "_COV;Bins;Bins;Covariance (#times10^{-76})");
}
TH2D* PlotUtils::GetInvCovarPlot(TMatrixDSym* cov, std::string name) {
return PlotUtils::GetCovarPlot(
cov, name + "_INVCOV",
name + "_INVCOV;Bins;Bins;Inv. Covariance (#times10^{-76})");
}
TH2D* PlotUtils::GetDecompCovarPlot(TMatrixDSym* cov, std::string name) {
return PlotUtils::GetCovarPlot(
cov, name + "_DECCOV",
name + "_DECCOV;Bins;Bins;Decomp Covariance (#times10^{-76})");
}
TH1D* PlotUtils::GetTH1DFromRootFile(std::string file, std::string name) {
if (name.empty()) {
std::vector tempfile = GeneralUtils::ParseToStr(file, ";");
file = tempfile[0];
name = tempfile[1];
}
TFile* rootHistFile = new TFile(file.c_str(), "READ");
TH1D* tempHist = (TH1D*)rootHistFile->Get(name.c_str())->Clone();
tempHist->SetDirectory(0);
rootHistFile->Close();
return tempHist;
}
TH2D* PlotUtils::GetTH2DFromRootFile(std::string file, std::string name) {
if (name.empty()) {
std::vector tempfile = GeneralUtils::ParseToStr(file, ";");
file = tempfile[0];
name = tempfile[1];
}
TFile* rootHistFile = new TFile(file.c_str(), "READ");
TH2D* tempHist = (TH2D*)rootHistFile->Get(name.c_str())->Clone();
tempHist->SetDirectory(0);
rootHistFile->Close();
delete rootHistFile;
return tempHist;
}
TH1* PlotUtils::GetTH1FromRootFile(std::string file, std::string name) {
if (name.empty()) {
std::vector tempfile = GeneralUtils::ParseToStr(file, ";");
file = tempfile[0];
name = tempfile[1];
}
TFile* rootHistFile = new TFile(file.c_str(), "READ");
if (!rootHistFile || rootHistFile->IsZombie()) {
THROW("Couldn't open root file: \"" << file << "\".");
}
TH1* tempHist = dynamic_cast(rootHistFile->Get(name.c_str())->Clone());
if (!tempHist) {
THROW("Couldn't retrieve: \"" << name << "\" from root file: \"" << file
<< "\".");
}
tempHist->SetDirectory(0);
rootHistFile->Close();
delete rootHistFile;
return tempHist;
}
/// Returns a vector of named TH1*s found in a single input file.
///
/// Expects a descriptor like: file.root[hist1|hist2|...]
std::vector PlotUtils::GetTH1sFromRootFile(
std::string const& descriptor) {
std::vector descriptors =
GeneralUtils::ParseToStr(descriptor, ",");
std::vector hists;
for (size_t d_it = 0; d_it < descriptors.size(); ++d_it) {
std::string& d = descriptors[d_it];
std::vector fname = GeneralUtils::ParseToStr(d, "[");
if (!fname.size() || !fname[0].length()) {
THROW("Couldn't find input file when attempting to parse : \""
<< d << "\". Expected input.root[hist1|hist2|...].");
}
if (fname[1][fname[1].length() - 1] == ']') {
fname[1] = fname[1].substr(0, fname[1].length() - 1);
}
std::vector histnames =
GeneralUtils::ParseToStr(fname[1], "|");
if (!histnames.size()) {
THROW(
"Couldn't find any histogram name specifiers when attempting to "
"parse "
": \""
<< fname[1] << "\". Expected hist1|hist2|...");
}
TFile* rootHistFile = new TFile(fname[0].c_str(), "READ");
if (!rootHistFile || rootHistFile->IsZombie()) {
THROW("Couldn't open root file: \"" << fname[0] << "\".");
}
for (size_t i = 0; i < histnames.size(); ++i) {
TH1* tempHist =
dynamic_cast(rootHistFile->Get(histnames[i].c_str())->Clone());
if (!tempHist) {
THROW("Couldn't retrieve: \"" << histnames[i] << "\" from root file: \""
<< fname[0] << "\".");
}
tempHist->SetDirectory(0);
hists.push_back(tempHist);
}
rootHistFile->Close();
}
return hists;
}
TH2D* PlotUtils::GetTH2DFromTextFile(std::string file) {
/// Contents should be
/// Low Edfe
return NULL;
}
void PlotUtils::AddNeutModeArray(TH1D* hist1[], TH1D* hist2[], double scaling) {
for (int i = 0; i < 60; i++) {
if (!hist2[i]) continue;
if (!hist1[i]) continue;
hist1[i]->Add(hist2[i], scaling);
}
return;
}
void PlotUtils::ScaleToData(TH1D* data, TH1D* mc, TH1I* mask) {
double scaleF = GetDataMCRatio(data, mc, mask);
mc->Scale(scaleF);
return;
}
void PlotUtils::MaskBins(TH1D* hist, TH1I* mask) {
for (int i = 0; i < hist->GetNbinsX(); i++) {
if (mask->GetBinContent(i + 1) <= 0.5) continue;
hist->SetBinContent(i + 1, 0.0);
hist->SetBinError(i + 1, 0.0);
LOG(REC) << "MaskBins: Set " << hist->GetName() << " Bin " << i + 1
<< " to 0.0 +- 0.0" << std::endl;
}
return;
}
void PlotUtils::MaskBins(TH2D* hist, TH2I* mask) {
for (int i = 0; i < hist->GetNbinsX(); i++) {
for (int j = 0; j < hist->GetNbinsY(); j++) {
if (mask->GetBinContent(i + 1, j + 1) <= 0.5) continue;
hist->SetBinContent(i + 1, j + 1, 0.0);
hist->SetBinError(i + 1, j + 1, 0.0);
LOG(REC) << "MaskBins: Set " << hist->GetName() << " Bin " << i + 1 << " "
<< j + 1 << " to 0.0 +- 0.0" << std::endl;
}
}
return;
}
double PlotUtils::GetDataMCRatio(TH1D* data, TH1D* mc, TH1I* mask) {
double rat = 1.0;
TH1D* newmc = (TH1D*)mc->Clone();
TH1D* newdt = (TH1D*)data->Clone();
if (mask) {
MaskBins(newmc, mask);
MaskBins(newdt, mask);
}
rat = newdt->Integral() / newmc->Integral();
return rat;
}
TH2D* PlotUtils::GetCorrelationPlot(TH2D* cov, std::string name) {
TH2D* cor = (TH2D*)cov->Clone();
cor->Reset();
for (int i = 0; i < cov->GetNbinsX(); i++) {
for (int j = 0; j < cov->GetNbinsY(); j++) {
if (cov->GetBinContent(i + 1, i + 1) != 0.0 and
cov->GetBinContent(j + 1, j + 1) != 0.0)
cor->SetBinContent(i + 1, j + 1,
cov->GetBinContent(i + 1, j + 1) /
(sqrt(cov->GetBinContent(i + 1, i + 1) *
cov->GetBinContent(j + 1, j + 1))));
}
}
if (!name.empty()) {
cor->SetNameTitle(name.c_str(), (name + ";;correlation").c_str());
}
cor->SetMinimum(-1);
cor->SetMaximum(1);
return cor;
}
TH2D* PlotUtils::GetDecompPlot(TH2D* cov, std::string name) {
TMatrixDSym* covarmat = new TMatrixDSym(cov->GetNbinsX());
for (int i = 0; i < cov->GetNbinsX(); i++)
for (int j = 0; j < cov->GetNbinsY(); j++)
(*covarmat)(i, j) = cov->GetBinContent(i + 1, j + 1);
TMatrixDSym* decompmat = StatUtils::GetDecomp(covarmat);
TH2D* dec = (TH2D*)cov->Clone();
for (int i = 0; i < cov->GetNbinsX(); i++)
for (int j = 0; j < cov->GetNbinsY(); j++)
dec->SetBinContent(i + 1, j + 1, (*decompmat)(i, j));
delete covarmat;
delete decompmat;
dec->SetNameTitle(name.c_str(), (name + ";;;decomposition").c_str());
return dec;
}
TH2D* PlotUtils::MergeIntoTH2D(TH1D* xhist, TH1D* yhist, std::string zname) {
std::vector xedges, yedges;
for (int i = 0; i < xhist->GetNbinsX() + 2; i++) {
xedges.push_back(xhist->GetXaxis()->GetBinLowEdge(i + 1));
}
for (int i = 0; i < yhist->GetNbinsX() + 2; i++) {
yedges.push_back(yhist->GetXaxis()->GetBinLowEdge(i + 1));
}
int nbinsx = xhist->GetNbinsX();
int nbinsy = yhist->GetNbinsX();
std::string name =
std::string(xhist->GetName()) + "_vs_" + std::string(yhist->GetName());
std::string titles = ";" + std::string(xhist->GetXaxis()->GetTitle()) + ";" +
std::string(yhist->GetXaxis()->GetTitle()) + ";" + zname;
TH2D* newplot = new TH2D(name.c_str(), (name + titles).c_str(), nbinsx,
&xedges[0], nbinsy, &yedges[0]);
return newplot;
}
//***************************************************
void PlotUtils::MatchEmptyBins(TH1D* data, TH1D* mc) {
//**************************************************
for (int i = 0; i < data->GetNbinsX(); i++) {
if (data->GetBinContent(i + 1) == 0.0 or data->GetBinError(i + 1) == 0.0)
mc->SetBinContent(i + 1, 0.0);
}
return;
}
//***************************************************
void PlotUtils::MatchEmptyBins(TH2D* data, TH2D* mc) {
//**************************************************
for (int i = 0; i < data->GetNbinsX(); i++) {
for (int j = 0; j < data->GetNbinsY(); j++) {
if (data->GetBinContent(i + 1, j + 1) == 0.0 or
data->GetBinError(i + 1, j + 1) == 0.0)
mc->SetBinContent(i + 1, j + 1, 0.0);
}
}
return;
}
//***************************************************
TH1D* PlotUtils::GetProjectionX(TH2D* hist, TH2I* mask) {
//***************************************************
TH2D* maskedhist = StatUtils::ApplyHistogramMasking(hist, mask);
TH1D* hist_X = maskedhist->ProjectionX();
delete maskedhist;
return hist_X;
}
//***************************************************
TH1D* PlotUtils::GetProjectionY(TH2D* hist, TH2I* mask) {
//***************************************************
TH2D* maskedhist = StatUtils::ApplyHistogramMasking(hist, mask);
TH1D* hist_Y = maskedhist->ProjectionY();
delete maskedhist;
return hist_Y;
}