Offline2ADST.cc

Go to the documentation of this file.

#include <utl/config.h>

#include "Offline2ADST.h"
#include "ConversionUtil.h"

#include <adst/RecEvent.h>
#include <adst/DetectorGeometry.h>
#include <adst/FileInfo.h>
#include <adst/GenShower.h>
#include <adst/GenParticle.h>

#include <det/Detector.h> // DON'T USE NAMESPACE det !!!!!!

#include <sdet/SDetector.h>

#include <fdet/FDetector.h>
#include <fdet/Pixel.h>
#include <fdet/Eye.h>
#include <fdet/Camera.h>
#include <fdet/Telescope.h>

#include <atm/Atmosphere.h>
#include <atm/InclinedAtmosphericProfile.h>
#include <atm/ProfileResult.h>
#include <atm/GOESDB.h>

#include <evt/Event.h>
#include <evt/ShowerSimData.h>
#include <evt/GenParticle.h>
#include <evt/RadioSimulation.h>
#include <fevt/FEvent.h>
#include <fevt/Eye.h>
#include <fevt/Telescope.h>

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

#include <utl/Branch.h>
#include <utl/ErrorLogger.h>
#include <utl/UTCDateTime.h>
#include <utl/TimeStamp.h>
#include <utl/ModifiedJulianDate.h>
#include <utl/MathConstants.h>
#include <utl/ReferenceEllipsoid.h>
#include <utl/AugerUnits.h>
#include <utl/Point.h>
#include <utl/Vector.h>
#include <utl/UTMPoint.h>
#include <utl/Transformation.h>
#include <utl/Trace.h>
#include <utl/AugerUnits.h>
#include <utl/MathConstants.h>
#include <utl/PhysicalConstants.h>
#include <utl/PhysicalFunctions.h>
#include <evt/VGaisserHillasParameter.h>
#include <utl/config.h>
#include <utl/Particle.h>
#include <utl/CoordinateSystemPtr.h>

#include <iostream>
#include <string>
#include <limits>
#include <unistd.h>

#include "SD2ADST.h"
#include "FD2ADST.h"
#include "Detector2ADST.h"
#include "ErrorPropagation.h"

#include <mdet/MDetector.h>
#include <mdet/Counter.h>
#include "MD2ADST.h"

#include <adst/RdEvent.h>
#include <adst/RdTrace.h>
#include <evt/ShowerRecData.h>
#include "RdFiller.h"

using namespace std;
using namespace otoa;
using namespace utl;


Offline2ADST::Offline2ADST(const otoa::Config& cfg) :
  fConfig(cfg)
{ }


bool
Offline2ADST::Convert(const evt::Event& event, RecEvent& recEvent)
{
  const auto& fDet = det::Detector::GetInstance().GetFDetector();
  const bool hasFdConfig = (fDet.EyesBegin() != fDet.EyesEnd());

  recEvent = RecEvent(); // Clear
  fStatus = otoa::Status(event); // reset conversion status
  fStatus.SetIsMCEvent(event.HasSimShower());

  FillGlobalInfo(event, recEvent);

  if (fConfig.WriteSD() && event.HasSEvent()) {
    INFO("Convert SD");
    SD2ADST(Config()).Convert(event, recEvent);
  }

  if (fConfig.WriteFD() && event.HasFEvent() && hasFdConfig) {
    INFO("Convert FD");
    FD2ADST(*this, fStatus).Convert(event, recEvent);
  }

  if (event.HasMEvent()) {
    INFO("Convert MD");
    MD2ADST(Config()).Convert(event, recEvent);
  }

  if (event.HasSimShower()) {
    bool saveProfiles = true;
    if (Config().DropUntriggeredMCProfiles())
      saveProfiles = recEvent.GetNEyes() > 0 || recEvent.GetSDEvent().GetRecLevel() > 0;
    FillSim(event, recEvent.GetGenShower(), saveProfiles);
  }

  if (!Config().DropUntriggeredDetector() ||
      recEvent.GetNEyes() > 0 ||
      recEvent.GetSDEvent().GetRecLevel() > 0) {
    INFO("Convert Detector");
    Detector2ADST(Config()).Convert(event, recEvent);
  }

  FillRadio(event, recEvent);

  FillDetectorGeometryIncremental();

  return true;
}


void
Offline2ADST::FillRadio(const evt::Event& event, RecEvent& recEvent)
{
  if (!event.HasREvent())
    return;

  const int storeRadioLevel = fConfig.StoreRadioLevel();
  if (storeRadioLevel < 0)
    return;

  INFO("Convert RD");

  const bool storeStationTraces = storeRadioLevel >= 1;
  const bool storeChannelTraces = storeRadioLevel >= 2;
  const bool storeExcludedStations = fConfig.StoreExcludedRdStations();

  const auto& rev = event.GetREvent();

  RdFiller rfiller;
  rfiller.AddRadioDetector(det::Detector::GetInstance(), recEvent.GetDetector());
  rfiller.FillRadioEventInfo(rev);
  if (event.HasRecShower() && event.GetRecShower().HasRRecShower())
    rfiller.FillRadioShower(event);
  rfiller.FillRadioStations(event, storeStationTraces, storeExcludedStations);
  rfiller.FillSimRadioStations(event);
  rfiller.FillRadioChannels(rev, storeChannelTraces, storeExcludedStations);

  recEvent.GetRdEvent() = rfiller.GetRdEvent();
}


void
Offline2ADST::FinishDetectorAndFileInfo()
{
  // Fill the FileInfo
  FileInfo& info = GetFileInfo();

  switch (fConfig.StoreSDTraces()) {
  case 2:
    info.SetHasSdAllTraces();
    // fall through
  case 1:
    info.SetHasVEMTraces();
    break;
  default:
    break;
  };

  // above code does not work if SdHorizontalReconstruction is used
  info.SetLDFType(eTabulated);

  if (fStatus.IsMCEvent())  // FIXME this should be done once per event
    info.SetHasMC();

  switch (fConfig.StoreFDTraces()) {
  case 2:
    info.SetHasFdSpotRecTraces();
    // fall through
  case 1:
    info.SetHasFdTraces();
    break;
  default:
    break;
  };

  const unsigned int hostLen = 80;
  char host[hostLen+1] = { '\0' };  // +1 = ensure zero termination when hostname is long/truncated
  gethostname(host, hostLen);

  info.SetOfflineVersion(OFFLINE_PACKAGE_VERSION); // from Offline utl/config.h
  info.SetHost(host);
  //info.SetHost(getenv("HOSTNAME"));
  info.SetUser(getenv("USER") ? getenv("USER") : "unknown");

  fwk::CentralConfig* const cc = fwk::CentralConfig::GetInstance();
  info.SetOfflineConfiguration(cc->GetConfig());

  ostringstream recEventVersion;
  recEventVersion << RecEvent::GetReleaseName() << ' ' << RecEvent::GetVersion()
                  << " (SVN: " << RecEvent::GetSVNVersion() << " rev "<< RecEvent::GetSVNRevision() << ')';
  info.SetRecEventVersion(recEventVersion.str());

  RdFiller::AddRadioInfo(&fFileInfo, fConfig.StoreRadioLevel());

  const CoordinateSystemPtr referenceCS = det::Detector::GetInstance().GetReferenceCoordinateSystem();

  const sdet::SDetector& sdDetector = det::Detector::GetInstance().GetSDetector();
  const DetectorGeometry::StationPosMap& adstStations = fDetectorGeometry.GetStations();

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

    if (station.IsDense())
      continue;

    const int id = station.GetId();

    // Fill only new stations
    // Note: It may be more reasonable to prefer newer information over older,
    // but an attempt to implement that was thwarted by an "overwriting" warning
    // in the DetectorGeometry.
    if (adstStations.find(id) != adstStations.end())
      continue;

    TVector3 pos;

    try {
      pos = ToTVector3(station.GetPosition(), referenceCS, meter);
    } catch (utl::NonExistentComponentException& e) {
      WARNING("Station position information cannot be read from SDetector. "
              "Probably this station is handled by the EventStationPositionsManager"
              "that is already shut down at this stage ... "
              "Station position will be set to 0,0,0.");
    }
    int gridType = 0;
    const int partners = station.GetNumberOfPartners();
    if (station.IsInGrid(sdet::SDetectorConstants::eInfill750))
      gridType |= eInfill;
    if (station.IsInGrid())
      gridType |= eRegularArray;
    if (partners == 1)
      gridType |= eDoublet;
    if (partners == 2)
      gridType |= eTriplet;

    fDetectorGeometry.SetStation(id, pos.X()/m, pos.Y()/m, pos.Z()/m, gridType);

    TString name = "unknown";
    try {
      name = station.GetName();
    } catch (utl::NonExistentComponentException& e) {
      WARNING("Station name information cannot be read from SDetector. "
              "Probably this station is handled by the EventStationPositionsManager"
              "that is already shut down at this stage ... "
              "Station name will be set to \"unknown\".");
    }
    fDetectorGeometry.SetStationName(id, name);
  }

  if (fConfig.StoreRadioLevel() >= 0)
    RdFiller::AddRadioDetectorGeometry(det::Detector::GetInstance(), &fDetectorGeometry);
}


void
Offline2ADST::FillGlobalInfo(const evt::Event& event, RecEvent& recEvent)
{
  recEvent.SetEventId(event.GetHeader().GetId());
  const auto& ts = event.GetHeader().GetTime();
  const int yymmdd = otoa::TimeStamp2YYMMDD(ts);
  const int hhmmss = otoa::TimeStamp2HHMMSS(ts);
  recEvent.SetYYMMDD(yymmdd);
  recEvent.SetHHMMSS(hhmmss);
  recEvent.SetMJD(ModifiedJulianDate(ts));
}


void
Offline2ADST::FillDetectorGeometryIncremental()
{
  const auto referenceCS = det::Detector::GetInstance().GetReferenceCoordinateSystem();
  const auto& fdDetector = det::Detector::GetInstance().GetFDetector();

  for (auto eyeiter = fdDetector.EyesBegin(fdet::FDetComponentSelector::eAll),
       end = fdDetector.EyesEnd(fdet::FDetComponentSelector::eAll); eyeiter != end; ++eyeiter) {

    const auto& eyeCS = eyeiter->GetEyeCoordinateSystem();
    const auto eyeId = eyeiter->GetId();

    if (!fDetectorGeometry.HasEye(eyeId)) {
      const auto pos = ToTVector3(eyeiter->GetPosition(), referenceCS, meter);
      const Vector vertEye(0,0,1, eyeCS);
      const double phiEye = vertEye.GetPhi(referenceCS);
      const double thetaEye = vertEye.GetTheta(referenceCS);
      const Vector backVec(1,0,0, eyeCS);
      double backwallAngle = atan2(backVec.GetY(referenceCS), backVec.GetX(referenceCS));
      if (backwallAngle > kTwoPi)
        backwallAngle -= kTwoPi;
      if (backwallAngle < 0)
        backwallAngle += kTwoPi;

      fDetectorGeometry.SetEye(eyeId,
                               pos.X()/m, pos.Y()/m, pos.Z()/m,
                               backwallAngle,
                               phiEye, thetaEye,
                               eyeiter->GetName(),
                               eyeiter->GetNameAbbr());
    }

    auto& eyeGeo = fDetectorGeometry.GetEye(eyeId);
    for (auto teliter = eyeiter->TelescopesBegin(), end = eyeiter->TelescopesEnd();
         teliter != end; ++teliter) {
      const int telId = teliter->GetId();

      if (!eyeGeo.HasTelescope(telId)) {
        // Add empty containers by default, waiting for the incremental filling
        map<TString, Double_t> adstTelPPhi;
        map<TString, Double_t> adstTelPElevation;
        map<TString, TelescopeGeometry::PixelList_t> pixelPhi;
        map<TString, TelescopeGeometry::PixelList_t> pixelOmega;
        fDetectorGeometry.SetTelescope(eyeId, telId,
                                       adstTelPPhi, adstTelPElevation,
                                       teliter->GetCamera().GetFADCTraceLength(),
                                       teliter->GetCamera().GetFADCBinSize(),
                                       pixelPhi, pixelOmega);
      }

      const unsigned int firstpixelid = teliter->GetFirstPixelId();
      const unsigned int lastpixelid = teliter->GetLastPixelId();

      const TString curPointingId = TString(teliter->GetTelescopePointingId());
      auto& telGeo = eyeGeo.GetTelescope(telId);
      auto& solidAngles = telGeo.GetSolidAngles();
      if (solidAngles.empty())
        for (unsigned int i = firstpixelid; i <= lastpixelid; ++i)
          solidAngles.push_back(teliter->GetPixel(i).GetSolidAngle());

      if (!telGeo.HasPointing(curPointingId)) {
        // Get the telescope pointing
        const Double_t elevation = kPi/2 - teliter->GetAxis().GetTheta(eyeCS);
        const Double_t azimuth = teliter->GetAxis().GetPhi(eyeCS);
        TelescopeGeometry::PixelList_t pixelPhi(lastpixelid - firstpixelid + 1);
        TelescopeGeometry::PixelList_t pixelOmega(lastpixelid - firstpixelid + 1);
        for (unsigned int i = firstpixelid; i <= lastpixelid; ++i) {
          const auto& detpixel = teliter->GetPixel(i);
          const auto& dir = detpixel.GetDirection();
          const double theta = dir.GetTheta(eyeCS) / utl::degree;
          double phi = dir.GetPhi(eyeCS) / utl::degree;
          if (phi < -60)
            phi += 360;

          pixelPhi[i-1] = phi;
          pixelOmega[i-1] = 90 - theta;
        }

        telGeo.AddPointing(curPointingId, azimuth, elevation, pixelPhi, pixelOmega);
      }

    }

  }

  // Now incrementally fill stations
  if (Verbosity() > 0)
    INFO("Updating SD stations");

  const auto& sdDetector = det::Detector::GetInstance().GetSDetector();

  const auto& adstStations = fDetectorGeometry.GetStations();

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

    if (station.IsDense())
      continue;

    const int id = station.GetId();

    // Fill only new stations
    // Note: It may be more reasonable to prefer newer information over older,
    // but an attempt to implement that was thwarted by an "overwriting" warning
    // in the DetectorGeometry.
    if (adstStations.find(id) != adstStations.end())
      continue;

    try {

      const auto pos = ToTVector3(station.GetPosition(), referenceCS, meter);
      int gridType = 0;
      const int partners = station.GetNumberOfPartners();
      if (station.IsInGrid(sdet::SDetectorConstants::eInfill750))
        gridType |= eInfill;
      if (station.IsInGrid())
        gridType |= eRegularArray;
      if (partners == 1)
        gridType |= eDoublet;
      if (partners == 2)
        gridType |= eTriplet;

      fDetectorGeometry.SetStation(id, pos.X()/m, pos.Y()/m, pos.Z()/m, gridType);

    } catch (NonExistentComponentException& e) {
      WARNING("Station position information cannot be read from SDetector.");
    }

    try {
      const auto name = station.GetName();
      fDetectorGeometry.SetStationName(id, name);
    } catch (NonExistentComponentException& e) {
      WARNING("Station name information cannot be read from SDetector.");
    }

    try {
      const double gpsSecond = station.GetCommissionTime().GetGPSSecond();
      fDetectorGeometry.SetStationCommissionTime(id, gpsSecond);
    } catch (NonExistentComponentException& e) {
      WARNING("Station commission time information cannot be read from SDetector.");
    }

  }

  // fill cloud on array pixelization if not yet done
  auto& cloudPixCenters = fDetectorGeometry.GetCloudPixelCenters();
  if (cloudPixCenters.empty()) {
    try {
      const auto& atm = det::Detector::GetInstance().GetAtmosphere();
      const auto& goesDB = atm.GetGOESDB();
      fDetectorGeometry.SetCloudPixelDimensions(goesDB.GetPixelWidthEasting(),
                                                goesDB.GetPixelWidthNorthing());
      for (unsigned int iPix = 0; iPix < goesDB.GetNumberOfPixels(); ++iPix) {
        const auto center = ToTVector3(goesDB.GetPixelCenter(iPix).GetPoint(referenceCS), referenceCS);
        cloudPixCenters.push_back(center);
      }
    } catch (utl::AugerException& augerException) {
      WARNING(augerException.GetExceptionName() + ", " + augerException.GetMessage());
    }
  }
}


void
Offline2ADST::FillSim(const evt::Event& event, GenShower& gen, const bool saveProfiles)
{
  if (!event.HasSimShower())
    return;

  const auto& simEvent= event.GetSimShower();
  const auto& referenceCS = det::Detector::GetInstance().GetReferenceCoordinateSystem();

  gen.SetEnergy(simEvent.GetEnergy()/eV);
  gen.SetMuonNumber(simEvent.GetMuonNumber());
  gen.SetMuonWeightScale(simEvent.GetMuonWeightScale());
  gen.SetCalorimetricEnergy(simEvent.GetCalorimetricEnergy() / eV);
  gen.SetElectromagneticEnergy(simEvent.GetElectromagneticEnergy() / eV);
  gen.SetNmu(simEvent.GetNmu());
  gen.SetXmaxMu(simEvent.GetXmaxMu() / (g/cm2));

  // Fill (radio) radiation energy
  if (simEvent.HasRadioSimulation()) {
    const auto& radiosim = simEvent.GetRadioSimulation();
    gen.SetRadiationEnergy(radiosim.GetRadiationEnergy());
  }

  if (simEvent.HasGeometry()) {
    const auto& core = simEvent.GetPosition();
    gen.SetCoreSiteCS(ToTVector3(core, referenceCS, meter));

    double northing = 0;
    double easting = 0;
    double altitude = 0;
    try {
      const utl::UTMPoint utm(core, ReferenceEllipsoid::GetWGS84());
      northing = utm.GetNorthing();
      easting = utm.GetEasting();
      altitude = utm.GetHeight();
    } catch (UTMPoint::UTMZoneException& zone_e) {
      ERROR("UTMZoneException: sd rec shower invalid");
      return;
    } catch (UTMPoint::UTMException& e) {
      ERROR("UTMException: sd rec shower invalid");
      return;
    }
    gen.SetCoreUTMCS(TVector3(easting, northing, altitude));

    if (simEvent.HasGeometry()) {

      const auto& axis = -simEvent.GetDirection();
      const auto& localCS = simEvent.GetLocalCoordinateSystem();
      gen.SetAxisCoreCS(ToTVector3(axis, localCS));
      gen.SetAxisSiteCS(ToTVector3(axis, referenceCS));

      gen.SetPrimary(simEvent.GetPrimaryParticle());

      gen.SetCoreTime(simEvent.GetTimeStamp().GetGPSSecond(),
                      simEvent.GetTimeStamp().GetGPSNanoSecond());

      // -------  angles
      otoa::err::CalculateEquatorialCoordinates equCalc;

      const auto& coreTime = simEvent.GetTimeStamp();
      const vector<double> inp({
        ModifiedJulianDate(coreTime) + 2400000.5,
        TimeStamp2GMST(coreTime),
        easting,  // x
        northing,  // y
        altitude,  // z
        axis.GetTheta(localCS),
        axis.GetPhi(localCS)
      });

      const auto equOut = equCalc(inp);
      gen.SetRightAscension(equOut[0]);
      gen.SetDeclination(equOut[1]);

      otoa::err::CalculateGalacticCoordinates galCalc;

      const auto galOut = galCalc(inp);
      gen.SetGalacticLongitude(galOut[0]);
      gen.SetGalacticLatitude(galOut[1]);
    }
  }

  // fill the profile MC info
  if (!simEvent.HasLongitudinalProfile())
    return;

  const TabulatedFunction& electronTrace = simEvent.GetLongitudinalProfile();

  if (!simEvent.HasdEdX())
    return;
  const TabulatedFunction& dEdXTrace = simEvent.GetdEdX();

  unsigned int size = electronTrace.GetNPoints();
  vector<double> depth(size);
  vector<double> electrons(size);
  vector<double> dEdX(size);

  for (unsigned int i = 0; i < size; ++i) {
    depth[i] = dEdXTrace[i].X() / (g/cm2);
    dEdX[i] = dEdXTrace[i].Y() / (PeV/(g/cm2));
    electrons[i] = electronTrace[i].Y();
  }

  // set distance to shower maximum
  const double distanceOfShowerMaximum = simEvent.GetDistanceOfShowerMaximum();
  // cm, sic!, this is a bug, but we keep it for backward compatibility -- DO NOT CHANGE
  // In GenShower::GetGeometricalDistanceToShowerMaximum() the value is converted to meter!
  gen.SetDistanceOfShowerMaximum(distanceOfShowerMaximum / cm);

  gen.SetMagneticFieldStrength(simEvent.GetMagneticFieldStrength());
  gen.SetMagneticFieldInclination(simEvent.GetMagneticFieldInclination());
  gen.SetMagneticFieldDeclination(simEvent.GetMagneticFieldDeclination());

  // maximum from generator
  if (simEvent.HasGHParameters()) {
    gen.SetXmaxGaisserHillas(simEvent.GetGHParameters().GetXMax() / (g/cm2));
    gen.SetNmaxGaisserHillas(simEvent.GetGHParameters().GetNMax() / (g/cm2));
  }

  // maximum from interpolation
  double xMaxConex = 0;
  double xMaxConexdEdX = 0;
  if (QuadraticMaximumInterpolation(depth, dEdX, xMaxConex, xMaxConexdEdX)) {
    gen.SetXmaxInterpolated(xMaxConex);
    gen.SetdEdXmaxInterpolated(xMaxConexdEdX);
  }

  // Get depths and distances to injection/first interaction
  const double X1 = simEvent.GetXFirst();
  gen.SetX1(X1 / (g/cm2));
  const double X0 = simEvent.GetXInject();
  gen.SetX0(X0 / (g/cm2));

  if (simEvent.HasGeometry()) {
    try {
      const double depthStepSize = 5*g/cm2;
      const auto& atmo = det::Detector::GetInstance().GetAtmosphere();
      atmo.InitSlantProfileModel(simEvent.GetPosition(), -simEvent.GetDirection(), depthStepSize);
      const atm::ProfileResult distanceVsSlantDepth = atmo.EvaluateDistanceVsSlantDepth();
      gen.SetD1(distanceVsSlantDepth.Y(X1));
      gen.SetD0(distanceVsSlantDepth.Y(X0));
    } catch (NonExistentComponentException& e) {
      WARNING("Can't load profile model. Core-to-depth distances not filled.");
    }
  }

  if (!saveProfiles)
    return;

  gen.SetDepth(depth);
  gen.SetElectrons(electrons);
  gen.SetEnergyDeposit(dEdX);

  // Fill information about leading particles
  auto& adstGenP = gen.GetParticleTree();
  const auto& partTree = simEvent.GetParticleTree();
  ConvertGenParticle(adstGenP, partTree);
}


void
Offline2ADST::ConvertGenParticle(GenParticle& adstGP, const evt::GenParticle& evtGP)
{
  adstGP.SetParentId(evtGP.GetParentId());
  adstGP.SetParentEnergy(evtGP.GetParentEnergy());
  adstGP.SetEnergyCM(evtGP.GetEnergyCM());
  adstGP.SetKinel(evtGP.GetKinel());
  adstGP.SetMultiplicity(evtGP.GetMultiplicity());
  adstGP.SetTargetMass(evtGP.GetTargetMass());
  adstGP.SetDepth(evtGP.GetDepth());
  adstGP.SetHeight(evtGP.GetHeight());

  unsigned int numberOfDaughterParticles = evtGP.GetDaughterParticles().size();
  adstGP.AddNDaughters(numberOfDaughterParticles);
  for (unsigned i = 0; i < numberOfDaughterParticles; ++i) {
    ConvertGenParticle(adstGP.GetDaughterParticles()[i], evtGP.GetDaughterParticles()[i]);
  }
}