diff --git a/sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc b/sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc index 4793171b9..6e836054c 100644 --- a/sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc +++ b/sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc @@ -116,19 +116,22 @@ void SystToolsEventWeight::produce(art::Event& e) { << sp->GetFullyQualifiedName() << "\n"; } - //==== syst_resp->size() is (Number of MCTruth) - //==== Each index corresponds to each of MCTruth + // syst_resp->size() is (Number of MCTruth) + // Each index corresponds to each of MCTruth int nMCTruthIndex = 0; if(fDebugMode) std::cout << "[SystToolsEventWeight::produce] syst_resp.size() (= Number of MCTruth) of this SystProvider = " << syst_resp->size() << std::endl; - //==== Looping over syst_resp is identical to looping over MCTruth + // Looping over syst_resp is identical to looping over MCTruth for(systtools::event_unit_response_t const& resp: *syst_resp) { - //==== resp.size() corresponds to number of knobs we altered; - //==== e.g., MaCCQE, MaCCRES, MvCCRE -> resp.size() = 3 + // resp.size() corresponds to number of knobs we altered; + // e.g., MaCCQE, MaCCRES, MvCCRE -> resp.size() = 3 if(fDebugMode){ std::cout << "[SystToolsEventWeight::produce] sp->GetSystMetaData().size() (expected) = " << sp->GetSystMetaData().size() << "\n" << "[SystToolsEventWeight::produce] resp.size() of this syst_resp (produced) = " << resp.size() << std::endl; } + + // Below check is not valid if there is a dependent dial + /* if(sp->GetSystMetaData().size()!=resp.size()){ throw cet::exception{ "SystToolsEventWeight" } << "sp->GetFullyQualifiedName() = " << sp->GetFullyQualifiedName() << "\n" @@ -136,22 +139,32 @@ void SystToolsEventWeight::produce(art::Event& e) { << resp.size() << " are produced. " << "Probably this particular event is not relevant to this systematic variation.\n"; } + */ sbn::evwgh::EventWeightMap mcwgh; for( auto const& r: resp){ - //==== responses.size() is the number of universes + // responses.size() is the number of universes systtools::SystParamHeader const &sph = fParamHeaderHelper.GetHeader( r.pid ); std::string prettyName = sph.prettyName; + if(sph.isResponselessParam) continue; + if(fDebugMode){ std::cout << "[SystToolsEventWeight::produce] pid of this resp = " << r.pid << "\n" << "[SystToolsEventWeight::produce] prettyName of this resp = " << prettyName << "\n" << "[SystToolsEventWeight::produce] paramVariations.size() of this resp (expected) = " << sph.paramVariations.size() << "\n" << "[SystToolsEventWeight::produce] responses.size() of this resp (produced) = " << r.responses.size() << std::endl; } - if(sph.paramVariations.size()!=r.responses.size()){ + + if(sph.isCorrection && (r.responses.size()!=1)){ + throw cet::exception{ "SystToolsEventWeight" } + << "prettyName of this resp = " << prettyName << "\n" + << "We expect to have 1 universes as it is a correction, but " + << r.responses.size() << " are produced\n"; + } + if(!sph.isCorrection && sph.paramVariations.size()!=r.responses.size()){ throw cet::exception{ "SystToolsEventWeight" } << "prettyName of this resp = " << prettyName << "\n" << "We expect to have " << sph.paramVariations.size() << " universes, but " @@ -159,8 +172,8 @@ void SystToolsEventWeight::produce(art::Event& e) { << "Probably this particular event is not relevant to this systematic variation.\n"; } - //==== r.responses : std::vector - //==== std::map > EventWeightMap + // r.responses : std::vector + // std::map > EventWeightMap mcwgh[sp->GetFullyQualifiedName()+"_"+prettyName].assign(r.responses.cbegin(), r.responses.cend()); } // END loop over knobs @@ -191,25 +204,92 @@ void SystToolsEventWeight::beginRun(art::Run& run) { << "sp->GetFullyQualifiedName() = " << sp->GetFullyQualifiedName() << "\n" << "sp->GetInstanceName() = " << sp->GetInstanceName() << "\n" << "Printing each SystParamHeader of this ISystProviderTool.."; - //==== Note: typedef std::vector SystMetaData; + // Note: typedef std::vector SystMetaData; auto const& smd = sp->GetSystMetaData(); + + // make a map of responsless-response params + std::map> map_resp_to_respless; for( auto &sph : smd ){ - if (fDebugMode) { - std::cout << " sph.prettyName = " << sph.prettyName << std::endl; + if(sph.isResponselessParam){ + auto it = map_resp_to_respless.find( sph.responseParamId ); + if( it != map_resp_to_respless.end() ){ + it->second.push_back( sph.systParamId ); + } + else{ + map_resp_to_respless[sph.responseParamId] = {}; + map_resp_to_respless[sph.responseParamId].push_back( sph.systParamId ); + } } + } + if (fDebugMode) { + for(const auto& it: map_resp_to_respless){ + const auto& sph = systtools::GetParam(smd, it.first); + std::cout << "Found a dependent dial: " << sph.prettyName << std::endl; + for(const auto& depdialid: it.second){ + const auto& sph_dep = systtools::GetParam(smd, depdialid); + std::cout << "- dep dial: " << sph_dep.prettyName << std::endl; + } + } + } - std::string rwmode = "multisigma"; - if(sph.isRandomlyThrown) rwmode = "multisim"; + for( auto &sph : smd ){ - if (fDebugMode) { - std::cout << " rwmode = " << rwmode << std::endl; + // responsless + if(sph.isResponselessParam){ + if (fDebugMode) { + std::cout << "Responsless dial found: " << sph.prettyName << ", thus skipping" << std::endl; + } + continue; } - std::vector widths { sph.paramVariations.begin(), sph.paramVariations.end() }; - sbn::evwgh::EventWeightParameterSet fParameterSet; - fParameterSet.AddParameter(sph.prettyName, std::move(widths)); - fParameterSet.Configure(sp->GetFullyQualifiedName()+"_"+sph.prettyName, rwmode, sph.paramVariations.size()); + std::string rwmode = ""; + + std::vector paramVars = sph.isCorrection ? std::vector(1, sph.centralParamValue) : sph.paramVariations; + + auto it = map_resp_to_respless.find( sph.systParamId ); + if(it!=map_resp_to_respless.end()){ + + for(const auto depdialid: it->second){ + const auto& sph_dep = systtools::GetParam(smd, depdialid); + std::vector paramVars_dep = sph_dep.isCorrection ? std::vector(1, sph_dep.centralParamValue) : sph_dep.paramVariations; + + if(rwmode=="") rwmode = sph_dep.isRandomlyThrown ? "multisim" : "multisigma"; + else{ + if(rwmode!= (sph_dep.isRandomlyThrown ? "multisim" : "multisigma")){ + throw cet::exception{ "SystToolsEventWeight" } + << sph.prettyName << " depends on other dials, but the remode are different between the deps dials\n"; + } + } + + std::vector widths_dep { paramVars_dep.begin(), paramVars_dep.end() }; + fParameterSet.AddParameter(sph_dep.prettyName, std::move(widths_dep)); + + } + + + + } + else{ + + rwmode = sph.isRandomlyThrown ? "multisim" : "multisigma"; + + std::vector widths { paramVars.begin(), paramVars.end() }; + + fParameterSet.AddParameter(sph.prettyName, std::move(widths)); + + } + + if(rwmode==""){ + throw cet::exception{ "SystToolsEventWeight" } + << "rwmode not set for " << sph.prettyName << "\n"; + } + + if (fDebugMode) { + std::cout << " sph.prettyName = " << sph.prettyName << ", rwmode = " << rwmode << " is added to the header" << std::endl; + } + + fParameterSet.Configure(sp->GetFullyQualifiedName()+"_"+sph.prettyName, rwmode, paramVars.size()); fParameterSet.FillKnobValues(); p->push_back( std::move(fParameterSet) ); diff --git a/ups/product_deps b/ups/product_deps index 8894b965c..0053be974 100644 --- a/ups/product_deps +++ b/ups/product_deps @@ -259,8 +259,8 @@ sbnanaobj v09_20_06_03 - sbndaq_artdaq_core v1_06_00of0 - sbndata v01_04 - sbnobj v09_16_00_02 - -systematicstools v01_02_00 - -nusystematics v01_02_10 - +systematicstools v01_02_01 - +nusystematics v01_02_11 - cetmodules v3_20_00 - only_for_build end_product_list ####################################