DataProviderAlg.h

Go to the documentation of this file.

// Class:       PointIdAlg
// Authors:     D.Stefan (Dorota.Stefan@ncbj.gov.pl),         from DUNE, CERN/NCBJ, since May 2016
//              R.Sulej (Robert.Sulej@cern.ch),               from DUNE, FNAL/NCBJ, since May 2016
//              P.Plonski,                                    from DUNE, WUT,       since May 2016
//
//
// Algorithm for making 2D image-like data from recob::Wire's. Used by CNN codes for training data
// preparation and application of trained models to new data. Also used by PMA to keep images of
// 2D projections used for the track validation.
//

#ifndef DataProviderAlg_h
#define DataProviderAlg_h

// Framework includes
#include "cetlib_except/exception.h"
#include "fhiclcpp/types/Atom.h"
#include "fhiclcpp/types/Comment.h"
#include "fhiclcpp/types/Name.h"
#include "fhiclcpp/types/Sequence.h"
#include "fhiclcpp/types/Table.h"

// LArSoft includes
#include "larcorealg/Geometry/fwd.h"
#include "lardataobj/RecoBase/Wire.h"
#include "larreco/Calorimetry/CalorimetryAlg.h"

#include "CLHEP/Random/JamesRandom.h" // for testing on noise, not used by any reco

// ROOT & C++
#include <iostream>
#include <memory>
#include <optional>
#include <vector>

namespace detinfo {
  class DetectorClocksData;
  class DetectorPropertiesData;
}

namespace img {
  class DataProviderAlg;
  struct DataProviderAlgView {
    unsigned int fNWires;
    unsigned int fNDrifts;
    unsigned int fNScaledDrifts;
    unsigned int fNCachedDrifts;
    std::vector<raw::ChannelID_t> fWireChannels;
    std::vector<std::vector<float>> fWireDriftData;
    std::vector<float> fLifetimeCorrFactors;
  };
}

class img::DataProviderAlg {
public:
  enum EDownscaleMode { kMax = 1, kMaxMean = 2, kMean = 3 };

  struct Config {
    using Name = fhicl::Name;
    using Comment = fhicl::Comment;

    fhicl::Table<calo::CalorimetryAlg::Config> CalorimetryAlg{
      Name("CalorimetryAlg"),
      Comment("Used to eliminate amplitude variation due to electron lifetime.")};

    fhicl::Atom<float> AdcMax{Name("AdcMax"), Comment("Saturation max value")};
    fhicl::Atom<float> AdcMin{Name("AdcMin"), Comment("Saturation min value")};
    fhicl::Atom<float> OutMax{Name("OutMax"), Comment("Output max value")};
    fhicl::Atom<float> OutMin{Name("OutMin"), Comment("Output min value")};

    fhicl::Atom<bool> CalibrateAmpl{Name("CalibrateAmpl"),
                                    Comment("Calibrate ADC values with CalAmpConstants")};

    fhicl::Atom<bool> CalibrateLifetime{Name("CalibrateLifetime"),
                                        Comment("Calibrate ADC values with the electron lifetime")};

    fhicl::Atom<unsigned int> DriftWindow{Name("DriftWindow"),
                                          Comment("Downsampling window (in drift ticks).")};

    fhicl::Atom<std::string> DownscaleFn{Name("DownscaleFn"), Comment("Downsampling function")};

    fhicl::Atom<bool> DownscaleFullView{
      Name("DownscaleFullView"),
      Comment("Downsample full view (faster / lower location precision)")};

    fhicl::Sequence<float> BlurKernel{Name("BlurKernel"), Comment("Blur kernel in wire direction")};

    fhicl::Atom<float> NoiseSigma{Name("NoiseSigma"), Comment("White noise sigma")};

    fhicl::Atom<float> CoherentSigma{Name("CoherentSigma"), Comment("Coherent noise sigma")};
  };

  DataProviderAlg(const fhicl::ParameterSet& pset)
    : DataProviderAlg(fhicl::Table<Config>(pset, {})())
  {}

  DataProviderAlg(const Config& config);

  virtual ~DataProviderAlg();

  bool setWireDriftData(const detinfo::DetectorClocksData& clock_data,
                        const detinfo::DetectorPropertiesData& det_prop,
                        const std::vector<recob::Wire>&
                          wires, // once per plane: setup ADC buffer, collect & downscale ADC's
                        unsigned int plane,
                        unsigned int tpc,
                        unsigned int cryo);

  std::vector<float> const& wireData(size_t widx) const { return fAlgView.fWireDriftData[widx]; }

  std::vector<std::vector<float>> getPatch(size_t wire,
                                           float drift,
                                           size_t patchSizeW,
                                           size_t patchSizeD) const
  {
    bool ok = false;
    std::vector<std::vector<float>> patch;
    if (fDownscaleFullView) {
      ok = patchFromDownsampledView(wire, drift, patchSizeW, patchSizeD, patch);
    }
    else {
      ok = patchFromOriginalView(wire, drift, patchSizeW, patchSizeD, patch);
    }

    if (ok) return patch;
    throw cet::exception("img::DataProviderAlg") << "Patch filling failed." << std::endl;
  }

  float getPixelOrZero(int wire, int drift) const
  {
    size_t didx = getDriftIndex(drift), widx = (size_t)wire;

    if ((widx < fAlgView.fWireDriftData.size()) && (didx < fAlgView.fNCachedDrifts)) {
      return fAlgView.fWireDriftData[widx][didx];
    }
    return 0;
  }

  double getAdcSum() const { return fAdcSumOverThr; }
  size_t getAdcArea() const { return fAdcAreaOverThr; }

  float poolMax(int wire, int drift, size_t r = 0) const;

  //float poolSum(int wire, int drift, size_t r = 0) const;

  unsigned int Cryo() const { return fCryo; }
  unsigned int TPC() const { return fTPC; }
  unsigned int Plane() const { return fPlane; }

  unsigned int NWires() const { return fAlgView.fNWires; }
  unsigned int NScaledDrifts() const { return fAlgView.fNScaledDrifts; }
  unsigned int NCachedDrifts() const { return fAlgView.fNCachedDrifts; }
  unsigned int DriftWindow() const { return fDriftWindow; }

  float ZeroLevel() const { return fAdcZero; }

  double LifetimeCorrection(detinfo::DetectorClocksData const& clock_data,
                            detinfo::DetectorPropertiesData const& det_prop,
                            double tick) const
  {
    return fCalorimetryAlg.LifetimeCorrection(clock_data, det_prop, tick);
  }

protected:
  DataProviderAlgView fAlgView;
  EDownscaleMode fDownscaleMode;
  //std::function<void (std::vector<float> &, std::vector<float> const &, size_t)> fnDownscale;

  size_t fDriftWindow;
  bool fDownscaleFullView;
  float fDriftWindowInv;

  std::vector<float> downscaleMax(std::size_t dst_size,
                                  std::vector<float> const& adc,
                                  size_t tick0) const;
  std::vector<float> downscaleMaxMean(std::size_t dst_size,
                                      std::vector<float> const& adc,
                                      size_t tick0) const;
  std::vector<float> downscaleMean(std::size_t dst_size,
                                   std::vector<float> const& adc,
                                   size_t tick0) const;
  std::vector<float> downscale(std::size_t dst_size,
                               std::vector<float> const& adc,
                               size_t tick0) const
  {
    switch (fDownscaleMode) {
    case img::DataProviderAlg::kMean: return downscaleMean(dst_size, adc, tick0);
    case img::DataProviderAlg::kMaxMean: return downscaleMaxMean(dst_size, adc, tick0);
    case img::DataProviderAlg::kMax: return downscaleMax(dst_size, adc, tick0);
    }
    throw cet::exception("img::DataProviderAlg") << "Downscale mode not supported." << std::endl;
  }

  size_t getDriftIndex(float drift) const
  {
    if (fDownscaleFullView)
      return (size_t)(drift * fDriftWindowInv);
    else
      return (size_t)drift;
  }

  std::optional<std::vector<float>> setWireData(std::vector<float> const& adc,
                                                size_t wireIdx) const;

  bool patchFromDownsampledView(size_t wire,
                                float drift,
                                size_t size_w,
                                size_t size_d,
                                std::vector<std::vector<float>>& patch) const;
  bool patchFromOriginalView(size_t wire,
                             float drift,
                             size_t size_w,
                             size_t size_d,
                             std::vector<std::vector<float>>& patch) const;

  virtual DataProviderAlgView resizeView(detinfo::DetectorClocksData const& clock_data,
                                         detinfo::DetectorPropertiesData const& det_prop,
                                         size_t wires,
                                         size_t drifts);

  // Calorimetry needed to equalize ADC amplitude along drift:
  calo::CalorimetryAlg fCalorimetryAlg;

  // Geometry and detector properties:
  geo::Geometry const* fGeometry;
  geo::WireReadoutGeom const* fWireReadoutGeom;

private:
  float scaleAdcSample(float val) const;
  void scaleAdcSamples(std::vector<float>& values) const;
  std::vector<float> fAmplCalibConst;
  bool fCalibrateAmpl, fCalibrateLifetime;
  unsigned int fCryo = 9999, fTPC = 9999, fPlane = 9999;
  float fAdcMax, fAdcMin, fAdcScale, fAdcOffset, fAdcZero;
  double fAdcSumOverThr, fAdcSumThr;
  size_t fAdcAreaOverThr;

  CLHEP::HepJamesRandom fRndEngine;

  void applyBlur();
  std::vector<float> fBlurKernel; // blur not applied if empty

  void addWhiteNoise();
  float fNoiseSigma; // noise not added if sigma=0

  void addCoherentNoise();
  float fCoherentSigma; // noise not added if sigma=0
};
// ------------------------------------------------------
// ------------------------------------------------------
// ------------------------------------------------------

#endif