T0MdcCalib.cxx

Go to the documentation of this file.

#include "MdcCalibAlg/T0MdcCalib.h"
 
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/IMessageSvc.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/StatusCode.h"
 
#include <cstring>
#include <iomanip>
#include <iostream>
 
T0MdcCalib::T0MdcCalib() {
  for ( int lay = 0; lay < MdcCalNLayer; lay++ )
  {
    if ( lay < 8 )
    {
      m_docaMin[lay] = 1.2; // mm
      m_docaMax[lay] = 4.8; // mm
    }
    else
    {
      m_docaMin[lay] = 1.6; // mm
      m_docaMax[lay] = 6.4; // mm
    }
  }
}
 
T0MdcCalib::~T0MdcCalib() {}
 
void T0MdcCalib::clear() {
  for ( int i = 0; i < MdcCalTotCell; i++ )
  {
    delete m_hleft[i];
    delete m_hright[i];
  }
  delete m_hLrResiSum;
  delete m_hLrResiSub;
  delete m_fdT0;
  delete m_fdResiWire;
 
  MdcCalib::clear();
}
 
void T0MdcCalib::initialize( TObjArray* hlist, IMdcGeomSvc* mdcGeomSvc,
                             IMdcCalibFunSvc* mdcFunSvc, IMdcUtilitySvc* mdcUtilitySvc ) {
  IMessageSvc* msgSvc;
  Gaudi::svcLocator()->service( "MessageSvc", msgSvc );
  MsgStream log( msgSvc, "T0MdcCalib" );
  log << MSG::INFO << "T0MdcCalib::initialize()" << endmsg;
 
  m_hlist         = hlist;
  m_mdcGeomSvc    = mdcGeomSvc;
  m_mdcFunSvc     = mdcFunSvc;
  m_mdcUtilitySvc = mdcUtilitySvc;
 
  m_vdr = 0.03;
 
  MdcCalib::initialize( m_hlist, m_mdcGeomSvc, m_mdcFunSvc, m_mdcUtilitySvc );
 
  int  i;
  int  nwire;
  int  lay;
  int  cel;
  char hname[200];
 
  m_fdT0 = new TFolder( "fdT0", "fdT0" );
  m_hlist->Add( m_fdT0 );
 
  m_fdResiWire = new TFolder( "ResiWire", "ResiWire" );
  m_hlist->Add( m_fdResiWire );
 
  nwire = m_mdcGeomSvc->getWireSize();
  for ( i = 0; i < nwire; i++ )
  {
    lay = m_mdcGeomSvc->Wire( i )->Layer();
    cel = m_mdcGeomSvc->Wire( i )->Cell();
 
    sprintf( hname, "Resi%04d_Lay%02d_Cell%03d_L", i, lay, cel );
    m_hleft[i] = new TH1F( hname, "", 400, -2.0, 2.0 );
    m_fdT0->Add( m_hleft[i] );
 
    sprintf( hname, "Resi%04d_Lay%02d_Cell%03d_R", i, lay, cel );
    m_hright[i] = new TH1F( hname, "", 400, -2.0, 2.0 );
    m_fdT0->Add( m_hright[i] );
  }
 
  m_hLrResiSum = new TH1F( "LrResiSum", "", 200, -0.5, 0.5 );
  m_fdResiWire->Add( m_hLrResiSum );
 
  m_hLrResiSub = new TH1F( "LrResiSub", "", 200, -0.5, 0.5 );
  m_fdResiWire->Add( m_hLrResiSub );
}
 
int T0MdcCalib::fillHist( MdcCalEvent* event ) {
  IMessageSvc* msgSvc;
  Gaudi::svcLocator()->service( "MessageSvc", msgSvc );
  MsgStream log( msgSvc, "T0MdcCalib" );
  log << MSG::DEBUG << "T0MdcCalib::fillHist()" << endmsg;
 
  MdcCalib::fillHist( event );
 
  // get EsTimeCol
  bool esCutFg = event->getEsCutFlag();
  if ( !esCutFg ) return -1;
 
  int i;
  int k;
  int ntrk;
  int nhit;
  int stat;
 
  int    lay;
  int    cel;
  int    wir;
  int    lr;
  double dmeas;
  double doca;
  double resi;
 
  bool fgHitLay[MdcCalNLayer];
 
  MdcCalRecTrk* rectrk;
  MdcCalRecHit* rechit;
 
  ntrk = event->getNTrk();
  if ( ( ntrk < m_param.nTrkCut[0] ) || ( ntrk > m_param.nTrkCut[1] ) ) return -1;
 
  for ( i = 0; i < ntrk; i++ )
  {
    rectrk = event->getRecTrk( i );
    nhit   = rectrk->getNHits();
 
    // dr cut
    double dr = rectrk->getDr();
    if ( fabs( dr ) > m_param.drCut ) continue;
 
    // dz cut
    double dz = rectrk->getDz();
    if ( fabs( dz ) > m_param.dzCut ) continue;
 
    // momentum cut
    double p = rectrk->getP();
    if ( ( fabs( p ) < m_param.pCut[0] ) || ( fabs( p ) > m_param.pCut[1] ) ) continue;
 
    for ( lay = 0; lay < MdcCalNLayer; lay++ ) fgHitLay[lay] = false;
    for ( k = 0; k < nhit; k++ )
    {
      rechit        = rectrk->getRecHit( k );
      lay           = rechit->getLayid();
      fgHitLay[lay] = true;
    }
 
    int nhitlay = 0;
    for ( lay = 0; lay < MdcCalNLayer; lay++ )
      if ( fgHitLay[lay] ) nhitlay++;
    if ( nhitlay < m_param.nHitLayCut ) continue;
 
    bool fgNoise = rectrk->getFgNoiseRatio();
    if ( m_param.noiseCut && ( !fgNoise ) ) continue;
 
    //    log << MSG::DEBUG << "nhits: " << nhit << endmsg;
    for ( k = 0; k < nhit; k++ )
    {
      rechit = rectrk->getRecHit( k );
      lay    = rechit->getLayid();
      cel    = rechit->getCellid();
      wir    = m_mdcGeomSvc->Wire( lay, cel )->Id();
      lr     = rechit->getLR();
      doca   = rechit->getDocaExc();
      dmeas  = rechit->getDmeas();
      resi   = rechit->getResiExc();
      stat   = rechit->getStat();
 
      if ( 1 != stat ) continue;
 
      if ( ( fabs( doca ) < m_docaMin[lay] ) || ( fabs( doca ) > m_docaMax[lay] ) ||
           ( fabs( resi ) > m_param.resiCut[lay] ) )
      { continue; }
 
      if ( 0 == lay )
      {
        if ( !fgHitLay[1] ) continue;
      }
      else if ( 42 == lay )
      {
        if ( !fgHitLay[41] ) continue;
      }
      else
      {
        if ( ( !fgHitLay[lay - 1] ) && ( !fgHitLay[lay + 1] ) ) continue;
      }
 
      if ( wir < ( m_mdcGeomSvc->getWireSize() ) )
      {
        if ( 0 == lr ) { m_hleft[wir]->Fill( resi ); }
        else if ( 1 == lr ) { m_hright[wir]->Fill( resi ); }
      }
      else { std::cout << "wireid: " << wir << std::endl; }
    }
  }
  return 1;
}
 
void T0MdcCalib::printCut() const { MdcCalib::printCut(); }
 
int T0MdcCalib::updateConst( MdcCalibConst* calconst ) {
  IMessageSvc* msgSvc;
  Gaudi::svcLocator()->service( "MessageSvc", msgSvc );
  MsgStream log( msgSvc, "T0MdcCalib" );
  log << MSG::INFO << "T0MdcCalib::updateConst()" << endmsg;
 
  MdcCalib::updateConst( calconst );
 
  int    i;
  int    lay;
  int    nwire;
  double t0;
  double delt0;
  double resiLrSum;
  double resiLrSub;
 
  Stat_t entry_l;
  double mean_l;
  //      double gmean_l;  // mean of fit with the gaussian distribution
  //      double sigma_l;
  //      double rms_l;
 
  Stat_t entry_r;
  double mean_r;
  //      double gmean_r;
  //      double sigma_r;
  //      double rms_r;
 
  nwire = m_mdcGeomSvc->getWireSize();
  std::cout << "totwire: " << nwire << std::endl;
  for ( i = 0; i < nwire; i++ )
  {
    lay = m_mdcGeomSvc->Wire( i )->Layer();
 
    if ( 1 == m_param.fgCalib[lay] )
    {
      entry_l = m_hleft[i]->GetEntries();
      mean_l  = m_hleft[i]->GetMean();
      //           m_hleft[i] -> Fit("gaus", "Q");
      //           gmean_l = m_hleft[i] -> GetFunction("gaus") -> GetParameter(1);
 
      entry_r = m_hright[i]->GetEntries();
      mean_r  = m_hright[i]->GetMean();
      //           m_hright[i] -> Fit("gaus", "Q");
      //           gmean_r = m_hright[i] -> GetFunction("gaus") -> GetParameter(1);
 
      delt0 = 0.5 * ( mean_l + mean_r ) / m_vdr;
    }
    else { delt0 = 0.0; }
    if ( 0 == m_param.fgWireCal[i] ) delt0 = 0.0;
 
    resiLrSum = 0.5 * ( mean_l + mean_r );
    resiLrSub = 0.5 * ( mean_l - mean_r );
    m_hLrResiSum->Fill( resiLrSum );
    m_hLrResiSub->Fill( resiLrSub );
 
    t0 = calconst->getT0( i );
    t0 += delt0;
 
    calconst->resetT0( i, t0 );
    calconst->resetDelT0( i, delt0 );
  }
  return 1;
}