RdMonitoring.cc

Go to the documentation of this file.

#include <cmath>
#include <complex>
#include <sstream>
#include <bitset>

#include "RdMonitoring.h"
#include <TROOT.h>

#include <TMinuit.h>

#include <TMath.h>
#include <cstdio>
#include <cstdlib>

#include <fwk/CentralConfig.h>
#include <fwk/LocalCoordinateSystem.h>
#include <fwk/CoordinateSystemRegistry.h>

#include <evt/Event.h>
#include <revt/Header.h>
#include <revt/EventTrigger.h>

#include <det/Detector.h>
#include <rdet/RDetector.h>
#include <rdet/Channel.h>
#include <rdet/Station.h>

#include <utl/Trace.h>
#include <utl/TimeStamp.h>
#include <utl/LeapSeconds.h>
#include <utl/TraceAlgorithm.h>
#include <utl/ErrorLogger.h>
#include <utl/Reader.h>
#include <utl/FFTDataContainerAlgorithm.h>
#include <utl/config.h>
#include <utl/AugerUnits.h>
#include <utl/ErrorLogger.h>
#include <utl/StringCompare.h>
#include <utl/UTCDateTime.h>

#include <utl/PhysicalConstants.h>
#include <utl/AxialVector.h>
#include <utl/Math.h>

#include <evt/Event.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerRRecData.h>
#include <revt/REvent.h>
#include <revt/Station.h>
#include <revt/Header.h>
#include <revt/StationRecData.h>

#include <CLHEP/Matrix/Matrix.h>
#include <CLHEP/Matrix/Vector.h>

using namespace std;
using namespace utl;
using namespace fwk;
using namespace revt;
using namespace evt;

namespace RdMonitoring {

  REvent* RdMonitoring::fgCurrentREvent;

  RdMonitoring::RdMonitoring() :
    fCurRunId(0),
    fCurEventId(0),
    fCurStationId(0),
    fCurChannelId(0),
    fCurGPSSecond(0),
    fPrevGPSSecond(0),
    fDownsample10sTriggers(true),
    fPathToStoreOutput(""),
    fMinFrequency(0),
    fMaxFrequency(0),
    fOutputFileName(""),
    fBeaconFrequencies(vector<double>()),
    fSpectrumNoiseFilterSize(31),
    fSpectrumPeakFilterSize(9),
    fSpectrumMinPeakSNR(0),
    fSpectrumFractionToSearch(0.01),
    fTraceNSigmaRMS(5)
  { }

  VModule::ResultFlag
  RdMonitoring::Init()
  {
    INFO("RdMonitoring::Init()");
    string tmpstring;
    
    Branch topBranch =
    CentralConfig::GetInstance()->GetTopBranch("RdMonitoring");
    topBranch.GetChild("Downsample10sTriggers").GetData(tmpstring);
    fDownsample10sTriggers = StringEquivalent(tmpstring, "yes");
    topBranch.GetChild("pathToStoreOutput").GetData(fPathToStoreOutput);
    topBranch.GetChild("minFrequency").GetData(fMinFrequency);
    topBranch.GetChild("maxFrequency").GetData(fMaxFrequency);
    topBranch.GetChild("OutputFileName").GetData(fOutputFileName);
    topBranch.GetChild("BeaconFrequencies").GetData(fBeaconFrequencies);

    topBranch.GetChild("SpectrumNoiseFilterSize").GetData(fSpectrumNoiseFilterSize);
    topBranch.GetChild("SpectrumPeakFilterSize").GetData(fSpectrumPeakFilterSize);
    topBranch.GetChild("SpectrumMinPeakSNR").GetData(fSpectrumMinPeakSNR);
    topBranch.GetChild("SpectrumFractionToSearch").GetData(fSpectrumFractionToSearch);

    topBranch.GetChild("TraceNSigmaRMS").GetData(fTraceNSigmaRMS);

    // read in the debug station list as vector, and then copy it to the set class variable
    vector<int> dummyStationIDlist;
    topBranch.GetChild("DebugStationIds").GetData(dummyStationIDlist);
    fDebugStationIds = set<int>(dummyStationIDlist.begin(), dummyStationIDlist.end());

    fMinFrequency *=1000; //to MHz
    fMaxFrequency *=1000; //to MHz

    // check xml settings for criteria
    if (fMinFrequency < 0 || fMinFrequency > 100) {
      WARNING("MinFrequency is negative or larger than 100 MHz! Set it to 0 MHz! \n");
      fMinFrequency = 0;
    }
    if (fMinFrequency >= fMaxFrequency) {
      WARNING("MinFrequency is not smaller than MaxFrequency! Set MinFrequency to 0 MHz and MaxFrequency to 100 MHz.\n");
      fMinFrequency = 0;
      fMaxFrequency = 100;
    }
    if (fMaxFrequency > 100) {
      WARNING("MaxFrequency is larger than 100 MHz! Set it to 100 MHz! \n");
      fMaxFrequency = 100;
    }

    ostringstream info;
    info << "Output files will be written as: \n" << fOutputFileName << '\n';
    INFO(info.str());

    std::ofstream mysqloutputfile;
    mysqloutputfile.open (fPathToStoreOutput+"mysql_"+fOutputFileName+".txt",std::ofstream::out | std::ofstream::trunc);
    mysqloutputfile.close();
    
    return eSuccess;
  }

  VModule::ResultFlag
  RdMonitoring::Run(evt::Event& event)
  {
    if (!event.HasREvent()) {
      ERROR("Read event without radio part");
      return eFailure;
    }

    // General event properties
    REvent& rEvent  = event.GetREvent();
    fgCurrentREvent = &event.GetREvent();
    fCurRunId       = rEvent.GetHeader().GetRunNumber();
    fCurEventId     = rEvent.GetHeader().GetId();
    fCurGPSSecond   = rEvent.GetHeader().GetTime().GetGPSSecond();
    //

    // following condition also allows for 5s inaccuracy in the timing of the 100 second trigger
    // If you want to use all data then set Downsample10sTriggers to false in the XML.
    if (fCurGPSSecond-fPrevGPSSecond > 95 || fDownsample10sTriggers==false) {
      
      // check if beacon frequencies have to be read from database, or are provided in XML file
      // in the first case, fBeaconFrequencies is empty, i.e. size()==0  
      // get Detector
      const det::Detector& Det = det::Detector::GetInstance();
      const rdet::RDetector& RDetector = Det.GetRDetector();

      std::vector<double> beaconFrequencies(fBeaconFrequencies);
      if(beaconFrequencies.size() < 1) {
        try {
          beaconFrequencies = RDetector.GetBeaconFrequencies(); // for example for 2013 returns: 0.058887 0.071191 0.068555 0.061523
        } catch (utl::DataNotFoundInDBException& err) {
          ERROR ("Error reading Beacon frequencies from database");
          return eFailure;
        }  
      }

      if (beaconFrequencies.size() > 100) {
        cout << "Too many active beacon frequencies extracted from the database!? Exiting... " << endl;
        return eFailure;
      }
      if (beaconFrequencies.size() < 1) {
        cout << "Can not read any active beacon frequencies from the database! Exiting ..." << endl;
        return eFailure;
      }
      //

      // Loop over stations
      for (REvent::StationIterator sIt = rEvent.StationsBegin(); sIt != rEvent.StationsEnd(); ++sIt) {
        fCurStationId     = sIt->GetId();

        // Loop over channels
        for (Station::ChannelIterator cIt = sIt->ChannelsBegin(); cIt != sIt->ChannelsEnd(); ++cIt) {
          if (!cIt->IsActive()) {
            continue;
          }
    
          if (rEvent.GetTrigger().IsPeriodicTrigger() ) {
            fCurChannelId=cIt->GetId();

            // Start calculating properties
            const revt::ChannelTimeSeries trace           = cIt->GetChannelTimeSeries();
            const revt::ChannelFrequencySpectrum spectrum = cIt->GetChannelFrequencySpectrum();


            // TRACE:
            //    * Mean
            //    * Mean-subtracted RMS
            //    * Min
            //    * Max
            double traceMean= TraceAlgorithm::Mean(trace, 0, trace.GetSize()-1) ;
            double traceRMS = TraceAlgorithm::StandardDeviation(trace, 0, trace.GetSize()-1) ;  //StandardDeviation() subtracts mean. RootMeanSquare() doesn't
            //double traceRMS = TraceAlgorithm::RootMeanSquare(trace, 0, trace.GetSize()-1) ;
            double traceMax = TraceAlgorithm::Max(trace, 0, trace.GetSize()-1)            ;
            double traceMin = TraceAlgorithm::Min(trace, 0, trace.GetSize()-1)            ;
            //


            for (set<int>::iterator it=fDebugStationIds.begin(); it != fDebugStationIds.end(); ++it) {
              if (*it==fCurStationId){
                std::ofstream mysqltraceoutputfile;
                mysqltraceoutputfile.open (fPathToStoreOutput+"trace_r"+std::to_string(fCurRunId)+"_e"+std::to_string(fCurEventId)+"_st"+std::to_string(fCurStationId)+"_ch"+std::to_string(fCurChannelId)+"_t"+std::to_string(fCurGPSSecond)+".txt",std::ios_base::app);
                cout << "\tStoring trace for   "
                                 << "\tst "       << fCurStationId   
                                 << "\tch "       << fCurChannelId   
                                 << "\tGPSSec="   << fCurGPSSecond    
                                 << "\tN_samples="<< trace.GetSize() << endl;
                float timeBinning = trace.GetBinning();
                for (unsigned int j = trace.GetStart(); j < trace.GetSize(); ++j) {
                   mysqltraceoutputfile << j*timeBinning<< "\t" << trace[j] << endl;
                }
                mysqltraceoutputfile.close();
              } // end if
            } // end for

            // SPECTRUM:
            //    * Power at [28,33,38,43,51,58,69,80] MHz: corresponding to the frequencies that drop in amplitude when a dipole antenna breaks
            //    * PowerInBand      : sum(spectrum amplitudes in band)*binsize (range set in xml file)
            //    * PowerOutsideBand
            //    * TraceOutliers
            //    * SpectrumNarrowBandPeaks (calculated in beacon section below)
            
            double PowerInBand      = 0;
            double PowerOutsideBand = 0;
            list<float> freqList    = {28,33,38,43,51,58,69,80}; //In MHz
            float PowerAt[8]        = { 0. };
            int i_freqList          = 0;
            
            for (float f : freqList) {
              PowerAt[i_freqList] = abs(spectrum[  round((f-cIt->GetFrequencyOfBin(0)/megahertz)/(spectrum.GetBinning()*1000))    ]);
              i_freqList++;
            }

            for (set<int>::iterator it=fDebugStationIds.begin(); it != fDebugStationIds.end(); ++it) {
              if (*it==fCurStationId){
                std::ofstream mysqlspectrumoutputfile;
                mysqlspectrumoutputfile.open (fPathToStoreOutput+"spectrum_r"+std::to_string(fCurRunId)+"_e"+std::to_string(fCurEventId)+"_st"+std::to_string(fCurStationId)+"_ch"+std::to_string(fCurChannelId)+"_t"+std::to_string(fCurGPSSecond)+".txt",std::ios_base::app);
                cout << "\tStoring spectrum for: "
                                 << "\tst "       << fCurStationId   
                                 << "\tch "       << fCurChannelId   
                                 << "\tGPSSec="   << fCurGPSSecond    
                                 << "\tN_samples="<< spectrum.GetSize() << endl;
                for (unsigned int j = spectrum.GetStart(); j < spectrum.GetSize(); ++j) {
                   mysqlspectrumoutputfile << cIt->GetFrequencyOfBin(j)/megahertz << "\t" << abs(spectrum[j]) << endl;
                }
                mysqlspectrumoutputfile.close();
              } //end if
            } //end for

            for (unsigned int j = spectrum.GetStart()+1; j < spectrum.GetSize(); ++j) {
              if (cIt->GetFrequencyOfBin(j)/megahertz > fMinFrequency && cIt->GetFrequencyOfBin(j)/megahertz < fMaxFrequency){
                PowerInBand +=  abs(spectrum[j]);
              }
              else {
                PowerOutsideBand +=  abs(spectrum[j]);
              }
            }
            PowerOutsideBand*=spectrum.GetBinning()*1000;   // PowerInBand[ADC*MHz]=SUM(spectrum amplitudes)*binsize
            PowerInBand     *=spectrum.GetBinning()*1000;
            //
         
            
            // BEACON
            //    * BeaconDetected        
            //    * BeaconSNR:        set noise window size around beacon in xml
            //    * (SpectrumNarrowBandPeaks: (used on channel level) amount of narrow-band peaks in the spectrum (excluding identified beacon lines))

            float f_step        = spectrum.GetBinning()*1000;                // in Mhz
            
            Trace<double> amplitudeSpectrum;
            for (ChannelFrequencySpectrum::SizeType i = 0; i < spectrum.GetSize(); ++i)
              amplitudeSpectrum.PushBack(abs(spectrum[i]));

            // WRITE BEACON INFO TO OUTPUT
            //

            std::ofstream mysqloutputfile;
            mysqloutputfile.open (fPathToStoreOutput+"mysql_"+fOutputFileName+".txt",std::ios_base::app);
            
            //
            // CONTINUES IN NEXT PART

            // RFI peak & beacon finding
            //  1) take highest value in spectrum --> peak candidate. 
            //  2) calculate snr around peak. 
            //  3) add bins around peak to SpectrumRejectedPositions to account for width of this already found peak so they are not identified on each bin. 
            //  4) take next highest value and repeat untill all highest values exhausted


            unsigned int SpectrumSize=spectrum.GetSize()-1; // max size is spectrum.GetSize()-1
            unsigned int SpectrumBinsToSearch=ceil(SpectrumSize*fSpectrumFractionToSearch);
            
            // create simple vectors of trace and spectum
            std::vector<float> UnsortedSpectrum;
            for (ChannelFrequencySpectrum::SizeType i = 0; i < SpectrumSize; ++i){
              UnsortedSpectrum.push_back(abs(spectrum[i+1]));   //skipping bin one as it is not part of the spectrum values.
            }

            std::vector<int> SpectrumSortedIndices;
            for(size_t i = 0; i < SpectrumSize; ++i) {
              SpectrumSortedIndices.push_back(i);   //indices will still start at zero!
            }
            std::sort(SpectrumSortedIndices.begin(),SpectrumSortedIndices.end(),   [&](int i1, int i2) { return UnsortedSpectrum[i1] > UnsortedSpectrum[i2]; }   );

            std::vector<int> _SpectrumPeakPos;
            std::vector<float> _SpectrumPeakHeight;
            for(size_t i = 0; i < SpectrumBinsToSearch; ++i){
              _SpectrumPeakPos.push_back(SpectrumSortedIndices[i]);
              _SpectrumPeakHeight.push_back(UnsortedSpectrum[SpectrumSortedIndices[i]]);
            }


            int SpectrumNarrowBandPeaks=0;

            std::string tempString="";
            std::vector<int> SpectrumRejectedPos;

            for(size_t i = 0; i < SpectrumBinsToSearch; ++i){
              if (  _SpectrumPeakPos[i]> trunc(fSpectrumNoiseFilterSize/2.0) && _SpectrumPeakPos[i]< SpectrumSize -ceil(fSpectrumNoiseFilterSize/2.0)   ) {
                if (( std::find(SpectrumRejectedPos.begin(), SpectrumRejectedPos.end(), _SpectrumPeakPos[i]) != SpectrumRejectedPos.end() )==false){
                  double noise = TraceAlgorithm::Median(amplitudeSpectrum,
                                1+_SpectrumPeakPos[i]-trunc(fSpectrumNoiseFilterSize/2.0),
                                1+_SpectrumPeakPos[i]+ceil(fSpectrumNoiseFilterSize/2.0)-1.); //amplitudeSpectrum starts at bin 0. for PeakPos bin 0 was ignored before. needs compensation now.
                  double snr  = (_SpectrumPeakHeight[i]*_SpectrumPeakHeight[i])/(noise*noise) ;
                  
                  if (snr>fSpectrumMinPeakSNR){
                    SpectrumNarrowBandPeaks++;
                  }

                  double CurrentFrequency=cIt->GetFrequencyOfBin(_SpectrumPeakPos[i]+1)*1000;

                  for (unsigned int BqIt = 0; BqIt < beaconFrequencies.size(); ++BqIt) {
                    if (TMath::Abs(CurrentFrequency - beaconFrequencies[BqIt]*1000.) < f_step) {
                      if (snr>fSpectrumMinPeakSNR){
                        SpectrumNarrowBandPeaks-=1; // do not count as NarrowBandPeak if it is a beacon
                      }
                    tempString+="INSERT INTO RdMonitoringBeacon (RdMonitoringChannel_id,BeaconFreq,BeaconSNR) VALUES (LAST_INSERT_ID(),"+std::to_string(beaconFrequencies[BqIt]*1000)+","+std::to_string(snr)+");\n";
                    break;
                    }
                  }

                  for (size_t j=_SpectrumPeakPos[i]-trunc(fSpectrumPeakFilterSize/2.0); j < _SpectrumPeakPos[i]+ceil(fSpectrumPeakFilterSize/2.0) ; j++){
                    SpectrumRejectedPos.push_back(j);
                  }
                
                }// end if non-rejected
              }// end if not near edge
            }// end peak candidate loop
            
            int TraceOutliers=0;
            for (unsigned int j = trace.GetStart(); j < trace.GetSize(); ++j) {
               if (TMath::Abs(trace[j]-traceMean) > traceRMS*fTraceNSigmaRMS) {
                  TraceOutliers++;
               }
            }
            //timedomainoutputfile.close();

            //
            //

            
            // WRITE CHANNEL INFO TO OUTPUT
            mysqloutputfile << "INSERT INTO RdMonitoringChannel (" <<
                                    "RunId, "         <<
                                    "StationId, "     <<
                                    "ChannelId, "     <<
                                    "GPSSecond, "     <<
                                    "traceMean, "      <<
                                    "traceRMS, "      <<
                                    "traceMin, "      <<
                                    "traceMax, "      <<
                                    "PowerAt28MHz, "  <<
                                    "PowerAt33MHz, "  <<
                                    "PowerAt38MHz, "  <<
                                    "PowerAt43MHz, "  <<
                                    "PowerAt51MHz, "  <<
                                    "PowerAt58MHz, "  <<
                                    "PowerAt69MHz, "  <<
                                    "PowerAt80MHz,"   <<
                                    "PowerInBand,"    <<
                                    "PowerOutsideband," <<
                                    "TraceOutliers,"  <<
                                    "SpectrumNarrowBandPeaks,"  <<
                                    "BadPeriodFlag"     <<                //initial flag=0: in post-processing the flag is set
                                    ") Values (" <<
                                        fCurRunId         << ","<< 
                                        fCurStationId     << "," <<
                                        fCurChannelId     << "," <<
                                        fCurGPSSecond     << "," <<
                                        traceMean          << "," <<
                                        traceRMS          << "," <<
                                        traceMin          << "," <<
                                        traceMax          << "," <<
                                        PowerAt[0]        << "," <<
                                        PowerAt[1]        << "," <<
                                        PowerAt[2]        << "," <<
                                        PowerAt[3]        << "," <<
                                        PowerAt[4]        << "," <<
                                        PowerAt[5]        << "," <<
                                        PowerAt[6]        << "," <<
                                        PowerAt[7]        << "," <<
                                        PowerInBand       << "," <<
                                        PowerOutsideBand  << "," <<
                                        TraceOutliers   << "," <<
                                        SpectrumNarrowBandPeaks   << "," <<
                                        0                 << ");\n"; 
            mysqloutputfile << tempString.c_str() ;
            mysqloutputfile.close();
            //


          } // end if periodic
        } // end channel loop        
      } // end station loop

      fPrevGPSSecond = fCurGPSSecond;
    } // end if periodic event
    return eSuccess;
  } // end event


  VModule::ResultFlag
  RdMonitoring::Finish()
  {
    INFO("Finish");
    return eSuccess;
  }

}