CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_2_7_hltpatch1/src/RecoVertex/KinematicFit/interface/ColinearityKinematicConstraintT.h

Go to the documentation of this file.
00001 #ifndef ColinearityKinematicConstraintT_H
00002 #define ColinearityKinematicConstraintT_H
00003 
00004 #include "RecoVertex/KinematicFitPrimitives/interface/MultiTrackKinematicConstraintT.h"
00005 #include "RecoVertex/KinematicFitPrimitives/interface/KinematicState.h"
00006 #include "DataFormats/CLHEP/interface/AlgebraicObjects.h"
00007 
00008 #include "RecoVertex/VertexPrimitives/interface/VertexException.h"
00009 
00010 
00019 namespace  colinearityKinematic {
00020   enum ConstraintDim {Phi=1, PhiTheta=2};
00021 }
00022 
00023 template<enum colinearityKinematic::ConstraintDim Dim>                                 
00024 class ColinearityKinematicConstraintT : public MultiTrackKinematicConstraintT<2,int(Dim)>{
00025 
00026 private:
00027   double a_1;
00028   double a_2;
00029 
00030   AlgebraicVector7 p1;
00031   AlgebraicVector7 p2;
00032 
00033   GlobalPoint point;
00034 
00035 public:
00036 
00037   typedef MultiTrackKinematicConstraintT<2,int(Dim)> super;
00038 
00039   ColinearityKinematicConstraintT(){}
00040   
00041   
00042   // initialize the constraint so it can precompute common qualtities to the three next call
00043   void init(const std::vector<KinematicState>& states,
00044             const GlobalPoint& ipoint,  const GlobalVector& fieldValue) {
00045     if(states.size()!=2) throw VertexException("ColinearityKinematicConstraint::<2 states passed");
00046 
00047     point = ipoint;
00048 
00049     a_1 = -states[0].particleCharge()*fieldValue.z();
00050     a_2 = -states[1].particleCharge()*fieldValue.z();
00051 
00052     p1 = states[0].kinematicParameters().vector();
00053     p2 = states[1].kinematicParameters().vector();
00054   }
00055 
00059   virtual int numberOfEquations() const {return Dim == colinearityKinematic::Phi ? 1 :2;}
00060   
00061   virtual ColinearityKinematicConstraintT<Dim> * clone()const
00062   {return new ColinearityKinematicConstraintT<Dim>(*this);}
00063 
00064 
00065 private:
00071   virtual void fillValue() const;
00072   
00073   
00079   virtual void fillParametersDerivative() const;
00080   
00086   virtual void fillPositionDerivative() const;
00087   
00088   
00089 };
00090 
00091 
00092 template<enum colinearityKinematic::ConstraintDim Dim>                                 
00093 void ColinearityKinematicConstraintT<Dim>::fillValue() const {
00094 
00095   typename super::valueType & vl = super::vl();
00096  
00097   double p1vx = p1(3) - a_1*(point.y() - p1(1));
00098   double p1vy = p1(4) + a_1*(point.x() - p1(0));
00099  
00100   double p2vx = p2(3) - a_2*(point.y() - p2(1));
00101   double p2vy = p2(4) + a_2*(point.x() - p2(0));
00102   
00103 
00104   // H_phi:
00105   vl(0)  = atan2(p1vy,p1vx) - atan2(p2vy,p2vx);
00106   if ( vl(0) >  M_PI ) vl(0) -= 2.0*M_PI;
00107   if ( vl(0) <= -M_PI ) vl(0) += 2.0*M_PI;
00108   // H_theta:
00109   if (Dim==colinearityKinematic::PhiTheta) {  
00110     double pt1  = sqrt(p1(3)*p1(3)+p1(4)*p1(4));
00111     double pt2  = sqrt(p2(3)*p2(3)+p2(4)*p2(4));
00112     vl(1)  = atan2(pt1,p1(5)) - atan2(pt2,p2(5));
00113     if ( vl(1) >  M_PI ) vl(1) -= 2.0*M_PI;
00114     if ( vl(1) <= -M_PI ) vl(1) += 2.0*M_PI;
00115   }
00116 }
00117 
00118 template<enum colinearityKinematic::ConstraintDim Dim>                                 
00119 void ColinearityKinematicConstraintT<Dim>::fillParametersDerivative() const {
00120 
00121 
00122   typename super::parametersDerivativeType & jac_d = super::jac_d(); 
00123 
00124   double p1vx = p1(3) - a_1*(point.y() - p1(1));
00125   double p1vy = p1(4) + a_1*(point.x() - p1(0));
00126   double k1 = 1.0/(p1vx*p1vx + p1vy*p1vy);
00127 
00128   double p2vx = p2(3) - a_2*(point.y() - p2(1));
00129   double p2vy = p2(4) + a_2*(point.x() - p2(0));
00130   double k2 = 1.0/(p2vx*p2vx + p2vy*p2vy);
00131 
00132   // H_phi:
00133 
00134   //x1 and x2 derivatives: 1st and 8th elements
00135   jac_d(0,0) =  -k1*a_1*p1vx;
00136   jac_d(0,7) =   k2*a_2*p2vx;
00137 
00138   //y1 and y2 derivatives: 2nd and 9th elements:
00139   jac_d(0,1) = -k1*a_1*p1vy;
00140   jac_d(0,8) =  k2*a_2*p2vy;
00141 
00142   //z1 and z2 components: 3d and 10th elmnets stay 0:
00143   jac_d(0,2)  = 0.; jac_d(0,9) = 0.;
00144 
00145   //px1 and px2 components: 4th and 11th elements:
00146   jac_d(0,3)  = -k1*p1vy;
00147   jac_d(0,10) =  k2*p2vy;
00148 
00149   //py1 and py2 components: 5th and 12 elements:
00150   jac_d(0,4)  =  k1*p1vx;
00151   jac_d(0,11) = -k2*p2vx;
00152 
00153 
00154   //pz1 and pz2 components: 6th and 13 elements:
00155   jac_d(0,5)  = 0.; jac_d(0,12) = 0.;
00156   //mass components: 7th and 14th elements:
00157   jac_d(0,6)  = 0.; jac_d(0,13) = 0.;
00158 
00159   if (Dim==colinearityKinematic::PhiTheta)  {
00160     double pt1 = sqrt(p1(3)*p1(3)+p1(4)*p1(4));
00161     double pTot1  = p1(3)*p1(3)+p1(4)*p1(4)+p1(5)*p1(5);
00162     double pt2  = sqrt(p2(3)*p2(3)+p2(4)*p2(4));
00163     double pTot2 = p2(3)*p2(3)+p2(4)*p2(4)+p2(5)*p2(5);
00164     
00165     // H_theta:
00166     //x1 and x2 derivatives: 1st and 8th elements
00167     jac_d(1,0) =  0.; jac_d(1,7) = 0.;
00168 
00169     //y1 and y2 derivatives: 2nd and 9th elements:
00170     jac_d(1,1) = 0.; jac_d(1,8) = 0.;
00171 
00172     //z1 and z2 components: 3d and 10th elmnets stay 0:
00173     jac_d(1,2) = 0.; jac_d(1,9) = 0.;
00174 
00175     jac_d(1,3)  =  p1(3) * (p1(5)/(pTot1*pt1));
00176     jac_d(1,10) =  p2(3) * (-p2(5)/(pTot2*pt2));
00177 
00178     //py1 and py2 components: 5th and 12 elements:
00179     jac_d(1,4)  =  p1(4)  * (p1(5)/(pTot1*pt1));
00180     jac_d(1,11) =  p2(4)  * (-p2(5)/(pTot2*pt2));
00181 
00182     //pz1 and pz2 components: 6th and 13 elements:
00183     jac_d(1,5)  = - pt1/pTot1;
00184     jac_d(1,12) =   pt2/pTot2;
00185 
00186     //mass components: 7th and 14th elements:
00187     jac_d(1,6)  = 0.; jac_d(1,13) = 0.;
00188   }
00189 }
00190 
00191 template<enum colinearityKinematic::ConstraintDim Dim>                                 
00192 void ColinearityKinematicConstraintT<Dim>::fillPositionDerivative() const
00193 {
00194 
00195   typename super::positionDerivativeType & jac_e = super::jac_e(); 
00196 
00197   double p1vx = p1(3) - a_1*(point.y() - p1(1));
00198   double p1vy = p1(4) + a_1*(point.x() - p1(0));
00199   double k1 = 1.0/(p1vx*p1vx + p1vy*p1vy);
00200 
00201   double p2vx = p2(3) - a_2*(point.y() - p2(1));
00202   double p2vy = p2(4) + a_2*(point.x() - p2(0));
00203   double k2 = 1.0/(p2vx*p2vx + p2vy*p2vy);
00204  
00205   // H_phi:
00206 
00207   // xv component
00208   jac_e(0,0) = k1*a_1*p1vx - k2*a_2*p2vx;
00209 
00210   //yv component
00211   jac_e(0,1) =  k1*a_1*p1vy - k2*a_2*p2vy;
00212 
00213   //zv component
00214   jac_e(0,2) = 0.;
00215 
00216   // H_theta: no correlation with vertex position
00217   if (Dim==colinearityKinematic::PhiTheta) {
00218     jac_e(1,0) = 0.;
00219     jac_e(1,1) = 0.;
00220     jac_e(1,2) = 0.;
00221   }
00222 }
00223 
00224 
00225 
00226 #endif