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

Go to the documentation of this file.

#include "MdcCalibAlg/MdcCalRecTrk.h"

#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/IMessageSvc.h"
#include "GaudiKernel/StatusCode.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/SmartDataPtr.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "MdcRawEvent/MdcDigi.h"
#include "McTruth/MdcMcHit.h"
#include "CLHEP/Units/PhysicalConstants.h"

MdcCalRecTrk::MdcCalRecTrk(int pid){
     m_pid = pid;
}

MdcCalRecTrk::~MdcCalRecTrk(){
     unsigned int i;
     for(i=0; i<m_rechit.size(); i++){
          delete m_rechit[i];
     }
     m_rechit.clear();
}

void MdcCalRecTrk::setRecTrk(RecMdcTrackCol::iterator it_trk){
     IMessageSvc *msgSvc;
     Gaudi::svcLocator()->service("MessageSvc", msgSvc);
     MsgStream log(msgSvc, "MdcCalRecTrk");
     log << MSG::DEBUG << "MdcCalRecTrk::setRecTrk()" << endreq;

     m_dr = (*it_trk)->helix(0);
     m_phi0 = (*it_trk)->helix(1);
     m_kappa = (*it_trk)->helix(2);
     m_dz = (*it_trk)->helix(3);
     m_tanl = (*it_trk)->helix(4);
     m_chisq = (*it_trk)->chi2();
     m_nhits = (*it_trk)->getNhits();
     m_helix = (*it_trk)->helix();
     m_helixerr = (*it_trk)->err();
     m_fgNoiseRatio = fgNoiseRatio(m_phi0);

     double mass;
     if(0 == m_pid) mass = 0.000511;
     else if(1 == m_pid) mass = 0.105658;
     else if(2 == m_pid) mass = 0.139570;
     else if(3 == m_pid) mass = 0.493677;
     else if(4 == m_pid) mass = 0.938272;
     else if(5 == m_pid) mass = 0.139570; // cosmic-ray
     else log << MSG::FATAL << "wrong particle type" << endreq;
     m_p4 = (*it_trk)->p4(mass);

     m_dr *= 10.0;              // cm -> mm
     m_dz *= 10.0;              // cm -> mm

     m_pt = 1.0 / m_kappa;
     m_p = m_pt * sqrt( m_tanl * m_tanl + 1.0 );

     HitRefVec gothits = (*it_trk) -> getVecHits();
     HitRefVec::iterator it_hit = gothits.begin();
     MdcCalRecHit* rechit;
     for(; it_hit != gothits.end(); it_hit++){
          rechit = new MdcCalRecHit();
          rechit->setRecHit(it_hit);
          m_rechit.push_back(rechit);
     }
}

void MdcCalRecTrk::setKalTrk(RecMdcKalTrackCol::iterator it_trk){
     IMessageSvc *msgSvc;
     Gaudi::svcLocator()->service("MessageSvc", msgSvc);
     MsgStream log(msgSvc, "MdcCalRecTrk");
     log << MSG::DEBUG << "MdcCalRecTrk::setKalTrk()" << endreq;

     if(0 == m_pid) RecMdcKalTrack::setPidType(RecMdcKalTrack::electron);
     else if(1 == m_pid) RecMdcKalTrack::setPidType(RecMdcKalTrack::muon);
     else if(2 == m_pid) RecMdcKalTrack::setPidType(RecMdcKalTrack::pion);
     else if(3 == m_pid) RecMdcKalTrack::setPidType(RecMdcKalTrack::kaon);
     else if(4 == m_pid) RecMdcKalTrack::setPidType(RecMdcKalTrack::proton);
     else if(5 == m_pid) RecMdcKalTrack::setPidType(RecMdcKalTrack::muon); // cosmic-ray
     else log << MSG::FATAL << "wrong particle type" << endreq;

     m_dr = (*it_trk)->dr();
     m_phi0 = (*it_trk)->fi0();
     m_kappa = (*it_trk)->kappa();
     m_dz = (*it_trk)->dz();
     m_tanl = (*it_trk)->tanl();
     m_chisq = (*it_trk)->chi2();
     m_fgNoiseRatio = fgNoiseRatio(m_phi0);

     m_p4.setPx(- sin(m_phi0) / fabs(m_kappa));
     m_p4.setPy(cos(m_phi0) / fabs(m_kappa));
     m_p4.setPz(m_tanl / fabs(m_kappa));
     double p3 = m_p4.mag();
     double mass;
     if(0 == m_pid) mass = 0.000511;
     else if(1 == m_pid) mass = 0.105658;
     else if(2 == m_pid) mass = 0.139570;
     else if(3 == m_pid) mass = 0.493677;
     else if(4 == m_pid) mass = 0.938272;
     else if(5 == m_pid) mass = 0.139570; // cosmic-ray
     m_p4.setE(sqrt(p3 * p3 + mass * mass));

     m_dr *= 10.0;              // cm -> mm
     m_dz *= 10.0;              // cm -> mm

     m_pt = 1.0 / m_kappa;
     m_p = m_pt * sqrt( m_tanl * m_tanl + 1.0 );

//      cout << setw(15) << m_p << setw(15) << m_dr << setw(15) << m_phi0
//        << setw(15) << m_kappa << setw(15) << m_dz << setw(15) << m_tanl << endl;

     HelixSegRefVec gothelixsegs = (*it_trk)->getVecHelixSegs();
     HelixSegRefVec::iterator it_hit = gothelixsegs.begin();
     MdcCalRecHit* rechit;

     int k = 0;
     for(; it_hit != gothelixsegs.end(); it_hit++){
          rechit = new MdcCalRecHit();
          rechit->setKalHit(it_hit);
          m_rechit.push_back(rechit);

          k++;
     }
     m_nhits = k;
//      cout<<"hits "<<m_nhits <<endl;
}

bool MdcCalRecTrk::fgNoiseRatio(double phi0){
     IMessageSvc* msgSvc;
     Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
     MsgStream log(msgSvc, "MdcCalib");

     IDataProviderSvc* eventSvc = NULL;
     Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
     SmartDataPtr<MdcDigiCol> mdcDigiCol(eventSvc,"/Event/Digi/MdcDigiCol");
     if (!mdcDigiCol) {
          log << MSG::FATAL << "Could not find event" << endreq;
     }

     int lay;
     int cel;
     bool hitCel[43][288];
     for(lay=0; lay<43; lay++){
          for(cel=0; cel<288; cel++){
               hitCel[lay][cel] = false;
          }
     }

     MdcDigiCol::iterator iter = mdcDigiCol->begin();
     unsigned fgOverFlow;
     int digiId = 0;
     Identifier id;
     for(; iter != mdcDigiCol->end(); iter++, digiId++) {
          MdcDigi *aDigi = (*iter);
          id = (aDigi)->identify();
          lay = MdcID::layer(id);
          cel = MdcID::wire(id);
          fgOverFlow = (aDigi) -> getOverflow();

          bool goodTQ = true;
          if ( ((fgOverFlow & 3) !=0 ) || ((fgOverFlow & 12) != 0) ||
               (aDigi->getTimeChannel() == 0x7FFFFFFF) ||
               (aDigi->getChargeChannel() == 0x7FFFFFFF) ) goodTQ = false;

          bool goodT = true;
          if ( ((fgOverFlow & 1) !=0 ) || (aDigi->getTimeChannel() == 0x7FFFFFFF) ) goodT = false;

//        if(!goodTQ) continue;
          if(!goodT) continue;

          hitCel[lay][cel] = true;
     }

     SmartDataPtr<RecMdcTrackCol> newtrkCol(eventSvc, "/Event/Recon/RecMdcTrackCol");
     if(!newtrkCol){
          log << MSG::ERROR << "Could not find RecMdcTrackCol" << endreq;
          return -1;
     }

     int nNoiRange = 4;         // number of cells
     double dphi = 1.0;
     RecMdcTrackCol::iterator it_trk = newtrkCol->begin();
     for(it_trk = newtrkCol->begin(); it_trk != newtrkCol->end(); it_trk++){
          HitRefVec gothits = (*it_trk) -> getVecHits();
          HitRefVec::iterator it_hit = gothits.begin();
          HepVector helix = (*it_trk)->helix();
          HepSymMatrix helixErr = (*it_trk)->err();
          double phi0Rec = (*it_trk)->helix(1);
          double delphi;
          if(phi0Rec>phi0) delphi = phi0Rec - phi0;
          else delphi = phi0 - phi0Rec;
          if(delphi > CLHEP::pi) delphi -= CLHEP::pi;
          if(delphi > (CLHEP::pi*0.17)) continue;

          int nhitLay1 = 0;     // lay0-3
          int nhitLay2 = 0;     // lay4-7
          int nhitLay3 = 0;     // lay8-19
          int nhitLay4 = 0;     // lay20-42
          int nhitT1 = 0;
          int nhitT2 = 0;
          int nhitT3 = 0;
          int nhitT4 = 0;
          for(lay=0; lay<8; lay++){
               double docamin = 0.7; // cm
               if(lay>7) docamin = 0.9;
               int celmin = -1;
               int ncel = m_mdcGeomSvc->Layer(lay)->NCell();
               for(cel=0; cel<ncel; cel++){
                    double wphi;
                    const MdcGeoWire* pWire = m_mdcGeomSvc -> Wire(lay, cel);
                    double xx = pWire->Backward().x();
                    double yy = pWire->Backward().y();
                    double rr = sqrt( (xx * xx) + (yy * yy) );
                    if( yy >= 0 )  wphi = acos(xx / rr);
                    else           wphi = CLHEP::twopi - acos(xx / rr);

                    double phiTrk = phi0 + CLHEP::halfpi;
                    if(phiTrk > CLHEP::twopi) phiTrk -= CLHEP::twopi;
                    if( !( (fabs(wphi-phiTrk) < dphi) || (fabs(wphi+CLHEP::twopi-phiTrk) < dphi) ||
                           (fabs(wphi-CLHEP::twopi-phiTrk) < dphi) ) ){
                         continue;
                    }

                    double doca = m_mdcUtilitySvc->doca(lay, cel, helix, helixErr);
//                  cout << setw(5) << lay << setw(5) << cel << setw(15) << doca << endl;
                    if(fabs(doca) < fabs(docamin)){
                         docamin = doca;
                         celmin = cel;
                    }
               }
               if(celmin > -1){
                    if(lay<4) nhitLay1++;
                    else if(lay<8) nhitLay2++;
                    else if(lay<20) nhitLay3++;
                    else nhitLay4++;
                    for(int ii=(-nNoiRange); ii<=nNoiRange; ii++){
                         if(0==ii) continue;
                         int icell = celmin + ii;
                         if(icell >= ncel) icell -= ncel;
                         if(icell < 0) icell += ncel;
//                       cout << "hit " << setw(5) << lay << setw(5) << celmin << setw(5) << icell << setw(5) << hitCel[lay][icell]<<endl;
                         if(hitCel[lay][icell]){
                              if(lay<4) nhitT1++;
                              else if(lay<8) nhitT2++;
                              else if(lay<20) nhitT3++;
                              else nhitT4++;
                         }
                    }
               }
          }
          if((nhitLay1<=0) || (nhitLay2<=0) || (nhitLay3<=0) || (nhitLay4<=0)) return false;
          double ratio1 = (double)nhitT1 / ((double)nhitLay1 * (double)nNoiRange*2.0);
          double ratio2 = (double)nhitT2 / ((double)nhitLay2 * (double)nNoiRange*2.0);
          double ratio3 = (double)nhitT3 / ((double)nhitLay3 * (double)nNoiRange*2.0);
          double ratio4 = (double)nhitT4 / ((double)nhitLay4 * (double)nNoiRange*2.0);

          if((ratio1>0.08) || (ratio2>0.08) || (ratio3>0.03) || (ratio4>0.03)) return false;
          else return true;
     }
     return false;
}

---