ZddReconAlg.cxx

Go to the documentation of this file.

#include "ZddReconAlg/ZddReconAlg.h"
#include "EvTimeEvent/RecEsTime.h"
#include "EventModel/EventHeader.h"
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/SmartDataPtr.h"
#include "ReconEvent/ReconEvent.h"
#include "ZddEvent/RecZddChannel.h"
#include "ZddEvent/ZddBoard.h"
 
using namespace std;
 
DECLARE_COMPONENT( ZddReconAlg )
/////////////////////////////////////////////////////////////////////////////
 
ZddReconAlg::ZddReconAlg( const std::string& name, ISvcLocator* pSvcLocator )
    : Algorithm( name, pSvcLocator ), m_errStat( false ) {
  declareProperty( "BreakWithError", m_errQuit = true );
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
StatusCode ZddReconAlg::initialize() {
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "in initialize()" << endmsg;
 
  return StatusCode::SUCCESS;
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
StatusCode ZddReconAlg::execute() {
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "in execute()" << endmsg;
 
  // Part 1: Get the event header, print out event and run number
  SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
  if ( !eventHeader )
  {
    log << MSG::FATAL << "Could not find Event Header" << endmsg;
    return StatusCode::FAILURE;
  }
  log << MSG::DEBUG << "Retrieved event: " << eventHeader->eventNumber()
      << "  run: " << eventHeader->runNumber() << endmsg;
 
  // Part 2: Register RecZddChannel to TDS
  DataObject* aRecEvent = 0;
  eventSvc()->findObject( "/Event/Recon", aRecEvent );
  if ( aRecEvent == 0 )
  {
    aRecEvent     = new ReconEvent();
    StatusCode sc = eventSvc()->registerObject( "/Event/Recon", aRecEvent );
    if ( sc.isFailure() )
    {
      log << MSG::FATAL << "Could not register ReconEvent" << endmsg;
      return StatusCode::FAILURE;
    }
  }
 
  RecZddChannelCol* recZddCol = new RecZddChannelCol;
  StatusCode sc = eventSvc()->registerObject( EventModel::Recon::RecZddChannelCol, recZddCol );
  if ( sc.isFailure() )
  {
    log << MSG::FATAL << "Could not register RecZddChannelCol" << endmsg;
    return StatusCode::FAILURE;
  }
 
  // Part 3: check the ZDD data
  // 3.1: some errors happened before, ignore the rest calculations
  if ( m_errStat ) return StatusCode::SUCCESS;
 
  // 3.2: check the status of the current ZDD event
  SmartDataPtr<Event::ZddEvent> zddEvt( eventSvc(), "/Event/ZddEvent" );
  int zddCheck = zddDataStat( zddEvt.ptr(), eventHeader->eventNumber() + 1 );
 
  if ( zddCheck != 0 )
  {
    if ( zddCheck < 0 )
    {
      // serious error happened, ignore all of the rest ZDD data
      m_errStat = true;
 
      if ( m_errQuit )
      { // whether to break the data processing
        return StatusCode::FAILURE;
      }
    }
    // else:
    // problems in the current event, only ignore this event
    return StatusCode::SUCCESS;
  }
 
  // Part 4: Get event T0
  double                     bes3_t0 = -10000.0;
  SmartDataPtr<RecEsTimeCol> evTimeCol( eventSvc(), "/Event/Recon/RecEsTimeCol" );
  if ( !evTimeCol || evTimeCol->size() == 0 )
  {
    log << MSG::WARNING << " Could not find RecEsTimeCol" << endmsg;
    // return StatusCode::SUCCESS;
  }
  else
  {
    RecEsTimeCol::iterator iter_evt = evTimeCol->begin();
    if ( iter_evt != evTimeCol->end() )
    {
      bes3_t0 = ( *iter_evt )->getTest();
      // std::cout << "BESIII t0: " << bes3_t0 << std::endl; //wangyadi
      bes3_t0 = 6400 - bes3_t0; // T_L1 - T_collision
    }
  }
 
  // Part 5: ZDD data processing
  const Event::ZddEvent::Channels& chs = zddEvt->channels();
  // cout << "ZDD has n channel: " << chs.size() << endl;
  //  loop the channels
  Event::ZddEvent::Channels::const_iterator end_ch = chs.end();
  for ( Event::ZddEvent::Channels::const_iterator it = chs.begin(); it != end_ch; ++it )
  {
    const ZddChannel*            pch   = *it;
    const ZddChannel::Fragments& frags = pch->fragments();
    // cout << "Channel " << pch->getChId() << "  has " << frags.size() << " fragments: " <<
    // endl;
 
    RecZddChannel* recZddCh = new RecZddChannel;
    recZddCh->setChannelId( pch->getChId() );
    recZddCh->setScanCode( pch->getScanCode() );
    double e_K = getEK( pch->getChId() );
 
    // loop the fragments, rebuild the original channel waveform into one buffer
    int           maxSamples = 800; // This value may change!
    unsigned char waveform[800];
    memset( waveform, 0, maxSamples );
    ZddChannel::Fragments::const_iterator end_fg     = frags.end();
    bool                                  quit_event = false;
    int                                   start      = 0;
    for ( ZddChannel::Fragments::const_iterator jt = frags.begin(); jt != end_fg; ++jt )
    {
      const ZddFragment& frag = *jt;
      // cout << "\tindex from " << frag.start_index << ",\t length " << frag.length << ":\t ";
      // for ( int i = 0; i < frag.length; ++i ) {
      //     cout << int(frag.sample[i]) << " ";
      // }
      // cout << endl;
      start = frag.start_index;
 
      // check for CAEN data corruption ( absent from 2013 onward ?? )
      if ( start + frag.length > maxSamples )
      {
        recZddCh->setScanCode( 2 * 10 + pch->getScanCode() );
        MsgStream log( msgSvc(), name() );
        log << MSG::ERROR << "ZDD BAD DATA: CAEN corruption problem" << endmsg;
        quit_event = true; // abandon this event...
        break;             // ...starting from this channel
      }
 
      for ( int i = 0; i < frag.length; ++i ) waveform[start++] = frag.sample[i];
    }
 
    // then reclusterize with a slightly higher threshold
    unsigned char threshold        = 20; // threshold for reclustering
    unsigned char rephaseThreshold = 40; // threshold for waveform rephasing
    unsigned char minSample = 255, maxSample = -1;
    int           maxTime   = -1;
    bool          closed    = true;               // no cluster is running yet
    int           phases[4] = { -1, -1, -1, -1 }; // a 4-channel histogram
    for ( int pt = 0; pt < maxSamples; pt++ )
    {
      bool notZero = waveform[pt] > 0;
      bool smaller = waveform[pt] < minSample;
      if ( notZero && smaller ) minSample = waveform[pt]; // find baseline for whole channel
      if ( waveform[pt] > threshold )
      {
        if ( closed )
        { // start a new cluster, initialize max and index
          closed    = false;
          maxSample = waveform[pt];
          maxTime   = pt;
        }
        else
        { // continue the old cluster, update max and index
          if ( waveform[pt] > maxSample )
          {
            maxSample = waveform[pt];
            maxTime   = pt;
          }
        }
      }
      else
      {
        if ( !closed )
        { // close the old cluster and store it
          closed       = true;
          double tNsec = 2. * maxTime;
          recZddCh->addFragment( tNsec,
                                 maxSample * e_K ); // time relative to start of FADC window
          // std::cout << "       ZDD E & T: " << int(maxSample) << ", " << maxTime <<
          // std::endl;
          if ( maxSample > rephaseThreshold )
          {                          // most peaks are at multiples of 8ns
            int phase = maxTime % 4; // only one of these bins will be most used
            phases[phase]++;         // but we do not know which one yet
          }
        }
      }
    }
    if ( !closed )
    { // close and store the last cluster if still running
      closed       = true;
      double tNsec = 2. * maxTime; // time relative to start of FADC window
      recZddCh->addFragment( tNsec, maxSample * e_K );
      // std::cout << "       ZDD E & T: " << int(maxSample) << ", " << maxTime << std::endl;
      if ( maxSample > rephaseThreshold )
      {
        int phase = maxTime % 4;
        phases[phase]++;
      }
    }
 
    // Compute channel phase
    int mostProb = -1;
    int chPhase  = -1;
    for ( int ph = 0; ph < 4; ph++ )
    {
      if ( phases[ph] > mostProb )
      {
        mostProb = phases[ph];
        chPhase  = ph;
      }
    }
 
    if ( chPhase == -1 )
    { // mark channel as not rephasable (no samples above rephaseThreshold)
      recZddCh->setScanCode( 3 * 10 + pch->getScanCode() );
      chPhase = -2;
    }
 
    recZddCh->setBaseLine( minSample );
    recZddCh->setPhase( chPhase );
    recZddCol->push_back( recZddCh );
 
    /* Previous code, left for reference
            // loop the fragments
            ZddChannel::Fragments::const_iterator end_fg = frags.end();
            for ( ZddChannel::Fragments::const_iterator jt = frags.begin(); jt != end_fg; ++jt)
       { const ZddFragment& frag = *jt; int efrag = -1, tfrag = -1; calFragEandT(frag, efrag,
       tfrag); recZdd->addFragment(bes3_t0-tfrag, efrag*e_K);
            }
    */
 
    if ( quit_event )
      break; // abandon the event entirely after a CAEN data corruption instance
 
  } // end of loop on channels
  return StatusCode::SUCCESS;
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
StatusCode ZddReconAlg::finalize() {
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "in finalize()" << endmsg;
 
  return StatusCode::SUCCESS;
}
 
double ZddReconAlg::getEK( int chId ) { return 1.0 / 16.0; }
 
int ZddReconAlg::zddDataStat( const Event::ZddEvent* zddEvt, int evtId ) {
  MsgStream log( msgSvc(), name() );
 
  if ( !zddEvt )
  {
    log << MSG::FATAL << "Could not find ZddEvent" << endmsg;
    return 1; /*-1;*/
  }
 
  // check the ZDD trigger number
  const std::vector<ZddBoard*>& boards = zddEvt->boards();
  if ( boards.size() != 2 )
  {
    log << MSG::FATAL << "incomplete ZDD data, no more ZDD data this run!" << endmsg;
    return -2;
  }
 
  if ( boards[0]->getCounter() != evtId || boards[1]->getCounter() != evtId )
  {
    log << MSG::FATAL << "missaligned ZDD triggers, no more ZDD data this run!" << endmsg;
    return -3;
  }
 
  /* This check is not useful anymore
      // check the Trigger Time Tag
      static unsigned int bd1_tt = boards[0]->getTimeTag() - 1026;
      static unsigned int bd2_tt = boards[1]->getTimeTag() - 1026;
 
      unsigned int tt0 = bd1_tt;
      unsigned int tt1 = bd2_tt;
      bd1_tt = boards[0]->getTimeTag();
      bd2_tt = boards[1]->getTimeTag();
 
      if ( tt0 > bd1_tt ) tt0 -= 0x80000000;
      if ( tt1 > bd2_tt ) tt1 -= 0x80000000;
 
      if ( (bd1_tt-tt0) < 1026 || (bd2_tt-tt1) < 1026 ) {
          log << MSG::INFO << "overlaped ZDD triggers" << endmsg;
          return 1;
      }
  */
  return 0;
}
 
// This code will be needed when we finally have a calibration
int ZddReconAlg::calFragEandT( const ZddFragment& frag, int& efrag, int& tfrag ) {
  int           minIndex = 0;
  unsigned char min = 255, max = 0;
  for ( int i = 0; i < frag.length; ++i )
  {
    unsigned char& sample = frag.sample[i];
    if ( sample < min )
    {
      min      = sample;
      minIndex = i;
    }
    if ( sample > max ) { max = sample; }
  }
 
  efrag = max - min;
  tfrag = 8160 - 2 * ( minIndex + frag.start_index );
 
  return 0;
}