diff --git a/bin/PROfit_setup.cxx b/bin/PROfit_setup.cxx index 164d6d5c..36e02514 100644 --- a/bin/PROfit_setup.cxx +++ b/bin/PROfit_setup.cxx @@ -14,6 +14,52 @@ static int find_spline_index(const PROsyst &systs, const PROconfig &config, cons return it == systs.spline_names.end() ? -1 : (int)std::distance(systs.spline_names.begin(), it); } +// File-name tag for an injected point, e.g. "dmsq3_sinsq2thmm0.25". Values use the short form +// unless it doesn't round-trip, so distinct points never share a cached binary. +static std::string inject_tag(const PROpt &options) { + std::string tag; + for(const auto &[name, value]: options.fake_data_osc_params) { // std::map: sorted by name + std::ostringstream v; v << value; + if(std::stof(v.str()) != value) { v.str(""); v << std::setprecision(std::numeric_limits::max_digits10) << value; } + tag += (tag.empty() ? "" : "_") + name + v.str(); + } + for(char &c: tag) if(!std::isalnum((unsigned char)c) && c != '.' && c != '-' && c != '_') c = '_'; + return tag; +} + +// --data + --inject with a data : process data like MC (no systs) and return the +// histograms at the injected point, ready to be cached like CreatePROdata's output. +static std::vector CreateInjectedPROdata(const PROconfig &dataconfig, const PROmodel &model, + const Eigen::VectorXf &fakedataparams, const PROpt &options) { + PROpeller dataprop; + std::vector> nosysts; + PROcess_CAFAna(dataconfig, nosysts, dataprop, options.noxrootd); + std::unique_ptr datamodel = get_model_from_string(dataconfig, dataprop); + for(const auto &[name, value]: options.fake_data_osc_params) { + if(std::find(datamodel->param_names.begin(), datamodel->param_names.end(), name) == datamodel->param_names.end()) { + log(L"%1% || --inject parameter %2% is not in the data model %3%.") % __func__ % name.c_str() % datamodel->param_names; + exit(EXIT_FAILURE); + } + } + // Map injected values by name; un-injected data-model params stay at their defaults. + Eigen::VectorXf dataparams = datamodel->default_val; + for(size_t j = 0; j < datamodel->nparams; ++j) { + auto it = std::find(model.param_names.begin(), model.param_names.end(), datamodel->param_names[j]); + if(it == model.param_names.end()) continue; + size_t k = std::distance(model.param_names.begin(), it); + float lin = model.is_log10[k] ? std::pow(10.0f, fakedataparams(k)) : fakedataparams(k); + dataparams(j) = datamodel->is_log10[j] ? std::log10(lin) : lin; + } + log(L"%1% || Filling data (%2% events) at injected point %3%") % __func__ % dataprop.NEvent() % dataparams; + PROsyst emptySyst; + std::vector alldata; + for(size_t io = 0; io < dataconfig.m_num_variables; ++io) { + Eigen::VectorXf vec = CollapseMatrix(dataconfig, FillSpectra(dataconfig, dataprop, emptySyst, *datamodel, dataparams, !options.eventbyevent, io).Spec(), io); + alldata.push_back(PROdata(vec, vec.array().sqrt())); + } + return alldata; +} + PROdata construct_data(std::vector &variable_data, bool use_real_data, const PROconfig &config, const PROpeller &prop, const PROmodel &model, const std::vector &variable_systs, const Eigen::VectorXf &fakedataparams, const PROpt &options) { PROdata data; std::uniform_int_distribution dseed(0, std::numeric_limits::max()); @@ -53,11 +99,16 @@ PROdata construct_data(std::vector &variable_data, bool use_real_data, exit(EXIT_FAILURE); } + const bool inject = !options.fake_data_osc_params.empty() && dataconfig.m_has_model; + if(!options.fake_data_osc_params.empty() && !dataconfig.m_has_model) + log(L"%1% || --inject is ignored: the data config has no .") % __func__; + // The injected point is in the name, so a binary is only reused at the same point. + if(inject) dataBinName = options.analysis_tag + "_data_inject_" + inject_tag(options) + ".bin"; if((*options.process_command) || (!std::filesystem::exists(dataBinName)) ){ log(L"%1% || Processing Data Spectrum and saving to binary output also: %2%") % __func__ % dataBinName.c_str(); //Process the CAF files to grab and fill spectrum directly - std::vector alldata = CreatePROdata(dataconfig); + std::vector alldata = inject ? CreateInjectedPROdata(dataconfig, model, fakedataparams, options) : CreatePROdata(dataconfig); PROdata::saveVector(dataconfig, alldata, dataBinName); data = alldata[config.i_prime]; //data.save(dataconfig,dataBinName); @@ -91,6 +142,12 @@ PROdata construct_data(std::vector &variable_data, bool use_real_data, } } + // Cache the Asimov injected spectrum; any Poisson throw is applied fresh each run. + if(inject && options.poisson_throw) + for(auto &d: variable_data) { + PROspec thrown = PROspec::PoissonVariation(PROspec(d.Spec(), d.Error()), dseed(PROseed::global_rng)); + d = PROdata(thrown.Spec(), thrown.Spec().array().sqrt()); + } }//if no data, use injected or fake data; else{ log(L"%1% || Going to get fake data set up for each variable.") % __func__ ; diff --git a/inc/PROconfig.h b/inc/PROconfig.h index 2fd499e8..2583bac9 100644 --- a/inc/PROconfig.h +++ b/inc/PROconfig.h @@ -581,6 +581,7 @@ namespace PROfit{ //Some model infomation std::string m_model_tag; + bool m_has_model = false; ///< A was parsed; lets --inject reweight a data config. std::vector m_model_rule_index; std::vector m_model_rule_names; std::vector m_model_parameter_index; diff --git a/src/PROconfig.cxx b/src/PROconfig.cxx index ceb2e03e..402dbfdb 100644 --- a/src/PROconfig.cxx +++ b/src/PROconfig.cxx @@ -2166,6 +2166,7 @@ int PROconfig::LoadFromXML(const std::string &filename){ pModel = doc.FirstChildElement("model"); if(pModel){ + m_has_model = true; // Read in how many bins this channel uses const std::vector expected_attrs = {"tag","name","index"}; @@ -2441,6 +2442,18 @@ int PROconfig::LoadFromXML(const std::string &filename){ } dataXml << "\n"; + // Opt-in data model (for --inject): a inside , or inherit_model="true". + tinyxml2::XMLElement* pDataModel = pData->FirstChildElement("model"); + if(!pDataModel && pData->BoolAttribute("inherit_model", false)) { + pDataModel = doc.FirstChildElement("model"); + if(!pDataModel) throw std::invalid_argument(" requires a top-level "); + } + if(pDataModel) { + tinyxml2::XMLPrinter printer; + pDataModel->Accept(&printer); + dataXml << printer.CStr() << "\n"; + } + // Serialize MCFile elements from inside the section tinyxml2::XMLElement* pDataMC = pData->FirstChildElement("MCFile"); while(pDataMC) {