EtsFixing.cxx

Go to the documentation of this file.

#include "EtsFixing.h"
#include "TFile.h"
#include "TTree.h"
 
static const int  NLastBackup    = 10;
static const long OneSecVar      = 2000000;
static const long ShiftThreshold = OneSecVar * 0.5;
 
// static int    _debug = 1;
// static TFile* _file = NULL;
// static TTree* _tree = NULL;
// static int    _count = 0;
// static int    _event;
// static int    _time;
// static int    _pileup;
// static long   _t1;
// static long   _rawt1;
// static long   _tdiff = 0;
// static float  _slope;
 
EtsFixing::EtsFixing( int preRefCount, double factor )
    : m_nPre( preRefCount )
    , m_factor( factor )
    , m_t0Sec( 0 )
    , m_t0SecLocal( 0 )
    , m_count( 0 )
    , m_refSlope( 0.0 )
    , m_critical( OneSecVar + 999 )
    , m_pileup( 0 ) {
  m_preEvt = new int[m_nPre];
  m_preSec = new int[m_nPre];
  m_preT1  = new long[m_nPre];
 
  m_preT1Sec  = new int[m_nPre];
  m_preT1Old  = new long[m_nPre];
  m_prePileup = new long[m_nPre];
 
  for ( int i = 0; i < m_nPre; ++i )
  {
    m_preEvt[i] = -999;
    m_preSec[i] = -999;
    m_preT1[i]  = -999;
 
    m_preT1Sec[i]  = -999;
    m_preT1Old[i]  = -999;
    m_prePileup[i] = -999;
  }
 
  m_lastEvt = new int[NLastBackup];
  m_lastT1  = new long[NLastBackup];
 
  for ( int i = 0; i < NLastBackup; ++i )
  {
    m_lastEvt[i] = -999;
    m_lastT1[i]  = -999;
  }
 
  // if (_debug) {
  //     _file = new TFile("fff.root", "RECREATE");
  //     _tree = new TTree("ff", "for test");
  //     _tree->Branch("count", &_count);
  //     _tree->Branch("event", &_event);
  //     _tree->Branch("time",  &_time);
  //     _tree->Branch("pileup",&_pileup);
  //     _tree->Branch("t1",    &_t1, "t1/L");
  //     _tree->Branch("rawt1", &_rawt1, "rawt1/L");
  //     _tree->Branch("tdiff", &_tdiff, "tdiff/L");
  //     _tree->Branch("slope", &_slope, "slope/F");
  // }
}
 
EtsFixing::~EtsFixing() {
  delete[] m_preEvt;
  delete[] m_preSec;
  delete[] m_preT1;
 
  delete[] m_lastEvt;
  delete[] m_lastT1;
 
  delete[] m_preT1Sec;
  delete[] m_preT1Old;
  delete[] m_prePileup;
 
  // if (_debug) {
  //     _file->Write();
  //     _file->Close();
  //     delete _file;
  // }
}
 
void EtsFixing::fixT1( Event::EventHeader* header ) {
  long t1Old = header->etsT1();
 
  if ( t1Old == 0 ) { return; }
 
  int  evtNo     = header->eventNumber();
  int  tnSec     = header->time() - m_t0Sec;
  long t1NanoSec = t1Old % OneSecVar;
  long t1New     = t1NanoSec; // tnSec*OneSecVar + t1NanoSec;
 
  header->setRawEtsT1( t1Old );
 
  if ( m_count != 0 )
  {
 
    int curIdx = tnSec % m_nPre;
 
    int preIdx   = curIdx;
    int secShift = 0;
    if ( tnSec != m_preSec[preIdx] )
    {
      preIdx   = ( preIdx + m_nPre - 1 ) % m_nPre;
      secShift = 1;
      while ( m_preEvt[preIdx] == -999 ||
              m_preSec[preIdx] < m_preSec[( preIdx + 1 ) % m_nPre] )
      {
        preIdx = ( preIdx + m_nPre - 1 ) % m_nPre;
        ++secShift;
      }
      if ( ( m_preSec[preIdx] + secShift ) <= ( tnSec - m_nPre ) )
      { secShift += m_nPre * ( ( tnSec - secShift - m_preSec[preIdx] ) / m_nPre ); }
      // if (_debug > 2) std::cout << tnSec << ", " << m_preSec[preIdx] << ", " << secShift <<
      // std::endl;
    }
    t1New += ( m_preT1[preIdx] / OneSecVar + secShift ) * OneSecVar;
    // if (_debug > 1) std::cout << __LINE__ << ": " << t1New << std::endl;
 
    /// fix with ideal slope
    int  eDiff    = evtNo - m_preEvt[preIdx];
    long tDiff    = t1New - m_preT1[preIdx];
    long t1Expect = m_preT1[preIdx] + m_refSlope * eDiff;
    long t1Shift  = t1New - t1Expect;
    if ( t1Shift > ShiftThreshold )
    {
      t1New -= OneSecVar;
      tDiff -= OneSecVar;
    }
    else if ( t1Shift < -ShiftThreshold )
    {
      t1New += OneSecVar;
      tDiff += OneSecVar;
    }
    // if (_debug) _tdiff = t1New - t1Expect;
    // if (_debug > 1) std::cout << __LINE__ << ": " << t1New << std::endl;
 
    /// adjust with event build time
    if ( tnSec - m_t0SecLocal > 3 )
    {
      long tXX = long( tnSec ) * OneSecVar - t1New -
                 long( tnSec - t1Old / OneSecVar ) * OneSecVar * ( m_factor - 1 );
      // if (_debug > 2) std::cout << "xxx: " << tnSec << ", " << m_critical << ", " << tXX <<
      // std::endl; expected range: (-m_critical, OneSecVar-m_critical], with a toleration of
      // 0.04s
      while ( tXX < -m_critical - OneSecVar * 0.04 )
      {
        t1New -= OneSecVar;
        tDiff -= OneSecVar;
        tXX += OneSecVar;
      }
      while ( tXX > OneSecVar * 1.04 - m_critical )
      {
        t1New += OneSecVar;
        tDiff += OneSecVar;
        tXX -= OneSecVar;
      }
    }
    else
    {
      if ( secShift == 1 )
      {
        m_critical = t1New - long( tnSec - 1 ) * OneSecVar;
        // if (_debug > 2) std::cout << "zzz: " << tnSec << ", " << t1New << ", " << t1Old <<
        // std::endl;
      }
    }
    // if (_debug > 1) std::cout << __LINE__ << ": " << t1New << std::endl;
 
    /// fix with neighboring event
    bool noNeighbor = true;
    int  nearestIdx = NLastBackup - 1;
    int  nearestDE  = 10000000;
    for ( int i = nearestIdx; i >= 0; --i )
    {
      int lidx   = ( m_count + i ) % NLastBackup;
      int lediff = evtNo - m_lastEvt[lidx];
      if ( abs( lediff ) < 20 )
      {
        long ltdiff = t1New - m_lastT1[lidx];
        /// the time difference should not be too large for very near (event id) events
        while ( ltdiff > ShiftThreshold )
        {
          t1New -= OneSecVar;
          tDiff -= OneSecVar;
          ltdiff -= OneSecVar;
        }
        while ( ltdiff < -ShiftThreshold )
        {
          t1New += OneSecVar;
          tDiff += OneSecVar;
          ltdiff += OneSecVar;
        }
        if ( lediff > 0 )
        {
          while ( ltdiff < 0 )
          {
            t1New += OneSecVar;
            tDiff += OneSecVar;
            ltdiff += OneSecVar;
          }
        }
        else
        {
          while ( ltdiff > 0 )
          {
            t1New -= OneSecVar;
            tDiff -= OneSecVar;
            ltdiff -= OneSecVar;
          }
        }
        noNeighbor = false;
        break;
      }
      else
      {
        if ( abs( lediff ) < nearestDE ) { nearestIdx = i; }
      }
    }
    // if (_debug > 1) std::cout << __LINE__ << ": " << t1New << std::endl;
 
    /// use nearest event if no neighboring event is found
    if ( noNeighbor )
    {
      int  lidx   = ( m_count + nearestIdx ) % NLastBackup;
      int  lediff = evtNo - m_lastEvt[lidx];
      long ltdiff = t1New - m_lastT1[lidx];
      /// in case of negative or abnormal slope
      if ( lediff > 0 )
      {
        while ( ltdiff * 1.0 / lediff < m_refSlope * 0.2 )
        {
          t1New += OneSecVar;
          tDiff += OneSecVar;
          ltdiff += OneSecVar;
          if ( ltdiff > OneSecVar && ltdiff * 1.0 / lediff > m_refSlope * 5 )
          {
            t1New -= OneSecVar;
            tDiff -= OneSecVar;
            ltdiff -= OneSecVar;
            break;
          }
        }
      }
      else
      {
        while ( ltdiff * 1.0 / lediff > m_refSlope * 5 )
        {
          t1New += OneSecVar;
          tDiff += OneSecVar;
          ltdiff += OneSecVar;
        }
        while ( ltdiff > 0 )
        {
          t1New -= OneSecVar;
          tDiff -= OneSecVar;
          ltdiff -= OneSecVar;
        }
        while ( ltdiff * 1.0 / lediff < m_refSlope * 0.2 )
        {
          t1New -= OneSecVar;
          tDiff -= OneSecVar;
          ltdiff -= OneSecVar;
        }
      }
    }
    // if (_debug > 1) std::cout << __LINE__ << ": " << t1New << std::endl;
 
    /// update slope with reliable events
    float curSlope = ( 1.0 * tDiff ) / eDiff;
    if ( curSlope > 0 )
    {
      if ( m_count > 30 )
      {
        float guard = curSlope / m_refSlope;
        if ( guard > 0.75 && guard < 1.5 ) { m_refSlope = m_refSlope * 0.9 + curSlope * 0.1; }
      }
      else { m_refSlope = ( m_refSlope * ( m_count - 1 ) + curSlope ) / m_count; }
    }
 
    /// update parameters for future event
    /// Note: DO NOT use the finalT1 here. t1New and finalT1 are not consistent
    if ( secShift == 1 )
    {
      // do not update these parameters when an event comes earlier more than 1 second
      m_preEvt[curIdx] = evtNo;
      m_preSec[curIdx] = tnSec;
      m_preT1[curIdx]  = t1New;
    }
    m_lastEvt[m_count % NLastBackup] = evtNo;
    m_lastT1[m_count % NLastBackup]  = t1New;
 
    /// Refined fixing of ETS with the factor (1.0000115 by default)
    if ( t1New != t1Old )
    {
      int t1Sec = t1New / OneSecVar;
      int fIdx  = t1Sec % m_nPre;
      if ( t1Sec == m_preT1Sec[fIdx] ) { m_pileup = m_prePileup[fIdx]; }
      else
      {
        for ( int i = 1; i < m_nPre; ++i )
        {
          int pT1Sec = t1Sec - i;
          int pFIdx  = ( pT1Sec + m_nPre ) % m_nPre;
          if ( pT1Sec == m_preT1Sec[pFIdx] )
          {
            int fTotalPileup = t1Sec - t1Old / OneSecVar;
            int pTotalPileup = m_preT1Sec[pFIdx] - m_preT1Old[pFIdx] / OneSecVar;
            if ( fTotalPileup != pTotalPileup )
            {
              m_pileup = m_prePileup[pFIdx] + long( fTotalPileup - pTotalPileup ) * OneSecVar;
            }
            else
            { // the GPS clock is recovered
              m_pileup = 0;
            }
            break;
          }
        }
        m_preT1Sec[fIdx]  = t1Sec;
        m_preT1Old[fIdx]  = t1Old;
        m_prePileup[fIdx] = m_pileup;
      }
      long finalT1 = t1New + ( m_factor - 1.0 ) * m_pileup;
      if ( finalT1 > 0 ) { header->setEtsT1( finalT1 ); }
      else { header->setEtsT1( 0 ); }
    }
 
    // if (_debug > 1) std::cout << "event " << evtNo << " : "
    //     << header->etsT1() << " - " << header->rawEtsT1()
    //     << " = " << long(header->etsT1()) - long(header->rawEtsT1())
    //     << ", " << m_refSlope
    //     << std::endl;
  }
  else
  { // first event
    int firstSec       = t1Old / OneSecVar;
    int firstIdx       = firstSec % m_nPre;
    m_t0Sec            = header->time() - firstSec;
    m_t0SecLocal       = firstSec;
    m_preEvt[firstIdx] = evtNo;
    m_preSec[firstIdx] = firstSec;
    m_preT1[firstIdx]  = t1Old;
 
    m_lastEvt[0] = evtNo;
    m_lastT1[0]  = t1Old;
 
    m_pileup              = 0;
    m_preT1Sec[firstIdx]  = firstSec;
    m_preT1Old[firstIdx]  = t1Old;
    m_prePileup[firstIdx] = 0;
  }
 
  ++m_count;
 
  // if (_debug) {
  //     _count = m_count;
  //     _event = header->eventNumber();
  //     _time  = header->time();
  //     _pileup = m_pileup/OneSecVar;
  //     _t1    = header->etsT1();
  //     _rawt1 = header->rawEtsT1();
  //     _slope = m_refSlope;
  //     _tree->Fill();
  // }
}