#include <MdcSegInfoSterO.h>
Inheritance diagram for MdcSegInfoSterO:
Public Member Functions | |
| double | arc () const |
| double | arc () const |
| int | calcStereo (MdcSeg *parentSeg, const TrkExchangePar &par, MdcSegWorks &segStuff, double Bz) |
| int | calcStereo (MdcSeg *parentSeg, const TrkRecoTrk &track, MdcSegWorks &segStuff) |
| int | calcStereo (MdcSeg *parentSeg, const TrkExchangePar &par, MdcSegWorks &segStuff, double Bz) |
| int | calcStereo (MdcSeg *parentSeg, const TrkRecoTrk &track, MdcSegWorks &segStuff) |
| double | ct () const |
| double | ct () const |
| const double * | errmat () const |
| const double * | errmat () const |
| const double * | inverr () const |
| const double * | inverr () const |
| MdcSegInfoSterO () | |
| MdcSegInfoSterO () | |
| double | par (int i) const |
| double | par (int i) const |
| bool | parIsAngle (int i) const |
| bool | parIsAngle (int i) const |
| void | plotSegInfo () const |
| void | plotSegInfo () const |
| double | z0 () const |
| double | z0 () const |
| ~MdcSegInfoSterO () | |
| ~MdcSegInfoSterO () | |
Protected Attributes | |
| double | _arc |
| double | _errmat [3] |
| double | _inverr [3] |
| double | _par0 |
| double | _par1 |
Private Member Functions | |
| MdcSegInfoSterO (const MdcSegInfoSterO &) | |
| MdcSegInfoSterO (const MdcSegInfoSterO &) | |
| MdcSegInfoSterO & | operator= (const MdcSegInfoSterO &) |
| MdcSegInfoSterO & | operator= (const MdcSegInfoSterO &) |
| int | zPosition (MdcHitUse &hitUse, const TrkExchangePar &par, MdcLine *span, int hitFit, double t0, double Bz) const |
| int | zPosition (MdcHitUse &hitUse, const TrkRecoTrk &track, MdcLine *span, int hitFit, double t0) const |
| int | zPosition (MdcHitUse &hitUse, const TrkExchangePar &par, MdcLine *span, int hitFit, double t0, double Bz) const |
| int | zPosition (MdcHitUse &hitUse, const TrkRecoTrk &track, MdcLine *span, int hitFit, double t0) const |
Private Attributes | |
| int | _debug |
|
|
00037 { };
|
|
|
00038 { };
|
|
|
|
|
|
00037 { };
|
|
|
00038 { };
|
|
|
|
|
|
00044 {return _arc;}
|
|
|
00044 {return _arc;}
|
|
||||||||||||||||||||
|
|
|
||||||||||||||||
|
|
|
||||||||||||||||||||
|
00065 {
00066 //----------------------------------------------------------------------------
00067 _debug = 0;//yzhang debug
00068 // phiAx = approximate phi of axial track @ stereo slayer
00069 // return 1 if calculation failed
00070 const MdcSuperLayer *slayer = parentSeg->superlayer();
00071
00072 double kap = 0.5 * par.omega();
00073 double d0 = par.d0();
00074 BesAngle phi0(par.phi0());
00075 BesAngle phiSeg(parentSeg->phi());
00076 double radius = slayer->rEnd();
00077 double dPhiZ = (-1)*slayer->delPhiinv() * slayer->zEnd();
00078 if (_debug >0){
00079 parentSeg->plotSeg();
00080 std::cout << "delPhi "<<slayer->delPhi()
00081 <<" delPhiInv "<<slayer->delPhiinv()
00082 <<" dPhiZ "<<dPhiZ
00083 <<" zEnd() "<<slayer->zEnd()
00084 <<" phiAx "<<segStuff.phiAx
00085 <<" phiSeg "<<phiSeg
00086 <<" phiDiff "<<phiSeg - segStuff.phiAx<< std::endl;//yzhang debug
00087 }
00088 bool lStraight = false;
00089 if (fabs(kap)<1.e-9) lStraight = true;
00090
00091 BesAngle phiDiff = phiSeg - segStuff.phiAx;
00092 double zApprox = phiDiff.rad() * -dPhiZ;
00093 if (_debug >0){
00094 std::cout<<"zApp "<<zApprox<<std::endl;//yzhang debug
00095 }
00096 double delRad = slayer->stDip() *
00097 (1. - zApprox * zApprox * segStuff.wirLen2inv);
00098 radius -= delRad;
00099 double delPhiArg;
00100 if (lStraight) {
00101 delPhiArg = delRad * d0 / (radius*radius);
00102 }else {
00103 delPhiArg = -delRad * kap;
00104 }
00105
00106
00107 //---------------------yzhang fix-------------------
00108 MdcLine span(12);
00109 int hitFit=0;
00110 bool szFaild = false;
00111 double arcTemp=0;
00112
00113 for (int ihit = 0; ihit < parentSeg->nHit(); ihit++){
00114 //*** Calc. z and correct for stereo dip
00115 //int szCode = zPosition(*(parentSeg->hit(ihit)), track, &span,hitFit,parentSeg->bunchTime());
00116 int szCode = zPosition(*(parentSeg->hit(ihit)), par, &span,hitFit,parentSeg->bunchTime(),Bz);
00117
00118 if (szCode > 0){
00119 span.sigma[hitFit] = 1;
00120 arcTemp += span.x[hitFit];
00121 hitFit ++;
00122 }else{
00123 szFaild = true;
00124 if (_debug >0){
00125 parentSeg->hit(ihit)->mdcHit()->print(std::cout);
00126 std::cout<<"MdcSegInfoSterO hit "<<ihit<<" z calc faild"<<std::endl;//yzhang debug
00127 }
00128 }
00129 }
00130 if (hitFit >0) span.fit(hitFit);
00131 //zhangy
00132
00133 // series expand for asin, dropping (delphiarg)**2 and higher terms
00134 segStuff.phiAx += delPhiArg +
00135 0.5 * delPhiArg * segStuff.phiArg * segStuff.phiArg;
00136 phiDiff = phiSeg - segStuff.phiAx;
00137 double z = phiDiff.rad() * -dPhiZ;
00138 if (_debug >0){
00139 std::cout<<"z "<<z<<std::endl;//yzhang debug
00140 }
00141 double ct, z0, arc;
00142 if (lStraight) { //straight track
00143 //*** Calculate parameters
00144 double arg = radius*radius - d0*d0;
00145 if (arg <= 0.0) return 1;
00146 double rD0Root = sqrt(arg);
00147 double rD0Rinv = 1. / rD0Root;
00148 double rinv = 1./radius;
00149 // ct
00150 double slope = parentSeg->slope();
00151 ct = dPhiZ * (rD0Root * slope + d0 * rinv) * rinv;
00152 // z0
00153 arc = TrkHelixUtils::fltToRad(par, radius);
00154 if (arc == 0.0) return 1;
00155 z0 = z - ct * arc;
00156
00157 if (_debug>0){
00158 std::cout << "in MdcSegInfoSterO : z0 "<<z0 <<" "<<_errmat[0]
00159 <<" ct "<<ct<<" "<<_errmat[1]
00160 <<" arc "<<arc<<" "<<_errmat[2]
00161 <<endl;//yzhang debug
00162 }
00163 // calculate errors
00164 double dctdm = dPhiZ * rD0Root * rinv;
00165 double dctdD = -dPhiZ * rinv * (slope * d0 * rD0Rinv - rinv);
00166 double dzdm = -arc * dPhiZ * rD0Root * rinv;
00167 double dzdphi = dPhiZ;
00168 double dzdphi0 = -dPhiZ;
00169 double dzdD = -dPhiZ + ct * d0 * rD0Rinv - arc * dctdD;
00170
00171 const double *inErr = parentSeg->errmat();
00172 //z0
00173 const HepMatrix trackErr = par.covariance();
00174 _errmat[0] = dzdm * dzdm * inErr[2] +
00175 2 * dzdm * dzdphi * inErr[1] +
00176 dzdphi * dzdphi * inErr[0] +
00177 dzdD * dzdD * trackErr(1,1) +
00178 dzdphi0 * dzdphi0 * trackErr(2,2) +
00179 2 * dzdD * dzdphi0 * trackErr(1,2);
00180 if (_errmat[0] < 0.) _errmat[0] = 0.;
00181
00182 //ct
00183 _errmat[2] = dctdm * dctdm * inErr[2] +
00184 dctdD * dctdD * trackErr(1,1);
00185 if (_errmat[2] < 0.) _errmat[2] = 0.;
00186
00187 // off-diag
00188 _errmat[1] = dzdm * dctdm * inErr[2] +
00189 dzdphi * dctdm * inErr[1] +
00190 dzdD * dctdD * trackErr(1,1) +
00191 dzdphi0 * dctdD * trackErr(1,2);
00192 }
00193 else { // curved track (treated as from origin)
00194 //*** Calculate parameters
00195 // ct
00196 double arg = 1. - kap*kap * radius*radius;
00197 if (arg < 0.0) return 1;
00198 double rKapRoot = sqrt(arg);
00199 double rKapRinv = 1. / rKapRoot;
00200 double ctFactor = -rKapRoot * -dPhiZ;
00201 double slopeAx = kap * rKapRinv;
00202 double slope = parentSeg->slope();
00203
00204 //std::cout<<"slopeAx="<<slopeAx<<" slopeSeg="<<slope
00205 // <<"\n diff "<<slopeAx - slope
00206 // <<" fac "<<ctFactor<<" rKapRoot "<<rKapRoot
00207 // <<" dPhiZ "<<dPhiZ<<std::endl;//yzhang debug
00208 //inner use szPozision
00209 //if (slayer->index()>1 || szFaild){//inner chamber use sz algorithm
00210 if (szFaild){ //use orinal algorithm
00211 //original
00212 ct = (slopeAx - slope) * ctFactor;
00213 arc = TrkHelixUtils::fltToRad(par, radius);
00214 if (arc == 0.0) return 1;
00215 z0 = z - ct * arc;
00216 }else{
00217 //yzhang use belle s/z
00218 if (!szFaild){
00219 arc = arcTemp/hitFit;
00220 ct = span.slope;
00221 z0 = span.intercept;
00222 }else{
00223 ct = 999;
00224 z0 = 999;
00225 arc = 999;
00226 }
00227 }
00228 if (_debug >0){//yzhang debug
00229 std::cout<<"--slayer NO. "<<slayer->index()<<std::endl;//yzhang debug
00230 std::cout<<"ori ct "<<(slopeAx - slope) * ctFactor
00231 <<" z0 "<<z - ct * TrkHelixUtils::fltToRad(par, radius)
00232 <<" arc "<<TrkHelixUtils::fltToRad(par, radius)<<std::endl;
00233 std::cout<<"fix ct "<<span.slope<<" z0 "<<span.intercept
00234 <<" arc "<<arcTemp/hitFit<<std::endl;
00235 std::cout<<"-------- "<<std::endl;//yzhang debug
00236 }
00237 // Calculate errors -- eliminate a bunch of terms when I get around to it
00238 double dctdm = dPhiZ * rKapRoot;
00239 double dctdkap = -dPhiZ * ( 1 + radius * radius * kap * rKapRinv * slope);
00240 double dzdm = arc * -dPhiZ * rKapRoot;
00241 double dzdphi = dPhiZ;
00242 double dzdkap = dPhiZ * radius * rKapRinv - arc * dctdkap -
00243 ct * ( radius * rKapRinv / kap - arc / kap);
00244 double dzdphi0 = -dPhiZ ;
00245
00246 const double *inErr = parentSeg->errmat();
00247 //z0
00248 const HepMatrix trackErr = par.covariance();
00249 _errmat[0] = dzdm * dzdm * inErr[2] +
00250 2 * dzdm * dzdphi * inErr[1] +
00251 dzdphi * dzdphi * inErr[0] +
00252 dzdkap * dzdkap * trackErr(3,3) +
00253 dzdphi0 * dzdphi0 * trackErr(2,2) +
00254 2 * dzdkap * dzdphi0 * trackErr(2,3);
00255 /*
00256 std::cout<<"dzdm "<<dzdm
00257 <<" inErr[2] "<<inErr[2]
00258 <<" inErr[1] "<<inErr[1]
00259 <<" inErr[0] "<<inErr[0]
00260 <<" dzdphi "<<dzdphi
00261 <<" dzdkap "<<dzdkap
00262 <<" dzdphi0 "<<dzdphi0
00263 <<" trackErr3,3 "<<trackErr(3,3)
00264 <<" trackErr2,2 "<<trackErr(2,2)
00265 <<" trackErr2,3 "<<trackErr(2,3)
00266 <<std::endl;
00267 */
00268 if (_errmat[0] < 0.) _errmat[0] = 0.;
00269
00270 //ct
00271 _errmat[2] = dctdm * dctdm * inErr[2] +
00272 dctdkap * dctdkap * trackErr(3,3);
00273 if (_errmat[2] < 0.) _errmat[2] = 0.;
00274
00275 //off-diag
00276 _errmat[1] = dzdm * dctdm * inErr[2] +
00277 dzdphi * dctdm * inErr[1] +
00278 dzdkap * dctdkap * trackErr(3,3) +
00279 dzdphi0 * dctdkap * trackErr(2,3);
00280
00281 } // end branch on straight/curved
00282
00283
00284 // Load parameters
00285 _par0 = z0;
00286 _par1 = ct;
00287 _arc = arc;
00288
00289 long error = mdcTwoInv(_errmat, _inverr);
00290 if (error) {
00291 std::cout << " ErrMsg(warning)" << "Failed to invert matrix -- MdcSegInfo::calcStereo" <<
00292 endl << _errmat[0] << " " << _errmat[1] << " " << _errmat[2]<< std::endl;
00293 }
00294 return 0;
00295
00296 }
|
|
||||||||||||||||
|
00055 {
00056 //----------------------------------------------------------------------------
00057 const TrkFit* tkFit = track.fitResult();
00058 if (tkFit == 0) return 1;
00059 TrkExchangePar par = tkFit->helix(0.0);
00060 double Bz = track.bField().bFieldZ();//??
00061 return calcStereo(parentSeg, par, segStuff,Bz);
00062 }
|
|
|
00041 {return _par1;}
|
|
|
00041 {return _par1;}
|
|
|
00039 {return _errmat;}
|
|
|
00039 {return _errmat;}
|
|
|
00040 {return _inverr;}
|
|
|
00040 {return _inverr;}
|
|
|
|
|
|
|
|
|
00043 {return (0 == i) ? _par0 : _par1;}
|
|
|
00043 {return (0 == i) ? _par0 : _par1;}
|
|
|
Implements MdcSegInfo. |
|
|
Implements MdcSegInfo. 00048 {
00049 //---------------------------------------------------------------------------
00050 assert (i >= 0 && i < 2);
00051 return false;
00052 }
|
|
|
|
|
|
00028 {
00029 //------------------------------------------------------------------------
00030 std::cout<<"seginfo: "<<_par0<<" "<<_par1<<" "<<_arc<< std::endl;
00031 }
|
|
|
00040 {return _par0;}
|
|
|
00040 {return _par0;}
|
|
||||||||||||||||||||||||||||
|
|
|
||||||||||||||||||||||||
|
|
|
||||||||||||||||||||||||||||
|
00298 {
00299
00300 const MdcHit & h = *(hitUse.mdcHit());
00301
00302 HepPoint3D fp ( h.wire()->getWestPoint()->x(),h.wire()->getWestPoint()->y(),h.wire()->getWestPoint()->z());
00303 HepPoint3D rp ( h.wire()->getEastPoint()->x(),h.wire()->getEastPoint()->y(),h.wire()->getEastPoint()->z());
00304
00305 HepVector3D X = 0.5 * (fp + rp);
00306 if(_debug>0){
00307 std::cout<<"---------- "<<std::endl;//yzhang debug
00308 //h.wire()->print(std::cout);
00309 std::cout<<"fp rp "<<fp<<" "<<rp<<std::endl;//yzhang debug
00310 std::cout<<"Xmid "<<X<<std::endl;//yzhang debug
00311 }
00312 HepVector3D xx = HepVector3D(X.x(), X.y(), 0.);
00313 //center of helix/circle
00314 //change to belle param
00315
00316 double d0 = -par.d0();
00317 double phi0 = (par.phi0()-pi/2);
00318 double _charge = -1;
00319 if((-1)*par.omega()/Bz > 0) _charge = 1;
00320
00321 double r;//??
00322 if (fabs(par.omega())< Constants::epsilon) r = 9999.;
00323 else r = 1/par.omega()*Bz;//???
00324
00325 double xc = cos(phi0) * (d0 + r);
00326 double yc = sin(phi0) * (d0 + r);
00327 HepPoint3D center (xc, yc, 0. );
00328
00329 HepVector3D yy = center - xx;
00330 HepVector3D ww = HepVector3D(yy.x(), yy.y(), 0.);
00331 double wwmag2 = ww.mag2();
00332 double wwmag = sqrt(wwmag2);
00333 //dirftDist
00334 double driftdist = fabs( h.driftDist(t0,hitUse.ambig()) );
00335 HepVector3D lr(driftdist/wwmag * ww.x(),
00336 driftdist/wwmag * ww.y(), 0.);
00337 if (_debug >0){
00338 std::cout<<"r1 d0 "<<r <<" "<<d0<<std::endl;//
00339 std::cout<<"lr "<<lr<<" LR "<<hitUse.ambig()
00340 <<" left "<<h.driftDist(0,1)
00341 <<" right "<<h.driftDist(0,-1)<<std::endl;
00342 }
00343
00344 //...Check left or right...
00345 HepPoint3D ORIGIN = HepPoint3D(0., 0., 0.);
00346 //ambig
00347 //-1 right -phi direction driftDist-=lr
00348 //, +1 left +phi directon driftDist+=lr
00349 if (_charge < 0){
00350 if (hitUse.ambig() == -1){
00351 lr = - lr;//left
00352 }else if (hitUse.ambig() == 0){
00353 lr = ORIGIN;//don't know
00354 }
00355 }else{
00356 if (hitUse.ambig() == 1){
00357 lr = - lr;//left
00358 }else if (hitUse.ambig() == 0){
00359 lr = ORIGIN;//don't know
00360 }
00361 }
00362 X += lr;
00363
00364 //...Prepare vectors...
00365 // HepPoint3D center = _helix->center();
00366 HepPoint3D tmp(-9999., -9999., 0.);
00367 HepVector3D x = HepVector3D(X.x(), X.y(), 0.);
00368 HepVector3D w = x - center;
00369 // //modified the next sentence because the direction are different from belle.
00370 HepVector3D V = h.wire()->zAxis();
00371 HepVector3D v = HepVector3D(V.x(), V.y(), 0.);
00372 double vmag2 = v.mag2();
00373 //double vmag = sqrt(vmag2);
00374
00375 //double r = _helix->curv();
00376 double wv = w.dot(v);
00377 // wv = abs(wv);
00378 double d2 = wv * wv - vmag2 * (w.mag2() - r * r);
00379 if (_debug >0){
00380 std::cout<<"X_fix "<<X<<" center "<<center<<std::endl;
00381 std::cout<<"V "<<V<<std::endl;//yzhang debug
00382 std::cout<<"w "<<w<<" v "<<v<<std::endl;
00383 }
00384 //...No crossing in R/Phi plane... This is too tight...
00385
00386 if (d2 < 0.) {
00387 //hitUse.position(tmp);
00388 if (_debug>0){
00389 std::cout << "in MdcSegInfoSterO !!! stereo: 0. > d2 = " << d2 << " "
00390 << hitUse.ambig() << std::endl;
00391 }
00392 return -1;
00393 }
00394 double d = sqrt(d2);
00395
00396 //...Cal. length to crossing points...
00397 double l[2];
00398 if (_debug >0){
00399 std::cout<<"wv "<<wv<<" d "<<d<<" vmag2 "<<vmag2<<std::endl;//yzhang debug
00400 }
00401 l[0] = (- wv + d) / vmag2;
00402 l[1] = (- wv - d) / vmag2;
00403
00404 //...Cal. z of crossing points...
00405 bool ok[2];
00406 ok[0] = true;
00407 ok[1] = true;
00408 double z[2];
00409 z[0] = X.z() + l[0] * V.z();
00410 z[1] = X.z() + l[1] * V.z();
00411 if (_debug >0){//yzhang debug
00412 std::cout << "X.z "<<X.z()<<" V.z "<<V.z()<<std::endl;
00413 std::cout << "l0, l1 = " << l[0] << ", " << l[1] << std::endl;
00414 std::cout << "z0, z1 = " << z[0] << ", " << z[1] << std::endl;
00415 }
00416 //...Check z position...
00417 if (hitUse.ambig() == 0) {
00418 if (z[0] > rp.z()+20.
00419 || z[0] < fp.z()-20.) ok[0] = false;
00420 if (z[1] > rp.z()+20.
00421 || z[1] < fp.z()-20.) ok[1] = false;
00422 }
00423 else {
00424 if (z[0] > rp.z()
00425 || z[0] < fp.z() ) ok[0] = false;
00426 if (z[1] > rp.z()
00427 || z[1] < fp.z() ) ok[1] = false;
00428 }
00429 if ((! ok[0]) && (! ok[1])) {
00430 //hitUse.position(tmp);
00431 return -2;
00432 }
00433
00434 //...Cal. xy position of crossing points...
00435 HepVector3D p[2];
00436 p[0] = x + l[0] * v;
00437 p[1] = x + l[1] * v;
00438
00439 if(_debug>0){
00440 std::cout<<__FILE__<<" "<<__LINE__<< " p0 "<<p[0].x()<<" "<<p[0].y()<< std::endl;
00441 std::cout<<__FILE__<<" "<<__LINE__<< " p1 "<<p[1].x()<<" "<<p[1].y()<< std::endl;
00442 std::cout<<__FILE__<<" "<<__LINE__<< " c "<<center.x()<<" "<<center.y()<<" "<<_charge<< std::endl;
00443 }
00444
00445 if (_charge * (center.x() * p[0].y() - center.y() * p[0].x()) < 0.)
00446 ok[0] = false;
00447 if (_charge * (center.x() * p[1].y() - center.y() * p[1].x()) < 0.)
00448 ok[1] = false;
00449
00450 if ((! ok[0]) && (! ok[1])){
00451 return -3;
00452 }
00453
00454 //...Which one is the best?... Study needed...
00455 unsigned best = 0;
00456 if (ok[1]) best = 1;
00457 //std::cout<<"in MdcSegInfoSterO Zbelle "<<z[best]<<std::endl;//yzhang debug
00458
00459 //...Cal. arc length...
00460 double cosdPhi = - center.dot((p[best] - center).unit()) / center.mag();
00461 double dPhi;
00462 if(fabs(cosdPhi)<=1.0) {
00463 dPhi = acos(cosdPhi);
00464 } else if (cosdPhi>1.0) {
00465 dPhi = 0.0;
00466 } else {
00467 dPhi = pi;
00468 }
00469
00470 //...Finish...
00471 tmp.setX(abs(r * dPhi));
00472 tmp.setY(z[best]);
00473 //hitUse.position(tmp);
00474
00475 //z
00476 span->y[hitFit] = z[best];
00477 //if (ok[0]) span->y[hitFit] = p[0].mag();//r
00478 //else if(ok[1]) span->y[hitFit] = p[1].mag();
00479 //else span->y[hitFit] = tmp.mag();
00480 span->x[hitFit] = abs(r * dPhi);
00481
00482 if (_debug >0){
00483 h.print(std::cout);
00484 std::cout<<"s z "<<span->x[hitFit]<<" "<<span->y[hitFit]<<std::endl;//yzhang debug
00485 }
00486 return 1;
00487 }
|
|
||||||||||||||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
1.3.9.1