RdFiller.cc

Go to the documentation of this file.

// Offline include
#include <evt/Event.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerRRecData.h>

#include <revt/REvent.h>
#include <revt/Station.h>
#include <revt/StationRecData.h>
#include <revt/StationSimData.h>
#include <revt/ChannelRecData.h>
#include <revt/Channel.h>
#include <revt/Header.h>

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

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

#include <utl/config.h>
#include <utl/CoordinateSystemPtr.h>
#include <utl/CoordinateSystem.h>
#include <utl/TimeStamp.h>
#include <utl/AugerUnits.h>
#include <utl/Point.h>
#include <utl/RadioGeometryUtilities.h>
#include <utl/UTMPoint.h>
#include <utl/AugerUnits.h>

// ADST includes
#include <RdTrace.h>
#include <Detector.h>
#include <DetectorGeometry.h>
#include <FileInfo.h>
#include "ConversionUtil.h"

// ROOT include
#include <TVector3.h>
#include <TF1.h>

#include "RdFiller.h"

using namespace utl;
using namespace otoa;
using namespace revt;
using std::vector;
using std::string;
using std::ostringstream;
using std::cerr;
using std::endl;
using std::complex;
using std::numeric_limits;
using std::out_of_range;


bool
RdFiller::FillRadioStations(const evt::Event& theev, bool saveTrace, bool storeExcludedStations)
{
  const revt::REvent& therev = theev.GetREvent();

  // increase RD RecLevel in case it is 0 and we have signal stations
  if ((therev.GetNumberOfSignalStations() > 0) &&
      (fRdEvent.GetRdRecLevel() == eNoRdEvent))
    fRdEvent.SetRdRecLevel(eHasRdSignalStations);

  for (const auto& station : therev.AllStationsRange()) {

    if (!storeExcludedStations && station.IsExcluded())
      continue;  // skip stations (do not store it)

    // Now we call (or create) the station
    // And we fill it with values from the member variables and parameter storage
    try {
      if (!fRdEvent.HasRdStation(station.GetId()))
        fRdEvent.AddRdStationWithId(station.GetId());

      RdRecStation& rstation = fRdEvent.GetRdStationById(station.GetId());

      if (station.HasRecData()) {
        const revt::StationRecData& recStation = station.GetRecData();

        // fill data from parameter storage
        const std::vector<StationRRecDataQuantities>& enumVector = recStation.GetEnumVector();
        for (const auto& index : enumVector) {
          if (recStation.HasParameter(index)) {
            rstation.SetParameter((index), recStation.GetParameter(index));
          } else {
            std::cout << "WARNING: RdFiller: Parameter " << index << " does  not exist" << std::endl;
          }
        }

        const auto& covarianceEnumVector = recStation.GetCovarianceEnumVector();
        int ncount = 0;
        for (const auto& index2 : covarianceEnumVector) {
          if (recStation.HasParameterCovariance(index2.first, index2.second)) {
            rstation.SetParameterCovariance(
              index2.first, index2.second,
              recStation.GetParameterCovariance(index2.first, index2.second)
            );
          } else
            ++ncount;
        }

        // for debugging purpose, will be removed later
        if (ncount > 0)
          WARNING(std::to_string(ncount) + string(" ParameterCovariance pairs does not exist"));

        rstation.SetHasPulse(station.GetRecData().GetPulseFound());
        rstation.SetSaturated(station.IsSaturated());
        rstation.SetExcluded(station.GetExcludedReason());
        rstation.SetRejectionStatus(station.GetRejectedReason());
      }

      // now store the trace if foreseen
      if (saveTrace) {
        const revt::StationTimeSeries& trace = station.GetStationTimeSeries();
        const revt::StationFrequencySpectrum& spectr = station.GetStationFrequencySpectrum();
        vector<Float_t> vectr[3]; // three-dimensional electric field
        vector<Float_t> vecsp[3]; // three-dimensional electric field

        RdTrace rtra;
        for (unsigned int i = 0; i < trace.GetSize(); ++i) {
          for (unsigned int j = 0; j < trace[i].GetSize(); ++j) {
            try {
              Float_t val = trace[i][j];
              val /= micro * volt / meter;
              vectr[j].push_back(val);
            } catch (std::out_of_range& e) {
              cerr << " Caught std::out_of_range\n"
                  << e.what() << endl;
            }
          }
        }

        // We do the same for the Spectrum
        for (unsigned int i = 0; i < spectr.GetSize(); ++i) {
          for (unsigned int j = 0; j < spectr[i].GetSize(); ++j) {
            try {
              Float_t val = abs(spectr[i][j]);
              val /= micro * volt / meter / megahertz;
              vecsp[j].push_back(val);
            } catch (std::out_of_range& e) {
              cerr << "Caught std::out_of_range\n"
                  << e.what() << endl;
            }
          }
        }

        for (int i = 0; i < 3; ++i) {
          rtra.SetTimeTrace(vectr[i]);
          rtra.SetSamplingRate(trace.GetBinning() / nanosecond);
          rtra.SetFreqRange(station.GetFrequencyOfBin(0) / megahertz,
                            station.GetFrequencyOfBin(spectr.GetSize() - 1) / megahertz);
          rtra.SetAbsoluteFreqSpectrum(vecsp[i]);
          rstation.SetRdTrace(rtra, i);
          vectr[i].clear();
          vecsp[i].clear();
        }

        // save trace in shower plane
        const evt::ShowerRRecData& recShower = theev.GetRecShower().GetRRecShower();
        if (recShower.HasAxis()) {

          Point core = recShower.GetCoordinateOrigin();
          if (recShower.HasCorePosition()) {
            core = recShower.GetCorePosition();
          }

          const CoordinateSystemPtr coreCS = fwk::LocalCoordinateSystem::Create(core);
          const utl::Vector& showerAxis = recShower.GetAxis();
          utl::RadioGeometryUtilities csTrans =
            utl::RadioGeometryUtilities(showerAxis, coreCS, recShower.GetMagneticFieldVector());
          const utl::TraceV3D traceShowerPlane = csTrans.GetTraceInShowerPlaneVxB(trace);
          vector<Float_t> vectraceShowerPlane[3];
          for (unsigned int i = 0; i < traceShowerPlane.GetSize(); ++i) {
            for (unsigned int j = 0; j < traceShowerPlane[i].GetSize(); ++j) {
              Float_t val = traceShowerPlane[i][j];
              val /= (micro * volt / meter);
              vectraceShowerPlane[j].push_back(val);
            }
          }
          for (int i = 0; i < 3; ++i) {
            rtra.SetTimeTrace(vectraceShowerPlane[i]);
            rtra.SetSamplingRate(traceShowerPlane.GetBinning() / nanosecond);
            rstation.SetRdTraceShowerPlane(rtra, i);
            vectraceShowerPlane[i].clear();
          }
        }

      }
    } catch (std::out_of_range& e) {
      WARNING(string("Caught std::out_of_range ") + e.what() + string(" \n"));
    }

  } // sIt
  return true;
}


bool
RdFiller::FillSimRadioStations(const evt::Event& theev)
{
  const revt::REvent& therev = theev.GetREvent();

  for (const auto& station : therev.AllStationsRange()) {
    // Now we call (or create) the station
    // And we fill it with values from the member variables and parameter storage
    try {
      if (!fRdEvent.HasRdSimStation(station.GetId()))
        fRdEvent.AddRdSimStationWithId(station.GetId());

      RdSimStation& rstation = fRdEvent.GetRdSimStationById(station.GetId());

      if (station.HasSimData()) {
        const revt::StationSimData& simStation = station.GetSimData();

        // fill data from parameter storage
        const std::vector<StationRRecDataQuantities>& enumVector = simStation.GetEnumVector();
        for (const auto& it : enumVector) {
          if (simStation.HasParameter(it)) {
            rstation.SetParameter((it), simStation.GetParameter(it));

          } else {
            std::cout << "WARNING: RdFiller: Parameter " << it << " does  not exist" << std::endl;
          }
        }

        const auto& covarianceEnumVector = simStation.GetCovarianceEnumVector();
        int ncount = 0;
        for (const auto& it2 : covarianceEnumVector) {
          if (simStation.HasParameterCovariance(it2.first, it2.second)) {
            rstation.SetParameterCovariance(
              it2.first, it2.second,
                          simStation.GetParameterCovariance(it2.first, it2.second)
            );
          } else
            ++ncount;
        }

        // for debugging purpose, will be removed later
        if (ncount > 0)
          WARNING(std::to_string(ncount) + string(" ParameterCovariance pairs does not exist"));
      }

      // now store the trace if foreseen

    } catch (std::out_of_range& e) {
      WARNING(string("Caught std::out_of_range ") + e.what() + string(" \n"));
    }

  } // sIt
  return true;
}


bool
RdFiller::FillRadioShower(const evt::Event& theev)
{
  //Shower parameters
  const CoordinateSystemPtr referenceCS = det::Detector::GetInstance().GetReferenceCoordinateSystem();
  const evt::ShowerRRecData& recShower = theev.GetRecShower().GetRRecShower();
  RdRecShower& theShower = fRdEvent.GetRdRecShower();

  try {
    const Point& coordinateOrigin = recShower.GetCoordinateOrigin();
    const CoordinateSystemPtr coordinateOriginCS = fwk::LocalCoordinateSystem::Create(coordinateOrigin);
    theShower.SetAxisSiteCS(ToTVector3(recShower.GetAxis(), referenceCS));
    theShower.SetAxisCoreCS(ToTVector3(recShower.GetAxis(), coordinateOriginCS));
    fRdEvent.SetRdRecLevel(eHasRdArrivalDirection); // there is an arrival direction available from the RD reconstruction

  } catch (utl::NonExistentComponentException& e) {
    WARNING(string("Caught utl::NonExistentComponentException") + e.what());
    // If the axis/core are not defined then they are filled with nan
    theShower.SetAxisSiteCS(TVector3(
      numeric_limits<double>::quiet_NaN(),
      numeric_limits<double>::quiet_NaN(),
      numeric_limits<double>::quiet_NaN()
    ));
    theShower.SetAxisCoreCS(TVector3(
      numeric_limits<double>::quiet_NaN(),
      numeric_limits<double>::quiet_NaN(),
      numeric_limits<double>::quiet_NaN())
    );
  }

  try {
    const utl::Point& core = recShower.GetCorePosition();
    double northing = 0;
    double easting = 0;
    double altitude = 0;

    try {
      utl::UTMPoint UTM(core, ReferenceEllipsoid::GetWGS84());
      northing = UTM.GetNorthing();
      easting = UTM.GetEasting();
      altitude = UTM.GetHeight();
      theShower.SetCoreUTMCS(TVector3(easting, northing, altitude));
    } catch (UTMPoint::UTMZoneException& zone_e) {
      ERROR ("UTMZoneException: radio rec shower invalid");
    } catch (UTMPoint::UTMException& e) {
      ERROR ("UTMException: radio rec shower invalid");
    }

    theShower.SetCoreSiteCS(ToTVector3(core, referenceCS, meter));

  } catch (utl::NonExistentComponentException& e) {
    WARNING(string("Caught utl::NonExistentComponentException") + e.what());
    // If the core are not defined then they are filled with nan
    theShower.SetCoreSiteCS(TVector3(
      numeric_limits<double>::quiet_NaN(),
      numeric_limits<double>::quiet_NaN(),
      numeric_limits<double>::quiet_NaN())
    );
  }

  const std::vector<ShowerRRecDataQuantities>& enumVector = recShower.GetEnumVector();
  for (auto it = enumVector.begin(); it != enumVector.end(); ++it) {
    if (recShower.HasParameter(*it)) {
      theShower.SetParameter((*it), recShower.GetParameter(*it));
    } else {
      std::cout << "WARNING: RdFiller: Parameter " <<  *it << " does not exist" << std::endl;
    }
  }

  const auto& covarianceEnumVector = recShower.GetCovarianceEnumVector();
  int ncount = 0;
  for (auto it2 = covarianceEnumVector.begin(); it2 != covarianceEnumVector.end(); ++it2) {
    if (recShower.HasParameterCovariance(it2->first, it2->second)) {
      theShower.SetParameterCovariance(it2->first, it2->second,recShower.GetParameterCovariance(it2->first,it2->second));
    } else {
      ++ncount;
    }
  }

  // for debugging purpose, will be removed later
  if (ncount > 0)
    WARNING(std::to_string(ncount) + string(" ParameterCovariance pairs does not exist"));

  // Set Errors of ShowerAxis in Spherical coordinates (has to get into the ADST because it is not trivial to compute the errors in spherical coordinates from a cartesian representation)
  if (recShower.HasParameter(revt::eShowerAxisX) &&
      recShower.HasParameterCovariance(eShowerAxisX, eShowerAxisX)) {  // assumes other parameter and covariances are already set
    theShower.SetRadiusError(recShower.GetRadiusError());
    theShower.SetAzimuthError(recShower.GetAzimuthError());
    theShower.SetZenithError(recShower.GetZenithError());
  }

  if (recShower.HasParameter(revt::ePreFitShowerAxisX) &&
      recShower.HasParameterCovariance(ePreFitShowerAxisX, ePreFitShowerAxisX)) {  // assumes other parameter and covariances are already set
    theShower.SetZenithPreFitError(recShower.GetZenithPreFitError());
    theShower.SetAzimuthPreFitError(recShower.GetAzimuthPreFitError());
  }

  // check whether electromagnetic energy has been reconstructed
  if (recShower.HasParameter(eReconstructedElectromagneticEnergy))
    fRdEvent.SetRdRecLevel(eHasRdElectromagneticEnergy);

  // check whether Xmax has been reconstructed
  if (recShower.HasParameter(eReconstructedXmax))
    fRdEvent.SetRdRecLevel(eHasRdDepthOfMaximum);

  return true;
}


bool
RdFiller::FillRadioChannels(const revt::REvent& therev, bool saveTrace, bool storeExcludedStations)
{
  const rdet::RDetector& rDetector = det::Detector::GetInstance().GetRDetector();

  for (const auto& station : therev.AllStationsRange()) {

    if (!storeExcludedStations && station.IsExcluded())
      continue;  // skip stations (do not store it)


    // Get the Station or Create a new one
    if (!fRdEvent.HasRdStation(station.GetId()))
      fRdEvent.AddRdStationWithId(station.GetId());

    try {
      RdRecStation& rstation = fRdEvent.GetRdStationById(station.GetId());

      // Here we save Radio data at the channel level
      for (const auto& channel : station.ChannelsRange()) {

        RdRecChannel rchannel; // create a new empty channel

        rchannel.SetId(channel.GetId());
        rchannel.SetADCSignalThreshold(channel.GetSignalThreshold());
        rchannel.SetADCNoiseThreshold(channel.GetNoiseThreshold());
        rchannel.SetSaturated(channel.IsSaturated());

        const rdet::Channel& cDetChannel =
          rDetector.GetStation(station.GetId()).GetChannel(channel.GetId());

        if (cDetChannel.GetChannelType() == "ScintillatorTop")
          rchannel.SetScintillatorTop(true);
        if (cDetChannel.GetChannelType() == "ScintillatorBottom")
          rchannel.SetScintillatorBottom(true);

        // now fill quantities form the ChannelRecData, in particular the parameter storage
        if (channel.HasRecData()) {
          const revt::ChannelRecData& recChannel = channel.GetRecData();

          // fill data from parameter storage
          const std::vector<ChannelRRecDataQuantities>& enumVector = recChannel.GetEnumVector();
          for (const auto& index : enumVector) {
            if (recChannel.HasParameter(index)) {
              rchannel.SetParameter((index),recChannel.GetParameter(index));
            } else {
              std::cout << "WARNING: RdFiller: Channel parameter " <<  index << " does  not exist" << std::endl;
            }
          }

          const auto& covarianceEnumVector = recChannel.GetCovarianceEnumVector();
          int ncount = 0;
          for (const auto& index2 : covarianceEnumVector) {
            if (recChannel.HasParameterCovariance(index2.first, index2.second)) {
              rchannel.SetParameterCovariance(index2.first, index2.second, recChannel.GetParameterCovariance(index2.first, index2.second));
            } else {
              ++ncount;
            }
          }

          // for debugging purpose, will be removed later
          if (ncount > 0)
            WARNING(std::to_string(ncount) + string(" ParameterCovariance pairs does not exist"));

        }

        // save channel traces if foreseen
        if (saveTrace) {
          std::vector<Float_t> timetrace;
          for (revt::ChannelTimeSeries::ConstIterator ctimt = channel.GetChannelTimeSeries().Begin();
               ctimt != channel.GetChannelTimeSeries().End(); ++ctimt) {
            timetrace.push_back(*ctimt / (micro * volt));
          }

          std::vector<Float_t> freqtrace;
          for (revt::ChannelFrequencySpectrum::ConstIterator cfreqit = channel.GetChannelFrequencySpectrum().Begin();
               cfreqit != channel.GetChannelFrequencySpectrum().End(); ++cfreqit) {
            freqtrace.push_back(abs(*cfreqit) / (micro * volt / megahertz));
          }

          // voltage
          rchannel.SetRdTimeTrace(timetrace);
          rchannel.SetRdAbsSpectrum(freqtrace);
          rchannel.GetRdTrace().SetFreqRange(
            channel.GetFrequencyOfBin(0) / megahertz,
            channel.GetFrequencyOfBin(channel.GetChannelFrequencySpectrum().GetSize()-1) / megahertz
          );
          rchannel.GetRdTrace().SetSamplingRate(channel.GetChannelTimeSeries().GetBinning() / nanosecond);

          // adc
          std::vector<Int_t> adctrace;
          const revt::ChannelADCTimeSeries& adctrace_offline = channel.GetChannelADCTimeSeries();
          for (unsigned int idx = 0; idx < adctrace_offline.GetSize(); ++idx)
            adctrace.push_back(adctrace_offline[idx]);

          RdADCTrace& rchannel_adc_trace = rchannel.GetRdADCTrace();
          rchannel_adc_trace.SetTimeTrace(adctrace);
          rchannel_adc_trace.SetSamplingRate(channel.GetChannelADCTimeSeries().GetBinning() / nanosecond);
          rchannel_adc_trace.SetMinVoltage(cDetChannel.GetMinVoltage());
          rchannel_adc_trace.SetMaxVoltage(cDetChannel.GetMaxVoltage());
          rchannel_adc_trace.SetBitDepth(cDetChannel.GetBitDepth());

        }

        rstation.GetRdChannelVector().push_back(rchannel); // store the channel in the RdRecStation object
      } // chanIt
    } catch (out_of_range& e) {
      WARNING(string("Caught std::out_of_range") + e.what());
    }
  }  // sIt
  return true;
}


bool
RdFiller::FillRadioEventInfo(const revt::REvent& therev)
{
  fRdEvent.SetRdEventId(therev.GetHeader().GetId());
  fRdEvent.SetRdRunNumber(therev.GetHeader().GetRunNumber());
  fRdEvent.SetRdEventTime(therev.GetHeader().GetTime().GetGPSSecond(), therev.GetHeader().GetTime().GetGPSNanoSecond());
  fRdEvent.SetRejectionStatus(therev.GetRejectionStatus());
  return true;
}


bool
RdFiller::AddRadioDetector(const det::Detector& offlineDet, Detector& adstDet)
{
  const rdet::RDetector& rdt = offlineDet.GetRDetector();
  const CoordinateSystemPtr referenceCS = offlineDet.GetReferenceCoordinateSystem();
  for (const auto& station : rdt.StationsRange()) {
    //const string& StationName = rdIt->GetName();
    const utl::Point& position = station.GetPosition();
    double x = position.GetX(referenceCS);
    double y = position.GetY(referenceCS);
    double z = position.GetZ(referenceCS);
    adstDet.SetRdStation(station.GetId(), x/m, y/m, z/m);
  }
  return true;
}


bool
RdFiller::AddRadioInfo(FileInfo* const info, const int radioSaveLevel)
{
  if (radioSaveLevel >= 1)
    info->SetHasRdStation();
  if (radioSaveLevel >= 2)
    info->SetHasRdChannel();
  return true;
}


bool
RdFiller::AddRadioDetectorGeometry(const det::Detector& offlineDet, DetectorGeometry* const adstGeo)
{
  const rdet::RDetector& rdt = offlineDet.GetRDetector();
  const CoordinateSystemPtr referenceCS = offlineDet.GetReferenceCoordinateSystem();
  for (const auto& station : rdt.StationsRange()) {
    const string& stationName = station.GetName();
    adstGeo->SetRdStationName(station.GetId(), stationName);
    double x = 0;
    double y = 0;
    double z = 0;
    try {
      const utl::Point& position = station.GetPosition();
      x = position.GetX(referenceCS);
      y = position.GetY(referenceCS);
      z = position.GetZ(referenceCS);
    } catch (NonExistentComponentException& e) {
      const string warn = "RdStation position information cannot be read from RDetector. "
                          "Station position will be set to 0,0,0";
      WARNING(warn);
    }
    adstGeo->SetRdStation(station.GetId(), x/m, y/m, z/m);
  }
  return true;
}