Detector2ADST.cc

Go to the documentation of this file.

#include "Detector2ADST.h"
#include "Config.h"
#include "ConversionUtil.h"

#include <det/Detector.h>  // DONT USE NAMESPACE det !!!!!!
#include <sdet/SDetector.h>
#include <fdet/FDetector.h>
#include <fdet/Eye.h>
#include <fdet/Telescope.h>
#include <fdet/Pixel.h>
#include <atm/OverallQualityDB.h>
#include <atm/ProfileResult.h>
#include <atm/MolecularDB.h>
#include <atm/GOESDB.h>
#include <atm/MolecularZone.h>
#include <atm/Atmosphere.h>
#include <atm/LidarDB.h>
#include <atm/LidarZone.h>

#include <fwk/CoordinateSystemRegistry.h>

#include <evt/Event.h>
#include <sevt/SEvent.h>
#include <fevt/Eye.h>
#include <fevt/Telescope.h>
#include <fevt/Pixel.h>
#include <fevt/FEvent.h>

#include <utl/CoordinateSystemPtr.h>
#include <utl/ErrorLogger.h>
#include <utl/Point.h>
#include <utl/PhysicalConstants.h>

#include <adst/Detector.h>
#include <adst/RecEvent.h>
#include <adst/FdLidarData.h>
#include <adst/FdCloudCameraData.h>
#include <adst/FdAerosols.h>

#include <cmath>
#include <sstream>
#include <iostream>
#include <unistd.h>

#include <mdet/MDetector.h>
#include <mevt/MEvent.h>

using namespace fwk;
using namespace fdet;
using namespace sdet;
using namespace fevt;
using namespace sevt;
using namespace utl;
using namespace atm;
using namespace std;


namespace otoa {

  Detector2ADST::Detector2ADST(const otoa::Config& config) :
    fConfig(config)
  { }


  void
  Detector2ADST::Convert(const evt::Event& event, RecEvent& recEvent)
  {
    FillDetector(event, recEvent.GetDetector());

    if (GetConfig().StoreLidarData())
      FillLidarData(recEvent);

    if (GetConfig().StoreCloudCameraData() == 2)
      FillCloudCameraData(event, recEvent.GetDetector());

    if (GetConfig().StoreGOESData() == 2)
      FillGOESData(recEvent.GetDetector());
  }


  void
  Detector2ADST::FillGOESData(Detector& detector)
  {
    try {
      const Atmosphere& atm = det::Detector::GetInstance().GetAtmosphere();
      const GOESDB& goesDB = atm.GetGOESDB();
      for (unsigned int iPix = 0; iPix < goesDB.GetNumberOfPixels(); ++iPix) {
        if (goesDB.HasData(iPix))
          detector.SetCloudFraction(iPix, goesDB.GetCloudProbability(iPix));
        else
          return;  // give up upon first HasNot()
      }
    } catch (utl::AugerException& augerException) {
      WARNING(augerException.GetExceptionName() + ", " +
              augerException.GetMessage());
      detector.ClearCloudMap();  // do not store half-filled cloud data
    }
  }


  void
  Detector2ADST::FillDetector(const evt::Event& theEvent, Detector& theDet)
  {
    const det::Detector& detector = det::Detector::GetInstance();

    // -------------> SD <-----------------

    TBits stationBits(1600);  // TBits enlarges automatically of needed!

    if (theEvent.HasSEvent() && fConfig.WriteSD()) {

      const sevt::SEvent& sEvent = theEvent.GetSEvent();

      const sdet::SDetector& sDet = detector.GetSDetector();
      const CoordinateSystemPtr cs = detector.GetReferenceCoordinateSystem();

      // loop on all stations
      for (sevt::SEvent::ConstStationIterator sIt = sEvent.StationsBegin();
           sIt != sEvent.StationsEnd(); ++sIt) {
        const int sId = sIt->GetId();
        try {
          if (sDet.GetStation(sId).IsDense()) {
            const TVector3 p = ToTVector3(sDet.GetStation(sId).GetPosition(), cs, meter);
            theDet.SetDenseStation(sId, p.X(), p.Y(), p.Z());
          }
        } catch (AugerException& e) {
          ERROR(e.GetMessage());
          continue;
        }

        stationBits.SetBitNumber(sId, true);
      }

    }

    //-------------> Md (AMIGA) <-----------------
    if (theEvent.HasMEvent()) {

      const mevt::MEvent& mEvent = theEvent.GetMEvent();

      const mdet::MDetector& mDet = det::Detector::GetInstance().GetMDetector();
      const CoordinateSystemPtr cs = det::Detector::GetInstance().GetReferenceCoordinateSystem();

      // loop on all stations
      for (mevt::MEvent::CounterConstIterator mIt = mEvent.CountersBegin();
           mIt != mEvent.CountersEnd(); ++mIt) {
        const int cId = mIt->GetId();
        try {
          if (mDet.GetCounter(cId).IsDense()) {
            const TVector3 p = ToTVector3(mDet.GetCounter(cId).GetPosition(), cs, meter);
            theDet.SetDenseStation(cId, p.X(), p.Y(), p.Z());
          }
        } catch (AugerException& e) {
          ERROR(e.GetMessage());
          continue;
        }

        stationBits.SetBitNumber(cId, true);
      }

    }

    // -------------> FD <-----------------

    TBits eyeBits = TBits(10);  // TBits enlarges automatically of needed!

    const bool hasFdConfig =
      detector.GetFDetector().EyesBegin() != detector.GetFDetector().EyesEnd();

    if (fConfig.WriteFD() && theEvent.HasFEvent() && hasFdConfig) {
      const Atmosphere& theAtm = det::Detector::GetInstance().GetAtmosphere();

      theDet.SetHasFdUpTime(detector.GetFDetector().HasFdUpTime());

      const fevt::FEvent& fevent = theEvent.GetFEvent();

      // loop on all eyes
      for (FEvent::ConstEyeIterator iEye = fevent.EyesBegin(fevt::ComponentSelector::eInDAQ);
           iEye != fevent.EyesEnd(fevt::ComponentSelector::eInDAQ); ++iEye) {

        const unsigned int thisEyeId = iEye->GetId();

        try {
          theDet.SetVAODAtReferenceHeight(
            thisEyeId,
            theAtm.GetVerticalAerosolOpticalDepth(
              thisEyeId,
              GetConfig().VAODReferenceHeight()
            )
          );
        } catch (AugerException& zone_e) {
           INFO("VAOD at reference not set!!!");
          // and do nothing about it
        }

        if (GetConfig().AerosolAttenuationLengthSampling() > 0)
          FillAerosols(thisEyeId, theDet);

        eyeBits.SetBitNumber(thisEyeId, true);

        const fdet::Eye& eyeDet = detector.GetFDetector().GetEye(thisEyeId);
        double sd_offset = eyeDet.GetSDTimeOffset().GetInterval();

        theDet.SetFDSDTimeOffset (thisEyeId, sd_offset);

        // telescopes

        TBits daqBits(10);
        for (fevt::Eye::ConstTelescopeIterator iTel = iEye->TelescopesBegin(fevt::ComponentSelector::eInDAQ);
             iTel != iEye->TelescopesEnd(fevt::ComponentSelector::eInDAQ); ++iTel) {

          const unsigned int thisTelId = (unsigned int)iTel->GetId();
          daqBits.SetBitNumber(thisTelId, true);

          const fdet::Telescope& telDet = eyeDet.GetTelescope(thisTelId);

          // pixels
          for (fevt::Telescope::ConstPixelIterator iPix = iTel->PixelsBegin(ComponentSelector::eExists);
               iPix != iTel->PixelsEnd(ComponentSelector::eExists); ++iPix) {

            const unsigned int pixelId = iPix->GetId();

            if (telDet.GetPixel(pixelId).GetChannelId() != pixelId)
              theDet.AddMismatchedPixel(thisEyeId, thisTelId, pixelId);

            if (iPix->GetStatus() == ComponentSelector::eDeSelected ||
                iPix->GetStatus() == ComponentSelector::eBadCalibration)
              theDet.AddBadPixel(thisEyeId, thisTelId, pixelId);
          }
        }

        daqBits.Compact();
        theDet.SetMirrorsInDAQ(thisEyeId, daqBits);

      }
    }
    eyeBits.Compact();
    stationBits.Compact();
    theDet.SetActiveStations(stationBits);
    theDet.SetActiveEyes(eyeBits);

    const fdet::FDetector& fdDet = det::Detector::GetInstance().GetFDetector();

    for (fdet::FDetector::EyeIterator eyeiter = fdDet.EyesBegin(fdet::FDetComponentSelector::eAll);
         eyeiter != fdDet.EyesEnd(fdet::FDetComponentSelector::eAll); ++eyeiter) {
      TBits corrBits(7);
      for (fdet::Eye::TelescopeIterator teliter = eyeiter->TelescopesBegin();
           teliter != eyeiter->TelescopesEnd(); ++teliter) {
        if (teliter->HasCorrectorRing())
          corrBits.SetBitNumber(teliter->GetId(), true);

        theDet.SetPointingId(eyeiter->GetId(), teliter->GetId(), teliter->GetTelescopePointingId());
      }

      theDet.SetCorrectorRingStatus(eyeiter->GetId(),corrBits);
    }

    // -------------> Atmosphere <-----------------
    int lastEventGPSSecond = 0;
    //double lastPressure = 0;
    //double lastTemperature = 0;
    if (GetConfig().UseWeatherStations() == 1) {

      const Atmosphere& theAtm = det::Detector::GetInstance().GetAtmosphere();
      const int newTime = theEvent.GetHeader().GetTime().GetGPSSecond();
      const int kTimeInterval = 30*60; //seconds
      if (abs(lastEventGPSSecond - newTime) >= kTimeInterval) {
        lastEventGPSSecond = theEvent.GetHeader().GetTime().GetGPSSecond();
        theDet.SetTemperaturePressureGPSSeconds(lastEventGPSSecond);
        const MolecularDB& molWeather = theAtm.GetMolecularDB(MolecularIds::eWeatherStation);
        for (MolecularDB::ZoneIterator zIt = molWeather.ZonesBegin();
             zIt != molWeather.ZonesEnd(); ++zIt) {
          for (MolecularZone::LayerIterator lIt = zIt->LayersBegin();
               lIt != zIt->LayersEnd(); ++lIt) {
            if (lIt->GetTemperature()) {
              theDet.SetPressure(lIt->GetPressure()/bar);
              theDet.SetTemperature(lIt->GetTemperature());
              theDet.SetHasWeatherStationData();
              //lastPressure = lIt->GetPressure()/bar;
              //lastTemperature = lIt->GetTemperature();
              break;
            }
          }
          if (theDet.HasWeatherStationData())
            break;
        }

        if (!theDet.HasWeatherStationData()) {
          const ProfileResult& temp = theAtm.EvaluateTemperatureVsHeight();
          const ProfileResult& pressure = theAtm.EvaluatePressureVsHeight();
          theDet.SetTemperature(temp.Y(1412*meter));
          theDet.SetPressure(pressure.Y(1412*meter)/bar);
          //lastPressure = pressure.Y(1412*meter)/bar;
          //lastTemperature = temp.Y(1412*meter);
        }
      }

    }

    det::Detector& det = det::Detector::GetInstance();
    if (GetConfig().ConnectToDatabases()) {
      try {
        const OverallQualityDB& overallDB = det.GetAtmosphere().GetOverallQualityDB();
        const double qDBHorizonalUniformity = overallDB.GetHorizontalUniformity();
        const double qDBCloudCoverage = overallDB.GetCloudCoverage();
        const double qDBMinCloudBase = overallDB.GetMinCloudBase();

        const atm::ProfileResult& depthProfile =
          det::Detector::GetInstance().GetAtmosphere().EvaluateDepthVsHeight();
        double maxHeight = depthProfile.MaxX();
        double minHeight = depthProfile.MinX();

        theDet.SetHasOverallQualityDatabase();
        theDet.SetQDBHorizontalUniformity(qDBHorizonalUniformity);
        theDet.SetQDBCloudCoverage(qDBCloudCoverage);
        const double h = qDBMinCloudBase*m;
        theDet.SetQDBMinCloudBaseHeight(h);
        if (h >= maxHeight)
          theDet.SetQDBMinCloudBaseDepth(depthProfile.Y(maxHeight)/g*cm2);
        else if (h <= minHeight)
          theDet.SetQDBMinCloudBaseDepth(depthProfile.Y(minHeight)/g*cm2);
        else
          theDet.SetQDBMinCloudBaseDepth(depthProfile.Y(h)/g*cm2);
      } catch (utl::AugerException& augerException) {
        WARNING(augerException.GetExceptionName() + ", " +
                augerException.GetMessage());
        theDet.SetQDBMinCloudBaseHeight(-1);
        theDet.SetQDBMinCloudBaseDepth(-1);
        theDet.SetQDBHorizontalUniformity(-1);
        theDet.SetQDBCloudCoverage(-1);
      }
    } else {
      theDet.SetQDBMinCloudBaseHeight(-1);
      theDet.SetQDBMinCloudBaseDepth(-1);
      theDet.SetQDBHorizontalUniformity(-1);
      theDet.SetQDBCloudCoverage(-1);
    }
  }


  void
  Detector2ADST::FillLidarData(RecEvent& theRecEvent)
  {
    const Atmosphere& theAtm = det::Detector::GetInstance().GetAtmosphere();
    try {
      const LidarDB& LidarDB = theAtm.GetLidarDB();
      for (LidarDB::ZoneIterator zIt = LidarDB.ZonesBegin();
           zIt != LidarDB.ZonesEnd(); ++zIt) {

        const string lidarZoneName = zIt->GetName();
        const double lidarCloudCoverage = zIt->GetCloudCoverage();
        double lidarCloudHeight = zIt->GetLowestCloudHeight();
        double lidarCloudThickness = zIt->GetLowestCloudThickness();
        double lidarCloudVAOD = zIt->GetLowestCloudVAOD();
//#warning update LidarRange once the information is available
        const double lidarRange = -1;

        //-------  Depth Conversion ---------------------------
        double cosZenith = -1;
        double lidarCloudDepth = -1;
        double lidarCloudSlantThickness = -1;

        int hottestEye = theRecEvent.GetHottestEyeId();
        if (hottestEye && theRecEvent.GetEye(hottestEye).GetRecLevel() > eHasAxis)
           cosZenith = theRecEvent.GetEye(hottestEye).GetFdRecShower().GetCosZenith();
        else
           if (theRecEvent.GetSDEvent().HasAxis())
             cosZenith = theRecEvent.GetSDEvent().GetSdRecShower().GetCosZenith();

        if (cosZenith > 0) {

          const atm::ProfileResult& depthProfile =
            det::Detector::GetInstance().GetAtmosphere().EvaluateDepthVsHeight();
          double maxHeight = depthProfile.MaxX();
          double minHeight = depthProfile.MinX();

          if (!lidarCloudCoverage) {
            // no clouds
            lidarCloudHeight = maxHeight;
            lidarCloudThickness = 0;
            lidarCloudVAOD = 0;
          }

          double lidarCloudMax = lidarCloudHeight+lidarCloudThickness;
          double lidarCloudSlantMax = -1.;

          if (lidarCloudHeight >= maxHeight)
            lidarCloudDepth = depthProfile.Y(maxHeight)/g*cm2;
          else if (lidarCloudHeight<= minHeight)
            lidarCloudDepth = depthProfile.Y(minHeight)/g*cm2;
          else
            lidarCloudDepth = depthProfile.Y(lidarCloudHeight)/g*cm2;

          if (lidarCloudMax >= maxHeight)
            lidarCloudSlantMax = depthProfile.Y(maxHeight)/g*cm2;
          else if (lidarCloudMax <= minHeight)
            lidarCloudSlantMax = depthProfile.Y(minHeight)/g*cm2;
          else
            lidarCloudSlantMax = depthProfile.Y(lidarCloudMax)/g*cm2;

          // simple transformation from vertical to slant depth
          double CloudDepthSlant = lidarCloudDepth/cosZenith;
          lidarCloudSlantMax /= cosZenith;
          if (lidarCloudSlantMax < CloudDepthSlant)
            lidarCloudSlantThickness = CloudDepthSlant - lidarCloudSlantMax;
        }

        int lidarZoneID = 0;
        if (lidarZoneName == "Los Leones")
          lidarZoneID = 1;
        else if (lidarZoneName == "Los Morados")
          lidarZoneID = 2;
        else if (lidarZoneName == "Loma Amarilla")
          lidarZoneID = 3;
        else if (lidarZoneName == "Coihueco")
          lidarZoneID = 4;

        FdLidarData thisLidarData(lidarCloudCoverage, lidarCloudHeight,
                                  lidarCloudDepth, lidarCloudThickness,
                                  lidarCloudSlantThickness, lidarCloudVAOD,
                                  lidarRange);
        theRecEvent.GetDetector().SetLidarData(lidarZoneID, thisLidarData);

      }

    } catch(AugerException& e) {
      WARNING("Exception retrieving lidar info.");
    }
  }


  void
  Detector2ADST::FillCloudCameraData(const evt::Event& theEvent,
                                     Detector& adstDetector)
  {
    const FDetector& theFDet = det::Detector::GetInstance().GetFDetector();

    if (!theEvent.HasFEvent())
      return;

    //loop over eyes
    for (fevt::FEvent::ConstEyeIterator iEye = theEvent.GetFEvent().EyesBegin(ComponentSelector::eHasData);
         iEye != theEvent.GetFEvent().EyesEnd(ComponentSelector::eHasData); ++iEye) {

      const unsigned int thisEyeId = iEye->GetId();
      FdCloudCameraData thisCloudCameraData;

      //loop over telescopes
      for (fevt::Eye::ConstTelescopeIterator iTel = iEye->TelescopesBegin(ComponentSelector::eHasData);
           iTel != iEye->TelescopesEnd(ComponentSelector::eHasData); ++iTel) {

        const unsigned int thisTelId = iTel->GetId();
        const fdet::Telescope& detTel = theFDet.GetEye(thisEyeId).GetTelescope(thisTelId);

        //loop over det::pixels
        for (unsigned int pixelId = detTel.GetLastPixelId(); //inverse order due to resize in FdCloudCameraData
             pixelId >= detTel.GetFirstPixelId(); --pixelId) {

          const fdet::Pixel& iPix = detTel.GetPixel(pixelId);
          if (!iPix.HasCloudFraction()) //maybe thisCloudCameraData should be cleared here
            break; //next tel if one pixel is empty

          const float cloudFraction = iPix.GetCloudFraction();

          thisCloudCameraData.SetCloudFraction(thisTelId, pixelId, cloudFraction);
        }
      }
      if (!thisCloudCameraData.IsEmpty())
        adstDetector.SetCloudCameraData(thisEyeId, thisCloudCameraData);
    }
  }


  void
  Detector2ADST::FillAerosols(const unsigned int thisEyeId, Detector& theDet)
  {
    det::Detector& offlDet = det::Detector::GetInstance();
    const Atmosphere& theAtm = offlDet.GetAtmosphere();
    const fdet::Eye& detEye = offlDet.GetFDetector().GetEye(thisEyeId);

    CoordinateSystemPtr eyeCS = detEye.GetEyeCoordinateSystem();
    const Point& eyepos = detEye.GetPosition();

    ReferenceEllipsoid wgs84(ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84));

    const double eyeHeight = (wgs84.PointToLatitudeLongitudeHeight(eyepos)).get<2>();
    const double maxHeightRelEye = 20000*meter - eyeHeight;

    const double wavelength = GetConfig().LaserWavelength();
    const double stepSize = GetConfig().AerosolAttenuationLengthSampling();

    vector<double> heights;
    vector<double> attenuationLengths;
    double heightRelEye = 0;
    while (heightRelEye < maxHeightRelEye) {
      const Point samplePoint(0, 0, heightRelEye, eyeCS);

      const double attLength =
        theAtm.GetMieAttenuationLength(samplePoint, wavelength);

      heights.push_back(heightRelEye + eyeHeight); // save abs height
      attenuationLengths.push_back(attLength);

      heightRelEye += stepSize;
    }

    FdAerosols aerosolData;
    aerosolData.SetMieAttenuationLength(attenuationLengths);
    aerosolData.SetMieAttenuationHeight(heights);

    theDet.SetAerosolData(thisEyeId, aerosolData);
  }

}