TimeOfFlightReco_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// Class:       TimeOfFlightReco
// Plugin Type: producer (art v2_12_01)
// File:        TimeOfFlightReco_module.cc
//
// Generated at Mon Jun  3 11:11:10 2019 by Aidan Medcalf using cetskelgen
// from cetlib version v3_06_00.
//
// Description: ToF reco from BeamlineDigits
////////////////////////////////////////////////////////////////////////

// framework
#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "canvas/Persistency/Common/Ptr.h"
#include "canvas/Utilities/InputTag.h"
#include "fhiclcpp/ParameterSet.h"

// nova
#include "BeamlineRecoBase/BeamlineDigit.h"
#include "BeamlineRecoBase/ToF.h"
#include "BeamlineReco/ToFClusterAlg.h"
#include "BeamlineUtils/BeamlineGeometry.h"
#include "Utilities/AssociationUtil.h"

// stl
#include <iostream>
#include <memory>

namespace beamlinereco {

  class TimeOfFlightReco : public art::EDProducer {
  public:
    explicit TimeOfFlightReco(fhicl::ParameterSet const &p);

    void reconfigure(const fhicl::ParameterSet &p);
    void produce(art::Event &e) override;

  private:

    art::InputTag fBeamlineToFDataLabel;
    double        fHitClusterWindow;
    bool          fUseSiPMForDS;
    std::string   fClusterTimeMethod;

    double calculateTime(ToFPulseCluster c, std::string method) const;
    double calculateTimeOfFlight(double ustime, double dstime) const;
  };

}

// Calculate single time for a cluster
double beamlinereco::TimeOfFlightReco::calculateTime(ToFPulseCluster c, std::string method) const {
  if (c.pulses.size() == 0) {
    return std::numeric_limits<double>::quiet_NaN();
  }

  double t = std::numeric_limits<double>::quiet_NaN();

  if (method == "Mean") {
    t = 0;
    for (auto p : c.pulses) { t += p->StartTimeInNanoSec(); }
    t = t / (double(c.pulses.size()));

  } else if (method == "First") {
    auto d = *std::min_element(c.pulses.begin(), c.pulses.end(),
      [](art::Ptr<brb::BeamlineDigit> p1, art::Ptr<brb::BeamlineDigit> p2){
        return p1->StartTimeInNanoSec() < p2->StartTimeInNanoSec();
      });
    t = d->StartTimeInNanoSec();

  } else {
    throw art::Exception(art::errors::Configuration)
      << "Unknown time of flight cluster time calculation method: " << method << "." << std::endl;
  }

  return t;
}

// Calculate time of flight
double beamlinereco::TimeOfFlightReco::calculateTimeOfFlight(double ustime, double dstime) const {
  return dstime - ustime;
}

beamlinereco::TimeOfFlightReco::TimeOfFlightReco(fhicl::ParameterSet const & p)
{
  produces<std::vector<brb::ToF>>();
  produces<art::Assns<brb::ToF, brb::BeamlineDigit>>();

  reconfigure(p);
}

void beamlinereco::TimeOfFlightReco::reconfigure(const fhicl::ParameterSet & p) {
  fBeamlineToFDataLabel = p.get<art::InputTag>("BeamlineToFDataLabel");
  fHitClusterWindow     = p.get<double>("HitClusterWindow", 10);
  fUseSiPMForDS         = p.get<bool>("UseSiPMForDS", false);
  fClusterTimeMethod    = p.get<std::string>("ClusterTimeMethod", "First");
}

void beamlinereco::TimeOfFlightReco::produce(art::Event & e)
{
  // Producer
  std::unique_ptr<std::vector<brb::ToF>>
    tof_reco(new std::vector<brb::ToF>);
  std::unique_ptr<art::Assns<brb::ToF, brb::BeamlineDigit>>
    tof_assn(new art::Assns<brb::ToF, brb::BeamlineDigit>);

  // Get ToF digits
  art::Handle<std::vector<brb::BeamlineDigit>> digitHandle;
  std::vector<art::Ptr<brb::BeamlineDigit>> digits;
  if (e.getByLabel(fBeamlineToFDataLabel, digitHandle))
    art::fill_ptr_vector(digits, digitHandle);

  // Find clusters
  ToFClusterAlg clusterAlg;
  clusterAlg.FindClusters(digits, fHitClusterWindow);

  std::vector<ToFPulseCluster> usclusters   = clusterAlg.GetClusters(beamlinegeo::ToFCounter::US);
  std::vector<ToFPulseCluster> dsclusters   = clusterAlg.GetClusters(beamlinegeo::ToFCounter::DS);
  std::vector<ToFPulseCluster> sipmclusters = clusterAlg.GetClusters(beamlinegeo::ToFCounter::DSSiPM);

  std::vector<ToFPulseCluster> dscounterclusters = dsclusters;
  if (fUseSiPMForDS) dscounterclusters = sipmclusters;

  // Make tofs
  for (auto usc : usclusters) {
    for (auto dsc : dscounterclusters) {
      double ustime = calculateTime(usc, fClusterTimeMethod);
      double dstime = calculateTime(dsc, fClusterTimeMethod);
      double time = calculateTimeOfFlight(ustime, dstime);
      if (time > 0) {
        brb::ToF tof;
        tof.SetTimestamps(std::make_pair(ustime, dstime));
        tof.SetTime(time);
        tof_reco->push_back(tof);
        util::CreateAssn(*this, e, *(tof_reco.get()), usc.pulses, *(tof_assn.get()), tof_reco->size()-1);
        util::CreateAssn(*this, e, *(tof_reco.get()), dsc.pulses, *(tof_assn.get()), tof_reco->size()-1);
      }
    }
  }

  e.put(std::move(tof_reco));
  e.put(std::move(tof_assn));
}

DEFINE_ART_MODULE(beamlinereco::TimeOfFlightReco)