00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include <iostream>
00014 #include <math.h>
00015 #include <float.h>
00016 #include "TrackUtil/Helix.h"
00017 #include "CLHEP/Matrix/Matrix.h"
00018 #include "GaudiKernel/StatusCode.h"
00019 #include "GaudiKernel/IInterface.h"
00020 #include "GaudiKernel/Kernel.h"
00021 #include "GaudiKernel/Service.h"
00022 #include "GaudiKernel/ISvcLocator.h"
00023 #include "GaudiKernel/SvcFactory.h"
00024 #include "GaudiKernel/IDataProviderSvc.h"
00025 #include "GaudiKernel/Bootstrap.h"
00026 #include "GaudiKernel/MsgStream.h"
00027 #include "GaudiKernel/SmartDataPtr.h"
00028 #include "GaudiKernel/IMessageSvc.h"
00029
00030
00031
00032
00033
00034
00035
00036 const double
00037 M_PI2 = 2. * M_PI;
00038
00039 const double
00040 M_PI4 = 4. * M_PI;
00041
00042 const double
00043 M_PI8 = 8. * M_PI;
00044
00045 const double Helix::ConstantAlpha = 333.564095;
00046
00047 Helix::Helix(const HepPoint3D & pivot,
00048 const HepVector & a,
00049 const HepSymMatrix & Ea)
00050 :
00051
00052 m_pivot(pivot),
00053 m_a(a),
00054 m_matrixValid(true),
00055 m_Ea(Ea) {
00056 StatusCode scmgn = Gaudi::svcLocator()->service("MagneticFieldSvc",m_pmgnIMF);
00057 if(scmgn!=StatusCode::SUCCESS) {
00058 std::cout<< "Unable to open Magnetic field service"<<std::endl;
00059 }
00060 m_bField = -10000*(m_pmgnIMF->getReferField());
00061 m_alpha = 10000. / 2.99792458 / m_bField;
00062
00063
00064 updateCache();
00065 }
00066
00067 Helix::Helix(const HepPoint3D & pivot,
00068 const HepVector & a)
00069 :
00070
00071 m_pivot(pivot),
00072 m_a(a),
00073 m_matrixValid(false),
00074 m_Ea(HepSymMatrix(5,0)) {
00075 StatusCode scmgn = Gaudi::svcLocator()->service("MagneticFieldSvc",m_pmgnIMF);
00076 if(scmgn!=StatusCode::SUCCESS) {
00077
00078 std::cout<< "Unable to open Magnetic field service"<<std::endl;
00079 }
00080 m_bField = -10000*(m_pmgnIMF->getReferField());
00081 m_alpha = 10000. / 2.99792458 / m_bField;
00082
00083
00084 updateCache();
00085 }
00086
00087 Helix::Helix(const HepPoint3D & position,
00088 const Hep3Vector & momentum,
00089 double charge)
00090 :
00091
00092 m_pivot(position),
00093 m_a(HepVector(5,0)),
00094 m_matrixValid(false),
00095 m_Ea(HepSymMatrix(5,0)) {
00096 StatusCode scmgn = Gaudi::svcLocator()->service("MagneticFieldSvc",m_pmgnIMF);
00097 if(scmgn!=StatusCode::SUCCESS) {
00098
00099 std::cout<< "Unable to open Magnetic field service"<<std::endl;
00100 }
00101 m_bField = -10000*(m_pmgnIMF->getReferField());
00102 m_alpha = 10000. / 2.99792458 / m_bField;
00103
00104 m_a[0] = 0.;
00105 m_a[1] = fmod(atan2(- momentum.x(), momentum.y())
00106 + M_PI4, M_PI2);
00107 m_a[3] = 0.;
00108 double perp(momentum.perp());
00109 if (perp != 0.0) {
00110 m_a[2] = charge / perp;
00111 m_a[4] = momentum.z() / perp;
00112 }
00113 else {
00114 m_a[2] = charge * (DBL_MAX);
00115 if (momentum.z() >= 0) {
00116 m_a[4] = (DBL_MAX);
00117 } else {
00118 m_a[4] = -(DBL_MAX);
00119 }
00120 }
00121
00122 updateCache();
00123 }
00124
00125 Helix::~Helix() {
00126 }
00127
00128 HepPoint3D
00129 Helix::x(double phi) const {
00130
00131
00132
00133
00134
00135
00136
00137
00138 double x = m_pivot.x() + m_ac[0] * m_cp + m_r * (m_cp - cos(m_ac[1] +phi));
00139 double y = m_pivot.y() + m_ac[0] * m_sp + m_r * (m_sp - sin(m_ac[1] +phi));
00140 double z = m_pivot.z() + m_ac[3] - m_r * m_ac[4] * phi;
00141
00142 return HepPoint3D(x, y, z);
00143 }
00144
00145 double *
00146 Helix::x(double phi, double p[3]) const {
00147
00148
00149
00150
00151
00152
00153
00154
00155 p[0] = m_pivot.x() + m_ac[0] * m_cp + m_r * (m_cp - cos(m_ac[1] + phi));
00156 p[1] = m_pivot.y() + m_ac[0] * m_sp + m_r * (m_sp - sin(m_ac[1] + phi));
00157 p[2] = m_pivot.z() + m_ac[3] - m_r * m_ac[4] * phi;
00158
00159 return p;
00160 }
00161
00162 HepPoint3D
00163 Helix::x(double phi, HepSymMatrix & Ex) const {
00164 double x = m_pivot.x() + m_ac[0] * m_cp + m_r * (m_cp - cos(m_ac[1] +phi));
00165 double y = m_pivot.y() + m_ac[0] * m_sp + m_r * (m_sp - sin(m_ac[1] +phi));
00166 double z = m_pivot.z() + m_ac[3] - m_r * m_ac[4] * phi;
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177 if (m_matrixValid) Ex = m_Ea.similarity(delXDelA(phi));
00178 else Ex = m_Ea;
00179
00180 return HepPoint3D(x, y, z);
00181 }
00182
00183 Hep3Vector
00184 Helix::momentum(double phi) const {
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195 double pt = fabs(m_pt);
00196 double px = - pt * sin(m_ac[1] + phi);
00197 double py = pt * cos(m_ac[1] + phi);
00198 double pz = pt * m_ac[4];
00199
00200 return Hep3Vector(px, py, pz);
00201 }
00202
00203 Hep3Vector
00204 Helix::momentum(double phi, HepSymMatrix & Em) const {
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215 double pt = fabs(m_pt);
00216 double px = - pt * sin(m_ac[1] + phi);
00217 double py = pt * cos(m_ac[1] + phi);
00218 double pz = pt * m_ac[4];
00219
00220 if (m_matrixValid) Em = m_Ea.similarity(delMDelA(phi));
00221 else Em = m_Ea;
00222
00223 return Hep3Vector(px, py, pz);
00224 }
00225
00226 HepLorentzVector
00227 Helix::momentum(double phi, double mass) const {
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239 double pt = fabs(m_pt);
00240 double px = - pt * sin(m_ac[1] + phi);
00241 double py = pt * cos(m_ac[1] + phi);
00242 double pz = pt * m_ac[4];
00243 double E = sqrt(pt*pt*(1.+m_ac[4]*m_ac[4])+mass*mass);
00244
00245 return HepLorentzVector(px, py, pz, E);
00246 }
00247
00248
00249 HepLorentzVector
00250 Helix::momentum(double phi, double mass, HepSymMatrix & Em) const {
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262 double pt = fabs(m_pt);
00263 double px = - pt * sin(m_ac[1] + phi);
00264 double py = pt * cos(m_ac[1] + phi);
00265 double pz = pt * m_ac[4];
00266 double E = sqrt(pt*pt*(1.+m_ac[4]*m_ac[4])+mass*mass);
00267
00268 if (m_matrixValid) Em = m_Ea.similarity(del4MDelA(phi,mass));
00269 else Em = m_Ea;
00270
00271 return HepLorentzVector(px, py, pz, E);
00272 }
00273
00274 HepLorentzVector
00275 Helix::momentum(double phi,
00276 double mass,
00277 HepPoint3D & x,
00278 HepSymMatrix & Emx) const {
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290 double pt = fabs(m_pt);
00291 double px = - pt * sin(m_ac[1] + phi);
00292 double py = pt * cos(m_ac[1] + phi);
00293 double pz = pt * m_ac[4];
00294 double E = sqrt(pt * pt * (1. + m_ac[4] * m_ac[4]) + mass * mass);
00295
00296 x.setX(m_pivot.x() + m_ac[0] * m_cp + m_r * (m_cp - cos(m_ac[1] + phi)));
00297 x.setY(m_pivot.y() + m_ac[0] * m_sp + m_r * (m_sp - sin(m_ac[1] + phi)));
00298 x.setZ(m_pivot.z() + m_ac[3] - m_r * m_ac[4] * phi);
00299
00300 if (m_matrixValid) Emx = m_Ea.similarity(del4MXDelA(phi,mass));
00301 else Emx = m_Ea;
00302
00303 return HepLorentzVector(px, py, pz, E);
00304 }
00305
00306
00307 const HepPoint3D &
00308 Helix::pivot(const HepPoint3D & newPivot) {
00309 const double & dr = m_ac[0];
00310 const double & phi0 = m_ac[1];
00311 const double & kappa = m_ac[2];
00312 const double & dz = m_ac[3];
00313 const double & tanl = m_ac[4];
00314
00315 double rdr = dr + m_r;
00316 double phi = fmod(phi0 + M_PI4, M_PI2);
00317 double csf0 = cos(phi);
00318 double snf0 = (1. - csf0) * (1. + csf0);
00319 snf0 = sqrt((snf0 > 0.) ? snf0 : 0.);
00320 if(phi > M_PI) snf0 = - snf0;
00321
00322 double xc = m_pivot.x() + rdr * csf0;
00323 double yc = m_pivot.y() + rdr * snf0;
00324 double csf, snf;
00325 if(m_r != 0.0) {
00326 csf = (xc - newPivot.x()) / m_r;
00327 snf = (yc - newPivot.y()) / m_r;
00328 double anrm = sqrt(csf * csf + snf * snf);
00329 if(anrm != 0.0) {
00330 csf /= anrm;
00331 snf /= anrm;
00332 phi = atan2(snf, csf);
00333 } else {
00334 csf = 1.0;
00335 snf = 0.0;
00336 phi = 0.0;
00337 }
00338 } else {
00339 csf = 1.0;
00340 snf = 0.0;
00341 phi = 0.0;
00342 }
00343 double phid = fmod(phi - phi0 + M_PI8, M_PI2);
00344 if(phid > M_PI) phid = phid - M_PI2;
00345 double drp = (m_pivot.x() + dr * csf0 + m_r * (csf0 - csf) - newPivot.x())
00346 * csf
00347 + (m_pivot.y() + dr * snf0 + m_r * (snf0 - snf) - newPivot.y()) * snf;
00348 double dzp = m_pivot.z() + dz - m_r * tanl * phid - newPivot.z();
00349
00350 HepVector ap(5);
00351 ap[0] = drp;
00352 ap[1] = fmod(phi + M_PI4, M_PI2);
00353 ap[2] = kappa;
00354 ap[3] = dzp;
00355 ap[4] = tanl;
00356
00357
00358 if (m_matrixValid) m_Ea = m_Ea.similarity(delApDelA(ap));
00359
00360 m_a = ap;
00361 m_pivot = newPivot;
00362
00363
00364 updateCache();
00365 return m_pivot;
00366 }
00367
00368 void
00369 Helix::set(const HepPoint3D & pivot,
00370 const HepVector & a,
00371 const HepSymMatrix & Ea) {
00372 m_pivot = pivot;
00373 m_a = a;
00374 m_Ea = Ea;
00375 m_matrixValid = true;
00376 updateCache();
00377 }
00378
00379 Helix &
00380 Helix::operator = (const Helix & i) {
00381 if (this == & i) return * this;
00382
00383 m_bField = i.m_bField;
00384 m_alpha = i.m_alpha;
00385 m_pivot = i.m_pivot;
00386 m_a = i.m_a;
00387 m_Ea = i.m_Ea;
00388 m_matrixValid = i.m_matrixValid;
00389
00390 m_center = i.m_center;
00391 m_cp = i.m_cp;
00392 m_sp = i.m_sp;
00393 m_pt = i.m_pt;
00394 m_r = i.m_r;
00395 m_ac[0] = i.m_ac[0];
00396 m_ac[1] = i.m_ac[1];
00397 m_ac[2] = i.m_ac[2];
00398 m_ac[3] = i.m_ac[3];
00399 m_ac[4] = i.m_ac[4];
00400
00401 return * this;
00402 }
00403
00404 void
00405 Helix::updateCache(void) {
00406
00407
00408
00409
00410
00411
00412
00413 m_ac[0] = m_a[0];
00414 m_ac[1] = m_a[1];
00415 m_ac[2] = m_a[2];
00416 m_ac[3] = m_a[3];
00417 m_ac[4] = m_a[4];
00418
00419 m_cp = cos(m_ac[1]);
00420 m_sp = sin(m_ac[1]);
00421 if (m_ac[2] != 0.0) {
00422 m_pt = 1. / m_ac[2];
00423 m_r = m_alpha / m_ac[2];
00424 }
00425 else {
00426 m_pt = (DBL_MAX);
00427 m_r = (DBL_MAX);
00428 }
00429
00430 double x = m_pivot.x() + (m_ac[0] + m_r) * m_cp;
00431 double y = m_pivot.y() + (m_ac[0] + m_r) * m_sp;
00432 m_center.setX(x);
00433 m_center.setY(y);
00434 m_center.setZ(0.);
00435 }
00436
00437 HepMatrix
00438 Helix::delApDelA(const HepVector & ap) const {
00439
00440
00441
00442
00443
00444 HepMatrix dApDA(5,5,0);
00445
00446 const double & dr = m_ac[0];
00447 const double & phi0 = m_ac[1];
00448 const double & cpa = m_ac[2];
00449 const double & dz = m_ac[3];
00450 const double & tnl = m_ac[4];
00451
00452 double drp = ap[0];
00453 double phi0p = ap[1];
00454 double cpap = ap[2];
00455 double dzp = ap[3];
00456 double tnlp = ap[4];
00457
00458 double rdr = m_r + dr;
00459 double rdrpr;
00460 if ((m_r + drp) != 0.0) {
00461 rdrpr = 1. / (m_r + drp);
00462 } else {
00463 rdrpr = (DBL_MAX);
00464 }
00465
00466
00467 double csfd = cos(phi0p - phi0);
00468 double snfd = sin(phi0p - phi0);
00469 double phid = fmod(phi0p - phi0 + M_PI8, M_PI2);
00470 if (phid > M_PI) phid = phid - M_PI2;
00471
00472 dApDA[0][0] = csfd;
00473 dApDA[0][1] = rdr*snfd;
00474 if(cpa!=0.0) {
00475 dApDA[0][2] = (m_r/cpa)*( 1.0 - csfd );
00476 } else {
00477 dApDA[0][2] = (DBL_MAX);
00478 }
00479
00480 dApDA[1][0] = - rdrpr*snfd;
00481 dApDA[1][1] = rdr*rdrpr*csfd;
00482 if(cpa!=0.0) {
00483 dApDA[1][2] = (m_r/cpa)*rdrpr*snfd;
00484 } else {
00485 dApDA[1][2] = (DBL_MAX);
00486 }
00487
00488 dApDA[2][2] = 1.0;
00489
00490 dApDA[3][0] = m_r*rdrpr*tnl*snfd;
00491 dApDA[3][1] = m_r*tnl*(1.0 - rdr*rdrpr*csfd);
00492 if(cpa!=0.0) {
00493 dApDA[3][2] = (m_r/cpa)*tnl*(phid - m_r*rdrpr*snfd);
00494 } else {
00495 dApDA[3][2] = (DBL_MAX);
00496 }
00497 dApDA[3][3] = 1.0;
00498 dApDA[3][4] = - m_r*phid;
00499
00500 dApDA[4][4] = 1.0;
00501
00502 return dApDA;
00503 }
00504
00505 HepMatrix
00506 Helix::delXDelA(double phi) const {
00507
00508
00509
00510
00511
00512
00513 HepMatrix dXDA(3,5,0);
00514
00515 const double & dr = m_ac[0];
00516 const double & phi0 = m_ac[1];
00517 const double & cpa = m_ac[2];
00518 const double & dz = m_ac[3];
00519 const double & tnl = m_ac[4];
00520
00521 double cosf0phi = cos(phi0 + phi);
00522 double sinf0phi = sin(phi0 + phi);
00523
00524 dXDA[0][0] = m_cp;
00525 dXDA[0][1] = - dr * m_sp + m_r * (- m_sp + sinf0phi);
00526 if(cpa!=0.0) {
00527 dXDA[0][2] = - (m_r / cpa) * (m_cp - cosf0phi);
00528 } else {
00529 dXDA[0][2] = (DBL_MAX);
00530 }
00531
00532
00533
00534 dXDA[1][0] = m_sp;
00535 dXDA[1][1] = dr * m_cp + m_r * (m_cp - cosf0phi);
00536 if(cpa!=0.0) {
00537 dXDA[1][2] = - (m_r / cpa) * (m_sp - sinf0phi);
00538 } else {
00539 dXDA[1][2] = (DBL_MAX);
00540 }
00541
00542
00543
00544
00545
00546 if(cpa!=0.0) {
00547 dXDA[2][2] = (m_r / cpa) * tnl * phi;
00548 } else {
00549 dXDA[2][2] = (DBL_MAX);
00550 }
00551 dXDA[2][3] = 1.0;
00552 dXDA[2][4] = - m_r * phi;
00553
00554 return dXDA;
00555 }
00556
00557
00558
00559 HepMatrix
00560 Helix::delMDelA(double phi) const {
00561
00562
00563
00564
00565
00566
00567 HepMatrix dMDA(3,5,0);
00568
00569 const double & phi0 = m_ac[1];
00570 const double & cpa = m_ac[2];
00571 const double & tnl = m_ac[4];
00572
00573 double cosf0phi = cos(phi0+phi);
00574 double sinf0phi = sin(phi0+phi);
00575
00576 double rho;
00577 if(cpa != 0.)rho = 1./cpa;
00578 else rho = (DBL_MAX);
00579
00580 double charge = 1.;
00581 if(cpa < 0.)charge = -1.;
00582
00583 dMDA[0][1] = -fabs(rho)*cosf0phi;
00584 dMDA[0][2] = charge*rho*rho*sinf0phi;
00585
00586 dMDA[1][1] = -fabs(rho)*sinf0phi;
00587 dMDA[1][2] = -charge*rho*rho*cosf0phi;
00588
00589 dMDA[2][2] = -charge*rho*rho*tnl;
00590 dMDA[2][4] = fabs(rho);
00591
00592 return dMDA;
00593 }
00594
00595
00596 HepMatrix
00597 Helix::del4MDelA(double phi, double mass) const {
00598
00599
00600
00601
00602
00603
00604 HepMatrix d4MDA(4,5,0);
00605
00606 double phi0 = m_ac[1];
00607 double cpa = m_ac[2];
00608 double tnl = m_ac[4];
00609
00610 double cosf0phi = cos(phi0+phi);
00611 double sinf0phi = sin(phi0+phi);
00612
00613 double rho;
00614 if(cpa != 0.)rho = 1./cpa;
00615 else rho = (DBL_MAX);
00616
00617 double charge = 1.;
00618 if(cpa < 0.)charge = -1.;
00619
00620 double E = sqrt(rho*rho*(1.+tnl*tnl)+mass*mass);
00621
00622 d4MDA[0][1] = -fabs(rho)*cosf0phi;
00623 d4MDA[0][2] = charge*rho*rho*sinf0phi;
00624
00625 d4MDA[1][1] = -fabs(rho)*sinf0phi;
00626 d4MDA[1][2] = -charge*rho*rho*cosf0phi;
00627
00628 d4MDA[2][2] = -charge*rho*rho*tnl;
00629 d4MDA[2][4] = fabs(rho);
00630
00631 if (cpa != 0.0 && E != 0.0) {
00632 d4MDA[3][2] = (-1.-tnl*tnl)/(cpa*cpa*cpa*E);
00633 d4MDA[3][4] = tnl/(cpa*cpa*E);
00634 } else {
00635 d4MDA[3][2] = (DBL_MAX);
00636 d4MDA[3][4] = (DBL_MAX);
00637 }
00638 return d4MDA;
00639 }
00640
00641
00642 HepMatrix
00643 Helix::del4MXDelA(double phi, double mass) const {
00644
00645
00646
00647
00648
00649
00650 HepMatrix d4MXDA(7,5,0);
00651
00652 const double & dr = m_ac[0];
00653 const double & phi0 = m_ac[1];
00654 const double & cpa = m_ac[2];
00655 const double & dz = m_ac[3];
00656 const double & tnl = m_ac[4];
00657
00658 double cosf0phi = cos(phi0+phi);
00659 double sinf0phi = sin(phi0+phi);
00660
00661 double rho;
00662 if(cpa != 0.)rho = 1./cpa;
00663 else rho = (DBL_MAX);
00664
00665 double charge = 1.;
00666 if(cpa < 0.)charge = -1.;
00667
00668 double E = sqrt(rho * rho * (1. + tnl * tnl) + mass * mass);
00669
00670 d4MXDA[0][1] = - fabs(rho) * cosf0phi;
00671 d4MXDA[0][2] = charge * rho * rho * sinf0phi;
00672
00673 d4MXDA[1][1] = - fabs(rho) * sinf0phi;
00674 d4MXDA[1][2] = - charge * rho * rho * cosf0phi;
00675
00676 d4MXDA[2][2] = - charge * rho * rho * tnl;
00677 d4MXDA[2][4] = fabs(rho);
00678
00679 if (cpa != 0.0 && E != 0.0) {
00680 d4MXDA[3][2] = (- 1. - tnl * tnl) / (cpa * cpa * cpa * E);
00681 d4MXDA[3][4] = tnl / (cpa * cpa * E);
00682 } else {
00683 d4MXDA[3][2] = (DBL_MAX);
00684 d4MXDA[3][4] = (DBL_MAX);
00685 }
00686
00687 d4MXDA[4][0] = m_cp;
00688 d4MXDA[4][1] = - dr * m_sp + m_r * (- m_sp + sinf0phi);
00689 if (cpa != 0.0) {
00690 d4MXDA[4][2] = - (m_r / cpa) * (m_cp - cosf0phi);
00691 } else {
00692 d4MXDA[4][2] = (DBL_MAX);
00693 }
00694
00695 d4MXDA[5][0] = m_sp;
00696 d4MXDA[5][1] = dr * m_cp + m_r * (m_cp - cosf0phi);
00697 if (cpa != 0.0) {
00698 d4MXDA[5][2] = - (m_r / cpa) * (m_sp - sinf0phi);
00699
00700 d4MXDA[6][2] = (m_r / cpa) * tnl * phi;
00701 } else {
00702 d4MXDA[5][2] = (DBL_MAX);
00703
00704 d4MXDA[6][2] = (DBL_MAX);
00705 }
00706
00707 d4MXDA[6][3] = 1.;
00708 d4MXDA[6][4] = - m_r * phi;
00709
00710 return d4MXDA;
00711 }
00712
00713 void
00714 Helix::ignoreErrorMatrix(void) {
00715 m_matrixValid = false;
00716 m_Ea *= 0.;
00717 }