CbmL1RichENNRingFinderParallel.cxx

Go to the documentation of this file.

/*
 *====================================================================
 *
 *  CBM Level 1 Reconstruction 
 *  
 *  Authors: I.Kisel,  S.Gorbunov
 *
 *  e-mail : ikisel@kip.uni-heidelberg.de 
 *
 *====================================================================
 *
 *  Standalone RICH ring finder based on the Elastic Neural Net
 *
 *====================================================================
 */
 
#define PRINT_TIMING
 
// CBM includes
#include "CbmL1RichENNRingFinderParallel.h"
 
#include "CbmRichHit.h"
#include "CbmRichRing.h"
 
#include "TClonesArray.h"
#include "TStopwatch.h"
 
#include "assert.h"
#include <algorithm>
#include <cmath>
#include <iostream>
#include <vector>
 
using std::cout;
using std::endl;
using std::fabs;
using std::ios;
using std::sqrt;
using std::vector;
 
CbmL1RichENNRingFinderParallel::CbmL1RichENNRingFinderParallel(Int_t verbose)
  : fRecoTime(0), fNEvents(0) {
  fVerbose = verbose;
#ifdef PRINT_TIMING
  TString name_tmp[NTimers] = {
 
    "All",
 
    "Ring finding",
    "Ring finding: Prepare data",
    "Ring finding: Init of param",
    "Ring finding: Hit selection",
    "Ring finding: Final fit",
    "Ring finding: Store ring",
 
    "Selection",
 
    "Sort",
    "Repack",
    "Copy data"};
  for (int i = 0; i < NTimers; i++) {
    fTimersNames[i] = name_tmp[i];
  }
#endif  // PRINT_TIMING
}
 
CbmL1RichENNRingFinderParallel::~CbmL1RichENNRingFinderParallel() {}
 
 
void CbmL1RichENNRingFinderParallel::Init() {}
 
Int_t CbmL1RichENNRingFinderParallel::DoFind(TClonesArray* HitArray,
                                             TClonesArray* /*ProjArray*/,
                                             TClonesArray* RingArray) {
  if (!RingArray || !HitArray) return 0;
 
  RingArray->Clear();
#ifdef PRINT_TIMING
  for (int i = 0; i < NTimers; i++) {
    fTimers[i].Reset();
  }
#endif  // PRINT_TIMING
 
#ifdef PRINT_TIMING
  GetTimer("Copy data").Start(0);
#endif  // PRINT_TIMING
  vector<ENNRingHit> Up;
  vector<ENNRingHit> Down;
 
  const Int_t nhits = HitArray->GetEntriesFast();
 
  for (register Int_t i = 0; i < nhits; ++i) {
    CbmRichHit* hit = L1_DYNAMIC_CAST<CbmRichHit*>(HitArray->At(i));
    if (!hit) continue;
    ENNRingHit tmp;
    tmp.x        = hit->GetX();
    tmp.y        = hit->GetY();
    tmp.outIndex = i;
    tmp.quality  = 0;
    if (tmp.y > 0.) {
      Up.push_back(tmp);
    } else {
      Down.push_back(tmp);
    }
  }
 
 
#ifdef PRINT_TIMING
  GetTimer("Copy data").Stop();
#endif  // PRINT_TIMING
 
#ifdef PRINT_TIMING
  GetTimer("Sort").Start(0);
#endif  // PRINT_TIMING
  sort(Up.begin(), Up.end(), ENNHit::Compare);
  sort(Down.begin(), Down.end(), ENNHit::Compare);
 
  // save local-out indices correspondece
  vector<Int_t> outIndicesUp;    // indices in HitArray indexed by index in Up
  vector<Int_t> outIndicesDown;  // indices in HitArray indexed by index in Down
  outIndicesUp.resize(Up.size());
  outIndicesDown.resize(Down.size());
  for (THitIndex k = 0; k < Up.size(); k++) {
    Up[k].localIndex = k;
    outIndicesUp[k]  = Up[k].outIndex;
  }
  for (THitIndex k = 0; k < Down.size(); k++) {
    Down[k].localIndex = k;
    outIndicesDown[k]  = Down[k].outIndex;
  }
 
#ifdef PRINT_TIMING
  GetTimer("Sort").Stop();
#endif  // PRINT_TIMING
 
#ifdef PRINT_TIMING
  GetTimer("Repack").Start(0);
#endif  // PRINT_TIMING
 
  // save local-out indices correspondece
  nsL1vector<ENNHitV>::TSimd UpV;
  nsL1vector<ENNHitV>::TSimd DownV;
  UpV.resize((Up.size() + fvecLen - 1) / fvecLen);
  DownV.resize((Down.size() + fvecLen - 1) / fvecLen);
  for (THitIndex k = 0; k < Up.size(); k++) {
    int k_4 = k % fvecLen, k_V = k / fvecLen;
    ENNHitV& hits = UpV[k_V];  // TODO change on ENNHitV
    hits.CopyHit(Up[k], k_4);
  }
  for (THitIndex k = 0; k < Down.size(); k++) {
    int k_4 = k % fvecLen, k_V = k / fvecLen;
    ENNHitV& hits = DownV[k_V];
    hits.CopyHit(Down[k], k_4);
  }
 
#ifdef PRINT_TIMING
  GetTimer("Repack").Stop();
#endif  // PRINT_TIMING
 
  vector<ENNRing> R;
 
  float HitSize         = .5;
  THitIndex MinRingHits = 5;
  float RMin            = 2.;
  float RMax            = 5.;
 
  TStopwatch timer;
 
  ENNRingFinder(Up.size(), UpV, R, HitSize, MinRingHits, RMin, RMax);
  ENNRingFinder(Down.size(), DownV, R, HitSize, MinRingHits, RMin, RMax);
 
  timer.Stop();
 
#ifdef PRINT_TIMING
  GetTimer("Copy data").Start(0);
#endif  // PRINT_TIMING
  Int_t NRings = 0;
  for (vector<ENNRing>::iterator i = R.begin(); i != R.end(); ++i) {
    if (!i->on) continue;
    new ((*RingArray)[NRings]) CbmRichRing();
    CbmRichRing* ring = L1_DYNAMIC_CAST<CbmRichRing*>(RingArray->At(NRings));
    NRings++;
    ring->SetCenterX(i->x);
    ring->SetCenterY(i->y);
    ring->SetRadius(i->r);
    ring->SetChi2(i->chi2);
    const THitIndex NHits = i->localIHits.size();
    for (THitIndex j = 0; j < NHits; j++) {
      if (i->y > 0.)
        ring->AddHit(outIndicesUp[i->localIHits[j]]);
      else
        ring->AddHit(outIndicesDown[i->localIHits[j]]);
    }
  }
#ifdef PRINT_TIMING
  GetTimer("Copy data").Stop();
#endif  // PRINT_TIMING
 
  fNEvents++;
  fRecoTime += timer.RealTime();
  cout.setf(ios::fixed);
  cout.setf(ios::showpoint);
  cout.precision(4);
  cout << "L1ENN Reco time stat/this [ms]  : " << fRecoTime * 1000. / fNEvents
       << "/" << timer.RealTime() * 1000. << endl;
 
#ifdef PRINT_TIMING
  static float timeStat[NTimers];
  static bool firstCall = 1;
  if (firstCall) {
    for (int i = 0; i < NTimers; i++) {
      timeStat[i] = 0;
    }
    firstCall = 0;
  }
 
  cout.width(30);
  cout << "Timing statEvents / curEvent[ms] | statEvents[%]  : " << endl;
  cout.flags(ios::left | ios::oct | ios::showpoint | ios::fixed);
  cout.precision(3);
  for (int i = 0; i < NTimers; i++) {
    timeStat[i] += fTimers[i].RealTime();
    cout.width(30);
    cout << fTimersNames[i] << " timer = " << timeStat[i] * 1000. / fNEvents
         << " / " << fTimers[i].RealTime() * 1000. << " | "
         << timeStat[i] / timeStat[1] * 100 << endl;
  }
#endif  // PRINT_TIMING
  return NRings;
}
 
 
void CbmL1RichENNRingFinderParallel::ENNRingFinder(
  const int NHits,
  nsL1vector<ENNHitV>::TSimd& HitsV,
  vector<ENNRing>& Rings,
  float HitSize,
  THitIndex MinRingHits,
  fvec /*RMin*/,
  fvec RMax) {
#ifdef PRINT_TIMING
  GetTimer("All").Start(0);
#endif  // PRINT_TIMING
  // INITIALIZATION
 
  const fvec MinRingHits_v      = MinRingHits;
  const fvec Rejection          = .5;
  const float ShadowSize        = HitSize / 4;
  const int StartHitMaxQuality  = 15;
  const int SearchHitMaxQuality = 100;  // TODO DELME
 
  const fvec R2MinCut = 3. * 3., R2MaxCut = 7. * 7.;
  //   const fvec R2Min = RMin*RMin, R2Max = RMax*RMax;
  const fvec HitSize2 = 2 * HitSize;
  //   const fvec HitSize4 = 4 * HitSize;
  const fvec HitSizeSq_v = HitSize * HitSize;
  const float HitSizeSq  = HitSizeSq_v[0];
  //   const fvec HitSizeSq4 = HitSize2 * HitSize2;
  const float AreaSize  = 2 * RMax[0] + HitSize;
  const float AreaSize2 = AreaSize * AreaSize;
 
  //   typedef vector<ENNRingHit>::iterator iH;
  //   typedef vector<ENNRingHit*>::iterator iP;
 
  THitIndex NInRings = Rings.size();
 
  // MAIN LOOP OVER HITS
#ifdef PRINT_TIMING
  GetTimer("Ring finding").Start(0);
#endif  // PRINT_TIMING
  nsL1vector<ENNSearchHitV>::TSimd SearchArea;
  nsL1vector<ENNHitV>::TSimd PickUpArea;
  const int MaxAreaHits = 200;  // TODO take into account NHits
 
  SearchArea.resize(MaxAreaHits, ENNSearchHitV());
  PickUpArea.resize(MaxAreaHits);
 
  // TODO 1
#if 0
  nsL1vector<ENNRingV>::TSimd rings_tmp; // simd hits for tmp store
  rings_tmp.resize(100); // TODO use NRings
  int nRings_tmp = 0;
#endif
 
  //  ENNRingHit* ileft = &(Hits[0]), *iright = ileft;//, i_main = ileft;
  int ileft = 0, iright = ileft;
  // int ileft[fvecLen] = {0, 0, 0, 0};
  // int iright[fvecLen] = {0, 0, 0, 0};
 
  THitIndex i_mains[fvecLen] = {0};
 
  THitIndex i_main_array
    [NHits];  // need for proceed in paralled almost independent areas
  for (THitIndex i = 0; i < NHits; i++) {
    i_main_array[i] = i;  // TODO
  }
 
 
  for (THitIndex ii_main = 0; ii_main < NHits;) {
 
#ifdef PRINT_TIMING
    GetTimer("Ring finding: Prepare data").Start(0);
#endif  // PRINT_TIMING
 
 
    fvec S0, S1, S2, S3, S4, S5, S6, S7;
    fvec X = 0, Y = 0, R = 0, R2 = 0;
 
    fvec validRing = (fvec(0) <= fvec(0));  // mask of the valid rings
 
    fvec SearchAreaSize = 0;  // number of hits to fit and search ring
    fvec PickUpAreaSize = 0;
 
    for (int i_4 = 0; (i_4 < fvecLen) && (ii_main < NHits); ii_main++) {
      const THitIndex i_main = i_main_array[ii_main];
      const int i_main_4 = i_main % fvecLen, i_main_V = i_main / fvecLen;
      ENNHitV* i = &HitsV[i_main_V];
      if (i->quality[i_main_4] >= StartHitMaxQuality)
        continue;  // already found hit
 
      i_mains[i_4] = i_main;
 
      float left  = i->x[i_main_4] - AreaSize;
      float right = i->x[i_main_4] + AreaSize;
 
      // while(                          (HitsV[ileft[i_4]/fvecLen] .x[ileft[i_4]%fvecLen] < left ) ) ++ileft[i_4];
      // while( (iright[i_4] < NHits) && (HitsV[iright[i_4]/fvecLen].x[ileft[i_4]%fvecLen] < right) ) ++iright[i_4];
      // for( int j = ileft[i_4]; j < iright[i_4]; ++j ){
      while ((HitsV[ileft / fvecLen].x[ileft % fvecLen] < left))
        ++ileft;  // TODO SIMDize
      while ((iright < NHits)
             && (HitsV[iright / fvecLen].x[ileft % fvecLen] < right))
        ++iright;
 
      for (int j = ileft; j < iright; ++j) {
        const int j_4 = j % fvecLen, j_V = j / fvecLen;
        ENNSearchHitV& sHit = SearchArea[int(SearchAreaSize[i_4])];
        sHit.CopyHit(HitsV[j_V], j_4, i_4);
 
        sHit.ly[i_4] = sHit.y[i_4] - i->y[i_main_4];
        //      if( fabs(sHit.ly) > AreaSize ) continue;
        sHit.lx[i_4] = sHit.x[i_4] - i->x[i_main_4];
        sHit.S0      = sHit.lx * sHit.lx;
        sHit.S1      = sHit.ly * sHit.ly;
        sHit.lr2     = sHit.S0 + sHit.S1;
        // if(( sHit.lr2 > AreaSize2 ) || ( j == i )) continue; // no difference in speed
        if (sHit.lr2[i_4] > AreaSize2) continue;
        if (ISUNLIKELY(j == i_main)) continue;
 
        if (sHit.quality[i_4]
            >= SearchHitMaxQuality) {  // CHECKME do we really need PickUpArea
          PickUpArea[static_cast<int>(PickUpAreaSize[i_4]++)].CopyHit(
            HitsV[j_V], j_4, i_4);
          continue;
        }
 
        SearchAreaSize[i_4]++;
      }  // hits
 
      // if(  NAreaHits+1 < MinRingHits  ) continue;
      if (SearchAreaSize[i_4] < MinRingHits - 1) continue;
      i_4++;
    }  // i_4
 
    ENNHitV iHit;
    int MaxSearchAreaSize = 0;
    int MaxPickUpAreaSize = 0;
    for (int i_4 = 0; i_4 < fvecLen; i_4++) {
      iHit.CopyHit(HitsV[i_mains[i_4] / fvecLen], i_mains[i_4] % fvecLen, i_4);
      MaxSearchAreaSize = (MaxSearchAreaSize < SearchAreaSize[i_4])
                            ? int(SearchAreaSize[i_4])
                            : MaxSearchAreaSize;
      MaxPickUpAreaSize = (MaxPickUpAreaSize < PickUpAreaSize[i_4])
                            ? int(PickUpAreaSize[i_4])
                            : MaxPickUpAreaSize;
    }
#ifdef PRINT_TIMING
    GetTimer("Ring finding: Prepare data").Stop();
    //      assert( MaxSearchAreaSize[i_4] <= MaxAreaHits );
#endif  // PRINT_TIMING
 
#ifdef PRINT_TIMING
    GetTimer("Ring finding: Init of param").Start(0);
#endif  // PRINT_TIMING
    for (int isa = 0; isa < MaxSearchAreaSize;
         isa++) {  // TODO don't work w\o this because of nan in wights
      ENNSearchHitV& sHit = SearchArea[isa];
      const fvec validHit = (fvec(isa) < SearchAreaSize) & validRing;
      sHit.lx             = if3(validHit, sHit.lx, 0);
      sHit.ly             = if3(validHit, sHit.ly, 0);
      sHit.lr2            = if3(validHit, sHit.lr2, 0);
    }
 
    // initialize hits in the search area
    fvec Dmax = 0.;
    S0 = S1 = S2 = S3 = S4 = 0.;
    for (int ih = 0; ih < MaxSearchAreaSize; ih++) {
      ENNSearchHitV& sHit = SearchArea[ih];
      const fvec validHit = (fvec(ih) < SearchAreaSize) & validRing;
 
      fvec& lr2 = sHit.lr2;
      fvec lr   = sqrt(lr2);
      Dmax      = if3((lr > Dmax) & validHit, lr, Dmax);
 
      sHit.S2 = sHit.lx * sHit.ly;
      sHit.S3 = sHit.lx * lr2;
      sHit.S4 = sHit.ly * lr2;
      sHit.C  = -lr * .5;
 
      const fvec w  = if3(validHit, if3(lr > fvec(1.E-4), 1. / lr, 1), 0);
      const fvec w2 = w * w;
      sHit.Cx       = w * sHit.lx;
      sHit.Cy       = w * sHit.ly;
      S0 += w2 * sHit.S0;
      S1 += w2 * sHit.S1;
      S2 += w2 * sHit.S2;
      S3 += w2 * sHit.S3;
      S4 += w2 * sHit.S4;
    }
 
    // end of initialization of the search area
#ifdef PRINT_TIMING
    GetTimer("Ring finding: Init of param").Stop();
    //      assert( MaxSearchAreaSize[i_4] <= MaxAreaHits );
#endif  // PRINT_TIMING
 
#ifdef PRINT_TIMING
    // loop for minimization of E and noise rejection
    GetTimer("Ring finding: Hit selection").Start(0);
#endif  // PRINT_TIMING
    do {
      // calculate parameters of the current ring
      fvec tmp = S0 * S1 - S2 * S2;
 
      //      if( fabs(tmp) < 1.E-10 ) break; // CHECKME
      tmp = 0.5 / tmp;
      X   = (S3 * S1 - S4 * S2) * tmp;
      Y   = (S4 * S0 - S3 * S2) * tmp;
      R2  = X * X + Y * Y;
      R   = sqrt(R2);
 
      const fvec Dcut = Dmax * Rejection;  // cut for noise hits
        //      float RingCut   = HitSize4 * R ; // closeness
      S0 = S1 = S2 = S3 = S4 = 0.0;
      Dmax                   = -1.;
 
      for (THitIndex ih = 0; ih < MaxSearchAreaSize; ih++) {
        const fvec validHit = (fvec(ih) < SearchAreaSize) & validRing;
        //        ENNHit *j = &(SearchArea[ih]);
        ENNSearchHitV& sHit = SearchArea[ih];
        const fvec dx       = sHit.lx - X;
        const fvec dy       = sHit.ly - Y;
        const fvec d        = fabs(sqrt(dx * dx + dy * dy) - R);
        const fvec dSq      = d * d;
        sHit.on_ring        = (d <= HitSize) & validHit;
        const fvec dp =
          if3(sHit.on_ring, -1, fabs(sHit.C + sHit.Cx * X + sHit.Cy * Y));
        Dmax = if3(((dp <= Dcut) & (dp > Dmax)), dp, Dmax);
 
        fvec w = if3((sHit.on_ring),
                     1. / (HitSizeSq_v + fabs(dSq)),
                     1. / (1.e-5 + fabs(dSq)));
        w      = if3((dp <= Dcut) & validHit, w, 0);
        S0 += w * sHit.S0;
        S1 += w * sHit.S1;
        S2 += w * sHit.S2;
        S3 += w * sHit.S3;
        S4 += w * sHit.S4;
      }
 
    } while (NotEmpty(Dmax > fvec(0.)));
 
 
#ifdef PRINT_TIMING
    GetTimer("Ring finding: Hit selection").Stop();
    // store the ring if it is found
 
    GetTimer("Ring finding: Final fit").Start(0);
#endif  // PRINT_TIMING
 
    fvec NRingHits = 1;
    {  // final fit of 3 parameters (X,Y,R)
      S0 = S1 = S2 = S3 = S4 = S5 = S6 = S7 = 0.0;
      for (THitIndex ih = 0; ih < MaxSearchAreaSize; ih++) {
        ENNSearchHitV& sHit = SearchArea[ih];
        const fvec w        = bool2int(SearchArea[ih].on_ring);
        S0 += w * sHit.S0;
        S1 += w * sHit.S1;
        S2 += w * sHit.S2;
        S3 += w * sHit.S3;
        S4 += w * sHit.S4;
        S5 += w * sHit.lx;
        S6 += w * sHit.ly;
        S7 += w * sHit.lr2;
        NRingHits += w;
      }
      const fvec s0 = S6 * S0 - S2 * S5;
      const fvec s1 = S0 * S1 - S2 * S2;
      const fvec s2 = S0 * S4 - S2 * S3;
      //        if( fabs(s0) < 1.E-6 || fabs(s1) < 1.E-6 ) continue; // CHECKME
      const fvec tmp = s1 * (S5 * S5 - NRingHits * S0) + s0 * s0;
      const fvec A   = ((S0 * S7 - S3 * S5) * s1 - s2 * s0) / tmp;
      Y              = (s2 + s0 * A) / s1 * 0.5;
      X              = (S3 + S5 * A - S2 * Y * 2) / S0 * 0.5;
      R2             = X * X + Y * Y - A;
      // if( R2 < 0 ) continue; // will be rejected letter by R2 > R2Min
      R = sqrt(R2);
    }  // end of the final fit
 
    validRing = !((NRingHits < MinRingHits_v) | (R2 > R2MaxCut)
                  | (R2 < R2MinCut));  // ghost suppresion // TODO constants
//    cout << validRing << endl;
#ifdef PRINT_TIMING
    GetTimer("Ring finding: Final fit").Stop();
#endif  // PRINT_TIMING
 
 
#ifdef PRINT_TIMING
    GetTimer("Ring finding: Store ring").Start(0);
#endif  // PRINT_TIMING
    if (ISUNLIKELY(Empty(validRing))) continue;
 
#if 0   // TODO 1
    {
      ENNRingV &ringV = rings_tmp[nRings_tmp++];
    
      ringV.localIHits.reserve(25);
      ringV.x = iHit.x + X;
      ringV.y = iHit.y + Y;
      ringV.r = R;
      ringV.localIHits.push_back(iHit.localIndex);
      ringV.NHits = 1;
      ringV.chi2  = 0;
      nsL1vector<fvec>::TSimd Shadow; // save hit indices of hits, which's quality will be changed
      Shadow.reserve(25);
      Shadow.push_back(iHit.localIndex);
      for( THitIndex ih = 0; ih < MaxSearchAreaSize; ih++){
        fvec validHit = ih < SearchAreaSize;
      
        ENNSearchHitV &sHit = SearchArea[ih];
        const fvec dx = sHit.x - ringV.x;
        const fvec dy = sHit.y - ringV.y;
        const fvec d = fabs( sqrt(dx*dx+dy*dy) - ringV.r );
        validHit = validHit & ( d <= HitSize );
        ringV.chi2 += d*d;
        ringV.localIHits.push_back( if3( validHit, sHit.localIndex, -1 ) );
        ringV.NHits += bool2int(validHit);
        validHit = validHit & ( d <= ShadowSize ); // TODO check *4
        if ( Empty (validHit) ) continue; // CHECKME
        Shadow.push_back( if3( validHit, sHit.localIndex, -1 ) );
      }
      for( int ipu = 0; ipu < MaxPickUpAreaSize; ipu++ ) {
        fvec validHit = ipu < PickUpAreaSize;
      
        ENNHitV &puHit = PickUpArea[ipu];
        const fvec dx = puHit.x - ringV.x;
        const fvec dy = puHit.y - ringV.y;
        const fvec d = fabs( sqrt(dx*dx+dy*dy) - ringV.r );
        validHit = validHit & ( d <= HitSize );
        if ( Empty (validHit) ) continue;
        ringV.chi2 += d*d;
        ringV.localIHits.push_back( if3( validHit, puHit.localIndex, -1 ) );
        ringV.NHits += bool2int(validHit);
        validHit = validHit & ( d <= ShadowSize ); // TODO check *4
        if ( Empty (validHit) ) continue; // CHECKME
        Shadow.push_back( if3( validHit, puHit.localIndex, -1 ) );
      }
 
      ringV.chi2 = ringV.chi2 / (( ringV.NHits - 3)*HitSizeSq);
      const fvec quality = ringV.NHits;
 
      
        // for( iP j = ringV.Hits.begin(); j != ringV.Hits.end(); ++j ){
        //   if( (*j)->quality<quality ) (*j)->quality = quality;
        // } 
   
        //quality *= ShadowOpacity;
      for( int i_4 = 0; (i_4 < fvecLen); i_4++) { 
        const int NShadow = Shadow.size();
        for( int is = 0; is < NShadow; is++ ) { // CHECKME change loops to speed up?
          cout << i_4 << Shadow[is] << endl;
           float ih_f = (Shadow[is])[i_4]-200;
          if (ih_f == -1) continue;
           int ih = static_cast<int>(ih_f);  // TODO ! problem in conversion...
           float ih_f2 =  static_cast<float>(ih);
          cout << ih_f << " " << ih << " " << ih_f2 << " " << ih%fvecLen << " " << ih/fvecLen << endl;
 
          const THitIndex ih_4 = ih%fvecLen;
          ENNHitV & hitV = HitsV[ih/fvecLen];
 
            //          hitV.quality[ih_4] = ( hitV.quality[ih_4] < quality[i_4] ) ? quality[i_4] : hitV.quality[ih_4];
            //        shHit->quality = if3( shHit->quality < quality, quality, shHit->quality );
        }
      } // i_4
    }
#endif  // 0
 
    for (int i_4 = 0; (i_4 < fvecLen); i_4++) {
      //      if( NRingHits < MinRingHits || R2 > R2Max || R2 < R2Min ) continue;
 
      if (/*ISUNLIKELY*/ (!validRing[i_4])) continue;
 
      ENNRing voidRing;
      Rings.push_back(voidRing);
      ENNRing& ring = Rings.back();
      ring.localIHits.reserve(15);
      ring.x = iHit.x[i_4] + X[i_4];
      ring.y = iHit.y[i_4] + Y[i_4];
      ring.r = R[i_4];
      ring.localIHits.push_back(static_cast<THitIndex>(iHit.localIndex[i_4]));
      ring.NHits = 1;
      ring.chi2  = 0;
      vector<THitIndex>
        Shadow;  // save hit indices of hits, which's quality will be changed
      Shadow.reserve(15);
      Shadow.push_back(static_cast<THitIndex>(iHit.localIndex[i_4]));
      for (THitIndex ih = 0; ih < SearchAreaSize[i_4]; ih++) {
        ENNSearchHitV& sHit = SearchArea[ih];
        const float dx      = sHit.x[i_4] - ring.x;
        const float dy      = sHit.y[i_4] - ring.y;
        const float d       = fabs(sqrt(dx * dx + dy * dy) - R[i_4]);
        if (ISLIKELY(d <= HitSize)) {
          ring.chi2 += d * d;
          ring.localIHits.push_back(int(sHit.localIndex[i_4]));
          ring.NHits++;
          if (d <= ShadowSize)
            Shadow.push_back(static_cast<THitIndex>(sHit.localIndex[i_4]));
        }
      }
      for (int ipu = 0; ipu < PickUpAreaSize[i_4]; ipu++) {
        ENNHitV& puHit = PickUpArea[ipu];
        const float dx = puHit.x[i_4] - ring.x;
        const float dy = puHit.y[i_4] - ring.y;
        const float d  = fabs(sqrt(dx * dx + dy * dy) - ring.r);
        if (ISLIKELY(d <= HitSize)) {
          ring.chi2 += d * d;
          ring.localIHits.push_back(
            static_cast<THitIndex>(puHit.localIndex[i_4]));
          ring.NHits++;
          if (d <= ShadowSize)
            Shadow.push_back(static_cast<THitIndex>(puHit.localIndex[i_4]));
        }
      }
      if (ISUNLIKELY(ring.NHits < MinRingHits)) {
        Rings.pop_back();
        continue;
      }
      ring.chi2   = ring.chi2 / ((ring.NHits - 3) * HitSizeSq);
      int quality = ring.NHits;
 
      // for( iP j = ring.Hits.begin(); j != ring.Hits.end(); ++j ){
      //   if( (*j)->quality<quality ) (*j)->quality = quality;
      // }
 
      //quality *= ShadowOpacity;
      const int NShadow = Shadow.size();
      for (int is = 0; is < NShadow; is++) {
        const THitIndex ih   = Shadow[is];
        const THitIndex ih_4 = ih % fvecLen;
        ENNHitV& hitV        = HitsV[ih / fvecLen];
 
        hitV.quality[ih_4] =
          (hitV.quality[ih_4] < quality) ? quality : hitV.quality[ih_4];
        //        shHit->quality = if3( shHit->quality < quality, quality, shHit->quality );
      }
    }  // i_4
 
 
#ifdef PRINT_TIMING
    GetTimer("Ring finding: Store ring").Stop();
#endif  // PRINT_TIMING
  }     // i_main
#ifdef PRINT_TIMING
  GetTimer("Ring finding").Stop();
 
  // SELECTION OF RINGS
  GetTimer("Selection").Start(0);
#endif  // PRINT_TIMING
  typedef vector<ENNRing>::iterator iR;
  iR Rbeg = Rings.begin() + NInRings;
  iR Rend = Rings.end();
 
//#define NEW_SELECTION
#ifdef NEW_SELECTION  // TODO optimize. at the moment just creates additional ghosts
 
  sort(Rings.begin() + NInRings, Rings.end(), ENNRing::CompareENNHRings);
 
  const int NHitsV = HitsV.size();
  for (int ih = 0; ih < NHitsV; ih++)
    HitsV[ih].quality = 0.;
 
  for (iR ir = Rbeg; ir != Rend; ++ir) {
    ir->on   = 0;
    ir->NOwn = ir->NHits;
  }
 
  for (iR i = Rbeg; i != Rend; ++i) {
    if (i->skip) continue;
    for (iR j = i + 1; j != Rend; ++j) {
      if (j->skip) continue;
 
      float dist = j->r + i->r + HitSize2[0];
      float distCentr =
        sqrt((j->x - i->x) * (j->x - i->x) + (j->y - i->y) * (j->y - i->y));
      if (distCentr > dist) continue;
      Int_t NOverlaped = 0;
 
      const THitIndex maxI = i->localIHits.size();
      const THitIndex maxJ = j->localIHits.size();
      for (THitIndex n = 0; n < maxI; n++)
        for (THitIndex m = 0; m < maxJ; ++m)
          if (i->localIHits[n] == j->localIHits[m]) NOverlaped++;
      ENNRing* BigRing = 0;
      if (NOverlaped > 0.7 * maxI) BigRing = &(*i);
      if (NOverlaped > 0.7 * maxJ) BigRing = &(*j);
      if (BigRing != 0) {
        std::vector<THitIndex> newIndices;
        for (THitIndex n = 0; n < maxI; n++) {
          bool IsNew = 1;
          for (THitIndex m = 0; m < maxJ; ++m)
            if (i->localIHits[n] == j->localIHits[m]) IsNew = 0;
          if (IsNew) newIndices.push_back(i->localIHits[n]);
        }
        if (maxI > maxJ) {
          j->x = i->x;
          j->y = i->y;
          j->r = i->r;
        }
        const THitIndex newISize = newIndices.size();
        for (THitIndex in = 0; in < newISize; in++)
          j->localIHits.push_back(newIndices[in]);
        i->skip = 1;
        i->on   = 0;
        break;
      }
 
    }  // j
  }    // i
 
 
  for (iR i = Rbeg; i != Rend; ++i) {
    if (!(i->skip)) {
      i->on = 1;
      float S0, S1, S2, S3, S4, S5, S6, S7, X, Y, R;
      S0 = S1 = S2 = S3 = S4 = S5 = S6 = S7 = 0.0;
 
 
      const THitIndex firstIh = i->localIHits[0];
      const ENNHitV& firstHit = HitsV[firstIh / fvecLen];
      const int firstIh_4     = firstIh % fvecLen;
      const THitIndex maxI    = i->localIHits.size();
 
      vector<ENNSearchHitV> shits;
      shits.resize(maxI);
      for (THitIndex iih = 0; iih < maxI; iih++) {
        const THitIndex ih  = i->localIHits[iih];
        const ENNHitV& hit  = HitsV[ih / fvecLen];
        const int ih_4      = ih % fvecLen;
        ENNSearchHitV& shit = shits[iih];
 
        shit.ly[0]  = hit.y[ih_4] - firstHit.y[firstIh_4];
        shit.lx[0]  = hit.x[ih_4] - firstHit.x[firstIh_4];
        shit.S0[0]  = shit.lx[0] * shit.lx[0];
        shit.S1[0]  = shit.ly[0] * shit.ly[0];
        shit.lr2[0] = shit.S0[0] + shit.S1[0];
        float lr2   = shit.lr2[0];
        float lr    = sqrt(lr2);
        shit.S2[0]  = shit.lx[0] * shit.ly[0];
        shit.S3[0]  = shit.lx[0] * lr2;
        shit.S4[0]  = shit.ly[0] * lr2;
        shit.C[0]   = -lr * 0.5;
        float w;
        if (lr > 1.E-4)
          w = 1. / lr;
        else
          w = 1.;
        shit.Cx[0] = w * shit.lx[0];
        shit.Cy[0] = w * shit.ly[0];
      }
      float Dmax = -1.;
 
      X               = i->x - firstHit.x[firstIh_4];
      Y               = i->y - firstHit.y[firstIh_4];
      R               = i->r;
      int search_stop = 0;
      do {  // final fit of 3 parameters (X,Y,R)
        float Dcut = Dmax * Rejection[0];
        int n      = 0;
        S0 = S1 = S2 = S3 = S4 = S5 = S6 = S7 = 0.0;
 
        for (THitIndex ih = 0; ih < maxI; ih++) {
          ENNSearchHitV& shit = shits[ih];
 
          float dx        = shit.lx[0] - X;
          float dy        = shit.ly[0] - Y;
          float d         = fabs(sqrt(dx * dx + dy * dy) - R);
          shit.on_ring[0] = (d <= HitSize2[0]);
          float w;
          if (shit.on_ring[0]) {
            n++;
            w = 1;
          } else {
            float dp = fabs(shit.C[0] + shit.Cx[0] * X + shit.Cy[0] * Y);
            if (dp > Dcut) continue;
            if (dp > Dmax) Dmax = dp;
            n++;
            w = 10. / (1.e-5 + fabs(d * d));
          }
 
          S0 += w * shit.S0[0];
          S1 += w * shit.S1[0];
          S2 += w * shit.S2[0];
          S3 += w * shit.S3[0];
          S4 += w * shit.S4[0];
          S5 += w * shit.lx[0];
          S6 += w * shit.ly[0];
          S7 += w * shit.lr2[0];
        }  // ih
 
        float s0 = S6 * S0 - S2 * S5;
        float s1 = S0 * S1 - S2 * S2;
        float s2 = S0 * S4 - S2 * S3;
 
        if (fabs(s0) < 1.E-6 || fabs(s1) < 1.E-6) continue;
        float tmp = s1 * (S5 * S5 - n * S0) + s0 * s0;
        float A   = ((S0 * S7 - S3 * S5) * s1 - s2 * s0) / tmp;
        Y         = (s2 + s0 * A) / s1 / 2;
        X         = (S3 + S5 * A - S2 * Y * 2) / S0 / 2;
        float R2  = X * X + Y * Y - A;
 
        if (R2 < 0) continue;
        R = sqrt(R2);
        if (Dmax <= 0) search_stop++;
      } while (search_stop < 2.);
 
      i->r = R;
      i->x = X + firstHit.x[firstIh_4];
      i->y = Y + firstHit.y[firstIh_4];
 
      if (R < 2.5 || R > 7.5) {
        i->on   = 0;
        i->skip = 1;
        continue;
      }
 
      std::vector<THitIndex> newHits;
      i->chi2 = 0;
 
      for (THitIndex iih = 0; iih < maxI; iih++) {
        const THitIndex ih  = i->localIHits[iih];
        const ENNHitV& hit  = HitsV[ih / fvecLen];
        ENNSearchHitV& shit = shits[iih];
        const int ih_4      = ih % fvecLen;
 
        float dx           = hit.x[ih_4] - i->x;
        float dy           = hit.y[ih_4] - i->y;
        float d            = fabs(sqrt(dx * dx + dy * dy) - i->r);
        shit.on_ring[ih_4] = (d <= HitSize);
        if (shit.on_ring[ih_4]) {
          newHits.push_back(ih);
          i->chi2 += d * d;
        }
      }
 
      i->localIHits.clear();
      i->localIHits = newHits;
      i->NHits      = i->localIHits.size();
      i->NOwn       = i->NHits;
 
      if (i->localIHits.size() < MinRingHits) {
        i->on   = 0;
        i->skip = 1;
        continue;
      }
      i->chi2 =
        i->chi2 / ((i->localIHits.size() - 3) * HitSize * HitSize);  //.3/.3;
    }
  }
 
  sort(Rings.begin() + NInRings, Rings.end(), ENNRing::CompareENNHRings);
  iR best;
  for (iR i = Rbeg; i != Rend; ++i) {
    i->on = 0;
    if (i->skip) continue;
    if (i->NOwn > 5) {
      best                 = i;
      const THitIndex maxI = i->localIHits.size();
      for (THitIndex n = 0; n < maxI; n++) {
        const THitIndex ih = i->localIHits[n];
        ENNHitV& hit       = HitsV[ih / fvecLen];
        const int ih_4     = ih % fvecLen;
        hit.quality[ih_4]  = 1;
      }
      for (iR j = i + 1; j != Rend; ++j) {
        if (i->skip) continue;
        float dist = j->r + best->r + HitSize2[0];
        if (fabs(j->x - best->x) > dist || fabs(j->y - best->y) > dist)
          continue;
        j->NOwn              = 0;
        const THitIndex maxJ = j->localIHits.size();
        for (THitIndex m = 0; m < maxJ; m++) {
          const THitIndex ihm = j->localIHits[m];
          const ENNHitV& hitm = HitsV[ihm / fvecLen];
          if (hitm.quality[ihm % fvecLen] == 0) j->NOwn++;
        }
      }
      i->on   = 1;
      i->skip = 1;
    } else
      i->skip = 1;
  }
#else   // NEW_SELECTION
 
  const int NHitsV = HitsV.size();
  for (int ih = 0; ih < NHitsV; ih++)
    HitsV[ih].quality = 0.;
  for (iR ir = Rbeg; ir != Rend; ++ir) {
    ir->on   = 0;
    ir->NOwn = ir->NHits;
    ir->skip = ((ir->NHits < MinRingHits) || (ir->NHits <= 6 && ir->chi2 > .3));
  }
 
  do {
    iR best           = Rend;
    THitIndex bestOwn = 0;
    float bestChi2    = 1.E20;
    for (iR ir = Rbeg; ir != Rend; ++ir) {
      // const bool skip = ir->skip || ( ( ir->NOwn  < 1.0*MinRingHits ) ||
      //                                 ( ir->NHits < 10 && ir->NOwn < 1.00*ir->NHits ) ||
      //                                 ( ir->NOwn  < .60*ir->NHits ) );
      if (ir->skip) continue;  // faster with if
      const bool skip = ((ir->NOwn < 1.0 * MinRingHits)
                         || (ir->NHits < 10 && ir->NOwn < 1.00 * ir->NHits)
                         || (ir->NOwn < .60 * ir->NHits));
      ir->skip        = skip;
      const bool isBetter =
        !skip
        & ((ir->NOwn > 1.2 * bestOwn)
           || (ir->NOwn >= 1. * bestOwn && ir->chi2 < bestChi2));
      bestOwn  = (isBetter) ? ir->NOwn : bestOwn;  //Hits;
      bestChi2 = (isBetter) ? ir->chi2 : bestChi2;
      best     = (isBetter) ? ir : best;
    }
    if (best == Rend) break;
    best->skip                = 1;
    best->on                  = 1;
    const THitIndex NHitsBest = best->localIHits.size();
    for (THitIndex iih = 0; iih < NHitsBest; iih++) {
      const THitIndex ih                        = best->localIHits[iih];
      HitsV[ih / fvecLen].quality[ih % fvecLen] = 1;
    }
    for (iR ir = Rbeg; ir != Rend; ++ir) {
      if (ir->skip) continue;
      float dist = ir->r + best->r + HitSize2[0];
      if (fabs(ir->x - best->x) > dist || fabs(ir->y - best->y) > dist)
        continue;
      ir->NOwn                 = 0;
      const THitIndex NHitsCur = ir->localIHits.size();
      for (THitIndex iih = 0; iih < NHitsCur; iih++) {
        const THitIndex ih = ir->localIHits[iih];
        ir->NOwn += (HitsV[ih / fvecLen].quality[ih % fvecLen] == 0);
      }
    }
  } while (1);
#endif  // else NEW_SELECTION
 
#ifdef PRINT_TIMING
  GetTimer("Selection").Stop();
  GetTimer("All").Stop();
#endif  // PRINT_TIMING
}
 
TStopwatch& CbmL1RichENNRingFinderParallel::GetTimer(TString t) {
  int i = 0;
  for (; (i < NTimers) && (fTimersNames[i] != t); i++)
    ;
  assert(i < NTimers);
  return fTimers[i];
}