/home/bes3soft/bes3soft/Boss/7.0.2/dist/7.0.2/Mdc/MdcCalibAlg/MdcCalibAlg-00-09-02/src/T0MdcCalib.cxx

Go to the documentation of this file.

#include "MdcCalibAlg/T0MdcCalib.h"

#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/IMessageSvc.h"
#include "GaudiKernel/StatusCode.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/Bootstrap.h"

#include <iostream>
#include <iomanip>
#include <cstring>

T0MdcCalib::T0MdcCalib(){
     for(int lay=0; lay<MdcCalNLayer; lay++){
          if(lay < 8){
               m_docaMin[lay] = 1.2; // mm
               m_docaMax[lay] = 4.8; // mm
          } else{
               m_docaMin[lay] = 1.6; // mm
               m_docaMax[lay] = 6.4; // mm
          }
     }
}

T0MdcCalib::~T0MdcCalib(){
}

void T0MdcCalib::clear(){
     for(int i=0; i<MdcCalTotCell; i++){
          delete m_hleft[i];
          delete m_hright[i];
     }
     delete m_hLrResiSum;
     delete m_hLrResiSub;
     delete m_fdT0;
     delete m_fdResiWire;

     MdcCalib::clear();
}

void T0MdcCalib::initialize(TObjArray* hlist, IMdcGeomSvc* mdcGeomSvc,
                            IMdcCalibFunSvc* mdcFunSvc, IMdcUtilitySvc* mdcUtilitySvc) {
     IMessageSvc* msgSvc;
     Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
     MsgStream log(msgSvc, "T0MdcCalib");
     log << MSG::INFO << "T0MdcCalib::initialize()" << endreq;

     m_hlist = hlist;
     m_mdcGeomSvc = mdcGeomSvc;
     m_mdcFunSvc = mdcFunSvc;
     m_mdcUtilitySvc = mdcUtilitySvc;

     m_vdr = 0.03;

     MdcCalib::initialize(m_hlist, m_mdcGeomSvc, m_mdcFunSvc, m_mdcUtilitySvc);

     int i;
     int nwire;
     int lay;
     int cel;
     char hname[200];

     m_fdT0 = new TFolder("fdT0", "fdT0");
     m_hlist -> Add(m_fdT0);

     m_fdResiWire = new TFolder("ResiWire", "ResiWire");
     m_hlist->Add(m_fdResiWire);

     nwire = m_mdcGeomSvc -> getWireSize();
     for(i=0; i<nwire; i++){
          lay = m_mdcGeomSvc -> Wire(i) -> Layer();
          cel = m_mdcGeomSvc -> Wire(i) -> Cell();

          sprintf(hname, "Resi%04d_Lay%02d_Cell%03d_L", i, lay, cel);
          m_hleft[i] = new TH1F(hname, "", 400, -2.0, 2.0);
          m_fdT0 -> Add(m_hleft[i]);

          sprintf(hname, "Resi%04d_Lay%02d_Cell%03d_R", i, lay, cel);
          m_hright[i] = new TH1F(hname, "", 400, -2.0, 2.0);
          m_fdT0 -> Add(m_hright[i]);
     }

     m_hLrResiSum = new TH1F("LrResiSum", "", 200, -0.5, 0.5);
     m_fdResiWire->Add(m_hLrResiSum);

     m_hLrResiSub = new TH1F("LrResiSub", "", 200, -0.5, 0.5);
     m_fdResiWire->Add(m_hLrResiSub);
}

int T0MdcCalib::fillHist(MdcCalEvent* event){
     IMessageSvc* msgSvc;
     Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
     MsgStream log(msgSvc, "T0MdcCalib");
     log << MSG::DEBUG << "T0MdcCalib::fillHist()" << endreq;

     MdcCalib::fillHist(event);

     // get EsTimeCol
     bool esCutFg = event->getEsCutFlag();
     if( ! esCutFg ) return -1;

     int i;
     int k;
     int ntrk;
     int nhit;
     int stat;

     int lay;
     int cel;
     int wir;
     int lr;
     double dmeas;
     double doca;
     double resi;

     bool fgHitLay[MdcCalNLayer];

     MdcCalRecTrk* rectrk;
     MdcCalRecHit* rechit;

     ntrk = event -> getNTrk();
     if((ntrk < m_param.nTrkCut[0]) || (ntrk > m_param.nTrkCut[1])) return -1;

     for(i=0; i<ntrk; i++){
          rectrk = event -> getRecTrk(i);
          nhit = rectrk -> getNHits();

          // dr cut
          double dr = rectrk->getDr();
          if(fabs(dr) > m_param.drCut) continue;

          // dz cut
          double 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();
               fgHitLay[lay] = true;
          }

          int nhitlay = 0;
          for(lay=0; lay<MdcCalNLayer; lay++) if(fgHitLay[lay]) nhitlay++;
          if(nhitlay < m_param.nHitLayCut) continue;

          bool fgNoise = rectrk->getFgNoiseRatio();
          if(m_param.noiseCut && (!fgNoise)) continue;

//        log << MSG::DEBUG << "nhits: " << nhit << endreq;
          for(k=0; k<nhit; k++){
               rechit = rectrk -> getRecHit(k);
               lay = rechit -> getLayid();
               cel = rechit -> getCellid();
               wir = m_mdcGeomSvc -> Wire(lay, cel) -> Id();
               lr = rechit -> getLR();
               doca = rechit -> getDocaExc();
               dmeas = rechit -> getDmeas();
               resi = rechit -> getResiExc();
               stat = rechit -> getStat();

               if(1 != stat) continue;

               if( (fabs(doca) < m_docaMin[lay]) ||
                   (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;
               }

               if(wir < (m_mdcGeomSvc->getWireSize()) ){
                    if( 0 == lr ){
                         m_hleft[wir] -> Fill(resi);
                    }else if( 1 == lr ){
                         m_hright[wir] -> Fill(resi);
                    }
               }else{
                    std::cout << "wireid: " << wir << std::endl;
               }
          }
     }
     return 1;
}

int T0MdcCalib::updateConst(MdcCalibConst* calconst){
     IMessageSvc* msgSvc;
     Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
     MsgStream log(msgSvc, "T0MdcCalib");
     log << MSG::INFO << "T0MdcCalib::updateConst()" << endreq;

     MdcCalib::updateConst(calconst);

     int i;
     int lay;
     int nwire;
     double t0;
     double delt0;
     double resiLrSum;
     double resiLrSub;

     Stat_t entry_l;
     double mean_l;
//      double gmean_l;  // mean of fit with the gaussian distribution
//      double sigma_l;
//      double rms_l;

     Stat_t entry_r;
     double mean_r;
//      double gmean_r;
//      double sigma_r;
//      double rms_r;

     nwire = m_mdcGeomSvc -> getWireSize();
     std::cout << "totwire: " << nwire << std::endl;
     for(i=0; i<nwire; i++){
          lay = m_mdcGeomSvc -> Wire(i) -> Layer();

          if(1 == m_param.fgCalib[lay]){
               entry_l = m_hleft[i] -> GetEntries();
               mean_l = m_hleft[i] -> GetMean();
//             m_hleft[i] -> Fit("gaus", "Q");
//             gmean_l = m_hleft[i] -> GetFunction("gaus") -> GetParameter(1);

               entry_r = m_hright[i] -> GetEntries();
               mean_r = m_hright[i] -> GetMean();
//             m_hright[i] -> Fit("gaus", "Q");
//             gmean_r = m_hright[i] -> GetFunction("gaus") -> GetParameter(1);

               delt0 = 0.5 * (mean_l + mean_r) / m_vdr;
          } else{
               delt0 = 0.0;
          }

          resiLrSum = 0.5 * (mean_l + mean_r);
          resiLrSub = 0.5 * (mean_l - mean_r);
          m_hLrResiSum->Fill(resiLrSum);
          m_hLrResiSub->Fill(resiLrSub);

          t0 = calconst -> getT0(i);
          t0 += delt0;

          calconst -> resetT0(i, t0);
          calconst -> resetDelT0(i, delt0);
     }
     return 1;
}

---