From f817b6deeaf72c66731086f8356e991ed2fd52e9 Mon Sep 17 00:00:00 2001 From: parkerw Date: Tue, 12 Nov 2024 21:26:42 +0000 Subject: [PATCH 1/6] add unit test that checks re-evaluations of llh aren't a problem if using a buffer region --- test/unit/BinnedNLLHTest.cpp | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/test/unit/BinnedNLLHTest.cpp b/test/unit/BinnedNLLHTest.cpp index 12a134f..7fc0912 100644 --- a/test/unit/BinnedNLLHTest.cpp +++ b/test/unit/BinnedNLLHTest.cpp @@ -101,6 +101,30 @@ TEST_CASE("Binned NLLH, 3 rates no systematics") REQUIRE_THAT(lh.Evaluate(), Catch::Matchers::WithinAbs(sumNorm + sumLogProb + constraint, 0.0001)); } + SECTION("Consistent Probability with buffer") + { + + BinnedNLLH lh_buffer; + lh_buffer.AddPdf(pdf1); + + // Add a generous buffer region + lh_buffer.SetBuffer("axis1",30,30); + lh_buffer.SetDataSet(&data); + lh_buffer.RegisterFitComponents(); + + ParameterDict params; + params["a"] = 1; + lh_buffer.SetParameters(params); + double llh1 = lh_buffer.Evaluate(); + + // Now let's re-evaluate, to check there's no issue with the shrinking being reapplied + double llh2 = lh_buffer.Evaluate(); + REQUIRE_THAT(llh1, Catch::Matchers::WithinAbs(llh2, 0.0001)); + + std::pair buffers = lh_buffer.GetBuffer("axis1"); + REQUIRE(buffers == std::pair(30,30)); + } + std::vector genRates(3, pow(10, 6)); BinnedNLLH lh2; lh2.SetBarlowBeeston(true); From 1fecdcd10f26cb205681f2bde2880788aa0e565d Mon Sep 17 00:00:00 2001 From: parkerw Date: Tue, 12 Nov 2024 21:27:48 +0000 Subject: [PATCH 2/6] don't re-shrink if using Direct or False norm fitting --- src/fitutil/BinnedEDManager.cpp | 36 +++++++++++++++++++++------------ 1 file changed, 23 insertions(+), 13 deletions(-) diff --git a/src/fitutil/BinnedEDManager.cpp b/src/fitutil/BinnedEDManager.cpp index e0fe0af..33b2bb3 100644 --- a/src/fitutil/BinnedEDManager.cpp +++ b/src/fitutil/BinnedEDManager.cpp @@ -144,21 +144,31 @@ BinnedEDManager::ApplyShrink(const BinnedEDShrinker& shrinker_){ for (size_t i = 0; i < fWorkingPdfs.size(); i++){ // Normalise if normalisation is a fittable param, but if indirect then track any change if (fAllowNormsFittable.at(i) == DIRECT) { - fWorkingPdfs[i] = shrinker_.ShrinkDist(fWorkingPdfs.at(i)); - fWorkingPdfs[i].Normalise(); + if( fWorkingPdfs.at(0).GetNBins() == fOriginalPdfs.at(0).GetNBins() ){ + fWorkingPdfs[i] = shrinker_.ShrinkDist(fWorkingPdfs.at(i)); + fWorkingPdfs[i].Normalise(); + } + else{ + return; + } } else if (fAllowNormsFittable.at(i) == FALSE) { - fWorkingPdfs[i] = shrinker_.ShrinkDist(fWorkingPdfs.at(i)); - const double integral_after = fWorkingPdfs[i].Integral(); - fNormalisations[i] = integral_after; - } else { - const double integral_before = fWorkingPdfs[i].Integral(); - fWorkingPdfs[i] = shrinker_.ShrinkDist(fWorkingPdfs.at(i)); - const double integral_after = fWorkingPdfs[i].Integral(); - if (integral_before == 0. && integral_after == 0.) { fNormalisations[i] = 0.; } - else { fNormalisations[i] *= integral_after/integral_before; } - } + if ( fWorkingPdfs.at(0).GetNBins() == fOriginalPdfs.at(0).GetNBins() ){ + fWorkingPdfs[i] = shrinker_.ShrinkDist(fWorkingPdfs.at(i)); + const double integral_after = fWorkingPdfs[i].Integral(); + fNormalisations[i] = integral_after; + } + else{ + return; + } + } else { + const double integral_before = fWorkingPdfs[i].Integral(); + fWorkingPdfs[i] = shrinker_.ShrinkDist(fWorkingPdfs.at(i)); + const double integral_after = fWorkingPdfs[i].Integral(); + if (integral_before == 0. && integral_after == 0.) { fNormalisations[i] = 0.; } + else { fNormalisations[i] *= integral_after/integral_before; } + } } - + } //////////////////////////////// From 700fef4f61f1e3031756a41a7437b6b5a5a5d68e Mon Sep 17 00:00:00 2001 From: parkerw Date: Tue, 10 Dec 2024 16:04:08 +0000 Subject: [PATCH 3/6] In SystMan AddPDFs, if using syst groups call corresponding AddPDF method --- src/teststat/BinnedNLLH.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/teststat/BinnedNLLH.cpp b/src/teststat/BinnedNLLH.cpp index 603b75e..c1deaaf 100644 --- a/src/teststat/BinnedNLLH.cpp +++ b/src/teststat/BinnedNLLH.cpp @@ -215,11 +215,11 @@ void BinnedNLLH::AddPdfs(const std::vector &pdfs, { if (norm_fitting_statuses != nullptr) { - AddPdf(pdfs.at(i), genrates_.at(i), norm_fitting_statuses->at(i)); + AddPdf(pdfs.at(i), sys_.at(i), genrates_.at(i), norm_fitting_statuses->at(i)); } else { - AddPdf(pdfs.at(i), genrates_.at(i)); + AddPdf(pdfs.at(i), sys_.at(i), genrates_.at(i)); } } } From fba48540a5cd07ce238aab46ce4ef513135059d2 Mon Sep 17 00:00:00 2001 From: parkerw Date: Tue, 10 Dec 2024 16:05:29 +0000 Subject: [PATCH 4/6] binned ed manager formatting --- src/fitutil/BinnedEDManager.cpp | 269 ++++++++++++++++++++------------ 1 file changed, 170 insertions(+), 99 deletions(-) diff --git a/src/fitutil/BinnedEDManager.cpp b/src/fitutil/BinnedEDManager.cpp index 33b2bb3..b3626ad 100644 --- a/src/fitutil/BinnedEDManager.cpp +++ b/src/fitutil/BinnedEDManager.cpp @@ -5,32 +5,41 @@ #include #include -unsigned -BinnedEDManager::GetNPdfs() const{ +unsigned +BinnedEDManager::GetNPdfs() const +{ return fOriginalPdfs.size(); } size_t -BinnedEDManager::GetNDims() const{ +BinnedEDManager::GetNDims() const +{ return fNDims; } -double -BinnedEDManager::Probability(const Event& data_) { +double +BinnedEDManager::Probability(const Event &data_) +{ ReassertNorms(true); double sum = 0; size_t j = 0; - for(size_t i = 0; i < fWorkingPdfs.size(); i++){ - if (fAllowNormsFittable.at(i) == DIRECT) { + for (size_t i = 0; i < fWorkingPdfs.size(); i++) + { + if (fAllowNormsFittable.at(i) == DIRECT) + { sum += fNormalisations.at(i) * fWorkingPdfs[i].Probability(data_); j++; - } else if (fAllowNormsFittable.at(i) == INDIRECT) { + } + else if (fAllowNormsFittable.at(i) == INDIRECT) + { sum += fFittableNorms.at(j) * fWorkingPdfs.at(i).Probability(data_); j++; - } else { - // In case where we don't auto-normalise the pdf, - // the normalisation is held within the pdf itself! - sum += fWorkingPdfs.at(i).Probability(data_); + } + else + { + // In case where we don't auto-normalise the pdf, + // the normalisation is held within the pdf itself! + sum += fWorkingPdfs.at(i).Probability(data_); } } @@ -38,137 +47,181 @@ BinnedEDManager::Probability(const Event& data_) { } double -BinnedEDManager::BinProbability(size_t bin_) { +BinnedEDManager::BinProbability(size_t bin_) +{ ReassertNorms(true); double sum = 0; size_t j = 0; - try{ - for(size_t i = 0; i < fWorkingPdfs.size(); i++){ - if (fAllowNormsFittable.at(i) == DIRECT) { + try + { + for (size_t i = 0; i < fWorkingPdfs.size(); i++) + { + if (fAllowNormsFittable.at(i) == DIRECT) + { sum += fNormalisations.at(i) * fWorkingPdfs.at(i).GetBinContent(bin_); j++; - } else if (fAllowNormsFittable.at(i) == INDIRECT) { + } + else if (fAllowNormsFittable.at(i) == INDIRECT) + { sum += fFittableNorms.at(j) * fWorkingPdfs.at(i).GetBinContent(bin_); j++; - } else { + } + else + { // In case where we don't auto-normalise the pdf, // the normalisation is held within the pdf itself! sum += fWorkingPdfs.at(i).GetBinContent(bin_); } } } - catch(const std::out_of_range&){ + catch (const std::out_of_range &) + { throw LogicError("BinnedEDManager:: Normalisation vector doesn't match pdf vector - are the normalisations set?"); } return sum; } - -void -BinnedEDManager::SetNormalisations(const std::vector& normalisations_){ +void BinnedEDManager::SetNormalisations(const std::vector &normalisations_) +{ if (normalisations_.size() != fOriginalPdfs.size()) throw LogicError("BinnedEDManager: number of norms doesn't match #pdfs"); fNormalisations = normalisations_; } -void -BinnedEDManager::ApplySystematics(const SystematicManager& sysMan_){ +void BinnedEDManager::ApplySystematics(const SystematicManager &sysMan_) +{ // If there are no systematics dont do anything // ( working pdfs = original pdfs from initialisation) - if(!sysMan_.GetNSystematics()) + if (!sysMan_.GetNSystematics()) return; - if (fAllNormsDirFittable) { - sysMan_.DistortEDs(fOriginalPdfs,fWorkingPdfs); - } else { + if (fAllNormsDirFittable) + { + sysMan_.DistortEDs(fOriginalPdfs, fWorkingPdfs); + } + else + { std::vector norms(fNormalisations.size(), 0); - sysMan_.DistortEDs(fOriginalPdfs,fWorkingPdfs, &norms); - for (size_t i = 0; i < norms.size(); i++) { - if (fAllowNormsFittable.at(i) == FALSE) { + sysMan_.DistortEDs(fOriginalPdfs, fWorkingPdfs, &norms); + for (size_t i = 0; i < norms.size(); i++) + { + if (fAllowNormsFittable.at(i) == FALSE) + { fNormalisations[i] = norms.at(i); - } else if (fAllowNormsFittable.at(i) == INDIRECT) { + } + else if (fAllowNormsFittable.at(i) == INDIRECT) + { fNormalisations[i] *= norms.at(i); } } } } -const BinnedED& -BinnedEDManager::GetOriginalPdf(size_t index_) const{ +const BinnedED & +BinnedEDManager::GetOriginalPdf(size_t index_) const +{ return fOriginalPdfs.at(index_); } -void -BinnedEDManager::AddPdf(const BinnedED& pdf_, const NormFittingStatus norm_fitting_status){ +void BinnedEDManager::AddPdf(const BinnedED &pdf_, const NormFittingStatus norm_fitting_status) +{ fOriginalPdfs.push_back(pdf_); fWorkingPdfs.push_back(pdf_); fNPdfs++; fNormalisations.resize(fOriginalPdfs.size(), 0); fAllowNormsFittable.push_back(norm_fitting_status); - if (norm_fitting_status) { fFittableNorms.push_back(0); } - if (norm_fitting_status != DIRECT) { fAllNormsDirFittable = false; } + if (norm_fitting_status) + { + fFittableNorms.push_back(0); + } + if (norm_fitting_status != DIRECT) + { + fAllNormsDirFittable = false; + } + RegisterParameters(); } -void -BinnedEDManager::AddPdfs(const std::vector& pdfs_, - const std::vector* norm_fitting_statuses){ - if(norm_fitting_statuses != nullptr && pdfs_.size() != norm_fitting_statuses->size()) { +void BinnedEDManager::AddPdfs(const std::vector &pdfs_, + const std::vector *norm_fitting_statuses) +{ + if (norm_fitting_statuses != nullptr && pdfs_.size() != norm_fitting_statuses->size()) + { throw DimensionError("BinnedEDManager: number of norm_fittable bools doesn't the number of pdfs"); } - for(size_t i = 0; i < pdfs_.size(); i++){ - if (norm_fitting_statuses != nullptr) { + for (size_t i = 0; i < pdfs_.size(); i++) + { + if (norm_fitting_statuses != nullptr) + { AddPdf(pdfs_.at(i), norm_fitting_statuses->at(i)); - } else { + } + else + { AddPdf(pdfs_.at(i)); } } RegisterParameters(); } -const std::vector& -BinnedEDManager::GetNormalisations() const{ +const std::vector & +BinnedEDManager::GetNormalisations() const +{ return fNormalisations; } -void -BinnedEDManager::ApplyShrink(const BinnedEDShrinker& shrinker_){ +void BinnedEDManager::ApplyShrink(const BinnedEDShrinker &shrinker_) +{ if (!shrinker_.GetBuffers().size()) return; - + // only shrink if not already shrunk! FIXME: more obvious behaviour if (!fWorkingPdfs.size()) return; - for (size_t i = 0; i < fWorkingPdfs.size(); i++){ + for (size_t i = 0; i < fWorkingPdfs.size(); i++) + { // Normalise if normalisation is a fittable param, but if indirect then track any change - if (fAllowNormsFittable.at(i) == DIRECT) { - if( fWorkingPdfs.at(0).GetNBins() == fOriginalPdfs.at(0).GetNBins() ){ + if (fAllowNormsFittable.at(i) == DIRECT) + { + if (fWorkingPdfs.at(i).GetNBins() == fOriginalPdfs.at(i).GetNBins()) + { fWorkingPdfs[i] = shrinker_.ShrinkDist(fWorkingPdfs.at(i)); fWorkingPdfs[i].Normalise(); } - else{ + else + { return; } - } else if (fAllowNormsFittable.at(i) == FALSE) { - if ( fWorkingPdfs.at(0).GetNBins() == fOriginalPdfs.at(0).GetNBins() ){ + } + else if (fAllowNormsFittable.at(i) == FALSE) + { + if (fWorkingPdfs.at(i).GetNBins() == fOriginalPdfs.at(i).GetNBins()) + { fWorkingPdfs[i] = shrinker_.ShrinkDist(fWorkingPdfs.at(i)); const double integral_after = fWorkingPdfs[i].Integral(); fNormalisations[i] = integral_after; } - else{ + else + { return; } - } else { - const double integral_before = fWorkingPdfs[i].Integral(); - fWorkingPdfs[i] = shrinker_.ShrinkDist(fWorkingPdfs.at(i)); - const double integral_after = fWorkingPdfs[i].Integral(); - if (integral_before == 0. && integral_after == 0.) { fNormalisations[i] = 0.; } - else { fNormalisations[i] *= integral_after/integral_before; } - } + } + else + { + const double integral_before = fWorkingPdfs[i].Integral(); + fWorkingPdfs[i] = shrinker_.ShrinkDist(fWorkingPdfs.at(i)); + const double integral_after = fWorkingPdfs[i].Integral(); + if (integral_before == 0. && integral_after == 0.) + { + fNormalisations[i] = 0.; + } + else + { + fNormalisations[i] *= integral_after / integral_before; + } + } } - } //////////////////////////////// @@ -176,80 +229,95 @@ BinnedEDManager::ApplyShrink(const BinnedEDShrinker& shrinker_){ //////////////////////////////// std::string -BinnedEDManager::GetName() const{ +BinnedEDManager::GetName() const +{ return fName; } -void -BinnedEDManager::SetName(const std::string& n_){ +void BinnedEDManager::SetName(const std::string &n_) +{ fName = n_; } -void -BinnedEDManager::RenameParameter(const std::string& old_, const std::string& new_){ +void BinnedEDManager::RenameParameter(const std::string &old_, const std::string &new_) +{ fParameterManager.RenameParameter(old_, new_); } -void -BinnedEDManager::SetParameter(const std::string& name_, double value_){ +void BinnedEDManager::SetParameter(const std::string &name_, double value_) +{ fParameterManager.SetParameter(name_, value_); } double -BinnedEDManager::GetParameter(const std::string& name_) const{ +BinnedEDManager::GetParameter(const std::string &name_) const +{ return fParameterManager.GetParameter(name_); } -void -BinnedEDManager::SetParameters(const ParameterDict& ps_){ +void BinnedEDManager::SetParameters(const ParameterDict &ps_) +{ fParameterManager.SetParameters(ps_); } ParameterDict -BinnedEDManager::GetParameters() const{ +BinnedEDManager::GetParameters() const +{ return fParameterManager.GetParameters(); } size_t -BinnedEDManager::GetParameterCount() const{ +BinnedEDManager::GetParameterCount() const +{ return fParameterManager.GetParameterCount(); } std::set -BinnedEDManager::GetParameterNames() const{ +BinnedEDManager::GetParameterNames() const +{ return fParameterManager.GetParameterNames(); } -void -BinnedEDManager::RegisterParameters(){ +void BinnedEDManager::RegisterParameters() +{ fParameterManager.Clear(); std::vector parameterNames; - if (fAllNormsDirFittable) { - for(size_t i = 0; i < fOriginalPdfs.size(); i++) { + if (fAllNormsDirFittable) + { + for (size_t i = 0; i < fOriginalPdfs.size(); i++) + { parameterNames.push_back(fOriginalPdfs.at(i).GetName()); } fParameterManager.AddContainer(fNormalisations, parameterNames); - } else { - for(size_t i = 0; i < fOriginalPdfs.size(); i++) { - if (fAllowNormsFittable.at(i)) { + } + else + { + for (size_t i = 0; i < fOriginalPdfs.size(); i++) + { + if (fAllowNormsFittable.at(i)) + { parameterNames.push_back(fOriginalPdfs.at(i).GetName()); } } - fParameterManager.AddContainer(fFittableNorms, parameterNames); + fParameterManager.AddContainer(fFittableNorms, parameterNames); } -} +} -void BinnedEDManager::ReassertNorms(bool calcing_binprob) { +void BinnedEDManager::ReassertNorms(bool calcing_binprob) +{ /* - * Normalisations that are fittable might have been changed - * in fFittableNorms, but not in fNormalisations! - * This method corrects any changes that may have occured, - * so that both vectors are consistent. - */ - if (!fAllNormsDirFittable) { + * Normalisations that are fittable might have been changed + * in fFittableNorms, but not in fNormalisations! + * This method corrects any changes that may have occured, + * so that both vectors are consistent. + */ + if (!fAllNormsDirFittable) + { size_t j = 0; - for (size_t i = 0; i < fNormalisations.size(); i++) { - if (fAllowNormsFittable.at(i) == ((calcing_binprob) ? DIRECT : INDIRECT)) { + for (size_t i = 0; i < fNormalisations.size(); i++) + { + if (fAllowNormsFittable.at(i) == ((calcing_binprob) ? DIRECT : INDIRECT)) + { fNormalisations[i] = fFittableNorms.at(j); j++; } @@ -257,15 +325,18 @@ void BinnedEDManager::ReassertNorms(bool calcing_binprob) { } } -void BinnedEDManager::AssertDimensions(const std::vector& observables) { +void BinnedEDManager::AssertDimensions(const std::vector &observables) +{ /* * Check whether each PDF in the manager actually has dimensions * matching that of the input vector. * Whenever this is not the case, the given PDF is marginalised onto * the input vector's axes. */ - for (auto &working_pdf : fWorkingPdfs) { - if (working_pdf.GetObservables() != observables) { + for (auto &working_pdf : fWorkingPdfs) + { + if (working_pdf.GetObservables() != observables) + { working_pdf = working_pdf.Marginalise(observables); } } From 7b8e2e417c533a947181aa0f3959760aa7dc121f Mon Sep 17 00:00:00 2001 From: parkerw Date: Tue, 10 Dec 2024 16:06:05 +0000 Subject: [PATCH 5/6] syst group applying reordering --- src/fitutil/SystematicManager.cpp | 51 +++++++++++-------------------- 1 file changed, 17 insertions(+), 34 deletions(-) diff --git a/src/fitutil/SystematicManager.cpp b/src/fitutil/SystematicManager.cpp index c0a3e86..740022d 100644 --- a/src/fitutil/SystematicManager.cpp +++ b/src/fitutil/SystematicManager.cpp @@ -55,10 +55,12 @@ SystematicManager::GetTotalResponse(const std::string &groupName_) const void SystematicManager::Add(Systematic *sys_, const std::string &groupName_) { - if (groupName_ == ""){ - std::cout << "Warning::Adding systematic with empty group name. Will apply to all events." << std::endl; - Add(sys_); - } + if (groupName_ == "") + { + std::cout << "Warning::Adding systematic with empty group name. Will apply to all events." << std::endl; + Add(sys_); + return; + } fGroups[groupName_].push_back(sys_); fNGroups = fGroups.size(); } @@ -151,6 +153,15 @@ void SystematicManager::DistortEDs(const std::vector &OrignalEDs_, if (fEDGroups.find(name) == fEDGroups.end()) continue; + // Copy EDs because WorkingED_s have to have the same binning as OrignalEDs_s. + WorkingEDs_[j] = OrignalEDs_.at(j); + + // If we haven't explicitly added the "" grooup to this dist's fEDGroups, apply that group first + if (std::find(fEDGroups.at(name).begin(), fEDGroups.at(name).end(), "") == fEDGroups.at(name).end() && fGroups.at("").size()) + { + WorkingEDs_[j].SetBinContents(GetTotalResponse("").operator()(WorkingEDs_.at(j).GetBinContents())); + } + // Apply everything else. for (size_t i = 0; i < fEDGroups.at(name).size(); ++i) { @@ -158,41 +169,13 @@ void SystematicManager::DistortEDs(const std::vector &OrignalEDs_, std::string groupName = fEDGroups.at(name).at(i); // If "" is empty the do nothing. if (groupName == "" && fGroups.at("").size() == 0) - { - WorkingEDs_[j] = OrignalEDs_.at(j); continue; - } if (!fGroups.at(groupName).size()) throw LogicError(Formatter() << "SystematicManager:: ED " << name << " has a systematic group of zero size acting on it"); - // Copy EDs because WorkingED_s have to have the same binning as OrignalEDs_s. - WorkingEDs_[j] = OrignalEDs_.at(j); - - // Logic overview: - // - Check whether fGroups has more that just the "" group in it. If - // not then apply the group as you would expect. - // - If groups have been added but there is nothing in "" then apply - // the groups on their own. - // - Else apply the "" group first and then apply the rest. - if (fGroups.size() > 1) - { - if (fGroups.at("").size() == 0) - { - WorkingEDs_[j].SetBinContents(GetTotalResponse(groupName).operator()(OrignalEDs_.at(j).GetBinContents())); - } - else - { - WorkingEDs_[j].SetBinContents( - GetTotalResponse(groupName).operator()( - GetTotalResponse("").operator()(OrignalEDs_.at(j).GetBinContents()))); - } - } - else - { - WorkingEDs_[j].SetBinContents( - GetTotalResponse(groupName).operator()(OrignalEDs_.at(j).GetBinContents())); - } + // Now apply this group to the WorkingEDs + WorkingEDs_[j].SetBinContents(GetTotalResponse(groupName).operator()(WorkingEDs_.at(j).GetBinContents())); } if (norms != nullptr) { From 25ab1453fc86e768c737a5d276d49192cf16b6c8 Mon Sep 17 00:00:00 2001 From: parkerw Date: Wed, 11 Dec 2024 12:07:48 +0000 Subject: [PATCH 6/6] add a bunch of unit tests --- test/unit/BinnedNLLHTest.cpp | 957 ++++++++++++++++++++++++++++++++++- 1 file changed, 955 insertions(+), 2 deletions(-) diff --git a/test/unit/BinnedNLLHTest.cpp b/test/unit/BinnedNLLHTest.cpp index 7fc0912..5c4adb9 100644 --- a/test/unit/BinnedNLLHTest.cpp +++ b/test/unit/BinnedNLLHTest.cpp @@ -4,6 +4,7 @@ #include #include #include +#include #include TEST_CASE("Binned NLLH, 3 rates no systematics") @@ -108,7 +109,7 @@ TEST_CASE("Binned NLLH, 3 rates no systematics") lh_buffer.AddPdf(pdf1); // Add a generous buffer region - lh_buffer.SetBuffer("axis1",30,30); + lh_buffer.SetBuffer("axis1", 30, 30); lh_buffer.SetDataSet(&data); lh_buffer.RegisterFitComponents(); @@ -122,7 +123,7 @@ TEST_CASE("Binned NLLH, 3 rates no systematics") REQUIRE_THAT(llh1, Catch::Matchers::WithinAbs(llh2, 0.0001)); std::pair buffers = lh_buffer.GetBuffer("axis1"); - REQUIRE(buffers == std::pair(30,30)); + REQUIRE(buffers == std::pair(30, 30)); } std::vector genRates(3, pow(10, 6)); @@ -288,3 +289,955 @@ TEST_CASE("Binned NLLH, 3 rates no systematics") REQUIRE_THAT(lh2.Evaluate(), Catch::Matchers::WithinAbs(sumNorm + sumLogProb + betaPen + constraint, 0.0001)); } } + +TEST_CASE("Binned NLLH, Grouping Systematics") +{ + + SECTION("Two PDFs, two systematics in \"\" Group") + { + + // Make BinnedEDs + AxisCollection ax; + // Two bins, each unit width + ax.AddAxis(BinAxis("xaxis", 0, 2, 2)); + BinnedED pdf1("pdf1", ax); + BinnedED pdf2("pdf2", ax); + BinnedED data("data", ax); + + // Some simple bin contents + double bin0_pdf1 = 4.; + double bin1_pdf1 = 6.; + double bin0_pdf2 = 3.; + double bin1_pdf2 = 2.; + double bin0_data = 3.; + double bin1_data = 4.; + pdf1.SetBinContent(0, bin0_pdf1); + pdf1.SetBinContent(1, bin1_pdf1); + pdf2.SetBinContent(0, bin0_pdf2); + pdf2.SetBinContent(1, bin1_pdf2); + data.SetBinContent(0, bin0_data); + data.SetBinContent(1, bin1_data); + + // Now make 2 simple shift systematics + Shift *shift1 = new Shift("shift1"); + shift1->RenameParameter("shift", "xshift1"); + shift1->SetShift(0.0); + ObsSet obs = {"xaxis"}; + shift1->SetAxes(ax); + shift1->SetTransformationObs(obs); + shift1->SetDistributionObs(obs); + shift1->Construct(); + + Shift *shift2 = new Shift("shift2"); + shift2->RenameParameter("shift", "xshift2"); + shift2->SetShift(0.0); + shift2->SetAxes(ax); + shift2->SetTransformationObs(obs); + shift2->SetDistributionObs(obs); + shift2->Construct(); + + // Make a LLH object + BinnedNLLH lh1; + lh1.SetDataDist(data); + lh1.AddSystematic(shift1, ""); + lh1.AddSystematic(shift2, ""); + std::vector group1vec = {""}; + std::vector group2vec = {""}; + lh1.AddPdf(pdf1, group1vec, 100); + lh1.AddPdf(pdf2, group2vec, 100); + lh1.RegisterFitComponents(); + + // And set some parameter values + double shiftval1 = 0.5; + double shiftval2 = 0.1; + double pdfnorm1 = 1.0; + double pdfnorm2 = 1.0; + ParameterDict parameterValues; + parameterValues["pdf1"] = pdfnorm1; + parameterValues["pdf2"] = pdfnorm2; + parameterValues["xshift1"] = shiftval1; + parameterValues["xshift2"] = shiftval2; + + lh1.SetParameters(parameterValues); + double llh = lh1.Evaluate(); + + // Now let's calculate the llh by hand, applying the systematics ourselves + // Bin 0 loses events being shifted into Bin 1: + double bin0_pdf1_shift1 = bin0_pdf1 - shiftval1 * bin0_pdf1; + // Bin 1 loses events being shifted into the overflow, but gains from Bin 0: + double bin1_pdf1_shift1 = bin1_pdf1 + (shiftval1 * bin0_pdf1) - (shiftval1 * bin1_pdf1); + // And we can repeat for the second shift: + double bin0_pdf1_shift2 = bin0_pdf1_shift1 - shiftval2 * bin0_pdf1_shift1; + double bin1_pdf1_shift2 = bin1_pdf1_shift1 + (shiftval2 * bin0_pdf1_shift1) - (shiftval2 * bin1_pdf1_shift1); + + // And we have the same for PDF 2: + double bin0_pdf2_shift1 = bin0_pdf2 - shiftval1 * bin0_pdf2; + double bin1_pdf2_shift1 = bin1_pdf2 + (shiftval1 * bin0_pdf2) - (shiftval1 * bin1_pdf2); + double bin0_pdf2_shift2 = bin0_pdf2_shift1 - shiftval2 * bin0_pdf2_shift1; + double bin1_pdf2_shift2 = bin1_pdf2_shift1 + (shiftval2 * bin0_pdf2_shift1) - (shiftval2 * bin1_pdf2_shift1); + + // Now we can sum the total MC bin contents: + double bin0_totmc = (pdfnorm1 * bin0_pdf1_shift2) + (pdfnorm2 * bin0_pdf2_shift2); + double bin1_totmc = (pdfnorm1 * bin1_pdf1_shift2) + (pdfnorm2 * bin1_pdf2_shift2); + + // And calculate the LLH! Summing -data*log(mc) + mc for each bin + double calcllh = (-bin0_data * log(bin0_totmc) + bin0_totmc) + (-bin1_data * log(bin1_totmc) + bin1_totmc); + + REQUIRE_THAT(llh, Catch::Matchers::WithinAbs(calcllh, 0.0001)); + // And compare to hand calculated value just to be sure: + REQUIRE_THAT(llh, Catch::Matchers::WithinAbs(-1.03259, 0.0001)); + } + + SECTION("Two PDFs, one systematic in \"\" Group, one in separate group") + { + + // Make BinnedEDs + AxisCollection ax; + // Two bins, each unit width + ax.AddAxis(BinAxis("xaxis", 0, 2, 2)); + BinnedED pdf1("pdf1", ax); + BinnedED pdf2("pdf2", ax); + BinnedED data("data", ax); + + // Some simple bin contents + double bin0_pdf1 = 4.; + double bin1_pdf1 = 6.; + double bin0_pdf2 = 3.; + double bin1_pdf2 = 2.; + double bin0_data = 3.; + double bin1_data = 4.; + pdf1.SetBinContent(0, bin0_pdf1); + pdf1.SetBinContent(1, bin1_pdf1); + pdf2.SetBinContent(0, bin0_pdf2); + pdf2.SetBinContent(1, bin1_pdf2); + data.SetBinContent(0, bin0_data); + data.SetBinContent(1, bin1_data); + + // Now make 2 simple shift systematics + Shift *shift1 = new Shift("shift1"); + shift1->RenameParameter("shift", "xshift1"); + shift1->SetShift(0.0); + ObsSet obs = {"xaxis"}; + shift1->SetAxes(ax); + shift1->SetTransformationObs(obs); + shift1->SetDistributionObs(obs); + shift1->Construct(); + + Shift *shift2 = new Shift("shift2"); + shift2->RenameParameter("shift", "xshift2"); + shift2->SetShift(0.0); + shift2->SetAxes(ax); + shift2->SetTransformationObs(obs); + shift2->SetDistributionObs(obs); + shift2->Construct(); + + // Make a LLH object + BinnedNLLH lh1; + lh1.SetDataDist(data); + lh1.AddSystematic(shift1, ""); + lh1.AddSystematic(shift2, "syst2"); + std::vector group1vec = {""}; + std::vector group2vec = {"syst2"}; + lh1.AddPdf(pdf1, group1vec, 100); + lh1.AddPdf(pdf2, group2vec, 100); + lh1.RegisterFitComponents(); + + // And set some parameter values + double shiftval1 = 0.5; + double shiftval2 = 0.1; + double pdfnorm1 = 1.0; + double pdfnorm2 = 1.0; + ParameterDict parameterValues; + parameterValues["pdf1"] = pdfnorm1; + parameterValues["pdf2"] = pdfnorm2; + parameterValues["xshift1"] = shiftval1; + parameterValues["xshift2"] = shiftval2; + + lh1.SetParameters(parameterValues); + double llh = lh1.Evaluate(); + + // Now let's calculate the llh by hand, applying the systematics ourselves + // Bin 0 loses events being shifted into Bin 1: + double bin0_pdf1_shift1 = bin0_pdf1 - shiftval1 * bin0_pdf1; + // Bin 1 loses events being shifted into the overflow, but gains from Bin 0: + double bin1_pdf1_shift1 = bin1_pdf1 + (shiftval1 * bin0_pdf1) - (shiftval1 * bin1_pdf1); + // But the second shift doesn''t apply to PDF 1 + + // And we have the same for PDF 2: + double bin0_pdf2_shift1 = bin0_pdf2 - shiftval1 * bin0_pdf2; + double bin1_pdf2_shift1 = bin1_pdf2 + (shiftval1 * bin0_pdf2) - (shiftval1 * bin1_pdf2); + double bin0_pdf2_shift2 = bin0_pdf2_shift1 - shiftval2 * bin0_pdf2_shift1; + double bin1_pdf2_shift2 = bin1_pdf2_shift1 + (shiftval2 * bin0_pdf2_shift1) - (shiftval2 * bin1_pdf2_shift1); + + // Now we can sum the total MC bin contents: + double bin0_totmc = (pdfnorm1 * bin0_pdf1_shift1) + (pdfnorm2 * bin0_pdf2_shift2); + double bin1_totmc = (pdfnorm1 * bin1_pdf1_shift1) + (pdfnorm2 * bin1_pdf2_shift2); + + // And calculate the LLH! Summing -data*log(mc) + mc for each bin + double calcllh = (-bin0_data * log(bin0_totmc) + bin0_totmc) + (-bin1_data * log(bin1_totmc) + bin1_totmc); + + REQUIRE_THAT(llh, Catch::Matchers::WithinAbs(calcllh, 0.0001)); + REQUIRE_THAT(llh, Catch::Matchers::WithinAbs(-0.88280, 0.0001)); + } + + SECTION("Two PDFs, two systematics in separate groups") + { + + // Make BinnedEDs + AxisCollection ax; + // Two bins, each unit width + ax.AddAxis(BinAxis("xaxis", 0, 2, 2)); + BinnedED pdf1("pdf1", ax); + BinnedED pdf2("pdf2", ax); + BinnedED data("data", ax); + + // Some simple bin contents + double bin0_pdf1 = 4.; + double bin1_pdf1 = 6.; + double bin0_pdf2 = 3.; + double bin1_pdf2 = 2.; + double bin0_data = 3.; + double bin1_data = 4.; + pdf1.SetBinContent(0, bin0_pdf1); + pdf1.SetBinContent(1, bin1_pdf1); + pdf2.SetBinContent(0, bin0_pdf2); + pdf2.SetBinContent(1, bin1_pdf2); + data.SetBinContent(0, bin0_data); + data.SetBinContent(1, bin1_data); + + // Now make 2 simple shift systematics + Shift *shift1 = new Shift("shift1"); + shift1->RenameParameter("shift", "xshift1"); + shift1->SetShift(0.0); + ObsSet obs = {"xaxis"}; + shift1->SetAxes(ax); + shift1->SetTransformationObs(obs); + shift1->SetDistributionObs(obs); + shift1->Construct(); + + Shift *shift2 = new Shift("shift2"); + shift2->RenameParameter("shift", "xshift2"); + shift2->SetShift(0.0); + shift2->SetAxes(ax); + shift2->SetTransformationObs(obs); + shift2->SetDistributionObs(obs); + shift2->Construct(); + + // Make a LLH object + BinnedNLLH lh1; + lh1.SetDataDist(data); + lh1.AddSystematic(shift1, "syst1"); + lh1.AddSystematic(shift2, "syst2"); + std::vector group1vec = {"syst1"}; + std::vector group2vec = {"syst2"}; + lh1.AddPdf(pdf1, group1vec, 100); + lh1.AddPdf(pdf2, group2vec, 100); + lh1.RegisterFitComponents(); + + // And set some parameter values + double shiftval1 = 0.5; + double shiftval2 = 0.1; + double pdfnorm1 = 1.0; + double pdfnorm2 = 1.0; + ParameterDict parameterValues; + parameterValues["pdf1"] = pdfnorm1; + parameterValues["pdf2"] = pdfnorm2; + parameterValues["xshift1"] = shiftval1; + parameterValues["xshift2"] = shiftval2; + + lh1.SetParameters(parameterValues); + double llh = lh1.Evaluate(); + + // Now let's calculate the llh by hand, applying the systematics ourselves + // Bin 0 loses events being shifted into Bin 1: + double bin0_pdf1_shift1 = bin0_pdf1 - shiftval1 * bin0_pdf1; + // Bin 1 loses events being shifted into the overflow, but gains from Bin 0: + double bin1_pdf1_shift1 = bin1_pdf1 + (shiftval1 * bin0_pdf1) - (shiftval1 * bin1_pdf1); + // But the second shift doesn''t apply to PDF 1 + + // And we have the same for PDF 2 and shift2: + double bin0_pdf2_shift2 = bin0_pdf2 - shiftval2 * bin0_pdf2; + double bin1_pdf2_shift2 = bin1_pdf2 + (shiftval2 * bin0_pdf2) - (shiftval2 * bin1_pdf2); + + // Now we can sum the total MC bin contents: + double bin0_totmc = (pdfnorm1 * bin0_pdf1_shift1) + (pdfnorm2 * bin0_pdf2_shift2); + double bin1_totmc = (pdfnorm1 * bin1_pdf1_shift1) + (pdfnorm2 * bin1_pdf2_shift2); + + // And calculate the LLH! Summing -data*log(mc) + mc for each bin + double calcllh = (-bin0_data * log(bin0_totmc) + bin0_totmc) + (-bin1_data * log(bin1_totmc) + bin1_totmc); + + REQUIRE_THAT(llh, Catch::Matchers::WithinAbs(calcllh, 0.0001)); + REQUIRE_THAT(llh, Catch::Matchers::WithinAbs(-0.68307, 0.0001)); + } + + SECTION("Two PDFs, one systematic in \"\" Group, one in separate group which contains both PDFs") + { + + // Make BinnedEDs + AxisCollection ax; + // Two bins, each unit width + ax.AddAxis(BinAxis("xaxis", 0, 2, 2)); + BinnedED pdf1("pdf1", ax); + BinnedED pdf2("pdf2", ax); + BinnedED data("data", ax); + + // Some simple bin contents + double bin0_pdf1 = 4.; + double bin1_pdf1 = 6.; + double bin0_pdf2 = 3.; + double bin1_pdf2 = 2.; + double bin0_data = 3.; + double bin1_data = 4.; + pdf1.SetBinContent(0, bin0_pdf1); + pdf1.SetBinContent(1, bin1_pdf1); + pdf2.SetBinContent(0, bin0_pdf2); + pdf2.SetBinContent(1, bin1_pdf2); + data.SetBinContent(0, bin0_data); + data.SetBinContent(1, bin1_data); + + // Now make 2 simple shift systematics + Shift *shift1 = new Shift("shift1"); + shift1->RenameParameter("shift", "xshift1"); + shift1->SetShift(0.0); + ObsSet obs = {"xaxis"}; + shift1->SetAxes(ax); + shift1->SetTransformationObs(obs); + shift1->SetDistributionObs(obs); + shift1->Construct(); + + Shift *shift2 = new Shift("shift2"); + shift2->RenameParameter("shift", "xshift2"); + shift2->SetShift(0.0); + shift2->SetAxes(ax); + shift2->SetTransformationObs(obs); + shift2->SetDistributionObs(obs); + shift2->Construct(); + + // Make a LLH object + BinnedNLLH lh1; + lh1.SetDataDist(data); + lh1.AddSystematic(shift1, ""); + lh1.AddSystematic(shift2, "syst2"); + std::vector group1vec = {"syst2"}; + std::vector group2vec = {"syst2"}; + lh1.AddPdf(pdf1, group1vec, 100); + lh1.AddPdf(pdf2, group2vec, 100); + lh1.RegisterFitComponents(); + + // And set some parameter values + double shiftval1 = 0.5; + double shiftval2 = 0.1; + double pdfnorm1 = 1.0; + double pdfnorm2 = 1.0; + ParameterDict parameterValues; + parameterValues["pdf1"] = pdfnorm1; + parameterValues["pdf2"] = pdfnorm2; + parameterValues["xshift1"] = shiftval1; + parameterValues["xshift2"] = shiftval2; + + lh1.SetParameters(parameterValues); + double llh = lh1.Evaluate(); + + // Now let's calculate the llh by hand, applying the systematics ourselves + // Bin 0 loses events being shifted into Bin 1: + double bin0_pdf1_shift1 = bin0_pdf1 - shiftval1 * bin0_pdf1; + // Bin 1 loses events being shifted into the overflow, but gains from Bin 0: + double bin1_pdf1_shift1 = bin1_pdf1 + (shiftval1 * bin0_pdf1) - (shiftval1 * bin1_pdf1); + // And we can repeat for the second shift: + double bin0_pdf1_shift2 = bin0_pdf1_shift1 - shiftval2 * bin0_pdf1_shift1; + double bin1_pdf1_shift2 = bin1_pdf1_shift1 + (shiftval2 * bin0_pdf1_shift1) - (shiftval2 * bin1_pdf1_shift1); + + // And we have the same for PDF 2 and shift2: + double bin0_pdf2_shift1 = bin0_pdf2 - shiftval1 * bin0_pdf2; + double bin1_pdf2_shift1 = bin1_pdf2 + (shiftval1 * bin0_pdf2) - (shiftval1 * bin1_pdf2); + double bin0_pdf2_shift2 = bin0_pdf2_shift1 - shiftval2 * bin0_pdf2_shift1; + double bin1_pdf2_shift2 = bin1_pdf2_shift1 + (shiftval2 * bin0_pdf2_shift1) - (shiftval2 * bin1_pdf2_shift1); + + // Now we can sum the total MC bin contents: + double bin0_totmc = (pdfnorm1 * bin0_pdf1_shift2) + (pdfnorm2 * bin0_pdf2_shift2); + double bin1_totmc = (pdfnorm1 * bin1_pdf1_shift2) + (pdfnorm2 * bin1_pdf2_shift2); + + // And calculate the LLH! Summing -data*log(mc) + mc for each bin + double calcllh = (-bin0_data * log(bin0_totmc) + bin0_totmc) + (-bin1_data * log(bin1_totmc) + bin1_totmc); + + REQUIRE_THAT(llh, Catch::Matchers::WithinAbs(calcllh, 0.0001)); + REQUIRE_THAT(llh, Catch::Matchers::WithinAbs(-1.03259, 0.0001)); + } + + SECTION("Three PDFs, two systematics with separate groups") + { + + // Make BinnedEDs + AxisCollection ax; + // Two bins, each unit width + ax.AddAxis(BinAxis("xaxis", 0, 2, 2)); + BinnedED pdf1("pdf1", ax); + BinnedED pdf2("pdf2", ax); + BinnedED pdf3("pdf3", ax); + BinnedED data("data", ax); + + // Some simple bin contents + double bin0_pdf1 = 4.; + double bin1_pdf1 = 6.; + double bin0_pdf2 = 3.; + double bin1_pdf2 = 2.; + double bin0_pdf3 = 2.; + double bin1_pdf3 = 4.; + double bin0_data = 3.; + double bin1_data = 4.; + pdf1.SetBinContent(0, bin0_pdf1); + pdf1.SetBinContent(1, bin1_pdf1); + pdf2.SetBinContent(0, bin0_pdf2); + pdf2.SetBinContent(1, bin1_pdf2); + pdf3.SetBinContent(0, bin0_pdf3); + pdf3.SetBinContent(1, bin1_pdf3); + data.SetBinContent(0, bin0_data); + data.SetBinContent(1, bin1_data); + + // Now make 2 simple shift systematics + Shift *shift1 = new Shift("shift1"); + shift1->RenameParameter("shift", "xshift1"); + shift1->SetShift(0.0); + ObsSet obs = {"xaxis"}; + shift1->SetAxes(ax); + shift1->SetTransformationObs(obs); + shift1->SetDistributionObs(obs); + shift1->Construct(); + + Shift *shift2 = new Shift("shift2"); + shift2->RenameParameter("shift", "xshift2"); + shift2->SetShift(0.0); + shift2->SetAxes(ax); + shift2->SetTransformationObs(obs); + shift2->SetDistributionObs(obs); + shift2->Construct(); + + // Make a LLH object + BinnedNLLH lh1; + lh1.SetDataDist(data); + lh1.AddSystematic(shift1, "syst1"); + lh1.AddSystematic(shift2, "syst2"); + std::vector group1vec = {"syst1"}; + std::vector group2vec = {"syst2"}; + std::vector group3vec = {"syst2"}; + lh1.AddPdf(pdf1, group1vec, 100); + lh1.AddPdf(pdf2, group2vec, 100); + lh1.AddPdf(pdf3, group3vec, 100); + lh1.RegisterFitComponents(); + + // And set some parameter values + double shiftval1 = 0.5; + double shiftval2 = 0.1; + double pdfnorm1 = 1.0; + double pdfnorm2 = 1.0; + double pdfnorm3 = 1.0; + ParameterDict parameterValues; + parameterValues["pdf1"] = pdfnorm1; + parameterValues["pdf2"] = pdfnorm2; + parameterValues["pdf3"] = pdfnorm3; + parameterValues["xshift1"] = shiftval1; + parameterValues["xshift2"] = shiftval2; + + lh1.SetParameters(parameterValues); + double llh = lh1.Evaluate(); + + // Now let's calculate the llh by hand, applying the systematics ourselves + // Bin 0 loses events being shifted into Bin 1: + double bin0_pdf1_shift1 = bin0_pdf1 - shiftval1 * bin0_pdf1; + // Bin 1 loses events being shifted into the overflow, but gains from Bin 0: + double bin1_pdf1_shift1 = bin1_pdf1 + (shiftval1 * bin0_pdf1) - (shiftval1 * bin1_pdf1); + // But the second shift doesn't apply to PDF 1 + + // And we have shift2 for PDFs 2 and 3: + double bin0_pdf2_shift2 = bin0_pdf2 - shiftval2 * bin0_pdf2; + double bin1_pdf2_shift2 = bin1_pdf2 + (shiftval2 * bin0_pdf2) - (shiftval2 * bin1_pdf2); + double bin0_pdf3_shift2 = bin0_pdf3 - shiftval2 * bin0_pdf3; + double bin1_pdf3_shift2 = bin1_pdf3 + (shiftval2 * bin0_pdf3) - (shiftval2 * bin1_pdf3); + + // Now we can sum the total MC bin contents: + double bin0_totmc = (pdfnorm1 * bin0_pdf1_shift1) + (pdfnorm2 * bin0_pdf2_shift2) + (pdfnorm3 * bin0_pdf3_shift2); + double bin1_totmc = (pdfnorm1 * bin1_pdf1_shift1) + (pdfnorm2 * bin1_pdf2_shift2) + (pdfnorm3 * bin1_pdf3_shift2); + + // And calculate the LLH! Summing -data*log(mc) + mc for each bin + double calcllh = (-bin0_data * log(bin0_totmc) + bin0_totmc) + (-bin1_data * log(bin1_totmc) + bin1_totmc); + + REQUIRE_THAT(llh, Catch::Matchers::WithinAbs(calcllh, 0.0001)); + REQUIRE_THAT(llh, Catch::Matchers::WithinAbs(2.22954, 0.0001)); + } + + SECTION("Three PDFs, two systematics in different groups, one in \"\" group") + { + + // Make BinnedEDs + AxisCollection ax; + // Two bins, each unit width + ax.AddAxis(BinAxis("xaxis", 0, 2, 2)); + BinnedED pdf1("pdf1", ax); + BinnedED pdf2("pdf2", ax); + BinnedED pdf3("pdf3", ax); + BinnedED data("data", ax); + + // Some simple bin contents + double bin0_pdf1 = 4.; + double bin1_pdf1 = 6.; + double bin0_pdf2 = 3.; + double bin1_pdf2 = 2.; + double bin0_pdf3 = 2.; + double bin1_pdf3 = 4.; + double bin0_data = 3.; + double bin1_data = 4.; + pdf1.SetBinContent(0, bin0_pdf1); + pdf1.SetBinContent(1, bin1_pdf1); + pdf2.SetBinContent(0, bin0_pdf2); + pdf2.SetBinContent(1, bin1_pdf2); + pdf3.SetBinContent(0, bin0_pdf3); + pdf3.SetBinContent(1, bin1_pdf3); + data.SetBinContent(0, bin0_data); + data.SetBinContent(1, bin1_data); + + // Now make 2 simple shift systematics + Shift *shift1 = new Shift("shift1"); + shift1->RenameParameter("shift", "xshift1"); + shift1->SetShift(0.0); + ObsSet obs = {"xaxis"}; + shift1->SetAxes(ax); + shift1->SetTransformationObs(obs); + shift1->SetDistributionObs(obs); + shift1->Construct(); + + Shift *shift2 = new Shift("shift2"); + shift2->RenameParameter("shift", "xshift2"); + shift2->SetShift(0.0); + shift2->SetAxes(ax); + shift2->SetTransformationObs(obs); + shift2->SetDistributionObs(obs); + shift2->Construct(); + + Shift *shift3 = new Shift("shift3"); + shift3->RenameParameter("shift", "xshift3"); + shift3->SetShift(0.0); + shift3->SetAxes(ax); + shift3->SetTransformationObs(obs); + shift3->SetDistributionObs(obs); + shift3->Construct(); + + // Make a LLH object + BinnedNLLH lh1; + lh1.SetDataDist(data); + lh1.AddSystematic(shift1, ""); + lh1.AddSystematic(shift2, "syst2"); + lh1.AddSystematic(shift3, "syst3"); + std::vector group1vec = {""}; + std::vector group2vec = {"syst2"}; + std::vector group3vec = {"syst3"}; + lh1.AddPdf(pdf1, group1vec, 100); + lh1.AddPdf(pdf2, group2vec, 100); + lh1.AddPdf(pdf3, group3vec, 100); + lh1.RegisterFitComponents(); + + // And set some parameter values + double shiftval1 = 0.5; + double shiftval2 = 0.1; + double shiftval3 = -0.3; + double pdfnorm1 = 1.0; + double pdfnorm2 = 1.0; + double pdfnorm3 = 1.0; + ParameterDict parameterValues; + parameterValues["pdf1"] = pdfnorm1; + parameterValues["pdf2"] = pdfnorm2; + parameterValues["pdf3"] = pdfnorm3; + parameterValues["xshift1"] = shiftval1; + parameterValues["xshift2"] = shiftval2; + parameterValues["xshift3"] = shiftval3; + + lh1.SetParameters(parameterValues); + double llh = lh1.Evaluate(); + + // Now let's calculate the llh by hand, applying the systematics ourselves + // Bin 0 loses events being shifted into Bin 1: + double bin0_pdf1_shift1 = bin0_pdf1 - shiftval1 * bin0_pdf1; + // Bin 1 loses events being shifted into the overflow, but gains from Bin 0: + double bin1_pdf1_shift1 = bin1_pdf1 + (shiftval1 * bin0_pdf1) - (shiftval1 * bin1_pdf1); + // But the second and third shifts don't apply to PDF 1 + + // And we have shift2 for PDF 2: + double bin0_pdf2_shift1 = bin0_pdf2 - shiftval1 * bin0_pdf2; + double bin1_pdf2_shift1 = bin1_pdf2 + (shiftval1 * bin0_pdf2) - (shiftval1 * bin1_pdf2); + double bin0_pdf2_shift2 = bin0_pdf2_shift1 - shiftval2 * bin0_pdf2_shift1; + double bin1_pdf2_shift2 = bin1_pdf2_shift1 + (shiftval2 * bin0_pdf2_shift1) - (shiftval2 * bin1_pdf2_shift1); + + // And we have shift3 for PDFs 3: + double bin0_pdf3_shift1 = bin0_pdf3 - shiftval1 * bin0_pdf3; + double bin1_pdf3_shift1 = bin1_pdf3 + (shiftval1 * bin0_pdf3) - (shiftval1 * bin1_pdf3); + // As shift 3 has a negative value, the shift goes the other way + double bin0_pdf3_shift3 = bin0_pdf3_shift1 - (shiftval3 * bin1_pdf3_shift1) + (shiftval3 * bin0_pdf3_shift1); + double bin1_pdf3_shift3 = bin1_pdf3_shift1 + shiftval3 * bin1_pdf3_shift1; + + // Now we can sum the total MC bin contents: + double bin0_totmc = (pdfnorm1 * bin0_pdf1_shift1) + (pdfnorm2 * bin0_pdf2_shift2) + (pdfnorm3 * bin0_pdf3_shift3); + double bin1_totmc = (pdfnorm1 * bin1_pdf1_shift1) + (pdfnorm2 * bin1_pdf2_shift2) + (pdfnorm3 * bin1_pdf3_shift3); + + // And calculate the LLH! Summing -data*log(mc) + mc for each bin + double calcllh = (-bin0_data * log(bin0_totmc) + bin0_totmc) + (-bin1_data * log(bin1_totmc) + bin1_totmc); + + REQUIRE_THAT(llh, Catch::Matchers::WithinAbs(calcllh, 0.0001)); + REQUIRE_THAT(llh, Catch::Matchers::WithinAbs(0.64667, 0.0001)); + } + + SECTION("Three PDFs, two systematics in \"\" group, one in a separate group") + { + + // Make BinnedEDs + AxisCollection ax; + // Two bins, each unit width + ax.AddAxis(BinAxis("xaxis", 0, 2, 2)); + BinnedED pdf1("pdf1", ax); + BinnedED pdf2("pdf2", ax); + BinnedED pdf3("pdf3", ax); + BinnedED data("data", ax); + + // Some simple bin contents + double bin0_pdf1 = 4.; + double bin1_pdf1 = 6.; + double bin0_pdf2 = 3.; + double bin1_pdf2 = 2.; + double bin0_pdf3 = 2.; + double bin1_pdf3 = 4.; + double bin0_data = 3.; + double bin1_data = 4.; + pdf1.SetBinContent(0, bin0_pdf1); + pdf1.SetBinContent(1, bin1_pdf1); + pdf2.SetBinContent(0, bin0_pdf2); + pdf2.SetBinContent(1, bin1_pdf2); + pdf3.SetBinContent(0, bin0_pdf3); + pdf3.SetBinContent(1, bin1_pdf3); + data.SetBinContent(0, bin0_data); + data.SetBinContent(1, bin1_data); + + // Now make 2 simple shift systematics + Shift *shift1 = new Shift("shift1"); + shift1->RenameParameter("shift", "xshift1"); + shift1->SetShift(0.0); + ObsSet obs = {"xaxis"}; + shift1->SetAxes(ax); + shift1->SetTransformationObs(obs); + shift1->SetDistributionObs(obs); + shift1->Construct(); + + Shift *shift2 = new Shift("shift2"); + shift2->RenameParameter("shift", "xshift2"); + shift2->SetShift(0.0); + shift2->SetAxes(ax); + shift2->SetTransformationObs(obs); + shift2->SetDistributionObs(obs); + shift2->Construct(); + + Shift *shift3 = new Shift("shift3"); + shift3->RenameParameter("shift", "xshift3"); + shift3->SetShift(0.0); + shift3->SetAxes(ax); + shift3->SetTransformationObs(obs); + shift3->SetDistributionObs(obs); + shift3->Construct(); + + // Make a LLH object + BinnedNLLH lh1; + lh1.SetDataDist(data); + lh1.AddSystematic(shift1, ""); + lh1.AddSystematic(shift2, ""); + lh1.AddSystematic(shift3, "syst3"); + std::vector group1vec = {""}; + std::vector group2vec = {""}; + std::vector group3vec = {"syst3"}; + lh1.AddPdf(pdf1, group1vec, 100); + lh1.AddPdf(pdf2, group2vec, 100); + lh1.AddPdf(pdf3, group3vec, 100); + lh1.RegisterFitComponents(); + + // And set some parameter values + double shiftval1 = 0.5; + double shiftval2 = 0.1; + double shiftval3 = -0.3; + double pdfnorm1 = 1.0; + double pdfnorm2 = 1.0; + double pdfnorm3 = 1.0; + ParameterDict parameterValues; + parameterValues["pdf1"] = pdfnorm1; + parameterValues["pdf2"] = pdfnorm2; + parameterValues["pdf3"] = pdfnorm3; + parameterValues["xshift1"] = shiftval1; + parameterValues["xshift2"] = shiftval2; + parameterValues["xshift3"] = shiftval3; + + lh1.SetParameters(parameterValues); + double llh = lh1.Evaluate(); + + // Now let's calculate the llh by hand, applying the systematics ourselves + // Bin 0 loses events being shifted into Bin 1: + double bin0_pdf1_shift1 = bin0_pdf1 - shiftval1 * bin0_pdf1; + // Bin 1 loses events being shifted into the overflow, but gains from Bin 0: + double bin1_pdf1_shift1 = bin1_pdf1 + (shiftval1 * bin0_pdf1) - (shiftval1 * bin1_pdf1); + // And do the same for shift 2: + double bin0_pdf1_shift2 = bin0_pdf1_shift1 - shiftval2 * bin0_pdf1_shift1; + double bin1_pdf1_shift2 = bin1_pdf1_shift1 + (shiftval2 * bin0_pdf1_shift1) - (shiftval2 * bin1_pdf1_shift1); + + // And similarly for PDF 2: + double bin0_pdf2_shift1 = bin0_pdf2 - shiftval1 * bin0_pdf2; + double bin1_pdf2_shift1 = bin1_pdf2 + (shiftval1 * bin0_pdf2) - (shiftval1 * bin1_pdf2); + double bin0_pdf2_shift2 = bin0_pdf2_shift1 - shiftval2 * bin0_pdf2_shift1; + double bin1_pdf2_shift2 = bin1_pdf2_shift1 + (shiftval2 * bin0_pdf2_shift1) - (shiftval2 * bin1_pdf2_shift1); + + // And we have shift3 for PDFs 3: + double bin0_pdf3_shift1 = bin0_pdf3 - shiftval1 * bin0_pdf3; + double bin1_pdf3_shift1 = bin1_pdf3 + (shiftval1 * bin0_pdf3) - (shiftval1 * bin1_pdf3); + double bin0_pdf3_shift2 = bin0_pdf3_shift1 - shiftval2 * bin0_pdf3_shift1; + double bin1_pdf3_shift2 = bin1_pdf3_shift1 + (shiftval2 * bin0_pdf3_shift1) - (shiftval2 * bin1_pdf3_shift1); + // As shift 3 has a negative value, the shift goes the other way + double bin0_pdf3_shift3 = bin0_pdf3_shift2 - (shiftval3 * bin1_pdf3_shift2) + (shiftval3 * bin0_pdf3_shift2); + double bin1_pdf3_shift3 = bin1_pdf3_shift2 + shiftval3 * bin1_pdf3_shift2; + + // Now we can sum the total MC bin contents: + double bin0_totmc = (pdfnorm1 * bin0_pdf1_shift2) + (pdfnorm2 * bin0_pdf2_shift2) + (pdfnorm3 * bin0_pdf3_shift3); + double bin1_totmc = (pdfnorm1 * bin1_pdf1_shift2) + (pdfnorm2 * bin1_pdf2_shift2) + (pdfnorm3 * bin1_pdf3_shift3); + + // And calculate the LLH! Summing -data*log(mc) + mc for each bin + double calcllh = (-bin0_data * log(bin0_totmc) + bin0_totmc) + (-bin1_data * log(bin1_totmc) + bin1_totmc); + + REQUIRE_THAT(llh, Catch::Matchers::WithinAbs(calcllh, 0.0001)); + REQUIRE_THAT(llh, Catch::Matchers::WithinAbs(0.27339, 0.0001)); + } + + SECTION("Three PDFs, three systematics all in separate groups") + { + + // Make BinnedEDs + AxisCollection ax; + // Two bins, each unit width + ax.AddAxis(BinAxis("xaxis", 0, 2, 2)); + BinnedED pdf1("pdf1", ax); + BinnedED pdf2("pdf2", ax); + BinnedED pdf3("pdf3", ax); + BinnedED data("data", ax); + + // Some simple bin contents + double bin0_pdf1 = 4.; + double bin1_pdf1 = 6.; + double bin0_pdf2 = 3.; + double bin1_pdf2 = 2.; + double bin0_pdf3 = 2.; + double bin1_pdf3 = 4.; + double bin0_data = 3.; + double bin1_data = 4.; + pdf1.SetBinContent(0, bin0_pdf1); + pdf1.SetBinContent(1, bin1_pdf1); + pdf2.SetBinContent(0, bin0_pdf2); + pdf2.SetBinContent(1, bin1_pdf2); + pdf3.SetBinContent(0, bin0_pdf3); + pdf3.SetBinContent(1, bin1_pdf3); + data.SetBinContent(0, bin0_data); + data.SetBinContent(1, bin1_data); + + // Now make 2 simple shift systematics + Shift *shift1 = new Shift("shift1"); + shift1->RenameParameter("shift", "xshift1"); + shift1->SetShift(0.0); + ObsSet obs = {"xaxis"}; + shift1->SetAxes(ax); + shift1->SetTransformationObs(obs); + shift1->SetDistributionObs(obs); + shift1->Construct(); + + Shift *shift2 = new Shift("shift2"); + shift2->RenameParameter("shift", "xshift2"); + shift2->SetShift(0.0); + shift2->SetAxes(ax); + shift2->SetTransformationObs(obs); + shift2->SetDistributionObs(obs); + shift2->Construct(); + + Shift *shift3 = new Shift("shift3"); + shift3->RenameParameter("shift", "xshift3"); + shift3->SetShift(0.0); + shift3->SetAxes(ax); + shift3->SetTransformationObs(obs); + shift3->SetDistributionObs(obs); + shift3->Construct(); + + // Make a LLH object + BinnedNLLH lh1; + lh1.SetDataDist(data); + lh1.AddSystematic(shift1, "syst1"); + lh1.AddSystematic(shift2, "syst2"); + lh1.AddSystematic(shift3, "syst3"); + std::vector group1vec = {"syst1"}; + std::vector group2vec = {"syst2"}; + std::vector group3vec = {"syst3"}; + lh1.AddPdf(pdf1, group1vec, 100); + lh1.AddPdf(pdf2, group2vec, 100); + lh1.AddPdf(pdf3, group3vec, 100); + lh1.RegisterFitComponents(); + + // And set some parameter values + double shiftval1 = 0.5; + double shiftval2 = 0.1; + double shiftval3 = -0.3; + double pdfnorm1 = 1.0; + double pdfnorm2 = 1.0; + double pdfnorm3 = 1.0; + ParameterDict parameterValues; + parameterValues["pdf1"] = pdfnorm1; + parameterValues["pdf2"] = pdfnorm2; + parameterValues["pdf3"] = pdfnorm3; + parameterValues["xshift1"] = shiftval1; + parameterValues["xshift2"] = shiftval2; + parameterValues["xshift3"] = shiftval3; + + lh1.SetParameters(parameterValues); + double llh = lh1.Evaluate(); + + // Now let's calculate the llh by hand, applying the systematics ourselves + // Bin 0 loses events being shifted into Bin 1: + double bin0_pdf1_shift1 = bin0_pdf1 - shiftval1 * bin0_pdf1; + // Bin 1 loses events being shifted into the overflow, but gains from Bin 0: + double bin1_pdf1_shift1 = bin1_pdf1 + (shiftval1 * bin0_pdf1) - (shiftval1 * bin1_pdf1); + + // And similarly for PDF 2: + double bin0_pdf2_shift2 = bin0_pdf2 - shiftval2 * bin0_pdf2; + double bin1_pdf2_shift2 = bin1_pdf2 + (shiftval2 * bin0_pdf2) - (shiftval2 * bin1_pdf2); + + // And we have shift3 for PDFs 3: + double bin0_pdf3_shift3 = bin0_pdf3 - (shiftval3 * bin1_pdf3) + (shiftval3 * bin0_pdf3); + double bin1_pdf3_shift3 = bin1_pdf3 + shiftval3 * bin1_pdf3; + + // Now we can sum the total MC bin contents: + double bin0_totmc = (pdfnorm1 * bin0_pdf1_shift1) + (pdfnorm2 * bin0_pdf2_shift2) + (pdfnorm3 * bin0_pdf3_shift3); + double bin1_totmc = (pdfnorm1 * bin1_pdf1_shift1) + (pdfnorm2 * bin1_pdf2_shift2) + (pdfnorm3 * bin1_pdf3_shift3); + + // And calculate the LLH! Summing -data*log(mc) + mc for each bin + double calcllh = (-bin0_data * log(bin0_totmc) + bin0_totmc) + (-bin1_data * log(bin1_totmc) + bin1_totmc); + + REQUIRE_THAT(llh, Catch::Matchers::WithinAbs(calcllh, 0.0001)); + REQUIRE_THAT(llh, Catch::Matchers::WithinAbs(2.06624, 0.0001)); + } + + SECTION("Three PDFs, one systematic in the \"\" group, one applies to one PDF, one applies to two PDFs") + { + + // Make BinnedEDs + AxisCollection ax; + // Two bins, each unit width + ax.AddAxis(BinAxis("xaxis", 0, 2, 2)); + BinnedED pdf1("pdf1", ax); + BinnedED pdf2("pdf2", ax); + BinnedED pdf3("pdf3", ax); + BinnedED data("data", ax); + + // Some simple bin contents + double bin0_pdf1 = 4.; + double bin1_pdf1 = 6.; + double bin0_pdf2 = 3.; + double bin1_pdf2 = 2.; + double bin0_pdf3 = 2.; + double bin1_pdf3 = 4.; + double bin0_data = 3.; + double bin1_data = 4.; + pdf1.SetBinContent(0, bin0_pdf1); + pdf1.SetBinContent(1, bin1_pdf1); + pdf2.SetBinContent(0, bin0_pdf2); + pdf2.SetBinContent(1, bin1_pdf2); + pdf3.SetBinContent(0, bin0_pdf3); + pdf3.SetBinContent(1, bin1_pdf3); + data.SetBinContent(0, bin0_data); + data.SetBinContent(1, bin1_data); + + // Now make 2 simple shift systematics + Shift *shift1 = new Shift("shift1"); + shift1->RenameParameter("shift", "xshift1"); + shift1->SetShift(0.0); + ObsSet obs = {"xaxis"}; + shift1->SetAxes(ax); + shift1->SetTransformationObs(obs); + shift1->SetDistributionObs(obs); + shift1->Construct(); + + Shift *shift2 = new Shift("shift2"); + shift2->RenameParameter("shift", "xshift2"); + shift2->SetShift(0.0); + shift2->SetAxes(ax); + shift2->SetTransformationObs(obs); + shift2->SetDistributionObs(obs); + shift2->Construct(); + + Shift *shift3 = new Shift("shift3"); + shift3->RenameParameter("shift", "xshift3"); + shift3->SetShift(0.0); + shift3->SetAxes(ax); + shift3->SetTransformationObs(obs); + shift3->SetDistributionObs(obs); + shift3->Construct(); + + // Make a LLH object + BinnedNLLH lh1; + lh1.SetDataDist(data); + lh1.AddSystematic(shift1, ""); + lh1.AddSystematic(shift2, "syst2"); + lh1.AddSystematic(shift3, "syst3"); + std::vector group1vec = {""}; + std::vector group2vec = {"syst2", "syst3"}; + std::vector group3vec = {"syst3"}; + lh1.AddPdf(pdf1, group1vec, 100); + lh1.AddPdf(pdf2, group2vec, 100); + lh1.AddPdf(pdf3, group3vec, 100); + lh1.RegisterFitComponents(); + + // And set some parameter values + double shiftval1 = 0.5; + double shiftval2 = 0.1; + double shiftval3 = -0.3; + double pdfnorm1 = 1.0; + double pdfnorm2 = 1.0; + double pdfnorm3 = 1.0; + ParameterDict parameterValues; + parameterValues["pdf1"] = pdfnorm1; + parameterValues["pdf2"] = pdfnorm2; + parameterValues["pdf3"] = pdfnorm3; + parameterValues["xshift1"] = shiftval1; + parameterValues["xshift2"] = shiftval2; + parameterValues["xshift3"] = shiftval3; + + lh1.SetParameters(parameterValues); + double llh = lh1.Evaluate(); + + // Now let's calculate the llh by hand, applying the systematics ourselves + // Bin 0 loses events being shifted into Bin 1: + double bin0_pdf1_shift1 = bin0_pdf1 - shiftval1 * bin0_pdf1; + // Bin 1 loses events being shifted into the overflow, but gains from Bin 0: + double bin1_pdf1_shift1 = bin1_pdf1 + (shiftval1 * bin0_pdf1) - (shiftval1 * bin1_pdf1); + + // And shifts 1, 2 and 3 for PDF 2: + double bin0_pdf2_shift1 = bin0_pdf2 - shiftval1 * bin0_pdf2; + double bin1_pdf2_shift1 = bin1_pdf2 + (shiftval1 * bin0_pdf2) - (shiftval1 * bin1_pdf2); + double bin0_pdf2_shift2 = bin0_pdf2_shift1 - shiftval2 * bin0_pdf2_shift1; + double bin1_pdf2_shift2 = bin1_pdf2_shift1 + (shiftval2 * bin0_pdf2_shift1) - (shiftval2 * bin1_pdf2_shift1); + // As shift 3 has a negative value, the shift goes the other way + double bin0_pdf2_shift3 = bin0_pdf2_shift2 - (shiftval3 * bin1_pdf2_shift2) + (shiftval3 * bin0_pdf2_shift2); + double bin1_pdf2_shift3 = bin1_pdf2_shift2 + shiftval3 * bin1_pdf2_shift2; + + // And we have shift 1 and 3 for PDFs 3: + double bin0_pdf3_shift1 = bin0_pdf3 - shiftval1 * bin0_pdf3; + double bin1_pdf3_shift1 = bin1_pdf3 + (shiftval1 * bin0_pdf3) - (shiftval1 * bin1_pdf3); + // As shift 3 has a negative value, the shift goes the other way + double bin0_pdf3_shift3 = bin0_pdf3_shift1 - (shiftval3 * bin1_pdf3_shift1) + (shiftval3 * bin0_pdf3_shift1); + double bin1_pdf3_shift3 = bin1_pdf3_shift1 + shiftval3 * bin1_pdf3_shift1; + + // Now we can sum the total MC bin contents: + double bin0_totmc = (pdfnorm1 * bin0_pdf1_shift1) + (pdfnorm2 * bin0_pdf2_shift3) + (pdfnorm3 * bin0_pdf3_shift3); + double bin1_totmc = (pdfnorm1 * bin1_pdf1_shift1) + (pdfnorm2 * bin1_pdf2_shift3) + (pdfnorm3 * bin1_pdf3_shift3); + + // And calculate the LLH! Summing -data*log(mc) + mc for each bin + double calcllh = (-bin0_data * log(bin0_totmc) + bin0_totmc) + (-bin1_data * log(bin1_totmc) + bin1_totmc); + + REQUIRE_THAT(llh, Catch::Matchers::WithinAbs(calcllh, 0.0001)); + REQUIRE_THAT(llh, Catch::Matchers::WithinAbs(0.371851, 0.0001)); + } +}