From 4b267f1a6e086b375168fcb12aa942a7798ec189 Mon Sep 17 00:00:00 2001 From: Junior Pena Date: Fri, 12 Sep 2025 11:54:28 -0400 Subject: [PATCH 01/16] Files created for cyclotron frequency reconstruction of cavity events: KTProcessedCavityEventData, KTCavityEventProcessing, and changes to KTROOTTreeTypeWriterEventAnalysis for writing the output to tree files. --- Source/Data/CMakeLists.txt | 2 + .../KTProcessedCavityEventData.cc | 54 ++++++++++++++ .../KTProcessedCavityEventData.hh | 40 ++++++++++ Source/EventAnalysis/CMakeLists.txt | 2 + .../EventAnalysis/KTCavityEventProcessing.cc | 47 ++++++++++++ .../EventAnalysis/KTCavityEventProcessing.hh | 65 +++++++++++++++++ .../KTROOTTreeTypeWriterEventAnalysis.cc | 73 +++++++++++++++++++ .../KTROOTTreeTypeWriterEventAnalysis.hh | 19 +++++ 8 files changed, 302 insertions(+) create mode 100644 Source/Data/EventAnalysis/KTProcessedCavityEventData.cc create mode 100644 Source/Data/EventAnalysis/KTProcessedCavityEventData.hh create mode 100644 Source/EventAnalysis/KTCavityEventProcessing.cc create mode 100644 Source/EventAnalysis/KTCavityEventProcessing.hh diff --git a/Source/Data/CMakeLists.txt b/Source/Data/CMakeLists.txt index 52fd90b83..ce4f5f879 100644 --- a/Source/Data/CMakeLists.txt +++ b/Source/Data/CMakeLists.txt @@ -30,6 +30,7 @@ set (DATA_HEADERFILES EventAnalysis/KTMultiTrackEventData.hh EventAnalysis/KTPowerFitData.hh EventAnalysis/KTProcessedMPTData.hh + EventAnalysis/KTProcessedCavityEventData.hh EventAnalysis/KTProcessedTrackData.hh EventAnalysis/KTRPTrackData.hh EventAnalysis/KTSequentialLineData.hh @@ -96,6 +97,7 @@ set (DATA_SOURCEFILES EventAnalysis/KTMultiTrackEventData.cc EventAnalysis/KTPowerFitData.cc EventAnalysis/KTProcessedMPTData.cc + EventAnalysis/KTProcessedCavityEventData.cc EventAnalysis/KTProcessedTrackData.cc EventAnalysis/KTRPTrackData.cc EventAnalysis/KTSequentialLineData.cc diff --git a/Source/Data/EventAnalysis/KTProcessedCavityEventData.cc b/Source/Data/EventAnalysis/KTProcessedCavityEventData.cc new file mode 100644 index 000000000..1ff72c9f5 --- /dev/null +++ b/Source/Data/EventAnalysis/KTProcessedCavityEventData.cc @@ -0,0 +1,54 @@ +/* + * KTProcessedCavityEventData.cc + * + * Created on: Sep 11, 2025 + * Author: juniorpe + */ + +#include "KTProcessedCavityEventData.hh" + +#include "KTLogger.hh" + +namespace Katydid +{ + + const std::string KTProcessedCavityEventData::sName("proc-cavity-event"); + + KTProcessedCavityEventData::KTProcessedCavityEventData() : + KTExtensibleData< KTProcessedCavityEventData >(), + fComponent(0), + fAcquisitionID(0.), + fEventID(0), + fTotalEventSequences(0), + fInitialCyclotronFrequency(0.) + { + } + + KTProcessedCavityEventData::KTProcessedCavityEventData(const KTProcessedCavityEventData& orig) : + KTExtensibleData< KTProcessedCavityEventData >(orig), + + fComponent(orig.fComponent), + fAcquisitionID(orig.fAcquisitionID), + fEventID(orig.fEventID), + fTotalEventSequences(orig.fTotalEventSequences), + fInitialCyclotronFrequency(orig.fInitialCyclotronFrequency) + { + } + + KTProcessedCavityEventData::~KTProcessedCavityEventData() + { + } + + KTProcessedCavityEventData& KTProcessedCavityEventData::operator=(const KTProcessedCavityEventData& rhs) + { + KTExtensibleData< KTProcessedCavityEventData >::operator=(rhs); + fComponent = rhs.fComponent; + fAcquisitionID = rhs.fAcquisitionID; + fEventID = rhs.fEventID; + fTotalEventSequences = rhs.fTotalEventSequences; + + fInitialCyclotronFrequency = rhs.fInitialCyclotronFrequency; + return *this; + } + +} diff --git a/Source/Data/EventAnalysis/KTProcessedCavityEventData.hh b/Source/Data/EventAnalysis/KTProcessedCavityEventData.hh new file mode 100644 index 000000000..1bcbd7c24 --- /dev/null +++ b/Source/Data/EventAnalysis/KTProcessedCavityEventData.hh @@ -0,0 +1,40 @@ +/* + * KTProcessedCavityEventData.hh + * + * Created on: Sep 11, 2025 + * Author: juniorpe + */ + +#ifndef KTPROCESSEDCAVITYEVENTDATA_HH_ +#define KTPROCESSEDCAVITYEVENTDATA_HH_ + +#include "KTData.hh" +#include "KTMultiTrackEventData.hh" +#include "KTMemberVariable.hh" + +namespace Katydid +{ + + class KTProcessedCavityEventData : public Nymph::KTExtensibleData< KTProcessedCavityEventData > + { + public: + KTProcessedCavityEventData(); + KTProcessedCavityEventData(const KTProcessedCavityEventData& orig); + virtual ~KTProcessedCavityEventData(); + + KTProcessedCavityEventData& operator=(const KTProcessedCavityEventData& rhs); + + MEMBERVARIABLE(unsigned, Component); + MEMBERVARIABLE(uint64_t, AcquisitionID); + MEMBERVARIABLE(unsigned, EventID); + MEMBERVARIABLE(unsigned, TotalEventSequences); + + MEMBERVARIABLE(double, InitialCyclotronFrequency); + + public: + static const std::string sName; + }; + +} +#endif + diff --git a/Source/EventAnalysis/CMakeLists.txt b/Source/EventAnalysis/CMakeLists.txt index bea0b58c0..d78bfb120 100644 --- a/Source/EventAnalysis/CMakeLists.txt +++ b/Source/EventAnalysis/CMakeLists.txt @@ -34,6 +34,7 @@ set (EVENTANALYSIS_HEADERFILES KTIterativeTrackClustering.hh KTMultiPeakEventBuilder.hh KTMultiPeakTrackProcessing.hh + KTCavityEventProcessing.hh KTMultiPeakTrackBuilder.hh KTMultiSliceClustering.hh KTOverlappingTrackClustering.hh @@ -73,6 +74,7 @@ set (EVENTANALYSIS_SOURCEFILES KTIterativeTrackClustering.cc KTMultiPeakEventBuilder.cc KTMultiPeakTrackProcessing.cc + KTCavityEventProcessing.cc KTMultiPeakTrackBuilder.cc KTMultiSliceClustering.cc KTOverlappingTrackClustering.cc diff --git a/Source/EventAnalysis/KTCavityEventProcessing.cc b/Source/EventAnalysis/KTCavityEventProcessing.cc new file mode 100644 index 000000000..d5a734321 --- /dev/null +++ b/Source/EventAnalysis/KTCavityEventProcessing.cc @@ -0,0 +1,47 @@ +/* + * KTCavityEventProcessing.cc + * + * Created on: Sep 11, 2025 + * Author: juniorpe + */ + +#include "KTCavityEventProcessing.hh" + +#include "KTLogger.hh" +#include "KTProcessedTrackData.hh" +#include "KTMultiTrackEventData.hh" + + +namespace Katydid +{ + KTLOGGER(evlog, "KTCavityEventProcessing"); + + // Register the processor + KT_REGISTER_PROCESSOR(KTCavityEventProcessing, "cavity-event-processing"); + + KTCavityEventProcessing::KTCavityEventProcessing(const std::string& name) : + KTProcessor(name), + fProcessedCavityEventSignal("proc-cavity-event", this), + fEventSlot("mt-event", this, &KTCavityEventProcessing::AnalyzeEvent, &fProcessedCavityEventSignal) + { + } + + KTCavityEventProcessing::~KTCavityEventProcessing() + { + } + + bool KTCavityEventProcessing::Configure(const scarab::param_node* node) + { + if (node == NULL) return false; + + return true; + } + + bool KTCavityEventProcessing::AnalyzeEvent( KTMultiTrackEventData& mtEventData ) + { + KTDEBUG(evlog, "Entered Analyze Event Function"); + + return true; + } + +} // namespace Katydid diff --git a/Source/EventAnalysis/KTCavityEventProcessing.hh b/Source/EventAnalysis/KTCavityEventProcessing.hh new file mode 100644 index 000000000..8a92f8cf2 --- /dev/null +++ b/Source/EventAnalysis/KTCavityEventProcessing.hh @@ -0,0 +1,65 @@ +/* + * KTCavityEventProcessing.hh + * + * Created on: Sep 11, 2025 + * Author: juniorpe + */ + +#ifndef KTCAVITYEVENTPROCESSING_HH_ +#define KTCAVITYEVENTPROCESSING_HH_ + +#include "KTProcessor.hh" +#include "KTData.hh" +#include "KTSlot.hh" + +namespace Katydid +{ + /* + @class KTCavityEventProcessing + @author J. I. Pena + @brief Assigns initial cyclotron frequency for events from symmetric-trap CCA cavity TE011 data + @details + + + Available configuration values: + (none) + + Slots: + - "mt-event": void (Nymph::KTDataPtr) -- Analyzes a multi-track-event; Requires KTMultiPeakEventData; Adds nothing + + Signals: + - "proc-cavity-event": void (Nymph::KTDataPtr) -- Emitted upon successful determination of cyclotron frequency; Guarantees KTProcessedCavityEventData + */ + + class KTMultiTrackEventData; + + class KTCavityEventProcessing : public Nymph::KTProcessor + { + public: + KTCavityEventProcessing(const std::string& name = "cavity-event-processing"); + virtual ~KTCavityEventProcessing(); + + bool Configure(const scarab::param_node* node); + + public: + bool AnalyzeEvent( KTMultiTrackEventData& mtEventData ); + + //*************** + // Signals + //*************** + + private: + Nymph::KTSignalData fProcessedCavityEventSignal; + + //*************** + // Slots + //*************** + + private: + Nymph::KTSlotDataOneType< KTMultiTrackEventData > fEventSlot; + + }; + +} + +#endif /* KTCAVITYEVENTPROCESSING_HH_ */ diff --git a/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.cc b/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.cc index d4d00992c..ce667104c 100644 --- a/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.cc +++ b/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.cc @@ -17,6 +17,7 @@ #include "KTMultiTrackEventData.hh" #include "KTPowerFitData.hh" #include "KTProcessedMPTData.hh" +#include "KTProcessedCavityEventData.hh" #include "KTProcessedTrackData.hh" #include "KTSliceHeader.hh" #include "KTSparseWaterfallCandidateData.hh" @@ -71,6 +72,7 @@ namespace Katydid fSequentialLineDataPtr(NULL), fProcessedTrackDataPtr(NULL), fProcessedMPTDataPtr(NULL), + fProcessedCavityEventData(), fMultiPeakTrackData(), fMultiTrackEventDataPtr(NULL), fMTEWithClassifierResultsDataPtr(NULL), @@ -95,6 +97,7 @@ namespace Katydid fWriter->RegisterSlot("swfc", this, &KTROOTTreeTypeWriterEventAnalysis::WriteSparseWaterfallCandidate); fWriter->RegisterSlot("seq-cand", this, &KTROOTTreeTypeWriterEventAnalysis::WriteSequentialLine); fWriter->RegisterSlot("processed-mpt", this, &KTROOTTreeTypeWriterEventAnalysis::WriteProcessedMPT); + fWriter->RegisterSlot("processed-cavity-event", this, &KTROOTTreeTypeWriterEventAnalysis::WriteProcessedCavityEvent); fWriter->RegisterSlot("proc-track", this, &KTROOTTreeTypeWriterEventAnalysis::WriteProcessedTrack); fWriter->RegisterSlot("mp-track", this, &KTROOTTreeTypeWriterEventAnalysis::WriteMultiPeakTrack); fWriter->RegisterSlot("mt-event", this, &KTROOTTreeTypeWriterEventAnalysis::WriteMultiTrackEvent); @@ -681,6 +684,76 @@ namespace Katydid return true; } + //**************************** + // Processed Cavity Event + //**************************** + + void KTROOTTreeTypeWriterEventAnalysis::WriteProcessedCavityEvent(Nymph::KTDataPtr data) + { + KTDEBUG(publog, "Attempting to write to processed cavity event root tree"); + KTProcessedCavityEventData& procCavityEventData = data->Of< KTProcessedCavityEventData >(); + + if (! fWriter->OpenAndVerifyFile()) return; + + if (fProcessedCavityEventTree == NULL) + { + if (! SetupProcessedCavityEventTree()) + { + KTERROR(publog, "Something went wrong while setting up the processed cavity mpt tree! Nothing was written."); + return; + } + } + + fProcessedCavityEventData.fComponent = procCavityEventData.GetComponent(); + fProcessedCavityEventData.fAcquisitionID = procCavityEventData.GetAcquisitionID(); + fProcessedCavityEventData.fEventID = procCavityEventData.GetEventID(); + fProcessedCavityEventData.fTotalEventSequences = procCavityEventData.GetTotalEventSequences(); + fProcessedCavityEventData.fInitialCyclotronFrequency = procCavityEventData.GetInitialCyclotronFrequency(); + + fProcessedCavityEventTree->Fill(); + + return; + } + + bool KTROOTTreeTypeWriterEventAnalysis::SetupProcessedCavityEventTree() + { + if( fWriter->GetAccumulate() ) + { + fWriter->GetFile()->GetObject( "processed-cavity-event", fProcessedCavityEventTree ); + + if (fProcessedCavityEventTree != NULL) + { + KTINFO(publog, "Tree already exists; will add to it"); + fWriter->AddTree( fProcessedCavityEventTree ); + + fProcessedCavityEventTree->SetBranchAddress( "Component", &fProcessedCavityEventData.fComponent ); + fProcessedCavityEventTree->SetBranchAddress( "AcquisitionID", &fProcessedCavityEventData.fAcquisitionID ); + fProcessedCavityEventTree->SetBranchAddress( "EventID", &fProcessedCavityEventData.fEventID ); + fProcessedCavityEventTree->SetBranchAddress( "TotalEventSequences", &fProcessedCavityEventData.fTotalEventSequences ); + fProcessedCavityEventTree->SetBranchAddress( "InitialCyclotronFrequency", &fProcessedCavityEventData.fInitialCyclotronFrequency ); + + return true; + } + } + + fProcessedCavityEventTree = new TTree("processed-cavity-event", "Processed Cavity Event"); + if( fProcessedCavityEventTree == NULL ) + { + KTERROR( publog, "Tree was not created!" ); + return false; + } + fWriter->AddTree( fProcessedCavityEventTree ); + + fProcessedCavityEventTree->Branch( "Component", &fProcessedCavityEventData.fComponent, "fComponent/i" ); + fProcessedCavityEventTree->Branch( "AcquisitionID", &fProcessedCavityEventData.fAcquisitionID, "fAcquisitionID/l" ); + fProcessedCavityEventTree->Branch( "EventID", &fProcessedCavityEventData.fEventID, "fEventID/i" ); + fProcessedCavityEventTree->Branch( "TotalEventSequences", &fProcessedCavityEventData.fTotalEventSequences, "fTotalEventSequences/i" ); + fProcessedCavityEventTree->Branch( "InitialCyclotronFrequency", &fProcessedCavityEventData.fInitialCyclotronFrequency, "fInitialCyclotronFrequency/d" ); + + return true; + } + + //************************** // Multi-Peak Track //************************** diff --git a/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh b/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh index 30348a4a5..03668b8dd 100644 --- a/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh +++ b/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh @@ -116,6 +116,15 @@ namespace Katydid UInt_t fUnknownEventTopology; }; + struct TProcessedCavityEventData + { + UInt_t fComponent; + ULong64_t fAcquisitionID; + UInt_t fEventID; + UInt_t fTotalEventSequences; + Double_t fInitialCyclotronFrequency; + }; + struct TLinearFitResult { UInt_t fFitNumber; @@ -198,6 +207,7 @@ namespace Katydid void WriteSparseWaterfallCandidate(Nymph::KTDataPtr data); void WriteSequentialLine(Nymph::KTDataPtr data); void WriteProcessedMPT(Nymph::KTDataPtr data); + void WriteProcessedCavityEvent(Nymph::KTDataPtr data); void WriteProcessedTrack(Nymph::KTDataPtr data); void WriteMultiPeakTrack(Nymph::KTDataPtr data); void WriteMultiTrackEvent(Nymph::KTDataPtr data); @@ -211,6 +221,7 @@ namespace Katydid TTree* GetSparseWaterfallCandidateTree() const; TTree* GetSequentialLineTree() const; TTree* GetProcessedMPTTree() const; + TTree* GetProcessedCavityEventTree() const; TTree* GetProcessedTrackTree() const; TTree* GetMultiPeakTrackTree() const; TTree* GetMultiTrackEventTree() const; @@ -224,6 +235,7 @@ namespace Katydid bool SetupSparseWaterfallCandidateTree(); bool SetupSequentialLineTree(); bool SetupProcessedMPTTree(); + bool SetupProcessedCavityEventTree(); bool SetupProcessedTrackTree(); bool SetupMultiPeakTrackTree(); bool SetupMultiTrackEventTree(); @@ -236,6 +248,7 @@ namespace Katydid TTree* fSparseWaterfallCandidateTree; TTree* fSequentialLineTree; TTree* fProcessedMPTTree; + TTree* fProcessedCavityEventTree; TTree* fProcessedTrackTree; TTree* fMultiPeakTrackTree; TTree* fMultiTrackEventTree; @@ -249,6 +262,7 @@ namespace Katydid TSequentialLineData* fSequentialLineDataPtr; Cicada::TProcessedTrackData* fProcessedTrackDataPtr; Cicada::TProcessedMPTData* fProcessedMPTDataPtr; + TProcessedCavityEventData fProcessedCavityEventData; TMultiPeakTrackData fMultiPeakTrackData; Cicada::TMultiTrackEventData* fMultiTrackEventDataPtr; Cicada::TMTEWithClassifierResultsData* fMTEWithClassifierResultsDataPtr; @@ -282,6 +296,11 @@ namespace Katydid return fProcessedMPTTree; } + inline TTree* KTROOTTreeTypeWriterEventAnalysis::GetProcessedCavityEventTree() const + { + return fProcessedCavityEventTree; + } + inline TTree* KTROOTTreeTypeWriterEventAnalysis::GetProcessedTrackTree() const { return fProcessedTrackTree; From 720f972f7bba7b0dbddd4a72e8c3fe390354e4c8 Mon Sep 17 00:00:00 2001 From: Junior Pena Date: Tue, 16 Sep 2025 11:14:04 -0400 Subject: [PATCH 02/16] Added code to cavity event processor that iterates through MPTs in events and sorts the bands within individual MPTs by order of start frequency. --- .../EventAnalysis/KTCavityEventProcessing.cc | 34 +++++++++++++++++++ 1 file changed, 34 insertions(+) diff --git a/Source/EventAnalysis/KTCavityEventProcessing.cc b/Source/EventAnalysis/KTCavityEventProcessing.cc index d5a734321..15bdfdd90 100644 --- a/Source/EventAnalysis/KTCavityEventProcessing.cc +++ b/Source/EventAnalysis/KTCavityEventProcessing.cc @@ -11,6 +11,8 @@ #include "KTProcessedTrackData.hh" #include "KTMultiTrackEventData.hh" +#include + namespace Katydid { @@ -41,6 +43,38 @@ namespace Katydid { KTDEBUG(evlog, "Entered Analyze Event Function"); + // initializing a map to store processed tracks in MPT groupings + std::map> mptsBySequence; + + for (TrackSetCIt trackIt = mtEventData.GetTracksBegin(); trackIt != mtEventData.GetTracksEnd(); ++trackIt) + { + const AllTrackData& track = *trackIt; + int eventSequenceID = track.fProcTrack.GetEventSequenceID(); + mptsBySequence[eventSequenceID].push_back(track); + } + + // iterating through MPTs and sorting bands by their start frequency used later for assigning track sideband order + for (auto& [eventSeqID, tracksInMPT] : mptsBySequence) + { + KTDEBUG(evlog, "MultiPeakTrack (Event Sequence ID: " << eventSeqID << ")"); + + std::multimap sortedMPTBands; + for (const auto& track : tracksInMPT) + { + sortedMPTBands.emplace(track.fProcTrack.GetStartFrequency(), &track); + } + + for (const auto& [startFreq, trackPtr] : sortedMPTBands) + { + KTDEBUG(evlog, "TrackID: " << trackPtr->fProcTrack.GetTrackID()); + KTDEBUG(evlog, "StartFrequency: " << startFreq); + KTDEBUG(evlog, "EndFrequency: " << trackPtr->fProcTrack.GetEndFrequency()); + KTDEBUG(evlog, "StartTimeInRunC: " << trackPtr->fProcTrack.GetStartTimeInRunC()); + KTDEBUG(evlog, "EndTimeInRunC: " << trackPtr->fProcTrack.GetEndTimeInRunC()); + } + + } + return true; } From 96179c3d3f546a4a84f46e01996be283b6f891be Mon Sep 17 00:00:00 2001 From: Junior Pena Date: Tue, 16 Sep 2025 17:24:35 -0400 Subject: [PATCH 03/16] Implemented cyclotron frequency reconstruction algorithm logic. Still need to go back and make 3 band track classification relative NUP threshold a configurable parameter. --- .../EventAnalysis/KTCavityEventProcessing.cc | 148 ++++++++++++++++-- 1 file changed, 131 insertions(+), 17 deletions(-) diff --git a/Source/EventAnalysis/KTCavityEventProcessing.cc b/Source/EventAnalysis/KTCavityEventProcessing.cc index 15bdfdd90..8d45a82f8 100644 --- a/Source/EventAnalysis/KTCavityEventProcessing.cc +++ b/Source/EventAnalysis/KTCavityEventProcessing.cc @@ -10,6 +10,7 @@ #include "KTLogger.hh" #include "KTProcessedTrackData.hh" #include "KTMultiTrackEventData.hh" +#include "KTProcessedCavityEventData.hh" #include @@ -41,11 +42,10 @@ namespace Katydid bool KTCavityEventProcessing::AnalyzeEvent( KTMultiTrackEventData& mtEventData ) { - KTDEBUG(evlog, "Entered Analyze Event Function"); + KTINFO(evlog, "Beginning start cyclotron frequency reconstruction of first multi-peak-tracks in events.") - // initializing a map to store processed tracks in MPT groupings + // Storing processed tracks in map by MPT groupings std::map> mptsBySequence; - for (TrackSetCIt trackIt = mtEventData.GetTracksBegin(); trackIt != mtEventData.GetTracksEnd(); ++trackIt) { const AllTrackData& track = *trackIt; @@ -53,28 +53,142 @@ namespace Katydid mptsBySequence[eventSequenceID].push_back(track); } - // iterating through MPTs and sorting bands by their start frequency used later for assigning track sideband order + // Initializing parameters used in cyclotron frequency reconstruction algorithm + double eventStartTime = mtEventData.GetStartTimeInRunC(); + double initialCyclotronFrequency = -1; + float trackClassRelPowThresh_3Bands = 0.7; + KTDEBUG(evlog, "Relative NUP threshold for classification of bands in 3 band multi-peak-tracks is " << trackClassRelPowThresh_3Bands); + + + // Iterating through MPTs for (auto& [eventSeqID, tracksInMPT] : mptsBySequence) { - KTDEBUG(evlog, "MultiPeakTrack (Event Sequence ID: " << eventSeqID << ")"); - - std::multimap sortedMPTBands; - for (const auto& track : tracksInMPT) + // Performing cyclotron reconstruction only on first MPT in event + if (eventSeqID == 0) { - sortedMPTBands.emplace(track.fProcTrack.GetStartFrequency(), &track); - } + KTDEBUG(evlog, "Looking at MultiPeakTrack with Event Sequence ID: " << eventSeqID); + + // Sorting bands in MPT by increasing order of start frequency used for track classification. + // Finding the maximum TotalTrackNUP out of all bands in MPT used for track classification of MPTs with 3 bands to decide between [-4, -2, 0] and [-2, 0, 2] topologies + double maxTotNUP = 0; + std::vector sortedMPTBands; + for (const auto& track : tracksInMPT) + { + sortedMPTBands.push_back(&track); + if (track.fProcTrack.GetTotalTrackNUP() >= maxTotNUP) {maxTotNUP=track.fProcTrack.GetTotalTrackNUP();} + } + std::sort(sortedMPTBands.begin(), sortedMPTBands.end(), + [](const AllTrackData* a, const AllTrackData* b) + { + return a->fProcTrack.GetStartFrequency() < b->fProcTrack.GetStartFrequency(); + }); + + KTDEBUG("Max Total NUP out of all bands in first track : " << maxTotNUP); + + /* + for (const AllTrackData* trackPtr : sortedMPTBands) + { + KTDEBUG(evlog, "TrackID: " << trackPtr->fProcTrack.GetTrackID()); + KTDEBUG(evlog, "TotalTrackNUP: " << trackPtr->fProcTrack.GetTotalTrackNUP()); + KTDEBUG(evlog, "StartFrequency: " << trackPtr->fProcTrack.GetStartFrequency()); + KTDEBUG(evlog, "EndFrequency: " << trackPtr->fProcTrack.GetEndFrequency()); + KTDEBUG(evlog, "StartTimeInRunC: " << trackPtr->fProcTrack.GetStartTimeInRunC()); + KTDEBUG(evlog, "EndTimeInRunC: " << trackPtr->fProcTrack.GetEndTimeInRunC()); + } + */ + + // Classifying bands in MPT. Number of bands determines possible topologies as follows + // Bands | Topologies + // 2 | [-2, 0] + // 3 | [-4, -2, 0] OR [-2, 0, 2] + // 4 | [-4, -2, 0, 2] + // Integers in topologies refer to sideband "order" (ie. -2 means second order lower sideband) + // Inital cyclotron frequency extracted from carrier with sideband order 0 by determine its position depending on topology + size_t numBands = sortedMPTBands.size(); + int carrierIndex = -1; + + if (numBands==0) + { + KTWARN(evlog, "MPT has no tracks; I cannot analyze this for cyclotron frequency. Aborting"); + return false; + } + else if (numBands==1) + { + carrierIndex = 0; + KTDEBUG(evlog, "MPT has 1 band!"); + } + else if (numBands==2) + { + carrierIndex = 1; + KTDEBUG(evlog, "MPT has 2 bands! Assumed band topology = [-2, 0].") + } + else if (numBands==3) + { + // For 3 band MPT 2 possibilities for topology. + // If the relative TotalTrackNUP of the band with the highest start frequency is larger than would be expected for a 2nd order upper sideband, it is assumed to be a carrier thus fixing the topology to [-4, -2, 0] + // Threshold for expectation of largest possible relative TotalTrackNUP for 2nd order sidebands initialized above (will be configurable parameter soon) + const AllTrackData* lastTrack = sortedMPTBands[numBands-1]; + if (maxTotNUP!=0) + { + double lastTrackRelTotalNUP = lastTrack->fProcTrack.GetTotalTrackNUP()/maxTotNUP; + if (lastTrackRelTotalNUP>trackClassRelPowThresh_3Bands) + { + carrierIndex = 2; + KTDEBUG(evlog, "MPT has 3 bands! Relative band threshold classification determined band topology = [-4, -2, 0].") + } + else + { + carrierIndex = 1; + KTDEBUG(evlog, "MPT has 3 bands! Relative band threshold classification determined band topology = [-2, 0, 2].") + + } + + } + else + { + KTWARN(evlog, "MPT has 3 bands; however, the band with the maximum total track NUP is 0 and I am unable to calculate the relative power of bands. Aborting"); + return false; + } + } + else if (numBands==4) + { + carrierIndex = 2; + KTDEBUG(evlog, "MPT has 4 bands! Assumed band topology = [-4, -2, 0, 2]."); + } + else if (numBands>=5) + { + KTWARN(evlog, "MPT has 5 or more bands; I don't know how to analyze this. Aborting"); + return false; + } + + // Calculating initial cyclotron frequency of event if correctly determined carrier location in MPT + if (carrierIndex!=-1) + { + const AllTrackData* carrier = sortedMPTBands[carrierIndex]; + initialCyclotronFrequency = carrier->fProcTrack.GetSlope()*eventStartTime + carrier->fProcTrack.GetIntercept(); + } + else + { + KTWARN(evlog, "Faulty index for carrier returned. Aborting") + return false; + } - for (const auto& [startFreq, trackPtr] : sortedMPTBands) - { - KTDEBUG(evlog, "TrackID: " << trackPtr->fProcTrack.GetTrackID()); - KTDEBUG(evlog, "StartFrequency: " << startFreq); - KTDEBUG(evlog, "EndFrequency: " << trackPtr->fProcTrack.GetEndFrequency()); - KTDEBUG(evlog, "StartTimeInRunC: " << trackPtr->fProcTrack.GetStartTimeInRunC()); - KTDEBUG(evlog, "EndTimeInRunC: " << trackPtr->fProcTrack.GetEndTimeInRunC()); } } + + // Create and fill new data object + KTProcessedCavityEventData& procEvent = mtEventData.Of(); + + procEvent.SetComponent(mtEventData.GetComponent()); + procEvent.SetAcquisitionID(mtEventData.GetAcquisitionID()); + procEvent.SetEventID(mtEventData.GetEventID()); + procEvent.SetTotalEventSequences(mtEventData.GetTotalEventSequences()); + + KTINFO(evlog, "Found initial cyclotron frequency: " << initialCyclotronFrequency) + procEvent.SetInitialCyclotronFrequency(initialCyclotronFrequency); + return true; } From 5d366ad9bc5ba1b0c67dbb0063a730ccb4166fd4 Mon Sep 17 00:00:00 2001 From: Junior Pena Date: Thu, 18 Sep 2025 15:38:09 -0400 Subject: [PATCH 04/16] Implemented logic to reconstruct the axial frequency of the first MPT. Getting break segmentation error in root writer, committing changes to see if there was random typo somewhere in writer files. --- .../KTProcessedCavityEventData.cc | 7 +- .../KTProcessedCavityEventData.hh | 1 + .../EventAnalysis/KTCavityEventProcessing.cc | 153 ++++++++++++++++-- .../EventAnalysis/KTCavityEventProcessing.hh | 3 + .../KTROOTTreeTypeWriterEventAnalysis.cc | 3 + .../KTROOTTreeTypeWriterEventAnalysis.hh | 1 + 6 files changed, 151 insertions(+), 17 deletions(-) diff --git a/Source/Data/EventAnalysis/KTProcessedCavityEventData.cc b/Source/Data/EventAnalysis/KTProcessedCavityEventData.cc index 1ff72c9f5..602fe40e6 100644 --- a/Source/Data/EventAnalysis/KTProcessedCavityEventData.cc +++ b/Source/Data/EventAnalysis/KTProcessedCavityEventData.cc @@ -20,7 +20,8 @@ namespace Katydid fAcquisitionID(0.), fEventID(0), fTotalEventSequences(0), - fInitialCyclotronFrequency(0.) + fInitialCyclotronFrequency(0.), + fFirstTrackAxialFrequency(0.) { } @@ -31,7 +32,8 @@ namespace Katydid fAcquisitionID(orig.fAcquisitionID), fEventID(orig.fEventID), fTotalEventSequences(orig.fTotalEventSequences), - fInitialCyclotronFrequency(orig.fInitialCyclotronFrequency) + fInitialCyclotronFrequency(orig.fInitialCyclotronFrequency), + fFirstTrackAxialFrequency(orig.fFirstTrackAxialFrequency) { } @@ -48,6 +50,7 @@ namespace Katydid fTotalEventSequences = rhs.fTotalEventSequences; fInitialCyclotronFrequency = rhs.fInitialCyclotronFrequency; + fFirstTrackAxialFrequency = rhs.fFirstTrackAxialFrequency; return *this; } diff --git a/Source/Data/EventAnalysis/KTProcessedCavityEventData.hh b/Source/Data/EventAnalysis/KTProcessedCavityEventData.hh index 1bcbd7c24..ae106c9c6 100644 --- a/Source/Data/EventAnalysis/KTProcessedCavityEventData.hh +++ b/Source/Data/EventAnalysis/KTProcessedCavityEventData.hh @@ -30,6 +30,7 @@ namespace Katydid MEMBERVARIABLE(unsigned, TotalEventSequences); MEMBERVARIABLE(double, InitialCyclotronFrequency); + MEMBERVARIABLE(double, FirstTrackAxialFrequency); public: static const std::string sName; diff --git a/Source/EventAnalysis/KTCavityEventProcessing.cc b/Source/EventAnalysis/KTCavityEventProcessing.cc index 8d45a82f8..033b5de61 100644 --- a/Source/EventAnalysis/KTCavityEventProcessing.cc +++ b/Source/EventAnalysis/KTCavityEventProcessing.cc @@ -13,6 +13,7 @@ #include "KTProcessedCavityEventData.hh" #include +#include namespace Katydid @@ -24,6 +25,7 @@ namespace Katydid KTCavityEventProcessing::KTCavityEventProcessing(const std::string& name) : KTProcessor(name), + fTrackClass3BandRelPowerThresh(0.), fProcessedCavityEventSignal("proc-cavity-event", this), fEventSlot("mt-event", this, &KTCavityEventProcessing::AnalyzeEvent, &fProcessedCavityEventSignal) { @@ -36,13 +38,24 @@ namespace Katydid bool KTCavityEventProcessing::Configure(const scarab::param_node* node) { if (node == NULL) return false; + + SetTrackClass3BandRelPowerThresh(node->get_value("3band-class-rel-power-thresh", GetTrackClass3BandRelPowerThresh())); return true; } bool KTCavityEventProcessing::AnalyzeEvent( KTMultiTrackEventData& mtEventData ) { - KTINFO(evlog, "Beginning start cyclotron frequency reconstruction of first multi-peak-tracks in events.") + // Create and fill new data object + KTProcessedCavityEventData& procEvent = mtEventData.Of(); + + procEvent.SetComponent(mtEventData.GetComponent()); + procEvent.SetAcquisitionID(mtEventData.GetAcquisitionID()); + procEvent.SetEventID(mtEventData.GetEventID()); + procEvent.SetTotalEventSequences(mtEventData.GetTotalEventSequences()); + procEvent.SetInitialCyclotronFrequency(-1.); + procEvent.SetFirstTrackAxialFrequency(-1.); + // Storing processed tracks in map by MPT groupings std::map> mptsBySequence; @@ -56,17 +69,18 @@ namespace Katydid // Initializing parameters used in cyclotron frequency reconstruction algorithm double eventStartTime = mtEventData.GetStartTimeInRunC(); double initialCyclotronFrequency = -1; - float trackClassRelPowThresh_3Bands = 0.7; - KTDEBUG(evlog, "Relative NUP threshold for classification of bands in 3 band multi-peak-tracks is " << trackClassRelPowThresh_3Bands); - + double firstTrackAxialFrequency = -1; // Iterating through MPTs for (auto& [eventSeqID, tracksInMPT] : mptsBySequence) { + // Performing cyclotron reconstruction only on first MPT in event if (eventSeqID == 0) { + KTINFO(evlog, "Beginning start cyclotron frequency reconstruction of first multi-peak-track in event.") KTDEBUG(evlog, "Looking at MultiPeakTrack with Event Sequence ID: " << eventSeqID); + KTDEBUG(evlog, "Relative NUP threshold for classification of bands in 3 band multi-peak-tracks is " << fTrackClass3BandRelPowerThresh); //float fTrackClass3BandRelPowerThresh = 0.7; // Tuned parameter for CCA simulation data // Sorting bands in MPT by increasing order of start frequency used for track classification. // Finding the maximum TotalTrackNUP out of all bands in MPT used for track classification of MPTs with 3 bands to decide between [-4, -2, 0] and [-2, 0, 2] topologies @@ -83,7 +97,7 @@ namespace Katydid return a->fProcTrack.GetStartFrequency() < b->fProcTrack.GetStartFrequency(); }); - KTDEBUG("Max Total NUP out of all bands in first track : " << maxTotNUP); + KTDEBUG("Max Total Track NUP out of all bands in first MPT : " << maxTotNUP); /* for (const AllTrackData* trackPtr : sortedMPTBands) @@ -131,7 +145,7 @@ namespace Katydid if (maxTotNUP!=0) { double lastTrackRelTotalNUP = lastTrack->fProcTrack.GetTotalTrackNUP()/maxTotNUP; - if (lastTrackRelTotalNUP>trackClassRelPowThresh_3Bands) + if (lastTrackRelTotalNUP>fTrackClass3BandRelPowerThresh) { carrierIndex = 2; KTDEBUG(evlog, "MPT has 3 bands! Relative band threshold classification determined band topology = [-4, -2, 0].") @@ -173,22 +187,131 @@ namespace Katydid return false; } - } - } + if(numBands>1) + { + KTINFO(evlog, "Beginning axial frequency reconstruction of first multi-peak-track in event. Attemping frequency distance calculation between bands in mpt.") + std::set timeStampSet; + std::vector frequencyDistances; + // Calculating all unique frequency distances between tracks at their start and end times + for (auto outerIt = tracksInMPT.begin(); outerIt != tracksInMPT.end(); ++outerIt) + { + double t1StartTime = outerIt->fProcTrack.GetStartTimeInRunC(); + double t1EndTime = outerIt->fProcTrack.GetEndTimeInRunC(); + double t1StartFreq = outerIt->fProcTrack.GetStartFrequency(); + double t1EndFreq = outerIt->fProcTrack.GetEndFrequency(); + double t1Slope = outerIt->fProcTrack.GetSlope(); + double t1Intercept = outerIt->fProcTrack.GetIntercept(); + + KTDEBUG(evlog, "Outer loop track id: " << outerIt->fProcTrack.GetTrackID()); - // Create and fill new data object - KTProcessedCavityEventData& procEvent = mtEventData.Of(); + // Looping through track start and end times to calculate frequency distance to other tracks at those time stamps + for (const double& timeStamp : {t1StartTime, t1EndTime}) + { + KTDEBUG(evlog, "Time stamp: " << timeStamp) + //Skip redundant time stamps + if (timeStampSet.count(timeStamp) > 0) + { + KTDEBUG(evlog, "Redundant time stamp encountered. Skipping frequency distance calculation!"); + continue; + } + //Store unique time stamps + timeStampSet.insert(timeStamp); + + double freqDist = 0; + double t1FreqAtTimeStamp = t1Slope*timeStamp + t1Intercept; + + for (auto innerIt = tracksInMPT.begin(); innerIt != tracksInMPT.end(); ++innerIt) + { + // Skip outer loop track + if (outerIt == innerIt) continue; + + double t2StartTime = innerIt->fProcTrack.GetStartTimeInRunC(); + double t2EndTime = innerIt->fProcTrack.GetEndTimeInRunC(); + double t2StartFreq = innerIt->fProcTrack.GetStartFrequency(); + double t2EndFreq = innerIt->fProcTrack.GetEndFrequency(); + double t2Slope = innerIt->fProcTrack.GetSlope(); + double t2Intercept = innerIt->fProcTrack.GetIntercept(); - procEvent.SetComponent(mtEventData.GetComponent()); - procEvent.SetAcquisitionID(mtEventData.GetAcquisitionID()); - procEvent.SetEventID(mtEventData.GetEventID()); - procEvent.SetTotalEventSequences(mtEventData.GetTotalEventSequences()); - KTINFO(evlog, "Found initial cyclotron frequency: " << initialCyclotronFrequency) + KTDEBUG(evlog, "Inner loop track id: " << innerIt->fProcTrack.GetTrackID()) + + // Check that time stamp is contained within other track lengths before calculating frequency distance + if (t2StartTime <= timeStamp && t2EndTime >= timeStamp) + { + double t2FreqAtTimeStamp = t2Slope*timeStamp + t2Intercept; + freqDist = std::abs(t1FreqAtTimeStamp - t2FreqAtTimeStamp); + frequencyDistances.push_back(freqDist); + KTDEBUG(evlog, "Successfully calculated frequency distance: " << freqDist) + } + else + { + KTDEBUG(evlog, "Time stamp not contained in inner loop track length. Aborting frequency distance calculation!") + continue; + } + } + + } + + } + + KTINFO(evlog, "Beginning calculation of average axial frequency from all frequency distances.") + // Calculating average axial frequency from frequency distances between tracks + auto minDist = std::min_element(frequencyDistances.begin(), frequencyDistances.end()); + std::vector separationOrder; + // Check if the vector is not empty and min is not zero to avoid division by zero + if (minDist != frequencyDistances.end() && *minDist != 0.0) + { + double min_val = *minDist; + // Divide all frequency distances by the minimum frequency distance and store value + for (double& val : frequencyDistances) { + separationOrder.push_back(std::round(val/min_val)); + } + } + else + { + KTDEBUG(evlog, "Error: Frequency distances vector is empty or minimum value is zero.") + } + + if (frequencyDistances.size() != separationOrder.size() || frequencyDistances.empty()) + { + KTDEBUG(evlog, "Error: Frequency distances vector empty or has different size to separation order vector.") + } + + KTDEBUG(evlog, "Separation order of frequency distances:") + for (std::size_t i = 0; i < separationOrder.size(); ++i) + { + KTDEBUG(evlog, separationOrder[i]); + } + + double sum = 0.0; + for (std::size_t i = 0; i < frequencyDistances.size(); ++i) + { + if (separationOrder[i] == 0.0) + { + KTDEBUG(evlog, "Error: Frequency distance division by 0 separation order.") + } + sum += frequencyDistances[i] / (2*separationOrder[i]);// Assuming only even order sidebands visible!!! + KTDEBUG(evlog, "Average axial frequency constribution " << i+1 << " : " << frequencyDistances[i] / (2*separationOrder[i])); + } + firstTrackAxialFrequency = sum/frequencyDistances.size(); + } + else + { + KTINFO(evlog, "First MPT in event has <= 1 band. Impossible to reconstruct axial frequency.") + } + + } + + } + + KTINFO(evlog, "Found first track initial cyclotron frequency: " << initialCyclotronFrequency) procEvent.SetInitialCyclotronFrequency(initialCyclotronFrequency); + KTINFO(evlog, "Found first track axial frequency: " << firstTrackAxialFrequency) + procEvent.SetFirstTrackAxialFrequency(firstTrackAxialFrequency); + return true; } diff --git a/Source/EventAnalysis/KTCavityEventProcessing.hh b/Source/EventAnalysis/KTCavityEventProcessing.hh index 8a92f8cf2..2759327bb 100644 --- a/Source/EventAnalysis/KTCavityEventProcessing.hh +++ b/Source/EventAnalysis/KTCavityEventProcessing.hh @@ -11,6 +11,7 @@ #include "KTProcessor.hh" #include "KTData.hh" #include "KTSlot.hh" +#include "KTMemberVariable.hh" namespace Katydid { @@ -41,6 +42,8 @@ namespace Katydid bool Configure(const scarab::param_node* node); + MEMBERVARIABLE(float, TrackClass3BandRelPowerThresh); + public: bool AnalyzeEvent( KTMultiTrackEventData& mtEventData ); diff --git a/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.cc b/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.cc index ce667104c..7a8f2e1d2 100644 --- a/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.cc +++ b/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.cc @@ -709,6 +709,7 @@ namespace Katydid fProcessedCavityEventData.fEventID = procCavityEventData.GetEventID(); fProcessedCavityEventData.fTotalEventSequences = procCavityEventData.GetTotalEventSequences(); fProcessedCavityEventData.fInitialCyclotronFrequency = procCavityEventData.GetInitialCyclotronFrequency(); + //fProcessedCavityEventData.fFirstTrackAxialFrequency = procCavityEventData.GetFirstTrackAxialFrequency(); fProcessedCavityEventTree->Fill(); @@ -731,6 +732,7 @@ namespace Katydid fProcessedCavityEventTree->SetBranchAddress( "EventID", &fProcessedCavityEventData.fEventID ); fProcessedCavityEventTree->SetBranchAddress( "TotalEventSequences", &fProcessedCavityEventData.fTotalEventSequences ); fProcessedCavityEventTree->SetBranchAddress( "InitialCyclotronFrequency", &fProcessedCavityEventData.fInitialCyclotronFrequency ); + //fProcessedCavityEventTree->SetBranchAddress( "FirstTrackAxialFrequency", &fProcessedCavityEventData.fFirstTrackAxialFrequency ); return true; } @@ -749,6 +751,7 @@ namespace Katydid fProcessedCavityEventTree->Branch( "EventID", &fProcessedCavityEventData.fEventID, "fEventID/i" ); fProcessedCavityEventTree->Branch( "TotalEventSequences", &fProcessedCavityEventData.fTotalEventSequences, "fTotalEventSequences/i" ); fProcessedCavityEventTree->Branch( "InitialCyclotronFrequency", &fProcessedCavityEventData.fInitialCyclotronFrequency, "fInitialCyclotronFrequency/d" ); + //fProcessedCavityEventTree->Branch( "FirstTrackAxialFrequency", &fProcessedCavityEventData.fFirstTrackAxialFrequency, "fFirstTrackAxialFrequency/d" ); return true; } diff --git a/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh b/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh index 03668b8dd..6ecf70b4a 100644 --- a/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh +++ b/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh @@ -123,6 +123,7 @@ namespace Katydid UInt_t fEventID; UInt_t fTotalEventSequences; Double_t fInitialCyclotronFrequency; + Double_t fFirstTrackAxialFrequency; }; struct TLinearFitResult From c8d815c0386caa35a1b545a5363344b23771a8ce Mon Sep 17 00:00:00 2001 From: Junior Pena Date: Thu, 18 Sep 2025 17:32:40 -0400 Subject: [PATCH 05/16] In previous commit I forgot to mention that I also implemented configurable parameter that is a relative power threshold to distinguish between topological classes of MPTs with 3 bands. This commit commented out writing of axial frequency to root file which caused a segmentation break. I will go back to see where this issue comes from in the next commit. --- Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.cc | 2 ++ Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.cc b/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.cc index 7a8f2e1d2..2a17535dc 100644 --- a/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.cc +++ b/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.cc @@ -711,6 +711,8 @@ namespace Katydid fProcessedCavityEventData.fInitialCyclotronFrequency = procCavityEventData.GetInitialCyclotronFrequency(); //fProcessedCavityEventData.fFirstTrackAxialFrequency = procCavityEventData.GetFirstTrackAxialFrequency(); + KTDEBUG(publog, "Writing first track initial cyclotron frequency: " << procCavityEventData.GetInitialCyclotronFrequency()); + fProcessedCavityEventTree->Fill(); return; diff --git a/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh b/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh index 6ecf70b4a..c9d8bc6c2 100644 --- a/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh +++ b/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh @@ -123,7 +123,7 @@ namespace Katydid UInt_t fEventID; UInt_t fTotalEventSequences; Double_t fInitialCyclotronFrequency; - Double_t fFirstTrackAxialFrequency; + //Double_t fFirstTrackAxialFrequency; }; struct TLinearFitResult From 2f283d58d8fc3c9e79c86827ac66717d3c91af1a Mon Sep 17 00:00:00 2001 From: Junior Pena Date: Thu, 18 Sep 2025 21:41:45 -0400 Subject: [PATCH 06/16] Writing error was a result of not having initialized ProcessedCavityEventTree member in constructor. It is working properly now after adding that line. --- .../ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.cc | 8 +++++--- .../ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh | 2 +- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.cc b/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.cc index 2a17535dc..bf5ef1d2e 100644 --- a/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.cc +++ b/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.cc @@ -63,6 +63,7 @@ namespace Katydid fProcessedTrackTree(NULL), fMultiPeakTrackTree(NULL), fMultiTrackEventTree(NULL), + fProcessedCavityEventTree(NULL), fMTEWithClassifierResultsTree(NULL), fLinearFitResultTree(NULL), fPowerFitDataTree(NULL), @@ -709,9 +710,10 @@ namespace Katydid fProcessedCavityEventData.fEventID = procCavityEventData.GetEventID(); fProcessedCavityEventData.fTotalEventSequences = procCavityEventData.GetTotalEventSequences(); fProcessedCavityEventData.fInitialCyclotronFrequency = procCavityEventData.GetInitialCyclotronFrequency(); - //fProcessedCavityEventData.fFirstTrackAxialFrequency = procCavityEventData.GetFirstTrackAxialFrequency(); + fProcessedCavityEventData.fFirstTrackAxialFrequency = procCavityEventData.GetFirstTrackAxialFrequency(); KTDEBUG(publog, "Writing first track initial cyclotron frequency: " << procCavityEventData.GetInitialCyclotronFrequency()); + KTDEBUG(publog, "Writing first track axial frequency: " << procCavityEventData.GetFirstTrackAxialFrequency()); fProcessedCavityEventTree->Fill(); @@ -734,7 +736,7 @@ namespace Katydid fProcessedCavityEventTree->SetBranchAddress( "EventID", &fProcessedCavityEventData.fEventID ); fProcessedCavityEventTree->SetBranchAddress( "TotalEventSequences", &fProcessedCavityEventData.fTotalEventSequences ); fProcessedCavityEventTree->SetBranchAddress( "InitialCyclotronFrequency", &fProcessedCavityEventData.fInitialCyclotronFrequency ); - //fProcessedCavityEventTree->SetBranchAddress( "FirstTrackAxialFrequency", &fProcessedCavityEventData.fFirstTrackAxialFrequency ); + fProcessedCavityEventTree->SetBranchAddress( "FirstTrackAxialFrequency", &fProcessedCavityEventData.fFirstTrackAxialFrequency ); return true; } @@ -753,7 +755,7 @@ namespace Katydid fProcessedCavityEventTree->Branch( "EventID", &fProcessedCavityEventData.fEventID, "fEventID/i" ); fProcessedCavityEventTree->Branch( "TotalEventSequences", &fProcessedCavityEventData.fTotalEventSequences, "fTotalEventSequences/i" ); fProcessedCavityEventTree->Branch( "InitialCyclotronFrequency", &fProcessedCavityEventData.fInitialCyclotronFrequency, "fInitialCyclotronFrequency/d" ); - //fProcessedCavityEventTree->Branch( "FirstTrackAxialFrequency", &fProcessedCavityEventData.fFirstTrackAxialFrequency, "fFirstTrackAxialFrequency/d" ); + fProcessedCavityEventTree->Branch( "FirstTrackAxialFrequency", &fProcessedCavityEventData.fFirstTrackAxialFrequency, "fFirstTrackAxialFrequency/d" ); return true; } diff --git a/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh b/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh index c9d8bc6c2..6ecf70b4a 100644 --- a/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh +++ b/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh @@ -123,7 +123,7 @@ namespace Katydid UInt_t fEventID; UInt_t fTotalEventSequences; Double_t fInitialCyclotronFrequency; - //Double_t fFirstTrackAxialFrequency; + Double_t fFirstTrackAxialFrequency; }; struct TLinearFitResult From 10f59df27f5352784acc6df79256b2232abfa6eb Mon Sep 17 00:00:00 2001 From: Junior Pena Date: Tue, 16 Dec 2025 14:09:46 -0500 Subject: [PATCH 07/16] Pulled out cycl. and ax. freq. reconstruction algorithms from inside AnalyzeEvent function. AnalyzeEvent now reconstructs frequencies for all multi-peak-tracks in an event but only saves the frequencies for the first mpt. --- .../KTProcessedCavityEventData.cc | 6 +- .../KTProcessedCavityEventData.hh | 2 +- .../EventAnalysis/KTCavityEventProcessing.cc | 455 ++++++++++-------- .../EventAnalysis/KTCavityEventProcessing.hh | 5 + .../KTROOTTreeTypeWriterEventAnalysis.cc | 8 +- .../KTROOTTreeTypeWriterEventAnalysis.hh | 2 +- 6 files changed, 259 insertions(+), 219 deletions(-) diff --git a/Source/Data/EventAnalysis/KTProcessedCavityEventData.cc b/Source/Data/EventAnalysis/KTProcessedCavityEventData.cc index 602fe40e6..88253b170 100644 --- a/Source/Data/EventAnalysis/KTProcessedCavityEventData.cc +++ b/Source/Data/EventAnalysis/KTProcessedCavityEventData.cc @@ -20,7 +20,7 @@ namespace Katydid fAcquisitionID(0.), fEventID(0), fTotalEventSequences(0), - fInitialCyclotronFrequency(0.), + fFirstTrackStartCyclotronFrequency(0.), fFirstTrackAxialFrequency(0.) { } @@ -32,7 +32,7 @@ namespace Katydid fAcquisitionID(orig.fAcquisitionID), fEventID(orig.fEventID), fTotalEventSequences(orig.fTotalEventSequences), - fInitialCyclotronFrequency(orig.fInitialCyclotronFrequency), + fFirstTrackStartCyclotronFrequency(orig.fFirstTrackStartCyclotronFrequency), fFirstTrackAxialFrequency(orig.fFirstTrackAxialFrequency) { } @@ -49,7 +49,7 @@ namespace Katydid fEventID = rhs.fEventID; fTotalEventSequences = rhs.fTotalEventSequences; - fInitialCyclotronFrequency = rhs.fInitialCyclotronFrequency; + fFirstTrackStartCyclotronFrequency = rhs.fFirstTrackStartCyclotronFrequency; fFirstTrackAxialFrequency = rhs.fFirstTrackAxialFrequency; return *this; } diff --git a/Source/Data/EventAnalysis/KTProcessedCavityEventData.hh b/Source/Data/EventAnalysis/KTProcessedCavityEventData.hh index ae106c9c6..d7f3f1b3b 100644 --- a/Source/Data/EventAnalysis/KTProcessedCavityEventData.hh +++ b/Source/Data/EventAnalysis/KTProcessedCavityEventData.hh @@ -29,7 +29,7 @@ namespace Katydid MEMBERVARIABLE(unsigned, EventID); MEMBERVARIABLE(unsigned, TotalEventSequences); - MEMBERVARIABLE(double, InitialCyclotronFrequency); + MEMBERVARIABLE(double, FirstTrackStartCyclotronFrequency); MEMBERVARIABLE(double, FirstTrackAxialFrequency); public: diff --git a/Source/EventAnalysis/KTCavityEventProcessing.cc b/Source/EventAnalysis/KTCavityEventProcessing.cc index 033b5de61..2ae9e29dd 100644 --- a/Source/EventAnalysis/KTCavityEventProcessing.cc +++ b/Source/EventAnalysis/KTCavityEventProcessing.cc @@ -53,7 +53,7 @@ namespace Katydid procEvent.SetAcquisitionID(mtEventData.GetAcquisitionID()); procEvent.SetEventID(mtEventData.GetEventID()); procEvent.SetTotalEventSequences(mtEventData.GetTotalEventSequences()); - procEvent.SetInitialCyclotronFrequency(-1.); + procEvent.SetFirstTrackStartCyclotronFrequency(-1.); procEvent.SetFirstTrackAxialFrequency(-1.); @@ -66,251 +66,286 @@ namespace Katydid mptsBySequence[eventSequenceID].push_back(track); } - // Initializing parameters used in cyclotron frequency reconstruction algorithm - double eventStartTime = mtEventData.GetStartTimeInRunC(); - double initialCyclotronFrequency = -1; - double firstTrackAxialFrequency = -1; + KTINFO(evlog, "Event Start Time : " << mtEventData.GetStartTimeInRunC()); // Iterating through MPTs for (auto& [eventSeqID, tracksInMPT] : mptsBySequence) { + KTDEBUG(evlog, "Looking at MultiPeakTrack with Event Sequence ID: " << eventSeqID); + + double reconstructedStartCyclotronFrequency = -1.0; + if (!ReconstructCyclotronFrequency(tracksInMPT, reconstructedStartCyclotronFrequency)) + { + KTWARN(evlog, "Cyclotron frequency reconstruction failed for MPT with Event Sequence ID: " << eventSeqID); + if(eventSeqID==0) {return false;} + continue; + } + + double reconstructedAxialFrequency = -1.0; + if (!ReconstructAxialFrequency(tracksInMPT, reconstructedAxialFrequency)) + { + KTWARN(evlog, "Axial frequency reconstruction failed."); + } - // Performing cyclotron reconstruction only on first MPT in event + // Storing reconstructed frequencies for only first multi-peak-track if (eventSeqID == 0) { - KTINFO(evlog, "Beginning start cyclotron frequency reconstruction of first multi-peak-track in event.") - KTDEBUG(evlog, "Looking at MultiPeakTrack with Event Sequence ID: " << eventSeqID); - KTDEBUG(evlog, "Relative NUP threshold for classification of bands in 3 band multi-peak-tracks is " << fTrackClass3BandRelPowerThresh); //float fTrackClass3BandRelPowerThresh = 0.7; // Tuned parameter for CCA simulation data + KTINFO(evlog, "Found first track initial cyclotron frequency: " << reconstructedStartCyclotronFrequency); + procEvent.SetFirstTrackStartCyclotronFrequency(reconstructedStartCyclotronFrequency); + + KTINFO(evlog, "Found first track axial frequency: " << reconstructedAxialFrequency); + procEvent.SetFirstTrackAxialFrequency(reconstructedAxialFrequency); + } + + } + + return true; + } + + bool KTCavityEventProcessing::ReconstructCyclotronFrequency(const std::vector& tracksInMPT, double& outStartCyclotronFrequency) const + { + outStartCyclotronFrequency = -1.0; + KTINFO(evlog, "Beginning start cyclotron frequency reconstruction of MPT."); + KTDEBUG(evlog, "Relative NUP threshold for classification of bands in 3 band multi-peak-tracks is " << this->fTrackClass3BandRelPowerThresh); //float fTrackClass3BandRelPowerThresh = 0.7; // Tuned parameter for CCA simulation data + + //Calculating MPT start time for start cyclotron frequency reconstruction + double startTime = tracksInMPT[0].fProcTrack.GetStartTimeInRunC(); + KTDEBUG(evlog, "Initial MPT start time guess : " << startTime); + for (const auto& track : tracksInMPT) + { + if(track.fProcTrack.GetStartTimeInRunC() < startTime) {startTime=track.fProcTrack.GetStartTimeInRunC();} - // Sorting bands in MPT by increasing order of start frequency used for track classification. - // Finding the maximum TotalTrackNUP out of all bands in MPT used for track classification of MPTs with 3 bands to decide between [-4, -2, 0] and [-2, 0, 2] topologies - double maxTotNUP = 0; - std::vector sortedMPTBands; - for (const auto& track : tracksInMPT) - { - sortedMPTBands.push_back(&track); - if (track.fProcTrack.GetTotalTrackNUP() >= maxTotNUP) {maxTotNUP=track.fProcTrack.GetTotalTrackNUP();} - } - std::sort(sortedMPTBands.begin(), sortedMPTBands.end(), - [](const AllTrackData* a, const AllTrackData* b) - { - return a->fProcTrack.GetStartFrequency() < b->fProcTrack.GetStartFrequency(); - }); + } + KTDEBUG(evlog, "MPT start time : " << startTime); - KTDEBUG("Max Total Track NUP out of all bands in first MPT : " << maxTotNUP); + // Sorting bands in MPT by increasing order of start frequency used for track classification. + // Finding the maximum TotalTrackNUP out of all bands in MPT used for track classification of MPTs with 3 bands to decide between [-4, -2, 0] and [-2, 0, 2] topologies + double maxTotNUP = 0; - /* - for (const AllTrackData* trackPtr : sortedMPTBands) - { - KTDEBUG(evlog, "TrackID: " << trackPtr->fProcTrack.GetTrackID()); - KTDEBUG(evlog, "TotalTrackNUP: " << trackPtr->fProcTrack.GetTotalTrackNUP()); - KTDEBUG(evlog, "StartFrequency: " << trackPtr->fProcTrack.GetStartFrequency()); - KTDEBUG(evlog, "EndFrequency: " << trackPtr->fProcTrack.GetEndFrequency()); - KTDEBUG(evlog, "StartTimeInRunC: " << trackPtr->fProcTrack.GetStartTimeInRunC()); - KTDEBUG(evlog, "EndTimeInRunC: " << trackPtr->fProcTrack.GetEndTimeInRunC()); - } - */ - - // Classifying bands in MPT. Number of bands determines possible topologies as follows - // Bands | Topologies - // 2 | [-2, 0] - // 3 | [-4, -2, 0] OR [-2, 0, 2] - // 4 | [-4, -2, 0, 2] - // Integers in topologies refer to sideband "order" (ie. -2 means second order lower sideband) - // Inital cyclotron frequency extracted from carrier with sideband order 0 by determine its position depending on topology - size_t numBands = sortedMPTBands.size(); - int carrierIndex = -1; - - if (numBands==0) - { - KTWARN(evlog, "MPT has no tracks; I cannot analyze this for cyclotron frequency. Aborting"); - return false; - } - else if (numBands==1) - { - carrierIndex = 0; - KTDEBUG(evlog, "MPT has 1 band!"); - } - else if (numBands==2) - { - carrierIndex = 1; - KTDEBUG(evlog, "MPT has 2 bands! Assumed band topology = [-2, 0].") - } - else if (numBands==3) - { - // For 3 band MPT 2 possibilities for topology. - // If the relative TotalTrackNUP of the band with the highest start frequency is larger than would be expected for a 2nd order upper sideband, it is assumed to be a carrier thus fixing the topology to [-4, -2, 0] - // Threshold for expectation of largest possible relative TotalTrackNUP for 2nd order sidebands initialized above (will be configurable parameter soon) - const AllTrackData* lastTrack = sortedMPTBands[numBands-1]; - if (maxTotNUP!=0) - { - double lastTrackRelTotalNUP = lastTrack->fProcTrack.GetTotalTrackNUP()/maxTotNUP; - if (lastTrackRelTotalNUP>fTrackClass3BandRelPowerThresh) - { - carrierIndex = 2; - KTDEBUG(evlog, "MPT has 3 bands! Relative band threshold classification determined band topology = [-4, -2, 0].") - } - else - { - carrierIndex = 1; - KTDEBUG(evlog, "MPT has 3 bands! Relative band threshold classification determined band topology = [-2, 0, 2].") + std::vector sortedMPTBands; + for (const auto& track : tracksInMPT) + { + sortedMPTBands.push_back(&track); + if (track.fProcTrack.GetTotalTrackNUP() >= maxTotNUP) {maxTotNUP=track.fProcTrack.GetTotalTrackNUP();} + } + std::sort(sortedMPTBands.begin(), sortedMPTBands.end(), + [](const AllTrackData* a, const AllTrackData* b) + { + return a->fProcTrack.GetStartFrequency() < b->fProcTrack.GetStartFrequency(); + }); - } - - } - else - { - KTWARN(evlog, "MPT has 3 bands; however, the band with the maximum total track NUP is 0 and I am unable to calculate the relative power of bands. Aborting"); - return false; - } - } - else if (numBands==4) - { - carrierIndex = 2; - KTDEBUG(evlog, "MPT has 4 bands! Assumed band topology = [-4, -2, 0, 2]."); - } - else if (numBands>=5) - { - KTWARN(evlog, "MPT has 5 or more bands; I don't know how to analyze this. Aborting"); - return false; - } + KTDEBUG("Max Total Track NUP out of all bands in first MPT : " << maxTotNUP); - // Calculating initial cyclotron frequency of event if correctly determined carrier location in MPT - if (carrierIndex!=-1) + /* + for (const AllTrackData* trackPtr : sortedMPTBands) + { + KTDEBUG(evlog, "TrackID: " << trackPtr->fProcTrack.GetTrackID()); + KTDEBUG(evlog, "TotalTrackNUP: " << trackPtr->fProcTrack.GetTotalTrackNUP()); + KTDEBUG(evlog, "StartFrequency: " << trackPtr->fProcTrack.GetStartFrequency()); + KTDEBUG(evlog, "EndFrequency: " << trackPtr->fProcTrack.GetEndFrequency()); + KTDEBUG(evlog, "StartTimeInRunC: " << trackPtr->fProcTrack.GetStartTimeInRunC()); + KTDEBUG(evlog, "EndTimeInRunC: " << trackPtr->fProcTrack.GetEndTimeInRunC()); + } + */ + + // Classifying bands in MPT. Number of bands determines possible topologies as follows + // Bands | Topologies + // 2 | [-2, 0] + // 3 | [-4, -2, 0] OR [-2, 0, 2] + // 4 | [-4, -2, 0, 2] + // Integers in topologies refer to sideband "order" (ie. -2 means second order lower sideband) + // Inital cyclotron frequency extracted from carrier with sideband order 0 by determine its position depending on topology + size_t numBands = sortedMPTBands.size(); + int carrierIndex = -1; + + if (numBands==0) + { + KTWARN(evlog, "MPT has no tracks; I cannot analyze this for cyclotron frequency. Aborting"); + return false; + } + else if (numBands==1) + { + carrierIndex = 0; + KTDEBUG(evlog, "MPT has 1 band!"); + } + else if (numBands==2) + { + carrierIndex = 1; + KTDEBUG(evlog, "MPT has 2 bands! Assumed band topology = [-2, 0]."); + } + else if (numBands==3) + { + // For 3 band MPT 2 possibilities for topology. + // If the relative TotalTrackNUP of the band with the highest start frequency is larger than would be expected for a 2nd order upper sideband, it is assumed to be a carrier thus fixing the topology to [-4, -2, 0] + // Threshold for expectation of largest possible relative TotalTrackNUP for 2nd order sidebands initialized above (will be configurable parameter soon) + const AllTrackData* lastTrack = sortedMPTBands[numBands-1]; + if (maxTotNUP!=0) + { + double lastTrackRelTotalNUP = lastTrack->fProcTrack.GetTotalTrackNUP()/maxTotNUP; + if (lastTrackRelTotalNUP> this->fTrackClass3BandRelPowerThresh) { - const AllTrackData* carrier = sortedMPTBands[carrierIndex]; - initialCyclotronFrequency = carrier->fProcTrack.GetSlope()*eventStartTime + carrier->fProcTrack.GetIntercept(); + carrierIndex = 2; + KTDEBUG(evlog, "MPT has 3 bands! Relative band threshold classification determined band topology = [-4, -2, 0]."); } else { - KTWARN(evlog, "Faulty index for carrier returned. Aborting") - return false; - } - + carrierIndex = 1; + KTDEBUG(evlog, "MPT has 3 bands! Relative band threshold classification determined band topology = [-2, 0, 2]."); - if(numBands>1) - { - KTINFO(evlog, "Beginning axial frequency reconstruction of first multi-peak-track in event. Attemping frequency distance calculation between bands in mpt.") - std::set timeStampSet; - std::vector frequencyDistances; + } + + } + else + { + KTWARN(evlog, "MPT has 3 bands; however, the band with the maximum total track NUP is 0 and I am unable to calculate the relative power of bands. Aborting"); + return false; + } + } + else if (numBands==4) + { + carrierIndex = 2; + KTDEBUG(evlog, "MPT has 4 bands! Assumed band topology = [-4, -2, 0, 2]."); + } + else if (numBands>=5) + { + KTWARN(evlog, "MPT has 5 or more bands; I don't know how to analyze this. Aborting"); + return false; + } - // Calculating all unique frequency distances between tracks at their start and end times - for (auto outerIt = tracksInMPT.begin(); outerIt != tracksInMPT.end(); ++outerIt) - { - double t1StartTime = outerIt->fProcTrack.GetStartTimeInRunC(); - double t1EndTime = outerIt->fProcTrack.GetEndTimeInRunC(); - double t1StartFreq = outerIt->fProcTrack.GetStartFrequency(); - double t1EndFreq = outerIt->fProcTrack.GetEndFrequency(); - double t1Slope = outerIt->fProcTrack.GetSlope(); - double t1Intercept = outerIt->fProcTrack.GetIntercept(); - - KTDEBUG(evlog, "Outer loop track id: " << outerIt->fProcTrack.GetTrackID()); + // Calculating initial cyclotron frequency of event if correctly determined carrier location in MPT + if (carrierIndex!=-1) + { + const AllTrackData* carrier = sortedMPTBands[carrierIndex]; + outStartCyclotronFrequency = carrier->fProcTrack.GetSlope()*startTime + carrier->fProcTrack.GetIntercept(); + KTDEBUG(evlog, "Found MPT start cyclotron frequency: " << outStartCyclotronFrequency); + } + else + { + KTWARN(evlog, "Faulty index for carrier returned. Aborting"); + return false; + } - // Looping through track start and end times to calculate frequency distance to other tracks at those time stamps - for (const double& timeStamp : {t1StartTime, t1EndTime}) - { - KTDEBUG(evlog, "Time stamp: " << timeStamp) - //Skip redundant time stamps - if (timeStampSet.count(timeStamp) > 0) - { - KTDEBUG(evlog, "Redundant time stamp encountered. Skipping frequency distance calculation!"); - continue; - } - //Store unique time stamps - timeStampSet.insert(timeStamp); - - double freqDist = 0; - double t1FreqAtTimeStamp = t1Slope*timeStamp + t1Intercept; - - for (auto innerIt = tracksInMPT.begin(); innerIt != tracksInMPT.end(); ++innerIt) - { - // Skip outer loop track - if (outerIt == innerIt) continue; - - double t2StartTime = innerIt->fProcTrack.GetStartTimeInRunC(); - double t2EndTime = innerIt->fProcTrack.GetEndTimeInRunC(); - double t2StartFreq = innerIt->fProcTrack.GetStartFrequency(); - double t2EndFreq = innerIt->fProcTrack.GetEndFrequency(); - double t2Slope = innerIt->fProcTrack.GetSlope(); - double t2Intercept = innerIt->fProcTrack.GetIntercept(); - - - KTDEBUG(evlog, "Inner loop track id: " << innerIt->fProcTrack.GetTrackID()) - - // Check that time stamp is contained within other track lengths before calculating frequency distance - if (t2StartTime <= timeStamp && t2EndTime >= timeStamp) - { - double t2FreqAtTimeStamp = t2Slope*timeStamp + t2Intercept; - freqDist = std::abs(t1FreqAtTimeStamp - t2FreqAtTimeStamp); - frequencyDistances.push_back(freqDist); - KTDEBUG(evlog, "Successfully calculated frequency distance: " << freqDist) - } - else - { - KTDEBUG(evlog, "Time stamp not contained in inner loop track length. Aborting frequency distance calculation!") - continue; - } - } + return true; + } - } - - } + bool KTCavityEventProcessing::ReconstructAxialFrequency(const std::vector& tracksInMPT, double& outAxialFrequency) const + { + outAxialFrequency = -1.0; + if(tracksInMPT.size()>1) + { + KTINFO(evlog, "Beginning axial frequency reconstruction of MPT."); + KTDEBUG(evlog, "Attemping frequency distance calculation between bands in MPT."); + std::set timeStampSet; + std::vector frequencyDistances; - KTINFO(evlog, "Beginning calculation of average axial frequency from all frequency distances.") - // Calculating average axial frequency from frequency distances between tracks - auto minDist = std::min_element(frequencyDistances.begin(), frequencyDistances.end()); - std::vector separationOrder; - // Check if the vector is not empty and min is not zero to avoid division by zero - if (minDist != frequencyDistances.end() && *minDist != 0.0) - { - double min_val = *minDist; - // Divide all frequency distances by the minimum frequency distance and store value - for (double& val : frequencyDistances) { - separationOrder.push_back(std::round(val/min_val)); - } - } - else + // Calculating all unique frequency distances between tracks at their start and end times + for (auto outerIt = tracksInMPT.begin(); outerIt != tracksInMPT.end(); ++outerIt) + { + double t1StartTime = outerIt->fProcTrack.GetStartTimeInRunC(); + double t1EndTime = outerIt->fProcTrack.GetEndTimeInRunC(); + double t1StartFreq = outerIt->fProcTrack.GetStartFrequency(); + double t1EndFreq = outerIt->fProcTrack.GetEndFrequency(); + double t1Slope = outerIt->fProcTrack.GetSlope(); + double t1Intercept = outerIt->fProcTrack.GetIntercept(); + + KTDEBUG(evlog, "Outer loop track id: " << outerIt->fProcTrack.GetTrackID()); + + // Looping through track start and end times to calculate frequency distance to other tracks at those time stamps + for (const double& timeStamp : {t1StartTime, t1EndTime}) + { + KTDEBUG(evlog, "Time stamp: " << timeStamp) + //Skip redundant time stamps + if (timeStampSet.count(timeStamp) > 0) { - KTDEBUG(evlog, "Error: Frequency distances vector is empty or minimum value is zero.") + KTDEBUG(evlog, "Redundant time stamp encountered. Skipping frequency distance calculation!"); + continue; } + //Store unique time stamps + timeStampSet.insert(timeStamp); - if (frequencyDistances.size() != separationOrder.size() || frequencyDistances.empty()) - { - KTDEBUG(evlog, "Error: Frequency distances vector empty or has different size to separation order vector.") - } + double freqDist = 0; + double t1FreqAtTimeStamp = t1Slope*timeStamp + t1Intercept; - KTDEBUG(evlog, "Separation order of frequency distances:") - for (std::size_t i = 0; i < separationOrder.size(); ++i) + for (auto innerIt = tracksInMPT.begin(); innerIt != tracksInMPT.end(); ++innerIt) { - KTDEBUG(evlog, separationOrder[i]); - } + // Skip outer loop track + if (outerIt == innerIt) continue; + + double t2StartTime = innerIt->fProcTrack.GetStartTimeInRunC(); + double t2EndTime = innerIt->fProcTrack.GetEndTimeInRunC(); + double t2StartFreq = innerIt->fProcTrack.GetStartFrequency(); + double t2EndFreq = innerIt->fProcTrack.GetEndFrequency(); + double t2Slope = innerIt->fProcTrack.GetSlope(); + double t2Intercept = innerIt->fProcTrack.GetIntercept(); - double sum = 0.0; - for (std::size_t i = 0; i < frequencyDistances.size(); ++i) - { - if (separationOrder[i] == 0.0) + + KTDEBUG(evlog, "Inner loop track id: " << innerIt->fProcTrack.GetTrackID()) + + // Check that time stamp is contained within other track lengths before calculating frequency distance + if (t2StartTime <= timeStamp && t2EndTime >= timeStamp) + { + double t2FreqAtTimeStamp = t2Slope*timeStamp + t2Intercept; + freqDist = std::abs(t1FreqAtTimeStamp - t2FreqAtTimeStamp); + frequencyDistances.push_back(freqDist); + KTDEBUG(evlog, "Successfully calculated frequency distance: " << freqDist) + } + else { - KTDEBUG(evlog, "Error: Frequency distance division by 0 separation order.") + KTDEBUG(evlog, "Time stamp not contained in inner loop track length. Aborting frequency distance calculation!"); + continue; } - sum += frequencyDistances[i] / (2*separationOrder[i]);// Assuming only even order sidebands visible!!! - KTDEBUG(evlog, "Average axial frequency constribution " << i+1 << " : " << frequencyDistances[i] / (2*separationOrder[i])); } - firstTrackAxialFrequency = sum/frequencyDistances.size(); - } - else - { - KTINFO(evlog, "First MPT in event has <= 1 band. Impossible to reconstruct axial frequency.") + } + + } + KTINFO(evlog, "Beginning calculation of average axial frequency from all frequency distances."); + // Calculating average axial frequency from frequency distances between tracks + auto minDist = std::min_element(frequencyDistances.begin(), frequencyDistances.end()); + std::vector separationOrder; + // Check if the vector is not empty and min is not zero to avoid division by zero + if (minDist != frequencyDistances.end() && *minDist != 0.0) + { + double min_val = *minDist; + // Divide all frequency distances by the minimum frequency distance and store value + for (double& val : frequencyDistances) { + separationOrder.push_back(std::round(val/min_val)); + } + } + else + { + KTDEBUG(evlog, "Error: Frequency distances vector is empty or minimum value is zero."); } - } + if (frequencyDistances.size() != separationOrder.size() || frequencyDistances.empty()) + { + KTDEBUG(evlog, "Error: Frequency distances vector empty or has different size to separation order vector."); + } - KTINFO(evlog, "Found first track initial cyclotron frequency: " << initialCyclotronFrequency) - procEvent.SetInitialCyclotronFrequency(initialCyclotronFrequency); + KTDEBUG(evlog, "Separation order of frequency distances:"); + for (std::size_t i = 0; i < separationOrder.size(); ++i) + { + KTDEBUG(evlog, separationOrder[i]); + } - KTINFO(evlog, "Found first track axial frequency: " << firstTrackAxialFrequency) - procEvent.SetFirstTrackAxialFrequency(firstTrackAxialFrequency); + double sum = 0.0; + for (std::size_t i = 0; i < frequencyDistances.size(); ++i) + { + if (separationOrder[i] == 0.0) + { + KTDEBUG(evlog, "Error: Frequency distance division by 0 separation order."); + } + sum += frequencyDistances[i] / (2*separationOrder[i]);// Assuming only even order sidebands visible!!! + KTDEBUG(evlog, "Average axial frequency constribution " << i+1 << " : " << frequencyDistances[i] / (2*separationOrder[i])); + } + outAxialFrequency = sum/frequencyDistances.size(); + } + else + { + KTINFO(evlog, "First MPT in event has <= 1 band. Impossible to reconstruct axial frequency."); + return false; + } return true; } diff --git a/Source/EventAnalysis/KTCavityEventProcessing.hh b/Source/EventAnalysis/KTCavityEventProcessing.hh index 2759327bb..368726dda 100644 --- a/Source/EventAnalysis/KTCavityEventProcessing.hh +++ b/Source/EventAnalysis/KTCavityEventProcessing.hh @@ -12,6 +12,7 @@ #include "KTData.hh" #include "KTSlot.hh" #include "KTMemberVariable.hh" +#include "KTProcessedTrackData.hh" namespace Katydid { @@ -46,6 +47,10 @@ namespace Katydid public: bool AnalyzeEvent( KTMultiTrackEventData& mtEventData ); + + bool ReconstructCyclotronFrequency(const std::vector& tracksInMPT, double& outStartCyclotronFrequency) const; + bool ReconstructAxialFrequency(const std::vector& tracksInMPT, double& outAxialFrequency) const; + //*************** // Signals diff --git a/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.cc b/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.cc index bf5ef1d2e..f10906eb9 100644 --- a/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.cc +++ b/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.cc @@ -709,10 +709,10 @@ namespace Katydid fProcessedCavityEventData.fAcquisitionID = procCavityEventData.GetAcquisitionID(); fProcessedCavityEventData.fEventID = procCavityEventData.GetEventID(); fProcessedCavityEventData.fTotalEventSequences = procCavityEventData.GetTotalEventSequences(); - fProcessedCavityEventData.fInitialCyclotronFrequency = procCavityEventData.GetInitialCyclotronFrequency(); + fProcessedCavityEventData.fFirstTrackStartCyclotronFrequency = procCavityEventData.GetFirstTrackStartCyclotronFrequency(); fProcessedCavityEventData.fFirstTrackAxialFrequency = procCavityEventData.GetFirstTrackAxialFrequency(); - KTDEBUG(publog, "Writing first track initial cyclotron frequency: " << procCavityEventData.GetInitialCyclotronFrequency()); + KTDEBUG(publog, "Writing first track initial cyclotron frequency: " << procCavityEventData.GetFirstTrackStartCyclotronFrequency()); KTDEBUG(publog, "Writing first track axial frequency: " << procCavityEventData.GetFirstTrackAxialFrequency()); fProcessedCavityEventTree->Fill(); @@ -735,7 +735,7 @@ namespace Katydid fProcessedCavityEventTree->SetBranchAddress( "AcquisitionID", &fProcessedCavityEventData.fAcquisitionID ); fProcessedCavityEventTree->SetBranchAddress( "EventID", &fProcessedCavityEventData.fEventID ); fProcessedCavityEventTree->SetBranchAddress( "TotalEventSequences", &fProcessedCavityEventData.fTotalEventSequences ); - fProcessedCavityEventTree->SetBranchAddress( "InitialCyclotronFrequency", &fProcessedCavityEventData.fInitialCyclotronFrequency ); + fProcessedCavityEventTree->SetBranchAddress( "FirstTrackStartCyclotronFrequency", &fProcessedCavityEventData.fFirstTrackStartCyclotronFrequency ); fProcessedCavityEventTree->SetBranchAddress( "FirstTrackAxialFrequency", &fProcessedCavityEventData.fFirstTrackAxialFrequency ); return true; @@ -754,7 +754,7 @@ namespace Katydid fProcessedCavityEventTree->Branch( "AcquisitionID", &fProcessedCavityEventData.fAcquisitionID, "fAcquisitionID/l" ); fProcessedCavityEventTree->Branch( "EventID", &fProcessedCavityEventData.fEventID, "fEventID/i" ); fProcessedCavityEventTree->Branch( "TotalEventSequences", &fProcessedCavityEventData.fTotalEventSequences, "fTotalEventSequences/i" ); - fProcessedCavityEventTree->Branch( "InitialCyclotronFrequency", &fProcessedCavityEventData.fInitialCyclotronFrequency, "fInitialCyclotronFrequency/d" ); + fProcessedCavityEventTree->Branch( "FirstTrackStartCyclotronFrequency", &fProcessedCavityEventData.fFirstTrackStartCyclotronFrequency, "fFirstTrackStartCyclotronFrequency/d" ); fProcessedCavityEventTree->Branch( "FirstTrackAxialFrequency", &fProcessedCavityEventData.fFirstTrackAxialFrequency, "fFirstTrackAxialFrequency/d" ); return true; diff --git a/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh b/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh index 6ecf70b4a..1d2c7a5a1 100644 --- a/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh +++ b/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh @@ -122,7 +122,7 @@ namespace Katydid ULong64_t fAcquisitionID; UInt_t fEventID; UInt_t fTotalEventSequences; - Double_t fInitialCyclotronFrequency; + Double_t fFirstTrackStartCyclotronFrequency; Double_t fFirstTrackAxialFrequency; }; From 6165e9932d84d727f7bba8775f97d0886c0f4a78 Mon Sep 17 00:00:00 2001 From: Junior Pena Date: Tue, 16 Dec 2025 15:32:25 -0500 Subject: [PATCH 08/16] Cleaned up guards for both reconstruction algorithms. --- .../EventAnalysis/KTCavityEventProcessing.cc | 215 +++++++++--------- 1 file changed, 112 insertions(+), 103 deletions(-) diff --git a/Source/EventAnalysis/KTCavityEventProcessing.cc b/Source/EventAnalysis/KTCavityEventProcessing.cc index 2ae9e29dd..69dbe9c73 100644 --- a/Source/EventAnalysis/KTCavityEventProcessing.cc +++ b/Source/EventAnalysis/KTCavityEventProcessing.cc @@ -105,6 +105,13 @@ namespace Katydid bool KTCavityEventProcessing::ReconstructCyclotronFrequency(const std::vector& tracksInMPT, double& outStartCyclotronFrequency) const { outStartCyclotronFrequency = -1.0; + + if (tracksInMPT.empty()) + { + KTWARN(evlog, "MPT has no tracks; cannot reconstruct cyclotron frequency."); + return false; + } + KTINFO(evlog, "Beginning start cyclotron frequency reconstruction of MPT."); KTDEBUG(evlog, "Relative NUP threshold for classification of bands in 3 band multi-peak-tracks is " << this->fTrackClass3BandRelPowerThresh); //float fTrackClass3BandRelPowerThresh = 0.7; // Tuned parameter for CCA simulation data @@ -158,12 +165,7 @@ namespace Katydid size_t numBands = sortedMPTBands.size(); int carrierIndex = -1; - if (numBands==0) - { - KTWARN(evlog, "MPT has no tracks; I cannot analyze this for cyclotron frequency. Aborting"); - return false; - } - else if (numBands==1) + if (numBands==1) { carrierIndex = 0; KTDEBUG(evlog, "MPT has 1 band!"); @@ -231,122 +233,129 @@ namespace Katydid bool KTCavityEventProcessing::ReconstructAxialFrequency(const std::vector& tracksInMPT, double& outAxialFrequency) const { outAxialFrequency = -1.0; - if(tracksInMPT.size()>1) + if (tracksInMPT.size() <= 1) { - KTINFO(evlog, "Beginning axial frequency reconstruction of MPT."); - KTDEBUG(evlog, "Attemping frequency distance calculation between bands in MPT."); - std::set timeStampSet; - std::vector frequencyDistances; + KTINFO(evlog, "MPT has <= 1 band; axial frequency reconstruction impossible."); + return false; + } + + + KTINFO(evlog, "Beginning axial frequency reconstruction of MPT."); + KTDEBUG(evlog, "Attemping frequency distance calculation between bands in MPT."); + std::set timeStampSet; + std::vector frequencyDistances; - // Calculating all unique frequency distances between tracks at their start and end times - for (auto outerIt = tracksInMPT.begin(); outerIt != tracksInMPT.end(); ++outerIt) + // Calculating all unique frequency distances between tracks at their start and end times + for (auto outerIt = tracksInMPT.begin(); outerIt != tracksInMPT.end(); ++outerIt) + { + double t1StartTime = outerIt->fProcTrack.GetStartTimeInRunC(); + double t1EndTime = outerIt->fProcTrack.GetEndTimeInRunC(); + double t1StartFreq = outerIt->fProcTrack.GetStartFrequency(); + double t1EndFreq = outerIt->fProcTrack.GetEndFrequency(); + double t1Slope = outerIt->fProcTrack.GetSlope(); + double t1Intercept = outerIt->fProcTrack.GetIntercept(); + + KTDEBUG(evlog, "Outer loop track id: " << outerIt->fProcTrack.GetTrackID()); + + // Looping through track start and end times to calculate frequency distance to other tracks at those time stamps + for (const double& timeStamp : {t1StartTime, t1EndTime}) { - double t1StartTime = outerIt->fProcTrack.GetStartTimeInRunC(); - double t1EndTime = outerIt->fProcTrack.GetEndTimeInRunC(); - double t1StartFreq = outerIt->fProcTrack.GetStartFrequency(); - double t1EndFreq = outerIt->fProcTrack.GetEndFrequency(); - double t1Slope = outerIt->fProcTrack.GetSlope(); - double t1Intercept = outerIt->fProcTrack.GetIntercept(); - - KTDEBUG(evlog, "Outer loop track id: " << outerIt->fProcTrack.GetTrackID()); + KTDEBUG(evlog, "Time stamp: " << timeStamp) + //Skip redundant time stamps + if (timeStampSet.count(timeStamp) > 0) + { + KTDEBUG(evlog, "Redundant time stamp encountered. Skipping frequency distance calculation!"); + continue; + } + //Store unique time stamps + timeStampSet.insert(timeStamp); - // Looping through track start and end times to calculate frequency distance to other tracks at those time stamps - for (const double& timeStamp : {t1StartTime, t1EndTime}) + double freqDist = 0; + double t1FreqAtTimeStamp = t1Slope*timeStamp + t1Intercept; + + for (auto innerIt = tracksInMPT.begin(); innerIt != tracksInMPT.end(); ++innerIt) { - KTDEBUG(evlog, "Time stamp: " << timeStamp) - //Skip redundant time stamps - if (timeStampSet.count(timeStamp) > 0) - { - KTDEBUG(evlog, "Redundant time stamp encountered. Skipping frequency distance calculation!"); - continue; - } - //Store unique time stamps - timeStampSet.insert(timeStamp); + // Skip outer loop track + if (outerIt == innerIt) continue; + + double t2StartTime = innerIt->fProcTrack.GetStartTimeInRunC(); + double t2EndTime = innerIt->fProcTrack.GetEndTimeInRunC(); + double t2StartFreq = innerIt->fProcTrack.GetStartFrequency(); + double t2EndFreq = innerIt->fProcTrack.GetEndFrequency(); + double t2Slope = innerIt->fProcTrack.GetSlope(); + double t2Intercept = innerIt->fProcTrack.GetIntercept(); - double freqDist = 0; - double t1FreqAtTimeStamp = t1Slope*timeStamp + t1Intercept; - for (auto innerIt = tracksInMPT.begin(); innerIt != tracksInMPT.end(); ++innerIt) + KTDEBUG(evlog, "Inner loop track id: " << innerIt->fProcTrack.GetTrackID()) + + // Check that time stamp is contained within other track lengths before calculating frequency distance + if (t2StartTime <= timeStamp && t2EndTime >= timeStamp) { - // Skip outer loop track - if (outerIt == innerIt) continue; - - double t2StartTime = innerIt->fProcTrack.GetStartTimeInRunC(); - double t2EndTime = innerIt->fProcTrack.GetEndTimeInRunC(); - double t2StartFreq = innerIt->fProcTrack.GetStartFrequency(); - double t2EndFreq = innerIt->fProcTrack.GetEndFrequency(); - double t2Slope = innerIt->fProcTrack.GetSlope(); - double t2Intercept = innerIt->fProcTrack.GetIntercept(); - - - KTDEBUG(evlog, "Inner loop track id: " << innerIt->fProcTrack.GetTrackID()) - - // Check that time stamp is contained within other track lengths before calculating frequency distance - if (t2StartTime <= timeStamp && t2EndTime >= timeStamp) - { - double t2FreqAtTimeStamp = t2Slope*timeStamp + t2Intercept; - freqDist = std::abs(t1FreqAtTimeStamp - t2FreqAtTimeStamp); - frequencyDistances.push_back(freqDist); - KTDEBUG(evlog, "Successfully calculated frequency distance: " << freqDist) - } - else - { - KTDEBUG(evlog, "Time stamp not contained in inner loop track length. Aborting frequency distance calculation!"); - continue; - } + double t2FreqAtTimeStamp = t2Slope*timeStamp + t2Intercept; + freqDist = std::abs(t1FreqAtTimeStamp - t2FreqAtTimeStamp); + frequencyDistances.push_back(freqDist); + KTDEBUG(evlog, "Successfully calculated frequency distance: " << freqDist) + } + else + { + KTDEBUG(evlog, "Time stamp not contained in inner loop track length. Aborting frequency distance calculation!"); + continue; } - - } - - } - - KTINFO(evlog, "Beginning calculation of average axial frequency from all frequency distances."); - // Calculating average axial frequency from frequency distances between tracks - auto minDist = std::min_element(frequencyDistances.begin(), frequencyDistances.end()); - std::vector separationOrder; - // Check if the vector is not empty and min is not zero to avoid division by zero - if (minDist != frequencyDistances.end() && *minDist != 0.0) - { - double min_val = *minDist; - // Divide all frequency distances by the minimum frequency distance and store value - for (double& val : frequencyDistances) { - separationOrder.push_back(std::round(val/min_val)); } - } - else - { - KTDEBUG(evlog, "Error: Frequency distances vector is empty or minimum value is zero."); - } - if (frequencyDistances.size() != separationOrder.size() || frequencyDistances.empty()) - { - KTDEBUG(evlog, "Error: Frequency distances vector empty or has different size to separation order vector."); } + + } - KTDEBUG(evlog, "Separation order of frequency distances:"); - for (std::size_t i = 0; i < separationOrder.size(); ++i) - { - KTDEBUG(evlog, separationOrder[i]); - } + if (frequencyDistances.empty()) + { + KTWARN(evlog, "No valid frequency distances found; cannot reconstruct axial frequency."); + return false; + } - double sum = 0.0; - for (std::size_t i = 0; i < frequencyDistances.size(); ++i) - { - if (separationOrder[i] == 0.0) - { - KTDEBUG(evlog, "Error: Frequency distance division by 0 separation order."); - } - sum += frequencyDistances[i] / (2*separationOrder[i]);// Assuming only even order sidebands visible!!! - KTDEBUG(evlog, "Average axial frequency constribution " << i+1 << " : " << frequencyDistances[i] / (2*separationOrder[i])); - } - outAxialFrequency = sum/frequencyDistances.size(); + KTINFO(evlog, "Beginning calculation of average axial frequency from all frequency distances."); + // Calculating average axial frequency from frequency distances between tracks + auto minDist = std::min_element(frequencyDistances.begin(), frequencyDistances.end()); + std::vector separationOrder; + // Check if the vector is not empty and min is not zero to avoid division by zero + if (minDist == frequencyDistances.end() || *minDist == 0.0) + { + KTWARN(evlog, "Cannot determine minimum frequency distance (empty or zero); cannot reconstruct axial frequency."); + return false; } - else + double min_val = *minDist; + // Divide all frequency distances by the minimum frequency distance and store value + for (double& val : frequencyDistances) { + separationOrder.push_back(std::round(val/min_val)); + } + + if (frequencyDistances.size() != separationOrder.size()) { - KTINFO(evlog, "First MPT in event has <= 1 band. Impossible to reconstruct axial frequency."); + KTWARN(evlog, "Error: Frequency distances vector has different size to separation order vector."); return false; } + KTDEBUG(evlog, "Separation order of frequency distances:"); + for (std::size_t i = 0; i < separationOrder.size(); ++i) + { + KTDEBUG(evlog, separationOrder[i]); + } + + double sum = 0.0; + for (std::size_t i = 0; i < frequencyDistances.size(); ++i) + { + if (separationOrder[i] == 0.0) + { + KTWARN(evlog, "Error: Frequency distance division by 0 separation order."); + return false; + } + sum += frequencyDistances[i] / (2*separationOrder[i]);// Assuming only even order sidebands visible!!! + KTDEBUG(evlog, "Average axial frequency constribution " << i+1 << " : " << frequencyDistances[i] / (2*separationOrder[i])); + } + + outAxialFrequency = sum/frequencyDistances.size(); + + return true; } From 4b0e15acb68329c6dbae23285dc2220c602d0e85 Mon Sep 17 00:00:00 2001 From: Junior Pena Date: Thu, 18 Dec 2025 17:02:10 -0500 Subject: [PATCH 09/16] Added vectors to KTProcessedCavityEventData to store cycl and ax frequencies for all MPTs in event. Included logic in KTCavityEventProcessing to store new data and modified writer to output the new data. --- .../KTProcessedCavityEventData.cc | 15 ++++++++++--- .../KTProcessedCavityEventData.hh | 22 +++++++++++++++++++ .../EventAnalysis/KTCavityEventProcessing.cc | 18 +++++---------- .../KTROOTTreeTypeWriterEventAnalysis.cc | 9 ++++++++ .../KTROOTTreeTypeWriterEventAnalysis.hh | 5 +++++ 5 files changed, 54 insertions(+), 15 deletions(-) diff --git a/Source/Data/EventAnalysis/KTProcessedCavityEventData.cc b/Source/Data/EventAnalysis/KTProcessedCavityEventData.cc index 88253b170..918675a8a 100644 --- a/Source/Data/EventAnalysis/KTProcessedCavityEventData.cc +++ b/Source/Data/EventAnalysis/KTProcessedCavityEventData.cc @@ -21,7 +21,11 @@ namespace Katydid fEventID(0), fTotalEventSequences(0), fFirstTrackStartCyclotronFrequency(0.), - fFirstTrackAxialFrequency(0.) + fFirstTrackAxialFrequency(0.), + fMPTEventSequenceID(), + fMPTStartCyclotronFrequency(), + fMPTAxialFrequency() + { } @@ -33,7 +37,10 @@ namespace Katydid fEventID(orig.fEventID), fTotalEventSequences(orig.fTotalEventSequences), fFirstTrackStartCyclotronFrequency(orig.fFirstTrackStartCyclotronFrequency), - fFirstTrackAxialFrequency(orig.fFirstTrackAxialFrequency) + fFirstTrackAxialFrequency(orig.fFirstTrackAxialFrequency), + fMPTEventSequenceID(orig.fMPTEventSequenceID), + fMPTStartCyclotronFrequency(orig.fMPTStartCyclotronFrequency), + fMPTAxialFrequency(orig.fMPTAxialFrequency) { } @@ -48,9 +55,11 @@ namespace Katydid fAcquisitionID = rhs.fAcquisitionID; fEventID = rhs.fEventID; fTotalEventSequences = rhs.fTotalEventSequences; - fFirstTrackStartCyclotronFrequency = rhs.fFirstTrackStartCyclotronFrequency; fFirstTrackAxialFrequency = rhs.fFirstTrackAxialFrequency; + fMPTEventSequenceID = rhs.fMPTEventSequenceID; + fMPTStartCyclotronFrequency = rhs.fMPTStartCyclotronFrequency; + fMPTAxialFrequency = rhs.fMPTAxialFrequency; return *this; } diff --git a/Source/Data/EventAnalysis/KTProcessedCavityEventData.hh b/Source/Data/EventAnalysis/KTProcessedCavityEventData.hh index d7f3f1b3b..ca29b72cf 100644 --- a/Source/Data/EventAnalysis/KTProcessedCavityEventData.hh +++ b/Source/Data/EventAnalysis/KTProcessedCavityEventData.hh @@ -12,6 +12,8 @@ #include "KTMultiTrackEventData.hh" #include "KTMemberVariable.hh" +#include + namespace Katydid { @@ -32,10 +34,30 @@ namespace Katydid MEMBERVARIABLE(double, FirstTrackStartCyclotronFrequency); MEMBERVARIABLE(double, FirstTrackAxialFrequency); + MEMBERVARIABLEREF(std::vector, MPTEventSequenceID); + MEMBERVARIABLEREF(std::vector, MPTStartCyclotronFrequency); + MEMBERVARIABLEREF(std::vector, MPTAxialFrequency); + void AddProcessedMPT(int seqID, double startCyclFreq, double axialFreq); + void ClearProcessedMPT(); + public: static const std::string sName; }; + inline void KTProcessedCavityEventData::AddProcessedMPT(int seqID, double startCyclFreq, double axialFreq) + { + fMPTEventSequenceID.push_back(seqID); + fMPTStartCyclotronFrequency.push_back(startCyclFreq); + fMPTAxialFrequency.push_back(axialFreq); + } + + inline void KTProcessedCavityEventData::ClearProcessedMPT() + { + fMPTEventSequenceID.clear(); + fMPTStartCyclotronFrequency.clear(); + fMPTAxialFrequency.clear(); + } + } #endif diff --git a/Source/EventAnalysis/KTCavityEventProcessing.cc b/Source/EventAnalysis/KTCavityEventProcessing.cc index 69dbe9c73..7f1172011 100644 --- a/Source/EventAnalysis/KTCavityEventProcessing.cc +++ b/Source/EventAnalysis/KTCavityEventProcessing.cc @@ -49,6 +49,8 @@ namespace Katydid // Create and fill new data object KTProcessedCavityEventData& procEvent = mtEventData.Of(); + procEvent.ClearProcessedMPT(); + procEvent.SetComponent(mtEventData.GetComponent()); procEvent.SetAcquisitionID(mtEventData.GetAcquisitionID()); procEvent.SetEventID(mtEventData.GetEventID()); @@ -97,6 +99,9 @@ namespace Katydid procEvent.SetFirstTrackAxialFrequency(reconstructedAxialFrequency); } + procEvent.AddProcessedMPT(eventSeqID, reconstructedStartCyclotronFrequency, reconstructedAxialFrequency); + + } return true; @@ -143,18 +148,6 @@ namespace Katydid KTDEBUG("Max Total Track NUP out of all bands in first MPT : " << maxTotNUP); - /* - for (const AllTrackData* trackPtr : sortedMPTBands) - { - KTDEBUG(evlog, "TrackID: " << trackPtr->fProcTrack.GetTrackID()); - KTDEBUG(evlog, "TotalTrackNUP: " << trackPtr->fProcTrack.GetTotalTrackNUP()); - KTDEBUG(evlog, "StartFrequency: " << trackPtr->fProcTrack.GetStartFrequency()); - KTDEBUG(evlog, "EndFrequency: " << trackPtr->fProcTrack.GetEndFrequency()); - KTDEBUG(evlog, "StartTimeInRunC: " << trackPtr->fProcTrack.GetStartTimeInRunC()); - KTDEBUG(evlog, "EndTimeInRunC: " << trackPtr->fProcTrack.GetEndTimeInRunC()); - } - */ - // Classifying bands in MPT. Number of bands determines possible topologies as follows // Bands | Topologies // 2 | [-2, 0] @@ -354,6 +347,7 @@ namespace Katydid } outAxialFrequency = sum/frequencyDistances.size(); + KTDEBUG(evlog, "Reconstructed Axial Frequency: " << outAxialFrequency); return true; diff --git a/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.cc b/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.cc index f10906eb9..f630fb788 100644 --- a/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.cc +++ b/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.cc @@ -711,6 +711,9 @@ namespace Katydid fProcessedCavityEventData.fTotalEventSequences = procCavityEventData.GetTotalEventSequences(); fProcessedCavityEventData.fFirstTrackStartCyclotronFrequency = procCavityEventData.GetFirstTrackStartCyclotronFrequency(); fProcessedCavityEventData.fFirstTrackAxialFrequency = procCavityEventData.GetFirstTrackAxialFrequency(); + fProcessedCavityEventData.fMPTEventSequenceID = procCavityEventData.GetMPTEventSequenceID(); + fProcessedCavityEventData.fMPTStartCyclotronFrequency = procCavityEventData.GetMPTStartCyclotronFrequency(); + fProcessedCavityEventData.fMPTAxialFrequency = procCavityEventData.GetMPTAxialFrequency(); KTDEBUG(publog, "Writing first track initial cyclotron frequency: " << procCavityEventData.GetFirstTrackStartCyclotronFrequency()); KTDEBUG(publog, "Writing first track axial frequency: " << procCavityEventData.GetFirstTrackAxialFrequency()); @@ -737,6 +740,9 @@ namespace Katydid fProcessedCavityEventTree->SetBranchAddress( "TotalEventSequences", &fProcessedCavityEventData.fTotalEventSequences ); fProcessedCavityEventTree->SetBranchAddress( "FirstTrackStartCyclotronFrequency", &fProcessedCavityEventData.fFirstTrackStartCyclotronFrequency ); fProcessedCavityEventTree->SetBranchAddress( "FirstTrackAxialFrequency", &fProcessedCavityEventData.fFirstTrackAxialFrequency ); + fProcessedCavityEventTree->SetBranchAddress( "MPTEventSequenceID", &fProcessedCavityEventData.fMPTEventSequenceID); + fProcessedCavityEventTree->SetBranchAddress( "MPTStartCyclotronFrequency", &fProcessedCavityEventData.fMPTStartCyclotronFrequency); + fProcessedCavityEventTree->SetBranchAddress( "MPTAxialFrequency", &fProcessedCavityEventData.fMPTAxialFrequency); return true; } @@ -756,6 +762,9 @@ namespace Katydid fProcessedCavityEventTree->Branch( "TotalEventSequences", &fProcessedCavityEventData.fTotalEventSequences, "fTotalEventSequences/i" ); fProcessedCavityEventTree->Branch( "FirstTrackStartCyclotronFrequency", &fProcessedCavityEventData.fFirstTrackStartCyclotronFrequency, "fFirstTrackStartCyclotronFrequency/d" ); fProcessedCavityEventTree->Branch( "FirstTrackAxialFrequency", &fProcessedCavityEventData.fFirstTrackAxialFrequency, "fFirstTrackAxialFrequency/d" ); + fProcessedCavityEventTree->Branch( "MPTEventSequenceID", "std::vector", &fProcessedCavityEventData.fMPTEventSequenceID); + fProcessedCavityEventTree->Branch( "MPTStartCyclotronFrequency", "std::vector", &fProcessedCavityEventData.fMPTStartCyclotronFrequency); + fProcessedCavityEventTree->Branch( "MPTAxialFrequency", "std::vector", &fProcessedCavityEventData.fMPTAxialFrequency); return true; } diff --git a/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh b/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh index 1d2c7a5a1..2fd7327ab 100644 --- a/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh +++ b/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh @@ -124,6 +124,11 @@ namespace Katydid UInt_t fTotalEventSequences; Double_t fFirstTrackStartCyclotronFrequency; Double_t fFirstTrackAxialFrequency; + + std::vector fMPTEventSequenceID; + std::vector fMPTStartCyclotronFrequency; + std::vector fMPTAxialFrequency; + }; struct TLinearFitResult From 41fa10b5560b0de084358f0fe4dfab806ca536e1 Mon Sep 17 00:00:00 2001 From: Junior Pena Date: Tue, 13 Jan 2026 16:51:34 -0500 Subject: [PATCH 10/16] Added data storage of individual track (ids, event seq ids, and sideband classifications) for users to visualize classification after reconstruction and tune 3 band classification parameter accordingly. The new data is stored as flatten vectors, as that was the easiest implementation. --- .../KTProcessedCavityEventData.cc | 13 ++++- .../KTProcessedCavityEventData.hh | 28 ++++++++++ .../EventAnalysis/KTCavityEventProcessing.cc | 56 +++++++++++++------ .../EventAnalysis/KTCavityEventProcessing.hh | 3 +- .../KTROOTTreeTypeWriterEventAnalysis.cc | 9 +++ .../KTROOTTreeTypeWriterEventAnalysis.hh | 4 ++ 6 files changed, 94 insertions(+), 19 deletions(-) diff --git a/Source/Data/EventAnalysis/KTProcessedCavityEventData.cc b/Source/Data/EventAnalysis/KTProcessedCavityEventData.cc index 918675a8a..f24b0bdd1 100644 --- a/Source/Data/EventAnalysis/KTProcessedCavityEventData.cc +++ b/Source/Data/EventAnalysis/KTProcessedCavityEventData.cc @@ -24,7 +24,10 @@ namespace Katydid fFirstTrackAxialFrequency(0.), fMPTEventSequenceID(), fMPTStartCyclotronFrequency(), - fMPTAxialFrequency() + fMPTAxialFrequency(), + fAllTrackIDs(), + fAllTrackEventSequenceIDs(), + fAllTrackBandClassifications() { } @@ -40,7 +43,10 @@ namespace Katydid fFirstTrackAxialFrequency(orig.fFirstTrackAxialFrequency), fMPTEventSequenceID(orig.fMPTEventSequenceID), fMPTStartCyclotronFrequency(orig.fMPTStartCyclotronFrequency), - fMPTAxialFrequency(orig.fMPTAxialFrequency) + fMPTAxialFrequency(orig.fMPTAxialFrequency), + fAllTrackIDs(orig.fAllTrackIDs), + fAllTrackEventSequenceIDs(orig.fAllTrackEventSequenceIDs), + fAllTrackBandClassifications(orig.fAllTrackBandClassifications) { } @@ -60,6 +66,9 @@ namespace Katydid fMPTEventSequenceID = rhs.fMPTEventSequenceID; fMPTStartCyclotronFrequency = rhs.fMPTStartCyclotronFrequency; fMPTAxialFrequency = rhs.fMPTAxialFrequency; + fAllTrackIDs = rhs.fAllTrackIDs; + fAllTrackEventSequenceIDs = rhs.fAllTrackEventSequenceIDs; + fAllTrackBandClassifications = rhs.fAllTrackBandClassifications; return *this; } diff --git a/Source/Data/EventAnalysis/KTProcessedCavityEventData.hh b/Source/Data/EventAnalysis/KTProcessedCavityEventData.hh index ca29b72cf..a2c2138ad 100644 --- a/Source/Data/EventAnalysis/KTProcessedCavityEventData.hh +++ b/Source/Data/EventAnalysis/KTProcessedCavityEventData.hh @@ -13,6 +13,7 @@ #include "KTMemberVariable.hh" #include +#include namespace Katydid { @@ -37,9 +38,17 @@ namespace Katydid MEMBERVARIABLEREF(std::vector, MPTEventSequenceID); MEMBERVARIABLEREF(std::vector, MPTStartCyclotronFrequency); MEMBERVARIABLEREF(std::vector, MPTAxialFrequency); + + MEMBERVARIABLEREF(std::vector, AllTrackIDs); + MEMBERVARIABLEREF(std::vector, AllTrackEventSequenceIDs); + MEMBERVARIABLEREF(std::vector, AllTrackBandClassifications); + void AddProcessedMPT(int seqID, double startCyclFreq, double axialFreq); void ClearProcessedMPT(); + void AddClassificationData(int trackID, int seqID, int bandClassification); + void AddClassificationData(const std::vector& tracksInMPT, int seqID, int bandClassification); + public: static const std::string sName; }; @@ -56,6 +65,25 @@ namespace Katydid fMPTEventSequenceID.clear(); fMPTStartCyclotronFrequency.clear(); fMPTAxialFrequency.clear(); + + fAllTrackIDs.clear(); + fAllTrackEventSequenceIDs.clear(); + fAllTrackBandClassifications.clear(); + } + + inline void KTProcessedCavityEventData::AddClassificationData(int trackID, int seqID, int bandClassification) + { + fAllTrackIDs.push_back(trackID); + fAllTrackEventSequenceIDs.push_back(seqID); + fAllTrackBandClassifications.push_back(bandClassification); + } + + inline void KTProcessedCavityEventData::AddClassificationData(const std::vector& tracksInMPT, int seqID, int bandClassification) + { + for (const auto& trk : tracksInMPT) + { + AddClassificationData(trk.fProcTrack.GetTrackID(), seqID, bandClassification); + } } } diff --git a/Source/EventAnalysis/KTCavityEventProcessing.cc b/Source/EventAnalysis/KTCavityEventProcessing.cc index 7f1172011..d663e7833 100644 --- a/Source/EventAnalysis/KTCavityEventProcessing.cc +++ b/Source/EventAnalysis/KTCavityEventProcessing.cc @@ -73,10 +73,12 @@ namespace Katydid // Iterating through MPTs for (auto& [eventSeqID, tracksInMPT] : mptsBySequence) { + + KTDEBUG(evlog, "Looking at MultiPeakTrack with Event Sequence ID: " << eventSeqID); double reconstructedStartCyclotronFrequency = -1.0; - if (!ReconstructCyclotronFrequency(tracksInMPT, reconstructedStartCyclotronFrequency)) + if (!ReconstructCyclotronFrequency(tracksInMPT, procEvent, reconstructedStartCyclotronFrequency)) { KTWARN(evlog, "Cyclotron frequency reconstruction failed for MPT with Event Sequence ID: " << eventSeqID); if(eventSeqID==0) {return false;} @@ -107,7 +109,7 @@ namespace Katydid return true; } - bool KTCavityEventProcessing::ReconstructCyclotronFrequency(const std::vector& tracksInMPT, double& outStartCyclotronFrequency) const + bool KTCavityEventProcessing::ReconstructCyclotronFrequency(const std::vector& tracksInMPT, KTProcessedCavityEventData& procEvent, double& outStartCyclotronFrequency) const { outStartCyclotronFrequency = -1.0; @@ -120,25 +122,21 @@ namespace Katydid KTINFO(evlog, "Beginning start cyclotron frequency reconstruction of MPT."); KTDEBUG(evlog, "Relative NUP threshold for classification of bands in 3 band multi-peak-tracks is " << this->fTrackClass3BandRelPowerThresh); //float fTrackClass3BandRelPowerThresh = 0.7; // Tuned parameter for CCA simulation data - //Calculating MPT start time for start cyclotron frequency reconstruction - double startTime = tracksInMPT[0].fProcTrack.GetStartTimeInRunC(); - KTDEBUG(evlog, "Initial MPT start time guess : " << startTime); - for (const auto& track : tracksInMPT) - { - if(track.fProcTrack.GetStartTimeInRunC() < startTime) {startTime=track.fProcTrack.GetStartTimeInRunC();} - - } - KTDEBUG(evlog, "MPT start time : " << startTime); - // Sorting bands in MPT by increasing order of start frequency used for track classification. // Finding the maximum TotalTrackNUP out of all bands in MPT used for track classification of MPTs with 3 bands to decide between [-4, -2, 0] and [-2, 0, 2] topologies + // Calculating MPT start time for start cyclotron frequency reconstruction double maxTotNUP = 0; - + double startTime = tracksInMPT.front().fProcTrack.GetStartTimeInRunC(); + KTDEBUG(evlog, "Initial MPT start time guess : " << startTime); std::vector sortedMPTBands; + size_t numBands = tracksInMPT.size(); + sortedMPTBands.reserve(numBands); + for (const auto& track : tracksInMPT) { sortedMPTBands.push_back(&track); - if (track.fProcTrack.GetTotalTrackNUP() >= maxTotNUP) {maxTotNUP=track.fProcTrack.GetTotalTrackNUP();} + if (track.fProcTrack.GetTotalTrackNUP() > maxTotNUP) {maxTotNUP=track.fProcTrack.GetTotalTrackNUP();} + if (track.fProcTrack.GetStartTimeInRunC() < startTime) {startTime=track.fProcTrack.GetStartTimeInRunC();} } std::sort(sortedMPTBands.begin(), sortedMPTBands.end(), [](const AllTrackData* a, const AllTrackData* b) @@ -146,7 +144,8 @@ namespace Katydid return a->fProcTrack.GetStartFrequency() < b->fProcTrack.GetStartFrequency(); }); - KTDEBUG("Max Total Track NUP out of all bands in first MPT : " << maxTotNUP); + KTDEBUG(evlog, "Max Total Track NUP out of all bands in MPT : " << maxTotNUP); + KTDEBUG(evlog, "MPT start time : " << startTime); // Classifying bands in MPT. Number of bands determines possible topologies as follows // Bands | Topologies @@ -155,17 +154,19 @@ namespace Katydid // 4 | [-4, -2, 0, 2] // Integers in topologies refer to sideband "order" (ie. -2 means second order lower sideband) // Inital cyclotron frequency extracted from carrier with sideband order 0 by determine its position depending on topology - size_t numBands = sortedMPTBands.size(); int carrierIndex = -1; + std::vector bandClassification; if (numBands==1) { carrierIndex = 0; + bandClassification = {0}; KTDEBUG(evlog, "MPT has 1 band!"); } else if (numBands==2) { carrierIndex = 1; + bandClassification = {-2, 0}; KTDEBUG(evlog, "MPT has 2 bands! Assumed band topology = [-2, 0]."); } else if (numBands==3) @@ -180,11 +181,13 @@ namespace Katydid if (lastTrackRelTotalNUP> this->fTrackClass3BandRelPowerThresh) { carrierIndex = 2; + bandClassification = {-4, -2, 0}; KTDEBUG(evlog, "MPT has 3 bands! Relative band threshold classification determined band topology = [-4, -2, 0]."); } else { carrierIndex = 1; + bandClassification = {-2, 0, 2}; KTDEBUG(evlog, "MPT has 3 bands! Relative band threshold classification determined band topology = [-2, 0, 2]."); } @@ -199,6 +202,7 @@ namespace Katydid else if (numBands==4) { carrierIndex = 2; + bandClassification = {-4, -2, 0, 2}; KTDEBUG(evlog, "MPT has 4 bands! Assumed band topology = [-4, -2, 0, 2]."); } else if (numBands>=5) @@ -220,6 +224,26 @@ namespace Katydid return false; } + if (bandClassification.size() != numBands) + { + KTWARN(evlog, "Band classification size (" << bandClassification.size()<< ") is different from number of bands found in MPT (" << numBands << "). Aborting"); + return false; + } + + const int seqID = sortedMPTBands.front()->fProcTrack.GetEventSequenceID(); + KTDEBUG(evlog, "Adding individual track sideband classification: (trackID, classification)"); + for (std::size_t i = 0; i < numBands; ++i) + { + const auto* track = sortedMPTBands[i]; + if (track->fProcTrack.GetEventSequenceID() != seqID) + { + KTWARN(evlog, "Sorted bands contain multiple EventSequenceIDs(expected " << seqID << "), something upstream is wrong. Aborting"); + return false; + } + procEvent.AddClassificationData(track->fProcTrack.GetTrackID(), seqID, bandClassification[i]); + KTDEBUG(evlog, "( " << track->fProcTrack.GetTrackID() << ", " << bandClassification[i] << " )"); + } + return true; } diff --git a/Source/EventAnalysis/KTCavityEventProcessing.hh b/Source/EventAnalysis/KTCavityEventProcessing.hh index 368726dda..0832b5e98 100644 --- a/Source/EventAnalysis/KTCavityEventProcessing.hh +++ b/Source/EventAnalysis/KTCavityEventProcessing.hh @@ -13,6 +13,7 @@ #include "KTSlot.hh" #include "KTMemberVariable.hh" #include "KTProcessedTrackData.hh" +#include "KTProcessedCavityEventData.hh" namespace Katydid { @@ -48,7 +49,7 @@ namespace Katydid public: bool AnalyzeEvent( KTMultiTrackEventData& mtEventData ); - bool ReconstructCyclotronFrequency(const std::vector& tracksInMPT, double& outStartCyclotronFrequency) const; + bool ReconstructCyclotronFrequency(const std::vector& tracksInMPT, KTProcessedCavityEventData& procEvent, double& outStartCyclotronFrequency) const; bool ReconstructAxialFrequency(const std::vector& tracksInMPT, double& outAxialFrequency) const; diff --git a/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.cc b/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.cc index f630fb788..a624ea028 100644 --- a/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.cc +++ b/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.cc @@ -714,6 +714,9 @@ namespace Katydid fProcessedCavityEventData.fMPTEventSequenceID = procCavityEventData.GetMPTEventSequenceID(); fProcessedCavityEventData.fMPTStartCyclotronFrequency = procCavityEventData.GetMPTStartCyclotronFrequency(); fProcessedCavityEventData.fMPTAxialFrequency = procCavityEventData.GetMPTAxialFrequency(); + fProcessedCavityEventData.fAllTrackIDs = procCavityEventData.GetAllTrackIDs(); + fProcessedCavityEventData.fAllTrackEventSequenceIDs = procCavityEventData.GetAllTrackEventSequenceIDs(); + fProcessedCavityEventData.fAllTrackBandClassifications = procCavityEventData.GetAllTrackBandClassifications(); KTDEBUG(publog, "Writing first track initial cyclotron frequency: " << procCavityEventData.GetFirstTrackStartCyclotronFrequency()); KTDEBUG(publog, "Writing first track axial frequency: " << procCavityEventData.GetFirstTrackAxialFrequency()); @@ -743,6 +746,9 @@ namespace Katydid fProcessedCavityEventTree->SetBranchAddress( "MPTEventSequenceID", &fProcessedCavityEventData.fMPTEventSequenceID); fProcessedCavityEventTree->SetBranchAddress( "MPTStartCyclotronFrequency", &fProcessedCavityEventData.fMPTStartCyclotronFrequency); fProcessedCavityEventTree->SetBranchAddress( "MPTAxialFrequency", &fProcessedCavityEventData.fMPTAxialFrequency); + fProcessedCavityEventTree->SetBranchAddress( "AllTrackIDs", &fProcessedCavityEventData.fAllTrackIDs); + fProcessedCavityEventTree->SetBranchAddress( "AllTrackEventSequenceIDs", &fProcessedCavityEventData.fAllTrackEventSequenceIDs); + fProcessedCavityEventTree->SetBranchAddress( "AllTrackBandClassifications", &fProcessedCavityEventData.fAllTrackBandClassifications); return true; } @@ -765,6 +771,9 @@ namespace Katydid fProcessedCavityEventTree->Branch( "MPTEventSequenceID", "std::vector", &fProcessedCavityEventData.fMPTEventSequenceID); fProcessedCavityEventTree->Branch( "MPTStartCyclotronFrequency", "std::vector", &fProcessedCavityEventData.fMPTStartCyclotronFrequency); fProcessedCavityEventTree->Branch( "MPTAxialFrequency", "std::vector", &fProcessedCavityEventData.fMPTAxialFrequency); + fProcessedCavityEventTree->Branch( "AllTrackIDs", "std::vector", &fProcessedCavityEventData.fAllTrackIDs); + fProcessedCavityEventTree->Branch( "AllTrackEventSequenceIDs", "std::vector", &fProcessedCavityEventData.fAllTrackEventSequenceIDs); + fProcessedCavityEventTree->Branch( "AllTrackBandClassifications", "std::vector", &fProcessedCavityEventData.fAllTrackBandClassifications); return true; } diff --git a/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh b/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh index 2fd7327ab..6fc5f148d 100644 --- a/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh +++ b/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh @@ -129,6 +129,10 @@ namespace Katydid std::vector fMPTStartCyclotronFrequency; std::vector fMPTAxialFrequency; + std::vector fAllTrackIDs; + std::vector fAllTrackEventSequenceIDs; + std::vector fAllTrackBandClassifications; + }; struct TLinearFitResult From 9eff712d352e93ee1f75850529bb7de224a6a8e5 Mon Sep 17 00:00:00 2001 From: Junior Pena Date: Fri, 30 Jan 2026 13:54:45 -0500 Subject: [PATCH 11/16] Added band classification data storage, so users do not have to extract data when only interested in first MPT of event. --- .../KTProcessedCavityEventData.cc | 6 +++ .../KTProcessedCavityEventData.hh | 37 +++++++++++++++---- .../EventAnalysis/KTCavityEventProcessing.cc | 23 ++++++++++-- .../KTROOTTreeTypeWriterEventAnalysis.cc | 6 +++ .../KTROOTTreeTypeWriterEventAnalysis.hh | 2 + 5 files changed, 63 insertions(+), 11 deletions(-) diff --git a/Source/Data/EventAnalysis/KTProcessedCavityEventData.cc b/Source/Data/EventAnalysis/KTProcessedCavityEventData.cc index f24b0bdd1..56ec664cf 100644 --- a/Source/Data/EventAnalysis/KTProcessedCavityEventData.cc +++ b/Source/Data/EventAnalysis/KTProcessedCavityEventData.cc @@ -22,6 +22,8 @@ namespace Katydid fTotalEventSequences(0), fFirstTrackStartCyclotronFrequency(0.), fFirstTrackAxialFrequency(0.), + fFirstTrackBandIDs(), + fFirstTrackBandClassifications(), fMPTEventSequenceID(), fMPTStartCyclotronFrequency(), fMPTAxialFrequency(), @@ -41,6 +43,8 @@ namespace Katydid fTotalEventSequences(orig.fTotalEventSequences), fFirstTrackStartCyclotronFrequency(orig.fFirstTrackStartCyclotronFrequency), fFirstTrackAxialFrequency(orig.fFirstTrackAxialFrequency), + fFirstTrackBandIDs(orig.fFirstTrackBandIDs), + fFirstTrackBandClassifications(orig.fFirstTrackBandClassifications), fMPTEventSequenceID(orig.fMPTEventSequenceID), fMPTStartCyclotronFrequency(orig.fMPTStartCyclotronFrequency), fMPTAxialFrequency(orig.fMPTAxialFrequency), @@ -63,6 +67,8 @@ namespace Katydid fTotalEventSequences = rhs.fTotalEventSequences; fFirstTrackStartCyclotronFrequency = rhs.fFirstTrackStartCyclotronFrequency; fFirstTrackAxialFrequency = rhs.fFirstTrackAxialFrequency; + fFirstTrackBandIDs = rhs.fFirstTrackBandIDs; + fFirstTrackBandClassifications = rhs.fFirstTrackBandClassifications; fMPTEventSequenceID = rhs.fMPTEventSequenceID; fMPTStartCyclotronFrequency = rhs.fMPTStartCyclotronFrequency; fMPTAxialFrequency = rhs.fMPTAxialFrequency; diff --git a/Source/Data/EventAnalysis/KTProcessedCavityEventData.hh b/Source/Data/EventAnalysis/KTProcessedCavityEventData.hh index a2c2138ad..eb3889d41 100644 --- a/Source/Data/EventAnalysis/KTProcessedCavityEventData.hh +++ b/Source/Data/EventAnalysis/KTProcessedCavityEventData.hh @@ -34,6 +34,8 @@ namespace Katydid MEMBERVARIABLE(double, FirstTrackStartCyclotronFrequency); MEMBERVARIABLE(double, FirstTrackAxialFrequency); + MEMBERVARIABLE(std::vector, FirstTrackBandIDs); + MEMBERVARIABLE(std::vector, FirstTrackBandClassifications); MEMBERVARIABLEREF(std::vector, MPTEventSequenceID); MEMBERVARIABLEREF(std::vector, MPTStartCyclotronFrequency); @@ -43,8 +45,12 @@ namespace Katydid MEMBERVARIABLEREF(std::vector, AllTrackEventSequenceIDs); MEMBERVARIABLEREF(std::vector, AllTrackBandClassifications); + void ClearProcessedEvent(); + + void AddFirstTrackBandInfo(int trackID, int classification); + void ClearFirstTrackBandInfo(); + void AddProcessedMPT(int seqID, double startCyclFreq, double axialFreq); - void ClearProcessedMPT(); void AddClassificationData(int trackID, int seqID, int bandClassification); void AddClassificationData(const std::vector& tracksInMPT, int seqID, int bandClassification); @@ -53,15 +59,11 @@ namespace Katydid static const std::string sName; }; - inline void KTProcessedCavityEventData::AddProcessedMPT(int seqID, double startCyclFreq, double axialFreq) + inline void KTProcessedCavityEventData::ClearProcessedEvent() { - fMPTEventSequenceID.push_back(seqID); - fMPTStartCyclotronFrequency.push_back(startCyclFreq); - fMPTAxialFrequency.push_back(axialFreq); - } + fFirstTrackBandIDs.clear(); + fFirstTrackBandClassifications.clear(); - inline void KTProcessedCavityEventData::ClearProcessedMPT() - { fMPTEventSequenceID.clear(); fMPTStartCyclotronFrequency.clear(); fMPTAxialFrequency.clear(); @@ -71,6 +73,25 @@ namespace Katydid fAllTrackBandClassifications.clear(); } + inline void KTProcessedCavityEventData::AddFirstTrackBandInfo(int trackID, int classification) + { + fFirstTrackBandIDs.push_back(trackID); + fFirstTrackBandClassifications.push_back(classification); + } + + inline void KTProcessedCavityEventData::ClearFirstTrackBandInfo() + { + fFirstTrackBandIDs.clear(); + fFirstTrackBandClassifications.clear(); + } + + inline void KTProcessedCavityEventData::AddProcessedMPT(int seqID, double startCyclFreq, double axialFreq) + { + fMPTEventSequenceID.push_back(seqID); + fMPTStartCyclotronFrequency.push_back(startCyclFreq); + fMPTAxialFrequency.push_back(axialFreq); + } + inline void KTProcessedCavityEventData::AddClassificationData(int trackID, int seqID, int bandClassification) { fAllTrackIDs.push_back(trackID); diff --git a/Source/EventAnalysis/KTCavityEventProcessing.cc b/Source/EventAnalysis/KTCavityEventProcessing.cc index d663e7833..73dd8cb6a 100644 --- a/Source/EventAnalysis/KTCavityEventProcessing.cc +++ b/Source/EventAnalysis/KTCavityEventProcessing.cc @@ -49,7 +49,7 @@ namespace Katydid // Create and fill new data object KTProcessedCavityEventData& procEvent = mtEventData.Of(); - procEvent.ClearProcessedMPT(); + procEvent.ClearProcessedEvent(); procEvent.SetComponent(mtEventData.GetComponent()); procEvent.SetAcquisitionID(mtEventData.GetAcquisitionID()); @@ -230,8 +230,16 @@ namespace Katydid return false; } + const int seqID = sortedMPTBands.front()->fProcTrack.GetEventSequenceID(); KTDEBUG(evlog, "Adding individual track sideband classification: (trackID, classification)"); + + // Clearing values for first MPT band classification + if (seqID==0) + { + procEvent.ClearFirstTrackBandInfo(); + } + for (std::size_t i = 0; i < numBands; ++i) { const auto* track = sortedMPTBands[i]; @@ -240,8 +248,17 @@ namespace Katydid KTWARN(evlog, "Sorted bands contain multiple EventSequenceIDs(expected " << seqID << "), something upstream is wrong. Aborting"); return false; } - procEvent.AddClassificationData(track->fProcTrack.GetTrackID(), seqID, bandClassification[i]); - KTDEBUG(evlog, "( " << track->fProcTrack.GetTrackID() << ", " << bandClassification[i] << " )"); + const int trackID = track->fProcTrack.GetTrackID(); + const int classification = bandClassification[i]; + + procEvent.AddClassificationData(trackID, seqID, classification); + KTDEBUG(evlog, "( " << trackID << ", " << classification << " )"); + + if (seqID == 0) + { + procEvent.AddFirstTrackBandInfo(trackID, classification); + KTINFO(evlog, "First MPT Band Classifications (trackID, bandClassification): (" << trackID << ", " << classification << " )"); + } } return true; diff --git a/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.cc b/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.cc index a624ea028..afb80471b 100644 --- a/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.cc +++ b/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.cc @@ -711,6 +711,8 @@ namespace Katydid fProcessedCavityEventData.fTotalEventSequences = procCavityEventData.GetTotalEventSequences(); fProcessedCavityEventData.fFirstTrackStartCyclotronFrequency = procCavityEventData.GetFirstTrackStartCyclotronFrequency(); fProcessedCavityEventData.fFirstTrackAxialFrequency = procCavityEventData.GetFirstTrackAxialFrequency(); + fProcessedCavityEventData.fFirstTrackBandIDs = procCavityEventData.GetFirstTrackBandIDs(); + fProcessedCavityEventData.fFirstTrackBandClassifications = procCavityEventData.GetFirstTrackBandClassifications(); fProcessedCavityEventData.fMPTEventSequenceID = procCavityEventData.GetMPTEventSequenceID(); fProcessedCavityEventData.fMPTStartCyclotronFrequency = procCavityEventData.GetMPTStartCyclotronFrequency(); fProcessedCavityEventData.fMPTAxialFrequency = procCavityEventData.GetMPTAxialFrequency(); @@ -743,6 +745,8 @@ namespace Katydid fProcessedCavityEventTree->SetBranchAddress( "TotalEventSequences", &fProcessedCavityEventData.fTotalEventSequences ); fProcessedCavityEventTree->SetBranchAddress( "FirstTrackStartCyclotronFrequency", &fProcessedCavityEventData.fFirstTrackStartCyclotronFrequency ); fProcessedCavityEventTree->SetBranchAddress( "FirstTrackAxialFrequency", &fProcessedCavityEventData.fFirstTrackAxialFrequency ); + fProcessedCavityEventTree->SetBranchAddress( "FirstTrackBandIDs", &fProcessedCavityEventData.fFirstTrackBandIDs ); + fProcessedCavityEventTree->SetBranchAddress( "FirstTrackBandClassifications", &fProcessedCavityEventData.fFirstTrackBandClassifications ); fProcessedCavityEventTree->SetBranchAddress( "MPTEventSequenceID", &fProcessedCavityEventData.fMPTEventSequenceID); fProcessedCavityEventTree->SetBranchAddress( "MPTStartCyclotronFrequency", &fProcessedCavityEventData.fMPTStartCyclotronFrequency); fProcessedCavityEventTree->SetBranchAddress( "MPTAxialFrequency", &fProcessedCavityEventData.fMPTAxialFrequency); @@ -768,6 +772,8 @@ namespace Katydid fProcessedCavityEventTree->Branch( "TotalEventSequences", &fProcessedCavityEventData.fTotalEventSequences, "fTotalEventSequences/i" ); fProcessedCavityEventTree->Branch( "FirstTrackStartCyclotronFrequency", &fProcessedCavityEventData.fFirstTrackStartCyclotronFrequency, "fFirstTrackStartCyclotronFrequency/d" ); fProcessedCavityEventTree->Branch( "FirstTrackAxialFrequency", &fProcessedCavityEventData.fFirstTrackAxialFrequency, "fFirstTrackAxialFrequency/d" ); + fProcessedCavityEventTree->Branch( "FirstTrackBandIDs", "std::vector", &fProcessedCavityEventData.fFirstTrackBandIDs); + fProcessedCavityEventTree->Branch( "FirstTrackBandClassifications", "std::vector", &fProcessedCavityEventData.fFirstTrackBandClassifications); fProcessedCavityEventTree->Branch( "MPTEventSequenceID", "std::vector", &fProcessedCavityEventData.fMPTEventSequenceID); fProcessedCavityEventTree->Branch( "MPTStartCyclotronFrequency", "std::vector", &fProcessedCavityEventData.fMPTStartCyclotronFrequency); fProcessedCavityEventTree->Branch( "MPTAxialFrequency", "std::vector", &fProcessedCavityEventData.fMPTAxialFrequency); diff --git a/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh b/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh index 6fc5f148d..b61aac6b4 100644 --- a/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh +++ b/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh @@ -124,6 +124,8 @@ namespace Katydid UInt_t fTotalEventSequences; Double_t fFirstTrackStartCyclotronFrequency; Double_t fFirstTrackAxialFrequency; + std::vector fFirstTrackBandIDs; + std::vector fFirstTrackBandClassifications; std::vector fMPTEventSequenceID; std::vector fMPTStartCyclotronFrequency; From 8f18af6a56f50c6edc777580e762e778252c0525 Mon Sep 17 00:00:00 2001 From: Junior Pena Date: Fri, 30 Jan 2026 14:08:53 -0500 Subject: [PATCH 12/16] Added comment on configurable parameter in KTCavityEventProcessing. --- Source/EventAnalysis/KTCavityEventProcessing.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/EventAnalysis/KTCavityEventProcessing.hh b/Source/EventAnalysis/KTCavityEventProcessing.hh index 0832b5e98..df8dcb516 100644 --- a/Source/EventAnalysis/KTCavityEventProcessing.hh +++ b/Source/EventAnalysis/KTCavityEventProcessing.hh @@ -25,7 +25,7 @@ namespace Katydid Available configuration values: - (none) + - "3band-class-rel-power-thresh": float -- Relative power threshold which is criteria in classifying bands in MPTs with 3 bands into either [-4,-2,0] or [-2,0,2] Slots: - "mt-event": void (Nymph::KTDataPtr) -- Analyzes a multi-track-event; Requires KTMultiPeakEventData; Adds nothing From f44a39662051d373c32de1459f6a0ccd699e219f Mon Sep 17 00:00:00 2001 From: Junior Pena Date: Mon, 4 May 2026 14:56:07 -0400 Subject: [PATCH 13/16] Updated variable and parameter name for configurable parameter in KTCavityEventProcessing. Also updated logic in 2 band MPT case to correctly classify 2 band MPTs where 2nd order lower sideband was missed. --- .../EventAnalysis/KTCavityEventProcessing.cc | 43 ++++++++++++++----- .../EventAnalysis/KTCavityEventProcessing.hh | 4 +- 2 files changed, 35 insertions(+), 12 deletions(-) diff --git a/Source/EventAnalysis/KTCavityEventProcessing.cc b/Source/EventAnalysis/KTCavityEventProcessing.cc index 73dd8cb6a..337974f67 100644 --- a/Source/EventAnalysis/KTCavityEventProcessing.cc +++ b/Source/EventAnalysis/KTCavityEventProcessing.cc @@ -25,7 +25,7 @@ namespace Katydid KTCavityEventProcessing::KTCavityEventProcessing(const std::string& name) : KTProcessor(name), - fTrackClass3BandRelPowerThresh(0.), + fMaxRelPowerThresh(0.), fProcessedCavityEventSignal("proc-cavity-event", this), fEventSlot("mt-event", this, &KTCavityEventProcessing::AnalyzeEvent, &fProcessedCavityEventSignal) { @@ -39,7 +39,7 @@ namespace Katydid { if (node == NULL) return false; - SetTrackClass3BandRelPowerThresh(node->get_value("3band-class-rel-power-thresh", GetTrackClass3BandRelPowerThresh())); + SetMaxRelPowerThresh(node->get_value("max-rel-power-thresh", GetMaxRelPowerThresh())); return true; } @@ -120,7 +120,7 @@ namespace Katydid } KTINFO(evlog, "Beginning start cyclotron frequency reconstruction of MPT."); - KTDEBUG(evlog, "Relative NUP threshold for classification of bands in 3 band multi-peak-tracks is " << this->fTrackClass3BandRelPowerThresh); //float fTrackClass3BandRelPowerThresh = 0.7; // Tuned parameter for CCA simulation data + KTDEBUG(evlog, "Relative NUP threshold for classification of bands in 3 band multi-peak-tracks is " << this->fMaxRelPowerThresh); //float fMaxRelPowerThresh = 0.1; // Tuned parameter for CCA simulation data // Sorting bands in MPT by increasing order of start frequency used for track classification. // Finding the maximum TotalTrackNUP out of all bands in MPT used for track classification of MPTs with 3 bands to decide between [-4, -2, 0] and [-2, 0, 2] topologies @@ -148,8 +148,8 @@ namespace Katydid KTDEBUG(evlog, "MPT start time : " << startTime); // Classifying bands in MPT. Number of bands determines possible topologies as follows - // Bands | Topologies - // 2 | [-2, 0] + // Bands | Topologies + // 2 | [-2, 0] OR [0, 2] // 3 | [-4, -2, 0] OR [-2, 0, 2] // 4 | [-4, -2, 0, 2] // Integers in topologies refer to sideband "order" (ie. -2 means second order lower sideband) @@ -165,20 +165,43 @@ namespace Katydid } else if (numBands==2) { - carrierIndex = 1; - bandClassification = {-2, 0}; - KTDEBUG(evlog, "MPT has 2 bands! Assumed band topology = [-2, 0]."); + // For 2 band MPT 2 possibilities for topology. + // If the relative TotalTrackNUP of the band with the highest start frequency is smaller than would be expected for a carrier, assume 2nd order lower sideband was missed and assign [0, 2] + // Threshold for expectation of largest possible relative TotalTrackNUP is the configurable parameter + const AllTrackData* lastTrack = sortedMPTBands[numBands-1]; + if (maxTotNUP!=0) + { + double lastTrackRelTotalNUP = lastTrack->fProcTrack.GetTotalTrackNUP()/maxTotNUP; + if (lastTrackRelTotalNUP> this->fMaxRelPowerThresh) + { + carrierIndex = 1; + bandClassification = {-2, 0}; + KTDEBUG(evlog, "MPT has 2 bands! Relative band threshold classification determined band topology = [-2, 0]."); + } + else + { + carrierIndex = 0; + bandClassification = {0, 2}; + KTDEBUG(evlog, "MPT has 2 bands! Relative band threshold classification determined band topology = [0, 2]."); + + } + + } + else + { + KTWARN(evlog, "MPT has 2 bands; however, the band with the maximum total track NUP is 0 and I am unable to calculate the relative power of bands. Aborting"); + return false; + } } else if (numBands==3) { // For 3 band MPT 2 possibilities for topology. // If the relative TotalTrackNUP of the band with the highest start frequency is larger than would be expected for a 2nd order upper sideband, it is assumed to be a carrier thus fixing the topology to [-4, -2, 0] - // Threshold for expectation of largest possible relative TotalTrackNUP for 2nd order sidebands initialized above (will be configurable parameter soon) const AllTrackData* lastTrack = sortedMPTBands[numBands-1]; if (maxTotNUP!=0) { double lastTrackRelTotalNUP = lastTrack->fProcTrack.GetTotalTrackNUP()/maxTotNUP; - if (lastTrackRelTotalNUP> this->fTrackClass3BandRelPowerThresh) + if (lastTrackRelTotalNUP> this->fMaxRelPowerThresh) { carrierIndex = 2; bandClassification = {-4, -2, 0}; diff --git a/Source/EventAnalysis/KTCavityEventProcessing.hh b/Source/EventAnalysis/KTCavityEventProcessing.hh index df8dcb516..601eb900f 100644 --- a/Source/EventAnalysis/KTCavityEventProcessing.hh +++ b/Source/EventAnalysis/KTCavityEventProcessing.hh @@ -25,7 +25,7 @@ namespace Katydid Available configuration values: - - "3band-class-rel-power-thresh": float -- Relative power threshold which is criteria in classifying bands in MPTs with 3 bands into either [-4,-2,0] or [-2,0,2] + - "max-rel-power-thresh": float -- Relative power threshold which is criteria in classifying bands in MPTs Slots: - "mt-event": void (Nymph::KTDataPtr) -- Analyzes a multi-track-event; Requires KTMultiPeakEventData; Adds nothing @@ -44,7 +44,7 @@ namespace Katydid bool Configure(const scarab::param_node* node); - MEMBERVARIABLE(float, TrackClass3BandRelPowerThresh); + MEMBERVARIABLE(float, MaxRelPowerThresh); public: bool AnalyzeEvent( KTMultiTrackEventData& mtEventData ); From 44ca300e61e9084a7bfb77b800b8551ed2711cff Mon Sep 17 00:00:00 2001 From: Junior Pena Date: Mon, 18 May 2026 16:04:17 -0400 Subject: [PATCH 14/16] Added example Katydid config file used on Locust-Kass simulated noiseless CCA data to validate the axial and cyclotron frequency reconstruction algorithms. --- .../KatydidNoiselessCCAFreqReco.yaml | 416 ++++++++++++++++++ 1 file changed, 416 insertions(+) create mode 100644 Examples/ConfigFiles/KatydidNoiselessCCAFreqReco.yaml diff --git a/Examples/ConfigFiles/KatydidNoiselessCCAFreqReco.yaml b/Examples/ConfigFiles/KatydidNoiselessCCAFreqReco.yaml new file mode 100644 index 000000000..b0f195c07 --- /dev/null +++ b/Examples/ConfigFiles/KatydidNoiselessCCAFreqReco.yaml @@ -0,0 +1,416 @@ +processor-toolbox: + + processors: + + - type: egg-processor + name: egg1 + - type: gaussian-noise-generator + name: noise-gen1 + - type: forward-fftw + name: fft1 + - type: convert-to-power + name: to-ps1 + + - type: egg-processor + name: egg2 + - type: gaussian-noise-generator + name: noise-gen2 + - type: forward-fftw + name: fft2 + - type: convert-to-power + name: to-ps2 + + - type: data-accumulator + name: acc + - type: gain-variation + name: gainvar + + - type: variable-spectrum-discriminator + name: discrim + - type: sequential-track-finder + name: seq-tr-clust + - type: overlapping-track-clustering + name: otc + - type: iterative-track-clustering + name: itc + - type: track-proc-ws + name: tr-proc-ws + - type: multi-peak-track-builder + name: mptb + - type: multi-peak-event-builder + name: mpeb + - type: cavity-event-processing + name: cavity-event-proc + + - type: apply-cut + name: ac1b + - type: apply-cut + name: event-nup-cut + - type: apply-cut + name: event-time-cut + + - type: basic-root-writer + name: brw + - type: root-tree-writer + name: trw + + connections: + + # Header processing + - signal: "egg1:header" + slot: "fft1:header" + - signal: "egg1:ts" + slot: "noise-gen1:slice" + - signal: "noise-gen1:slice" + slot: "fft1:ts-fftw" + + # First egg processing + - signal: "fft1:fft" + slot: "to-ps1:fs-fftw-to-psd" + - signal: "to-ps1:psd" + slot: "acc:ps" + - signal: "acc:ps-finished" + slot: "gainvar:ps-var" + # - signal: "acc:ps-finished" + # slot: "brw:psd" # Write accumulated PSD histogram + # - signal: "gainvar:gain-var" + # slot: "brw:gain-var" # Write gain-var histogram + - signal: "gainvar:gain-var" + slot: "discrim:gv" + + # Second egg processing + - signal: "egg2:header" + slot: "fft2:header" + - signal: "egg2:ts" + slot: "noise-gen2:slice" + - signal: "noise-gen2:slice" + slot: "fft2:ts-fftw" + - signal: "fft2:fft" + slot: "to-ps2:fs-fftw-to-psd" + - signal: "to-ps2:psd" + slot: "discrim:ps-pre" + - signal: "discrim:disc-1d" + slot: "seq-tr-clust:disc-1d" +# - signal: "discrim:disc-1d" +# slot: "trw:disc-1d" # Writing discriminated points to ROOT tree + - signal: "seq-tr-clust:seq-cand" + slot: "otc:seq-cand" + - signal: "otc:seq-cand" + slot: "itc:seq-cand" + - signal: "itc:seq-cand" + slot: "ac1b:apply" + - signal: "ac1b:pass" + slot: "tr-proc-ws:seq-cand" + - signal: "tr-proc-ws:track" + slot: "mptb:track" + - signal: "tr-proc-ws:track" + slot: "trw:proc-track" # Write processed tracks to a ROOT tree + - signal: "egg2:egg-done" + slot: "seq-tr-clust:done" + - signal: "seq-tr-clust:clustering-done" + slot: "otc:do-clustering" + - signal: "otc:clustering-done" + slot: "itc:do-clustering" + - signal: "itc:clustering-done" + slot: "mptb:do-clustering" + - signal: "mptb:mpt" + slot: "mpeb:mpt" + - signal: "mptb:mpt-done" + slot: "mpeb:do-clustering" + - signal: "mpeb:event" + slot: "event-nup-cut:apply" + - signal: "event-nup-cut:pass" + slot: "event-time-cut:apply" + - signal: "event-time-cut:pass" + slot: "trw:mt-event" + - signal: "event-time-cut:pass" + slot: "cavity-event-proc:mt-event" + - signal: "cavity-event-proc:proc-cavity-event" + slot: "trw:processed-cavity-event" + + + + run-queue: + - egg1 + - egg2 + + +egg1: + filename: "noElectronSignalFile.egg" + egg-reader: egg3 + number-of-slices: 0 + slice-size: 8192 + +noise-gen1: + mean: 0.0 + noise-floor-psd: 2.2e-13 # Noise power in W/Hz + seed: 12345 + +fft1: + transform-flag: ESTIMATE + +acc: + number-to-average: 0 + signal-interval: 0 + +gainvar: + normalize: false + min-frequency: 0 + max-frequency: 403e6 + fit-points: 50 + +egg2: + filename: "$eggfilename" + egg-reader: egg3 + number-of-slices: 0 + slice-size: 8192 + +noise-gen2: + mean: 0.0 + noise-floor-psd: 2.2e-13 # Noise power in W/Hz + seed: 12345 + +fft2: + transform-flag: ESTIMATE + +discrim: + min-frequency: 105.e6 # Here we chop the frequency range to (100 +/- + max-frequency: 168.e6 # 42.5) MHz, the active bandwidth of the RSA + #snr-threshold-power: 6.0 + sigma-threshold: 9.0 + normalize: true + neighborhood-radius: 2 + + +seq-tr-clust: + min-frequency: 0e6 + max-frequency: 403e6 + initial-slope : 0.1e6 # slope of a track with 1 point + slope-method: weighted-first-point-ref + n-slope-points: 10 + time-gap-tolerance: 0.1e-3 + frequency-acceptance: 56e3 + trimming-threshold: 6.0 # snr threshold for the first and last point of a candidate + min-points: 3 + min-slope: 0 + +otc: + max-track-width: 150e3 # For most tracks this only needs to be 150e3-ish. Only curved tracks profit from this number being bigger + +itc: + time-gap-tolerance: 1.0e-3 + frequency-acceptance: 100.0e3 + max-track-width: 100.0e3 + +tr-proc-ws: + min-slope: 0 + +mptb: + sideband-time-tol: 1e-3 + +mpeb: + jump-time-tol: 3.0e-3 + +cavity-event-proc: + max-rel-power-thresh: 0.1 + +ac1b: + seq-line-nup-cut: # can be applied after stf, after otc or after itc + min-total-nup: 0 + min-average-nup: 7.0 + time-or-bin-average: bin + wide-or-narrow: wide # wide is default + + +event-nup-cut: + ntracks-npoints-nup-cut: # cut events if the summed or average nup is below threshold + wide-or-narrow: wide # wide is default + time-or-bin-average: bin + + default-parameters: # default-parameters are thresholds for event first track nup cut + min-total-nup: 0 + min-average-nup: 7.454357052489394 + min-max-nup: 0 + + parameters: + - ft-npoints: 3 + ntracks: 1 + min-total-nup: 0 + min-average-nup: 15.74668644848628 + min-max-nup: 0 + + - ft-npoints: 3 + ntracks: 2 + min-total-nup: 0 + min-average-nup: 14.396403491431705 + min-max-nup: 0 + + - ft-npoints: 3 + ntracks: 3 + min-total-nup: 0 + min-average-nup: 13.032932650345844 + min-max-nup: 0 + + - ft-npoints: 3 + ntracks: 4 + min-total-nup: 0 + min-average-nup: 11.651179906121854 + min-max-nup: 0 + + - ft-npoints: 3 + ntracks: 5 + min-total-nup: 0 + min-average-nup: 10.242207772240628 + min-max-nup: 0 + + - ft-npoints: 3 + ntracks: 6 + min-total-nup: 0 + min-average-nup: 8.78782196327628 + min-max-nup: 0 + + - ft-npoints: 4 + ntracks: 1 + min-total-nup: 0 + min-average-nup: 12.735418234003228 + min-max-nup: 0 + + - ft-npoints: 4 + ntracks: 2 + min-total-nup: 0 + min-average-nup: 11.681651640175646 + min-max-nup: 0 + + - ft-npoints: 4 + ntracks: 3 + min-total-nup: 0 + min-average-nup: 10.610759932504259 + min-max-nup: 0 + + - ft-npoints: 4 + ntracks: 4 + min-total-nup: 0 + min-average-nup: 9.515056968514685 + min-max-nup: 0 + + - ft-npoints: 4 + ntracks: 5 + min-total-nup: 0 + min-average-nup: 8.379691045166915 + min-max-nup: 0 + + - ft-npoints: 5 + ntracks: 1 + min-total-nup: 0 + min-average-nup: 10.869025304065605 + min-max-nup: 0 + + - ft-npoints: 5 + ntracks: 2 + min-total-nup: 0 + min-average-nup: 9.987041877860296 + min-max-nup: 0 + + - ft-npoints: 5 + ntracks: 3 + min-total-nup: 0 + min-average-nup: 9.083012737896947 + min-max-nup: 0 + + - ft-npoints: 5 + ntracks: 4 + min-total-nup: 0 + min-average-nup: 8.144983909889586 + min-max-nup: 0 + + - ft-npoints: 6 + ntracks: 1 + min-total-nup: 0 + min-average-nup: 9.831510394306793 + min-max-nup: 0 + + - ft-npoints: 6 + ntracks: 2 + min-total-nup: 0 + min-average-nup: 9.148730213081482 + min-max-nup: 0 + + - ft-npoints: 6 + ntracks: 3 + min-total-nup: 0 + min-average-nup: 8.436488780835134 + min-max-nup: 0 + + - ft-npoints: 6 + ntracks: 4 + min-total-nup: 0 + min-average-nup: 7.67794591793375 + min-max-nup: 0 + + - ft-npoints: 7 + ntracks: 1 + min-total-nup: 0 + min-average-nup: 8.893879615160312 + min-max-nup: 0 + + - ft-npoints: 7 + ntracks: 2 + min-total-nup: 0 + min-average-nup: 8.294285562892895 + min-max-nup: 0 + + - ft-npoints: 7 + ntracks: 3 + min-total-nup: 0 + min-average-nup: 7.658861173418163 + min-max-nup: 0 + + - ft-npoints: 8 + ntracks: 1 + min-total-nup: 0 + min-average-nup: 8.277213651025455 + min-max-nup: 0 + + - ft-npoints: 8 + ntracks: 2 + min-total-nup: 0 + min-average-nup: 7.722303890626092 + min-max-nup: 0 + + - ft-npoints: 9 + ntracks: 1 + min-total-nup: 0 + min-average-nup: 7.9962037692806796 + min-max-nup: 0 + + - ft-npoints: 9 + ntracks: 2 + min-total-nup: 0 + min-average-nup: 7.46339335340889 + min-max-nup: 0 + + - ft-npoints: 10 + ntracks: 1 + min-total-nup: 0 + min-average-nup: 7.725369787485035 + min-max-nup: 0 + + - ft-npoints: 11 + ntracks: 1 + min-total-nup: 0 + min-average-nup: 7.454357052489394 + min-max-nup: 0 + +event-time-cut: + event-time-in-acq-cut: + min-time: 0.0e-3 + max-time: 1.0e-1 + +brw: + output-file: "$katydidbasicfilename" + file-flag: recreate + +trw: + output-file: "$katydidoutputfilename" + file-flag: recreate + From 0dcee67b1b7c0f1905063e5afe6e34fcc502f3f6 Mon Sep 17 00:00:00 2001 From: Junior Pena Date: Fri, 22 May 2026 20:40:02 -0400 Subject: [PATCH 15/16] Added a more precise description of the KTCavityEventProcessing processor in brief and details. --- Source/EventAnalysis/KTCavityEventProcessing.hh | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/Source/EventAnalysis/KTCavityEventProcessing.hh b/Source/EventAnalysis/KTCavityEventProcessing.hh index 601eb900f..d7a307ccc 100644 --- a/Source/EventAnalysis/KTCavityEventProcessing.hh +++ b/Source/EventAnalysis/KTCavityEventProcessing.hh @@ -20,9 +20,14 @@ namespace Katydid /* @class KTCavityEventProcessing @author J. I. Pena - @brief Assigns initial cyclotron frequency for events from symmetric-trap CCA cavity TE011 data - @details + @brief Reconstructs axial and cyclotron frequency for all multi-peak-tracks in an event + @details + Iterates through all MPTs in an event by sequence ID + Reconstructs axial frequency by calculating frequency distance between different bands in MPTs at distinct band start and end times + Assumes an axially symmetric trap and cavity mode map such that sideband separation is 2*n*f_ax + Reconstructs cyclotron frequency by classifying MPT topology based on relative power of bands to the band in the MPT with the highest power + Relative power threshold parameter used 2 decide between 2 classes in the cases of 2 and 3 band MPTs Available configuration values: - "max-rel-power-thresh": float -- Relative power threshold which is criteria in classifying bands in MPTs From 5bc23a6c37699b93c5b51670d0312e7c00fc0d86 Mon Sep 17 00:00:00 2001 From: Junior Pena Date: Fri, 22 May 2026 20:52:32 -0400 Subject: [PATCH 16/16] Corrected a typo in KTCavityEventProcessing.hh comment --- Source/EventAnalysis/KTCavityEventProcessing.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/EventAnalysis/KTCavityEventProcessing.hh b/Source/EventAnalysis/KTCavityEventProcessing.hh index d7a307ccc..a56aa663b 100644 --- a/Source/EventAnalysis/KTCavityEventProcessing.hh +++ b/Source/EventAnalysis/KTCavityEventProcessing.hh @@ -33,7 +33,7 @@ namespace Katydid - "max-rel-power-thresh": float -- Relative power threshold which is criteria in classifying bands in MPTs Slots: - - "mt-event": void (Nymph::KTDataPtr) -- Analyzes a multi-track-event; Requires KTMultiPeakEventData; Adds nothing + - "mt-event": void (Nymph::KTDataPtr) -- Analyzes a multi-track-event; Requires KTMultiTrackEventData; Adds nothing Signals: - "proc-cavity-event": void (Nymph::KTDataPtr) -- Emitted upon successful determination of cyclotron frequency; Guarantees KTProcessedCavityEventData