/Home/aooliver/hepsoft/lcdroot/LCDRootApps/PhUtil/src/LCDVToplTRK.cxx

Go to the documentation of this file.
00001 // ----------------------------------------------------------------------------
00002 // $Id: LCDVToplTRK.cxx,v 1.4 2001/10/08 17:54:11 toshi Exp $
00003 // ----------------------------------------------------------------------------
00004 //
00005 #include "LCDVToplTRK.h"
00006 
00007 ////////////////////////////////////////////////////////////
00008 //                                              
00009 // LCDVToplTRK
00010 //                                                  
00011 ////////////////////////////////////////////////////////////
00012 
00013 ClassImp(LCDVToplTRK)
00014 
00015 //__________________________________________________________________________
00016 //Constructor&deconstructor
00017 LCDVToplTRK::LCDVToplTRK() {
00018   Init();
00019 }
00020 
00021 LCDVToplTRK::LCDVToplTRK(Int_t tkid, Double_t* tkpar, Double_t* e_tkpar, 
00022                          Double_t chrg,
00023                          Double_t mfld, Int_t mc_ind,const TVector3& ipv) {
00024   Init(tkid,tkpar,e_tkpar,chrg,mfld,mc_ind,ipv);
00025 }
00026 
00027 LCDVToplTRK::LCDVToplTRK(TObjArray* track_list, 
00028                          Int_t trk_id,
00029                          const TVector3& ipv) {
00030   
00031   LCDTrack* trk=(LCDTrack*)(track_list->At(trk_id));
00032   Double_t* trk_parm=trk->GetTrackParameters();
00033   Double_t  emtrk[15];
00034   Int_t i=0,j=0,k=0;
00035   for(i=0; i < 5 ; i++) {
00036     for(j=0; j <= i ; j++ ) {
00037       emtrk[k++] = trk->GetErrorMatrixElement(i,j);
00038     }
00039   }
00040   Init(trk_id,trk_parm,emtrk,trk->GetCharge(),trk->GetMagneticField(),
00041              trk->GetParticle(),ipv);
00042 }
00043 
00044 //Initialize.
00045 void LCDVToplTRK::Init(){
00046   m_nvrt=0;
00047   Int_t i=0;
00048   for (i=0 ; i < 100 ; i++) {
00049     m_vertex_list[i]=0;
00050   }
00051   m_track_id=-1;
00052 
00053   for(i=0; i < 5 ; i++) {
00054     m_trackPar[i] = 0.0;
00055   }
00056   Int_t k=0;
00057   for(k=0; k < 15 ; k++) {
00058     m_trackErrorMatrix[k++] = 0;
00059   }
00060 
00061   m_charge=0;
00062   m_magneticfield=0;
00063   m_index = -1;
00064 
00065   m_trkval = 0;
00066   
00067   m_rmax=2.0;
00068   Int_t j;
00069   for (j=0 ; j < 3 ; j++) {
00070     m_err_tsi[j]=1E10;
00071     m_err_eta[j]=1E10;
00072   }
00073 
00074   m_v0flag=0;
00075   m_mcflav=0;
00076 }
00077 
00078 void LCDVToplTRK::Init(Int_t tkid,Double_t* tkpar, Double_t* e_tkpar, 
00079                        Double_t chrg, 
00080                        Double_t mfld,Int_t mc_ind,
00081                        const TVector3& ipv){
00082 
00083   Init();
00084   SetUp(tkid,tkpar,e_tkpar,chrg,mfld,mc_ind,ipv);
00085 }
00086 
00087 void LCDVToplTRK::SetUp(Int_t tkid, Double_t* tkpar, Double_t* e_tkpar, 
00088                         Double_t chrg,
00089                         Double_t mfld, Int_t mc_ind, const TVector3& ipv) {
00090   m_track_id=tkid;
00091   Int_t i=0;
00092   for(i=0; i < 5 ; i++) {
00093     m_trackPar[i] = tkpar[i];
00094   }
00095   Int_t k=0;
00096   for(k=0; k < 15 ; k++) {
00097     m_trackErrorMatrix[k] = e_tkpar[k];
00098   }
00099   m_charge       =chrg;
00100   m_magneticfield=mfld;
00101   m_index        =mc_ind;
00102 
00103   Double_t s1i=TMath::Sqrt(GetErrorTsi(0.0     ,ipv));
00104   Double_t s2i=TMath::Sqrt(GetErrorEta(0.0     ,ipv));
00105   Double_t s1m=TMath::Sqrt(GetErrorTsi(m_rmax/2.0,ipv));
00106   Double_t s2m=TMath::Sqrt(GetErrorEta(m_rmax/2.0,ipv));
00107   Double_t s1f=TMath::Sqrt(GetErrorTsi(m_rmax    ,ipv));
00108   Double_t s2f=TMath::Sqrt(GetErrorEta(m_rmax    ,ipv));
00109   m_err_tsi[0] = s1i;
00110   m_err_tsi[2] = (s1f + s1i)/2.0 - s1m;
00111   m_err_tsi[1] = (s1m - s1i) - m_err_tsi[2];
00112   m_err_eta[0] = s2i;
00113   m_err_eta[2] = (s2f + s2i)/2.0 - s2m;
00114   m_err_eta[1] = (s2m - s2i) - m_err_eta[2];
00115 }
00116 
00117 Double_t LCDVToplTRK::GetPositionX(Double_t l) {
00118   Double_t d0   = m_trackPar[0];
00119   Double_t phi0 = m_trackPar[1];
00120   Double_t omega= m_trackPar[2];
00121   Double_t phi  = phi0+omega*l;
00122 
00123   return 1.0/omega*TMath::Sin(phi) - (1.0/omega + d0)*TMath::Sin(phi0);
00124 }
00125 
00126 Double_t LCDVToplTRK::GetPositionY(Double_t l) {
00127   Double_t d0   = m_trackPar[0];
00128   Double_t phi0 = m_trackPar[1];
00129   Double_t omega= m_trackPar[2];
00130   Double_t phi  = phi0+omega*l;
00131 
00132   return -1.0/omega*TMath::Cos(phi) + (1.0/omega + d0)*TMath::Cos(phi0);
00133 }
00134 
00135 Double_t LCDVToplTRK::GetPositionZ(Double_t l) {
00136   Double_t z0   = m_trackPar[3];
00137   Double_t tanl = m_trackPar[4];
00138 
00139   return z0 + l*tanl;
00140 }
00141 
00142 TVector3 LCDVToplTRK::GetPositionVector(Double_t l) {
00143   return TVector3(GetPositionX(l),GetPositionY(l),GetPositionZ(l));
00144 }
00145 
00146 Double_t LCDVToplTRK::GetMomentumPx(Double_t l) {
00147   Double_t phi0 = m_trackPar[1];
00148   Double_t omega= m_trackPar[2];
00149   Double_t phi  = phi0+omega*l;
00150   Double_t pt   = m_magneticfield/333.567/TMath::Abs(m_trackPar[2]);
00151 
00152   return pt*TMath::Cos(phi);
00153 }
00154 
00155 Double_t LCDVToplTRK::GetMomentumPy(Double_t l) {
00156   Double_t phi0 = m_trackPar[1];
00157   Double_t omega= m_trackPar[2];
00158   Double_t phi  = phi0+omega*l;
00159   Double_t pt   = m_magneticfield/333.567/TMath::Abs(m_trackPar[2]);
00160 
00161   return pt*TMath::Sin(phi);
00162 }
00163 
00164 Double_t LCDVToplTRK::GetMomentumPz(Double_t l) {
00165   Double_t omega= m_trackPar[2];
00166   Double_t tanl = m_trackPar[4];
00167   Double_t pt   = m_magneticfield/333.567/TMath::Abs(omega);
00168 
00169   return pt*tanl;
00170 }
00171 
00172 TVector3 LCDVToplTRK::GetMomentumVector(Double_t l) {
00173   return TVector3(GetMomentumPx(l),GetMomentumPy(l),GetMomentumPz(l));
00174 }
00175 
00176 Double_t LCDVToplTRK::GetTsi(Double_t l, const TVector3& ipv) {
00177   Double_t bz    = GetMagneticField();
00178   Double_t phi0  = GetTrackParameter(1);
00179   Double_t omega = GetTrackParameter(2);
00180   Double_t phi   = phi0 + omega*l;
00181 
00182   TVector3 n(TMath::Sin(phi),-TMath::Cos(phi),0.0);
00183   if (bz < 0) { n*=-1; }
00184   TVector3 dx=GetPositionVector(l) - ipv;
00185 
00186   return n*dx;
00187 }
00188 
00189 Double_t LCDVToplTRK::GetEta(Double_t l, const TVector3& ipv) {
00190   Double_t bz    = GetMagneticField();
00191   Double_t phi0  = GetTrackParameter(1);
00192   Double_t omega = GetTrackParameter(2);
00193   Double_t tanl  = GetTrackParameter(4);
00194   Double_t cosl  = 1.0/TMath::Sqrt(1.0 + tanl*tanl);
00195   Double_t phi   = phi0 + omega*l;
00196 
00197   TVector3 q(-TMath::Cos(phi)*tanl*cosl,-TMath::Sin(phi)*tanl*cosl,cosl);
00198   if (bz < 0) { q *= -1; }
00199   TVector3 dx= GetPositionVector(l) - ipv;
00200 
00201   return q*dx;
00202 }
00203 
00204 Double_t LCDVToplTRK::GetTrackChiSqContrib(const TVector3& xvrt,
00205                                            const TVector3& ipv,
00206                                            Int_t f_debug) {
00207 
00208   Double_t l=0;
00209   TVector2 x0(GetPositionX(l)-ipv[0],GetPositionY(l)-ipv[1]);
00210   Double_t omega = GetTrackParameter(2);
00211   Double_t phi   = GetTrackParameter(1) - TMath::Pi()/2.0;
00212   TVector2 x0p   = x0.Rotate(-phi);
00213   TVector2 x0pr  = (xvrt-ipv).XYvector().Rotate(-phi);
00214   Double_t y0pr  = x0pr.Y() - x0p.Y();
00215   Double_t lp    = y0pr;
00216   // Double_t lp    = GetLPOCA2(xvrt);
00217 
00218   Double_t z0p   = GetPositionZ(l);
00219   Double_t tanl  = GetTrackParameter(4);
00220   Double_t s1s   = GetErrorTsiApprox(lp);
00221   Double_t s2s   = GetErrorEtaApprox(lp)*(1.0 + tanl*tanl);
00222 
00223   Double_t dx0pr = x0pr.X() - (x0p.X() - omega/2.0*y0pr*y0pr);
00224   Double_t dz0pr = xvrt[2]  - (z0p     + tanl     *lp       );
00225 
00226   if (f_debug) {
00227     printf("   x0pvrt =%f z0pvrt=%f\n",x0pr.X(),xvrt[2]);
00228     printf("   x0ptk  =%f z0ptk =%f\n",
00229            x0p.X() - omega/2.0*y0pr*y0pr,z0p+tanl*lp);
00230     printf("   x0p.X()=%f y0pr=%f omega=%f z0p=%f tanl=%f\n",
00231            x0p.X(),y0pr,omega,z0p,tanl);
00232   }
00233 
00234   return dx0pr/s1s*dx0pr + dz0pr/s2s*dz0pr;
00235 }
00236 
00237 Double_t LCDVToplTRK::GetTrackProbabirity(const TVector3& xvrt,
00238                                           const TVector3& ipv) {
00239 
00240   Double_t  u    = GetTrackChiSqContrib(xvrt,ipv);
00241   Double_t  prfu = 0.0;
00242   // protect against floating overflow
00243   if (u < 100.0) {
00244     prfu = TMath::Exp(-0.5*u);
00245   }
00246   return prfu;
00247 }
00248 
00249 Double_t LCDVToplTRK::GetDistance(Double_t l, const TVector3& xvrt) {
00250   return (GetPositionVector(l) - xvrt).Mag();
00251 }
00252 
00253 Double_t LCDVToplTRK::GetDistance(Double_t l, Double_t x, Double_t y, Double_t z) {
00254   TVector3 x0(x,y,z);
00255   return GetDistance(l,x0);
00256 }
00257 
00258 Double_t LCDVToplTRK::GetDistance2D(Double_t l, Double_t x, Double_t y) {
00259   Double_t x0=GetPositionX(l);
00260   Double_t y0=GetPositionY(l);
00261   return TMath::Sqrt((x0 - x)*(x0 - x) + (y0 - y)*(y0 - y));
00262 }
00263 
00264 Double_t LCDVToplTRK::GetTrdi(const TVector3& ipv, const TVector3& xvrt) {
00265   Double_t l0=GetLPOCA2(xvrt);
00266   Double_t l=GetLminTrdi(l0,ipv,xvrt);
00267   return GetTrdi(l,ipv,xvrt);
00268 }
00269 
00270 Double_t LCDVToplTRK::GetTrdi(Double_t l, const TVector3& ipv, const TVector3& xvrt) {
00271   TVector3 ipv2xvrt=xvrt - ipv;
00272   TVector3 ipv2xtrk=GetPositionVector(l) - ipv;
00273   return ipv2xtrk.Perp(ipv2xvrt);
00274 }
00275 
00276 Double_t LCDVToplTRK::GetLodi(const TVector3& ipv, const TVector3& xvrt) {
00277   Double_t l0=GetLPOCA2(xvrt);
00278   Double_t l=GetLminTrdi(l0,ipv,xvrt);
00279   return GetLodi(l,ipv,xvrt);
00280 }
00281 
00282 Double_t LCDVToplTRK::GetLodi(Double_t l, const TVector3& ipv, const TVector3& xvrt) {
00283   TVector3 ipv2xvrt=xvrt - ipv;
00284   TVector3 ipv2xtrk=GetPositionVector(l) - ipv;
00285   Double_t lodi=0;
00286   Double_t aipv2xvrt=ipv2xvrt.Mag();
00287   if (aipv2xvrt > 0) lodi=(ipv2xvrt*ipv2xtrk)/aipv2xvrt;
00288   return lodi;
00289 }
00290 
00291 void LCDVToplTRK::GetTrdiLodi(Double_t l, const TVector3& ipv, const TVector3& xvt,
00292                                   Double_t* trdi, Double_t* lodi) {
00293   *lodi=0;
00294   *trdi=0;
00295 
00296   TVector3 xtk=GetPositionVector(l);
00297   TVector3 ptk=GetMomentumVector(l);
00298 
00299   TVector3 vx=xvt - ipv;
00300 
00301   //Double_t anta=vx.Angle(ptk);
00302 
00303   TVector3 it2ip=xtk-ipv;
00304   Double_t aa=it2ip*vx;
00305   Double_t dd=it2ip*ptk;
00306   Double_t bb=ptk*vx;
00307   Double_t cc=-vx.Mag2();
00308   Double_t ee=ptk.Mag2();
00309 
00310   Double_t lambda=0;
00311   Double_t mu    =0;
00312   if (cc+bb*bb/ee == 0.0) { // no solution, track, axis parallel.
00313     *lodi  =-1000.0;
00314     lambda= 0.0;
00315     mu    = -dd/ee;
00316   } else {
00317     lambda = (dd*bb/ee - aa)/(cc + bb*bb/ee);
00318     mu     = (bb*lambda - dd)/ee;
00319   }
00320 
00321   // xu = point of normal hit on vertex axis
00322   // xt = point of normal hit on track vector
00323 
00324   TVector3 xu=ipv+lambda*vx ;
00325   TVector3 xt=xtk+mu    *ptk;
00326 
00327   TVector3 xu2ip=xu - ipv;
00328   TVector3 xu2xt=xu - xt;
00329   *lodi = xu2ip.Mag();
00330   //if we're behing IP...
00331   if (xu2ip*vx < 0.0) { *lodi=-*lodi; }
00332   *trdi=xu2xt.Mag();
00333 
00334   if (cc+bb*bb/ee == 0.0) {//no solution, track, axis parallel: reset lodi
00335     *lodi=-1000.0;
00336   }
00337 }
00338 
00339 void LCDVToplTRK::GetTrdiLodiAnta(Double_t l, const TVector3& ipv, const TVector3& xvt,
00340                                   Double_t* trdi, 
00341                                   Double_t* lodi, 
00342                                   Double_t* anta) {
00343   *lodi=0;
00344   *trdi=0;
00345   *anta=0;
00346 
00347   TVector3 xtk=GetPositionVector(l);
00348   TVector3 ptk=GetMomentumVector(l);
00349 
00350   TVector3 vx=xvt - ipv;
00351 
00352   *anta=vx.Angle(ptk);
00353 
00354   TVector3 it2ip=xtk-ipv;
00355   Double_t aa=it2ip*vx;
00356   Double_t dd=it2ip*ptk;
00357   Double_t bb=ptk*vx;
00358   Double_t cc=-vx.Mag2();
00359   Double_t ee=ptk.Mag2();
00360 
00361   Double_t lambda=0;
00362   Double_t mu    =0;
00363   if (cc+bb*bb/ee == 0.0) { // no solution, track, axis parallel.
00364     *lodi  =-1000.0;
00365     lambda= 0.0;
00366     mu    = -dd/ee;
00367   } else {
00368     lambda = (dd*bb/ee - aa)/(cc + bb*bb/ee);
00369     mu     = (bb*lambda - dd)/ee;
00370   }
00371 
00372   // xu = point of normal hit on vertex axis
00373   // xt = point of normal hit on track vector
00374 
00375   TVector3 xu=ipv+lambda*vx ;
00376   TVector3 xt=xtk+mu    *ptk;
00377 
00378   TVector3 xu2ip=xu - ipv;
00379   TVector3 xu2xt=xu - xt;
00380   *lodi = xu2ip.Mag();
00381   //if we're behing IP...
00382   if (xu2ip*vx < 0.0) { *lodi=-*lodi; }
00383   *trdi=xu2xt.Mag();
00384 
00385   if (cc+bb*bb/ee == 0.0) {//no solution, track, axis parallel: reset lodi
00386     *lodi=-1000.0;
00387   }
00388 }
00389 
00390 Double_t LCDVToplTRK::GetErrorMatrixElement(Int_t i, Int_t j) {
00391   Int_t ij[5][5]= { { 0, 1, 3, 6,10},
00392                     { 1, 2, 4, 7,11},
00393                     { 3, 4, 5, 8,12},
00394                     { 6, 7, 8, 9,13},
00395                     {10,11,12,13,14}};
00396 
00397   return m_trackErrorMatrix[ij[i][j]];
00398 }
00399 
00400 void LCDVToplTRK::GetErrorMatrix(TMatrixD& trackErrorMatrix){
00401   Int_t ij[5][5]= { { 0, 1, 3, 6,10},
00402                     { 1, 2, 4, 7,11},
00403                     { 3, 4, 5, 8,12},
00404                     { 6, 7, 8, 9,13},
00405                     {10,11,12,13,14}};
00406   Int_t ip=0,jp=0;
00407   for (ip=0 ; ip < 5 ; ip++) {
00408     for (jp=0 ; jp < 5 ; jp++) {
00409       trackErrorMatrix(ip,jp) = m_trackErrorMatrix[ij[ip][jp]];
00410     }
00411   }
00412 }
00413 
00414 void LCDVToplTRK::GetErrorMatrixPosition(Double_t l,TMatrixD& var2) {
00415   TMatrixD var0(5,5); GetErrorMatrix(var0);
00416   TMatrixD dxdaT(5,3); GetdxdaT(l,dxdaT);
00417   TMatrixD var1(var0,TMatrixD::kMult,dxdaT);
00418   TMatrixD dxda(TMatrixD::kTransposed,dxdaT);
00419   var2.Mult(dxda,var1);
00420 }
00421 
00422 void LCDVToplTRK::GetErrorMatrixMomentum(Double_t l,TMatrixD& var2) {
00423   TMatrixD var0(5,5); GetErrorMatrix(var0);
00424   TMatrixD dpdaT(5,3); GetdpdaT(l,dpdaT);
00425   TMatrixD var1(var0,TMatrixD::kMult,dpdaT);
00426   TMatrixD dpda(TMatrixD::kTransposed,dpdaT);
00427   var2.Mult(dpda,var1);
00428 }
00429 
00430 Double_t LCDVToplTRK::GetErrorTsi(Double_t l, const TVector3& ipv) {
00431   Double_t bz       = GetMagneticField();  
00432   if (bz != 0) bz /= TMath::Abs(bz);
00433   Double_t d0       = GetTrackParameter(0);
00434   Double_t phi0     = GetTrackParameter(1);
00435   Double_t omega    = GetTrackParameter(2);
00436   Double_t omegal   = omega*l;
00437   Double_t phi      = phi0 + omegal;
00438   Double_t csomegal = TMath::Cos(omegal);
00439   Double_t snomegal = TMath::Sin(omegal);
00440   Double_t csphi    = TMath::Cos(phi);
00441   Double_t snphi    = TMath::Sin(phi);
00442 
00443   TMatrixD dtsida(5,1);
00444   dtsida(0,0)=-bz*csomegal;
00445   dtsida(1,0)=-bz*(ipv[0]*csphi + ipv[1]*snphi);
00446   dtsida(2,0)=-bz*((csomegal - 1.0)/omega/omega + (1.0/omega + d0)*l*snomegal
00447                    - ipv[0]*l*csphi - ipv[1]*l*snphi);
00448   dtsida(3,0)= 0.0;
00449   dtsida(4,0)= 0.0;
00450   TMatrixD trkparmerror(5,5); GetErrorMatrix(trkparmerror);
00451   TMatrixD tmp0(trkparmerror,TMatrixD::kMult,dtsida);
00452   TMatrixD tmp1(dtsida,TMatrixD::kTransposeMult,tmp0);
00453 
00454   return tmp1(0,0);
00455 }
00456 
00457 Double_t LCDVToplTRK::GetErrorEta(Double_t l, const TVector3& ipv) {
00458   Double_t bz       = GetMagneticField(); 
00459   if (bz != 0) bz /= TMath::Abs(bz);
00460   Double_t d0       = GetTrackParameter(0);
00461   Double_t phi0     = GetTrackParameter(1);
00462   Double_t omega    = GetTrackParameter(2);
00463   Double_t z0       = GetTrackParameter(3);
00464   Double_t tanl     = GetTrackParameter(4);
00465   Double_t cosl     = 1.0/TMath::Sqrt(1.0 + tanl*tanl);
00466   Double_t sinl     = tanl*cosl;
00467   Double_t omegal   = omega*l;
00468   Double_t phi      = phi0 + omegal;
00469   Double_t csomegal = TMath::Cos(omegal);
00470   Double_t snomegal = TMath::Sin(omegal);
00471   Double_t csphi    = TMath::Cos(phi);
00472   Double_t snphi    = TMath::Sin(phi);
00473 
00474   TMatrixD detadaT(5,1);
00475   detadaT(0,0)=-bz*sinl*snomegal;
00476   detadaT(1,0)= bz*sinl*(-ipv[0]*snphi + ipv[1]*csphi);
00477   detadaT(2,0)= bz*sinl*(snomegal/omega/omega - (1.0/omega + d0)*l*csomegal
00478                          - ipv[0]*l*snphi - ipv[1]*l*csphi);
00479   detadaT(3,0)= bz*cosl;
00480   detadaT(4,0)= bz/cosl*(-(1.0/omega + d0)*snomegal + l + ipv[0]*csphi
00481                          +ipv[1]*snphi - (z0 - ipv[2])*tanl);
00482   TMatrixD trkparmerror(5,5); GetErrorMatrix(trkparmerror);
00483   TMatrixD tmp0(trkparmerror,TMatrixD::kMult,detadaT);
00484   TMatrixD tmp1(detadaT,TMatrixD::kTransposeMult,tmp0);
00485 
00486   return tmp1(0,0);
00487 }
00488 
00489 Double_t LCDVToplTRK::GetErrorTsiApprox(Double_t l) {
00490   Double_t ll=TMath::Min(l,m_rmax);
00491   ll=TMath::Max(ll,0.0);
00492   return TMath::Power(m_err_tsi[0] + m_err_tsi[1]*ll + m_err_tsi[2]*ll*ll,2);
00493 }
00494 
00495 Double_t LCDVToplTRK::GetErrorEtaApprox(Double_t l) {
00496   Double_t ll=TMath::Min(l,m_rmax);
00497   ll=TMath::Max(ll,0.0);
00498   return TMath::Power(m_err_eta[0] + m_err_eta[1]*ll + m_err_eta[2]*ll*ll,2);
00499 }
00500 
00501 Double_t LCDVToplTRK::GetErrorDistance(Double_t l, const TVector3& xvrt) {
00502     return GetErrorDistance(l,xvrt[0],xvrt[1],xvrt[2]);
00503 }
00504 
00505 Double_t LCDVToplTRK::GetErrorDistance(Double_t l,
00506                                        Double_t x, Double_t y, Double_t z) {
00507   TMatrixD dddaT(5,1); GetdddaT(l,x,y,z,dddaT);
00508   TMatrixD em(5,5); GetErrorMatrix(em);
00509   TMatrixD tmp0(em,TMatrixD::kMult,dddaT);
00510   TMatrixD tmp1(dddaT,TMatrixD::kTransposeMult,tmp0);
00511 
00512   return tmp1(0,0);  
00513 }
00514 
00515 Double_t LCDVToplTRK::GetErrorDistance(Double_t l, const TVector3& x_from, 
00516                                        TMatrixD& ex_from) {
00517   Double_t svrt=GetErrorDistance(l,x_from);
00518 
00519   TVector3 xvrt=GetPositionVector(l);
00520   TVector3 dx  =xvrt - x_from;
00521   Double_t d   =dx.Mag();
00522   TMatrixD dddxipT(3,1);
00523   if (d > 0) {
00524     dddxipT(0,0)=-dx[0]/d;
00525     dddxipT(1,0)=-dx[1]/d;
00526     dddxipT(2,0)=-dx[2]/d;
00527   } else {
00528     dddxipT(0,0)=0.0;
00529     dddxipT(1,0)=0.0;
00530     dddxipT(2,0)=0.0;
00531   }
00532   TMatrixD tmp0(ex_from,TMatrixD::kMult,dddxipT);
00533   TMatrixD tmp1(dddxipT,TMatrixD::kTransposeMult,tmp0);
00534 
00535   return tmp1(0,0)+svrt;
00536   
00537 }
00538 
00539 Double_t LCDVToplTRK::GetErrorDistance(Double_t l,
00540                                        Double_t x, Double_t y, Double_t z,
00541                                        TMatrixD& ex_from) {
00542   TVector3 xv(x,y,z);
00543   return GetErrorDistance(l,xv,ex_from);
00544 }
00545 
00546 Double_t LCDVToplTRK::GetErrorDistance2D(Double_t l, Double_t x, Double_t y) {
00547   TMatrixD dddaT(5,1); GetdddaT2D(l,x,y,dddaT);
00548   TMatrixD em(5,5); GetErrorMatrix(em);
00549   TMatrixD tmp0(em,TMatrixD::kMult,dddaT);
00550   TMatrixD tmp1(dddaT,TMatrixD::kTransposeMult,tmp0);
00551 
00552   return tmp1(0,0);  
00553 }
00554 
00555 Double_t LCDVToplTRK::GetErrorDistance2D(Double_t l, Double_t x, Double_t y,
00556                                          const TMatrixD& ex_from) {
00557   Double_t svrt=GetErrorDistance2D(l,x,y);
00558 
00559   TVector2 dx(GetPositionX(l) - x, GetPositionY(l) - y);
00560   Double_t d   =dx.Mod();
00561   TMatrixD dddxipT(3,1);
00562   if (d > 0) {
00563     dddxipT(0,0)=-dx.X()/d;
00564     dddxipT(1,0)=-dx.Y()/d;
00565     dddxipT(2,0)= 0.0;
00566   } else {
00567     dddxipT(0,0)=0.0;
00568     dddxipT(1,0)=0.0;
00569     dddxipT(2,0)=0.0;
00570   }
00571   TMatrixD tmp0(ex_from,TMatrixD::kMult,dddxipT);
00572   TMatrixD tmp1(dddxipT,TMatrixD::kTransposeMult,tmp0);
00573 
00574   return tmp1(0,0)+svrt;
00575 }
00576 
00577 Double_t LCDVToplTRK::GetErrorDistanceZ(Double_t l, Double_t z_from) {
00578   Double_t zp=GetPositionZ(l);
00579   TMatrixD dx(3,1);
00580   Double_t d=TMath::Abs(zp - z_from);
00581   dx(0,0) = 0.0;
00582   dx(1,0) = 0.0;
00583   if (d > 0) {
00584     dx(2,0) = -(z_from - zp)/d;
00585   } else {
00586     dx(2,0) = 0.0;
00587   }
00588   TMatrixD dxdaT(5,3); GetdxdaT(l,dxdaT);
00589   TMatrixD dddaT(dxdaT,TMatrixD::kMult,dx);
00590   TMatrixD em(5,5); GetErrorMatrix(em);
00591   TMatrixD tmp0(em,TMatrixD::kMult,dddaT);
00592   TMatrixD tmp1(dddaT,TMatrixD::kTransposeMult,tmp0);
00593   return tmp1(0,0);
00594 }
00595 
00596 Double_t LCDVToplTRK::GetErrorDistanceZ(Double_t l, Double_t z_from,
00597                                          const TMatrixD& ex_from) {
00598   Double_t svrt= GetErrorDistanceZ(l,z_from);
00599 
00600   Double_t zp  = GetPositionZ(l);
00601   Double_t d   = TMath::Abs(zp - z_from);
00602   TMatrixD dddxipT(3,1);
00603   if (d > 0) {
00604     dddxipT(0,0)= 0.0;
00605     dddxipT(1,0)= 0.0;
00606     dddxipT(2,0)=-(z_from - zp)/d;
00607   } else {
00608     dddxipT(0,0)=0.0;
00609     dddxipT(1,0)=0.0;
00610     dddxipT(2,0)=0.0;
00611   }
00612   TMatrixD tmp0(ex_from,TMatrixD::kMult,dddxipT);
00613   TMatrixD tmp1(dddxipT,TMatrixD::kTransposeMult,tmp0);
00614 
00615   return tmp1(0,0)+svrt;
00616 }
00617 
00618 void LCDVToplTRK::GetdxdaT(Double_t l, TMatrixD& dxdaT) {
00619   Double_t d0    = m_trackPar[0];
00620   Double_t phi0  = m_trackPar[1];
00621   Double_t omega = m_trackPar[2];
00622   Double_t phi   = phi0 + omega*l;
00623   Double_t csphi0= TMath::Cos(phi0);
00624   Double_t snphi0= TMath::Sin(phi0);
00625   Double_t csphi = TMath::Cos(phi);
00626   Double_t snphi = TMath::Sin(phi);
00627 
00628   dxdaT(0,0) =-snphi0;
00629   dxdaT(1,0) = 1.0/omega*csphi - (1.0/omega + d0)*csphi0;
00630   dxdaT(2,0) = (1.0/omega*(-snphi + snphi0) + csphi*l)/omega;
00631   dxdaT(3,0) = 0.0;
00632   dxdaT(4,0) = 0.0;
00633 
00634   dxdaT(0,1) = csphi0;
00635   dxdaT(1,1) = 1.0/omega*snphi - (1.0/omega + d0)*snphi0;
00636   dxdaT(2,1) = (1.0/omega*( csphi - csphi0) + snphi*l)/omega;
00637   dxdaT(3,1) = 0.0;
00638   dxdaT(4,1) = 0.0;
00639 
00640   dxdaT(0,2) = 0.0;
00641   dxdaT(1,2) = 0.0;
00642   dxdaT(2,2) = 0.0;
00643   dxdaT(3,2) = 1.0;
00644   dxdaT(4,2) = l;
00645 }
00646 
00647 void LCDVToplTRK::GetdpdaT(Double_t l, TMatrixD& dpdaT) {
00648   Double_t phi0  = m_trackPar[1];
00649   Double_t omega = m_trackPar[2];
00650   Double_t tanl  = m_trackPar[4];
00651   Double_t phi   = phi0 + omega*l;
00652   Double_t csphi = TMath::Cos(phi);
00653   Double_t snphi = TMath::Sin(phi);
00654   Double_t r2p   = TMath::Abs(m_magneticfield)/333.567;
00655 
00656   dpdaT(0,0) = 0.0;
00657   dpdaT(1,0) =-r2p/omega*snphi;
00658   dpdaT(2,0) =-r2p/omega/omega*csphi - r2p/omega*snphi*l;
00659   dpdaT(3,0) = 0.0;
00660   dpdaT(4,0) = 0.0;
00661 
00662   dpdaT(0,1) = 0.0;
00663   dpdaT(1,1) = r2p/omega*csphi;
00664   dpdaT(2,1) =-r2p/omega/omega*snphi + r2p/omega*csphi*l;
00665   dpdaT(3,1) = 0.0;
00666   dpdaT(4,1) = 0.0;
00667 
00668   dpdaT(0,2) = 0.0;
00669   dpdaT(1,2) = 0.0;
00670   dpdaT(2,2) =-r2p/omega/omega*tanl;
00671   dpdaT(3,2) = 0.0;
00672   dpdaT(4,2) = r2p/omega;
00673 }
00674 
00675 void LCDVToplTRK::GetdddaT(Double_t l, 
00676                            Double_t x, Double_t y, Double_t z,
00677                            TMatrixD& ddda) {
00678   Double_t xp=GetPositionX(l);
00679   Double_t yp=GetPositionY(l);
00680   Double_t zp=GetPositionZ(l);
00681   Double_t d =GetDistance(l,x,y,z);
00682   TMatrixD dx(3,1);
00683   dx(0,0) = -(x - xp)/d;
00684   dx(1,0) = -(y - yp)/d;
00685   dx(2,0) = -(z - zp)/d;
00686   TMatrixD dxdaT(5,3); GetdxdaT(l,dxdaT);
00687   ddda.Mult(dxdaT,dx);
00688 }
00689 
00690 void LCDVToplTRK::GetdddaT2D(Double_t l, 
00691                              Double_t x, Double_t y,
00692                              TMatrixD& ddda) {
00693   Double_t xp=GetPositionX(l);
00694   Double_t yp=GetPositionY(l);
00695   Double_t d =GetDistance2D(l,x,y);
00696   TMatrixD dx(3,1);
00697   dx(0,0) = -(x - xp)/d;
00698   dx(1,0) = -(y - yp)/d;
00699   dx(2,0) =  0.0;
00700   TMatrixD dxdaT(5,3); GetdxdaT(l,dxdaT);
00701   ddda.Mult(dxdaT,dx);
00702 }
00703 
00704 void LCDVToplTRK::Getd2xd2l(Double_t l, TMatrixD& d2xd2l) {
00705   Double_t phi0  = m_trackPar[1];
00706   Double_t omega = m_trackPar[2];
00707   Double_t phi   = phi0 + omega*l;
00708   
00709   d2xd2l(0,0)=-omega*TMath::Sin(phi);
00710   d2xd2l(1,0)= omega*TMath::Cos(phi);
00711   d2xd2l(2,0)= 0.0;
00712 }
00713 
00714 void LCDVToplTRK::Getdxdl(Double_t l,TMatrixD& dxdl) {
00715   Double_t phi0=GetTrackParameter(1);
00716   Double_t omega=GetTrackParameter(2);
00717   Double_t tanl=GetTrackParameter(4);
00718   Double_t phi=phi0+omega*l;
00719 
00720   dxdl(0,0)=TMath::Cos(phi);
00721   dxdl(1,0)=TMath::Sin(phi);
00722   dxdl(2,0)=tanl;
00723 }
00724 
00725 TVector3 LCDVToplTRK::GetdxdlVector(Double_t l) {
00726   Double_t phi0=GetTrackParameter(1);
00727   Double_t omega=GetTrackParameter(2);
00728   Double_t tanl=GetTrackParameter(4);
00729   Double_t phi=phi0+omega*l;
00730 
00731   return TVector3(TMath::Cos(phi),TMath::Sin(phi),tanl);
00732 }
00733 
00734 Double_t LCDVToplTRK::GetLPOCA2(const TVector3& x_given) {
00735   Double_t d0   =m_trackPar[0];
00736   Double_t phi0 =m_trackPar[1];
00737   Double_t omega=m_trackPar[2];
00738   TVector3 xc(-(d0+1.0/omega)*TMath::Sin(phi0),
00739               +(d0+1.0/omega)*TMath::Cos(phi0),
00740               0.0);
00741   TVector3 xxc=x_given-xc;
00742   xxc[2]=0.0;
00743   TVector3 nxxc=xxc.Unit();
00744   TVector3 nxc=-xc.Unit();
00745   TVector3 nxcXnxxc=nxc.Cross(nxxc);
00746   Double_t nxcXnxxcZ=nxcXnxxc[2];
00747   if      (nxcXnxxcZ < -1.0) { nxcXnxxcZ=-1.0; }
00748   else if (nxcXnxxcZ > +1.0) { nxcXnxxcZ=+1.0; }
00749 
00750   return TMath::ASin(nxcXnxxcZ)/omega;
00751 }
00752 
00753 // Get Closest Approach Point to beam axis...
00754 
00755 Double_t LCDVToplTRK::CalcPOCA3(Double_t l0,const TVector3& x_given) {
00756   const Int_t MAXTRY=50;
00757   const Double_t DLEN2TH=1E-10;
00758 
00759   Double_t l=l0;
00760 
00761   Double_t l_best=0.0,ddlen2dl_best=0.0,d2dlen2d2l_best=0.0;
00762   Double_t dlen2_c;
00763   Double_t dlen2_p=1e10;
00764   Double_t a;
00765   Double_t f=1.0;
00766   for (Int_t i=0; i < MAXTRY ; i++) {
00767     dlen2_c=dlen2(l,x_given);
00768     if ( TMath::Abs(dlen2_c - dlen2_p) < DLEN2TH && i != 0) {
00769       break;
00770     } else if ( i >=  MAXTRY-1 ) {
00771       l = l_best;
00772       break;
00773     } else if ( dlen2_c > dlen2_p && i != 0 ) {
00774       f /= 2.0;
00775       l  = l_best;
00776       l -=  ddlen2dl_best/d2dlen2d2l_best*f;
00777     } else {
00778       a= d2dlen2d2l(l,x_given);
00779       if (a != 0) {
00780         dlen2_p         = dlen2_c;
00781         f               = 1.0;
00782         l_best          = l;
00783         ddlen2dl_best   = ddlen2dl(l,x_given);
00784         d2dlen2d2l_best = a;
00785       } else {
00786         f /= 2.0;
00787         l  = l_best;
00788       }
00789       l -=  ddlen2dl_best/d2dlen2d2l_best*f;
00790     }
00791   }
00792 
00793   return l;
00794 }
00795 
00796 Double_t LCDVToplTRK::CalcPOCA3(const TVector3& x_given) {
00797   Double_t l=GetLPOCA2(x_given);
00798 
00799   TVector3 a=x_given-GetPositionVector(l);
00800   TVector3 b=GetMomentumVector(l);
00801   Double_t tanl=GetTrackParameter(4);
00802   l += a*b/b.Mag()*TMath::Sqrt(1.0/(1.0 + tanl*tanl));
00803 
00804   return CalcPOCA3(l,x_given);
00805 }
00806 
00807 //void LCDVToplTRK::AddVertex(TObject* vtx) {
00808 void LCDVToplTRK::AddVertex(LCDVToplVRT* vtx) {
00809   if (m_nvrt < 99) {
00810     m_vertex_list[m_nvrt++]=vtx;
00811   } else {
00812     fprintf(stderr,"LCDVToplTRK::AddVertex  # of vertices is over 100!\n");
00813   }
00814 }
00815 
00816 // void LCDVToplTRK::AddVertexAt(TObject* vtx,Int_t i) {
00817 void LCDVToplTRK::AddVertexAt(LCDVToplVRT* vtx,Int_t i) {
00818   if (i < m_nvrt+1) {
00819     if (m_vertex_list[i] == 0) m_nvrt++;
00820     m_vertex_list[i]=vtx;
00821   } else {
00822     fprintf(stderr,
00823             "LCDVToplTRK::AddVertex  vertex index is over 100! idx=%d\n nv=%d",
00824             i,m_nvrt);
00825   }
00826 }
00827 
00828 // void LCDVToplTRK::RemoveVertex(TObject* vrt){
00829 void LCDVToplTRK::RemoveVertex(LCDVToplVRT* vrt){
00830   for (Int_t i=0 ; i < m_nvrt ; i++) {
00831     if(m_vertex_list[i] == vrt) {
00832       m_vertex_list[i] = 0;
00833     }
00834   }
00835   CompressVerticesList();
00836 }
00837 
00838 void LCDVToplTRK::RemoveVertexAt(Int_t ivrt){
00839   if (ivrt < m_nvrt) {
00840     m_vertex_list[ivrt]=0;
00841     CompressVerticesList();
00842   } else {
00843     fprintf(stderr,
00844             "LCDVToplTRK::RemoveVertexAt ver. idx is over nv! idx=%d nv=%d\n",
00845             ivrt,m_nvrt);
00846   }
00847 }
00848 
00849 void LCDVToplTRK::CompressVerticesList() {
00850   Int_t nzero=0;
00851   for (Int_t i=0 ; i < m_nvrt ; i++) {
00852     if(m_vertex_list[i] == 0) {
00853       nzero++;
00854       continue;
00855     }
00856     if (nzero == 0) continue;
00857     m_vertex_list[i-nzero]=m_vertex_list[i];
00858   }
00859   m_nvrt -= nzero;
00860 }
00861 
00862 Double_t LCDVToplTRK::dlen2(Double_t l, const TVector3& x_given) {
00863   return (x_given - GetPositionVector(l)).Mag2();
00864 }
00865 
00866 Double_t LCDVToplTRK::ddlen2dl(Double_t l,const TVector3& x_given)
00867 {
00868   Double_t d0   =m_trackPar[0];
00869   Double_t phi0 =m_trackPar[1];
00870   Double_t omega=m_trackPar[2];
00871   Double_t z0   =m_trackPar[3];
00872   Double_t tanl =m_trackPar[4];
00873   Double_t x0p  =x_given[0];
00874   Double_t y0p  =x_given[1];
00875   Double_t z0p  =x_given[2];
00876 
00877   return 2*(d0 + 1.0/omega)*TMath::Sin(omega*l) + 2.0*z0*tanl + 2.0*l*tanl*tanl
00878     -2*x0p*TMath::Cos(phi0+omega*l) - 2*y0p*TMath::Sin(phi0+omega*l)-2*z0p*tanl;
00879 }
00880 
00881 Double_t LCDVToplTRK::d2dlen2d2l(Double_t l,const TVector3& x_given) 
00882 {
00883   Double_t d0   =m_trackPar[0];
00884   Double_t phi0 =m_trackPar[1];
00885   Double_t omega=m_trackPar[2];
00886   Double_t tanl =m_trackPar[4];
00887   Double_t x0p  =x_given[0];
00888   Double_t y0p  =x_given[1];
00889 
00890   return 2*omega*(d0 + 1.0/omega)*TMath::Cos(omega*l) + 2.0*tanl*tanl
00891     + 2*omega*x0p*TMath::Sin(phi0+omega*l)-2*omega*y0p*TMath::Cos(phi0+omega*l);
00892 }
00893 
00894 // Get Closest Approach Point to given axis...
00895 
00896 Double_t LCDVToplTRK::GetLminTrdiApprox(const TVector3& xfrom,const TVector3& xto) {
00897   Double_t l0=GetLPOCA2(xto);
00898   return GetLminTrdiApprox(l0,xfrom,xto);
00899 }
00900 
00901 Double_t LCDVToplTRK::GetLminTrdiApprox(Double_t l,
00902                                         const TVector3& ipv, 
00903                                         const TVector3& xvrt) {
00904   TVector3 ipv2xvrt=xvrt - ipv;
00905   TVector3 xtk=GetPositionVector(l);
00906   TVector3 ptk=GetMomentumVector(l);
00907   TVector3 ipv2xtk=xtk - ipv;
00908 
00909   Double_t a=ipv2xtk*ipv2xvrt;
00910   Double_t d=ipv2xtk*ptk;
00911   Double_t b=ptk*ipv2xvrt;
00912   Double_t c=-ipv2xvrt.Mag2();
00913   Double_t e= ptk.Mag2();
00914   Double_t lambda=0;
00915   Double_t mu=0;
00916 
00917   if (c+b*b/e == 0.0) { // no solution, track, axis parallel.
00918     lambda= 0.0;
00919     mu    = -d/e;
00920   } else {
00921     lambda = (d*b/e - a)/(c + b*b/e);
00922     mu     = (b*lambda - d)/e;
00923   }
00924 
00925   return l+mu*ptk.Pt();
00926 }
00927 
00928 Double_t LCDVToplTRK::GetLminTrdi(const TVector3& xfrom,const TVector3& xto) {
00929   Double_t l0=GetLPOCA2(xto);
00930   return GetLminTrdi(l0,xfrom,xto);
00931 }
00932 
00933 Double_t LCDVToplTRK::GetLminTrdi(Double_t l0,const TVector3& xfrom,const TVector3& xto) {
00934   const Int_t MAXTRY=50;
00935   const Double_t TRDITH=1E-10;
00936 
00937   Double_t l=l0;
00938 
00939   Double_t l_best=0.0;
00940   Double_t dtrdidl_best=0.0;
00941   Double_t d2trdid2l_best=0.0;
00942   Double_t trdi_c=1e10;
00943   Double_t trdi_p=1e10;
00944   Double_t a=0.0;
00945   Double_t f=1.0;
00946   for (Int_t i=0; i < MAXTRY ; i++) {
00947     trdi_c=TMath::Power(GetTrdi(l,xfrom,xto),2);
00948     if ( TMath::Abs(trdi_c - trdi_p) < TRDITH && i != 0) {
00949       break;
00950     } else if ( i >=  MAXTRY-1 ) {
00951       l = l_best;
00952       break;
00953     } else if ( trdi_c > trdi_p && i != 0 ) {
00954       f /= 2.0;
00955       l  = l_best;
00956       l -=  dtrdidl_best/d2trdid2l_best*f;
00957     } else {
00958       a= Getd2Trdi2d2l(l,xfrom,xto);
00959       if (a != 0) {
00960         trdi_p         = trdi_c;
00961         f              = 1.0;
00962         l_best         = l;
00963         dtrdidl_best   = GetdTrdi2dl(l,xfrom,xto);
00964         d2trdid2l_best = a;
00965       } else {
00966         f /= 2.0;
00967         l  = l_best;
00968       }
00969       l -=  dtrdidl_best/d2trdid2l_best*f;
00970     }
00971   }
00972 
00973   return l;
00974 }
00975 
00976 Double_t LCDVToplTRK::GetdTrdi2dl(Double_t l, const TVector3& xfrom, const TVector3& xto) {
00977   TVector3 n       = xto - xfrom; n.Unit();
00978   TVector3 dxdl    = GetdxdlVector(l);
00979   TVector3 x       = GetPositionVector(l);
00980   TVector3 xfrom2x = x - xfrom;
00981   return 2*(xfrom2x*dxdl) - 2*(x*n)*(n*dxdl);
00982 }
00983 
00984 Double_t LCDVToplTRK::Getd2Trdi2d2l(Double_t l, const TVector3& xfrom, const TVector3& xto) {
00985   Double_t d0      = GetTrackParameter(0);
00986   Double_t phi0    = GetTrackParameter(1);
00987   Double_t omega   = GetTrackParameter(2);
00988   Double_t tanl    = GetTrackParameter(4);
00989   Double_t omegal  = omega*l;
00990   Double_t phi     = phi0 + omegal;
00991   Double_t csphi   = TMath::Cos(phi);
00992   Double_t snphi   = TMath::Sin(phi);
00993   Double_t csomegal= TMath::Cos(omegal);
00994   Double_t snomegal= TMath::Sin(omegal);
00995   TVector3 n       = xto - xfrom; n.Unit();
00996   TVector3 x       = GetPositionVector(l);
00997   TVector3 xfrom2x = x - xfrom;
00998 
00999   return 2*(d0 + 1.0/omega)*omega*csomegal + 2*xfrom[0]*omega*snphi
01000     - 2*xfrom[1]*omega*csphi + 2*tanl*tanl 
01001     - 2*n[0]*n[0]*csphi*csphi + 2*n[0]*n[0]*omega*xfrom2x[0]*snphi 
01002     - 2*n[1]*n[1]*snphi*snphi - 2*n[1]*n[1]*omega*xfrom2x[1]*csphi
01003     - 2*n[2]*n[2]*tanl*tanl
01004     + 2*n[0]*n[1]*(d0 + 1.0/omega)*omega*snomegal
01005     + 2*n[0]*n[1]*(-xfrom[1]*snphi + xfrom[0]*snphi)
01006     - 2*n[0]*n[2]*tanl*csphi + 2*n[0]*n[2]*xfrom2x[2]*omega*snphi
01007     - 2*n[0]*n[2]*snphi*tanl
01008     - 2*n[1]*n[2]*tanl*snphi - 2*n[1]*n[2]*xfrom2x[2]*omega*csphi
01009     - 2*n[1]*n[2]*snphi*tanl;
01010 }
01011 
01012 Double_t LCDVToplTRK::GetSignedDistance(const TVector3& xfrom, 
01013                                         const TVector3& pjet) {
01014   Double_t l=CalcPOCA3(xfrom);
01015   return GetSignedDistance(l,xfrom,pjet);
01016 }
01017 
01018 Double_t LCDVToplTRK::GetSignedDistance(Double_t l, 
01019                                         const TVector3& xfrom, 
01020                                         const TVector3& pjet) {
01021   TVector3 dx(GetPositionVector(l) - xfrom);
01022   Double_t dlen=dx.Mag();
01023   return (dx*pjet > 0) ? dlen : -dlen;
01024 }
01025 
01026 Double_t LCDVToplTRK::GetSignedDistance2D(Double_t l, const TVector3& xfrom,
01027                                           const TVector3& pjet) {
01028   TVector3 dx(GetPositionVector(l) - xfrom);
01029   Double_t dlen=dx.Perp();
01030   return (dx*pjet > 0) ? dlen : -dlen;
01031 }
01032 
01033 Double_t LCDVToplTRK::GetDistanceNormByError(const TVector3& xfrom) {
01034   Double_t l=CalcPOCA3(xfrom);
01035   return GetDistanceNormByError(l,xfrom);
01036 }
01037 
01038 Double_t LCDVToplTRK::GetDistanceNormByError(Double_t l, const TVector3& xfrom) {
01039   Double_t dist=GetDistance(l,xfrom);
01040   Double_t sdst=TMath::Max(1e-20,TMath::Sqrt(GetErrorDistance(l,xfrom)));
01041   return dist/sdst;
01042 }
01043 
01044 Double_t LCDVToplTRK::GetDistanceNormByError(const TVector3& xfrom,
01045                                              TMatrixD& ex_from) {
01046   Double_t l=CalcPOCA3(xfrom);
01047   return GetDistanceNormByError(l,xfrom,ex_from);
01048 }
01049 
01050 Double_t LCDVToplTRK::GetDistanceNormByError(Double_t l, const TVector3& xfrom,
01051                                              TMatrixD& ex_from) {
01052   Double_t dist=GetDistance(l,xfrom);
01053   Double_t sdst=
01054     TMath::Max(1e-20,TMath::Sqrt(GetErrorDistance(l,xfrom,ex_from)));
01055   return dist/sdst;
01056 }
01057 
01058 Double_t LCDVToplTRK::GetSignedDistanceNormByError(const TVector3& xfrom,
01059                                                    const TVector3& pjet) {
01060   Double_t l=CalcPOCA3(xfrom);
01061   return GetSignedDistanceNormByError(l,xfrom,pjet);
01062 }
01063 
01064 Double_t LCDVToplTRK::GetSignedDistanceNormByError(Double_t l,
01065                                                    const TVector3& xfrom,
01066                                                    const TVector3& pjet) {
01067   Double_t dist=GetSignedDistance(l,xfrom,pjet);
01068   Double_t sdst=TMath::Max(1e-20,TMath::Sqrt(GetErrorDistance(l,xfrom)));
01069   return dist/sdst;
01070 }
01071 
01072 Double_t LCDVToplTRK::GetSignedDistanceNormByError(const TVector3& xfrom,
01073                                                    const TVector3& pjet,
01074                                                    TMatrixD& ex_from) {
01075   Double_t l=CalcPOCA3(xfrom);
01076   return GetSignedDistanceNormByError(l,xfrom,pjet,ex_from);
01077 }
01078 
01079 Double_t LCDVToplTRK::GetSignedDistanceNormByError(Double_t l, const TVector3& xfrom,
01080                                                    const TVector3& pjet,
01081                                                    TMatrixD& ex_from) {
01082   Double_t dist=GetSignedDistance(l,xfrom,pjet);
01083   Double_t sdst=
01084     TMath::Max(1e-20,TMath::Sqrt(GetErrorDistance(l,xfrom,ex_from)));
01085   return dist/sdst;
01086 }
01087 
01088 Double_t LCDVToplTRK::GetSignedDistance2DNormByError(Double_t l, 
01089                                                      const TVector3& xfrom,
01090                                                      const TVector3& pjet,
01091                                                      const TMatrixD& ex_from) {
01092   Double_t dist=GetSignedDistance2D(l,xfrom,pjet);
01093   Double_t sdst=
01094     TMath::Max(1e-20,TMath::Sqrt(GetErrorDistance2D(l,xfrom.X(),xfrom.Y(),
01095                                                     ex_from)));
01096   return dist/sdst;
01097 }
01098 

Generated on Tue Jul 18 18:33:59 2006 for LCDROOT by  doxygen 1.4.6