00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #include "MdcxReco/MdcxHits.h"
00023 #include "CLHEP/Alist/AIterator.h"
00024 #include "MdcData/MdcHit.h"
00025 #include "GaudiKernel/Bootstrap.h"
00026 #include "GaudiKernel/SmartDataPtr.h"
00027 #include "GaudiKernel/IService.h"
00028 #include "GaudiKernel/ISvcLocator.h"
00029 #include "GaudiKernel/IDataProviderSvc.h"
00030 #include "EventModel/Event.h"
00031 #include "Identifier/MdcID.h"
00032 #include "MdcRawEvent/MdcDigi.h"
00033 #include "MdcxReco/MdcxHel.h"
00034 #include "RawEvent/RawDataUtil.h"
00035 #include "MdcGeom/Constants.h"
00036 #include "MdcGeom/EntranceAngle.h"
00037 #include <iomanip>
00038
00039 using std::endl;
00040 using std::ostream;
00041 #ifdef MDCXRECO_RESLAYER
00042 extern int g_resLayer;
00043 #endif
00044 const MdcCalibFunSvc* MdcxHit::m_mdcCalibFunSvc = 0;
00045 bool MdcxHit::m_countPropTime = true;
00046 const MdcDetector* MdcxHit::m_gm= 0;
00047
00048 MdcxHit::MdcxHit(const MdcDigi *pdcdatum, float c0, float cresol)
00049 :_mdcHit(0), _mdcDigi(pdcdatum), _c0(c0), _cresol(cresol)
00050 {
00051 process();
00052 }
00053
00054 MdcxHit::MdcxHit(const MdcHit *pdchhit, float c0, float cresol)
00055 :_mdcHit(pdchhit),_mdcDigi(pdchhit->digi()), _c0(c0), _cresol(cresol)
00056 {
00057 process();
00058 }
00059
00060 void MdcxHit::setMdcCalibFunSvc(const MdcCalibFunSvc* calibSvc) {
00061 m_mdcCalibFunSvc = calibSvc;
00062 }
00063
00064 void MdcxHit::setCountPropTime(bool countPropTime) {
00065 m_countPropTime = countPropTime;
00066 }
00067
00068 void MdcxHit::setMdcDetector(const MdcDetector* gm) {
00069 m_gm = gm;
00070 }
00071
00072 void MdcxHit::process() {
00073 assert(m_gm);
00074 assert(m_mdcCalibFunSvc);
00075 Identifier _id = (*_mdcDigi).identify();
00076 _layernumber=MdcID::layer(_id);
00077 _wirenumber=MdcID::wire(_id);
00078 _superlayer=(_layernumber)/4;
00079 _iTdc = _mdcDigi->getTimeChannel();
00080 _iAdc = _mdcDigi->getChargeChannel();
00081 _t = RawDataUtil::MdcTime(_iTdc);
00082 _q = RawDataUtil::MdcCharge(_iAdc);
00083 _T0Walk = m_mdcCalibFunSvc->getT0(_layernumber,_wirenumber)
00084 + m_mdcCalibFunSvc->getTimeWalk(_layernumber,_iAdc);
00085
00086 const MdcLayer* layerPtr= m_gm->Layer(_layernumber);
00087 _L = layerPtr->zLength();
00088 _x = layerPtr->xWire(_wirenumber);
00089 _y = layerPtr->yWire(_wirenumber);
00090 _r = sqrt(_x*_x + _y*_y);
00091 _s = layerPtr->stereo();
00092 _p = layerPtr->phiOffset() + _wirenumber*layerPtr->dPhi();
00093
00094
00095 double tst = _s;
00096 double pw = atan2(_y, _x);
00097 _pw = pw;
00098 _wx = -tst*sin(pw);
00099 _wy = tst*cos(pw);
00100 _wz = (tst*tst < 1.0) ? sqrt(1.0-tst*tst) : 0.0;
00101 _sp = sin(_p);
00102 _cp = cos(_p);
00103 _consterr = 0;
00104 _e = _cresol;
00105 _d = d();
00106
00107 _v = (_t-_c0 > 0.0) ? _d/(_t-_c0) : 0.0013;
00108 _xneg = _x + _d*_sp;
00109 _yneg = _y - _d*_cp;
00110 _xpos = _x - _d*_sp;
00111 _ypos = _y + _d*_cp;
00112 usedonhel = 0;
00113 _type = Constants::viewOfsLayer[_superlayer];
00114 }
00115
00116 float MdcxHit::tcor(float z, float tof, float tzero)const {
00117
00118 double tprop = 0.;
00119 if (m_countPropTime){ tprop = m_mdcCalibFunSvc->getTprop(_layernumber,z*10); }
00120 double dt = (_t - _T0Walk -_c0 - tprop - tof - tzero);
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132 return dt;
00133 }
00134
00135 float MdcxHit::d(float z, float tof, float tzero, int wamb, float entranceAngle)const {
00136
00137
00138
00139 float t_corr = tcor(z,tof,tzero);
00140 double eAngle = EntranceAngle(entranceAngle);
00141
00142
00143 if (fabs(z)>150. || fabs(t_corr)>1500. || fabs(eAngle)>999){
00144 return 9999.;
00145 }
00146
00147
00148 int lrCalib = 2;
00149 if (wamb==1)lrCalib = 0;
00150 else if (wamb==-1) lrCalib = 1;
00151
00152 double driftD = 0.1 * m_mdcCalibFunSvc->driftTimeToDist(t_corr,_layernumber,_wirenumber,lrCalib,eAngle);
00153
00154
00155
00156 if (fabs(driftD)<=0.0001) driftD = 0.001;
00157 return driftD;
00158 }
00159
00160 float MdcxHit::d(MdcxHel &hel) const {
00161 float doca=hel.Doca(*this);
00162 return d(hel.Doca_Zh(),hel.Doca_Tof(),hel.T0(),
00163 hel.Doca_Wamb(),hel.Doca_Eang());
00164 }
00165
00166 float MdcxHit::e(float dd) const{
00167
00168 float errTemp = fabs(getSigma(dd));
00169
00170 float errMin = 1.0e-7;
00171 return errTemp<errMin?errMin:errTemp;
00172 }
00173
00174 float MdcxHit::pull(MdcxHel &hel) const {
00175
00176 float doca=hel.Doca(*this); if(hel.Mode() == 0)doca=fabs(doca);
00177 return (d(hel.Doca_Zh(),hel.Doca_Tof(),hel.T0(),
00178 hel.Doca_Wamb(),hel.Doca_Eang())-doca)/e(doca);
00179
00180 }
00181
00182 float MdcxHit::residual(MdcxHel &hel) const {
00183
00184 float doca = hel.Doca(_wx, _wy, _wz, _x, _y);
00185 if(hel.Mode() == 0) doca = fabs(doca);
00186 doca += v()*hel.T0();
00187 return d(hel.Doca_Zh(),hel.Doca_Tof(),hel.T0(),
00188 hel.Doca_Wamb(),hel.Doca_Eang())-doca;
00189
00190 }
00191
00192 std::vector<float> MdcxHit::derivatives(MdcxHel &hel) const {
00193
00194 std::vector<float> deriv = hel.derivatives(*this);
00195 float dtemp = d(hel.Doca_Zh(), hel.Doca_Tof(), hel.T0(),
00196 hel.Doca_Wamb(), hel.Doca_Eang());
00197
00198 deriv[0] = dtemp - deriv[0];
00199
00200 float ewire = e(dtemp);
00201 for (uint32_t i = 0; i < deriv.size(); i++) deriv[i] /= ewire;
00202 return deriv;
00203 }
00204
00205 void MdcxHit::print(ostream &o, int i) const {
00206 o << "(" << Layer() << "," << WireNo() << ",id "<< getDigi()->getTrackIndex()<<") " ;
00207 }
00208
00209 void MdcxHit::printAll(ostream &o, int i) const {
00210 o << " hit" << i << " (" << Layer() << "," << WireNo() << ") ";
00211 if (getMdcHit()) o << " phi "<< getMdcHit()->phi();
00212 o << "dd " << d() << " dde "
00213 << e() << " rt " << t() << endl;
00214 }
00215
00216 double MdcxHit::getSigma(float driftDist, int ambig, double entranceAngle,
00217 double dipAngle, double z) const {
00218 int lrCalib = 2;
00219 if (ambig != 0) lrCalib = (ambig == 1) ? 0 : 1;
00220 double eAngle = EntranceAngle(entranceAngle);
00221
00222 double sig = 0.1 * m_mdcCalibFunSvc->getSigma(_layernumber, lrCalib,
00223 driftDist*10., eAngle, tan(dipAngle), z*10., _iAdc);
00224
00225
00226 #ifdef MDCXRECO_RESLAYER
00227 if (Layer() == g_resLayer) {
00228
00229 return 9999.;
00230 }
00231 #endif
00232 return sig;
00233 }