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; };