From 5ed443ee366d7504ecdd910b58228f9c0ceecd2f Mon Sep 17 00:00:00 2001 From: ckierans Date: Mon, 13 Apr 2026 22:11:49 -0400 Subject: [PATCH 1/4] 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 1d24feb8..0e8b99ac 100644 --- a/src/MModuleDepthCalibration.cxx +++ b/src/MModuleDepthCalibration.cxx @@ -224,202 +224,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 b1f4548d772080f756d15d2409fd264ea3059817 Mon Sep 17 00:00:00 2001 From: ckierans Date: Tue, 14 Apr 2026 01:04:54 -0400 Subject: [PATCH 2/4] Included option to save QA-flagged events --- include/MGUIOptionsEventSaver.h | 3 +++ include/MModuleEventSaver.h | 8 ++++++++ include/MReadOutAssembly.h | 2 ++ src/MGUIOptionsEventSaver.cxx | 5 +++++ src/MModuleDepthCalibration.cxx | 1 - src/MModuleEventSaver.cxx | 10 ++++++++++ src/MReadOutAssembly.cxx | 27 +++++++++++++++++++++++++-- 7 files changed, 53 insertions(+), 3 deletions(-) diff --git a/include/MGUIOptionsEventSaver.h b/include/MGUIOptionsEventSaver.h index fe466302..a6f79050 100644 --- a/include/MGUIOptionsEventSaver.h +++ b/include/MGUIOptionsEventSaver.h @@ -81,6 +81,9 @@ class MGUIOptionsEventSaver : public MGUIOptions //! Checkbutton to save or reject bad events TGCheckButton* m_SaveBadEvents; + //! Checkbutton to save or reject quality flag events + TGCheckButton* m_SavePoorQualityEvents; + //! Checkbutton to save veto events TGCheckButton* m_SaveVetoEvents; diff --git a/include/MModuleEventSaver.h b/include/MModuleEventSaver.h index e3f1d952..44064268 100644 --- a/include/MModuleEventSaver.h +++ b/include/MModuleEventSaver.h @@ -63,6 +63,11 @@ class MModuleEventSaver : public MModule //! Set whether the Bad events should be saved void SetSaveBadEvents(bool SaveBadEvents) { m_SaveBadEvents = SaveBadEvents; } + //! Return true if the Poor Quality events should be saved + bool GetSavePoorQualityEvents() const { return m_SavePoorQualityEvents; } + //! Set whether the Poor Quality events should be saved + void SetSavePoorQualityEvents(bool SavePoorQualityEvents) { m_SavePoorQualityEvents = SavePoorQualityEvents; } + //! Return true if the Veto events should be saved bool GetSaveVetoEvents() const { return m_SaveVetoEvents; } //! Set whether the Veto events should be saved @@ -193,6 +198,9 @@ class MModuleEventSaver : public MModule //! Save bad events bool m_SaveBadEvents; + //! Save poor quality events + bool m_SavePoorQualityEvents; + //! Save Veto events bool m_SaveVetoEvents; diff --git a/include/MReadOutAssembly.h b/include/MReadOutAssembly.h index 9d5005f5..1da4a236 100644 --- a/include/MReadOutAssembly.h +++ b/include/MReadOutAssembly.h @@ -249,6 +249,8 @@ class MReadOutAssembly : public MReadOutSequence bool IsGood() const; //! Returns true if any of the "bad" or "Error" flags has been set bool IsBad() const; + //! Returns true if any of the Quality flags have been set + bool IsPoorQuality() const; //! Set a specific analysis progress void SetAnalysisProgress(uint64_t Progress) { m_AnalysisProgress |= Progress; } diff --git a/src/MGUIOptionsEventSaver.cxx b/src/MGUIOptionsEventSaver.cxx index ab913a43..832203e4 100644 --- a/src/MGUIOptionsEventSaver.cxx +++ b/src/MGUIOptionsEventSaver.cxx @@ -108,6 +108,10 @@ void MGUIOptionsEventSaver::Create() m_SaveBadEvents->SetOn(dynamic_cast(m_Module)->GetSaveBadEvents()); GeneralFrame->AddFrame(m_SaveBadEvents, FirstLabelLayout); + m_SavePoorQualityEvents = new TGCheckButton(GeneralFrame, "Save events with quality flag (QA)", 1); + m_SavePoorQualityEvents->SetOn(dynamic_cast(m_Module)->GetSavePoorQualityEvents()); + GeneralFrame->AddFrame(m_SavePoorQualityEvents, FirstLabelLayout); + m_SaveVetoEvents = new TGCheckButton(GeneralFrame, "Save guard ring and shield veto events (Veto)", 1); m_SaveVetoEvents->SetOn(dynamic_cast(m_Module)->GetSaveVetoEvents()); GeneralFrame->AddFrame(m_SaveVetoEvents, TightButtonLayout); @@ -229,6 +233,7 @@ bool MGUIOptionsEventSaver::OnApply() dynamic_cast(m_Module)->SetFileName(m_FileSelector->GetFileName()); dynamic_cast(m_Module)->SetSaveBadEvents(m_SaveBadEvents->IsOn()); + dynamic_cast(m_Module)->SetSavePoorQualityEvents(m_SavePoorQualityEvents->IsOn()); dynamic_cast(m_Module)->SetSaveVetoEvents(m_SaveVetoEvents->IsOn()); dynamic_cast(m_Module)->SetAddTimeTag(m_AddTimeTag->IsOn()); dynamic_cast(m_Module)->SetSplitFile(m_SplitFile->IsOn()); diff --git a/src/MModuleDepthCalibration.cxx b/src/MModuleDepthCalibration.cxx index 0e8b99ac..5d46150d 100644 --- a/src/MModuleDepthCalibration.cxx +++ b/src/MModuleDepthCalibration.cxx @@ -232,7 +232,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"<IsBad() == true) return true; } + if (m_SavePoorQualityEvents == false) { + if (Event->IsPoorQuality() == true) return true; + } + if (m_SaveVetoEvents == false) { if (Event->IsVeto() == true) return true; } @@ -415,6 +420,10 @@ bool MModuleEventSaver::ReadXmlConfiguration(MXmlNode* Node) if (SaveBadEventsNode != 0) { m_SaveBadEvents = SaveBadEventsNode->GetValueAsBoolean(); } + MXmlNode* SavePoorQualityEventsNode = Node->GetNode("SavePoorQualityEvents"); + if (SavePoorQualityEventsNode != 0) { + m_SavePoorQualityEvents = SavePoorQualityEventsNode->GetValueAsBoolean(); + } MXmlNode* SaveVetoEventsNode = Node->GetNode("SaveVetoEvents"); if (SaveVetoEventsNode != 0) { m_SaveVetoEvents = SaveVetoEventsNode->GetValueAsBoolean(); @@ -480,6 +489,7 @@ MXmlNode* MModuleEventSaver::CreateXmlConfiguration() new MXmlNode(Node, "FileName", m_FileName); new MXmlNode(Node, "Mode", m_Mode); new MXmlNode(Node, "SaveBadEvents", m_SaveBadEvents); + new MXmlNode(Node, "SavePoorQualityEvents", m_SavePoorQualityEvents); new MXmlNode(Node, "SaveVetoEvents", m_SaveVetoEvents); new MXmlNode(Node, "AddTimeTag", m_AddTimeTag); new MXmlNode(Node, "SplitFile", m_SplitFile); diff --git a/src/MReadOutAssembly.cxx b/src/MReadOutAssembly.cxx index d4084d83..905e327f 100644 --- a/src/MReadOutAssembly.cxx +++ b/src/MReadOutAssembly.cxx @@ -618,7 +618,13 @@ void MReadOutAssembly::StreamEvta(ostream& S) } for (unsigned int h = 0; h < m_Hits.size(); ++h) { - m_Hits[h]->StreamEvta(S); + // 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); + } } S<<"CC NStripHits "< Date: Thu, 7 May 2026 21:30:41 -0400 Subject: [PATCH 3/4] GR hits now passed thru revan --- include/MHit.h | 144 +++++++-------- include/MModuleDepthCalibration.h | 1 + src/MHit.cxx | 87 +++++----- src/MModuleDepthCalibration.cxx | 164 +++++++++++------- ...MModuleStripPairingMultiRoundChiSquare.cxx | 3 - src/MReadOutAssembly.cxx | 8 +- 6 files changed, 227 insertions(+), 180 deletions(-) diff --git a/include/MHit.h b/include/MHit.h index 4d79b11f..06fb015b 100644 --- a/include/MHit.h +++ b/include/MHit.h @@ -58,15 +58,15 @@ class MHit //! Return the energy double GetEnergy() const { return m_Energy; } - //! Set the LVenergy - void SetLVEnergy(double LVEnergy) { m_LVEnergy = LVEnergy; } - //! Return the LVenergy - double GetLVEnergy() const { return m_LVEnergy; } + //! Set the LVenergy + void SetLVEnergy(double LVEnergy) { m_LVEnergy = LVEnergy; } + //! Return the LVenergy + double GetLVEnergy() const { return m_LVEnergy; } - //! Set the HVenergy - void SetHVEnergy(double HVEnergy) { m_HVEnergy = HVEnergy; } - //! Return the HVenergy - double GetHVEnergy() const { return m_HVEnergy; } + //! Set the HVenergy + void SetHVEnergy(double HVEnergy) { m_HVEnergy = HVEnergy; } + //! Return the HVenergy + double GetHVEnergy() const { return m_HVEnergy; } //! Set the energy resolution void SetEnergyResolution(double EnergyResolution) { m_EnergyResolution = EnergyResolution; } @@ -89,53 +89,55 @@ class MHit //! Remove a strip hit void RemoveStripHit(MStripHit* StripHit); - //! set cross talk flag - void SetCrossTalkFlag(bool PossibleCrossTalk) {m_PossibleCrossTalk = PossibleCrossTalk;} - //! get cross talk flag value - bool GetCrossTalkFlag() const { return m_PossibleCrossTalk; } + //! set cross talk flag + void SetCrossTalkFlag(bool PossibleCrossTalk) {m_PossibleCrossTalk = PossibleCrossTalk;} + //! get cross talk flag value + bool GetCrossTalkFlag() const { return m_PossibleCrossTalk; } - //! Set guard ring hit flag - void SetGuardRingHitFlag(bool GuardRingHit) {m_GuardRingHit = GuardRingHit;} - //! Get guard ring hit flag - bool GetGuardRingHitFlag() const { return m_GuardRingHit; } - - //! set charge loss flag - void SetChargeLossFlag(bool PossibleChargeLoss) {m_PossibleChargeLoss = PossibleChargeLoss;} - //! get charge loss flag value - bool GetChargeLossFlag() const { return m_PossibleChargeLoss; } - - //! set x strip hit multiple times flag - void SetStripHitMultipleTimesX(bool stripHitMultipleTimesX) {m_StripHitMultipleTimesX = stripHitMultipleTimesX;} - //! get m_StripHitMultipleTimesX - bool GetStripHitMultipleTimesX() const { return m_StripHitMultipleTimesX; } - //! set y strip hit multiple times flag - void SetStripHitMultipleTimesY(bool stripHitMultipleTimesY) {m_StripHitMultipleTimesY = stripHitMultipleTimesY;} - //! get m_StripHitMultipleTimesY - bool GetStripHitMultipleTimesY() const { return m_StripHitMultipleTimesY; } - - //! set charge sharing flag for LV side - void SetChargeSharingLV(bool chargeSharingLV) {m_ChargeSharingLV = chargeSharingLV; } - //! get m_ChargeSharingLV - bool GetChargeSharingLV() const { return m_ChargeSharingLV; } - //! set charge sharing flag for HV side - void SetChargeSharingHV(bool chargeSharingHV) {m_ChargeSharingHV = chargeSharingHV; } - //! get m_ChargeSharingHV - bool GetChargeSharingHV() const { return m_ChargeSharingHV; } - //! Keeping general charge sharing flags because Greedy strip pairing relies on it - void SetChargeSharing(bool chargeSharing) {m_ChargeSharing = chargeSharing; } - //! get m_ChargeSharing - bool GetChargeSharing() const { return m_ChargeSharing; } - //! set m_NoDepth - void SetNoDepth(bool X = true) { m_NoDepth = X;} - //! get m_NoDepth - bool GetNoDepth(void) const { return m_NoDepth; } - //! set m_IsNonDominantNeighborStrip - void SetIsNondominantNeighborStrip(bool X = true) {m_IsNonDominantNeighborStrip = X;} - //! get m_IsNonDominantNeighborStrip - bool GetIsNondominantNeighborStrip(void) const {return m_IsNonDominantNeighborStrip;} + //! Set guard ring hit flag + void SetGuardRingHitFlag(bool GuardRingHit) {m_GuardRingHit = GuardRingHit;} + //! Get guard ring hit flag + bool GetGuardRingHitFlag() const { return m_GuardRingHit; } + + //! set charge loss flag + void SetChargeLossFlag(bool PossibleChargeLoss) {m_PossibleChargeLoss = PossibleChargeLoss;} + //! get charge loss flag value + bool GetChargeLossFlag() const { return m_PossibleChargeLoss; } + + //! set x strip hit multiple times flag + void SetStripHitMultipleTimesX(bool stripHitMultipleTimesX) {m_StripHitMultipleTimesX = stripHitMultipleTimesX;} + //! get m_StripHitMultipleTimesX + bool GetStripHitMultipleTimesX() const { return m_StripHitMultipleTimesX; } + //! set y strip hit multiple times flag + void SetStripHitMultipleTimesY(bool stripHitMultipleTimesY) {m_StripHitMultipleTimesY = stripHitMultipleTimesY;} + //! get m_StripHitMultipleTimesY + bool GetStripHitMultipleTimesY() const { return m_StripHitMultipleTimesY; } + + //! set charge sharing flag for LV side + void SetChargeSharingLV(bool chargeSharingLV) {m_ChargeSharingLV = chargeSharingLV; } + //! get m_ChargeSharingLV + bool GetChargeSharingLV() const { return m_ChargeSharingLV; } + //! set charge sharing flag for HV side + void SetChargeSharingHV(bool chargeSharingHV) {m_ChargeSharingHV = chargeSharingHV; } + //! get m_ChargeSharingHV + bool GetChargeSharingHV() const { return m_ChargeSharingHV; } + //! Keeping general charge sharing flags because Greedy strip pairing relies on it + void SetChargeSharing(bool chargeSharing) {m_ChargeSharing = chargeSharing; } + //! get m_ChargeSharing + bool GetChargeSharing() const { return m_ChargeSharing; } + + //! set m_NoDepth + void SetNoDepth(bool X = true) { m_NoDepth = X;} + //! get m_NoDepth + bool GetNoDepth(void) const { return m_NoDepth; } + + //! set m_IsNonDominantNeighborStrip + void SetIsNondominantNeighborStrip(bool X = true) {m_IsNonDominantNeighborStrip = X;} + //! get m_IsNonDominantNeighborStrip + bool GetIsNondominantNeighborStrip(void) const {return m_IsNonDominantNeighborStrip;} - //! Set the origins from the simulations (take care of duplicates) - void AddOrigins(vector Origins); + //! Set the origins from the simulations (take care of duplicates) + void AddOrigins(vector Origins); //! Get the origins from the simulation vector GetOrigins() const { return m_Origins; } @@ -186,28 +188,28 @@ class MHit //! List of strip hits vector m_StripHits; - //! Flag: possible cross talk - bool m_PossibleCrossTalk; - //! Flag: possible charge loss - bool m_PossibleChargeLoss; - //! Flag: hit containing guard ring strip - bool m_GuardRingHit; + //! Flag: possible cross talk + bool m_PossibleCrossTalk; + //! Flag: possible charge loss + bool m_PossibleChargeLoss; + //! Flag: hit containing guard ring strip + bool m_GuardRingHit; - //! true if hit contains strip that was hit multiple times on X - bool m_StripHitMultipleTimesX = false; - //! true if hit contains strip that was hit multiple times on Y - bool m_StripHitMultipleTimesY = false; + //! true if hit contains strip that was hit multiple times on X + bool m_StripHitMultipleTimesX = false; + //! true if hit contains strip that was hit multiple times on Y + bool m_StripHitMultipleTimesY = false; - //! true if hit contains charge sharing - bool m_ChargeSharing; - bool m_ChargeSharingLV; - bool m_ChargeSharingHV; + //! true if hit contains charge sharing + bool m_ChargeSharing; + bool m_ChargeSharingLV; + bool m_ChargeSharingHV; - //! true if depth is invalid, either because the pixel depth was uncalibrated, the hit was mapped too far out of the detector,or there was no timing data - bool m_NoDepth; + //! true if depth is invalid, either because the pixel depth was uncalibrated, the hit was mapped too far out of the detector,or there was no timing data + bool m_NoDepth; - //! true if hit was made from a charge sharing event using a neighbor strip that had the lowere energy fraction - bool m_IsNonDominantNeighborStrip; + //! true if hit was made from a charge sharing event using a neighbor strip that had the lowere energy fraction + bool m_IsNonDominantNeighborStrip; //! Origin IAs from simulations vector m_Origins; diff --git a/include/MModuleDepthCalibration.h b/include/MModuleDepthCalibration.h index 8489a8ac..6fd3169d 100644 --- a/include/MModuleDepthCalibration.h +++ b/include/MModuleDepthCalibration.h @@ -169,6 +169,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 24043d6c..ded5c7ac 100644 --- a/src/MHit.cxx +++ b/src/MHit.cxx @@ -165,54 +165,63 @@ void MHit::StreamEvta(ostream& S) { //! Stream the content to an ASCII file - // Assemble the origin information; - vector Origins; + if (m_GuardRingHit == false) { + // Assemble the origin information; + vector Origins; - // Fix the origins: only those existing both on x and y strips count - vector xOrigins; - vector yOrigins; - for (unsigned int s = 0; s < GetNStripHits(); ++s) { - //GetStripHit(s)->StreamRoa(cout); - vector NewOrigins = GetStripHit(s)->GetOrigins(); - if (GetStripHit(s)->IsLowVoltageStrip() == true) { - for (int o: NewOrigins) { - xOrigins.push_back(o); - } - } else { - for (int o: NewOrigins) { - yOrigins.push_back(o); + // Fix the origins: only those existing both on x and y strips count + vector xOrigins; + vector yOrigins; + for (unsigned int s = 0; s < GetNStripHits(); ++s) { + //GetStripHit(s)->StreamRoa(cout); + vector NewOrigins = GetStripHit(s)->GetOrigins(); + if (GetStripHit(s)->IsLowVoltageStrip() == true) { + for (int o: NewOrigins) { + xOrigins.push_back(o); + } + } else { + for (int o: NewOrigins) { + yOrigins.push_back(o); + } } } - } - sort(xOrigins.begin(), xOrigins.end()); - xOrigins.erase(unique(xOrigins.begin(), xOrigins.end()), xOrigins.end()); - sort(yOrigins.begin(), yOrigins.end()); - yOrigins.erase(unique(yOrigins.begin(), yOrigins.end()), yOrigins.end()); + sort(xOrigins.begin(), xOrigins.end()); + xOrigins.erase(unique(xOrigins.begin(), xOrigins.end()), xOrigins.end()); + sort(yOrigins.begin(), yOrigins.end()); + yOrigins.erase(unique(yOrigins.begin(), yOrigins.end()), yOrigins.end()); - set_intersection(xOrigins.begin(), xOrigins.end(), - yOrigins.begin(), yOrigins.end(), - std::back_inserter(Origins)); - - if ((xOrigins.size() != 0 || yOrigins.size() != 0) && Origins.size() == 0) { - // This is the case when the strip pairing got screwed up completely, and the hits got mixed - // In this case we need to keep the mixing: - for (unsigned int s = 0; s < GetNStripHits(); ++s) { - vector NewOrigins = GetStripHit(s)->GetOrigins(); - for (int o: NewOrigins) { - Origins.push_back(o); + set_intersection(xOrigins.begin(), xOrigins.end(), + yOrigins.begin(), yOrigins.end(), + std::back_inserter(Origins)); + + if ((xOrigins.size() != 0 || yOrigins.size() != 0) && Origins.size() == 0) { + // This is the case when the strip pairing got screwed up completely, and the hits got mixed + // In this case we need to keep the mixing: + for (unsigned int s = 0; s < GetNStripHits(); ++s) { + vector NewOrigins = GetStripHit(s)->GetOrigins(); + for (int o: NewOrigins) { + Origins.push_back(o); + } } + sort(Origins.begin(), Origins.end()); + Origins.erase(unique(Origins.begin(), Origins.end()), Origins.end()); } - sort(Origins.begin(), Origins.end()); - Origins.erase(unique(Origins.begin(), Origins.end()), Origins.end()); - } - S<<"HT 3;"< 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 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 == true) { - DetID = 11; - } - + + unsigned int DetID; 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(); - - if (g_Verbosity >= c_Info) { - cout << "Found detector " << det_name << " corresponding to DetID=" << 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 { - if (g_Verbosity >= c_Error) { - cout<<"ERROR in MModuleDepthCalibration::Initialize: Found a duplicate detector: "<GetSensitiveVolume(0)->GetName()<<", type: "<GetTypeName()<GetName(); + + // Search first for GeD Strip3D detectors + if (DetectorName.BeginsWith("GeD_") == true) { + DetectorName.RemoveAllInPlace("GeD_"); + cout<<"Found GeD name"<GetTypeName() == "Strip3D") { + cout<<"DetNames LOOP"<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(); + + if (g_Verbosity >= c_Info) { + cout << "Found detector " << det->GetName() << " corresponding to DetID = " << 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; } + + // Add to list of of Detectors + m_DetectorIDs.push_back(DetID); + // Map Detector to Detector ID + m_Detectors[DetID] = det; + } else { - if (g_Verbosity >= c_Error) { - cout<<"ERROR in MModuleDepthCalibration::Initialize: Found a Strip3D detector with "<GetNSensitiveVolumes()<<" Sensitive Volumes."<= c_Error) cout<< m_XmlTag <<"GeD Detector Volume " << det->GetName() << " is not Strip3D! " << endl; + return false; + } + + } else if (DetectorName.BeginsWith("GuardRingDetector") == true) { + + if ((det->GetTypeName() == "Simple") || (det->GetTypeName() == "Scintillator" )) { + if (DetectorName == "GuardRingDetector") { + DetID = 0; + } else { + DetectorName.RemoveAllInPlace("GuardRingDetector_GeD_"); + DetID = DetectorName.ToInt(); } + cout<<"Det ID: "<= c_Error) cout<< m_XmlTag << ": No GuardRing Detectors were found."<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(); @@ -259,7 +306,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; @@ -274,20 +321,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 @@ -300,15 +346,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); @@ -337,8 +383,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()); @@ -408,7 +452,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); @@ -1093,6 +1136,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 56681015..82eb5256 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 905e327f..cbf9cfde 100644 --- a/src/MReadOutAssembly.cxx +++ b/src/MReadOutAssembly.cxx @@ -618,13 +618,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/4] Include GR and XE hits in evta --- src/MHit.cxx | 98 +++++++++++++-------- src/MModuleDepthCalibration.cxx | 147 +++++++++++++++++--------------- 2 files changed, 141 insertions(+), 104 deletions(-) diff --git a/src/MHit.cxx b/src/MHit.cxx index ded5c7ac..2e465f12 100644 --- a/src/MHit.cxx +++ b/src/MHit.cxx @@ -82,6 +82,20 @@ void MHit::Clear() m_StripHits.clear(); m_Origins.clear(); + + m_HitQuality = g_DoubleNotDefined; + + m_PossibleCrossTalk = false; + m_PossibleChargeLoss = false; + m_GuardRingHit = false; + m_StripHitMultipleTimesX = false; + m_StripHitMultipleTimesY = false; + m_ChargeSharing = false; + m_ChargeSharingLV = false; + m_ChargeSharingHV = false; + m_NoDepth = false; + m_IsNonDominantNeighborStrip = false; + } @@ -165,7 +179,7 @@ void MHit::StreamEvta(ostream& S) { //! Stream the content to an ASCII file - if (m_GuardRingHit == false) { + if ((m_GuardRingHit == false) && (m_NoDepth == false)) { // Assemble the origin information; vector Origins; @@ -216,13 +230,16 @@ 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; - } - */ + /* + 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; + } + */ } diff --git a/src/MModuleDepthCalibration.cxx b/src/MModuleDepthCalibration.cxx index 10d58894..ab04efe1 100644 --- a/src/MModuleDepthCalibration.cxx +++ b/src/MModuleDepthCalibration.cxx @@ -112,29 +112,23 @@ MModuleDepthCalibration::~MModuleDepthCalibration() bool MModuleDepthCalibration::Initialize() { - // The detectors need to be in the same order as DetIDs. - // ie DetID=0 should be the 0th detector in m_Detectors, DetID=1 should the 1st, etc. + // Look through the Geometry and get the names and thicknesses of all the detectors. 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 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; MDDetector* det = DetList[i]; - cout<<"Found detector in DetList: "<GetSensitiveVolume(0)->GetName()<<", type: "<GetTypeName()<GetName(); // Search first for GeD Strip3D detectors if (DetectorName.BeginsWith("GeD_") == true) { DetectorName.RemoveAllInPlace("GeD_"); - cout<<"Found GeD name"<GetTypeName() == "Strip3D") { - cout<<"DetNames LOOP"<GetStructuralSize().GetZ()); MDStrip3D* strip = dynamic_cast(det); m_XPitches[DetID] = strip->GetPitchX(); @@ -151,7 +145,7 @@ bool MModuleDepthCalibration::Initialize() cout << "Y strip pitch: " << m_YPitches[DetID] << endl; } - // Add to list of of Detectors + // Add to list of Detectors m_DetectorIDs.push_back(DetID); // Map Detector to Detector ID m_Detectors[DetID] = det; @@ -161,21 +155,25 @@ bool MModuleDepthCalibration::Initialize() return false; } + // Search for Guard Rings } else if (DetectorName.BeginsWith("GuardRingDetector") == true) { + // Need to check for Simple and Scintillator detector types for different COSI mass models if ((det->GetTypeName() == "Simple") || (det->GetTypeName() == "Scintillator" )) { - if (DetectorName == "GuardRingDetector") { + // For single detector mass models, the GR detector had no number, set it to 0 + if (DetectorName == "GuardRingDetector") { DetID = 0; } else { DetectorName.RemoveAllInPlace("GuardRingDetector_GeD_"); DetID = DetectorName.ToInt(); } - cout<<"Det ID: "<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 @@ -270,19 +264,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. @@ -298,7 +288,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) { @@ -321,7 +314,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; @@ -362,7 +354,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"); @@ -385,64 +377,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); @@ -454,14 +455,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; + }