Reference Manual

Go to the documentation of this file.
00001 #include "TrackingTools/AnalyticalJacobians/interface/AnalyticalCurvilinearJacobian.h"
00002 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
00003 
00004 AnalyticalCurvilinearJacobian::AnalyticalCurvilinearJacobian 
00005 (const GlobalTrajectoryParameters& globalParameters,
00006  const GlobalPoint& x, 
00007  const GlobalVector& p, 
00008  const double& s) :  theJacobian(AlgebraicMatrixID())
00009 {
00010   //
00011   // helix: calculate full jacobian
00012   //
00013   if ( s*s*fabs(globalParameters.transverseCurvature())>1.e-5 ) { 
00014     GlobalPoint xStart = globalParameters.position();
00015     GlobalVector h  = globalParameters.magneticFieldInInverseGeV(xStart);
00016     computeFullJacobian(globalParameters,x,p,h,s);
00017   }
00018   //
00019   // straight line approximation, error in RPhi about 0.1um
00020   //
00021   else
00022     computeStraightLineJacobian(globalParameters,x,p,s);
00023    //dbg::dbg_trace(1,"ACJ1", globalParameters.vector(),x,p,s,theJacobian);
00024 }
00025 
00026 
00027 AnalyticalCurvilinearJacobian::AnalyticalCurvilinearJacobian
00028 (const GlobalTrajectoryParameters& globalParameters,
00029  const GlobalPoint& x, 
00030  const GlobalVector& p, 
00031  const GlobalVector& h, // h is the magnetic Field in Inverse GeV
00032  const double& s) :  theJacobian(AlgebraicMatrixID())
00033 {
00034   //
00035   // helix: calculate full jacobian
00036   //
00037   if ( s*s*fabs(globalParameters.transverseCurvature())>1.e-5 )
00038     computeFullJacobian(globalParameters,x,p,h,s);
00039   //
00040   // straight line approximation, error in RPhi about 0.1um
00041   //
00042   else
00043     computeStraightLineJacobian(globalParameters,x,p,s);
00044   
00045   //dbg::dbg_trace(1,"ACJ2", globalParameters.vector(),x,p,s,theJacobian);
00046 }
00047 
00048 // #ifdef TRPRFN_SSE
00049 #if defined(USE_SSEVECT) && !defined(TRPRFN_SCALAR)
00050 #include "AnalyticalCurvilinearJacobianSSE.icc"
00051 #else
00052 
00053 void
00054 AnalyticalCurvilinearJacobian::computeFullJacobian
00055 (const GlobalTrajectoryParameters& globalParameters,
00056  const GlobalPoint& x, 
00057  const GlobalVector& p, 
00058  const GlobalVector& h, 
00059  const double& s)
00060 {    
00061   //GlobalVector p1 = fts.momentum().unit();
00062   GlobalVector p1 = globalParameters.momentum().unit();
00063   GlobalVector p2 = p.unit();
00064   //GlobalPoint xStart = fts.position();
00065   GlobalPoint xStart = globalParameters.position();
00066   GlobalVector dx = xStart - x;
00067   //GlobalVector h  = MagneticField::inInverseGeV(xStart);
00068   // Martijn: field is now given as parameter.. GlobalVector h  = globalParameters.magneticFieldInInverseGeV(xStart);
00069 
00070   //double qbp = fts.signedInverseMomentum();
00071   double qbp = globalParameters.signedInverseMomentum();
00072   double absS = s;
00073   
00074   // calculate transport matrix
00075   // Origin: TRPRFN
00076   double t11 = p1.x(); double t12 = p1.y(); double t13 = p1.z();
00077   double t21 = p2.x(); double t22 = p2.y(); double t23 = p2.z();
00078   double cosl0 = p1.perp(); double cosl1 = 1./p2.perp();
00079   //AlgebraicMatrix a(5,5,1);
00080   // define average magnetic field and gradient 
00081   // at initial point - inlike TRPRFN
00082   GlobalVector hn = h.unit();
00083   double qp = -h.mag();
00084 //   double q = -h.mag()*qbp;
00085   double q = qp*qbp;
00086   double theta = q*absS; double sint = sin(theta); double cost = cos(theta);
00087   double hn1 = hn.x(); double hn2 = hn.y(); double hn3 = hn.z();
00088   double dx1 = dx.x(); double dx2 = dx.y(); double dx3 = dx.z();
00089   double gamma = hn1*t21 + hn2*t22 + hn3*t23;
00090   double an1 = hn2*t23 - hn3*t22;
00091   double an2 = hn3*t21 - hn1*t23;
00092   double an3 = hn1*t22 - hn2*t21;
00093   double au = 1./sqrt(t11*t11 + t12*t12);
00094   double u11 = -au*t12; double u12 = au*t11;
00095   double v11 = -t13*u12; double v12 = t13*u11; double v13 = t11*u12 - t12*u11;
00096   au = 1./sqrt(t21*t21 + t22*t22);
00097   double u21 = -au*t22; double u22 = au*t21;
00098   double v21 = -t23*u22; double v22 = t23*u21; double v23 = t21*u22 - t22*u21;
00099   // now prepare the transport matrix
00100   // pp only needed in high-p case (WA)
00101 //   double pp = 1./qbp;
00103 // moved up (where -h.mag() is needed()
00104 //   double qp = q*pp;
00105   double anv = -(hn1*u21 + hn2*u22          );
00106   double anu =  (hn1*v21 + hn2*v22 + hn3*v23);
00107   double omcost = 1. - cost; double tmsint = theta - sint;
00108   
00109   double hu1 =         - hn3*u12;
00110   double hu2 = hn3*u11;
00111   double hu3 = hn1*u12 - hn2*u11;
00112   
00113   double hv1 = hn2*v13 - hn3*v12;
00114   double hv2 = hn3*v11 - hn1*v13;
00115   double hv3 = hn1*v12 - hn2*v11;
00116   
00117   //   1/p - doesn't change since |p1| = |p2|
00118   
00119   //   lambda
00120   
00121   theJacobian(1,0) = -qp*anv*(t21*dx1 + t22*dx2 + t23*dx3);
00122   
00123   theJacobian(1,1) = cost*(v11*v21 + v12*v22 + v13*v23) +
00124                       sint*(hv1*v21 + hv2*v22 + hv3*v23) +
00125                     omcost*(hn1*v11 + hn2*v12 + hn3*v13) *
00126                            (hn1*v21 + hn2*v22 + hn3*v23) +
00127                 anv*(-sint*(v11*t21 + v12*t22 + v13*t23) +
00128                     omcost*(v11*an1 + v12*an2 + v13*an3) -
00129               tmsint*gamma*(hn1*v11 + hn2*v12 + hn3*v13) );
00130 
00131   theJacobian(1,2) = cost*(u11*v21 + u12*v22          ) +
00132                       sint*(hu1*v21 + hu2*v22 + hu3*v23) +
00133                     omcost*(hn1*u11 + hn2*u12          ) *
00134                            (hn1*v21 + hn2*v22 + hn3*v23) +
00135                 anv*(-sint*(u11*t21 + u12*t22          ) +
00136                     omcost*(u11*an1 + u12*an2          ) -
00137               tmsint*gamma*(hn1*u11 + hn2*u12          ) );
00138   theJacobian(1,2) *= cosl0;
00139 
00140   theJacobian(1,3) = -q*anv*(u11*t21 + u12*t22          );
00141 
00142   theJacobian(1,4) = -q*anv*(v11*t21 + v12*t22 + v13*t23);
00143 
00144   //   phi
00145 
00146   theJacobian(2,0) = -qp*anu*(t21*dx1 + t22*dx2 + t23*dx3)*cosl1;
00147 
00148   theJacobian(2,1) = cost*(v11*u21 + v12*u22          ) +
00149                       sint*(hv1*u21 + hv2*u22          ) +
00150                     omcost*(hn1*v11 + hn2*v12 + hn3*v13) *
00151                            (hn1*u21 + hn2*u22          ) +
00152                 anu*(-sint*(v11*t21 + v12*t22 + v13*t23) +
00153                     omcost*(v11*an1 + v12*an2 + v13*an3) -
00154               tmsint*gamma*(hn1*v11 + hn2*v12 + hn3*v13) );
00155   theJacobian(2,1) *= cosl1;
00156 
00157   theJacobian(2,2) = cost*(u11*u21 + u12*u22          ) +
00158                       sint*(hu1*u21 + hu2*u22          ) +
00159                     omcost*(hn1*u11 + hn2*u12          ) *
00160                            (hn1*u21 + hn2*u22          ) +
00161                 anu*(-sint*(u11*t21 + u12*t22          ) +
00162                     omcost*(u11*an1 + u12*an2          ) -
00163               tmsint*gamma*(hn1*u11 + hn2*u12          ) );
00164   theJacobian(2,2) *= cosl1*cosl0;
00165 
00166   theJacobian(2,3) = -q*anu*(u11*t21 + u12*t22          )*cosl1;
00167 
00168   theJacobian(2,4) = -q*anu*(v11*t21 + v12*t22 + v13*t23)*cosl1;
00169 
00170   //   yt
00171 
00172   //double cutCriterion = abs(s/fts.momentum().mag());
00173   double cutCriterion = fabs(s/globalParameters.momentum().mag());
00174   const double limit = 5.; // valid for propagations with effectively float precision
00175 
00176   if (cutCriterion > limit) {
00177     double pp = 1./qbp;
00178     theJacobian(3,0) = pp*(u21*dx1 + u22*dx2            );
00179     theJacobian(4,0) = pp*(v21*dx1 + v22*dx2 + v23*dx3);
00180   }
00181   else {
00182     double hp11 = hn2*t13 - hn3*t12;
00183     double hp12 = hn3*t11 - hn1*t13;
00184     double hp13 = hn1*t12 - hn2*t11;
00185     double temp1 = hp11*u21 + hp12*u22;
00186     double s2 = s*s;
00187     double secondOrder41 = 0.5 * qp * temp1 * s2;
00188     double ghnmp1 = gamma*hn1 - t11;
00189     double ghnmp2 = gamma*hn2 - t12;
00190     double ghnmp3 = gamma*hn3 - t13;
00191     double temp2 = ghnmp1*u21 + ghnmp2*u22;
00192     double s3 = s2 * s;
00193     double s4 = s3 * s;
00194     double h1 = h.mag();
00195     double h2 = h1 * h1;
00196     double h3 = h2 * h1;
00197     double qbp2 = qbp * qbp;
00198     //                           s*qp*s* (qp*s *qbp)
00199     double thirdOrder41 = 1./3 * h2 * s3 * qbp * temp2;
00200     //                           -qp * s * qbp  * above
00201     double fourthOrder41 = 1./8 * h3 * s4 * qbp2 * temp1;
00202     theJacobian(3,0) = secondOrder41 + (thirdOrder41 + fourthOrder41);
00203 
00204     double temp3 = hp11*v21 + hp12*v22 + hp13*v23;
00205     double secondOrder51 = 0.5 * qp * temp3 * s2;
00206     double temp4 = ghnmp1*v21 + ghnmp2*v22 + ghnmp3*v23;
00207     double thirdOrder51 = 1./3 * h2 * s3 * qbp * temp4;
00208     double fourthOrder51 = 1./8 * h3 * s4 * qbp2 * temp3;
00209     theJacobian(4,0) = secondOrder51 + (thirdOrder51 + fourthOrder51);
00210   }
00211 
00212   theJacobian(3,1) = (sint*(v11*u21 + v12*u22          ) +
00213                      omcost*(hv1*u21 + hv2*u22          ) +
00214                      tmsint*(hn1*u21 + hn2*u22          ) *
00215                             (hn1*v11 + hn2*v12 + hn3*v13))/q;
00216 
00217   theJacobian(3,2) = (sint*(u11*u21 + u12*u22          ) +
00218                      omcost*(hu1*u21 + hu2*u22          ) +
00219                      tmsint*(hn1*u21 + hn2*u22          ) *
00220                             (hn1*u11 + hn2*u12          ))*cosl0/q;
00221 
00222   theJacobian(3,3) = (u11*u21 + u12*u22          );
00223   
00224   theJacobian(3,4) = (v11*u21 + v12*u22          );
00225 
00226   //   zt
00227 
00228   theJacobian(4,1) = (sint*(v11*v21 + v12*v22 + v13*v23) +
00229                      omcost*(hv1*v21 + hv2*v22 + hv3*v23) +
00230                      tmsint*(hn1*v21 + hn2*v22 + hn3*v23) *
00231                             (hn1*v11 + hn2*v12 + hn3*v13))/q;
00232 
00233   theJacobian(4,2) = (sint*(u11*v21 + u12*v22          ) +
00234                      omcost*(hu1*v21 + hu2*v22 + hu3*v23) +
00235                      tmsint*(hn1*v21 + hn2*v22 + hn3*v23) *
00236                             (hn1*u11 + hn2*u12          ))*cosl0/q;
00237 
00238   theJacobian(4,3) = (u11*v21 + u12*v22          );
00239 
00240   theJacobian(4,4) = (v11*v21 + v12*v22 + v13*v23);
00241   // end of TRPRFN
00242 }
00243 
00244 #endif
00245 
00246 void AnalyticalCurvilinearJacobian::computeInfinitesimalJacobian 
00247 (const GlobalTrajectoryParameters& globalParameters,
00248  const GlobalPoint&, 
00249  const GlobalVector& p, 
00250  const GlobalVector& h, 
00251  const double& s) {
00252   /*
00253    * origin  TRPROP
00254    *
00255    C *** ERROR PROPAGATION ALONG A PARTICLE TRAJECTORY IN A MAGNETIC FIELD
00256    C     ROUTINE ASSUMES THAT IN THE INTERVAL (X1,X2) THE QUANTITIES 1/P
00257    C     AND (HX,HY,HZ) ARE RATHER CONSTANT. DELTA(PHI) MUST NOT BE TOO LARGE
00258    C
00259    C     Authors: A. Haas and W. Wittek
00260    C
00261    
00262   */
00263   
00264   
00265   double qbp = globalParameters.signedInverseMomentum();
00266   double absS = s;
00267   
00268   // average momentum
00269   GlobalVector tn = (globalParameters.momentum()+p).unit(); 
00270   double sinl = tn.z(); 
00271   double cosl = std::sqrt(1.-sinl*sinl); 
00272   double cosl1 = 1./cosl;
00273   double tgl=sinl*cosl1;
00274   double sinp = tn.y()*cosl1;
00275   double cosp = tn.x()*cosl1;
00276 
00277   // define average magnetic field and gradient 
00278   // at initial point - inlike TRPROP
00279   double b0= h.x()*cosp+h.y()*sinp;
00280   double b2=-h.x()*sinp+h.y()*cosp;
00281   double b3=-b0*sinl+h.z()*cosl;
00282 
00283   theJacobian(3,2)=absS*cosl;
00284   theJacobian(4,1)=absS;
00285 
00286 
00287   theJacobian(1,0) =  absS*b2;
00288   //if ( qbp<0) theJacobian(1,0) = -theJacobian(1,0);
00289   theJacobian(1,2) = -b0*(absS*qbp);
00290   theJacobian(1,3) =  b3*(b2*qbp*(absS*qbp));
00291   theJacobian(1,4) = -b2*(b2*qbp*(absS*qbp));
00292   
00293   theJacobian(2,0) = -absS*b3*cosl1;
00294   // if ( qbp<0) theJacobian(2,0) = -theJacobian(2,0);
00295   theJacobian(2,1) = b0*(absS*qbp)*cosl1*cosl1;
00296   theJacobian(2,2) = 1.+tgl*b2*(absS*qbp);
00297   theJacobian(2,3) = -b3*(b3*qbp*(absS*qbp)*cosl1);
00298   theJacobian(2,4) =  b2*(b3*qbp*(absS*qbp)*cosl1);
00299   
00300   theJacobian(3,4) = -b3*tgl*(absS*qbp);
00301   theJacobian(4,3) =  b3*tgl*(absS*qbp);
00302 
00303 
00304 }
00305 
00306 void
00307 AnalyticalCurvilinearJacobian::computeStraightLineJacobian
00308 (const GlobalTrajectoryParameters& globalParameters,
00309  const GlobalPoint& x, const GlobalVector& p, const double& s)
00310 {
00311   //
00312   // matrix: elements =1 on diagonal and =0 are already set
00313   // in initialisation
00314   //
00315   GlobalVector p1 = globalParameters.momentum().unit();
00316   double cosl0 = p1.perp();
00317   theJacobian(3,2) = cosl0 * s;
00318   theJacobian(4,1) = s;
00319 }
/data/refman/pasoursint/CMSSW_6_1_2_SLHC4_patch1/src/TrackingTools/AnalyticalJacobians/src/AnalyticalCurvilinearJacobian.cc
