From 706c2e56de4a0b5036db08cc48feaaf4aeb289e9 Mon Sep 17 00:00:00 2001 From: ckierans Date: Mon, 13 Apr 2026 22:11:49 -0400 Subject: [PATCH 1/6] Pass through GR hits --- src/MModuleDepthCalibration.cxx | 312 ++++++++++++++++---------------- 1 file changed, 153 insertions(+), 159 deletions(-) diff --git a/src/MModuleDepthCalibration.cxx b/src/MModuleDepthCalibration.cxx index 3d5d3874..f629b2ed 100644 --- a/src/MModuleDepthCalibration.cxx +++ b/src/MModuleDepthCalibration.cxx @@ -177,202 +177,196 @@ void MModuleDepthCalibration::CreateExpos() bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) { - - 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 ){ - // Each event represents one photon. It contains Hits, representing interaction sites. - // H is a pointer to an instance of the MHit class. Each Hit has activated strips, represented by - // instances of the MStripHit class. - MHit* H = Event->GetHit(i); - - int Grade = GetHitGrade(H); - - // Handle different grades differently - // GRADE=-1 is an error. Break from the loop and continue. - if (Grade < 0){ - H->SetNoDepth(); - Event->SetDepthCalibrationError("Error in depth calibration"); - if (Grade == -1) { - ++m_ErrorSH; - } else if (Grade == -2) { - ++m_ErrorNullSH; - } else if (Grade == -3) { - ++m_ErrorNoE; - } - } else if (Grade > 4) { // GRADE=5 is some complicated geometry with multiple hits on a single strip. GRADE=6 means not all strips are adjacent. - H->SetNoDepth(); - Event->SetDepthCalibrationError("Multiple hits on single strip"); - if (Grade==5) { - ++m_Error5; - } else if (Grade==6) { - ++m_Error6; - } - } else { // If the Grade is 0-4, we can handle it. - - // 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) { - MStripHit* SH = H->GetStripHit(j); - if (SH->IsLowVoltageStrip()) LVStrips.push_back(SH); else HVStrips.push_back(SH); - } - double LVEnergyFraction; - double HVEnergyFraction; - MStripHit* LVSH = GetDominantStrip(LVStrips, LVEnergyFraction); - MStripHit* HVSH = GetDominantStrip(HVStrips, HVEnergyFraction); - - double CTD_s = 0.0; - - //now try and get z position - int DetID = LVSH->GetDetectorID(); - int LVStripID = LVSH->GetStripID(); - int HVStripID = HVSH->GetStripID(); - int PixelCode = 10000*DetID + 100*LVStripID + HVStripID; + for (unsigned int i = 0; i < Event->GetNHits(); ++i ){ + // H is a pointer to an instance of the MHit class. Each Hit has activated strips, represented by + // instances of the MStripHit class. + MHit* H = Event->GetHit(i); + + // Skip the depth calibration for Hits that have a GR Strip Hit + if (H->GetGuardRingHitFlag() == true) { + cout<<"Found GR Hit"<SetNoDepth(); + Event->SetDepthCalibrationError("Error in depth calibration"); + if (Grade == -1) { + ++m_ErrorSH; + } else if (Grade == -2) { + ++m_ErrorNullSH; + } else if (Grade == -3) { + ++m_ErrorNoE; + } + } else if (Grade > 4) { // GRADE=5 is some complicated geometry with multiple hits on a single strip. GRADE=6 means not all strips are adjacent. + H->SetNoDepth(); + Event->SetDepthCalibrationError("Multiple hits on single strip"); + if (Grade==5) { + ++m_Error5; + } else if (Grade==6) { + ++m_Error6; + } + } else { // If the Grade is 0-4, we can handle it. - 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()); + // 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; - // 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]; + for (unsigned int j = 0; j < H->GetNStripHits(); ++j) { + MStripHit* SH = H->GetStripHit(j); + if (SH->IsLowVoltageStrip()) LVStrips.push_back(SH); else HVStrips.push_back(SH); + } - } + double LVEnergyFraction; + double HVEnergyFraction; + MStripHit* LVSH = GetDominantStrip(LVStrips, LVEnergyFraction); + MStripHit* HVSH = GetDominantStrip(HVStrips, HVEnergyFraction); + + double CTD_s = 0.0; + + //now try and get z position + int DetID = LVSH->GetDetectorID(); + int LVStripID = LVSH->GetStripID(); + int HVStripID = HVSH->GetStripID(); + int PixelCode = 10000*DetID + 100*LVStripID + 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; + + 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. + // 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); + double Xsigma = m_YPitches[DetID]/sqrt(12.0); + double Ysigma = m_XPitches[DetID]/sqrt(12.0); + double Zsigma = m_Thicknesses[DetID]/sqrt(12.0); - vector* Coeffs = GetPixelCoeffs(PixelCode); + vector* Coeffs = GetPixelCoeffs(PixelCode); - vector CTDVec = GetCTD(DetID, Grade); - vector DepthVec = GetDepth(DetID); + 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(); + double LVTiming = LVSH->GetTiming(); + double HVTiming = HVSH->GetTiming(); - // If there aren't coefficients loaded, then report a depth calibration error. - if( Coeffs == nullptr ){ - // Set the bad flag for depth + // If there aren't coefficients loaded, then report a depth calibration error. + if( Coeffs == nullptr ){ + // Set the bad flag for depth + H->SetNoDepth(); + Event->SetDepthCalibrationError("No calibration coefficients"); + ++m_Error1; + } else if (CTDVec.size() == 0) { + if (g_Verbosity >= c_Error) cout << m_XmlTag << "Empty CTD vector" << endl; H->SetNoDepth(); Event->SetDepthCalibrationError("No calibration coefficients"); - ++m_Error1; - } else if (CTDVec.size() == 0) { - if (g_Verbosity >= c_Error) cout << m_XmlTag << "Empty CTD vector" << endl; - H->SetNoDepth(); - Event->SetDepthCalibrationError("No calibration coefficients"); - } else if (DepthVec.size() == 0) { - 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)) { - ++m_Error3; - H->SetNoDepth(); - Event->SetDepthCalibrationError("No timing"); - } else { + } else if (DepthVec.size() == 0) { + 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)) { + ++m_Error3; + H->SetNoDepth(); + Event->SetDepthCalibrationError("No timing"); + } else { - // If there are coefficients and timing information is loaded, try calculating the CTD and depth - double CTD = (HVTiming - LVTiming); + // 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 + // 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()); + 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))) { - H->SetNoDepth(); - Event->SetDepthCalibrationError("Out of Range"); - ++m_Error2; - } + //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))) { + H->SetNoDepth(); + Event->SetDepthCalibrationError("Out of Range"); + ++m_Error2; + } - // If the CTD is in range, calculate the depth - // Rather than plugging CTD into a spline to get depth, use the depth-CTD relation to calculate a probability-weighted depth value. - // This way we can avoid problems like non-monotonicity or assigning depth to events "outside" the detector - // Note that this requires that we don't massively overestimate the timing noise - else { - // Calculate the probability given timing noise of CTD_s corresponding to the values of depth in DepthVec - // Utlize symmetry of the normal distribution. - vector prob_dist = norm_pdf(CTDVec, CTD_s, noise/2.355); + // If the CTD is in range, calculate the depth + // Rather than plugging CTD into a spline to get depth, use the depth-CTD relation to calculate a probability-weighted depth value. + // This way we can avoid problems like non-monotonicity or assigning depth to events "outside" the detector + // Note that this requires that we don't massively overestimate the timing noise + else { + // Calculate the probability given timing noise of CTD_s corresponding to the values of depth in DepthVec + // Utlize symmetry of the normal distribution. + vector prob_dist = norm_pdf(CTDVec, CTD_s, noise/2.355); - // Weight the depth by probability - double prob_sum = 0.0; - for (unsigned int k=0; k < prob_dist.size(); ++k) { - prob_sum += prob_dist[k]; - } - double weighted_depth = 0.0; - - for (unsigned int k = 0; k < DepthVec.size(); ++k) { - weighted_depth += prob_dist[k] * DepthVec[k]; - } + // Weight the depth by probability + double prob_sum = 0.0; + for (unsigned int k=0; k < prob_dist.size(); ++k) { + prob_sum += prob_dist[k]; + } + double weighted_depth = 0.0; + + for (unsigned int k = 0; k < DepthVec.size(); ++k) { + weighted_depth += prob_dist[k] * DepthVec[k]; + } - // Calculate the expectation value of the depth - double mean_depth = weighted_depth/prob_sum; + // Calculate the expectation value of the depth + double mean_depth = weighted_depth/prob_sum; - // Calculate the standard deviation of the depth - double depth_var = 0.0; + // Calculate the standard deviation of the depth + double depth_var = 0.0; - for (unsigned int k=0; k < DepthVec.size(); ++k) { - depth_var += prob_dist[k] * pow(DepthVec[k] - mean_depth, 2.0); - } + for (unsigned int k=0; k < DepthVec.size(); ++k) { + depth_var += prob_dist[k] * pow(DepthVec[k] - mean_depth, 2.0); + } - Zsigma = sqrt(depth_var/prob_sum); - Zpos = mean_depth; - // Zpos = mean_depth - (m_Thicknesses[DetID]/2.0); + 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->HasStripPairingError()==false) { - if (HasExpos() == true) { - m_ExpoDepthCalibration->AddDepth(DetID, Zpos); - } + // Add the depth to the GUI histogram. + if (Event->HasStripPairingError()==false) { + if (HasExpos() == true) { + m_ExpoDepthCalibration->AddDepth(DetID, Zpos); } - m_NoError+=1; } + m_NoError+=1; } + } - if (g_Verbosity >= c_Info) cout << m_XmlTag << "Strip ID :" << LVStripID << " " << HVStripID << endl << "Hit position: "<< Xpos << " " << Ypos << " " << Zpos << endl; + if (g_Verbosity >= c_Info) cout << m_XmlTag << "Strip ID :" << LVStripID << " " << HVStripID << endl << "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); + 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. - MVector PositionResolution(Xsigma, Ysigma, Zsigma); - MVector GlobalResolution = ((m_Detectors[DetID]->GetSensitiveVolume(0)->GetPositionInWorldVolume(PositionResolution)) - (m_Detectors[DetID]->GetSensitiveVolume(0)->GetPositionInWorldVolume(LocalOrigin))).Abs(); + // Make sure XYZ resolution are correctly mapped to the global coord system. + MVector PositionResolution(Xsigma, Ysigma, Zsigma); + MVector GlobalResolution = ((m_Detectors[DetID]->GetSensitiveVolume(0)->GetPositionInWorldVolume(PositionResolution)) - (m_Detectors[DetID]->GetSensitiveVolume(0)->GetPositionInWorldVolume(LocalOrigin))).Abs(); - H->SetPosition(GlobalPosition); + H->SetPosition(GlobalPosition); - H->SetPositionResolution(GlobalResolution); + H->SetPositionResolution(GlobalResolution); - } } } From adfaf1c6817fe67934f836467c879cfcad3fd4be Mon Sep 17 00:00:00 2001 From: ckierans Date: Tue, 14 Apr 2026 01:04:54 -0400 Subject: [PATCH 2/6] Rebased with current develop/em Included option to save QA-flagged events --- src/MModuleDepthCalibration.cxx | 1 - src/MReadOutAssembly.cxx | 1 - 2 files changed, 2 deletions(-) diff --git a/src/MModuleDepthCalibration.cxx b/src/MModuleDepthCalibration.cxx index f629b2ed..a89f0a2d 100644 --- a/src/MModuleDepthCalibration.cxx +++ b/src/MModuleDepthCalibration.cxx @@ -185,7 +185,6 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) // Skip the depth calibration for Hits that have a GR Strip Hit if (H->GetGuardRingHitFlag() == true) { - cout<<"Found GR Hit"<GetGuardRingHitFlag() == true) { From c24de0822ee7b3de4f7948c246646d0fafd0e25e Mon Sep 17 00:00:00 2001 From: ckierans Date: Thu, 7 May 2026 21:30:41 -0400 Subject: [PATCH 3/6] GR hits now passed thru revan Continue rebase with develop/em --- include/MModuleDepthCalibration.h | 1 + src/MHit.cxx | 90 +++++++++++-------- src/MModuleDepthCalibration.cxx | 78 +++++++++++----- ...MModuleStripPairingMultiRoundChiSquare.cxx | 3 - src/MReadOutAssembly.cxx | 8 +- 5 files changed, 107 insertions(+), 73 deletions(-) diff --git a/include/MModuleDepthCalibration.h b/include/MModuleDepthCalibration.h index 0926a61f..6b6ed684 100644 --- a/include/MModuleDepthCalibration.h +++ b/include/MModuleDepthCalibration.h @@ -184,6 +184,7 @@ class MModuleDepthCalibration : public MModule uint64_t m_ErrorNullSH; uint64_t m_ErrorNoE; unordered_map m_Detectors; + unordered_map m_GRDetectors; vector m_DetectorIDs; MModuleEnergyCalibration* m_EnergyCalibration; MGUIExpoDepthCalibration* m_ExpoDepthCalibration; diff --git a/src/MHit.cxx b/src/MHit.cxx index 260352c4..a6ae19e9 100644 --- a/src/MHit.cxx +++ b/src/MHit.cxx @@ -193,54 +193,66 @@ void MHit::StreamEvta(ostream& S) { // Stream the hit in MEGAlib's EVTA format - // Assemble the origin information - vector Origins; - - // Only origins existing on both low-voltage and high-voltage strips count - vector LVOrigins; - vector HVOrigins; - for (unsigned int s = 0; s < GetNStripHits(); ++s) { - MStripHit* StripHit = m_StripHits[s]; - vector NewOrigins = StripHit->GetOrigins(); - if (StripHit->IsLowVoltageStrip() == true) { - for (int o: NewOrigins) { - LVOrigins.push_back(o); - } - } else { - for (int o: NewOrigins) { - HVOrigins.push_back(o); - } - } - } - - sort(LVOrigins.begin(), LVOrigins.end()); - LVOrigins.erase(unique(LVOrigins.begin(), LVOrigins.end()), LVOrigins.end()); - sort(HVOrigins.begin(), HVOrigins.end()); - HVOrigins.erase(unique(HVOrigins.begin(), HVOrigins.end()), HVOrigins.end()); + if ((m_GuardRingHit == false) && (m_NoDepth == false)) { - set_intersection(LVOrigins.begin(), LVOrigins.end(), - HVOrigins.begin(), HVOrigins.end(), - std::back_inserter(Origins)); + // Assemble the origin information + vector Origins; - if ((LVOrigins.size() != 0 || HVOrigins.size() != 0) && Origins.size() == 0) { - // If strip pairing mixed the hits completely, keep the mixed origin information + // Only origins existing on both low-voltage and high-voltage strips count + vector LVOrigins; + vector HVOrigins; for (unsigned int s = 0; s < GetNStripHits(); ++s) { MStripHit* StripHit = m_StripHits[s]; vector NewOrigins = StripHit->GetOrigins(); - for (int o: NewOrigins) { - Origins.push_back(o); + if (StripHit->IsLowVoltageStrip() == true) { + for (int o: NewOrigins) { + LVOrigins.push_back(o); + } + } else { + for (int o: NewOrigins) { + HVOrigins.push_back(o); + } } } - sort(Origins.begin(), Origins.end()); - Origins.erase(unique(Origins.begin(), Origins.end()), Origins.end()); - } - S<<"HT 3;"< NewOrigins = StripHit->GetOrigins(); + for (int o: NewOrigins) { + Origins.push_back(o); + } + } + sort(Origins.begin(), Origins.end()); + Origins.erase(unique(Origins.begin(), Origins.end()), Origins.end()); + } + + S<<"HT 3;"<GetHit(i); - - // Skip the depth calibration for Hits that have a GR Strip Hit + + //Initalize position variables: + double Xpos = 0; + double Ypos = 0; + double Zpos = 0; + double Xsigma = 0; + double Ysigma = 0; + double Zsigma = 0; + + + // First, check for GR Hits if (H->GetGuardRingHitFlag() == true) { + // Skip the depth calibration for Hits that have a GR Strip Hit. Define position + // to be anywhere within the GR volume to pass through to revan + + int GRDetID = H->GetStripHit(0)->GetDetectorID(); + // Find unique/random position within the GR volume to assign as the hit position, as revan expects + MVector GRPosition = m_Geometry->GetDetector(m_GRDetectors[GRDetID]->GetName())->GetSensitiveVolume(0)->GetRandomPositionExclusivelyInside(); + + Xpos = GRPosition[0]; + Ypos = GRPosition[1]; + Zpos = GRPosition[2]; + + if (g_Verbosity >= c_Info) cout << m_XmlTag << "GR Hit :" << GRDetID << ", " << "set hit position: "<< Xpos << " " << Ypos << " " << Zpos << endl; + + MVector GlobalPositionGR = m_GRDetectors[GRDetID]->GetSensitiveVolume(0)->GetPositionInWorldVolume(GRPosition); + MVector GRResolution(0.1, 0.1, 0.1); + + H->SetPosition(GlobalPositionGR); + H->SetPositionResolution(GRResolution); + + // Skip past the rest of the calibration and move to next Hit continue; - } + + } + // If not a GR Hit, perform the depth/position calibration... - int Grade = GetHitGrade(H); - // Handle different grades differently + //TODO Rename Grade variables + // The event "Grade" is the sub-pixel region determined via charge sharing + int Grade = GetHitGrade(H); // GRADE=-1 is an error. Break from the loop and continue. if (Grade < 0){ H->SetNoDepth(); @@ -212,7 +242,7 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) } } else { // If the Grade is 0-4, we can handle it. - // 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; @@ -227,20 +257,19 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) MStripHit* LVSH = GetDominantStrip(LVStrips, LVEnergyFraction); MStripHit* HVSH = GetDominantStrip(HVStrips, HVEnergyFraction); - double CTD_s = 0.0; - - //now try and get z position int DetID = LVSH->GetDetectorID(); int LVStripID = LVSH->GetStripID(); int HVStripID = HVSH->GetStripID(); int PixelCode = 10000*DetID + 100*LVStripID + HVStripID; - //Define the X/Y positions based on the detector pitch and number of strip hits + + // TODO: Calculate X and Y positions more rigorously using charge sharing. + + // 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; + Xpos = m_YPitches[DetID]*((double)LVStripID - ((m_NYStrips[DetID]-1)/2.0)); + Ypos = m_XPitches[DetID]*((double)HVStripID - ((m_NXStrips[DetID]-1)/2.0)); 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 @@ -253,15 +282,15 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) Ypos = inter[1]; } + Xsigma = m_YPitches[DetID]/sqrt(12.0); + Ysigma = m_XPitches[DetID]/sqrt(12.0); + Zsigma = m_Thicknesses[DetID]/sqrt(12.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); - double Zsigma = m_Thicknesses[DetID]/sqrt(12.0); + // Now try and get z position + double CTD_s = 0.0; vector* Coeffs = GetPixelCoeffs(PixelCode); - vector CTDVec = GetCTD(DetID, Grade); vector DepthVec = GetDepth(DetID); @@ -290,8 +319,6 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) // 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()); @@ -361,7 +388,6 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) MVector GlobalResolution = ((m_Detectors[DetID]->GetSensitiveVolume(0)->GetPositionInWorldVolume(PositionResolution)) - (m_Detectors[DetID]->GetSensitiveVolume(0)->GetPositionInWorldVolume(LocalOrigin))).Abs(); H->SetPosition(GlobalPosition); - H->SetPositionResolution(GlobalResolution); @@ -507,6 +533,9 @@ bool MModuleDepthCalibration::LoadDetectorDimensions(MDGeometryQuest* Geometry) } +///////////////////////////////////////////////////////////////////////////////// + + bool MModuleDepthCalibration::LoadCoeffsFile(MString FileName) { // Read in the stretch and offset file, which should have a header line with information on the measurements: @@ -1139,6 +1168,7 @@ void MModuleDepthCalibration::Finalize() m_XPitches.clear(); m_YPitches.clear(); m_Detectors.clear(); + m_GRDetectors.clear(); m_CTDMap.clear(); m_DepthGrid.clear(); m_SplineMap.clear(); diff --git a/src/MModuleStripPairingMultiRoundChiSquare.cxx b/src/MModuleStripPairingMultiRoundChiSquare.cxx index 3b6b00c4..8e9913fa 100644 --- a/src/MModuleStripPairingMultiRoundChiSquare.cxx +++ b/src/MModuleStripPairingMultiRoundChiSquare.cxx @@ -833,9 +833,6 @@ bool MModuleStripPairingMultiRoundChiSquare::AnalyzeEvent(MReadOutAssembly* Even Event->GetHit(h)->SetGuardRingHitFlag(true); } } - if (Event->GetHit(h)->GetGuardRingHitFlag() == true) { - Event->SetStripPairing_QualityFlag("GR Hit: Detector ID " + to_string(d) + " and Energy " + to_string(Event->GetHit(h)->GetEnergy())); - } } } // End Detector loop diff --git a/src/MReadOutAssembly.cxx b/src/MReadOutAssembly.cxx index 3ec5170f..8c5e6953 100644 --- a/src/MReadOutAssembly.cxx +++ b/src/MReadOutAssembly.cxx @@ -632,13 +632,7 @@ void MReadOutAssembly::StreamEvta(ostream& S) } for (unsigned int h = 0; h < m_Hits.size(); ++h) { - // Don't print Guard Ring hits as normal strip hits as they don't have positions defined - // the corresponding energy is saved in the StripPairing QA message - if (m_Hits[h]->GetGuardRingHitFlag() == true) { - continue; - } else { - m_Hits[h]->StreamEvta(S); - } + m_Hits[h]->StreamEvta(S); } S<<"CC NStripHits "< Date: Mon, 18 May 2026 23:18:31 -0400 Subject: [PATCH 4/6] Include GR and XE hits in evta Continue rebasing --- src/MHit.cxx | 61 +++++++++++++++- src/MModuleDepthCalibration.cxx | 126 +++++++++++++++++--------------- 2 files changed, 125 insertions(+), 62 deletions(-) diff --git a/src/MHit.cxx b/src/MHit.cxx index a6ae19e9..726ad246 100644 --- a/src/MHit.cxx +++ b/src/MHit.cxx @@ -244,11 +244,10 @@ void MHit::StreamEvta(ostream& S) } S< tokens = Line.Tokenize(" "); + if( tokens.size() >= 5 ){ + m_Position.SetX( tokens.at(1).ToDouble() ); + m_Position.SetY( tokens.at(2).ToDouble() ); + m_Position.SetZ( tokens.at(3).ToDouble() ); + m_Energy = tokens.at(4).ToDouble(); + return true; + } else { + return false; + } + } else { + return false; + } + */ + +} + + +//////////////////////////////////////////////////////////////////////////////// + + +//! Set the origins from the simulations (take care of duplicates) +void MHit::AddOrigins(vector Origins) +>>>>>>> 5a54288 (Include GR and XE hits in evta) { // Parse a hit in Nuclearizer's DAT format diff --git a/src/MModuleDepthCalibration.cxx b/src/MModuleDepthCalibration.cxx index ce3b2587..9f0106f4 100644 --- a/src/MModuleDepthCalibration.cxx +++ b/src/MModuleDepthCalibration.cxx @@ -116,6 +116,7 @@ bool MModuleDepthCalibration::Initialize() return false; } + // Check for issues with geometry or file loading, and return appropriate errors if (m_DetectorIDs.size() == 0) { cout<<"No Strip3D detectors were found."<GetNHits(); ++i ){ + for (unsigned int i = 0; i < Event->GetNHits(); ++i ) { // H is a pointer to an instance of the MHit class. Each Hit has activated strips, represented by // instances of the MStripHit class. MHit* H = Event->GetHit(i); + int DetID = H->GetStripHit(0)->GetDetectorID(); //Initalize position variables: - double Xpos = 0; - double Ypos = 0; - double Zpos = 0; - double Xsigma = 0; - double Ysigma = 0; - double Zsigma = 0; + double Xpos = 0.0, Ypos = 0.0, Zpos = 0.0; + double Xsigma = 0.0, Ysigma = 0.0, Zsigma = 0.0; // First, check for GR Hits if (H->GetGuardRingHitFlag() == true) { - // Skip the depth calibration for Hits that have a GR Strip Hit. Define position - // to be anywhere within the GR volume to pass through to revan + // For GR Hit, define position to be anywhere within the GR volume to pass through to revan int GRDetID = H->GetStripHit(0)->GetDetectorID(); // Find unique/random position within the GR volume to assign as the hit position, as revan expects @@ -206,19 +203,15 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) if (g_Verbosity >= c_Info) cout << m_XmlTag << "GR Hit :" << GRDetID << ", " << "set hit position: "<< Xpos << " " << Ypos << " " << Zpos << endl; MVector GlobalPositionGR = m_GRDetectors[GRDetID]->GetSensitiveVolume(0)->GetPositionInWorldVolume(GRPosition); - MVector GRResolution(0.1, 0.1, 0.1); - H->SetPosition(GlobalPositionGR); - H->SetPositionResolution(GRResolution); // Skip past the rest of the calibration and move to next Hit continue; - } - // If not a GR Hit, perform the depth/position calibration... + } // If not a GR Hit, perform the depth/position calibration... - //TODO Rename Grade variables + //TODO: Rename Grade variables // The event "Grade" is the sub-pixel region determined via charge sharing int Grade = GetHitGrade(H); // GRADE=-1 is an error. Break from the loop and continue. @@ -234,7 +227,10 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) } } else if (Grade > 4) { // GRADE=5 is some complicated geometry with multiple hits on a single strip. GRADE=6 means not all strips are adjacent. H->SetNoDepth(); - Event->SetDepthCalibrationError("Multiple hits on single strip"); + if (Event->HasDepthCalibrationError() == false) { + // Check if the DepthCalibrationError is already define for the Event before duplciating it + Event->SetDepthCalibrationError("Multiple hits on single strip"); + } if (Grade==5) { ++m_Error5; } else if (Grade==6) { @@ -257,7 +253,6 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) MStripHit* LVSH = GetDominantStrip(LVStrips, LVEnergyFraction); MStripHit* HVSH = GetDominantStrip(HVStrips, HVEnergyFraction); - int DetID = LVSH->GetDetectorID(); int LVStripID = LVSH->GetStripID(); int HVStripID = HVSH->GetStripID(); int PixelCode = 10000*DetID + 100*LVStripID + HVStripID; @@ -298,7 +293,7 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) double HVTiming = HVSH->GetTiming(); // If there aren't coefficients loaded, then report a depth calibration error. - if( Coeffs == nullptr ){ + if (Coeffs == nullptr) { // Set the bad flag for depth H->SetNoDepth(); Event->SetDepthCalibrationError("No calibration coefficients"); @@ -321,64 +316,73 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) double CTD = (HVTiming - LVTiming); 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 CTDmin = * std::min_element(CTDVec.begin(), CTDVec.end()); + double CTDmax = * std::max_element(CTDVec.begin(), CTDVec.end()); 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 the CTD is not crazy out of range, stick it on the edge of the detector and report a QA flag: + if ((CTD_s < CTDmin) && (CTD_s > (CTDmin - 5.0*noise))) { + CTD_s = CTDmin; + Event->SetDepthCalibrationError("Forced to edge"); + } else if ((CTD_s > CTDmax) && (CTD_s < (CTDmax + 5.0*noise))) { + CTD_s = CTDmax; + Event->SetDepthCalibrationError("Forced to edge"); + } else if ((CTD_s < (CTDmin - 5.0*noise)) || (CTD_s > (CTDmax + 5.0*noise))) { + // If the CTD is >5 sigma away from the max/min CTD, then don't calibrate this hit and return an error H->SetNoDepth(); Event->SetDepthCalibrationError("Out of Range"); ++m_Error2; } - // If the CTD is in range, calculate the depth + // Now, calculate the depth // Rather than plugging CTD into a spline to get depth, use the depth-CTD relation to calculate a probability-weighted depth value. // This way we can avoid problems like non-monotonicity or assigning depth to events "outside" the detector // Note that this requires that we don't massively overestimate the timing noise - else { - // Calculate the probability given timing noise of CTD_s corresponding to the values of depth in DepthVec - // Utlize symmetry of the normal distribution. - vector prob_dist = norm_pdf(CTDVec, CTD_s, noise/2.355); + // Calculate the probability given timing noise of CTD_s corresponding to the values of depth in DepthVec + // Utlize symmetry of the normal distribution. + vector prob_dist = norm_pdf(CTDVec, CTD_s, noise/2.355); - // Weight the depth by probability - double prob_sum = 0.0; - for (unsigned int k=0; k < prob_dist.size(); ++k) { - prob_sum += prob_dist[k]; - } - double weighted_depth = 0.0; + // Weight the depth by probability + double prob_sum = 0.0; + for (unsigned int k=0; k < prob_dist.size(); ++k) { + prob_sum += prob_dist[k]; + } + double weighted_depth = 0.0; - for (unsigned int k = 0; k < DepthVec.size(); ++k) { - weighted_depth += prob_dist[k] * DepthVec[k]; - } + 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 expectation value of the depth + double mean_depth = weighted_depth/prob_sum; - // Calculate the standard deviation of the depth - double depth_var = 0.0; + // Calculate the standard deviation of the depth + double depth_var = 0.0; - for (unsigned int k=0; k < DepthVec.size(); ++k) { - depth_var += prob_dist[k] * pow(DepthVec[k] - mean_depth, 2.0); - } + for (unsigned int k=0; k < DepthVec.size(); ++k) { + depth_var += prob_dist[k] * pow(DepthVec[k] - mean_depth, 2.0); + } - Zsigma = sqrt(depth_var/prob_sum); - Zpos = mean_depth; - // Zpos = mean_depth - (m_Thicknesses[DetID]/2.0); + 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->HasStripPairingError()==false) { - if (HasExpos() == true) { - m_ExpoDepthCalibration->AddDepth(DetID, Zpos); - } + // Add the depth to the GUI histogram. + if (Event->HasStripPairingError()==false) { + if (HasExpos() == true) { + m_ExpoDepthCalibration->AddDepth(DetID, Zpos); } - m_NoError+=1; } + + m_NoError+=1; + } + + if (g_Verbosity >= c_Info) cout << m_XmlTag << "Strip ID :" << LVStripID << " " << HVStripID << endl << "Hit position: "<< Xpos << " " << Ypos << " " << Zpos << endl; - if (g_Verbosity >= c_Info) cout << m_XmlTag << "Strip ID :" << LVStripID << " " << HVStripID << endl << "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); @@ -390,14 +394,20 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) H->SetPosition(GlobalPosition); H->SetPositionResolution(GlobalResolution); - - + // For events that have NoDepth, these can be passed through revan with an XE flag and any position within the detector + if (H->GetNoDepth() == true) { + DetID = H->GetStripHit(0)->GetDetectorID(); + MVector XEPosition = m_Geometry->GetDetector(m_Detectors[DetID]->GetName())->GetSensitiveVolume(0)->GetRandomPositionExclusivelyInside(); + + MVector GlobalPositionXE = m_Detectors[DetID]->GetSensitiveVolume(0)->GetPositionInWorldVolume(XEPosition); + H->SetPosition(GlobalPositionXE); } - } + } + Event->SetAnalysisProgress(MAssembly::c_DepthCorrection | MAssembly::c_PositionDetermiation); - return true; + } From 1acb1e5f4d23ab964bee0b84c2a7f236894d3407 Mon Sep 17 00:00:00 2001 From: ckierans Date: Tue, 11 Aug 2026 01:20:54 -0400 Subject: [PATCH 5/6] Rebase error, now compiles --- src/MHit.cxx | 54 --------------------------------- src/MModuleDepthCalibration.cxx | 37 +++++++++++++++------- 2 files changed, 26 insertions(+), 65 deletions(-) diff --git a/src/MHit.cxx b/src/MHit.cxx index 726ad246..3af65968 100644 --- a/src/MHit.cxx +++ b/src/MHit.cxx @@ -259,60 +259,6 @@ void MHit::StreamEvta(ostream& S) bool MHit::Parse(MString& Line, int Version) - - //check that line begins with HT - const char* line = Line.Data(); - - if( line[0] == 'H' && line[1] == 'T' ) { - float X,Y,Z,E; - sscanf(&line[3],"%f %f %f %f",&X, &Y, &Z, &E); - m_Position.SetX(X); - m_Position.SetY(Y); - m_Position.SetZ(Z); - m_Energy = E; - return true; - - } else if (line[0] == 'G' && line[1] == 'R' ) { - float X,Y,Z,E; - sscanf(&line[2],"%f %f %f %f",&X, &Y, &Z, &E); - m_Position.SetX(X); - m_Position.SetY(Y); - m_Position.SetZ(Z); - m_Energy = E; - m_GuardRingHit = true; - return true; - - } else { - return false; - } - - - /* - if( Line.BeginsWith("HT") ){ - vector tokens = Line.Tokenize(" "); - if( tokens.size() >= 5 ){ - m_Position.SetX( tokens.at(1).ToDouble() ); - m_Position.SetY( tokens.at(2).ToDouble() ); - m_Position.SetZ( tokens.at(3).ToDouble() ); - m_Energy = tokens.at(4).ToDouble(); - return true; - } else { - return false; - } - } else { - return false; - } - */ - -} - - -//////////////////////////////////////////////////////////////////////////////// - - -//! Set the origins from the simulations (take care of duplicates) -void MHit::AddOrigins(vector Origins) ->>>>>>> 5a54288 (Include GR and XE hits in evta) { // Parse a hit in Nuclearizer's DAT format diff --git a/src/MModuleDepthCalibration.cxx b/src/MModuleDepthCalibration.cxx index 9f0106f4..96be3bbe 100644 --- a/src/MModuleDepthCalibration.cxx +++ b/src/MModuleDepthCalibration.cxx @@ -116,7 +116,6 @@ bool MModuleDepthCalibration::Initialize() return false; } - // Check for issues with geometry or file loading, and return appropriate errors if (m_DetectorIDs.size() == 0) { cout<<"No Strip3D detectors were found."<GetGuardRingHitFlag() == true) { // For GR Hit, define position to be anywhere within the GR volume to pass through to revan - int GRDetID = H->GetStripHit(0)->GetDetectorID(); +// int GRDetID = H->GetStripHit(0)->GetDetectorID(); // Find unique/random position within the GR volume to assign as the hit position, as revan expects - MVector GRPosition = m_Geometry->GetDetector(m_GRDetectors[GRDetID]->GetName())->GetSensitiveVolume(0)->GetRandomPositionExclusivelyInside(); +// MVector GRPosition = m_Geometry->GetDetector(m_GRDetectors[GRDetID]->GetName())->GetSensitiveVolume(0)->GetRandomPositionExclusivelyInside(); - Xpos = GRPosition[0]; - Ypos = GRPosition[1]; - Zpos = GRPosition[2]; +// Xpos = GRPosition[0]; +// Ypos = GRPosition[1]; +// Zpos = GRPosition[2]; - if (g_Verbosity >= c_Info) cout << m_XmlTag << "GR Hit :" << GRDetID << ", " << "set hit position: "<< Xpos << " " << Ypos << " " << Zpos << endl; +// if (g_Verbosity >= c_Info) cout << m_XmlTag << "GR Hit :" << GRDetID << ", " << "set hit position: "<< Xpos << " " << Ypos << " " << Zpos << endl; - MVector GlobalPositionGR = m_GRDetectors[GRDetID]->GetSensitiveVolume(0)->GetPositionInWorldVolume(GRPosition); - H->SetPosition(GlobalPositionGR); +// MVector GlobalPositionGR = m_GRDetectors[GRDetID]->GetSensitiveVolume(0)->GetPositionInWorldVolume(GRPosition); +// H->SetPosition(GlobalPositionGR); // Skip past the rest of the calibration and move to next Hit - continue; + // continue; } // If not a GR Hit, perform the depth/position calibration... @@ -536,9 +535,25 @@ bool MModuleDepthCalibration::LoadDetectorDimensions(MDGeometryQuest* Geometry) cout<<"ERROR in MModuleDepthCalibration::Initialize: Found a Strip3D detector with "<GetNSensitiveVolumes()<<" Sensitive Volumes."<GetTypeName() == "Simple") { +// if (det->GetNSensitiveVolumes() == 1) { +// MString DetectorName = det->GetName(); +// string DetName = DetectorName.GetString(); + + // Check that the DetID agrees with the naming scheme GeD_X +// if (DetectorName.BeginsWith("GuardRingDetector_GeD_") == true) { +// DetectorName.RemoveAllInPlace("GuardRingDetector_GeD_"); // The number after GeD is the COSI detector ID +// if (DetID != (DetectorName.ToUnsignedInt()-1)) { // The GR detector ID is +1 compared to the GeD detector for the same DetID. +// if (g_Verbosity >= c_Error) { +// cout << "ERROR in MModuleDepthCalibration::Initialize: Non-matching DetID="< Date: Tue, 11 Aug 2026 09:35:51 -0400 Subject: [PATCH 6/6] Correcting GRHit position definition --- src/MModuleDepthCalibration.cxx | 84 ++++++++++++++++++--------------- 1 file changed, 46 insertions(+), 38 deletions(-) diff --git a/src/MModuleDepthCalibration.cxx b/src/MModuleDepthCalibration.cxx index 96be3bbe..e0e38868 100644 --- a/src/MModuleDepthCalibration.cxx +++ b/src/MModuleDepthCalibration.cxx @@ -117,7 +117,7 @@ bool MModuleDepthCalibration::Initialize() } if (m_DetectorIDs.size() == 0) { - cout<<"No Strip3D detectors were found."<= c_Error) cout << m_XmlTag << ": No Strip3D detectors were found"<= c_Info) cout << m_XmlTag << ": !!! Mask Metrology Enabled !!!" << endl; m_MaskMetrologyFileIsLoaded = LoadMaskMetrologyFile(m_MaskMetrologyFileName); if (m_MaskMetrologyFileIsLoaded == false) { - if (g_Verbosity >= c_Error) cout << m_XmlTag << "Unable to open Metrology file" << endl; + if (g_Verbosity >= c_Error) cout << m_XmlTag << ": Unable to open Metrology file" << endl; return false; } } @@ -191,21 +191,28 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) if (H->GetGuardRingHitFlag() == true) { // For GR Hit, define position to be anywhere within the GR volume to pass through to revan -// int GRDetID = H->GetStripHit(0)->GetDetectorID(); - // Find unique/random position within the GR volume to assign as the hit position, as revan expects -// MVector GRPosition = m_Geometry->GetDetector(m_GRDetectors[GRDetID]->GetName())->GetSensitiveVolume(0)->GetRandomPositionExclusivelyInside(); + int GRDetID = H->GetStripHit(0)->GetDetectorID(); -// Xpos = GRPosition[0]; -// Ypos = GRPosition[1]; -// Zpos = GRPosition[2]; + // Find unique/random position within the GR volume to assign as the hit position, as revan expects + if (m_GRDetectors[GRDetID]->GetName().BeginsWith("GuardRing") == true) { + MVector GRPosition = m_Geometry->GetDetector(m_GRDetectors[GRDetID]->GetName())->GetSensitiveVolume(0)->GetRandomPositionExclusivelyInside(); + + Xpos = GRPosition[0]; + Ypos = GRPosition[1]; + Zpos = GRPosition[2]; + + if (g_Verbosity >= c_Info) cout << m_XmlTag << ": GR Hit: Det ID " << GRDetID << ", " << "set hit position: "<< Xpos << " " << Ypos << " " << Zpos << endl; -// if (g_Verbosity >= c_Info) cout << m_XmlTag << "GR Hit :" << GRDetID << ", " << "set hit position: "<< Xpos << " " << Ypos << " " << Zpos << endl; + MVector GlobalPositionGR = m_GRDetectors[GRDetID]->GetSensitiveVolume(0)->GetPositionInWorldVolume(GRPosition); + H->SetPosition(GlobalPositionGR); -// MVector GlobalPositionGR = m_GRDetectors[GRDetID]->GetSensitiveVolume(0)->GetPositionInWorldVolume(GRPosition); -// H->SetPosition(GlobalPositionGR); + } else { + if (g_Verbosity >= c_Error) cout << m_XmlTag << ": Could not find GuardRing volume for position determination" << endl; + } + // Skip past the rest of the calibration and move to next Hit - // continue; + continue; } // If not a GR Hit, perform the depth/position calibration... @@ -298,11 +305,11 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) Event->SetDepthCalibrationError("No calibration coefficients"); ++m_Error1; } else if (CTDVec.size() == 0) { - if (g_Verbosity >= c_Error) cout << m_XmlTag << "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) { - if (g_Verbosity >= c_Error) cout << m_XmlTag << "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)) { @@ -378,7 +385,7 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) } - if (g_Verbosity >= c_Info) cout << m_XmlTag << "Strip ID :" << LVStripID << " " << HVStripID << endl << "Hit position: "<< Xpos << " " << Ypos << " " << Zpos << endl; + if (g_Verbosity >= c_Info) cout << m_XmlTag << ": Strip ID: " << LVStripID << " " << HVStripID << endl << "Hit position: "<< Xpos << " " << Ypos << " " << Zpos << endl; } @@ -393,8 +400,8 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) H->SetPosition(GlobalPosition); H->SetPositionResolution(GlobalResolution); - // For events that have NoDepth, these can be passed through revan with an XE flag and any position within the detector - if (H->GetNoDepth() == true) { + // For non-GR hits that have NoDepth, these can be passed through revan with an XE flag and any position within the detector + if (H->GetNoDepth() == true && H->GetGuardRingHitFlag() == false) { DetID = H->GetStripHit(0)->GetDetectorID(); MVector XEPosition = m_Geometry->GetDetector(m_Detectors[DetID]->GetName())->GetSensitiveVolume(0)->GetRandomPositionExclusivelyInside(); @@ -535,23 +542,24 @@ bool MModuleDepthCalibration::LoadDetectorDimensions(MDGeometryQuest* Geometry) cout<<"ERROR in MModuleDepthCalibration::Initialize: Found a Strip3D detector with "<GetNSensitiveVolumes()<<" Sensitive Volumes."<GetTypeName() == "Simple") { -// if (det->GetNSensitiveVolumes() == 1) { -// MString DetectorName = det->GetName(); -// string DetName = DetectorName.GetString(); - - // Check that the DetID agrees with the naming scheme GeD_X -// if (DetectorName.BeginsWith("GuardRingDetector_GeD_") == true) { -// DetectorName.RemoveAllInPlace("GuardRingDetector_GeD_"); // The number after GeD is the COSI detector ID -// if (DetID != (DetectorName.ToUnsignedInt()-1)) { // The GR detector ID is +1 compared to the GeD detector for the same DetID. -// if (g_Verbosity >= c_Error) { -// cout << "ERROR in MModuleDepthCalibration::Initialize: Non-matching DetID="<GetTypeName() == "Simple" || det->GetTypeName() == "Scintillator") { + if (det->GetNSensitiveVolumes() == 1) { + MString DetectorName = det->GetName(); + string DetName = DetectorName.GetString(); + // Check that the DetID agrees with the naming scheme GuardRingDetector_GeD_ + if (DetectorName.BeginsWith("GuardRingDetector_GeD_") == true) { + DetectorName.RemoveAllInPlace("GuardRingDetector_GeD_"); // The number after GeD is the COSI detector ID + if (DetID != (DetectorName.ToUnsignedInt()-1)) { // The GR detector ID is +1 compared to the GeD detector for the same DetID. + if (g_Verbosity >= c_Error) { + cout << "ERROR in MModuleDepthCalibration::Initialize: Non-matching DetID="< 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; + cout << m_XmlTag << ": The stretch and offset were calculated for " << m_Coeffs_Energy << " keV." << endl; } } else { std::vector Tokens = Line.Tokenize(","); @@ -807,7 +815,7 @@ int MModuleDepthCalibration::GetHitGrade(MHit* H){ } if (H->GetNStripHits() == 0) { // Error if no strip hits listed. Bad grade is returned - if (g_Verbosity >= c_Error) cout << m_XmlTag << "ERROR in MModuleDepthCalibration: HIT WITH NO STRIP HITS" << endl; + if (g_Verbosity >= c_Error) cout << m_XmlTag << ": ERROR in MModuleDepthCalibration: HIT WITH NO STRIP HITS" << endl; return -1; } @@ -819,11 +827,11 @@ int MModuleDepthCalibration::GetHitGrade(MHit* H){ for (unsigned int j = 0; j < H->GetNStripHits(); ++j) { MStripHit* SH = H->GetStripHit(j); if (SH == nullptr ) { - if (g_Verbosity >= c_Error) cout << m_XmlTag << "ERROR in MModuleDepthCalibration: Depth Calibration: got NULL strip hit :( " << endl; + if (g_Verbosity >= c_Error) cout << m_XmlTag << ": ERROR in MModuleDepthCalibration: Depth Calibration: got NULL strip hit :( " << endl; return -1; } if (SH->GetEnergy() == 0 ) { - if (g_Verbosity >= c_Error) cout << m_XmlTag << "ERROR in MModuleDepthCalibration: Depth Calibration: got strip without energy :( " << endl; + if (g_Verbosity >= c_Error) cout << m_XmlTag << ": ERROR in MModuleDepthCalibration: Depth Calibration: got strip without energy :( " << endl; return -1; } if (SH->IsLowVoltageStrip()) {