EtsFixing.cxx

Go to the documentation of this file.

#include "ResetEtsAlg/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();
    //}
}