diff --git a/app/PrepareNuwroEvents.cxx b/app/PrepareNuwroEvents.cxx
index c6daed4..2b44f8f 100644
--- a/app/PrepareNuwroEvents.cxx
+++ b/app/PrepareNuwroEvents.cxx
@@ -1,492 +1,492 @@
#include "event1.h"
#include <stdio.h>
#include <stdlib.h>
// Hopefully we don't need these as they're included above.
// #include "params_all.h"
// #include "params.h"
#include "FitLogger.h"
#include "PlotUtils.h"
#include "TFile.h"
#include "TH1D.h"
#include "TTree.h"
void printInputCommands(char *argv[]) {
std::cout << "[USAGE]: " << argv[0]
<< " [-h] [-f] [-F <FluxRootFile>,<FluxHistName>] [-o output.root] "
"inputfile.root [file2.root ...]"
<< std::endl
<< "\t-h : Print this message." << std::endl
<< "\t-f : Pass -f argument to '$ hadd' invocation." << std::endl
<< "\t-F : Read input flux from input descriptor." << std::endl
<< "\t-o : Write full output to a new file." << std::endl
<< std::endl;
};
void CreateRateHistograms(std::string inputs, bool force_out);
void HaddNuwroFiles(std::vector<std::string> &inputs, bool force_out);
bool outputNewFile = false;
std::string ofile = "";
bool haveFluxInputs = false;
struct FluxInputBlob {
FluxInputBlob(std::string _File, std::string _Hist, int _PDG,
double _Fraction)
: File(_File), Hist(_Hist), PDG(_PDG), Fraction(_Fraction) {}
std::string File;
std::string Hist;
int PDG;
double Fraction;
};
std::vector<FluxInputBlob> FluxInputs;
bool haddedFiles = false;
TH1D *F2D(TH1F *f) {
Double_t *bins = new Double_t[f->GetXaxis()->GetNbins() + 1];
for (Int_t bi_it = 0; bi_it < f->GetXaxis()->GetNbins(); ++bi_it) {
bins[bi_it] = f->GetXaxis()->GetBinLowEdge(bi_it + 1);
}
bins[f->GetXaxis()->GetNbins()] =
f->GetXaxis()->GetBinUpEdge(f->GetXaxis()->GetNbins());
TH1D *d = new TH1D((std::string(f->GetName()) + "_f").c_str(), f->GetTitle(),
f->GetXaxis()->GetNbins(), bins);
std::cout << "Converted TH1F with " << f->GetXaxis()->GetXbins()
<< " bins : " << std::endl;
for (Int_t bi_it = 0; bi_it < f->GetXaxis()->GetNbins() + 2; ++bi_it) {
d->SetBinContent(bi_it, f->GetBinContent(bi_it));
d->SetBinError(bi_it, f->GetBinError(bi_it));
std::cout << "\tF " << f->GetXaxis()->GetBinLowEdge(bi_it) << "--[" << bi_it
<< "]--" << f->GetXaxis()->GetBinUpEdge(bi_it) << ": "
<< f->GetBinContent(bi_it) << " D "
<< d->GetXaxis()->GetBinLowEdge(bi_it) << "--[" << bi_it << "]--"
<< d->GetXaxis()->GetBinUpEdge(bi_it) << ": "
<< d->GetBinContent(bi_it) << std::endl;
}
delete bins;
return d;
}
TH1D *GetTH1DFromFile(std::string const &rootFile,
std::string const &histName) {
TFile *inpFile = new TFile(rootFile.c_str(), "READ");
if (!inpFile || !inpFile->IsOpen()) {
NUIS_ABORT("Cannot open input root file: " << rootFile
- << " to read input histo.");
+ << " to read input histo.");
}
TH1D *histD = dynamic_cast<TH1D *>(inpFile->Get(histName.c_str()));
if (!histD) {
TH1F *histF = dynamic_cast<TH1F *>(inpFile->Get(histName.c_str()));
if (!histF) {
- NUIS_ABORT("Cannot find TH1D/F: " << histName << " in root file: " << rootFile
- << ".");
+ NUIS_ABORT("Cannot find TH1D/F: " << histName << " in root file: "
+ << rootFile << ".");
}
histD = F2D(histF);
} else {
histD = static_cast<TH1D *>(histD->Clone());
}
histD->SetDirectory(NULL);
inpFile->Close();
return histD;
}
//*******************************
int main(int argc, char *argv[]) {
//*******************************
// If No Arguments print commands
if (argc == 1) {
printInputCommands(argv);
return 0;
}
std::vector<std::string> inputfiles;
bool force_output = false;
// Get Inputs
for (int i = 1; i < argc; ++i) {
if (!std::strcmp(argv[i], "-h")) {
printInputCommands(argv);
return 0;
} else if (!std::strcmp(argv[i], "-f")) {
force_output = true;
} else if (!std::strcmp(argv[i], "-o")) {
outputNewFile = true;
ofile = argv[++i];
} else if (!std::strcmp(argv[i], "-F")) {
std::string inpLine = argv[++i];
std::vector<std::string> fluxInputDescriptor =
GeneralUtils::ParseToStr(inpLine, ",");
if ((fluxInputDescriptor.size() != 2) &&
(fluxInputDescriptor.size() != 3) &&
(fluxInputDescriptor.size() != 4)) {
NUIS_ABORT("Received -F argument with option: \""
- << inpLine
- << "\", was expecting "
- "<FluxRootFile>,<FluxHistName>[,PDG[,speciesFraction]].");
+ << inpLine
+ << "\", was expecting "
+ "<FluxRootFile>,<FluxHistName>[,PDG[,speciesFraction]].");
}
haveFluxInputs = true;
FluxInputs.push_back(
FluxInputBlob(fluxInputDescriptor[0], fluxInputDescriptor[1],
(fluxInputDescriptor.size() > 2)
? GeneralUtils::StrToInt(fluxInputDescriptor[2])
: 14,
(fluxInputDescriptor.size() > 3)
? GeneralUtils::StrToDbl(fluxInputDescriptor[3])
: 1));
if (!FluxInputs.back().File.length() ||
!FluxInputs.back().Hist.length()) {
NUIS_ABORT("Received -F argument with option: \""
- << inpLine
- << "\", was expecting "
- "<FluxRootFile>,<FluxHistName>[,PDG[,speciesFraction]].");
+ << inpLine
+ << "\", was expecting "
+ "<FluxRootFile>,<FluxHistName>[,PDG[,speciesFraction]].");
}
} else {
inputfiles.push_back(std::string(argv[i]));
}
}
// If one input file just create flux histograms
if (inputfiles.size() > (UInt_t)1) {
HaddNuwroFiles(inputfiles, force_output);
} else if (inputfiles.size() < (UInt_t)1) {
printInputCommands(argv);
}
CreateRateHistograms(inputfiles[0], force_output);
NUIS_LOG(FIT, "Finished NUWRO Prep.");
};
//*******************************
void CreateRateHistograms(std::string inputs, bool force_out) {
//*******************************
// Open root file
TFile *outRootFile = 0;
TTree *nuwrotree = 0;
if (!haddedFiles &&
outputNewFile) { // we need to make the new file and clone the tree.
TFile *inpFile = new TFile(inputs.c_str(), "READ");
if (!inpFile || !inpFile->IsOpen()) {
NUIS_ABORT("Cannot open input root file: " << inputs);
}
TTree *inpTree = dynamic_cast<TTree *>(inpFile->Get("treeout"));
if (!inpTree) {
NUIS_ABORT("Cannot find TTree \"treeout\" in input root file: "
- << inputs.c_str());
+ << inputs.c_str());
}
outRootFile = new TFile(ofile.c_str(), force_out ? "RECREATE" : "CREATE");
if (!outRootFile || !outRootFile->IsOpen()) {
NUIS_ABORT("Couldn't open root file: "
- << ofile << " for writing, does it already exist?");
+ << ofile << " for writing, does it already exist?");
}
nuwrotree = inpTree->CloneTree(-1, "fast");
nuwrotree->SetDirectory(outRootFile);
nuwrotree->Write(nuwrotree->GetName());
} else {
outRootFile = new TFile(inputs.c_str(), "UPDATE");
if (!outRootFile || !outRootFile->IsOpen()) {
NUIS_ABORT("Cannot open input root file: " << inputs);
}
nuwrotree = dynamic_cast<TTree *>(outRootFile->Get("treeout"));
if (!nuwrotree) {
NUIS_ABORT("Cannot find TTree \"treeout\" in input root file: "
- << inputs.c_str());
+ << inputs.c_str());
}
}
// Get Flux Histogram
event *evt = new event();
nuwrotree->SetBranchAddress("e", &evt);
nuwrotree->GetEntry(0);
int fluxtype = evt->par.beam_type;
std::map<int, TH1D *> fluxlist;
std::map<int, TH1D *> eventlist;
std::vector<int> allpdg;
std::map<int, int> nevtlist;
std::map<int, double> intxseclist;
// Did the input file have a mono-energetic flux?
bool isMono = false;
nevtlist[0] = 0.0;
intxseclist[0] = 0.0;
allpdg.push_back(0);
NUIS_LOG(FIT, "Nuwro fluxtype = " << fluxtype);
if (haveFluxInputs) {
double totalFraction = 0;
for (size_t flux_it = 0; flux_it < FluxInputs.size(); ++flux_it) {
FluxInputBlob &fb = FluxInputs[flux_it];
int pdg = fb.PDG;
TH1D *fluxHist = GetTH1DFromFile(fb.File, fb.Hist);
double pctg = fb.Fraction;
totalFraction += pctg;
double Elow = fluxHist->GetXaxis()->GetBinLowEdge(1);
double Ehigh = fluxHist->GetXaxis()->GetBinLowEdge(
fluxHist->GetXaxis()->GetNbins() + 1);
NUIS_LOG(FIT, "Adding new nuwro flux "
- << "pdg: " << pdg << " pctg: " << pctg << " Elow: " << Elow
- << " Ehigh: " << Ehigh);
+ << "pdg: " << pdg << " pctg: " << pctg
+ << " Elow: " << Elow << " Ehigh: " << Ehigh);
// Sort total flux plot
if (!fluxlist[0]) {
// Setup total flux
fluxlist[0] = (TH1D *)fluxHist->Clone();
fluxlist[0]->SetNameTitle("FluxHist", "FluxHist");
// Prep empty total events
eventlist[0] = (TH1D *)fluxHist->Clone();
eventlist[0]->SetNameTitle("EvtHist", "EvtHist");
eventlist[0]->Reset();
} else {
// Add up each new plot
fluxlist[0]->Add(fluxHist);
}
fluxHist->SetNameTitle(Form("nuwro_pdg%i_pct%f_Flux", pdg, pctg),
Form("nuwro_pdg%i_pct%f_Flux", pdg, pctg));
TH1D *eventplot = (TH1D *)fluxHist->Clone();
eventplot->SetNameTitle(Form("nuwro_pdg%i_pct%f_Evt", pdg, pctg),
Form("nuwro_pdg%i_pct%f_Evt", pdg, pctg));
eventplot->Reset();
fluxlist[pdg] = (TH1D *)fluxHist->Clone();
eventlist[pdg] = eventplot;
nevtlist[pdg] = 0;
intxseclist[pdg] = 0.0;
allpdg.push_back(pdg);
delete fluxHist;
}
if (fabs(totalFraction - 1) > 1E-5) {
- NUIS_ABORT(FTL, "Total species fraction for input flux histos = "
- << totalFraction << ", expected to sum to 1.");
+ NUIS_ABORT("Total species fraction for input flux histos = "
+ << totalFraction << ", expected to sum to 1.");
}
} else if (fluxtype == 0) {
std::string fluxstring = evt->par.beam_energy;
std::vector<double> fluxvals = GeneralUtils::ParseToDbl(fluxstring, " ");
int pdg = evt->par.beam_particle;
double Elow = double(fluxvals[0]) / 1000.0;
double Ehigh = double(fluxvals[1]) / 1000.0;
TH1D *fluxplot = NULL;
if (Elow > Ehigh)
isMono = true;
// For files produced with a flux distribution
if (!isMono) {
NUIS_LOG(FIT, "Adding new nuwro flux "
- << "pdg: " << pdg << " Elow: " << Elow
- << " Ehigh: " << Ehigh);
+ << "pdg: " << pdg << " Elow: " << Elow
+ << " Ehigh: " << Ehigh);
fluxplot =
new TH1D("fluxplot", "fluxplot", fluxvals.size() - 4, Elow, Ehigh);
for (uint j = 2; j < fluxvals.size(); j++) {
NUIS_LOG(DEB, j << " " << fluxvals[j]);
fluxplot->SetBinContent(j - 1, fluxvals[j]);
}
} else { // For monoenergetic fluxes
NUIS_LOG(FIT, "Adding mono-energetic nuwro flux "
- << "pdg: " << pdg << " E: " << Elow);
+ << "pdg: " << pdg << " E: " << Elow);
fluxplot = new TH1D("fluxplot", "fluxplot", 100, 0, Elow * 2);
fluxplot->SetBinContent(fluxplot->FindBin(Elow), 1);
}
// Setup total flux
fluxlist[0] = (TH1D *)fluxplot->Clone();
fluxlist[0]->SetNameTitle("FluxHist", "FluxHist");
// Prep empty total events
eventlist[0] = (TH1D *)fluxplot->Clone();
eventlist[0]->SetNameTitle("EvtHist", "EvtHist");
eventlist[0]->Reset();
fluxplot->SetNameTitle(Form("nuwro_pdg%i_Flux", pdg),
Form("nuwro_pdg%i_Flux", pdg));
TH1D *eventplot = (TH1D *)fluxplot->Clone();
eventplot->SetNameTitle(Form("nuwro_pdg%i_Evt", pdg),
Form("nuwro_pdg%i_Evt", pdg));
eventplot->Reset();
fluxlist[pdg] = fluxplot;
eventlist[pdg] = eventplot;
nevtlist[pdg] = 0;
intxseclist[pdg] = 0.0;
allpdg.push_back(pdg);
} else if (fluxtype == 1) {
std::string fluxstring = evt->par.beam_content;
std::vector<std::string> fluxlines =
GeneralUtils::ParseToStr(fluxstring, "\n");
for (uint i = 0; i < fluxlines.size(); i++) {
std::vector<double> fluxvals =
GeneralUtils::ParseToDbl(fluxlines[i], " ");
int pdg = int(fluxvals[0]);
double pctg = double(fluxvals[1]) / 100.0;
double Elow = double(fluxvals[2]) / 1000.0;
double Ehigh = double(fluxvals[3]) / 1000.0;
NUIS_LOG(FIT, "Adding new nuwro flux "
- << "pdg: " << pdg << " pctg: " << pctg << " Elow: " << Elow
- << " Ehigh: " << Ehigh);
+ << "pdg: " << pdg << " pctg: " << pctg
+ << " Elow: " << Elow << " Ehigh: " << Ehigh);
TH1D *fluxplot =
new TH1D("fluxplot", "fluxplot", fluxvals.size() - 4, Elow, Ehigh);
for (uint j = 4; j < fluxvals.size(); j++) {
fluxplot->SetBinContent(j + 1, fluxvals[j]);
}
// Sort total flux plot
if (!fluxlist[0]) {
// Setup total flux
fluxlist[0] = (TH1D *)fluxplot->Clone();
fluxlist[0]->SetNameTitle("FluxHist", "FluxHist");
// Prep empty total events
eventlist[0] = (TH1D *)fluxplot->Clone();
eventlist[0]->SetNameTitle("EvtHist", "EvtHist");
eventlist[0]->Reset();
} else {
// Add up each new plot
fluxlist[0]->Add(fluxplot);
}
fluxplot->SetNameTitle(Form("nuwro_pdg%i_pct%f_Flux", pdg, pctg),
Form("nuwro_pdg%i_pct%f_Flux", pdg, pctg));
TH1D *eventplot = (TH1D *)fluxplot->Clone();
eventplot->SetNameTitle(Form("nuwro_pdg%i_pct%f_Evt", pdg, pctg),
Form("nuwro_pdg%i_pct%f_Evt", pdg, pctg));
eventplot->Reset();
fluxlist[pdg] = fluxplot;
eventlist[pdg] = eventplot;
nevtlist[pdg] = 0;
intxseclist[pdg] = 0.0;
allpdg.push_back(pdg);
}
}
// Start main event loop to fill plots
int nevents = nuwrotree->GetEntries();
double Enu = 0.0;
double TotXSec = 0.0;
// double totaleventmode = 0.0;
// double totalevents = 0.0;
int pdg = 0;
int countwidth = nevents / 50.0;
countwidth = countwidth ? countwidth : 1;
for (int i = 0; i < nevents; i++) {
nuwrotree->GetEntry(i);
// Get Variables
Enu = evt->in[0].t / 1000.0;
TotXSec = evt->weight;
pdg = evt->in[0].pdg;
eventlist[0]->Fill(Enu);
eventlist[pdg]->Fill(Enu);
nevtlist[0] += 1;
nevtlist[pdg] += 1;
intxseclist[0] += TotXSec;
intxseclist[pdg] += TotXSec;
if (i % countwidth == 0) {
NUIS_LOG(FIT, "Processed " << i << " events "
- << " (" << int(i * 100.0 / nevents) << "%)"
- << " : E, W, PDG = " << Enu << ", " << TotXSec
- << ", " << pdg)
+ << " (" << int(i * 100.0 / nevents) << "%)"
+ << " : E, W, PDG = " << Enu << ", " << TotXSec
+ << ", " << pdg)
}
}
TH1D *zeroevents = (TH1D *)eventlist[0]->Clone();
outRootFile->cd();
// Loop over eventlist
for (uint i = 0; i < allpdg.size(); i++) {
int pdg = allpdg[i];
double AvgXSec = intxseclist[0] * 1E38 / double(nevtlist[0]);
NUIS_LOG(FIT, pdg << " Avg XSec = " << AvgXSec);
NUIS_LOG(FIT, pdg << " nevents = " << double(nevtlist[pdg]));
if (!isMono) {
// Convert events to PDF
eventlist[pdg]->Scale(1.0 / zeroevents->Integral("width"));
// Multiply by total predicted event rate
eventlist[pdg]->Scale(fluxlist[0]->Integral("width") * AvgXSec);
} else {
// If a mono-energetic flux was used, width should not be used
// The output is (now) forced to be flux = 1, evtrt = xsec (in 1E38 * nb
// cm^2)
eventlist[pdg]->Scale(1.0 / zeroevents->Integral());
eventlist[pdg]->Scale(fluxlist[0]->Integral() * AvgXSec);
}
// Save everything
fluxlist[pdg]->Write("", TObject::kOverwrite);
eventlist[pdg]->Write("", TObject::kOverwrite);
}
// Tidy up
outRootFile->Close();
fluxlist.clear();
eventlist.clear();
// Exit Program
return;
}
//*******************************
void HaddNuwroFiles(std::vector<std::string> &inputs, bool force_out) {
//*******************************
// Get output file name
std::string outputname = inputs[0];
// Make command line string
std::string cmd = "hadd ";
if (outputNewFile) {
cmd += ofile + " ";
outputname = ofile;
} else if (force_out) {
cmd += "-f ";
}
for (UInt_t i = 0; i < inputs.size(); i++) {
cmd += inputs[i] + " ";
}
NUIS_LOG(FIT, " Running HADD from PrepareNuwro: " << cmd);
// Start HADD
system(cmd.c_str());
// Return name of output file
inputs.clear();
inputs.push_back(outputname);
haddedFiles = true;
return;
}
diff --git a/cmake/GENIESetup.cmake b/cmake/GENIESetup.cmake
index 43d8096..2b67e48 100644
--- a/cmake/GENIESetup.cmake
+++ b/cmake/GENIESetup.cmake
@@ -1,252 +1,251 @@
# 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 <http://www.gnu.org/licenses/>.
################################################################################
# TODO
# check system for libxml2
# check whether we need the includes
# check if we can use a subset of the GENIE libraries
include(${CMAKE_SOURCE_DIR}/cmake/parseConfigApp.cmake)
################################################################################
# Check Dependencies
################################################################################
################################# GENIE ######################################
if(GENIE STREQUAL "")
cmessage(FATAL_ERROR "Variable GENIE is not defined. "
"The location of a pre-built GENIE install must be defined either as"
" $ cmake -DGENIE=/path/to/GENIE or as an environment variable"
" $ export GENIE=/path/to/GENIE")
endif()
if (BUILD_GEVGEN)
cmessage(STATUS "Building custom gevgen")
LIST(APPEND EXTRA_CXX_FLAGS -D__GEVGEN_ENABLED__)
endif()
execute_process(COMMAND genie-config --version
OUTPUT_VARIABLE GENIE_VER OUTPUT_STRIP_TRAILING_WHITESPACE)
cmessage(STATUS "genie_ver: ${GENIE_VER}")
if(GENIE_VER VERSION_GREATER 3.0.0)
set(GENIE_POST_R3 1)
string(REPLACE "." "" GENIE_VERSION ${GENIE_VER})
cmessage(STATUS "set genie_post_r3")
endif()
if(NOT GENIE_POST_R3)
LIST(APPEND EXTRA_CXX_FLAGS -DGENIE_PRE_R3)
cmessage(STATUS "setting genie_pre_r3 ${EXTRA_CXX_FLAGS}")
if(GENIE_EMPMEC_REWEIGHT)
cmessage(STATUS "Enable EMPMEC dials")
LIST(APPEND EXTRA_CXX_FLAGS -D__GENIE_EMP_MECRW_ENABLED)
endif()
else()
cmessage(STATUS "Enable EMPMEC dials")
LIST(APPEND EXTRA_CXX_FLAGS -D__GENIE_EMP_MECRW_ENABLED)
endif()
execute_process (COMMAND genie-config
--topsrcdir OUTPUT_VARIABLE GENIE_INCLUDES_DIR OUTPUT_STRIP_TRAILING_WHITESPACE)
#Allows for external override in the case where genie-config lies.
if(NOT DEFINED GENIE_LIB_DIR OR GENIE_LIB_DIR STREQUAL "")
#This looks like it should call libdir, but it strips the argument with -L from the response of --libs
GETLIBDIRS(genie-config --libs GENIE_LIB_DIR)
endif()
GETLIBS(genie-config --libs GENIE_LIBS)
cmessage(STATUS "GENIE version : ${GENIE_VERSION}")
cmessage(STATUS "GENIE libdir : ${GENIE_LIB_DIR}")
cmessage(STATUS "GENIE libs : ${GENIE_LIBS}")
string(REGEX MATCH "ReinSeghal" WASMATCHED ${GENIE_LIBS})
if(WASMATCHED AND GENIE_VERSION STREQUAL "210")
set(GENIE_SEHGAL ${GENIE_LIBS})
STRING(REPLACE "ReinSeghal" "ReinSehgal" GENIE_LIBS ${GENIE_SEHGAL})
cmessage(DEBUG "Fixed inconsistency in library naming: ${GENIE_LIBS}")
endif()
if(NOT USE_REWEIGHT)
SET(USING_GENIE_RW FALSE)
elseif(NOT GENIE_POST_R3)
LIST(FIND GENIE_LIBS GReWeight FOUND_GENIE_RW)
if(FOUND_GENIE_RW EQUAL -1)
cmessage(DEBUG "Did NOT find ReWeight library. Here are libs: ${GENIE_LIBS}")
SET(USING_GENIE_RW FALSE)
else()
SET(USING_GENIE_RW TRUE)
endif()
elseif(DEFINED GENIE_REWEIGHT AND NOT GENIE_REWEIGHT STREQUAL "")
LIST(FIND GENIE_LIBS GRwFwk FOUND_GENIE_RW)
if(FOUND_GENIE_RW EQUAL -1)
LIST(APPEND GENIE_LIBS GRwClc GRwFwk GRwIO)
cmessage(DEBUG "Force added ReWeight library. Here are libs: ${GENIE_LIBS}")
SET(USING_GENIE_RW TRUE)
else()
SET(USING_GENIE_RW FALSE)
endif()
endif()
if(USING_GENIE_RW)
cmessage(STATUS "Using GENIE ReWeight library.")
else()
cmessage(STATUS "Building without GENIE ReWeight support.")
- LIST(APPEND EXTRA_CXX_FLAGS -D__NO_GENIE_REWEIGHT__)
endif()
LIST(APPEND GENIE_LIBS -Wl,--end-group )
LIST(REVERSE GENIE_LIBS)
LIST(APPEND GENIE_LIBS -Wl,--start-group -Wl,--no-as-needed )
LIST(REVERSE GENIE_LIBS)
cmessage(DEBUG "GENIE_LIBS: ${GENIE_LIBS}")
################################ LHAPDF ######################################
if(LHAPDF_LIB STREQUAL "")
cmessage(FATAL_ERROR "Variable LHAPDF_LIB is not defined. The location of a pre-built lhapdf install must be defined either as $ cmake -DLHAPDF_LIB=/path/to/LHAPDF_libraries or as an environment variable $ export LHAPDF_LIB=/path/to/LHAPDF_libraries")
endif()
if(LHAPDF_INC STREQUAL "")
cmessage(FATAL_ERROR "Variable LHAPDF_INC is not defined. The location of a pre-built lhapdf install must be defined either as $ cmake -DLHAPDF_INC=/path/to/LHAPDF_includes or as an environment variable $ export LHAPDF_INC=/path/to/LHAPDF_includes")
endif()
if(LHAPATH STREQUAL "")
cmessage(FATAL_ERROR "Variable LHAPATH is not defined. The location of a the LHAPATH directory must be defined either as $ cmake -DLHAPATH=/path/to/LHAPATH or as an environment variable $ export LHAPATH=/path/to/LHAPATH")
endif()
################################ LIBXML ######################################
if(LIBXML2_LIB STREQUAL "")
GETLIBDIR(xml2-config --libs LIBXML2_LIB IGNORE_EMPTY_RESPONSE)
if(LIBXML2_LIB STREQUAL "")
message(WARNING "Variable LIBXML2_LIB is not defined, as xml2-config was found and didn't report a library include path, it is likely that libxml2.so can be found in the standard system location, lets hope so.")
endif()
endif()
if(LIBXML2_INC STREQUAL "")
GETINCDIR(xml2-config --cflags LIBXML2_INC)
if(LIBXML2_INC STREQUAL "")
message(FATAL_ERROR "Variable LIBXML2_INC is not defined and could not be found with xml2-config. The location of a pre-built libxml2 install must be defined either as $ cmake -DLIBXML2_INC=/path/to/LIBXML2_includes or as an environment variable $ export LIBXML2_INC=/path/to/LIBXML2_includes")
endif()
endif()
############################### log4cpp ######################################
if(LOG4CPP_LIB STREQUAL "")
GETLIBDIR(log4cpp-config --libs LOG4CPP_LIB)
if(LOG4CPP_LIB STREQUAL "")
message(FATAL_ERROR "Variable LOG4CPP_LIB is not defined and could not be found with log4cpp-config. The location of a pre-built log4cpp install must be defined either as $ cmake -DLOG4CPP_LIB=/path/to/LOG4CPP_libraries or as an environment variable $ export LOG4CPP_LIB=/path/to/LOG4CPP_libraries")
endif()
endif()
if(LOG4CPP_INC STREQUAL "")
GETINCDIR(log4cpp-config --cflags LOG4CPP_INC)
if(LOG4CPP_INC STREQUAL "")
message(FATAL_ERROR "Variable LOG4CPP_INC is not defined and could not be found with log4cpp-config. The location of a pre-built log4cpp install must be defined either as $ cmake -DLOG4CPP_INC=/path/to/LOG4CPP_includes or as an environment variable $ export LOG4CPP_INC=/path/to/LOG4CPP_includes")
endif()
endif()
################################################################################
# Set the compiler defines
LIST(APPEND EXTRA_CXX_FLAGS -D__GENIE_ENABLED__ -D__GENIE_VERSION__=${GENIE_VERSION})
LIST(APPEND EXTRA_LIBS ${GENIE_LIBS})
############################### GSL ######################################
if(GENIE_POST_R3)
if(GSL_LIB STREQUAL "")
GETLIBDIR(gsl-config --libs GSL_LIB)
if(GSL_LIB STREQUAL "")
message(FATAL_ERROR "Variable GSL_LIB is not defined and could not be found with gsl-config. The location of a pre-built gsl install must be defined either as $ cmake -DGSL_LIB=/path/to/GSL_libraries or as an environment variable $ export GSL_LIB=/path/to/GSL_libraries")
endif()
endif()
if(GSL_INC STREQUAL "")
GETINCDIR(gsl-config --cflags GSL_INC)
if(GSL_INC STREQUAL "")
message(FATAL_ERROR "Variable GSL_INC is not defined and could not be found with gsl-config. The location of a pre-built gsl install must be defined either as $ cmake -DGSL_INC=/path/to/GSL_includes or as an environment variable $ export GSL_INC=/path/to/GSL_includes")
endif()
endif()
GETLIBS(gsl-config --libs GSL_LIB_LIST)
if(USING_GENIE_RW AND GENIE_REWEIGHT STREQUAL "")
message(FATAL_ERROR "Variable GENIE_REWEIGHT is not defined. When using GENIE v3+, we require the reweight product to be built and accessible via the environment variable GENIE_REWEIGHT")
endif()
endif()
################################################################################
LIST(APPEND EXTRA_LIBS LHAPDF xml2 log4cpp)
LIST(APPEND EXTRA_LINK_DIRS
${GENIE_LIB_DIR}
${LHAPDF_LIB}
${LOG4CPP_LIB})
# Append only if we have found GENIE ReWeight
if(NOT GENIE_POST_R3)
LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES
${GENIE_INCLUDES_DIR}
${GENIE_INCLUDES_DIR}/GHEP
${GENIE_INCLUDES_DIR}/Ntuple)
if(USING_GENIE_RW)
LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${GENIE_INCLUDES_DIR}/ReWeight)
endif()
LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES
${GENIE_INCLUDES_DIR}/Apps
${GENIE_INCLUDES_DIR}/FluxDrivers
${GENIE_INCLUDES_DIR}/EVGDrivers
${LHAPDF_INC}
${LIBXML2_INC}
${LOG4CPP_INC})
else()
LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES
${GENIE_INCLUDES_DIR})
if(USING_GENIE_RW)
LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${GENIE_REWEIGHT}/src)
endif()
LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${GSL_INC}
${LHAPDF_INC}
${LIBXML2_INC}
${LOG4CPP_INC})
if(USING_GENIE_RW)
LIST(APPEND EXTRA_LINK_DIRS
${GENIE_REWEIGHT}/lib)
endif()
LIST(APPEND EXTRA_LINK_DIRS
${GSL_LIB}
)
LIST(APPEND EXTRA_LIBS ${GSL_LIB_LIST})
endif()
if(USE_PYTHIA8)
set(NEED_PYTHIA8 TRUE)
set(NEED_ROOTPYTHIA8 TRUE)
else()
set(NEED_PYTHIA6 TRUE)
set(NEED_ROOTPYTHIA6 TRUE)
endif()
set(NEED_ROOTEVEGEN TRUE)
SET(USE_GENIE TRUE CACHE BOOL "Whether to enable GENIE (reweight) support. Requires external libraries. <FALSE>" FORCE)
diff --git a/cmake/NEUTSetup.cmake b/cmake/NEUTSetup.cmake
index 3cd85c6..eb4dd92 100644
--- a/cmake/NEUTSetup.cmake
+++ b/cmake/NEUTSetup.cmake
@@ -1,145 +1,209 @@
# 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 <http://www.gnu.org/licenses/>.
################################################################################
-if(NEUT_ROOT STREQUAL "")
- cmessage(FATAL_ERROR "Variable NEUT_ROOT is not defined. Please export environment variable NEUT_ROOT or configure with -DNEUT_ROOT=/path/to/NEUT. This must be set to point to a prebuilt NEUT instance.")
-endif()
+include(cmake/parseConfigApp.cmake)
-if(CERN STREQUAL "")
- cmessage(FATAL_ERROR "Variable CERN is not defined. Please export environment variable CERN or configure with -DCERN=/path/to/CERNLIB. This must be set to point to a prebuilt CERNLIB instance.")
-endif()
+find_program(NEUTCONFIG NAMES neut-config)
-if(CERN_LEVEL STREQUAL "")
- cmessage(FATAL_ERROR "Variable CERN_LEVEL is not defined. Please export environment variable CERN_LEVEL or configure with -DCERN_LEVEL=XXXX (likely to be 2005).")
-endif()
+LIST(APPEND EXTRA_CXX_FLAGS -DNEED_FILL_NEUT_COMMONS)
-if(${NEUT_VERSION} VERSION_LESS 5.4.0)
- set(NEUT_LIB_DIR ${NEUT_ROOT}/lib/Linux_pc)
+SET(HAVENEUTCONFIG FALSE)
+# We are dealing with shiny NEUT
+if(NOT "${NEUTCONFIG}" STREQUAL "NEUTCONFIG-NOTFOUND")
+ SET(HAVENEUTCONFIG TRUE)
+ cmessage(STATUS "Found neut-config, using it to determine configuration.")
else()
- set(NEUT_LIB_DIR ${NEUT_ROOT}/lib)
+ cmessage(STATUS "Failed to find neut-config, assuming older NEUT build.")
endif()
-set(NEUT_CLASS ${NEUT_ROOT}/src/neutclass)
-
-LIST(APPEND EXTRA_CXX_FLAGS -D__NEUT_ENABLED__ -DNEUT_VERSION=${NEUT_VERSION})
-
-LIST(APPEND EXTRA_CXX_FLAGS
- -I${NEUT_ROOT}/include
- -I${NEUT_ROOT}/src/neutclass
- -I${NEUT_ROOT}/src/reweight)
-
-LIST(APPEND EXTRA_LINK_DIRS
- ${NEUT_LIB_DIR}
- ${CERN}/${CERN_LEVEL}/lib
- ${NEUT_ROOT}/src/reweight)
-
-if(${NEUT_VERSION} VERSION_EQUAL 5.4.2)
- LIST(APPEND EXTRA_LIBS
- -Wl,--as-needed
- NReWeight
- -Wl,--start-group
- neutcore_5.4.2
- nuccorspl_5.4.2 #typo in NEUT, may hopefully disappear
- nuceff_5.4.2
- partnuck_5.4.2
- skmcsvc_5.4.2
- tauola_5.4.2
- HT2p2h_5.4.0
- N1p1h_5.4.0
- -Wl,--end-group
- jetset74
- pdflib804
- mathlib
- packlib
- pawlib)
-
- LIST(APPEND EXTRA_CXX_FLAGS -DNEUT_COMMON_QEAV)
-elseif(${NEUT_VERSION} VERSION_EQUAL 5.4.0)
- LIST(APPEND EXTRA_LIBS
- -Wl,--as-needed
- NReWeight
- -Wl,--start-group
- neutcore_5.4.0
- nuccorspl_5.4.0 #typo in NEUT, may hopefully disappear
- nuceff_5.4.0
- partnuck_5.4.0
- skmcsvc_5.4.0
- tauola_5.4.0
- HT2p2h_5.4.0
- N1p1h_5.4.0
- specfunc_5.4.0
- radcorr_5.4.0
- gfortran
- -Wl,--end-group
- jetset74
- pdflib804
- mathlib
- packlib
- pawlib)
-else()
- LIST(APPEND EXTRA_LIBS
- -Wl,--as-needed
- NReWeight
- -Wl,--start-group
- neutcore
- nuccorrspl
- nuceff
- partnuck
- skmcsvc
- tauola
- -Wl,--end-group
- jetset74
- pdflib804
- mathlib
- packlib
- pawlib)
-endif()
+if(HAVENEUTCONFIG)
+ execute_process (COMMAND neut-config
+ --version OUTPUT_VARIABLE NEUT_VERSION
+ OUTPUT_STRIP_TRAILING_WHITESPACE)
+
+ execute_process (COMMAND neut-config
+ --incdir OUTPUT_VARIABLE NEUT_INCLUDE_DIRS
+ OUTPUT_STRIP_TRAILING_WHITESPACE)
+
+ execute_process (COMMAND neut-config
+ --libdir OUTPUT_VARIABLE NEUT_LINK_DIRS
+ OUTPUT_STRIP_TRAILING_WHITESPACE)
+
+ GETLIBDIRS(neut-config --cernflags CERN_LIB_DIR)
+ LIST(APPEND NEUT_LINK_DIRS ${CERN_LIB_DIR})
+ GETLIBS(neut-config --cernflags CERN_LIBS)
+
+ if(USE_REWEIGHT)
+ execute_process (COMMAND neut-config
+ --rwlibflags OUTPUT_VARIABLE NEUT_RWLIBS
+ OUTPUT_STRIP_TRAILING_WHITESPACE)
+ GETLIBS(neut-config --rwlibflags NEUT_RWLIBS)
+ LIST(APPEND NEUT_LIBS ${NEUT_RWLIBS})
+ else()
+ GETLIBS(neut-config --libflags NEUT_GENLIBS)
+ GETLIBS(neut-config --iolibflags NEUT_LIBS)
+ LIST(APPEND NEUT_LIBS ${NEUT_IOLIBS})
+ LIST(APPEND NEUT_LIBS ${NEUT_GENLIBS})
+ endif()
+
+ LIST(APPEND NEUT_LIBS ${CERN_LIBS};gfortran)
+ LIST(APPEND EXTRA_LIBS ${NEUT_LIBS})
+
+ PrefixList(NEUT_INCLUDE_DIRS "-I" ${NEUT_INCLUDE_DIRS})
+ LIST(APPEND EXTRA_CXX_FLAGS ${NEUT_INCLUDE_DIRS} -D__NEUT_ENABLED__ -DNEUT_VERSION=${NEUT_VERSION})
+
+ LIST(APPEND EXTRA_LINK_DIRS ${NEUT_LINK_DIRS})
+
+ cmessage(STATUS "NEUT")
+ cmessage(STATUS " Version : ${NEUT_VERSION}")
+ cmessage(STATUS " Flags : ${NEUT_CXX_FLAGS}")
+ cmessage(STATUS " Includes : ${NEUT_INCLUDE_DIRS}")
+ cmessage(STATUS " Link Dirs : ${NEUT_LINK_DIRS}")
+ cmessage(STATUS " Libs : ${NEUT_LIBS}")
+
+
+else() # Everything better be set up already
+
+ if(NEUT_ROOT STREQUAL "")
+ cmessage(FATAL_ERROR "Variable NEUT_ROOT is not defined. Please export environment variable NEUT_ROOT or configure with -DNEUT_ROOT=/path/to/NEUT. This must be set to point to a prebuilt NEUT instance.")
+ endif()
+
+ if(CERN STREQUAL "")
+ cmessage(FATAL_ERROR "Variable CERN is not defined. Please export environment variable CERN or configure with -DCERN=/path/to/CERNLIB. This must be set to point to a prebuilt CERNLIB instance.")
+ endif()
+
+ if(CERN_LEVEL STREQUAL "")
+ cmessage(FATAL_ERROR "Variable CERN_LEVEL is not defined. Please export environment variable CERN_LEVEL or configure with -DCERN_LEVEL=XXXX (likely to be 2005).")
+ endif()
+
+ if(${NEUT_VERSION} VERSION_LESS 5.4.0)
+ set(NEUT_LIB_DIR ${NEUT_ROOT}/lib/Linux_pc)
+ else()
+ set(NEUT_LIB_DIR ${NEUT_ROOT}/lib)
+ endif()
+
+ set(NEUT_CLASS ${NEUT_ROOT}/src/neutclass)
+
+ LIST(APPEND EXTRA_CXX_FLAGS -D__NEUT_ENABLED__ -DNEUT_VERSION=${NEUT_VERSION})
+
+ LIST(APPEND EXTRA_CXX_FLAGS
+ -I${NEUT_ROOT}/include
+ -I${NEUT_ROOT}/src/neutclass
+ -I${NEUT_ROOT}/src/reweight)
+
+ LIST(APPEND EXTRA_LINK_DIRS
+ ${NEUT_LIB_DIR}
+ ${CERN}/${CERN_LEVEL}/lib
+ ${NEUT_ROOT}/src/reweight)
+
+ if(${NEUT_VERSION} VERSION_EQUAL 5.4.2)
+ LIST(APPEND EXTRA_LIBS
+ -Wl,--as-needed
+ NReWeight
+ -Wl,--start-group
+ neutcore_5.4.2
+ nuccorspl_5.4.2 #typo in NEUT, may hopefully disappear
+ nuceff_5.4.2
+ partnuck_5.4.2
+ skmcsvc_5.4.2
+ tauola_5.4.2
+ HT2p2h_5.4.0
+ N1p1h_5.4.0
+ -Wl,--end-group
+ jetset74
+ pdflib804
+ mathlib
+ packlib
+ pawlib)
+
+ LIST(APPEND EXTRA_CXX_FLAGS -DNEUT_COMMON_QEAV)
+ elseif(${NEUT_VERSION} VERSION_EQUAL 5.4.0)
+ LIST(APPEND EXTRA_LIBS
+ -Wl,--as-needed
+ NReWeight
+ -Wl,--start-group
+ neutcore_5.4.0
+ nuccorspl_5.4.0 #typo in NEUT, may hopefully disappear
+ nuceff_5.4.0
+ partnuck_5.4.0
+ skmcsvc_5.4.0
+ tauola_5.4.0
+ HT2p2h_5.4.0
+ N1p1h_5.4.0
+ specfunc_5.4.0
+ radcorr_5.4.0
+ gfortran
+ -Wl,--end-group
+ jetset74
+ pdflib804
+ mathlib
+ packlib
+ pawlib)
+ else()
+ LIST(APPEND EXTRA_LIBS
+ -Wl,--as-needed
+ NReWeight
+ -Wl,--start-group
+ neutcore
+ nuccorrspl
+ nuceff
+ partnuck
+ skmcsvc
+ tauola
+ -Wl,--end-group
+ jetset74
+ pdflib804
+ mathlib
+ packlib
+ pawlib)
+ endif()
+
+ set(NEUT_ROOT_LIBS)
-set(NEUT_ROOT_LIBS)
+ LIST(APPEND NEUT_ROOT_LIBS
+ ${NEUT_CLASS}/neutctrl.so
+ ${NEUT_CLASS}/neutfsivert.so)
-LIST(APPEND NEUT_ROOT_LIBS
- ${NEUT_CLASS}/neutctrl.so
- ${NEUT_CLASS}/neutfsivert.so)
+ # Check for new versions of NEUT with NUCLEON FSI
+ if(EXISTS "${NEUT_CLASS}/neutnucfsistep.so")
+ set(NEUT_NUCFSI 1)
+ LIST(APPEND EXTRA_CXX_FLAGS -DNEUT_NUCFSI_ENABLED)
-# Check for new versions of NEUT with NUCLEON FSI
-if(EXISTS "${NEUT_CLASS}/neutnucfsistep.so")
- set(NEUT_NUCFSI 1)
- LIST(APPEND EXTRA_CXX_FLAGS -DNEUT_NUCFSI_ENABLED)
+ LIST(APPEND NEUT_ROOT_LIBS
+ ${NEUT_CLASS}/neutnucfsistep.so
+ ${NEUT_CLASS}/neutnucfsivert.so
+ )
+ endif()
+
+ if(${NEUT_VERSION} VERSION_LESS 5.4.0)
+ LIST(APPEND NEUT_ROOT_LIBS
+ ${NEUT_CLASS}/neutrootTreeSingleton.so)
+ endif()
LIST(APPEND NEUT_ROOT_LIBS
- ${NEUT_CLASS}/neutnucfsistep.so
- ${NEUT_CLASS}/neutnucfsivert.so
+ ${NEUT_CLASS}/neutvtx.so
+ ${NEUT_CLASS}/neutfsipart.so
+ ${NEUT_CLASS}/neutpart.so
+ ${NEUT_CLASS}/neutvect.so
)
-endif()
-if(${NEUT_VERSION} VERSION_LESS 5.4.0)
- LIST(APPEND NEUT_ROOT_LIBS
- ${NEUT_CLASS}/neutrootTreeSingleton.so)
+ foreach(OBJ ${NEUT_ROOT_LIBS})
+ LIST(APPEND EXTRA_SHAREDOBJS ${OBJ})
+ endforeach()
endif()
-
-LIST(APPEND NEUT_ROOT_LIBS
- ${NEUT_CLASS}/neutvtx.so
- ${NEUT_CLASS}/neutfsipart.so
- ${NEUT_CLASS}/neutpart.so
- ${NEUT_CLASS}/neutvect.so
- )
-
-foreach(OBJ ${NEUT_ROOT_LIBS})
- LIST(APPEND EXTRA_SHAREDOBJS ${OBJ})
-endforeach()
diff --git a/cmake/NuWroSetup.cmake b/cmake/NuWroSetup.cmake
index f649bff..af9d455 100644
--- a/cmake/NuWroSetup.cmake
+++ b/cmake/NuWroSetup.cmake
@@ -1,112 +1,111 @@
# 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 <http://www.gnu.org/licenses/>.
################################################################################
if(NOT NUWRO_INPUT_FILE STREQUAL "")
if(NOT EXISTS ${NUWRO_INPUT_FILE})
cmessage(FATAL_ERROR "Expected -DNUWRO_INPUT_FILE to point to a valid input file. Cannot find: '${NUWRO_INPUT_FILE}'")
endif()
if(CMAKE_BUILD_TYPE MATCHES DEBUG)
BuildROOTProject(NuWro_event1 ${NUWRO_INPUT_FILE} "event,vec,vect,particle,flags,params,line" STATIC)
SET(ROOTLIBNAME "libNuWro_event1.a")
else(CMAKE_BUILD_TYPE MATCHES RELEASE)
BuildROOTProject(NuWro_event1 ${NUWRO_INPUT_FILE} "event,vec,vect,particle,flags,params,line" SHARED)
SET(ROOTLIBNAME "libNuWro_event1.so")
endif()
ADD_CUSTOM_TARGET(NuWro_event1HeaderLink ALL
COMMAND ${CMAKE_COMMAND} -E create_symlink
${CMAKE_BINARY_DIR}/NuWro_event1/event.h
${CMAKE_BINARY_DIR}/NuWro_event1/event1.h
DEPENDS NuWro_event1)
LIST(APPEND EXTRA_CXX_FLAGS -D__NUWRO_ENABLED__)
LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${CMAKE_BINARY_DIR}/NuWro_event1)
LIST(APPEND EXTRA_LINK_DIRS ${CMAKE_BINARY_DIR})
LIST(APPEND EXTRA_LIBS NuWro_event1)
LIST(APPEND PROJECTWIDE_EXTRA_DEPENDENCIES NuWro_event1HeaderLink )
install(TARGETS NuWro_event1 DESTINATION lib)
SET(NUWRO_BUILT_FROM_FILE TRUE)
else()
if(NUWRO STREQUAL "")
cmessage(FATAL_ERROR "Variable NUWRO is not defined. "
"This must be set to point to a prebuilt NuWro instance.")
endif()
# If you are using a version of NuWro without reweighting use this to compile.
if(USE_NuWro_RW)
if(NUWRO_INC STREQUAL "")
cmessage(FATAL_ERROR "Variable NUWRO_INC is not defined. "
"This must be set to point to an installed NuWro instance.")
endif()
LIST(APPEND EXTRA_CXX_FLAGS -D__NUWRO_ENABLED__ -D__NUWRO_REWEIGHT_ENABLED__)
if(USE_NuWro_SRW_Event)
LIST(APPEND EXTRA_CXX_FLAGS -D__USE_NUWRO_SRW_EVENTS__)
endif()
LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES
${NUWRO}/src
${NUWRO}/src/reweight
${NUWRO_INC}/nuwro)
LIST(APPEND EXTRA_LINK_DIRS ${NUWRO}/build/${CMAKE_SYSTEM_NAME}/lib)
LIST(APPEND EXTRA_LIBS reweight event)
else ()
LIST(APPEND EXTRA_CXX_FLAGS -D__NUWRO_ENABLED__)
LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${NUWRO}/src)
if(NOT EXISTS ${NUWRO}/bin/event1.so)
if(EXISTS ${NUWRO}/build/${CMAKE_SYSTEM_NAME}/lib)
if(NUWRO_INC STREQUAL "")
cmessage(FATAL_ERROR "Variable NUWRO_INC is not defined. "
"This must be set to point to an installed NuWro instance.")
endif()
LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${NUWRO_INC}/nuwro)
LIST(APPEND EXTRA_LINK_DIRS ${NUWRO}/build/${CMAKE_SYSTEM_NAME}/lib)
- LIST(APPEND EXTRA_LIBS event)
+ LIST(APPEND EXTRA_LIBS -Wl,--no-as-needed event -Wl,--as-needed)
else()
cmessage(FATAL_ERROR "Expected to find the NuWro event library in: ${NUWRO}/bin/event1.so, or if using NuWro with reweight support: ${NUWRO}/build/${CMAKE_SYSTEM_NAME}/lib/libevent.a. Is NuWro built?")
endif()
else()
- LIST(APPEND EXTRA_SHAREDOBJS ${NUWRO}/bin/event1.so)
+ LIST(APPEND EXTRA_SHAREDOBJS -Wl,--no-as-needed ${NUWRO}/bin/event1.so -Wl,--as-needed)
endif()
endif()
set(NEED_PYTHIA6 TRUE)
set(NEED_ROOTPYTHIA6 TRUE)
endif()
-
diff --git a/cmake/ReweightEnginesSetup.cmake b/cmake/ReweightEnginesSetup.cmake
index f72d2ac..31dc580 100644
--- a/cmake/ReweightEnginesSetup.cmake
+++ b/cmake/ReweightEnginesSetup.cmake
@@ -1,93 +1,97 @@
# 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 <http://www.gnu.org/licenses/>.
################################################################################
+if(NOT USE_REWEIGHT)
+ LIST(APPEND EXTRA_CXX_FLAGS -D__NO_REWEIGHT__)
+endif()
+
################################## T2K ######################################
if(USE_T2K)
include(${CMAKE_SOURCE_DIR}/cmake/T2KSetup.cmake)
cmessage(STATUS "Using T2K Reweight engine.")
set(USE_T2K TRUE CACHE BOOL "Whether to enable T2KReWeight support. Requires external libraries. <FALSE>" FORCE)
endif()
################################## NIWG ######################################
if(USE_NIWG)
include(${CMAKE_SOURCE_DIR}/cmake/NIWGSetup.cmake)
cmessage(STATUS "Using NIWG Reweight engine.")
set(USE_NIWG TRUE CACHE BOOL "Whether to enable (T2K) NIWG ReWeight support. Requires external libraries. <FALSE>" FORCE)
endif()
################################## MINERvA ######################################
if(USE_MINERvA_RW)
include(${CMAKE_SOURCE_DIR}/cmake/MINERvASetup.cmake)
cmessage(STATUS "Using MINERvA Reweight engine.")
set(USE_MINERvA_RW TRUE CACHE BOOL "Whether to enable MINERvA ReWeight support. <FALSE>" FORCE)
endif()
################################## NEUT ######################################
if(USE_NEUT)
include(${CMAKE_SOURCE_DIR}/cmake/NEUTSetup.cmake)
cmessage(STATUS "Using NEUT Reweight engine.")
set(USE_NEUT TRUE CACHE BOOL "Whether to enable NEUT (reweight) support. Requires external libraries. <FALSE>" FORCE)
endif()
################################# NuWro ######################################
if(USE_NuWro)
include(${CMAKE_SOURCE_DIR}/cmake/NuWroSetup.cmake)
cmessage(STATUS "Using NuWro Reweight engine.")
set(USE_NuWro TRUE CACHE BOOL "Whether to enable NuWro support. <FALSE>" FORCE)
endif()
################################## GENIE #####################################
if(USE_GENIE)
include(${CMAKE_SOURCE_DIR}/cmake/GENIESetup.cmake)
cmessage(STATUS "Using GENIE.")
set(USE_GENIE TRUE CACHE BOOL "Whether to enable GENIE support. Requires external libraries. <FALSE>" FORCE)
endif()
################################################################################
################################ NOvARwgt ####################################
if(USE_NOvARwgt)
include(${CMAKE_SOURCE_DIR}/cmake/NOvARwgtSetup.cmake)
cmessage(STATUS "Using NOvARwgt Reweight engine.")
set(USE_NOvARwgt TRUE CACHE BOOL "Whether to enable NOvARwgt (reweight) support. Requires external libraries. <FALSE>" FORCE)
endif()
################################################################################
################################ Prob3++ ####################################
include(${CMAKE_SOURCE_DIR}/cmake/Prob3++Setup.cmake)
################################################################################
cmessage(STATUS "Reweight engine include directories: ${RWENGINE_INCLUDE_DIRECTORIES}")
if(NEED_ROOTEVEGEN)
cmessage(STATUS "Require ROOT eve generation libraries")
LIST(REVERSE ROOT_LIBS)
LIST(APPEND ROOT_LIBS
Gui
Ged
Geom
TreePlayer
EG
Eve)
LIST(REVERSE ROOT_LIBS)
endif()
if(NEED_ROOTPYTHIA6)
cmessage(STATUS "Require ROOT Pythia6 libraries")
LIST(APPEND ROOT_LIBS
EGPythia6)
endif()
LIST(APPEND EXTRA_LIBS ${ROOT_LIBS})
diff --git a/cmake/parseConfigApp.cmake b/cmake/parseConfigApp.cmake
index 2f72ce5..76238c0 100644
--- a/cmake/parseConfigApp.cmake
+++ b/cmake/parseConfigApp.cmake
@@ -1,148 +1,185 @@
function(GETFIRSTMATCHINGDELMIMMEDDIR DELIM CONFIGAPP ARG DIR_OUT FAILURE_IS_NOT_ERROR)
- cmessage(DEBUG "CONFIGAPP: ${CONFIGAPP}, ARG: ${ARG} DIR_OUT: ${DIR_OUT}, FAILURE_IS_NOT_ERROR: ${FAILURE_IS_NOT_ERROR}")
-
if(DELIM STREQUAL "")
cmessage(FATAL_ERROR "GETFIRSTMATCHINGDELMIMMEDDIR Passed no delimiter. This is a build configuration bug in NUISANCE, please report to the developers.")
endif()
if(CONFIGAPP STREQUAL "")
cmessage(FATAL_ERROR "GETFIRSTMATCHINGDELMIMMEDDIR Passed no configuration application. This is a build configuration bug in NUISANCE, please report to the developers.")
endif()
SET(CONFIGAPP_LOCATION "CONFIGAPP_LOCATION-NOTFOUND")
find_program(CONFIGAPP_LOCATION ${CONFIGAPP})
if(NOT CONFIGAPP_LOCATION STREQUAL "CONFIGAPP_LOCATION-NOTFOUND")
execute_process (COMMAND ${CONFIGAPP_LOCATION}
${ARG} OUTPUT_VARIABLE CONFIGAPP_RESPONSE_RAW OUTPUT_STRIP_TRAILING_WHITESPACE)
cmessage(DEBUG "${CONFIGAPP_LOCATION} ${ARG} responded with: \"${CONFIGAPP_RESPONSE_RAW}\"")
if(CONFIGAPP_RESPONSE_RAW STREQUAL "")
if(FAILURE_IS_NOT_ERROR)
cmessage(DEBUG "\"${CONFIGAPP_LOCATION} ${ARG}\" produced no output and was expected to.")
set(${DIR_OUT} "" PARENT_SCOPE)
else()
cmessage(FATAL_ERROR "\"${CONFIGAPP_LOCATION} ${ARG}\" produced no output and was required to.")
endif()
else()
string(REGEX MATCH "${DELIM}\([^ ]+\)" PARSE_CONFIGAPP_RESPONSE_MATCH ${CONFIGAPP_RESPONSE_RAW})
if(NOT PARSE_CONFIGAPP_RESPONSE_MATCH)
if(FAILURE_IS_NOT_ERROR)
cmessage(DEBUG "Couldn't find ${DELIM} flag, found: \"${CONFIGAPP_RESPONSE_RAW}\"")
set(${CMAKE_MATCH_1} "")
else()
cmessage(FATAL_ERROR "Expected to be able to parse the result of ${CONFIGAPP} ${ARG} to a lib directory, but couldn't find a ${DELIM} flag, found: \"${CONFIGAPP_RESPONSE_RAW}\"")
endif()
endif()
set(${DIR_OUT} ${CMAKE_MATCH_1} PARENT_SCOPE)
endif()
else()
cmessage(FATAL_ERROR "[ERROR]: Failed to find dependency configuration application: \"${CONFIGAPP}\"")
endif()
endfunction()
#Uselike GETLIBDIR(gsl-config --libs GSL_LIB_DIR)
function(GETLIBDIR CONFIGAPP ARG LIBDIR_OUT)
if(ARGN)
set(FAILURE_IS_NOT_ERROR TRUE)
else()
set(FAILURE_IS_NOT_ERROR FALSE)
endif()
GETFIRSTMATCHINGDELMIMMEDDIR(
"-L"
${CONFIGAPP}
${ARG}
MATCHING_DIR
${FAILURE_IS_NOT_ERROR})
set(${LIBDIR_OUT} ${MATCHING_DIR} PARENT_SCOPE)
endfunction()
#Uselike GETINCDIR(gsl-config --cflags GSL_INC_DIR)
function(GETINCDIR CONFIGAPP ARG INCDIR_OUT)
if(ARGN)
set(FAILURE_IS_NOT_ERROR TRUE)
else()
set(FAILURE_IS_NOT_ERROR FALSE)
endif()
GETFIRSTMATCHINGDELMIMMEDDIR(
"-I"
${CONFIGAPP}
${ARG}
MATCHING_DIR
${FAILURE_IS_NOT_ERROR})
set(${INCDIR_OUT} ${MATCHING_DIR} PARENT_SCOPE)
endfunction()
function(GETALLMATCHINGDELMIMMEDDIR DELIM CONFIGAPP ARG LIST_OUT)
if(DELIM STREQUAL "")
cmessage(FATAL_ERROR "GETALLMATCHINGDELMIMMEDDIR Passed no delimiter. This is a build configuration bug in NUISANCE, please report to the developers.")
endif()
if(CONFIGAPP STREQUAL "")
cmessage(FATAL_ERROR "GETALLMATCHINGDELMIMMEDDIR Passed no configuration application. This is a build configuration bug in NUISANCE, please report to the developers.")
endif()
- #cmessage(WARNING "DELIM: ${DELIM}, CONFIGAPP: ${CONFIGAPP}, ARG: ${ARG}, LIST_OUT: ${LIST_OUT}")
-
SET(CONFIGAPP_LOCATION "CONFIGAPP_LOCATION-NOTFOUND")
find_program(CONFIGAPP_LOCATION ${CONFIGAPP})
if(NOT CONFIGAPP_LOCATION STREQUAL "CONFIGAPP_LOCATION-NOTFOUND")
execute_process (COMMAND ${CONFIGAPP_LOCATION}
${ARG} OUTPUT_VARIABLE CONFIGAPP_RESPONSE_RAW OUTPUT_STRIP_TRAILING_WHITESPACE)
string(REPLACE " " ";" CONFIGAPP_RESPONSE_LIST "${CONFIGAPP_RESPONSE_RAW}")
- # cmessage(WARNING "CONFIGAPP_RESPONSE_RAW: ${CONFIGAPP_RESPONSE_RAW}, CONFIGAPP_RESPONSE_LIST: ${CONFIGAPP_RESPONSE_LIST}")
-
set(LIST_BUILD)
foreach(I ${CONFIGAPP_RESPONSE_LIST})
if(I)
string(REGEX MATCH "^${DELIM}" WASMATCHED ${I})
if(WASMATCHED)
string(REPLACE "${DELIM}" "" I_STRIPPED "${I}")
LIST(APPEND LIST_BUILD ${I_STRIPPED})
endif()
endif()
endforeach()
set(${LIST_OUT} ${LIST_BUILD} PARENT_SCOPE)
else()
cmessage(FATAL_ERROR "[ERROR]: Failed to find dependency configuration application: \"${CONFIGAPP}\"")
endif()
endfunction()
#Uselike GETLIBDIRS(gsl-config --libs GSL_LIB_DIR)
function(GETLIBDIRS CONFIGAPP ARG LIBDIR_OUT)
GETALLMATCHINGDELMIMMEDDIR(
"-L"
${CONFIGAPP}
${ARG}
MATCHING_DIR)
set(${LIBDIR_OUT} ${MATCHING_DIR} PARENT_SCOPE)
endfunction()
#Uselike GETINCDIRS(gsl-config --cflags GSL_INC_DIR)
function(GETINCDIRS CONFIGAPP ARG INCDIR_OUT)
GETALLMATCHINGDELMIMMEDDIR(
"-I"
${CONFIGAPP}
${ARG}
MATCHING_DIR)
set(${INCDIR_OUT} ${MATCHING_DIR} PARENT_SCOPE)
endfunction()
#Uselike GETLIBS(gsl-config --libs GSL_LIB_DIR)
function(GETLIBS CONFIGAPP ARG LIBLIST_OUT)
- #cmessage(WARNING "LIBLIST_OUT: ${LIBLIST_OUT}")
GETALLMATCHINGDELMIMMEDDIR(
"-l"
${CONFIGAPP}
${ARG}
MATCHING_ITEMS)
set(${LIBLIST_OUT} ${MATCHING_ITEMS} PARENT_SCOPE)
- #cmessage(WARNING "LIBLIST_OUT: ${LIBLIST_OUT}: ${${LIBLIST_OUT}}")
+endfunction()
+
+function(BuildFlagString OUTPUTLIST DELIMITER)
+ LIST(APPEND INPUT_LIST ${ARGN})
+ if(NOT DELIMITER STREQUAL "")
+ string(STRIP ${DELIMITER} DELIMITER)
+ endif()
+ if(NOT "${INPUT_LIST}" STREQUAL "")
+ string(REPLACE ";" " ${DELIMITER}" LIST_STR "${DELIMITER}${INPUT_LIST}")
+ string(STRIP "${LIST_STR}" LIST_STR)
+ set(${OUTPUTLIST} ${LIST_STR} PARENT_SCOPE)
+ else()
+ set(${OUTPUTLIST} PARENT_SCOPE)
+ endif()
+endfunction()
+
+#Here we want to look for linker flags in a list of libraries and escape them
+function(BuildLibraryFlagString OUTPUTLIST)
+ set(INPUT_LIST ${ARGN})
+ BuildFlagString(NEED_SCRUB_LINKOPTS "-l" ${INPUT_LIST})
+ string(REPLACE "-l-" "-" LIST_STR ${NEED_SCRUB_LINKOPTS})
+ set(${OUTPUTLIST} ${LIST_STR} PARENT_SCOPE)
+endfunction()
+
+function(CatStringsIfNotEmpty OUTPUT_STRING)
+ set(BUILD_STR)
+ foreach(I ${ARGN})
+ if(NOT I STREQUAL "")
+ set(BUILD_STR "${BUILD_STR} ${I}")
+ endif()
+ endforeach()
+ if(NOT BUILD_STR STREQUAL "")
+ string(STRIP "${BUILD_STR}" BUILD_STR)
+ endif()
+ set(${OUTPUT_STRING} ${BUILD_STR} PARENT_SCOPE)
+endfunction()
+
+function(PrefixList OUTPUT_LIST PREFIX)
+ set(${OUTPUT_LIST})
+ foreach(I ${ARGN})
+ if(NOT I STREQUAL "")
+ LIST(APPEND ${OUTPUT_LIST} "${PREFIX}${I}")
+ endif()
+ endforeach()
+ set(${OUTPUT_LIST} ${${OUTPUT_LIST}} PARENT_SCOPE)
endfunction()
diff --git a/src/FitBase/Measurement2D.cxx b/src/FitBase/Measurement2D.cxx
index 37316bc..0cc89f3 100644
--- a/src/FitBase/Measurement2D.cxx
+++ b/src/FitBase/Measurement2D.cxx
@@ -1,1982 +1,2018 @@
// 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#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;
+ fResidualHist = 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;
+
+ delete fResidualHist;
}
//********************************************************************
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
NUIS_LOG(SAM, "Finalising Sample Settings: " << fName);
if ((fSettings.GetS("originalname").find("Evt") != std::string::npos)) {
fIsRawEvents = true;
NUIS_LOG(SAM, "Found event rate measurement but using poisson likelihoods.");
}
if (fSettings.GetS("originalname").find("Enu") != std::string::npos) {
fIsEnu1D = true;
NUIS_LOG(SAM, "::" << fName << "::");
NUIS_LOG(SAM, "Found XSec Enu measurement, applying flux integrated scaling, "
<< "not flux averaged!");
}
if (fIsEnu1D && fIsRawEvents) {
NUIS_ERR(FTL, "Found 2D Enu XSec distribution AND fIsRawEvents, is this "
"really correct?!");
NUIS_ERR(FTL, "Check experiment constructor for " << fName
<< " and correct this!");
NUIS_ABORT("I live in " << __FILE__ << ":" << __LINE__);
}
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) {
NUIS_ERR(FTL, "Norm error for class " << fName << " is 0.0!");
NUIS_ABORT("If you want to use it please add fNormError=VAL");
}
}
}
void Measurement2D::CreateDataHistogram(int dimx, double *binx, int dimy,
double *biny) {
if (fDataHist)
delete fDataHist;
NUIS_LOG(SAM, "Creating Data Histogram dim : " << dimx << " " << dimy);
fDataHist = new TH2D((fSettings.GetName() + "_data").c_str(),
(fSettings.GetFullTitles()).c_str(), dimx - 1, binx,
dimy - 1, biny);
}
void Measurement2D::SetDataFromTextFile(std::string data, std::string binx,
std::string biny) {
// Get the data hist
fDataHist = PlotUtils::GetTH2DFromTextFile(data, binx, biny);
// Set the name properly
fDataHist->SetName((fSettings.GetName() + "_data").c_str());
fDataHist->SetTitle(fSettings.GetFullTitles().c_str());
}
void Measurement2D::SetDataFromRootFile(std::string datfile,
std::string histname) {
NUIS_LOG(SAM, "Reading data from root file: " << datfile << ";" << histname);
fDataHist = PlotUtils::GetTH2DFromRootFile(datfile, histname);
fDataHist->SetNameTitle((fSettings.GetName() + "_data").c_str(),
(fSettings.GetFullTitles()).c_str());
}
void Measurement2D::SetDataValuesFromTextFile(std::string datfile, TH2D *hist) {
NUIS_LOG(SAM, "Setting data values from text file");
if (!hist)
hist = fDataHist;
// Read TH2D From textfile
TH2D *valhist = (TH2D *)hist->Clone();
valhist->Reset();
PlotUtils::Set2DHistFromText(datfile, valhist, 1.0, true);
NUIS_LOG(SAM, " -> Filling values from read hist.");
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));
}
}
NUIS_LOG(SAM, " --> Done");
}
void Measurement2D::SetDataErrorsFromTextFile(std::string datfile, TH2D *hist) {
NUIS_LOG(SAM, "Setting data errors from text file");
if (!hist)
hist = fDataHist;
// Read TH2D From textfile
TH2D *valhist = (TH2D *)hist->Clone();
valhist->Reset();
PlotUtils::Set2DHistFromText(datfile, valhist, 1.0);
// Fill Errors
NUIS_LOG(SAM, " -> Filling errors from read hist.");
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));
}
}
NUIS_LOG(SAM, " --> Done");
}
void Measurement2D::SetMapValuesFromText(std::string dataFile) {
TH2D *hist = fDataHist;
std::vector<double> edgex;
std::vector<double> 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]);
NUIS_LOG(SAM, "Reading map from: " << dataFile);
PlotUtils::Set2DHistFromText(dataFile, fMapHist, 1.0);
}
//********************************************************************
void Measurement2D::SetPoissonErrors() {
//********************************************************************
if (!fDataHist) {
NUIS_ERR(FTL, "Need a data hist to setup possion errors! ");
NUIS_ABORT("Setup Data First!");
}
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) {
NUIS_LOG(SAM, "Setting diagonal covariance for: " << data->GetName());
fFullCovar = StatUtils::MakeDiagonalCovarMatrix(data);
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
} else {
NUIS_ABORT("No data input provided to set diagonal covar from!");
}
// 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();
}
NUIS_LOG(SAM, "Reading covariance from text file: " << covfile << " " << dim);
fFullCovar = StatUtils::GetCovarFromTextFile(covfile, dim);
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
//********************************************************************
void Measurement2D::SetCovarFromRootFile(std::string covfile,
std::string histname) {
//********************************************************************
NUIS_LOG(SAM,
"Reading covariance from text file: " << covfile << ";" << histname);
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();
}
NUIS_LOG(SAM, "Reading inverted covariance from text file: " << covfile);
covar = StatUtils::GetCovarFromTextFile(covfile, dim);
fFullCovar = StatUtils::GetInvert(covar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
//********************************************************************
void Measurement2D::SetCovarInvertFromRootFile(std::string covfile,
std::string histname) {
//********************************************************************
NUIS_LOG(SAM, "Reading inverted covariance from text file: " << covfile << ";"
<< histname);
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();
NUIS_LOG(SAM,
"Reading data correlations from text file: " << covfile << ";" << dim);
TMatrixDSym *correlation = StatUtils::GetCovarFromTextFile(covfile, dim);
if (!fDataHist) {
NUIS_ABORT("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");
}
// 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) {
//********************************************************************
NUIS_LOG(SAM, "Reading data correlations from text file: " << covfile << ";"
<< histname);
TMatrixDSym *correlation = StatUtils::GetCovarFromRootFile(covfile, histname);
if (!fDataHist) {
NUIS_ABORT("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");
}
// 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();
}
NUIS_LOG(SAM, "Reading cholesky from text file: " << covfile << " " << dim);
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) {
//********************************************************************
NUIS_LOG(SAM, "Reading cholesky decomp from root file: " << covfile << ";"
<< histname);
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;
NUIS_LOG(SAM, "Reading bin mask from file: " << maskfile);
// 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(), std::ifstream::in);
if (!mask.is_open()) {
NUIS_LOG(FTL, " Cannot find mask file.");
throw;
}
while (std::getline(mask >> std::ws, line, '\n')) {
std::vector<int> entries = GeneralUtils::ParseToInt(line, " ");
// Skip lines with poorly formatted lines
if (entries.size() < 2) {
NUIS_LOG(WRN, "Measurement2D::SetBinMask(), couldn't parse line: " << line);
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() {
//********************************************************************
NUIS_LOG(SAM, "Finalising Measurement: " << fName);
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) {
NUIS_ABORT("No data has been setup inside " << fName << " constructor!");
}
// 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) {
NUIS_ERR(FTL, "I found a negative fScaleFactor in " << __FILE__ << ":"
<< __LINE__);
NUIS_ERR(FTL, "fScaleFactor = " << fScaleFactor);
NUIS_ABORT("EXITING");
}
// 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");
}
if (!fMapHist)
fMapHist = StatUtils::GenerateMap(fDataHist);
}
//********************************************************************
void Measurement2D::SetFitOptions(std::string opt) {
//********************************************************************
// Do nothing if default given
if (opt == "DEFAULT")
return;
// CHECK Conflicting Fit Options
std::vector<std::string> fit_option_allow =
GeneralUtils::ParseToStr(fAllowedTypes, "/");
for (UInt_t i = 0; i < fit_option_allow.size(); i++) {
std::vector<std::string> 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) {
NUIS_ABORT("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.");
}
}
}
// Check all options are allowed
std::vector<std::string> 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) {
NUIS_ERR(FTL, "ERROR: Fit Option '"
<< fit_options_input.at(i)
<< "' Provided is not allowed for this measurement.");
NUIS_ERR(FTL, "Fit Options should be provided as a '/' seperated list "
"(e.g. FREE/DIAG/NORM)");
NUIS_ABORT("Available options for " << fName << " are '" << fAllowedTypes
<< "'");
}
}
// 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;
NUIS_ERR(FTL, "No other LIKELIHOODS properly supported!");
NUIS_ABORT("Try to use a chi2!");
} 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(); ++xBin) {
for (int yBin = 0; yBin < fDataHist->GetNbinsY(); ++yBin) {
if (fDataHist->GetBinError(xBin + 1, yBin + 1) != 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);
+ fMaskHist, fResidualHist);
+ if (fChi2LessBinHist) {
+ TH2I *binmask = new TH2I("mask", "", fDataHist->GetNbinsX(), 0,
+ fDataHist->GetNbinsX(), fDataHist->GetNbinsY(),
+ 0, fDataHist->GetNbinsY());
+ binmask->SetDirectory(NULL);
+ for (int xi = 0; xi < fDataHist->GetNbinsX(); ++xi) {
+ for (int yi = 0; yi < fDataHist->GetNbinsY(); ++yi) {
+ binmask->Reset();
+ binmask->SetBinContent(xi + 1, yi + 1, 1);
+ fChi2LessBinHist->SetBinContent(
+ xi + 1, yi + 1,
+ StatUtils::GetChi2FromCov(fDataHist, fMCHist, covar, fMapHist,
+ binmask));
+ }
+ }
+ delete binmask;
+ }
}
}
// Add a normal penalty term
if (fAddNormPen) {
chi2 +=
(1 - (fCurrentNorm)) * (1 - (fCurrentNorm)) / (fNormError * fNormError);
NUIS_LOG(REC, "Norm penalty = " << (1 - (fCurrentNorm)) * (1 - (fCurrentNorm)) /
(fNormError * fNormError));
}
// 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;
NUIS_LOG(SAM, "Setting fake data from : " << fFakeDataInput);
// 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();
}
+ if (fIsChi2 && !fIsDiag) {
+ fResidualHist = (TH2D *)fMCHist->Clone((fName + "_RESIDUAL").c_str());
+ fResidualHist->GetYaxis()->SetTitle("#Delta#chi^{2}");
+ fResidualHist->Reset();
+
+ fChi2LessBinHist = (TH2D *)fMCHist->Clone((fName + "_Chi2NMinusOne").c_str());
+ fChi2LessBinHist->GetYaxis()->SetTitle("Total #chi^{2} without bin_{i}");
+ fChi2LessBinHist->Reset();
+
+ (void)GetLikelihood();
+
+ fResidualHist->Write();
+ fChi2LessBinHist->Write();
+ }
+
// // 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;
// }
/// 2D VERSION
// If null pointer return
if (!fMCHist and !fDataHist) {
NUIS_LOG(SAM, fName << "Incomplete histogram set!");
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 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 drawSliceMC = (drawOpt.find("CANVSLICEMC") != std::string::npos);
bool drawWeighted =
(drawOpt.find("WEIGHTS") != std::string::npos && fMCWeighted);
if (FitPar::Config().GetParB("EventManager")) {
drawFlux = false;
drawEvents = false;
}
// Save standard plots
if (drawData) {
GetDataList().at(0)->Write();
// Generate a simple map
if (!fMapHist)
fMapHist = StatUtils::GenerateMap(fDataHist);
// Convert to 1D Lists
TH1D *data_1D = StatUtils::MapToTH1D(fDataHist, fMapHist);
data_1D->Write();
delete data_1D;
}
if (drawNormal) {
GetMCList().at(0)->Write();
if (!fMapHist)
fMapHist = StatUtils::GenerateMap(fDataHist);
TH1D *mc_1D = StatUtils::MapToTH1D(fMCHist, fMapHist);
mc_1D->SetLineColor(kRed);
mc_1D->Write();
delete mc_1D;
}
// Write Weighted Histogram
if (drawWeighted)
fMCWeighted->Write();
if (drawCov) {
TH2D(*fFullCovar).Write((fName + "_COV").c_str());
}
if (drawOpt.find("INVCOV") != std::string::npos) {
TH2D(*covar).Write((fName + "_INVCOV").c_str());
}
// 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
TH1I *mask_1D = StatUtils::MapToMask(fMaskHist, fMapHist);
if (mask_1D) {
mask_1D->Write();
TMatrixDSym *calc_cov =
StatUtils::ApplyInvertedMatrixMasking(covar, mask_1D);
TH1D *data_1D = StatUtils::MapToTH1D(fDataHist, fMapHist);
TH1D *mc_1D = StatUtils::MapToTH1D(fMCHist, fMapHist);
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;
}
}
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 || fIsProjFitY || 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);
// This is not the correct way of doing it
// 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 (drawSliceMC) {
TCanvas *c1 = new TCanvas((fName + "_MC_CANV_Y").c_str(),
(fName + "_MC_CANV_Y").c_str(), 1024, 1024);
c1->Divide(2, int(fDataHist->GetNbinsY() / 3. + 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.0, 0.0, 1.0, 1.0);
leg->AddEntry(fDataHist, (fName + " Data").c_str(), "lep");
leg->AddEntry(fMCHist, (fName + " MC").c_str(), "l");
leg->AddEntry(mcShape, (fName + " Shape").c_str(), "l");
leg->SetLineColor(0);
leg->SetLineStyle(0);
leg->SetFillColor(0);
leg->SetLineStyle(0);
leg->Draw("SAME");
// Make Y slices
for (int i = 1; i < fDataHist->GetNbinsY() + 1; i++) {
c1->cd(i + 1);
TH1D *fDataHist_SliceY = PlotUtils::GetSliceY(fDataHist, i);
fDataHist_SliceY->Draw("E1");
TH1D *fMCHist_SliceY = PlotUtils::GetSliceY(fMCHist, i);
fMCHist_SliceY->Draw("SAME");
TH1D *mcShape_SliceY = PlotUtils::GetSliceY(mcShape, i);
mcShape_SliceY->Draw("SAME");
mcShape_SliceY->SetLineStyle(mcShape->GetLineStyle());
}
c1->Write();
delete c1;
delete leg;
TCanvas *c2 = new TCanvas((fName + "_MC_CANV_X").c_str(),
(fName + "_MC_CANV_X").c_str(), 1024, 1024);
c2->Divide(2, int(fDataHist->GetNbinsX() / 3. + 1));
mcShape = (TH2D *)fMCHist->Clone((fName + "_MC_SHAPE").c_str());
shapeScale = fDataHist->Integral("width") / fMCHist->Integral("width");
mcShape->Scale(shapeScale);
mcShape->SetLineStyle(7);
c2->cd(1);
TLegend *leg2 = new TLegend(0.0, 0.0, 1.0, 1.0);
leg2->AddEntry(fDataHist, (fName + " Data").c_str(), "lep");
leg2->AddEntry(fMCHist, (fName + " MC").c_str(), "l");
leg2->AddEntry(mcShape, (fName + " Shape").c_str(), "l");
leg2->SetLineColor(0);
leg2->SetLineStyle(0);
leg2->SetFillColor(0);
leg2->SetLineStyle(0);
leg2->Draw("SAME");
// Make Y slices
for (int i = 1; i < fDataHist->GetNbinsX() + 1; i++) {
c2->cd(i + 1);
TH1D *fDataHist_SliceX = PlotUtils::GetSliceX(fDataHist, i);
fDataHist_SliceX->Draw("E1");
TH1D *fMCHist_SliceX = PlotUtils::GetSliceX(fMCHist, i);
fMCHist_SliceX->Draw("SAME");
TH1D *mcShape_SliceX = PlotUtils::GetSliceX(mcShape, i);
mcShape_SliceX->Draw("SAME");
mcShape_SliceX->SetLineStyle(mcShape->GetLineStyle());
}
c2->Write();
delete c2;
delete leg2;
}
// Write Extra Histograms
AutoWriteExtraTH1();
WriteExtraHistograms();
// Returning
NUIS_LOG(SAM, "Written Histograms: " << fName);
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;
NUIS_LOG(SAM, "Found event rate measurement but fIsRawEvents == false!");
NUIS_LOG(SAM, "Overriding this and setting fIsRawEvents == true!");
}
fIsEnu = false;
if ((fName.find("XSec") != std::string::npos) &&
(fName.find("Enu") != std::string::npos)) {
fIsEnu = true;
NUIS_LOG(SAM, "::" << fName << "::");
NUIS_LOG(SAM, "Found XSec Enu measurement, applying flux integrated scaling, "
"not flux averaged!");
if (FitPar::Config().GetParB("EventManager")) {
NUIS_ERR(FTL, "Enu Measurements do not yet work with the Event Manager!");
NUIS_ERR(FTL, "If you want decent flux unfolded results please run in "
"series mode (-q EventManager=0)");
sleep(2);
throw;
}
}
if (fIsEnu && fIsRawEvents) {
NUIS_ERR(FTL, "Found 1D Enu XSec distribution AND fIsRawEvents, is this "
"really correct?!");
NUIS_ERR(FTL, "Check experiment constructor for " << fName
<< " and correct this!");
NUIS_ABORT("I live in " << __FILE__ << ":" << __LINE__);
}
// 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) {
NUIS_ERR(FTL, "Error! " << dataFile << " is not a .root file");
NUIS_ERR(FTL, "Currently only .root file reading is supported (MiniBooNE "
"CC1pi+ 2D), but implementing .txt should be dirt easy");
NUIS_ABORT("See me at " << __FILE__ << ":" << __LINE__);
} 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(), std::ifstream::in);
fDataHist = new TH2D((fName + "_data").c_str(), (fName + fPlotTitles).c_str(),
fNDataPointsX - 1, fXBins, fNDataPointsY - 1, fYBins);
if (data.is_open()) {
NUIS_LOG(SAM, "Reading data from: " << dataFile.c_str());
}
while (std::getline(data >> std::ws, line, '\n')) {
int xBin = 0;
// Loop over entries and insert them into the histogram
std::vector<double> entries = GeneralUtils::ParseToDbl(line, " ");
for (std::vector<double>::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(), std::ifstream::in);
if (error.is_open()) {
NUIS_LOG(SAM, "Reading errors from: " << errorFile.c_str());
}
while (std::getline(error >> std::ws, line, '\n')) {
int xBin = 0;
// Loop over entries and insert them into the histogram
std::vector<double> entries = GeneralUtils::ParseToDbl(line, " ");
for (std::vector<double>::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);
NUIS_LOG(SAM, "Reading data from: " << dataFile);
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 *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());
// 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 {
NUIS_ERR(WRN, " Incorrect thrown_covariance option in parameters.");
}
// 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(), std::ifstream::in);
this->covar = new TMatrixDSym(dim);
fFullCovar = new TMatrixDSym(dim);
if (covar.is_open()) {
NUIS_LOG(SAM, "Reading covariance matrix from file: " << covarFile);
}
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<double> entries = GeneralUtils::ParseToDbl(line, " ");
for (std::vector<double>::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(), std::ifstream::in);
TMatrixD *newcov = new TMatrixD(dim, dim);
if (covarread.is_open()) {
NUIS_LOG(SAM, "Reading covariance matrix from file: " << covarFile);
}
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<double> entries = GeneralUtils::ParseToDbl(line, " ");
for (std::vector<double>::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/FitBase/Measurement2D.h b/src/FitBase/Measurement2D.h
index 2ab738c..00d32f8 100644
--- a/src/FitBase/Measurement2D.h
+++ b/src/FitBase/Measurement2D.h
@@ -1,641 +1,644 @@
// 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#ifndef MEASUREMENT_2D_HXX_SEEN
#define MEASUREMENT_2D_HXX_SEEN
/*!
* \addtogroup FitBase
* @{
*/
#include <math.h>
#include <stdlib.h>
#include <deque>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <numeric>
#include <sstream>
#include <string>
// ROOT includes
#include <TArrayF.h>
#include <TDecompSVD.h>
#include <TGraph.h>
#include <TGraphErrors.h>
#include <TH1D.h>
#include <TH2D.h>
#include <TMatrixDSym.h>
#include <TROOT.h>
#include <TSystem.h>
// External data fit includes
#include "FitEvent.h"
#include "FitUtils.h"
#include "MeasurementBase.h"
#include "PlotUtils.h"
#include "SignalDef.h"
#include "StatUtils.h"
#include "MeasurementVariableBox2D.h"
//********************************************************************
//! 2D Measurement base class. Histogram handling is done in this base layer.
class Measurement2D : public MeasurementBase {
//********************************************************************
public:
/*
Constructor/Deconstuctor
*/
//! Default Constructor
Measurement2D();
//! Default Destructor
virtual ~Measurement2D();
/*
Setup Functions
*/
/// \brief Setup all configs once initialised
///
/// Should be called after all configs have been setup inside fSettings container.
/// Handles the processing of inputs and setting up of types.
/// Replaces the old 'SetupMeasurement' function.
void FinaliseSampleSettings();
/// \brief Creates the 2D data distribution given the binning provided.
virtual void CreateDataHistogram(int dimx, double* binx, int dimy, double* biny);
/// \brief Set Data Histogram from a list of contents in a text file
///
/// Assumes the format: \n
/// x_low_1 y_low_1 cont_11 err_11 \n
/// x_low_1 y_low_2 cont_12 err_12 \n
/// x_low_2 y_low_1 cont_21 err_21 \n
/// x_low_2 y_low_2 cont_22 err_22 \n
/// x_low_2 y_low_3 cont_23 err_23 \n
/// x_low_3 y_low_2 cont_32 err_32 \n
virtual void SetDataFromTextFile(std::string data, std::string binx, std::string biny);
/// \brief Set Data Histogram from a TH2D in a file
///
/// - datfile = Full path to data file
/// - histname = Name of histogram
///
/// If histname not given it assumes that datfile
/// is in the format: \n
/// 'file.root;histname'
virtual void SetDataFromRootFile(std::string datfile, std::string histname="");
/// \brief Set data values from a 2D array in text file
///
/// \warning requires DATA HISTOGRAM TO BE SET FIRST
///
/// Assumes form: \n
/// cont_11 cont_12 ... cont_1N \n
/// cont_21 cont_22 ... cont_2N \n
/// ... ... ... ... \n
/// cont_N1 cont_N2 ... cont_NN \n
virtual void SetDataValuesFromTextFile(std::string datfile, TH2D* hist = NULL);
/// \brief Set data errors from a 2D array in text file
///
/// \warning requires DATA HISTOGRAM TO BE SET FIRST
///
/// Assumes form: \n
/// errs_11 errs_12 ... errs_1N \n
/// errs_21 errs_22 ... errs_2N \n
/// ... ... ... ... \n
/// errs_N1 errs_N2 ... errs_NN \n
virtual void SetDataErrorsFromTextFile(std::string datfile, TH2D* hist = NULL);
/// \brief Set data bin errors to sqrt(entries)
///
/// \warning REQUIRES DATA HISTOGRAM TO BE SET FIRST
///
/// Sets the data errors as the sqrt of the bin contents
/// Should be use for counting experiments
virtual void SetPoissonErrors();
/// \brief Make diagonal covariance from data
///
/// \warning If no histogram passed, data must be setup first!
/// Setup the covariance inputs by taking the data histogram
/// errors and setting up a diagonal covariance matrix.
///
/// If no data is supplied, fDataHist is used if already set.
virtual void SetCovarFromDiagonal(TH2D* data = NULL);
/// \brief Read the data covariance from a text file.
///
/// Inputfile should have the format: \n
/// covariance_11 covariance_12 covariance_13 ... \n
/// covariance_21 covariance_22 covariance_23 ... \n
/// ... ... ... ... \n
///
/// If no dimensions are given, it is assumed from the number
/// entries in the first line of covfile.
virtual void SetCovarFromTextFile(std::string covfile, int dim = -1);
/// \brief Read the data covariance from a ROOT file.
///
/// - covfile specifies the full path to the file
/// - histname specifies the name of the covariance object. Both TMatrixDSym and TH2D are supported.
///
/// If no histogram name is given the inhistfile value
/// is automatically parsed with ; so that: \n
/// mycovfile.root;myhistname \n
/// will also work.
virtual void SetCovarFromRootFile(std::string covfile, std::string histname="");
/// \brief Read the inverted data covariance from a text file.
///
/// Inputfile should have similar format to that shown
/// in SetCovarFromTextFile.
///
/// If no dimensions are given, it is assumed from the number
/// entries in the first line of covfile.
virtual void SetCovarInvertFromTextFile(std::string covfile, int dim = -1);
/// \brief Read the inverted data covariance from a ROOT file.
///
/// Inputfile should have similar format to that shown
/// in SetCovarFromRootFile.
///
/// If no histogram name is given the inhistfile value
/// is automatically parsed with ; so that: \n
/// mycovfile.root;myhistname \n
/// will also work.
virtual void SetCovarInvertFromRootFile(std::string covfile, std::string histname="");
/// \brief Read the data correlations from a text file.
///
/// \warning REQUIRES DATA HISTOGRAM TO BE SET FIRST
///
/// Inputfile should have similar format to that shown
/// in SetCovarFromTextFile.
///
/// If no dimensions are given, it is assumed from the number
/// entries in the first line of covfile.
virtual void SetCorrelationFromTextFile(std::string covfile, int dim = -1);
/// \brief Read the data correlations from a ROOT file.
///
/// \warning REQUIRES DATA TO BE SET FIRST
///
/// Inputfile should have similar format to that shown
/// in SetCovarFromRootFile.
///
/// If no histogram name is given the inhistfile value
/// is automatically parsed with ; so that: \n
/// mycovfile.root;myhistname \n
/// will also work.
virtual void SetCorrelationFromRootFile(std::string covfile, std::string histname="");
/// \brief Read the cholesky decomposed covariance from a text file and turn it into a covariance
///
/// Inputfile should have similar format to that shown
/// in SetCovarFromTextFile.
///
/// If no dimensions are given, it is assumed from the number
/// entries in the first line of covfile.
virtual void SetCholDecompFromTextFile(std::string covfile, int dim = -1);
/// \brief Read the cholesky decomposed covariance from a ROOT file and turn it into a covariance
///
/// Inputfile should have similar format to that shown
/// in SetCovarFromRootFile.
///
/// If no histogram name is given the inhistfile value
/// is automatically parsed with ; so that: \n
/// mycovfile.root;myhistname \n
/// will also work.
virtual void SetCholDecompFromRootFile(std::string covfile, std::string histname="");
/// \brief Read the map values from a text file
///
/// \warning Requires DATA hist to be set beforehand.
- /// Format should be a 2D array of mappings.
+ /// Format should be a 2D array of mappings.
/// -1 indicates empty bins. \n
/// e.g.: \n
/// 1 2 3 4 5 \n
/// -1 6 7 8 9 \n
/// -1 -1 10 11 -1 \n
virtual void SetMapValuesFromText(std::string dataFile);
/// \brief Scale the data by some scale factor
virtual void ScaleData(double scale);
/// \brief Scale the data error bars by some scale factor
virtual void ScaleDataErrors(double scale);
/// \brief Scale the covariaince and its invert/decomp by some scale factor.
virtual void ScaleCovar(double scale);
/// \brief Setup a bin masking histogram and apply masking to data
///
/// \warning REQUIRES DATA HISTOGRAM TO BE SET FIRST
///
/// Reads in a list of bins in a text file to be masked. Format is: \n
/// bin_index_x_1 bin_index_y_1 1 \n
/// bin_index_x_2 bin_index_y_2 1 \n
/// bin_index_x_3 bin_index_y_3 1 \n
///
/// If 0 is given then a bin entry will NOT be masked. So for example: \n\n
/// 1 1 1 \n
/// 2 0 1 \n
/// 3 4 0 \n
/// 4 0 1 \n\n
/// Will mask only the (1,1), (2,0), and (4,0) bins.
///
/// Masking can be turned on by specifiying the MASK option when creating a sample.
/// When this is passed NUISANCE will look in the following locations for the mask file:
/// - FitPar::Config().GetParS(fName + ".mask")
/// - "data/masks/" + fName + ".mask";
virtual void SetBinMask(std::string maskfile);
/// \brief Set the current fit options from a string.
///
/// This is called twice for each sample, once to set the default
/// and once to set the current setting (if anything other than default given)
///
/// For this to work properly it requires the default and allowed types to be
/// set correctly. These should be specified as a string listing options.
///
/// To split up options so that NUISANCE can automatically detect ones that
/// are conflicting. Any options seperated with the '/' symbol are non conflicting
/// and can be given together, whereas any seperated with the ',' symbol cannot
/// be specified by the end user at the same time.
///
/// Default Type Examples:
/// - DIAG/FIX = Default option will be a diagonal covariance, with FIXED norm.
/// - MASK/SHAPE = Default option will be a masked hist, with SHAPE always on.
///
/// Allowed Type examples:
/// - 'FULL/DIAG/NORM/MASK' = Any of these options can be specified.
/// - 'FULL,FREE,SHAPE/MASK/NORM' = User can give either FULL, FREE, or SHAPE as on option.
/// MASK and NORM can also be included as options.
virtual void SetFitOptions(std::string opt);
/// \brief Final constructor setup
/// \warning Should be called right at the end of the constructor.
///
/// Contains a series of checks to ensure the data and inputs have been setup.
/// Also creates the MC histograms needed for fitting.
void FinaliseMeasurement();
/*
Reconfigure
*/
/// \brief Create a Measurement1D box
///
/// Creates a new 1D variable box containing just fXVar.
///
/// This box is the bare minimum required by the JointFCN when
/// running fast reconfigures during a routine.
/// If for some reason a sample needs extra variables to be saved then
/// it should override this function creating its own MeasurementVariableBox
/// that contains the extra variables.
virtual MeasurementVariableBox* CreateBox() {return new MeasurementVariableBox2D();};
/// \brief Reset all MC histograms
///
/// Resets all standard histograms and those registered to auto
/// process to zero.
///
/// If extra histograms are not included in auto processing, then they must be reset
/// by overriding this function and doing it manually if required.
virtual void ResetAll(void);
/// \brief Fill MC Histograms from XVar, YVar
///
/// Fill standard histograms using fXVar, fYVar, Weight read from the variable box.
///
/// WARNING : Any extra MC histograms need to be filled by overriding this function,
/// even if they have been set to auto process.
virtual void FillHistograms(void);
// \brief Convert event rates to final histogram
///
/// Apply standard scaling procedure to standard mc histograms to convert from
/// raw events to xsec prediction.
///
/// If any distributions have been set to auto process
/// that is done during this function call, and a differential xsec is assumed.
/// If that is not the case this function must be overriden.
virtual void ScaleEvents(void);
/// \brief Scale MC by a factor=1/norm
///
/// Apply a simple normalisation scaling if the option FREE or a norm_parameter
/// has been specified in the NUISANCE routine.
virtual void ApplyNormScale(double norm);
/*
Statistical Functions
*/
/// \brief Get Number of degrees of freedom
///
/// Returns the number bins inside the data histogram accounting for
/// any bin masking applied.
virtual int GetNDOF();
/// \brief Return Data/MC Likelihood at current state
///
/// Returns the likelihood of the data given the current MC prediction.
/// Diferent likelihoods definitions are used depending on the FitOptions.
virtual double GetLikelihood(void);
/*
Fake Data
*/
/// \brief Set the fake data values from either a file, or MC
///
/// - Setting from a file "path": \n
/// When reading from a file the full path must be given to a standard
/// nuisance output. The standard MC histogram should have a name that matches
/// this sample for it to be read in.
/// \n\n
/// - Setting from "MC": \n
/// If the MC option is given the current MC prediction is used as fake data.
virtual void SetFakeDataValues(std::string fakeOption);
/// \brief Reset fake data back to starting fake data
///
/// Reset the fake data back to original fake data (Reset back to before
/// ThrowCovariance was first called)
virtual void ResetFakeData(void);
/// \brief Reset fake data back to original data
///
/// Reset the data histogram back to the true original dataset for this sample
/// before any fake data was defined.
virtual void ResetData(void);
/// \brief Generate fake data by throwing the covariance.
///
/// Can be used on fake MC data or just the original dataset.
/// Call ResetFakeData or ResetData to return to values before the throw.
virtual void ThrowCovariance(void);
- /// \brief Throw the data by its assigned errors and assign this to MC
- ///
- /// Used when creating data toys by assign the MC to this thrown data
- /// so that the likelihood is calculated between data and thrown data
+ /// \brief Throw the data by its assigned errors and assign this to MC
+ ///
+ /// Used when creating data toys by assign the MC to this thrown data
+ /// so that the likelihood is calculated between data and thrown data
virtual void ThrowDataToy(void);
/*
Access Functions
*/
/// \brief Returns nicely formatted MC Histogram
///
/// Format options can also be given in the samplesettings:
/// - linecolor
/// - linestyle
/// - linewidth
/// - fillcolor
/// - fillstyle
///
/// So to have a sample line colored differently in the xml cardfile put: \n
/// <sample name="MiniBooNE_CCQE_XSec_1DQ2_nu" input="NEUT:input.root"
/// linecolor="2" linestyle="7" linewidth="2" />
virtual TH2D* GetMCHistogram(void);
/// \brief Returns nicely formatted data Histogram
///
/// Format options can also be given in the samplesettings:
/// - datacolor
/// - datastyle
/// - datawidth
///
/// So to have a sample data colored differently in the xml cardfile put: \n
/// <sample name="MiniBooNE_CCQE_XSec_1DQ2_nu" input="NEUT:input.root"
/// datacolor="2" datastyle="7" datawidth="2" />
virtual TH2D* GetDataHistogram(void);
/// \brief Returns a list of all MC histograms.
///
/// Override this if you have extra histograms that need to be
/// accessed outside of the Measurement1D class.
inline virtual std::vector<TH1*> GetMCList(void) {
return std::vector<TH1*>(1, GetMCHistogram());
}
/// \brief Returns a list of all Data histograms.
///
/// Override this if you have extra histograms that need to be
/// accessed outside of the Measurement1D class.
inline virtual std::vector<TH1*> GetDataList(void) {
return std::vector<TH1*>(1, GetDataHistogram());
}
/// \brief Returns a list of all Mask histograms.
///
/// Override this if you have extra histograms that need to be
/// accessed outside of the Measurement1D class.
inline virtual std::vector<TH1*> GetMaskList(void) {
return std::vector<TH1*>(1, fMaskHist);
};
/// \brief Returns a list of all Fine histograms.
///
/// Override this if you have extra histograms that need to be
/// accessed outside of the Measurement1D class.
inline virtual std::vector<TH1*> GetFineList(void) {
return std::vector<TH1*>(1, fMCFine);
};
/*
Write Functions
*/
/// \brief Save the current state to the current TFile directory \n
///
/// Data/MC are both saved by default.
/// A range of other histograms can be saved by setting the
/// config option 'drawopts'.
///
/// Possible options: \n
/// - FINE = Write Fine Histogram \n
/// - WEIGHTS = Write Weighted MC Histogram (before scaling) \n
/// - FLUX = Write Flux Histogram from MC Input \n
/// - EVT = Write Event Histogram from MC Input \n
/// - XSEC = Write XSec Histogram from MC Input \n
/// - MASK = Write Mask Histogram \n
/// - COV = Write Covariance Histogram \n
/// - INVCOV = Write Inverted Covariance Histogram \n
/// - DECMOP = Write Decomp. Covariance Histogram \n
/// - RESIDUAL= Write Resudial Histograms \n
/// - RATIO = Write Data/MC Ratio Histograms \n
/// - SHAPE = Write MC Shape Histograms norm. to Data \n
/// - CANVMC = Build MC Canvas Showing Data, MC, Shape \n
/// - MODES = Write PDG Stack \n
/// - CANVPDG = Build MC Canvas Showing Data, PDGStack \n
///
/// So to save a range of these in parameters/config.xml set: \n
/// <config drawopts='FINE/COV/SHAPE/RATIO' />
virtual void Write(std::string drawOpt);
//////// OLD FUNCTIONS ////////////
//! Intial setup of common measurement variables. Parse input files, types,
//! etc.
virtual void SetupMeasurement(std::string input, std::string type,
FitWeight* rw, std::string fkdt);
//! Setup the default mc Hist given a data histogram
virtual void SetupDefaultHist();
//! Set the data values and errors from two files
virtual void SetDataValues(std::string dataFile, double dataNorm,
std::string errorFile, double errorNorm);
virtual void SetDataValues(std::string dataFile, std::string TH2Dname);
//! Set the data values only from a text file
virtual void SetDataValuesFromText(std::string dataFile, double norm);
//! Read a covariance matrix from a file (Default name "covar" in file)
virtual void SetCovarMatrix(std::string covarFile);
//! Set the covariance matrix from a text file
virtual void SetCovarMatrixFromText(std::string covarFile, int dim);
//! Set the covariance matrix from a text file containing the cholesky
//! fDecomposition
virtual void SetCovarMatrixFromChol(std::string covarFile, int dim);
protected:
// The data histograms
TH2D* fDataHist; //!< default data histogram (use in chi2 calculations)
TH2D* fDataOrig; //!< histogram to store original data before throws.
TH2D* fDataTrue; //!< histogram to store true dataset
TH1D* fDataHist_X; //!< Projections onto X of the fDataHist
TH1D* fDataHist_Y; //!< Projections onto Y of the fDataHist
// The MC histograms
TH2D* fMCHist; //!< MC Histogram (used in chi2 calculations)
TH2D* fMCFine; //!< Finely binned MC Histogram
TH2D* fMCHist_PDG[61]; //!< MC Histograms for each interaction mode
TH1D* fMCHist_X; //!< Projections onto X of the fMCHist
TH1D* fMCHist_Y; //!< Projections onto Y of the fMCHist
TH2D* fMCWeighted; //!< Raw Event Weights
TH2D* fMCStat;
TH2I* fMaskHist; //!< mask histogram for the data
TH2I* fMapHist; //!< map histogram used to convert 2D to 1D distributions
+ TH2D *fResidualHist;
+ TH2D *fChi2LessBinHist;
+
bool fIsFakeData; //!< is current data actually fake
std::string fakeDataFile; //!< MC fake data input file
std::string fPlotTitles; //!< X and Y plot titles.
std::string fFitType;
std::string fDefaultTypes; //!< Default Fit Options
std::string fAllowedTypes; //!< Any allowed Fit Options
TMatrixDSym* covar; //!< inverted covariance matrix
TMatrixDSym* fFullCovar; //!< covariance matrix
TMatrixDSym* fDecomp; //!< fDecomposed covariance matrix
TMatrixDSym* fCorrel; //!< correlation matrix
double fCovDet; //!< covariance deteriminant
double fNormError; //!< Normalisation on the error on the data
double fLikelihood; //!< Likelihood value
Double_t* fXBins; //!< X Bin Edges
Double_t* fYBins; //!< Y Bin Edges
Int_t fNDataPointsX; //!< Number of X data points
Int_t fNDataPointsY; //!< NUmber of Y data points
// Fit specific flags
bool fIsShape; //!< Flag: Perform shape-only fit
bool fIsFree; //!< Flag: Perform normalisation free fit
bool fIsDiag; //!< Flag: Only use diagonal bin errors in stats
bool fIsMask; //!< Flag: Apply bin masking
bool fIsRawEvents; //!< Flag: Only event rates in histograms
bool fIsEnu; //!< Needs Enu Unfolding
bool fIsChi2SVD; //!< Flag: Chi2 SVD Method (DO NOT USE)
bool fAddNormPen; //!< Flag: Add normalisation penalty to fi
bool fIsProjFitX; //!< Flag: Use 1D projections onto X and Y to calculate the
//!Chi2 Method. If flagged X will be used to set the rate.
bool fIsProjFitY; //!< Flag: Use 1D projections onto X and Y to calculate the
//!Chi2 Method. If flagged Y will be used to set the rate.
bool fIsFix; //!< Flag: Fixed Histogram Norm
bool fIsFull; //!< Flag; Use Full Covar
bool fIsDifXSec; //!< Flag: Differential XSec
bool fIsEnu1D; //!< Flag: Flux Unfolded XSec
bool fIsChi2; //!< Flag; Use Chi2 over LL
TrueModeStack* fMCHist_Modes; ///< Optional True Mode Stack
TMatrixDSym* fCovar; ///< New FullCovar
TMatrixDSym* fInvert; ///< New covar
// Fake Data
std::string fFakeDataInput; ///< Input fake data file path
TFile* fFakeDataFile; ///< Input fake data file
// Arrays for data entries
Double_t* fDataValues; ///< REMOVE data bin contents
Double_t* fDataErrors; ///< REMOVE data bin errors
};
/*! @} */
#endif
diff --git a/src/MINERvA/MINERvA_CC0pi_XSec_2D_nu.cxx b/src/MINERvA/MINERvA_CC0pi_XSec_2D_nu.cxx
index 7990f68..fe5eb18 100644
--- a/src/MINERvA/MINERvA_CC0pi_XSec_2D_nu.cxx
+++ b/src/MINERvA/MINERvA_CC0pi_XSec_2D_nu.cxx
@@ -1,245 +1,284 @@
// 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
/*
Authors: Adrian Orea (v1 2017)
Clarence Wret (v2 2018)
*/
#include "MINERvA_CC0pi_XSec_2D_nu.h"
#include "MINERvA_SignalDef.h"
//********************************************************************
void MINERvA_CC0pi_XSec_2D_nu::SetupDataSettings() {
//********************************************************************
// Set Distribution
// See header file for enum and some descriptions
std::string name = fSettings.GetS("name");
// Have one distribution as of summer 2018
if (!name.compare("MINERvA_CC0pi_XSec_2Dptpz_nu"))
fDist = kPtPz;
// Define what files to use from the dist
std::string datafile = "";
std::string corrfile = "";
std::string titles = "";
std::string distdescript = "";
std::string histname = "";
switch (fDist) {
case (kPtPz):
datafile = "MINERvA/CC0pi_2D/cov_ptpl_2D_qelike.root";
corrfile = "MINERvA/CC0pi_2D/cov_ptpl_2D_qelike.root";
titles = "MINERvA CC0#pi #nu_{#mu} p_{t} p_{z};p_{z} (GeV);p_{t} "
"(GeV);d^{2}#sigma/dP_{t}dP_{z} (cm^{2}/GeV^{2}/nucleon)";
distdescript = "MINERvA_CC0pi_XSec_2Dptpz_nu sample";
histname = "pt_pl_cross_section";
break;
default:
NUIS_ABORT("Unknown Analysis Distribution : " << name);
}
fSettings.SetTitle(GeneralUtils::ParseToStr(titles, ";")[0]);
fSettings.SetXTitle(GeneralUtils::ParseToStr(titles, ";")[1]);
fSettings.SetYTitle(GeneralUtils::ParseToStr(titles, ";")[2]);
fSettings.SetZTitle(GeneralUtils::ParseToStr(titles, ";")[3]);
// Sample overview ---------------------------------------------------
std::string descrip = distdescript + "\n"
"Target: CH \n"
"Flux: MINERvA Low Energy FHC numu \n"
"Signal: CC-0pi \n";
fSettings.SetDescription(descrip);
// The input ROOT file
fSettings.SetDataInput(FitPar::GetDataBase() + datafile);
fSettings.SetCovarInput(FitPar::GetDataBase() + corrfile);
// Set the data file
SetDataValues(fSettings.GetDataInput(), histname);
}
//********************************************************************
MINERvA_CC0pi_XSec_2D_nu::MINERvA_CC0pi_XSec_2D_nu(nuiskey samplekey) {
//********************************************************************
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/FULL,DIAG/MASK", "FIX/FULL");
fSettings.SetEnuRange(0.0, 100.0);
fSettings.DefineAllowedTargets("C,H");
fSettings.DefineAllowedSpecies("numu");
SetupDataSettings();
FinaliseSampleSettings();
fScaleFactor =
(GetEventHistogram()->Integral("width") * 1E-38 / (fNEvents + 0.)) /
this->TotalIntegratedFlux();
TMatrixDSym *tempmat = StatUtils::GetCovarFromRootFile(
fSettings.GetCovarInput(), "TotalCovariance");
fFullCovar = tempmat;
// Decomposition is stable for entries that aren't E-xxx
double ScalingFactor = 1E38 * 1E38;
(*fFullCovar) *= ScalingFactor;
// Just check that the data error and covariance are the same
for (int i = 0; i < fFullCovar->GetNrows(); ++i) {
for (int j = 0; j < fFullCovar->GetNcols(); ++j) {
// Get the global bin
int xbin1, ybin1, zbin1;
fDataHist->GetBinXYZ(i, xbin1, ybin1, zbin1);
double xlo1 = fDataHist->GetXaxis()->GetBinLowEdge(xbin1);
double xhi1 = fDataHist->GetXaxis()->GetBinLowEdge(xbin1 + 1);
double ylo1 = fDataHist->GetYaxis()->GetBinLowEdge(ybin1);
double yhi1 = fDataHist->GetYaxis()->GetBinLowEdge(ybin1 + 1);
if (xlo1 < fDataHist->GetXaxis()->GetBinLowEdge(1) ||
ylo1 < fDataHist->GetYaxis()->GetBinLowEdge(1) ||
xhi1 > fDataHist->GetXaxis()->GetBinLowEdge(
fDataHist->GetXaxis()->GetNbins() + 1) ||
yhi1 > fDataHist->GetYaxis()->GetBinLowEdge(
fDataHist->GetYaxis()->GetNbins() + 1))
continue;
double data_error = fDataHist->GetBinError(xbin1, ybin1);
double cov_error = sqrt((*fFullCovar)(i, i) / ScalingFactor);
if (fabs(data_error - cov_error) > 1E-5) {
std::cerr << "Error on data is different to that of covariance"
<< std::endl;
NUIS_ERR(FTL, "Data error: " << data_error);
NUIS_ERR(FTL, "Cov error: " << cov_error);
NUIS_ERR(FTL, "Data/Cov: " << data_error / cov_error);
NUIS_ERR(FTL, "Data-Cov: " << data_error - cov_error);
NUIS_ERR(FTL, "For x: " << xlo1 << "-" << xhi1);
NUIS_ABORT("For y: " << ylo1 << "-" << yhi1);
}
}
}
// Now can make the inverted covariance
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
// Now scale back
(*fFullCovar) *= 1.0 / ScalingFactor;
(*covar) *= ScalingFactor;
(*fDecomp) *= ScalingFactor;
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
//********************************************************************
void MINERvA_CC0pi_XSec_2D_nu::FillEventVariables(FitEvent *event) {
//********************************************************************
// Checking to see if there is a Muon
if (event->NumFSParticle(13) == 0)
return;
// Get the muon kinematics
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
Double_t px = Pmu.X() / 1000;
Double_t py = Pmu.Y() / 1000;
Double_t pt = sqrt(px * px + py * py);
// y-axis is transverse momentum for both measurements
fYVar = pt;
// Now we set the x-axis
switch (fDist) {
case kPtPz: {
// Don't want to assume the event generators all have neutrino coming along
// z pz is muon momentum projected onto the neutrino direction
Double_t pz = Pmu.Vect().Dot(Pnu.Vect() * (1.0 / Pnu.Vect().Mag())) / 1000.;
// Set Hist Variables
fXVar = pz;
break;
}
default:
NUIS_ABORT("DIST NOT FOUND : " << fDist);
break;
}
};
//********************************************************************
bool MINERvA_CC0pi_XSec_2D_nu::isSignal(FitEvent *event) {
//********************************************************************
return SignalDef::isCC0pi_MINERvAPTPZ(event, 14, EnuMin, EnuMax);
};
//********************************************************************
// Custom likelihood calculator because binning of covariance matrix
double MINERvA_CC0pi_XSec_2D_nu::GetLikelihood() {
//********************************************************************
// The calculated chi2
double chi2 = 0.0;
// Support shape comparisons
double scaleF = fDataHist->Integral() / fMCHist->Integral();
if (fIsShape) {
fMCHist->Scale(scaleF);
fMCFine->Scale(scaleF);
}
// Even though this chi2 calculation looks ugly it is _EXACTLY_ what MINERvA
// used for their measurement Can be prettified in due time but for now keep
int nbinsx = fMCHist->GetNbinsX();
int nbinsy = fMCHist->GetNbinsY();
Int_t Nbins = nbinsx * nbinsy;
// Loop over the covariance matrix bins
for (int i = 0; i < Nbins; ++i) {
- int xbin = i % nbinsx + 1;
- int ybin = i / nbinsx + 1;
+ int xbin = (i % nbinsx) + 1;
+ int ybin = (i / nbinsx) + 1;
double datax = fDataHist->GetBinContent(xbin, ybin);
double mcx = fMCHist->GetBinContent(xbin, ybin);
+ double chi2_bin = 0;
for (int j = 0; j < Nbins; ++j) {
- int xbin2 = j % nbinsx + 1;
- int ybin2 = j / nbinsx + 1;
+ int xbin2 = (j % nbinsx) + 1;
+ int ybin2 = (j / nbinsx) + 1;
double datay = fDataHist->GetBinContent(xbin2, ybin2);
double mcy = fMCHist->GetBinContent(xbin2, ybin2);
double chi2_xy = (datax - mcx) * (*covar)(i, j) * (datay - mcy);
- chi2 += chi2_xy;
+ chi2_bin += chi2_xy;
+ }
+ if (fResidualHist) {
+ fResidualHist->SetBinContent(xbin, ybin, chi2_bin);
+ }
+ chi2 += chi2_bin;
+ }
+
+ if (fChi2LessBinHist) {
+ for (int igbin = 0; igbin < Nbins; ++igbin) {
+ int igxbin = (igbin % nbinsx) + 1;
+ int igybin = (igbin / nbinsx) + 1;
+ double tchi2 = 0;
+ for (int i = 0; i < Nbins; ++i) {
+ int xbin = (i % nbinsx) + 1;
+ int ybin = (i / nbinsx) + 1;
+ if ((xbin == igxbin) && (ybin == igybin)) {
+ continue;
+ }
+ double datax = fDataHist->GetBinContent(xbin, ybin);
+ double mcx = fMCHist->GetBinContent(xbin, ybin);
+ double chi2_bin = 0;
+ for (int j = 0; j < Nbins; ++j) {
+ int xbin2 = (j % nbinsx) + 1;
+ int ybin2 = (j / nbinsx) + 1;
+ if ((xbin2 == igxbin) && (ybin2 == igybin)) {
+ continue;
+ }
+ double datay = fDataHist->GetBinContent(xbin2, ybin2);
+ double mcy = fMCHist->GetBinContent(xbin2, ybin2);
+
+ double chi2_xy = (datax - mcx) * (*covar)(i, j) * (datay - mcy);
+ chi2_bin += chi2_xy;
+ }
+ tchi2 += chi2_bin;
+ }
+
+ fChi2LessBinHist->SetBinContent(igxbin, igybin, tchi2);
}
}
// Normalisation penalty term if included
if (fAddNormPen) {
chi2 +=
(1 - (fCurrentNorm)) * (1 - (fCurrentNorm)) / (fNormError * fNormError);
- NUIS_LOG(REC, "Norm penalty = " << (1 - (fCurrentNorm)) * (1 - (fCurrentNorm)) /
- (fNormError * fNormError));
+ NUIS_LOG(REC, "Norm penalty = " << (1 - (fCurrentNorm)) *
+ (1 - (fCurrentNorm)) /
+ (fNormError * fNormError));
}
// 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;
};
diff --git a/src/Reweight/GENIEWeightEngine.cxx b/src/Reweight/GENIEWeightEngine.cxx
index 97c060e..827a00f 100644
--- a/src/Reweight/GENIEWeightEngine.cxx
+++ b/src/Reweight/GENIEWeightEngine.cxx
@@ -1,478 +1,478 @@
#include "GENIEWeightEngine.h"
#ifdef __GENIE_ENABLED__
#ifdef GENIE_PRE_R3
#include "Algorithm/AlgConfigPool.h"
#include "EVGCore/EventRecord.h"
#include "GHEP/GHepRecord.h"
#include "Ntuple/NtpMCEventRecord.h"
-#ifndef __NO_GENIE_REWEIGHT__
+#ifndef __NO_REWEIGHT__
#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"
#include "ReWeight/GSystUncertainty.h"
#ifdef __GENIE_EMP_MECRW_ENABLED
#include "ReWeight/GReWeightXSecEmpiricalMEC.h"
#endif
#endif
#if __GENIE_VERSION__ >= 212
#include "ReWeight/GReWeightNuXSecCCQEaxial.h"
#endif
#else // GENIE v3
#include "Framework/Algorithm/AlgConfigPool.h"
#include "Framework/EventGen/EventRecord.h"
#include "Framework/GHEP/GHepParticle.h"
#include "Framework/GHEP/GHepRecord.h"
#include "Framework/Ntuple/NtpMCEventRecord.h"
#include "Framework/Utils/AppInit.h"
#include "Framework/Utils/RunOpt.h"
using namespace genie;
-#ifndef __NO_GENIE_REWEIGHT__
+#ifndef __NO_REWEIGHT__
#include "RwCalculators/GReWeightAGKY.h"
#include "RwCalculators/GReWeightDISNuclMod.h"
#include "RwCalculators/GReWeightFGM.h"
#include "RwCalculators/GReWeightFZone.h"
#include "RwCalculators/GReWeightINuke.h"
#include "RwCalculators/GReWeightNonResonanceBkg.h"
#include "RwCalculators/GReWeightNuXSecCCQE.h"
#include "RwCalculators/GReWeightNuXSecCCQEaxial.h"
#include "RwCalculators/GReWeightNuXSecCCQEvec.h"
#include "RwCalculators/GReWeightNuXSecCCRES.h"
#include "RwCalculators/GReWeightNuXSecCOH.h"
#include "RwCalculators/GReWeightNuXSecDIS.h"
#include "RwCalculators/GReWeightNuXSecNC.h"
#include "RwCalculators/GReWeightNuXSecNCEL.h"
#include "RwCalculators/GReWeightNuXSecNCRES.h"
#include "RwCalculators/GReWeightResonanceDecay.h"
#include "RwCalculators/GReWeightXSecEmpiricalMEC.h"
#include "RwFramework/GSystUncertainty.h"
using namespace genie::rew;
#endif
#endif
#endif
#include "FitLogger.h"
GENIEWeightEngine::GENIEWeightEngine(std::string name) {
#ifdef __GENIE_ENABLED__
-#ifndef __NO_GENIE_REWEIGHT__
+#ifndef __NO_REWEIGHT__
// Setup the NEUT Reweight engien
fCalcName = name;
NUIS_LOG(DEB, "Setting up GENIE RW : " << fCalcName);
// Create RW Engine suppressing cout
StopTalking();
fGenieRW = new genie::rew::GReWeight();
// Get List of Vetos (Just for debugging)
std::string rw_engine_list =
FitPar::Config().GetParS("FitWeight_fGenieRW_veto");
bool xsec_ncel = rw_engine_list.find("xsec_ncel") == std::string::npos;
bool xsec_ccqe = rw_engine_list.find("xsec_ccqe") == std::string::npos;
bool xsec_coh = rw_engine_list.find("xsec_coh") == std::string::npos;
bool xsec_nnres = rw_engine_list.find("xsec_nonresbkg") == std::string::npos;
bool xsec_nudis = rw_engine_list.find("nuclear_dis") == std::string::npos;
bool xsec_resdec =
rw_engine_list.find("hadro_res_decay") == std::string::npos;
bool xsec_fzone = rw_engine_list.find("hadro_intranuke") == std::string::npos;
bool xsec_intra = rw_engine_list.find("hadro_fzone") == std::string::npos;
bool xsec_agky = rw_engine_list.find("hadro_agky") == std::string::npos;
bool xsec_qevec = rw_engine_list.find("xsec_ccqe_vec") == std::string::npos;
bool xsec_dis = rw_engine_list.find("xsec_dis") == std::string::npos;
bool xsec_nc = rw_engine_list.find("xsec_nc") == std::string::npos;
bool xsec_ccres = rw_engine_list.find("xsec_ccres") == std::string::npos;
bool xsec_ncres = rw_engine_list.find("xsec_ncres") == std::string::npos;
bool xsec_nucqe = rw_engine_list.find("nuclear_qe") == std::string::npos;
bool xsec_qeaxial =
rw_engine_list.find("xsec_ccqe_axial") == std::string::npos;
#ifdef __GENIE_EMP_MECRW_ENABLED
bool xsec_empMEC = rw_engine_list.find("xsec_empMEC") == std::string::npos;
#endif
#ifndef GENIE_PRE_R3
genie::RunOpt *grunopt = genie::RunOpt::Instance();
grunopt->EnableBareXSecPreCalc(true);
grunopt->SetEventGeneratorList(Config::GetParS("GENIEEventGeneratorList"));
if(!Config::HasPar("GENIETune")){
NUIS_ABORT("GENIE tune was not specified, this is required when reweighting GENIE V3+ events. Add a config parameter like: <config GENIETune=\"G18_10a_02_11a\" /> to your nuisance card.");
}
grunopt->SetTuneName(Config::GetParS("GENIETune"));
grunopt->BuildTune();
std::string genv =
std::string(getenv("GENIE")) + "/config/Messenger_laconic.xml";
genie::utils::app_init::MesgThresholds(genv);
#endif
// Now actually add the RW Calcs
if (xsec_ncel)
fGenieRW->AdoptWghtCalc("xsec_ncel", new genie::rew::GReWeightNuXSecNCEL);
if (xsec_ccqe) {
fGenieRW->AdoptWghtCalc("xsec_ccqe", new genie::rew::GReWeightNuXSecCCQE);
// (dynamic_cast<GReWeightNuXSecCCQE*> (fGenieRW->WghtCalc("xsec_ccqe")))
// ->SetXSecModel( FitPar::Config().GetParS("GENIEXSecModelCCQE") );
}
#ifdef __GENIE_EMP_MECRW_ENABLED
if (xsec_empMEC) {
fGenieRW->AdoptWghtCalc("xsec_empMEC",
new genie::rew::GReWeightXSecEmpiricalMEC);
}
#endif
if (xsec_coh) {
fGenieRW->AdoptWghtCalc("xsec_coh", new genie::rew::GReWeightNuXSecCOH());
// (dynamic_cast<GReWeightNuXSecCOH*> (fGenieRW->WghtCalc("xsec_coh")))
// ->SetXSecModel( FitPar::Config().GetParS("GENIEXSecModelCOH") );
}
if (xsec_nnres)
fGenieRW->AdoptWghtCalc("xsec_nonresbkg",
new genie::rew::GReWeightNonResonanceBkg);
if (xsec_nudis)
fGenieRW->AdoptWghtCalc("nuclear_dis", new genie::rew::GReWeightDISNuclMod);
if (xsec_resdec)
fGenieRW->AdoptWghtCalc("hadro_res_decay",
new genie::rew::GReWeightResonanceDecay);
if (xsec_fzone)
fGenieRW->AdoptWghtCalc("hadro_fzone", new genie::rew::GReWeightFZone);
if (xsec_intra)
fGenieRW->AdoptWghtCalc("hadro_intranuke", new genie::rew::GReWeightINuke);
if (xsec_agky)
fGenieRW->AdoptWghtCalc("hadro_agky", new genie::rew::GReWeightAGKY);
if (xsec_qevec)
fGenieRW->AdoptWghtCalc("xsec_ccqe_vec",
new genie::rew::GReWeightNuXSecCCQEvec);
#if __GENIE_VERSION__ >= 212
if (xsec_qeaxial)
fGenieRW->AdoptWghtCalc("xsec_ccqe_axial",
new genie::rew::GReWeightNuXSecCCQEaxial);
#endif
if (xsec_dis)
fGenieRW->AdoptWghtCalc("xsec_dis", new genie::rew::GReWeightNuXSecDIS);
if (xsec_nc)
fGenieRW->AdoptWghtCalc("xsec_nc", new genie::rew::GReWeightNuXSecNC);
if (xsec_ccres) {
#if __GENIE_VERSION__ < 213
fGenieRW->AdoptWghtCalc("xsec_ccres", new genie::rew::GReWeightNuXSecCCRES);
#else
fGenieRW->AdoptWghtCalc(
"xsec_ccres",
new genie::rew::GReWeightNuXSecCCRES(
FitPar::Config().GetParS("GENIEXSecModelCCRES"), "Default"));
#endif
}
if (xsec_ncres)
fGenieRW->AdoptWghtCalc("xsec_ncres", new genie::rew::GReWeightNuXSecNCRES);
if (xsec_nucqe)
fGenieRW->AdoptWghtCalc("nuclear_qe", new genie::rew::GReWeightFGM);
// Set the CCQE reweighting style
GReWeightNuXSecCCQE *rwccqe =
dynamic_cast<GReWeightNuXSecCCQE *>(fGenieRW->WghtCalc("xsec_ccqe"));
// For MaCCQE reweighting
std::string ccqetype = FitPar::Config().GetParS("GENIEWeightEngine_CCQEMode");
if (ccqetype == "kModeMa") {
NUIS_LOG(DEB, "Setting GENIE ReWeight CCQE to kModeMa");
rwccqe->SetMode(GReWeightNuXSecCCQE::kModeMa);
} else if (ccqetype == "kModeNormAndMaShape") {
NUIS_LOG(DEB, "Setting GENIE ReWeight CCQE to kModeNormAndMaShape");
rwccqe->SetMode(GReWeightNuXSecCCQE::kModeNormAndMaShape);
// For z-expansion reweighting
} else if (ccqetype == "kModeZExp") {
NUIS_LOG(DEB, "Setting GENIE ReWeight CCQE to kModeZExp");
rwccqe->SetMode(GReWeightNuXSecCCQE::kModeZExp);
} else {
NUIS_ERR(FTL, "Did not find specified GENIE ReWeight CCQE mode");
NUIS_ABORT("You provided: " << ccqetype << " in parameters/config.xml");
}
#if (__GENIE_VERSION__ >= 212) and \
(__GENIE_VERSION__ <= \
300) // This doesn't currently work as is for GENIE v3, but the reweighting
// in v3 supposedly does similar checks anyway.
// Check the UserPhysicsOptions too!
AlgConfigPool *Pool = genie::AlgConfigPool::Instance();
Registry *full = Pool->GlobalParameterList();
std::string name_ax = full->GetAlg("AxialFormFactorModel").name;
std::string config_ax = full->GetAlg("AxialFormFactorModel").config;
if (name_ax == "genie::DipoleAxialFormFactorModel" &&
ccqetype == "kModeZExp") {
NUIS_ERR(
FTL,
"Trying to run Z Expansion reweighting with Llewelyn-Smith model.");
NUIS_ERR(FTL, "Please check your "
<< std::getenv("GENIE")
<< "/config/UserPhysicsOptions.xml to match generated");
NUIS_ERR(FTL, "You're telling me " << name_ax << "/" << config_ax);
NUIS_ABORT("Also check your "
<< std::getenv("NUISANCE")
<< "/parameters/config.xml GENIEWeightEngine_CCQEMode: "
<< ccqetype);
}
if (name_ax == "genie::ZExpAxialFormFactorModel" && ccqetype != "kModeZExp") {
NUIS_ERR(
FTL,
"Trying to run Llewelyn-Smith reweighting with Z Expansion model.");
NUIS_ERR(FTL, "Please change your "
<< std::getenv("GENIE")
<< "/config/UserPhysicsOptions.xml to match generated");
NUIS_ERR(FTL, "You're telling me " << name_ax << "/" << config_ax);
NUIS_ABORT("Also check your "
<< std::getenv("NUISANCE")
<< "/parameters/config.xml GENIEWeightEngine_CCQEMode: "
<< ccqetype);
}
std::string name_qelcc =
full->GetAlg("XSecModel@genie::EventGenerator/QEL-CC").name;
std::string config_qelcc =
full->GetAlg("XSecModel@genie::EventGenerator/QEL-CC").config;
if (config_qelcc == "Default" && ccqetype == "kModeZExp") {
NUIS_ERR(
FTL,
"Trying to run Z Expansion reweighting with Llewelyn-Smith model.");
NUIS_ERR(FTL, "Please change your "
<< std::getenv("GENIE")
<< "/config/UserPhysicsOptions.xml to match generated");
NUIS_ERR(FTL, "You're telling me " << name_qelcc << "/" << config_qelcc);
NUIS_ABORT("Also check your "
<< std::getenv("NUISANCE")
<< "/parameters/config.xml GENIEWeightEngine_CCQEMode: "
<< ccqetype);
}
if (config_qelcc == "ZExp" && ccqetype != "kModeZExp") {
NUIS_ERR(
FTL,
"Trying to run Llewelyn-Smith reweighting with Z Expansion model.");
NUIS_ERR(FTL, "Please change your "
<< std::getenv("GENIE")
<< "/config/UserPhysicsOptions.xml to match generated");
NUIS_ERR(FTL, "You're telling me " << name_qelcc << "/" << config_qelcc);
NUIS_ABORT("Also check your "
<< std::getenv("NUISANCE")
<< "/parameters/config.xml GENIEWeightEngine_CCQEMode: "
<< ccqetype);
}
#endif
if (xsec_ccres) {
// Default to include shape and normalization changes for CCRES (can be
// changed downstream if desired)
GReWeightNuXSecCCRES *rwccres =
dynamic_cast<GReWeightNuXSecCCRES *>(fGenieRW->WghtCalc("xsec_ccres"));
std::string marestype =
FitPar::Config().GetParS("GENIEWeightEngine_CCRESMode");
if (!marestype.compare("kModeNormAndMaMvShape")) {
rwccres->SetMode(GReWeightNuXSecCCRES::kModeNormAndMaMvShape);
} else if (!marestype.compare("kModeMaMv")) {
rwccres->SetMode(GReWeightNuXSecCCRES::kModeMaMv);
} else {
NUIS_ABORT("Unkown MARES Mode in GENIE Weight Engine : " << marestype);
}
}
if (xsec_ncres) {
// Default to include shape and normalization changes for NCRES (can be
// changed downstream if desired)
GReWeightNuXSecNCRES *rwncres =
dynamic_cast<GReWeightNuXSecNCRES *>(fGenieRW->WghtCalc("xsec_ncres"));
rwncres->SetMode(GReWeightNuXSecNCRES::kModeMaMv);
}
if (xsec_dis) {
// Default to include shape and normalization changes for DIS (can be
// changed downstream if desired)
GReWeightNuXSecDIS *rwdis =
dynamic_cast<GReWeightNuXSecDIS *>(fGenieRW->WghtCalc("xsec_dis"));
rwdis->SetMode(GReWeightNuXSecDIS::kModeABCV12u);
// Set Abs Twk Config
fIsAbsTwk = (FitPar::Config().GetParB("setabstwk"));
}
// allow cout again
StartTalking();
#else
NUIS_ERR(FTL,
"GENIE ReWeight is __NOT ENABLED__ in GENIE and you're trying to "
"run NUISANCE with it enabled");
NUIS_ERR(FTL, "Check your genie-config --libs for reweighting");
NUIS_ERR(FTL, "If not present you need to recompile GENIE");
NUIS_ABORT("If present you need to contact NUISANCE authors");
#endif
#endif
};
void GENIEWeightEngine::IncludeDial(std::string name, double startval) {
#ifdef __GENIE_ENABLED__
-#ifndef __NO_GENIE_REWEIGHT__
+#ifndef __NO_REWEIGHT__
// Get First enum
int nuisenum = Reweight::ConvDial(name, kGENIE);
// Check ZExp sillyness in GENIE
// If ZExpansion parameters are used we need to set a different mode in GENIE
// ReWeight... GENIE doesn't have a setter either...
#if (__GENIE_VERSION__ >= 212) and (__GENIE_VERSION__ <= 300)
std::string ccqetype = FitPar::Config().GetParS("GENIEWeightEngine_CCQEMode");
if (ccqetype != "kModeZExp" &&
(name == "ZExpA1CCQE" || name == "ZExpA2CCQE" || name == "ZExpA3CCQE" ||
name == "ZExpA4CCQE")) {
NUIS_ERR(FTL, "Found a Z-expansion parameter in GENIE although the GENIE "
"ReWeighting engine is set to use Llewelyn-Smith and MaQE!");
NUIS_ABORT("Change your GENIE UserPhysicsOptions.xml in "
<< std::getenv("GENIE")
<< "/config/UserPhysicsOptions.xml to match requirements");
}
if ((ccqetype != "kModeMa" && ccqetype != "kModeMaNormAndMaShape") &&
(name == "MaCCQE")) {
NUIS_ERR(FTL, "Found MaCCQE parameter in GENIE although the GENIE "
"ReWeighting engine is set to not use this!");
NUIS_ABORT("Change your GENIE UserPhysicsOptions.xml in "
<< std::getenv("GENIE")
<< "/config/UserPhysicsOptions.xml to match requirements");
}
#endif
// Setup Maps
fEnumIndex[nuisenum]; // = std::vector<size_t>(0);
fNameIndex[name]; // = std::vector<size_t>(0);
// Split by commas
std::vector<std::string> allnames = GeneralUtils::ParseToStr(name, ",");
for (uint i = 0; i < allnames.size(); i++) {
std::string singlename = allnames[i];
// Get RW
genie::rew::GSyst_t rwsyst = GSyst::FromString(singlename);
// Fill Maps
int index = fValues.size();
fValues.push_back(0.0);
fGENIESysts.push_back(rwsyst);
// Initialize dial
NUIS_LOG(DEB, "Registering " << singlename << " from " << name);
fGenieRW->Systematics().Init(fGENIESysts[index]);
// If Absolute
if (fIsAbsTwk) {
GSystUncertainty::Instance()->SetUncertainty(rwsyst, 1.0, 1.0);
}
// Setup index
fEnumIndex[nuisenum].push_back(index);
fNameIndex[name].push_back(index);
}
// Set Value if given
if (startval != -999.9) {
SetDialValue(nuisenum, startval);
}
#endif
#endif
}
void GENIEWeightEngine::SetDialValue(int nuisenum, double val) {
#ifdef __GENIE_ENABLED__
-#ifndef __NO_GENIE_REWEIGHT__
+#ifndef __NO_REWEIGHT__
std::vector<size_t> indices = fEnumIndex[nuisenum];
for (uint i = 0; i < indices.size(); i++) {
fValues[indices[i]] = val;
fGenieRW->Systematics().Set(fGENIESysts[indices[i]], val);
}
#endif
#endif
}
void GENIEWeightEngine::SetDialValue(std::string name, double val) {
#ifdef __GENIE_ENABLED__
-#ifndef __NO_GENIE_REWEIGHT__
+#ifndef __NO_REWEIGHT__
std::vector<size_t> indices = fNameIndex[name];
for (uint i = 0; i < indices.size(); i++) {
fValues[indices[i]] = val;
fGenieRW->Systematics().Set(fGENIESysts[indices[i]], val);
}
#endif
#endif
}
void GENIEWeightEngine::Reconfigure(bool silent) {
#ifdef __GENIE_ENABLED__
-#ifndef __NO_GENIE_REWEIGHT__
+#ifndef __NO_REWEIGHT__
// Hush now...
if (silent)
StopTalking();
// Reconf
fGenieRW->Reconfigure();
fGenieRW->Print();
// Shout again
if (silent)
StartTalking();
#endif
#endif
}
double GENIEWeightEngine::CalcWeight(BaseFitEvt *evt) {
double rw_weight = 1.0;
#ifdef __GENIE_ENABLED__
-#ifndef __NO_GENIE_REWEIGHT__
+#ifndef __NO_REWEIGHT__
// Make nom weight
if (!evt) {
NUIS_ABORT("evt not found : " << evt);
}
// Skip Non GENIE
if (evt->fType != kGENIE)
return 1.0;
if (!(evt->genie_event)) {
NUIS_ABORT("evt->genie_event not found!" << evt->genie_event);
}
if (!(evt->genie_event->event)) {
NUIS_ABORT("evt->genie_event->event GHepRecord not found!"
<< (evt->genie_event->event));
}
if (!fGenieRW) {
NUIS_ABORT("GENIE RW Not Found!" << fGenieRW);
}
rw_weight = fGenieRW->CalcWeight(*(evt->genie_event->event));
// std::cout << "Returning GENIE Weight for electron scattering = " <<
// rw_weight << std::endl;
// if (rw_weight != 1.0 )std::cout << "mode=" << evt->Mode << " rw_weight = "
// << rw_weight << std::endl;
#endif
#endif
// Return rw_weight
return rw_weight;
}
diff --git a/src/Reweight/GENIEWeightEngine.h b/src/Reweight/GENIEWeightEngine.h
index 4c39f68..b2e6dc5 100644
--- a/src/Reweight/GENIEWeightEngine.h
+++ b/src/Reweight/GENIEWeightEngine.h
@@ -1,46 +1,46 @@
#ifndef WEIGHT_ENGINE_GENIE_H
#define WEIGHT_ENGINE_GENIE_H
#ifdef __GENIE_ENABLED__
#ifdef GENIE_PRE_R3
-#ifndef __NO_GENIE_REWEIGHT__
+#ifndef __NO_REWEIGHT__
#include "ReWeight/GSyst.h"
#include "ReWeight/GReWeight.h"
#endif
#else
-#ifndef __NO_GENIE_REWEIGHT__
+#ifndef __NO_REWEIGHT__
#include "RwFramework/GSyst.h"
#include "RwFramework/GReWeight.h"
using namespace genie;
using namespace genie::rew;
#endif
#endif
#endif
#include "GeneratorUtils.h"
#include "WeightEngineBase.h"
#include "FitWeight.h"
class GENIEWeightEngine : public WeightEngineBase {
public:
GENIEWeightEngine(std::string name);
~GENIEWeightEngine() {};
void IncludeDial(std::string name, double startval);
void SetDialValue(int rwenum, double val);
void SetDialValue(std::string name, double val);
void Reconfigure(bool silent = false);
double CalcWeight(BaseFitEvt* evt);
inline bool NeedsEventReWeight() { return true; };
#ifdef __GENIE_ENABLED__
-#ifndef __NO_GENIE_REWEIGHT__
+#ifndef __NO_REWEIGHT__
std::vector<genie::rew::GSyst_t> fGENIESysts;
genie::rew::GReWeight* fGenieRW; //!< Genie RW Object
#endif
#endif
};
#endif
diff --git a/src/Reweight/NEUTWeightEngine.cxx b/src/Reweight/NEUTWeightEngine.cxx
index 88b515f..966a09b 100644
--- a/src/Reweight/NEUTWeightEngine.cxx
+++ b/src/Reweight/NEUTWeightEngine.cxx
@@ -1,179 +1,178 @@
#include "NEUTWeightEngine.h"
NEUTWeightEngine::NEUTWeightEngine(std::string name) {
-#ifdef __NEUT_ENABLED__
+#if defined(__NEUT_ENABLED__) and !defined(__NO_REWEIGHT__)
// Setup the NEUT Reweight engien
fCalcName = name;
NUIS_LOG(FIT, "Setting up NEUT RW : " << fCalcName);
// Create RW Engine suppressing cout
StopTalking();
fNeutRW = new neut::rew::NReWeight();
TDirectory *olddir = gDirectory;
// get list of vetoed calc engines (just for debug really)
std::string rw_engine_list =
FitPar::Config().GetParS("FitWeight_fNeutRW_veto");
bool xsec_ccqe = rw_engine_list.find("xsec_ccqe") == std::string::npos;
bool xsec_res = rw_engine_list.find("xsec_res") == std::string::npos;
bool xsec_ccres = rw_engine_list.find("xsec_ccres") == std::string::npos;
bool xsec_coh = rw_engine_list.find("xsec_coh") == std::string::npos;
bool xsec_dis = rw_engine_list.find("xsec_dis") == std::string::npos;
bool xsec_ncel = rw_engine_list.find("xsec_ncel") == std::string::npos;
bool xsec_nc = rw_engine_list.find("xsec_nc") == std::string::npos;
bool xsec_ncres = rw_engine_list.find("xsec_ncres") == std::string::npos;
bool nucl_casc = rw_engine_list.find("nucl_casc") == std::string::npos;
bool nucl_piless = rw_engine_list.find("nucl_piless") == std::string::npos;
// Activate each calc engine
if (xsec_ccqe)
fNeutRW->AdoptWghtCalc("xsec_ccqe", new neut::rew::NReWeightNuXSecCCQE);
if (xsec_res)
fNeutRW->AdoptWghtCalc("xsec_res", new neut::rew::NReWeightNuXSecRES);
if (xsec_ccres)
fNeutRW->AdoptWghtCalc("xsec_ccres", new neut::rew::NReWeightNuXSecCCRES);
if (xsec_coh)
fNeutRW->AdoptWghtCalc("xsec_coh", new neut::rew::NReWeightNuXSecCOH);
if (xsec_dis)
fNeutRW->AdoptWghtCalc("xsec_dis", new neut::rew::NReWeightNuXSecDIS);
if (xsec_ncel)
fNeutRW->AdoptWghtCalc("xsec_ncel", new neut::rew::NReWeightNuXSecNCEL);
if (xsec_nc)
fNeutRW->AdoptWghtCalc("xsec_nc", new neut::rew::NReWeightNuXSecNC);
if (xsec_ncres)
fNeutRW->AdoptWghtCalc("xsec_ncres", new neut::rew::NReWeightNuXSecNCRES);
if (nucl_casc)
fNeutRW->AdoptWghtCalc("nucl_casc", new neut::rew::NReWeightCasc);
if (nucl_piless)
fNeutRW->AdoptWghtCalc("nucl_piless", new neut::rew::NReWeightNuclPiless);
fNeutRW->Reconfigure();
olddir->cd();
// Set Abs Twk Config
fIsAbsTwk = (FitPar::Config().GetParB("setabstwk"));
// allow cout again
StartTalking();
#else
NUIS_ABORT("NEUT RW NOT ENABLED!");
#endif
};
void NEUTWeightEngine::IncludeDial(std::string name, double startval) {
-#ifdef __NEUT_ENABLED__
+#if defined(__NEUT_ENABLED__) and !defined(__NO_REWEIGHT__)
// Get First enum
int nuisenum = Reweight::ConvDial(name, kNEUT);
// Setup Maps
fEnumIndex[nuisenum]; // = std::vector<size_t>(0);
fNameIndex[name]; // = std::vector<size_t>(0);
// Split by commas
std::vector<std::string> allnames = GeneralUtils::ParseToStr(name, ",");
for (uint i = 0; i < allnames.size(); i++) {
std::string singlename = allnames[i];
// Get Syst
neut::rew::NSyst_t gensyst = NSyst::FromString(singlename);
// Fill Maps
int index = fValues.size();
fValues.push_back(0.0);
fNEUTSysts.push_back(gensyst);
// Initialize dial
NUIS_LOG(FIT, "Registering " << singlename << " dial.");
fNeutRW->Systematics().Init(fNEUTSysts[index]);
// If Absolute
if (fIsAbsTwk) {
NSystUncertainty::Instance()->SetUncertainty(fNEUTSysts[index], 1.0, 1.0);
}
// Setup index
fEnumIndex[nuisenum].push_back(index);
fNameIndex[name].push_back(index);
}
// Set Value if given
if (startval != -999.9) {
SetDialValue(nuisenum, startval);
}
-
#endif
}
void NEUTWeightEngine::SetDialValue(int nuisenum, double val) {
-#ifdef __NEUT_ENABLED__
+#if defined(__NEUT_ENABLED__) and !defined(__NO_REWEIGHT__)
std::vector<size_t> indices = fEnumIndex[nuisenum];
for (uint i = 0; i < indices.size(); i++) {
fValues[indices[i]] = val;
- NUIS_LOG(FIT, "Setting Dial Value for " << nuisenum << " " << i << " "
- << indices[i] << " "
- << fValues[indices[i]]
- << " Enum: " << fNEUTSysts[indices[i]]);
+ NUIS_LOG(FIT, "Setting Dial Value for "
+ << nuisenum << " " << i << " " << indices[i] << " "
+ << fValues[indices[i]]
+ << " Enum: " << fNEUTSysts[indices[i]]);
fNeutRW->Systematics().Set(fNEUTSysts[indices[i]], val);
}
#endif
}
void NEUTWeightEngine::SetDialValue(std::string name, double val) {
-#ifdef __NEUT_ENABLED__
+#if defined(__NEUT_ENABLED__) and !defined(__NO_REWEIGHT__)
std::vector<size_t> indices = fNameIndex[name];
for (uint i = 0; i < indices.size(); i++) {
fValues[indices[i]] = val;
NUIS_LOG(FIT, "Setting Dial Value for "
- << name << " = " << i << " " << indices[i] << " "
- << fValues[indices[i]]
- << " Enum: " << fNEUTSysts[indices[i]]);
+ << name << " = " << i << " " << indices[i] << " "
+ << fValues[indices[i]]
+ << " Enum: " << fNEUTSysts[indices[i]]);
fNeutRW->Systematics().Set(fNEUTSysts[indices[i]], val);
}
#endif
}
void NEUTWeightEngine::Reconfigure(bool silent) {
-#ifdef __NEUT_ENABLED__
+#if defined(__NEUT_ENABLED__) and !defined(__NO_REWEIGHT__)
// Hush now...
if (silent)
StopTalking();
// Reconf
fNeutRW->Reconfigure();
// if (LOG_LEVEL(DEB)){
fNeutRW->Print();
// }
// Shout again
if (silent)
StartTalking();
#endif
}
double NEUTWeightEngine::CalcWeight(BaseFitEvt *evt) {
double rw_weight = 1.0;
-#ifdef __NEUT_ENABLED__
+#if defined(__NEUT_ENABLED__) and !defined(__NO_REWEIGHT__)
// Skip Non NEUT
if (evt->fType != kNEUT)
return 1.0;
// Hush now
StopTalking();
// Fill NEUT Common blocks
NEUTUtils::FillNeutCommons(evt->fNeutVect);
// Call Weight calculation
rw_weight = fNeutRW->CalcWeight();
// Speak Now
StartTalking();
#endif
// Return rw_weight
return rw_weight;
}
diff --git a/src/Reweight/NEUTWeightEngine.h b/src/Reweight/NEUTWeightEngine.h
index a4ad648..7cc7bb0 100644
--- a/src/Reweight/NEUTWeightEngine.h
+++ b/src/Reweight/NEUTWeightEngine.h
@@ -1,53 +1,55 @@
#ifndef WEIGHT_ENGINE_NEUT_H
#define WEIGHT_ENGINE_NEUT_H
#include "FitLogger.h"
#ifdef __NEUT_ENABLED__
+#ifndef __NO_REWEIGHT__
+#include "NEUTInputHandler.h"
#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"
-#include "NEUTInputHandler.h"
+#endif
#endif
-
+#include "FitWeight.h"
#include "GeneratorUtils.h"
#include "WeightEngineBase.h"
-#include "FitWeight.h"
class NEUTWeightEngine : public WeightEngineBase {
public:
- NEUTWeightEngine(std::string name);
- ~NEUTWeightEngine() {};
+ NEUTWeightEngine(std::string name);
+ ~NEUTWeightEngine(){};
- void IncludeDial(std::string name, double startval);
+ void IncludeDial(std::string name, double startval);
- void SetDialValue(std::string name, double val);
- void SetDialValue(int nuisenum, double val);
+ void SetDialValue(std::string name, double val);
+ void SetDialValue(int nuisenum, double val);
- void Reconfigure(bool silent = false);
+ void Reconfigure(bool silent = false);
- double CalcWeight(BaseFitEvt* evt);
-
- inline bool NeedsEventReWeight() { return true; };
+ double CalcWeight(BaseFitEvt *evt);
+ inline bool NeedsEventReWeight() { return true; };
#ifdef __NEUT_ENABLED__
- std::vector<neut::rew::NSyst_t> fNEUTSysts;
- neut::rew::NReWeight* fNeutRW;
+#ifndef __NO_REWEIGHT__
+ std::vector<neut::rew::NSyst_t> fNEUTSysts;
+ neut::rew::NReWeight *fNeutRW;
+#endif
#endif
};
#endif
diff --git a/src/Reweight/WeightUtils.cxx b/src/Reweight/WeightUtils.cxx
index 8429a6a..b381b80 100644
--- a/src/Reweight/WeightUtils.cxx
+++ b/src/Reweight/WeightUtils.cxx
@@ -1,592 +1,595 @@
#include "WeightUtils.h"
#include "FitLogger.h"
+#ifndef __NO_REWEIGHT__
#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 "NIWGSyst.h"
#endif
#ifdef __NEUT_ENABLED__
#include "NReWeight.h"
#include "NSyst.h"
#endif
#ifdef __NUWRO_REWEIGHT_ENABLED__
#include "NuwroReWeight.h"
#include "NuwroSyst.h"
#endif
#ifdef __GENIE_ENABLED__
#ifdef GENIE_PRE_R3
#ifndef __NO_GENIE_REWEIGHT__
#include "ReWeight/GReWeight.h"
#include "ReWeight/GSyst.h"
#endif
#else
using namespace genie;
#ifndef __NO_GENIE_REWEIGHT__
#include "RwFramework/GReWeight.h"
#include "RwFramework/GSyst.h"
using namespace genie::rew;
#endif
#endif
#endif
#ifdef __NOVA_ENABLED__
#include "NOvARwgtEngine.h"
#endif
+#endif // end of no reweight
+
#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;
std::ifstream card(
(GeneralUtils::GetTopLevelDir() + "/parameters/dial_conversion.card")
.c_str(),
std::ifstream::in);
while (std::getline(card >> std::ws, line, '\n')) {
std::vector<std::string> 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(),
std::ifstream::in);
while (std::getline(card >> std::ws, line, '\n')) {
std::vector<std::string> 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;
NUIS_LOG(FIT, "AbsToSigma(" << name << ") = " << val << " -> " << conv_val);
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 if (!type.compare("nova_parameter"))
return kNOvARWGT;
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";
}
case kNOvARWGT: {
return "nova_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
NUIS_LOG(DEB, "Getting dial enum " << type << " " << name);
// Select Types
switch (type) {
// NEUT DIAL TYPE
case kNEUT: {
-#ifdef __NEUT_ENABLED__
+#if defined(__NEUT_ENABLED__) && !defined(__NO_REWEIGHT__)
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__
+#if defined(__NIWG_ENABLED__) && !defined(__NO_REWEIGHT__)
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__
+#if defined(__NUWRO_REWEIGHT_ENABLED__) && !defined(__NO_REWEIGHT__)
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: {
-#if defined(__GENIE_ENABLED__) && !defined(__NO_GENIE_REWEIGHT__)
+#if defined(__GENIE_ENABLED__) && !defined(__NO_REWEIGHT__)
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__
+#if defined(__T2KREW_ENABLED__) && !defined(__NO_REWEIGHT__)
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());
NUIS_LOG(FTL, "Getting mode num " << mode_num);
if (!mode_num) {
NUIS_ABORT("Attempting to parse dial name: \""
<< name << "\" as a mode norm dial but failed.");
}
this_enum = 60 + mode_num + offset;
break;
}
}
// If Not Enabled
if (this_enum == Reweight::kNoTypeFound) {
NUIS_ERR(FTL, "RW Engine not supported for " << FitBase::ConvDialType(type));
NUIS_ABORT("Check dial " << name);
}
// If Not Found
if (this_enum == Reweight::kNoDialFound) {
NUIS_ABORT("Dial " << name << " not found.");
}
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 > kNOvARWGT ? Reweight::kNoDialFound : t;
}
int Reweight::RemoveDialType(int type) { return (type % 1000); }
int Reweight::NEUTEnumFromName(std::string const &name) {
-#ifdef __NEUT_ENABLED__
+#if defined(__NEUT_ENABLED__) && !defined(__NO_REWEIGHT__)
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__
+#if defined(__NIWG_ENABLED__) && !defined(__NO_REWEIGHT__)
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__
+#if defined(__NUWRO_REWEIGHT_ENABLED__) && !defined(__NO_REWEIGHT__)
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) {
-#if defined(__GENIE_ENABLED__) && !defined(__NO_GENIE_REWEIGHT__)
+#if defined(__GENIE_ENABLED__) && !defined(__NO_REWEIGHT__)
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__
+#if defined(__T2KREW_ENABLED__) && !defined(__NO_REWEIGHT__)
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;
#ifdef __NOVA_ENABLED__
case kNOvARWGT:
genenum = NOvARwgtEngine::GetWeightGeneratorIndex(name);
break;
#endif
default:
genenum = Reweight::kNoTypeFound;
break;
}
// Throw if required.
if (exceptions) {
// If Not Enabled
if (genenum == Reweight::kGeneratorNotBuilt) {
NUIS_ERR(FTL,
"RW Engine not supported for " << FitBase::ConvDialType(type));
NUIS_ABORT("Check dial " << name);
}
// If no type enabled
if (genenum == Reweight::kNoTypeFound) {
NUIS_ABORT("Type mismatch inside ConvDialEnum");
}
// If Not Found
if (genenum == Reweight::kNoDialFound) {
NUIS_ABORT("Dial " << name << " not found.");
}
}
// 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];
}
}
NUIS_LOG(FIT, "Cannot find dial with enum = " << nuisenum);
return "";
}
diff --git a/src/Statistical/StatUtils.cxx b/src/Statistical/StatUtils.cxx
index c4fa68b..586cbed 100644
--- a/src/Statistical/StatUtils.cxx
+++ b/src/Statistical/StatUtils.cxx
@@ -1,1406 +1,1448 @@
// 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "StatUtils.h"
#include "GeneralUtils.h"
#include "NuisConfig.h"
#include "TH1D.h"
//*******************************************************************
Double_t StatUtils::GetChi2FromDiag(TH1D *data, TH1D *mc, TH1I *mask) {
//*******************************************************************
Double_t Chi2 = 0.0;
- TH1D *calc_data = (TH1D *)data->Clone();
- TH1D *calc_mc = (TH1D *)mc->Clone();
+ TH1D *calc_data = (TH1D *)data->Clone("calc_data");
+ calc_data->SetDirectory(NULL);
+ TH1D *calc_mc = (TH1D *)mc->Clone("calc_mc");
+ calc_mc->SetDirectory(NULL);
// Add MC Error to data if required
if (FitPar::Config().GetParB("addmcerror")) {
for (int i = 0; i < calc_data->GetNbinsX(); i++) {
double dterr = calc_data->GetBinError(i + 1);
double mcerr = calc_mc->GetBinError(i + 1);
if (dterr > 0.0) {
calc_data->SetBinError(i + 1, sqrt(dterr * dterr + mcerr * mcerr));
}
}
}
// Apply masking if required
if (mask) {
calc_data = ApplyHistogramMasking(data, mask);
+ calc_data->SetDirectory(NULL);
calc_mc = ApplyHistogramMasking(mc, mask);
+ calc_mc->SetDirectory(NULL);
}
// Iterate over bins in X
for (int i = 0; i < calc_data->GetNbinsX(); i++) {
// Ignore bins with zero data or zero bin error
if (calc_data->GetBinError(i + 1) <= 0.0 ||
calc_data->GetBinContent(i + 1) == 0.0)
continue;
// Take mc data difference
double diff =
calc_data->GetBinContent(i + 1) - calc_mc->GetBinContent(i + 1);
double err = calc_data->GetBinError(i + 1);
Chi2 += (diff * diff) / (err * err);
}
// cleanup
delete calc_data;
delete calc_mc;
return Chi2;
};
//*******************************************************************
Double_t StatUtils::GetChi2FromDiag(TH2D *data, TH2D *mc, TH2I *map,
TH2I *mask) {
//*******************************************************************
// Generate a simple map
if (!map)
map = GenerateMap(data);
// Convert to 1D Histograms
TH1D *data_1D = MapToTH1D(data, map);
TH1D *mc_1D = MapToTH1D(mc, map);
TH1I *mask_1D = MapToMask(mask, map);
// Calculate 1D chi2 from 1D Plots
Double_t Chi2 = StatUtils::GetChi2FromDiag(data_1D, mc_1D, mask_1D);
// CleanUp
delete data_1D;
delete mc_1D;
delete mask_1D;
return Chi2;
};
//*******************************************************************
Double_t StatUtils::GetChi2FromCov(TH1D *data, TH1D *mc, TMatrixDSym *invcov,
TH1I *mask, double data_scale,
- double covar_scale) {
+ double covar_scale, TH1D *outchi2perbin) {
//*******************************************************************
Double_t Chi2 = 0.0;
TMatrixDSym *calc_cov = (TMatrixDSym *)invcov->Clone();
TH1D *calc_data = (TH1D *)data->Clone();
TH1D *calc_mc = (TH1D *)mc->Clone();
// If a mask if applied we need to apply it before the matrix is inverted
if (mask) {
calc_cov = ApplyInvertedMatrixMasking(invcov, mask);
calc_data = ApplyHistogramMasking(data, mask);
calc_mc = ApplyHistogramMasking(mc, mask);
}
// Add MC Error to data if required
if (FitPar::Config().GetParB("statutils.addmcerror")) {
// Make temp cov
TMatrixDSym *newcov = StatUtils::GetInvert(calc_cov);
// Add MC err to diag
for (int i = 0; i < calc_data->GetNbinsX(); i++) {
double mcerr = calc_mc->GetBinError(i + 1) * sqrt(covar_scale);
double oldval = (*newcov)(i, i);
NUIS_LOG(FIT,
"Adding cov stat " << mcerr * mcerr << " to " << (*newcov)(i, i));
(*newcov)(i, i) = oldval + mcerr * mcerr;
}
// Reset the calc_cov to new invert
delete calc_cov;
calc_cov = GetInvert(newcov);
// Delete the tempcov
delete newcov;
}
calc_data->Scale(data_scale);
calc_mc->Scale(data_scale);
(*calc_cov) *= covar_scale;
// iterate over bins in X (i,j)
NUIS_LOG(DEB, "START Chi2 Calculation=================");
for (int i = 0; i < calc_data->GetNbinsX(); i++) {
+ double ibin_contrib = 0;
NUIS_LOG(DEB,
"[CHI2] i = " << i << " ["
<< calc_data->GetXaxis()->GetBinLowEdge(i + 1) << " -- "
<< calc_data->GetXaxis()->GetBinUpEdge(i + 1) << "].");
for (int j = 0; j < calc_data->GetNbinsX(); j++) {
NUIS_LOG(DEB, "[CHI2]\t j = "
<< i << " [" << calc_data->GetXaxis()->GetBinLowEdge(j + 1)
<< " -- " << calc_data->GetXaxis()->GetBinUpEdge(j + 1)
<< "].");
if ((calc_data->GetBinContent(i + 1) != 0 ||
calc_mc->GetBinContent(i + 1) != 0) &&
((*calc_cov)(i, j) != 0)) {
NUIS_LOG(DEB,
"[CHI2]\t\t Chi2 contribution (i,j) = (" << i << "," << j << ")");
NUIS_LOG(DEB, "[CHI2]\t\t Data - MC(i) = "
<< calc_data->GetBinContent(i + 1) << " - "
<< calc_mc->GetBinContent(i + 1) << " = "
<< (calc_data->GetBinContent(i + 1) -
calc_mc->GetBinContent(i + 1)));
NUIS_LOG(DEB, "[CHI2]\t\t Data - MC(j) = "
<< calc_data->GetBinContent(j + 1) << " - "
<< calc_mc->GetBinContent(j + 1) << " = "
<< (calc_data->GetBinContent(j + 1) -
calc_mc->GetBinContent(j + 1)));
NUIS_LOG(DEB, "[CHI2]\t\t Covar = " << (*calc_cov)(i, j));
NUIS_LOG(DEB,
"[CHI2]\t\t Cont chi2 = " << ((calc_data->GetBinContent(i + 1) -
calc_mc->GetBinContent(i + 1)) *
(*calc_cov)(i, j) *
(calc_data->GetBinContent(j + 1) -
calc_mc->GetBinContent(j + 1)))
<< " " << Chi2);
-
- Chi2 +=
+ double bin_cont =
((calc_data->GetBinContent(i + 1) - calc_mc->GetBinContent(i + 1)) *
(*calc_cov)(i, j) *
(calc_data->GetBinContent(j + 1) - calc_mc->GetBinContent(j + 1)));
+ Chi2 += bin_cont;
+ ibin_contrib += bin_cont;
} else {
NUIS_LOG(DEB, "Skipping chi2 contribution (i,j) = ("
<< i << "," << j
<< "), Data = " << calc_data->GetBinContent(i + 1)
<< ", MC = " << calc_mc->GetBinContent(i + 1)
<< ", Cov = " << (*calc_cov)(i, j));
Chi2 += 0.;
}
}
+ if (outchi2perbin) {
+ outchi2perbin->SetBinContent(i + 1, ibin_contrib);
+ }
}
// Cleanup
delete calc_cov;
delete calc_data;
delete calc_mc;
return Chi2;
}
//*******************************************************************
Double_t StatUtils::GetChi2FromCov(TH2D *data, TH2D *mc, TMatrixDSym *invcov,
- TH2I *map, TH2I *mask) {
+ TH2I *map, TH2I *mask, TH2D *outchi2perbin) {
//*******************************************************************
// Generate a simple map
if (!map) {
map = StatUtils::GenerateMap(data);
}
// Convert to 1D Histograms
TH1D *data_1D = MapToTH1D(data, map);
TH1D *mc_1D = MapToTH1D(mc, map);
TH1I *mask_1D = MapToMask(mask, map);
+ TH1D *outchi2perbin_1D = NULL;
+ TH1D *outchi2perbin_map_1D = NULL;
+
+ if (outchi2perbin) {
+ outchi2perbin_1D = MapToTH1D(outchi2perbin, map);
+ for (Int_t xbi_it = 0; xbi_it < outchi2perbin->GetXaxis()->GetNbins();
+ ++xbi_it) {
+ for (Int_t ybi_it = 0; ybi_it < outchi2perbin->GetYaxis()->GetNbins();
+ ++ybi_it) {
+ int gbin = outchi2perbin->GetBin(xbi_it + 1, ybi_it + 1);
+ // std::cout << " gbin " << gbin << " corresponds to "
+ // << " x: " << (xbi_it + 1) << ", y: " << (ybi_it + 1)
+ // << std::endl;
+ outchi2perbin->SetBinContent(xbi_it + 1, ybi_it + 1, gbin);
+ }
+ }
+ outchi2perbin_map_1D = MapToTH1D(outchi2perbin, map);
+ }
// Calculate 1D chi2 from 1D Plots
- Double_t Chi2 = StatUtils::GetChi2FromCov(data_1D, mc_1D, invcov, mask_1D);
+ Double_t Chi2 = StatUtils::GetChi2FromCov(data_1D, mc_1D, invcov, mask_1D, 1,
+ 1E76, outchi2perbin_1D);
+ if (outchi2perbin) {
+ for (int xbi_it = 0; xbi_it < outchi2perbin_1D->GetXaxis()->GetNbins();
+ ++xbi_it) {
+ // std::cout << " adding chi2 "
+ // << outchi2perbin_1D->GetBinContent(xbi_it + 1)
+ // << " from 1d bin " << (xbi_it + 1) << " to gbin "
+ // << outchi2perbin_map_1D->GetBinContent(xbi_it + 1) << std::endl;
+ outchi2perbin->SetBinContent(
+ outchi2perbin_map_1D->GetBinContent(xbi_it + 1),
+ outchi2perbin_1D->GetBinContent(xbi_it + 1));
+ }
+ }
// CleanUp
delete data_1D;
delete mc_1D;
delete mask_1D;
+ delete outchi2perbin_1D;
+ delete outchi2perbin_map_1D;
return Chi2;
}
//*******************************************************************
Double_t StatUtils::GetChi2FromSVD(TH1D *data, TH1D *mc, TMatrixDSym *cov,
TH1I *mask) {
//*******************************************************************
Double_t Chi2 = 0.0;
TMatrixDSym *calc_cov = (TMatrixDSym *)cov->Clone();
TH1D *calc_data = (TH1D *)data->Clone();
TH1D *calc_mc = (TH1D *)mc->Clone();
// If a mask if applied we need to apply it before the matrix is inverted
if (mask) {
calc_cov = StatUtils::ApplyMatrixMasking(cov, mask);
calc_data = StatUtils::ApplyHistogramMasking(data, mask);
calc_mc = StatUtils::ApplyHistogramMasking(mc, mask);
}
// Decompose matrix
TDecompSVD LU = TDecompSVD((*calc_cov));
LU.Decompose();
TMatrixDSym *cov_U =
new TMatrixDSym(calc_data->GetNbinsX(), LU.GetU().GetMatrixArray(), "");
TVectorD *cov_S = new TVectorD(LU.GetSig());
// Apply basis rotation before adding up chi2
Double_t rotated_difference = 0.0;
for (int i = 0; i < calc_data->GetNbinsX(); i++) {
rotated_difference = 0.0;
// Rotate basis of Data - MC
for (int j = 0; j < calc_data->GetNbinsY(); j++)
rotated_difference +=
(calc_data->GetBinContent(j + 1) - calc_mc->GetBinContent(j + 1)) *
(*cov_U)(j, i);
// Divide by rotated error cov_S
Chi2 += rotated_difference * rotated_difference * 1E76 / (*cov_S)(i);
}
// Cleanup
delete calc_cov;
delete calc_data;
delete calc_mc;
delete cov_U;
delete cov_S;
return Chi2;
}
//*******************************************************************
Double_t StatUtils::GetChi2FromSVD(TH2D *data, TH2D *mc, TMatrixDSym *cov,
TH2I *map, TH2I *mask) {
//*******************************************************************
// Generate a simple map
if (!map)
map = StatUtils::GenerateMap(data);
// Convert to 1D Histograms
TH1D *data_1D = MapToTH1D(data, map);
TH1D *mc_1D = MapToTH1D(mc, map);
TH1I *mask_1D = MapToMask(mask, map);
// Calculate from 1D
Double_t Chi2 = StatUtils::GetChi2FromSVD(data_1D, mc_1D, cov, mask_1D);
// CleanUp
delete data_1D;
delete mc_1D;
delete mask_1D;
return Chi2;
}
//*******************************************************************
double StatUtils::GetChi2FromEventRate(TH1D *data, TH1D *mc, TH1I *mask) {
//*******************************************************************
// If just an event rate, for chi2 just use Poission Likelihood to calculate
// the chi2 component
double chi2 = 0.0;
TH1D *calc_data = (TH1D *)data->Clone();
TH1D *calc_mc = (TH1D *)mc->Clone();
// Apply masking if required
if (mask) {
calc_data = ApplyHistogramMasking(data, mask);
calc_mc = ApplyHistogramMasking(mc, mask);
}
// Iterate over bins in X
for (int i = 0; i < calc_data->GetNbinsX(); i++) {
double dt = calc_data->GetBinContent(i + 1);
double mc = calc_mc->GetBinContent(i + 1);
if (mc <= 0)
continue;
if (dt <= 0) {
// Only add difference
chi2 += 2 * (mc - dt);
} else {
// Do the chi2 for Poisson distributions
chi2 += 2 * (mc - dt + (dt * log(dt / mc)));
}
/*
LOG(REC)<<"Evt Chi2 cont = "<<i<<" "
<<mc<<" "<<dt<<" "
<<2 * (mc - dt + (dt+0.) * log((dt+0.) / (mc+0.)))
<<" "<<Chi2<<std::endl;
*/
}
// cleanup
delete calc_data;
delete calc_mc;
return chi2;
}
//*******************************************************************
Double_t StatUtils::GetChi2FromEventRate(TH2D *data, TH2D *mc, TH2I *map,
TH2I *mask) {
//*******************************************************************
// Generate a simple map
if (!map)
map = StatUtils::GenerateMap(data);
// Convert to 1D Histograms
TH1D *data_1D = MapToTH1D(data, map);
TH1D *mc_1D = MapToTH1D(mc, map);
TH1I *mask_1D = MapToMask(mask, map);
// Calculate from 1D
Double_t Chi2 = StatUtils::GetChi2FromEventRate(data_1D, mc_1D, mask_1D);
// CleanUp
delete data_1D;
delete mc_1D;
delete mask_1D;
return Chi2;
}
//*******************************************************************
Double_t StatUtils::GetLikelihoodFromDiag(TH1D *data, TH1D *mc, TH1I *mask) {
//*******************************************************************
// Currently just a placeholder!
(void)data;
(void)mc;
(void)mask;
return 0.0;
};
//*******************************************************************
Double_t StatUtils::GetLikelihoodFromDiag(TH2D *data, TH2D *mc, TH2I *map,
TH2I *mask) {
//*******************************************************************
// Generate a simple map
if (!map)
map = StatUtils::GenerateMap(data);
// Convert to 1D Histograms
TH1D *data_1D = MapToTH1D(data, map);
TH1D *mc_1D = MapToTH1D(mc, map);
TH1I *mask_1D = MapToMask(mask, map);
// Calculate from 1D
Double_t MLE = StatUtils::GetLikelihoodFromDiag(data_1D, mc_1D, mask_1D);
// CleanUp
delete data_1D;
delete mc_1D;
delete mask_1D;
return MLE;
};
//*******************************************************************
Double_t StatUtils::GetLikelihoodFromCov(TH1D *data, TH1D *mc,
TMatrixDSym *invcov, TH1I *mask) {
//*******************************************************************
// Currently just a placeholder !
(void)data;
(void)mc;
(void)invcov;
(void)mask;
return 0.0;
};
//*******************************************************************
Double_t StatUtils::GetLikelihoodFromCov(TH2D *data, TH2D *mc,
TMatrixDSym *invcov, TH2I *map,
TH2I *mask) {
//*******************************************************************
// Generate a simple map
if (!map)
map = StatUtils::GenerateMap(data);
// Convert to 1D Histograms
TH1D *data_1D = MapToTH1D(data, map);
TH1D *mc_1D = MapToTH1D(mc, map);
TH1I *mask_1D = MapToMask(mask, map);
// Calculate from 1D
Double_t MLE =
StatUtils::GetLikelihoodFromCov(data_1D, mc_1D, invcov, mask_1D);
// CleanUp
delete data_1D;
delete mc_1D;
delete mask_1D;
return MLE;
};
//*******************************************************************
Double_t StatUtils::GetLikelihoodFromSVD(TH1D *data, TH1D *mc, TMatrixDSym *cov,
TH1I *mask) {
//*******************************************************************
// Currently just a placeholder!
(void)data;
(void)mc;
(void)cov;
(void)mask;
return 0.0;
};
//*******************************************************************
Double_t StatUtils::GetLikelihoodFromSVD(TH2D *data, TH2D *mc, TMatrixDSym *cov,
TH2I *map, TH2I *mask) {
//*******************************************************************
// Generate a simple map
if (!map)
map = StatUtils::GenerateMap(data);
// Convert to 1D Histograms
TH1D *data_1D = MapToTH1D(data, map);
TH1D *mc_1D = MapToTH1D(mc, map);
TH1I *mask_1D = MapToMask(mask, map);
// Calculate from 1D
Double_t MLE = StatUtils::GetLikelihoodFromSVD(data_1D, mc_1D, cov, mask_1D);
// CleanUp
delete data_1D;
delete mc_1D;
delete mask_1D;
return MLE;
};
//*******************************************************************
Double_t StatUtils::GetLikelihoodFromEventRate(TH1D *data, TH1D *mc,
TH1I *mask) {
//*******************************************************************
// Currently just a placeholder!
(void)data;
(void)mc;
(void)mask;
return 0.0;
};
//*******************************************************************
Double_t StatUtils::GetLikelihoodFromEventRate(TH2D *data, TH2D *mc, TH2I *map,
TH2I *mask) {
//*******************************************************************
// Generate a simple map
if (!map)
map = StatUtils::GenerateMap(data);
// Convert to 1D Histograms
TH1D *data_1D = MapToTH1D(data, map);
TH1D *mc_1D = MapToTH1D(mc, map);
TH1I *mask_1D = MapToMask(mask, map);
// Calculate from 1D
Double_t MLE = StatUtils::GetChi2FromEventRate(data_1D, mc_1D, mask_1D);
// CleanUp
delete data_1D;
delete mc_1D;
delete mask_1D;
return MLE;
};
//*******************************************************************
Int_t StatUtils::GetNDOF(TH1D *hist, TH1I *mask) {
//*******************************************************************
TH1D *calc_hist = (TH1D *)hist->Clone();
// If a mask is provided we need to apply it before getting NDOF
if (mask) {
calc_hist = StatUtils::ApplyHistogramMasking(hist, mask);
}
// NDOF is defined as total number of bins with non-zero errors
Int_t NDOF = 0;
for (int i = 0; i < calc_hist->GetNbinsX(); i++) {
if (calc_hist->GetBinError(i + 1) > 0.0)
NDOF++;
}
delete calc_hist;
return NDOF;
};
//*******************************************************************
Int_t StatUtils::GetNDOF(TH2D *hist, TH2I *map, TH2I *mask) {
//*******************************************************************
Int_t NDOF = 0;
if (!map)
map = StatUtils::GenerateMap(hist);
for (int i = 0; i < hist->GetNbinsX(); i++) {
for (int j = 0; j < hist->GetNbinsY(); j++) {
if (mask->GetBinContent(i + 1, j + 1))
continue;
if (map->GetBinContent(i + 1, j + 1) <= 0)
continue;
NDOF++;
}
}
return NDOF;
};
//*******************************************************************
TH1D *StatUtils::ThrowHistogram(TH1D *hist, TMatrixDSym *cov, bool throwdiag,
TH1I *mask) {
//*******************************************************************
TH1D *calc_hist =
(TH1D *)hist->Clone((std::string(hist->GetName()) + "_THROW").c_str());
TMatrixDSym *calc_cov = (TMatrixDSym *)cov->Clone();
Double_t correl_val = 0.0;
// If a mask if applied we need to apply it before the matrix is decomposed
if (mask) {
calc_cov = ApplyMatrixMasking(cov, mask);
calc_hist = ApplyHistogramMasking(calc_hist, mask);
}
// If a covariance is provided we need a preset random vector and a decomp
std::vector<Double_t> rand_val;
TMatrixDSym *decomp_cov = NULL;
if (cov) {
for (int i = 0; i < hist->GetNbinsX(); i++) {
rand_val.push_back(gRandom->Gaus(0.0, 1.0));
}
// Decomp the matrix
decomp_cov = StatUtils::GetDecomp(calc_cov);
}
// iterate over bins
for (int i = 0; i < hist->GetNbinsX(); i++) {
// By Default the errors on the histogram are thrown uncorrelated to the
// other errors
/*
if (throwdiag) {
calc_hist->SetBinContent(i + 1, (calc_hist->GetBinContent(i + 1) + \
gRandom->Gaus(0.0, 1.0) * calc_hist->GetBinError(i + 1)) );
}
*/
// If a covariance is provided that is also thrown
if (cov) {
correl_val = 0.0;
for (int j = 0; j < hist->GetNbinsX(); j++) {
correl_val += rand_val[j] * (*decomp_cov)(j, i);
}
calc_hist->SetBinContent(
i + 1, (calc_hist->GetBinContent(i + 1) + correl_val * 1E-38));
}
}
delete calc_cov;
delete decomp_cov;
// return this new thrown data
return calc_hist;
};
//*******************************************************************
TH2D *StatUtils::ThrowHistogram(TH2D *hist, TMatrixDSym *cov, TH2I *map,
bool throwdiag, TH2I *mask) {
//*******************************************************************
// PLACEHOLDER!!!!!!!!!
// Currently no support for throwing 2D Histograms from a covariance
(void)hist;
(void)cov;
(void)map;
(void)throwdiag;
(void)mask;
// /todo
// Sort maps if required
// Throw the covariance for a 1D plot
// Unmap back to 2D Histogram
return hist;
}
//*******************************************************************
TH1D *StatUtils::ApplyHistogramMasking(TH1D *hist, TH1I *mask) {
//*******************************************************************
if (!mask)
return ((TH1D *)hist->Clone());
// This masking is only sufficient for chi2 calculations, and will have dodgy
// bin edges.
// Get New Bin Count
Int_t NBins = 0;
for (int i = 0; i < hist->GetNbinsX(); i++) {
if (mask->GetBinContent(i + 1))
continue;
NBins++;
}
// Make new hist
std::string newmaskname = std::string(hist->GetName()) + "_MSKD";
TH1D *calc_hist =
new TH1D(newmaskname.c_str(), newmaskname.c_str(), NBins, 0, NBins);
// fill new hist
int binindex = 0;
for (int i = 0; i < hist->GetNbinsX(); i++) {
if (mask->GetBinContent(i + 1)) {
NUIS_LOG(REC, "Applying mask to bin " << i + 1 << " " << hist->GetName());
continue;
}
calc_hist->SetBinContent(binindex + 1, hist->GetBinContent(i + 1));
calc_hist->SetBinError(binindex + 1, hist->GetBinError(i + 1));
binindex++;
}
return calc_hist;
};
//*******************************************************************
TH2D *StatUtils::ApplyHistogramMasking(TH2D *hist, TH2I *mask) {
//*******************************************************************
TH2D *newhist = (TH2D *)hist->Clone();
if (!mask)
return newhist;
for (int i = 0; i < hist->GetNbinsX(); i++) {
for (int j = 0; j < hist->GetNbinsY(); j++) {
if (mask->GetBinContent(i + 1, j + 1) > 0) {
newhist->SetBinContent(i + 1, j + 1, 0.0);
newhist->SetBinContent(i + 1, j + 1, 0.0);
}
}
}
return newhist;
}
//*******************************************************************
TMatrixDSym *StatUtils::ApplyMatrixMasking(TMatrixDSym *mat, TH1I *mask) {
//*******************************************************************
if (!mask)
return (TMatrixDSym *)(mat->Clone());
// Get New Bin Count
Int_t NBins = 0;
for (int i = 0; i < mask->GetNbinsX(); i++) {
if (mask->GetBinContent(i + 1))
continue;
NBins++;
}
// make new matrix
TMatrixDSym *calc_mat = new TMatrixDSym(NBins);
int col, row;
// Need to mask out bins in the current matrix
row = 0;
for (int i = 0; i < mask->GetNbinsX(); i++) {
col = 0;
// skip if masked
if (mask->GetBinContent(i + 1) > 0.5)
continue;
for (int j = 0; j < mask->GetNbinsX(); j++) {
// skip if masked
if (mask->GetBinContent(j + 1) > 0.5)
continue;
(*calc_mat)(row, col) = (*mat)(i, j);
col++;
}
row++;
}
return calc_mat;
};
//*******************************************************************
TMatrixDSym *StatUtils::ApplyMatrixMasking(TMatrixDSym *mat, TH2D *data,
TH2I *mask, TH2I *map) {
//*******************************************************************
if (!map)
map = StatUtils::GenerateMap(data);
TH1I *mask_1D = StatUtils::MapToMask(mask, map);
TMatrixDSym *newmat = StatUtils::ApplyMatrixMasking(mat, mask_1D);
delete mask_1D;
return newmat;
}
//*******************************************************************
TMatrixDSym *StatUtils::ApplyInvertedMatrixMasking(TMatrixDSym *mat,
TH1I *mask) {
//*******************************************************************
TMatrixDSym *new_mat = GetInvert(mat);
TMatrixDSym *masked_mat = ApplyMatrixMasking(new_mat, mask);
TMatrixDSym *inverted_mat = GetInvert(masked_mat);
delete masked_mat;
delete new_mat;
return inverted_mat;
};
//*******************************************************************
TMatrixDSym *StatUtils::ApplyInvertedMatrixMasking(TMatrixDSym *mat, TH2D *data,
TH2I *mask, TH2I *map) {
//*******************************************************************
if (!map)
map = StatUtils::GenerateMap(data);
TH1I *mask_1D = StatUtils::MapToMask(mask, map);
TMatrixDSym *newmat = ApplyInvertedMatrixMasking(mat, mask_1D);
delete mask_1D;
return newmat;
}
//*******************************************************************
TMatrixDSym *StatUtils::GetInvert(TMatrixDSym *mat) {
//*******************************************************************
TMatrixDSym *new_mat = (TMatrixDSym *)mat->Clone();
// Check for diagonal
bool non_diagonal = false;
for (int i = 0; i < new_mat->GetNrows(); i++) {
for (int j = 0; j < new_mat->GetNrows(); j++) {
if (i == j)
continue;
if ((*new_mat)(i, j) != 0.0) {
non_diagonal = true;
break;
}
}
}
// If diag, just flip the diag
if (!non_diagonal or new_mat->GetNrows() == 1) {
for (int i = 0; i < new_mat->GetNrows(); i++) {
if ((*new_mat)(i, i) != 0.0)
(*new_mat)(i, i) = 1.0 / (*new_mat)(i, i);
else
(*new_mat)(i, i) = 0.0;
}
return new_mat;
}
// Invert full matrix
TDecompSVD LU = TDecompSVD((*new_mat));
new_mat =
new TMatrixDSym(new_mat->GetNrows(), LU.Invert().GetMatrixArray(), "");
return new_mat;
}
//*******************************************************************
TMatrixDSym *StatUtils::GetDecomp(TMatrixDSym *mat) {
//*******************************************************************
TMatrixDSym *new_mat = (TMatrixDSym *)mat->Clone();
int nrows = new_mat->GetNrows();
// Check for diagonal
bool diagonal = true;
for (int i = 0; i < nrows; i++) {
for (int j = 0; j < nrows; j++) {
if (i == j)
continue;
if ((*new_mat)(i, j) != 0.0) {
diagonal = false;
break;
}
}
}
// If diag, just flip the diag
if (diagonal or nrows == 1) {
for (int i = 0; i < nrows; i++) {
if ((*new_mat)(i, i) > 0.0)
(*new_mat)(i, i) = sqrt((*new_mat)(i, i));
else
(*new_mat)(i, i) = 0.0;
}
return new_mat;
}
TDecompChol LU = TDecompChol(*new_mat);
LU.Decompose();
delete new_mat;
TMatrixDSym *dec_mat = new TMatrixDSym(nrows, LU.GetU().GetMatrixArray(), "");
return dec_mat;
}
//*******************************************************************
void StatUtils::ForceNormIntoCovar(TMatrixDSym *&mat, TH1D *hist, double norm) {
//*******************************************************************
if (!mat)
mat = MakeDiagonalCovarMatrix(hist);
int nbins = mat->GetNrows();
TMatrixDSym *new_mat = new TMatrixDSym(nbins);
for (int i = 0; i < nbins; i++) {
for (int j = 0; j < nbins; j++) {
double valx = hist->GetBinContent(i + 1) * 1E38;
double valy = hist->GetBinContent(j + 1) * 1E38;
(*new_mat)(i, j) = (*mat)(i, j) + norm * norm * valx * valy;
}
}
// Swap the two
delete mat;
mat = new_mat;
return;
};
//*******************************************************************
void StatUtils::ForceNormIntoCovar(TMatrixDSym *mat, TH2D *data, double norm,
TH2I *map) {
//*******************************************************************
if (!map)
map = StatUtils::GenerateMap(data);
TH1D *data_1D = MapToTH1D(data, map);
StatUtils::ForceNormIntoCovar(mat, data_1D, norm);
delete data_1D;
return;
}
//*******************************************************************
TMatrixDSym *StatUtils::MakeDiagonalCovarMatrix(TH1D *data, double scaleF) {
//*******************************************************************
TMatrixDSym *newmat = new TMatrixDSym(data->GetNbinsX());
for (int i = 0; i < data->GetNbinsX(); i++) {
(*newmat)(i, i) =
data->GetBinError(i + 1) * data->GetBinError(i + 1) * scaleF * scaleF;
}
return newmat;
}
//*******************************************************************
TMatrixDSym *StatUtils::MakeDiagonalCovarMatrix(TH2D *data, TH2I *map,
double scaleF) {
//*******************************************************************
if (!map)
map = StatUtils::GenerateMap(data);
TH1D *data_1D = MapToTH1D(data, map);
return StatUtils::MakeDiagonalCovarMatrix(data_1D, scaleF);
};
//*******************************************************************
void StatUtils::SetDataErrorFromCov(TH1D *DataHist, TMatrixDSym *cov,
double scale, bool ErrorCheck) {
//*******************************************************************
// Check
if (ErrorCheck) {
if (cov->GetNrows() != DataHist->GetNbinsX()) {
NUIS_ERR(
FTL,
"Nrows in cov don't match nbins in DataHist for SetDataErrorFromCov");
NUIS_ERR(FTL, "Nrows = " << cov->GetNrows());
NUIS_ABORT("Nbins = " << DataHist->GetNbinsX());
}
}
// Set bin errors form cov diag
// Check if the errors are set
bool ErrorsSet = false;
for (int i = 0; i < DataHist->GetNbinsX(); i++) {
if (ErrorsSet == true)
break;
if (DataHist->GetBinError(i + 1) != 0 && DataHist->GetBinContent(i + 1) > 0)
ErrorsSet = true;
}
// Now loop over
if (ErrorsSet && ErrorCheck) {
for (int i = 0; i < DataHist->GetNbinsX(); i++) {
double DataHisterr = DataHist->GetBinError(i + 1);
double coverr = sqrt((*cov)(i, i)) * scale;
// Check that the errors are within 1% of eachother
if (fabs(DataHisterr - coverr) / DataHisterr > 0.01) {
NUIS_ERR(WRN, "Data error does not match covariance error for bin "
<< i + 1 << " ("
<< DataHist->GetXaxis()->GetBinLowEdge(i + 1) << "-"
<< DataHist->GetXaxis()->GetBinLowEdge(i + 2) << ")");
NUIS_ERR(WRN, "Data error: " << DataHisterr);
NUIS_ERR(WRN, "Cov error: " << coverr);
}
}
// Else blindly trust the covariance
} else {
for (int i = 0; i < DataHist->GetNbinsX(); i++) {
DataHist->SetBinError(i + 1, sqrt((*cov)(i, i)) * scale);
}
}
return;
}
//*******************************************************************
void StatUtils::SetDataErrorFromCov(TH2D *data, TMatrixDSym *cov, TH2I *map,
double scale, bool ErrorCheck) {
//*******************************************************************
// Check
if (ErrorCheck) {
if (cov->GetNrows() != data->GetNbinsX() * data->GetNbinsY()) {
NUIS_ERR(
FTL,
"Nrows in cov don't match nbins in data for SetDataNUIS_ERRorFromCov");
NUIS_ERR(FTL, "Nrows = " << cov->GetNrows());
NUIS_ABORT("Nbins = " << data->GetNbinsX());
}
}
// Set bin errors form cov diag
// Check if the errors are set
bool ErrorsSet = false;
for (int i = 0; i < data->GetNbinsX(); i++) {
for (int j = 0; j < data->GetNbinsX(); j++) {
if (ErrorsSet == true)
break;
if (data->GetBinError(i + 1, j + 1) != 0)
ErrorsSet = true;
}
}
// Create map if required
if (!map)
map = StatUtils::GenerateMap(data);
// Set Bin Errors from cov diag
int count = 0;
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)
continue;
// If we have errors on our histogram the map is good
count = map->GetBinContent(i + 1, j + 1) - 1;
double dataerr = data->GetBinError(i + 1, j + 1);
double coverr = sqrt((*cov)(count, count)) * scale;
// Check that the errors are within 1% of eachother
if (ErrorsSet && ErrorCheck) {
if (fabs(dataerr - coverr) / dataerr > 0.01) {
NUIS_ERR(WRN, "Data error does not match covariance error for bin "
<< i + 1 << " ("
<< data->GetXaxis()->GetBinLowEdge(i + 1) << "-"
<< data->GetXaxis()->GetBinLowEdge(i + 2) << ")");
NUIS_ERR(WRN, "Data error: " << dataerr);
NUIS_ERR(WRN, "Cov error: " << coverr);
}
} else {
data->SetBinError(i + 1, j + 1, sqrt((*cov)(count, count)) * scale);
}
}
}
// Delete the map now that we don't need it
// Woops, it's needed elsewhere there! (Grrrrr)
// map->Delete();
}
TMatrixDSym *StatUtils::ExtractShapeOnlyCovar(TMatrixDSym *full_covar,
TH1 *data_hist,
double data_scale) {
int nbins = full_covar->GetNrows();
TMatrixDSym *shape_covar = new TMatrixDSym(nbins);
// Check nobody is being silly
if (data_hist->GetNbinsX() != nbins) {
NUIS_ERR(WRN, "Inconsistent matrix and data histogram passed to "
"StatUtils::ExtractShapeOnlyCovar!");
NUIS_ERR(WRN, "data_hist has " << data_hist->GetNbinsX() << " matrix has "
<< nbins);
int err_bins = data_hist->GetNbinsX();
if (nbins > err_bins)
err_bins = nbins;
for (int i = 0; i < err_bins; ++i) {
NUIS_ERR(WRN, "Matrix diag. = " << (*full_covar)(i, i) << " data = "
<< data_hist->GetBinContent(i + 1));
}
return NULL;
}
double total_data = 0;
double total_covar = 0;
// Initial loop to calculate some constants
for (int i = 0; i < nbins; ++i) {
total_data += data_hist->GetBinContent(i + 1) * data_scale;
for (int j = 0; j < nbins; ++j) {
total_covar += (*full_covar)(i, j);
}
}
if (total_data == 0 || total_covar == 0) {
NUIS_ERR(WRN, "Stupid matrix or data histogram passed to "
"StatUtils::ExtractShapeOnlyCovar! Ignoring...");
return NULL;
}
NUIS_LOG(SAM, "Norm error = " << sqrt(total_covar) / total_data);
// Now loop over and calculate the shape-only matrix
for (int i = 0; i < nbins; ++i) {
double data_i = data_hist->GetBinContent(i + 1) * data_scale;
for (int j = 0; j < nbins; ++j) {
double data_j = data_hist->GetBinContent(j + 1) * data_scale;
double norm_term =
data_i * data_j * total_covar / total_data / total_data;
double mix_sum1 = 0;
double mix_sum2 = 0;
for (int k = 0; k < nbins; ++k) {
mix_sum1 += (*full_covar)(k, j);
mix_sum2 += (*full_covar)(i, k);
}
double mix_term1 =
data_i * (mix_sum1 / total_data -
total_covar * data_j / total_data / total_data);
double mix_term2 =
data_j * (mix_sum2 / total_data -
total_covar * data_i / total_data / total_data);
(*shape_covar)(i, j) =
(*full_covar)(i, j) - mix_term1 - mix_term2 - norm_term;
}
}
return shape_covar;
}
//*******************************************************************
TH2I *StatUtils::GenerateMap(TH2D *hist) {
//*******************************************************************
std::string maptitle = std::string(hist->GetName()) + "_MAP";
TH2I *map =
new TH2I(maptitle.c_str(), maptitle.c_str(), hist->GetNbinsX(), 0,
hist->GetNbinsX(), hist->GetNbinsY(), 0, hist->GetNbinsY());
Int_t index = 1;
for (int i = 0; i < hist->GetNbinsX(); i++) {
for (int j = 0; j < hist->GetNbinsY(); j++) {
if (hist->GetBinContent(i + 1, j + 1) > 0) {
map->SetBinContent(i + 1, j + 1, index);
index++;
} else {
map->SetBinContent(i + 1, j + 1, 0);
}
}
}
return map;
}
//*******************************************************************
TH1D *StatUtils::MapToTH1D(TH2D *hist, TH2I *map) {
//*******************************************************************
if (!hist)
return NULL;
// Get N bins for 1D plot
Int_t Nbins = map->GetMaximum();
std::string name1D = std::string(hist->GetName()) + "_1D";
// Make new 1D Hist
TH1D *newhist = new TH1D(name1D.c_str(), name1D.c_str(), Nbins, 0, Nbins);
// map bin contents
for (int i = 0; i < map->GetNbinsX(); i++) {
for (int j = 0; j < map->GetNbinsY(); j++) {
if (map->GetBinContent(i + 1, j + 1) == 0)
continue;
newhist->SetBinContent(map->GetBinContent(i + 1, j + 1),
hist->GetBinContent(i + 1, j + 1));
newhist->SetBinError(map->GetBinContent(i + 1, j + 1),
hist->GetBinError(i + 1, j + 1));
}
}
// return
return newhist;
}
//*******************************************************************
TH1I *StatUtils::MapToMask(TH2I *hist, TH2I *map) {
//*******************************************************************
TH1I *newhist = NULL;
if (!hist)
return newhist;
// Get N bins for 1D plot
Int_t Nbins = map->GetMaximum();
std::string name1D = std::string(hist->GetName()) + "_1D";
// Make new 1D Hist
newhist = new TH1I(name1D.c_str(), name1D.c_str(), Nbins, 0, Nbins);
// map bin contents
for (int i = 0; i < map->GetNbinsX(); i++) {
for (int j = 0; j < map->GetNbinsY(); j++) {
if (map->GetBinContent(i + 1, j + 1) == 0)
continue;
newhist->SetBinContent(map->GetBinContent(i + 1, j + 1),
hist->GetBinContent(i + 1, j + 1));
}
}
// return
return newhist;
}
TMatrixDSym *StatUtils::GetCovarFromCorrel(TMatrixDSym *correl, TH1D *data) {
int nbins = correl->GetNrows();
TMatrixDSym *covar = new TMatrixDSym(nbins);
for (int i = 0; i < nbins; i++) {
for (int j = 0; j < nbins; j++) {
(*covar)(i, j) =
(*correl)(i, j) * data->GetBinError(i + 1) * data->GetBinError(j + 1);
}
}
return covar;
}
//*******************************************************************
TMatrixD *StatUtils::GetMatrixFromTextFile(std::string covfile, int dimx,
int dimy) {
//*******************************************************************
// Determine dim
if (dimx == -1 and dimy == -1) {
std::string line;
std::ifstream covar(covfile.c_str(), std::ifstream::in);
int row = 0;
while (std::getline(covar >> std::ws, line, '\n')) {
int column = 0;
std::vector<double> entries = GeneralUtils::ParseToDbl(line, " ");
if (entries.size() <= 1) {
NUIS_ERR(WRN, "StatUtils::GetMatrixFromTextFile, matrix only has <= 1 "
"entries on this line: "
<< row);
}
for (std::vector<double>::iterator iter = entries.begin();
iter != entries.end(); iter++) {
column++;
if (column > dimx)
dimx = column;
}
row++;
if (row > dimy)
dimy = row;
}
}
// Or assume symmetric
if (dimx != -1 and dimy == -1) {
dimy = dimx;
}
assert(dimy != -1 && " matrix dimy not set.");
// Make new matrix
TMatrixD *mat = new TMatrixD(dimx, dimy);
std::string line;
std::ifstream covar(covfile.c_str(), std::ifstream::in);
int row = 0;
while (std::getline(covar >> std::ws, line, '\n')) {
int column = 0;
std::vector<double> entries = GeneralUtils::ParseToDbl(line, " ");
if (entries.size() <= 1) {
NUIS_ERR(WRN, "StatUtils::GetMatrixFromTextFile, matrix only has <= 1 "
"entries on this line: "
<< row);
}
for (std::vector<double>::iterator iter = entries.begin();
iter != entries.end(); iter++) {
// Check Rows
// assert(row > mat->GetNrows() && " covar rows doesn't match matrix
// rows.");
// assert(column > mat->GetNcols() && " covar cols doesn't match matrix
// cols.");
// Fill Matrix
(*mat)(row, column) = (*iter);
column++;
}
row++;
}
return mat;
}
//*******************************************************************
TMatrixD *StatUtils::GetMatrixFromRootFile(std::string covfile,
std::string histname) {
//*******************************************************************
std::string inputfile = covfile + ";" + histname;
std::vector<std::string> splitfile = GeneralUtils::ParseToStr(inputfile, ";");
if (splitfile.size() < 2) {
NUIS_ABORT("No object name given!");
}
// Get file
TFile *tempfile = new TFile(splitfile[0].c_str(), "READ");
// Get Object
TObject *obj = tempfile->Get(splitfile[1].c_str());
if (!obj) {
NUIS_ABORT("Object " << splitfile[1] << " doesn't exist!");
}
// Try casting
TMatrixD *mat = dynamic_cast<TMatrixD *>(obj);
if (mat) {
TMatrixD *newmat = (TMatrixD *)mat->Clone();
delete mat;
tempfile->Close();
return newmat;
}
TMatrixDSym *matsym = dynamic_cast<TMatrixDSym *>(obj);
if (matsym) {
TMatrixD *newmat = new TMatrixD(matsym->GetNrows(), matsym->GetNrows());
for (int i = 0; i < matsym->GetNrows(); i++) {
for (int j = 0; j < matsym->GetNrows(); j++) {
(*newmat)(i, j) = (*matsym)(i, j);
}
}
delete matsym;
tempfile->Close();
return newmat;
}
TH2D *mathist = dynamic_cast<TH2D *>(obj);
if (mathist) {
TMatrixD *newmat = new TMatrixD(mathist->GetNbinsX(), mathist->GetNbinsX());
for (int i = 0; i < mathist->GetNbinsX(); i++) {
for (int j = 0; j < mathist->GetNbinsX(); j++) {
(*newmat)(i, j) = mathist->GetBinContent(i + 1, j + 1);
}
}
delete mathist;
tempfile->Close();
return newmat;
}
return NULL;
}
//*******************************************************************
TMatrixDSym *StatUtils::GetCovarFromTextFile(std::string covfile, int dim) {
//*******************************************************************
// Delete TempMat
TMatrixD *tempmat = GetMatrixFromTextFile(covfile, dim, dim);
// Make a symmetric covariance
TMatrixDSym *newmat = new TMatrixDSym(tempmat->GetNrows());
for (int i = 0; i < tempmat->GetNrows(); i++) {
for (int j = 0; j < tempmat->GetNrows(); j++) {
(*newmat)(i, j) = (*tempmat)(i, j);
}
}
delete tempmat;
return newmat;
}
//*******************************************************************
TMatrixDSym *StatUtils::GetCovarFromRootFile(std::string covfile,
std::string histname) {
//*******************************************************************
TMatrixD *tempmat = GetMatrixFromRootFile(covfile, histname);
TMatrixDSym *newmat = new TMatrixDSym(tempmat->GetNrows());
for (int i = 0; i < tempmat->GetNrows(); i++) {
for (int j = 0; j < tempmat->GetNrows(); j++) {
(*newmat)(i, j) = (*tempmat)(i, j);
}
}
delete tempmat;
return newmat;
}
diff --git a/src/Statistical/StatUtils.h b/src/Statistical/StatUtils.h
index 435e6a7..d2f82c9 100644
--- a/src/Statistical/StatUtils.h
+++ b/src/Statistical/StatUtils.h
@@ -1,255 +1,255 @@
// 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#ifndef STATUTILS_H
#define STATUTILS_H
// C Includes
#include <stdlib.h>
#include <numeric>
#include <math.h>
#include <string>
#include <iostream>
#include <sstream>
#include <iomanip>
#include <deque>
#include "assert.h"
// Root Includes
#include "TH1D.h"
#include "TH2I.h"
#include "TH2D.h"
#include "TFile.h"
#include "TMatrixDSym.h"
#include "TDecompSVD.h"
#include "TMath.h"
#include "TRandom3.h"
#include "TDecompChol.h"
#include "TGraphErrors.h"
// Fit Includes
#include "FitLogger.h"
/*!
* \addtogroup Utils
* @{
*/
//! Functions for handling statistics calculations
namespace StatUtils{
/*
Chi2 Functions
*/
//! Get Chi2 using diagonal bin errors from the histogram. Masking applied before calculation if mask provided.
Double_t GetChi2FromDiag(TH1D* data, TH1D* mc, TH1I* mask=NULL);
//! Get Chi2 using diagonal bin errors from the histogram.
//! Plots converted to 1D histograms before using 1D calculation.
Double_t GetChi2FromDiag(TH2D* data, TH2D* mc, TH2I* map=NULL, TH2I* mask=NULL);
//! Get Chi2 using an inverted covariance for the data
- Double_t GetChi2FromCov( TH1D* data, TH1D* mc, TMatrixDSym* invcov, TH1I* mask=NULL, double data_scale=1, double covar_scale=1E76);
+ Double_t GetChi2FromCov( TH1D* data, TH1D* mc, TMatrixDSym* invcov, TH1I* mask=NULL, double data_scale=1, double covar_scale=1E76, TH1D *outchi2perbin=NULL);
//! Get Chi2 using an inverted covariance for the data
//! Plots converted to 1D histograms before using 1D calculation.
- Double_t GetChi2FromCov( TH2D* data, TH2D* mc, TMatrixDSym* invcov, TH2I* map=NULL, TH2I* mask=NULL);
+ Double_t GetChi2FromCov( TH2D* data, TH2D* mc, TMatrixDSym* invcov, TH2I* map=NULL, TH2I* mask=NULL, TH2D *outchi2perbin=NULL);
//! Get Chi2 using an SVD method on the covariance before calculation.
//! Method suggested by Rex at MiniBooNE. Shown that it doesn't actually work.
Double_t GetChi2FromSVD( TH1D* data, TH1D* mc, TMatrixDSym* cov, TH1I* mask=NULL);
//! Get Chi2 using an SVD method on the covariance before calculation.
//! Method suggested by Rex at MiniBooNE. Shown that it doesn't actually work.
//! Plots converted to 1D histograms before using 1D calculation.
Double_t GetChi2FromSVD( TH2D* data, TH2D* mc, TMatrixDSym* cov, TH2I* map=NULL, TH2I* mask=NULL);
//! Get Chi2 using only the raw event rates given in each bin using a -2LL method.
Double_t GetChi2FromEventRate(TH1D* data, TH1D* mc, TH1I* mask=NULL);
//! Get Chi2 using only the raw event rates given in each bin using a -2LL method.
//! Plots converted to 1D histograms before using 1D calculation.
Double_t GetChi2FromEventRate(TH2D* data, TH2D* mc, TH2I* map=NULL, TH2I* mask=NULL);
// Likelihood Functions
//! Placeholder for 1D binned likelihood method
Double_t GetLikelihoodFromDiag(TH1D* data, TH1D* mc, TH1I* mask=NULL);
//! Placeholder for 2D binned likelihood method
Double_t GetLikelihoodFromDiag(TH2D* data, TH2D* mc, TH2I* map=NULL, TH2I* mask=NULL);
//! Placeholder for 1D binned likelihood method
Double_t GetLikelihoodFromCov( TH1D* data, TH1D* mc, TMatrixDSym* invcov, TH1I* mask=NULL);
//! Placeholder for 2D binned likelihood method
Double_t GetLikelihoodFromCov( TH2D* data, TH2D* mc, TMatrixDSym* invcov, TH2I* map=NULL, TH2I* mask=NULL);
//! Placeholder for 1D binned likelihood method
Double_t GetLikelihoodFromSVD( TH1D* data, TH1D* mc, TMatrixDSym* cov, TH1I* mask=NULL);
//! Placeholder for 2D binned likelihood method
Double_t GetLikelihoodFromSVD( TH2D* data, TH2D* mc, TMatrixDSym* cov, TH2I* map=NULL, TH2I* mask=NULL);
//! Placeholder for 1D binned likelihood method
Double_t GetLikelihoodFromEventRate(TH1D* data, TH1D* mc, TH1I* mask=NULL);
//! Placeholder for 2D binned likelihood method
Double_t GetLikelihoodFromEventRate(TH2D* data, TH2D* mc, TH2I* map=NULL, TH2I* mask=NULL);
/*
NDOF Functions
*/
//! Return 1D Histogram NDOF considering masking and empty bins
Int_t GetNDOF(TH1D* hist, TH1I* mask=NULL);
//! Return 2D Histogram NDOF considering masking and empty bins
Int_t GetNDOF(TH2D* hist, TH2I* map=NULL, TH2I* mask=NULL);
/*
Fake Data Functions
*/
//! Given a full covariance for a 1D data set throw the decomposition to generate fake data.
//! throwdiag determines whether diagonal statistical errors are thrown.
//! If no covariance is provided only statistical errors are thrown.
TH1D* ThrowHistogram(TH1D* hist, TMatrixDSym* cov, bool throwdiag=true, TH1I* mask=NULL);
//! Given a full covariance for a 2D data set throw the decomposition to generate fake data.
//! Plots are converted to 1D histograms and the 1D ThrowHistogram is used, before being converted back to 2D histograms.
TH2D* ThrowHistogram(TH2D* hist, TMatrixDSym* cov, TH2I* map=NULL, bool throwdiag=true, TH2I* mask=NULL);
/*
Masking Functions
*/
//! Given a mask histogram, mask out any bins in hist with non zero entries in mask.
TH1D* ApplyHistogramMasking(TH1D* hist, TH1I* mask);
//! Given a mask histogram, mask out any bins in hist with non zero entries in mask.
TH2D* ApplyHistogramMasking(TH2D* hist, TH2I* mask);
//! Given a mask histogram apply the masking procedure to each of the rows/columns in a covariance, before recalculating its inverse.
TMatrixDSym* ApplyInvertedMatrixMasking(TMatrixDSym* mat, TH1I* mask);
//! Given a mask histogram apply the masking procedure to each of the rows/columns in a covariance, before recalculating its inverse.
//! Converts to 1D data before using the 1D ApplyInvertedMatrixMasking function and converting back to 2D.
TMatrixDSym* ApplyInvertedMatrixMasking(TMatrixDSym* mat, TH2D* data, TH2I* mask, TH2I* map=NULL);
//! Given a mask histogram apply the masking procedure to each of the rows/columns in a covariance
TMatrixDSym* ApplyMatrixMasking(TMatrixDSym* mat, TH1I* mask);
//! Given a mask histogram apply the masking procedure to each of the rows/columns in a covariance
//! Converts to 1D data before using the 1D ApplyInvertedMatrixMasking function and converting back to 2D.
TMatrixDSym* ApplyMatrixMasking(TMatrixDSym* mat, TH2D* data, TH2I* mask, TH2I* map=NULL);
/*
Covariance Handling Functions
*/
//! Return inverted matrix of TMatrixDSym
TMatrixDSym* GetInvert(TMatrixDSym* mat);
//! Return Cholesky Decomposed matrix of TMatrixDSym
TMatrixDSym* GetDecomp(TMatrixDSym* mat);
//! Return full covariances
TMatrixDSym* GetCovarFromCorrel(TMatrixDSym* correl, TH1D* data);
//! Given a normalisation factor for a dataset add in a new normalisation term to the covariance.
void ForceNormIntoCovar(TMatrixDSym*& mat, TH1D* data, double norm);
//! Given a normalisation factor for a dataset add in a new normalisation term to the covariance.
//! Convertes 2D to 1D, before using 1D ForceNormIntoCovar
void ForceNormIntoCovar(TMatrixDSym* mat, TH2D* data, double norm, TH2I* map=NULL);
//! Given a dataset generate an uncorrelated covariance matrix using the bin errors.
TMatrixDSym* MakeDiagonalCovarMatrix(TH1D* data, double scaleF=1E38);
//! Given a dataset generate an uncorrelated covariance matrix using the bin errors.
TMatrixDSym* MakeDiagonalCovarMatrix(TH2D* data, TH2I* map=NULL, double scaleF=1E38);
//! Given a covariance set the errors in each bin on the data from the covariance diagonals.
void SetDataErrorFromCov(TH1D* data, TMatrixDSym* cov, double scale=1.0, bool ErrorCheck = true);
//! Given a covariance set the errors in each bin on the data from the covariance diagonals.
void SetDataErrorFromCov(TH2D* data, TMatrixDSym* cov, TH2I* map=NULL, double scale=1.0, bool ErrorCheck = true);
//! Given a covariance, extracts the shape-only matrix using the method from the MiniBooNE TN
TMatrixDSym* ExtractShapeOnlyCovar(TMatrixDSym* full_covar, TH1* data_hist, double data_scale=1E38);
/*
Mapping Functions
*/
//! If no map is provided for the 2D histogram, generate one by counting up through the bins along x and y.
TH2I* GenerateMap(TH2D* hist);
//! Apply a map to a 2D histogram converting it into a 1D histogram.
TH1D* MapToTH1D(TH2D* hist, TH2I* map);
//! Apply a map to a 2D mask convering it into a 1D mask.
TH1I* MapToMask(TH2I* hist, TH2I* map);
/// \brief Read TMatrixD from a text file
///
/// - covfile = full path to text file
/// - dimx = x dimensions of matrix
/// - dimy = y dimensions of matrix
///
/// Format of textfile should be: \n
/// cov_11 cov_12 ... cov_1N \n
/// cov_21 cov_22 ... cov_2N \n
/// ... ... ... ... \n
/// cov_N1 ... ... cov_NN \n
///
/// If no dimensions are given, dimx and dimy are determined from rows/columns
/// inside textfile.
///
/// If only dimx is given a symmetric matrix is assumed.
TMatrixD* GetMatrixFromTextFile(std::string covfile, int dimx=-1, int dimy=-1);
/// \brief Read TMatrixD from a ROOT file
///
/// - covfile = full path to root file (+';histogram')
/// - histname = histogram name
///
/// If no histogram name is given function assumes it has been appended
/// covfile path as: \n
- /// 'covfile.root;histname'
+ /// 'covfile.root;histname'
///
/// histname can point to a TMatrixD object, a TMatrixDSym object, or
/// a TH2D object.
TMatrixD* GetMatrixFromRootFile(std::string covfile, std::string histname="");
/// \brief Calls GetMatrixFromTextFile and turns it into a TMatrixDSym
TMatrixDSym* GetCovarFromTextFile(std::string covfile, int dim);
/// \brief Calls GetMatrixFromRootFile and turns it into a TMatrixDSym
TMatrixDSym* GetCovarFromRootFile(std::string covfile, std::string histname);
};
/*! @} */
#endif
diff --git a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddat_nu.cxx b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddat_nu.cxx
index 0cc4155..956401e 100644
--- a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddat_nu.cxx
+++ b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddat_nu.cxx
@@ -1,87 +1,90 @@
// 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 <http://www.gnu.org/licenses/>.
-*******************************************************************************/
+ * 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 <http://www.gnu.org/licenses/>.
+ *******************************************************************************/
#include "T2K_SignalDef.h"
#include "T2K_CC0pinp_STV_XSec_1Ddat_nu.h"
-
-//********************************************************************
-T2K_CC0pinp_STV_XSec_1Ddat_nu::T2K_CC0pinp_STV_XSec_1Ddat_nu(nuiskey samplekey) {
//********************************************************************
+T2K_CC0pinp_STV_XSec_1Ddat_nu::T2K_CC0pinp_STV_XSec_1Ddat_nu(
+ nuiskey samplekey) {
+ //********************************************************************
// Sample overview ---------------------------------------------------
- std::string descrip = "T2K_CC0pinp_STV_XSec_1Ddat_nu sample. \n" \
- "Target: CH \n" \
- "Flux: T2K 2.5 degree off-axis (ND280) \n" \
- "Signal: CC0piNp (N>=1) with 450M eV < p_p < 1 GeV \n" \
- " p_mu > 250 MeV \n" \
- " cth_p > 0.6 \n" \
- " cth_mu > -0.4 \n";
+ std::string descrip = "T2K_CC0pinp_STV_XSec_1Ddat_nu sample. \n"
+ "Target: CH \n"
+ "Flux: T2K 2.5 degree off-axis (ND280) \n"
+ "Signal: CC0piNp (N>=1) with 450M eV < p_p < 1 GeV \n"
+ " p_mu > 250 MeV \n"
+ " cth_p > 0.6 \n"
+ " cth_mu > -0.4 \n"
+ "https://doi.org/10.1103/PhysRevD.98.032003 \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("#delta#it{#alpha}_{T} (GeV c^{-1})");
- fSettings.SetYTitle("#frac{d#sigma}{d#delta#it{#alpha}_{T}} (cm^{2} nucleon^{-1} rads^{-1})");
+ fSettings.SetYTitle(
+ "#frac{d#sigma}{d#delta#it{#alpha}_{T}} (cm^{2} nucleon^{-1} rads^{-1})");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG");
fSettings.SetEnuRange(0.0, 50.0);
fSettings.DefineAllowedTargets("C,H");
// CCQELike plot information
fSettings.SetTitle("T2K_CC0pinp_STV_XSec_1Ddat_nu");
- // fSettings.SetDataInput( GeneralUtils::GetTopLevelDir() + "/data/T2K/T2K_CC0pinp_STV_XSec_1Ddat_nu.dat");
+ // fSettings.SetDataInput( GeneralUtils::GetTopLevelDir() +
+ // "/data/T2K/T2K_CC0pinp_STV_XSec_1Ddat_nu.dat");
- fSettings.SetDataInput( FitPar::GetDataBase() + "/T2K/CC0pi/STV/datResults.root;Result" );
- fSettings.SetCovarInput( FitPar::GetDataBase() + "/T2K/CC0pi/STV/datResults.root;Correlation_Matrix" );
+ fSettings.SetDataInput(FitPar::GetDataBase() +
+ "/T2K/CC0pi/STV/datResults.root;Result");
+ fSettings.SetCovarInput(FitPar::GetDataBase() +
+ "/T2K/CC0pi/STV/datResults.root;Correlation_Matrix");
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = GetEventHistogram()->Integral("width") * double(1E-38) /
- (double(fNEvents) * TotalIntegratedFlux("width"));
+ (double(fNEvents) * TotalIntegratedFlux("width"));
// Plot Setup -------------------------------------------------------
- //SetDataFromTextFile( fSettings.GetDataInput() );
- //ScaleData(1E-38);
- //SetCovarFromDiagonal();
-
- SetDataFromRootFile( fSettings.GetDataInput() );
- SetCorrelationFromRootFile(fSettings.GetCovarInput() );
+ // SetDataFromTextFile( fSettings.GetDataInput() );
+ // ScaleData(1E-38);
+ // SetCovarFromDiagonal();
+ SetDataFromRootFile(fSettings.GetDataInput());
+ SetCorrelationFromRootFile(fSettings.GetCovarInput());
// Final setup ---------------------------------------------------
FinaliseMeasurement();
-
};
void T2K_CC0pinp_STV_XSec_1Ddat_nu::FillEventVariables(FitEvent *event) {
fXVar = FitUtils::Get_STV_dalphat(event, 14, true);
return;
};
//********************************************************************
bool T2K_CC0pinp_STV_XSec_1Ddat_nu::isSignal(FitEvent *event)
//********************************************************************
{
return SignalDef::isT2K_CC0pi_STV(event, EnuMin, EnuMax);
}
diff --git a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddphit_nu.cxx b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddphit_nu.cxx
index e14e2b6..cdc60dd 100644
--- a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddphit_nu.cxx
+++ b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddphit_nu.cxx
@@ -1,88 +1,90 @@
// 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 <http://www.gnu.org/licenses/>.
-*******************************************************************************/
+ * 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 <http://www.gnu.org/licenses/>.
+ *******************************************************************************/
#include "T2K_SignalDef.h"
#include "T2K_CC0pinp_STV_XSec_1Ddphit_nu.h"
-
-//********************************************************************
-T2K_CC0pinp_STV_XSec_1Ddphit_nu::T2K_CC0pinp_STV_XSec_1Ddphit_nu(nuiskey samplekey) {
//********************************************************************
+T2K_CC0pinp_STV_XSec_1Ddphit_nu::T2K_CC0pinp_STV_XSec_1Ddphit_nu(
+ nuiskey samplekey) {
+ //********************************************************************
// Sample overview ---------------------------------------------------
- std::string descrip = "T2K_CC0pinp_STV_XSec_1Ddphit_nu sample. \n" \
- "Target: CH \n" \
- "Flux: T2K 2.5 degree off-axis (ND280) \n" \
- "Signal: CC0piNp (N>=1) with 450M eV < p_p < 1 GeV \n"\
- " p_mu > 250 MeV \n"\
- " cth_p > 0.6 \n"\
- " cth_mu > -0.4 \n";\
+ std::string descrip = "T2K_CC0pinp_STV_XSec_1Ddphit_nu sample. \n"
+ "Target: CH \n"
+ "Flux: T2K 2.5 degree off-axis (ND280) \n"
+ "Signal: CC0piNp (N>=1) with 450M eV < p_p < 1 GeV \n"
+ " p_mu > 250 MeV \n"
+ " cth_p > 0.6 \n"
+ " cth_mu > -0.4 \n"
+ "https://doi.org/10.1103/PhysRevD.98.032003 \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("#delta#it{#phi}_{T} (GeV c^{-1})");
- fSettings.SetYTitle("#frac{d#sigma}{d#delta#it{#phi}_{T}} (cm^{2} nucleon^{-1} rads^{-1})");
+ fSettings.SetYTitle(
+ "#frac{d#sigma}{d#delta#it{#phi}_{T}} (cm^{2} nucleon^{-1} rads^{-1})");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG");
fSettings.SetEnuRange(0.0, 50.0);
fSettings.DefineAllowedTargets("C,H");
// CCQELike plot information
fSettings.SetTitle("T2K_CC0pinp_STV_XSec_1Ddphit_nu");
- //fSettings.SetDataInput( GeneralUtils::GetTopLevelDir() + "/data/T2K/T2K_CC0pinp_STV_XSec_1Ddphit_nu.dat");
-
- fSettings.SetDataInput( FitPar::GetDataBase() + "/T2K/CC0pi/STV/dphitResults.root;Result" );
- fSettings.SetCovarInput( FitPar::GetDataBase() + "/T2K/CC0pi/STV/dphitResults.root;Correlation_Matrix" );
+ // fSettings.SetDataInput( GeneralUtils::GetTopLevelDir() +
+ // "/data/T2K/T2K_CC0pinp_STV_XSec_1Ddphit_nu.dat");
+
+ fSettings.SetDataInput(FitPar::GetDataBase() +
+ "/T2K/CC0pi/STV/dphitResults.root;Result");
+ fSettings.SetCovarInput(
+ FitPar::GetDataBase() +
+ "/T2K/CC0pi/STV/dphitResults.root;Correlation_Matrix");
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = GetEventHistogram()->Integral("width") * double(1E-38) /
- (double(fNEvents) * TotalIntegratedFlux("width"));
+ (double(fNEvents) * TotalIntegratedFlux("width"));
// Plot Setup -------------------------------------------------------
- //SetDataFromTextFile( fSettings.GetDataInput() );
- //ScaleData(1E-38);
- //SetCovarFromDiagonal();
+ // SetDataFromTextFile( fSettings.GetDataInput() );
+ // ScaleData(1E-38);
+ // SetCovarFromDiagonal();
- SetDataFromRootFile( fSettings.GetDataInput() );
- SetCorrelationFromRootFile(fSettings.GetCovarInput() );
+ SetDataFromRootFile(fSettings.GetDataInput());
+ SetCorrelationFromRootFile(fSettings.GetCovarInput());
// Final setup ---------------------------------------------------
FinaliseMeasurement();
-
};
-
-
-
void T2K_CC0pinp_STV_XSec_1Ddphit_nu::FillEventVariables(FitEvent *event) {
fXVar = FitUtils::Get_STV_dphit(event, 14, true);
return;
};
//********************************************************************
bool T2K_CC0pinp_STV_XSec_1Ddphit_nu::isSignal(FitEvent *event)
//********************************************************************
{
return SignalDef::isT2K_CC0pi_STV(event, EnuMin, EnuMax);
}
diff --git a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.cxx b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.cxx
index 0b862ed..809db8a 100644
--- a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.cxx
+++ b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.cxx
@@ -1,89 +1,90 @@
// 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 <http://www.gnu.org/licenses/>.
-*******************************************************************************/
+ * 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 <http://www.gnu.org/licenses/>.
+ *******************************************************************************/
#include "T2K_SignalDef.h"
#include "T2K_CC0pinp_STV_XSec_1Ddpt_nu.h"
-
-//********************************************************************
-T2K_CC0pinp_STV_XSec_1Ddpt_nu::T2K_CC0pinp_STV_XSec_1Ddpt_nu(nuiskey samplekey) {
//********************************************************************
+T2K_CC0pinp_STV_XSec_1Ddpt_nu::T2K_CC0pinp_STV_XSec_1Ddpt_nu(
+ nuiskey samplekey) {
+ //********************************************************************
// Sample overview ---------------------------------------------------
- std::string descrip = "T2K_CC0pinp_STV_XSec_1Ddpt_nu sample. \n" \
- "Target: CH \n" \
- "Flux: T2K 2.5 degree off-axis (ND280) \n" \
- "Signal: CC0piNp (N>=1) with 450M eV < p_p < 1 GeV \n"\
- " p_mu > 250 MeV \n"\
- " cth_p > 0.6 \n"\
- " cth_mu > -0.4 \n";\
+ std::string descrip = "T2K_CC0pinp_STV_XSec_1Ddpt_nu sample. \n"
+ "Target: CH \n"
+ "Flux: T2K 2.5 degree off-axis (ND280) \n"
+ "Signal: CC0piNp (N>=1) with 450M eV < p_p < 1 GeV \n"
+ " p_mu > 250 MeV \n"
+ " cth_p > 0.6 \n"
+ " cth_mu > -0.4 \n"
+ "https://doi.org/10.1103/PhysRevD.98.032003 \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("#delta#it{p}_{T} (GeV c^{-1})");
- fSettings.SetYTitle("#frac{d#sigma}{d#delta#it{p}_{T}} (cm^{2} nucleon^{-1} GeV^{-1} c)");
+ fSettings.SetYTitle(
+ "#frac{d#sigma}{d#delta#it{p}_{T}} (cm^{2} nucleon^{-1} GeV^{-1} c)");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG");
fSettings.SetEnuRange(0.0, 50.0);
fSettings.DefineAllowedTargets("C,H");
// CCQELike plot information
fSettings.SetTitle("T2K_CC0pinp_STV_XSec_1Ddpt_nu");
- //fSettings.SetDataInput( GeneralUtils::GetTopLevelDir() + "/data/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.dat");
-
- fSettings.SetDataInput( FitPar::GetDataBase() + "/T2K/CC0pi/STV/dptResults.root;Result" );
- fSettings.SetCovarInput( FitPar::GetDataBase() + "/T2K/CC0pi/STV/dptResults.root;Correlation_Matrix" );
+ // fSettings.SetDataInput( GeneralUtils::GetTopLevelDir() +
+ // "/data/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.dat");
+ fSettings.SetDataInput(FitPar::GetDataBase() +
+ "/T2K/CC0pi/STV/dptResults.root;Result");
+ fSettings.SetCovarInput(FitPar::GetDataBase() +
+ "/T2K/CC0pi/STV/dptResults.root;Correlation_Matrix");
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = GetEventHistogram()->Integral("width") * double(1E-38) /
- (double(fNEvents) * TotalIntegratedFlux("width"));
+ (double(fNEvents) * TotalIntegratedFlux("width"));
// Plot Setup -------------------------------------------------------
- //SetDataFromTextFile( fSettings.GetDataInput() );
- //SetCovarFromDiagonal();
- //ScaleData(1E-38);
-
-
- SetDataFromRootFile( fSettings.GetDataInput() );
- SetCorrelationFromRootFile(fSettings.GetCovarInput() );
- //SetCovarianceFromRootFile(fSettings.GetCovarInput() );
+ // SetDataFromTextFile( fSettings.GetDataInput() );
+ // SetCovarFromDiagonal();
+ // ScaleData(1E-38);
+ SetDataFromRootFile(fSettings.GetDataInput());
+ SetCorrelationFromRootFile(fSettings.GetCovarInput());
+ // SetCovarianceFromRootFile(fSettings.GetCovarInput() );
// Final setup ---------------------------------------------------
FinaliseMeasurement();
-
};
void T2K_CC0pinp_STV_XSec_1Ddpt_nu::FillEventVariables(FitEvent *event) {
fXVar = FitUtils::Get_STV_dpt(event, 14, true) / 1000.0;
return;
};
//********************************************************************
bool T2K_CC0pinp_STV_XSec_1Ddpt_nu::isSignal(FitEvent *event)
//********************************************************************
{
return SignalDef::isT2K_CC0pi_STV(event, EnuMin, EnuMax);
}