diff --git a/include/MReadOutAssembly.h b/include/MReadOutAssembly.h index a0aa8da..77ac046 100644 --- a/include/MReadOutAssembly.h +++ b/include/MReadOutAssembly.h @@ -55,14 +55,8 @@ class MReadOutAssembly : public MReadOutSequence //! Delete Hits void DeleteHits(); - - /* - //! Set the ID of this event - void SetID(unsigned long ID) { m_ID = ID; } - //! Return the ID of this event - unsigned long GetID() const { return m_ID; } - */ - + + //! TODO Scrub all clock/time variables for COSI SMEX //! Set the Frame Counter of this event void SetFC(unsigned int FC) { m_FC = FC; } //! Return the Frame Counter of this event @@ -76,12 +70,6 @@ class MReadOutAssembly : public MReadOutSequence void SetCL(uint64_t CL) { m_CL = CL;} uint64_t GetCL() const { return m_CL;} - /* - //! Set and get the Time of this event - void SetTime(MTime Time) { m_Time = Time; } - MTime GetTime() const { return m_Time; } - */ - //! Set and get the Modified Julian Date of this event void SetMJD(double MJD) { m_MJD = MJD; } double GetMJD() const { return m_MJD; } @@ -90,24 +78,25 @@ class MReadOutAssembly : public MReadOutSequence void SetTimeUTC(const MTime& TimeUTC) { m_EventTimeUTC = TimeUTC; } MTime GetTimeUTC() const { return m_EventTimeUTC; } + //! TODO Scrub Aspect in nuclearizer for COSI SMEX //! Set the aspect void SetAspect(MAspect* Aspect) { if (m_Aspect != 0) delete m_Aspect; m_Aspect = Aspect; } //! Get the aspect - will be zero if the aspect has not been set! MAspect* GetAspect() { return m_Aspect; } - //! Set and get simulation aspect information - void SetGalacticPointingXAxisTheta(double theta){ m_GalacticPointingXAxisTheta = theta; } - void SetGalacticPointingXAxisPhi(double phi){ m_GalacticPointingXAxisPhi = phi; } - void SetGalacticPointingZAxisTheta(double theta){ m_GalacticPointingZAxisTheta = theta; } - void SetGalacticPointingZAxisPhi(double phi){ m_GalacticPointingZAxisPhi = phi; } + //! Set and get simulation aspect information + void SetGalacticPointingXAxisTheta(double theta){ m_GalacticPointingXAxisTheta = theta; } + void SetGalacticPointingXAxisPhi(double phi){ m_GalacticPointingXAxisPhi = phi; } + void SetGalacticPointingZAxisTheta(double theta){ m_GalacticPointingZAxisTheta = theta; } + void SetGalacticPointingZAxisPhi(double phi){ m_GalacticPointingZAxisPhi = phi; } - double GetGalacticPointingXAxisTheta(){ if (m_HasSimAspectInfo){return m_GalacticPointingXAxisTheta;} else{return 0;}} - double GetGalacticPointingXAxisPhi(){ if (m_HasSimAspectInfo){return m_GalacticPointingXAxisPhi;} else{return 0;}} - double GetGalacticPointingZAxisTheta(){ if (m_HasSimAspectInfo){return m_GalacticPointingZAxisTheta;} else{return 0;}} - double GetGalacticPointingZAxisPhi(){ if (m_HasSimAspectInfo){return m_GalacticPointingZAxisPhi;} else{return 0;}} + double GetGalacticPointingXAxisTheta(){ if (m_HasSimAspectInfo){return m_GalacticPointingXAxisTheta;} else{return 0;}} + double GetGalacticPointingXAxisPhi(){ if (m_HasSimAspectInfo){return m_GalacticPointingXAxisPhi;} else{return 0;}} + double GetGalacticPointingZAxisTheta(){ if (m_HasSimAspectInfo){return m_GalacticPointingZAxisTheta;} else{return 0;}} + double GetGalacticPointingZAxisPhi(){ if (m_HasSimAspectInfo){return m_GalacticPointingZAxisPhi;} else{return 0;}} - void SetSimAspectInfo(bool TF){ m_HasSimAspectInfo = TF; } - bool HasSimAspectInfo(){ return m_HasSimAspectInfo; } + void SetSimAspectInfo(bool TF){ m_HasSimAspectInfo = TF; } + bool HasSimAspectInfo(){ return m_HasSimAspectInfo; } //! Find out if the event contains strip hits in a given detector @@ -128,11 +117,6 @@ class MReadOutAssembly : public MReadOutSequence //! Return the trigger flag bool GetTrigger() const { return m_Trigger; } - //! Set the aspect good flag - void SetAspectGood(bool AspectGood = true) { m_AspectGood = AspectGood; } - //! Return the aspect good flag - bool GetAspectGood() const { return m_AspectGood; } - //! Return the number of strip hits unsigned int GetNStripHits() const { return m_StripHits.size(); } //! Return strip hit i @@ -143,6 +127,7 @@ class MReadOutAssembly : public MReadOutSequence void RemoveStripHit(unsigned int i); //! Return the number of T Only strip hits + //! TODO Is this a hold-over from balloon days? unsigned int GetNStripHitsTOnly() const { return m_StripHitsTOnly.size(); } //! Return strip hit i MStripHit* GetStripHitTOnly(unsigned int i); @@ -160,7 +145,6 @@ class MReadOutAssembly : public MReadOutSequence //! Remove a crystal hit void RemoveCrystalHit(unsigned int i); - //! Return the number of guardring hits unsigned int GetNGuardringHits() const { return m_GuardringHits.size(); } //! Return guardring hit i @@ -179,13 +163,13 @@ class MReadOutAssembly : public MReadOutSequence //! Return the number of simulation hits // TODO: Remove - part of m_SimEvent - unsigned int GetNHitsSim() const { return m_HitsSim.size(); } +// unsigned int GetNHitsSim() const { return m_HitsSim.size(); } //! Return simulation hit i // TODO: Remove - part of m_SimEvent - MHit* GetHitSim(unsigned int i); +// MHit* GetHitSim(unsigned int i); //! Move hits to simulation hits list // TODO: Why ?? - void MoveHitsToSim() {m_HitsSim = m_Hits; m_Hits.clear();} +// void MoveHitsToSim() {m_HitsSim = m_Hits; m_Hits.clear();} /* //! Return the number of simulation interactions @@ -226,86 +210,68 @@ class MReadOutAssembly : public MReadOutSequence list& GetDEECrystalHitListReference() { return m_DEECrystalHits; } - //! Return the number of read outs - //unsigned int GetNReadOuts() const { return m_ReadOuts.size(); } - //! Return read out i - throws an exception of the index is not found - //MReadOut& GetReadOut(unsigned int i); - //! Add a read out - // void AddReadOut(MReadOut& ReadOut) {} - //! Remove a read out - does do nothing if the index is not found - //void RemoveReadOut(unsigned int i); + //Track BD Flags + //! Set the energy-calibration-Error flag + void SetEnergyCalibrationError(bool Flag = true, MString Text = "") { m_EnergyCalibrationError = Flag; m_EnergyCalibrationErrorString = Text; } + //! Get the energy-calibration-Error flag + bool IsEnergyCalibrationError() const { return m_EnergyCalibrationError; } + + //! Set the strip-pairing-Error flag + void SetStripPairingError(bool Flag = true, MString Text = "") { m_StripPairingError = Flag; m_StripPairingErrorString = Text; } + //! Get the strip-pairing-Error flag + bool IsStripPairingError() const { return m_StripPairingError; } + + //! Set the depth-calibration-Error flag + void SetDepthCalibrationError(bool Flag = true, MString Text = "") { m_DepthCalibrationError = Flag; m_DepthCalibrationErrorString = Text; } + //! Get the depth-calibration-Error flag + bool IsDepthCalibrationError() const { return m_DepthCalibrationError; } + + //! Set the event-reconstruction-Error flag + void SetEventReconstructionError(bool Flag = true, MString Text = "") { m_EventReconstructionError = Flag; m_EventReconstructionErrorString = Text; } + //! Get the event-reconstruction-Error flag + bool IsEventReconstructionError() const { return m_EventReconstructionError; } + + // Track Quality Flags + + //! Set the Strip Hit Below Threshold quality flag + void SetStripHitBelowThreshold_QualityFlag(bool Flag = true, MString Text = ""){ m_StripHitBelowThreshold_QualityFlag = Flag; m_StripHitBelowThresholdString_QualityFlag = Text;} + //! Get the Strip Hit Below Threshold quality flag + bool IsStripHitBelowThreshold_QualityFlag() const { return m_StripHitBelowThreshold_QualityFlag; } + + //! Track the energy and number of strip hits removed through the threshold cut + void AddStripHitBelowThreshold(double Energy) { m_StripHitBelowThreshold_Energy = m_StripHitBelowThreshold_Energy + Energy; m_StripHitBelowThreshold_Number = m_StripHitBelowThreshold_Number + 1; } + //! Get total energy removed through threshold cut + double GetStripHitBelowThreshold_Energy() const { return m_StripHitBelowThreshold_Energy; } + //! Get total number of strip hits removed through threshold cut + int GetStripHitBelowThreshold_Number() const { return m_StripHitBelowThreshold_Number;} + + //! Set the Reduced Chi^2 used in MultiRoundChiSquare module + void SetReducedChiSquare(double ReducedChiSquare) { m_ReducedChiSquare = ReducedChiSquare; } + //! Return the Reduced Chi^2 + double GetReducedChiSquare() const { return m_ReducedChiSquare; } + //! Set the Quality of this Event used in Greedy Strip pairing module + //! TODO Change name of this variable to be more descriptive + void SetEventQuality(double EventQuality){ m_EventQuality = EventQuality; } + //!Return the Quality of this Event + double GetEventQuality() const { return m_EventQuality; } - //! Set the aspect-incomplete flag - void SetAspectIncomplete(bool Flag = true, MString Text = "") { m_AspectIncomplete = Flag; m_AspectIncompleteString = Text; } - //! Get the aspect-incomplete flag - bool IsAspectIncomplete() const { return m_AspectIncomplete; } - - //! Set the time-incomplete flag - void SetTimeIncomplete(bool Flag = true, MString Text = "") { m_TimeIncomplete = Flag; m_TimeIncompleteString = Text; } - //! Get the time-incomplete flag - bool IsTimeIncomplete() const { return m_TimeIncomplete; } - - //! Set the energy-calibration-incomplete flag for strips without a calibration - void SetEnergyCalibrationIncomplete_BadStrip(bool Flag = true, MString Text = "") { m_EnergyCalibrationIncomplete_BadStrip = Flag; m_EnergyCalibrationIncomplete_BadStripString = Text; } - //! Get the energy-calibration-incomplete flag for strips without a calibration - bool IsEnergyCalibrationIncomplete_BadStrip() const { return m_EnergyCalibrationIncomplete_BadStrip; } - - //! Set the energy-calibration-incomplete flag - void SetEnergyCalibrationIncomplete(bool Flag = true, MString Text = "") { m_EnergyCalibrationIncomplete = Flag; m_EnergyCalibrationIncompleteString = Text; } - //! Get the energy-calibration-incomplete flag - bool IsEnergyCalibrationIncomplete() const { return m_EnergyCalibrationIncomplete; } - - //! Set the energy resolution calibration incomplete flag - void SetEnergyResolutionCalibrationIncomplete(bool Flag = true, MString Text = "") { m_EnergyResolutionCalibrationIncomplete = Flag; m_EnergyResolutionCalibrationIncompleteString = Text;} - //! Get the energy resolution calibration incomplete flag - bool IsEnergyResolutionCalibrationIncomplete() const { return m_EnergyResolutionCalibrationIncomplete; } - - //! set the Strip Hit Below Threshold flag - void SetStripHitBelowThreshold(bool Flag = true, MString Text = ""){ - m_StripHitBelowThreshold = Flag; m_StripHitBelowThreshold = Text;} - //! get the Strip Hit Below Threshold flag - bool IsStripHitBelowThreshold() const { return m_StripHitBelowThreshold; } - - - //! Set the strip-pairing-incomplete flag - void SetStripPairingIncomplete(bool Flag = true, MString Text = "") { m_StripPairingIncomplete = Flag; m_StripPairingIncompleteString = Text; } - //! Get the strip-pairing-incomplete flag - bool IsStripPairingIncomplete() const { return m_StripPairingIncomplete; } - //! Set the LLD Event flag - void SetLLDEvent(bool Flag = true, MString Text = "") { m_LLDEvent = Flag; m_LLDEventString = Text; } - //! Get the LLD Event flag - bool IsLLDEvent() const {return m_LLDEvent;} + // Track Vetos - //! Set the depth-calibration-incomplete flag - void SetDepthCalibrationIncomplete(bool Flag = true, MString Text = "") { m_DepthCalibrationIncomplete = Flag; m_DepthCalibrationIncompleteString = Text; } - //! Get the depth-calibration-incomplete flag - bool IsDepthCalibrationIncomplete() const { return m_DepthCalibrationIncomplete; } + //! Returns true if any of the "veto" flags have been set + bool IsVeto() const; - //! Set the depth-calibration out of range flag - void SetDepthCalibration_OutofRange(bool Flag = true, MString Text = "") {m_DepthCalibration_OutofRange = Flag; m_DepthCalibration_OutofRangeString = Text; } - //! Get the depth calibration out of range flag - bool IsDepthCalibration_OutofRange() const { return m_DepthCalibration_OutofRange; } //! Set the filtered-out flag void SetFilteredOut(bool Flag = true) { m_FilteredOut = Flag; } //! Get the filgtered-out flag bool IsFilteredOut() const { return m_FilteredOut; } - //! Set the event-reconstruction-incomplete flag - void SetEventReconstructionIncomplete(bool Flag = true, MString Text = "") { m_EventReconstructionIncomplete = Flag; m_EventReconstructionIncompleteString = Text; } - //! Get the event-reconstruction-incomplete flag - bool IsEventReconstructionIncomplete() const { return m_EventReconstructionIncomplete; } - - - //! Returns true if any of the "veto" flags have been set - bool IsVeto() const; - - //! Returns true if none of the "bad" or "incomplete" flags has been set and the event has not been filtered out or rejected + //! Returns true if none of the "bad" or "Error" flags has been set and the event has not been filtered out or rejected bool IsGood() const; - //! Returns true if any of the "bad" or "incomplete" flags has been set + //! Returns true if any of the "bad" or "Error" flags has been set bool IsBad() const; //! Set a specific analysis progress @@ -315,11 +281,6 @@ class MReadOutAssembly : public MReadOutSequence //! Return the analysis progress flag uint64_t GetAnalysisProgress() const { return m_AnalysisProgress; } - //! Set the Quality of this Event - void SetEventQuality(double EventQuality){ m_EventQuality = EventQuality; } - //!Return the Quality of this Event - double GetEventQuality() const { return m_EventQuality; } - //! Parse some content from a line bool Parse(MString& Line, int Version = 1); @@ -345,11 +306,6 @@ class MReadOutAssembly : public MReadOutSequence //! Get the MTime corresponding to absolute UTC time MTime GetAbsoluteTime() const {return m_EventTimeUTC; } - //! Set the Reduced Chi^2 - void SetReducedChiSquare(double ReducedChiSquare) { m_ReducedChiSquare = ReducedChiSquare; } - //! Return the Reduced Chi^2 - double GetReducedChiSquare() const { return m_ReducedChiSquare; } - // protected methods: protected: //MReadOutAssembly() {}; @@ -366,8 +322,6 @@ class MReadOutAssembly : public MReadOutSequence // private members: private: - //! ID of this event - // unsigned long m_ID; // in base class //! Frame Counter of this event unsigned int m_FC; @@ -386,17 +340,14 @@ class MReadOutAssembly : public MReadOutSequence //! The aspect information - will be zero if not set! MAspect* m_Aspect; - //Added by Clio: - //! The aspect information from the simulation, only used in DEE - // (Simulation aspect information doesn't have everything in Aspect packet) - double m_GalacticPointingXAxisTheta; - double m_GalacticPointingXAxisPhi; - double m_GalacticPointingZAxisTheta; - double m_GalacticPointingZAxisPhi; - bool m_HasSimAspectInfo; - - //! Quality of this event - double m_EventQuality; + //Added by Clio: + //! The aspect information from the simulation, only used in DEE + // (Simulation aspect information doesn't have everything in Aspect packet) + double m_GalacticPointingXAxisTheta; + double m_GalacticPointingXAxisPhi; + double m_GalacticPointingZAxisTheta; + double m_GalacticPointingZAxisPhi; + bool m_HasSimAspectInfo; //! Guard ring veto flag bool m_GuardRingVeto; @@ -411,7 +362,7 @@ class MReadOutAssembly : public MReadOutSequence bool m_AspectGood; //! Whether event contains strip hits in given detector - bool m_InDetector[12]; + bool m_InDetector[16]; //! List of strip hits vector m_StripHits; @@ -445,40 +396,31 @@ class MReadOutAssembly : public MReadOutSequence //! The physical event from event reconstruction MPhysicalEvent* m_PhysicalEvent; - //! Reduced Chi^2 of the Strip Paired Event - double m_ReducedChiSquare; - - // Flags indicating the quality of the event - bool m_AspectIncomplete; - MString m_AspectIncompleteString; - - bool m_TimeIncomplete; - MString m_TimeIncompleteString; - - bool m_EnergyCalibrationIncomplete_BadStrip; - MString m_EnergyCalibrationIncomplete_BadStripString; - bool m_EnergyCalibrationIncomplete; - MString m_EnergyCalibrationIncompleteString; - - bool m_EnergyResolutionCalibrationIncomplete; - MString m_EnergyResolutionCalibrationIncompleteString; + // Flags indicating bad events: + bool m_EnergyCalibrationError; + MString m_EnergyCalibrationErrorString; - bool m_StripPairingIncomplete; - MString m_StripPairingIncompleteString; + bool m_StripPairingError; + MString m_StripPairingErrorString; - bool m_LLDEvent; - MString m_LLDEventString; - - bool m_DepthCalibrationIncomplete; - MString m_DepthCalibrationIncompleteString; - bool m_DepthCalibration_OutofRange; - MString m_DepthCalibration_OutofRangeString; + bool m_DepthCalibrationError; + MString m_DepthCalibrationErrorString; + + bool m_EventReconstructionError; + MString m_EventReconstructionErrorString; + + // Flags indicating the quality of the event: quality warning, but not to be filtered out + bool m_StripHitBelowThreshold_QualityFlag; + MString m_StripHitBelowThresholdString_QualityFlag; + double m_StripHitBelowThreshold_Energy; + int m_StripHitBelowThreshold_Number; - bool m_StripHitBelowThreshold; - MString m_StripHitBelowThresholdString; + //! Reduced Chi^2 of the Strip Paired Event + double m_ReducedChiSquare; - bool m_EventReconstructionIncomplete; - MString m_EventReconstructionIncompleteString; + //! Quality of this event in Greedy strip pairing + //! TODO change variable name... + double m_EventQuality; //! True if event has been filtered out bool m_FilteredOut; diff --git a/src/MModuleDepthCalibration.cxx b/src/MModuleDepthCalibration.cxx index 447f739..f3a7d7d 100644 --- a/src/MModuleDepthCalibration.cxx +++ b/src/MModuleDepthCalibration.cxx @@ -301,7 +301,7 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) } */ - //Event->SetDepthCalibrationIncomplete(); //AWL1x1 + //Event->SetDepthCalibrationError(); //AWL1x1 } else if( (XStrips.size() == 2) && (YStrips.size() == 2) ){ //in this case use depth from dominant strips but use weighted X and Y positions @@ -328,7 +328,7 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) GlobalPosition = m_DetectorVolumes[DetID]->GetPositionInWorldVolume(LocalPosition); H->SetPosition(GlobalPosition); H->SetPositionResolution(PositionResolution); - //Event->SetDepthCalibrationIncomplete(); //AWL1x1 + //Event->SetDepthCalibrationError(); //AWL1x1 } else { //set too many SH bad flag @@ -339,24 +339,24 @@ bool MModuleDepthCalibration::AnalyzeEvent(MReadOutAssembly* Event) //good ++m_NoError; } else if( Error == 1 ){ - Event->SetDepthCalibrationIncomplete(); + Event->SetDepthCalibrationError(); ++m_Error1; } else if( Error == 2 ){ - Event->SetDepthCalibrationIncomplete(); + Event->SetDepthCalibrationError(); ++m_Error2; } else if( Error == 3){ //Hits that were missing timing information //EHist->Fill(H->GetEnergy()); //don't set the globally bad flag - //Event->SetDepthCalibrationIncomplete(); - //Event->SetDepthCalibrationIncomplete(); //AWL1x1 + //Event->SetDepthCalibrationError(); + //Event->SetDepthCalibrationError(); //AWL1x1 ++m_Error3; } else if( Error == 4){ //hit was bad because of StripHitMultipleTimes flag from strip pairing ++m_Error4; - Event->SetDepthCalibrationIncomplete(); + Event->SetDepthCalibrationError(); } else if( Error == -1){ - Event->SetDepthCalibrationIncomplete(); + Event->SetDepthCalibrationError(); ++m_ErrorSH; } diff --git a/src/MModuleDepthCalibration2024.cxx b/src/MModuleDepthCalibration2024.cxx index 46de19a..b586cad 100644 --- a/src/MModuleDepthCalibration2024.cxx +++ b/src/MModuleDepthCalibration2024.cxx @@ -140,19 +140,26 @@ bool MModuleDepthCalibration2024::Initialize() m_YPitches[DetID] = strip->GetPitchY(); m_NXStrips[DetID] = strip->GetNStripsX(); m_NYStrips[DetID] = strip->GetNStripsY(); - 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; + + 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 { - cout<<"ERROR in MModuleDepthCalibration2024::Initialize: Found a duplicate detector: "<= c_Error) { + cout<<"ERROR in MModuleDepthCalibration2024::Initialize: Found a duplicate detector: "<GetNSensitiveVolumes()<<" Sensitive Volumes."<= c_Error) { + cout<<"ERROR in MModuleDepthCalibration2024::Initialize: Found a Strip3D detector with "<GetNSensitiveVolumes()<<" Sensitive Volumes."<GetGuardRingVeto()==true) { - Event->SetDepthCalibrationIncomplete(); + Event->SetDepthCalibrationError(true, "GR Veto"); return false; } else { @@ -224,7 +231,7 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) // GRADE=-1 is an error. Break from the loop and continue. if ( Grade < 0 ){ H->SetNoDepth(); - Event->SetDepthCalibrationIncomplete(); + Event->SetDepthCalibrationError(); if (Grade == -1) { ++m_ErrorSH; } else if (Grade == -2) { @@ -234,7 +241,7 @@ bool MModuleDepthCalibration2024::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->SetDepthCalibrationIncomplete(); + Event->SetDepthCalibrationError(true, "Multiple hits on single strip"); if (Grade==5) { ++m_Error5; } else if (Grade==6) { @@ -285,24 +292,24 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) double LVTiming = LVSH->GetTiming(); double HVTiming = HVSH->GetTiming(); - // If there aren't coefficients loaded, then calibration is incomplete. + // If there aren't coefficients loaded, then report a depth calibration error. if( Coeffs == nullptr ){ //set the bad flag for depth H->SetNoDepth(); - Event->SetDepthCalibrationIncomplete(); + Event->SetDepthCalibrationError(true, "No calibration coefficients"); ++m_Error1; } else if (CTDVec.size() == 0) { cout << "Empty CTD vector" << endl; H->SetNoDepth(); - Event->SetDepthCalibrationIncomplete(); + Event->SetDepthCalibrationError(true, "No calibration coefficients"); } else if (DepthVec.size() == 0) { cout << "Empty Depth vector" << endl; H->SetNoDepth(); - Event->SetDepthCalibrationIncomplete(); + Event->SetDepthCalibrationError(true, "No calibration coefficients"); } else if ((LVTiming < 1.0E-6) || (HVTiming < 1.0E-6)) { ++m_Error3; H->SetNoDepth(); - Event->SetDepthCalibrationIncomplete(); + Event->SetDepthCalibrationError(true, "No timing"); } else { // If there are coefficients and timing information is loaded, try calculating the CTD and depth @@ -320,7 +327,7 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) //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->SetDepthCalibrationIncomplete(); + Event->SetDepthCalibrationError(true, "Out of Range"); ++m_Error2; } @@ -355,7 +362,7 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) Zpos = mean_depth; // Add the depth to the GUI histogram. - if (Event->IsStripPairingIncomplete()==false) { + if (Event->IsStripPairingError()==false) { if (HasExpos() == true) { m_ExpoDepthCalibration->AddDepth(DetID, Zpos); } @@ -475,7 +482,9 @@ std::vector* MModuleDepthCalibration2024::GetPixelCoeffs(int PixelCode) if( m_Coeffs.count(PixelCode) > 0 ){ return &m_Coeffs[PixelCode]; } else { - cout << "MModuleDepthCalibration2024::GetPixelCoeffs: cannot get stretch and offset; pixel code " << PixelCode << " not found." << endl; + if (g_Verbosity >= c_Warning) { + cout << "MModuleDepthCalibration2024::GetPixelCoeffs: cannot get stretch and offset; pixel code " << PixelCode << " not found." << endl; + } return nullptr; } } else { diff --git a/src/MModuleDepthCalibrationB.cxx b/src/MModuleDepthCalibrationB.cxx index 3b32494..c8e7892 100644 --- a/src/MModuleDepthCalibrationB.cxx +++ b/src/MModuleDepthCalibrationB.cxx @@ -252,7 +252,7 @@ bool MModuleDepthCalibrationB::AnalyzeEvent(MReadOutAssembly* Event) NewHits.push_back(NH); } -// Event->SetDepthCalibrationIncomplete(); //for testing purposes!!! filtering out events that are not 1x1 strip AWL1x1 +// Event->SetDepthCalibrationError(); //for testing purposes!!! filtering out events that are not 1x1 strip AWL1x1 } else if( (XStrips.size() == 2) && (YStrips.size() == 2) ){ //in this case use depth from dominant strips but use weighted X and Y positions @@ -278,7 +278,7 @@ bool MModuleDepthCalibrationB::AnalyzeEvent(MReadOutAssembly* Event) //GlobalPosition = m_Geometry->GetGlobalPosition( LocalPosition, m_DetectorNames[DetID]); GlobalPosition = m_DetectorVolumes[DetID]->GetPositionInWorldVolume(LocalPosition); H->SetPosition(GlobalPosition); H->SetPositionResolution(PositionResolution); -// Event->SetDepthCalibrationIncomplete(); //for testing purposes!!! filtering out events that are not 1x1 strip AWL1X1 +// Event->SetDepthCalibrationError(); //for testing purposes!!! filtering out events that are not 1x1 strip AWL1X1 } else { //set too many SH bad flag @@ -289,24 +289,24 @@ bool MModuleDepthCalibrationB::AnalyzeEvent(MReadOutAssembly* Event) //good ++m_NoError; } else if( Error == 1 ){ - Event->SetDepthCalibrationIncomplete(); + Event->SetDepthCalibrationError(); ++m_Error1; } else if( Error == 2 ){ - Event->SetDepthCalibrationIncomplete(); + Event->SetDepthCalibrationError(); ++m_Error2; } else if( Error == 3){ //Hits that were missing timing information //EHist->Fill(H->GetEnergy()); //don't set the globally bad flag - //Event->SetDepthCalibrationIncomplete(); -// Event->SetDepthCalibrationIncomplete(); //for testing purposes!!! filtering out events that are not 1x1 strip AWL1X1 + //Event->SetDepthCalibrationError(); +// Event->SetDepthCalibrationError(); //for testing purposes!!! filtering out events that are not 1x1 strip AWL1X1 ++m_Error3; } else if( Error == 4){ //hit was bad because of StripHitMultipleTimes flag from strip pairing ++m_Error4; - Event->SetDepthCalibrationIncomplete(); + Event->SetDepthCalibrationError(); } else if( Error == -1){ - Event->SetDepthCalibrationIncomplete(); + Event->SetDepthCalibrationError(); ++m_ErrorSH; } diff --git a/src/MModuleEnergyCalibrationUniversal.cxx b/src/MModuleEnergyCalibrationUniversal.cxx index 5cf37bf..8f16fd9 100644 --- a/src/MModuleEnergyCalibrationUniversal.cxx +++ b/src/MModuleEnergyCalibrationUniversal.cxx @@ -63,38 +63,38 @@ ClassImp(MModuleEnergyCalibrationUniversal) MModuleEnergyCalibrationUniversal::MModuleEnergyCalibrationUniversal() : MModule() { // Construct an instance of MModuleEnergyCalibrationUniversal - + // Set all module relevant information - + // Set the module name --- has to be unique m_Name = "Universal energy calibrator"; - + // Set the XML tag --- has to be unique --- no spaces allowed m_XmlTag = "EnergyCalibrationUniversal"; - + // Set all modules, which have to be done before this module AddPreceedingModuleType(MAssembly::c_EventLoader); // AddPreceedingModuleType(MAssembly::c_TACcut); - + // Set all types this modules handles AddModuleType(MAssembly::c_EnergyCalibration); - + // Set all modules, which can follow this module AddSucceedingModuleType(MAssembly::c_TACcut); - + // Set if this module has an options GUI m_HasOptionsGUI = true; - + // Allow the use of multiple threads and instances m_AllowMultiThreading = true; m_AllowMultipleInstances = true; - // + // Initiate the Slow Threshold Cut variables m_SlowThresholdCutMode = MSlowThresholdCutModes::e_Ignore; - m_SlowThresholdCutFixedValue = 35; + m_SlowThresholdCutFixedValue = 15; m_SlowThresholdCutFileName = ""; -} +} //////////////////////////////////////////////////////////////////////////////// @@ -111,12 +111,14 @@ MModuleEnergyCalibrationUniversal::~MModuleEnergyCalibrationUniversal() void MModuleEnergyCalibrationUniversal::CreateExpos() { // If they are already created, return - if (m_Expos.size() != 0) return; - + if (m_Expos.size() != 0) { + return; + } + // Set the histogram display m_ExpoEnergyCalibration = new MGUIExpoEnergyCalibration(this); m_ExpoEnergyCalibration->SetEnergyHistogramParameters(200, 0, 2000); - m_Expos.push_back(m_ExpoEnergyCalibration); + m_Expos.push_back(m_ExpoEnergyCalibration); } @@ -126,11 +128,13 @@ void MModuleEnergyCalibrationUniversal::CreateExpos() bool MModuleEnergyCalibrationUniversal::Initialize() { // Initialize the module - + //Parse the Energy Calibration file MParser Parser; if (Parser.Open(m_FileName, MFile::c_Read) == false) { - if (g_Verbosity >= c_Error) cout<= c_Error) { + cout << m_XmlTag << ": Unable to open calibration file " << m_FileName << endl; + } return false; } //Parse the slow Threshold file @@ -139,30 +143,27 @@ bool MModuleEnergyCalibrationUniversal::Initialize() if (Parser_Threshold.Open(m_SlowThresholdCutFileName, MFile::c_Read) == false) { if (g_Verbosity >= c_Error) cout< Tokens = Line.Tokenize(","); // for each line, Create tokens seperated by commas if (Tokens.size() == 6) { int IndexOffset = Tokens.size() % 6; //index counter - int DetID = Tokens[1+IndexOffset].ToInt(); // Detector ID - MString Side = Tokens[2+IndexOffset].ToString(); // side is a string, either 'l' or 'h' - int StripID = Tokens[3+IndexOffset].ToInt(); // stripID - int ThresholdADC = Tokens[4+IndexOffset].ToInt(); // energy threshold in ADC - double ThresholdKeVFile = Tokens[5+IndexOffset].ToDouble(); //energy threshold in keV - + int DetID = Tokens[1 + IndexOffset].ToInt(); // Detector ID + MString Side = Tokens[2 + IndexOffset].ToString(); // side is a string, either 'l' or 'h' + int StripID = Tokens[3 + IndexOffset].ToInt(); // stripID + //int ThresholdADC = Tokens[4 + IndexOffset].ToInt(); // energy threshold in ADC + double ThresholdKeVFile = Tokens[5 + IndexOffset].ToDouble(); //energy threshold in keV + MReadOutElementDoubleStrip R; R.SetDetectorID(DetID); R.SetStripID(StripID); R.IsLowVoltageStrip(Side == "l"); - + // map detectorID, strip number, and voltage side to the threshold (keV) m_ThresholdMap[R] = ThresholdKeVFile; - } - } } } @@ -175,11 +176,13 @@ bool MModuleEnergyCalibrationUniversal::Initialize() //tokenize ecal file for (unsigned int i = 0; i < Parser.GetNLines(); ++i) { unsigned int NTokens = Parser.GetTokenizerAt(i)->GetNTokens(); - if (NTokens < 2) continue; + if (NTokens < 2) { + continue; + } if (Parser.GetTokenizerAt(i)->IsTokenAt(0, "CP") == true || - Parser.GetTokenizerAt(i)->IsTokenAt(0, "CM") == true || Parser.GetTokenizerAt(i)->IsTokenAt(0,"CR") == true) { + Parser.GetTokenizerAt(i)->IsTokenAt(0, "CM") == true || Parser.GetTokenizerAt(i)->IsTokenAt(0, "CR") == true) { if (Parser.GetTokenizerAt(i)->IsTokenAt(1, "dss") == true) { - + // input token values to map MReadOutElementDoubleStrip R; R.SetDetectorID(Parser.GetTokenizerAt(i)->GetTokenAtAsUnsignedInt(2)); @@ -193,90 +196,96 @@ bool MModuleEnergyCalibrationUniversal::Initialize() CR_ROEToLine[R] = i; } } else { - if (g_Verbosity >= c_Error) cout<GetTokenAt(1)<<")"<= c_Error) { + cout << m_XmlTag << ": Line parser: Unknown read-out element (" << Parser.GetTokenizerAt(i)->GetTokenAt(1) << ")" << endl; + } return false; } } } - for (auto CM: CM_ROEToLine) { + for (auto CM : CM_ROEToLine) { // If we have at least three data points, we store the calibration - + if (CP_ROEToLine.find(CM.first) != CP_ROEToLine.end()) { unsigned int i = CP_ROEToLine[CM.first]; if (Parser.GetTokenizerAt(i)->IsTokenAt(5, "pakw") == false) { - if (g_Verbosity >= c_Warning) cout<GetTokenAt(5)<= c_Warning) { + cout << m_XmlTag << ": Unknown calibration point descriptor found: " << Parser.GetTokenizerAt(i)->GetTokenAt(5) << endl; + } continue; } } else { - if (g_Verbosity >= c_Warning) cout<= c_Warning) { + cout << m_XmlTag << ": No good calibration for the following strip found: " << CM.first << endl; + } continue; } - + // Read the calibrator, i.e. read the fit function from the .ecal Melinator file. - + unsigned int Pos = 5; MString CalibratorType = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsString(Pos); CalibratorType.ToLower(); - + // Below inclusion of poly1zero written by J. Beechert on 2019/11/15 if (CalibratorType == "poly1zero") { double a0 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); - + //From the fit parameters I just extracted from the .ecal file, I can define a function - TF1* melinatorfit = new TF1("poly1zero","0. + [0]*x", 0., 8191.); + TF1* melinatorfit = new TF1("poly1zero", "0. + [0]*x", 0., 8191.); melinatorfit->FixParameter(0, a0); - + //Define the map by saving the fit function I just created as a map to the current ReadOutElement m_Calibration[CM.first] = melinatorfit; - - } - + + } + // Below inclusion of poly1 and poly2 written by J. Beechert on 2019/10/24 else if (CalibratorType == "poly1") { double a0 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); double a1 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); - + //From the fit parameters I just extracted from the .ecal file, I can define a function - TF1* melinatorfit = new TF1("poly1","[0] + [1]*x", 0., 8191.); + TF1* melinatorfit = new TF1("poly1", "[0] + [1]*x", 0., 8191.); melinatorfit->FixParameter(0, a0); melinatorfit->FixParameter(1, a1); - + //Define the map by saving the fit function I just created as a map to the current ReadOutElement m_Calibration[CM.first] = melinatorfit; - + } else if (CalibratorType == "poly2") { double a0 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); double a1 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); double a2 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); - + //From the fit parameters I just extracted from the .ecal file, I can define a function - TF1* melinatorfit = new TF1("poly2","[0] + [1]*x + [2]*x^2", 0., 8191.); + TF1* melinatorfit = new TF1("poly2", "[0] + [1]*x + [2]*x^2", 0., 8191.); melinatorfit->FixParameter(0, a0); melinatorfit->FixParameter(1, a1); melinatorfit->FixParameter(2, a2); - + //Define the map by saving the fit function I just created as a map to the current ReadOutElement m_Calibration[CM.first] = melinatorfit; - - } - //Eventually, I'll be including other possible fits, but for now, we've just include poly3 and poly4 - else if (CalibratorType == "poly3") { + + } + //Eventually, I'll be including other possible fits, but for now, we've just include poly3 and poly4 + else if (CalibratorType == "poly3") { double a0 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); double a1 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); double a2 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); double a3 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); - + //From the fit parameters I just extracted from the .ecal file, I can define a function - TF1* melinatorfit = new TF1("poly3","[0] + [1]*x + [2]*x^2 + [3]*x^3", 0., 8191.); + TF1* melinatorfit = new TF1("poly3", "[0] + [1]*x + [2]*x^2 + [3]*x^3", 0., 8191.); melinatorfit->FixParameter(0, a0); melinatorfit->FixParameter(1, a1); melinatorfit->FixParameter(2, a2); melinatorfit->FixParameter(3, a3); - + //Define the map by saving the fit function I just created as a map to the current ReadOutElement m_Calibration[CM.first] = melinatorfit; - + } else if (CalibratorType == "poly4") { double a0 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); double a1 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); @@ -285,7 +294,7 @@ bool MModuleEnergyCalibrationUniversal::Initialize() double a4 = Parser.GetTokenizerAt(CM.second)->GetTokenAtAsDouble(++Pos); //From the fit parameters I just extracted from the .ecal file, I can define a function - TF1 * melinatorfit = new TF1("poly4","[0] + [1]*x + [2]*x^2 + [3]*x^3 + [4]*x^4", 0., 8191.); + TF1* melinatorfit = new TF1("poly4", "[0] + [1]*x + [2]*x^2 + [3]*x^3 + [4]*x^4", 0., 8191.); melinatorfit->FixParameter(0, a0); melinatorfit->FixParameter(1, a1); melinatorfit->FixParameter(2, a2); @@ -296,12 +305,14 @@ bool MModuleEnergyCalibrationUniversal::Initialize() m_Calibration[CM.first] = melinatorfit; } else { - if (g_Verbosity >= c_Error) cout<= c_Error) { + cout << m_XmlTag << ": Line parser: Unknown calibrator type (" << CalibratorType << ") for strip" << CM.first << endl; + } continue; } } - - for (auto CR: CR_ROEToLine) { + + for (auto CR : CR_ROEToLine) { unsigned int Pos = 5; MString CalibratorType = Parser.GetTokenizerAt(CR.second)->GetTokenAtAsString(Pos); @@ -309,97 +320,120 @@ bool MModuleEnergyCalibrationUniversal::Initialize() if (CalibratorType == "p1" || CalibratorType == "poly1") { double f0 = Parser.GetTokenizerAt(CR.second)->GetTokenAtAsDouble(++Pos); double f1 = Parser.GetTokenizerAt(CR.second)->GetTokenAtAsDouble(++Pos); - TF1* resolutionfit = new TF1("P1","([0]+[1]*x) / 2.355",0.,2000.); - resolutionfit->FixParameter(0,f0); - resolutionfit->FixParameter(1,f1); + TF1* resolutionfit = new TF1("P1", "([0]+[1]*x) / 2.355", 0., 2000.); + resolutionfit->FixParameter(0, f0); + resolutionfit->FixParameter(1, f1); m_ResolutionCalibration[CR.first] = resolutionfit; - } - else if (CalibratorType == "p2" || CalibratorType == "poly2") { + } else if (CalibratorType == "p2" || CalibratorType == "poly2") { double f0 = Parser.GetTokenizerAt(CR.second)->GetTokenAtAsDouble(++Pos); double f1 = Parser.GetTokenizerAt(CR.second)->GetTokenAtAsDouble(++Pos); double f2 = Parser.GetTokenizerAt(CR.second)->GetTokenAtAsDouble(++Pos); - TF1* resolutionfit = new TF1("P2","([0]+[1]*x+[2]*x*x) / 2.355",0.,2000.); - resolutionfit->FixParameter(0,f0); - resolutionfit->FixParameter(1,f1); - resolutionfit->FixParameter(2,f2); + TF1* resolutionfit = new TF1("P2", "([0]+[1]*x+[2]*x*x) / 2.355", 0., 2000.); + resolutionfit->FixParameter(0, f0); + resolutionfit->FixParameter(1, f1); + resolutionfit->FixParameter(2, f2); m_ResolutionCalibration[CR.first] = resolutionfit; } else { - if (g_Verbosity >= c_Error) cout<= c_Error) { + cout << m_XmlTag << ": Line parser: Unknown resolution calibrator type (" << CalibratorType << ") for strip" << CR.first << endl; + } continue; } } return MModule::Initialize(); - } //////////////////////////////////////////////////////////////////////////////// -bool MModuleEnergyCalibrationUniversal::AnalyzeEvent(MReadOutAssembly* Event) +bool MModuleEnergyCalibrationUniversal::AnalyzeEvent(MReadOutAssembly* Event) { // Main data analysis routine, which updates the event to a new level, i.e. takes the raw ADC value from the .roa file loaded through nuclearizer and converts it into energy units. - - for (unsigned int i = 0; i < Event->GetNStripHits(); ++i) { + + for (unsigned int i = 0; i < Event->GetNStripHits();) { + MStripHit* SH = Event->GetStripHit(i); MReadOutElementDoubleStrip R = *dynamic_cast(SH->GetReadOutElement()); - + TF1* Fit = m_Calibration[R]; TF1* FitRes = m_ResolutionCalibration[R]; - double ADCMod, newADC; if (Fit == nullptr) { - if (g_Verbosity >= c_Error) cout<SetEnergyCalibrationIncomplete_BadStrip(true); - } else { + if (g_Verbosity >= c_Error) { + cout << m_XmlTag << ": Error: Energy-fit not found for read-out element " << R << endl; + } + Event->SetEnergyCalibrationError(true); + ++i; // iterate to next SH + continue; + + } else { double Energy = 0; // declare energy variable + double Threshold = 0; // 0 means ignore Energy = Fit->Eval(SH->GetADCUnits()); - double Threshold = 0; // 0 means ignore + + if (Energy < 0) { + Energy = 0; // If the calibrated energy is less than 0, force it to be 0. + } if (m_SlowThresholdCutMode == MSlowThresholdCutModes::e_Fixed) { //check if user input threshold is enabled (one value applied to all strips) Threshold = m_SlowThresholdCutFixedValue; } else if (m_SlowThresholdCutMode == MSlowThresholdCutModes::e_File) { //check if threshold file is enabled (unique value applied to each strip) - double Threshold_map = m_ThresholdMap[R]; // if file enabled, declare value from map - - if (Threshold_map == 0) { - if (g_Verbosity >= c_Error) cout<= c_Error) { + cout << m_XmlTag << ": Error: Threshold not found for read-out element " << R << endl; + } + const double DefaultSlowThreshold = 15.0; + Threshold = DefaultSlowThreshold; // set default threshold if threshold not found } else { - Threshold = Threshold_map; // set threshold variable to value found in map + Threshold = ThresholdFromFile; // set threshold variable to value found in map } } - - //! Remove SH for any energy value below the established threshold + + //! Remove SH for any energy value below the established threshold (0 is default) if (Energy < Threshold) { + if (g_Verbosity >= c_Warning) { + cout << m_XmlTag << ": Strip Hit below threshold, deleting SH with Energy " << Energy << " keV " << endl; + } + Event->AddStripHitBelowThreshold(Energy); + Event->SetStripHitBelowThreshold_QualityFlag(true, to_string(Event->GetStripHitBelowThreshold_Number()) + " SH(s) with total energy " + to_string(Event->GetStripHitBelowThreshold_Energy())); Event->RemoveStripHit(i); - if (g_Verbosity >= c_Error) cout<SetStripHitBelowThreshold(true); - return false; - } - - SH->SetEnergy(Energy); - if (FitRes == nullptr) { - if (g_Verbosity >= c_Error) cout<SetEnergyResolutionCalibrationIncomplete(true); + delete SH; + continue; // continue to next SH without iterating i } else { - double EnergyResolution = FitRes->Eval(Energy); - SH->SetEnergyResolution(EnergyResolution); - } - if (R.IsLowVoltageStrip() == true) { // check voltage side - if (HasExpos() == true) { - m_ExpoEnergyCalibration->AddEnergy(Energy); + + // if the energy isn't filtered out with the threshold, then assign the energy to the SH + ++i; // iterate to next SH + SH->SetEnergy(Energy); + + if (FitRes == nullptr) { + if (g_Verbosity >= c_Error) { + cout << m_XmlTag << ": Error: Energy Resolution fit not found for read-out element " << R << endl; + } + // There is not expected to be a time in which the energy resolution calibration is not defined when the energy calibration itself is. Therefore, don't need a seperate BD flag for this. + } else { + double EnergyResolution = FitRes->Eval(Energy); + SH->SetEnergyResolution(EnergyResolution); + } + if (R.IsLowVoltageStrip() == true) { // check voltage side to plot only LV hits to the Expo histogram + if (HasExpos() == true) { + m_ExpoEnergyCalibration->AddEnergy(Energy); + } + } + if (g_Verbosity >= c_Info) { + cout << m_XmlTag << ": Energy: " << SH->GetADCUnits() << " adc --> " << Energy << " keV" << endl; } } - if (g_Verbosity >= c_Info) cout<GetADCUnits()<<" adu --> "<SetAnalysisProgress(MAssembly::c_EnergyCalibration); - + return true; } @@ -419,7 +453,7 @@ double MModuleEnergyCalibrationUniversal::GetEnergy(MReadOutElementDoubleStrip R Energy = 0.0; } } else { - cout<GetX(Energy); } else { - cout<Create(); gClient->WaitForUnmap(Options); @@ -479,7 +517,7 @@ void MModuleEnergyCalibrationUniversal::ShowOptionsGUI() bool MModuleEnergyCalibrationUniversal::ReadXmlConfiguration(MXmlNode* Node) { //! Read the configuration data from an XML node - + MXmlNode* FileNameNode = Node->GetNode("FileName"); if (FileNameNode != nullptr) { m_FileName = FileNameNode->GetValue(); @@ -497,7 +535,7 @@ bool MModuleEnergyCalibrationUniversal::ReadXmlConfiguration(MXmlNode* Node) if (SlowThresholdCutFileNameNode != nullptr) { m_SlowThresholdCutFileName = SlowThresholdCutFileNameNode->GetValue(); } - + return true; } @@ -505,10 +543,10 @@ bool MModuleEnergyCalibrationUniversal::ReadXmlConfiguration(MXmlNode* Node) //////////////////////////////////////////////////////////////////////////////// -MXmlNode* MModuleEnergyCalibrationUniversal::CreateXmlConfiguration() +MXmlNode* MModuleEnergyCalibrationUniversal::CreateXmlConfiguration() { //! Create an XML node tree from the configuration - + MXmlNode* Node = new MXmlNode(0, m_XmlTag); new MXmlNode(Node, "FileName", m_FileName); new MXmlNode(Node, "SlowThresholdCutMode", static_cast(m_SlowThresholdCutMode)); @@ -527,12 +565,12 @@ double MModuleEnergyCalibrationUniversal::LookupEnergyResolution(MStripHit* SH, MReadOutElementDoubleStrip* ROE = dynamic_cast(SH->GetReadOutElement()); if (ROE == nullptr) { - cout<Eval(Energy); @@ -540,6 +578,5 @@ double MModuleEnergyCalibrationUniversal::LookupEnergyResolution(MStripHit* SH, } - // MModuleEnergyCalibrationUniversal.cxx: the end... //////////////////////////////////////////////////////////////////////////////// diff --git a/src/MModuleRevan.cxx b/src/MModuleRevan.cxx index 56f6d8b..5104b6a 100644 --- a/src/MModuleRevan.cxx +++ b/src/MModuleRevan.cxx @@ -185,10 +185,10 @@ bool MModuleRevan::AnalyzeEvent(MReadOutAssembly* Event) MPhysicalEvent* PE = BestRawEvent->GetPhysicalEvent(); Event->SetPhysicalEvent(PE); } else { - Event->SetEventReconstructionIncomplete(true); + Event->SetEventReconstructionError(true); } } else { - Event->SetEventReconstructionIncomplete(true); + Event->SetEventReconstructionError(true); } Event->SetAnalysisProgress(MAssembly::c_EventReconstruction); diff --git a/src/MModuleStripPairingChiSquare.cxx b/src/MModuleStripPairingChiSquare.cxx index 7237a74..4611215 100644 --- a/src/MModuleStripPairingChiSquare.cxx +++ b/src/MModuleStripPairingChiSquare.cxx @@ -190,7 +190,7 @@ bool MModuleStripPairingChiSquare::AnalyzeEvent(MReadOutAssembly* Event) unsigned int MaxCombinations = 5; if (Event->GetNStripHits() == 0) { - Event->SetStripPairingIncomplete(true, "No strip hits"); + Event->SetStripPairingError(true, "No strip hits"); Event->SetAnalysisProgress(MAssembly::c_StripPairing); return false; } @@ -230,7 +230,7 @@ bool MModuleStripPairingChiSquare::AnalyzeEvent(MReadOutAssembly* Event) for (unsigned int d = 0; d < StripHits.size(); ++d) { // Detector loop for (unsigned int side = 0; side <=1; ++side) { // side loop if (StripHits[d][side].size() > MaxStripHits) { - Event->SetStripPairingIncomplete(true, "More than 6 hit strIps on one side"); + Event->SetStripPairingError(true, "More than 6 hit strIps on one side"); Event->SetAnalysisProgress(MAssembly::c_StripPairing); return false; } @@ -260,7 +260,7 @@ bool MModuleStripPairingChiSquare::AnalyzeEvent(MReadOutAssembly* Event) for (unsigned int d = 0; d < StripHits.size(); ++d) { // Detector loop if (StripHits[d][0].size() == 0 || StripHits[d][1].size() == 0) { - Event->SetStripPairingIncomplete(true, "One detector side has not strip hits"); + Event->SetStripPairingError(true, "One detector side has not strip hits"); Event->SetAnalysisProgress(MAssembly::c_StripPairing); return false; } @@ -470,7 +470,7 @@ bool MModuleStripPairingChiSquare::AnalyzeEvent(MReadOutAssembly* Event) // Now create hits: if (BestChiSquare == numeric_limits::max()) { - Event->SetStripPairingIncomplete(true, "Pairing did not find a single match"); + Event->SetStripPairingError(true, "Pairing did not find a single match"); Event->SetAnalysisProgress(MAssembly::c_StripPairing); return false; } @@ -543,7 +543,7 @@ bool MModuleStripPairingChiSquare::AnalyzeEvent(MReadOutAssembly* Event) } if (EnergyTotal > max(XEnergyTotal, YEnergyTotal) + 2.5*max(XEnergyRes, YEnergyRes) || EnergyTotal < min(XEnergyTotal, YEnergyTotal) - 2.5*max(XEnergyRes, YEnergyRes)) { - Event->SetStripPairingIncomplete(true, "Strips not pairable wihin 2.5 sigma of measure denergy"); + Event->SetStripPairingError(true, "Strips not pairable wihin 2.5 sigma of measure denergy"); Event->SetAnalysisProgress(MAssembly::c_StripPairing); return false; } diff --git a/src/MModuleStripPairingGreedy.cxx b/src/MModuleStripPairingGreedy.cxx index a697a94..f55708c 100644 --- a/src/MModuleStripPairingGreedy.cxx +++ b/src/MModuleStripPairingGreedy.cxx @@ -333,7 +333,7 @@ bool MModuleStripPairingGreedy::AnalyzeEvent(MReadOutAssembly* Event){ //flag hits without enough strips for (int det=0; det<12; det++){ if (notEnoughStrips[det] == 1 || notEnoughStrips[det] == 2){ - Event->SetStripPairingIncomplete(true,"bad number of strip hits"); + Event->SetStripPairingError(true,"bad number of strip hits"); break; } } @@ -341,19 +341,19 @@ bool MModuleStripPairingGreedy::AnalyzeEvent(MReadOutAssembly* Event){ for (unsigned int h = 0; h < Event->GetNHits(); h++){ if (Event->GetHit(h)->GetStripHitMultipleTimesX() == true || Event->GetHit(h)->GetStripHitMultipleTimesY() == true){ - Event->SetStripPairingIncomplete(true,"multiple hits per strip"); + Event->SetStripPairingError(true,"multiple hits per strip"); } /* if (Event->GetHit(h)->GetChargeSharing() == true){ - Event->SetStripPairingIncomplete(true,"charge sharing"); + Event->SetStripPairingError(true,"charge sharing"); } */ /* int detID = Event->GetHit(h)->GetStripHit(0)->GetDetectorID(); if (chi_sq[detID] <= 25){ - Event->SetStripPairingIncomplete(true,"good pairing, multiple hits per strip"); + Event->SetStripPairingError(true,"good pairing, multiple hits per strip"); } else { - Event->SetStripPairingIncomplete(true,"bad pairing, multiple hits per strip"); + Event->SetStripPairingError(true,"bad pairing, multiple hits per strip"); } */ } @@ -399,12 +399,12 @@ bool MModuleStripPairingGreedy::AnalyzeEvent(MReadOutAssembly* Event){ */ // Difference must be more than 10 keV for cross talk + 2 sigma energy resolution on *both* sides // if (Event->GetHit(h)->GetStripHitMultipleTimes() == true){ -// Event->SetStripPairingIncomplete(true,"multiple hits per strip"); +// Event->SetStripPairingError(true,"multiple hits per strip"); // } if (Difference > 2*pUncertainty + 2*nUncertainty + 20) { if (Event->GetHit(h)->GetStripHitMultipleTimesX() == false && Event->GetHit(h)->GetStripHitMultipleTimesY() == false){ - Event->SetStripPairingIncomplete(true, "bad pairing"); + Event->SetStripPairingError(true, "bad pairing"); } /* if (nHits[0] != 0 && nHits[1] != 0){ PrintXYStripsHitOrig(); @@ -427,7 +427,7 @@ bool MModuleStripPairingGreedy::AnalyzeEvent(MReadOutAssembly* Event){ // } // else{ -// Event->SetStripPairingIncomplete(true,"multiple hits per strip"); +// Event->SetStripPairingError(true,"multiple hits per strip"); // } @@ -438,7 +438,7 @@ bool MModuleStripPairingGreedy::AnalyzeEvent(MReadOutAssembly* Event){ } //flag events with charge sharing // if (Event->GetHit(h)->GetChargeSharing()){ -// Event->SetStripPairingIncomplete(true,"charge sharing"); +// Event->SetStripPairingError(true,"charge sharing"); // } } // } @@ -455,7 +455,7 @@ bool MModuleStripPairingGreedy::AnalyzeEvent(MReadOutAssembly* Event){ } /* - if (Event->GetNHits()==0 && Event->IsStripPairingIncomplete() == false){ + if (Event->GetNHits()==0 && Event->IsStripPairingError() == false){ for (unsigned int sh=0; shGetNStripHits(); sh++){ MStripHit* striphit = Event->GetStripHit(sh); cout << striphit->GetDetectorID() << '\t'; @@ -468,7 +468,7 @@ bool MModuleStripPairingGreedy::AnalyzeEvent(MReadOutAssembly* Event){ */ if (m_StripPairingFailed != "") { - Event->SetStripPairingIncomplete(true, m_StripPairingFailed); + Event->SetStripPairingError(true, m_StripPairingFailed); } Event->SetAnalysisProgress(MAssembly::c_StripPairing); diff --git a/src/MModuleStripPairingMultiRoundChiSquare.cxx b/src/MModuleStripPairingMultiRoundChiSquare.cxx index a3e42f1..2f9ca65 100644 --- a/src/MModuleStripPairingMultiRoundChiSquare.cxx +++ b/src/MModuleStripPairingMultiRoundChiSquare.cxx @@ -272,14 +272,14 @@ bool MModuleStripPairingMultiRoundChiSquare::EventSelection(MReadOutAssembly* Ev for (unsigned int d = 0; d < StripHits.size(); ++d) { // Detector loop for (unsigned int side = 0; side <= 1; ++side) { // Side loop if (StripHits[d][side].size() > MaxStripHits) { - Event->SetStripPairingIncomplete(true, "More than 6 hit strips on one side"); + Event->SetStripPairingError(true, "More than 6 hit strips on one side"); Event->SetAnalysisProgress(MAssembly::c_StripPairing); return false; } // Check if one side of the detector has no strip hits if (StripHits[d][side].size() == 0) { - Event->SetStripPairingIncomplete(true, "One detector side has no strip hits"); + Event->SetStripPairingError(true, "One detector side has no strip hits"); Event->SetAnalysisProgress(MAssembly::c_StripPairing); return false; } @@ -649,7 +649,7 @@ bool MModuleStripPairingMultiRoundChiSquare::CreateHits(unsigned int d, MReadOut // If both HV and LV have multiple hits per strip, can't pair else if (AllAdjacentLV == false and AllAdjacentHV == false) { - Event->SetStripPairingIncomplete(true, "Strips not pairable. Multiple hits per strip on LV and HV sides"); + Event->SetStripPairingError(true, "Strips not pairable. Multiple hits per strip on LV and HV sides"); Event->SetAnalysisProgress(MAssembly::c_StripPairing); return false; } @@ -660,7 +660,7 @@ bool MModuleStripPairingMultiRoundChiSquare::CreateHits(unsigned int d, MReadOut // One last quality selection based on total event energies if ((EnergyTotal > max(LVEnergyTotal, HVEnergyTotal) + 2.5 * max(LVEnergyResTotal, HVEnergyResTotal) || EnergyTotal < min(LVEnergyTotal, HVEnergyTotal) - 2.5 * max(LVEnergyResTotal, HVEnergyResTotal))) { - Event->SetStripPairingIncomplete(true, "Strips not pairable wihin 2.5 sigma of measured energy"); + Event->SetStripPairingError(true, "Strips not pairable wihin 2.5 sigma of measured energy"); Event->SetAnalysisProgress(MAssembly::c_StripPairing); return false; } @@ -695,7 +695,7 @@ bool MModuleStripPairingMultiRoundChiSquare::AnalyzeEvent(MReadOutAssembly* Even // Check if there are actually any strip hits if (Event->GetNStripHits() == 0) { - Event->SetStripPairingIncomplete(true, "No strip hits"); + Event->SetStripPairingError(true, "No strip hits"); Event->SetAnalysisProgress(MAssembly::c_StripPairing); return false; } @@ -759,13 +759,13 @@ bool MModuleStripPairingMultiRoundChiSquare::AnalyzeEvent(MReadOutAssembly* Even } // Check if chi^2 was ever actually updated if (BestChiSquare == numeric_limits::max()) { - Event->SetStripPairingIncomplete(true, "Pairing did not find a single match"); + Event->SetStripPairingError(true, "Pairing did not find a single match"); Event->SetAnalysisProgress(MAssembly::c_StripPairing); return false; } // Flag events with a reduced chi square > 25 else if (BestChiSquare > 25) { - Event->SetStripPairingIncomplete(true, "Best reduced chi square is not below 25"); + Event->SetStripPairingError(true, "Best reduced chi square is not below 25"); Event->SetAnalysisProgress(MAssembly::c_StripPairing); return false; } diff --git a/src/MReadOutAssembly.cxx b/src/MReadOutAssembly.cxx index 1d0558a..a78af52 100644 --- a/src/MReadOutAssembly.cxx +++ b/src/MReadOutAssembly.cxx @@ -53,7 +53,7 @@ MReadOutAssembly::MReadOutAssembly() : MReadOutSequence(), m_EventTimeUTC(0) m_PhysicalEvent = nullptr; m_SimEvent = nullptr; m_Aspect = nullptr; - m_HasSimAspectInfo = false; + m_HasSimAspectInfo = false; Clear(); } @@ -76,6 +76,11 @@ MReadOutAssembly::~MReadOutAssembly() } m_StripHitsTOnly.clear(); + // Delete all DEE Strip hits + m_DEEStripHitsLV.clear(); + m_DEEStripHitsHV.clear(); + m_DEECrystalHits.clear(); + // Delete all crystal hits for (unsigned int h = 0; h < m_CrystalHits.size(); ++h) { delete m_CrystalHits[h]; @@ -94,17 +99,13 @@ MReadOutAssembly::~MReadOutAssembly() } m_HitsSim.clear(); - // Delete all guardring hits for (unsigned int h = 0; h < m_GuardringHits.size(); ++h) { delete m_GuardringHits[h]; } m_GuardringHits.clear(); - m_DEEStripHitsLV.clear(); - m_DEEStripHitsHV.clear(); - m_DEECrystalHits.clear(); - + // Delete all Events delete m_SimEvent; delete m_PhysicalEvent; delete m_Aspect; @@ -127,7 +128,6 @@ void MReadOutAssembly::Clear() m_Time = 0; m_EventTimeUTC = 0; m_MJD = 0.0; - m_ReducedChiSquare = -1; m_ShieldVeto = false; m_GuardRingVeto = false; @@ -172,26 +172,22 @@ void MReadOutAssembly::Clear() } m_GuardringHits.clear(); - m_AspectIncomplete = false; - m_AspectIncompleteString = ""; - m_TimeIncomplete = false; - m_TimeIncompleteString = ""; - m_EnergyCalibrationIncomplete_BadStrip = false; - m_EnergyCalibrationIncomplete_BadStripString = ""; - m_EnergyCalibrationIncomplete = false; - m_EnergyCalibrationIncompleteString = ""; - m_EnergyResolutionCalibrationIncomplete = false; - m_EnergyResolutionCalibrationIncompleteString = ""; - m_StripPairingIncomplete = false; - m_StripPairingIncompleteString = ""; - m_LLDEvent = false; - m_LLDEventString = ""; - m_DepthCalibrationIncomplete = false; - m_DepthCalibrationIncompleteString = ""; - m_DepthCalibration_OutofRange = false; - m_DepthCalibration_OutofRangeString = ""; - m_EventReconstructionIncomplete = false; - m_EventReconstructionIncompleteString = ""; + // Delete all event flags and associated variables + m_EnergyCalibrationError = false; + m_EnergyCalibrationErrorString = ""; + m_StripPairingError = false; + m_StripPairingErrorString = ""; + m_DepthCalibrationError = false; + m_DepthCalibrationErrorString = ""; + m_EventReconstructionError = false; + m_EventReconstructionErrorString = ""; + + m_ReducedChiSquare = -1; + + m_StripHitBelowThreshold_QualityFlag = false; + m_StripHitBelowThresholdString_QualityFlag = ""; + m_StripHitBelowThreshold_Energy = 0; + m_StripHitBelowThreshold_Number = 0; m_FilteredOut = false; @@ -207,7 +203,7 @@ void MReadOutAssembly::Clear() delete m_Aspect; m_Aspect = nullptr; -} +}// //////////////////////////////////////////////////////////////////////////////// @@ -672,56 +668,32 @@ void MReadOutAssembly::StreamTra(ostream& S) void MReadOutAssembly::StreamBDFlags(ostream& S) { - // Stream the BD flags + // Stream the BD and QA flags - if (m_AspectIncomplete == true) { - S<<"BD AspectIncomplete"; - if (m_AspectIncompleteString != "") S<<" ("<