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

Go to the documentation of this file.

//-----------------------------------------------------------------------------
// $Id: TMDCWire.cxx,v 1.14 2011/08/02 12:31:31 liucy Exp $
//-----------------------------------------------------------------------------
// Filename : TMDCWire.cc
// Section  : Tracking MDC
// Owner    : Yoshi Iwasaki
// Email    : yoshihito.iwasaki@kek.jp
//-----------------------------------------------------------------------------
// Description : A class to represent a wire in MDC.
//               See http://bsunsrv1.kek.jp/~yiwasaki/tracking/
//-----------------------------------------------------------------------------

#define TMDCWire_INLINE_DEFINE_HERE
#include "CLHEP/String/Strings.h"
//#include "belle.h"
//#include MDC_H
//#include "MdcRecGeo/MdcRecGeo.h"
#include "MdcGeomSvc/MdcGeomSvc.h"
#include "TrkReco/TMDC.h"
#include "TrkReco/TMDCWire.h"
#include "TrkReco/TMDCTsf.h"
#include "TrkReco/TMDCUtil.h"

TMDC *
TMDCWire::_cdc = 0;

/*zsl
extern "C" void
calcdc_sag3_(int *, float *, float[3], float *, float *, float *);
*/

TMDCWire::TMDCWire(const MdcGeoWire * w, TMDCLayer * l)
: _forwardPosition(w->Backward()*0.1), 
  _backwardPosition(w->Forward()*0.1), 
  _tension(w->Tension()),
  _direction(w->Backward() -w->Forward()) { 
//: _forwardPosition(w->Forward()*0.1),
//  _backwardPosition(w->Backward()*0.1),
//  _direction(w->Forward() - w->Backward()) {
    _state = 0;
    _geo = w;
    _id = w->Id();
    _localId = w->Cell();
    _layerId = w->Layer();
    _layer = l;
    _tsf = 0;
    _hit = 0;
    _xyPosition = 0.5 * (_forwardPosition + _backwardPosition);
    _xyPosition.setZ(0.);
    _direction = _direction.unit();
}

TMDCWire::~TMDCWire() {
}

void
TMDCWire::dump(const std::string & msg, const std::string & pre) const {
    std::cout << pre;
    std::cout << "w " << _id;
    std::cout << ",local " << _localId;
    std::cout << ",layer " << _layerId;
    std::cout << ",super layer " << superLayerId();
    std::cout << ",local layer " << localLayerId();
    std::cout << std::endl;
    if (msg.find("neighbor") != std::string::npos || msg.find("detail") != std::string::npos) {
        for (unsigned i = 0; i < 6; i++)
            if (neighbor(i))
                neighbor(i)->dump("", pre + itostring(i) + "   ");
    }    
    if (msg.find("tsf") != std::string::npos || msg.find("detail") != std::string::npos) {
        if (_tsf) _tsf->dump("", pre + "    ");
    }
}

const HepPoint3D 
TMDCWire::xyPosition(double z) const {
    double A=47.35E-6/_tension;
    HepPoint3D x0 =(z - _backwardPosition.z()) /_direction.z() * _direction + _backwardPosition;

    double zinit = x0.z();
    HepVector3D wire;
    wire.setX(_forwardPosition.x()-_backwardPosition.x());
    wire.setY(_forwardPosition.y()-_backwardPosition.y());
    wire.setZ(_forwardPosition.z()-_backwardPosition.z());
    double lxz=sqrt(wire.z()*wire.z()+wire.x()*wire.x());
    double Zp = (zinit - _backwardPosition.z())*lxz/wire.z();
    double x=_direction.x()*(zinit-_backwardPosition.z())/_direction.z() + _backwardPosition.x();
    double y=(A*(Zp-lxz)+wire.y()/lxz)*Zp+_backwardPosition.y();
    x0.setX(x);
    x0.setY(y);
    return x0;
}

const TMDCWire * const
TMDCWire::neighbor(unsigned i) const {
    if (_cdc == 0) _cdc = TMDC::getTMDC();
    static const unsigned nLayers = _cdc->nLayers();

    //...Expensive...
    unsigned l = layerId();
    int local = (int) localId();
    float phi0 = _cdc->layer(l)->offset();
    int nWir = (int) _cdc->layer(l)->nWires();
    float phi = phi0 + local*2*pi/nWir;
    unsigned s = superLayerId();    
//    cout<<"i:"<<i<<"   layerId:"<<l<<"  local:"<<local<<"  nWire:"<<nWir<<"   superLyrId:"<<s<<endl;

    if (i == WireInnerLeft || i == WireInnerRight) {
        if (l == 0) return 0;
        --l;
        if (_cdc->layer(l)->superLayerId() != s) return 0;
        phi0 = _cdc->layer(l)->offset();
        nWir = _cdc->layer(l)->nWires();
        local = (int)((phi-phi0)/(2*pi/nWir)); 

//      cout<<"inner local"<<local<<", "<<local+1<<endl;
        
//      cout<<(phi-phi0)/(2*pi/nWir)<<"  local:"<<local<<endl;
//      cout<<"InnerLeft phi:"<<phi0 + local*2*pi/nWir
//              <<"   InnerRight phi:"<<phi0 + (local+1)*2*pi/nWir 
//              <<"   phi:"<<phi<<endl;
//      if(local<0 || local+1>nWir-1) 
//              cout<<"TMDCWire::neighbor()---Pay attation the cellNo !!!"<<endl
//                  <<l+1<<"th layer: InnerLeftRight,   local:"<<local<<"    local+1:"<<local+1<<endl;  
// LR changed to bes3
        if (i == WireInnerLeft) return _cdc->wire(l, local);
        else                    return _cdc->wire(l, local + 1);
//      if (i == WireInnerRight) return _cdc->wire(l, local);
//      else                     return _cdc->wire(l, local + 1);
        
    }
    else if (i == WireLeft || i == WireRight) {
//        cout<<"Left phi:"<<phi0 + (local-1)*2*pi/nWir  
//            <<"   Right phi:"<<phi0 + (local+1)*2*pi/nWir 
//            <<"   phi:"<<phi<<endl;
//      if(local-1<0 || local+1>nWir-1) 
//              cout<<"TMDCWire::neighbor()---Pay attation the Cell No !!!"<<endl
//                <<l<<"th layer: LeftRight,   local-1:"<<local-1<<"    local+1:"<<local+1<<endl;       
        if (i == WireLeft) return _cdc->wire(l, local - 1);
        else                return _cdc->wire(l, local + 1);
    }
    else if (i == WireOuterLeft || i == WireOuterRight) {
        if (l == nLayers - 1) return 0;
        ++l;
        if (_cdc->layer(l)->superLayerId() != s) return 0;
        phi0 = _cdc->layer(l)->offset();
        nWir = _cdc->layer(l)->nWires();
        local = (int)((phi-phi0)/(2*pi/nWir));
 
//      cout<<"outer local"<<local<<", "<<local+1<<endl;
                
//      cout<<(phi-phi0)/(2*pi/nWir)<<"  local:"<<local<<endl;          
//        cout<<"OuterLeft phi:"<<phi0 + local*2*pi/nWir 
//                <<"   OuterRight phi:"<<phi0 + (local+1)*2*pi/nWir 
//              <<"   phi:"<<phi<<endl;
//        if(local<0 || local+1>nWir-1) 
//              cout<<"TMDCWire::neighbor()---Pay attation the Cell No !!!"<<endl
//               <<l-1<<"th layer:  OuterLeftRight,   local:"<<local<<"    local+1:"<<local+1<<endl;    
        if (i == WireOuterLeft) return _cdc->wire(l, local);
        else                     return _cdc->wire(l, local + 1);
    }
    return 0;
}

int
TMDCWire::localIdForPlus(void) const {
/*
  if(_layerId >= 0 && _layerId <= 5){
    if(_localId == 63)
      return -1;
    else return _localId;
  }

  if(_layerId >= 6 && _layerId <= 8){
    if(_localId == 79)
      return -1;
    else return _localId;
  }

  if(_layerId >= 9 && _layerId <= 14){
    if(_localId == 95)
      return -1;
    else return _localId;
  }
  
  if(_layerId >= 15 && _layerId <= 17){
    if(_localId == 127)
      return -1;
    else return _localId;
  }

  if(_layerId >= 18 && _layerId <= 22){
    if(_localId == 143)
      return -1;
    else return _localId;
  }

  if(_layerId >= 23 && _layerId <= 26){
    if(_localId == 159)
      return -1;
    else return _localId;
  }
  
  if(_layerId >= 27 && _layerId <= 31){
    if(_localId == 175)
      return -1;
    else return _localId;
  }
  
  if(_layerId >= 32 && _layerId <= 35){
    if(_localId == 207)
      return -1;
    else return _localId;
  }

  if(_layerId >= 36 && _layerId <= 40){
    if(_localId == 239)
      return -1;
    else return _localId;
  }

  if(_layerId >= 41 && _layerId <= 44){
    if(_localId == 255)
      return -1;
    else return _localId;
  }

  if(_layerId >= 45){
    if(_localId == 287)
      return -1;
    else return _localId;
  }

  return -1;
*/
//Liuqg 060915
  if(_layerId == 0){
    if(_localId == 39)
      return -1;
    else return _localId;
  }

  if(_layerId == 1){
    if(_localId == 43)
      return -1;
    else return _localId;
  }

  if(_layerId == 2){
    if(_localId == 47)
      return -1;
    else return _localId;
  }

  if(_layerId == 3){
    if(_localId == 55)
      return -1;
    else return _localId;
  }

  if(_layerId == 4){
    if(_localId == 63)
      return -1;
    else return _localId;
  }

  if(_layerId == 5){
    if(_localId == 71)
      return -1;
    else return _localId;
  }

  if(_layerId == 6 || _layerId == 7){
    if(_localId == 79)
      return -1;
    else return _localId;
  }

  if(_layerId == 8 || _layerId == 9){
    if(_localId == 75)
      return -1;
    else return _localId;
  }

  if(_layerId == 10 || _layerId == 11){
    if(_localId == 87)
      return -1;
    else return _localId;
  }

  if(_layerId == 12 || _layerId == 13){
    if(_localId == 99)
      return -1;
    else return _localId;
  }

  if(_layerId == 14 || _layerId == 15){
    if(_localId == 111)
      return -1;
    else return _localId;
  }

  if(_layerId == 16 || _layerId == 17){
    if(_localId == 127)
      return -1;
    else return _localId;
  }

  if(_layerId == 18 || _layerId == 19){
    if(_localId == 139)
      return -1;
    else return _localId;
  }

  if(_layerId >= 20 && _layerId <= 23){
    if(_localId == 159)
      return -1;
    else return _localId;
  }

  if(_layerId >= 24 && _layerId <= 27){
    if(_localId == 175)
      return -1;
    else return _localId;
  }
  
  if(_layerId >= 28 && _layerId <= 31){
    if(_localId == 207)
      return -1;
    else return _localId;
  }

  if(_layerId >= 32 && _layerId <= 35){
    if(_localId == 239)
      return -1;
    else return _localId;
  }

  if(_layerId >= 36 && _layerId <= 39){
    if(_localId == 255)
      return -1;
    else return _localId;
  }
  
  if(_layerId >= 40 && _layerId <= 42){
    if(_localId == 287)
      return -1;
    else return _localId;
  }

  return -1;
}

int
TMDCWire::localIdForMinus(void) const {
/*
  if(_layerId >= 0 && _layerId <= 5){
    if(_localId == 0)
      return 64;
    else return _localId;
  }

  if(_layerId >= 6 && _layerId <= 8){
    if(_localId == 0)
      return 80;
    else return _localId;
  }

  if(_layerId >= 9 && _layerId <= 14){
    if(_localId == 0)
      return 96;
    else return _localId;
  }
  
  if(_layerId >= 15 && _layerId <= 17){
    if(_localId == 0)
      return 128;
    else return _localId;
  }

  if(_layerId >= 18 && _layerId <= 22){
    if(_localId == 0)
      return 144;
    else return _localId;
  }

  if(_layerId >= 23 && _layerId <= 26){
    if(_localId == 0)
      return 160;
    else return _localId;
  }
  
  if(_layerId >= 27 && _layerId <= 31){
    if(_localId == 0)
      return 192;
    else return _localId;
  }
  
  if(_layerId >= 32 && _layerId <= 35){
    if(_localId == 0)
      return 208;
    else return _localId;
  }

  if(_layerId >= 36 && _layerId <= 40){
    if(_localId == 0)
      return 240;
    else return _localId;
  }

  if(_layerId >= 41 && _layerId <= 44){
    if(_localId == 0)
      return 256;
    else return _localId;
  }

  if(_layerId >= 45){
    if(_localId == 0)
      return 288;
    else return _localId;
  }
*/
//Liuqg 060915
  if(_layerId == 0){
    if(_localId == 0)
      return 40;
    else return _localId;
  }

  if(_layerId == 1){
    if(_localId == 0)
      return 44;
    else return _localId;
  }

  if(_layerId == 2){
    if(_localId == 0)
      return 48;
    else return _localId;
  }

  if(_layerId == 3){
    if(_localId == 0)
      return 56;
    else return _localId;
  }

  if(_layerId == 4){
    if(_localId == 0)
      return 64;
    else return _localId;
  }

  if(_layerId == 5){
    if(_localId == 0)
      return 72;
    else return _localId;
  }

  if(_layerId == 6 || _layerId == 7){
    if(_localId == 0)
      return 80;
    else return _localId;
  }

  if(_layerId == 8 || _layerId == 9){
    if(_localId == 0)
      return 76;
    else return _localId;
  }

  if(_layerId == 10 || _layerId == 11){
    if(_localId == 0)
      return 88;
    else return _localId;
  }

  if(_layerId == 12 || _layerId == 13){
    if(_localId == 0)
      return 100;
    else return _localId;
  }

  if(_layerId == 14 || _layerId == 15){
    if(_localId == 0)
      return 112;
    else return _localId;
  }

  if(_layerId == 16 || _layerId == 17){
    if(_localId == 0)
      return 128;
    else return _localId;
  }

  if(_layerId == 18 || _layerId == 19){
    if(_localId == 0)
      return 140;
    else return _localId;
  }

  if(_layerId >= 20 && _layerId <= 23){
    if(_localId == 0)
      return 160;
    else return _localId;
  }

  if(_layerId >= 24 && _layerId <= 27){
    if(_localId == 0)
      return 176;
    else return _localId;
  }

  if(_layerId >= 28 && _layerId <= 31){
    if(_localId == 0)
      return 208;
    else return _localId;
  }

  if(_layerId >= 32 && _layerId <= 35){
    if(_localId == 0)
      return 240;
    else return _localId;
  }

  if(_layerId >= 36 && _layerId <= 39){
    if(_localId == 0)
      return 256;
    else return _localId;
  }

  if(_layerId >= 40 && _layerId <= 42){
    if(_localId == 0)
      return 288;
    else return _localId;
  }
  return -1;
}

void
TMDCWire::wirePosition(float z,
                       HepPoint3D & xy,
                       HepPoint3D & back,
                       HepVector3D & dir) const {
    back = _backwardPosition;

    //...Check z position...
    if (! (z > _backwardPosition.z() && z < _forwardPosition.z())) {
        xy = _xyPosition;
        dir = _direction;
    }

    //...Setup...
    int wireID = id();
    float wirePosition[3] = {0., 0., 0.};
    float dydz = 0;
    float ybSag = 0;
    float yfSag = 0;

    //...Get corrections...
//zsl    calcdc_sag3_(& wireID, & z, wirePosition, & dydz, & ybSag, & yfSag);   
        
    //...Wire position...
    xy.setX((double) wirePosition[0]);
    xy.setY((double) wirePosition[1]);
    xy.setZ((double) wirePosition[2]);
    back.setY((double) ybSag);
    HepVector3D v_aux(_forwardPosition.x() - _backwardPosition.x(),
                   yfSag - ybSag,
                   _forwardPosition.z() - _backwardPosition.z());
    dir = v_aux.unit();

    return;
}

int
TMDCWire::localIdDifference(const TMDCWire & a) const {

#ifdef TRKRECO_DEBUG_DETAIL
    if (superLayerId() != a.superLayerId()) {
        std::cout << "TMDCWire::localIdDifference !!!";
        std::cout << "super layer assumption violation" << std::endl;
    }
#endif

    int diff = int(a.localId()) - int(localId());
    unsigned nWires = layer()->nWires();
    if (diff > 0) {
        int difR = nWires - diff;
        if (diff < difR) return diff;
        else return - difR;
    }
    else {
        int difR = nWires + diff;
        if (- diff < difR) return diff;
        else return difR;
    }
}

---