RdGalacticDatasetMaker.cc

Go to the documentation of this file.

#include "RdGalacticDatasetMaker.h"

// framework
#include <fwk/CentralConfig.h>

// utilities
#include <utl/ErrorLogger.h>
#include <utl/TimeStamp.h>
#include <utl/TraceAlgorithm.h>
#include <utl/UTCDateTime.h>
#include <utl/Minou.h>
// event
#include <evt/Event.h>

// revent
#include <revt/REvent.h>
#include <revt/Header.h>
#include <revt/Station.h>

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

// general
#include <stdio.h>
#include <stdlib.h>
#include <iomanip>

// for folder creation
#include <string.h>
#include <errno.h>
#include <sys/stat.h> 
#include <sys/types.h> 

// boost
#include <boost/algorithm/string.hpp>

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


namespace RdGalacticDatasetMaker {

  REvent* RdGalacticDatasetMaker::fgCurrentREvent;
  
  RdGalacticDatasetMaker::RdGalacticDatasetMaker() :
    fInfoLevel(0),
    fCurRunId(0),
    fCurEventId(0),
    fCurStationId(0),
    fCurChannelId(0),
    fADCthreshold(0),
    fSTDthreshold(0),
    ffreqBinWidth(0),
    ftimeBinWidth(0),
    fdebugData(0),
    fdayStack(false),
    ffilteredDebugData(false),
    fchannels(vector<int>()),
    fMinFrequency(0),
    fMaxFrequency(0),
    fstartBin(0),
    ftraceCutBinLength(0),
    fsamplingFrequency4Header(0),
    fimpedanceFile(false)
    {}

  VModule::ResultFlag
  RdGalacticDatasetMaker::Init()
  {
    INFO("RdGalacticDatasetMaker::Init()");
    string tmpstring;
    ostringstream info;

    Branch topBranch =
    CentralConfig::GetInstance()->GetTopBranch("RdGalacticDatasetMaker");
    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);
    topBranch.GetChild("freqBinWidth").GetData(ffreqBinWidth);
    topBranch.GetChild("timeBinWidth").GetData(ftimeBinWidth);
    topBranch.GetChild("ADCthreshold").GetData(fADCthreshold);
    topBranch.GetChild("STDthreshold").GetData(fSTDthreshold);
    topBranch.GetChild("debugData").GetData(fdebugData);
    string fdayStackStr = "";
    topBranch.GetChild("dayStack").GetData(fdayStackStr);
    fdayStack = (fdayStackStr=="yes");
    string filteredDebugDataStr = "";
    topBranch.GetChild("filteredDebugData").GetData(filteredDebugDataStr);
    ffilteredDebugData = (filteredDebugDataStr=="yes");
    topBranch.GetChild("channels").GetData(fchannels);
    topBranch.GetChild("MinFrequency").GetData(fMinFrequency);
    topBranch.GetChild("MaxFrequency").GetData(fMaxFrequency);
    topBranch.GetChild("startBin").GetData(fstartBin);
    topBranch.GetChild("traceCutBinLength").GetData(ftraceCutBinLength);
    topBranch.GetChild("samplingFrequency").GetData(fsamplingFrequency4Header);
    string fimpedanceFileStr = "";
    topBranch.GetChild("impedanceFile").GetData(fimpedanceFileStr);
    fimpedanceFile = (fimpedanceFileStr=="yes");
    
    vector<int> StationIDlist_dummy;
    Branch stationsB = topBranch.GetChild("stations");
    for (Branch stationB = stationsB.GetFirstChild(); stationB; stationB = stationB.GetNextSibling()) {
      StationIDlist_dummy.push_back(stoi(stationB.GetAttributes().find("id")->second));
    }
    
    fStationIds = vector<int>(StationIDlist_dummy.begin(), StationIDlist_dummy.end());
    
    for (unsigned int i = 0; i < fStationIds.size(); ++i) {
      info << "Station: " << fStationIds[i] << endl;
    }
    INFODebug(info);
    info.str("");
    
    if (fimpedanceFile == true) {
      (impedanceTopBranch = topBranch.GetChild("RImpedanceProfile"));
    } else {
      WARNING( "No impedance file selected. Impedance will be set to 1 Ohms for all frequencies.");
    }
  
    info << "Settings: Frequency bin width: " << ffreqBinWidth << "MHz" << "; Time bin width: "
    << ftimeBinWidth << " hours" << "; ADC threshold: " << fADCthreshold << " ADC" << "; STD threshold: "
    << fSTDthreshold << " ADC" << "; DayStack is " << fdayStack << "; MinFrequency: " << fMinFrequency << " MHz"
    << " MaxFrequency: " << fMaxFrequency << " MHz" << "; Start bin: " << fstartBin << "; Total bins to take from time trace: " 
    << ftraceCutBinLength << "; Sampling frequency (for the header in the output files: "
    << fsamplingFrequency4Header << endl;
    INFOFinal(info.str());
    info.str("");
  
    // create time and frequency bin vector based on bin widths
    // create time bin vector
    for (float currentTimeBin = 0; currentTimeBin <= 24; currentTimeBin += ftimeBinWidth){
      timeBins.push_back(currentTimeBin);
    }
    // create freq bin vector
    for (float currentFreqBin = fMinFrequency; currentFreqBin <= fMaxFrequency; currentFreqBin += ffreqBinWidth){
      freqBins.push_back(currentFreqBin);
    }
    
    // init dataset vectors
    vector<vector<vector<vector<float>>>> AllChannelsAndStationsDatasets_dummy(fchannels.size(), vector<vector<vector<float>>>
    (fStationIds.size(), vector<vector<float>>(timeBins.size()-1, vector<float>(freqBins.size()+1)) ));
    
    AllChannelsAndStationsDatasets = AllChannelsAndStationsDatasets_dummy;
    vector<vector<float> > stationTraceCounter_dummy(fStationIds.size(), vector<float>(2,0));
    stationTraceCounter = stationTraceCounter_dummy;

    vector<vector<vector<vector<float>>>> (0, vector<vector<vector<float>>>
    (0, vector<vector<float>>(0, vector<float>(0)) )).swap(AllChannelsAndStationsDatasets_dummy); // deallocating memory
    
    // fill the Time columns
    for (unsigned int chan = 0; chan < AllChannelsAndStationsDatasets.size(); ++chan) {
      for (unsigned int st = 0; st < AllChannelsAndStationsDatasets[chan].size(); ++st) {
        for (unsigned int row = 0; row < AllChannelsAndStationsDatasets[chan][st].size(); ++row) {
          AllChannelsAndStationsDatasets[chan][st][row][0] = timeBins[row];
        }
      } //end station
    } //end chan

    info << "Dataset vectors with dimensions " << AllChannelsAndStationsDatasets.size() << "x"
    << AllChannelsAndStationsDatasets[0].size() << "x" << AllChannelsAndStationsDatasets[0][0].size() 
    << "x" << AllChannelsAndStationsDatasets[0][0][0].size() << "  (channels x stations x Time bins x frequency bins) has been created.";
    INFOFinal(info.str());
    info.str("");
    
    // Create dataset directory 
    if (mkdir("datasets", 0777) == -1){
        WARNING("Directory for datasets cannot be created, maybe it exists already?");
      //  cerr << "Error, directory for datasets :  " << strerror(errno) << endl;
    } else {
      info << "Directory for datasets has been created";
      INFO(info.str());
    } 
    info.str("");

    // Create debug directory 
    if (fdebugData >= 1){
      if (mkdir("debug", 0777) == -1){
        WARNING("Directory for debug data cannot be created, maybe it exists already?");
    //    cerr << "Error, directory for debug : " << strerror(errno) << endl; 
      } else {
        info << "Directory for debug has been created";
        INFO(info.str());
      }
    }
    info.str("");

    for (unsigned int i=0; i<fStationIds.size(); ++i) {
        // initialize save to file
      for (unsigned int j=0; j< fchannels.size(); ++j) {
        if (fdebugData >= 1){
          
          ofstream tracePropertiesDebugFile; // initialize std and peakmax spectrum trace to file
          tracePropertiesDebugFile.open 
          ("./debug/tracePropertiesDebugFile_"+to_string(fStationIds[i])+"_channel_"+to_string(fchannels[j])+".txt", ios::out | ios::app);
          
          tracePropertiesDebugFile  << "event" << "  " << "station" <<"  " << "channel" << "  gpsTime" 
          << "  UTChour" << "  sidTime"  << "  STD" << "  absMaxPeakADC" <<endl;
          tracePropertiesDebugFile.close();
          
        } // end  fdebugData >=1
        if (fdebugData == 2){
          
          ofstream ADCtraceFile;
          ADCtraceFile.open 
          ("./debug/ADCtracesDebugFile_"+to_string(fStationIds[i])+"_channel_"+to_string(fchannels[j])+".txt", ios::out | ios::app);
          
          ADCtraceFile  << "event" << "  " << "station" <<"  " << "channel" << "  gpsTime" << "  UTChour" << "  sidTime" ;
          for (int k=fstartBin; k<ftraceCutBinLength; ++k) {
            ADCtraceFile << "  " << k;
          }
          ADCtraceFile << endl;
          ADCtraceFile.close();
          
          // initialize save  to file
          ofstream spectrumTraceFile;
          spectrumTraceFile.open 
          ("./debug/spectraDebugFile_"+to_string(fStationIds[i])+"_channel_"+to_string(fchannels[j])+".txt", ios::out | ios::app);
          
          spectrumTraceFile  << "event" << "  " << "station" <<"  " << "channel" << "  gpsTime" << "  UTChour" << "  sidTime";
          for (int k=fstartBin; k<(ftraceCutBinLength/2+1); ++k) {
            spectrumTraceFile << "  " << k*fsamplingFrequency4Header/ftraceCutBinLength;
          }
          spectrumTraceFile << endl;
          spectrumTraceFile.close();
          
        } // end  fdebugData == 2
      } // end channel file creation loop
    // end of station file creation loop to file
    }
    return eSuccess;
  } //end init


  VModule::ResultFlag
  RdGalacticDatasetMaker::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();
    fHHMMSS = TimeStamp2HHMMSS(fCurGPSSecond);
    fCurUTChour = HHMMSS2hour(fHHMMSS);

    // Loop over stations
    for (REvent::StationIterator sIt = rEvent.StationsBegin(); sIt != rEvent.StationsEnd(); ++sIt) {
      fCurStationId = sIt->GetId();
      Station& station = *sIt;
      
      // init filter flag
      bool filterPass = true;

      // init off all temp vectors
      vector<vector<float>> powerDepositedInFreqBins(fchannels.size(), vector<float>(freqBins.size()-1, 0.0));
      vector<vector<float>> cachedSpectra(fchannels.size(), vector<float>());
      vector<vector<float>> cachedADCtraces(fchannels.size(), vector<float>());
      vector<vector<float>> cachedSTD(fchannels.size(), vector<float>());
      vector<vector<float>> cachedTraceADCmax(fchannels.size(), vector<float>());
        
      // Loop over channels
      const rdet::Station& detStation = det::Detector::GetInstance().GetRDetector().GetStation(station); // rdStation
      for (Station::ChannelIterator cIt = station.ChannelsBegin(); cIt != station.ChannelsEnd(); ++cIt) {
        Channel& channel = *cIt;
        
        // Check if channel exists
        if (!channel.IsActive()) {
          continue;
        }

        ostringstream info;

        fCurChannelId=channel.GetId();
        int channelIndex = getChannelIndex(fchannels,fCurChannelId); //find the channel index in the selected channels
        const rdet::Channel& rdetChannel = detStation.GetChannel(fCurChannelId);

        // Start
        

        // check length before
          
        // get time traces
        revt::ChannelADCTimeSeries ADCtrace = cIt->GetChannelADCTimeSeries();
        revt::ChannelTimeSeries& trace      = cIt->GetChannelTimeSeries();
        
        // time length of the trace
        float T = trace.GetSize()*trace.GetBinning();

        info << "BEFORE cutting: trace size: " << trace.GetSize() << ";  ADC trace size: " << ADCtrace.GetSize()
        << "; Spectrum should be half+1 of that" << "; The length of the time trace is: "<< T << "ns" <<endl; 
        INFODebug(info);
        info.str("");
        
        // if different length than default, do the cut
        if (ftraceCutBinLength !=2048) {
          revt::ChannelTimeSeries tempTimeSeries;
          revt::ChannelADCTimeSeries tempADCTimeSeries;

          tempTimeSeries.SetBinning(trace.GetBinning());
          tempADCTimeSeries.SetBinning(ADCtrace.GetBinning());
            
          unsigned int endBin = fstartBin + ftraceCutBinLength;
          for (unsigned int i = fstartBin; i < endBin; ++i) {
            tempTimeSeries.PushBack(trace[i]);
            tempADCTimeSeries.PushBack(ADCtrace[i]);
          }
          trace = tempTimeSeries;
          ADCtrace = tempADCTimeSeries;
          // request new instance
          trace  = cIt->GetChannelTimeSeries();
          // the length of the time trace should change accordingly
          T = trace.GetSize()*trace.GetBinning();
        }
        
        // get spectrum
        const revt::ChannelFrequencySpectrum spectrum = channel.GetChannelFrequencySpectrum();
        // get sampling frequency
        const double samplingfreq = rdetChannel.GetSamplingFrequency()*1000;
        
        info << "AFTER cutting: trace size: " << trace.GetSize() << "; ADC trace size: " << ADCtrace.GetSize()
        << "; Spectrum size: " << spectrum.GetSize() << "; The length of the time trace is: "<< T << "ns" <<endl; 
        INFODebug(info);
        info.str("");
        
        info << "Trace size: " << trace.GetSize() << "; Trace bin length (ns): " <<  trace.GetBinning() 
        << "Spectrum bin size (GHz): " <<  spectrum.GetBinning() << "; Trace time length (nano sec): " << T << endl;
        info << "GPS time:  " << fCurGPSSecond << ";  LST time: " << GPSSecToLSTHourAuger(fCurGPSSecond)  
        << ";  HHMMSS: " << fHHMMSS << ";  Converted to UTC hour: "  << fCurUTChour  << ";  Spectrum.GetSize(): " 
        << spectrum.GetSize()  << ";  Sampling rate from XML: " << samplingfreq << "MHz" << endl;
        INFOIntermediate(info);
        info.str("");
          
        // get statistics
        int ADCtraceAbsMax;
        double ADCtraceSTD = TraceAlgorithm::StandardDeviation(ADCtrace, 0, trace.GetSize()-1);
        double ADCtraceMax = TraceAlgorithm::Max(ADCtrace, 0, trace.GetSize()-1);
        double ADCtraceMin = TraceAlgorithm::Min(ADCtrace, 0, trace.GetSize()-1);
        
        if (ADCtraceMax >= abs(ADCtraceMin)){
          ADCtraceAbsMax = ADCtraceMax;
        }
        if (ADCtraceMax < abs(ADCtraceMin)){
          ADCtraceAbsMax = abs(ADCtraceMin);
        }
        
        // update filter flag
        if (filterPass == true){
          if ( (ADCtraceSTD < fSTDthreshold) & (ADCtraceAbsMax < fADCthreshold) ){
            filterPass = true;
          } else {
            filterPass = false;
          }
        }

        info << "Threshold settings and quantity values in the Channel " << to_string(fCurChannelId)
        << ";  STD threshold: " << fSTDthreshold << ";  ADC threshold: " << fADCthreshold << "\n" 
        << "Properties of current trace::: abs Max: " << ADCtraceAbsMax <<";  STD: " << ADCtraceSTD 
        << ";  Min: " << ADCtraceMin << ";  Max" << ADCtraceMax << ";  Filterpass flag: " << filterPass << endl;
        INFOIntermediate(info);
        info.str("");
        
        //bin integration
        float Z;
        for (unsigned j=0; j < freqBins.size()-1; ++j){
          for (unsigned int i = 0; i < spectrum.GetSize(); ++i){
            if ((channel.GetFrequencyOfBin(i)*1000 >= freqBins[j]) && (channel.GetFrequencyOfBin(i)*1000 < freqBins[j+1] )){
              if (fimpedanceFile == false){
                Z = 1;
                WARNING("No impedance from file! Setting Z = 1 for all frequencies.");
              } else {
                Z = GetImpedanceAt(impedanceTopBranch, fCurChannelId, channel.GetFrequencyOfBin(i)*1000);
              }
              
              info << "Interpolated impedance for frequency " << to_string(channel.GetFrequencyOfBin(i)*1000) 
              << "in channel " << to_string(fCurChannelId) << " is " << to_string(Z) << " This was done in bin interval: " 
              << freqBins[j] << " to " << freqBins[j+1] << "[MHz]" << endl;
              INFODebug(info);
              info.str("");
              
              // power integration, the units in T, spectrum.GetBinning and spectrum*spectrum cancel each other
              powerDepositedInFreqBins[channelIndex][j] += 2*(1/Z)*(1/T)*abs(spectrum[i]*spectrum[i])*spectrum.GetBinning();
            
            }
          }
        } // end bin integration
  
        // fill debug vectors
        if (fdebugData >= 1) {
          cachedSTD[channelIndex].push_back(ADCtraceSTD); 
          cachedTraceADCmax[channelIndex].push_back(ADCtraceAbsMax);
        } // end of debugData >= 1
        if (fdebugData == 2) {
          for (unsigned int j = spectrum.GetStart(); j < spectrum.GetSize(); ++j) {
            cachedSpectra[channelIndex].push_back(abs(spectrum[j]));
          }
          for (unsigned int j = ADCtrace.GetStart(); j < ADCtrace.GetSize(); ++j) {
            cachedADCtraces[channelIndex].push_back(ADCtrace[j]);
          }
        } // end of debugData == 2
      } // end channel loop 

      // flag for if dump filtered debug data or not
      bool filterPass4debugData = true;
      if (ffilteredDebugData == true){
        if (filterPass == true){
          filterPass4debugData = true;
        } else {
          filterPass4debugData = false;
        }
      }
          
      // debug data dump
      if (filterPass4debugData == true) {
        for (unsigned int channelIndex=0; channelIndex < fchannels.size(); ++channelIndex) {
          if (fdebugData >= 1){
          // debuggin tool, save all STD and ADC peak height with event,station and channel number to to file
            ofstream tracePropertiesDebugFile;
            tracePropertiesDebugFile.open 
            ("./debug/tracePropertiesDebugFile_"+to_string(fCurStationId)+"_channel_"+to_string(fchannels[channelIndex])+".txt", ios::out | ios::app);
            
            tracePropertiesDebugFile  << fCurEventId << "  " << fCurStationId <<"  " << fCurChannelId 
            << "  " << fCurGPSSecond << "  " << fCurUTChour << "  " << GPSSecToLSTHourAuger(fCurGPSSecond) 
            << "  " << cachedSTD[channelIndex][0] << "  " << cachedTraceADCmax[channelIndex][0] << endl;
            tracePropertiesDebugFile.close();
            
          } //end if debug >= 1
          if (fdebugData == 2){
          // debuggin tool, save all spectra with event,station and channel number to to file
            ofstream spectrumTraceFile;
            spectrumTraceFile.open
            ("./debug/spectraDebugFile_"+to_string(fCurStationId)+"_channel_"+to_string(fchannels[channelIndex])+".txt", ios::out | ios::app);
            
            spectrumTraceFile  << fCurEventId << "  " << fCurStationId <<"  " << fCurChannelId << "  " 
            << fCurGPSSecond << "  " << fCurUTChour << "  " << GPSSecToLSTHourAuger(fCurGPSSecond);
            for (unsigned int j = 0; j < cachedSpectra[channelIndex].size(); ++j) {
               spectrumTraceFile <<  "  " << cachedSpectra[channelIndex][j];
            }
            spectrumTraceFile << "\n";
            spectrumTraceFile.close();
            // end of saving to file
          
            // debuggin tool, save all ADC time traces with event,station and channel number to to file
            ofstream ADCtraceFile;
            ADCtraceFile.open 
            ("./debug/ADCtracesDebugFile_"+to_string(fCurStationId)+"_channel_"+to_string(fchannels[channelIndex])+".txt", ios::out | ios::app);
            
            ADCtraceFile  << fCurEventId << "  " << fCurStationId <<"  "<< fCurChannelId << "  " 
            << fCurGPSSecond << "  " << fCurUTChour  << "  " << GPSSecToLSTHourAuger(fCurGPSSecond);
            for (unsigned int j = 0; j < cachedADCtraces[channelIndex].size(); j++) {
               ADCtraceFile <<  "  " << cachedADCtraces[channelIndex][j] ;
            }
            ADCtraceFile << "\n";
            ADCtraceFile.close();
            
          } //end if debug == 2
        } // end channel loop
      } // end if pass

      
      // put the vector to the correct sidereal time if the trace passed the threshold filters
      if (filterPass == true){
        float currentTime;
        if (fdayStack == true){
          currentTime = fCurUTChour; // "normal" hour time
        } else {
          currentTime = GPSSecToLSTHourAuger(fCurGPSSecond); // LST hour time
        }
        
        int n;
        for (unsigned int chan = 0; chan < AllChannelsAndStationsDatasets.size(); ++chan) {
          for (unsigned int st = 0; st < AllChannelsAndStationsDatasets[0].size(); ++st) {
            if (fCurStationId == fStationIds[st]){
              for (unsigned int row = 0; row < timeBins.size()-1; ++row) {
                if ((currentTime >= timeBins[row]) & (currentTime < timeBins[row+1] )) {
                  AllChannelsAndStationsDatasets[chan][st][row][1] += 1; // this is the weight column
                  n = AllChannelsAndStationsDatasets[chan][st][row][1]; // it easier to pass this to the function "online average"
                  
                  // the row is starting from 2nd columns, since the zeroth is lst time and the first is the weight of the row
                  transform(AllChannelsAndStationsDatasets[chan][st][row].begin()+2, AllChannelsAndStationsDatasets[chan][st][row].end(), 
                  powerDepositedInFreqBins[chan].begin(), AllChannelsAndStationsDatasets[chan][st][row].begin()+2, [n](float v1, float v2)
                  {return onlineMean(v1,v2,n);});
                  
                } // end if is in Time bin
              } // end row
            } // end if station
          } // end station loop
        } // end channel loop
      } // end filterPass
      
      // just some statistics, to see at the end how many traces in each stations were and how many has passed
      for (unsigned int st = 0; st < stationTraceCounter.size(); ++st){
        if (fCurStationId == fStationIds[st]){
          stationTraceCounter[st][1] += 1;
          if (filterPass == true){
            stationTraceCounter[st][0] += 1;
          }
        }
      } //end stats
      
    } // end station loop
    
    ostringstream info;
    info << "Output of the cached dataset vector after reading the event: " << endl;
    info << "*********************" << endl;
    for (unsigned int chan = 0; chan < AllChannelsAndStationsDatasets.size(); ++chan) {
      info << " Channel number:" << fchannels[chan] << " " << endl;
      for (unsigned int st = 0; st < AllChannelsAndStationsDatasets[chan].size(); ++st) {
        info << " Station number:" << fStationIds[st] << " " << endl;
        for (unsigned int row = 0; row < AllChannelsAndStationsDatasets[chan][st].size(); ++row) {
          for (unsigned int col = 0; col < AllChannelsAndStationsDatasets[chan][st][row].size(); ++col) {
            info << AllChannelsAndStationsDatasets[chan][st][row][col] << "  ";
          } // end col loop
          info << endl;
        } // end row loop
      } //end station
    } //end channel loop
    info << "*********************" << endl;
    INFODebug(info);
    info.str("");

    return eSuccess;
      
  } // end event

  VModule::ResultFlag
  RdGalacticDatasetMaker::Finish()
  {

    ostringstream info;
    
    info << "Output of the whole final dataset vector: \n";
    info << "*********************" << endl;
      for (unsigned int chan = 0; chan < AllChannelsAndStationsDatasets.size(); ++chan) {
        info << "\n";
        info << " Channel number:" << fchannels[chan] << " " << endl;
        for (unsigned int st = 0; st < AllChannelsAndStationsDatasets[chan].size(); ++st) {
          info << " Station number:" << fStationIds[st] << " " << endl;
          for (unsigned int row = 0; row < AllChannelsAndStationsDatasets[chan][st].size(); ++row){
            for (unsigned int col = 0; col < AllChannelsAndStationsDatasets[chan][st][row].size(); ++col) {
              info << AllChannelsAndStationsDatasets[chan][st][row][col] << "  ";
            } // end col loop
            info << endl;
          } // end row loop
        } //end station
      } //end channel loop
    info << "*********************" << endl;
    INFOIntermediate(info);
    info.str("");
    
    // initialize average dataset file header, this is now 3d vector (channel X rows(Time bins) X cols (freq bins)
    vector<vector<vector<float>>> 
    AllChannelsAndStationsDatasets_average(fchannels.size(), vector<vector<float>>(timeBins.size()-1, vector<float>(freqBins.size()+1)));
    
    // fill the Time columns
    for (unsigned int chan = 0; chan < AllChannelsAndStationsDatasets_average.size(); ++chan) {
      for (unsigned int row = 0; row < AllChannelsAndStationsDatasets_average[chan].size(); ++row) {
        AllChannelsAndStationsDatasets_average[chan][row][0] = timeBins[row]; // fill the Time bin times in the zeroth col
      }
    }
    
    // make average from station vectors
    vector<float> tempVector;
    int st;
    for (unsigned int chan = 0; chan < AllChannelsAndStationsDatasets.size(); ++chan) {
      st = 0;
      tempVector.clear();
      for (unsigned int row=0; row < AllChannelsAndStationsDatasets[chan][st].size(); ++row) {
        for (unsigned int col=1; col < AllChannelsAndStationsDatasets[chan][st][row].size(); ++col) {
          for (unsigned int st=0; st < AllChannelsAndStationsDatasets[chan].size(); ++st) {
            if ((AllChannelsAndStationsDatasets[chan][st][row][col] != 0) & (col!=1)) {
              tempVector.push_back(AllChannelsAndStationsDatasets[chan][st][row][col]);
              if (tempVector.size() != 0){
                float tempAverage = accumulate(tempVector.begin(), tempVector.end(), 0.0) / tempVector.size();
                AllChannelsAndStationsDatasets_average[chan][row][col] = tempAverage;
                tempVector.clear();
              } else {
                AllChannelsAndStationsDatasets_average[chan][row][col] = 0;
              }
            } //end if for freq cols
            if (col == 1){
              
              info << "Weight of the rows in each station and channel: \n";
              info << "weight test: station:" << fStationIds[st] << "channel: "<< fchannels[chan] <<   " row:"<<row<< "  " 
              << AllChannelsAndStationsDatasets[chan][st][row][col] << endl; 
              INFODebug(info);
              info.str("");
              
              // data
              AllChannelsAndStationsDatasets_average[chan][row][1] +=  AllChannelsAndStationsDatasets[chan][st][row][col];
              
            } // end if for weight
          } // end station loop
        } // end cols loop
      } //end row loop
    } // end channel loop
    
    // save  average from station and channels to file
    INFO("Saving average of station datasets for each channel.");
    
    string timeHeaderString;
    if (fdayStack == true){
      timeHeaderString = "UTChour"; // "normal" hour time
    } else {
      timeHeaderString = "LSTtime"; // LST hour time
    }
    for (unsigned int chan = 0; chan < AllChannelsAndStationsDatasets_average.size(); ++chan) {
      ofstream AllChannelsAndStationsDatasets_average_FILE;
      AllChannelsAndStationsDatasets_average_FILE.open("./datasets/dataSetFile_Ch"+to_string(fchannels[chan])+"_average.txt", ios::out | ios::app);
      AllChannelsAndStationsDatasets_average_FILE << std::setprecision(6) << fixed;
      AllChannelsAndStationsDatasets_average_FILE  << timeHeaderString << setw(16) << "weight";
      for (unsigned int i=0; i<freqBins.size()-1; ++i){
        AllChannelsAndStationsDatasets_average_FILE << setw(16) << freqBins[i];
      }
      AllChannelsAndStationsDatasets_average_FILE << endl;
    // end of saving to file
      cout << "\tStoring dataset. " << endl;
      AllChannelsAndStationsDatasets_average_FILE << std::setprecision(9) << scientific;
      for (unsigned int row = 0; row < AllChannelsAndStationsDatasets_average[chan].size(); ++row) {
        for (unsigned int col = 0; col < AllChannelsAndStationsDatasets_average[chan][row].size(); ++col) {
          AllChannelsAndStationsDatasets_average_FILE << AllChannelsAndStationsDatasets_average[chan][row][col] <<  "\t" ;
        }
        AllChannelsAndStationsDatasets_average_FILE << endl;
      }
      AllChannelsAndStationsDatasets_average_FILE.close();
    }


    info << "Output of the final averaged datasets for each channel: \n";
    info << "*********************" << endl;
    info <<"**********************" << endl;
    // test of 3d vector
    for (unsigned int chan = 0; chan < AllChannelsAndStationsDatasets_average.size(); ++chan) {
      info <<"***channel:  " << fchannels[chan] << "  ***" << endl;
      for (unsigned int row=0; row < AllChannelsAndStationsDatasets_average[chan].size(); ++row) {
        for (unsigned int col=0; col < AllChannelsAndStationsDatasets_average[chan][row].size(); ++col) {
          info << AllChannelsAndStationsDatasets_average[chan][row][col] << "  ";
        }
        info << endl;
      }
    }
    info <<  "\n" <<"**********************" << "\n";
    info <<  "\n" <<"**********************" << "\n";
    INFOIntermediate(info);
    info.str("");
    
    // save all channel and station datasets separately 
    INFO("Saving station and channel datasets (separately).");
    // initialize dataset file header
    for (unsigned chan=0; chan < AllChannelsAndStationsDatasets.size(); ++chan) {
      for (unsigned st=0; st< AllChannelsAndStationsDatasets[chan].size(); ++st) {
        ofstream dataSetFileOutput;
        dataSetFileOutput.open("./datasets/dataSetFile_Station"+to_string(fStationIds[st])+"_Ch"+to_string(fchannels[chan])+".txt", ios::out | ios::app);
        dataSetFileOutput << std::setprecision(6) << fixed;
        dataSetFileOutput  << timeHeaderString << setw(16) << "weight";
        for (unsigned int u=0; u<freqBins.size()-1; ++u) {
            dataSetFileOutput  << setw(16) << freqBins[u];
        }
        dataSetFileOutput << endl;
        dataSetFileOutput << std::setprecision(9) << scientific;
        for (unsigned int row = 0; row < AllChannelsAndStationsDatasets[chan][st].size(); ++row) {
          for (unsigned int col = 0; col < AllChannelsAndStationsDatasets[chan][st][row].size(); ++col) {
            dataSetFileOutput << AllChannelsAndStationsDatasets[chan][st][row][col] << "  ";
          } // end col loop
          dataSetFileOutput << endl;
        } // end row loop
      } //end station
    } //end channel loop
    //
      
    // just some statistics, to see at the end how many traces in each stations were and how many has passed
    INFO("How many traces in each station passed the filtering conditions from the total number of traces that were in the station:");
      for (unsigned int st = 0; st < stationTraceCounter.size(); ++st) {
        cout << "Station: " << fStationIds[st] << "  ::: " << stationTraceCounter[st][0] << "/" << stationTraceCounter[st][1] << endl;
      }

    return eSuccess;
    
  }

  float 
  RdGalacticDatasetMaker::onlineMean(float oldMean, float nextValue, int n)
  {
    return oldMean + (nextValue-oldMean)/n;
  }

  int 
  RdGalacticDatasetMaker::getChannelIndex(vector<int> vector, int element)
  { 
    auto it = find(vector.begin(), vector.end(), element);
    int index;
    // If element was found get the index of element in vector  
    if (it != vector.end()) { 
      index = distance(vector.begin(), it); 
    } else {
      ERROR("Channel number not found! \n"); //If the element is not 
      return eFailure;
    }
    return index;
  } 

  // copied from RdChannelGalacticConstantsGenerator
  // function to calculate LST
  double
  RdGalacticDatasetMaker::GPSSecToLSTHourAuger(const int gps)
  {
    const double J2000 = 2451545.0;
    const float Dtohr = 0.066666666666666;
    const time_t utc = gps + 315964800 - 18; // GPS UTC offset, leap seconds are added to utc, so they have to be substracted from GPS time

    const tm* const t = gmtime(&utc);
    int year = t->tm_year + 1900;
    int month = t->tm_mon + 1;
    const int day = t->tm_mday;
    const int hour = t->tm_hour;
    const int min = t->tm_min;
    const int sec = t->tm_sec;

//#warning can you use the utl::ModifiedJulianDate here?
    //Date to Julian Date
    if (month <= 2) {
      --year;
      month += 12;
    }
    const int a = int(floor(year / 100.));
    const int b = 2 - a + a/4;
    double jd = floor(365.25*(year + 4716)) + floor(30.6001*(month + 1)) + day + b - 1524.5;

    const double fd = double(60*60*hour + 60*min + sec) / 86400;
    jd = jd + fd;

    static double sidereal_coeff[4] = { 280.46061837, 360.98564736629, 0.000387933, 38710000 };
    const double dt = jd - J2000;
    const double dt0 = dt / 36525.;
    const double theta =
      sidereal_coeff[0] + dt * sidereal_coeff[1] +
      dt0 * dt0 * (sidereal_coeff[2] - dt0 / sidereal_coeff[3]); // is in degrees
    // original J. Meeus: longitude is positive westwards  *lst = mod(theta-longitude,360.) * kSTC::Dtohr;
    // with the new convention for longitudes (positive eastwards) :
    // the PAO longitude -69.25 is taken from GAP-2001-038
    const double lst = fmod(theta + (360 - 69.25), 360.) * Dtohr;
    return lst;
  }

  float 
  RdGalacticDatasetMaker::HHMMSS2hour(float ffHHMMSS)
  {
    float fhour = floor(ffHHMMSS/10000);
    float fminute = floor((ffHHMMSS - fhour*10000)/100);
    float fsecond = ffHHMMSS - fhour*10000 - fminute*100;
    return fhour + fminute/60 + fsecond/3600;
  }

  unsigned int
  RdGalacticDatasetMaker::TimeStamp2HHMMSS(const utl::TimeStamp& timest)
  {
    const UTCDateTime t = timest;
    return t.GetHour()*10000 + t.GetMinute()*100 + t.GetSecond();
  }
  
  double
  RdGalacticDatasetMaker::GetImpedanceAt(Branch topBranch, int channelID, float frequencyMHz)
  {
    TabulatedFunction returnData;
    Branch mainBranch = topBranch.GetFirstChild();
   // Branch mainBranch = topBranch;
    Branch impBranch;
    vector<float> frequency;
    vector<float> impedance;
    for (; mainBranch; mainBranch = mainBranch.GetNextSibling()){
      if (stoi(mainBranch.GetAttributes().begin()->second) == channelID){
        impBranch = mainBranch.GetChild("impedance");
        impBranch.GetChild("x").GetData(frequency);
        impBranch.GetChild("y").GetData(impedance);
        for (unsigned int i = 0; i < frequency.size(); ++i)
          returnData.PushBack(frequency[i], impedance[i]);
        return returnData.Y(frequencyMHz/1000);
      }
    }
     ERROR("Impedance for the channel not found!"); //If the element is not 
     return eFailure;
   }
}