AlgSubShowerSRList.cxx

Go to the documentation of this file.

// $Id: AlgSubShowerSRList.cxx,v 1.32 2009/02/13 15:53:37 asousa Exp $
//
// AlgSubShowerSRList.cxx
//
// AlgSubShowerSRList is a concrete SubShowerSRList Algorithm class.
//

#include <cassert>
#include <fstream>
#include <vector>
#include "Algorithm/AlgFactory.h"
#include "Algorithm/AlgHandle.h"
#include "Algorithm/AlgConfig.h"
#include "Candidate/CandContext.h"
#include "CandSubShowerSR/AlgSubShowerSR.h"
#include "CandSubShowerSR/AlgSubShowerSRList.h"
#include "CandSubShowerSR/CandSubShowerSR.h"
#include "CandSubShowerSR/CandSubShowerSRHandle.h"
#include "CandSubShowerSR/CandSubShowerSRList.h"
#include "CandSubShowerSR/CandSubShowerSRListHandle.h"
#include "Conventions/PlaneView.h"
#include "CandSubShowerSR/ClusterType.h"
#include "MessageService/MsgService.h"
#include "MinosObjectMap/MomNavigator.h"
#include "Navigation/NavKey.h"
#include "Navigation/NavSet.h"
#include "RecoBase/CandSliceHandle.h"
#include "RecoBase/CandTrackHandle.h"
#include "RecoBase/CandClusterHandle.h"
#include "RecoBase/CandSliceListHandle.h"
#include "RecoBase/CandTrackListHandle.h"
#include "RecoBase/CandClusterListHandle.h"
#include "Validity/VldContext.h"
#include "UgliGeometry/UgliGeomHandle.h"
#include "Calibrator/Calibrator.h"
#include "UgliGeometry/UgliGeomHandle.h"
#include "TSpectrum.h"
#include "TMath.h"
#include "TGraphErrors.h"
#include "TCanvas.h"
#include "TStyle.h"

#define MAXHOUGHITER 5
#define STRIPWIDTHINMETRES 0.041
#define PITCHINMETRES 0.06
#define FIRSTNDSPECPLANE 121

ClassImp(AlgSubShowerSRList)

//______________________________________________________________________
CVSID("$Id: AlgSubShowerSRList.cxx,v 1.32 2009/02/13 15:53:37 asousa Exp $");
//______________________________________________________________________

AlgSubShowerSRList::AlgSubShowerSRList() :
  fLongWindowMipCut(0.0),fLongWindowGapPlaneCut(4),
  fMinHoughPlanes(3),fMinHoughPH(1.0), fBestHoughMaxPHFrac(0.75),
  fCleanUpTimeCut(50.0),fCleanUpPlaneCut(2),fCleanUpStripCut(2.5),
  fFormHalo(1),fHaloMaxPlaneGap(2),fHaloNPeakFindingBins(64),
  fMinStripPE(2)
{
  vtxHist = new TH1F("vtxHist","vtxHist",1,-4.,4.);
  houghHist = new TH2F("houghHist","Hough Histogram",1,-1,1,1,-1,1);
  phHoughHist = new TH2F("phHoughHist","PH Weighted Hough Histogram",1,-1,1,1,-1,1);
  ssPol1 = new TF1("ssPol1","pol1",0,500);
}

//______________________________________________________________________
AlgSubShowerSRList::~AlgSubShowerSRList()
{
  delete vtxHist;
  delete houghHist;
  delete phHoughHist;
  delete ssPol1;
}

//______________________________________________________________________
void AlgSubShowerSRList::RunAlg(AlgConfig &ac, CandHandle &ch,CandContext &cx)
{

  MSG("SubShowerSR", Msg::kDebug) << "Starting AlgSubShowerSRList::RunAlg()" 
                                  << endl;
  assert(cx.GetDataIn());
  if (!(cx.GetDataIn()->InheritsFrom("TObjArray"))) {
    return;
  }
  
  const CandSliceListHandle *slicelist = 0;
  const CandClusterListHandle *clusterlist = 0;
  const CandTrackListHandle *tracklist = 0;
  const TObjArray *cxin = dynamic_cast<const TObjArray *>(cx.GetDataIn());
  for (Int_t i=0; i<=cxin->GetLast(); i++) {
    TObject *tobj = cxin->At(i);
    MSG("SubShowerSR", Msg::kDebug) << "cxin->At(" << i << ")"<< endl;
    if (tobj->InheritsFrom("CandSliceListHandle")) {
      slicelist = dynamic_cast<CandSliceListHandle*>(tobj);
      MSG("SubShowerSR", Msg::kDebug) << "Got SliceList" << endl;
    }
    if (tobj->InheritsFrom("CandClusterListHandle")) {
      clusterlist = dynamic_cast<CandClusterListHandle*>(tobj);
      MSG("SubShowerSR", Msg::kDebug) << "Got ClusterList" << endl;
    }
    if (tobj->InheritsFrom("CandTrackListHandle")) {
      tracklist = dynamic_cast<CandTrackListHandle*>(tobj);
      MSG("SubShowerSR", Msg::kDebug) << "Got TrackList" << endl;
    }
  }
  if (!slicelist || !clusterlist || !tracklist) {
    MSG("SubShowerSR",Msg::kInfo) << "CandSliceListHandle or " 
                                  << "CandClusterListHandle or " 
                                  << "CandTrackListHandle missing\n";
    return;
  }

  MSG("SubShowerSR", Msg::kDebug) << "Creating CandContext" << endl;
  // Create Candcontext
  CandContext cxx(this,cx.GetMom());

  //get config for AlgSubShowerSRList
  const char *charSubShowerSRAlgConfig = 0;
  ac.Get("SubShowerSRAlgConfig",charSubShowerSRAlgConfig);

  Int_t minNTrkOnlyPlanesToUseTrackStrips = ac.GetInt("MinNTrkOnlyPlanesToUseTrackStrips");
  Int_t maxPlanesFromVtxToUseTrackStrips  = ac.GetInt("MaxPlanesFromVertexToUseTrackStrips");
  Double_t minMIPsToUseTrackStrips        = ac.GetDouble("MinMIPsToUseTrackStrips");

  fLongWindowMipCut      = ac.GetDouble("LongWindowMipCut");
  fLongWindowGapPlaneCut = ac.GetInt("LongWindowGapPlaneCut");
  fMinHoughPlanes        = ac.GetInt("MinHoughPlanes");
  fMinHoughPH            = ac.GetDouble("MinHoughPH");
  fBestHoughMaxPHFrac    = ac.GetDouble("BestHoughMaxPHFrac");
  fCleanUpTimeCut        = ac.GetDouble("CleanUpTimeCut");
  fCleanUpPlaneCut       = ac.GetInt("CleanUpPlaneCut");
  fCleanUpStripCut       = ac.GetDouble("CleanUpStripCut");
  fFormHalo              = ac.GetInt("FormHalo");
  fHaloMaxPlaneGap       = ac.GetInt("HaloMaxPlaneGap");
  fHaloNPeakFindingBins  = ac.GetInt("HaloNPeakFindingBins");
  fMinStripPE            = ac.GetDouble("MinStripPE");
  
  //Get singleton instance of AlgFactory
  AlgFactory &af = AlgFactory::GetInstance();
  AlgHandle ah = af.GetAlgHandle("AlgSubShowerSR",charSubShowerSRAlgConfig);
                                                                   
  const CandRecord *candrec = cx.GetCandRecord();
  assert(candrec);
  const VldContext *vldcptr = candrec->GetVldContext();
  assert(vldcptr);
  VldContext vldc = *vldcptr;                                         
  UgliGeomHandle ugh(vldc);
  Int_t detector = 0;
  if(vldc.GetDetector()==Detector::kFar) detector = 1;

  Calibrator& calibrator = Calibrator::Instance();

  //keep track of all unassigned strips in all slices and try to add them later
  TObjArray unassignedStrips;

  //loop over all slices
  Int_t sliceCounter = 0;
  CandSliceHandleItr sliceItr(slicelist->GetDaughterIterator());
  while (CandSliceHandle *slice = sliceItr()) {

    MSG("SubShowerSR",Msg::kVerbose) << "Slice " <<  sliceCounter << endl;

    //also keep track of cluster planes and track planes
    Int_t cluPln[500] = {0};
    Int_t trkPln[500] = {0};
    Int_t nTrkOnlyPlanes = 0;
    if(tracklist!=0) {
      //loop through all tracks in track list
      CandTrackHandleItr trackItr(tracklist->GetDaughterIterator());
      while (CandTrackHandle *track = trackItr()){
        //if this track is in the current slice then...
        if (*track->GetCandSlice()==*slice) {
          for(Int_t i=track->GetBegPlane();i<=track->GetEndPlane();i++){
            trkPln[i] += 1;   //if trkPln > 1 => multiple tracks in plane
            nTrkOnlyPlanes += 1;
          }
        }
      }
    }
    if(clusterlist!=0) {
      //loop through all clusters in cluster list
      CandClusterHandleItr clusterItr(clusterlist->GetDaughterIterator());
      while (CandClusterHandle *cluster = clusterItr()){
        //if this cluster is in the current slice then...
        if (*cluster->GetCandSlice()==*slice) {
          for(Int_t i=cluster->GetBegPlane();i<=cluster->GetEndPlane();i+=2){
            cluPln[i] += 1;   //if cluPln > 1 => multiple clusters in plane
            if(trkPln[i]>0) nTrkOnlyPlanes-=1;
            if(detector==1){//account for SM gap
              if(i==247) i = 250; 
              if(i==248) i = 251;
            }
          }
        }
      }
    }

    //for each slice, make an Object Array to hold all the 
    //non track strips in the slice
    TObjArray allUStrips;
    TObjArray allVStrips;
    
    //loop through the strips in the slice
    CandStripHandleItr stripslcItr(slice->GetDaughterIterator());    
    MSG("SubShowerSR",Msg::kVerbose) << "About to loop through slice strips" 
                                     << endl;
    
    while (CandStripHandle *stp =stripslcItr()) { //while over slice stp

      // Only consider strips above the required minimum
      if (stp->GetCharge(CalDigitType::kPE) < fMinStripPE) continue;

      bool gotStrip = false;
      MSG("SubShowerSR",Msg::kVerbose) << "PlaneView " 
                                       << stp->GetPlaneView() << endl;
      
      CandClusterHandleItr clusterItr(clusterlist->GetDaughterIterator());
      Int_t clusterCounter = 0;
      while(CandClusterHandle *cluster = clusterItr()){ 
        MSG("SubShowerSR",Msg::kVerbose) 
          << "Looking for stp match in cluster " << clusterCounter << endl;
        if(*cluster->GetCandSlice()==*slice){
          CandStripHandleItr stripcluItr(cluster->GetDaughterIterator());
          while (CandStripHandle *stpclu =stripcluItr()){
            if(*stp==*stpclu){
              gotStrip = true;
              //add strip to object array then break out of loop
              if(stp->GetPlaneView()==PlaneView::kU ||
                 stp->GetPlaneView()==PlaneView::kX) {
                allUStrips.Add(stp);
                MSG("SubShowerSR",Msg::kVerbose) << "Found match! U strip" 
                                                 << endl;             
              }
              else if(stp->GetPlaneView()==PlaneView::kV ||
                      stp->GetPlaneView()==PlaneView::kY) {
                allVStrips.Add(stp);
                MSG("SubShowerSR",Msg::kVerbose) << "Found match! V strip" 
                                                 << endl;
              }
              else MSG("SubShowerSR",Msg::kError) 
                << "Unknown PlaneView! Not adding strip to initial ObjArray" 
                << endl;
              break;
            }
          }
        }
        if(gotStrip) break;
        clusterCounter++;
      }
      
      //if this is an ND strip in spectrometer, don't add to array
      if(!gotStrip && detector==0 && stp->GetPlane()>=FIRSTNDSPECPLANE) gotStrip=true;

      if(!gotStrip){  //if !gotStrip
        bool inTrack = false;
        Int_t nPlanesFromBegPlane = 0;
        if(tracklist!=0) {
          MSG("SubShowerSR",Msg::kVerbose) 
            << "No matched strip found in clusters, looking in tracks"
            << endl;
          
          //loop through all tracks in track list
          CandTrackHandleItr trackItr(tracklist->GetDaughterIterator());
          Int_t trackCounter = 0;
          while (CandTrackHandle *track = trackItr()){
            MSG("SubShowerSR",Msg::kVerbose) 
              << "Looking for match of slc stp with trk stp in track "
              << trackCounter << endl;
            //if this track is in the current slice then...
            if (*track->GetCandSlice()==*slice) {
              //...loop over all strips in this track
              CandStripHandleItr striptrkItr(track->GetDaughterIterator());
              while (CandStripHandle *stptrk =striptrkItr()){
                //if this strip is in this track then...
                if (*stp == *stptrk) {
                  MSG("SubShowerSR",Msg::kVerbose) 
                    << "Found match! This strip is in a track!" << endl;
                  //...break and move on to the next strip in the slice
                  inTrack = true;
                  nPlanesFromBegPlane = TMath::Abs(stp->GetPlane()-track->GetBegPlane());
                  break;
                }
              }
            }
            if(inTrack) break;
            trackCounter++;
          }
        }
        
        MSG("SubShowerSR",Msg::kVerbose) << "nTrkOnlyPlanes = " 
                                         << nTrkOnlyPlanes << endl;
        if(inTrack && nTrkOnlyPlanes<minNTrkOnlyPlanesToUseTrackStrips) {
          MSG("SubShowerSR",Msg::kDebug) 
            << "nTrkOnlyPlanes<" << minNTrkOnlyPlanesToUseTrackStrips 
            << " - going to include this strip " 
            << "even though it's in a track." << endl;
          inTrack = false; //ignore tracks that are buried in showers
        }
        else if(inTrack && nPlanesFromBegPlane <
                maxPlanesFromVtxToUseTrackStrips) { //near vertex
          if(calibrator.GetMIP(stp->GetCharge(CalDigitType::kSigCorr)) >
             minMIPsToUseTrackStrips) {
            MSG("SubShowerSR",Msg::kDebug) 
              << "nPlanesFromBegPlane < "
              << maxPlanesFromVtxToUseTrackStrips 
              << " for this strip and ph = "
              << calibrator.GetMIP(stp->GetCharge(CalDigitType::kSigCorr))
              << " mips. "
              << "Going to include it even though it's in a track." 
              << endl;
            inTrack = false;
          }
        }

        //then: add it to our array
        if(!inTrack){
          MSG("SubShowerSR",Msg::kVerbose) << "Strip not found in track! ";
          if(stp->GetPlaneView()==PlaneView::kU ||
             stp->GetPlaneView()==PlaneView::kX) {
            allUStrips.Add(stp);
            MSG("SubShowerSR",Msg::kVerbose) << "Adding U strip" << endl;
          }
          else if(stp->GetPlaneView()==PlaneView::kV ||
                  stp->GetPlaneView()==PlaneView::kY) {
            allVStrips.Add(stp);
            MSG("SubShowerSR",Msg::kVerbose) << "Adding V strip" << endl;
          }
          else MSG("SubShowerSR",Msg::kError) 
            << "Unknown PlaneView! Not adding strip to initial ObjArray" 
            << endl;
        }
      }
    }

    MSG("SubShowerSR",Msg::kDebug) << "Found total of " 
                                   << allUStrips.GetLast()+1 << " U strips" 
                                   << " and " << allVStrips.GetLast()+1 
                                   << " V strips" << endl;
    
    //In here put code to pull out the strips belonging to the SubShowers:
    bool keepGoingU = true;
    bool keepGoingV = true;
    while(keepGoingU || keepGoingV){
      if(keepGoingU) 
        MSG("SubShowerSR",Msg::kDebug) << "Still have " 
                                       << allUStrips.GetLast()+1
                                       << " U strips to process" << endl;
      
      else if(keepGoingV) 
        MSG("SubShowerSR",Msg::kDebug) << "Still have " 
                                       << allVStrips.GetLast()+1
                                       << " V strips to process" << endl;

      if(allUStrips.GetLast()+1<=0) keepGoingU = false;
      if(allVStrips.GetLast()+1<=0) keepGoingV = false;

      TObjArray newSubShower;
      if(keepGoingU){ //do U first
        this->FindCluster(&allUStrips,&newSubShower,detector);
        this->TimeTest(&newSubShower);
        while(this->TestOverLap(&newSubShower,ch,slice));
        MSG("SubShowerSR",Msg::kDebug) << "Added " << newSubShower.GetLast()+1 
                                       << " U strips to newSubShower" << endl;
      }
      else if(keepGoingV) { //then V strips
        TObjArray protoSubShower;
        this->FindCluster(&allVStrips,&newSubShower,detector);
        this->TimeTest(&newSubShower);
        while(this->TestOverLap(&newSubShower,ch,slice));
        MSG("SubShowerSR",Msg::kDebug) << "Added " << newSubShower.GetLast()+1 
                                       << " V strips to newSubShower" << endl;
     }

      //Everytime you have a newSubShower with all the strips 
      //you want to add do this
      cxx.SetDataIn(&newSubShower);
      MSG("SubShowerSR",Msg::kDebug) << "forming SubShowerSR\n";
      CandSubShowerSRHandle subshowersrhandle = 
        CandSubShowerSR::MakeCandidate(ah,cxx);
      // Set the strip PE cut - this gets passed through the subshower
      // to the shower, so that it is accessible from the truthhelper
      // when it fills the purity and completeness variables
      subshowersrhandle.SetMinStripPE(fMinStripPE);

      if(subshowersrhandle.GetNStrip()!=0) {
        
        MSG("SubShowerSR",Msg::kDebug) << "Formed SubShowerSR with " 
                                       << subshowersrhandle.GetNStrip() 
                                       << " strips" << endl;
        
        subshowersrhandle.SetCandSlice(slice);
        ch.AddDaughterLink(subshowersrhandle);
        
        CandStripHandleItr stripItr(subshowersrhandle.GetDaughterIterator());

        std::vector<Int_t> tmpUArray(allUStrips.GetLast()+1,0);
        for(int i=0;i<allUStrips.GetLast()+1;i++) tmpUArray[i] = 0;
        std::vector<Int_t> tmpVArray(allVStrips.GetLast()+1,0);
        for(int i=0;i<allVStrips.GetLast()+1;i++) tmpVArray[i] = 0;     
        
        while(CandStripHandle *stp = stripItr()){
          if(stp->GetPlaneView()==PlaneView::kU ||
             stp->GetPlaneView()==PlaneView::kX){
            for (int i=0; i<=allUStrips.GetLast(); i++) {
              TObject *tobj = allUStrips.At(i);
              assert(tobj->InheritsFrom("CandStripHandle"));
              CandStripHandle *striphandle = 
                dynamic_cast<CandStripHandle*>(tobj);
              if(*stp==*striphandle) {
                tmpUArray[i] = 1;
                break;
              }
            }
          }
          else if(stp->GetPlaneView()==PlaneView::kV ||
                  stp->GetPlaneView()==PlaneView::kY){
            for (int i=0; i<=allVStrips.GetLast(); i++) {
              TObject *tobj = allVStrips.At(i);
              assert(tobj->InheritsFrom("CandStripHandle"));
              CandStripHandle *striphandle = 
                dynamic_cast<CandStripHandle*>(tobj);
              if(*stp==*striphandle) {
                tmpVArray[i] = 1;
                break;
              }
            }
          }
        }
        
        //remove used strips from arrays
        for(int i=0;i<allUStrips.GetLast()+1;i++) {
          if(tmpUArray[i]==1) allUStrips.RemoveAt(i);
        }
        allUStrips.Compress();
        if(allUStrips.GetLast()==-1) keepGoingU = false;
        
        for(int i=0;i<allVStrips.GetLast()+1;i++) {
          if(tmpVArray[i]==1) allVStrips.RemoveAt(i);
        }
        allVStrips.Compress();
        if(allVStrips.GetLast()==-1) keepGoingV = false;
      
      }
      else {
        MSG("SubShowerSR",Msg::kDebug) 
          << "SubShowerSR formed with 0 strips - not adding to list" 
          << endl;
        if(keepGoingU) {
          MSG("SubShowerSR",Msg::kDebug) << "Running CleanUp() on U Strips" 
                                         << endl;       
          keepGoingU = CleanUp(&allUStrips,ch,detector,slice);
          if(!keepGoingU) {
            for(int i=0;i<allUStrips.GetLast()+1;i++) {
              unassignedStrips.Add(allUStrips.At(i));
            }
          }
        }
        else {
          MSG("SubShowerSR",Msg::kDebug) << "Running CleanUp() on V Strips" 
                                         << endl;
          keepGoingV = CleanUp(&allVStrips,ch,detector,slice);
          if(!keepGoingV) {
            for(int i=0;i<allVStrips.GetLast()+1;i++) {
              unassignedStrips.Add(allVStrips.At(i));
            }
          }
        }
      }
    }
    sliceCounter++;
  }

  //slice-wise cleanup
  if(fFormHalo&&FormHalo(&unassignedStrips,ch,cxx,ah,slicelist,detector)) 
    MSG("SubShowerSR",Msg::kDebug) << "Halo SubShowers formed" << endl;  

}

//******************************************************************
Bool_t AlgSubShowerSRList::FindCluster(TObjArray *allStrips,
                                       TObjArray *newSubShower,Int_t det)
{

  Int_t nstp = allStrips->GetLast()+1;
  if(nstp==0){
    MSG("SubShower", Msg::kWarning) << "No strips present!" << endl;
    return false;
  }

  //Get useful info from candidates for clustering
  //Read in info:
  std::vector<Int_t> plane(nstp,0);
  std::vector<Int_t> strip(nstp,0);
  std::vector<Double_t> ph(nstp,0.0);
  std::vector<Double_t> pherr(nstp,0.0);
  std::vector<Double_t> z(nstp,0.0);
  std::vector<Double_t> tpos(nstp,0.0);

  //neighbouring strip info
  std::vector<Double_t> TransGradientPlus(nstp,0.0);
  std::vector<Double_t> TransGradientMinus(nstp,0.0);
  std::vector<Double_t> TransGradientErrorPlus(nstp,0.0);
  std::vector<Double_t> TransGradientErrorMinus(nstp,0.0);

  CandStripHandle *usefulStrip = 0;
  int begPlane = 500;
  int endPlane = 0;
  
  Int_t doubleCounter = 0;
  Int_t alreadyGot[500][192] = {{0}};
  Calibrator& calibrator=Calibrator::Instance();
  for (Int_t i=0; i<=allStrips->GetLast(); i++) {
    TObject *tobj = allStrips->At(i);
    assert(tobj->InheritsFrom("CandStripHandle"));
    CandStripHandle *stp = dynamic_cast<CandStripHandle*>(tobj);
    if(alreadyGot[stp->GetPlane()][stp->GetStrip()]==1) {
      for(int j=0;j<i-doubleCounter;j++){
        if(plane[j]==stp->GetPlane() && strip[j]==stp->GetStrip()){
          ph[j] += calibrator.GetMIP(stp->GetCharge(CalDigitType::kSigCorr));
          TransGradientPlus[j] = ph[j];
          TransGradientMinus[j] = ph[j];
          break;
        }
      }
      doubleCounter+=1;
      nstp -= 1;
    }
    else {
      alreadyGot[stp->GetPlane()][stp->GetStrip()] = 1;
      plane[i-doubleCounter] = stp->GetPlane();
      if (plane[i-doubleCounter]<begPlane) begPlane = plane[i-doubleCounter];
      if (plane[i-doubleCounter]>endPlane) endPlane = plane[i-doubleCounter];
      strip[i-doubleCounter] = stp->GetStrip();
      
      ph[i-doubleCounter] = 
        calibrator.GetMIP(stp->GetCharge(CalDigitType::kSigCorr));
      pherr[i-doubleCounter] = ph[i-doubleCounter] /
        TMath::Sqrt(stp->GetCharge(CalDigitType::kPE)); //pe stats err
      
      z[i-doubleCounter] = stp->GetZPos();
      tpos[i-doubleCounter] = stp->GetTPos();
      
      TransGradientPlus[i-doubleCounter] = ph[i-doubleCounter];
      TransGradientErrorPlus[i-doubleCounter] = 
        TMath::Power(pherr[i-doubleCounter],2);
      
      TransGradientMinus[i-doubleCounter] = ph[i-doubleCounter];
      TransGradientErrorMinus[i-doubleCounter] = 
        TMath::Power(pherr[i-doubleCounter],2);
      
      usefulStrip = stp;
    }
  }

  //find transverse neighbour strip index and PH gradient for each strip
  for(Int_t i=0;i<nstp;i++){
    for(Int_t j=0;j<nstp;j++){
      if (plane[j]==plane[i]){
        if(strip[j]==strip[i]-1) {
          TransGradientMinus[i] = TransGradientMinus[i] - ph[j];
          TransGradientErrorMinus[i] = TMath::Sqrt(TransGradientErrorMinus[i]+ 
                                                   TMath::Power(pherr[j],2));
        }
        if(strip[j]==strip[i]+1) {
          TransGradientPlus[i] = TransGradientPlus[i] - ph[j];
          TransGradientErrorPlus[i] = TMath::Sqrt(TransGradientErrorPlus[i]+ 
                                                  TMath::Power(pherr[j],2));
          TransGradientErrorMinus[i] = 0;
          TransGradientErrorPlus[i] = 0;
        }
      }
    }
  }

  //find per plane PH
  int npls = (endPlane-begPlane)/2+1;
  if(det==1&&(begPlane-249)*(endPlane-249)<0) npls = (endPlane-begPlane-1)/2+1;
  std::vector<Int_t> allpl(npls,0);
  std::vector<Double_t> allplph(npls,0.0);
  std::vector<Int_t> allstp(npls,0);

  for(Int_t i=0;i<npls;i++){
    allpl[i] = 2*i+begPlane;
    if(det==1&&(begPlane-249)*(endPlane-249)<0&&(2*i+begPlane)>249)
      allpl[i] = 2*i+begPlane+1;
    allplph[i] = 0.;
    allstp[i] = 0;
    for(Int_t j=0;j<nstp;j++){
      if (plane[j]==allpl[i]){
        allplph[i] += ph[j];
        allstp[i] += 1;
      }
    }
  }
  UgliGeomHandle ugh(*usefulStrip->GetVldContext());

  //Start Clustering: 
  //Find contiguous longitudinal windows
  MSG("SubShowerSR",Msg::kDebug) <<"start longitudinal clustering"<<endl;

  std::vector<int> zeroPlanes;  //vector of planes with PH = 0
  int start = begPlane-2;
  if(det==1&&(begPlane-2-249)*(begPlane-249)<0) start = begPlane-3;
  zeroPlanes.push_back(start-2);
  zeroPlanes.push_back(start);
  for(Int_t i=0;i<npls;i++){
    if (allplph[i]<=fLongWindowMipCut) zeroPlanes.push_back(allpl[i]);
  }
  int over = endPlane+2;
  if(det==1&&(endPlane+2-249)*(endPlane-249)<0) over = endPlane+3;
  zeroPlanes.push_back(over);
  zeroPlanes.push_back(over+2);
  
  MSG("SubShowerSR",Msg::kVerbose) << "Number of Zero Planes = " 
                                   << zeroPlanes.size()
                                   << " First Zero Plane = " << start 
                                   << " Last Zero Plane =  "<< over << endl;
  
  //Find "Zero PH" windows

  //position of first plane in a longitudinal gap window:
  std::vector<int> firstGapPlane;
  //position of last plane in a longitudinal gap window:
  std::vector<int> lastGapPlane;
  int begGap = 0;
  int endGap = -5;
  //filling gap window vectors:
  for (UInt_t i = 0; i<zeroPlanes.size(); i++){
    bool foundgap = false;
    if(zeroPlanes.at(i)>endGap) {
      begGap = zeroPlanes.at(i);
      endGap = zeroPlanes.at(i);
      int j = 0;
      while ((i+j+1)<zeroPlanes.size() && 
             zeroPlanes.at(i+j+1)-zeroPlanes.at(i+j) < 
             fLongWindowGapPlaneCut) {
        foundgap = true;
        endGap = zeroPlanes.at(i+j+1);
        j++;
      }
      if (foundgap){
        firstGapPlane.push_back(begGap);
        lastGapPlane.push_back(endGap);
      }
    }
  }

  //Find highest strip density longitudinal window
  Int_t maxWinBegPlane = 0;  //first boundary plane of max window
  Int_t maxWinEndPlane = 0;
  Float_t maxDensity = 0;
  for (UInt_t i = 0; i<firstGapPlane.size(); i++){
    if((i+1)<firstGapPlane.size()){
      Float_t density = 0;
      Int_t bpln = lastGapPlane.at(i);
      Int_t epln = firstGapPlane.at(i+1);
      for(int j=0;j<npls;j++){
        if(allpl[j]>=bpln && allpl[j]<=epln){
          density += Float_t(allstp[j]);
        }
      }
      density/=Float_t(epln-bpln+2);
      MSG("SubShowerSR",Msg::kDebug) << "Proto-long window [" << bpln
                                     << "," << epln << "] density = " 
                                     << density << endl;
      if(density > maxDensity) {
        maxWinBegPlane = bpln;
        maxWinEndPlane = epln;
        maxDensity = density;
      }
    }
  }

  MSG("SubShowerSR",Msg::kDebug) << "Using Longitudinal Window: "
                                 << maxWinBegPlane << "-"
                                 << maxWinEndPlane << endl;
  
  Int_t maxPlane = -1;
  Double_t maxPlanePH = 0;
  for(int i=0;i<npls;i++){
    if(allpl[i]<maxWinEndPlane && allpl[i]>maxWinBegPlane) {
      if(allplph[i]>maxPlanePH) {
        maxPlane = allpl[i];
        maxPlanePH = allplph[i];
      }
    }
  }
  MSG("SubShowerSR",Msg::kVerbose) 
    << "Maximum PH plane within Long. window = " 
    << maxPlane << endl;

  Int_t numPreMaxPlanes = maxPlane - maxWinBegPlane;
  Int_t numPostMaxPlanes = maxWinEndPlane - maxPlane;
  Int_t numLongPlanes = numPostMaxPlanes + numPreMaxPlanes - 1;

  //Start Transverse Clustering
  
  //window info:
  std::vector<window> winvec;  //vector containing all transverse windows

  //Loop over all planes in long. window
  //Find transverse windows per plane

  for (int i=0;i<numLongPlanes;i++){  //loop over numLongPlanes
 
    //find next appropriate plane to look at starting from max plane:
    int pl = -1;
    if(numPreMaxPlanes>numPostMaxPlanes){
      if (i<numPreMaxPlanes) pl = maxPlane-i;
      else pl = maxPlane + i - numPreMaxPlanes + 1; 
      //add 1 to prevent re-doing max plane
    }
    else{
      if (i<numPostMaxPlanes) pl = maxPlane+i;
      else pl = maxPlane - i + numPostMaxPlanes - 1;
      //subtract 1 to prevent re-doing max plane
    }

    std::vector<Int_t> winu;  // vector of upper dip
    std::vector<Int_t> wind;  // vector of lower dip
    std::vector<Double_t> winu_tpos;  // vector of tpos of upper dip
    std::vector<Double_t> wind_tpos;  // vector of tpos of lower dip
    std::vector<Int_t> wintu; // vector of upper dip type
    std::vector<Int_t> wintd; // vector of lower dip type

    //find dips and categorize them
    MSG("SubShowerSR",Msg::kVerbose) << "Find Upper/Lower boundaries in plane = "
                                     << pl << endl;
    for(Int_t i=0;i<nstp;i++){
      if(plane[i]==pl){
        //using neighbouring strip info from above:
        //(only called dips if delta PH is significant given pe stats)
        if(TransGradientPlus[i]<0 /*-(TransGradientErrorPlus[i])*/ &&
           TransGradientMinus[i]<0 /*-(TransGradientErrorMinus[i])*/){
          winu.push_back(strip[i]);
          winu_tpos.push_back(tpos[i]);
          wintu.push_back(0);      // type-0 (this is a dip) 
          MSG("SubShowerSR",Msg::kVerbose) 
            << "up dip strip = " << strip[i] << " ph = " << ph[i] 
            << " tgp = " << TransGradientPlus[i]
            << " tgm = " << TransGradientMinus[i] << endl;
        }
        else if(TransGradientPlus[i]==ph[i] && 
                TransGradientMinus[i]<=ph[i]){
          winu.push_back(strip[i]);
          winu_tpos.push_back(tpos[i]);
          MSG("SubShowerSR",Msg::kVerbose) 
            << "up edge/iso strip = " << strip[i] << " ph = " << ph[i] 
            << " tgp = " << TransGradientPlus[i]
            << " tgm = " << TransGradientMinus[i] << endl;
          if (TransGradientMinus[i]<ph[i]) 
            wintu.push_back(1);     // type-1 (this is an edge strip) 
          else wintu.push_back(2);  // type-2 (this is an isolated strip)
        }
        if(TransGradientPlus[i]<0 /*-(TransGradientErrorPlus[i])*/ &&
           TransGradientMinus[i]<0 /*-(TransGradientErrorMinus[i])*/ ){
          wind.push_back(strip[i]);
          wind_tpos.push_back(tpos[i]);
          wintd.push_back(0);      // type-0
          MSG("SubShowerSR",Msg::kVerbose) 
            << "down dip strip = " << strip[i] << " ph = " << ph[i]
            << " tgp = " << TransGradientPlus[i]
            << " tgm = " << TransGradientMinus[i] << endl;
        }
        else if(TransGradientMinus[i]==ph[i] && 
                TransGradientPlus[i]<=ph[i]){
          wind.push_back(strip[i]);
          wind_tpos.push_back(tpos[i]);
          MSG("SubShowerSR",Msg::kVerbose) 
            << "down edge/iso strip = "  << strip[i] << " ph = " << ph[i]
            << " tgp = "  << TransGradientPlus[i]
            << " tgm = "  << TransGradientMinus[i] << endl;
          if (TransGradientPlus[i]<ph[i]) 
            wintd.push_back(1);     //type-1
          else wintd.push_back(2);  //type-2
        }
      }
    }

    //use dips to form windows                                   
    UInt_t wins = winu.size();
    if(wins!=wind.size()) {
      MSG("SubShowerSR",Msg::kError) 
        << "Number of Upper Transverse Window Boundaries does not equal " 
        << "number of Lower Transverse Window Boundaries... "
        << "Leaving FindCluster()" << endl;
      MSG("SubShowerSR",Msg::kError) << "winu.size() = " << winu.size() << endl;
      MSG("SubShowerSR",Msg::kError) << "wind.size() = " << wind.size() << endl;
      MSG("SubShowerSR",Msg::kError) << "wins = " << wins << endl;
      
      return false;
    }

    //temporarily hold window info
    std::vector<Int_t> win0(wins,0);             // strip of upper dip
    std::vector<Int_t> win1(wins,0);             // strip of lower dip
    std::vector<Float_t> win0_tpos(wins,0.0);      // tpos of strip of upper dip
    std::vector<Float_t> win1_tpos(wins,0.0);      // tpos of strip of lower dip
    std::vector<Float_t> win2(wins,0.0);           // PH of window
    std::vector<Int_t> win3(wins,0);             // ID of window:
                                                 //   0-both dips; 
                                                 //   1 or 10-one dip,one gap; 
                                                 //   11-both gaps; 
                                                 //   22-isolated strip
    std::vector<Float_t> win4(wins,0.0);           // Centroid of window in tpos
    std::vector<Float_t> win5(wins,0.0);           // PH weighted tpos of window
    std::vector<Float_t> win6(wins,0.0);           // PH weighted width of window

    //match up dips to form windows
    for (UInt_t w = 0; w<wins; w++){
      win0[w] = -1;
      win1[w] = -1;
      win0_tpos[w] = 0;
      win1_tpos[w] = 0;
      win2[w] = 0;
      win3[w] = -1;
      win4[w] = 0;
      win5[w] = 0;
      win6[w] = 0;
      int winw = 200;
      for (UInt_t s = 0; s<wind.size(); s++){
        int diff = winu.at(w)-wind.at(s);
        if ( ( (diff>0  && wintu.at(w)<2  && wintd.at(s)<2) ||
               (diff==0 && wintu.at(w)==2 && wintd.at(s)==2) ) &&
             diff<winw) {
          winw = diff;
          win0[w] = winu.at(w);
          win1[w] = wind.at(s);
          win0_tpos[w] = winu_tpos.at(w);
          win1_tpos[w] = wind_tpos.at(s);
          win3[w] = wintu.at(w)*10+wintd.at(s);
        }
      }
    }

    //now get PH of windows
    //and push_back windows into vectors
    for (UInt_t w = 0; w<wins; w++){            
      Double_t biggestStripPH = 0;
      Int_t n = 0;
      for(Int_t i=0;i<nstp;i++){
        if (plane[i]==pl && strip[i]<=win0[w] && strip[i]>=win1[w]) {
          win2[w] += ph[i];
          if(ph[i]>biggestStripPH) {
            win4[w] = tpos[i];
            biggestStripPH = ph[i];
          }
          win5[w] += ph[i]*tpos[i];
          win6[w] += ph[i]*tpos[i]*tpos[i];
          n++;
        }
      }
      win5[w] /= win2[w];
      if(n!=1){
        win6[w] /= win2[w];
        win6[w] -= win5[w]*win5[w];
        if(win6[w]>0) win6[w] = TMath::Sqrt(win6[w]);
        else win6[w] = STRIPWIDTHINMETRES/TMath::Sqrt(12.); 
     }
      else win6[w] = STRIPWIDTHINMETRES/TMath::Sqrt(12.);
      
      window win;
      win.plane = pl;
      win.upper = win0[w];
      win.lower = win1[w];
      win.upper_tpos = win0_tpos[w];
      win.lower_tpos = win1_tpos[w];
      win.ph = win2[w];
      win.type = win3[w];
      win.centroid = win4[w];
      win.phpos = win5[w];
      win.phwidth = win6[w];
      winvec.push_back(win);

    }
  }

  //Perform clustering on windows:

  // if this longitudinal window only contains a single 
  // active plane then just return largest window
  if(numLongPlanes/2<2){ 
    std::vector<window>::iterator winIter = winvec.begin();
    std::vector<window>::iterator winEnd = winvec.end();
    Int_t best_upper_strip = 0;
    Int_t best_lower_strip = 200;
    Double_t best_pulse_height = 0;
    while(winIter!=winEnd){
      if(winIter->ph>best_pulse_height){
        best_lower_strip = winIter->lower;
        best_upper_strip = winIter->upper;
        best_pulse_height = winIter->ph;
      }
      winIter++;
    }
    for (Int_t i=0; i<=allStrips->GetLast(); i++) {
      TObject *tobj = allStrips->At(i);
      assert(tobj->InheritsFrom("CandStripHandle"));
      CandStripHandle *stp = dynamic_cast<CandStripHandle*>(tobj);
      if(stp->GetPlane()>=maxWinBegPlane && stp->GetPlane()<=maxWinEndPlane &&
         stp->GetStrip()>=best_lower_strip && stp->GetStrip()<=best_upper_strip){
        newSubShower->Add(stp);
      }
    }
    return true;
  }

  //call transverse clustering algorithms
  std::vector<window> clusterWinVec = TransCluster(winvec,det);
  
  //no windows added to cluster:
  if(clusterWinVec.size()==0) return false;

  //Clean up proto-subshower
  //Check for longitudinal gaps within the proto-subshower
  //Remove windows after first gap on both sides of shower max plane
  maxPlane = -1;  //recalculate maxPlane
  maxPlanePH = 0.;
  Double_t PlanePH[500] = {0.};
  Int_t non0Win = 0;
  std::vector<window>::iterator winIter = clusterWinVec.begin();
  std::vector<window>::iterator winEnd = clusterWinVec.end();
  while(winIter!=winEnd){
    if(winIter->plane>0){
      PlanePH[winIter->plane] += winIter->ph;
      if(PlanePH[winIter->plane]>0) non0Win++;
      if(PlanePH[winIter->plane]>maxPlanePH) {
        maxPlane = winIter->plane;
        maxPlanePH = PlanePH[maxPlane];
      }
    }
    winIter++;
  }
  
  if(non0Win==0){
    MSG("SubShowerSR",Msg::kDebug) <<"No window solutions found"<<endl;
    return false;
  }

  Int_t maxLongPlane = -1;
  Int_t minLongPlane = -1;
  Int_t counter = 0;
  Int_t delta = 2;
  while((maxPlane+counter)>=0&&(PlanePH[maxPlane+counter]>0||(maxPlane+counter)==249)) {
    maxLongPlane = maxPlane+counter;
    delta = 2;
    if(det==1&&(maxLongPlane-249)*(maxLongPlane+delta-249)<0)
      delta = 3;
    counter+=delta;
  }
  counter = 0;
  while((maxPlane+counter)>=0&&(PlanePH[maxPlane+counter]>0||(maxPlane+counter)==249)) {
    minLongPlane = maxPlane+counter;
    delta = 2;
    if(det==1&&(minLongPlane-249)*(minLongPlane-delta-249)<0)
      delta = 3;
    counter-=delta;
  }
  MSG("SubShowerSR",Msg::kDebug) <<"New Longitudinal Window: "<< minLongPlane
                                 << "-" << maxLongPlane << endl;

  for (Int_t i=0; i<=allStrips->GetLast(); i++) {
    TObject *tobj = allStrips->At(i);
    assert(tobj->InheritsFrom("CandStripHandle"));
    CandStripHandle *stp = dynamic_cast<CandStripHandle*>(tobj);
    if(stp->GetPlane()>=minLongPlane && stp->GetPlane()<=maxLongPlane){
      winIter = clusterWinVec.begin();
      while(winIter!=winEnd){
        if(winIter->plane==stp->GetPlane()){
          Int_t strp = stp->GetStrip();
          if(strp<=winIter->upper && strp>=winIter->lower) {
            newSubShower->Add(stp);
            break;
          }
        }
        winIter++;
      }
    }
  }

  return true;
}

//******************************************************************
std::vector<window> AlgSubShowerSRList::TransCluster(std::vector<window> win,
                                                     Int_t det)
{
  std::vector<window> cluster;

  Int_t maxPlane = -1;  //calculate maxPlane
  Double_t maxPlanePH = 0;
  Double_t PlanePH[500] = {0};
  std::vector<window>::iterator winIter = win.begin();
  std::vector<window>::iterator winEnd = win.end();
  while(winIter!=winEnd){
    PlanePH[winIter->plane] += winIter->ph;
    if(PlanePH[winIter->plane]>maxPlanePH) {
      maxPlane = winIter->plane;
      maxPlanePH = PlanePH[maxPlane];
    }
    winIter++;
  }

  Int_t maxLongPlane = -1;
  Int_t minLongPlane = -1;
  Int_t counter = 0;
  Int_t delta = 2;
  while(PlanePH[maxPlane+counter]>0||(maxPlane+counter)==249) {
    maxLongPlane = maxPlane+counter;
    delta = 2;
    if(det==1&&(maxLongPlane-249)*(maxLongPlane+delta-249)<0)
      delta = 3;
    counter+=delta;
  }
  counter = 0;
  while(PlanePH[maxPlane+counter]<=0||(maxPlane+counter)==249) {
    minLongPlane = maxPlane+counter;
    delta = 2;
    if(det==1&&(minLongPlane-249)*(minLongPlane-delta-249)<0)
      delta = 3;
    counter-=delta;
  }
  Int_t numPreMaxPlanes = maxPlane-minLongPlane+1;
  Int_t numPostMaxPlanes = maxLongPlane-maxPlane+1;
  Int_t numLongPlanes = numPostMaxPlanes + numPreMaxPlanes - 1;

  //match up windows plane by plane
  Double_t shwmid = 0;    // centroid of cluster transversely
  Double_t shwwid = 0;    // trans width of widest window in proto-cluster
  Double_t maxwid = 0;
  Double_t clusterph = 0; // total PH of cluster
  Double_t refp = 0;      // reference plane  }  essentially pivot point for 
  Double_t refs = 0;      // reference strip  }  cluster
  
  //long. window info:  
  std::vector<Int_t> winpl(numLongPlanes,0);          // plane number 
  std::vector<Int_t> winpl0(numLongPlanes,0);         // upper bound
  std::vector<Int_t> winpl1(numLongPlanes,0);         // lower bound
  std::vector<Float_t> winpl0_tpos(numLongPlanes,0.0);  // upper bound in tpos
  std::vector<Float_t> winpl1_tpos(numLongPlanes,0.0);  // lower bound in tpos
  std::vector<Double_t> winpl2(numLongPlanes,0.0);      // PH
  std::vector<Int_t> winpl3(numLongPlanes,0);         // window type
  std::vector<Double_t> winpl4(numLongPlanes,0.0);      // width
  std::vector<Double_t> winpl5(numLongPlanes,0.0);      // centroid

  Double_t slopen = 0;
  Double_t eslopen = 0;
  Bool_t trustslopen = false;

  std::vector<window> houghcluster = HoughTransCluster(win,det);

  if(houghcluster.size()==0) { //houghcluster.size()=0 try usual tracking
    for (int d = 0;d<numLongPlanes;d++){
      int pl = -1;
      if(numPreMaxPlanes>numPostMaxPlanes){
        if (d<numPreMaxPlanes) pl = maxPlane-d;
        else pl = maxPlane+d-numPreMaxPlanes+1;
      }
      else{
        if (d<numPostMaxPlanes) pl = maxPlane+d;
        else pl = maxPlane-d+numPostMaxPlanes-1;
      }
      winpl[d] = 0;
      winpl0[d] = -1;
      winpl1[d] = -1;
      winpl0_tpos[d] = -10;
      winpl1_tpos[d] = -10;
      winpl2[d] = 0;
      winpl3[d] = 0;
      winpl4[d] = -1;
      winpl5[d] = -1;
      // PH for current best matched window:
      double weight = 0.;
      //deviation of best matched window from prediction of 
      //centroid for this plane:
      double weight1 = 10000.;
      //loop keeps track of # windows already checked on this plane:
      int loop = 0;
      // counts number of windows with a good match:
      int within = 0;
      //prediction of centroid for this plane (not using slope in info):
      double mid0 = 0;
      //prediction of centroid for this plane using slope info where available:
      double mid = 0;
      
      //loop over all in windows
      winIter = win.begin();
      while (winIter!=winEnd){
        if (winIter->plane==pl){  //if it's in current plane:
          loop++;
          double dev = shwwid;
          if (dev<(winIter->upper_tpos-winIter->lower_tpos+STRIPWIDTHINMETRES))
            dev = winIter->upper_tpos-winIter->lower_tpos+STRIPWIDTHINMETRES;
          if (dev<2.*STRIPWIDTHINMETRES) dev = 2.*STRIPWIDTHINMETRES;
          if (shwmid!=0) {
            mid0 = shwmid;
            mid = shwmid;
          }
          else {
            mid = (winIter->upper+winIter->lower)/2.;
            mid0 = mid;
          }
          
          // check if this window is consistent with previously selected 
          // windows and has a higher PH than current best match for this
          // plane or else if it's shower max plane
          if((shwmid==0 ||
              (TMath::Abs(winIter->upper_tpos-mid)<dev &&
               TMath::Abs(winIter->lower_tpos-mid)<dev))
             && winIter->ph>weight){
            
            //loop through planes in long. window again to do a local 
            //scan to ensure best match is found.
            //step forward and backward once each if possible
          
            int loopp = 0;
            int loopm = 0;
            int withinp = 0;
            int withinm = 0;
            //double midp0 = 0;
            //double midm0 = 0;
            double midp = 0;
            double midm = 0;
            
            // step forward:
            std::vector<window>::iterator winIter1 = win.begin();
            while(winIter1!=winEnd){
              if ((((det==1&&(pl-249)*(pl+2-249)>0)||det==0)
                   &&winIter1->plane==pl+2)
                  ||(det==1&&(pl-249)*(pl+2-249)<0
                     &&winIter1->plane==pl+3)){   //go to next same view plane
                loopp++;
                if(shwmid!=0){  //not shower max plane
                  midp = shwmid;
                  if(TMath::Abs(winIter1->upper_tpos-midp)<dev/2.
                     ||TMath::Abs(winIter1->lower_tpos-midp)<dev/2.){  //good match
                    withinp++;
                  }
                }
                else if(shwmid==0){  //in shower max plane
                  midp = (winIter->upper_tpos+winIter->lower_tpos)/2.;
                  if(TMath::Abs(winIter1->upper_tpos-midp)<dev/2.
                     ||TMath::Abs(winIter1->lower_tpos-midp)<dev/2.){  //good match
                    withinp++;  
                  }
                }
              }
              
              //step backwards:
              if ((((det==1&&(pl-249)*(pl-2-249)>0)||det==0)
                   &&winIter1->plane==pl-2)
                  ||(det==1&&(pl-249)*(pl-2-249)<0
                     &&winIter1->plane==pl-3)){   
                //go to previous same view plane                
                loopm++;
                if(shwmid!=0){   // not shower max plane
                  midm = shwmid;
                  if(TMath::Abs(winIter1->upper_tpos-midm)<dev
                     ||TMath::Abs(winIter1->lower_tpos-midm)<dev){  //good match
                    withinm++;
                  }
                }
                else if(shwmid==0){  //shower max plane
                  midm = (winIter->upper_tpos+winIter->lower_tpos)/2.;
                  if(TMath::Abs(winIter1->upper_tpos-midm)<dev
                     ||TMath::Abs(winIter1->lower_tpos-midm)<dev){  //good match
                    withinm++;
                  }
                }
              }
              winIter1++;
            }
            if((withinp>0 && withinm>0) || (loopp==0 && withinm>0) ||
               (loopm==0 && withinp>0)){  
              //success! current window is a good solution
              within++;
              // write out solution: 
              winpl[d] = winIter->plane;
              winpl0[d] = winIter->upper;
              winpl1[d] = winIter->lower;
              winpl0_tpos[d] = winIter->upper_tpos;
              winpl1_tpos[d] = winIter->lower_tpos;
              winpl2[d] = winIter->ph;
              winpl3[d] = winIter->type;
              weight = winIter->ph;
              weight1 = (winIter->upper_tpos+winIter->lower_tpos)/2.-mid0;
              //these will be overwritten if there is another window 
              //which is also a good match, but has more PH, 
              //ensuring we get best possible match
            }
          }
        }
        winIter++;
      }
      
      if(loop*within>0){   //got a matched window
        if (shwmid==0) {   // shower max plane
          refp = winpl[d];  // set reference point
          refs = (winpl0_tpos[d]+winpl1_tpos[d])/2.;
        }
        //determine widest window:
        if(shwwid<(winpl0_tpos[d]-winpl1_tpos[d]+STRIPWIDTHINMETRES)) 
          shwwid = winpl0_tpos[d]-winpl1_tpos[d]+STRIPWIDTHINMETRES; 
        // average centroid of event so far:
        shwmid = (shwmid*clusterph+(winpl0_tpos[d]+winpl1_tpos[d]) * 
                  winpl2[d]/2.)/(clusterph+winpl2[d]); 
        clusterph+=winpl2[d];
        winpl4[d] = shwwid;
        winpl5[d] = shwmid;
      }
    }
  
    //Fitting the boundaries of windows matched
    Double_t slopetn = 0;
    Double_t slopebn = 0;
    Double_t eslopetn = 0;
    Double_t eslopebn = 0;
    Double_t chi2tn = 0;
    Double_t chi2bn = 0;
    int nt = 0;
    int nb = 0;

    for (int d = 0;d<numLongPlanes;d++){
      int pl = -1;
      if(numPreMaxPlanes>numPostMaxPlanes){
        if (d<numPreMaxPlanes) pl = maxPlane-d;
        else pl = maxPlane+d-numPreMaxPlanes+1;
      }
      else{
        if (d<numPostMaxPlanes) pl = maxPlane+d;
        else pl = maxPlane-d+numPostMaxPlanes-1;
      }
      if(winpl2[d]>0) {nt++; nb++;}
      if (maxwid<winpl4[d]) maxwid = winpl4[d];
      if (maxwid<2.*STRIPWIDTHINMETRES) maxwid = 2.*STRIPWIDTHINMETRES;
    }
    if (nt>2 && nb>2) {
      TGraphErrors grat(nt);
      TGraphErrors grab(nb);
      int t = 0;
      int b = 0;
      for (int d = 0;d<numLongPlanes;d++){
        int pl = -1;
        if(numPreMaxPlanes>numPostMaxPlanes){
          if (d<numPreMaxPlanes) pl = maxPlane-d;
          else pl = maxPlane+d-numPreMaxPlanes+1;
        }
        else{
          if (d<numPostMaxPlanes) pl = maxPlane+d;
          else pl = maxPlane-d+numPostMaxPlanes-1;
        }
        if(winpl2[d]>0) {
          grat.SetPoint(t,winpl[d],winpl0_tpos[d]);
          grat.SetPointError(t,0,1./TMath::Sqrt(winpl2[d]));
          t++;
          grab.SetPoint(b,winpl[d],winpl1_tpos[d]);
          grab.SetPointError(b,0,1./TMath::Sqrt(winpl2[d]));
          b++;
        }
      }

      grat.Fit("ssPol1","NNQ");
      slopetn = ssPol1->GetParameter(1);
      eslopetn = ssPol1->GetParError(1);
      chi2tn = ssPol1->GetChisquare();
      grab.Fit("ssPol1","NNQ");
      slopebn = ssPol1->GetParameter(1);
      eslopebn = ssPol1->GetParError(1);
      chi2bn = ssPol1->GetChisquare();
      
      if ((slopetn+slopebn)!=0. && slopetn*slopebn!=0.){
        if(chi2tn/nt<100. && chi2bn/nb<100. &&
           ((((eslopetn/TMath::Abs(slopetn)<0.1 && eslopebn/TMath::Abs(slopebn)<0.1
               &&chi2tn/nt>0.1 && chi2bn/nb>0.1)) &&
             (TMath::Abs(slopetn-slopebn)/TMath::Abs(slopetn+slopebn)<0.25)) ||
            TMath::Abs(slopetn-slopebn)==0)) {
          slopen = (slopetn+slopebn)/2.;
          eslopen = TMath::Sqrt(TMath::Power(eslopetn,2)+TMath::Power(eslopebn,2))/2.;
          if(TMath::Abs(slopetn-slopebn)>eslopen) eslopen = TMath::Abs(slopetn-slopebn);
          trustslopen = true;
        }
      }
      MSG("SubShowerSR",Msg::kDebug) <<"fit window "<<slopetn<<" "<<slopebn
                                     <<" "<<slopen<<" "<<eslopetn<<" "
                                     <<eslopebn<<" "<<eslopen<<" "
                                     <<chi2tn/nt<<" "<<chi2bn/nb<<" "
                                     <<nt<<" "<<nb<<endl;
    }
    
    //check if solution is good:
    double usewid = maxwid;
    for (Int_t d = 0;d<numLongPlanes;d++){
      int pl = -1;
      if(numPreMaxPlanes>numPostMaxPlanes){
        if (d<numPreMaxPlanes) pl = maxPlane-d;
        else pl = maxPlane+d-numPreMaxPlanes+1;
      }
      else{
        if (d<numPostMaxPlanes) pl = maxPlane+d;
        else pl = maxPlane-d+numPostMaxPlanes-1;
      }
    
      if(d>0 && (winpl0[d]>=0||winpl1[d]>=0) ){ //if there is a solution
        if(!trustslopen && //don't trust new slope
           //window is too far away from centroid
           (winpl0_tpos[d]-shwmid)*(winpl1_tpos[d]-shwmid) > 
           (maxwid/2.)*(maxwid/2.) ){ 
          MSG("SubShowerSR",Msg::kDebug) << "removed plane " << winpl[d] 
                                         << " " << winpl0[d] << " " 
                                         << winpl1[d] << " " 
                                         << winpl0_tpos[d] << " " 
                                         << winpl1_tpos[d] << endl;
          shwmid = (shwmid*clusterph-(winpl0_tpos[d]+winpl1_tpos[d]) * 
                    winpl2[d]/2.)/(clusterph-winpl2[d]);
          winpl0[d] = -1;
          winpl1[d] = -1;
          winpl0_tpos[d] = -10;
          winpl1_tpos[d] = -10;
          winpl2[d] = 0;
          winpl3[d] = 0;
          winpl4[d] = -1;
          winpl5[d] = -1;
        }
        else if(trustslopen){ 
          //this time account for slope on window boundaries:
          //try to find another window as a reference point       
          //to compare the window in this plane with
          int i = 0;  
          while(winpl2[i]==0&&i<=d) i++;
          if(i==d){
            while(winpl2[i]==0&&i<numLongPlanes) i++;
          }
          if(usewid<TMath::Abs(winpl[d]-winpl[i])*eslopen)  
            //find appropriate cut value on window width
            usewid = TMath::Abs(winpl[d]-winpl[i])*eslopen; 
          //expected centroid for plane:
          double expmid = (winpl0_tpos[i] + 
                           winpl1_tpos[i])/2.+(winpl[d]-winpl[i])*slopen; 
          if((winpl0_tpos[d]-expmid)*(winpl1_tpos[d]-expmid) >
             (usewid/2.)*(usewid/2.)){ 
            MSG("SubShowerSR",Msg::kDebug) <<"removed plane using slope "
                                           <<winpl[d]<<" "<<winpl0[d]<<" "
                                           <<winpl1[d]<<" with pl "
                                           <<maxPlane<<" "<<winpl[i]
                                           <<" "<<expmid<<" "<<usewid<<endl;
            shwmid = (shwmid*clusterph-(winpl0_tpos[d]+winpl1_tpos[d]) * 
                      winpl2[d]/2.)/(clusterph-winpl2[d]);
            clusterph-=winpl2[d];
            winpl0[d] = -1;
            winpl1[d] = -1;
            winpl0_tpos[d] = -10;
            winpl1_tpos[d] = -10;
            winpl2[d] = 0;
            winpl3[d] = 0;
            winpl4[d] = -1;
            winpl5[d] = -1;
          }
        }
      }
    }
  }
  else { //use hough solution as a starting point:
    for (int d = 0;d<numLongPlanes;d++){  
      int pl = -1;
      if(numPreMaxPlanes>numPostMaxPlanes){
        if (d<numPreMaxPlanes) pl = maxPlane-d;
        else pl = maxPlane+d-numPreMaxPlanes+1;
      }
      else{
        if (d<numPostMaxPlanes) pl = maxPlane+d;
        else pl = maxPlane-d+numPostMaxPlanes-1;
      }
      winpl[d] = 0;
      winpl0[d] = -1;
      winpl1[d] = -1;
      winpl0_tpos[d] = -10;
      winpl1_tpos[d] = -10;
      winpl2[d] = 0;
      winpl3[d] = 0;
      winpl4[d] = -1;
      winpl5[d] = -1;

      //loop over all in windows
      std::vector<window>::iterator hcIter = houghcluster.begin();
      std::vector<window>::iterator hcEnd = houghcluster.end();      
      while (hcIter!=hcEnd){
        if (hcIter->plane==pl){  //if it's in current plane:
          winpl[d] = hcIter->plane;
          winpl0[d] = hcIter->upper;
          winpl1[d] = hcIter->lower;
          winpl0_tpos[d] = hcIter->upper_tpos;
          winpl1_tpos[d] = hcIter->lower_tpos;
          winpl2[d] = hcIter->ph;
          winpl3[d] = hcIter->type;
                
          if (shwmid==0) {   // shower max plane
            refp = winpl[d];  // set reference point
            refs = (winpl0_tpos[d]+winpl1_tpos[d])/2.;
          }
          //determine widest window:
          if(shwwid<(winpl0_tpos[d]-winpl1_tpos[d]+STRIPWIDTHINMETRES)) 
            shwwid = winpl0_tpos[d]-winpl1_tpos[d]+STRIPWIDTHINMETRES; 
          // average centroid of event so far:
          shwmid = (shwmid*clusterph+(winpl0_tpos[d]+winpl1_tpos[d]) * 
                    winpl2[d]/2.)/(clusterph+winpl2[d]); 
          // calculate average distance to centroid:
          clusterph += winpl2[d];
          winpl4[d] = shwwid;
          winpl5[d] = shwmid;
          if (maxwid<winpl4[d]) maxwid = winpl4[d];
          if (maxwid<2.*STRIPWIDTHINMETRES) maxwid = 2.*STRIPWIDTHINMETRES;
        }
        hcIter++;
      }
    }
    std::vector<Double_t> MCV = BestHough(houghcluster,false);
    if(MCV[5]!=0){
      slopen = MCV[5];
      eslopen = MCV[8];
      trustslopen = true;
    }
    else if(MCV[0]!=0){
      slopen = MCV[0];
      eslopen = MCV[3];
      trustslopen = true;
    }
    else {
      slopen = 0;
      eslopen = 0;
      trustslopen = false;
    }
  }
  
  //Try to fill planes with no windows
  
  //First find centroid and width for +/- 1 planes
  
  for (Int_t d = 0;d<numLongPlanes;d++){
    int pl = -1;
    if(numPreMaxPlanes>numPostMaxPlanes){
      if (d<numPreMaxPlanes) pl = maxPlane-d;
      else pl = maxPlane+d-numPreMaxPlanes+1;
    }
    else{
      if (d<numPostMaxPlanes) pl = maxPlane+d;
      else pl = maxPlane-d+numPostMaxPlanes-1;
    }
    
    double minwei = 200.;   //minimum distance to centroid (weight)
    double minwei_tpos = 10.;   //minimum distance to centroid (weight)
    int newwinpl0 = -1;     //new window upper and lower strip bounds
    int newwinpl1 = -1;
    Float_t newwinpl0_tpos = -10; //new window upp and low strip tpos bounds
    Float_t newwinpl1_tpos = -10;
    double newwinpl2 = 0.;
    int newwinpl3 = 0;

    if(winpl0[d]==-1&&winpl1[d]==-1){  //does this plane have a solution?
      
      int p = 0;
      int m = 0;
      if (det==0||(det==1&&(pl-249)*(pl+2-249)>0)) p = 2;
      else p = 3;
      if (det==0||(det==1&&(pl-249)*(pl-2-249)>0)) m = 2;
      else m = 3;
      
      double pmid = -1;   //centroid for next plane (plus)
      double mmid = -1;   //centroid for prev plane (minus)
      double pwid = -1;   //width
      double mwid = -1;   
      int pwt = -1;      // window type  plus
      int mwt = -1;     //               minus
      bool useslopep = false; 
      bool useslopem = false;
      
      for (Int_t c = 0;c<numLongPlanes;c++){
        if (winpl[c]==pl+p) {
          pmid = winpl5[c];
          if (trustslopen) useslopep = true;
          pwid = winpl4[c];
          pwt = winpl3[c];
        }
        if (winpl[c]==pl-m) {
          mmid = winpl5[c];      
          if (trustslopen) useslopem = true;
          mwid = winpl4[c];
          mwt = winpl3[c];
        }
      }
      
      //Find compatible windows
      
      double usewid = maxwid;
      if(maxwid>0){
        if((mmid==-1&&pmid==-1)||
           (!useslopem&&mmid>=0&&TMath::Abs(shwmid-mmid)>maxwid)||
           (!useslopep&&pmid>=0&&TMath::Abs(shwmid-pmid)>maxwid)||
           (useslopem&&mmid>=0&&
              (TMath::Abs(shwmid-mmid-(maxPlane-pl+m)*slopen)>maxwid)||
            (TMath::Abs(shwmid-mmid-(maxPlane-pl+m)*slopen) > 
             TMath::Abs(maxPlane-pl+m)*eslopen))||
           (useslopep&&pmid>=0&&
              (TMath::Abs(shwmid-pmid-(maxPlane-pl-p)*slopen)>maxwid)||
            (TMath::Abs(shwmid-pmid-(maxPlane-pl-p)*slopen) > 
             TMath::Abs(maxPlane-pl-p)*eslopen))){
          minwei = 200.;
          minwei_tpos = 10.;
          newwinpl0 = -1;
          newwinpl1 = -1;
          newwinpl0_tpos = -10;
          newwinpl1_tpos = -10;
          newwinpl2 = 0.;
          winIter = win.begin();
          while(winIter!=winEnd){
            if(winIter->plane==pl){
              if (//shwmid>=0&&
                  ((winIter->upper_tpos-shwmid) * 
                   (winIter->lower_tpos-shwmid)<=
                   (maxwid/2.)*(maxwid/2.)) && 
                  TMath::Abs((winIter->upper_tpos + 
                              winIter->lower_tpos)/2.-shwmid)<minwei){
                minwei = TMath::Abs((winIter->upper_tpos + 
                                     winIter->lower_tpos)/2.-shwmid);
                newwinpl0 = winIter->upper;
                newwinpl1 = winIter->lower;
                newwinpl0_tpos = winIter->upper_tpos;
                newwinpl1_tpos = winIter->lower_tpos;
                newwinpl2 = winIter->ph;
                newwinpl3 = winIter->type;
              }
            }
            winIter++;
          }
        }
        else{
          minwei = 200.;
          minwei_tpos = 10.;
          newwinpl0 = -1;
          newwinpl1 = -1;
          newwinpl0_tpos = -10.;
          newwinpl1_tpos = -10.;
          newwinpl2 = 0.;
          winIter = win.begin();
          while(winIter!=winEnd){
            if(winIter->plane==pl){
              if (!useslopem&&mmid>=0
                  &&((winIter->upper_tpos-mmid) * 
                     (winIter->lower_tpos-mmid) <=
                     (maxwid/2.)*(maxwid/2.)) &&
                  TMath::Abs((winIter->upper_tpos + 
                              winIter->lower_tpos)/2.-mmid)<minwei){
                minwei = TMath::Abs((winIter->upper_tpos + 
                                     winIter->lower_tpos)/2.-mmid);
                newwinpl0 = winIter->upper;
                newwinpl1 = winIter->lower;
                newwinpl0_tpos = winIter->upper_tpos;
                newwinpl1_tpos = winIter->lower_tpos;
                newwinpl2 = winIter->ph;
                newwinpl3 = winIter->type;
              }
              else if (!useslopep&&pmid>=0
                       &&((winIter->upper_tpos-pmid) * 
                          (winIter->lower_tpos-pmid) <= 
                          (maxwid/2.)*(maxwid/2.))
                       && TMath::Abs((winIter->upper_tpos + 
                                      winIter->lower_tpos)/2.-pmid)<minwei){
                minwei = TMath::Abs((winIter->upper_tpos + 
                                     winIter->lower_tpos)/2.-pmid);
                newwinpl0 = winIter->upper;
                newwinpl1 = winIter->lower;
                newwinpl0_tpos = winIter->upper_tpos;
                newwinpl1_tpos = winIter->lower_tpos;
                newwinpl2 = winIter->ph;
                newwinpl3 = winIter->type;
              }
              else if (useslopem&&mmid>=0){
                if(usewid<m*eslopen) usewid = m*eslopen;
                else usewid = maxwid;
                if((winIter->upper_tpos-mmid-m*slopen) * 
                   (winIter->lower_tpos-mmid-m*slopen) <= 
                   (usewid/2.)*(usewid/2.)
                   &&TMath::Abs((winIter->upper_tpos + 
                                 winIter->lower_tpos)/2.-mmid-m*slopen)<minwei){
                  minwei = TMath::Abs((winIter->upper_tpos + 
                                       winIter->lower_tpos)/2.-mmid-m*slopen);
                  newwinpl0 = winIter->upper;
                  newwinpl1 = winIter->lower;
                  newwinpl0_tpos = winIter->upper_tpos;
                  newwinpl1_tpos = winIter->lower_tpos;
                  newwinpl2 = winIter->ph;
                  newwinpl3 = winIter->type;
                }
              }
              else if (useslopep&&pmid>=0){
                if(usewid<p*eslopen) usewid = p*eslopen;
                else usewid = maxwid; 
                if((winIter->upper_tpos-pmid+p*slopen) * 
                   (winIter->lower_tpos-pmid+p*slopen) <= 
                   (usewid/2.)*(usewid/2.)
                   &&TMath::Abs((winIter->upper_tpos + 
                                 winIter->lower_tpos)/2.-pmid+p*slopen)<minwei){
                  minwei = TMath::Abs((winIter->upper_tpos + 
                                       winIter->lower_tpos)/2.-pmid+p*slopen);
                  newwinpl0 = winIter->upper;
                  newwinpl1 = winIter->lower;
                  newwinpl0_tpos = winIter->upper_tpos;
                  newwinpl1_tpos = winIter->lower_tpos;
                  newwinpl2 = winIter->ph;
                  newwinpl3 = winIter->type;
                }
              }
            }
            winIter++;
          }
        }
      }
      if(newwinpl0>=0) {
        winpl0[d] = newwinpl0;
        winpl1[d] = newwinpl1;
        winpl0_tpos[d] = newwinpl0_tpos;
        winpl1_tpos[d] = newwinpl1_tpos;
        winpl2[d] = newwinpl2;
        winpl3[d] = newwinpl3;
        shwmid = (shwmid*clusterph + 
                  (winpl0_tpos[d]+winpl1_tpos[d]) * 
                  winpl2[d]/2.)/(clusterph+winpl2[d]);
        clusterph+=winpl2[d];
        MSG("SubShowerSR",Msg::kDebug) <<"find missed windows "<<pl
                                       <<" "<<newwinpl0<<" "
                                       <<newwinpl1<<" "<<minwei<<endl;
      }
    }
    
    window winSolution;
    winSolution.plane = pl;
    winSolution.upper = winpl0[d];
    winSolution.lower = winpl1[d];
    winSolution.upper_tpos = winpl0_tpos[d];
    winSolution.lower_tpos = winpl1_tpos[d];
    winSolution.ph = winpl2[d];
    winSolution.type = winpl3[d];
    winSolution.centroid = 0.;
    winSolution.phpos = 0.;
    winSolution.phwidth = 0.;
    cluster.push_back(winSolution);
  }

  return cluster;
}


//******************************************************************
std::vector<window> AlgSubShowerSRList::HoughTransCluster(std::vector<window> win,
                                                          Int_t det)
{

  MSG("SubShowerSR",Msg::kDebug) << "In HoughTransCluster()" << endl;
  
  std::vector<window> cluster;
  std::vector<window> nonCluster;
  std::vector<window> emptyCluster;
  
  //copy all windows into nonCluster:
  //also get best vertex estimate in first plane
  Int_t vertexPlane = 500;
  Double_t vertexPlanePH = 0;
  Double_t vertexTPos = 0;
  Double_t vertexTPosWidth = 0;

  //use these to calculate the density of hits in the long. window
  //this can be used to decide which Hough method to use
  Double_t totStrips = 0; 
  Double_t maxTpos = -5.0;
  Double_t minTpos = 5.0;
  Int_t lastPlane = 0;

  std::vector<window>::iterator wIt = win.begin();
  std::vector<window>::iterator wEn = win.end();
  while(wIt!=wEn){
    nonCluster.push_back(*wIt);
    if(wIt->plane<vertexPlane) {
      if(wIt->ph>1.0) {
        vertexPlane = wIt->plane;
        vertexPlanePH = wIt->ph;
        vertexTPos = wIt->phpos;
        vertexTPosWidth = wIt->phwidth;
        if(vertexTPosWidth<STRIPWIDTHINMETRES) vertexTPosWidth=STRIPWIDTHINMETRES;
      }
      else if(vertexPlanePH<1.0){
        vertexPlane = wIt->plane;
        vertexPlanePH = wIt->ph;
        vertexTPos = wIt->phpos;
        vertexTPosWidth = wIt->phwidth;
        if(vertexTPosWidth<STRIPWIDTHINMETRES) vertexTPosWidth=STRIPWIDTHINMETRES;
      }
    }
    else if(wIt->plane==vertexPlane){
      if(wIt->ph>vertexPlanePH){
        vertexPlanePH = wIt->ph;
        vertexTPos = wIt->phpos;
        vertexTPosWidth = wIt->phwidth;
        if(vertexTPosWidth<STRIPWIDTHINMETRES) vertexTPosWidth=STRIPWIDTHINMETRES;
      }
    }
    if(wIt->ph>1.0) {
      if(wIt->plane>lastPlane) lastPlane = wIt->plane;
      if(wIt->upper_tpos>maxTpos) maxTpos = wIt->upper_tpos;
      if(wIt->lower_tpos<minTpos) minTpos = wIt->lower_tpos;      
      totStrips += Double_t(wIt->upper - wIt->lower + 1);
    }
    wIt++;
  }

  //keep track of some quantities for quality checking at the end:
  Int_t   lowPlane  = 500;
  Int_t   uppPlane  = 0;
  Double_t nStrips  = 0;
  Double_t   avePH  = 0;
  Double_t aveWidth = 0;
  
  //keep looping until we have all windows
  Bool_t keepGoing = true;
  Int_t nLoops = 0;
  while(keepGoing && nLoops<MAXHOUGHITER) {
    
    MSG("SubShowerSR",Msg::kDebug) << "Running BestHough(): Loop number " 
                                   << nLoops << endl;

    //remember number of strips in cluster after last loop
    //once this value does not change, then stop looping
    Double_t oldNStrips = nStrips;

    //Get Hough m,c:
    std::vector<Double_t> MCV(14,0.0);
    Int_t badGradient[4] = {0};

    //If we are on the first loop, use all windows to get angle
    //If we are not on first loop, just use previously 
    //selected windows to get better determination of angle
    if(cluster.size()==0){
      MSG("SubShowerSR",Msg::kDebug) << "cluster.size()=0" << endl;
      MCV = BestHough(win,false);

      //verify that best hough passes reasonably close to vertex:
      if(true){
        MSG("SubShowerSR",Msg::kDebug) << "vertexPlane = " 
                                       << vertexPlane << endl;
        MSG("SubShowerSR",Msg::kDebug) << "vertexTPos = " 
                                       << vertexTPos << endl;
        MSG("SubShowerSR",Msg::kDebug) << "vertexTPosWidth = " 
                                       << vertexTPosWidth << endl;
        
        double vertexY_simple = (vertexPlane - MCV[2])*MCV[0] + MCV[1];
        double vertexY_simple_peak = (vertexPlane - MCV[2])*MCV[10] + MCV[11];
        
        MSG("SubShowerSR",Msg::kDebug) << "vertexY_simple = " 
                                       << vertexY_simple << endl;
        
        if( vertexY_simple < vertexTPos - vertexTPosWidth ||
            vertexY_simple > vertexTPos + vertexTPosWidth ||
            MCV[0]==0 ) badGradient[0] = 1;
        if( vertexY_simple_peak < vertexTPos - vertexTPosWidth || 
            vertexY_simple_peak > vertexTPos + vertexTPosWidth || 
            MCV[10]==0 ) badGradient[1] = 1;
    
        double vertexY_ph      = (vertexPlane - MCV[7])*MCV[5] + MCV[6];
        double vertexY_ph_peak = (vertexPlane - MCV[7])*MCV[12] + MCV[13];
        MSG("SubShowerSR",Msg::kDebug) << "vertexY_ph = " 
                                       << vertexY_ph << endl;
        if( vertexY_ph < vertexTPos - vertexTPosWidth || 
            vertexY_ph > vertexTPos + vertexTPosWidth ||
            MCV[5]==0 ) badGradient[2] = 1;
        if( vertexY_ph_peak < vertexTPos - vertexTPosWidth || 
            vertexY_ph_peak > vertexTPos + vertexTPosWidth ||
            MCV[12]==0 ) badGradient[3] = 1;
      }
    }
    else {
      MSG("SubShowerSR",Msg::kDebug) << "cluster.size()=" 
                                     << cluster.size() << endl;
      MCV = BestHough(cluster,false);
    }

    MSG("SubShowerSR",Msg::kDebug) 
      << "\nHough gradient, intercept, plane vertex: " 
      << MCV[0] << "+/-" << MCV[3] << " " 
      << MCV[1] << "+/-" << MCV[4] << " " << MCV[2]
      << "\nPH Hough gradient, intercept, plane vertex: " 
      << MCV[5] << "+/-" << MCV[8] << " " 
      << MCV[6] << "+/-" << MCV[9] << " " << MCV[7]
      << "\n Hough Max Bin Gradient and Intercept Values: "
      << MCV[10] << " " << MCV[11] << " "
      << MCV[12] << " " << MCV[13] << endl;
    
    //if we've learnt nothing from first run of BestHough
    //stop now and try other clustering algorithm
    if(badGradient[0]==1 && badGradient[1]==1 && 
       badGradient[2]==1 && badGradient[3]==1 && 
       cluster.size()==0) return emptyCluster;

    cluster.clear();
    avePH = 0;
    aveWidth = 0;
    nStrips = 0;
    lowPlane = 500;
    uppPlane = 0;
    std::vector<int> tempcluster[4];
    std::vector<window>::iterator winIter = nonCluster.begin();
    std::vector<window>::iterator winEnd  = nonCluster.end();
    Int_t counter = 0;
    while(winIter!=winEnd){
      double y[4] = {};
      double y_upp[4] = {};
      double y_low[4] = {};
      y[0] = (winIter->plane - MCV[2])*MCV[0]  + MCV[1];
      y[1] = (winIter->plane - MCV[2])*MCV[10] + MCV[11];
      y[2] = (winIter->plane - MCV[7])*MCV[5]  + MCV[6];
      y[3] = (winIter->plane - MCV[7])*MCV[12] + MCV[13];
      
      if(cluster.size()==0) {
        for(int ii=0;ii<4;ii++) {
          y_upp[ii] = y[ii] + STRIPWIDTHINMETRES; 
          y_low[ii] = y[ii] - STRIPWIDTHINMETRES;
        }
      }
      else {
        for(int ii=0;ii<4;ii++) {
          y_upp[ii] = y[ii] + MCV[4+int(ii/2)*5];
          y_low[ii] = y[ii] - MCV[4+int(ii/2)*5];
        }
      }
      
      for(int ii=0;ii<4;ii++){
        if( ( winIter->upper_tpos>=y[ii] && winIter->lower_tpos<=y[ii] ) ||
            ( winIter->phpos<=y_upp[ii]  && winIter->phpos>=y_low[ii]  ) ) {
          if(badGradient[ii]==0) {
            tempcluster[ii].push_back(counter);
          }
        }
      }
      winIter++;
      counter++;
    }
    
    Int_t best_cluster = -1;
    UInt_t max_size = 0;
    for(int ii=0;ii<4;ii++){
      if(tempcluster[ii].size()>max_size) {
        max_size = tempcluster[ii].size();
        best_cluster = ii;
      }
    }
    if(best_cluster>=0 && best_cluster<=3) {
      std::vector<window>::iterator winIter = nonCluster.begin();
      std::vector<window>::iterator winEnd  = nonCluster.end();
      std::vector<int>::iterator tempIter = tempcluster[best_cluster].begin();
      std::vector<int>::iterator tempEnd  = tempcluster[best_cluster].end();
      counter = 0;
      while(winIter!=winEnd){
        if(tempIter!=tempEnd && (*tempIter) == counter) {
          cluster.push_back(*winIter);
          avePH += winIter->ph;
          aveWidth += winIter->upper_tpos - winIter->lower_tpos;
          nStrips += Double_t(winIter->upper - winIter->lower + 1);
          if(winIter->plane<lowPlane) lowPlane = winIter->plane;
          if(winIter->plane>uppPlane) uppPlane = winIter->plane;
          tempIter++;
        }
        winIter++;
        counter++;
      }
    }
    
    MSG("SubShowerSR",Msg::kDebug) << "Proto-SubShower currently has " 
                                   << (uppPlane-lowPlane+2)/2
                                   << " planes and " << nStrips 
                                   << " strips" << endl;

    if(nStrips == oldNStrips) keepGoing = false;
    nLoops+=1;
  }

  if(nStrips<=0) return emptyCluster;
    
  MSG("SubShowerSR",Msg::kDebug) 
    << "Returning Proto-SubShower from HoughTransCluster()" << endl;
  
  if( det==1 && uppPlane>249 && lowPlane<249 ) return cluster;
  avePH/=nStrips;
  if( ((uppPlane - lowPlane + 2)/2)<fMinHoughPlanes && 
      avePH<fMinHoughPH ) return emptyCluster;
  
  return cluster;

}

//******************************************************************
Bool_t AlgSubShowerSRList::TestOverLap(TObjArray *newSubShower,CandHandle &ch,
                                       const CandSliceHandle *slice)
{
  //Going to check that proto-subshower does not go through an existing subshower  
  //if it does, take largest, plane-wise contiguous set of strips

  MSG("SubShowerSR", Msg::kDebug) << "Running TestOverLap()" << endl;

  //if no strips return false. Empty subshower will be caught later.
  Int_t totStp = newSubShower->GetLast()+1;  
  if(totStp<=1) return false;

  assert(ch.InheritsFrom("CandSubShowerSRListHandle"));
  CandSubShowerSRListHandle &csslh = dynamic_cast<CandSubShowerSRListHandle &>(ch);
  //if this is the first subshower, just return false:
  if(csslh.GetNDaughters()==0) return false;

  //Get plane view:
  TObject *firstObj = newSubShower->At(0);
  assert(firstObj->InheritsFrom("CandStripHandle"));
  CandStripHandle *firstStp = dynamic_cast<CandStripHandle*>(firstObj);
  PlaneView::PlaneView_t pv = firstStp->GetPlaneView();
  MSG("SubShowerSR", Msg::kDebug) << "PlaneView = " << pv << endl;

  //loop through subshowers in view + energy order
  //so that newSubShower is tested against largest subshower first
  CandSubShowerSRHandleItr subshowerItr(csslh.GetDaughterIterator());
  CandSubShowerSRHandleKeyFunc *engKF = subshowerItr.CreateKeyFunc();
  engKF->SetFun(CandSubShowerSRHandle::KeyFromViewEnergy);
  subshowerItr.GetSet()->AdoptSortKeyFunc(engKF);
  while (CandSubShowerSRHandle *subshower = subshowerItr()) {
    if((*subshower->GetCandSlice())!=(*slice)) continue;
    if(subshower->GetPlaneView()!=pv) continue;
    if(subshower->GetClusterID()==ClusterType::kXTalk ||
       subshower->GetClusterID()==ClusterType::kHalo) continue;
    //Get vertex and angle of subshower:
    Double_t tpos_vtx = 0;
    Double_t z_vtx = subshower->GetVtxZ();
    if((subshower->GetPlaneView()==PlaneView::kU || 
        subshower->GetPlaneView()==PlaneView::kX)) 
      tpos_vtx = subshower->GetVtxU();
    if((subshower->GetPlaneView()==PlaneView::kV || 
        subshower->GetPlaneView()==PlaneView::kY)) {
      tpos_vtx = subshower->GetVtxV();
    }
    Double_t slope = subshower->GetSlope();
    Int_t beg_plane = subshower->GetBegPlane();
    Int_t end_plane = subshower->GetEndPlane();

    //for calculating asymmetry:
    Double_t n_upp = 0;
    Double_t n_low = 0;
    Double_t ph_upp = 0;
    Double_t ph_low = 0;
    Double_t n_plane_upp[500] = {0};
    Double_t n_plane_low[500] = {0};
    Double_t ph_plane_upp[500] = {0};
    Double_t ph_plane_low[500] = {0};

    Int_t firstPlane = 500;
    Int_t lastPlane = 0;

    //calculate number of strips above and below axis of subshower
    for(int i=0;i<totStp;i++){
      CandStripHandle *stp = dynamic_cast<CandStripHandle*>(newSubShower->At(i));
      if(firstPlane>stp->GetPlane()) firstPlane = stp->GetPlane();
      if(lastPlane<stp->GetPlane()) lastPlane = stp->GetPlane();
      if(stp->GetTPos() > (slope*(stp->GetZPos()-z_vtx) + tpos_vtx)) {
        n_plane_upp[stp->GetPlane()] += 1;
        ph_plane_upp[stp->GetPlane()] += stp->GetCharge(CalDigitType::kSigCorr);
        if(stp->GetPlane()>=beg_plane && stp->GetPlane()<=end_plane) {
          n_upp += 1;
          ph_upp += stp->GetCharge(CalDigitType::kSigCorr);
        }
      }
      else {
        n_plane_low[stp->GetPlane()] += 1;
        ph_plane_low[stp->GetPlane()] += stp->GetCharge(CalDigitType::kSigCorr);
        if(stp->GetPlane()>=beg_plane && stp->GetPlane()<=end_plane) {  
          n_low += 1;
          ph_low += stp->GetCharge(CalDigitType::kSigCorr);
        }
      }
    }

    //if asymmetry < 95% (should be 100%!) in the plane 
    //range of the current subshower then only keep the 
    //larger part of the new subshower:
    Double_t asymmetry = 1;
    Double_t asymmetryPH = 1;
    if( n_upp + n_low > 1 ) asymmetry = TMath::Abs(n_upp-n_low)/(n_upp+n_low);
    if( ph_upp + ph_low > 1 ) asymmetryPH = TMath::Abs(ph_upp-ph_low)/(ph_upp+ph_low);
    MSG("SubShowerSR", Msg::kDebug) << "Asymmetry = " << asymmetry
                                    << " AsymmetryPH = " << asymmetryPH
                                    << " TotStp = " << totStp
                                    << " n_upp = " << n_upp
                                    << " n_low = " << n_low
                                    << " ph_upp = " << ph_upp
                                    << " ph_low = " << ph_low
                                    << endl;
    if(asymmetry < 0.95) {
      //we know now that this proto-subshower is cutting through another subshower
      //find the greatest number of contiguous planes above or below slope
      //and take those strips only

      Int_t bestAsymFirstPlane = firstPlane;
      Int_t bestAsymLastPlane = firstPlane;
      Int_t AsymFirstPlane = firstPlane;
      Int_t currentSense = 0;
      Int_t bestSense = 0;

      //loop through all planes and find longest segments
      for(int i=firstPlane;i<=lastPlane;i+=2){
        //find longest set of contiguous planes with same asym sign:
        if(n_plane_upp[i]>n_plane_low[i] || 
           ( n_plane_upp[i]==n_plane_low[i] && 
             ph_plane_upp[i]>ph_plane_low[i]) ) {
          if(currentSense==0) {
            AsymFirstPlane = i;
            currentSense = 1;
          }
          else if(currentSense>0) currentSense += 1;
          else if(currentSense<0) { //sense has changed:
            if(TMath::Abs(currentSense)>TMath::Abs(bestSense)){
              bestSense = currentSense;
              bestAsymFirstPlane = AsymFirstPlane;
              bestAsymLastPlane = i-1;
            }
            currentSense = 1;
            AsymFirstPlane = i;
          }
        }
        else if(n_plane_low[i]>n_plane_upp[i] || 
                ( n_plane_upp[i]==n_plane_low[i] && 
                  ph_plane_low[i]>ph_plane_upp[i]) ){
          if(currentSense==0) {
            AsymFirstPlane = i;
            currentSense = -1;
          }
          else if(currentSense<0) currentSense -= 1;
          else if(currentSense>0) {
            if(TMath::Abs(currentSense)>TMath::Abs(bestSense)){
              bestSense = currentSense;
              bestAsymFirstPlane = AsymFirstPlane;
              bestAsymLastPlane = i-1;
            }
            currentSense = -1;
            AsymFirstPlane = i;
          }
        }
      }

      if(TMath::Abs(currentSense)>TMath::Abs(bestSense)){
        bestSense = currentSense;
        bestAsymFirstPlane = AsymFirstPlane;
        bestAsymLastPlane = lastPlane;
      }

      for(int i=0;i<totStp;i++){
        CandStripHandle *stp = dynamic_cast<CandStripHandle*>(newSubShower->At(i));
        if(stp->GetPlane()<bestAsymFirstPlane || stp->GetPlane()>bestAsymLastPlane){
          newSubShower->RemoveAt(i);
        }
      }
      newSubShower->Compress();
      if(newSubShower->GetLast()+1<totStp) {
        return true;
      }
    }
  }
  return false;
}

//******************************************************************
Bool_t AlgSubShowerSRList::TimeTest(TObjArray *newSubShower,TObjArray *allStrips)
{
  //Going to check that proto-subshower does not have hits with extreme times
  MSG("SubShowerSR", Msg::kDebug) << "Running TestTime()" << endl;

  //if no strips return false. Empty subshower will be caught later.
  Int_t totStp = newSubShower->GetLast()+1;
  if(totStp<=1) return false;
  
  //Look for >50ns gaps in time structure:
  std::vector<Double_t> stripTime(totStp,0.0);
  std::vector<Float_t> stripCharge(totStp,0.0);
  for(int i=0;i<totStp;i++){
    CandStripHandle *stp = dynamic_cast<CandStripHandle*>(newSubShower->At(i));
    Int_t nEarlier = 0;
    Int_t nRepeats = 0;
    for(int j=0;j<totStp;j++){
      if(i==j) continue;
      CandStripHandle *stp2 = dynamic_cast<CandStripHandle*>(newSubShower->At(j));    
      if(stp->GetTime()>stp2->GetTime()) nEarlier+=1;
      else if(stp->GetTime()==stp2->GetTime()) {
        if(stp->GetPlane()>stp2->GetPlane()) nEarlier+=1;
        else if(stp->GetPlane()==stp2->GetPlane()) {
          if(stp->GetStrip()>stp2->GetStrip()) nEarlier+=1;
          else if(stp->GetStrip()==stp2->GetStrip()) {
            if(stp->GetPlaneView()>stp2->GetPlaneView()) nEarlier+=1;
            else if(stp->GetPlaneView()==stp2->GetPlaneView()) nRepeats += 1;         
          }
        }
      }
    }
    //Fill vectors with time and charge information in time order (earliest first)
    stripTime[nEarlier] = stp->GetTime()*1e9;
    stripCharge[nEarlier] = stp->GetCharge(CalDigitType::kSigCorr);
    for(int j=1;j<=nRepeats;j++) {
      stripTime[nEarlier+j] = stp->GetTime()*1e9;
      stripCharge[nEarlier+j] = stp->GetCharge(CalDigitType::kSigCorr);
    }
  }

  //Loop through time vector to find gaps greater than 50ns
  //count the number of gaps
  Int_t nGaps = 0;
  for(int i=0;i<totStp-1;i++){
    if(stripTime[i+1] - stripTime[i] > fCleanUpTimeCut) {
      nGaps+=1;
    }
  }
  //if no gaps, return now
  if(nGaps==0) return false;
  //otherwise find the largest chunk of charge between the gaps
  std::vector<Float_t> totalCharge(nGaps+1,0.0);
  nGaps = 0;
  for(int i=0;i<totStp-1;i++) {
    totalCharge[nGaps] += stripCharge[i];
    if(stripTime[i+1] - stripTime[i] > fCleanUpTimeCut) {
      nGaps+=1;
    }
  }
  //don't forget to add in the last strip!
  totalCharge[nGaps] += stripCharge[totStp-1];

  //find the biggest chunk:
  Int_t biggestChunk = -1;
  Double_t biggestCharge = 0;
  for(int i=0;i<nGaps+1;i++) {
    if(totalCharge[i] > biggestCharge){
      biggestChunk = i;
      biggestCharge = totalCharge[i];
    }
  }

  //now go back to time vector to get the min, max times to use
  Double_t minTime = -1;
  Double_t maxTime = 0;
  nGaps = 0;
  for(int i=0;i<totStp-1;i++){
    if(nGaps>biggestChunk) break;
    if(nGaps == biggestChunk) {
      if(minTime<0) minTime = stripTime[i];
      maxTime = stripTime[i];
    }
    if(stripTime[i+1] - stripTime[i] > fCleanUpTimeCut) {
      nGaps+=1;
    }
  }

  //loop through strip list again and remove strips outside of time range
  for(int i=0;i<totStp;i++){
    CandStripHandle *stp = dynamic_cast<CandStripHandle*>(newSubShower->At(i));
    if(stp->GetTime()*1e9<minTime - fCleanUpTimeCut || 
       stp->GetTime()*1e9>maxTime + fCleanUpTimeCut) {
      newSubShower->RemoveAt(i);
      if(allStrips) allStrips->Add(stp);
      MSG("SubShowerSR",Msg::kDebug) << "Removing Strip From SubShower" << endl;
    }
  }
  newSubShower->Compress();

  if(newSubShower->GetLast()+1<totStp) return true;
  return false;

}

//******************************************************************
Bool_t AlgSubShowerSRList::CleanUp(TObjArray *allStrips,CandHandle &ch,
                                   Int_t det,const CandSliceHandle *cslice)
{
  MSG("SubShowerSR", Msg::kDebug) << "Running CleanUp()" << endl;

  Int_t totUnassigned = allStrips->GetLast()+1;
  if(totUnassigned==0) return false;

  Calibrator& calibrator=Calibrator::Instance();

  assert(ch.InheritsFrom("CandSubShowerSRListHandle"));
  CandSubShowerSRListHandle &csslh = 
    dynamic_cast<CandSubShowerSRListHandle &>(ch);

  //Get planeview for this cleanup
  TObject *firstObj = allStrips->At(0);
  assert(firstObj->InheritsFrom("CandStripHandle"));
  CandStripHandle *firstStp = dynamic_cast<CandStripHandle*>(firstObj);
  PlaneView::PlaneView_t pv = firstStp->GetPlaneView();
  MSG("SubShowerSR", Msg::kDebug) << "PlaneView = " << pv << endl;

  //Any changes made?
  Bool_t anyChanges = false;
  //once strips have been assigned do not try to assign again:
  std::vector<Bool_t> useStrip(totUnassigned,0);
  std::vector<Bool_t> assignedStrip(totUnassigned,0);
  
  for(int i=0;i<totUnassigned;i++) {
    assignedStrip[i] = false;
    //if cslice defined, check that strip is part of slice:
    if(cslice){
      useStrip[i] = false;
      TObject *tobj = allStrips->At(i);
      CandStripHandle *unassigned = dynamic_cast<CandStripHandle*>(tobj);
      CandStripHandleItr stpSliceItr(cslice->GetDaughterIterator());
      while(CandStripHandle *stp = stpSliceItr()) {
        if(*stp==*unassigned) {
          useStrip[i] = true;
          break;
        }
      }
    }
    else useStrip[i] = true;
  }
  MSG("SubShowerSR", Msg::kDebug) << "Checking " << totUnassigned 
                                  << " strips" << endl;

  //going to check if any isolated strips are nearby within +/- 2 strip
  for (Int_t i=0; i<totUnassigned; i++) {
    
    MSG("SubShowerSR", Msg::kVerbose) << "Checking unassigned strip " 
                                      << i << endl;

    TObject *tobj = allStrips->At(i);
    if(!useStrip[i]) continue;
    MSG("SubShowerSR", Msg::kVerbose) << "Can still use unassigned strip " 
                                      << i << endl;

    CandStripHandle *unassigned = dynamic_cast<CandStripHandle*>(tobj);
    std::vector<Float_t> StripDistance(csslh.GetNDaughters(),0.0);
    Float_t theCharge = 
      calibrator.GetMIP(unassigned->GetCharge(CalDigitType::kSigCorr));

    //loop through existing subshowers
    Int_t ssCount = 0;
    CandSubShowerSRHandleItr subshowerItr(csslh.GetDaughterIterator());
    while (CandSubShowerSRHandle *subshower = subshowerItr()) {       

      StripDistance[ssCount] = 99999; //set to 999999 as a flag

      //increment early because of continue statements, later use ssCount-1
      ssCount++; 

      //if a slice is given in argument list
      //demand that slices are the same
      if(cslice){
        if(*subshower->GetCandSlice()!=*cslice) continue;
      }

      //demand that any hit added to subshower is within +/- 50 ns
      if(TMath::Abs(subshower->GetAveTime() - 
                    unassigned->GetTime())*1.0e9 > fCleanUpTimeCut) continue;
      
      //only look at the right planeview subshowers
      if(subshower->GetPlaneView()!=pv) continue; 
      
      MSG("SubShowerSR", Msg::kVerbose) << "SubShower " << ssCount-1 
                                        << " in the same planeview" << endl;

      // if plane of strip is more than fCleanUpPlaneCut planes away
      // don't try to add it to this subshower 
      // (also beware of SM boundary in FarDet)
      Int_t thePlane = unassigned->GetPlane();
      if(thePlane<subshower->GetBegPlane()-fCleanUpPlaneCut || 
         thePlane>subshower->GetEndPlane()+fCleanUpPlaneCut || 
         (det==1 && 
          ( (subshower->GetBegPlane()>249 && thePlane<249) || 
            (subshower->GetEndPlane()<249 && thePlane>249) ))) continue;

      MSG("SubShowerSR", Msg::kVerbose) << "SubShower " << ssCount-1 
                                        << " in the same plane range" << endl;
 
      //loop through strips in subshower
      //check for proximity to other strips in subshower
      //(on a plane-by-plane basis)
      Bool_t checkProximity = false;
      Float_t theStrip_tpos = unassigned->GetTPos();
      CandStripHandleItr stripItr(subshower->GetDaughterIterator());
      while(CandStripHandle *stp = stripItr()){
        //check we're on the same plane 
        //or on the begin plane if the unassigned plane < begin plane
        //or on the end plane if the unassigned plane > end plane
        if((stp->GetPlane() == thePlane) ||
           (thePlane < subshower->GetBegPlane() && 
            stp->GetPlane() == subshower->GetBegPlane()) ||
            (thePlane > subshower->GetEndPlane() && 
            stp->GetPlane() == subshower->GetEndPlane())) {
          //allow +/-fCleanUpStripCut strips:
          if(TMath::Abs(theStrip_tpos-stp->GetTPos()) <=
             fCleanUpStripCut*STRIPWIDTHINMETRES) { 
            MSG("SubShowerSR",Msg::kVerbose) 
              << "TMath::Abs(theStrip_tpos-stp->GetTPos())=" 
              << TMath::Abs(theStrip_tpos-stp->GetTPos()) 
              << " stp->Plane()=" << stp->GetPlane()
              << " stp->Strip()=" << stp->GetStrip()
              << endl;
            checkProximity = true;
          }
        }
      }
      if(!checkProximity) continue;  //not close enough to this subshower

      MSG("SubShowerSR", Msg::kVerbose) 
        << "SubShower " << ssCount-1 
        << " has strips in the same tpos range" 
        << endl;

      //calculate pH weighted distance from cluster centre
      //get subshower vertex and angle:
      Double_t tposVtx = 0;
      if(pv==PlaneView::kU || 
         pv==PlaneView::kX) tposVtx = subshower->GetVtxU();
      else if(pv==PlaneView::kV || 
              pv==PlaneView::kY) tposVtx = subshower->GetVtxV();      
      Double_t zposVtx = subshower->GetVtxZ();
      Float_t slope = subshower->GetSlope(); 
      Float_t avgdev = subshower->GetAvgDev();
      //get strip position and charge:
      Float_t tpos = unassigned->GetTPos();
      Float_t zpos = unassigned->GetZPos();
      if(avgdev!=0.) { 
        StripDistance[ssCount-1] = 
          (TMath::Abs(tpos - (tposVtx + slope*(zpos - zposVtx))) * 
           TMath::Cos(TMath::ATan(slope)))*theCharge/avgdev; 
      }
      else {
        MSG("SubShowerSR",Msg::kVerbose) << "SubShower " << ssCount-1 
                                         << " has AvgDev = 0!" << endl;
        StripDistance[ssCount-1] = 99998;
      }

      MSG("SubShowerSR", Msg::kVerbose) 
        << "Unassigned strip " << i 
        << " and SubShower " << ssCount-1 
        << " are a good match with dev(strip)/avgdev(subshower) = " 
        << StripDistance[ssCount-1] << endl;
    }

    //make decision here about closest subshower
    Float_t minSD = 99999;
    Int_t bestSS = -1;
    for(int j=0;j<csslh.GetNDaughters();j++){
      if(StripDistance[j]<minSD) {
        minSD = StripDistance[j];
        bestSS = j;
      }
    }
    if(bestSS!=-1){
      MSG("SubShowerSR", Msg::kVerbose) << "Best match for strip " << i 
                                        << ": SubShower " 
                                        << bestSS << " Adding plane: " 
                                        << unassigned->GetPlane() 
                                        << " strip: " 
                                        << unassigned->GetStrip() << endl; 
      
      CandSubShowerSRHandleItr subshowerItr(csslh.GetDaughterIterator());
      ssCount = 0;
      while(CandSubShowerSRHandle *cssh = subshowerItr()){
        if(ssCount==bestSS){
          MSG("SubShowerSR",Msg::kVerbose) << "Before: cssh->GetNStrip() = " 
                                           <<  cssh->GetNStrip() 
                                           << " cssh->GetEnergy() = " 
                                           << cssh->GetEnergy() << endl; 
                                         
          cssh->AddDaughterLink(*unassigned);
          Double_t mip2gev = cssh->GetEnergy() / 
            calibrator.GetMIP(cssh->GetCharge(CalStripType::kSigCorr));
          cssh->SetEnergy(cssh->GetEnergy()+mip2gev*theCharge);
          
          MSG("SubShowerSR",Msg::kVerbose) << "After: cssh->GetNStrip() = " 
                                           << cssh->GetNStrip() 
                                           << " cssh->GetEnergy() = " 
                                           << cssh->GetEnergy() << endl;
          
          assignedStrip[i] = true;
          anyChanges = true;
          break;
        }
        ssCount+=1;
      }
    }
    else {
      MSG("SubShowerSR", Msg::kVerbose) 
        << "No matching subshower for unassigned strip " 
        << i << endl;
    }
  }
  
  //remove assigned strips
  Int_t assigned = 0;
  for(int i=0;i<totUnassigned;i++) {
    if(assignedStrip[i]) {
      allStrips->RemoveAt(i);
      assigned+=1;
    }  
  }
  allStrips->Compress();
  
  if(anyChanges) {
    MSG("SubShowerSR", Msg::kDebug) 
      << assigned << "/" << totUnassigned
      << " strips assigned to SubShowers" << endl;
  }
  else {
    MSG("SubShowerSR", Msg::kDebug) 
      << "No Changes made to SubShowers" << endl;
  }
  return anyChanges;
  
}

//******************************************************************
Bool_t AlgSubShowerSRList::FormHalo(TObjArray *allStrips,CandHandle &ch,
                                    CandContext &cxx,AlgHandle &ah,
                                    const CandSliceListHandle *slicelist,
                                    Int_t det)
{

  MSG("SubShowerSR", Msg::kDebug) << "Running FormHalo()" << endl;
  
  assert(ch.InheritsFrom("CandSubShowerSRListHandle"));

  Int_t totUnassigned = allStrips->GetLast()+1;
  if(totUnassigned==0) return false;
  
  MSG("SubShowerSR", Msg::kDebug) << "Forming slice-wise, " 
                                  << "plane-wise Halo clusters from " 
                                  << totUnassigned 
                                  << " unassigned strips" << endl;
  
  //loop over slices
  Int_t nslice = 0;
  CandSliceHandleItr sliceItr(slicelist->GetDaughterIterator());
  while (CandSliceHandle *slice = sliceItr()) {
    MSG("SubShowerSR",Msg::kDebug) << "Slice " << nslice << endl;
    
    //get valid plane ranges for this slice:
    CandSubShowerSRListHandle &subshowerlist = 
      dynamic_cast<CandSubShowerSRListHandle&>(ch);
    CandSubShowerSRHandleItr subshowerItr(subshowerlist.GetDaughterIterator());
    Int_t plane[500] = {0};
    Int_t firstPlane = 500;
    Int_t lastPlane = 0;
    while (CandSubShowerSRHandle *subshower = subshowerItr()) {
      if (*subshower->GetCandSlice()==*slice) { 
        for(int i=subshower->GetBegPlane();i<=subshower->GetEndPlane();i++){
          plane[i] = 1;
          if(i<firstPlane) firstPlane = i;
          if(i>lastPlane) lastPlane = i;
        }
      }
    }

    //loop over valid plane ranges and look for gaps > fHaloMaxPlaneGap
    while(firstPlane<=lastPlane){
      Int_t begPlane = 500;
      Int_t endPlane = 0;
      Int_t nGaps = 0;
      for(int i=firstPlane;i<=lastPlane;i++){
        if(plane[i]==0) {
          if(begPlane!=500) nGaps+=1;
        }
        else {
          if(begPlane==500) begPlane = i;
          endPlane = i;
        }
        plane[i] = 0;
        if(nGaps>fHaloMaxPlaneGap) {
          break;
        }
      }
      firstPlane = endPlane+1;

      //use plane range between begPlane and endPlane to form subshowers
      //loop over views:
      for(Int_t pv=0;pv<2;pv++){

        MSG("SubShowerSR",Msg::kDebug) << "View : " << pv << endl;

        //first loop through subshower list again and histogram tposvtx
        //in order to look for peaks
        vtxHist->Reset();
        vtxHist->SetBins(fHaloNPeakFindingBins+2,
                 -4. - 8./Float_t(fHaloNPeakFindingBins),
                 +4. + 8./Float_t(fHaloNPeakFindingBins));

        subshowerItr.ResetFirst();
        while (CandSubShowerSRHandle *subshower = subshowerItr()) {
          if (*subshower->GetCandSlice()==*slice && 
              subshower->GetClusterID()!=ClusterType::kXTalk) {
            if(pv==0 && 
               (subshower->GetPlaneView()==PlaneView::kV || 
                subshower->GetPlaneView()==PlaneView::kY)) continue;
            else if(pv==1 && 
                    (subshower->GetPlaneView()==PlaneView::kU || 
                     subshower->GetPlaneView()==PlaneView::kX)) continue;
            if(subshower->GetBegPlane()>=begPlane && 
               subshower->GetBegPlane()<=endPlane){
              if(pv==0) {
                MSG("SubShowerSR",Msg::kDebug) 
                  << "Filling vtxHist with Subshower "
                  << " with tposvtv = " << subshower->GetVtxU() 
                  << " and energy = " << subshower->GetEnergy() 
                  << " and slope = " << subshower->GetSlope()
                  << " and AvgDev = " << subshower->GetAvgDev()
                  << endl;
                CandStripHandleItr stripssItr(subshower->GetDaughterIterator());
                while (CandStripHandle *stp = stripssItr()) {
                  //fill with max charge for this tpos
                  Float_t val = stp->GetCharge(CalDigitType::kSigCorr) - 
                    vtxHist->GetBinContent(vtxHist->FindBin(stp->GetTPos()));
                  if(val>0) vtxHist->Fill(stp->GetTPos(),val);
                }
              }
              else if(pv==1) {
                MSG("SubShowerSR",Msg::kDebug) 
                  << "Filling vtxHist with Subshower "
                  << " with tposvtv = " << subshower->GetVtxV() 
                  << " and energy = " << subshower->GetEnergy() 
                  << " and slope = " << subshower->GetSlope()
                  << " and AvgDev = " << subshower->GetAvgDev()
                  << endl;
                CandStripHandleItr stripssItr(subshower->GetDaughterIterator());
                while (CandStripHandle *stp = stripssItr()) {
                  //fill with max charge for this tpos
                  Float_t val = stp->GetCharge(CalDigitType::kSigCorr) - 
                    vtxHist->GetBinContent(vtxHist->FindBin(stp->GetTPos()));
                  if(val>0) vtxHist->Fill(stp->GetTPos(),val);            
                }
              }
            }
          }
        }

        if(false){
          TCanvas *can = new TCanvas();
          can->cd(1);
          vtxHist->Draw();
          char nomus[256];
          sprintf(nomus,"peakPlot_FH_view%i_slice%i.eps",pv,nslice);
          can->Print(nomus);
          delete can;
        }
        
        //width of shower in units of number of bins
        //assume shower typically has a width of ~4 strips
        Double_t sigma = (4.*0.0417/8.)*(fHaloNPeakFindingBins);
        //don't let width be <=1 bin to prevent
        //against clipping error messages
        if(sigma<=1.0) sigma = 1.01;

        Int_t nPeaks = 0;
        Float_t peakPos[10] = {};
        Float_t peakVal[10] = {};
        Float_t tmp_peakPos = vtxHist->GetBinCenter(vtxHist->GetMaximumBin());
        Float_t tmp_peakVal = vtxHist->GetMaximum();
        for(int bins = 1;bins<=vtxHist->GetNbinsX();bins++){
          if(vtxHist->GetBinContent(bins)>0) nPeaks++;
        }
        if(nPeaks==1) {
          peakPos[0] = tmp_peakPos;
          peakVal[0] = tmp_peakVal;
        }
        else if(nPeaks>0) {
          TSpectrum spec(10);
          Float_t *source = new Float_t[vtxHist->GetNbinsX()];
          Float_t *dest = new Float_t[vtxHist->GetNbinsX()];
          for(int h_loop=0;h_loop<vtxHist->GetNbinsX();h_loop++){
            source[h_loop] = vtxHist->GetBinContent(h_loop+1);
            dest[h_loop] = 0;
          }
#if ROOT_VERSION_CODE >= ROOT_VERSION(5,9,1)
          nPeaks = spec.SearchHighRes(source,dest,
                                      vtxHist->GetNbinsX(),
                                      sigma,5,0,1,0,1);
#else
          nPeaks = spec.Search1HighRes(source,dest,
                                      vtxHist->GetNbinsX(),
                                      sigma,5,0,1,0,1);
#endif
          Float_t *posX = spec.GetPositionX();
          for(Int_t p_loop=0;p_loop<nPeaks;p_loop++){
            Int_t bin = 1+Int_t(posX[p_loop]+0.5);
            peakPos[p_loop] = vtxHist->GetBinCenter(bin);
            peakVal[p_loop] = vtxHist->GetBinContent(bin);
          }
          delete [] source;
          delete [] dest;
          //spec.Search(vtxHist,sigma);
          //nPeaks = spec.GetNPeaks();
        }

        //if no peaks are found
        //then there's no reconstructed subshower 
        //in this slice + plane range + view
        if(nPeaks==0) {
          MSG("SubShowerSR",Msg::kDebug) << "No peaks found" << endl;
          //so take everything and put it in a subshower
          CandStripHandleItr stripslcItr(slice->GetDaughterIterator());
          TObjArray newSubShower;
          while (CandStripHandle *stp = stripslcItr()) {
            if(pv==0 && (stp->GetPlaneView()==PlaneView::kV || 
                         stp->GetPlaneView()==PlaneView::kY)) continue;
            else if(pv==1 && (stp->GetPlaneView()==PlaneView::kU || 
                              stp->GetPlaneView()==PlaneView::kX)) continue;
            //is this strip in the subshower plane range?
            //allow +/- 2 planes in this window
            if(stp->GetPlane()>endPlane+2 || 
               stp->GetPlane()<begPlane-2) continue;
            //loop over unassigned strips, and see if it's there
            for (Int_t i=0; i<allStrips->GetLast()+1; i++) {
              TObject *tobj = allStrips->At(i);
              CandStripHandle *unassigned = 
                dynamic_cast<CandStripHandle*>(tobj);
              if(*unassigned == *stp) {
                newSubShower.Add(unassigned);
                allStrips->RemoveAt(i);
              }
            }
            allStrips->Compress();
          }
          //make a new subshower with the unassigned strips in this slice
          if(newSubShower.GetLast()+1>0){
            TimeTest(&newSubShower,allStrips);
            cxx.SetDataIn(&newSubShower);
            MSG("SubShowerSR",Msg::kDebug) << "forming SubShowerSR\n";
            CandSubShowerSRHandle subshowersrhandle = 
              CandSubShowerSR::MakeCandidate(ah,cxx);
            // Set the strip PE cut - this gets passed through the subshower
            // to the shower, so that it is accessible from the truthhelper
            // when it fills the purity and completeness variables
            subshowersrhandle.SetMinStripPE(fMinStripPE);
            
            if(subshowersrhandle.GetNStrip()!=0) {
              MSG("SubShowerSR",Msg::kDebug) << "Formed SubShowerSR " 
                                             << ch.GetNDaughters()
                                             << " with " 
                                             << subshowersrhandle.GetNStrip()
                                             << " strips" << endl;
              subshowersrhandle.SetClusterID(ClusterType::kHalo);
              subshowersrhandle.SetCandSlice(slice);
              ch.AddDaughterLink(subshowersrhandle);          
            }
          }
          continue;
        }
      
        //have one or more peak => maybe one or 
        //more shower in this plane range
        //first find acceptable limits in tpos for a halo
        //subshower to be formed:
        std::vector<Float_t> peakPosLow(nPeaks,0.0);
        std::vector<Float_t> peakPosUpp(nPeaks,0.0);
        MSG("SubShowerSR",Msg::kDebug) << "Number of Peaks found = "
                                       << nPeaks 
                                       << " with positions, heights "
                                       << "and bounds: "  << endl;
        for (int i=0;i<nPeaks;i++){
          if(i==0) peakPosLow[i] = -4;
          else peakPosLow[i] = peakPosUpp[i-1];
          if(i==nPeaks-1) peakPosUpp[i] = 4;
          else {
            if(peakVal[i]+peakVal[i+1] == 0) {
              MSG("SubShowerSR",Msg::kError) << "peakVal[" << i << "]=" << peakVal[i] 
                                             << " peakVal[" << i+1 << "]=" << peakVal[i+1]
                                             << " - Setting nPeaks=1!" << endl;
              nPeaks=1;
              peakPosUpp[0] = 4;
              break;
            }
            peakPosUpp[i] = peakPos[i] + 
              (peakVal[i])*(peakPos[i+1] - peakPos[i]) / 
              ((peakVal[i]) + (peakVal[i+1]));
          }
          MSG("SubShowerSR",Msg::kDebug) << peakPos[i] << " " 
                                         << peakVal[i] << " ["
                                         << peakPosLow[i] << "," 
                                         << peakPosUpp[i] << "]" 
                                         << endl;
        }

        //loop through peaks and assign strips to form halo subshowers
        //according to acceptable tpos and plane ranges
        for(int i=0;i<nPeaks;i++){
          //loop through slice strips
          CandStripHandleItr stripslcItr(slice->GetDaughterIterator());
          TObjArray newSubShower;
        
          //loop over strips in slices in this view:
          // (have to loop through slice strips and not just
          // unassigned strips because you cannot get slice number
          // from strip handle)
          while (CandStripHandle *stp = stripslcItr()) {
          
            //only look at one view at a time:
            if(pv==0 && (stp->GetPlaneView()==PlaneView::kV || 
                         stp->GetPlaneView()==PlaneView::kY)) continue;
            else if(pv==1 && (stp->GetPlaneView()==PlaneView::kU || 
                              stp->GetPlaneView()==PlaneView::kX)) continue;
          
            //is this strip in the subshower plane range?
            //allow +/- 2 planes in this window
            if(stp->GetPlane()>endPlane+2 || 
               stp->GetPlane()<begPlane-2) continue;
            
            //is this strip in the tpos range
            if(stp->GetTPos()>peakPosUpp[i] ||
               stp->GetTPos()<=peakPosLow[i]) continue;

            //loop over unassigned strips, and see if it's there
            for (Int_t i=0; i<allStrips->GetLast()+1; i++) {
              TObject *tobj = allStrips->At(i);
              CandStripHandle *unassigned = 
                dynamic_cast<CandStripHandle*>(tobj);
              if(*unassigned == *stp) {
                newSubShower.Add(unassigned);
                allStrips->RemoveAt(i);
              }
            }
            allStrips->Compress();
          }
          
          //make a new subshower with the unassigned strips in this slice
          if(newSubShower.GetLast()+1>0){  
            TimeTest(&newSubShower,allStrips);
            cxx.SetDataIn(&newSubShower);
            MSG("SubShowerSR",Msg::kDebug) << "forming SubShowerSR\n";
            CandSubShowerSRHandle subshowersrhandle = 
              CandSubShowerSR::MakeCandidate(ah,cxx);
            // Set the strip PE cut - this gets passed through the subshower
            // to the shower, so that it is accessible from the truthhelper
            // when it fills the purity and completeness variables
            subshowersrhandle.SetMinStripPE(fMinStripPE);
            
            if(subshowersrhandle.GetNStrip()!=0) {
              MSG("SubShowerSR",Msg::kDebug) << "Formed SubShowerSR with " 
                                             << subshowersrhandle.GetNStrip()
                                             << " strips" << endl;
              subshowersrhandle.SetClusterID(ClusterType::kHalo);
              subshowersrhandle.SetCandSlice(slice);
              ch.AddDaughterLink(subshowersrhandle); 
            }
          }
        }
      }
    }
    //slice-wise clean up
    Int_t numberOfCleanUps = 0;
    while(CleanUp(allStrips,ch,det,slice)) {
      numberOfCleanUps+=1;
      //prevent against infinite loop! (which should NEVER happen)
      if(numberOfCleanUps>1000) break;
    }
    nslice++;
  }
  //snarl-wide clean up
  Int_t numberOfCleanUps = 0;
  while(CleanUp(allStrips,ch,det)) {
    numberOfCleanUps+=1;
    //prevent against infinite loop! (which should NEVER happen)
    if(numberOfCleanUps>1000) break;
  }
  return true;
}

//******************************************************************
std::vector<Double_t> AlgSubShowerSRList::BestHough(std::vector<window> winVec,
                                                    Bool_t makeEPS)
{

  std::vector<window>::iterator winIterBeg = winVec.begin();
  std::vector<window>::iterator winIterEnd = winVec.end();
  std::vector<Double_t> MCV(14,0.0);
  for(int i=0;i<14;i++) MCV[i] = 0.0;

  //get extrema of windows:
  Int_t MaxPlane = 0;
  Int_t MinPlane = 500;
  Double_t MaxStrip = -5.;
  Double_t MinStrip = 5.;
  //total PH:
  Double_t totPH = 0;
  //PH of window with max PH:
  Double_t maxWinPH = 0;
  //plane with maximum PH, and value of max PH:
  Int_t maxPlanePos = -1;
  Double_t maxPlanePosPH = 0;
  //plane by plane sum PH:
  Double_t PlanePH[500] = {0};
  while(winIterBeg!=winIterEnd){
    if(MaxPlane<winIterBeg->plane) MaxPlane = winIterBeg->plane;
    if(MinPlane>winIterBeg->plane) MinPlane = winIterBeg->plane;
    if(MaxStrip<winIterBeg->upper_tpos) MaxStrip = winIterBeg->upper_tpos;
    if(MinStrip>winIterBeg->lower_tpos) MinStrip = winIterBeg->lower_tpos;
    if(maxWinPH<winIterBeg->ph) maxWinPH = winIterBeg->ph;
    totPH+=winIterBeg->ph;
    PlanePH[winIterBeg->plane] += winIterBeg->ph;
    if(PlanePH[winIterBeg->plane]>maxPlanePosPH) {
      maxPlanePos = winIterBeg->plane;
      maxPlanePosPH = PlanePH[maxPlanePos];
    }
    winIterBeg++;
  }
  
  //window tpos position with max PH, and value of max PH
  Double_t maxWindowPos = 0.;
  Double_t maxWindowPH = 0.;
  winIterBeg = winVec.begin();
  while(winIterBeg!=winIterEnd){
    if(winIterBeg->plane==maxPlanePos){
      if(winIterBeg->ph>maxWindowPH){
        maxWindowPos = winIterBeg->phpos;
        maxWindowPH = winIterBeg->ph;
      }
    }
    winIterBeg++;
  }

  if(winVec.size()<=1 || MinPlane==MaxPlane){
    MCV[0] = 0;
    MCV[1] = maxWindowPos;
    MCV[2] = maxPlanePos;
    MCV[3] = 0;
    MCV[4] = 0;
    MCV[5] = 0;
    MCV[6] = maxWindowPos;
    MCV[7] = maxPlanePos;
    MCV[8] = 0;
    MCV[9] = 0;
    MCV[10] = 0;
    MCV[11] = maxWindowPos;
    MCV[12] = 0;
    MCV[13] = maxWindowPos;
    return MCV;
  }

  Double_t TPosWin = MaxStrip-MinStrip+STRIPWIDTHINMETRES;

  MSG("SubShowerSR",Msg::kDebug) << "Hough extrema: plane " << MinPlane 
                                 << " to " << MaxPlane
                                 << " strip tpos " << MinStrip 
                                 << " to " << MaxStrip << endl;
  
  //define "x=0" to be at MinPlane-2, 
  //so y-intercept is strip value at MinPlane-2
  //gradient is in units of Delta(tpos) per plane

  houghHist->Reset();
  houghHist->SetBins(Int_t(TPosWin/STRIPWIDTHINMETRES) + 
                     MaxPlane-MinPlane+2,
                     -2.*TPosWin/Float_t(MaxPlane-MinPlane+2),
                     2.*TPosWin/Float_t(MaxPlane-MinPlane+2),
                     Int_t(2.*TPosWin/STRIPWIDTHINMETRES),
                     MinStrip-TPosWin,
                     MaxStrip+TPosWin);

  phHoughHist->Reset();
  phHoughHist->SetBins(Int_t(TPosWin/STRIPWIDTHINMETRES) + 
                       MaxPlane-MinPlane+2,
                       -2.*TPosWin/Float_t(MaxPlane-MinPlane+2),
                       2.*TPosWin/Float_t(MaxPlane-MinPlane+2),
                       Int_t(2.*TPosWin/STRIPWIDTHINMETRES),
                       MinStrip-TPosWin,
                       MaxStrip+TPosWin);
  
  winIterBeg = winVec.begin();
  while(winIterBeg!=winIterEnd){
    Int_t plane = winIterBeg->plane;
    Double_t strip = winIterBeg->phpos;
    Double_t ph = winIterBeg->ph;
    Int_t lastBin = -1;
    for(Double_t c=MinStrip-TPosWin; c<=MaxStrip+TPosWin;
        c+=STRIPWIDTHINMETRES){
      Double_t m = (strip - c)/Float_t(plane-MinPlane+2);
      MSG("SubShowerSR",Msg::kVerbose) << "Hough entries: " << m << " " 
                                       << c << " " << ph << endl;      
      if(houghHist->FindBin(m,c)!=lastBin){
        houghHist->Fill(m,c,1);
        phHoughHist->Fill(m,c,TMath::Sqrt(ph));
        lastBin = houghHist->FindBin(m,c);
      }
    }
    winIterBeg++;
  }

  Int_t maxBin = houghHist->GetMaximumBin();
  Float_t maxVal = houghHist->GetMaximum();
  Int_t xbin = maxBin%(houghHist->GetNbinsX()+2);
  Int_t ybin = maxBin/(houghHist->GetNbinsX()+2);
  Int_t phMaxBin = phHoughHist->GetMaximumBin();
  Float_t phMaxVal = phHoughHist->GetMaximum();
  Int_t phXbin = phMaxBin%(phHoughHist->GetNbinsX()+2);
  Int_t phYbin = phMaxBin/(phHoughHist->GetNbinsX()+2);

  Double_t totx = 0;
  Double_t sumx=0,sumx2=0;
  Double_t ptotx = 0;
  Double_t psumx=0,psumx2=0;

  Double_t toty=0;
  Double_t sumy=0,sumy2=0;
  Double_t ptoty = 0;
  Double_t psumy=0,psumy2=0;

  for(int xx = 1;xx<=houghHist->GetNbinsX();xx++){
    for(int yy = 1;yy<=houghHist->GetNbinsY();yy++){

      if(houghHist->GetBinContent(xx,yy)>maxVal*fBestHoughMaxPHFrac) {
        
        totx += houghHist->GetBinContent(xx,yy);
        sumx += ( houghHist->GetBinContent(xx,yy) * 
                  houghHist->GetXaxis()->GetBinCenter(xx)  );
        sumx2 += ( houghHist->GetBinContent(xx,yy) * 
                   houghHist->GetXaxis()->GetBinCenter(xx)  *
                   houghHist->GetXaxis()->GetBinCenter(xx)  );
        
        toty += houghHist->GetBinContent(xx,yy);
        sumy += ( houghHist->GetBinContent(xx,yy) * 
                  houghHist->GetYaxis()->GetBinCenter(yy)  );
        sumy2 += ( houghHist->GetBinContent(xx,yy) * 
                   houghHist->GetYaxis()->GetBinCenter(yy)  *
                   houghHist->GetYaxis()->GetBinCenter(yy)  );
      }

      if(phHoughHist->GetBinContent(xx,yy)>phMaxVal*fBestHoughMaxPHFrac) {

        ptotx += phHoughHist->GetBinContent(xx,yy);
        psumx += ( phHoughHist->GetBinContent(xx,yy) * 
                   phHoughHist->GetXaxis()->GetBinCenter(xx)  );
        psumx2 += ( phHoughHist->GetBinContent(xx,yy) * 
                    phHoughHist->GetXaxis()->GetBinCenter(xx)  *
                    phHoughHist->GetXaxis()->GetBinCenter(xx)  );
        
        ptoty += phHoughHist->GetBinContent(xx,yy);
        psumy += ( phHoughHist->GetBinContent(xx,yy) * 
                   phHoughHist->GetYaxis()->GetBinCenter(yy)  );
        psumy2 += ( phHoughHist->GetBinContent(xx,yy) * 
                    phHoughHist->GetYaxis()->GetBinCenter(yy)  *
                    phHoughHist->GetYaxis()->GetBinCenter(yy)  );
      }

    }
  }

  MCV[2] = MinPlane - 2;
  if(totx>0) {
    MCV[0] = sumx/totx;
    MCV[3] = (sumx2/totx - TMath::Power(MCV[0],2))/totx;
    if(MCV[3]>0) MCV[3] = TMath::Sqrt(MCV[3]);
    else MCV[3] = 0;
  }
  if(toty>0) {
    MCV[1] = sumy/toty;
    MCV[4] = (sumy2/toty - TMath::Power(MCV[1],2))/toty;
    if(MCV[4]>0) MCV[4] = TMath::Sqrt(MCV[4]);
    else MCV[4] = 0;
  }
  if(MCV[3]<houghHist->GetXaxis()->GetBinWidth(1))
    MCV[3] = houghHist->GetXaxis()->GetBinWidth(1);
  if(MCV[4]<houghHist->GetYaxis()->GetBinWidth(1))
    MCV[4] = houghHist->GetYaxis()->GetBinWidth(1);

  MCV[7] = MinPlane - 2;
  if(ptotx>0) {
    MCV[5] = psumx/ptotx;
    MCV[8] = (psumx2/ptotx - TMath::Power(MCV[5],2))/ptotx;
    if(MCV[8]>0) MCV[8] = TMath::Sqrt(MCV[8]);
    else MCV[8] = 0;
  }
  if(ptoty>0) {
    MCV[6] = psumy/ptoty;
    MCV[9] = (psumy2/ptoty - TMath::Power(MCV[6],2))/ptoty;
    if(MCV[9]>0) MCV[9] = TMath::Sqrt(MCV[9]);
    else MCV[9] = 0;  
  }
  if(MCV[8]<phHoughHist->GetXaxis()->GetBinWidth(1))
    MCV[8] = phHoughHist->GetXaxis()->GetBinWidth(1);  
  if(MCV[9]<phHoughHist->GetYaxis()->GetBinWidth(1))
    MCV[9] = phHoughHist->GetYaxis()->GetBinWidth(1);

  MCV[10] = houghHist->GetXaxis()->GetBinCenter(xbin);
  MCV[11] = houghHist->GetYaxis()->GetBinCenter(ybin);
  MCV[12] = phHoughHist->GetXaxis()->GetBinCenter(phXbin);
  MCV[13] = phHoughHist->GetYaxis()->GetBinCenter(phYbin);

  if(makeEPS) {
    gStyle->SetOptStat(0);
    TCanvas *can = new TCanvas();
    can->Divide(2,2);
    can->cd(1);
    houghHist->Draw("COLZ");
    can->cd(2);
    phHoughHist->Draw("COLZ");
    can->cd(3);
    houghHist->Draw("LEGO2");
    can->cd(4);
    phHoughHist->Draw("LEGO2");
    Int_t countr = 0;
    Bool_t noPlot = true;
    while(noPlot){
      char name[256];
      sprintf(name,"best_%i.root",countr);
      std::ifstream Test(name);
      if(!Test) { //if files does not exist:
        can->Print(name);
        noPlot = false;
      }
      countr+=1;
    }
    delete can;
  }

  MSG("SubShowerSR",Msg::kDebug) 
    << "Before Cuts in BestHough:"
    << "\nHough gradient, intercept, plane vertex: " 
    << MCV[0] << "+/-" << MCV[3] << " " 
    << MCV[1] << "+/-" << MCV[4] << " " << MCV[2]
    << "\nPH Hough gradient, intercept, plane vertex: " 
    << MCV[5] << "+/-" << MCV[8] << " " 
    << MCV[6] << "+/-" << MCV[9] << " " << MCV[7] 
    << "\n Hough Max Bin Gradient and Intercept Values: "
    << MCV[10] << " " << MCV[11] << " "
    << MCV[12] << " " << MCV[13] << endl;

  //check that there is a coincidence of 3 or more windows
  //and gradient is bigger than error
  if(maxVal < 3 || 
     (TMath::Abs(MCV[0])-MCV[3]<0 && 
      TMath::Abs(MCV[10])-MCV[3]<0)) {
    MCV[0] = 0;
    MCV[1] = maxWindowPos;
    MCV[2] = maxPlanePos;
    MCV[3] = 0;
    MCV[4] = 0; 
    MCV[10] = 0;
    MCV[11] = maxWindowPos;
  }
  
  if((phMaxVal - maxWindowPH)/(totPH - maxWindowPH) < 0.1 || 
     (TMath::Abs(MCV[5])-MCV[8]<0 ) &&
     (TMath::Abs(MCV[12])-MCV[8]<0 )) {
    MCV[5] = 0;
    MCV[6] = maxWindowPos;
    MCV[7] = maxPlanePos;
    MCV[8] = 0;
    MCV[9] = 0;
    MCV[12] = 0;
    MCV[13] = maxWindowPos;
  }
  return MCV;
  
}

// ----------------------------------------------------------
NavKey AlgSubShowerSRList::StripKeyFromPlane(const CandStripHandle *strip)
{
  return const_cast<CandStripHandle *>(strip)->GetPlane();
}

// ----------------------------------------------------------
void AlgSubShowerSRList::Trace(const char * /* c */) const
{
}

---