diff --git a/Source/Data/CMakeLists.txt b/Source/Data/CMakeLists.txt index 52fd90b83..e2318efad 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/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 diff --git a/Source/Data/EventAnalysis/KTProcessedCavityMPTData.cc b/Source/Data/EventAnalysis/KTProcessedCavityMPTData.cc new file mode 100644 index 000000000..2798feba2 --- /dev/null +++ b/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; + } + +} diff --git a/Source/Data/EventAnalysis/KTProcessedCavityMPTData.hh b/Source/Data/EventAnalysis/KTProcessedCavityMPTData.hh new file mode 100644 index 000000000..39587d022 --- /dev/null +++ b/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 + diff --git a/Source/EventAnalysis/CMakeLists.txt b/Source/EventAnalysis/CMakeLists.txt index bea0b58c0..b2914d078 100644 --- a/Source/EventAnalysis/CMakeLists.txt +++ b/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 diff --git a/Source/EventAnalysis/KTCavityMultiPeakTrackProcessing.cc b/Source/EventAnalysis/KTCavityMultiPeakTrackProcessing.cc new file mode 100644 index 000000000..345e50cb7 --- /dev/null +++ b/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 +#include +#include + +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 timeStampSet; + std::vector 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 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 diff --git a/Source/EventAnalysis/KTCavityMultiPeakTrackProcessing.hh b/Source/EventAnalysis/KTCavityMultiPeakTrackProcessing.hh new file mode 100644 index 000000000..683e04752 --- /dev/null +++ b/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_ */ diff --git a/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.cc b/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.cc index d4d00992c..92f80c8d6 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 "KTProcessedCavityMPTData.hh" #include "KTProcessedTrackData.hh" #include "KTSliceHeader.hh" #include "KTSparseWaterfallCandidateData.hh" @@ -59,6 +60,7 @@ namespace Katydid fSparseWaterfallCandidateTree(NULL), fSequentialLineTree(NULL), fProcessedMPTTree(NULL), + fProcessedCavityMPTTree(NULL), fProcessedTrackTree(NULL), fMultiPeakTrackTree(NULL), fMultiTrackEventTree(NULL), @@ -71,6 +73,7 @@ namespace Katydid fSequentialLineDataPtr(NULL), fProcessedTrackDataPtr(NULL), fProcessedMPTDataPtr(NULL), + fProcessedCavityMPTData(), fMultiPeakTrackData(), fMultiTrackEventDataPtr(NULL), fMTEWithClassifierResultsDataPtr(NULL), @@ -95,6 +98,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-mpt", this, &KTROOTTreeTypeWriterEventAnalysis::WriteProcessedCavityMPT); fWriter->RegisterSlot("proc-track", this, &KTROOTTreeTypeWriterEventAnalysis::WriteProcessedTrack); fWriter->RegisterSlot("mp-track", this, &KTROOTTreeTypeWriterEventAnalysis::WriteMultiPeakTrack); fWriter->RegisterSlot("mt-event", this, &KTROOTTreeTypeWriterEventAnalysis::WriteMultiTrackEvent); @@ -681,6 +685,67 @@ namespace Katydid return true; } + //**************************** + // Processed Cavity Multi-Peak Track + //**************************** + + void KTROOTTreeTypeWriterEventAnalysis::WriteProcessedCavityMPT(Nymph::KTDataPtr data) + { + KTDEBUG(publog, "Attempting to write to processed cavity mpt root tree"); + KTProcessedCavityMPTData& procCavityMPTData = data->Of< KTProcessedCavityMPTData >(); + + if (! fWriter->OpenAndVerifyFile()) return; + + if (fProcessedCavityMPTTree == NULL) + { + if (! SetupProcessedCavityMPTTree()) + { + KTERROR(publog, "Something went wrong while setting up the processed cavity mpt tree! Nothing was written."); + return; + } + } + + fProcessedCavityMPTData.fComponent = procCavityMPTData.GetComponent(); + fProcessedCavityMPTData.fAxialFrequency = procCavityMPTData.GetAxialFrequency(); + + fProcessedCavityMPTTree->Fill(); + + return; + } + + bool KTROOTTreeTypeWriterEventAnalysis::SetupProcessedCavityMPTTree() + { + if( fWriter->GetAccumulate() ) + { + fWriter->GetFile()->GetObject( "processed-cavity-mpt", fProcessedCavityMPTTree ); + + if (fProcessedCavityMPTTree != NULL) + { + KTINFO(publog, "Tree already exists; will add to it"); + fWriter->AddTree( fProcessedCavityMPTTree ); + + fProcessedCavityMPTTree->SetBranchAddress( "Component", &fProcessedCavityMPTData.fComponent ); + fProcessedCavityMPTTree->SetBranchAddress( "AxialFrequency", &fProcessedCavityMPTData.fAxialFrequency ); + + return true; + } + } + + fProcessedCavityMPTTree = new TTree("processed-cavity-mpt", "Processed Cavity MPT"); + if( fProcessedCavityMPTTree == NULL ) + { + KTERROR( publog, "Tree was not created!" ); + return false; + } + fWriter->AddTree( fProcessedCavityMPTTree ); + + fProcessedCavityMPTTree->Branch( "Component", &fProcessedCavityMPTData.fComponent, "fComponent/i" ); + fProcessedCavityMPTTree->Branch( "AxialFrequency", &fProcessedCavityMPTData.fAxialFrequency, "fAxialFrequency/d" ); + + return true; + } + + //************************** // Multi-Peak Track //************************** diff --git a/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh b/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh index 30348a4a5..e962db7ec 100644 --- a/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh +++ b/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh @@ -116,6 +116,12 @@ namespace Katydid UInt_t fUnknownEventTopology; }; + struct TProcessedCavityMPTData + { + UInt_t fComponent; + Double_t fAxialFrequency; + }; + struct TLinearFitResult { UInt_t fFitNumber; @@ -198,6 +204,7 @@ namespace Katydid void WriteSparseWaterfallCandidate(Nymph::KTDataPtr data); void WriteSequentialLine(Nymph::KTDataPtr data); void WriteProcessedMPT(Nymph::KTDataPtr data); + void WriteProcessedCavityMPT(Nymph::KTDataPtr data); void WriteProcessedTrack(Nymph::KTDataPtr data); void WriteMultiPeakTrack(Nymph::KTDataPtr data); void WriteMultiTrackEvent(Nymph::KTDataPtr data); @@ -211,6 +218,7 @@ namespace Katydid TTree* GetSparseWaterfallCandidateTree() const; TTree* GetSequentialLineTree() const; TTree* GetProcessedMPTTree() const; + TTree* GetProcessedCavityMPTTree() const; TTree* GetProcessedTrackTree() const; TTree* GetMultiPeakTrackTree() const; TTree* GetMultiTrackEventTree() const; @@ -224,6 +232,7 @@ namespace Katydid bool SetupSparseWaterfallCandidateTree(); bool SetupSequentialLineTree(); bool SetupProcessedMPTTree(); + bool SetupProcessedCavityMPTTree(); bool SetupProcessedTrackTree(); bool SetupMultiPeakTrackTree(); bool SetupMultiTrackEventTree(); @@ -236,6 +245,7 @@ namespace Katydid TTree* fSparseWaterfallCandidateTree; TTree* fSequentialLineTree; TTree* fProcessedMPTTree; + TTree* fProcessedCavityMPTTree; TTree* fProcessedTrackTree; TTree* fMultiPeakTrackTree; TTree* fMultiTrackEventTree; @@ -249,6 +259,7 @@ namespace Katydid TSequentialLineData* fSequentialLineDataPtr; Cicada::TProcessedTrackData* fProcessedTrackDataPtr; Cicada::TProcessedMPTData* fProcessedMPTDataPtr; + TProcessedCavityMPTData fProcessedCavityMPTData; TMultiPeakTrackData fMultiPeakTrackData; Cicada::TMultiTrackEventData* fMultiTrackEventDataPtr; Cicada::TMTEWithClassifierResultsData* fMTEWithClassifierResultsDataPtr; @@ -282,6 +293,11 @@ namespace Katydid return fProcessedMPTTree; } + inline TTree* KTROOTTreeTypeWriterEventAnalysis::GetProcessedCavityMPTTree() const + { + return fProcessedCavityMPTTree; + } + inline TTree* KTROOTTreeTypeWriterEventAnalysis::GetProcessedTrackTree() const { return fProcessedTrackTree;