From 51e5f7c4e5d83239d71f06993e2b0b431fa55e62 Mon Sep 17 00:00:00 2001 From: edbarriga Date: Mon, 23 Mar 2026 10:00:28 -0400 Subject: [PATCH 1/2] Rework of vecps_plotter --- .../vecps_plotter/vecps_plotter.cc | 1338 +++++++++++++---- 1 file changed, 1015 insertions(+), 323 deletions(-) diff --git a/src/programs/AmplitudeAnalysis/vecps_plotter/vecps_plotter.cc b/src/programs/AmplitudeAnalysis/vecps_plotter/vecps_plotter.cc index b4efa15f4..b10dec0d7 100644 --- a/src/programs/AmplitudeAnalysis/vecps_plotter/vecps_plotter.cc +++ b/src/programs/AmplitudeAnalysis/vecps_plotter/vecps_plotter.cc @@ -2,6 +2,7 @@ #include #include #include +#include #include "AMPTOOLS_AMPS/BreitWigner.h" #include "AMPTOOLS_AMPS/ComplexCoeff.h" @@ -50,367 +51,1029 @@ void atiSetup() { AmpToolsInterface::registerDataReader(ROOTDataReaderTEM()); } -using namespace std; +// ...oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... +// Enum Classes to Simplify Boolean Logic +// ...oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... -int main(int argc, char* argv[]) { - // ************************ - // usage - // ************************ +// Classify the sum of the reaction according to reflectivity +enum class ReflSum { + Neither, // something like an incoherent background sum + Positive, + Negative +}; - cout << endl - << " *** Viewing Results Using AmpPlotter and writing root histograms " - "*** " - << endl - << endl; +// Specify the group of amplitudes to be plotted +enum class AmpPlotGroup{ + AllAmps, // for the total fit result + AllReflAmps, // for the total refl + JlFamily, // sum over the m projections + SingleAmp // one amplitude at a time +}; - if (argc < 2) { - cout << "Usage:" << endl << endl; - cout << "\tvecps_plotter -o " - << endl - << endl; - return 0; - } +// Specify the group of sums to be plotted +enum class SumPlotGroup{ + AllSums, // for the total fit results + BothReflSums, // Incoherent sum of both reflectivities + PosReflSum, // one coherent sum + NegReflSum, // one coherent sum + NoReflSum // This could be broken down for multiple non-refl sums +}; +// ...oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... +// Structures to Gather Relevant Info Together +// ...oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... + +// Store the information from the command line arguments and configuration file +// and allows to easily change the default values of the options +// TODO: add an option for only GUI mode (aka no root file output) +struct ProgramOptions{ + std::string resultsName; + std::string outName = "vecps_plot.root"; bool showGui = false; - string outName = "vecps_plot.root"; - string resultsName(argv[1]); - for (int i = 2; i < argc; i++) { - string arg(argv[i]); + bool mcGen = false; + bool bkgFile = false; + bool sFile = false; + bool skipPlots = false; + bool txtOutput = false; + bool csvOutput = false; +}; - if (arg == "-g") { - showGui = true; - } - if (arg == "-o") { - outName = argv[++i]; - } - if (arg == "-h") { - cout << endl << " Usage for: " << argv[0] << endl << endl; - cout << "\t -o \t output file path" << endl; - cout << "\t -g \t show GUI" << endl; - exit(1); - } +// Store the information needed to identify a wave +struct WaveKey { + std::string ampName; // full amp name (e.g. 1S+1) + std::string jFamily; // Jl family (e.g. 1S) + std::string m; // m projection (e.g. +1) + bool isNumeric; // non-Jl amps do not start with a number +}; + +// Store all the information that identifies an amplitude +struct AmplitudeKey { + std::string fullName; // unparsed AmpTools full name + std::string reactionName; + std::string sumName; + ReflSum reflSum; // Pos/neg refl or neither + WaveKey key; + + // index in uniqueAmplitudes() + size_t ampIndex = std::numeric_limits::max(); + // index in uniqueSums() + size_t sumIndex = std::numeric_limits::max(); +}; + +// Store all the relation between the unique sums and unique amps +// could be extended to global (between reactions) relations +struct Sisterhood{ + // maps to look up the indices of the unique amplitudes and sums + // the key is the string/enum name of the amplitude/sum + std::map> jlIndexMap; + std::map> sumIndexMap; + + // Relationship between sums and each amplitude (index-based) + std::vector> ampToSums; + std::vector> sumToAmps; + // Note: amplitudes that do not have refl are not in the following map + std::map> reflToAmps; + + // Relations for amplitudes based on the full name + std::vector> ampToFull; // uniqueAmp to full indices + std::vector> sumToFull; // uniqueSum to full indices + std::map> reflToFull; + std::map, std::vector> ampReflToFull; + std::map, std::vector> jlReflToFull; +}; + +// Store the information that is static when a reaction is read +// this structure avoids having to parse the string information again +// and it becomes a look up table for the indices of the amplitudes and +// sums when plotting +struct AmplitudeRegistry { + // Collection of all of the amplitudes in the fit file + std::vector ampsList; + + // Parsed unique amplitudes + std::vector uniqueAmps; + size_t nUnqAmps; + // Parsed unique sums + std::vector uniqueSums; + size_t nUnqSums; + // containers to look up the indices of the unique amplitudes and sums + // and their relation to other sums/amps + Sisterhood relations; +}; + +// This stores the specific configuration for a plotting +struct PlotConfig{ + AmpPlotGroup ampGroup; + SumPlotGroup sumGroup; + + // Store the indices to be enable + std::vector ampIndices; + std::vector sumIndices; +}; + +struct PlotLabels{ + std::string human; // human friendly output + std::string histName; // better for saving in a ROOT file + // TODO: add csv format + // std::string csvFormat; // better for saving in a CSV file +}; + +// ...oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... +// Helper Functions to Fill the Structures, etc +// ...oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... + +// Parse the name of the sum to classify it according to reflectivity +// Note: this follow standard naming conventions for the vec_ps amplitudes +// this means that: +// [ImagNegSign & RealPosSign] = Negative Reflectivity +// [RealNegSign & ImagPosSign] = Positive Reflectivity +// The 3rd option is "Neither", but it can be easily expanded if needed +// dont forget to update the ReflSum enum if you do! +// General naming convention for amps: +// https://halldweb.jlab.org/doc-private/DocDB/ShowDocument?docid=7036 +ReflSum classifyReflSum(const std::string& sumName){ + if (sumName.find("RealNegSign") != std::string::npos || + sumName.find("ImagPosSign") != std::string::npos){ + return ReflSum::Positive; } + else if (sumName.find("RealPosSign") != std::string::npos || + sumName.find("ImagNegSign") != std::string::npos){ + return ReflSum::Negative; + } + else { + return ReflSum::Neither; + } +} - // ************************ - // parse the command line parameters - // ************************ +std::string reflToString(ReflSum r){ + switch(r){ + case ReflSum::Positive: return "PosRefl"; + case ReflSum::Negative: return "NegRefl"; + default: return "Neither"; + } +} - cout << "Fit results file name = " << resultsName << endl; - cout << "Output file name = " << outName << endl << endl; +// Store in a vector the different parts of the full amplitude name, given +// the common delimiter used in AmpTools (::) +std::vector breakDownName(const std::string& fullAmpName, + const std::string& delim){ + std::vector brokenDownName; + size_t start = 0; size_t end; - // ************************ - // load the results and display the configuration info - // ************************ - cout << "Initializing AmpTools interface..." << endl; - atiSetup(); - cout << "Interface initialized." << endl; + while((end = fullAmpName.find(delim, start)) != std::string::npos){ + brokenDownName.push_back(fullAmpName.substr(start, end-start)); + start = end + delim.size(); + } - cout << "Loading Fit results" << endl; - FitResults results(resultsName); - if (!results.valid()) { - cout << "Invalid fit results in file: " << resultsName << endl; - exit(1); + brokenDownName.push_back(fullAmpName.substr(start)); + return brokenDownName; +} + +WaveKey getJlm(const std::string& ampName){ + WaveKey key; + key.ampName = ampName; + // For the standard format "Jl.m" (e.g. "1S+1") + // it splits the m-projection and the Jl family + key.isNumeric = std::isdigit(ampName[0]); + if (!ampName.empty() && key.isNumeric){ + key.jFamily = ampName.substr(0, 2); + key.m = ampName.substr(2); } - cout << "Fit results loaded" << endl; - // ************************ - // set up the plot generator - // ************************ - cout << "Plotgen results" << endl; - vecps_PlotGen plotGen(results, PlotGenerator::kNoGenMC); - cout << " Initialized ati and PlotGen" << endl; - - // ************************ - // set up an output ROOT file to store histograms - // ************************ - // Loop over different polarization files - // the default will remain to just use a single file - // ************************ - size_t nReactions = results.reactionList().size(); - for (int polFile = 0; polFile < nReactions; polFile++) { - if (polFile > 0) break; // TEMPORARY - REMOVE TO ENABLE ALL POL FILES - string reactionName = results.reactionList()[polFile]; - outName = reactionName + ".root"; - - TFile* plotfile = new TFile(outName.c_str(), "recreate"); - TH1::AddDirectory(kFALSE); - - plotGen.enableReaction(reactionName); - vector sums = plotGen.uniqueSums(); - vector amps = plotGen.uniqueAmplitudes(); - cout << "Reaction " << reactionName << " enabled with " << sums.size() - << " sums and " << amps.size() << " amplitudes" << endl; - - const vector amphistname = {"1S-1", "1S+0", "1S+1", "1S", - "1P-1", "1P+0", "1P+1", "1P", - "1D-1", "1D+0", "1D+1", "1D"}; - const vector reflname = {"PosRefl", "NegRefl"}; - - // loop over sum configurations (one for each of the individual - // contributions [when irefl < reflname.size()], and the combined sum - // of all [when irefl == reflname.size()]) - for (unsigned int irefl = 0; irefl <= reflname.size(); irefl++) { - // loop over desired amplitudes (one for each of the individual - // contributions [when iamp < amphistname.size()], and the combined - // sum of all [when iamp == amphistname.size()]) - for (unsigned int iamp = 0; iamp <= amphistname.size(); iamp++) { - // turn all ampltiudes by default - for (unsigned int jamp = 0; jamp < amps.size(); jamp++) { - plotGen.enableAmp(jamp); - } + else { + key.jFamily = ""; + key.m = ""; + } + return key; +} - // turn off unwanted amplitudes and sums - if (iamp < amphistname.size()) { - string locampname = amphistname[iamp]; - - // turn on all sums by default - for (unsigned int i = 0; i < sums.size(); i++) - plotGen.enableSum(i); - - // turn off unwanted sums for reflectivity (based on - // naturality) cout<<"refl = "< tokens = breakDownName(fullAmpName, "::"); + + // Expect exactly 3 tokens: + // [reaction] :: [sum] :: [ampName] + if (tokens.size() != 3){ + std::cerr << "Warning: Unexpected amplitude name format: " + << fullAmpName << "\n"; + ampKey.reactionName = ""; + ampKey.sumName = ""; + ampKey.key.jFamily = ""; + ampKey.key.m = ""; + } + else { + ampKey.reactionName = tokens[0]; + ampKey.sumName = tokens[1]; + ampKey.key = getJlm(tokens[2]); + } + ampKey.reflSum = classifyReflSum(ampKey.sumName); + + return ampKey; +} + +// Just a reminder that the registry is meant to be once per file, and +// it serves as a look up table for everything related to the amps +AmplitudeRegistry buildRegistry(std::unique_ptr& plotResults){ + AmplitudeRegistry registry; + + std::vector uniqueAmps = plotResults->uniqueAmplitudes(); + registry.nUnqAmps = uniqueAmps.size(); + std::vector uniqueSums = plotResults->uniqueSums(); + registry.nUnqSums = uniqueSums.size(); + + // Look up maps for parsed info + std::map ampNameToIndex; + std::map sumNameToIndex; + + std::vector uniqueAmpsParsed; + for (size_t i = 0; i < uniqueAmps.size(); ++i){ + ampNameToIndex[uniqueAmps[i]] = i; + WaveKey key = getJlm(uniqueAmps[i]); + uniqueAmpsParsed.push_back(key); + if(!key.jFamily.empty()) // non-numeric amps have and empty jFamily + registry.relations.jlIndexMap[key.jFamily].push_back(i); + } + registry.uniqueAmps = uniqueAmpsParsed; + + std::vector uniqueSumsParsed; + for (size_t i = 0; i < uniqueSums.size(); ++i){ + sumNameToIndex[uniqueSums[i]] = i; + ReflSum refl = classifyReflSum(uniqueSums[i]); + uniqueSumsParsed.push_back(refl); + registry.relations.sumIndexMap[refl].push_back(i); + } + registry.uniqueSums = uniqueSumsParsed; + + std::vector parsedAmps; + for (const auto& amp : plotResults->fullAmplitudes()){ + AmplitudeKey ampKey = parseAmp(amp); + ampKey.ampIndex = ampNameToIndex[ampKey.key.ampName]; + ampKey.sumIndex = sumNameToIndex[ampKey.sumName]; + ampKey.reflSum = classifyReflSum(ampKey.sumName); + parsedAmps.push_back(ampKey); + } + registry.ampsList = parsedAmps; + + // Sisterhood vectors + std::vector> ampToSums(registry.nUnqAmps); + std::vector> sumToAmps(registry.nUnqSums); + std::map> reflToAmps; + + // Sisterhood vectors for full reference + std::vector> ampToFull(registry.nUnqAmps); + std::vector> sumToFull(registry.nUnqSums); + std::map> reflToFull; + std::map, std::vector> ampReflToFull; + std::map, std::vector> jlReflToFull; + + for(size_t iFull = 0; iFull < registry.ampsList.size(); ++iFull){ + const auto& amp = registry.ampsList[iFull]; + size_t ampIndex = amp.ampIndex; + size_t sumIndex = amp.sumIndex; + + auto& amp2Sum = ampToSums[ampIndex]; + if(std::find(amp2Sum.begin(), amp2Sum.end(), sumIndex) == amp2Sum.end()) + amp2Sum.push_back(sumIndex); + + auto& sum2Amp = sumToAmps[sumIndex]; + if(std::find(sum2Amp.begin(), sum2Amp.end(), ampIndex) == sum2Amp.end()) + sum2Amp.push_back(ampIndex); + + auto& refl2Amp = reflToAmps[amp.reflSum]; + if(std::find(refl2Amp.begin(), refl2Amp.end(), ampIndex) == refl2Amp.end()) + refl2Amp.push_back(ampIndex); + + ampToFull[amp.ampIndex].push_back(iFull); + sumToFull[amp.sumIndex].push_back(iFull); + reflToFull[amp.reflSum].push_back(iFull); + + ampReflToFull[{ampIndex, amp.reflSum}].push_back(iFull); + if(!amp.key.jFamily.empty()) + jlReflToFull[{amp.key.jFamily, amp.reflSum}].push_back(iFull); + } - // turn off unwanted amplitudes - for (unsigned int jamp = 0; jamp < amps.size(); jamp++) { - if (amps[jamp].find(locampname.data()) == - std::string::npos) { - plotGen.disableAmp(jamp); - // cout<<"disable amplitude "<& amps, + const std::vector& sums, + const Sisterhood& rel){ + if(amps.empty() || sums.empty()) + return false; + + std::unordered_set sumSet(sums.begin(), sums.end()); + for(size_t a : amps){ + for(size_t s : rel.ampToSums[a]){ + if(sumSet.count(s)) + return true; + } + } + return false; +} + +// The vector that this function outputs contains all the configurations that +// are allowed by the enum classes defined above with the caveat that we +// screen for "physically meaningful" combinations in isNonEmptyConfig() +std::vector generatePlotConfigs(const AmplitudeRegistry& reg){ + std::vector configs; + std::cout << reg.uniqueAmps.size() << " unique amps and " + << reg.uniqueSums.size() << " unique sums." << std::endl; + + const std::vector ampGroups = { + AmpPlotGroup::AllAmps, + AmpPlotGroup::AllReflAmps, + AmpPlotGroup::SingleAmp, + AmpPlotGroup::JlFamily + }; + + const std::vector sumGroups = { + SumPlotGroup::AllSums, + SumPlotGroup::PosReflSum, + SumPlotGroup::NegReflSum, + SumPlotGroup::NoReflSum, + SumPlotGroup::BothReflSums + }; + + std::vector allAmpIndices(reg.nUnqAmps); + std::iota(allAmpIndices.begin(), allAmpIndices.end(), 0); + + std::vector allSumIndices(reg.nUnqSums); + std::iota(allSumIndices.begin(), allSumIndices.end(), 0); + + // Extract the relation from the registry + const auto& rel = reg.relations; + + // These vector is meant to be used with AllReflAmps. Note that generally + // both Positive and Negative share the same wave. But there are situations + // where only one refl is present. So, it's safer to loop over both and + // then remove duplicates + std::vector reflAmps; + if(rel.reflToAmps.count(ReflSum::Positive)) + reflAmps.insert(reflAmps.end(), + rel.reflToAmps.at(ReflSum::Positive).begin(), + rel.reflToAmps.at(ReflSum::Positive).end()); + + if(rel.reflToAmps.count(ReflSum::Negative)) + reflAmps.insert(reflAmps.end(), + rel.reflToAmps.at(ReflSum::Negative).begin(), + rel.reflToAmps.at(ReflSum::Negative).end()); + + // Remove duplicates + std::sort(reflAmps.begin(), reflAmps.end()); + reflAmps.erase(std::unique(reflAmps.begin(), reflAmps.end()), reflAmps.end()); + + // To build all possible combinations we need to loop according to the + // following logic: (AmpPlotGroup × SumPlotGroup) × (ampSets × sumSets) + // the first group, (AmpPlotGroup × SumPlotGroup), tell us what type of + // plot we want; while the second group, (ampSets × sumSets), tell us which + // amplitudes and sums to include in the plot + for(auto aGroup : ampGroups){ + + std::vector> ampSets; + + // build ampSets + switch(aGroup){ + case AmpPlotGroup::AllAmps: + ampSets.push_back(allAmpIndices); + break; + + case AmpPlotGroup::AllReflAmps: + ampSets.push_back(reflAmps); + break; + + case AmpPlotGroup::SingleAmp: + for(size_t i : allAmpIndices) + ampSets.push_back({i}); + break; + + case AmpPlotGroup::JlFamily: + for(const auto& [jl, indices] : rel.jlIndexMap) + ampSets.push_back(indices); + break; + } + + for(auto sGroup : sumGroups){ + // AmpPlotGroup::AllAmps can only be with SumPlotGroup::AllSums + if ((aGroup == AmpPlotGroup::AllAmps) != + (sGroup == SumPlotGroup::AllSums)) + continue; + + // AmpPlotGroup::AllReflAmps can only be with + // SumPlotGroup::PosReflSum, SumPlotGroup::NegReflSum, or + // SumPlotGroup::BothReflSums + if (aGroup == AmpPlotGroup::AllReflAmps) { + if (!(sGroup == SumPlotGroup::PosReflSum || + sGroup == SumPlotGroup::NegReflSum || + sGroup == SumPlotGroup::BothReflSums)) + continue; + } + + std::vector> sumSets; + + // build sumSets + switch(sGroup){ + case SumPlotGroup::AllSums: + sumSets.push_back(allSumIndices); + break; + + case SumPlotGroup::PosReflSum: + if(rel.sumIndexMap.count(ReflSum::Positive)) + sumSets.push_back(rel.sumIndexMap.at(ReflSum::Positive)); + break; + + case SumPlotGroup::NegReflSum: + if(rel.sumIndexMap.count(ReflSum::Negative)) + sumSets.push_back(rel.sumIndexMap.at(ReflSum::Negative)); + break; + + case SumPlotGroup::NoReflSum: + if(rel.sumIndexMap.count(ReflSum::Neither)) + sumSets.push_back(rel.sumIndexMap.at(ReflSum::Neither)); + break; + + case SumPlotGroup::BothReflSums:{ + std::vector posNeg; + + if(rel.sumIndexMap.count(ReflSum::Positive)){ + posNeg.insert(posNeg.end(), + rel.sumIndexMap.at(ReflSum::Positive).begin(), + rel.sumIndexMap.at(ReflSum::Positive).end()); + } + if(rel.sumIndexMap.count(ReflSum::Negative)){ + posNeg.insert(posNeg.end(), + rel.sumIndexMap.at(ReflSum::Negative).begin(), + rel.sumIndexMap.at(ReflSum::Negative).end()); } + sumSets.push_back(posNeg); + break; } + } - cout << "Looping over input data" << endl; - // loop over data, accMC, and genMC - for (unsigned int iplot = 0; iplot < PlotGenerator::kNumTypes; - iplot++) { - if (iplot == PlotGenerator::kGenMC || - iplot == PlotGenerator::kBkgnd) - continue; - bool singleData = - irefl == reflname.size() && iamp == amphistname.size(); - if (iplot == PlotGenerator::kData && !singleData) - continue; // only plot data once - - // loop over different variables - for (unsigned int ivar = 0; - ivar < VecPsPlotGenerator::kNumHists; ivar++) { - // set unique histogram name for each plot (could put in - // directories...) - string histname = ""; - if (ivar == VecPsPlotGenerator::kVecPsMass) - histname += "MVecPs"; - else if (ivar == VecPsPlotGenerator::kCosTheta) - histname += "CosTheta"; - else if (ivar == VecPsPlotGenerator::kPhi) - histname += "Phi"; - else if (ivar == VecPsPlotGenerator::kCosThetaH) - histname += "CosTheta_H"; - else if (ivar == VecPsPlotGenerator::kPhiH) - histname += "Phi_H"; - else if (ivar == VecPsPlotGenerator::kProd_Ang) - histname += "Prod_Ang"; - else if (ivar == VecPsPlotGenerator::kt) - histname += "t"; - else if (ivar == VecPsPlotGenerator::kRecoilMass) - histname += "MRecoil"; - else if (ivar == VecPsPlotGenerator::kProtonPsMass) - histname += "MProtonPs"; - else if (ivar == VecPsPlotGenerator::kRecoilPsMass) - histname += "MRecoilPs"; - else if (ivar == VecPsPlotGenerator::kLambda) - histname += "Lambda"; - else if (ivar == VecPsPlotGenerator::kDalitz) - histname += "Dalitz"; - else if (ivar == VecPsPlotGenerator::kPhi_ProdVsPhi) - histname += "Phi_ProdVsPhi"; - else if (ivar == VecPsPlotGenerator::kPhiOffsetVsPhi) - histname += "PhiOffsetVsPhi"; - else - continue; - - if (iplot == PlotGenerator::kData) histname += "dat"; - if (iplot == PlotGenerator::kBkgnd) histname += "bkgnd"; - if (iplot == PlotGenerator::kAccMC) histname += "acc"; - if (iplot == PlotGenerator::kGenMC) histname += "gen"; - - if (irefl < reflname.size()) { - // get name of sum for naming histogram - string sumName = reflname[irefl]; - histname += "_"; - histname += sumName; - } - if (iamp < amphistname.size()) { - // get name of amp for naming histogram - histname += "_"; - histname += amphistname[iamp]; - } - - Histogram* hist = - plotGen.projection(ivar, reactionName, iplot); - TH1* thist = hist->toRoot(); - thist->SetName(histname.c_str()); - plotfile->cd(); - thist->Write(); - } + // Use the sets to build amp x sum combinations + for(const auto& amps : ampSets){ + for(const auto& sums : sumSets){ + if(!isNonEmptyConfig(amps, sums, rel)) continue; + configs.push_back({aGroup, sGroup, amps, sums}); } } } + } + return configs; +} + +// Rather than reseting the active amps/sums, +// this function checks the active status for them and changes what's needed. +void applyConfig(const PlotConfig& newConfig, const AmplitudeRegistry& reg, + std::unique_ptr& plotResults){ + // Initialize a vector with all entries as false and flip to true the + // indices for the amps that need to be active + std::vector newAmps(reg.nUnqAmps,false); + for(auto a : newConfig.ampIndices) newAmps[a] = true; + + for(size_t i=0;iisAmpEnabled(i) && !newAmps[i]) + plotResults->disableAmp(i); + if(!plotResults->isAmpEnabled(i) && newAmps[i]) + plotResults->enableAmp(i); + } - plotfile->Close(); + // Initialize a vector with all entries as false and flip to true the + // indices for the sums that need to be active + std::vector newSums(reg.nUnqSums,false); + for(auto s : newConfig.sumIndices) newSums[s] = true; - // The next calculations only need to happen for the first file - // and we do not need to repeat for each polarization - if (polFile > 0) continue; - - // model parameters - cout << "Checking Parameters" << endl; + for(size_t i=0;iisSumEnabled(i) && !newSums[i]) + plotResults->disableSum(i); + if(!plotResults->isSumEnabled(i) && newSums[i]) + plotResults->enableSum(i); + } +} - // parameters to check - vector pars; +PlotLabels makeLabel(const PlotConfig& cfg, const AmplitudeRegistry& reg){ + PlotLabels labels; + std::string amp; + std::string sum; - pars.push_back("dsratio"); + switch(cfg.ampGroup){ + case AmpPlotGroup::AllAmps: + amp = "All"; + break; - // file for writing parameters (later switch to putting in ROOT file) - ofstream outfile; - outfile.open("vecps_fitPars.txt"); + case AmpPlotGroup::AllReflAmps: + amp = "AllRefl"; + break; - for (unsigned int i = 0; i < pars.size(); i++) { - double parValue = results.parValue(pars[i]); - double parError = results.parError(pars[i]); - outfile << parValue << "\t" << parError << "\t" << endl; + case AmpPlotGroup::SingleAmp:{ + const auto& ampKey = reg.uniqueAmps[cfg.ampIndices[0]]; + amp = ampKey.ampName; + break; } - outfile << "TOTAL EVENTS = " << results.intensity().first << " +- " - << results.intensity().second << endl; - vector fullamps = plotGen.fullAmplitudes(); - for (unsigned int i = 0; i < fullamps.size(); i++) { - vector useamp; - useamp.push_back(fullamps[i]); - outfile << "FIT FRACTION " << fullamps[i] << " = " - << results.intensity(useamp).first / - results.intensity().first - << " +- " - << results.intensity(useamp).second / - results.intensity().first - << endl; + case AmpPlotGroup::JlFamily:{ + const auto& ampKey = reg.uniqueAmps[cfg.ampIndices[0]]; + amp = ampKey.jFamily; + break; } + } - const int nAmps = amphistname.size(); - vector ampsumPosRefl[nAmps]; - vector ampsumNegRefl[nAmps]; - vector > phaseDiffNames; + switch(cfg.sumGroup){ + case SumPlotGroup::AllSums: + sum = "Full"; + break; - for (unsigned int i = 0; i < fullamps.size(); i++) { - // combine amplitudes with names defined above - for (int iamp = 0; iamp < nAmps; iamp++) { - string locampname = amphistname[iamp]; + case SumPlotGroup::PosReflSum: + sum = "posRefl"; + break; - if (fullamps[i].find("::" + locampname) == std::string::npos) - continue; - // cout<& plotResults, + const std::string& reactionName, + const TFile& plotFile, bool& singleData, + const PlotConfig& cfg, const PlotLabels& cfgLabel){ + // std::cout << "Creating plots for: " << cfgLabel.human << std::endl; + // Loop over data types + for (unsigned int iType = 0; iType < PlotGenerator::kNumTypes; ++iType){ + if (!options.bkgFile && iType == PlotGenerator::kBkgnd) continue; + + if (iType == PlotGenerator::kData){ + if(!singleData) singleData = true; + else if(singleData) continue; // only plot data once + } + + std::string dataType; + switch (iType) { + case PlotGenerator::kData: dataType = "dat"; break; + case PlotGenerator::kAccMC: dataType = "acc"; break; + case PlotGenerator::kGenMC: dataType = "gen"; break; + case PlotGenerator::kBkgnd: dataType = "bkgnd"; break; + default: continue; + } + + // Loop over all histograms defined in VecPsPlotGenerator + for (unsigned int ivar = 0; ivar < VecPsPlotGenerator::kNumHists; ++ivar){ + Histogram* hist = plotResults->projection(ivar, reactionName, iType); + if (!hist) continue; // safety skip + TH1* histRoot = hist->toRoot(); + + std::string histName = ""; + // Assign proper names to each histogram. Note that the names are + // defined in VecPsPlotGenerator + // TODO: add a parser for the variables similar to reflToString to + // VecPsPlotGenerator to make this more automatic and + // avoid the need to edit this file when adding new variables + switch (ivar) { + case VecPsPlotGenerator::kVecPsMass: histName = "MVecPs"; break; + case VecPsPlotGenerator::kCosTheta: histName = "CosTheta"; break; + case VecPsPlotGenerator::kPhi: histName = "Phi"; break; + case VecPsPlotGenerator::kCosThetaH: histName = "CosTheta_H"; break; + case VecPsPlotGenerator::kPhiH: histName = "Phi_H"; break; + case VecPsPlotGenerator::kProd_Ang: histName = "Prod_Ang"; break; + case VecPsPlotGenerator::kt: histName = "t"; break; + case VecPsPlotGenerator::kRecoilMass: histName = "MRecoil"; break; + case VecPsPlotGenerator::kProtonPsMass: histName = "MProtonPs"; break; + case VecPsPlotGenerator::kRecoilPsMass: histName = "MRecoilPs"; break; + case VecPsPlotGenerator::kLambda: histName = "Lambda"; break; + case VecPsPlotGenerator::kDalitz: histName = "Dalitz"; break; + case VecPsPlotGenerator::kPhi_ProdVsPhi: histName = "Phi_ProdVsPhi"; break; + case VecPsPlotGenerator::kPhiOffsetVsPhi: histName = "PhiOffsetVsPhi"; break; + // Add more variables here if needed + default: continue; } + + if(!(iType==PlotGenerator::kData)) + histName += "_" + dataType + cfgLabel.histName; + histRoot->SetName(histName.c_str()); + histRoot->Write(); } + } + +} + +// This function output a semi-human readable text file with the fit fractions +// and phase differences between amplitudes. It has some patterns to help +// parsing, but csv file might be preferred for that purpose. +// Note: the default behavior is the acceptance corrected intensities +// TODO: include things like d/s ratios and other parameters +void writeAmpInfo(const AmplitudeRegistry& reg, + std::unique_ptr& results, + std::ofstream& out){ + out.open("test.txt"); + // This is the total number of events for the acceptance correction + auto total = results->intensity(); + double totalVal = total.first; + double totalErr = total.second; + + out << "################################################\n"; + out << "# All Values are acceptance corrected #\n"; + out << "################################################\n"; + out << "Total Events = " << totalVal << " +- " << totalErr << "\n"; + + // Coherent sums + for(size_t iAmp = 0; iAmp < reg.nUnqAmps; ++iAmp){ + const auto& wave = reg.uniqueAmps[iAmp]; + + // Coherent sums for unique amp based in refl + for(auto refl : {ReflSum::Positive, ReflSum::Negative, ReflSum::Neither}){ + auto key = std::make_pair(iAmp, refl); + if(reg.relations.ampReflToFull.count(key) == 0) + continue; + const auto& fullIndices = reg.relations.ampReflToFull.at(key); + if(fullIndices.empty()) + continue; + // Build full name list + std::vector fullNames; + for(auto idx : fullIndices) + fullNames.push_back(reg.ampsList[idx].fullName); - for (int i = 0; i < nAmps; i++) { - if (ampsumPosRefl[i].empty()) continue; - outfile << "FIT FRACTION (coherent sum) PosRefl " << amphistname[i] - << " = " - << results.intensity(ampsumPosRefl[i]).first / - results.intensity().first - << " +- " - << results.intensity(ampsumPosRefl[i]).second / - results.intensity().first - << endl; - outfile << "FIT FRACTION (coherent sum) NegRefl " << amphistname[i] - << " = " - << results.intensity(ampsumNegRefl[i]).first / - results.intensity().first - << " +- " - << results.intensity(ampsumNegRefl[i]).second / - results.intensity().first - << endl; + auto val = results->intensity(fullNames); + + out << "FIT FRACTION (coherent sum) " + << reflToString(refl) << " " + << wave.ampName << " = " + << val.first / totalVal << " +- " + << val.second / totalVal << "\n"; } + } + + // Coherent sums for Jl families + for(const auto& [key, fullIndices] : reg.relations.jlReflToFull){ + auto [jlFamily, refl] = key; + if(fullIndices.empty()) continue; + + std::vector fullNames; + for(auto idx : fullIndices) + fullNames.push_back(reg.ampsList[idx].fullName); + + auto val = results->intensity(fullNames); + + out << "FIT FRACTION (coherent sum) " + << reflToString(refl) << " " + << jlFamily << " = " + << val.first / totalVal << " +- " + << val.second / totalVal << "\n"; + } + + // Individual amplitudes + for(const auto& ampKey : reg.ampsList){ + + std::vector singleAmp = { ampKey.fullName }; + + auto val = results->intensity(singleAmp); + + out << "FIT FRACTION " << ampKey.fullName << " = " + << val.first / totalVal << " +- " + << val.second / totalVal << "\n"; + } + + out << "################################################\n"; + out << "################################################\n"; + out << "# Phase Differences #\n"; + out << "################################################\n"; + out << "################################################\n"; + + // In general, the phase differences are the same for all + // reactions. This is because the current use case is that + // multiple "reactions" are actually the same physical reaction + // but for a different data set like different run periods or + // beam pol. angle + const auto& firstReaction = results->reactionList()[0]; + + for(size_t i = 0; i < reg.ampsList.size(); ++i){ + const auto& amp1 = reg.ampsList[i]; + + // Only amplitudes in the first reaction to avoid duplicates + if(amp1.reactionName != firstReaction) continue; + if(amp1.reflSum == ReflSum::Neither) continue; + + for(size_t j = i+1; j < reg.ampsList.size(); ++j){ + const auto& amp2 = reg.ampsList[j]; - cout << "Computing phase differences" << endl; - for (unsigned int i = 0; i < phaseDiffNames.size(); i++) { - pair phaseDiff = results.phaseDiff( - phaseDiffNames[i].first, phaseDiffNames[i].second); - outfile << "PHASE DIFF " << phaseDiffNames[i].first << " " - << phaseDiffNames[i].second << " " << phaseDiff.first << " " - << phaseDiff.second << endl; + if(amp2.reactionName != firstReaction) continue; + if(amp2.reflSum == ReflSum::Neither) continue; + + // Only compare amplitudes in the same coherent sum + if(amp1.sumIndex != amp2.sumIndex) continue; + + auto phaseDiff = results->phaseDiff(amp1.fullName, amp2.fullName); + + out << "PHASE DIFF " + << amp1.fullName << " " + << amp2.fullName << " = " + << phaseDiff.first << " +- " + << phaseDiff.second << "\n"; } } +} - // covariance matrix - vector > covMatrix; - covMatrix = results.errorMatrix(); +void writeAmpInfoCSV(const AmplitudeRegistry& reg, + std::unique_ptr& results, + std::ofstream& out){ + out.open("fit_results.csv"); - // ************************ - // start the GUI - // ************************ + // This is the total number of events for the acceptance correction + auto total = results->intensity(); + double totalVal = total.first; + double totalErr = total.second; + + out << "Type,Name,Value,Error\n"; + out << "Total,Acceptance Corrected Events," << totalVal << "," + << totalErr << "\n"; + + // Coherent sums + for(size_t iAmp = 0; iAmp < reg.nUnqAmps; ++iAmp){ + const auto& wave = reg.uniqueAmps[iAmp]; + + // Coherent sums for unique amp based in refl + for(auto refl : {ReflSum::Positive, ReflSum::Negative, ReflSum::Neither}){ + auto key = std::make_pair(iAmp, refl); + if(reg.relations.ampReflToFull.count(key) == 0) continue; + const auto& fullIndices = reg.relations.ampReflToFull.at(key); + if(fullIndices.empty()) continue; + + std::vector fullNames; + for(auto idx : fullIndices) + fullNames.push_back(reg.ampsList[idx].fullName); + + auto val = results->intensity(fullNames); + + out << "CoherentSum," + << reflToString(refl) << " " << wave.ampName << "," + << val.first / totalVal << "," + << val.second / totalVal << "\n"; + } + } + + // Coherent sums for Jl families + for(const auto& [key, fullIndices] : reg.relations.jlReflToFull){ + auto [jlFamily, refl] = key; + if(fullIndices.empty()) continue; - if (showGui) { - cout << ">> Plot generator ready, starting GUI..." << endl; + std::vector fullNames; + for(auto idx : fullIndices) + fullNames.push_back(reg.ampsList[idx].fullName); + + auto val = results->intensity(fullNames); + + out << "CoherentSum," + << reflToString(refl) << " " << jlFamily << "," + << val.first / totalVal << "," + << val.second / totalVal << "\n"; + } + + // Individual amplitudes + for(const auto& ampKey : reg.ampsList){ + std::vector singleAmp = { ampKey.fullName }; + auto val = results->intensity(singleAmp); + + out << "IndividualAmp," + << ampKey.fullName << "," + << val.first / totalVal << "," + << val.second / totalVal << "\n"; + } + + // Phase differences + const auto& firstReaction = results->reactionList()[0]; + + for(size_t i = 0; i < reg.ampsList.size(); ++i){ + const auto& amp1 = reg.ampsList[i]; + if(amp1.reactionName != firstReaction) continue; + if(amp1.reflSum == ReflSum::Neither) continue; + + for(size_t j = i+1; j < reg.ampsList.size(); ++j){ + const auto& amp2 = reg.ampsList[j]; + if(amp2.reactionName != firstReaction) continue; + if(amp2.reflSum == ReflSum::Neither) continue; + if(amp1.sumIndex != amp2.sumIndex) continue; + + auto phaseDiff = results->phaseDiff(amp1.fullName, amp2.fullName); + + out << "PhaseDiff," + << amp1.fullName << " vs " << amp2.fullName << "," + << phaseDiff.first << "," + << phaseDiff.second << "\n"; + } + } +} + +// ...oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... +// Main Functions +// ...oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... + +bool parseInputs(int argc, char* argv[], ProgramOptions& options){ + // CLI parsing + if (argc < 2) { + std::cerr << "Usage: \n"; + std::cerr << "\tvecps_plotter " + "-o \n"<< std::endl; + return false; + } + + for (int i = 1; i < argc; i++){ + std::string arg(argv[i]); + if (arg == "-f"){ + options.resultsName = std::string(argv[++i]); + } + else if (arg == "-g"){ + options.showGui = true; + } + else if (arg == "-o"){ + options.outName = std::string(argv[++i]); + } + else if (arg == "-mcGen"){ + options.mcGen = true; + } + else if (arg == "-bkgFile"){ + options.bkgFile = true; + } + else if (arg == "-sFile"){ + options.sFile = true; + } + else if (arg == "-skipPlots"){ + options.skipPlots = true; + } + else if (arg == "-txt"){ + options.txtOutput = true; + } + else if (arg == "-csv"){ + options.csvOutput = true; + } + else if (arg == "-h"){ + std::cout << "\n Usage for: " << argv[0] << "\n\n"; + std::cout << "\t -f \t input fit results file path\n"; + std::cout << "\t -o \t output file path\n"; + std::cout << "\t -g \t show GUI. Default off\n"; + std::cout << "\t -mcGen\t output genMC. Default off\n"; + std::cout << "\t -sFile\t use a single reaction file set. Default off\n"; + std::cout << "\t -skipPlots\t skip generating plots. Default off\n"; + std::cout << "\t -txt\t output text file. Default off\n"; + std::cout << "\t -csv\t output CSV file. Default off\n"; + std::cout << "\t -h\t show this help message\n" + << std::endl; + return false; + } + else { + std::cerr << "Error: Unknown option '" << arg << "'\n"; + std::cerr << "Run with -h for help." << std::endl; + return false; + } + } + + // report used options + std::cout << "Fit results file name = " << options.resultsName << "\n"; + if(!options.skipPlots) + std::cout << "Output file name = " << options.outName << "\n"; + std::cout << "Show GUI? = " + << (options.showGui ? "Yes" : "No") << "\n"; + std::cout << "MC Gen? = " + << (options.mcGen ? "Yes" : "No") << "\n"; + std::cout << "Single File? = " + << (options.sFile ? "Yes" : "No") << "\n"; + std::cout << "Text Output? = " + << (options.txtOutput ? "Yes" : "No") << "\n"; + std::cout << "CSV Output? = " + << (options.csvOutput ? "Yes" : "No") << "\n" + << std::endl; + + return true; +} + +bool initializeAmpTools(ProgramOptions& options, + std::unique_ptr& results, + std::unique_ptr& plotResults, + AmplitudeRegistry& registry){ + // load the results and display the configuration info + std::cout << "Initializing AmpTools interface...\n"; + atiSetup(); + std::cout << "Interface initialized\n"; + std::cout << "Loading Fit results\n"; + results = std::make_unique(options.resultsName); + if (!results->valid()){ + std::cerr << "Invalid fit results in file: " + << options.resultsName << std::endl; + return false; + } + std::cout << "Fit results loaded\n"; + // set up the plot generator + if (options.mcGen){ + std::cout << "Plot results (MC gen on) \n"; + plotResults = + std::make_unique(*results, PlotGenerator::kDefault); + } + else { + std::cout << "Plot results (MC gen off) \n"; + plotResults = + std::make_unique(*results, PlotGenerator::kNoGenMC); + } + std::cout << "Initialized ati and PlotGen" << std::endl; + + // Build amplitude registry (once per file) + registry = buildRegistry(plotResults); + + return true; +} + +void analyzeFitFile(ProgramOptions& options, + std::unique_ptr& results, + std::unique_ptr& plotResults, + AmplitudeRegistry& registry){ + if(options.skipPlots){ + std::cout << "Skipping plot generation.\n"; + return; + } + cout << "\n*** Viewing Results Using AmpPlotter***\n\n"; + // Loop over reactions + // User will have to add the different root file outputs manually + auto reactionList = results->reactionList(); + for (size_t fileSet = 0; fileSet < reactionList.size(); ++fileSet){ + if (options.sFile && fileSet > 0) break; + + std::string reactionName = reactionList[fileSet]; + std::cout << "Processing reaction " << reactionName << " (" + << fileSet + 1 << "/" << reactionList.size() << ")\n"; + + if (!options.sFile){ + options.outName = reactionName + ".root"; + std::cout << "Output file: " << options.outName << std::endl; + } + + TFile* plotFile = new TFile(options.outName.c_str(), "RECREATE"); + TH1::AddDirectory(kFALSE); // avoid automatic ownership for GUI use + + plotResults->enableReaction(reactionName); + // Generate plot configurations + auto plotConfigs = generatePlotConfigs(registry); + // Loop over plot configurations + bool first = true; + bool singleData = false; + for (const auto& cfg : plotConfigs){ + + if (!first){ + applyConfig(cfg, registry, plotResults); + } + else{ + // On first iteration, just enable the relevant amps/sums + for (auto i : cfg.ampIndices) plotResults->enableAmp(i); + for (auto i : cfg.sumIndices) plotResults->enableSum(i); + first = false; + } + + PlotLabels cfgLabel = makeLabel(cfg, registry); + generatePlots(options, plotResults, reactionName, + *plotFile, singleData,cfg, cfgLabel); + } + plotFile->Close(); + } +} + +void startGUI(ProgramOptions& options, + std::unique_ptr& results, + std::unique_ptr& plotResults){ + if(options.skipPlots){ + std::cout << "No plots = No GUI\n"; + return; + } + string reactionName = results->reactionList()[0]; + plotResults->enableReaction(reactionName); + // start the GUI + if (options.showGui){ + std::cout << ">> Plot generator ready, starting GUI..." << std::endl; int dummy_argc = 0; char* dummy_argv[] = {}; @@ -424,15 +1087,44 @@ int main(int argc, char* argv[]) { gStyle->SetFrameFillColor(10); gStyle->SetFrameFillStyle(1001); - cout << " Initialized App " << endl; - PlotFactory factory(plotGen); - cout << " Created Plot Factory " << endl; + std::cout << " Initialized App " << std::endl; + PlotFactory factory(*plotResults); + std::cout << " Created Plot Factory " << std::endl; PlotterMainWindow mainFrame(gClient->GetRoot(), factory); - cout << " Main frame created " << endl; + std::cout << " Main frame created " << std::endl; app.Run(); - cout << " App running" << endl; + std::cout << " App running" << std::endl; } - return 0; } + +int main(int argc, char* argv[]){ + + ProgramOptions options; + if (!parseInputs(argc, argv, options)){ + std::cout << "Error parsing inputs" << std::endl; + return 1; // Parsing failed + } + + std::unique_ptr results; + std::unique_ptr plotResults; + AmplitudeRegistry registry; + if (!initializeAmpTools(options, results, plotResults, registry)){ + std::cout << "Error initializing AmpTools" << std::endl; + return 1; // Parsing failed + } + + analyzeFitFile(options, results, plotResults,registry); + ofstream parsOutputFile; + writeAmpInfo(registry,results,parsOutputFile); + parsOutputFile.close(); + + ofstream csvOutputFile; + writeAmpInfoCSV(registry,results,csvOutputFile); + csvOutputFile.close(); + + startGUI(options, results, plotResults); + + return 0; +} \ No newline at end of file From 7f5f9ec26adefe7f9a3ee59a66ebf68e80ffa903 Mon Sep 17 00:00:00 2001 From: edbarriga Date: Mon, 30 Mar 2026 15:05:43 -0400 Subject: [PATCH 2/2] Now VecPsPlotGenerator will give the names to vecpsplotter, and parameters will now be written in txt-based outputs --- .../AMPTOOLS_DATAIO/VecPsPlotGenerator.cc | 2 + .../AMPTOOLS_DATAIO/VecPsPlotGenerator.h | 27 ++- .../vecps_plotter/vecps_plotter.cc | 174 +++++++++--------- 3 files changed, 117 insertions(+), 86 deletions(-) diff --git a/src/libraries/AMPTOOLS_DATAIO/VecPsPlotGenerator.cc b/src/libraries/AMPTOOLS_DATAIO/VecPsPlotGenerator.cc index 7bfac7e70..15ab13dd5 100644 --- a/src/libraries/AMPTOOLS_DATAIO/VecPsPlotGenerator.cc +++ b/src/libraries/AMPTOOLS_DATAIO/VecPsPlotGenerator.cc @@ -59,6 +59,8 @@ void VecPsPlotGenerator::createHistograms( ) { bookHistogram(kPhiH, new Histogram1D(50, -PI, PI, "Phi_H", ";#phi_H [rad.]")); bookHistogram(kProd_Ang, new Histogram1D(50, -PI, PI, "Prod_Ang", ";Prod_Ang (#Phi) [rad.]")); + bookHistogram(kProdOffset, new Histogram1D(50, -PI, PI, "ProdOffset", + ";Prod_Ang (#Phi) Uncorrected [rad.]" ) ); bookHistogram(kt, new Histogram1D(100, 0, 2.0, "t", ";-t")); bookHistogram(kRecoilMass, new Histogram1D(100, 0.9, 1.9, "MRecoil", ";Invariant Mass of Recoil [GeV]")); diff --git a/src/libraries/AMPTOOLS_DATAIO/VecPsPlotGenerator.h b/src/libraries/AMPTOOLS_DATAIO/VecPsPlotGenerator.h index 8ad988655..4f54c2a5b 100644 --- a/src/libraries/AMPTOOLS_DATAIO/VecPsPlotGenerator.h +++ b/src/libraries/AMPTOOLS_DATAIO/VecPsPlotGenerator.h @@ -17,14 +17,33 @@ class VecPsPlotGenerator : public PlotGenerator public: // create an index for different histograms - enum { kVecPsMass = 0, kCosTheta = 1, kPhi = 2, kCosThetaH = 3, kPhiH = 4, - kProd_Ang = 5, kt = 6, kRecoilMass = 7, kProtonPsMass = 8, - kRecoilPsMass = 9, kLambda = 10, kDalitz = 11, kPhi_ProdVsPhi = 12, - kPhiOffsetVsPhi = 13, kNumHists}; + enum Hist_index{ kVecPsMass, kCosTheta, kPhi, kCosThetaH, kPhiH, + kProd_Ang, kProdOffset, kt, kRecoilMass, kProtonPsMass, + kRecoilPsMass, kLambda, kDalitz, kPhi_ProdVsPhi, + kPhiOffsetVsPhi, kNumHists}; VecPsPlotGenerator( const FitResults& results, Option opt); VecPsPlotGenerator( const FitResults& results ); VecPsPlotGenerator( ); + + static std::string numToString(Hist_index kHistName){ + switch(kHistName){ + case VecPsPlotGenerator::kVecPsMass: return "MVecPs"; + case VecPsPlotGenerator::kCosTheta: return "CosTheta"; + case VecPsPlotGenerator::kPhi: return "Phi"; + case VecPsPlotGenerator::kCosThetaH: return "CosTheta_H"; + case VecPsPlotGenerator::kPhiH: return "Phi_H"; + case VecPsPlotGenerator::kProd_Ang: return "Prod_Ang"; + case VecPsPlotGenerator::kt: return "t"; + case VecPsPlotGenerator::kRecoilPsMass: return "MRecoilPs"; + case VecPsPlotGenerator::kLambda: return "Lambda"; + case VecPsPlotGenerator::kProdOffset: return "ProdOffset"; + case VecPsPlotGenerator::kPhi_ProdVsPhi: return "Phi_ProdVsPhi"; + case VecPsPlotGenerator::kPhiOffsetVsPhi: return "PhiOffsetVsPhi"; + // Add more variables here if needed + default: return "Unknown"; + } +} private: diff --git a/src/programs/AmplitudeAnalysis/vecps_plotter/vecps_plotter.cc b/src/programs/AmplitudeAnalysis/vecps_plotter/vecps_plotter.cc index b10dec0d7..8d1c44e99 100644 --- a/src/programs/AmplitudeAnalysis/vecps_plotter/vecps_plotter.cc +++ b/src/programs/AmplitudeAnalysis/vecps_plotter/vecps_plotter.cc @@ -647,35 +647,19 @@ void generatePlots(ProgramOptions& options, std::unique_ptr& plot } // Loop over all histograms defined in VecPsPlotGenerator - for (unsigned int ivar = 0; ivar < VecPsPlotGenerator::kNumHists; ++ivar){ + // The loop variable should be defined using the first entry in the + // enum, in this case kVecPsMass (change if needed). The last one is + // kNumHists and this one should not change + for (VecPsPlotGenerator::Hist_index ivar = VecPsPlotGenerator::kVecPsMass; + ivar < VecPsPlotGenerator::kNumHists; + ivar = static_cast(ivar + 1)){ Histogram* hist = plotResults->projection(ivar, reactionName, iType); if (!hist) continue; // safety skip TH1* histRoot = hist->toRoot(); - std::string histName = ""; - // Assign proper names to each histogram. Note that the names are - // defined in VecPsPlotGenerator - // TODO: add a parser for the variables similar to reflToString to - // VecPsPlotGenerator to make this more automatic and - // avoid the need to edit this file when adding new variables - switch (ivar) { - case VecPsPlotGenerator::kVecPsMass: histName = "MVecPs"; break; - case VecPsPlotGenerator::kCosTheta: histName = "CosTheta"; break; - case VecPsPlotGenerator::kPhi: histName = "Phi"; break; - case VecPsPlotGenerator::kCosThetaH: histName = "CosTheta_H"; break; - case VecPsPlotGenerator::kPhiH: histName = "Phi_H"; break; - case VecPsPlotGenerator::kProd_Ang: histName = "Prod_Ang"; break; - case VecPsPlotGenerator::kt: histName = "t"; break; - case VecPsPlotGenerator::kRecoilMass: histName = "MRecoil"; break; - case VecPsPlotGenerator::kProtonPsMass: histName = "MProtonPs"; break; - case VecPsPlotGenerator::kRecoilPsMass: histName = "MRecoilPs"; break; - case VecPsPlotGenerator::kLambda: histName = "Lambda"; break; - case VecPsPlotGenerator::kDalitz: histName = "Dalitz"; break; - case VecPsPlotGenerator::kPhi_ProdVsPhi: histName = "Phi_ProdVsPhi"; break; - case VecPsPlotGenerator::kPhiOffsetVsPhi: histName = "PhiOffsetVsPhi"; break; - // Add more variables here if needed - default: continue; - } + // Assign proper names to each histogram. Note that the names + // and the translation table are defined in VecPsPlotGenerator + std::string histName = VecPsPlotGenerator::numToString(ivar); if(!(iType==PlotGenerator::kData)) histName += "_" + dataType + cfgLabel.histName; @@ -690,20 +674,19 @@ void generatePlots(ProgramOptions& options, std::unique_ptr& plot // and phase differences between amplitudes. It has some patterns to help // parsing, but csv file might be preferred for that purpose. // Note: the default behavior is the acceptance corrected intensities -// TODO: include things like d/s ratios and other parameters void writeAmpInfo(const AmplitudeRegistry& reg, - std::unique_ptr& results, - std::ofstream& out){ - out.open("test.txt"); + std::unique_ptr& results){ + + std::ofstream parsOutputFile("params.txt"); // This is the total number of events for the acceptance correction auto total = results->intensity(); double totalVal = total.first; double totalErr = total.second; - out << "################################################\n"; - out << "# All Values are acceptance corrected #\n"; - out << "################################################\n"; - out << "Total Events = " << totalVal << " +- " << totalErr << "\n"; + parsOutputFile << "################################################\n"; + parsOutputFile << "# All Values are acceptance corrected #\n"; + parsOutputFile << "################################################\n"; + parsOutputFile << "Total Events = " << totalVal << " +- " << totalErr << "\n"; // Coherent sums for(size_t iAmp = 0; iAmp < reg.nUnqAmps; ++iAmp){ @@ -724,11 +707,11 @@ void writeAmpInfo(const AmplitudeRegistry& reg, auto val = results->intensity(fullNames); - out << "FIT FRACTION (coherent sum) " - << reflToString(refl) << " " - << wave.ampName << " = " - << val.first / totalVal << " +- " - << val.second / totalVal << "\n"; + parsOutputFile << "FIT FRACTION (coherent sum) " + << reflToString(refl) << " " + << wave.ampName << " = " + << val.first / totalVal << " +- " + << val.second / totalVal << "\n"; } } @@ -743,11 +726,11 @@ void writeAmpInfo(const AmplitudeRegistry& reg, auto val = results->intensity(fullNames); - out << "FIT FRACTION (coherent sum) " - << reflToString(refl) << " " - << jlFamily << " = " - << val.first / totalVal << " +- " - << val.second / totalVal << "\n"; + parsOutputFile << "FIT FRACTION (coherent sum) " + << reflToString(refl) << " " + << jlFamily << " = " + << val.first / totalVal << " +- " + << val.second / totalVal << "\n"; } // Individual amplitudes @@ -757,16 +740,16 @@ void writeAmpInfo(const AmplitudeRegistry& reg, auto val = results->intensity(singleAmp); - out << "FIT FRACTION " << ampKey.fullName << " = " - << val.first / totalVal << " +- " - << val.second / totalVal << "\n"; + parsOutputFile << "FIT FRACTION " << ampKey.fullName << " = " + << val.first / totalVal << " +- " + << val.second / totalVal << "\n"; } - out << "################################################\n"; - out << "################################################\n"; - out << "# Phase Differences #\n"; - out << "################################################\n"; - out << "################################################\n"; + parsOutputFile << "################################################\n"; + parsOutputFile << "################################################\n"; + parsOutputFile << "# Phase Differences #\n"; + parsOutputFile << "################################################\n"; + parsOutputFile << "################################################\n"; // In general, the phase differences are the same for all // reactions. This is because the current use case is that @@ -793,27 +776,43 @@ void writeAmpInfo(const AmplitudeRegistry& reg, auto phaseDiff = results->phaseDiff(amp1.fullName, amp2.fullName); - out << "PHASE DIFF " - << amp1.fullName << " " - << amp2.fullName << " = " - << phaseDiff.first << " +- " - << phaseDiff.second << "\n"; + parsOutputFile << "PHASE DIFF " + << amp1.fullName << " vs " + << amp2.fullName << " = " + << phaseDiff.first << " +- " + << phaseDiff.second << "\n"; } } + + parsOutputFile << "################################################\n"; + parsOutputFile << "################################################\n"; + parsOutputFile << "# Amp Parameters #\n"; + parsOutputFile << "################################################\n"; + parsOutputFile << "################################################\n"; + + // currently plots all parameters, including production amplitudes + // a helper function could be added to filter the parameters of interest + std::vector parNameList = results->parNameList(); + for(const auto& parName : parNameList){ + double parVal = results->parValue(parName); + double parErr = results->parError(parName); + parsOutputFile << parName << ": " + << parVal << " +- " + << parErr << "\n"; + } + parsOutputFile.close(); } void writeAmpInfoCSV(const AmplitudeRegistry& reg, - std::unique_ptr& results, - std::ofstream& out){ - out.open("fit_results.csv"); - + std::unique_ptr& results){ + std::ofstream csvOutputFile("params.csv"); // This is the total number of events for the acceptance correction auto total = results->intensity(); double totalVal = total.first; double totalErr = total.second; - out << "Type,Name,Value,Error\n"; - out << "Total,Acceptance Corrected Events," << totalVal << "," + csvOutputFile << "Type,Name,Value,Error\n"; + csvOutputFile << "Total,Acceptance Corrected Events," << totalVal << "," << totalErr << "\n"; // Coherent sums @@ -833,7 +832,7 @@ void writeAmpInfoCSV(const AmplitudeRegistry& reg, auto val = results->intensity(fullNames); - out << "CoherentSum," + csvOutputFile << "CoherentSum," << reflToString(refl) << " " << wave.ampName << "," << val.first / totalVal << "," << val.second / totalVal << "\n"; @@ -851,7 +850,7 @@ void writeAmpInfoCSV(const AmplitudeRegistry& reg, auto val = results->intensity(fullNames); - out << "CoherentSum," + csvOutputFile << "CoherentSum," << reflToString(refl) << " " << jlFamily << "," << val.first / totalVal << "," << val.second / totalVal << "\n"; @@ -862,7 +861,7 @@ void writeAmpInfoCSV(const AmplitudeRegistry& reg, std::vector singleAmp = { ampKey.fullName }; auto val = results->intensity(singleAmp); - out << "IndividualAmp," + csvOutputFile << "IndividualAmp," << ampKey.fullName << "," << val.first / totalVal << "," << val.second / totalVal << "\n"; @@ -884,12 +883,25 @@ void writeAmpInfoCSV(const AmplitudeRegistry& reg, auto phaseDiff = results->phaseDiff(amp1.fullName, amp2.fullName); - out << "PhaseDiff," + csvOutputFile << "PhaseDiff," << amp1.fullName << " vs " << amp2.fullName << "," << phaseDiff.first << "," << phaseDiff.second << "\n"; } } + + // currently plots all parameters, including production amplitudes + // a helper function could be added to filter the parameters of interest + std::vector parNameList = results->parNameList(); + for(const auto& parName : parNameList){ + double parVal = results->parValue(parName); + double parErr = results->parError(parName); + csvOutputFile << "Parameter," + << parName << "," + << parVal << "," + << parErr << "\n"; + } + csvOutputFile.close(); } // ...oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... @@ -1063,14 +1075,12 @@ void analyzeFitFile(ProgramOptions& options, } void startGUI(ProgramOptions& options, - std::unique_ptr& results, std::unique_ptr& plotResults){ if(options.skipPlots){ std::cout << "No plots = No GUI\n"; return; } - string reactionName = results->reactionList()[0]; - plotResults->enableReaction(reactionName); + auto plotResults_local = std::move(plotResults); // start the GUI if (options.showGui){ std::cout << ">> Plot generator ready, starting GUI..." << std::endl; @@ -1078,6 +1088,7 @@ void startGUI(ProgramOptions& options, int dummy_argc = 0; char* dummy_argv[] = {}; TApplication app("app", &dummy_argc, dummy_argv); + app.SetReturnFromRun(true); gStyle->SetFillColor(10); gStyle->SetCanvasColor(10); @@ -1088,13 +1099,14 @@ void startGUI(ProgramOptions& options, gStyle->SetFrameFillStyle(1001); std::cout << " Initialized App " << std::endl; - PlotFactory factory(*plotResults); + auto factory = new PlotFactory(*plotResults_local); std::cout << " Created Plot Factory " << std::endl; - PlotterMainWindow mainFrame(gClient->GetRoot(), factory); + auto mainFrame = new PlotterMainWindow(gClient->GetRoot(), *factory); std::cout << " Main frame created " << std::endl; - - app.Run(); + mainFrame->Connect("CloseWindow()", "TApplication", gApplication, "Terminate()"); std::cout << " App running" << std::endl; + app.Run(); + std::cout << "GUI closed cleanly\n"; } } @@ -1115,16 +1127,14 @@ int main(int argc, char* argv[]){ return 1; // Parsing failed } - analyzeFitFile(options, results, plotResults,registry); - ofstream parsOutputFile; - writeAmpInfo(registry,results,parsOutputFile); - parsOutputFile.close(); - - ofstream csvOutputFile; - writeAmpInfoCSV(registry,results,csvOutputFile); - csvOutputFile.close(); + analyzeFitFile(options, results, plotResults, registry); + + writeAmpInfo(registry,results); + + writeAmpInfoCSV(registry,results); - startGUI(options, results, plotResults); + // GUI needs to be started last + startGUI(options, plotResults); return 0; } \ No newline at end of file