MD2ADST.cc

Go to the documentation of this file.

#include <utl/config.h>

#include "MD2ADST.h"
#include "Config.h"
#include "ConversionUtil.h"

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

#include <det/Detector.h>

#include <mdet/MDetector.h>
#include <mdet/Counter.h>
#include <mdet/Scintillator.h>

#include <evt/Event.h>
#include <evt/ShowerSimData.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerSRecData.h>

#include <mevt/MEvent.h>
#include <mevt/Counter.h>
#include <mevt/CounterSimData.h>

#include <utl/Branch.h>
#include <utl/ErrorLogger.h>
#include <utl/UTCDateTime.h>
#include <utl/TimeStamp.h>
#include <utl/CoordinateSystemPtr.h>
#include <utl/ReferenceEllipsoid.h>
#include <utl/Point.h>
#include <utl/Vector.h>
#include <utl/UTMPoint.h>
#include <utl/Particle.h>
#include <utl/AugerUnits.h>
#include <utl/MathConstants.h>
#include <utl/PhysicalConstants.h>
#include <utl/PhysicalFunctions.h>
#include <utl/Is.h>

#include <iostream>
#include <string>

using utl::m;
using utl::meter;
using std::cout;
using std::endl;


namespace otoa {

  MD2ADST::MD2ADST(const otoa::Config& config) :
    fConfig(config),
    fDetector(&det::Detector::GetInstance()),
    fMdDetector(&fDetector->GetMDetector()),
    fReferenceCS(fDetector->GetReferenceCoordinateSystem())
  { }


  void
  MD2ADST::Convert(const evt::Event& inEvent, RecEvent& outEvent)
    const
  {
    const auto& inMdEvent = inEvent.GetMEvent();
    auto& outMdEvent = outEvent.GetMDEvent();

    // Fill event level data
    FillMEvent(inEvent, outMdEvent);

    // Fill shower data
    utl::CoordinateSystemPtr recShowerPlaneCS;
    if (inEvent.HasRecShower() && inEvent.GetRecShower().HasMRecShower()) {

      const auto& inMdShower = inEvent.GetRecShower().GetMRecShower();
      auto& outMdShower = outMdEvent.GetMdRecShower();

      // Fill the reconstruction flags
      outMdShower.SetLdfReconstructed(inMdShower.IsLdfReconstructed());
      outMdShower.SetBetaFixed(inMdShower.IsBetaFixed());
      outMdShower.SetCoreFixedLdf(inMdShower.IsCoreFixedLdf());
      outMdShower.SetGeometryReconstructed(inMdShower.IsGeometryReconstructed());
      outMdShower.SetCoreFixedGeo(inMdShower.IsCoreFixedGeo());
      outMdShower.SetCurvatureFixed(inMdShower.IsCurvatureFixed());

      FillRecShower(inMdShower, outMdShower);
      recShowerPlaneCS = GetShowerPlaneCS(inMdShower);

    } else {

      WARNING("Distances from the MD counters to the core will not be set because the event was not reconstructed with the MD.");

    }

    // Add counters
    for (auto inCounter = inMdEvent.CountersBegin(); inCounter != inMdEvent.CountersEnd(); ++inCounter)
      outMdEvent.AddCounter(MakeCounter(*inCounter, recShowerPlaneCS));

    // Add simulated data
    if (inEvent.HasSimShower()) {
      const auto& showerPlaneCS = inEvent.GetSimShower().GetShowerCoordinateSystem();
      for (auto inCounter = inMdEvent.CountersBegin(); inCounter != inMdEvent.CountersEnd(); ++inCounter) {
        if (inCounter->HasSimData())
          outMdEvent.AddSimCounter(MakeSimCounter(*inCounter, showerPlaneCS));
      }
    }
  }


  void
  MD2ADST::FillMEvent(const evt::Event& inEvent, MDEvent& outMdEvent)
    const
  {
    // Setting the number of candidate counters
    const auto& inMdEvent = inEvent.GetMEvent();
    const unsigned int nCandidates = inMdEvent.GetNumberOfCandidateCounters();
    outMdEvent.SetNumberOfCandidates(nCandidates);

    if (inEvent.HasRecShower() && inEvent.GetRecShower().HasMRecShower()) {
      const auto& inMdShower = inEvent.GetRecShower().GetMRecShower();
      if (inMdShower.HasMLDF()) {
        outMdEvent.SetRecLevel(eHasMdLDF);
      }
    } else if (nCandidates > 0) {
      outMdEvent.SetRecLevel(eHasTriggeredMdStations);
    } else {
      outMdEvent.SetRecLevel(eNoMdEvent);
    }
  }


  void
  MD2ADST::FillRecShower(const evt::ShowerMRecData& inMdShower, MdRecShower& outMdShower)
    const
  {
    // Fill the shower axis (in the core coordinate system)
    const auto& corePosition = inMdShower.GetCorePosition();
    const auto& coreCS = fwk::LocalCoordinateSystem::Create(corePosition);
    const auto& axis1 = inMdShower.GetAxis();
    TVector3 axis2(0,0,1);
    axis2.SetTheta(axis1.GetTheta(coreCS));
    axis2.SetPhi(axis1.GetPhi(coreCS));
    outMdShower.SetAxisCoreCS(axis2);
    outMdShower.SetMdSdAngle(inMdShower.GetMdSdAngle());

    // Fill the error info of the geometrical reconstruction
    outMdShower.SetZenithError(inMdShower.GetThetaError());
    outMdShower.SetAzimuthError(inMdShower.GetPhiError());
    outMdShower.SetZenithAzimuthCorrelation(inMdShower.GetThetaPhiCorrelation());
    outMdShower.SetAngleChi2(inMdShower.GetAngleChi2(), inMdShower.GetAngleNdof());
    outMdShower.SetTimeResidualMean(inMdShower.GetTimeResidualMean());
    outMdShower.SetTimeResidualSpread(inMdShower.GetTimeResidualSpread());

    const double radius = inMdShower.GetCurvature();
    const double dradius = inMdShower.GetCurvatureError();
    outMdShower.SetCurvatureRadius(radius, dradius);

    // Fill the LDF
    if (inMdShower.HasMLDF()) {
      const auto& ldfTable = inMdShower.GetMLDF();
      const auto n = ldfTable.GetNPoints();
      std::vector<double> rhos;
      rhos.reserve(n);
      std::vector<double> values;
      values.reserve(n);
      for (size_t i = 0; i < n; i += 2) {
        const double x = ldfTable[i].X();
        const double y = ldfTable[i].Y();
        rhos.push_back(x);
        values.push_back(y);
        if (x > 1500)
          i += 4;
      }
      outMdShower.SetTabulatedValues(rhos, values);

      // Fill LDF parameters
      outMdShower.SetReferenceDistance(inMdShower.GetReferenceDistance());
      outMdShower.SetNMuRef(inMdShower.GetNMuRef(), inMdShower.GetNMuRefError());
      outMdShower.SetNMuRefSystematics(inMdShower.GetNMuRefSystematics());
      outMdShower.SetBeta(inMdShower.GetBeta(), inMdShower.GetBetaError());
      outMdShower.SetBetaSystematics(inMdShower.GetBetaSystematics());

      // Fill goodness of fit parameters
      outMdShower.SetMLDFNdof(inMdShower.GetMLDFNdof());
      outMdShower.SetMLDFLikelihood(inMdShower.GetMLDFLikelihood());
      outMdShower.SetMLDFChi2(inMdShower.GetMLDFChi2());

      // Fill the core position
      outMdShower.SetCoreSiteCS(ToTVector3(corePosition, fReferenceCS, meter));

      //Fill the core times
      outMdShower.SetCoreTime(inMdShower.GetCoreTime().GetGPSSecond(),
                              inMdShower.GetCoreTime().GetGPSNanoSecond());

      // Fill the core error
      const utl::Vector& coreError = inMdShower.GetCoreError();
      outMdShower.SetCoreNorthingError(coreError.GetY(fReferenceCS));
      outMdShower.SetCoreEastingError(coreError.GetX(fReferenceCS));
      outMdShower.SetCoreNorthingEastingCorrelation(inMdShower.GetCorrelationXY());
    }
  }


  MdRecCounter
  MD2ADST::MakeCounter(const mevt::Counter& inCounter, const utl::CoordinateSystemPtr& recShowerPlaneCS)
    const
  {
    MdRecCounter outCounter;

    // set ids
    const unsigned int counterId = inCounter.GetId();
    outCounter.SetId(counterId);

    const auto& detCounter = fMdDetector->GetCounter(counterId);

    // get the counter for future master data lookups
    outCounter.SetSdPartnerId(detCounter.GetAssociatedTankId());

    // map the status of counter to outCounter
    if (inCounter.IsCandidate()) {
      outCounter.SetCandidate();
    } else if (inCounter.IsSilent()) {
      outCounter.SetSilent();
    } else if (inCounter.IsRejected()) {
      outCounter.SetRejected();
    }

    // set the reconstructed data flag
    const bool recDataFlag = inCounter.HasRecData();
    outCounter.SetRecDataFlag(recDataFlag);
    if (!recDataFlag) {  // if not reconstructed data
      return outCounter;
    }

    // set dense flag
    try {
      if (detCounter.IsDense())
        outCounter.SetIsDense();
    } catch (utl::AugerException& e) {
      ERROR(e.GetMessage());
      exit(-1);
    }

    outCounter.SetNumberOfChannelsOn(inCounter.GetNumberOfChannelsOn());

    outCounter.SetNumberOfMuons(inCounter.GetNumberOfEstimatedMuons());
    outCounter.SetNumberOfMuonsErrorHigh(inCounter.GetNumberOfMuonsErrorHigh());
    outCounter.SetNumberOfMuonsErrorLow(inCounter.GetNumberOfMuonsErrorLow());
    outCounter.SetNumberOfMuonsLowLimit(inCounter.GetNumberOfMuonsLowLimit());

    outCounter.SetMeanMuons(inCounter.GetMeanMuons());
    outCounter.SetMeanMuonsErrorHigh(inCounter.GetMeanMuonsErrorHigh());
    outCounter.SetMeanMuonsErrorLow(inCounter.GetMeanMuonsErrorLow());
    outCounter.SetMeanMuonsLowLimit(inCounter.GetMeanMuonsLowLimit());

    outCounter.SetActiveArea(inCounter.GetActiveArea());

    outCounter.SetMuonDensity(inCounter.GetMuonDensity());
    outCounter.SetMuonDensityErrorHigh(inCounter.GetMuonDensityErrorHigh());
    outCounter.SetMuonDensityErrorLow(inCounter.GetMuonDensityErrorLow());

    outCounter.SetMeanMuonDensity(inCounter.GetMeanMuonDensity());
    outCounter.SetMeanMuonDensityErrorHigh(inCounter.GetMeanMuonDensityErrorHigh());
    outCounter.SetMeanMuonDensityErrorLow(inCounter.GetMeanMuonDensityErrorLow());

    outCounter.SetLDFResidual(inCounter.GetLDFResidual());
    outCounter.SetTimeResidual(inCounter.GetTimeResidual());
    outCounter.SetPlaneFrontDelay(inCounter.GetPlaneFrontDelay());

    // ADC LG
    outCounter.SetADCCalibratedLG(inCounter.IsADCCalibratedLG());
    if (outCounter.IsADCCalibratedLG()) {
      outCounter.SetNumberOfMuonsLG(inCounter.GetNumberOfEstimatedMuonsLG());
      outCounter.SetNumberOfMuonsErrorHighLG(inCounter.GetNumberOfMuonsErrorHighLG());
      outCounter.SetNumberOfMuonsErrorLowLG(inCounter.GetNumberOfMuonsErrorLowLG());
      outCounter.SetActiveAreaLG(inCounter.GetActiveAreaLG());
      outCounter.SetMuonDensityLG(inCounter.GetMuonDensityLG());
      outCounter.SetMuonDensityErrorHighLG(inCounter.GetMuonDensityErrorHighLG());
      outCounter.SetMuonDensityErrorLowLG(inCounter.GetMuonDensityErrorLowLG());
    }

    // ADC HG
    outCounter.SetADCCalibratedHG(inCounter.IsADCCalibratedHG());
    if (outCounter.IsADCCalibratedHG()) {
      outCounter.SetNumberOfMuonsHG(inCounter.GetNumberOfEstimatedMuonsHG());
      outCounter.SetNumberOfMuonsErrorHighHG(inCounter.GetNumberOfMuonsErrorHighHG());
      outCounter.SetNumberOfMuonsErrorLowHG(inCounter.GetNumberOfMuonsErrorLowHG());
      outCounter.SetActiveAreaHG(inCounter.GetActiveAreaHG());
      outCounter.SetMuonDensityHG(inCounter.GetMuonDensityHG());
      outCounter.SetMuonDensityErrorHighHG(inCounter.GetMuonDensityErrorHighHG());
      outCounter.SetMuonDensityErrorLowHG(inCounter.GetMuonDensityErrorLowHG());
    }

    // Fill the counter position
    const auto& counterPosition = detCounter.GetPosition();
    outCounter.SetPosition(ToTVector3(counterPosition, fReferenceCS, meter));

    // Fill the position in the shower plane reference system
    if (recShowerPlaneCS) {  // valid check of a boost::shared_ptr!
      const double rho = counterPosition.GetRho(recShowerPlaneCS) / m;
      const double azimuth = counterPosition.GetPhi(recShowerPlaneCS) / utl::deg;
      const double delta = counterPosition.GetZ(recShowerPlaneCS) / m;
      outCounter.SetSPDistance(rho, 0);   // errors are set to zero...
      outCounter.SetSPAzimuth(azimuth, 0);
      outCounter.SetSPDelta(delta, 0);
    }

    if (inCounter.HasT50()) {
      const auto signalT50 = inCounter.GetSignalT50();
      outCounter.SetT50(signalT50.GetGPSSecond(), signalT50.GetGPSNanoSecond());
      outCounter.SetT50Error(inCounter.GetT50Error());
      outCounter.SetT502SdStart(inCounter.GetT502SdStart());
    }

    // add the modules
    for (auto inModule = inCounter.ModulesBegin(); inModule != inCounter.ModulesEnd(); ++inModule) {
      const unsigned int moduleId = inModule->GetId();
      const MdRecModule outModule = MakeModule(*inModule, detCounter.GetModule(moduleId), recShowerPlaneCS);
      outCounter.AddModule(outModule);
    }

    return outCounter;
  }


  MdRecModule
  MD2ADST::MakeModule(const mevt::Module& inModule, const mdet::Module& mdetModule, const utl::CoordinateSystemPtr& recShowerPlaneCS)
    const
  {
    MdRecModule outModule;
    const unsigned int moduleId = inModule.GetId();
    outModule.SetId(moduleId);

    // set flags
    if (inModule.IsCandidate()) {
      outModule.SetCandidate();
    } else if (inModule.IsSilent()) {
      outModule.SetSilent();
    } else if (inModule.IsRejected()) {
      outModule.SetRejected();
    }

    // return if no reconstructed data
    if (!inModule.HasRecData())
      return outModule;

    // beware some data is at the module and other at the module reconstructed data levels
    const auto& inRecData = inModule.GetRecData();

    // Set the master data
    outModule.SetActiveArea(inRecData.GetActiveArea());
    outModule.SetSegmentation(inRecData.GetSegmentation());

    double sampleNanos = 0;
    int numberOfSamples = 0;

    if (!mdetModule.IsSiPM()) {
      sampleNanos = mdetModule.GetFrontEnd().GetMeanSampleRatePeriod();
      numberOfSamples = mdetModule.GetFrontEnd().GetBufferLength();
    } else {
      sampleNanos = mdetModule.GetFrontEndSiPM().GetMeanSampleRatePeriod();
      numberOfSamples = mdetModule.GetFrontEndSiPM().GetBufferLength();
    }

    outModule.SetSampleTime(sampleNanos);
    outModule.SetNumberOfSamples(numberOfSamples);
    //outModule.SetSaturationFlag(inRecData.SaturationFlag());

    // Fill the channelsOn trace
    if (inRecData.HasChannelsOn()) {
      const auto& inChannelsOn = inRecData.GetChannelsOn();
      outModule.GetChannelsOn().assign(inChannelsOn.Begin(), inChannelsOn.End());
      outModule.SetChannelsOnBinning(inChannelsOn.GetBinning());
    }

    if (inRecData.HasChannelsOnStartTime()) {
      const auto& t = inRecData.GetChannelsOnStartTime();
      outModule.SetChannelsOnStartTime(t.GetGPSSecond(), t.GetGPSNanoSecond());
    }

    // Fill the pattern matches vector
    outModule.GetPatternMatchTimes() = inRecData.GetPatternMatchTimes();

    // Fill the nmuVsTime and muVsTime traces
    if (inRecData.HasMeanMuonsVsTime()) {
      const auto& inMeanMuonsVsTime = inRecData.GetMeanMuonsVsTime();
      outModule.SetMeanMuonsVsTime(inMeanMuonsVsTime.Begin(), inMeanMuonsVsTime.End());
    }
    if (inRecData.HasNumberOfMuonsVsTime()) {
      const auto& inNumberOfMuonsVsTime = inRecData.GetNumberOfMuonsVsTime();
      outModule.SetNumberOfMuonsVsTime(inNumberOfMuonsVsTime.Begin(), inNumberOfMuonsVsTime.End());
    }

    // Fill the number of muons data
    outModule.SetNumberOfEstimatedMuons(inRecData.GetNumberOfEstimatedMuons());
    outModule.SetNumberOfMuonsErrorHigh(inRecData.GetNumberOfMuonsErrorHigh());
    outModule.SetNumberOfMuonsErrorLow(inRecData.GetNumberOfMuonsErrorLow());
    outModule.SetNumberOfMuonsLowLimit(inRecData.GetNumberOfMuonsLowLimit());

    outModule.SetMeanMuons(inRecData.GetMeanMuons());
    outModule.SetMeanMuonsErrorHigh(inRecData.GetMeanMuonsErrorHigh());
    outModule.SetMeanMuonsErrorLow(inRecData.GetMeanMuonsErrorLow());
    outModule.SetMeanMuonsLowLimit(inRecData.GetMeanMuonsLowLimit());

    outModule.SetMuonDensity(inRecData.GetMuonDensity());
    outModule.SetMuonDensityErrorHigh(inRecData.GetMuonDensityErrorHigh());
    outModule.SetMuonDensityErrorLow(inRecData.GetMuonDensityErrorLow());

    outModule.SetMeanMuonDensity(inRecData.GetMeanMuonDensity());
    outModule.SetMeanMuonDensityErrorHigh(inRecData.GetMeanMuonDensityErrorHigh());
    outModule.SetMeanMuonDensityErrorLow(inRecData.GetMeanMuonDensityErrorLow());

    outModule.SetLDFResidual(inRecData.GetLDFResidual());

    // Fill the module position
    const auto& modulePosition = mdetModule.GetPosition();
    outModule.SetPosition(ToTVector3(modulePosition, fReferenceCS, meter));

    // Fill the position in the shower plane reference system
    if (recShowerPlaneCS) {
      const double rho = modulePosition.GetRho(recShowerPlaneCS) / m;
      const double azimuth = modulePosition.GetPhi(recShowerPlaneCS) / utl::deg;
      const double delta = modulePosition.GetZ(recShowerPlaneCS) / m;
      outModule.SetSPDistance(rho, 0);   // errors are set to zero...
      outModule.SetSPAzimuth(azimuth, 0);
      outModule.SetSPDelta(delta, 0);
    }

    // Store integrator traces
    if (fConfig.StoreMDTraces()) {  // if config set to save traces
      if (inModule.HasDynodeTrace()) {
        const auto& t = inModule.GetDynodeTrace();
        outModule.SetDynodeTrace(std::vector<unsigned short int>(t.Begin(), t.End()));
      }

      // SiPM Integrator
      if (inModule.HasIntegratorATrace()) { // if has integrator LG trace
        const auto& t = inModule.GetIntegratorATrace();
        outModule.SetIntegratorATrace(std::vector<unsigned short int>(t.Begin(), t.End()));
      }

      if (inModule.HasIntegratorBTrace()) { // if has integrator HG trace
        const auto& t = inModule.GetIntegratorBTrace();
        outModule.SetIntegratorBTrace(std::vector<unsigned short int>(t.Begin(), t.End()));
      }
    }

    // store integrator reconstruction

    if (inRecData.IsADCCalibratedLG()) {
      outModule.SetMeanChargeMuonLG(inRecData.GetMeanChargeMuonLG());
      outModule.SetStdDevChargeMuonLG(inRecData.GetStdDevChargeMuonLG());

      outModule.SetNumberOfEstimatedMuonsADCLG(inRecData.GetNumberOfEstimatedMuonsADCLG());
      outModule.SetNumberOfMuonsADCErrorHighLG(inRecData.GetNumberOfMuonsADCErrorHighLG());
      outModule.SetNumberOfMuonsADCErrorLowLG(inRecData.GetNumberOfMuonsADCErrorLowLG());

      outModule.SetMuonDensityLG(inRecData.GetMuonDensityLG());
      outModule.SetMuonDensityErrorHighLG(inRecData.GetMuonDensityErrorHighLG());
      outModule.SetMuonDensityErrorLowLG(inRecData.GetMuonDensityErrorLowLG());
      outModule.SetTotalChargeLG(inRecData.GetTotalChargeLG());
    }

    if (inRecData.IsADCCalibratedHG()) {
      outModule.SetMeanChargeMuonHG(inRecData.GetMeanChargeMuonHG());
      outModule.SetStdDevChargeMuonHG(inRecData.GetStdDevChargeMuonHG());

      outModule.SetNumberOfEstimatedMuonsADCHG(inRecData.GetNumberOfEstimatedMuonsADCHG());
      outModule.SetNumberOfMuonsADCErrorHighHG(inRecData.GetNumberOfMuonsADCErrorHighHG());
      outModule.SetNumberOfMuonsADCErrorLowHG(inRecData.GetNumberOfMuonsADCErrorLowHG());

      outModule.SetMuonDensityHG(inRecData.GetMuonDensityHG());
      outModule.SetMuonDensityErrorHighHG(inRecData.GetMuonDensityErrorHighHG());
      outModule.SetMuonDensityErrorLowHG(inRecData.GetMuonDensityErrorLowHG());
      outModule.SetTotalChargeHG(inRecData.GetTotalChargeHG());
    }

    // add channels
    for (auto inChannel = inModule.ChannelsBegin(); inChannel != inModule.ChannelsEnd(); ++inChannel) {
      const auto outChannel = MakeChannel(*inChannel, mdetModule);
      outModule.AddChannel(outChannel);
    }

    return outModule;
  }


  MdRecChannel
  MD2ADST::MakeChannel(const mevt::Channel& inChannel, const mdet::Module& detModule)
    const
  {
    MdRecChannel outChannel;

    // set ids
    const unsigned int channelId = inChannel.GetId();
    outChannel.SetId(channelId);
    const unsigned int scintillatorId = detModule.ChannelToScintillatorId(channelId);
    outChannel.SetScintillatorId(scintillatorId);
    const unsigned int pixelId = detModule.ChannelToPixelId(channelId);
    outChannel.SetPixelId(pixelId);

    // set channel mask
    outChannel.SetMask(inChannel.IsMasked());

    // set trace info
    if (inChannel.HasTrace()) {
      const auto& trace = inChannel.GetTrace();
      outChannel.SetNumberOfBins(trace.GetSize());
      outChannel.SetBinning(trace.GetBinning());

      if (inChannel.HasTraceStartTime()) {
        const auto& t = inChannel.GetTraceStartTime();
        outChannel.SetTraceStartTime(t.GetGPSSecond(), t.GetGPSNanoSecond());
      }

      if (fConfig.StoreMDTraces()) {
        // if config set to save traces
        outChannel.SetTrace(std::vector<char>(trace.Begin(), trace.End()));
      }

    }

    if (inChannel.HasRecData()) {
      const auto& channelRecData = inChannel.GetRecData();
      // Fill the pattern matches vector
      outChannel.GetPatternMatchTimes().assign(channelRecData.PatternMatchBinsBegin(), channelRecData.PatternMatchBinsEnd());
    }

    return outChannel;
  }


  MdSimCounter
  MD2ADST::MakeSimCounter(const mevt::Counter& inCounter, const utl::CoordinateSystemPtr& showerPlaneCS)
    const
  {
    MdSimCounter outCounter;

    const auto& sPos = fMdDetector->GetCounter(inCounter.GetId()).GetPosition();
    const auto& detCounter = fMdDetector->GetCounter(inCounter.GetId());
    const utl::CoordinateSystemPtr& counterCS = detCounter.GetLocalCoordinateSystem();
    const auto& cSimData = inCounter.GetSimData();

    outCounter.SetId(inCounter.GetId());
    outCounter.SetSdPartnerId(detCounter.GetAssociatedTankId());
    outCounter.SetInsideRMinFlag(cSimData.IsInsideMinRadius());

    outCounter.SetSPDelta(sPos.GetZ(showerPlaneCS) / m);
    outCounter.SetSPDistance(sPos.GetRho(showerPlaneCS) / m);
    outCounter.SetSPAzimuth(sPos.GetPhi(showerPlaneCS) / utl::deg);

    for (auto mIt = inCounter.ModulesBegin(); mIt != inCounter.ModulesEnd(); ++mIt) {
      const int mId = mIt->GetId();
      const auto& detModule = detCounter.GetModule(mId);

      // get area of each module
      const double areaPerModule = mIt->HasRecData() ?
        mIt->GetRecData().GetActiveArea() : detModule.GetArea();

      outCounter.SetAreaByModule(areaPerModule, mIt->GetId());

      for (auto chIt = mIt->ChannelsBegin(); chIt != mIt->ChannelsEnd(); ++chIt) {
        const int chId = chIt->GetId();
        const int scId = detModule.ChannelToScintillatorId(chId);
        if (mIt->HasScintillator(scId)) {
          const auto& scintillator = mIt->GetScintillator(scId);
          if (scintillator.HasSimData()) {

            MdSimScintillator myscint;
            const auto& scintSimData = scintillator.GetSimData();
            const int injmu = scintSimData.GetNumberOfInjectedMuons();
            const int ccmu = scintSimData.GetNumberOfCornerClippingMuons();
            const int electrons = scintSimData.GetNumberOfElectrons();
            myscint.SetNumberOfInjectedMuons(injmu);
            myscint.SetNumberOfCornerClippingMuons(ccmu);
            myscint.SetNumberOfElectrons(electrons);
            myscint.SetId(scId);
            myscint.SetChannelId(chId);
            myscint.SetPixelId(detModule.ChannelToPixelId(chId));
            myscint.SetModuleId(mId);

            for (auto pIt = scintSimData.SPEPulsesBegin(); pIt != scintSimData.SPEPulsesEnd(); ++pIt) {
              myscint.AddSPEPulse(pIt->fMu, pIt->fSigma, pIt->fAmplitude, pIt->fDestPixelId);
            }
            for (auto pIt = scintSimData.PEPulsesBegin(); pIt != scintSimData.PEPulsesEnd(); ++pIt) {
              myscint.AddPEPulse(pIt->fT0, pIt->fAmplitude1, pIt->fAmplitude2, pIt->fAmplitude3, pIt->fRiseTime,
                                 pIt->fFallTime1, pIt->fFallTime2, pIt->fFallTime3, pIt->fDestPixelId);
            }

            double sumKinEnergy = 0;
            double sumKinEnergyMuons = 0;
            {
              const auto& pTypes = fConfig.StoreMDParticleTypes();
              for (auto pIt = scintSimData.ParticlesBegin(); pIt != scintSimData.ParticlesEnd(); ++pIt) {
                const double kinEnergy = pIt->GetKineticEnergy();
                sumKinEnergy += kinEnergy;
                const auto pType = pIt->GetType();
                const bool isMuon = (pType == utl::Particle::eMuon || pType == utl::Particle::eAntiMuon);
                if (isMuon)
                  sumKinEnergyMuons += kinEnergy;
                if (fConfig.StoreMDInjectedParticles() &&
                    (pTypes.empty() || utl::Is(pType).In(pTypes))) {
                  AddParticle(*pIt, detCounter, mId, scId, myscint);
                }
              }
            }

            myscint.SetKineticEnergy(sumKinEnergy);
            myscint.SetMuonKineticEnergy(sumKinEnergyMuons);

            myscint.SetEnergyDeposit(scintSimData.GetEnergyDeposit());
            myscint.SetEnergyDepositMuons(scintSimData.GetEnergyDepositMuons());

            outCounter.AddSimScintillator(myscint);
          }
        }
      }
    }

    return outCounter;
  }


  utl::CoordinateSystemPtr
  MD2ADST::GetShowerPlaneCS(const evt::ShowerMRecData& mRecShower)
    const
  {
    //const auto& rCore = mRecShower.GetCorePosition();  // core from MD rec. If no MD geometry reco, then it is core from SD rec
    const auto& rCore = mRecShower.GetReferenceCorePosition();  // core from SD reco
    const auto& coreCS = fwk::LocalCoordinateSystem::Create(rCore);

    //const utl::Vector& rAxis = mRecShower.GetAxis();  // axis from MD geometry reco. Currently deprecated so it outputs an empty vector
    const auto& rAxis = mRecShower.GetReferenceAxis();  // axis from SD reco
    if (!rAxis)
      throw utl::NonExistentComponentException("No axis is provided to compute shower plane CS.");

    const double theta = rAxis.GetTheta(coreCS);
    const double phi = rAxis.GetPhi(coreCS);

    const auto& recShowerPlaneCS =
      utl::CoordinateSystem::RotationY(theta, utl::CoordinateSystem::RotationZ(phi, coreCS));
    return recShowerPlaneCS;
  }


  void
  MD2ADST::AddParticle(const utl::Particle& p, const mdet::Counter& detCounter, const int mId, const int scintId, MdSimScintillator& myScint)
    const
  {
    const utl::CoordinateSystemPtr& counterCS = detCounter.GetLocalCoordinateSystem();
    const auto& detModule = detCounter.GetModule(mId);
    const utl::CoordinateSystemPtr& moduleCS = detModule.GetLocalCoordinateSystem();
    const auto& detScint = detModule.GetScintillator(scintId);
    const utl::CoordinateSystemPtr& scintCS = detScint.GetLocalCoordinateSystem();

    const utl::Point& pos = p.GetPosition();
    const auto pc = pos.GetCoordinates(counterCS);
    const TVector3 positionStationCS(pc.get<0>(), pc.get<1>(), pc.get<2>());
    const auto pm = pos.GetCoordinates(moduleCS);
    const TVector3 positionModuleCS(pm.get<0>(), pm.get<1>(), pm.get<2>());
    const auto ps = pos.GetCoordinates(scintCS);
    const TVector3 positionScintillatorCS(ps.get<0>(), ps.get<1>(), ps.get<2>());
    const utl::Vector& direction = p.GetDirection();
    const double zenithStationCS = direction.GetTheta(counterCS);
    const double azimuthStationCS = direction.GetPhi(counterCS);
    const double azimuthModuleCS = direction.GetPhi(moduleCS);
    const double azimuthScintCS = direction.GetPhi(scintCS);

    //const utl::Particle::Source source = partIt->GetSource();
    const short int source = p.GetSource();
    const double kineticEnergy = p.GetKineticEnergy();
    const int type = p.GetType();

    myScint.AddParticle(zenithStationCS, azimuthStationCS,
                        azimuthModuleCS, azimuthScintCS,
                        positionStationCS, positionModuleCS, positionScintillatorCS,
                        kineticEnergy, source, type);
  }

}