From e1d121e76301428e5bccc2754a774b08e55f157d Mon Sep 17 00:00:00 2001 From: Gregory J Ottino Date: Tue, 26 May 2026 14:23:01 -0400 Subject: [PATCH] added Systematic ratio divider macro --- .../SystematicDetermination.C | 236 ++++++++++++++++++ 1 file changed, 236 insertions(+) create mode 100644 LightFlavorRatios/systematic_uncertainty/SystematicDetermination.C diff --git a/LightFlavorRatios/systematic_uncertainty/SystematicDetermination.C b/LightFlavorRatios/systematic_uncertainty/SystematicDetermination.C new file mode 100644 index 0000000000..10b08a485d --- /dev/null +++ b/LightFlavorRatios/systematic_uncertainty/SystematicDetermination.C @@ -0,0 +1,236 @@ +#include +#include + +#include "TCanvas.h" +#include "TFile.h" +#include "TH1F.h" +#include "TH2F.h" +#include "TStyle.h" +#include +#include +#include "TMatrixD.h" + + +namespace { + +TMatrixD BuildSystematicCorrelationMatrix(const std::vector& baseline, + const std::vector& syst, + bool systIsRelative = true, + bool fullyCorrelated = true) +{ + const int n = static_cast(baseline.size()); + if (n == 0 || static_cast(syst.size()) != n) { + throw std::runtime_error("Input vectors must have same non-zero size."); + } + + std::vector delta(n, 0.0); // absolute systematic shift per bin + for (int i = 0; i < n; ++i) { + delta[i] = systIsRelative ? baseline[i] * syst[i] : syst[i]; + } + + TMatrixD cov(n, n); + cov.Zero(); + + if (fullyCorrelated) { + // One systematic source shared across bins: Cov_ij = delta_i *delta_j + for (int i = 0; i < n; ++i) + for (int j = 0; j < n; ++j) + cov(i, j) = delta[i] * delta[j]; + } else { + // Independent bin-by-bin systematic nuisances + for (int i = 0; i < n; ++i) + cov(i, i) = delta[i] * delta[i]; + } + + TMatrixD corr(n, n); + corr.Zero(); + + for (int i = 0; i < n; ++i) { + for (int j = 0; j < n; ++j) { + const double denom = std::sqrt(cov(i, i) * cov(j, j)); + corr(i, j) = (denom > 0.0) ? cov(i, j) / denom : 0.0; + } + } + + return corr; + } + + +bool HaveMatchingBinning(const TH1F* a, const TH1F* b) +{ + if (!a || !b) { + return false; + } + + if (a->GetNbinsX() != b->GetNbinsX()) { + return false; + } + + const int nBins = a->GetNbinsX(); + for (int i = 1; i <= nBins + 1; ++i) { + if (std::fabs(a->GetXaxis()->GetBinLowEdge(i) - b->GetXaxis()->GetBinLowEdge(i)) > 1e-12) { + return false; + } + } + + return true; +} + +void FillRatio(const TH1F* numerator, const TH1F* denominator, TH1F* ratio) +{ + const int nBins = ratio->GetNbinsX(); + for (int i = 1; i <= nBins; ++i) { + const double n = numerator->GetBinContent(i); + const double dn = numerator->GetBinError(i); + const double d = denominator->GetBinContent(i); + const double dd = denominator->GetBinError(i); + + if (d == 0.0 || n == 0.0) { + ratio->SetBinContent(i, 0.0); + ratio->SetBinError(i, 0.0); + continue; + } + + const double r = n / d; + const double relErr2 = (dn / n) * (dn / n) + (dd / d) * (dd / d); + ratio->SetBinContent(i, r); + ratio->SetBinError(i, std::fabs(r) * std::sqrt(relErr2)); + } +} + +} // namespace + +void SystematicDetermination(const char* inputFileName = "/sphenix/user/rosstom/analysis/LightFlavorRatios/systematics/geoAccEffSystematics.root", + const char* baselineParticle1HistName = "h_pT_base_KS0", + const char* baselineParticle2HistName = "h_pT_base_Xi", + const char* variationParticle1HistName = "h_pT_base_KS0_max", + const char* variationParticle2HistName = "h_pT_base_Xi_max", + const char* outputFileName = "SystematicDetermination_output.root") +{ + TFile* inputFile = TFile::Open(inputFileName, "READ"); + if (!inputFile || inputFile->IsZombie()) { + std::cerr << "Error: could not open input file: " << inputFileName << std::endl; + return; + } + + TH1F* hBaseP1 = dynamic_cast(inputFile->Get(baselineParticle1HistName)); + TH1F* hBaseP2 = dynamic_cast(inputFile->Get(baselineParticle2HistName)); + TH1F* hVarP1 = dynamic_cast(inputFile->Get(variationParticle1HistName)); + TH1F* hVarP2 = dynamic_cast(inputFile->Get(variationParticle2HistName)); + + if (!hBaseP1 || !hBaseP2 || !hVarP1 || !hVarP2) { + std::cerr << "Error: one or more requested TH1F histograms were not found." << std::endl; + inputFile->Close(); + return; + } + + if (!HaveMatchingBinning(hBaseP1, hBaseP2) || + !HaveMatchingBinning(hBaseP1, hVarP1) || + !HaveMatchingBinning(hBaseP1, hVarP2)) { + std::cerr << "Error: input histograms do not have matching binning." << std::endl; + inputFile->Close(); + return; + } + + TH1F* hRatioBaseline = dynamic_cast(hBaseP2->Clone("hRatioBaseline")); + TH1F* hRatioVariation = dynamic_cast(hVarP2->Clone("hRatioVariation")); + TH1F* hDifference = dynamic_cast(hBaseP2->Clone("hDifference")); + + hRatioBaseline->Reset("ICES"); + hRatioVariation->Reset("ICES"); + hDifference->Reset("ICES"); + + hRatioBaseline->SetTitle("Baseline ratio (particle2 / particle1)"); + hRatioVariation->SetTitle("Variation ratio (particle2 / particle1)"); + hDifference->SetTitle("Relative Difference: (variation - baseline) / baseline"); + + FillRatio(hBaseP2, hBaseP1, hRatioBaseline); + FillRatio(hVarP2, hVarP1, hRatioVariation); + + const int nBins = hRatioBaseline->GetNbinsX(); + for (int i = 1; i <= nBins; ++i) { + const double b = hRatioBaseline->GetBinContent(i); + const double vb = hRatioBaseline->GetBinError(i); + const double v = hRatioVariation->GetBinContent(i); + const double vv = hRatioVariation->GetBinError(i); + + if (b == 0.0) { + hDifference->SetBinContent(i, 0.0); + hDifference->SetBinError(i, 0.0); + continue; + } + + const double diff = (v - b)/b; + //const double diffErr = std::sqrt((vv / b) * (vv / b) + (v * vb / (b * b)) * (v * vb / (b * b))); + + hDifference->SetBinContent(i, diff); + //hDifference->SetBinError(i, diffErr); + } + + TH2F* hCorrelationMatrix = new TH2F("hSystematicCorrelation", "Bin-to-bin correlation of systematic uncertainty", nBins, 0.5, nBins + 0.5, nBins, 0.5, nBins + 0.5); + hCorrelationMatrix->GetXaxis()->SetTitle("Bin i"); + hCorrelationMatrix->GetYaxis()->SetTitle("Bin j"); + + for (int i = 1; i <= nBins; ++i) { + const double di = hDifference->GetBinContent(i); + for (int j = 1; j <= nBins; ++j) { + const double dj = hDifference->GetBinContent(j); + + double corr = 0.0; + if (i == j) { + corr = 1.0; + } + + if (di != 0.0 && dj != 0.0) { + corr = (di * dj) / (std::fabs(di) * std::fabs(dj)); + } + + hCorrelationMatrix->SetBinContent(i, j, corr); + } + } + + gStyle->SetOptStat(0); + + TCanvas* cRatioBaseline = new TCanvas("cRatioBaseline", "Baseline Ratio", 800, 600); + hRatioBaseline->SetLineColor(kBlue + 1); + hRatioBaseline->SetMarkerColor(kBlue + 1); + hRatioBaseline->SetMarkerStyle(20); + hRatioBaseline->Draw("E1"); + + TCanvas* cRatioVariation = new TCanvas("cRatioVariation", "Variation Ratio", 800, 600); + hRatioVariation->SetLineColor(kRed + 1); + hRatioVariation->SetMarkerColor(kRed + 1); + hRatioVariation->SetMarkerStyle(20); + hRatioVariation->Draw("E1"); + + TCanvas* cRelativeDifference = new TCanvas("cRelativeDifference", "Relative Difference", 800, 600); + hDifference->SetLineColor(kBlack); + hDifference->SetMarkerColor(kBlack); + hDifference->SetMarkerStyle(21); + hDifference->Draw("E1"); + + TCanvas* cCorrelation = new TCanvas("cCorrelation", "Systematic Correlation Matrix", 800, 700); + hCorrelationMatrix->Draw("COLZ TEXT"); + + TFile* outputFile = TFile::Open(outputFileName, "RECREATE"); + if (!outputFile || outputFile->IsZombie()) { + std::cerr << "Error: could not create output file: " << outputFileName << std::endl; + inputFile->Close(); + return; + } + + hRatioBaseline->Write(); + hRatioVariation->Write(); + hDifference->Write(); + hCorrelationMatrix->Write(); + + cRatioBaseline->Write(); + cRatioVariation->Write(); + cRelativeDifference->Write(); + cCorrelation->Write(); + + outputFile->Close(); + inputFile->Close(); + + std::cout << "Wrote output to: " << outputFileName << std::endl; +}