////////////////////////////////////////////////////////////////////////
// $Id: HoughTrackSR.cxx,v 1.17 2007/01/15 20:02:44 rhatcher Exp $
//
// HoughTrackSR
//
// Author:  R. Lee 2002.01.18
//
////////////////////////////////////////////////////////////////////////

#include <iostream>
#include <string>
#include <cmath>
#include "TMath.h"

#include "CandTrackSR/HoughTrackSR.h"
#include "MessageService/MsgService.h"
#include "CandTrackSR/TrackClusterSR.h"
#include "RecoBase/LinearFit.h"

ClassImp(HoughTrackSR)

//______________________________________________________________________
CVSID("$Id: HoughTrackSR.cxx,v 1.17 2007/01/15 20:02:44 rhatcher Exp $");

HoughTrackSR::HoughTrackSR(TObjArray clusterlist)
{
  fSlopeLim[0] = -20.;
  fSlopeLim[1] =  21.;
  fInterLim[0] = -10.;
  fInterLim[1] =  11.;
  fMaxBin=40;
  fPeakMin = 4;
  fPeakMinFrac = .25;
  fPeakMinFracZoom = .75;
  fMinInterBinSize = .02;
  fClusterList = new TObjArray(clusterlist);
  fFillHough = kFALSE;
}

HoughTrackSR::HoughTrackSR(const HoughTrackSR &oldHT) : TObject(oldHT)
{
  fSlopeLim[0] = oldHT.fSlopeLim[0];
  fSlopeLim[1] = oldHT.fSlopeLim[1];
  fInterLim[0] = oldHT.fInterLim[0];
  fInterLim[1] = oldHT.fInterLim[1];
  fPeakMin = oldHT.fPeakMin;
  fPeakMinFrac = oldHT.fPeakMinFrac;
  fPeakMinFracZoom = oldHT.fPeakMinFracZoom;
  fMinInterBinSize = oldHT.fMinInterBinSize;
  fClusterList = new TObjArray(*(oldHT.fClusterList));
  fFillHough = oldHT.fFillHough;

  fXZmean[0] = oldHT.fXZmean[0];
  fXZmean[1] = oldHT.fXZmean[1];
  fXZrms[0] = oldHT.fXZrms[0];
  fXZrms[1] = oldHT.fXZrms[1];
  fDSlope = oldHT.fDSlope;
  fDInter = oldHT.fDInter;
  fMaxBin = oldHT.fMaxBin;
  fNBinMax = oldHT.fNBinMax;
  fNBinCut = oldHT.fNBinCut;

  for (Int_t im=0; im<fNBin; im++) {
    for (Int_t ic=0; ic<fNBin; ic++) {
      fHough[im][ic] = oldHT.fHough[im][ic];
      fMaxHoughBin[im][ic][0] = oldHT.fMaxHoughBin[im][ic][0];
      fMaxHoughBin[im][ic][1] = oldHT.fMaxHoughBin[im][ic][1];
    }
  }

  fNPeak = oldHT.fNPeak;
  fNMaxima = oldHT.fNMaxima;

  for (Int_t i=0; i<fNMaxima; i++) {
    fMPeak.push_back(oldHT.fMPeak[i]);
    fCPeak.push_back(oldHT.fCPeak[i]);
    fIPeak.push_back(oldHT.fIPeak[i]);
    fSlopePeak.push_back(oldHT.fSlopePeak[i]);
  }

}

HoughTrackSR::HoughTrackSR()
{
  fSlopeLim[0] = -20.;
  fSlopeLim[1] =  21.;
  fInterLim[0] = -10.;
  fInterLim[1] =  11.;
  fMaxBin=40;
  fPeakMin = 4;
  fPeakMinFrac = .5;
  fPeakMinFracZoom = .75;
  fMinInterBinSize = .02;
  fFillHough = kFALSE;
}

HoughTrackSR::~HoughTrackSR()
{
  if (fClusterList) {
    delete fClusterList;
  }
}

void HoughTrackSR::SetClusterList(TObjArray clusterlist)
{
  fClusterList = new TObjArray(clusterlist);
}

void HoughTrackSR::FillHough()
{
  fFillHough = kTRUE;

  fDSlope = fSlopeLim[1]-fSlopeLim[0];
  fDInter = fInterLim[1]-fInterLim[0];
  for (int im=0; im<fMaxBin; im++) {
    for (int ic=0; ic<fMaxBin; ic++) {
      fHough[im][ic] = 0;
    }
  }
  fMPeak.erase(fMPeak.begin(),fMPeak.end());
  fCPeak.erase(fCPeak.begin(),fCPeak.end());
  fIPeak.erase(fIPeak.begin(),fIPeak.end());
  fSlopePeak.erase(fSlopePeak.begin(),fSlopePeak.end());
  for (int i=0; i<=fClusterList->GetLast(); i++) {
    TrackClusterSR *cluster = dynamic_cast<TrackClusterSR*>(fClusterList->At(i));
    if (cluster) {
      for (int im=0; im<fMaxBin; im++) {
        Double_t slope = (Double_t)(im+.5)/(Double_t)fMaxBin*fDSlope+fSlopeLim[0];
        Double_t inter = cluster->GetTPos()-fXZmean[0]-slope*(cluster->GetZPos()-fXZmean[1]);
	int ic=-1;
	if(fDInter!=0){
	  ic = (int)((inter-fInterLim[0])/fDInter*fMaxBin+.5);
	}
        if (ic>=0 && ic<fMaxBin) {
          fHough[im][ic]++;
        }
      }
    }
  }
// find maximum bin
  fNBinMax = 0;
  fNPeak = 0;
  fNMaxima = 0;
  for (int im=0; im<fMaxBin; im++) {
    for (int ic=0; ic<fMaxBin; ic++) {
      if (fHough[im][ic]>fNBinMax) {
        fNBinMax = fHough[im][ic];
      }
      Int_t nmax = fHough[im][ic];
      Int_t nmin=10000;
      Int_t immax = 0;
      Int_t icmax = 0;
      for (int im1=-3; im1<=3; im1++) {
        for (int ic1=-3; ic1<=3; ic1++) {
          if (im+im1>=0 && im+im1<fNBin && ic+ic1>=0 && ic+ic1<fMaxBin && fHough[im+im1][ic+ic1]>nmax) {
            nmax = fHough[im+im1][ic+ic1];
            immax = im1;
            icmax = ic1;
          }
	  if (im+im1>=0 && im+im1<fMaxBin && ic+ic1>=0 && ic+ic1<fMaxBin && fHough[im+im1][ic+ic1]<nmin) nmin=fHough[im+im1][ic+ic1];
	}
      }
      if((nmax-nmin)<fPeakMin){
	immax=-40;
	icmax=-40;
      }
      fMaxHoughBin[im][ic][0] = immax;
      fMaxHoughBin[im][ic][1] = icmax;
    }
  }
  /*
  MSG("HoughTrackSR", Msg::kInfo)  << " MAX -----";
  for (int ic=0; ic<fMaxBin; ic++) {
	  MSG("HoughTrackSR", Msg::kInfo)  << "---";
  }
  MSG("HoughTrackSR", Msg::kInfo)  << "\n";
  for (int im=0; im<fMaxBin; im++) {
    printf("%2d | ",im);
    for (int ic=0; ic<fMaxBin; ic++) {
      printf("%2d/%2d ",fMaxHoughBin[im][ic][0],fMaxHoughBin[im][ic][1] );
    }
    MSG("HoughTrackSR", Msg::kInfo)  << "\n";
  }
  */

// find peak(s)
  Int_t nmuon = 0;
  Int_t imaxpeak = 0;
  for (int im=0; im<fMaxBin; im++) {
    for (int ic=0; ic<fMaxBin; ic++) {
      if (fHough[im][ic]>=fPeakMin && fMaxHoughBin[im][ic][0]==0 && fMaxHoughBin[im][ic][1]==0) {
        int ifound = 0;
        for (int imaxima=0; imaxima<fNMaxima; imaxima++) {
          if (abs(fMPeak[imaxima]-im)<=3 && abs(fCPeak[imaxima]-ic)<=3) {
            if (!ifound) {
              ifound = fIPeak[imaxima];
            }
            else if (ifound!=fIPeak[imaxima]) {
              for (int jmaxima=0; jmaxima<fNMaxima; jmaxima++) {
                if (fIPeak[jmaxima]==fIPeak[imaxima]) {
                  fIPeak[jmaxima] = ifound;
                }
              }
            }
          }
        }
 
        fMPeak.push_back(im);
        fCPeak.push_back(ic);
        if (ifound) {
          fIPeak.push_back(ifound);
        }
        else {
          nmuon++;
          fIPeak.push_back(nmuon);
        }
        if (!fNMaxima || fHough[im][ic]==fNBinMax) {
          imaxpeak = fNMaxima;
        }
        fNMaxima++;
      }
    }
  }
// calculate slope in hough space
  for (int i=0; i<fNMaxima; i++) {
    int ic = fCPeak[i];
    Double_t xhough[7] = {-3.,-2.,-1.,0.,1.,2.,3.};
    Double_t yhough[7] = {0.,0.,0.,0.,0.,0.,0.};
    Double_t whough[7] = {0.,0.,0.,0.,0.,0.,0.};
    Int_t ncexist=0;
    for (int ic1=-3; ic1<=3; ic1++) {
      if (ic+ic1>=0 && ic+ic1<fMaxBin) {
        for (int im1=0; im1<fMaxBin; im1++) {
          yhough[ic1+3] += (Double_t)(im1*fHough[im1][ic+ic1]);
          whough[ic1+3] += (Double_t)(fHough[im1][ic+ic1]);
        }
      }
      if (whough[ic1+3]>0.) {
        yhough[ic1+3] /= whough[ic1+3];
        ncexist++;
      }
    }
    if (ncexist>1) {
      Double_t parm[2];
      LinearFit::Weighted(7,xhough,yhough,whough,parm);
      fSlopePeak.push_back(parm[1]);
    } else {
      fSlopePeak.push_back(99999.); // no nearby peaks, set to non physical number so that it does not match any other peaks
    }
  }



// match local maxima, and remove those below threshold
  for (int i=0; i<fNMaxima; i++) {
    for (int j=0; j<fNMaxima; j++) {
      if (i!=j) {
	int  avec = (fCPeak[j] + fCPeak[i])/2;
	int  avem = (fMPeak[j] + fMPeak[i])/2;
	Double_t peakisum = 0;
	Double_t peakjsum = 0;
	Double_t valleysum = 0;
	for (int ij=-1;ij<=1;ij++){
	  for (int ik=-1;ik<=1;ik++){
	    if(fMPeak[i]+ij>=0 && fMPeak[i]+ij<fMaxBin && fCPeak[i]+ik>=0 && fCPeak[i]+ik<fMaxBin)  peakisum += fHough[fMPeak[i]+ij][fCPeak[i]+ik];	
	    if(fMPeak[j]+ij>=0 && fMPeak[j]+ij<fMaxBin && fCPeak[j]+ik>=0 && fCPeak[j]+ik<fMaxBin)   peakjsum +=fHough[fMPeak[j]+ij][fCPeak[j]+ik];
	    if(avem+ij>=0 && avem+ij<fMaxBin && avec+ik>=0 && avec+ik<fMaxBin) valleysum += fHough[avem+ij][avec+ik];
	  }
	}
	Double_t minpeak  = min(peakisum,peakjsum);
	Double_t threshold = 2 * valleysum + fPeakMin;
	if(valleysum>10) threshold = valleysum + 2.* TMath::Sqrt(valleysum);
	//	cout << i << " " << j << " " << avem << " " << avec << " " << peakisum << " " << peakjsum << " " << valleysum << endl;
	if(minpeak<threshold){
	  //	  cout << "remove " << endl;
	  if(peakisum<peakjsum)fIPeak[i]=-1;
	  else fIPeak[j]=-1;
	}

	
      }
    }
  }

  // remove all peaks that are shorter than 20% of longest
  for (int i=0; i<fNMaxima; i++) {
    if(fHough[fMPeak[i]][fCPeak[i]]<0.2*fNBinMax)fIPeak[i]=-1;
    
  }
// require slopes to be within +-3 bins of maximum peak (this only seems to make sense for cosmics)
/*
  for (int i=0; i<fNMaxima; i++) {
    if (i!=imaxpeak && abs(fMPeak[i]-fMPeak[imaxpeak])>3) {
      fIPeak[i] = -1;
    }
  }
*/

  int muonmult[fNBin*fNBin+1];
  for (int i=0; i<=fNBin*fNBin; i++) {
    muonmult[i] = 0;
  }
  for (int i=0; i<fNMaxima; i++) {
    if (fIPeak[i]>=0) {
      muonmult[fIPeak[i]] = 1;
    }
  }
  for (int i=0; i<=fMaxBin*fMaxBin; i++) {
    fNPeak += muonmult[i];
  } 
}

Int_t HoughTrackSR::GetNPeak() const
{
  return fNPeak;
}

void HoughTrackSR::Print(Option_t* /* option */) const
{
	MSG("HoughTrackSR", Msg::kInfo)  << "MEAN " << fXZmean[0] << " " << fXZmean[1] << "\n";
	MSG("HoughTrackSR", Msg::kInfo)  << "RMS  " << fXZrms[0] << " " << fXZrms[1]  << "\n";

	MSG("HoughTrackSR", Msg::kInfo)  << "slope limits: " << fSlopeLim[0] << " / " << fSlopeLim[1] << "\n";
	MSG("HoughTrackSR", Msg::kInfo)  << "intercept limits: " << fInterLim[0] << " / " << fInterLim[1] << "\n";
	MSG("HoughTrackSR", Msg::kInfo)  << "\n";

	MSG("HoughTrackSR", Msg::kInfo)  << "   | ";
  for (int ic=0; ic<fMaxBin; ic++) {
    printf("%2d ",ic);
  }
  MSG("HoughTrackSR", Msg::kInfo)  << "\n";
  MSG("HoughTrackSR", Msg::kInfo)  << "-----";
  for (int ic=0; ic<fMaxBin; ic++) {
	  MSG("HoughTrackSR", Msg::kInfo)  << "---";
  }
  MSG("HoughTrackSR", Msg::kInfo)  << "\n";
  for (int im=0; im<fMaxBin; im++) {
    printf("%2d | ",im);
    for (int ic=0; ic<fMaxBin; ic++) {
      printf("%2d ",fHough[im][ic]);
    }
    MSG("HoughTrackSR", Msg::kInfo)  << "\n";
  }
  MSG("HoughTrackSR", Msg::kInfo)  << "\n\n";

  MSG("HoughTrackSR", Msg::kInfo)  << "# of found muons = " << fNPeak << "\n";
  MSG("HoughTrackSR", Msg::kInfo)  << "# of found local maxima = " << fNMaxima << "\n";

  for (int i=0; i<fNMaxima; i++) {
	  MSG("HoughTrackSR", Msg::kInfo)  << "  " << i << " " << fMPeak[i] << " " << fCPeak[i] << " " << fIPeak[i] << " " << fSlopePeak[i] << " " << fHough[fMPeak[i]][fCPeak[i]] << "\n";
    if (fIPeak[i]>=0) {
		MSG("HoughTrackSR", Msg::kInfo)  << "  PEAK: " << GetSlope(i) << " " << GetInter(i) << "\n";
    }
  }
  MSG("HoughTrackSR", Msg::kInfo)  << "\n\n";

}

TObjArray HoughTrackSR::Iterate()
{
  TObjArray tracklist;
  if (fFillHough) {
    Double_t SlopeLim[2] = {fSlopeLim[0],fSlopeLim[1]};
    Double_t InterLim[2] = {fInterLim[0],fInterLim[1]};
    if (fNPeak>1) {
      Int_t imuon=0;
      Double_t newslopelim[2] = { 0, 0 }, newinterlim[2] = { 0, 0 };
      for (Int_t i=0; i<fNMaxima; i++) {
        if (fIPeak[i]>=0) {
          Int_t nentry = fHough[fMPeak[i]][fCPeak[i]];
          fNBinCut = max(fPeakMin,(Int_t)(fPeakMinFracZoom*nentry));
          Int_t impeak = (Int_t)(GetSlope(i)+.5);
          Int_t icpeak = (Int_t)(GetInter(i)+.5);
          Int_t minbin[2] = {impeak-5,icpeak-5};
          Int_t maxbin[2] = {impeak+5,icpeak+5};
	  /*
          Int_t found=1;
          for (Int_t ic1=icpeak-1; ic1>=0 && found; ic1--) {
            found = 0;
            Int_t im1 = (Int_t)((ic1-icpeak)*fSlopePeak[i]+.5)+impeak;
            for (Int_t im2=-3; im2<=3; im2++) {
              if (im1+im2>=0 && im1+im2<fMaxBin && fHough[im1+im2][ic1]>=fNBinCut) {
                found = 1;
                minbin[1] = ic1;
                minbin[0] = min(minbin[0],im1+im2);
                maxbin[0] = max(maxbin[0],im1+im2);
              }
            }
          }
          found = 1;
          for (Int_t ic1=icpeak+1; ic1<fMaxBin && found; ic1++) {
            found = 0;
            Int_t im1 = (Int_t)((ic1-icpeak)*fSlopePeak[i]+.5)+impeak;
            for (Int_t im2=-3; im2<=3; im2++) {
              if (im1+im2>=0 && im1+im2<fMaxBin && fHough[im1+im2][ic1]>=fNBinCut) {
                found = 1;
                maxbin[1] = ic1;
                minbin[0] = min(minbin[0],im1+im2);
                maxbin[0] = max(maxbin[0],im1+im2);
              }
            }
          }
	  */
	  // expand intercept bin range to minimum allowed if we have specified
	  // a range that is smaller than allowed by fMinInterBinSize
          Int_t dMinInterBin = (Int_t)(fMinInterBinSize*(Double_t)(fMaxBin*fMaxBin)/fDInter);
          while (dMinInterBin<fMaxBin && maxbin[1]-minbin[1]+1<dMinInterBin) {
            minbin[1]--;
            if (minbin[1]<0) minbin[1]=0;
            if (maxbin[1]-minbin[1]<dMinInterBin) {
              maxbin[1]++;
              if (maxbin[1]>=fMaxBin) maxbin[1]=fMaxBin-1;
            }
          }

          if (!imuon) {
	    // set new limits on original track space
	    if (fMPeak[i]-minbin[0]<fMaxBin/4) minbin[0] = max(0,fMPeak[i]-fMaxBin/4);
	    if (fCPeak[i]-minbin[1]<fMaxBin/4) minbin[1] = max(0,fCPeak[i]-fMaxBin/4);
	    if (maxbin[0]-fMPeak[i]<fMaxBin/4) maxbin[0] = min(fMaxBin-1,fMPeak[i]+fMaxBin/4);
	    if (maxbin[1]-fCPeak[i]<fMaxBin/4) maxbin[1] = min(fMaxBin-1,fCPeak[i]+fMaxBin/4);
	    minbin[0] = min(minbin[0],fMPeak[i]);
	    minbin[1] = min(minbin[1],fCPeak[i]);
	    maxbin[0] = max(maxbin[0],fMPeak[i]);
	    maxbin[1] = max(maxbin[1],fCPeak[i]);

            newslopelim[0] = SlopeLim[0]+(Double_t)minbin[0]/(Double_t)fMaxBin*fDSlope;
            newslopelim[1] = SlopeLim[0]+(Double_t)maxbin[0]/(Double_t)fMaxBin*fDSlope;
            newinterlim[0] = InterLim[0]+(Double_t)minbin[1]/(Double_t)fMaxBin*fDInter;
            newinterlim[1] = InterLim[0]+(Double_t)maxbin[1]/(Double_t)fMaxBin*fDInter;
          }
          else {
	    // create new hough track space centered on peak
            HoughTrackSR *newhough = new HoughTrackSR(*fClusterList);
            newhough->SetX0(fXZmean[0]);
            newhough->SetZ0(fXZmean[1]);
            newhough->SetXRMS(fXZrms[0]);
            newhough->SetZRMS(fXZrms[1]);
            newhough->SetPeakMin(fPeakMin);
            newhough->SetPeakMinFrac(fPeakMinFrac);

            newhough->SetSlopeMin(SlopeLim[0]+(Double_t)minbin[0]/(Double_t)fMaxBin*fDSlope);
            newhough->SetSlopeMax(SlopeLim[0]+(Double_t)maxbin[0]/(Double_t)fMaxBin*fDSlope);
            newhough->SetInterceptMin(InterLim[0]+(Double_t)minbin[1]/(Double_t)fMaxBin*fDInter);
            newhough->SetInterceptMax(InterLim[0]+(Double_t)maxbin[1]/(Double_t)fMaxBin*fDInter);
            newhough->FillHough();
            tracklist.Add(newhough);
          }
          imuon++;
        }
      }

      // set new space limits and refill this track
      fSlopeLim[0] = newslopelim[0];
      fSlopeLim[1] = newslopelim[1];
      fInterLim[0] = newinterlim[0];
      fInterLim[1] = newinterlim[1];
      FillHough();
    }
    else if (fNPeak==1) {
      Int_t i;
      for (i=0; i<fNMaxima && fIPeak[i]<0; i++);

      // we now determine what bin count we use to define new boundaries of
      // slope/int space
      fNBinCut = max(fPeakMin,(Int_t)(fPeakMinFrac*fHough[fMPeak[i]][fCPeak[i]]));
      Int_t minbin[2] = {fMaxBin,fMaxBin};
      Int_t maxbin[2] = {0,0};
      for (Int_t im=0; im<fMaxBin; im++) {
        for (Int_t ic=0; ic<fMaxBin; ic++) {
          if (fHough[im][ic]>=fNBinCut) {
            minbin[0] = min(minbin[0],im);
            minbin[1] = min(minbin[1],ic);
            maxbin[0] = max(maxbin[0],im);
            maxbin[1] = max(maxbin[1],ic);
          }
        }
      }
      Int_t dinter = (Int_t)(fMinInterBinSize*(Double_t)(fMaxBin*fMaxBin)/fDInter);
      while (dinter<fMaxBin && maxbin[1]-minbin[1]+1<dinter) {
        minbin[1]--;
        if (minbin[1]<0) minbin[1]=0;
        if (maxbin[1]-minbin[1]<dinter) {
          maxbin[1]++;
          if (maxbin[1]>=fMaxBin) maxbin[1]=fMaxBin-1;
        }
      }

      // don't allow excessive scale change in one iteration
      if (fMPeak[i]-minbin[0]<fMaxBin/4) minbin[0] = max(0,fMPeak[i]-fMaxBin/4);
      if (fCPeak[i]-minbin[1]<fMaxBin/4) minbin[1] = max(0,fCPeak[i]-fMaxBin/4);
      if (maxbin[0]-fMPeak[i]<fMaxBin/4) maxbin[0] = min(fMaxBin-1,fMPeak[i]+fMaxBin/4);
      if (maxbin[1]-fCPeak[i]<fMaxBin/4) maxbin[1] = min(fMaxBin-1,fCPeak[i]+fMaxBin/4);
      minbin[0] = min(minbin[0],fMPeak[i]);
      minbin[1] = min(minbin[1],fCPeak[i]);
      maxbin[0] = max(maxbin[0],fMPeak[i]);
      maxbin[1] = max(maxbin[1],fCPeak[i]);
      fSlopeLim[1] = fSlopeLim[0]+(Double_t)maxbin[0]/(Double_t)(fMaxBin-1)*fDSlope;
      fSlopeLim[0] += (Double_t)minbin[0]/(Double_t)(fMaxBin-1)*fDSlope;
      fInterLim[1] = fInterLim[0]+(Double_t)maxbin[1]/(Double_t)(fMaxBin-1)*fDInter;
      fInterLim[0] += (Double_t)minbin[1]/(Double_t)(fMaxBin-1)*fDInter;
      FillHough();
    }
  }
  else {
    int nrms = 1;
    fInterLim[0] = min(-nrms*fXZrms[0]-fSlopeLim[1]*nrms*fXZrms[1],
                       -nrms*fXZrms[0]+fSlopeLim[0]*nrms*fXZrms[1]);
    fInterLim[1] = max(nrms*fXZrms[0]+fSlopeLim[1]*nrms*fXZrms[1],
                       nrms*fXZrms[0]-fSlopeLim[0]*nrms*fXZrms[1]);
    FillHough();
  }


  return tracklist;

}

Int_t HoughTrackSR::GetNBin() const
{
  return fMaxBin;
}

void HoughTrackSR::SetPeakMin(Int_t nbin)
{
  fPeakMin = nbin;
}

void HoughTrackSR::SetPeakMinFrac(Double_t dvar)
{
  fPeakMinFrac = dvar;
}

void HoughTrackSR::SetPeakMinFracZoom(Double_t dvar)
{
  fPeakMinFracZoom = dvar;
}

void HoughTrackSR::SetMinInterBinSize(Double_t dvar)
{
  fMinInterBinSize = dvar;
}

Double_t HoughTrackSR::GetSlopeMin() const
{
  return fSlopeLim[0];
}

Double_t HoughTrackSR::GetSlopeMax() const
{
  return fSlopeLim[1];
}

Double_t HoughTrackSR::GetInterceptMin() const
{
  return fInterLim[0];
}

Double_t HoughTrackSR::GetInterceptMax() const
{
  return fInterLim[1];
}

void HoughTrackSR::SetX0(Double_t dvar)
{
  fXZmean[0] = dvar;
}

void HoughTrackSR::SetZ0(Double_t dvar)
{
  fXZmean[1] = dvar;
}

void HoughTrackSR::SetXRMS(Double_t dvar)
{
  fXZrms[0] = dvar;
}

void HoughTrackSR::SetZRMS(Double_t dvar)
{
  fXZrms[1] = dvar;
}


void HoughTrackSR::SetSlopeMin(Double_t dvar)
{
  fSlopeLim[0] = dvar;
}

void HoughTrackSR::SetSlopeMax(Double_t dvar)
{
  fSlopeLim[1] = dvar;
}

void HoughTrackSR::SetInterceptMin(Double_t dvar)
{
  fInterLim[0] = dvar;
}

void HoughTrackSR::SetInterceptMax(Double_t dvar)
{
  fInterLim[1] = dvar;
}

Double_t HoughTrackSR::GetPeakSlope(Int_t i) const
{
  if (fNPeak==0) return -9999.;
  if (i<0 || i>=fNPeak) i = 0;
  Int_t j;
  Int_t n=0;
  for (j=0; j<fNMaxima && n<=i; j++) {
    if (fIPeak[j]>=0) n++;
  }
  j--;
  Double_t slope = fMPeak[j]*fDSlope/fMaxBin+fSlopeLim[0];
  return slope;
}

Double_t HoughTrackSR::GetPeakIntercept(Int_t i) const
{
  if (fNPeak==0) return -9999.;
  if (i<0 || i>=fNPeak) i = 0;
  Int_t j;
  Int_t n=0;
  for (j=0; j<fNMaxima && n<=i; j++) {
    if (fIPeak[j]>=0) n++;
  }
  j--;

  Double_t slope = fMPeak[j]*fDSlope/fMaxBin+fSlopeLim[0];
  Double_t inter = fCPeak[j]*fDInter/fMaxBin+fInterLim[0];

  return inter-slope*fXZmean[1]+fXZmean[0]; 
}

Int_t HoughTrackSR::GetPeakBin(Int_t i) const
{
  if (fNPeak==0) return -9999;
  if (i<0 || i>=fNPeak) i = 0;
  Int_t j;
  Int_t n=0;
  for (j=0; j<fNMaxima && n<=i; j++) {
    if (fIPeak[j]>=0) n++;
  }
  j--;
  return fHough[fMPeak[j]][fCPeak[j]];
}

Int_t HoughTrackSR::GetLargestPeakIndex() const
{
  Int_t nmaxbin = 0;
  Int_t imax = -9999;
  Int_t ipeak = 0;
  for (Int_t i=0; i<fNMaxima; i++) {
    if (fIPeak[i]>=0) {
      if (fHough[fMPeak[i]][fCPeak[i]]>nmaxbin) {
        nmaxbin = fHough[fMPeak[i]][fCPeak[i]];
        imax = ipeak;
      }
      ipeak++;
    }
  }
  return imax;
}

Double_t HoughTrackSR::GetSlope(Int_t i) const
{
  Int_t impeak = fMPeak[i];
  Int_t icpeak = fCPeak[i];
  Int_t nentry = fHough[impeak][icpeak];
  Int_t found=1;
  Int_t nskip=0;
  Int_t nfound = 0;
  Double_t mpeak = 0.;
  Double_t cpeak = 0.;
  for (Int_t ic1=icpeak-1; ic1>=0 && nskip<2; ic1--) {
    found = 0;
    Int_t im1 = (Int_t)((ic1-icpeak)*fSlopePeak[i]+.5)+impeak;
    for (Int_t im2=-3; im2<=3; im2++) {
      if (im1+im2>=0 && im1+im2<fMaxBin && fHough[im1+im2][ic1]==nentry) {
        found = 1;
        nfound++;
        mpeak += (Double_t)(im1+im2);
        cpeak += (Double_t)(ic1);
      }
    }
    if (!found) {
      nskip++;
    }
    else {
      nskip = 0;
    }
  }
  found=1;
  nskip=0;
  for (Int_t ic1=icpeak; ic1<fMaxBin && nskip<2; ic1++) {
    found = 0;
    Int_t im1 = (Int_t)((ic1-icpeak)*fSlopePeak[i]+.5)+impeak;
    for (Int_t im2=-3; im2<=3; im2++) {
      if (im1+im2>=0 && im1+im2<fMaxBin && fHough[im1+im2][ic1]==nentry) {
        found = 1;
        nfound++;
        mpeak += (Double_t)(im1+im2);
        cpeak += (Double_t)(ic1);
      }
    }
    if (!found) {
      nskip++;
    }
    else {
      nskip = 0;
    }
  }
  if (!nfound) return -1.;
  mpeak /= (Double_t)nfound;
  cpeak /= (Double_t)nfound;
  return mpeak;
}

Double_t HoughTrackSR::GetInter(Int_t i) const
{
  Int_t impeak = fMPeak[i];
  Int_t icpeak = fCPeak[i];
  Int_t nentry = fHough[impeak][icpeak];
  Int_t found=1;
  Int_t nskip=0;
  Int_t nfound = 0;
  Double_t mpeak = 0.;
  Double_t cpeak = 0.;
  for (Int_t ic1=icpeak-1; ic1>=0 && nskip<2; ic1--) {
    found = 0;
    Int_t im1 = (Int_t)((ic1-icpeak)*fSlopePeak[i]+.5)+impeak;
    for (Int_t im2=-3; im2<=3; im2++) {
      if (im1+im2>=0 && im1+im2<fMaxBin && fHough[im1+im2][ic1]==nentry) {
        found = 1;
        nfound++;
        mpeak += (Double_t)(im1+im2);
        cpeak += (Double_t)(ic1);
      }
    }
    if (!found) {
      nskip++;
    }
    else {
      nskip = 0;
    }
  }
  found=1;
  nskip=0;
  for (Int_t ic1=icpeak; ic1<fMaxBin && nskip<2; ic1++) {
    found = 0;
    Int_t im1 = (Int_t)((ic1-icpeak)*fSlopePeak[i]+.5)+impeak;
    for (Int_t im2=-3; im2<=3; im2++) {
      if (im1+im2>=0 && im1+im2<fMaxBin && fHough[im1+im2][ic1]==nentry) {
        found = 1;
        nfound++;
        mpeak += (Double_t)(im1+im2);
        cpeak += (Double_t)(ic1);
      }
    }
    if (!found) {
      nskip++;
    }
    else {
      nskip = 0;
    }
  }
  if (!nfound) return -1.;
  mpeak /= (Double_t)nfound;
  cpeak /= (Double_t)nfound;
  return cpeak;
}