ValidationTests/SdCalibration/Validatrix.cc

Go to the documentation of this file.

#include <fwk/VModule.h>
#include <sdet/Station.h>
#include <evt/Event.h>
#include <sevt/SEvent.h>
#include <sevt/SortCriteria.h>
#include <sevt/Header.h>
#include <sevt/Station.h>
#include <sevt/StationRecData.h>
#include <sevt/PMT.h>
#include <sevt/PMTCalibData.h>
#include <sevt/PMTRecData.h>
#include <utl/ErrorLogger.h>
#include <utl/Math.h>
#include <utl/SVector.h>

#include <cmath>
#include <vector>
#include <fstream>

using namespace evt;
using namespace sevt;
using namespace utl;
using namespace fwk;
using namespace std;


namespace utl {

  template<size_t n, typename T>
  inline
  double
  abs(const SVector<n, T>& v)
  {
    double sum2 = 0;
    for (unsigned int i = 0; i < n; ++i)
      sum2 += Sqr(v[i]);
    return sqrt(sum2);
  }

}


namespace Validatrix {

  template<typename T>
  class AbsoluteDifference {
  public:
    AbsoluteDifference(const double tolerance = 0)
      : fTolerance(tolerance) { }

    static double GetDistance(const T& x, const T& y)
    { return abs(x - y); }

    bool IsCompatible(const T& x, const T& y) const
    { return GetDistance(x, y) <= fTolerance; }

  private:
    double fTolerance;
  };


  template<typename T>
  class RelativeDifference {
  public:
    RelativeDifference(const double tolerance = 0)
      : fTolerance(tolerance) { }

    static
    double
    GetDistance(const T& x, const T& y)
    {
      const double avg = 0.5 * (abs(x) + abs(y));
      const double diff = abs(x - y);
      return avg ? diff/avg : diff;
    }

    bool IsCompatible(const T& x, const T& y) const
    { return GetDistance(x, y) <= fTolerance; }

  private:
    double fTolerance;
  };


  template<typename T, template<typename> class TolerancePolicy = AbsoluteDifference>
  class NamedType {
  public:
    NamedType(const string& name,
              const TolerancePolicy<T>& tolerance = TolerancePolicy<T>(),
              const T& value = T())
      : fName(name), fType(value), fTolerance(tolerance) { }

    operator T&() { return fType; }
    operator const T&() const { return fType; }

    T& operator=(const T& t) { return fType = t; }

    bool
    IsCompatible(const NamedType& nt, const int line = 0)
      const
    {
      if (!fTolerance.IsCompatible((const T&)(fType), (const T&)(nt.fType))) {
        ostringstream err;
        if (line)
          err << "In line " << line << ": ";
        err << fName << ' '
            << (const T&)(fType) << " not compatible with " << (const T&)(nt.fType)
            << ", distance "
            << fTolerance.GetDistance((const T&)(fType), (const T&)(nt.fType));
        ERROR(err);
        return false;
      }
      return true;
    }

  private:
    const string fName;
    T fType;
    const TolerancePolicy<T> fTolerance;
  };


  template<class T, template<typename> class TolerancePolicy = AbsoluteDifference>
  class NamedClass : public T {
  public:
    NamedClass(const string& name,
               const TolerancePolicy<T>& tolerance = TolerancePolicy<T>(),
               const T& value = T())
      : T(value), fName(name), fTolerance(tolerance) { }

    bool
    IsCompatible(const NamedClass& nc, const int line = 0)
      const
    {
      if (!fTolerance.IsCompatible((const T&)(*this), (const T&)(nc))) {
        ostringstream err;
        if (line)
          err << "In line " << line << ": ";
        err << fName << ' '
            << (const T&)(*this) << " not compatible with " << (const T&)(nc)
            << ", distance "
            << fTolerance.GetDistance((const T&)(*this), (const T&)(nc));
        ERROR(err);
        return false;
      }
      return true;
    }

  private:
    const string fName;
    const TolerancePolicy<T> fTolerance;
  };


  const double eps = 1e-5;


  class CalibrationDigest {
  private:
    typedef utl::SVector<3> SVector3;

  public:
    CalibrationDigest() :
      fEventId("EventId"),
      fStationId("StationId"),
      fHighGainSaturation("HighGainSaturation"),
      fLowGainSaturation("LowGainSaturation"),
      fTriggerAlgorithm("TriggerAlgorithm"),
      fVEMCharge("VEMCharge", RelativeDifference<SVector3>(eps), SVector3(0.)),
      fVEMPeak("VEMPeak", RelativeDifference<SVector3>(eps), SVector3(0.)),
      fMuonChargeSlope("MuonChargeSlope", RelativeDifference<SVector3>(eps), SVector3(0.)),
      fMuonPulseDecayTime("MuonPulseDecayTime", RelativeDifference<SVector3>(eps), SVector3(0.)),
      fPeakAmplitude("PeakAmplitude", RelativeDifference<SVector3>(eps), SVector3(0.)),
      fRiseTime("RiseTime", RelativeDifference<SVector3>(eps), SVector3(0.)),
      fFallTime("FallTime", RelativeDifference<SVector3>(eps), SVector3(0.)),
      fT50("T50", RelativeDifference<SVector3>(eps), SVector3(0.)),
      fTotalCharge("TotalCharge", RelativeDifference<SVector3>(eps), SVector3(0.)),
      fAreaOverPeak("AreaOverPeak", RelativeDifference<SVector3>(eps), SVector3(0.)),
      fStationPeakAmplitude("StationPeakAmplitude", RelativeDifference<double>(eps)),
      fSignalStartSlot("SignalStartSlot"),
      fSignalEndSlot("SignalEndSlot"),
      fSignalStartTime("SignalStartTime", AbsoluteDifference<SVector<2, int>>(), SVector<2, int>(0)),
      fTotalSignal("TotalSignal", RelativeDifference<double>(eps)),
      fShapeParameter("ShapeParameter", RelativeDifference<double>(eps)),
      fStationRiseTime("StationRiseTime", RelativeDifference<double>(eps)),
      fStationFallTime("StationFallTime", RelativeDifference<double>(eps)),
      fStationT50("StationT50", RelativeDifference<double>(eps))
    { }

    bool
    IsCompatible(const CalibrationDigest& calib, const int line = 0)
      const
    {
      return fEventId.IsCompatible(calib.fEventId, line) &&
             fStationId.IsCompatible(calib.fStationId, line) &&
             fHighGainSaturation.IsCompatible(calib.fHighGainSaturation, line) &&
             fLowGainSaturation.IsCompatible(calib.fLowGainSaturation, line) &&
             fTriggerAlgorithm.IsCompatible(calib.fTriggerAlgorithm, line) &&
             fVEMCharge.IsCompatible(calib.fVEMCharge, line) &&
             fVEMPeak.IsCompatible(calib.fVEMPeak, line) &&
             fMuonChargeSlope.IsCompatible(calib.fMuonChargeSlope, line) &&
             fMuonPulseDecayTime.IsCompatible(calib.fMuonPulseDecayTime, line) &&
             fPeakAmplitude.IsCompatible(calib.fPeakAmplitude, line) &&
             fRiseTime.IsCompatible(calib.fRiseTime, line) &&
             fFallTime.IsCompatible(calib.fFallTime, line) &&
             fT50.IsCompatible(calib.fT50, line) &&
             fTotalCharge.IsCompatible(calib.fTotalCharge, line) &&
             fAreaOverPeak.IsCompatible(calib.fAreaOverPeak, line) &&
             fStationPeakAmplitude.IsCompatible(calib.fStationPeakAmplitude, line) &&
             fSignalStartSlot.IsCompatible(calib.fSignalStartSlot, line) &&
             fSignalEndSlot.IsCompatible(calib.fSignalEndSlot, line) &&
             fSignalStartTime.IsCompatible(calib.fSignalStartTime, line) &&
             fTotalSignal.IsCompatible(calib.fTotalSignal, line) &&
             fShapeParameter.IsCompatible(calib.fShapeParameter, line) &&
             fStationRiseTime.IsCompatible(calib.fStationRiseTime, line) &&
             fStationFallTime.IsCompatible(calib.fStationFallTime, line) &&
             fStationT50.IsCompatible(calib.fStationT50, line);
    }

    NamedType<int> fEventId;
    NamedType<int> fStationId;
    NamedType<bool> fHighGainSaturation;
    NamedType<bool> fLowGainSaturation;
    NamedType<int> fTriggerAlgorithm;
    NamedClass<SVector3, RelativeDifference> fVEMCharge;
    NamedClass<SVector3, RelativeDifference> fVEMPeak;
    NamedClass<SVector3, RelativeDifference> fMuonChargeSlope;
    NamedClass<SVector3, RelativeDifference> fMuonPulseDecayTime;
    NamedClass<SVector3, RelativeDifference> fPeakAmplitude;
    NamedClass<SVector3, RelativeDifference> fRiseTime;
    NamedClass<SVector3, RelativeDifference> fFallTime;
    NamedClass<SVector3, RelativeDifference> fT50;
    NamedClass<SVector3, RelativeDifference> fTotalCharge;
    NamedClass<SVector3, RelativeDifference> fAreaOverPeak;
    NamedType<double, RelativeDifference> fStationPeakAmplitude;
    NamedType<int> fSignalStartSlot;
    NamedType<int> fSignalEndSlot;
    NamedClass<SVector<2, int>> fSignalStartTime;
    NamedType<double, RelativeDifference> fTotalSignal;
    NamedType<double, RelativeDifference> fShapeParameter;
    NamedType<double, RelativeDifference> fStationRiseTime;
    NamedType<double, RelativeDifference> fStationFallTime;
    NamedType<double, RelativeDifference> fStationT50;
  };


  ostream&
  operator<<(ostream& os, const CalibrationDigest& calib)
  {
    return os << calib.fEventId << ' '
              << calib.fStationId << ' '
              << calib.fHighGainSaturation << ' '
              << calib.fLowGainSaturation << ' '
              << calib.fTriggerAlgorithm << ' '
              << (const SVector<3>&)(calib.fVEMCharge) << ' '
              << (const SVector<3>&)(calib.fVEMPeak) << ' '
              << (const SVector<3>&)(calib.fMuonChargeSlope) << ' '
              << (const SVector<3>&)(calib.fMuonPulseDecayTime) << ' '
              << (const SVector<3>&)(calib.fPeakAmplitude) << ' '
              << (const SVector<3>&)(calib.fRiseTime) << ' '
              << (const SVector<3>&)(calib.fFallTime) << ' '
              << (const SVector<3>&)(calib.fT50) << ' '
              << (const SVector<3>&)(calib.fTotalCharge) << ' '
              << (const SVector<3>&)(calib.fAreaOverPeak) << ' '
              << calib.fStationPeakAmplitude << ' '
              << calib.fSignalStartSlot << ' '
              << calib.fSignalEndSlot << ' '
              << (const SVector<2, int>&)(calib.fSignalStartTime) << ' '
              << calib.fTotalSignal << ' '
              << calib.fShapeParameter << ' '
              << calib.fStationRiseTime << ' '
              << calib.fStationFallTime << ' '
              << calib.fStationT50;
  }


  istream&
  operator>>(istream& is, CalibrationDigest& calib)
  {
    return is >> calib.fEventId
              >> calib.fStationId
              >> calib.fHighGainSaturation
              >> calib.fLowGainSaturation
              >> calib.fTriggerAlgorithm
              >> (SVector<3>&)(calib.fVEMCharge)
              >> (SVector<3>&)(calib.fVEMPeak)
              >> (SVector<3>&)(calib.fMuonChargeSlope)
              >> (SVector<3>&)(calib.fMuonPulseDecayTime)
              >> (SVector<3>&)(calib.fPeakAmplitude)
              >> (SVector<3>&)(calib.fRiseTime)
              >> (SVector<3>&)(calib.fFallTime)
              >> (SVector<3>&)(calib.fT50)
              >> (SVector<3>&)(calib.fTotalCharge)
              >> (SVector<3>&)(calib.fAreaOverPeak)
              >> calib.fStationPeakAmplitude
              >> calib.fSignalStartSlot
              >> calib.fSignalEndSlot
              >> (SVector<2, int>&)(calib.fSignalStartTime)
              >> calib.fTotalSignal
              >> calib.fShapeParameter
              >> calib.fStationRiseTime
              >> calib.fStationFallTime
              >> calib.fStationT50;
  }


  template<typename T>
  inline
  ostream&
  operator<<(ostream& os, const vector<T>& vec)
  {
    for (const auto& v : vec)
      os << v << '\n';
    return os;
  }


  template<typename T>
  istream&
  operator>>(istream& is, vector<T>& vec)
  {
    int lineNo = 0;
    istringstream lineis;
    T t;
    string line;
    while (getline(is, line)) {
      ++lineNo;
      if (line.length() && line[0] != '#') {
        lineis.clear();
        lineis.str(line);
        if (lineis >> t)
          vec.push_back(t);
        else if (!lineis.eof()) {
          ostringstream err;
          err << "Read error in line " << lineNo;
          ERROR(err);
        }
      }
    }
    return is;
  }


  template<typename T>
  inline
  bool
  IsCompatible(const vector<T>& x, const vector<T>& y)
  {
    if (x.size() != y.size()) {
      ERROR("Sizes do not match!");
      return false;
    }

    int line = 0;
    for (typename vector<T>::const_iterator xIt = x.begin(), yIt = y.begin();
         xIt != x.end(); ++xIt, ++yIt)
      if (!xIt->IsCompatible(*yIt, ++line))
        return false;

    return true;
  }


  inline
  ostream&
  operator<<(ostream& os, const vector<int>& vec)
  {
    for (const auto v : vec)
      os << v << '\n';
    return os;
  }


  class Validatrix : public fwk::VModule {
  public:
    virtual ~Validatrix() = default;

    fwk::VModule::ResultFlag Init() override;
    fwk::VModule::ResultFlag Run(evt::Event& event) override;
    fwk::VModule::ResultFlag Finish() override;

  private:
    unsigned int fLineNo = 0;
    bool fStatus = true;
    vector<CalibrationDigest> fReferenceValues;
    vector<CalibrationDigest> fCurrentValues;

    // set this to 'true' to get muon histograms
    static const bool fDumpHistograms = true;
    ofstream fMuonBaseHistoFile;
    ofstream fMuonChargeHistoFile;
    ofstream fMuonPeakHistoFile;
    ofstream fMuonShapeHistoFile;

    REGISTER_MODULE("Validatrix", Validatrix);
  };


  VModule::ResultFlag
  Validatrix::Init()
  {
    fLineNo = 0;
    fStatus = true;
    fReferenceValues.clear();
    fCurrentValues.clear();

    const string referenceFilename = "SdCalibrationValues.ref";
    ifstream reference(referenceFilename);

    if (!reference || !reference.is_open()) {
      ERROR("Failed to open '" + referenceFilename + "' for reading!");
      return eFailure;
    }

    if (!(reference >> fReferenceValues) && reference.eof()) {
      ostringstream info;
      info << "Read " << fReferenceValues.size() << " reference items.";
      INFO(info);
    } else {
      ERROR("Error while reading '" + referenceFilename + "'...");
      return eFailure;
    }

    if (fDumpHistograms) {
      fMuonBaseHistoFile.open("MuonBaseHisto.txt");
      fMuonChargeHistoFile.open("MuonChargeHisto.txt");
      fMuonPeakHistoFile.open("MuonPeakHisto.txt");
      fMuonShapeHistoFile.open("MuonShapeHisto.txt");
    }

    return eSuccess;
  }


  VModule::ResultFlag
  Validatrix::Run(Event& event)
  {
    if (!event.HasSEvent())
      return eSuccess;

    event.GetSEvent().SortStations(ByIncreasingId());
    const auto& sEvent = event.GetSEvent();

    const int firstPMT = sdet::Station::GetFirstPMTId();

    const int eventId = sEvent.GetHeader().GetId();

    for (const auto& s : sEvent.StationsRange()) {
      CalibrationDigest calib;
      calib.fEventId = eventId;
      const int sId = s.GetId();
      calib.fStationId = sId;
      calib.fHighGainSaturation = s.IsHighGainSaturation();
      calib.fLowGainSaturation = s.IsLowGainSaturation();
      for (int p = 0; p < 3; ++p) {
        const int pmtId = p + firstPMT;
        const auto& pmt = s.GetPMT(pmtId);
        if (!pmt.HasCalibData() && pmt.HasRecData()) {
          ostringstream err;
          err << "Station " << sId << ", PMT " << pmtId << " "
                 "has PMTRecData but no PMTCalibData!";
          ERROR(err);
          return eFailure;
        }
        if (fDumpHistograms && pmt.HasCalibData()) {
          // formatted for gnuplot "index first:last" plot flag
          fMuonBaseHistoFile << pmt.GetCalibData().GetMuonBaseHisto() << "\n\n";
          fMuonChargeHistoFile << pmt.GetCalibData().GetMuonChargeHisto() << "\n\n";
          fMuonPeakHistoFile << pmt.GetCalibData().GetMuonPeakHisto() << "\n\n";
          fMuonShapeHistoFile << pmt.GetCalibData().GetMuonShapeHisto() << "\n\n";
        }
        if (!pmt.HasRecData())
          calib.fVEMCharge[p] = calib.fVEMPeak[p] = -1;
        else {
          const PMTRecData& pmtRec = pmt.GetRecData();
          calib.fVEMCharge[p] = pmtRec.GetVEMCharge();
          calib.fVEMPeak[p] = pmtRec.GetVEMPeak();
          calib.fMuonChargeSlope[p] = pmtRec.GetMuonChargeSlope();
          calib.fMuonPulseDecayTime[p] = pmtRec.GetMuonPulseDecayTime();
          calib.fPeakAmplitude[p] = pmtRec.GetPeakAmplitude();
          calib.fRiseTime[p] = pmtRec.GetRiseTime();
          calib.fFallTime[p] = pmtRec.GetFallTime();
          calib.fT50[p] = pmtRec.GetT50();
          calib.fTotalCharge[p] = pmtRec.GetTotalCharge();
          calib.fAreaOverPeak[p] = pmtRec.GetAreaOverPeak();
        }
      }
      if (s.HasRecData()) {
        const StationRecData& sRec = s.GetRecData();
        calib.fStationPeakAmplitude = sRec.GetPeakAmplitude();
        calib.fSignalStartSlot = sRec.GetSignalStartSlot();
        calib.fSignalEndSlot = sRec.GetSignalEndSlot();
        calib.fSignalStartTime[0] = sRec.GetSignalStartTime().GetGPSSecond();
        calib.fSignalStartTime[1] = sRec.GetSignalStartTime().GetGPSNanoSecond();
        calib.fTotalSignal = sRec.GetTotalSignal();
        calib.fShapeParameter = sRec.GetShapeParameter();
        calib.fStationRiseTime = sRec.GetRiseTime();
        calib.fStationFallTime = sRec.GetFallTime();
        calib.fStationT50 = sRec.GetT50();
      }

      if (fStatus) {
        if (fReferenceValues.size() <= fLineNo) {
          ERROR("Running out of reference values!");
          fStatus = false;
        } else
          fStatus = fReferenceValues.at(fLineNo).IsCompatible(calib, fLineNo+1);
      }

      ++fLineNo;
      fCurrentValues.push_back(calib);
    }

    return eSuccess;
  }


  VModule::ResultFlag
  Validatrix::Finish()
  {
    const string currentFilename = "SdCalibrationValues.ref.new";
    ofstream currentFile(currentFilename);

    if (!currentFile || !currentFile.is_open()) {
      ERROR("Cannot open '" + currentFilename + "' for writing!");
      return eFailure;
    }

    if (!(currentFile << fCurrentValues)) {
      ERROR("Cannot write current values to '" + currentFilename + "'!");
      return eFailure;
    }

    if (fReferenceValues.size() != fCurrentValues.size()) {
      ostringstream err;
      err << "Sizes of the current (" << fCurrentValues.size() << ") "
             "and reference (" << fReferenceValues.size() << ") values do not agree!";
      ERROR(err);
      return eFailure;
    }

    if (fDumpHistograms) {
      fMuonBaseHistoFile.close();
      fMuonChargeHistoFile.close();
      fMuonPeakHistoFile.close();
      fMuonShapeHistoFile.close();
    }

    return fStatus ? eSuccess : eFailure;
  }

}