diff --git a/apps/TrappingCorrectionAm241.cxx b/apps/TrappingCorrectionAm241.cxx index dadf4894..e650d358 100644 --- a/apps/TrappingCorrectionAm241.cxx +++ b/apps/TrappingCorrectionAm241.cxx @@ -274,13 +274,16 @@ bool TrappingCorrectionAm241::Analyze() // Store the HV and LV input files vector FileNames; FileNames.push_back(m_HVFileName); + cout<<"HV file names stored"< IllumSide; IllumSide.push_back(MString("HV")); IllumSide.push_back(MString("LV")); + cout<<"HV and LV integers mapped"< +#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 = 1500; + +int g_HVStrips = 64; +int g_LVStrips = 64; + +double g_CsPhotopeak = 661.7; + +const int NCTDBins = 50; +// 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* GenerateCTDFunction(double CTDFitMin, double CTDFitMax, double CTDGuess); + 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; + + double CTDFitMin = g_MinCTD; + double CTDFitMax = g_MaxCTD; + + 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)) { + + 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 (FullCTDHist == nullptr) { + // char name[128]; sprintf(name, "CTD_Detector%d_bin%d", DetID, CTDBin); + // FullCTDHist = new TH1D(name, name, (g_MaxCTD - g_MinCTD) / 2, g_MinCTD, g_MaxCTD); + // FullDetCTDHistograms[CTDBin][DetID] = FullCTDHist; + // } + 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) { + + // double CTDGuess = FullCTDHist->GetBinCenter(FullCTDHist->GetMaximumBin()); + // TF1* CTDFunction = GenerateCTDFunction(CTDFitMin, CTDFitMax, CTDGuess); + // TFitResultPtr CTDFit = FullCTDHist->Fit(CTDFunction, "SQ", "", CTDFitMin, CTDFitMax); + + 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 + + + // ofstream CTDFitFile(DetID + MString("_CTDFitResult_.txt")); + // streambuf* coutbuf = cout.rdbuf(); + // cout.rdbuf(CTDFitFile.rdbuf()); + // cout << "CTD Fit Results for Detector " << DetID << " in CTD bin " << c << endl; + // if (CTDFit >= 0) { + // CTDFit->Print(); + // } + // cout.rdbuf(coutbuf); + // CTDFitFile.close(); + + + 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 CTDFile(m_OutFile + MString("_Det") + DetID + MString("_CTDbin_") + c + MString("_CTDHist_Illum.root"), "recreate"); + // TCanvas* CTDCanvas = new TCanvas(); + // CTDCanvas->cd(); + // FullCTDHist->Draw("hist"); + // CTDFunction->Draw("same"); + // FullCTDHist->Write(); + // CTDFile.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; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +// TF1* TrappingCorrectionCs137::GenerateCTDFunction(double CTDFitMin, double CTDFitMax, double CTDGuess) +// { +// // Exponentially modified gaussian +// TF1* CTDFunction = new TF1("CTDFunction", "[0]*([1]/2)*exp(([1]/2)*(([1]*[3]*[3]) - 2*[4]*(x-[2])))*erfc((([1]*[3]*[3]) - [4]*(x-[2]))/([3]*sqrt(2)))", CTDFitMin, CTDFitMax); + +// CTDFunction->SetParName(0, "Norm"); +// CTDFunction->SetParName(1, "Lambda"); +// CTDFunction->SetParName(2, "Mu"); +// CTDFunction->SetParName(3, "Sigma"); +// CTDFunction->SetParName(4, "Flip"); + +// CTDFunction->SetParameter("Norm", 1000); +// CTDFunction->SetParameter("Lambda", 0.05); +// CTDFunction->SetParameter("Sigma", 12); +// CTDFunction->SetParameter("Mu", CTDGuess); + + +// CTDFunction->SetParLimits(0, 0, 1e8); +// CTDFunction->SetParLimits(1, 0.01, 1); +// CTDFunction->SetParLimits(2, CTDFitMin, CTDFitMax); +// CTDFunction->SetParLimits(3, 6, 30); + +// return CTDFunction; +// } + + +//////////////////////////////////////////////////////////////////////////////// + + +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!"< +#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 energy histogram + void AddEnergy(double Energy); + + // protected methods: + protected: + + + // protected members: + protected: + + // private members: + private: + //! Energy canvas + TRootEmbeddedCanvas* m_EnergyCanvas; + //! Energy histogram + TH1D* m_Energy; + + + +#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/MModuleTrappingCorrection.h b/include/MModuleTrappingCorrection.h new file mode 100644 index 00000000..32887ce2 --- /dev/null +++ b/include/MModuleTrappingCorrection.h @@ -0,0 +1,158 @@ +/* + * 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" + +// 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_SimCCEFile = FileName; } + + //! Get filename for SimCCE file + MString GetSimCCEFileName() const { return m_SimCCEFile; } + + //! Finalize the module + virtual void Finalize(); + + + //! Read the XML configuration + bool ReadXmlConfiguration(MXmlNode* Node); + + //! Create the XML configuration + MXmlNode* CreateXmlConfiguration(); + + + // 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); + + //! Determine the Grade (geometry of charge sharing) of the Hit + int GetHitGrade(MHit* H); + + //! 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, 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: + + unordered_map> m_SimCCE; + double m_SimCCE_Energy; + MString m_SimCCEFile; + + unordered_map m_Detectors; + vector m_DetectorIDs; + MModuleEnergyCalibration* m_EnergyCalibration; + MGUIExpoTrappingCorrection* m_ExpoTrappingCorrection; + + bool m_SimCCEFileIsLoaded; + + double m_ParamA_HV; + double m_ParamB_HV; + double m_ParamC_HV; + double m_ParamA_LV; + double m_ParamB_LV; + double m_ParamC_LV; + std::vector m_Depths; + std::vector m_CCEs_HV; + std::vector m_CCEs_LV; + + + + // private members: + private: + + //! Updated GUI to display the energy histogram + MGUIExpoPlotSpectrum* m_ExpoSpectrum; + + TF1* GeneratePhotopeakFunction(); + + +#ifdef ___CLING___ + public: + ClassDef(MModuleTrappingCorrection, 0) // no description +#endif + +}; + +#endif + + +//////////////////////////////////////////////////////////////////////////////// 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..e5e55b2d --- /dev/null +++ b/src/MGUIExpoTrappingCorrection.cxx @@ -0,0 +1,189 @@ +/* + * 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 + +// MEGAlib libs: +#include "MStreams.h" + + + +//////////////////////////////////////////////////////////////////////////////// + + +#ifdef ___CLING___ +ClassImp(MGUIExpoTrappingCorrection) +#endif + + +//////////////////////////////////////////////////////////////////////////////// + + +MGUIExpoTrappingCorrection::MGUIExpoTrappingCorrection(MModule* Module) : MGUIExpo(Module) +{ + // standard constructor + + // Set the new title of the tab here: + m_TabTitle = "Trapping Correction"; + + // Add all histograms and canvases below + m_Energy = new TH1D("", "Spectrum combined hits", 200, 0, 1000); + m_Energy->SetXTitle("Energy [keV]"); + m_Energy->SetYTitle("counts"); + m_Energy->SetFillColor(kAzure+7); + + m_EnergyCanvas = 0; + + // use hierarchical cleaning + SetCleanup(kDeepCleanup); +} + + +//////////////////////////////////////////////////////////////////////////////// + + +MGUIExpoTrappingCorrection::~MGUIExpoTrappingCorrection() +{ + // kDeepCleanup is activated +} + + +//////////////////////////////////////////////////////////////////////////////// + + +void MGUIExpoTrappingCorrection::Reset() +{ + //! Reset the data in the UI + + m_Mutex.Lock(); + + m_Energy->Reset(); + + m_Mutex.UnLock(); +} + + +//////////////////////////////////////////////////////////////////////////////// + + +void MGUIExpoTrappingCorrection::SetEnergyHistogramParameters(int NBins, double Min, double Max) +{ + // Set the energy histogram parameters + + m_Mutex.Lock(); + + m_Energy->SetBins(NBins, Min, Max); + + m_Mutex.UnLock(); +} + + +//////////////////////////////////////////////////////////////////////////////// + + +void MGUIExpoTrappingCorrection::AddEnergy(double Energy) +{ + // Add data to the energy histogram + + m_Mutex.Lock(); + + m_Energy->Fill(Energy); + + m_Mutex.UnLock(); +} + + +//////////////////////////////////////////////////////////////////////////////// + + +void MGUIExpoTrappingCorrection::Export(const MString& FileName) +{ + // Add data to the energy histogram + + m_Mutex.Lock(); + + m_EnergyCanvas->GetCanvas()->SaveAs(FileName); + + m_Mutex.UnLock(); +} + + +//////////////////////////////////////////////////////////////////////////////// + + +void MGUIExpoTrappingCorrection::Create() +{ + // Add the GUI options here + + // Do not create it twice! + if (m_IsCreated == true) return; + + m_Mutex.Lock(); + + TGLayoutHints* CanvasLayout = new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX | kLHintsExpandY, + 2, 2, 2, 2); + + TGHorizontalFrame* HFrame = new TGHorizontalFrame(this); + AddFrame(HFrame, CanvasLayout); + + + m_EnergyCanvas = new TRootEmbeddedCanvas("Energy", HFrame, 100, 100); + HFrame->AddFrame(m_EnergyCanvas, CanvasLayout); + + m_EnergyCanvas->GetCanvas()->cd(); + m_EnergyCanvas->GetCanvas()->SetGridy(); + m_EnergyCanvas->GetCanvas()->SetGridx(); + m_Energy->Draw(); + m_EnergyCanvas->GetCanvas()->Update(); + + m_IsCreated = true; + + m_Mutex.UnLock(); +} + + +//////////////////////////////////////////////////////////////////////////////// + + +void MGUIExpoTrappingCorrection::Update() +{ + //! Update the frame + + m_Mutex.Lock(); + + if (m_EnergyCanvas != 0) { + m_EnergyCanvas->GetCanvas()->Modified(); + m_EnergyCanvas->GetCanvas()->Update(); + } + + m_Mutex.UnLock(); +} + + +// MGUIExpoTrappingCorrection: the end... +//////////////////////////////////////////////////////////////////////////////// diff --git a/src/MGUIOptionsLoaderMeasurementsHDF.cxx b/src/MGUIOptionsLoaderMeasurementsHDF.cxx index a763f1bc..5126b5b8 100644 --- a/src/MGUIOptionsLoaderMeasurementsHDF.cxx +++ b/src/MGUIOptionsLoaderMeasurementsHDF.cxx @@ -85,12 +85,6 @@ 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(); @@ -138,7 +132,6 @@ bool MGUIOptionsLoaderMeasurementsHDF::OnApply() dynamic_cast(m_Module)->SetFileName(m_FileSelectorHDF->GetFileName()); dynamic_cast(m_Module)->SetLoadContinuationFiles(m_LoadContinuationFiles->IsOn()); dynamic_cast(m_Module)->SetFileNameStripMap(m_FileSelectorStripMap->GetFileName()); - dynamic_cast(m_Module)->SetIncludeNearestNeighbor(m_IncludeNearestNeighbor->IsOn()); return true; } 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/MModuleTrappingCorrection.cxx b/src/MModuleTrappingCorrection.cxx new file mode 100644 index 00000000..75d16306 --- /dev/null +++ b/src/MModuleTrappingCorrection.cxx @@ -0,0 +1,684 @@ +/* + * MModuleTrappingCorrection.cxx + * + * + * Copyright (C) 2008-2008 by Andreas Zoglauer, Alex Lowell, + * Sophie Haight, Carolyn Kierans. + * 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. + * + */ + + +//////////////////////////////////////////////////////////////////////////////// +// +// MModuleTrappingCorrection +// +//////////////////////////////////////////////////////////////////////////////// + + +// Include the header: +#include "MModuleTrappingCorrection.h" + +// Standard libs: + +// ROOT libs: +#include "TMath.h" +#include "TGClient.h" +#include "TH1.h" + +// MEGAlib libs: +#include "MString.h" + +// Nuclearizer libs: +#include "MGUIOptionsTrappingCorrection.h" +#include "MGUIExpoTrappingCorrection.h" +#include "MGUIExpoPlotSpectrum.h" +#include "MModuleEnergyCalibration.h" + +//////////////////////////////////////////////////////////////////////////////// + + +#ifdef ___CLING___ +ClassImp(MModuleTrappingCorrection) +#endif + + +//////////////////////////////////////////////////////////////////////////////// + + +MModuleTrappingCorrection::MModuleTrappingCorrection() : MModule() +{ + // Construct an instance of MModuleTrappingCorrection + + // Set all module relevant information + + // Set the module name --- has to be unique + m_Name = "Trapping Correction"; // - correcting energies for charge trapping (by Sophie); + + // Set the XML tag --- has to be unique --- no spaces allowed + m_XmlTag = "TrappingCorrection"; + + // Set all modules, which have to be done before this module + AddPreceedingModuleType(MAssembly::c_EnergyCalibration, true); + AddPreceedingModuleType(MAssembly::c_StripPairing, true); + AddPreceedingModuleType(MAssembly::c_TACcut, true); + AddPreceedingModuleType(MAssembly::c_EnergyCalibration, true); + AddPreceedingModuleType(MAssembly::c_DepthCorrection, true); + + // Set all types this modules handles + AddModuleType(MAssembly::c_TrappingCorrection); + + // Set all modules, which can follow this module + AddSucceedingModuleType(MAssembly::c_NoRestriction); + + // Set if this module has an options GUI + // If true, overwrite ShowOptionsGUI() with the call to the GUI! + m_HasOptionsGUI = true; + + // Allow the use of multiple threads and instances + m_AllowMultiThreading = true; + m_AllowMultipleInstances = false; + + +} + + +//////////////////////////////////////////////////////////////////////////////// + + +MModuleTrappingCorrection::~MModuleTrappingCorrection() +{ + // Delete this instance of MModuleTrappingCorrection +} + + +//////////////////////////////////////////////////////////////////////////////// + + +bool MModuleTrappingCorrection::Initialize() +{ + + // The detectors need to be in the same order as DetIDs. + // ie DetID=0 should be the 0th detector in m_Detectors, DetID=1 should the 1st, etc. + vector DetList = m_Geometry->GetDetectorList(); + + // Look through the Geometry and get the names of all the detectors. + 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. + unsigned int DetID = i; + + 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); + + if (g_Verbosity >= c_Info) { + cout << "Found detector " << det_name << " corresponding to DetID=" << DetID << "." << endl; + } + m_DetectorIDs.push_back(DetID); + m_Detectors[DetID] = det; + } else { + if (g_Verbosity >= c_Error) { + cout<<"ERROR in MModuleTrappingCorrection::Initialize: Found a duplicate detector: "<GetAvailableModuleByXmlTag("EnergyCalibration"); + if (m_EnergyCalibration == nullptr) { + cout << "MModuleTrappingCorrection: couldn't resolve pointer to Energy Calibration Module... need access to this module for energy resolution lookup!" << endl; + return false; + } + + return MModule::Initialize(); +} + + +//////////////////////////////////////////////////////////////////////////////// + +void MModuleTrappingCorrection::CreateExpos() +{ + // Create all expos + + if (HasExpos() == true) return; + + // Set the histogram display + m_ExpoSpectrum = new MGUIExpoPlotSpectrum(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 + + // Event->SetTrappingCorrectionError("GR Veto"); + return false; + + } else { + + for (unsigned int i = 0; i < Event->GetNHits(); ++i ){ + // Each event represents one photon. It contains Hits, representing interaction sites. + // H is a pointer to an instance of the MHit class. Each Hit has activated strips, represented by + // instances of the MStripHit class. + MHit* H = Event->GetHit(i); + + int Grade = GetHitGrade(H); + + // Handle different grades differently + // Get the position from the depth cal. If error is thrown, record and no depth. + + // GRADE=-1 is an error. Break from the loop and continue. + if (Grade < 0){ + H->SetNoDepth(); + } else if (Grade > 4) { // GRADE=5 is some complicated geometry with multiple hits on a single strip. GRADE=6 means not all strips are adjacent. + H->SetNoDepth(); + } else { // If the Grade is 0-4, we can handle it. + + // Take a Hit and separate its activated X- and Y-strips into separate vectors. + 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); + } + + // Get the dominant strip for the hit and its energy fraction for both LV and HV sides + double LVEnergyFraction; + double HVEnergyFraction; + MStripHit* LVSH = GetDominantStrip(LVStrips, LVEnergyFraction); + MStripHit* HVSH = GetDominantStrip(HVStrips, HVEnergyFraction); + + // Get the position value (assumed from event/hit context H) + double depth_val = H->GetPosition().GetZ(); + // double depth_val = static_cast(Zpos); + + // Correct the Low Voltage side energy if the hit pointer exists + if (LVSH != nullptr) { + double rawLVEnergy = LVSH->GetEnergy(); + double correctedLVEnergy = GetSimBasedCorrectedEnergy(depth_val, rawLVEnergy, m_CCEs_LV, m_ParamA_LV, m_ParamB_LV, m_ParamC_LV); + LVSH->SetEnergy(correctedLVEnergy); + + if (HasExpos() == true) { + m_ExpoSpectrum->AddEnergyFinal(correctedLVEnergy, LVSH->IsNearestNeighbor(), LVSH->IsLowVoltageStrip()); + } + } + + // Correct the High Voltage side energy if the hit pointer exists + if (HVSH != nullptr) { + double rawHVEnergy = HVSH->GetEnergy(); + double correctedHVEnergy = GetSimBasedCorrectedEnergy(depth_val, rawHVEnergy, m_CCEs_HV, m_ParamA_HV, m_ParamB_HV, m_ParamC_HV); + HVSH->SetEnergy(correctedHVEnergy); + + 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; + } + + // Locate the histograms dynamically from ROOT's global memory map + TH1D* histLV = (TH1D*) gDirectory->Get("EnergyHistogramLVFinal"); + TH1D* histHV = (TH1D*) gDirectory->Get("EnergyHistogramHVFinal"); + + if (histLV == nullptr) histLV = (TH1D*) gDirectory->Get("m_EnergyHistogramLVFinal"); + if (histHV == nullptr) histHV = (TH1D*) gDirectory->Get("m_EnergyHistogramHVFinal"); + + // Perform the photopeak fit for the LV peak + if (histLV != nullptr && histLV->GetEntries() > 0) { + TF1* fitFuncLV = GeneratePhotopeakFunction(); + + histLV->Fit(fitFuncLV, "RQ"); + + double mu = fitFuncLV->GetParameter("x0 (Mu)"); + double fwhm = 2.35482 * fitFuncLV->GetParameter("Sigma Gauss"); + + if (g_Verbosity >= c_Info) { + cout << m_XmlTag << " --- LV FINAL SPECTRUM FIT ---" << endl; + cout << " Centroid (Mu): " << mu << " keV | FWHM: " << fwhm << " keV" << endl; + } + delete fitFuncLV; + } + + // Perform the photopeak fit for the HV peak + if (histHV != nullptr && histHV->GetEntries() > 0) { + TF1* fitFuncHV = GeneratePhotopeakFunction(); + + histHV->Fit(fitFuncHV, "RQ"); + + double mu = fitFuncHV->GetParameter("x0 (Mu)"); + double fwhm = 2.35482 * fitFuncHV->GetParameter("Sigma Gauss"); + + if (g_Verbosity >= c_Info) { + cout << m_XmlTag << " --- HV FINAL SPECTRUM FIT ---" << endl; + cout << " Centroid (Mu): " << mu << " keV | FWHM: " << fwhm << " keV" << endl; + } + delete fitFuncHV; + } + + return; +} +///////////////////////////////////////////////////////////////////////////////// + +MStripHit* MModuleTrappingCorrection::GetDominantStrip(vector& Strips, double& EnergyFraction) +{ + double MaxEnergy = -numeric_limits::max(); // AZ: When both energies are zero (which shouldn't happen) we still pick one + 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; +} + + +///////////////////////////////////////////////////////////////////////////////// + + +bool MModuleTrappingCorrection::LoadSimCCEFile(MString FileName) +{ + MFile SimCCEFile; + if (SimCCEFile.Open(FileName) == false) { + cout << "ERROR in MModuleTrappingCorrection::LoadSimCCEFile: failed to open file." << endl; + return false; + } + + // Clear existing array data before loading new files + m_Depths.clear(); + m_CCEs_HV.clear(); + m_CCEs_LV.clear(); + + MString Line; + int ValidLineCount = 0; + + while (SimCCEFile.ReadLine(Line)) { + // Skip comment lines + if (Line.BeginsWith('#') == true) { + continue; + } + + std::vector Tokens = Line.Tokenize(","); + + // Skip empty lines safely + if (Tokens.size() == 0) { + continue; + } + + if (ValidLineCount == 0) { + // Read parameters A, B, and C from the first line + if (Tokens.size() == 6) { + m_ParamA_HV = Tokens[0].ToDouble(); + m_ParamB_HV = Tokens[1].ToDouble(); + m_ParamC_HV = Tokens[2].ToDouble(); + m_ParamA_LV = Tokens[3].ToDouble(); + m_ParamB_LV = Tokens[4].ToDouble(); + m_ParamC_LV = Tokens[5].ToDouble(); + ValidLineCount++; + } else { + cout << "ERROR in LoadSimCCEFile: Expected 6 parameters (A,B,C for HV, LV) on the first line." << endl; + SimCCEFile.Close(); + return false; + } + } + else if (ValidLineCount == 1) { + // Skip the second line which is the column header (z_depth_mm,CCE_HV) + ValidLineCount++; + } + else { + // Read the rest of the rows into your depth and CCE HV arrays + if (Tokens.size() == 3) { + m_Depths.push_back(Tokens[0].ToDouble()); + m_CCEs_HV.push_back(Tokens[1].ToDouble()); + m_CCEs_LV.push_back(Tokens[2].ToDouble()); + ValidLineCount++; + } + } + } + + SimCCEFile.Close(); + + // Print summary to console if verbose logging is enabled + if (g_Verbosity >= c_Info) { + cout << m_XmlTag << "Loaded HV parameters: A=" << m_ParamA_HV + << ", B=" << m_ParamB_HV << ", C=" << m_ParamC_HV << endl; + cout << m_XmlTag << "Loaded " << m_Depths.size() << " data points into arrays." << endl; + } + + return true; +} + + +///////////////////////////////////////////////////////////////////////////////// + +double MModuleTrappingCorrection::GetSimBasedCorrectedEnergy(double depth_val, double uncorrected_energy, const std::vector& sim_cce_sorted, double paramA, double paramB, double paramC) { + + // Look up the simulation CCE baseline using the interpolate function + double cce_base = Interpolate(depth_val, m_Depths, sim_cce_sorted); + + // Evaluate the physical trapping function model using class global popt variables + double expected_centroid_scaled = paramA * (1.0 - paramB * (1.0 - cce_base)) * (1.0 - paramC * (1.0 - cce_base)); + + // 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 continuour 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); +} + +///////////////////////////////////////////////////////////////////////////////// + + +int MModuleTrappingCorrection::GetHitGrade(MHit* H){ + // Function for choosing which Depth-to-CTD relation to use for a given event. + // At time of writing, intention is to choose a CTD based on sub-pixel region determined via charge sharing (Event "grade"). + // 5 possible grades, and one Error Grade, -1. GRADE 4 is as yet uncategorized complicated geometry. GRADE 5 means multiple, presumably separated strip hits. + + //organize x and y strips into vectors + if (H == nullptr) { + return -1; + } + if (H->GetNStripHits() == 0) { + // Error if no strip hits listed. Bad grade is returned + if (g_Verbosity >= c_Error) cout << m_XmlTag << "ERROR in MModuleTrappingCorrection: HIT WITH NO STRIP HITS" << endl; + return -1; + } + + // Take a Hit and separate its activated p and n strips into separate vectors. + std::vector LVStrips; + std::vector HVStrips; + vector LVStripIDs; + vector HVStripIDs; + for (unsigned int j = 0; j < H->GetNStripHits(); ++j) { + MStripHit* SH = H->GetStripHit(j); + if (SH == nullptr ) { + if (g_Verbosity >= c_Error) cout << m_XmlTag << "ERROR in MModuleTrappingCorrection: Trapping Correction: got NULL strip hit :( " << endl; + return -1; + } + if (SH->GetEnergy() == 0 ) { + if (g_Verbosity >= c_Error) cout << m_XmlTag << "ERROR in MModuleTrappingCorrection: Trapping Correction: got strip without energy :( " << endl; + return -1; + } + if (SH->IsLowVoltageStrip()) { + LVStrips.push_back(SH); + LVStripIDs.push_back(SH->GetStripID()); + } + else { + HVStrips.push_back(SH); + HVStripIDs.push_back(SH->GetStripID()); + } + } + + // If the same strip has multiple hits, this is a bad grade. + bool MultiHitX = H->GetStripHitMultipleTimesLV(); + bool MultiHitY = H->GetStripHitMultipleTimesHV(); + if (MultiHitX || MultiHitY) { + return 5; + } + + if (LVStrips.size()>0 && HVStrips.size()>0) { + int HVmin = * std::min_element(HVStripIDs.begin(), HVStripIDs.end()); + int HVmax = * std::max_element(HVStripIDs.begin(), HVStripIDs.end()); + + int LVmin = * std::min_element(LVStripIDs.begin(), LVStripIDs.end()); + int LVmax = * std::max_element(LVStripIDs.begin(), LVStripIDs.end()); + + // If the strip hits are not all adjacent, it's a bad grade. + if ( ((HVmax - HVmin) >= (HVStrips.size())) || ((LVmax - LVmin) >= (LVStrips.size())) ) { + return 6; + } + } + else{ + return -1; + } + + int return_value; + // If 1 strip on each side, GRADE=0 + // This represents the center of the pixel + if ( ((LVStrips.size() == 1) && (HVStrips.size() == 1)) || ((LVStrips.size() == 3) && (HVStrips.size() == 3)) ) { + return_value = 0; + } + // If 2 hits on N side and 1 on P, GRADE=1 + // This represents the middle of the edges of the pixel + else if ( (LVStrips.size() == 1) && (HVStrips.size() == 2) ) { + return_value = 1; + } + + // If 2 hits on P and 1 on N, GRADE=2 + // This represents the middle of the edges of the pixel + else if ( (LVStrips.size() == 2) && (HVStrips.size() == 1) ) { + return_value = 2; + } + + // If 2 strip hits on both sides, GRADE=3 + // This represents the corners the pixel + else if ( (LVStrips.size() == 2) && (HVStrips.size() == 2) ) { + return_value = 3; + } + + // If 3 hits on N side and 1 on P, GRADE=0 + // This represents the middle of the pixel, near the p (LV) side of the detector. + else if ( (LVStrips.size() == 1) && (HVStrips.size() == 3) ) { + return_value = 0; + } + + // If 3 hits on P and 1 on N, GRADE=0 + // This represents the middle of the pixel, near the n (HV) side of the detector. + else if ( (LVStrips.size() == 3) && (HVStrips.size() == 1) ) { + return_value = 0; + } + + // If 3 hits on N side and 2 on P, GRADE=0 + // This represents the middle of the edge of the pixel, near the p (LV) side of the detector. + else if ( (LVStrips.size() == 2) && (HVStrips.size() == 3) ) { + return_value = 2; + } + + // If 3 hits on P and 2 on N, GRADE=0 + // This represents the middle of the edge of the pixel, near the n (HV) side of the detector. + else if ( (LVStrips.size() == 3) && (HVStrips.size() == 2) ) { + return_value = 1; + } + + else { + // If more complicated than the above cases, return 4 for now. + // TODO: Handle more complicated charge distributions. + return_value = 4; + } + + return return_value; +} + +///////////////////////////////////////////////////////////////////////////////// + + +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_SimCCEFile = 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_SimCCEFile); + + 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; +} + + +//////////////////////////////////////////////////////////////////////////////// + + + +// MModuleTrappingCorrection.cxx: the end... +////////////////////////////////////////////////////////////////////////////////