From f7c71fa4ae6facc05781302bc19768969e6a35e2 Mon Sep 17 00:00:00 2001 From: Field Date: Tue, 14 Jul 2026 16:03:18 -0700 Subject: [PATCH 01/15] moved calculation of z from ctd to a function --- include/MModuleDepthCalibration.h | 5 ++ src/MModuleDepthCalibration.cxx | 140 +++++++++++++++++++++--------- 2 files changed, 103 insertions(+), 42 deletions(-) diff --git a/include/MModuleDepthCalibration.h b/include/MModuleDepthCalibration.h index 0926a61f..3e69ab14 100644 --- a/include/MModuleDepthCalibration.h +++ b/include/MModuleDepthCalibration.h @@ -21,6 +21,7 @@ #include #include #include +#include // ROOT libs: @@ -124,6 +125,10 @@ class MModuleDepthCalibration : public MModule // protected methods: protected: + + //! Returns the z as a function of ctd, given some ctd and some detector + std::tuple CalculateZfromCTD(double CTDvalue, double noise, int DetID,int Grade, bool sean_weighting); + //! Returns the strip with most energy from vector Strips, also gives back the energy fraction MStripHit* GetDominantStrip(std::vector& Strips, double& EnergyFraction); diff --git a/src/MModuleDepthCalibration.cxx b/src/MModuleDepthCalibration.cxx index adcb9044..7f761245 100644 --- a/src/MModuleDepthCalibration.cxx +++ b/src/MModuleDepthCalibration.cxx @@ -65,7 +65,7 @@ MModuleDepthCalibration::MModuleDepthCalibration() : MModule() AddPreceedingModuleType(MAssembly::c_EnergyCalibration, true); AddPreceedingModuleType(MAssembly::c_StripPairing, true); AddPreceedingModuleType(MAssembly::c_TACcut, true); -// AddPreceedingModuleType(MAssembly::c_CrosstalkCorrection, false); // Soft requirement + // AddPreceedingModuleType(MAssembly::c_CrosstalkCorrection, false); // Soft requirement // Set all types this modules handles AddModuleType(MAssembly::c_DepthCorrection); @@ -196,6 +196,7 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) // Handle different grades differently // GRADE=-1 is an error. Break from the loop and continue. + // TODO in what circumstances do we get errors? if (Grade < 0){ H->SetNoDepth(); Event->SetDepthCalibrationError("Error in depth calibration"); @@ -226,12 +227,13 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) if (SH->IsLowVoltageStrip()) LVStrips.push_back(SH); else HVStrips.push_back(SH); } + // TODO classify based on charge sharing degree. > 10% == charge sharing event! double LVEnergyFraction; double HVEnergyFraction; MStripHit* LVSH = GetDominantStrip(LVStrips, LVEnergyFraction); MStripHit* HVSH = GetDominantStrip(HVStrips, HVEnergyFraction); - double CTD_s = 0.0; + double rawCTD_s = 0.0; //now try and get z position int DetID = LVSH->GetDetectorID(); @@ -275,6 +277,7 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) double HVTiming = HVSH->GetTiming(); // If there aren't coefficients loaded, then report a depth calibration error. + // TODO check adjacent strips if so if( Coeffs == nullptr ){ // Set the bad flag for depth H->SetNoDepth(); @@ -295,58 +298,37 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) } else { // If there are coefficients and timing information is loaded, try calculating the CTD and depth - double CTD = (HVTiming - LVTiming); + // TODO FR start here + double rawCTD = (HVTiming - LVTiming); - // Confirmed that this matches SP's python code. - CTD_s = (CTD - Coeffs->at(1))/(Coeffs->at(0)); //apply inverse stretch and offset + rawCTD_s = (rawCTD - Coeffs->at(1))/(Coeffs->at(0)); //apply inverse stretch and offset double Xmin = * std::min_element(CTDVec.begin(), CTDVec.end()); double Xmax = * std::max_element(CTDVec.begin(), CTDVec.end()); - double noise = GetTimingNoiseFWHM(PixelCode, H->GetEnergy()); + double noise = GetTimingNoiseFWHM(PixelCode, H->GetEnergy()); + // TODO make the energy dependence correct.... or at least something! //if the CTD is out of range, check if we should reject the event. - if ((CTD_s < (Xmin - 2.0*noise)) || (CTD_s > (Xmax + 2.0*noise))) { + // TODO -- nope, check consistency with adjacent strips + // ALSO TODO would we want to reject tht whole event? What if this is the last hit and we still have the energy? + // also, how is slow timing dealt with? I think a lot of last-hits could be out of range for this reason... + if ((rawCTD_s < (Xmin - 2.0*noise)) || (rawCTD_s > (Xmax + 2.0*noise))) { H->SetNoDepth(); Event->SetDepthCalibrationError("Out of Range"); ++m_Error2; } // If the CTD is in range, calculate the depth - // Rather than plugging CTD into a spline to get depth, use the depth-CTD relation to calculate a probability-weighted depth value. - // This way we can avoid problems like non-monotonicity or assigning depth to events "outside" the detector - // Note that this requires that we don't massively overestimate the timing noise else { - // Calculate the probability given timing noise of CTD_s corresponding to the values of depth in DepthVec - // Utlize symmetry of the normal distribution. - vector prob_dist = norm_pdf(CTDVec, CTD_s, noise/2.355); - - // Weight the depth by probability - double prob_sum = 0.0; - for (unsigned int k=0; k < prob_dist.size(); ++k) { - prob_sum += prob_dist[k]; - } - double weighted_depth = 0.0; - - for (unsigned int k = 0; k < DepthVec.size(); ++k) { - weighted_depth += prob_dist[k] * DepthVec[k]; - } - - // Calculate the expectation value of the depth - double mean_depth = weighted_depth/prob_sum; - - // Calculate the standard deviation of the depth - double depth_var = 0.0; - - for (unsigned int k=0; k < DepthVec.size(); ++k) { - depth_var += prob_dist[k] * pow(DepthVec[k] - mean_depth, 2.0); - } - - Zsigma = sqrt(depth_var/prob_sum); - Zpos = mean_depth; - // Zpos = mean_depth - (m_Thicknesses[DetID]/2.0); - - // Add the depth to the GUI histogram. + // FR TODO the last boolean is for sean's weighting method; make it a flag + auto [rawZpos, rawZsigma] = CalculateZfromCTD(rawCTD_s, noise,DetID, Grade, false); + + // TODO depth correction loop! + Zpos = rawZpos; + Zsigma = rawZsigma; + + // Add the depth to the GUI histogram. if (Event->HasStripPairingError()==false) { if (HasExpos() == true) { m_ExpoDepthCalibration->AddDepth(DetID, Zpos); @@ -370,8 +352,6 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) H->SetPositionResolution(GlobalResolution); - - } } } @@ -381,6 +361,82 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) return true; } +///////////////////////////////////////////////////////////////////////////////// + +// TODO noise needs to be broken down into strip noise and calculated as a function of energy per strip... +std::tuple MModuleDepthCalibration::CalculateZfromCTD(double CTDvalue, double noise, int DetID,int Grade, bool sean_weighting=false) +{ + vector CTDVec = GetCTD(DetID, Grade); + vector DepthVec = GetDepth(DetID); + + + // Rather than plugging CTD into a spline to get depth, use the depth-CTD relation to calculate a probability-weighted depth value. + // This way we can avoid problems like non-monotonicity or assigning depth to events "outside" the detector + // Note that this requires that we don't massively overestimate the timing noise + if (sean_weighting){ + vector prob_dist = norm_pdf(CTDVec, CTDvalue, noise/2.355); + + // Weight the depth by probability + double prob_sum = 0.0; + for (unsigned int k=0; k < prob_dist.size(); ++k) { + prob_sum += prob_dist[k]; + } + double weighted_depth = 0.0; + + for (unsigned int k = 0; k < DepthVec.size(); ++k) { + weighted_depth += prob_dist[k] * DepthVec[k]; + } + + // Calculate the expectation value of the depth + double mean_depth = weighted_depth/prob_sum; + + // Calculate the standard deviation of the depth + double depth_var = 0.0; + + for (unsigned int k=0; k < DepthVec.size(); ++k) { + depth_var += prob_dist[k] * pow(DepthVec[k] - mean_depth, 2.0); + } + return std::make_tuple(mean_depth,sqrt(depth_var/prob_sum)); + } + // otherwise, use the standard appropach with no rounding off + // if out of bounds, return boundary + if (CTDvalue <= CTDVec.front()) { + double CTD_high = CTDvalue + noise/2.355; + auto it = std::upper_bound(CTDVec.begin(), CTDVec.end(), CTD_high); + unsigned int i = std::distance(CTDVec.begin(), it); + double fraction = (CTD_high - CTDVec[i - 1]) / (CTDVec[i] - CTDVec[i - 1]); + double depth_high = DepthVec[i - 1] + fraction * (DepthVec[i] - DepthVec[i - 1]); + return std::make_tuple(DepthVec.front(),depth_high-DepthVec.front()); + } + if (CTDvalue >= CTDVec.back()) { + double CTD_low = CTDvalue - noise/2.355; + auto it = std::upper_bound(CTDVec.begin(), CTDVec.end(), CTD_low); + unsigned int i = std::distance(CTDVec.begin(), it); + double fraction = (CTD_low - CTDVec[i - 1]) / (CTDVec[i] - CTDVec[i - 1]); + double depth_low = DepthVec[i - 1] + fraction * (DepthVec[i] - DepthVec[i - 1]); + return std::make_tuple(DepthVec.back(),depth_low-DepthVec.back()); + } + + // if not out of bounds, extrapolate and calculate errors.... + auto it = std::upper_bound(CTDVec.begin(), CTDVec.end(), CTDvalue); + unsigned int i = std::distance(CTDVec.begin(), it); + double fraction = (CTDvalue - CTDVec[i - 1]) / (CTDVec[i] - CTDVec[i - 1]); + double depth = DepthVec[i - 1] + fraction * (DepthVec[i] - DepthVec[i - 1]); + + double CTD_low = CTDvalue - noise/2.355; + double CTD_high = CTDvalue + noise/2.355; + it = std::upper_bound(CTDVec.begin(), CTDVec.end(), CTD_low); + i = std::distance(CTDVec.begin(), it); + fraction = (CTD_low - CTDVec[i - 1]) / (CTDVec[i] - CTDVec[i - 1]); + double depth_low = DepthVec[i - 1] + fraction * (DepthVec[i] - DepthVec[i - 1]); + it = std::upper_bound(CTDVec.begin(), CTDVec.end(), CTD_high); + i = std::distance(CTDVec.begin(), it); + fraction = (CTD_high - CTDVec[i - 1]) / (CTDVec[i] - CTDVec[i - 1]); + double depth_high = DepthVec[i - 1] + fraction * (DepthVec[i] - DepthVec[i - 1]); + + return std::make_tuple(depth, (depth_high - depth_low) / 2.); +} + ///////////////////////////////////////////////////////////////////////////////// From 17a3c38638ad359d18d5f7d66719b0ebef3a0a00 Mon Sep 17 00:00:00 2001 From: Field Date: Tue, 14 Jul 2026 17:19:42 -0700 Subject: [PATCH 02/15] bug fix :) --- src/MModuleDepthCalibration.cxx | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/src/MModuleDepthCalibration.cxx b/src/MModuleDepthCalibration.cxx index 7f761245..3b0168c2 100644 --- a/src/MModuleDepthCalibration.cxx +++ b/src/MModuleDepthCalibration.cxx @@ -364,7 +364,7 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) ///////////////////////////////////////////////////////////////////////////////// // TODO noise needs to be broken down into strip noise and calculated as a function of energy per strip... -std::tuple MModuleDepthCalibration::CalculateZfromCTD(double CTDvalue, double noise, int DetID,int Grade, bool sean_weighting=false) +std::tuple MModuleDepthCalibration::CalculateZfromCTD(double CTDvalue, double noise, int DetID,int Grade, bool sean_weighting) { vector CTDVec = GetCTD(DetID, Grade); vector DepthVec = GetDepth(DetID); @@ -402,6 +402,9 @@ std::tuple MModuleDepthCalibration::CalculateZfromCTD(double CTD // if out of bounds, return boundary if (CTDvalue <= CTDVec.front()) { double CTD_high = CTDvalue + noise/2.355; + if (CTD_high <= CTDVec.front()) { + return std::make_tuple(DepthVec.front(), 0.0); + } auto it = std::upper_bound(CTDVec.begin(), CTDVec.end(), CTD_high); unsigned int i = std::distance(CTDVec.begin(), it); double fraction = (CTD_high - CTDVec[i - 1]) / (CTDVec[i] - CTDVec[i - 1]); @@ -410,6 +413,9 @@ std::tuple MModuleDepthCalibration::CalculateZfromCTD(double CTD } if (CTDvalue >= CTDVec.back()) { double CTD_low = CTDvalue - noise/2.355; + if (CTD_low >= CTDVec.back()) { + return std::make_tuple(DepthVec.back(), 0.0); + } auto it = std::upper_bound(CTDVec.begin(), CTDVec.end(), CTD_low); unsigned int i = std::distance(CTDVec.begin(), it); double fraction = (CTD_low - CTDVec[i - 1]) / (CTDVec[i] - CTDVec[i - 1]); @@ -423,8 +429,8 @@ std::tuple MModuleDepthCalibration::CalculateZfromCTD(double CTD double fraction = (CTDvalue - CTDVec[i - 1]) / (CTDVec[i] - CTDVec[i - 1]); double depth = DepthVec[i - 1] + fraction * (DepthVec[i] - DepthVec[i - 1]); - double CTD_low = CTDvalue - noise/2.355; - double CTD_high = CTDvalue + noise/2.355; + double CTD_low = std::max(CTDvalue - noise/2.355,CTDVec.front()); + double CTD_high = std::min(CTDvalue + noise/2.355,CTDVec.back()); it = std::upper_bound(CTDVec.begin(), CTDVec.end(), CTD_low); i = std::distance(CTDVec.begin(), it); fraction = (CTD_low - CTDVec[i - 1]) / (CTDVec[i] - CTDVec[i - 1]); From fb8ae0de2b4bb9c819c7a549eafa934d4b1c4059 Mon Sep 17 00:00:00 2001 From: Field Date: Tue, 14 Jul 2026 17:45:21 -0700 Subject: [PATCH 03/15] turning sean method on for now as the default --- src/MModuleDepthCalibration.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/MModuleDepthCalibration.cxx b/src/MModuleDepthCalibration.cxx index 3b0168c2..1550ecbe 100644 --- a/src/MModuleDepthCalibration.cxx +++ b/src/MModuleDepthCalibration.cxx @@ -322,7 +322,7 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) // If the CTD is in range, calculate the depth else { // FR TODO the last boolean is for sean's weighting method; make it a flag - auto [rawZpos, rawZsigma] = CalculateZfromCTD(rawCTD_s, noise,DetID, Grade, false); + auto [rawZpos, rawZsigma] = CalculateZfromCTD(rawCTD_s, noise,DetID, Grade, true); // TODO depth correction loop! Zpos = rawZpos; From dad676a1707c0333b940ffd9dd09818e66f2f081 Mon Sep 17 00:00:00 2001 From: Field Date: Tue, 14 Jul 2026 18:03:21 -0700 Subject: [PATCH 04/15] made a plan to make a plan ! --- src/MModuleDepthCalibration.cxx | 27 +++++++++++++++++++++++++-- 1 file changed, 25 insertions(+), 2 deletions(-) diff --git a/src/MModuleDepthCalibration.cxx b/src/MModuleDepthCalibration.cxx index 1550ecbe..9a7b22af 100644 --- a/src/MModuleDepthCalibration.cxx +++ b/src/MModuleDepthCalibration.cxx @@ -227,7 +227,6 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) if (SH->IsLowVoltageStrip()) LVStrips.push_back(SH); else HVStrips.push_back(SH); } - // TODO classify based on charge sharing degree. > 10% == charge sharing event! double LVEnergyFraction; double HVEnergyFraction; MStripHit* LVSH = GetDominantStrip(LVStrips, LVEnergyFraction); @@ -262,6 +261,7 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) // TODO: Calculate X and Y positions more rigorously using charge sharing. + // FR note: i think this actually has to go later in the loop double Xsigma = m_YPitches[DetID]/sqrt(12.0); double Ysigma = m_XPitches[DetID]/sqrt(12.0); @@ -307,7 +307,7 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) double Xmax = * std::max_element(CTDVec.begin(), CTDVec.end()); double noise = GetTimingNoiseFWHM(PixelCode, H->GetEnergy()); - // TODO make the energy dependence correct.... or at least something! + // TODO make the energy dependence correct... and make it depend on teh two strip energies, as it should!! //if the CTD is out of range, check if we should reject the event. // TODO -- nope, check consistency with adjacent strips @@ -325,6 +325,29 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) auto [rawZpos, rawZsigma] = CalculateZfromCTD(rawCTD_s, noise,DetID, Grade, true); // TODO depth correction loop! + // + // step 1 -- check the HV side: + // -- how many strips share > 10% of the total energy (or are over slow threshold, maybe?) need this info for next steps + // -- check dTAC between adjacent strips with charge sharing and also strips relative to their low-eneryg neighbors + // -- correct the HV timing asic jitter bug, if needed, to make everything consistent + // -- TODO FLAG that HV timing asic jitter bug correction was used, if needed + // -- if charge sharing, calculate the corrected timing value for each strip in teh absense of charge sharing, + // -- and also calculate the HV tac as the weighted average, with its own uncertainty + // -- note, rawZpos is used in the above calculations! + // + // step 2 -- check the LV side: + // -- how many strips share > 10% of the total energy (or are over the slow threshold, maybe?) need this info for the next steps + // -- check dTAC between adjacent strips with charge sharing and also strips relative to their low-energy neighbors + // -- deal with the zombie bump! + // -- check for consistency between all adjacent strip pairs. If any strips with significant energy are consistent with bump: flag and correct! + // -- if charge sharing, calculate the corrected timing value for each strip in the absence of charge sharing, + // -- and then also calculate the corrected LV tac as the weighted average, with its own uncertainty (they should be consistent) + // -- note, rawZpos is used in teh above calculations + // + // step 3 -- calculate new corrected CTD, and from that calculate new corrected Z + // + // bonus points: implement x and y localization based on info in step 1 and 2 with charge sharing :) + // bonus points -- can re-calculated depth between 1 and 2 if significant changes with HV correction would impact ZB correction! Zpos = rawZpos; Zsigma = rawZsigma; From bda4fb846cfd52763f0ba6d085e0ac5c15bfb8f6 Mon Sep 17 00:00:00 2001 From: Field Date: Mon, 20 Jul 2026 14:12:38 -0700 Subject: [PATCH 05/15] added variables and functionality to enable reading in charge sharing correction config file --- include/MGUIOptionsDepthCalibration.h | 3 + include/MModuleDepthCalibration.h | 24 ++++++++ src/MGUIOptionsDepthCalibration.cxx | 9 ++- src/MModuleDepthCalibration.cxx | 89 ++++++++++++++++++++++++++- 4 files changed, 123 insertions(+), 2 deletions(-) diff --git a/include/MGUIOptionsDepthCalibration.h b/include/MGUIOptionsDepthCalibration.h index 62ca990d..8911d929 100644 --- a/include/MGUIOptionsDepthCalibration.h +++ b/include/MGUIOptionsDepthCalibration.h @@ -71,6 +71,9 @@ class MGUIOptionsDepthCalibration : public MGUIOptions //! Select which coefficients file (stretching factors and offsets) to load MGUIEFileSelector* m_CoeffsFileSelector; + //! Select which coefficients file (charge sharing correction) to load + MGUIEFileSelector* m_DtacCoeffsFileSelector; + //! Select spline file to load, splines will convert CTD->Depth MGUIEFileSelector* m_SplinesFileSelector; diff --git a/include/MModuleDepthCalibration.h b/include/MModuleDepthCalibration.h index 3e69ab14..1563b629 100644 --- a/include/MModuleDepthCalibration.h +++ b/include/MModuleDepthCalibration.h @@ -63,6 +63,12 @@ class MModuleDepthCalibration : public MModule //! Show the options GUI virtual void ShowOptionsGUI(); + //! Set filename for charge sharing coefficients file + void SetDtacCoeffsFileName( const MString& FileName) { m_DtacCoeffsFileName = FileName; } + //! Get filename for coefficients file + MString GetDtacCoeffsFileName() const { return m_DtacCoeffsFileName; } + + //! Set filename for coefficients file void SetCoeffsFileName( const MString& FileName) { m_CoeffsFileName = FileName; } //! Get filename for coefficients file @@ -92,11 +98,21 @@ class MModuleDepthCalibration : public MModule //! Load the detector and strip dimensions from the geometry object bool LoadDetectorDimensions(MDGeometryQuest* Geometry); + //! Load in the specified dTAC coefficients file + bool LoadDtacCoeffsFile(MString FName); + //! Load in the specified coefficients file bool LoadCoeffsFile(MString FName); + //! Set the charge sharing correction depth calibration coefficients + void SetDtacCoeffs( unordered_map> DtacCoeffs ) { m_DtacCoeffs = DtacCoeffs; } + //! Set the depth calibration coefficients void SetCoeffs( unordered_map> Coeffs ) { m_Coeffs = Coeffs; } + + //! Get the charge sharing correction calibration coefficients + unordered_map> GetDtacCoeffs() { return m_DtacCoeffs; } + //! Get the depth calibration coefficients unordered_map> GetCoeffs() { return m_Coeffs; } @@ -149,6 +165,9 @@ class MModuleDepthCalibration : public MModule //! Determine the Grade (geometry of charge sharing) of the Hit int GetHitGrade(MHit* H); + //! Return the charge-sharing coefficients for a strip + vector* GetDtacCoeffs(int StripPairCode); + //! Return the coefficients for a pixel vector* GetPixelCoeffs(int PixelCode); @@ -169,8 +188,12 @@ class MModuleDepthCalibration : public MModule // protected members: protected: + unordered_map> m_DtacCoeffs; // maps StripPairID to a vector of coefficients... + unordered_map> m_LVDtacPolyCoeffs; // maps DetID to the LV coefficients for dTAC vs CS map + unordered_map> m_HVDtacPolyCoeffs; // maps DetID to the HV coefficients for dTAC vs CS map unordered_map> m_Coeffs; double m_Coeffs_Energy; + MString m_DtacCoeffsFileName; MString m_CoeffsFileName; MString m_SplinesFile; unordered_map m_Thicknesses; @@ -199,6 +222,7 @@ class MModuleDepthCalibration : public MModule unordered_map> m_SplineMap; bool m_SplinesFileIsLoaded; bool m_CoeffsFileIsLoaded; + bool m_DtacCoeffsFileIsLoaded; //! The Mask Metrology file name MString m_MaskMetrologyFileName; diff --git a/src/MGUIOptionsDepthCalibration.cxx b/src/MGUIOptionsDepthCalibration.cxx index 7da46225..097f9fcb 100644 --- a/src/MGUIOptionsDepthCalibration.cxx +++ b/src/MGUIOptionsDepthCalibration.cxx @@ -66,7 +66,7 @@ void MGUIOptionsDepthCalibration::Create() { PreCreate(); - m_CoeffsFileSelector = new MGUIEFileSelector(m_OptionsFrame, "Select a coefficients file:", + m_CoeffsFileSelector = new MGUIEFileSelector(m_OptionsFrame, "Select a depth calcoefficients file:", dynamic_cast(m_Module)->GetCoeffsFileName()); m_CoeffsFileSelector->SetFileType("coeffs", "*.csv"); TGLayoutHints* LabelLayout = new TGLayoutHints(kLHintsTop | kLHintsCenterX | kLHintsExpandX, 10, 10, 10, 10); @@ -78,6 +78,12 @@ void MGUIOptionsDepthCalibration::Create() // TGLayoutHints* Label2Layout = new TGLayoutHints(kLHintsTop | kLHintsCenterX | kLHintsExpandX, 10, 10, 10, 10); m_OptionsFrame->AddFrame(m_SplinesFileSelector, LabelLayout); + m_DtacCoeffsFileSelector = new MGUIEFileSelector(m_OptionsFrame, "Select a dtac (charge sharing) coefficients file:", + dynamic_cast(m_Module)->GetDtacCoeffsFileName()); + m_DtacCoeffsFileSelector->SetFileType("dtacCoeffs", "*.csv"); +// TGLayoutHints* Label2Layout = new TGLayoutHints(kLHintsTop | kLHintsCenterX | kLHintsExpandX, 10, 10, 10, 10); + m_OptionsFrame->AddFrame(m_DtacCoeffsFileSelector, LabelLayout); + m_MaskMetModeCB = new TGCheckButton(m_OptionsFrame, "Enable mask metrology correction and read calibration from file:", c_MetrologyFile); m_MaskMetModeCB->SetState((dynamic_cast(m_Module)->GetMaskMetrologyCorrectionEnable() == true) ? kButtonDown : kButtonUp); m_MaskMetModeCB->Associate(this); @@ -159,6 +165,7 @@ bool MGUIOptionsDepthCalibration::OnApply() { // Modify this to store the data in the module! + dynamic_cast(m_Module)->SetDtacCoeffsFileName(m_DtacCoeffsFileSelector->GetFileName()); dynamic_cast(m_Module)->SetCoeffsFileName(m_CoeffsFileSelector->GetFileName()); dynamic_cast(m_Module)->SetSplinesFileName(m_SplinesFileSelector->GetFileName()); diff --git a/src/MModuleDepthCalibration.cxx b/src/MModuleDepthCalibration.cxx index 9a7b22af..6277653c 100644 --- a/src/MModuleDepthCalibration.cxx +++ b/src/MModuleDepthCalibration.cxx @@ -127,6 +127,10 @@ bool MModuleDepthCalibration::Initialize() if (m_CoeffsFileIsLoaded == false) { return false; } + m_DtacCoeffsFileIsLoaded = LoadDtacCoeffsFile(m_DtacCoeffsFileName); + if (m_DtacCoeffsFileIsLoaded == false) { + return false; + } m_SplinesFileIsLoaded = LoadSplinesFile(m_SplinesFile); if (m_SplinesFileIsLoaded == false) { return false; @@ -327,6 +331,7 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) // TODO depth correction loop! // // step 1 -- check the HV side: + int StripPairCode = 10000*DetID + HVStripID; // note, it is the lower strip ID always (eg, 15 if sharing between strip 15 and 16) // -- how many strips share > 10% of the total energy (or are over slow threshold, maybe?) need this info for next steps // -- check dTAC between adjacent strips with charge sharing and also strips relative to their low-eneryg neighbors // -- correct the HV timing asic jitter bug, if needed, to make everything consistent @@ -336,6 +341,7 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) // -- note, rawZpos is used in the above calculations! // // step 2 -- check the LV side: + StripPairCode = 10000*DetID + 100*LVStripID; // note, it is the lower Strip ID always! // -- how many strips share > 10% of the total energy (or are over the slow threshold, maybe?) need this info for the next steps // -- check dTAC between adjacent strips with charge sharing and also strips relative to their low-energy neighbors // -- deal with the zombie bump! @@ -575,6 +581,58 @@ bool MModuleDepthCalibration::LoadDetectorDimensions(MDGeometryQuest* Geometry) return true; } +bool MModuleDepthCalibration::LoadDtacCoeffsFile(MString FileName) +{ + // Read in the dTAC coefficients file, which gives the coefficients needed for per-strip charge sharing correction + // it should have a header line with the following info: + // TODO info here once finalized! + // it should contain for each pixel + // TODO final form here!! + + // TODO replace this! temp fix while waiting for actual config file... + for (int DetID = 0; DetID < 1; DetID++){ + vector poly_coeffs; + poly_coeffs.push_back(139.209); poly_coeffs.push_back(-106.849);// these are the coefficients, where we'll have a polynomial dTac = coeffs[0]*(f-0.5) - coeffs[1]*(f-0.5)^3 + m_LVDtacPolyCoeffs[DetID] = poly_coeffs; + m_HVDtacPolyCoeffs[DetID] = poly_coeffs; // in principle they would be different + + vector coeffs; // the stretch and offset, currently set to the same values for all strips (which are almost certainly wrong) + coeffs.push_back(1.033454449); coeffs.push_back(-1.996517705);coeffs.push_back(6.373983606); // stretch, offset, dTacSigma. Need to figure out how to deal with dTac sigma in an energy-dependent way + for (int LVStripID = 0; LVStripID < 63; LVStripID++){ // up to 62 since these are pairs + int StripPairCode = 10000*DetID + 100*LVStripID; + m_DtacCoeffs[StripPairCode] = coeffs; + } + for (int HVStripID = 0; HVStripID < 63; HVStripID++){ + int StripPairCode = 10000*DetID + HVStripID; + m_DtacCoeffs[StripPairCode] = coeffs; + } + } + return true; + + //read in file + MFile DtacCoeffsFile; + std::vector HeaderTokens; + if (DtacCoeffsFile.Open(FileName) == false) { + cout << "ERROR in MModuleDepthCalibration::LoadDtacCoeffsFile: failed to open dtac coefficients file." < Tokens = Line.Tokenize(","); + vector coeffs; + + } + } + DtacCoeffsFile.Close(); + return true; +} bool MModuleDepthCalibration::LoadCoeffsFile(MString FileName) { @@ -592,7 +650,7 @@ bool MModuleDepthCalibration::LoadCoeffsFile(MString FileName) MString Line; while (CoeffsFile.ReadLine(Line) == true) { if (Line.BeginsWith('#') == true) { - std::vector Tokens = Line.Tokenize(" "); + std::vector Tokens = Line.Tokenize(" "); // TODO why is it like this? Did i not put commas in the header of the depth cal? m_Coeffs_Energy = Tokens[5].ToDouble(); if (g_Verbosity >= c_Info) { cout << m_XmlTag << "The stretch and offset were calculated for " << m_Coeffs_Energy << " keV." << endl; @@ -623,6 +681,28 @@ bool MModuleDepthCalibration::LoadCoeffsFile(MString FileName) ///////////////////////////////////////////////////////////////////////////////// +std::vector* MModuleDepthCalibration::GetDtacCoeffs(int StripPairCode) +{ + // Check to see if the charge sharing coefficients have been loaded. If so, try to get the coefficients for the specified strip pair. + if (m_DtacCoeffsFileIsLoaded == true) { + if (m_DtacCoeffs.count(StripPairCode) > 0) { + return &m_Coeffs[StripPairCode]; + } else { + if (g_Verbosity >= c_Warning) { + cout << "MModuleDepthCalibration::GetDtacCoeffs: cannot get charge sharing coefficients; strip pair code " << StripPairCode << " not found." << endl; + } + return nullptr; + } + } else { + cout << "MModuleDepthCalibration::GetDtacCoeffs: cannot get charge sharing coefficients; file has not yet been loaded." << endl; + return nullptr; + } + +} + +///////////////////////////////////////////////////////////////////////////////// + + std::vector* MModuleDepthCalibration::GetPixelCoeffs(int PixelCode) { // Check to see if the stretch and offset have been loaded. If so, try to get the coefficients for the specified pixel. @@ -1142,6 +1222,11 @@ bool MModuleDepthCalibration::ReadXmlConfiguration(MXmlNode* Node) m_CoeffsFileName = CoeffsFileNameNode->GetValue(); } + MXmlNode* DtacCoeffsFileNameNode = Node->GetNode("DtacCoeffsFileName"); + if (DtacCoeffsFileNameNode != nullptr) { + m_DtacCoeffsFileName = DtacCoeffsFileNameNode->GetValue(); + } + MXmlNode* SplinesFileNameNode = Node->GetNode("SplinesFileName"); if (SplinesFileNameNode != nullptr) { m_SplinesFile = SplinesFileNameNode->GetValue(); @@ -1174,6 +1259,7 @@ MXmlNode* MModuleDepthCalibration::CreateXmlConfiguration() MXmlNode* Node = new MXmlNode(0,m_XmlTag); new MXmlNode(Node, "CoeffsFileName", m_CoeffsFileName); + new MXmlNode(Node, "DtacCoeffsFileName", m_DtacCoeffsFileName); new MXmlNode(Node, "SplinesFileName", m_SplinesFile); new MXmlNode(Node, "MaskMetrology", (bool)m_MaskMetrologyEnabled); new MXmlNode(Node, "MaskMetrologyFileName", m_MaskMetrologyFileName); @@ -1202,6 +1288,7 @@ void MModuleDepthCalibration::Finalize() // Clean up maps and vectors m_Coeffs.clear(); + m_DtacCoeffs.clear(); m_Thicknesses.clear(); m_NXStrips.clear(); m_NYStrips.clear(); From 953bf8d80e1b4d569acda74c06fca0d6d34f5c34 Mon Sep 17 00:00:00 2001 From: Field Date: Mon, 20 Jul 2026 18:14:41 -0700 Subject: [PATCH 06/15] updating to enable mulitple parameters and polynomial coefficients in the charge sharing conf, for different depths --- include/MGUIOptionsDepthCalibration.h | 2 +- include/MModuleDepthCalibration.h | 62 ++++++---- src/MGUIOptionsDepthCalibration.cxx | 10 +- src/MModuleDepthCalibration.cxx | 162 ++++++++++++++++++++------ 4 files changed, 176 insertions(+), 60 deletions(-) diff --git a/include/MGUIOptionsDepthCalibration.h b/include/MGUIOptionsDepthCalibration.h index 8911d929..dc82d192 100644 --- a/include/MGUIOptionsDepthCalibration.h +++ b/include/MGUIOptionsDepthCalibration.h @@ -72,7 +72,7 @@ class MGUIOptionsDepthCalibration : public MGUIOptions MGUIEFileSelector* m_CoeffsFileSelector; //! Select which coefficients file (charge sharing correction) to load - MGUIEFileSelector* m_DtacCoeffsFileSelector; + MGUIEFileSelector* m_ChargeSharingConfigFileSelector; //! Select spline file to load, splines will convert CTD->Depth MGUIEFileSelector* m_SplinesFileSelector; diff --git a/include/MModuleDepthCalibration.h b/include/MModuleDepthCalibration.h index 1563b629..ad88cba3 100644 --- a/include/MModuleDepthCalibration.h +++ b/include/MModuleDepthCalibration.h @@ -64,12 +64,11 @@ class MModuleDepthCalibration : public MModule virtual void ShowOptionsGUI(); //! Set filename for charge sharing coefficients file - void SetDtacCoeffsFileName( const MString& FileName) { m_DtacCoeffsFileName = FileName; } + void SetChargeSharingConfigFileName( const MString& FileName) { m_ChargeSharingConfigFileName = FileName; } //! Get filename for coefficients file - MString GetDtacCoeffsFileName() const { return m_DtacCoeffsFileName; } + MString GetChargeSharingConfigFileName() const { return m_ChargeSharingConfigFileName; } - - //! Set filename for coefficients file + //! Set filename for ctd-z calibration coefficients file void SetCoeffsFileName( const MString& FileName) { m_CoeffsFileName = FileName; } //! Get filename for coefficients file MString GetCoeffsFileName() const { return m_CoeffsFileName; } @@ -98,20 +97,34 @@ class MModuleDepthCalibration : public MModule //! Load the detector and strip dimensions from the geometry object bool LoadDetectorDimensions(MDGeometryQuest* Geometry); - //! Load in the specified dTAC coefficients file - bool LoadDtacCoeffsFile(MString FName); + //! Load in the specified charge sharing calibration coefficients file + bool LoadChargeSharingConfigFile(MString FName); - //! Load in the specified coefficients file + //! Load in the specified ctd-z coefficients file bool LoadCoeffsFile(MString FName); - //! Set the charge sharing correction depth calibration coefficients - void SetDtacCoeffs( unordered_map> DtacCoeffs ) { m_DtacCoeffs = DtacCoeffs; } + //! Set the vector of z values for each detector (which depths are in the config for the charge sharing correction for each detector) + void SetChargeSharingDepths( unordered_map> ChargeSharingDepths) {m_ChargeSharingDepths = ChargeSharingDepths;} + + //! Set the charge sharing correction depth calibration coefficient, one vector per strip-pair (map-int) per depth (the vector) + void SetChargeSharingCoeffs( unordered_map>> ChargeSharingCoeffs ) { m_ChargeSharingCoeffs = ChargeSharingCoeffs; } + + //! Set the coefficients of the polynomial describing the charge sharing correction depth calibration coefficient, one vector per detector (map-int) per depth (the vector) + void SetChargeSharingPolyCoeffsHV( unordered_map>> ChargeSharingPolyCoeffsHV ) { m_ChargeSharingPolyCoeffsHV = ChargeSharingPolyCoeffsHV; } + void SetChargeSharingPolyCoeffsLV( unordered_map>> ChargeSharingPolyCoeffsLV ) { m_ChargeSharingPolyCoeffsLV = ChargeSharingPolyCoeffsLV; } //! Set the depth calibration coefficients void SetCoeffs( unordered_map> Coeffs ) { m_Coeffs = Coeffs; } + //! Get the depths for the charge sharing correction calibration coefficients for a detector + unordered_map> GetChargeSharingDepths() { return m_ChargeSharingDepths; } + //! Get the charge sharing correction calibration coefficients - unordered_map> GetDtacCoeffs() { return m_DtacCoeffs; } + unordered_map>> GetChargeSharingCoeffs() { return m_ChargeSharingCoeffs; } + + //! Get the coefficients of the polynomial for the ccharge sharing correction calibration + unordered_map>> GetChargeSharingPolyCoeffsHV() { return m_ChargeSharingPolyCoeffsHV; } + unordered_map>> GetChargeSharingPolyCoeffsLV() { return m_ChargeSharingPolyCoeffsLV; } //! Get the depth calibration coefficients unordered_map> GetCoeffs() { return m_Coeffs; } @@ -165,10 +178,14 @@ class MModuleDepthCalibration : public MModule //! Determine the Grade (geometry of charge sharing) of the Hit int GetHitGrade(MHit* H); - //! Return the charge-sharing coefficients for a strip - vector* GetDtacCoeffs(int StripPairCode); + //! Return the coefficients of the dTac polynomial for a detector at a depth + vector* GetChargeSharingPolyCoeffsHV(int DetID,double z); + vector* GetChargeSharingPolyCoeffsLV(int DetID,double z); + + //! Return the charge-sharing stretch/offset coefficients for a strip-pair at a depth + vector* GetChargeSharingCoeffs(int StripPairCode,double z); - //! Return the coefficients for a pixel + //! Return the ctd-z stretch/offset coefficients for a pixel vector* GetPixelCoeffs(int PixelCode); //! Load the metrology mask file @@ -177,6 +194,7 @@ class MModuleDepthCalibration : public MModule //! Get the x, y position of the intersection of two strips based on the Metrology Mask vector GetStripIntersection(MReadOutElementDoubleStrip LVStrip, MReadOutElementDoubleStrip HVStrip); + // TODO this should require strip energy, not hit energy... //! Get the timing FWHM noise for the specified pixel and Energy double GetTimingNoiseFWHM(int PixelCode, double Energy); @@ -186,14 +204,16 @@ class MModuleDepthCalibration : public MModule // protected members: - protected: - - unordered_map> m_DtacCoeffs; // maps StripPairID to a vector of coefficients... - unordered_map> m_LVDtacPolyCoeffs; // maps DetID to the LV coefficients for dTAC vs CS map - unordered_map> m_HVDtacPolyCoeffs; // maps DetID to the HV coefficients for dTAC vs CS map - unordered_map> m_Coeffs; + protected: + + unordered_map> m_ChargeSharingDepths; // maps DetID to the depths for which the charge sharing polynomial correction and charge sharing coefficients per-strip-pair were calculated + unordered_map>> m_ChargeSharingCoeffs; // maps StripPairID to a vector of coefficients, for a vector of depths (needs interpolation) + vector m_InterpolatedCoeffs; // holder for interpolated charge sharing coeffs between different depths + unordered_map>> m_ChargeSharingPolyCoeffsLV; // maps DetID to the LV coefficients for dTAC vs CS map and charge sharing correction vs CS map, for a vector of depths + unordered_map>> m_ChargeSharingPolyCoeffsHV; // maps DetID to the HV coefficients for dTAC vs CS map and charge sharing correction vs CS map, for a given depth + unordered_map> m_Coeffs; // maps pix id to a vector of coefficients... double m_Coeffs_Energy; - MString m_DtacCoeffsFileName; + MString m_ChargeSharingConfigFileName; MString m_CoeffsFileName; MString m_SplinesFile; unordered_map m_Thicknesses; @@ -222,7 +242,7 @@ class MModuleDepthCalibration : public MModule unordered_map> m_SplineMap; bool m_SplinesFileIsLoaded; bool m_CoeffsFileIsLoaded; - bool m_DtacCoeffsFileIsLoaded; + bool m_ChargeSharingConfigFileIsLoaded; //! The Mask Metrology file name MString m_MaskMetrologyFileName; diff --git a/src/MGUIOptionsDepthCalibration.cxx b/src/MGUIOptionsDepthCalibration.cxx index 097f9fcb..81c3ceec 100644 --- a/src/MGUIOptionsDepthCalibration.cxx +++ b/src/MGUIOptionsDepthCalibration.cxx @@ -78,11 +78,11 @@ void MGUIOptionsDepthCalibration::Create() // TGLayoutHints* Label2Layout = new TGLayoutHints(kLHintsTop | kLHintsCenterX | kLHintsExpandX, 10, 10, 10, 10); m_OptionsFrame->AddFrame(m_SplinesFileSelector, LabelLayout); - m_DtacCoeffsFileSelector = new MGUIEFileSelector(m_OptionsFrame, "Select a dtac (charge sharing) coefficients file:", - dynamic_cast(m_Module)->GetDtacCoeffsFileName()); - m_DtacCoeffsFileSelector->SetFileType("dtacCoeffs", "*.csv"); + m_ChargeSharingConfigFileSelector = new MGUIEFileSelector(m_OptionsFrame, "Select a dtac (charge sharing) coefficients file:", + dynamic_cast(m_Module)->GetChargeSharingConfigFileName()); + m_ChargeSharingConfigFileSelector->SetFileType("Charge Sharing Config", "*.csv"); // TGLayoutHints* Label2Layout = new TGLayoutHints(kLHintsTop | kLHintsCenterX | kLHintsExpandX, 10, 10, 10, 10); - m_OptionsFrame->AddFrame(m_DtacCoeffsFileSelector, LabelLayout); + m_OptionsFrame->AddFrame(m_ChargeSharingConfigFileSelector, LabelLayout); m_MaskMetModeCB = new TGCheckButton(m_OptionsFrame, "Enable mask metrology correction and read calibration from file:", c_MetrologyFile); m_MaskMetModeCB->SetState((dynamic_cast(m_Module)->GetMaskMetrologyCorrectionEnable() == true) ? kButtonDown : kButtonUp); @@ -165,7 +165,7 @@ bool MGUIOptionsDepthCalibration::OnApply() { // Modify this to store the data in the module! - dynamic_cast(m_Module)->SetDtacCoeffsFileName(m_DtacCoeffsFileSelector->GetFileName()); + dynamic_cast(m_Module)->SetChargeSharingConfigFileName(m_ChargeSharingConfigFileSelector->GetFileName()); dynamic_cast(m_Module)->SetCoeffsFileName(m_CoeffsFileSelector->GetFileName()); dynamic_cast(m_Module)->SetSplinesFileName(m_SplinesFileSelector->GetFileName()); diff --git a/src/MModuleDepthCalibration.cxx b/src/MModuleDepthCalibration.cxx index 6277653c..c95ed410 100644 --- a/src/MModuleDepthCalibration.cxx +++ b/src/MModuleDepthCalibration.cxx @@ -127,8 +127,8 @@ bool MModuleDepthCalibration::Initialize() if (m_CoeffsFileIsLoaded == false) { return false; } - m_DtacCoeffsFileIsLoaded = LoadDtacCoeffsFile(m_DtacCoeffsFileName); - if (m_DtacCoeffsFileIsLoaded == false) { + m_ChargeSharingConfigFileIsLoaded = LoadChargeSharingConfigFile(m_ChargeSharingConfigFileName); + if (m_ChargeSharingConfigFileIsLoaded == false) { return false; } m_SplinesFileIsLoaded = LoadSplinesFile(m_SplinesFile); @@ -508,6 +508,7 @@ double MModuleDepthCalibration::GetTimingNoiseFWHM(int PixelCode, double Energy) // Placeholder for determining the timing noise with energy, and possibly even on a pixel-by-pixel basis. // Should follow 1/E relation // TODO: Determine real energy dependence and implement it here. + // TODO: should be a function of strip, not pixel double noiseFWHM = 0.0; if (m_CoeffsFileIsLoaded == true) { noiseFWHM = m_Coeffs[PixelCode][2] * m_Coeffs_Energy/Energy; @@ -581,7 +582,7 @@ bool MModuleDepthCalibration::LoadDetectorDimensions(MDGeometryQuest* Geometry) return true; } -bool MModuleDepthCalibration::LoadDtacCoeffsFile(MString FileName) +bool MModuleDepthCalibration::LoadChargeSharingConfigFile(MString FileName) { // Read in the dTAC coefficients file, which gives the coefficients needed for per-strip charge sharing correction // it should have a header line with the following info: @@ -591,34 +592,44 @@ bool MModuleDepthCalibration::LoadDtacCoeffsFile(MString FileName) // TODO replace this! temp fix while waiting for actual config file... for (int DetID = 0; DetID < 1; DetID++){ - vector poly_coeffs; - poly_coeffs.push_back(139.209); poly_coeffs.push_back(-106.849);// these are the coefficients, where we'll have a polynomial dTac = coeffs[0]*(f-0.5) - coeffs[1]*(f-0.5)^3 - m_LVDtacPolyCoeffs[DetID] = poly_coeffs; - m_HVDtacPolyCoeffs[DetID] = poly_coeffs; // in principle they would be different - - vector coeffs; // the stretch and offset, currently set to the same values for all strips (which are almost certainly wrong) - coeffs.push_back(1.033454449); coeffs.push_back(-1.996517705);coeffs.push_back(6.373983606); // stretch, offset, dTacSigma. Need to figure out how to deal with dTac sigma in an energy-dependent way - for (int LVStripID = 0; LVStripID < 63; LVStripID++){ // up to 62 since these are pairs - int StripPairCode = 10000*DetID + 100*LVStripID; - m_DtacCoeffs[StripPairCode] = coeffs; - } - for (int HVStripID = 0; HVStripID < 63; HVStripID++){ - int StripPairCode = 10000*DetID + HVStripID; - m_DtacCoeffs[StripPairCode] = coeffs; + + // fill m_ChargeSharingDepths for each detector + vector depths; + for (double i = -7.; i < 7.4; i = i + 1.) depths.push_back(i); + m_ChargeSharingDepths[DetID] = depths; + + for(int z = 0; z < depths.size(); z++){ + // fill m_ChargeSharingPolyCoeffs (HV and LV) for each detector and each depth + vector poly_coeffs; + poly_coeffs.push_back(139.209); poly_coeffs.push_back(-106.849);// these are the coefficients, where we'll have a polynomial dTac = coeffs[0]*(f-0.5) - coeffs[1]*(f-0.5)^3 + m_ChargeSharingPolyCoeffsHV[DetID].push_back(poly_coeffs); + m_ChargeSharingPolyCoeffsLV[DetID].push_back(poly_coeffs); // in principle they would be different + + // fill m_ChargeSharingConfig for each detector / depth / strip pair + vector coeffs; // the stretch and offset, currently set to the same values for all strips (which are almost certainly wrong) + coeffs.push_back(1.033454449); coeffs.push_back(-1.996517705);coeffs.push_back(6.373983606); // stretch, offset, dTacSigma. Need to figure out how to deal with dTac sigma in an energy-dependent way + for (int LVStripID = 0; LVStripID < 63; LVStripID++){ // up to 62 since these are pairs + int StripPairCode = 10000*DetID + 100*LVStripID; + m_ChargeSharingCoeffs[StripPairCode].push_back(coeffs); + } + for (int HVStripID = 0; HVStripID < 63; HVStripID++){ + int StripPairCode = 10000*DetID + HVStripID; + m_ChargeSharingCoeffs[StripPairCode].push_back(coeffs); + } } } return true; //read in file - MFile DtacCoeffsFile; + MFile ChargeSharingConfigFile; std::vector HeaderTokens; - if (DtacCoeffsFile.Open(FileName) == false) { - cout << "ERROR in MModuleDepthCalibration::LoadDtacCoeffsFile: failed to open dtac coefficients file." <* MModuleDepthCalibration::GetDtacCoeffs(int StripPairCode) +std::vector* MModuleDepthCalibration::GetChargeSharingCoeffs(int StripPairCode, double z) { + int DetID = StripPairCode / 10000; // Check to see if the charge sharing coefficients have been loaded. If so, try to get the coefficients for the specified strip pair. - if (m_DtacCoeffsFileIsLoaded == true) { - if (m_DtacCoeffs.count(StripPairCode) > 0) { - return &m_Coeffs[StripPairCode]; + if (m_ChargeSharingConfigFileIsLoaded == true) { + if (m_ChargeSharingCoeffs.count(StripPairCode) > 0) { + + // if we only sampled one depth, or we're beyond the depth range, just return the closest coefficients + if (z <= m_ChargeSharingDepths[DetID].front()) return &m_ChargeSharingCoeffs[StripPairCode].at(0); + if (z >= m_ChargeSharingDepths[DetID].back()) return &m_ChargeSharingCoeffs[StripPairCode].at(m_ChargeSharingDepths[DetID].size()-1); + + // otherwise, interpolate + for (unsigned int i = 0; i < m_ChargeSharingDepths[DetID].size() - 1; i++){ + if (z >= m_ChargeSharingDepths[DetID].at(i) && z < m_ChargeSharingDepths[DetID].at(i + 1)) { + double f = (z - m_ChargeSharingDepths[DetID][i]) / (m_ChargeSharingDepths[DetID][i + 1] - m_ChargeSharingDepths[DetID][i]); + m_InterpolatedCoeffs.clear(); + for (int j = 0; j < m_ChargeSharingCoeffs[StripPairCode].at(i).size(); j++) m_InterpolatedCoeffs.push_back((1.0 - f) * m_ChargeSharingCoeffs[StripPairCode][i][j] + f * m_ChargeSharingCoeffs[StripPairCode][i + 1][j]); + return &m_InterpolatedCoeffs; + } + } + } else { + if (g_Verbosity >= c_Warning) { + cout << "MModuleDepthCalibration::GetChargeSharingCoeffs: cannot get charge sharing coefficients; strip pair code " << StripPairCode << " not found." << endl; + } + return nullptr; + } + } else { + cout << "MModuleDepthCalibration::GetChargeSharingCoeffs: cannot get charge sharing coefficients; file has not yet been loaded." << endl; + return nullptr; + } + +} +///////////////////////////////////////////////////////////////////////////////// + + +std::vector* MModuleDepthCalibration::GetChargeSharingPolyCoeffsLV(int DetID, double z) +{ + // Check to see if the charge sharing coefficients have been loaded. If so, try to get the coefficients for the specified strip pair. + if (m_ChargeSharingConfigFileIsLoaded == true) { + if (m_ChargeSharingPolyCoeffsLV.count(DetID) > 0) { + + // if we only sampled one depth, or we're beyond the depth range, just return the closest coefficients + if (z <= m_ChargeSharingDepths[DetID].front()) return &m_ChargeSharingPolyCoeffsLV[DetID].at(0); + if (z >= m_ChargeSharingDepths[DetID].back()) return &m_ChargeSharingPolyCoeffsLV[DetID].at(m_ChargeSharingDepths[DetID].size()-1); + + // otherwise, interpolate + for (unsigned int i = 0; i < m_ChargeSharingDepths[DetID].size() - 1; i++){ + if (z >= m_ChargeSharingDepths[DetID].at(i) && z < m_ChargeSharingDepths[DetID].at(i + 1)) { + double f = (z - m_ChargeSharingDepths[DetID][i]) / (m_ChargeSharingDepths[DetID][i + 1] - m_ChargeSharingDepths[DetID][i]); + m_InterpolatedCoeffs.clear(); + for (int j = 0; j < m_ChargeSharingPolyCoeffsLV[DetID].at(i).size(); j++) m_InterpolatedCoeffs.push_back((1.0 - f) * m_ChargeSharingPolyCoeffsLV[DetID][i][j] + f * m_ChargeSharingPolyCoeffsLV[DetID][i + 1][j]); + return &m_InterpolatedCoeffs; + } + } + } else { + if (g_Verbosity >= c_Warning) { + cout << "MModuleDepthCalibration::GetChargeSharingPolyCoeffsLV: cannot get charge sharing polynomial coefficients; detector id code " << DetID << " not found." << endl; + } + return nullptr; + } + } else { + cout << "MModuleDepthCalibration::GetChargeSharingPolyCoeffsLV: cannot get charge sharing coefficients; file has not yet been loaded." << endl; + return nullptr; + } + +} +///////////////////////////////////////////////////////////////////////////////// + + +std::vector* MModuleDepthCalibration::GetChargeSharingPolyCoeffsHV(int DetID, double z) +{ + // Check to see if the charge sharing coefficients have been loaded. If so, try to get the coefficients for the specified strip pair. + if (m_ChargeSharingConfigFileIsLoaded == true) { + if (m_ChargeSharingPolyCoeffsHV.count(DetID) > 0) { + + // if we only sampled one depth, or we're beyond the depth range, just return the closest coefficients + if (z <= m_ChargeSharingDepths[DetID].front()) return &m_ChargeSharingPolyCoeffsHV[DetID].at(0); + if (z >= m_ChargeSharingDepths[DetID].back()) return &m_ChargeSharingPolyCoeffsHV[DetID].at(m_ChargeSharingDepths[DetID].size()-1); + + // otherwise, interpolate + for (unsigned int i = 0; i < m_ChargeSharingDepths[DetID].size() - 1; i++){ + if (z >= m_ChargeSharingDepths[DetID].at(i) && z < m_ChargeSharingDepths[DetID].at(i + 1)) { + double f = (z - m_ChargeSharingDepths[DetID][i]) / (m_ChargeSharingDepths[DetID][i + 1] - m_ChargeSharingDepths[DetID][i]); + m_InterpolatedCoeffs.clear(); + for (int j = 0; j < m_ChargeSharingPolyCoeffsHV[DetID].at(i).size(); j++) m_InterpolatedCoeffs.push_back((1.0 - f) * m_ChargeSharingPolyCoeffsHV[DetID][i][j] + f * m_ChargeSharingPolyCoeffsHV[DetID][i + 1][j]); + return &m_InterpolatedCoeffs; + } + } } else { if (g_Verbosity >= c_Warning) { - cout << "MModuleDepthCalibration::GetDtacCoeffs: cannot get charge sharing coefficients; strip pair code " << StripPairCode << " not found." << endl; + cout << "MModuleDepthCalibration::GetChargeSharingPolyCoeffsHV: cannot get charge sharing polynomial coefficients; detector id code " << DetID << " not found." << endl; } return nullptr; } } else { - cout << "MModuleDepthCalibration::GetDtacCoeffs: cannot get charge sharing coefficients; file has not yet been loaded." << endl; + cout << "MModuleDepthCalibration::GetChargeSharingPolyCoeffsHV: cannot get charge sharing coefficients; file has not yet been loaded." << endl; return nullptr; } @@ -1222,9 +1315,9 @@ bool MModuleDepthCalibration::ReadXmlConfiguration(MXmlNode* Node) m_CoeffsFileName = CoeffsFileNameNode->GetValue(); } - MXmlNode* DtacCoeffsFileNameNode = Node->GetNode("DtacCoeffsFileName"); - if (DtacCoeffsFileNameNode != nullptr) { - m_DtacCoeffsFileName = DtacCoeffsFileNameNode->GetValue(); + MXmlNode* ChargeSharingConfigFileNameNode = Node->GetNode("ChargeSharingConfigFileName"); + if (ChargeSharingConfigFileNameNode != nullptr) { + m_ChargeSharingConfigFileName = ChargeSharingConfigFileNameNode->GetValue(); } MXmlNode* SplinesFileNameNode = Node->GetNode("SplinesFileName"); @@ -1259,7 +1352,7 @@ MXmlNode* MModuleDepthCalibration::CreateXmlConfiguration() MXmlNode* Node = new MXmlNode(0,m_XmlTag); new MXmlNode(Node, "CoeffsFileName", m_CoeffsFileName); - new MXmlNode(Node, "DtacCoeffsFileName", m_DtacCoeffsFileName); + new MXmlNode(Node, "ChargeSharingConfigFileName", m_ChargeSharingConfigFileName); new MXmlNode(Node, "SplinesFileName", m_SplinesFile); new MXmlNode(Node, "MaskMetrology", (bool)m_MaskMetrologyEnabled); new MXmlNode(Node, "MaskMetrologyFileName", m_MaskMetrologyFileName); @@ -1288,7 +1381,10 @@ void MModuleDepthCalibration::Finalize() // Clean up maps and vectors m_Coeffs.clear(); - m_DtacCoeffs.clear(); + m_ChargeSharingCoeffs.clear(); + m_ChargeSharingPolyCoeffsHV.clear(); + m_ChargeSharingPolyCoeffsLV.clear(); + m_ChargeSharingDepths.clear(); m_Thicknesses.clear(); m_NXStrips.clear(); m_NYStrips.clear(); From a98b7d379a10e8d5f49a7bab3b43c1d4606e152f Mon Sep 17 00:00:00 2001 From: Field Date: Mon, 20 Jul 2026 22:56:57 -0700 Subject: [PATCH 07/15] i somehow made more zombie bump, but committing for benchmarking purposes --- include/MModuleDepthCalibration.h | 10 ++- src/MModuleDepthCalibration.cxx | 109 ++++++++++++++++++++++++++---- 2 files changed, 103 insertions(+), 16 deletions(-) diff --git a/include/MModuleDepthCalibration.h b/include/MModuleDepthCalibration.h index ad88cba3..62b1a3a2 100644 --- a/include/MModuleDepthCalibration.h +++ b/include/MModuleDepthCalibration.h @@ -158,6 +158,9 @@ class MModuleDepthCalibration : public MModule //! Returns the z as a function of ctd, given some ctd and some detector std::tuple CalculateZfromCTD(double CTDvalue, double noise, int DetID,int Grade, bool sean_weighting); + //! Returns the strip with the specified strip ID + MStripHit* GetStrip(std::vector& Strips, int StripID); + //! Returns the strip with most energy from vector Strips, also gives back the energy fraction MStripHit* GetDominantStrip(std::vector& Strips, double& EnergyFraction); @@ -205,7 +208,7 @@ class MModuleDepthCalibration : public MModule // protected members: protected: - + unordered_map> m_ChargeSharingDepths; // maps DetID to the depths for which the charge sharing polynomial correction and charge sharing coefficients per-strip-pair were calculated unordered_map>> m_ChargeSharingCoeffs; // maps StripPairID to a vector of coefficients, for a vector of depths (needs interpolation) vector m_InterpolatedCoeffs; // holder for interpolated charge sharing coeffs between different depths @@ -231,6 +234,7 @@ class MModuleDepthCalibration : public MModule uint64_t m_ErrorSH; uint64_t m_ErrorNullSH; uint64_t m_ErrorNoE; + uint64_t m_ZombieBump; unordered_map m_Detectors; vector m_DetectorIDs; MModuleEnergyCalibration* m_EnergyCalibration; @@ -257,7 +261,9 @@ class MModuleDepthCalibration : public MModule // boolean for use with the card cage at UCSD since it tags all events as detector 11 bool m_UCSDOverride; - + // variable to describe the fraction of charge sharing that we define as a single strip + // TODO -- should this be a variable in the config? + double m_SingleStripChargeSharing = 0.9; // private members: private: diff --git a/src/MModuleDepthCalibration.cxx b/src/MModuleDepthCalibration.cxx index c95ed410..33d5e93f 100644 --- a/src/MModuleDepthCalibration.cxx +++ b/src/MModuleDepthCalibration.cxx @@ -324,14 +324,15 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) } // If the CTD is in range, calculate the depth + else { // FR TODO the last boolean is for sean's weighting method; make it a flag auto [rawZpos, rawZsigma] = CalculateZfromCTD(rawCTD_s, noise,DetID, Grade, true); - + // TODO depth correction loop! // // step 1 -- check the HV side: - int StripPairCode = 10000*DetID + HVStripID; // note, it is the lower strip ID always (eg, 15 if sharing between strip 15 and 16) + //int StripPairCode = 10000*DetID + HVStripID; // note, it is the lower strip ID always (eg, 15 if sharing between strip 15 and 16) // -- how many strips share > 10% of the total energy (or are over slow threshold, maybe?) need this info for next steps // -- check dTAC between adjacent strips with charge sharing and also strips relative to their low-eneryg neighbors // -- correct the HV timing asic jitter bug, if needed, to make everything consistent @@ -339,9 +340,58 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) // -- if charge sharing, calculate the corrected timing value for each strip in teh absense of charge sharing, // -- and also calculate the HV tac as the weighted average, with its own uncertainty // -- note, rawZpos is used in the above calculations! - // + // step 2 -- check the LV side: - StripPairCode = 10000*DetID + 100*LVStripID; // note, it is the lower Strip ID always! + bool ValidatedTiming = false; // should put a flag here eventually TODO; if NN do not have fast timing + bool ZombieBump = false; + vector CorrectedTiming; + vector CorrectedTimingUncertainty; + bool CorrectedChargeSharingLV = false; + // also should put a flag here TODO (is there another flag for charge sharing?) + + if (LVEnergyFraction > m_SingleStripChargeSharing){// if we have one obvious main strip, we are not going to be correcting charge sharing but just checking for zombie bump + // compare with the neighbors, if possible + for (int neighbor = 0; neighbor < 2; neighbor++){// 0 for left neighbor, 1 for right + int pm = 2*neighbor - 1; // -1 for neighbor 0 (neighbor is left); + 1 for neighbor == 1 (right neighbor, which is nominal for the convention StripPairID = left StripID of pair + int NeighborStripID = LVStripID + pm; + MStripHit* NSH = GetStrip(LVStrips, NeighborStripID); + if (NeighborStripID >=0 && NeighborStripID <=63 && NSH){ // neighbor is not a guard ring strip, NSH exists (not a null pointer) TODO and NSH has fast timing! + double dTacData = (LVSH->GetTiming() - NSH->GetTiming())*pm; // always the left strip - right strip; neighbor on left means pm = -1 -> NSH - LVSH timing + double fracData = (NSH->GetEnergy()/(LVSH->GetEnergy() + NSH->GetEnergy())*pm) + 1 - neighbor; // always the fraction on the right stripHit; nominally NSH for right neighbor + int StripPairCode = 10000*DetID + 100*(LVStripID - 1 + neighbor); //LVStripID -1 + 0 = LVStripID -1 (left neighbor); or LVStripID -1 + 1 = LVStripID (LVStrip is the StripID when we consider right negihbor) + int x = (fracData - 0.5); // x is LVEnergyFraction - 0.5 + vector* CSPolyCoeffs = GetChargeSharingPolyCoeffsLV(DetID,rawZpos); // TODO need to actually check and fill the variable that checks the length of this, and check that it's right when loading + vector* CSCoeffs = GetChargeSharingCoeffs(StripPairCode,rawZpos); + if (CSCoeffs && CSPolyCoeffs){ + double dTacExpect = (CSPolyCoeffs->at(0)*x + CSPolyCoeffs->at(1)*x*x*x)*CSCoeffs->at(0) + CSCoeffs->at(1);// TODO update if not cubic polynomial + // TODO we should display (dTacData - dTacExpect)*pm to keep an eye on the prevalence of the bump + if (rawZpos > -5 && dTacData < 500 && dTacData > -500 && (dTacData - dTacExpect)*pm > 2*noise){// zombie bump! TODO update bump criteria, fix noise, remove Slow Timing check once checked earlier + ZombieBump = true; + ValidatedTiming = false; + if ((dTacData - dTacExpect)*pm > -1*pm*dTacExpect) { // zombie bump with good neighbor. Can we not always do this in this case, though, since we know what the timing should be? + CorrectedTiming.push_back(LVTiming-(dTacData - dTacExpect)*pm); + CorrectedTimingUncertainty.push_back(noise*2); // to do: quantify and make into something real + } + } else { + if (dTacData < 500 && dTacData > -500 && !ZombieBump && CSCoeffs && CSPolyCoeffs) ValidatedTiming = true; + } + } + } + } + } else { // charge sharing correction + CorrectedChargeSharingLV = true; + } + if (ZombieBump) m_ZombieBump++; + if (CorrectedTiming.size() > 0){// if we have a correction. Also, implemented weighting! TODO an + double correctionSum = 0; + for (unsigned int j = 0; j < CorrectedTiming.size(); j++) { + // TODO check that they are consistent and drop one if not.... + correctionSum += CorrectedTiming.at(j); + } + LVTiming = correctionSum / CorrectedTiming.size(); + } + rawCTD = (HVTiming - LVTiming); + rawCTD_s = (rawCTD - Coeffs->at(1))/(Coeffs->at(0)); // -- how many strips share > 10% of the total energy (or are over the slow threshold, maybe?) need this info for the next steps // -- check dTAC between adjacent strips with charge sharing and also strips relative to their low-energy neighbors // -- deal with the zombie bump! @@ -349,21 +399,32 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) // -- if charge sharing, calculate the corrected timing value for each strip in the absence of charge sharing, // -- and then also calculate the corrected LV tac as the weighted average, with its own uncertainty (they should be consistent) // -- note, rawZpos is used in teh above calculations - // + // step 3 -- calculate new corrected CTD, and from that calculate new corrected Z // // bonus points: implement x and y localization based on info in step 1 and 2 with charge sharing :) // bonus points -- can re-calculated depth between 1 and 2 if significant changes with HV correction would impact ZB correction! - Zpos = rawZpos; - Zsigma = rawZsigma; - - // Add the depth to the GUI histogram. - if (Event->HasStripPairingError()==false) { - if (HasExpos() == true) { - m_ExpoDepthCalibration->AddDepth(DetID, Zpos); - } + + + if ((rawCTD_s < (Xmin - 2.0*noise)) || (rawCTD_s > (Xmax + 2.0*noise))) { + H->SetNoDepth(); + Event->SetDepthCalibrationError("Out of Range"); + ++m_Error2; } - m_NoError+=1; + // If the CTD is in range, calculate the depth + else { + auto [rawZpos, rawZsigma] = CalculateZfromCTD(rawCTD_s, noise,DetID, Grade, true); + Zpos = rawZpos; + Zsigma = rawZsigma; + + // Add the depth to the GUI histogram. + if (Event->HasStripPairingError()==false) { + if (HasExpos() == true) { + m_ExpoDepthCalibration->AddDepth(DetID, Zpos); + } + } + m_NoError+=1; + } } } @@ -500,6 +561,25 @@ MStripHit* MModuleDepthCalibration::GetDominantStrip(vector& Strips, } + + +///////////////////////////////////////////////////////////////////////////////// + + +MStripHit* MModuleDepthCalibration::GetStrip(vector& Strips, int StripID) +{ + MStripHit* MaxStrip = nullptr; + + // Iterate through strip hits and get the strip with highest energy + for (const auto SH : Strips) { + if (SH->GetStripID() == StripID) return SH; + } + return MaxStrip; +} + + + + ///////////////////////////////////////////////////////////////////////////////// @@ -1378,6 +1458,7 @@ void MModuleDepthCalibration::Finalize() cout << "Number of hits with no strip hits on one or both sides: " << m_ErrorSH << endl; cout << "Number of hits with null strip hits: " << m_ErrorNullSH << endl; cout << "Number of hits 0 energy on a strip hit: " << m_ErrorNoE << endl; + cout << "Number of hits with zombie bump:" << m_ZombieBump << endl; // Clean up maps and vectors m_Coeffs.clear(); From 046276994bd82a3df90c4bf266e229bfc6bfd9c6 Mon Sep 17 00:00:00 2001 From: Field Date: Thu, 23 Jul 2026 20:10:24 -0700 Subject: [PATCH 08/15] added a very basic bump correction routine; added a zombie bump Expo to the displays; reverted to non-weighted depth reco for now because of problem with putting everyting at z = 0; added before and after correction plots to the Depth Expo --- include/MGUIExpoDepthCalibration.h | 32 ++- include/MGUIExpoPlotTacDiff.h | 88 +++++++ include/MModuleDepthCalibration.h | 12 +- src/MGUIExpoDepthCalibration.cxx | 261 ++++++++++++++++++-- src/MGUIExpoPlotTacDiff.cxx | 368 +++++++++++++++++++++++++++++ src/MModuleDepthCalibration.cxx | 222 +++++++++++------ 6 files changed, 874 insertions(+), 109 deletions(-) create mode 100644 include/MGUIExpoPlotTacDiff.h create mode 100644 src/MGUIExpoPlotTacDiff.cxx diff --git a/include/MGUIExpoDepthCalibration.h b/include/MGUIExpoDepthCalibration.h index 5429ca67..a393def7 100644 --- a/include/MGUIExpoDepthCalibration.h +++ b/include/MGUIExpoDepthCalibration.h @@ -28,6 +28,9 @@ #include #include #include +#include +#include +#include // MEGAlib libs: #include "MGlobal.h" @@ -83,21 +86,44 @@ class MGUIExpoDepthCalibration : public MGUIExpo //! 0 1 2 3 //! 4 5 6 7 //! 8 9 10 11 - void AddDepth(unsigned int DetID, double Depth); + void AddDepth(unsigned int DetID, int LVStrip, int HVStrip, double Depth); + void AddRawDepth(unsigned int DetID, int LVStrip, int HVStrip, double Depth); + void OnStripSelectionChanged(); + void RebuildDisplayHistograms(); + virtual bool ProcessMessage(long Message, long Parameter1, long Parameter2); + // protected methods: protected: + + // protected members: + protected: + std::map m_RawDepthPerStrip; + std::map m_DepthPerStrip; + TGComboBox* m_SideSelector; + TGNumberEntry* m_StripMinEntry; + TGNumberEntry* m_StripMaxEntry; - // protected members: - protected: + int m_SelectedSide; // 0=LV, 1=HV + int m_StripMin; + int m_StripMax; + + int GetStripKey(int DetID, int Side, int StripID) { return DetID * 1000 + Side * 100 + StripID; } // private members: private: + TLegend* m_Legend = nullptr; + TGTextButton* m_UpdateSelectionButton; + + static const int c_UpdateSelection = 1001; + + //! Depth canvas unordered_map m_DepthCanvases; //! Depth vs detector ID histogram unordered_map m_DepthHistograms; + unordered_map m_RawDepthHistograms; //! Detectors in x direction unsigned int m_NColumns; diff --git a/include/MGUIExpoPlotTacDiff.h b/include/MGUIExpoPlotTacDiff.h new file mode 100644 index 00000000..bb6fa8ca --- /dev/null +++ b/include/MGUIExpoPlotTacDiff.h @@ -0,0 +1,88 @@ +#ifndef __MGUIExpoPlotTacDiff__ +#define __MGUIExpoPlotTacDiff__ + +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include "MGUIExpo.h" +#include "MModule.h" + +class MGUIExpoPlotTacDiff : public MGUIExpo +{ +public: + MGUIExpoPlotTacDiff(MModule* Module); + virtual ~MGUIExpoPlotTacDiff(); + + virtual void Reset(); + virtual void Create(); + virtual void Update(); + virtual void Export(const MString& FileName); + + void SetHistogramParameters(unsigned int DetID, unsigned int NBinsDepth, double DepthMin, double DepthMax, + unsigned int NBinsDtac, double DtacMin, double DtacMax, + unsigned int NBinsFrac, double FracMin, double FracMax); + void AddData(int StripPairCode, double Depth, double dTac, double Fraction, double dTacAlt,double Energy); + + //void OnDetectorSelected(Int_t DetID); + //void OnSideSelected(Int_t Side); + //void OnStripSelected(); + //void OnEnergyRangeChanged(); + virtual bool ProcessMessage(long Message, long Parameter1, long Parameter2); + void OnUpdateSelection(); + + +protected: + void RedrawPlots(); + + int m_SelectedDetector; + int m_SelectedSide; + int m_SelectedStrip; + int m_SelectedStripPairCode; + + TGComboBox* m_DetectorSelector; + TGComboBox* m_SideSelector; + TGNumberEntry* m_StripEntry; + TGNumberEntry* m_EnergyMinEntry; + TGNumberEntry* m_EnergyMaxEntry; + double m_EnergyMin; + double m_EnergyMax; + + TRootEmbeddedCanvas* m_DtacVsDepthCanvas; + TRootEmbeddedCanvas* m_DtacVsFracCanvas; + TRootEmbeddedCanvas* m_DtacVsDtacCanvas; + + // Key: StripPairCode = 10000*DetID + 100*LVStripID (for LV) or 10000*DetID + HVStripID (for HV) + std::map m_DtacVsDepthHistograms; + std::map m_DtacVsFracHistograms; + std::map m_DtacVsDtacHistograms; + + std::map m_NBinsDepth; + std::map m_DepthMin; + std::map m_DepthMax; + std::map m_NBinsDtac; + std::map m_DtacMin; + std::map m_DtacMax; + std::map m_NBinsFrac; + std::map m_FracMin; + std::map m_FracMax; + + std::vector m_DetIDs; +private: + static const int c_UpdateSelection = 1001; + TGTextButton* m_UpdateSelectionButton; + +#ifdef ___CLING___ +public: + ClassDef(MGUIExpoPlotTacDiff, 0) +#endif +}; + +#endif diff --git a/include/MModuleDepthCalibration.h b/include/MModuleDepthCalibration.h index 62b1a3a2..526dddec 100644 --- a/include/MModuleDepthCalibration.h +++ b/include/MModuleDepthCalibration.h @@ -32,6 +32,7 @@ #include "MDStrip3D.h" #include "MDShapeBRIK.h" #include "MGUIExpoDepthCalibration.h" +#include "MGUIExpoPlotTacDiff.h" // Forward declarations: @@ -182,11 +183,11 @@ class MModuleDepthCalibration : public MModule int GetHitGrade(MHit* H); //! Return the coefficients of the dTac polynomial for a detector at a depth - vector* GetChargeSharingPolyCoeffsHV(int DetID,double z); - vector* GetChargeSharingPolyCoeffsLV(int DetID,double z); + vector GetChargeSharingPolyCoeffsHV(int DetID,double z); + vector GetChargeSharingPolyCoeffsLV(int DetID,double z); //! Return the charge-sharing stretch/offset coefficients for a strip-pair at a depth - vector* GetChargeSharingCoeffs(int StripPairCode,double z); + vector GetChargeSharingCoeffs(int StripPairCode,double z); //! Return the ctd-z stretch/offset coefficients for a pixel vector* GetPixelCoeffs(int PixelCode); @@ -211,7 +212,6 @@ class MModuleDepthCalibration : public MModule unordered_map> m_ChargeSharingDepths; // maps DetID to the depths for which the charge sharing polynomial correction and charge sharing coefficients per-strip-pair were calculated unordered_map>> m_ChargeSharingCoeffs; // maps StripPairID to a vector of coefficients, for a vector of depths (needs interpolation) - vector m_InterpolatedCoeffs; // holder for interpolated charge sharing coeffs between different depths unordered_map>> m_ChargeSharingPolyCoeffsLV; // maps DetID to the LV coefficients for dTAC vs CS map and charge sharing correction vs CS map, for a vector of depths unordered_map>> m_ChargeSharingPolyCoeffsHV; // maps DetID to the HV coefficients for dTAC vs CS map and charge sharing correction vs CS map, for a given depth unordered_map> m_Coeffs; // maps pix id to a vector of coefficients... @@ -235,10 +235,14 @@ class MModuleDepthCalibration : public MModule uint64_t m_ErrorNullSH; uint64_t m_ErrorNoE; uint64_t m_ZombieBump; + uint64_t m_ChargeSharingLV; + uint64_t m_ChargeSharingHV; unordered_map m_Detectors; vector m_DetectorIDs; MModuleEnergyCalibration* m_EnergyCalibration; MGUIExpoDepthCalibration* m_ExpoDepthCalibration; + MGUIExpoPlotTacDiff* m_ExpoPlotTacDiff; + // The CTD Map maps each detector (int) to a 2D array of CTD values. unordered_map>> m_CTDMap; diff --git a/src/MGUIExpoDepthCalibration.cxx b/src/MGUIExpoDepthCalibration.cxx index c4a0751e..cfac4742 100644 --- a/src/MGUIExpoDepthCalibration.cxx +++ b/src/MGUIExpoDepthCalibration.cxx @@ -49,8 +49,13 @@ MGUIExpoDepthCalibration::MGUIExpoDepthCalibration(MModule* Module) : MGUIExpo(M // standard constructor // Set the new title of the tab here: - m_TabTitle = "Depth Calibration"; - + m_TabTitle = "Depth"; + m_SideSelector = nullptr; + m_StripMinEntry = nullptr; + m_StripMaxEntry = nullptr; + m_SelectedSide = 0; + m_StripMin = 0; + m_StripMax = 63; // Set the histogram arrangment // SetDepthHistogramArrangement(1, 1); @@ -76,9 +81,10 @@ void MGUIExpoDepthCalibration::Reset() //! Reset the data in the UI m_Mutex.Lock(); - for (auto H: m_DepthHistograms) { - (H.second)->Reset(); - } + for (auto H: m_DepthHistograms) (H.second)->Reset(); + for (auto H: m_RawDepthHistograms) (H.second)->Reset(); + for (auto& H : m_RawDepthPerStrip) H.second->Reset(); + for (auto& H : m_DepthPerStrip) H.second->Reset(); m_Mutex.UnLock(); } @@ -115,10 +121,25 @@ void MGUIExpoDepthCalibration::SetDepthHistogramArrangement(vector TH1D* Depth = new TH1D("", "Depth", m_NBins[DetID], m_Min[DetID], m_Max[DetID]); Depth->SetXTitle("Depth [cm]"); Depth->SetYTitle("counts"); - Depth->SetFillColor(kAzure+7); - + Depth->SetFillColorAlpha(kAzure+7,0.5); + Depth->SetFillStyle(3001); m_DepthHistograms[DetID] = Depth; - // m_DepthCanvases[DetID] = 0; + + TH1D* RawDepth = new TH1D("", "Depth", m_NBins[DetID], m_Min[DetID], m_Max[DetID]); + RawDepth->SetXTitle("Depth [cm]"); + RawDepth->SetYTitle("counts"); + RawDepth->SetFillColorAlpha(kRed,0.5); + Depth->SetFillStyle(3001); + m_RawDepthHistograms[DetID] = RawDepth;// m_DepthCanvases[DetID] = 0; + + // Create per-strip histograms + for (int side = 0; side < 2; side++) { + for (int strip = 0; strip < 64; strip++) { + int key = GetStripKey(DetID, side, strip); + m_RawDepthPerStrip[key] = new TH1D("", "", m_NBins[DetID], m_Min[DetID], m_Max[DetID]); + m_DepthPerStrip[key] = new TH1D("", "", m_NBins[DetID], m_Min[DetID], m_Max[DetID]); + } + } ++column; } @@ -141,20 +162,25 @@ void MGUIExpoDepthCalibration::SetDepthHistogramArrangement(vector void MGUIExpoDepthCalibration::SetDepthHistogramParameters(unsigned int DetID, unsigned int NBins, double DepthMin, double DepthMax) { - // Set the energy histogram parameters - m_Mutex.Lock(); m_NBins[DetID] = NBins; m_Min[DetID] = DepthMin; m_Max[DetID] = DepthMax; - TH1D* H = m_DepthHistograms[DetID]; - H->SetBins(NBins, DepthMin, DepthMax); + + // Only update bins if histograms already exist + if (m_DepthHistograms.count(DetID) > 0 && m_DepthHistograms[DetID] != nullptr) { + m_DepthHistograms[DetID]->SetBins(NBins, DepthMin, DepthMax); + } + if (m_RawDepthHistograms.count(DetID) > 0 && m_RawDepthHistograms[DetID] != nullptr) { + m_RawDepthHistograms[DetID]->SetBins(NBins, DepthMin, DepthMax); + } m_Mutex.UnLock(); } + //////////////////////////////////////////////////////////////////////////////// @@ -172,22 +198,95 @@ void MGUIExpoDepthCalibration::SetDepthHistogramName(unsigned int DetID, MString } +//////////////////////////////////////////////////////////////////////////////// + +void MGUIExpoDepthCalibration::AddRawDepth(unsigned int DetID, int LVStrip, int HVStrip, double Depth) +{ + m_Mutex.Lock(); + int lvKey = GetStripKey(DetID, 0, LVStrip); + int hvKey = GetStripKey(DetID, 1, HVStrip); + if (m_RawDepthPerStrip.count(lvKey) > 0) m_RawDepthPerStrip[lvKey]->Fill(Depth); + if (m_RawDepthPerStrip.count(hvKey) > 0) m_RawDepthPerStrip[hvKey]->Fill(Depth); + m_Mutex.UnLock(); +} + //////////////////////////////////////////////////////////////////////////////// -void MGUIExpoDepthCalibration::AddDepth(unsigned int DetID, double Depth) +void MGUIExpoDepthCalibration::AddDepth(unsigned int DetID, int LVStrip, int HVStrip, double Depth) { - // Add data to the energy histogram + m_Mutex.Lock(); + int lvKey = GetStripKey(DetID, 0, LVStrip); + int hvKey = GetStripKey(DetID, 1, HVStrip); + if (m_DepthPerStrip.count(lvKey) > 0) m_DepthPerStrip[lvKey]->Fill(Depth); + if (m_DepthPerStrip.count(hvKey) > 0) m_DepthPerStrip[hvKey]->Fill(Depth); + m_Mutex.UnLock(); +} + +//////////////////////////////////////////////////////////////////////////////// +void MGUIExpoDepthCalibration::OnStripSelectionChanged() +{ m_Mutex.Lock(); - if (m_DepthHistograms.find(DetID) != m_DepthHistograms.end()) { - m_DepthHistograms[DetID]->Fill(Depth); + m_SelectedSide = m_SideSelector->GetSelected(); + m_StripMin = (int)m_StripMinEntry->GetNumber(); + m_StripMax = (int)m_StripMaxEntry->GetNumber(); + if (m_StripMin < 0) m_StripMin = 0; + if (m_StripMax > 63) m_StripMax = 63; + if (m_StripMin > m_StripMax) m_StripMin = m_StripMax; + + RebuildDisplayHistograms(); + + // Rescale Y-axis + double Max = 0; + for (const auto& pair : m_DepthHistograms) { + TH1D* H = pair.second; + for (int bx = 2; bx < H->GetNbinsX(); ++bx) { + if (Max < H->GetBinContent(bx)) Max = H->GetBinContent(bx); + } + } + for (const auto& pair : m_RawDepthHistograms) { + TH1D* H = pair.second; + for (int bx = 2; bx < H->GetNbinsX(); ++bx) { + if (Max < H->GetBinContent(bx)) Max = H->GetBinContent(bx); + } + } + Max *= 1.1; + if (Max == 0) Max = 1.0; + + for (const auto& pair : m_DepthHistograms) pair.second->SetMaximum(Max); + for (const auto& pair : m_RawDepthHistograms) pair.second->SetMaximum(Max); + + // Redraw canvases + for (auto& C : m_DepthCanvases) { + C.second->GetCanvas()->Modified(); + C.second->GetCanvas()->Update(); } m_Mutex.UnLock(); } +//////////////////////////////////////////////////////////////////////////////// + +void MGUIExpoDepthCalibration::RebuildDisplayHistograms() +{ + cout << "RebuildDisplayHistograms called" << endl; + for (auto& pair : m_DepthHistograms) pair.second->Reset(); + for (auto& pair : m_RawDepthHistograms) pair.second->Reset(); + + for (auto& pair : m_DepthHistograms) { + unsigned int DetID = pair.first; + for (int s = m_StripMin; s <= m_StripMax; s++) { + int key = GetStripKey(DetID, m_SelectedSide, s); + if (m_DepthPerStrip.count(key) > 0) + m_DepthHistograms[DetID]->Add(m_DepthPerStrip[key]); + cout << " added "<< key << ": " << m_DepthPerStrip[key]->GetEntries() << " entries" << endl; + if (m_RawDepthPerStrip.count(key) > 0) + m_RawDepthHistograms[DetID]->Add(m_RawDepthPerStrip[key]); + } + } +} //////////////////////////////////////////////////////////////////////////////// @@ -200,13 +299,64 @@ void MGUIExpoDepthCalibration::Create() if (m_IsCreated == true) return; m_Mutex.Lock(); - - TGLayoutHints* CanvasLayout = new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX | kLHintsExpandY, 2, 2, 2, 2); - for (unsigned int y = 0; y < m_DetectorMap.size(); ++y) { - TGHorizontalFrame* HFrame = new TGHorizontalFrame(this); - AddFrame(HFrame, CanvasLayout); + TGLayoutHints* ExpandLayout = new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX | kLHintsExpandY, 2, 2, 2, 2); + TGLayoutHints* RowLayout = new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX, 2, 2, 2, 2); + TGLayoutHints* LabelLayout = new TGLayoutHints(kLHintsLeft | kLHintsCenterY, 2, 5, 0, 0); + TGLayoutHints* EntryLayout = new TGLayoutHints(kLHintsLeft | kLHintsCenterY, 0, 2, 0, 0); + + TGHorizontalFrame* MainHFrame = new TGHorizontalFrame(this); + AddFrame(MainHFrame, ExpandLayout); + + // === Left: Controls === + TGVerticalFrame* ControlFrame = new TGVerticalFrame(MainHFrame, 180, 400); + MainHFrame->AddFrame(ControlFrame, new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandY, 5, 5, 5, 5)); + + // Strip selection: Side dropdown + min + max + TGHorizontalFrame* StripRow = new TGHorizontalFrame(ControlFrame); + ControlFrame->AddFrame(StripRow, RowLayout); + + m_SideSelector = new TGComboBox(StripRow); + m_SideSelector->AddEntry("LV", 0); + m_SideSelector->AddEntry("HV", 1); + m_SideSelector->Select(0); + m_SideSelector->Resize(45, 20); + //m_SideSelector->Connect("Selected(Int_t)", "MGUIExpoDepthCalibration", this, "OnStripSelectionChanged()"); + StripRow->AddFrame(m_SideSelector, EntryLayout); + + m_StripMinEntry = new TGNumberEntry(StripRow, 0, 3, -1, + TGNumberFormat::kNESInteger, + TGNumberFormat::kNEANonNegative, + TGNumberFormat::kNELLimitMinMax, 0, 63); + //m_StripMinEntry->Connect("ValueSet(Long_t)", "MGUIExpoDepthCalibration", this, "OnStripSelectionChanged()"); + //m_StripMinEntry->GetNumberEntry()->Connect("ReturnPressed()", "MGUIExpoDepthCalibration", this, "OnStripSelectionChanged()"); + m_StripMinEntry->Resize(45, 20); + StripRow->AddFrame(m_StripMinEntry, EntryLayout); + + m_StripMaxEntry = new TGNumberEntry(StripRow, 63, 3, -1, + TGNumberFormat::kNESInteger, + TGNumberFormat::kNEANonNegative, + TGNumberFormat::kNELLimitMinMax, 0, 63); + //m_StripMaxEntry->Connect("ValueSet(Long_t)", "MGUIExpoDepthCalibration", this, "OnStripSelectionChanged()"); + //m_StripMaxEntry->GetNumberEntry()->Connect("ReturnPressed()", "MGUIExpoDepthCalibration", this, "OnStripSelectionChanged()"); + m_StripMaxEntry->Resize(45, 20); + StripRow->AddFrame(m_StripMaxEntry, EntryLayout); + + // Update + TGTextButton* m_UpdateSelectionButton = new TGTextButton(ControlFrame, "Update Plot Selection", c_UpdateSelection); + m_UpdateSelectionButton->Associate(this); + ControlFrame->AddFrame(m_UpdateSelectionButton, RowLayout); + + // === Right: Plots === + TGVerticalFrame* PlotFrame = new TGVerticalFrame(MainHFrame); + MainHFrame->AddFrame(PlotFrame, ExpandLayout); + + TGLayoutHints* CanvasLayout = new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX | kLHintsExpandY, 2, 2, 2, 2); + for (unsigned int y = 0; y < m_DetectorMap.size(); ++y) { + TGHorizontalFrame* HFrame = new TGHorizontalFrame(PlotFrame); + PlotFrame->AddFrame(HFrame, CanvasLayout); + for (unsigned int x = 0; x < m_DetectorMap[y].size(); ++x) { unsigned int DetID = m_DetectorMap[y][x]; TRootEmbeddedCanvas* DepthCanvas = new TRootEmbeddedCanvas("Depth", HFrame, 100, 100); @@ -214,13 +364,30 @@ void MGUIExpoDepthCalibration::Create() m_DepthCanvases[DetID] = DepthCanvas; DepthCanvas->GetCanvas()->cd(); - m_DepthHistograms[DetID]->Draw("colz"); + m_RawDepthHistograms[DetID]->Draw("colz"); + m_DepthHistograms[DetID]->Draw("SAME"); + + // Add legend + if (x == 0 && m_DepthHistograms[0] != nullptr) { + + m_Legend = new TLegend(0.15, 0.75, 0.35, 0.88); + m_Legend->SetTextSize(0.03); + m_Legend->AddEntry(m_RawDepthHistograms[DetID], "Raw Depth", "f"); + m_Legend->AddEntry(m_DepthHistograms[DetID],"Final Depth", "f"); + m_Legend->Draw(); + } + DepthCanvas->GetCanvas()->Update(); } } + m_SelectedSide = m_SideSelector->GetSelected(); + m_StripMin = (int)m_StripMinEntry->GetNumber(); + m_StripMax = (int)m_StripMaxEntry->GetNumber(); m_IsCreated = true; + + m_Mutex.UnLock(); } @@ -234,6 +401,13 @@ void MGUIExpoDepthCalibration::Update() m_Mutex.Lock(); + if (m_SideSelector) m_SideSelector->SetEnabled(false); + if (m_StripMinEntry) m_StripMinEntry->SetState(false); + if (m_StripMaxEntry) m_StripMaxEntry->SetState(false); + + + RebuildDisplayHistograms(); + double Max = 0; // for (auto H : m_DepthHistograms) { for ( const auto dethistpair : m_DepthHistograms ){ @@ -244,11 +418,17 @@ void MGUIExpoDepthCalibration::Update() } } } - Max *= 1.1; - for ( const auto dethistpair : m_DepthHistograms ){ + for ( const auto dethistpair : m_RawDepthHistograms ){ TH1D* H = dethistpair.second; - H->SetMaximum(Max); + for (int bx = 2; bx < H->GetNbinsX(); ++bx) { // Skip first and last + if (Max < H->GetBinContent(bx)) { + Max = H->GetBinContent(bx); + } + } } + Max *= 1.1; + for (const auto& pair : m_DepthHistograms) pair.second->SetMaximum(Max); + for (const auto& pair : m_RawDepthHistograms) pair.second->SetMaximum(Max); for (auto C : m_DepthCanvases) { @@ -256,10 +436,38 @@ void MGUIExpoDepthCalibration::Update() (C.second)->GetCanvas()->Update(); } + if (m_SideSelector) m_SideSelector->SetEnabled(true); + if (m_StripMinEntry) m_StripMinEntry->SetState(true); + if (m_StripMaxEntry) m_StripMaxEntry->SetState(true); + m_Mutex.UnLock(); } +//////////////////////////////////////////////////////////////////////////////// +bool MGUIExpoDepthCalibration::ProcessMessage(long Message, long Parameter1, long Parameter2) +{ + switch (GET_MSG(Message)) { + case kC_COMMAND: + switch (GET_SUBMSG(Message)) { + case kCM_BUTTON: + switch (Parameter1) { + case c_UpdateSelection: + OnStripSelectionChanged(); + return true; + default: + break; + } + break; + default: + break; + } + break; + default: + break; + } + return true; +} //////////////////////////////////////////////////////////////////////////////// @@ -275,7 +483,8 @@ void MGUIExpoDepthCalibration::Export(const MString& FileName) for (unsigned int x = 0; x < m_DetectorMap[y].size(); ++x) { unsigned int DetID = m_DetectorMap[y][x]; P->cd((x+1) + m_NColumns*y); - m_DepthHistograms[DetID]->DrawCopy("colz"); + m_RawDepthHistograms[DetID]->DrawCopy("colz"); + m_DepthHistograms[DetID]->DrawCopy("SAME"); } } P->SaveAs(FileName); diff --git a/src/MGUIExpoPlotTacDiff.cxx b/src/MGUIExpoPlotTacDiff.cxx new file mode 100644 index 00000000..300ca631 --- /dev/null +++ b/src/MGUIExpoPlotTacDiff.cxx @@ -0,0 +1,368 @@ +#include "MGUIExpoPlotTacDiff.h" + +#include + +#include +#include +#include +#include +#include + +#include "MStreams.h" + +using namespace std; + +#ifdef ___CLING___ +ClassImp(MGUIExpoPlotTacDiff) +#endif + +//////////////////////////////////////////////////////////////////////////////// + +MGUIExpoPlotTacDiff::MGUIExpoPlotTacDiff(MModule* Module) : MGUIExpo(Module) +{ + m_TabTitle = "Zombie"; + m_SelectedDetector = -1; + m_SelectedSide = 0; + m_SelectedStrip = 0; + m_SelectedStripPairCode = 0; + m_DetectorSelector = nullptr; + m_SideSelector = nullptr; + m_StripEntry = nullptr; + m_DtacVsDepthCanvas = nullptr; + m_DtacVsFracCanvas = nullptr; + m_DtacVsDtacCanvas = nullptr; + m_EnergyMinEntry = nullptr; + m_EnergyMaxEntry = nullptr; + m_EnergyMin = 0; + m_EnergyMax = 1e9; // default: accept all + SetCleanup(kDeepCleanup); +} + +//////////////////////////////////////////////////////////////////////////////// + +MGUIExpoPlotTacDiff::~MGUIExpoPlotTacDiff() +{ +} + +//////////////////////////////////////////////////////////////////////////////// + +void MGUIExpoPlotTacDiff::Reset() +{ + m_Mutex.Lock(); + for (auto& pair : m_DtacVsDepthHistograms) pair.second->Reset(); + for (auto& pair : m_DtacVsFracHistograms) pair.second->Reset(); + for (auto& pair : m_DtacVsDtacHistograms) pair.second->Reset(); + m_Mutex.UnLock(); +} + +//////////////////////////////////////////////////////////////////////////////// + +void MGUIExpoPlotTacDiff::SetHistogramParameters(unsigned int DetID, unsigned int NBinsDepth, double DepthMin, double DepthMax, + unsigned int NBinsDtac, double DtacMin, double DtacMax, + unsigned int NBinsFrac, double FracMin, double FracMax) +{ + m_Mutex.Lock(); + + m_NBinsDepth[DetID] = NBinsDepth; + m_DepthMin[DetID] = DepthMin; + m_DepthMax[DetID] = DepthMax; + m_NBinsDtac[DetID] = NBinsDtac; + m_DtacMin[DetID] = DtacMin; + m_DtacMax[DetID] = DtacMax; + m_NBinsFrac[DetID] = NBinsFrac; + m_FracMin[DetID] = FracMin; + m_FracMax[DetID] = FracMax; + + bool found = false; + for (auto id : m_DetIDs) { if (id == DetID) { found = true; break; } } + if (!found) m_DetIDs.push_back(DetID); + + // LV strips: code = 10000*DetID + 100*StripID + 99 + for (int s = 0; s < 63; s++) { + int key = 10000 * DetID + 100 * s + 99; + if (m_DtacVsDepthHistograms.count(key) == 0) { + m_DtacVsDepthHistograms[key] = new TH2D("", TString::Format("Det %d LV %d: dTAC vs Depth", DetID, s), + NBinsDepth, DepthMin, DepthMax, NBinsDtac, DtacMin, DtacMax); + m_DtacVsDepthHistograms[key]->SetXTitle("Depth [cm]"); + m_DtacVsDepthHistograms[key]->SetYTitle("dTAC [ns]"); + + m_DtacVsFracHistograms[key] = new TH2D("", TString::Format("Det %d LV %d: dTAC vs Frac", DetID, s), + NBinsFrac, FracMin, FracMax, NBinsDtac, DtacMin, DtacMax); + m_DtacVsFracHistograms[key]->SetXTitle("Fraction"); + m_DtacVsFracHistograms[key]->SetYTitle("dTAC [ns]"); + + m_DtacVsDtacHistograms[key] = new TH2D("", TString::Format("Det %d LV %d: dTAC vs dTAC", DetID, s), + NBinsDtac, DtacMin, DtacMax, NBinsDtac, DtacMin, DtacMax); + m_DtacVsDtacHistograms[key]->SetXTitle("dTAC [ns]"); + m_DtacVsDtacHistograms[key]->SetYTitle("dTAC alt [ns]"); + } + } + + // HV strips: code = 10000*DetID + StripID + for (int s = 0; s < 63; s++) { + int key = 10000 * DetID + s + 9900; + if (m_DtacVsDepthHistograms.count(key) == 0) { + m_DtacVsDepthHistograms[key] = new TH2D("", TString::Format("Det %d HV %d: dTAC vs Depth", DetID, s), + NBinsDepth, DepthMin, DepthMax, NBinsDtac, DtacMin, DtacMax); + m_DtacVsDepthHistograms[key]->SetXTitle("Depth [cm]"); + m_DtacVsDepthHistograms[key]->SetYTitle("dTAC [ns]"); + + m_DtacVsFracHistograms[key] = new TH2D("", TString::Format("Det %d HV %d: dTAC vs Frac", DetID, s), + NBinsFrac, FracMin, FracMax, NBinsDtac, DtacMin, DtacMax); + m_DtacVsFracHistograms[key]->SetXTitle("Fraction"); + m_DtacVsFracHistograms[key]->SetYTitle("dTAC [ns]"); + + m_DtacVsDtacHistograms[key] = new TH2D("", TString::Format("Det %d HV %d: dTAC vs dTAC", DetID, s), + NBinsDtac, DtacMin, DtacMax, NBinsDtac, DtacMin, DtacMax); + m_DtacVsDtacHistograms[key]->SetXTitle("dTAC [ns]"); + m_DtacVsDtacHistograms[key]->SetYTitle("dTAC alt [ns]"); + } + } + + m_Mutex.UnLock(); + cout << "MGUIExpoPlotTacDiff: created histograms, total count = " << m_DtacVsDepthHistograms.size() << endl; +} + +//////////////////////////////////////////////////////////////////////////////// + +void MGUIExpoPlotTacDiff::AddData(int StripPairCode, double Depth, double dTac, double Fraction, double dTacAlt,double Energy) +{ + m_Mutex.Lock(); + if (m_DtacVsDepthHistograms.count(StripPairCode) > 0) { + if (!std::isnan(Depth) && !std::isnan(dTac)) + m_DtacVsDepthHistograms[StripPairCode]->Fill(Depth, dTac); + if (!std::isnan(Fraction) && !std::isnan(dTac)) + m_DtacVsFracHistograms[StripPairCode]->Fill(Fraction, dTac); + if (!std::isnan(dTac) && !std::isnan(dTacAlt)) + m_DtacVsDtacHistograms[StripPairCode]->Fill(dTac, dTacAlt); + }else { + cout << "MGUIExpoPlotTacDiff::AddData: StripPairCode " << StripPairCode << " not found in histograms!" << endl; + } + m_Mutex.UnLock(); +} + +//////////////////////////////////////////////////////////////////////////////// + +void MGUIExpoPlotTacDiff::Create() +{ + if (m_IsCreated == true) return; + if (m_DtacVsDepthHistograms.empty()) return; + + m_Mutex.Lock(); + + m_SelectedDetector = m_DetIDs[0]; + m_SelectedSide = 0; + m_SelectedStrip = 0; + m_SelectedStripPairCode = 10000 * m_SelectedDetector + 100 * m_SelectedStrip + 99; + + TGLayoutHints* ExpandLayout = new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX | kLHintsExpandY, 2, 2, 2, 2); + TGLayoutHints* WidgetLayout = new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX, 5, 5, 5, 5); + + // Main vertical frame (two rows) + TGVerticalFrame* MainVFrame = new TGVerticalFrame(this); + AddFrame(MainVFrame, ExpandLayout); + + // === Top row: two canvases side by side === + TGHorizontalFrame* TopRow = new TGHorizontalFrame(MainVFrame); + MainVFrame->AddFrame(TopRow, ExpandLayout); + + m_DtacVsDepthCanvas = new TRootEmbeddedCanvas("DtacVsDepth", TopRow, 400, 300); + TopRow->AddFrame(m_DtacVsDepthCanvas, ExpandLayout); + + m_DtacVsFracCanvas = new TRootEmbeddedCanvas("DtacVsFrac", TopRow, 400, 300); + TopRow->AddFrame(m_DtacVsFracCanvas, ExpandLayout); + + // === Bottom row: controls left, square plot right === + TGHorizontalFrame* BottomRow = new TGHorizontalFrame(MainVFrame); + MainVFrame->AddFrame(BottomRow, ExpandLayout); + +// === Controls === + TGVerticalFrame* ControlFrame = new TGVerticalFrame(BottomRow, 180, 300); + BottomRow->AddFrame(ControlFrame, new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandY, 5, 5, 5, 5)); + + TGLayoutHints* RowLayout = new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX, 2, 2, 2, 2); + TGLayoutHints* LabelLayout = new TGLayoutHints(kLHintsLeft | kLHintsCenterY, 2, 5, 0, 0); + TGLayoutHints* EntryLayout = new TGLayoutHints(kLHintsLeft | kLHintsCenterY, 0, 2, 0, 0); + + // Detector row + TGHorizontalFrame* DetRow = new TGHorizontalFrame(ControlFrame); + ControlFrame->AddFrame(DetRow, RowLayout); + TGLabel* DetLabel = new TGLabel(DetRow, "Detector"); + DetRow->AddFrame(DetLabel, LabelLayout); + m_DetectorSelector = new TGComboBox(DetRow); + for (auto id : m_DetIDs) { + m_DetectorSelector->AddEntry(TString::Format("%d", id), id); + } + m_DetectorSelector->Select(m_SelectedDetector); + m_DetectorSelector->Resize(50, 20); + //m_DetectorSelector->Connect("Selected(Int_t)", "MGUIExpoPlotTacDiff", this, "OnDetectorSelected(Int_t)"); + DetRow->AddFrame(m_DetectorSelector, EntryLayout); + + // Side + Strip row + TGHorizontalFrame* StripRow = new TGHorizontalFrame(ControlFrame); + ControlFrame->AddFrame(StripRow, RowLayout); + m_SideSelector = new TGComboBox(StripRow); + m_SideSelector->AddEntry("LV", 0); + m_SideSelector->AddEntry("HV", 1); + m_SideSelector->Select(0); + m_SideSelector->Resize(45, 20); + //m_SideSelector->Connect("Selected(Int_t)", "MGUIExpoPlotTacDiff", this, "OnSideSelected(Int_t)"); + StripRow->AddFrame(m_SideSelector, EntryLayout); + m_StripEntry = new TGNumberEntry(StripRow, 0, 3, -1, + TGNumberFormat::kNESInteger, + TGNumberFormat::kNEANonNegative, + TGNumberFormat::kNELLimitMinMax, 0, 62); + //m_StripEntry->Connect("ValueSet(Long_t)", "MGUIExpoPlotTacDiff", this, "OnStripSelected()"); + m_StripEntry->GetNumberEntry()->Connect("ReturnPressed()", "MGUIExpoPlotTacDiff", this, "OnStripSelected()"); + m_StripEntry->Resize(50, 20); + StripRow->AddFrame(m_StripEntry, EntryLayout); + + // Energy min row + TGHorizontalFrame* EMinRow = new TGHorizontalFrame(ControlFrame); + ControlFrame->AddFrame(EMinRow, RowLayout); + TGLabel* EMinLabel = new TGLabel(EMinRow, "Min E"); + EMinRow->AddFrame(EMinLabel, LabelLayout); + m_EnergyMinEntry = new TGNumberEntry(EMinRow, 0, 5, -1, + TGNumberFormat::kNESRealTwo, + TGNumberFormat::kNEANonNegative); + //m_EnergyMinEntry->Connect("ValueSet(Long_t)", "MGUIExpoPlotTacDiff", this, "OnEnergyRangeChanged()"); + m_EnergyMinEntry->GetNumberEntry()->Connect("ReturnPressed()", "MGUIExpoPlotTacDiff", this, "OnEnergyRangeChanged()"); + m_EnergyMinEntry->Resize(70, 20); + EMinRow->AddFrame(m_EnergyMinEntry, EntryLayout); + + // Energy max row + TGHorizontalFrame* EMaxRow = new TGHorizontalFrame(ControlFrame); + ControlFrame->AddFrame(EMaxRow, RowLayout); + TGLabel* EMaxLabel = new TGLabel(EMaxRow, "Max E"); + EMaxRow->AddFrame(EMaxLabel, LabelLayout); + m_EnergyMaxEntry = new TGNumberEntry(EMaxRow, 10000, 5, -1, + TGNumberFormat::kNESRealTwo, + TGNumberFormat::kNEANonNegative); + //m_EnergyMaxEntry->Connect("ValueSet(Long_t)", "MGUIExpoPlotTacDiff", this, "OnEnergyRangeChanged()"); + m_EnergyMaxEntry->GetNumberEntry()->Connect("ReturnPressed()", "MGUIExpoPlotTacDiff", this, "OnEnergyRangeChanged()"); + m_EnergyMaxEntry->Resize(70, 20); + EMaxRow->AddFrame(m_EnergyMaxEntry, EntryLayout); + + // === UPDATE BUTTON === + m_UpdateSelectionButton = new TGTextButton(ControlFrame, "Update Plot Selection", c_UpdateSelection); + m_UpdateSelectionButton->Associate(this); + ControlFrame->AddFrame(m_UpdateSelectionButton, RowLayout); + + // Square canvas + m_DtacVsDtacCanvas = new TRootEmbeddedCanvas("DtacVsDtac", BottomRow, 400, 400); + BottomRow->AddFrame(m_DtacVsDtacCanvas, ExpandLayout); + + RedrawPlots(); + + m_IsCreated = true; + m_Mutex.UnLock(); +} + +//////////////////////////////////////////////////////////////////////////////// + +bool MGUIExpoPlotTacDiff::ProcessMessage(long Message, long Parameter1, long Parameter2) +{ + switch (GET_MSG(Message)) { + case kC_COMMAND: + switch (GET_SUBMSG(Message)) { + case kCM_BUTTON: + switch (Parameter1) { + case c_UpdateSelection: + OnUpdateSelection(); + return true; + default: + break; + } + break; + default: + break; + } + break; + default: + break; + } + return true; +} + +//////////////////////////////////////////////////////////////////////////////// + +void MGUIExpoPlotTacDiff::OnUpdateSelection() +{ + m_Mutex.Lock(); + + // Read current widget values + m_SelectedDetector = m_DetectorSelector->GetSelected(); + m_SelectedSide = m_SideSelector->GetSelected(); + m_SelectedStrip = (int)m_StripEntry->GetNumber(); + m_EnergyMin = m_EnergyMinEntry->GetNumber(); + m_EnergyMax = m_EnergyMaxEntry->GetNumber(); + + RedrawPlots(); + + m_Mutex.UnLock(); +} +//////////////////////////////////////////////////////////////////////////////// + +void MGUIExpoPlotTacDiff::RedrawPlots() +{ + if (m_SelectedSide == 0) { + m_SelectedStripPairCode = 10000 * m_SelectedDetector + 100 * m_SelectedStrip + 99; + } else { + m_SelectedStripPairCode = 10000 * m_SelectedDetector + m_SelectedStrip + 9900; + } + + int key = m_SelectedStripPairCode; + if (m_DtacVsDepthHistograms.count(key) > 0) { + m_DtacVsDepthCanvas->GetCanvas()->cd(); + m_DtacVsDepthHistograms[key]->Draw("colz"); + m_DtacVsDepthCanvas->GetCanvas()->Modified(); + m_DtacVsDepthCanvas->GetCanvas()->Update(); + + m_DtacVsFracCanvas->GetCanvas()->cd(); + m_DtacVsFracHistograms[key]->Draw("colz"); + m_DtacVsFracCanvas->GetCanvas()->Modified(); + m_DtacVsFracCanvas->GetCanvas()->Update(); + + m_DtacVsDtacCanvas->GetCanvas()->cd(); + m_DtacVsDtacHistograms[key]->Draw("colz"); + m_DtacVsDtacCanvas->GetCanvas()->Modified(); + m_DtacVsDtacCanvas->GetCanvas()->Update(); + } +} + +//////////////////////////////////////////////////////////////////////////////// + + +void MGUIExpoPlotTacDiff::Update() +{ + m_Mutex.Lock(); + if (m_IsCreated) { + RedrawPlots(); + } + m_Mutex.UnLock(); +} + +//////////////////////////////////////////////////////////////////////////////// + +void MGUIExpoPlotTacDiff::Export(const MString& FileName) +{ + m_Mutex.Lock(); + + int key = m_SelectedStripPairCode; + TCanvas* P = new TCanvas("", "", 1200, 800); + P->Divide(2, 2); + + if (m_DtacVsDepthHistograms.count(key) > 0) { + P->cd(1); + m_DtacVsDepthHistograms[key]->DrawCopy("colz"); + P->cd(2); + m_DtacVsFracHistograms[key]->DrawCopy("colz"); + P->cd(3); + m_DtacVsDtacHistograms[key]->DrawCopy("colz"); + } + + P->SaveAs(FileName); + delete P; + + m_Mutex.UnLock(); +} diff --git a/src/MModuleDepthCalibration.cxx b/src/MModuleDepthCalibration.cxx index 33d5e93f..c2a79d43 100644 --- a/src/MModuleDepthCalibration.cxx +++ b/src/MModuleDepthCalibration.cxx @@ -163,19 +163,35 @@ void MModuleDepthCalibration::CreateExpos() // Create all expos if (HasExpos() == true) return; + + cout << "CreateExpos: m_DetectorIDs.size() = " << m_DetectorIDs.size() << endl; + + // dTAC expo for bump diagnostics + m_ExpoPlotTacDiff = new MGUIExpoPlotTacDiff(this); + for (unsigned int i = 0; i < m_DetectorIDs.size(); ++i) { + unsigned int DetID = m_DetectorIDs[i]; + double thickness = m_Thicknesses[DetID]; + m_ExpoPlotTacDiff->SetHistogramParameters(DetID, + 120, -thickness/2.0, thickness/2.0, // depth bins + 200, -100, 100, // dTac bins + 100, 0, 1); // fraction bins + } + m_Expos.push_back(m_ExpoPlotTacDiff); - // Set the histogram display + // Depth calibration expo m_ExpoDepthCalibration = new MGUIExpoDepthCalibration(this); - m_ExpoDepthCalibration->SetDepthHistogramArrangement(&m_DetectorIDs); - for (unsigned int i = 0; i < m_DetectorIDs.size(); ++i){ + // Set parameters first (so m_NBins/m_Min/m_Max are available) + for (unsigned int i = 0; i < m_DetectorIDs.size(); ++i) { unsigned int DetID = m_DetectorIDs[i]; double thickness = m_Thicknesses[DetID]; - m_ExpoDepthCalibration->SetDepthHistogramParameters(DetID, 120, -thickness/2.0,thickness/2.0); + cout << " DetID=" << DetID << " thickness=" << thickness << endl; + m_ExpoDepthCalibration->SetDepthHistogramParameters(DetID, 120, -thickness/2.0, thickness/2.0); } + // Then create arrangement (which uses m_NBins/m_Min/m_Max to create histograms) + m_ExpoDepthCalibration->SetDepthHistogramArrangement(&m_DetectorIDs); m_Expos.push_back(m_ExpoDepthCalibration); } - ///////////////////////////////////////////////////////////////////////////////// @@ -195,6 +211,7 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) // 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); + double HitEnergy = H->GetEnergy(); int Grade = GetHitGrade(H); @@ -327,11 +344,20 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) else { // FR TODO the last boolean is for sean's weighting method; make it a flag - auto [rawZpos, rawZsigma] = CalculateZfromCTD(rawCTD_s, noise,DetID, Grade, true); - + auto [rawZpos, rawZsigma] = CalculateZfromCTD(rawCTD_s, noise,DetID, Grade, false); // true (sean weighting) + // add the raw depth to the raw depth histogram + if (Event->HasStripPairingError()==false) { + if (HasExpos() == true) { + m_ExpoDepthCalibration->AddRawDepth(DetID, LVStripID, HVStripID, rawZpos); + } + } + // TODO depth correction loop! // // step 1 -- check the HV side: + bool ChargeSharingHV = false; + vector CorrectedHVTiming; + vector CorrectedHVTimingUncertainty; //int StripPairCode = 10000*DetID + HVStripID; // note, it is the lower strip ID always (eg, 15 if sharing between strip 15 and 16) // -- how many strips share > 10% of the total energy (or are over slow threshold, maybe?) need this info for next steps // -- check dTAC between adjacent strips with charge sharing and also strips relative to their low-eneryg neighbors @@ -342,12 +368,11 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) // -- note, rawZpos is used in the above calculations! // step 2 -- check the LV side: - bool ValidatedTiming = false; // should put a flag here eventually TODO; if NN do not have fast timing - bool ZombieBump = false; - vector CorrectedTiming; - vector CorrectedTimingUncertainty; - bool CorrectedChargeSharingLV = false; - // also should put a flag here TODO (is there another flag for charge sharing?) + bool ValidatedLVTiming = false; // should put a flag here eventually TODO; if NN do not have fast timing + bool ZombieBump = false;// should put a flag here TODO + vector CorrectedLVTiming; + vector CorrectedLVTimingUncertainty; + bool ChargeSharingLV = false;// also should put a flag here TODO (is there another flag for charge sharing?) if (LVEnergyFraction > m_SingleStripChargeSharing){// if we have one obvious main strip, we are not going to be correcting charge sharing but just checking for zombie bump // compare with the neighbors, if possible @@ -355,43 +380,59 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) int pm = 2*neighbor - 1; // -1 for neighbor 0 (neighbor is left); + 1 for neighbor == 1 (right neighbor, which is nominal for the convention StripPairID = left StripID of pair int NeighborStripID = LVStripID + pm; MStripHit* NSH = GetStrip(LVStrips, NeighborStripID); - if (NeighborStripID >=0 && NeighborStripID <=63 && NSH){ // neighbor is not a guard ring strip, NSH exists (not a null pointer) TODO and NSH has fast timing! + if (NeighborStripID >=0 && NeighborStripID <=63 && NSH && NSH->HasFastTiming()){ // neighbor is not a guard ring strip, NSH exists (not a null pointer) and has fast timing! double dTacData = (LVSH->GetTiming() - NSH->GetTiming())*pm; // always the left strip - right strip; neighbor on left means pm = -1 -> NSH - LVSH timing double fracData = (NSH->GetEnergy()/(LVSH->GetEnergy() + NSH->GetEnergy())*pm) + 1 - neighbor; // always the fraction on the right stripHit; nominally NSH for right neighbor - int StripPairCode = 10000*DetID + 100*(LVStripID - 1 + neighbor); //LVStripID -1 + 0 = LVStripID -1 (left neighbor); or LVStripID -1 + 1 = LVStripID (LVStrip is the StripID when we consider right negihbor) - int x = (fracData - 0.5); // x is LVEnergyFraction - 0.5 - vector* CSPolyCoeffs = GetChargeSharingPolyCoeffsLV(DetID,rawZpos); // TODO need to actually check and fill the variable that checks the length of this, and check that it's right when loading - vector* CSCoeffs = GetChargeSharingCoeffs(StripPairCode,rawZpos); - if (CSCoeffs && CSPolyCoeffs){ - double dTacExpect = (CSPolyCoeffs->at(0)*x + CSPolyCoeffs->at(1)*x*x*x)*CSCoeffs->at(0) + CSCoeffs->at(1);// TODO update if not cubic polynomial + int StripPairCode = 10000*DetID + 100*(LVStripID - 1 + neighbor) + 99; //LVStripID -1 + 0 = LVStripID -1 (left neighbor); or LVStripID -1 + 1 = LVStripID (LVStrip is the StripID when we consider right negihbor) + double x = (fracData - 0.5); // x is LVEnergyFraction - 0.5, ie the parameter of Isidro's polynomials, which are forced to go through (0.5, 0) + vector CSPolyCoeffs = GetChargeSharingPolyCoeffsLV(DetID,rawZpos); // TODO need to actually check and fill the variable that checks the length of this, and check that it's right when loading + vector CSCoeffs = GetChargeSharingCoeffs(StripPairCode,rawZpos); + + // add to the expo + if (HasExpos() == true) { + m_ExpoPlotTacDiff->AddData(StripPairCode, rawZpos, dTacData, fracData, std::nan(""),HitEnergy); + } + + if (!CSCoeffs.empty() && !CSPolyCoeffs.empty()){ + double dTacExpect = (CSPolyCoeffs.at(0)*x + CSPolyCoeffs.at(1)*x*x*x)*CSCoeffs.at(0) + CSCoeffs.at(1);// TODO update if not cubic polynomial // TODO we should display (dTacData - dTacExpect)*pm to keep an eye on the prevalence of the bump - if (rawZpos > -5 && dTacData < 500 && dTacData > -500 && (dTacData - dTacExpect)*pm > 2*noise){// zombie bump! TODO update bump criteria, fix noise, remove Slow Timing check once checked earlier + if (rawZpos > -0.5 && (dTacData - dTacExpect)*pm > 2*noise){// zombie bump! TODO update bump criteria, fix noise ZombieBump = true; - ValidatedTiming = false; + ValidatedLVTiming = false; + cout << "LV Strip: "< -1*pm*dTacExpect) { // zombie bump with good neighbor. Can we not always do this in this case, though, since we know what the timing should be? - CorrectedTiming.push_back(LVTiming-(dTacData - dTacExpect)*pm); - CorrectedTimingUncertainty.push_back(noise*2); // to do: quantify and make into something real + CorrectedLVTiming.push_back(LVTiming-60); // needs to be a config + CorrectedLVTimingUncertainty.push_back(noise*2); // to do: quantify and make into something real + } else { + CorrectedLVTiming.push_back(0); // needs to be a config + CorrectedLVTimingUncertainty.push_back(0); // to do: quantify and make into something real } } else { - if (dTacData < 500 && dTacData > -500 && !ZombieBump && CSCoeffs && CSPolyCoeffs) ValidatedTiming = true; + if (!ZombieBump && !CSCoeffs.empty() && !CSPolyCoeffs.empty()) ValidatedLVTiming = true; } } } } } else { // charge sharing correction - CorrectedChargeSharingLV = true; + ChargeSharingLV = true; } - if (ZombieBump) m_ZombieBump++; - if (CorrectedTiming.size() > 0){// if we have a correction. Also, implemented weighting! TODO an + // correct the timing + if (CorrectedLVTiming.size() > 0){// if we have a correction. Also, implemented weighting! TODO an double correctionSum = 0; - for (unsigned int j = 0; j < CorrectedTiming.size(); j++) { + for (unsigned int j = 0; j < CorrectedLVTiming.size(); j++) { // TODO check that they are consistent and drop one if not.... - correctionSum += CorrectedTiming.at(j); + correctionSum += CorrectedLVTiming.at(j); } - LVTiming = correctionSum / CorrectedTiming.size(); + LVTiming = correctionSum / CorrectedLVTiming.size(); } - rawCTD = (HVTiming - LVTiming); - rawCTD_s = (rawCTD - Coeffs->at(1))/(Coeffs->at(0)); + + // update counters + if (ZombieBump) m_ZombieBump++; + if (ChargeSharingLV) m_ChargeSharingLV++; + if (ChargeSharingHV) m_ChargeSharingHV++; + + double CTD = (HVTiming - LVTiming); + double CTD_s = (CTD - Coeffs->at(1))/(Coeffs->at(0)); // -- how many strips share > 10% of the total energy (or are over the slow threshold, maybe?) need this info for the next steps // -- check dTAC between adjacent strips with charge sharing and also strips relative to their low-energy neighbors // -- deal with the zombie bump! @@ -400,27 +441,24 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) // -- and then also calculate the corrected LV tac as the weighted average, with its own uncertainty (they should be consistent) // -- note, rawZpos is used in teh above calculations - // step 3 -- calculate new corrected CTD, and from that calculate new corrected Z - // + // step 3 // bonus points: implement x and y localization based on info in step 1 and 2 with charge sharing :) - // bonus points -- can re-calculated depth between 1 and 2 if significant changes with HV correction would impact ZB correction! - if ((rawCTD_s < (Xmin - 2.0*noise)) || (rawCTD_s > (Xmax + 2.0*noise))) { + if ((CTD_s < (Xmin - 2.0*noise)) || (CTD_s > (Xmax + 2.0*noise))) { H->SetNoDepth(); Event->SetDepthCalibrationError("Out of Range"); ++m_Error2; } + // If the CTD is in range, calculate the depth else { - auto [rawZpos, rawZsigma] = CalculateZfromCTD(rawCTD_s, noise,DetID, Grade, true); - Zpos = rawZpos; - Zsigma = rawZsigma; + auto [Zpos, Zsigma] = CalculateZfromCTD(CTD_s, noise,DetID, Grade, false); // Add the depth to the GUI histogram. if (Event->HasStripPairingError()==false) { if (HasExpos() == true) { - m_ExpoDepthCalibration->AddDepth(DetID, Zpos); + m_ExpoDepthCalibration->AddDepth(DetID, LVStripID, HVStripID, Zpos); } } m_NoError+=1; @@ -458,6 +496,11 @@ std::tuple MModuleDepthCalibration::CalculateZfromCTD(double CTD { vector CTDVec = GetCTD(DetID, Grade); vector DepthVec = GetDepth(DetID); + + if (CTDVec.empty() || DepthVec.empty()) { + cout << "NO CTD Vector for this detector and GRADE!!!!" < MModuleDepthCalibration::CalculateZfromCTD(double CTD } // otherwise, use the standard appropach with no rounding off // if out of bounds, return boundary - if (CTDvalue <= CTDVec.front()) { + double Vecmin = CTDVec.front(); + double Vecmax = CTDVec.back(); + if (Vecmin > Vecmax) cout << "CTD vec in descenting order!! front: "<< Vecmin << " back: "<= Vecmax) { + return std::make_tuple(DepthVec.front(), DepthVec.back() - DepthVec.front()); + } auto it = std::upper_bound(CTDVec.begin(), CTDVec.end(), CTD_high); unsigned int i = std::distance(CTDVec.begin(), it); + if (i == 0) return std::make_tuple(DepthVec.front(), 0.0); + if (i >= CTDVec.size()) return std::make_tuple(DepthVec.front(), DepthVec.back() - DepthVec.front()); double fraction = (CTD_high - CTDVec[i - 1]) / (CTDVec[i] - CTDVec[i - 1]); double depth_high = DepthVec[i - 1] + fraction * (DepthVec[i] - DepthVec[i - 1]); return std::make_tuple(DepthVec.front(),depth_high-DepthVec.front()); } - if (CTDvalue >= CTDVec.back()) { + if (CTDvalue >= Vecmax) { double CTD_low = CTDvalue - noise/2.355; if (CTD_low >= CTDVec.back()) { return std::make_tuple(DepthVec.back(), 0.0); } + if (CTD_low <= Vecmin) { + return std::make_tuple(DepthVec.back(), DepthVec.back() - DepthVec.front()); + } auto it = std::upper_bound(CTDVec.begin(), CTDVec.end(), CTD_low); unsigned int i = std::distance(CTDVec.begin(), it); + if (i == 0) return std::make_tuple(DepthVec.back(), DepthVec.back() - DepthVec.front()); + if (i >= CTDVec.size()) return std::make_tuple(DepthVec.back(), 0.0); double fraction = (CTD_low - CTDVec[i - 1]) / (CTDVec[i] - CTDVec[i - 1]); double depth_low = DepthVec[i - 1] + fraction * (DepthVec[i] - DepthVec[i - 1]); return std::make_tuple(DepthVec.back(),depth_low-DepthVec.back()); @@ -516,20 +573,29 @@ std::tuple MModuleDepthCalibration::CalculateZfromCTD(double CTD // if not out of bounds, extrapolate and calculate errors.... auto it = std::upper_bound(CTDVec.begin(), CTDVec.end(), CTDvalue); unsigned int i = std::distance(CTDVec.begin(), it); + if (i == 0) return std::make_tuple(DepthVec.front(), 0.0); + if (i >= CTDVec.size()) return std::make_tuple(DepthVec.back(), 0.0); double fraction = (CTDvalue - CTDVec[i - 1]) / (CTDVec[i] - CTDVec[i - 1]); double depth = DepthVec[i - 1] + fraction * (DepthVec[i] - DepthVec[i - 1]); double CTD_low = std::max(CTDvalue - noise/2.355,CTDVec.front()); double CTD_high = std::min(CTDvalue + noise/2.355,CTDVec.back()); + + double depth_low = depth; + double depth_high = depth; + it = std::upper_bound(CTDVec.begin(), CTDVec.end(), CTD_low); i = std::distance(CTDVec.begin(), it); - fraction = (CTD_low - CTDVec[i - 1]) / (CTDVec[i] - CTDVec[i - 1]); - double depth_low = DepthVec[i - 1] + fraction * (DepthVec[i] - DepthVec[i - 1]); + if (i > 0 && i < CTDVec.size()) { + fraction = (CTD_low - CTDVec[i - 1]) / (CTDVec[i] - CTDVec[i - 1]); + depth_low = DepthVec[i - 1] + fraction * (DepthVec[i] - DepthVec[i - 1]); + }// give an else TODO it = std::upper_bound(CTDVec.begin(), CTDVec.end(), CTD_high); i = std::distance(CTDVec.begin(), it); - fraction = (CTD_high - CTDVec[i - 1]) / (CTDVec[i] - CTDVec[i - 1]); - double depth_high = DepthVec[i - 1] + fraction * (DepthVec[i] - DepthVec[i - 1]); - + if (i > 0 && i < CTDVec.size()) { + fraction = (CTD_high - CTDVec[i - 1]) / (CTDVec[i] - CTDVec[i - 1]); + depth_high = DepthVec[i - 1] + fraction * (DepthVec[i] - DepthVec[i - 1]); + } return std::make_tuple(depth, (depth_high - depth_low) / 2.); } @@ -675,7 +741,7 @@ bool MModuleDepthCalibration::LoadChargeSharingConfigFile(MString FileName) // fill m_ChargeSharingDepths for each detector vector depths; - for (double i = -7.; i < 7.4; i = i + 1.) depths.push_back(i); + for (double i = -7.; i < 7.4; i = i + 1.) depths.push_back(i/10); // should be 100% in cm m_ChargeSharingDepths[DetID] = depths; for(int z = 0; z < depths.size(); z++){ @@ -689,11 +755,11 @@ bool MModuleDepthCalibration::LoadChargeSharingConfigFile(MString FileName) vector coeffs; // the stretch and offset, currently set to the same values for all strips (which are almost certainly wrong) coeffs.push_back(1.033454449); coeffs.push_back(-1.996517705);coeffs.push_back(6.373983606); // stretch, offset, dTacSigma. Need to figure out how to deal with dTac sigma in an energy-dependent way for (int LVStripID = 0; LVStripID < 63; LVStripID++){ // up to 62 since these are pairs - int StripPairCode = 10000*DetID + 100*LVStripID; + int StripPairCode = 10000*DetID + 100*LVStripID + 99; m_ChargeSharingCoeffs[StripPairCode].push_back(coeffs); } for (int HVStripID = 0; HVStripID < 63; HVStripID++){ - int StripPairCode = 10000*DetID + HVStripID; + int StripPairCode = 10000*DetID + HVStripID + 9900; m_ChargeSharingCoeffs[StripPairCode].push_back(coeffs); } } @@ -772,7 +838,7 @@ bool MModuleDepthCalibration::LoadCoeffsFile(MString FileName) ///////////////////////////////////////////////////////////////////////////////// -std::vector* MModuleDepthCalibration::GetChargeSharingCoeffs(int StripPairCode, double z) +std::vector MModuleDepthCalibration::GetChargeSharingCoeffs(int StripPairCode, double z) { int DetID = StripPairCode / 10000; // Check to see if the charge sharing coefficients have been loaded. If so, try to get the coefficients for the specified strip pair. @@ -780,95 +846,95 @@ std::vector* MModuleDepthCalibration::GetChargeSharingCoeffs(int StripPa if (m_ChargeSharingCoeffs.count(StripPairCode) > 0) { // if we only sampled one depth, or we're beyond the depth range, just return the closest coefficients - if (z <= m_ChargeSharingDepths[DetID].front()) return &m_ChargeSharingCoeffs[StripPairCode].at(0); - if (z >= m_ChargeSharingDepths[DetID].back()) return &m_ChargeSharingCoeffs[StripPairCode].at(m_ChargeSharingDepths[DetID].size()-1); + if (z <= m_ChargeSharingDepths[DetID].front()) return m_ChargeSharingCoeffs[StripPairCode].at(0); + if (z >= m_ChargeSharingDepths[DetID].back()) return m_ChargeSharingCoeffs[StripPairCode].at(m_ChargeSharingDepths[DetID].size()-1); // otherwise, interpolate for (unsigned int i = 0; i < m_ChargeSharingDepths[DetID].size() - 1; i++){ if (z >= m_ChargeSharingDepths[DetID].at(i) && z < m_ChargeSharingDepths[DetID].at(i + 1)) { double f = (z - m_ChargeSharingDepths[DetID][i]) / (m_ChargeSharingDepths[DetID][i + 1] - m_ChargeSharingDepths[DetID][i]); - m_InterpolatedCoeffs.clear(); - for (int j = 0; j < m_ChargeSharingCoeffs[StripPairCode].at(i).size(); j++) m_InterpolatedCoeffs.push_back((1.0 - f) * m_ChargeSharingCoeffs[StripPairCode][i][j] + f * m_ChargeSharingCoeffs[StripPairCode][i + 1][j]); - return &m_InterpolatedCoeffs; + vector result; + for (unsigned int j = 0; j < m_ChargeSharingCoeffs[StripPairCode].at(i).size(); j++) result.push_back((1.0 - f) * m_ChargeSharingCoeffs[StripPairCode][i][j] + f * m_ChargeSharingCoeffs[StripPairCode][i + 1][j]); + return result; } } } else { if (g_Verbosity >= c_Warning) { cout << "MModuleDepthCalibration::GetChargeSharingCoeffs: cannot get charge sharing coefficients; strip pair code " << StripPairCode << " not found." << endl; } - return nullptr; + return {}; } } else { cout << "MModuleDepthCalibration::GetChargeSharingCoeffs: cannot get charge sharing coefficients; file has not yet been loaded." << endl; - return nullptr; + return {}; } } ///////////////////////////////////////////////////////////////////////////////// -std::vector* MModuleDepthCalibration::GetChargeSharingPolyCoeffsLV(int DetID, double z) +std::vector MModuleDepthCalibration::GetChargeSharingPolyCoeffsLV(int DetID, double z) { // Check to see if the charge sharing coefficients have been loaded. If so, try to get the coefficients for the specified strip pair. if (m_ChargeSharingConfigFileIsLoaded == true) { if (m_ChargeSharingPolyCoeffsLV.count(DetID) > 0) { // if we only sampled one depth, or we're beyond the depth range, just return the closest coefficients - if (z <= m_ChargeSharingDepths[DetID].front()) return &m_ChargeSharingPolyCoeffsLV[DetID].at(0); - if (z >= m_ChargeSharingDepths[DetID].back()) return &m_ChargeSharingPolyCoeffsLV[DetID].at(m_ChargeSharingDepths[DetID].size()-1); + if (z <= m_ChargeSharingDepths[DetID].front()) return m_ChargeSharingPolyCoeffsLV[DetID].at(0); + if (z >= m_ChargeSharingDepths[DetID].back()) return m_ChargeSharingPolyCoeffsLV[DetID].at(m_ChargeSharingDepths[DetID].size()-1); // otherwise, interpolate for (unsigned int i = 0; i < m_ChargeSharingDepths[DetID].size() - 1; i++){ if (z >= m_ChargeSharingDepths[DetID].at(i) && z < m_ChargeSharingDepths[DetID].at(i + 1)) { double f = (z - m_ChargeSharingDepths[DetID][i]) / (m_ChargeSharingDepths[DetID][i + 1] - m_ChargeSharingDepths[DetID][i]); - m_InterpolatedCoeffs.clear(); - for (int j = 0; j < m_ChargeSharingPolyCoeffsLV[DetID].at(i).size(); j++) m_InterpolatedCoeffs.push_back((1.0 - f) * m_ChargeSharingPolyCoeffsLV[DetID][i][j] + f * m_ChargeSharingPolyCoeffsLV[DetID][i + 1][j]); - return &m_InterpolatedCoeffs; + vector result; + for (int j = 0; j < m_ChargeSharingPolyCoeffsLV[DetID].at(i).size(); j++) result.push_back((1.0 - f) * m_ChargeSharingPolyCoeffsLV[DetID][i][j] + f * m_ChargeSharingPolyCoeffsLV[DetID][i + 1][j]); + return result; } } } else { if (g_Verbosity >= c_Warning) { cout << "MModuleDepthCalibration::GetChargeSharingPolyCoeffsLV: cannot get charge sharing polynomial coefficients; detector id code " << DetID << " not found." << endl; } - return nullptr; + return {}; } } else { cout << "MModuleDepthCalibration::GetChargeSharingPolyCoeffsLV: cannot get charge sharing coefficients; file has not yet been loaded." << endl; - return nullptr; + return {}; } } ///////////////////////////////////////////////////////////////////////////////// -std::vector* MModuleDepthCalibration::GetChargeSharingPolyCoeffsHV(int DetID, double z) +std::vector MModuleDepthCalibration::GetChargeSharingPolyCoeffsHV(int DetID, double z) { // Check to see if the charge sharing coefficients have been loaded. If so, try to get the coefficients for the specified strip pair. if (m_ChargeSharingConfigFileIsLoaded == true) { if (m_ChargeSharingPolyCoeffsHV.count(DetID) > 0) { // if we only sampled one depth, or we're beyond the depth range, just return the closest coefficients - if (z <= m_ChargeSharingDepths[DetID].front()) return &m_ChargeSharingPolyCoeffsHV[DetID].at(0); - if (z >= m_ChargeSharingDepths[DetID].back()) return &m_ChargeSharingPolyCoeffsHV[DetID].at(m_ChargeSharingDepths[DetID].size()-1); + if (z <= m_ChargeSharingDepths[DetID].front()) return m_ChargeSharingPolyCoeffsHV[DetID].at(0); + if (z >= m_ChargeSharingDepths[DetID].back()) return m_ChargeSharingPolyCoeffsHV[DetID].at(m_ChargeSharingDepths[DetID].size()-1); // otherwise, interpolate for (unsigned int i = 0; i < m_ChargeSharingDepths[DetID].size() - 1; i++){ if (z >= m_ChargeSharingDepths[DetID].at(i) && z < m_ChargeSharingDepths[DetID].at(i + 1)) { double f = (z - m_ChargeSharingDepths[DetID][i]) / (m_ChargeSharingDepths[DetID][i + 1] - m_ChargeSharingDepths[DetID][i]); - m_InterpolatedCoeffs.clear(); - for (int j = 0; j < m_ChargeSharingPolyCoeffsHV[DetID].at(i).size(); j++) m_InterpolatedCoeffs.push_back((1.0 - f) * m_ChargeSharingPolyCoeffsHV[DetID][i][j] + f * m_ChargeSharingPolyCoeffsHV[DetID][i + 1][j]); - return &m_InterpolatedCoeffs; + vector result; + for (int j = 0; j < m_ChargeSharingPolyCoeffsHV[DetID].at(i).size(); j++) result.push_back((1.0 - f) * m_ChargeSharingPolyCoeffsHV[DetID][i][j] + f * m_ChargeSharingPolyCoeffsHV[DetID][i + 1][j]); + return result; } } } else { if (g_Verbosity >= c_Warning) { cout << "MModuleDepthCalibration::GetChargeSharingPolyCoeffsHV: cannot get charge sharing polynomial coefficients; detector id code " << DetID << " not found." << endl; } - return nullptr; + return {}; } } else { cout << "MModuleDepthCalibration::GetChargeSharingPolyCoeffsHV: cannot get charge sharing coefficients; file has not yet been loaded." << endl; - return nullptr; + return {}; } } @@ -939,6 +1005,7 @@ bool MModuleDepthCalibration::LoadSplinesFile(MString FileName) if (DepthVec.size() > 0) { Result &= AddDepthCTD(DepthVec, CTDArr, DetID, m_DepthGrid, m_CTDMap, m_SplineMap, 1000); + cout << "loaded spline for detector"< 0) { Result &= AddDepthCTD(DepthVec, CTDArr, DetID, m_DepthGrid, m_CTDMap, m_SplineMap, 1000); + cout << "loaded spline for detector"< Date: Thu, 23 Jul 2026 21:27:52 -0700 Subject: [PATCH 09/15] better labels for plots in dTAC Expo --- src/MGUIExpoPlotTacDiff.cxx | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/MGUIExpoPlotTacDiff.cxx b/src/MGUIExpoPlotTacDiff.cxx index 300ca631..32b5a46e 100644 --- a/src/MGUIExpoPlotTacDiff.cxx +++ b/src/MGUIExpoPlotTacDiff.cxx @@ -77,19 +77,19 @@ void MGUIExpoPlotTacDiff::SetHistogramParameters(unsigned int DetID, unsigned in for (auto id : m_DetIDs) { if (id == DetID) { found = true; break; } } if (!found) m_DetIDs.push_back(DetID); - // LV strips: code = 10000*DetID + 100*StripID + 99 + // LV strips: code = 10000*DetID + 100*StripID + 99 for (int s = 0; s < 63; s++) { int key = 10000 * DetID + 100 * s + 99; if (m_DtacVsDepthHistograms.count(key) == 0) { m_DtacVsDepthHistograms[key] = new TH2D("", TString::Format("Det %d LV %d: dTAC vs Depth", DetID, s), NBinsDepth, DepthMin, DepthMax, NBinsDtac, DtacMin, DtacMax); m_DtacVsDepthHistograms[key]->SetXTitle("Depth [cm]"); - m_DtacVsDepthHistograms[key]->SetYTitle("dTAC [ns]"); + m_DtacVsDepthHistograms[key]->SetYTitle(TString::Format("TAC %s - TAC %s [ns]",s,s+1)); - m_DtacVsFracHistograms[key] = new TH2D("", TString::Format("Det %d LV %d: dTAC vs Frac", DetID, s), + m_DtacVsFracHistograms[key] = new TH2D("", TString::Format("Det %d Charge Shared between LV %d and %d", DetID, s,s+1), NBinsFrac, FracMin, FracMax, NBinsDtac, DtacMin, DtacMax); - m_DtacVsFracHistograms[key]->SetXTitle("Fraction"); - m_DtacVsFracHistograms[key]->SetYTitle("dTAC [ns]"); + m_DtacVsFracHistograms[key]->SetXTitle(TString::Format("Charge Sharing Fraction: (strip %d) / (strip %d + strip %d",s,s,s+1)); + m_DtacVsFracHistograms[key]->SetYTitle(TString::Format("TAC %s - TAC %s [ns]",s,s+1)); m_DtacVsDtacHistograms[key] = new TH2D("", TString::Format("Det %d LV %d: dTAC vs dTAC", DetID, s), NBinsDtac, DtacMin, DtacMax, NBinsDtac, DtacMin, DtacMax); @@ -105,12 +105,12 @@ void MGUIExpoPlotTacDiff::SetHistogramParameters(unsigned int DetID, unsigned in m_DtacVsDepthHistograms[key] = new TH2D("", TString::Format("Det %d HV %d: dTAC vs Depth", DetID, s), NBinsDepth, DepthMin, DepthMax, NBinsDtac, DtacMin, DtacMax); m_DtacVsDepthHistograms[key]->SetXTitle("Depth [cm]"); - m_DtacVsDepthHistograms[key]->SetYTitle("dTAC [ns]"); + m_DtacVsDepthHistograms[key]->SetYTitle(TString::Format("TAC %s - TAC %s [ns]",s,s+1)); - m_DtacVsFracHistograms[key] = new TH2D("", TString::Format("Det %d HV %d: dTAC vs Frac", DetID, s), + m_DtacVsFracHistograms[key] = new TH2D("", TString::Format("Det %d Charge Shared between LV %s and %s", DetID, s,s+1), NBinsFrac, FracMin, FracMax, NBinsDtac, DtacMin, DtacMax); - m_DtacVsFracHistograms[key]->SetXTitle("Fraction"); - m_DtacVsFracHistograms[key]->SetYTitle("dTAC [ns]"); + m_DtacVsFracHistograms[key]->SetXTitle(TString::Format("Charge Sharing Fraction: (strip %d) / (strip %d + strip %d",s,s,s+1)); + m_DtacVsFracHistograms[key]->SetYTitle(TString::Format("TAC %s - TAC %s [ns]",s,s+1)); m_DtacVsDtacHistograms[key] = new TH2D("", TString::Format("Det %d HV %d: dTAC vs dTAC", DetID, s), NBinsDtac, DtacMin, DtacMax, NBinsDtac, DtacMin, DtacMax); From dda405ea6bda2d8f43a91bf69c8d65032072bffe Mon Sep 17 00:00:00 2001 From: Field Date: Thu, 23 Jul 2026 21:58:56 -0700 Subject: [PATCH 10/15] added the charge sharing correction polynomial coefficient classes (still need to get actual correction coefficients --- include/MModuleDepthCalibration.h | 16 +++++- src/MModuleDepthCalibration.cxx | 89 +++++++++++++++++++++++++++++-- 2 files changed, 99 insertions(+), 6 deletions(-) diff --git a/include/MModuleDepthCalibration.h b/include/MModuleDepthCalibration.h index 526dddec..91cf93a3 100644 --- a/include/MModuleDepthCalibration.h +++ b/include/MModuleDepthCalibration.h @@ -114,6 +114,10 @@ class MModuleDepthCalibration : public MModule void SetChargeSharingPolyCoeffsHV( unordered_map>> ChargeSharingPolyCoeffsHV ) { m_ChargeSharingPolyCoeffsHV = ChargeSharingPolyCoeffsHV; } void SetChargeSharingPolyCoeffsLV( unordered_map>> ChargeSharingPolyCoeffsLV ) { m_ChargeSharingPolyCoeffsLV = ChargeSharingPolyCoeffsLV; } + //! Set the coefficients of the polynomial describing the charge sharing correction, one vector per detector (map-int) per depth (the vector) + void SetChargeSharingCorrectionCoeffsHV( unordered_map>> ChargeSharingCorrectionCoeffsHV ) { m_ChargeSharingCorrectionCoeffsHV = ChargeSharingCorrectionCoeffsHV; } + void SetChargeSharingCorrectionCoeffsLV( unordered_map>> ChargeSharingCorrectionCoeffsLV ) { m_ChargeSharingCorrectionCoeffsLV = ChargeSharingCorrectionCoeffsLV; } + //! Set the depth calibration coefficients void SetCoeffs( unordered_map> Coeffs ) { m_Coeffs = Coeffs; } @@ -123,10 +127,14 @@ class MModuleDepthCalibration : public MModule //! Get the charge sharing correction calibration coefficients unordered_map>> GetChargeSharingCoeffs() { return m_ChargeSharingCoeffs; } - //! Get the coefficients of the polynomial for the ccharge sharing correction calibration + //! Get the coefficients of the polynomial for the ccharge sharing expectation / bump calibration unordered_map>> GetChargeSharingPolyCoeffsHV() { return m_ChargeSharingPolyCoeffsHV; } unordered_map>> GetChargeSharingPolyCoeffsLV() { return m_ChargeSharingPolyCoeffsLV; } + //! Get the coefficients of the polynomial for the charge sharing correction calibration + unordered_map>> GetChargeSharingCorrectionCoeffsHV() { return m_ChargeSharingCorrectionCoeffsHV; } + unordered_map>> GetChargeSharingCorrectionCoeffsLV() { return m_ChargeSharingCorrectionCoeffsLV; } + //! Get the depth calibration coefficients unordered_map> GetCoeffs() { return m_Coeffs; } @@ -182,6 +190,10 @@ class MModuleDepthCalibration : public MModule //! Determine the Grade (geometry of charge sharing) of the Hit int GetHitGrade(MHit* H); + //! Return the coefficients of the charge sharing correction polynomial for a detector at a depth + vector GetChargeSharingCorrectionCoeffsHV(int DetID,double z); + vector GetChargeSharingCorrectionCoeffsLV(int DetID,double z); + //! Return the coefficients of the dTac polynomial for a detector at a depth vector GetChargeSharingPolyCoeffsHV(int DetID,double z); vector GetChargeSharingPolyCoeffsLV(int DetID,double z); @@ -212,6 +224,8 @@ class MModuleDepthCalibration : public MModule unordered_map> m_ChargeSharingDepths; // maps DetID to the depths for which the charge sharing polynomial correction and charge sharing coefficients per-strip-pair were calculated unordered_map>> m_ChargeSharingCoeffs; // maps StripPairID to a vector of coefficients, for a vector of depths (needs interpolation) + unordered_map>> m_ChargeSharingCorrectionCoeffsLV; // maps DetID to the LV coefficients for dTAC vs CS map and charge sharing correction vs CS map, for a vector of depths + unordered_map>> m_ChargeSharingCorrectionCoeffsHV; // maps DetID to the HV coefficients for dTAC vs CS map and charge sharing correction vs CS map, for a given depth unordered_map>> m_ChargeSharingPolyCoeffsLV; // maps DetID to the LV coefficients for dTAC vs CS map and charge sharing correction vs CS map, for a vector of depths unordered_map>> m_ChargeSharingPolyCoeffsHV; // maps DetID to the HV coefficients for dTAC vs CS map and charge sharing correction vs CS map, for a given depth unordered_map> m_Coeffs; // maps pix id to a vector of coefficients... diff --git a/src/MModuleDepthCalibration.cxx b/src/MModuleDepthCalibration.cxx index c2a79d43..f57b5558 100644 --- a/src/MModuleDepthCalibration.cxx +++ b/src/MModuleDepthCalibration.cxx @@ -352,13 +352,13 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) } } - // TODO depth correction loop! - // - // step 1 -- check the HV side: - bool ChargeSharingHV = false; + // === check and correct the HV side TAC + bool ChargeSharingHV = false; // TODO flag + bool TacJitterCorrectHV = false; // TODO flag; check if problem in FM ASICS vector CorrectedHVTiming; vector CorrectedHVTimingUncertainty; - //int StripPairCode = 10000*DetID + HVStripID; // note, it is the lower strip ID always (eg, 15 if sharing between strip 15 and 16) + int StripPairCode = 10000*DetID + HVStripID + 9900; // note, it is the lower strip ID always (eg, 15 if sharing between strip 15 and 16) + // TODO actually fill the vectors and correct the TAC // -- how many strips share > 10% of the total energy (or are over slow threshold, maybe?) need this info for next steps // -- check dTAC between adjacent strips with charge sharing and also strips relative to their low-eneryg neighbors // -- correct the HV timing asic jitter bug, if needed, to make everything consistent @@ -366,6 +366,9 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) // -- if charge sharing, calculate the corrected timing value for each strip in teh absense of charge sharing, // -- and also calculate the HV tac as the weighted average, with its own uncertainty // -- note, rawZpos is used in the above calculations! + double rawCTD = (HVTiming - LVTiming); + rawCTD_s = (rawCTD - Coeffs->at(1))/(Coeffs->at(0)); //apply inverse stretch and offset + std::tie(rawZpos, rawZsigma) = CalculateZfromCTD(rawCTD_s, noise,DetID, Grade, false); // true (sean weighting) // step 2 -- check the LV side: bool ValidatedLVTiming = false; // should put a flag here eventually TODO; if NN do not have fast timing @@ -750,6 +753,14 @@ bool MModuleDepthCalibration::LoadChargeSharingConfigFile(MString FileName) poly_coeffs.push_back(139.209); poly_coeffs.push_back(-106.849);// these are the coefficients, where we'll have a polynomial dTac = coeffs[0]*(f-0.5) - coeffs[1]*(f-0.5)^3 m_ChargeSharingPolyCoeffsHV[DetID].push_back(poly_coeffs); m_ChargeSharingPolyCoeffsLV[DetID].push_back(poly_coeffs); // in principle they would be different + + vector correction_coeffs; + correction_coeffs.push_back(0); // TODO get actual coeffs from Isidro!!! + correction_coeffs.push_back(0); + correction_coeffs.push_back(0); + correction_coeffs.push_back(0); + m_ChargeSharingCorrectionCoeffsLV[DetID].push_back(correction_coeffs); + m_ChargeSharingCorrectionCoeffsHV[DetID].push_back(correction_coeffs); // fill m_ChargeSharingConfig for each detector / depth / strip pair vector coeffs; // the stretch and offset, currently set to the same values for all strips (which are almost certainly wrong) @@ -937,6 +948,74 @@ std::vector MModuleDepthCalibration::GetChargeSharingPolyCoeffsHV(int De return {}; } +} +///////////////////////////////////////////////////////////////////////////////// + + +std::vector MModuleDepthCalibration::GetChargeSharingCorrectionCoeffsLV(int DetID, double z) +{ + // Check to see if the charge sharing coefficients have been loaded. If so, try to get the coefficients for the specified strip pair. + if (m_ChargeSharingConfigFileIsLoaded == true) { + if (m_ChargeSharingCorrectionCoeffsLV.count(DetID) > 0) { + + // if we only sampled one depth, or we're beyond the depth range, just return the closest coefficients + if (z <= m_ChargeSharingDepths[DetID].front()) return m_ChargeSharingCorrectionCoeffsLV[DetID].at(0); + if (z >= m_ChargeSharingDepths[DetID].back()) return m_ChargeSharingCorrectionCoeffsLV[DetID].at(m_ChargeSharingDepths[DetID].size()-1); + + // otherwise, interpolate + for (unsigned int i = 0; i < m_ChargeSharingDepths[DetID].size() - 1; i++){ + if (z >= m_ChargeSharingDepths[DetID].at(i) && z < m_ChargeSharingDepths[DetID].at(i + 1)) { + double f = (z - m_ChargeSharingDepths[DetID][i]) / (m_ChargeSharingDepths[DetID][i + 1] - m_ChargeSharingDepths[DetID][i]); + vector result; + for (int j = 0; j < m_ChargeSharingCorrectionCoeffsLV[DetID].at(i).size(); j++) result.push_back((1.0 - f) * m_ChargeSharingCorrectionCoeffsLV[DetID][i][j] + f * m_ChargeSharingCorrectionCoeffsLV[DetID][i + 1][j]); + return result; + } + } + } else { + if (g_Verbosity >= c_Warning) { + cout << "MModuleDepthCalibration::GetChargeSharingCorrectionCoeffsLV: cannot get charge sharing correction coefficients; detector id code " << DetID << " not found." << endl; + } + return {}; + } + } else { + cout << "MModuleDepthCalibration::GetChargeSharingCorrectionCoeffsLV: cannot get charge sharing correction coefficients; file has not yet been loaded." << endl; + return {}; + } + +} +///////////////////////////////////////////////////////////////////////////////// + + +std::vector MModuleDepthCalibration::GetChargeSharingCorrectionCoeffsHV(int DetID, double z) +{ + // Check to see if the charge sharing coefficients have been loaded. If so, try to get the coefficients for the specified strip pair. + if (m_ChargeSharingConfigFileIsLoaded == true) { + if (m_ChargeSharingCorrectionCoeffsHV.count(DetID) > 0) { + + // if we only sampled one depth, or we're beyond the depth range, just return the closest coefficients + if (z <= m_ChargeSharingDepths[DetID].front()) return m_ChargeSharingCorrectionCoeffsHV[DetID].at(0); + if (z >= m_ChargeSharingDepths[DetID].back()) return m_ChargeSharingCorrectionCoeffsHV[DetID].at(m_ChargeSharingDepths[DetID].size()-1); + + // otherwise, interpolate + for (unsigned int i = 0; i < m_ChargeSharingDepths[DetID].size() - 1; i++){ + if (z >= m_ChargeSharingDepths[DetID].at(i) && z < m_ChargeSharingDepths[DetID].at(i + 1)) { + double f = (z - m_ChargeSharingDepths[DetID][i]) / (m_ChargeSharingDepths[DetID][i + 1] - m_ChargeSharingDepths[DetID][i]); + vector result; + for (int j = 0; j < m_ChargeSharingCorrectionCoeffsHV[DetID].at(i).size(); j++) result.push_back((1.0 - f) * m_ChargeSharingCorrectionCoeffsHV[DetID][i][j] + f * m_ChargeSharingCorrectionCoeffsHV[DetID][i + 1][j]); + return result; + } + } + } else { + if (g_Verbosity >= c_Warning) { + cout << "MModuleDepthCalibration::GetChargeSharingCorrectionCoeffsHV: cannot get charge sharing correction coefficients; detector id code " << DetID << " not found." << endl; + } + return {}; + } + } else { + cout << "MModuleDepthCalibration::GetChargeSharingCorrectionCoeffsHV: cannot get charge sharing coefficients; file has not yet been loaded." << endl; + return {}; + } + } ///////////////////////////////////////////////////////////////////////////////// From 149b39a71270bcbf2b0928eb6f61249daf117786 Mon Sep 17 00:00:00 2001 From: Field Date: Thu, 23 Jul 2026 22:39:04 -0700 Subject: [PATCH 11/15] added functionality for charge sharing correction coefficients --- src/MModuleDepthCalibration.cxx | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/src/MModuleDepthCalibration.cxx b/src/MModuleDepthCalibration.cxx index f57b5558..f5473015 100644 --- a/src/MModuleDepthCalibration.cxx +++ b/src/MModuleDepthCalibration.cxx @@ -177,6 +177,7 @@ void MModuleDepthCalibration::CreateExpos() 100, 0, 1); // fraction bins } m_Expos.push_back(m_ExpoPlotTacDiff); + cout << "added dTAC expo"< 4) { // GRADE=5 is some complicated geometry with multiple hits on a single strip. GRADE=6 means not all strips are adjacent. H->SetNoDepth(); Event->SetDepthCalibrationError("Multiple hits on single strip"); @@ -236,6 +238,7 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) } else if (Grade==6) { ++m_Error6; } + continue; } else { // If the Grade is 0-4, we can handle it. // Calculate the position. If error is thrown, record and no depth. @@ -253,7 +256,6 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) MStripHit* LVSH = GetDominantStrip(LVStrips, LVEnergyFraction); MStripHit* HVSH = GetDominantStrip(HVStrips, HVEnergyFraction); - double rawCTD_s = 0.0; //now try and get z position int DetID = LVSH->GetDetectorID(); @@ -304,25 +306,29 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) H->SetNoDepth(); Event->SetDepthCalibrationError("No calibration coefficients"); ++m_Error1; + continue; } else if (CTDVec.size() == 0) { if (g_Verbosity >= c_Error) cout << m_XmlTag << "Empty CTD vector" << endl; H->SetNoDepth(); Event->SetDepthCalibrationError("No calibration coefficients"); + continue; } else if (DepthVec.size() == 0) { if (g_Verbosity >= c_Error) cout << m_XmlTag << "Empty Depth vector" << endl; H->SetNoDepth(); Event->SetDepthCalibrationError("No calibration coefficients"); + continue; } else if ((LVTiming < 1.0E-6) || (HVTiming < 1.0E-6)) { ++m_Error3; H->SetNoDepth(); Event->SetDepthCalibrationError("No timing"); + continue; } else { // If there are coefficients and timing information is loaded, try calculating the CTD and depth // TODO FR start here double rawCTD = (HVTiming - LVTiming); - rawCTD_s = (rawCTD - Coeffs->at(1))/(Coeffs->at(0)); //apply inverse stretch and offset + double rawCTD_s = (rawCTD - Coeffs->at(1))/(Coeffs->at(0)); //apply inverse stretch and offset double Xmin = * std::min_element(CTDVec.begin(), CTDVec.end()); double Xmax = * std::max_element(CTDVec.begin(), CTDVec.end()); @@ -355,9 +361,9 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) // === check and correct the HV side TAC bool ChargeSharingHV = false; // TODO flag bool TacJitterCorrectHV = false; // TODO flag; check if problem in FM ASICS - vector CorrectedHVTiming; - vector CorrectedHVTimingUncertainty; - int StripPairCode = 10000*DetID + HVStripID + 9900; // note, it is the lower strip ID always (eg, 15 if sharing between strip 15 and 16) + //vector CorrectedHVTiming; + //vector CorrectedHVTimingUncertainty; + //int StripPairCode = 10000*DetID + HVStripID + 9900; // note, it is the lower strip ID always (eg, 15 if sharing between strip 15 and 16) // TODO actually fill the vectors and correct the TAC // -- how many strips share > 10% of the total energy (or are over slow threshold, maybe?) need this info for next steps // -- check dTAC between adjacent strips with charge sharing and also strips relative to their low-eneryg neighbors @@ -366,9 +372,9 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) // -- if charge sharing, calculate the corrected timing value for each strip in teh absense of charge sharing, // -- and also calculate the HV tac as the weighted average, with its own uncertainty // -- note, rawZpos is used in the above calculations! - double rawCTD = (HVTiming - LVTiming); - rawCTD_s = (rawCTD - Coeffs->at(1))/(Coeffs->at(0)); //apply inverse stretch and offset - std::tie(rawZpos, rawZsigma) = CalculateZfromCTD(rawCTD_s, noise,DetID, Grade, false); // true (sean weighting) + //rawCTD = (HVTiming - LVTiming); + //rawCTD_s = (rawCTD - Coeffs->at(1))/(Coeffs->at(0)); //apply inverse stretch and offset + //auto [HVZpos, HVZsigma] = CalculateZfromCTD(rawCTD_s, noise,DetID, Grade, false); // true (sean weighting) // step 2 -- check the LV side: bool ValidatedLVTiming = false; // should put a flag here eventually TODO; if NN do not have fast timing From 22a835faffb03c1e33262a8bba4b76d50a5c2f00 Mon Sep 17 00:00:00 2001 From: Field Date: Thu, 23 Jul 2026 22:39:23 -0700 Subject: [PATCH 12/15] fixed a very stupid seg fault in src/MGUIExpoPlotSpectrum.cxx --- src/MGUIExpoPlotSpectrum.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/MGUIExpoPlotSpectrum.cxx b/src/MGUIExpoPlotSpectrum.cxx index 35ba3acc..2210f598 100644 --- a/src/MGUIExpoPlotSpectrum.cxx +++ b/src/MGUIExpoPlotSpectrum.cxx @@ -62,7 +62,7 @@ MGUIExpoPlotSpectrum::MGUIExpoPlotSpectrum(MModule* Module) : MGUIExpo(Module) // Set the new title of the tab here: if (Module != nullptr) { - m_TabTitle = "Energy Spectrum (" + Module->GetName() + ")"; + m_TabTitle = "ESpec (" + Module->GetName() + ")"; } else { m_TabTitle = "Energy Spectrum"; } From 3c988e98abd5c210b53a82f2f3a791de251f7820 Mon Sep 17 00:00:00 2001 From: Field Date: Fri, 24 Jul 2026 12:50:39 -0700 Subject: [PATCH 13/15] added all charge sharing back in for now --- src/MModuleDepthCalibration.cxx | 104 +++++++++++++++++--------------- 1 file changed, 54 insertions(+), 50 deletions(-) diff --git a/src/MModuleDepthCalibration.cxx b/src/MModuleDepthCalibration.cxx index f5473015..44b90f0f 100644 --- a/src/MModuleDepthCalibration.cxx +++ b/src/MModuleDepthCalibration.cxx @@ -173,7 +173,7 @@ void MModuleDepthCalibration::CreateExpos() double thickness = m_Thicknesses[DetID]; m_ExpoPlotTacDiff->SetHistogramParameters(DetID, 120, -thickness/2.0, thickness/2.0, // depth bins - 200, -100, 100, // dTac bins + 100, -150, 150, // dTac bins 100, 0, 1); // fraction bins } m_Expos.push_back(m_ExpoPlotTacDiff); @@ -361,9 +361,9 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) // === check and correct the HV side TAC bool ChargeSharingHV = false; // TODO flag bool TacJitterCorrectHV = false; // TODO flag; check if problem in FM ASICS - //vector CorrectedHVTiming; - //vector CorrectedHVTimingUncertainty; - //int StripPairCode = 10000*DetID + HVStripID + 9900; // note, it is the lower strip ID always (eg, 15 if sharing between strip 15 and 16) + vector CorrectedHVTiming; + vector CorrectedHVTimingUncertainty; + int StripPairCode = 10000*DetID + HVStripID + 9900; // note, it is the lower strip ID always (eg, 15 if sharing between strip 15 and 16) // TODO actually fill the vectors and correct the TAC // -- how many strips share > 10% of the total energy (or are over slow threshold, maybe?) need this info for next steps // -- check dTAC between adjacent strips with charge sharing and also strips relative to their low-eneryg neighbors @@ -372,64 +372,68 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) // -- if charge sharing, calculate the corrected timing value for each strip in teh absense of charge sharing, // -- and also calculate the HV tac as the weighted average, with its own uncertainty // -- note, rawZpos is used in the above calculations! - //rawCTD = (HVTiming - LVTiming); - //rawCTD_s = (rawCTD - Coeffs->at(1))/(Coeffs->at(0)); //apply inverse stretch and offset - //auto [HVZpos, HVZsigma] = CalculateZfromCTD(rawCTD_s, noise,DetID, Grade, false); // true (sean weighting) + rawCTD = (HVTiming - LVTiming); + rawCTD_s = (rawCTD - Coeffs->at(1))/(Coeffs->at(0)); //apply inverse stretch and offset + auto [HVZpos, HVZsigma] = CalculateZfromCTD(rawCTD_s, noise,DetID, Grade, false); // true (sean weighting) - // step 2 -- check the LV side: + + // == check and correct the LV side: + // -------------------------------- bool ValidatedLVTiming = false; // should put a flag here eventually TODO; if NN do not have fast timing bool ZombieBump = false;// should put a flag here TODO vector CorrectedLVTiming; vector CorrectedLVTimingUncertainty; - bool ChargeSharingLV = false;// also should put a flag here TODO (is there another flag for charge sharing?) - - if (LVEnergyFraction > m_SingleStripChargeSharing){// if we have one obvious main strip, we are not going to be correcting charge sharing but just checking for zombie bump - // compare with the neighbors, if possible - for (int neighbor = 0; neighbor < 2; neighbor++){// 0 for left neighbor, 1 for right - int pm = 2*neighbor - 1; // -1 for neighbor 0 (neighbor is left); + 1 for neighbor == 1 (right neighbor, which is nominal for the convention StripPairID = left StripID of pair - int NeighborStripID = LVStripID + pm; - MStripHit* NSH = GetStrip(LVStrips, NeighborStripID); - if (NeighborStripID >=0 && NeighborStripID <=63 && NSH && NSH->HasFastTiming()){ // neighbor is not a guard ring strip, NSH exists (not a null pointer) and has fast timing! - double dTacData = (LVSH->GetTiming() - NSH->GetTiming())*pm; // always the left strip - right strip; neighbor on left means pm = -1 -> NSH - LVSH timing - double fracData = (NSH->GetEnergy()/(LVSH->GetEnergy() + NSH->GetEnergy())*pm) + 1 - neighbor; // always the fraction on the right stripHit; nominally NSH for right neighbor - int StripPairCode = 10000*DetID + 100*(LVStripID - 1 + neighbor) + 99; //LVStripID -1 + 0 = LVStripID -1 (left neighbor); or LVStripID -1 + 1 = LVStripID (LVStrip is the StripID when we consider right negihbor) - double x = (fracData - 0.5); // x is LVEnergyFraction - 0.5, ie the parameter of Isidro's polynomials, which are forced to go through (0.5, 0) - vector CSPolyCoeffs = GetChargeSharingPolyCoeffsLV(DetID,rawZpos); // TODO need to actually check and fill the variable that checks the length of this, and check that it's right when loading - vector CSCoeffs = GetChargeSharingCoeffs(StripPairCode,rawZpos); - - // add to the expo - if (HasExpos() == true) { - m_ExpoPlotTacDiff->AddData(StripPairCode, rawZpos, dTacData, fracData, std::nan(""),HitEnergy); - } - - if (!CSCoeffs.empty() && !CSPolyCoeffs.empty()){ - double dTacExpect = (CSPolyCoeffs.at(0)*x + CSPolyCoeffs.at(1)*x*x*x)*CSCoeffs.at(0) + CSCoeffs.at(1);// TODO update if not cubic polynomial - // TODO we should display (dTacData - dTacExpect)*pm to keep an eye on the prevalence of the bump - if (rawZpos > -0.5 && (dTacData - dTacExpect)*pm > 2*noise){// zombie bump! TODO update bump criteria, fix noise - ZombieBump = true; - ValidatedLVTiming = false; - cout << "LV Strip: "< -1*pm*dTacExpect) { // zombie bump with good neighbor. Can we not always do this in this case, though, since we know what the timing should be? - CorrectedLVTiming.push_back(LVTiming-60); // needs to be a config - CorrectedLVTimingUncertainty.push_back(noise*2); // to do: quantify and make into something real - } else { - CorrectedLVTiming.push_back(0); // needs to be a config - CorrectedLVTimingUncertainty.push_back(0); // to do: quantify and make into something real - } + bool ChargeSharingLV = LVEnergyFraction > m_SingleStripChargeSharing;// also should put a flag here TODO (is there another flag for charge sharing?) + double MaxEnergy = LVSH->GetEnergy(); + + // compare with the neighbors -- calculate parameters from data + for (int neighbor = 0; neighbor < 2; neighbor++){// 0 for left neighbor, 1 for right + int pm = 2*neighbor - 1; // -1 for neighbor 0 (neighbor is left); + 1 for neighbor == 1 (right neighbor, which is nominal for the convention StripPairID = left StripID of pair + int NeighborStripID = LVStripID + pm; + MStripHit* NSH = GetStrip(LVStrips, NeighborStripID); + if (NSH && NSH->HasFastTiming()){ // NSH exists (not a null pointer) and has fast timing! + double dTacData = (LVSH->GetTiming() - NSH->GetTiming())*pm; // always the left strip - right strip; neighbor on left means pm = -1 -> NSH - LVSH timing + double fracData = (NSH->GetEnergy()/(LVSH->GetEnergy() + NSH->GetEnergy())*pm) + 1 - neighbor; // always the fraction on the right stripHit; nominally NSH for right neighbor + int StripPairCode = 10000*DetID + 100*(LVStripID - 1 + neighbor) + 99; //LVStripID -1 + 0 = LVStripID -1 (left neighbor); or LVStripID -1 + 1 = LVStripID (LVStrip is the StripID when we consider right negihbor) + if (HasExpos() == true) m_ExpoPlotTacDiff->AddData(StripPairCode, rawZpos, dTacData, fracData, std::nan(""),HitEnergy); // TODO add logic to do this later so we can get dTac dTac + + // get the expectation + vector CSPolyCoeffs = GetChargeSharingPolyCoeffsLV(DetID,rawZpos); // TODO need to actually check and fill the variable that checks the length of this, and check that it's right when loading + vector CSCoeffs = GetChargeSharingCoeffs(StripPairCode,rawZpos); + if (!CSCoeffs.empty() && !CSPolyCoeffs.empty()){ + double x = (fracData - 0.5); // x is LVEnergyFraction - 0.5, ie the parameter of Isidro's polynomials, which are forced to go through (0.5, 0) + double dTacExpect = (CSPolyCoeffs.at(0)*x + CSPolyCoeffs.at(1)*x*x*x)*CSCoeffs.at(0) + CSCoeffs.at(1);// TODO update if not cubic polynomial + + // check for zombie bump + if (rawZpos > -0.5 && (dTacData - dTacExpect)*pm > 25){// zombie bump! TODO update bump criteria, fix noise (should be noise, not 25, but the noise is all wrong + ZombieBump = true; + ValidatedLVTiming = false; + cout << "Zombie Bump!!! LV Strip: "< 80) { // zombie bump with good neighbor. Can we not always do this in this case, though, since we know what the timing should be? TODO 100 should be calibrated + CorrectedLVTiming.push_back(LVTiming-150); // TODO needs to be a config + CorrectedLVTimingUncertainty.push_back(noise*2); // TODO: quantify and make into something real } else { - if (!ZombieBump && !CSCoeffs.empty() && !CSPolyCoeffs.empty()) ValidatedLVTiming = true; - } - } - } + CorrectedLVTiming.push_back(0); // we will not use it in the calculation later + CorrectedLVTimingUncertainty.push_back(0); + } + } else { + if (!ZombieBump && !CSCoeffs.empty() && !CSPolyCoeffs.empty()) ValidatedLVTiming = true; + } + + // now do the charge sharing correction + // TODO charge sharing correction + // TODO push bakc to corrected LV timing goes here!!! + } } - } else { // charge sharing correction - ChargeSharingLV = true; - } + // TODO loop to check NNN if high charge sharing on respective neighbor + } + // correct the timing if (CorrectedLVTiming.size() > 0){// if we have a correction. Also, implemented weighting! TODO an double correctionSum = 0; for (unsigned int j = 0; j < CorrectedLVTiming.size(); j++) { // TODO check that they are consistent and drop one if not.... + // TODO figure out actual logic with the charge sharing correction included, in combination with zombie bump. I think this involves calculating both the expectation for neighbor and main strip in the case of charge sharing and zombie bump (?) correctionSum += CorrectedLVTiming.at(j); } LVTiming = correctionSum / CorrectedLVTiming.size(); From 57ba33e48b306645123b1d568d1318623c3138da Mon Sep 17 00:00:00 2001 From: Field Date: Fri, 24 Jul 2026 12:50:58 -0700 Subject: [PATCH 14/15] fixed plot axis labels on dTAC display --- src/MGUIExpoPlotTacDiff.cxx | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/MGUIExpoPlotTacDiff.cxx b/src/MGUIExpoPlotTacDiff.cxx index 32b5a46e..00b7e64b 100644 --- a/src/MGUIExpoPlotTacDiff.cxx +++ b/src/MGUIExpoPlotTacDiff.cxx @@ -84,12 +84,12 @@ void MGUIExpoPlotTacDiff::SetHistogramParameters(unsigned int DetID, unsigned in m_DtacVsDepthHistograms[key] = new TH2D("", TString::Format("Det %d LV %d: dTAC vs Depth", DetID, s), NBinsDepth, DepthMin, DepthMax, NBinsDtac, DtacMin, DtacMax); m_DtacVsDepthHistograms[key]->SetXTitle("Depth [cm]"); - m_DtacVsDepthHistograms[key]->SetYTitle(TString::Format("TAC %s - TAC %s [ns]",s,s+1)); + m_DtacVsDepthHistograms[key]->SetYTitle(TString::Format("TAC %d - TAC %d [ns]",s,s+1)); m_DtacVsFracHistograms[key] = new TH2D("", TString::Format("Det %d Charge Shared between LV %d and %d", DetID, s,s+1), NBinsFrac, FracMin, FracMax, NBinsDtac, DtacMin, DtacMax); - m_DtacVsFracHistograms[key]->SetXTitle(TString::Format("Charge Sharing Fraction: (strip %d) / (strip %d + strip %d",s,s,s+1)); - m_DtacVsFracHistograms[key]->SetYTitle(TString::Format("TAC %s - TAC %s [ns]",s,s+1)); + m_DtacVsFracHistograms[key]->SetXTitle(TString::Format("Charge Sharing Fraction: (strip %d) / (strip %d + strip %d)",s,s,s+1)); + m_DtacVsFracHistograms[key]->SetYTitle(TString::Format("TAC %d - TAC %d [ns]",s,s+1)); m_DtacVsDtacHistograms[key] = new TH2D("", TString::Format("Det %d LV %d: dTAC vs dTAC", DetID, s), NBinsDtac, DtacMin, DtacMax, NBinsDtac, DtacMin, DtacMax); @@ -105,12 +105,12 @@ void MGUIExpoPlotTacDiff::SetHistogramParameters(unsigned int DetID, unsigned in m_DtacVsDepthHistograms[key] = new TH2D("", TString::Format("Det %d HV %d: dTAC vs Depth", DetID, s), NBinsDepth, DepthMin, DepthMax, NBinsDtac, DtacMin, DtacMax); m_DtacVsDepthHistograms[key]->SetXTitle("Depth [cm]"); - m_DtacVsDepthHistograms[key]->SetYTitle(TString::Format("TAC %s - TAC %s [ns]",s,s+1)); + m_DtacVsDepthHistograms[key]->SetYTitle(TString::Format("TAC %d - TAC %d [ns]",s,s+1)); - m_DtacVsFracHistograms[key] = new TH2D("", TString::Format("Det %d Charge Shared between LV %s and %s", DetID, s,s+1), + m_DtacVsFracHistograms[key] = new TH2D("", TString::Format("Det %d Charge Shared between LV %d and %d", DetID, s,s+1), NBinsFrac, FracMin, FracMax, NBinsDtac, DtacMin, DtacMax); - m_DtacVsFracHistograms[key]->SetXTitle(TString::Format("Charge Sharing Fraction: (strip %d) / (strip %d + strip %d",s,s,s+1)); - m_DtacVsFracHistograms[key]->SetYTitle(TString::Format("TAC %s - TAC %s [ns]",s,s+1)); + m_DtacVsFracHistograms[key]->SetXTitle(TString::Format("Charge Sharing Fraction: (strip %d) / (strip %d + strip %d)",s,s,s+1)); + m_DtacVsFracHistograms[key]->SetYTitle(TString::Format("TAC %d - TAC %d [ns]",s,s+1)); m_DtacVsDtacHistograms[key] = new TH2D("", TString::Format("Det %d HV %d: dTAC vs dTAC", DetID, s), NBinsDtac, DtacMin, DtacMax, NBinsDtac, DtacMin, DtacMax); @@ -213,7 +213,7 @@ void MGUIExpoPlotTacDiff::Create() TGNumberFormat::kNEANonNegative, TGNumberFormat::kNELLimitMinMax, 0, 62); //m_StripEntry->Connect("ValueSet(Long_t)", "MGUIExpoPlotTacDiff", this, "OnStripSelected()"); - m_StripEntry->GetNumberEntry()->Connect("ReturnPressed()", "MGUIExpoPlotTacDiff", this, "OnStripSelected()"); + //m_StripEntry->GetNumberEntry()->Connect("ReturnPressed()", "MGUIExpoPlotTacDiff", this, "OnStripSelected()"); m_StripEntry->Resize(50, 20); StripRow->AddFrame(m_StripEntry, EntryLayout); @@ -226,7 +226,7 @@ void MGUIExpoPlotTacDiff::Create() TGNumberFormat::kNESRealTwo, TGNumberFormat::kNEANonNegative); //m_EnergyMinEntry->Connect("ValueSet(Long_t)", "MGUIExpoPlotTacDiff", this, "OnEnergyRangeChanged()"); - m_EnergyMinEntry->GetNumberEntry()->Connect("ReturnPressed()", "MGUIExpoPlotTacDiff", this, "OnEnergyRangeChanged()"); + //m_EnergyMinEntry->GetNumberEntry()->Connect("ReturnPressed()", "MGUIExpoPlotTacDiff", this, "OnEnergyRangeChanged()"); m_EnergyMinEntry->Resize(70, 20); EMinRow->AddFrame(m_EnergyMinEntry, EntryLayout); @@ -239,7 +239,7 @@ void MGUIExpoPlotTacDiff::Create() TGNumberFormat::kNESRealTwo, TGNumberFormat::kNEANonNegative); //m_EnergyMaxEntry->Connect("ValueSet(Long_t)", "MGUIExpoPlotTacDiff", this, "OnEnergyRangeChanged()"); - m_EnergyMaxEntry->GetNumberEntry()->Connect("ReturnPressed()", "MGUIExpoPlotTacDiff", this, "OnEnergyRangeChanged()"); + //m_EnergyMaxEntry->GetNumberEntry()->Connect("ReturnPressed()", "MGUIExpoPlotTacDiff", this, "OnEnergyRangeChanged()"); m_EnergyMaxEntry->Resize(70, 20); EMaxRow->AddFrame(m_EnergyMaxEntry, EntryLayout); From 5a7b946b20ce24630d4720805d0c74b3681dac76 Mon Sep 17 00:00:00 2001 From: Field Date: Fri, 24 Jul 2026 12:53:16 -0700 Subject: [PATCH 15/15] suppress verbosity from debugging src/MGUIExpoDepthCalibration.cxx --- src/MGUIExpoDepthCalibration.cxx | 1 - 1 file changed, 1 deletion(-) diff --git a/src/MGUIExpoDepthCalibration.cxx b/src/MGUIExpoDepthCalibration.cxx index cfac4742..755efec5 100644 --- a/src/MGUIExpoDepthCalibration.cxx +++ b/src/MGUIExpoDepthCalibration.cxx @@ -281,7 +281,6 @@ void MGUIExpoDepthCalibration::RebuildDisplayHistograms() int key = GetStripKey(DetID, m_SelectedSide, s); if (m_DepthPerStrip.count(key) > 0) m_DepthHistograms[DetID]->Add(m_DepthPerStrip[key]); - cout << " added "<< key << ": " << m_DepthPerStrip[key]->GetEntries() << " entries" << endl; if (m_RawDepthPerStrip.count(key) > 0) m_RawDepthHistograms[DetID]->Add(m_RawDepthPerStrip[key]); }