00001 #include <assert.h>
00002 #include <math.h>
00003
00004 #include "MdcGeom/Constants.h"
00005 #include "MdcGeom/BesAngle.h"
00006 #include "CLHEP/Geometry/Point3D.h"
00007 #include "CLHEP/Vector/ThreeVector.h"
00008 #include "CLHEP/Matrix/SymMatrix.h"
00009 #include "TrkBase/NeutTraj.h"
00010 #include "TrkBase/TrkVisitor.h"
00011 #include "MdcRecoUtil/DifNumber.h"
00012 #include "MdcRecoUtil/DifPoint.h"
00013 #include "MdcRecoUtil/DifVector.h"
00014
00015
00016 NeutTraj::NeutTraj(const NeutParams& np, double lowlim,double hilim,
00017 const HepPoint3D& refpoint) :
00018 TrkSimpTraj(np.parameter(), np.covariance(), lowlim,hilim)
00019 {
00020 }
00021
00022 NeutTraj::NeutTraj( const NeutTraj& n )
00023 : TrkSimpTraj(n.parameters()->parameter(), n.parameters()->covariance(),
00024 n.lowRange(),n.hiRange())
00025 {
00026 }
00027
00028
00029 NeutTraj*
00030 NeutTraj::clone() const
00031 {
00032 return new NeutTraj(*this);
00033 }
00034
00035 NeutTraj&
00036 NeutTraj::operator=(const NeutTraj& n)
00037 {
00038 if( &n != this ){
00039 for(int iend=0;iend<2;iend++)
00040 flightrange[iend] = n.flightrange[iend];
00041 _dtparams = *n._np();
00042 }
00043 return *this;
00044 }
00045
00046
00047 NeutTraj::~NeutTraj() {
00048 }
00049
00050
00051 double
00052 NeutTraj::x( const double& f ) const {
00053
00054 return cosDip()*f*cos(phi0()) - d0()*sin(phi0());
00055 }
00056
00057
00058 double
00059 NeutTraj::y( const double& f ) const {
00060
00061 return cosDip()*f*sin(phi0()) + d0()*cos(phi0());
00062 }
00063
00064
00065 double
00066 NeutTraj::z(const double& f) const {
00067
00068 return (z0() + f*sinDip());
00069 }
00070
00071
00072 HepPoint3D
00073 NeutTraj::position( double f) const {
00074
00075 double cdip = cosDip();
00076 double sphi0 = sin(phi0());
00077 double cphi0 = cos(phi0());
00078
00079 return HepPoint3D(cdip * f * cphi0 - d0()*sphi0,
00080 cdip * f * sphi0 + d0()*cphi0,
00081 z0() + f * tanDip() * cdip);
00082 }
00083
00084
00085 Hep3Vector
00086 NeutTraj::direction( double f) const {
00087
00088 return Hep3Vector ( cos(phi0())*cosDip(), sin(phi0())*cosDip(),
00089 sinDip() );
00090 }
00091
00092 Hep3Vector
00093 NeutTraj::delDirect( double fltLen ) const{
00094
00095 return Hep3Vector(0.0, 0.0, 0.0);
00096 }
00097
00098
00099 double
00100 NeutTraj::distTo1stError(double , double tol, int) const {
00101
00102 return 9999.;
00103 }
00104
00105
00106 double
00107 NeutTraj::distTo2ndError(double , double tol, int) const {
00108
00109 return 9999.;
00110 }
00111
00112 void
00113 NeutTraj::getInfo(double fltLen, HepPoint3D& pos, Hep3Vector& dir,
00114 Hep3Vector& delDir) const
00115 {
00116
00117 pos = position(fltLen);
00118 dir = direction(fltLen);
00119 delDir = delDirect(fltLen);
00120 }
00121
00122 void
00123 NeutTraj::getInfo( double fltLen, HepPoint3D& pos, Hep3Vector& dir ) const
00124 {
00125
00126 pos = position(fltLen);
00127 dir = direction(fltLen);
00128 return;
00129 }
00130
00131 void
00132 NeutTraj::getDFInfo(double flt, DifPoint& pos, DifVector& dir,
00133 DifVector& delDir) const
00134 {
00135
00136
00137
00138
00139
00140
00141 DifNumber phi0Df(phi0(), NeutParams::_phi0+1, NeutParams::_nneutprm);
00142 DifNumber d0Df(d0(), NeutParams::_d0+1, NeutParams::_nneutprm);
00143 DifNumber z0Df(z0(), NeutParams::_z0+1, NeutParams::_nneutprm);
00144 DifNumber tanDipDf(tanDip(), NeutParams::_tanDip+1, NeutParams::_nneutprm);
00145 DifNumber pDf(p(), NeutParams::_p+1, NeutParams::_nneutprm);
00146
00147
00148
00149 DifNumber s0Df(flt, NeutParams::_s0+1, NeutParams::_nneutprm);
00150 phi0Df.setIndepPar( parameters() );
00151 d0Df.setIndepPar( parameters() );
00152 z0Df.setIndepPar( parameters() );
00153 tanDipDf.setIndepPar( parameters() );
00154 pDf.setIndepPar( parameters() );
00155 s0Df.setIndepPar( parameters() );
00156 DifNumber dipDf = atan(tanDipDf);
00157
00158 DifNumber sDip, cDip;
00159 dipDf.cosAndSin(cDip, sDip);
00160 DifNumber sinPhi0, cosPhi0;
00161 phi0Df.cosAndSin(cosPhi0, sinPhi0);
00162
00163 DifNumber x = cDip*s0Df*cosPhi0 - d0Df*sinPhi0;
00164 DifNumber y = cDip*s0Df*sinPhi0 + d0Df*cosPhi0;
00165 DifNumber z = z0Df + s0Df*sDip ;
00166 pos = DifPoint(x, y, z);
00167 dir = DifVector( cosPhi0*cDip, sinPhi0*cDip, sDip );
00168 delDir = DifVector(0., 0., 0.);
00169 }
00170
00171 HepMatrix
00172 NeutTraj::derivDeflect(double fltlen,deflectDirection idirect) const
00173 {
00174
00175
00176
00177
00178
00179
00180
00181
00182 HepMatrix ddflct(NeutParams::_nneutprm,1);
00183
00184
00185
00186 double tand = tanDip();
00187 double cosd = cosDip();
00188
00189
00190
00191 switch (idirect) {
00192 case theta1:
00193 ddflct[NeutParams::_p][0] = 0.;
00194 ddflct[NeutParams::_tanDip][0] = 1.0/(cosd*cosd);
00195 ddflct[NeutParams::_d0][0] = 0.;
00196 ddflct[NeutParams::_phi0][0] = 0.;
00197 ddflct[NeutParams::_s0][0] = 0.;
00198 ddflct[NeutParams::_z0][0] = translen(fltlen) * (-1. - (tand*tand));
00199 break;
00200 case theta2:
00201 ddflct[NeutParams::_p][0] = 0;
00202 ddflct[NeutParams::_tanDip][0] = 0;
00203 ddflct[NeutParams::_d0][0] = -translen(fltlen)/cosd;
00204 ddflct[NeutParams::_phi0][0] = 1./cosd;
00205 ddflct[NeutParams::_z0][0] = 0.;
00206 ddflct[NeutParams::_s0][0] = 0.;
00207 break;
00208 }
00209 return ddflct;
00210 }
00211
00212 HepMatrix
00213 NeutTraj::derivDisplace(double fltlen,deflectDirection idirect) const
00214 {
00215
00216
00217
00218
00219
00220
00221
00222
00223 HepMatrix ddflct(NeutParams::_nneutprm,1);
00224
00225
00226
00227 double cosd = cosDip();
00228
00229
00230
00231 switch (idirect) {
00232 case theta1:
00233 ddflct[NeutParams::_p][0] = 0.;
00234 ddflct[NeutParams::_tanDip][0] = 0.0;
00235 ddflct[NeutParams::_d0][0] = 0.;
00236 ddflct[NeutParams::_phi0][0] = 0.;
00237 ddflct[NeutParams::_s0][0] = 0.;
00238 ddflct[NeutParams::_z0][0] = 1.0/cosd;
00239 break;
00240 case theta2:
00241 ddflct[NeutParams::_p][0] = 0;
00242 ddflct[NeutParams::_tanDip][0] = 0;
00243 ddflct[NeutParams::_d0][0] = 1.0;
00244 ddflct[NeutParams::_phi0][0] = 0.0;
00245 ddflct[NeutParams::_z0][0] = 0.;
00246 ddflct[NeutParams::_s0][0] = 0.;
00247 break;
00248 }
00249 return ddflct;
00250 }
00251
00252 HepMatrix
00253 NeutTraj::derivPFract(double fltlen) const
00254 {
00255
00256
00257
00258
00259
00260
00261
00262
00263 HepMatrix dmomfrac(NeutParams::_nneutprm,1);
00264
00265
00266 dmomfrac[NeutParams::_p][0] = p();
00267 dmomfrac[NeutParams::_tanDip][0] = 0.0;
00268 dmomfrac[NeutParams::_d0][0] = 0.0;
00269 dmomfrac[NeutParams::_phi0][0] = 0.0;
00270 dmomfrac[NeutParams::_z0][0] = 0.0;
00271 dmomfrac[NeutParams::_s0][0] = 0.0;
00272 return dmomfrac;
00273 }
00274
00275
00276 double
00277 NeutTraj::phi0() const
00278 {
00279 return BesAngle(_np()->phi0()).rad();
00280 }
00281
00282 void
00283 NeutTraj::paramFunc(const HepPoint3D& oldpoint,const HepPoint3D& newpoint,
00284 const HepVector& oldvec,const HepSymMatrix& oldcov,
00285 HepVector& newvec,HepSymMatrix& newcov,
00286 double fltlen)
00287 {
00288
00289
00290 HepVector parvec(oldvec);
00291
00292 newvec = parvec;
00293
00294 double delx = newpoint.x()-oldpoint.x();
00295 double dely = newpoint.y()-oldpoint.y();
00296 double delz = newpoint.z()-oldpoint.z();
00297
00298 double cos0 = cos(parvec[NeutParams::_phi0]);
00299 double sin0 = sin(parvec[NeutParams::_phi0]);
00300 double perp = delx*sin0-dely*cos0;
00301 double para = delx*cos0+dely*sin0;
00302 double tand = parvec[NeutParams::_tanDip];
00303
00304
00305
00306 newvec[NeutParams::_d0] = parvec[NeutParams::_d0] + perp;
00307
00308 newvec[NeutParams::_phi0] = parvec[NeutParams::_phi0];
00309
00310 newvec[NeutParams::_z0] +=
00311 parvec[NeutParams::_tanDip]*(delx/cos0)*(1.+(sin0/cos0)) - delz;
00312
00313 HepMatrix covrot(NeutParams::_nneutprm,NeutParams::_nneutprm,0);
00314
00315
00316 covrot[NeutParams::_p][NeutParams::_p] = 1.0;
00317
00318 covrot[NeutParams::_tanDip][NeutParams::_tanDip] = 1.0;
00319
00320 covrot[ NeutParams::_d0][ NeutParams::_d0] = 1.;
00321 covrot[ NeutParams::_d0][ NeutParams::_phi0] = para;
00322
00323 covrot[ NeutParams::_phi0][ NeutParams::_p] = para;
00324 covrot[ NeutParams::_phi0][ NeutParams::_phi0] = 1.;
00325
00326 covrot[ NeutParams::_z0][ NeutParams::_phi0] = tand*-2.*perp;
00327 covrot[ NeutParams::_z0][ NeutParams::_tanDip] =
00328 (delx/cos0)*(1.+(sin0/cos0));
00329 covrot[ NeutParams::_z0][ NeutParams::_z0] = 1.0;
00330 covrot[ NeutParams::_s0][ NeutParams::_s0] = 1.0;
00331
00332
00333 newcov = oldcov.similarity(covrot);
00334
00335 }
00336
00337 void
00338 NeutTraj::invertParams(TrkParams* newparams, std::vector<bool>& flags) const
00339 {
00340 assert(1==0);
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364 }
00365
00366 void
00367 NeutTraj::visitAccept(TrkVisitor* vis) const
00368 {
00369
00370 vis->trkVisitNeutTraj(this);
00371 }