SdTopDownSignalSelector.cc

Go to the documentation of this file.

#include "SdTopDownSignalSelector.h"

#include <algorithm>

#include <utl/ErrorLogger.h>
#include <utl/Branch.h>
#include <utl/Accumulator.h>

#include <fwk/CentralConfig.h>
#include <fwk/LocalCoordinateSystem.h>

#include <evt/Event.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerSRecData.h>

#include <sevt/StationTriggerData.h>
#include <sevt/EventTrigger.h>
#include <sevt/SortCriteria.h>

#include <sdet/SDetector.h>

#include <TCdasUtil.h>

using namespace std;
using namespace sevt;
using namespace utl;
using namespace fwk;


static bool bFirstIteration;


namespace SdTopDownSignalSelectorUGR {

  Point SdTopDownSignalSelectorUGR::fBarycenter;
  TimeStamp SdTopDownSignalSelectorUGR::fBaryTime;


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


  VModule::ResultFlag
  SdTopDownSignalSelectorUGR::Init()
  {
    // Initialize your module here. This method
    INFO("SdTopDownSignalSelectorUGR::Init: INFO:");

    Branch topB =
      CentralConfig::GetInstance()->GetTopBranch("SdTopDownSignalSelectorUGR");

    topB.GetChild("Verbose").GetData(fVerbose);
    fVerbose = (fVerbose < 0) ? 0 : fVerbose;

    topB.GetChild("SignalThreshold").GetData(fSignalThreshold);

    topB.GetChild("Nmin").GetData(fMinSlots);
    if (fMinSlots < 1) {
      cerr << "WARNING! You specified a non valid Nmin (" << fMinSlots << "). Should be at least 1. Setting it to default (10 Slots)" << endl;
      fMinSlots = 10;
    }

    topB.GetChild("Smin").GetData(fMinSignal);
    if (fMinSignal <= 0) {
      cerr << "WARNING! You specified a non valid Smin (" << fMinSignal << "). Should be bigger than 0. Setting it to default (1 VEMs)" << endl;
      fMinSignal = 1;
    }

    topB.GetChild("N0").GetData(fMingap);
    if (fMingap < 0) {
      cerr << "WARNING! You specified a negative N0 (" << fMingap << "). Setting it to 0 slots." << endl;
      fMingap = 0;
    }

    topB.GetChild("RejectionProcedure").GetData(fRejectionProcedure);

    topB.GetChild("MinStationSignal").GetData(fMinStSignal);
    if (fMinStSignal < 0) {
      cerr << "WARNING! You specified a non valid MinStationSignal (" << fMinStSignal << "). Should be at least 0. Setting it to 0." << endl;
      fMinStSignal = 0;
    }

    topB.GetChild("MinNumberOfStationForSignalCut").GetData(fMinNumberStSignalCut);
    if (fMinNumberStSignalCut < 0) {
      cerr << "WARNING! You specified a non valid MinNumberOfStationForSignalCut (" << fMinNumberStSignalCut << "). Should be at least 0. Setting it to 0." << endl;
      fMinNumberStSignalCut = 0;
    }

    topB.GetChild("AnodeSignalFactor").GetData(fLowFactor);
    if (fLowFactor <= 0) {
      cerr << "WARNING! You specified a non valid AnodeSignalFactor (" << fLowFactor << "). Should be bigger than 0. Setting it to 1." << endl;
      fLowFactor = 1;
    }

    // Get Top Down Parameters
    // Pre Segment identification
    topB.GetChild("PreSignalThreshold").GetData(fPreSignalThreshold);

    topB.GetChild("PreNmin").GetData(fPreMinSlots);
    if (fPreMinSlots < 1) {
      cerr << "WARNING! You specified a non valid Nmin (" << fPreMinSlots << "). Should be at least 1. Setting it to default (3 Slots)" << endl;
      fPreMinSlots = 3;
    }

    topB.GetChild("PreN0").GetData(fPreMingap);
    if (fPreMingap < 0) {
      cerr << "WARNING! You specified a non valid GAP (" << fPreMingap << "). Should be at least 0. Setting it to default (3 Slots)" << endl;
      fPreMingap = 3;
    }

    topB.GetChild("PreSmin").GetData(fPreMinSignal);
    if (fPreMinSignal <= 0) {
      cerr << "WARNING! You specified a non valid PreSmin (" << fPreMinSignal << "). Should be bigger than 0. Setting it to default (0.8 VEMs)" << endl;
      fPreMinSignal = 0.8;
    }

    // Top Down Parameters
    topB.GetChild("RejectNonT4Events").GetData(fRejectNonT4Events);
    topB.GetChild("MinNumberStations").GetData(fNMinStations);
    topB.GetChild("MaxNumberStations").GetData(fNMaxStations);
    topB.GetChild("MaxRejectedSignals").GetData(fNMaxRemovedStations);
    topB.GetChild("BaryPower").GetData(kBaryPower);
    topB.GetChild("TimeResidualTolerance").GetData(kTolRes);
    topB.GetChild("MTimeResidualTolerance").GetData(kTolResM);
    topB.GetChild("IsolatedDistance1").GetData(kIsoDist1);
    topB.GetChild("IsolatedDistance2").GetData(kIsoDist2);
    topB.GetChild("Curv0").GetData(kCurv0);
    topB.GetChild("AreaMin").GetData(kAreaMin);
    topB.GetChild("AreaMax").GetData(kAreaMax);
    topB.GetChild("Ana3RatMin").GetData(kAna3RatMin);
    topB.GetChild("DistMax").GetData(kDistMax);
    topB.GetChild("IsoTime1").GetData(kIsoTime1);
    topB.GetChild("IsoTime2").GetData(kIsoTime2);

    // Get SdCalibrator Parameters

    Branch topCB =
      CentralConfig::GetInstance()->GetTopBranch("SdCalibrator");
    if (topCB) {
      topCB.GetChild("riseTimeFractions").GetData(fRiseTimeFractions);
      topCB.GetChild("fallTimeFractions").GetData(fFallTimeFractions);
      if (fRiseTimeFractions.first  < 0 || fRiseTimeFractions.first  > 1 ||
          fRiseTimeFractions.second < 0 || fRiseTimeFractions.second > 1 ||
          fFallTimeFractions.first  < 0 || fFallTimeFractions.first  > 1 ||
          fFallTimeFractions.second < 0 || fFallTimeFractions.second > 1) {

        ERROR("Rise/fall time fractions must be within [0, 1]");
        return eFailure;
      }

      if (fRiseTimeFractions.first >= fRiseTimeFractions.second ||
          fFallTimeFractions.first >= fFallTimeFractions.second ||
          fRiseTimeFractions.first >= fFallTimeFractions.second) {

        ERROR("Rise/fall time definition is not in the ascending order");
        return eFailure;
      }
    } else {
      WARNING("Rise/Fall Time Fractions can not be readed. Did you use the SdTraceCalibratorOG?\n Using default values");
      fRiseTimeFractions.first = 0.1;
      fRiseTimeFractions.second = 0.5;
      fFallTimeFractions.first = 0.5;
      fFallTimeFractions.second = 0.9;
    }

    if (fVerbose) {
      cout << " Accidental Signal Rejector parameters: (See GAP 2005-074)\n"
           "\t          Verbosity: " << fVerbose << "\n"
           "\t   Signal Threshold: " << fSignalThreshold << "\n"
           "\t   Min Slots (Nmin): " << fMinSlots << "\n"
           "\t Min. Signal (Smin): " << fMinSignal << "\n"
           "\t      Min. GAP (N0): " << fMingap << "\n"
           "\tRejection Procedure: " << fRejectionProcedure << "\n"
           "\tAnode Signal factor: " << fLowFactor << endl;

      if (fMinStSignal)
        cout << "\t   Reject St. below: " << fMinStSignal << " VEMs\n"
             "\t  Min. St. Sig. Cut: " << fMinNumberStSignalCut << endl;

      if (fRejectionProcedure == 4)
        cout << "\t PreSignalThreshold: " << fPreSignalThreshold << " VEMs\n"
             "\t           Pre Nmin: " << fPreMinSlots << "\n"
             "\t             Pre N0: " << fPreMingap << "\n"
             "\t           Pre Smin: " << fPreMinSignal << "\n"
             "\t Reject non T4 evt.: " << fRejectNonT4Events << "\n"
             "\t    Min. number St.: " << fNMinStations << "\n"
             "\t    Max. number St.: " << fNMaxStations << "\n"
             "\t  Max Rejected Sig.: " << fNMaxRemovedStations << "\n"
             "\t         Bary Power: " << kBaryPower << "\n"
             "\t Time Residual Tol.: " << kTolRes << "\n"
             "\t  M. Time Res. Tol.: " << kTolResM << "\n"
             "\t   Isolated Dist. 1: " << kIsoDist1 << "\n"
             "\t   Isolated Dist. 2: " << kIsoDist2 << "\n"
             "\t        Curvature 0: " << kCurv0 << "\n"
             "\t          Area Min.: " << kAreaMin << "\n"
             "\t          Area Max.: " << kAreaMax << "\n"
             "\t   Ana. 3 Rat. Min.: " << kAna3RatMin << "\n"
             "\t         Dist. Max.: " << kDistMax << "\n"
             "\t    Isolated Time 1: " << kIsoTime1 << "\n"
             "\t    Isolated Time 2: " << kIsoTime2 << endl;
      cout << endl;
    }

    return eSuccess;
  }


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

    if (!event.HasSEvent()) {
      INFO("Event has not SEvent");
      return eContinueLoop;
    }

    fSEvent = &event.GetSEvent();

    if (fRejectionProcedure == 4.) { // Do we really need this? The SdEventSelector should do this.
      cout << " Removing non grid stations." << endl;
      RemoveEA();
      RemoveDoublets();
      RemoveIsolatedStations();
      fAllSegments.clear();
    }

    // Sort by Increasing signal. important for cut station with low signal.
    fSEvent->SortStations(sevt::ByIncreasingSignal());

    if (fVerbose > 2) {
      cout << "   Sorted Stations by signal:" << endl;
      for (sevt::SEvent::CandidateStationIterator sIt = fSEvent->CandidateStationsBegin();
           sIt != fSEvent->CandidateStationsEnd(); ++sIt)
        cout << "    St: " << sIt->GetId() << "\t Signal: " << sIt->GetRecData().GetTotalSignal() << " VEMs" << endl;
    }

    int candidateStations = 0;

    if (fVerbose)
      cout << " Station pre-processing... " << endl;

    for (sevt::SEvent::CandidateStationIterator sIt = fSEvent->CandidateStationsBegin();
         sIt != fSEvent->CandidateStationsEnd(); ++sIt) {

      // Safety filters for accidental uses of the module before the SdCalibrator
      // It should never happens in Candidate Stations

      if (!sIt->HasTriggerData()) {
        if (fVerbose)
          cout << "  St: " << sIt->GetId() << " -> Skipped. No Trigger Data." << endl;
        sIt->SetRejected(sevt::StationConstants::eNoTrigger);
        continue;
      }
      const StationTriggerData& trig = sIt->GetTriggerData();
      if (trig.GetErrorCode()) {
        if (fVerbose)
          cout << "  St: " << sIt->GetId() << " -> Skipped. Error by code:" << trig.GetErrorCode() << endl;
        sIt->SetRejected(sevt::StationConstants::eErrorCode);
        continue;
      }
      if (trig.IsRandom()) {
        if (fVerbose)
          cout << "  St: " << sIt->GetId() << " -> Skipped. Random Trigger." << endl;
        sIt->SetRejected(sevt::StationConstants::eRandom);
        continue;
      }
      if (!sIt->HasCalibData()) {
        if (fVerbose)
          cout << "  St: " << sIt->GetId() << " -> Skipped. No Calib Data." << endl;
        sIt->SetRejected(sevt::StationConstants::eNoCalibData);
        continue;
      }
      // exclude FADCs with Patrick's data
      if (sIt->GetCalibData().GetVersion() > 32000) {
        if (fVerbose)
          cout << "  St: " << sIt->GetId() << " -> Skipped. Calib vervison > 32000." << endl;
        sIt->SetRejected(sevt::StationConstants::eBadCalib);
        continue;
      }
      if (!sIt->HasGPSData()) {
        if (fVerbose)
          cout << "  St: " << sIt->GetId() << " -> Skipped. No GPS Data." << endl;
        sIt->SetRejected(sevt::StationConstants::eNoGPSData);
        continue;
      }
      if (!sIt->HasVEMTrace()) {
        if (fVerbose)
          cout << "  St: " << sIt->GetId() << " -> Skipped. No VEM signal." << endl;
        sIt->SetRejected(sevt::StationConstants::eNoRecData);
        continue;
      }
      ++candidateStations;
    }

    if (fVerbose > 1)
      cout << "  Event with " << candidateStations << " candidate stations" << endl;

    // trivial selection
    if (fRejectionProcedure == 4 && fRejectNonT4Events && candidateStations < fNMinStations) {
      INFO("Event was not selected.");
      return eContinueLoop;
    }

    if (fVerbose)
      cout << " Station processing... " << endl;

    // Iterate over candidate stations
    for (sevt::SEvent::CandidateStationIterator sIt = fSEvent->CandidateStationsBegin();
         sIt != fSEvent->CandidateStationsEnd(); ++sIt) {

      if (fVerbose)
        cout << "\n Processing station: " << sIt->GetId() << endl;

      if (sIt->GetRecData().GetTotalSignal() < fMinStSignal && candidateStations > fMinNumberStSignalCut) {
        if (fVerbose)
          cout << "  Rejected. Signal " << sIt->GetRecData().GetTotalSignal() << " VEMs below the minimun " << fMinStSignal << " VEMs" << endl;
        sIt->SetRejected(sevt::StationConstants::eRandom);
        --candidateStations;
        continue;
      }

      if (fRejectionProcedure != 4 || candidateStations >= fNMaxStations) {
        // Search for different segments
        GeoSegmentCollection segments(FoundSegments(*sIt));

        if (segments.size() < 1) { // Delete station
          sIt->SetRejected(sevt::StationConstants::eRandom);
          if (fVerbose)
            cout << "  Station " << sIt->GetId() << " rejected. No valid segments. Flagged as random." << endl;
          --candidateStations;
          continue;
        }
        SelectMainSegment(segments);

        // Set New station values
        const int start = segments[0].fStart;
        const int stop = segments[0].fStop + 1;
        if (fVerbose)
          cout << "  Selected Segment [StartBin, StopBin): [" << start << " , " << stop << ")" << endl;
        if (fRejectionProcedure == 4) {
          fAllSegments.push_back(segments[0]);
          if (fVerbose)
            cout << "      Segment Stored for Top Down selection " << endl;
        }
        if (!sIt->HasRecData())
          sIt->MakeRecData();
        UpdateStationValues(*sIt, start, stop);
      } else {
        GeoSegmentCollection BasicSegments(FoundBasicSegments(*sIt));
        for (unsigned int seg = 0; seg < BasicSegments.size(); ++seg)
          fAllSegments.push_back(BasicSegments[seg]);
      }

    } // Station iterator


    if (fRejectionProcedure == 4.) {

      if (fVerbose)
        cout << "  Using Top Down signal selector...";
      if (fVerbose > 2) {
        cout << " ...for the segments:\n"
             "   (StId, Start, Stop, Time, Signal)" << endl;
        for (unsigned int seg = 0; seg < fAllSegments.size(); ++seg)
          cout << "    (" << fAllSegments[seg].fStEvt->GetId() << " , "
               << fAllSegments[seg].fStart << " , "
               << fAllSegments[seg].fStop << " , "
               << fAllSegments[seg].fStartTime.GetGPSNanoSecond() << " , "
               << fAllSegments[seg].fSignal << endl;
      }
      if (fVerbose)
        cout << endl;

      // trivial selection
      if (fRejectNonT4Events && candidateStations < fNMinStations) {
        INFO("Event was not selected.");
        return eContinueLoop;
      }

      const bool isFD = fSEvent->GetTrigger().IsFD();
      const bool isSdT4 = Select();  // Top Down procedure here!!.

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

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

      if (isSdT4) {
        EstimateCore();
        const Vector axis(-fU, -fV, cos(EstimatedZenith("rad")), LocalCoordinateSystem::Create(fBarycenter));
        shSRec.SetAxis(axis);

        // shSRec.SetAngleErrors(axis.GetCoordinates(fgLocalCS),
        //                        Vector::Triple(fit->sigmaU2, fit->sigmaUV, fit->sigmaV2));

        shSRec.SetCurvature(0, 0);
        shSRec.SetBarycenter(fBarycenter);
        shSRec.SetCorePosition(fBarycenter);
        shSRec.SetCoreTime(fBaryTime + TimeInterval(fT0), 0); //kSpeedOfLight),
        // TimeInterval(sqrt(fit->sigmact02)));

        // double timeMean3d;
        // double timeSpread3d;
        // double chi2;
        // CalculateTimeResidual3D(axis, fT0/kSpeedOfLight, timeMean3d, timeSpread3d, chi2);

        // shSRec.SetTimeResidualMean(timeMean3d);
        // shSRec.SetTimeResidualSpread(timeSpread3d);
        // const int angleNDOF = nStations - 3;
        // shSRec.SetAngleChi2(chi2, angleNDOF);

        //  shSRec.SetLDFRecStage(fit->stage);
      }

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

    return eSuccess;
  }


  VModule::ResultFlag
  SdTopDownSignalSelectorUGR::Finish()
  {
    // Put any termination or cleanup code here.
    // This method is called once at the end of the run.

    INFO("SdTopDownSignalSelectorUGR::Finish");

    return eSuccess;
  }


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

  bool
  SdTopDownSignalSelectorUGR::Select()
  {
    sort(fAllSegments.begin(), fAllSegments.end(), ByIncreasingSignal());

    const int ns = fAllSegments.size();

    int k = (ns - fNMinStations) < fNMaxRemovedStations ? ns - fNMinStations : fNMaxRemovedStations;

    if (fVerbose > 4)
      cout << " Debuged info: Number of segments: " << ns << "\n"
           "               Minimun of stations: " << fNMinStations << "\n"
           "               MaxRemoved Stations: " << fNMaxRemovedStations << "\n"
           "               Index k: " << k << endl;

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

    // "top down" selection
    bFirstIteration = true;
    for (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 (int ik = 0; ik < i; ++ik)
          cout << "   Select: fIndexes[" << ik << "] = " << fIndexes[ik] << endl;
      }
      if (TryRemovingStations(i)) {
        UpdateStations();

        if (fVerbose > 2) {
          cout << " Selected segments:\n"
               "   (StId, Start, Stop, Time, Signal)" << endl;
          for (unsigned int seg = 0; seg < fAllSegments.size(); ++seg) {
            if (!fAllSegments[seg].fCandidate)
              continue;
            cout << "    (" << fAllSegments[seg].fStEvt->GetId() << " , "
                 << fAllSegments[seg].fStart << " , "
                 << fAllSegments[seg].fStop << " , "
                 << fAllSegments[seg].fStartTime.GetGPSNanoSecond() << " , "
                 << fAllSegments[seg].fSignal << endl;
          }
        }

        if (fVerbose > 0) {
          cout << " Event accepted with " << i << " rejected. Theta ~ "
               << EstimatedZenith("deg") << " deg." << endl;
        }
        return true;
      }
    }

    bFirstIteration = false;
    for (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 (int ik = 0; ik < i; ++ik)
          cout << "   Select: fIndexes[" << ik << "] = " << fIndexes[ik] << endl;
      }
      if (TryRemovingStations(i)) {
        UpdateStations();

        if (fVerbose > 2) {
          cout << " Selected segments:\n"
               "   (StId, Start, Stop, Time, Signal)" << endl;
          for (unsigned int seg = 0; seg < fAllSegments.size(); ++seg) {
            if (!fAllSegments[seg].fCandidate)
              continue;
            cout << "    (" << fAllSegments[seg].fStEvt->GetId() << " , "
                 << fAllSegments[seg].fStart << " , "
                 << fAllSegments[seg].fStop << " , "
                 << fAllSegments[seg].fStartTime.GetGPSNanoSecond() << " , "
                 << fAllSegments[seg].fSignal << endl;
          }
        }

        if (fVerbose > 0) {
          cout << " Event accepted with " << i << " rejected. Theta ~ "
               << EstimatedZenith("deg") << " deg." << endl;
        }
        return true;
      }
    }

    return false;
  }


  void
  SdTopDownSignalSelectorUGR::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 IsInGrid
  // by the corresponding one
  void
  SdTopDownSignalSelectorUGR::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
  SdTopDownSignalSelectorUGR::RemoveIsolatedStations()
  {
    const sdet::SDetector& sdetector = det::Detector::GetInstance().GetSDetector();

    for (sevt::SEvent::CandidateStationIterator sIt_1 = fSEvent->CandidateStationsBegin();
         sIt_1 != fSEvent->CandidateStationsEnd(); ++sIt_1) {
      const utl::Point fStPoint1 = sdetector.GetStation(*sIt_1).GetPosition();
      int n1 = 0;
      int n2 = 0;
      const 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);
    }
  }


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

  bool
  SdTopDownSignalSelectorUGR::PatternRecovering()
  {
    unsigned int ns = 0;
    for (unsigned int is = 0; is < fAllSegments.size(); ++is)
      if (fAllSegments[is].fCandidate)
        ++ns;

    if (ns == 3) {
      if (fNMinStations == 3)
        return Recover3Stations(); // Note that 3 segments could not imply 3 stations
      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
  SdTopDownSignalSelectorUGR::TryRemovingStations(const int n)
  {
    if (n < 2) {
      for (unsigned int is = 0; is < fIndexes[0]; ++is) {
        fAllSegments[is].fCandidate = false;
        if (PatternRecovering())
          return true;
        fAllSegments[is].fCandidate = true;
      }
    } else {
      for (unsigned int is = n - 1; is < fIndexes[n - 1]; ++is) {
        fAllSegments[is].fCandidate = false;
        fIndexes[n - 2] = is;
        if (TryRemovingStations(n - 1))
          return true;
        fAllSegments[is].fCandidate = true;
      }
    }
    return false;
  }


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

    Vector barySum(0, 0, 0, siteCS);
    TimeInterval timeSum = 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;

    const unsigned int size = fAllSegments.size();
    double sumw = 0;
    for (unsigned int is = 0; is < size; ++is) {
      if (!fAllSegments[is].fCandidate)
        continue;
      const sevt::Station& station = *fAllSegments[is].fStEvt;

      const double weight = exp(kBaryPower * log(fAllSegments[is].fSignal));
      barySum += weight * (sdetector.GetStation(station).GetPosition() - siteOrigin);
      timeSum += weight *
                 (fAllSegments[is].fStartTime - eventTime).GetInterval();

      sumw += weight;
    }
    if (sumw == 0) {
      if (fVerbose > 2)
        cout << "    ...FAILED." << endl;
      return false;
    }
    barySum /= sumw;
    fBarycenter = siteOrigin + barySum;
    fBaryTime = eventTime + timeSum;

    timeSum /= sumw;

    fTCore = timeSum.GetInterval();
    fXCore = barySum.GetX(siteCS);
    fYCore = barySum.GetY(siteCS);
    fZCore = barySum.GetZ(siteCS);

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

    return true;
  }


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


  double
  SdTopDownSignalSelectorUGR::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 = fAllSegments[is].fStDet->GetPosition().GetX(siteCS);
    const double dy = fAllSegments[is].fStDet->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 * (fAllSegments[is].fStDet->GetPosition().GetZ(siteCS) - fZCore) - curv * d2 / 2.) / kSpeedOfLight;
  }


  bool
  SdTopDownSignalSelectorUGR::IsGoodTimeConfig(const TimeCorrectionType correctionType)
  {
    const CoordinateSystemPtr siteCS =
      det::Detector::GetInstance().GetSiteCoordinateSystem();
    unsigned int ns = 0;

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

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

      const double dx = fAllSegments[is].fStDet->GetPosition().GetX(siteCS) - fXCore;
      const double dy = fAllSegments[is].fStDet->GetPosition().GetY(siteCS) - fYCore;

      double dt;
      if (correctionType == SdTopDownSignalSelectorUGR::eAltitude)
        dt = CorrectTimeAltitude(is);
      else if (correctionType == SdTopDownSignalSelectorUGR::eAltitudeCurved)
        dt = CorrectTimeAltitudeCurv(is);
      else if (correctionType == SdTopDownSignalSelectorUGR::eNone)
        dt = 0;
      else {
        dt = 0;
        INFO("Timing Correction Type not recognized. No correction done.");
      }
      dt += (fBaryTime - fAllSegments[is].fStartTime).GetInterval();
      const double 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;
    }

    const double rxx = sumw * sumy2 - sumy * sumy;
    const double ryy = sumw * sumx2 - sumx * sumx;
    const double rxy = sumx * sumy - sumw * sumxy;
    const double rx = sumxy * sumy - sumx * sumy2;
    const double ry = sumxy * sumx - sumy * sumx2;
    const double rr = sumx2 * sumy2 - sumxy * sumxy;
    const double 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;
    const 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 (unsigned int is = 0; is < size; ++is) {
      if (!fAllSegments[is].fCandidate || !fAllSegments[is].fStEvt->HasRecData())
        continue; // trying patterns...
      const double dx = fAllSegments[is].fStDet->GetPosition().GetX(siteCS) - fXCore;
      const double dy = fAllSegments[is].fStDet->GetPosition().GetY(siteCS) - fYCore;
      const double dt =
        (fBaryTime - fAllSegments[is].fStartTime).GetInterval() - (fT0 - fU * dx / kSpeedOfLight - fV * dy / kSpeedOfLight);
      sumdt2 += dt * dt;

      if (fVerbose > 3)
        cout << "    St: " << fAllSegments[is].fStEvt->GetId()
             << "\tResidual: " << dt << " fabs: " << fabs(dt) << "\tLimit: " << (ns - 2)*kTolResM* max(costh, 0.2) << " --> ";

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

      if (fVerbose > 3)
        cout << " Segment ACCEPTED." << endl;
      fAllSegments[is].fTimeResidual = dt;
      fResiduals.push_back(dt);
    }

    if (sqrt(sumdt2 / (ns - 3)) > (ns - 2)*kTolRes * max(costh, 0.2)) {
      if (fVerbose > 3)
        cout << "    Mean RMS: " << sqrt(sumdt2 / (ns - 3)) << "\tMaxValue: " << (ns - 2)*kTolRes* max(costh, 0.2) << " --> REJECTED " << endl;
      return false;
    }

    const 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
  SdTopDownSignalSelectorUGR::IsGoodSpaceConfig()
  {
    const CoordinateSystemPtr siteCS =
      det::Detector::GetInstance().GetSiteCoordinateSystem();
    const unsigned int ns = fAllSegments.size();
    const 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 (!fAllSegments[i].fCandidate)
        continue;
      const double dx = fAllSegments[i].fStDet->GetPosition().GetX(siteCS) - fXCore;
      const double dy = fAllSegments[i].fStDet->GetPosition().GetY(siteCS) - fYCore;
      const double dist2 = pow(dx, 2) + pow(dy, 2) - pow(fU * dx + fV * dy, 2);
      if (fVerbose > 3)
        cout << "    St: " << fAllSegments[i].fStEvt->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
  SdTopDownSignalSelectorUGR::IsAligned()
  {
    const CoordinateSystemPtr siteCS =
      det::Detector::GetInstance().GetSiteCoordinateSystem();
    const unsigned int ns = fAllSegments.size();
    double xmean = 0;
    double ymean = 0;
    for (unsigned int is = 0; is < ns; ++is) {
      if (!fAllSegments[is].fCandidate)
        continue;
      xmean += fAllSegments[is].fStDet->GetPosition().GetX(siteCS);
      ymean += fAllSegments[is].fStDet->GetPosition().GetY(siteCS);
    }
    xmean /= ns;
    ymean /= ns;
    vector<double> dx(ns);
    vector<double> dy(ns);
    for (unsigned int is = 0; is < ns; ++is) {
      if (!fAllSegments[is].fCandidate)
        continue;
      dx[is] = fAllSegments[is].fStDet->GetPosition().GetX(siteCS) - xmean;
      dy[is] = fAllSegments[is].fStDet->GetPosition().GetY(siteCS) - ymean;
    }
    double ixx = 0;
    double ixy = 0;
    double iyy = 0;
    for (unsigned int is = 0; is < ns; ++is) {
      if (!fAllSegments[is].fCandidate)
        continue;
      ixx += dx[is] * dx[is];
      iyy += dy[is] * dy[is];
      ixy += dx[is] * dy[is];
    }
    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 (!fAllSegments[is].fCandidate)
        continue;
      lg[is] = dx[is] * cos(zeta) + dy[is] * sin(zeta);
      tr[is] = dx[is] * sin(zeta) - dy[is] * cos(zeta);
    }
    double longi = 0;
    double trans = 0;
    for (unsigned int is = 0; is < ns; ++is) {
      if (!fAllSegments[is].fCandidate)
        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
  SdTopDownSignalSelectorUGR::EstimatedZenith(const string& s)
  {
    if (s == "deg")
      return asin(sqrt(fU * fU + fV * fV)) / degree;
    return asin(sqrt(fU * fU + fV * fV));
  }


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


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

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

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

    // area and maximal distance
    const double st1X = sdetector.GetStation(*fAllSegments[index[0]].fStEvt).GetPosition().GetX(siteCS);
    const double st1Y = sdetector.GetStation(*fAllSegments[index[0]].fStEvt).GetPosition().GetY(siteCS);
    const double st2X = sdetector.GetStation(*fAllSegments[index[1]].fStEvt).GetPosition().GetX(siteCS);
    const double st2Y = sdetector.GetStation(*fAllSegments[index[1]].fStEvt).GetPosition().GetY(siteCS);
    const double st3X = sdetector.GetStation(*fAllSegments[index[2]].fStEvt).GetPosition().GetX(siteCS);
    const double st3Y = sdetector.GetStation(*fAllSegments[index[2]].fStEvt).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(deter / 2) > kAreaMax || fabs(deter / 2) < kAreaMin)
      return false;
    // cut on distance between stations
    if (dist2 > pow(kDistMax, 2))
      return false;
    //--- computation of direction

    const double dt1 = (fAllSegments[index[1]].fStartTime -
                        fAllSegments[index[0]].fStartTime).GetInterval();
    const double dt2 = (fAllSegments[index[2]].fStartTime -
                        fAllSegments[index[0]].fStartTime).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;
    const double sin2th = fU * fU + fV * fV;
    if (sin2th > 1) {
      return false; // no SINTHMAX condition...
    }
    return true;
  }


  // Segment Identification

  GeoSegmentCollection
  SdTopDownSignalSelectorUGR::FoundSegments(sevt::Station& st)
  {
    // For station detector
    const sdet::Station& dStation =
      det::Detector::GetInstance().GetSDetector().GetStation(st);

    // For start time of the signal
    const double fadcBinSize = dStation.GetFADCBinSize();
    const StationGPSData& gpsData = st.GetGPSData();
    const TimeStamp gpsTime(gpsData.GetSecond(), gpsData.GetCorrectedNanosecond());
    const TimeStamp traceTime =
      gpsTime - TimeInterval(dStation.GetFADCTraceLength() * fadcBinSize);

    GeoSegmentCollection segments;

    const utl::TraceD& vemChargeTrace = GetVEMChargeTrace(st);
    const utl::TraceD& vemPeakTrace = st.GetVEMTrace();

    //sevt::Station::Signal subsegment(0,0,0,0);
    GeoSegment subsegment(0, 0, 0, 0., 0., 0., traceTime, &st, &dStation, true);

    double lowFactor = 1;
    if (st.IsHighGainSaturation()) {
      lowFactor *= fLowFactor;
      if (fVerbose)
        cout << "  Station with High gain trace saturated. Using factor " << lowFactor << " for signal threshold." << endl;
    }

    if (fVerbose > 2)
      cout << "\n\t\t Scaning for segments..." << endl;

    // Uncomment next lines scan only betwen start signal bin and stop signal bin.
    // sevt::StationRecData & stRec = st.GetRecData();
    // for(int bin=stRec.GetSignalStartSlot(); bin<=stRec.GetSignalEndSlot(); ++bin ){

    // Search for elemental subsegments
    for (utl::TraceD::SizeType bin = vemChargeTrace.GetStart(); bin < vemChargeTrace.GetStop(); ++bin) {
      if (fVerbose > 4)
        cout << "\t\t\t Bin: " << bin
             << "\t Value: " << vemPeakTrace[bin]
             << " VEM Peak\t Value: " << vemChargeTrace[bin] << " VEM Charge" << endl;
      if (vemChargeTrace[bin] < (fSignalThreshold * lowFactor) && subsegment.fBinsOverThresh == 0)
        continue;
      else if (vemChargeTrace[bin] < (fSignalThreshold * lowFactor) && subsegment.fBinsOverThresh != 0) {
        if (fVerbose > 2)
          cout << "\t\t Adding new subsegment (start, stop, charge, bins): ("
               << subsegment.fStart << " , "
               << subsegment.fStop  << " , "
               << subsegment.fCharge  << " , "
               << subsegment.fAoP  << " , "
               << subsegment.fBinsOverThresh  << ")" << endl;

        //      if((subsegment.fBinsOverThresh>=fMinSlots) && (subsegment.fCharge>=fMinSignal))
        segments.push_back(subsegment); // Store one segment
        // Init subsegment
        subsegment.fBinsOverThresh = 0;
        subsegment.fCharge = 0;
        subsegment.fAoP = 0;
        subsegment.fSignal = 0;
        subsegment.fStart = 0;
        subsegment.fStop = 0;
        subsegment.fStartTime = 0;
      } else if (vemChargeTrace[bin] >= (fSignalThreshold * lowFactor) && subsegment.fBinsOverThresh == 0) {
        subsegment.fStart = bin;
        subsegment.fStop = bin;
        subsegment.fCharge += vemChargeTrace[bin];
        subsegment.fSignal += vemPeakTrace[bin];
        ++subsegment.fBinsOverThresh;
        subsegment.fStartTime = traceTime + TimeInterval((bin - 0.5) * fadcBinSize);
        if (fVerbose >= 4)
          cout << "\t\t\t New segment starting: "
               << subsegment.fStart << " , "
               << subsegment.fStop << " , "
               << subsegment.fCharge << " , "
               << subsegment.fBinsOverThresh << endl;
      } else if (vemChargeTrace[bin] >= (fSignalThreshold * lowFactor) && subsegment.fBinsOverThresh != 0) {
        subsegment.fStop = bin;
        subsegment.fCharge += vemChargeTrace[bin];
        ++subsegment.fBinsOverThresh;
        subsegment.fSignal += vemPeakTrace[bin];
        if (fVerbose >= 4)
          cout << "\t\t\t Updating segment: "
               <<  subsegment.fStart << " , "
               <<  subsegment.fStop << " , "
               <<  subsegment.fCharge << " , "
               <<  subsegment.fBinsOverThresh << endl;
      }
    }

    if (segments.size() <= 1)
      return segments; // Only one segment or no segments

    if (fVerbose > 2)
      cout << "\n\t\t Merging subsegments..." << endl;
    if (fVerbose > 3)
      cout << "\t\t\t(Start, Stop, Charge, BinsOverThreshold)" << endl;

    // Merge subsegments into segments
    for (unsigned int seg = 0; seg < segments.size() - 1; ++seg) {
      //if (seg == segments.size()-2) break;
      if ((segments[seg].fStop + fMingap + segments[seg].fBinsOverThresh) > segments[seg + 1].fStart) { // Merge next segment
        if (fVerbose > 3)
          cout << "\t\t\t Merging: ("
               <<  segments[seg].fStart << " , "
               <<  segments[seg].fStop << " , "
               <<  segments[seg].fCharge << " , "
               <<  segments[seg].fBinsOverThresh << ") + ("
               <<  segments[seg + 1].fStart << " , "
               <<  segments[seg + 1].fStop << " , "
               <<  segments[seg + 1].fCharge << " , "
               <<  segments[seg + 1].fBinsOverThresh << ") = (";

        segments[seg].fStop = segments[seg + 1].fStop;
        segments[seg].fBinsOverThresh += segments[seg + 1].fBinsOverThresh;
        segments[seg].fCharge += segments[seg + 1].fCharge;
        segments[seg].fSignal += segments[seg + 1].fSignal;

        if (fVerbose > 3)
          cout <<  segments[seg].fStart << " , "
               <<  segments[seg].fStop << " , "
               <<  segments[seg].fCharge << " , "
               <<  segments[seg].fBinsOverThresh << ")" << endl;

        segments.erase(segments.begin() + seg + 1);
        --seg;
      }
    }

    if (fVerbose > 2)
      cout << "\n\t\t Deleting segments below threshold... " << endl;
    if (fVerbose > 3)
      cout << "\t\t\t(Start, Stop, Charge, BinsOverThreshold)" << endl;
    // Check for acceptable segments
    for (unsigned int seg = 0; seg < segments.size(); ++seg) {
      if (segments[seg].fBinsOverThresh >= fMinSlots && segments[seg].fCharge > fMinSignal) { // Good segment
        // Compute AoP
        double peak = 0;
        for (int bin = segments[seg].fStart; bin <= segments[seg].fStop; ++bin)
          peak = vemPeakTrace[bin] > peak ? vemPeakTrace[bin] : peak;
        segments[seg].fAoP = segments[seg].fCharge / peak;
        continue;
      }

      if (fVerbose > 3)
        cout << "\t\t\t Deleting segment: ("
             <<  segments[seg].fStart << " , "
             <<  segments[seg].fStop << " , "
             <<  segments[seg].fCharge << " , "
             <<  segments[seg].fBinsOverThresh << ")" << endl;

      segments.erase(segments.begin() + seg);
      --seg;
    }

    if (fVerbose > 1) {
      cout << "\n\t Selected segments: (Start, Stop, VEMCharge, AoP, VEMPeakSignal, BinsOverTh, SignalStartTime)" << endl;
      for (unsigned int seg = 0; seg < segments.size(); ++seg) {
        cout << "\t\t ("
             <<  segments[seg].fStart << " , "
             <<  segments[seg].fStop << " , "
             <<  segments[seg].fCharge << " , "
             <<  segments[seg].fAoP << " , "
             <<  segments[seg].fSignal << " , "
             <<  segments[seg].fBinsOverThresh << " , "
             <<  segments[seg].fStartTime.GetGPSNanoSecond() << ")" << endl;
      }
    }

    return segments;
  }


  bool
  ByIncreasingSignal::operator()(const GeoSegment& i, const GeoSegment& j)
    const
  {
    return i.fCharge < j.fCharge;
  }


  bool
  ByDecreasingSignal::operator()(const GeoSegment& i, const GeoSegment& j)
    const
  {
    cout << " Segment i: " << i.fCharge << " " << i.fCandidate << endl
         << " Segment j: " << j.fCharge << " " << j.fCandidate << endl
         << endl;
    if (!j.fCandidate)
      return true;  //
    else if (!i.fCandidate)
      return false;
    else
      return i.fCharge > j.fCharge;
  }


  bool
  ByDecreasingSignalQ::operator()(const GeoSegment& i, const GeoSegment& j)
    const
  {
    return i.fCharge * i.fBinsOverThresh > j.fCharge * j.fBinsOverThresh;
  }


  bool
  ByDecreasingAoP::operator()(const GeoSegment& i, const GeoSegment& j)
    const
  {
    return i.fAoP > j.fAoP;
  }


  bool
  ByDecreasingAoPQ::operator()(const GeoSegment& i, const GeoSegment& j)
    const
  {
    return i.fAoP * i.fBinsOverThresh > j.fAoP * j.fBinsOverThresh;
  }


  bool
  SdTopDownSignalSelectorUGR::SelectMainSegment(GeoSegmentCollection& segments)
  {
    // For the moment we select the segment with biggest signal.
    // sort(segments.begin(), segments.end(), ByDecreasingSignal());
    if (fVerbose > 1)
      cout << "\n Sorted segments by ";

    switch (fRejectionProcedure) {
    case 0:
      cout << "Signal";
      sort(segments.begin(), segments.end(), ByDecreasingSignal());
      break;
    case 1:
      cout << "Signal quality";
      sort(segments.begin(), segments.end(), ByDecreasingSignalQ());
      break;
    case 2:
      cout << "AoP";
      sort(segments.begin(), segments.end(), ByDecreasingAoP());
      break;
    case 3:
      cout << "AoP quality";
      sort(segments.begin(), segments.end(), ByDecreasingAoPQ());
      break;
    case 4:
      cout << "Signal quality (default for Top Down)";
      sort(segments.begin(), segments.end(), ByDecreasingSignalQ());
      break;
    default :
      cout << " Unrecognized Rejection Procedure. Setting to default: Signal Quality";
      fRejectionProcedure = 1;
      sort(segments.begin(), segments.end(), ByDecreasingSignalQ());
      break;
    }

    if (fVerbose > 1) {
      cout << ": " << endl;
      for (unsigned int seg = 0; seg < segments.size(); ++seg) {
        cout << "\t ("
             << segments[seg].fStart << " , "
             << segments[seg].fStop << " , "
             << segments[seg].fCharge << " , "
             << segments[seg].fAoP << " , "
             << segments[seg].fBinsOverThresh << ")" << endl;
      }
    }
    return true;
  }


  void
  SdTopDownSignalSelectorUGR::UpdateStations()
  {
    // For various segments in the same station, reject the ones with bigger Time Residual
    if (fVerbose > 2) {
      cout << " Will check duplicity of stations for the candidates out of the following segments\n"
           "\t (StId, Start, Stop, Charge, AoP, Cand, BinsOverThreshold, TimeResidual)" << endl;
      for (unsigned int i = 0; i < fAllSegments.size(); ++i) {
        cout << "\t ("
             << fAllSegments[i].fStEvt->GetId() << " , "
             << fAllSegments[i].fStart << " , "
             << fAllSegments[i].fStop << " , "
             << fAllSegments[i].fCharge << " , "
             << fAllSegments[i].fAoP << " , "
             << fAllSegments[i].fCandidate << " , "
             << fAllSegments[i].fBinsOverThresh << " , "
             << fAllSegments[i].fTimeResidual << ")" << endl;
      }
    }

    for (unsigned int i = 0; i < fAllSegments.size(); ++i) {
      if (!fAllSegments[i].fCandidate)
        continue;
      for (unsigned int j = i + 1; j < fAllSegments.size(); ++j) {
        if (!fAllSegments[j].fCandidate)
          continue;
        if (fAllSegments[i].fStEvt->GetId() == fAllSegments[j].fStEvt->GetId()) {
          if (abs(fAllSegments[i].fTimeResidual) < abs(fAllSegments[j].fTimeResidual)) {
            if (fVerbose > 2)
              cout << "   Duplicate segment with bigger residual time rejected. Station " << fAllSegments[j].fStEvt->GetId() << endl;
            fAllSegments[j].fCandidate = false;
          } else {
            if (fVerbose > 2)
              cout << "   Duplicate segment with bigger residual time rejected. Station " << fAllSegments[i].fStEvt->GetId() << endl;
            fAllSegments[i].fCandidate = false;
          }
        }
      }
    }

    // Update Stations
    if (fVerbose > 2) {
      cout << " Duplicate candidate segments checked and cleared (only one candidate per station)\n"
           "\t (StId, Start, Stop, Charge, AoP, Cand, BinsOverThreshold, TimeResidual)" << endl;
      for (unsigned int i = 0; i < fAllSegments.size(); ++i) {
        cout << "\t ("
             << fAllSegments[i].fStEvt->GetId() << " , "
             << fAllSegments[i].fStart << " , "
             << fAllSegments[i].fStop << " , "
             << fAllSegments[i].fCharge << " , "
             << fAllSegments[i].fAoP << " , "
             << fAllSegments[i].fCandidate << " , "
             << fAllSegments[i].fBinsOverThresh << " , "
             << fAllSegments[i].fTimeResidual << ")" << endl;
      }
    }

    // First we flag each rejected segment and delete it from the list.
    GeoSegmentCollection rejectedSegments;
    for (unsigned int seg = 0; seg < fAllSegments.size(); ++seg) {
      if (fAllSegments[seg].fCandidate)
        continue;
      UpdateStationValues(*fAllSegments[seg].fStEvt, fAllSegments[seg].fStart, fAllSegments[seg].fStop);
      //~ fAllSegments[seg].fStEvt->SetRejected(sevt::StationConstants::eOutOfTime);
      rejectedSegments.push_back(fAllSegments[seg]);
      fAllSegments.erase(fAllSegments.begin() + seg);
      --seg;
    }
    // Finally, if the rejected candidate's station doesn't have a candidate segment at all, we set the station to rejected.
    for (unsigned int r = 0; r < rejectedSegments.size(); ++r) {
      bool rejectCurrSt = true;
      for (unsigned int seg = 0; seg < fAllSegments.size() && rejectCurrSt; ++seg) {
        if (rejectedSegments[r].fStEvt == fAllSegments[seg].fStEvt)
          rejectCurrSt = false;
      }
      if (rejectCurrSt)
        rejectedSegments[r].fStEvt->SetRejected(sevt::StationConstants::eOutOfTime);
    }

    if (fVerbose > 2) {
      cout << " Duplicate segments deleted. Final segments to update\n"
           "\t (StId, Start, Stop, Charge, AoP, Cand, BinsOverThreshold, TimeResidual)" << endl;
      for (unsigned int i = 0; i < fAllSegments.size(); ++i) {
        cout << "\t ("
             << fAllSegments[i].fStEvt->GetId() << " , "
             << fAllSegments[i].fStart << " , "
             << fAllSegments[i].fStop << " , "
             << fAllSegments[i].fCharge << " , "
             << fAllSegments[i].fAoP << " , "
             << fAllSegments[i].fCandidate << " , "
             << fAllSegments[i].fBinsOverThresh << " , "
             << fAllSegments[i].fTimeResidual << ")" << endl;
      }
    }

    // Now we update the segments with the final criteria for a segment definition and update the station values.
    for (unsigned int seg = 0; seg < fAllSegments.size(); ++seg) {
      if (!fAllSegments[seg].fCandidate)
        continue;
      UpdateSegmentValues(seg);
      UpdateStationValues(*fAllSegments[seg].fStEvt, fAllSegments[seg].fStart, fAllSegments[seg].fStop);
      fAllSegments[seg].fStEvt->SetCandidate();
    }

    if (fVerbose > 2) {
      cout << " updated segments:\n"
           "\t (StId, Start, Stop, Charge, AoP, Cand, BinsOverThreshold, TimeResidual)" << endl;
      for (unsigned int i = 0; i < fAllSegments.size(); ++i) {
        cout << "\t ("
             <<  fAllSegments[i].fStEvt->GetId() << " , "
             <<  fAllSegments[i].fStart << " , "
             <<  fAllSegments[i].fStop << " , "
             <<  fAllSegments[i].fCharge << " , "
             <<  fAllSegments[i].fAoP << " , "
             <<  fAllSegments[i].fCandidate << " , "
             <<  fAllSegments[i].fBinsOverThresh << " , "
             <<  fAllSegments[i].fTimeResidual << ")" << endl;
      }
    }
  }


  bool
  SdTopDownSignalSelectorUGR::UpdateStationValues(sevt::Station& st, const int start, const int stop)
  {
    sevt::StationRecData& stRec = st.GetRecData();
    stRec.SetSignalStartSlot(start);
    stRec.SetSignalEndSlot(stop);

    {
      const TraceD& vemTrace = st.GetVEMTrace();
      double peak = 0;
      for (int i = start; i < stop; ++i) {
        const double val = vemTrace[i];
        if (val > peak)
          peak = val;
      }
      stRec.SetPeakAmplitude(peak);
    }

    const sdet::Station& dStation =
      det::Detector::GetInstance().GetSDetector().GetStation(st);

    // bool highGainSaturation = false;
    // bool lowGainSaturation = false;
    int nPMTs = 0;
    double totalCharge = 0;
    utl::Accumulator::SampleStandardDeviationN shapeStat;
    utl::Accumulator::SampleStandardDeviationN riseStat;
    utl::Accumulator::SampleStandardDeviationN fallStat;
    utl::Accumulator::SampleStandardDeviationN t50Stat;
    for (Station::PMTIterator pmtIt = st.PMTsBegin();
         pmtIt != st.PMTsEnd(); ++pmtIt) {
      if (pmtIt->HasCalibData() && pmtIt->GetCalibData().IsTubeOk()) {
        PMTRecData& pmtRec = pmtIt->GetRecData();
        // if (pmtRec.GetFADCSaturatedBins(PMTConstants::eHighGain)) // set but not used
        //   highGainSaturation = true;
        // if (pmtRec.GetFADCSaturatedBins(PMTConstants::eLowGain))
        // set but not used
        //   lowGainSaturation = true;
        const TraceD& vemTrace = pmtRec.GetVEMTrace();
        double charge = 0;
        for (int i = start; i < stop; ++i)
          charge += vemTrace[i];
        ComputeShapeRiseFallPeak(pmtRec, dStation.GetFADCBinSize(), start, start, stop, charge);
        charge *= pmtRec.GetVEMPeak() / pmtRec.GetVEMCharge();
        pmtRec.SetTotalCharge(charge);
        totalCharge += charge;
        shapeStat(pmtRec.GetShapeParameter());
        riseStat(pmtRec.GetRiseTime());
        fallStat(pmtRec.GetFallTime());
        t50Stat(pmtRec.GetT50());
        const double peak = pmtRec.GetPeakAmplitude();
        if (peak)
          pmtRec.SetAreaOverPeak(charge / peak);
        ++nPMTs;
      }
    }
    totalCharge /= nPMTs;

    // this is done on the pmt vem traces due to individual pmt vem peak/charge values
    stRec.SetTotalSignal(totalCharge);
    // if (totalCharge <= 0)
    //   return false;

    const StationGPSData& gpsData = st.GetGPSData();

    // timing of the trace END
    const TimeStamp gpsTime(gpsData.GetSecond(), gpsData.GetCorrectedNanosecond());

    const double fadcBinSize = dStation.GetFADCBinSize();

    // timing of the trace BEGINNING
    const TimeStamp traceTime = gpsTime -
                                TimeInterval(dStation.GetFADCTraceLength() * fadcBinSize);
    st.SetTraceStartTime(traceTime);

    // timing of the SIGNAL START
    const TimeStamp signalTime = traceTime +
                                 TimeInterval((start - 0.5) * fadcBinSize);
    stRec.SetSignalStartTime(signalTime);

    return true;
  }


  // Copy paste from SdCalibratorOG
  // set rise/fall time for PMT trace and return shape parameter
  void
  SdTopDownSignalSelectorUGR::ComputeShapeRiseFallPeak(PMTRecData& pmtRec,
      const double binSize,
      const unsigned int startBin,
      const unsigned int startIntegration,
      const unsigned int endIntegration,
      const double traceIntegral)
  const
  {
    // sorry for not using the nice TraceAlgorithms, but all this values can
    // be calculated in single trace pass

    if (traceIntegral <= 0)
      return;

    const double riseStartSentry = fRiseTimeFractions.first * traceIntegral;
    const double riseEndSentry = fRiseTimeFractions.second * traceIntegral;
    const double fallStartSentry = fFallTimeFractions.first * traceIntegral;
    const double fallEndSentry = fFallTimeFractions.second * traceIntegral;
    const double binTiming = binSize;
    const unsigned int shapeSentry =
      startIntegration + (unsigned int)(600 * nanosecond / binTiming);
    const double t50Sentry = 0.5 * traceIntegral;
    double riseStartBin = 0;
    double riseEndBin = 0;
    double fallStartBin = 0;
    double fallEndBin = 0;
    double t50Bin = 0;
    double sumEarly = 0;

    double peakAmplitude = 0;
    double runningSum = 0;
    double oldSum = 0;

    const TraceD& trace = pmtRec.GetVEMTrace();

    for (unsigned int i = startIntegration; i < endIntegration; ++i) {

      const double binValue = trace[i];
      runningSum += binValue;

      if (peakAmplitude < binValue)
        peakAmplitude = binValue;

      if (!sumEarly && i >= shapeSentry)
        sumEarly = oldSum;

      if (!riseStartBin && runningSum > riseStartSentry)
        riseStartBin = i - (runningSum - riseStartSentry) / (runningSum - oldSum);

      if (!riseEndBin && runningSum > riseEndSentry)
        riseEndBin = i - (runningSum - riseEndSentry) / (runningSum - oldSum);

      if (!fallStartBin && runningSum > fallStartSentry)
        fallStartBin = i - (runningSum - fallStartSentry) / (runningSum - oldSum);

      if (!fallEndBin && runningSum > fallEndSentry)
        fallEndBin = i - (runningSum - fallEndSentry) / (runningSum - oldSum);

      if (!t50Bin && runningSum > t50Sentry)
        t50Bin = i - (runningSum - t50Sentry) / (runningSum - oldSum);

      oldSum = runningSum;

    }

    pmtRec.SetPeakAmplitude(peakAmplitude);
    pmtRec.SetRiseTime(binTiming * (riseEndBin - riseStartBin), 0);
    pmtRec.SetFallTime(binTiming * (fallEndBin - fallStartBin), 0);
    pmtRec.SetT50(binTiming * (t50Bin - startBin));
    const double sumLate = runningSum - sumEarly;
    if (shapeSentry < endIntegration && sumLate > 0)
      pmtRec.SetShapeParameter(sumEarly / sumLate);
  }


  TraceD
  SdTopDownSignalSelectorUGR::GetVEMChargeTrace(const sevt::Station& station)
  const
  {

    const sdet::Station& dStation =
      det::Detector::GetInstance().GetSDetector().GetStation(station);

    vector<const PMT*> validPMTs;

    TraceD vemTrace(dStation.GetFADCTraceLength(), dStation.GetFADCBinSize(),
                    sevt::StationConstants::eMySource0);

    const int traceLength = vemTrace.GetSize();
    for (int i = 0; i < traceLength; ++i)
      vemTrace[i] = 0;

    for (Station::ConstPMTIterator pmtIt = station.PMTsBegin();
         pmtIt != station.PMTsEnd(); ++pmtIt)
      if (pmtIt->HasCalibData() && pmtIt->GetCalibData().IsTubeOk() &&
          pmtIt->HasRecData() && pmtIt->GetRecData().HasVEMTrace())
        validPMTs.push_back(&(*pmtIt));

    if (!validPMTs.empty()) {
      const int nPMTs = validPMTs.size();
      for (vector<const PMT*>::const_iterator pmtIt = validPMTs.begin();
           pmtIt != validPMTs.end(); ++pmtIt) {
        const PMTRecData& pmtRec = (*pmtIt)->GetRecData();
        const TraceD& pmtTrace = pmtRec.GetVEMTrace();
        const double peakChargeFactor = pmtRec.GetVEMPeak() / pmtRec.GetVEMCharge();
        for (int i = 0; i < traceLength; ++i)
          vemTrace[i] += (pmtTrace[i] * peakChargeFactor) / nPMTs;
      }
    }

    return vemTrace;
  }


  GeoSegmentCollection
  SdTopDownSignalSelectorUGR::FoundBasicSegments(sevt::Station& st)
  {
    // For station detector
    const sdet::Station& dStation =
      det::Detector::GetInstance().GetSDetector().GetStation(st);

    // For start time of the signal
    const double fadcBinSize = dStation.GetFADCBinSize();
    const StationGPSData& gpsData = st.GetGPSData();
    const TimeStamp gpsTime(gpsData.GetSecond(), gpsData.GetCorrectedNanosecond());
    const TimeStamp traceTime = gpsTime -
                                TimeInterval(dStation.GetFADCTraceLength() * fadcBinSize);

    GeoSegmentCollection segments;

    const utl::TraceD& vemChargeTrace = GetVEMChargeTrace(st);
    const utl::TraceD& vemPeakTrace = st.GetVEMTrace();

    //sevt::Station::Signal subsegment(0,0,0,0);
    GeoSegment subsegment(0, 0, 0, 0., 0., 0., traceTime, &st, &dStation, true);

    double lowFactor = 1;
    if (st.IsHighGainSaturation()) {
      lowFactor *= fLowFactor;
      if (fVerbose)
        cout << "  Station with High gain trace saturated. Using factor " << lowFactor << " for signal threshold." << endl;
    }

    if (fVerbose > 2)
      cout << "\n\t\t Scaning for segments..." << endl;

    // Search for elemental subsegments
    for (utl::TraceD::SizeType bin = vemChargeTrace.GetStart(); bin < vemChargeTrace.GetStop(); ++bin) {
      if (fVerbose > 4)
        cout << "\t\t\t Bin: " << bin
             << "\t Value: " << vemPeakTrace[bin]
             << " VEM Peak\t Value: " << vemChargeTrace[bin] << " VEM Charge" << endl;
      if (vemChargeTrace[bin] < fPreSignalThreshold * lowFactor && subsegment.fBinsOverThresh == 0)
        continue;
      else if (vemChargeTrace[bin] < fPreSignalThreshold * lowFactor && subsegment.fBinsOverThresh != 0) {
        if (fVerbose > 2)
          cout << "\t\t Adding new subsegment (start, stop, charge, bins): ("
               << subsegment.fStart << " , "
               << subsegment.fStop  << " , "
               << subsegment.fCharge  << " , "
               << subsegment.fAoP  << " , "
               << subsegment.fBinsOverThresh  << ")" << endl;

        //      if((subsegment.fBinsOverThresh>=fMinSlots) && (subsegment.fCharge>=fMinSignal))
        segments.push_back(subsegment); // Store one segment
        // Init subsegment
        subsegment.fBinsOverThresh = 0;
        subsegment.fCharge = 0;
        subsegment.fAoP = 0;
        subsegment.fSignal = 0;
        subsegment.fStart = 0;
        subsegment.fStop = 0;
        subsegment.fStartTime = 0;
      } else if (vemChargeTrace[bin] >= fPreSignalThreshold * lowFactor && subsegment.fBinsOverThresh == 0) {
        subsegment.fStart = bin;
        subsegment.fStop = bin;
        subsegment.fCharge += vemChargeTrace[bin];
        subsegment.fSignal += vemPeakTrace[bin];
        ++subsegment.fBinsOverThresh;
        subsegment.fStartTime = traceTime + TimeInterval((bin - 0.5) * fadcBinSize);
        if (fVerbose >= 4)
          cout << "\t\t\t New segment starting: "
               << subsegment.fStart << " , "
               << subsegment.fStop << " , "
               << subsegment.fCharge << " , "
               << subsegment.fBinsOverThresh << endl;
      } else if (vemChargeTrace[bin] >= fPreSignalThreshold * lowFactor && subsegment.fBinsOverThresh != 0) {
        subsegment.fStop = bin;
        subsegment.fCharge += vemChargeTrace[bin];
        ++subsegment.fBinsOverThresh;
        subsegment.fSignal += vemPeakTrace[bin];
        if (fVerbose >= 4)
          cout << "\t\t\t Updating segment: "
               << subsegment.fStart << " , "
               << subsegment.fStop << " , "
               << subsegment.fCharge << " , "
               << subsegment.fBinsOverThresh << endl;
      }
    }

    if (segments.size() <= 1)
      return segments; // Only one segment or no segments

    if (fVerbose > 2)
      cout << "\n\t\t Merging subsegments..." << endl;
    if (fVerbose > 3)
      cout << "\t\t\t(Start, Stop, Charge, BinsOverThreshold)" << endl;

    // Merge subsegments into segments
    for (unsigned int seg = 0; seg < segments.size() - 1; ++seg) {
      // if(seg==(segments.size()-2)) break;
      if ((segments[seg].fStop + fPreMingap + segments[seg].fBinsOverThresh) > segments[seg + 1].fStart) { // Merge next segment
        if (fVerbose > 3)
          cout << "\t\t\t Merging: ("
               << segments[seg].fStart << " , "
               << segments[seg].fStop << " , "
               << segments[seg].fCharge << " , "
               << segments[seg].fBinsOverThresh << ") + ("
               << segments[seg + 1].fStart << " , "
               << segments[seg + 1].fStop << " , "
               << segments[seg + 1].fCharge << " , "
               << segments[seg + 1].fBinsOverThresh << ") = (";

        segments[seg].fStop = segments[seg + 1].fStop;
        segments[seg].fBinsOverThresh += segments[seg + 1].fBinsOverThresh;
        segments[seg].fCharge += segments[seg + 1].fCharge;
        segments[seg].fSignal += segments[seg + 1].fSignal;

        if (fVerbose > 3)
          cout << segments[seg].fStart << " , "
               << segments[seg].fStop << " , "
               << segments[seg].fCharge << " , "
               << segments[seg].fBinsOverThresh << ")" << endl;

        segments.erase(segments.begin() + seg + 1);
        --seg;
      }
    }

    if (fVerbose > 2)
      cout << "\n\t\t Deleting Pre segments below threshold... " << endl;
    if (fVerbose > 3)
      cout << "\t\t\t(Start, Stop, Charge, BinsOverThreshold)" << endl;
    // Check for acceptable segments
    for (unsigned int seg = 0; seg < segments.size(); ++seg) {
      if (segments[seg].fBinsOverThresh >= fPreMinSlots && segments[seg].fCharge > fPreMinSignal) { // Good segment
        // Compute AoP
        double peak = 0;
        for (int bin = segments[seg].fStart; bin <= segments[seg].fStop; ++bin)
          peak = vemPeakTrace[bin] > peak ? vemPeakTrace[bin] : peak;
        segments[seg].fAoP = segments[seg].fCharge / peak;
        continue;
      }

      if (fVerbose > 3)
        cout << "\t\t\t Deleting segment: ("
             << segments[seg].fStart << " , "
             << segments[seg].fStop << " , "
             << segments[seg].fCharge << " , "
             << segments[seg].fBinsOverThresh << ")" << endl;

      segments.erase(segments.begin() + seg);
      --seg;
    }

    if (fVerbose > 1) {
      cout << "\n\t Selected segments: (Start, Stop, VEMCharge, AoP, VEMPeakSignal, BinsOverTh, SignalStartTime)" << endl;
      for (unsigned int seg = 0; seg < segments.size(); ++seg) {
        cout << "\t\t ("
             << segments[seg].fStart << " , "
             << segments[seg].fStop << " , "
             << segments[seg].fCharge << " , "
             << segments[seg].fAoP << " , "
             << segments[seg].fSignal << " , "
             << segments[seg].fBinsOverThresh << " , "
             << segments[seg].fStartTime.GetGPSNanoSecond() << ")" << endl;
      }
    }

    return segments;
  }


  void
  SdTopDownSignalSelectorUGR::UpdateSegmentValues(const int seg)
  {
    sevt::Station& st = *(fAllSegments[seg].fStEvt);

    if (fVerbose > 2)
      cout << "\n\t\t Updating segment for station: " << st.GetId()
           << "\n\t\t (Start, Stop): (" << fAllSegments[seg].fStart
           << " , " << fAllSegments[seg].fStop << endl;

    GeoSegmentCollection segments;

    const utl::TraceD& vemChargeTrace = GetVEMChargeTrace(st);
    utl::TraceD& vemPeakTrace = st.GetVEMTrace();

    //sevt::Station::Signal subsegment(0,0,0,0);
    GeoSegment subsegment(0, 0, 0, 0., 0., 0., 0, NULL, NULL, true);

    double lowFactor = 1;
    if (st.IsHighGainSaturation()) {
      if (fVerbose)
        cout << "  Station with High gain trace saturated. Using factor " << lowFactor << " for signal threshold." << endl;
      lowFactor *= fLowFactor;
    }

    if (fVerbose > 2)
      cout << "\n\t\t Scaning for segments..." << endl;

    // Search for elemental subsegments
    for (utl::TraceD::SizeType bin = fAllSegments[seg].fStop; bin < vemChargeTrace.GetStop(); ++bin) {

      if ((fAllSegments[seg].fStop + fMingap + fAllSegments[seg].fBinsOverThresh) < int(bin) &&
          !segments.size() && !subsegment.fStart) {
        if (fVerbose > 2)
          cout << "\t\tNo more segments found " << endl;
        return;
      }

      if (fVerbose > 4)
        cout << "\t\t\t Bin: " << bin
             << "\t Value: " << vemPeakTrace[bin]
             << " VEM Peak\t Value: " << vemChargeTrace[bin] << " VEM Charge" << endl;
      if (vemChargeTrace[bin] < fSignalThreshold * lowFactor && subsegment.fBinsOverThresh == 0)
        continue;
      else if (vemChargeTrace[bin] < fSignalThreshold * lowFactor && subsegment.fBinsOverThresh != 0) {
        if (fVerbose > 2)
          cout << "\t\t Adding new subsegment (start, stop, charge, bins): ("
               << subsegment.fStart << " , "
               << subsegment.fStop  << " , "
               << subsegment.fCharge  << " , "
               << subsegment.fAoP  << " , "
               << subsegment.fBinsOverThresh  << ")" << endl;

        segments.push_back(subsegment); // Store one segment
        // Init subsegment
        subsegment.fBinsOverThresh = 0;
        subsegment.fCharge = 0;
        subsegment.fAoP = 0;
        subsegment.fSignal = 0;
        subsegment.fStart = 0;
        subsegment.fStop = 0;
        subsegment.fStartTime = 0;
      } else if (vemChargeTrace[bin] >= fSignalThreshold * lowFactor && subsegment.fBinsOverThresh == 0) {
        subsegment.fStart = bin;
        subsegment.fStop = bin;
        subsegment.fCharge += vemChargeTrace[bin];
        subsegment.fSignal += vemPeakTrace[bin];
        ++subsegment.fBinsOverThresh;

        //If first subsegment founded is far enought return
        if ((fAllSegments[seg].fStop + fMingap + fAllSegments[seg].fBinsOverThresh) < int(bin)) {
          if (fVerbose > 2)
            cout << "\t\tNo more segments found " << endl;
          return;
        }

        if (fVerbose >= 4)
          cout << "\t\t\t New segment starting: "
               <<  subsegment.fStart << " , "
               <<  subsegment.fStop << " , "
               <<  subsegment.fCharge << " , "
               <<  subsegment.fBinsOverThresh << endl;
      } else if (vemChargeTrace[bin] >= fSignalThreshold * lowFactor && subsegment.fBinsOverThresh != 0) {
        subsegment.fStop = bin;
        subsegment.fCharge += vemChargeTrace[bin];
        ++subsegment.fBinsOverThresh;
        subsegment.fSignal += vemPeakTrace[bin];
        if (fVerbose >= 4)
          cout << "\t\t\t Updating segment: "
               << subsegment.fStart << " , "
               << subsegment.fStop << " , "
               << subsegment.fCharge << " , "
               << subsegment.fBinsOverThresh << endl;
      }
    }

    // Merge subsegments into segments
    for (unsigned int i = 0; i < segments.size(); ++i) {
      if ((fAllSegments[seg].fStop + fMingap + fAllSegments[seg].fBinsOverThresh) > segments[i].fStart) { // Merge next segment
        if (fVerbose > 3)
          cout << "\t\t\t Merging: ("
               << fAllSegments[seg].fStart << " , "
               << fAllSegments[seg].fStop << " , "
               << fAllSegments[seg].fCharge << " , "
               << fAllSegments[seg].fBinsOverThresh << ") + ("
               << segments[i].fStart << " , "
               << segments[i].fStop << " , "
               << segments[i].fCharge << " , "
               << segments[i].fBinsOverThresh << ") = (";

        fAllSegments[seg].fStop = segments[i].fStop;
        fAllSegments[seg].fBinsOverThresh += segments[i].fBinsOverThresh;
        fAllSegments[seg].fCharge += segments[i].fCharge;
        fAllSegments[seg].fSignal += segments[i].fSignal;

        if (fVerbose > 3)
          cout << fAllSegments[seg].fStart << " , "
               << fAllSegments[seg].fStop << " , "
               << fAllSegments[seg].fCharge << " , "
               << fAllSegments[seg].fBinsOverThresh << ")" << endl;
      }
    }

    // Recompute  AoP
    {
      double peak = 0;
      for (int bin = fAllSegments[seg].fStart; bin <= fAllSegments[seg].fStop; ++bin)
        peak = vemPeakTrace[bin] > peak ? vemPeakTrace[bin] : peak;
      fAllSegments[seg].fAoP = fAllSegments[seg].fCharge / peak;
    }

    if (fVerbose > 3)
      cout << "\t\t\t Merging done! " << endl;
  }

}