diff --git a/apps/TrappingCorrectionAm241.cxx b/apps/TrappingCorrectionAm241.cxx index dadf4894..d2283cb5 100644 --- a/apps/TrappingCorrectionAm241.cxx +++ b/apps/TrappingCorrectionAm241.cxx @@ -344,7 +344,6 @@ bool TrappingCorrectionAm241::Analyze() for (unsigned int f = 0; fSetFileNameStripMap(m_StripMapFile); Loader->SetFileName(File); @@ -363,21 +361,18 @@ bool TrappingCorrectionAm241::Analyze() S->SetModule(Loader, MNumber); ++MNumber; - cout<<"Creating TAC calibrator"<SetTACCalFileName(m_TACCalFile); TACCalibrator->SetTACCutFileName(m_TACCutFile); S->SetModule(TACCalibrator, MNumber); ++MNumber; - cout<<"Creating energy calibrator"<SetFileName(m_EcalFile); //EnergyCalibrator->EnablePreampTempCorrection(false); S->SetModule(EnergyCalibrator, MNumber); ++MNumber; - cout<<"Creating Event filter"<SetMaximumTotalEnergy(m_MaxEnergy*2); // Multiply by 2 because this is the event-level energy, i.e. sum over both sides S->SetModule(EventFilter, MNumber); ++MNumber; - - cout<<"Creating strip pairing"<SetModule(Pairing, MNumber); - cout<<"Initializing Loader"<Initialize() == false) return false; - cout<<"Initializing TAC calibrator"<Initialize() == false) return false; - cout<<"Initializing Energy calibrator"<Initialize() == false) return false; - cout<<"Initializing Event filter"<Initialize() == false) return false; - cout<<"Initializing Pairing"<Initialize() == false) return false; bool IsFinished = false; MReadOutAssembly* Event = new MReadOutAssembly(); // Pass Events through each module. Once calibrated, add the Event to the histograms - cout<<"Analyzing..."<Clear(); @@ -475,7 +463,7 @@ bool TrappingCorrectionAm241::Analyze() MStripHit* HVSH = GetDominantStrip(HVStrips, HVEnergyFraction); MStripHit* LVSH = GetDominantStrip(LVStrips, LVEnergyFraction); - if ((LVSH->HasCalibratedTiming()==true) && (HVSH->HasCalibratedTiming()==true)) { + if ((LVSH->HasCalibratedTiming()==true) && (HVSH->HasCalibratedTiming()==true)&& (LVSH != nullptr) && (HVSH != nullptr)) { double CTD = LVSH->GetTiming() - HVSH->GetTiming(); diff --git a/apps/TrappingCorrectionCs137.cxx b/apps/TrappingCorrectionCs137.cxx new file mode 100644 index 00000000..11f092ca --- /dev/null +++ b/apps/TrappingCorrectionCs137.cxx @@ -0,0 +1,830 @@ +/* + * TrappingCorrectionCs137.cxx + * + * + * Copyright (C) by Sophie Haight. + * All rights reserved. + * + * + * This code implementation is the intellectual property of + * Sophie Haight. + * + * By copying, distributing or modifying the Program (or any work + * based on the Program) you indicate your acceptance of this statement, + * and all its terms. + * + */ + +// Standard +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +using namespace std; + +// ROOT +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// MEGAlib +#include "MGlobal.h" +#include "MFile.h" +#include "MReadOutElementDoubleStrip.h" +#include "MFileReadOuts.h" +#include "MReadOutAssembly.h" +#include "MStripHit.h" +#include "MReadOutSequence.h" +#include "MSupervisor.h" +#include "MModuleLoaderMeasurementsHDF.h" +#include "MModuleEnergyCalibration.h" +#include "MModuleEventFilter.h" +#include "MModuleStripPairingMultiRoundChiSquare.h" +#include "MModuleStripPairingChiSquare.h" +#include "MModuleTACcut.h" +#include "MAssembly.h" + + +double g_MinCTD = -250; +double g_MaxCTD = 250; +int g_MinCounts = 1000; +int g_HVStrips = 64; +int g_LVStrips = 64; + +double g_CsPhotopeak = 661.7; + +const int NCTDBins = 10; +// We need NCTDBins + 1 edges to define the boundaries of NCTDBins +double g_CTDBinEdges[NCTDBins + 1]; + +// Run this initialization function ONCE at the start of your program (e.g., in main or class constructor) +void InitializeCTDBins() { + double center = (g_MaxCTD + g_MinCTD) / 2.0; + double halfWidth = (g_MaxCTD - g_MinCTD) / 2.0; + + for (int i = 0; i <= NCTDBins; ++i) { + // Map linear fraction from -1.0 (at i=0) to +1.0 (at i=NCTDBins) + double fraction = -1.0 + 2.0 * double(i) / double(NCTDBins); + + // Sinusoidal transformation: creates a higher density of points near the center + // If you prefer an even steeper density difference, you can use: pow(fraction, 3) + double nonLinearFraction = sin(fraction * M_PI / 2.0); + + // Calculate the actual CTD boundary coordinate + g_CTDBinEdges[i] = center + halfWidth * nonLinearFraction; + } +} + +int GetCTDBin(double CTD) { + // Hard bounds check + if (CTD < g_MinCTD || CTD >= g_MaxCTD) return -1; + + // Perform binary search to find the first edge that is strictly greater than our CTD value + auto it = std::upper_bound(g_CTDBinEdges, g_CTDBinEdges + NCTDBins + 1, CTD); + + // The bin index is simply the distance from the beginning boundary minus 1 + int bin = std::distance(g_CTDBinEdges, it) - 1; + + // Guard against edge cases at the absolute maximum limit + if (bin >= NCTDBins) bin = NCTDBins - 1; + if (bin < 0) bin = 0; + + return bin; +} + +//////////////////////////////////////////////////////////////////////////////// + + +//! A standalone program based on MEGAlib and ROOT +class TrappingCorrectionCs137 +{ +public: + //! Default constructor + TrappingCorrectionCs137(); + //! Default destructor + ~TrappingCorrectionCs137(); + + //! Parse the command line + bool ParseCommandLine(int argc, char** argv); + //! Analyze what ever needs to be analyzed... + bool Analyze(); + //! Interrupt the analysis + void Interrupt() { m_Interrupt = true; } + + //! Produce functions for fitting + + TF1* GeneratePhotopeakFunction(); + + MStripHit* GetDominantStrip(vector& Strips, double& EnergyFraction); + + private: + //! True, if the analysis needs to be interrupted + bool m_Interrupt; + //! The input file name + MString m_FileName; + MString m_EcalFile; + MString m_TACCalFile; + MString m_TACCutFile; + MString m_StripMapFile; + //! output file names + MString m_OutFile; + //! option to do a pixel-by-pixel calibration (instead of detector-by-detector) + bool m_PixelCorrect; + bool m_MultiRoundStripPairing; + bool m_ExcludeNN; + bool m_ContinueHDF5; + + double m_MinEnergy; + double m_MaxEnergy; + +}; + +//////////////////////////////////////////////////////////////////////////////// + + +//! Default constructor +TrappingCorrectionCs137::TrappingCorrectionCs137() : m_Interrupt(false) +{ + gStyle->SetPalette(1, 0); +} + + +//////////////////////////////////////////////////////////////////////////////// + + +//! Default destructor +TrappingCorrectionCs137::~TrappingCorrectionCs137() +{ + // Intentionally left blank +} + + +//////////////////////////////////////////////////////////////////////////////// + + +//! Parse the command line +bool TrappingCorrectionCs137::ParseCommandLine(int argc, char** argv) +{ + ostringstream Usage; + Usage<"< i+1) && (argv[i+1][0] != '-' || isalpha(argv[i+1][1]) == 0))){ + cout<<"Error: Option "<>> FullDetEndpoints; + + // Store the input files + vector FileNames; + FileNames.push_back(m_FileName); + cout << "file name stored" << endl; + + MString InputFile = FileNames[0]; + cout << " input file stored as: " << InputFile << endl; + vector HDFNames; + + TH1::SetDefaultSumw2(); + + // Detector-level maps organized by [CTDBin][DetID] + map> FullDetCTDHistograms; + map> FullDetHVEnergyHistograms; + map> FullDetLVEnergyHistograms; + + // Read in the input files and make a list of hdf5 files to calibrate + if ((InputFile.GetSubString(InputFile.Length() - 4)) == "hdf5") { + HDFNames.push_back(InputFile); + cout << "hdf names loaded correctly" << endl; + } else if ((InputFile.GetSubString(InputFile.Length() - 3)) == "txt") { + cout << "Reading input file " << InputFile << endl; + cout << "WARNING: When passing a list of files, ensure that you have chosen the correct HDF5 continuous reading mode. Use the --nocontinue option to suppress continuous file reading." << endl; + MFile F; + if (F.Open(InputFile) == false) { + cout << "Error: Failed to open input file." << endl; + } else { + MString Line; + while (F.ReadLine(Line)) { + MString Trimmed = Line.Trim(); + if ((Trimmed != "")) { + if (F.Exists(Trimmed) == true) { + HDFNames.push_back(Trimmed); + } else { + cout << "Error: Could not find file " << Trimmed << endl; + } + } + } + } + } else { + cout << "Error: Unrecognized file format: " << InputFile << endl; + } + + // Analyze all the data and fill in the histograms + for (unsigned int f = 0; f < HDFNames.size(); ++f) { + + MString File = HDFNames[f]; + cout << "Beginning analysis of file " << File << endl; + + // Create and initialize nuclearizer modules + MSupervisor* S = MSupervisor::GetSupervisor(); + + MModuleLoaderMeasurementsHDF* Loader; + MModuleTACcut* TACCalibrator; + MModuleEnergyCalibration* EnergyCalibrator; + MModuleEventFilter* EventFilter; + + unsigned int MNumber = 0; + cout << "Creating HDF5 loader" << endl; + Loader = new MModuleLoaderMeasurementsHDF(); + Loader->SetFileNameStripMap(m_StripMapFile); + Loader->SetFileName(File); + Loader->SetLoadContinuationFiles(m_ContinueHDF5); + S->SetModule(Loader, MNumber); + ++MNumber; + + cout << "Creating TAC calibrator" << endl; + TACCalibrator = new MModuleTACcut(); + TACCalibrator->SetTACCalFileName(m_TACCalFile); + TACCalibrator->SetTACCutFileName(m_TACCutFile); + S->SetModule(TACCalibrator, MNumber); + ++MNumber; + + cout << "Creating energy calibrator" << endl; + EnergyCalibrator = new MModuleEnergyCalibration(); + EnergyCalibrator->SetFileName(m_EcalFile); + S->SetModule(EnergyCalibrator, MNumber); + ++MNumber; + + cout << "Creating Event filter" << endl; + EventFilter = new MModuleEventFilter(); + EventFilter->SetMinimumLVStrips(1); + EventFilter->SetMaximumLVStrips(3); + EventFilter->SetMinimumHVStrips(1); + EventFilter->SetMaximumHVStrips(3); + EventFilter->SetMinimumHits(0); + EventFilter->SetMaximumHits(100); + EventFilter->SetMinimumTotalEnergy(m_MinEnergy); + EventFilter->SetMaximumTotalEnergy(m_MaxEnergy * 2); + S->SetModule(EventFilter, MNumber); + ++MNumber; + + cout << "Creating strip pairing" << endl; + MModule* Pairing; + if (m_MultiRoundStripPairing == true) { + Pairing = new MModuleStripPairingMultiRoundChiSquare(); + } else { + Pairing = new MModuleStripPairingChiSquare(); + } + S->SetModule(Pairing, MNumber); + + cout<<"Initializing Loader"<Initialize() == false) return false; + cout<<"Initializing TAC calibrator"<Initialize() == false) return false; + cout<<"Initializing Energy calibrator"<Initialize() == false) return false; + cout<<"Initializing Event filter"<Initialize() == false) return false; + cout<<"Initializing Pairing"<Initialize() == false) return false; + + bool IsFinished = false; + MReadOutAssembly* Event = new MReadOutAssembly(); + cout<<"Modules initialized, starting event loop"<Clear(); + if (Loader->IsReady()) { + + Loader->AnalyzeEvent(Event); + TACCalibrator->AnalyzeEvent(Event); + EnergyCalibrator->AnalyzeEvent(Event); + bool Unfiltered = EventFilter->AnalyzeEvent(Event); + + if (Unfiltered == true) { + + Pairing->AnalyzeEvent(Event); + + if ((Event->HasAnalysisProgress(MAssembly::c_StripPairing) == true) && (Unfiltered == true)) { + + for (unsigned int h = 0; h < Event->GetNHits(); ++h) { + double HVEnergy = 0.0; + double LVEnergy = 0.0; + vector HVStrips; + vector LVStrips; + + MHit* H = Event->GetHit(h); + int DetID = H->GetStripHit(0)->GetDetectorID(); + + for (unsigned int sh = 0; sh < H->GetNStripHits(); ++sh) { + MStripHit* SH = H->GetStripHit(sh); + + if ((m_ExcludeNN == false) || ((m_ExcludeNN == true) && (SH->IsNearestNeighbor() == false))) { + if (SH->IsLowVoltageStrip() == true) { + LVEnergy += SH->GetEnergy(); + LVStrips.push_back(SH); + } else { + HVEnergy += SH->GetEnergy(); + HVStrips.push_back(SH); + } + } + } + + if ((HVStrips.size() > 0) && (LVStrips.size() > 0)) { + + double HVEnergyFraction = 0; + double LVEnergyFraction = 0; + MStripHit* HVSH = GetDominantStrip(HVStrips, HVEnergyFraction); + MStripHit* LVSH = GetDominantStrip(LVStrips, LVEnergyFraction); + + if ((LVSH->HasCalibratedTiming() == true) && (HVSH->HasCalibratedTiming() == true) && (LVSH != nullptr) && (HVSH != nullptr)) { + + double CTD = LVSH->GetTiming() - HVSH->GetTiming(); + int CTDBin = GetCTDBin(CTD); + if (CTDBin < 0) continue; + + // TH1D* FullCTDHist = FullDetCTDHistograms[CTDBin][DetID]; + TH1D* FullHVHist = FullDetHVEnergyHistograms[CTDBin][DetID]; + TH1D* FullLVHist = FullDetLVEnergyHistograms[CTDBin][DetID]; + + if (FullHVHist == nullptr) { + char name[128]; sprintf(name, "HV_Detector%d_bin%d", DetID, CTDBin); + FullHVHist = new TH1D(name, name, (m_MaxEnergy - m_MinEnergy) * 2, m_MinEnergy, m_MaxEnergy); + FullDetHVEnergyHistograms[CTDBin][DetID] = FullHVHist; + } + if (FullLVHist == nullptr) { + char name[128]; sprintf(name, "LV_Detector%d_bin%d", DetID, CTDBin); + FullLVHist = new TH1D(name, name, (m_MaxEnergy - m_MinEnergy) * 2, m_MinEnergy, m_MaxEnergy); + FullDetLVEnergyHistograms[CTDBin][DetID] = FullLVHist; + } + + // FullCTDHist->Fill(CTD); + FullHVHist->Fill(HVEnergy); + FullLVHist->Fill(LVEnergy); + } + } + } + } + } + } + IsFinished = Loader->IsFinished(); + } + } + + // Place this outside/before your CTD bin and Detector loops! + ofstream MasterFitFile; + MasterFitFile.open(m_OutFile + MString("_All_CTDBin_FitResults.txt")); + MasterFitFile << "======================================================================" << endl; + MasterFitFile << "MASTER PHOTOPEAK FIT LOG FOR ALL CTD BINS AND DETECTORS" << endl; + MasterFitFile << "======================================================================" << endl << endl; + + // Do function fitting and recording for full detector outputs + for (int c = 0; c < NCTDBins; ++c) { + + cout << "Processing CTD bin " << c << endl; + + for (auto const& [DetID, FullHVHist] : FullDetHVEnergyHistograms[c]) { + + TH1D* HVHist = FullDetHVEnergyHistograms[c][DetID]; + TH1D* LVHist = FullDetLVEnergyHistograms[c][DetID]; + + if (HVHist->Integral() > g_MinCounts) { + + TF1* PhotopeakFunctionHV = GeneratePhotopeakFunction(); + TFitResultPtr HVFit = HVHist->Fit(PhotopeakFunctionHV, "S", "", 645, 675); + + TF1* PhotopeakFunctionLV = GeneratePhotopeakFunction(); + TFitResultPtr LVFit = LVHist->Fit(PhotopeakFunctionLV, "S", "", 645, 675); + + if ((HVFit >= 0)) { + + // Clear or initialize the vector for this specific bin and detector + FullDetEndpoints[c][DetID].clear(); + + // Parameter(2) is Mu for your CTD function, Parameter(1) is Mu for the Photopeaks + // FullDetEndpoints[c][DetID].push_back(CTDFit->Parameter(2)); // Index 0: CTD Centroid + FullDetEndpoints[c][DetID].push_back(HVFit->Parameter(1)); // Index 1: HV Photopeak Mu + FullDetEndpoints[c][DetID].push_back(HVFit->ParError(1)); // Index 2: HV Photopeak Mu Error + FullDetEndpoints[c][DetID].push_back(LVFit->Parameter(1)); // Index 2: LV Photopeak Mu + FullDetEndpoints[c][DetID].push_back(LVFit->ParError(1)); // Index 2: HV Photopeak Mu Error + + + MasterFitFile << "------------------------------------------------------------" << endl; + MasterFitFile << " DETECTOR ID: " << DetID << " | CTD BIN INDEX: " << c << endl; + MasterFitFile << "------------------------------------------------------------" << endl; + + // Redirect cout to our master file stream + std::streambuf* coutbuf = cout.rdbuf(); + cout.rdbuf(MasterFitFile.rdbuf()); + + // Passing "V" forces ROOT to print out all parameters, errors, and Chi2 configurations + HVFit->Print("V"); + + // Restore normal terminal routing + cout.rdbuf(coutbuf); + + MasterFitFile << endl << endl; // Add spacing between different bin entries + + // ofstream LVFitFile(DetID + MString("_CTDbin_") + c + MString("_LVEnergyFitResult_.txt")); + // coutbuf = cout.rdbuf(); + // cout.rdbuf(LVFitFile.rdbuf()); + // if (LVFit >= 0) { + // LVFit->Print(); + // } + // cout.rdbuf(coutbuf); + // LVFitFile.close(); + + TFile HVHistFile(m_OutFile + MString("_Det") + DetID + MString("_CTDbin_") + c + MString("_HVEnergyHist_Illum.root"), "recreate"); + TCanvas* HVHistCanvas = new TCanvas(); + HVHistCanvas->cd(); + HVHist->Draw("Hist"); + PhotopeakFunctionHV->Draw("same"); + HVHistCanvas->Write(); + HVHistFile.Close(); + + TFile LVHistFile(m_OutFile + MString("_Det") + DetID + MString("_CTDbin_") + c + MString("_LVEnergyHist_Illum.root"), "recreate"); + TCanvas* LVHistCanvas = new TCanvas(); + LVHistCanvas->cd(); + LVHist->Draw("Hist"); + PhotopeakFunctionLV->Draw("same"); + LVHistCanvas->Write(); + LVHistFile.Close(); + + } else { + cout << "Fits failed for CTD bin " << c << " Detector " << DetID << endl; + } + } else { + cout << "Fewer than " << g_MinCounts << " counts in CTD bin " << c << " Detector " << DetID << endl; + } + } + } + + // Place this at the absolute end of your Analyze() function + MasterFitFile.close(); + cout << "Master fit results log saved successfully." << endl; + + // Setup parameter file + ofstream OutputCalFile; + OutputCalFile.open(m_OutFile + MString("_parameters.txt")); + + // Updated header matching your request + OutputCalFile << "Det_ID" << '\t' + << "CTD_Bin" << '\t' + << "CTD_BinMidpoint_ns" << '\t' + << "HV_Centroid_keV" << '\t' + <<"HV_Centroid_error_keV" << '\t' + << "LV_Centroid_keV" << '\t' + <<"LV_Centroid_error_keV" << '\t' << endl; + cout << "Parameter file set up" << endl; + + // Loop systematically over each CTD bin first + for (int c = 0; c < NCTDBins; ++c) { + + cout << "Writing output tracking parameters for CTD bin: " << c << endl; + + // Calculate the exact midpoint of this specific variable CTD bin + double ctdBinMidpoint = (g_CTDBinEdges[c] + g_CTDBinEdges[c + 1]) / 2.0; + + // Loop over the detectors found inside this specific CTD bin + for (auto const& [DetID, FitsVec] : FullDetEndpoints[c]) { + + double HVCentroid = 0.0; + double HVCentroidError = 0.0; + double LVCentroid = 0.0; + double LVCentroidError = 0.0; + + // Ensure all 3 parameters (CTD mu, HV mu, LV mu) were successfully saved + if (FitsVec.size() >= 1) { + HVCentroid = FitsVec[0]; + HVCentroidError = FitsVec[1]; + LVCentroid = FitsVec[2]; + LVCentroidError = FitsVec[3]; + } + + // Write row entries corresponding purely to the detector level measurements + OutputCalFile << DetID << '\t' + << c << '\t' + << ctdBinMidpoint << '\t' + << HVCentroid << '\t' + << HVCentroidError << '\t' + << LVCentroid << '\t' + << LVCentroidError << '\t' << endl; + } + } + + OutputCalFile.close(); + cout << "Parameters file saved successfully." << endl; + watch.Stop(); + cout << "total time (s): " << watch.CpuTime() << endl; + + return true; +} +//////////////////////////////////////////////////////////////////////////////// + + +TF1* TrappingCorrectionCs137::GeneratePhotopeakFunction() +{ + // Component 1: Core Gaussian + // exp(-(x-x0)^2 / (2*sigma^2)) + MString gaussStr = "exp(-(x-[1])^2 / (2*[2]^2))"; + + // Component 2: Exponential Tail + Shelf + // BoverA * exp(gamma*(x-x0)) * 0.5 * erfc((x-x0)/(sigma*sigma_ratio*sqrt(2))) + MString expTailStr = "[3] * exp([4]*(x-[1])) * 0.5 * erfc((x-[1])/([2]*[5]*sqrt(2)))"; + + // Component 3: Linear Tail + Shelf + // BoverA * CoverB * (1 + D*(x-x0)) * 0.5 * erfc((x-x0)/(sigma*sigma_ratio*sqrt(2))) + MString linTailStr = "[3] * [6] * (1 + [7]*(x-[1])) * 0.5 * erfc((x-[1])/([2]*[5]*sqrt(2)))"; + + // Combine components with an overall normalization scaling factor [0] + MString fullFormula = "[0] * (" + gaussStr + " + " + expTailStr + " + " + linTailStr + ")"; + + // Instantiate TF1 over your expected fit window + TF1* PhotopeakFunction = new TF1("PhotopeakFunction", fullFormula.Data(), 645, 675); + + // Set Parameter Names + PhotopeakFunction->SetParName(0, "Amplitude"); + PhotopeakFunction->SetParName(1, "x0 (Mu)"); + PhotopeakFunction->SetParName(2, "Sigma Gauss"); + PhotopeakFunction->SetParName(3, "BoverA"); + PhotopeakFunction->SetParName(4, "Gamma"); + PhotopeakFunction->SetParName(5, "Sigma Ratio"); + PhotopeakFunction->SetParName(6, "CoverB"); + PhotopeakFunction->SetParName(7, "D (Lin Slope)"); + + // Provide initial sensible guesses for a Cs137 photopeak + PhotopeakFunction->SetParameter("Amplitude", 1000); + PhotopeakFunction->SetParameter("x0 (Mu)", 661.7); + PhotopeakFunction->SetParameter("Sigma Gauss", 2.0); + PhotopeakFunction->SetParameter("BoverA", 0.05); + PhotopeakFunction->SetParameter("Gamma", 0.5); + PhotopeakFunction->SetParameter("Sigma Ratio", 0.85); + PhotopeakFunction->SetParameter("CoverB", 0.13); + PhotopeakFunction->SetParameter("D (Lin Slope)", 0.028); + + // Set boundary limits to stabilize convergence + PhotopeakFunction->SetParLimits(0, 1, 1e8); + PhotopeakFunction->SetParLimits(1, 645, 675); // Keeps peak centered around 662 keV + PhotopeakFunction->SetParLimits(2, 0.5, 10); // Prevents sigma from blowing up or hitting zero + PhotopeakFunction->SetParLimits(3, 0.0, 1.0); // Tail shouldn't be larger than the main peak + PhotopeakFunction->SetParLimits(4, 0.001, 2.0); // Standard range for exponential decay factor + PhotopeakFunction->SetParLimits(5, 0.1, 5.0); // Ratio of shelf width to peak width + PhotopeakFunction->SetParLimits(6, 0.0, 5.0); + PhotopeakFunction->SetParLimits(7, -1.0, 1.0); + // // Gaussian with a low-E shelf + // TF1* PhotopeakFunction = new TF1("PhotopeakFunction", "gaus(0) + [0]*[3]*(1 - erf((x-[1])/(sqrt(2)*[2])))", 620, 680); + + // PhotopeakFunction->SetParName(0, "Gauss norm"); + // PhotopeakFunction->SetParName(1, "Mu"); + // PhotopeakFunction->SetParName(2, "Sigma"); + // PhotopeakFunction->SetParName(3, "Shelf norm"); + + // PhotopeakFunction->SetParameter("Gauss norm", 1000); + // PhotopeakFunction->SetParameter("Mu", 661.7); + // PhotopeakFunction->SetParameter("Sigma", 2); + // PhotopeakFunction->SetParameter("Shelf norm", 0.05); + + // PhotopeakFunction->SetParLimits(0, 10, 1e8); + // PhotopeakFunction->SetParLimits(1, 652, 672); + // PhotopeakFunction->SetParLimits(2, 1.0, 10); + // PhotopeakFunction->SetParLimits(3, 0, 0.1); + + return PhotopeakFunction; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +TrappingCorrectionCs137* g_Prg = 0; +int g_NInterruptCatches = 1; + +MStripHit* TrappingCorrectionCs137::GetDominantStrip(vector& Strips, double& EnergyFraction) +{ + double MaxEnergy = -numeric_limits::max(); + double TotalEnergy = 0.0; + MStripHit* MaxStrip = nullptr; + + // Iterate through strip hits and get the strip with highest energy + for (const auto SH : Strips) { + double Energy = SH->GetEnergy(); + TotalEnergy += Energy; + if (Energy > MaxEnergy) { + MaxStrip = SH; + MaxEnergy = Energy; + } + } + if (TotalEnergy == 0) { + EnergyFraction = 0; + } else { + EnergyFraction = MaxEnergy/TotalEnergy; + } + return MaxStrip; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +//! Called when an interrupt signal is flagged +//! All catched signals lead to a well defined exit of the program +void CatchSignal(int a) +{ + if (g_Prg != 0 && g_NInterruptCatches-- > 0) { + cout<<"Catched signal Ctrl-C (ID="<Interrupt(); + } else { + abort(); + } +} + + +//////////////////////////////////////////////////////////////////////////////// + + +//! Main program +int main(int argc, char** argv) +{ + // Catch a user interupt for graceful shutdown + signal(SIGINT, CatchSignal); + + // Initialize global MEGALIB variables, especially mgui, etc. + MGlobal::Initialize("Standalone", "a standalone example program"); + + TApplication TrappingCorrectionApp("TrappingCorrectionApp", 0, 0); + + InitializeCTDBins(); + + g_Prg = new TrappingCorrectionCs137(); + + if (g_Prg->ParseCommandLine(argc, argv) == false) { + cerr<<"Error during parsing of command line!"<Analyze() == false) { + cerr<<"Error during analysis!"< m_DataBufferLVInitial; diff --git a/include/MGUIExpoTrappingCorrection.h b/include/MGUIExpoTrappingCorrection.h new file mode 100644 index 00000000..023ef66c --- /dev/null +++ b/include/MGUIExpoTrappingCorrection.h @@ -0,0 +1,136 @@ +/* + * MGUIExpoTrappingCorrection.h + * + * Copyright (C) by Andreas Zoglauer. + * All rights reserved. + * + * Please see the source-file for the copyright-notice. + * + */ + + +#ifndef __MGUIExpoTrappingCorrection__ +#define __MGUIExpoTrappingCorrection__ + + +//////////////////////////////////////////////////////////////////////////////// + + +// ROOT libs: +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// MEGAlib libs: +#include "MGlobal.h" +#include "MGUIERBList.h" + +// NuSTAR libs +#include "MGUIExpo.h" + +// Forward declarations: + + +//////////////////////////////////////////////////////////////////////////////// + + +class MGUIExpoTrappingCorrection : public MGUIExpo +{ + // public Session: + public: + //! Default constructor + MGUIExpoTrappingCorrection(MModule* Module); + //! Default destructor + virtual ~MGUIExpoTrappingCorrection(); + + //! The creation part which gets overwritten + virtual void Create(); + + //! Update the frame + virtual void Update(); + + //! Reset the data in the UI + virtual void Reset(); + + //! Export the data in the UI + virtual void Export(const MString& FileName); + + //! Set the energy histogram parameters + void SetEnergyHistogramParameters(int NBins, double Min, double Max); + + //! Add data to the uncorrected energy histogram + void AddEnergyInitial(double Energy, bool IsNearestNeighbor, bool IsLV); + + //! Add data to the corrected energy histogram + void AddEnergyFinal(double Energy, bool IsNearestNeighbor, bool IsLV); + + // //! Add data to the energy histogram + // void AddEnergy(double Energy); + + // Callback slot for the Apply button + void OnApply(); + + // Getters for Finalize analysis + TH1D* GetEnergyHistogramLVInitial() const { return m_EnergyLVInitial; } + TH1D* GetEnergyHistogramLVFinal() const { return m_EnergyLVFinal; } + TH1D* GetEnergyHistogramHVInitial() const { return m_EnergyHVInitial; } + TH1D* GetEnergyHistogramHVFinal() const { return m_EnergyHVFinal; } + + // protected methods: + protected: + + + // protected members: + protected: + + // private members: + private: + // //! Energy canvas + // TRootEmbeddedCanvas* m_EnergyCanvas; + // //! Energy histogram + // TH1D* m_Energy; + + //! Energy Histograms + TH1D* m_EnergyLVInitial; + TH1D* m_EnergyLVFinal; + TH1D* m_EnergyHVInitial; + TH1D* m_EnergyHVFinal; + + //! Canvases & Display elements + TRootEmbeddedCanvas* m_CanvasLV; + TRootEmbeddedCanvas* m_CanvasHV; + TLegend* m_LegendLV; + TLegend* m_LegendHV; + + //! Interactive Control Widgets + TGNumberEntry* m_EntryNBins; + TGNumberEntry* m_EntryMinEnergy; + TGNumberEntry* m_EntryMaxEnergy; + TGTextButton* m_ButtonApply; + + TGCheckButton* m_CheckLogY; + + +#ifdef ___CLING___ + public: + ClassDef(MGUIExpoTrappingCorrection, 1) // basic class for dialog windows +#endif + +}; + +#endif + + +//////////////////////////////////////////////////////////////////////////////// diff --git a/include/MGUIOptionsTrappingCorrection.h b/include/MGUIOptionsTrappingCorrection.h new file mode 100644 index 00000000..e5ff733c --- /dev/null +++ b/include/MGUIOptionsTrappingCorrection.h @@ -0,0 +1,87 @@ +/* + * MGUIOptionsTrappingCorrection.h + * + * Copyright (C) by Andreas Zoglauer. + * All rights reserved. + * + * Please see the source-file for the copyright-notice.:q + * + */ + + +#ifndef __MGUIOptionsTrappingCorrection__ +#define __MGUIOptionsTrappingCorrection__ + + +//////////////////////////////////////////////////////////////////////////////// + + +// ROOT libs: +#include +#include +#include +#include +#include +#include +#include +#include + +// MEGAlib libs: +#include "MGlobal.h" +#include "MGUIEFileSelector.h" +#include "MGUIOptions.h" +#include "MGUIERBList.h" +#include "MGUIEEntry.h" + +// Nuclearizer libs: +#include "MModule.h" + + +// Forward declarations: + + +//////////////////////////////////////////////////////////////////////////////// + + +//! The user interface for the universal energy calibration +class MGUIOptionsTrappingCorrection : public MGUIOptions +{ + // public Session: + public: + //! Default constructor + MGUIOptionsTrappingCorrection(MModule* Module); + //! Default destructor + virtual ~MGUIOptionsTrappingCorrection(); + + //! The creation part which gets overwritten + virtual void Create(); + + // protected methods: + protected: + + //! Actions after the Apply or OK button has been pressed + virtual bool OnApply(); + + + + // protected members: + protected: + + // private members: + private: + + //! Select which file to load + MGUIEFileSelector* m_SimCCEFileSelector; + + +#ifdef ___CLING___ + public: + ClassDef(MGUIOptionsTrappingCorrection, 1) +#endif + +}; + +#endif + + +//////////////////////////////////////////////////////////////////////////////// diff --git a/include/MHit.h b/include/MHit.h index a44ac6b6..1dd21e38 100644 --- a/include/MHit.h +++ b/include/MHit.h @@ -73,6 +73,14 @@ class MHit //! Return the position resolution of the hit MVector GetPositionResolution() const { return m_PositionResolution; } + // Local Position: + + //! Set the local/raw position of the hit + void SetLocalPosition(const MVector& Position) { m_LocalPosition = Position; } + + //! Return the local/raw position of the hit + MVector GetLocalPosition() const { return m_LocalPosition; } + // Energy: @@ -178,6 +186,10 @@ class MHit //! Position resolution of the hit MVector m_PositionResolution; + + //! Position of the hit in local detector coordinates + MVector m_LocalPosition; + //! Energy of the hit double m_Energy; //! Energy resolution of the hit diff --git a/include/MModuleDepthCalibration.h b/include/MModuleDepthCalibration.h index 0926a61f..196678b0 100644 --- a/include/MModuleDepthCalibration.h +++ b/include/MModuleDepthCalibration.h @@ -119,13 +119,20 @@ class MModuleDepthCalibration : public MModule //! Create the XML configuration MXmlNode* CreateXmlConfiguration(); + //! Returns the strip with most energy from vector Strips, also gives back the energy fraction + MStripHit* GetDominantStrip(std::vector& Strips, double& EnergyFraction); + + + //! Determine the Grade (geometry of charge sharing) of the Hit + int GetHitGrade(MHit* H); + + const std::vector& GetDetectorIDs() const { return m_DetectorIDs; } + //! Finalize void Finalize(); // protected methods: protected: - //! Returns the strip with most energy from vector Strips, also gives back the energy fraction - MStripHit* GetDominantStrip(std::vector& Strips, double& EnergyFraction); //! Retrieve the appropriate Depth values given the DetID vector GetDepth(int DetID); @@ -141,9 +148,6 @@ class MModuleDepthCalibration : public MModule //! Adds a Depth-to-CTD relation bool AddDepthCTD(vector Depth, vector> CTDArr, int DetID, unordered_map>& DepthGrid, unordered_map>>& CTDMap, unordered_map>& SplineMap, unsigned int NPoints); - //! Determine the Grade (geometry of charge sharing) of the Hit - int GetHitGrade(MHit* H); - //! Return the coefficients for a pixel vector* GetPixelCoeffs(int PixelCode); diff --git a/include/MModuleTrappingCorrection.h b/include/MModuleTrappingCorrection.h new file mode 100644 index 00000000..b7ad1271 --- /dev/null +++ b/include/MModuleTrappingCorrection.h @@ -0,0 +1,162 @@ +/* + * MModuleTrappingCorrection.h + * + * Copyright (C) 2008-2008 by Andreas Zoglauer. + * All rights reserved. + * + * Please see the source-file for the copyright-notice. + * + */ + + +#ifndef __MModuleTrappingCorrection__ +#define __MModuleTrappingCorrection__ + + +//////////////////////////////////////////////////////////////////////////////// + + +// Standard libs: +#include +#include +#include +#include + +// ROOT libs: + +// MEGAlib libs: +#include "MGlobal.h" +#include "MModule.h" +#include "MGUIEEntry.h" + + + +// Nuclearizer libs: +#include "MGUIExpoPlotSpectrum.h" +#include "MModuleEnergyCalibration.h" +#include "MGUIExpoTrappingCorrection.h" +#include "MGUIOptionsTrappingCorrection.h" +#include "MModuleDepthCalibration.h" + +// Forward declarations: + + +//////////////////////////////////////////////////////////////////////////////// + + +class MModuleTrappingCorrection : public MModule +{ + // public interface: + public: + //! Default constructor + MModuleTrappingCorrection(); + //! Default destructor + virtual ~MModuleTrappingCorrection(); + + //! Create a new object of this class + virtual MModuleTrappingCorrection* Clone() { return new MModuleTrappingCorrection(); } + + //! Initialize the module + virtual bool Initialize(); + + //! Create the expos + virtual void CreateExpos(); + + //! Main data analysis routine, which updates the event to a new level + virtual bool AnalyzeEvent(MReadOutAssembly* Event); + + //! Show the options GUI + virtual void ShowOptionsGUI(); + + //! Set filename for SimCCE file + void SetSimCCEFileName( const MString& FileName) { m_SimCCEFileName = FileName; } + + //! Get filename for SimCCE file + MString GetSimCCEFileName() const { return m_SimCCEFileName; } + + //! Finalize the module + virtual void Finalize(); + + //! Read the XML configuration + bool ReadXmlConfiguration(MXmlNode* Node); + + //! Create the XML configuration + MXmlNode* CreateXmlConfiguration(); + + // Getters to retrieve the calculated values after Finalize() runs + double GetDirectFWHM_LV() const { return m_DirectFWHM_LV; } + double GetDirectFWHM_HV() const { return m_DirectFWHM_HV; } + + + // protected methods: + protected: + + //! Load in the specified SimCCE file + bool LoadSimCCEFile(MString FName); + + //! Get the Sim-based corrected energy given the CTD value, uncorrected energy, and the sorted Sim CCE values + double GetSimBasedCorrectedEnergy(double ctd_val, double uncorrected_energy, const std::vector& sim_cce_sorted_e, const std::vector& sim_cce_sorted_h, double paramA, double paramB, double paramC); + + //! Interpolate a value given x, xp, and fp + double Interpolate(double x, const std::vector& xp, const std::vector& fp); + + + + // private methods + private: + + + // protected members: + protected: + + double m_SimCCE_Energy; + MString m_SimCCEFileName; + + // unordered_map m_Detectors; + vector m_DetectorIDs; + MModuleEnergyCalibration* m_EnergyCalibration; + MGUIExpoTrappingCorrection* m_ExpoTrappingCorrection; + + bool m_SimCCEFileIsLoaded; + + double m_ParamA_HV; + double m_ParamA_LV; + double m_ParamB; + double m_ParamC; + std::vector m_Depths; + std::vector m_CCEs_HV_e; + std::vector m_CCEs_HV_h; + std::vector m_CCEs_LV_e; + std::vector m_CCEs_LV_h; + + + + // private members: + private: + + MModuleDepthCalibration* m_DepthCalibration = nullptr; + + //! Updated GUI to display the energy histogram + MGUIExpoTrappingCorrection* m_ExpoSpectrum; + + TF1* GeneratePhotopeakFunction(); + + double CalculateDirectFWHM(TH1D* hist); + + double m_DirectFWHM_LV = 0.0; + double m_DirectFWHM_HV = 0.0; + + TGCheckButton* m_LogYButton; + + +#ifdef ___CLING___ + public: + ClassDef(MModuleTrappingCorrection, 0) // no description +#endif + +}; + +#endif + + +//////////////////////////////////////////////////////////////////////////////// diff --git a/resource/dee/dummy_trapping_parameters_singledetector.csv b/resource/dee/dummy_trapping_parameters_singledetector.csv new file mode 100644 index 00000000..a4ff0a52 --- /dev/null +++ b/resource/dee/dummy_trapping_parameters_singledetector.csv @@ -0,0 +1,154 @@ +# A_HV,A_LV,B,C +1.0009264205609762,1.00084829470918,1.288689028046797,0.3010529785474152 +# z_depth_cm,e_CCE_HV,h_CCE_HV,e_CCE_LV,h_CCE_LV +-0.75,0.9960754005722672,1.0,0.9995680702824706,1.0 +-0.74,0.9961283655288202,1.0,0.9995997923338564,0.9999868765995326 +-0.73,0.9961799143100644,1.0,0.9996279267722878,0.9999704487457096 +-0.72,0.9962307077978106,1.0,0.999652945652101,0.999951566155109 +-0.71,0.9962799906985604,1.0,0.9996737160428892,0.9999320226738372 +-0.7,0.9963278574240018,1.0,0.9996934479141384,0.9999112518241762 +-0.6900000000000001,0.9963751576793192,1.0,0.9997102530485034,0.9998898200838444 +-0.6799999999999999,0.9964219858762003,0.999999905589124,0.9997259252524616,0.9998655559549224 +-0.67,0.996467775544522,0.999999905589124,0.9997390483630052,0.9998428024332484 +-0.6599999999999999,0.9965127155076576,0.999999905589124,0.9997533988148228,0.9998179718266084 +-0.65,0.996555200766925,0.9999994335347432,0.9997637840102168,0.9997949350660752 +-0.64,0.9965989133781268,0.9999994335347432,0.9997751133142831,0.9997692547428584 +-0.63,0.9966405489322088,0.9999993391238672,0.999784932044474,0.9997448962009834 +-0.62,0.9966822788979778,0.9999991503021148,0.9997933346116564,0.9997195935295784 +-0.61,0.9967226871001252,0.9999991503021148,0.9998014539462374,0.9996941020322672 +-0.6,0.9967628120672112,0.9999990558912386,0.9998102341568889,0.99966936583858 +-0.5900000000000001,0.9968020873291114,0.999997167673716,0.9998177870262664,0.9996436855153632 +-0.58,0.9968404184741392,0.999997167673716,0.9998235460891668,0.999618099605099 +-0.5700000000000001,0.9968790328542289,0.999997167673716,0.9998298716172704,0.9995930801725528 +-0.5599999999999999,0.9969165142940712,0.99999707326284,0.9998377077192496,0.999567588675242 +-0.55,0.996954184557288,0.99999707326284,0.999842239440876,0.9995420027649778 +-0.54,0.996991477173756,0.9999966012084592,0.999848187325511,0.9995169833324316 +-0.53,0.9970278256733512,0.999996506797583,0.9998520581710671,0.9994913974221676 +-0.52,0.9970642685846338,0.9999963179758308,0.9998577228231004,0.9994659059248566 +-0.51,0.997099484143982,0.9999962235649548,0.9998614992577888,0.9994411697311696 +-0.5,0.997133944409832,0.999994240936556,0.9998669750880878,0.9994155838209052 +-0.49000000000000005,0.9971688767341186,0.999994240936556,0.9998706571119091,0.9993901867365474 +-0.48,0.9972029593532196,0.9999941465256798,0.9998763217639424,0.9993659226076254 +-0.47000000000000003,0.9972367587372588,0.9999936744712992,0.9998783988030212,0.9993401478714552 +-0.45999999999999996,0.99727084135636,0.9999936744712992,0.9998830249355152,0.9993145619611912 +-0.45,0.99730322456509,0.9999933912386708,0.999886046083266,0.9992912419617992 +-0.44000000000000006,0.9973361742439442,0.9999932968277946,0.9998896336962204,0.9992652783997232 +-0.43,0.9973690295111108,0.999991408610272,0.9998926548439714,0.99924129750966 +-0.42000000000000004,0.9974000909562196,0.999991408610272,0.9998956759917224,0.9992161836641612 +-0.41,0.9974331350467608,0.9999913141993956,0.9998986971394734,0.999191069818662 +-0.4,0.997464385315244,0.9999908421450152,0.999902001519826,0.999167938645176 +-0.39,0.9974948802902291,0.9999905589123866,0.99990492825671,0.999142447147865 +-0.38,0.9975268858522104,0.9999905589123866,0.9999059667762494,0.999118843909614 +-0.37,0.9975576640622574,0.9999904645015104,0.999908799102266,0.999093918890021 +-0.36,0.9975877813904936,0.9999903700906344,0.9999116314282824,0.999069843587005 +-0.35,0.997617332248606,0.9999898980362538,0.9999145581651664,0.999045296219224 +-0.33999999999999997,0.997646977518406,0.9999896148036254,0.9999155022738384,0.9990207488514428 +-0.32999999999999996,0.997677566905078,0.9999896148036254,0.9999183345998552,0.9989975232650038 +-0.32,0.9977074009982526,0.9999876321752268,0.9999211669258716,0.9989727870713168 +-0.31,0.9977362909745544,0.9999876321752268,0.9999239992518882,0.99894993913669 +-0.3,0.9977654641859178,0.9999870657099698,0.999924754538826,0.9989252029430028 +-0.29,0.9977938821037832,0.9999867824773414,0.9999273036322408,0.998902543834282 +-0.27999999999999997,0.997822394433336,0.9999867824773414,0.9999299471365232,0.998878090879454 +-0.27,0.997850717939514,0.9999847054380664,0.9999326850516724,0.9988551485318738 +-0.26,0.997878947034005,0.9999841389728096,0.9999335347494774,0.9988315452936228 +-0.25,0.997907176128496,0.9999838557401812,0.9999343844472824,0.99880850853309 +-0.24,0.9979341778710524,0.9999838557401812,0.9999371223624318,0.998785377359604 +-0.22999999999999998,0.9979613684369836,0.9999818731117824,0.9999394826341124,0.9987617741213528 +-0.22000000000000003,0.9979886534146016,0.999981401057402,0.999941087618855,0.998739775903303 +-0.21000000000000002,0.9980156551571582,0.9999811178247736,0.9999429758361996,0.998715606187334 +-0.2,0.9980431289581512,0.9999810234138974,0.9999434478905356,0.9986931359045192 +-0.19,0.9980701307007076,0.9999809290030212,0.9999459025730832,0.9986692494274092 +-0.18,0.9980973212666384,0.9999804569486404,0.9999468466817554,0.9986469679705 +-0.16999999999999998,0.9981228124221988,0.999980079305136,0.9999491125425686,0.9986243088617792 +-0.16,0.9981487756361956,0.9999781910876132,0.9999516616359836,0.9986010832753404 +-0.15,0.9981753053203158,0.9999780966767372,0.9999525113337884,0.9985788962313844 +-0.13999999999999999,0.9982020238278104,0.9999774358006044,0.9999527945663902,0.9985551985801804 +-0.13,0.998226476454811,0.999977246978852,0.9999534554424608,0.9985332947750832 +-0.12,0.9982527229038692,0.9999751699395772,0.9999560989467428,0.9985111077311274 +-0.11000000000000001,0.9982785917061788,0.9999746034743202,0.9999563821793446,0.9984875989058294 +-0.1,0.9983034219799284,0.999974320241692,0.9999590256836268,0.9984655062748266 +-0.09,0.9983288187238014,0.9999724320241692,0.9999611027227056,0.9984432248179176 +-0.08,0.9983543098793618,0.9999718655589124,0.9999618580096432,0.9984207545351028 +-0.06999999999999999,0.9983788569180494,0.999971487915408,0.9999620468313778,0.9983985674911467 +-0.06,0.9984038760151736,0.999971487915408,0.9999628021183156,0.9983752474917548 +-0.05,0.9984295559941084,0.9999709214501512,0.99996469033566,0.9983531548607518 +-0.04,0.9984528756808618,0.9999705438066464,0.9999656344443322,0.9983317231204198 +-0.03,0.998477989189673,0.9999686555891238,0.9999658232660666,0.9983083087080749 +-0.02,0.9985017809348627,0.9999680891238673,0.9999679947160126,0.998286593728884 +-0.01,0.9985266112086122,0.9999677114803625,0.999968655592083,0.9982652564015052 +0.0,0.9985504029538018,0.9999657288519638,0.9999706382202947,0.9982417475762072 +0.01,0.9985743835223662,0.9999651623867069,0.9999712990963652,0.9982202214229224 +0.02,0.9985982696792428,0.9999638406344412,0.9999714879180998,0.9981987896825903 +0.03,0.9986223446594944,0.9999628021148036,0.9999717711507012,0.9981756585091044 +0.04,0.9986454755228732,0.9999620468277944,0.9999723376159048,0.9981542267687724 +0.05,0.9986693616797502,0.9999619524169184,0.999973376135444,0.9981327950284404 +0.06,0.9986927757781908,0.9999613859516616,0.999974603476718,0.9981107968103904 +0.06999999999999999,0.9987164731116932,0.999960158610272,0.9999751699419211,0.9980885153534818 +0.08,0.9987380933880756,0.9999589312688822,0.9999753587636556,0.9980668947872436 +0.09,0.9987626404267635,0.9999582703927492,0.9999774358027346,0.9980461239375829 +0.1,0.998784921584957,0.9999561933534744,0.999978002267938,0.998022709525238 +0.11000000000000001,0.9988089965652084,0.9999553436555892,0.9999790407874772,0.998001561023765 +0.12,0.9988307112532784,0.9999533610271902,0.9999802681287512,0.9979806013481982 +0.13,0.9988522371179738,0.9999525113293052,0.9999808345939544,0.9979586975431012 +0.13999999999999999,0.9988764065099124,0.9999524169184292,0.999981023415689,0.9979365104991452 +0.15,0.9988983100213568,0.99995166163142,0.999981117826556,0.9979156452365314 +0.16,0.9989203079444886,0.999949584592145,0.9999814010591578,0.9978946855609646 +0.16999999999999998,0.9989429667494312,0.999948829305136,0.999981967524361,0.9978725929299616 +0.18,0.9989647758491884,0.9999466578549848,0.999983666919971,0.9978508779507708 +0.19,0.998986018478822,0.9999458081570995,0.99998404456344,0.9978296350363448 +0.2,0.999009243753888,0.999943825528701,0.9999846110286432,0.9978090530125898 +0.21000000000000002,0.9990302975601472,0.9999429758308158,0.9999847998503776,0.9977862050779628 +0.22000000000000003,0.9990513513664064,0.9999424093655588,0.9999856495481826,0.997765434228302 +0.22999999999999998,0.999074104583036,0.9999401435045318,0.9999868768894564,0.9977440969009232 +0.24,0.9990952528009824,0.9999394826283988,0.9999874433546598,0.9977237037030744 +0.25,0.9991166842539904,0.9999372167673716,0.9999876321763944,0.9977020831368364 +0.26,0.9991369827667516,0.9999363670694864,0.9999877265872616,0.9976801793317396 +0.27,0.9991602080418176,0.9999338179758308,0.9999896148046058,0.9976594084820788 +0.27999999999999997,0.9991806009662658,0.9999334403323265,0.9999898036263404,0.9976386376324176 +0.29,0.999201277125776,0.9999325906344412,0.9999902756806764,0.997617772369804 +0.3,0.9992220476969732,0.9999306080060424,0.9999904645024108,0.997595207674036 +0.31,0.9992421573863596,0.9999296638972812,0.999990558913278,0.9975750977150462 +0.32,0.999265099426364,0.999926831570997,0.999990558913278,0.9975539492135732 +0.32999999999999996,0.9992851147040632,0.999924848942598,0.9999907477350124,0.9975335560157244 +0.33999999999999997,0.9993056020401988,0.999923999244713,0.9999912197893484,0.9975122186883454 +0.35,0.9993258061412724,0.9999231495468278,0.9999913142002156,0.9974904092962016 +0.36,0.9993458214189718,0.9999208836858008,0.999991408611083,0.9974698272724468 +0.37,0.9993659311083584,0.9999185234138972,0.9999915030219504,0.9974491508357388 +0.38,0.999387551384741,0.9999172960725076,0.9999933912392946,0.9974290408767488 +0.39,0.9994081331325636,0.999914746978852,0.9999935800610292,0.9974076091364172 +0.4,0.9994276763518268,0.9999143693353476,0.9999935800610292,0.997385799744273 +0.41,0.9994475972178388,0.9999116314199396,0.9999941465262324,0.9973654065464244 +0.42000000000000004,0.9994672348487889,0.9999106873111784,0.9999941465262324,0.9973447301097164 +0.43,0.9994872501264882,0.9999078549848942,0.9999943353479668,0.9973246201507264 +0.44000000000000006,0.999506415699002,0.9999052114803626,0.9999943353479668,0.9973034716492536 +0.45,0.9995258645065777,0.9999043617824774,0.9999962235653111,0.9972815678441568 +0.45999999999999996,0.9995450300790916,0.9999015294561934,0.9999963179761784,0.997261741124026 +0.47000000000000003,0.999564101239918,0.9998989803625378,0.9999964123870456,0.997241159100271 +0.48,0.9995831724007448,0.999896336858006,0.9999965067979129,0.9972206714894692 +0.49000000000000005,0.9996022435615716,0.9998952983383684,0.999996978852249,0.9972009391822914 +0.5,0.9996211258990236,0.9998923716012084,0.999996978852249,0.9971795074419594 +0.51,0.9996401026481628,0.9998890672205438,0.999996978852249,0.9971578868757216 +0.52,0.999658984985615,0.9998860460725076,0.9999971676739834,0.9971380601555908 +0.53,0.9996773952646308,0.999885101963746,0.9999971676739834,0.9971175725447888 +0.54,0.9996957111319592,0.9998821752265862,0.9999981117826556,0.997096990521034 +0.55,0.999714121410975,0.9998773602719034,0.9999990558913278,0.997077258213856 +0.5599999999999999,0.9997316819848056,0.9998757552870092,0.999999150302195,0.9970560152994304 +0.5700000000000001,0.9997500922638216,0.9998726397280968,0.9999992447130622,0.9970343947331924 +0.58,0.9997674640142772,0.9998674471299092,0.9999992447130622,0.9970146624260146 +0.5900000000000001,0.9997830419426754,0.9998658421450152,0.9999993391239296,0.9969942692281656 +0.6,0.9998008857515676,0.9998605551359516,0.9999993391239296,0.99697397044327 +0.61,0.9998173133851508,0.9998573451661632,0.9999998111782656,0.996954049310186 +0.62,0.9998351571940431,0.999853663141994,0.9999998111782656,0.9969334672864312 +0.63,0.9998495077705069,0.999847998489426,0.9999998111782656,0.9969120355460992 +0.64,0.9998661242274646,0.9998442220543808,0.9999998111782656,0.9968916423482505 +0.65,0.9998827406844224,0.9998384629909366,0.9999999055891328,0.9968720044540256 +0.6599999999999999,0.9998965247907624,0.999834592145015,1.0,0.9968516112561766 +0.67,0.9999117250724114,0.999828077794562,1.0,0.9968314068842338 +0.6799999999999999,0.9999250371203152,0.9998212802114804,1.0,0.996811391338197 +0.6900000000000001,0.9999396709318406,0.9998156155589124,1.0,0.996791186966254 +0.7,0.9999515668044354,0.9998088179758308,1.0,0.996769755225922 +0.71,0.9999628962069066,0.9998000377643504,1.0,0.9967493620280732 +0.72,0.9999747920795016,0.9997908799093655,1.0,0.9967296297208956 +0.73,0.9999852717767874,0.99978125,1.0,0.9967094253489528 +0.74,0.9999932967702044,0.9997709592145017,1.0,0.9966895042158688 +0.75,1.0,0.9997587802114803,1.0,0.996669488669832 diff --git a/src/MAssembly.cxx b/src/MAssembly.cxx index e6b43581..499a49e5 100644 --- a/src/MAssembly.cxx +++ b/src/MAssembly.cxx @@ -66,6 +66,7 @@ using namespace std; #include "MModuleLoaderMeasurementsL0.h" #include "MModuleEnergyCalibration.h" #include "MModuleDepthCalibration.h" +#include "MModuleTrappingCorrection.h" #include "MModuleStripPairingMultiRoundChiSquare.h" #include "MModuleStripPairingChiSquare.h" #include "MModuleEventFilter.h" @@ -135,6 +136,7 @@ MAssembly::MAssembly() m_Supervisor->AddAvailableModule(new MModuleStripPairingMultiRoundChiSquare()); m_Supervisor->AddAvailableModule(new MModuleStripPairingChiSquare()); m_Supervisor->AddAvailableModule(new MModuleDepthCalibration()); + m_Supervisor->AddAvailableModule(new MModuleTrappingCorrection()); m_Supervisor->AddAvailableModule(new MModuleEventSaver()); m_Supervisor->AddAvailableModule(new MModuleSaverMeasurementsL0()); diff --git a/src/MGUIExpoTrappingCorrection.cxx b/src/MGUIExpoTrappingCorrection.cxx new file mode 100644 index 00000000..a935a0e1 --- /dev/null +++ b/src/MGUIExpoTrappingCorrection.cxx @@ -0,0 +1,358 @@ +/* + * MGUIExpoTrappingCorrection.cxx + * + * + * Copyright (C) by Andreas Zoglauer + * All rights reserved. + * + * + * This code implementation is the intellectual property of + * Andreas Zoglauer. + * + * By copying, distributing or modifying the Program (or any work + * based on the Program) you indicate your acceptance of this statement, + * and all its terms. + * + */ + + +// Include the header: +#include "MGUIExpoTrappingCorrection.h" + +// Standard libs: + +// ROOT libs: +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// MEGAlib libs: +#include "MStreams.h" + + + +//////////////////////////////////////////////////////////////////////////////// + + +#ifdef ___CLING___ +ClassImp(MGUIExpoTrappingCorrection) +#endif + +//////////////////////////////////////////////////////////////////////////////// + +MGUIExpoTrappingCorrection::MGUIExpoTrappingCorrection(MModule* Module) : MGUIExpo(Module) +{ + m_TabTitle = "Trapping Correction"; + + // LV Histograms + m_EnergyLVInitial = new TH1D("EnergyLVInitial", "LV Spectrum (Uncorrected vs Corrected)", 200, 620, 700); + m_EnergyLVInitial->SetXTitle("Energy [keV]"); + m_EnergyLVInitial->SetYTitle("Counts"); + m_EnergyLVInitial->GetYaxis()->SetNoExponent(kTRUE); + m_EnergyLVInitial->SetLineColor(kGray+2); + m_EnergyLVInitial->SetLineWidth(2); + m_EnergyLVInitial->SetLineStyle(2); + + m_EnergyLVFinal = new TH1D("EnergyLVFinal", "LV Spectrum (Uncorrected vs Corrected)", 200, 620, 700); + m_EnergyLVFinal->SetXTitle("Energy [keV]"); + m_EnergyLVFinal->SetYTitle("Counts"); + m_EnergyLVFinal->GetYaxis()->SetNoExponent(kTRUE); + m_EnergyLVFinal->SetLineColor(kAzure+1); + m_EnergyLVFinal->SetLineWidth(2); + m_EnergyLVFinal->SetFillColorAlpha(kAzure-9, 0.35); + + // HV Histograms + m_EnergyHVInitial = new TH1D("EnergyHVInitial", "HV Spectrum (Uncorrected vs Corrected)", 200, 620, 700); + m_EnergyHVInitial->SetXTitle("Energy [keV]"); + m_EnergyHVInitial->SetYTitle("Counts"); + m_EnergyHVInitial->GetYaxis()->SetNoExponent(kTRUE); + m_EnergyHVInitial->SetLineColor(kGray+2); + m_EnergyHVInitial->SetLineWidth(2); + m_EnergyHVInitial->SetLineStyle(2); + + m_EnergyHVFinal = new TH1D("EnergyHVFinal", "HV Spectrum (Uncorrected vs Corrected)", 200, 620, 700); + m_EnergyHVFinal->SetXTitle("Energy [keV]"); + m_EnergyHVFinal->SetYTitle("Counts"); + m_EnergyHVFinal->GetYaxis()->SetNoExponent(kTRUE); + m_EnergyHVFinal->SetLineColor(kOrange+7); + m_EnergyHVFinal->SetLineWidth(2); + m_EnergyHVFinal->SetFillColorAlpha(kOrange-9, 0.35); + + // Initialize canvases and buttons + m_CanvasLV = nullptr; + m_CanvasHV = nullptr; + m_LegendLV = nullptr; + m_LegendHV = nullptr; + + m_EntryNBins = nullptr; + m_EntryMinEnergy = nullptr; + m_EntryMaxEnergy = nullptr; + m_CheckLogY = nullptr; + m_ButtonApply = nullptr; + + SetCleanup(kDeepCleanup); +} + +//////////////////////////////////////////////////////////////////////////////// + +MGUIExpoTrappingCorrection::~MGUIExpoTrappingCorrection() +{ + // Memory cleanup handled by ROOT's kDeepCleanup +} + +//////////////////////////////////////////////////////////////////////////////// + +void MGUIExpoTrappingCorrection::Create() +{ + if (m_IsCreated == true) return; + + m_Mutex.Lock(); + + // Create labels and buttons on GUI + // Top frame + TGHorizontalFrame* ControlFrame = new TGHorizontalFrame(this, 800, 30); + + // Binning entry + TGLabel* LabelBins = new TGLabel(ControlFrame, "Bins:"); + ControlFrame->AddFrame(LabelBins, new TGLayoutHints(kLHintsLeft | kLHintsCenterY, 5, 2, 2, 2)); + + m_EntryNBins = new TGNumberEntry(ControlFrame, 200, 5, -1, TGNumberFormat::kNESInteger, TGNumberFormat::kNEAPositive); + ControlFrame->AddFrame(m_EntryNBins, new TGLayoutHints(kLHintsLeft | kLHintsCenterY, 2, 10, 2, 2)); + + // Min Energy entry + TGLabel* LabelMin = new TGLabel(ControlFrame, "Min Energy [keV]:"); + ControlFrame->AddFrame(LabelMin, new TGLayoutHints(kLHintsLeft | kLHintsCenterY, 5, 2, 2, 2)); + + m_EntryMinEnergy = new TGNumberEntry(ControlFrame, 0, 6, -1, TGNumberFormat::kNESRealOne); + ControlFrame->AddFrame(m_EntryMinEnergy, new TGLayoutHints(kLHintsLeft | kLHintsCenterY, 2, 10, 2, 2)); + + // Max Energy entry + TGLabel* LabelMax = new TGLabel(ControlFrame, "Max Energy [keV]:"); + ControlFrame->AddFrame(LabelMax, new TGLayoutHints(kLHintsLeft | kLHintsCenterY, 5, 2, 2, 2)); + + m_EntryMaxEnergy = new TGNumberEntry(ControlFrame, 1000, 6, -1, TGNumberFormat::kNESRealOne); + ControlFrame->AddFrame(m_EntryMaxEnergy, new TGLayoutHints(kLHintsLeft | kLHintsCenterY, 2, 10, 2, 2)); + + // Log Y Checkbox + m_CheckLogY = new TGCheckButton(ControlFrame, "Log Y Scale"); + ControlFrame->AddFrame(m_CheckLogY, new TGLayoutHints(kLHintsLeft | kLHintsCenterY, 10, 10, 2, 2)); + m_CheckLogY->Connect("Clicked()", "MGUIExpoTrappingCorrection", this, "OnApply()"); + + // Apply Button + m_ButtonApply = new TGTextButton(ControlFrame, " Apply Range "); + ControlFrame->AddFrame(m_ButtonApply, new TGLayoutHints(kLHintsLeft | kLHintsCenterY, 5, 5, 2, 2)); + m_ButtonApply->Connect("Clicked()", "MGUIExpoTrappingCorrection", this, "OnApply()"); + + AddFrame(ControlFrame, new TGLayoutHints(kLHintsTop | kLHintsExpandX, 5, 5, 5, 2)); + + // Main Canvas Frame (Side-by-Side Plots) + TGLayoutHints* CanvasLayout = new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX | kLHintsExpandY, 2, 2, 2, 2); + TGHorizontalFrame* HFrame = new TGHorizontalFrame(this); + AddFrame(HFrame, CanvasLayout); + + // LV Canvas setup + m_CanvasLV = new TRootEmbeddedCanvas("CanvasLV", HFrame, 100, 100); + HFrame->AddFrame(m_CanvasLV, CanvasLayout); + + m_CanvasLV->GetCanvas()->cd(); + m_CanvasLV->GetCanvas()->SetGridx(); + m_CanvasLV->GetCanvas()->SetGridy(); + + m_EnergyLVInitial->Draw("HIST"); + m_EnergyLVFinal->Draw("HIST SAME"); + + m_LegendLV = new TLegend(0.55, 0.72, 0.88, 0.88); + m_LegendLV->AddEntry(m_EnergyLVInitial, "LV Uncorrected", "l"); + m_LegendLV->AddEntry(m_EnergyLVFinal, "LV Corrected", "f"); + m_LegendLV->Draw(); + + // HV Canvas setup + m_CanvasHV = new TRootEmbeddedCanvas("CanvasHV", HFrame, 100, 100); + HFrame->AddFrame(m_CanvasHV, CanvasLayout); + + m_CanvasHV->GetCanvas()->cd(); + m_CanvasHV->GetCanvas()->SetGridx(); + m_CanvasHV->GetCanvas()->SetGridy(); + + m_EnergyHVInitial->Draw("HIST"); + m_EnergyHVFinal->Draw("HIST SAME"); + + m_LegendHV = new TLegend(0.55, 0.72, 0.88, 0.88); + m_LegendHV->AddEntry(m_EnergyHVInitial, "HV Uncorrected", "l"); + m_LegendHV->AddEntry(m_EnergyHVFinal, "HV Corrected", "f"); + m_LegendHV->Draw(); + + MapSubwindows(); + Layout(); + + // Signal that the canvas has been created so we know when we can start to fill it + m_IsCreated = true; + + m_Mutex.UnLock(); +} + +//////////////////////////////////////////////////////////////////////////////// + +void MGUIExpoTrappingCorrection::OnApply() +{ + //Create function has already been run + if (m_IsCreated == false) return; + if (m_EntryNBins == nullptr || m_EntryMinEnergy == nullptr || m_EntryMaxEnergy == nullptr) return; + + int nBins = m_EntryNBins->GetIntNumber(); + double minE = m_EntryMinEnergy->GetNumber(); + double maxE = m_EntryMaxEnergy->GetNumber(); + + if (maxE <= minE || nBins <= 0) return; + + SetEnergyHistogramParameters(nBins, minE, maxE); + Update(); +} + +//////////////////////////////////////////////////////////////////////////////// + +void MGUIExpoTrappingCorrection::SetEnergyHistogramParameters(int NBins, double Min, double Max) +{ + m_Mutex.Lock(); + + // 1. Zoom the X-axis view without altering underlying data/binning + if (m_EnergyLVInitial != nullptr) m_EnergyLVInitial->GetXaxis()->SetRangeUser(Min, Max); + if (m_EnergyLVFinal != nullptr) m_EnergyLVFinal->GetXaxis()->SetRangeUser(Min, Max); + if (m_EnergyHVInitial != nullptr) m_EnergyHVInitial->GetXaxis()->SetRangeUser(Min, Max); + if (m_EnergyHVFinal != nullptr) m_EnergyHVFinal->GetXaxis()->SetRangeUser(Min, Max); + + // 2. Update entry boxes + if (m_EntryNBins != nullptr) m_EntryNBins->SetIntNumber(NBins); + if (m_EntryMinEnergy != nullptr) m_EntryMinEnergy->SetNumber(Min); + if (m_EntryMaxEnergy != nullptr) m_EntryMaxEnergy->SetNumber(Max); + + m_Mutex.UnLock(); +} +//////////////////////////////////////////////////////////////////////////////// + +void MGUIExpoTrappingCorrection::Update() +{ + if (m_IsCreated == false) return; + + m_Mutex.Lock(); + + bool isLog = (m_CheckLogY != nullptr && m_CheckLogY->IsOn()); + + // --- Update LV Canvas --- + if (m_CanvasLV != nullptr && m_CanvasLV->GetCanvas() != nullptr) { + TCanvas* canvasLV = m_CanvasLV->GetCanvas(); + canvasLV->cd(); + canvasLV->SetLogy(isLog ? 1 : 0); + + double maxLVInitial = m_EnergyLVInitial->GetBinContent(m_EnergyLVInitial->GetMaximumBin()); + double maxLVFinal = m_EnergyLVFinal->GetBinContent(m_EnergyLVFinal->GetMaximumBin()); + double realMaxLV = std::max(maxLVInitial, maxLVFinal); + + if (isLog) { + m_EnergyLVInitial->SetMinimum(0.1); + m_EnergyLVInitial->SetMaximum(realMaxLV > 0 ? realMaxLV * 5.0 : 10.0); + } else { + m_EnergyLVInitial->SetMinimum(-1111); + m_EnergyLVInitial->SetMaximum(realMaxLV > 0 ? realMaxLV * 1.15 : 10.0); + } + + canvasLV->Modified(); + canvasLV->Update(); + } + + // --- Update HV Canvas --- + if (m_CanvasHV != nullptr && m_CanvasHV->GetCanvas() != nullptr) { + TCanvas* canvasHV = m_CanvasHV->GetCanvas(); + canvasHV->cd(); + canvasHV->SetLogy(isLog ? 1 : 0); + + double maxHVInitial = m_EnergyHVInitial->GetBinContent(m_EnergyHVInitial->GetMaximumBin()); + double maxHVFinal = m_EnergyHVFinal->GetBinContent(m_EnergyHVFinal->GetMaximumBin()); + double realMaxHV = std::max(maxHVInitial, maxHVFinal); + + if (isLog) { + m_EnergyHVInitial->SetMinimum(0.1); + m_EnergyHVInitial->SetMaximum(realMaxHV > 0 ? realMaxHV * 5.0 : 10.0); + } else { + m_EnergyHVInitial->SetMinimum(-1111); + m_EnergyHVInitial->SetMaximum(realMaxHV > 0 ? realMaxHV * 1.15 : 10.0); + } + + canvasHV->Modified(); + canvasHV->Update(); + } + + m_Mutex.UnLock(); +} + +//////////////////////////////////////////////////////////////////////////////// + +void MGUIExpoTrappingCorrection::Reset() +{ + m_Mutex.Lock(); + + if (m_EnergyLVInitial != nullptr) m_EnergyLVInitial->Reset(); + if (m_EnergyLVFinal != nullptr) m_EnergyLVFinal->Reset(); + if (m_EnergyHVInitial != nullptr) m_EnergyHVInitial->Reset(); + if (m_EnergyHVFinal != nullptr) m_EnergyHVFinal->Reset(); + + m_Mutex.UnLock(); +} + +//////////////////////////////////////////////////////////////////////////////// + +void MGUIExpoTrappingCorrection::AddEnergyInitial(double Energy, bool IsNearestNeighbor, bool IsLV) +{ + m_Mutex.Lock(); + + if (IsLV == true) { + if (m_EnergyLVInitial != nullptr) m_EnergyLVInitial->Fill(Energy); + } else { + if (m_EnergyHVInitial != nullptr) m_EnergyHVInitial->Fill(Energy); + } + + m_Mutex.UnLock(); +} + +//////////////////////////////////////////////////////////////////////////////// + +void MGUIExpoTrappingCorrection::AddEnergyFinal(double Energy, bool IsNearestNeighbor, bool IsLV) +{ + m_Mutex.Lock(); + + if (IsLV == true) { + if (m_EnergyLVFinal != nullptr) m_EnergyLVFinal->Fill(Energy); + } else { + if (m_EnergyHVFinal != nullptr) m_EnergyHVFinal->Fill(Energy); + } + + m_Mutex.UnLock(); +} + +//////////////////////////////////////////////////////////////////////////////// + +void MGUIExpoTrappingCorrection::Export(const MString& FileName) +{ + m_Mutex.Lock(); + + if (m_CanvasLV != nullptr && m_CanvasLV->GetCanvas() != nullptr) { + MString nameLV = FileName; + nameLV.ReplaceAll(".png", "_LV.png"); + m_CanvasLV->GetCanvas()->SaveAs(nameLV); + } + + if (m_CanvasHV != nullptr && m_CanvasHV->GetCanvas() != nullptr) { + MString nameHV = FileName; + nameHV.ReplaceAll(".png", "_HV.png"); + m_CanvasHV->GetCanvas()->SaveAs(nameHV); + } + + m_Mutex.UnLock(); +} \ No newline at end of file diff --git a/src/MGUIOptionsLoaderMeasurementsHDF.cxx b/src/MGUIOptionsLoaderMeasurementsHDF.cxx index a763f1bc..eedb6add 100644 --- a/src/MGUIOptionsLoaderMeasurementsHDF.cxx +++ b/src/MGUIOptionsLoaderMeasurementsHDF.cxx @@ -85,14 +85,13 @@ void MGUIOptionsLoaderMeasurementsHDF::Create() dynamic_cast(m_Module)->GetFileNameStripMap()); m_FileSelectorStripMap->SetFileType("Strip map file", "*.map"); m_OptionsFrame->AddFrame(m_FileSelectorStripMap, LabelLayout); - + // Nearest neighbor checkbox m_IncludeNearestNeighbor = new TGCheckButton(m_OptionsFrame, "Include Nearest Neighbors"); m_IncludeNearestNeighbor->SetOn(dynamic_cast(m_Module)->GetIncludeNearestNeighbor()); m_IncludeNearestNeighbor->Associate(this); m_OptionsFrame->AddFrame(m_IncludeNearestNeighbor, LabelLayout); - PostCreate(); } diff --git a/src/MGUIOptionsTrappingCorrection.cxx b/src/MGUIOptionsTrappingCorrection.cxx new file mode 100644 index 00000000..117ba531 --- /dev/null +++ b/src/MGUIOptionsTrappingCorrection.cxx @@ -0,0 +1,99 @@ +/* + * MGUIOptionsTrappingCorrection.cxx + * + * + * Copyright (C) by Andreas Zoglauer. + * All rights reserved. + * + * + * This code implementation is the intellectual property of + * Andreas Zoglauer. + * + * By copying, distributing or modifying the Program (or any work + * based on the Program) you indicate your acceptance of this statement, + * and all its terms. + * + */ + + +// Include the header: +#include "MGUIOptionsTrappingCorrection.h" + +// Standard libs: + +// ROOT libs: +#include +#include +#include +#include + +// MEGAlib libs: +#include "MStreams.h" +#include "MModule.h" +#include "MModuleTrappingCorrection.h" + + +//////////////////////////////////////////////////////////////////////////////// + + +#ifdef ___CLING___ +ClassImp(MGUIOptionsTrappingCorrection) +#endif + + +//////////////////////////////////////////////////////////////////////////////// + + +MGUIOptionsTrappingCorrection::MGUIOptionsTrappingCorrection(MModule* Module) + : MGUIOptions(Module) +{ + // standard constructor +} + + +//////////////////////////////////////////////////////////////////////////////// + + +MGUIOptionsTrappingCorrection::~MGUIOptionsTrappingCorrection() +{ + // kDeepCleanup is activated +} + + +//////////////////////////////////////////////////////////////////////////////// + + +void MGUIOptionsTrappingCorrection::Create() +{ + PreCreate(); + + m_SimCCEFileSelector = new MGUIEFileSelector(m_OptionsFrame, "Select a trapping parameter file:", + dynamic_cast(m_Module)->GetSimCCEFileName()); + m_SimCCEFileSelector->SetFileType("trapping parameters", "*.csv"); + TGLayoutHints* LabelLayout = new TGLayoutHints(kLHintsTop | kLHintsCenterX | kLHintsExpandX, 10, 10, 10, 10); + m_OptionsFrame->AddFrame(m_SimCCEFileSelector, LabelLayout); + + TGLayoutHints* RBLayout = new TGLayoutHints(kLHintsLeft | kLHintsTop, 40, 10, 2, 0); + TGLayoutHints* RBOptionLayout = new TGLayoutHints(kLHintsLeft | kLHintsTop, 60, 10, 2, 0); + TGLayoutHints* RBOptionStretchLayout = new TGLayoutHints(kLHintsLeft | kLHintsTop | kLHintsExpandX, 60, 10, 2, 0); + + + PostCreate(); +} + + +//////////////////////////////////////////////////////////////////////////////// + + +bool MGUIOptionsTrappingCorrection::OnApply() +{ + // Modify this to store the data in the module! + + dynamic_cast(m_Module)->SetSimCCEFileName(m_SimCCEFileSelector->GetFileName()); + + return true; +} + + +// MGUIOptionsTrappingCorrection: the end... +//////////////////////////////////////////////////////////////////////////////// diff --git a/src/MHit.cxx b/src/MHit.cxx index 260352c4..8e4e4cca 100644 --- a/src/MHit.cxx +++ b/src/MHit.cxx @@ -74,6 +74,8 @@ void MHit::Clear() m_Position = g_VectorNotDefined; m_Energy = g_DoubleNotDefined; + m_LocalPosition = g_VectorNotDefined; + m_LVEnergy = g_DoubleNotDefined; m_HVEnergy = g_DoubleNotDefined; diff --git a/src/MModuleDepthCalibration.cxx b/src/MModuleDepthCalibration.cxx index adcb9044..3b28d6e2 100644 --- a/src/MModuleDepthCalibration.cxx +++ b/src/MModuleDepthCalibration.cxx @@ -365,11 +365,13 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) // Make sure XYZ resolution are correctly mapped to the global coord system. MVector PositionResolution(Xsigma, Ysigma, Zsigma); MVector GlobalResolution = ((m_Detectors[DetID]->GetSensitiveVolume(0)->GetPositionInWorldVolume(PositionResolution)) - (m_Detectors[DetID]->GetSensitiveVolume(0)->GetPositionInWorldVolume(LocalOrigin))).Abs(); - + H->SetPosition(GlobalPosition); H->SetPositionResolution(GlobalResolution); + H->SetLocalPosition(LocalPosition); + } @@ -440,23 +442,45 @@ bool MModuleDepthCalibration::LoadDetectorDimensions(MDGeometryQuest* Geometry) // ie DetID=0 should be the 0th detector in m_Detectors, DetID=1 should the 1st, etc. vector DetList = Geometry->GetDetectorList(); - // Look through the Geometry and get the names and thicknesses of all the detectors. + // Look through the Geometry and get the names and thicknesses of all Strip3D detectors. + vector DetectorNames; + unsigned int DetID = 0; + for (unsigned int i = 0; i < DetList.size(); ++i) { // For now, DetID is in order of detectors, which puts contraints on how the geometry file should be written. // If using the card cage at UCSD, default to DetID=11. - unsigned int DetID = i; if (m_UCSDOverride == true) { DetID = 11; } MDDetector* det = DetList[i]; - vector DetectorNames; if (det->GetTypeName() == "Strip3D") { if (det->GetNSensitiveVolumes() == 1) { MDVolume* vol = det->GetSensitiveVolume(0); - string det_name = vol->GetName().GetString(); - if (find(DetectorNames.begin(), DetectorNames.end(), det_name) == DetectorNames.end()) { - DetectorNames.push_back(det_name); + MString DetectorName = det->GetName(); + string DetName = DetectorName.GetString(); + + // Check that the DetID agrees with the naming scheme GeD_X + if (DetectorName.BeginsWith("GeD_") == true) { + DetectorName.RemoveAllInPlace("GeD_"); // The number after GeD is the COSI detector ID + if (DetID != DetectorName.ToUnsignedInt()) { + if (g_Verbosity >= c_Error) { + cout << "ERROR in MModuleDepthCalibration::Initialize: Non-matching DetID="<GetName() == "COSI-SMEX-Payload"){ + return false; + } + } + } else if (Geometry->GetName() == "COSI-SMEX-Payload") { + if (g_Verbosity >= c_Error) { + cout << "ERROR in MModuleDepthCalibration::Initialize: COSI-SMEX-Payload expects all Strip3D detectors to follow the name scheme GeD_X"<GetStructuralSize().GetZ()); MDStrip3D* strip = dynamic_cast(det); m_XPitches[DetID] = strip->GetPitchX(); @@ -465,7 +489,7 @@ bool MModuleDepthCalibration::LoadDetectorDimensions(MDGeometryQuest* Geometry) m_NYStrips[DetID] = strip->GetNStripsY(); if (g_Verbosity >= c_Info) { - cout << "Found detector " << det_name << " corresponding to DetID=" << DetID << "." << endl; + cout << "Found detector " << DetName << " corresponding to DetID=" << DetID << "." << endl; cout << "Detector thickness: " << m_Thicknesses[DetID] << endl; cout << "Number of X strips: " << m_NXStrips[DetID] << endl; cout << "Number of Y strips: " << m_NYStrips[DetID] << endl; @@ -474,9 +498,10 @@ bool MModuleDepthCalibration::LoadDetectorDimensions(MDGeometryQuest* Geometry) } m_DetectorIDs.push_back(DetID); m_Detectors[DetID] = det; + DetID += 1; } else { if (g_Verbosity >= c_Error) { - cout<<"ERROR in MModuleDepthCalibration::Initialize: Found a duplicate detector: "<GetAvailableModuleByXmlTag("EnergyCalibration"); + if (m_EnergyCalibration == nullptr) { + if (g_Verbosity >= c_Error) { + cout << "ERROR in MModuleTrappingCorrection::Initialize: couldn't resolve pointer to Energy Calibration Module... need access to this module for energy resolution lookup!" << endl; + } + return false; + } + + m_DepthCalibration = (MModuleDepthCalibration*) S->GetAvailableModuleByXmlTag("DepthCalibration"); + if (m_DepthCalibration == nullptr) { + if (g_Verbosity >= c_Error) { + cout << "ERROR in MModuleTrappingCorrection::Initialize: couldn't resolve pointer to Depth Calibration Module... need access to this module for depth resolution lookup!" << endl; + } + return false; + } + +// TO DO: remove this check once we successfully process multiple detectors + m_DetectorIDs = m_DepthCalibration-> GetDetectorIDs(); + + if (m_DetectorIDs.empty()) { + if (g_Verbosity >= c_Error) { + cout << "ERROR in MModuleTrappingCorrection::Initialize: Depth Calibration has no registered detector IDs!" << endl; + } + return false; + } + + return MModule::Initialize(); +} + + +//////////////////////////////////////////////////////////////////////////////// + +void MModuleTrappingCorrection::CreateExpos() +{ + // Create all expos + + if (HasExpos() == true) return; + + // Set the histogram display using the new double-canvas GUI class + m_ExpoSpectrum = new MGUIExpoTrappingCorrection(this); + m_ExpoSpectrum->SetEnergyHistogramParameters(200, 0, 2000); + m_Expos.push_back(m_ExpoSpectrum); +} + + +///////////////////////////////////////////////////////////////////////////////// + + +bool MModuleTrappingCorrection::AnalyzeEvent(MReadOutAssembly* Event) +{ + if (Event->GetGuardRingVeto() == true) { + // Right now we cannot use events w GR veto + return false; + } else { + for (unsigned int i = 0; i < Event->GetNHits(); ++i) { + MHit* H = Event->GetHit(i); + + int Grade = m_DepthCalibration->GetHitGrade(H); + + if (Grade < 0 || Grade > 4) { + H->SetNoDepth(); + } else { // If Grade is 0-4, proceed with analysis + + vector LVStrips; + vector HVStrips; + + for (unsigned int j = 0; j < H->GetNStripHits(); ++j) { + MStripHit* SH = H->GetStripHit(j); + if (SH->IsLowVoltageStrip()) LVStrips.push_back(SH); else HVStrips.push_back(SH); + } + + double LVEnergyFraction; + double HVEnergyFraction; + MStripHit* LVSH = m_DepthCalibration->GetDominantStrip(LVStrips, LVEnergyFraction); + MStripHit* HVSH = m_DepthCalibration->GetDominantStrip(HVStrips, HVEnergyFraction); + + // Local Z depth position + double depth_val = H->GetLocalPosition().GetZ(); + + // --- Low Voltage Side --- + if (LVSH != nullptr) { + double rawLVEnergy = LVSH->GetEnergy(); + + // --- 1. FILL UNCORRECTED (RAW) --- + if (HasExpos() == true) { + m_ExpoSpectrum->AddEnergyInitial(rawLVEnergy, LVSH->IsNearestNeighbor(), LVSH->IsLowVoltageStrip()); + } + + // --- 2. CALCULATE CORRECTION --- + double correctedLVEnergy = GetSimBasedCorrectedEnergy(depth_val, rawLVEnergy, m_CCEs_LV_e, m_CCEs_LV_h, m_ParamA_LV, m_ParamB, m_ParamC); + LVSH->SetEnergy(correctedLVEnergy); + + // --- 3. FILL CORRECTED (FINAL) --- + if (HasExpos() == true) { + m_ExpoSpectrum->AddEnergyFinal(correctedLVEnergy, LVSH->IsNearestNeighbor(), LVSH->IsLowVoltageStrip()); + } + } + + // --- High Voltage Side --- + if (HVSH != nullptr) { + double rawHVEnergy = HVSH->GetEnergy(); + + // 1. Record UNCORRECTED (raw) HV energy to expo spectrum + if (HasExpos() == true) { + m_ExpoSpectrum->AddEnergyInitial(rawHVEnergy, HVSH->IsNearestNeighbor(), HVSH->IsLowVoltageStrip()); + } + + // 2. Compute trapping correction + double correctedHVEnergy = GetSimBasedCorrectedEnergy(depth_val, rawHVEnergy, m_CCEs_HV_e, m_CCEs_HV_h, m_ParamA_HV, m_ParamB, m_ParamC); + HVSH->SetEnergy(correctedHVEnergy); + + // 3. Record CORRECTED (final) HV energy to expo spectrum + if (HasExpos() == true) { + m_ExpoSpectrum->AddEnergyFinal(correctedHVEnergy, HVSH->IsNearestNeighbor(), HVSH->IsLowVoltageStrip()); + } + } + } + } + } + + Event->SetAnalysisProgress(MAssembly::c_TrappingCorrection); + return true; +} + +///////////////////////////////////////////////////////////////////////////////// + +void MModuleTrappingCorrection::Finalize() +{ + MModule::Finalize(); + + if (m_ExpoSpectrum == nullptr) { + cout << "ERROR in MModuleTrappingCorrection::Finalize: Expo plot spectrum is null." << endl; + return; + } + + if (g_Verbosity >= c_Info) { + cout << "INFO: Finalizing Trapping Correction Module..." << endl; + } + + // Retrieve both uncorrected and corrected histograms + TH1D* histLVInit = m_ExpoSpectrum->GetEnergyHistogramLVInitial(); + TH1D* histLVFinal = m_ExpoSpectrum->GetEnergyHistogramLVFinal(); + + TH1D* histHVInit = m_ExpoSpectrum->GetEnergyHistogramHVInitial(); + TH1D* histHVFinal = m_ExpoSpectrum->GetEnergyHistogramHVFinal(); + + // Helper lambda to perform fit, output results, and return {mu, fwhm} + auto FitAndPrintSpectrum = [&](TH1D* hist, const string& titleLabel) -> std::pair { + if (hist == nullptr || hist->GetEntries() <= 0) { + cout << "WARNING: " << titleLabel << " histogram is null or has 0 entries." << endl; + return std::make_pair(0.0, 0.0); + } + + double directFWHM = CalculateDirectFWHM(hist); + + TF1* fitFunc = GeneratePhotopeakFunction(); + fitFunc->SetParameter("Amplitude", hist->GetBinContent(hist->GetMaximumBin())); + hist->Fit(fitFunc, "RQ"); + + double mu = fitFunc->GetParameter("x0 (Mu)"); + double fwhm = 2.35482 * fitFunc->GetParameter("Sigma Gauss"); + + cout << "\n" << m_XmlTag << " --- " << titleLabel << " ---" << endl; + cout << " Centroid (Mu) : " << mu << " keV" << endl; + cout << " Fitted Gaussian FWHM : " << fwhm << " keV" << endl; + cout << " Direct Histogram FWHM: " << directFWHM << " keV" << endl; + + delete fitFunc; + return std::make_pair(mu, fwhm); + }; + + cout << "\n========================================================" << endl; + cout << " TRAPPING CORRECTION SPECTRUM FIT RESULTS " << endl; + cout << "========================================================" << endl; + + // --- Fit Uncorrected and Corrected Spectra --- + std::pair lv_raw = FitAndPrintSpectrum(histLVInit, "LV UNCORRECTED (RAW) SPECTRUM"); + std::pair lv_corr = FitAndPrintSpectrum(histLVFinal, "LV CORRECTED SPECTRUM"); + + std::pair hv_raw = FitAndPrintSpectrum(histHVInit, "HV UNCORRECTED (RAW) SPECTRUM"); + std::pair hv_corr = FitAndPrintSpectrum(histHVFinal, "HV CORRECTED SPECTRUM"); + + // Summary Comparison Output + cout << "\n========================================================" << endl; + cout << " SUMMARY COMPARISON " << endl; + cout << "========================================================" << endl; + cout << " LV Side: " << endl; + cout << " Raw Centroid: " << lv_raw.first << " keV | Corrected Centroid: " << lv_corr.first << " keV" << endl; + cout << " Raw FWHM : " << lv_raw.second << " keV | Corrected FWHM : " << lv_corr.second << " keV" << endl; + cout << " HV Side: " << endl; + cout << " Raw Centroid: " << hv_raw.first << " keV | Corrected Centroid: " << hv_corr.first << " keV" << endl; + cout << " Raw FWHM : " << hv_raw.second << " keV | Corrected FWHM : " << hv_corr.second << " keV" << endl; + cout << "========================================================\n" << endl; + + return; +} + +///////////////////////////////////////////////////////////////////////////////// + + +bool MModuleTrappingCorrection::LoadSimCCEFile(MString FileName) +{ + MFile SimCCEFile; + if (SimCCEFile.Open(FileName) == false) { + if (g_Verbosity >= c_Error) { + cout << "ERROR in MModuleTrappingCorrection::LoadSimCCEFile: failed to open file " << FileName << endl; + } + return false; + } + + m_Depths.clear(); + m_CCEs_HV_e.clear(); + m_CCEs_HV_h.clear(); + m_CCEs_LV_e.clear(); + m_CCEs_LV_h.clear(); + + MString Line; + bool ParsedHeaderParameters = false; + + while (SimCCEFile.ReadLine(Line)) { + + // Skip empty lines or pure comment lines + if (Line.IsEmpty() == true || Line.BeginsWith('#') == true) { + continue; + } + + std::vector Tokens = Line.Tokenize(","); + + // Read parameters (A_HV, A_LV, B, C) + if (ParsedHeaderParameters == false) { + if (Tokens.size() == 4) { + m_ParamA_HV = Tokens[0].Strip().ToDouble(); + m_ParamA_LV = Tokens[1].Strip().ToDouble(); + m_ParamB = Tokens[2].Strip().ToDouble(); + m_ParamC = Tokens[3].Strip().ToDouble(); + + ParsedHeaderParameters = true; + } else { + if (g_Verbosity >= c_Error) { + cout << "ERROR in LoadSimCCEFile: Expected 4 parameters (A_HV, A_LV, B, C) on first data line, found " + << Tokens.size() << " tokens." << endl; + } + SimCCEFile.Close(); + return false; + } + } + // Read CCE curves + else { + if (Tokens.size() == 5) { + m_Depths.push_back(Tokens[0].Strip().ToDouble()); + m_CCEs_HV_e.push_back(Tokens[1].Strip().ToDouble()); + m_CCEs_HV_h.push_back(Tokens[2].Strip().ToDouble()); + m_CCEs_LV_e.push_back(Tokens[3].Strip().ToDouble()); + m_CCEs_LV_h.push_back(Tokens[4].Strip().ToDouble()); + } + } + } + + SimCCEFile.Close(); + + if (ParsedHeaderParameters == false || m_Depths.size() == 0) { + if (g_Verbosity >= c_Error) { + cout << "ERROR in LoadSimCCEFile: No valid CCE data points were loaded!" << endl; + } + return false; + } + + // Console output on successful load + if (g_Verbosity >= c_Info) { + cout << m_XmlTag << "Loaded CCE simulation parameters from " << FileName << ":" << endl; + cout << m_XmlTag << " A_HV = " << m_ParamA_HV << ", A_LV = " << m_ParamA_LV + << ", B = " << m_ParamB << ", C = " << m_ParamC << endl; + cout << m_XmlTag << " Loaded " << m_Depths.size() << " depth grid points." << endl; + } + + return true; +} + +///////////////////////////////////////////////////////////////////////////////// + +double MModuleTrappingCorrection::GetSimBasedCorrectedEnergy(double depth_val, double uncorrected_energy, const std::vector& sim_cce_sorted_e, const std::vector& sim_cce_sorted_h, double paramA, double paramB, double paramC) { + + // Look up the simulation CCE baseline using the interpolate function + + double cce_base_e = Interpolate(depth_val, m_Depths, sim_cce_sorted_e); + double cce_base_h = Interpolate(depth_val, m_Depths, sim_cce_sorted_h); + + + // Evaluate the physical trapping function model using class global popt variables + double expected_centroid_scaled = paramA * (1.0 - paramB * (1.0 - cce_base_e)) * (1.0 - paramC * (1.0 - cce_base_h)); + + // Prevent division-by-zero or non-physical negative values + if (expected_centroid_scaled <= 0.0) { + return uncorrected_energy; + } + + //Reconstruct the true un-trapped energy + return uncorrected_energy / expected_centroid_scaled; +} + + +///////////////////////////////////////////////////////////////////////////////// + + +double MModuleTrappingCorrection::Interpolate(double x, const std::vector& xp, const std::vector& fp) { + // need an interpolation function to get continuous CCE values from the discrete simulation data + if (xp.empty()) return 0.0; + if (x <= xp.front()) return fp.front(); + if (x >= xp.back()) return fp.back(); + + // Find the first element which is greater than or equal to x + auto it = std::lower_bound(xp.begin(), xp.end(), x); + size_t idx = std::distance(xp.begin(), it); + + // Linear interpolation formula + double x0 = xp[idx - 1]; + double x1 = xp[idx]; + double y0 = fp[idx - 1]; + double y1 = fp[idx]; + + return y0 + (x - x0) * (y1 - y0) / (x1 - x0); +} + +///////////////////////////////////////////////////////////////////////////////// + + +void MModuleTrappingCorrection::ShowOptionsGUI() +{ + // Show the options GUI - or do nothing + MGUIOptionsTrappingCorrection* Options = new MGUIOptionsTrappingCorrection(this); + Options->Create(); + gClient->WaitForUnmap(Options); +} + + +///////////////////////////////////////////////////////////////////////////////// + + +bool MModuleTrappingCorrection::ReadXmlConfiguration(MXmlNode* Node) +{ + //! Read the configuration data from an XML node + + MXmlNode* SimCCEFileNameNode = Node->GetNode("SimCCEFileName"); + if (SimCCEFileNameNode != nullptr) { + m_SimCCEFileName = SimCCEFileNameNode->GetValue(); + } + + return true; +} + + +///////////////////////////////////////////////////////////////////////////////// + +MXmlNode* MModuleTrappingCorrection::CreateXmlConfiguration() +{ + //! Create an XML node tree from the configuration + + MXmlNode* Node = new MXmlNode(0,m_XmlTag); + new MXmlNode(Node, "SimCCEFileName", m_SimCCEFileName); + + return Node; +} + +////////////////////////////////////////////////////////////////////////////// + + +TF1* MModuleTrappingCorrection::GeneratePhotopeakFunction() +{ + // Component 1: Core Gaussian + // exp(-(x-x0)^2 / (2*sigma^2)) + MString gaussStr = "exp(-(x-[1])^2 / (2*[2]^2))"; + + // Component 2: Exponential Tail + Shelf + // BoverA * exp(gamma*(x-x0)) * 0.5 * erfc((x-x0)/(sigma*sigma_ratio*sqrt(2))) + MString expTailStr = "[3] * exp([4]*(x-[1])) * 0.5 * erfc((x-[1])/([2]*[5]*sqrt(2)))"; + + // Component 3: Linear Tail + Shelf + // BoverA * CoverB * (1 + D*(x-x0)) * 0.5 * erfc((x-x0)/(sigma*sigma_ratio*sqrt(2))) + MString linTailStr = "[3] * [6] * (1 + [7]*(x-[1])) * 0.5 * erfc((x-[1])/([2]*[5]*sqrt(2)))"; + + // Combine components with an overall normalization scaling factor [0] + MString fullFormula = "[0] * (" + gaussStr + " + " + expTailStr + " + " + linTailStr + ")"; + + // Instantiate TF1 over your expected fit window + TF1* PhotopeakFunction = new TF1("PhotopeakFunction", fullFormula.Data(), 645, 675); + + // Set Parameter Names + PhotopeakFunction->SetParName(0, "Amplitude"); + PhotopeakFunction->SetParName(1, "x0 (Mu)"); + PhotopeakFunction->SetParName(2, "Sigma Gauss"); + PhotopeakFunction->SetParName(3, "BoverA"); + PhotopeakFunction->SetParName(4, "Gamma"); + PhotopeakFunction->SetParName(5, "Sigma Ratio"); + PhotopeakFunction->SetParName(6, "CoverB"); + PhotopeakFunction->SetParName(7, "D (Lin Slope)"); + + // Provide initial sensible guesses for a Cs137 photopeak + PhotopeakFunction->SetParameter("Amplitude", 1000); + PhotopeakFunction->SetParameter("x0 (Mu)", 661.7); + PhotopeakFunction->SetParameter("Sigma Gauss", 2.0); + PhotopeakFunction->SetParameter("BoverA", 0.05); + PhotopeakFunction->SetParameter("Gamma", 0.5); + PhotopeakFunction->SetParameter("Sigma Ratio", 0.85); + PhotopeakFunction->SetParameter("CoverB", 0.13); + PhotopeakFunction->SetParameter("D (Lin Slope)", 0.028); + + // Set boundary limits to stabilize convergence + PhotopeakFunction->SetParLimits(0, 1, 1e8); + PhotopeakFunction->SetParLimits(1, 645, 675); // Keeps peak centered around 662 keV + PhotopeakFunction->SetParLimits(2, 0.5, 10); // Prevents sigma from blowing up or hitting zero + PhotopeakFunction->SetParLimits(3, 0.0, 1.0); // Tail shouldn't be larger than the main peak + PhotopeakFunction->SetParLimits(4, 0.001, 2.0); // Standard range for exponential decay factor + PhotopeakFunction->SetParLimits(5, 0.1, 5.0); // Ratio of shelf width to peak width + PhotopeakFunction->SetParLimits(6, 0.0, 5.0); + PhotopeakFunction->SetParLimits(7, -1.0, 1.0); + + return PhotopeakFunction; +} + + +//////////////////////////////////////////////////////////////////////////////// +double MModuleTrappingCorrection::CalculateDirectFWHM(TH1D* hist) +{ + if (hist == nullptr || hist->GetEntries() == 0) return 0.0; + + //Locate the peak bin + int maxBin = hist->GetMaximumBin(); + + if (g_Verbosity >= c_Info){ + cout << "INFO: Maximum bin located at: " << maxBin << " with content: " << hist->GetBinContent(maxBin) << endl; + } + + double peakHeight = hist->GetBinContent(maxBin); + if (peakHeight <= 0.0) return 0.0; + + //Define a local search window (depends on binning you set on the GUI) + // Set the binning in the GUI to 0.1 keV before finalize function is called + const int searchWindowBins = 150; + + int minSearchBin = std::max(1, maxBin - searchWindowBins); + int maxSearchBin = std::min(hist->GetNbinsX(), maxBin + searchWindowBins); + + //Estimate local background level from the window edges + double bgLeft = hist->GetBinContent(minSearchBin); + double bgRight = hist->GetBinContent(maxSearchBin); + double localBG = 0.0; // assume the background is zero for now + //(bgLeft + bgRight) / 2.0; + + // Calculate net peak height above background + double netPeakHeight = peakHeight - localBG; + if (netPeakHeight <= 0.0) return 0.0; + + // Target level is half-maximum relative to local background + double targetHalfMax = localBG + (netPeakHeight / 2.0); + + //Search left within the restricted window + double xLeft = -1.0; + for (int b = maxBin; b >= minSearchBin; --b) { + if (hist->GetBinContent(b) <= targetHalfMax) { + double x1 = hist->GetBinCenter(b); + double y1 = hist->GetBinContent(b); + double x2 = hist->GetBinCenter(b + 1); + double y2 = hist->GetBinContent(b + 1); + + // Linear interpolation between adjacent bins + xLeft = (y2 != y1) ? x1 + (targetHalfMax - y1) * (x2 - x1) / (y2 - y1) : x1; + break; + } + } + + //Search right within the restricted window + double xRight = -1.0; + for (int b = maxBin; b <= maxSearchBin; ++b) { + if (hist->GetBinContent(b) <= targetHalfMax) { + double x1 = hist->GetBinCenter(b - 1); + double y1 = hist->GetBinContent(b - 1); + double x2 = hist->GetBinCenter(b); + double y2 = hist->GetBinContent(b); + + // Linear interpolation between adjacent bins + xRight = (y2 != y1) ? x1 + (targetHalfMax - y1) * (x2 - x1) / (y2 - y1) : x2; + break; + } + } + + // Ensure valid crossing points were found on both sides + if (xLeft < 0.0 || xRight < 0.0) { + if (g_Verbosity >= c_Warning) { + cout << "WARNING in CalculateDirectFWHM: Peak did not cross half-maximum inside local window." << endl; + } + return 0.0; + } + + return (xRight - xLeft); +} +//////////////////////////////////////////////////////////////////////////////// + + +// MModuleTrappingCorrection.cxx: the end... +////////////////////////////////////////////////////////////////////////////////