From d50dc1e431a7e7efda581610dd830f5413c04cec Mon Sep 17 00:00:00 2001 From: lgiannes Date: Wed, 30 Aug 2023 13:24:15 +0200 Subject: [PATCH 001/131] try --- src/Fitter/src/MinimizerBase.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Fitter/src/MinimizerBase.cpp b/src/Fitter/src/MinimizerBase.cpp index 8f24593f6..86775eca6 100644 --- a/src/Fitter/src/MinimizerBase.cpp +++ b/src/Fitter/src/MinimizerBase.cpp @@ -23,7 +23,7 @@ void MinimizerBase::readConfigImpl(){ bool showParametersOnFitMonitor = getLikelihood().getShowParametersOnFitMonitor(); showParametersOnFitMonitor = GenericToolbox::Json::fetchValue(_config_, "showParametersOnFitMonitor", showParametersOnFitMonitor); getLikelihood().setShowParametersOnFitMonitor(showParametersOnFitMonitor); -// +/// bool maxNbParametersPerLineOnMonitor = getLikelihood().getMaxNbParametersPerLineOnMonitor(); maxNbParametersPerLineOnMonitor = GenericToolbox::Json::fetchValue(_config_, "maxNbParametersPerLineOnMonitor", maxNbParametersPerLineOnMonitor); getLikelihood().setMaxNbParametersPerLineOnMonitor(maxNbParametersPerLineOnMonitor); From 1a01d00d626e8b29c902e965e63d954b56acd410 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Thu, 28 Sep 2023 16:54:45 +0200 Subject: [PATCH 002/131] added gundamMarginalise first version --- src/Applications/CMakeLists.txt | 4 +- src/Applications/src/gundamMarginalise.cxx | 275 +++++++++++++++++++++ 2 files changed, 278 insertions(+), 1 deletion(-) create mode 100644 src/Applications/src/gundamMarginalise.cxx diff --git a/src/Applications/CMakeLists.txt b/src/Applications/CMakeLists.txt index 5e7a2e29d..f74250138 100644 --- a/src/Applications/CMakeLists.txt +++ b/src/Applications/CMakeLists.txt @@ -10,7 +10,8 @@ set( APPLICATION_LIST gundamConfigCompare gundamPlotExtractor gundamConfigUnfolder - ) + gundamMarginalise +) if( ENABLE_DEV_MODE ) @@ -79,6 +80,7 @@ target_link_libraries( gundamConfigUnfolder GundamUtils ) target_link_libraries( gundamConfigCompare GundamUtils ) target_link_libraries( gundamPlotExtractor GundamUtils ) target_link_libraries( gundamCalcXsec GundamPropagator ) +target_link_libraries( gundamMarginalise GundamPropagator ) if( WITH_GUNDAM_ROOT_APP ) #target_sources( gundamRoot PRIVATE G__GundamRootDict.cxx ) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx new file mode 100644 index 000000000..e3d077d38 --- /dev/null +++ b/src/Applications/src/gundamMarginalise.cxx @@ -0,0 +1,275 @@ +// +// Created by Lorenzo Giannessi on 05.09.23. +// + +#include "GundamGlobals.h" +#include "GundamApp.h" +#include "GundamUtils.h" +#include "Propagator.h" +#include "ConfigUtils.h" + +#ifdef GUNDAM_USING_CACHE_MANAGER +#include "CacheManager.h" +#endif + +#include "Logger.h" +#include "CmdLineParser.h" +#include "GenericToolbox.h" +#include "GenericToolbox.Json.h" +#include "GenericToolbox.Root.h" + +#include +#include "TDirectory.h" +#include "TH1D.h" +#include "TH2D.h" + +#include +#include + + +LoggerInit([]{ + Logger::getUserHeader() << "[" << FILENAME << "]"; +}); + + +int main(int argc, char** argv){ + + GundamApp app{"contours marginalisation tool"}; + + // -------------------------- + // Read Command Line Args: + // -------------------------- + CmdLineParser clParser; + + clParser.addDummyOption("Main options:"); + + // I need a config file where the list of parameters to marginalise over are specified (look at the XSec config file to get inspired) + // TODO: understand the format/syntax of this file. This should be a required file. + clParser.addOption("configFile", {"-c"}, "Specify the parameters to marginalise over"); + + // (I think) I need the output file from a fitter to use as input here + clParser.addOption("fitterOutputFile", {"-f"}, "Specify the fitter output file"); + + // Think carefully what do you need to put in the output file + // 1. Marginalised covariance matrix + // what else? + clParser.addOption("outputFile", {"-o", "--out-file"}, "Specify the Marginaliser output file"); + + clParser.addOption("nbThreads", {"-t", "--nb-threads"}, "Specify nb of parallel threads"); + clParser.addOption("nToys", {"-n"}, "Specify number of toys"); + clParser.addOption("randomSeed", {"-s", "--seed"}, "Set random seed"); + + clParser.addDummyOption("Trigger options:"); + clParser.addTriggerOption("dryRun", {"-d", "--dry-run"}, "Only overrides fitter config and print it."); + + // I probably don't need this + //clParser.addTriggerOption("useBfAsXsec", {"--use-bf-as-xsec"}, "Use best-fit as x-sec value instead of mean of toys."); + + LogInfo << "Usage: " << std::endl; + LogInfo << clParser.getConfigSummary() << std::endl << std::endl; + + clParser.parseCmdLine(argc, argv); + + LogThrowIf(clParser.isNoOptionTriggered(), "No option was provided."); + + LogInfo << "Provided arguments: " << std::endl; + LogInfo << clParser.getValueSummary() << std::endl << std::endl; + + + // Sanity checks + LogThrowIf(not clParser.isOptionTriggered("configFile"), "Marginaliser config file not provided."); + LogThrowIf(not clParser.isOptionTriggered("fitterOutputFile"), "Did not provide the output fitter file."); + // Do I need to throw toys actually??? YES: throwing toys in your case is throwing according to the post-fit covariance matrix. + // You do it to compute the weights that eventually go in the computation of the marginalised likelihood. + LogThrowIf(not clParser.isOptionTriggered("nToys"), "Did not provide number of toys."); + + + // Global parameters + gRandom = new TRandom3(0); // Initialize with a UUID + if( clParser.isOptionTriggered("randomSeed") ){ + LogAlert << "Using user-specified random seed: " << clParser.getOptionVal("randomSeed") << std::endl; + gRandom->SetSeed(clParser.getOptionVal("randomSeed")); + } + else{ + ULong_t seed = time(nullptr); + LogInfo << "Using \"time(nullptr)\" random seed: " << seed << std::endl; + gRandom->SetSeed(seed); + } + + GundamGlobals::setNbThreads(clParser.getOptionVal("nbThreads", 1)); + LogInfo << "Running the fitter with " << GundamGlobals::getNbThreads() << " parallel threads." << std::endl; + + // Reading fitter file + LogInfo << "Opening fitter output file: " << clParser.getOptionVal("fitterOutputFile") << std::endl; + auto fitterFile = std::unique_ptr( TFile::Open( clParser.getOptionVal("fitterOutputFile").c_str() ) ); + LogThrowIf( fitterFile == nullptr, "Could not open fitter output file." ); + + using namespace GundamUtils; + ObjectReader::throwIfNotFound = true; + + nlohmann::json fitterConfig; + ObjectReader::readObject(fitterFile.get(), {{"gundam/config_TNamed"}, {"gundamFitter/unfoldedConfig_TNamed"}}, [&](TNamed* config_){ + fitterConfig = GenericToolbox::Json::readConfigJsonStr( config_->GetTitle() ); + }); + ConfigUtils::ConfigHandler cHandler{ fitterConfig }; + + // Disabling defined samples: + LogInfo << "Removing defined samples..." << std::endl; + ConfigUtils::applyOverrides( + cHandler.getConfig(), + GenericToolbox::Json::readConfigJsonStr(R"({"fitterEngineConfig":{"propagatorConfig":{"fitSampleSetConfig":{"fitSampleList":[]}}}})") + ); + + // Disabling defined plots: + LogInfo << "Removing defined plots..." << std::endl; + ConfigUtils::applyOverrides( + cHandler.getConfig(), + GenericToolbox::Json::readConfigJsonStr(R"({"fitterEngineConfig":{"propagatorConfig":{"plotGeneratorConfig":{}}}})") + ); + + // Defining signal samples + nlohmann::json margConfig{ ConfigUtils::readConfigFile( clParser.getOptionVal("configFile") ) }; + cHandler.override( margConfig ); + LogInfo << "Override done." << std::endl; + + if( clParser.isOptionTriggered("dryRun") ){ + std::cout << cHandler.toString() << std::endl; + + LogAlert << "Exiting as dry-run is set." << std::endl; + return EXIT_SUCCESS; + } + + auto configPropagator = GenericToolbox::Json::fetchValuePath( cHandler.getConfig(), "fitterEngineConfig/propagatorConfig" ); + + // Create a propagator object + Propagator propagator; + + // Read the whole fitter config with the override parameters + propagator.readConfig( configPropagator ); + + // We are only interested in our MC. Data has already been used to get the post-fit error/values + propagator.setLoadAsimovData( true ); + + // Disabling eigen decomposed parameters + propagator.setEnableEigenToOrigInPropagate( false ); + + // Load post-fit parameters as "prior" so we can reset the weight to this point when throwing toys + ObjectReader::readObject( fitterFile.get(), "FitterEngine/postFit/parState_TNamed", [&](TNamed* parState_){ + propagator.injectParameterValues( GenericToolbox::Json::readConfigJsonStr( parState_->GetTitle() ) ); + for( auto& parSet : propagator.getParameterSetsList() ){ + if( not parSet.isEnabled() ){ continue; } + for( auto& par : parSet.getParameterList() ){ + if( not par.isEnabled() ){ continue; } + par.setPriorValue( par.getParameterValue() ); + } + } + }); + + // Load the post-fit covariance matrix + ObjectReader::readObject( + fitterFile.get(), "FitterEngine/postFit/Hesse/hessian/postfitCovarianceOriginal_TH2D", + [&](TH2D* hCovPostFit_){ + propagator.setGlobalCovarianceMatrix(std::make_shared(hCovPostFit_->GetNbinsX(), hCovPostFit_->GetNbinsX())); + for( int iBin = 0 ; iBin < hCovPostFit_->GetNbinsX() ; iBin++ ){ + for( int jBin = 0 ; jBin < hCovPostFit_->GetNbinsX() ; jBin++ ){ + (*propagator.getGlobalCovarianceMatrix())[iBin][jBin] = hCovPostFit_->GetBinContent(1 + iBin, 1 + jBin); + } + } + }); + + + // Sample binning using parameterSetName + for( auto& sample : propagator.getFitSampleSet().getFitSampleList() ){ + auto associatedParSet = GenericToolbox::Json::fetchValue(sample.getConfig(), "parameterSetName"); + + // Looking for parSet + auto foundDialCollection = std::find_if( + propagator.getDialCollections().begin(), propagator.getDialCollections().end(), + [&](const DialCollection& dialCollection_){ + auto* parSetPtr{dialCollection_.getSupervisedParameterSet()}; + if( parSetPtr == nullptr ){ return false; } + return ( parSetPtr->getName() == associatedParSet ); + }); + LogThrowIf( + foundDialCollection == propagator.getDialCollections().end(), + "Could not find " << associatedParSet << " among fit dial collections: " + << GenericToolbox::iterableToString(propagator.getDialCollections(), [](const DialCollection& dialCollection_){ + return dialCollection_.getTitle(); + } + )); + + LogThrowIf(foundDialCollection->getDialBinSet().isEmpty(), "Could not find binning"); + sample.setBinningFilePath( foundDialCollection->getDialBinSet().getFilePath() ); + } + + // Load everything + propagator.initialize(); + + // Creating output file + std::string outFilePath{}; + if( clParser.isOptionTriggered("outputFile") ){ outFilePath = clParser.getOptionVal("outputFile"); } + else { + // appendixDict["optionName"] = "Appendix" + // this list insure all appendices will appear in the same order + std::vector> appendixDict{ + {"configFile", "%s"}, + {"fitterOutputFile", "Fit_%s"}, + {"nToys", "nToys_%s"}, + {"randomSeed", "Seed_%s"}, + }; + outFilePath = "marg_" + GundamUtils::generateFileName(clParser, appendixDict) + ".root"; + + std::string outFolder{GenericToolbox::Json::fetchValue(margConfig, "outputFolder", "./")}; + outFilePath = GenericToolbox::joinPath(outFolder, outFilePath); + } + + auto* marginalisationDir{ GenericToolbox::mkdirTFile(app.getOutfilePtr(), "marginalisation") }; + LogInfo << "Creating throws tree" << std::endl; + auto* margThrowTree = new TTree("margThrow", "margThrow"); + margThrowTree->SetDirectory( GenericToolbox::mkdirTFile(marginalisationDir, "throws") ); // temp saves will be done here + + + int nToys{ clParser.getOptionVal("nToys") }; + + + ////////////////////////////////////// + // THROWS LOOP + ///////////////////////////////////// + std::stringstream ss; ss << LogWarning.getPrefixString() << "Generating " << nToys << " toys..."; + for( int iToy = 0 ; iToy < nToys ; iToy++ ){ + + // loading... + GenericToolbox::displayProgressBar( iToy+1, nToys, ss.str() ); + + // Do the throwing: + // TODO: ask Adrien: when I do the throwing at this stage, the parameters to "be thrown" are already set? should I already set before the parameters to marginalise over? + propagator.throwParametersFromGlobalCovariance(); + //propagator.propagateParametersOnSamples(); + propagator.updateLlhCache(); + double LLH_value = propagator.getLlhBuffer(); + // get the prior covariance matrix for a subset of parameters (to marginalise over them) + + + // for now set it to false. Still need to understand/implement this + bool enableStatThrowInToys{false}; + if( enableStatThrowInToys ){ +// for( auto& xsec : crossSectionDataList ){ +// if( enableEventMcThrow ){ +// // Take into account the finite amount of event in MC +// xsec.samplePtr->getMcContainer().throwEventMcError(); +// } +// // Asimov bin content -> toy data +// xsec.samplePtr->getMcContainer().throwStatError(); +// } + } + + // i don't know what is this for, so I comment it for now +// writeBinDataFct(); + + // Write the branches + margThrowTree->Fill(); + } + + + GundamGlobals::getParallelWorker().reset(); +} From cc0425b625c986e902735a33183619dedcb3dc7d Mon Sep 17 00:00:00 2001 From: lgiannes Date: Thu, 28 Sep 2023 18:15:25 +0200 Subject: [PATCH 003/131] CMakeLists.txt modified for manual yaml lib finding --- CMakeLists.txt | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 09c45d129..0a1613d92 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -190,9 +190,13 @@ endif() # YAML-CPP cmessage( WARNING "Looking for YAML install..." ) -find_package( yaml-cpp REQUIRED HINTS ${YAMLCPP_DIR} ) -if(NOT yaml-cpp_FOUND) - cmessage(FATAL_ERROR "yaml-cpp library not found.") +if( "${YAML_CPP_INCLUDE_DIR} " STREQUAL " ") + find_package( yaml-cpp REQUIRED HINTS ${YAMLCPP_DIR} ) + if(NOT yaml-cpp_FOUND) + cmessage(FATAL_ERROR "yaml-cpp library not found.") + endif() +else() + add_definitions( -DCMDLINEPARSER_YAML_CPP_ENABLED=1 ) endif() include_directories( ${YAML_CPP_INCLUDE_DIR} ) cmessage( STATUS "Custom yaml-cpp installation: ${YAMLCPP_DIR}") From 0c2d5f566c9da9d5891f51b96e76d44cff29c2c2 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Mon, 2 Oct 2023 20:08:09 +0200 Subject: [PATCH 004/131] create TTree with parameters and LLH, prior, gLLH --- src/Applications/src/gundamMarginalise.cxx | 159 ++++++++++++++------- 1 file changed, 111 insertions(+), 48 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index e3d077d38..a97e64670 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -44,14 +44,13 @@ int main(int argc, char** argv){ clParser.addDummyOption("Main options:"); // I need a config file where the list of parameters to marginalise over are specified (look at the XSec config file to get inspired) - // TODO: understand the format/syntax of this file. This should be a required file. clParser.addOption("configFile", {"-c"}, "Specify the parameters to marginalise over"); // (I think) I need the output file from a fitter to use as input here clParser.addOption("fitterOutputFile", {"-f"}, "Specify the fitter output file"); // Think carefully what do you need to put in the output file - // 1. Marginalised covariance matrix + // 1. Marginalised covariance matrix: gotten from teh outoutFitter file // what else? clParser.addOption("outputFile", {"-o", "--out-file"}, "Specify the Marginaliser output file"); @@ -107,27 +106,31 @@ int main(int argc, char** argv){ using namespace GundamUtils; ObjectReader::throwIfNotFound = true; + + nlohmann::json fitterConfig; ObjectReader::readObject(fitterFile.get(), {{"gundam/config_TNamed"}, {"gundamFitter/unfoldedConfig_TNamed"}}, [&](TNamed* config_){ fitterConfig = GenericToolbox::Json::readConfigJsonStr( config_->GetTitle() ); }); ConfigUtils::ConfigHandler cHandler{ fitterConfig }; - // Disabling defined samples: - LogInfo << "Removing defined samples..." << std::endl; - ConfigUtils::applyOverrides( - cHandler.getConfig(), - GenericToolbox::Json::readConfigJsonStr(R"({"fitterEngineConfig":{"propagatorConfig":{"fitSampleSetConfig":{"fitSampleList":[]}}}})") - ); - - // Disabling defined plots: - LogInfo << "Removing defined plots..." << std::endl; - ConfigUtils::applyOverrides( - cHandler.getConfig(), - GenericToolbox::Json::readConfigJsonStr(R"({"fitterEngineConfig":{"propagatorConfig":{"plotGeneratorConfig":{}}}})") - ); - - // Defining signal samples +// // Disabling defined samples: +// LogInfo << "Removing defined samples..." << std::endl; +// ConfigUtils::applyOverrides( +// cHandler.getConfig(), +// GenericToolbox::Json::readConfigJsonStr(R"({"fitterEngineConfig":{"propagatorConfig":{"fitSampleSetConfig":{"fitSampleList":[]}}}})") +// ); + +// // Disabling defined plots: +// LogInfo << "Removing defined plots..." << std::endl; +// ConfigUtils::applyOverrides( +// cHandler.getConfig(), +// GenericToolbox::Json::readConfigJsonStr(R"({"fitterEngineConfig":{"propagatorConfig":{"plotGeneratorConfig":{}}}})") +// ); + + + + // Reading marginaliser config file nlohmann::json margConfig{ ConfigUtils::readConfigFile( clParser.getOptionVal("configFile") ) }; cHandler.override( margConfig ); LogInfo << "Override done." << std::endl; @@ -144,7 +147,7 @@ int main(int argc, char** argv){ // Create a propagator object Propagator propagator; - // Read the whole fitter config with the override parameters +// // Read the whole fitter config with the override parameters propagator.readConfig( configPropagator ); // We are only interested in our MC. Data has already been used to get the post-fit error/values @@ -153,13 +156,31 @@ int main(int argc, char** argv){ // Disabling eigen decomposed parameters propagator.setEnableEigenToOrigInPropagate( false ); + // Load everything + propagator.initialize(); + + propagator.getParameterSetsList(); + for( auto& parSet : propagator.getParameterSetsList() ){ + if( not parSet.isEnabled() ){ continue; } + LogInfo < "<( fitterFile.get(), "FitterEngine/postFit/parState_TNamed", [&](TNamed* parState_){ propagator.injectParameterValues( GenericToolbox::Json::readConfigJsonStr( parState_->GetTitle() ) ); for( auto& parSet : propagator.getParameterSetsList() ){ if( not parSet.isEnabled() ){ continue; } + LogInfo<< parSet.getName()< "<(sample.getConfig(), "parameterSetName"); +// +// // Looking for parSet +// auto foundDialCollection = std::find_if( +// propagator.getDialCollections().begin(), propagator.getDialCollections().end(), +// [&](const DialCollection& dialCollection_){ +// auto* parSetPtr{dialCollection_.getSupervisedParameterSet()}; +// if( parSetPtr == nullptr ){ return false; } +// return ( parSetPtr->getName() == associatedParSet ); +// }); +// LogThrowIf( +// foundDialCollection == propagator.getDialCollections().end(), +// "Could not find " << associatedParSet << " among fit dial collections: " +// << GenericToolbox::iterableToString(propagator.getDialCollections(), [](const DialCollection& dialCollection_){ +// return dialCollection_.getTitle(); +// } +// )); +// +// LogThrowIf(foundDialCollection->getDialBinSet().isEmpty(), "Could not find binning"); +// sample.setBinningFilePath( foundDialCollection->getDialBinSet().getFilePath() ); +// } - // Sample binning using parameterSetName - for( auto& sample : propagator.getFitSampleSet().getFitSampleList() ){ - auto associatedParSet = GenericToolbox::Json::fetchValue(sample.getConfig(), "parameterSetName"); - - // Looking for parSet - auto foundDialCollection = std::find_if( - propagator.getDialCollections().begin(), propagator.getDialCollections().end(), - [&](const DialCollection& dialCollection_){ - auto* parSetPtr{dialCollection_.getSupervisedParameterSet()}; - if( parSetPtr == nullptr ){ return false; } - return ( parSetPtr->getName() == associatedParSet ); - }); - LogThrowIf( - foundDialCollection == propagator.getDialCollections().end(), - "Could not find " << associatedParSet << " among fit dial collections: " - << GenericToolbox::iterableToString(propagator.getDialCollections(), [](const DialCollection& dialCollection_){ - return dialCollection_.getTitle(); - } - )); - - LogThrowIf(foundDialCollection->getDialBinSet().isEmpty(), "Could not find binning"); - sample.setBinningFilePath( foundDialCollection->getDialBinSet().getFilePath() ); - } - - // Load everything - propagator.initialize(); // Creating output file std::string outFilePath{}; @@ -223,10 +241,32 @@ int main(int argc, char** argv){ outFilePath = GenericToolbox::joinPath(outFolder, outFilePath); } + app.setCmdLinePtr( &clParser ); + app.setConfigString( ConfigUtils::ConfigHandler{margConfig}.toString() ); + app.openOutputFile( outFilePath ); + app.writeAppInfo(); + auto* marginalisationDir{ GenericToolbox::mkdirTFile(app.getOutfilePtr(), "marginalisation") }; LogInfo << "Creating throws tree" << std::endl; auto* margThrowTree = new TTree("margThrow", "margThrow"); margThrowTree->SetDirectory( GenericToolbox::mkdirTFile(marginalisationDir, "throws") ); // temp saves will be done here + // make a TTree with the following branches, to be filled with the throws + // 1. marginalised parameters drew: vector + // 2. non-marginalised parameters drew: vector + // 3. LLH value: double + // 4. "g" value: the chi square value as extracted from the covariance matrix: double + // 5. value of the priors for the marginalised parameters: vector + + std::vector eta; + std::vector theta; + std::vector prior; + double LLH, gLLH; + margThrowTree->Branch("marginalisedParameters", &eta); + margThrowTree->Branch("nonMarginalisedParameters", &theta); + margThrowTree->Branch("prior", &prior); + margThrowTree->Branch("LLH", &LLH); + margThrowTree->Branch("gLLH", &gLLH); + int nToys{ clParser.getOptionVal("nToys") }; @@ -236,18 +276,39 @@ int main(int argc, char** argv){ // THROWS LOOP ///////////////////////////////////// std::stringstream ss; ss << LogWarning.getPrefixString() << "Generating " << nToys << " toys..."; + + for( int iToy = 0 ; iToy < nToys ; iToy++ ){ // loading... GenericToolbox::displayProgressBar( iToy+1, nToys, ss.str() ); // Do the throwing: - // TODO: ask Adrien: when I do the throwing at this stage, the parameters to "be thrown" are already set? should I already set before the parameters to marginalise over? propagator.throwParametersFromGlobalCovariance(); - //propagator.propagateParametersOnSamples(); + //propagator.propagateParametersOnSamples(); // Adrien says this is not necessary propagator.updateLlhCache(); - double LLH_value = propagator.getLlhBuffer(); - // get the prior covariance matrix for a subset of parameters (to marginalise over them) + LLH = propagator.getLlhBuffer(); + + + + for( auto& parSet : propagator.getParameterSetsList() ){ + if( not parSet.isEnabled() ){ continue; } +// LogInfo<< parSet.getName()< "<getParameterList().at(0).getParameterValue(); +// +// } // for now set it to false. Still need to understand/implement this @@ -270,6 +331,8 @@ int main(int argc, char** argv){ margThrowTree->Fill(); } + margThrowTree->Write(); + GundamGlobals::getParallelWorker().reset(); } From 465d02dfb706346019e138a252b1b7f36431dd0c Mon Sep 17 00:00:00 2001 From: lgiannes Date: Tue, 3 Oct 2023 11:57:54 +0200 Subject: [PATCH 005/131] updated TTree structure --- src/Applications/src/gundamMarginalise.cxx | 24 ++++++++++++---------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index a97e64670..56ac57ae5 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -177,10 +177,10 @@ int main(int argc, char** argv){ propagator.injectParameterValues( GenericToolbox::Json::readConfigJsonStr( parState_->GetTitle() ) ); for( auto& parSet : propagator.getParameterSetsList() ){ if( not parSet.isEnabled() ){ continue; } - LogInfo<< parSet.getName()< "< "< - std::vector eta; - std::vector theta; + std::vector parameters; + std::vector margThis; std::vector prior; double LLH, gLLH; - margThrowTree->Branch("marginalisedParameters", &eta); - margThrowTree->Branch("nonMarginalisedParameters", &theta); + margThrowTree->Branch("Parameters", ¶meters); + margThrowTree->Branch("Marginalise", &margThis); margThrowTree->Branch("prior", &prior); margThrowTree->Branch("LLH", &LLH); margThrowTree->Branch("gLLH", &gLLH); @@ -288,18 +288,20 @@ int main(int argc, char** argv){ //propagator.propagateParametersOnSamples(); // Adrien says this is not necessary propagator.updateLlhCache(); LLH = propagator.getLlhBuffer(); + gLLH = 0; - - + parameters.clear(); + margThis.clear(); + prior.clear(); for( auto& parSet : propagator.getParameterSetsList() ){ if( not parSet.isEnabled() ){ continue; } // LogInfo<< parSet.getName()< "< Date: Tue, 31 Oct 2023 10:53:44 +0100 Subject: [PATCH 006/131] added a function in Propagator.cpp to save the weights while parameters are thrown --- src/Applications/src/gundamMarginalise.cxx | 66 ++++++++++--- src/Propagator/include/Propagator.h | 3 +- src/Propagator/src/Propagator.cpp | 105 +++++++++++++++++++++ 3 files changed, 161 insertions(+), 13 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 56ac57ae5..d8dd1fc8c 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -57,6 +57,8 @@ int main(int argc, char** argv){ clParser.addOption("nbThreads", {"-t", "--nb-threads"}, "Specify nb of parallel threads"); clParser.addOption("nToys", {"-n"}, "Specify number of toys"); clParser.addOption("randomSeed", {"-s", "--seed"}, "Set random seed"); + clParser.addTriggerOption("usingGpu", {"--gpu"}, "Use GPU parallelization"); + clParser.addDummyOption("Trigger options:"); clParser.addTriggerOption("dryRun", {"-d", "--dry-run"}, "Only overrides fitter config and print it."); @@ -135,6 +137,11 @@ int main(int argc, char** argv){ cHandler.override( margConfig ); LogInfo << "Override done." << std::endl; + // read the parameters to include in the TTree + + // read the parameters to marginalise over + + if( clParser.isOptionTriggered("dryRun") ){ std::cout << cHandler.toString() << std::endl; @@ -185,6 +192,12 @@ int main(int argc, char** argv){ } } }); + // also save the value of the LLH at the best fit point: + propagator.propagateParametersOnSamples(); + propagator.updateLlhCache(); + double bestFitLLH = propagator.getLlhBuffer(); + LogInfo<<"Best fit LLH: "<( @@ -260,51 +273,80 @@ int main(int argc, char** argv){ std::vector parameters; std::vector margThis; std::vector prior; - double LLH, gLLH; + std::vector weightsChiSquare; +// weightsChiSquare.reserve(1000); + double LLH, gLLH, priorSum, LLHwrtBestFit; margThrowTree->Branch("Parameters", ¶meters); margThrowTree->Branch("Marginalise", &margThis); margThrowTree->Branch("prior", &prior); margThrowTree->Branch("LLH", &LLH); + margThrowTree->Branch("LLHwrtBestFit", &LLHwrtBestFit); margThrowTree->Branch("gLLH", &gLLH); - + margThrowTree->Branch("priorSum", &priorSum); + margThrowTree->Branch("weightsChiSquare", &weightsChiSquare); int nToys{ clParser.getOptionVal("nToys") }; + // Get parameters to be marginalised + //std::string marginalisedParameters{GenericToolbox::Json::fetchValue(margConfig, "marginalisedParameters", "")}; + + + ////////////////////////////////////// // THROWS LOOP ///////////////////////////////////// std::stringstream ss; ss << LogWarning.getPrefixString() << "Generating " << nToys << " toys..."; + LogInfo<<"Prior information: "< type: "< "< gLLH: "< toy data // xsec.samplePtr->getMcContainer().throwStatError(); // } - } +// } // i don't know what is this for, so I comment it for now // writeBinDataFct(); diff --git a/src/Propagator/include/Propagator.h b/src/Propagator/include/Propagator.h index 6fcb44d39..8a8665af0 100644 --- a/src/Propagator/include/Propagator.h +++ b/src/Propagator/include/Propagator.h @@ -80,8 +80,9 @@ class Propagator : public JsonBaseClass { [[nodiscard]] nlohmann::json exportParameterInjectorConfig() const; void injectParameterValues(const nlohmann::json &config_); void throwParametersFromGlobalCovariance(); + void throwParametersFromGlobalCovariance(std::vector &weightsChiSquare); - // Logger related + // Logger related static void muteLogger(); static void unmuteLogger(); diff --git a/src/Propagator/src/Propagator.cpp b/src/Propagator/src/Propagator.cpp index c4598b39c..75c8a0892 100644 --- a/src/Propagator/src/Propagator.cpp +++ b/src/Propagator/src/Propagator.cpp @@ -867,7 +867,112 @@ void Propagator::throwParametersFromGlobalCovariance(){ } +void Propagator::throwParametersFromGlobalCovariance(std::vector &weightsChiSquare){ + // check that weightsChiSquare is an empty vector + LogThrowIf( not weightsChiSquare.empty(), "ERROR: argument weightsChiSquare is not empty" ); + + if( _strippedCovarianceMatrix_ == nullptr ){ + LogInfo << "Creating stripped global covariance matrix..." << std::endl; + LogThrowIf( _globalCovarianceMatrix_ == nullptr, "Global covariance matrix not set." ); + int nStripped{0}; + for( int iDiag = 0 ; iDiag < _globalCovarianceMatrix_->GetNrows() ; iDiag++ ){ + if( (*_globalCovarianceMatrix_)[iDiag][iDiag] != 0 ){ nStripped++; } + } + + LogInfo << "Stripped global covariance matrix is " << nStripped << "x" << nStripped << std::endl; + _strippedCovarianceMatrix_ = std::make_shared(nStripped, nStripped); + int iStrippedBin{-1}; + for( int iBin = 0 ; iBin < _globalCovarianceMatrix_->GetNrows() ; iBin++ ){ + if( (*_globalCovarianceMatrix_)[iBin][iBin] == 0 ){ continue; } + iStrippedBin++; + int jStrippedBin{-1}; + for( int jBin = 0 ; jBin < _globalCovarianceMatrix_->GetNrows() ; jBin++ ){ + if( (*_globalCovarianceMatrix_)[jBin][jBin] == 0 ){ continue; } + jStrippedBin++; + (*_strippedCovarianceMatrix_)[iStrippedBin][jStrippedBin] = (*_globalCovarianceMatrix_)[iBin][jBin]; + } + } + + _strippedParameterList_.reserve( nStripped ); + for( auto& parSet : _parameterSetList_ ){ + if( not parSet.isEnabled() ) continue; + for( auto& par : parSet.getParameterList() ){ + if( not par.isEnabled() ) continue; + _strippedParameterList_.emplace_back(&par); + } + } + LogThrowIf( _strippedParameterList_.size() != nStripped, "Enabled parameters list don't correspond to the matrix" ); + } + + if( _choleskyMatrix_ == nullptr ){ + LogInfo << "Generating global cholesky matrix" << std::endl; + _choleskyMatrix_ = std::shared_ptr( + GenericToolbox::getCholeskyMatrix(_strippedCovarianceMatrix_.get()) + ); + } + + bool keepThrowing{true}; +// int throwNb{0}; + + while( keepThrowing ){ +// throwNb++; + bool rethrow{false}; + std::vector throws,weights; + GenericToolbox::throwCorrelatedParameters(_choleskyMatrix_.get(),throws, weights); + if(throws.size() != weights.size()){ + LogInfo<<"WARNING: throws.size() != weights.size() "<< throws.size()<GetNrows() ; iPar++ ){ + _strippedParameterList_[iPar]->setParameterValue( + _strippedParameterList_[iPar]->getPriorValue() + + throws[iPar] + ); + weightsChiSquare.push_back(weights[iPar]); + + if( _reThrowParSetIfOutOfBounds_ ){ + if( not _strippedParameterList_[iPar]->isValueWithinBounds() ){ + // re-do the throwing +// LogDebug << "Not within bounds: " << _strippedParameterList_[iPar]->getSummary() << std::endl; + rethrow = true; + } + } + } + + // Making sure eigen decomposed parameters get the conversion done + for( auto& parSet : _parameterSetList_ ){ + if( not parSet.isEnabled() ) continue; + if( parSet.isUseEigenDecompInFit() ){ + parSet.propagateOriginalToEigen(); + + // also check the bounds of real parameter space + if( _reThrowParSetIfOutOfBounds_ ){ + for( auto& par : parSet.getEigenParameterList() ){ + if( not par.isEnabled() ) continue; + if( not par.isValueWithinBounds() ){ + // re-do the throwing + rethrow = true; + break; + } + } + } + } + } + + + if( rethrow ){ + // wrap back to the while loop +// LogDebug << "RE-THROW #" << throwNb << std::endl; + continue; + } + + // reached this point: all parameters are within bounds + keepThrowing = false; + } + + + +} // Protected void Propagator::initializeThreads() { reweightMcEventsFct = [this](int iThread){ From ef0c45132b58cd2d5f5740aecda67e3fa73513bb Mon Sep 17 00:00:00 2001 From: lgiannes Date: Tue, 31 Oct 2023 11:52:53 +0100 Subject: [PATCH 007/131] moved function throwParametersFromGlobalCovariance(weightschiSquare) from Propagator to ParametersManager --- src/Applications/src/gundamMarginalise.cxx | 18 +-- .../include/ParametersManager.h | 1 + .../src/ParametersManager.cpp | 107 ++++++++++++++++++ src/Propagator/src/Propagator.cpp | 106 ----------------- 4 files changed, 117 insertions(+), 115 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index d8dd1fc8c..ede6aa5ef 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -166,8 +166,8 @@ int main(int argc, char** argv){ // Load everything propagator.initialize(); - propagator.getParameterSetsList(); - for( auto& parSet : propagator.getParameterSetsList() ){ + propagator.getParametersManager().getParameterSetsList(); + for( auto& parSet : propagator.getParametersManager().getParameterSetsList() ){ if( not parSet.isEnabled() ){ continue; } LogInfo <( fitterFile.get(), "FitterEngine/postFit/parState_TNamed", [&](TNamed* parState_){ - propagator.injectParameterValues( GenericToolbox::Json::readConfigJsonStr( parState_->GetTitle() ) ); - for( auto& parSet : propagator.getParameterSetsList() ){ + propagator.getParametersManager().injectParameterValues( GenericToolbox::Json::readConfigJsonStr( parState_->GetTitle() ) ); + for( auto& parSet : propagator.getParametersManager().getParameterSetsList() ){ if( not parSet.isEnabled() ){ continue; } // LogInfo<< parSet.getName()<( fitterFile.get(), "FitterEngine/postFit/Hesse/hessian/postfitCovarianceOriginal_TH2D", [&](TH2D* hCovPostFit_){ - propagator.setGlobalCovarianceMatrix(std::make_shared(hCovPostFit_->GetNbinsX(), hCovPostFit_->GetNbinsX())); + propagator.getParametersManager().setGlobalCovarianceMatrix(std::make_shared(hCovPostFit_->GetNbinsX(), hCovPostFit_->GetNbinsX())); for( int iBin = 0 ; iBin < hCovPostFit_->GetNbinsX() ; iBin++ ){ for( int jBin = 0 ; jBin < hCovPostFit_->GetNbinsX() ; jBin++ ){ - (*propagator.getGlobalCovarianceMatrix())[iBin][jBin] = hCovPostFit_->GetBinContent(1 + iBin, 1 + jBin); + (*propagator.getParametersManager().getGlobalCovarianceMatrix())[iBin][jBin] = hCovPostFit_->GetBinContent(1 + iBin, 1 + jBin); } } }); @@ -300,7 +300,7 @@ int main(int argc, char** argv){ std::stringstream ss; ss << LogWarning.getPrefixString() << "Generating " << nToys << " toys..."; LogInfo<<"Prior information: "< &weightsChiSquare); ParameterSet* getFitParameterSetPtr(const std::string& name_); // Logger related diff --git a/src/ParametersManager/src/ParametersManager.cpp b/src/ParametersManager/src/ParametersManager.cpp index 509efc5b1..8af4f28a5 100644 --- a/src/ParametersManager/src/ParametersManager.cpp +++ b/src/ParametersManager/src/ParametersManager.cpp @@ -277,6 +277,113 @@ void ParametersManager::throwParametersFromGlobalCovariance(bool quietVerbose_){ // reached this point: all parameters are within bounds keepThrowing = false; } +} + +void ParametersManager::throwParametersFromGlobalCovariance(std::vector &weightsChiSquare){ + + // check that weightsChiSquare is an empty vector + LogThrowIf( not weightsChiSquare.empty(), "ERROR: argument weightsChiSquare is not empty" ); + + if( _strippedCovarianceMatrix_ == nullptr ){ + LogInfo << "Creating stripped global covariance matrix..." << std::endl; + LogThrowIf( _globalCovarianceMatrix_ == nullptr, "Global covariance matrix not set." ); + int nStripped{0}; + for( int iDiag = 0 ; iDiag < _globalCovarianceMatrix_->GetNrows() ; iDiag++ ){ + if( (*_globalCovarianceMatrix_)[iDiag][iDiag] != 0 ){ nStripped++; } + } + + LogInfo << "Stripped global covariance matrix is " << nStripped << "x" << nStripped << std::endl; + _strippedCovarianceMatrix_ = std::make_shared(nStripped, nStripped); + int iStrippedBin{-1}; + for( int iBin = 0 ; iBin < _globalCovarianceMatrix_->GetNrows() ; iBin++ ){ + if( (*_globalCovarianceMatrix_)[iBin][iBin] == 0 ){ continue; } + iStrippedBin++; + int jStrippedBin{-1}; + for( int jBin = 0 ; jBin < _globalCovarianceMatrix_->GetNrows() ; jBin++ ){ + if( (*_globalCovarianceMatrix_)[jBin][jBin] == 0 ){ continue; } + jStrippedBin++; + (*_strippedCovarianceMatrix_)[iStrippedBin][jStrippedBin] = (*_globalCovarianceMatrix_)[iBin][jBin]; + } + } + + _strippedParameterList_.reserve( nStripped ); + for( auto& parSet : _parameterSetList_ ){ + if( not parSet.isEnabled() ) continue; + for( auto& par : parSet.getParameterList() ){ + if( not par.isEnabled() ) continue; + _strippedParameterList_.emplace_back(&par); + } + } + LogThrowIf( _strippedParameterList_.size() != nStripped, "Enabled parameters list don't correspond to the matrix" ); + } + + if( _choleskyMatrix_ == nullptr ){ + LogInfo << "Generating global cholesky matrix" << std::endl; + _choleskyMatrix_ = std::shared_ptr( + GenericToolbox::getCholeskyMatrix(_strippedCovarianceMatrix_.get()) + ); + } + + bool keepThrowing{true}; +// int throwNb{0}; + + while( keepThrowing ){ +// throwNb++; + bool rethrow{false}; + std::vector throws,weights; + GenericToolbox::throwCorrelatedParameters(_choleskyMatrix_.get(),throws, weights); + if(throws.size() != weights.size()){ + LogInfo<<"WARNING: throws.size() != weights.size() "<< throws.size()<GetNrows() ; iPar++ ){ + _strippedParameterList_[iPar]->setParameterValue( + _strippedParameterList_[iPar]->getPriorValue() + + throws[iPar] + ); + weightsChiSquare.push_back(weights[iPar]); + + if( _reThrowParSetIfOutOfBounds_ ){ + if( not _strippedParameterList_[iPar]->isValueWithinBounds() ){ + // re-do the throwing +// LogDebug << "Not within bounds: " << _strippedParameterList_[iPar]->getSummary() << std::endl; + rethrow = true; + } + } + } + + // Making sure eigen decomposed parameters get the conversion done + for( auto& parSet : _parameterSetList_ ){ + if( not parSet.isEnabled() ) continue; + if( parSet.isUseEigenDecompInFit() ){ + parSet.propagateOriginalToEigen(); + + // also check the bounds of real parameter space + if( _reThrowParSetIfOutOfBounds_ ){ + for( auto& par : parSet.getEigenParameterList() ){ + if( not par.isEnabled() ) continue; + if( not par.isValueWithinBounds() ){ + // re-do the throwing + rethrow = true; + break; + } + } + } + } + } + + + if( rethrow ){ + // wrap back to the while loop +// LogDebug << "RE-THROW #" << throwNb << std::endl; + continue; + } + + // reached this point: all parameters are within bounds + keepThrowing = false; + } + + + } void ParametersManager::injectParameterValues(const nlohmann::json &config_) { LogWarning << "Injecting parameters..." << std::endl; diff --git a/src/Propagator/src/Propagator.cpp b/src/Propagator/src/Propagator.cpp index 5fa816727..518d56bf9 100644 --- a/src/Propagator/src/Propagator.cpp +++ b/src/Propagator/src/Propagator.cpp @@ -614,112 +614,6 @@ std::string Propagator::getSampleBreakdownTableStr() const{ std::stringstream ss; ss << t.generateTableString(); return ss.str(); -} -void Propagator::throwParametersFromGlobalCovariance(std::vector &weightsChiSquare){ - - // check that weightsChiSquare is an empty vector - LogThrowIf( not weightsChiSquare.empty(), "ERROR: argument weightsChiSquare is not empty" ); - - if( _strippedCovarianceMatrix_ == nullptr ){ - LogInfo << "Creating stripped global covariance matrix..." << std::endl; - LogThrowIf( _globalCovarianceMatrix_ == nullptr, "Global covariance matrix not set." ); - int nStripped{0}; - for( int iDiag = 0 ; iDiag < _globalCovarianceMatrix_->GetNrows() ; iDiag++ ){ - if( (*_globalCovarianceMatrix_)[iDiag][iDiag] != 0 ){ nStripped++; } - } - - LogInfo << "Stripped global covariance matrix is " << nStripped << "x" << nStripped << std::endl; - _strippedCovarianceMatrix_ = std::make_shared(nStripped, nStripped); - int iStrippedBin{-1}; - for( int iBin = 0 ; iBin < _globalCovarianceMatrix_->GetNrows() ; iBin++ ){ - if( (*_globalCovarianceMatrix_)[iBin][iBin] == 0 ){ continue; } - iStrippedBin++; - int jStrippedBin{-1}; - for( int jBin = 0 ; jBin < _globalCovarianceMatrix_->GetNrows() ; jBin++ ){ - if( (*_globalCovarianceMatrix_)[jBin][jBin] == 0 ){ continue; } - jStrippedBin++; - (*_strippedCovarianceMatrix_)[iStrippedBin][jStrippedBin] = (*_globalCovarianceMatrix_)[iBin][jBin]; - } - } - - _strippedParameterList_.reserve( nStripped ); - for( auto& parSet : _parameterSetList_ ){ - if( not parSet.isEnabled() ) continue; - for( auto& par : parSet.getParameterList() ){ - if( not par.isEnabled() ) continue; - _strippedParameterList_.emplace_back(&par); - } - } - LogThrowIf( _strippedParameterList_.size() != nStripped, "Enabled parameters list don't correspond to the matrix" ); - } - - if( _choleskyMatrix_ == nullptr ){ - LogInfo << "Generating global cholesky matrix" << std::endl; - _choleskyMatrix_ = std::shared_ptr( - GenericToolbox::getCholeskyMatrix(_strippedCovarianceMatrix_.get()) - ); - } - - bool keepThrowing{true}; -// int throwNb{0}; - - while( keepThrowing ){ -// throwNb++; - bool rethrow{false}; - std::vector throws,weights; - GenericToolbox::throwCorrelatedParameters(_choleskyMatrix_.get(),throws, weights); - if(throws.size() != weights.size()){ - LogInfo<<"WARNING: throws.size() != weights.size() "<< throws.size()<GetNrows() ; iPar++ ){ - _strippedParameterList_[iPar]->setParameterValue( - _strippedParameterList_[iPar]->getPriorValue() - + throws[iPar] - ); - weightsChiSquare.push_back(weights[iPar]); - - if( _reThrowParSetIfOutOfBounds_ ){ - if( not _strippedParameterList_[iPar]->isValueWithinBounds() ){ - // re-do the throwing -// LogDebug << "Not within bounds: " << _strippedParameterList_[iPar]->getSummary() << std::endl; - rethrow = true; - } - } - } - - // Making sure eigen decomposed parameters get the conversion done - for( auto& parSet : _parameterSetList_ ){ - if( not parSet.isEnabled() ) continue; - if( parSet.isUseEigenDecompInFit() ){ - parSet.propagateOriginalToEigen(); - - // also check the bounds of real parameter space - if( _reThrowParSetIfOutOfBounds_ ){ - for( auto& par : parSet.getEigenParameterList() ){ - if( not par.isEnabled() ) continue; - if( not par.isValueWithinBounds() ){ - // re-do the throwing - rethrow = true; - break; - } - } - } - } - } - - - if( rethrow ){ - // wrap back to the while loop -// LogDebug << "RE-THROW #" << throwNb << std::endl; - continue; - } - - // reached this point: all parameters are within bounds - keepThrowing = false; - } - - - } // Protected void Propagator::initializeThreads() { From 9f7a080b64a52e4175c59a2e0f3f126737a1a3ed Mon Sep 17 00:00:00 2001 From: lgiannes Date: Tue, 31 Oct 2023 11:57:20 +0100 Subject: [PATCH 008/131] fixed according to the update of Toolbox::ParallelWorker --- src/Applications/src/gundamMarginalise.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index ede6aa5ef..a9aed5bf7 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -97,8 +97,8 @@ int main(int argc, char** argv){ gRandom->SetSeed(seed); } - GundamGlobals::setNbThreads(clParser.getOptionVal("nbThreads", 1)); - LogInfo << "Running the fitter with " << GundamGlobals::getNbThreads() << " parallel threads." << std::endl; + GundamGlobals::getParallelWorker().setNThreads(clParser.getOptionVal("nbThreads", 1)); + LogInfo << "Running the fitter with " << GundamGlobals::getParallelWorker().getNbThreads() << " parallel threads." << std::endl; // Reading fitter file LogInfo << "Opening fitter output file: " << clParser.getOptionVal("fitterOutputFile") << std::endl; From 05f18a306192bd7b155db44fbe6c05274a8524eb Mon Sep 17 00:00:00 2001 From: lgiannes Date: Thu, 2 Nov 2023 11:38:46 +0100 Subject: [PATCH 009/131] added function in generictoolbox.root --- src/Applications/src/gundamMarginalise.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index a9aed5bf7..e9f9344ba 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -378,5 +378,5 @@ int main(int argc, char** argv){ margThrowTree->Write(); - GundamGlobals::getParallelWorker().reset(); + //GundamGlobals::getParallelWorker().reset(); } From 510329ff757130ac310fb4a13f9a8f0b4d0f5464 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Fri, 10 Nov 2023 18:38:51 +0100 Subject: [PATCH 010/131] print out swrt(det) of the cov. matrix --- src/Applications/src/gundamMarginalise.cxx | 90 +++++++++++++++++++--- 1 file changed, 79 insertions(+), 11 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index e9f9344ba..c0fe4160e 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -210,6 +210,10 @@ int main(int argc, char** argv){ } } }); + + // Compute the determinant + double det = propagator.getParametersManager().getGlobalCovarianceMatrix()->Determinant(); + // // Sample binning using parameterSetName // for( auto& sample : propagator.getFitSampleSet().getFitSampleList() ){ // auto associatedParSet = GenericToolbox::Json::fetchValue(sample.getConfig(), "parameterSetName"); @@ -289,8 +293,60 @@ int main(int argc, char** argv){ int nToys{ clParser.getOptionVal("nToys") }; // Get parameters to be marginalised - //std::string marginalisedParameters{GenericToolbox::Json::fetchValue(margConfig, "marginalisedParameters", "")}; + std::vector marginalisedParameters; + std::vector marginalisedParameterSets; + marginalisedParameters = GenericToolbox::Json::fetchValue>(margConfig, "parameterList"); + marginalisedParameterSets = GenericToolbox::Json::fetchValue>(margConfig, "parameterSetList"); + LogInfo<<"Marginalised parameters: "< type: " << par.getPriorType() << " mu=" << par.getPriorValue() + << " sigma= " << par.getStdDevValue() << " limits: " << par.getMinValue() << " - " + << par.getMaxValue() << " limits (phys): " << par.getMinPhysical() << " - " + << par.getMaxPhysical() << " limits (mirr): " << par.getMinMirror() << " - " + << par.getMaxMirror() <<" --- marg? "< type: "< "<Fill(); + + //debug +// if(iToy==0 || iToy==1){ +// LogInfo<<"DEBUG: Toy "<Write(); + det = 1.0; + TMatrixD eigenVectors = (*propagator.getParametersManager().getGlobalCovarianceMatrix()); + TVectorD eigenValues(parameters.size()); + eigenVectors.EigenVectors(eigenValues); + LogInfo<<"Eigenvalues: "< Date: Fri, 10 Nov 2023 18:40:25 +0100 Subject: [PATCH 011/131] added marginalise flag to Parameter class --- src/ParametersManager/include/Parameter.h | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/ParametersManager/include/Parameter.h b/src/ParametersManager/include/Parameter.h index 04fb42364..01db602aa 100644 --- a/src/ParametersManager/include/Parameter.h +++ b/src/ParametersManager/include/Parameter.h @@ -60,6 +60,7 @@ class Parameter : public JsonBaseClass { void setParameterDefinitionConfig(const nlohmann::json &config_); void setOwner(const ParameterSet *owner_){ _owner_ = owner_; } void setPriorType(PriorType::PriorType priorType){ _priorType_ = priorType; } + void setMarginalised(bool isMarginalised){ _isMarginalised_ = isMarginalised; } // Getters [[nodiscard]] bool isFree() const{ return _isFree_; } @@ -83,6 +84,7 @@ class Parameter : public JsonBaseClass { [[nodiscard]] const nlohmann::json &getDialDefinitionsList() const{ return _dialDefinitionsList_; } [[nodiscard]] const ParameterSet *getOwner() const{ return _owner_; } [[nodiscard]] PriorType::PriorType getPriorType() const{ return _priorType_; } + [[nodiscard]] bool isMarginalised() const{ return _isMarginalised_; } // Core void setValueAtPrior(); @@ -116,6 +118,7 @@ class Parameter : public JsonBaseClass { std::string _dialsWorkingDirectory_{"."}; nlohmann::json _parameterConfig_{}; nlohmann::json _dialDefinitionsList_{}; + bool _isMarginalised_{false}; // Internals const ParameterSet* _owner_{nullptr}; From 4008c6fad08b1c5be701fe5e4b5c89fd1efbcacf Mon Sep 17 00:00:00 2001 From: lgiannes Date: Fri, 10 Nov 2023 18:44:49 +0100 Subject: [PATCH 012/131] submodule generic toolbox --- submodules/cpp-generic-toolbox | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/submodules/cpp-generic-toolbox b/submodules/cpp-generic-toolbox index d7acf7f06..fa7710379 160000 --- a/submodules/cpp-generic-toolbox +++ b/submodules/cpp-generic-toolbox @@ -1 +1 @@ -Subproject commit d7acf7f06a1fe69b554cd1bf1555502dd62c3b66 +Subproject commit fa7710379eb8a2c8f079c13c768c01e677493c88 From 6278fd1b0b4ff9f8f32eb1ada80e0bb1ddd6566d Mon Sep 17 00:00:00 2001 From: lgiannes Date: Tue, 14 Nov 2023 09:47:33 +0100 Subject: [PATCH 013/131] debug parameterSetList --- src/Applications/src/gundamMarginalise.cxx | 35 +++++++++++++++------- 1 file changed, 24 insertions(+), 11 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index c0fe4160e..043743912 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -307,23 +307,31 @@ int main(int argc, char** argv){ if (not parSet.isEnabled()) { LogInfo << "Parameter set " << parSet.getName() << " is disabled" << std::endl; continue; - }else{ + } else { bool setMatches = false; + LogInfo << "Set: " << parSet.getName().c_str(); for (int i = 0; i < marginalisedParameterSets.size(); i++) { + if (0 == std::strcmp(parSet.getName().c_str(), marginalisedParameterSets[i].c_str())) { setMatches = (true); break; } else { setMatches = (false); + } } - if (setMatches){ + if (setMatches) { + LogInfo << " will be marginalized out. \n"; // loop over the parameters in the set and set the corresponding margThis to true for (auto &par: parSet.getParameterList()) { par.setMarginalised(true); } + continue; // skip the single params thing and go to the next ParameterSet + } else { + LogInfo << " will not be marginalized out. \n"; } } + // Do the same for single parameters for (auto &par: parSet.getParameterList()) { if (not par.isEnabled()) { LogInfo << "Parameter " << par.getName() << " is disabled" << std::endl; @@ -340,14 +348,19 @@ int main(int argc, char** argv){ } } par.setMarginalised(matches); - LogInfo << par.getFullTitle() << " -> type: " << par.getPriorType() << " mu=" << par.getPriorValue() - << " sigma= " << par.getStdDevValue() << " limits: " << par.getMinValue() << " - " - << par.getMaxValue() << " limits (phys): " << par.getMinPhysical() << " - " - << par.getMaxPhysical() << " limits (mirr): " << par.getMinMirror() << " - " - << par.getMaxMirror() <<" --- marg? "< type: " << par.getPriorType() << " mu=" << par.getPriorValue() +// << " sigma= " << par.getStdDevValue() << " limits: " << par.getMinValue() << " - " +// << par.getMaxValue() << " limits (phys): " << par.getMinPhysical() << " - " +// << par.getMaxPhysical() << " limits (mirr): " << par.getMinMirror() << " - " +// << par.getMaxMirror() <<" --- marg? "< Date: Tue, 14 Nov 2023 11:14:05 +0100 Subject: [PATCH 014/131] check if margConfig is a list --- src/Applications/src/gundamMarginalise.cxx | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 043743912..e4bce930f 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -114,6 +114,8 @@ int main(int argc, char** argv){ ObjectReader::readObject(fitterFile.get(), {{"gundam/config_TNamed"}, {"gundamFitter/unfoldedConfig_TNamed"}}, [&](TNamed* config_){ fitterConfig = GenericToolbox::Json::readConfigJsonStr( config_->GetTitle() ); }); + // Check if the config is an array (baobab compatibility) + ConfigUtils::ConfigHandler cHandler{ fitterConfig }; // // Disabling defined samples: @@ -137,6 +139,8 @@ int main(int argc, char** argv){ cHandler.override( margConfig ); LogInfo << "Override done." << std::endl; + + // read the parameters to include in the TTree // read the parameters to marginalise over @@ -295,8 +299,15 @@ int main(int argc, char** argv){ // Get parameters to be marginalised std::vector marginalisedParameters; std::vector marginalisedParameterSets; - marginalisedParameters = GenericToolbox::Json::fetchValue>(margConfig, "parameterList"); - marginalisedParameterSets = GenericToolbox::Json::fetchValue>(margConfig, "parameterSetList"); + if(margConfig.size()==1){ + // broken json library puts everything in the same level + marginalisedParameters = GenericToolbox::Json::fetchValue>(margConfig[0], "parameterList"); + marginalisedParameterSets = GenericToolbox::Json::fetchValue>(margConfig[0], "parameterSetList"); + }else{ + // normal behavior + marginalisedParameters = GenericToolbox::Json::fetchValue>(margConfig, "parameterList"); + marginalisedParameterSets = GenericToolbox::Json::fetchValue>(margConfig, "parameterSetList"); + } LogInfo<<"Marginalised parameters: "< Date: Tue, 14 Nov 2023 11:34:30 +0100 Subject: [PATCH 015/131] check if margConfig is a list --- src/Applications/src/gundamMarginalise.cxx | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index e4bce930f..44e0bced8 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -138,10 +138,14 @@ int main(int argc, char** argv){ nlohmann::json margConfig{ ConfigUtils::readConfigFile( clParser.getOptionVal("configFile") ) }; cHandler.override( margConfig ); LogInfo << "Override done." << std::endl; + if(margConfig.size()==1) { + // broken json library puts everything in the same level + // normal behavior + margConfig = margConfig[0]; + } - - // read the parameters to include in the TTree + // read the parameters to include in the TTree // read the parameters to marginalise over @@ -299,15 +303,9 @@ int main(int argc, char** argv){ // Get parameters to be marginalised std::vector marginalisedParameters; std::vector marginalisedParameterSets; - if(margConfig.size()==1){ - // broken json library puts everything in the same level - marginalisedParameters = GenericToolbox::Json::fetchValue>(margConfig[0], "parameterList"); - marginalisedParameterSets = GenericToolbox::Json::fetchValue>(margConfig[0], "parameterSetList"); - }else{ - // normal behavior marginalisedParameters = GenericToolbox::Json::fetchValue>(margConfig, "parameterList"); marginalisedParameterSets = GenericToolbox::Json::fetchValue>(margConfig, "parameterSetList"); - } + LogInfo<<"Marginalised parameters: "< Date: Tue, 14 Nov 2023 11:35:59 +0100 Subject: [PATCH 016/131] also print gLLH --- src/Applications/src/gundamMarginalise.cxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 44e0bced8..995f0f072 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -392,7 +392,6 @@ int main(int argc, char** argv){ //propagator.propagateParametersOnSamples(); // Probably not necessary (what's that for?) propagator.updateLlhCache(); LLH = propagator.getLlhBuffer(); - LogInfo<<"LLH: "< gLLH: "< Date: Tue, 14 Nov 2023 12:17:15 +0100 Subject: [PATCH 017/131] minor debug --- src/Applications/src/gundamMarginalise.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 995f0f072..2efed6733 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -136,13 +136,13 @@ int main(int argc, char** argv){ // Reading marginaliser config file nlohmann::json margConfig{ ConfigUtils::readConfigFile( clParser.getOptionVal("configFile") ) }; - cHandler.override( margConfig ); - LogInfo << "Override done." << std::endl; if(margConfig.size()==1) { // broken json library puts everything in the same level // normal behavior margConfig = margConfig[0]; } + cHandler.override( margConfig ); + LogInfo << "Override done." << std::endl; // read the parameters to include in the TTree From a04ae7275391b54bc6da761c5166003c6b235f8c Mon Sep 17 00:00:00 2001 From: lgiannes Date: Wed, 15 Nov 2023 16:45:03 +0100 Subject: [PATCH 018/131] implementation of profiling result in output TTree --- src/Applications/src/gundamMarginalise.cxx | 60 ++++++++++++++-------- 1 file changed, 40 insertions(+), 20 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 2efed6733..566d03b35 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -186,8 +186,9 @@ int main(int argc, char** argv){ } - + std::vector parametersBestFit; // Load post-fit parameters as "prior" so we can reset the weight to this point when throwing toys + // also save the values in a vector so we can use them to compute the LLH at the best fit point ObjectReader::readObject( fitterFile.get(), "FitterEngine/postFit/parState_TNamed", [&](TNamed* parState_){ propagator.getParametersManager().injectParameterValues( GenericToolbox::Json::readConfigJsonStr( parState_->GetTitle() ) ); for( auto& parSet : propagator.getParametersManager().getParameterSetsList() ){ @@ -196,6 +197,7 @@ int main(int argc, char** argv){ for( auto& par : parSet.getParameterList() ){ if( not par.isEnabled() ){ continue; } // LogInfo<<" "< "<("outputFile"); } + else { + // appendixDict["optionName"] = "Appendix" + // this list insure all appendices will appear in the same order + std::vector> appendixDict{ + {"configFile", "%s"}, + {"fitterOutputFile", "Fit_%s"}, + {"nToys", "nToys_%s"}, + {"randomSeed", "Seed_%s"}, + }; + outFilePath = GundamUtils::generateFileName(clParser, appendixDict) + ".root"; + + std::string outFolder{GenericToolbox::Json::fetchValue(margConfig, "outputFolder", "./")}; + outFilePath = GenericToolbox::joinPath(outFolder, outFilePath); + } // Load the post-fit covariance matrix ObjectReader::readObject( @@ -219,8 +238,6 @@ int main(int argc, char** argv){ } }); - // Compute the determinant - double det = propagator.getParametersManager().getGlobalCovarianceMatrix()->Determinant(); // // Sample binning using parameterSetName // for( auto& sample : propagator.getFitSampleSet().getFitSampleList() ){ @@ -248,29 +265,32 @@ int main(int argc, char** argv){ - // Creating output file - std::string outFilePath{}; - if( clParser.isOptionTriggered("outputFile") ){ outFilePath = clParser.getOptionVal("outputFile"); } - else { - // appendixDict["optionName"] = "Appendix" - // this list insure all appendices will appear in the same order - std::vector> appendixDict{ - {"configFile", "%s"}, - {"fitterOutputFile", "Fit_%s"}, - {"nToys", "nToys_%s"}, - {"randomSeed", "Seed_%s"}, - }; - outFilePath = "marg_" + GundamUtils::generateFileName(clParser, appendixDict) + ".root"; - std::string outFolder{GenericToolbox::Json::fetchValue(margConfig, "outputFolder", "./")}; - outFilePath = GenericToolbox::joinPath(outFolder, outFilePath); - } app.setCmdLinePtr( &clParser ); app.setConfigString( ConfigUtils::ConfigHandler{margConfig}.toString() ); app.openOutputFile( outFilePath ); app.writeAppInfo(); + // write the post fit covariance matrix in the output file + TDirectory* postFitInfo = app.getOutfilePtr()->mkdir("postFitInfo"); + postFitInfo->cd(); + ObjectReader::readObject( + fitterFile.get(), "FitterEngine/postFit/Hesse/hessian/postfitCovarianceOriginal_TH2D", + [&](TH2D* hCovPostFit_) { + hCovPostFit_->SetName("postFitCovarianceOriginal_TH2D"); + hCovPostFit_->Write();// save the TH2D cov. matrix also in the output file + }); + // write the best fit parameters vector to the output file + // TODO: implement using TH1D: (FitterEngine/postFit/Migrad/errors//values/postFitErrors_TH1D) + ObjectReader::readObject( fitterFile.get(), "FitterEngine/postFit/parState_TNamed", [&](TNamed* parState_){ + parState_->SetName("postFitParameters_TNamed"); + parState_->Write(); + }); + TVectorD bestFitParameters (parametersBestFit.size(), parametersBestFit.data()); + bestFitParameters.Write("bestFitParameters_TVectorD"); + app.getOutfilePtr()->cd(); + auto* marginalisationDir{ GenericToolbox::mkdirTFile(app.getOutfilePtr(), "marginalisation") }; LogInfo << "Creating throws tree" << std::endl; auto* margThrowTree = new TTree("margThrow", "margThrow"); @@ -457,7 +477,7 @@ int main(int argc, char** argv){ margThrowTree->Write(); - det = 1.0; + double det = 1.0; TMatrixD eigenVectors = (*propagator.getParametersManager().getGlobalCovarianceMatrix()); TVectorD eigenValues(parameters.size()); eigenVectors.EigenVectors(eigenValues); From c4ddfa6b331a58d1d44c94ff7a198f2e7928d935 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Wed, 15 Nov 2023 17:21:08 +0100 Subject: [PATCH 019/131] less verbose --- src/Applications/src/gundamMarginalise.cxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 566d03b35..db5e0578a 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -435,7 +435,8 @@ int main(int argc, char** argv){ iPar++; } } - LogInfo<<"LLH: "< gLLH: "< Date: Wed, 22 Nov 2023 11:55:32 +0100 Subject: [PATCH 020/131] less verbose --- src/Applications/src/gundamMarginalise.cxx | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index db5e0578a..9ec4b759e 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -187,6 +187,7 @@ int main(int argc, char** argv){ std::vector parametersBestFit; + std::vector parameterFullTitles; // Load post-fit parameters as "prior" so we can reset the weight to this point when throwing toys // also save the values in a vector so we can use them to compute the LLH at the best fit point ObjectReader::readObject( fitterFile.get(), "FitterEngine/postFit/parState_TNamed", [&](TNamed* parState_){ @@ -198,6 +199,7 @@ int main(int argc, char** argv){ if( not par.isEnabled() ){ continue; } // LogInfo<<" "< "<WriteObject(¶meterFullTitles, "parameterFullTitles"); app.getOutfilePtr()->cd(); + auto* marginalisationDir{ GenericToolbox::mkdirTFile(app.getOutfilePtr(), "marginalisation") }; LogInfo << "Creating throws tree" << std::endl; auto* margThrowTree = new TTree("margThrow", "margThrow"); From 5728fc75acb2d2e009fd3a102122da056ef8e358 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Mon, 11 Dec 2023 15:58:54 +0100 Subject: [PATCH 021/131] dummy test --- src/Applications/src/gundamMarginalise.cxx | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 9ec4b759e..ebfcdba31 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -34,6 +34,8 @@ LoggerInit([]{ int main(int argc, char** argv){ + std::cout<<"Hello world"< margThis; std::vector prior; std::vector weightsChiSquare; + std::vector dialResponse; // weightsChiSquare.reserve(1000); double LLH, gLLH, priorSum, LLHwrtBestFit; margThrowTree->Branch("Parameters", ¶meters); @@ -320,6 +323,7 @@ int main(int argc, char** argv){ margThrowTree->Branch("gLLH", &gLLH); margThrowTree->Branch("priorSum", &priorSum); margThrowTree->Branch("weightsChiSquare", &weightsChiSquare); + margThrowTree->Branch("dialResponse", &dialResponse); int nToys{ clParser.getOptionVal("nToys") }; @@ -413,7 +417,7 @@ int main(int argc, char** argv){ weightsChiSquare.clear(); // Do the throwing: propagator.getParametersManager().throwParametersFromGlobalCovariance(weightsChiSquare); - //propagator.propagateParametersOnSamples(); // Probably not necessary (what's that for?) + //propagator.propagateParametersOnSamples(); // Not necessary (included in updateLlhCache()) propagator.updateLlhCache(); LLH = propagator.getLlhBuffer(); LLHwrtBestFit = LLH - bestFitLLH; From b2a09d903b6d3abddb54adca077f7d4119cb3db3 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Mon, 11 Dec 2023 16:03:23 +0100 Subject: [PATCH 022/131] debug: print "dial response is out" --- src/Applications/src/gundamMarginalise.cxx | 5 ++--- .../DialEngine/src/DialResponseSupervisor.cpp | 12 ++++++++++-- 2 files changed, 12 insertions(+), 5 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index ebfcdba31..3f62d0516 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -34,7 +34,7 @@ LoggerInit([]{ int main(int argc, char** argv){ - std::cout<<"Hello world"< gLLH: "< double DialResponseSupervisor::process(double reponse_) const { // apply cap? - if ( not std::isnan(_minResponse_) and reponse_ < _minResponse_ ){ return _minResponse_; } - else if( not std::isnan(_maxResponse_) and reponse_ > _maxResponse_ ){ return _maxResponse_; } + // print out info if out of range + if ( not std::isnan(_minResponse_) and reponse_ < _minResponse_ ){ + LogWarning << "Response " << reponse_ << " is below the minimum response " << _minResponse_ << std::endl; + return _minResponse_; + } + else if( not std::isnan(_maxResponse_) and reponse_ > _maxResponse_ ){ + LogWarning << "Response " << reponse_ << " is above the maximum response " << _maxResponse_ << std::endl; + return _maxResponse_; + } // else? From a686686076a6e5d883183b99e6bd628ec3b57052 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Mon, 11 Dec 2023 16:14:31 +0100 Subject: [PATCH 023/131] also print out difference --- src/Applications/src/gundamMarginalise.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 3f62d0516..ae956567e 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -443,7 +443,7 @@ int main(int argc, char** argv){ iPar++; } } - LogInfo<<"LLH: "< gLLH: "< Date: Mon, 11 Dec 2023 16:26:59 +0100 Subject: [PATCH 024/131] always print our dial responses --- src/Applications/src/gundamMarginalise.cxx | 3 +++ .../DialEngine/src/DialResponseSupervisor.cpp | 6 +++--- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index ae956567e..7bde5a562 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -444,6 +444,9 @@ int main(int argc, char** argv){ } } LogInfo<<"LLH: "< 10){ + LogInfo<<"WARNING: BIG THROW!!"< gLLH: "< _maxResponse_ ){ - LogWarning << "Response " << reponse_ << " is above the maximum response " << _maxResponse_ << std::endl; return _maxResponse_; } From 8391e41379d72cd98cbd0a04d22e30dc42a126a6 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Mon, 11 Dec 2023 17:31:55 +0100 Subject: [PATCH 025/131] inject problematic parvector --- src/Applications/src/gundamMarginalise.cxx | 31 +++++++++++++++++-- .../DialEngine/src/DialResponseSupervisor.cpp | 1 - 2 files changed, 29 insertions(+), 3 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 7bde5a562..252109cf2 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -31,6 +31,8 @@ LoggerInit([]{ Logger::getUserHeader() << "[" << FILENAME << "]"; }); +double getParameterValueFromTextFile(std::string fileName, std::string parameterName); + int main(int argc, char** argv){ @@ -202,7 +204,7 @@ int main(int argc, char** argv){ // LogInfo<<" "< "<> parName >> parValue; + if (parName == parameterName) { + return parValue; + } + } + LogInfo << "Parameter " << parameterName << " not found in file " << fileName << std::endl; + return -999; +} \ No newline at end of file diff --git a/src/DialDictionary/DialEngine/src/DialResponseSupervisor.cpp b/src/DialDictionary/DialEngine/src/DialResponseSupervisor.cpp index ffa1969d4..9d37a356a 100644 --- a/src/DialDictionary/DialEngine/src/DialResponseSupervisor.cpp +++ b/src/DialDictionary/DialEngine/src/DialResponseSupervisor.cpp @@ -11,7 +11,6 @@ double DialResponseSupervisor::process(double reponse_) const { // apply cap? // print out info - LogWarning << "Response " << reponse_ << " limits: " << getSummary() << std::endl; if ( not std::isnan(_minResponse_) and reponse_ < _minResponse_ ){ return _minResponse_; From bc58bfc469040f7e5915c516cce03211ef396852 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Mon, 11 Dec 2023 19:15:59 +0100 Subject: [PATCH 026/131] file to inject parameters --- src/Applications/src/gundamMarginalise.cxx | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 252109cf2..f27fcdc3a 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -62,7 +62,7 @@ int main(int argc, char** argv){ clParser.addOption("nToys", {"-n"}, "Specify number of toys"); clParser.addOption("randomSeed", {"-s", "--seed"}, "Set random seed"); clParser.addTriggerOption("usingGpu", {"--gpu"}, "Use GPU parallelization"); - + clParser.addOption("parInject", {"--parameters-inject"}, "Input txt file for injecting params"); clParser.addDummyOption("Trigger options:"); clParser.addTriggerOption("dryRun", {"-d", "--dry-run"}, "Only overrides fitter config and print it."); @@ -160,6 +160,12 @@ int main(int argc, char** argv){ LogAlert << "Exiting as dry-run is set." << std::endl; return EXIT_SUCCESS; } + bool injectParamsManually = false; + std::string parInjectFile; + if( clParser.isOptionTriggered("parInject") ){ + parInjectFile = clParser.getOptionVal("parInject"); + injectParamsManually = true; + } auto configPropagator = GenericToolbox::Json::fetchValuePath( cHandler.getConfig(), "fitterEngineConfig/propagatorConfig" ); @@ -204,7 +210,11 @@ int main(int argc, char** argv){ // LogInfo<<" "< "< Date: Mon, 11 Dec 2023 19:26:51 +0100 Subject: [PATCH 027/131] uncommented throw params from global covariance --- src/Applications/src/gundamMarginalise.cxx | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index f27fcdc3a..abd946666 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -165,6 +165,9 @@ int main(int argc, char** argv){ if( clParser.isOptionTriggered("parInject") ){ parInjectFile = clParser.getOptionVal("parInject"); injectParamsManually = true; + LogInfo << "Injecting parameters from file: " << parInjectFile << std::endl; + }else{ + LogInfo << "Injecting best fit parameters as prior" << std::endl; } auto configPropagator = GenericToolbox::Json::fetchValuePath( cHandler.getConfig(), "fitterEngineConfig/propagatorConfig" ); @@ -428,7 +431,7 @@ int main(int argc, char** argv){ // reset weights vector weightsChiSquare.clear(); // Do the throwing: - //propagator.getParametersManager().throwParametersFromGlobalCovariance(weightsChiSquare); + propagator.getParametersManager().throwParametersFromGlobalCovariance(weightsChiSquare); //propagator.propagateParametersOnSamples(); // Not necessary (included in updateLlhCache()) propagator.updateLlhCache(); LLH = propagator.getLlhBuffer(); From 41184c80d536f6eb5f19b0b271b72019b320f9a9 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Mon, 11 Dec 2023 19:46:25 +0100 Subject: [PATCH 028/131] debugged parameter injector --- src/Applications/src/gundamMarginalise.cxx | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index abd946666..35a2bb681 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -520,7 +520,7 @@ int main(int argc, char** argv){ -double getParameterValueFromTextFile(std::string fileName, std::string parameterName) { +double getParameterValueFromTextFile(std::string fileName="LargeWeight_parVector.txt", std::string parameterName="Parameter Flux Systematics/#98") { std::ifstream file(fileName); std::string line; std::string parName; @@ -534,11 +534,12 @@ double getParameterValueFromTextFile(std::string fileName, std::string parameter std::getline(iss, title, ':'); std::getline(iss, value); - iss >> parName >> parValue; - if (parName == parameterName) { + + parValue = std::stod(value); + if ("Parameter "+title == parameterName) { return parValue; } } - LogInfo << "Parameter " << parameterName << " not found in file " << fileName << std::endl; + cout << "Parameter \"" << parameterName << "\" not found in file " << fileName << std::endl; return -999; } \ No newline at end of file From 08d0b73a65d0d7eca0bca156d901525e80c375ee Mon Sep 17 00:00:00 2001 From: lgiannes Date: Mon, 11 Dec 2023 19:49:15 +0100 Subject: [PATCH 029/131] debugged parameter injector --- src/Applications/src/gundamMarginalise.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 35a2bb681..6cbd06804 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -540,6 +540,6 @@ double getParameterValueFromTextFile(std::string fileName="LargeWeight_parVector return parValue; } } - cout << "Parameter \"" << parameterName << "\" not found in file " << fileName << std::endl; + LogInfo << "Parameter \"" << parameterName << "\" not found in file " << fileName << std::endl; return -999; } \ No newline at end of file From 1eec9f0007847824bba04dcdbcc7c9e252c83538 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Mon, 11 Dec 2023 19:56:24 +0100 Subject: [PATCH 030/131] debugged parameter injector --- src/Applications/src/gundamMarginalise.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 6cbd06804..d1914832a 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -536,7 +536,7 @@ double getParameterValueFromTextFile(std::string fileName="LargeWeight_parVector parValue = std::stod(value); - if ("Parameter "+title == parameterName) { + if (title == parameterName) { return parValue; } } From c33a2a6f64bff8282957b447e5bd52cd62ee7f1f Mon Sep 17 00:00:00 2001 From: lgiannes Date: Mon, 11 Dec 2023 20:12:54 +0100 Subject: [PATCH 031/131] debugged parameter injector --- src/Applications/src/gundamMarginalise.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index d1914832a..1d42574f3 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -458,7 +458,7 @@ int main(int argc, char** argv){ iPar++; } } - LogInfo<<"LLH: "< 10){ LogInfo<<"WARNING: BIG THROW!!"< Date: Mon, 11 Dec 2023 20:36:30 +0100 Subject: [PATCH 032/131] debugged parameter injector --- src/Applications/src/gundamMarginalise.cxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 1d42574f3..f280ccfc0 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -216,7 +216,8 @@ int main(int argc, char** argv){ if (not injectParamsManually) { par.setPriorValue(par.getParameterValue()); }else{ - par.setPriorValue(getParameterValueFromTextFile(parInjectFile,par.getFullTitle()) ); + par.setPriorValue(par.getParameterValue()); + par.setParameterValue(getParameterValueFromTextFile(parInjectFile,par.getFullTitle()) ); } } } From 36361f283874d4a010a42a5962c4f829a4aa14c5 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Mon, 11 Dec 2023 21:45:55 +0100 Subject: [PATCH 033/131] debugged parameter injector --- src/Applications/src/gundamMarginalise.cxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index f280ccfc0..23035b1bc 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -217,7 +217,8 @@ int main(int argc, char** argv){ par.setPriorValue(par.getParameterValue()); }else{ par.setPriorValue(par.getParameterValue()); - par.setParameterValue(getParameterValueFromTextFile(parInjectFile,par.getFullTitle()) ); + double parValue = getParameterValueFromTextFile(parInjectFile,par.getFullTitle()); + par.setParameterValue( parValue ); } } } From b81a156508439509c92cd3fb9d534a92b478b8e2 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Mon, 11 Dec 2023 21:46:34 +0100 Subject: [PATCH 034/131] debugged parameter injector --- src/Applications/src/gundamMarginalise.cxx | 1 + 1 file changed, 1 insertion(+) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 23035b1bc..19eba17ba 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -218,6 +218,7 @@ int main(int argc, char** argv){ }else{ par.setPriorValue(par.getParameterValue()); double parValue = getParameterValueFromTextFile(parInjectFile,par.getFullTitle()); + LogInfo<<"Setting "< "< Date: Mon, 11 Dec 2023 21:57:44 +0100 Subject: [PATCH 035/131] bug fix --- src/Applications/src/gundamMarginalise.cxx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 19eba17ba..5f840933b 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -218,7 +218,7 @@ int main(int argc, char** argv){ }else{ par.setPriorValue(par.getParameterValue()); double parValue = getParameterValueFromTextFile(parInjectFile,par.getFullTitle()); - LogInfo<<"Setting "< "< "< 10){ + LogInfo<<"LLH: "< 10){ LogInfo<<"WARNING: BIG THROW!!"< Date: Mon, 11 Dec 2023 23:00:48 +0100 Subject: [PATCH 036/131] bug fix --- src/Applications/src/gundamMarginalise.cxx | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 5f840933b..894d07820 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -217,6 +217,7 @@ int main(int argc, char** argv){ par.setPriorValue(par.getParameterValue()); }else{ par.setPriorValue(par.getParameterValue()); + // TODO: this is probably wrong, do it in "injectParameterValues" instead double parValue = getParameterValueFromTextFile(parInjectFile,par.getFullTitle()); LogInfo<<"Setting "< "< Date: Tue, 12 Dec 2023 10:36:38 +0100 Subject: [PATCH 037/131] debug. RUN WITH MAX 10 TOYS!! --- src/Applications/src/gundamMarginalise.cxx | 26 ++++++++-------------- 1 file changed, 9 insertions(+), 17 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 894d07820..c638edafc 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -213,15 +213,7 @@ int main(int argc, char** argv){ // LogInfo<<" "< "< "< Date: Tue, 12 Dec 2023 10:39:51 +0100 Subject: [PATCH 038/131] debug. RUN WITH MAX 10 TOYS!! --- src/Applications/src/gundamMarginalise.cxx | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index c638edafc..a928b3c50 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -430,12 +430,14 @@ int main(int argc, char** argv){ propagator.getParametersManager().throwParametersFromGlobalCovariance(weightsChiSquare); //propagator.propagateParametersOnSamples(); // Not necessary (included in updateLlhCache()) - for( auto& parSet : propagator.getParametersManager().getParameterSetsList() ){ - if( not parSet.isEnabled() ){ continue; } - for( auto& par : parSet.getParameterList() ){ - if( not par.isEnabled() ){ continue; } - par.setParameterValue(getParameterValueFromTextFile(parInjectFile, par.getFullTitle())); - LogInfo<<"Setting: "< Date: Tue, 12 Dec 2023 11:04:59 +0100 Subject: [PATCH 039/131] debug. RUN WITH MAX 10 TOYS!! --- src/Applications/src/gundamMarginalise.cxx | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index a928b3c50..b568bafc7 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -438,8 +438,10 @@ int main(int argc, char** argv){ par.setParameterValue(getParameterValueFromTextFile(parInjectFile, par.getFullTitle())); LogInfo << "Setting: " << par.getFullTitle() << " to " << par.getParameterValue() << std::endl; } + propagator.getParametersManager().injectParameterValues( parSet.exportInjectorConfig() ); } } + LogInfo<<"Computing LH..."< Date: Tue, 12 Dec 2023 11:06:34 +0100 Subject: [PATCH 040/131] debug. RUN WITH MAX 10 TOYS!! --- src/Applications/src/gundamMarginalise.cxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index b568bafc7..68b4f8eff 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -436,8 +436,9 @@ int main(int argc, char** argv){ for (auto &par: parSet.getParameterList()) { if (not par.isEnabled()) { continue; } par.setParameterValue(getParameterValueFromTextFile(parInjectFile, par.getFullTitle())); - LogInfo << "Setting: " << par.getFullTitle() << " to " << par.getParameterValue() << std::endl; + //LogInfo << "Setting: " << par.getFullTitle() << " to " << par.getParameterValue() << std::endl; } + LogInfo << "Injecting parameterSet values for parSet "< Date: Tue, 12 Dec 2023 11:23:00 +0100 Subject: [PATCH 041/131] debug. RUN WITH MAX 10 TOYS!! --- src/Applications/src/gundamMarginalise.cxx | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 68b4f8eff..a78ebc0ac 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -438,14 +438,13 @@ int main(int argc, char** argv){ par.setParameterValue(getParameterValueFromTextFile(parInjectFile, par.getFullTitle())); //LogInfo << "Setting: " << par.getFullTitle() << " to " << par.getParameterValue() << std::endl; } - LogInfo << "Injecting parameterSet values for parSet "< 10){ - LogInfo<<"WARNING: BIG THROW!!"< 10){ +// LogInfo<<"WARNING: BIG THROW!!"< gLLH: "< Date: Tue, 12 Dec 2023 11:25:02 +0100 Subject: [PATCH 042/131] debug. RUN WITH MAX 10 TOYS!! --- src/Applications/src/gundamMarginalise.cxx | 1 + 1 file changed, 1 insertion(+) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index a78ebc0ac..8dff46e92 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -431,6 +431,7 @@ int main(int argc, char** argv){ //propagator.propagateParametersOnSamples(); // Not necessary (included in updateLlhCache()) if(injectParamsManually) { + LogInfo<< "Injecting parameters from file: " << parInjectFile << std::endl; for (auto &parSet: propagator.getParametersManager().getParameterSetsList()) { if (not parSet.isEnabled()) { continue; } for (auto &par: parSet.getParameterList()) { From c47c4ba693d9b7545b0356b78666678e9baf0291 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Tue, 12 Dec 2023 11:30:12 +0100 Subject: [PATCH 043/131] debug (getEffectiveParameterList). RUN WITH MAX 10 TOYS!! --- src/Applications/src/gundamMarginalise.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 8dff46e92..e48dac962 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -434,7 +434,7 @@ int main(int argc, char** argv){ LogInfo<< "Injecting parameters from file: " << parInjectFile << std::endl; for (auto &parSet: propagator.getParametersManager().getParameterSetsList()) { if (not parSet.isEnabled()) { continue; } - for (auto &par: parSet.getParameterList()) { + for (auto &par: parSet.getEffectiveParameterList()) { if (not par.isEnabled()) { continue; } par.setParameterValue(getParameterValueFromTextFile(parInjectFile, par.getFullTitle())); //LogInfo << "Setting: " << par.getFullTitle() << " to " << par.getParameterValue() << std::endl; From 5ae550d5aca59f98ea3099d5b725b5dd6e0350ae Mon Sep 17 00:00:00 2001 From: lgiannes Date: Tue, 12 Dec 2023 11:35:35 +0100 Subject: [PATCH 044/131] debug (getEffectiveParameterList). RUN WITH MAX 10 TOYS!! --- src/Applications/src/gundamMarginalise.cxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index e48dac962..5be0984d3 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -434,11 +434,12 @@ int main(int argc, char** argv){ LogInfo<< "Injecting parameters from file: " << parInjectFile << std::endl; for (auto &parSet: propagator.getParametersManager().getParameterSetsList()) { if (not parSet.isEnabled()) { continue; } - for (auto &par: parSet.getEffectiveParameterList()) { + for (auto &par: parSet.getParameterList()) { if (not par.isEnabled()) { continue; } par.setParameterValue(getParameterValueFromTextFile(parInjectFile, par.getFullTitle())); //LogInfo << "Setting: " << par.getFullTitle() << " to " << par.getParameterValue() << std::endl; } + parSet.propagateOriginalToEigen(); } } From 5bc2ccdb6eee035d58dc13b988e105978a069bd9 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Tue, 12 Dec 2023 11:46:46 +0100 Subject: [PATCH 045/131] debug (getEffectiveParameterList). RUN WITH MAX 10 TOYS!! --- src/Applications/src/gundamMarginalise.cxx | 1 - 1 file changed, 1 deletion(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 5be0984d3..8dff46e92 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -439,7 +439,6 @@ int main(int argc, char** argv){ par.setParameterValue(getParameterValueFromTextFile(parInjectFile, par.getFullTitle())); //LogInfo << "Setting: " << par.getFullTitle() << " to " << par.getParameterValue() << std::endl; } - parSet.propagateOriginalToEigen(); } } From f7b38cac8d7fc899d812515adf4c789fd5e6053e Mon Sep 17 00:00:00 2001 From: lgiannes Date: Tue, 12 Dec 2023 12:02:44 +0100 Subject: [PATCH 046/131] debug (getEffectiveParameterList). RUN WITH MAX 10 TOYS!! --- src/Applications/src/gundamMarginalise.cxx | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 8dff46e92..4e81f6855 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -439,6 +439,9 @@ int main(int argc, char** argv){ par.setParameterValue(getParameterValueFromTextFile(parInjectFile, par.getFullTitle())); //LogInfo << "Setting: " << par.getFullTitle() << " to " << par.getParameterValue() << std::endl; } + if (parSet.isUseEigenDecompInFit()){ + parSet.propagateEigenToOriginal(); + } } } From 1fc80c4ca6583296096be03c599ef702b6dbe152 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Tue, 12 Dec 2023 12:11:58 +0100 Subject: [PATCH 047/131] debug (getEffectiveParameterList). RUN WITH MAX 10 TOYS!! --- src/Applications/src/gundamMarginalise.cxx | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 4e81f6855..f0d890a02 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -432,13 +432,22 @@ int main(int argc, char** argv){ if(injectParamsManually) { LogInfo<< "Injecting parameters from file: " << parInjectFile << std::endl; + + std::vector parPtrList{}; + // save a list of parameter pointers (NOT IN EIGEN SPACE) + for( auto& parSet : propagator.getParametersManager().getParameterSetsList() ){ + for( auto& par : parSet.getParameterList() ){ + if( par.isEnabled() and not par.isFixed() ) { + parPtrList.emplace_back(&par); + // set the parameter value from the injector file + parPtrList.back()->setParameterValue(getParameterValueFromTextFile(parInjectFile, par.getFullTitle())); + } + } + } + + // If is in eigen space, propagateOriginalToEigen for (auto &parSet: propagator.getParametersManager().getParameterSetsList()) { if (not parSet.isEnabled()) { continue; } - for (auto &par: parSet.getParameterList()) { - if (not par.isEnabled()) { continue; } - par.setParameterValue(getParameterValueFromTextFile(parInjectFile, par.getFullTitle())); - //LogInfo << "Setting: " << par.getFullTitle() << " to " << par.getParameterValue() << std::endl; - } if (parSet.isUseEigenDecompInFit()){ parSet.propagateEigenToOriginal(); } From d010a620a1699ba6b3799d93f48e993c9463b347 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Tue, 12 Dec 2023 13:16:52 +0100 Subject: [PATCH 048/131] debug (getEffectiveParameterList). RUN WITH MAX 10 TOYS!! --- src/Applications/src/gundamMarginalise.cxx | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index f0d890a02..b2160abd8 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -441,6 +441,7 @@ int main(int argc, char** argv){ parPtrList.emplace_back(&par); // set the parameter value from the injector file parPtrList.back()->setParameterValue(getParameterValueFromTextFile(parInjectFile, par.getFullTitle())); + LogInfo<<"Setting: "< Date: Tue, 12 Dec 2023 13:18:05 +0100 Subject: [PATCH 049/131] debug (getEffectiveParameterList). RUN WITH MAX 10 TOYS!! --- src/Applications/src/gundamMarginalise.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index b2160abd8..8a0b7b6bf 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -427,7 +427,7 @@ int main(int argc, char** argv){ // reset weights vector weightsChiSquare.clear(); // Do the throwing: - propagator.getParametersManager().throwParametersFromGlobalCovariance(weightsChiSquare); + //propagator.getParametersManager().throwParametersFromGlobalCovariance(weightsChiSquare); //propagator.propagateParametersOnSamples(); // Not necessary (included in updateLlhCache()) if(injectParamsManually) { From a4e9544a77839166052810a0a023e9b59bab0592 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Tue, 12 Dec 2023 13:22:19 +0100 Subject: [PATCH 050/131] debug (getEffectiveParameterList). RUN WITH MAX 10 TOYS!! --- src/Applications/src/gundamMarginalise.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 8a0b7b6bf..b2160abd8 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -427,7 +427,7 @@ int main(int argc, char** argv){ // reset weights vector weightsChiSquare.clear(); // Do the throwing: - //propagator.getParametersManager().throwParametersFromGlobalCovariance(weightsChiSquare); + propagator.getParametersManager().throwParametersFromGlobalCovariance(weightsChiSquare); //propagator.propagateParametersOnSamples(); // Not necessary (included in updateLlhCache()) if(injectParamsManually) { From 355c9784144476c1d167f82747afde2df8e2f9c7 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Tue, 12 Dec 2023 13:38:56 +0100 Subject: [PATCH 051/131] debug (getEffectiveParameterList). RUN WITH MAX 10 TOYS!! --- src/Applications/src/gundamMarginalise.cxx | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index b2160abd8..f504fb5f5 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -433,6 +433,15 @@ int main(int argc, char** argv){ if(injectParamsManually) { LogInfo<< "Injecting parameters from file: " << parInjectFile << std::endl; + + + // If is in eigen space, propagateOriginalToEigen + for (auto &parSet: propagator.getParametersManager().getParameterSetsList()) { + if (not parSet.isEnabled()) { continue; } + if (parSet.isUseEigenDecompInFit()){ + parSet.propagateEigenToOriginal(); + } + } std::vector parPtrList{}; // save a list of parameter pointers (NOT IN EIGEN SPACE) for( auto& parSet : propagator.getParametersManager().getParameterSetsList() ){ @@ -442,18 +451,11 @@ int main(int argc, char** argv){ // set the parameter value from the injector file parPtrList.back()->setParameterValue(getParameterValueFromTextFile(parInjectFile, par.getFullTitle())); LogInfo<<"Setting: "< Date: Tue, 12 Dec 2023 13:49:59 +0100 Subject: [PATCH 052/131] debug (getEffectiveParameterList). RUN WITH MAX 10 TOYS!! --- src/Applications/src/gundamMarginalise.cxx | 39 +++++++++++++--------- 1 file changed, 24 insertions(+), 15 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index f504fb5f5..3ca135585 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -433,8 +433,30 @@ int main(int argc, char** argv){ if(injectParamsManually) { LogInfo<< "Injecting parameters from file: " << parInjectFile << std::endl; - - + // count the number of parameters to be injected + int nStripped{0}; + for( auto& parSet : propagator.getParametersManager().getParameterSetsList() ){ + if( not parSet.isEnabled() ) continue; + for( auto& par : parSet.getParameterList() ){ + if( not par.isEnabled() ) continue; + nStripped++; + } + } + // allocate a vector of parameter pointers + std::vector strippedParameterList.reserve( nStripped ); + for( auto& parSet : propagator.getParametersManager().getParameterSetsList() ){ + if( not parSet.isEnabled() ) continue; + for( auto& par : parSet.getParameterList() ){ + if( not par.isEnabled() ) continue; + strippedParameterList.emplace_back(&par); + } + } + // change their values + for( int iPar = 0 ; iPar < nStripped ; iPar++ ) { + strippedParameterList[iPar]->setParameterValue( + getParameterValueFromTextFile(parInjectFile, strippedParameterList[iPar]->getFullTitle()) + ); + } // If is in eigen space, propagateOriginalToEigen for (auto &parSet: propagator.getParametersManager().getParameterSetsList()) { if (not parSet.isEnabled()) { continue; } @@ -442,19 +464,6 @@ int main(int argc, char** argv){ parSet.propagateEigenToOriginal(); } } - std::vector parPtrList{}; - // save a list of parameter pointers (NOT IN EIGEN SPACE) - for( auto& parSet : propagator.getParametersManager().getParameterSetsList() ){ - for( auto& par : parSet.getParameterList() ){ - if( par.isEnabled() and not par.isFixed() ) { - parPtrList.emplace_back(&par); - // set the parameter value from the injector file - parPtrList.back()->setParameterValue(getParameterValueFromTextFile(parInjectFile, par.getFullTitle())); - LogInfo<<"Setting: "< Date: Tue, 12 Dec 2023 13:53:12 +0100 Subject: [PATCH 053/131] debug (getEffectiveParameterList). RUN WITH MAX 10 TOYS!! --- src/Applications/src/gundamMarginalise.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 3ca135585..3af6ce99d 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -427,7 +427,7 @@ int main(int argc, char** argv){ // reset weights vector weightsChiSquare.clear(); // Do the throwing: - propagator.getParametersManager().throwParametersFromGlobalCovariance(weightsChiSquare); + //propagator.getParametersManager().throwParametersFromGlobalCovariance(weightsChiSquare); //propagator.propagateParametersOnSamples(); // Not necessary (included in updateLlhCache()) if(injectParamsManually) { @@ -451,7 +451,7 @@ int main(int argc, char** argv){ strippedParameterList.emplace_back(&par); } } - // change their values + // change the parameter values for( int iPar = 0 ; iPar < nStripped ; iPar++ ) { strippedParameterList[iPar]->setParameterValue( getParameterValueFromTextFile(parInjectFile, strippedParameterList[iPar]->getFullTitle()) From 90d9a3416e61e4b676741711af54c316269122c5 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Tue, 12 Dec 2023 13:54:07 +0100 Subject: [PATCH 054/131] debug (getEffectiveParameterList). RUN WITH MAX 10 TOYS!! --- src/Applications/src/gundamMarginalise.cxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 3af6ce99d..8f1282334 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -443,7 +443,8 @@ int main(int argc, char** argv){ } } // allocate a vector of parameter pointers - std::vector strippedParameterList.reserve( nStripped ); + std::vector strippedParameterList; + strippedParameterList.reserve( nStripped ); for( auto& parSet : propagator.getParametersManager().getParameterSetsList() ){ if( not parSet.isEnabled() ) continue; for( auto& par : parSet.getParameterList() ){ From bc595072d5d6726ceddac7bfb3901a5a914341c3 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Tue, 12 Dec 2023 13:57:12 +0100 Subject: [PATCH 055/131] debug (getEffectiveParameterList). RUN WITH MAX 10 TOYS!! --- src/Applications/src/gundamMarginalise.cxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 8f1282334..be31dd82a 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -427,7 +427,7 @@ int main(int argc, char** argv){ // reset weights vector weightsChiSquare.clear(); // Do the throwing: - //propagator.getParametersManager().throwParametersFromGlobalCovariance(weightsChiSquare); + propagator.getParametersManager().throwParametersFromGlobalCovariance(weightsChiSquare); //propagator.propagateParametersOnSamples(); // Not necessary (included in updateLlhCache()) if(injectParamsManually) { @@ -458,6 +458,7 @@ int main(int argc, char** argv){ getParameterValueFromTextFile(parInjectFile, strippedParameterList[iPar]->getFullTitle()) ); } + // print out the parameter values // If is in eigen space, propagateOriginalToEigen for (auto &parSet: propagator.getParametersManager().getParameterSetsList()) { if (not parSet.isEnabled()) { continue; } From 66efaf2359ca8b82976ed03c9c0053eb3bcd85f4 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Tue, 12 Dec 2023 18:20:22 +0100 Subject: [PATCH 056/131] print out prior info --- src/Applications/src/gundamMarginalise.cxx | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index be31dd82a..7d4ff859b 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -397,9 +397,21 @@ int main(int argc, char** argv){ } par.setMarginalised(matches); if(par.isMarginalised()){ - LogInfo << "Parameter " << par.getFullTitle() << " will be marginalized out. \n"; + LogInfo << "Parameter " << par.getFullTitle() + << " -> type: " << par.getPriorType() << " mu=" << par.getPriorValue() + << " sigma= " << par.getStdDevValue() << " limits: " << par.getMinValue() << " - " + << par.getMaxValue() << " limits (phys): " << par.getMinPhysical() << " - " + << par.getMaxPhysical() << " limits (mirr): " << par.getMinMirror() << " - " + << par.getMaxMirror() <<" --- marg? "< type: " << par.getPriorType() << " mu=" << par.getPriorValue() + << " sigma= " << par.getStdDevValue() << " limits: " << par.getMinValue() << " - " + << par.getMaxValue() << " limits (phys): " << par.getMinPhysical() << " - " + << par.getMaxPhysical() << " limits (mirr): " << par.getMinMirror() << " - " + << par.getMaxMirror() <<" --- marg? "< Date: Tue, 12 Dec 2023 18:27:45 +0100 Subject: [PATCH 057/131] print out prior info --- src/Applications/src/gundamMarginalise.cxx | 13 ++++--------- 1 file changed, 4 insertions(+), 9 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 7d4ff859b..79a071596 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -352,11 +352,11 @@ int main(int argc, char** argv){ // Also display the prior information LogInfo<<"Prior information: "< type: " << par.getPriorType() << " mu=" << par.getPriorValue() << " sigma= " << par.getStdDevValue() << " limits: " << par.getMinValue() << " - " - << par.getMaxValue() << " limits (phys): " << par.getMinPhysical() << " - " - << par.getMaxPhysical() << " limits (mirr): " << par.getMinMirror() << " - " - << par.getMaxMirror() <<" --- marg? "< type: " << par.getPriorType() << " mu=" << par.getPriorValue() << " sigma= " << par.getStdDevValue() << " limits: " << par.getMinValue() << " - " - << par.getMaxValue() << " limits (phys): " << par.getMinPhysical() << " - " - << par.getMaxPhysical() << " limits (mirr): " << par.getMinMirror() << " - " - << par.getMaxMirror() <<" --- marg? "< Date: Tue, 12 Dec 2023 18:34:50 +0100 Subject: [PATCH 058/131] print out prior info --- src/Applications/src/gundamMarginalise.cxx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 79a071596..19a608898 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -380,7 +380,7 @@ int main(int argc, char** argv){ } // Do the same for single parameters for (auto &par: parSet.getParameterList()) { - if (not par.isEnabled() and not setMatches) { + if (not par.isEnabled()) { LogInfo << "Parameter " << par.getName() << " is disabled" << std::endl; continue; @@ -394,19 +394,19 @@ int main(int argc, char** argv){ matches = (false); } } - par.setMarginalised(matches); + par.setMarginalised(matches or setMatches); if(par.isMarginalised()){ LogInfo << "Parameter " << par.getFullTitle() << " -> type: " << par.getPriorType() << " mu=" << par.getPriorValue() << " sigma= " << par.getStdDevValue() << " limits: " << par.getMinValue() << " - " << par.getMaxValue() < type: " << par.getPriorType() << " mu=" << par.getPriorValue() << " sigma= " << par.getStdDevValue() << " limits: " << par.getMinValue() << " - " << par.getMaxValue() << par.isMarginalised() - <<" will NOT be marginalized out. \n"; + <<" will NOT be marg. out\n"; } } } From 1b943f9b2ebc381410045d567d87ca2b617e6c97 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Wed, 13 Dec 2023 16:06:57 +0100 Subject: [PATCH 059/131] check param value, does it change? --- src/Applications/src/gundamMarginalise.cxx | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 19a608898..ee953afd6 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -399,14 +399,14 @@ int main(int argc, char** argv){ LogInfo << "Parameter " << par.getFullTitle() << " -> type: " << par.getPriorType() << " mu=" << par.getPriorValue() << " sigma= " << par.getStdDevValue() << " limits: " << par.getMinValue() << " - " - << par.getMaxValue() < will be marg. out\n"; }else{ LogInfo << "Parameter " << par.getFullTitle() << " -> type: " << par.getPriorType() << " mu=" << par.getPriorValue() << " sigma= " << par.getStdDevValue() << " limits: " << par.getMinValue() << " - " - << par.getMaxValue() << par.isMarginalised() - <<" will NOT be marg. out\n"; + << par.getMaxValue() + <<" -> will NOT be marg. out\n"; } } } @@ -457,6 +457,9 @@ int main(int argc, char** argv){ for( auto& par : parSet.getParameterList() ){ if( not par.isEnabled() ) continue; strippedParameterList.emplace_back(&par); + if (par.getFullTitle()=="Parameter Flux Systematics/#98"){ + LogInfo<<"Parameter "<setParameterValue( getParameterValueFromTextFile(parInjectFile, strippedParameterList[iPar]->getFullTitle()) ); + if (strippedParameterList[iPar]->getFullTitle()=="Parameter Flux Systematics/#98"){ + LogInfo<<"Parameter "< Date: Wed, 13 Dec 2023 16:10:19 +0100 Subject: [PATCH 060/131] check param value, does it change? --- src/Applications/src/gundamMarginalise.cxx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index ee953afd6..da8006b26 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -468,7 +468,7 @@ int main(int argc, char** argv){ getParameterValueFromTextFile(parInjectFile, strippedParameterList[iPar]->getFullTitle()) ); if (strippedParameterList[iPar]->getFullTitle()=="Parameter Flux Systematics/#98"){ - LogInfo<<"Parameter "<getFullTitle()<<" has value "<getParameterValue()< gLLH: "<Write(); From dea72b80256f0f979cd856f8e4b201d7c4319b03 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Wed, 13 Dec 2023 16:16:28 +0100 Subject: [PATCH 061/131] check param value, does it change? --- src/Applications/src/gundamMarginalise.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index da8006b26..2cbaf6b5e 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -457,7 +457,7 @@ int main(int argc, char** argv){ for( auto& par : parSet.getParameterList() ){ if( not par.isEnabled() ) continue; strippedParameterList.emplace_back(&par); - if (par.getFullTitle()=="Parameter Flux Systematics/#98"){ + if (par.getFullTitle()=="ND280 Detector Systematics/#551"){ LogInfo<<"Parameter "< Date: Wed, 13 Dec 2023 16:19:59 +0100 Subject: [PATCH 062/131] check param value, does it change? --- src/Applications/src/gundamMarginalise.cxx | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 2cbaf6b5e..bd386ecdd 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -457,9 +457,8 @@ int main(int argc, char** argv){ for( auto& par : parSet.getParameterList() ){ if( not par.isEnabled() ) continue; strippedParameterList.emplace_back(&par); - if (par.getFullTitle()=="ND280 Detector Systematics/#551"){ - LogInfo<<"Parameter "<setParameterValue( getParameterValueFromTextFile(parInjectFile, strippedParameterList[iPar]->getFullTitle()) ); - if (strippedParameterList[iPar]->getFullTitle()=="Parameter Flux Systematics/#98"){ - LogInfo<<"Parameter "<getFullTitle()<<" has value "<getParameterValue()<getFullTitle()<<" has value "<getParameterValue()< Date: Wed, 13 Dec 2023 16:28:35 +0100 Subject: [PATCH 063/131] check param value, does it change? --- src/Applications/src/gundamMarginalise.cxx | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index bd386ecdd..8d3b7d369 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -457,8 +457,10 @@ int main(int argc, char** argv){ for( auto& par : parSet.getParameterList() ){ if( not par.isEnabled() ) continue; strippedParameterList.emplace_back(&par); - LogInfo<setParameterValue( getParameterValueFromTextFile(parInjectFile, strippedParameterList[iPar]->getFullTitle()) ); - LogInfo<getFullTitle()<<" has value "<getParameterValue()<getFullTitle() == "ND280 Detector Systematics/#551" or + strippedParameterList[iPar]->getFullTitle() == "Cross-Section (binned) Systematics/#0_EB_bin_C_nu") { + LogInfo << strippedParameterList[iPar]->getFullTitle() << " has value " + << strippedParameterList[iPar]->getParameterValue() << std::endl; + } } // print out the parameter values // If is in eigen space, propagateOriginalToEigen @@ -477,7 +483,10 @@ int main(int argc, char** argv){ } for( auto& par : parSet.getParameterList() ){ if( not par.isEnabled() ) continue; - LogInfo< Date: Wed, 13 Dec 2023 16:37:04 +0100 Subject: [PATCH 064/131] check param value, does it change? --- src/Applications/src/gundamMarginalise.cxx | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 8d3b7d369..f4dab70e5 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -458,7 +458,8 @@ int main(int argc, char** argv){ if( not par.isEnabled() ) continue; strippedParameterList.emplace_back(&par); if (par.getFullTitle()=="ND280 Detector Systematics/#551" or - par.getFullTitle()=="Cross-Section (binned) Systematics/#0_EB_bin_C_nu") { + par.getFullTitle()=="Cross-Section (binned) Systematics/#0_EB_bin_C_nu" or + par.getFullTitle()=="Flux Systematics/#8") { LogInfo << par.getFullTitle() << " has value " << par.getParameterValue() << std::endl; } } @@ -469,7 +470,8 @@ int main(int argc, char** argv){ getParameterValueFromTextFile(parInjectFile, strippedParameterList[iPar]->getFullTitle()) ); if (strippedParameterList[iPar]->getFullTitle() == "ND280 Detector Systematics/#551" or - strippedParameterList[iPar]->getFullTitle() == "Cross-Section (binned) Systematics/#0_EB_bin_C_nu") { + strippedParameterList[iPar]->getFullTitle() == "Cross-Section (binned) Systematics/#0_EB_bin_C_nu" or + strippedParameterList[iPar]->getFullTitle() == "Flux Systematics/#8") { LogInfo << strippedParameterList[iPar]->getFullTitle() << " has value " << strippedParameterList[iPar]->getParameterValue() << std::endl; } @@ -480,11 +482,13 @@ int main(int argc, char** argv){ if (not parSet.isEnabled()) { continue; } if (parSet.isUseEigenDecompInFit()){ parSet.propagateEigenToOriginal(); + LogInfo<< "propagating to original space: "< Date: Wed, 13 Dec 2023 16:57:55 +0100 Subject: [PATCH 065/131] check param value, does it change? --- src/Applications/src/gundamMarginalise.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index f4dab70e5..c434f4aab 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -481,8 +481,8 @@ int main(int argc, char** argv){ for (auto &parSet: propagator.getParametersManager().getParameterSetsList()) { if (not parSet.isEnabled()) { continue; } if (parSet.isUseEigenDecompInFit()){ - parSet.propagateEigenToOriginal(); - LogInfo<< "propagating to original space: "< Date: Wed, 13 Dec 2023 17:02:42 +0100 Subject: [PATCH 066/131] check param value, does it change? --- src/Applications/src/gundamMarginalise.cxx | 21 --------------------- 1 file changed, 21 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index c434f4aab..91bf028fb 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -457,11 +457,6 @@ int main(int argc, char** argv){ for( auto& par : parSet.getParameterList() ){ if( not par.isEnabled() ) continue; strippedParameterList.emplace_back(&par); - if (par.getFullTitle()=="ND280 Detector Systematics/#551" or - par.getFullTitle()=="Cross-Section (binned) Systematics/#0_EB_bin_C_nu" or - par.getFullTitle()=="Flux Systematics/#8") { - LogInfo << par.getFullTitle() << " has value " << par.getParameterValue() << std::endl; - } } } // change the parameter values @@ -469,12 +464,6 @@ int main(int argc, char** argv){ strippedParameterList[iPar]->setParameterValue( getParameterValueFromTextFile(parInjectFile, strippedParameterList[iPar]->getFullTitle()) ); - if (strippedParameterList[iPar]->getFullTitle() == "ND280 Detector Systematics/#551" or - strippedParameterList[iPar]->getFullTitle() == "Cross-Section (binned) Systematics/#0_EB_bin_C_nu" or - strippedParameterList[iPar]->getFullTitle() == "Flux Systematics/#8") { - LogInfo << strippedParameterList[iPar]->getFullTitle() << " has value " - << strippedParameterList[iPar]->getParameterValue() << std::endl; - } } // print out the parameter values // If is in eigen space, propagateOriginalToEigen @@ -482,16 +471,6 @@ int main(int argc, char** argv){ if (not parSet.isEnabled()) { continue; } if (parSet.isUseEigenDecompInFit()){ parSet.propagateOriginalToEigen(); - LogInfo<< "propagating to eigen space: "< Date: Wed, 13 Dec 2023 17:15:14 +0100 Subject: [PATCH 067/131] check param value, does it change? --- src/Applications/src/gundamMarginalise.cxx | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 91bf028fb..5d565e35c 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -434,7 +434,7 @@ int main(int argc, char** argv){ // reset weights vector weightsChiSquare.clear(); // Do the throwing: - propagator.getParametersManager().throwParametersFromGlobalCovariance(weightsChiSquare); + //propagator.getParametersManager().throwParametersFromGlobalCovariance(weightsChiSquare); //propagator.propagateParametersOnSamples(); // Not necessary (included in updateLlhCache()) if(injectParamsManually) { @@ -461,9 +461,11 @@ int main(int argc, char** argv){ } // change the parameter values for( int iPar = 0 ; iPar < nStripped ; iPar++ ) { + double epsilon = gRandom->Gaus(0, 0.001); strippedParameterList[iPar]->setParameterValue( - getParameterValueFromTextFile(parInjectFile, strippedParameterList[iPar]->getFullTitle()) + epsilon + getParameterValueFromTextFile(parInjectFile, strippedParameterList[iPar]->getFullTitle()) ); + } // print out the parameter values // If is in eigen space, propagateOriginalToEigen From beff0c6fa30fb95b9a2e5448d20a0982ba27f57b Mon Sep 17 00:00:00 2001 From: lgiannes Date: Wed, 13 Dec 2023 17:21:36 +0100 Subject: [PATCH 068/131] check param value, does it change? --- src/Applications/src/gundamMarginalise.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 5d565e35c..1a2279a74 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -461,7 +461,7 @@ int main(int argc, char** argv){ } // change the parameter values for( int iPar = 0 ; iPar < nStripped ; iPar++ ) { - double epsilon = gRandom->Gaus(0, 0.001); + double epsilon = gRandom->Gaus(0, 0.000001); strippedParameterList[iPar]->setParameterValue( epsilon + getParameterValueFromTextFile(parInjectFile, strippedParameterList[iPar]->getFullTitle()) ); From 599fdcf1fabbf700821161983c98762f76d0cca5 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Wed, 13 Dec 2023 17:39:07 +0100 Subject: [PATCH 069/131] check param value, does it change? --- src/Applications/src/gundamMarginalise.cxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 1a2279a74..3d5751a66 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -461,7 +461,8 @@ int main(int argc, char** argv){ } // change the parameter values for( int iPar = 0 ; iPar < nStripped ; iPar++ ) { - double epsilon = gRandom->Gaus(0, 0.000001); + double sigma = strippedParameterList[iPar]->getStdDevValue(); + double epsilon = gRandom->Gaus(0, sigma/10.); strippedParameterList[iPar]->setParameterValue( epsilon + getParameterValueFromTextFile(parInjectFile, strippedParameterList[iPar]->getFullTitle()) ); From fc9690ce36d91d2d09bc2f30d4d7e6bdc90266e3 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Wed, 13 Dec 2023 17:39:33 +0100 Subject: [PATCH 070/131] check param value, does it change? --- src/Applications/src/gundamMarginalise.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 3d5751a66..0b00e9a65 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -434,7 +434,7 @@ int main(int argc, char** argv){ // reset weights vector weightsChiSquare.clear(); // Do the throwing: - //propagator.getParametersManager().throwParametersFromGlobalCovariance(weightsChiSquare); + propagator.getParametersManager().throwParametersFromGlobalCovariance(weightsChiSquare); //propagator.propagateParametersOnSamples(); // Not necessary (included in updateLlhCache()) if(injectParamsManually) { From 8465bee04d7e9ed19b272f7c5408e12851f8098f Mon Sep 17 00:00:00 2001 From: lgiannes Date: Wed, 13 Dec 2023 17:48:22 +0100 Subject: [PATCH 071/131] check param value, does it change? --- src/Applications/src/gundamMarginalise.cxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 0b00e9a65..30bd6d085 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -462,7 +462,8 @@ int main(int argc, char** argv){ // change the parameter values for( int iPar = 0 ; iPar < nStripped ; iPar++ ) { double sigma = strippedParameterList[iPar]->getStdDevValue(); - double epsilon = gRandom->Gaus(0, sigma/10.); + double epsilon = gRandom->Gaus(0, sigma/50.); + LogInfo<getFullTitle()<<" e: "<setParameterValue( epsilon + getParameterValueFromTextFile(parInjectFile, strippedParameterList[iPar]->getFullTitle()) ); From 5c57711c101340582c6f88200b217b2c59179811 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Wed, 13 Dec 2023 17:52:48 +0100 Subject: [PATCH 072/131] check param value, does it change? --- src/Applications/src/gundamMarginalise.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 30bd6d085..66f4ae159 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -463,7 +463,7 @@ int main(int argc, char** argv){ for( int iPar = 0 ; iPar < nStripped ; iPar++ ) { double sigma = strippedParameterList[iPar]->getStdDevValue(); double epsilon = gRandom->Gaus(0, sigma/50.); - LogInfo<getFullTitle()<<" e: "<getFullTitle()<<" e: "<setParameterValue( epsilon + getParameterValueFromTextFile(parInjectFile, strippedParameterList[iPar]->getFullTitle()) ); From 3969d7c7fe29ec1fc74c650849772d0ed4286c5d Mon Sep 17 00:00:00 2001 From: lgiannes Date: Wed, 13 Dec 2023 18:01:30 +0100 Subject: [PATCH 073/131] check param value, does it change? --- src/Applications/src/gundamMarginalise.cxx | 21 ++++++++++++++++----- 1 file changed, 16 insertions(+), 5 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 66f4ae159..82d30d595 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -423,8 +423,8 @@ int main(int argc, char** argv){ ///////////////////////////////////// std::stringstream ss; ss << LogWarning.getPrefixString() << "Generating " << nToys << " toys..."; - - + double LLH_sum{0};// needed when injecting parameters manually + double injectedLLH{0};// needed when injecting parameters manually for( int iToy = 0 ; iToy < nToys ; iToy++ ){ @@ -463,6 +463,7 @@ int main(int argc, char** argv){ for( int iPar = 0 ; iPar < nStripped ; iPar++ ) { double sigma = strippedParameterList[iPar]->getStdDevValue(); double epsilon = gRandom->Gaus(0, sigma/50.); + if (iToy==0) epsilon = 0; //LogInfo<getFullTitle()<<" e: "<setParameterValue( epsilon + getParameterValueFromTextFile(parInjectFile, strippedParameterList[iPar]->getFullTitle()) @@ -480,10 +481,14 @@ int main(int argc, char** argv){ }// end if(injectParamsManually) - LogInfo<<"Computing LH..."<Write(); double det = 1.0; From 8acc351654b5e10c382651c465ecd033225ebe97 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Wed, 13 Dec 2023 18:10:42 +0100 Subject: [PATCH 074/131] check param value, does it change? --- src/Applications/src/gundamMarginalise.cxx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 82d30d595..749767969 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -424,7 +424,7 @@ int main(int argc, char** argv){ std::stringstream ss; ss << LogWarning.getPrefixString() << "Generating " << nToys << " toys..."; double LLH_sum{0};// needed when injecting parameters manually - double injectedLLH{0};// needed when injecting parameters manually + double injectedLLH{-1};// needed when injecting parameters manually for( int iToy = 0 ; iToy < nToys ; iToy++ ){ @@ -462,7 +462,7 @@ int main(int argc, char** argv){ // change the parameter values for( int iPar = 0 ; iPar < nStripped ; iPar++ ) { double sigma = strippedParameterList[iPar]->getStdDevValue(); - double epsilon = gRandom->Gaus(0, sigma/50.); + double epsilon = gRandom->Gaus(0, sigma/100.); if (iToy==0) epsilon = 0; //LogInfo<getFullTitle()<<" e: "<setParameterValue( @@ -485,7 +485,7 @@ int main(int argc, char** argv){ propagator.updateLlhCache(); LLH = propagator.getLlhBuffer(); LLH_sum += LLH; - if(iToy==0){ + if(iToy==0 and injectParamsManually){ injectedLLH = LLH; } //LogInfo<<"LLH: "< Date: Wed, 13 Dec 2023 18:16:36 +0100 Subject: [PATCH 075/131] check param value, does it change? --- src/Applications/src/gundamMarginalise.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 749767969..b0f972494 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -462,7 +462,7 @@ int main(int argc, char** argv){ // change the parameter values for( int iPar = 0 ; iPar < nStripped ; iPar++ ) { double sigma = strippedParameterList[iPar]->getStdDevValue(); - double epsilon = gRandom->Gaus(0, sigma/100.); + double epsilon = gRandom->Gaus(0, sigma/200.); if (iToy==0) epsilon = 0; //LogInfo<getFullTitle()<<" e: "<setParameterValue( From 9d9dc7d4b9d393f75fadeb4fa42ad6958f5e07d7 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Wed, 13 Dec 2023 18:57:23 +0100 Subject: [PATCH 076/131] check param value, does it change? --- src/Applications/src/gundamMarginalise.cxx | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index b0f972494..714e48294 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -425,6 +425,9 @@ int main(int argc, char** argv){ double LLH_sum{0};// needed when injecting parameters manually double injectedLLH{-1};// needed when injecting parameters manually + if (injectParamsManually){ + LogInfo<< "Injecting parameters from file: " << parInjectFile << std::endl; + } for( int iToy = 0 ; iToy < nToys ; iToy++ ){ @@ -438,7 +441,6 @@ int main(int argc, char** argv){ //propagator.propagateParametersOnSamples(); // Not necessary (included in updateLlhCache()) if(injectParamsManually) { - LogInfo<< "Injecting parameters from file: " << parInjectFile << std::endl; // count the number of parameters to be injected int nStripped{0}; From 67060f4c4e47302e47a7f1857775ceb5e91f31f4 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Wed, 13 Dec 2023 19:26:41 +0100 Subject: [PATCH 077/131] check param value, does it change? --- src/Applications/src/gundamMarginalise.cxx | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 714e48294..ef0d94eb5 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -462,16 +462,20 @@ int main(int argc, char** argv){ } } // change the parameter values + double epsilonNorm=0; for( int iPar = 0 ; iPar < nStripped ; iPar++ ) { double sigma = strippedParameterList[iPar]->getStdDevValue(); - double epsilon = gRandom->Gaus(0, sigma/200.); + double epsilon = gRandom->Gaus(0, sigma/nStripped); + epsilonNorm += epsilon*epsilon; if (iToy==0) epsilon = 0; //LogInfo<getFullTitle()<<" e: "<setParameterValue( epsilon + getParameterValueFromTextFile(parInjectFile, strippedParameterList[iPar]->getFullTitle()) ); - } + epsilonNorm = sqrt(epsilonNorm); + LogInfo<<"epsilonNorm: "< Date: Wed, 13 Dec 2023 19:42:52 +0100 Subject: [PATCH 078/131] check param value, does it change? --- src/Applications/src/gundamMarginalise.cxx | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index ef0d94eb5..0bbd45329 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -50,12 +50,8 @@ int main(int argc, char** argv){ // I need a config file where the list of parameters to marginalise over are specified (look at the XSec config file to get inspired) clParser.addOption("configFile", {"-c"}, "Specify the parameters to marginalise over"); - // (I think) I need the output file from a fitter to use as input here clParser.addOption("fitterOutputFile", {"-f"}, "Specify the fitter output file"); - // Think carefully what do you need to put in the output file - // 1. Marginalised covariance matrix: gotten from teh outoutFitter file - // what else? clParser.addOption("outputFile", {"-o", "--out-file"}, "Specify the Marginaliser output file"); clParser.addOption("nbThreads", {"-t", "--nb-threads"}, "Specify nb of parallel threads"); @@ -234,6 +230,7 @@ int main(int argc, char** argv){ {"fitterOutputFile", "Fit_%s"}, {"nToys", "nToys_%s"}, {"randomSeed", "Seed_%s"}, + {"parInject", "parInject_%s" }; outFilePath = GundamUtils::generateFileName(clParser, appendixDict) + ".root"; From 7ae5d2ee8dd86bd0c6a15716b8bb90c373af7581 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Wed, 13 Dec 2023 19:44:13 +0100 Subject: [PATCH 079/131] check param value, does it change? --- src/Applications/src/gundamMarginalise.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 0bbd45329..2fed4d5a3 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -230,7 +230,7 @@ int main(int argc, char** argv){ {"fitterOutputFile", "Fit_%s"}, {"nToys", "nToys_%s"}, {"randomSeed", "Seed_%s"}, - {"parInject", "parInject_%s" + {"parInject", "parInject_%s"} }; outFilePath = GundamUtils::generateFileName(clParser, appendixDict) + ".root"; From d25f9fe3d5795f1871a8b9306d7bb6b5184391e8 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Wed, 13 Dec 2023 19:49:19 +0100 Subject: [PATCH 080/131] check param value, does it change? --- src/Applications/src/gundamMarginalise.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 2fed4d5a3..0e2445b20 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -462,7 +462,7 @@ int main(int argc, char** argv){ double epsilonNorm=0; for( int iPar = 0 ; iPar < nStripped ; iPar++ ) { double sigma = strippedParameterList[iPar]->getStdDevValue(); - double epsilon = gRandom->Gaus(0, sigma/nStripped); + double epsilon = gRandom->Gaus(0, sigma/(nStripped/10.)); epsilonNorm += epsilon*epsilon; if (iToy==0) epsilon = 0; //LogInfo<getFullTitle()<<" e: "< Date: Wed, 13 Dec 2023 19:58:12 +0100 Subject: [PATCH 081/131] check param value, does it change? --- src/Applications/src/gundamMarginalise.cxx | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 0e2445b20..c49562dc7 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -422,6 +422,7 @@ int main(int argc, char** argv){ double LLH_sum{0};// needed when injecting parameters manually double injectedLLH{-1};// needed when injecting parameters manually + double epsilonNormAverage{0};// needed when injecting parameters manually if (injectParamsManually){ LogInfo<< "Injecting parameters from file: " << parInjectFile << std::endl; } @@ -462,7 +463,7 @@ int main(int argc, char** argv){ double epsilonNorm=0; for( int iPar = 0 ; iPar < nStripped ; iPar++ ) { double sigma = strippedParameterList[iPar]->getStdDevValue(); - double epsilon = gRandom->Gaus(0, sigma/(nStripped/10.)); + double epsilon = gRandom->Gaus(0, sigma/(nStripped)); epsilonNorm += epsilon*epsilon; if (iToy==0) epsilon = 0; //LogInfo<getFullTitle()<<" e: "<Write(); From 4ccbcb986b28d11f95606eae424665b8c80d541d Mon Sep 17 00:00:00 2001 From: lgiannes Date: Mon, 18 Dec 2023 13:27:27 +0100 Subject: [PATCH 082/131] another override implementation of throwCorrelatedParameters --- src/Applications/src/gundamMarginalise.cxx | 13 +------------ 1 file changed, 1 insertion(+), 12 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index c49562dc7..907cc364f 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -534,18 +534,7 @@ int main(int argc, char** argv){ // } - // for now set it to false. Still need to understand/implement this - bool enableStatThrowInToys{false}; -// if( enableStatThrowInToys ){ -// for( auto& xsec : crossSectionDataList ){ -// if( enableEventMcThrow ){ -// // Take into account the finite amount of event in MC -// xsec.samplePtr->getMcContainer().throwEventMcError(); -// } -// // Asimov bin content -> toy data -// xsec.samplePtr->getMcContainer().throwStatError(); -// } -// } + // i don't know what is this for, so I comment it for now // writeBinDataFct(); From 671505f0dadcd09e6927d0f9d1ea774c4b025d31 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Mon, 18 Dec 2023 14:19:39 +0100 Subject: [PATCH 083/131] update toolbox --- submodules/cpp-generic-toolbox | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/submodules/cpp-generic-toolbox b/submodules/cpp-generic-toolbox index fa7710379..6a671d02c 160000 --- a/submodules/cpp-generic-toolbox +++ b/submodules/cpp-generic-toolbox @@ -1 +1 @@ -Subproject commit fa7710379eb8a2c8f079c13c768c01e677493c88 +Subproject commit 6a671d02c610af65f28f3105b0a349eb311d0599 From edd0119d844b1223f6f78a9c1f840c76e7862bc9 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Tue, 19 Dec 2023 12:28:58 +0100 Subject: [PATCH 084/131] implemented pedestal+gaus sampling distribution --- src/Applications/src/gundamMarginalise.cxx | 69 +++++++++++-------- .../include/ParametersManager.h | 1 + .../src/ParametersManager.cpp | 17 +++-- 3 files changed, 55 insertions(+), 32 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 907cc364f..642fad651 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -59,6 +59,7 @@ int main(int argc, char** argv){ clParser.addOption("randomSeed", {"-s", "--seed"}, "Set random seed"); clParser.addTriggerOption("usingGpu", {"--gpu"}, "Use GPU parallelization"); clParser.addOption("parInject", {"--parameters-inject"}, "Input txt file for injecting params"); + clParser.addOption("pedestal", {"--pedestal"}, "Add pedestal to the sampling distribution (percents)"); clParser.addDummyOption("Trigger options:"); clParser.addTriggerOption("dryRun", {"-d", "--dry-run"}, "Only overrides fitter config and print it."); @@ -165,6 +166,16 @@ int main(int argc, char** argv){ }else{ LogInfo << "Injecting best fit parameters as prior" << std::endl; } + double pedestalEntity = 0; + bool usePedestal = false; + double pedestalLeftEdge = -5, pedestalRightEdge = 5; + if(clParser.isOptionTriggered("pedestal")){ + usePedestal = true; + pedestalEntity = clParser.getOptionVal("pedestal")/100.; + LogInfo << "Using Gauss+pedestal sampling: "<< pedestalEntity*100 << "% of the throws will be drawn from a uniform distribution."<( cHandler.getConfig(), "fitterEngineConfig/propagatorConfig" ); @@ -230,7 +241,8 @@ int main(int argc, char** argv){ {"fitterOutputFile", "Fit_%s"}, {"nToys", "nToys_%s"}, {"randomSeed", "Seed_%s"}, - {"parInject", "parInject_%s"} + {"parInject", "parInject_%s"}, + {"pedestal", "ped_%s_percent"} }; outFilePath = GundamUtils::generateFileName(clParser, appendixDict) + ".root"; @@ -299,7 +311,8 @@ int main(int argc, char** argv){ parState_->SetName("postFitParameters_TNamed"); parState_->Write(); }); - TVectorD bestFitParameters (parametersBestFit.size(), parametersBestFit.data()); + int nParameters = parametersBestFit.size(); + TVectorD bestFitParameters (nParameters, parametersBestFit.data()); bestFitParameters.Write("bestFitParameters_TVectorD"); app.getOutfilePtr()->WriteObject(¶meterFullTitles, "parameterFullTitles"); app.getOutfilePtr()->cd(); @@ -323,6 +336,7 @@ int main(int argc, char** argv){ std::vector dialResponse; // weightsChiSquare.reserve(1000); double LLH, gLLH, priorSum, LLHwrtBestFit; + double totalWeight; // use in the pedestal case, just as a placeholder margThrowTree->Branch("Parameters", ¶meters); margThrowTree->Branch("Marginalise", &margThis); margThrowTree->Branch("prior", &prior); @@ -415,11 +429,7 @@ int main(int argc, char** argv){ // << par.getMaxMirror() <<" --- marg? "< 10){ -// LogInfo<<"WARNING: BIG THROW!!"< 10){ + LogInfo<<"WARNING: BIG THROW!!"<getParameterList().at(0).getParameterValue(); -// // } - - - - // i don't know what is this for, so I comment it for now -// writeBinDataFct(); - // Write the branches margThrowTree->Fill(); - //debug // if(iToy==0 || iToy==1){ // LogInfo<<"DEBUG: Toy "< &weightsChiSquare); + void throwParametersFromGlobalCovariance(std::vector &weightsChiSquare, double pedestalEntity, double pedestalLeftEdge, double pedestalRightEdge); ParameterSet* getFitParameterSetPtr(const std::string& name_); // Logger related diff --git a/src/ParametersManager/src/ParametersManager.cpp b/src/ParametersManager/src/ParametersManager.cpp index 8af4f28a5..576ee2f79 100644 --- a/src/ParametersManager/src/ParametersManager.cpp +++ b/src/ParametersManager/src/ParametersManager.cpp @@ -280,6 +280,12 @@ void ParametersManager::throwParametersFromGlobalCovariance(bool quietVerbose_){ } void ParametersManager::throwParametersFromGlobalCovariance(std::vector &weightsChiSquare){ + throwParametersFromGlobalCovariance(weightsChiSquare,0,0,0); +}// end of function + +void ParametersManager::throwParametersFromGlobalCovariance(std::vector &weightsChiSquare, + double pedestalEntity, double pedestalLeftEdge, double pedestalRightEdge + ){ // check that weightsChiSquare is an empty vector LogThrowIf( not weightsChiSquare.empty(), "ERROR: argument weightsChiSquare is not empty" ); @@ -331,7 +337,12 @@ void ParametersManager::throwParametersFromGlobalCovariance(std::vector // throwNb++; bool rethrow{false}; std::vector throws,weights; - GenericToolbox::throwCorrelatedParameters(_choleskyMatrix_.get(),throws, weights); + if(pedestalEntity==0){ + GenericToolbox::throwCorrelatedParameters(_choleskyMatrix_.get(),throws, weights); + }else{ + GenericToolbox::throwCorrelatedParameters(_choleskyMatrix_.get(),throws, weights, + pedestalEntity,pedestalLeftEdge,pedestalRightEdge); + } if(throws.size() != weights.size()){ LogInfo<<"WARNING: throws.size() != weights.size() "<< throws.size()< // reached this point: all parameters are within bounds keepThrowing = false; } +}// end of function - - -} void ParametersManager::injectParameterValues(const nlohmann::json &config_) { LogWarning << "Injecting parameters..." << std::endl; From 231ef4ca4a38ada69474c503c5e83b2a8122ad2e Mon Sep 17 00:00:00 2001 From: lgiannes Date: Tue, 19 Dec 2023 12:40:32 +0100 Subject: [PATCH 085/131] -- --- submodules/cpp-generic-toolbox | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/submodules/cpp-generic-toolbox b/submodules/cpp-generic-toolbox index 6a671d02c..3f01f2f91 160000 --- a/submodules/cpp-generic-toolbox +++ b/submodules/cpp-generic-toolbox @@ -1 +1 @@ -Subproject commit 6a671d02c610af65f28f3105b0a349eb311d0599 +Subproject commit 3f01f2f919492930073fbd4bdc52101303f64a10 From b253bdff0c12ec5ee9be472f15b05dabacf0cbe1 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Fri, 22 Dec 2023 12:27:38 +0100 Subject: [PATCH 086/131] save prior values in the output file --- src/Applications/src/gundamMarginalise.cxx | 63 ++++++---------------- 1 file changed, 17 insertions(+), 46 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 642fad651..f6a6546c9 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -195,19 +195,11 @@ int main(int argc, char** argv){ propagator.initialize(); propagator.getParametersManager().getParameterSetsList(); - for( auto& parSet : propagator.getParametersManager().getParameterSetsList() ){ - if( not parSet.isEnabled() ){ continue; } - LogInfo < "< parametersBestFit; - std::vector parameterFullTitles; + std::vector parameterNames; + std::vector bestFitValues; + std::vector priorValues; + std::vector priorSigmas; // Load post-fit parameters as "prior" so we can reset the weight to this point when throwing toys // also save the values in a vector so we can use them to compute the LLH at the best fit point ObjectReader::readObject( fitterFile.get(), "FitterEngine/postFit/parState_TNamed", [&](TNamed* parState_){ @@ -218,8 +210,11 @@ int main(int argc, char** argv){ for( auto& par : parSet.getParameterList() ){ if( not par.isEnabled() ){ continue; } // LogInfo<<" "< "<(sample.getConfig(), "parameterSetName"); -// -// // Looking for parSet -// auto foundDialCollection = std::find_if( -// propagator.getDialCollections().begin(), propagator.getDialCollections().end(), -// [&](const DialCollection& dialCollection_){ -// auto* parSetPtr{dialCollection_.getSupervisedParameterSet()}; -// if( parSetPtr == nullptr ){ return false; } -// return ( parSetPtr->getName() == associatedParSet ); -// }); -// LogThrowIf( -// foundDialCollection == propagator.getDialCollections().end(), -// "Could not find " << associatedParSet << " among fit dial collections: " -// << GenericToolbox::iterableToString(propagator.getDialCollections(), [](const DialCollection& dialCollection_){ -// return dialCollection_.getTitle(); -// } -// )); -// -// LogThrowIf(foundDialCollection->getDialBinSet().isEmpty(), "Could not find binning"); -// sample.setBinningFilePath( foundDialCollection->getDialBinSet().getFilePath() ); -// } - - - - app.setCmdLinePtr( &clParser ); app.setConfigString( ConfigUtils::ConfigHandler{margConfig}.toString() ); @@ -311,13 +279,16 @@ int main(int argc, char** argv){ parState_->SetName("postFitParameters_TNamed"); parState_->Write(); }); - int nParameters = parametersBestFit.size(); - TVectorD bestFitParameters (nParameters, parametersBestFit.data()); - bestFitParameters.Write("bestFitParameters_TVectorD"); - app.getOutfilePtr()->WriteObject(¶meterFullTitles, "parameterFullTitles"); + unsigned int nParameters = parameterNames.size(); + TVectorD bestFitParameters_TVectorD(nParameters,bestFitValues.data()); + TVectorD priorParameters_TVectorD(nParameters,priorValues.data()); + TVectorD priorSigmas_TVectorD(nParameters,priorSigmas.data()); + app.getOutfilePtr()->WriteObject(¶meterNames, "parameterFullTitles"); + app.getOutfilePtr()->WriteObject(&bestFitParameters_TVectorD, "bestFitParameters_TVectorD"); + app.getOutfilePtr()->WriteObject(&priorParameters_TVectorD, "priorParameters_TVectorD"); + app.getOutfilePtr()->WriteObject(&priorSigmas_TVectorD, "priorSigmas_TVectorD"); app.getOutfilePtr()->cd(); - auto* marginalisationDir{ GenericToolbox::mkdirTFile(app.getOutfilePtr(), "marginalisation") }; LogInfo << "Creating throws tree" << std::endl; auto* margThrowTree = new TTree("margThrow", "margThrow"); From 9ff9372dd53b75bfc2bea1664c1d7c1cddd2067a Mon Sep 17 00:00:00 2001 From: lgiannes Date: Fri, 5 Jan 2024 12:49:56 +0100 Subject: [PATCH 087/131] debug for pedestal --- src/Applications/src/gundamMarginalise.cxx | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index f6a6546c9..06d4723b8 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -274,7 +274,6 @@ int main(int argc, char** argv){ hCovPostFit_->Write();// save the TH2D cov. matrix also in the output file }); // write the best fit parameters vector to the output file - // TODO: implement using TH1D: (FitterEngine/postFit/Migrad/errors//values/postFitErrors_TH1D) ObjectReader::readObject( fitterFile.get(), "FitterEngine/postFit/parState_TNamed", [&](TNamed* parState_){ parState_->SetName("postFitParameters_TNamed"); parState_->Write(); @@ -479,7 +478,6 @@ int main(int argc, char** argv){ // make the LH a probability distribution (but still work with the log) LLH /= 2.0; //LLH -= bestFitLLH/2.0; - // TODO: figure out normalization factor here. //LLH -= nParameters*log(1.0/2.0*TMath::Pi()); LLH_sum += LLH; From 80cdbf5549c75cb2542df403a5ee2157f9a04a0c Mon Sep 17 00:00:00 2001 From: lgiannes Date: Fri, 5 Jan 2024 17:07:37 +0100 Subject: [PATCH 088/131] debug for pedestal (see if it changes) --- src/Applications/src/gundamMarginalise.cxx | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 06d4723b8..f2f2dd0a6 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -168,7 +168,7 @@ int main(int argc, char** argv){ } double pedestalEntity = 0; bool usePedestal = false; - double pedestalLeftEdge = -5, pedestalRightEdge = 5; + double pedestalLeftEdge = -3, pedestalRightEdge = -pedestalLeftEdge; if(clParser.isOptionTriggered("pedestal")){ usePedestal = true; pedestalEntity = clParser.getOptionVal("pedestal")/100.; @@ -512,10 +512,9 @@ int main(int argc, char** argv){ } } - LogInfo<<"LLH: "< 10){ - LogInfo<<"WARNING: BIG THROW!!"< Date: Wed, 24 Jan 2024 16:11:19 +0100 Subject: [PATCH 089/131] no change --- src/Applications/src/gundamMarginalise.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index f2f2dd0a6..7fee8812e 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -512,7 +512,7 @@ int main(int argc, char** argv){ } } - LogInfo<<"LLH: "< Date: Thu, 25 Jan 2024 17:14:28 +0100 Subject: [PATCH 090/131] implemented TTree with throws in the PTheta format. [NEED TO CHECK SORTING] --- src/Applications/src/gundamMarginalise.cxx | 88 +++++++++++++++------- 1 file changed, 62 insertions(+), 26 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 7fee8812e..3609a5cb1 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -60,10 +60,12 @@ int main(int argc, char** argv){ clParser.addTriggerOption("usingGpu", {"--gpu"}, "Use GPU parallelization"); clParser.addOption("parInject", {"--parameters-inject"}, "Input txt file for injecting params"); clParser.addOption("pedestal", {"--pedestal"}, "Add pedestal to the sampling distribution (percents)"); + clParser.addOption("weightCap", {"--weight-cap"}, "Cap the weight of the throws (default: 1.e8)", 1); clParser.addDummyOption("Trigger options:"); clParser.addTriggerOption("dryRun", {"-d", "--dry-run"}, "Only overrides fitter config and print it."); + // I probably don't need this //clParser.addTriggerOption("useBfAsXsec", {"--use-bf-as-xsec"}, "Use best-fit as x-sec value instead of mean of toys."); @@ -176,6 +178,13 @@ int main(int argc, char** argv){ LogInfo <<"Using default pedestal interval: ["<("weightCap"); + LogInfo << "Using weight cap: "<< weightCap << std::endl; + }else{ + LogInfo << "Using default weight cap: "<< weightCap << std::endl; + } auto configPropagator = GenericToolbox::Json::fetchValuePath( cHandler.getConfig(), "fitterEngineConfig/propagatorConfig" ); @@ -292,6 +301,10 @@ int main(int argc, char** argv){ LogInfo << "Creating throws tree" << std::endl; auto* margThrowTree = new TTree("margThrow", "margThrow"); margThrowTree->SetDirectory( GenericToolbox::mkdirTFile(marginalisationDir, "throws") ); // temp saves will be done here + auto* ThrowsPThetaFormat = new TTree("PThetaThrows", "PThetaThrows"); + + + // make a TTree with the following branches, to be filled with the throws // 1. marginalised parameters drew: vector // 2. non-marginalised parameters drew: vector @@ -299,6 +312,7 @@ int main(int argc, char** argv){ // 4. "g" value: the chi square value as extracted from the covariance matrix: double // 5. value of the priors for the marginalised parameters: vector + // branches for margThrowTree std::vector parameters; std::vector margThis; std::vector prior; @@ -317,6 +331,12 @@ int main(int argc, char** argv){ margThrowTree->Branch("weightsChiSquare", &weightsChiSquare); margThrowTree->Branch("dialResponse", &dialResponse); + // branches for ThrowsPThetaFormat + std::vector survivingParameterValues; + double LhOverGauss; + ThrowsPThetaFormat->Branch("ParValues", &survivingParameterValues); + ThrowsPThetaFormat->Branch("weight", &LhOverGauss); + int nToys{ clParser.getOptionVal("nToys") }; @@ -418,10 +438,12 @@ int main(int argc, char** argv){ // reset weights vector weightsChiSquare.clear(); - // Do the throwing: + // Do the throwing if (usePedestal){ + // if the pedestal option is enabled, a uniform distribution is added to the gaussian sampling distribution propagator.getParametersManager().throwParametersFromGlobalCovariance(weightsChiSquare, pedestalEntity, pedestalLeftEdge, pedestalRightEdge); }else { + // standard case: throw according to the covariance matrix propagator.getParametersManager().throwParametersFromGlobalCovariance(weightsChiSquare); } @@ -476,9 +498,10 @@ int main(int argc, char** argv){ propagator.updateLlhCache(); LLH = propagator.getLlhBuffer(); // make the LH a probability distribution (but still work with the log) - LLH /= 2.0; - //LLH -= bestFitLLH/2.0; - //LLH -= nParameters*log(1.0/2.0*TMath::Pi()); + // This is an approximation, it works only in case of gaussian LH + LLH /= -2.0; + LLH += nParameters/2.*log(1.0/(2.0*TMath::Pi())); + LLH = -LLH; LLH_sum += LLH; if(iToy==0 and injectParamsManually){ @@ -494,7 +517,6 @@ int main(int argc, char** argv){ parameters.clear(); margThis.clear(); - prior.clear(); int iPar=0; for( auto& parSet : propagator.getParametersManager().getParameterSetsList() ) { if (not parSet.isEnabled()) { continue; } @@ -511,18 +533,31 @@ int main(int argc, char** argv){ iPar++; } } + //debug + LogInfo<<"LLH: "<getParameterList().at(0).getParameterValue(); -// } + if ( LLH-gLLH > log(weightCap)){ + LogInfo<<"Throw "<Write(); - - double det = 1.0; - TMatrixD eigenVectors = (*propagator.getParametersManager().getGlobalCovarianceMatrix()); - TVectorD eigenValues(parameters.size()); - eigenVectors.EigenVectors(eigenValues); - //LogInfo<<"Eigenvalues: "<Write(); + + // Compute the determinant of the covariance matrix +// double det = 1.0; +// TMatrixD eigenVectors = (*propagator.getParametersManager().getGlobalCovarianceMatrix()); +// TVectorD eigenValues(parameters.size()); +// eigenVectors.EigenVectors(eigenValues); +// //LogInfo<<"Eigenvalues: "< Date: Mon, 29 Jan 2024 23:33:06 +0100 Subject: [PATCH 091/131] (maybe) fixed a bug: uncleared prior --- src/Applications/src/gundamMarginalise.cxx | 1 + 1 file changed, 1 insertion(+) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 3609a5cb1..a74d824a7 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -517,6 +517,7 @@ int main(int argc, char** argv){ parameters.clear(); margThis.clear(); + prior.clear(); int iPar=0; for( auto& parSet : propagator.getParametersManager().getParameterSetsList() ) { if (not parSet.isEnabled()) { continue; } From 876bb87c90644bf9df732a52b985666ac7f137dd Mon Sep 17 00:00:00 2001 From: lgiannes Date: Tue, 30 Jan 2024 09:52:33 +0100 Subject: [PATCH 092/131] removed re-throw: too much time-consuming --- src/Applications/src/gundamMarginalise.cxx | 31 +++++++++------------- 1 file changed, 13 insertions(+), 18 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index a74d824a7..dfbb0b611 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -179,6 +179,7 @@ int main(int argc, char** argv){ LogInfo <<". User defined interval to be implemented"<("weightCap"); LogInfo << "Using weight cap: "<< weightCap << std::endl; @@ -518,6 +519,7 @@ int main(int argc, char** argv){ parameters.clear(); margThis.clear(); prior.clear(); + survivingParameterValues.clear(); int iPar=0; for( auto& parSet : propagator.getParametersManager().getParameterSetsList() ) { if (not parSet.isEnabled()) { continue; } @@ -530,31 +532,21 @@ int main(int argc, char** argv){ prior.push_back(par.getDistanceFromNominal() * par.getDistanceFromNominal()); priorSum += prior.back(); gLLH += weightsChiSquare[iPar]; + survivingParameterValues.push_back(par.getParameterValue()); + //LogInfo< log(weightCap)) { + LogInfo << "Throw " << iToy << " rejected: LLH-gLLH = " << LLH - gLLH << std::endl; + countBigThrows++; + } + //debug LogInfo<<"LLH: "< log(weightCap)){ - LogInfo<<"Throw "< log(weightCap): "< Date: Sun, 4 Feb 2024 16:27:02 +0100 Subject: [PATCH 093/131] added marg_param_list --- src/Applications/src/gundamMarginalise.cxx | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index dfbb0b611..60ec94712 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -25,6 +25,7 @@ #include #include +#include LoggerInit([]{ @@ -349,6 +350,9 @@ int main(int argc, char** argv){ LogInfo<<"Marginalised parameters: "<Add(new TObjString(par.getFullTitle().c_str())); + } if(par.isMarginalised()){ LogInfo << "Parameter " << par.getFullTitle() << " -> type: " << par.getPriorType() << " mu=" << par.getPriorValue() @@ -413,6 +420,8 @@ int main(int argc, char** argv){ } } + // write the list of parameters that will not be marginalised to the output file + app.getOutfilePtr()->WriteObject(marg_param_list, "marg_param_list"); // LogInfo << par.getFullTitle() << " -> type: " << par.getPriorType() << " mu=" << par.getPriorValue() // << " sigma= " << par.getStdDevValue() << " limits: " << par.getMinValue() << " - " // << par.getMaxValue() << " limits (phys): " << par.getMinPhysical() << " - " @@ -532,7 +541,8 @@ int main(int argc, char** argv){ prior.push_back(par.getDistanceFromNominal() * par.getDistanceFromNominal()); priorSum += prior.back(); gLLH += weightsChiSquare[iPar]; - survivingParameterValues.push_back(par.getParameterValue()); + if(!par.isMarginalised()) + survivingParameterValues.push_back(par.getParameterValue()); //LogInfo< Date: Sun, 4 Feb 2024 16:35:08 +0100 Subject: [PATCH 094/131] added cmakelists --- CMakeLists.txt | 54 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 54 insertions(+) create mode 100644 CMakeLists.txt diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 000000000..c132db471 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,54 @@ +# +# GUNDAM project CMake main file +# + +# cmake_minimum_required() should be called prior to this top-level project() +# 3.5 minimum a priori. Taking a lower version as min will make recent CMake +# version complain about the deprecation version older than 3.5. +# Might require higher version for specific features. +cmake_minimum_required( VERSION 3.5 FATAL_ERROR ) + +# Define project +project( GUNDAM ) + + +# colored messages +include( ${CMAKE_SOURCE_DIR}/cmake/utils/cmessage.cmake ) + +# define cmake options +include( ${CMAKE_SOURCE_DIR}/cmake/options.cmake ) + +# git version checking +include( ${CMAKE_SOURCE_DIR}/cmake/version.cmake ) + +# checking dependencies +include( ${CMAKE_SOURCE_DIR}/cmake/dependencies.cmake ) + +# setup compiler options +include( ${CMAKE_SOURCE_DIR}/cmake/compiler.cmake ) + +# setting up submodule libraries +include( ${CMAKE_SOURCE_DIR}/cmake/submodules.cmake ) + +# defining GUNDAM modules +include( ${CMAKE_SOURCE_DIR}/cmake/modules.cmake ) + +# defining GUNDAM applications +include( ${CMAKE_SOURCE_DIR}/cmake/applications.cmake ) + +# defining GUNDAM tests +include( ${CMAKE_SOURCE_DIR}/cmake/tests.cmake ) + + +message("") +cmessage( WARNING "Identified GUNDAM version: ${GUNDAM_FULL_VERSION_STR}" ) + +# A command that works to compile gundam from any sub-directory. You +# can get this aliased to "gundam-build" by sourcing +# +# . "$(git rev-parse --show-toplevel)/cmake/scripts/gundam-setup.sh" +# +# That runs: +# Local Variables: +# compile-command:"$(git rev-parse --show-toplevel)/cmake/gundam-build.sh" +# End: From f128f85fcd14a0caac8cec99698b57b5ae43e126 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Mon, 5 Feb 2024 16:17:56 +0100 Subject: [PATCH 095/131] hide error in JointProbability --- src/SamplesManager/src/JointProbability.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/SamplesManager/src/JointProbability.cpp b/src/SamplesManager/src/JointProbability.cpp index 9ae731beb..f8c83d753 100644 --- a/src/SamplesManager/src/JointProbability.cpp +++ b/src/SamplesManager/src/JointProbability.cpp @@ -211,7 +211,7 @@ namespace JointProbability{ << GET_VAR_NAME_VALUE(nomMC->GetBinError(bin_)) << std::endl << GET_VAR_NAME_VALUE(predMC->GetBinError(bin_)) << std::endl << GET_VAR_NAME_VALUE(predMC->GetBinContent(bin_)) << std::endl; - LogThrow("Bad chisq for bin"); + //LogThrow("Bad chisq for bin"); } if(verboseLevel>=3){ From 74b71420fbaa679384b353eab5a886fed455893d Mon Sep 17 00:00:00 2001 From: lgiannes Date: Mon, 5 Feb 2024 16:41:28 +0100 Subject: [PATCH 096/131] restored JointProbability.cpp --- src/SamplesManager/src/JointProbability.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/SamplesManager/src/JointProbability.cpp b/src/SamplesManager/src/JointProbability.cpp index f8c83d753..9ae731beb 100644 --- a/src/SamplesManager/src/JointProbability.cpp +++ b/src/SamplesManager/src/JointProbability.cpp @@ -211,7 +211,7 @@ namespace JointProbability{ << GET_VAR_NAME_VALUE(nomMC->GetBinError(bin_)) << std::endl << GET_VAR_NAME_VALUE(predMC->GetBinError(bin_)) << std::endl << GET_VAR_NAME_VALUE(predMC->GetBinContent(bin_)) << std::endl; - //LogThrow("Bad chisq for bin"); + LogThrow("Bad chisq for bin"); } if(verboseLevel>=3){ From 12977a2188c98d0b863ad0cbdb8845c74c55ce7e Mon Sep 17 00:00:00 2001 From: lgiannes Date: Tue, 6 Feb 2024 16:10:18 +0100 Subject: [PATCH 097/131] commented "LogThrow" in JointProbability.cpp --- src/SamplesManager/src/JointProbability.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/SamplesManager/src/JointProbability.cpp b/src/SamplesManager/src/JointProbability.cpp index 9ae731beb..17919238c 100644 --- a/src/SamplesManager/src/JointProbability.cpp +++ b/src/SamplesManager/src/JointProbability.cpp @@ -211,7 +211,7 @@ namespace JointProbability{ << GET_VAR_NAME_VALUE(nomMC->GetBinError(bin_)) << std::endl << GET_VAR_NAME_VALUE(predMC->GetBinError(bin_)) << std::endl << GET_VAR_NAME_VALUE(predMC->GetBinContent(bin_)) << std::endl; - LogThrow("Bad chisq for bin"); + // LogThrow("Bad chisq for bin"); } if(verboseLevel>=3){ From e3d0436d2fc4f1786a8a8fd486bf13887ddcb1eb Mon Sep 17 00:00:00 2001 From: lgiannes Date: Tue, 6 Feb 2024 18:49:26 +0100 Subject: [PATCH 098/131] changed throw in error in SampleSet.cpp --- src/SamplesManager/src/SampleSet.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/SamplesManager/src/SampleSet.cpp b/src/SamplesManager/src/SampleSet.cpp index d37ed0633..c6c6c140e 100644 --- a/src/SamplesManager/src/SampleSet.cpp +++ b/src/SamplesManager/src/SampleSet.cpp @@ -104,7 +104,9 @@ double SampleSet::evalLikelihood(){ double llh = 0.; for( auto& sample : _fitSampleList_ ){ llh += this->evalLikelihood(sample); - LogThrowIf(std::isnan(llh) or std::isinf(llh), sample.getName() << ": reportde likelihood is invalid:" << llh); + if(std::isnan(llh) or std::isinf(llh)) + LogError(sample.getName() << ": reported likelihood is invalid:" << llh); +// LogThrowIf(std::isnan(llh) or std::isinf(llh), sample.getName() << ": reported likelihood is invalid:" << llh); } return llh; } From 097254eac3c465c44dc299953e37bc6579198de2 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Tue, 6 Feb 2024 18:52:10 +0100 Subject: [PATCH 099/131] debug in SampleSet.cpp --- src/SamplesManager/src/SampleSet.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/SamplesManager/src/SampleSet.cpp b/src/SamplesManager/src/SampleSet.cpp index c6c6c140e..3a3d5c0f7 100644 --- a/src/SamplesManager/src/SampleSet.cpp +++ b/src/SamplesManager/src/SampleSet.cpp @@ -104,8 +104,7 @@ double SampleSet::evalLikelihood(){ double llh = 0.; for( auto& sample : _fitSampleList_ ){ llh += this->evalLikelihood(sample); - if(std::isnan(llh) or std::isinf(llh)) - LogError(sample.getName() << ": reported likelihood is invalid:" << llh); + LogErrorIf(std::isnan(llh) or std::isinf(llh)) << sample.getName() << ": reported likelihood is invalid:" << llh; // LogThrowIf(std::isnan(llh) or std::isinf(llh), sample.getName() << ": reported likelihood is invalid:" << llh); } return llh; From d91a21580c578baadd19bb8d59228984261c4fc6 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Wed, 7 Feb 2024 13:35:28 +0100 Subject: [PATCH 100/131] SUPER verbose version HANDLE WITH CARE (Jointprobability.cpp) --- src/SamplesManager/src/JointProbability.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/SamplesManager/src/JointProbability.cpp b/src/SamplesManager/src/JointProbability.cpp index 17919238c..c9d43f758 100644 --- a/src/SamplesManager/src/JointProbability.cpp +++ b/src/SamplesManager/src/JointProbability.cpp @@ -199,7 +199,7 @@ namespace JointProbability{ } if (not std::isfinite(chisq)) [[unlikely]] { - LogError << "Infinite chi2: " << chisq << std::endl + LogInfo << "Infinite chi2: " << chisq << std::endl << " bin " << bin_ << GET_VAR_NAME_VALUE(predVal) << std::endl << GET_VAR_NAME_VALUE(dataVal) << std::endl From 339910738d870ee2f287ffa8a4a9d7697e9a4891 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Wed, 7 Feb 2024 16:09:41 +0100 Subject: [PATCH 101/131] SUPER verbose version HANDLE WITH CARE (Jointprobability.cpp) --- src/SamplesManager/src/JointProbability.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/SamplesManager/src/JointProbability.cpp b/src/SamplesManager/src/JointProbability.cpp index c9d43f758..78053325f 100644 --- a/src/SamplesManager/src/JointProbability.cpp +++ b/src/SamplesManager/src/JointProbability.cpp @@ -198,7 +198,7 @@ namespace JointProbability{ } } - if (not std::isfinite(chisq)) [[unlikely]] { + //if (not std::isfinite(chisq)) [[unlikely]] { LogInfo << "Infinite chi2: " << chisq << std::endl << " bin " << bin_ << GET_VAR_NAME_VALUE(predVal) << std::endl @@ -212,7 +212,7 @@ namespace JointProbability{ << GET_VAR_NAME_VALUE(predMC->GetBinError(bin_)) << std::endl << GET_VAR_NAME_VALUE(predMC->GetBinContent(bin_)) << std::endl; // LogThrow("Bad chisq for bin"); - } + //} if(verboseLevel>=3){ LogTrace << "Bin #" << bin_ << ": chisq(" << chisq << ") / predVal(" << predVal << ") / dataVal(" << dataVal << ")" << std::endl; From 2bfefb3adea35b5ac9b57fb69c1c6e0aec6b66f3 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Wed, 7 Feb 2024 21:56:39 +0100 Subject: [PATCH 102/131] throw disable parameters as prior --- src/Applications/src/gundamMarginalise.cxx | 9 +++++++-- src/SamplesManager/src/JointProbability.cpp | 6 +++--- 2 files changed, 10 insertions(+), 5 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 60ec94712..7675daba5 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -37,7 +37,6 @@ double getParameterValueFromTextFile(std::string fileName, std::string parameter int main(int argc, char** argv){ - std::cout<<"Hello world 11.12.23"<Gaus(mu,sigma)); + } // LogInfo<<" "< "<GetBinError(bin_)) << std::endl << GET_VAR_NAME_VALUE(predMC->GetBinContent(bin_)) << std::endl; // LogThrow("Bad chisq for bin"); - //} + } if(verboseLevel>=3){ LogTrace << "Bin #" << bin_ << ": chisq(" << chisq << ") / predVal(" << predVal << ") / dataVal(" << dataVal << ")" << std::endl; From 63cec4b0ecdd675c9d76203e68c2422951113106 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Thu, 8 Feb 2024 15:59:38 +0100 Subject: [PATCH 103/131] included parameters thrown as prior (unconstrained, used in PTheta) --- src/Applications/src/gundamMarginalise.cxx | 45 +++++++++++++++++----- 1 file changed, 35 insertions(+), 10 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 7675daba5..87751b140 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -383,7 +383,7 @@ int main(int argc, char** argv){ LogInfo << " will not be marginalized out. \n"; } } - // Do the same for single parameters + // Do the same for single parameters for (auto &par: parSet.getParameterList()) { if (not par.isEnabled()) { LogInfo << "Parameter " << par.getName() << " is disabled" << std::endl; @@ -418,6 +418,15 @@ int main(int argc, char** argv){ } } } + // now deal with xsec disabled parameters (unconstrained by the ND, thrown according to prior) + for( auto& parSet : propagator.getParametersManager().getParameterSetsList() ) { + if (parSet.getName().find("Cross-section") != std::string::npos) { continue; } + // I want to throw single parameters from the xsec systematics, I don't care about the rest + for (auto &par: parSet.getParameterList()) { + if (par.isEnabled()) continue; + marg_param_list->Add(new TObjString(par.getFullTitle().c_str())); + } + } // write the list of parameters that will not be marginalised to the output file app.getOutfilePtr()->WriteObject(marg_param_list, "marg_param_list"); @@ -533,20 +542,14 @@ int main(int argc, char** argv){ if (not parSet.isEnabled()) { continue; } // LogInfo<< parSet.getName()<Gaus(mu,sigma)); - } + if (not par.isEnabled()) continue; // LogInfo<<" "< "<Gaus(mu, sigma); + } else { + LogInfo << " -> I don't know this prior. Prior type: " << par.getPriorType() << std::endl; + } + margThis.push_back(par.isMarginalised()); + parameters.push_back(value); + if(not par.isMarginalised()) + survivingParameterValues.push_back(value); + } + } // Write the ttrees margThrowTree->Fill(); ThrowsPThetaFormat->Fill(); From 4070d0f1948730e460b1fd43e049cafa9e61d470 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Wed, 14 Feb 2024 10:18:10 +0100 Subject: [PATCH 104/131] do not add to the margThrowTree the prior parameters --- src/Applications/src/gundamMarginalise.cxx | 1 - 1 file changed, 1 deletion(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 87751b140..e0e47ecbb 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -582,7 +582,6 @@ int main(int argc, char** argv){ LogInfo << " -> I don't know this prior. Prior type: " << par.getPriorType() << std::endl; } margThis.push_back(par.isMarginalised()); - parameters.push_back(value); if(not par.isMarginalised()) survivingParameterValues.push_back(value); } From 1124d72cf084cce0f888441a78211e1103227a35 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Wed, 14 Feb 2024 11:25:09 +0100 Subject: [PATCH 105/131] improved handling of unconstrained parameters to be thrown as priors --- src/Applications/src/gundamMarginalise.cxx | 51 ++++++++++++++-------- 1 file changed, 32 insertions(+), 19 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index e0e47ecbb..0e45a0d93 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -33,7 +33,7 @@ LoggerInit([]{ }); double getParameterValueFromTextFile(std::string fileName, std::string parameterName); - +bool isParUnconstrained(std::string parName, std::vector unconstrainedParNames); // TODO: implement this as parameterManager method int main(int argc, char** argv){ @@ -341,6 +341,9 @@ int main(int argc, char** argv){ int nToys{ clParser.getOptionVal("nToys") }; + // Get unconstrained parameters (thrown as prior in the PThetaThrowsTree) + std::vector unconstrainedParameters; + unconstrainedParameters = GenericToolbox::Json::fetchValue>(margConfig, "unconstrainedParameterList"); // Get parameters to be marginalised std::vector marginalisedParameters; std::vector marginalisedParameterSets; @@ -418,13 +421,13 @@ int main(int argc, char** argv){ } } } - // now deal with xsec disabled parameters (unconstrained by the ND, thrown according to prior) + // now deal with unconstrained parameters (unconstrained by the ND, thrown according to prior) for( auto& parSet : propagator.getParametersManager().getParameterSetsList() ) { - if (parSet.getName().find("Cross-section") != std::string::npos) { continue; } - // I want to throw single parameters from the xsec systematics, I don't care about the rest for (auto &par: parSet.getParameterList()) { if (par.isEnabled()) continue; - marg_param_list->Add(new TObjString(par.getFullTitle().c_str())); + if (isParUnconstrained(par.getFullTitle(), unconstrainedParameters)) { + marg_param_list->Add(new TObjString(par.getFullTitle().c_str())); + } } } @@ -564,26 +567,26 @@ int main(int argc, char** argv){ //debug LogInfo<<"LLH: "<Gaus(mu, sigma); - } else { - LogInfo << " -> I don't know this prior. Prior type: " << par.getPriorType() << std::endl; + double mu = par.getPriorValue(); + double sigma = par.getStdDevValue(); + double value = 0; + if (par.getPriorType() == 0) { + value = gRandom->Gaus(mu, sigma); + } else { + LogInfo << " -> I don't know this prior. Prior type: " << par.getPriorType() << std::endl; + } + margThis.push_back(par.isMarginalised()); + if (not par.isMarginalised()) + survivingParameterValues.push_back(value); } - margThis.push_back(par.isMarginalised()); - if(not par.isMarginalised()) - survivingParameterValues.push_back(value); } } // Write the ttrees @@ -652,4 +655,14 @@ double getParameterValueFromTextFile(std::string fileName="LargeWeight_parVector } LogInfo << "Parameter \"" << parameterName << "\" not found in file " << fileName << std::endl; return -999; +} + + +bool isParUnconstrained(std::string parName, std::vector unconstrainedParNames){ + for (auto &par: unconstrainedParNames){ + if (parName == par){ + return true; + } + } + return false; } \ No newline at end of file From e4684cf517a26f6af232ab6b632db2cf417e99c3 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Thu, 15 Feb 2024 13:49:25 +0100 Subject: [PATCH 106/131] Removed unconstrained parameters. Throws (I will do it in a separate macro) --- src/Applications/src/gundamMarginalise.cxx | 53 +--------------------- 1 file changed, 1 insertion(+), 52 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 0e45a0d93..4377018af 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -33,7 +33,6 @@ LoggerInit([]{ }); double getParameterValueFromTextFile(std::string fileName, std::string parameterName); -bool isParUnconstrained(std::string parName, std::vector unconstrainedParNames); // TODO: implement this as parameterManager method int main(int argc, char** argv){ @@ -147,12 +146,6 @@ int main(int argc, char** argv){ cHandler.override( margConfig ); LogInfo << "Override done." << std::endl; - - // read the parameters to include in the TTree - - // read the parameters to marginalise over - - if( clParser.isOptionTriggered("dryRun") ){ std::cout << cHandler.toString() << std::endl; @@ -341,9 +334,6 @@ int main(int argc, char** argv){ int nToys{ clParser.getOptionVal("nToys") }; - // Get unconstrained parameters (thrown as prior in the PThetaThrowsTree) - std::vector unconstrainedParameters; - unconstrainedParameters = GenericToolbox::Json::fetchValue>(margConfig, "unconstrainedParameterList"); // Get parameters to be marginalised std::vector marginalisedParameters; std::vector marginalisedParameterSets; @@ -367,13 +357,11 @@ int main(int argc, char** argv){ } else { LogInfo << "Set: " << parSet.getName().c_str(); for (int i = 0; i < marginalisedParameterSets.size(); i++) { - if (0 == std::strcmp(parSet.getName().c_str(), marginalisedParameterSets[i].c_str())) { setMatches = (true); break; } else { setMatches = (false); - } } if (setMatches) { @@ -421,15 +409,6 @@ int main(int argc, char** argv){ } } } - // now deal with unconstrained parameters (unconstrained by the ND, thrown according to prior) - for( auto& parSet : propagator.getParametersManager().getParameterSetsList() ) { - for (auto &par: parSet.getParameterList()) { - if (par.isEnabled()) continue; - if (isParUnconstrained(par.getFullTitle(), unconstrainedParameters)) { - marg_param_list->Add(new TObjString(par.getFullTitle().c_str())); - } - } - } // write the list of parameters that will not be marginalised to the output file app.getOutfilePtr()->WriteObject(marg_param_list, "marg_param_list"); @@ -567,28 +546,7 @@ int main(int argc, char** argv){ //debug LogInfo<<"LLH: "<Gaus(mu, sigma); - } else { - LogInfo << " -> I don't know this prior. Prior type: " << par.getPriorType() << std::endl; - } - margThis.push_back(par.isMarginalised()); - if (not par.isMarginalised()) - survivingParameterValues.push_back(value); - } - } - } + // Write the ttrees margThrowTree->Fill(); ThrowsPThetaFormat->Fill(); @@ -657,12 +615,3 @@ double getParameterValueFromTextFile(std::string fileName="LargeWeight_parVector return -999; } - -bool isParUnconstrained(std::string parName, std::vector unconstrainedParNames){ - for (auto &par: unconstrainedParNames){ - if (parName == par){ - return true; - } - } - return false; -} \ No newline at end of file From 0b0fa22ba444d2f356eeb9212c2c4a0ef8027eb0 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Tue, 5 Mar 2024 15:38:05 +0100 Subject: [PATCH 107/131] debug message --- src/Applications/src/gundamMarginalise.cxx | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 4377018af..9b795fcf7 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -379,7 +379,6 @@ int main(int argc, char** argv){ if (not par.isEnabled()) { LogInfo << "Parameter " << par.getName() << " is disabled" << std::endl; continue; - } bool matches = false; for (int i = 0; i < marginalisedParameters.size(); i++) { @@ -418,6 +417,10 @@ int main(int argc, char** argv){ // << par.getMaxPhysical() << " limits (mirr): " << par.getMinMirror() << " - " // << par.getMaxMirror() <<" --- marg? "<Print(); + LogInfo<<"\n\n"; // initializing variables for "parametersInject" mode double LLH_sum{0};// needed when injecting parameters manually From ab4ea78ebe08e173d4b12ea776743cc4086a725f Mon Sep 17 00:00:00 2001 From: lgiannes Date: Tue, 5 Mar 2024 15:50:09 +0100 Subject: [PATCH 108/131] adjusted, info on number of parameters --- src/Applications/src/gundamMarginalise.cxx | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 9b795fcf7..a4a8820b8 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -418,9 +418,7 @@ int main(int argc, char** argv){ // << par.getMaxMirror() <<" --- marg? "<Print(); - LogInfo<<"\n\n"; + LogInfo<<"Number of parameter in marg_param_list: "<GetEntries()< Date: Sun, 24 Mar 2024 23:41:23 +0900 Subject: [PATCH 109/131] implementation of sampling from multivariate t-student distribution (first) --- src/Applications/src/gundamMarginalise.cxx | 49 ++++++-- .../include/ParametersManager.h | 1 + .../src/ParametersManager.cpp | 105 ++++++++++++++++++ 3 files changed, 148 insertions(+), 7 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index a4a8820b8..160cd6c13 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -60,6 +60,7 @@ int main(int argc, char** argv){ clParser.addOption("parInject", {"--parameters-inject"}, "Input txt file for injecting params"); clParser.addOption("pedestal", {"--pedestal"}, "Add pedestal to the sampling distribution (percents)"); clParser.addOption("weightCap", {"--weight-cap"}, "Cap the weight of the throws (default: 1.e8)", 1); + clParser.addOption("tStudent", {"--t-student"}, "ndf for tStudent throw (default: gaussian throws)", 1); clParser.addDummyOption("Trigger options:"); clParser.addTriggerOption("dryRun", {"-d", "--dry-run"}, "Only overrides fitter config and print it."); @@ -179,6 +180,20 @@ int main(int argc, char** argv){ }else{ LogInfo << "Using default weight cap: "<< weightCap << std::endl; } + bool TStudent = false; + double tStudentNu = -1; + if(clParser.isOptionTriggered("tStudent")){ + tStudentNu = clParser.getOptionVal("tStudent"); + LogInfo << "Throwing according to a multivariate t-student with ndf: "<< tStudentNu << std::endl; + } + if(tStudentNu==-1){ + TStudent = false; + }else if(tStudentNu<=2){ + LogError<< "Not possible to use ndf for t-student smaller than 2!"<( cHandler.getConfig(), "fitterEngineConfig/propagatorConfig" ); @@ -240,7 +255,8 @@ int main(int argc, char** argv){ {"nToys", "nToys_%s"}, {"randomSeed", "Seed_%s"}, {"parInject", "parInject_%s"}, - {"pedestal", "ped_%s_percent"} + {"pedestal", "ped_%s_percent"}, + {"tStudent", "tStudentNu_%s"} }; outFilePath = GundamUtils::generateFileName(clParser, appendixDict) + ".root"; @@ -428,6 +444,7 @@ int main(int argc, char** argv){ LogInfo<< "Injecting parameters from file: " << parInjectFile << std::endl; } + double weightSum = 0, weightSquareSum = 0, ESS = 0; std::stringstream ss; ss << LogWarning.getPrefixString() << "Generating " << nToys << " toys..."; ////////////////////////////////////// // THROWS LOOP @@ -443,11 +460,16 @@ int main(int argc, char** argv){ if (usePedestal){ // if the pedestal option is enabled, a uniform distribution is added to the gaussian sampling distribution propagator.getParametersManager().throwParametersFromGlobalCovariance(weightsChiSquare, pedestalEntity, pedestalLeftEdge, pedestalRightEdge); - }else { - // standard case: throw according to the covariance matrix - propagator.getParametersManager().throwParametersFromGlobalCovariance(weightsChiSquare); + }else{ + if(!TStudent) { + // standard case: throw according to the covariance matrix + propagator.getParametersManager().throwParametersFromGlobalCovariance(weightsChiSquare); + }else{ + propagator.getParametersManager().throwParametersFromTStudent(weightsChiSquare,tStudentNu); + } } + if(injectParamsManually) { // count the number of parameters to be injected int nStripped{0}; @@ -541,13 +563,14 @@ int main(int argc, char** argv){ } LhOverGauss = exp(LLH-gLLH); if ( LLH-gLLH > log(weightCap)) { - LogInfo << "Throw " << iToy << " rejected: LLH-gLLH = " << LLH - gLLH << std::endl; + LogInfo << "Throw " << iToy << " over weight cap: LLH-gLLH = " << LLH - gLLH << std::endl; countBigThrows++; } //debug LogInfo<<"LLH: "<Fill(); ThrowsPThetaFormat->Fill(); @@ -562,7 +585,19 @@ int main(int argc, char** argv){ }// end of main throws loop LogInfo<<"weight cap: "< log(weightCap): "< log(weightCap): "<WriteObject(&ESS_writable,"ESS"); + app.getOutfilePtr()->WriteObject(&weightSum_writable,"weight_sum"); + app.getOutfilePtr()->WriteObject(&weightSquareSum_writable,"weight2_sum"); + double averageLLH = LLH_sum/nToys; epsilonNormAverage /= nToys; diff --git a/src/ParametersManager/include/ParametersManager.h b/src/ParametersManager/include/ParametersManager.h index 2cc4cbe1f..935666f59 100644 --- a/src/ParametersManager/include/ParametersManager.h +++ b/src/ParametersManager/include/ParametersManager.h @@ -50,6 +50,7 @@ class ParametersManager : public JsonBaseClass { void throwParametersFromGlobalCovariance(bool quietVerbose_ = true); void throwParametersFromGlobalCovariance(std::vector &weightsChiSquare); void throwParametersFromGlobalCovariance(std::vector &weightsChiSquare, double pedestalEntity, double pedestalLeftEdge, double pedestalRightEdge); + void throwParametersFromTStudent(std::vector &weightsChiSquare,double nu_); ParameterSet* getFitParameterSetPtr(const std::string& name_); // Logger related diff --git a/src/ParametersManager/src/ParametersManager.cpp b/src/ParametersManager/src/ParametersManager.cpp index 576ee2f79..649610711 100644 --- a/src/ParametersManager/src/ParametersManager.cpp +++ b/src/ParametersManager/src/ParametersManager.cpp @@ -394,6 +394,111 @@ void ParametersManager::throwParametersFromGlobalCovariance(std::vector } }// end of function +void ParametersManager::throwParametersFromTStudent(std::vector &weightsChiSquare,double nu_){ + // check that weightsChiSquare is an empty vector + LogThrowIf( not weightsChiSquare.empty(), "ERROR: argument weightsChiSquare is not empty" ); + + if( _strippedCovarianceMatrix_ == nullptr ){ + LogInfo << "Creating stripped global covariance matrix..." << std::endl; + LogThrowIf( _globalCovarianceMatrix_ == nullptr, "Global covariance matrix not set." ); + int nStripped{0}; + for( int iDiag = 0 ; iDiag < _globalCovarianceMatrix_->GetNrows() ; iDiag++ ){ + if( (*_globalCovarianceMatrix_)[iDiag][iDiag] != 0 ){ nStripped++; } + } + + LogInfo << "Stripped global covariance matrix is " << nStripped << "x" << nStripped << std::endl; + _strippedCovarianceMatrix_ = std::make_shared(nStripped, nStripped); + int iStrippedBin{-1}; + for( int iBin = 0 ; iBin < _globalCovarianceMatrix_->GetNrows() ; iBin++ ){ + if( (*_globalCovarianceMatrix_)[iBin][iBin] == 0 ){ continue; } + iStrippedBin++; + int jStrippedBin{-1}; + for( int jBin = 0 ; jBin < _globalCovarianceMatrix_->GetNrows() ; jBin++ ){ + if( (*_globalCovarianceMatrix_)[jBin][jBin] == 0 ){ continue; } + jStrippedBin++; + (*_strippedCovarianceMatrix_)[iStrippedBin][jStrippedBin] = (*_globalCovarianceMatrix_)[iBin][jBin]; + } + } + + _strippedParameterList_.reserve( nStripped ); + for( auto& parSet : _parameterSetList_ ){ + if( not parSet.isEnabled() ) continue; + for( auto& par : parSet.getParameterList() ){ + if( not par.isEnabled() ) continue; + _strippedParameterList_.emplace_back(&par); + } + } + LogThrowIf( _strippedParameterList_.size() != nStripped, "Enabled parameters list don't correspond to the matrix" ); + } + + if( _choleskyMatrix_ == nullptr ){ + LogInfo << "Generating global cholesky matrix" << std::endl; + _choleskyMatrix_ = std::shared_ptr( + GenericToolbox::getCholeskyMatrix(_strippedCovarianceMatrix_.get()) + ); + } + + bool keepThrowing{true}; +// int throwNb{0}; + + while( keepThrowing ){ +// throwNb++; + bool rethrow{false}; + std::vector throws,weights; + // calling Toolbox function to throw random parameters + GenericToolbox::throwTStudentParameters(_choleskyMatrix_.get(),nu_,throws, weights); + /////// + if(throws.size() != weights.size()){ + LogInfo<<"WARNING: throws.size() != weights.size() "<< throws.size()<GetNrows() ; iPar++ ){ + _strippedParameterList_[iPar]->setParameterValue( + _strippedParameterList_[iPar]->getPriorValue() + + throws[iPar] + ); + weightsChiSquare.push_back(weights[iPar]); + + if( _reThrowParSetIfOutOfBounds_ ){ + if( not _strippedParameterList_[iPar]->isValueWithinBounds() ){ + // re-do the throwing +// LogDebug << "Not within bounds: " << _strippedParameterList_[iPar]->getSummary() << std::endl; + rethrow = true; + } + } + } + + // Making sure eigen decomposed parameters get the conversion done + for( auto& parSet : _parameterSetList_ ){ + if( not parSet.isEnabled() ) continue; + if( parSet.isUseEigenDecompInFit() ){ + parSet.propagateOriginalToEigen(); + + // also check the bounds of real parameter space + if( _reThrowParSetIfOutOfBounds_ ){ + for( auto& par : parSet.getEigenParameterList() ){ + if( not par.isEnabled() ) continue; + if( not par.isValueWithinBounds() ){ + // re-do the throwing + rethrow = true; + break; + } + } + } + } + } + + + if( rethrow ){ + // wrap back to the while loop +// LogDebug << "RE-THROW #" << throwNb << std::endl; + continue; + } + + // reached this point: all parameters are within bounds + keepThrowing = false; + } +} + void ParametersManager::injectParameterValues(const nlohmann::json &config_) { LogWarning << "Injecting parameters..." << std::endl; From 6ccb15035bcb2cf8cf8f203f1a47dd0c582de319 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Mon, 25 Mar 2024 15:48:30 +0900 Subject: [PATCH 110/131] counterattack to double limit hit (first attempt) --- src/Applications/src/gundamMarginalise.cxx | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 160cd6c13..3be98a4a6 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -445,6 +445,7 @@ int main(int argc, char** argv){ } double weightSum = 0, weightSquareSum = 0, ESS = 0; + int weightSumE50 = 0, weightSquareSumE50 = 0; std::stringstream ss; ss << LogWarning.getPrefixString() << "Generating " << nToys << " toys..."; ////////////////////////////////////// // THROWS LOOP @@ -571,6 +572,14 @@ int main(int argc, char** argv){ weightSum += LhOverGauss; weightSquareSum += LhOverGauss*LhOverGauss; + while(weightSum>1.e50){ + weightSum /= 1.e50; + weightSumE50++; + } + while(weightSquareSum>1.e50){ + weightSquareSum /= 1.e50; + weightSquareSumE50++; + } // Write the ttrees margThrowTree->Fill(); ThrowsPThetaFormat->Fill(); @@ -584,13 +593,15 @@ int main(int argc, char** argv){ // } }// end of main throws loop - LogInfo<<"weight cap: "< log(weightCap): "< Date: Sun, 31 Mar 2024 16:54:26 +0200 Subject: [PATCH 111/131] SampleSet.cpp: reduce log --- src/SamplesManager/src/SampleSet.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/SamplesManager/src/SampleSet.cpp b/src/SamplesManager/src/SampleSet.cpp index 3a3d5c0f7..8c03ce9e3 100644 --- a/src/SamplesManager/src/SampleSet.cpp +++ b/src/SamplesManager/src/SampleSet.cpp @@ -104,7 +104,7 @@ double SampleSet::evalLikelihood(){ double llh = 0.; for( auto& sample : _fitSampleList_ ){ llh += this->evalLikelihood(sample); - LogErrorIf(std::isnan(llh) or std::isinf(llh)) << sample.getName() << ": reported likelihood is invalid:" << llh; +// LogErrorIf(std::isnan(llh) or std::isinf(llh)) << sample.getName() << ": reported likelihood is invalid:" << llh; // LogThrowIf(std::isnan(llh) or std::isinf(llh), sample.getName() << ": reported likelihood is invalid:" << llh); } return llh; From f0a1bc297e0047ced3da74eac81ef8981efe1536 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Tue, 2 Apr 2024 12:23:53 +0200 Subject: [PATCH 112/131] SampleSet.cpp: LogThrowIf instead of LogErrorIf --- src/SamplesManager/src/JointProbability.cpp | 2 +- src/SamplesManager/src/SampleSet.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/SamplesManager/src/JointProbability.cpp b/src/SamplesManager/src/JointProbability.cpp index 17919238c..9ae731beb 100644 --- a/src/SamplesManager/src/JointProbability.cpp +++ b/src/SamplesManager/src/JointProbability.cpp @@ -211,7 +211,7 @@ namespace JointProbability{ << GET_VAR_NAME_VALUE(nomMC->GetBinError(bin_)) << std::endl << GET_VAR_NAME_VALUE(predMC->GetBinError(bin_)) << std::endl << GET_VAR_NAME_VALUE(predMC->GetBinContent(bin_)) << std::endl; - // LogThrow("Bad chisq for bin"); + LogThrow("Bad chisq for bin"); } if(verboseLevel>=3){ diff --git a/src/SamplesManager/src/SampleSet.cpp b/src/SamplesManager/src/SampleSet.cpp index 8c03ce9e3..221fee2dc 100644 --- a/src/SamplesManager/src/SampleSet.cpp +++ b/src/SamplesManager/src/SampleSet.cpp @@ -105,7 +105,7 @@ double SampleSet::evalLikelihood(){ for( auto& sample : _fitSampleList_ ){ llh += this->evalLikelihood(sample); // LogErrorIf(std::isnan(llh) or std::isinf(llh)) << sample.getName() << ": reported likelihood is invalid:" << llh; -// LogThrowIf(std::isnan(llh) or std::isinf(llh), sample.getName() << ": reported likelihood is invalid:" << llh); + LogThrowIf(std::isnan(llh) or std::isinf(llh), sample.getName() << ": reported likelihood is invalid:" << llh); } return llh; } From aad3f4ada53790aa1420c464f173dfee9f19d81a Mon Sep 17 00:00:00 2001 From: lgiannes Date: Tue, 2 Apr 2024 12:48:36 +0200 Subject: [PATCH 113/131] SampleSet.cpp: LogThrowIf instead of LogErrorIf --- src/SamplesManager/src/JointProbability.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/SamplesManager/src/JointProbability.cpp b/src/SamplesManager/src/JointProbability.cpp index 9ae731beb..f8c83d753 100644 --- a/src/SamplesManager/src/JointProbability.cpp +++ b/src/SamplesManager/src/JointProbability.cpp @@ -211,7 +211,7 @@ namespace JointProbability{ << GET_VAR_NAME_VALUE(nomMC->GetBinError(bin_)) << std::endl << GET_VAR_NAME_VALUE(predMC->GetBinError(bin_)) << std::endl << GET_VAR_NAME_VALUE(predMC->GetBinContent(bin_)) << std::endl; - LogThrow("Bad chisq for bin"); + //LogThrow("Bad chisq for bin"); } if(verboseLevel>=3){ From 14d717fcee996af3dbde2a9ad248b2f462997094 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Tue, 2 Apr 2024 14:52:18 +0200 Subject: [PATCH 114/131] JointProbability.cpp: remove error if chisq=infinity --- src/SamplesManager/src/JointProbability.cpp | 24 ++++++++++----------- 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/src/SamplesManager/src/JointProbability.cpp b/src/SamplesManager/src/JointProbability.cpp index f8c83d753..2ee882fcb 100644 --- a/src/SamplesManager/src/JointProbability.cpp +++ b/src/SamplesManager/src/JointProbability.cpp @@ -199,18 +199,18 @@ namespace JointProbability{ } if (not std::isfinite(chisq)) [[unlikely]] { - LogError << "Infinite chi2: " << chisq << std::endl - << " bin " << bin_ - << GET_VAR_NAME_VALUE(predVal) << std::endl - << GET_VAR_NAME_VALUE(dataVal) << std::endl - << GET_VAR_NAME_VALUE(newmc) << std::endl - << GET_VAR_NAME_VALUE(stat) << std::endl - << GET_VAR_NAME_VALUE(penalty) << std::endl - << GET_VAR_NAME_VALUE(mcuncert) << std::endl - << GET_VAR_NAME_VALUE(nomMC->GetBinContent(bin_)) << std::endl - << GET_VAR_NAME_VALUE(nomMC->GetBinError(bin_)) << std::endl - << GET_VAR_NAME_VALUE(predMC->GetBinError(bin_)) << std::endl - << GET_VAR_NAME_VALUE(predMC->GetBinContent(bin_)) << std::endl; +// LogError << "Infinite chi2: " << chisq << std::endl +// << " bin " << bin_ +// << GET_VAR_NAME_VALUE(predVal) << std::endl +// << GET_VAR_NAME_VALUE(dataVal) << std::endl +// << GET_VAR_NAME_VALUE(newmc) << std::endl +// << GET_VAR_NAME_VALUE(stat) << std::endl +// << GET_VAR_NAME_VALUE(penalty) << std::endl +// << GET_VAR_NAME_VALUE(mcuncert) << std::endl +// << GET_VAR_NAME_VALUE(nomMC->GetBinContent(bin_)) << std::endl +// << GET_VAR_NAME_VALUE(nomMC->GetBinError(bin_)) << std::endl +// << GET_VAR_NAME_VALUE(predMC->GetBinError(bin_)) << std::endl +// << GET_VAR_NAME_VALUE(predMC->GetBinContent(bin_)) << std::endl; //LogThrow("Bad chisq for bin"); } From 9255435b49eaff041be34f156744aac3acae83f0 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Tue, 2 Apr 2024 15:04:34 +0200 Subject: [PATCH 115/131] JointProbability.cpp --- src/SamplesManager/src/JointProbability.cpp | 26 ++++++++++----------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/src/SamplesManager/src/JointProbability.cpp b/src/SamplesManager/src/JointProbability.cpp index 2ee882fcb..9ae731beb 100644 --- a/src/SamplesManager/src/JointProbability.cpp +++ b/src/SamplesManager/src/JointProbability.cpp @@ -199,19 +199,19 @@ namespace JointProbability{ } if (not std::isfinite(chisq)) [[unlikely]] { -// LogError << "Infinite chi2: " << chisq << std::endl -// << " bin " << bin_ -// << GET_VAR_NAME_VALUE(predVal) << std::endl -// << GET_VAR_NAME_VALUE(dataVal) << std::endl -// << GET_VAR_NAME_VALUE(newmc) << std::endl -// << GET_VAR_NAME_VALUE(stat) << std::endl -// << GET_VAR_NAME_VALUE(penalty) << std::endl -// << GET_VAR_NAME_VALUE(mcuncert) << std::endl -// << GET_VAR_NAME_VALUE(nomMC->GetBinContent(bin_)) << std::endl -// << GET_VAR_NAME_VALUE(nomMC->GetBinError(bin_)) << std::endl -// << GET_VAR_NAME_VALUE(predMC->GetBinError(bin_)) << std::endl -// << GET_VAR_NAME_VALUE(predMC->GetBinContent(bin_)) << std::endl; - //LogThrow("Bad chisq for bin"); + LogError << "Infinite chi2: " << chisq << std::endl + << " bin " << bin_ + << GET_VAR_NAME_VALUE(predVal) << std::endl + << GET_VAR_NAME_VALUE(dataVal) << std::endl + << GET_VAR_NAME_VALUE(newmc) << std::endl + << GET_VAR_NAME_VALUE(stat) << std::endl + << GET_VAR_NAME_VALUE(penalty) << std::endl + << GET_VAR_NAME_VALUE(mcuncert) << std::endl + << GET_VAR_NAME_VALUE(nomMC->GetBinContent(bin_)) << std::endl + << GET_VAR_NAME_VALUE(nomMC->GetBinError(bin_)) << std::endl + << GET_VAR_NAME_VALUE(predMC->GetBinError(bin_)) << std::endl + << GET_VAR_NAME_VALUE(predMC->GetBinContent(bin_)) << std::endl; + LogThrow("Bad chisq for bin"); } if(verboseLevel>=3){ From aab661468c1886c504c1f9d89e9d6fe686ad58b3 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Wed, 3 Apr 2024 18:49:12 +0200 Subject: [PATCH 116/131] remove all error/throws in JointProbability.cpp and SampleSet.cpp --- src/SamplesManager/src/JointProbability.cpp | 26 ++++++++++----------- src/SamplesManager/src/SampleSet.cpp | 2 +- 2 files changed, 14 insertions(+), 14 deletions(-) diff --git a/src/SamplesManager/src/JointProbability.cpp b/src/SamplesManager/src/JointProbability.cpp index 9ae731beb..11f81e98b 100644 --- a/src/SamplesManager/src/JointProbability.cpp +++ b/src/SamplesManager/src/JointProbability.cpp @@ -199,19 +199,19 @@ namespace JointProbability{ } if (not std::isfinite(chisq)) [[unlikely]] { - LogError << "Infinite chi2: " << chisq << std::endl - << " bin " << bin_ - << GET_VAR_NAME_VALUE(predVal) << std::endl - << GET_VAR_NAME_VALUE(dataVal) << std::endl - << GET_VAR_NAME_VALUE(newmc) << std::endl - << GET_VAR_NAME_VALUE(stat) << std::endl - << GET_VAR_NAME_VALUE(penalty) << std::endl - << GET_VAR_NAME_VALUE(mcuncert) << std::endl - << GET_VAR_NAME_VALUE(nomMC->GetBinContent(bin_)) << std::endl - << GET_VAR_NAME_VALUE(nomMC->GetBinError(bin_)) << std::endl - << GET_VAR_NAME_VALUE(predMC->GetBinError(bin_)) << std::endl - << GET_VAR_NAME_VALUE(predMC->GetBinContent(bin_)) << std::endl; - LogThrow("Bad chisq for bin"); +// LogError << "Infinite chi2: " << chisq << std::endl +// << " bin " << bin_ +// << GET_VAR_NAME_VALUE(predVal) << std::endl +// << GET_VAR_NAME_VALUE(dataVal) << std::endl +// << GET_VAR_NAME_VALUE(newmc) << std::endl +// << GET_VAR_NAME_VALUE(stat) << std::endl +// << GET_VAR_NAME_VALUE(penalty) << std::endl +// << GET_VAR_NAME_VALUE(mcuncert) << std::endl +// << GET_VAR_NAME_VALUE(nomMC->GetBinContent(bin_)) << std::endl +// << GET_VAR_NAME_VALUE(nomMC->GetBinError(bin_)) << std::endl +// << GET_VAR_NAME_VALUE(predMC->GetBinError(bin_)) << std::endl +// << GET_VAR_NAME_VALUE(predMC->GetBinContent(bin_)) << std::endl; +// LogThrow("Bad chisq for bin"); } if(verboseLevel>=3){ diff --git a/src/SamplesManager/src/SampleSet.cpp b/src/SamplesManager/src/SampleSet.cpp index 221fee2dc..d3e0776f9 100644 --- a/src/SamplesManager/src/SampleSet.cpp +++ b/src/SamplesManager/src/SampleSet.cpp @@ -105,7 +105,7 @@ double SampleSet::evalLikelihood(){ for( auto& sample : _fitSampleList_ ){ llh += this->evalLikelihood(sample); // LogErrorIf(std::isnan(llh) or std::isinf(llh)) << sample.getName() << ": reported likelihood is invalid:" << llh; - LogThrowIf(std::isnan(llh) or std::isinf(llh), sample.getName() << ": reported likelihood is invalid:" << llh); +// LogThrowIf(std::isnan(llh) or std::isinf(llh), sample.getName() << ": reported likelihood is invalid:" << llh); } return llh; } From 330786e386fc91cffa34833372504038b1e463dd Mon Sep 17 00:00:00 2001 From: lgiannes Date: Thu, 4 Apr 2024 10:58:35 +0200 Subject: [PATCH 117/131] gitmodules --- .gitmodules | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/.gitmodules b/.gitmodules index 430ed230c..7542918c0 100644 --- a/.gitmodules +++ b/.gitmodules @@ -2,6 +2,10 @@ path = submodules/simple-cpp-logger url = https://github.com/gundam-organization/simple-cpp-logger.git branch = main +[submodule "submodules/cpp-generic-toolbox"] + path = submodules/cpp-generic-toolbox + url = https://github.com/gundam-organization/cpp-generic-toolbox.git + branch = main [submodule "submodules/simple-cpp-cmd-line-parser"] path = submodules/simple-cpp-cmd-line-parser url = https://github.com/gundam-organization/simple-cpp-cmd-line-parser.git From 24564b94858971f9b2f125dd847668f73e0d5317 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Thu, 4 Apr 2024 12:18:10 +0200 Subject: [PATCH 118/131] added sampling function to GundamUtils.h --- src/Utils/include/GundamUtils.h | 102 ++++++++++++++++++++++++++++++++ 1 file changed, 102 insertions(+) diff --git a/src/Utils/include/GundamUtils.h b/src/Utils/include/GundamUtils.h index 2e3fce89f..5a123a5d1 100644 --- a/src/Utils/include/GundamUtils.h +++ b/src/Utils/include/GundamUtils.h @@ -92,6 +92,108 @@ namespace GundamUtils { }; + inline void throwCorrelatedParameters(TMatrixD* choleskyCovMatrix_, std::vector& thrownParListOut_){ + if( choleskyCovMatrix_ == nullptr ) return; + if( thrownParListOut_.size() != choleskyCovMatrix_->GetNcols() ){ + thrownParListOut_.resize(choleskyCovMatrix_->GetNcols(), 0); + } + TVectorD thrownParVec(choleskyCovMatrix_->GetNcols()); + for( int iPar = 0 ; iPar < choleskyCovMatrix_->GetNcols() ; iPar++ ){ + thrownParVec[iPar] = gRandom->Gaus(); + } + thrownParVec *= (*choleskyCovMatrix_); + for( int iPar = 0 ; iPar < choleskyCovMatrix_->GetNcols() ; iPar++ ){ + thrownParListOut_.at(iPar) = thrownParVec[iPar]; + } + }// end of function throwCorrelatedParameters(TMatrixD* choleskyCovMatrix_, std::vector& thrownParListOut_) + inline std::vector throwCorrelatedParameters(TMatrixD* choleskyCovMatrix_){ + std::vector out; + throwCorrelatedParameters(choleskyCovMatrix_, out); + return out; + }// end of function throwCorrelatedParameters(TMatrixD* choleskyCovMatrix_) + inline void throwCorrelatedParameters( TMatrixD* choleskyCovMatrix_, std::vector& thrownParListOut_, std::vector& weights, + double pedestalEntity, double pedestalLeftEdge, double pedestalRightEdge + ){ + + double pi = TMath::Pi(); + double NormalizingFactor = 1.0 / (TMath::Sqrt(2.0 * pi)); + double pedestalRange = pedestalRightEdge - pedestalLeftEdge; + if( choleskyCovMatrix_ == nullptr ) return; + if( thrownParListOut_.size() != choleskyCovMatrix_->GetNcols() ){ + thrownParListOut_.resize(choleskyCovMatrix_->GetNcols(), 0); + } + if( weights.size() != choleskyCovMatrix_->GetNcols() ){ + weights.resize(choleskyCovMatrix_->GetNcols(), 0); + } + + TVectorD thrownParVec(choleskyCovMatrix_->GetNcols()); + double choice = gRandom->Uniform(0,1); + if (choice>pedestalEntity) { + for (int iPar = 0; iPar < choleskyCovMatrix_->GetNcols(); iPar++) { + thrownParVec[iPar] = gRandom->Gaus(); + if (thrownParVec[iPar]>pedestalLeftEdge and thrownParVec[iPar]GetNcols(); iPar++) { + thrownParVec[iPar] = gRandom->Uniform(pedestalLeftEdge, pedestalRightEdge); + weights.at(iPar) = -TMath::Log( + pedestalEntity*1.0/pedestalRange + (1.0-pedestalEntity) * NormalizingFactor * TMath::Exp(-0.500 * thrownParVec[iPar] * thrownParVec[iPar]) + ); + } + } + thrownParVec *= (*choleskyCovMatrix_); + for( int iPar = 0 ; iPar < choleskyCovMatrix_->GetNcols() ; iPar++ ){ + thrownParListOut_.at(iPar) = thrownParVec[iPar]; + } + }// end of function throwCorrelatedParameters(TMatrixD* choleskyCovMatrix_, std::vector& thrownParListOut_, std::vector& weights, double pedestalEntity, double pedestalLeftEdge, double pedestalRightEdge) + + inline void throwCorrelatedParameters(TMatrixD* choleskyCovMatrix_, std::vector& thrownParListOut_, std::vector& weights){ + throwCorrelatedParameters(choleskyCovMatrix_, thrownParListOut_, weights, 0, 0, 0); + }// end of function throwCorrelatedParameters(TMatrixD* choleskyCovMatrix_, std::vector& thrownParListOut_, std::vector& weights) + + inline void throwTStudentParameters(TMatrixD* choleskyCovMatrix_, double nu_, std::vector& thrownParListOut_, std::vector& weights){ + // Simple sanity check + if( choleskyCovMatrix_ == nullptr ) return; + if( thrownParListOut_.size() != choleskyCovMatrix_->GetNcols() ){ + thrownParListOut_.resize(choleskyCovMatrix_->GetNcols(), 0); + } + if( weights.size() != choleskyCovMatrix_->GetNcols() ){ + weights.resize(choleskyCovMatrix_->GetNcols(), 0); + } + // Throw N independent normal distributions + TVectorD thrownParVec(choleskyCovMatrix_->GetNcols()); + for( int iPar = 0 ; iPar < choleskyCovMatrix_->GetNcols() ; iPar++ ){ + thrownParVec[iPar] = gRandom->Gaus(); + } + // Multiply by cov matrix to obtain the gaussian part of the t-student (expanded as the covariance matrix says) + TVectorD thrownParVecExpanded = (*choleskyCovMatrix_)*thrownParVec; + std::vector chiSquareForStudentT(choleskyCovMatrix_->GetNcols()); + int p = 1; // because only THEN I sum over all dimensions. At this level it's all independent single-dim variables + for( int iPar = 0 ; iPar < choleskyCovMatrix_->GetNcols() ; iPar++ ){ + // Throw a chisquare with nu_ degrees of freedom + double chiSquareProb = gRandom->Uniform(0,1); + chiSquareForStudentT.at(iPar) = TMath::ChisquareQuantile(chiSquareProb, nu_); + // Build the t-student throw by multiplying the multivariate gaussian (with input cov. matrix) and the chisquare + thrownParListOut_.at(iPar) = thrownParVecExpanded[iPar] * sqrt(nu_/chiSquareForStudentT.at(iPar)); + // Fill the weights vector now (according to the pdf of the t-student multivariate distribution) + double logFactor1 = TMath::LnGamma(0.5*(nu_+p)) - TMath::LnGamma(0.5*nu_); + double logDenominator = + (0.5*p)*TMath::Log(nu_*TMath::Pi()); + double normalizedTStudentThrow = thrownParVec[iPar] * sqrt(nu_/chiSquareForStudentT.at(iPar)); + double logFactor2 = -0.5*(nu_+p)*TMath::Log( 1 + 1/nu_*normalizedTStudentThrow*normalizedTStudentThrow ); + weights.at(iPar) = -(logFactor1) +(logDenominator) - logFactor2; + //std::cout<& thrownParListOut_, std::vector& weights) + } #endif //GUNDAM_GUNDAMUTILS_H From 2cfa4b03923d4551b3a986e639d2375e1b109d2d Mon Sep 17 00:00:00 2001 From: lgiannes Date: Fri, 5 Apr 2024 11:47:23 +0200 Subject: [PATCH 119/131] gundamMarginalise.cxx compliant with gundam version 1.9.0 --- src/Applications/src/gundamMarginalise.cxx | 70 +++++++++++----------- 1 file changed, 35 insertions(+), 35 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 3be98a4a6..ddd8809ea 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -22,6 +22,8 @@ #include "TDirectory.h" #include "TH1D.h" #include "TH2D.h" +#include "../../Fitter/Engine/include/FitterEngine.h" +#include "../../StatisticalInference/Likelihood/include/LikelihoodInterface.h" #include #include @@ -113,7 +115,7 @@ int main(int argc, char** argv){ - nlohmann::json fitterConfig; + JsonType fitterConfig; ObjectReader::readObject(fitterFile.get(), {{"gundam/config_TNamed"}, {"gundamFitter/unfoldedConfig_TNamed"}}, [&](TNamed* config_){ fitterConfig = GenericToolbox::Json::readConfigJsonStr( config_->GetTitle() ); }); @@ -121,6 +123,7 @@ int main(int argc, char** argv){ ConfigUtils::ConfigHandler cHandler{ fitterConfig }; + // // Disabling defined samples: // LogInfo << "Removing defined samples..." << std::endl; // ConfigUtils::applyOverrides( @@ -144,7 +147,7 @@ int main(int argc, char** argv){ // normal behavior margConfig = margConfig[0]; } - cHandler.override( margConfig ); + cHandler.override( (JsonType)margConfig ); LogInfo << "Override done." << std::endl; if( clParser.isOptionTriggered("dryRun") ){ @@ -180,37 +183,39 @@ int main(int argc, char** argv){ }else{ LogInfo << "Using default weight cap: "<< weightCap << std::endl; } - bool TStudent = false; + bool tStudent = false; double tStudentNu = -1; if(clParser.isOptionTriggered("tStudent")){ tStudentNu = clParser.getOptionVal("tStudent"); LogInfo << "Throwing according to a multivariate t-student with ndf: "<< tStudentNu << std::endl; } if(tStudentNu==-1){ - TStudent = false; + tStudent = false; }else if(tStudentNu<=2){ LogError<< "Not possible to use ndf for t-student smaller than 2!"<( cHandler.getConfig(), "fitterEngineConfig/propagatorConfig" ); - // Create a propagator object - Propagator propagator; + // Initialize the fitterEngine + LogInfo << "FitterEngine setup..." << std::endl; + FitterEngine fitter{nullptr}; + fitter.readConfig(GenericToolbox::Json::fetchSubEntry((JsonType)cHandler.getConfig(), {"fitterEngineConfig"})); -// // Read the whole fitter config with the override parameters - propagator.readConfig( configPropagator ); + DataSetManager& dataSetManager{fitter.getLikelihoodInterface().getDataSetManager()}; // We are only interested in our MC. Data has already been used to get the post-fit error/values - propagator.setLoadAsimovData( true ); + dataSetManager.getPropagator().setLoadAsimovData( true ); // Disabling eigen decomposed parameters - propagator.setEnableEigenToOrigInPropagate( false ); + dataSetManager.getPropagator().setEnableEigenToOrigInPropagate( false ); // Load everything - propagator.initialize(); + dataSetManager.initialize(); + Propagator& propagator{dataSetManager.getPropagator()}; propagator.getParametersManager().getParameterSetsList(); @@ -237,11 +242,6 @@ int main(int argc, char** argv){ } } }); - // also save the value of the LLH at the best fit point: - propagator.propagateParametersOnSamples(); - propagator.updateLlhCache(); - double bestFitLLH = propagator.getLlhBuffer(); - LogInfo<<"LLH (computed on injected parameters): "<( fitterFile.get(), "FitterEngine/postFit/Hesse/hessian/postfitCovarianceOriginal_TH2D", @@ -279,10 +282,16 @@ int main(int argc, char** argv){ app.setCmdLinePtr( &clParser ); - app.setConfigString( ConfigUtils::ConfigHandler{margConfig}.toString() ); + app.setConfigString( ConfigUtils::ConfigHandler{(JsonType)margConfig}.toString() ); app.openOutputFile( outFilePath ); app.writeAppInfo(); + + // also save the value of the LLH at the best fit point: + propagator.propagateParameters(); + fitter.getLikelihoodInterface().propagateAndEvalLikelihood(); + double bestFitLLH = fitter.getLikelihoodInterface().getBuffer().totalLikelihood; + LogInfo<<"LLH (computed on injected parameters): "<mkdir("postFitInfo"); postFitInfo->cd(); @@ -462,7 +471,7 @@ int main(int argc, char** argv){ // if the pedestal option is enabled, a uniform distribution is added to the gaussian sampling distribution propagator.getParametersManager().throwParametersFromGlobalCovariance(weightsChiSquare, pedestalEntity, pedestalLeftEdge, pedestalRightEdge); }else{ - if(!TStudent) { + if(!tStudent) { // standard case: throw according to the covariance matrix propagator.getParametersManager().throwParametersFromGlobalCovariance(weightsChiSquare); }else{ @@ -511,7 +520,7 @@ int main(int argc, char** argv){ // If is in eigen space, propagateOriginalToEigen for (auto &parSet: propagator.getParametersManager().getParameterSetsList()) { if (not parSet.isEnabled()) { continue; } - if (parSet.isUseEigenDecompInFit()){ + if (parSet.isEnableEigenDecomp()){ parSet.propagateOriginalToEigen(); } } @@ -519,8 +528,8 @@ int main(int argc, char** argv){ }// end if(injectParamsManually) //LogInfo<<"Computing LH..."<Fill(); ThrowsPThetaFormat->Fill(); - //debug -// if(iToy==0 || iToy==1){ -// LogInfo<<"DEBUG: Toy "< log(weightCap): "< Date: Fri, 5 Apr 2024 11:48:11 +0200 Subject: [PATCH 120/131] sampler in GundamUtils (not anymore in generic-toolbox) --- src/ParametersManager/include/Parameter.h | 11 +---------- src/ParametersManager/src/ParametersManager.cpp | 10 +++++----- src/Propagator/src/Propagator.cpp | 2 +- src/Utils/include/GundamUtils.h | 1 + 4 files changed, 8 insertions(+), 16 deletions(-) diff --git a/src/ParametersManager/include/Parameter.h b/src/ParametersManager/include/Parameter.h index dc2880d2e..03c3a38a3 100644 --- a/src/ParametersManager/include/Parameter.h +++ b/src/ParametersManager/include/Parameter.h @@ -81,12 +81,8 @@ class Parameter : public JsonBaseClass { [[nodiscard]] const std::string &getName() const{ return _name_; } [[nodiscard]] const JsonType &getDialDefinitionsList() const{ return _dialDefinitionsList_; } [[nodiscard]] const ParameterSet *getOwner() const{ return _owner_; } -<<<<<<< HEAD - [[nodiscard]] PriorType::PriorType getPriorType() const{ return _priorType_; } - [[nodiscard]] bool isMarginalised() const{ return _isMarginalised_; } -======= [[nodiscard]] PriorType getPriorType() const{ return _priorType_; } ->>>>>>> origin/main + [[nodiscard]] bool isMarginalised() const{ return _isMarginalised_; } // Core void setValueAtPrior(){ _parameterValue_ = _priorValue_; } @@ -118,14 +114,9 @@ class Parameter : public JsonBaseClass { double _stepSize_{std::nan("unset")}; std::string _name_{}; std::string _dialsWorkingDirectory_{"."}; -<<<<<<< HEAD - nlohmann::json _parameterConfig_{}; - nlohmann::json _dialDefinitionsList_{}; bool _isMarginalised_{false}; -======= JsonType _parameterConfig_{}; JsonType _dialDefinitionsList_{}; ->>>>>>> origin/main // Internals const ParameterSet* _owner_{nullptr}; diff --git a/src/ParametersManager/src/ParametersManager.cpp b/src/ParametersManager/src/ParametersManager.cpp index e64a4345d..dfb2d1785 100644 --- a/src/ParametersManager/src/ParametersManager.cpp +++ b/src/ParametersManager/src/ParametersManager.cpp @@ -8,6 +8,7 @@ #include "GenericToolbox.Utils.h" #include "GenericToolbox.Json.h" #include "Logger.h" +#include "GundamUtils.h" #include @@ -217,7 +218,7 @@ void ParametersManager::throwParametersFromGlobalCovariance(bool quietVerbose_){ while( keepThrowing ){ throwNb++; bool rethrow{false}; - auto throws = GenericToolbox::throwCorrelatedParameters(_choleskyMatrix_.get()); + auto throws = GundamUtils::throwCorrelatedParameters(_choleskyMatrix_.get()); for( int iPar = 0 ; iPar < _choleskyMatrix_->GetNrows() ; iPar++ ){ auto* parPtr = _strippedParameterList_[iPar]; parPtr->setParameterValue( parPtr->getPriorValue() + throws[iPar] ); @@ -290,7 +291,6 @@ void ParametersManager::throwParametersFromGlobalCovariance(bool quietVerbose_){ } } -<<<<<<< HEAD void ParametersManager::throwParametersFromGlobalCovariance(std::vector &weightsChiSquare){ throwParametersFromGlobalCovariance(weightsChiSquare,0,0,0); }// end of function @@ -350,9 +350,9 @@ void ParametersManager::throwParametersFromGlobalCovariance(std::vector bool rethrow{false}; std::vector throws,weights; if(pedestalEntity==0){ - GenericToolbox::throwCorrelatedParameters(_choleskyMatrix_.get(),throws, weights); + GundamUtils::throwCorrelatedParameters(_choleskyMatrix_.get(),throws, weights); }else{ - GenericToolbox::throwCorrelatedParameters(_choleskyMatrix_.get(),throws, weights, + GundamUtils::throwCorrelatedParameters(_choleskyMatrix_.get(),throws, weights, pedestalEntity,pedestalLeftEdge,pedestalRightEdge); } if(throws.size() != weights.size()){ @@ -458,7 +458,7 @@ void ParametersManager::throwParametersFromTStudent(std::vector &weights bool rethrow{false}; std::vector throws,weights; // calling Toolbox function to throw random parameters - GenericToolbox::throwTStudentParameters(_choleskyMatrix_.get(),nu_,throws, weights); + GundamUtils::throwTStudentParameters(_choleskyMatrix_.get(),nu_,throws, weights); /////// if(throws.size() != weights.size()){ LogInfo<<"WARNING: throws.size() != weights.size() "<< throws.size()< Date: Fri, 5 Apr 2024 11:50:12 +0200 Subject: [PATCH 121/131] gundamMarginalise compliant with Gundam 1.9.0 --- src/Applications/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Applications/CMakeLists.txt b/src/Applications/CMakeLists.txt index 9e86d63d7..5adb32fcf 100644 --- a/src/Applications/CMakeLists.txt +++ b/src/Applications/CMakeLists.txt @@ -82,7 +82,7 @@ target_link_libraries( gundamConfigUnfolder GundamUtils ) target_link_libraries( gundamConfigCompare GundamUtils ) target_link_libraries( gundamPlotExtractor GundamUtils ) target_link_libraries( gundamCalcXsec GundamPropagator ) -target_link_libraries( gundamMarginalise GundamPropagator ) +target_link_libraries( gundamMarginalise GundamPropagator GundamFitter) if( WITH_GUNDAM_ROOT_APP ) #target_sources( gundamRoot PRIVATE G__GundamRootDict.cxx ) #target_link_libraries( gundamRoot GundamPropagator ) From 419eaf86367aa34187aa165b491a89aa381026be Mon Sep 17 00:00:00 2001 From: lgiannes Date: Fri, 5 Apr 2024 14:20:41 +0200 Subject: [PATCH 122/131] 1.9.0 --- src/Applications/src/gundamCalcXsec.cxx | 3 +++ src/Applications/src/gundamMarginalise.cxx | 6 ++++-- 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/src/Applications/src/gundamCalcXsec.cxx b/src/Applications/src/gundamCalcXsec.cxx index 6ec400c9a..48ef4fc8e 100644 --- a/src/Applications/src/gundamCalcXsec.cxx +++ b/src/Applications/src/gundamCalcXsec.cxx @@ -665,6 +665,8 @@ int main(int argc, char** argv){ xsec.histogram.GetXaxis()->SetLabelSize(0.02); xsec.histogram.GetYaxis()->SetTitle( GenericToolbox::Json::fetchValue(xsec.samplePtr->getConfig(), "yAxis", "#delta#sigma").c_str() ); + + GenericToolbox::writeInTFile( GenericToolbox::mkdirTFile(calcXsecDir, "histograms"), &xsec.histogram, GenericToolbox::generateCleanBranchName( xsec.samplePtr->getName() ) @@ -743,6 +745,7 @@ int main(int argc, char** argv){ } } + LogInfo << "Generating canvases..." << std::endl; propagator.getPlotGenerator().generateCanvas( propagator.getPlotGenerator().getHistHolderList(0), GenericToolbox::mkdirTFile(calcXsecDir, "plots/canvas") diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index ddd8809ea..af0877d7a 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -527,9 +527,11 @@ int main(int argc, char** argv){ }// end if(injectParamsManually) - //LogInfo<<"Computing LH..."< Date: Fri, 5 Apr 2024 14:50:02 +0200 Subject: [PATCH 123/131] gundamMarginalise compliant with Gundam 1.9.0 - bug fix in link libraries --- src/Applications/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/src/Applications/CMakeLists.txt b/src/Applications/CMakeLists.txt index db6924f64..83c603dc1 100644 --- a/src/Applications/CMakeLists.txt +++ b/src/Applications/CMakeLists.txt @@ -81,6 +81,7 @@ target_link_libraries( gundamFitPlot GundamUtils ) target_link_libraries( gundamConfigUnfolder GundamUtils ) target_link_libraries( gundamConfigCompare GundamUtils ) target_link_libraries( gundamPlotExtractor GundamUtils ) +target_link_libraries( gundamMarginalise GundamPropagator GundamFitter GundamUtils ) if( WITH_GUNDAM_ROOT_APP ) #target_sources( gundamRoot PRIVATE G__GundamRootDict.cxx ) From 3882c6975839960139c2358e99af4be956f5b678 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Fri, 5 Apr 2024 15:53:12 +0200 Subject: [PATCH 124/131] reset parameters at the end of each LH computation --- src/Applications/src/gundamMarginalise.cxx | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index af0877d7a..979eeeec7 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -575,12 +575,13 @@ int main(int argc, char** argv){ if ( LLH-gLLH > log(weightCap)) { LogInfo << "Throw " << iToy << " over weight cap: LLH-gLLH = " << LLH - gLLH << std::endl; countBigThrows++; + }else{ + weightSum += LhOverGauss; + weightSquareSum += LhOverGauss*LhOverGauss; } //debug //LogInfo<<"LLH: "<1.e50){ weightSum /= 1.e50; weightSumE50++; @@ -593,6 +594,15 @@ int main(int argc, char** argv){ margThrowTree->Fill(); ThrowsPThetaFormat->Fill(); + // reset the parameters to the best fit values + for( auto& parSet : propagator.getParametersManager().getParameterSetsList() ){ + if( not parSet.isEnabled() ) continue; + for( auto& par : parSet.getParameterList() ){ + if( not par.isEnabled() ) continue; + par.setParameterValue(par.getPriorValue()); + } + } + }// end of main throws loop From 29e39511577d0f9fcc8045bda85ad777912ac135 Mon Sep 17 00:00:00 2001 From: lgiannes Date: Fri, 5 Apr 2024 16:01:28 +0200 Subject: [PATCH 125/131] set verbosity properly --- src/Applications/src/gundamMarginalise.cxx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 979eeeec7..51839fd69 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -527,11 +527,11 @@ int main(int argc, char** argv){ }// end if(injectParamsManually) - LogInfo<<"Computing LH... "; + // LogInfo<<"Computing LH... "; fitter.getLikelihoodInterface().propagateAndEvalLikelihood(); - LogInfo<<"Done. "; + // LogInfo<<"Done. "; LLH = fitter.getLikelihoodInterface().getBuffer().totalLikelihood; - LogInfo<<"LLH: "<1.e50){ weightSum /= 1.e50; From 3489e096a2c1524d906b25a91cec65975ee28331 Mon Sep 17 00:00:00 2001 From: Adrien Blanchet Date: Tue, 24 Sep 2024 15:05:39 +0200 Subject: [PATCH 126/131] restoring state of DialResponseSupervisor.cpp --- .../DialEngine/src/DialResponseSupervisor.cpp | 11 ++--------- 1 file changed, 2 insertions(+), 9 deletions(-) diff --git a/src/DialDictionary/DialEngine/src/DialResponseSupervisor.cpp b/src/DialDictionary/DialEngine/src/DialResponseSupervisor.cpp index 9d37a356a..321af24ac 100644 --- a/src/DialDictionary/DialEngine/src/DialResponseSupervisor.cpp +++ b/src/DialDictionary/DialEngine/src/DialResponseSupervisor.cpp @@ -3,21 +3,14 @@ // #include "DialResponseSupervisor.h" -#include "Logger.h" #include double DialResponseSupervisor::process(double reponse_) const { // apply cap? - // print out info - - if ( not std::isnan(_minResponse_) and reponse_ < _minResponse_ ){ - return _minResponse_; - } - else if( not std::isnan(_maxResponse_) and reponse_ > _maxResponse_ ){ - return _maxResponse_; - } + if ( not std::isnan(_minResponse_) and reponse_ < _minResponse_ ){ return _minResponse_; } + else if( not std::isnan(_maxResponse_) and reponse_ > _maxResponse_ ){ return _maxResponse_; } // else? From 168ee72ea6df793a585cd05a1d4e7e52e13a0314 Mon Sep 17 00:00:00 2001 From: Adrien Blanchet Date: Tue, 24 Sep 2024 15:05:51 +0200 Subject: [PATCH 127/131] link lib to the fitter --- src/Applications/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Applications/CMakeLists.txt b/src/Applications/CMakeLists.txt index 83c603dc1..26a93a5eb 100644 --- a/src/Applications/CMakeLists.txt +++ b/src/Applications/CMakeLists.txt @@ -74,6 +74,7 @@ endforeach() # Dependencies target_link_libraries( gundamFitter GundamFitter ) target_link_libraries( gundamCalcXsec GundamFitter ) # using the fitter engine to parse back the config file +target_link_libraries( gundamMarginalise GundamFitter ) target_link_libraries( gundamFitReader GundamUtils ) target_link_libraries( gundamInputZipper GundamUtils ) target_link_libraries( gundamFitCompare GundamUtils ) @@ -81,7 +82,6 @@ target_link_libraries( gundamFitPlot GundamUtils ) target_link_libraries( gundamConfigUnfolder GundamUtils ) target_link_libraries( gundamConfigCompare GundamUtils ) target_link_libraries( gundamPlotExtractor GundamUtils ) -target_link_libraries( gundamMarginalise GundamPropagator GundamFitter GundamUtils ) if( WITH_GUNDAM_ROOT_APP ) #target_sources( gundamRoot PRIVATE G__GundamRootDict.cxx ) From 1440120a45e0078b97c11253b87bf3f6f56b9203 Mon Sep 17 00:00:00 2001 From: Adrien Blanchet Date: Tue, 24 Sep 2024 15:08:18 +0200 Subject: [PATCH 128/131] resync generictoolbox --- submodules/cpp-generic-toolbox | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/submodules/cpp-generic-toolbox b/submodules/cpp-generic-toolbox index d4d0bdb8d..a161bb691 160000 --- a/submodules/cpp-generic-toolbox +++ b/submodules/cpp-generic-toolbox @@ -1 +1 @@ -Subproject commit d4d0bdb8d7583dc69377993a66287fb97065413f +Subproject commit a161bb6912843faf4d83ca20172a480c8a750c30 From bedd3f068328b770974a5ab35005e4115d856215 Mon Sep 17 00:00:00 2001 From: Adrien Blanchet Date: Tue, 24 Sep 2024 15:14:05 +0200 Subject: [PATCH 129/131] getting rid of outdated option --- .../src/ParametersManager.cpp | 24 +++++++------------ 1 file changed, 8 insertions(+), 16 deletions(-) diff --git a/src/ParametersManager/src/ParametersManager.cpp b/src/ParametersManager/src/ParametersManager.cpp index 3df964670..70386c144 100644 --- a/src/ParametersManager/src/ParametersManager.cpp +++ b/src/ParametersManager/src/ParametersManager.cpp @@ -376,24 +376,21 @@ void ParametersManager::throwParametersFromGlobalCovariance(std::vector ); weightsChiSquare.push_back(weights[iPar]); - if( _reThrowParSetIfOutOfBounds_ ){ - if( not _strippedParameterList_[iPar]->isValueWithinBounds() ){ - // re-do the throwing + if( not _strippedParameterList_[iPar]->isValueWithinBounds() ){ + // re-do the throwing // LogDebug << "Not within bounds: " << _strippedParameterList_[iPar]->getSummary() << std::endl; - rethrow = true; - } + rethrow = true; } } // Making sure eigen decomposed parameters get the conversion done for( auto& parSet : _parameterSetList_ ){ if( not parSet.isEnabled() ) continue; - if( parSet.isUseEigenDecompInFit() ){ + if( parSet.isEnableEigenDecomp() ){ parSet.propagateOriginalToEigen(); // also check the bounds of real parameter space - if( _reThrowParSetIfOutOfBounds_ ){ - for( auto& par : parSet.getEigenParameterList() ){ + for( auto& par : parSet.getEigenParameterList() ){ if( not par.isEnabled() ) continue; if( not par.isValueWithinBounds() ){ // re-do the throwing @@ -401,7 +398,6 @@ void ParametersManager::throwParametersFromGlobalCovariance(std::vector break; } } - } } } @@ -481,24 +477,21 @@ void ParametersManager::throwParametersFromTStudent(std::vector &weights ); weightsChiSquare.push_back(weights[iPar]); - if( _reThrowParSetIfOutOfBounds_ ){ - if( not _strippedParameterList_[iPar]->isValueWithinBounds() ){ + if( not _strippedParameterList_[iPar]->isValueWithinBounds() ){ // re-do the throwing // LogDebug << "Not within bounds: " << _strippedParameterList_[iPar]->getSummary() << std::endl; rethrow = true; } - } } // Making sure eigen decomposed parameters get the conversion done for( auto& parSet : _parameterSetList_ ){ if( not parSet.isEnabled() ) continue; - if( parSet.isUseEigenDecompInFit() ){ + if( parSet.isEnableEigenDecomp() ){ parSet.propagateOriginalToEigen(); // also check the bounds of real parameter space - if( _reThrowParSetIfOutOfBounds_ ){ - for( auto& par : parSet.getEigenParameterList() ){ + for( auto& par : parSet.getEigenParameterList() ){ if( not par.isEnabled() ) continue; if( not par.isValueWithinBounds() ){ // re-do the throwing @@ -506,7 +499,6 @@ void ParametersManager::throwParametersFromTStudent(std::vector &weights break; } } - } } } From d7c808a8478b32868452333eadd30b7af82876cb Mon Sep 17 00:00:00 2001 From: Adrien Blanchet Date: Tue, 24 Sep 2024 15:18:15 +0200 Subject: [PATCH 130/131] restoring header declarations --- src/ParametersManager/include/ParametersManager.h | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/ParametersManager/include/ParametersManager.h b/src/ParametersManager/include/ParametersManager.h index 0bcb62ea5..384c46d73 100644 --- a/src/ParametersManager/include/ParametersManager.h +++ b/src/ParametersManager/include/ParametersManager.h @@ -49,6 +49,9 @@ class ParametersManager : public JsonBaseClass { void throwParameters(); void throwParametersFromParSetCovariance(); void throwParametersFromGlobalCovariance(bool quietVerbose_ = true); + void throwParametersFromGlobalCovariance(std::vector &weightsChiSquare); + void throwParametersFromGlobalCovariance(std::vector &weightsChiSquare, double pedestalEntity, double pedestalLeftEdge, double pedestalRightEdge); + void throwParametersFromTStudent(std::vector &weightsChiSquare,double nu_); void initializeStrippedGlobalCov(); ParameterSet* getFitParameterSetPtr(const std::string& name_); From 64e8eb9c2c58ef9b369f02ae14e26ab09ecb3d5e Mon Sep 17 00:00:00 2001 From: Adrien Blanchet Date: Tue, 24 Sep 2024 15:29:13 +0200 Subject: [PATCH 131/131] adapting with the new structure of the propagator --- src/Applications/src/gundamMarginalise.cxx | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/src/Applications/src/gundamMarginalise.cxx b/src/Applications/src/gundamMarginalise.cxx index 51839fd69..86dfd99e7 100644 --- a/src/Applications/src/gundamMarginalise.cxx +++ b/src/Applications/src/gundamMarginalise.cxx @@ -102,8 +102,8 @@ int main(int argc, char** argv){ gRandom->SetSeed(seed); } - GundamGlobals::getParallelWorker().setNThreads(clParser.getOptionVal("nbThreads", 1)); - LogInfo << "Running the fitter with " << GundamGlobals::getParallelWorker().getNbThreads() << " parallel threads." << std::endl; + GundamGlobals::setNumberOfThreads(clParser.getOptionVal("nbThreads", 1)); + LogInfo << "Running the fitter with " << GundamGlobals::getNumberOfThreads() << " parallel threads." << std::endl; // Reading fitter file LogInfo << "Opening fitter output file: " << clParser.getOptionVal("fitterOutputFile") << std::endl; @@ -198,24 +198,22 @@ int main(int argc, char** argv){ tStudent = true; } - auto configPropagator = GenericToolbox::Json::fetchValuePath( cHandler.getConfig(), "fitterEngineConfig/propagatorConfig" ); + auto configPropagator = GenericToolbox::Json::fetchValue( cHandler.getConfig(), "fitterEngineConfig/propagatorConfig" ); // Initialize the fitterEngine LogInfo << "FitterEngine setup..." << std::endl; FitterEngine fitter{nullptr}; - fitter.readConfig(GenericToolbox::Json::fetchSubEntry((JsonType)cHandler.getConfig(), {"fitterEngineConfig"})); - - DataSetManager& dataSetManager{fitter.getLikelihoodInterface().getDataSetManager()}; + fitter.readConfig(GenericToolbox::Json::fetchValue( cHandler.getConfig(), "fitterEngineConfig")); // We are only interested in our MC. Data has already been used to get the post-fit error/values - dataSetManager.getPropagator().setLoadAsimovData( true ); + fitter.getLikelihoodInterface().setForceAsimovData( true ); // Disabling eigen decomposed parameters - dataSetManager.getPropagator().setEnableEigenToOrigInPropagate( false ); + fitter.getLikelihoodInterface().getModelPropagator().setEnableEigenToOrigInPropagate( false ); // Load everything - dataSetManager.initialize(); - Propagator& propagator{dataSetManager.getPropagator()}; + fitter.getLikelihoodInterface().initialize(); + Propagator& propagator{fitter.getLikelihoodInterface().getModelPropagator()}; propagator.getParametersManager().getParameterSetsList();