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 + 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..56ec664cf --- /dev/null +++ b/Source/Data/EventAnalysis/KTProcessedCavityEventData.cc @@ -0,0 +1,81 @@ +/* + * 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), + fFirstTrackStartCyclotronFrequency(0.), + fFirstTrackAxialFrequency(0.), + fFirstTrackBandIDs(), + fFirstTrackBandClassifications(), + fMPTEventSequenceID(), + fMPTStartCyclotronFrequency(), + fMPTAxialFrequency(), + fAllTrackIDs(), + fAllTrackEventSequenceIDs(), + fAllTrackBandClassifications() + + { + } + + KTProcessedCavityEventData::KTProcessedCavityEventData(const KTProcessedCavityEventData& orig) : + KTExtensibleData< KTProcessedCavityEventData >(orig), + + fComponent(orig.fComponent), + fAcquisitionID(orig.fAcquisitionID), + fEventID(orig.fEventID), + fTotalEventSequences(orig.fTotalEventSequences), + fFirstTrackStartCyclotronFrequency(orig.fFirstTrackStartCyclotronFrequency), + fFirstTrackAxialFrequency(orig.fFirstTrackAxialFrequency), + fFirstTrackBandIDs(orig.fFirstTrackBandIDs), + fFirstTrackBandClassifications(orig.fFirstTrackBandClassifications), + fMPTEventSequenceID(orig.fMPTEventSequenceID), + fMPTStartCyclotronFrequency(orig.fMPTStartCyclotronFrequency), + fMPTAxialFrequency(orig.fMPTAxialFrequency), + fAllTrackIDs(orig.fAllTrackIDs), + fAllTrackEventSequenceIDs(orig.fAllTrackEventSequenceIDs), + fAllTrackBandClassifications(orig.fAllTrackBandClassifications) + { + } + + KTProcessedCavityEventData::~KTProcessedCavityEventData() + { + } + + KTProcessedCavityEventData& KTProcessedCavityEventData::operator=(const KTProcessedCavityEventData& rhs) + { + KTExtensibleData< KTProcessedCavityEventData >::operator=(rhs); + fComponent = rhs.fComponent; + fAcquisitionID = rhs.fAcquisitionID; + fEventID = rhs.fEventID; + fTotalEventSequences = rhs.fTotalEventSequences; + fFirstTrackStartCyclotronFrequency = rhs.fFirstTrackStartCyclotronFrequency; + fFirstTrackAxialFrequency = rhs.fFirstTrackAxialFrequency; + fFirstTrackBandIDs = rhs.fFirstTrackBandIDs; + fFirstTrackBandClassifications = rhs.fFirstTrackBandClassifications; + 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 new file mode 100644 index 000000000..eb3889d41 --- /dev/null +++ b/Source/Data/EventAnalysis/KTProcessedCavityEventData.hh @@ -0,0 +1,112 @@ +/* + * KTProcessedCavityEventData.hh + * + * Created on: Sep 11, 2025 + * Author: juniorpe + */ + +#ifndef KTPROCESSEDCAVITYEVENTDATA_HH_ +#define KTPROCESSEDCAVITYEVENTDATA_HH_ + +#include "KTData.hh" +#include "KTMultiTrackEventData.hh" +#include "KTMemberVariable.hh" + +#include +#include + +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, FirstTrackStartCyclotronFrequency); + MEMBERVARIABLE(double, FirstTrackAxialFrequency); + MEMBERVARIABLE(std::vector, FirstTrackBandIDs); + MEMBERVARIABLE(std::vector, FirstTrackBandClassifications); + + 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 ClearProcessedEvent(); + + void AddFirstTrackBandInfo(int trackID, int classification); + void ClearFirstTrackBandInfo(); + + void AddProcessedMPT(int seqID, double startCyclFreq, double axialFreq); + + void AddClassificationData(int trackID, int seqID, int bandClassification); + void AddClassificationData(const std::vector& tracksInMPT, int seqID, int bandClassification); + + public: + static const std::string sName; + }; + + inline void KTProcessedCavityEventData::ClearProcessedEvent() + { + fFirstTrackBandIDs.clear(); + fFirstTrackBandClassifications.clear(); + + fMPTEventSequenceID.clear(); + fMPTStartCyclotronFrequency.clear(); + fMPTAxialFrequency.clear(); + + fAllTrackIDs.clear(); + fAllTrackEventSequenceIDs.clear(); + 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); + 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); + } + } + +} +#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..337974f67 --- /dev/null +++ b/Source/EventAnalysis/KTCavityEventProcessing.cc @@ -0,0 +1,420 @@ +/* + * KTCavityEventProcessing.cc + * + * Created on: Sep 11, 2025 + * Author: juniorpe + */ + +#include "KTCavityEventProcessing.hh" + +#include "KTLogger.hh" +#include "KTProcessedTrackData.hh" +#include "KTMultiTrackEventData.hh" +#include "KTProcessedCavityEventData.hh" + +#include +#include + + +namespace Katydid +{ + KTLOGGER(evlog, "KTCavityEventProcessing"); + + // Register the processor + KT_REGISTER_PROCESSOR(KTCavityEventProcessing, "cavity-event-processing"); + + KTCavityEventProcessing::KTCavityEventProcessing(const std::string& name) : + KTProcessor(name), + fMaxRelPowerThresh(0.), + 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; + + SetMaxRelPowerThresh(node->get_value("max-rel-power-thresh", GetMaxRelPowerThresh())); + + return true; + } + + bool KTCavityEventProcessing::AnalyzeEvent( KTMultiTrackEventData& mtEventData ) + { + // Create and fill new data object + KTProcessedCavityEventData& procEvent = mtEventData.Of(); + + procEvent.ClearProcessedEvent(); + + procEvent.SetComponent(mtEventData.GetComponent()); + procEvent.SetAcquisitionID(mtEventData.GetAcquisitionID()); + procEvent.SetEventID(mtEventData.GetEventID()); + procEvent.SetTotalEventSequences(mtEventData.GetTotalEventSequences()); + procEvent.SetFirstTrackStartCyclotronFrequency(-1.); + procEvent.SetFirstTrackAxialFrequency(-1.); + + + // Storing processed tracks in map by 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); + } + + 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, procEvent, 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."); + } + + // Storing reconstructed frequencies for only first multi-peak-track + if (eventSeqID == 0) + { + KTINFO(evlog, "Found first track initial cyclotron frequency: " << reconstructedStartCyclotronFrequency); + procEvent.SetFirstTrackStartCyclotronFrequency(reconstructedStartCyclotronFrequency); + + KTINFO(evlog, "Found first track axial frequency: " << reconstructedAxialFrequency); + procEvent.SetFirstTrackAxialFrequency(reconstructedAxialFrequency); + } + + procEvent.AddProcessedMPT(eventSeqID, reconstructedStartCyclotronFrequency, reconstructedAxialFrequency); + + + } + + return true; + } + + bool KTCavityEventProcessing::ReconstructCyclotronFrequency(const std::vector& tracksInMPT, KTProcessedCavityEventData& procEvent, 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->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 + // 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.GetStartTimeInRunC() < startTime) {startTime=track.fProcTrack.GetStartTimeInRunC();} + } + std::sort(sortedMPTBands.begin(), sortedMPTBands.end(), + [](const AllTrackData* a, const AllTrackData* b) + { + return a->fProcTrack.GetStartFrequency() < b->fProcTrack.GetStartFrequency(); + }); + + 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 + // 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) + // Inital cyclotron frequency extracted from carrier with sideband order 0 by determine its position depending on topology + int carrierIndex = -1; + std::vector bandClassification; + + if (numBands==1) + { + carrierIndex = 0; + bandClassification = {0}; + KTDEBUG(evlog, "MPT has 1 band!"); + } + else if (numBands==2) + { + // 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] + const AllTrackData* lastTrack = sortedMPTBands[numBands-1]; + if (maxTotNUP!=0) + { + double lastTrackRelTotalNUP = lastTrack->fProcTrack.GetTotalTrackNUP()/maxTotNUP; + if (lastTrackRelTotalNUP> this->fMaxRelPowerThresh) + { + 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]."); + + } + + } + 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; + bandClassification = {-4, -2, 0, 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]; + 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; + } + + 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)"); + + // 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]; + if (track->fProcTrack.GetEventSequenceID() != seqID) + { + KTWARN(evlog, "Sorted bands contain multiple EventSequenceIDs(expected " << seqID << "), something upstream is wrong. Aborting"); + return false; + } + 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; + } + + bool KTCavityEventProcessing::ReconstructAxialFrequency(const std::vector& tracksInMPT, double& outAxialFrequency) const + { + outAxialFrequency = -1.0; + if (tracksInMPT.size() <= 1) + { + 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) + { + 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 = 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; + } + } + + } + + } + + if (frequencyDistances.empty()) + { + KTWARN(evlog, "No valid frequency distances found; cannot reconstruct axial frequency."); + return false; + } + + 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; + } + 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()) + { + 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(); + KTDEBUG(evlog, "Reconstructed Axial Frequency: " << outAxialFrequency); + + + return true; + } + +} // namespace Katydid diff --git a/Source/EventAnalysis/KTCavityEventProcessing.hh b/Source/EventAnalysis/KTCavityEventProcessing.hh new file mode 100644 index 000000000..a56aa663b --- /dev/null +++ b/Source/EventAnalysis/KTCavityEventProcessing.hh @@ -0,0 +1,79 @@ +/* + * KTCavityEventProcessing.hh + * + * Created on: Sep 11, 2025 + * Author: juniorpe + */ + +#ifndef KTCAVITYEVENTPROCESSING_HH_ +#define KTCAVITYEVENTPROCESSING_HH_ + +#include "KTProcessor.hh" +#include "KTData.hh" +#include "KTSlot.hh" +#include "KTMemberVariable.hh" +#include "KTProcessedTrackData.hh" +#include "KTProcessedCavityEventData.hh" + +namespace Katydid +{ + /* + @class KTCavityEventProcessing + @author J. I. Pena + @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 + + Slots: + - "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 + */ + + 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); + + MEMBERVARIABLE(float, MaxRelPowerThresh); + + public: + bool AnalyzeEvent( KTMultiTrackEventData& mtEventData ); + + bool ReconstructCyclotronFrequency(const std::vector& tracksInMPT, KTProcessedCavityEventData& procEvent, double& outStartCyclotronFrequency) const; + bool ReconstructAxialFrequency(const std::vector& tracksInMPT, double& outAxialFrequency) const; + + + //*************** + // 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..afb80471b 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" @@ -62,6 +63,7 @@ namespace Katydid fProcessedTrackTree(NULL), fMultiPeakTrackTree(NULL), fMultiTrackEventTree(NULL), + fProcessedCavityEventTree(NULL), fMTEWithClassifierResultsTree(NULL), fLinearFitResultTree(NULL), fPowerFitDataTree(NULL), @@ -71,6 +73,7 @@ namespace Katydid fSequentialLineDataPtr(NULL), fProcessedTrackDataPtr(NULL), fProcessedMPTDataPtr(NULL), + fProcessedCavityEventData(), 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-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 +685,106 @@ 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.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(); + 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()); + + 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( "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); + fProcessedCavityEventTree->SetBranchAddress( "AllTrackIDs", &fProcessedCavityEventData.fAllTrackIDs); + fProcessedCavityEventTree->SetBranchAddress( "AllTrackEventSequenceIDs", &fProcessedCavityEventData.fAllTrackEventSequenceIDs); + fProcessedCavityEventTree->SetBranchAddress( "AllTrackBandClassifications", &fProcessedCavityEventData.fAllTrackBandClassifications); + + 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( "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); + fProcessedCavityEventTree->Branch( "AllTrackIDs", "std::vector", &fProcessedCavityEventData.fAllTrackIDs); + fProcessedCavityEventTree->Branch( "AllTrackEventSequenceIDs", "std::vector", &fProcessedCavityEventData.fAllTrackEventSequenceIDs); + fProcessedCavityEventTree->Branch( "AllTrackBandClassifications", "std::vector", &fProcessedCavityEventData.fAllTrackBandClassifications); + + return true; + } + + //************************** // Multi-Peak Track //************************** diff --git a/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh b/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh index 30348a4a5..b61aac6b4 100644 --- a/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh +++ b/Source/IO/ROOTTreeWriter/KTROOTTreeTypeWriterEventAnalysis.hh @@ -116,6 +116,27 @@ namespace Katydid UInt_t fUnknownEventTopology; }; + struct TProcessedCavityEventData + { + UInt_t fComponent; + ULong64_t fAcquisitionID; + UInt_t fEventID; + UInt_t fTotalEventSequences; + Double_t fFirstTrackStartCyclotronFrequency; + Double_t fFirstTrackAxialFrequency; + std::vector fFirstTrackBandIDs; + std::vector fFirstTrackBandClassifications; + + std::vector fMPTEventSequenceID; + std::vector fMPTStartCyclotronFrequency; + std::vector fMPTAxialFrequency; + + std::vector fAllTrackIDs; + std::vector fAllTrackEventSequenceIDs; + std::vector fAllTrackBandClassifications; + + }; + struct TLinearFitResult { UInt_t fFitNumber; @@ -198,6 +219,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 +233,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 +247,7 @@ namespace Katydid bool SetupSparseWaterfallCandidateTree(); bool SetupSequentialLineTree(); bool SetupProcessedMPTTree(); + bool SetupProcessedCavityEventTree(); bool SetupProcessedTrackTree(); bool SetupMultiPeakTrackTree(); bool SetupMultiTrackEventTree(); @@ -236,6 +260,7 @@ namespace Katydid TTree* fSparseWaterfallCandidateTree; TTree* fSequentialLineTree; TTree* fProcessedMPTTree; + TTree* fProcessedCavityEventTree; TTree* fProcessedTrackTree; TTree* fMultiPeakTrackTree; TTree* fMultiTrackEventTree; @@ -249,6 +274,7 @@ namespace Katydid TSequentialLineData* fSequentialLineDataPtr; Cicada::TProcessedTrackData* fProcessedTrackDataPtr; Cicada::TProcessedMPTData* fProcessedMPTDataPtr; + TProcessedCavityEventData fProcessedCavityEventData; TMultiPeakTrackData fMultiPeakTrackData; Cicada::TMultiTrackEventData* fMultiTrackEventDataPtr; Cicada::TMTEWithClassifierResultsData* fMTEWithClassifierResultsDataPtr; @@ -282,6 +308,11 @@ namespace Katydid return fProcessedMPTTree; } + inline TTree* KTROOTTreeTypeWriterEventAnalysis::GetProcessedCavityEventTree() const + { + return fProcessedCavityEventTree; + } + inline TTree* KTROOTTreeTypeWriterEventAnalysis::GetProcessedTrackTree() const { return fProcessedTrackTree;