RecDataWriter.cc

Go to the documentation of this file.

#include "RecDataWriter.h"

#include "Offline2ADST.h"

#include "FOVCalculator.h"

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

#include <det/Detector.h>

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

#include <sdet/SDetector.h>
#include <sdet/Station.h>

#include <evt/Event.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerSimData.h>
#include <evt/ShowerFRecData.h>
#include <evt/ShowerSRecData.h>
#include <evt/VGaisserHillasParameter.h>
#include <evt/GaisserHillas4Parameter.h>
#include <evt/GaisserHillas6Parameter.h>
#include <evt/ShowerSimData.h>

#include <fevt/FEvent.h>
#include <fevt/Header.h>
#include <fevt/Eye.h>
#include <fevt/Telescope.h>
#include <fevt/EyeHeader.h>
#include <fevt/EyeRecData.h>
#include <fevt/TelescopeSimData.h>
#include <fevt/TelescopeRecData.h>
#include <fevt/Pixel.h>
#include <fevt/PixelRecData.h>

#include <sevt/SEvent.h>
#include <sevt/Station.h>
#include <sevt/Header.h>
#include <sevt/StationRecData.h>
#include <sevt/StationSimData.h>
#include <sevt/StationTriggerData.h>
#include <sevt/SortCriteria.h>
#include <sevt/PMT.h>
#include <sevt/PMTRecData.h>
#include <sevt/PMTCalibData.h>
#include <sevt/StationConstants.h>

#include <revt/REvent.h>
#include <revt/Station.h>
#include <revt/Channel.h>
#include <revt/EventTrigger.h>
#include <revt/Header.h>
#include <revt/StationGPSData.h>
#include <revt/StationRecData.h>
#include <revt/StationSimData.h>
#include <revt/StationTriggerData.h>

#include <rdet/RDetector.h>

#include <utl/Reader.h>
#include <utl/ErrorLogger.h>
#include <utl/TimeStamp.h>
#include <utl/ModifiedJulianDate.h>
#include <utl/UTCDateTime.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 <utl/config.h>
#include <utl/Particle.h>
#include <utl/CoordinateSystemPtr.h>
#include <utl/SaveCurrentTDirectory.h>
#include <utl/String.h>

#include <atm/ProfileResult.h>
#include <atm/MolecularDB.h>
#include <atm/MolecularZone.h>
#include <atm/Atmosphere.h>
#include <atm/LidarDB.h>
#include <atm/LidarZone.h>

#ifdef INCLUDE_HAS
#  include <has/MuonProfile.h>
#  include <has/VTankResponse.h>
#  include <has/TankResponseFactory.h>
#  include <has/EMComponent.h>
#  include <has/HASUtilities.h>
#endif

#include <boost/tuple/tuple.hpp>

#include <adst/Detector.h>
#include <adst/FileInfo.h>
#include <adst/SDEvent.h>
#include <adst/FDEvent.h>
#include <adst/RecEvent.h>
#include <adst/FdRecApertureLight.h>
#include <adst/FdGenApertureLight.h>
#include <adst/FdRecPixel.h>
#include <adst/FdGenShower.h>
#include <adst/FdRecLevel.h>
#include <adst/SdRecLevel.h>
#include <adst/DetectorGeometry.h>
#include <adst/FdRecLevel.h>
#include <adst/FdLidarData.h>
#include <adst/FdAerosols.h>
#include <adst/Traces.h>

#include <adst/AsciiSDConverter.h>
#include <adst/AsciiHybridConverter.h>
#include <adst/AsciiHybridOnlineConverter.h>

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

#include <cmath>
#include <vector>
#include <sstream>
#include <iostream>
#include <set>

using namespace otoa;
using namespace fwk;
using namespace fdet;
using namespace sdet;
using namespace fevt;
using namespace sevt;
using namespace utl;
using namespace atm;
using namespace std;
#ifdef INCLUDE_HAS
using namespace HASTools;
using namespace HASUtilities;
#endif


namespace RecDataWriterNG {

  RecDataWriter::RecDataWriter() :
    fADSTConverter(new otoa::Offline2ADST)
  { }


  RecDataWriter::~RecDataWriter()
  {
    delete fRecEvent;
    delete fADSTConverter;
  }


  VModule::ResultFlag
  RecDataWriter::Init()
  {
    if (!ReadConfig())
      return eFailure;

    const auto& cfg = fADSTConverter->Config();

    // info output
    ostringstream info;
    info << "Parameters:\n"
            "             output file: " << fRootOutputFileName << "\n"
            "  min ASCII FD rec level: " << fAsciiMinFdRecLevel << "\n"
            "  min ASCII SD rec level: " << fAsciiMinSdRecLevel << "\n"
            "          save FD traces: " << cfg.StoreFDTraces() << "\n"
            "          save SD traces: " << cfg.StoreSDTraces() << "\n"
            "          save MC traces: " << cfg.StoreMCTraces() << "\n"
            "     use weather station: " << cfg.UseWeatherStations() << "\n"
            "    connect to databases: " << cfg.ConnectToDatabases() << "\n"
            "        store Lidar data: " << cfg.StoreLidarData() << "\n"
            " store cloud camera data: " << cfg.StoreCloudCameraData() << "\n"
            "   store GOES cloud data: " << cfg.StoreGOESData() << "\n"
            "        store all pixels: " << cfg.StoreAllPixels() << "\n"
            "  untriggered MC profile: " << (cfg.DropUntriggeredMCProfiles() ? "drop"  :"keep") << "\n"
            "     untriggered MC tels: " << (cfg.DropUntriggeredMCTelescopes() ? "drop"  :"keep") << "\n"
            "    untriggered detector: " << (cfg.DropUntriggeredDetector() ? "drop"  :"keep") << "\n"
            "  re-bin sim. tel traces: ";
    if (cfg.GetRebinSimTelescopeTraces() > 0)
      info << cfg.GetRebinSimTelescopeTraces()/ns << " ns\n";
    else
      info << "NO\n";
    info << "    store photoelectrons: " << cfg.StoreSDPETimeDistribution() << "\n"
            " Mie att. prof. sampling: " << cfg.AerosolAttenuationLengthSampling()/m << " m\n"
            "        Laser wavelength: " << cfg.LaserWavelength()/nanometer << " nm\n"
            "         store particles: " << cfg.StoreSDParticles() << "\n"
            "        Radio Save level: " << cfg.StoreRadioLevel() << "\n"
            "  save MD traces (AMIGA): " << cfg.StoreMDTraces();
    INFO(info);

    if (cfg.AerosolAttenuationLengthSampling() > 0 && cfg.AerosolAttenuationLengthSampling() != 200*m) {
      ERROR("Aerosol attenuation length sampling was not 200 m; this will invoke interpolation "
            "which is not properly implemented, so the aerosol data stored in the ADST will not "
            "be accurate. Please check your configuration and use 200 m!");
    }

    // This radio part just checks that the radio save level is good
    switch (cfg.StoreRadioLevel()) {
    case -1:
      INFO("Radio will not be saved");
      break;
    case 0:
      INFO("Only Radio shower and station information will be saved");
      break;
    case 1:
      INFO("Radio data will be saved at the station level");
      break;
    case 2:
      INFO("Radio data will be saved at the Channel level");
      break;
    default:
      INFO("Correct Values for SaveRadio field are -1/0/1/2! "
           "Defaulting to saving the radio at the shower level! "
           "Check RecDataWriter.xml");
      break;
    }

    return eSuccess;
  }


  VModule::ResultFlag
  RecDataWriter::Run(evt::Event& event)
  {
    // initialize RecEvent and RecEventFile first
    if (!PrepareRecEventAndRootFile())  // checks fWriteROOT!
      return eFailure;

    // initialize AsciiConverter first
    if (!fAsciiConverter && fWriteASCII && !PrepareAsciiFile())
      return eFailure;

    if (!fWriteASCII && !fWriteROOT)
      return eSuccess;

    const bool okay = fADSTConverter->Convert(event, *fRecEvent);
    if (!okay)
      return eSuccess; // skip

    if (fWriteROOT && fRecEventFile)
      fRecEventFile->WriteEvent();

    if (fWriteASCII && fAsciiConverter)
      fAsciiConverter->Convert(*fRecEvent);

    ++fNWritten;

    switch (fRecEvent->GetNEyes()) {
    case 2:
      ++fNWrittenStereo;
      break;
    case 3:
      ++fNWrittenTriple;
      break;
    case 4:
      ++fNWrittenQuadruple;
      break;
    default:
      break;
    }

    return eSuccess;
  }


  VModule::ResultFlag
  RecDataWriter::Finish()
  {
    if (fRecEventFile) {
      fADSTConverter->FinishDetectorAndFileInfo();

      if (fVerbosity > -1)
        INFO("writing detector");
      fRecEventFile->WriteDetectorGeometry(fADSTConverter->GetDetectorGeometry());

      if (fVerbosity > -1)
        INFO("writing file info");
      fRecEventFile->WriteFileInfo(fADSTConverter->GetFileInfo());

      fRecEventFile->Close();
      delete fRecEventFile;
      fRecEventFile = nullptr;

      ostringstream info;
      info << fNWritten << " event" << String::Plural(fNWritten) << " written to ADST ROOT file: "
           << fNWrittenStereo << " stereo | "
           << fNWrittenTriple << " triple | "
           << fNWrittenQuadruple << " quadruple";
      INFO(info);
    }

    // finish ASCII writing
    if (fAsciiConverter) {
      delete fAsciiConverter;
      fAsciiConverter = nullptr;
    }

    return eSuccess;
  }


  bool
  RecDataWriter::ReadConfig()
  {
    const auto cc = CentralConfig::GetInstance();
    auto topB = cc->GetTopBranch("RecDataWriter");

    if (!topB) {
      ERROR("Could not find configuration branch RecDataWriter");
      return eFailure;
    }

    auto& adstCfg = fADSTConverter->Config();

    if (!ReadRootConfig(topB) || !ReadAsciiConfig(topB))
      return false;

    topB.GetChild("verbosity").GetData(fVerbosity);
    adstCfg.SetVerbosity(fVerbosity);

    // configuration of WHAT is written to ADST
    adstCfg.SetWriteFD(topB.GetChild("writeFD").Get<bool>());
    adstCfg.SetWriteSD(topB.GetChild("writeSD").Get<bool>());

    // "Mostly" FD stuff
    adstCfg.SetUseWeatherStations(topB.GetChild("useWeatherStations").Get<bool>());
    adstCfg.SetConnectToDatabases(topB.GetChild("connectToDatabases").Get<bool>());
    adstCfg.SetStoreLidarData(topB.GetChild("storeLidarData").Get<bool>());
    adstCfg.SetStoreCloudCameraData(topB.GetChild("storeCloudCameraData").Get<int>());
    adstCfg.SetStoreGOESData(topB.GetChild("storeGOESData").Get<int>());
    adstCfg.SetStoreAllPixels(topB.GetChild("storeAllPixels").Get<bool>());
    adstCfg.SetStoreFDTraces(topB.GetChild("saveFDTraces").Get<int>());
    adstCfg.SetStoreLCEfficiency(topB.GetChild("saveLCEfficiency").Get<bool>());
    adstCfg.SetVAODReferenceHeight(topB.GetChild("vaodReferenceHeight").Get<double>());
    adstCfg.SetFieldOfViewOption(topB.GetChild("calculateFOVs").Get<int>());
    adstCfg.SetAerosolAttenuationLengthSampling(topB.GetChild("aerosolAttenuationLengthSampling").Get<double>());
    adstCfg.SetLaserWavelength(topB.GetChild("laserWavelength").Get<double>());

    // "mostly" SD stuff
    adstCfg.SetStoreSDTraces(topB.GetChild("saveSDTraces").Get<int>());
    adstCfg.SetStoreCalibHistos(topB.GetChild("storeCalibHistos").Get<int>());
    adstCfg.SetStoreSDTracesMinEnergy(topB.GetChild("minEnergyForTraces").Get<double>());
    adstCfg.SetStoreMPDsMinEnergy(topB.GetChild("minEnergyForStationsMPDs").Get<double>());
    adstCfg.SetStoreSDPETimeDistribution(topB.GetChild("savePETimeDistribution").Get<bool>());
    adstCfg.SetStoreMCTraces(topB.GetChild("saveMCTraces").Get<bool>());
    adstCfg.SetStoreSDParticles(topB.GetChild("saveParticles").Get<bool>());

    auto saveParticleTypesB = topB.GetChild("saveParticleTypes");
    if (saveParticleTypesB) {
      const auto p = topB.GetChild("saveParticleTypes").Get<vector<int>>();
      adstCfg.SetStoreSDParticleTypes(std::set<int>(p.begin(), p.end()));
    }

    // MC stuff
    adstCfg.SetDropUntriggeredMCProfiles(topB.GetChild("dropUntriggeredMCProfiles").Get<bool>());

    adstCfg.SetDropUntriggeredMCTelescopes(topB.GetChild("dropUntriggeredMCTelescopes").Get<bool>());

    adstCfg.SetDropUntriggeredDetector(topB.GetChild("dropUntriggeredDetector").Get<bool>());

    adstCfg.SetRebinSimTelescopeTraces(topB.GetChild("rebinSimTelescopesTraces").Get<double>());

    // Radio stuff
    adstCfg.SetStoreRadioLevel(topB.GetChild("SaveRadio").Get<int>());
    adstCfg.SetStoreExcludedRdStations(topB.GetChild("StoreExcludedRdStations").Get<bool>());

    // Muon detector stuff
    adstCfg.SetStoreMDTraces(topB.GetChild("saveMDTraces").Get<bool>());
    adstCfg.SetStoreMDInjectedParticles(topB.GetChild("saveMDInjectedParticles").Get<bool>());
    auto saveMDParticleTypesB = topB.GetChild("saveMDParticleTypes");
    if (saveMDParticleTypesB) {
      const auto p = saveMDParticleTypesB.Get<vector<int>>();
      adstCfg.SetStoreMDParticleTypes(std::set<int>(p.begin(), p.end()));
    }

    // FIXME move HAS stuff to otoa
  #ifdef INCLUDE_HAS
    const auto hasB = cc->GetTopBranch("SdHorizontalReconstruction");

    if (hasB) {
      const auto muonProfile = hasB.GetChild("MuonProfile").Get<string>();
      auto muonMapBranch = cc->GetTopBranch("MuonProfile").GetChild(muonProfile);
      fMuonProfile = new HASTools::MuonProfile(muonMapBranch);
      const string tmp = "MuonProfile";
      fMuonMapType = muonProfile.substr(0, muonProfile.length() - tmp.length());

      const auto tankResponse = hasB.GetChild("TankResponse").Get<string>();
      auto tankResponseBranch = cc->GetTopBranch("TankResponse").GetChild(tankResponse);
      fTankResponse = &HASTools::TankResponseFactory::GetInstance().GetTankResponse(tankResponseBranch);

      const auto emComponent = hasB.GetChild("EMComponent").Get<string>();
      auto EMComponentBranch = cc->GetTopBranch("EMComponent").GetChild(emComponent);
      fEMComponent = new HASTools::EMComponent(EMComponentBranch);

      const auto muonMapB = topB.GetChild("muonMap");
      muonMapB.GetChild("nrOfPointsPerContour").GetData(fNoOfPoints);
      muonMapB.GetChild("nrOfContours").GetData(fNoOfContours);
      muonMapB.GetChild("contourSignalFactor").GetData(fContourSignalFactor);
    } else {
      INFO("Configuration for SdHorizontalReconstruction not found, no special treatment for horizontal events.");
    }
  #endif
    return true;
  }


  bool
  RecDataWriter::ReadRootConfig(utl::Branch& topB)
  {
    auto outputB = topB.GetChild("rootOutput");
    if (!outputB)
      return false;
    outputB.GetChild("outputFileName").GetData(fRootOutputFileName);
    INFO("ROOT output to " + fRootOutputFileName);
    const auto writeMode = outputB.GetChild("outputFileMode").Get<string>();
    if (writeMode == "eWrite") {
      fRootOutputFileMode = RecEventFile::eWrite;
      fWriteROOT = true;
    } else if (writeMode == "eAppend")  {
      fRootOutputFileMode = RecEventFile::eAppend;
      fWriteROOT = true;
    } else if (writeMode == "eNone")  {
      fRootOutputFileMode = RecEventFile::eWrite;
      fWriteROOT = false;
    } else {
      ostringstream err;
      err << "Unrecognized file mode \"" << writeMode << "\"";
      ERROR(err);
      throw utl::NonExistentComponentException();
    }

    return true;
  }


  bool
  RecDataWriter::ReadAsciiConfig(utl::Branch& topB)
  {
    auto aOutputB = topB.GetChild("asciiOutput");
    if (!aOutputB)
      return false;

    const auto aWriteMode = aOutputB.GetChild("outputFileMode").Get<string>();
    if (aWriteMode == "eWrite") {
      fAsciiOutputFileMode = RecEventFile::eWrite;
      fWriteASCII = true;
    } else if (aWriteMode == "eAppend")  {
      fAsciiOutputFileMode = RecEventFile::eAppend;
      fWriteASCII = true;
    } else if (aWriteMode == "eNone") {
      fAsciiOutputFileMode = RecEventFile::eWrite;
      fWriteASCII = false;
    } else {
      ostringstream err;
      err << "Unrecognized ASCII file mode \"" << aWriteMode << "\"";
      ERROR(err);
      throw utl::NonExistentComponentException();
    }
    aOutputB.GetChild("mode").GetData(fAsciiDataMode);

    if (fWriteASCII) {
      aOutputB.GetChild("minimumFDRecLevel").GetData(fAsciiMinFdRecLevel);
      aOutputB.GetChild("minimumSDRecLevel").GetData(fAsciiMinSdRecLevel);
      aOutputB.GetChild("outputFileName").GetData(fAsciiOutputFileName);
      INFO("Ascii output to " + fAsciiOutputFileName);
    }

    return true;
  }


  bool
  RecDataWriter::PrepareRecEventAndRootFile()
  {
    const SaveCurrentTDirectory save;
    if (!fRecEvent)
      fRecEvent = new RecEvent();
    else
      *fRecEvent = RecEvent();

    if (fWriteROOT && !fRecEventFile) {
      // ---- open the RecEventFile and link the buffers ------
      fRecEventFile = new RecEventFile(fRootOutputFileName, fRootOutputFileMode);
      if (!fRecEventFile)
        return false;
      fRecEventFile->SetBuffers(&fRecEvent);
    }

    return true;
  }


  bool
  RecDataWriter::PrepareAsciiFile()
  {
    const auto openMode = ios::out | (fAsciiOutputFileMode == RecEventFile::eAppend ? ios::app : ios::trunc);
    const auto asciiOut = new ofstream(fAsciiOutputFileName, openMode);
    if (!asciiOut->is_open()) {
      ERROR("Could not open ascii output file");
      return false;
    }

    if (fAsciiDataMode == "SD") {
      const auto cv = new ADST::AsciiSDConverter(asciiOut);
      cv->SetMinSdRecLevel(fAsciiMinSdRecLevel);
      fAsciiConverter = (ADST::AsciiConverter*)cv;
    } else if (fAsciiDataMode == "Hybrid") {
      const auto cv = new ADST::AsciiHybridConverter(asciiOut);
      ((ADST::AsciiFDConverter*)cv)->SetMinFdRecLevel(fAsciiMinFdRecLevel);
      fAsciiConverter = (ADST::AsciiConverter*)cv;
    } else if (fAsciiDataMode == "HybridOnline") {
      const auto cv = new ADST::AsciiHybridOnlineConverter(asciiOut);
      ((ADST::AsciiFDConverter*)cv)->SetMinFdRecLevel(fAsciiMinFdRecLevel);
      fAsciiConverter = (ADST::AsciiConverter*)cv;
    }
    // transfer ownership of output to the converter
    fAsciiConverter->SetOwnershipOfStream(true);

    if (fAsciiOutputFileMode == RecEventFile::eWrite)
      fAsciiConverter->WriteFileHeader();

    return true;
  }

}