From 6d926ee4576d2adc36dd7eb4ee3199ec475345c2 Mon Sep 17 00:00:00 2001 From: ckierans Date: Tue, 5 Aug 2025 10:13:21 -0400 Subject: [PATCH 01/22] ADD: Including Mask Metrology file in Depth Calibration, correction still needs work --- include/MGUIOptionsDepthCalibration2024.h | 7 + include/MModuleDepthCalibration2024.h | 42 +++++- src/MGUIOptionsDepthCalibration2024.cxx | 54 +++++-- src/MModuleDepthCalibration2024.cxx | 167 +++++++++++++++++++++- 4 files changed, 254 insertions(+), 16 deletions(-) diff --git a/include/MGUIOptionsDepthCalibration2024.h b/include/MGUIOptionsDepthCalibration2024.h index 15c72ed6..6c14df24 100644 --- a/include/MGUIOptionsDepthCalibration2024.h +++ b/include/MGUIOptionsDepthCalibration2024.h @@ -74,6 +74,13 @@ class MGUIOptionsDepthCalibration2024 : public MGUIOptions //! Select spline file to load, splines will convert CTD->Depth MGUIEFileSelector* m_SplinesFileSelector; + //! Select mask metrology file to load. This gives the translation and rotation for each strip in the detector frame + MGUIEFileSelector* m_MaskMetrologyFileSelector; + + int m_UseMaskMetCorr; + TGCheckButton* m_MaskMetModeCB; + enum ButtonIDs {c_MetFile}; + //! Check button if working with the Card Cage at UCSD TGCheckButton* m_UCSDOverride; diff --git a/include/MModuleDepthCalibration2024.h b/include/MModuleDepthCalibration2024.h index 21f48291..d0cb51c5 100644 --- a/include/MModuleDepthCalibration2024.h +++ b/include/MModuleDepthCalibration2024.h @@ -1,4 +1,4 @@ -/* +/*:q * MModuleDepthCalibration2024.h * * Copyright (C) 2008-2008 by Andreas Zoglauer. @@ -72,6 +72,16 @@ class MModuleDepthCalibration2024 : public MModule //! Get filename for CTD->Depth splines MString GetSplinesFileName() const {return m_SplinesFile;} + //! Enable/Disable Preamp Temp Correction + void EnableMaskMetrologyCorrection(bool X) {m_MaskMetrologyEnabled = X;} + //! Get coincidence merging true/false + bool GetMaskMetrologyCorrection() const { return m_MaskMetrologyEnabled; } + + //! Set filename for mask metrology + void SetMaskMetrologyFileName( const MString& FileName) {m_MaskMetrologyFile = FileName;} + //! Get filename for CTD->Depth splines + MString GetMaskMetrologyFileName() const {return m_MaskMetrologyFile;} + //! Set whether the data came from the card cage at UCSD void SetUCSDOverride( bool Override ) {m_UCSDOverride = Override;} //! Get whether the data came from the card cage at UCSD @@ -91,22 +101,40 @@ class MModuleDepthCalibration2024 : public MModule protected: //! Returns the strip with most energy from vector Strips, also gives back the energy fraction MStripHit* GetDominantStrip(std::vector& Strips, double& EnergyFraction); + //! Retrieve the appropriate Depth values given the DetID - vector GetDepth(int DetID); + vector GetDepth(int DetID); + //! Retrieve the appropriate CTD values given the DetID and Grade vector GetCTD(int DetID, int Grade); + //! Normal distribution vector norm_pdf(vector x, double mu, double sigma); - //! Adds a Depth-to-CTD relation - bool AddDepthCTD(vector depthvec, vector> ctdarr, int DetID, unordered_map>& DepthGrid, unordered_map>>& CTDMap); + + //! Adds a Depth-to-CTD relation + bool AddDepthCTD(vector depthvec, vector> ctdarr, int DetID, unordered_map>& DepthGrid, unordered_map>>& CTDMap); + //! Determine the Grade (geometry of charge sharing) of the Hit int GetHitGrade(MHit* H); + //! Load in the specified coefficients file bool LoadCoeffsFile(MString FName); + //! Return the coefficients for a pixel vector* GetPixelCoeffs(int PixelCode); + //! Load the splines file bool LoadSplinesFile(MString FName); + + //! Mask Metrology Correction + bool m_MaskMetrologyEnabled; + + //! Load the metrology mask file + bool LoadMaskMetrologyFile(MString FName); + + //! Get the x, y position of the intersection of two strips based on the Metrology Mask + vector GetStripIntersection(MReadOutElementDoubleStrip LVStrip, MReadOutElementDoubleStrip HVStrip); + //! Get the timing FWHM noise for the specified pixel and Energy double GetTimingNoiseFWHM(int PixelCode, double Energy); @@ -147,6 +175,12 @@ class MModuleDepthCalibration2024 : public MModule unordered_map> m_DepthGrid; bool m_SplinesFileIsLoaded; bool m_CoeffsFileIsLoaded; + bool m_MaskMetrologyFileIsLoaded; + + // The Mask Metrology + MString m_MaskMetrologyFile; + map> m_MaskMetrology; + // boolean for use with the card cage at UCSD since it tags all events as detector 11 bool m_UCSDOverride; diff --git a/src/MGUIOptionsDepthCalibration2024.cxx b/src/MGUIOptionsDepthCalibration2024.cxx index cd6c07ec..9ead8652 100644 --- a/src/MGUIOptionsDepthCalibration2024.cxx +++ b/src/MGUIOptionsDepthCalibration2024.cxx @@ -75,8 +75,30 @@ void MGUIOptionsDepthCalibration2024::Create() m_SplinesFileSelector = new MGUIEFileSelector(m_OptionsFrame, "Select a splines file:", dynamic_cast(m_Module)->GetSplinesFileName()); m_SplinesFileSelector->SetFileType("splines", "*.csv"); - TGLayoutHints* Label2Layout = new TGLayoutHints(kLHintsTop | kLHintsCenterX | kLHintsExpandX, 10, 10, 10, 10); - m_OptionsFrame->AddFrame(m_SplinesFileSelector, Label2Layout); +// TGLayoutHints* Label2Layout = new TGLayoutHints(kLHintsTop | kLHintsCenterX | kLHintsExpandX, 10, 10, 10, 10); + m_OptionsFrame->AddFrame(m_SplinesFileSelector, LabelLayout); + + m_MaskMetModeCB = new TGCheckButton(m_OptionsFrame, "Enable mask metrology correction and read calibration from file:", c_MetFile); + m_MaskMetModeCB->SetState((dynamic_cast(m_Module)->GetMaskMetrologyCorrection() == 1) ? kButtonDown : kButtonUp); + m_MaskMetModeCB->Associate(this); + m_OptionsFrame->AddFrame(m_MaskMetModeCB, LabelLayout); + + m_UseMaskMetCorr = dynamic_cast(m_Module)->GetMaskMetrologyCorrection(); + +// TGLayoutHints* FileLabelLayout = new TGLayoutHints(kLHintsTop | kLHintsExpandX, m_FontScaler*65 + 21*m_FontScaler, m_FontScaler*65, 0, 2*m_FontScaler); + + m_MaskMetrologyFileSelector = new MGUIEFileSelector(m_OptionsFrame, "", dynamic_cast(m_Module)->GetMaskMetrologyFileName()); + m_MaskMetrologyFileSelector->SetFileType("metrology", "*.csv"); + m_OptionsFrame->AddFrame(m_MaskMetrologyFileSelector, LabelLayout); + + + if (m_UseMaskMetCorr) { + m_MaskMetrologyFileSelector->SetEnabled(true); + } else { + m_MaskMetrologyFileSelector->SetEnabled(false); + } + + m_UCSDOverride = new TGCheckButton(m_OptionsFrame, "Check this box if you're using the card cage at UCSD", 1); m_UCSDOverride->SetOn(dynamic_cast(m_Module)->GetUCSDOverride()); @@ -97,16 +119,28 @@ bool MGUIOptionsDepthCalibration2024::ProcessMessage(long Message, long Paramete bool Status = true; switch (GET_MSG(Message)) { - case kC_COMMAND: - switch (GET_SUBMSG(Message)) { - case kCM_BUTTON: - break; - default: + case kC_COMMAND: + switch (GET_SUBMSG(Message)) { + case kCM_BUTTON: + break; + case kCM_CHECKBUTTON: + switch (Parameter1) { + case c_MetFile: + if (m_MaskMetModeCB->GetState() == kButtonDown) { + m_UseMaskMetCorr = 1; + m_MaskMetrologyFileSelector->SetEnabled(true); + } else if (m_MaskMetModeCB->GetState() == kButtonUp) { + m_UseMaskMetCorr = 0; + m_MaskMetrologyFileSelector->SetEnabled(false); + } break; } - break; default: break; + } + break; + default: + break; } if (Status == false) { @@ -127,6 +161,10 @@ bool MGUIOptionsDepthCalibration2024::OnApply() dynamic_cast(m_Module)->SetCoeffsFileName(m_CoeffsFileSelector->GetFileName()); dynamic_cast(m_Module)->SetSplinesFileName(m_SplinesFileSelector->GetFileName()); + + if (dynamic_cast(m_Module)->GetMaskMetrologyCorrection() != m_UseMaskMetCorr) dynamic_cast(m_Module)->EnableMaskMetrologyCorrection(m_UseMaskMetCorr); + + dynamic_cast(m_Module)->SetMaskMetrologyFileName(m_MaskMetrologyFileSelector->GetFileName()); dynamic_cast(m_Module)->SetUCSDOverride(m_UCSDOverride->IsOn()); return true; diff --git a/src/MModuleDepthCalibration2024.cxx b/src/MModuleDepthCalibration2024.cxx index 135db6dc..011e959b 100644 --- a/src/MModuleDepthCalibration2024.cxx +++ b/src/MModuleDepthCalibration2024.cxx @@ -30,6 +30,7 @@ // Standard libs: // ROOT libs: +#include "TMath.h" #include "TGClient.h" #include "TH1.h" @@ -171,6 +172,15 @@ bool MModuleDepthCalibration2024::Initialize() return false; } + if (m_MaskMetrologyEnabled == true) { + cout << "\n !!! Mask Metrology Enabled !!! \n" << endl; + m_MaskMetrologyFileIsLoaded = LoadMaskMetrologyFile(m_MaskMetrologyFile); + if (m_MaskMetrologyFile == false) { + cout << "Unable to open Metrology file" << endl; + return false; + } + } + MSupervisor* S = MSupervisor::GetSupervisor(); m_EnergyCalibration = (MModuleEnergyCalibrationUniversal*) S->GetAvailableModuleByXmlTag("EnergyCalibrationUniversal"); if (m_EnergyCalibration == nullptr) { @@ -202,6 +212,9 @@ void MModuleDepthCalibration2024::CreateExpos() } +///////////////////////////////////////////////////////////////////////////////// + + bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) { @@ -268,9 +281,25 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) // TODO: Calculate X and Y positions more rigorously using charge sharing. // Somewhat confusing notation: HVStrips run parallel to X-axis, so we calculate X position with LVStrips. - double Xpos = m_YPitches[DetID]*((m_NYStrips[DetID]/2.0) - ((double)LVStripID)); - double Ypos = m_XPitches[DetID]*((m_NXStrips[DetID]/2.0) - ((double)HVStripID)); - double Zpos = 0.0; + double Xpos, Ypos, Zpos; + if ( m_MaskMetrologyEnabled ) { + // If we are applying the mask metrology correction, first define two new readout elements to help determine the intersection of these two strips + MReadOutElementDoubleStrip R_LV = *dynamic_cast(LVSH->GetReadOutElement()); + MReadOutElementDoubleStrip R_HV = *dynamic_cast(HVSH->GetReadOutElement()); + + vector inter = GetStripIntersection(R_LV, R_HV); + Xpos = inter[0]; + Ypos = inter[1]; + //Xpos = m_YPitches[DetID]*((m_NYStrips[DetID]/2.0) - ((double)LVStripID)); + //Ypos = m_XPitches[DetID]*((m_NXStrips[DetID]/2.0) - ((double)HVStripID)); + + Zpos = 0.0; + } else { + Xpos = m_YPitches[DetID]*((m_NYStrips[DetID]/2.0) - ((double)LVStripID)); + Ypos = m_XPitches[DetID]*((m_NXStrips[DetID]/2.0) - ((double)HVStripID)); + Zpos = 0.0; + } + double Xsigma = m_YPitches[DetID]/sqrt(12.0); double Ysigma = m_XPitches[DetID]/sqrt(12.0); @@ -364,6 +393,9 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) } } +// cout << "Strip ID :" << LVStripID << " " << HVStripID << endl; +// cout << "Hit position: "<< Xpos << " " << Ypos << " " << Zpos << endl; + LocalPosition.SetXYZ(Xpos, Ypos, Zpos); LocalOrigin.SetXYZ(0.0,0.0,0.0); GlobalPosition = m_Detectors[DetID]->GetSensitiveVolume(0)->GetPositionInWorldVolume(LocalPosition); @@ -387,6 +419,10 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) return true; } + +///////////////////////////////////////////////////////////////////////////////// + + MStripHit* MModuleDepthCalibration2024::GetDominantStrip(vector& Strips, double& EnergyFraction) { double MaxEnergy = -numeric_limits::max(); // AZ: When both energies are zero (which shouldn't happen) we still pick one @@ -410,6 +446,10 @@ MStripHit* MModuleDepthCalibration2024::GetDominantStrip(vector& Str return MaxStrip; } + +///////////////////////////////////////////////////////////////////////////////// + + double MModuleDepthCalibration2024::GetTimingNoiseFWHM(int PixelCode, double Energy) { // Placeholder for determining the timing noise with energy, and possibly even on a pixel-by-pixel basis. @@ -428,6 +468,10 @@ double MModuleDepthCalibration2024::GetTimingNoiseFWHM(int PixelCode, double Ene return noiseFWHM; } + +///////////////////////////////////////////////////////////////////////////////// + + bool MModuleDepthCalibration2024::LoadCoeffsFile(MString FName) { // Read in the stretch and offset file, which should have a header line with information on the measurements: @@ -468,6 +512,10 @@ bool MModuleDepthCalibration2024::LoadCoeffsFile(MString FName) } + +///////////////////////////////////////////////////////////////////////////////// + + std::vector* MModuleDepthCalibration2024::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. @@ -475,7 +523,7 @@ std::vector* MModuleDepthCalibration2024::GetPixelCoeffs(int PixelCode) if( m_Coeffs.count(PixelCode) > 0 ){ return &m_Coeffs[PixelCode]; } else { - cout << "MModuleDepthCalibration2024::GetPixelCoeffs: cannot get stretch and offset; pixel code " << PixelCode << " not found." << endl; + // cout << "MModuleDepthCalibration2024::GetPixelCoeffs: cannot get stretch and offset; pixel code " << PixelCode << " not found." << endl; return nullptr; } } else { @@ -485,6 +533,10 @@ std::vector* MModuleDepthCalibration2024::GetPixelCoeffs(int PixelCode) } + +///////////////////////////////////////////////////////////////////////////////// + + vector MModuleDepthCalibration2024::norm_pdf(vector x, double mu, double sigma) { vector result; @@ -557,6 +609,83 @@ bool MModuleDepthCalibration2024::LoadSplinesFile(MString FName) } + +///////////////////////////////////////////////////////////////////////////////// + + +bool MModuleDepthCalibration2024::LoadMaskMetrologyFile(MString FName) +{ + //Read the Mask Metrology File + // Det ID, Strip ID (0-63), Side (l,h), x_mm, y_mm, z_mm, roll_deg, pitch_deg, yaw_deg + MFile F; + if( F.Open(FName) == false ){ + cout << "ERROR in MModuleDepthCalibration2024::LoadMaskMetrologyFile: failed to open metrology file." << endl; + return false; + } else { + MString Line; + while( F.ReadLine( Line ) ){ + if ( Line.BeginsWith('#') ){ + } + else { + std::vector Tokens = Line.Tokenize(","); + if( Tokens.size() == 9 ){ + //Define the readout element to track det ID, strip ID, and lv/hv + MReadOutElementDoubleStrip R; + R.SetDetectorID(Tokens[0].ToInt()); + R.SetStripID(Tokens[1].ToInt()); + R.IsLowVoltageStrip((Tokens[2].ToString() == "p") || (Tokens[2].ToString() == "l")); + double Strip_MetX = Tokens[3].ToDouble()/10; //convert to cm + double Strip_MetY = Tokens[4].ToDouble()/10; //convert to cm + double Strip_MetZ = Tokens[5].ToDouble()/10; //convert to cm + double Strip_Roll = Tokens[6].ToDouble(); + double Strip_Pitch = Tokens[7].ToDouble(); + double Strip_Yaw = Tokens[8].ToDouble(); + vector maskmet; + maskmet.push_back(Strip_MetX); maskmet.push_back(Strip_MetY); maskmet.push_back(Strip_MetZ); + maskmet.push_back(Strip_Roll); maskmet.push_back(Strip_Pitch); maskmet.push_back(Strip_Yaw); + + //make the map that defines the metrology info for each readout element + m_MaskMetrology[R] = maskmet; + } + } + } + + F.Close(); + + } + + return true; + +} + + +///////////////////////////////////////////////////////////////////////////////// + + +vector MModuleDepthCalibration2024::GetStripIntersection(MReadOutElementDoubleStrip R_LVStrip, MReadOutElementDoubleStrip R_HVStrip){ +// Function to find the intersection between two strips based on the Mask Metrology + + vector lv_strip_met = m_MaskMetrology[R_LVStrip]; + vector hv_strip_met = m_MaskMetrology[R_HVStrip]; + + //Find the x position of two lines represented by the dominate strips: + // LVstrip is centered at (x,y,z) = (lv_strip_met[0], lv_strip_met[1], lv_strip_met[2]) + // and is approximately parallel to the y axis, but rotated at angle lv_strip_met[3] + // around the z axis of the detector + // HVstrip is centered at (x,y,z) = (hv_strip_met[0], hv_strip_met[1], hv_strip_met[2]) + // and is approximately parallel to the x axis, but rotated at angle hv_strip_met[3] + // around the z axis of the detector + double x_intercept = (hv_strip_met[0]*tan(hv_strip_met[3]*TMath::DegToRad()) - lv_strip_met[0]/tan(lv_strip_met[3]*TMath::DegToRad()) - lv_strip_met[1] + hv_strip_met[1])/(tan(hv_strip_met[3]*TMath::DegToRad())-1/tan(lv_strip_met[3]*TMath::DegToRad())); + + double y_intercept = (x_intercept - hv_strip_met[0])*tan(hv_strip_met[3]*TMath::DegToRad()) + hv_strip_met[1]; + + return {x_intercept, y_intercept}; +} + + +///////////////////////////////////////////////////////////////////////////////// + + int MModuleDepthCalibration2024::GetHitGrade(MHit* H){ // Function for choosing which Depth-to-CTD relation to use for a given event. // At time of writing, intention is to choose a CTD based on sub-pixel region determined via charge sharing (Event "grade"). @@ -672,6 +801,10 @@ int MModuleDepthCalibration2024::GetHitGrade(MHit* H){ return return_value; } + +///////////////////////////////////////////////////////////////////////////////// + + bool MModuleDepthCalibration2024::AddDepthCTD(vector depthvec, vector> ctdarr, int DetID, unordered_map>& DepthGrid, unordered_map>>& CTDMap){ // Saves a CTD array, basically allowing for multiple CTDs as a function of depth @@ -714,6 +847,9 @@ bool MModuleDepthCalibration2024::AddDepthCTD(vector depthvec, vector MModuleDepthCalibration2024::GetCTD(int DetID, int Grade) { // Retrieves the appropriate CTD vector given the Detector ID and Event Grade passed @@ -738,6 +874,10 @@ vector MModuleDepthCalibration2024::GetCTD(int DetID, int Grade) } } + +///////////////////////////////////////////////////////////////////////////////// + + vector MModuleDepthCalibration2024::GetDepth(int DetID) { // Retrieves the appropriate CTD vector given the Detector ID and Event Grade passed @@ -757,6 +897,10 @@ vector MModuleDepthCalibration2024::GetDepth(int DetID) } } + +///////////////////////////////////////////////////////////////////////////////// + + void MModuleDepthCalibration2024::ShowOptionsGUI() { // Show the options GUI - or do nothing @@ -766,6 +910,9 @@ void MModuleDepthCalibration2024::ShowOptionsGUI() } +///////////////////////////////////////////////////////////////////////////////// + + bool MModuleDepthCalibration2024::ReadXmlConfiguration(MXmlNode* Node) { //! Read the configuration data from an XML node @@ -780,6 +927,16 @@ bool MModuleDepthCalibration2024::ReadXmlConfiguration(MXmlNode* Node) m_SplinesFile = SplinesFileNameNode->GetValue(); } + MXmlNode* MasKMetrologyNode = Node->GetNode("MaskMetrology"); + if( MasKMetrologyNode != NULL ){ + m_MaskMetrologyEnabled = (bool) MasKMetrologyNode->GetValueAsInt(); + } + + MXmlNode* MaskMetrologyFileNameNode = Node->GetNode("MaskMetrologyFileName"); + if (MaskMetrologyFileNameNode != 0) { + m_MaskMetrologyFile = MaskMetrologyFileNameNode->GetValue(); + } + MXmlNode* UCSDOverrideNode = Node->GetNode("UCSDOverride"); if( UCSDOverrideNode != NULL ){ m_UCSDOverride = (bool) UCSDOverrideNode->GetValueAsInt(); @@ -798,6 +955,8 @@ MXmlNode* MModuleDepthCalibration2024::CreateXmlConfiguration() MXmlNode* Node = new MXmlNode(0,m_XmlTag); new MXmlNode(Node, "CoeffsFileName", m_CoeffsFile); new MXmlNode(Node, "SplinesFileName", m_SplinesFile); + new MXmlNode(Node, "MaskMetrology", (unsigned int) m_MaskMetrologyEnabled); + new MXmlNode(Node, "MaskMetrologyFileName", m_MaskMetrologyFile); new MXmlNode(Node, "UCSDOverride",(unsigned int) m_UCSDOverride); return Node; From 1f2163f303a9ce2688f734755de91369c150ac8d Mon Sep 17 00:00:00 2001 From: ckierans Date: Tue, 2 Sep 2025 09:31:07 -0400 Subject: [PATCH 02/22] CHG: Updated x,y definition with mask rotation, added print statements --- src/MModuleDepthCalibration2024.cxx | 20 +++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/src/MModuleDepthCalibration2024.cxx b/src/MModuleDepthCalibration2024.cxx index 011e959b..baccd430 100644 --- a/src/MModuleDepthCalibration2024.cxx +++ b/src/MModuleDepthCalibration2024.cxx @@ -284,15 +284,21 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) double Xpos, Ypos, Zpos; if ( m_MaskMetrologyEnabled ) { // If we are applying the mask metrology correction, first define two new readout elements to help determine the intersection of these two strips + int HVStripID_flipped = 63 - HVStripID; + HVSH->SetStripID(HVStripID_flipped); MReadOutElementDoubleStrip R_LV = *dynamic_cast(LVSH->GetReadOutElement()); MReadOutElementDoubleStrip R_HV = *dynamic_cast(HVSH->GetReadOutElement()); vector inter = GetStripIntersection(R_LV, R_HV); Xpos = inter[0]; Ypos = inter[1]; - //Xpos = m_YPitches[DetID]*((m_NYStrips[DetID]/2.0) - ((double)LVStripID)); - //Ypos = m_XPitches[DetID]*((m_NXStrips[DetID]/2.0) - ((double)HVStripID)); - + + cout<<"StripLV: "< MModuleDepthCalibration2024::norm_pdf(vector x, double mu return result; } + +///////////////////////////////////////////////////////////////////////////////// + + bool MModuleDepthCalibration2024::LoadSplinesFile(MString FName) { //when invert flag is set to true, the splines returned are CTD->Depth @@ -675,9 +685,9 @@ vector MModuleDepthCalibration2024::GetStripIntersection(MReadOutElement // HVstrip is centered at (x,y,z) = (hv_strip_met[0], hv_strip_met[1], hv_strip_met[2]) // and is approximately parallel to the x axis, but rotated at angle hv_strip_met[3] // around the z axis of the detector - double x_intercept = (hv_strip_met[0]*tan(hv_strip_met[3]*TMath::DegToRad()) - lv_strip_met[0]/tan(lv_strip_met[3]*TMath::DegToRad()) - lv_strip_met[1] + hv_strip_met[1])/(tan(hv_strip_met[3]*TMath::DegToRad())-1/tan(lv_strip_met[3]*TMath::DegToRad())); + double x_intercept = (hv_strip_met[0]*tan(-hv_strip_met[5]*TMath::DegToRad()) - lv_strip_met[0]/tan(lv_strip_met[5]*TMath::DegToRad()) - lv_strip_met[1] + hv_strip_met[1])/(tan(-hv_strip_met[5]*TMath::DegToRad())-1/tan(lv_strip_met[5]*TMath::DegToRad())); - double y_intercept = (x_intercept - hv_strip_met[0])*tan(hv_strip_met[3]*TMath::DegToRad()) + hv_strip_met[1]; + double y_intercept = (x_intercept - hv_strip_met[0])*tan(hv_strip_met[5]*TMath::DegToRad()) + hv_strip_met[1]; return {x_intercept, y_intercept}; } From 5577c2ec4d742f61d9057738ae338fd8f7c7cdba Mon Sep 17 00:00:00 2001 From: ckierans Date: Mon, 27 Oct 2025 12:37:53 -0400 Subject: [PATCH 03/22] Updated metrology with correct understanding of rotation angles, now works for 2 detectors based on latest files from Aldo --- src/MModuleDepthCalibration2024.cxx | 38 +++++++++++++++-------------- 1 file changed, 20 insertions(+), 18 deletions(-) diff --git a/src/MModuleDepthCalibration2024.cxx b/src/MModuleDepthCalibration2024.cxx index baccd430..a6c3bd70 100644 --- a/src/MModuleDepthCalibration2024.cxx +++ b/src/MModuleDepthCalibration2024.cxx @@ -219,7 +219,8 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) { if (Event->GetGuardRingVeto()==true) { - + //TODO: Handle events with GR vetos + Event->SetDepthCalibrationIncomplete(); return false; @@ -280,27 +281,28 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) int PixelCode = 10000*DetID + 100*LVStripID + HVStripID; // TODO: Calculate X and Y positions more rigorously using charge sharing. - // Somewhat confusing notation: HVStrips run parallel to X-axis, so we calculate X position with LVStrips. double Xpos, Ypos, Zpos; + if ( m_MaskMetrologyEnabled ) { // If we are applying the mask metrology correction, first define two new readout elements to help determine the intersection of these two strips + + + //TODO: Aldo's file has the HV strip ID flipped (as evident in the x,y position reported in the metrology files), so this is a placeholder until the files are fixed. int HVStripID_flipped = 63 - HVStripID; - HVSH->SetStripID(HVStripID_flipped); + MStripHit* HVSH_flipped = HVSH; + HVSH_flipped->SetStripID(HVStripID_flipped); MReadOutElementDoubleStrip R_LV = *dynamic_cast(LVSH->GetReadOutElement()); - MReadOutElementDoubleStrip R_HV = *dynamic_cast(HVSH->GetReadOutElement()); - + MReadOutElementDoubleStrip R_HV = *dynamic_cast(HVSH_flipped->GetReadOutElement()); + + //find the intercept of the two dominate strips based on the mask metrology vector inter = GetStripIntersection(R_LV, R_HV); Xpos = inter[0]; Ypos = inter[1]; - cout<<"StripLV: "< MModuleDepthCalibration2024::GetStripIntersection(MReadOutElement //Find the x position of two lines represented by the dominate strips: // LVstrip is centered at (x,y,z) = (lv_strip_met[0], lv_strip_met[1], lv_strip_met[2]) - // and is approximately parallel to the y axis, but rotated at angle lv_strip_met[3] + // and is approximately parallel to the y axis, but rotated at angle lv_strip_met[5] // around the z axis of the detector // HVstrip is centered at (x,y,z) = (hv_strip_met[0], hv_strip_met[1], hv_strip_met[2]) - // and is approximately parallel to the x axis, but rotated at angle hv_strip_met[3] + // and is approximately parallel to the x axis, but rotated at angle (hv_strip_met[5] - pi/2) // around the z axis of the detector - double x_intercept = (hv_strip_met[0]*tan(-hv_strip_met[5]*TMath::DegToRad()) - lv_strip_met[0]/tan(lv_strip_met[5]*TMath::DegToRad()) - lv_strip_met[1] + hv_strip_met[1])/(tan(-hv_strip_met[5]*TMath::DegToRad())-1/tan(lv_strip_met[5]*TMath::DegToRad())); + double x_intercept = (hv_strip_met[0]*tan((hv_strip_met[5]-90)*TMath::DegToRad()) - lv_strip_met[0]/tan(lv_strip_met[5]*TMath::DegToRad()) - lv_strip_met[1] + hv_strip_met[1])/(tan((hv_strip_met[5]-90)*TMath::DegToRad())-1/tan(lv_strip_met[5]*TMath::DegToRad())); - double y_intercept = (x_intercept - hv_strip_met[0])*tan(hv_strip_met[5]*TMath::DegToRad()) + hv_strip_met[1]; + double y_intercept = (x_intercept - hv_strip_met[0])*tan((hv_strip_met[5]-90)*TMath::DegToRad()) + hv_strip_met[1]; return {x_intercept, y_intercept}; } From d05aa6a98c9e8e218beb24eec33040592bbe001f Mon Sep 17 00:00:00 2001 From: julianmgerber <95662611+julianmgerber@users.noreply.github.com> Date: Tue, 18 Nov 2025 10:00:39 -0800 Subject: [PATCH 04/22] Strip Pairing Chi Square Algorithm Updated (#61) * Merging select files from StripPairingTesting into StripPairingCleanUp so I can do pull request * Clean up comments. Remove long debugging statements for readability * Calling charge trapping correction function in chi square calculation and hit population * Cleaning up comments * Reformatting strip pairing module * Update copywrite. Change name of module. * Fixing documentation in header files * Change how strip pairing chi square is read out to .dat file * Adding comments for all the nested vectors * Adding new version to StreamDat that reads out LV and HV energy of each hit * Adding Sean's new strip pairing expos to chi square version of strip pairing * Changing all instances of RedChiSquare to ReducedChiSquare for clarity * Flags for multiple hits on a single strip are now defined on the MHit level * Initializing m_StripHitMultipleTimesX and m_StripHitMultipleTimesY to false * Changing StreamDAT version from Version 1 to Version 3 * Rebased and fixed merge conflict with strip pairing --------- Co-authored-by: Julian Gerber --- include/MReadOutAssembly.h | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/include/MReadOutAssembly.h b/include/MReadOutAssembly.h index 1a8b0aa6..23db7b7c 100644 --- a/include/MReadOutAssembly.h +++ b/include/MReadOutAssembly.h @@ -454,8 +454,7 @@ class MReadOutAssembly : public MReadOutSequence bool m_StripHitBelowThreshold; MString m_StripHitBelowThresholdString; - - + //! True if event has been filtered out bool m_FilteredOut; From 6ba6e6bf15f93ee1a4afc3b210912e4669182166 Mon Sep 17 00:00:00 2001 From: ckierans Date: Tue, 9 Dec 2025 22:59:42 -0500 Subject: [PATCH 05/22] Resolving merge conflict with Sean's recent PR --- include/MGUIOptionsDepthCalibration2024.h | 7 + include/MModuleDepthCalibration2024.h | 42 +++++- src/MGUIOptionsDepthCalibration2024.cxx | 54 ++++++- src/MModuleDepthCalibration2024.cxx | 170 +++++++++++++++++++++- 4 files changed, 256 insertions(+), 17 deletions(-) diff --git a/include/MGUIOptionsDepthCalibration2024.h b/include/MGUIOptionsDepthCalibration2024.h index 15c72ed6..6c14df24 100644 --- a/include/MGUIOptionsDepthCalibration2024.h +++ b/include/MGUIOptionsDepthCalibration2024.h @@ -74,6 +74,13 @@ class MGUIOptionsDepthCalibration2024 : public MGUIOptions //! Select spline file to load, splines will convert CTD->Depth MGUIEFileSelector* m_SplinesFileSelector; + //! Select mask metrology file to load. This gives the translation and rotation for each strip in the detector frame + MGUIEFileSelector* m_MaskMetrologyFileSelector; + + int m_UseMaskMetCorr; + TGCheckButton* m_MaskMetModeCB; + enum ButtonIDs {c_MetFile}; + //! Check button if working with the Card Cage at UCSD TGCheckButton* m_UCSDOverride; diff --git a/include/MModuleDepthCalibration2024.h b/include/MModuleDepthCalibration2024.h index 17237b4f..6a697761 100644 --- a/include/MModuleDepthCalibration2024.h +++ b/include/MModuleDepthCalibration2024.h @@ -1,4 +1,4 @@ -/* +/*:q * MModuleDepthCalibration2024.h * * Copyright (C) 2008-2008 by Andreas Zoglauer. @@ -72,6 +72,16 @@ class MModuleDepthCalibration2024 : public MModule //! Get filename for CTD->Depth splines MString GetSplinesFileName() const {return m_SplinesFile;} + //! Enable/Disable Preamp Temp Correction + void EnableMaskMetrologyCorrection(bool X) {m_MaskMetrologyEnabled = X;} + //! Get coincidence merging true/false + bool GetMaskMetrologyCorrection() const { return m_MaskMetrologyEnabled; } + + //! Set filename for mask metrology + void SetMaskMetrologyFileName( const MString& FileName) {m_MaskMetrologyFile = FileName;} + //! Get filename for CTD->Depth splines + MString GetMaskMetrologyFileName() const {return m_MaskMetrologyFile;} + //! Set whether the data came from the card cage at UCSD void SetUCSDOverride( bool Override ) {m_UCSDOverride = Override;} //! Get whether the data came from the card cage at UCSD @@ -91,24 +101,42 @@ class MModuleDepthCalibration2024 : public MModule protected: //! Returns the strip with most energy from vector Strips, also gives back the energy fraction MStripHit* GetDominantStrip(std::vector& Strips, double& EnergyFraction); + //! Retrieve the appropriate Depth values given the DetID - vector GetDepth(int DetID); + vector GetDepth(int DetID); + //! Retrieve the appropriate CTD values given the DetID and Grade vector GetCTD(int DetID, int Grade); + //! Retrieve the appropriate depth-to-CTD spline given the DetID and Grade TSpline3* GetSpline(int DetID, int Grade); //! Normal distribution vector norm_pdf(vector x, double mu, double sigma); - //! Adds a Depth-to-CTD relation - bool AddDepthCTD(vector Depth, vector> CTDArr, int DetID, unordered_map>& DepthGrid, unordered_map>>& CTDMap, unordered_map>& SplineMap, unsigned int NPoints); + + //! Adds a Depth-to-CTD relation + bool AddDepthCTD(vector Depth, vector> CTDArr, int DetID, unordered_map>& DepthGrid, unordered_map>>& CTDMap, unordered_map>& SplineMap, unsigned int NPoints); + //! Determine the Grade (geometry of charge sharing) of the Hit int GetHitGrade(MHit* H); + //! Load in the specified coefficients file bool LoadCoeffsFile(MString FName); + //! Return the coefficients for a pixel vector* GetPixelCoeffs(int PixelCode); + //! Load the splines file bool LoadSplinesFile(MString FName); + + //! Mask Metrology Correction + bool m_MaskMetrologyEnabled; + + //! Load the metrology mask file + bool LoadMaskMetrologyFile(MString FName); + + //! Get the x, y position of the intersection of two strips based on the Metrology Mask + vector GetStripIntersection(MReadOutElementDoubleStrip LVStrip, MReadOutElementDoubleStrip HVStrip); + //! Get the timing FWHM noise for the specified pixel and Energy double GetTimingNoiseFWHM(int PixelCode, double Energy); @@ -150,6 +178,12 @@ class MModuleDepthCalibration2024 : public MModule unordered_map> m_SplineMap; bool m_SplinesFileIsLoaded; bool m_CoeffsFileIsLoaded; + bool m_MaskMetrologyFileIsLoaded; + + // The Mask Metrology + MString m_MaskMetrologyFile; + map> m_MaskMetrology; + // boolean for use with the card cage at UCSD since it tags all events as detector 11 bool m_UCSDOverride; diff --git a/src/MGUIOptionsDepthCalibration2024.cxx b/src/MGUIOptionsDepthCalibration2024.cxx index cd6c07ec..9ead8652 100644 --- a/src/MGUIOptionsDepthCalibration2024.cxx +++ b/src/MGUIOptionsDepthCalibration2024.cxx @@ -75,8 +75,30 @@ void MGUIOptionsDepthCalibration2024::Create() m_SplinesFileSelector = new MGUIEFileSelector(m_OptionsFrame, "Select a splines file:", dynamic_cast(m_Module)->GetSplinesFileName()); m_SplinesFileSelector->SetFileType("splines", "*.csv"); - TGLayoutHints* Label2Layout = new TGLayoutHints(kLHintsTop | kLHintsCenterX | kLHintsExpandX, 10, 10, 10, 10); - m_OptionsFrame->AddFrame(m_SplinesFileSelector, Label2Layout); +// TGLayoutHints* Label2Layout = new TGLayoutHints(kLHintsTop | kLHintsCenterX | kLHintsExpandX, 10, 10, 10, 10); + m_OptionsFrame->AddFrame(m_SplinesFileSelector, LabelLayout); + + m_MaskMetModeCB = new TGCheckButton(m_OptionsFrame, "Enable mask metrology correction and read calibration from file:", c_MetFile); + m_MaskMetModeCB->SetState((dynamic_cast(m_Module)->GetMaskMetrologyCorrection() == 1) ? kButtonDown : kButtonUp); + m_MaskMetModeCB->Associate(this); + m_OptionsFrame->AddFrame(m_MaskMetModeCB, LabelLayout); + + m_UseMaskMetCorr = dynamic_cast(m_Module)->GetMaskMetrologyCorrection(); + +// TGLayoutHints* FileLabelLayout = new TGLayoutHints(kLHintsTop | kLHintsExpandX, m_FontScaler*65 + 21*m_FontScaler, m_FontScaler*65, 0, 2*m_FontScaler); + + m_MaskMetrologyFileSelector = new MGUIEFileSelector(m_OptionsFrame, "", dynamic_cast(m_Module)->GetMaskMetrologyFileName()); + m_MaskMetrologyFileSelector->SetFileType("metrology", "*.csv"); + m_OptionsFrame->AddFrame(m_MaskMetrologyFileSelector, LabelLayout); + + + if (m_UseMaskMetCorr) { + m_MaskMetrologyFileSelector->SetEnabled(true); + } else { + m_MaskMetrologyFileSelector->SetEnabled(false); + } + + m_UCSDOverride = new TGCheckButton(m_OptionsFrame, "Check this box if you're using the card cage at UCSD", 1); m_UCSDOverride->SetOn(dynamic_cast(m_Module)->GetUCSDOverride()); @@ -97,16 +119,28 @@ bool MGUIOptionsDepthCalibration2024::ProcessMessage(long Message, long Paramete bool Status = true; switch (GET_MSG(Message)) { - case kC_COMMAND: - switch (GET_SUBMSG(Message)) { - case kCM_BUTTON: - break; - default: + case kC_COMMAND: + switch (GET_SUBMSG(Message)) { + case kCM_BUTTON: + break; + case kCM_CHECKBUTTON: + switch (Parameter1) { + case c_MetFile: + if (m_MaskMetModeCB->GetState() == kButtonDown) { + m_UseMaskMetCorr = 1; + m_MaskMetrologyFileSelector->SetEnabled(true); + } else if (m_MaskMetModeCB->GetState() == kButtonUp) { + m_UseMaskMetCorr = 0; + m_MaskMetrologyFileSelector->SetEnabled(false); + } break; } - break; default: break; + } + break; + default: + break; } if (Status == false) { @@ -127,6 +161,10 @@ bool MGUIOptionsDepthCalibration2024::OnApply() dynamic_cast(m_Module)->SetCoeffsFileName(m_CoeffsFileSelector->GetFileName()); dynamic_cast(m_Module)->SetSplinesFileName(m_SplinesFileSelector->GetFileName()); + + if (dynamic_cast(m_Module)->GetMaskMetrologyCorrection() != m_UseMaskMetCorr) dynamic_cast(m_Module)->EnableMaskMetrologyCorrection(m_UseMaskMetCorr); + + dynamic_cast(m_Module)->SetMaskMetrologyFileName(m_MaskMetrologyFileSelector->GetFileName()); dynamic_cast(m_Module)->SetUCSDOverride(m_UCSDOverride->IsOn()); return true; diff --git a/src/MModuleDepthCalibration2024.cxx b/src/MModuleDepthCalibration2024.cxx index 46de19ae..1f26bca6 100644 --- a/src/MModuleDepthCalibration2024.cxx +++ b/src/MModuleDepthCalibration2024.cxx @@ -30,6 +30,7 @@ // Standard libs: // ROOT libs: +#include "TMath.h" #include "TGClient.h" #include "TH1.h" @@ -171,6 +172,15 @@ bool MModuleDepthCalibration2024::Initialize() return false; } + if (m_MaskMetrologyEnabled == true) { + cout << "\n !!! Mask Metrology Enabled !!! \n" << endl; + m_MaskMetrologyFileIsLoaded = LoadMaskMetrologyFile(m_MaskMetrologyFile); + if (m_MaskMetrologyFile == false) { + cout << "Unable to open Metrology file" << endl; + return false; + } + } + MSupervisor* S = MSupervisor::GetSupervisor(); m_EnergyCalibration = (MModuleEnergyCalibrationUniversal*) S->GetAvailableModuleByXmlTag("EnergyCalibrationUniversal"); if (m_EnergyCalibration == nullptr) { @@ -202,6 +212,9 @@ void MModuleDepthCalibration2024::CreateExpos() } +///////////////////////////////////////////////////////////////////////////////// + + bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) { @@ -268,9 +281,25 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) // TODO: Calculate X and Y positions more rigorously using charge sharing. // Somewhat confusing notation: HVStrips run parallel to X-axis, so we calculate X position with LVStrips. - double Xpos = m_YPitches[DetID]*((m_NYStrips[DetID]/2.0) - ((double)LVStripID)); - double Ypos = m_XPitches[DetID]*((m_NXStrips[DetID]/2.0) - ((double)HVStripID)); - double Zpos = 0.0; + double Xpos, Ypos, Zpos; + if ( m_MaskMetrologyEnabled ) { + // If we are applying the mask metrology correction, first define two new readout elements to help determine the intersection of these two strips + MReadOutElementDoubleStrip R_LV = *dynamic_cast(LVSH->GetReadOutElement()); + MReadOutElementDoubleStrip R_HV = *dynamic_cast(HVSH->GetReadOutElement()); + + vector inter = GetStripIntersection(R_LV, R_HV); + Xpos = inter[0]; + Ypos = inter[1]; + //Xpos = m_YPitches[DetID]*((m_NYStrips[DetID]/2.0) - ((double)LVStripID)); + //Ypos = m_XPitches[DetID]*((m_NXStrips[DetID]/2.0) - ((double)HVStripID)); + + Zpos = 0.0; + } else { + Xpos = m_YPitches[DetID]*((m_NYStrips[DetID]/2.0) - ((double)LVStripID)); + Ypos = m_XPitches[DetID]*((m_NXStrips[DetID]/2.0) - ((double)HVStripID)); + Zpos = 0.0; + } + double Xsigma = m_YPitches[DetID]/sqrt(12.0); double Ysigma = m_XPitches[DetID]/sqrt(12.0); @@ -364,6 +393,9 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) } } +// cout << "Strip ID :" << LVStripID << " " << HVStripID << endl; +// cout << "Hit position: "<< Xpos << " " << Ypos << " " << Zpos << endl; + LocalPosition.SetXYZ(Xpos, Ypos, Zpos); LocalOrigin.SetXYZ(0.0,0.0,0.0); GlobalPosition = m_Detectors[DetID]->GetSensitiveVolume(0)->GetPositionInWorldVolume(LocalPosition); @@ -387,6 +419,10 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) return true; } + +///////////////////////////////////////////////////////////////////////////////// + + MStripHit* MModuleDepthCalibration2024::GetDominantStrip(vector& Strips, double& EnergyFraction) { double MaxEnergy = -numeric_limits::max(); // AZ: When both energies are zero (which shouldn't happen) we still pick one @@ -410,6 +446,10 @@ MStripHit* MModuleDepthCalibration2024::GetDominantStrip(vector& Str return MaxStrip; } + +///////////////////////////////////////////////////////////////////////////////// + + double MModuleDepthCalibration2024::GetTimingNoiseFWHM(int PixelCode, double Energy) { // Placeholder for determining the timing noise with energy, and possibly even on a pixel-by-pixel basis. @@ -428,6 +468,10 @@ double MModuleDepthCalibration2024::GetTimingNoiseFWHM(int PixelCode, double Ene return noiseFWHM; } + +///////////////////////////////////////////////////////////////////////////////// + + bool MModuleDepthCalibration2024::LoadCoeffsFile(MString FName) { // Read in the stretch and offset file, which should have a header line with information on the measurements: @@ -468,6 +512,10 @@ bool MModuleDepthCalibration2024::LoadCoeffsFile(MString FName) } + +///////////////////////////////////////////////////////////////////////////////// + + std::vector* MModuleDepthCalibration2024::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. @@ -475,7 +523,7 @@ std::vector* MModuleDepthCalibration2024::GetPixelCoeffs(int PixelCode) if( m_Coeffs.count(PixelCode) > 0 ){ return &m_Coeffs[PixelCode]; } else { - cout << "MModuleDepthCalibration2024::GetPixelCoeffs: cannot get stretch and offset; pixel code " << PixelCode << " not found." << endl; + // cout << "MModuleDepthCalibration2024::GetPixelCoeffs: cannot get stretch and offset; pixel code " << PixelCode << " not found." << endl; return nullptr; } } else { @@ -485,6 +533,10 @@ std::vector* MModuleDepthCalibration2024::GetPixelCoeffs(int PixelCode) } + +///////////////////////////////////////////////////////////////////////////////// + + vector MModuleDepthCalibration2024::norm_pdf(vector x, double mu, double sigma) { vector result; @@ -544,6 +596,83 @@ bool MModuleDepthCalibration2024::LoadSplinesFile(MString FName) } + +///////////////////////////////////////////////////////////////////////////////// + + +bool MModuleDepthCalibration2024::LoadMaskMetrologyFile(MString FName) +{ + //Read the Mask Metrology File + // Det ID, Strip ID (0-63), Side (l,h), x_mm, y_mm, z_mm, roll_deg, pitch_deg, yaw_deg + MFile F; + if( F.Open(FName) == false ){ + cout << "ERROR in MModuleDepthCalibration2024::LoadMaskMetrologyFile: failed to open metrology file." << endl; + return false; + } else { + MString Line; + while( F.ReadLine( Line ) ){ + if ( Line.BeginsWith('#') ){ + } + else { + std::vector Tokens = Line.Tokenize(","); + if( Tokens.size() == 9 ){ + //Define the readout element to track det ID, strip ID, and lv/hv + MReadOutElementDoubleStrip R; + R.SetDetectorID(Tokens[0].ToInt()); + R.SetStripID(Tokens[1].ToInt()); + R.IsLowVoltageStrip((Tokens[2].ToString() == "p") || (Tokens[2].ToString() == "l")); + double Strip_MetX = Tokens[3].ToDouble()/10; //convert to cm + double Strip_MetY = Tokens[4].ToDouble()/10; //convert to cm + double Strip_MetZ = Tokens[5].ToDouble()/10; //convert to cm + double Strip_Roll = Tokens[6].ToDouble(); + double Strip_Pitch = Tokens[7].ToDouble(); + double Strip_Yaw = Tokens[8].ToDouble(); + vector maskmet; + maskmet.push_back(Strip_MetX); maskmet.push_back(Strip_MetY); maskmet.push_back(Strip_MetZ); + maskmet.push_back(Strip_Roll); maskmet.push_back(Strip_Pitch); maskmet.push_back(Strip_Yaw); + + //make the map that defines the metrology info for each readout element + m_MaskMetrology[R] = maskmet; + } + } + } + + F.Close(); + + } + + return true; + +} + + +///////////////////////////////////////////////////////////////////////////////// + + +vector MModuleDepthCalibration2024::GetStripIntersection(MReadOutElementDoubleStrip R_LVStrip, MReadOutElementDoubleStrip R_HVStrip){ +// Function to find the intersection between two strips based on the Mask Metrology + + vector lv_strip_met = m_MaskMetrology[R_LVStrip]; + vector hv_strip_met = m_MaskMetrology[R_HVStrip]; + + //Find the x position of two lines represented by the dominate strips: + // LVstrip is centered at (x,y,z) = (lv_strip_met[0], lv_strip_met[1], lv_strip_met[2]) + // and is approximately parallel to the y axis, but rotated at angle lv_strip_met[3] + // around the z axis of the detector + // HVstrip is centered at (x,y,z) = (hv_strip_met[0], hv_strip_met[1], hv_strip_met[2]) + // and is approximately parallel to the x axis, but rotated at angle hv_strip_met[3] + // around the z axis of the detector + double x_intercept = (hv_strip_met[0]*tan(hv_strip_met[3]*TMath::DegToRad()) - lv_strip_met[0]/tan(lv_strip_met[3]*TMath::DegToRad()) - lv_strip_met[1] + hv_strip_met[1])/(tan(hv_strip_met[3]*TMath::DegToRad())-1/tan(lv_strip_met[3]*TMath::DegToRad())); + + double y_intercept = (x_intercept - hv_strip_met[0])*tan(hv_strip_met[3]*TMath::DegToRad()) + hv_strip_met[1]; + + return {x_intercept, y_intercept}; +} + + +///////////////////////////////////////////////////////////////////////////////// + + int MModuleDepthCalibration2024::GetHitGrade(MHit* H){ // Function for choosing which Depth-to-CTD relation to use for a given event. // At time of writing, intention is to choose a CTD based on sub-pixel region determined via charge sharing (Event "grade"). @@ -658,9 +787,11 @@ int MModuleDepthCalibration2024::GetHitGrade(MHit* H){ return return_value; } + +///////////////////////////////////////////////////////////////////////////////// + bool MModuleDepthCalibration2024::AddDepthCTD(vector Depth, vector> CTDArr, int DetID, unordered_map>& DepthGrid, unordered_map>>& CTDMap, unordered_map>& SplineMap, unsigned int NPoints) { - // Saves a CTD array, basically allowing for multiple CTDs as a function of depth // Depth: list of simulated depth values // CTDArr: vector of vectors of simulated CTD values. Each vector of CTDs must be the same length as Depth @@ -757,6 +888,9 @@ bool MModuleDepthCalibration2024::AddDepthCTD(vector Depth, vector MModuleDepthCalibration2024::GetCTD(int DetID, int Grade) { // Retrieves the appropriate CTD vector given the Detector ID and Event Grade passed @@ -780,6 +914,10 @@ vector MModuleDepthCalibration2024::GetCTD(int DetID, int Grade) } } + +///////////////////////////////////////////////////////////////////////////////// + + vector MModuleDepthCalibration2024::GetDepth(int DetID) { // Retrieves the appropriate CTD vector given the Detector ID and Event Grade passed @@ -800,6 +938,9 @@ vector MModuleDepthCalibration2024::GetDepth(int DetID) } +///////////////////////////////////////////////////////////////////////////////// + + TSpline3* MModuleDepthCalibration2024::GetSpline(int DetID, int Grade) { // Retrieves the appropriate depth->CTD spline given the Detector ID and Event Grade passed @@ -824,6 +965,10 @@ TSpline3* MModuleDepthCalibration2024::GetSpline(int DetID, int Grade) } } + +///////////////////////////////////////////////////////////////////////////////// + + void MModuleDepthCalibration2024::ShowOptionsGUI() { // Show the options GUI - or do nothing @@ -833,6 +978,9 @@ void MModuleDepthCalibration2024::ShowOptionsGUI() } +///////////////////////////////////////////////////////////////////////////////// + + bool MModuleDepthCalibration2024::ReadXmlConfiguration(MXmlNode* Node) { //! Read the configuration data from an XML node @@ -847,6 +995,16 @@ bool MModuleDepthCalibration2024::ReadXmlConfiguration(MXmlNode* Node) m_SplinesFile = SplinesFileNameNode->GetValue(); } + MXmlNode* MasKMetrologyNode = Node->GetNode("MaskMetrology"); + if( MasKMetrologyNode != NULL ){ + m_MaskMetrologyEnabled = (bool) MasKMetrologyNode->GetValueAsInt(); + } + + MXmlNode* MaskMetrologyFileNameNode = Node->GetNode("MaskMetrologyFileName"); + if (MaskMetrologyFileNameNode != 0) { + m_MaskMetrologyFile = MaskMetrologyFileNameNode->GetValue(); + } + MXmlNode* UCSDOverrideNode = Node->GetNode("UCSDOverride"); if( UCSDOverrideNode != NULL ){ m_UCSDOverride = (bool) UCSDOverrideNode->GetValueAsInt(); @@ -865,6 +1023,8 @@ MXmlNode* MModuleDepthCalibration2024::CreateXmlConfiguration() MXmlNode* Node = new MXmlNode(0,m_XmlTag); new MXmlNode(Node, "CoeffsFileName", m_CoeffsFile); new MXmlNode(Node, "SplinesFileName", m_SplinesFile); + new MXmlNode(Node, "MaskMetrology", (unsigned int) m_MaskMetrologyEnabled); + new MXmlNode(Node, "MaskMetrologyFileName", m_MaskMetrologyFile); new MXmlNode(Node, "UCSDOverride",(unsigned int) m_UCSDOverride); return Node; From beb5c322b6c6814afb4a2add8a47516e772a3afd Mon Sep 17 00:00:00 2001 From: ckierans Date: Tue, 2 Sep 2025 09:31:07 -0400 Subject: [PATCH 06/22] CHG: Updated x,y definition with mask rotation, added print statements --- src/MModuleDepthCalibration2024.cxx | 20 +++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/src/MModuleDepthCalibration2024.cxx b/src/MModuleDepthCalibration2024.cxx index 1f26bca6..d306c8d4 100644 --- a/src/MModuleDepthCalibration2024.cxx +++ b/src/MModuleDepthCalibration2024.cxx @@ -284,15 +284,21 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) double Xpos, Ypos, Zpos; if ( m_MaskMetrologyEnabled ) { // If we are applying the mask metrology correction, first define two new readout elements to help determine the intersection of these two strips + int HVStripID_flipped = 63 - HVStripID; + HVSH->SetStripID(HVStripID_flipped); MReadOutElementDoubleStrip R_LV = *dynamic_cast(LVSH->GetReadOutElement()); MReadOutElementDoubleStrip R_HV = *dynamic_cast(HVSH->GetReadOutElement()); vector inter = GetStripIntersection(R_LV, R_HV); Xpos = inter[0]; Ypos = inter[1]; - //Xpos = m_YPitches[DetID]*((m_NYStrips[DetID]/2.0) - ((double)LVStripID)); - //Ypos = m_XPitches[DetID]*((m_NXStrips[DetID]/2.0) - ((double)HVStripID)); - + + cout<<"StripLV: "< MModuleDepthCalibration2024::norm_pdf(vector x, double mu return result; } + +///////////////////////////////////////////////////////////////////////////////// + + bool MModuleDepthCalibration2024::LoadSplinesFile(MString FName) { // Input spline files should have the following format: @@ -662,9 +672,9 @@ vector MModuleDepthCalibration2024::GetStripIntersection(MReadOutElement // HVstrip is centered at (x,y,z) = (hv_strip_met[0], hv_strip_met[1], hv_strip_met[2]) // and is approximately parallel to the x axis, but rotated at angle hv_strip_met[3] // around the z axis of the detector - double x_intercept = (hv_strip_met[0]*tan(hv_strip_met[3]*TMath::DegToRad()) - lv_strip_met[0]/tan(lv_strip_met[3]*TMath::DegToRad()) - lv_strip_met[1] + hv_strip_met[1])/(tan(hv_strip_met[3]*TMath::DegToRad())-1/tan(lv_strip_met[3]*TMath::DegToRad())); + double x_intercept = (hv_strip_met[0]*tan(-hv_strip_met[5]*TMath::DegToRad()) - lv_strip_met[0]/tan(lv_strip_met[5]*TMath::DegToRad()) - lv_strip_met[1] + hv_strip_met[1])/(tan(-hv_strip_met[5]*TMath::DegToRad())-1/tan(lv_strip_met[5]*TMath::DegToRad())); - double y_intercept = (x_intercept - hv_strip_met[0])*tan(hv_strip_met[3]*TMath::DegToRad()) + hv_strip_met[1]; + double y_intercept = (x_intercept - hv_strip_met[0])*tan(hv_strip_met[5]*TMath::DegToRad()) + hv_strip_met[1]; return {x_intercept, y_intercept}; } From c138b17eec8e114368000d2b4fbe86dbf42003c2 Mon Sep 17 00:00:00 2001 From: ckierans Date: Mon, 27 Oct 2025 12:37:53 -0400 Subject: [PATCH 07/22] Updated metrology with correct understanding of rotation angles, now works for 2 detectors based on latest files from Aldo --- src/MModuleDepthCalibration2024.cxx | 38 +++++++++++++++-------------- 1 file changed, 20 insertions(+), 18 deletions(-) diff --git a/src/MModuleDepthCalibration2024.cxx b/src/MModuleDepthCalibration2024.cxx index d306c8d4..a0db0c83 100644 --- a/src/MModuleDepthCalibration2024.cxx +++ b/src/MModuleDepthCalibration2024.cxx @@ -219,7 +219,8 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) { if (Event->GetGuardRingVeto()==true) { - + //TODO: Handle events with GR vetos + Event->SetDepthCalibrationIncomplete(); return false; @@ -280,27 +281,28 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) int PixelCode = 10000*DetID + 100*LVStripID + HVStripID; // TODO: Calculate X and Y positions more rigorously using charge sharing. - // Somewhat confusing notation: HVStrips run parallel to X-axis, so we calculate X position with LVStrips. double Xpos, Ypos, Zpos; + if ( m_MaskMetrologyEnabled ) { // If we are applying the mask metrology correction, first define two new readout elements to help determine the intersection of these two strips + + + //TODO: Aldo's file has the HV strip ID flipped (as evident in the x,y position reported in the metrology files), so this is a placeholder until the files are fixed. int HVStripID_flipped = 63 - HVStripID; - HVSH->SetStripID(HVStripID_flipped); + MStripHit* HVSH_flipped = HVSH; + HVSH_flipped->SetStripID(HVStripID_flipped); MReadOutElementDoubleStrip R_LV = *dynamic_cast(LVSH->GetReadOutElement()); - MReadOutElementDoubleStrip R_HV = *dynamic_cast(HVSH->GetReadOutElement()); - + MReadOutElementDoubleStrip R_HV = *dynamic_cast(HVSH_flipped->GetReadOutElement()); + + //find the intercept of the two dominate strips based on the mask metrology vector inter = GetStripIntersection(R_LV, R_HV); Xpos = inter[0]; Ypos = inter[1]; - cout<<"StripLV: "< MModuleDepthCalibration2024::GetStripIntersection(MReadOutElement //Find the x position of two lines represented by the dominate strips: // LVstrip is centered at (x,y,z) = (lv_strip_met[0], lv_strip_met[1], lv_strip_met[2]) - // and is approximately parallel to the y axis, but rotated at angle lv_strip_met[3] + // and is approximately parallel to the y axis, but rotated at angle lv_strip_met[5] // around the z axis of the detector // HVstrip is centered at (x,y,z) = (hv_strip_met[0], hv_strip_met[1], hv_strip_met[2]) - // and is approximately parallel to the x axis, but rotated at angle hv_strip_met[3] + // and is approximately parallel to the x axis, but rotated at angle (hv_strip_met[5] - pi/2) // around the z axis of the detector - double x_intercept = (hv_strip_met[0]*tan(-hv_strip_met[5]*TMath::DegToRad()) - lv_strip_met[0]/tan(lv_strip_met[5]*TMath::DegToRad()) - lv_strip_met[1] + hv_strip_met[1])/(tan(-hv_strip_met[5]*TMath::DegToRad())-1/tan(lv_strip_met[5]*TMath::DegToRad())); + double x_intercept = (hv_strip_met[0]*tan((hv_strip_met[5]-90)*TMath::DegToRad()) - lv_strip_met[0]/tan(lv_strip_met[5]*TMath::DegToRad()) - lv_strip_met[1] + hv_strip_met[1])/(tan((hv_strip_met[5]-90)*TMath::DegToRad())-1/tan(lv_strip_met[5]*TMath::DegToRad())); - double y_intercept = (x_intercept - hv_strip_met[0])*tan(hv_strip_met[5]*TMath::DegToRad()) + hv_strip_met[1]; + double y_intercept = (x_intercept - hv_strip_met[0])*tan((hv_strip_met[5]-90)*TMath::DegToRad()) + hv_strip_met[1]; return {x_intercept, y_intercept}; } From 9410061c264edf45c67ccbab882ceb4a28e37166 Mon Sep 17 00:00:00 2001 From: julianmgerber <95662611+julianmgerber@users.noreply.github.com> Date: Tue, 18 Nov 2025 10:00:39 -0800 Subject: [PATCH 08/22] Strip Pairing Chi Square Algorithm Updated (#61) * Merging select files from StripPairingTesting into StripPairingCleanUp so I can do pull request * Clean up comments. Remove long debugging statements for readability * Calling charge trapping correction function in chi square calculation and hit population * Cleaning up comments * Reformatting strip pairing module * Update copywrite. Change name of module. * Fixing documentation in header files * Change how strip pairing chi square is read out to .dat file * Adding comments for all the nested vectors * Adding new version to StreamDat that reads out LV and HV energy of each hit * Adding Sean's new strip pairing expos to chi square version of strip pairing * Changing all instances of RedChiSquare to ReducedChiSquare for clarity * Flags for multiple hits on a single strip are now defined on the MHit level * Initializing m_StripHitMultipleTimesX and m_StripHitMultipleTimesY to false * Changing StreamDAT version from Version 1 to Version 3 * Rebased and fixed merge conflict with strip pairing --------- Co-authored-by: Julian Gerber --- include/MReadOutAssembly.h | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/include/MReadOutAssembly.h b/include/MReadOutAssembly.h index 1a8b0aa6..23db7b7c 100644 --- a/include/MReadOutAssembly.h +++ b/include/MReadOutAssembly.h @@ -454,8 +454,7 @@ class MReadOutAssembly : public MReadOutSequence bool m_StripHitBelowThreshold; MString m_StripHitBelowThresholdString; - - + //! True if event has been filtered out bool m_FilteredOut; From f858369c9e908fd82aeb68913e8a0cc18334169b Mon Sep 17 00:00:00 2001 From: ckierans Date: Tue, 9 Dec 2025 23:01:53 -0500 Subject: [PATCH 09/22] Resolving merge conflict with Sean's recent PR --- src/MModuleDepthCalibration2024.cxx | 16 +++++----------- 1 file changed, 5 insertions(+), 11 deletions(-) diff --git a/src/MModuleDepthCalibration2024.cxx b/src/MModuleDepthCalibration2024.cxx index a0db0c83..e2dd9da1 100644 --- a/src/MModuleDepthCalibration2024.cxx +++ b/src/MModuleDepthCalibration2024.cxx @@ -61,7 +61,6 @@ MModuleDepthCalibration2024::MModuleDepthCalibration2024() : MModule() m_XmlTag = "DepthCalibration2024"; // Set all modules, which have to be done before this module - AddPreceedingModuleType(MAssembly::c_EventLoader, true); AddPreceedingModuleType(MAssembly::c_EnergyCalibration, true); AddPreceedingModuleType(MAssembly::c_StripPairing, true); AddPreceedingModuleType(MAssembly::c_TACcut, true); @@ -285,14 +284,8 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) if ( m_MaskMetrologyEnabled ) { // If we are applying the mask metrology correction, first define two new readout elements to help determine the intersection of these two strips - - - //TODO: Aldo's file has the HV strip ID flipped (as evident in the x,y position reported in the metrology files), so this is a placeholder until the files are fixed. - int HVStripID_flipped = 63 - HVStripID; - MStripHit* HVSH_flipped = HVSH; - HVSH_flipped->SetStripID(HVStripID_flipped); MReadOutElementDoubleStrip R_LV = *dynamic_cast(LVSH->GetReadOutElement()); - MReadOutElementDoubleStrip R_HV = *dynamic_cast(HVSH_flipped->GetReadOutElement()); + MReadOutElementDoubleStrip R_HV = *dynamic_cast(HVSH->GetReadOutElement()); //find the intercept of the two dominate strips based on the mask metrology vector inter = GetStripIntersection(R_LV, R_HV); @@ -300,11 +293,11 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) Ypos = inter[1]; Zpos = 0.0; + } else { //If no metrology is used, define the strip positions based on the detector pitch and number of strip hits - - Xpos = m_YPitches[DetID]*((m_NYStrips[DetID]/2.0) - ((double)LVStripID)); - Ypos = m_XPitches[DetID]*((m_NXStrips[DetID]/2.0) - ((double)HVStripID)); + Xpos = m_YPitches[DetID]*(((m_NYStrips[DetID]-1)/2.0) - ((double)LVStripID)); + Ypos = m_XPitches[DetID]*((double)HVStripID - ((m_NXStrips[DetID]-1)/2.0)); Zpos = 0.0; } @@ -390,6 +383,7 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) Zsigma = sqrt(depth_var/prob_sum); Zpos = mean_depth; + //Zpos = mean_depth - (m_Thicknesses[DetID]/2.0); // Add the depth to the GUI histogram. if (Event->IsStripPairingIncomplete()==false) { From 108275f9c9fc273333026c67d8b18acc6ac9ab96 Mon Sep 17 00:00:00 2001 From: ckierans Date: Tue, 9 Dec 2025 23:19:44 -0500 Subject: [PATCH 10/22] Resolving merge conflict with Sean's recent PR --- include/MGUIOptionsDepthCalibration2024.h | 2 +- include/MModuleDepthCalibration2024.h | 2 +- src/MGUIOptionsDepthCalibration2024.cxx | 8 +- src/MModuleDepthCalibration2024.cxx | 312 +++++++++++----------- 4 files changed, 169 insertions(+), 155 deletions(-) diff --git a/include/MGUIOptionsDepthCalibration2024.h b/include/MGUIOptionsDepthCalibration2024.h index 6c14df24..6641fddf 100644 --- a/include/MGUIOptionsDepthCalibration2024.h +++ b/include/MGUIOptionsDepthCalibration2024.h @@ -77,7 +77,7 @@ class MGUIOptionsDepthCalibration2024 : public MGUIOptions //! Select mask metrology file to load. This gives the translation and rotation for each strip in the detector frame MGUIEFileSelector* m_MaskMetrologyFileSelector; - int m_UseMaskMetCorr; + bool m_UseMaskMetCorr; TGCheckButton* m_MaskMetModeCB; enum ButtonIDs {c_MetFile}; diff --git a/include/MModuleDepthCalibration2024.h b/include/MModuleDepthCalibration2024.h index 6a697761..b2794d47 100644 --- a/include/MModuleDepthCalibration2024.h +++ b/include/MModuleDepthCalibration2024.h @@ -1,4 +1,4 @@ -/*:q +/* * MModuleDepthCalibration2024.h * * Copyright (C) 2008-2008 by Andreas Zoglauer. diff --git a/src/MGUIOptionsDepthCalibration2024.cxx b/src/MGUIOptionsDepthCalibration2024.cxx index 9ead8652..569ce96b 100644 --- a/src/MGUIOptionsDepthCalibration2024.cxx +++ b/src/MGUIOptionsDepthCalibration2024.cxx @@ -88,11 +88,11 @@ void MGUIOptionsDepthCalibration2024::Create() // TGLayoutHints* FileLabelLayout = new TGLayoutHints(kLHintsTop | kLHintsExpandX, m_FontScaler*65 + 21*m_FontScaler, m_FontScaler*65, 0, 2*m_FontScaler); m_MaskMetrologyFileSelector = new MGUIEFileSelector(m_OptionsFrame, "", dynamic_cast(m_Module)->GetMaskMetrologyFileName()); - m_MaskMetrologyFileSelector->SetFileType("metrology", "*.csv"); + m_MaskMetrologyFileSelector->SetFileType("metrology", "*.metrology.csv"); m_OptionsFrame->AddFrame(m_MaskMetrologyFileSelector, LabelLayout); - if (m_UseMaskMetCorr) { + if (m_UseMaskMetCorr == true) { m_MaskMetrologyFileSelector->SetEnabled(true); } else { m_MaskMetrologyFileSelector->SetEnabled(false); @@ -127,10 +127,10 @@ bool MGUIOptionsDepthCalibration2024::ProcessMessage(long Message, long Paramete switch (Parameter1) { case c_MetFile: if (m_MaskMetModeCB->GetState() == kButtonDown) { - m_UseMaskMetCorr = 1; + m_UseMaskMetCorr = true; m_MaskMetrologyFileSelector->SetEnabled(true); } else if (m_MaskMetModeCB->GetState() == kButtonUp) { - m_UseMaskMetCorr = 0; + m_UseMaskMetCorr = false; m_MaskMetrologyFileSelector->SetEnabled(false); } break; diff --git a/src/MModuleDepthCalibration2024.cxx b/src/MModuleDepthCalibration2024.cxx index e2dd9da1..ece42791 100644 --- a/src/MModuleDepthCalibration2024.cxx +++ b/src/MModuleDepthCalibration2024.cxx @@ -118,11 +118,11 @@ bool MModuleDepthCalibration2024::Initialize() vector DetList = m_Geometry->GetDetectorList(); // Look through the Geometry and get the names and thicknesses of all the detectors. - for(unsigned int i = 0; i < DetList.size(); ++i){ + for (unsigned int i = 0; i < DetList.size(); ++i) { // For now, DetID is in order of detectors, which puts contraints on how the geometry file should be written. // If using the card cage at UCSD, default to DetID=11. unsigned int DetID = i; - if ( m_UCSDOverride ){ + if (m_UCSDOverride == true) { DetID = 11; } @@ -134,7 +134,7 @@ bool MModuleDepthCalibration2024::Initialize() string det_name = vol->GetName().GetString(); if (find(DetectorNames.begin(), DetectorNames.end(), det_name) == DetectorNames.end()) { DetectorNames.push_back(det_name); - m_Thicknesses[DetID] = 2*(det->GetStructuralSize().GetZ()); + m_Thicknesses[DetID] = 2 * (det->GetStructuralSize().GetZ()); MDStrip3D* strip = dynamic_cast(det); m_XPitches[DetID] = strip->GetPitchX(); m_YPitches[DetID] = strip->GetPitchY(); @@ -202,10 +202,10 @@ void MModuleDepthCalibration2024::CreateExpos() // Set the histogram display m_ExpoDepthCalibration = new MGUIExpoDepthCalibration2024(this); m_ExpoDepthCalibration->SetDepthHistogramArrangement(&m_DetectorIDs); - for(unsigned int i = 0; i < m_DetectorIDs.size(); ++i){ + 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.,thickness/2.); + m_ExpoDepthCalibration->SetDepthHistogramParameters(DetID, 120, -thickness/2.0,thickness/2.0); } m_Expos.push_back(m_ExpoDepthCalibration); } @@ -217,7 +217,7 @@ void MModuleDepthCalibration2024::CreateExpos() bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) { - if (Event->GetGuardRingVeto()==true) { + if (Event->GetGuardRingVeto() == true) { //TODO: Handle events with GR vetos Event->SetDepthCalibrationIncomplete(); @@ -225,7 +225,7 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) } else { - for( unsigned int i = 0; i < Event->GetNHits(); ++i ){ + for (unsigned int i = 0; i < Event->GetNHits(); ++i ){ // Each event represents one photon. It contains Hits, representing interaction sites. // H is a pointer to an instance of the MHit class. Each Hit has activated strips, represented by // instances of the MStripHit class. @@ -235,7 +235,7 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) // Handle different grades differently // GRADE=-1 is an error. Break from the loop and continue. - if ( Grade < 0 ){ + if (Grade < 0){ H->SetNoDepth(); Event->SetDepthCalibrationIncomplete(); if (Grade == -1) { @@ -255,15 +255,14 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) } } else { // If the Grade is 0-4, we can handle it. - MVector LocalPosition, PositionResolution, GlobalPosition, GlobalResolution, LocalOrigin; - // Calculate the position. If error is thrown, record and no depth. // Take a Hit and separate its activated X- and Y-strips into separate vectors. vector LVStrips; vector HVStrips; - for( unsigned int j = 0; j < H->GetNStripHits(); ++j){ + + for (unsigned int j = 0; j < H->GetNStripHits(); ++j){ MStripHit* SH = H->GetStripHit(j); - if( SH->IsLowVoltageStrip() ) LVStrips.push_back(SH); else HVStrips.push_back(SH); + if (SH->IsLowVoltageStrip()) LVStrips.push_back(SH); else HVStrips.push_back(SH); } double LVEnergyFraction; @@ -279,28 +278,29 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) int HVStripID = HVSH->GetStripID(); int PixelCode = 10000*DetID + 100*LVStripID + HVStripID; - // TODO: Calculate X and Y positions more rigorously using charge sharing. - double Xpos, Ypos, Zpos; + //Define the X/Y positions based on the detector pitch and number of strip hits + // LV strip 0 is in -ve X direction, HV strip 0 is in -ve Y direction. + // Confusingly, the strips parallel to the Y axis determines the X position, and the "X strips" determine the Y position + double Xpos = m_YPitches[DetID]*((double)LVStripID - ((m_NYStrips[DetID]-1)/2.0)); + double Ypos = m_XPitches[DetID]*((double)HVStripID - ((m_NXStrips[DetID]-1)/2.0)); + double Zpos = 0.0; - if ( m_MaskMetrologyEnabled ) { + + // TODO: Confirm X and Y implementation below with Aldo's new metrology files. His old files swapped these. + if (m_MaskMetrologyEnabled == true) { // If we are applying the mask metrology correction, first define two new readout elements to help determine the intersection of these two strips MReadOutElementDoubleStrip R_LV = *dynamic_cast(LVSH->GetReadOutElement()); MReadOutElementDoubleStrip R_HV = *dynamic_cast(HVSH->GetReadOutElement()); - //find the intercept of the two dominate strips based on the mask metrology + //find the intercept of the two dominate strips based on the mask metrology, and update Xpos and Ypos vector inter = GetStripIntersection(R_LV, R_HV); Xpos = inter[0]; Ypos = inter[1]; - Zpos = 0.0; + } - } else { - //If no metrology is used, define the strip positions based on the detector pitch and number of strip hits - Xpos = m_YPitches[DetID]*(((m_NYStrips[DetID]-1)/2.0) - ((double)LVStripID)); - Ypos = m_XPitches[DetID]*((double)HVStripID - ((m_NXStrips[DetID]-1)/2.0)); - Zpos = 0.0; - } + // TODO: Calculate X and Y positions more rigorously using charge sharing. double Xsigma = m_YPitches[DetID]/sqrt(12.0); double Ysigma = m_XPitches[DetID]/sqrt(12.0); @@ -316,11 +316,12 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) double HVTiming = HVSH->GetTiming(); // If there aren't coefficients loaded, then calibration is incomplete. - if( Coeffs == nullptr ){ + if (Coeffs == nullptr){ //set the bad flag for depth H->SetNoDepth(); Event->SetDepthCalibrationIncomplete(); ++m_Error1; + } else if (CTDVec.size() == 0) { cout << "Empty CTD vector" << endl; H->SetNoDepth(); @@ -336,7 +337,6 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) } else { // If there are coefficients and timing information is loaded, try calculating the CTD and depth - double CTD = (HVTiming - LVTiming); // Confirmed that this matches SP's python code. @@ -348,7 +348,7 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) double noise = GetTimingNoiseFWHM(PixelCode, H->GetEnergy()); //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)) ){ + if ((CTD_s < (Xmin - 2.0*noise)) || (CTD_s > (Xmax + 2.0*noise))) { H->SetNoDepth(); Event->SetDepthCalibrationIncomplete(); ++m_Error2; @@ -365,20 +365,23 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) // Weight the depth by probability double prob_sum = 0.0; - for( unsigned int k=0; k < prob_dist.size(); ++k ){ + 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]; + + 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; kGetSensitiveVolume(0)->GetPositionInWorldVolume(LocalPosition); + MVector LocalPosition(Xpos, Ypos, Zpos); + MVector LocalOrigin(0.0, 0.0, 0.0); + MVector GlobalPosition = m_Detectors[DetID]->GetSensitiveVolume(0)->GetPositionInWorldVolume(LocalPosition); // Make sure XYZ resolution are correctly mapped to the global coord system. - PositionResolution.SetXYZ(Xsigma, Ysigma, Zsigma); - GlobalResolution = ((m_Detectors[DetID]->GetSensitiveVolume(0)->GetPositionInWorldVolume(PositionResolution)) - (m_Detectors[DetID]->GetSensitiveVolume(0)->GetPositionInWorldVolume(LocalOrigin))).Abs(); + MVector PositionResolution(Xsigma, Ysigma, Zsigma); + MVector GlobalResolution = ((m_Detectors[DetID]->GetSensitiveVolume(0)->GetPositionInWorldVolume(PositionResolution)) - (m_Detectors[DetID]->GetSensitiveVolume(0)->GetPositionInWorldVolume(LocalOrigin))).Abs(); H->SetPosition(GlobalPosition); @@ -458,13 +461,12 @@ double MModuleDepthCalibration2024::GetTimingNoiseFWHM(int PixelCode, double Ene // Should follow 1/E relation // TODO: Determine real energy dependence and implement it here. double noiseFWHM = 0.0; - if ( m_Coeffs_Energy != 0 ){ + if (m_Coeffs_Energy != 0) { noiseFWHM = m_Coeffs[PixelCode][2] * m_Coeffs_Energy/Energy; - if ( noiseFWHM < 3.0*2.355 ){ + if (noiseFWHM < 3.0*2.355) { noiseFWHM = 3.0*2.355; } - } - else { + } else { noiseFWHM = 6.0*2.355; } return noiseFWHM; @@ -474,42 +476,43 @@ double MModuleDepthCalibration2024::GetTimingNoiseFWHM(int PixelCode, double Ene ///////////////////////////////////////////////////////////////////////////////// -bool MModuleDepthCalibration2024::LoadCoeffsFile(MString FName) +bool MModuleDepthCalibration2024::LoadCoeffsFile(MString FileName) { // Read in the stretch and offset file, which should have a header line with information on the measurements: // ### 800 V 80 K 59.5 keV // And which should contain for each pixel: // Pixel code (10000*det + 100*LVStrip + HVStrip), Stretch, Offset, Timing/CTD noise, Chi2 for the CTD fit (for diagnostics mainly) - MFile F; - if( F.Open(FName) == false ){ + + MFile CoeffsFile; + if (CoeffsFile.Open(FileName) == false) { cout << "ERROR in MModuleDepthCalibration2024::LoadCoeffsFile: failed to open coefficients file." << endl; return false; - } else { - MString Line; - while( F.ReadLine( Line ) ){ - if ( Line.BeginsWith('#') ){ - std::vector Tokens = Line.Tokenize(" "); - m_Coeffs_Energy = Tokens[5].ToDouble(); - cout << "The stretch and offset were calculated for " << m_Coeffs_Energy << " keV." << endl; - } - else { - std::vector Tokens = Line.Tokenize(","); - if( Tokens.size() == 5 ){ - int PixelCode = Tokens[0].ToInt(); - double Stretch = Tokens[1].ToDouble(); - double Offset = Tokens[2].ToDouble(); - double CTD_FWHM = Tokens[3].ToDouble() * 2.355; - double Chi2 = Tokens[4].ToDouble(); - // Previous iteration of depth calibration read in "Scale" instead of ctd resolution. - vector coeffs; - coeffs.push_back(Stretch); coeffs.push_back(Offset); coeffs.push_back(CTD_FWHM); coeffs.push_back(Chi2); - m_Coeffs[PixelCode] = coeffs; - } + } + + MString Line; + while (CoeffsFile.ReadLine(Line)) { + if (Line.BeginsWith('#') == true) { + std::vector Tokens = Line.Tokenize(" "); + m_Coeffs_Energy = Tokens[5].ToDouble(); + cout << "The stretch and offset were calculated for " << m_Coeffs_Energy << " keV." << endl; + } else { + std::vector Tokens = Line.Tokenize(","); + if (Tokens.size() == 5) { + int PixelCode = Tokens[0].ToInt(); + double Stretch = Tokens[1].ToDouble(); + double Offset = Tokens[2].ToDouble(); + double CTD_FWHM = Tokens[3].ToDouble() * 2.355; + double Chi2 = Tokens[4].ToDouble(); + // Previous iteration of depth calibration read in "Scale" instead of ctd resolution. + vector coeffs; + coeffs.push_back(Stretch); coeffs.push_back(Offset); coeffs.push_back(CTD_FWHM); coeffs.push_back(Chi2); + m_Coeffs[PixelCode] = coeffs; } } - F.Close(); } + CoeffsFile.Close(); + return true; } @@ -521,11 +524,11 @@ bool MModuleDepthCalibration2024::LoadCoeffsFile(MString FName) std::vector* MModuleDepthCalibration2024::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. - if( m_CoeffsFileIsLoaded ){ - if( m_Coeffs.count(PixelCode) > 0 ){ + if (m_CoeffsFileIsLoaded == true) { + if (m_Coeffs.count(PixelCode) > 0) { return &m_Coeffs[PixelCode]; } else { - // cout << "MModuleDepthCalibration2024::GetPixelCoeffs: cannot get stretch and offset; pixel code " << PixelCode << " not found." << endl; + cout << "MModuleDepthCalibration2024::GetPixelCoeffs: cannot get stretch and offset; pixel code " << PixelCode << " not found." << endl; return nullptr; } } else { @@ -542,7 +545,7 @@ std::vector* MModuleDepthCalibration2024::GetPixelCoeffs(int PixelCode) vector MModuleDepthCalibration2024::norm_pdf(vector x, double mu, double sigma) { vector result; - for( unsigned int i=0; i MModuleDepthCalibration2024::norm_pdf(vector x, double mu ///////////////////////////////////////////////////////////////////////////////// -bool MModuleDepthCalibration2024::LoadSplinesFile(MString FName) +bool MModuleDepthCalibration2024::LoadSplinesFile(MString FileName) { // Input spline files should have the following format: // ### DetID, HV, Temperature, Photopeak Energy (TODO: More? Fewer?) // depth, ctd0, ctd1, ctd2.... (Basically, allow for CTDs for different subpixel regions) // '' '' '' - MFile F; - if( F.Open(FName) == false ){ + + MFile SplineFile; + if (SplineFile.Open(FileName) == false) { cout << "ERROR in MModuleDepthCalibration2024::LoadSplinesFile: failed to open splines file." << endl; return false; } + vector DepthVec; vector> CTDArr; + bool Result = true; - MString line; + MString Line; int DetID = 0; - while (F.ReadLine(line)) { - if (line.Length() != 0) { - if (line.BeginsWith("#")) { + while (SplineFile.ReadLine(Line)) { + if (Line.Length() != 0) { + if (Line.BeginsWith("#") == true) { // If we've reached a new ctd spline then record the previous one in the m_SplineMaps and start a new one. - vector tokens = line.Tokenize(" "); + vector tokens = Line.Tokenize(" "); + if (DepthVec.size() > 0) { Result &= AddDepthCTD(DepthVec, CTDArr, DetID, m_DepthGrid, m_CTDMap, m_SplineMap, 1000); } + DepthVec.clear(); CTDArr.clear(); DetID = tokens[1].ToInt(); + } else { - vector tokens = line.Tokenize(","); + + vector tokens = Line.Tokenize(","); DepthVec.push_back(tokens[0].ToDouble()); + // Multiple CTDs allowed. for (unsigned int i = 0; i < (tokens.size() - 1); ++i) { while (i>=CTDArr.size()) { @@ -593,6 +604,9 @@ bool MModuleDepthCalibration2024::LoadSplinesFile(MString FName) } } } + + SplineFile.Close(); + //make last spline if (DepthVec.size() > 0) { Result &= AddDepthCTD(DepthVec, CTDArr, DetID, m_DepthGrid, m_CTDMap, m_SplineMap, 1000); @@ -606,47 +620,48 @@ bool MModuleDepthCalibration2024::LoadSplinesFile(MString FName) ///////////////////////////////////////////////////////////////////////////////// -bool MModuleDepthCalibration2024::LoadMaskMetrologyFile(MString FName) +bool MModuleDepthCalibration2024::LoadMaskMetrologyFile(MString FileName) { //Read the Mask Metrology File // Det ID, Side (l,h), Strip ID (0-63), x_mm, y_mm, z_mm, roll_deg, pitch_deg, yaw_deg - MFile F; - if( F.Open(FName) == false ){ + MFile MetrologyFile; + if (MetrologyFile.Open(FileName) == false) { cout << "ERROR in MModuleDepthCalibration2024::LoadMaskMetrologyFile: failed to open metrology file." << endl; return false; - } else { - MString Line; - while( F.ReadLine( Line ) ){ - if ( Line.BeginsWith('#') ){ - } - else { - std::vector Tokens = Line.Tokenize(","); - if( Tokens.size() == 9 ){ - //Define the readout element to track det ID, strip ID, and lv/hv - MReadOutElementDoubleStrip R; - R.SetDetectorID(Tokens[0].ToInt()); - R.IsLowVoltageStrip((Tokens[1].ToString() == "p") || (Tokens[1].ToString() == "l")); - R.SetStripID(Tokens[2].ToInt()); - double Strip_MetX = Tokens[3].ToDouble()/10; //convert to cm - double Strip_MetY = Tokens[4].ToDouble()/10; //convert to cm - double Strip_MetZ = Tokens[5].ToDouble()/10; //convert to cm - double Strip_Roll = Tokens[6].ToDouble(); - double Strip_Pitch = Tokens[7].ToDouble(); - double Strip_Yaw = Tokens[8].ToDouble(); - vector maskmet; - maskmet.push_back(Strip_MetX); maskmet.push_back(Strip_MetY); maskmet.push_back(Strip_MetZ); - maskmet.push_back(Strip_Roll); maskmet.push_back(Strip_Pitch); maskmet.push_back(Strip_Yaw); - - //make the map that defines the metrology info for each readout element - m_MaskMetrology[R] = maskmet; - } + } + + MString Line; + while (MetrologyFile.ReadLine(Line)) { + if (Line.BeginsWith('#') == true) continue; + else { + std::vector Tokens = Line.Tokenize(","); + if (Tokens.size() == 9) { + //Define the readout element to track det ID, strip ID, and lv/hv + MReadOutElementDoubleStrip R; + R.SetDetectorID(Tokens[0].ToInt()); + R.IsLowVoltageStrip((Tokens[1].ToString() == "p") || (Tokens[1].ToString() == "l")); + R.SetStripID(Tokens[2].ToInt()); + double Strip_MetX = Tokens[3].ToDouble()/10; //convert to cm + double Strip_MetY = Tokens[4].ToDouble()/10; //convert to cm + double Strip_MetZ = Tokens[5].ToDouble()/10; //convert to cm + double Strip_Roll = Tokens[6].ToDouble(); + double Strip_Pitch = Tokens[7].ToDouble(); + double Strip_Yaw = Tokens[8].ToDouble(); + vector maskmet; + maskmet.push_back(Strip_MetX); maskmet.push_back(Strip_MetY); maskmet.push_back(Strip_MetZ); + maskmet.push_back(Strip_Roll); maskmet.push_back(Strip_Pitch); maskmet.push_back(Strip_Yaw); + + //make the map that defines the metrology info for each readout element + m_MaskMetrology[R] = maskmet; + } else { + cout << "ERROR in MModuleDepthCalibration2024::LoadMaskMetrologyFile: incorrect number of tokens in the file." << endl; + return false; } } - - F.Close(); - } + MetrologyFile.Close(); + return true; } @@ -658,8 +673,8 @@ bool MModuleDepthCalibration2024::LoadMaskMetrologyFile(MString FName) vector MModuleDepthCalibration2024::GetStripIntersection(MReadOutElementDoubleStrip R_LVStrip, MReadOutElementDoubleStrip R_HVStrip){ // Function to find the intersection between two strips based on the Mask Metrology - vector lv_strip_met = m_MaskMetrology[R_LVStrip]; - vector hv_strip_met = m_MaskMetrology[R_HVStrip]; + vector LVStripMet = m_MaskMetrology[R_LVStrip]; + vector HVStripMet = m_MaskMetrology[R_HVStrip]; //Find the x position of two lines represented by the dominate strips: // LVstrip is centered at (x,y,z) = (lv_strip_met[0], lv_strip_met[1], lv_strip_met[2]) @@ -668,11 +683,11 @@ vector MModuleDepthCalibration2024::GetStripIntersection(MReadOutElement // HVstrip is centered at (x,y,z) = (hv_strip_met[0], hv_strip_met[1], hv_strip_met[2]) // and is approximately parallel to the x axis, but rotated at angle (hv_strip_met[5] - pi/2) // around the z axis of the detector - double x_intercept = (hv_strip_met[0]*tan((hv_strip_met[5]-90)*TMath::DegToRad()) - lv_strip_met[0]/tan(lv_strip_met[5]*TMath::DegToRad()) - lv_strip_met[1] + hv_strip_met[1])/(tan((hv_strip_met[5]-90)*TMath::DegToRad())-1/tan(lv_strip_met[5]*TMath::DegToRad())); + double XIntercept = (HVStripMet[0]*tan((HVStripMet[5]-90)*TMath::DegToRad()) - LVStripMet[0]/tan(LVStripMet[5]*TMath::DegToRad()) - LVStripMet[1] + HVStripMet[1])/(tan((HVStripMet[5]-90)*TMath::DegToRad())-1/tan(LVStripMet[5]*TMath::DegToRad())); - double y_intercept = (x_intercept - hv_strip_met[0])*tan((hv_strip_met[5]-90)*TMath::DegToRad()) + hv_strip_met[1]; + double YIntercept = (XIntercept - HVStripMet[0])*tan((HVStripMet[5]-90)*TMath::DegToRad()) + HVStripMet[1]; - return {x_intercept, y_intercept}; + return {XIntercept, YIntercept}; } @@ -685,10 +700,10 @@ int MModuleDepthCalibration2024::GetHitGrade(MHit* H){ // 5 possible grades, and one Error Grade, -1. GRADE 4 is as yet uncategorized complicated geometry. GRADE 5 means multiple, presumably separated strip hits. //organize x and y strips into vectors - if( H == NULL ){ + if (H == nullptr) { return -1; } - if( H->GetNStripHits() == 0 ){ + if (H->GetNStripHits() == 0) { // Error if no strip hits listed. Bad grade is returned cout << "ERROR in MModuleDepthCalibration2024: HIT WITH NO STRIP HITS" << endl; return -1; @@ -699,11 +714,11 @@ int MModuleDepthCalibration2024::GetHitGrade(MHit* H){ std::vector HVStrips; vector LVStripIDs; vector HVStripIDs; - for( unsigned int j = 0; j < H->GetNStripHits(); ++j){ + for (unsigned int j = 0; j < H->GetNStripHits(); ++j) { MStripHit* SH = H->GetStripHit(j); - if( SH == NULL ) { cout << "ERROR in MModuleDepthCalibration2024: Depth Calibration: got NULL strip hit :( " << endl; return -1;} - if( SH->GetEnergy() == 0 ) { cout << "ERROR in MModuleDepthCalibration2024: Depth Calibration: got strip without energy :( " << endl; return -1;} - if( SH->IsLowVoltageStrip() ){ + if (SH == nullptr ) { cout << "ERROR in MModuleDepthCalibration2024: Depth Calibration: got NULL strip hit :( " << endl; return -1;} + if (SH->GetEnergy() == 0 ) { cout << "ERROR in MModuleDepthCalibration2024: Depth Calibration: got strip without energy :( " << endl; return -1;} + if (SH->IsLowVoltageStrip()) { LVStrips.push_back(SH); LVStripIDs.push_back(SH->GetStripID()); } @@ -716,11 +731,11 @@ int MModuleDepthCalibration2024::GetHitGrade(MHit* H){ // if the same strip has multiple hits, this is a bad grade. bool MultiHitX = H->GetStripHitMultipleTimesX(); bool MultiHitY = H->GetStripHitMultipleTimesY(); - if( MultiHitX || MultiHitY ){ + if (MultiHitX || MultiHitY) { return 5; } - if( LVStrips.size()>0 && HVStrips.size()>0 ){ + if (LVStrips.size()>0 && HVStrips.size()>0) { int HVmin = * std::min_element(HVStripIDs.begin(), HVStripIDs.end()); int HVmax = * std::max_element(HVStripIDs.begin(), HVStripIDs.end()); @@ -728,7 +743,7 @@ int MModuleDepthCalibration2024::GetHitGrade(MHit* H){ int LVmax = * std::max_element(LVStripIDs.begin(), LVStripIDs.end()); // If the strip hits are not all adjacent, it's a bad grade. - if ( ((HVmax - HVmin) >= (HVStrips.size())) || ((LVmax - LVmin) >= (LVStrips.size())) ){ + if ( ((HVmax - HVmin) >= (HVStrips.size())) || ((LVmax - LVmin) >= (LVStrips.size())) ) { return 6; } } @@ -739,48 +754,48 @@ int MModuleDepthCalibration2024::GetHitGrade(MHit* H){ int return_value; // If 1 strip on each side, GRADE=0 // This represents the center of the pixel - if( ((LVStrips.size() == 1) && (HVStrips.size() == 1)) || ((LVStrips.size() == 3) && (HVStrips.size() == 3)) ){ + if ( ((LVStrips.size() == 1) && (HVStrips.size() == 1)) || ((LVStrips.size() == 3) && (HVStrips.size() == 3)) ) { return_value = 0; } // If 2 hits on N side and 1 on P, GRADE=1 // This represents the middle of the edges of the pixel - else if( (LVStrips.size() == 1) && (HVStrips.size() == 2) ){ + else if ( (LVStrips.size() == 1) && (HVStrips.size() == 2) ) { return_value = 1; } // If 2 hits on P and 1 on N, GRADE=2 // This represents the middle of the edges of the pixel - else if( (LVStrips.size() == 2) && (HVStrips.size() == 1) ){ + else if ( (LVStrips.size() == 2) && (HVStrips.size() == 1) ) { return_value = 2; } // If 2 strip hits on both sides, GRADE=3 // This represents the corners the pixel - else if( (LVStrips.size() == 2) && (HVStrips.size() == 2) ){ + else if ( (LVStrips.size() == 2) && (HVStrips.size() == 2) ) { return_value = 3; } // If 3 hits on N side and 1 on P, GRADE=0 // This represents the middle of the pixel, near the p (LV) side of the detector. - else if( (LVStrips.size() == 1) && (HVStrips.size() == 3) ){ + else if ( (LVStrips.size() == 1) && (HVStrips.size() == 3) ) { return_value = 0; } // If 3 hits on P and 1 on N, GRADE=0 // This represents the middle of the pixel, near the n (HV) side of the detector. - else if( (LVStrips.size() == 3) && (HVStrips.size() == 1) ){ + else if ( (LVStrips.size() == 3) && (HVStrips.size() == 1) ) { return_value = 0; } // If 3 hits on N side and 2 on P, GRADE=0 // This represents the middle of the edge of the pixel, near the p (LV) side of the detector. - else if( (LVStrips.size() == 2) && (HVStrips.size() == 3) ){ + else if ( (LVStrips.size() == 2) && (HVStrips.size() == 3) ) { return_value = 2; } // If 3 hits on P and 2 on N, GRADE=0 // This represents the middle of the edge of the pixel, near the n (HV) side of the detector. - else if( (LVStrips.size() == 3) && (HVStrips.size() == 2) ){ + else if ( (LVStrips.size() == 3) && (HVStrips.size() == 2) ) { return_value = 1; } @@ -901,17 +916,16 @@ vector MModuleDepthCalibration2024::GetCTD(int DetID, int Grade) { // Retrieves the appropriate CTD vector given the Detector ID and Event Grade passed - if( !m_SplinesFileIsLoaded ){ + if (!m_SplinesFileIsLoaded) { cout << "MModuleDepthCalibration2024::GetCTD: cannot return Depth to CTD relation because the file was not loaded." << endl; return vector (); } // If there is a CTD array for the given detector, return it. // If the Grade is larger than the number of CTD vectors stored, then just return Grade 0 vector. - if( m_CTDMap.count(DetID) > 0 ){ + if (m_CTDMap.count(DetID) > 0) { if ( ((int)m_CTDMap[DetID].size()) > Grade) { return (m_CTDMap[DetID][Grade]); - } - else { + } else { return (m_CTDMap[DetID][0]); } } else { @@ -928,14 +942,14 @@ vector MModuleDepthCalibration2024::GetDepth(int DetID) { // Retrieves the appropriate CTD vector given the Detector ID and Event Grade passed - if( !m_SplinesFileIsLoaded ){ + if (!m_SplinesFileIsLoaded) { cout << "MModuleDepthCalibration2024::GetDepth: cannot return Depth grid because the file was not loaded." << endl; return vector (); } // If there is a CTD array for the given detector, return it. // If the Grade is larger than the number of CTD vectors stored, then just return Grade 0 vector. - if( m_DepthGrid.count(DetID) > 0 ){ + if (m_DepthGrid.count(DetID) > 0){ return m_DepthGrid[DetID]; } else { cout << "MModuleDepthCalibration2024::GetDepth: No Depth grid is loaded for Det " << DetID << "." << endl; @@ -992,28 +1006,28 @@ bool MModuleDepthCalibration2024::ReadXmlConfiguration(MXmlNode* Node) //! Read the configuration data from an XML node MXmlNode* CoeffsFileNameNode = Node->GetNode("CoeffsFileName"); - if (CoeffsFileNameNode != 0) { + if (CoeffsFileNameNode != nullptr) { m_CoeffsFile = CoeffsFileNameNode->GetValue(); } MXmlNode* SplinesFileNameNode = Node->GetNode("SplinesFileName"); - if (SplinesFileNameNode != 0) { + if (SplinesFileNameNode != nullptr) { m_SplinesFile = SplinesFileNameNode->GetValue(); } MXmlNode* MasKMetrologyNode = Node->GetNode("MaskMetrology"); - if( MasKMetrologyNode != NULL ){ - m_MaskMetrologyEnabled = (bool) MasKMetrologyNode->GetValueAsInt(); + if (MasKMetrologyNode != nullptr) { + m_MaskMetrologyEnabled = (bool) MasKMetrologyNode->GetValueAsBoolean(); } MXmlNode* MaskMetrologyFileNameNode = Node->GetNode("MaskMetrologyFileName"); - if (MaskMetrologyFileNameNode != 0) { + if (MaskMetrologyFileNameNode != nullptr) { m_MaskMetrologyFile = MaskMetrologyFileNameNode->GetValue(); } MXmlNode* UCSDOverrideNode = Node->GetNode("UCSDOverride"); - if( UCSDOverrideNode != NULL ){ - m_UCSDOverride = (bool) UCSDOverrideNode->GetValueAsInt(); + if (UCSDOverrideNode != nullptr) { + m_UCSDOverride = (bool) UCSDOverrideNode->GetValueAsBoolean(); } return true; @@ -1029,9 +1043,9 @@ MXmlNode* MModuleDepthCalibration2024::CreateXmlConfiguration() MXmlNode* Node = new MXmlNode(0,m_XmlTag); new MXmlNode(Node, "CoeffsFileName", m_CoeffsFile); new MXmlNode(Node, "SplinesFileName", m_SplinesFile); - new MXmlNode(Node, "MaskMetrology", (unsigned int) m_MaskMetrologyEnabled); + new MXmlNode(Node, "MaskMetrology", (bool)m_MaskMetrologyEnabled); new MXmlNode(Node, "MaskMetrologyFileName", m_MaskMetrologyFile); - new MXmlNode(Node, "UCSDOverride",(unsigned int) m_UCSDOverride); + new MXmlNode(Node, "UCSDOverride", (bool)m_UCSDOverride); return Node; } From e0861b295cae7e05534d2bd7a6060d08b92c7a71 Mon Sep 17 00:00:00 2001 From: Felix Hagemann Date: Sat, 22 Nov 2025 10:54:37 -0800 Subject: [PATCH 11/22] [DEE] Improve strip ID handling for guard ring events --- src/MSubModuleChargeTransport.cxx | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/MSubModuleChargeTransport.cxx b/src/MSubModuleChargeTransport.cxx index 0afd6406..e2b41061 100644 --- a/src/MSubModuleChargeTransport.cxx +++ b/src/MSubModuleChargeTransport.cxx @@ -101,7 +101,8 @@ bool MSubModuleChargeTransport::AnalyzeEvent(MReadOutAssembly* Event) // Calculate LV strip ID by rounding down intentionally to avoid truncation towards zero int ID = static_cast(std::floor((Pos.X() + 7.4/2.0) / (7.4/64.0))); - if (ID >= 0 && ID < 64) { + // Check for strip ID and if the position is within the allowed strip length or on the guard ring + if (ID >= 0 && ID < 64 && std::abs(Pos.Y()) <= 7.4/2.0 && std::hypot(Pos.X(), Pos.Y()) <= 4.712) { SH.m_ROE.SetStripID(ID); SH.m_IsGuardRing = false; } else { @@ -118,7 +119,8 @@ bool MSubModuleChargeTransport::AnalyzeEvent(MReadOutAssembly* Event) // Calculate HV strip ID by rounding down intentionally to avoid truncation towards zero int ID = static_cast(std::floor((Pos.Y() + 7.4/2.0) / (7.4/64.0))); - if (ID >= 0 && ID < 64) { + // Check for strip ID and if the position is within the allowed strip length or on the guard ring + if (ID >= 0 && ID < 64 && std::abs(Pos.X()) <= 7.4/2.0 && std::hypot(Pos.X(), Pos.Y()) <= 4.712) { SH.m_ROE.SetStripID(ID); SH.m_IsGuardRing = false; } else { From 7c3d1d616c5a8d7ef1741cda1155bc86bfdd5efc Mon Sep 17 00:00:00 2001 From: Felix Hagemann Date: Mon, 24 Nov 2025 19:12:37 -0800 Subject: [PATCH 12/22] Add TODOs --- src/MSubModuleChargeTransport.cxx | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/MSubModuleChargeTransport.cxx b/src/MSubModuleChargeTransport.cxx index e2b41061..edf770db 100644 --- a/src/MSubModuleChargeTransport.cxx +++ b/src/MSubModuleChargeTransport.cxx @@ -99,9 +99,11 @@ bool MSubModuleChargeTransport::AnalyzeEvent(MReadOutAssembly* Event) MVector Pos = SH.m_SimulatedPositionInDetector; // Calculate LV strip ID by rounding down intentionally to avoid truncation towards zero + // TODO: Confirm the correct strip pitch based on SMEX detector models int ID = static_cast(std::floor((Pos.X() + 7.4/2.0) / (7.4/64.0))); // Check for strip ID and if the position is within the allowed strip length or on the guard ring + // TODO: Confirm the correct boundary of the guard ring based on SMEX detector models if (ID >= 0 && ID < 64 && std::abs(Pos.Y()) <= 7.4/2.0 && std::hypot(Pos.X(), Pos.Y()) <= 4.712) { SH.m_ROE.SetStripID(ID); SH.m_IsGuardRing = false; @@ -117,9 +119,11 @@ bool MSubModuleChargeTransport::AnalyzeEvent(MReadOutAssembly* Event) MVector Pos = SH.m_SimulatedPositionInDetector; // Calculate HV strip ID by rounding down intentionally to avoid truncation towards zero + // TODO: Confirm the correct strip pitch based on SMEX detector models int ID = static_cast(std::floor((Pos.Y() + 7.4/2.0) / (7.4/64.0))); // Check for strip ID and if the position is within the allowed strip length or on the guard ring + // TODO: Confirm the correct boundary of the guard ring based on SMEX detector models if (ID >= 0 && ID < 64 && std::abs(Pos.X()) <= 7.4/2.0 && std::hypot(Pos.X(), Pos.Y()) <= 4.712) { SH.m_ROE.SetStripID(ID); SH.m_IsGuardRing = false; From 3d3e0e3bbd6c2340c5dd34fe757c547eb4f1fa00 Mon Sep 17 00:00:00 2001 From: Andreas Zoglauer Date: Fri, 24 Oct 2025 16:44:29 -0500 Subject: [PATCH 13/22] ADD: How to handle ToDo's --- CodingConventions.md | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/CodingConventions.md b/CodingConventions.md index 00a13ef1..3798a748 100644 --- a/CodingConventions.md +++ b/CodingConventions.md @@ -22,6 +22,19 @@ To format a file do the following: 4. Do your pull request / check it in +## How to handle future ToDo's in the code + +1. Mark them clearly in a way that IDE's can parse them. Choose one of: + ``` + // TODO: + or + // FIXME: + or + // HACK: + ``` +2. Open an issue for each of them to keep track of them + + ## New classes For new classes, copy and modify an existing one, or use the MModuleTemplate class as tempplate. From 3eeacffe99ed4c275cabba71696402cc8f1d4a86 Mon Sep 17 00:00:00 2001 From: Andreas Zoglauer Date: Fri, 24 Oct 2025 17:08:31 -0500 Subject: [PATCH 14/22] CHG: Updated README.md --- README.md | 72 +++++++++++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 70 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index cf0ffb6f..cf335fbd 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,70 @@ -# nuclearizer -COSI's calibration tool +# COSI's calibration tool nuclearizer + +## What is nuclearizer? + +Nuclearizer is the detectior calibration tool of the gamma-ray telescope COSI, the Compton Spectrometer and Imager. +It handles (at least partly) the detector calibration parameter generation, the appliation of the calibration to measured data, and it as a detector effects engine to make simulation look like real data. +Nuclearizer is written in C++ and based on MEGAlib's detector calibration framework Fretalon. + +## Short history + +The ancestor of Nuclearizer is the calibration tool of the gamma-ray telescope MEGA, which was called MEGAlyze. +It inherited the data flow, the class organization, and the fact that it was based on CERN's ROOT library. +The development of nuclearizer started the year before COSI's second Fort Sumner balloon flight in 2008. +The code base was significanly rewritten and reorganized for version 2 of nuclearizer which was developped for the 2014 Antarctic balloon flight opd COSI. This is the first version based on MEGAlib's Fretalon calibration framework. +The current version is written for the COSI SMEX mission, which resulted in more code cleanups and further code reorganizations + +## Installing nuclearizer + +1. Nuclearizer is based on MEGAlib. Thus the first step is to install MEGAlib: https://github.com/zoglauer/MEGAlib +2. Then set the nuclearizer environment variable. For example in bash/zsh do: + ``` + export NUCLEARIZER=/path/to/nuclearziers/main/directory + ``` +3. Then switch to the nuclearzier main directory and do + ``` + make -j10 + ``` + Where 10 is the number of cores you want to use for compilation +4. Launch nuclearzier via + ``` + nuclearizer + ``` + +## Contributing + +Everyboidy is welcome to contribute. +We follow the standard fork-pull-request workflow. + +Please take a look at the file CodingConventions.md for the coding conventions and other programming tips. + +### Debugging + +In case you need to debug nuclearizer do the following: +1. First you have to compile MEGALib and nuclearizer in debug mode: + ``` + cd $MEGALIB + bash configure --deb=on --opt=off + make -j10 + cd $NUCLEARIZER + make clean + make -j 10 + ``` +2. Then start nuclearizer in debug mode + ``` + debug nuclearizer [you can add command line options here - it automatically runs] + ``` + This works both on macOS and Linux +3. Then when nuclearizer crashes, use the standard gdb / lldb syntax to figure out what went wrong: + ``` + bt (for backtrace) + up (to move to the crash in nuclearizer/megalib code) + p (fro priny the content of a variable) + ``` + +### Issues + +Feel free to report all issues in GitHub's issue tracker. + + + From 5038c55518183d3baa88725e641b70bd6370b072 Mon Sep 17 00:00:00 2001 From: Andreas Zoglauer Date: Fri, 24 Oct 2025 17:20:31 -0500 Subject: [PATCH 15/22] CHG: Add historic contributors --- README.md | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/README.md b/README.md index cf335fbd..16210f64 100644 --- a/README.md +++ b/README.md @@ -66,5 +66,20 @@ In case you need to debug nuclearizer do the following: Feel free to report all issues in GitHub's issue tracker. +## Acknowledgements + +We want to thank all the contributors not listed in GitHub contributor list: + ++ Mark Bandstra ++ Eric Bellm ++ Alan Chiu ++ Daniel Perez ++ Carolyn Kierans ++ Clio Sleator ++ Alex Lowell ++ Hadar Lazar ++ Jacqueline Beechert ++ plus everybody listed on GitHub + From 2fb1580944a9dba08eac040120fb4b2ea3b67830 Mon Sep 17 00:00:00 2001 From: ckierans Date: Tue, 9 Dec 2025 23:46:50 -0500 Subject: [PATCH 16/22] Resolving merge conflict with Sean's recent PR --- include/MModuleDepthCalibration2024.h | 10 +- src/MModuleDepthCalibration2024.cxx | 211 ++++++++++++++++---------- 2 files changed, 138 insertions(+), 83 deletions(-) diff --git a/include/MModuleDepthCalibration2024.h b/include/MModuleDepthCalibration2024.h index d0cb51c5..17f0f671 100644 --- a/include/MModuleDepthCalibration2024.h +++ b/include/MModuleDepthCalibration2024.h @@ -107,12 +107,13 @@ class MModuleDepthCalibration2024 : public MModule //! Retrieve the appropriate CTD values given the DetID and Grade vector GetCTD(int DetID, int Grade); - + + //! Retrieve the appropriate depth-to-CTD spline given the DetID and Grade + TSpline3* GetSpline(int DetID, int Grade); //! Normal distribution vector norm_pdf(vector x, double mu, double sigma); - - //! Adds a Depth-to-CTD relation - bool AddDepthCTD(vector depthvec, vector> ctdarr, int DetID, unordered_map>& DepthGrid, unordered_map>>& CTDMap); + //! Adds a Depth-to-CTD relation + bool AddDepthCTD(vector Depth, vector> CTDArr, int DetID, unordered_map>& DepthGrid, unordered_map>>& CTDMap, unordered_map>& SplineMap, unsigned int NPoints); //! Determine the Grade (geometry of charge sharing) of the Hit int GetHitGrade(MHit* H); @@ -173,6 +174,7 @@ class MModuleDepthCalibration2024 : public MModule // The CTD Map maps each detector (int) to a 2D array of CTD values. unordered_map>> m_CTDMap; unordered_map> m_DepthGrid; + unordered_map> m_SplineMap; bool m_SplinesFileIsLoaded; bool m_CoeffsFileIsLoaded; bool m_MaskMetrologyFileIsLoaded; diff --git a/src/MModuleDepthCalibration2024.cxx b/src/MModuleDepthCalibration2024.cxx index a6c3bd70..2e4a6d02 100644 --- a/src/MModuleDepthCalibration2024.cxx +++ b/src/MModuleDepthCalibration2024.cxx @@ -260,8 +260,8 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) // Calculate the position. If error is thrown, record and no depth. // Take a Hit and separate its activated X- and Y-strips into separate vectors. - std::vector LVStrips; - std::vector HVStrips; + vector LVStrips; + vector HVStrips; for( unsigned int j = 0; j < H->GetNStripHits(); ++j){ MStripHit* SH = H->GetStripHit(j); if( SH->IsLowVoltageStrip() ) LVStrips.push_back(SH); else HVStrips.push_back(SH); @@ -315,6 +315,9 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) vector* Coeffs = GetPixelCoeffs(PixelCode); + vector CTDVec = GetCTD(DetID, Grade); + vector DepthVec = GetDepth(DetID); + // TODO: For Card Cage, may need to add noise double LVTiming = LVSH->GetTiming(); double HVTiming = HVSH->GetTiming(); @@ -325,34 +328,29 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) H->SetNoDepth(); Event->SetDepthCalibrationIncomplete(); ++m_Error1; - } - // If there isn't timing information, set no depth. - // Alex's old comments suggest assigning the event to the middle of the detector and the position resolution to be large. - else if( (LVTiming < 1.0E-6) || (HVTiming < 1.0E-6) ){ + } else if (CTDVec.size() == 0) { + cout << "Empty CTD vector" << endl; + H->SetNoDepth(); + Event->SetDepthCalibrationIncomplete(); + } else if (DepthVec.size() == 0) { + cout << "Empty Depth vector" << endl; + H->SetNoDepth(); + Event->SetDepthCalibrationIncomplete(); + } else if ((LVTiming < 1.0E-6) || (HVTiming < 1.0E-6)) { ++m_Error3; H->SetNoDepth(); Event->SetDepthCalibrationIncomplete(); - } - - // If there are coefficients and timing information is loaded, try calculating the CTD and depth - else { - - vector ctdvec = GetCTD(DetID, Grade); - vector depthvec = GetDepth(DetID); - - if ( ctdvec.size() == 0){ - cout << "Empty CTD vector" << endl; - H->SetNoDepth(); - Event->SetDepthCalibrationIncomplete(); - } + } else { + + // If there are coefficients and timing information is loaded, try calculating the CTD and depth double CTD = (HVTiming - LVTiming); // Confirmed that this matches SP's python code. CTD_s = (CTD - 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 Xmin = * std::min_element(CTDVec.begin(), CTDVec.end()); + double Xmax = * std::max_element(CTDVec.begin(), CTDVec.end()); double noise = GetTimingNoiseFWHM(PixelCode, H->GetEnergy()); @@ -365,9 +363,9 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) // If the CTD is in range, calculate the depth else { - // Calculate the probability given timing noise of CTD_s corresponding to the values of depth in depthvec + // 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); + vector prob_dist = norm_pdf(CTDVec, CTD_s, noise/2.355); // Weight the depth by probability double prob_sum = 0.0; @@ -376,20 +374,20 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) } //double prob_sum = std::accumulate(prob_dist.begin(), prob_dist.end(), 0); double weighted_depth = 0.0; - for( unsigned int k = 0; k < depthvec.size(); ++k ){ - weighted_depth += prob_dist[k]*depthvec[k]; + 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; kIsStripPairingIncomplete()==false) { @@ -562,24 +560,17 @@ vector MModuleDepthCalibration2024::norm_pdf(vector x, double mu bool MModuleDepthCalibration2024::LoadSplinesFile(MString FName) { - //when invert flag is set to true, the splines returned are CTD->Depth - // Previously saved cathode and anode timing in addition to CTD. This may be redundant, commenting out for now. // Input spline files should have the following format: // ### DetID, HV, Temperature, Photopeak Energy (TODO: More? Fewer?) // depth, ctd0, ctd1, ctd2.... (Basically, allow for CTDs for different subpixel regions) // '' '' '' - MFile F; + MFile F; if( F.Open(FName) == false ){ cout << "ERROR in MModuleDepthCalibration2024::LoadSplinesFile: failed to open splines file." << endl; return false; } - // vector depthvec, ctdvec, anovec, catvec; - vector depthvec; - vector> ctdarr; - for( unsigned int i=0; i < 5; ++i ){ - vector temp_vec; - ctdarr.push_back(temp_vec); - } + vector DepthVec; + vector> CTDArr; bool Result = true; MString line; int DetID = 0; @@ -588,33 +579,28 @@ bool MModuleDepthCalibration2024::LoadSplinesFile(MString FName) if (line.BeginsWith("#")) { // If we've reached a new ctd spline then record the previous one in the m_SplineMaps and start a new one. vector tokens = line.Tokenize(" "); - if (depthvec.size() > 0) { - Result &= AddDepthCTD(depthvec, ctdarr, DetID, m_DepthGrid, m_CTDMap); - } - depthvec.clear(); ctdarr.clear(); - for (unsigned int i=0; i < 5; ++i) { - vector temp_vec; - ctdarr.push_back(temp_vec); + if (DepthVec.size() > 0) { + Result &= AddDepthCTD(DepthVec, CTDArr, DetID, m_DepthGrid, m_CTDMap, m_SplineMap, 1000); } + DepthVec.clear(); CTDArr.clear(); DetID = tokens[1].ToInt(); } else { vector tokens = line.Tokenize(","); - depthvec.push_back(tokens[0].ToDouble()); - + DepthVec.push_back(tokens[0].ToDouble()); // Multiple CTDs allowed. for (unsigned int i = 0; i < (tokens.size() - 1); ++i) { - ctdarr[i].push_back(tokens[1+i].ToDouble()); - } - // Fill in the higher grades with the GRADE=0 CTD if there are none listed in the file. - for (unsigned int i=tokens.size()-1; i<5; ++i) { - ctdarr[i].push_back(tokens[1].ToDouble()); + while (i>=CTDArr.size()) { + vector TempVec; + CTDArr.push_back(TempVec); + } + CTDArr[i].push_back(tokens[1+i].ToDouble()); } } } } //make last spline - if (depthvec.size() > 0) { - Result &= AddDepthCTD(depthvec, ctdarr, DetID, m_DepthGrid, m_CTDMap); + if (DepthVec.size() > 0) { + Result &= AddDepthCTD(DepthVec, CTDArr, DetID, m_DepthGrid, m_CTDMap, m_SplineMap, 1000); } return Result; @@ -817,44 +803,91 @@ int MModuleDepthCalibration2024::GetHitGrade(MHit* H){ ///////////////////////////////////////////////////////////////////////////////// -bool MModuleDepthCalibration2024::AddDepthCTD(vector depthvec, vector> ctdarr, int DetID, unordered_map>& DepthGrid, unordered_map>>& CTDMap){ - +bool MModuleDepthCalibration2024::AddDepthCTD(vector Depth, vector> CTDArr, int DetID, unordered_map>& DepthGrid, unordered_map>>& CTDMap, unordered_map>& SplineMap, unsigned int NPoints) +{ // Saves a CTD array, basically allowing for multiple CTDs as a function of depth - // depthvec: list of simulated depth values - // ctdarr: vector of vectors of simulated CTD values. Each vector of CTDs must be the same length as depthvec + // Depth: list of simulated depth values + // CTDArr: vector of vectors of simulated CTD values. Each vector of CTDs must be the same length as Depth // DetID: An integer which specifies which detector. // CTDMap: unordered map into which the array of CTDs should be placed // TODO: Possible energy dependence of CTD? - // TODO: Depth values need to be evenly spaced. Check this when reading the files in. // Check to make sure things look right. // First check that the CTDs all have the right length. - for (unsigned int i = 0; i < ctdarr.size(); ++i) { - if ((ctdarr[i].size() != depthvec.size()) && (ctdarr[i].size() > 0)) { + + for (unsigned int i = 0; i < CTDArr.size(); ++i) { + if ((CTDArr[i].size() != Depth.size()) && (CTDArr[i].size() > 0)) { cout<<"ERROR in MModuleDepthCalibration2024::AddDepthCTD: The number of values in the CTD list is not equal to the number of depth values."< 0.01) { + double MaxDepth = * std::max_element(Depth.begin(), Depth.end()); + double MinDepth = * std::min_element(Depth.begin(), Depth.end()); + if (fabs((MaxDepth-MinDepth) - m_Thicknesses[DetID]) > 0.01) { cout<<"ERROR in MModuleDepthCalibration2024::AddDepthCTD: The thickness of detector "<0) { + cout<<"MModuleDepthCalibration2024::AddSplines: Splines already added for DetID "< TempVec; + SplineMap[DetID] = TempVec; + } + + double dx_low = Depth[1] - Depth[0]; + double newx_low = Depth[0] - (dx_low/2.0); + Depth.insert(Depth.begin(), newx_low); + + size_t N = Depth.size(); + double dx_high = Depth[N-1] - Depth[N-2]; + double newx_high = Depth[N-1] + (dx_high/2.0); + Depth.push_back(newx_high); + + if (NPoints < N+1) { + NPoints = N+1; + } + + vector NewDepth; + for (unsigned int i=0; i> NewCTDArr; + for (unsigned int i=0; i CTD = CTDArr[i]; + vector NewCTD; - //Now make sure the values for the depth start with 0.0. - if (mindepth != 0.0) { - cout<<"MModuleDepthCalibration2024::AddDepthCTD: The minimum depth is not zero. Editing the depth vector."<Eval(NewDepth[d])); + } + NewCTDArr.push_back(NewCTD); } - CTDMap[DetID] = ctdarr; - DepthGrid[DetID] = depthvec; + CTDMap[DetID] = NewCTDArr; return true; } @@ -877,7 +910,7 @@ vector MModuleDepthCalibration2024::GetCTD(int DetID, int Grade) return (m_CTDMap[DetID][Grade]); } else { - cout << "MModuleDepthCalibration2024::GetCTD: No CTD map is loaded for Grade " << Grade << ". Returning Grade 0 CTD." << endl; + // cout << "MModuleDepthCalibration2024::GetCTD: No CTD map is loaded for Grade " << Grade << ". Returning Grade 0 CTD." << endl; return (m_CTDMap[DetID][0]); } } else { @@ -912,6 +945,31 @@ vector MModuleDepthCalibration2024::GetDepth(int DetID) ///////////////////////////////////////////////////////////////////////////////// +TSpline3* MModuleDepthCalibration2024::GetSpline(int DetID, int Grade) +{ + // Retrieves the appropriate depth->CTD spline given the Detector ID and Event Grade passed + + if( !m_SplinesFileIsLoaded ){ + cout << "MModuleDepthCalibration2024::GetSpline: cannot return Depth to CTD spline because the file was not loaded." << endl; + return nullptr; + } + + // If there is a spline for the given detector, return it. + // If the Grade is larger than the number of splines stored, then just return Grade 0 spline. + if( m_SplineMap.count(DetID) > 0 ){ + if ( ((int)m_SplineMap[DetID].size()) > Grade) { + return (m_SplineMap[DetID][Grade]); + } + else { + // cout << "MModuleDepthCalibration2024::GetSpline: No spline is loaded for Grade " << Grade << ". Returning Grade 0 spline." << endl; + return (m_SplineMap[DetID][0]); + } + } else { + cout << "MModuleDepthCalibration2024::GetSpline: No spline is loaded for Det " << DetID << "." << endl; + return nullptr; + } +} +>>>>>>> 1899e35 (CHG: Update spline handling to ensure the full detector is sampled and the depth grid is evenly spaced) void MModuleDepthCalibration2024::ShowOptionsGUI() { @@ -991,11 +1049,6 @@ void MModuleDepthCalibration2024::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; - /* - TFile* rootF = new TFile("EHist.root","recreate"); - rootF->WriteTObject( EHist ); - rootF->Close(); - */ } From 230f1907c94cd5c758065a210a5ef8eb10c09be1 Mon Sep 17 00:00:00 2001 From: Sean Pike Date: Tue, 23 Sep 2025 15:12:06 -0700 Subject: [PATCH 17/22] Cleanup: added comments and cleaned up extraneous print statements --- src/MModuleDepthCalibration2024.cxx | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/src/MModuleDepthCalibration2024.cxx b/src/MModuleDepthCalibration2024.cxx index 2e4a6d02..f331765e 100644 --- a/src/MModuleDepthCalibration2024.cxx +++ b/src/MModuleDepthCalibration2024.cxx @@ -362,6 +362,9 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) } // 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. @@ -372,7 +375,6 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) for( unsigned int k=0; k < prob_dist.size(); ++k ){ prob_sum += prob_dist[k]; } - //double prob_sum = std::accumulate(prob_dist.begin(), prob_dist.end(), 0); double weighted_depth = 0.0; for( unsigned int k = 0; k < DepthVec.size(); ++k ){ weighted_depth += prob_dist[k]*DepthVec[k]; @@ -548,7 +550,6 @@ vector MModuleDepthCalibration2024::norm_pdf(vector x, double mu vector result; for( unsigned int i=0; i Depth, vector Depth, vector 0.01) { @@ -831,14 +834,16 @@ bool MModuleDepthCalibration2024::AddDepthCTD(vector Depth, vector0) { - cout<<"MModuleDepthCalibration2024::AddSplines: Splines already added for DetID "< TempVec; SplineMap[DetID] = TempVec; } + // Add a new point to the depth vector on the bottom and top of the detector to make sure we're sampling the full depth double dx_low = Depth[1] - Depth[0]; double newx_low = Depth[0] - (dx_low/2.0); Depth.insert(Depth.begin(), newx_low); @@ -848,16 +853,20 @@ bool MModuleDepthCalibration2024::AddDepthCTD(vector Depth, vector NewDepth; for (unsigned int i=0; i> NewCTDArr; for (unsigned int i=0; i CTD = CTDArr[i]; @@ -910,7 +919,6 @@ vector MModuleDepthCalibration2024::GetCTD(int DetID, int Grade) return (m_CTDMap[DetID][Grade]); } else { - // cout << "MModuleDepthCalibration2024::GetCTD: No CTD map is loaded for Grade " << Grade << ". Returning Grade 0 CTD." << endl; return (m_CTDMap[DetID][0]); } } else { @@ -961,7 +969,6 @@ TSpline3* MModuleDepthCalibration2024::GetSpline(int DetID, int Grade) return (m_SplineMap[DetID][Grade]); } else { - // cout << "MModuleDepthCalibration2024::GetSpline: No spline is loaded for Grade " << Grade << ". Returning Grade 0 spline." << endl; return (m_SplineMap[DetID][0]); } } else { From 649f74c0b92aab8467351131c7a980be01c6a2e1 Mon Sep 17 00:00:00 2001 From: Sean Pike Date: Tue, 23 Sep 2025 15:35:09 -0700 Subject: [PATCH 18/22] BUG: Clear variables on Finalize --- src/MModuleDepthCalibration2024.cxx | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/src/MModuleDepthCalibration2024.cxx b/src/MModuleDepthCalibration2024.cxx index f331765e..17dd195a 100644 --- a/src/MModuleDepthCalibration2024.cxx +++ b/src/MModuleDepthCalibration2024.cxx @@ -1057,6 +1057,19 @@ void MModuleDepthCalibration2024::Finalize() cout << "Number of hits with null strip hits: " << m_ErrorNullSH << endl; cout << "Number of hits 0 energy on a strip hit: " << m_ErrorNoE << endl; + // Clean up maps and vectors + m_Coeffs.clear(); + m_Thicknesses.clear(); + m_NXStrips.clear(); + m_NYStrips.clear(); + m_XPitches.clear(); + m_YPitches.clear(); + m_Detectors.clear(); + m_CTDMap.clear(); + m_DepthGrid.clear(); + m_SplineMap.clear(); + m_DetectorIDs.clear(); + } From f5c13f0411077b93cd15d96965e9f67d9cf8fc6f Mon Sep 17 00:00:00 2001 From: Felix Hagemann Date: Tue, 9 Dec 2025 09:10:57 -0800 Subject: [PATCH 19/22] [DEE] Read detector dimensions from the geometry (#76) * Set the geometry in SMEX DEE SubModules * Read detector dimensions from the geometry * Read detector radius from geometry --- include/MDEEStripHit.h | 2 +- include/MSubModuleChargeTransport.h | 18 ++++++ src/MModuleDEESMEX.cxx | 3 + src/MSubModuleChargeTransport.cxx | 96 ++++++++++++++++++++++++++--- 4 files changed, 110 insertions(+), 9 deletions(-) diff --git a/include/MDEEStripHit.h b/include/MDEEStripHit.h index b57d3da4..44aff194 100644 --- a/include/MDEEStripHit.h +++ b/include/MDEEStripHit.h @@ -60,7 +60,7 @@ struct MDEEStripHit double m_SimulatedRelativeDepth; //! The energy from simulations double m_SimulatedEnergy; // original: m_EnergyOrig - //! The list of origin IDs form the simulation + //! The list of origin IDs from the simulation list m_SimulatedOrigins; // original: m_Origins //! True if this is a guard ring bool m_SimulatedIsGuardRing; // original: m_IsGuardRing diff --git a/include/MSubModuleChargeTransport.h b/include/MSubModuleChargeTransport.h index 96e84790..9a7cbf67 100644 --- a/include/MSubModuleChargeTransport.h +++ b/include/MSubModuleChargeTransport.h @@ -23,6 +23,7 @@ // MEGAlib libs: #include "MGlobal.h" #include "MSubModule.h" +#include "MDStrip3D.h" // Forward declarations: @@ -50,6 +51,9 @@ class MSubModuleChargeTransport : public MSubModule //! Default destructor virtual ~MSubModuleChargeTransport(); + //! Set geometry + void SetGeometry(MDGeometryQuest* Geometry) { m_Geometry = Geometry; } + //! Initialize the module virtual bool Initialize(); @@ -77,6 +81,20 @@ class MSubModuleChargeTransport : public MSubModule // protected members: protected: + //! The geometry + MDGeometryQuest* m_Geometry; + + //! The detector dimensions + unordered_map m_Thicknesses; + unordered_map m_NXStrips; + unordered_map m_NYStrips; + unordered_map m_XPitches; + unordered_map m_YPitches; + unordered_map m_XWidths; + unordered_map m_YWidths; + unordered_map m_Radii; + unordered_map m_Detectors; + vector m_DetectorIDs; // private members: diff --git a/src/MModuleDEESMEX.cxx b/src/MModuleDEESMEX.cxx index 4858d822..d9d00223 100644 --- a/src/MModuleDEESMEX.cxx +++ b/src/MModuleDEESMEX.cxx @@ -90,6 +90,9 @@ MModuleDEESMEX::~MModuleDEESMEX() bool MModuleDEESMEX::Initialize() { + // Set the geometry to the SubModules using it + m_ChargeTransport.SetGeometry(m_Geometry); + // Initialize the module // Each Initialize() should handle its own error messaging diff --git a/src/MSubModuleChargeTransport.cxx b/src/MSubModuleChargeTransport.cxx index edf770db..4e67ba38 100644 --- a/src/MSubModuleChargeTransport.cxx +++ b/src/MSubModuleChargeTransport.cxx @@ -32,6 +32,8 @@ // MEGAlib libs: #include "MSubModule.h" +#include "MDShapeIntersection.h" +#include "MDShapeTUBS.h" //////////////////////////////////////////////////////////////////////////////// @@ -67,7 +69,70 @@ MSubModuleChargeTransport::~MSubModuleChargeTransport() bool MSubModuleChargeTransport::Initialize() { - // Initialize the module + + // The detectors need to be in the same order as DetIDs. + // ie DetID=0 should be the 0th detector in m_Detectors, DetID=1 should the 1st, etc. + vector DetList = m_Geometry->GetDetectorList(); + + // Look through the Geometry and get the names and thicknesses of all the detectors. + for(unsigned int i = 0; i < DetList.size(); ++i){ + + unsigned int DetID = i; + + MDDetector* det = DetList[i]; + vector DetectorNames; + if (det->GetTypeName() == "Strip3D") { + if (det->GetNSensitiveVolumes() == 1) { + MDVolume* vol = det->GetSensitiveVolume(0); + string det_name = vol->GetName().GetString(); + if (find(DetectorNames.begin(), DetectorNames.end(), det_name) == DetectorNames.end()) { + DetectorNames.push_back(det_name); + m_Thicknesses[DetID] = 2*(det->GetStructuralSize().GetZ()); + MDStrip3D* strip = dynamic_cast(det); + m_XPitches[DetID] = strip->GetPitchX(); + m_YPitches[DetID] = strip->GetPitchY(); + m_NXStrips[DetID] = strip->GetNStripsX(); + m_NYStrips[DetID] = strip->GetNStripsY(); + m_XWidths[DetID] = strip->GetWidthX(); + m_YWidths[DetID] = strip->GetWidthY(); + + // Read the detector radius from the geometry (assuming it is the second shape of an intersection) + if (vol->GetShape()->GetType() == "Intersection" && dynamic_cast(vol->GetShape())->GetShapeB()->GetType() == "TUBS") { + m_Radii[DetID] = dynamic_cast(dynamic_cast(vol->GetShape())->GetShapeB())->GetRmax(); + } else { + // If that does not exist, set the detector radius to a value high enough to not have an impact + m_Radii[DetID] = m_XWidths[DetID] + m_YWidths[DetID]; + if (g_Verbosity >= c_Info) { + cout << m_Name << ": No bounding tube volume found for this detector" << endl; + } + } + + if (g_Verbosity >= c_Info) { + cout << "Found detector " << det_name << " corresponding to DetID=" << DetID << "." << endl; + cout << "Detector width (X): " << m_XWidths[DetID] << endl; + cout << "Detector width (Y): " << m_YWidths[DetID] << endl; + cout << "Detector radius (R): " << m_Radii[DetID] << endl; + cout << "Detector thickness: " << m_Thicknesses[DetID] << endl; + cout << "Number of X strips: " << m_NXStrips[DetID] << endl; + cout << "Number of Y strips: " << m_NYStrips[DetID] << endl; + cout << "X strip pitch: " << m_XPitches[DetID] << endl; + cout << "Y strip pitch: " << m_YPitches[DetID] << endl; + } + m_DetectorIDs.push_back(DetID); + m_Detectors[DetID] = det; + } else { + cout << "ERROR in MSubModuleChargeTransport::Initialize: Found a duplicate detector: " << det_name << endl; + } + } else { + cout << "ERROR in MSubModuleChargeTransport::Initialize: Found a Strip3D detector with " << det->GetNSensitiveVolumes() << " Sensitive Volumes." << endl; + } + } + } + + if (m_DetectorIDs.size() == 0) { + cout<<"No Strip3D detectors were found."<(std::floor((Pos.X() + 7.4/2.0) / (7.4/64.0))); + int ID = static_cast(std::floor((Pos.X() + XWidth/2.0) / XPitch)); // Check for strip ID and if the position is within the allowed strip length or on the guard ring // TODO: Confirm the correct boundary of the guard ring based on SMEX detector models - if (ID >= 0 && ID < 64 && std::abs(Pos.Y()) <= 7.4/2.0 && std::hypot(Pos.X(), Pos.Y()) <= 4.712) { + if (ID >= 0 && ID < NXStrips && std::abs(Pos.Y()) <= YWidth/2.0 && std::hypot(Pos.X(), Pos.Y()) <= Radius) { SH.m_ROE.SetStripID(ID); SH.m_IsGuardRing = false; } else { - SH.m_ROE.SetStripID(64); + SH.m_ROE.SetStripID(NXStrips); SH.m_IsGuardRing = true; } SH.m_Energy = SH.m_SimulatedEnergy; @@ -118,17 +190,25 @@ bool MSubModuleChargeTransport::AnalyzeEvent(MReadOutAssembly* Event) for (MDEEStripHit& SH: HVHits) { MVector Pos = SH.m_SimulatedPositionInDetector; + // Determine detector and strip dimensions from the geometry + unsigned int DetID = SH.m_ROE.GetDetectorID(); + double XWidth = m_XWidths[DetID]; + double YWidth = m_YWidths[DetID]; + double YPitch = m_YPitches[DetID]; + double Radius = m_Radii[DetID]; + int NYStrips = m_NYStrips[DetID]; + // Calculate HV strip ID by rounding down intentionally to avoid truncation towards zero // TODO: Confirm the correct strip pitch based on SMEX detector models - int ID = static_cast(std::floor((Pos.Y() + 7.4/2.0) / (7.4/64.0))); + int ID = static_cast(std::floor((Pos.Y() + YWidth/2.0) / YPitch)); // Check for strip ID and if the position is within the allowed strip length or on the guard ring // TODO: Confirm the correct boundary of the guard ring based on SMEX detector models - if (ID >= 0 && ID < 64 && std::abs(Pos.X()) <= 7.4/2.0 && std::hypot(Pos.X(), Pos.Y()) <= 4.712) { + if (ID >= 0 && ID < NYStrips && std::abs(Pos.X()) <= XWidth/2.0 && std::hypot(Pos.X(), Pos.Y()) <= Radius) { SH.m_ROE.SetStripID(ID); SH.m_IsGuardRing = false; } else { - SH.m_ROE.SetStripID(64); + SH.m_ROE.SetStripID(NYStrips); SH.m_IsGuardRing = true; } SH.m_Energy = SH.m_SimulatedEnergy; From 22f32141a51a7a8674683341e5fc83392263755a Mon Sep 17 00:00:00 2001 From: ckierans Date: Tue, 9 Dec 2025 23:49:29 -0500 Subject: [PATCH 20/22] Resolving merge conflict with Sean's recent PR --- src/MModuleDepthCalibration2024.cxx | 39 +++++++++++++---------------- 1 file changed, 18 insertions(+), 21 deletions(-) diff --git a/src/MModuleDepthCalibration2024.cxx b/src/MModuleDepthCalibration2024.cxx index 17dd195a..c8750d3e 100644 --- a/src/MModuleDepthCalibration2024.cxx +++ b/src/MModuleDepthCalibration2024.cxx @@ -281,28 +281,21 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) int PixelCode = 10000*DetID + 100*LVStripID + HVStripID; // TODO: Calculate X and Y positions more rigorously using charge sharing. + // Somewhat confusing notation: HVStrips run parallel to X-axis, so we calculate X position with LVStrips. double Xpos, Ypos, Zpos; - if ( m_MaskMetrologyEnabled ) { // If we are applying the mask metrology correction, first define two new readout elements to help determine the intersection of these two strips - - - //TODO: Aldo's file has the HV strip ID flipped (as evident in the x,y position reported in the metrology files), so this is a placeholder until the files are fixed. - int HVStripID_flipped = 63 - HVStripID; - MStripHit* HVSH_flipped = HVSH; - HVSH_flipped->SetStripID(HVStripID_flipped); MReadOutElementDoubleStrip R_LV = *dynamic_cast(LVSH->GetReadOutElement()); - MReadOutElementDoubleStrip R_HV = *dynamic_cast(HVSH_flipped->GetReadOutElement()); - - //find the intercept of the two dominate strips based on the mask metrology + MReadOutElementDoubleStrip R_HV = *dynamic_cast(HVSH->GetReadOutElement()); + vector inter = GetStripIntersection(R_LV, R_HV); Xpos = inter[0]; Ypos = inter[1]; - + //Xpos = m_YPitches[DetID]*((m_NYStrips[DetID]/2.0) - ((double)LVStripID)); + //Ypos = m_XPitches[DetID]*((m_NXStrips[DetID]/2.0) - ((double)HVStripID)); + Zpos = 0.0; } else { - //If no metrology is used, define the strip positions based on the detector pitch and number of strip hits - Xpos = m_YPitches[DetID]*((m_NYStrips[DetID]/2.0) - ((double)LVStripID)); Ypos = m_XPitches[DetID]*((m_NXStrips[DetID]/2.0) - ((double)HVStripID)); Zpos = 0.0; @@ -615,7 +608,7 @@ bool MModuleDepthCalibration2024::LoadSplinesFile(MString FName) bool MModuleDepthCalibration2024::LoadMaskMetrologyFile(MString FName) { //Read the Mask Metrology File - // Det ID, Side (l,h), Strip ID (0-63), x_mm, y_mm, z_mm, roll_deg, pitch_deg, yaw_deg + // Det ID, Strip ID (0-63), Side (l,h), x_mm, y_mm, z_mm, roll_deg, pitch_deg, yaw_deg MFile F; if( F.Open(FName) == false ){ cout << "ERROR in MModuleDepthCalibration2024::LoadMaskMetrologyFile: failed to open metrology file." << endl; @@ -631,8 +624,8 @@ bool MModuleDepthCalibration2024::LoadMaskMetrologyFile(MString FName) //Define the readout element to track det ID, strip ID, and lv/hv MReadOutElementDoubleStrip R; R.SetDetectorID(Tokens[0].ToInt()); - R.IsLowVoltageStrip((Tokens[1].ToString() == "p") || (Tokens[1].ToString() == "l")); - R.SetStripID(Tokens[2].ToInt()); + R.SetStripID(Tokens[1].ToInt()); + R.IsLowVoltageStrip((Tokens[2].ToString() == "p") || (Tokens[2].ToString() == "l")); double Strip_MetX = Tokens[3].ToDouble()/10; //convert to cm double Strip_MetY = Tokens[4].ToDouble()/10; //convert to cm double Strip_MetZ = Tokens[5].ToDouble()/10; //convert to cm @@ -669,14 +662,14 @@ vector MModuleDepthCalibration2024::GetStripIntersection(MReadOutElement //Find the x position of two lines represented by the dominate strips: // LVstrip is centered at (x,y,z) = (lv_strip_met[0], lv_strip_met[1], lv_strip_met[2]) - // and is approximately parallel to the y axis, but rotated at angle lv_strip_met[5] + // and is approximately parallel to the y axis, but rotated at angle lv_strip_met[3] // around the z axis of the detector // HVstrip is centered at (x,y,z) = (hv_strip_met[0], hv_strip_met[1], hv_strip_met[2]) - // and is approximately parallel to the x axis, but rotated at angle (hv_strip_met[5] - pi/2) + // and is approximately parallel to the x axis, but rotated at angle hv_strip_met[3] // around the z axis of the detector - double x_intercept = (hv_strip_met[0]*tan((hv_strip_met[5]-90)*TMath::DegToRad()) - lv_strip_met[0]/tan(lv_strip_met[5]*TMath::DegToRad()) - lv_strip_met[1] + hv_strip_met[1])/(tan((hv_strip_met[5]-90)*TMath::DegToRad())-1/tan(lv_strip_met[5]*TMath::DegToRad())); + double x_intercept = (hv_strip_met[0]*tan(hv_strip_met[3]*TMath::DegToRad()) - lv_strip_met[0]/tan(lv_strip_met[3]*TMath::DegToRad()) - lv_strip_met[1] + hv_strip_met[1])/(tan(hv_strip_met[3]*TMath::DegToRad())-1/tan(lv_strip_met[3]*TMath::DegToRad())); - double y_intercept = (x_intercept - hv_strip_met[0])*tan((hv_strip_met[5]-90)*TMath::DegToRad()) + hv_strip_met[1]; + double y_intercept = (x_intercept - hv_strip_met[0])*tan(hv_strip_met[3]*TMath::DegToRad()) + hv_strip_met[1]; return {x_intercept, y_intercept}; } @@ -800,7 +793,7 @@ int MModuleDepthCalibration2024::GetHitGrade(MHit* H){ } -///////////////////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////////// bool MModuleDepthCalibration2024::AddDepthCTD(vector Depth, vector> CTDArr, int DetID, unordered_map>& DepthGrid, unordered_map>>& CTDMap, unordered_map>& SplineMap, unsigned int NPoints) @@ -978,6 +971,10 @@ TSpline3* MModuleDepthCalibration2024::GetSpline(int DetID, int Grade) } >>>>>>> 1899e35 (CHG: Update spline handling to ensure the full detector is sampled and the depth grid is evenly spaced) + +///////////////////////////////////////////////////////////////////////////////// + + void MModuleDepthCalibration2024::ShowOptionsGUI() { // Show the options GUI - or do nothing From a5a5de2a69df1b7857a7e237a400c891db22f04c Mon Sep 17 00:00:00 2001 From: ckierans Date: Tue, 2 Sep 2025 09:31:07 -0400 Subject: [PATCH 21/22] CHG: Updated x,y definition with mask rotation, added print statements --- src/MModuleDepthCalibration2024.cxx | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/src/MModuleDepthCalibration2024.cxx b/src/MModuleDepthCalibration2024.cxx index c8750d3e..26392a9a 100644 --- a/src/MModuleDepthCalibration2024.cxx +++ b/src/MModuleDepthCalibration2024.cxx @@ -285,15 +285,21 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) double Xpos, Ypos, Zpos; if ( m_MaskMetrologyEnabled ) { // If we are applying the mask metrology correction, first define two new readout elements to help determine the intersection of these two strips + int HVStripID_flipped = 63 - HVStripID; + HVSH->SetStripID(HVStripID_flipped); MReadOutElementDoubleStrip R_LV = *dynamic_cast(LVSH->GetReadOutElement()); MReadOutElementDoubleStrip R_HV = *dynamic_cast(HVSH->GetReadOutElement()); vector inter = GetStripIntersection(R_LV, R_HV); Xpos = inter[0]; Ypos = inter[1]; - //Xpos = m_YPitches[DetID]*((m_NYStrips[DetID]/2.0) - ((double)LVStripID)); - //Ypos = m_XPitches[DetID]*((m_NXStrips[DetID]/2.0) - ((double)HVStripID)); - + + cout<<"StripLV: "< MModuleDepthCalibration2024::GetStripIntersection(MReadOutElement // HVstrip is centered at (x,y,z) = (hv_strip_met[0], hv_strip_met[1], hv_strip_met[2]) // and is approximately parallel to the x axis, but rotated at angle hv_strip_met[3] // around the z axis of the detector - double x_intercept = (hv_strip_met[0]*tan(hv_strip_met[3]*TMath::DegToRad()) - lv_strip_met[0]/tan(lv_strip_met[3]*TMath::DegToRad()) - lv_strip_met[1] + hv_strip_met[1])/(tan(hv_strip_met[3]*TMath::DegToRad())-1/tan(lv_strip_met[3]*TMath::DegToRad())); + double x_intercept = (hv_strip_met[0]*tan(-hv_strip_met[5]*TMath::DegToRad()) - lv_strip_met[0]/tan(lv_strip_met[5]*TMath::DegToRad()) - lv_strip_met[1] + hv_strip_met[1])/(tan(-hv_strip_met[5]*TMath::DegToRad())-1/tan(lv_strip_met[5]*TMath::DegToRad())); - double y_intercept = (x_intercept - hv_strip_met[0])*tan(hv_strip_met[3]*TMath::DegToRad()) + hv_strip_met[1]; + double y_intercept = (x_intercept - hv_strip_met[0])*tan(hv_strip_met[5]*TMath::DegToRad()) + hv_strip_met[1]; return {x_intercept, y_intercept}; } From 92ed5b72a1bad6d33606ef42f69484daad112259 Mon Sep 17 00:00:00 2001 From: ckierans Date: Wed, 10 Dec 2025 00:09:56 -0500 Subject: [PATCH 22/22] Removing a vim mistype --- include/MModuleDepthCalibration2024.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/MModuleDepthCalibration2024.h b/include/MModuleDepthCalibration2024.h index 6a697761..b2794d47 100644 --- a/include/MModuleDepthCalibration2024.h +++ b/include/MModuleDepthCalibration2024.h @@ -1,4 +1,4 @@ -/*:q +/* * MModuleDepthCalibration2024.h * * Copyright (C) 2008-2008 by Andreas Zoglauer.