CgemCalibFunSvc.cxx

Go to the documentation of this file.

#include "CgemCalibFunSvc/CgemCalibFunSvc.h"
#include "GaudiKernel/Kernel.h"
#include "GaudiKernel/IInterface.h"
#include "GaudiKernel/StatusCode.h"
#include "GaudiKernel/SvcFactory.h"
#include "GaudiKernel/MsgStream.h"
 
#include "GaudiKernel/IIncidentSvc.h"
#include "GaudiKernel/Incident.h"
#include "GaudiKernel/IIncidentListener.h"
 
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/Bootstrap.h"
 
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/SmartDataPtr.h"
#include "GaudiKernel/DataSvc.h"
#include "GaudiKernel/IJobOptionsSvc.h"
 
#include "CalibData/CalibModel.h"
// #include "CalibData/Cgem/CgemCalibData.h"
// #include "CalibData/Cgem/CgemDataConst.h"
#include "EventModel/EventHeader.h"
#include "GaudiKernel/SmartDataPtr.h"
 
#include <iomanip>
#include <iostream>
#include <fstream>
#include <cmath>
#include <cstring>
 
#include "TFile.h"
#include "TTree.h"
 
using namespace std;
 
CgemCalibFunSvc::CgemCalibFunSvc( const string& name, ISvcLocator* svcloc) :Service (name, svcloc){
     declareProperty("PrintInfo", m_printInfo=false);
     declareProperty("TimeFitFile", m_timeFitFile);
     declareProperty("TimeFitFilePath", m_timeFitFilePath);
     declareProperty("TimeWalkCalibFile", m_timeWalkFileName);
     declareProperty("LUTFile", m_fileLUT);
}
 
CgemCalibFunSvc::~CgemCalibFunSvc(){
}
 
StatusCode CgemCalibFunSvc::queryInterface(const InterfaceID& riid, void** ppvInterface){
     if( IID_ICgemCalibFunSvc.versionMatch(riid) ){
      *ppvInterface = static_cast<ICgemCalibFunSvc*> (this);
     } else{
      return Service::queryInterface(riid, ppvInterface);
     }
     return StatusCode::SUCCESS;
}
 
StatusCode CgemCalibFunSvc::initialize(){
     MsgStream log(messageService(), name());
     log << MSG::INFO << "CgemCalibFunSvc::initialize()" << endreq;
 
     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);
     // }
 
     m_getRunNo = false;
     static IJobOptionsSvc* jobSvc = 0;
     if ( jobSvc == 0 ) {
      sc = service("JobOptionsSvc", jobSvc);
      if ( sc.isFailure() ) {
           std::cout << "Can't get the JobOptionsSvc @ DistBoss::GetPropertyValue()" << std::endl;
           return sc;
      }
     }
     const vector<const Property*>* properties = jobSvc->getProperties("ReadCosmicRayData");
     if ( properties == NULL ) {
      log << MSG::WARNING << "can not get runNo from ReadCosmicRayData" << endreq;
     } else{
      for ( unsigned int i = 0; i < properties->size(); ++i ) {
           if ( properties->at(i)->name() == "runNo" ) {
            cout << "name of property: " << properties->at(i)->name() << endl;
            string strRnd = properties->at(i)->toString();
            sscanf(strRnd.c_str(), "%d", &m_run);
            // cout << "runNo: " << m_run << endl;
            break;
           }
      }
      m_getRunNo = true;
     }
     if( !initCalibConst() ){
          log << MSG::ERROR << "can not initialize CGEM calibration constants" << endreq;
      return StatusCode::FAILURE;
     }
 
     sc = service("CalibDataSvc", m_pCalDataSvc, true);
     if( sc == StatusCode::SUCCESS ){
          log << MSG::INFO << "Retrieve IDataProviderSvc" << endreq;
     }else{
          log << MSG::FATAL << "can not get IDataProviderSvc" << endreq;
     }
 
     m_readLUT = true;
     m_cgemLUTReader = new CgemLUTReader(m_fileLUT.c_str());
     cout << "LUT file: " << m_fileLUT << endl;
     if(!m_cgemLUTReader->ReadLUT()){
      m_readLUT = false;
      log << MSG::FATAL << "can not read LUT " << m_fileLUT << endl;
     }
 
     return StatusCode::SUCCESS;
}
 
StatusCode CgemCalibFunSvc::finalize(){
     MsgStream log(messageService(), name());
     log << MSG::INFO << "CgemCalibFunSvc::finalize()" << endreq;
 
     return StatusCode::SUCCESS;
}
 
void CgemCalibFunSvc::handle(const Incident& inc){
     MsgStream log( messageService(), name() );
     log << MSG::DEBUG << "handle: " << inc.type() << endreq;
 
     if ( inc.type() == "NewRun" ){
          log << MSG::DEBUG << "NewRun" << endreq;
 
          if( ! initCalibConst() ){
               log << MSG::ERROR 
               << "can not initilize Cgem Calib Constants" << endl
               << "  Please insert the following statement "
               << "in your \"jobOption.txt\" "
               << "before the include file of Cgem Reconstruction: "
               << endl << "        "
               << "#include \"$CALIBSVCROOT/share/job-CalibData.txt\""
               << endreq;
          }
     } 
}
 
bool CgemCalibFunSvc::initCalibConst(){
     MsgStream log(messageService(), name());
     log << MSG::INFO << "read calib const from TCDS" << endreq;
 
     // IDataProviderSvc* eventSvc = NULL;
     // Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
     // SmartDataPtr<Event::EventHeader> eventHeader(eventSvc,"/Event/EventHeader");
     // if (!eventHeader) {
     //       log << MSG::FATAL << "Could not find Event Header" << endreq;
     //       return( StatusCode::FAILURE);
     // }
     // m_run = eventHeader->runNumber();
 
     // clear calibconst
 
     // string fullPath = "/Calib/CgemCal";
     // SmartDataPtr<CalibData::CgemCalibData> calConst(m_pCalDataSvc, fullPath);
     // if( ! calConst ){
     //       log << MSG::ERROR << "can not get CgemCalibConst via SmartPtr" 
     //           << endreq;
     //       return false;
     // }
 
     char fname[500];
     if(!m_timeFitFile.empty()){
      sprintf(fname, "%s", m_timeFitFile.c_str());
     } else if(m_getRunNo && m_timeFitFile.empty()){
      sprintf(fname, "%s/timeFit_Run%d.txt", m_timeFitFilePath.c_str(), m_run);
     } else{
      cout << "ERROR: can not get runNo and CGEM calibration file" << endl;
      return false;
     }
 
     cout << "Open CGEM calibration file: " << fname << endl;
     ifstream fconst(fname);
     if(!fconst){
          log << MSG::ERROR << "can not find calibration constants. Please check runNo or calibration file" << endreq;
      return false;
     }
 
     string str;
     getline(fconst, str);
     for(int layer=0; layer<NLAYER; layer++){
      if(0==layer){
           for(int i=0; i<5; i++) fconst >> str;
           for(int xv=0; xv<NXV; xv++){
            fconst >> str >> m_tRising[layer][xv][0] >> m_tFalling[layer][xv][0];
            m_tRising[layer][xv][1] = m_tRising[layer][xv][0];
            m_tFalling[layer][xv][1] = m_tFalling[layer][xv][0];
           }
      } else{
           for(int i=0; i<7; i++) fconst >> str;
           for(int xv=0; xv<NXV; xv++){
            fconst >> str >> m_tRising[layer][xv][0] >> m_tFalling[layer][xv][0]
               >> m_tRising[layer][xv][1] >> m_tFalling[layer][xv][1];
           }
      }
     }
     fconst.close();
 
     if(m_printInfo){
      cout << "print CGEM time fit calibration constants" << endl;
      for(int layer=0; layer<NLAYER; layer++){
           cout << "layer " << layer << endl;
           for(int xv=0; xv<NXV; xv++){
            for(int sheet=0; sheet<NSHEET; sheet++){
             if((0==layer) && (1==sheet)) continue;
             cout << setw(15) << m_tRising[layer][xv][sheet] << setw(15) << m_tFalling[layer][xv][sheet];
            }
           }
           cout << endl;
      }
     }
 
     // read time walk calibration constants
     Float_t tw, charge, threshold;
     TFile ftw(m_timeWalkFileName.c_str());
     cout << "time walk file: " << m_timeWalkFileName << endl;
     TTree* trTw = (TTree*)ftw.Get("tree");
     trTw->SetBranchAddress("time_correction", &tw);
     trTw->SetBranchAddress("charge", &charge);
     trTw->SetBranchAddress("threshold", &threshold);
     int nEntries = trTw->GetEntries();
     for(int i=0; i<nEntries; i++){
      trTw->GetEntry(i);
      m_tw.push_back((double)tw);
      m_charge.push_back((double)charge);
      m_threshold.push_back((double)threshold);
 
      if(0==i){
           m_valThreshold.push_back(threshold);
           m_thIdFrom.push_back(i);
      } else{
           bool pushNewThreshold = true;
           int vecSize = m_valThreshold.size();
           for(int k=0; k<vecSize; k++){
            if(fabs(threshold-m_valThreshold[k])<0.1){
             pushNewThreshold = false;
             break;
            }
           }
           if(pushNewThreshold){
            m_valThreshold.push_back(threshold);
            m_thIdEnd.push_back(i-1);
            m_thIdFrom.push_back(i);
           }
      }
     }
     ftw.Close();
     m_thIdEnd.push_back(nEntries-1);
     if(m_printInfo){
      cout << "Print time walk calibration info" << endl;
      // for(int i=0; i<m_tw.size(); i++){
      //      cout << "time walk: " << setw(15)<<m_tw[i]<<setw(15)<<m_charge[i]
      //        <<setw(15)<<m_threshold[i]<<endl;
      // }
      for(int i=0; i<m_valThreshold.size(); i++){
           cout << setw(10) << m_valThreshold[i] << setw(10) << m_thIdFrom[i]
            << setw(10) << m_thIdEnd[i] << setw(15)<<m_charge[m_thIdFrom[i]]
            << setw(15) << m_charge[m_thIdEnd[i]] << endl;
      }
     }
 
     return true;
}
 
double CgemCalibFunSvc::getSigma(int layer, int xvFlag, int readoutMode, double angle, double Q, double T) const{
     double sigma = 0.13;   // unit: mm
     return sigma;
}
 
double CgemCalibFunSvc::getTimeRising(int layer, int xvFlag, int sheet, int stripID, double Q, double z) const{
     return m_tRising[layer][xvFlag][sheet];
}
 
double CgemCalibFunSvc::getTimeFalling(int layer, int xvFlag, int sheet, int stripID, double Q, double z) const{
     return m_tFalling[layer][xvFlag][sheet];
}
 
double CgemCalibFunSvc::getTimeWalk(int layer, int xvFlag, int sheet, int stripID, double Q) const{
     double tw = 0.0;
     if(!m_readLUT){
      cout << "ERROR: can not calculate time walk with LUT, return tw=0." << endl;
      return tw;
     }
     double threshold = m_cgemLUTReader->Get_thr_T_fC(layer, sheet, xvFlag, stripID);
     if(m_printInfo){
      cout << setw(8)<<layer <<setw(8)<<sheet<<setw(8)<<xvFlag<<setw(8)<<stripID
           <<setw(15)<<threshold<<endl;
     }
 
     tw = getTimeWalk(Q, threshold);
     return tw;
}
 
double CgemCalibFunSvc::getTimeWalk(double Q, double threshold) const{
     MsgStream log(messageService(), name());
 
     double tw = 0.0;
     int thresholdId = -1;
     double minval = 1000.;
     for(int i=0; i<m_valThreshold.size(); i++){
      double val = fabs(threshold-m_valThreshold[i]);
      if(val<minval){
           minval = val;
           thresholdId = i;
      }
     }
 
     int binMin = m_thIdFrom[thresholdId];
     int binMax = m_thIdEnd[thresholdId];
     double qMin = m_charge[binMin];
     double qMax = m_charge[binMax];
 
     if(minval>1.0){
      log << MSG::WARNING << "can not find time walk calibration constants for threshold="
          << threshold << " , use the parameters with threshold=" << m_valThreshold[thresholdId] << endreq;
      if(m_printInfo){
           cout << "threshold_binId = " << thresholdId << "   qMin = " << qMin
            << "    qMax = " << qMax << endl;
      }
     }
 
     if(Q < 0.0){
      log << MSG::WARNING << "charge<0, calculate time walk with charge=0" << endreq;
      tw = (0.0-m_charge[binMin])*(m_tw[binMin]-m_tw[binMin+1])/(m_charge[binMin]-m_charge[binMin+1]) + m_tw[binMin];
     } else if(Q < qMin){
      log << MSG::WARNING << "charge < minimumCharge in time walk calibration function " << qMin << endreq;
      tw = (Q-m_charge[binMin])*(m_tw[binMin]-m_tw[binMin+1])/(m_charge[binMin]-m_charge[binMin+1]) + m_tw[binMin];
     } else if(Q > qMax){
      log << MSG::WARNING << "charge > maximumCharge in time walk calibration function " << qMax << endreq;
      if(m_printInfo) cout << "binMin = " << binMin << "   binMax = " << binMax << endl;
      tw = (Q-m_charge[binMax])*(m_tw[binMax]-m_tw[binMax-1])/(m_charge[binMax]-m_charge[binMax-1]) + m_tw[binMax];
     } else{
      for(int i=binMin; i<binMax; i++){
           if((Q>=m_charge[i]) && (Q<m_charge[i+1])){
            tw = (m_tw[i+1]-m_tw[i])*(Q-m_charge[i])/(m_charge[i+1]-m_charge[i]) + m_tw[i];
           }
      }
     }
 
     return tw;
}