diff --git a/app/CMakeLists.txt b/app/CMakeLists.txt
index 3983c0a..c260ba6 100644
--- a/app/CMakeLists.txt
+++ b/app/CMakeLists.txt
@@ -1,178 +1,192 @@
# Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
################################################################################
# This file is part of NUISANCE.
#
# NUISANCE is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# NUISANCE is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with NUISANCE. If not, see .
################################################################################
set(TARGETS_TO_BUILD)
if(USE_MINIMIZER)
add_executable(nuismin nuismin.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuismin)
target_link_libraries(nuismin ${MODULETargets})
target_link_libraries(nuismin ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(nuismin ${ROOT_LIBS})
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(nuismin PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
add_executable(nuissplines nuissplines.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuissplines)
target_link_libraries(nuissplines ${MODULETargets})
target_link_libraries(nuissplines ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(nuissplines ${ROOT_LIBS})
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(nuissplines PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
endif()
include_directories(${RWENGINE_INCLUDE_DIRECTORIES})
include_directories(${CMAKE_SOURCE_DIR}/src/Routines)
include_directories(${CMAKE_SOURCE_DIR}/src/InputHandler)
include_directories(${CMAKE_SOURCE_DIR}/src/Genie)
include_directories(${CMAKE_SOURCE_DIR}/src/FitBase)
include_directories(${CMAKE_SOURCE_DIR}/src/Utils)
include_directories(${CMAKE_SOURCE_DIR}/src/Splines)
include_directories(${CMAKE_SOURCE_DIR}/src/Reweight)
include_directories(${CMAKE_SOURCE_DIR}/src/FCN)
include_directories(${CMAKE_SOURCE_DIR}/src/MCStudies)
include_directories(${CMAKE_SOURCE_DIR}/src/Smearceptance)
include_directories(${EXP_INCLUDE_DIRECTORIES})
if (USE_NuWro AND NOT NUWRO_BUILT_FROM_FILE)
add_executable(nuwro_nuisance nuwro_NUISANCE.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuwro_nuisance)
target_link_libraries(nuwro_nuisance ${MODULETargets})
target_link_libraries(nuwro_nuisance ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(nuwro_nuisance ${ROOT_LIBS})
include_directories(${CMAKE_SOURCE_DIR}/src/FitBase)
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(nuwro_nuisance PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
endif()
if (USE_NEUT)
add_executable(neut_nuisance neut_NUISANCE.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};neut_nuisance)
target_link_libraries(neut_nuisance ${MODULETargets})
target_link_libraries(neut_nuisance ${CMAKE_DEPENDLIB_FLAGS})
target_link_libraries(neut_nuisance ${ROOT_LIBS})
include_directories(${CMAKE_SOURCE_DIR}/src/FitBase)
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(neut_nuisance PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
endif()
if (BUILD_GEVGEN)
add_executable(gevgen_nuisance gEvGen_NUISANCE.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};gevgen_nuisance)
target_link_libraries(gevgen_nuisance ${MODULETargets})
target_link_libraries(gevgen_nuisance ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(gevgen_nuisance ${ROOT_LIBS})
include_directories(${CMAKE_SOURCE_DIR}/src/FitBase)
include_directories(${GENIE_INCLUDES}/Apps)
include_directories(${GENIE_INCLUDES}/FluxDrivers)
include_directories(${GENIE_INCLUDES}/EVGDrivers)
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(gevgen_nuisance PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
+
+
+ add_executable(gevgen_nuisance_mixed gEvGen_NUISANCE_MIXED.cxx)
+ set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};gevgen_nuisance_mixed)
+ target_link_libraries(gevgen_nuisance_mixed ${MODULETargets})
+ target_link_libraries(gevgen_nuisance_mixed ${CMAKE_DEPENDLIB_FLAGS})
+ # target_link_libraries(gevgen_nuisance_mixed ${ROOT_LIBS})
+ include_directories(${CMAKE_SOURCE_DIR}/src/FitBase)
+ include_directories(${GENIE_INCLUDES}/Apps)
+ include_directories(${GENIE_INCLUDES}/FluxDrivers)
+ include_directories(${GENIE_INCLUDES}/EVGDrivers)
+ if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
+ set_target_properties(gevgen_nuisance_mixed PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
+ endif()
endif()
if (USE_GiBUU)
add_executable(DumpGiBUUEvents DumpGiBUUEvents.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};DumpGiBUUEvents)
target_link_libraries(DumpGiBUUEvents ${MODULETargets})
target_link_libraries(DumpGiBUUEvents ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(DumpGiBUUEvents ${ROOT_LIBS})
include_directories(${CMAKE_SOURCE_DIR}/src/FitBase)
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(DumpGiBUUEvents PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
endif()
add_executable(nuiscomp nuiscomp.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuiscomp)
target_link_libraries(nuiscomp ${MODULETargets})
target_link_libraries(nuiscomp ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(nuiscomp ${ROOT_LIBS})
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(nuiscomp PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
add_executable(nuisflat nuisflat.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuisflat)
target_link_libraries(nuisflat ${MODULETargets})
target_link_libraries(nuisflat ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(nuisflat ${ROOT_LIBS})
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(nuisflat PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
add_executable(nuissmear nuissmear.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuissmear)
target_link_libraries(nuissmear ${MODULETargets})
target_link_libraries(nuissmear ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(nuissmear ${ROOT_LIBS})
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(nuissmear PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
add_executable(nuissyst nuissyst.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuissyst)
target_link_libraries(nuissyst ${MODULETargets})
target_link_libraries(nuissyst ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(nuissyst ${ROOT_LIBS})
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(nuissyst PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
if(USE_GENIE)
add_executable(PrepareGENIE PrepareGENIE.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};PrepareGENIE)
target_link_libraries(PrepareGENIE ${MODULETargets})
target_link_libraries(PrepareGENIE ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(PrepareGENIE ${ROOT_LIBS})
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(PrepareGENIE PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
endif()
if(USE_NEUT)
add_executable(PrepareNEUT PrepareNEUT.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};PrepareNEUT)
target_link_libraries(PrepareNEUT ${MODULETargets})
target_link_libraries(PrepareNEUT ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(PrepareNEUT ${ROOT_LIBS})
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(PrepareNEUT PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
endif()
# PREPARE NUWRO
# Commented out for the time being until it is finished..
if(USE_NuWro)
add_executable(PrepareNuwro PrepareNuwroEvents.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};PrepareNuwro)
target_link_libraries(PrepareNuwro ${MODULETargets})
target_link_libraries(PrepareNuwro ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(PrepareNuwro ${ROOT_LIBS})
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(PrepareNuwro PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
endif()
install(TARGETS ${TARGETS_TO_BUILD} DESTINATION bin)
diff --git a/src/FCN/JointFCN.cxx b/src/FCN/JointFCN.cxx
index a4136e2..31023b2 100755
--- a/src/FCN/JointFCN.cxx
+++ b/src/FCN/JointFCN.cxx
@@ -1,1087 +1,1010 @@
#include "JointFCN.h"
#include
#include "FitUtils.h"
-// //***************************************************
-// JointFCN::JointFCN(std::string cardfile, TFile* outfile) {
-// //***************************************************
-
-// fOutputDir = gDirectory;
-// FitPar::Config().out = outfile;
-
-// fCardFile = cardfile;
-
-// LoadSamples(fCardFile);
-
-// fCurIter = 0;
-// fMCFilled = false;
-
-// fOutputDir->cd();
-
-// fIterationTree = NULL;
-// fDialVals = NULL;
-// fNDials = 0;
-
-// fUsingEventManager = FitPar::Config().GetParB("EventManager");
-// fOutputDir->cd();
-// };
+//***************************************************
JointFCN::JointFCN(TFile* outfile) {
+//***************************************************
+
fOutputDir = gDirectory;
if (outfile) FitPar::Config().out = outfile;
std::vector samplekeys = Config::QueryKeys("sample");
LoadSamples(samplekeys);
std::vector covarkeys = Config::QueryKeys("covar");
LoadPulls(covarkeys);
fCurIter = 0;
fMCFilled = false;
- fIterationTree = NULL;
+ fIterationTree = false;
fDialVals = NULL;
fNDials = 0;
fUsingEventManager = FitPar::Config().GetParB("EventManager");
fOutputDir->cd();
}
+//***************************************************
JointFCN::JointFCN(std::vector samplekeys, TFile* outfile) {
+//***************************************************
+
fOutputDir = gDirectory;
if (outfile) FitPar::Config().out = outfile;
LoadSamples(samplekeys);
fCurIter = 0;
fMCFilled = false;
fOutputDir->cd();
- fIterationTree = NULL;
+ fIterationTree = false;
fDialVals = NULL;
fNDials = 0;
fUsingEventManager = FitPar::Config().GetParB("EventManager");
fOutputDir->cd();
}
//***************************************************
JointFCN::~JointFCN() {
//***************************************************
// Delete Samples
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
delete exp;
}
for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) {
ParamPull* pull = *iter;
delete pull;
}
// Sort Tree
if (fIterationTree) DestroyIterationTree();
if (fDialVals) delete fDialVals;
if (fSampleLikes) delete fSampleLikes;
};
//***************************************************
void JointFCN::CreateIterationTree(std::string name, FitWeight* rw) {
- //***************************************************
-
- LOG(FIT) << " Creating new iteration tree! " << std::endl;
- if (fIterationTree && !name.compare(fIterationTree->GetName())) {
- fIterationTree->Reset();
- return;
- }
-
- fIterationTree = new TTree(name.c_str(), name.c_str());
-
- // Setup Main Branches
- fIterationTree->Branch("total_likelihood", &fLikelihood,
- "total_likelihood/D");
- fIterationTree->Branch("total_ndof", &fNDOF, "total_ndof/I");
+//***************************************************
- // Setup Sample Arrays
- int ninputs = fSamples.size() + fPulls.size();
- fSampleLikes = new double[ninputs];
- fSampleNDOF = new int[ninputs];
- fNDOF = GetNDOF();
+ LOG(FIT) << " Creating new iteration container! " << std::endl;
+ DestroyIterationTree();
+ fIterationTreeName = name;
- // Setup Sample Branches
- int count = 0;
- for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
+ // Add sample likelihoods and ndof
+ for (MeasListConstIter iter = fSamples.begin();
+ iter != fSamples.end();
iter++) {
- MeasurementBase* exp = *iter;
+ MeasurementBase* exp = *iter;
std::string name = exp->GetName();
- std::string liketag = name + "_likelihood";
- std::string ndoftag = name + "_ndof";
- fIterationTree->Branch(liketag.c_str(), &fSampleLikes[count],
- (liketag + "/D").c_str());
- fIterationTree->Branch(ndoftag.c_str(), &fSampleNDOF[count],
- (ndoftag + "/D").c_str());
+ std::string liketag = name + "_likelihood";
+ fNameValues.push_back(liketag);
+ fCurrentValues.push_back(0.0);
- count++;
+ std::string ndoftag = name + "_ndof";
+ fNameValues.push_back(ndoftag);
+ fCurrentValues.push_back(0.0);
}
- for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) {
- ParamPull* pull = *iter;
+ // Add Pull terms
+ for (PullListConstIter iter = fPulls.begin();
+ iter != fPulls.end(); iter++) {
+ ParamPull* pull = *iter;
std::string name = pull->GetName();
+
std::string liketag = name + "_likelihood";
+ fNameValues.push_back(liketag);
+ fCurrentValues.push_back(0.0);
+
std::string ndoftag = name + "_ndof";
+ fNameValues.push_back(ndoftag);
+ fCurrentValues.push_back(0.0);
+ }
- fIterationTree->Branch(liketag.c_str(), &fSampleLikes[count],
- (liketag + "/D").c_str());
- fIterationTree->Branch(ndoftag.c_str(), &fSampleNDOF[count],
- (ndoftag + "/D").c_str());
+ // Add Likelihoods
+ fNameValues.push_back("total_likelihood");
+ fCurrentValues.push_back(0.0);
- count++;
- }
+ fNameValues.push_back("total_ndof");
+ fCurrentValues.push_back(0.0);
- // Add Dial Branches
- std::vector dials = rw->GetDialNames();
- fNDials = dials.size();
- fDialVals = new double[fNDials];
+ // Setup Containers
+ fSampleN = fSamples.size() + fPulls.size();
+ fSampleLikes = new double[fSampleN];
+ fSampleNDOF = new int[fSampleN];
- for (int i = 0; i < fNDials; i++) {
- fIterationTree->Branch(dials[i].c_str(), &fDialVals[i],
- (dials[i] + "/D").c_str());
+ // Add Dials
+ std::vector dials = rw->GetDialNames();
+ for (size_t i = 0; i < dials.size(); i++){
+ fNameValues.push_back( dials[i] );
+ fCurrentValues.push_back( 0.0 );
}
+ fNDials = dials.size();
+ fDialVals = new double[fNDials];
+
+ // Set IterationTree Flag
+ fIterationTree = true;
+
}
//***************************************************
void JointFCN::DestroyIterationTree() {
- //***************************************************
+//***************************************************
+
+ fIterationCount.clear();
+ fCurrentValues.clear();
+ fNameValues.clear();
+ fIterationValues.clear();
- if (!fIterationTree) {
- delete fIterationTree;
- }
}
//***************************************************
void JointFCN::WriteIterationTree() {
- //***************************************************
+//***************************************************
+ LOG(FIT) << "Writing iteration tree" << std::endl;
+
+ // Make a new TTree
+ TTree* itree = new TTree(fIterationTreeName.c_str(),
+ fIterationTreeName.c_str());
+
+ double* vals = new double[fNameValues.size()];
+ int count = 0;
+
+ itree->Branch("iteration",&count,"Iteration/I");
+ for (int i = 0; i < fNameValues.size(); i++) {
+ itree->Branch( fNameValues[i].c_str(),
+ &vals[i],
+ (fNameValues[i] + "/D").c_str() );
+ }
+
+ // Fill Iterations
+ for (size_t i = 0; i < fIterationValues.size(); i++){
+ std::vector itervals = fIterationValues[i];
+
+ // Fill iteration state
+ count = fIterationCount[i];
+ for (size_t j = 0; j < itervals.size(); j++){
+ vals[j] = itervals[j];
+ }
- if (!fIterationTree) {
- ERR(FTL) << "Can't save empty iteration tree!" << std::endl;
- throw;
+ // Save to TTree
+ itree->Fill();
}
- fIterationTree->Write();
+
+ // Write to file
+ itree->Write();
}
//***************************************************
void JointFCN::FillIterationTree(FitWeight* rw) {
- //***************************************************
+//***************************************************
- if (!fIterationTree) {
- ERR(FTL) << "Trying to fill iteration_tree when it is NULL!" << std::endl;
- throw;
+ // Loop over samples count
+ int count = 0;
+ for (int i = 0; i < fSampleN; i++){
+ fCurrentValues[count++] = fSampleLikes[i];
+ fCurrentValues[count++] = double(fSampleNDOF[i]);
}
+ // Fill Totals
+ fCurrentValues[count++] = fLikelihood;
+ fCurrentValues[count++] = double(fNDOF);
+
+ // Loop Over Parameter Counts
rw->GetAllDials(fDialVals, fNDials);
- fIterationTree->Fill();
+ for (int i = 0; i < fNDials; i++){
+ fCurrentValues[count++] = double(fDialVals[i]);
+ }
+
+ // Push Back Into Container
+ fIterationCount.push_back( fCurIter );
+ fIterationValues.push_back(fCurrentValues);
}
//***************************************************
double JointFCN::DoEval(const double* x) {
//***************************************************
// WEIGHT ENGINE
fDialChanged = FitBase::GetRW()->HasRWDialChanged(x);
FitBase::GetRW()->UpdateWeightEngine(x);
if (fDialChanged) {
FitBase::GetRW()->Reconfigure();
FitBase::EvtManager().ResetWeightFlags();
}
if (LOG_LEVEL(REC)) {
FitBase::GetRW()->Print();
}
// SORT SAMPLES
ReconfigureSamples();
// GET TEST STAT
fLikelihood = GetLikelihood();
-
+ fNDOF = GetNDOF();
+
// PRINT PROGRESS
LOG(FIT) << "Current Stat (iter. " << this->fCurIter << ") = " << fLikelihood
<< std::endl;
// UPDATE TREE
if (fIterationTree) FillIterationTree(FitBase::GetRW());
return fLikelihood;
}
//***************************************************
int JointFCN::GetNDOF() {
//***************************************************
int totaldof = 0;
int count = 0;
// Total number of Free bins in each MC prediction
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
int dof = exp->GetNDOF();
// Save Seperate DOF
if (fIterationTree) {
fSampleNDOF[count] = dof;
}
// Add to total
totaldof += dof;
count++;
}
// Loop over pulls
for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) {
ParamPull* pull = *iter;
double dof = pull->GetLikelihood();
// Save seperate DOF
if (fIterationTree) {
fSampleNDOF[count] = dof;
}
// Add to total
totaldof += dof;
count++;
}
// Set Data Variable
- fNDOF = totaldof;
+ fSampleNDOF[count] = totaldof;
return totaldof;
}
//***************************************************
double JointFCN::GetLikelihood() {
//***************************************************
LOG(MIN) << std::left << std::setw(43) << "Getting likelihoods..."
<< " : "
<< "-2logL" << std::endl;
// Loop and add up likelihoods in an uncorrelated way
double like = 0.0;
int count = 0;
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
double newlike = exp->GetLikelihood();
int ndof = exp->GetNDOF();
// Save seperate likelihoods
if (fIterationTree) {
fSampleLikes[count] = newlike;
}
LOG(MIN) << "-> " << std::left << std::setw(40) << exp->GetName() << " : "
<< newlike << "/" << ndof << std::endl;
// Add Weight Scaling
// like *= FitBase::GetRW()->GetSampleLikelihoodWeight(exp->GetName());
// Add to total
like += newlike;
count++;
}
// Loop over pulls
for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) {
ParamPull* pull = *iter;
double newlike = pull->GetLikelihood();
// Save seperate likelihoods
if (fIterationTree) {
fSampleLikes[count] = newlike;
}
// Add to total
like += newlike;
count++;
}
// Set Data Variable
fLikelihood = like;
return like;
};
void JointFCN::LoadSamples(std::vector samplekeys) {
LOG(MIN) << "Loading Samples : " << samplekeys.size() << std::endl;
for (size_t i = 0; i < samplekeys.size(); i++) {
nuiskey key = samplekeys[i];
// Get Sample Options
std::string samplename = key.GetS("name");
std::string samplefile = key.GetS("input");
std::string sampletype = key.GetS("type");
std::string fakeData = "";
LOG(MIN) << "Loading Sample : " << samplename << std::endl;
fOutputDir->cd();
MeasurementBase* NewLoadedSample = SampleUtils::CreateSample(key);
if (!NewLoadedSample) {
ERR(FTL) << "Could not load sample provided: " << samplename << std::endl;
ERR(FTL) << "Check spelling with that in src/FCN/SampleList.cxx"
<< std::endl;
throw;
} else {
fSamples.push_back(NewLoadedSample);
}
}
}
+//***************************************************
void JointFCN::LoadPulls(std::vector pullkeys) {
+//***************************************************
for (size_t i = 0; i < pullkeys.size(); i++) {
nuiskey key = pullkeys[i];
std::string pullname = key.GetS("name");
std::string pullfile = key.GetS("input");
std::string pulltype = key.GetS("type");
fOutputDir->cd();
- std::cout << "Creating Pull Term : " << std::endl;
- sleep(1);
fPulls.push_back(new ParamPull(pullname, pullfile, pulltype));
}
-
- // // Sample Inputs
- // if (!samplevect[0].compare("covar") || !samplevect[0].compare("pulls")
- // ||
- // !samplevect[0].compare("throws")) {
- // // Get all inputs
- // std::string name = samplevect[1];
- // std::string files = samplevect[2];
- // std::string type = "DEFAULT";
-
- // if (samplevect.size() > 3) {
- // type = samplevect[3];
- // } else if (!samplevect[0].compare("pull")) {
- // type = "GAUSPULL";
- // } else if (!samplevect[0].compare("throws")) {
- // type = "GAUSTHROWS";
- // }
-
- // // Create Pull Class
- // LOG(MIN) << "Loading up pull term: " << name << " << " << files << "
- // ("
- // << type << ")" << std::endl;
- // std::string fakeData = "";
- // fOutputDir->cd();
- // fPulls.push_back(new ParamPull(name, files, type));
- // }
}
-// //***************************************************
-// void JointFCN::LoadSamples(std::string cardinput)
-// //***************************************************
-// {
-// LOG(MIN) << "Initializing Samples" << std::endl;
-
-// // Read the card file here and load objects
-// std::string line;
-// std::ifstream card(cardinput.c_str(), ifstream::in);
-
-// // Make sure they are created in correct working DIR
-// fOutputDir->cd();
-
-// while (std::getline(card >> std::ws, line, '\n')) {
-// // Skip Empties
-// if (line.c_str()[0] == '#') continue;
-// if (line.empty()) continue;
-
-// // Parse line
-// std::vector samplevect = GeneralUtils::ParseToStr(line, "
-// ");
-
-// // Sample Inputs
-// if (!samplevect[0].compare("sample")) {
-// // Get all inputs
-// std::string name = samplevect[1];
-// std::string files = samplevect[2];
-// std::string type = "DEFAULT";
-// if (samplevect.size() > 3) type = samplevect[3];
-
-// // Create Sample Class
-// LOG(MIN) << "Loading up sample: " << name << " << " << files << " ("
-// << type << ")" << std::endl;
-// std::string fakeData = "";
-// fOutputDir->cd();
-// MeasurementBase* NewLoadedSample = SampleUtils::CreateSample(name,
-// files, type,
-// fakeData, FitBase::GetRW());
-
-// if (!NewLoadedSample) {
-// ERR(FTL) << "Could not load sample provided: " << name << std::endl;
-// ERR(FTL) << "Check spelling with that in src/FCN/SampleList.cxx"
-// << std::endl;
-// throw;
-// } else {
-// fSamples.push_back(NewLoadedSample);
-// }
-// }
-
-// // Sample Inputs
-// if (!samplevect[0].compare("covar") || !samplevect[0].compare("pulls") ||
-// !samplevect[0].compare("throws")) {
-// // Get all inputs
-// std::string name = samplevect[1];
-// std::string files = samplevect[2];
-// std::string type = "DEFAULT";
-
-// if (samplevect.size() > 3) {
-// type = samplevect[3];
-// } else if (!samplevect[0].compare("pull")) {
-// type = "GAUSPULL";
-// } else if (!samplevect[0].compare("throws")) {
-// type = "GAUSTHROWS";
-// }
-
-// // Create Pull Class
-// LOG(MIN) << "Loading up pull term: " << name << " << " << files << " ("
-// << type << ")" << std::endl;
-// std::string fakeData = "";
-// fOutputDir->cd();
-// fPulls.push_back(new ParamPull(name, files, type));
-// }
-// }
-// card.close();
-// };
-
//***************************************************
void JointFCN::ReconfigureSamples(bool fullconfig) {
- //***************************************************
+//***************************************************
int starttime = time(NULL);
LOG(REC) << "Starting Reconfigure iter. " << this->fCurIter << std::endl;
// std::cout << fUsingEventManager << " " << fullconfig << " " << fMCFilled <<
// std::endl;
// Event Manager Reconf
if (fUsingEventManager) {
if (!fullconfig and fMCFilled)
ReconfigureFastUsingManager();
else
ReconfigureUsingManager();
} else {
// Loop over all Measurement Classes
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
// If RW Either do signal or full reconfigure.
if (fDialChanged or !fMCFilled or fullconfig) {
if (!fullconfig and fMCFilled)
exp->ReconfigureFast();
else
exp->Reconfigure();
// If RW Not needed just do normalisation
} else {
exp->Renormalise();
}
}
}
// Loop over pulls and update
for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) {
ParamPull* pull = *iter;
pull->Reconfigure();
}
fMCFilled = true;
LOG(MIN) << "Finished Reconfigure iter. " << fCurIter << " in "
<< time(NULL) - starttime << "s" << std::endl;
fCurIter++;
}
//***************************************************
void JointFCN::ReconfigureSignal() {
- //***************************************************
- this->ReconfigureSamples(false);
+//***************************************************
+ ReconfigureSamples(false);
}
//***************************************************
void JointFCN::ReconfigureAllEvents() {
//***************************************************
FitBase::GetRW()->Reconfigure();
FitBase::EvtManager().ResetWeightFlags();
ReconfigureSamples(true);
}
std::vector JointFCN::GetInputList() {
std::vector InputList;
fIsAllSplines = true;
MeasListConstIter iterSam = fSamples.begin();
for (; iterSam != fSamples.end(); iterSam++) {
MeasurementBase* exp = (*iterSam);
std::vector subsamples = exp->GetSubSamples();
for (size_t i = 0; i < subsamples.size(); i++) {
InputHandlerBase* inp = subsamples[i]->GetInput();
if (std::find(InputList.begin(), InputList.end(), inp) ==
InputList.end()) {
if (subsamples[i]->GetInput()->GetType() != kSPLINEPARAMETER)
fIsAllSplines = false;
InputList.push_back(subsamples[i]->GetInput());
}
}
}
return InputList;
}
std::vector JointFCN::GetSubSampleList() {
std::vector SampleList;
MeasListConstIter iterSam = fSamples.begin();
for (; iterSam != fSamples.end(); iterSam++) {
MeasurementBase* exp = (*iterSam);
std::vector subsamples = exp->GetSubSamples();
for (size_t i = 0; i < subsamples.size(); i++) {
SampleList.push_back(subsamples[i]);
}
}
return SampleList;
}
//***************************************************
void JointFCN::ReconfigureUsingManager() {
- //***************************************************
+//***************************************************
// 'Slow' Event Manager Reconfigure
LOG(REC) << "Event Manager Reconfigure" << std::endl;
int timestart = time(NULL);
// Reset all samples
MeasListConstIter iterSam = fSamples.begin();
for (; iterSam != fSamples.end(); iterSam++) {
MeasurementBase* exp = (*iterSam);
exp->ResetAll();
}
// If we are siving signal, reset all containers.
bool savesignal = (FitPar::Config().GetParB("SignalReconfigures"));
if (savesignal) {
// Reset all of our event signal vectors
fSignalEventBoxes.clear();
fSignalEventFlags.clear();
fSampleSignalFlags.clear();
fSignalEventSplines.clear();
}
// Make sure we have a list of inputs
if (fInputList.empty()) {
fInputList = GetInputList();
fSubSampleList = GetSubSampleList();
}
// If all inputs are splines make sure the readers are told
// they need to be reconfigured.
std::vector::iterator inp_iter = fInputList.begin();
if (fIsAllSplines) {
for (; inp_iter != fInputList.end(); inp_iter++) {
InputHandlerBase* curinput = (*inp_iter);
// Tell reader in each BaseEvent it needs a Reconfigure next weight calc.
BaseFitEvt* curevent = curinput->FirstBaseEvent();
if (curevent->fSplineRead) {
curevent->fSplineRead->SetNeedsReconfigure(true);
}
}
}
// MAIN INPUT LOOP ====================
int fillcount = 0;
int inputcount = 0;
inp_iter = fInputList.begin();
// Loop over each input in manager
for (; inp_iter != fInputList.end(); inp_iter++) {
InputHandlerBase* curinput = (*inp_iter);
// Get event information
FitEvent* curevent = curinput->FirstNuisanceEvent();
curinput->CreateCache();
int i = 0;
int nevents = curinput->GetNEvents();
int countwidth = nevents / 5;
// Start event loop iterating until we get a NULL pointer.
while (curevent) {
// Get Event Weight
curevent->RWWeight = FitBase::GetRW()->CalcWeight(curevent);
curevent->Weight = curevent->RWWeight * curevent->InputWeight;
double rwweight = curevent->Weight;
// std::cout << "RWWeight = " << curevent->RWWeight << " " <<
// curevent->InputWeight << std::endl;
// Logging
// std::cout << CHECKLOG(1) << std::endl;
if (LOGGING(REC)) {
if (i % countwidth == 0) {
QLOG(REC, curinput->GetName()
- << " : Processed " << i << " events. [M, W] = ["
- << curevent->Mode << ", " << rwweight << "]");
+ << " : Processed " << i << " events. [M, W] = ["
+ << curevent->Mode << ", " << rwweight << "]");
}
}
// Setup flag for if signal found in at least one sample
bool foundsignal = false;
// Create a new signal bitset for this event
std::vector signalbitset(fSubSampleList.size());
// Create a new signal box vector for this event
std::vector signalboxes;
// Start measurement iterator
size_t measitercount = 0;
std::vector::iterator meas_iter =
- fSubSampleList.begin();
+ fSubSampleList.begin();
// Loop over all subsamples (sub in JointMeas)
for (; meas_iter != fSubSampleList.end(); meas_iter++) {
MeasurementBase* curmeas = (*meas_iter);
// Compare input pointers, to current input, skip if not.
// Pointer tells us if it matches without doing ID checks.
if (curinput != curmeas->GetInput()) {
if (savesignal) {
// Set bit to 0 as definitely not signal
signalbitset[measitercount] = 0;
}
// Count up what measurement we are on.
measitercount++;
// Skip sample as input not signal.
continue;
}
// Fill events for matching inputs.
MeasurementVariableBox* box = curmeas->FillVariableBox(curevent);
bool signal = curmeas->isSignal(curevent);
curmeas->SetSignal(signal);
curmeas->FillHistograms(curevent->Weight);
// If its Signal tally up fills
if (signal) {
fillcount++;
}
// If we are saving signal/splines fill the bitset
if (savesignal) {
signalbitset[measitercount] = signal;
}
// If signal save a clone of the event box for use later.
if (savesignal and signal) {
foundsignal = true;
signalboxes.push_back(box->CloneSignalBox());
}
// Keep track of Measurement we are on.
measitercount++;
}
// Once we've filled the measurements, if saving signal
// push back if any sample flagged this event as signal
if (savesignal) {
fSignalEventFlags.push_back(foundsignal);
}
// Save the vector of signal boxes for this event
if (savesignal and foundsignal) {
fSignalEventBoxes.push_back(signalboxes);
fSampleSignalFlags.push_back(signalbitset);
}
// If all inputs are splines we can save the spline coefficients
// for fast in memory reconfigures later.
if (fIsAllSplines and savesignal and foundsignal) {
// Make temp vector to push back with
std::vector coeff;
for (size_t l = 0; l < (UInt_t)curevent->fSplineRead->GetNPar(); l++) {
coeff.push_back(curevent->fSplineCoeff[l]);
}
// Push back to signal event splines. Kept in sync with
// fSignalEventBoxes size.
// int splinecount = fSignalEventSplines.size();
fSignalEventSplines.push_back(coeff);
// if (splinecount % 1000 == 0) {
// std::cout << "Pushed Back Coeff " << splinecount << " : ";
// for (size_t l = 0; l < fSignalEventSplines[splinecount].size(); l++)
// {
// std::cout << " " << fSignalEventSplines[splinecount][l];
// }
// std::cout << std::endl;
// }
}
// Clean up vectors once done with this event
signalboxes.clear();
signalbitset.clear();
// Iterate to the next event.
curevent = curinput->NextNuisanceEvent();
i++;
}
curinput->RemoveCache();
// Keep track of what input we are on.
inputcount++;
}
// End of Event Loop ===============================
// Now event loop is finished loop over all Measurements
// Converting Binned events to XSec Distributions
iterSam = fSamples.begin();
for (; iterSam != fSamples.end(); iterSam++) {
MeasurementBase* exp = (*iterSam);
exp->ConvertEventRates();
}
// Print out statements on approximate memory usage for profiling.
LOG(REC) << "Filled " << fillcount << " signal events." << std::endl;
if (savesignal) {
int mem =
- ( // sizeof(fSignalEventBoxes) +
- // fSignalEventBoxes.size() * sizeof(fSignalEventBoxes.at(0)) +
- sizeof(MeasurementVariableBox1D) * fillcount) *
- 1E-6;
+ ( // sizeof(fSignalEventBoxes) +
+ // fSignalEventBoxes.size() * sizeof(fSignalEventBoxes.at(0)) +
+ sizeof(MeasurementVariableBox1D) * fillcount) *
+ 1E-6;
LOG(REC) << " -> Saved " << fillcount
<< " signal boxes for faster access. (~" << mem << " MB)"
<< std::endl;
if (fIsAllSplines and !fSignalEventSplines.empty()) {
int splmem = sizeof(float) * fSignalEventSplines.size() *
fSignalEventSplines[0].size() * 1E-6;
LOG(REC) << " -> Saved " << fillcount << " " << fSignalEventSplines.size()
<< " spline sets into memory. (~" << splmem << " MB)"
<< std::endl;
}
}
LOG(REC) << "Time taken ReconfigureUsingManager() : "
<< time(NULL) - timestart << std::endl;
// Check SignalReconfigures works for all samples
- if (savesignal){
+ if (savesignal) {
double likefull = GetLikelihood();
ReconfigureFastUsingManager();
double likefast = GetLikelihood();
-
+
if (fabs(likefull - likefast) > 0.0001)
- {
- ERROR(FTL,"Fast and Full Likelihoods DIFFER! : " << likefull << " : " << likefast);
- ERROR(FTL,"This means some samples you are using are not setup to use SignalReconfigures=1");
- ERROR(FTL,"Please turn OFF signal reconfigures.");
- throw;
- } else {
+ {
+ ERROR(FTL, "Fast and Full Likelihoods DIFFER! : " << likefull << " : " << likefast);
+ ERROR(FTL, "This means some samples you are using are not setup to use SignalReconfigures=1");
+ ERROR(FTL, "Please turn OFF signal reconfigures.");
+ throw;
+ } else {
LOG(FIT) << "Likelihoods for FULL and FAST match. Will use FAST next time." << std::endl;
}
}
// End of reconfigure
return;
};
//***************************************************
void JointFCN::ReconfigureFastUsingManager() {
- //***************************************************
+//***************************************************
LOG(FIT) << " -> Doing FAST using manager" << std::endl;
// Get Start time for profilling
int timestart = time(NULL);
// Reset all samples
MeasListConstIter iterSam = fSamples.begin();
for (; iterSam != fSamples.end(); iterSam++) {
MeasurementBase* exp = (*iterSam);
exp->ResetAll();
}
// Check for saved variables if not do a full reconfigure.
if (fSignalEventFlags.empty()) {
ERR(WRN) << "Signal Flags Empty! Using normal manager." << std::endl;
ReconfigureUsingManager();
return;
}
bool fFillNuisanceEvent =
- FitPar::Config().GetParB("FullEventOnSignalReconfigure");
+ FitPar::Config().GetParB("FullEventOnSignalReconfigure");
// Setup fast vector iterators.
std::vector::iterator inpsig_iter = fSignalEventFlags.begin();
std::vector >::iterator box_iter =
- fSignalEventBoxes.begin();
+ fSignalEventBoxes.begin();
std::vector >::iterator spline_iter =
- fSignalEventSplines.begin();
+ fSignalEventSplines.begin();
std::vector >::iterator samsig_iter =
- fSampleSignalFlags.begin();
+ fSampleSignalFlags.begin();
int splinecount = 0;
// Setup stuff for logging
int fillcount = 0;
int nevents = fSignalEventFlags.size();
int countwidth = nevents / 20;
// If All Splines tell splines they need a reconfigure.
std::vector::iterator inp_iter = fInputList.begin();
if (fIsAllSplines) {
LOG(REC) << "All Spline Inputs so using fast spline loop." << std::endl;
for (; inp_iter != fInputList.end(); inp_iter++) {
InputHandlerBase* curinput = (*inp_iter);
// Tell each fSplineRead in BaseFitEvent to reconf next weight calc
BaseFitEvt* curevent = curinput->FirstBaseEvent();
if (curevent->fSplineRead)
curevent->fSplineRead->SetNeedsReconfigure(true);
}
}
// Loop over all possible spline inputs
double* coreeventweights = new double[fSignalEventBoxes.size()];
splinecount = 0;
inp_iter = fInputList.begin();
inpsig_iter = fSignalEventFlags.begin();
spline_iter = fSignalEventSplines.begin();
// Loop over all signal flags
// For each valid signal flag add one to splinecount
// Get Splines from that count and add to weight
// Add splinecount
int sigcount = 0;
splinecount = 0;
// #pragma omp parallel for shared(splinecount,sigcount)
for (uint iinput = 0; iinput < fInputList.size(); iinput++) {
InputHandlerBase* curinput = fInputList[iinput];
BaseFitEvt* curevent = curinput->FirstBaseEvent();
for (int i = 0; i < curinput->GetNEvents(); i++) {
double rwweight = 0.0;
if (fSignalEventFlags[sigcount]) {
// Get Event Info
if (!fIsAllSplines) {
if (fFillNuisanceEvent)
curinput->GetNuisanceEvent(i);
else
curevent = curinput->GetBaseEvent(i);
} else {
curevent->fSplineCoeff = &fSignalEventSplines[splinecount][0];
}
curevent->RWWeight = FitBase::GetRW()->CalcWeight(curevent);
curevent->Weight = curevent->RWWeight * curevent->InputWeight;
rwweight = curevent->Weight;
coreeventweights[splinecount] = rwweight;
if (splinecount % countwidth == 0) {
LOG(REC) << "Processed " << splinecount
<< " event weights. W = " << rwweight << std::endl;
}
// #pragma omp atomic
splinecount++;
}
// #pragma omp atomic
sigcount++;
}
}
LOG(SAM) << "Processed event weights." << std::endl;
// #pragma omp barrier
// Reset Iterators
inpsig_iter = fSignalEventFlags.begin();
spline_iter = fSignalEventSplines.begin();
box_iter = fSignalEventBoxes.begin();
samsig_iter = fSampleSignalFlags.begin();
int nsplineweights = splinecount;
splinecount = 0;
// Start of Fast Event Loop ============================
// Start input iterators
// Loop over number of inputs
for (int ispline = 0; ispline < nsplineweights; ispline++) {
double rwweight = coreeventweights[ispline];
// Get iterators for this event
std::vector::iterator subsamsig_iter = (*samsig_iter).begin();
std::vector::iterator subbox_iter =
- (*box_iter).begin();
+ (*box_iter).begin();
// Loop over all sub measurements.
std::vector::iterator meas_iter = fSubSampleList.begin();
for (; meas_iter != fSubSampleList.end(); meas_iter++, subsamsig_iter++) {
MeasurementBase* curmeas = (*meas_iter);
// If event flagged as signal for this sample fill from the box.
if (*subsamsig_iter) {
curmeas->SetSignal(true);
curmeas->FillHistogramsFromBox((*subbox_iter), rwweight);
// Move onto next box if there is one.
subbox_iter++;
fillcount++;
}
}
if (ispline % countwidth == 0) {
LOG(REC) << "Filled " << ispline << " sample weights." << std::endl;
}
// Iterate over the main signal event containers.
samsig_iter++;
box_iter++;
spline_iter++;
splinecount++;
}
// End of Fast Event Loop ===================
LOG(SAM) << "Filled sample distributions." << std::endl;
// Now loop over all Measurements
// Convert Binned events
iterSam = fSamples.begin();
for (; iterSam != fSamples.end(); iterSam++) {
MeasurementBase* exp = (*iterSam);
exp->ConvertEventRates();
}
// Cleanup coreeventweights
if (fIsAllSplines) {
delete coreeventweights;
}
// Print some reconfigure profiling.
LOG(REC) << "Filled " << fillcount << " signal events." << std::endl;
LOG(REC) << "Time taken ReconfigureFastUsingManager() : "
<< time(NULL) - timestart << std::endl;
}
//***************************************************
void JointFCN::Write() {
- //***************************************************
+//***************************************************
// Save a likelihood/ndof plot
LOG(MIN) << "Writing likelihood plot.." << std::endl;
std::vector likes;
std::vector ndofs;
std::vector names;
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
- iter++){
+ iter++) {
MeasurementBase* exp = *iter;
double like = exp->GetLikelihood();
double ndof = exp->GetNDOF();
std::string name = exp->GetName();
likes.push_back(like);
ndofs.push_back(ndof);
names.push_back(name);
}
- TH1D likehist = TH1D("likelihood_hist","likelihood_hist",
- likes.size(), 0.0, double(likes.size()));
- TH1D ndofhist = TH1D("ndof_hist","ndof_hist",
- ndofs.size(), 0.0, double(ndofs.size()));
- TH1D divhist = TH1D("likedivndof_hist","likedivndof_hist",
- likes.size(), 0.0, double(likes.size()));
- for (size_t i = 0; i < likehist.GetNbinsX(); i++){
- likehist.SetBinContent(i+1, likes[i]);
- ndofhist.SetBinContent(i+1, ndofs[i]);
- if (ndofs[i] != 0.0){
- divhist.SetBinContent(i+1, likes[i]/ndofs[i]);
+ TH1D likehist = TH1D("likelihood_hist", "likelihood_hist",
+ likes.size(), 0.0, double(likes.size()));
+ TH1D ndofhist = TH1D("ndof_hist", "ndof_hist",
+ ndofs.size(), 0.0, double(ndofs.size()));
+ TH1D divhist = TH1D("likedivndof_hist", "likedivndof_hist",
+ likes.size(), 0.0, double(likes.size()));
+ for (size_t i = 0; i < likehist.GetNbinsX(); i++) {
+ likehist.SetBinContent(i + 1, likes[i]);
+ ndofhist.SetBinContent(i + 1, ndofs[i]);
+ if (ndofs[i] != 0.0) {
+ divhist.SetBinContent(i + 1, likes[i] / ndofs[i]);
}
- likehist.GetXaxis()->SetBinLabel(i+1, names[i].c_str());
- ndofhist.GetXaxis()->SetBinLabel(i+1, names[i].c_str());
- divhist.GetXaxis()->SetBinLabel(i+1, names[i].c_str());
+ likehist.GetXaxis()->SetBinLabel(i + 1, names[i].c_str());
+ ndofhist.GetXaxis()->SetBinLabel(i + 1, names[i].c_str());
+ divhist.GetXaxis()->SetBinLabel(i + 1, names[i].c_str());
}
likehist.Write();
ndofhist.Write();
divhist.Write();
// Loop over individual experiments and call Write
LOG(MIN) << "Writing each of the data classes..." << std::endl;
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
exp->Write();
}
// Save Pull Terms
for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) {
ParamPull* pull = *iter;
pull->Write();
}
if (FitPar::Config().GetParB("EventManager")) {
// Get list of inputs
std::map fInputs =
- FitBase::EvtManager().GetInputs();
+ FitBase::EvtManager().GetInputs();
std::map::const_iterator iterInp;
for (iterInp = fInputs.begin(); iterInp != fInputs.end(); iterInp++) {
InputHandlerBase* input = (iterInp->second);
input->GetFluxHistogram()->Write();
input->GetXSecHistogram()->Write();
input->GetEventHistogram()->Write();
}
}
};
//***************************************************
void JointFCN::SetFakeData(std::string fakeinput) {
- //***************************************************
+//***************************************************
LOG(MIN) << "Setting fake data from " << fakeinput << std::endl;
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
exp->SetFakeDataValues(fakeinput);
}
return;
}
//***************************************************
void JointFCN::ThrowDataToy() {
- //***************************************************
+//***************************************************
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
exp->ThrowDataToy();
}
return;
}
diff --git a/src/FCN/JointFCN.h b/src/FCN/JointFCN.h
index 361a75a..bf87424 100755
--- a/src/FCN/JointFCN.h
+++ b/src/FCN/JointFCN.h
@@ -1,165 +1,176 @@
#ifndef _JOINT_FCN_H_
#define _JOINT_FCN_H_
/*!
* \addtogroup FCN
* @{
*/
#include
#include
#include
#include
// ROOT headers
#include "TTree.h"
#include "TH1D.h"
#include "TH2D.h"
#include
#include "TGraphErrors.h"
#include
#include
#include "FitUtils.h"
// External fitter headers
#include "MeasurementBase.h"
#include "SampleList.h"
#define GetCurrentDir getcwd
#include "EventManager.h"
#include "Math/Functor.h"
#include "ParamPull.h"
#include "NuisConfig.h"
#include "NuisKey.h"
#include "MeasurementVariableBox.h"
#include "MeasurementVariableBox1D.h"
using namespace FitUtils;
using namespace FitBase;
//! Main FCN Class which ROOT's joint function needs to evaulate the chi2 at each stage of the fit.
class JointFCN
{
public:
//! Constructor
//! cardfile = Path to input card file listing samples
JointFCN(std::vector samplekeys, TFile* outfile=NULL);
JointFCN(TFile* outfile=NULL); // Loads from global config
//! Destructor
~JointFCN();
//! Create sample list from cardfile
void LoadSamples(std::vector samplekeys);
void LoadPulls(std::vector pullkeys);
//! Main Likelihood evaluation FCN
double DoEval(const double *x);
//! Func Wrapper for ROOT
inline double operator() (const std::vector & x) {
double* x_array = new double[x.size()];
return this->DoEval(x_array);
};
//! Func Wrapper for ROOT
inline double operator() (const double *x) {
return this->DoEval(x);
};
//! Create a TTree to save all dial value iterations for this FCN
void CreateIterationTree(std::string name, FitWeight* rw);
//! Fills dial values and sample likelihoods into tree
void FillIterationTree(FitWeight* rw);
//! Writes TTree to fOutput directory
void WriteIterationTree();
//! Deletes TTree
void DestroyIterationTree();
//! Get Degrees of Freedom for samples (NBins)
int GetNDOF();
//! Return NDOF wrapper
inline unsigned int NDim() {return this->GetNDOF();};
//! Reconfigure samples where we force all events to be looped over.
void ReconfigureAllEvents() ;
//! Call Reconfigure on samples.
//! Choice of reconfigure type depends on whether dials have changed
//! and the MC has been filled.
void ReconfigureSamples(bool fullconfig = false);
//! Call reconfigure on only signal events (defaults to all events if CurIter=0)
void ReconfigureSignal();
//! Gets likelihood for all samples in FCN (assuming uncorrelated)
double GetLikelihood();
//! Returns list of pointers to the samples
inline std::list GetSampleList(){ return fSamples; }
//! Return list of pointers to all the pulls
inline std::list GetPullList(){ return fPulls; };
//! Write all samples to output DIR
void Write();
//! Set Fake data from file/MC
void SetFakeData(std::string fakeinput);
//! Reconfigure looping over duplicate inputs
void ReconfigureUsingManager();
//! Reconfigure Fast looping over duplicate inputs
void ReconfigureFastUsingManager();
/// Throws data according to current stats
void ThrowDataToy();
std::vector GetSubSampleList();
std::vector GetInputList();
private:
//! Append the experiments to include in the fit to this list
std::list fSamples;
//! Append parameter pull terms to include penalties in the fit to this list
std::list fPulls;
TDirectory *fOutputDir; //!< Directory to save contents
std::string fCardFile; //!< Input Card text file
bool fDialChanged; //!< Flag for if RW dials changed
UInt_t fCurIter; //!< Counter for how many times reconfigure called
bool fMCFilled; //!< Check MC has at least been filled once
- TTree* fIterationTree; //!< Tree to save RW values on each function call
+ bool fIterationTree; //!< Tree to save RW values on each function call
int fNDials; //!< Number of RW Dials in FitWeight
double* fDialVals; //!< Current Values of RW Dials
double fLikelihood; //!< Current likelihood for joint sample likelihood
double fNDOF; //!< Total NDOF
double* fSampleLikes; //!< Likelihoods for each individual measurement in list
int * fSampleNDOF; //!< NDOF for each individual measurement in list
bool fUsingEventManager; //!< Flag for doing joint comparisons
std::vector< std::vector > fSignalEventSplines;
std::vector< std::vector > fSignalEventBoxes;
std::vector< bool > fSignalEventFlags;
std::vector< std::vector > fSampleSignalFlags;
std::vector fInputList;
std::vector fSubSampleList;
bool fIsAllSplines;
+
+ std::vector< int > fIterationCount;
+ std::vector< double > fCurrentValues;
+ std::vector< std::string > fNameValues;
+ std::vector< std::vector > fIterationValues;
+ int fSampleN;
+ std::string fIterationTreeName;
+
+
+
+
};
/*! @} */
#endif // _JOINT_FCN_H_
diff --git a/src/FitBase/ParamPull.cxx b/src/FitBase/ParamPull.cxx
index 4a0c6bc..8f9cb6e 100644
--- a/src/FitBase/ParamPull.cxx
+++ b/src/FitBase/ParamPull.cxx
@@ -1,810 +1,844 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#include "ParamPull.h"
//*******************************************************************************
ParamPull::ParamPull(std::string name, std::string inputfile, std::string type, std::string dials) {
//*******************************************************************************
fMinHist = NULL;
fMaxHist = NULL;
fTypeHist = NULL;
fDialSelection = dials;
+ fLimitHist = NULL;
fName = name;
fInput = inputfile;
fType = type;
// Set the pull type
SetType(fType);
std::cout << fType << std::endl;
// Setup Histograms from input file
SetupHistograms(fInput);
};
//*******************************************************************************
void ParamPull::SetType(std::string type) {
//*******************************************************************************
fType = type;
// Assume Default if empty
if (type.empty() || type == "DEFAULT") {
ERR(WRN) << "No type specified for ParmPull class " << fName << std::endl;
ERR(WRN) << "Assuming GAUSTHROW/GAUSPULL" << std::endl;
type = "GAUSTHROW/GAUSPULL";
}
// Set Dial options
if (type.find("FRAC") != std::string::npos) {
fDialOptions = "FRAC";
fPlotTitles = ";; Fractional RW Value";
} else if (type.find("ABS") != std::string::npos) {
fDialOptions = "ABS";
fPlotTitles = ";; ABS RW Value";
} else {
fDialOptions = "";
fPlotTitles = ";; RW Value";
}
// Parse throw types
if (type.find("GAUSPULL") != std::string::npos) fCalcType = kGausPull;
else fCalcType = kNoPull;
if (type.find("GAUSTHROW") != std::string::npos) fThrowType = kGausThrow;
+ else if (type.find("FLATTHROW") != std::string::npos) fThrowType = kFlatThrow;
else fThrowType = kNoThrow;
// Extra check to see if throws or pulls are turned off
if (type.find("NOPULL") != std::string::npos) fCalcType = kNoPull;
if (type.find("NOTHROW") != std::string::npos) fThrowType = kNoThrow;
}
//*******************************************************************************
void ParamPull::SetupHistograms(std::string input) {
//*******************************************************************************
// Extract Types from Input
fFileType = "";
const int nfiletypes = 4;
const std::string filetypes[nfiletypes] = {"FIT", "ROOT", "TXT", "DIAL"};
for (int i = 0; i < nfiletypes; i++) {
std::string tempTypes = filetypes[i] + ":";
if (input.find(tempTypes) != std::string::npos) {
fFileType = filetypes[i];
input.replace(input.find(tempTypes), tempTypes.size(), "");
break;
}
}
// Read Files
if (!fFileType.compare("FIT")) ReadFitFile(input);
else if (!fFileType.compare("ROOT")) ReadRootFile(input);
else if (!fFileType.compare("VECT")) ReadVectFile(input);
else if (!fFileType.compare("DIAL")) ReadDialInput(input);
else {
ERR(FTL) << "Unknown ParamPull Type: " << input << std::endl;
throw;
}
// Check Dials are all good
CheckDialsValid();
// Setup MC Histogram
fMCHist = (TH1D*) fDataHist->Clone();
fMCHist->Reset();
fMCHist->SetNameTitle( (fName + "_MC").c_str(),
(fName + " MC" + fPlotTitles).c_str() );
// If no Covar input make an uncorrelated one
if (!fCovar) {
fCovar = StatUtils::MakeDiagonalCovarMatrix(fDataHist, 1.0);
}
// If no types or limits are provided give them a default option
if (!fMinHist) {
fMinHist = (TH1D*) fDataHist->Clone();
fMinHist->SetNameTitle( (fName + "_min").c_str(),
(fName + " min" + fPlotTitles).c_str() );
for (int i = 0; i < fMinHist->GetNbinsX(); i++) {
// TODO (P.Stowell) Change this to a NULL system where limits are actually free!
fMinHist->SetBinContent(i + 1, fDataHist->GetBinContent(i + 1) - 1E6);
}
}
if (!fMaxHist) {
fMaxHist = (TH1D*) fDataHist->Clone();
fMaxHist->SetNameTitle( (fName + "_min").c_str(),
(fName + " min" + fPlotTitles).c_str() );
for (int i = 0; i < fMaxHist->GetNbinsX(); i++) {
fMaxHist->SetBinContent(i + 1, fDataHist->GetBinContent(i + 1) - 1E6);
}
}
// Set types from state, or to unknown
if (!fTypeHist) {
int deftype = -1;
if (fType.find("T2K") != std::string::npos) { deftype = kT2K; }
else if (fType.find("NEUT") != std::string::npos) { deftype = kNEUT; }
else if (fType.find("NIWG") != std::string::npos) { deftype = kNIWG; }
else if (fType.find("GENIE") != std::string::npos) { deftype = kGENIE; }
else if (fType.find("NORM") != std::string::npos) { deftype = kNORM; }
else if (fType.find("NUWRO") != std::string::npos) { deftype = kNUWRO; }
fTypeHist = new TH1I( (fName + "_type").c_str(),
(fName + " type" + fPlotTitles).c_str(),
fDataHist->GetNbinsX(), 0, fDataHist->GetNbinsX() );
for (int i = 0; i < fTypeHist->GetNbinsX(); i++) {
fTypeHist->SetBinContent(i + 1, deftype );
}
}
// Sort Covariances
fInvCovar = StatUtils::GetInvert(fCovar);
fDecomp = StatUtils::GetDecomp(fCovar);
// Create DataTrue for Throws
fDataTrue = (TH1D*) fDataHist->Clone();
fDataTrue->SetNameTitle( (fName + "_truedata").c_str(),
(fName + " truedata" + fPlotTitles).c_str() );
fDataOrig = (TH1D*) fDataHist->Clone();
fDataOrig->SetNameTitle( (fName + "_origdata").c_str(),
(fName + " origdata" + fPlotTitles).c_str() );
// Select only dials we want
if (!fDialSelection.empty()) {
(*fDataHist) = RemoveBinsNotInString(*fDataHist, fDialSelection);
}
}
//*******************************************************************************
TH1D ParamPull::RemoveBinsNotInString(TH1D hist, std::string mystr) {
//*******************************************************************************
// Make list of allowed bins
std::vector allowedbins;
for (int i = 0; i < hist.GetNbinsX(); i++) {
std::string syst = std::string(hist.GetXaxis()->GetBinLabel(i + 1));
if (mystr.find(syst) != std::string::npos) {
allowedbins.push_back(syst);
}
}
// Make new histogram
UInt_t nbins = allowedbins.size();
TH1D newhist = TH1D( hist.GetName(), hist.GetTitle(),
(Int_t)nbins, 0.0, (Double_t)nbins);
// Setup bins
for (UInt_t i = 0; i < nbins; i++) {
// Set Labels
newhist.GetXaxis()->SetBinLabel(i + 1, allowedbins[i].c_str());
// Copy Values
for (Int_t j = 0; j < hist.GetNbinsX(); j++) {
if (!allowedbins[i].compare(hist.GetXaxis()->GetBinLabel(j + 1))) {
newhist.SetBinContent(i + 1, hist.GetBinContent(j + 1) );
newhist.SetBinError(i + 1, hist.GetBinError(j + 1) );
}
}
}
return newhist;
}
//*******************************************************************************
TH1I ParamPull::RemoveBinsNotInString(TH1I hist, std::string mystr) {
//*******************************************************************************
// Make list of allowed bins
std::vector allowedbins;
for (int i = 0; i < hist.GetNbinsX(); i++) {
std::string syst = std::string(hist.GetXaxis()->GetBinLabel(i + 1));
if (mystr.find(syst) != std::string::npos) {
allowedbins.push_back(syst);
}
}
// Make new histogram
UInt_t nbins = allowedbins.size();
TH1I newhist = TH1I( hist.GetName(), hist.GetTitle(),
(Int_t)nbins, 0.0, (Int_t)nbins);
// Setup bins
for (UInt_t i = 0; i < nbins; i++) {
// Set Labels
newhist.GetXaxis()->SetBinLabel(i + 1, allowedbins[i].c_str());
// Copy Values
for (Int_t j = 0; j < hist.GetNbinsX(); j++) {
if (!allowedbins[i].compare(hist.GetXaxis()->GetBinLabel(j + 1))) {
newhist.SetBinContent(i + 1, hist.GetBinContent(j + 1) );
newhist.SetBinError(i + 1, hist.GetBinError(j + 1) );
}
}
}
return newhist;
}
//*******************************************************************************
void ParamPull::ReadFitFile(std::string input) {
//*******************************************************************************
TFile* tempfile = new TFile(input.c_str(), "READ");
// Read Data
fDataHist = (TH1D*) tempfile->Get("fit_dials_free");
if (!fDataHist) {
ERR(FTL) << "Can't find TH1D hist fit_dials in " << fName << std::endl;
ERR(FTL) << "File Entries:" << std::endl;
tempfile->ls();
throw;
}
fDataHist->SetDirectory(0);
fDataHist->SetNameTitle( (fName + "_data").c_str(),
(fName + " data" + fPlotTitles).c_str() );
// Read Covar
TH2D* tempcov = (TH2D*) tempfile->Get("covariance_free");
if (!tempcov) {
ERR(FTL) << "Can't find TH2D covariance_free in " << fName << std::endl;
ERR(FTL) << "File Entries:" << std::endl;
tempfile->ls();
throw;
}
// Setup Covar
int nbins = fDataHist->GetNbinsX();
fCovar = new TMatrixDSym( nbins );
for (int i = 0; i < nbins; i++) {
for (int j = 0; j < nbins; j++) {
(*fCovar)(i, j) = tempcov->GetBinContent(i + 1, j + 1);
}
}
return;
}
//*******************************************************************************
void ParamPull::ReadRootFile(std::string input) {
//*******************************************************************************
std::vector inputlist = GeneralUtils::ParseToStr(input, ";");
// Check all given
if (inputlist.size() < 2) {
ERR(FTL) << "Covar supplied in 'ROOT' format should have 3 semi-colon seperated entries!" << std::endl
<< "ROOT:filename;histname[;covarname]" << std::endl;
ERR(FTL) << "histname = TH1D, covarname = TH2D" << std::endl;
throw;
}
// Get Entries
std::string filename = inputlist[0];
LOG(DEB) << filename << std::endl;
std::string histname = inputlist[1];
LOG(DEB) << histname << std::endl;
LOG(DEB) << input << std::endl;
// Read File
TFile* tempfile = new TFile(filename.c_str(), "READ");
if (tempfile->IsZombie()) {
ERR(FTL) << "Can't find file in " << fName << std::endl;
ERR(FTL) << "location = " << filename << std::endl;
throw;
}
// Read Hist
fDataHist = (TH1D*) tempfile->Get(histname.c_str());
if (!fDataHist) {
ERR(FTL) << "Can't find TH1D hist " << histname << " in " << fName << std::endl;
ERR(FTL) << "File Entries:" << std::endl;
tempfile->ls();
throw;
}
fDataHist->SetDirectory(0);
fDataHist->SetNameTitle( (fName + "_data").c_str(),
(fName + " data" + fPlotTitles).c_str() );
LOG(DEB) << "READING COVAR" << std::endl;
// Read Covar
if (inputlist.size() > 2) {
std::string covarname = inputlist[2];
LOG(DEB) << "COVARNAME = " << covarname << std::endl;
TH2D* tempcov = (TH2D*) tempfile->Get(covarname.c_str());
if (!tempcov) {
ERR(FTL) << "Can't find TH2D covar " << covarname << " in " << fName << std::endl;
ERR(FTL) << "File Entries:" << std::endl;
tempfile->ls();
throw;
}
// Setup Covar
int nbins = fDataHist->GetNbinsX();
fCovar = new TMatrixDSym( nbins );
for (int i = 0; i < nbins; i++) {
for (int j = 0; j < nbins; j++) {
(*fCovar)(i, j) = tempcov->GetBinContent(i + 1, j + 1);
}
}
// Uncorrelated
} else {
LOG(SAM) << "No Covar provided so using diagonal errors for "
<< fName << std::endl;
fCovar = NULL;
}
}
//*******************************************************************************
void ParamPull::ReadVectFile(std::string input) {
//*******************************************************************************
std::vector inputlist = GeneralUtils::ParseToStr(input, ";");
if (inputlist.size() < 4) {
ERR(FTL) << "Need 3 inputs for vector input in " << fName << std::endl;
ERR(FTL) << "Inputs: " << input << std::endl;
throw;
}
// Open File
std::string rootname = inputlist[0];
TFile* tempfile = new TFile(rootname.c_str(), "READ");
if (tempfile->IsZombie()) {
ERR(FTL) << "Can't find file in " << fName << std::endl;
ERR(FTL) << "location = " << rootname << std::endl;
throw;
}
// Get Name
std::string tagname = inputlist[1];
// TVector dialtags = tempfile->Get(tagname.c_str());
// if (!dialtags){
// ERR(FTL) << "Can't find list of dial names!" << std::endl;
// }
// Get Values
std::string valuename = inputlist[2];
TVectorD* dialvals = (TVectorD*)tempfile->Get(valuename.c_str());
if (!dialvals) {
ERR(FTL) << "Can't find dial values" << std::endl;
}
// Get Matrix
std::string matrixname = inputlist[3];
TMatrixD* matrixvals = (TMatrixD*)tempfile->Get(matrixname.c_str());
if (!matrixvals) {
ERR(FTL) << "Can't find matirx values" << std::endl;
}
// Get Types
if (inputlist.size() > 4) {
std::string typesname = inputlist[3];
}
// Get Minimum
if (inputlist.size() > 5) {
std::string minname = inputlist[4];
}
// Get Maximum
if (inputlist.size() > 6) {
std::string maxname = inputlist[5];
}
}
//*******************************************************************************
void ParamPull::ReadDialInput(std::string input) {
//*******************************************************************************
std::vector inputlist = GeneralUtils::ParseToStr(input, ";");
if (inputlist.size() < 3) {
ERR(FTL) << "Need 3 inputs for dial input in " << fName << std::endl;
ERR(FTL) << "Inputs: " << input << std::endl;
throw;
}
std::vector inputvals = GeneralUtils::ParseToDbl(input, ";");
std::string dialname = inputlist[0];
double val = inputvals[1];
double err = inputvals[2];
fDataHist = new TH1D( (fName + "_data").c_str(),
(fName + "_data" + fPlotTitles).c_str(), 1, 0, 1);
fDataHist->SetBinContent(1, val);
fDataHist->SetBinError(1, err);
fDataHist->GetXaxis()->SetBinLabel(1, dialname.c_str());
+
+ fLimitHist = new TH1D( (fName + "_limits").c_str(),
+ (fName + "_limits" + fPlotTitles).c_str(), 1, 0, 1);
+ fLimitHist->Reset();
+ if (inputvals.size() > 4) {
+ fLimitHist->SetBinContent(1, (inputvals[3] + inputvals[4]) / 2.0);
+ fLimitHist->SetBinError(1, (inputvals[4] - inputvals[3]) / 2.0);
+ }
+
fCovar = NULL;
}
//*******************************************************************************
std::map ParamPull::GetAllDials() {
//*******************************************************************************
std::map dialtypemap;
for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
std::string name = fDataHist->GetXaxis()->GetBinLabel(i + 1);
int type = fTypeHist->GetBinContent(i + 1);
dialtypemap[name] = type;
}
return dialtypemap;
}
//*******************************************************************************
bool ParamPull::CheckDialsValid() {
//*******************************************************************************
return true;
std::string helpstring = "";
for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
std::string name = std::string(fDataHist->GetXaxis()->GetBinLabel(i + 1));
// If dial exists its all good
if (FitBase::GetRW()->DialIncluded(name)) continue;
// If it doesn't but its a sample norm also continue
if (name.find("_norm") != std::string::npos) {
ERR(WRN) << "Norm dial included in covar but not set in FitWeight." << std::endl;
ERR(WRN) << "Assuming its a sample norm and skipping..." << std::endl;
}
// Dial unknown so print a help statement
std::ostringstream tempstr;
tempstr << "unknown_parameter " << name << " "
<< fDataHist->GetBinContent(i + 1) << " "
<< fDataHist->GetBinContent(i + 1) - fDataHist->GetBinError(i + 1) << " "
<< fDataHist->GetBinContent(i + 1) + fDataHist->GetBinError(i + 1) << " "
<< fDataHist->GetBinError(i + 1) << " ";
if (!fType.empty()) tempstr << fType << std::endl;
else tempstr << "FREE" << std::endl;
helpstring += tempstr.str();
}
// Show statement before failing
if (!helpstring.empty()) {
ERR(WRN) << "Dial(s) included in covar but not set in FitWeight." << std::endl
<< "ParamPulls needs to know how you want it to be treated." << std::endl
<< "Include the following lines into your card:" << std::endl;
ERR(WRN) << helpstring << std::endl;
throw;
return false;
} else {
return true;
}
}
//*******************************************************************************
void ParamPull::Reconfigure() {
//*******************************************************************************
FitWeight* rw = FitBase::GetRW();
// Get Dial Names that are valid
std::vector namevec = rw->GetDialNames();
std::vector valuevec = rw->GetDialValues();
// Set Bin Values from RW
for (UInt_t i = 0; i < namevec.size(); i++) {
// Loop over bins and check name matches
std::string syst = namevec.at(i);
double systval = valuevec.at(i);
std::vector allsyst = GeneralUtils::ParseToStr(syst, ",");
// Proper Reconf using RW
for (int j = 0; j < fMCHist->GetNbinsX(); j++) {
// Search for the name of this bin in the corrent dial
std::string binname = std::string(fMCHist->GetXaxis()->GetBinLabel(j + 1) );
// Check Full Name
if (!syst.compare(binname.c_str())) {
fMCHist->SetBinContent(j + 1, systval);
break;
}
std::vector splitbinname = GeneralUtils::ParseToStr(binname, ",");
for (size_t l = 0; l < splitbinname.size(); l++) {
std::string singlebinname = splitbinname[l];
for (size_t k = 0; k < allsyst.size(); k++) {
if (!allsyst[k].compare(singlebinname.c_str())) {
fMCHist->SetBinContent(j + 1, systval);
}
}
}
}
}
return;
};
//*******************************************************************************
void ParamPull::ResetToy(void) {
//*******************************************************************************
if (fDataHist) delete fDataHist;
LOG(DEB) << "Resetting toy" << std::endl;
LOG(DEB) << fDataTrue << std::endl;
fDataHist = (TH1D*)fDataTrue->Clone();
LOG(DEB) << "Setting name" << std::endl;
fDataHist->SetNameTitle( (fName + "_data").c_str(),
(fName + " data" + fPlotTitles).c_str() );
}
//*******************************************************************************
void ParamPull::SetFakeData(std::string fakeinput) {
//*******************************************************************************
// Set from MC Setting
if (!fakeinput.compare("MC")) {
// Copy MC into data
if (fDataHist) delete fDataHist;
fDataHist = (TH1D*)fMCHist->Clone();
fDataHist->SetNameTitle( (fName + "_data").c_str(),
(fName + " fakedata" + fPlotTitles).c_str() );
// Copy original data errors
for (int i = 0; i < fDataOrig->GetNbinsX(); i++) {
fDataHist->SetBinError(i + 1, fDataOrig->GetBinError(i + 1) );
}
// Make True Toy Central Value Hist
fDataTrue = (TH1D*) fDataHist->Clone();
fDataTrue->SetNameTitle( (fName + "_truedata").c_str(),
(fName + " truedata" + fPlotTitles).c_str() );
} else {
ERR(FTL) << "Trying to set fake data for ParamPulls not from MC!" << std::endl;
ERR(FTL) << "Not currently implemented.." << std::endl;
throw;
}
}
//*******************************************************************************
void ParamPull::RemoveFakeData() {
//*******************************************************************************
delete fDataHist;
fDataHist = (TH1D*)fDataOrig->Clone();
fDataHist->SetNameTitle( (fName + "_data").c_str(),
(fName + " data" + fPlotTitles).c_str() );
fDataTrue = (TH1D*) fDataHist->Clone();
fDataTrue->SetNameTitle( (fName + "_truedata").c_str(),
(fName + " truedata" + fPlotTitles).c_str() );
}
//*******************************************************************************
double ParamPull::GetLikelihood() {
//*******************************************************************************
double like = 0.0;
switch (fCalcType) {
// Gaussian Calculation with correlations
case kGausPull:
like = StatUtils::GetChi2FromCov(fDataHist, fMCHist, fInvCovar, NULL);
like *= 1E-76;
break;
// Default says this has no pull
case kNoThrow:
default:
like = 0.0;
break;
}
LOG(DEB) << "Likelihood = " << like << " " << fCalcType << std::endl;
return like;
};
//*******************************************************************************
int ParamPull::GetNDOF() {
//*******************************************************************************
int ndof = 0;
if (fCalcType != kNoThrow) {
- ndof += fDataHist->GetNbinsX();
+ ndof = fDataHist->GetNbinsX();
}
return ndof;
};
//*******************************************************************************
void ParamPull::ThrowCovariance() {
//*******************************************************************************
// Reset toy for throw
ResetToy();
- LOG(DEB) << "Toy Reset " << std::endl;
+ LOG(FIT) << "Creating new toy dataset" << std::endl;
// Generate random Gaussian throws
std::vector randthrows;
for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
double randtemp = 0.0;
switch (fThrowType) {
// Gaussian Throws
case kGausThrow:
randtemp = gRandom->Gaus(0.0, 1.0);
break;
+ // Uniform Throws
+ case kFlatThrow:
+ randtemp = gRandom->Uniform(0.0, 1.0);
+ if (fLimitHist) {
+ randtemp = fLimitHist->GetBinContent(i + 1) + fLimitHist->GetBinError(i + 1) * ( randtemp * 2 - 1 );
+ }
+ break;
+
// No Throws (DEFAULT)
default:
break;
}
randthrows.push_back(randtemp);
}
// Create Bin Modifications
double totalres = 0.0;
for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
// Calc Bin Mod
double binmod = 0.0;
- for (int j = 0; j < fDataHist->GetNbinsX(); j++) {
- // std::cout << "DECOMP " << j << " " << i << " " << randthrows.at(j) << std::endl;
- binmod += (*fDecomp)(j, i) * randthrows.at(j);
+
+ if (fThrowType == kGausThrow) {
+ for (int j = 0; j < fDataHist->GetNbinsX(); j++) {
+ binmod += (*fDecomp)(j, i) * randthrows.at(j);
+ }
+ } else if (fThrowType == kFlatThrow) {
+ binmod = randthrows.at(i) - fDataHist->GetBinContent(i + 1);
}
+
// Add up fraction dif
totalres += binmod;
// Add to current data
fDataHist->SetBinContent(i + 1, fDataHist->GetBinContent(i + 1) + binmod);
}
// Rename
fDataHist->SetNameTitle( (fName + "_data").c_str(),
(fName + " toydata" + fPlotTitles).c_str() );
- // Print Status
- LOG(REC) << "Created new toy histogram. Total Dif = "
- << totalres << std::endl;
+ // Check Limits
+ if (fLimitHist) {
+ for (int i = 0; i < fLimitHist->GetNbinsX(); i++) {
+ if (fLimitHist->GetBinError(i + 1) == 0.0) continue;
+ if (fDataHist->GetBinContent(i + 1) > fLimitHist->GetBinContent(i + 1) + fLimitHist->GetBinError(i + 1) ||
+ fDataHist->GetBinContent(i + 1) < fLimitHist->GetBinContent(i + 1) - fLimitHist->GetBinError(i + 1)) {
+ this->ThrowCovariance();
+ }
+ }
+ }
+
return;
};
//*******************************************************************************
TH2D ParamPull::GetCovar() {
//*******************************************************************************
TH2D tempCov = TH2D(*fInvCovar);
for (int i = 0; i < tempCov.GetNbinsX(); i++) {
tempCov.GetXaxis()->SetBinLabel(i + 1, fDataHist->GetXaxis()->GetBinLabel(i + 1));
tempCov.GetYaxis()->SetBinLabel(i + 1, fDataHist->GetXaxis()->GetBinLabel(i + 1));
}
tempCov.SetNameTitle((fName + "_INVCOV").c_str(),
(fName + " InvertedCovariance;Dials;Dials").c_str());
return tempCov;
}
//*******************************************************************************
TH2D ParamPull::GetFullCovar() {
//*******************************************************************************
TH2D tempCov = TH2D(*fCovar);
for (int i = 0; i < tempCov.GetNbinsX(); i++) {
tempCov.GetXaxis()->SetBinLabel(i + 1, fDataHist->GetXaxis()->GetBinLabel(i + 1));
tempCov.GetYaxis()->SetBinLabel(i + 1, fDataHist->GetXaxis()->GetBinLabel(i + 1));
}
tempCov.SetNameTitle((fName + "_COV").c_str(),
(fName + " Covariance;Dials;Dials").c_str());
return tempCov;
}
//*******************************************************************************
TH2D ParamPull::GetDecompCovar() {
//*******************************************************************************
TH2D tempCov = TH2D(*fCovar);
for (int i = 0; i < tempCov.GetNbinsX(); i++) {
tempCov.GetXaxis()->SetBinLabel(i + 1, fDataHist->GetXaxis()->GetBinLabel(i + 1));
tempCov.GetYaxis()->SetBinLabel(i + 1, fDataHist->GetXaxis()->GetBinLabel(i + 1));
}
tempCov.SetNameTitle((fName + "_DEC").c_str(),
(fName + " Decomposition;Dials;Dials").c_str());
return tempCov;
}
//*******************************************************************************
void ParamPull::Write(std::string writeoptt) {
//*******************************************************************************
fDataHist->Write();
fMCHist->Write();
-
+ if (fLimitHist) {
+ fLimitHist->Write();
+ }
GetCovar().Write();
GetFullCovar().Write();
GetDecompCovar().Write();
return;
};
diff --git a/src/FitBase/ParamPull.h b/src/FitBase/ParamPull.h
index 22e707d..7da9eac 100644
--- a/src/FitBase/ParamPull.h
+++ b/src/FitBase/ParamPull.h
@@ -1,188 +1,193 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#ifndef PARAM_PULL_H_SEEN
#define PARAM_PULL_H_SEEN
/*!
* \addtogroup FitBase
* @{
*/
#include
#include
#include
#include
#include
#include
#include
#include
// ROOT includes
#include
#include
#include
#include
#include
#include
#include
// Fit Includes
#include "PlotUtils.h"
#include "StatUtils.h"
#include "FitWeight.h"
#include "FitLogger.h"
#include "EventManager.h"
#include "TVector.h"
using namespace std;
//! Enums to allow Non Guassian Pulls in Future
enum FitPullTypes {
kUnknownPull = -1,
kNoPull = 0,
kGausPull = 1
};
enum FitThrowTypes {
kUnknownThrow = -1,
kNoThrow = 0,
- kGausThrow = 1
+ kGausThrow = 1,
+ kFlatThrow = 2
};
//! Used to produce gaussian penalty terms in the fit.
class ParamPull {
public:
//! Default Constructor
ParamPull(std::string name, std::string inputfile, std::string type, std::string dials="");
//! Default destructor
virtual ~ParamPull(void) {};
// Set dial types (DEFAULT,ABS,FRAC)
void SetType(std::string type);
// Setup Histogram inputs (from previous fit file, or ROOT file)
void SetupHistograms(std::string input);
// Read a previous NUISANCE file
void ReadFitFile(std::string input);
// Read a ROOT file with any histograms in
void ReadRootFile(std::string input);
// Read Text file
void ReadVectFile(std::string input);
// Read a single dial central value and error
void ReadDialInput(std::string input);
//! Reconfigure function reads in current RW engine dials and sets their value to MC
void Reconfigure(void);
//! Get likelihood given the current values
double GetLikelihood(void);
//! Get NDOF if used in likelihoods
int GetNDOF(void);
// Get Covariance Matrices as plots
TH2D GetCovar(void);
TH2D GetFullCovar(void);
TH2D GetDecompCovar(void);
// Get Covariance Matrices
inline TMatrixDSym GetCovarMatrix (void) const { return *fInvCovar; };
inline TMatrixDSym GetFullCovarMatrix (void) const { return *fCovar; };
inline TMatrixDSym GetDecompCovarMatrix (void) const { return *fDecomp; };
//! Save the histograms
void Write(std::string writeopt="");
//! Throw the dial values using the current covariance. Useful for parameter throws.
void ThrowCovariance(void);
//! Compare dials to RW
bool CheckDialsValid(void);
//! Reset toy data back to original data
void ResetToy(void);
//! Read fake data from MC
void SetFakeData(std::string fakeinput);
//! Reset fake data back to original data (before fake data or throws)
void RemoveFakeData();
// Get Functions
inline std::string GetName (void) const { return fName; };
inline std::string GetInput (void) const { return fInput; };
inline std::string GetType (void) const { return fType; };
inline std::string GetFileType (void) const { return fFileType; };
inline std::string GetDialOptions (void) const { return fDialOptions; };
std::map GetAllDials();
inline TH1D GetDataHist (void) const { return *fDataHist; };
inline TH1D GetDataTrue (void) const { return *fDataTrue; };
inline TH1D GetDataOrig (void) const { return *fDataOrig; };
inline TH1D GetMCHist (void) const { return *fMCHist; };
inline TH1D GetMaxHist (void) const { return *fMaxHist; };
inline TH1D GetMinHist (void) const { return *fMinHist; };
inline TH1I GetDialTypes (void) const { return *fTypeHist; };
+ inline TH1D GetLimitHist (void) const { return *fLimitHist; };
+
private:
TH1D RemoveBinsNotInString(TH1D hist, std::string mystr);
TH1I RemoveBinsNotInString(TH1I hist, std::string mystr);
std::string fName; //!< Pull Name
std::string fInput; //!< Pull input string
std::string fType; //!< Pull options type
std::string fFileType; //!< Pull input file types
std::string fPlotTitles; //! Axis format
std::string fDialOptions; //!< Dial handling options
std::string fDialSelection; //!< Dial Selection
TH1D* fMCHist; //!< Current MC Histogram
TH1D* fDataHist; //!< Current data Histogram
TH1D* fDataTrue; //!< True Data (before histogram throws)
TH1D* fDataOrig; //!< Orig Data (without toys or fake data)
TH1D* fMaxHist; //!< Maximum limit on MC/Data
TH1D* fMinHist; //!< Maximum limit on MC/Data
TH1I* fTypeHist; //!< Dial Types
int fCalcType; //!< Method to calculate likelihood
int fThrowType; //!< Method to calculate throws
TMatrixDSym* fCovar; //!< Covariance
TMatrixDSym* fInvCovar; //!< Inverted Covariance
TMatrixDSym* fDecomp; //!< Decomposition
+
+ TH1D* fLimitHist;
};
// Class TypeDefs
typedef std::list::const_iterator PullListConstIter;
typedef std::list::iterator PullListIter;
typedef std::vector::const_iterator PullVectConstIter;
typedef std::vector::iterator PullVectIter;
/*! @} */
#endif
diff --git a/src/Routines/SystematicRoutines.cxx b/src/Routines/SystematicRoutines.cxx
index 4c846a5..d33b35e 100755
--- a/src/Routines/SystematicRoutines.cxx
+++ b/src/Routines/SystematicRoutines.cxx
@@ -1,1516 +1,1517 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#include "SystematicRoutines.h"
void SystematicRoutines::Init(){
fInputFile = "";
fInputRootFile = NULL;
fOutputFile = "";
fOutputRootFile = NULL;
fCovar = fCovarFree = NULL;
fCorrel = fCorrelFree = NULL;
fDecomp = fDecompFree = NULL;
fStrategy = "ErrorBands";
fRoutines.clear();
fRoutines.push_back("ErrorBands");
fCardFile = "";
fFakeDataInput = "";
fSampleFCN = NULL;
fAllowedRoutines = ("ErrorBands,PlotLimits");
};
SystematicRoutines::~SystematicRoutines(){
};
SystematicRoutines::SystematicRoutines(int argc, char* argv[]){
// Initialise Defaults
Init();
nuisconfig configuration = Config::Get();
// Default containers
std::string cardfile = "";
std::string maxevents = "-1";
int errorcount = 0;
int verbocount = 0;
std::vector xmlcmds;
std::vector configargs;
fNThrows = 250;
fStartThrows = 0;
fThrowString = "";
// Make easier to handle arguments.
std::vector args = GeneralUtils::LoadCharToVectStr(argc, argv);
ParserUtils::ParseArgument(args, "-c", fCardFile, true);
ParserUtils::ParseArgument(args, "-o", fOutputFile, false, false);
ParserUtils::ParseArgument(args, "-n", maxevents, false, false);
ParserUtils::ParseArgument(args, "-f", fStrategy, false, false);
ParserUtils::ParseArgument(args, "-d", fFakeDataInput, false, false);
ParserUtils::ParseArgument(args, "-s", fStartThrows, false, false);
ParserUtils::ParseArgument(args, "-t", fNThrows, false, false);
ParserUtils::ParseArgument(args, "-p", fThrowString, false, false);
ParserUtils::ParseArgument(args, "-i", xmlcmds);
ParserUtils::ParseArgument(args, "-q", configargs);
ParserUtils::ParseCounter(args, "e", errorcount);
ParserUtils::ParseCounter(args, "v", verbocount);
ParserUtils::CheckBadArguments(args);
// Add extra defaults if none given
if (fCardFile.empty() and xmlcmds.empty()) {
ERR(FTL) << "No input supplied!" << std::endl;
throw;
}
if (fOutputFile.empty() and !fCardFile.empty()) {
fOutputFile = fCardFile + ".root";
ERR(WRN) << "No output supplied so saving it to: " << fOutputFile << std::endl;
} else if (fOutputFile.empty()) {
ERR(FTL) << "No output file or cardfile supplied!" << std::endl;
throw;
}
// Configuration Setup =============================
// Check no comp key is available
if (Config::Get().GetNodes("nuiscomp").empty()) {
fCompKey = Config::Get().CreateNode("nuiscomp");
} else {
fCompKey = Config::Get().GetNodes("nuiscomp")[0];
}
if (!fCardFile.empty()) fCompKey.AddS("cardfile", fCardFile);
if (!fOutputFile.empty()) fCompKey.AddS("outputfile", fOutputFile);
if (!fStrategy.empty()) fCompKey.AddS("strategy", fStrategy);
// Load XML Cardfile
configuration.LoadConfig( fCompKey.GetS("cardfile"), "");
// Add CMD XML Structs
for (size_t i = 0; i < xmlcmds.size(); i++) {
configuration.AddXMLLine(xmlcmds[i]);
}
// Add Config Args
for (size_t i = 0; i < configargs.size(); i++) {
configuration.OverrideConfig(configargs[i]);
}
if (maxevents.compare("-1")){
configuration.OverrideConfig("MAXEVENTS=" + maxevents);
}
// Finish configuration XML
configuration.FinaliseConfig(fCompKey.GetS("outputfile") + ".xml");
// Add Error Verbo Lines
verbocount += Config::Get().GetParI("VERBOSITY");
errorcount += Config::Get().GetParI("ERROR");
std::cout << "[ NUISANCE ]: Setting VERBOSITY=" << verbocount << std::endl;
std::cout << "[ NUISANCE ]: Setting ERROR=" << errorcount << std::endl;
SETVERBOSITY(verbocount);
// Proper Setup
if (fStrategy.find("ErrorBands") != std::string::npos ||
fStrategy.find("MergeErrors") != std::string::npos){
fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
}
// fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
SetupSystematicsFromXML();
SetupCovariance();
SetupRWEngine();
SetupFCN();
GetCovarFromFCN();
// Run();
return;
};
void SystematicRoutines::SetupSystematicsFromXML(){
LOG(FIT) << "Setting up nuismin" << std::endl;
// Setup Parameters ------------------------------------------
std::vector parkeys = Config::QueryKeys("parameter");
if (!parkeys.empty()) {
LOG(FIT) << "Number of parameters : " << parkeys.size() << std::endl;
}
for (size_t i = 0; i < parkeys.size(); i++) {
nuiskey key = parkeys.at(i);
// Check for type,name,nom
if (!key.Has("type")) {
ERR(FTL) << "No type given for parameter " << i << std::endl;
throw;
} else if (!key.Has("name")) {
ERR(FTL) << "No name given for parameter " << i << std::endl;
throw;
} else if (!key.Has("nominal")) {
ERR(FTL) << "No nominal given for parameter " << i << std::endl;
throw;
}
// Get Inputs
std::string partype = key.GetS("type");
std::string parname = key.GetS("name");
double parnom = key.GetD("nominal");
double parlow = parnom - 1;
double parhigh = parnom + 1;
double parstep = 1;
// Override if state not given
if (!key.Has("state")){
key.SetS("state","FIX");
}
std::string parstate = key.GetS("state");
// Extra limits
if (key.Has("low")) {
parlow = key.GetD("low");
parhigh = key.GetD("high");
parstep = key.GetD("step");
LOG(FIT) << "Read " << partype << " : "
<< parname << " = "
<< parnom << " : "
<< parlow << " < p < " << parhigh
<< " : " << parstate << std::endl;
} else {
LOG(FIT) << "Read " << partype << " : "
<< parname << " = "
<< parnom << " : "
<< parstate << std::endl;
}
// Run Parameter Conversion if needed
if (parstate.find("ABS") != std::string::npos) {
parnom = FitBase::RWAbsToSigma( partype, parname, parnom );
parlow = FitBase::RWAbsToSigma( partype, parname, parlow );
parhigh = FitBase::RWAbsToSigma( partype, parname, parhigh );
parstep = FitBase::RWAbsToSigma( partype, parname, parstep );
} else if (parstate.find("FRAC") != std::string::npos) {
parnom = FitBase::RWFracToSigma( partype, parname, parnom );
parlow = FitBase::RWFracToSigma( partype, parname, parlow );
parhigh = FitBase::RWFracToSigma( partype, parname, parhigh );
parstep = FitBase::RWFracToSigma( partype, parname, parstep );
}
// Push into vectors
fParams.push_back(parname);
fTypeVals[parname] = FitBase::ConvDialType(partype);;
fStartVals[parname] = parnom;
fCurVals[parname] = parnom;
fErrorVals[parname] = 0.0;
fStateVals[parname] = parstate;
bool fixstate = parstate.find("FIX") != std::string::npos;
fFixVals[parname] = fixstate;
fStartFixVals[parname] = fFixVals[parname];
fMinVals[parname] = parlow;
fMaxVals[parname] = parhigh;
fStepVals[parname] = parstep;
}
// Setup Samples ----------------------------------------------
std::vector samplekeys = Config::QueryKeys("sample");
if (!samplekeys.empty()) {
LOG(FIT) << "Number of samples : " << samplekeys.size() << std::endl;
}
for (size_t i = 0; i < samplekeys.size(); i++) {
nuiskey key = samplekeys.at(i);
// Get Sample Options
std::string samplename = key.GetS("name");
std::string samplefile = key.GetS("input");
std::string sampletype =
key.Has("type") ? key.GetS("type") : "DEFAULT";
double samplenorm =
key.Has("norm") ? key.GetD("norm") : 1.0;
// Print out
LOG(FIT) << "Read sample info " << i << " : "
<< samplename << std::endl
<< "\t\t input -> " << samplefile << std::endl
<< "\t\t state -> " << sampletype << std::endl
<< "\t\t norm -> " << samplenorm << std::endl;
// If FREE add to parameters otherwise continue
if (sampletype.find("FREE") == std::string::npos) {
continue;
}
// Form norm dial from samplename + sampletype + "_norm";
std::string normname = samplename + "_norm";
// Check normname not already present
if (fTypeVals.find(normname) != fTypeVals.end()) {
continue;
}
// Add new norm dial to list if its passed above checks
fParams.push_back(normname);
fTypeVals[normname] = kNORM;
fStateVals[normname] = sampletype;
fCurVals[normname] = samplenorm;
fErrorVals[normname] = 0.0;
fMinVals[normname] = 0.1;
fMaxVals[normname] = 10.0;
fStepVals[normname] = 0.5;
bool state = sampletype.find("FREE") == std::string::npos;
fFixVals[normname] = state;
fStartFixVals[normname] = state;
}
// Setup Fake Parameters -----------------------------
std::vector fakekeys = Config::QueryKeys("fakeparameter");
if (!fakekeys.empty()) {
LOG(FIT) << "Number of fake parameters : " << fakekeys.size() << std::endl;
}
for (size_t i = 0; i < fakekeys.size(); i++) {
nuiskey key = fakekeys.at(i);
// Check for type,name,nom
if (!key.Has("name")) {
ERR(FTL) << "No name given for fakeparameter " << i << std::endl;
throw;
} else if (!key.Has("nom")) {
ERR(FTL) << "No nominal given for fakeparameter " << i << std::endl;
throw;
}
// Get Inputs
std::string parname = key.GetS("name");
double parnom = key.GetD("nom");
// Push into vectors
fFakeVals[parname] = parnom;
}
}
void SystematicRoutines::ReadCard(std::string cardfile){
// Read cardlines into vector
std::vector cardlines = GeneralUtils::ParseFileToStr(cardfile,"\n");
FitPar::Config().cardLines = cardlines;
// Read Samples first (norm params can be overridden)
int linecount = 0;
for (std::vector::iterator iter = cardlines.begin();
iter != cardlines.end(); iter++){
std::string line = (*iter);
linecount++;
// Skip Empties
if (line.empty()) continue;
if (line.c_str()[0] == '#') continue;
// Read Valid Samples
int samstatus = ReadSamples(line);
// Show line if bad to help user
if (samstatus == kErrorStatus) {
ERR(FTL) << "Bad Input in cardfile " << fCardFile
<< " at line " << linecount << "!" << std::endl;
LOG(FIT) << line << std::endl;
throw;
}
}
// Read Parameters second
linecount = 0;
for (std::vector::iterator iter = cardlines.begin();
iter != cardlines.end(); iter++){
std::string line = (*iter);
linecount++;
// Skip Empties
if (line.empty()) continue;
if (line.c_str()[0] == '#') continue;
// Try Parameter Reads
int parstatus = ReadParameters(line);
int fakstatus = ReadFakeDataPars(line);
// Show line if bad to help user
if (parstatus == kErrorStatus ||
fakstatus == kErrorStatus ){
ERR(FTL) << "Bad Parameter Input in cardfile " << fCardFile
<< " at line " << linecount << "!" << std::endl;
LOG(FIT) << line << std::endl;
throw;
}
}
return;
};
int SystematicRoutines::ReadParameters(std::string parstring){
std::string inputspec = "RW Dial Inputs Syntax \n"
"free input w/ limits: TYPE NAME START MIN MAX STEP [STATE] \n"
"fix input: TYPE NAME VALUE [STATE] \n"
"free input w/o limits: TYPE NAME START FREE,[STATE] \n"
"Allowed Types: \n"
"neut_parameter,niwg_parameter,t2k_parameter,"
"nuwro_parameter,gibuu_parameter";
// Check sample input
if (parstring.find("parameter") == std::string::npos) return kGoodStatus;
// Parse inputs
std::vector strvct = GeneralUtils::ParseToStr(parstring, " ");
// Skip if comment or parameter somewhere later in line
if (strvct[0].c_str()[0] == '#' ||
strvct[0].find("parameter") == std::string::npos){
return kGoodStatus;
}
// Check length
if (strvct.size() < 3){
ERR(FTL) << "Input rw dials need to provide at least 3 inputs." << std::endl;
std::cout << inputspec << std::endl;
return kErrorStatus;
}
// Setup default inputs
std::string partype = strvct[0];
std::string parname = strvct[1];
double parval = GeneralUtils::StrToDbl(strvct[2]);
double minval = parval - 1.0;
double maxval = parval + 1.0;
double stepval = 1.0;
std::string state = "FIX"; //[DEFAULT]
// Check Type
if (FitBase::ConvDialType(partype) == kUNKNOWN){
ERR(FTL) << "Unknown parameter type! " << partype << std::endl;
std::cout << inputspec << std::endl;
return kErrorStatus;
}
// Check Parameter Name
if (FitBase::GetDialEnum(partype, parname) == -1){
ERR(FTL) << "Bad RW parameter name! " << partype << " " << parname << std::endl;
std::cout << inputspec << std::endl;
return kErrorStatus;
}
// Option Extra (No Limits)
if (strvct.size() == 4){
state = strvct[3];
}
// Check for weirder inputs
if (strvct.size() > 4 && strvct.size() < 6){
ERR(FTL) << "Provided incomplete limits for " << parname << std::endl;
std::cout << inputspec << std::endl;
return kErrorStatus;
}
// Option Extra (With limits and steps)
if (strvct.size() >= 6){
minval = GeneralUtils::StrToDbl(strvct[3]);
maxval = GeneralUtils::StrToDbl(strvct[4]);
stepval = GeneralUtils::StrToDbl(strvct[5]);
}
// Option Extra (dial state after limits)
if (strvct.size() == 7){
state = strvct[6];
}
// Run Parameter Conversion if needed
if (state.find("ABS") != std::string::npos){
parval = FitBase::RWAbsToSigma( partype, parname, parval );
minval = FitBase::RWAbsToSigma( partype, parname, minval );
maxval = FitBase::RWAbsToSigma( partype, parname, maxval );
stepval = FitBase::RWAbsToSigma( partype, parname, stepval );
} else if (state.find("FRAC") != std::string::npos){
parval = FitBase::RWFracToSigma( partype, parname, parval );
minval = FitBase::RWFracToSigma( partype, parname, minval );
maxval = FitBase::RWFracToSigma( partype, parname, maxval );
stepval = FitBase::RWFracToSigma( partype, parname, stepval );
}
// Check no repeat params
if (std::find(fParams.begin(), fParams.end(), parname) != fParams.end()){
ERR(FTL) << "Duplicate parameter names given for " << parname << std::endl;
throw;
}
// Setup Containers
fParams.push_back(parname);
fTypeVals[parname] = FitBase::ConvDialType(partype);
fStartVals[parname] = parval;
fCurVals[parname] = fStartVals[parname];
fErrorVals[parname] = 0.0;
fStateVals[parname] = state;
bool fixstate = state.find("FIX") != std::string::npos;
fFixVals[parname] = fixstate;
fStartFixVals[parname] = fFixVals[parname];
fMinVals[parname] = minval;
fMaxVals[parname] = maxval;
fStepVals[parname] = stepval;
// Print the parameter
LOG(MIN) << "Read Parameter " << parname << " " << parval << " "
<< minval << " " << maxval << " "
<< stepval << " " << state << std::endl;
// Tell reader its all good
return kGoodStatus;
}
//*******************************************
int SystematicRoutines::ReadFakeDataPars(std::string parstring){
//******************************************
std::string inputspec = "Fake Data Dial Inputs Syntax \n"
"fake value: fake_parameter NAME VALUE \n"
"Name should match dialnames given in actual dial specification.";
// Check sample input
if (parstring.find("fake_parameter") == std::string::npos)
return kGoodStatus;
// Parse inputs
std::vector strvct = GeneralUtils::ParseToStr(parstring, " ");
// Skip if comment or parameter somewhere later in line
if (strvct[0].c_str()[0] == '#' ||
strvct[0] == "fake_parameter"){
return kGoodStatus;
}
// Check length
if (strvct.size() < 3){
ERR(FTL) << "Fake dials need to provide at least 3 inputs." << std::endl;
std::cout << inputspec << std::endl;
return kErrorStatus;
}
// Read Inputs
std::string parname = strvct[1];
double parval = GeneralUtils::StrToDbl(strvct[2]);
// Setup Container
fFakeVals[parname] = parval;
// Print the fake parameter
LOG(MIN) << "Read Fake Parameter " << parname << " " << parval << std::endl;
// Tell reader its all good
return kGoodStatus;
}
//******************************************
int SystematicRoutines::ReadSamples(std::string samstring){
//******************************************
const static std::string inputspec =
"\tsample :inputfile.root [OPTS] "
"[norm]\nsample_name: Name "
"of sample to include. e.g. MiniBooNE_CCQE_XSec_1DQ2_nu\ninput_type: The "
"input event format. e.g. NEUT, GENIE, EVSPLN, ...\nOPTS: Additional, "
"optional sample options.\nnorm: Additional, optional sample "
"normalisation factor.";
// Check sample input
if (samstring.find("sample") == std::string::npos)
return kGoodStatus;
// Parse inputs
std::vector strvct = GeneralUtils::ParseToStr(samstring, " ");
// Skip if comment or parameter somewhere later in line
if (strvct[0].c_str()[0] == '#' ||
strvct[0] != "sample"){
return kGoodStatus;
}
// Check length
if (strvct.size() < 3){
ERR(FTL) << "Sample need to provide at least 3 inputs." << std::endl;
return kErrorStatus;
}
// Setup default inputs
std::string samname = strvct[1];
std::string samfile = strvct[2];
if (samfile == "FIX") {
ERR(FTL) << "Input filename was \"FIX\", this line is probably malformed "
"in the input card file. Line:\'"
<< samstring << "\'" << std::endl;
ERR(FTL) << "Expect sample lines to look like:\n\t" << inputspec
<< std::endl;
throw;
}
std::string samtype = "DEFAULT";
double samnorm = 1.0;
// Optional Type
if (strvct.size() > 3) {
samtype = strvct[3];
// samname += "_"+samtype;
// Also get rid of the / and replace it with underscore because it might not be supported character
// while (samname.find("/") != std::string::npos) {
// samname.replace(samname.find("/"), 1, std::string("_"));
// }
}
// Optional Norm
if (strvct.size() > 4) samnorm = GeneralUtils::StrToDbl(strvct[4]);
// Add Sample Names as Norm Dials
std::string normname = samname + "_norm";
// Check no repeat params
if (std::find(fParams.begin(), fParams.end(), normname) != fParams.end()){
ERR(FTL) << "Duplicate samples given for " << samname << std::endl;
throw;
}
fParams.push_back(normname);
fTypeVals[normname] = kNORM;
fStartVals[normname] = samnorm;
fCurVals[normname] = fStartVals[normname];
fErrorVals[normname] = 0.0;
fMinVals[normname] = 0.1;
fMaxVals[normname] = 10.0;
fStepVals[normname] = 0.5;
bool state = samtype.find("FREE") == std::string::npos;
fFixVals[normname] = state;
fStartFixVals[normname] = state;
// Print read in
LOG(MIN) << "Read sample " << samname << " "
<< samfile << " " << samtype << " "
<< samnorm << std::endl;
// Tell reader its all good
return kGoodStatus;
}
/*
Setup Functions
*/
//*************************************
void SystematicRoutines::SetupRWEngine(){
//*************************************
for (UInt_t i = 0; i < fParams.size(); i++){
std::string name = fParams[i];
FitBase::GetRW() -> IncludeDial(name, fTypeVals.at(name) );
}
UpdateRWEngine(fStartVals);
return;
}
//*************************************
void SystematicRoutines::SetupFCN(){
//*************************************
LOG(FIT)<<"Making the jointFCN"<Reconfigure();
fSampleFCN->ReconfigureAllEvents();
fSampleFCN->SetFakeData("MC");
UpdateRWEngine(fCurVals);
LOG(FIT)<<"Set all data to fake MC predictions."<SetFakeData(fFakeDataInput);
}
return;
}
//*****************************************
void SystematicRoutines::GetCovarFromFCN(){
//*****************************************
LOG(FIT) << "Loading ParamPull objects from FCN to build covar" << std::endl;
// Make helperstring
std::ostringstream helperstr;
// Keep track of what is being thrown
std::map dialthrowhandle;
// Get Covariance Objects from FCN
std::list inputpulls = fSampleFCN->GetPullList();
for (PullListConstIter iter = inputpulls.begin();
iter != inputpulls.end(); iter++){
ParamPull* pull = (*iter);
if (pull->GetType().find("THROW")){
fInputThrows.push_back(pull);
fInputCovar.push_back(pull->GetFullCovarMatrix());
fInputDials.push_back(pull->GetDataHist());
LOG(FIT) << "Read ParamPull: " << pull->GetName() << " " << pull->GetType() << std::endl;
}
TH1D dialhist = pull->GetDataHist();
TH1D minhist = pull->GetMinHist();
TH1D maxhist = pull->GetMaxHist();
TH1I typehist = pull->GetDialTypes();
for (int i = 0; i < dialhist.GetNbinsX(); i++){
std::string name = std::string(dialhist.GetXaxis()->GetBinLabel(i+1));
dialthrowhandle[name] = pull->GetName();
if (fCurVals.find(name) == fCurVals.end()){
// Add to Containers
fParams.push_back(name);
fCurVals[name] = dialhist.GetBinContent(i+1);
fStartVals[name] = dialhist.GetBinContent(i+1);
fMinVals[name] = minhist.GetBinContent(i+1);
fMaxVals[name] = maxhist.GetBinContent(i+1);
fStepVals[name] = 1.0;
fFixVals[name] = false;
fStartFixVals[name] = false;
fTypeVals[name] = typehist.GetBinContent(i+1);
fStateVals[name] = "FREE" + pull->GetType();
// Maker Helper
helperstr << std::string(16, ' ' ) << FitBase::ConvDialType(fTypeVals[name]) << " "
<< name << " " << fMinVals[name] << " "
<< fMaxVals[name] << " " << fStepVals[name] << " " << fStateVals[name]
<< std::endl;
}
}
}
// Check if no throws given
if (fInputThrows.empty()){
ERR(WRN) << "No covariances given to nuissyst" << std::endl;
ERR(WRN) << "Pushing back an uncorrelated gaussian throw error for each free parameter using step size" << std::endl;
for (UInt_t i = 0; i < fParams.size(); i++){
std::string syst = fParams[i];
if (fFixVals[syst]) continue;
// Make Terms
std::string name = syst + "_pull";
std::ostringstream pullterm;
pullterm << "DIAL:" << syst << ";"
<< fStartVals[syst] << ";"
<< fStepVals[syst];
std::string type = "GAUSTHROW/NEUT";
// Push Back Pulls
ParamPull* pull = new ParamPull( name, pullterm.str(), type );
fInputThrows.push_back(pull);
fInputCovar.push_back(pull->GetFullCovarMatrix());
fInputDials.push_back(pull->GetDataHist());
// Print Whats added
ERR(WRN) << "Added ParamPull : " << name << " " << pullterm.str() << " " << type << std::endl;
// Add helper string for future fits
helperstr << std::string(16, ' ' ) << "covar " << name << " " << pullterm.str() << " " << type << std::endl;
// Keep Track of Throws
dialthrowhandle[syst] = pull->GetName();
}
}
// Print Helper String
if (!helperstr.str().empty()){
LOG(FIT) << "To remove these statements in future studies, add the lines below to your card:" << std::endl;
// Can't use the logger properly because this can be multi-line. Use cout and added spaces to look better!
std::cout << helperstr.str();
sleep(2);
}
// Print Throw State
for (UInt_t i = 0; i < fParams.size(); i++){
std::string syst = fParams[i];
if (dialthrowhandle.find(syst) != dialthrowhandle.end()){
LOG(FIT) << "Dial " << i << ". " << setw(40) << syst << " = THROWING with " << dialthrowhandle[syst] << std::endl;
} else {
LOG(FIT) << "Dial " << i << ". " << setw(40) << syst << " = FIXED" << std::endl;
}
}
// Pause anyway
sleep(1);
return;
}
/*
Fitting Functions
*/
//*************************************
void SystematicRoutines::UpdateRWEngine(std::map& updateVals){
//*************************************
for (UInt_t i = 0; i < fParams.size(); i++){
std::string name = fParams[i];
if (updateVals.find(name) == updateVals.end()) continue;
FitBase::GetRW()->SetDialValue(name,updateVals.at(name));
}
FitBase::GetRW()->Reconfigure();
return;
}
//*************************************
void SystematicRoutines::PrintState(){
//*************************************
LOG(FIT)<<"------------"<GetLikelihood();
LOG(FIT) << std::left << std::setw(46) << "Likelihood for JointFCN: " << like << std::endl;
LOG(FIT)<<"------------"<cd();
SaveCurrentState();
}
//*************************************
void SystematicRoutines::SaveCurrentState(std::string subdir){
//*************************************
LOG(FIT)<<"Saving current full FCN predictions" <mkdir(subdir.c_str());
newdir->cd();
}
FitBase::GetRW()->Reconfigure();
fSampleFCN->ReconfigureAllEvents();
fSampleFCN->Write();
// Change back to current DIR
curdir->cd();
return;
}
//*************************************
void SystematicRoutines::SaveNominal(){
//*************************************
if (!fOutputRootFile)
fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
fOutputRootFile->cd();
LOG(FIT)<<"Saving Nominal Predictions (be cautious with this)" <Reconfigure();
SaveCurrentState("nominal");
};
//*************************************
void SystematicRoutines::SavePrefit(){
//*************************************
if (!fOutputRootFile)
fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
fOutputRootFile->cd();
LOG(FIT)<<"Saving Prefit Predictions"< 0){
fCovarFree = new TH2D("covariance_free",
"covariance_free",
NFREE,0,NFREE,
NFREE,0,NFREE);
}
// Set Bin Labels
int countall = 0;
int countfree = 0;
for (UInt_t i = 0; i < fParams.size(); i++){
fCovar->GetXaxis()->SetBinLabel(countall+1,fParams[i].c_str());
fCovar->GetYaxis()->SetBinLabel(countall+1,fParams[i].c_str());
countall++;
if (!fFixVals[fParams[i]] and NFREE > 0){
fCovarFree->GetXaxis()->SetBinLabel(countfree+1,fParams[i].c_str());
fCovarFree->GetYaxis()->SetBinLabel(countfree+1,fParams[i].c_str());
countfree++;
}
}
fCorrel = PlotUtils::GetCorrelationPlot(fCovar,"correlation");
fDecomp = PlotUtils::GetDecompPlot(fCovar,"decomposition");
if (NFREE > 0)fCorrelFree = PlotUtils::GetCorrelationPlot(fCovarFree, "correlation_free");
if (NFREE > 0)fDecompFree = PlotUtils::GetDecompPlot(fCovarFree,"decomposition_free");
return;
};
//*************************************
void SystematicRoutines::ThrowCovariance(bool uniformly){
//*************************************
// Set fThrownVals to all values in currentVals
for (UInt_t i = 0; i < fParams.size(); i++){
std::string name = fParams.at(i);
fThrownVals[name] = fCurVals[name];
}
for (PullListConstIter iter = fInputThrows.begin();
iter != fInputThrows.end(); iter++){
ParamPull* pull = *iter;
pull->ThrowCovariance();
TH1D dialhist = pull->GetDataHist();
for (int i = 0; i < dialhist.GetNbinsX(); i++){
std::string name = std::string(dialhist.GetXaxis()->GetBinLabel(i+1));
if (fCurVals.find(name) != fCurVals.end()){
fThrownVals[name] = dialhist.GetBinContent(i+1);
}
}
// Reset throw incase pulls are calculated.
pull->ResetToy();
}
return;
};
//*************************************
void SystematicRoutines::PlotLimits(){
//*************************************
std::cout << "Plotting Limits" << std::endl;
if (!fOutputRootFile)
fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
TDirectory* limfolder = (TDirectory*) fOutputRootFile->mkdir("Limits");
limfolder->cd();
// Set all parameters at their starting values
for (UInt_t i = 0; i < fParams.size(); i++){
fCurVals[fParams[i]] = fStartVals[fParams[i]];
}
TDirectory* nomfolder = (TDirectory*) limfolder->mkdir("nominal");
nomfolder->cd();
UpdateRWEngine(fCurVals);
fSampleFCN->ReconfigureAllEvents();
fSampleFCN->Write();
limfolder->cd();
std::vector allfolders;
// Loop through each parameter
for (UInt_t i = 0; i < fParams.size(); i++){
std::string syst = fParams[i];
std::cout << "Starting Param " << syst << std::endl;
if (fFixVals[syst]) continue;
// Loop Downwards
while (fCurVals[syst] > fMinVals[syst]){
fCurVals[syst] = fCurVals[syst] - fStepVals[syst];
// Check Limit
if (fCurVals[syst] < fMinVals[syst])
fCurVals[syst] = fMinVals[syst];
// Check folder exists
std::string curvalstring = std::string( Form( (syst + "_%f").c_str(), fCurVals[syst] ) );
if (std::find(allfolders.begin(), allfolders.end(), curvalstring) != allfolders.end())
break;
// Make new folder for variation
TDirectory* minfolder = (TDirectory*) limfolder->mkdir(Form( (syst + "_%f").c_str(), fCurVals[syst] ) );
minfolder->cd();
allfolders.push_back(curvalstring);
// Update Iterations
double *vals = FitUtils::GetArrayFromMap( fParams, fCurVals );
fSampleFCN->DoEval( vals );
delete vals;
// Save to folder
fSampleFCN->Write();
}
// Reset before next loop
fCurVals[syst] = fStartVals[syst];
// Loop Upwards now
while (fCurVals[syst] < fMaxVals[syst]){
fCurVals[syst] = fCurVals[syst] + fStepVals[syst];
// Check Limit
if (fCurVals[syst] > fMaxVals[syst])
fCurVals[syst] = fMaxVals[syst];
// Check folder exists
std::string curvalstring = std::string( Form( (syst + "_%f").c_str(), fCurVals[syst] ) );
if (std::find(allfolders.begin(), allfolders.end(), curvalstring) != allfolders.end())
break;
// Make new folder
TDirectory* maxfolder = (TDirectory*) limfolder->mkdir(Form( (syst + "_%f").c_str(), fCurVals[syst] ) );
maxfolder->cd();
allfolders.push_back(curvalstring);
// Update Iterations
double *vals = FitUtils::GetArrayFromMap( fParams, fCurVals );
fSampleFCN->DoEval( vals );
delete vals;
// Save to file
fSampleFCN->Write();
}
// Reset before leaving
fCurVals[syst] = fStartVals[syst];
UpdateRWEngine(fCurVals);
}
return;
}
//*************************************
void SystematicRoutines::Run(){
//*************************************
std::cout << "Running routines "<< std::endl;
fRoutines = GeneralUtils::ParseToStr(fStrategy,",");
for (UInt_t i = 0; i < fRoutines.size(); i++){
std::string routine = fRoutines.at(i);
int fitstate = kFitUnfinished;
LOG(FIT)<<"Running Routine: "<cd();
int nthrows = fNThrows;
int startthrows = fStartThrows;
int endthrows = startthrows + nthrows;
if (nthrows < 0) nthrows = endthrows;
if (startthrows < 0) startthrows = 0;
if (endthrows < 0) endthrows = startthrows + nthrows;
int seed = (gRandom->Uniform(0.0,1.0)*100000 + 100000000*(startthrows + endthrows) + time(NULL))/35;
gRandom->SetSeed(seed);
LOG(FIT) << "Using Seed : " << seed << std::endl;
LOG(FIT) << "nthrows = " << nthrows << std::endl;
LOG(FIT) << "startthrows = " << startthrows << std::endl;
LOG(FIT) << "endthrows = " << endthrows << std::endl;
UpdateRWEngine(fCurVals);
fSampleFCN->ReconfigureAllEvents();
if (startthrows == 0){
LOG(FIT) << "Making nominal " << std::endl;
TDirectory* nominal = (TDirectory*) tempfile->mkdir("nominal");
nominal->cd();
fSampleFCN->Write();
}
LOG(SAM) << "nthrows = " << nthrows << std::endl;
LOG(SAM) << "startthrows = " << startthrows << std::endl;
LOG(SAM) << "endthrows = " << endthrows << std::endl;
TTree* parameterTree = new TTree("throws","throws");
double chi2;
for (UInt_t i = 0; i < fParams.size(); i++)
parameterTree->Branch(fParams[i].c_str(), &fThrownVals[fParams[i]], (fParams[i] + "/D").c_str());
parameterTree->Branch("chi2",&chi2,"chi2/D");
fSampleFCN->CreateIterationTree("error_iterations", FitBase::GetRW());
// Would anybody actually want to do uniform throws of any parameter??
bool uniformly = FitPar::Config().GetParB("error_uniform");
// Run Throws and save
for (Int_t i = 0; i < endthrows+1; i++){
LOG(FIT) << "Loop " << i << std::endl;
ThrowCovariance(uniformly);
if (i < startthrows) continue;
if (i == 0) continue;
LOG(FIT) << "Throw " << i << " ================================" << std::endl;
// Generate Random Parameter Throw
// ThrowCovariance(uniformly);
TDirectory* throwfolder = (TDirectory*)tempfile->mkdir(Form("throw_%i",i));
throwfolder->cd();
// Run Eval
double *vals = FitUtils::GetArrayFromMap( fParams, fThrownVals );
chi2 = fSampleFCN->DoEval( vals );
delete vals;
// Save the FCN
fSampleFCN->Write();
parameterTree->Fill();
}
+ tempfile->cd();
fSampleFCN->WriteIterationTree();
tempfile->Close();
}
void SystematicRoutines::MergeThrows(){
fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
fOutputRootFile->cd();
// Make a container folder
TDirectory* errorDIR = (TDirectory*) fOutputRootFile->mkdir("error_bands");
errorDIR->cd();
TDirectory* outnominal = (TDirectory*) fOutputRootFile->mkdir("nominal_throw");
outnominal->cd();
// Split Input Files
if (!fThrowString.empty()) fThrowList = GeneralUtils::ParseToStr(fThrowString,",");
// Add default if no throwlist given
if (fThrowList.size() < 1) fThrowList.push_back( fOutputFile + ".throws.root" );
/// Save location of file containing nominal
std::string nominalfile;
bool nominalfound;
// Loop over files and check they exist.
for (uint i = 0; i < fThrowList.size(); i++){
std::string file = fThrowList[i];
bool found = false;
// normal
std::string newfile = file;
TFile* throwfile = new TFile(file.c_str(),"READ");
if (throwfile and !throwfile->IsZombie()){
found = true;
}
// normal.throws.root
if (!found){
newfile = file + ".throws.root";
throwfile = new TFile((file + ".throws.root").c_str(),"READ");
if (throwfile and !throwfile->IsZombie()) {
found = true;
}
}
// If its found save to throwlist, else save empty.
// Also search for nominal
if (found){
fThrowList[i] = newfile;
LOG(FIT) << "Throws File :" << newfile << std::endl;
// Find input which contains nominal
if (throwfile->Get("nominal")){
nominalfound = true;
nominalfile = newfile;
}
throwfile->Close();
} else {
fThrowList[i] = "";
}
delete throwfile;
}
// Make sure we have a nominal file
if (!nominalfound or nominalfile.empty()){
ERR(FTL) << "No nominal found when mergining! Exiting!" << std::endl;
throw;
}
// Get the nominal throws file
TFile* tempfile = new TFile((nominalfile).c_str(),"READ");
tempfile->cd();
TDirectory* nominal = (TDirectory*)tempfile->Get("nominal");
// int nthrows = FitPar::Config().GetParI("error_throws");
bool uniformly = FitPar::Config().GetParB("error_uniform");
// Check percentage of bad files is okay.
int badfilecount = 0;
for (uint i = 0; i < fThrowList.size(); i++){
if (!fThrowList[i].empty()){
LOG(FIT) << "Loading Throws From File " << i << " : "
<< fThrowList[i] << std::endl;
} else {
badfilecount++;
}
}
// Check we have at least one good file
if ((uint)badfilecount == fThrowList.size()){
ERR(FTL) << "Found no good throw files for MergeThrows" << std::endl;
throw;
} else if (badfilecount > fThrowList.size()*0.25){
ERR(WRN) << "Over 25% of your throw files are dodgy. Please check this is okay!" << std::endl;
ERR(WRN) << "Will continue for the time being..." << std::endl;
sleep(5);
}
// Now go through the keys in the temporary file and look for TH1D, and TH2D plots
TIter next(nominal->GetListOfKeys());
TKey *key;
while ((key = (TKey*)next())) {
TClass *cl = gROOT->GetClass(key->GetClassName());
if (!cl->InheritsFrom("TH1D") and !cl->InheritsFrom("TH2D")) continue;
TH1* baseplot = (TH1D*)key->ReadObj();
std::string plotname = std::string(baseplot->GetName());
LOG(FIT) << "Creating error bands for " << plotname;
if (LOG_LEVEL(FIT)){
if (!uniformly) std::cout << " : Using COVARIANCE Throws! " << std::endl;
else std::cout << " : Using UNIFORM THROWS!!! " << std::endl;
}
int nbins = 0;
if (cl->InheritsFrom("TH1D")) nbins = ((TH1D*)baseplot)->GetNbinsX();
else nbins = ((TH1D*)baseplot)->GetNbinsX()* ((TH1D*)baseplot)->GetNbinsY();
// Setup TProfile with RMS option
TProfile* tprof = new TProfile((plotname + "_prof").c_str(),(plotname + "_prof").c_str(),nbins, 0, nbins, "S");
// Setup The TTREE
double* bincontents;
bincontents = new double[nbins];
double* binlowest;
binlowest = new double[nbins];
double* binhighest;
binhighest = new double[nbins];
errorDIR->cd();
TTree* bintree = new TTree((plotname + "_tree").c_str(), (plotname + "_tree").c_str());
for (Int_t i = 0; i < nbins; i++){
bincontents[i] = 0.0;
binhighest[i] = 0.0;
binlowest[i] = 0.0;
bintree->Branch(Form("content_%i",i),&bincontents[i],Form("content_%i/D",i));
}
// Make new throw plot
TH1* newplot;
// Run Throw Merging.
for (UInt_t i = 0; i < fThrowList.size(); i++){
TFile* throwfile = new TFile(fThrowList[i].c_str(), "READ");
// Loop over all throws in a folder
TIter nextthrow(throwfile->GetListOfKeys());
TKey *throwkey;
while ((throwkey = (TKey*)nextthrow())) {
// Skip non throw folders
if (std::string(throwkey->GetName()).find("throw_") == std::string::npos) continue;
// Get Throw DIR
TDirectory* throwdir = (TDirectory*)throwkey->ReadObj();
// Get Plot From Throw
newplot = (TH1*)throwdir->Get(plotname.c_str());
if (!newplot) continue;
// Loop Over Plot
for (Int_t j = 0; j < nbins; j++){
tprof->Fill(j+0.5, newplot->GetBinContent(j+1));
bincontents[j] = newplot->GetBinContent(j+1);
if (bincontents[j] < binlowest[j] or i == 0) binlowest[j] = bincontents[j];
if (bincontents[j] > binhighest[j] or i == 0) binhighest[j] = bincontents[j];
}
errorDIR->cd();
bintree->Fill();
}
}
errorDIR->cd();
if (uniformly){
LOG(FIT) << "Uniformly Calculating Plot Errors!" << std::endl;
}
TH1* statplot = (TH1*) baseplot->Clone();
for (Int_t j = 0; j < nbins; j++){
if (!uniformly){
// if ((baseplot->GetBinError(j+1)/baseplot->GetBinContent(j+1)) < 1.0) {
// baseplot->SetBinError(j+1,sqrt(pow(tprof->GetBinError(j+1),2) + pow(baseplot->GetBinError(j+1),2)));
// } else {
baseplot->SetBinContent(j+1,tprof->GetBinContent(j+1));
baseplot->SetBinError(j+1,tprof->GetBinError(j+1));
// }
} else {
baseplot->SetBinContent(j+1, 0.0);//(binlowest[j] + binhighest[j]) / 2.0);
baseplot->SetBinError(j+1, 0.0); //(binhighest[j] - binlowest[j])/2.0);
}
}
errorDIR->cd();
baseplot->Write();
tprof->Write();
bintree->Write();
outnominal->cd();
for (int i = 0; i < nbins; i++){
baseplot->SetBinError(i+1, sqrt(pow(statplot->GetBinError(i+1),2) + pow(baseplot->GetBinError(i+1),2)));
}
baseplot->Write();
delete statplot;
delete baseplot;
delete tprof;
delete bintree;
delete [] bincontents;
}
return;
};