XtInteMdcCalib.cxx

Go to the documentation of this file.

#include "MdcCalibAlg/XtInteMdcCalib.h"
 
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/IMessageSvc.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/StatusCode.h"
 
#include <cmath>
#include <cstring>
#include <fstream>
#include <iomanip>
#include <iostream>
 
#include "TF1.h"
 
using namespace std;
 
XtInteMdcCalib::XtInteMdcCalib() {
  m_fgIni       = false;
  m_nMaxGrPoint = 5000;
  for ( int lay = 0; lay < MdcCalNLayer; lay++ )
  {
    m_tbinWid[lay][0] = 5.0;
    m_tbinWid[lay][1] = 10.0;
    m_tbinWid[lay][2] = 20.0;
 
    m_tbinLim[lay][0] = -10.0;
    m_tbinLim[lay][1] = 30.0;
    if ( lay < 8 ) m_tbinLim[lay][2] = 210.0;
    else m_tbinLim[lay][2] = 350.0;
    m_tbinLim[lay][3] = 990.0;
  }
}
 
XtInteMdcCalib::~XtInteMdcCalib() {}
 
void XtInteMdcCalib::clear() {
  for ( int lay = 0; lay < MdcCalNLayer; lay++ )
  {
    for ( int iEntr = 0; iEntr < NENTR; iEntr++ )
    {
      for ( int lr = 0; lr < 2; lr++ )
      {
        delete m_pfNear[lay][iEntr][lr];
        delete m_pfMid[lay][iEntr][lr];
        delete m_pfFar[lay][iEntr][lr];
        delete m_grXt[lay][iEntr][lr];
      }
    }
  }
  delete m_fdPf;
  MdcCalib::clear();
}
 
void XtInteMdcCalib::initialize( TObjArray* hlist, IMdcGeomSvc* mdcGeomSvc,
                                 IMdcCalibFunSvc* mdcFunSvc, IMdcUtilitySvc* mdcUtilitySvc ) {
  IMessageSvc* msgSvc;
  Gaudi::svcLocator()->service( "MessageSvc", msgSvc );
  MsgStream log( msgSvc, "XtInteMdcCalib" );
  log << MSG::INFO << "XtInteMdcCalib::initialize()" << endmsg;
 
  m_hlist         = hlist;
  m_mdcGeomSvc    = mdcGeomSvc;
  m_mdcFunSvc     = mdcFunSvc;
  m_mdcUtilitySvc = mdcUtilitySvc;
 
  MdcCalib::initialize( m_hlist, m_mdcGeomSvc, m_mdcFunSvc, m_mdcUtilitySvc );
 
  m_fdPf = new TFolder( "fdProfile", "fdProfile" );
  m_hlist->Add( m_fdPf );
 
  char hname[200];
  for ( int lay = 0; lay < MdcCalNLayer; lay++ )
  {
    for ( int iEntr = 0; iEntr < NENTR; iEntr++ )
    {
      for ( int lr = 0; lr < 2; lr++ )
      {
        sprintf( hname, "xt%02d_%02d_%d_gr", lay, iEntr, lr );
        m_grXt[lay][iEntr][lr] = new TGraph();
        m_grXt[lay][iEntr][lr]->SetName( hname );
        m_fdPf->Add( m_grXt[lay][iEntr][lr] );
 
        int xbinN =
            (int)( ( m_tbinLim[lay][1] - m_tbinLim[lay][0] ) / m_tbinWid[lay][0] + 0.5 );
        sprintf( hname, "xt%02d_%02d_%d_near", lay, iEntr, lr );
        m_pfNear[lay][iEntr][lr] =
            new TProfile( hname, hname, xbinN, m_tbinLim[lay][0], m_tbinLim[lay][1] );
        m_fdPf->Add( m_pfNear[lay][iEntr][lr] );
 
        int xbinM =
            (int)( ( m_tbinLim[lay][2] - m_tbinLim[lay][1] ) / m_tbinWid[lay][1] + 0.5 );
        sprintf( hname, "xt%02d_%02d_%d_mid", lay, iEntr, lr );
        m_pfMid[lay][iEntr][lr] =
            new TProfile( hname, hname, xbinM, m_tbinLim[lay][1], m_tbinLim[lay][2] );
        m_fdPf->Add( m_pfMid[lay][iEntr][lr] );
 
        int xbinF =
            (int)( ( m_tbinLim[lay][3] - m_tbinLim[lay][2] ) / m_tbinWid[lay][2] + 0.5 );
        sprintf( hname, "xt%02d_%02d_%d_far", lay, iEntr, lr );
        m_pfFar[lay][iEntr][lr] =
            new TProfile( hname, hname, xbinF, m_tbinLim[lay][2], m_tbinLim[lay][3] );
        m_fdPf->Add( m_pfFar[lay][iEntr][lr] );
 
        //          if((0==iEntr)&&(0==lr))
        //          cout << setw(5) << lay
        //           << setw(5) << xbinN << setw(10) << m_tbinLim[lay][0] <<
        // setw(10) << m_tbinLim[lay][1]
        //           << setw(5) << xbinM << setw(10) << m_tbinLim[lay][1] <<
        // setw(10) << m_tbinLim[lay][2]
        //           << setw(5) << xbinF << setw(10) << m_tbinLim[lay][2] <<
        // setw(10) << m_tbinLim[lay][3] << endl;
      }
    }
  }
}
 
int XtInteMdcCalib::fillHist( MdcCalEvent* event ) {
  IMessageSvc* msgSvc;
  Gaudi::svcLocator()->service( "MessageSvc", msgSvc );
  MsgStream log( msgSvc, "XtInteMdcCalib" );
  log << MSG::DEBUG << "XtInteMdcCalib::fillHist()" << endmsg;
 
  MdcCalib::fillHist( event );
 
  // get EsTimeCol
  bool esCutFg = event->getEsCutFlag();
  if ( !esCutFg ) return -1;
 
  int i;
  int k;
  int lay;
  int iLR;
  int iEntr;
 
  double dr;
  double dz;
  double doca;
  double tdr;
  double resi;
  double entrance;
 
  int  nhitlay;
  bool fgHitLay[MdcCalNLayer];
  bool fgTrk;
 
  if ( !m_fgIni )
  {
    for ( lay = 0; lay < MdcCalNLayer; lay++ )
    {
      if ( lay < 8 ) m_docaMax[lay] = m_param.maxDocaInner;
      else m_docaMax[lay] = m_param.maxDocaOuter;
    }
    m_fgIni = true;
  }
 
  MdcCalRecTrk* rectrk;
  MdcCalRecHit* rechit;
 
  int nhit;
  int ntrk = event->getNTrk();
  for ( i = 0; i < ntrk; i++ )
  {
    fgTrk  = true;
    rectrk = event->getRecTrk( i );
    nhit   = rectrk->getNHits();
 
    // dr cut
    dr = rectrk->getDr();
    if ( fabs( dr ) > m_param.drCut ) continue;
 
    // dz cut
    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;
 
    // select cosmic track with y<0
    double phi0   = rectrk->getPhi0();
    double phiTrk = phi0 + CLHEP::halfpi;
    if ( phiTrk > CLHEP::twopi ) phiTrk -= CLHEP::twopi;
    if ( m_param.cosmicDwTrk && ( phiTrk < CLHEP::pi ) ) continue;
    //    cout << "cosmicDwTrk " << m_param.cosmicDwTrk << setw(15) << phiTrk*180./3.14 <<
    // endl;
 
    for ( lay = 0; lay < MdcCalNLayer; lay++ ) { fgHitLay[lay] = false; }
 
    for ( k = 0; k < nhit; k++ )
    {
      rechit        = rectrk->getRecHit( k );
      lay           = rechit->getLayid();
      doca          = rechit->getDocaInc();
      resi          = rechit->getResiInc();
      fgHitLay[lay] = true;
 
      //           if( (fabs(doca) > m_docaMax[lay]) ||
      //           (fabs(resi) > m_param.resiCut[lay]) ){
      //            fgTrk = false;
      //           }
    }
    if ( !fgTrk ) continue;
 
    nhitlay = 0;
    for ( lay = 0; lay < MdcCalNLayer; lay++ )
    {
      if ( fgHitLay[lay] ) nhitlay++;
    }
    if ( nhitlay < m_param.nHitLayCut ) continue;
 
    for ( k = 0; k < nhit; k++ )
    {
      rechit   = rectrk->getRecHit( k );
      lay      = rechit->getLayid();
      doca     = rechit->getDocaInc();
      resi     = rechit->getResiInc();
      iLR      = rechit->getLR();
      entrance = rechit->getEntra();
      tdr      = rechit->getTdrift();
 
      if ( ( 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;
      }
 
      int nPM = 1;
      if ( m_param.fgCombinePlusMinue ) nPM = 2;
      for ( int iPM = 0; iPM < nPM; iPM++ )
      {
        if ( 0 == iPM ) iEntr = m_mdcFunSvc->getXtEntrIndex( entrance );
        else iEntr = m_mdcFunSvc->getXtEntrIndex( -1.0 * entrance );
 
        int npoint = m_grXt[lay][iEntr][iLR]->GetN();
        if ( npoint < m_nMaxGrPoint )
          m_grXt[lay][iEntr][iLR]->SetPoint( npoint, tdr, fabs( doca ) );
 
        if ( tdr < m_tbinLim[lay][0] ) continue;
        else if ( tdr < m_tbinLim[lay][1] )
          m_pfNear[lay][iEntr][iLR]->Fill( tdr, fabs( doca ) );
        else if ( tdr < m_tbinLim[lay][2] )
          m_pfMid[lay][iEntr][iLR]->Fill( tdr, fabs( doca ) );
        else if ( tdr < m_tbinLim[lay][3] )
          m_pfFar[lay][iEntr][iLR]->Fill( tdr, fabs( doca ) );
      }
    } // hit loop
  }   // track loop
  return 1;
}
 
void XtInteMdcCalib::printCut() const { MdcCalib::printCut(); }
 
int XtInteMdcCalib::updateConst( MdcCalibConst* calconst ) {
  IMessageSvc* msgSvc;
  Gaudi::svcLocator()->service( "MessageSvc", msgSvc );
  MsgStream log( msgSvc, "XtInteMdcCalib" );
  log << MSG::INFO << "XtInteMdcCalib::updateConst()" << endmsg;
 
  MdcCalib::updateConst( calconst );
}