PreXtMdcCalib.cxx

Go to the documentation of this file.

#include "MdcCalibAlg/PreXtMdcCalib.h"
 
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/IMessageSvc.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/StatusCode.h"
 
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/PropertyMgr.h"
#include "GaudiKernel/SmartDataPtr.h"
 
#include "EvTimeEvent/RecEsTime.h"
#include "EventModel/Event.h"
#include "Identifier/Identifier.h"
#include "Identifier/MdcID.h"
#include "MdcRawEvent/MdcDigi.h"
#include "TStyle.h"
 
#include <cmath>
#include <cstring>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <string>
 
#include "TF1.h"
 
using namespace Event;
 
PreXtMdcCalib::PreXtMdcCalib() { m_fgIniTm = false; }
 
PreXtMdcCalib::~PreXtMdcCalib() {}
 
void PreXtMdcCalib::clear() {
  for ( int lay = 0; lay < MdcCalNLayer; lay++ )
  {
    delete m_grXt[lay];
    delete m_htrec[lay];
  }
  delete m_haxis;
  delete m_fdPreXt;
}
 
void PreXtMdcCalib::initialize( TObjArray* hlist, IMdcGeomSvc* mdcGeomSvc,
                                IMdcCalibFunSvc* mdcFunSvc, IMdcUtilitySvc* mdcUtilitySvc ) {
  IMessageSvc* msgSvc;
  Gaudi::svcLocator()->service( "MessageSvc", msgSvc );
  MsgStream log( msgSvc, "PreXtMdcCalib" );
  log << MSG::INFO << "PreXtMdcCalib::initialize()" << endmsg;
 
  m_hlist         = hlist;
  m_mdcGeomSvc    = mdcGeomSvc;
  m_mdcFunSvc     = mdcFunSvc;
  m_mdcUtilitySvc = mdcUtilitySvc;
 
  int lay;
 
  m_nWire  = m_mdcGeomSvc->getWireSize();
  m_nLayer = m_mdcGeomSvc->getLayerSize();
 
  m_fdPreXt = new TFolder( "PreXt", "PreXt" );
  m_hlist->Add( m_fdPreXt );
 
  m_fdNhit = new TFolder( "XtNhit", "XtNhit" );
  m_hlist->Add( m_fdNhit );
 
  m_haxis = new TH2F( "axis", "", 50, 0, 300, 50, 0, 9 );
  m_haxis->SetStats( 0 );
  m_fdPreXt->Add( m_haxis );
 
  char hname[200];
  for ( lay = 0; lay < MdcCalNLayer; lay++ )
  {
    sprintf( hname, "trec%02d", lay );
    m_htrec[lay] = new TH1F( hname, "", 310, -20, 600 );
    m_fdPreXt->Add( m_htrec[lay] );
  }
 
  m_nhitTot = new TH1F( "nhitTot", "", 43, -0.5, 42.5 );
  m_fdNhit->Add( m_nhitTot );
 
  for ( lay = 0; lay < MdcCalNLayer; lay++ )
  {
    sprintf( hname, "nhitBin%02d", lay );
    m_nhitBin[lay] = new TH1F( hname, "", 40, 5.0, 405.0 );
    m_fdNhit->Add( m_nhitBin[lay] );
  }
 
  /* integrated time Spectrum for determination X-T */
  int bin;
  m_nXtBin    = 40;
  double twid = 10.0; // ns
  for ( bin = 0; bin < m_nXtBin; bin++ ) m_tbin[bin] = (double)( bin + 1 ) * twid;
 
  for ( lay = 0; lay < MdcCalNLayer; lay++ )
  {
    m_nTot[lay] = 0;
    for ( bin = 0; bin < m_nXtBin; bin++ ) { m_nEntries[lay][bin] = 0; }
  }
}
 
int PreXtMdcCalib::fillHist( MdcCalEvent* event ) {
  IMessageSvc* msgSvc;
  Gaudi::svcLocator()->service( "MessageSvc", msgSvc );
  MsgStream log( msgSvc, "PreXtMdcCalib" );
  log << MSG::DEBUG << "PreXtMdcCalib::fillHist()" << endmsg;
 
  int        bin;
  int        lay;
  int        digiId;
  unsigned   fgOverFlow;
  double     tdc;
  double     adc;
  double     traw;
  double     trec;
  Identifier id;
 
  IDataProviderSvc* eventSvc = NULL;
  Gaudi::svcLocator()->service( "EventDataSvc", eventSvc );
 
  double xtpar[8];
  if ( !m_fgIniTm )
  {
    for ( lay = 0; lay < MdcCalNLayer; lay++ )
    {
      m_t0[lay] = m_mdcFunSvc->getT0( lay, 0 );
      m_mdcFunSvc->getXtpar( lay, 0, 0, xtpar );
      m_tm[lay] = xtpar[6];
    }
    m_fgIniTm = true;
  }
 
  // get EsTimeCol
  double                     tes        = -9999.0;
  int                        esTimeflag = -1;
  SmartDataPtr<RecEsTimeCol> aevtimeCol( eventSvc, "/Event/Recon/RecEsTimeCol" );
  if ( ( !aevtimeCol ) || ( aevtimeCol->size() == 0 ) )
  {
    tes        = -9999.0;
    esTimeflag = -1;
  }
  else
  {
    RecEsTimeCol::iterator iter_evt = aevtimeCol->begin();
    for ( ; iter_evt != aevtimeCol->end(); iter_evt++ )
    {
      tes        = ( *iter_evt )->getTest();
      esTimeflag = ( *iter_evt )->getStat();
    }
  }
  bool flagTes = false;
  for ( int iEs = 0; iEs < m_param.nEsFlag; iEs++ )
  {
    if ( esTimeflag == m_param.esFlag[iEs] )
    {
      flagTes = true;
      break;
    }
  }
  if ( -1 == esTimeflag ) return -1;
  // if( (!flagTes) || (tes < m_param.tesMin) || (tes > m_param.tesMax) ) return -1;
  if ( ( m_param.esCut ) &&
       ( ( !flagTes ) || ( tes < m_param.tesMin ) || ( tes > m_param.tesMax ) ) )
    return -1;
 
  // retrieve Mdc digi
  SmartDataPtr<MdcDigiCol> mdcDigiCol( eventSvc, "/Event/Digi/MdcDigiCol" );
  if ( !mdcDigiCol ) { log << MSG::FATAL << "Could not find event" << endmsg; }
 
  MdcDigiCol::iterator iter = mdcDigiCol->begin();
  digiId                    = 0;
  for ( ; iter != mdcDigiCol->end(); iter++, digiId++ )
  {
    MdcDigi* aDigi = ( *iter );
    id             = ( aDigi )->identify();
 
    lay = MdcID::layer( id );
 
    tdc        = ( aDigi )->getTimeChannel();
    adc        = ( aDigi )->getChargeChannel();
    fgOverFlow = ( aDigi )->getOverflow();
    traw       = tdc * MdcCalTdcCnv;
    trec       = traw - tes - m_t0[lay];
 
    //    if ( !((fgOverFlow == 0)||(fgOverFlow ==12)) ||
    //         (aDigi->getTimeChannel() == 0x7FFFFFFF) ||
    //         (aDigi->getChargeChannel() == 0x7FFFFFFF) ) continue;
    if ( ( ( fgOverFlow & 1 ) != 0 ) || ( ( fgOverFlow & 12 ) != 0 ) ||
         ( aDigi->getTimeChannel() == 0x7FFFFFFF ) ||
         ( aDigi->getChargeChannel() == 0x7FFFFFFF ) )
      continue;
 
    /* integrated time Spectrum for determination X-T */
    if ( trec > 0 )
    {
      if ( trec < m_tm[lay] )
      {
        m_nTot[lay]++;
        m_htrec[lay]->Fill( trec );
        m_nhitTot->Fill( lay );
      }
      for ( bin = 0; bin < m_nXtBin; bin++ )
      {
        if ( trec < m_tbin[bin] )
        {
          m_nEntries[lay][bin]++;
          m_nhitBin[lay]->Fill( m_tbin[bin] );
        }
      }
    }
  }
  return 1;
}
 
void PreXtMdcCalib::printCut() const {}
 
int PreXtMdcCalib::updateConst( MdcCalibConst* calconst ) {
  IMessageSvc* msgSvc;
  Gaudi::svcLocator()->service( "MessageSvc", msgSvc );
  MsgStream log( msgSvc, "PreXtMdcCalib" );
  log << MSG::DEBUG << "PreXtMdcCalib::updateConst()" << endmsg;
 
  int lay;
  int bin;
  int iLR;
  int iEntr;
  int ord;
 
  double layRad;
  double ncel;
  double pi = 3.141592653;
 
  double dm;
  double dist[40];
  double xtpar[6];
  char   hname[200];
 
  TF1* funXt = new TF1( "funXt", xtfun, 0, 300, 6 );
  funXt->FixParameter( 0, 0.0 );
  funXt->SetParameter( 1, 0.03 );
  funXt->SetParameter( 2, 0.0 );
  funXt->SetParameter( 3, 0.0 );
  funXt->SetParameter( 4, 0.0 );
  funXt->SetParameter( 5, 0.0 );
 
  ofstream fxtlog( "preXtpar.dat" );
  for ( lay = 0; lay < MdcCalNLayer; lay++ )
  {
    sprintf( hname, "grPreXt%02d", lay );
    m_grXt[lay] = new TGraph();
    m_grXt[lay]->SetName( hname );
    m_grXt[lay]->SetMarkerStyle( 20 );
    m_fdPreXt->Add( m_grXt[lay] );
 
    layRad = m_mdcGeomSvc->Layer( lay )->Radius();
    ncel   = m_mdcGeomSvc->Layer( lay )->NCell();
    dm     = pi * layRad / ncel;
 
    fxtlog << "layer " << lay << endl;
    for ( bin = 0; bin < m_nXtBin; bin++ )
    {
      dist[bin] = dm * m_nEntries[lay][bin] / m_nTot[lay];
      m_grXt[lay]->SetPoint( bin, m_tbin[bin], dist[bin] );
      fxtlog << setw( 4 ) << bin << setw( 15 ) << m_tbin[bin] << setw( 15 ) << dist[bin]
             << setw( 15 ) << dm << setw( 10 ) << m_nEntries[lay][bin] << setw( 10 )
             << m_nTot[lay] << endl;
 
      if ( m_tbin[bin] >= m_tm[lay] ) break;
    }
 
    if ( 1 == m_param.fgCalib[lay] )
    {
      m_grXt[lay]->Fit( funXt, "Q", "", 0.0, m_tm[lay] );
      for ( ord = 0; ord < 6; ord++ ) { xtpar[ord] = funXt->GetParameter( ord ); }
 
      for ( iEntr = 0; iEntr < MdcCalNENTRXT; iEntr++ )
      {
        for ( iLR = 0; iLR < MdcCalLR; iLR++ )
        {
          for ( ord = 0; ord < 6; ord++ )
          { calconst->resetXtpar( lay, iEntr, iLR, ord, xtpar[ord] ); }
        }
      }
    }
    else
    {
      for ( ord = 0; ord < 6; ord++ ) xtpar[ord] = calconst->getXtpar( lay, 0, 0, ord );
    }
 
    for ( ord = 0; ord < 6; ord++ ) fxtlog << setw( 14 ) << xtpar[ord];
    fxtlog << setw( 10 ) << m_tm[lay] << "  0" << endl;
  } // end of layer loop
  fxtlog.close();
  cout << "preXt.dat was written." << endl;
 
  delete funXt;
  return 1;
}
 
Double_t PreXtMdcCalib::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;
}