Cluster.cxx

Go to the documentation of this file.

// class header
#include "lardataobj/RecoBase/Cluster.h"

// C/C++ standard libraries
#include <algorithm> // std::min(), std::max()
#include <cmath>     // std::sqrt
#include <iomanip>   // std::setw(), ...

namespace {
  template <typename T>
  inline T sqr(T v)
  {
    return v * v;
  }
} // local namespace

namespace recob {

  const Cluster::SentryArgument_t Cluster::Sentry;

  //----------------------------------------------------------------------
  Cluster::Cluster()
    : fNHits(0)
    , fEndWires{0., 0.}
    , fSigmaEndWires{0., 0.}
    , fEndTicks{0., 0.}
    , fSigmaEndTicks{0., 0.}
    , fEndCharges{0., 0.}
    , fAngles{0., 0.}
    , fOpeningAngles{0., 0.}
    , fChargeSum{0., 0.}
    , fChargeStdDev{0., 0.}
    , fChargeAverage{0., 0.}
    , fMultipleHitDensity(0.)
    , fWidth(0.)
    , fID(InvalidID)
    , fView(geo::kUnknown)
    , fPlaneID()
  {} // Cluster::Cluster()

  //----------------------------------------------------------------------
  Cluster::Cluster(float start_wire,
                   float sigma_start_wire,
                   float start_tick,
                   float sigma_start_tick,
                   float start_charge,
                   float start_angle,
                   float start_opening,
                   float end_wire,
                   float sigma_end_wire,
                   float end_tick,
                   float sigma_end_tick,
                   float end_charge,
                   float end_angle,
                   float end_opening,
                   float integral,
                   float integral_stddev,
                   float summedADC,
                   float summedADC_stddev,
                   unsigned int n_hits,
                   float multiple_hit_density,
                   float width,
                   ID_t ID,
                   geo::View_t view,
                   geo::PlaneID const& plane,
                   SentryArgument_t /* sentry = Sentry */
                   )
    : fNHits(n_hits)
    , fEndWires{start_wire, end_wire}
    , fSigmaEndWires{sigma_start_wire, sigma_end_wire}
    , fEndTicks{start_tick, end_tick}
    , fSigmaEndTicks{sigma_start_tick, sigma_end_tick}
    , fEndCharges{start_charge, end_charge}
    , fAngles{start_angle, end_angle}
    , fOpeningAngles{start_opening, end_opening}
    , fChargeSum{integral, summedADC}
    , fChargeStdDev{integral_stddev, summedADC_stddev}
    , fChargeAverage{}
    , fMultipleHitDensity(multiple_hit_density)
    , fWidth(width)
    , fID(ID)
    , fView(view)
    , fPlaneID(plane)
  {

    for (unsigned int mode = cmFirstMode; mode < NChargeModes; ++mode)
      fChargeAverage[mode] = (fNHits > 0) ? fChargeSum[mode] / fNHits : 0.;

  } // Cluster::Cluster(float...)

#if 0
  // FIXME DELME
  //----------------------------------------------------------------------
  //  Addition operator.
  //
  // The two clusters must have the same view and must lay on the same plane.
  // If one of the clusters has an invalid plane, the sum will inherit the
  // other's plane. If both are invalid, sum will also have an invalid plane.
  //
  Cluster Cluster::operator + (Cluster const& other) const {

    // throw exception if the clusters are not from the same view
    if (other.View() != View()) return {};

    if (other.hasPlane() && hasPlane() && (other.Plane() != Plane())) return {};

    const unsigned int n_hits = NHits() + other.NHits();
    double charge_stddev[2];
    for (unsigned int mode = cmFirstMode; mode < NChargeModes; ++mode) {

      // this assumes that the definition of the std dev is unbiased...
      const double this_variance = sqr(ChargeStdDev(mode)) * (NHits()-1.)/NHits();
      const double other_variance
        = sqr(other.ChargeStdDev(mode)) * (other.NHits()-1.) / other.NHits();

      const double e2 = (
        (sqr(NHits()) + sqr(other.NHits())) / sqr(n_hits)
          * sqr(ChargeAverage(mode) - other.ChargeAverage(mode))
        + NHits() * this_variance
        + other.NHits() * other_variance
        ) / (n_hits - 1.); // unbiased

      charge_stddev[mode] = std::sqrt(e2);
    } // for charge mode

    return Cluster (
      std::min(StartWire(), other.StartWire()), // start_wire
      std::min(StartTick(), other.StartTick()), // start_tick
      0.,                                       // start_charge
      0.,                                       // start_angle
      0.,                                       // start_opening
      std::max(EndWire(), other.EndWire()),     // end_wire
      std::max(EndTick(), other.EndTick()),     // end_tick
      0.,                                       // end_charge
      0.,                                       // end_angle
      0.,                                       // end_opening
      Integral() + other.Integral(),            // integral
      charge_stddev[cmFit],                     // integral_stddev
      SummedADC() + other.SummedADC(),          // summedADC
      charge_stddev[cmADC],                     // summedADC_stddev
      n_hits,                                   // n_hits
      0.,                                       // multiple_hit_density
      0.,                                       // width
      InvalidID,                                // ID
      View(),                                   // view
      hasPlane()? Plane(): other.Plane()        // plane
      );

  } // Cluster::operator+ ()

#endif // 0

  //----------------------------------------------------------------------
  // ostream operator.
  //
  std::ostream& operator<<(std::ostream& o, Cluster const& c)
  {
    o << std::setiosflags(std::ios::fixed) << std::setprecision(2);
    o << "Cluster ID " << std::setw(5) << std::right << c.ID() << " : Cryo = " << std::setw(3)
      << std::right << c.Plane().Cryostat << " TPC = " << std::setw(3) << std::right
      << c.Plane().TPC << " Plane = " << std::setw(3) << std::right << c.Plane().Plane
      << " View = " << std::setw(3) << std::right << c.View() << " StartWire = " << std::setw(7)
      << std::right << c.StartWire() << " EndWire = " << std::setw(7) << std::right << c.EndWire()
      << " StartTime = " << std::setw(9) << std::right << c.StartTick()
      << " EndTime = " << std::setw(9) << std::right << c.EndTick()
      << " N hits =      " << std::setw(5) << std::right << c.NHits()
      << " Width =       " << std::setw(5) << std::right << c.Width()
      << " Charge(fit) = " << std::setw(10) << std::right << c.Integral()
      << " Charge(ADC) = " << std::setw(10) << std::right << c.SummedADC();
    return o;
  } // operator<< (ostream, Cluster)

  //----------------------------------------------------------------------
  // < operator.
  //
  bool operator<(Cluster const& a, Cluster const& b)
  {

    if (a.hasPlane() && b.hasPlane() && a.Plane() != b.Plane()) return a.Plane() < b.Plane();
    if (a.View() != b.View()) return a.View() < b.View();
    if (a.ID() != b.ID()) return a.ID() < b.ID();
    if (a.StartTick() != b.StartTick()) return a.StartTick() < b.StartTick();
    if (a.EndTick() != b.EndTick()) return a.EndTick() < b.EndTick();

    return false; // they are equal enough
  }               // operator < (Cluster, Cluster)

} // namespace