00001
00002
00003 #include "MdcData/MdcHitOnTrack.h"
00004 #include "MdcData/MdcHit.h"
00005 #include "MdcGeom/Constants.h"
00006 #include "MdcGeom/BesAngle.h"
00007
00008
00009 #include "MdcGeom/MdcLayer.h"
00010 #include "MdcGeom/MdcSWire.h"
00011 #include "TrkBase/TrkPoca.h"
00012 #include "TrkBase/TrkDifTraj.h"
00013 #include "TrkBase/TrkRep.h"
00014 #include "TrkBase/TrkSimpTraj.h"
00015 #include "TrkBase/TrkRecoTrk.h"
00016 #include "MdcGeom/EntranceAngle.h"
00017 using std::cout;
00018 using std::endl;
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047 MdcHitOnTrack::MdcHitOnTrack(const TrkFundHit& fundHit,
00048 const MdcHit& baseHit,
00049 int ambig, double t0)
00050 : TrkHitOnTrk(&fundHit,10.e-4),
00051 _ambig(ambig),
00052 _hitTraj(baseHit.hitTraj()),
00053 _fitTime(baseHit.rawTime()-t0*1e-9),
00054 _dHit(&baseHit)
00055 {
00056
00057
00058
00059
00060
00061 _drift[0] = +9999;
00062 _drift[1] = -9999;
00063 setHitResid(-21212121.0);
00064 setHitRms( 0.02 );
00065 setHitLen(0.5 * baseHit.layer()->zLength());
00066 setFltLen(0.);
00067 _startLen = hitTraj()->lowRange() - 5.;
00068 _endLen = hitTraj()->hiRange() + 5.;
00069 }
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099 MdcHitOnTrack::MdcHitOnTrack(const MdcHitOnTrack &hot,TrkRep *newRep,
00100 const TrkDifTraj* trkTraj, const MdcHit *hb)
00101 : TrkHitOnTrk(hot,newRep,trkTraj)
00102 {
00103 _ambig = hot._ambig;
00104 _hitTraj = hot._hitTraj;
00105 _fitTime = hot._fitTime;
00106 _drift[0] = hot._drift[0];
00107 _drift[1] = hot._drift[1];
00108 _startLen = hot._startLen;
00109 _endLen = hot._endLen;
00110 _dHit = (hb==0?hot._dHit:hb);
00111 }
00112
00113 MdcHitOnTrack::~MdcHitOnTrack()
00114 { ; }
00115
00116 void
00117 MdcHitOnTrack::setT0(double t0)
00118 {
00119 _fitTime= _dHit->rawTime()-t0*1e-9;
00120 }
00121
00122 double
00123 MdcHitOnTrack::dcaToWire() const
00124 {
00125 double dca = -9999.;
00126 if ( getParentRep() == 0 ) {
00127
00128 return dca;
00129 }
00130
00131
00132 if (isActive()) {
00133 dca = resid()+drift();
00134 } else {
00135 TrkPoca poca(getParentRep()->traj(), fltLen(), *hitTraj(), hitLen(),
00136 _tolerance);
00137 if (poca.status().success()) dca = poca.doca();
00138 }
00139 return dca;
00140 }
00141
00142 bool
00143 MdcHitOnTrack::updateAmbiguity(double dca)
00144 {
00145 if (dca < 0 && ambig() >= 0) {
00146 setAmbig(-1); return isActive();
00147 } else if (dca > 0 && ambig() <= 0) {
00148 setAmbig(1); return isActive();
00149 } else {
00150 return false;
00151 }
00152 }
00153
00154 const MdcHitOnTrack*
00155 MdcHitOnTrack::mdcHitOnTrack() const
00156 {
00157 return this;
00158 }
00159
00160 double
00161 MdcHitOnTrack::entranceAngleHit() const
00162 {
00163 static Hep3Vector dir;
00164 static HepPoint3D pos;
00165 if (getParentRep() == 0) return 0.;
00166 getParentRep()->traj().getInfo(fltLen(), pos, dir);
00167
00168 return BesAngle(dir.phi() - pos.phi());
00169 }
00170
00171 double
00172 MdcHitOnTrack::entranceAngle() const
00173 {
00174 static Hep3Vector dir;
00175 static HepPoint3D pos;
00176 if (getParentRep() == 0) return 0.;
00177 getParentRep()->traj().getInfo(fltLen(), pos, dir);
00178
00179 return entranceAngle(pos, dir);
00180 }
00181
00182 double
00183 MdcHitOnTrack::entranceAngle(const HepPoint3D pos, const Hep3Vector dir) const
00184 {
00185 double angle = EntranceAngle(dir.phi() - _dHit->phi(pos.z()));
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217 return angle;
00218 }
00219
00220 unsigned
00221 MdcHitOnTrack::layerNumber() const
00222 {
00223 return layernumber();
00224 }
00225
00226 double
00227 MdcHitOnTrack::dipAngle() const
00228 {
00229 return getParentRep()==0?0:Constants::pi/2-getParentRep()->traj().direction(fltLen()).theta();
00230 }
00231
00232 TrkErrCode
00233 MdcHitOnTrack::updateMeasurement(const TrkDifTraj* traj, bool maintainAmb)
00234 {
00235 TrkErrCode status=updatePoca(traj,maintainAmb);
00236 if (status.failure()) {
00237 #ifdef MDCPATREC_DEBUG
00238 std::cout<<" ErrMsg(warning) " << "MdcHitOnTrack::updateMeasurement failed " << status << std::endl;
00239 #endif
00240 return status;
00241 }
00242 assert (_poca!=0);
00243 double dca=_poca->doca();
00244 bool forceIteration = (maintainAmb&&ambig()!=0)?false:updateAmbiguity(dca);
00245
00246
00247
00248 assert(ambig()!=0);
00249
00250 if (isBeyondEndflange()) setUsability(false);
00251 if (forceIteration || !driftCurrent() ) {
00252 updateCorrections();
00253 forceIteration=true;
00254 }
00255 setHitResid(dca-drift());
00256 return !forceIteration?status:
00257 TrkErrCode(TrkErrCode::succeed, 11, "Ambiguity flipped");
00258 }
00259
00260 void
00261 MdcHitOnTrack::updateCorrections()
00262 {
00263 const TrkRep* tkRep = getParentRep();
00264 assert(tkRep != 0);
00265 static HepPoint3D pos; static Hep3Vector dir;
00266 _trkTraj->getInfo(fltLen(), pos, dir);
00267
00268
00269 int wireAmb = wireAmbig();
00270 double tof = tkRep->arrivalTime(fltLen());
00271
00272 assert( ambig() == -1 || ambig() == 1 );
00273 double eAngle = entranceAngle(pos, dir);
00274 double dAngle = Constants::pi/2 - dir.theta();
00275 double z = pos.z();
00276
00277
00278
00279
00280
00281
00282
00283 double dist = _dHit->driftDist(tof, wireAmb, eAngle, dAngle, z);
00284
00285
00286 _fitTime = _dHit->driftTime(tof,z);
00287
00288 _drift[ambig()<0?0:1] = ambig() * dist;
00289 assert( driftCurrent() );
00290
00291 double newSig = _dHit->sigma(dist, wireAmb, eAngle, dAngle, z);
00292
00293
00294 setHitRms(newSig);
00295 }
00296
00297 double
00298 MdcHitOnTrack::driftVelocity() const
00299 {
00300 const TrkRep* tkRep = getParentRep();
00301 assert(tkRep != 0);
00302 double tof = tkRep->arrivalTime(fltLen());
00303 static HepPoint3D pos; static Hep3Vector dir;
00304 _trkTraj->getInfo(fltLen(), pos, dir);
00305 double eAngle = entranceAngle(pos, dir);
00306 double dAngle = Constants::pi/2 - dir.theta();
00307 double z = pos.z();
00308
00309
00310 int wireAmb= fabs(eAngle)<Constants::pi/2?_ambig:-_ambig;
00311
00312 static const double epsilon = 0.5e-9;
00313 double dist1 = _dHit->driftDist(tof+epsilon, wireAmb, eAngle, dAngle, z);
00314 double dist2 = _dHit->driftDist(tof-epsilon, wireAmb, eAngle, dAngle, z);
00315
00316 return (dist2-dist1)/(2*epsilon);
00317 }
00318
00319 bool
00320 MdcHitOnTrack::timeResid(double &t, double &tErr) const
00321 {
00322 double v = driftVelocity();
00323 if (v <= 0) return false;
00324 t = (fabs(drift())-fabs(dcaToWire()))/v;
00325 tErr= hitRms()/v;
00326 return true;
00327 }
00328
00329 bool
00330 MdcHitOnTrack::timeAbsolute(double &t, double &tErr) const
00331 {
00332 double tresid(-1.0);
00333 if(timeResid(tresid,tErr)){
00334
00335 t = tresid + getParentRep()->parentTrack()->trackT0();
00336 return true;
00337 } else
00338 return false;
00339 }
00340
00341 TrkEnums::TrkViewInfo
00342 MdcHitOnTrack::whatView() const
00343 {
00344 return _dHit->whatView();
00345 }
00346
00347 int
00348 MdcHitOnTrack::layernumber() const
00349 {
00350 return _dHit->layernumber();
00351 }
00352
00353 const MdcLayer*
00354 MdcHitOnTrack::layer() const
00355 {
00356 return _dHit->layer();
00357 }
00358
00359 int
00360 MdcHitOnTrack::wire() const
00361 {
00362 return _dHit->wirenumber();
00363 }
00364
00365 int
00366 MdcHitOnTrack::whichView() const
00367 {
00368 return _dHit->whichView();
00369 }
00370
00371 double
00372 MdcHitOnTrack::rawTime() const
00373 {
00374 return _dHit->rawTime();
00375 }
00376
00377 double
00378 MdcHitOnTrack::charge() const
00379 {
00380 return _dHit->charge();
00381 }
00382
00383 const Trajectory*
00384 MdcHitOnTrack::hitTraj() const
00385 {
00386 return _hitTraj;
00387 }
00388
00389 const MdcHit*
00390 MdcHitOnTrack::mdcHit() const
00391 {
00392 return 0;
00393 }
00394
00395
00396
00397
00398 void
00399 MdcHitOnTrack::changeBase(MdcHit* newBase)
00400 {
00401 _dHit = newBase;
00402 }
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415 int MdcHitOnTrack::wireAmbig() const {
00416
00417
00418
00419 const TrkRep* tkRep = getParentRep();
00420 static Hep3Vector dir;
00421 static HepPoint3D pos;
00422 if (getParentRep() == 0) return 0.;
00423 getParentRep()->traj().getInfo(fltLen(), pos, dir);
00424
00425 double wireAmb = ambig();
00426 if (mdcHit()->isCosmicFit()){
00427 HepPoint3D poca = tkRep->position(0.);
00428 if ( pos.y() > poca.y()){
00429 wireAmb = -1.*_ambig;
00430
00431 }
00432 }
00433
00434 return wireAmb;
00435 }