Feature/axial freq recon

Source/Data/CMakeLists.txt

@@ -30,6 +30,7 @@ set (DATA_HEADERFILES
     EventAnalysis/KTMultiTrackEventData.hh
     EventAnalysis/KTPowerFitData.hh
     EventAnalysis/KTProcessedMPTData.hh
+    EventAnalysis/KTProcessedCavityMPTData.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/KTProcessedCavityMPTData.cc
     EventAnalysis/KTProcessedTrackData.cc
     EventAnalysis/KTRPTrackData.cc
     EventAnalysis/KTSequentialLineData.cc

Source/Data/EventAnalysis/KTProcessedCavityMPTData.cc

@@ -0,0 +1,44 @@
+/*
+ * KTProcessedCavityMPTData.cc
+ *
+ *  Created on: Aug 19, 2025
+ *      Author: juniorpe
+ */
+
+#include "KTProcessedCavityMPTData.hh"
+
+#include "KTLogger.hh"
+
+namespace Katydid
+{
+
+    const std::string KTProcessedCavityMPTData::sName("proc-cavity-mpt");
+
+    KTProcessedCavityMPTData::KTProcessedCavityMPTData() :
+            KTExtensibleData< KTProcessedCavityMPTData >(),
+            fComponent(0),
+            fAxialFrequency(0.)
+    {
+    }
+
+    KTProcessedCavityMPTData::KTProcessedCavityMPTData(const KTProcessedCavityMPTData& orig) :
+            KTExtensibleData< KTProcessedCavityMPTData >(orig),
+
+            fComponent(orig.fComponent),
+            fAxialFrequency(orig.fAxialFrequency)
+    {
+    }
+
+    KTProcessedCavityMPTData::~KTProcessedCavityMPTData()
+    {
+    }
+
+    KTProcessedCavityMPTData& KTProcessedCavityMPTData::operator=(const KTProcessedCavityMPTData& rhs)
+    {
+        KTExtensibleData< KTProcessedCavityMPTData >::operator=(rhs);
+        fComponent = rhs.fComponent;
+        fAxialFrequency = rhs.fAxialFrequency;
+        return *this;
+    }
+
+}

Source/Data/EventAnalysis/KTProcessedCavityMPTData.hh

@@ -0,0 +1,36 @@
+/*
+ * KTProcessedCavityMPTData.hh
+ *
+ *  Created on: Aug 19, 2025
+ *      Author: juniorpe
+ */
+
+#ifndef KTPROCESSEDCAVITYMPTDATA_HH_
+#define KTPROCESSEDCAVITYMPTDATA_HH_
+
+#include "KTData.hh"
+#include "KTProcessedTrackData.hh"
+#include "KTMemberVariable.hh"
+
+namespace Katydid
+{
+
+    class KTProcessedCavityMPTData : public Nymph::KTExtensibleData< KTProcessedCavityMPTData >
+    {
+        public:
+            KTProcessedCavityMPTData();
+            KTProcessedCavityMPTData(const KTProcessedCavityMPTData& orig);
+            virtual ~KTProcessedCavityMPTData();
+
+            KTProcessedCavityMPTData& operator=(const KTProcessedCavityMPTData& rhs);
+
+            MEMBERVARIABLE(unsigned, Component);
+            MEMBERVARIABLE(double, AxialFrequency);
+
+        public:
+            static const std::string sName;
+    };
+
+}
+#endif
+

Source/EventAnalysis/CMakeLists.txt

@@ -34,6 +34,7 @@ set (EVENTANALYSIS_HEADERFILES
     KTIterativeTrackClustering.hh
     KTMultiPeakEventBuilder.hh
     KTMultiPeakTrackProcessing.hh
+    KTCavityMultiPeakTrackProcessing.hh
     KTMultiPeakTrackBuilder.hh
     KTMultiSliceClustering.hh
     KTOverlappingTrackClustering.hh
@@ -73,6 +74,7 @@ set (EVENTANALYSIS_SOURCEFILES
     KTIterativeTrackClustering.cc
     KTMultiPeakEventBuilder.cc
     KTMultiPeakTrackProcessing.cc
+    KTCavityMultiPeakTrackProcessing.cc
     KTMultiPeakTrackBuilder.cc
     KTMultiSliceClustering.cc
     KTOverlappingTrackClustering.cc

Source/EventAnalysis/KTCavityMultiPeakTrackProcessing.cc

@@ -0,0 +1,183 @@
+/*
+ * KTCavityMultiPeakTrackProcessing.cc
+ *
+ *  Created on: Aug 19, 2025
+ *      Author: juniorpe
+ */
+
+#include "KTCavityMultiPeakTrackProcessing.hh"
+
+#include "KTLogger.hh"
+#include "KTProcessedTrackData.hh"
+#include "KTProcessedCavityMPTData.hh"
+#include "KTMultiTrackEventData.hh"
+
+#include <cmath>
+#include <set>
+#include <algorithm>
+
+namespace Katydid
+{
+    KTLOGGER(evlog, "KTCavityMultiPeakTrackProcessing");
+
+    // Register the processor
+    KT_REGISTER_PROCESSOR(KTCavityMultiPeakTrackProcessing, "cavity-mpt-processing");
+
+    KTCavityMultiPeakTrackProcessing::KTCavityMultiPeakTrackProcessing(const std::string& name) :
+            KTProcessor(name),
+            fProcessedCavityMPTSignal("proc-cavity-mpt", this),
+            fMPTSlot("mpt", this, &KTCavityMultiPeakTrackProcessing::AnalyzeMPT, &fProcessedCavityMPTSignal)
+    {
+    }
+
+    KTCavityMultiPeakTrackProcessing::~KTCavityMultiPeakTrackProcessing()
+    {
+    }
+
+    bool KTCavityMultiPeakTrackProcessing::Configure(const scarab::param_node* node)
+    {
+        if (node == NULL) return false;
+
+        return true;
+    }
+
+    bool KTCavityMultiPeakTrackProcessing::AnalyzeMPT( KTMultiPeakTrackData& mptData )
+    {
+        TrackSetCItSet allTracks = mptData.GetMPTrack().fTrackRefs;
+
+        // Determine multiplicity
+        int mult = mptData.GetMultiplicity();
+        KTDEBUG(evlog, "Determined multiplicity = " << mult);
+
+        if(mult <= 1)
+        {
+            KTWARN(evlog, "MPT only has one or no tracks; I cannot analyze this for axial frequency. Aborting");
+            return false;
+        }
+
+        // Create and fill new data object
+
+        KTProcessedCavityMPTData& procData = mptData.Of< KTProcessedCavityMPTData >();
+
+        procData.SetComponent( mptData.GetComponent() );
+        procData.SetAxialFrequency( 0. );
+
+        double axialFreq = 0.;
+        std::set<double> timeStampSet;
+        std::vector<double> frequencyDistances;
+
+        KTINFO(evlog, "Beginning axial frequency reconstruction. Attemping frequency distance calculation between bands in mpt.")
+        // Calculating all unique frequency distances between tracks at their start and end times
+        for (auto outerIt = allTracks.begin(); outerIt != allTracks.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, "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 = allTracks.begin(); innerIt != allTracks.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;
+                    }
+                }
+
+            }
+
+        }
+
+        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<double> 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 
+        {
+            KTWARN(evlog, "Frequency distances vector is empty or minimum value is zero. Aborting!")
+            return false;
+        }
+
+        if (frequencyDistances.size() != separationOrder.size() || frequencyDistances.empty()) 
+        {
+            KTWARN(evlog, "Frequency distances vector empty or has different size to separation order vector. Aborting!")
+            return false;
+        }
+
+        KTDEBUG(evlog, "Separation order of frequency distances:")
+        for (std::size_t i = 0; i < separationOrder.size(); ++i)
+        {
+            KTDEBUG(evlog, i+1 << " : " << separationOrder[i]);
+        }
+
+        double sum = 0.0;
+        for (std::size_t i = 0; i < frequencyDistances.size(); ++i) 
+        {
+            if (separationOrder[i] == 0.0) 
+            {
+                KTWARN(evlog, "Frequency distance division by 0 separation order. Aborting!")
+                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]));
+        }
+        axialFreq = sum/frequencyDistances.size();
+
+        // Set axial frequency
+        KTINFO(evlog, "Found axial frequency: " << axialFreq);
+        procData.SetAxialFrequency( axialFreq );
+
+        return true;
+    }
+
+} // namespace Katydid

Source/EventAnalysis/KTCavityMultiPeakTrackProcessing.hh

@@ -0,0 +1,66 @@
+/*
+ * KTCavityMultiPeakTrackProcessing.hh
+ *
+ *  Created on: Aug 19, 2025
+ *      Author: juniorpe
+ */
+
+#ifndef KTCAVITYMULTIPEAKTRACKPROCESSING_HH_
+#define KTCAVITYMULTIPEAKTRACKPROCESSING_HH_
+
+#include "KTProcessor.hh"
+#include "KTData.hh"
+#include "KTSlot.hh"
+
+namespace Katydid
+{
+    /*
+     @class KTCavityMultiPeakTrackProcessing
+     @author J. I. Pena
+     @brief Assigns axial frequency to MPT structure for symmetric-trap cavity TE011 data
+     @details
+     Iterates through bands in MPT to calculate frequency distance to other bands at distinct start and end times.
+     Assumes that trap is symmetric and carrier is present so that so sidebands are 2n(f_a) apart.
+
+     Available configuration values:
+     (none)
+
+     Slots:
+     - "mpt": void (Nymph::KTDataPtr) -- Analyzes a multi-peak track with > 1 band; Requires KTMultiPeakTrackData; Adds nothing
+
+     Signals:
+     - "proc-cavity-mpt": void (Nymph::KTDataPtr) -- Emitted upon successful determination of axial frequency; Guarantees KTProcessedCavityMPTData
+    */
+
+    class KTMultiPeakTrackData;
+
+    class KTCavityMultiPeakTrackProcessing : public Nymph::KTProcessor
+    {
+        public:
+            KTCavityMultiPeakTrackProcessing(const std::string& name = "cavity-mpt-processing");
+            virtual ~KTCavityMultiPeakTrackProcessing();
+
+            bool Configure(const scarab::param_node* node);
+
+        public:
+            bool AnalyzeMPT( KTMultiPeakTrackData& mptData );
+
+            //***************
+            // Signals
+            //***************
+
+        private:
+            Nymph::KTSignalData fProcessedCavityMPTSignal;
+
+            //***************
+            // Slots
+            //***************
+
+        private:
+            Nymph::KTSlotDataOneType< KTMultiPeakTrackData > fMPTSlot;
+
+    };
+
+}
+
+#endif /* KTCAVITYMULTIPEAKTRACKPROCESSING_HH_ */