beamtime/cern-oct11/go4/RICH/CbmRichRingFinderHoughImpl.cxx (r4864/r3566)

Go to the documentation of this file.

#include "CbmRichRingFinderHoughImpl.h"
#include "CbmRichRingLight.h"
#include "CbmRichRingFitterCOP.h"

#include <iostream>
#include <cmath>
#include <set>
#include <algorithm>
#include <iostream>

using std::cout;
using std::endl;
using std::vector;
using std::set;
using std::sort;

CbmRichRingFinderHoughImpl::CbmRichRingFinderHoughImpl():
   fNofParts(0),

   fMaxDistance(0.f),
   fMinDistance(0.f),
   fMinDistanceSq(0.f),
   fMaxDistanceSq(0.f),

   fMinRadius(0.f),
   fMaxRadius(0.f),

   fDx(0.f),
   fDy(0.f),
   fDr(0.f),
   fNofBinsX(0),
   fNofBinsY(0),
   fNofBinsXY(0),

   fHTCut(0),

   fNofBinsR(0),
   fHTCutR(0),

   fMinNofHitsInArea(0),

   fRmsCoeffCOP(0.f),
   fMaxCutCOP(0.f),

   fCurMinX(0.f),
   fCurMinY(0.f),

   fData(),
   fHist(),
   fHistR(),
   fHitInd(),
   fFoundRings(),
   fFitCOP(NULL)
{
   SetParameters();

   fHist.resize(fNofBinsXY);
   fHistR.resize(fNofBinsR);
   fHitInd.resize(fNofParts);

   fFitCOP = new CbmRichRingFitterCOP();
}

CbmRichRingFinderHoughImpl::~CbmRichRingFinderHoughImpl()
{
        if (NULL != fFitCOP) delete fFitCOP;
}

void CbmRichRingFinderHoughImpl::DoFind(
                const vector<CbmRichHitLight>& hits)
{
        for (int i = 0 ; i < fFoundRings.size(); i++){
                if (fFoundRings[i] != NULL) delete fFoundRings[i];
        }
        fFoundRings.clear();
        fData.clear();

        if (hits.size() < 1) return;

        if (fData.size() > MAX_NOF_HITS) {
           cout << endl << endl << "-E- CbmRichRingFinderHoughImpl::DoFind"
                 << ". Number of hits is more than " << MAX_NOF_HITS << endl << endl;
           return;
        }

        for (int iH = 0; iH < hits.size(); iH++){
                CbmRichHoughHit tempPoint;
                tempPoint.fHit.fX = hits[iH].fX;
                tempPoint.fHit.fY = hits[iH].fY;
                tempPoint.fX2plusY2 = hits[iH].fX * hits[iH].fX + hits[iH].fY * hits[iH].fY;
                tempPoint.fId = iH;
                tempPoint.fIsUsed = false;
                fData.push_back(tempPoint);
        }


        std::sort(fData.begin(), fData.end(), CbmRichHoughHitCmpUp());
        HoughTransformReconstruction();
        fData.clear();
}

void CbmRichRingFinderHoughImpl::SetParameters()
{
        fMaxDistance = 115.;
        fMinDistance = 10.;
        fMinDistanceSq = fMinDistance*fMinDistance;
        fMaxDistanceSq = fMaxDistance*fMaxDistance;

        fMinRadius = 20.;
        fMaxRadius = 65.;

        fHTCut = 7;
        fHTCutR = 5;
        fMinNofHitsInArea = 4;

        fNofBinsX = 25;
        fNofBinsY = 25;
        fNofBinsR = 32;

        fRmsCoeffCOP = 3.;
        fMaxCutCOP = 10.;

        fNofParts = 1;
   fDx = 2.f * fMaxDistance / (float)fNofBinsX;
   fDy = 2.f * fMaxDistance / (float)fNofBinsY;
   fDr = fMaxRadius / (float)fNofBinsR;
   fNofBinsXY = fNofBinsX * fNofBinsY;
}

void CbmRichRingFinderHoughImpl::HoughTransformReconstruction()
{
   int indmin, indmax;
   unsigned int size = fData.size();

   fCurMinX = fData[0].fHit.fX - fMaxDistance;
   fCurMinY = fData[0].fHit.fY - fMaxDistance;

   DefineLocalAreaAndHits(fData[0].fHit.fX, fData[0].fHit.fY , &indmin, &indmax);
   HoughTransform(indmin, indmax);
   FindPeak(indmin, indmax);

   DefineLocalAreaAndHits(fData[0].fHit.fX, fData[0].fHit.fY , &indmin, &indmax);
   HoughTransform(indmin, indmax);
   FindPeak(indmin, indmax);
}

void CbmRichRingFinderHoughImpl::DefineLocalAreaAndHits(
      float x0,
      float y0,
      int *indmin,
          int *indmax)
{
   CbmRichHoughHit mpnt;
   vector<CbmRichHoughHit>::iterator itmin, itmax;

        //find all hits which are in the corridor
        mpnt.fHit.fX = x0 - fMaxDistance;
        itmin = std::lower_bound(fData.begin(), fData.end(), mpnt, CbmRichHoughHitCmpUp());

        mpnt.fHit.fX = x0 + fMaxDistance;
        itmax = std::lower_bound(fData.begin(), fData.end(), mpnt, CbmRichHoughHitCmpUp()) - 1;

        *indmin = itmin - fData.begin();
        *indmax = itmax - fData.begin();

        int arSize = *indmax - *indmin + 1;
        if (arSize <= fMinNofHitsInArea) return;

        for (unsigned short i = 0; i < fNofParts; i++){
                fHitInd[i].clear();
                fHitInd[i].reserve( (*indmax-*indmin) / fNofParts);
        }

        unsigned short indmin1 = (unsigned short)(*indmin);
        unsigned short indmax1 = (unsigned short)(*indmax);

   for (unsigned short i = indmin1; i <= indmax1; i++) {
      if (fData[i].fIsUsed == true) continue;
      float ry = y0 - fData[i].fHit.fY;
      if (fabs(ry) > fMaxDistance) continue;
      float rx = x0 - fData[i].fHit.fX;
      float d = rx * rx + ry * ry;
      if (d > fMaxDistanceSq) continue;
           fHitInd[i % fNofParts].push_back(i);
   }

        for (unsigned short j = 0; j < fNofBinsXY; j++){
                fHist[j] = 0;
        }

        for (unsigned short k = 0; k < fNofBinsR; k++) {
                fHistR[k] = 0;
        }
}

void CbmRichRingFinderHoughImpl::HoughTransform(
      unsigned short int indmin,
      unsigned short int indmax)
{
    for (int iPart = 0; iPart < fNofParts; iPart++){
         HoughTransformGroup(indmin, indmax, iPart);
    }//iPart
}

void CbmRichRingFinderHoughImpl::HoughTransformGroup(
      unsigned short int indmin,
                unsigned short int indmax,
                int iPart)
{
   unsigned short nofHits = fHitInd[iPart].size();
   float xcs, ycs; // xcs = xc - fCurMinX
   float dx = 1.0f/fDx, dy = 1.0f/fDy, dr = 1.0f/fDr;

   vector<CbmRichHoughHit> data;
   data.resize(nofHits);
   for (int i = 0; i < nofHits; i++){
      data[i] = fData[ fHitInd[iPart][i] ];
   }

   typedef vector<CbmRichHoughHit>::iterator iH;

        for (iH iHit1 = data.begin(); iHit1 != data.end(); iHit1++) {
           float iH1X = iHit1->fHit.fX;
                float iH1Y = iHit1->fHit.fY;

      for (iH iHit2 = iHit1 + 1; iHit2 != data.end(); iHit2++) {
         float iH2X = iHit2->fHit.fX;
         float iH2Y = iHit2->fHit.fY;

         float rx0 = iH1X - iH2X;//rx12
         float ry0 = iH1Y- iH2Y;//ry12
         float r12 = rx0 * rx0 + ry0 * ry0;
         if (r12 < fMinDistanceSq || r12 > fMaxDistanceSq)      continue;

         float t10 = iHit1->fX2plusY2 - iHit2->fX2plusY2;
         for (iH iHit3 = iHit2 + 1; iHit3 != data.end(); iHit3++) {
            float iH3X = iHit3->fHit.fX;
            float iH3Y = iHit3->fHit.fY;

            float rx1 = iH1X - iH3X;//rx13
            float ry1 = iH1Y - iH3Y;//ry13
            float r13 = rx1 * rx1 + ry1 * ry1;
            if (r13 < fMinDistanceSq || r13 > fMaxDistanceSq)continue;

            float rx2 = iH2X - iH3X;//rx23
            float ry2 = iH2Y - iH3Y;//ry23
            float r23 = rx2 * rx2 + ry2 * ry2;
            if (r23     < fMinDistanceSq || r23 > fMaxDistanceSq)continue;

            float det = rx2*ry0 - rx0*ry2;
            if (det == 0.0f) continue;
            float t19 = 0.5f / det;
            float t5 = iHit2->fX2plusY2 - iHit3->fX2plusY2;

            float xc = (t5 * ry0 - t10 * ry2) * t19;
            xcs = xc - fCurMinX;
            int intX = int( xcs *dx);
            if (intX < 0 || intX >= fNofBinsX ) continue;

            float yc = (t10 * rx2 - t5 * rx0) * t19;
            ycs = yc - fCurMinY;
            int intY = int( ycs *dy);
            if (intY < 0 || intY >= fNofBinsY ) continue;

            //radius calculation
            float t6 = iH1X - xc;
            float t7 = iH1Y - yc;
            //if (t6 > fMaxRadius || t7 > fMaxRadius) continue;
            float r = sqrt(t6 * t6 + t7 * t7);
            //if (r < fMinRadius) continue;
            int intR = int(r *dr);
            if (intR < 0 || intR >= fNofBinsR) continue;

            fHist[intX * fNofBinsX + intY]++;
            fHistR[intR]++;
         }// iHit1
      }// iHit2
   }// iHit3
        //hitIndPart.clear();
}

void CbmRichRingFinderHoughImpl::FindPeak(
      int indmin,
      int indmax)
{
   // Find MAX bin R and compare it with CUT
   int maxBinR = -1, maxR = -1;
   unsigned int size = fHistR.size();
   for (unsigned int k = 0; k < size; k++){
      if (fHistR[k] > maxBinR){
         maxBinR = fHistR[k];
         maxR = k;
      }
   }
        if (maxBinR < fHTCutR) return;

   // Find MAX bin XY and compare it with CUT
        int maxBinXY = -1, maxXY = -1;
   size = fHist.size();
   for (unsigned int k = 0; k < size; k++){
      if (fHist[k] > maxBinXY){
         maxBinXY = fHist[k];
         maxXY = k;
      }
   }
   if (maxBinXY < fHTCut) return;

        CbmRichRingLight* ring1 = new CbmRichRingLight();

   // Find Preliminary Xc, Yc, R
   float xc, yc, r;
   float rx, ry, dr;
        xc = (maxXY / fNofBinsX + 0.5f)* fDx + fCurMinX;
        yc = (maxXY % fNofBinsX + 0.5f)* fDy + fCurMinY;
        r = (maxR + 0.5f) * fDr;
        for (int j = indmin; j < indmax + 1; j++) {
                rx = fData[j].fHit.fX - xc;
                ry = fData[j].fHit.fY - yc;

                dr = fabs(sqrt(rx * rx + ry * ry) - r);
                if (dr > 6.) continue;
                ring1->AddHit(fData[j].fHit, fData[j].fId);
        }

        if (ring1->GetNofHits() < 5) {
           delete ring1;
           return;
        }

        fFitCOP->DoFit(ring1);
        float drCOPCut = fRmsCoeffCOP * sqrt(ring1->GetChi2() / ring1->GetNofHits());
        if (drCOPCut > fMaxCutCOP)      drCOPCut = fMaxCutCOP;

        xc = ring1->GetCenterX();
        yc = ring1->GetCenterY();
        r = ring1->GetRadius();

        delete ring1;

        CbmRichRingLight* ring2 = new CbmRichRingLight();
        for (int j = indmin; j < indmax + 1; j++) {
                rx = fData[j].fHit.fX - xc;
                ry = fData[j].fHit.fY - yc;

                dr = fabs(sqrt(rx * rx + ry * ry) - r);
                if (dr > drCOPCut) continue;
                fData[j].fIsUsed = true;
                ring2->AddHit(fData[j].fHit, fData[j].fId);
        }
        if (ring2->GetNofHits() < 5) {
           delete ring2;
           return;
        }

   fFitCOP->DoFit(ring2);
   fFoundRings.push_back(ring2);
        ring2=NULL;
        delete ring2;
}