FitterInputHandler.hh

 
#ifndef __RAT_FitterInputHandler__
#define __RAT_FitterInputHandler__
 
#include <RAT/DB.hh>
#include <RAT/DS/EV.hh>
#include <RAT/Log.hh>
#include <algorithm>
#include <string>
 
#include "RAT/DS/DigitPMT.hh"
#include "RAT/DS/FitResult.hh"
#include "RAT/DS/PMT.hh"
 
namespace RAT {
 
class FitterInputHandler {
 public:
  enum class Mode {
    kPMT = 0,
    kDigitPMT = 1,
    kWaveformAnalysis = 2,
  };
 
  Mode mode;
  std::string wfm_ana_name;
  uint64_t hit_cleaning_mask;
  std::string vertex_seed, direction_seed, energy_seed;
 
  FitterInputHandler() : FitterInputHandler(""){};
  FitterInputHandler(const std::string& index) { Configure(index); };
 
  void Configure(const std::string& index) {
    DBLinkPtr tbl = DB::Get()->GetLink("FIT_COMMON", index);
    mode = static_cast<Mode>(tbl->GetI("mode"));
    if (mode == Mode::kWaveformAnalysis) wfm_ana_name = tbl->GetS("waveform_analyzer");
    hit_cleaning_mask = static_cast<uint64_t>(tbl->GetI("hit_cleaning_mask"));
    vertex_seed = tbl->GetS("vertex_seed");
    direction_seed = tbl->GetS("direction_seed");
    energy_seed = tbl->GetS("energy_seed");
  }
 
  DS::FitResult* FindFitResult(const std::string& name) const {
    if (!ev) Log::Die("FitterInputHandler: Trying to acccess event info without registering the event.");
    std::vector<DS::FitResult*> fit_results = ev->GetFitResults();
    for (auto& fit_result : fit_results) {
      if (fit_result->GetFullName() == name) return fit_result;
    }
    return nullptr;
  }
 
  TVector3 GetSeedPosition(const std::string& _fitter_name = "") const {
    std::string fitter_name = _fitter_name;
    if (_fitter_name.empty()) {
      fitter_name = vertex_seed;
    }
    DS::FitResult* fit_result = FindFitResult(fitter_name);
    if (!ValidSeedPosition(fitter_name)) {
      warn << "The Requested Seed Position (from " << fitter_name << ") is not valid." << newline;
      return TVector3(0, 0, 0);
    }
    return fit_result->GetPosition();
  }
 
  bool ValidSeedPosition(const std::string& _fitter_name = "") const {
    std::string fitter_name = _fitter_name;
    if (_fitter_name.empty()) {
      fitter_name = vertex_seed;
    }
    DS::FitResult* fit_result = FindFitResult(fitter_name);
    if (!fit_result) return false;
    return (fit_result->GetEnablePosition() && fit_result->GetValidPosition());
  }
 
  double GetSeedTime(const std::string& _fitter_name = "") const {
    std::string fitter_name = _fitter_name;
    if (_fitter_name.empty()) {
      fitter_name = vertex_seed;
    }
    DS::FitResult* fit_result = FindFitResult(fitter_name);
    if (!ValidSeedTime(fitter_name)) {
      warn << "The Requested Seed Time (from " << fitter_name << ") is not valid." << newline;
      return 0;
    }
    return fit_result->GetTime();
  }
 
  bool ValidSeedTime(const std::string& _fitter_name = "") const {
    std::string fitter_name = _fitter_name;
    if (_fitter_name.empty()) {
      fitter_name = vertex_seed;
    }
    DS::FitResult* fit_result = FindFitResult(fitter_name);
    if (!fit_result) return false;
    return (fit_result->GetEnableTime() && fit_result->GetValidTime());
  }
 
  TVector3 GetSeedDirection(const std::string& _fitter_name = "") const {
    std::string fitter_name = _fitter_name;
    if (_fitter_name.empty()) {
      fitter_name = direction_seed;
    }
    DS::FitResult* fit_result = FindFitResult(fitter_name);
    if (!ValidSeedDirection(fitter_name)) {
      warn << "The Requested Seed Direction (from " << fitter_name << ") is not valid." << newline;
      return TVector3(0, 0, 0);
    }
    return fit_result->GetDirection();
  }
 
  bool ValidSeedDirection(const std::string& _fitter_name = "") const {
    std::string fitter_name = _fitter_name;
    if (_fitter_name.empty()) {
      fitter_name = direction_seed;
    }
    DS::FitResult* fit_result = FindFitResult(fitter_name);
    if (!fit_result) return false;
    return (fit_result->GetEnableDirection() && fit_result->GetValidDirection());
  }
 
  double GetSeedEnergy(const std::string& _fitter_name = "") const {
    std::string fitter_name = _fitter_name;
    if (_fitter_name.empty()) {
      fitter_name = energy_seed;
    }
    DS::FitResult* fit_result = FindFitResult(fitter_name);
    if (!ValidSeedEnergy(fitter_name)) {
      warn << "The Requested Seed Energy (from " << fitter_name << ") is not valid." << newline;
      return 0;
    }
    return fit_result->GetEnergy();
  }
 
  bool ValidSeedEnergy(const std::string& _fitter_name = "") const {
    std::string fitter_name = _fitter_name;
    if (_fitter_name.empty()) {
      fitter_name = energy_seed;
    }
    DS::FitResult* fit_result = FindFitResult(fitter_name);
    if (!fit_result) return false;
    return (fit_result->GetEnableEnergy() && fit_result->GetValidEnergy());
  }
 
  void RegisterEvent(DS::EV* _ev) {
    ev = _ev;
    hitPMTChannels.clear();
    // NOTE: class implementation assumes GetAllPMTIDs and GetAllDigitPMTIDs returns _sorted_ results.
    // This is true since ev->digitpmt and ev->pmt are both std::maps.
    hitPMTChannels = mode == Mode::kPMT ? ev->GetAllPMTIDs() : ev->GetAllDigitPMTIDs();
    // Only use channels that pass the chosen hit cleaning bit
    if (hit_cleaning_mask != 0) {
      std::vector<Int_t> hitPMTChannelsCleaned;
      for (int id : hitPMTChannels) {
        DS::DigitPMT* digitpmt = ev->GetOrCreateDigitPMT(id);
        uint64_t mask = digitpmt->GetHitCleaningMask();
        // passes if none of the bits specified in the cuts are set in the mask
        if ((mask & hit_cleaning_mask) == 0) {
          hitPMTChannelsCleaned.push_back(id);
        }
      }
      hitPMTChannels = hitPMTChannelsCleaned;
    }
  }
 
  const std::vector<Int_t>& GetAllHitPMTIDs() const {
    if (!ev) Log::Die("FitterInputHandler: Trying to acccess event info without registering the event.");
    return hitPMTChannels;
  }
 
  size_t GetNHits() const {
    if (!ev) Log::Die("FitterInputHandler: Trying to acccess event info without registering the event.");
    return hitPMTChannels.size();
  }
 
  double GetCharge(Int_t id) const {
    if (!ev) Log::Die("FitterInputHandler: Trying to acccess event info without registering the event.");
    if (!std::binary_search(hitPMTChannels.begin(), hitPMTChannels.end(), id))
      Log::Die("FitterInputHandler: Trying to access a channel with no hit registered!");
 
    switch (mode) {
      case Mode::kPMT:
        return ev->GetOrCreatePMT(id)->GetCharge();
      case Mode::kDigitPMT:
        return ev->GetOrCreateDigitPMT(id)->GetDigitizedCharge();
      case Mode::kWaveformAnalysis: {
        std::vector<double> charges = GetCharges(id);
        double charge = 0;
        for (size_t ic = 0; ic < charges.size(); ic++) charge += charges[ic];
        return charge;
      }
      default:
        Log::Die("INVALID TYPE! Should never reach here.");
    }
  }
 
  std::vector<double> GetCharges(Int_t id) const {
    if (mode != Mode::kWaveformAnalysis) return std::vector<double>{GetCharge(id)};
    if (!ev) Log::Die("FitterInputHandler: Trying to acccess event info without registering the event.");
    if (!std::binary_search(hitPMTChannels.begin(), hitPMTChannels.end(), id))
      Log::Die("FitterInputHandler: Trying to access a channel with no hit registered!");
    DS::DigitPMT* digitpmt = ev->GetOrCreateDigitPMT(id);
    std::vector<std::string> fitterNames = digitpmt->GetFitterNames();
    if (std::find(fitterNames.begin(), fitterNames.end(), wfm_ana_name) == fitterNames.end()) {
      info << "FitResult not found for pmt id " << id << " " << wfm_ana_name << newline;
    }
    return digitpmt->GetOrCreateWaveformAnalysisResult(wfm_ana_name)->getCharges();
  }
 
  double GetTime(Int_t id) const {
    if (!ev) Log::Die("FitterInputHandler: Trying to acccess event info without registering the event.");
    if (!std::binary_search(hitPMTChannels.begin(), hitPMTChannels.end(), id))
      Log::Die("FitterInputHandler: Trying to access a channel with no hit registered!");
 
    switch (mode) {
      case Mode::kPMT:
        return ev->GetOrCreatePMT(id)->GetTime();
      case Mode::kDigitPMT:
        return ev->GetOrCreateDigitPMT(id)->GetDigitizedTime();
      case Mode::kWaveformAnalysis: {
        DS::DigitPMT* digitpmt = ev->GetOrCreateDigitPMT(id);
        std::vector<std::string> fitterNames = digitpmt->GetFitterNames();
        if (std::find(fitterNames.begin(), fitterNames.end(), wfm_ana_name) == fitterNames.end()) {
          info << "FitResult not found for pmt id " << id << " " << wfm_ana_name << newline;
        }
        return digitpmt->GetOrCreateWaveformAnalysisResult(wfm_ana_name)->getTime(0);
      }
      default:
        Log::Die("INVALID TYPE! Should never reach here.");
    }
  }
 
  std::vector<double> GetTimes(Int_t id) const {
    if (mode != Mode::kWaveformAnalysis) return std::vector<double>{GetTime(id)};
    if (!ev) Log::Die("FitterInputHandler: Trying to acccess event info without registering the event.");
    if (!std::binary_search(hitPMTChannels.begin(), hitPMTChannels.end(), id))
      Log::Die("FitterInputHandler: Trying to access a channel with no hit registered!");
    DS::DigitPMT* digitpmt = ev->GetOrCreateDigitPMT(id);
    std::vector<std::string> fitterNames = digitpmt->GetFitterNames();
    if (std::find(fitterNames.begin(), fitterNames.end(), wfm_ana_name) == fitterNames.end()) {
      info << "FitResult not found for pmt id " << id << " " << wfm_ana_name << newline;
    }
    return digitpmt->GetOrCreateWaveformAnalysisResult(wfm_ana_name)->getTimes();
  }
 
  unsigned int GetNPEs(Int_t id) const {
    if (!ev) Log::Die("FitterInputHandler: Trying to acccess event info without registering the event.");
    if (!std::binary_search(hitPMTChannels.begin(), hitPMTChannels.end(), id))
      Log::Die("FitterInputHandler: Trying to access a channel with no hit registered!");
 
    switch (mode) {
      case Mode::kPMT:
        return 1;  // no nhit information
      case Mode::kDigitPMT:
        return ev->GetOrCreateDigitPMT(id)->GetNCrossings();  // approximate
      case Mode::kWaveformAnalysis: {
        DS::DigitPMT* digitpmt = ev->GetOrCreateDigitPMT(id);
        std::vector<std::string> fitterNames = digitpmt->GetFitterNames();
        if (std::find(fitterNames.begin(), fitterNames.end(), wfm_ana_name) == fitterNames.end()) {
          info << "FitResult not found for pmt id " << id << " " << wfm_ana_name << newline;
        }
        return digitpmt->GetOrCreateWaveformAnalysisResult(wfm_ana_name)->getNPEs();
      }
      default:
        Log::Die("INVALID TYPE! Should never reach here.");
    }
  }
 
 protected:
  DS::EV* ev = nullptr;
  std::vector<Int_t> hitPMTChannels;
};
 
}  // namespace RAT
 
#endif