diff --git a/include/MGUIOptionsDepthCalibration2024.h b/include/MGUIOptionsDepthCalibration2024.h index 15c72ed..8eaf646 100644 --- a/include/MGUIOptionsDepthCalibration2024.h +++ b/include/MGUIOptionsDepthCalibration2024.h @@ -74,9 +74,18 @@ class MGUIOptionsDepthCalibration2024 : public MGUIOptions //! Select spline file to load, splines will convert CTD->Depth MGUIEFileSelector* m_SplinesFileSelector; + //! Use metrology correction + bool m_UseMaskMetCorr; + //! Check button to use metrology correction + TGCheckButton* m_MaskMetModeCB; + //! Select mask metrology file to load. This gives the translation and rotation for each strip in the detector frame. + MGUIEFileSelector* m_MaskMetrologyFileSelector; + //! Check button if working with the Card Cage at UCSD TGCheckButton* m_UCSDOverride; + //! IDs of check buttons + enum ButtonIDs {c_MetrologyFile}; #ifdef ___CLING___ public: diff --git a/include/MModuleDepthCalibration2024.h b/include/MModuleDepthCalibration2024.h index 17237b4..25d58d6 100644 --- a/include/MModuleDepthCalibration2024.h +++ b/include/MModuleDepthCalibration2024.h @@ -63,19 +63,30 @@ class MModuleDepthCalibration2024 : public MModule virtual void ShowOptionsGUI(); //! Set filename for coefficients file - void SetCoeffsFileName( const MString& FileName) {m_CoeffsFile = FileName;} + void SetCoeffsFileName( const MString& FileName) { m_CoeffsFile = FileName; } //! Get filename for coefficients file - MString GetCoeffsFileName() const {return m_CoeffsFile;} + MString GetCoeffsFileName() const { return m_CoeffsFile; } //! Set filename for CTD->Depth splines - void SetSplinesFileName( const MString& FileName) {m_SplinesFile = FileName;} + void SetSplinesFileName( const MString& FileName) { m_SplinesFile = FileName; } //! Get filename for CTD->Depth splines MString GetSplinesFileName() const {return m_SplinesFile;} + //! Enable/Disable Mask Metrology Correction + void SetMaskMetrologyCorrectionEnable(bool X) { m_MaskMetrologyEnabled = X; } + //! Get enable/disable status of mask metrology correction + bool GetMaskMetrologyCorrectionEnable() 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; } + + //TODO Remove UCSD code here and place within it's own branch //! Set whether the data came from the card cage at UCSD - void SetUCSDOverride( bool Override ) {m_UCSDOverride = Override;} + void SetUCSDOverride( bool Override ) { m_UCSDOverride = Override; } //! Get whether the data came from the card cage at UCSD - bool GetUCSDOverride() const {return m_UCSDOverride;} + bool GetUCSDOverride() const { return m_UCSDOverride; } //! Read the XML configuration @@ -91,24 +102,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 +179,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/include/MReadOutAssembly.h b/include/MReadOutAssembly.h index a8d501b..e4f299a 100644 --- a/include/MReadOutAssembly.h +++ b/include/MReadOutAssembly.h @@ -426,7 +426,7 @@ class MReadOutAssembly : public MReadOutSequence //! TODO change variable name or remove Greedy approach all together double m_EventQuality; - //! True if event has been filtered out + //! True if event has been filtered out bool m_FilteredOut; //! The analysis progress diff --git a/src/MGUIOptionsDepthCalibration2024.cxx b/src/MGUIOptionsDepthCalibration2024.cxx index cd6c07e..ff2490f 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_MetrologyFile); + m_MaskMetModeCB->SetState((dynamic_cast(m_Module)->GetMaskMetrologyCorrectionEnable() == true) ? kButtonDown : kButtonUp); + m_MaskMetModeCB->Associate(this); + m_OptionsFrame->AddFrame(m_MaskMetModeCB, LabelLayout); + + m_UseMaskMetCorr = dynamic_cast(m_Module)->GetMaskMetrologyCorrectionEnable(); + + 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", "*.metrology.csv"); + m_OptionsFrame->AddFrame(m_MaskMetrologyFileSelector, LabelLayout); + + + if (m_UseMaskMetCorr == true) { + 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_MetrologyFile: + if (m_MaskMetModeCB->GetState() == kButtonDown) { + m_UseMaskMetCorr = true; + m_MaskMetrologyFileSelector->SetEnabled(true); + } else if (m_MaskMetModeCB->GetState() == kButtonUp) { + m_UseMaskMetCorr = false; + 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)->GetMaskMetrologyCorrectionEnable() != m_UseMaskMetCorr) dynamic_cast(m_Module)->SetMaskMetrologyCorrectionEnable(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 d129969..8ec2114 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" @@ -60,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); @@ -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(); @@ -178,6 +178,15 @@ bool MModuleDepthCalibration2024::Initialize() return false; } + if (m_MaskMetrologyEnabled == true) { + if (g_Verbosity >= c_Info) cout << m_XmlTag << ": !!! Mask Metrology Enabled !!!" << endl; + m_MaskMetrologyFileIsLoaded = LoadMaskMetrologyFile(m_MaskMetrologyFile); + if (m_MaskMetrologyFile == false) { + if (g_Verbosity >= c_Error) cout << m_XmlTag << "Unable to open Metrology file" << endl; + return false; + } + } + MSupervisor* S = MSupervisor::GetSupervisor(); m_EnergyCalibration = (MModuleEnergyCalibrationUniversal*) S->GetAvailableModuleByXmlTag("EnergyCalibrationUniversal"); if (m_EnergyCalibration == nullptr) { @@ -200,26 +209,30 @@ 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); } +///////////////////////////////////////////////////////////////////////////////// + + bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) { - if (Event->GetGuardRingVeto()==true) { + if (Event->GetGuardRingVeto() == true) { + //TODO: Handle events with GR vetos Event->SetDepthCalibrationError("GR Veto"); return false; } 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. @@ -229,7 +242,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->SetDepthCalibrationError("Error in depth calibration"); if (Grade == -1) { @@ -249,15 +262,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; @@ -273,12 +285,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. - // 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)); + //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; + // 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, and update Xpos and Ypos + vector inter = GetStripIntersection(R_LV, R_HV); + Xpos = inter[0]; + Ypos = inter[1]; + + } + + + // 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); double Zsigma = m_Thicknesses[DetID]/sqrt(12.0); @@ -299,11 +328,11 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) Event->SetDepthCalibrationError("No calibration coefficients"); ++m_Error1; } else if (CTDVec.size() == 0) { - cout << "Empty CTD vector" << endl; + if (g_Verbosity >= c_Error) cout << m_XmlTag << "Empty CTD vector" << endl; H->SetNoDepth(); Event->SetDepthCalibrationError("No calibration coefficients"); } else if (DepthVec.size() == 0) { - cout << "Empty Depth vector" << endl; + if (g_Verbosity >= c_Error) cout << m_XmlTag << "Empty Depth vector" << endl; H->SetNoDepth(); Event->SetDepthCalibrationError("No calibration coefficients"); } else if ((LVTiming < 1.0E-6) || (HVTiming < 1.0E-6)) { @@ -313,7 +342,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. @@ -325,7 +353,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->SetDepthCalibrationError("Out of Range"); ++m_Error2; @@ -342,24 +370,28 @@ 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; kHasStripPairingError()==false) { @@ -371,13 +403,16 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) } } - LocalPosition.SetXYZ(Xpos, Ypos, Zpos); - LocalOrigin.SetXYZ(0.0,0.0,0.0); - GlobalPosition = m_Detectors[DetID]->GetSensitiveVolume(0)->GetPositionInWorldVolume(LocalPosition); + //cout << "Strip ID :" << LVStripID << " " << HVStripID << endl; + //cout << "Hit position: "<< Xpos << " " << Ypos << " " << Zpos << endl; + + 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); @@ -394,6 +429,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 @@ -417,69 +456,81 @@ 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. // 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; } -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(); + if (g_Verbosity >= c_Info) cout << m_XmlTag << "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; } + +///////////////////////////////////////////////////////////////////////////////// + + 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 { if (g_Verbosity >= c_Warning) { @@ -494,45 +545,59 @@ std::vector* MModuleDepthCalibration2024::GetPixelCoeffs(int PixelCode) } + +///////////////////////////////////////////////////////////////////////////////// + + vector MModuleDepthCalibration2024::norm_pdf(vector x, double mu, double sigma) { vector result; - for( unsigned int i=0; i 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()) { @@ -544,6 +609,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); @@ -553,18 +621,111 @@ bool MModuleDepthCalibration2024::LoadSplinesFile(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 MetrologyFile; + if (MetrologyFile.Open(FileName) == false) { + cout << "ERROR in MModuleDepthCalibration2024::LoadMaskMetrologyFile: failed to open metrology file." << endl; + return false; + } + + 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; + } + } + } + + MetrologyFile.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 + + // 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] + // 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) + // around the z axis of the detector + + vector LVStripMet = m_MaskMetrology[R_LVStrip]; + vector HVStripMet = m_MaskMetrology[R_HVStrip]; + int DetID = R_LVStrip.GetDetectorID(); + + // Check for division by zero and return standard intersection without mask rotation in this case. + // These values will be the same for every detector since the rotation is the same for each strip within a detector. + double denominator1 = tan(LVStripMet[5]*TMath::DegToRad()); + double denominator2 = tan((HVStripMet[5]-90)*TMath::DegToRad())-1/tan(LVStripMet[5]*TMath::DegToRad()); + if (denominator1 == 0.0 || denominator2 == 0.0) { + if (g_Verbosity >= c_Error) cout << m_XmlTag << ": Strip Intersection gives divide by zero - returning unrotated hit position" << endl; + double Xpos = m_YPitches[DetID]*((double)R_LVStrip.GetStripID() - ((m_NYStrips[DetID]-1)/2.0)); + double Ypos = m_XPitches[DetID]*((double)R_HVStrip.GetStripID() - ((m_NXStrips[DetID]-1)/2.0)); + return {Xpos, Ypos}; + } + + // Calculate XIntercept + 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())); + + // Solve for YIntercept + double YIntercept = (XIntercept - HVStripMet[0])*tan((HVStripMet[5]-90)*TMath::DegToRad()) + HVStripMet[1]; + + return {XIntercept, YIntercept}; +} + + +///////////////////////////////////////////////////////////////////////////////// + + 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"). // 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; + if (g_Verbosity >= c_Error) cout << m_XmlTag << "ERROR in MModuleDepthCalibration2024: HIT WITH NO STRIP HITS" << endl; return -1; } @@ -573,11 +734,17 @@ 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 ) { + if (g_Verbosity >= c_Error) cout << m_XmlTag << "ERROR in MModuleDepthCalibration2024: Depth Calibration: got NULL strip hit :( " << endl; + return -1; + } + if (SH->GetEnergy() == 0 ) { + if (g_Verbosity >= c_Error) cout << m_XmlTag << "ERROR in MModuleDepthCalibration2024: Depth Calibration: got strip without energy :( " << endl; + return -1; + } + if (SH->IsLowVoltageStrip()) { LVStrips.push_back(SH); LVStripIDs.push_back(SH->GetStripID()); } @@ -590,11 +757,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()); @@ -602,7 +769,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; } } @@ -613,48 +780,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; } @@ -667,9 +834,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 @@ -766,21 +935,23 @@ 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 - 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 { @@ -789,18 +960,22 @@ 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 - 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; @@ -809,6 +984,8 @@ 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 @@ -833,6 +1010,10 @@ TSpline3* MModuleDepthCalibration2024::GetSpline(int DetID, int Grade) } } + +///////////////////////////////////////////////////////////////////////////////// + + void MModuleDepthCalibration2024::ShowOptionsGUI() { // Show the options GUI - or do nothing @@ -842,23 +1023,36 @@ void MModuleDepthCalibration2024::ShowOptionsGUI() } +///////////////////////////////////////////////////////////////////////////////// + + 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 != nullptr) { + m_MaskMetrologyEnabled = (bool) MasKMetrologyNode->GetValueAsBoolean(); + } + + MXmlNode* MaskMetrologyFileNameNode = Node->GetNode("MaskMetrologyFileName"); + 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; @@ -874,7 +1068,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, "UCSDOverride",(unsigned int) m_UCSDOverride); + new MXmlNode(Node, "MaskMetrology", (bool)m_MaskMetrologyEnabled); + new MXmlNode(Node, "MaskMetrologyFileName", m_MaskMetrologyFile); + new MXmlNode(Node, "UCSDOverride", (bool)m_UCSDOverride); return Node; }