Reference Manual

Go to the documentation of this file.
00001 #include "RecoVertex/KinematicFit/interface/SimplePointingConstraint.h"
00002 #include "RecoVertex/VertexPrimitives/interface/VertexException.h"
00003 
00004 
00005 std::pair<AlgebraicVector, AlgebraicVector> SimplePointingConstraint::value(const AlgebraicVector& exPoint) const
00006 {
00007  if(exPoint.num_row() ==0 ) throw VertexException("PointingKinematicConstraint::value requested for zero Linearization point");
00008 
00009 //security check for extended cartesian parametrization 
00010  int inSize = exPoint.num_row(); 
00011  if((inSize%7) !=0) throw VertexException("PointingKinematicConstraint::linearization point has a wrong dimension");
00012  int nStates = inSize/7;
00013  if(nStates != 1) throw VertexException("PointingKinematicConstraint::Current version does not support the multistate refit");
00014  
00015  AlgebraicVector lPar = exPoint;
00016  AlgebraicVector vl(2,0);
00017  
00018 //vector of values 1x2  for given particle
00019  AlgebraicVector lValue = makeValue(lPar).first;
00020  vl(1) =lValue(1);
00021  vl(2) =lValue(2);
00022  return std::pair<AlgebraicVector,AlgebraicVector>(vl,lPar); 
00023 }
00024 
00025 std::pair<AlgebraicMatrix, AlgebraicVector> SimplePointingConstraint::derivative(const AlgebraicVector& exPoint) const
00026 {
00027  if(exPoint.num_row() ==0 ) throw VertexException("PointingKinematicConstraint::value requested for zero Linearization point");
00028 
00029 //security check for extended cartesian parametrization 
00030  int inSize = exPoint.num_row(); 
00031  if((inSize%7) !=0) throw VertexException("PointingKinematicConstraint::linearization point has a wrong dimension");
00032  int nStates = inSize/7;
00033  if(nStates != 1) throw VertexException("PointingKinematicConstraint::Current version does not support the multistate refit");
00034  AlgebraicVector lPar = exPoint;
00035 
00036 //2x7 derivative matrix for given particle
00037  AlgebraicMatrix lDeriv = makeDerivative(lPar).first;
00038  AlgebraicMatrix dr(2,7,0);
00039  dr.sub(1,1,lDeriv);
00040  return std::pair<AlgebraicMatrix,AlgebraicVector>(dr,lPar);
00041 }
00042 
00043 std::pair<AlgebraicMatrix, AlgebraicVector> SimplePointingConstraint::derivative(const std::vector<RefCountedKinematicParticle> par) const
00044 {
00045  int nStates = par.size();
00046  if(nStates == 0) throw VertexException("PointingKinematicConstraint::Empty vector of particles passed");
00047  if(nStates != 1) throw VertexException("PointingKinematicConstraint::Current version does not support the multistate refit");
00048  
00049  AlgebraicMatrix dr(2,7,0);
00050  AlgebraicVector lPoint = asHepVector<7>(par.front()->currentState().kinematicParameters().vector());
00051 
00052 //2x7 derivative matrix for given state  
00053  AlgebraicMatrix lDeriv = makeDerivative(lPoint).first;
00054  dr.sub(1,1,lDeriv);
00055 // cout<<"Derivative returned: "<<dr<<endl;
00056 // cout<<"For the value: "<<lPoint<<endl;
00057  return std::pair<AlgebraicMatrix,AlgebraicVector>(dr,lPoint);
00058 }
00059 
00060 std::pair<AlgebraicVector, AlgebraicVector> SimplePointingConstraint::value(const std::vector<RefCountedKinematicParticle> par) const
00061 { 
00062  int nStates = par.size();
00063  if(nStates == 0) throw VertexException("PointingKinematicConstraint::Empty vector of particles passed");
00064  if(nStates != 1) throw VertexException("PointingKinematicConstraint::Current version does not support the multistate refit");
00065  AlgebraicVector vl(2,0);
00066  AlgebraicVector lPoint = asHepVector<7>(par.front()->currentState().kinematicParameters().vector());
00067  vl(1) = makeValue(lPoint).first(1);
00068  vl(2) = makeValue(lPoint).first(2);
00069 // cout<<"Value returned: "<<vl<<endl;
00070 // cout<<"For the point: "<<lPoint<<endl;
00071  
00072  return std::pair<AlgebraicVector,AlgebraicVector>(vl,lPoint);
00073 }
00074  
00075 AlgebraicVector SimplePointingConstraint::deviations(int nStates) const
00076 {return AlgebraicVector(7*nStates,0);}
00077  
00078 int SimplePointingConstraint::numberOfEquations() const
00079 {return 2;}
00080  
00081 std::pair<AlgebraicVector,AlgebraicVector> SimplePointingConstraint::makeValue(const AlgebraicVector& exPoint)const 
00082 { 
00083 // cout<<"Make value called"<<endl;
00084  AlgebraicVector vl(2,0);
00085  AlgebraicVector point = exPoint;
00086  double dx = point(1) - refPoint.x();
00087  double dy = point(2) - refPoint.y();
00088  double dz = point(3) - refPoint.z();
00089  double px = point(4);
00090  double py = point(5); 
00091  double pz = point(6);
00092 
00093 
00094 //half angle solution: sin((alpha - betha)/2)
00095  double cos_phi_p = px/sqrt(px*px + py*py);
00096  double cos_phi_x = dx/sqrt(dx*dx + dy*dy);
00097 // cout<<"mom cos phi"<<cos_phi_p<<endl;
00098 // cout<<"x cos phi"<<cos_phi_x<<endl;
00099  
00100  
00101  double cos_theta_p = sqrt(px*px + py*py)/sqrt(px*px + py*py + pz*pz); 
00102  double cos_theta_x = sqrt(dx*dx + dy*dy)/sqrt(dx*dx + dy*dy + dz*dz);
00103  
00104  float feq = sqrt((1-cos_phi_p)*(1+cos_phi_x)) - sqrt((1+cos_phi_p)*(1-cos_phi_x));
00105  float seq = sqrt((1-cos_theta_p)*(1+cos_theta_x)) - sqrt((1+cos_theta_p)*(1-cos_theta_x));
00106  
00107 // cout<<"First factor: "<<feq/2<<endl;
00108 // cout<<"Second factor: "<<seq/2<<endl;
00109  
00110  vl(1) = feq/2;
00111  vl(2) = seq/2;
00112 
00113 // cout<<"Value "<<vl<<endl;
00114 //half angle corrected
00115 // vl(1) = (sin_x/(1+cos_x)) - (sin_p/(1+cos_p));
00116 // vl(2) = (sin_xt/(1+cos_xt)) - (sin_pt/(1+cos_pt));
00117  return std::pair<AlgebraicVector,AlgebraicVector>(vl,point);
00118 }
00119 
00120 std::pair<AlgebraicMatrix, AlgebraicVector> SimplePointingConstraint::makeDerivative(const AlgebraicVector& exPoint) const
00121 { 
00122  AlgebraicMatrix dr(2,7,0);
00123  AlgebraicVector point = exPoint;
00124  double dx = point(1) - refPoint.x();
00125  double dy = point(2) - refPoint.y();
00126  double dz = point(3) - refPoint.z();
00127  double px = point(4);
00128  double py = point(5); 
00129  double pz = point(6);
00130  
00131 
00132 //half angle solution
00133 //d/dx_i
00134  dr(1,1) = (sqrt((1 + dx/sqrt(pow(dx,2) + pow(dy,2)))*(1 - px/sqrt(pow(px,2) + pow(py,2)))) - 
00135             sqrt((1 - dx/sqrt(pow(dx,2) + pow(dy,2)))*(1 + px/sqrt(pow(px,2) + pow(py,2)))))/2.;
00136                                                                   
00137  dr(1,2) = (((-(pow(dx,2)/pow(pow(dx,2) + pow(dy,2),1.5)) + 1/sqrt(pow(dx,2) + pow(dy,2)))*
00138            (1 - px/sqrt(pow(px,2) + pow(py,2))))/
00139            (2.*sqrt((1 + dx/sqrt(pow(dx,2) + pow(dy,2)))*(1 - px/sqrt(pow(px,2) + pow(py,2))))) - 
00140            ((pow(dx,2)/pow(pow(dx,2) + pow(dy,2),1.5) - 1/sqrt(pow(dx,2) + pow(dy,2)))*
00141            (1 + px/sqrt(pow(px,2) + pow(py,2))))/
00142            (2.*sqrt((1 - dx/sqrt(pow(dx,2) + pow(dy,2)))*(1 + px/sqrt(pow(px,2) + pow(py,2))))))/2.;
00143  
00144  
00145  dr(1,3) = 0;
00146  
00147 //d/dp_i  
00148 //debug: x->p index xhange in denominator
00149  dr(1,4) = (-(dx*dy*(1 - px/sqrt(pow(px,2) + pow(py,2))))/
00150            (2.*pow(pow(dx,2) + pow(dy,2),1.5)*
00151            sqrt((1 + dx/sqrt(pow(dx,2) + pow(dy,2)))*(1 - px/sqrt(pow(px,2) + pow(py,2))))) - 
00152            (dx*dy*(1 + px/sqrt(pow(px,2) + pow(py,2))))/
00153            (2.*pow(pow(dx,2) + pow(dy,2),1.5)*
00154            sqrt((1 - dx/sqrt(pow(dx,2) + pow(dy,2)))*(1 + px/sqrt(pow(px,2) + pow(py,2))))))/2.;
00155         
00156                            
00157  dr(1,5) = (((1 + dx/sqrt(pow(dx,2) + pow(dy,2)))*px*py)/
00158            (2.*pow(pow(px,2) + pow(py,2),1.5)*
00159            sqrt((1 + dx/sqrt(pow(dx,2) + pow(dy,2)))*(1 - px/sqrt(pow(px,2) + pow(py,2))))) + 
00160            ((1 - dx/sqrt(pow(dx,2) + pow(dy,2)))*px*py)/
00161            (2.*pow(pow(px,2) + pow(py,2),1.5)*
00162            sqrt((1 - dx/sqrt(pow(dx,2) + pow(dy,2)))*(1 + px/sqrt(pow(px,2) + pow(py,2))))))/2.;
00163  
00164  
00165  
00166  dr(1,6) = 0;
00167  dr(1,7) = 0; 
00168 
00169 //2nd equation
00170 //d/dx_i
00171 
00172  dr(2,1) =(((-((dx*sqrt(pow(dx,2) + pow(dy,2)))/pow(pow(dx,2) + pow(dy,2) + pow(dz,2),1.5)) + 
00173           dx/(sqrt(pow(dx,2) + pow(dy,2))*sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2))))*
00174           (1 - sqrt(pow(px,2) + pow(py,2))/sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))/
00175           (2.*sqrt((1 + sqrt(pow(dx,2) + pow(dy,2))/sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2)))*
00176           (1 - sqrt(pow(px,2) + pow(py,2))/sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))) - 
00177           (((dx*sqrt(pow(dx,2) + pow(dy,2)))/pow(pow(dx,2) + pow(dy,2) + pow(dz,2),1.5) - 
00178           dx/(sqrt(pow(dx,2) + pow(dy,2))*sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2))))*
00179           (1 + sqrt(pow(px,2) + pow(py,2))/sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))/
00180           (2.*sqrt((1 - sqrt(pow(dx,2) + pow(dy,2))/sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2)))*
00181           (1 + sqrt(pow(px,2) + pow(py,2))/sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))))/2.;
00182            
00183  
00184  dr(2,2) = (((-((dy*sqrt(pow(dx,2) + pow(dy,2)))/pow(pow(dx,2) + pow(dy,2) + pow(dz,2),1.5)) + 
00185            dy/(sqrt(pow(dx,2) + pow(dy,2))*sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2))))*
00186            (1 - sqrt(pow(px,2) + pow(py,2))/sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))/
00187            (2.*sqrt((1 + sqrt(pow(dx,2) + pow(dy,2))/sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2)))*
00188            (1 - sqrt(pow(px,2) + pow(py,2))/sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))) - 
00189            (((dy*sqrt(pow(dx,2) + pow(dy,2)))/pow(pow(dx,2) + pow(dy,2) + pow(dz,2),1.5) - 
00190            dy/(sqrt(pow(dx,2) + pow(dy,2))*sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2))))*
00191            (1 + sqrt(pow(px,2) + pow(py,2))/sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))/
00192            (2.*sqrt((1 - sqrt(pow(dx,2) + pow(dy,2))/sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2)))*
00193            (1 + sqrt(pow(px,2) + pow(py,2))/sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))))/2.;
00194            
00195         
00196  dr(2,3) = (-(sqrt(pow(dx,2) + pow(dy,2))*dz*(1 - sqrt(pow(px,2) + pow(py,2))/
00197            sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))/
00198            (2.*pow(pow(dx,2) + pow(dy,2) + pow(dz,2),1.5)*
00199            sqrt((1 + sqrt(pow(dx,2) + pow(dy,2))/sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2)))*
00200            (1 - sqrt(pow(px,2) + pow(py,2))/sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))) - 
00201            (sqrt(pow(dx,2) + pow(dy,2))*dz*(1 + 
00202            sqrt(pow(px,2) + pow(py,2))/sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))/
00203            (2.*pow(pow(dx,2) + pow(dy,2) + pow(dz,2),1.5)*
00204            sqrt((1 - sqrt(pow(dx,2) + pow(dy,2))/sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2)))*
00205            (1 + sqrt(pow(px,2) + pow(py,2))/sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))))/2.;
00206            
00207 
00208  
00209 //d/dp_i 
00210 //debug: x->p index xhange in denominator
00211 
00212  dr(2,4) = (((1 + sqrt(pow(dx,2) + pow(dy,2))/sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2)))*
00213            ((px*sqrt(pow(px,2) + pow(py,2)))/pow(pow(px,2) + pow(py,2) + pow(pz,2),1.5) - 
00214            px/(sqrt(pow(px,2) + pow(py,2))*sqrt(pow(px,2) + pow(py,2) + pow(pz,2)))))/
00215            (2.*sqrt((1 + sqrt(pow(dx,2) + pow(dy,2))/sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2)))*
00216            (1 - sqrt(pow(px,2) + pow(py,2))/sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))) - 
00217            ((1 - sqrt(pow(dx,2) + pow(dy,2))/sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2)))*
00218            (-((px*sqrt(pow(px,2) + pow(py,2)))/pow(pow(px,2) + pow(py,2) + pow(pz,2),1.5)) + 
00219            px/(sqrt(pow(px,2) + pow(py,2))*sqrt(pow(px,2) + pow(py,2) + pow(pz,2)))))/
00220            (2.*sqrt((1 - sqrt(pow(dx,2) + pow(dy,2))/sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2)))*
00221            (1 + sqrt(pow(px,2) + pow(py,2))/sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))))/2.;
00222            
00223  
00224  dr(2,5) = (((1 + sqrt(pow(dx,2) + pow(dy,2))/sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2)))*
00225            ((py*sqrt(pow(px,2) + pow(py,2)))/pow(pow(px,2) + pow(py,2) + pow(pz,2),1.5) - 
00226            py/(sqrt(pow(px,2) + pow(py,2))*sqrt(pow(px,2) + pow(py,2) + pow(pz,2)))))/
00227            (2.*sqrt((1 + sqrt(pow(dx,2) + pow(dy,2))/sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2)))*
00228            (1 - sqrt(pow(px,2) + pow(py,2))/sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))) - 
00229            ((1 - sqrt(pow(dx,2) + pow(dy,2))/sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2)))*
00230            (-((py*sqrt(pow(px,2) + pow(py,2)))/pow(pow(px,2) + pow(py,2) + pow(pz,2),1.5)) + 
00231            py/(sqrt(pow(px,2) + pow(py,2))*sqrt(pow(px,2) + pow(py,2) + pow(pz,2)))))/
00232            (2.*sqrt((1 - sqrt(pow(dx,2) + pow(dy,2))/sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2)))*
00233            (1 + sqrt(pow(px,2) + pow(py,2))/sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))))/2.;
00234            
00235  
00236  dr(2,6) = (((1 + sqrt(pow(dx,2) + pow(dy,2))/sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2)))*
00237            sqrt(pow(px,2) + pow(py,2))*pz)/
00238            (2.*pow(pow(px,2) + pow(py,2) + pow(pz,2),1.5)*
00239            sqrt((1 + sqrt(pow(dx,2) + pow(dy,2))/sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2)))*
00240            (1 - sqrt(pow(px,2) + pow(py,2))/sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))) + 
00241            ((1 - sqrt(pow(dx,2) + pow(dy,2))/sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2)))*
00242             sqrt(pow(px,2) + pow(py,2))*pz)/
00243            (2.*pow(pow(px,2) + pow(py,2) + pow(pz,2),1.5)*
00244            sqrt((1 - sqrt(pow(dx,2) + pow(dy,2))/sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2)))*
00245            (1 + sqrt(pow(px,2) + pow(py,2))/sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))))/2.;
00246            
00247  
00248  dr(2,7) = 0;
00249  
00250 // cout<<"derivative matrix "<<dr<<endl;
00251  return std::pair<AlgebraicMatrix,AlgebraicVector>(dr,point); 
00252 }
/data/refman/pasoursint/CMSSW_5_3_3/src/RecoVertex/KinematicFit/src/SimplePointingConstraint.cc
