/home/bes3soft/bes3soft/Boss/7.0.2/dist/7.0.2/Reconstruction/TrkReco/TrkReco-00-08-59-patch4-slc6tag/src/TMDC.cxx

Go to the documentation of this file.

//-----------------------------------------------------------------------------
// $Id: TMDC.cxx,v 1.16 2012/05/28 05:16:29 maoh Exp $
//-----------------------------------------------------------------------------
// Filename : TMDC.cc
// Section  : Tracking MDC
// Owner    : Yoshi Iwasaki
// Email    : yoshihito.iwasaki@kek.jp
//-----------------------------------------------------------------------------
// Description : A class to represent MDC.
//               See http://bsunsrv1.kek.jp/~yiwasaki/tracking/
//-----------------------------------------------------------------------------

#ifdef TRKRECO_DEBUG_DETAIL
#ifndef TRKRECO_DEBUG
#define TRKRECO_DEBUG
#endif
#endif
#include <iostream>
#include "CLHEP/String/Strings.h"
//#include "belle.h"
#include "TrkReco/TMDC.h"
#include "TrkReco/TMDCUtil.h"
#include "TrkReco/TMDCWire.h"
#include "TrkReco/TMDCLayer.h"
#include "TrkReco/TMDCWireHit.h"
#include "TrkReco/TMDCWireHitMC.h"
#include "TrkReco/TMLink.h"
#include "TrkReco/TTrack.h"
#include "TrkReco/TTrackHEP.h"
//#include MDC_H
//#include BELLETDF_H
//#include "MdcRecGeo/MdcRecGeo.h"
#include "MdcGeomSvc/MdcGeomSvc.h"
#include "MdcTables/MdcTables.h"
//#include "tables/bestdf.h"

//#include "GaudiKernel/Algorithm.h"
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/AlgFactory.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/IMessageSvc.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/SmartDataPtr.h"
#include "GaudiKernel/PropertyMgr.h"
#include "GaudiKernel/IJobOptionsSvc.h"
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/StatusCode.h"
#include "EventModel/Event.h"
//#include "MDCRawEvent/MdcDigi.h"
//#include "TOFRawEvent/TofDigi.h"
//#include "EMCRawEvent/CalDigi.h"
//#include "McTruth/McKine.h"
 
#include <vector>
#include <iostream>
 
using namespace std;
using namespace Event;

TMDC * GMDC = 0;
/*zsl
extern "C" {
    void calcdc_driftdist_(int *,
                           int *,
                           int *,
                           float[3],
                           float[3],
                           float *,
                           float *,
                           float *);
    void calcdc_tof2_(int *, float *, float *, float *);
};
*/

std::string
TMDC::name(void) const {
    return "TMDC";
}

std::string
TMDC::version(void) const {
    return "2.18";
}

TMDC *
TMDC::_cdc = 0;

TMDC *
TMDC::getTMDC(const std::string & version) {
    if (! _cdc) _cdc = new TMDC(version);
    return _cdc;
}

TMDC *
TMDC::getTMDC(void) {
    if (! _cdc) _cdc = new TMDC("1.0");
    return _cdc;
}

TMDC::TMDC(const std::string & version) 
: _debugLevel(0),
  _fudgeFactor(1.),
  _cdcVersion(version), 
//  _cdcVersion(-2000),
//  _nWires(BsCouTab(GEOMDC_WIRE)),
//  _nSuperLayers(BsCouTab(GEOMDC_SUPER)),
//  _nLayers(BsCouTab(GEOMDC_LAYER)),
//  _geo(MdcRecGeo::getMdcRecGeo()),
//  _nWires(_geo->getWireSize()),
//  _nSuperLayers(_geo->getSuperLayerSize()),
//  _nLayers(_geo->getLayerSize()),
  _nWires(6796),
  _nSuperLayers(11),
  _nLayers(43),
  _newCdc(version == "-2000") {
    
    GMDC = this;

    //...Get MDC information...
    std::cout << "TMDC ... MDC version = " << _cdcVersion << std::endl;

    //...Check GEOMDC_WIRE...
/*    if (BsCouTab(GEOMDC_LAYER) == 0) {
        std::cout << "TMDC !!! Bank GEOMDC_LAYER not found" << std::endl;
        std::cout << "     !!! mission impossible" << std::endl;
        std::cout << "     !!! TrkReco will crash!" << std::endl;
        std::cout << "     !!! Please stop execution" << std::endl;
        return;
    }
*/

  IMdcGeomSvc* mdcGeomSvc;
  StatusCode sc = Gaudi::svcLocator()->service("MdcGeomSvc", mdcGeomSvc);
  if (sc ==  StatusCode::SUCCESS) {
//   mdcGeomSvc->Dump();
    
//   _nWires = mdcGeomSvc->getWireSize();
//   _nSuperLayers = mdcGeomSvc->getSuperLayerSize();
//   _nLayers = mdcGeomSvc->getLayerSize();
   
    //...Create super-layers...
    for (unsigned i = 0; i < _nSuperLayers; i++)
        _superLayers.append(new AList<TMDCLayer>);

    //...Loop over all wires layers...
    unsigned nWires = 0;
    for (unsigned i = 0; i < _nLayers; i++) {
//      struct geocdc_layer * l =
//          (struct geocdc_layer *)
//          BsGetEnt(GEOMDC_LAYER, i + 1, BBS_No_Index);
            
        const MdcGeoLayer * l = mdcGeomSvc->Layer(i);
        TMDCLayer * layer =  new TMDCLayer(l);
        _layers.append(layer);
        _superLayers[layer->superLayerId()]->append(layer);

        //...Loop over all wires in a layer...
        for (unsigned j = 0; j < l->NCell(); j++) {
//      for (unsigned j = 0; j < l->m_div; j++) {
//          struct geocdc_wire * w =
//              (struct geocdc_wire *)
//              BsGetEnt(GEOMDC_WIRE, nWires + 1, BBS_No_Index);
            const MdcGeoWire * w = mdcGeomSvc->Wire(l->Wirst()+j);
            TMDCWire * tw = new TMDCWire(w, layer);
            _wires.append(tw);
            layer->append(tw);

            ++nWires;
        }
    }

  } else {
//    return StatusCode::FAILURE;
    return;
  }
 
//  return StatusCode::SUCCESS;
}

void
TMDC::dump(const std::string & msg) const {
//      if (   msg.index("name")    != -1
//      || msg.index("version") != -1
//      || msg.index("detail")  != -1
//      || msg == "") {
  if (   msg.find("name")    != std::string::npos
        || msg.find("version") != std::string::npos
        || msg.find("detail")  != std::string::npos
        || msg == "") {
        std::cout << name() << "(" << version() << ") ";
    }
    if (msg.find("detail") != std::string::npos || msg.find("state") != std::string::npos) {
        std::cout << "Debug Level=" << _debugLevel;
    }
    std::cout << std::endl;

    std::string tab("        ");

    if (msg == "" || msg.find("geometry") != std::string::npos) {

        //...Get information..."
        unsigned nLayer = _layers.length();
        std::cout << "    version    : " << version() << std::endl;
        std::cout << "    cdc version: " << cdcVersion() << std::endl;
        std::cout << "    # of wires : " << _wires.length() << std::endl;
        std::cout << "    # of layers: " << nLayer << std::endl;
        std::cout << "    super layer information" << std::endl;
        std::cout << "        # of super layers = " << _superLayers.length()
                  << std::endl;

        std::cout << "    layer information" << std::endl;
        for (unsigned i = 0; i < _nLayers; i++)
            _layers[i]->dump("", tab);

        std::cout << "    wire information" << std::endl;
        for (unsigned i = 0; i < _nWires; i++)
            (_wires[i])->dump("neighbor", tab);

        return;
    }
    if (msg.find("hits") != std::string::npos) {
        std::cout << "    hits     : " << _hits.length() << std::endl;
        for (unsigned i = 0; i < _hits.length(); i++) {
            _hits[i]->dump("state mc", tab);
        }
    }
    if (msg.find("axialHits") != std::string::npos) {
        std::cout << "    hits     : " << _axialHits.length() << std::endl;
        for (unsigned i = 0; i < _axialHits.length(); i++) {
            _axialHits[i]->dump("state mc", tab);
        }
    }
    if (msg.find("stereoHits") != std::string::npos) {
        std::cout << "    hits     : " << _stereoHits.length() << std::endl;
        for (unsigned i = 0; i < _stereoHits.length(); i++) {
            _stereoHits[i]->dump("state mc", tab);
        }
    }
}

const TMDCWire * const
TMDC::wire(unsigned layerId, int localId) const {
    if (layerId < _nLayers)
        return _layers[layerId]->wire(localId);

    return NULL;
}

const TMDCWire * const
TMDC::wire(const HepPoint3D & p) const {
    float r = p.mag();
    float phi = p.phi();
    return wire(r, phi);
}

const TMDCWire * const
TMDC::wire(float r, float p) const {

    std::cout << "r,phi = " << r << "," << p << std::endl;

    unsigned id = 25;
    bool ok = false;
    const TMDCLayer * l;
    while (! ok) {
        l = layer(id);
        if (! l) return NULL;
        
        const MdcGeoLayer * geo = l->geocdc();
        if (geo->Radius()/10 + geo->RCSiz2()/10 < r) ++id;
        else if (geo->Radius()/10 - geo->RCSiz1()/10 > r) --id;
        else ok = true;
    }
//    float dPhi = 2. * M_PI / float(l->nWires());
//    if (l->geocdc()->m_offset > 0.) p -= dPhi / 2.;
    float dPhi = 2. * M_PI / float(l->nWires());   
    p -= l->geocdc()->Offset()/dPhi; 
    unsigned localId = unsigned(phi(p) / dPhi);
    return l->wire(localId);
}

void
TMDC::clear(void) {
    unsigned i = 0;
    while (TMDCWire * w = _wires[i++])
        w->clear();

    _hitWires.removeAll();
    _axialHits.removeAll();
    _stereoHits.removeAll();
    HepAListDeleteAll(_hits);
    HepAListDeleteAll(_hitsMC);
    HepAListDeleteAll(_badHits);

    TUpdater::clear();
}

void
TMDC::fastClear(void) {
    unsigned i = 0;
    while (TMDCWire * w = _hitWires[i++])
        w->clear();
    i = 0;
    while (TMDCWireHit * h = _badHits[i++])
        ((TMDCWire *) h->wire())->clear();

    _hitWires.removeAll();
    _axialHits.removeAll();
    _stereoHits.removeAll();
    HepAListDeleteAll(_hits);
    HepAListDeleteAll(_hitsMC);
    HepAListDeleteAll(_badHits);

    TUpdater::clear();
}

void
TMDC::update(bool mcAnalysis) {
    //...Already updated?...
//zsl    if (updated()) return;

    //...Clear old information...
    fastClear();

    //...Loop over RECMDC_WIRHIT bank...
//    unsigned nReccdc = BsCouTab(RECMDC_WIRHIT);
     unsigned nReccdc = MdcRecWirhitCol::getMdcRecWirhitCol()->size();
#ifdef TRKRECO_DEBUG
    cout<<"size of MdcRecWirhit:"<<nReccdc<<endl;   
#endif
    for (unsigned i = 0; i < nReccdc; i++) {
//      struct reccdc_wirhit * h =
//          (struct reccdc_wirhit *)
//          BsGetEnt(RECMDC_WIRHIT, i + 1, BBS_No_Index);
        MdcRec_wirhit* h = &(*MdcRecWirhitCol::getMdcRecWirhitCol())[i]; 
        //...Check validity...
//zsl   if (! (h->stat & WireHitFindingValid)) continue;

        //...Obtain a pointer to GEOMDC...
//      struct geocdc_wire * g =
//          (struct geocdc_wire *)
//          BsGetEnt(GEOMDC_WIRE, h->m_geo, BBS_No_Index);
        const MdcGeoWire* g = h->geo;
//      if (_newCdc)
//          if (g->m_ID < 256)
//              continue;

        //...Get a pointer to a TMDCWire...
//      TMDCWire * w = _wires[g->m_ID - 1];
        TMDCWire * w = _wires[g->Id()];   
        _hitWires.append(w);

        //...Create TMDCWireHit...
        TMDCWireHit * hit = new TMDCWireHit(w, h, _fudgeFactor);
        _hits.append(hit);

        //... update the wirehit pointer in the TMDCWire  (by zang shilei)
        w->hit(hit);
       //if (w->hit()) cout<<"test .......... right~~~~~"<<endl; 
        
        //...Axial or stereo...
        if (w->axial()) _axialHits.append(hit);
        else            _stereoHits.append(hit);
    }
            
    //...Hit classfication...
    classification();

    //...MC information...
    if (mcAnalysis) updateMC();

    //...Update information...
    TUpdater::update();
}

void
TMDC::updateMC(void) {

    //...Create TTrackHEP...
    TTrackHEP::update();

    //...Loop over DATMDC_MCWIRHIT bank...
    unsigned n = 0;
//    for (unsigned i = 0; i < BsCouTab(DATMDC_MCWIRHIT); i++) {
    for (unsigned i = 0; i < MdcDatMcwirhitCol::getMdcDatMcwirhitCol()->size(); i++) {   
//        struct datcdc_mcwirhit * h =
//          (struct datcdc_mcwirhit *)
//          BsGetEnt(DATMDC_MCWIRHIT, i + 1, BBS_No_Index);
        MdcDat_mcwirhit* h = &(*MdcDatMcwirhitCol::getMdcDatMcwirhitCol())[i];
        //...Get a pointer to RECMDC_WIRHIT...
//      reccdc_wirhit * whp =
//          (reccdc_wirhit *) BsGetEnt(RECMDC_WIRHIT, h->m_dat, BBS_No_Index);
        MdcRec_wirhit* whp = h->dat->rec;
        //...Get TrkReco objects...
        TMDCWireHit * wh = 0;
        TMDCWire * w = 0;
        if (whp) {
            if (whp->stat & WireHitFindingValid) {
                unsigned n = _hits.length();
                unsigned j = (whp->id < n) ? whp->id : n;
                while (j) {
                    --j;
                    if (_hits[j]->reccdc() == whp) {
                        wh = _hits[j];
                        w = _wires[wh->wire()->id()];
                        break;
                    }
                }
            }
        }
        if (! w) {
//          geocdc_wire * g =
//              (geocdc_wire *) BsGetEnt(GEOMDC_WIRE, h->m_geo, BBS_No_Index);
//          w = _wires[g->m_ID - 1];
          const MdcGeoWire* g = h->geo;
        }

        //...Create TMDCWireHitMC...
        TMDCWireHitMC * hit = new TMDCWireHitMC(w, wh, h);
        _hitsMC.append(hit);
        w->hit(hit);
        if (wh) wh->mc(hit);

        //...TTrackHEP...
//      TTrackHEP * hep = TTrackHEP::list()[h->m_hep - 1];
        TTrackHEP * hep = TTrackHEP::list()[h->hep->id];     
        hit->_hep = hep;
        if (hep) hep->_hits.append(hit);
        else {
            std::cout << "TMDC::updateMC !!! mission impossible" << std::endl;
            std::cout << "                   This error will cause TrkReco crush";
            std::cout << std::endl;
#ifdef TRKRECO_DEBUG_DETAIL
            //cout << "    h->m_hep, h->m_hep -1 = " << h->m_hep;
            //std::cout << ", " << h->m_hep - 1 << std::endl;
            std::cout << "    TTrackHEP list length = ";
            std::cout << TTrackHEP::list().length() << std::endl;
//          BsShwDat(GEN_HEPEVT);
//          BsShwDat(DATMDC_MCWIRHIT);
#endif
        }
    }
}

void
TMDC::classification(void) {
    unsigned n = _hits.length();

    for (unsigned i = 0; i < n; i++) {
        TMDCWireHit * h = _hits[i];
        const TMDCWire * w = h->wire();
        unsigned state = h->state();

        //...Cache pointers to a neighbor...
        const TMDCWire * neighbor[6];
        for (unsigned j = 0; j < 6; j++) neighbor[j] = w->neighbor(j);

        //output to test ...
//      cout<<w->localId()<<endl;
//      if(neighbor[0]) cout<<"inner local:"<<neighbor[0]->localId()<<", "<<neighbor[1]->localId()<<endl;
//      if(neighbor[5]) cout<<"outer local:"<<neighbor[4]->localId()<<", "<<neighbor[5]->localId()<<endl;   
                
        //...Decide hit pattern...
        unsigned pattern = 0;
        for (unsigned j = 0; j < 6; j++) {
            if (neighbor[j])
                if (neighbor[j]->hit())
                    pattern += (1 << j);
        }
        state |= (pattern << WireHitNeighborHit);

        //...Check isolation...
        const TMDCWireHit * hr1 = neighbor[2]->hit();
        const TMDCWireHit * hl1 = neighbor[3]->hit();
        if ((hr1 == 0) && (hl1 == 0)) {
            state |= WireHitIsolated;
        }
        else {
            const TMDCWireHit * hr2 = neighbor[2]->neighbor(2)->hit();
            const TMDCWireHit * hl2 = neighbor[3]->neighbor(3)->hit();
            if ((hr2 == 0) && (hr1 != 0) && (hl1 == 0) ||
                (hl2 == 0) && (hl1 != 0) && (hr1 == 0))
                state |= WireHitIsolated;
        }

        //...Check continuation...
        unsigned superLayer = w->superLayerId();
        bool previous = false;
        bool next = false;
        if (neighbor[0] == 0) previous = true;
        else {
            if ((neighbor[0]->hit()) || neighbor[1]->hit())
                previous = true;
        }
        if (neighbor[5] == 0) next = true;
        else {
            if ((neighbor[4]->hit()) || neighbor[5]->hit())
                next = true;
        }
        // if (previous && next) state |= WireHitContinuous;
        if (previous || next) state |= WireHitContinuous;

        //...Solve LR locally...
        if ((pattern == 34) || (pattern == 42) ||
            (pattern == 40) || (pattern == 10) ||
            (pattern == 35) || (pattern == 50))
            state |= WireHitPatternRight;
        else if ((pattern == 17) || (pattern == 21) ||
                 (pattern == 20) || (pattern ==  5) ||
                 (pattern == 19) || (pattern == 49))
            state |= WireHitPatternLeft;

        //...Store it...
        h->state(state);
    }
}

const AList<TMDCWireHit> &
TMDC::axialHits(unsigned mask) const {
    if (! mask) return _axialHits;
    else if (mask == WireHitFindingValid) return _axialHits;
    std::cout << "TMDC::axialHits !!! unsupported mask given" << std::endl;
    return _axialHits;
}

const AList<TMDCWireHit> &
TMDC::stereoHits(unsigned mask) const {
    if (! mask) return _stereoHits;
    else if (mask == WireHitFindingValid) return _stereoHits;
    std::cout << "TMDC::stereoHits !!! unsupported mask given" << std::endl;
    return _stereoHits;
}

const AList<TMDCWireHit> &
TMDC::hits(unsigned mask) const {
    if (! mask) return _hits;
    else if (mask == WireHitFindingValid) return _hits;
    std::cout << "TMDC::hits !!! unsupported mask given" << std::endl;
    return _hits;
}

const AList<TMDCWireHit> &
TMDC::badHits(void) {
    if (! updated()) update();
    if (_badHits.length()) return _badHits;

    //...Loop over RECMDC_WIRHIT bank...
//    unsigned nReccdc = BsCouTab(RECMDC_WIRHIT);
    unsigned nReccdc = MdcRecWirhitCol::getMdcRecWirhitCol()->size();
    for (unsigned i = 0; i < nReccdc; i++) {
//      struct reccdc_wirhit * h =
//          (struct reccdc_wirhit *)
//          BsGetEnt(RECMDC_WIRHIT, i + 1, BBS_No_Index);
        MdcRec_wirhit* h = &(*MdcRecWirhitCol::getMdcRecWirhitCol())[i];
        
        //...Check validity...
        if (h->stat & WireHitFindingValid) continue;

        //...Obtain a pointer to GEOMDC...
//      geocdc_wire * g =
//          (geocdc_wire *) BsGetEnt(GEOMDC_WIRE, h->m_geo, BBS_No_Index);
        const MdcGeoWire* g = h->geo;
        
        //...Get a pointer to a TMDCWire...
//      TMDCWire * w = _wires[g->m_ID - 1];
        TMDCWire * w = _wires[g->Id()];
        
        //...Create TMDCWireHit...
        _badHits.append(new TMDCWireHit(w, h, _fudgeFactor));
    }

    return _badHits;
}

std::string
TMDC::wireName(unsigned wireId) {
    if (superLayerId(wireId) % 2)
        return itostring(layerId(wireId)) + "=" +
            itostring(localId(wireId));
    return itostring(layerId(wireId)) + "-" + itostring(localId(wireId));
}

std::string
TMDC::wireName(const MdcGeoWire * const w) {
    if (! w) return std::string("no such a wire");
//    unsigned id = w->m_ID - 1;
    unsigned id = w->Id();
    return wireName(id);
}

//std::string
//TMDC::wireName(const Geocdc_wire * const w) {
//    if (! w) return std::string("no such a wire");
//    unsigned id = w->get_ID() - 1;
//    return wireName(id);
//}

unsigned
TMDC::layerId(unsigned id) {
//    int length = _wires.length();
//    if(id < length) return _wires[id]->geocdc()->Layer();
//    return 9999;    

    if (id < 40) return 0;
    else if (id < 84) return 1;
    else if (id < 132) return 2;
    else if (id < 188) return 3;

    else if (id < 252) return 4;
    else if (id < 324) return 5;
    else if (id < 404) return 6;
    else if (id < 484) return 7;

    else if (id < 560) return 8;
    else if (id < 636) return 9;
    else if (id < 724) return 10;
    else if (id < 812) return 11;

    else if (id < 912) return 12;
    else if (id < 1012) return 13;
    else if (id < 1124) return 14;
    else if (id < 1236) return 15;

    else if (id < 1364) return 16;
    else if (id < 1492) return 17;
    else if (id < 1632) return 18;
    else if (id < 1772) return 19;

    else if (id < 1932) return 20;
    else if (id < 2092) return 21;
    else if (id < 2252) return 22;
    else if (id < 2412) return 23;

    else if (id < 2604) return 24;
    else if (id < 2796) return 25;
    else if (id < 2988) return 26;
    else if (id < 3180) return 27;

    else if (id < 3388) return 28;
    else if (id < 3596) return 29;
    else if (id < 3804) return 30;
    else if (id < 4012) return 31;

    else if (id < 4252) return 32;
    else if (id < 4492) return 33;
    else if (id < 4732) return 34;
    else if (id < 4972) return 35;

    else if (id < 5228) return 36;
    else if (id < 5484) return 37;
    else if (id < 5740) return 38;
    else if (id < 5996) return 39;

    else if (id < 6284) return 40;
    else if (id < 6572) return 41;
    else if (id < 6860) return 42;

    return 9999;    
}

unsigned
TMDC::layerId(const MdcGeoWire * const w) {
    if (! w) return 9999;
    //zsl    unsigned id = w->Id();
    //    return layerId(id);
    return w->Layer();
}

//unsigned
//TMDC::layerId(const Geocdc_wire * const w) {
//    if (! w) return 9999;
//    unsigned id = w->get_ID() - 1;
//    return layerId(id);
//}

unsigned
TMDC::localId(unsigned id) {
    if (id < 40) return id;
    else if (id < 84) return id-40;
    else if (id < 132) return id-84;
    else if (id < 188) return id-132;

    else if (id < 252) return id-188;
    else if (id < 324) return id-252;
    else if (id < 404) return id-324;
    else if (id < 484) return id-404;

    else if (id < 560) return id-484;
    else if (id < 636) return id-560;
    else if (id < 724) return id-636;
    else if (id < 812) return id-724;

    else if (id < 912) return id-812;
    else if (id < 1012) return id-912;
    else if (id < 1124) return id-1012;
    else if (id < 1236) return id-1124;

    else if (id < 1364) return id-1236;
    else if (id < 1492) return id-1364;
    else if (id < 1632) return id-1492;
    else if (id < 1772) return id-1632;

    else if (id < 1932) return id-1772;
    else if (id < 2092) return id-1932;
    else if (id < 2252) return id-2092;
    else if (id < 2412) return id-2252;

    else if (id < 2604) return id-2412;
    else if (id < 2796) return id-2604;
    else if (id < 2988) return id-2796;
    else if (id < 3180) return id-2988;

    else if (id < 3388) return id-3180;
    else if (id < 3596) return id-3388;
    else if (id < 3804) return id-3596;
    else if (id < 4012) return id-3804;

    else if (id < 4252) return id-4012;
    else if (id < 4492) return id-4252;
    else if (id < 4732) return id-4492;
    else if (id < 4972) return id-4732;

    else if (id < 5228) return id-4972;
    else if (id < 5484) return id-5228;
    else if (id < 5740) return id-5484;
    else if (id < 5996) return id-5740;

    else if (id < 6284) return id-5996;
    else if (id < 6572) return id-6284;
    else if (id < 6860) return id-6572;

    return 9999;    
 
/*    if (id < 384) return id % 64;
    else if (id < 624) return (id - 384) % 80;
    else if (id < 1200) return (id - 624) % 96;
    else if (id < 1584) return (id - 1200) % 128;
    else if (id < 2304) return (id - 1584) % 144;
    else if (id < 2944) return (id - 2304) % 160;
    else if (id < 3904) return (id - 2944) % 176;
    else if (id < 4736) return (id - 3904) % 208;
    else if (id < 5936) return (id - 4736) % 240;
    else if (id < 6960) return (id - 5936) % 256;
    else if (id < 8400) return (id - 6960) % 288;

    return 9999;
*/
}

unsigned
TMDC::localId(const MdcGeoWire * const w) {
    if (! w) return 9999;
//    unsigned id = w->Id();
//    return localId(id);
    return w->Cell();    
}

//unsigned
//TMDC::localId(const Geocdc_wire * const w) {
//    if (! w) return 9999;
//    unsigned id = w->get_ID() - 1;
//    return localId(id);
//}

unsigned
TMDC::superLayerId(unsigned id) {
    if (id < 188) return 0;
    else if (id < 484) return 1;
    else if (id < 812) return 2;
    else if (id < 1236) return 3;
    else if (id < 1772) return 4;
    else if (id < 2412) return 5;
    else if (id < 3180) return 6;
    else if (id < 4012) return 7;
    else if (id < 4972) return 8;
    else if (id < 5996) return 9;
    else if (id < 6860) return 10;
 
    return 9999;    
/*    if (id < 384) return 0;
    else if (id < 624) return 1;
    else if (id < 1200) return 2;
    else if (id < 1584) return 3;
    else if (id < 2304) return 4;
    else if (id < 2944) return 5;
    else if (id < 3904) return 6;
    else if (id < 4736) return 7;
    else if (id < 5936) return 8;
    else if (id < 6960) return 9;
    else if (id < 8400) return 10;

    return 9999;
*/
}

unsigned
TMDC::superLayerId(const MdcGeoWire * const w) {
    if (! w) return 9999;
    unsigned id = w->Id();
    return superLayerId(id);
}

//unsigned
//TMDC::superLayerId(const Geocdc_wire * const w) {
//    if (! w) return 9999;
//    unsigned id = w->get_ID() - 1;
//    return superLayerId(id);
//}

unsigned
TMDC::superLayerId(const MdcGeoLayer * const l) {
    if (! l) return 9999;
    unsigned id = l->Id();

    if(id < 44) return id/4;
    return 9999;     
/*    if (id < 6) return 0;
    else if (id < 9) return 1;
    else if (id < 15) return 2;
    else if (id < 18) return 3;
    else if (id < 23) return 4;
    else if (id < 27) return 5;
    else if (id < 32) return 6;
    else if (id < 36) return 7;
    else if (id < 41) return 8;
    else if (id < 45) return 9;
    else if (id < 50) return 10;

    return 9999;
*/
}

unsigned
TMDC::localLayerId(unsigned id) {
    if(id < 6860) return layerId(id)%4;
    return 9999;
    
/*    if (id < 64) return 0;
    else if (id < 128) return 1;
    else if (id < 192) return 2;
    else if (id < 256) return 3;
    else if (id < 320) return 4;
    else if (id < 384) return 5;

    else if (id < 464) return 0;
    else if (id < 544) return 1;
    else if (id < 624) return 2;

    else if (id < 720) return 0;
    else if (id < 816) return 1;
    else if (id < 912) return 2;
    else if (id < 1008) return 3;
    else if (id < 1104) return 4;
    else if (id < 1200) return 5;

    else if (id < 1328) return 0;
    else if (id < 1456) return 1;
    else if (id < 1584) return 2;

    else if (id < 1728) return 0;
    else if (id < 1872) return 1;
    else if (id < 2016) return 2;
    else if (id < 2160) return 3;
    else if (id < 2304) return 4;

    else if (id < 2464) return 0;
    else if (id < 2624) return 1;
    else if (id < 2784) return 2;
    else if (id < 2944) return 3;

    else if (id < 3136) return 0;
    else if (id < 3328) return 1;
    else if (id < 3520) return 2;
    else if (id < 3712) return 3;
    else if (id < 3904) return 4;

    else if (id < 4112) return 0;
    else if (id < 4320) return 1;
    else if (id < 4528) return 2;
    else if (id < 4736) return 3;

    else if (id < 4976) return 0;
    else if (id < 5216) return 1;
    else if (id < 5456) return 2;
    else if (id < 5696) return 3;
    else if (id < 5936) return 4;

    else if (id < 6192) return 0;
    else if (id < 6448) return 1;
    else if (id < 6704) return 2;
    else if (id < 6960) return 3;

    else if (id < 7248) return 0;
    else if (id < 7536) return 1;
    else if (id < 7824) return 2;
    else if (id < 8112) return 3;
    else if (id < 8400) return 4;

    return 9999;    
*/
}

unsigned
TMDC::localLayerId(const MdcGeoWire * const w) {
    if (! w) return 9999;
    return w->Layer()%4;
    
//    unsigned id = w->Id();
//    return localLayerId(id);
}

//unsigned
//TMDC::localLayerId(const Geocdc_wire * const w) {
//    if (! w) return 9999;
//    unsigned id = w->get_ID() - 1;
//    return localLayerId(id);
//}

unsigned
TMDC::localLayerId(const MdcGeoLayer * const l) {
    if (! l) return 9999;
    unsigned id = l->Id();
    if (id < 44) return id%4;      
    return 9999;
    
/*    if (id < 1) return 0;
    else if (id < 2) return 1;
    else if (id < 3) return 2;
    else if (id < 4) return 3;
    else if (id < 5) return 4;
    else if (id < 6) return 5;

    else if (id < 7) return 0;
    else if (id < 8) return 1;
    else if (id < 9) return 2;

    else if (id < 10) return 0;
    else if (id < 11) return 1;
    else if (id < 12) return 2;
    else if (id < 13) return 3;
    else if (id < 14) return 4;
    else if (id < 15) return 5;

    else if (id < 16) return 0;
    else if (id < 17) return 1;
    else if (id < 18) return 2;

    else if (id < 19) return 0;
    else if (id < 20) return 1;
    else if (id < 21) return 2;
    else if (id < 22) return 3;
    else if (id < 23) return 4;

    else if (id < 24) return 0;
    else if (id < 25) return 1;
    else if (id < 26) return 2;
    else if (id < 27) return 3;

    else if (id < 28) return 0;
    else if (id < 29) return 1;
    else if (id < 30) return 2;
    else if (id < 31) return 3;
    else if (id < 32) return 4;

    else if (id < 33) return 0;
    else if (id < 34) return 1;
    else if (id < 35) return 2;
    else if (id < 36) return 3;

    else if (id < 37) return 0;
    else if (id < 38) return 1;
    else if (id < 39) return 2;
    else if (id < 40) return 3;
    else if (id < 41) return 4;

    else if (id < 42) return 0;
    else if (id < 43) return 1;
    else if (id < 44) return 2;
    else if (id < 45) return 3;

    else if (id < 46) return 0;
    else if (id < 47) return 1;
    else if (id < 48) return 2;
    else if (id < 49) return 3;
    else if (id < 50) return 4;

    return 9999;    
*/
}

unsigned
TMDC::axialStereoLayerId(const MdcGeoLayer * const l) {
    if (! l) return 9999;
    unsigned id = l->Id();

    if (id < 8) return id;
    else if (id < 20) return id-8;
    else if (id < 36) return id-12;
    else if (id < 43) return id-24;

    return 9999;    
}

unsigned
TMDC::layerId(unsigned as, unsigned id) {

    //...Axial case...
    if (as == 0) {
        if (id < 12) return id + 8;
        else if (id < 19) return id + 24;
        return 9999;
    }

    //...Stereo case...
    if (id < 8) return id;
    else if (id < 24) return id + 12;
    return 9999;
}

std::string
TMDC::wireName(const MdcRec_wirhit & h) {
//    geocdc_wire * g = (geocdc_wire *) BsGetEnt(GEOMDC_WIRE,
//                                             h.m_geo,
//                                             BBS_No_Index);
      const MdcGeoWire* g = h.geo;
      return wireName(g);
}

void
TMDC::driftDistance(TMLink & l,
                    const TTrack & t,
                    unsigned flag,
                    float t0Offset) {

    //...No correction...
    if (flag == 0) {
        if (l.hit()) {
            l.drift(0, l.hit()->drift(0));
            l.drift(1, l.hit()->drift(1));
            l.dDrift(0, l.hit()->dDrift(0));
            l.dDrift(1, l.hit()->dDrift(1));
        }
        else {
            l.drift(0, 0.);
            l.drift(1, 0.);
            l.dDrift(0, 0.);
            l.dDrift(1, 0.);
        }

        return;
    }

    //...TOF correction...
    float tof = 0.;
    if (flag && 1) {
        int imass = 3;
        float tl = t.helix().a()[4];
        float f = sqrt(1. + tl * tl);
        float s = fabs(t.helix().curv()) * fabs(l.dPhi()) * f;
        float p = f / fabs(t.helix().a()[2]);
//zsl   calcdc_tof2_(& imass, & p, & s, & tof); 
    }

    //...T0 correction....
    if (! (flag && 2)) t0Offset = 0.;

    //...Propagation corrections...
    const TMDCWireHit & h = * l.hit();
    int wire = h.wire()->id();
    HepVector3D tp = t.helix().momentum(l.dPhi());
    float p[3] = {tp.x(), tp.y(), tp.z()};
    const HepPoint3D & onWire = l.positionOnWire();
    float x[3] = {onWire.x(), onWire.y(), onWire.z()};
    float time = h.reccdc()->tdc + t0Offset - tof;
    float dist;
    float edist;
    int prop = (flag & 4);

    //...Calculation with left side...
    int side = -1;
    if (side == 0) side = -1;
//zsl    calcdc_driftdist_(& prop,
//                    & wire,
//                    & side,
//                    p,
//                    x,
//                    & time,
//                    & dist,
//                    & edist);
    l.drift(0, dist);
    l.dDrift(0, edist);

    //...Calculation with left side...
    side = 1;
//zsl    calcdc_driftdist_(& prop,
//                    & wire,
//                    & side,
//                    p,
//                    x,
//                    & time,
//                    & dist,
//                    & edist);
    l.drift(1, dist);
    l.dDrift(1, edist);

    //...tan(lambda) correction...
    if (flag && 8) {
        float tanl = abs(p[2]) / tp.perp();
        float c;
        if ((tanl >= 0.0) && (tanl < 0.5))      c = -0.48 * tanl + 1.3;
        else if ((tanl >= 0.5) && (tanl < 1.0)) c = -0.28 * tanl + 1.2;
        else if ((tanl >= 1.0) && (tanl < 1.5)) c = -0.16 * tanl + 1.08;
        else                                    c =  0.84;

        l.dDrift(0, l.dDrift(0) * c);
        l.dDrift(1, l.dDrift(1) * c);
    }
}

---