Cluster.cxx

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// \file    Cluster.cxx
// \brief   Cluster of CellHits. Base of other reconstruction classes.
// \version $Id: Cluster.cxx,v 1.31 2012-10-26 18:35:05 bckhouse Exp $
// \author  Christopher Backhouse - bckhouse@caltech.edu
// \date    $Date: 2012-10-26 18:35:05 $
////////////////////////////////////////////////////////////////////////
#include <cfloat>
#include <iomanip>

#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

#include "RecoBase/Cluster.h"
#include "Calibrator/Calibrator.h"
#include "Geometry/Geometry.h"
#include "GeometryObjects/CellGeo.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/RecoHit.h"
#include "RecoBase/WeightedHit.h"

namespace rb
{
  //------------------------------------------------------------
  Cluster::Cluster()
    : fView(geo::kXorY)
    , fID(0)
    , fNoiseCluster(false)
    , fPrecalcTotalGeV(-1)
  {
  }

  //------------------------------------------------------------
  Cluster::Cluster(const art::PtrVector<CellHit>& cells, int id)
    : fView(geo::kXorY)
    , fID(id)
    , fNoiseCluster(false)
    , fPrecalcTotalGeV(-1)
  {
    for(unsigned int i = 0; i < cells.size(); ++i) Add(cells[i]);
  }

  //------------------------------------------------------------
  Cluster::Cluster(geo::View_t view, int id)
    : fView(view)
    , fID(id)
    , fNoiseCluster(false)
    , fPrecalcTotalGeV(-1)
  {
  }

  //------------------------------------------------------------
  Cluster::Cluster(geo::View_t view, const art::PtrVector<CellHit>& cells, int id)
    : fView(view)
    , fID(id)
    , fNoiseCluster(false)
    , fPrecalcTotalGeV(-1)
  {
    for(unsigned int i = 0; i < cells.size(); ++i) Add(cells[i]);

    assert((fView == geo::kX && fYCell.empty()) ||
     (fView == geo::kY && fXCell.empty()));
  }

  //------------------------------------------------------------
  Cluster::Cluster(geo::View_t view, const std::vector<art::Ptr<CellHit> >& cells, int id)
    : fView(view)
    , fID(id)
    , fNoiseCluster(false)
    , fPrecalcTotalGeV(-1)
  {
    for(unsigned int i = 0; i < cells.size(); ++i) Add(cells[i]);

    assert((fView == geo::kX && fYCell.empty()) ||
     (fView == geo::kY && fXCell.empty()));
  }

  //------------------------------------------------------------
  Cluster::~Cluster()
  {
  }

  //------------------------------------------------------------
  void Cluster::Add(const art::Ptr<rb::CellHit>& cell, double weight)
  {
    fPrecalcTotalGeV = -1; // Invalidate any cached value

    assert(weight >= 0 && weight <= 1);

    // Are we optimizing (omitting) the weights storage?
    const bool optWeights = (weight == 1 && fXWeights.empty() && fYWeights.empty());
    if(!optWeights) EnsureWeightAlloc();

    if(cell->View() == geo::kX){
      assert(fView != geo::kY);
      fXCell.push_back(cell);
      if(!optWeights) fXWeights.push_back(weight);
    }
    else{
      assert(fView != geo::kX);
      fYCell.push_back(cell);
      if(!optWeights) fYWeights.push_back(weight);
    }
  }

  //------------------------------------------------------------
  void Cluster::Add(const art::PtrVector<rb::CellHit>& cells,
                    const std::vector<double>& weights)
  {
    const unsigned int N = cells.size();

    assert(weights.empty() || weights.size() == N);

    for(unsigned int n = 0; n < N; ++n){
      assert(weights.empty() || (weights[n] >= 0 && weights[n] <= 1));
      Add(cells[n], weights.empty() ? 1 : weights[n]);
    }
  }

  //------------------------------------------------------------
  double Cluster::W(const rb::CellHit* hit) const
  {
    geo::View_t otherView = geo::kX;
    if(hit->View() == geo::kX) otherView = geo::kY;

    // TB horizontal modules (Y-view) have readout on opposite end as ND/FD
    art::ServiceHandle<geo::Geometry> geom;
    if(geom->DetId() == novadaq::cnv::kTESTBEAM && hit->View() == geo::kY) return -MeanV(otherView);

    return MeanV(otherView);
  }

  //------------------------------------------------------------
  unsigned int Cluster::NCell(geo::View_t view) const
  {
    switch(view){
    case geo::kX: return NXCell();
    case geo::kY: return NYCell();
    case geo::kXorY: return NCell();
    default: assert(0 && "Unknown view");
    }
  }

  //------------------------------------------------------------
  art::Ptr<rb::CellHit> Cluster::Cell(geo::View_t view,
              unsigned int viewIdx) const
  {
    switch(view){
    case geo::kX: return XCell(viewIdx);
    case geo::kY: return YCell(viewIdx);
    case geo::kXorY: return Cell(viewIdx);
    default: assert(0 && "Unknown view");
    }
  }

  //------------------------------------------------------------
  art::Ptr<rb::CellHit> Cluster::XCell(unsigned int xIdx) const
  {
    assert(xIdx < fXCell.size());

    return fXCell[xIdx];
  }

  //------------------------------------------------------------
  art::Ptr<rb::CellHit> Cluster::YCell(unsigned int yIdx) const
  {
    assert(yIdx < fYCell.size());

    return fYCell[yIdx];
  }

  //------------------------------------------------------------
  art::Ptr<rb::CellHit> Cluster::Cell(unsigned int globalIdx) const
  {
    if(globalIdx < NXCell()) return XCell(globalIdx);
    return YCell(globalIdx-NXCell());
  }

  //------------------------------------------------------------
  art::PtrVector<rb::CellHit> Cluster::AllCells() const
  {
    art::PtrVector<rb::CellHit> all;
    for(size_t c = 0; c < fXCell.size(); ++c) all.push_back(fXCell[c]);
    for(size_t c = 0; c < fYCell.size(); ++c) all.push_back(fYCell[c]);

    return all;
  }

  //------------------------------------------------------------
  std::vector<geo::OfflineChan> Cluster::OfflineChans() const
  {
    std::vector<geo::OfflineChan> ret;
    ret.reserve(NCell());
    for(unsigned int i = 0; i < NCell(); ++i)
      ret.push_back(Cell(i)->OfflineChan());
    return ret;
  }

  //------------------------------------------------------------
  std::vector<rb::WeightedHit> Cluster::WeightedHits() const
  {
    std::vector<rb::WeightedHit> all;
    for(unsigned int cellIdx = 0; cellIdx < NCell(); ++cellIdx)
      all.push_back(rb::WeightedHit(Cell(cellIdx), Weight(cellIdx)));
    return all;
  }

  //------------------------------------------------------------
  double Cluster::Weight(unsigned int globalIdx) const
  {
    if(fXWeights.empty() && fYWeights.empty()) return 1;

    if(globalIdx < NXCell()) return Weight(geo::kX, globalIdx);
    return Weight(geo::kY, globalIdx-NXCell());
  }

  //------------------------------------------------------------
  double Cluster::Weight(geo::View_t view, unsigned int viewIdx) const
  {
    if(fXWeights.empty() && fYWeights.empty()) return 1;

    if(view == geo::kX){
      assert(viewIdx < fXWeights.size());
      return fXWeights[viewIdx];
    }
    else{
      assert(viewIdx < fYWeights.size());
      return fYWeights[viewIdx];
    }
  }

  //------------------------------------------------------------
  rb::Cluster Cluster::Exclude(const rb::Cluster* excl) const
  {
    rb::Cluster ret(fView, fID);
    ret.SetNoise(fNoiseCluster);

    const unsigned int I = NCell();
    const unsigned int J = excl->NCell();
    for(unsigned int i = 0; i < I; ++i){
      const art::Ptr<rb::CellHit> chit = Cell(i);
      bool clash = false;
      for(unsigned int j = 0; j < J; ++j){
        if(excl->Cell(j) == chit){
          clash = true;
          break;
        }
      } // end for j
      if(!clash) ret.Add(chit, Weight(i));
    }

    if(ret.NCell() == 0)
      mf::LogWarning("Cluster") << "Exclude created empty cluster";

    return ret;
  }

  //------------------------------------------------------------
  rb::RecoHit Cluster::RecoHit(const art::Ptr<rb::CellHit>& chit) const
  {
    art::ServiceHandle<calib::Calibrator> cal;
    rb::RecoHit rb(cal->MakeRecoHit(*chit, this->W(chit.get())));
    return rb;
  }

  //------------------------------------------------------------
  rb::RecoHit Cluster::RecoHit(geo::View_t view, unsigned int viewIdx) const
  {
    return RecoHit(Cell(view, viewIdx));
  }

  //------------------------------------------------------------
  rb::RecoHit Cluster::RecoHit(unsigned int globalIdx) const
  {
    return RecoHit(Cell(globalIdx));
  }

  //------------------------------------------------------------
  void Cluster::Clear()
  {
    fPrecalcTotalGeV = -1; // Invalidate any cached value

    fXCell.clear();
    fYCell.clear();
    fXWeights.clear();
    fYWeights.clear();
  }

  //-----------------------------------------------------------
  void Cluster::RemoveHit(const art::Ptr<rb::CellHit> hit)
  {
    fPrecalcTotalGeV = -1; // Invalidate any cached value

    if(hit->View() == geo::kX){
      // try to find the hit in the xcells
      for(unsigned int i = 0;i<NXCell();++i){
        if(fXCell[i] == hit){
          fXCell.erase(fXCell.begin()+i);
          if(!fXWeights.empty()) fXWeights.erase(fXWeights.begin()+i);
          break;
        }
      }
    }
    else if(hit->View() == geo::kY){
     // try to find the hit in the ycells
      for(unsigned int i = 0;i<NYCell();++i){
        if(fYCell[i] == hit){
          fYCell.erase(fYCell.begin()+i);
          if(!fYWeights.empty()) fYWeights.erase(fYWeights.begin()+i);
          break;
        }
      }
    }
  }

  //-----------------------------------------------------------
  void Cluster::RemoveHit(unsigned int globalIdx)
  {
    fPrecalcTotalGeV = -1; // Invalidate any cached value

    assert(globalIdx < NCell());
    if(globalIdx < NXCell()) fXCell.erase(fXCell.begin()+globalIdx);
    else fYCell.erase(fYCell.begin()+globalIdx-NXCell());
  }

  //------------------------------------------------------------
  void Cluster::SetWeight(unsigned int globalIdx, double weight)
  {
    fPrecalcTotalGeV = -1; // Invalidate any cached value

    if(globalIdx < NXCell())
      SetWeight(geo::kX, globalIdx, weight);
    else
      SetWeight(geo::kY, globalIdx-NXCell(), weight);
  }

  //------------------------------------------------------------
  void Cluster::SetWeight(geo::View_t view, unsigned int viewIdx,
                          double weight)
  {
    fPrecalcTotalGeV = -1; // Invalidate any cached value

    assert(weight >= 0 && weight <= 1);

    // Are we optimizing (omitting) the weights storage?
    const bool optWeights = (weight == 1 && fXWeights.empty() && fYWeights.empty());
    if(!optWeights) EnsureWeightAlloc();

    if(view == geo::kX){
      assert(viewIdx < fXCell.size());
      if(!optWeights) fXWeights[viewIdx] = weight;
    }
    else{
      assert(viewIdx < fYCell.size());
      if(!optWeights) fYWeights[viewIdx] = weight;
    }
  }

  //------------------------------------------------------------
  double Cluster::TotalADC() const
  {
    double ret = 0;
    for(unsigned int i = 0; i < NCell(); ++i)
      ret += Cell(i)->ADC()*Weight(i);
    return ret;
  }

  //------------------------------------------------------------
  double Cluster::TotalPE() const
  {
    double ret = 0;
    for(unsigned int i = 0; i < NCell(); ++i)
      ret += Cell(i)->PE()*Weight(i);
    return ret;
  }

  //------------------------------------------------------------
  double Cluster::TotalGeV(EEnergyCalcScheme escheme) const
  {
    if(escheme == kUsePrecalcEnergy){
      if(fPrecalcTotalGeV >= 0) return fPrecalcTotalGeV;

      mf::LogWarning("rb::Cluster::TotalGeV()") << "Passed kUsePrecalcEnergy but no energy was cached in the Cluster. Falling back to recomputing." << std::endl;
    }

    // Optimization that's only safe for literal Clusters
    if(typeid(*this) == typeid(rb::Cluster)) return TotalGeVFastClusterOnly();

    // For derived classes, fall through to the generic implementation

    double ret = 0;
    for(unsigned int i = 0; i < NCell(); ++i){
      const rb::RecoHit rhit = RecoHit(i);
      if(rhit.IsCalibrated()) ret += rhit.GeV()*Weight(i);
    }
    return ret;
  }

  //------------------------------------------------------------
  double Cluster::TotalGeVFastClusterOnly() const
  {
    // This is only safe for objects that return the same W for all hits in the
    // same view. We only know that's the case for Cluster. It's definitely not
    // correct for Prong and Track.
    assert(typeid(*this) == typeid(rb::Cluster));

    art::ServiceHandle<calib::Calibrator> cal;

    double gev = 0;

    // We just use the mean position for the W estimate in all cases.
    // Precalculate once up front, instead of NCell() times as we would
    // otherwise.
    const TVector3 mean = MeanXYZ();

    // X and Y cells have the mean of the opposite view as their W
    for(unsigned int i = 0; i < NXCell(); ++i){
      const rb::RecoHit rhit(cal->MakeRecoHit(*(XCell(i)), mean.Y()));
      if(rhit.IsCalibrated()) gev += rhit.GeV()*Weight(geo::kX, i);
    }
    for(unsigned int i = 0; i < NYCell(); ++i){
      const rb::RecoHit rhit(cal->MakeRecoHit(*(YCell(i)), mean.X()));
      if(rhit.IsCalibrated()) gev += rhit.GeV()*Weight(geo::kY, i);
    }

    return gev;
  }

  //------------------------------------------------------------
  double Cluster::TotalWeight() const
  {
    double ret = 0;
    for(unsigned int cellIdx = 0; cellIdx < NCell(); ++cellIdx)
      ret += Weight(cellIdx);
    return ret;
  }

  //------------------------------------------------------------
  double Cluster::CalorimetricEnergy(EEnergyCalcScheme escheme) const
  {
    art::ServiceHandle<calib::Calibrator> cal;
    return cal->GetGeVToCalorimetricScale()*TotalGeV(escheme);
  }

  //------------------------------------------------------------
  TVector3 Cluster::MinXYZ() const
  {
    TVector3 lo, junk;
    MinMaxMeanXYZ(lo, junk, junk);
    return lo;
  }

  //------------------------------------------------------------
  double Cluster::MinV(geo::View_t view) const
  {
    assert(view == geo::kX || view == geo::kY);
    if(view == geo::kX) return MinXYZ().X();
    return MinXYZ().Y();
  }

  //------------------------------------------------------------
  unsigned int Cluster::MinPlane(geo::View_t view) const
  {
    assert(NCell(view) > 0);
    unsigned int ret = UINT_MAX;
    for(unsigned int i = 0; i < NCell(view); ++i)
      if(Cell(view, i)->Plane() < ret) ret = Cell(view, i)->Plane();
    return ret;
  }

  //------------------------------------------------------------
  unsigned int Cluster::MinCell(geo::View_t view) const
  {
    assert(NCell(view) > 0);
    unsigned int ret = UINT_MAX;
    for(unsigned int i = 0; i < NCell(view); ++i)
      if(Cell(view, i)->Cell() < ret) ret = Cell(view, i)->Cell();
    return ret;
  }

  //------------------------------------------------------------
  double Cluster::MinTNS() const
  {
    assert(NCell() > 0);
    double ret = DBL_MAX;
    for(unsigned int i = 0; i < NCell(); ++i)
      if(Cell(i)->TNS() < ret) ret = Cell(i)->TNS();
    return ret;
  }

  //------------------------------------------------------------
  TVector3 Cluster::MaxXYZ() const
  {
    TVector3 hi, junk;
    MinMaxMeanXYZ(junk, hi, junk);
    return hi;
  }

  //------------------------------------------------------------
  double Cluster::MaxV(geo::View_t view) const
  {
    assert(view == geo::kX || view == geo::kY);
    if(view == geo::kX) return MaxXYZ().X();
    return MaxXYZ().Y();
  }

  //------------------------------------------------------------
  unsigned int Cluster::MaxPlane(geo::View_t view) const
  {
    assert(NCell(view) > 0);
    unsigned int ret = 0;
    for(unsigned int i = 0; i < NCell(view); ++i)
      if(Cell(view, i)->Plane() > ret) ret = Cell(view, i)->Plane();
    return ret;
  }

  //------------------------------------------------------------
  unsigned int Cluster::MaxCell(geo::View_t view) const
  {
    assert(NCell(view) > 0);
    unsigned int ret = 0;
    for(unsigned int i = 0; i < NCell(view); ++i)
      if(Cell(view, i)->Cell() > ret) ret = Cell(view, i)->Cell();
    return ret;
  }

  //------------------------------------------------------------
  double Cluster::MaxTNS() const
  {
    assert(NCell() > 0);
    double ret = -DBL_MAX;
    for(unsigned int i = 0; i < NCell(); ++i)
      if(Cell(i)->TNS() > ret) ret = Cell(i)->TNS();
    return ret;
  }

  //------------------------------------------------------------
  TVector3 Cluster::MeanXYZ(rb::AveragingScheme scheme) const
  {
    TVector3 mean, junk;
    MinMaxMeanXYZ(junk, junk, mean, scheme);
    return mean;
  }

  //------------------------------------------------------------
  double Cluster::MeanV(geo::View_t view, rb::AveragingScheme scheme) const
  {
    assert(view == geo::kX || view == geo::kY);
    if(view == geo::kX) return MeanXYZ(scheme).X();
    return MeanXYZ(scheme).Y();
  }

  //------------------------------------------------------------
  double Cluster::MeanTNS(rb::AveragingScheme scheme) const
  {
    assert(NCell() > 0);
    double ret = 0;
    double totW = 0;
    for(unsigned int i = 0; i < NCell(); ++i){
      double w = Weight(i);
      if(scheme == kByEnergy) w *= Cell(i)->PE();

      ret += Cell(i)->TNS()*w;
      totW += w;
    }
    return ret/totW;
  }

  //------------------------------------------------------------
  unsigned int Cluster::ExtentCell(geo::View_t view) const
  {
    return MaxCell(view)-MinCell(view)+1;
  }

  //------------------------------------------------------------
  bool Cluster::operator<(const Cluster& other) const
  {
    return this->MeanZ() < other.MeanZ();
  }

  //------------------------------------------------------------
  void Cluster::MinMaxMeanXYZ(TVector3& lo, TVector3& hi, TVector3& mean,
                              rb::AveragingScheme scheme) const
  {
    // If this becomes a bottleneck we can cache the results

    /// \todo Won't be quite right in the FD when cells aren't at right angles?

    assert(NCell() > 0);

    art::ServiceHandle<geo::Geometry> geom;

    double Wx = 0;
    double Wy = 0;
    for(int n = 0; n < 3; ++n){
      lo(n) = DBL_MAX;
      hi(n) = -DBL_MAX;
      mean(n) = 0;
    }

    for(unsigned int i = 0; i < NCell(); ++i){
      const art::Ptr<rb::CellHit>& h = Cell(i);
      double xyz[3];
      geom->Plane(h->Plane())->Cell(h->Cell())->GetCenter(xyz);
      const double x = xyz[0];
      const double y = xyz[1];
      const double z = xyz[2];
      const geo::View_t view = geom->Plane(h->Plane())->View();
      double weight = Weight(i);
      if(scheme == kByEnergy) weight *= h->PE();

      if(view == geo::kX){
        if(x < lo.X()) lo.SetX(x);
        if(x > hi.X()) hi.SetX(x);
        mean.SetX(mean.X()+x*weight);
        Wx += weight;
      }
      if(view == geo::kY){
        if(y < lo.Y()) lo.SetY(y);
        if(y > hi.Y()) hi.SetY(y);
        mean.SetY(mean.Y()+y*weight);
        Wy += weight;
      }
      if(z < lo.Z()) lo.SetZ(z);
      if(z > hi.Z()) hi.SetZ(z);
      mean.SetZ(mean.Z()+z*weight);
    }

    if(Wx) mean.SetX(mean.X()/Wx);
    if(Wy) mean.SetY(mean.Y()/Wy);
    if(Wx || Wy) mean.SetZ(mean.Z()/(Wx+Wy));
  }

  //------------------------------------------------------------
  int Cluster::MostContiguousPlanes(geo::View_t view) const
  {
    std::set<int> planes;

    for(unsigned int i = 0; i < NCell(view); ++i)
      planes.insert(Cell(view, i)->Plane());

    if(planes.empty()) return 0;

    // In 3D expect a hit every plane, in 2D only every other plane.
    //
    // This isn't right if the planes aren't alternating in views, but
    // thankfully we stopped building detectors like that.
    const int expGap = (view == geo::kXorY) ? 1 : 2;

    int prevPlane = -999;
    int curRun = 1;
    int longestRun = 0;
    for(auto plane: planes){
      if(plane - prevPlane == expGap){
        ++curRun;
        longestRun = std::max(longestRun, curRun);
      }
      else{
        curRun = 1;
      }
      prevPlane = plane;
    } // end for plane

    return longestRun;
  }

  //------------------------------------------------------------
  int Cluster::MostMissingPlanes(geo::View_t view) const
  {
    std::set<int> planes;

    for(unsigned int i = 0; i < NCell(view); ++i)
      planes.insert(Cell(view, i)->Plane());

    if(planes.empty()) return 0;

    int maxDiff = 1;
    int prevPlane = *planes.begin();
    for(int plane: planes){
      maxDiff = std::max(maxDiff, plane-prevPlane);
      prevPlane = plane;
    }

    // One plane difference as expected in 3D isn't a gap
    if(view == geo::kXorY) return maxDiff-1;

    // Two plane difference in 2D isn't a gap. Beyond that, every two planes
    // difference only indicates one missing plane of hits.
    return (maxDiff-2)/2;
  }

  //------------------------------------------------------------
  int Cluster::NMissingPlanes(geo::View_t view) const
  {
    if(view == geo::kXorY) view = View();

    const int stride = (view == geo::kXorY) ? 1 : 2;
    
    std::set<int> planes;

    for(unsigned int i = 0; i < NCell(view); ++i)
      planes.insert(Cell(view, i)->Plane());

    if(planes.empty()) return 0;
    
    return (*planes.rbegin() - *planes.begin())/stride - (int)planes.size() + 1; 

  }


  //------------------------------------------------------------
  void Cluster::EnsureWeightAlloc()
  {
    if(fXWeights.size() != fXCell.size()) fXWeights.resize(fXCell.size(), 1);
    if(fYWeights.size() != fYCell.size()) fYWeights.resize(fYCell.size(), 1);
  }

  //------------------------------------------------------------
  bool standard_compare(const art::Ptr<rb::CellHit>& a,
                        const art::Ptr<rb::CellHit>& b)
  {
    if (a->fTDC[0]!=b->fTDC[0]) return a->fTDC[0]<b->fTDC[0];
    if (a->Plane()!=b->Plane()) return a->Plane()<b->Plane();
    if (a->Cell() !=b->Cell())  return a->Cell() <b->Cell();
    if (a->fADC[0]!=b->fADC[0]) return a->fADC[0]<b->fADC[0];
    if (a->ID()   !=b->ID())    return a->ID()   <b->ID();
    return false; // Should never get here...
  }

  //------------------------------------------------------------
  void Cluster::StandardSort()
  {
    //
    // Unweighted clusters are an easy sorting job
    //
    if (fXWeights.size()==0 && fYWeights.size()==0) {
      sort(fXCell.begin(), fXCell.end(), standard_compare);
      sort(fYCell.begin(), fYCell.end(), standard_compare);
      return;
    }
    //
    // \todo Sorting weighted hits requires more complicated code
    //
    std::cerr << __FILE__ << ":" << __LINE__ 
              << " Sort for weighted hits not implemented"
              << std::endl;
    abort();
  }

  //------------------------------------------------------------
  void Cluster::SavePrecalcTotalGeV(ESaveGeVMode savemode)
  {
    const double gev = TotalGeV(kRecomputeEnergy);
    if(fPrecalcTotalGeV < 0 || savemode == kResetTotalGeV){
      fPrecalcTotalGeV = gev;
    }
    else{
      if(fabs(gev-fPrecalcTotalGeV) > .01){
        mf::LogError("rb::Cluster") << "SaveTotalGeV() already called. Old value was " << fPrecalcTotalGeV << ", now " << gev
                                    << ". If you intended to reset the stored value, you must pass kResetTotalGeV. Aborting." << std::endl;
        abort();
      }
      else{
        mf::LogWarning("rb::Cluster") << "Called SaveTotalGeV() twice for the same cluster, which is inefficient." << std::endl;
      }
    }
  }

  //------------------------------------------------------------
  std::ostream& operator<< (std::ostream& o, const Cluster& c)
  {
    o << std::setiosflags(std::ios::fixed) << std::setprecision(2);
    o << "View = "    << std::setw(1) << std::right << c.View()
      << " ID = "     << std::setw(1) << std::right << c.ID()
      << " #Cells = " << std::setw(3) << std::right << c.NCell()
      << " Is3D = "   << std::setw(1) << std::right << c.Is3D()
      << " TotalADC = " << std::setw(5) << std::right << c.TotalADC();

    return o;
  }

} // end namespace
////////////////////////////////////////////////////////////////////////