/******************************************************************************
*  $Id: CbmRichRingSelect.cxx,v 1.1 2006/09/13 14:53:31 hoehne Exp $
*
*  Class  : CbmRichRingSelect
*  Description : Abstract base class for concrete RICH ring selection algorithm:
*                to be run after ring-track assign for fake-ring rejection
*
*  Author : Simeon Lebedev
*  E-mail : salebedev@jinr.ru
*
*******************************************************************************
*  $Log: CbmRichRingSelect.cxx,v $
*  Revision 1.1  2006/09/13 14:53:31  hoehne
*  initial version
*
*  
*
*******************************************************************************/



#include "CbmRichRingSelect.h"
#include "CbmTrackParam.h"

//------------------------------------------------------------------
Double_t CbmRichRingSelect::GetChi2(CbmRichRing* ring){
    Double_t chi2 = 0;
    Int_t nHits = ring->GetNoOfHits();
    if (nHits < 3) {
    
    }
    CbmRichHit* hitRing;
    for(Int_t iH = 0; iH < nHits; iH++){
      hitRing = (CbmRichHit*)ring->GetHit(iH);
      chi2 += pow( sqrt((ring->GetCenterX() - hitRing->X())*
                         (ring->GetCenterX() - hitRing->X()) +
		         (ring->GetCenterY() - hitRing->Y())*
			 (ring->GetCenterY() - hitRing->Y()) )
			 - ring->GetRadius(), 2 );
    }
    if (nHits < 4) {
      //cout << "-W- CbmRichRingSelect: Rings has less than 3 hits " << nHits << endl;
      return 999. ;
      }
   
    return chi2 = chi2 / (Double_t)(nHits - 3);
}
//------------------------------------------------------------------

//------------------------------------------------------------------
Double_t CbmRichRingSelect::GetRadPos(CbmRichRing* ring){
   Double_t radPos;
   if (ring->GetCenterY() > 0){
		radPos = sqrt((ring->GetCenterX() - 0)*
                        (ring->GetCenterX() - 0) +
		         		   (ring->GetCenterY() - 110)*
			 			      (ring->GetCenterY() - 110));
   } else {
      radPos = sqrt((ring->GetCenterX() - 0)*
                        (ring->GetCenterX() - 0) +
		         		   (ring->GetCenterY() + 110)*
			 			      (ring->GetCenterY() + 110));
   }

   return radPos;
}
//------------------------------------------------------------------

//------------------------------------------------------------------
Double_t CbmRichRingSelect::GetTBSum(CbmRichRing* ring){
   Int_t histTB[20];
   CbmRichHit* hitRing;
   Int_t nHits = ring->GetNoOfHits();
   for (Int_t iH = 0; iH < 20; iH++) histTB[iH] = 0;
      for(Int_t iH = 0; iH < nHits; iH++){
         hitRing = (CbmRichHit*)ring->GetHit(iH);
         Double_t r = pow((hitRing->X() - ring->GetCenterX())*
                        (hitRing->X() - ring->GetCenterX()) +
                        (hitRing->Y() - ring->GetCenterY())*
                        (hitRing->Y() - ring->GetCenterY()), 0.5);

         if (r > 3.5 && r < 6.5) {
            Double_t ir ;
            modf(r / (6.5 / 20.0), &ir);
            if (ir < 20 && ir > 0 ) histTB[ (Int_t) ir ]++;
         }
      }
   Int_t ii;
   Int_t max = 0;
   for (Int_t i = 0; i < 20; i++)
      if (histTB[i] >= max) {
         max = histTB[i];
         ii = i;
      }
   UInt_t maxsum = histTB[ii];
   if (ii+1 < 20) maxsum += histTB[ii+1];
   if (ii-1 > 0) maxsum += histTB[ii-1];
 
   return maxsum;
}
//------------------------------------------------------------------

//------------------------------------------------------------------
Double_t CbmRichRingSelect::GetAngle(CbmRichRing* ring){
  
   vector<Double_t> alpha;
   vector<Double_t> phi;
	phi.clear();
	alpha.clear();
   Double_t Pi = 3.1415;
   CbmRichHit * hitRing;
   Int_t nHits = ring->GetNoOfHits();
   Double_t xRing = ring->GetCenterX();
   Double_t yRing = ring->GetCenterY();
  
   if (nHits < 2) return 999.;
   for(Int_t iHit = 0; iHit < nHits; iHit++){
      hitRing = (CbmRichHit*)ring->GetHit(iHit);
      if (!hitRing) {
         cout << "-E- Double_t CbmRichRingSelect::GetAngle(CbmRichRing* ring)"<< iHit <<endl;
         continue;
      }
      
      if ((hitRing->Y()-yRing) == 0) {
          cout << " -W- CbmRichRingSelect 0 in angle determination, y " << endl;
	  return 999.;
	  }
      if ((hitRing->X()-xRing) == 0) {
          cout << " -W- CbmRichRingSelect 0 in angle determination, x " << endl;
	  return 999.;
	  }
	  
      if( hitRing->X() > xRing && hitRing->Y() > yRing ){
	      alpha.push_back(atan(fabs((hitRing->X()-xRing)/
                                    (hitRing->Y()-yRing))));
      }
      if( hitRing->X() < xRing && hitRing->Y() > yRing ){
        alpha.push_back(Pi - atan(fabs((hitRing->X()-xRing)/
                                       (hitRing->Y()-yRing))));
      }
      if( hitRing->X() < xRing && hitRing->Y() < yRing ){
        alpha.push_back(Pi + atan(fabs((hitRing->X()-xRing)/
                                       (hitRing->Y()-yRing))));
      }
      if( hitRing->X() > xRing && hitRing->Y() < yRing ){
        alpha.push_back(2*Pi - atan(fabs((hitRing->X()-xRing)/
                                         (hitRing->Y()-yRing))));
      }
    }
    
    sort(alpha.begin(),alpha.end());

    for(Int_t i = 0; i < nHits-1; i++) phi.push_back(alpha[i+1] - alpha[i]);
    phi.push_back(2*Pi - alpha[nHits-1] + alpha[0]);
    sort(phi.begin(),phi.end());
    
    Double_t bigestAngle = phi[nHits-1]+phi[nHits-2]+phi[nHits-3];

    return bigestAngle;
}
//------------------------------------------------------------------

//------------------------------------------------------------------
Double_t CbmRichRingSelect::GetTrackDist(CbmRichRing* ring, TClonesArray *projAr){
    Double_t xRing = ring->GetCenterX();
    Double_t yRing = ring->GetCenterY();
    if (ring->GetTrackID() == -1) return 999;
    CbmTrackParam* proj= (CbmTrackParam*)projAr->At(ring->GetTrackID());
    Double_t xProj = proj->GetX();
    Double_t yProj = proj->GetY();
    if (xProj == 0 && yProj == 0 ) return 999;    // no projection to photodetector plane
    Double_t distRingTrack = TMath::Sqrt( (xRing-xProj)*(xRing-xProj) +
                                (yRing-yProj)*(yRing-yProj) );
				
   
    return distRingTrack;
}
//------------------------------------------------------------------

//------------------------------------------------------------------
Double_t CbmRichRingSelect::GetNofHits(CbmRichRing* ring){
   Int_t nHits = ring->GetNoOfHits();
 
   return nHits;
}
//------------------------------------------------------------------


ClassImp(CbmRichRingSelect)