diff --git a/Source/SpectrumAnalysis/CMakeLists.txt b/Source/SpectrumAnalysis/CMakeLists.txt index 48accce6a..5049f6df7 100644 --- a/Source/SpectrumAnalysis/CMakeLists.txt +++ b/Source/SpectrumAnalysis/CMakeLists.txt @@ -18,6 +18,7 @@ set (SPECTRUMANALYSIS_HEADERFILES KTHoughTransform.hh KTMergeKDTree.hh KTNNFilter.hh + KTPhasedAggregator.hh KTSequentialTrackFinder.hh KTSpectrumDiscriminator.hh KTSpectrogramStriper.hh @@ -43,6 +44,7 @@ set (SPECTRUMANALYSIS_SOURCEFILES KTHoughTransform.cc KTMergeKDTree.cc KTNNFilter.cc + KTPhasedAggregator.cc KTSequentialTrackFinder.cc KTSpectrumDiscriminator.cc KTSpectrogramStriper.cc diff --git a/Source/SpectrumAnalysis/KTChannelAggregator.cc b/Source/SpectrumAnalysis/KTChannelAggregator.cc index 965057a57..b7a7ceac7 100644 --- a/Source/SpectrumAnalysis/KTChannelAggregator.cc +++ b/Source/SpectrumAnalysis/KTChannelAggregator.cc @@ -246,6 +246,7 @@ namespace Katydid double gridLocationY = 0; double gridLocationZ = 0; newAggFreqData.GetGridPoint(gridPointNumber, gridLocationX, gridLocationY, gridLocationZ); + std::cout << "nComponents: " << nComponents << ", nRings: " << fNRings << std::endl; for (unsigned iComponent = 0; iComponent < nComponents; ++iComponent) { // Arbitarily assign 0 to the first channel and progresively add 2pi/N for the rest of the channels in increasing order diff --git a/Source/SpectrumAnalysis/KTPhasedAggregator.cc b/Source/SpectrumAnalysis/KTPhasedAggregator.cc new file mode 100644 index 000000000..eb05140db --- /dev/null +++ b/Source/SpectrumAnalysis/KTPhasedAggregator.cc @@ -0,0 +1,309 @@ +/* + * KTPhasedAggregator.cc + * + * Created on: Dec 13, 2024 + * Author: J. K. Gaison + * Based on KTChannelAggregator.cc by P. T. Surukuchi + */ + +#include "KTPhasedAggregator.hh" +#include "KTLogger.hh" + +#include +using namespace boost::algorithm; + +#include + +namespace Katydid +{ + KTLOGGER(agglog, "KTPhasedAggregator"); + + // Register the processor + KT_REGISTER_PROCESSOR(KTPhasedAggregator, "phased-aggregator"); + + KTPhasedAggregator::KTPhasedAggregator(const std::string& name) : + KTProcessor(name), + fSummedFrequencyData("agg-fft", this), + fPhaseChFrequencySumSlot("fft", this, &KTPhasedAggregator::SumChannelVoltageWithPhase, &fSummedFrequencyData), + fAxialSumSlot("ax-agg-fft", this, &KTPhasedAggregator::SumChannelVoltageWithPhase, &fSummedFrequencyData), + fActiveRadius(0), + fNGrid(0), + fWavelength(0), + fIsGridDefined(false), + fIsUserDefinedGrid(false), + fIsPartialRing(false), + fPartialRingMultiplicity(), + fSummationMinFreq(0e6), + fSummationMaxFreq(200e6), + fUseAntiSpiralPhaseShifts(false), + fAntiSpiralPhaseShifts(), + fPhase_str("0."), + fScale_str("1."), + fNRings(1), + fPhases({0}), + fScales({1.}) + { + } + + KTPhasedAggregator::~KTPhasedAggregator() + { + } + + bool KTPhasedAggregator::Configure(const scarab::param_node* node) + { + if (node != NULL) + { + fNGrid = node->get_value< unsigned >("grid-size", fNGrid); + fIsUserDefinedGrid = node->get_value< bool >("use-grid-file", fIsUserDefinedGrid); + fIsPartialRing= node->get_value< bool >("partial-ring", fIsPartialRing); + fPartialRingMultiplicity= node->get_value< unsigned >("partial-ring-Multiplicity", fPartialRingMultiplicity); + fUserDefinedGridFile = node->get_value< >("grid-file", fUserDefinedGridFile); + fActiveRadius = node->get_value< double >("active-radius", fActiveRadius); + fWavelength = node->get_value< double >("wavelength", fWavelength); + fSummationMinFreq= node->get_value< double >("min-freq", fSummationMinFreq); + fSummationMaxFreq= node->get_value< double >("max-freq", fSummationMaxFreq); + fNRings = node->get_value< unsigned >("n-rings", fNRings); + fPhase_str = node->get_value< std::string >("phases", fPhase_str); + std::vector sPhase; + boost::split(sPhase, fPhase_str, boost::is_any_of(", "), boost::token_compress_on); + fPhases.resize(sPhase.size()); + fScale_str = node->get_value< std::string >("scales", fScale_str); + std::vector sScale; + boost::split(sScale, fScale_str, boost::is_any_of(", "), boost::token_compress_on); + if(sPhase.size() != sScale.size()) + { + KTERROR(agglog,"The imported phases and scales are not the same size."); + } + fPhases.resize(sPhase.size()); + fScales.resize(sScale.size()); + for(int i=0; iget_value< bool>("use-antispiral-phase-shifts", fUseAntiSpiralPhaseShifts); + } + return true; + } + + bool KTPhasedAggregator::ApplyPhaseShift(double &realVal, double &imagVal, double phase) + { + double tempRealVal = realVal; + double tempImagVal = imagVal; + realVal = tempRealVal * cos(phase) - tempImagVal * sin(phase); + imagVal = tempRealVal * sin(phase) + tempImagVal * cos(phase); + return true; + } + + double KTPhasedAggregator::GetPhaseShift(double xPosition, double yPosition, double wavelength, double channelAngle) const + { + // X position based on the angle of the channel + double xChannel = fActiveRadius * cos(channelAngle); + // X position based on the angle of the channel + double yChannel = fActiveRadius * sin(channelAngle); + // Distance of the input point from the input channel + double pointDistance = pow(pow(xChannel - xPosition, 2) + pow(yChannel - yPosition, 2), 0.5); + // Phase of the input signal based on the input point, channel location and the wavelength + return 2.0 * KTMath::Pi() * pointDistance / wavelength; + } + + double KTPhasedAggregator::GetAntiSpiralPhaseShift(double xPosition, double yPosition, double wavelength, double channelAngle) const + { + // X position based on the angle of the channel + double xChannel = fActiveRadius * cos(channelAngle); + // X position based on the angle of the channel + double yChannel = fActiveRadius * sin(channelAngle); + // Angle based on the deltaX and deltaY + return atan2(yChannel-yPosition,xChannel-xPosition); + } + + bool KTPhasedAggregator::GetGridLocation(unsigned gridNumber, unsigned gridSize, double &gridLocation) + { + if (gridNumber >= gridSize) return false; + gridLocation = fActiveRadius * (((2.0 * gridNumber + 1.0) / gridSize) - 1); + return true; + } + + bool KTPhasedAggregator::GenerateAntiSpiralPhaseShifts(unsigned channelCount) + { + for(unsigned i=0;i channelPhaseShift=std::make_pair(i,phaseShift); + fAntiSpiralPhaseShifts.insert(channelPhaseShift); + } + return true; + } + + bool KTPhasedAggregator::SumChannelVoltageWithPhase(KTFrequencySpectrumDataFFTW& fftwData) + { + KTAggregatedFrequencySpectrumDataFFTW& newAggFreqData = fftwData.Of< KTAggregatedFrequencySpectrumDataFFTW >().SetNComponents(0); + return PerformPhaseSummation(fftwData, newAggFreqData); + } + + bool KTPhasedAggregator::SumChannelVoltageWithPhase(KTAxialAggregatedFrequencySpectrumDataFFTW& fftwData) + { + KTAggregatedFrequencySpectrumDataFFTW& newAggFreqData = fftwData.Of< KTAggregatedFrequencySpectrumDataFFTW >().SetNComponents(0); + return PerformPhaseSummation(fftwData, newAggFreqData); + } + unsigned KTPhasedAggregator::DefineGrid(KTAggregatedFrequencySpectrumDataFFTW &newAggFreqData) + { + unsigned nTotalGridPoints=0; + for (unsigned iRing = 0; iRing < fNRings; ++iRing) + { + if(fIsUserDefinedGrid) + { + if( !ends_with(fUserDefinedGridFile, ".txt") ) + { + //if(!ends_with(fTextFileName.c_str(),".txt")) + KTDEBUG(agglog,"`grid-file` file must end in .txt"); + return -1; + } + double gridLocationX; + double gridLocationY; + std::fstream textFile(fUserDefinedGridFile.c_str(),std::ios::in); + if (textFile.fail()) + { + KTDEBUG(agglog,"`grid-file` cannot be opened"); + return -1; + } + while(!textFile.eof()) + { + std::string lineContent; + while(std::getline(textFile,lineContent)) + { + if (lineContent.find('#')!=std::string::npos) continue; + std::string token; + std::stringstream ss(lineContent); + unsigned wordCount=0; + while (ss >> token) + { + if( wordCount==0 ) gridLocationX = std::stod(token); + else if( wordCount==1 ) gridLocationY = std::stod(token); + else + { + KTDEBUG(agglog,"`grid-file` cannot have more than 2 columns"); + return -1; + } + ++wordCount; + } + // Check to make sure that the grid point is within the active detector volume, skip otherwise + if((pow(gridLocationX, 2) + pow( gridLocationY,2 ) )> pow(fActiveRadius, 2)) continue; + newAggFreqData.SetNComponents(nTotalGridPoints+1); + newAggFreqData.SetGridPoint(nTotalGridPoints, gridLocationX, gridLocationY, iRing); + ++nTotalGridPoints; + } + } + } + else + { + // Loop over the grid points and fill the values + for (unsigned iGridX = 0; iGridX < fNGrid; ++iGridX) + { + double gridLocationX = 0; + GetGridLocation(iGridX, fNGrid, gridLocationX); + for (unsigned iGridY = 0; iGridY < fNGrid; ++iGridY) + { + double gridLocationY = 0; + GetGridLocation(iGridY, fNGrid, gridLocationY); + // Check to make sure that the grid point is within the active detector volume, skip otherwise + if( (pow(gridLocationX,2)+pow(gridLocationY,2))>pow(fActiveRadius,2) ) continue; + newAggFreqData.SetNComponents(nTotalGridPoints+1); + newAggFreqData.SetGridPoint(nTotalGridPoints, gridLocationX, gridLocationY, iRing); + ++nTotalGridPoints; + } + } + } + } + return nTotalGridPoints; + } + bool KTPhasedAggregator::PerformPhaseSummation(KTFrequencySpectrumDataFFTWCore& fftwData,KTAggregatedFrequencySpectrumDataFFTW &newAggFreqData) + { + const KTFrequencySpectrumFFTW* freqSpectrum = fftwData.GetSpectrumFFTW(0); + unsigned nTimeBins = freqSpectrum->GetNTimeBins(); + // Get the number of frequency bins from the first component of fftwData + unsigned nFreqBins = freqSpectrum->GetNFrequencyBins(); + unsigned nTotalComponents = fftwData.GetNComponents(); // Get number of components + if( fIsPartialRing ) + { + if( fPartialRingMultiplicity%2!=0 || fPartialRingMultiplicity<=0 ) return false; // The partial ring case should only work if the multiplicity is even + nTotalComponents*=fPartialRingMultiplicity; + } + if( nTotalComponents%fNRings!=0 ) + { + KTERROR(agglog,"The number of rings has to be an integer multiple of total components"); + } + unsigned nComponents = nTotalComponents/fNRings;// Get number of components + + GenerateAntiSpiralPhaseShifts(nComponents); + double maxValue = 0.0; + double maxGridLocationX = 0.0; + double maxGridLocationY = 0.0; + // Setting up the active radius of the KTAggregatedFrequencySpectrumDataFFTW object to maintain consistency + // This doesn't need to be done if there is a way to provide config values to data objects + newAggFreqData.SetActiveRadius(fActiveRadius); + // Set the number of rings present + newAggFreqData.SetNAxialPositions(fNRings); + unsigned nTotalGridPoints = DefineGrid(newAggFreqData); + if(nTotalGridPoints<=0) return false; + unsigned gridPointsPerRing=nTotalGridPoints/fNRings; + // Loop over the grid points and rings and fill the values + for (unsigned iRing = 0; iRing < fNRings; ++iRing) + { + // Loop over all grid points and find the one that gives the highest value + for (unsigned iGrid = 0; iGrid < gridPointsPerRing; ++iGrid) + { // Loop over the grid points + unsigned gridPointNumber=iGrid+gridPointsPerRing*iRing; + KTFrequencySpectrumFFTW* newFreqSpectrum = new KTFrequencySpectrumFFTW(nFreqBins, freqSpectrum->GetRangeMin(), freqSpectrum->GetRangeMax()); + // Empty values in the frequency spectrum, not sure if this is needed but there were some issues when this was not done for the power spectrum + NullFreqSpectrum(*newFreqSpectrum); + double gridLocationX = 0; + double gridLocationY = 0; + double gridLocationZ = 0; + newAggFreqData.GetGridPoint(gridPointNumber, gridLocationX, gridLocationY, gridLocationZ); + for (unsigned iComponent = 0; iComponent < nComponents; ++iComponent) + { + // Arbitarily assign 0 to the first channel and progresively add 2pi/N for the rest of the channels in increasing order + double PhaseShift = fPhases[fNRings*iRing + iComponent]; + double ScaleCoefficient = fScales[fNRings*iRing + iComponent]; + // Just being redundantly cautious, the phaseShifts are already zerors but checking to make sure anyway + freqSpectrum = fftwData.GetSpectrumFFTW(iComponent+iRing*nComponents); + double maxVoltage = 0.0; + unsigned maxFrequencyBin = 0; + //Loop over the frequency bins + for (unsigned iFreqBin = 0; iFreqBin < nFreqBins; ++iFreqBin) + { + if( newFreqSpectrum->GetBinCenter(iFreqBin)GetBinCenter(iFreqBin)>fSummationMaxFreq ) continue; + double realVal = freqSpectrum->GetReal(iFreqBin); + double imagVal = freqSpectrum->GetImag(iFreqBin); + ApplyPhaseShift(realVal, imagVal, PhaseShift); + double summedRealVal = ScaleCoefficient*realVal + newFreqSpectrum->GetReal(iFreqBin); + double summedImagVal = ScaleCoefficient*imagVal + newFreqSpectrum->GetImag(iFreqBin); + (*newFreqSpectrum)(iFreqBin) = std::complex{ summedRealVal, summedImagVal }; + } // End of loop over freq bins + } // End of loop over all comps + newFreqSpectrum->SetNTimeBins(nTimeBins); + newAggFreqData.SetSpectrum(newFreqSpectrum,gridPointNumber); + + double maxVoltageFreq = 0.0; + //Loop over all the freq bins and get the highest value and save to the aggregated frequency data + for (unsigned iFreqBin = 0; iFreqBin < nFreqBins; ++iFreqBin) + { + if (newFreqSpectrum->GetAbs(iFreqBin) > maxVoltageFreq) + { + maxVoltageFreq = newFreqSpectrum->GetAbs(iFreqBin); + } + } // end of freqeuncy bin loops + newAggFreqData.SetSummedGridVoltage(gridPointNumber, maxVoltageFreq); + } // End of grid + }// End of loop over all rings + KTDEBUG(agglog,"Channel summation performed over "<< fNRings<<" rings and "< -- relative phases of separate channels in radians [given as , , ...] + - "scales": vector -- scaling coefficients for separate channels, set individual channels to 0.0 if not included [given as , , ...] + + Slots: + - "fft": void (Nymph::KTDataPtr) -- Adds channels voltages using FFTW-phase information for appropriate phase addition; Requires KTFrequencySpectrumDataFFTW; Adds summation of the channel results; Emits signal "fft" + + Signals: + - "fft": void (Nymph::KTDataPtr) -- Emitted upon summation of all channels; Guarantees KTFrequencySpectrumDataFFTW + */ + + class KTPhasedAggregator : public Nymph::KTProcessor + { + public: + KTPhasedAggregator(const std::string& name = "phased-aggregator"); + virtual ~KTPhasedAggregator(); + + bool Configure(const scarab::param_node* node); + + // in meters, should not be hard-coded + MEMBERVARIABLE(double, ActiveRadius); + + // Get the grid size assuming a square grid + // Set default to 30 currently + MEMBERVARIABLE(unsigned, NGrid); + + // PTS:: for 18.6 keV electrons, this somehow needs to come from the data file or config file + MEMBERVARIABLE(double, Wavelength); + + //For exception handling to make sure the grid is defined before the spectra are assigned. + MEMBERVARIABLE(bool, IsGridDefined); + + // A boolean value to check whether the grid should be defined by the user + MEMBERVARIABLE(bool, IsUserDefinedGrid); + + // A boolean to run channel aggreagator when only a partial ring is defined + MEMBERVARIABLE(bool, IsPartialRing); + + // If IsPartialRing is true, the total number of channels is the product of the number of channels extracted from the egg file and the PartialRingMultiplicity + MEMBERVARIABLE(unsigned, PartialRingMultiplicity) + + // The text file to be used for the user-defined grid + MEMBERVARIABLE(std::string, UserDefinedGridFile); + + //The minimum frequency value above which the channel aggregated spectrum is calculated + MEMBERVARIABLE(double, SummationMinFreq); + + //The maximum frequency value below which the channel aggregated spectrum is calculated + MEMBERVARIABLE(double, SummationMaxFreq); + + // Number of axial rings/subarrays + MEMBERVARIABLE(unsigned, NRings); + + // Phase shifts and scale factors for complex signal summation + // variables read in as strings from config file and cast into vector + MEMBERVARIABLE(std::string, Phase_str); + MEMBERVARIABLE(std::string, Scale_str); + MEMBERVARIABLE(std::vector, Phases); + MEMBERVARIABLE(std::vector, Scales); + + //AN electron undergoiing cyclotron motion has a spiral motion and not all receving channels are in phase. + //If selected this option will make sure that there is a relative phase-shift applied + MEMBERVARIABLE(bool,UseAntiSpiralPhaseShifts); + + virtual bool SumChannelVoltageWithPhase(KTFrequencySpectrumDataFFTW& fftwData); + virtual bool SumChannelVoltageWithPhase(KTAxialAggregatedFrequencySpectrumDataFFTW& fftwData); + + protected: + + ///map that stores antispiral phase shifts + std::map fAntiSpiralPhaseShifts; + + /// Define the grid based (one per ring) + /*If fIsUserDefinedGrid it true, the grid is defined based on text file defined by fUserDefinedGridFile. + * If fIsUserDefinedGrid is false, defines a square grid with fNGrid*fNGrid points + * Returns total number of grid points defined + */ + unsigned DefineGrid(KTAggregatedFrequencySpectrumDataFFTW &newAggFreqData); + + /// Returns the phase shift based on a given point, angle of the channel and the wavelength + double GetPhaseShift(double xPosition, double yPosition, double wavelength, double channelAngle) const; + + /// Returns the phase shift based on a given point, angle of the channel and the wavelength + double GetAntiSpiralPhaseShift(double xPosition, double yPosition, double wavelength, double channelAngle) const; + + /// Get location of the point in the grid based on the given grid number and the size of the grid. + /* Returns true if the assigment went well, false if there was some mistake + */ + bool GetGridLocation(unsigned gridNumber, unsigned gridSize, double &gridLocation); + + /// Apply shift phase to the supplied points based on the phase provided + bool ApplyPhaseShift(double &realVal, double &imagVal, double phase); + + /// Generate antispiral phase shifts and save in fAntiSpiralPhaseShifts vector to be applied to channels + bool GenerateAntiSpiralPhaseShifts(unsigned channelCount); + + /// Convert frquency to wavlength + double ConvertFrequencyToWavelength(double frequency); + + virtual bool PerformPhaseSummation(KTFrequencySpectrumDataFFTWCore& fftwData,KTAggregatedFrequencySpectrumDataFFTW& newAggFreqData); + protected: + //PTS: This needs fixing, currently just setting each element to 0. But why does it have to be done to begin with. + // Perhaps there is some function in the utilities to do this ? + bool NullFreqSpectrum(KTFrequencySpectrumFFTW &freqSpectrum); + + //*************** + // Signals + //*************** + + protected: + Nymph::KTSignalData fSummedFrequencyData; + + //*************** + // Slots + //*************** + + protected: + Nymph::KTSlotDataOneType< KTFrequencySpectrumDataFFTW > fPhaseChFrequencySumSlot; + Nymph::KTSlotDataOneType< KTAxialAggregatedFrequencySpectrumDataFFTW> fAxialSumSlot; + }; + + inline bool KTPhasedAggregator::NullFreqSpectrum(KTFrequencySpectrumFFTW &freqSpectrum) + { + for (unsigned i = 0; i < freqSpectrum.size(); ++i) + { + freqSpectrum.SetRect(i, 0.0, 0.0); + } + return true; + } +} + +#endif /* KTPHASEDAGGREGATOR_HH_ */