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/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/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/MHit.cxx b/src/MHit.cxx index 24043d6c..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,54 +179,66 @@ void MHit::StreamEvta(ostream& S) { //! Stream the content to an ASCII file - // Assemble the origin information; - vector Origins; + if ((m_GuardRingHit == false) && (m_NoDepth == 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;"< 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 1d24feb8..ab04efe1 100644 --- a/src/MModuleDepthCalibration.cxx +++ b/src/MModuleDepthCalibration.cxx @@ -1,6 +1,4 @@ -/* - * MModuleDepthCalibration.cxx - * + /* * * Copyright (C) 2008-2008 by Andreas Zoglauer, Alex Lowell, * Sean Pike, Carolyn Kierans. @@ -114,62 +112,77 @@ 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. - vector DetList = m_Geometry->GetDetectorList(); - // Look through the Geometry and get the names and thicknesses of all the detectors. + vector DetList = m_Geometry->GetDetectorList(); + 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: "<GetName(); + + // Search first for GeD Strip3D detectors + if (DetectorName.BeginsWith("GeD_") == true) { + DetectorName.RemoveAllInPlace("GeD_"); + DetID = DetectorName.ToInt(); + + // Save the strip pitch and thickness of each detector + if (det->GetTypeName() == "Strip3D") { + 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->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 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; + } + + // 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" )) { + // 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(); } + // Map GR Detector to Detector ID + m_GRDetectors[DetID] = det; } } } + // Check for issues with geometry or file loading, and return appropriate errors if (m_DetectorIDs.size() == 0) { cout<<"No Strip3D detectors were found."<= c_Error) cout<< m_XmlTag << ": No GuardRing Detectors were found."<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); + 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(); - // 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(); + //Initalize position variables: + 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) { + // 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(); + + 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); + H->SetPosition(GlobalPositionGR); + + // Skip past the rest of the calibration and move to next Hit + continue; + + } // If not a GR Hit, perform the depth/position calibration... + + + //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(); + 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(); + 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) { - ++m_Error6; - } - } else { // If the Grade is 0-4, we can handle it. + } + 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); - } + // Take a Hit and separate its activated X- and Y-strips into separate vectors. + vector LVStrips; + vector HVStrips; - double LVEnergyFraction; - double HVEnergyFraction; - MStripHit* LVSH = GetDominantStrip(LVStrips, LVEnergyFraction); - MStripHit* HVSH = GetDominantStrip(HVStrips, HVEnergyFraction); + 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 CTD_s = 0.0; + double LVEnergyFraction; + double HVEnergyFraction; + MStripHit* LVSH = GetDominantStrip(LVStrips, LVEnergyFraction); + MStripHit* HVSH = GetDominantStrip(HVStrips, HVEnergyFraction); - //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; + 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()); + // TODO: Calculate X and Y positions more rigorously using charge sharing. - // 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]; + // 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 + 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 + 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. + Xsigma = m_YPitches[DetID]/sqrt(12.0); + Ysigma = m_XPitches[DetID]/sqrt(12.0); + 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); + // Now try and get z position + double CTD_s = 0.0; - vector CTDVec = GetCTD(DetID, Grade); - vector DepthVec = GetDepth(DetID); + vector* Coeffs = GetPixelCoeffs(PixelCode); + vector CTDVec = GetCTD(DetID, Grade); + vector DepthVec = GetDepth(DetID); - // TODO: For Card Cage, may need to add noise - double LVTiming = LVSH->GetTiming(); - double HVTiming = HVSH->GetTiming(); + 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"); + } else if (DepthVec.size() == 0) { + if (g_Verbosity >= c_Error) cout << m_XmlTag << "Empty Depth 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 ((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); + CTD_s = (CTD - Coeffs->at(1))/(Coeffs->at(0)); //apply inverse stretch and offset + + 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 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; + } - // Confirmed that this matches SP's python code. - CTD_s = (CTD - Coeffs->at(1))/(Coeffs->at(0)); //apply inverse stretch and offset + // 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 + // 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; - double Xmin = * std::min_element(CTDVec.begin(), CTDVec.end()); - double Xmax = * std::max_element(CTDVec.begin(), CTDVec.end()); + for (unsigned int k = 0; k < DepthVec.size(); ++k) { + weighted_depth += prob_dist[k] * DepthVec[k]; + } - double noise = GetTimingNoiseFWHM(PixelCode, H->GetEnergy()); + // Calculate the expectation value of the depth + double mean_depth = weighted_depth/prob_sum; - //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; - } + // Calculate the standard deviation of the depth + double depth_var = 0.0; - // If the CTD is in range, calculate the depth - // Rather than plugging CTD into a spline to get depth, use the depth-CTD relation to calculate a probability-weighted depth value. - // This way we can avoid problems like non-monotonicity or assigning depth to events "outside" the detector - // Note that this requires that we don't massively overestimate the timing noise - else { - // Calculate the probability given timing noise of CTD_s corresponding to the values of depth in DepthVec - // Utlize symmetry of the normal distribution. - vector prob_dist = norm_pdf(CTDVec, CTD_s, noise/2.355); - - // Weight the depth by probability - double prob_sum = 0.0; - for (unsigned int k=0; k < prob_dist.size(); ++k) { - prob_sum += prob_dist[k]; - } - double weighted_depth = 0.0; - - for (unsigned int k = 0; k < DepthVec.size(); ++k) { - weighted_depth += prob_dist[k] * DepthVec[k]; - } - - // Calculate the expectation value of the depth - double mean_depth = weighted_depth/prob_sum; - - // Calculate the standard deviation of the depth - double depth_var = 0.0; - - for (unsigned int k=0; k < DepthVec.size(); ++k) { - depth_var += prob_dist[k] * pow(DepthVec[k] - mean_depth, 2.0); - } - - Zsigma = sqrt(depth_var/prob_sum); - Zpos = mean_depth; - // Zpos = mean_depth - (m_Thicknesses[DetID]/2.0); - - // Add the depth to the GUI histogram. - if (Event->HasStripPairingError()==false) { - if (HasExpos() == true) { - m_ExpoDepthCalibration->AddDepth(DetID, Zpos); - } - } - m_NoError+=1; - } + for (unsigned int k=0; k < DepthVec.size(); ++k) { + depth_var += prob_dist[k] * pow(DepthVec[k] - mean_depth, 2.0); } - 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); + Zsigma = sqrt(depth_var/prob_sum); + Zpos = mean_depth; + // Zpos = mean_depth - (m_Thicknesses[DetID]/2.0); - // 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(); + // Add the depth to the GUI histogram. + if (Event->HasStripPairingError()==false) { + if (HasExpos() == true) { + m_ExpoDepthCalibration->AddDepth(DetID, Zpos); + } + } + + m_NoError+=1; - H->SetPosition(GlobalPosition); - - H->SetPositionResolution(GlobalResolution); + } + + 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); + // 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->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; + } @@ -1100,6 +1143,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/MModuleEventSaver.cxx b/src/MModuleEventSaver.cxx index dc7671c3..c1ff7e08 100644 --- a/src/MModuleEventSaver.cxx +++ b/src/MModuleEventSaver.cxx @@ -75,6 +75,7 @@ MModuleEventSaver::MModuleEventSaver() : MModule() m_InternalFileName = ""; m_Zip = false; m_SaveBadEvents = true; + m_SavePoorQualityEvents = true; m_SaveVetoEvents = true; m_AddTimeTag = false; @@ -345,6 +346,10 @@ bool MModuleEventSaver::AnalyzeEvent(MReadOutAssembly* Event) if (Event->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/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 d4084d83..cbf9cfde 100644 --- a/src/MReadOutAssembly.cxx +++ b/src/MReadOutAssembly.cxx @@ -804,7 +804,24 @@ bool MReadOutAssembly::IsBad() const return false; } - + + +//////////////////////////////////////////////////////////////////////////////// + + +bool MReadOutAssembly::IsPoorQuality() const +{ + //! Returns true if none of the Quality flag has been set + + // Let's not filter out the strips below threshold events since these aren't less quality + //if (m_StripHitBelowThreshold_QualityFlag == true) return true; + if (m_StripPairing_QualityFlag == true) return true; + + if (m_FilteredOut == true) return true; + + return false; +} + //////////////////////////////////////////////////////////////////////////////