BOSS_7.0.2: /home/bes3soft/bes3soft/Boss/7.0.2/dist/7.0.2/Mdc/MdcCalibAlg/MdcCalibAlg-00-09-02/src/MdcCalibAlg.cxx Source File

00001 // MdcCalibAlg
00002 // 2006/01/09   Wu Linghui   IHEP
00003 
00004 #include "GaudiKernel/MsgStream.h"
00005 #include "GaudiKernel/StatusCode.h"
00006 
00007 #include "MdcCalibAlg/MdcCalibAlg.h"
00008 #include "EventModel/EventHeader.h"
00009 
00010 #include <iostream>
00011 #include <fstream>
00012 #include <cstring>
00013 
00014 #include "TStyle.h"
00015 
00016 using namespace Event;
00017 
00019 
00020 MdcCalibAlg::MdcCalibAlg(const std::string& name, ISvcLocator* pSvcLocator) :
00021      Algorithm(name, pSvcLocator), m_mdcCalFlg(0), m_flgKalFit(0), m_evtType(0){
00022 
00023      m_histname = "NULL";
00024      m_configFile = "NULL";
00025      m_nEvtDisp = 1000;
00026      m_distCalib = false;
00027      m_nEvt = 0;
00028 
00029      // declare properties
00030      declareProperty("MdcCalFlg", m_mdcCalFlg);
00031      declareProperty("MdcKalmanFlg", m_flgKalFit);
00032      declareProperty("Event", m_evtType);
00033      declareProperty("HistOutput", m_histname);
00034      declareProperty("ConfigFile", m_configFile);
00035      declareProperty("WirePosCorFile", m_wpcFile);
00036      declareProperty("NumEvtDisplay", m_nEvtDisp);
00037      declareProperty("DistCalib", m_distCalib);
00038      declareProperty("Ecm", m_ecm = 3.686);
00039 
00040      m_initCalConstFlg = false;
00041 }
00042 
00043 MdcCalibAlg::~MdcCalibAlg(){
00044      delete m_constmgr;
00045      std::cout << "m_constmgr deleted" << std::endl;
00046 
00047      delete m_mdccalib;
00048      std::cout << "delete m_mdccalib" << std::endl;
00049 
00050      delete m_mdcevt;
00051      std::cout << "m_mdcevt deleted" << std::endl;
00052 
00053      delete m_hlist;
00054      std::cout << "m_hlist deleted" << std::endl;
00055 
00056      delete m_calconst;
00057      std::cout << "m_calconst deleted" << std::endl;
00058 
00059      if(!m_distCalib){
00060           std::ofstream fend("endcal.out");
00061           fend << "MdcCalib end." << std::endl;
00062           fend.close();
00063      }
00064 
00065      std::cout << "MdcCalibAlg End." << std::endl;
00066 }
00067 
00068 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
00069 StatusCode MdcCalibAlg::initialize() {
00070      MsgStream log( msgSvc(), name() );
00071      log << MSG::INFO << "MdcCalibAlg initialze() ..." << endreq;
00072      log << MSG::INFO << "MdcCalFlg = " << m_mdcCalFlg << endreq;
00073 
00074      if("NULL"==m_histname){
00075           log << MSG::ERROR << "not defined histogram file." << endreq;
00076           return StatusCode::FAILURE;
00077      }
00078      if("NULL"==m_configFile){
00079           log << MSG::ERROR << "not defined MdcCalibConfig file." << endreq;
00080           return StatusCode::FAILURE;
00081      }
00082 
00083      StatusCode sc = service("MdcGeomSvc", m_mdcGeomSvc);
00084      if(sc != StatusCode::SUCCESS){
00085           log << MSG::ERROR << "can not use MdcGeomSvc" << endreq;
00086      }
00087 
00088      sc = service("MdcCalibFunSvc", m_mdcFunSvc);
00089      if( sc != StatusCode::SUCCESS ){
00090           log << MSG::FATAL << "can not use MdcCalibFunSvc" << endreq;
00091      }
00092 
00093      sc = service ("MdcUtilitySvc", m_mdcUtilitySvc);
00094      if ( sc.isFailure() ){
00095           log << MSG::FATAL << "Could not load MdcUtilitySvc!" << endreq;
00096      }
00097 
00098      string estver = getenv("ESTIMEALGROOT");
00099      cout << "EsTimeAlg_ver: "<< estver << endl;
00100 
00101      // initialize m_param
00102      initParam();
00103 
00104      // read mdc config parameters
00105      int i;
00106      int lay;
00107      std::string strconfig;
00108      std::string strcomment;
00109      std::string strtmp;
00110      std::string strTes;
00111      int fgTes[50];
00112      for(i=0; i<50; i++) fgTes[i] = -999;
00113      ifstream fconfig( m_configFile.c_str() );
00114      if( ! fconfig.is_open() ){
00115           log << MSG::ERROR << "can not open config file " << m_configFile
00116               << ". Use defalt config parameters" << endreq;
00117           return StatusCode::FAILURE;
00118      } else {
00119           log << MSG::INFO << "Open config file " << m_configFile << endreq;
00120           while( fconfig >> strconfig ){
00121                if('#' == strconfig[0]){
00122                     getline(fconfig, strcomment);
00123                } else if("FillNtuple" == strconfig){
00124                     fconfig >> m_param.fillNtuple;
00125                } else if("nEvtNtuple" == strconfig){
00126                     fconfig >> m_param.nEvtNtuple;
00127                }else if("FlagCalDetEffi" == strconfig){
00128                     fconfig >> m_param.fgCalDetEffi;
00129                }else if("BoostPar" == strconfig){
00130                     fconfig >> m_param.boostPar[0] >> m_param.boostPar[1] >> m_param.boostPar[2];
00131                } else if("EsTimeFlag" == strconfig){
00132                     getline(fconfig, strTes);
00133                     int n = sscanf(strTes.c_str(), "%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d",
00134                                    fgTes+0, fgTes+1, fgTes+2, fgTes+3, fgTes+4,
00135                                    fgTes+5, fgTes+6, fgTes+7, fgTes+8, fgTes+9,
00136                                    fgTes+10, fgTes+11, fgTes+12, fgTes+13, fgTes+14,
00137                                    fgTes+15, fgTes+16, fgTes+17, fgTes+18, fgTes+19);
00138                     for(i=0; i<n; i++) m_param.esFlag[i] = fgTes[i];
00139                     m_param.nEsFlag = n;
00140                } else if("TimeShift" == strconfig){
00141                     fconfig >> m_param.timeShift;
00142                } else if("TesMin" == strconfig){
00143                     fconfig >> m_param.tesMin;
00144                } else if("TesMax" == strconfig){
00145                     fconfig >> m_param.tesMax;
00146                } else if("FlagIniCalConst" == strconfig){
00147                     fconfig >> m_param.fgIniCalConst;
00148                } else if("FlagUpdateTmInPreT0" == strconfig){
00149                     fconfig >> m_param.preT0SetTm;
00150                } else if("InitT0" == strconfig){
00151                     fconfig >> m_param.initT0;
00152                } else if("T0Shift" == strconfig){
00153                     fconfig >> m_param.t0Shift;
00154                } else if("TminFitChindf" == strconfig){
00155                     fconfig >> m_param.tminFitChindf;
00156                } else if("TmaxFitChindf" == strconfig){
00157                     fconfig >> m_param.tmaxFitChindf;
00158                } else if("InitSigma" == strconfig){
00159                     fconfig >> m_param.initSigma;
00160                } else if("UpdateSigma" == strconfig){
00161                     fconfig >> m_param.calSigma;
00162                } else if("ResidualType" == strconfig){
00163                     fconfig >> m_param.resiType;
00164                } else if("FixXtC0" == strconfig){
00165                     fconfig >> m_param.fixXtC0;
00166                } else if("FixXtEdge" == strconfig){
00167                     fconfig >> m_param.fixXtEdge;
00168                } else if("FlagAdjacLayerCut" == strconfig){
00169                     fconfig >> m_param.fgAdjacLayerCut;
00170                } else if("FlagBoundLayerCut" == strconfig){
00171                     fconfig >> m_param.fgBoundLayerCut;
00172                } else if("hitStatCut" == strconfig){
00173                     fconfig >> m_param.hitStatCut;
00174                } else if("nTrkCut" == strconfig){
00175                     fconfig >> m_param.nTrkCut[0] >> m_param.nTrkCut[1];
00176                } else if("nHitLayCut" == strconfig){
00177                     fconfig >> m_param.nHitLayCut;
00178                } else if("nHitCut" == strconfig){
00179                     fconfig >> m_param.nHitCut;
00180                } else if("noiseCut" == strconfig){
00181                     fconfig >> m_param.noiseCut;
00182                } else if("cosThetaCut" == strconfig){
00183                     fconfig >> m_param.costheCut[0] >> m_param.costheCut[1];
00184                } else if("DrCut" == strconfig){
00185                     fconfig >> m_param.drCut;
00186                } else if("DzCut" == strconfig){
00187                     fconfig >> m_param.dzCut;
00188                }else if("MaximalDocaInner" == strconfig){
00189                     fconfig >> m_param.maxDocaInner;
00190                }else if("MaximalDocaOuter" == strconfig){
00191                     fconfig >> m_param.maxDocaOuter;
00192                } else if("RawTimeFitRange" == strconfig){
00193                     for(lay=0; lay<MdcCalNLayer; lay++){
00194                          fconfig >> strtmp
00195                                  >> m_param.fgCalib[lay]
00196                                  >> m_param.tminFitRange[lay][0]
00197                                  >> m_param.tminFitRange[lay][1]
00198                                  >> m_param.tmaxFitRange[lay][0]
00199                                  >> m_param.tmaxFitRange[lay][1]
00200                                  >> m_param.initTm[lay]
00201                                  >> m_param.resiCut[lay]
00202                                  >> m_param.qmin[lay]
00203                                  >> m_param.qmax[lay];
00204                     }
00205                }
00206           }
00207           fconfig.close();
00208      }
00209      // set single wire position calibration file
00210      m_param.wpcFile = m_wpcFile;
00211 
00212      // set event type
00213      m_param.particle = m_evtType;
00214 
00215      // set Ecm
00216      m_param.ecm = m_ecm;
00217 
00218      cout << "EsFlag for Mdc Calibration: ";
00219      for(int iEs=0; iEs<m_param.nEsFlag; iEs++) cout << setw(6) << m_param.esFlag[iEs];
00220      cout << endl;
00221 
00222      m_hlist    = new TObjArray(0);
00223      m_constmgr = new MdcCalConstMgr();
00224      m_calconst = new MdcCalibConst();
00225 
00226      m_mdcevt   = new MdcCalEvent();
00227      m_mdcevt -> setParam(m_param);
00228      m_mdcevt -> setGeomSvc(m_mdcGeomSvc);
00229      m_mdcevt -> setUtilSvc(m_mdcUtilitySvc);
00230 
00231      if( 0 == m_mdcCalFlg ){
00232           m_mdccalib  = new IniMdcCalib();
00233      }else if( 1 == m_mdcCalFlg ){
00234           m_mdccalib  = new PreXtMdcCalib();
00235      }else if( 2 == m_mdcCalFlg ){
00236           m_mdccalib  = new PreT0MdcCalib();
00237      }else if( 3 == m_mdcCalFlg ){
00238           m_mdccalib  = new XtMdcCalib();
00239      }else if( 4 == m_mdcCalFlg ){
00240           m_mdccalib  = new GrXtMdcCalib();
00241      }else if( 9 == m_mdcCalFlg ){
00242           m_mdccalib  = new XtInteMdcCalib();
00243      }else if( 5 == m_mdcCalFlg ){
00244           m_mdccalib  = new T0MdcCalib();
00245      }else if( 6 == m_mdcCalFlg ){
00246           m_mdccalib  = new WrMdcCalib();
00247      }else if( 7 == m_mdcCalFlg ){
00248           m_mdccalib  = new Wr2dMdcCalib();
00249      }else if( 8 == m_mdcCalFlg ){
00250           m_mdccalib  = new QtMdcCalib();
00251      }
00252      // initialize and read the calibraion data from TCDS
00253      m_constmgr -> init(m_mdcGeomSvc, m_mdcFunSvc);
00254 
00255      m_mdccalib -> setParam(m_param);
00256      m_mdccalib -> initialize(m_hlist, m_mdcGeomSvc, m_mdcFunSvc, m_mdcUtilitySvc);
00257 
00258 //      if(0 == m_mdcCalFlg) m_calconst->initCalibConst();
00259 
00260      return StatusCode::SUCCESS;
00261 }
00262 
00263 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
00264 StatusCode MdcCalibAlg::execute() {
00265      MsgStream log(msgSvc(), name());
00266      log << MSG::DEBUG << "MdcCalibAlg execute()" << endreq;
00267 
00268      // display event number for debug
00269      SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00270      if (!eventHeader) {
00271           log << MSG::FATAL << "Could not find Event Header" << endreq;
00272           return( StatusCode::FAILURE);
00273      }
00274      int iRun = eventHeader->runNumber();
00275      int iEvt = eventHeader->eventNumber();
00276      if(0 == (m_nEvt % m_nEvtDisp))
00277           std::cout << "Run " << iRun << ", Event " << iEvt << ",   Total Event number " << m_nEvt << endl;
00278 
00279      m_mdcevt->setEvtNoOnline(iEvt);
00280      m_mdcevt->setEvtNoOffline(m_nEvt);
00281 
00282      if( ! m_initCalConstFlg ){
00283 //        if((0 == m_mdcCalFlg) && (0 == m_param.fgIniCalConst)){
00284 //             m_calconst->initCalibConst();
00285 //        } else{
00286 //             m_constmgr -> rdConstTcds(m_calconst);
00287 //        }
00288           m_constmgr -> rdConstTcds(m_calconst);
00289           m_initCalConstFlg = true;
00290      }
00291 
00292      if(m_mdcCalFlg > 1){
00293           if(m_flgKalFit)  m_mdcevt -> setKalEvent();
00294           else             m_mdcevt -> setRecEvent();
00295      }
00296 
00297      // fill histograms
00298      m_mdccalib -> fillHist(m_mdcevt);
00299      m_mdcevt   -> clear();
00300      m_nEvt++;
00301 
00302      return StatusCode::SUCCESS;
00303 }
00304 
00305 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
00306 StatusCode MdcCalibAlg::finalize() {
00307      MsgStream log(msgSvc(), name());
00308      log << MSG::INFO << "MdcCalibAlg finalize()" << endreq;
00309 
00310      gStyle -> SetCanvasBorderMode(0);
00311      gStyle -> SetCanvasColor(10);
00312      gStyle -> SetOptFit(0011);
00313      gStyle -> SetStatColor(10);
00314      gStyle -> SetStatBorderSize(1);
00315      gStyle -> SetStatFont(42);
00316      gStyle -> SetStatFontSize(0.04);
00317      gStyle -> SetStatX(0.9);
00318      gStyle -> SetStatY(0.9);
00319      gStyle -> SetPadColor(10);
00320      gStyle -> SetFuncColor(2);
00321 
00322      // execute calibration and write new calibration data file
00323      if(!m_distCalib){
00324           m_mdccalib -> updateConst(m_calconst);
00325           m_constmgr -> wrtConst( m_calconst );
00326      }
00327 
00328      // write histogram file
00329      TFile fhist(m_histname.c_str(), "recreate");
00330      m_hlist -> Write();
00331      fhist.Close();
00332      std::cout << m_histname << " was written" << std::endl;
00333 
00334      m_mdccalib->clear();
00335 
00336      return StatusCode::SUCCESS;
00337 }
00338 
00339 void MdcCalibAlg::initParam(){
00340      m_param.fillNtuple = 0;
00341      m_param.nEvtNtuple = 10000;
00342      m_param.fgCalDetEffi = 0;
00343      m_param.ecm = 3.686;
00344      m_param.boostPar[0] = 0.011;
00345      m_param.boostPar[1] = 0.0;
00346      m_param.boostPar[2] = 0.0;
00347      m_param.particle = 0;
00348      m_param.nEsFlag = 0;
00349      for(int i=0; i<50; i++) m_param.esFlag[i] = -999;
00350      m_param.timeShift = 0.0;
00351      m_param.tesMin = 0.0;
00352      m_param.tesMax = 9999.0;
00353      m_param.fgIniCalConst = 2;
00354      m_param.preT0SetTm = true;
00355      m_param.initT0 = 50.0;
00356      m_param.timeShift = 0.0;
00357      m_param.t0Shift = 0.0;
00358      m_param.tminFitChindf = 20.0;
00359      m_param.tmaxFitChindf = 20.0;
00360      m_param.initSigma = 0.16;  // mm
00361      m_param.resiType = 1;
00362      m_param.resiType = 0;
00363      m_param.fixXtC0 = 0;
00364      m_param.fixXtEdge = 0;
00365      m_param.fgAdjacLayerCut = 0;
00366      m_param.fgBoundLayerCut = 0;
00367      m_param.hitStatCut = 0;
00368      m_param.nTrkCut[0] = 2;
00369      m_param.nTrkCut[1] = 2;
00370      m_param.nHitLayCut = 35;
00371      m_param.nHitCut = 0;
00372      m_param.noiseCut = false;
00373      m_param.costheCut[0] = -1.0;
00374      m_param.costheCut[1] = 1.0;
00375      m_param.drCut = 50.0;      // mm
00376      m_param.dzCut = 300.0;     // mm
00377      m_param.maxDocaInner = 8.0;        // mm
00378      m_param.maxDocaOuter = 12.0;       // mm
00379 
00380      for(int lay=0; lay<MdcCalNLayer; lay++){
00381           m_param.fgCalib[lay] = 1;
00382           m_param.tminFitRange[lay][0] = 0.0;
00383           m_param.tminFitRange[lay][1] = 100.0;
00384           m_param.tmaxFitRange[lay][0] = 120.0;
00385           m_param.tmaxFitRange[lay][1] = 400.0;
00386           m_param.initTm[lay] = 300.0;
00387           m_param.resiCut[lay] = 1.2;
00388           m_param.qmin[lay] = 400;
00389           m_param.qmax[lay] = 1000;
00390      }
00391 }