From 582ae8161a3750a76cff0fed6cec61263d2fffde Mon Sep 17 00:00:00 2001 From: Pawel Plesniak Date: Sat, 24 Jan 2026 17:29:45 +0100 Subject: [PATCH 1/2] New tools --- .../MCTruthToDetectorResponseDifference.C | 120 +++++ DataSummaries/SignalFractionTable.C | 121 +++++ Plotting/plotMWDResults_noADCClock.C | 442 ++++++++++++++++++ 3 files changed, 683 insertions(+) create mode 100644 DataSummaries/MCTruthToDetectorResponseDifference.C create mode 100644 DataSummaries/SignalFractionTable.C create mode 100644 Plotting/plotMWDResults_noADCClock.C diff --git a/DataSummaries/MCTruthToDetectorResponseDifference.C b/DataSummaries/MCTruthToDetectorResponseDifference.C new file mode 100644 index 0000000..6a0d1ec --- /dev/null +++ b/DataSummaries/MCTruthToDetectorResponseDifference.C @@ -0,0 +1,120 @@ +// Generates the table of ratio of MC truth normalization to detector response normalization differences +// Usage example - $ root -l -q 'MCTruthToDetectorResponseDifference.C({{1.679e10, 3.02e9}, {0.0, 0.0}, {1.583e11, 2.28e10}}, {{4.275e11, 8.56e10}, {0.0, 0.0}, {7.660e13, 1.180e13}})' +// Original author - Pawel Plesniak + + +void customErrorHandler(int level, Bool_t abort, const char* location, const char* message) { + /* + Description + Define a custom error handler that won't print the stack trace but will print an error message and exit. + */ + std::cerr << message << std::endl; + if (level > kInfo) + exit(1); +}; + +std::string doubleToStringScientific(double value, int sigFigs) { + /* + Description + Converts a double to a string with a defined number of significant figures + + Arguments + value - value to convert to string + sigFigs - number of significant figures to use + + Variables + out - output string stream used to convert 'value' to a std::string + */ + std::ostringstream out; + out << std::scientific << std::setprecision(sigFigs - 1) << value; + return out.str(); +}; + +std::string doubleToStringFixed(double value, int sigFigs) { + /* + Description + Converts a double to a string with a defined number of significant figures + + Arguments + value - value to convert to string + sigFigs - number of significant figures to use + + Variables + out - output string stream used to convert 'value' to a std::string + */ + std::ostringstream out; + out << std::fixed << std::setprecision(sigFigs - 1) << value; + return out.str(); +}; + +void MCTruthToDetectorResponseDifference(std::vector> mcTruthNormalizations, std::vector> detectorResponseNormalizations, double nPOTs = 2.69e16) { + /* + Description + Generates the table describing the ratio of MC truth normalization to detector response normalization differences + + Arguments + mcTruthNormalizations - vector of vectors containing the MC truth normalizations and associated uncertainties for 347, 844, and 1809 keV + detectorResponseNormalizations - vector of vectors containing the detector response normalizations and associated uncertainties for 347, 844, and 1809 keV + + Variables + nSF - number of significant figures to use in the table + correctionNameColumnWidth - width of the correction name column + signal347ColumnWidth - width of the 347 keV signal column + signal844ColumnWidth - width of the 844 keV signal column + signal1809ColumnWidth - width of the 1809 keV signal column + fullWidth - total width of the table + i, j - iterators + ratioMCToDetector - vector of vectors containing the ratio of MC truth normalization to detector response normalization and associated uncertainties for 347, 844, and 1809 keV + */ + + // Sanity check + if (mcTruthNormalizations.size() != 3 || detectorResponseNormalizations.size() != 3) + Fatal("MCTruthToDetectorResponseDifference", "Incorrect format of input normalizations, expected as {{value347, uncertainty347}, {value844, uncertainty844}, {value1809, uncertainty1809}}"); + if (mcTruthNormalizations == detectorResponseNormalizations) + Fatal("MCTruthToDetectorResponseDifference", "Input normalizations are identical, no difference to calculate"); + std::vector> zeroVector = {{0, 0}, {0, 0}, {0, 0}}; + if (mcTruthNormalizations == zeroVector || detectorResponseNormalizations == zeroVector) + Fatal("MCTruthToDetectorResponseDifference", "Input normalizations cannot be zero"); + + // Calculate the ratio of MC truth normalization to detector response normalization and its errors + std::vector> ratioMCToDetector = {{0, 0}, {0, 0}, {0, 0}}; + for (int i = 0; i < 3; i++) { + ratioMCToDetector[i][0] = detectorResponseNormalizations[i][0] / mcTruthNormalizations[i][0]; + ratioMCToDetector[i][1] = ratioMCToDetector[i][0] * std::sqrt(std::pow(mcTruthNormalizations[i][1]/mcTruthNormalizations[i][0], 2) + std::pow(detectorResponseNormalizations[i][1]/detectorResponseNormalizations[i][0], 2)); + }; + + + // Define the signal photon energy strings + std::string e347 = "347 keV"; + std::string e844 = "844 keV"; + std::string e1809 = "1809 keV"; + + // Define table formatting + const int nSF = 4; + const int correctionNameColumnWidth = 20, signalColumnWidth = 30; + const int fullWidth = correctionNameColumnWidth + signalColumnWidth * 3; + + // Print the title line and rules + std::string title = "Normalization ratio (MC truth / Detector response)"; + std::cout << std::string(fullWidth, '=') << std::endl; // Title line + std::cout << std::string((fullWidth - static_cast(title.size())) / 2, ' ') << title << std::endl; + std::cout << std::string(fullWidth, '=') << std::endl; // Title line + std::cout << std::setw(correctionNameColumnWidth) << std::left << "" << std::setw(signalColumnWidth) << std::left << "347 keV" << std::setw(signalColumnWidth) << std::left << "844 keV" << std::setw(signalColumnWidth) << std::left << "1809 keV" << std::endl; + std::cout << std::string(fullWidth, '-') << std::endl; // Section line + std::cout << std::setw(correctionNameColumnWidth) << std::left << "MC Truth";; + for (int j = 0; j < 3; j++) + std::cout << std::setw(signalColumnWidth) << std::left << doubleToStringScientific(mcTruthNormalizations[j][0], nSF) + " ± " + doubleToStringScientific(mcTruthNormalizations[j][1], nSF); + std::cout << std::endl; + std::cout << std::setw(correctionNameColumnWidth) << std::left << "Detector Response";; + for (int j = 0; j < 3; j++) + std::cout << std::setw(signalColumnWidth) << std::left << doubleToStringScientific(detectorResponseNormalizations[j][0], nSF) + " ± " + doubleToStringScientific(detectorResponseNormalizations[j][1], nSF); + std::cout << std::endl; + std::cout << std::string(fullWidth, '-') << std::endl; // Section line + std::cout << std::setw(correctionNameColumnWidth) << std::left << "Ratio";; + for (int j = 0; j < 3; j++) + std::cout << std::setw(signalColumnWidth) << std::left << doubleToStringFixed(ratioMCToDetector[j][0], nSF) + " ± " + doubleToStringFixed(ratioMCToDetector[j][1], nSF); + std::cout << std::endl; + std::cout << std::string(fullWidth, '=') << std::endl; // End line + std::cout << std::endl; // Buffer line + return; +}; \ No newline at end of file diff --git a/DataSummaries/SignalFractionTable.C b/DataSummaries/SignalFractionTable.C new file mode 100644 index 0000000..f286510 --- /dev/null +++ b/DataSummaries/SignalFractionTable.C @@ -0,0 +1,121 @@ +// Generates the table of signal photons from the total measured signal spectrum +// Usage example - $ root -l -q 'SignalFractionTable.C(2082. 0. 21777)' +// Original author - Pawel Plesniak + +void customErrorHandler(int level, Bool_t abort, const char* location, const char* message) { + /* + Description + Define a custom error handler that won't print the stack trace but will print an error message and exit. + */ + std::cerr << message << std::endl; + if (level > kInfo) + exit(1); +}; + +std::string doubleToString(double value, int sigFigs) { + /* + Description + Converts a double to a string with a defined number of significant figures + + Arguments + value - value to convert to string + sigFigs - number of significant figures to use + + Variables + out - output string stream used to convert 'value' to a std::string + */ + std::ostringstream out; + out << std::fixed << std::setprecision(sigFigs - 1) << value; + return out.str(); +}; + +void SignalFractionTable(double n347, double n844, double n1809) { + /* + Description + Generates the table describing the number of signal photons measured in the STM energy range + + Arguments + n347 - total number of measured photons 347 keV STM energy range (use value from plot generated by plotMWDResults.C) + n844 - total number of measured photons 844 keV STM energy range (use value from plot generated by plotMWDResults.C) + n1809 - total number of measured photons 1809 keV STM energy range (use value from plot generated by plotMWDResults.C) + + Variables + */ + + // Determine the total signal count uncertainties + const double u347 = std::sqrt(n347); + const double u844 = std::sqrt(n844); + const double u1809 = std::sqrt(n1809); + + // Determine the normalized total signal count uncertainties + const double uf347 = u347 / n347; + const double uf844 = u844 / n844; + const double uf1809 = u1809 / n1809; + + + // Define the signal to total ratios + const double st347 = 0.2712; + const double st844 = 0.6987; + const double st1809 = 0.6032; + + // Define the signal to total ratio uncertainties + const double ust347 = 0.0541; + const double ust844 = 0.2064; + const double ust1809 = 0.0715; + + // Determine the normalized signal to total ratio uncertainties + const double ufst347 = ust347 / st347; + const double ufst844 = ust844 / st844; + const double ufst1809 = ust1809 / st1809; + + + // Determine the signal photon counts + const double nSignal347 = n347 * st347; + const double nSignal844 = n844 * st844; + const double nSignal1809 = n1809 * st1809; + + // Determine the signal photon count uncertainties + const double uSignal347 = nSignal347 * std::sqrt(std::pow(uf347, 2) + std::pow(ufst347, 2));; + const double uSignal844 = nSignal844 * std::sqrt(std::pow(uf844, 2) + std::pow(ufst844, 2));; + const double uSignal1809 = nSignal1809 * std::sqrt(std::pow(uf1809, 2) + std::pow(ufst1809, 2));; + + + // Define the signal photon energy strings + std::string e347 = "347 keV"; + std::string e844 = "844 keV"; + std::string e1809 = "1809 keV"; + + // Define measured photon strings + std::string t347 = doubleToString(n347, 1) + " ± " + doubleToString(u347, 1); + std::string t844 = doubleToString(n844, 1) + " ± " + doubleToString(u844, 1); + std::string t1809 = doubleToString(n1809, 1) + " ± " + doubleToString(u1809, 1); + + // Define signal to total ratio strings + std::string ratio347 = doubleToString(st347, 5) + " ± " + doubleToString(ust347, 5); + std::string ratio844 = doubleToString(st844, 5) + " ± " + doubleToString(ust844, 5); + std::string ratio1809 = doubleToString(st1809, 5) + " ± " + doubleToString(ust1809, 5); + + // Define signal photon strings + std::string sig347 = doubleToString(nSignal347, 1) + " ± " + doubleToString(uSignal347, 1); + std::string sig844 = doubleToString(nSignal844, 1) + " ± " + doubleToString(uSignal844, 1); + std::string sig1809 = doubleToString(nSignal1809, 1) + " ± " + doubleToString(uSignal1809, 1); + + + // Define column widths + const int titleWidth = 40, signalWidth = 30, fullWidth = titleWidth + signalWidth * 3; + + // Print the table + std::string tableTitle = "Signal Photons Summary Table"; + std::cout << std::string(fullWidth, '=') << std::endl; + std::cout << std::string((fullWidth - tableTitle.size())/2, ' ') << tableTitle << std::endl; + std::cout << std::string(fullWidth, '=') << std::endl; + std::cout << std::setw(titleWidth) << std::left << " " << std::setw(signalWidth) << std::left << e347 << std::setw(signalWidth) << std::left << e844 << std::setw(signalWidth) << std::left << e1809 << std::endl; + std::cout << std::string(fullWidth, '-') << std::endl; + std::cout << std::setw(titleWidth) << std::left << "Total Measured Photons: " << std::setw(signalWidth) << std::left << t347 << std::setw(signalWidth) << std::left << t844 << std::setw(signalWidth) << std::left << t1809 << std::endl; + std::cout << std::setw(titleWidth) << std::left << "Signal to Total ratio: " << std::setw(signalWidth) << std::left << ratio347 << std::setw(signalWidth) << std::left << ratio844 << std::setw(signalWidth) << std::left << ratio1809 << std::endl; + std::cout << std::string(fullWidth, '-') << std::endl; + std::cout << std::setw(titleWidth) << std::left << "Signal Photons: " << std::setw(signalWidth) << std::left << sig347 << std::setw(signalWidth) << std::left << sig844 << std::setw(signalWidth) << std::left << sig1809 << std::endl; + std::cout << std::string(fullWidth, '=') << std::endl; + + return; +}; \ No newline at end of file diff --git a/Plotting/plotMWDResults_noADCClock.C b/Plotting/plotMWDResults_noADCClock.C new file mode 100644 index 0000000..a76730f --- /dev/null +++ b/Plotting/plotMWDResults_noADCClock.C @@ -0,0 +1,442 @@ +// Generates the analysis path of the MWD algorithm +// Usage example - $ root -l -q 'plotMWDResults.C("FinalData/DigNoZS.root", "MWDSpectra/ttree", 2000.0, 50, 0.1, 500, 5, 0.5, 0.8, 1.1)' +// Original author - Pawel Plesniak + +#include + +void customErrorHandler(int level, Bool_t abort, const char* location, const char* message) { + /* + Description + Define a custom error handler that won't print the stack trace but will print an error message and exit. + */ + std::cerr << message << std::endl; + if (level > kInfo) + exit(1); +}; + +void collectData(const std::string fileName, const std::string treeName, std::vector ×, std::vector &energies) { + /* + Description + Collects all the required data from virtual detector TTrees + + Arguments + fileName - as documented in function "plotMWDResults" + treeName - as documented in function "plotMWDResults" + times - as documented in function "plotMWDResults" + energies - as documented in function "plotMWDResults" + + Variables + file - ROOT TFile interface + branches - ROOT TFile interface to TTree branches + branchNames - vector of branch names + dataTime - time from file + dataE - energy from file + entries - number of entries in the TTree + */ + std::cout << "Processing file " << fileName << std::endl; + + // Get the branch + TFile *file = new TFile(fileName.c_str()); + if (!file || file->IsZombie()) { + Fatal("collectData", "Failed to open the file."); + }; + TTree *tree = (TTree*)file->Get(treeName.c_str()); + if (!tree) + Fatal("collectData", "Requested tree does not exist in the file."); + + // Get the list of branches to check if they exist + TObjArray *branches = tree->GetListOfBranches(); + std::vector branchNames; + for (int i = 0; i < branches->GetEntries(); ++i) { + TBranch *branch = dynamic_cast(branches->At(i)); + if (branch) + branchNames.push_back(branch->GetName()); + }; + + // If the branches exist, assign them to the appropriate variables + double dataE; + double dataTime; + if (std::find(branchNames.begin(), branchNames.end(), "time") != branchNames.end()) + tree->SetBranchAddress("time", &dataTime); + else + Fatal("collectData", "Requested branch 'time' does not exist in the file."); + if (std::find(branchNames.begin(), branchNames.end(), "E") != branchNames.end()) + tree->SetBranchAddress("E", &dataE); + else + Fatal("collectData", "Requested branch 'E' does not exist in the file."); + + // Get the number of entries + int entries = tree->GetEntries(); + + // Collect the data + for (int i = 0; i < entries; i++) { + // Get the corresponding entry + tree->GetEntry(i); + + // Collect the data + times.push_back(dataTime); + energies.push_back(dataE); + }; + // Check if data has been collected + if (times.size() == 0 || energies.size() == 0) + Fatal("collectData", "No data was collected from this file"); + if (times.size() != energies.size()) + Fatal("collectData", "Unequal number of data points were collected"); + + // Clear + file->Close(); + delete file; + std::cout << "Finished processing file " << fileName << "\n" << std::endl; + return; +}; + + +std::string convertBinWidthToStr(double binWidth) { + /* + Description + Converts the bin width as a double to a string and returns it + + Arguments + binWidth - bin width + + Variables + stream - input string stream to read the bin width + binWidthStr - bin width as a str + */ + // Set up the stream to read in the bin width + std::stringstream stream; + + // Set the precision and read it in. If the bin width is an integer, set 0 decimal places, otherwise set 3 + stream << std::fixed << std::setprecision(((binWidth - (int)binWidth) < std::numeric_limits::epsilon()) ? 0 : 3) << binWidth; + + // Convert the stream object to a string + std::string binWidthStr = stream.str(); + + // Clear + stream.str(""); + std::cout << std::defaultfloat; + return binWidthStr; +}; + +void plot(std::vector times, std::vector energies, std::vector binWidths, double ERed, double EMin, bool convertkeVToMeV, const double signalAcceptance, bool highResolution, bool timeCuts, bool flashCut347) { + /* + Description + Plots the full, reduced, and signal spectra + + Arguments + times - as documented in function "plotMWDResults" + energies - as documented in function "plotMWDResults" + binWidths - as documented in function "plotMWDResults" + ERed - as documented in function "plotMWDResults" + EMin - as documented in function "plotMWDResults" + convertkeVToMeV - as documented in function "plotMWDResults" + signalAcceptance - as documented in function "plotMWDResults" + highResolution - as documented in function "plotMWDResults" + timeCuts - as documented in function "plotMWDResults" + flashCut347 - as documented in function "plotMWDResults" + + Variables + eMax - maximum energy in plotting data + eMin - minomum energy in plotting data + eRangeFull - energy range of data to be plotted + eRangeRed - energy range of the reduced spectrum + e347 - energy of 347 keV signal + e844 - energy of 844 keV signal + e1809 - energy of 1809 keV signal + eRange347 - energy range of 347 keV signal + eRange844 - energy range of 844 keV signal + eRange1809 - energy range of 1809 keV signal + eMin347 - minimum of 347 keV signal plot + eMin844 - minimum of 844 keV signal plot + eMin1809 - minimum of 1809 keV signal plot + eMax347 - maximum of 347 keV signal plot + eMax844 - maximum of 844 keV signal plot + eMax1809 - maximum of 1809 keV signal plot + tMicrospill - microspill duration in ns + tMin347 - minimum time for 347keV signal acceptance + tMax347 - maximum time for 347keV signal acceptance + tMod347 - time modulus for 347keV signal acceptance + tMin844 - minimum time for 844keV signal acceptance + tMax844 - maximum time for 844keV signal acceptance + tMod844 - time modulus for 844keV signal acceptance + tMin1809 - minimum time for 1809keV signal acceptance + tMax1809 - maximum time for 1809keV signal acceptance + tMod1809 - time modulus for 1809keV signal acceptance + vars - temporary buffer to store data if it requires unit conversion from keV to MeV + order - order of plot generation histograms + eMinMax - energy ranges of histograms as {minimum, maximum} + eRanges - energy ranges of the histograms as {maximum - minimum} + nOrder - number of plot categories + nBins - bin numbers in histograms + binWidthStrs - bin widths as strings for histograms + unit - energy unit as either keV or MeV + sigRange - controls when the signal specific plot title extensions are added + fileNames - vector of plot file names + titles - vector of plot titles + px - number of x pixels for canvases + py - number of y pixels for canvases + canvases - vector of all canvases + plot - order iterator + lowResAxisTextSize - low resolution axis text size + c - canvases iterator + hists - vector of histograms + nEntries - number of data points to use with histograms + e - energy buffer vector + t - time buffer vector + lx1 - defines an x coordinate for stat boxes and legend + ly1 - defines a y coordinate for stat boxes and legend + lx2 - defines the other x coordinate for stat boxes and legend + ly2 - defines the other y coordinate for stat boxes and legend + count - buffer vector for histogram bin checks + + */ + std::cout << std::string("Generating plots with ") + (highResolution ? "high resolution, " : "low resolution, ") + (timeCuts ? "time cuts, " : "no time cuts, ") + (convertkeVToMeV ? "in MeV" : "in MeV") << std::endl; + + // Define the energy parameters in keV + double eMax = *std::max_element(energies.begin(), energies.end()) + 1; + double eMin = EMin < std::numeric_limits::epsilon() ? *std::min_element(energies.begin(), energies.end()) : EMin; + if (eMin < 0) { + std::cout << "==========Warning==========" << std::endl; + std::cout << "Negative minimum energy of " << *std::min_element(energies.begin(), energies.end()) << " keV" << std::endl; + std::cout << "===========================" << std::endl; + }; + double eRangeFull = eMax - eMin; + double eRangeRed = ERed - eMin; + double e347 = 347, eRange347 = e347 * signalAcceptance, eMin347 = e347 - eRange347, eMax347 = e347 + eRange347; + double e844 = 844, eRange844 = e844 * signalAcceptance, eMin844 = e844 - eRange844, eMax844 = e844 + eRange844; + double e1809 = 1809, eRange1809 = e1809 * signalAcceptance, eMin1809 = e1809 - eRange1809, eMax1809 = e1809 + eRange1809; + + // Double the energy range to include both sides of the acceptance + eRange347 *=2; + eRange844 *=2; + eRange1809 *=2; + + // Define the time parameters in ns + const double tMicroSpill = 1695; + const double tMin347 = flashCut347 ? 200 : 300, tMax347 = 700, tMod347 = tMicroSpill; + const double tMin844 = 492000, tMax844 = 1330000, tMod844 = tMax844; + const double tMin1809 = 500, tMax1809 = 1600, tMod1809 = tMicroSpill; + + // Convert everything to keV if needed + double* vars[] = {&e347, &eMin347, &eMax347, &eRange347, + &e844, &eMin844, &eMax844, &eRange844, + &e1809, &eMin1809, &eMax1809, &eRange1809, + &eMin, &eMax, &ERed, &eRangeFull, &eRangeRed}; + if (convertkeVToMeV) { + double conversion = 1000; + std::transform(std::begin(vars), std::end(vars), std::begin(vars), [](double* x)-> double*{ *x /= 1000; return x;}); + std::transform(energies.begin(), energies.end(), energies.begin(), [conversion](double x) { return x / conversion; }); + std::transform(binWidths.begin(), binWidths.end(), binWidths.begin(), [conversion](double x) { return x / conversion; }); + }; + std::vector order = {"Full", "Reduced", "347", "844", "1809"}; + std::vector> eMinMax = {{eMin, eMax}, {eMin, ERed}, {eMin347, eMax347}, {eMin844, eMax844}, {eMin1809, eMax1809}}; // Energy ranges for TH1Ds + std::vector eRanges = {eRangeFull, eRangeRed, eRange347, eRange844, eRange1809}; + const int nOrder = order.size(); + std::vector nBins(nOrder); + std::transform(eRanges.begin(), eRanges.end(), binWidths.begin(), nBins.begin(), [](double E, double BW) {if (BW == 0) throw std::runtime_error("Illegal zero bin width, exiting."); return E / BW;}); + + // Convert non str types to str types + std::vector binWidthStrs; + for (double binWidth : binWidths) + binWidthStrs.push_back(convertBinWidthToStr(binWidth)); + std::string unit = convertkeVToMeV ? "MeV" : "keV"; + + // Construct file names + const int sigRange = 1; + std::vector fileNames; + std::string signalExtension = std::string("") + (timeCuts ? ".time-cut" : ".no-time-cut") + (flashCut347 ? ".347-flash" : ""); + for (int i = 0; i < nOrder; i++) + fileNames.push_back("MWD." + order[i] + "." + unit + "." + (highResolution ? "high" : "low") + (i > sigRange ? signalExtension : "") + ".png"); + + // Construct plot titles + std::vector titles; + for (int i = 0; i < nOrder; i++) + titles.push_back("Reconstructed signal energy; Energy [" + unit + "]; Count / " + binWidthStrs[i] + " " + unit); + + // Set up the TCanvases + const int px = (highResolution ? 1500 : 750), py = (highResolution ? 1000 : 500); + std::vector canvases; + for (std::string plot : order) + canvases.push_back(new TCanvas(("c" + plot).c_str(), ("c" + plot).c_str(), px, py)); + + // Apply TCanvas formatting + const double lowResAxisTextSize = 0.06; + gStyle->SetTitleFontSize(lowResAxisTextSize); + gStyle->SetTitleAlign(33); + gStyle->SetTitleX(0.75); + if (!highResolution) { + gStyle->SetTitleX(0.99); + for (TCanvas* c : canvases) { + c->SetBottomMargin(0.14); + c->SetLeftMargin(0.175); + c->Update(); + }; + }; + + // Set up TH1Ds + std::vector hists; + for (int i = 0; i < nOrder; i++) + hists.emplace_back(new TH1D((order[i] + " spectrum").c_str(), (order[i] + " spectrum").c_str(), nBins[i], eMinMax[i][0], eMinMax[i][1])); + + // Populate the TH1Ds + const int nEntries = energies.size(); + double e = 0.0, t = 0.0; + for (uint i = 0; i < nEntries; i++) { + e = energies[i]; + t = times[i]; + if (eMin < e && e < eMax) + hists[0]->Fill(e); + if (eMin < e && e < ERed) + hists[1]->Fill(e); + if ((!timeCuts && eMin347 < e && e < eMax347) || (timeCuts && tMin347 < t && t < tMax347 && eMin347 < e && e < eMax347)) + hists[2]->Fill(e); + if ((!timeCuts && eMin844 < e && e < eMax844) || (timeCuts && tMin844 < t && t < tMax844 && eMin844 < e && e < eMax844)) + hists[3]->Fill(e); + if ((!timeCuts && eMin1809 < e && e < eMax1809) || (timeCuts && tMin1809 < t && t < tMax1809 && eMin1809 < e && e < eMax1809)) + hists[4]->Fill(e); + }; + + // Set up the stat boxes + const double lx1 = (highResolution ? 0.7 : 0.6); + const double ly1 = lx1; + const double lx2 = 0.9; + const double ly2 = 0.9; + gStyle->SetStatX(0.9); + gStyle->SetStatY(0.9); + gStyle->SetStatW(lx2 - lx1); + gStyle->SetStatH(ly2 - ly1); + + // Check the hists for overflow and underflow bins + int count = 0; + for (int i = 0; i < nOrder; i++) { + count = hists[i]->GetBinContent(0); + if (count) { + std::cout << "In plot " << fileNames[i] << ", number of underflow bins is " << count << std::endl; + std::cout << "This histogram's lower range is: " << hists[i]->GetXaxis()->GetXmin() << std::endl; + Fatal("makePlot", "Generated plot has values in the underflow bins"); + }; + count = hists[i]->GetBinContent(hists[i]->GetNbinsX() + 1); + if (count) { + std::cout << "In plot " << fileNames[i] << ", number of overflow bins is " << count << std::endl; + std::cout << "This histogram's upper range is: " << hists[i]->GetXaxis()->GetXmax() << std::endl; + Fatal("makePlot", "Generated plot has values in the overflow bins"); + }; + count = hists[i]->GetEntries() - hists[i]->Integral(1, hists[i]->GetNbinsX()); + if (count) { // scaling affects integral but not entries + std::cout << "In plot " << fileNames[i] << ", difference between entries and integral is " << count << std::endl; + Fatal("makePlot", "Generated plot has values in the underflow/overflow bins"); + }; + }; + + // Draw and save the unstacked plots + for (int i = 0; i < nOrder; i++) { + canvases[i]->cd(); + count = hists[i]->GetEntries(); + if (count == 0) { + std::cout << fileNames[i] << " has no entries, not saving" << std::endl; + continue; + }; + hists[i]->SetTitle(titles[i].c_str()); + hists[i]->SetMinimum(0); // Sets the y minimum + hists[i]->Draw("HIST"); + if (!highResolution) { + hists[i]->GetXaxis()->SetLabelSize(lowResAxisTextSize); + hists[i]->GetXaxis()->SetTitleSize(lowResAxisTextSize); + hists[i]->GetYaxis()->SetLabelSize(lowResAxisTextSize); + hists[i]->GetYaxis()->SetTitleSize(lowResAxisTextSize); + }; + canvases[i]->SaveAs(fileNames[i].c_str()); + }; + + // Clean up from the unstacked plots + for (int i = nOrder - 1; i > -1; i--) { + hists[i]->Delete(); + canvases[i]->Close(); + }; + + // Done + std::cout << "Finished\n" << std::endl; + return; +}; + +void makePlots(std::vector ×, std::vector &energies, double ERed, double EMin, const double signalAcceptance, std::vector &binWidths) { + /* + Description + Sets up loops for function "plot" + + Arguments + times - as documented in function "plotMWDResults" + energies - as documented in function "plotMWDResults" + ERed - as documented in function "plotMWDResults" + EMin - as documented in function "plotMWDResults" + binWidths - as documented in function "plotMWDResults" + + Variables + bools - vector of allowed boolean values + highResolution - selects whether plot is generated in high or low resolution + convertkevToMeV - selects whether plot is generated in keV or MeV + */ + + std::vector bools = {true, false}; + for (bool highResolution : bools) { + for (bool convertkeVToMeV : bools) { + plot(times, energies, binWidths, ERed, EMin, convertkeVToMeV, signalAcceptance, highResolution, false, false); + plot(times, energies, binWidths, ERed, EMin, convertkeVToMeV, signalAcceptance, highResolution, true, false); + plot(times, energies, binWidths, ERed, EMin, convertkeVToMeV, signalAcceptance, highResolution, true, true); + }; + }; + + return; +}; + + +void plotMWDResults_noADCClock(std::string fileName, std::string treeName, double ERed = 2000.0, double EMin = 0.0, const double signalAcceptance = 0.1, double binWidthFull = 500, double binWidthRed = 5, double binWidth347 = 5, double binWidth844 = 5, double binWidth1809 = 10) { + /* + Description + Plots spectra measured by the detectors using the MWD algorithm. Names the files as + MWD.....<347-flash/>.png + such that + - determines the range of the plot, with "full" plotting the full spectrum, "red" plotting the spectrum up to ERed, and "signal" corresponding to each signal in the STM energy range + - units of the energy (x) axis + - whether the plot is in high or low resolution + - defines whether or not the relevant time cut has been applied + <347-flash/> - if present, uses the 347keV flash cut instead of the off-flash cut + + Arguments + fileName - name of file containing MWD derived spectra in a ROOT file as a relative path to cwd + treeName - name of tree containing all the data + + Optional arguments + ERed - reduced spectrum max energy [MeV] + eMin - minimum energy of plots + signalAcceptance - for the signal plots, this controls the width of the plotted region as a multiple of the signal energy + binWidthFull - full energy spectrum bin width + binWidthRed - reduced energy spectrum bin width + binWidth347 - bin width of the 347keV signal plot + binWidth844 - bin width of the 844keV signal plot + binWidth1809 - bin width of the 1809keV signal plot + + Variables + times - vector of times converted from ADC clock tick times to real times from fileName + energies - vector of energies measured by detector [keV] + binWidths - vector of bin widths + */ + + // Update global parameters + SetErrorHandler(customErrorHandler); + gROOT->SetBatch(kTRUE); + + // Set up the data storage variables + std::vector times, energies; + + // Collect the data from the file + collectData(fileName, treeName, times, energies); + + // Gather the bin widths together + std::vector binWidths = {binWidthFull, binWidthRed, binWidth347, binWidth844, binWidth1809}; + + // Generate the plots + makePlots(times, energies, ERed, EMin, signalAcceptance, binWidths); + + return; +}; From ffbbc2049d0f679e3e9007b6213e67ebc75aa13f Mon Sep 17 00:00:00 2001 From: Pawel Plesniak Date: Thu, 19 Feb 2026 11:08:32 +0100 Subject: [PATCH 2/2] WIP --- DataSummaries/CountMuCapPerMeasuredPhoton.C | 7 +- DataSummaries/ExpectedPhotonCount.C | 197 +++++++ .../NormalizationComparisonToExpectedCount.C | 184 ++++++ DataSummaries/SignalFractionTable.C | 4 +- Plotting/plotDigitizationStages.C | 16 +- ...WDResults_noADCClock.C => plotHPGeEDeps.C} | 16 +- Plotting/shifting_check.py | 523 ++++++++++++++++++ 7 files changed, 930 insertions(+), 17 deletions(-) create mode 100644 DataSummaries/ExpectedPhotonCount.C create mode 100644 DataSummaries/NormalizationComparisonToExpectedCount.C rename Plotting/{plotMWDResults_noADCClock.C => plotHPGeEDeps.C} (95%) create mode 100644 Plotting/shifting_check.py diff --git a/DataSummaries/CountMuCapPerMeasuredPhoton.C b/DataSummaries/CountMuCapPerMeasuredPhoton.C index 7993a76..b2e5b06 100644 --- a/DataSummaries/CountMuCapPerMeasuredPhoton.C +++ b/DataSummaries/CountMuCapPerMeasuredPhoton.C @@ -100,6 +100,7 @@ void plot(std::vector> capturedMuons, const unsigned long lo muonCaptureCount.push_back(nExpectedMuonCaptures); muonCaptureUncertainty.push_back(uExpectedMuonCaptures); }; + if (capturedMuons[0][0] > std::numeric_limits::epsilon()) { normalizationSource.push_back("347 keV signal"); muonCaptureCount.push_back(capturedMuons[0][0]); @@ -336,6 +337,10 @@ void CountMuCapPerMeasuredPhoton(bool makePlot = false, std::vector signalColumnWidths = {signal347ColumnWidth, signal844ColumnWidth, signal1809ColumnWidth}; + std::string tableTitle = "Normalized muon capture count per measured signal photon"; + std::cout << std::endl; // Buffer line + std::cout << std::string(fullWidth, '=') << std::endl; + std::cout << std::string((fullWidth - tableTitle.size())/2, ' ') << tableTitle << std::endl; std::cout << std::string(fullWidth, '-') << std::endl; // Title line std::cout << std::setw(correctionNameColumnWidth) << std::left << "Correction factor"; for (i = 0; i < nOrder; i++) @@ -372,7 +377,7 @@ void CountMuCapPerMeasuredPhoton(bool makePlot = false, std::vector kInfo) + exit(1); +}; + +std::string doubleToStringScientific(double value, int sigFigs) { + /* + Description + Converts a double to a string with a defined number of significant figures + + Arguments + value - value to convert to string + sigFigs - number of significant figures to use + + Variables + out - output string stream used to convert 'value' to a std::string + */ + std::ostringstream out; + out << std::scientific << std::setprecision(sigFigs - 1) << value; + return out.str(); +}; + +std::string doubleToStringFixed(double value, int sigFigs) { + /* + Description + Converts a double to a string with a defined number of significant figures + + Arguments + value - value to convert to string + sigFigs - number of significant figures to use + + Variables + out - output string stream used to convert 'value' to a std::string + */ + std::ostringstream out; + out << std::fixed << std::setprecision(sigFigs - 1) << value; + return out.str(); +}; + +void ExpectedPhotonCount(double nPOTs = 2.69e16) { + /* + Description + Generates the table describing the expected number of signal photons for each energy + + Arguments + nPOTs - number of Protons on Target simulated including resampling + + Variables + nSF - number of significant figures to use in the table + correctionNameColumnWidth - width of the correction name column + signal347ColumnWidth - width of the 347 keV signal column + signal844ColumnWidth - width of the 844 keV signal column + signal1809ColumnWidth - width of the 1809 keV signal column + fullWidth - total width of the table + i, j - iterators + ratioMCToDetector - vector of vectors containing the ratio of MC truth normalization to detector response normalization and associated uncertainties for 347, 844, and 1809 keV + */ + + + // Update global parameters + SetErrorHandler(customErrorHandler); + gROOT->SetBatch(kTRUE); + + // Define the table formatting + const int correctionNameColumnWidth = 40, signal347ColumnWidth = 30, signal844ColumnWidth = 30, signal1809ColumnWidth = 30, fullWidth = correctionNameColumnWidth + signal347ColumnWidth + signal844ColumnWidth + signal1809ColumnWidth; + const int nSF = 4; + + // Define the signal order + std::vector order = {"347", "844", "1809"}; + const int nOrder = order.size(); + + // Declare iterator variables + int i = 0, j = 0; + + // Define the parameters that contribute to the number of measured signal photons per muon capture, defined as {value, uncertainty} for each signal photon + std::vector correctionFactorNames = { + // "Probability of final state", + "Absorber acceptance", + "Detector acceptance", + "Path attenuation", + "GEANT rate correction", + "Energy window acceptance", + "Clipping factor", + "Geometric acceptance", + "Time cut acceptance" + }; + + // Define the correction factors and their associated uncertianties + const int nCorrectionFactors = correctionFactorNames.size(); + // 347 corr uncert 844 corr uncert 1809 corr uncert + std::vector> pFinalState = {{1.31, 0.013}, {0.093, 0.007}, {0.51, 0.05} }; + std::vector> absorberAcceptance = {{0.87, 0}, {1, 0}, {1, 0} }; + std::vector> detectorAcceptance = {{0.628, 1.528e-4}, {0.288, 1.432e-4}, {0.179, 1.212e-4} }; + std::vector> pathAttenuation = {{1, 0}, {1, 0}, {1, 0} }; + std::vector> geantRateCorrection = {{1, 0}, {0.259, 0}, {1.0, 0} }; + std::vector> signalInEnergyWindow = {{0.67, 0}, {1, 0}, {1, 0} }; + std::vector> clippingFactor = {{0.85, 0}, {1, 0}, {0.85, 0} }; + std::vector> geometricAcceptance = {{3.25e-9, 0}, {3.25e-9, 0}, {3.25e-9, 0} }; + std::vector> timeCutAcceptance = {{0.9976, 0}, {0.6298, 0}, {0.68, 0} }; + + + // Calculate the ratio of MC truth normalization to expected number of captured muons and its errors + std::vector> expectation = {{0, 0}, {0, 0}, {0, 0}}; + const double pMuonStopMDC2020 = 1432535.0 / (2e8 * (4e8 / 869305)); // Based on the MDC2020 workflow + const double uMuonStopMDC2020 = std::sqrt(1432535) / (2e8 * (4e8 / 869305)); // Based on the MDC2020 workflow + const double pMuonCapture = 0.61; + const double uMuonCapture = 0.001; + const double nExpectedMuonCaptures = nPOTs * pMuonStopMDC2020 * pMuonCapture; + const double uExpectedMuonCaptures = nExpectedMuonCaptures * std::sqrt(std::pow(uMuonStopMDC2020/pMuonStopMDC2020, 2) + std::pow(uMuonCapture/pMuonCapture, 2)); + + for (int i = 0; i < 3; i++) { + expectation[i][0] = pFinalState[i][0] * nExpectedMuonCaptures; + expectation[i][1] = expectation[i][0] * std::sqrt(std::pow(pFinalState[i][1]/pFinalState[i][0], 2) + std::pow(uExpectedMuonCaptures/nExpectedMuonCaptures, 2)); + }; + + // Store all the correction factors in a single variable + std::vector>> correctionFactors = {absorberAcceptance, detectorAcceptance, pathAttenuation, geantRateCorrection, signalInEnergyWindow, clippingFactor, geometricAcceptance, timeCutAcceptance}; + + // Incorporate the effect of the correction factors into the number of measured signal photons per muon capture + std::vector> signalPhotonCount(expectation.begin(), expectation.end()); // Initialize the number of measured signal photons per muon capture to the expected number of signal photons per muon capture before corrections + for (i = 0; i < nOrder; i++) { + signalPhotonCount[i][1] = std::pow(expectation[i][1]/expectation[i][0], 2); // Initialize the quadratic sum of the relative uncertainty to the relative uncertainty on the expected number of signal photons per muon capture + }; + + // Incorporate the effect of the correction factors into the number of measured photons per muon capture + for (std::vector> correctionFactor : correctionFactors) { + for (i = 0; i < nOrder; i++) { + // Update the rate + signalPhotonCount[i][0] *= correctionFactor[i][0]; + + // Update the quadratic sum of the uncertainty + if (correctionFactor[i][0] < std::numeric_limits::epsilon()) + Fatal("CountMuCapPerMeasuredPhoton", "Declared correction factor is zero, this would mean we don't get any signal. Did you mean to set this to unity? Fix this!"); + signalPhotonCount[i][1] += std::pow(correctionFactor[i][1]/correctionFactor[i][0], 2); + }; + }; + + // Finalize the uncertainty + for (i = 0; i < nOrder; i++) + signalPhotonCount[i][1] = signalPhotonCount[i][0] * std::sqrt(signalPhotonCount[i][1]); + + + // Print the title line and rules + const std::vector signalColumnWidths = {signal347ColumnWidth, signal844ColumnWidth, signal1809ColumnWidth}; + std::string tableTitle = "Expected signal photon count"; + std::cout << "Probability of a POT producing a captured muon: " << pMuonStopMDC2020 << " ± " << uMuonStopMDC2020 << std::endl; + std::cout << std::endl; // Buffer line + std::cout << std::string(fullWidth, '=') << std::endl; + std::cout << std::string((fullWidth - tableTitle.size())/2, ' ') << tableTitle << std::endl; + std::cout << std::string(fullWidth, '-') << std::endl; // Title line + std::cout << std::setw(correctionNameColumnWidth) << std::left << "nPOTs"; + for (i = 0; i < nOrder; i++) + std::cout << std::setw(signalColumnWidths[i]) << std::left << doubleToStringScientific(nPOTs, nSF); + std::cout << std::endl; + std::cout << std::setw(correctionNameColumnWidth) << std::left << "nCapturedMuons"; + for (i = 0; i < nOrder; i++) + std::cout << std::setw(signalColumnWidths[i]) << std::left << doubleToStringScientific(nExpectedMuonCaptures, nSF) + " ± " + doubleToStringScientific(uExpectedMuonCaptures, nSF); + std::cout << std::endl; + std::cout << std::setw(correctionNameColumnWidth) << std::left << "Expected generated signal photons"; + for (i = 0; i < nOrder; i++) + std::cout << std::setw(signalColumnWidths[i]) << std::left << doubleToStringScientific(expectation[i][0], nSF) + " ± " + doubleToStringScientific(expectation[i][1], nSF); + std::cout << std::endl; + std::cout << std::string(fullWidth, '-') << std::endl; // Section line + + // Print the correction factors and bottom rule + for (i = 0; i < nCorrectionFactors; i++) { + std::cout << std::setw(correctionNameColumnWidth) << std::left << correctionFactorNames[i]; + if (correctionFactorNames[i] == "Geometric acceptance") { + for (j = 0; j < nOrder; j++) + std::cout << std::setw(signalColumnWidths[j]) << std::left << doubleToStringScientific(correctionFactors[i][j][0], nSF) + " ± " + doubleToStringScientific(correctionFactors[i][j][1], nSF); + } + else { + for (j = 0; j < nOrder; j++) + std::cout << std::setw(signalColumnWidths[j]) << std::left << doubleToStringFixed(correctionFactors[i][j][0], nSF) + " ± " + doubleToStringFixed(correctionFactors[i][j][1], nSF); + } + std::cout << std::endl; + }; + std::cout << std::string(fullWidth, '-') << std::endl; // Section line + + // Print the normalized quantities and bottom rule + std::cout << std::setw(correctionNameColumnWidth) << std::left << "Expected signal photons"; + for (i = 0; i < nOrder; i++) + std::cout << std::setw(signalColumnWidths[i]) << std::left << doubleToStringScientific(signalPhotonCount[i][0], nSF) + " ± " + doubleToStringScientific(signalPhotonCount[i][1], nSF); + std::cout << std::endl; + std::cout << std::string(fullWidth, '=') << std::endl; // End line + std::cout << std::endl; // Buffer line +}; \ No newline at end of file diff --git a/DataSummaries/NormalizationComparisonToExpectedCount.C b/DataSummaries/NormalizationComparisonToExpectedCount.C new file mode 100644 index 0000000..644d2b1 --- /dev/null +++ b/DataSummaries/NormalizationComparisonToExpectedCount.C @@ -0,0 +1,184 @@ +// Generates the table of ratio of any normalization type to the expected based on the POT count +// Usage example - $ root -l -q 'NormalizationComparisonToExpectedCount.C({{1.679e10, 3.02e9}, {0.0, 0.0}, {1.583e11, 2.28e10}}, 2.69e16)' +// Original author - Pawel Plesniak + + +void customErrorHandler(int level, Bool_t abort, const char* location, const char* message) { + /* + Description + Define a custom error handler that won't print the stack trace but will print an error message and exit. + */ + std::cerr << message << std::endl; + if (level > kInfo) + exit(1); +}; + +std::string doubleToStringScientific(double value, int sigFigs) { + /* + Description + Converts a double to a string with a defined number of significant figures + + Arguments + value - value to convert to string + sigFigs - number of significant figures to use + + Variables + out - output string stream used to convert 'value' to a std::string + */ + std::ostringstream out; + out << std::scientific << std::setprecision(sigFigs - 1) << value; + return out.str(); +}; + +std::string doubleToStringFixed(double value, int sigFigs) { + /* + Description + Converts a double to a string with a defined number of significant figures + + Arguments + value - value to convert to string + sigFigs - number of significant figures to use + + Variables + out - output string stream used to convert 'value' to a std::string + */ + std::ostringstream out; + out << std::fixed << std::setprecision(sigFigs - 1) << value; + return out.str(); +}; + +void NormalizationComparisonToExpectedCount(std::vector> normalizations, const double nPOTs = 2.69e16) { + /* + Description + Generates the table describing the ratio of a normalization method to the expected number of captured muons + + Arguments + normalizations - vector of vectors containing the normalizations and associated uncertainties for 347, 844, and 1809 keV + nPOTs - number of Protons on Target simulated including resampling + + Variables + nSF - number of significant figures to use in the table + correctionNameColumnWidth - width of the correction name column + signal347ColumnWidth - width of the 347 keV signal column + signal844ColumnWidth - width of the 844 keV signal column + signal1809ColumnWidth - width of the 1809 keV signal column + fullWidth - total width of the table + i, j - iterators + expectation - vector of vectors containing the expected number of captured muons and associated uncertainties for 347, 844, and 1809 keV + signalProbability - vector of vectors containing the probability of final state for 347, 844, and 1809 keV + */ + + // Update global parameters + SetErrorHandler(customErrorHandler); + gROOT->SetBatch(kTRUE); + + // Define the signal order + std::vector order = {"347", "844", "1809"}; + const int nOrder = order.size(); + + // Declare iterator variables + int i = 0, j = 0; + + // Define the parameters that contribute to the number of measured signal photons per muon capture, defined as {value, uncertainty} for each signal photon + std::vector correctionFactorNames = { + // "Probability of final state", + "Absorber acceptance", + "Detector acceptance", + "Path attenuation", + "GEANT rate correction", + "Energy window acceptance", + "Clipping factor", + "Geometric acceptance", + "Time cut acceptance" + }; + + // Define the correction factors and their associated uncertianties + const int nCorrectionFactors = correctionFactorNames.size(); + // 347 corr uncert 844 corr uncert 1809 corr uncert + std::vector> pFinalState = {{1.31, 0.013}, {0.093, 0.007}, {0.51, 0.05} }; + std::vector> absorberAcceptance = {{0.87, 0}, {1, 0}, {1, 0} }; + std::vector> detectorAcceptance = {{0.628, 1.528e-4}, {0.288, 1.432e-4}, {0.179, 1.212e-4} }; + std::vector> pathAttenuation = {{1, 0}, {1, 0}, {1, 0} }; + std::vector> geantRateCorrection = {{1, 0}, {0.259, 0}, {1.0, 0} }; + std::vector> signalInEnergyWindow = {{0.67, 0}, {1, 0}, {1, 0} }; + std::vector> clippingFactor = {{0.85, 0}, {1, 0}, {0.85, 0} }; + std::vector> geometricAcceptance = {{3.25e-9, 0}, {3.25e-9, 0}, {3.25e-9, 0} }; + std::vector> timeCutAcceptance = {{0.9976, 0}, {0.6298, 0}, {0.68, 0} }; + + + // Calculate the ratio of MC truth normalization to expected number of captured muons and its errors + std::vector> expectation = {{0, 0}, {0, 0}, {0, 0}}; + const double pMuonStopMDC2020 = 1432535.0 / (2e8 * (4e8 / 869305)); // Based on the MDC2020 workflow + const double uMuonStopMDC2020 = std::sqrt(1432535) / (2e8 * (4e8 / 869305)); // Based on the MDC2020 workflow + const double pMuonCapture = 0.61; + const double uMuonCapture = 0.001; + const double nExpectedMuonCaptures = nPOTs * pMuonStopMDC2020 * pMuonCapture; + const double uExpectedMuonCaptures = nExpectedMuonCaptures * std::sqrt(std::pow(uMuonStopMDC2020/pMuonStopMDC2020, 2) + std::pow(uMuonCapture/pMuonCapture, 2)); + + for (int i = 0; i < 3; i++) { + expectation[i][0] = pFinalState[i][0] * nExpectedMuonCaptures; + expectation[i][1] = expectation[i][0] * std::sqrt(std::pow(pFinalState[i][1]/pFinalState[i][0], 2) + std::pow(uExpectedMuonCaptures/nExpectedMuonCaptures, 2)); + }; + + // Store all the correction factors in a single variable + std::vector>> correctionFactors = {absorberAcceptance, detectorAcceptance, pathAttenuation, geantRateCorrection, signalInEnergyWindow, clippingFactor, geometricAcceptance, timeCutAcceptance}; + + // Incorporate the effect of the correction factors into the number of measured signal photons per muon capture + std::vector> signalPhotonCount(expectation.begin(), expectation.end()); // Initialize the number of measured signal photons per muon capture to the expected number of signal photons per muon capture before corrections + for (i = 0; i < nOrder; i++) { + signalPhotonCount[i][1] = std::pow(expectation[i][1]/expectation[i][0], 2); // Initialize the quadratic sum of the relative uncertainty to the relative uncertainty on the expected number of signal photons per muon capture + }; + + // Incorporate the effect of the correction factors into the number of measured photons per muon capture + for (std::vector> correctionFactor : correctionFactors) { + for (i = 0; i < nOrder; i++) { + // Update the rate + signalPhotonCount[i][0] *= correctionFactor[i][0]; + + // Update the quadratic sum of the uncertainty + if (correctionFactor[i][0] < std::numeric_limits::epsilon()) + Fatal("CountMuCapPerMeasuredPhoton", "Declared correction factor is zero, this would mean we don't get any signal. Did you mean to set this to unity? Fix this!"); + signalPhotonCount[i][1] += std::pow(correctionFactor[i][1]/correctionFactor[i][0], 2); + }; + }; + + // Finalize the uncertainty + for (i = 0; i < nOrder; i++) + signalPhotonCount[i][1] = signalPhotonCount[i][0] * std::sqrt(signalPhotonCount[i][1]); + + + // Define the signal photon energy strings + std::string e347 = "347 keV"; + std::string e844 = "844 keV"; + std::string e1809 = "1809 keV"; + + // Define table formatting + const int nSF = 4; + const int correctionNameColumnWidth = 20, signalColumnWidth = 30; + const int fullWidth = correctionNameColumnWidth + signalColumnWidth * 3; + + // Print the title line and rules + std::string title = "Expected to normalization signal photon count comparison for " + doubleToStringScientific(nPOTs, nSF) + " POTs"; + std::cout << std::endl; // Buffer line + std::cout << std::string(fullWidth, '=') << std::endl; // Title line + std::cout << std::string((fullWidth - static_cast(title.size())) / 2, ' ') << title << std::endl; + std::cout << std::string(fullWidth, '=') << std::endl; // Title line + std::cout << std::setw(correctionNameColumnWidth) << std::left << "" << std::setw(signalColumnWidth) << std::left << "347 keV" << std::setw(signalColumnWidth) << std::left << "844 keV" << std::setw(signalColumnWidth) << std::left << "1809 keV" << std::endl; + std::cout << std::string(fullWidth, '-') << std::endl; // Section line + std::cout << std::setw(correctionNameColumnWidth) << std::left << "Expected values";; + for (int j = 0; j < 3; j++) + std::cout << std::setw(signalColumnWidth) << std::left << doubleToStringScientific(signalPhotonCount[j][0], nSF) + " ± " + doubleToStringScientific(signalPhotonCount[j][1], nSF); + std::cout << std::endl; + std::cout << std::setw(correctionNameColumnWidth) << std::left << "Normalization";; + for (int j = 0; j < 3; j++) + std::cout << std::setw(signalColumnWidth) << std::left << doubleToStringScientific(normalizations[j][0], nSF) + " ± " + doubleToStringScientific(normalizations[j][1], nSF); + std::cout << std::endl; + std::cout << std::string(fullWidth, '-') << std::endl; // Section line + std::cout << std::setw(correctionNameColumnWidth) << std::left << "Ratio";; + for (int j = 0; j < 3; j++) + std::cout << std::setw(signalColumnWidth) << std::left << doubleToStringFixed(normalizations[j][0]/signalPhotonCount[j][0], nSF) + " ± " + doubleToStringFixed(std::sqrt(std::pow(normalizations[j][1]/normalizations[j][0], 2) + std::pow(signalPhotonCount[j][1]/signalPhotonCount[j][0], 2)), nSF); + std::cout << std::endl; + std::cout << std::string(fullWidth, '=') << std::endl; // End line + std::cout << std::endl; // Buffer line + return; +}; \ No newline at end of file diff --git a/DataSummaries/SignalFractionTable.C b/DataSummaries/SignalFractionTable.C index f286510..85f6fe8 100644 --- a/DataSummaries/SignalFractionTable.C +++ b/DataSummaries/SignalFractionTable.C @@ -1,5 +1,5 @@ // Generates the table of signal photons from the total measured signal spectrum -// Usage example - $ root -l -q 'SignalFractionTable.C(2082. 0. 21777)' +// Usage example - $ root -l -q 'SignalFractionTable.C(2082, 0, 21777)' // Original author - Pawel Plesniak void customErrorHandler(int level, Bool_t abort, const char* location, const char* message) { @@ -105,6 +105,7 @@ void SignalFractionTable(double n347, double n844, double n1809) { const int titleWidth = 40, signalWidth = 30, fullWidth = titleWidth + signalWidth * 3; // Print the table + std::cout << std::endl; // Buffer line std::string tableTitle = "Signal Photons Summary Table"; std::cout << std::string(fullWidth, '=') << std::endl; std::cout << std::string((fullWidth - tableTitle.size())/2, ' ') << tableTitle << std::endl; @@ -116,6 +117,7 @@ void SignalFractionTable(double n347, double n844, double n1809) { std::cout << std::string(fullWidth, '-') << std::endl; std::cout << std::setw(titleWidth) << std::left << "Signal Photons: " << std::setw(signalWidth) << std::left << sig347 << std::setw(signalWidth) << std::left << sig844 << std::setw(signalWidth) << std::left << sig1809 << std::endl; std::cout << std::string(fullWidth, '=') << std::endl; + std::cout << std::endl; // Buffer line return; }; \ No newline at end of file diff --git a/Plotting/plotDigitizationStages.C b/Plotting/plotDigitizationStages.C index 0e79cd6..aa8103f 100644 --- a/Plotting/plotDigitizationStages.C +++ b/Plotting/plotDigitizationStages.C @@ -86,7 +86,7 @@ void collectData(const std::string fileName, const std::string treeName, std::ve int entries = tree->GetEntries(); // Collect the data - for (int i = 0; i < entries; i++) { + for (int i = 1000; i < 2000; i++) { // Get the corresponding entry tree->GetEntry(i); @@ -154,14 +154,14 @@ void plot(std::vector &chargeCollected, std::vector &chargeDecay // Construct the time vector for plotting std::vector time; const double tADC = 3.125; - for (int i = 0; i < N; i++) time.push_back(i * tADC); + for (int i = 1000; i < (N+1000); i++) time.push_back(i * tADC); double* timeA = time.data(); // Generate the plot for the collected charge if (!chargeCollected.empty()) { TCanvas* cChargeCollected = new TCanvas("cChargeCollected", "chargeCollected", 800, 600); - gPad->SetLeftMargin(0.14); - gPad->SetRightMargin(0.02); + // gPad->SetLeftMargin(0.14); + // gPad->SetRightMargin(0.02); TGraph *gChargeCollected = new TGraph(N, timeA, chargeCollectedA); gChargeCollected->SetTitle("Charge collected;Time [ns];Charge [e]"); gChargeCollected->Draw("APL"); @@ -174,8 +174,8 @@ void plot(std::vector &chargeCollected, std::vector &chargeDecay // Generate the plot for the decayed charge if (!chargeDecayed.empty()) { TCanvas* cChargeDecayed = new TCanvas("cChargeDecayed", "chargeDecayed", 800, 600); - gPad->SetLeftMargin(0.14); - gPad->SetRightMargin(0.02); + // gPad->SetLeftMargin(0.14); + // gPad->SetRightMargin(0.02); TGraph *gChargeDecayed = new TGraph(N, timeA, chargeDecayedA); gChargeDecayed->SetTitle("Charge deacyed;Time [ns];Charge [e]"); gChargeDecayed->Draw("APL"); @@ -188,8 +188,8 @@ void plot(std::vector &chargeCollected, std::vector &chargeDecay // Generate the plot for the digitized waveform if (!ADCs.empty()) { TCanvas* cADCs = new TCanvas("cADCs", "ADCs", 800, 600); - gPad->SetLeftMargin(0.12); - gPad->SetRightMargin(0.02); + // gPad->SetLeftMargin(0.12); + // gPad->SetRightMargin(0.02); TGraph *gADCs = new TGraph(N, timeA, ADCsA); gADCs->SetTitle("Digitized waveform;Time [ns];ADC [arb. unit]"); gADCs->Draw("APL"); diff --git a/Plotting/plotMWDResults_noADCClock.C b/Plotting/plotHPGeEDeps.C similarity index 95% rename from Plotting/plotMWDResults_noADCClock.C rename to Plotting/plotHPGeEDeps.C index a76730f..21c177f 100644 --- a/Plotting/plotMWDResults_noADCClock.C +++ b/Plotting/plotHPGeEDeps.C @@ -1,5 +1,5 @@ -// Generates the analysis path of the MWD algorithm -// Usage example - $ root -l -q 'plotMWDResults.C("FinalData/DigNoZS.root", "MWDSpectra/ttree", 2000.0, 50, 0.1, 500, 5, 0.5, 0.8, 1.1)' +// Plots the spectra associated with the MC truth information of the energy depositions inside the HPGe detector +// Usage example - $ root -l -q 'plotHPGeEDeps.C("FinalData/DigNoZS.root", "MWDSpectra/ttree", 2000.0, 50, 0.1, 500, 5, 0.5, 0.8, 1.1)' // Original author - Pawel Plesniak #include @@ -75,7 +75,7 @@ void collectData(const std::string fileName, const std::string treeName, std::ve // Collect the data times.push_back(dataTime); - energies.push_back(dataE); + energies.push_back(dataE * 1000); // Convert from MeV to keV }; // Check if data has been collected if (times.size() == 0 || energies.size() == 0) @@ -200,6 +200,8 @@ void plot(std::vector times, std::vector energies, std::vector times, std::vector energies, std::vectorFill(e); if (eMin < e && e < ERed) hists[1]->Fill(e); - if ((!timeCuts && eMin347 < e && e < eMax347) || (timeCuts && tMin347 < t && t < tMax347 && eMin347 < e && e < eMax347)) + if ((!timeCuts && eMin347 < e && e < eMax347) || (timeCuts && tMin347 < fmod(t, tMod347) && fmod(t, tMod347) < tMax347 && eMin347 < e && e < eMax347)) hists[2]->Fill(e); - if ((!timeCuts && eMin844 < e && e < eMax844) || (timeCuts && tMin844 < t && t < tMax844 && eMin844 < e && e < eMax844)) + if ((!timeCuts && eMin844 < e && e < eMax844) || (timeCuts && tMin844 < fmod(t, tMod844) && fmod(t, tMod844) < tMax844 && eMin844 < e && e < eMax844)) hists[3]->Fill(e); - if ((!timeCuts && eMin1809 < e && e < eMax1809) || (timeCuts && tMin1809 < t && t < tMax1809 && eMin1809 < e && e < eMax1809)) + if ((!timeCuts && eMin1809 < e && e < eMax1809) || (timeCuts && tMin1809 < fmod(t, tMod1809) && fmod(t, tMod1809) < tMax1809 && eMin1809 < e && e < eMax1809)) hists[4]->Fill(e); }; @@ -390,7 +392,7 @@ void makePlots(std::vector ×, std::vector &energies, double }; -void plotMWDResults_noADCClock(std::string fileName, std::string treeName, double ERed = 2000.0, double EMin = 0.0, const double signalAcceptance = 0.1, double binWidthFull = 500, double binWidthRed = 5, double binWidth347 = 5, double binWidth844 = 5, double binWidth1809 = 10) { +void plotHPGeEDeps(std::string fileName, std::string treeName, double ERed = 2000.0, double EMin = 0.0, const double signalAcceptance = 0.1, double binWidthFull = 500, double binWidthRed = 5, double binWidth347 = 5, double binWidth844 = 5, double binWidth1809 = 10) { /* Description Plots spectra measured by the detectors using the MWD algorithm. Names the files as diff --git a/Plotting/shifting_check.py b/Plotting/shifting_check.py new file mode 100644 index 0000000..ce86161 --- /dev/null +++ b/Plotting/shifting_check.py @@ -0,0 +1,523 @@ +""" +Plot particle distributions at shifting stage. +""" + +# TODOs +# Add the stats boxes to the plots +# Fix the legend style + +import array +import glob +import os + +import ROOT +import pandas as pd +import uproot +import numpy as np + +INPUT_DIRS = { + "Ele": "InputDatasetBiasing/output_data/S1/Ele/StepPointMCs/", + "Mu": "InputDatasetBiasing/output_data/S1/Mu/StepPointMCs/", +} + +TBRANCH_NAMES = ["virtualdetectorId", "x", "y", "Ekin"] + +# Plot settings +N_BINS = 100 +INDIVIDUAL_PLOT_PDGIDS = [11, -11, 13, -13, 22, 2112] + +# Data handling variables +TREE_NAME = "ROOTDump/ttree" +ROOT.gROOT.SetBatch(True) # pylint: disable=no-member +ROOT.gStyle.SetOptStat(0) # pylint: disable=no-member +cLow = ROOT.TCanvas("cLow", "", 800, 600) # pylint: disable=no-member +cHigh = ROOT.TCanvas("cHigh", "", 1600, 1200) # pylint: disable=no-member + +# Boring stuff to make the code more legible +BOOL_VALUES = [True, False] + +# The S1 simulation of EleBeamCat and MuBeamCat have effective POT counts, define the ratio of these: +MU_ELE_RATIO = 9.95e12 / 2.37e11 + +# Define the titles for the TH1D plots +TH1D_PLOT_TITLES = { + "VD101": "Energy at VD101", + "HPGe Absorber": "Energy at HPGe Absorber Region", + "LaBr Absorber": "Energy at LaBr Absorber Region", + "HPGe SSC": "Energy at HPGe SSC Aperture Region", + "LaBr SSC": "Energy at LaBr SSC Aperture Region", +} + +# Define the legend entries and titles for overlap plots +OVERLAP_PLOT_LEGEND_ENTRIES_AND_TITLE = { + "HPGe SSC vs LaBr SSC": { + "legend": { + "Ele": "HPGe in SSC Aperture", + "LaBr": "LaBr in SSC Aperture", + }, + "title": "SSC Aperture energy distributions", + }, + "HPGe Absorber vs LaBr Absorber": { + "legend": { + "Ele": "HPGe in Absorber Aperture", + "LaBr": "LaBr in Absorber Aperture", + }, + "title": "Absorber Aperture energy distributions", + }, + "VD101 vs HPGe Absorber": { + "legend": { + "Ele": "VD101 Data", + "LaBr": "HPGe in Absorber Aperture", + }, + "title": "VD101 vs HPGe Absorber energy distributions", + }, + "HPGe Absorber vs HPGe SSC": { + "legend": { + "Ele": "HPGe in Absorber Aperture", + "LaBr": "HPGe in SSC Aperture", + }, + "title": "HPGe Absorber vs SSC energy distributions", + }, + "VD101 vs LaBr Absorber": { + "legend": { + "Ele": "VD101 Data", + "LaBr": "LaBr in Absorber Aperture", + }, + "title": "VD101 vs LaBr Absorber energy distributions", + }, + "LaBr Absorber vs SSC": { + "legend": { + "Ele": "LaBr in Absorber Aperture", + "LaBr": "LaBr in SSC Aperture", + }, + "title": "LaBr Absorber vs SSC energy distributions", + }, + "VD101 vs HPGe SSC": { + "legend": { + "Ele": "VD101 Data", + "LaBr": "HPGe in SSC Aperture", + }, + "title": "VD101 vs HPGe SSC energy distributions", + }, + "VD101 vs LaBr SSC": { + "legend": { + "Ele": "VD101 Data", + "LaBr": "LaBr in SSC Aperture", + }, + "title": "VD101 vs LaBr SSC energy distributions", + }, +} + + +def parse_dir(data_dir: str) -> pd.DataFrame: + """ + Parse all ROOT files in the specified directory and return a DataFrame + of particles going through FILTER_VDID. + """ + all_dfs = [] + files = sorted(glob.glob(f"{data_dir}/*.root")) + num_files = len(files) + if num_files == 0: + raise ValueError(f"No ROOT files found in directory: {data_dir}") + counter = 1 + for input_file in files: + print(f"Reading file {counter}/{num_files}:", input_file) + counter += 1 + with uproot.open(input_file) as f: + arrs = f[TREE_NAME].arrays(TBRANCH_NAMES, library="np") + if arrs is None or any(len(arrs[param]) == 0 for param in TBRANCH_NAMES): + print( + "Empty arrs:", + [param for param in TBRANCH_NAMES if len(arrs[param]) == 0], + ) + raise ValueError( + f"No entries found in table {TREE_NAME} in file: " f"{input_file}" + ) + save_arrs = {} + for param in TBRANCH_NAMES: + save_arrs[param] = arrs[param] + df_tmp = pd.DataFrame(save_arrs).copy() + if df_tmp.empty: + raise ValueError( + f"DataFrame is empty after loading table {TREE_NAME} in " + f"file: {input_file}" + ) + df_tmp = df_tmp[df_tmp["virtualdetectorId"] == 101] + all_dfs.append(df_tmp) + # break # For debugging, remove later + return pd.concat(all_dfs, ignore_index=True) + + +def read_data_from_dir() -> pd.DataFrame: + """Read data from all ROOT files in the specified directory.""" + + ele_dir = INPUT_DIRS["Ele"] + print("Reading ele files from directory:", ele_dir) + read_ele_df = parse_dir(ele_dir) + if read_ele_df.empty: + raise ValueError(f"No data found in directory: {ele_dir}") + + mu_dir = INPUT_DIRS["Mu"] + print("Reading mu files from directory:", mu_dir) + read_mu_df = parse_dir(mu_dir) + if read_mu_df.empty: + raise ValueError(f"No data found in directory: {mu_dir}") + + return read_ele_df, read_mu_df + + +def convert_plot_title_to_file_name(title_key: str) -> str: + """Convert plot title to a valid file name.""" + return title_key.replace(" ", "_").replace("-", "_").lower() + + +def plot_th1d( + ele_data: pd.DataFrame, + mu_data: pd.DataFrame, + title_key: str, + eMin: float = 0.0, + eMax: float = 2.0, +) -> None: + """Plot TH1D histograms for electron and muon data.""" + + # Get energy values + energies_ele = ele_data["Ekin"].values + energies_mu = mu_data["Ekin"].values + + # Check if there are any energies to plot + if len(energies_ele) + len(energies_mu) == 0: + print(f"Warning: No energy data available for plot {title_key}, skipping.") + return + + # Define the bin width + bin_wid = (eMax - eMin) / N_BINS + + # Create histograms + h_ele = ROOT.TH1D( # pylint: disable=no-member + "hEle", TH1D_PLOT_TITLES[title_key], N_BINS, eMin, eMax + ) + for energy in energies_ele: + h_ele.Fill(energy) + h_mu = ROOT.TH1D( # pylint: disable=no-member + "hMu", TH1D_PLOT_TITLES[title_key], N_BINS, eMin, eMax + ) + for energy in energies_mu: + h_mu.Fill(energy) + + # Make these plots reflect the same number of POTs + h_ele.Scale(MU_ELE_RATIO) + + # Check if there is data to merge + if h_ele.GetEntries() + h_mu.GetEntries() == 0: + print(f"Warning: Both histograms for plot {title_key} are empty, skipping.") + return + + # Define the file name + file_name = convert_plot_title_to_file_name(TH1D_PLOT_TITLES[title_key]) + + # Create a merged histogram + h_total = h_ele.Clone("hTotal") + h_total.Add(h_mu) + + # Define the plot title and axis labels + main_title = TH1D_PLOT_TITLES[title_key] + x_title = "Kinetic Energy (MeV)" + y_title = f"Counts / {bin_wid:.2f} MeV" + h_total.SetTitle(f"{main_title};{x_title};{y_title}") + + # Plot and save the low res canvas + cLow.cd() + h_total.Draw("HIST") + cLow.SaveAs(f"{file_name}_th1d_low.png") + + # Plot and save the high res canvas + cHigh.cd() + h_total.Draw("HIST") + cHigh.SaveAs(f"{file_name}_th1d_high.png") + + +def plot_th1d_overlap( + ele_hpge_data: pd.DataFrame, + mu_hpge_data: pd.DataFrame, + ele_labr_data: pd.DataFrame, + mu_labr_data: pd.DataFrame, + title_key: str, + eMin: float = 0.0, + eMax: float = 2.0, +) -> None: + """Plot TH1D histograms for electron and muon data.""" + + # Get energy values + energies_ele_hpge = ele_hpge_data["Ekin"].values + energies_mu_hpge = mu_hpge_data["Ekin"].values + energies_ele_labr = ele_labr_data["Ekin"].values + energies_mu_labr = mu_labr_data["Ekin"].values + + # Check if there are any energies to plot + if ( + len(energies_ele_hpge) + + len(energies_mu_hpge) + + len(energies_ele_labr) + + len(energies_mu_labr) + == 0 + ): + print(f"Warning: No energy data available for plot {title_key}, skipping.") + return + + # Define the bin width + bin_wid = (eMax - eMin) / N_BINS + + # Create histograms for HPGe data + h_ele_hpge = ROOT.TH1D( # pylint: disable=no-member + "hEleHPGe", + OVERLAP_PLOT_LEGEND_ENTRIES_AND_TITLE[title_key]["title"], + N_BINS, + eMin, + eMax, + ) + for energy in energies_ele_hpge: + h_ele_hpge.Fill(energy) + h_mu_hpge = ROOT.TH1D( # pylint: disable=no-member + "hMuHPGe", + OVERLAP_PLOT_LEGEND_ENTRIES_AND_TITLE[title_key]["title"], + N_BINS, + eMin, + eMax, + ) + for energy in energies_mu_hpge: + h_mu_hpge.Fill(energy) + + # Create histograms for LaBr data + h_ele_labr = ROOT.TH1D( # pylint: disable=no-member + "hEleLabr", + OVERLAP_PLOT_LEGEND_ENTRIES_AND_TITLE[title_key]["title"], + N_BINS, + eMin, + eMax, + ) + for energy in energies_ele_labr: + h_ele_labr.Fill(energy) + h_mu_labr = ROOT.TH1D( # pylint: disable=no-member + "hMuLabr", + OVERLAP_PLOT_LEGEND_ENTRIES_AND_TITLE[title_key]["title"], + N_BINS, + eMin, + eMax, + ) + for energy in energies_mu_labr: + h_mu_labr.Fill(energy) + + # Make these plots reflect the same number of POTs + h_ele_hpge.Scale(MU_ELE_RATIO) + h_ele_labr.Scale(MU_ELE_RATIO) + + # Check if there is data to merge + if ( + h_ele_hpge.GetEntries() + + h_mu_hpge.GetEntries() + + h_ele_labr.GetEntries() + + h_mu_labr.GetEntries() + == 0 + ): + print(f"Warning: Both histograms for plot {title_key} are empty, skipping.") + return + + # Define the file name + file_name = convert_plot_title_to_file_name( + OVERLAP_PLOT_LEGEND_ENTRIES_AND_TITLE[title_key]["title"] + ) + + # Create a merged histograms + h_total_hpge = h_ele_hpge.Clone("hTotalHPGe") + h_total_hpge.Add(h_mu_hpge) + h_total_labr = h_ele_labr.Clone("hTotalLabr") + h_total_labr.Add(h_mu_labr) + h_total_hpge.SetLineColor(ROOT.kRed) # pylint: disable=no-member + h_total_labr.SetLineColor(ROOT.kBlue) # pylint: disable=no-member + h_total_hpge.SetLineWidth(2) + h_total_labr.SetLineWidth(2) + + # Define the THStack + hs = ROOT.THStack("hs", "") + hs.Add(h_total_hpge) + hs.Add(h_total_labr) + + # Define the legend + legend = ROOT.TLegend(0.65, 0.7, 0.85, 0.85) # pylint: disable=no-member + legend.SetBorderSize(0) + legend.SetFillStyle(0) + legend.AddEntry( + h_total_hpge, + OVERLAP_PLOT_LEGEND_ENTRIES_AND_TITLE[title_key]["legend"]["Ele"], + "l", + ) + legend.AddEntry( + h_total_labr, + OVERLAP_PLOT_LEGEND_ENTRIES_AND_TITLE[title_key]["legend"]["LaBr"], + "l", + ) + + # Define the plot title and axis labels + main_title = OVERLAP_PLOT_LEGEND_ENTRIES_AND_TITLE[title_key]["title"] + x_title = "Kinetic Energy (MeV)" + y_title = f"Counts / {bin_wid:.2f} MeV" + hs.SetTitle(f"{main_title};{x_title};{y_title}") + + # Plot and save the low res canvas + cLow.cd() + hs.Draw("nostack HIST") + legend.Draw() + cLow.SaveAs(f"{file_name}_th1d_low.png") + + # Plot and save the high res canvas + cHigh.cd() + hs.Draw("nostack HIST") + legend.Draw() + cHigh.SaveAs(f"{file_name}_th1d_high.png") + + return + + +ele_df, mu_df = read_data_from_dir() + +# Filter data within ±25mm of the aperture center +hpge_ssc_aperture_x = -3945.0 +labr_ssc_aperture_x = -3865.0 +aperture_tolerance = 25.0 + +ele_hpge_absorber_df = ele_df[ + (ele_df["x"] >= hpge_ssc_aperture_x - aperture_tolerance) + & (ele_df["x"] <= hpge_ssc_aperture_x + aperture_tolerance) + & (ele_df["y"] >= -aperture_tolerance) + & (ele_df["y"] <= aperture_tolerance) +] +ele_labr_absorber_df = ele_df[ + (ele_df["x"] >= labr_ssc_aperture_x - aperture_tolerance) + & (ele_df["x"] <= labr_ssc_aperture_x + aperture_tolerance) + & (ele_df["y"] >= -aperture_tolerance) + & (ele_df["y"] <= aperture_tolerance) +] + +mu_hpge_absorber_df = mu_df[ + (mu_df["x"] >= hpge_ssc_aperture_x - aperture_tolerance) + & (mu_df["x"] <= hpge_ssc_aperture_x + aperture_tolerance) + & (mu_df["y"] >= -aperture_tolerance) + & (mu_df["y"] <= aperture_tolerance) +] +mu_labr_absorber_df = mu_df[ + (mu_df["x"] >= labr_ssc_aperture_x - aperture_tolerance) + & (mu_df["x"] <= labr_ssc_aperture_x + aperture_tolerance) + & (mu_df["y"] >= -aperture_tolerance) + & (mu_df["y"] <= aperture_tolerance) +] + + +# Filter data within radius of aperture centers +aperture_radius = np.sqrt(50) / np.pi + +ele_hpge_ssc_df = ele_df[ + (ele_df["x"] - hpge_ssc_aperture_x) ** 2 + ele_df["y"] ** 2 <= aperture_radius**2 +] +ele_labr_ssc_df = ele_df[ + (ele_df["x"] - labr_ssc_aperture_x) ** 2 + ele_df["y"] ** 2 <= aperture_radius**2 +] + +mu_hpge_ssc_df = mu_df[ + (mu_df["x"] - hpge_ssc_aperture_x) ** 2 + mu_df["y"] ** 2 <= aperture_radius**2 +] +mu_labr_ssc_df = mu_df[ + (mu_df["x"] - labr_ssc_aperture_x) ** 2 + mu_df["y"] ** 2 <= aperture_radius**2 +] + +# Generate energy distribution plots at VD101 +# plot_th1d( +# ele_df, +# mu_df, +# "VD101", +# ) + +# # Generate energy distribution plots at HPGe absorber +# plot_th1d( +# ele_hpge_absorber_df, +# mu_hpge_absorber_df, +# "HPGe Absorber", +# ) + +# # Generate energy distribution plots at LaBr absorber +# plot_th1d( +# ele_labr_absorber_df, +# mu_labr_absorber_df, +# "LaBr Absorber", +# ) + +# # Generate energy distribution plots at HPGe SSC aperture +# plot_th1d( +# ele_hpge_ssc_df, +# mu_hpge_ssc_df, +# "HPGe SSC", +# ) + +# # Generate energy distribution plots at LaBr SSC aperture +# plot_th1d( +# ele_labr_ssc_df, +# mu_labr_ssc_df, +# "LaBr SSC", +# ) + + +# Generate overlap plots comparing HPGe and LaBr SSC data +plot_th1d_overlap( + ele_hpge_ssc_df, + mu_hpge_ssc_df, + ele_labr_ssc_df, + mu_labr_ssc_df, + "HPGe SSC vs LaBr SSC", +) + +# plot_th1d_overlap( +# ele_hpge_absorber_df, +# mu_hpge_absorber_df, +# ele_labr_absorber_df, +# mu_labr_absorber_df, +# "HPGe Absorber vs LaBr Absorber", +# ) + +# plot_th1d_overlap( +# ele_hpge_absorber_df, +# mu_hpge_absorber_df, +# ele_hpge_ssc_df, +# mu_hpge_ssc_df, +# "HPGe Absorber vs HPGe SSC", +# ) + +# plot_th1d_overlap( +# ele_labr_absorber_df, +# mu_labr_absorber_df, +# ele_labr_ssc_df, +# mu_labr_ssc_df, +# "LaBr Absorber vs LaBr SSC", +# ) + +# plot_th1d_overlap( +# ele_df, mu_df, ele_hpge_absorber_df, mu_hpge_absorber_df, "VD101 vs HPGe Absorber" +# ) + +# plot_th1d_overlap( +# ele_df, mu_df, ele_labr_absorber_df, mu_labr_absorber_df, "VD101 vs LaBr Absorber" +# ) + +# plot_th1d_overlap( +# ele_df, +# mu_df, +# ele_hpge_ssc_df, +# mu_hpge_ssc_df, +# "VD101 vs HPGe SSC", +# ) + +# plot_th1d_overlap( +# ele_df, +# mu_df, +# ele_labr_ssc_df, +# mu_labr_ssc_df, +# "VD101 vs LaBr SSC", +# )