AltAlgSliceList.cxx

Go to the documentation of this file.

// $Id: AltAlgSliceList.cxx,v 1.25 2007/11/11 08:07:50 rhatcher Exp $
//
// AltAlgSliceList.cxx
//
//    -- Event Slicing Algorithm
//
// Costas Andreopoulos <C.V.Andreopoulos@rl.ac.uk>
// CCLRC, Rutherford Appleton Laboratory
// July 01, 2003

#include <map>
#include <vector>
#include <sstream>
#include <numeric>
#include <cassert>
#include <cmath>
#include <algorithm>
#include <functional>

#include "AltAlgSliceList.h"
#include "AltWrapperStlVecStripHandle.h"

#include "Algorithm/AlgFactory.h"
#include "Algorithm/AlgHandle.h"
#include "Algorithm/AlgConfig.h"
#include "Candidate/CandContext.h"
#include "CandDigit/CandDigitHandle.h"
#include "RecoBase/CandSlice.h"
#include "RecoBase/CandSliceHandle.h"
#include "RecoBase/CandSliceList.h"
#include "RecoBase/CandSliceListHandle.h"
#include "Conventions/Detector.h"
#include "LeakChecker/Lea.h"
#include "MessageService/MsgService.h"
#include "MinosObjectMap/MomNavigator.h"
#include "Navigation/NavKey.h"
#include "Navigation/NavSet.h"
#include "Navigation/XxxItr.h"
#include "RawData/RawDigit.h"
#include "RawData/RawHeader.h"
#include "RawData/RawRecord.h"
#include "RawData/RawChannelId.h"
#include "RawData/RawDigitDataBlock.h"
#include "RecoBase/CandSliceHandle.h"
#include "RecoBase/CandStripHandle.h"
#include "RecoBase/CandStripListHandle.h"
#include "UgliGeometry/UgliGeomHandle.h"
#include "Validity/VldContext.h"

#include <TCanvas.h>
#include <TPad.h>
#include <TLine.h>
#include <TMath.h>
#include <TH1D.h>
#include <TH2D.h>
#include <TH3D.h>

ClassImp(AltAlgSliceList)

//____________________________________________________________________________
CVSID("$Id: AltAlgSliceList.cxx,v 1.25 2007/11/11 08:07:50 rhatcher Exp $");

//____________________________________________________________________________
// --STL Functors (unary/binary predicates, unary/binary functions etc...)
// --Sorting Functions and Functors to be used with the Navigation Tools

struct min_time_distance: 
           public binary_function<CandStripHandle *, CandStripHandle *, bool>
{
   CandStripHandle * _ref;
   min_time_distance(CandStripHandle * ref) : _ref(ref) {}

   bool operator() (CandStripHandle * a, CandStripHandle * b) {
           return  fabs( a->GetCorrBegTime() - _ref->GetCorrBegTime() ) <
                   fabs( b->GetCorrBegTime() - _ref->GetCorrBegTime() );
   }
};
struct is_in_z_window_noref:  public unary_function<CandStripHandle *, bool>
{
   double _pl_min, _pl_max;
   is_in_z_window_noref(double pl_min, double pl_max) :
                                         _pl_min(pl_min), _pl_max(pl_max) {}
   bool operator()(CandStripHandle * csh) {
           return  csh->GetPlane() >= _pl_min && csh->GetPlane() <= _pl_max;
   }
};
struct is_in_t_window_noref:  public unary_function<CandStripHandle *, bool>
{
   double _t_min, _t_max;
   is_in_t_window_noref(double t_min, double t_max) :
                                             _t_min(t_min), _t_max(t_max) {}
   bool operator()(CandStripHandle * csh) {
               return  csh->GetCorrBegTime() >= _t_min && csh->GetCorrBegTime() <= _t_max;
   }
};
struct is_in_tz_window:  public unary_function<CandStripHandle *, bool>
{
   CandStripHandle * _ref;
   double            _dt, _dpln;
   is_in_tz_window(CandStripHandle * ref, double dt, double dpln) :
                                            _ref(ref), _dt(dt), _dpln(dpln) {}
   bool operator()(CandStripHandle * csh) {
     return  fabs( csh->GetCorrBegTime()  - _ref->GetCorrBegTime() )  <= _dt  &&
             abs(  csh->GetPlane() - _ref->GetPlane() ) <= _dpln;
   }
};
struct is_in_tztpos_window:  public unary_function<CandStripHandle *, bool>
{
   CandStripHandle * _ref;
   double            _dt, _dpln, _tpos;
   is_in_tztpos_window(
                 CandStripHandle * ref, double dt, double dpln, double tpos) :
                               _ref(ref), _dt(dt), _dpln(dpln), _tpos(tpos) {}
   bool operator()(CandStripHandle * csh) {
     if( csh->GetPlaneView() != _ref->GetPlaneView() ) return false;
     return  fabs(csh->GetCorrBegTime()  - _ref->GetCorrBegTime())  <= _dt   &&
             abs( csh->GetPlane() - _ref->GetPlane()) <= _dpln &&
             fabs( csh->GetTPos()  - _ref->GetTPos())  <= _tpos;
   }
};
struct min_t:
          public binary_function<CandStripHandle *, CandStripHandle *, bool> {
   bool operator()(CandStripHandle * a, CandStripHandle * b) {
                return a->GetCorrBegTime() < b->GetCorrBegTime();
   }
};
struct min_plane:
          public binary_function<CandStripHandle *, CandStripHandle *, bool> {
   bool operator()(CandStripHandle * a, CandStripHandle * b) {
                return a->GetPlane() < b->GetPlane();
   }
};
struct peak_before : public binary_function<TimeSlice_t, TimeSlice_t, bool>  {
   bool operator()(TimeSlice_t a, TimeSlice_t b) { return a.tpeak < b.tpeak; }
};
struct zero_charge : public unary_function<TimeSlice_t, bool> {
   bool operator()(TimeSlice_t a) { return a.q <= 0; }
};
struct sum_q :     public binary_function<double, CandStripHandle *, double> {
   double operator()(double sum, CandStripHandle * current) {
          return (sum + current->GetCharge());
   }
};
struct sum_qtpos : public binary_function<double, CandStripHandle *, double> {
   double operator()(double sum, CandStripHandle * current) {
          return (sum + (current->GetCharge() * current->GetTPos()));
   }
};
struct sum_qz :    public binary_function<double, CandStripHandle *, double> {
   double operator()(double sum, CandStripHandle * current) {
          return (sum + (current->GetCharge() * current->GetZPos()));
   }
};
struct sum_qtime : public binary_function<double, CandStripHandle *, double> {
   double operator()(double sum, CandStripHandle * current) {
          return (sum + (current->GetCharge() * current->GetCorrBegTime()));
   }
};
struct is_u_view : public unary_function<CandStripHandle *, bool> {
   bool operator()(CandStripHandle * csh) {
          return csh->GetPlaneView() == PlaneView::kU;
   }
};
struct is_v_view : public unary_function<CandStripHandle *, bool> {
   bool operator()(CandStripHandle * csh) {
          return csh->GetPlaneView() == PlaneView::kV;
   }
};
struct same_strip : public unary_function<CandStripHandle *, bool> {
   CandStripHandle * _ref;
   same_strip(CandStripHandle * ref) : _ref (ref) { }
   bool operator()(CandStripHandle * csh) {
          return csh->GetUidInt() == _ref->GetUidInt();
   }
};
class PlaneCandStripHandleKeyFunctor : public CandStripHandleKeyFunctor {

public:
   PlaneCandStripHandleKeyFunctor(int pln, float Q, const char * comp):
                           fPln(pln), fQ(Q), fComp(std::string(comp)) { }
  ~PlaneCandStripHandleKeyFunctor() { }
   virtual NavKey operator() (const CandStripHandle * csh) const {
          if(fComp.compare("le") == 0)
                 return csh->GetPlane() <= fPln && csh->GetCharge() > fQ;
          else if(fComp.compare("g")  == 0)
                 return csh->GetPlane() >  fPln && csh->GetCharge() > fQ;
          else
                 return false;
   }
   virtual CandStripHandleKeyFunctor * Clone() const
                    { return new PlaneCandStripHandleKeyFunctor(*this); }
private:
     int fPln;
     float fQ;
     std::string fComp;
};

static NavKey StripKeyFromTime(const CandStripHandle * csh) {
        return csh->GetCorrBegTime();
}
static NavKey SelectNonZeroQ(const CandStripHandle * csh) {
        return (csh->GetCharge() > 0);
}
//____________________________________________________________________________
AltAlgSliceList::AltAlgSliceList()
{
  LEA_CTOR;

  internalInit(); // initialize private data members

  //-- Msg formatting objects
  fFmtStp  = MsgFormat("%3d");
  fFmtPln  = MsgFormat("%3d");
  fFmtSlc  = MsgFormat("%3d");
  fFmtTime = MsgFormat("%11.10f");
  fFmtQ    = MsgFormat("%6.3f");
}
//____________________________________________________________________________
AltAlgSliceList::~AltAlgSliceList()
{
  LEA_DTOR;
}
//____________________________________________________________________________
void AltAlgSliceList::Trace(const char * /* c */) const
{

}
//____________________________________________________________________________
void AltAlgSliceList::RunAlg(AlgConfig & ac, CandHandle & ch,CandContext & cx)
{
  MSG("AltAlg",Msg::kInfo)
         << "**** Begin of AltAlgSliceList::RunAlg() with parameters" << endl;

  //cout.setf(ios_base::fmtflags(0));
  ac.Print();

  assert(cx.GetDataIn());
  assert(cx.GetDataIn()->InheritsFrom("CandStripListHandle"));

  //-- get a handle to input CandStripList & and create an iterator

  const CandStripListHandle * cslh =
                   dynamic_cast<const CandStripListHandle *> (cx.GetDataIn());

  CandStripHandleItr cshItr   ( cslh->GetDaughterIterator() );
  CandStripHandleItr cshItr_u ( cslh->GetDaughterIterator() );
  CandStripHandleItr cshItr_d ( cslh->GetDaughterIterator() );
  assert( cshItr.GetSet()->Size() > 0 );

  //-- this check will be removed - the algorithm is in principle applicable
  //   to far detector as well

  if(cslh->GetVldContext()->GetDetector() != Detector::kNear) {
    MSG("AltAlg", Msg::kFatal)
             << "DATA COMING FROM A WRONG DETECTOR!! QUITING AT ONCE" << endl;
    return;
  }

  //====================== get algorithm configuration =======================

  getAlgorithmConfiguration(ac);

  assert(fTimeResolution       > 0);
  assert(fPkfNofMergedTimeBins > 0);
  assert(fZDiffBetweenPeaks    > 0);
  assert(fTimeDiffBetweenPeaks > 0);

  //====== Use Navigation tools to suitably sort/select the NavSets ==========
  //======      accessible by my CandStripHandleItr iterators       ==========

  // cshItr            can be used to loop over a time-sorted collection of
  //                   all hit strips with Q>0

  CandStripHandleKeyFunc * key_func_1 =  cshItr.CreateKeyFunc();
  key_func_1->SetFun(SelectNonZeroQ);
  cshItr.GetSet()->AdoptSelectKeyFunc(key_func_1);
  key_func_1=0;
  CandStripHandleKeyFunc * key_func_2 =  cshItr.CreateKeyFunc();
  key_func_2->SetFun(StripKeyFromTime);
  cshItr.GetSet()->AdoptSortKeyFunc(key_func_2);
  key_func_2=0;

  printStripList(cshItr, "Input CandStripList: Q > 0, time-sorted");

  // cshItr_u          Can be used to loop over a time-sorted collection of
  //                   all hit strips with Q>0 in the upstream detector
  //                   This is used to exclude the mu-spectrometer hits in ND
  // cshItr_d          Like cshItr_u but for the downstream detector
  //                   This is used to select only mu-spectrometer hits in ND

  PlaneCandStripHandleKeyFunctor  plane_functor_u(kMuSpec1stPlane-1, 0, "le");
  PlaneCandStripHandleKeyFunctor  plane_functor_d(kMuSpec1stPlane-1, 0, "g");

  CandStripHandleKeyFunc * key_func_u1 =  cshItr_u.CreateKeyFunc();
  key_func_u1->SetFun(plane_functor_u);
  cshItr_u.GetSet()->AdoptSelectKeyFunc(key_func_u1);
  key_func_u1=0;
  CandStripHandleKeyFunc * key_func_u2 =  cshItr_u.CreateKeyFunc();
  key_func_u2->SetFun(StripKeyFromTime);
  cshItr_u.GetSet()->AdoptSortKeyFunc(key_func_u2);
  key_func_u2=0;

  std::ostringstream sstream_u;
  sstream_u  << "Input CandStripList: Pln <= "<< kMuSpec1stPlane-1
             << " && Q > 0, time-sorted";
  printStripList(cshItr_u, sstream_u.str().c_str());

  CandStripHandleKeyFunc * key_func_d1 =  cshItr_d.CreateKeyFunc();
  key_func_d1->SetFun(plane_functor_d);
  cshItr_d.GetSet()->AdoptSelectKeyFunc(key_func_d1);
  key_func_d1=0;
  CandStripHandleKeyFunc * key_func_d2 =  cshItr_d.CreateKeyFunc();
  key_func_d2->SetFun(StripKeyFromTime);
  cshItr_d.GetSet()->AdoptSortKeyFunc(key_func_d2);
  key_func_d2=0;

  std::ostringstream sstream_d;
  sstream_d  << "Input CandStripList: Pln > " << kMuSpec1stPlane-1
             << " && Q > 0, time-sorted";
  printStripList(cshItr_d, sstream_d.str().c_str());


  //=========================  Get spill time-profile  =======================

  //-- get earliest & latest time in the snarl

  getSnarlTimeBoundaries(&cshItr, ftmin, ftmax);

  MSG("AltAlg", Msg::kDebug)
       << "Spill time boundaries: tmin = "
       << fFmtTime( ftmin ) << " - tmax = " << fFmtTime( ftmax ) << endl;

  //-- group raw times in bins given by the near det time resolution

  fNTimeBins   = (int) ( (ftmax-ftmin+15*fTimeResolution) / (fTimeResolution) );

  fTimeProfile = new TH1D("fTimeProfile", "time profile", fNTimeBins,
                             ftmin-5*fTimeResolution, ftmax+10*fTimeResolution);

  MSG("AltAlg", Msg::kDebug)
      << "tprof[min] = " <<  fFmtTime( ftmin -  5*fTimeResolution ) << " "
      << "tprof[max] = " <<  fFmtTime( ftmax + 15*fTimeResolution ) << " "
      << "nbins = " << fNTimeBins << endl;

  // For the time profile use only strips in the upstream detector
  // (the cshItr_u CandStripHandleItr points only to those CandStripHandles -
  //  see above how the Sort/Select KeyFucs are applied using NavigationTools)
  //
  // This is used to exclude hits in the muon spectrometer.
  // These strips can not be included in the subsequent 3-D clustering
  // (in time, tpos, z) space *as they are* since they are not demuxed.
  //
  // Muon spectrometer hits are treated separately and are assigned to
  // event slices *after* they have been formed (see later)

  cshItr_u.ResetFirst();
  while( CandStripHandle * striph = cshItr_u.Ptr() ) {

     (fPkfWeightProfileWithCharge) ?
                 fTimeProfile->Fill( striph->GetCorrBegTime(), striph->GetCharge() )  :
           fTimeProfile->Fill( striph->GetCorrBegTime() )  ;

     cshItr_u.Next();
  }
  fTimeProfile->Rebin(fPkfNofMergedTimeBins);
  fTimeProfileMax = fTimeProfile->GetMaximum();

  //========================= Compute Slice 'Seeds'  =========================

  MSG("AltAlg", Msg::kInfo) << "--> Finding slice seeds " << endl;

  std::vector<TimeSlice_t> slice_seeds = getSliceSeeds(
                                          &cshItr_u, fPkfRecursivePeakSearch);

  MSG("AltAlg", Msg::kInfo)
        << "--> " << slice_seeds.size() << " slice seeds were found " << endl;

  //========================= Fill Slice 'Seed' maps =========================

  MSG("AltAlg", Msg::kInfo)
            << "--> Filling map<int, vector<CandStripHandle *> > event_slices"
            << endl;

  std::map<int, std::vector<CandStripHandle *> >
                          event_slices = fillSliceSeeds(cshItr_u, slice_seeds);

  //-- print / draw the event slices (time structure & UZ, VZ views)
  printSlices(event_slices, "rough slice seeds");
  if(fGrfxDebugGraphics & 1) display(event_slices,
                        "before space/time clustering - rough time splitting");

  //======================= Dissolve very small slices =======================

  //-- handle the most obvious cases where a 'fake' slice-seed was created by
  //   the recursive peak finder -- usually it must be very small or very
  //   scattered slices -- dissolve them and append the orphan strips to the
  //   best available slice...

  if( fRefinementDissolving ) {

     MSG("AltAlg", Msg::kInfo) << "--> Slice filtering" << endl;

     initSliceFiltering(event_slices);

     printSlices(event_slices, "after slice filtering");
     if(fGrfxDebugGraphics & 2) display(event_slices, "after slice filtering");

  } //if !skipping refinement

  //============================= Merge slices ===============================

  //-- handle the most obvious cases where an event is split between two or
  //   more slices...

  if( fRefinementMerging ) {

     MSG("AltAlg", Msg::kInfo)
       << "--> Merging slices where it seems that an event was split" << endl;

     sliceMerger(event_slices);

     printSlices(event_slices, "after slice merging");
     if(fGrfxDebugGraphics & 4) display(event_slices, "after slice merging");
  }

  //============== Re-order "misplaced" strips : 3-D clustering ==============

  //-- compute 3-D clusters in UZT & VZT space with a k-Means
  //   clustering algorithm

  if( fRefinementKMeansClustering ) {

     MSG("AltAlg", Msg::kInfo)
          << "--> strip 2 slice reassignments based on k-Means 3-D clustering"
          << endl;

     eventClustering(event_slices);

     printSlices(event_slices, "after 3-D k-Means clustering");
     if(fGrfxDebugGraphics & 8)
                 display(event_slices, "after k_Means space/time clustering");
  }

  //============== Look for substructure within existing slices ==============

  if( fRefinementMSTClustering ) {

     MSG("AltAlg", Msg::kInfo)
      << "--> Look for big slices that could be (better) splitted in N others"
      << endl;

     //-- need to code this part

     printSlices(event_slices, "after attempt to split slices");
     if(fGrfxDebugGraphics & 16)
                               display(event_slices, "after MST clustering");
  }

  //=========== Assign muon spectrometer hit strips to the slices ==========

  // this is relevant to the near detector only

  if(cslh->GetVldContext()->GetDetector() == Detector::kNear) {

    MSG("AltAlg", Msg::kInfo)
                 << "--> Assign ND Muon Spectrometer hits in slices" << endl;

    //-- hits in the muon spectrometer will be assigned to the closest
    //   (t-wise) slice in which hit strips just before the spectrometer
    //   have been assigned to...

    assignMuonSpectrHits2Slices(cshItr_d, event_slices);
  }

  MSG("AltAlg", Msg::kInfo)
      << "--> " << event_slices.size()
      << ((event_slices.size() == 1) ? " slice was" : " slices were")
      << " found in this spill" << endl;

  printSlices(event_slices, "final reconstructed slices");
  if(fGrfxDebugGraphics & 32) display(event_slices, "final slices");

  //========================= Build CandSlice's ============================

  MSG("AltAlg", Msg::kInfo) << "--> Building CandSlices" << endl;

  buildCandidates(event_slices, ch, cx);

  delete fTimeProfile;

  MSG("AltAlg", Msg::kInfo) << "<-- all CandSlices created!" << endl;
}
//____________________________________________________________________________
void AltAlgSliceList::display(
                  std::map<int, std::vector<CandStripHandle *> > event_slices,
                                                         const char * comment)
{
  MSG("AltAlg", Msg::kInfo)
                 << "--> Debug graphics requested - start drawing..." << endl;

  TCanvas * c = new TCanvas("c", comment, 20, 20, 800, 800);

  eventDisplay(c,event_slices);
  delete c;

  MSG("AltAlg", Msg::kInfo)
          << "<-- Graphics displayed - continuing with slice reco..." << endl;
}
//____________________________________________________________________________
void AltAlgSliceList::eventDisplay(TCanvas * c,
                  std::map<int, std::vector<CandStripHandle *> > event_slices)
{
// Quick debugging 'tool' for tracking down obvious mistakes almost 'online'.
// For thorough testing I run my debugging macros on the std output ntuples

  //-- create the canvas pads & the display histograms

  TPad * padT   = new TPad("padT","strip times",0.01,0.55,0.99,0.99,0);
  TPad * padUZ  = new TPad("padUZ","uz view",0.01,0.01,0.49,0.54,0);
  TPad * padVZ  = new TPad("padVZ","vz view",0.51,0.01,0.99,0.54,0);

  std::vector<TH1D *> time_profileVec;
  std::vector<TH2D *> hVviewVec;
  std::vector<TH2D *> hUviewVec ;

  TH1D * time_profile = 0;
  TH2D * hVview = 0;
  TH2D * hUview = 0;

  const int kNColors   = 16;
  int colors[kNColors] = { 1,40,2,30,3,4,5,6,7,8,9,11,25,42,50,28 };

  int icolor = 0;
  int islice = 0;

  //-- beautify...

  c->SetBorderMode(0);
  c->SetFillColor(0);
  c->Draw();
  padT->Draw();
  padUZ->Draw();
  padVZ->Draw();
  padT->SetBorderMode(0);
  padUZ->SetBorderMode(0);
  padVZ->SetBorderMode(0);

  for(std::map<int, std::vector<CandStripHandle *> >::iterator iter =
             event_slices.begin(); iter != event_slices.end(); ++iter) {

     std::vector<CandStripHandle *> slice_strips = (*iter).second;

     //-- instantiate...

     time_profile = new TH1D("","tb - time stamps",
              fNTimeBins,ftmin-5*fTimeResolution,ftmax+10*fTimeResolution);
     hVview = new TH2D("","VZ-view",200,0,200,100,-4,4);
     hUview = new TH2D("","UZ-view",200,0,200,100,-4,4);

     time_profile->SetMaximum( fTimeProfileMax );

     if(fGrfxTimeProfileLogView) {
        time_profile->SetMinimum( 0.9 );
        padT->SetLogy();
     }

     //-- fill...

     for(std::vector<CandStripHandle *>::iterator strip_iter =
         slice_strips.begin(); strip_iter != slice_strips.end(); ++strip_iter) {

        //time_profile->Fill( (*strip_iter)->GetCorrBegTime() );

        (fPkfWeightProfileWithCharge) ?
           time_profile->Fill( (*strip_iter)->GetCorrBegTime(),
                               (*strip_iter)->GetCharge() )  :
                 time_profile->Fill( (*strip_iter)->GetCorrBegTime() )  ;

        if( (*strip_iter)->GetPlaneView() == PlaneView::kU)
             hUview->Fill((*strip_iter)->GetPlane(), (*strip_iter)->GetTPos());

        if( (*strip_iter)->GetPlaneView() == PlaneView::kV)
             hVview->Fill((*strip_iter)->GetPlane(), (*strip_iter)->GetTPos());
     }

     //-- format & draw...

     padT->cd();
     time_profile->SetStats(kFALSE);
     time_profile->SetFillColor(colors[icolor]);
     time_profile->SetLineColor(colors[icolor]);
     if(islice==0) time_profile->Draw();
     else          time_profile->Draw("SAME");

     padUZ->cd();
     hUview->SetMarkerStyle(20);
     hUview->SetMarkerSize(0.7);
     hUview->SetMarkerColor(colors[icolor]);
     hUview->SetStats(kFALSE);
     if(islice==0) hUview->Draw("P");
     else          hUview->Draw("PSAME");

     padVZ->cd();
     hVview->SetMarkerStyle(20);
     hVview->SetMarkerSize(0.7);
     hVview->SetMarkerColor(colors[icolor]);
     hVview->SetStats(kFALSE);
     if(islice==0) hVview->Draw("P");
     else          hVview->Draw("PSAME");

     c->Update();

     icolor = ((icolor == kNColors - 1) ? ( 0 ) : ( icolor+1 ) );
     islice++;

     time_profileVec.push_back(time_profile);
     hVviewVec.push_back(hVview);
     hUviewVec.push_back(hUview);
  }

  c->SaveAs("slices.gif");

  //-- delete all allocated memory...

  MSG("AltAlg",Msg::kVerbose)
                << "<< Dellacolating memory for graphics objects..." << endl;
  delete padT;
  delete padUZ;
  delete padVZ;

  for(std::vector<TH2D *>::iterator iterTH2 = hVviewVec.begin();
                  iterTH2 != hVviewVec.end(); ++iterTH2)   delete (*iterTH2);
  for(std::vector<TH2D *>::iterator iterTH2 = hUviewVec.begin();
                  iterTH2 != hUviewVec.end(); ++iterTH2)   delete (*iterTH2);
  for(std::vector<TH1D *>::iterator iterTH1 = time_profileVec.begin();
            iterTH1 != time_profileVec.end(); ++iterTH1)   delete (*iterTH1);

  MSG("AltAlg",Msg::kVerbose) << "... Done!>>" << endl;
}
//____________________________________________________________________________
void AltAlgSliceList::eventDisplaySingleSlice(TCanvas * c, int sliceKey,
                  std::map<int, std::vector<CandStripHandle *> > event_slices)
{
  //-- create the canvas pads & the display histograms

  TPad *  padT   = new TPad("padT","strip times",0.01,0.55,0.99,0.99,0);
  TPad *  padUZ  = new TPad("padUZ","uz view",0.01,0.01,0.49,0.54,0);
  TPad *  padVZ  = new TPad("padVZ","vz view",0.51,0.01,0.99,0.54,0);

  TH1D *  time_profileSlc = new TH1D("time_profileSlc","tb - time stamps",
            fNTimeBins,ftmin-5*fTimeResolution,ftmax+10*fTimeResolution);
  TH1D *  time_profileAll = new TH1D("time_profileSlc","tb - time stamps",
            fNTimeBins,ftmin-5*fTimeResolution,ftmax+10*fTimeResolution);

  TH2D *  hVviewSlc = new TH2D("hVviewSlc","VZ-view",200,0,200,100,-4,4);
  TH2D *  hVviewAll = new TH2D("hVviewAll","VZ-view",200,0,200,100,-4,4);
  TH2D *  hUviewSlc = new TH2D("hUviewSlc","UZ-view",200,0,200,100,-4,4);
  TH2D *  hUviewAll = new TH2D("hUviewAll","UZ-view",200,0,200,100,-4,4);

  //-- beautify...

  c->SetBorderMode(0);
  c->SetFillColor(0);
  c->Draw();
  padT->Draw();
  padUZ->Draw();
  padVZ->Draw();
  padT->SetBorderMode(0);
  padUZ->SetBorderMode(0);
  padVZ->SetBorderMode(0);

  //-- fill...

  for(std::map<int, std::vector<CandStripHandle *> >::iterator iter =
            event_slices.begin(); iter != event_slices.end(); ++iter) {

     std::vector<CandStripHandle *> slice_strips = (*iter).second;

     for(std::vector<CandStripHandle *>::iterator strip_iter =
        slice_strips.begin(); strip_iter != slice_strips.end(); ++strip_iter) {

        //time_profileAll->Fill( (*strip_iter)->GetCorrBegTime() );

        (fPkfWeightProfileWithCharge) ?
                 time_profileAll->Fill( (*strip_iter)->GetCorrBegTime(),
                                        (*strip_iter)->GetCharge() )  :
                 time_profileAll->Fill( (*strip_iter)->GetCorrBegTime() )  ;

        if( (*strip_iter)->GetPlaneView() == PlaneView::kU)
            hUviewAll->Fill((*strip_iter)->GetPlane(), (*strip_iter)->GetTPos());

        if( (*strip_iter)->GetPlaneView() == PlaneView::kV)
            hVviewAll->Fill((*strip_iter)->GetPlane(), (*strip_iter)->GetTPos());

        if(sliceKey == (*iter).first) {

            //time_profileSlc->Fill( (*strip_iter)->GetCorrBegTime() );

           (fPkfWeightProfileWithCharge) ?
                 time_profileSlc->Fill( (*strip_iter)->GetCorrBegTime(),
                                        (*strip_iter)->GetCharge()  )  :
                 time_profileSlc->Fill( (*strip_iter)->GetCorrBegTime() )  ;

            if( (*strip_iter)->GetPlaneView() == PlaneView::kU)
                 hUviewSlc->Fill((*strip_iter)->GetPlane(), (*strip_iter)->GetTPos());

            if( (*strip_iter)->GetPlaneView() == PlaneView::kV)
                 hVviewSlc->Fill((*strip_iter)->GetPlane(), (*strip_iter)->GetTPos());
        }
     }

     padT->cd();
     time_profileAll->SetStats(kFALSE);
     time_profileAll->SetFillColor(8);
     time_profileAll->SetLineColor(8);
     time_profileSlc->SetStats(kFALSE);
     time_profileSlc->SetFillColor(2);
     time_profileSlc->SetLineColor(2);
     time_profileAll->Draw();
     time_profileSlc->Draw("SAME");

     padUZ->cd();
     hUviewAll->SetMarkerStyle(20);
     hUviewAll->SetMarkerSize(0.7);
     hUviewAll->SetMarkerColor(8);
     hUviewAll->SetStats(kFALSE);
     hUviewSlc->SetMarkerStyle(20);
     hUviewSlc->SetMarkerSize(0.7);
     hUviewSlc->SetMarkerColor(2);
     hUviewSlc->SetStats(kFALSE);
     hUviewAll->Draw("P");
     hUviewSlc->Draw("PSAME");

     padVZ->cd();
     hVviewAll->SetMarkerStyle(20);
     hVviewAll->SetMarkerSize(0.7);
     hVviewAll->SetMarkerColor(8);
     hVviewAll->SetStats(kFALSE);
     hVviewSlc->SetMarkerStyle(20);
     hVviewSlc->SetMarkerSize(0.7);
     hVviewSlc->SetMarkerColor(2);
     hVviewSlc->SetStats(kFALSE);
     hVviewAll->Draw("P");
     hVviewSlc->Draw("PSAME");

     c->Update();
  }

  delete padT;
  delete padUZ;
  delete padVZ;
  delete time_profileAll;
  delete hVviewAll;
  delete hUviewAll;
  delete time_profileSlc;
  delete hVviewSlc;
  delete hUviewSlc;
}
//____________________________________________________________________________
void AltAlgSliceList::plot3DClusters(TCanvas * c,
                  std::map<int, std::vector<CandStripHandle *> > event_slices)
{
  TPad *  padU  = new TPad("padU","z-tpos-t @ u view",0.01,0.01,0.49,0.99,0);
  TPad *  padV  = new TPad("padV","z-tpos-t @ v view",0.51,0.01,0.99,0.99,0);

  std::vector<TH3D *> uview;
  std::vector<TH3D *> vview;

  TH3D *  hZTposT_u = 0;
  TH3D *  hZTposT_v = 0;

  int colors[26] = {
         1,40,2,30,3,4,5,6,7,8,9,11,25,42,50,28,1,40,2,30,3,4,5,6,7,8};
  int icolor = 0;

  double clight = 300000000;

  //-- beautify canvas

  c->SetBorderMode(0);
  c->SetFillColor(0);
  c->Draw();
  padU->Draw();
  padV->Draw();
  padU->SetBorderMode(0);
  padV->SetBorderMode(0);

  for(std::map<int, std::vector<CandStripHandle *> >::iterator iter =
            event_slices.begin(); iter != event_slices.end(); ++iter) {


     hZTposT_u = new TH3D("hZTposT_u","",200,0,20,100,-4,4,
                                       fNTimeBins,ftmin*clight,ftmax*clight);

     hZTposT_v = new TH3D("hZTposT_v","",200,0,20,100,-4,4,
                                       fNTimeBins,ftmin*clight,ftmax*clight);

     std::vector<CandStripHandle *> slice_strips = (*iter).second;

     for(std::vector<CandStripHandle *>::iterator strip_iter =
          slice_strips.begin(); strip_iter != slice_strips.end(); ++strip_iter) {

        if( (*strip_iter)->GetPlaneView() == PlaneView::kU)
            hZTposT_u->Fill( (*strip_iter)->GetZPos(), (*strip_iter)->GetTPos(),
                              clight*(*strip_iter)->GetCorrBegTime());

        if( (*strip_iter)->GetPlaneView() == PlaneView::kV)
            hZTposT_v->Fill( (*strip_iter)->GetZPos(), (*strip_iter)->GetTPos(),
                             clight*(*strip_iter)->GetCorrBegTime());
     }

     uview.push_back(hZTposT_u);
     vview.push_back(hZTposT_v);

     hZTposT_u->SetMarkerStyle(20);
     hZTposT_u->SetMarkerSize(0.7);
     hZTposT_u->SetMarkerColor(colors[icolor]);
     hZTposT_u->SetStats(kFALSE);

     hZTposT_v->SetMarkerStyle(20);
     hZTposT_v->SetMarkerSize(0.7);
     hZTposT_v->SetMarkerColor(colors[icolor]);
     hZTposT_v->SetStats(kFALSE);

     if(uview.size() == 1) {
         padU->cd();
         hZTposT_u->Draw("P");
         padV->cd();
         hZTposT_v->Draw("P");
     } else {
         padU->cd();
         hZTposT_u->Draw("PSAME");
         padV->cd();
         hZTposT_v->Draw("PSAME");
     }
     c->Update();
     icolor++;
  }


  c->Update();

  for(std::vector<TH3D *>::iterator iterTH3 = uview.begin();
                        iterTH3 != uview.end(); ++iterTH3) delete (*iterTH3);
  for(std::vector<TH3D *>::iterator iterTH3 = vview.begin();
                        iterTH3 != vview.end(); ++iterTH3) delete (*iterTH3);
}
//____________________________________________________________________________
void AltAlgSliceList::getSnarlTimeBoundaries(CandStripHandleItr * cshItr,
                                          double & tmin, double & tmax) const
{
// Get the minimum & the maximum time in the NavSet that the cshItr
// CandStripHandleItr iterates on. The NavSet must be sorted in time (in
// ascending order)
//
  CandStripHandle * striph = 0;

  cshItr->ResetFirst();
  striph = cshItr->Ptr();
  tmin = striph->GetCorrBegTime();

  cshItr->ResetLast();
  striph = cshItr->Ptr();
  tmax = striph->GetCorrBegTime();

  cshItr->ResetFirst();
}
//____________________________________________________________________________
std::vector<TimeSlice_t> AltAlgSliceList::getSliceSeeds(
      CandStripHandleItr * cshItr, bool recursive_mode, PeakFinderConf_t conf)
{
// Recosntruct "slice seeds" ro be refined later.
// The a vector of time boundaries corresponding to the seed slices.
// The vec element i corresponds to the ith seed slice.
// The slice seeds are sorted in pulse height peak time (in ascending order).
//

  MSG("AltAlg", Msg::kDebug)
    << "[-] Peak Finder ------- Nesting level: " << ++fPeakFinderNestingLevel
    << endl;

  assert(fPeakFinderNestingLevel < 50); // crash if things go very wrong...

  std::vector<TimeSlice_t> slice_seeds;

  MSG("AltAlg", Msg::kDebug)  << " |-->[-] Getting gross slice seeds" << endl;

  peakFinder(slice_seeds, conf);
  updateSliceSeedInfo(cshItr, slice_seeds);

  MSG("AltAlg", Msg::kDebug)
                  << "      |--> Seeds found= " << slice_seeds.size() << endl;

  printSliceSeeds(slice_seeds, "initial, very gross slice seeds");

  //-- Check the time-window of these slice-seeds.
  //   In some cases the time-window might be just too large, much bigger than
  //   the typical 'duration' of an event and therefore it collects hit strips
  //   from nearby slices where the overall activity was below the initial
  //   peak-search threshold
  //   In this case just confine the slice seed to an acceptable time-window
  //   around the found pulse height peak, set PeakFinderConf_t to kLowActivity
  //   and look for smaller peaks below the initial threshold.

  //   This step could be by-passed in conjunction with making the peak-finder
  //   very sensitive in the first place (low threshold - small max number of
  //   successive empty bins). However, the two approaches might not be equivalent.
  //   So, I keep both and I need to fine-tune the algorithm externally.

  if(recursive_mode) {

     // by calling the next function, essentially getSliceSeeds() will start
     // calling itself, until the duration  of all currently existing slice seeds
     // has been reduced to the maximum allowed duration and all the 'de-allocated'
     // time has been searched for smaller peaks...

     reduceTimeWindows(cshItr, slice_seeds);
  }

  fPeakFinderNestingLevel--;

  return slice_seeds;
}
//____________________________________________________________________________
void AltAlgSliceList::peakFinder(
          std::vector<TimeSlice_t> & slice_seeds, PeakFinderConf_t conf) const
{
// This method is a gross peak finder that goes through a specified time-profile
// in the time window [tmin, tmax] and identifies 'peaks' (i.e. time intervals
// [t_slice_min, t_slice_max]) where 'maxSuccessiveEmptyBins' were encountered
// (bins with content < 'peakFinderThreshold') after at least one bin with
// non-empty contents was found.
//

  MSG("AltAlg",Msg::kDebug)
     << "      |--> Time-profile peak search with PeakFinderConf_t conf = "
     << asString(conf) << endl;

  int cur_bin, start_bin=0;

  int    n_successive_zeros  = 0;     // zero = ( bin content < threshold )
  int    n_non_zero_values   = 0;     // non-zero = (bin content >= threshold)
  bool   non_zero_bin_found  = false; // <-- in current slice
  bool   slice_is_not_over   = true;
  bool   new_slice_starts    = true;
  Axis_t content_sum         = 0;

  //-- These values depend on PeakFinderConf_t or the algorithm stage
  //   when the peak finder was called

  int     peak_finder_threshold     = 0; // current PeakFinderConf_t value
  int     max_successive_empty_bins = 0; // current PeakFinderConf_t value
  double  tmin                      = 0; // time window: minimum
  double  tmax                      = 0; // time window: maximum
  TH1D *  time_profile              = 0; // t-prof. the peak-finder will run on

  switch(conf) {

  case kMuSpectrometer:
      peak_finder_threshold     = fPkfMuSpecPeakThreshold;
      max_successive_empty_bins = fPkfMuSpecNSuccessiveEmptyBins;
      time_profile              = fSubsetTimeProfile;
      tmin                      = fSubsettmin;
      tmax                      = fSubsettmax;
      break;
  case kLowActivity:
      peak_finder_threshold     = fPkfLowPeakThreshold;
      max_successive_empty_bins = fPkfLowNSuccessiveEmptyBins;
      time_profile              = fSubsetTimeProfile;
      tmin                      = fSubsettmin;
      tmax                      = fSubsettmax;
      break;
  case kDefault:
  default:
      peak_finder_threshold     = fPkfPeakThreshold;
      max_successive_empty_bins = fPkfNSuccessiveEmptyBins;
      time_profile              = fTimeProfile;
      tmin                      = ftmin;
      tmax                      = ftmax;
      break;
  }

  MsgFormat fmtInt("%4d");

  // Compute threshold:
  //     (# of strips / time bin) * (number of merged bins in time profile)
  int threshold = peak_finder_threshold * fPkfNofMergedTimeBins;

  // Find the minimum and maximum bin for this search...
  // Add/subtract a small dt from tmin/tmax so that we do not widen
  // the search range by a bin because of round-off errors if the t value
  // is exactly on the bin-boundary - this would lead to non-exclusive
  // slice-seeds.
  int min_bin  = time_profile->FindBin( (Axis_t) (tmin + fTimeResolution/100.) );
  int max_bin  = time_profile->FindBin( (Axis_t) (tmax - fTimeResolution/100.) );

  MSG("AltAlg",Msg::kDebug)
    << "      |--> For this peak search: "
    << "tmin = " << fFmtTime(tmin) << " bin[" << fmtInt(min_bin) << "] "
    << "tmax = " << fFmtTime(tmax) << " bin[" << fmtInt(max_bin) << "] "
    << endl;

  MSG("AltAlg",Msg::kVerbose)
       << "      |-->[-] Printing time-profile contents & peak finder actions"
       << endl;

  // Loop over the time frames for the time window under examination
  for(cur_bin=min_bin; cur_bin <= max_bin; cur_bin++) {

     // init a new slice seed & mark its starting time bin
     if(new_slice_starts) {

         new_slice_starts = false;
         content_sum      = 0;

         // skip leading zeros...
         for(start_bin = cur_bin; start_bin <= max_bin; start_bin++) {
                  bool skip = time_profile->GetBinContent(start_bin) == 0;

            if(!skip) break;
            else  printTimeBin(start_bin, time_profile, "skipped!");
         }
         cur_bin = start_bin;
     }

     // Keep on adding non-empty (>threshold) time frames to this slice seed

     content_sum += time_profile->GetBinContent(cur_bin);

     if( time_profile->GetBinContent(cur_bin) <= threshold )
         n_successive_zeros++;
     else {
         non_zero_bin_found   = true;
         n_successive_zeros   = 0;
         n_non_zero_values++;
     }

     std::ostringstream sstream;
     sstream  << " thr: "  << fmtInt( threshold )
              << " - N(0's) = " << fmtInt( n_successive_zeros )
              << " - N(1's) = " << fmtInt( n_non_zero_values );

     printTimeBin(cur_bin, time_profile, sstream.str().c_str());

     // If
     // a) we have not reached the maximum allowed number of successive
     //    time frames with pulse height < threshold, or
     // b) we have not included any time frame with pulse height > threshold

     if(n_successive_zeros < max_successive_empty_bins || !non_zero_bin_found)

         slice_is_not_over = true;  // then keep on expanding the slice seed
     else
         slice_is_not_over = false; // else stop expanding the slice seed now!

     // There is no life after tmax:
     //  --the current slice seed is forced to stop growing if we reach
     //    the 'end of the time window'
     //  --this only happens if there was at least one entry in the bins
     //    associated with this slice-seed
     if(cur_bin == max_bin && content_sum > 0) slice_is_not_over = false;

     if(!slice_is_not_over) {

       MSG("AltAlg",Msg::kVerbose)
             << "          **** Peak included - could stop at any time"
             << " - adding any trailing bins below threshold"<< endl;

        int end_bin;
        for(end_bin = cur_bin+1; end_bin <= max_bin; end_bin++) {
          bool expand = time_profile->GetBinContent(end_bin) > 0 &&
                         time_profile->GetBinContent(end_bin) <= threshold;

          if(!expand) break;
          else printTimeBin(end_bin, time_profile,
                                        "trailing bin appended to slice!");
        }

        MSG("AltAlg",Msg::kVerbose)
            << "          **** Stop! One more slice seed is added" << endl;

        cur_bin = end_bin-1;

        //--- reset for next slice seed within this spill
        n_successive_zeros  =  0;
        n_non_zero_values   =  0;
        non_zero_bin_found  =  false;
        slice_is_not_over   =  true;
        new_slice_starts    =  true;

        //--- add time info for the added seed

        TimeSlice_t current_slice_seed;

        current_slice_seed.tmin =
                        (double)   time_profile->GetBinLowEdge(start_bin);
        current_slice_seed.tmax =
                        (double)   time_profile->GetBinLowEdge(end_bin);

        // with the new approach that was causing overlaps...
        //(double) ( time_profile->GetBinLowEdge(end_bin) +
        //           time_profile->GetBinWidth(end_bin)  );

        // this will be update later (and then all_sync = true)
        current_slice_seed.tpeak    = -1;
        current_slice_seed.qt       = -1;
        current_slice_seed.q        = -1;
        current_slice_seed.all_sync = false;
        slice_seeds.push_back(current_slice_seed);
     }

  } // loop over time frames
}
//____________________________________________________________________________
void AltAlgSliceList::updateSliceSeedInfo(CandStripHandleItr * cshItr,
                                 std::vector<TimeSlice_t> & slice_seeds) const
{
// Compute <tpeak> = Sum_{i} (charge(i)*time(i)) / Sum_{i} (charge(i)) for all
// strips within slice seed's [tmin, tmax]
//
  unsigned int islice = 0;
  MsgFormat fmtInt("%3d");

  std::vector<TimeSlice_t>::iterator slice_iter;

  for(slice_iter = slice_seeds.begin();
                            slice_iter != slice_seeds.end(); ++slice_iter) {

      MSG("AltAlg", Msg::kDebug)
             << "    |-->[-] Updating information for slice-seed: "
             << fmtInt(islice++) << endl;

      if(! (*slice_iter).all_sync)
                            updateSingleSliceSeedInfo( cshItr, *slice_iter );
  }
  // Eliminate slice-seeds with 0 charge that might have made it to
  // this point doe to relaxed criteria or rearrangements in previous steps

  removeNullSeeds(slice_seeds);
}
//____________________________________________________________________________
void AltAlgSliceList::updateSingleSliceSeedInfo(
                        CandStripHandleItr * cshItr, TimeSlice_t & seed) const
{
// Compute <tpeak> = Sum_{i} (charge(i)*time(i)) / Sum_{i} (charge(i)) for all
// strips within slice seed's [tmin, tmax]
//
  //cshItr->ResetFirst();

  //-- reset previous selections
  cshItr->GetSet()->Slice();

  //-- select strips in this time-window
  cshItr->GetSet()->Slice( seed.tmin, seed.tmax  );

  MSG("AltAlg", Msg::kDebug)
         << "         |-->[-] t = ["
         << fFmtTime(seed.tmin) << ", " << fFmtTime(seed.tmax ) << "]" << endl;

  MSG("AltAlg", Msg::kDebug)
         << "              |--> selected "
         << cshItr->GetSet()->SizeSelect() << " out of "
         << cshItr->GetSet()->Size()       << " entries" << endl;

  seed.q  = 0.;
  seed.qt = 0.;

  cshItr->ResetFirst();
  while( CandStripHandle * striph = cshItr->Ptr() ) {

      seed.q  += ( striph->GetCharge() );
      seed.qt += ( striph->GetCharge() * striph->GetCorrBegTime()   );

      printStrip(striph);

      cshItr->Next();
  }

  if( seed.q == 0) resetSliceSeed(seed);
  else             seed.tpeak = seed.qt / seed.q;

  seed.all_sync = true;

  printSliceSeed(seed);
}
//____________________________________________________________________________
void AltAlgSliceList::removeNullSeeds(
                                 std::vector<TimeSlice_t> & slice_seeds) const
{
  MSG("AltAlg", Msg::kDebug)
          << "    |-->[-] removing any 'null' slice-seed" << endl;
  MSG("AltAlg", Msg::kDebug)
          << "         |--> init. size = " << slice_seeds.size() << endl;

  std::vector<TimeSlice_t>::iterator slice_iter_new_end =
           remove_if( slice_seeds.begin(), slice_seeds.end(), zero_charge() );
  slice_seeds.erase( slice_iter_new_end, slice_seeds.end() );

  MSG("AltAlg", Msg::kDebug)
          << "         |--> fin. size = " << slice_seeds.size() << endl;
}
//____________________________________________________________________________
void AltAlgSliceList::reduceTimeWindows(CandStripHandleItr * cshItr,
                                       std::vector<TimeSlice_t> & slice_seeds)
{
// The 'recursive' peak finding approach is actually implemented here.
// After the first pass, definite peaks were found and now we are looking
// for smaller peaks in the 'time' between these peaks... by calling the
// method ...that has called this method (and so on, until the full space
// is serached for peaks).
//
  MsgFormat fmtInt("%3d");

  // Maximum time interval that can be assigned to any single slice-seed.
  // bin width = time resolution
  // maximum duration in time resolution units =
  //  # bins before peak + # bins after peak + the bin that contains the peak

  double max_duration = fTimeResolution *
                          (fPkfTimeWindowBefPeak + fPkfTimeWindowAftPeak + 1);

  MSG("AltAlg", Msg::kDebug)
      << "Narrowing seeds: duration limited to "
      << fFmtTime(max_duration)
      << " (-" << fPkfTimeWindowBefPeak << ", "
      << "+"   << fPkfTimeWindowAftPeak << ") bins around peak's bin"
      << endl;

  // caution: in this loop I perform insert() / erase() operations on the
  //          "std::vector<TimeSlice_t> slice_seeds"...
  //          Do not use STL iterators in this case because the condition
  //          "iter != slice_seeds.end()" gets invalidated by these operations.

  for(unsigned int islice = 0; islice < slice_seeds.size(); islice++) {

      double dt = slice_seeds[islice].tmax - slice_seeds[islice].tmin;

      if(dt > max_duration) {

          bool search_before_peak = false;
          bool search_after_peak  = false;

          MSG("AltAlg", Msg::kDebug)
             << "[-] Got a wide slice seed: Dt = " << fFmtTime( dt )
             << " > (maxDuration = ) " << fFmtTime(max_duration)
             << endl;

          //-- reduce the slice-seed width

          MSG("AltAlg", Msg::kDebug)
              << " |-->[-] initial:  t(slc-seed="    << fmtInt( islice ) << ")"
              << " = ["  << fFmtTime( slice_seeds[islice].tmin )
              << ", "    << fFmtTime( slice_seeds[islice].tmax )
              << "] with pulse height peak at t = "
              << fFmtTime( slice_seeds[islice].tpeak )
              << endl;

          int peak_bin = fTimeProfile->FindBin(
                                 (Axis_t) slice_seeds[islice].tpeak );

          int min_bin  = peak_bin - fPkfTimeWindowBefPeak;
          int max_bin  = peak_bin + fPkfTimeWindowAftPeak;

          double tmin_new = (double)   fTimeProfile->GetBinLowEdge(min_bin);
          double tmax_new = (double) ( fTimeProfile->GetBinLowEdge(max_bin) +
                                       fTimeProfile->GetBinWidth(max_bin)  );

          double tmin_old = slice_seeds[islice].tmin;
          double tmax_old = slice_seeds[islice].tmax;

          //-- check whether the 'after peak' time boundary should be shifted

          if( tmax_old > tmax_new ) {
             search_after_peak            = true;
             slice_seeds[islice].tmax     = tmax_new;
             slice_seeds[islice].all_sync = false;
          }

          //-- check whether the 'before peak' time boundary should be shifted

          if( tmin_old < tmin_new ) {
             search_before_peak           = true;
             slice_seeds[islice].tmin     = tmin_new;
             slice_seeds[islice].all_sync = false;
          }

          MSG("AltAlg", Msg::kDebug)
              << " |-->[-] updated:  t(slc-seed="    << fmtInt( islice ) << ")"
              << " = ["  << fFmtTime( slice_seeds[islice].tmin )
              << ", "    << fFmtTime( slice_seeds[islice].tmax )
              << "] with pulse height peak at t = "
              << fFmtTime( slice_seeds[islice].tpeak ) << endl;


          //-- tmin/tmax changed --keep all TimeSlice_t info in sync
          MSG("AltAlg", Msg::kDebug)
              << "      |-[-] sync TimeSlice info" << endl;

          if(!slice_seeds[islice].all_sync)
                      updateSingleSliceSeedInfo(cshItr, slice_seeds[islice]);

          //-- change PeakFinderConf_t to kLowActivity and look for new
          //   slice seeds in 'de-allocated' time-space

          if( search_after_peak ) {

            MSG("AltAlg", Msg::kDebug)
               << "      |-->[-] Looking for small peaks (after the main peak) "
               << "in [" << fFmtTime(tmax_new)
               << ", "   << fFmtTime(tmax_old) << "] time window" << endl;

            fSubsettmin = tmax_new;
            fSubsettmax = tmax_old;

            std::vector<TimeSlice_t> new_seeds = findSmallerPeaks(cshItr);
            if(new_seeds.size() > 0) mergeSeeds(new_seeds, slice_seeds);
          }

          if( search_before_peak ) {

            MSG("AltAlg", Msg::kDebug)
              << "      |-->[-] Looking for small peaks (before the main peak) "
              << "in [" << fFmtTime(tmin_old)
              << ", "   << fFmtTime(tmin_new) << "] time window" << endl;

            fSubsettmin = tmin_old;
            fSubsettmax = tmin_new;

            std::vector<TimeSlice_t> new_seeds = findSmallerPeaks(cshItr);
            if(new_seeds.size() > 0) mergeSeeds(new_seeds, slice_seeds);
          }

      } //  if dt > max_duration

      //islice++;
  } // slice seed iter

}
//____________________________________________________________________________
std::vector<TimeSlice_t> AltAlgSliceList::findSmallerPeaks(
                                                  CandStripHandleItr * cshItr)
{
// Runs the peak-finder (in a kLowActivity configuration state) at a time
// window within spill (the time-window de-allocated when the duration of all
// major peaks was limited) so as to search for smaller peaks.
// The time window must have been defined prior to calling this method
// (by setting the fSubsettmin, fSubsettmax private data members).
//
  std::vector<TimeSlice_t> new_seeds;

  //-- select strips in this time-window

  double tmin = fSubsettmin; // copy here because these two data members change
  double tmax = fSubsettmax; // in nested calls of the peak-finder

  MSG("AltAlg", Msg::kDebug)
     << "           |--> cshItr.GetSet()->Slice("
     << fFmtTime(tmin) << ", " << fFmtTime(tmax) << ")" << endl;

  cshItr->GetSet()->Slice();       // reset previous
  cshItr->GetSet()->Slice(tmin, tmax);

  MSG("AltAlg", Msg::kDebug)
     << "           |--> selected "
     << cshItr->GetSet()->SizeSelect() << " out of " << cshItr->GetSet()->Size()
     << " entries" << endl;

  if( cshItr->GetSet()->SizeSelect() > 0) {

       // watch out for memory leaks when running the peak finder recursively
       if(fSubsetTimeProfile) delete fSubsetTimeProfile;

       fSubsetTimeProfile = createSubsetTimeProfile(cshItr);

       MSG("AltAlg", Msg::kDebug)
            << "---------------------------- PEAK FINDER (NESTED LEVEL)/ "
            << "---------------------------- " << endl;

       //-- search for new slice-seeds
       new_seeds = getSliceSeeds(cshItr, fPkfRecursivePeakSearch, kLowActivity);

       MSG("AltAlg", Msg::kDebug)
            << "------------------------------------------- /PEAK FINDER "
            << "------------------------------------------- " << endl;

       MSG("AltAlg", Msg::kDebug) << "           |--> search in the ["
          << fFmtTime(tmin) << ", " << fFmtTime(tmax)
          << "] time window gave " << new_seeds.size() << " more seeds" << endl;

       if(fSubsetTimeProfile) delete fSubsetTimeProfile;
       fSubsetTimeProfile = 0;

  } else

       MSG("AltAlg", Msg::kDebug)
            << "           |--> cshItr.GetSet()->SizeSelect() = 0, "
            << "terminating recursive calls for this window" << endl;

  return new_seeds;
}
//____________________________________________________________________________
void AltAlgSliceList::mergeSeeds(std::vector<TimeSlice_t> new_seeds,
                              std::vector<TimeSlice_t> & full_seed_list) const
{
  std::vector<TimeSlice_t>::iterator time_slice_iter;

   //-- add new slices
   MSG("AltAlg", Msg::kVerbose)
            << "           |--> adding new seeds - initial number: "
            << full_seed_list.size() << endl;

   for(time_slice_iter = new_seeds.begin();
              time_slice_iter != new_seeds.end(); ++time_slice_iter) {

              full_seed_list.push_back( *time_slice_iter );
   }
   MSG("AltAlg", Msg::kVerbose)
            << "           |--> seeds added - new number: "
            << full_seed_list.size() << endl;

  //printSliceSeeds(new_seeds, "new seeds that were added in this step");
  //printSliceSeeds(full_seed_list, "full list of current slice seeds");
}
//____________________________________________________________________________
std::map<int,std::vector<CandStripHandle *> > AltAlgSliceList::fillSliceSeeds(
             CandStripHandleItr cshItr, std::vector<TimeSlice_t> & slice_seeds)
{
// This method is called when the construction of the slice seeds has finished
// and fills the std::map<int,std::vector<CandStripHandle *> > structure that
// will be used for holding the reconstructed slices in the subsequent steps
// of the event slicing algorithm.
// This change is necessary because it will no longer be adequate to describe
// the slices with mutually exclusive [tmin, tmax] time-intervals.
// During slice-refinement, the reconstructed slices might pick up strips from
// other slices so as to minimize a combined temporal-spatial distance and
// therefore *hard* time boundaries might not be respected.
//
  MSG("AltAlg", Msg::kDebug) <<  "* Constructing rough slices map"  << endl;

  // 'narrowing' might completely strip a scattered seed
  removeNullSeeds(slice_seeds);

  sort( slice_seeds.begin(), slice_seeds.end(), peak_before() );

  printSliceSeeds(slice_seeds, "final slice seeds -- sorted in peak time");

  MSG("AltAlg", Msg::kDebug)
         << " |--> selected " << cshItr.GetSet()->SizeSelect()
         << " out of "        << cshItr.GetSet()->Size()
         << " entries" << endl;

  int islice = 0;
  MsgFormat fmtInt("%3d");

  std::map<int, std::vector<CandStripHandle *> >  event_slices;

  std::vector<TimeSlice_t>::const_iterator slice_seeds_iter;

  for(slice_seeds_iter = slice_seeds.begin();
                slice_seeds_iter != slice_seeds.end(); ++slice_seeds_iter ) {

     MSG("AltAlg", Msg::kVerbose)
         << "Assigning strips to slice seed: " << fmtInt(islice) << endl;;
     MSG("AltAlg", Msg::kVerbose)
         << "[-] cshItr.GetSet()->Slice("
         <<   "t = " << fFmtTime( (*slice_seeds_iter).tmin )
         << ", t = " << fFmtTime( (*slice_seeds_iter).tmax ) << ")"
         << " - dt<slice> = "
         << fFmtTime( (*slice_seeds_iter).tmax - (*slice_seeds_iter).tmin )
         << endl;

     //-- reset previous selections
     cshItr.GetSet()->Slice();

     //-- select strips in this time-window
     cshItr.GetSet()->Slice( (*slice_seeds_iter).tmin,
                             (*slice_seeds_iter).tmax  );

     MSG("AltAlg", Msg::kVerbose)
        << " |--> selected " << cshItr.GetSet()->SizeSelect() << " out of "
        << cshItr.GetSet()->Size() << " entries" << endl;

     std::vector<CandStripHandle *> slice_strips;

     MSG("AltAlg", Msg::kVerbose)
         << "[-] while( CandStripHandle * striph = cshItr.Ptr() )" << endl;

     cshItr.ResetFirst();
     while( CandStripHandle * striph = cshItr.Ptr() ) {
        slice_strips.push_back(striph);

        MSG("AltAlg", Msg::kVerbose)
           << " |--> slice_strips.push_back( "
           << "pl = "     << fFmtPln(  striph->GetPlane()   )
           << ", str = "  << fFmtStp(  striph->GetStrip()   )
           << ", Q = "    << fFmtQ(    striph->GetCharge()  )
           << ", view = " << PlaneView::AsString( striph->GetPlaneView() )
           << ", t = "    << fFmtTime( striph->GetCorrBegTime() )
           << " )" << endl;

        cshItr.Next();
     }

     MSG("AltAlg", Msg::kVerbose)
         << "[-] add pair: event_slices.insert( "
         << "std::map<int, std::vector<CandStripHandle *> >::"
         << "value_type(" << islice << ", " << &slice_strips << ") );"
         << endl;

     event_slices.insert(
       std::map<int, std::vector<CandStripHandle *> >::value_type(
                                                      islice, slice_strips) );
     islice++;
  }

  MSG("AltAlg", Msg::kDebug)
               << "* Constructing of rough reco slices map completed" << endl;

  return event_slices;
}
//____________________________________________________________________________
void AltAlgSliceList::sortSlicesInTime(
          std::map<int, std::vector<CandStripHandle *> > & event_slices) const
{
// Sorts CandStripHandles associated with each slice is ascending time order
//
  MSG("AltAlg", Msg::kDebug) << "* Sorting Slices in time " << endl;

  std::map<int, std::vector<CandStripHandle *> >::iterator slice_iter;

  for(slice_iter = event_slices.begin();
                              slice_iter != event_slices.end(); ++slice_iter)
     sort( (*slice_iter).second.begin(), (*slice_iter).second.end(), min_t());
}
//____________________________________________________________________________
void AltAlgSliceList::sortSlicesInPlaneNo(
          std::map<int, std::vector<CandStripHandle *> > & event_slices) const
{
// Sorts CandStripHandles associated with each slice is ascending plane number
//
  MSG("AltAlg", Msg::kDebug) << "* Sorting Slices in plane number " << endl;

  std::map<int, std::vector<CandStripHandle *> >::iterator slice_iter;

  for(slice_iter = event_slices.begin();
                               slice_iter != event_slices.end(); ++slice_iter)
          sort( (*slice_iter).second.begin(),
                                     (*slice_iter).second.end(), min_plane());
}
//____________________________________________________________________________
void AltAlgSliceList::printSliceSeeds(
     const std::vector<TimeSlice_t> & slice_seeds, const char * comment) const
{
// Method for printing the slice seeds found by the time-profile peak-finder
//
  MSG("AltAlg", Msg::kDebug) << "[-] Printing: " << comment << endl;

  int islice=0;
  MsgFormat fmtInt("%3d");

  std::vector<TimeSlice_t>::const_iterator time_slice_iter;

  for(time_slice_iter = slice_seeds.begin();
               time_slice_iter != slice_seeds.end(); ++time_slice_iter) {

      if( (*time_slice_iter).tpeak == -1 )
        MSG("AltAlg", Msg::kDebug)
             << " |--> t(slc-seed ="    << fmtInt( islice++ ) << ")"
             << " = ["  << fFmtTime( (*time_slice_iter).tmin )
             << ", "    << fFmtTime( (*time_slice_iter).tmax )
             << "] with pulse height peak not defined yet" << endl;
      else
         MSG("AltAlg", Msg::kDebug)
             << " |--> t(slc-seed ="    << fmtInt( islice++ ) << ")"
             << " = ["  << fFmtTime( (*time_slice_iter).tmin )
             << ", "    << fFmtTime( (*time_slice_iter).tmax )
             << "] with pulse height peak at t = "
             << fFmtTime( (*time_slice_iter).tpeak ) << endl;
  }
}
//____________________________________________________________________________
void AltAlgSliceList::printSliceSeed(const TimeSlice_t & seed) const
{
   MSG("AltAlg", Msg::kDebug)
          << "              |--> Seed info: t = ["
          << fFmtTime(seed.tmin) << ", "    << fFmtTime(seed.tmax)
          << "] with peak at t = " << fFmtTime(seed.tpeak)
          << " and Qtot = " << fFmtQ(seed.q) << endl;
}
//____________________________________________________________________________
void AltAlgSliceList::resetSliceSeed(TimeSlice_t & seed) const
{
  seed.tmin  = 0;  seed.tpeak = 0;  seed.tmax  = 0;
  seed.q     = 0;  seed.qt    = 0;  seed.all_sync = true;
}
//____________________________________________________________________________
void AltAlgSliceList::printSlices(
          const std::map<int, std::vector<CandStripHandle *> > & event_slices,
                                                   const char * comment) const
{
// Method for printing the reconstructed event slices
//
  MSG("AltAlg", Msg::kDebug) << "* Printing Reco Slices: " << comment << endl;

  std::map<int, std::vector<CandStripHandle *> >::const_iterator slice_iter;

  for(slice_iter = event_slices.begin();
               slice_iter != event_slices.end(); ++slice_iter)
                                                       printSlice(slice_iter);

  MSG("AltAlg", Msg::kDebug) << "* Printing Reco Slices completed" << endl;
}
//____________________________________________________________________________
void AltAlgSliceList::printSlice(
   std::map<int, std::vector<CandStripHandle *> >::const_iterator slice) const
{
// Method for printing the reconstructed event slices
//
  int nv = count_if(
                 slice->second.begin(), slice->second.end(), is_v_view());
  int nu = count_if(
                 slice->second.begin(), slice->second.end(), is_u_view());

  MSG("AltAlg", Msg::kDebug)
     << "[-] Slice: " << slice->first
     << " - (N(U) = " << nu << " -- N(V) = " << nv << ")" << endl;

  const std::vector<CandStripHandle *> & slice_strips = slice->second;
  std::vector<CandStripHandle *>::const_iterator strip_iter;

  for(strip_iter = slice_strips.begin();
                          strip_iter != slice_strips.end(); ++strip_iter) {

     MSG("AltAlg", Msg::kDebug)
        << " |--> Pl = " << fFmtPln(  (*strip_iter)->GetPlane()   )
        << " Str = "     << fFmtStp(  (*strip_iter)->GetStrip()   )
        << " Q = "       << fFmtQ(    (*strip_iter)->GetCharge()  )
        << " tb = "      << fFmtTime( (*strip_iter)->GetCorrBegTime() )
        << " view = " << PlaneView::AsString( (*strip_iter)->GetPlaneView() )
        << endl;
  } // strip_iter
}
//____________________________________________________________________________
void AltAlgSliceList::printSlice(
                          const std::vector<CandStripHandle *> & slice) const
{
  int nv = count_if(slice.begin(), slice.end(), is_v_view());
  int nu = count_if(slice.begin(), slice.end(), is_u_view());

  MSG("AltAlg", Msg::kDebug)
     << "[-] printing slice whose id was not specified or a strip collection"
     << " - (N(U) = " << nu << " -- N(V) = " << nv << ")" << endl;

  std::vector<CandStripHandle *>::const_iterator strip_iter;

  for(strip_iter = slice.begin(); strip_iter != slice.end(); ++strip_iter) {

     MSG("AltAlg", Msg::kDebug)
        << " |--> Pl = " << fFmtPln(  (*strip_iter)->GetPlane()   )
        << " Str = "     << fFmtStp(  (*strip_iter)->GetStrip()   )
        << " Q = "       << fFmtQ(    (*strip_iter)->GetCharge()  )
        << " tb = "      << fFmtTime( (*strip_iter)->GetCorrBegTime() )
        << " view = " << PlaneView::AsString( (*strip_iter)->GetPlaneView() )
        << endl;
  } // strip_iter
}
//____________________________________________________________________________
void AltAlgSliceList::printStripList(
                        CandStripHandleItr cshItr, const char * comment) const
{
  MSG("AltAlg", Msg::kVerbose) << comment << endl;

  MSG("AltAlg", Msg::kVerbose)
      << " |--> selected " << cshItr.GetSet()->SizeSelect()
      << " out of "        << cshItr.GetSet()->Size()
      << " entries" << endl;

  cshItr.ResetFirst();
  while( CandStripHandle * striph = cshItr.Ptr() ) {

     MSG("AltAlg", Msg::kVerbose)
        << " |--> Pl = " << fFmtPln(  striph->GetPlane()   )
        << " Str = "     << fFmtStp(  striph->GetStrip()   )
        << " Q = "       << fFmtQ(    striph->GetCharge()  )
        << " tb = "      << fFmtTime( striph->GetCorrBegTime() )
        << " view = " << PlaneView::AsString( striph->GetPlaneView() )
        << endl;

     cshItr.Next();
  }
}
//____________________________________________________________________________
void AltAlgSliceList::printStrip(const CandStripHandle * strip) const
{
   MSG("AltAlg", Msg::kVerbose)
       << "              |--> Pl = " << fFmtPln(  strip->GetPlane()   )
       << " Str = "     << fFmtStp(  strip->GetStrip()   )
       << " Q = "       << fFmtQ(    strip->GetCharge()  )
       << " tb = "      << fFmtTime( strip->GetCorrBegTime() )
       << " view = " << PlaneView::AsString( strip->GetPlaneView() )
       << endl;
}
//____________________________________________________________________________
void AltAlgSliceList::printTimeBin(
               int bin, const TH1D * time_profile, const char * comment) const
{
   MsgFormat fmtInt("%4d");

   MSG("AltAlg",Msg::kVerbose)
       << "          |--> bin: "  << fmtInt( bin )
       << ":[" << fFmtTime( (double)    time_profile->GetBinLowEdge(bin) )
       << ", " << fFmtTime( (double)  ( time_profile->GetBinLowEdge(bin) +
                                        time_profile->GetBinWidth(bin)) )
       << "], cnt: "  << fmtInt( int( time_profile->GetBinContent(bin) ))
       << " / " << comment << endl;
}
//____________________________________________________________________________
std::string AltAlgSliceList::asString(WhichCandSlice_t which_slices) const
{
  switch( which_slices ) {
  case (eNone):
      return std::string("isolated strip - no really good match found");
      break;
  case (eSame):
      return std::string("no action - this strip is @ home");
      break;
  case (eOther):
      return std::string("a single other candidate was found");
      break;
  case (eMany):
      return std::string("many 'candidates' - need to resolve the ambiguity");
      break;
  default:
      return std::string(" --- ");
  }
}
//____________________________________________________________________________
std::string AltAlgSliceList::asString(PeakFinderConf_t conf) const
{
  switch( conf ) {
  case (kDefault):        return std::string("[default config]");      break;
  case (kLowActivity):    return std::string("[low activity config]"); break;
  case (kMuSpectrometer): return std::string("[mu-spectr. config]");   break;
  default:                return std::string(" --- ");
  }
}
//____________________________________________________________________________
void AltAlgSliceList::Merge(const std::pair<int, int> & slices_to_merge,
          std::map<int, std::vector<CandStripHandle *> > & event_slices) const
{
// Method for handling the map<int, vector<CandStripHandle *> > slice map.
// Gets a pair of slice ids that have to be merged and applies this operation
// the map. All CandStripHandles are appended to the first slice of the pair
// and the other is deleted. At the end, the map is sorted in time.
//
  MSG("AltAlg", Msg::kVerbose) << "* Begin Slice Merging" << endl;

  //printSlices(event_slices,"");

  int target_slice = slices_to_merge.first;
  int source_slice = slices_to_merge.second;

  std::map<int, std::vector<CandStripHandle *> >::iterator slice_iter;

  slice_iter = event_slices.find(source_slice);

  std::vector<CandStripHandle *> source_strips = slice_iter->second;

  std::vector<CandStripHandle *>::iterator strip_iter;

  for(strip_iter = source_strips.begin();
                           strip_iter != source_strips.end(); ++strip_iter) {

     MSG("AltAlg", Msg::kVerbose)
             << "   --- Strip: " << fFmtStp( (*strip_iter)->GetStrip() )
             << " from plane: "  << fFmtPln( (*strip_iter)->GetPlane() )
             << " to be removed from slice: " << fFmtSlc( source_slice )
             << " and appended in slice: " << fFmtSlc( target_slice ) << endl;

     //-- Add cand strip handle to target slice
     event_slices[target_slice].push_back(*strip_iter);

  } // loop over source slice strips

  MSG("AltAlg", Msg::kVerbose)
                       << "   ---- Removing slice: " << source_slice << endl;
  event_slices.erase(slice_iter); // remove source slice

  MSG("AltAlg", Msg::kVerbose) << "* Slice Merging Completed" << endl;
}
//____________________________________________________________________________
void AltAlgSliceList::Split(
              int slice, std::vector<CandStripHandle *>::size_type n_separator,
          std::map<int, std::vector<CandStripHandle *> > & event_slices) const
{
// A specific case of Split(). ** The generic case is coded next **
// In this specific case, the original slice is split in just two parts and
// its vector of strips is somehow sorted to reflect this splitting.
// The first n_separator elements are assigned to the first sub-slice and the
// remaining size()-n_separator elements are assigned to the second one.
//
  MSG("AltAlg", Msg::kVerbose) << "[-] splitting slice: " << slice << endl;

  vector<CandStripHandle *>::iterator strip_iter = event_slices[slice].begin();

  advance(strip_iter, n_separator);

  std::vector<CandStripHandle *> new_slice_strips_1(n_separator);
  copy(event_slices[slice].begin(),  strip_iter,  new_slice_strips_1.begin());
  event_slices.insert(
      std::map<int, std::vector<CandStripHandle *> >::value_type(
                     nextAvailableSliceId(event_slices), new_slice_strips_1 ));

  std::vector<CandStripHandle *> new_slice_strips_2(
                                    event_slices[slice].size() - n_separator);
  copy( strip_iter, event_slices[slice].end(), new_slice_strips_2.begin());
  event_slices.insert(
      std::map<int, std::vector<CandStripHandle *> >::value_type(
                    nextAvailableSliceId(event_slices), new_slice_strips_2 ));

  event_slices.erase(slice);
}
//____________________________________________________________________________
void AltAlgSliceList::Split(
              int slice, std::vector< std::vector<CandStripHandle *> > groups,
          std::map<int, std::vector<CandStripHandle *> > & event_slices) const
{
// Method for updating the slices map to reflect the splitting of an existing
// slice. New slices are created to accomodate the strips of each elemement of
// the std::vector< std::vector<CandStripHandle *> > "groups" STL-vector...
// These elements define a partition of the strips corresponding to the input
// slice identifies.
// After creating the new slices, the original slice is deleted.
//
  MSG("AltAlg", Msg::kVerbose) << "[-] splitting slice: " << slice << endl;

  std::vector< std::vector<CandStripHandle *> >::iterator group_iter;

  for(group_iter = groups.begin(); group_iter != groups.end(); ++group_iter)

     event_slices.insert(
          std::map<int, std::vector<CandStripHandle *> >::value_type(
                           nextAvailableSliceId(event_slices), *group_iter) );

  event_slices.erase(slice);
}
//____________________________________________________________________________
int AltAlgSliceList::nextAvailableSliceId(
    const std::map<int, std::vector<CandStripHandle *> > & event_slices) const
{
  std::map<int, std::vector<CandStripHandle *> >::const_iterator slice_iter;

  slice_iter = --( event_slices.end() );

  return ( 1 + (int) (*slice_iter).first );
}
//____________________________________________________________________________
void AltAlgSliceList::initSliceFiltering(
          std::map<int, std::vector<CandStripHandle *> > & event_slices) const
{

  MSG("AltAlg", Msg::kDebug)
     << "[-] Dissolving 'mini'-slices & creating an orphan strip list" << endl;

  std::vector<CandStripHandle *> orphan_strips_list
                                            = dissolveMiniSlices(event_slices);

  if(orphan_strips_list.size() > 0) {

     giveOrphanStripsForAdoption(orphan_strips_list, event_slices);
     sortSlicesInTime(event_slices);
  }
}
//____________________________________________________________________________
std::vector<CandStripHandle *> AltAlgSliceList::dissolveMiniSlices(
          std::map<int, std::vector<CandStripHandle *> > & event_slices) const
{
// Removes bad slices from the slices-map and creates a list of orphan
// strips to be appended in nearby slices.
//
  MSG("AltAlg", Msg::kDebug)
                 << " |-->[-] Trying to eliminate some 'fake' slices" << endl;

  MSG("AltAlg", Msg::kDebug)
            << "      |---> Current number of slices: " << event_slices.size()
            << endl;

  std::vector<int>* slices_to_delete = findSlicesToBeEliminated(event_slices);

  MSG("AltAlg", Msg::kVerbose)
         << "      |---[-] Eliminate slices & copy strips to list of orphans"
         << endl;

  std::vector<CandStripHandle *> orphan_strips;
  std::vector<CandStripHandle *>::iterator strip_iter;

  std::vector<int>::iterator slice_id;

  for(slice_id = slices_to_delete->begin();
                           slice_id != slices_to_delete->end(); ++slice_id) {

      for(strip_iter = event_slices[*slice_id].begin();
              strip_iter != event_slices[*slice_id].end(); ++strip_iter)
                                      orphan_strips.push_back( *strip_iter );

       event_slices.erase( *slice_id );

       MSG("AltAlg", Msg::kVerbose)
                 << "           |--> Deleted slice : " << *slice_id  << endl;
  }

  MSG("AltAlg", Msg::kDebug)
        << "      |---> Number of slices *after* eliminating the fake ones: "
        << event_slices.size() << endl;

  delete slices_to_delete;

  return orphan_strips;
}
//____________________________________________________________________________
std::vector<int> * AltAlgSliceList::findSlicesToBeEliminated(
    const std::map<int, std::vector<CandStripHandle *> > & event_slices) const
{
  std::vector<int> * slices_to_delete = new std::vector<int>;

  std::map<int, std::vector<CandStripHandle *> >::const_iterator slice_iter;

  for(slice_iter = event_slices.begin();
                            slice_iter != event_slices.end(); ++slice_iter) {

     MSG("AltAlg", Msg::kVerbose)
         << "      |---[-] Checking slice: " << (*slice_iter).first  << endl;

     if( sliceShouldBeEliminated(slice_iter) )
                           slices_to_delete->push_back( (*slice_iter).first );
  }

  return slices_to_delete; // the calling function should delete this
}
//____________________________________________________________________________
bool AltAlgSliceList::sliceShouldBeEliminated(
   std::map<int, std::vector<CandStripHandle *> >::const_iterator slice) const
{
  bool should_be_eliminated = false;

  //**** apply the list of criteria for slice-elimination

  std::vector<bool> check_list;
  std::vector<bool>::iterator check;

  check_list.push_back( smallNumberOfHitStrips (slice)  );
  check_list.push_back( smallAmountOfCharge    (slice)  );
  check_list.push_back( noStripsInOneView      (slice)  );

  for(check = check_list.begin(); check != check_list.end(); ++check)
                      should_be_eliminated = should_be_eliminated || (*check);

  return should_be_eliminated;
}
//____________________________________________________________________________
bool AltAlgSliceList::smallNumberOfHitStrips(
   std::map<int, std::vector<CandStripHandle *> >::const_iterator slice) const
{
// check for small number of hit strips in the slice
//
  if( (int) (*slice).second.size() < fMinNoHitStrips  )  {

     MSG("AltAlg", Msg::kVerbose)
           << "           |--> Slice: " << (*slice).first
           << " will be eliminated --> too few hit strips: "
           << (*slice).second.size() << " < (MinStrips=) " << fMinNoHitStrips
           << endl;

  } else return false;

  return true;
}
//____________________________________________________________________________
bool AltAlgSliceList::smallAmountOfCharge(
   std::map<int, std::vector<CandStripHandle *> >::const_iterator slice) const
{
// check for small amount of charge in the slice
//

// double totalCharge = 0;
// for(std::vector<CandStripHandle *>::iterator strip_iter =
//   (*iter).second.begin();strip_iter != (*iter).second.end(); ++strip_iter)
//                       totalCharge += ( (*strip_iter)->GetCharge() );

  double totalCharge = accumulate( (*slice).second.begin(),
                                   (*slice).second.end(),  0.0, sum_q());

  if( totalCharge < fMinCharge )  {

      MSG("AltAlg", Msg::kVerbose)
            << "           |--> Slice: " << (*slice).first
            << " will be eliminated --> too small charge: "
            << totalCharge << " < (MinCharge=) " << fMinCharge << endl;

  } else return false;

  return true;
}
//____________________________________________________________________________
bool AltAlgSliceList::noStripsInOneView(
   std::map<int, std::vector<CandStripHandle *> >::const_iterator slice) const
{
// check for strips in both views in the slice
//
  int nv = count_if(
                 (*slice).second.begin(), (*slice).second.end(), is_v_view());
  int nu = count_if(
                 (*slice).second.begin(), (*slice).second.end(), is_u_view());

  if(nv == 0 || nu == 0) {

      MSG("AltAlg", Msg::kVerbose)
           << "           |--> Slice: " << (*slice).first
           << " will be eliminated --> 0 strips in one view: "
           << "  > N(U) = " << nu << " -- N(V) = " << nv << endl;

  } else return false;

  return true;
}
//____________________________________________________________________________
void AltAlgSliceList::giveOrphanStripsForAdoption(
                                 std::vector<CandStripHandle *> orphan_strips,
          std::map<int, std::vector<CandStripHandle *> > & event_slices) const
{
// Appends the CandStripHandles found in the orphan_strips STL vector
// to the reco. slice with the minimum distance in the z-time-tpos
// 3-D space...
//
  //-- assume that centroids do not change too much by the addition
  //   of few orhan strips, so compute them just once...

  MSG("AltAlg", Msg::kVerbose)
          << " |-->[-] Appending orphan strips to closest reco slice" << endl;

  std::vector<CandStripHandle *>::iterator strip_iter;
  for(strip_iter = orphan_strips.begin();
                           strip_iter != orphan_strips.end(); ++strip_iter) {

     MSG("AltAlg", Msg::kVerbose)
       << "      |--[-] Looking a home for orphan strip "
       << (*strip_iter)->GetUidInt()
       << " (pl = "  << fFmtPln(  (*strip_iter)->GetPlane() )
       << ", str = " << fFmtStp(  (*strip_iter)->GetStrip() )
       << ", t = "   << fFmtTime( (*strip_iter)->GetCorrBegTime() )
       << ")" << endl;

     int islice = findBestSliceToAdoptAStrip(
                       fOrphanStripsQWeight, fOrphanStripsPlaneWindow,
                       fOrphanStripsTimeWindow, *strip_iter, &event_slices );

     event_slices[islice].push_back( *strip_iter );

     MSG("AltAlg", Msg::kVerbose)
                     << "          |--> Adopted by slice "  << islice << endl;
  } // loop over orphan strips
}
//____________________________________________________________________________
int AltAlgSliceList::findBestSliceToAdoptAStrip(
            bool qweight, int dpln, double dt, CandStripHandle * orphan_strip,
          std::map<int, std::vector<CandStripHandle *> > * event_slices) const
{
// Finds the best slice (amongst the already existing ones) to adopt the input
// strip from the list of orphan ones.
// The best slice is the one which has the biggest number of hit strips (or
// the biggest amount of charge, if qweight is turned on) inside a temporal
// and spatial window centered in the orphan strip.
// If none is found, the method relaxes the criteria (by increasing the window
// size) and calls itself again.
//
  int    best_slice = -1;
  double nq_max = 0, nq;

  std::map<int, std::vector<CandStripHandle *> >::iterator slice_iter;

  for(slice_iter = event_slices->begin();
                           slice_iter != event_slices->end(); ++slice_iter) {

      //-- no weighting by charge - just count number of hit strips
      if(qweight == 0)
        nq = (double) count_if(
             (*slice_iter).second.begin(), (*slice_iter).second.end(),
                                    is_in_tz_window(orphan_strip, dt, dpln) );
      //-- weighting by charge
      else {

        std::vector<CandStripHandle *>::iterator part_iter = partition(
           (*slice_iter).second.begin(), (*slice_iter).second.end(),
                                    is_in_tz_window(orphan_strip, dt, dpln) );
        nq = accumulate(
                       (*slice_iter).second.begin(), part_iter, 0.0, sum_q());
      }

      if(nq > nq_max) { nq_max = nq; best_slice = (*slice_iter).first; }

      //MSG("AltAlg", Msg::kVerbose) << "SLICE " << (*slice_iter).first
      //   << " ---> nq = " << nq << endl;
  }// slice_iter

  // No slice was found - increase the window size and try again...
  // Obviously, this can go many levels deep... If, in the next call, none is
  // found as well, it will relax the criteria some more and call itself again
  // and then possibly again... until the 'best' slice is found...
  // In this case the 'best' might not be good enough, but for now all we
  // want is to assign all orphan strips to the best available reco slice and
  // we will check later if they have substructure and need to be further
  // sliced...

  if(best_slice == -1) {
     MSG("AltAlg", Msg::kVerbose)
       << "          |--> No slice was found to adopt strip - relaxing limits"
       << " (dpln < " << fFmtPln(  min(2*dpln, 80) )
       << " - dt << " << fFmtTime( 2*dt   ) << ") & re-trying" << endl;

     return findBestSliceToAdoptAStrip(
                  qweight, min(2*dpln, 80), 2*dt, orphan_strip, event_slices);
  }
  return best_slice;
}
//____________________________________________________________________________
void AltAlgSliceList::sliceMerger(
          std::map<int, std::vector<CandStripHandle *> > & event_slices) const
{
// Sometimes it might be possible that a real event is split between 2 or more
// reconstructed slices...
// So, if slices are very close in time and span successive (or overlapping)
// areas of the detector, just merge them.
//
  MSG("AltAlg", Msg::kVerbose)
                          << "[-] Trying to spot slices to be merged" << endl;

  pair<int, int> slices_to_merge;

  while( (slices_to_merge = findSlicesToMerge(event_slices)).first != -1 )
                                         Merge(slices_to_merge, event_slices);
}
//____________________________________________________________________________
std::pair<int, int> AltAlgSliceList::findSlicesToMerge(
          std::map<int, std::vector<CandStripHandle *> > & event_slices) const
{
// Returns the first pair of slices that should be merged...
//
  std::map<int, std::vector<CandStripHandle *> >::iterator slice_iter_0;
  std::map<int, std::vector<CandStripHandle *> >::iterator slice_iter_1;

  pair<int, int> slices_to_merge(-1,-1);

  // The following 2-ble loop, loops over all *different* slice-pairs once.
  // Each pair will be checked against the 'merging' criteria.

  for(slice_iter_0 = event_slices.begin();
                        slice_iter_0 != event_slices.end(); ++slice_iter_0) {

      for( ++(slice_iter_1 = slice_iter_0);
                        slice_iter_1 != event_slices.end(); ++slice_iter_1) {

          MSG("AltAlg", Msg::kVerbose) << " |--[-] examining slices "
                         << fFmtSlc( (*slice_iter_0).first ) << " and "
                         << fFmtSlc( (*slice_iter_1).first ) << endl;

           bool should_be_merged = checkIfSlicesShouldBeMerged(
                         (*slice_iter_0).second, (*slice_iter_1).second );

           if(should_be_merged) {
                slices_to_merge.first  = (*slice_iter_0).first;
                slices_to_merge.second = (*slice_iter_1).first;

                // return this pair *now*. If two slices are merged, then
                // previous checks are potentially invalidated.
                // So merge the slices & then start checking from scratch.
                return slices_to_merge;
           }

      } // slice_iter_1
  } // slice_iter_0

  MSG("AltAlg", Msg::kVerbose)
                 << "* Finished trying to spot clusters to be merged" << endl;

  return slices_to_merge;
}
//____________________________________________________________________________
bool AltAlgSliceList::checkIfSlicesShouldBeMerged(
                               std::vector<CandStripHandle *> & slice_0,
                               std::vector<CandStripHandle *> & slice_1) const
{
// Applies all the slice merging criteria...

  std::vector<bool> check_list;
  std::vector<bool>::iterator check;

  check_list.push_back( checkSliceDistanceInTimeAndSpace( slice_0, slice_1 ) );
  // --- push_back() more checks later...
  //

  bool merge_them = true;
            for(check = check_list.begin();
                  check != check_list.end(); ++check)
                                         merge_them = merge_them && (*check);
  return merge_them;
}
//____________________________________________________________________________
bool AltAlgSliceList::checkSliceDistanceInTimeAndZ(
                               std::vector<CandStripHandle *> & slice_0,
                               std::vector<CandStripHandle *> & slice_1) const
{
// A simple test to check whether two slices should be merged on the basis of
// small time & z difference

  // check the time difference between the pulse height peaks
  if( checkPeakTimeDifference(slice_0, slice_1) ) {

    // check the z-difference between the pulse height peaks
    if( checkPeakZDifference(slice_0, slice_1) ) {

        // check the z-diff. between the latest hits of the earliest slice
        // and the earliest hits of the latest slice....
        if( checkLeadingTrailingEdgeDistance(slice_0, slice_1) ) {

            MSG("AltAlg", Msg::kVerbose)
                              << "     |--[*] Slices will be merged" << endl;
            return true;
        }
    }
  }
  return false;
}
//____________________________________________________________________________
bool AltAlgSliceList::checkSliceDistanceInTimeAndSpace(
                               std::vector<CandStripHandle *> & slice_0,
                               std::vector<CandStripHandle *> & slice_1) const
{
// A simple test to check whether two slices should be merged on the basis of
// small time & z,u,v difference

  // check the time difference between the pulse height peaks
  if( checkPeakTimeDifference(slice_0, slice_1) ) {

    // check the z-difference between the pulse height peaks
    if( checkPeakZDifference(slice_0, slice_1) ) {

        // check the z-diff. between the latest hits of the earliest slice
        // and the earliest hits of the latest slice....
        if( checkLeadingTrailingEdgeDistance(slice_0, slice_1) ) {

            // check uv difference
            if( checkPeakUVDifference(slice_0, slice_1) ) {
                  MSG("AltAlg", Msg::kVerbose)
                              << "     |--[*] Slices will be merged" << endl;
                  return true;
            }
        }
    }
  }
  return false;
}
//____________________________________________________________________________
bool AltAlgSliceList::checkPeakTimeDifference(
                         const std::vector<CandStripHandle *> & slice_0,
                         const std::vector<CandStripHandle *> & slice_1) const
{
// Checks the difference between the charge weighted average times of the two
// input slices.
// Returns true if dt < fTimeDiffBetweenPeaks

  double q0  = accumulate(slice_0.begin(), slice_0.end(), 0.0, sum_q());
  double qt0 = accumulate(slice_0.begin(), slice_0.end(), 0.0, sum_qtime());

  double q1  = accumulate(slice_1.begin(), slice_1.end(), 0.0, sum_q());
  double qt1 = accumulate(slice_1.begin(), slice_1.end(), 0.0, sum_qtime());

  //remove this eventually
  if(q0 == 0 || q1 == 0) {
    MSG("AltAlg", Msg::kWarning) << "Zero charge in one view" << endl;
    return false;
  }

  assert(q0 > 0); // slices that are striped out of their strips should
  assert(q1 > 0); // be removed from the slice map

  double t0 = qt0 / q0;
  double t1 = qt1 / q1;

  MSG("AltAlg", Msg::kVerbose)
            << "     |--> peak time difference = " << fFmtTime(t0-t1) << endl;

  return ( fabs(t0-t1) < fTimeDiffBetweenPeaks );
}
//____________________________________________________________________________
bool AltAlgSliceList::checkPeakZDifference(
                         const std::vector<CandStripHandle *> & slice_0,
                         const std::vector<CandStripHandle *> & slice_1) const
{
// Checks the difference between the charged-weighted longitudinal position
// of the input slices.
// Returns true if dz < fZDiffBetweenPeaks

  double q0  = accumulate(slice_0.begin(), slice_0.end(), 0.0, sum_q());
  double qz0 = accumulate(slice_0.begin(), slice_0.end(), 0.0, sum_qz());

  double q1  = accumulate(slice_1.begin(), slice_1.end(), 0.0, sum_q());
  double qz1 = accumulate(slice_1.begin(), slice_1.end(), 0.0, sum_qz());

  //remove this eventually
  if(q0 == 0 || q1 == 0) {
    MSG("AltAlg", Msg::kWarning) << "Zero charge in one view" << endl;
    return false;
  }

  assert(q0 > 0); // slices that are striped out of their strips should
  assert(q1 > 0); // be removed from the slice map

  double z0 = qz0 / q0;
  double z1 = qz1 / q1;

  MSG("AltAlg", Msg::kVerbose) << "     |--> d<z> = " << z0-z1 << endl;

  return ( fabs(z0-z1) < fZDiffBetweenPeaks );
}
//____________________________________________________________________________
bool AltAlgSliceList::checkPeakUVDifference(
                         const std::vector<CandStripHandle *> & slice_0,
                         const std::vector<CandStripHandle *> & slice_1) const
{
// Checks the difference between the charged-weighted transverse position
// of the input slices.
// Returns true if du < fUVDiffBetweenPeaks && dv < fUVDiffBetweenPeaks

  //-- u-v partition
  std::vector<CandStripHandle *> slice_0_cp( slice_0.size() );
  std::vector<CandStripHandle *> slice_1_cp( slice_1.size() );

  copy( slice_0.begin(), slice_0.end(), slice_0_cp.begin() );
  copy( slice_1.begin(), slice_1.end(), slice_1_cp.begin() );

  std::vector<CandStripHandle *>::iterator part_iter_0 = partition(
                        slice_0_cp.begin(), slice_0_cp.end(), is_v_view() );
  std::vector<CandStripHandle *>::iterator part_iter_1 = partition(
                        slice_1_cp.begin(), slice_1_cp.end(), is_v_view() );

  // slice 0
  double q_v_0 = accumulate(slice_0_cp.begin(), part_iter_0, 0.0, sum_q());
  double qtpos_v_0 = accumulate(
                        slice_0_cp.begin(), part_iter_0, 0.0, sum_qtpos());

  double q_u_0 = accumulate(part_iter_0, slice_0_cp.end(),   0.0, sum_q());
  double qtpos_u_0 = accumulate(
                        part_iter_0, slice_0_cp.end(),   0.0, sum_qtpos());

  // slice 1
  double q_v_1 = accumulate(slice_1_cp.begin(), part_iter_1, 0.0, sum_q());
  double qtpos_v_1 = accumulate(
                        slice_1_cp.begin(), part_iter_1, 0.0, sum_qtpos());

  double q_u_1 = accumulate(part_iter_1, slice_1_cp.end(),   0.0, sum_q());
  double qtpos_u_1 = accumulate(
                        part_iter_1, slice_1_cp.end(),   0.0, sum_qtpos());

  // remove this eventually
  if(q_v_0 == 0 || q_u_0 == 0 || q_v_1 == 0 || q_u_1 == 0) {
    MSG("AltAlg", Msg::kWarning) << "Zero charge in one view" << endl;
    return false;
  }
  assert(q_v_0 > 0); // slices that are striped out of their strips should
  assert(q_u_0 > 0); // be removed from the slice map
  assert(q_v_1 > 0);
  assert(q_u_1 > 0);

  double u0 = qtpos_u_0 / q_u_0;
  double v0 = qtpos_v_0 / q_v_0;
  double u1 = qtpos_u_1 / q_u_1;
  double v1 = qtpos_v_1 / q_v_1;

  MSG("AltAlg", Msg::kVerbose)
         << "     |--> d<u> = " << u0-u1 << " and d<v> = " << v0 - v1 << endl;

  return ( fabs(v0-v1) < fUVDiffBetweenPeaks &&
                                          fabs(u0-u1) < fUVDiffBetweenPeaks );
}
//____________________________________________________________________________
bool AltAlgSliceList::checkLeadingTrailingEdgeDistance(
                               std::vector<CandStripHandle *> & slice_0,
                               std::vector<CandStripHandle *> & slice_1) const
{
// Checks the longitudinal (z) distance between the leading edge of the latest
// slice and the trailing edge of the earliest one...
// Returns true if the distance is smaller than fZDiffBetweenEnds
// This is hoped to prevent slice breaking when use in conjunction with
// temporal cuts.
//
  double q0  = accumulate(slice_0.begin(), slice_0.end(), 0.0, sum_q());
  double qt0 = accumulate(slice_0.begin(), slice_0.end(), 0.0, sum_qtime());

  double q1  = accumulate(slice_1.begin(), slice_1.end(), 0.0, sum_q());
  double qt1 = accumulate(slice_1.begin(), slice_1.end(), 0.0, sum_qtime());

  // remove this eventually
  if(q0 == 0 || q1 == 0) {
    MSG("AltAlg", Msg::kWarning) << "Zero charge in one view" << endl;
    return false;
  }

  assert(q0 > 0); // slices that are striped out of their strips should
  assert(q1 > 0); // be removed from the slice map

  double t0 = qt0 / q0;
  double t1 = qt1 / q1;

  double z_earliest, z_latest;

  if(t0 < t1) {
      z_earliest = averageSliceZ( slice_0, t0, maxSliceTime(slice_0) );
      z_latest   = averageSliceZ( slice_1, minSliceTime(slice_1), t1 );
  } else {
      z_earliest = averageSliceZ( slice_1, t1, maxSliceTime(slice_1) );
      z_latest   = averageSliceZ( slice_0, minSliceTime(slice_0), t0 );
  }

  MSG("AltAlg", Msg::kVerbose)
     << "     |--> leading - trailing edge z difference = "
                                             << z_earliest - z_latest << endl;

  if(fabs(z_earliest-z_latest) < fZDiffBetweenEnds)  return true;
  else return false;
}
//____________________________________________________________________________
double AltAlgSliceList::averageSliceZ(
       std::vector<CandStripHandle *> & slice, double tmin, double tmax) const
{
// Computes a (charge weighted) average longitudinal position of the input
// strips for which time e [tmin, tmax]
//
  std::vector<CandStripHandle *>::iterator strip_iter = partition(
                slice.begin(), slice.end(), is_in_t_window_noref(tmin, tmax));

  double q  = accumulate(slice.begin(), strip_iter, 0.0, sum_q() );
  double qz = accumulate(slice.begin(), strip_iter, 0.0, sum_qz());

  // remove this eventually
  if(q==0) {
    MSG("AltAlg", Msg::kWarning) << "Zero charge in one view" << endl;
    return -999;
  }

  assert(q>0);

  return qz/q;
}
//____________________________________________________________________________
double AltAlgSliceList::minSliceTime(
                                 std::vector<CandStripHandle *> & slice) const
{
// Returns the earliest strip time in the input vector of strips

  sort(slice.begin(), slice.end(), min_t());

  std::vector<CandStripHandle *>::iterator strip_iter = slice.begin();

  return (*strip_iter)->GetCorrBegTime();
}
//____________________________________________________________________________
double AltAlgSliceList::maxSliceTime(
                                 std::vector<CandStripHandle *> & slice) const
{
// Returns the latest strip time in the input vector of strips

  sort(slice.begin(), slice.end(), min_t());

  std::vector<CandStripHandle *>::iterator strip_iter = --(slice.end());

  return (*strip_iter)->GetCorrBegTime();
}
//____________________________________________________________________________
void AltAlgSliceList::eventClustering(
                std::map<int, std::vector<CandStripHandle *> > & event_slices)
{
// Main method for k-Means - based 3-D event clustering. Initiates single
// k-Means clustering method steps until the algorithm has converged.
//
  MSG("AltAlg", Msg::kDebug)
          << "[-] Applying 3-D (c*time, tpos, z) k-Means clustering"  << endl;

  bool converged   = false;
  fkMeansIteration = 0;

  while(!converged) {

      MSG("AltAlg", Msg::kDebug)
                   << " |--[-] k-Means clustering has not converged yet - "
                   << "next iteration:" << endl;

      converged = singleKMeansIteration(event_slices);
  }

  // refilter to catch/remove any slice that has become too 'thin'...
  initSliceFiltering(event_slices);
}
//____________________________________________________________________________
bool AltAlgSliceList::singleKMeansIteration(
                std::map<int, std::vector<CandStripHandle *> > & event_slices)
{
// Performs a single iteration of a K-Means clustering algorithm on the
// std::map<int, std::vector<CandStripHandle *> > of reco. slices
// Clustering is performed in time-z-tpos 3D space and strips are
// re-arranged between existing slices as appropriate.
// Returns "true" if the algorithm has converged.
//
  bool is_convergent = true;

  // The outer loop can go on for ever if for some reason something goes wrong.
  // Plan your escape route, if such an ill case appears... Force convergence:
  if( ++fkMeansIteration > 30 ) return true;

  fUpdateCentroids = true; // make known that centroids are not updated

  std::map<int, Centroid_t>  centroids;

  std::map<int, std::vector<CandStripHandle *> >::iterator slice_iter;
  for(slice_iter = event_slices.begin();
                            slice_iter != event_slices.end(); ++slice_iter) {

     int source = slice_iter->first; // source slice for following strips

     // check if all strips feel like home in this slice...
     for(std::vector<CandStripHandle *>::size_type istrip = 0;
                              istrip < slice_iter->second.size(); istrip++) {

         CandStripHandle * current_strip = slice_iter->second[istrip];

         pair<int, CandStripHandle *> slc_strip(source, current_strip);

         // try to update centroids - if it fails because it has to re-filter
         // the slices, filter them and 'admit' you have not converged yet...

         if ( !updateCentroids(centroids, event_slices) ) return false;

         // find out which is the best slice, amongst the existing ones
         // to host this strip...

         int target = findBestSliceToHostStrip(
                                         slc_strip, centroids, event_slices);

         if( source != target ) {
            //-- k-Means clustering is not convergent yet...
            is_convergent = false;

            //-- we rellocate a strip - recompute centroids if needed again
            fUpdateCentroids = true;

            //-- rellocate strip
            rellocateStrip(source, target, current_strip, event_slices);
         }

     } // strips
  } // slice_iter

  return is_convergent;
}
//____________________________________________________________________________
bool AltAlgSliceList::needToFilterStrips(
    const std::map<int, std::vector<CandStripHandle *> > & event_slices) const
{
// Somethimes it is possible that strip exchanges during the 3-D clustering
// will (almost or completely) strip out a slice from strips.
// In this case, computing slice centroids will fail because Sum{Q}=0 in
// at least one of the views.
// If this is happening then these small slices have to be filtered out.
// In this method I check if this is necessary at *this* point...

  std::map<int, std::vector<CandStripHandle *> >::const_iterator slice_iter;

  for(slice_iter = event_slices.begin();
                            slice_iter != event_slices.end(); ++slice_iter) {

      std::vector<CandStripHandle *> slice_strips = slice_iter->second;

      if( ZeroChargeInAView(slice_strips) ) {

        MSG("AltAlg", Msg::kVerbose)
           << " |--[-] Oops! At least one of slices was striped out during"
           << " strip exchanges -- run the filter to clear them out" << endl;

        return true;
      }
  }//slices

  return false;
}
//____________________________________________________________________________
bool AltAlgSliceList::ZeroChargeInAView(
                           const std::vector<CandStripHandle *> & slice) const
{
  //-- u-v partition
  std::vector<CandStripHandle *> slice_cp( slice.size() );

  copy( slice.begin(), slice.end(), slice_cp.begin() );

  std::vector<CandStripHandle *>::iterator part_iter = partition(
                             slice_cp.begin(), slice_cp.end(), is_v_view() );

  //-- accumulate
  double sum_q_v = accumulate(slice_cp.begin(), part_iter, 0.0, sum_q());
  double sum_q_u = accumulate(part_iter, slice_cp.end(),   0.0, sum_q());

  //--check
  if(sum_q_v <= 0 || sum_q_u <= 0) return true;

  return false;
}
//____________________________________________________________________________
bool AltAlgSliceList::updateCentroids(
                                        std::map<int, Centroid_t> & centroids,
                std::map<int, std::vector<CandStripHandle *> > & event_slices)
{
// Try to update centroids... if its needed... and if you can...
// Returns false if it had to re-filter the slices.

  if(fUpdateCentroids) {

     if(needToFilterStrips(event_slices)) {
        initSliceFiltering(event_slices);
        return false;
     } else {
        centroids = computeSlicesCentroid(event_slices);
        fUpdateCentroids = false;
     }
  }
  return true;
}
//____________________________________________________________________________
std::map<int, Centroid_t> AltAlgSliceList::computeSlicesCentroid(
    const std::map<int, std::vector<CandStripHandle *> > & event_slices) const
{
// Method for computing the centroids of the input reconstructed slices.
//
  std::map<int, Centroid_t> clusters_centroid;

  std::map<int, std::vector<CandStripHandle *> >::const_iterator slice_iter;

  for(slice_iter = event_slices.begin();
                            slice_iter != event_slices.end(); ++slice_iter) {

      std::vector<CandStripHandle *> slice_strips = slice_iter->second;

      Centroid_t centroid = computeSliceCentroid(slice_strips);

      clusters_centroid.insert(
            std::map<int, Centroid_t>::value_type(
                                              slice_iter->first, centroid) );
  }

  return clusters_centroid;
}
//____________________________________________________________________________
Centroid_t AltAlgSliceList::computeSliceCentroid(
                           const std::vector<CandStripHandle *> & slice) const
{
// Method for computing the centroids of the input reconstructed slice.
//
  Centroid_t centroid;

  //-- u-v partition
  std::vector<CandStripHandle *> slice_cp( slice.size() );

  copy( slice.begin(), slice.end(), slice_cp.begin() );

  std::vector<CandStripHandle *>::iterator part_iter = partition(
                             slice_cp.begin(), slice_cp.end(), is_v_view() );

  //-- accumulate
  double sum_q_v = accumulate(
                          slice_cp.begin(), part_iter, 0.0, sum_q());
  double sum_qtpos_v = accumulate(
                      slice_cp.begin(), part_iter, 0.0, sum_qtpos());
  double sum_qz_v = accumulate(
                         slice_cp.begin(), part_iter, 0.0, sum_qz());
  double sum_qtime_v = accumulate(
                      slice_cp.begin(), part_iter, 0.0, sum_qtime());

  double sum_q_u = accumulate(
                          part_iter, slice_cp.end(),   0.0, sum_q());
  double sum_qtpos_u = accumulate(
                      part_iter, slice_cp.end(),   0.0, sum_qtpos());
  double sum_qz_u = accumulate(
                         part_iter, slice_cp.end(),   0.0, sum_qz());
  double sum_qtime_u = accumulate(
                      part_iter, slice_cp.end(),   0.0, sum_qtime());

  assert(sum_q_u > 0);
  assert(sum_q_v > 0);

  centroid.u  = sum_qtpos_u    / sum_q_u;
  centroid.v  = sum_qtpos_v    / sum_q_v;
  centroid.zu = sum_qz_u       / sum_q_u;
  centroid.zv = sum_qz_v       / sum_q_v;
  centroid.tu = sum_qtime_u    / sum_q_u;
  centroid.tv = sum_qtime_v    / sum_q_v;

  return centroid;
}
//____________________________________________________________________________
int AltAlgSliceList::findBestSliceToHostStrip(
                         const pair<int, CandStripHandle *> & slc_strip,
                            const std::map<int, Centroid_t> & centroids,
    const std::map<int, std::vector<CandStripHandle *> > & event_slices) const
{
// Find which of the existing reconstructed slices is the best one for hosting
// the input strip.

  int target=0;

  std::vector<int> alt_slices = checkForAlternativeSlices(
                                                     slc_strip, event_slices);

  WhichCandSlice_t which_slice = characterizeCandidateSlices(
                                                 slc_strip.first, alt_slices);
  MSG("AltAlg", Msg::kVerbose)
               << "     |--[-] slice = " << fFmtSlc(slc_strip.first) << " -->"
               << " plane = "  << fFmtPln( slc_strip.second->GetPlane() )
               << ", strip = " << fFmtStp( slc_strip.second->GetStrip() )
               << " : " << asString(which_slice) << endl;

  switch( which_slice ) {
     case (eNone):
           // There is no really good slice candidate... This happens
           // with fairly isolated hit strips... Not much you can do
           // here & at this stage. Assign the strip to the cluster
           // with the closest centroid...
           target = getClosestCentroid(slc_strip.second, centroids);
           break;
     case (eSame):
           // There is a single slice candidate and it is the same
           // with the one that already owns the strip. No rellocation.
           target = slc_strip.first;
           break;
     case (eOther):
           // There is a single candidate that is different from the
           // one that already owns the strip. Mark for rellocation.
           target = alt_slices[0];
           break;
     case (eMany):
           target = selectBestCandidateSlice(
                                         slc_strip, alt_slices, event_slices);
     break;
  }

  return target;
}
//____________________________________________________________________________
std::vector<int> AltAlgSliceList::checkForAlternativeSlices(
                         const pair<int, CandStripHandle *> & slc_strip,
    const std::map<int, std::vector<CandStripHandle *> > & event_slices) const
{
// Make a list with all slices that could possibly host the input strip

  std::vector<int> alt_slices;
  std::map<int, std::vector<CandStripHandle *> >::const_iterator slice_iter;

  int source = slc_strip.first;
  CandStripHandle * current_strip = slc_strip.second;

  // find slices that have hit strips in a time, z, tpos window centered at the
  // current strip...

  for(slice_iter = event_slices.begin();
                            slice_iter != event_slices.end(); ++slice_iter) {

      int n = count_if(
               slice_iter->second.begin(), slice_iter->second.end(),
               is_in_tztpos_window(current_strip, fKMeansTimeWindow,
                                   fKMeansPlaneWindow, fKMeansTPosWindow) );
      if( slice_iter->first == source ) n--;

      if(n > 0) alt_slices.push_back( slice_iter->first );
  }

  // if none is found, relax the criteria and check in a t, z window centered
  // at the current strip...

  if( alt_slices.size() == 0 ) {

       for(slice_iter = event_slices.begin();
                            slice_iter != event_slices.end(); ++slice_iter) {

            int n = count_if(
                         slice_iter->second.begin(), slice_iter->second.end(),
                         is_in_tz_window(current_strip,
                                    fKMeansTimeWindow, 2*fKMeansPlaneWindow) );

            if( slice_iter->first == source ) n--;

            if(n > 0) alt_slices.push_back( slice_iter->first );
       }
  }

  // return whatever was found - if it is still empty it will be handled later
  return alt_slices;
}
//____________________________________________________________________________
WhichCandSlice_t AltAlgSliceList::characterizeCandidateSlices(
                        int source, const std::vector<int> & alt_slices) const
{
  if( alt_slices.size() == 0 )     return eNone;
  else if ( alt_slices.size() == 1 ) {

       if(alt_slices[0] == source) return eSame;
       else                        return eOther;

  } else                           return eMany;

  return eSame;
}
//____________________________________________________________________________
int AltAlgSliceList::getClosestCentroid(
   CandStripHandle * strip, const std::map<int, Centroid_t> & centroids) const
{
  std::multimap<double, int, less<double> > slices_dt_sorted;
  std::multimap<double, int, less<double> > slices_dz_sorted;

  std::map<int, Centroid_t>::const_iterator cd_iter;
  for(cd_iter = centroids.begin(); cd_iter != centroids.end(); ++cd_iter) {

     double dt = ((*cd_iter).second.tu+(*cd_iter).second.tv)/2. -
                                                             strip->GetCorrBegTime();

     double dz = ((*cd_iter).second.zu+(*cd_iter).second.zv)/2. -
                                                             strip->GetZPos();

     slices_dt_sorted.insert(std::multimap<double, int>::value_type(
                                                 fabs(dt), (*cd_iter).first));
     slices_dz_sorted.insert(std::multimap<double, int>::value_type(
                                                 fabs(dz), (*cd_iter).first));
  }

  // Use this & similar information to decide what is the 'closest' centroid...
  // For now just:

  MSG("AltAlg", Msg::kVerbose)
       << "         |---> slice with 'closest' centroid: "
                               << (*(slices_dt_sorted.begin())).second << endl;

  return (*(slices_dt_sorted.begin())).second;
}
//____________________________________________________________________________
int AltAlgSliceList::selectBestCandidateSlice(
                          const pair<int, CandStripHandle *> & slc_strip,
                                    const std::vector<int> & cand_slices,
     const std::map<int, std::vector<CandStripHandle *> > & event_slices) const
{
// Selects the "best" slice to host a strip, out of a list of good candidates.
// What follows is a first approach

  int source = slc_strip.first;
  CandStripHandle * strip = slc_strip.second;

  //-- check how many of these candidate slices have other strips in a t,tpos
  //   window (centered at the strip) at the two succesive planes of the same,
  //   orientation (1 upstream & 1 downstream)

  std::vector<int> slice_list;
  std::vector<int>::const_iterator slice_id_iter;
  std::map<int, std::vector<CandStripHandle *> >::const_iterator slice_iter;

  for(slice_id_iter = cand_slices.begin();
                       slice_id_iter != cand_slices.end(); ++slice_id_iter) {

      slice_iter = event_slices.find(*slice_id_iter);

      int n = count_if( slice_iter->second.begin(), slice_iter->second.end(),
                        is_in_tztpos_window(strip, fKMeansTimeWindow,
                                                     2, fKMeansTPosWindow) );

      if( *slice_id_iter == source ) n--; // do not count yourself...

      if(n > 0) slice_list.push_back( *slice_id_iter );
  }

  //-- if there was only one then return this as the best candidate
  if(slice_list.size() == 1) return slice_list[0];

  else {
    //-- select the slice which has more strips in the list of the N closest
    //   neighbors...

    std::multimap<double, int, less<double> > slices_dist_sorted =
                           MakeDistanceMap(strip, cand_slices, event_slices);

    //-- find the slice with the more occurences in the N closest positions

    int best_cand_slice = moreOccurences(slices_dist_sorted, 10);

    MSG("AltAlg", Msg::kVerbose)
       << "         |---> best candidate slice: " << best_cand_slice << endl;

    return best_cand_slice;
  }

  MSG("AltAlg", Msg::kVerbose)
       << "         |---> keep strip in slice: " << source << endl;

  return source; // if  something go wrong - just leave the strip in its place
}
//____________________________________________________________________________
std::multimap<double, int, less<double> > AltAlgSliceList::MakeDistanceMap(
          CandStripHandle * strip, const std::vector<int> & cand_slices,
    const std::map<int, std::vector<CandStripHandle *> > & event_slices) const
{
// Create a map of 'distance' -> 'slice_id', sorted in distance.
// The 'distance' is computed with respect to the input CandStripHandle.

  std::multimap<double, int, less<double> > slices_dist_sorted;

  std::vector<int>::const_iterator slice_id_iter;
  std::vector<CandStripHandle *>::const_iterator strip_iter;
  std::map<int, std::vector<CandStripHandle *> >::const_iterator slice_iter;

  for(slice_id_iter = cand_slices.begin();
                       slice_id_iter != cand_slices.end(); ++slice_id_iter) {

     slice_iter = event_slices.find(*slice_id_iter);

     for(strip_iter = slice_iter->second.begin();
                      strip_iter != slice_iter->second.end(); ++strip_iter) {

        double cdt = 3.0e8 * ( (*strip_iter)->GetCorrBegTime() -
                                                          strip->GetCorrBegTime() );
        double dz  = (*strip_iter)->GetZPos() - strip->GetZPos();

        double d2 = cdt*cdt+dz*dz; // neglect tpos for now - keep 2 views

        slices_dist_sorted.insert(
                 std::multimap<double, int>::value_type(d2, *slice_id_iter));
     } // slice strips
  }// slices

  return slices_dist_sorted;
}
//____________________________________________________________________________
int AltAlgSliceList::moreOccurences(
                            std::multimap<double, int> slice_ids, int n) const
{
// Finds the slice with the more occurences in the n first slots of the multimap
// The 'double' argument of the multimap is a measure of 'distance' in some
// variable. It therefore represents a "distance" -> "slice id" mapping.
// The multimap<Key, Data, Compare, Alloc> must have Compare = less<double>
// It is a multimap, rather than a map, because, in principle, >1 slices might
// have the same key (generalized distance from something).
// If n exceeds the multimap size it is set equal to it.
// If more than one slice is found then it returns the one of them which has
// the entry with the smallest distance.

  std::vector<int> slc_ids;
  std::vector<int> unique_slc_ids;

  std::map<double, int>::iterator id_iter;
  std::map<double, int>::iterator id_iter_last = slice_ids.begin();
  advance(id_iter_last, min(n, (int) slice_ids.size()) );

  for(id_iter = slice_ids.begin(); id_iter != id_iter_last; ++id_iter) {

       slc_ids.push_back( (*id_iter).second );

       if(! count_if( unique_slc_ids.begin(),  unique_slc_ids.end(),
                 bind2nd(equal_to<int>(), (*id_iter).second) ) )
                                unique_slc_ids.push_back( (*id_iter).second );
  }

  std::map<int, int, greater<int> > occurences;
  std::vector<int>::iterator id;

  for(id = unique_slc_ids.begin(); id != unique_slc_ids.end(); ++id)
            occurences.insert( std::map<int, int>::value_type(
                        count_if( slc_ids.begin(),  slc_ids.end(),
                                       bind2nd(equal_to<int>(), *id) ), *id));

  std::multimap<int, int>::iterator max_occurences =
                          max_element( occurences.begin(), occurences.end() );

  return max_occurences->second;
}
//____________________________________________________________________________
void AltAlgSliceList::rellocateStrip(
                        int source, int target, CandStripHandle * strip,
          std::map<int, std::vector<CandStripHandle *> > & event_slices) const
{
// Rellocates strip from slice: source to slice: target
//
  MSG("AltAlg", Msg::kVerbose)
           << "         |--[-] strip rellocation is attempted now..." << endl;

  MSG("AltAlg", Msg::kVerbose)
             << "             |--> source: " << fFmtSlc( source )
             << " [elements = " << event_slices[source].size() << "]"
             << " --> target: " << fFmtSlc( target )
             << " [elements = " << event_slices[target].size() << "]" << endl;

  event_slices[target].push_back( strip );

  event_slices[source].erase( find_if( event_slices[source].begin(),
                                       event_slices[source].end(),
                                       same_strip(strip) ));

  MSG("AltAlg", Msg::kVerbose)
         << "             |--> strip rellocation done..."
         << " [src slice elements] = " << event_slices[source].size()
         << " [tgt slice elements] = " << event_slices[target].size() << endl;
}
//____________________________________________________________________________
void AltAlgSliceList::sliceSplitter(
            std::map<int, std::vector<CandStripHandle *> > & /*event_slices*/)
{
// Minimal Spanning Tree  code goes here...
// Utility methods go below...
}
//____________________________________________________________________________
void AltAlgSliceList::assignMuonSpectrHits2Slices(CandStripHandleItr cshItr,
                std::map<int, std::vector<CandStripHandle *> > & event_slices)
{
// Method for assigning muon spectrometer strips to upstream detector slices

  // Get slices in the muon spectrometer

  MSG("AltAlg", Msg::kVerbose)
                     << "[-] Slicing the muon spectrometer activity" << endl;

  std::map<int, std::vector<CandStripHandle *> > mu_slices =
                                                     getMuSpecSlices(cshItr);
  printSlices(mu_slices, "muon spectrometer slices");

  // Try to append muon spectrometer slices to existing slice (if one exists
  // quite close in space & time). Otherwise, form new slices or append to
  // a rubbish slice made exclusively from unmatched mu-spectr. strips...
  // to slices.

  MSG("AltAlg", Msg::kVerbose)
     << "[-] examine mu-spectrometer & upstream detector slices for matches"
     << endl;

  sortSlicesInTime(event_slices);

  sliceMatcher(mu_slices, event_slices);
}
//____________________________________________________________________________
std::map<int, std::vector<CandStripHandle *> >
                             AltAlgSliceList::getMuSpecSlices(
                                                    CandStripHandleItr cshItr)
{
// Finding slices in the muon spectrometer

  //-- Get the time-profile of muon spectrometer activity:
  //   CandStripHandleItr cshItr points only to strips in pl > MaxPlane which
  //   defines the max plane before the muon spectrometer
  //   Use the same time-resolution as in the fTimeResolution

  MSG("AltAlg", Msg::kVerbose)
          << " |--> Getting the muon spectrometer strip time profile" << endl;

  fSubsetTimeProfile = createSubsetTimeProfile(&cshItr);

  fSubsettmin = ftmin;
  fSubsettmax = ftmax;

  //-- Use the peak-finder to determine time slices in the profile of the muon
  //   spectrometer activity

  MSG("AltAlg", Msg::kVerbose)
                  << " |--> Finding time slices in muon spectrometer" << endl;

  std::vector<TimeSlice_t> mu_seeds = getSliceSeeds(
                                             &cshItr, false, kMuSpectrometer);

  std::map<int, std::vector<CandStripHandle *> >
                                 mu_slices = fillSliceSeeds(cshItr, mu_seeds);

  MSG("AltAlg", Msg::kVerbose)
                    << " |--> Muon spectrometer slice seed filtering" << endl;

  //initSliceFiltering(mu_slices);

  sortSlicesInTime(mu_slices);

  delete fSubsetTimeProfile;
  fSubsetTimeProfile = 0;

  return mu_slices;
}
//____________________________________________________________________________
void AltAlgSliceList::sliceMatcher(
          std::map<int, std::vector<CandStripHandle *> > & mu_slices,
          std::map<int, std::vector<CandStripHandle *> > & event_slices) const

{
// Match muon spectrometer slices with upstream detector slices

  // Find existing slices with 'activity' just before the muon spectrometer

  MSG("AltAlg", Msg::kVerbose)
      << " |--[-] Finding existing slices with activity just upstream of the "
      << "muon spectrometer" << endl;

  std::vector<int> slice_ids = slicesActiveUpstreamOfMuSpec(event_slices);

  // Examine all muon spectrometer slices for matches with the above slices.

  std::map<int, std::vector<CandStripHandle *> >::iterator mu_slice_iter;

  std::vector<CandStripHandle *> rubbish;

  for(mu_slice_iter = mu_slices.begin();
                         mu_slice_iter != mu_slices.end(); ++mu_slice_iter) {

     std::vector<CandStripHandle *> & mu_slice = (*mu_slice_iter).second;

     MSG("AltAlg", Msg::kVerbose)
         << " |--[-] muon-spectr. slice: " << (*mu_slice_iter).first << endl;

     // Meaningfull to proceed only if muon spectrometer slice has hit strips
     // in the upstream part of the muon spectrometer

     if( sliceHasStripsInUpstreamMuSpectrometer(mu_slice) ) {

        // Check for match... If a match is found the strips will have already
        // be appended to the right slice... if no match is found take care of
        // the unmatched strips...

        if( ! findMatchForMuSpecSlice(mu_slice, slice_ids, event_slices) ) {

          if(fMuSpecSuppressUnmatchedSlices) {

             // do not create a separate slice for every unmatched mu-spec slc.
             MSG("AltAlg", Msg::kVerbose)
                << "     |--> NO match - adding strips to rubbish slc" << endl;
             addMuSpecStripsToExistingSlice(mu_slice, rubbish);

          } else {
             // do create a separate slice for every unmatched mu-spec slc.
             MSG("AltAlg", Msg::kVerbose)
                          << "     |--> NO match - creating new slice" << endl;
             newSliceToHostMuSpecStrips(mu_slice, event_slices);

          } // if suppress creation of unmatched slices
       } // if no match was found
     } // if mu spectrometer slice has activity in its upstream section
  }// loop over muon spectrometer slices

  // If it chosen to suppress the creation of separate slices for every single
  // unmatched mu-spectrometer slice, the rubbish strips are added as a single
  // (rubbish) slice...

  if(fMuSpecSuppressUnmatchedSlices)
                            newSliceToHostMuSpecStrips(rubbish, event_slices);
}
//____________________________________________________________________________
std::vector<int> AltAlgSliceList::slicesActiveUpstreamOfMuSpec(
    const std::map<int, std::vector<CandStripHandle *> > & event_slices) const
{
// Find out which upstream slices have activity in the part just upstream of
// the muon spectrometer (in the most downstream part of the forward detector)

  std::ostringstream sstream;
  sstream   << "     |--[-] Using AlgConf params: 'activity' is > "
            << fMuSpecNHitStripsBefSpectr  << " hit strips in the "
            << fMuSpecNPlnBefSpectr        << " planes before the mu-spectr.";

  MSG("AltAlg", Msg::kDebug) << sstream.str().c_str() << endl;

  std::vector<int> slice_ids;

  int min_plane_inclusive = kMuSpec1stPlane - fMuSpecNPlnBefSpectr;
  int max_plane_inclusive = kMuSpec1stPlane - 1;

  std::map<int, std::vector<CandStripHandle *> >::const_iterator slice_iter;
  for(slice_iter = event_slices.begin();
                           slice_iter != event_slices.end(); ++slice_iter) {

     int nhits = count_if(
                    (*slice_iter).second.begin(), (*slice_iter).second.end(),
             is_in_z_window_noref(min_plane_inclusive, max_plane_inclusive));

     MSG("AltAlg", Msg::kVerbose)
           << "         |--> slice " << (*slice_iter).first<< " has "<< nhits
           << " hit strips in planes [" << min_plane_inclusive << ", "
           << max_plane_inclusive << "]" << endl;

     if(nhits > fMuSpecNHitStripsBefSpectr) {

         slice_ids.push_back( (*slice_iter).first );

         MSG("AltAlg", Msg::kVerbose)
            << "         |--> *** use this one to when checking for matches "
            << " with mu-spectrometer activity " << endl;
     }
  }  // slice_iter

  return slice_ids;
}
//____________________________________________________________________________
bool AltAlgSliceList::findMatchForMuSpecSlice(
  std::vector<CandStripHandle *> mu_strips, std::vector<int> slice_ids,
         std::map<int, std::vector<CandStripHandle *> > & event_slices) const
{
  bool match_found = false;

  // Since we are at this point, it is likely to have a match with an upstream
  // spectrometer slice. Find the closest one...
  // The closest is decided in terms of temporal distance between the
  // most downstream section of the upstream slice and the
  // most upstream section of the downstream slice...

  std::vector<int>::iterator slice_id_iter;

  std::map<double, int, less<double> > closest_slices;

  double tmu  = averageTimeAtBeginOfDownstreamSlice(mu_strips);

  MSG("AltAlg", Msg::kVerbose)
        << "     |--> t (mu. spec begin) = " << fFmtTime(tmu) << endl;
  MSG("AltAlg", Msg::kVerbose)
        << "     |--[-] checking for matches " << endl;

  for(slice_id_iter = slice_ids.begin();
                         slice_id_iter != slice_ids.end(); ++slice_id_iter) {

      double tslc = averageTimeAtEndOfUpstreamSlice(
                                                event_slices[*slice_id_iter] );

      MSG("AltAlg", Msg::kVerbose)
         << "         |--> tend(slc = " << fFmtSlc(*slice_id_iter) << ") = "
         << fFmtTime(tslc) << " -- dtbeg(mu.slc - slc) = "
         << fFmtTime(tmu-tslc) << endl;

      if(tmu  - tslc < fMuSpecTimeAftUpstrActivity &&
         tslc - tmu  < fMuSpecTimeBefUpstrActivity)
               closest_slices.insert(std::map<double, int>::value_type(
                                             fabs(tslc-tmu), *slice_id_iter));
  }

  if(closest_slices.size() >= 1) {

      // if at least one good candidate was found append the strips to the
      // closest one - this can be modified to use tpos info...

     std::map<double, int, less<double> >::iterator closest;

     closest = closest_slices.begin();

     MSG("AltAlg", Msg::kVerbose)
           << "     |--> mu-slice is appended to slice " << (*closest).second
           << " -- the closest out of " << closest_slices.size() << " slices"
           << endl;

      match_found = true;
      addMuSpecStripsToExistingSlice(
                                  mu_strips, event_slices[(*closest).second]);
  }

  return match_found;
}
//____________________________________________________________________________
bool AltAlgSliceList::sliceHasStripsInUpstreamMuSpectrometer(
                      const std::vector<CandStripHandle *> & mu_strips) const
{
// Check if the muon spectrometer slice has strips in the first few muon
// spectrometer planes

  int nhits = count_if( mu_strips.begin(), mu_strips.end(),
                  is_in_z_window_noref( kMuSpec1stPlane,
                                        kMuSpec1stPlane+fMuSpecNUpstrPlanes));
  MSG("AltAlg", Msg::kVerbose)
           << "     |--> slice has " << nhits << " strips (out of "
           << mu_strips.size() << " in total) in the forward "
           << " mu-spectr. region: planes = [" << kMuSpec1stPlane << ", "
           << kMuSpec1stPlane+fMuSpecNUpstrPlanes << "]" << endl;

  if(nhits == 0) {
      MSG("AltAlg", Msg::kVerbose) << "     |--> *** skipping" << endl;
      return false;
  }

  return true;
}
//____________________________________________________________________________
void AltAlgSliceList::addMuSpecStripsToExistingSlice(
                             const std::vector<CandStripHandle *> & mu_strips,
                          std::vector<CandStripHandle *> & slice_strips) const
{
  std::vector<CandStripHandle *>::const_iterator mu_strip_iter;
  for(mu_strip_iter = mu_strips.begin();
           mu_strip_iter != mu_strips.end(); ++mu_strip_iter)
                                     slice_strips.push_back( *mu_strip_iter );

  sort(slice_strips.begin(), slice_strips.end(), min_t());
}
//____________________________________________________________________________
void AltAlgSliceList::newSliceToHostMuSpecStrips(
                             const std::vector<CandStripHandle *> & mu_strips,
          std::map<int, std::vector<CandStripHandle *> > & event_slices) const
{
  int slice_id = nextAvailableSliceId(event_slices);

  event_slices.insert(
       std::map<int, std::vector<CandStripHandle *> >::value_type(
                                                       slice_id, mu_strips) );

  sort(event_slices[slice_id].begin(), event_slices[slice_id].end(), min_t());
}
//____________________________________________________________________________
TH1D * AltAlgSliceList::createSubsetTimeProfile(CandStripHandleItr * cshItr)
{
// Creates the time-profile for a subset of hit strips.
// The selection of the strip subset (i.e. muon spectrimeter hits only etc)
// is *encapsulated* in the CandStripHandleItr (a proper selection on its
// NavSet using NavTools must have preceded the call to this function).
// The time profile histogram is instantiated by cloning and reseting the
// spill time-profile. This guarantees that the time-profile of all subsets
// have the same bin size (corresponding to the same time resolution)
// The returned object is owned by the calling function.
//
  TH1D * subset_time_profile = (TH1D *) fTimeProfile->Clone();
  subset_time_profile->Reset("ICE"); // "ICE"--check that I remember correctly

  cshItr->ResetFirst();
  while( CandStripHandle * striph = cshItr->Ptr() ) {

     (fPkfWeightProfileWithCharge) ?
          subset_time_profile->Fill( striph->GetCorrBegTime(),
                                     striph->GetCharge()  )  :
          subset_time_profile->Fill( striph->GetCorrBegTime() )  ;

     cshItr->Next();
  }

  return subset_time_profile;
}
//____________________________________________________________________________
double AltAlgSliceList::averageTimeAtBeginOfDownstreamSlice(
                            std::vector<CandStripHandle *> & mu_strips) const
{
  std::vector<CandStripHandle *>::iterator mu_end = partition(
                  mu_strips.begin(), mu_strips.end(),
                  is_in_z_window_noref( kMuSpec1stPlane,
                                        kMuSpec1stPlane+fMuSpecNUpstrPlanes));

  return averageTime( mu_strips, mu_end );
}
//____________________________________________________________________________
double AltAlgSliceList::averageTimeAtEndOfUpstreamSlice(
                            std::vector<CandStripHandle *> & slc_strips) const
{
  std::vector<CandStripHandle *>::iterator slc_end = partition(
                 slc_strips.begin(), slc_strips.end(),
                 is_in_z_window_noref( kMuSpec1stPlane - fMuSpecNPlnBefSpectr,
                                       kMuSpec1stPlane - 1));

  return averageTime(slc_strips, slc_end );
}
//____________________________________________________________________________
double AltAlgSliceList::averageTime(
                     const std::vector<CandStripHandle *> & strips,
                     std::vector<CandStripHandle *>::const_iterator end) const
{
  double qt = accumulate( strips.begin(), end, 0.0, sum_qtime() );
  double q  = accumulate( strips.begin(), end, 0.0, sum_q()     );

  //eventually assert(q!=0)

  if(q==0) {
      MSG("AltAlg", Msg::kWarning)
           << "Charge = " << q << " -- This should not be happenning" << endl;
      return 0;
  } else  return qt/q;
}
//____________________________________________________________________________
double AltAlgSliceList::getStripTime(
         CandStripHandle * csh, StripTime_t st, StripEnd::StripEnd_t se) const
{
  switch ( st ) {
  case ( eTime )        : return csh->GetTime(se);       break;
  case ( eBegTime )     : return csh->GetBegTime(se);    break;
  case ( eEndTime )     : return csh->GetEndTime(se);    break;
  case ( eCorrBegTime ) : return csh->GetCorrBegTime();  break;
  default               : return 0;
  }
}
//____________________________________________________________________________
void AltAlgSliceList::getAlgorithmConfiguration(const AlgConfig & ac)
{
  internalInit(); // initialize internal data members

  fGrfxDebugGraphics              =  ac.GetInt("Grfx_DebugGraphics");
  fGrfxTimeProfileLogView         =  ac.GetInt("Grfx_TimeProfileLogView");  
  fPkfPeakThreshold               =  ac.GetInt("Pkf_PeakThreshold");  
  fPkfNSuccessiveEmptyBins        =  ac.GetInt("Pkf_NSuccessiveEmptyBins");     
  fPkfMuSpecPeakThreshold         =  ac.GetInt("Pkf_MuSpecPeakThreshold");    
  fPkfMuSpecNSuccessiveEmptyBins  =  ac.GetInt("Pkf_MuSpecNSuccessiveEmptyBins");   
  fPkfLowPeakThreshold            =  ac.GetInt("Pkf_LowPeakThreshold");  
  fPkfLowNSuccessiveEmptyBins     =  ac.GetInt("Pkf_LowNSuccessiveEmptyBins");   
  fPkfNofMergedTimeBins           =  ac.GetInt("Pkf_NofMergedTimeBins");  
  fPkfWeightProfileWithCharge     = (ac.GetInt("Pkf_WeightProfileWithCharge") == 1);    
  fPkfRecursivePeakSearch         = (ac.GetInt("Pkf_RecursivePeakSearch") == 1);
  fPkfTimeWindowBefPeak           =  ac.GetInt("Pkf_TimeWindowBefPeak");
  fPkfTimeWindowAftPeak           =  ac.GetInt("Pkf_TimeWindowAftPeak");
  fRefinementDissolving           = (ac.GetInt("Refinement_Dissolving") == 1);
  fRefinementMerging              = (ac.GetInt("Refinement_Merging") == 1);
  fRefinementKMeansClustering     = (ac.GetInt("Refinement_KMeansClustering") == 1);
  fRefinementMSTClustering        = (ac.GetInt("Refinement_MSTClustering") == 1);
  fOrphanStripsQWeight            = (ac.GetInt("OrphanStrips_QWeight") == 1);
  fOrphanStripsPlaneWindow        =  ac.GetInt("OrphanStrips_PlaneWindow");
  fOrphanStripsTimeWindow         =  ac.GetDouble("OrphanStrips_TimeWindow");    
  fKMeansTimeWindow               =  ac.GetDouble("KMeans_TimeWindow");
  fKMeansPlaneWindow              =  ac.GetInt("KMeans_PlaneWindow");
  fKMeansTPosWindow               =  ac.GetDouble("KMeans_TPosWindow");  
  fMuSpecNUpstrPlanes             =  ac.GetInt("MuSpec_NUpstrPlanes");
  fMuSpecTimeAftUpstrActivity     =  ac.GetDouble("MuSpec_TimeAftUpstrActivity");
  fMuSpecTimeBefUpstrActivity     =  ac.GetDouble("MuSpec_TimeBefUpstrActivity");
  fMuSpecNHitStripsBefSpectr      =  ac.GetInt("MuSpec_HitStripsBefSpectr");
  fMuSpecNPlnBefSpectr            =  ac.GetInt("MuSpec_NPlnBefSpectr");
  fMuSpecSuppressUnmatchedSlices  = (ac.GetInt("MuSpec_SuppressUnmatchedSlices") == 1);
  fTimeResolution                 =  ac.GetDouble("TimeResolution");
  fMinNoHitStrips                 =  ac.GetInt("MinHitStripsInSlice");  
  fTimeDiffBetweenPeaks           =  ac.GetDouble("TimeDiffBetweenPeaks");  
  fZDiffBetweenPeaks              =  ac.GetDouble("ZDiffBetweenPeaks");  
  fZDiffBetweenEnds               =  ac.GetDouble("ZDiffBetweenEnds");  
  fUVDiffBetweenPeaks             =  ac.GetDouble("UVDiffBetweenPeaks");
  fMinCharge                      =  ac.GetDouble("MinChargeInSlice");  
}
//____________________________________________________________________________
void AltAlgSliceList::internalInit(void)
{
// Initialize private data members

  //--  ordinary private data members
  
  fNTimeBins                           = 0;
  fUpdateCentroids                     = false;         
  fTimeProfileMax                      = 0;
  fTimeProfile                         = 0;
  ftmin                                = 0;
  ftmax                                = 0;
  fSubsetTimeProfile                   = 0;
  fSubsettmin                          = 0;
  fSubsettmax                          = 0;
  fPeakFinderNestingLevel              = 0;
  fkMeansIteration                     = 0;

  //-- private data members that are set by external AlgConf object

  fGrfxDebugGraphics                   = 0;
  fGrfxTimeProfileLogView              = 0;
  fPkfPeakThreshold                    = 0;
  fPkfNSuccessiveEmptyBins             = 0;
  fPkfMuSpecPeakThreshold              = 0;
  fPkfMuSpecNSuccessiveEmptyBins       = 0;
  fPkfLowPeakThreshold                 = 0;
  fPkfLowNSuccessiveEmptyBins          = 0;
  fPkfNofMergedTimeBins                = 0;   
  fPkfTimeWindowBefPeak                = 0;  
  fPkfTimeWindowAftPeak                = 0;  
  fOrphanStripsPlaneWindow             = 0;
  fKMeansPlaneWindow                   = 0;
  fMuSpecNUpstrPlanes                  = 0;
  fMuSpecNHitStripsBefSpectr           = 0;
  fMuSpecNPlnBefSpectr                 = 0;
  fMinNoHitStrips                      = 0;            
  fPkfWeightProfileWithCharge          = false;
  fPkfRecursivePeakSearch              = false;
  fRefinementDissolving                = false;
  fRefinementMerging                   = false;      
  fRefinementKMeansClustering          = false;
  fRefinementMSTClustering             = false;
  fMuSpecSuppressUnmatchedSlices       = false;
  fOrphanStripsQWeight                 = false;   
  fOrphanStripsTimeWindow              = 0.0; 
  fKMeansTimeWindow                    = 0.0;
  fKMeansTPosWindow                    = 0.0;      
  fMuSpecTimeAftUpstrActivity          = 0.0;
  fMuSpecTimeBefUpstrActivity          = 0.0;
  fTimeResolution                      = 0.0;
  fTimeDiffBetweenPeaks                = 0.0;
  fZDiffBetweenPeaks                   = 0.0;
  fZDiffBetweenEnds                    = 0.0;
  fUVDiffBetweenPeaks                  = 0.0;  
  fMinCharge                           = 0.0;         
}
//____________________________________________________________________________
void AltAlgSliceList::buildCandidates(
                std::map<int, std::vector<CandStripHandle *> > & event_slices,
                                            CandHandle & ch, CandContext & cx) 
{
  sortSlicesInTime(event_slices);

  //-- Get an AlgFactory - ask it for a handle to AltAlgSlice algorithm
  //   and create a candidate context to supply the algorithm with

  AlgFactory & af = AlgFactory::GetInstance();
  AlgHandle    ah = af.GetAlgHandle("AltAlgSlice","default");
  
  CandContext ccx(this, cx.GetMom());

  std::map<int, std::vector<CandStripHandle *> >::const_iterator slc_iter;
  
  for( slc_iter = event_slices.begin(); 
                                 slc_iter != event_slices.end(); ++slc_iter) {
 
     MSG("AltAlg", Msg::kDebug) 
        << " Slice: " << fFmtSlc( (*slc_iter).first )
        << " / Running MakeCandidate on AltAlgSlice & Adding Doughter Link" 
        << endl;

     AltWrapperStlVecStripHandle vec_wrapper( (*slc_iter).second );
     ccx.SetDataIn( & vec_wrapper); 
     CandSliceHandle slchandle = CandSlice::MakeCandidate(ah,ccx);
     ch.AddDaughterLink(slchandle);     
  }
}
//____________________________________________________________________________

---