EstTofCaliSvc.cxx

Go to the documentation of this file.

//********************************************************
//
//  Class: EstTofCaliSvc
//  Sun Shengsen <sunss@ihep.ac.cn>
//
//  Description:
//     This service is used to get Tof Calibrarion constants
//     from DataBase and do the Tof data reconstruction.
//
//*********************************************************
#include "GaudiKernel/DataSvc.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/IIncidentListener.h"
#include "GaudiKernel/IIncidentSvc.h"
#include "GaudiKernel/IInterface.h"
#include "GaudiKernel/Incident.h"
#include "GaudiKernel/Kernel.h"
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/Service.h"
#include "GaudiKernel/SmartDataPtr.h"
#include "GaudiKernel/StatusCode.h"
#include "TMath.h"
#include <fstream>
#include <iostream>
#include <math.h>
 
#include "CalibData/CalibModel.h"
#include "CalibData/Tof/TofCalibData.h"
#include "CalibData/Tof/bTofCalibBase.h"
 
#include "EstTofCaliSvc/EstBTofCal.h"
#include "EstTofCaliSvc/EstBTofCommonCal.h"
#include "EstTofCaliSvc/EstETofCal.h"
#include "EstTofCaliSvc/EstEtfCal.h"
#include "EstTofCaliSvc/EstTofInfoCal.h"
 
#include "EstTofCaliSvc.h"
 
using namespace std;
 
// static double radius_Inner =   81.35;
// static double radius_Outer =   87.05;
DECLARE_COMPONENT( EstTofCaliSvc )
EstTofCaliSvc::EstTofCaliSvc( const std::string& name, ISvcLocator* svcloc )
    : base_class( name, svcloc ) {
  declareProperty( "Run", m_run = 1 );
}
 
StatusCode EstTofCaliSvc::initialize() {
  m_hasbeeninitialized = false;
  m_st                 = false;
 
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << name() << ": Start of run initialisation" << endmsg;
 
  StatusCode sc = Service::initialize();
  if ( sc.isFailure() ) return sc;
 
  IIncidentSvc* incsvc;
  sc           = service( "IncidentSvc", incsvc );
  int priority = 100;
  if ( sc.isSuccess() ) { incsvc->addListener( this, "NewRun", priority ); }
 
  StatusCode scc;
 
  log << MSG::INFO << "setProperties()" << endmsg;
  scc = service( "CalibDataSvc", m_pCalibDataSvc, true );
  if ( !scc.isSuccess() )
  {
    log << MSG::ERROR << "Could not get IDataProviderSvc interface of CalibXmlCnvSvc"
        << endmsg;
    return scc;
  }
  else
  { log << MSG::DEBUG << "Retrieved IDataProviderSvc interface of CalibXmlCnvSvc" << endmsg; }
  // Get properties from the JobOptionsSvc
  // scc = setProperties();
 
  return sc;
}
 
StatusCode EstTofCaliSvc::finalize() {
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << name() << ": End of Run" << endmsg;
  return StatusCode::SUCCESS;
}
 
EstTofCaliSvc::~EstTofCaliSvc() {
  for ( auto it : fBTofCal )
  {
    if ( it ) delete it;
  };
  fBTofCal.clear();
 
  for ( auto it : fETofCal )
  {
    if ( it ) delete it;
  };
  fETofCal.clear();
 
  if ( fEtfCal.size() != 0 )
  {
    for ( auto it : fEtfCal )
    {
      if ( it ) delete it;
    };
    fEtfCal.clear();
  }
 
  for ( auto it : fBTofCommonCal )
  {
    if ( it ) delete it;
  };
  fBTofCommonCal.clear();
 
  for ( auto it : fTofInfoCal )
  {
    if ( it ) delete it;
  };
  fTofInfoCal.clear();
  return;
}
 
StatusCode EstTofCaliSvc::FillfromDatabase() {
 
  MsgStream   log( msgSvc(), name() );
  std::string fullPath = "/Calib/EstTofCal";
  log << MSG::INFO << " Tof calib fullPath = " << fullPath << endmsg;
  SmartDataPtr<CalibData::TofCalibData> test( m_pCalibDataSvc, fullPath );
  if ( !test )
  {
    log << MSG::FATAL << "EstTofCaliSvc could not find TofCalibData in TCDS!!" << endmsg;
    m_st = false;
    return StatusCode::FAILURE;
  }
 
  m_st = true;
 
  for ( auto it : fBTofCal )
  {
    if ( it ) delete it;
  };
  fBTofCal.clear();
 
  for ( auto it : fETofCal )
  {
    if ( it ) delete it;
  };
  fETofCal.clear();
 
  if ( fEtfCal.size() != 0 )
  {
    for ( auto it : fEtfCal )
    {
      if ( it ) delete it;
    };
    fEtfCal.clear();
  }
 
  for ( auto it : fBTofCommonCal )
  {
    if ( it ) delete it;
  };
  fBTofCommonCal.clear();
 
  for ( auto it : fTofInfoCal )
  {
    if ( it ) delete it;
  };
  fTofInfoCal.clear();
 
  unsigned int ibSize    = test->getBTofSize();
  unsigned int ibcomSize = test->getBTofComSize();
  unsigned int ieSize    = test->getETofSize();
  unsigned int ietfSize  = test->getEtfSize();
  ietfSize               = ietfSize / 12;
  unsigned int iinfoSize = test->getTofInfoSize();
 
  for ( int ib = 0; ib < ibSize; ib++ )
  {
    EstBTofCal* btof = new EstBTofCal;
    for ( int j = 0; j < EstNBarPar; j++ )
    {
      btof->setP1( j, test->getBTofPleft( ib, j ) );
      btof->setP2( j, test->getBTofPright( ib, j ) );
    }
    fBTofCal.push_back( btof );
  }
 
  for ( int ie = 0; ie < ieSize; ie++ )
  {
    EstETofCal* etof = new EstETofCal;
    for ( int i = 0; i < EstNEndPar; i++ ) { etof->setP( i, test->getETofP( ie, i ) ); }
    fETofCal.push_back( etof );
  }
 
  for ( int ietf = 0; ietf < ietfSize; ietf++ )
  {
    for ( int jetf = 0; jetf < 12; jetf++ )
    {
      EstEtfCal* etf = new EstEtfCal;
      for ( int k = 0; k < EstNEtfPar; k++ )
      {
        etf->setP( k, test->getEtfPcombine( ietf, jetf, k ) );
        etf->setP1( k, test->getEtfPleft( ietf, jetf, k ) );
        etf->setP2( k, test->getEtfPright( ietf, jetf, k ) );
      }
      fEtfCal.push_back( etf );
    }
  }
 
  for ( int ibcom = 0; ibcom < ibcomSize; ibcom++ )
  {
    EstBTofCommonCal* bcomtof = new EstBTofCommonCal;
    for ( int i = 0; i < EstNBarOffset; i++ )
    { bcomtof->setOffset( i, test->getBTofOffset( ibcom, i ) ); }
    fBTofCommonCal.push_back( bcomtof );
  }
 
  for ( int iinfo = 0; iinfo < iinfoSize; iinfo++ )
  {
    EstTofInfoCal* tofinfo = new EstTofInfoCal;
    tofinfo->setRunFrom( test->getRunFrom( iinfo ) );
    tofinfo->setRunTo( test->getRunTo( iinfo ) );
    tofinfo->setEventFrom( test->getEventFrom( iinfo ) );
    tofinfo->setEventTo( test->getEventTo( iinfo ) );
    fTofInfoCal.push_back( tofinfo );
  }
 
  return StatusCode::SUCCESS;
}
 
StatusCode EstTofCaliSvc::chooseConstants( int run, int event ) {
 
  MsgStream log( msgSvc(), name() );
 
  m_sequence  = 0;
  bool filled = false;
 
  if ( fTofInfoCal[0]->getRunFrom() == -1 )
  {
    if ( fTofInfoCal.size() == 1 )
    {
      log << MSG::INFO
          << "EstTofCaliSvc::chooseConstants() -- Run From is equal to -1! The ONLY TOF "
             "calibration constants are used!"
          << endmsg;
      return StatusCode::SUCCESS;
    }
    else
    {
      log << MSG::ERROR
          << "EstTofCaliSvc::chooseConstants() -- Run From is equal to -1! The NUMBER of "
             "calibration constants are NOT "
             "equal to 1!"
          << endmsg;
      return StatusCode::FAILURE;
    }
  }
 
  for ( int i = 0; i < fTofInfoCal.size(); i++ )
  {
    auto element = fTofInfoCal[i];
    if ( element->getRunTo() != -1 && element->getRunTo() < element->getRunFrom() )
    {
      log << MSG::ERROR << "EstTofCaliSvc::chooseConstants() -- The " << i
          << "th calibration constatns is ABNORMAL! Run From is LARGER than RUN To!" << endmsg;
      return StatusCode::FAILURE;
    }
 
    if ( element->getRunFrom() == element->getRunTo() && element->getEventFrom() != -1 &&
         element->getEventTo() != -1 && element->getEventFrom() > element->getEventTo() )
    {
      log << MSG::ERROR << "EstTofCaliSvc::chooseConstants() -- The " << i
          << "th calibration constatns is ABNORMAL! Event From is LARGER than Event To!"
          << endmsg;
      return StatusCode::FAILURE;
    }
  }
 
  unsigned int inumber = 0;
  for ( ; inumber < fTofInfoCal.size(); inumber++ )
  {
    auto element   = fTofInfoCal[inumber];
    int  runFrom   = element->getRunFrom();
    int  runTo     = element->getRunTo();
    int  eventFrom = element->getEventFrom();
    int  eventTo   = element->getEventTo();
    if ( ( run == runFrom ) && ( ( eventFrom == -1 ) || ( event >= eventFrom ) ) )
    {
      if ( ( run < runTo ) ||
           ( ( run == runTo ) && ( ( eventTo == -1 ) || ( event <= eventTo ) ) ) )
      {
        filled = true;
        break;
      }
    }
    if ( run > runFrom )
    {
      if ( ( run < runTo ) ||
           ( ( run == runTo ) && ( ( eventTo == -1 ) || ( event <= eventTo ) ) ) )
      {
        filled = true;
        break;
      }
    }
  }
 
  if ( filled ) { m_sequence = inumber; }
  else
  {
    log << MSG::ERROR << "EstTofCaliSvc::chooseConstants() -- The event with run number "
        << run << " and event number " << event
        << " is NOT suitable for this group of calibration constants" << endmsg;
    return StatusCode::FAILURE;
  }
 
  return StatusCode::SUCCESS;
}
 
const double EstTofCaliSvc::BTCorr1( double ADC, double z, unsigned id ) {
  MsgStream log( msgSvc(), name() );
 
  double p1[EstNBarPar];
  for ( int i = 0; i < static_cast<int>( EstNBarPar ); i++ )
  { p1[i] = fBTofCal[176 * m_sequence + id]->getP1( i ); }
 
  log << MSG::DEBUG << "BTCorr1   id =" << id << "  1=" << p1[0] << " 2=" << p1[1]
      << " 3=" << p1[2] << " 4=" << p1[3] << " 5=" << p1[4] << " 6=" << p1[5] << " 7=" << p1[6]
      << " 8=" << p1[7] << " 9=" << p1[8] << endmsg;
 
  double tcorr1 = p1[0] + p1[1] / TMath::Sqrt( ADC ) + p1[2] * z / TMath::Sqrt( ADC ) +
                  p1[3] / ADC + p1[4] * z + p1[5] * z * z + p1[6] * z * z * z;
 
  return tcorr1;
}
 
const double EstTofCaliSvc::BTCorr2( double ADC, double z, unsigned id ) {
  MsgStream log( msgSvc(), name() );
 
  double p2[EstNBarPar];
  for ( int i = 0; i < static_cast<int>( EstNBarPar ); i++ )
  { p2[i] = fBTofCal[176 * m_sequence + id]->getP2( i ); }
 
  log << MSG::DEBUG << "BTCorr2   id =" << id << "  1=" << p2[0] << " 2=" << p2[1]
      << " 3=" << p2[2] << " 4=" << p2[3] << " 5=" << p2[4] << " 6=" << p2[5] << " 7=" << p2[6]
      << " 8=" << p2[7] << " 9=" << p2[8] << endmsg;
 
  double tcorr2 = p2[0] + p2[1] / TMath::Sqrt( ADC ) + p2[2] * z / TMath::Sqrt( ADC ) +
                  p2[3] / ADC + p2[4] * z + p2[5] * z * z + p2[6] * z * z * z;
 
  return tcorr2;
}
 
const double EstTofCaliSvc::TOffset() {
  double toffset = fBTofCommonCal[m_sequence]->getOffset( 0 );
  return toffset;
}
 
const double EstTofCaliSvc::BTime1( double ADC, double TDC, double z, unsigned id ) {
  MsgStream log( msgSvc(), name() );
 
  double tcorr1 = BTCorr1( ADC, z, id );
 
  double tcorr3 = fBTofCommonCal[m_sequence]->getOffset( 0 );
 
  log << MSG::DEBUG << "BTime1   id =" << id << "  tcorr=" << tcorr1
      << "  total offset=" << tcorr3 << endmsg;
 
  double time = TDC - tcorr1 - tcorr3;
 
  if ( time < 0. )
  {
    log << MSG::WARNING << "TofCaliSvc::BTime1() -- A minus time is given : " << time
        << " [Input] ADC = " << ADC << " TDC = " << TDC << " z = " << z << " ID = " << id
        << endmsg;
  }
  return time;
}
 
const double EstTofCaliSvc::BTime2( double ADC, double TDC, double z, unsigned id ) {
  MsgStream log( msgSvc(), name() );
 
  double tcorr1 = BTCorr2( ADC, z, id );
 
  double tcorr3 = fBTofCommonCal[m_sequence]->getOffset( 0 );
 
  log << MSG::DEBUG << "BTime2   id =" << id << "  tcorr=" << tcorr1
      << "  total offset=" << tcorr3 << endmsg;
 
  double time = TDC - tcorr1 - tcorr3;
 
  if ( time < 0. )
  {
    log << MSG::WARNING << "TofCaliSvc::BTime2() -- A minus time is given : " << time
        << " [Input] ADC = " << ADC << " TDC = " << TDC << " z = " << z << " ID = " << id
        << endmsg;
  }
  return time;
}
 
const double EstTofCaliSvc::ETime( double ADC, double TDC, double rHit, unsigned id ) {
  MsgStream log( msgSvc(), name() );
 
  double p[EstNEndPar];
  for ( int i = 0; i < static_cast<int>( EstNEndPar ); i++ )
  { p[i] = fETofCal[96 * m_sequence + id]->getP( i ); }
 
  log << MSG::DEBUG << "ETime   id =" << id << "  1=" << p[0] << " 2=" << p[1] << " 3=" << p[2]
      << " 4=" << p[3] << " 5=" << p[4] << " 6=" << p[5] << " 7=" << p[6] << endmsg;
 
  double time = TDC - ( p[0] + p[1] / TMath::Sqrt( ADC ) + p[2] / ADC + p[3] * ADC +
                        p[4] * rHit + p[5] * rHit * rHit + p[6] * rHit * rHit * rHit );
 
  if ( time < 0. )
  {
    log << MSG::WARNING << "TofCaliSvc::ETime() -- A minus time is given : " << time
        << " [Input] ADC = " << ADC << " TDC = " << TDC << " r = " << rHit << " ID = " << id
        << endmsg;
  }
 
  return time;
}
 
const double EstTofCaliSvc::EtfTime( double ADC1, double ADC2, double TDC1, double TDC2,
                                     unsigned id, unsigned int strip ) {
  MsgStream log( msgSvc(), name() );
 
  double q = ( ADC1 + ADC2 ) / 2.0;
  double t = ( TDC1 + TDC2 ) / 2.0;
 
  double p[EstNEtfPar];
  for ( int i = 0; i < static_cast<int>( EstNEtfPar ); i++ )
  { p[i] = fEtfCal[72 * 12 * m_sequence + 12 * id + strip]->getP( i ); }
 
  double tcorr = p[0] + p[1] / TMath::Sqrt( q ) + p[2] / q + p[3] * q + p[4] * q * q +
                 p[5] * q * q * q + p[6] * q * q * q * q;
  double time = t - tcorr;
 
  log << MSG::DEBUG << "EtfTime   module =" << id << " strip=" << strip << "  1=" << p[0]
      << " 2=" << p[1] << " 3=" << p[2] << " 4=" << p[3] << " 5=" << p[4] << " 6=" << p[5]
      << " 7=" << p[6] << " 8=" << p[7] << " 9=" << p[8] << " tcorr=" << tcorr
      << " time=" << time << endmsg;
  if ( time < 0. )
  {
    log << MSG::WARNING << "TofCaliSvc::EtfTime() -- A minus time is given : " << time
        << " [Input] ADC1 = " << ADC1 << " ADC2 = " << ADC2 << " TDC1 = " << TDC1
        << " TDC2 = " << TDC2 << " ID = " << id << " Strip = " << strip << endmsg;
  }
 
  return time;
}
 
const double EstTofCaliSvc::EtfTime1( double ADC, double TDC, double z, unsigned id,
                                      unsigned int strip ) {
  MsgStream log( msgSvc(), name() );
 
  double p1[EstNEtfPar];
  for ( int i = 0; i < static_cast<int>( EstNEtfPar ); i++ )
  { p1[i] = fEtfCal[72 * 12 * m_sequence + 12 * id + strip]->getP1( i ); }
 
  double tcorr1 = p1[0] + p1[1] / TMath::Sqrt( ADC ) + p1[2] * ADC + p1[3] * ADC * ADC +
                  p1[4] * ADC * ADC * ADC + p1[5] * ADC * ADC * ADC * ADC +
                  p1[6] * z / TMath::Sqrt( ADC ) + p1[7] / ADC + p1[8] * z / ADC +
                  p1[9] * z * ADC + p1[10] * z + p1[11] * z * z + p1[12] * z * z * z;
  double time1 = TDC - tcorr1;
 
  log << MSG::DEBUG << "EtfTime1   module =" << id << " strip=" << strip << "  1=" << p1[0]
      << " 2=" << p1[1] << " 3=" << p1[2] << " 4=" << p1[3] << " 5=" << p1[4] << " 6=" << p1[5]
      << " 7=" << p1[6] << " 8=" << p1[7] << " 9=" << p1[8] << " tcorr=" << tcorr1
      << " time1=" << time1 << endmsg;
  if ( time1 < 0. )
  {
    log << MSG::WARNING << "TofCaliSvc::EtfTime1() -- A minus time is given : " << time
        << " [Input] ADC = " << ADC << " TDC = " << TDC << " z = " << z << " ID = " << id
        << " Strip = " << strip << endmsg;
  }
 
  return time1;
}
 
const double EstTofCaliSvc::EtfTime2( double ADC, double TDC, double z, unsigned id,
                                      unsigned int strip ) {
  MsgStream log( msgSvc(), name() );
 
  double p2[EstNEtfPar];
  for ( int i = 0; i < static_cast<int>( EstNEtfPar ); i++ )
  { p2[i] = fEtfCal[72 * 12 * m_sequence + 12 * id + strip]->getP2( i ); }
 
  double tcorr2 = p2[0] + p2[1] / TMath::Sqrt( ADC ) + p2[2] * ADC + p2[3] * ADC * ADC +
                  p2[4] * ADC * ADC * ADC + p2[5] * ADC * ADC * ADC * ADC +
                  p2[6] * z / TMath::Sqrt( ADC ) + p2[7] / ADC + p2[8] * z / ADC +
                  p2[9] * z * ADC + p2[10] * z + p2[11] * z * z + p2[12] * z * z * z;
  double time2 = TDC - tcorr2;
 
  log << MSG::DEBUG << "EtfTime2   module =" << id << " strip=" << strip << "  1=" << p2[0]
      << " 2=" << p2[1] << " 3=" << p2[2] << " 4=" << p2[3] << " 5=" << p2[4] << " 6=" << p2[5]
      << " 7=" << p2[6] << " 8=" << p2[7] << " 9=" << p2[8] << " tcorr=" << tcorr2
      << " time2=" << time2 << endmsg;
  if ( time2 < 0. )
  {
    log << MSG::WARNING << "TofCaliSvc::EtfTime2() -- A minus time is given : " << time
        << " [Input] ADC = " << ADC << " TDC = " << TDC << " z = " << z << " ID = " << id
        << " Strip = " << strip << endmsg;
  }
 
  return time2;
}
 
const double EstTofCaliSvc::EtfTime( double TDC1, double TDC2, unsigned id,
                                     unsigned int strip ) {
  MsgStream log( msgSvc(), name() );
 
  double time = -99;
  double t, tcorr;
 
  if ( TDC1 > 0 && TDC2 > 0 )
  {
    t     = ( TDC1 + TDC2 ) / 2.0;
    tcorr = fEtfCal[72 * 12 * m_sequence + 12 * id + strip]->getP( 0 );
    time  = t - tcorr;
    log << MSG::DEBUG << "EtfTime   module =" << id << " strip=" << strip << " t1=" << TDC1
        << " t2=" << TDC2 << "  combine const=" << tcorr << " time=" << time << endmsg;
  }
  else if ( TDC1 > 0 && TDC2 < 0 )
  {
    t     = TDC1;
    tcorr = fEtfCal[72 * 12 * m_sequence + 12 * id + strip]->getP1( 0 );
    time  = t - tcorr;
    log << MSG::DEBUG << "EtfTime   module =" << id << " strip=" << strip << " t1=" << TDC1
        << "  left const=" << tcorr << " time=" << time << endmsg;
  }
  else if ( TDC1 < 0 && TDC2 > 0 )
  {
    t     = TDC2;
    tcorr = fEtfCal[72 * 12 * m_sequence + 12 * id + strip]->getP2( 0 );
    time  = t - tcorr;
    log << MSG::DEBUG << "EtfTime   module =" << id << " strip=" << strip << " t2=" << TDC2
        << "  right const=" << tcorr << " time=" << time << endmsg;
  }
 
  if ( time < 0. )
  {
    log << MSG::WARNING << "TofCaliSvc::EtfTime() -- A minus time is given : " << time
        << " [Input]  TDC1 = " << TDC1 << " TDC2 = " << TDC2 << " ID = " << id
        << " Strip = " << strip << endmsg;
  }
 
  return time;
}
 
EstBTofCal* EstTofCaliSvc::BTof( unsigned id ) const {
  return fBTofCal[176 * m_sequence + id];
}
 
EstETofCal* EstTofCaliSvc::ETof( unsigned id ) const { return fETofCal[96 * m_sequence + id]; }
 
EstEtfCal* EstTofCaliSvc::EtfTof( unsigned id ) const {
  return fEtfCal[72 * 12 * m_sequence + id];
}
 
EstBTofCommonCal* EstTofCaliSvc::BTofCommon() const { return fBTofCommonCal[m_sequence]; }
 
EstTofInfoCal* EstTofCaliSvc::TofInfo() const { return fTofInfoCal[m_sequence]; }
 
void EstTofCaliSvc::Dump() {
  std::cout << "Now We can get the TOF Calibtration Service" << std::endl;
  std::cout << "Barrel TOF Counter Number = " << fBTofCal.size() << std::endl;
  std::cout << "Endcap TOF Counter Number = " << fETofCal.size() << std::endl;
 
  cout << "The O  Barrel TOF PMT1 p1[0]  is"
       << " " << BTof( 0 )->getP1( 0 ) << endl;
  cout << "The 0  Barrel TOF PMT1 p1[1]  is"
       << " " << BTof( 0 )->getP1( 1 ) << endl;
}
 
void EstTofCaliSvc::handle( const Incident& inc ) {
  MsgStream log( msgSvc(), name() );
  log << MSG::DEBUG << "handle: " << inc.type() << endmsg;
 
  if ( inc.type() == "NewRun" )
  {
    log << MSG::DEBUG << "New Run" << endmsg;
    StatusCode sc = FillfromDatabase();
    if ( sc.isSuccess() ) { m_hasbeeninitialized = true; }
  }
}