TopDownSelector.cc

Go to the documentation of this file.

#include "TopDownSelector.h"
#include <utl/ErrorLogger.h>

#include <evt/Event.h>
#include <evt/ShowerSRecData.h>
#include <evt/ShowerRecData.h>
#include <sevt/EventTrigger.h>
#include <sevt/SEvent.h>
#include <sevt/Station.h>
#include <sevt/StationConstants.h>
#include <sevt/SortCriteria.h>

#include <fwk/CentralConfig.h>
#include <det/Detector.h>
#include <sdet/Station.h>
#include <sdet/SDetector.h>
#include <utl/CoordinateSystemPtr.h>
#include <utl/Branch.h>

#include <TCdasUtil.h>

using namespace sdet;
using namespace std;
using namespace sevt;
using namespace evt;
using namespace fwk;
using namespace utl;
using namespace TopDownSelectorNS;


static bool bFirstIteration;


TopDownSelector::TopDownSelector() :
  fSEvent(0),
  fVerbose(0),
  fRejectNonT4Events(false),
  fU(0),
  fV(0),
  fXCore(0),
  fYCore(0),
  fZCore(0),
  fTCore(0),
  fNMinStations(3),
  fNMaxRemovedStations(4)
{ }


VModule::ResultFlag
TopDownSelector::Init()
{
  INFO("TopDownSelector::Init()");

  Branch tb = CentralConfig::GetInstance()->GetTopBranch("TopDownSelector");
  tb.GetChild("RejectNonT4Events").GetData(fRejectNonT4Events);

  tb.GetChild("MinNumberStations").GetData(fNMinStations);
  tb.GetChild("MaxRejectedStations").GetData(fNMaxRemovedStations);

  tb.GetChild("Verbose").GetData(fVerbose);
  tb.GetChild("BaryPower").GetData(kBaryPower);
  tb.GetChild("TimeResidualTolerance").GetData(kTolRes);
  tb.GetChild("MTimeResidualTolerance").GetData(kTolResM);
  tb.GetChild("IsolatedDistance1").GetData(kIsoDist1);
  tb.GetChild("IsolatedDistance2").GetData(kIsoDist2);
  tb.GetChild("Curv0").GetData(kCurv0);
  tb.GetChild("AreaMin").GetData(kAreaMin);
  tb.GetChild("AreaMax").GetData(kAreaMax);
  tb.GetChild("Ana3RatMin").GetData(kAna3RatMin);
  tb.GetChild("DistMax").GetData(kDistMax);
  tb.GetChild("IsoTime1").GetData(kIsoTime1);
  tb.GetChild("IsoTime2").GetData(kIsoTime2);

  return eSuccess;
}


VModule::ResultFlag
TopDownSelector::Run(evt::Event& event)
{
  INFO("TopDownSelector::Run()");

  if (!event.HasSEvent()) {
    INFO("TopDownSelector::Run() => Event has not SEvent");
    return eContinueLoop;
  }

  fSEvent = &event.GetSEvent();

  const bool isFD = fSEvent->GetTrigger().IsFD();
  const bool isSdT4 = Select();

  if (!event.HasRecShower())
    event.MakeRecShower();
  if (!event.GetRecShower().HasSRecShower())
    event.GetRecShower().MakeSRecShower();

  ShowerSRecData& shSRec = event.GetRecShower().GetSRecShower();
  shSRec.SetT4Trigger((isFD * ShowerSRecData::eT4_FD) |
                      (isSdT4 * ShowerSRecData::eT4_4C1));

  if (fRejectNonT4Events && !isSdT4) {
    INFO("Event was not selected.");
    return eContinueLoop;
  }
  return eSuccess;
}


VModule::ResultFlag
TopDownSelector::Finish()
{
  INFO("TopDownSelector::Finish()");
  return eSuccess;
}


// ==============================================
// driving function
// ==============================================

bool
TopDownSelector::Select()
{
  InitializeEvent();

  const unsigned int ns = fStations.size();
  const unsigned int k = ((ns - fNMinStations) < fNMaxRemovedStations) ? ns - fNMinStations : fNMaxRemovedStations;

  // trivial selection
  if (ns < fNMinStations)
    return false;

  if (PatternRecovering()) {
    if (fVerbose > 0)
      cout << " Event Accepted with all stations. Theta ~ "
           << EstimatedZenith("deg") << " deg." << endl;
    return true;
  }

  // "top down" selection
  bFirstIteration = true;
  for (unsigned int i = 1; i <= k; ++i) {
    fIndexes.clear();
    fIndexes.resize(i);
    fIndexes[i-1] = ns;
    if (fVerbose > 1) {
      cout << "  Select: Trying with " << i << " less stations over "<< ns << endl;
      for (unsigned int ik = 0; ik < i; ++ik)
        cout << "   Select: fIndexes[" << ik << "] = "<< fIndexes[ik] << endl;
    }
    if (TryRemovingStations(i)) {
      for (unsigned int j = 0; j < fStations.size(); ++j) {
        if (!fStations[j].candidate)
          fStations[j].event->SetRejected(sevt::StationConstants::eOutOfTime);
      }
      if (fVerbose > 0) {
        cout << " Event accepted with " << i << " rejected. Theta ~ "
             << EstimatedZenith("deg") << " deg." << endl;
      }
      return true;
    }
  }

  bFirstIteration = false;
  for (unsigned int i = 1; i <= k; ++i) {
    fIndexes.clear();
    fIndexes.resize(i);
    fIndexes[i-1] = ns;
    if (fVerbose > 1) {
      cout << "  Select: Trying with " << i << " less stations over "<< ns << endl;
      for (unsigned int ik = 0; ik < i; ++ik)
        cout << "   Select: fIndexes[" << ik << "] = "<< fIndexes[ik] << endl;
    }
    if (TryRemovingStations(i)) {
      for (unsigned int j = 0; j < fStations.size(); ++j) {
        if (!fStations[j].candidate)
          fStations[j].event->SetRejected(sevt::StationConstants::eOutOfTime);
      }
      if (fVerbose > 0) {
        cout << " Event accepted with " << i << " rejected. Theta ~ "
             << EstimatedZenith("deg") << " deg." << endl;
      }
      return true;
    }
  }

  return false;
}


void
TopDownSelector::InitializeEvent()
{
  RemoveEA();
  RemoveDoublets();
  RemoveIsolatedStations();

  fStations.clear();
  const sdet::SDetector& sdetector = det::Detector::GetInstance().GetSDetector();
  fSEvent->SortStations(sevt::ByIncreasingSignal());
  for (sevt::SEvent::CandidateStationIterator sIt = fSEvent->CandidateStationsBegin();
       sIt != fSEvent->CandidateStationsEnd(); ++sIt) {
    StationPairStatus stationPair;
    stationPair.detector = &sdetector.GetStation(*sIt);
    stationPair.event = &(*sIt);
    stationPair.candidate = true;
    fStations.push_back(stationPair);
  }
}


void
TopDownSelector::RemoveEA()
{
  for (sevt::SEvent::CandidateStationIterator sIt = fSEvent->CandidateStationsBegin();
       sIt != fSEvent->CandidateStationsEnd(); ++sIt) {
    if (sIt->GetId() < 100) {
      sIt->SetRejected(sevt::StationConstants::eEngineeringArray);
    }
  }
}


// This one actually removes all off-grid stations (in-fill, triplets, doublets),
// not only doublets. If we want to remove only doublets we need to change IsInFill
// by the corresponding one
void
TopDownSelector::RemoveDoublets()
{
  const sdet::SDetector& sdetector = det::Detector::GetInstance().GetSDetector();

  for (sevt::SEvent::CandidateStationIterator sIt = fSEvent->CandidateStationsBegin();
       sIt != fSEvent->CandidateStationsEnd(); ++sIt) {
    const sdet::Station& station = sdetector.GetStation(*sIt);
    if (!station.IsInGrid()) {
      sIt->SetRejected(sevt::StationConstants::eOffGrid);
    }
  }
}


void
TopDownSelector::RemoveIsolatedStations()
{
  const sdet::SDetector& sdetector = det::Detector::GetInstance().GetSDetector();

  for (sevt::SEvent::CandidateStationIterator sIt_1 = fSEvent->CandidateStationsBegin();
       sIt_1 != fSEvent->CandidateStationsEnd(); ++sIt_1) {
    utl::Point fStPoint1 = sdetector.GetStation(*sIt_1).GetPosition();
    int n1 = 0;
    int n2 = 0;
    TimeStamp t1 = sIt_1->GetRecData().GetSignalStartTime();

    for (sevt::SEvent::CandidateStationIterator sIt_2 = fSEvent->CandidateStationsBegin();
         sIt_2 != fSEvent->CandidateStationsEnd(); ++sIt_2) {
      if (sIt_1->GetId() == sIt_2->GetId())
        continue;

      const utl::Point fStPoint2 = sdetector.GetStation(*sIt_2).GetPosition();
      const double d = ( fStPoint1 - fStPoint2 ).GetMag();
      const TimeStamp t2 = sIt_2->GetRecData().GetSignalStartTime();
      const double dt = fabs((t1-t2).GetInterval());

      if (d < kIsoDist1 && dt < kIsoTime1)
        ++n1;
      if (d < kIsoDist2 && dt < kIsoTime2)
        ++n2;
    }

    if (n1 == 0)
      sIt_1->SetRejected(sevt::StationConstants::eLonely);

    if (n1 && n2 < 2)
      sIt_1->SetRejected(sevt::StationConstants::eLonely);
  }
}


/*void
TopDownSelector::RemoveIsolatedStations()
{
  const sdet::SDetector& sdetector = det::Detector::GetInstance().GetSDetector();
  // discarding very isolated stations
  unsigned int size = fStations.size();
  if (size>=40) kIsoDist1=3400.; // => gain time in case of huge events...
  for (sevt::SEvent::CandidateStationIterator sIt_1 = fSEvent->CandidateStationsBegin(); sIt_1 != fSEvent->CandidateStationsEnd(); ++sIt_1) {
    int n1 = 0;
    int n2 = 0;
    for (sevt::SEvent::CandidateStationIterator sIt_2 = fSEvent->CandidateStationsBegin(); sIt_2 != fSEvent->CandidateStationsEnd(); ++sIt_2) {
      const double d = (sdetector.GetStation(*sIt_1).GetPosition() - sdetector.GetStation(*sIt_2).GetPosition()).GetMag();
      if ( d<kIsoDist1 ) n1++;
      if ( d<kIsoDist2 ) n2++;
    }
    if (n1==0) sIt_1->SetRejected(sevt::StationConstants::eLonely);
    if (n1 && n2<2) sIt_1->SetRejected(sevt::StationConstants::eLonely);
  }
}
*/


// ==============================================
// selection algorithm
// ==============================================

bool
TopDownSelector::PatternRecovering()
{
  unsigned int ns = 0;
  for (unsigned int is = 0; is < fStations.size(); ++is)
    if (fStations[is].candidate)
      ++ns;

  if (ns == 3) {
    if (fNMinStations == 3)
      return Recover3Stations();
    else
      return false;
  }
  if (!EstimateCore())
    return false;
  fIsLastTimingIteration = false;
  if (!IsGoodTimeConfig(eNone))
    return false;
  if (!IsGoodSpaceConfig())
    return false;
  fIsLastTimingIteration = true;
  if (!IsGoodTimeConfig(eAltitude))
    return false;
  if (IsAligned())
    return false;
  return true;
}


bool
TopDownSelector::TryRemovingStations(int n)
{
  if (n < 2) {
    for (unsigned int is = 0; is < fIndexes[0]; ++is) {
      fStations[is].candidate = false;
      if (PatternRecovering())
        return true;
      fStations[is].candidate = true;
    }
  } else {
    for (unsigned int is = n-1; is < fIndexes[n-1]; ++is) {
      fStations[is].candidate = false;
      fIndexes[n-2] = is;
      if (TryRemovingStations(n-1))
        return true;
      fStations[is].candidate = true;
    }
  }
  return false;
}


bool
TopDownSelector::EstimateCore()
{
  const CoordinateSystemPtr siteCS =
    det::Detector::GetInstance().GetSiteCoordinateSystem();
  const sdet::SDetector& sdetector = det::Detector::GetInstance().GetSDetector();

  fBarycenter = Vector(0,0,0, siteCS);
  fBaryTime = 0;

  // the absolute points in time and space
  const Point siteOrigin(0,0,0, siteCS);
  const TimeStamp& eventTime = fSEvent->GetTrigger().GetTime();

  if (fVerbose > 2)
    cout << "   Estimating core..." << endl;

  unsigned int size = fStations.size();
  double sumw = 0;
  for (unsigned int is = 0; is < size; ++is) {
    if (!fStations[is].candidate)
      continue;
    sevt::Station& station = *fStations[is].event;

    const double weight = exp(kBaryPower*log(station.GetRecData().GetTotalSignal()));
    fBarycenter += weight * (sdetector.GetStation(station).GetPosition() - siteOrigin);
    fBaryTime += weight *
      (eventTime - station.GetRecData().GetSignalStartTime()).GetInterval();

    sumw += weight;
  }
  if (sumw == 0) {
    if (fVerbose > 2)
      cout << "    ...FAILED." << endl;
    return false;
  }
  fBarycenter /= sumw;
  fBaryTime /= sumw;
  fTCore = fBaryTime;
  fXCore = fBarycenter.GetX(siteCS);
  fYCore = fBarycenter.GetY(siteCS);
  fZCore = fBarycenter.GetZ(siteCS);

  if (fVerbose > 2)
    cout << "    ... Core(X,Y,Z,t): ("
         << fXCore << " , "
         << fYCore << " , "
         << fZCore << " , "
         << fTCore << ")" << endl;

  return true;
}


double
TopDownSelector::CorrectTimeAltitude(const int is)
{
  const CoordinateSystemPtr siteCS =
    det::Detector::GetInstance().GetSiteCoordinateSystem();
  const double costh = sqrt(max(0., 1 - fU*fU - fV*fV));
  return costh * (fStations[is].detector->GetPosition().GetZ(siteCS) - fZCore) / kSpeedOfLight;
}


double
TopDownSelector::CorrectTimeAltitudeCurv(const int is)
{
  // this is not used right now and probably shouldn't be used (HD)
  const CoordinateSystemPtr siteCS =
    det::Detector::GetInstance().GetSiteCoordinateSystem();
  const double costh = sqrt(max(0., 1 - fU*fU - fV*fV));
  const double dx = fStations[is].detector->GetPosition().GetX(siteCS);
  const double dy = fStations[is].detector->GetPosition().GetY(siteCS);
  const double d2 = dx*dx + dy*dy - pow(fU*dx + fV*dy, 2);
  // WARNING: model of the curvature seems too simple (HD)
  const double curv = kCurv0 * costh;
  return (costh * (fStations[is].detector->GetPosition().GetZ(siteCS) - fZCore) - 0.5*curv*d2) / kSpeedOfLight;
}


bool
TopDownSelector::IsGoodTimeConfig(TimeCorrectionType correctionType)
{
  const CoordinateSystemPtr siteCS =
    det::Detector::GetInstance().GetSiteCoordinateSystem();
  unsigned int is, ns = 0;
  double wgt, dx, dy, dt, sumw, sumx, sumy, sumx2, sumy2, sumxy, sumt, sumxt, sumyt;
  double rxx, ryy, rxy, rx, ry, rr, deter;

  if (fVerbose > 2)
    cout << "   Checking for Good Time Configuration. Correction Type:  " << correctionType << " ..." << endl;

  const TimeStamp& eventTime = fSEvent->GetTrigger().GetTime();
  unsigned int size = fStations.size();
  sumw = sumx = sumy = sumx2 = sumy2 = sumxy = sumt = sumxt = sumyt = 0;
  for (is = 0; is < size; ++is) {
    if (!fStations[is].candidate)
      continue;

    dx = fStations[is].detector->GetPosition().GetX(siteCS) - fXCore;
    dy = fStations[is].detector->GetPosition().GetY(siteCS) - fYCore;

    if (correctionType == TopDownSelector::eAltitude)
      dt = CorrectTimeAltitude(is);
    else if (correctionType == TopDownSelector::eAltitudeCurved)
      dt = CorrectTimeAltitudeCurv(is);
    else if (correctionType == TopDownSelector::eNone)
      dt = 0;
    else {
      dt = 0;
      INFO("Timing Correction Type not recognized. No correction done.");
    }
    dt +=
      (eventTime-fStations[is].event->GetRecData().GetSignalStartTime()).GetInterval() -
      fBaryTime;
    wgt = 1;
    sumw += wgt;
    sumx += wgt * dx;
    sumy += wgt * dy;
    sumx2 += wgt * dx*dx;
    sumy2 += wgt * dy*dy;
    sumxy += wgt * dx * dy;
    sumt += wgt * dt;
    sumxt += wgt * dt * dx;
    sumyt += wgt * dt * dy;
    ++ns;
  }

  rxx = sumw * sumy2 - sumy*sumy;
  ryy = sumw * sumx2 - sumx*sumx;
  rxy = sumx * sumy - sumw * sumxy;
  rx = sumxy * sumy - sumx * sumy2;
  ry = sumxy * sumx - sumy * sumx2;
  rr = sumx2 * sumy2 - sumxy*sumxy;
  deter = sumx2 * rxx + sumxy * rxy + sumx * rx;

  if (deter == 0 || std::isnan(deter))
    return false;
  fU = (rxx * sumxt + rxy * sumyt + rx * sumt) / deter;
  fV = (rxy * sumxt + ryy * sumyt + ry * sumt) / deter;
  fT0 = (rx * sumxt + ry * sumyt + rr * sumt) / deter;
  fU *= -kSpeedOfLight;
  fV *= -kSpeedOfLight;

  if (!fIsLastTimingIteration)
    return true;
  if (fU*fU + fV*fV > 1)
    return false;

  if (ns == 3)
    return true;

  double sumdt2 = 0;
  double costh = sqrt(max(0., 1 - fU*fU - fV*fV));

  if (fVerbose > 3)
    cout << "    Weight sum: " << sumw
         << "\t Residual Max: " << (ns-2)*kTolResM*max(costh, 0.2) << endl;

  fResiduals.clear();
  for (is = 0; is < size; ++is) {
    if (!fStations[is].candidate || !fStations[is].event->HasRecData())
      continue; // trying patterns...
    const double dx = fStations[is].detector->GetPosition().GetX(siteCS) - fXCore;
    const double dy = fStations[is].detector->GetPosition().GetY(siteCS) - fYCore;
    const double dt =
      (eventTime-fStations[is].event->GetRecData().GetSignalStartTime()).GetInterval() -
      fBaryTime - (fT0 -fU*dx/kCdasUtil::CSPEED - fV*dy/kCdasUtil::CSPEED);
    sumdt2 += dt*dt;

    if (fVerbose > 3)
      cout << "    St: " << fStations[is].event->GetId()
           << "\tResidual: " << dt << " --> ";

    if (fabs(dt) > (ns - 2)*kTolResM*max(costh, 0.2)){
      if (fVerbose > 3)
        cout << " Pattern REJECTED." << endl;
      return false;
    }

    if (fVerbose > 3)
      cout << " St. ACCEPTED." << endl;
    fResiduals.push_back(dt);
  }

  if (sqrt(sumdt2 / (ns - 3)) > (ns - 2)*kTolRes*max(costh, 0.2)) {
    if (fVerbose > 3)
      cout << "    REJECTED by tolerance on mean RMS of residuals" << endl;
    return false;
  }

  double sin2th = fU*fU + fV*fV;
  if (sin2th >= 1 && bFirstIteration) {
    if (fVerbose > 3)
      cout << "    REJECTED (firstIt && sin2th > 1)" << endl;
    return false;
  }
  if (sin2th >= 1) {
    fU /= sqrt(sin2th);
    fV /= sqrt(sin2th);
  }

  return true;
}


bool
TopDownSelector::IsGoodSpaceConfig()
{
  const CoordinateSystemPtr siteCS =
    det::Detector::GetInstance().GetSiteCoordinateSystem();
  unsigned int ns = fStations.size();
  double dist2_max_proj = pow(1300., 2.) * (ns - 2);

  if (fVerbose > 2)
    cout << "   Checking Space Configuration..." << endl;

  if (fVerbose > 3)
    cout << "    Maximum Proj. Dist squared: " << dist2_max_proj << endl;

  for (int i = ns-1; i >= 0; --i) {
    if (!fStations[i].candidate)
      continue;
    const double dx = fStations[i].detector->GetPosition().GetX(siteCS) - fXCore;
    const double dy = fStations[i].detector->GetPosition().GetY(siteCS) - fYCore;
    const double dist2 = pow(dx, 2.) + pow(dy, 2.) - pow(fU*dx + fV*dy, 2.);
    if (fVerbose > 3)
      cout << "    St: " << fStations[i].event->GetId()
           << "\tDist square: " << dist2 << " --> ";

    if (dist2 > dist2_max_proj) {
      if (fVerbose > 3)
        cout << " Pattern REJECTED." << endl;
      return false;
    }
    if (fVerbose > 3)
      cout << " St. ACCEPTED." << endl;
  }
  return true;
}


bool
TopDownSelector::IsAligned()
{
  const CoordinateSystemPtr siteCS =
    det::Detector::GetInstance().GetSiteCoordinateSystem();
  unsigned int ns = fStations.size();
  double xmean = 0, ymean = 0;
  for (unsigned int is = 0; is < ns; ++is) {
    if (!fStations[is].candidate)
      continue;
    xmean += fStations[is].detector->GetPosition().GetX(siteCS) / double(ns);
    ymean += fStations[is].detector->GetPosition().GetY(siteCS) / double(ns);
  }
  vector<double> dx(ns);
  vector<double> dy(ns);
  for (unsigned int is = 0; is < ns; ++is) {
    if (!fStations[is].candidate)
      continue;
    dx[is] = fStations[is].detector->GetPosition().GetX(siteCS) - xmean;
    dy[is] = fStations[is].detector->GetPosition().GetY(siteCS) - ymean;
  }
  double ixx = 0, ixy = 0, iyy = 0;
  for (unsigned int is = 0; is < ns; ++is) {
    if (!fStations[is].candidate)
      continue;
    ixx += dx[is]*dx[is];
    iyy += dy[is]*dy[is];
    ixy += dx[is]*dy[is];
  }
  const double zeta = 0.5*atan2(2*ixy, ixx - iyy);
  vector<double> lg(ns);
  vector<double> tr(ns);
  for (unsigned int is = 0; is < ns; ++is) {
    if (!fStations[is].candidate)
      continue;
    lg[is] = dx[is]*cos(zeta) + dy[is]*sin(zeta);
    tr[is] = dx[is]*sin(zeta) - dy[is]*cos(zeta);
  }
  double longi = 0, trans = 0;
  for (unsigned int is = 0; is < ns; ++is) {
    if (!fStations[is].candidate)
      continue;
    longi += lg[is]*lg[is];
    trans += tr[is]*tr[is];
  }
  longi = sqrt(longi);
  trans = sqrt(trans);
  if (longi < 300)
    return true;

  if (trans < 300)
    return true;

  return false;
}


double
TopDownSelector::EstimatedZenith(const string& s)
{
  if (s == "deg")
    return asin(sqrt(fU*fU + fV*fV)) / degree;
  return asin(sqrt(fU*fU + fV*fV));
}


double
TopDownSelector::EstimatedAzimuth(const string& s)
{
  if (s == "deg")
    return atan2(fV, fU) / degree;
  return atan2(fV, fU);
}


bool
TopDownSelector::Recover3Stations()
{
  const CoordinateSystemPtr siteCS =
    det::Detector::GetInstance().GetSiteCoordinateSystem();
  const sdet::SDetector& sdetector = det::Detector::GetInstance().GetSDetector();

  int index[3];
  unsigned int ns = fStations.size();
  for (unsigned int i = 0, j = 0; i < ns; ++i) {
    if (!fStations[i].candidate)
      continue;

    index[j] = i;
    ++j;
  }

  // area and maximal distance
  const double st1X = sdetector.GetStation(*fStations[index[0]].event).GetPosition().GetX(siteCS);
  const double st1Y = sdetector.GetStation(*fStations[index[0]].event).GetPosition().GetY(siteCS);
  const double st2X = sdetector.GetStation(*fStations[index[1]].event).GetPosition().GetX(siteCS);
  const double st2Y = sdetector.GetStation(*fStations[index[1]].event).GetPosition().GetY(siteCS);
  const double st3X = sdetector.GetStation(*fStations[index[2]].event).GetPosition().GetX(siteCS);
  const double st3Y = sdetector.GetStation(*fStations[index[2]].event).GetPosition().GetY(siteCS);
  const double deter =
    st1X * st2Y - st1Y * st2X +
    st2X * st3Y - st2Y * st3X +
    st3X * st1Y - st3Y * st1X;
  double dist2 = pow(st1X - st2X, 2.) + pow(st1Y - st2Y, 2.);
  dist2 = max(dist2, pow(st2X - st3X, 2.) + pow(st2Y - st3Y, 2.));
  dist2 = max(dist2, pow(st3X - st1X, 2.) + pow(st3Y - st1Y, 2.));

  // geometrical extension : cut on area of the triangle
  if (fabs(0.5*deter) > kAreaMax || fabs(0.5*deter) < kAreaMin) {
    return false;
  }
  // cut on distance between stations
  if (dist2 > pow(kDistMax, 2.)) {
    return false;
  }
  //--- computation of direction

  double dt1 = (fStations[index[1]].event->GetRecData().GetSignalStartTime() -
                fStations[index[0]].event->GetRecData().GetSignalStartTime()).GetInterval();
  double dt2 = (fStations[index[2]].event->GetRecData().GetSignalStartTime() -
                fStations[index[0]].event->GetRecData().GetSignalStartTime()).GetInterval();

  fU = kSpeedOfLight * ((st1Y - st3Y)*dt1 - (st1Y - st2Y)*dt2) / deter;
  fV = kSpeedOfLight * ((st3X - st1X)*dt1 - (st2X - st1X)*dt2) / deter;

  EstimateCore();

  fIsLastTimingIteration = false;
  if (!IsGoodTimeConfig(eAltitude)) {
    return false;
  }

  //  fIsLastTimingIteration=true; if (!TopDownIsGoodTimeConfig(2)) return false;
  double sin2th = fU*fU + fV*fV;
  if (sin2th > 1) {
    return false; // no SINTHMAX condition...
  }
  return true;
}