/* $Id: CbmRichRingFinderTrackCircleFit.cxx,v 1.1 2006/01/19 10:56:12 hoehne Exp $*/

/* History of CVS commits:
 *
 * $Log: CbmRichRingFinderTrackCircleFit.cxx,v $
 * Revision 1.1  2006/01/19 10:56:12  hoehne
 * Rename CbmRichRingFinder.cxx to CbmRichRingFinderTrackCircleFit.cxx
 *
 *
 *
 */

// -------------------------------------------------------------------------
// -----         CbmRichRingFinderTrackCircleFit source file           -----
// -----      Created 01/07/04 A.Soloviev  <solovjev@cv.jinr.ru>       -----
// -------------------------------------------------------------------------


#include "CbmRichRingFinder.h"
#include "CbmRichRingFinderImp.h"
#include "CbmRichHit.h"
#include "CbmRichRing.h"
#include "CbmTrackParam.h"

#include <iostream>
#include <algorithm>
using namespace std;

// -----   Default constructor   -------------------------------------------
CbmRichRingFinder::CbmRichRingFinder()
  : CbmTask(),
    fManager ( CbmRootManager::Instance() ),
    fThreshold ( 3 ),
    fQuiet ( kFALSE ),
    fNEvent ( 0 )
{}
// -------------------------------------------------------------------------


// -----   Standard constructor   ------------------------------------------
CbmRichRingFinder::CbmRichRingFinder ( const char * name, const char * title )
  : CbmTask ( name ),
    fManager ( CbmRootManager::Instance() ),
    fThreshold ( 3 ),
    fQuiet ( kFALSE ),
    fNEvent ( 0 )
{}
// -------------------------------------------------------------------------


// -----   Destructor   ----------------------------------------------------
CbmRichRingFinder::~CbmRichRingFinder()
{
//   if ( fManager ) fManager->Write();
}
// -------------------------------------------------------------------------


// -----   Initialization   ------------------------------------------------
InitStatus CbmRichRingFinder::Init()
{
  if ( !fManager ) return kERROR;

  fHitCollection = 
    (TClonesArray*)fManager->GetObject ( "RichHit" );

  fProjectionCollection =
    (TClonesArray*)fManager->GetObject ( "RichRingProjection" );

  fRingCollection = new TClonesArray ( "CbmRichRing" );
  fManager->Register ( "RichRing", "RICH", fRingCollection ,kTRUE);

  return kSUCCESS;
}
// -------------------------------------------------------------------------


// -----   Execution of Task   ---------------------------------------------
void CbmRichRingFinder::Exec ( Option_t * option )
{
  ++fNEvent;
  cout << "Event no. " << fNEvent << endl;
  
  Double_t Radius = 6.;
  Double_t RadiusMax = 8.;

  static CbmRichRingFinderImp sigma;
  sigma.clear();

  const Int_t nhits = fHitCollection->GetEntriesFast();
  const Int_t nrings = fProjectionCollection->GetEntriesFast();

  bool robust = true; // robust fitting?

  for ( register Int_t i = 0; i < nhits; ++i ) {
    CbmRichHit * hit = (CbmRichHit*)fHitCollection->At ( i );
    if ( hit ) sigma.add_data_point ( hit->X(), hit->Y() );
  }

  if ( nrings ) { // have ring Projections...
    for ( register Int_t i = 0; i < nrings; ++i ) {
      CbmTrackParam * Projection = (CbmTrackParam*)fProjectionCollection->At ( i );
      if ( Projection ) {
	sigma.set_Projection ( Projection->GetX(),
			     Projection->GetY(),
			     Radius );

	sigma.set_threshold ( fThreshold );
	bool fake = 
	  !sigma.minimize ( robust ) ||
	  sigma.ring_radius() > 1.2 * RadiusMax;

	if ( !fake ) { // reject duplicates...
	  for ( Int_t j = 0; j < fRingCollection->GetEntriesFast(); ++j ) {
	    CbmRichRing * r = (CbmRichRing*)fRingCollection->At ( j );
	    if ( r ) {
	      double dx = r->GetCenterX() - sigma.ring_center_x();
	      double dy = r->GetCenterY() - sigma.ring_center_y();
	      double dist = sqrt ( (dx*=dx) + (dy*=dy) );
	      bool criteria4reject = 
		dist < min ( r->GetRadius(), sigma.ring_radius() );
	      if ( criteria4reject ) {
		if ( r->GetUncertainty() > sigma.value() ) {
		  r->SetCenterX ( sigma.ring_center_x() );
		  r->SetCenterY ( sigma.ring_center_y() );
		  r->SetRadius ( sigma.ring_radius() );
		  r->SetTrackID ( Projection->GetTrackID() );
		  r->SetUncertainty ( sigma.value() );
		}
		fake = true; // this ring is already taken
		break; // nothing else to do!
	      }
	    }
	  }
	}

	if ( !fake ) {
	  CbmRichRing ring ( sigma.ring_center_x(), 
			     sigma.ring_center_y(), 
			     sigma.ring_radius() );
	  ring.SetTrackID ( Projection->GetTrackID() );
	  ring.SetUncertainty ( sigma.value() );
	  new ( (*fRingCollection)[fRingCollection->GetEntriesFast()] ) 
	    CbmRichRing ( ring );

	  if ( !fQuiet ) 
	    cout << "CbmRichRingFinder::Exec --[" 
		 << fRingCollection->GetEntriesFast() 
		 << "]--> " << sigma << endl;
	}
      }
    }
  } else { // have NO ring Projections...
    sigma.find_Projection();

    sigma.minimize ( robust );

    new ( (*fRingCollection)[fRingCollection->GetEntriesFast()] ) 
      CbmRichRing ( sigma.ring_center_x(), 
		    sigma.ring_center_y(), 
		    sigma.ring_radius() );

    if ( !fQuiet ) 
      cout << "CbmRichRingFinder::Exec ----> " << sigma << endl;
  }

  cout << "Number of input hits: " 
       << fHitCollection->GetEntriesFast() << endl;
  cout << "Number of input Projections: "
       << fProjectionCollection->GetEntriesFast() << endl;
  cout << "Number of output rings: " 
       << fRingCollection->GetEntriesFast() << endl;

}
// -------------------------------------------------------------------------


// -----   Finish Task   ---------------------------------------------------
void CbmRichRingFinder::Finish()
{
  if ( fRingCollection ) fRingCollection->Clear();
  if ( fProjectionCollection ) fProjectionCollection->Clear();
  if ( fHitCollection ) fHitCollection->Clear();
}
// -------------------------------------------------------------------------



void CbmRichRingFinder::SetThreshold ( const Double_t & threshold )
{
  fThreshold = threshold;
}

void CbmRichRingFinder::SetQuiet ( Bool_t quiet )
{
  fQuiet = quiet;
}

ClassImp(CbmRichRingFinder)