From 706c4d8895a1458164feed5a31de6924a19d4dcf Mon Sep 17 00:00:00 2001 From: Jaesung Kim Date: Fri, 2 Jun 2023 12:25:54 -0500 Subject: [PATCH 1/5] Update SystToolsEventWeight for dependent dials --- .../App/SystToolsEventWeight_module.cc | 102 +++++++++++++++--- 1 file changed, 87 insertions(+), 15 deletions(-) diff --git a/sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc b/sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc index 4793171b9..78fee3c44 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,12 +139,13 @@ 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; @@ -159,8 +163,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,24 +195,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(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; + } + } + } + for( auto &sph : smd ){ if (fDebugMode) { std::cout << " sph.prettyName = " << sph.prettyName << std::endl; } - std::string rwmode = "multisigma"; - if(sph.isRandomlyThrown) rwmode = "multisim"; + // responsless + if(sph.isResponselessParam){ + if (fDebugMode) { + std::cout << "Responsless dial found: " << sph.prettyName << ", thus skipping" << std::endl; + } + continue; + } + + sbn::evwgh::EventWeightParameterSet fParameterSet; + std::string rwmode = ""; + + 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); + + 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 { sph_dep.paramVariations.begin(), sph_dep.paramVariations.end() }; + fParameterSet.AddParameter(sph_dep.prettyName, std::move(widths_dep)); + + } + + + + } + else{ + + rwmode = sph.isRandomlyThrown ? "multisim" : "multisigma"; + + std::vector widths { sph.paramVariations.begin(), sph.paramVariations.end() }; + + sbn::evwgh::EventWeightParameterSet fParameterSet; + fParameterSet.AddParameter(sph.prettyName, std::move(widths)); + + } + + if(rwmode==""){ + throw cet::exception{ "SystToolsEventWeight" } + << "rwmode not set for " << sph.prettyName << "\n"; + } if (fDebugMode) { std::cout << " rwmode = " << rwmode << std::endl; } - 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()); fParameterSet.FillKnobValues(); From f0c3c49448d9aeaadc0df911305103d5bc2a3556 Mon Sep 17 00:00:00 2001 From: Jaesung Kim Date: Mon, 12 Jun 2023 09:38:25 -0500 Subject: [PATCH 2/5] Duplicated declaration of fParameterSet removed --- sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc b/sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc index 78fee3c44..ed3c870c0 100644 --- a/sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc +++ b/sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc @@ -149,6 +149,8 @@ void SystToolsEventWeight::produce(art::Event& e) { 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" @@ -267,7 +269,6 @@ void SystToolsEventWeight::beginRun(art::Run& run) { std::vector widths { sph.paramVariations.begin(), sph.paramVariations.end() }; - sbn::evwgh::EventWeightParameterSet fParameterSet; fParameterSet.AddParameter(sph.prettyName, std::move(widths)); } From 12aa7494968f529df72a2e4e9f740b07d2ee474b Mon Sep 17 00:00:00 2001 From: Jaesung Kim Date: Sat, 24 Jun 2023 22:19:53 -0500 Subject: [PATCH 3/5] Update logger for SystToolsEventWeight_module.cc --- sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc b/sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc index ed3c870c0..62f03b6ce 100644 --- a/sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc +++ b/sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc @@ -226,9 +226,6 @@ void SystToolsEventWeight::beginRun(art::Run& run) { } for( auto &sph : smd ){ - if (fDebugMode) { - std::cout << " sph.prettyName = " << sph.prettyName << std::endl; - } // responsless if(sph.isResponselessParam){ @@ -279,7 +276,7 @@ void SystToolsEventWeight::beginRun(art::Run& run) { } if (fDebugMode) { - std::cout << " rwmode = " << rwmode << std::endl; + std::cout << " sph.prettyName = " << sph.prettyName << ", rwmode = " << rwmode << " is added to the header" << std::endl; } fParameterSet.Configure(sp->GetFullyQualifiedName()+"_"+sph.prettyName, rwmode, sph.paramVariations.size()); From 58c6946968495db9184c2d4777d64fe361436bd0 Mon Sep 17 00:00:00 2001 From: Jaesung Kim Date: Thu, 27 Jul 2023 13:24:07 -0500 Subject: [PATCH 4/5] Update SystToolsEventWeight module to treat isCorrection=true paramheader --- .../App/SystToolsEventWeight_module.cc | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) diff --git a/sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc b/sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc index 62f03b6ce..6e836054c 100644 --- a/sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc +++ b/sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc @@ -157,7 +157,14 @@ void SystToolsEventWeight::produce(art::Event& e) { << "[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 " @@ -238,11 +245,14 @@ void SystToolsEventWeight::beginRun(art::Run& run) { sbn::evwgh::EventWeightParameterSet fParameterSet; 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{ @@ -252,7 +262,7 @@ void SystToolsEventWeight::beginRun(art::Run& run) { } } - std::vector widths_dep { sph_dep.paramVariations.begin(), sph_dep.paramVariations.end() }; + std::vector widths_dep { paramVars_dep.begin(), paramVars_dep.end() }; fParameterSet.AddParameter(sph_dep.prettyName, std::move(widths_dep)); } @@ -264,7 +274,7 @@ void SystToolsEventWeight::beginRun(art::Run& run) { rwmode = sph.isRandomlyThrown ? "multisim" : "multisigma"; - std::vector widths { sph.paramVariations.begin(), sph.paramVariations.end() }; + std::vector widths { paramVars.begin(), paramVars.end() }; fParameterSet.AddParameter(sph.prettyName, std::move(widths)); @@ -279,7 +289,7 @@ void SystToolsEventWeight::beginRun(art::Run& run) { std::cout << " sph.prettyName = " << sph.prettyName << ", rwmode = " << rwmode << " is added to the header" << std::endl; } - fParameterSet.Configure(sp->GetFullyQualifiedName()+"_"+sph.prettyName, rwmode, sph.paramVariations.size()); + fParameterSet.Configure(sp->GetFullyQualifiedName()+"_"+sph.prettyName, rwmode, paramVars.size()); fParameterSet.FillKnobValues(); p->push_back( std::move(fParameterSet) ); From 43893b9b84e91cf3f2aa6676f1ab7cded4317d8a Mon Sep 17 00:00:00 2001 From: Jaesung Kim Date: Thu, 15 Feb 2024 22:53:55 -0600 Subject: [PATCH 5/5] Update systematicstools/nusystematics versions --- ups/product_deps | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ups/product_deps b/ups/product_deps index 8d2b04768..3c0df5cca 100644 --- a/ups/product_deps +++ b/ups/product_deps @@ -259,8 +259,8 @@ sbnanaobj v09_20_06 - sbndaq_artdaq_core v1_06_00of0 - sbndata v01_04 - sbnobj v09_16_00 - -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 ####################################