Skip to content

Add file reading/writing for flattened dchi2 dist #442

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
nathanielerowe wants to merge 8 commits into
base: project-SBN-dev
Choose a base branch
from
Open

Add file reading/writing for flattened dchi2 dist #442

Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
bf00d8b
Add file reading/writing for flattened dchi2 dist
nathanielerowe Jun 23, 2026
9b7e1d0
Fix dtype
nathanielerowe Jun 23, 2026
8a16452
add command line arg
nathanielerowe Jun 23, 2026
ea75a5c
fix missing var
nathanielerowe Jun 23, 2026
98da064
fitx bug with multithreading
nathanielerowe Jun 25, 2026
79c98d2
Combine with existing FC files
nathanielerowe Jul 15, 2026
1b8b59e
Add catches for reusing files
nathanielerowe Aug 6, 2026
dad86ba
cleaning
nathanielerowe Aug 7, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
344 changes: 264 additions & 80 deletions bin/PROfit.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,154 @@

using namespace PROfit;

std::vector<fc_out> readAnalysisFromFile(const std::string& filename,
const PROmodel &model, // To parse log10 parameters
const std::vector<PROsyst>& variable_systs, // To map spline names to Eigen indices
size_t i_prime, // Current system config index
bool gof_mode_requested) { // Requested runtime mode
std::vector<fc_out> entries;

// 1. Open the file
TFile* inFile = TFile::Open(filename.c_str(), "READ");
if (!inFile || inFile->IsZombie()) {
log<LOG_ERROR>(L"%1% || Error: Could not open file %2%") % __func__ % filename.c_str();
return entries;
}

// 2. Get the TTree
TTree* tree = nullptr;
inFile->GetObject("tree", tree);
if (!tree) {
log<LOG_ERROR>(L"%1% || Error: Could not find TTree 'tree' in %2%") % __func__ % filename.c_str();
inFile->Close();
delete inFile;
return entries;
}

// 3. Inspect branch existence
TBranch* b_osc = tree->GetBranch("chi2_osc");
TBranch* b_sosc = tree->GetBranch("best_systs_osc");

if (!b_osc || !b_sosc) {
if (!gof_mode_requested) {
log<LOG_ERROR>(L"%1% || ERROR: The reused file '%2%' is missing oscillation branches ('chi2_osc' / 'best_systs_osc').")
% __func__ % filename.c_str();
log<LOG_ERROR>(L"%1% || Cannot compute standard Feldman-Cousins p-values (--pval) without oscillation fit data.") % __func__;
log<LOG_ERROR>(L"%1% || Please re-run fresh throws without --gof to generate a complete distribution.") % __func__;
inFile->Close();
delete inFile;
return entries; // Returns empty vector to halt main execution
}
}

// 4. Set up local reading buffers
float chi2_osc = 0;
float chi2_syst = 0;
std::map<std::string, float>* p_best_systs_osc = nullptr;
std::map<std::string, float>* p_best_systs = nullptr;
std::map<std::string, float>* p_syst_throw = nullptr;

// Allocate memory for physical parameter values matching model size
std::vector<float> phys_param_vals(model.nparams, 0.0f);

// 5. Bind tree branches to our buffers
if (b_osc) tree->SetBranchAddress("chi2_osc", &chi2_osc);
tree->SetBranchAddress("chi2_syst", &chi2_syst);
if (b_sosc) tree->SetBranchAddress("best_systs_osc", &p_best_systs_osc);
tree->SetBranchAddress("best_systs", &p_best_systs);
tree->SetBranchAddress("syst_throw", &p_syst_throw);

// Bind physical parameter branches using the model's exact parameter names
for (size_t i = 0; i < model.nparams; ++i) {
std::string branch_name = "best_" + model.param_names[i];
if (tree->GetBranch(branch_name.c_str())) {
tree->SetBranchAddress(branch_name.c_str(), &phys_param_vals[i]);
} else {
log<LOG_WARNING>(L"%1% || Warning: Branch %2% not found in file.") % __func__ % branch_name.c_str();
phys_param_vals[i] = 0.0f;
}
}

// Get the systematic size details
size_t n_splines = variable_systs[i_prime].GetNSplines();

// 6. Loop over all records, reconstruct fc_out structures, and validate oscillation values
Long64_t nEntries = tree->GetEntries();
entries.reserve(nEntries);

bool non_zero_osc_found = false;

for (Long64_t entry = 0; entry < nEntries; ++entry) {
// Reset buffers to sentinel value before reading entry
chi2_osc = -999.0f;
chi2_syst = -999.0f;

tree->GetEntry(entry);

// Only count oscillation fits as valid if the branch existed AND yielded a physical chi2 (chi2 >= 0)
if (b_osc && chi2_osc >= 0.0f) {
non_zero_osc_found = true;
}

fc_out fco;
fco.chi2_osc = chi2_osc;
fco.chi2_syst = chi2_syst;

// --- Reconstruct best_phys_osc (and apply inverse conversion if log-scaled) ---
fco.best_phys_osc.resize(model.nparams);
for (size_t i = 0; i < model.nparams; ++i) {
float val = phys_param_vals[i];
if (model.is_log10[i]) {
fco.best_phys_osc(i) = (val > 0.0f) ? std::log10(val) : 0.0f;
} else {
fco.best_phys_osc(i) = val;
}
}

// --- Reconstruct Eigen Vectors from maps ---
fco.best_fit_osc.setZero(n_splines);
fco.best_fit_syst.setZero(n_splines);
fco.syst_throw.setZero(n_splines);

for (size_t i = 0; i < n_splines; ++i) {
const std::string& spline_name = variable_systs[i_prime].spline_names[i];

if (p_best_systs_osc && p_best_systs_osc->count(spline_name)) {
fco.best_fit_osc(i) = (*p_best_systs_osc)[spline_name];
}
if (p_best_systs && p_best_systs->count(spline_name)) {
fco.best_fit_syst(i) = (*p_best_systs)[spline_name];
}
if (p_syst_throw && p_syst_throw->count(spline_name)) {
fco.syst_throw(i) = (*p_syst_throw)[spline_name];
}
}

entries.push_back(fco);
}

// Clean up ROOT IO
inFile->Close();
delete inFile;

// 7. Validate oscillation data contents against requested mode
if (!gof_mode_requested && !non_zero_osc_found) {
log<LOG_ERROR>(L"%1% || ERROR: The reused file '%2%' contains only zeroed oscillation fits (chi2_osc == 0).")
% __func__ % filename.c_str();
log<LOG_ERROR>(L"%1% || This occurs when reusing a file generated with --gof.") % __func__;
log<LOG_ERROR>(L"%1% || Cannot compute standard Feldman-Cousins p-values (--pval) from a GOF-only distribution.") % __func__;
log<LOG_ERROR>(L"%1% || Please re-run fresh throws without --gof to generate a full distribution.") % __func__;
entries.clear();
return entries;
}

if (gof_mode_requested && non_zero_osc_found) {
log<LOG_INFO>(L"%1% || Reusing full Feldman-Cousins distribution for GOF evaluation (extracting chi2_syst).") % __func__;
}

return entries;
}

// Unique key for DetVar propeller maps (names can be reused across sections).
static std::string DetVarKey(const PROconfig& config, size_t file_index) {
const auto& dv = config.m_detvar_files[file_index];
Expand Down Expand Up @@ -305,6 +453,7 @@ int main(int argc, char* argv[])
PlotBounds pbounds;
size_t nuniv;
bool gof_pvalue = false;
bool reuse_dist = false;
bool pvalue = false;


Expand Down Expand Up @@ -408,6 +557,7 @@ int main(int argc, char* argv[])
CLI::App *profc_command = app.add_subcommand("fc", "Run Feldman-Cousins for this injected signal");
profc_command->add_option("-u,--universes", nuniv, "Number of Feldman Cousins universes to throw")->default_val(1000);
profc_command->add_flag("--gof", gof_pvalue, "Get GOF pvalue");
profc_command->add_flag("--reuse", reuse_dist, "Reuse existing chi2_null.root file for pvalue calculation.");
profc_command->add_flag("--pval", pvalue, "Get FC pvalue")->excludes("--gof");

//PROglobal
Expand Down Expand Up @@ -2483,34 +2633,67 @@ int main(int argc, char* argv[])
fc_PB_configs.push_back({int(nuniv/FCthreads), "Thread " + std::to_string(i)});
}
MultiPROgressBar fc_progress(fc_PB_configs);
fc_progress.initialize_display();
fc_progress.start_display_thread();

bool gen_null_dist = true;
std::string FC_file = analysis_tag+"_FC.root";
if(reuse_dist){
TFile* file = TFile::Open(FC_file.c_str(), "READ");
if (!file) {
gen_null_dist = true;
}
else {
gen_null_dist = false;
}
}
if(gen_null_dist){
fc_progress.initialize_display();
fc_progress.start_display_thread();
}

for(size_t i = 0; i < FCthreads; i++) {
dchi2s.emplace_back();
outs.emplace_back();
fc_args args{todo + (i >= addone), &dchi2s.back(), &outs.back(), config, prop, variable_systs[config.i_prime], chi2, fakeDataParams, L, scanFitConfig,(*myseed.getThreadSeeds())[i], (int)i, !eventbyevent, gof_mode};
std::vector<float> flattened_dchi2s;
if(gen_null_dist){
for(size_t i = 0; i < FCthreads; i++) {
dchi2s.emplace_back();
outs.emplace_back();
fc_args args{todo + (i >= addone), &dchi2s.back(), &outs.back(), config, prop, variable_systs[config.i_prime], chi2, fakeDataParams, L, scanFitConfig,(*myseed.getThreadSeeds())[i], (int)i, !eventbyevent, gof_mode};


threads.emplace_back([args, &fc_progress]() {
PROfit::fc_worker(args, std::ref(fc_progress));
});
}
for(auto&& t: threads) {
t.join();
threads.emplace_back([args, &fc_progress]() {
PROfit::fc_worker(args, std::ref(fc_progress));
});
}
for(auto&& t: threads) {
t.join();
}
fc_progress.finish_all();
}
fc_progress.finish_all();
else if(!gen_null_dist){
fc_progress.finish_all();
std::vector<fc_out> flat_outs = readAnalysisFromFile(FC_file, *model, variable_systs, config.i_prime, gof_pvalue);
nuniv = flat_outs.size();
if(flat_outs.empty()) {
log<LOG_ERROR>(L"%1% || Aborting execution due to invalid or unreadable Feldman-Cousins distribution.") % __func__;
return 1;
}

std::vector<float> flattened_dchi2s;
for(const auto& v: dchi2s) for(const auto& dchi2: v) flattened_dchi2s.push_back(dchi2);
outs.push_back(flat_outs);
}
// Unify array flattening across fresh generation and reused distributions
std::vector<float> flattened_syst_chi2;
for(const auto &out: outs) {
for(const auto &fco: out) {
float dchi2 = gof_pvalue ? fco.chi2_syst : (fco.chi2_syst - fco.chi2_osc);
flattened_dchi2s.push_back(dchi2);
flattened_syst_chi2.push_back(fco.chi2_syst);
}
}
std::sort(flattened_dchi2s.begin(), flattened_dchi2s.end());
log<LOG_INFO>(L"%1% || 90%% Feldman-Cousins delta chi2 after throwing %2% universes is %3%")
std::sort(flattened_syst_chi2.begin(), flattened_syst_chi2.end());

log<LOG_INFO>(L"%1% || 90%% Feldman-Cousins delta chi2 after throwing %2% universes is %3%")
% __func__ % nuniv % flattened_dchi2s[0.9*flattened_dchi2s.size()];

if(gof_pvalue) {
std::vector<float> flattened_syst_chi2;
for(const auto &out : outs) for(const auto &fco : out) flattened_syst_chi2.push_back(fco.chi2_syst);
std::sort(flattened_syst_chi2.begin(), flattened_syst_chi2.end());

log<LOG_ERROR>(L"%1% || All: %2% ") % __func__ % flattened_syst_chi2;
log<LOG_ERROR>(L"%1% || chi: %2% ") % __func__ % global_chi2;
auto it = std::lower_bound(flattened_syst_chi2.begin(), flattened_syst_chi2.end(), global_chi2);
Expand All @@ -2532,68 +2715,69 @@ int main(int argc, char* argv[])
log<LOG_ERROR>(L"%1% || FC Corrected pval after throwing %2% universes is %3%") % __func__ % nuniv % pvalFC ;
}

{
TFile fout((final_output_tag+"_FC.root").c_str(), "RECREATE");
fout.cd();
float chi2_osc, chi2_syst;
// One float per physics parameter — plain branches named "best_<param_name>".
// Vector kept alive for the full lifetime of the TTree.
std::vector<float> best_phys_vals(model->nparams, 0.0f);
std::map<std::string, float> best_systs_osc, best_systs, syst_throw;
TTree tree("tree", "tree");
tree.Branch("chi2_osc", &chi2_osc);
tree.Branch("chi2_syst", &chi2_syst);
for(size_t i = 0; i < model->nparams; ++i)
tree.Branch(("best_" + model->param_names[i]).c_str(), &best_phys_vals[i]);
tree.Branch("best_systs_osc", &best_systs_osc);
tree.Branch("best_systs", &best_systs);
tree.Branch("syst_throw", &syst_throw);

for(const auto &out: outs) {
for(const auto &fco: out) {
chi2_osc = fco.chi2_osc;
chi2_syst = fco.chi2_syst;
for(size_t i = 0; i < model->nparams; ++i) {
float raw = fco.best_phys_osc.size() > (Eigen::Index)i ? fco.best_phys_osc(i) : 0.0f;
best_phys_vals[i] = model->is_log10[i] ? std::pow(10.0f, raw) : raw;
}
for(size_t i = 0; i < variable_systs[config.i_prime].GetNSplines(); ++i) {
if(!gof_pvalue) best_systs_osc[variable_systs[config.i_prime].spline_names[i]] = fco.best_fit_osc(i);
best_systs[variable_systs[config.i_prime].spline_names[i]] = fco.best_fit_syst(i);
syst_throw[variable_systs[config.i_prime].spline_names[i]] = fco.syst_throw(i);
}
tree.Fill();
}
}

tree.Write();
}
{
ofstream fcout(final_output_tag+"_FC.csv");
fcout << "chi2_osc,chi2_syst";
for(const std::string &name: model->param_names)
fcout << ",best_" << name;
for(const std::string &name: variable_systs[config.i_prime].spline_names)
fcout << ",best_" << name << "_osc,best_" << name << "," << name << "_throw";
fcout << "\r\n";

for(const auto &out: outs) {
for(const auto &fco: out) {
fcout << fco.chi2_osc << "," << fco.chi2_syst;
for(size_t i = 0; i < model->nparams; ++i) {
float raw = fco.best_phys_osc.size() > (Eigen::Index)i ? fco.best_phys_osc(i) : 0.0f;
float val = model->is_log10[i] ? std::pow(10.0f, raw) : raw;
fcout << "," << val;
}
for(size_t i = 0; i < variable_systs[config.i_prime].GetNSplines(); ++i)
fcout << "," << (gof_pvalue ? 0 : fco.best_fit_osc(i)) << "," << fco.best_fit_syst(i) << "," << fco.syst_throw(i);
fcout << "\r\n";
}
}
}
//if(save_null_dist){
{
TFile fout((analysis_tag+"_FC.root").c_str(), "RECREATE");
fout.cd();
float chi2_osc, chi2_syst;
// One float per physics parameter — plain branches named "best_<param_name>".
// Vector kept alive for the full lifetime of the TTree.
std::vector<float> best_phys_vals(model->nparams, 0.0f);
std::map<std::string, float> best_systs_osc, best_systs, syst_throw;
TTree tree("tree", "tree");
tree.Branch("chi2_osc", &chi2_osc);
tree.Branch("chi2_syst", &chi2_syst);
for(size_t i = 0; i < model->nparams; ++i)
tree.Branch(("best_" + model->param_names[i]).c_str(), &best_phys_vals[i]);
tree.Branch("best_systs_osc", &best_systs_osc);
tree.Branch("best_systs", &best_systs);
tree.Branch("syst_throw", &syst_throw);

for(const auto &out: outs) {
for(const auto &fco: out) {
chi2_osc = fco.chi2_osc;
chi2_syst = fco.chi2_syst;
for(size_t i = 0; i < model->nparams; ++i) {
float raw = fco.best_phys_osc.size() > (Eigen::Index)i ? fco.best_phys_osc(i) : 0.0f;
best_phys_vals[i] = model->is_log10[i] ? std::pow(10.0f, raw) : raw;
}
for(size_t i = 0; i < variable_systs[config.i_prime].GetNSplines(); ++i) {
if(!gof_pvalue) best_systs_osc[variable_systs[config.i_prime].spline_names[i]] = fco.best_fit_osc(i);
best_systs[variable_systs[config.i_prime].spline_names[i]] = fco.best_fit_syst(i);
syst_throw[variable_systs[config.i_prime].spline_names[i]] = fco.syst_throw(i);
}
tree.Fill();
}
}

tree.Write();
}
{
ofstream fcout(analysis_tag+"_FC.csv");
fcout << "chi2_osc,chi2_syst";
for(const std::string &name: model->param_names)
fcout << ",best_" << name;
for(const std::string &name: variable_systs[config.i_prime].spline_names)
fcout << ",best_" << name << "_osc,best_" << name << "," << name << "_throw";
fcout << "\r\n";

for(const auto &out: outs) {
for(const auto &fco: out) {
fcout << fco.chi2_osc << "," << fco.chi2_syst;
for(size_t i = 0; i < model->nparams; ++i) {
float raw = fco.best_phys_osc.size() > (Eigen::Index)i ? fco.best_phys_osc(i) : 0.0f;
float val = model->is_log10[i] ? std::pow(10.0f, raw) : raw;
fcout << "," << val;
}
for(size_t i = 0; i < variable_systs[config.i_prime].GetNSplines(); ++i)
fcout << "," << (gof_pvalue ? 0 : fco.best_fit_osc(i)) << "," << fco.best_fit_syst(i) << "," << fco.syst_throw(i);
fcout << "\r\n";
}
}
}
//}
}


//***********************************************************************
//***********************************************************************
//******************** global **************************
Expand Down
Loading