GrXtMdcCalib.cxx

Go to the documentation of this file.

#include "MdcCalibAlg/GrXtMdcCalib.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;
 
double GrXtMdcCalib::DMAX;
double GrXtMdcCalib::TMAX;
 
GrXtMdcCalib::GrXtMdcCalib() {
  m_maxNhit = 5000;
  m_fgIni   = false;
}
 
GrXtMdcCalib::~GrXtMdcCalib() {}
 
void GrXtMdcCalib::clear() {
  cout << "~GrXtMdcCalib" << endl;
  for ( int lay = 0; lay < MdcCalNLayer; lay++ )
  {
    for ( int iEntr = 0; iEntr < MdcCalNENTRXT; iEntr++ )
    {
      for ( int iLR = 0; iLR < MdcCalLR; iLR++ ) { delete m_grxt[lay][iEntr][iLR]; }
    }
  }
  delete m_haxis;
  delete m_fdXt;
 
  MdcCalib::clear();
}
 
void GrXtMdcCalib::initialize( TObjArray* hlist, IMdcGeomSvc* mdcGeomSvc,
                               IMdcCalibFunSvc* mdcFunSvc, IMdcUtilitySvc* mdcUtilitySvc ) {
  IMessageSvc* msgSvc;
  Gaudi::svcLocator()->service( "MessageSvc", msgSvc );
  MsgStream log( msgSvc, "GrXtMdcCalib" );
  log << MSG::INFO << "GrXtMdcCalib::initialize()" << endmsg;
 
  m_hlist         = hlist;
  m_mdcGeomSvc    = mdcGeomSvc;
  m_mdcFunSvc     = mdcFunSvc;
  m_mdcUtilitySvc = mdcUtilitySvc;
 
  MdcCalib::initialize( m_hlist, m_mdcGeomSvc, m_mdcFunSvc, m_mdcUtilitySvc );
 
  int  lay;
  int  iLR;
  int  iEntr;
  char hname[200];
 
  m_fdXt = new TFolder( "fdXtGr", "fdXtGr" );
  m_hlist->Add( m_fdXt );
 
  m_haxis = new TH2F( "axis", "", 50, 0, 300, 50, 0, 9 );
  m_haxis->SetStats( 0 );
  m_fdXt->Add( m_haxis );
 
  for ( lay = 0; lay < MdcCalNLayer; lay++ )
  {
    for ( iEntr = 0; iEntr < MdcCalNENTRXT; iEntr++ )
    {
      for ( iLR = 0; iLR < MdcCalLR; iLR++ )
      {
        m_nhit[lay][iEntr][iLR] = 0;
 
        sprintf( hname, "grXt%02d_%02d_lr%01d", lay, iEntr, iLR );
        m_grxt[lay][iEntr][iLR] = new TGraph();
        m_grxt[lay][iEntr][iLR]->SetName( hname );
        m_grxt[lay][iEntr][iLR]->SetMarkerStyle( 10 );
        m_grxt[lay][iEntr][iLR]->SetLineColor( 10 );
        m_fdXt->Add( m_grxt[lay][iEntr][iLR] );
      }
    }
  }
}
 
int GrXtMdcCalib::fillHist( MdcCalEvent* event ) {
  IMessageSvc* msgSvc;
  Gaudi::svcLocator()->service( "MessageSvc", msgSvc );
  MsgStream log( msgSvc, "GrXtMdcCalib" );
  log << MSG::DEBUG << "GrXtMdcCalib::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;
 
    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;
      }
 
      iEntr = m_mdcFunSvc->getXtEntrIndex( entrance );
 
      if ( iLR < 2 )
      {
        if ( m_nhit[lay][iEntr][iLR] < m_maxNhit )
        {
          m_grxt[lay][iEntr][iLR]->SetPoint( m_nhit[lay][iEntr][iLR], tdr, fabs( doca ) );
          m_nhit[lay][iEntr][iLR]++;
        }
      }
 
      if ( m_nhit[lay][iEntr][2] < m_maxNhit )
      {
        m_grxt[lay][iEntr][2]->SetPoint( m_nhit[lay][iEntr][2], tdr, fabs( doca ) );
        m_nhit[lay][iEntr][2]++;
      }
    } // hit loop
  }   // track loop
  return 1;
}
 
void GrXtMdcCalib::printCut() const { MdcCalib::printCut(); }
 
int GrXtMdcCalib::updateConst( MdcCalibConst* calconst ) {
  IMessageSvc* msgSvc;
  Gaudi::svcLocator()->service( "MessageSvc", msgSvc );
  MsgStream log( msgSvc, "GrXtMdcCalib" );
  log << MSG::INFO << "GrXtMdcCalib::updateConst()" << endmsg;
 
  MdcCalib::updateConst( calconst );
 
  int    ord;
  double xtpar[MdcCalNLayer][MdcCalNENTRXT][MdcCalLR][8];
  TF1*   fxtDr = new TF1( "fxtDr", xtfun, 0, 300, 6 );
  TF1*   fxtEd = new TF1( "fxtEd", xtedge, 150, 500, 1 );
  if ( 1 == m_param.fixXtC0 ) fxtDr->FixParameter( 0, 0 );
 
  for ( int lay = 0; lay < MdcCalNLayer; lay++ )
  {
    for ( int iEntr = 0; iEntr < MdcCalNENTRXT; iEntr++ )
    {
      for ( int lr = 0; lr < MdcCalLR; lr++ )
      {
        m_fgFit[lay][iEntr][lr] = false;
        if ( 0 == m_param.fgCalib[lay] ) continue;
 
        if ( m_nhit[lay][iEntr][lr] > 1000 )
        {
          TMAX = calconst->getXtpar( lay, iEntr, lr, 6 );
 
          m_grxt[lay][iEntr][lr]->Fit( "fxtDr", "Q+", "", 0, TMAX );
 
          for ( ord = 0; ord < 6; ord++ )
          { xtpar[lay][iEntr][lr][ord] = fxtDr->GetParameter( ord ); }
          xtpar[lay][iEntr][lr][6] = TMAX;
 
          DMAX = 0.0;
          for ( ord = 0; ord < 6; ord++ )
            DMAX += xtpar[lay][iEntr][lr][ord] * pow( TMAX, ord );
 
          m_grxt[lay][iEntr][lr]->Fit( "fxtEd", "Q+", "", TMAX, TMAX + 300 );
          xtpar[lay][iEntr][lr][7] = fxtEd->GetParameter( 0 );
          if ( xtpar[lay][iEntr][lr][7] < 0.0 ) xtpar[lay][iEntr][lr][7] = 0.0;
          m_fgFit[lay][iEntr][lr] = true;
 
          //             for(ord=0; ord<8; ord++){
          //                  calconst -> resetXtpar(lay, iEntr, lr, ord, xtpar[ord]);
          //             }
        }
 
      } // end of lr loop
    }   // end of entrance angle loop
  }     // end of layer loop
 
  ofstream fxtlog( "xtlog" );
  for ( int lay = 0; lay < MdcCalNLayer; lay++ )
  {
    for ( int iEntr = 0; iEntr < MdcCalNENTRXT; iEntr++ )
    {
      for ( int lr = 0; lr < MdcCalLR; lr++ )
      {
        fxtlog << setw( 3 ) << lay << setw( 3 ) << iEntr << setw( 3 ) << lr;
 
        int fgUpdate = -1;
        if ( m_fgFit[lay][iEntr][lr] )
        {
          fgUpdate = 1;
          for ( ord = 0; ord < 8; ord++ )
            calconst->resetXtpar( lay, iEntr, lr, ord, xtpar[lay][iEntr][lr][ord] );
        }
        else
        {
          int iEntrNew = findXtEntr( lay, iEntr, lr );
          if ( -1 != iEntrNew )
          {
            fgUpdate = 2;
            for ( ord = 0; ord < 8; ord++ )
            { calconst->resetXtpar( lay, iEntr, lr, ord, xtpar[lay][iEntrNew][lr][ord] ); }
          }
        }
        fxtlog << setw( 3 ) << fgUpdate;
        for ( ord = 0; ord < 8; ord++ )
        {
          double par = calconst->getXtpar( lay, iEntr, lr, ord );
          fxtlog << setw( 14 ) << par;
        }
        fxtlog << endl;
      }
    }
  }
  fxtlog.close();
 
  cout << "Xt update finished. File xtlog was written." << endl;
  delete fxtDr;
  delete fxtEd;
}
 
Double_t GrXtMdcCalib::xtfun( Double_t* x, Double_t* par ) {
  Double_t val = 0.0;
  for ( Int_t ord = 0; ord < 6; ord++ ) { val += par[ord] * pow( x[0], ord ); }
  return val;
}
 
Double_t GrXtMdcCalib::xtedge( Double_t* x, Double_t* par ) {
  double val = DMAX + ( x[0] - TMAX ) * par[0];
  return val;
}
 
int GrXtMdcCalib::findXtEntr( int lay, int iEntr, int lr ) const {
  int id0    = 8;
  int id1    = 9;
  int idmax  = 17;
  int entrId = -1;
  if ( iEntr <= id0 )
  {
    int id = -1;
    for ( int i = iEntr; i <= id0; i++ )
    {
      if ( m_fgFit[lay][i][lr] )
      {
        id = i;
        break;
      }
    }
    if ( -1 != id ) entrId = id;
    else
    {
      for ( int i = iEntr; i >= 0; i-- )
      {
        if ( m_fgFit[lay][i][lr] )
        {
          id = i;
          break;
        }
      }
      entrId = id;
    }
  }
  else
  {
    int id = -1;
    for ( int i = iEntr; i >= id1; i-- )
    {
      if ( m_fgFit[lay][i][lr] )
      {
        id = i;
        break;
      }
    }
    if ( -1 != id ) entrId = id;
    else
    {
      for ( int i = iEntr; i < idmax; i++ )
      {
        if ( m_fgFit[lay][i][lr] )
        {
          id = i;
          break;
        }
      }
      entrId = id;
    }
  }
  if ( -1 == entrId )
  {
    cout << "find EntrId error "
         << "layer " << lay << "  iEntr " << iEntr << "  lr " << lr << endl;
  }
 
  return entrId;
}