From bf00d8bdcc0b69b11a943b74d0faea42d9db5159 Mon Sep 17 00:00:00 2001 From: Nathaniel Date: Tue, 23 Jun 2026 12:15:22 -0500 Subject: [PATCH 1/8] Add file reading/writing for flattened dchi2 dist --- bin/PROfit.cxx | 125 ++++++++++++++++++++++++++++++++++++++++++------- 1 file changed, 107 insertions(+), 18 deletions(-) diff --git a/bin/PROfit.cxx b/bin/PROfit.cxx index f49d1c5e..75a9d9bc 100644 --- a/bin/PROfit.cxx +++ b/bin/PROfit.cxx @@ -49,6 +49,77 @@ using namespace PROfit; + +void writeVectorToFile(const std::string& filename, const std::vector& dataVector) { + // 1. Create/Overwrite the specified ROOT file + TFile* outFile = TFile::Open(filename.c_str(), "RECREATE"); + if (!outFile || outFile->IsZombie()) { + std::cerr << "Error: Could not create file " << filename << std::endl; + return; + } + + // 2. Create the TTree + TTree* tree = new TTree("myTree", "A tree with a vector"); + + // 3. ROOT needs a pointer to the vector. + // Since dataVector is const&, we make a local pointer to a non-const copy + // or simply point to a local copy to ensure safety. + std::vector localCopy = dataVector; + std::vector* vecPtr = &localCopy; + + // 4. Create the branch pointing to our pointer address + tree->Branch("floatVector", &vecPtr); + + // 5. Fill the tree with the vector data (creates 1 entry containing the entire vector) + tree->Fill(); + + // 6. Write and clean up + outFile->Write(); + outFile->Close(); + delete outFile; + + std::cout << "Successfully wrote vector of size " << dataVector.size() << " to " << filename << std::endl; +} + +std::vector readVectorFromFile(const std::string& filename) { + std::vector resultVector; + + // 1. Open the specified file + TFile* inFile = TFile::Open(filename.c_str(), "READ"); + if (!inFile || inFile->IsZombie()) { + std::cerr << "Error: Could not open file " << filename << std::endl; + return resultVector; // Returns empty vector + } + + // 2. Get the TTree + TTree* tree = nullptr; + inFile->GetObject("myTree", tree); + if (!tree) { + std::cerr << "Error: Could not find TTree 'myTree' in " << filename << std::endl; + inFile->Close(); + delete inFile; + return resultVector; + } + + // 3. Set up the pointer for ROOT tracking + std::vector* vecPtr = nullptr; + tree->SetBranchAddress("floatVector", &vecPtr); + + // 4. Read the first entry (assuming 1 entry containing your saved vector) + if (tree->GetEntries() > 0) { + tree->GetEntry(0); + if (vecPtr) { + resultVector = *vecPtr; // Copy the data to our return vector + } + } + + // 5. Clean up ROOT memory + inFile->Close(); + delete inFile; + + return resultVector; +} + // 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]; @@ -2486,27 +2557,45 @@ int main(int argc, char* argv[]) 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}; - - - threads.emplace_back([args, &fc_progress]() { - PROfit::fc_worker(args, std::ref(fc_progress)); - }); - } - for(auto&& t: threads) { - t.join(); - } - fc_progress.finish_all(); + bool save_null_dist = true; + if(pvalue){ + std::string pvalue_file = "chi2_null.root"; + TFile* file = TFile::Open(filename.c_str(), "READ"); + if (!file) { + save_null_dist = true; + } + else { + save_null_dist = false; + } + } std::vector flattened_dchi2s; - for(const auto& v: dchi2s) for(const auto& dchi2: v) flattened_dchi2s.push_back(dchi2); - std::sort(flattened_dchi2s.begin(), flattened_dchi2s.end()); - log(L"%1% || 90%% Feldman-Cousins delta chi2 after throwing %2% universes is %3%") + if(save_null_dist || !pvalue){ + 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(); + } + fc_progress.finish_all(); + + for(const auto& v: dchi2s) for(const auto& dchi2: v) flattened_dchi2s.push_back(dchi2); + std::sort(flattened_dchi2s.begin(), flattened_dchi2s.end()); + writeVectorToFile(pvalue_file, flattened_dchi2s); + } + else if(!save_null_dist){ + flattened_dchi2s = readVectorFromFile(pvalue_file); + } + log(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 flattened_syst_chi2; for(const auto &out : outs) for(const auto &fco : out) flattened_syst_chi2.push_back(fco.chi2_syst); From 9b7e1d0516c61175708b1ba0a65c5974314b3004 Mon Sep 17 00:00:00 2001 From: Nathaniel Date: Tue, 23 Jun 2026 13:00:49 -0500 Subject: [PATCH 2/8] Fix dtype --- bin/PROfit.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/bin/PROfit.cxx b/bin/PROfit.cxx index 75a9d9bc..7fa87331 100644 --- a/bin/PROfit.cxx +++ b/bin/PROfit.cxx @@ -2588,10 +2588,10 @@ int main(int argc, char* argv[]) for(const auto& v: dchi2s) for(const auto& dchi2: v) flattened_dchi2s.push_back(dchi2); std::sort(flattened_dchi2s.begin(), flattened_dchi2s.end()); - writeVectorToFile(pvalue_file, flattened_dchi2s); + writeVectorToFile(&pvalue_file, flattened_dchi2s); } else if(!save_null_dist){ - flattened_dchi2s = readVectorFromFile(pvalue_file); + flattened_dchi2s = readVectorFromFile(&pvalue_file); } log(L"%1% || 90%% Feldman-Cousins delta chi2 after throwing %2% universes is %3%") % __func__ % nuniv % flattened_dchi2s[0.9*flattened_dchi2s.size()]; From 8a16452e497304e74803f833400b8cfe57c5480a Mon Sep 17 00:00:00 2001 From: Nathaniel Date: Tue, 23 Jun 2026 16:19:03 -0500 Subject: [PATCH 3/8] add command line arg --- bin/PROfit.cxx | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/bin/PROfit.cxx b/bin/PROfit.cxx index 7fa87331..f7f79dad 100644 --- a/bin/PROfit.cxx +++ b/bin/PROfit.cxx @@ -479,6 +479,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 @@ -2558,9 +2559,9 @@ int main(int argc, char* argv[]) fc_progress.start_display_thread(); bool save_null_dist = true; - if(pvalue){ - std::string pvalue_file = "chi2_null.root"; - TFile* file = TFile::Open(filename.c_str(), "READ"); + std::string pvalue_file = "chi2_null.root"; + if(reuse_dist){ + TFile* file = TFile::Open(pvalue_file.c_str(), "READ"); if (!file) { save_null_dist = true; } @@ -2570,7 +2571,7 @@ int main(int argc, char* argv[]) } std::vector flattened_dchi2s; - if(save_null_dist || !pvalue){ + if(save_null_dist){ for(size_t i = 0; i < FCthreads; i++) { dchi2s.emplace_back(); outs.emplace_back(); @@ -2588,10 +2589,10 @@ int main(int argc, char* argv[]) for(const auto& v: dchi2s) for(const auto& dchi2: v) flattened_dchi2s.push_back(dchi2); std::sort(flattened_dchi2s.begin(), flattened_dchi2s.end()); - writeVectorToFile(&pvalue_file, flattened_dchi2s); + writeVectorToFile(pvalue_file, flattened_dchi2s); } else if(!save_null_dist){ - flattened_dchi2s = readVectorFromFile(&pvalue_file); + flattened_dchi2s = readVectorFromFile(pvalue_file); } log(L"%1% || 90%% Feldman-Cousins delta chi2 after throwing %2% universes is %3%") % __func__ % nuniv % flattened_dchi2s[0.9*flattened_dchi2s.size()]; From ea75a5c0bd57cd5ce24bc05f1321a095045e1dd5 Mon Sep 17 00:00:00 2001 From: Nathaniel Date: Tue, 23 Jun 2026 16:44:17 -0500 Subject: [PATCH 4/8] fix missing var --- bin/PROfit.cxx | 1 + 1 file changed, 1 insertion(+) diff --git a/bin/PROfit.cxx b/bin/PROfit.cxx index f7f79dad..905dce4a 100644 --- a/bin/PROfit.cxx +++ b/bin/PROfit.cxx @@ -376,6 +376,7 @@ int main(int argc, char* argv[]) PlotBounds pbounds; size_t nuniv; bool gof_pvalue = false; + bool reuse_dist = false; bool pvalue = false; From 98da06422c28e50b8df62352a836959b820b6935 Mon Sep 17 00:00:00 2001 From: Nathaniel Date: Wed, 24 Jun 2026 23:36:25 -0500 Subject: [PATCH 5/8] fitx bug with multithreading --- bin/PROfit.cxx | 128 +++++++++++++++++++++++++------------------------ 1 file changed, 66 insertions(+), 62 deletions(-) diff --git a/bin/PROfit.cxx b/bin/PROfit.cxx index 905dce4a..392f6cc9 100644 --- a/bin/PROfit.cxx +++ b/bin/PROfit.cxx @@ -75,8 +75,8 @@ void writeVectorToFile(const std::string& filename, const std::vector& da // 6. Write and clean up outFile->Write(); - outFile->Close(); - delete outFile; + //outFile->Close(); + //delete outFile; std::cout << "Successfully wrote vector of size " << dataVector.size() << " to " << filename << std::endl; } @@ -96,8 +96,8 @@ std::vector readVectorFromFile(const std::string& filename) { inFile->GetObject("myTree", tree); if (!tree) { std::cerr << "Error: Could not find TTree 'myTree' in " << filename << std::endl; - inFile->Close(); - delete inFile; + //inFile->Close(); + //delete inFile; return resultVector; } @@ -114,8 +114,8 @@ std::vector readVectorFromFile(const std::string& filename) { } // 5. Clean up ROOT memory - inFile->Close(); - delete inFile; + //inFile->Close(); + //delete inFile; return resultVector; } @@ -2556,9 +2556,6 @@ 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 save_null_dist = true; std::string pvalue_file = "chi2_null.root"; if(reuse_dist){ @@ -2570,6 +2567,10 @@ int main(int argc, char* argv[]) save_null_dist = false; } } + if(save_null_dist){ + fc_progress.initialize_display(); + fc_progress.start_display_thread(); + } std::vector flattened_dchi2s; if(save_null_dist){ @@ -2593,6 +2594,7 @@ int main(int argc, char* argv[]) writeVectorToFile(pvalue_file, flattened_dchi2s); } else if(!save_null_dist){ + fc_progress.finish_all(); flattened_dchi2s = readVectorFromFile(pvalue_file); } log(L"%1% || 90%% Feldman-Cousins delta chi2 after throwing %2% universes is %3%") @@ -2623,68 +2625,69 @@ int main(int argc, char* argv[]) log(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_". - // Vector kept alive for the full lifetime of the TTree. - std::vector best_phys_vals(model->nparams, 0.0f); - std::map 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); + if(save_null_dist){ + { + 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_". + // Vector kept alive for the full lifetime of the TTree. + std::vector best_phys_vals(model->nparams, 0.0f); + std::map 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.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; + 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"; } - 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 ************************** @@ -3218,6 +3221,7 @@ int main(int argc, char* argv[]) } log(L"%1% || ########################################################") % __func__; } + log(L"%1% || Before timing info") % __func__; if(global_fit_out.is_open()) global_fit_out.close(); delete metric; From 79c98d234e725a57592f9f283df8a84355d9f623 Mon Sep 17 00:00:00 2001 From: Nathaniel Date: Wed, 15 Jul 2026 15:29:42 -0500 Subject: [PATCH 6/8] Combine with existing FC files --- bin/PROfit.cxx | 272 +++++++++++++++++++++++++++++++++++-------------- 1 file changed, 198 insertions(+), 74 deletions(-) diff --git a/bin/PROfit.cxx b/bin/PROfit.cxx index 392f6cc9..ca8a8685 100644 --- a/bin/PROfit.cxx +++ b/bin/PROfit.cxx @@ -50,7 +50,7 @@ using namespace PROfit; -void writeVectorToFile(const std::string& filename, const std::vector& dataVector) { +/*void writeVectorToFile(const std::string& filename, const std::vector& dataVector) { // 1. Create/Overwrite the specified ROOT file TFile* outFile = TFile::Open(filename.c_str(), "RECREATE"); if (!outFile || outFile->IsZombie()) { @@ -79,9 +79,115 @@ void writeVectorToFile(const std::string& filename, const std::vector& da //delete outFile; std::cout << "Successfully wrote vector of size " << dataVector.size() << " to " << filename << std::endl; +}*/ + +std::vector readAnalysisFromFile(const std::string& filename, + const PROmodel &model, // To parse log10 parameters + const std::vector& variable_systs, // To map spline names to Eigen indices + size_t i_prime) { // Current system config index + std::vector entries; + + // 1. Open the file + TFile* inFile = TFile::Open(filename.c_str(), "READ"); + if (!inFile || inFile->IsZombie()) { + std::cerr << "Error: Could not open file " << filename << std::endl; + return entries; + } + + // 2. Get the TTree + TTree* tree = nullptr; + inFile->GetObject("tree", tree); + if (!tree) { + std::cerr << "Error: Could not find TTree 'tree' in " << filename << std::endl; + inFile->Close(); + delete inFile; + return entries; + } + + // 3. Set up local reading buffers + float chi2_osc = 0; + float chi2_syst = 0; + std::map* p_best_systs_osc = nullptr; + std::map* p_best_systs = nullptr; + std::map* p_syst_throw = nullptr; + + // Allocate memory for physical parameter values matching model size + std::vector phys_param_vals(model.nparams, 0.0f); + + // 4. Bind tree branches to our buffers + tree->SetBranchAddress("chi2_osc", &chi2_osc); + tree->SetBranchAddress("chi2_syst", &chi2_syst); + 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 { + std::cerr << "Warning: Branch " << branch_name << " not found in file." << std::endl; + phys_param_vals[i] = 0.0f; + } + } + + // Get the systematic size details + size_t n_splines = variable_systs[i_prime].GetNSplines(); + + // 5. Loop over all records and reconstruct fc_out structures + Long64_t nEntries = tree->GetEntries(); + entries.reserve(nEntries); + + for (Long64_t entry = 0; entry < nEntries; ++entry) { + tree->GetEntry(entry); + + 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 the parameter was log-scaled, the file has 10^val. We retrieve val = log10(stored_val) + 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; + + return entries; } -std::vector readVectorFromFile(const std::string& filename) { +/*std::vector readVectorFromFile(const std::string& filename) { std::vector resultVector; // 1. Open the specified file @@ -118,7 +224,7 @@ std::vector readVectorFromFile(const std::string& filename) { //delete inFile; return resultVector; -} +}*/ // Unique key for DetVar propeller maps (names can be reused across sections). static std::string DetVarKey(const PROconfig& config, size_t file_index) { @@ -2556,24 +2662,24 @@ 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); - bool save_null_dist = true; - std::string pvalue_file = "chi2_null.root"; + bool gen_null_dist = true; + std::string FC_file = analysis_tag+"_FC.root"; if(reuse_dist){ - TFile* file = TFile::Open(pvalue_file.c_str(), "READ"); + TFile* file = TFile::Open(FC_file.c_str(), "READ"); if (!file) { - save_null_dist = true; + gen_null_dist = true; } else { - save_null_dist = false; + gen_null_dist = false; } } - if(save_null_dist){ + if(gen_null_dist){ fc_progress.initialize_display(); fc_progress.start_display_thread(); } std::vector flattened_dchi2s; - if(save_null_dist){ + if(gen_null_dist){ for(size_t i = 0; i < FCthreads; i++) { dchi2s.emplace_back(); outs.emplace_back(); @@ -2591,11 +2697,29 @@ int main(int argc, char* argv[]) for(const auto& v: dchi2s) for(const auto& dchi2: v) flattened_dchi2s.push_back(dchi2); std::sort(flattened_dchi2s.begin(), flattened_dchi2s.end()); - writeVectorToFile(pvalue_file, flattened_dchi2s); + //writeVectorToFile(FC_file, flattened_dchi2s); } - else if(!save_null_dist){ + else if(!gen_null_dist){ fc_progress.finish_all(); - flattened_dchi2s = readVectorFromFile(pvalue_file); + std::vector flat_outs = readAnalysisFromFile(FC_file, *model, variable_systs, config.i_prime); + // The following code assumes non-flat, so we can just make it fake non-flat + outs.push_back(flat_outs); + float dchi2; + float chi2_osc; + float chi2_syst; + for(const auto &out: outs) { + for(const auto &fco: out) { + chi2_osc = fco.chi2_osc; + chi2_syst = fco.chi2_syst; + if(gof_pvalue){ + dchi2 = chi2_syst; + } + else{ + dchi2 = chi2_syst - chi2_osc; + } + flattened_dchi2s.push_back(dchi2); + } + } } log(L"%1% || 90%% Feldman-Cousins delta chi2 after throwing %2% universes is %3%") % __func__ % nuniv % flattened_dchi2s[0.9*flattened_dchi2s.size()]; @@ -2625,67 +2749,67 @@ int main(int argc, char* argv[]) log(L"%1% || FC Corrected pval after throwing %2% universes is %3%") % __func__ % nuniv % pvalFC ; } - if(save_null_dist){ - { - 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_". - // Vector kept alive for the full lifetime of the TTree. - std::vector best_phys_vals(model->nparams, 0.0f); - std::map 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_". + // Vector kept alive for the full lifetime of the TTree. + std::vector best_phys_vals(model->nparams, 0.0f); + std::map 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"; + } + } + } + //} } //*********************************************************************** From 1b8b59eddeeca480b669cc5ddb1d4f59ccc67f6c Mon Sep 17 00:00:00 2001 From: nathanielerowe <70993723+nathanielerowe@users.noreply.github.com> Date: Thu, 6 Aug 2026 14:17:42 -0500 Subject: [PATCH 7/8] Add catches for reusing files --- bin/PROfit.cxx | 145 +++++++++++++++++++++++++------------------------ 1 file changed, 75 insertions(+), 70 deletions(-) diff --git a/bin/PROfit.cxx b/bin/PROfit.cxx index ca8a8685..d47347e4 100644 --- a/bin/PROfit.cxx +++ b/bin/PROfit.cxx @@ -49,48 +49,17 @@ using namespace PROfit; - -/*void writeVectorToFile(const std::string& filename, const std::vector& dataVector) { - // 1. Create/Overwrite the specified ROOT file - TFile* outFile = TFile::Open(filename.c_str(), "RECREATE"); - if (!outFile || outFile->IsZombie()) { - std::cerr << "Error: Could not create file " << filename << std::endl; - return; - } - - // 2. Create the TTree - TTree* tree = new TTree("myTree", "A tree with a vector"); - - // 3. ROOT needs a pointer to the vector. - // Since dataVector is const&, we make a local pointer to a non-const copy - // or simply point to a local copy to ensure safety. - std::vector localCopy = dataVector; - std::vector* vecPtr = &localCopy; - - // 4. Create the branch pointing to our pointer address - tree->Branch("floatVector", &vecPtr); - - // 5. Fill the tree with the vector data (creates 1 entry containing the entire vector) - tree->Fill(); - - // 6. Write and clean up - outFile->Write(); - //outFile->Close(); - //delete outFile; - - std::cout << "Successfully wrote vector of size " << dataVector.size() << " to " << filename << std::endl; -}*/ - std::vector readAnalysisFromFile(const std::string& filename, const PROmodel &model, // To parse log10 parameters const std::vector& variable_systs, // To map spline names to Eigen indices - size_t i_prime) { // Current system config index + size_t i_prime, // Current system config index + bool gof_mode_requested) { // Requested runtime mode std::vector entries; // 1. Open the file TFile* inFile = TFile::Open(filename.c_str(), "READ"); if (!inFile || inFile->IsZombie()) { - std::cerr << "Error: Could not open file " << filename << std::endl; + log(L"%1% || Error: Could not open file %2%") % __func__ % filename.c_str(); return entries; } @@ -98,13 +67,29 @@ std::vector readAnalysisFromFile(const std::string& filename, TTree* tree = nullptr; inFile->GetObject("tree", tree); if (!tree) { - std::cerr << "Error: Could not find TTree 'tree' in " << filename << std::endl; + log(L"%1% || Error: Could not find TTree 'tree' in %2%") % __func__ % filename.c_str(); inFile->Close(); delete inFile; return entries; } - // 3. Set up local reading buffers + // 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(L"%1% || ERROR: The reused file '%2%' is missing oscillation branches ('chi2_osc' / 'best_systs_osc').") + % __func__ % filename.c_str(); + log(L"%1% || Cannot compute standard Feldman-Cousins p-values (--pval) without oscillation fit data.") % __func__; + log(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* p_best_systs_osc = nullptr; @@ -114,10 +99,10 @@ std::vector readAnalysisFromFile(const std::string& filename, // Allocate memory for physical parameter values matching model size std::vector phys_param_vals(model.nparams, 0.0f); - // 4. Bind tree branches to our buffers - tree->SetBranchAddress("chi2_osc", &chi2_osc); + // 5. Bind tree branches to our buffers + if (b_osc) tree->SetBranchAddress("chi2_osc", &chi2_osc); tree->SetBranchAddress("chi2_syst", &chi2_syst); - tree->SetBranchAddress("best_systs_osc", &p_best_systs_osc); + 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); @@ -127,7 +112,7 @@ std::vector readAnalysisFromFile(const std::string& filename, if (tree->GetBranch(branch_name.c_str())) { tree->SetBranchAddress(branch_name.c_str(), &phys_param_vals[i]); } else { - std::cerr << "Warning: Branch " << branch_name << " not found in file." << std::endl; + log(L"%1% || Warning: Branch %2% not found in file.") % __func__ % branch_name.c_str(); phys_param_vals[i] = 0.0f; } } @@ -135,13 +120,24 @@ std::vector readAnalysisFromFile(const std::string& filename, // Get the systematic size details size_t n_splines = variable_systs[i_prime].GetNSplines(); - // 5. Loop over all records and reconstruct fc_out structures + // 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; @@ -150,7 +146,6 @@ std::vector readAnalysisFromFile(const std::string& filename, fco.best_phys_osc.resize(model.nparams); for (size_t i = 0; i < model.nparams; ++i) { float val = phys_param_vals[i]; - // If the parameter was log-scaled, the file has 10^val. We retrieve val = log10(stored_val) if (model.is_log10[i]) { fco.best_phys_osc(i) = (val > 0.0f) ? std::log10(val) : 0.0f; } else { @@ -184,6 +179,21 @@ std::vector readAnalysisFromFile(const std::string& filename, inFile->Close(); delete inFile; + // 7. Validate oscillation data contents against requested mode + if (!gof_mode_requested && !non_zero_osc_found) { + log(L"%1% || ERROR: The reused file '%2%' contains only zeroed oscillation fits (chi2_osc == 0).") + % __func__ % filename.c_str(); + log(L"%1% || This occurs when reusing a file generated with --gof.") % __func__; + log(L"%1% || Cannot compute standard Feldman-Cousins p-values (--pval) from a GOF-only distribution.") % __func__; + log(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(L"%1% || Reusing full Feldman-Cousins distribution for GOF evaluation (extracting chi2_syst).") % __func__; + } + return entries; } @@ -2694,40 +2704,35 @@ int main(int argc, char* argv[]) t.join(); } fc_progress.finish_all(); - - for(const auto& v: dchi2s) for(const auto& dchi2: v) flattened_dchi2s.push_back(dchi2); - std::sort(flattened_dchi2s.begin(), flattened_dchi2s.end()); - //writeVectorToFile(FC_file, flattened_dchi2s); - } - else if(!gen_null_dist){ + } + else if(!gen_null_dist){ fc_progress.finish_all(); - std::vector flat_outs = readAnalysisFromFile(FC_file, *model, variable_systs, config.i_prime); - // The following code assumes non-flat, so we can just make it fake non-flat + std::vector flat_outs = readAnalysisFromFile(FC_file, *model, variable_systs, config.i_prime, gof_pvalue); + nuniv = flat_outs.size(); + if(flat_outs.empty()) { + log(L"%1% || Aborting execution due to invalid or unreadable Feldman-Cousins distribution.") % __func__; + return 1; + } + outs.push_back(flat_outs); - float dchi2; - float chi2_osc; - float chi2_syst; - for(const auto &out: outs) { - for(const auto &fco: out) { - chi2_osc = fco.chi2_osc; - chi2_syst = fco.chi2_syst; - if(gof_pvalue){ - dchi2 = chi2_syst; - } - else{ - dchi2 = chi2_syst - chi2_osc; - } - flattened_dchi2s.push_back(dchi2); - } + } + // Unify array flattening across fresh generation and reused distributions + std::vector 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); } - } - log(L"%1% || 90%% Feldman-Cousins delta chi2 after throwing %2% universes is %3%") + } + std::sort(flattened_dchi2s.begin(), flattened_dchi2s.end()); + std::sort(flattened_syst_chi2.begin(), flattened_syst_chi2.end()); + + log(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 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(L"%1% || All: %2% ") % __func__ % flattened_syst_chi2; log(L"%1% || chi: %2% ") % __func__ % global_chi2; auto it = std::lower_bound(flattened_syst_chi2.begin(), flattened_syst_chi2.end(), global_chi2); From dad86ba047cc1cab64af52fee74eb428e6a5be13 Mon Sep 17 00:00:00 2001 From: nathanielerowe <70993723+nathanielerowe@users.noreply.github.com> Date: Fri, 7 Aug 2026 12:33:28 -0500 Subject: [PATCH 8/8] cleaning --- bin/PROfit.cxx | 40 ---------------------------------------- 1 file changed, 40 deletions(-) diff --git a/bin/PROfit.cxx b/bin/PROfit.cxx index d47347e4..c263834c 100644 --- a/bin/PROfit.cxx +++ b/bin/PROfit.cxx @@ -197,45 +197,6 @@ std::vector readAnalysisFromFile(const std::string& filename, return entries; } -/*std::vector readVectorFromFile(const std::string& filename) { - std::vector resultVector; - - // 1. Open the specified file - TFile* inFile = TFile::Open(filename.c_str(), "READ"); - if (!inFile || inFile->IsZombie()) { - std::cerr << "Error: Could not open file " << filename << std::endl; - return resultVector; // Returns empty vector - } - - // 2. Get the TTree - TTree* tree = nullptr; - inFile->GetObject("myTree", tree); - if (!tree) { - std::cerr << "Error: Could not find TTree 'myTree' in " << filename << std::endl; - //inFile->Close(); - //delete inFile; - return resultVector; - } - - // 3. Set up the pointer for ROOT tracking - std::vector* vecPtr = nullptr; - tree->SetBranchAddress("floatVector", &vecPtr); - - // 4. Read the first entry (assuming 1 entry containing your saved vector) - if (tree->GetEntries() > 0) { - tree->GetEntry(0); - if (vecPtr) { - resultVector = *vecPtr; // Copy the data to our return vector - } - } - - // 5. Clean up ROOT memory - //inFile->Close(); - //delete inFile; - - return resultVector; -}*/ - // 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]; @@ -3350,7 +3311,6 @@ int main(int argc, char* argv[]) } log(L"%1% || ########################################################") % __func__; } - log(L"%1% || Before timing info") % __func__; if(global_fit_out.is_open()) global_fit_out.close(); delete metric;