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

#include <iostream>
#include <string>
#include <cmath>

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

ClassImp(HoughViewSR)

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

//-------------------------------------------------------------------
HoughViewSR::HoughViewSR(TObjArray clusterlist)
{
  // set default values for hough fitter parameters.
  // these get reset by AlgTrackSRList
  fPeakMin = 4;
  fPeakMinFrac = .25;
  fPeakMinFracZoom = .75;
  fMinInterBinSize = 0.02;
  fClusterList = new TObjArray(clusterlist);
  Init();
}

//-------------------------------------------------------------------
HoughViewSR::HoughViewSR()
{
 // set default values for hough fitter parameters.
  // these get reset by AlgTrackSRList 
  fPeakMin = 4;
  fPeakMinFrac = .5;
  fPeakMinFracZoom = .75;
  fMinInterBinSize = 0.02;
}

//-------------------------------------------------------------------
HoughViewSR::~HoughViewSR()
{
  if (fClusterList) {
    delete fClusterList;
  }
  for (Int_t i=0; i<=fHoughTrackList.GetLast(); i++) {
    HoughTrackSR *htrack = dynamic_cast<HoughTrackSR*>(fHoughTrackList.At(i));
    if (htrack) delete htrack;
  }
}
//-------------------------------------------------------------------
void HoughViewSR::SetClusterList(TObjArray clusterlist)
{
  fClusterList = new TObjArray(clusterlist);
  Init();

}

//-------------------------------------------------------------------
void HoughViewSR::SetClusterList(TObjArray *  clusterlist,  Double_t minPulseHeight)
{
  fClusterList = new TObjArray(1,0);

 for (Int_t i=0; i<=clusterlist->GetLast(); i++) {
    TrackClusterSR *cluster = dynamic_cast<TrackClusterSR*>(clusterlist->At(i)); 
    if(cluster->GetCharge()>minPulseHeight){
      
      //loop over the cluster list and look to see if there is a cluster from
      //this clusters plane already in the list.  if so, check to see if the 
      //current cluster has a larger signal and if it does replace the one in the
      //list with the current cluster
      fClusterList->AddLast(cluster);
      
      MSG("TrackSR", Msg::kDebug) << "plane " << cluster->GetPlane() << " cluster " 
				  <<  cluster->GetMaxStrip() << " - " << cluster->GetMinStrip() << endl;
    }
  }
  
  MSG("TrackSR", Msg::kDebug) << fClusterList->GetLast()+1 
			      << " clusters in the HoughView" << endl;
  Init();


}



//-------------------------------------------------------------------
void HoughViewSR::Init()
{
  fXZmean[0] = 0.;
  fXZmean[1] = 0.;
  fXZrms[0] = 0.;
  fXZrms[1] = 0.;
 // reduce requirement for min peak height for < 5 track clusters 
  if(fClusterList->GetLast()<4)fPeakMin=3;
  int ncluster=0;
  for (int i=0; i<=fClusterList->GetLast(); ++i) {
    TrackClusterSR *cluster = dynamic_cast<TrackClusterSR*>(fClusterList->At(i));
    
    Double_t tpos = cluster->GetTPos();
    Double_t zpos = cluster->GetZPos();
    ++ncluster;
    fXZmean[0] += tpos;
    fXZrms[0] += tpos*tpos;
    fXZmean[1] += zpos;
    fXZrms[1] += zpos*zpos;
  }
  if (ncluster) {
    fXZmean[0] /= (Double_t)ncluster;
    fXZmean[1] /= (Double_t)ncluster;
    fXZrms[0] /= (Double_t)ncluster;
    fXZrms[1] /= (Double_t)ncluster;
    fXZrms[0] -= fXZmean[0]*fXZmean[0];
    if (fXZrms[0]>0.) {
      fXZrms[0] = sqrt(fXZrms[0]);
    }
    else {
      fXZrms[0] = 0.;
    }
    fXZrms[1] -= fXZmean[1]*fXZmean[1];
    if (fXZrms[1]>0.) {
      fXZrms[1] = sqrt(fXZrms[1]);
    }
    else {
      fXZrms[1] = 0.;
    }
  }
  HoughTrackSR *houghtrack = new HoughTrackSR(*fClusterList);
  houghtrack->SetX0(fXZmean[0]);
  houghtrack->SetZ0(fXZmean[1]);
  houghtrack->SetXRMS(fXZrms[0]);
  houghtrack->SetZRMS(fXZrms[1]);
  houghtrack->SetPeakMin(fPeakMin);
  houghtrack->SetPeakMinFrac(fPeakMinFrac);
  houghtrack->SetPeakMinFracZoom(fPeakMinFracZoom);
  houghtrack->SetMinInterBinSize(fPeakMinFracZoom);
  houghtrack->SetMaxBin(40);
  //  if(fClusterList->GetLast()<10) houghtrack->SetMaxBin(20);
  houghtrack->Iterate();
  fHoughTrackList.Add(houghtrack);
  
}

//-------------------------------------------------------------------
void HoughViewSR::Print(Option_t* /* option */) const
{
  for (Int_t i=0; i<=fHoughTrackList.GetLast(); i++) {
    HoughTrackSR *houghtrack = dynamic_cast<HoughTrackSR*>(fHoughTrackList.At(i));
    cout << "**************************** TRACK " << i << "/" << fHoughTrackList.GetLast() << " ****************************\n\n";
    houghtrack->Print();
  }
}

//-------------------------------------------------------------------
Int_t HoughViewSR::Iterate(Int_t imode)
{

  //imode<=0: max # iterations = 20
  //imode>0: max # iterations = imode
  if (imode<=0) {
    Int_t niterate=0;
    Int_t icount=0;
    do {
      MSG("TrackSR",Msg::kDebug) << "ITERATING " << icount << "\n";
      niterate = IterateSingle();
      ++icount;
    } while (niterate && icount<20);
    return icount;
  }
  else {
    Int_t niterate=0;
    for (Int_t i=0; i<imode; ++i) {
      niterate = IterateSingle();
    }
    return niterate;
  }
}

//-------------------------------------------------------------------
Int_t HoughViewSR::IterateSingle()
{
  MSG("TrackSR",Msg::kDebug) << "begin IterateSingle\n";
  TObjArray tracklist;
  Int_t niterate=0;
 
  for (Int_t i=0; i<=fHoughTrackList.GetLast(); ++i) {
    MSG("TrackSR",Msg::kDebug) << "TRACK " << i << "/" << fHoughTrackList.GetLast() << "\n";
    HoughTrackSR *houghtrack = dynamic_cast<HoughTrackSR*>(fHoughTrackList.At(i));
    MSG("TrackSR",Msg::kDebug) << "OLD\n"; 
    //            houghtrack->Print();
    HoughTrackSR *oldHT = 0;
    if (houghtrack->GetNPeak()>0 && houghtrack->GetInterceptMax()-houghtrack->GetInterceptMin()>fMinInterBinSize*(Double_t)houghtrack->GetNBin()) {
      oldHT = new HoughTrackSR(*houghtrack);
      TObjArray newlist(houghtrack->Iterate());
      for (Int_t i=0; i<=newlist.GetLast(); ++i) {
        HoughTrackSR *htrack = dynamic_cast<HoughTrackSR*>(newlist.At(i));
        MSG("TrackSR",Msg::kDebug) << "newly created " << i << "\n";
	//  	htrack->Print();
        if (htrack->GetNPeak()>0) {
	  fHoughTrackList.Add(htrack);
          MSG("TrackSR",Msg::kDebug) << "adding to list\n";
        }
        else {
          delete htrack;
        }
      }
      MSG("TrackSR",Msg::kDebug) << "NEW\n";
      //   houghtrack->Print();
      if (houghtrack->GetNPeak()>0) {
        ++niterate;
        delete oldHT;
      }
      else if (newlist.GetLast()<0) {
        delete houghtrack;
        houghtrack = oldHT;
        MSG("TrackSR",Msg::kDebug) << "NO PEAKS, REVERTING BACK\n";
      }
      else {
        delete oldHT;
      }
    }
    if (houghtrack->GetNPeak()>0) {
      tracklist.Add(houghtrack);
    }
    else {
      delete houghtrack;
    }
    MSG("TrackSR",Msg::kDebug) << "tracklist size " << tracklist.GetLast() << "\n";
  }
  MSG("TrackSR",Msg::kDebug) << "clearing houghtracklist\n";
  fHoughTrackList.Clear();
  fHoughTrackList.Compress();
  for (Int_t i=0; i<=tracklist.GetLast(); ++i) {
  MSG("TrackSR",Msg::kDebug) << "adding to houghtracklist " << i << "/" << tracklist.GetLast() << "\n";
    HoughTrackSR *htrack = dynamic_cast<HoughTrackSR*>(tracklist.At(i));
    fHoughTrackList.Add(htrack);
  }
  MSG("TrackSR",Msg::kDebug) << "end iteratesingle\n";
  return niterate;
}

//-------------------------------------------------------------------
void HoughViewSR::SetPeakMin(Int_t nbin)
{
  fPeakMin = nbin;
}

//-------------------------------------------------------------------
void HoughViewSR::SetPeakMinFrac(Double_t frac)
{
  fPeakMinFrac = frac;
}

//-------------------------------------------------------------------
void HoughViewSR::SetPeakMinFracZoom(Double_t frac)
{
  fPeakMinFracZoom = frac;
}

//-------------------------------------------------------------------
void HoughViewSR::SetMinInterBinSize(Double_t frac)
{
  fMinInterBinSize = frac;
}

//-------------------------------------------------------------------
Double_t HoughViewSR::GetX0()
{
  return fXZmean[0];
}

//-------------------------------------------------------------------
Double_t HoughViewSR::GetZ0()
{
  return fXZmean[1];
}

//-------------------------------------------------------------------
Int_t HoughViewSR::GetNTrack() const
{
  return fHoughTrackList.GetLast()+1;
}

//-------------------------------------------------------------------
Int_t HoughViewSR::GetNPeak() const
{
  Int_t npeak = 0;
  for (Int_t i=0; i<=fHoughTrackList.GetLast(); ++i) {
    const HoughTrackSR *houghtrack = dynamic_cast<const HoughTrackSR*>(fHoughTrackList.At(i));
    npeak += houghtrack->GetNPeak();
  }
  return npeak;
}

//-------------------------------------------------------------------
const HoughTrackSR *HoughViewSR::GetHoughTrack(Int_t itrack)
{
  if (itrack<=fHoughTrackList.GetLast()) {
    const HoughTrackSR *houghtrack = dynamic_cast<const HoughTrackSR*>(fHoughTrackList.At(itrack));
    return houghtrack;
  }
  else {
    return 0;
  }
}

//-------------------------------------------------------------------
Int_t HoughViewSR::GetLargestTrackIndex() const
{
  Int_t itrack = -9999;
  Int_t nmaxbin = 0;
  for (Int_t i=0; i<=fHoughTrackList.GetLast(); ++i) {
    const HoughTrackSR *houghtrack = dynamic_cast<const HoughTrackSR*>(fHoughTrackList.At(i));
    if (houghtrack->GetPeakBin(houghtrack->GetLargestPeakIndex())>nmaxbin) {
      nmaxbin = houghtrack->GetPeakBin(houghtrack->GetLargestPeakIndex());
      itrack = i;
    }
  }
  return itrack;
}