Reference Manual

Go to the documentation of this file.
00001 
00002 #include "Alignment/TwoBodyDecay/interface/TwoBodyDecayDerivatives.h"
00003 #include "Alignment/TwoBodyDecay/interface/TwoBodyDecayModel.h"
00004 
00005 #include "FWCore/Utilities/interface/Exception.h"
00006 
00007 #include <algorithm>
00008 
00009 TwoBodyDecayDerivatives::TwoBodyDecayDerivatives( double mPrimary, double mSecondary ) : 
00010   thePrimaryMass( mPrimary ), theSecondaryMass( mSecondary ) {}
00011 
00012 
00013 TwoBodyDecayDerivatives::~TwoBodyDecayDerivatives() {}
00014 
00015 
00016 const std::pair< AlgebraicMatrix, AlgebraicMatrix >
00017 TwoBodyDecayDerivatives::derivatives( const TwoBodyDecay & tbd )  const
00018 {
00019   return derivatives( tbd.decayParameters() );
00020 }
00021 
00022 const std::pair< AlgebraicMatrix, AlgebraicMatrix >
00023 TwoBodyDecayDerivatives::derivatives( const TwoBodyDecayParameters & param )  const
00024 {
00025   // get the derivatives with respect to all parameters
00026   std::pair< AlgebraicMatrix, AlgebraicMatrix > dqsdpx = this->dqsdpx( param );
00027   std::pair< AlgebraicMatrix, AlgebraicMatrix > dqsdpy = this->dqsdpy( param );
00028   std::pair< AlgebraicMatrix, AlgebraicMatrix > dqsdpz = this->dqsdpz( param );
00029   std::pair< AlgebraicMatrix, AlgebraicMatrix > dqsdtheta = this->dqsdtheta( param );
00030   std::pair< AlgebraicMatrix, AlgebraicMatrix > dqsdphi = this->dqsdphi( param );
00031   std::pair< AlgebraicMatrix, AlgebraicMatrix > dqsdm = this->dqsdm( param );
00032 
00033   AlgebraicMatrix dqplusdz( 3, dimension );
00034   dqplusdz.sub( 1, px, dqsdpx.first );
00035   dqplusdz.sub( 1, py, dqsdpy.first );
00036   dqplusdz.sub( 1, pz, dqsdpz.first );
00037   dqplusdz.sub( 1, theta, dqsdtheta.first );
00038   dqplusdz.sub( 1, phi, dqsdphi.first );
00039   dqplusdz.sub( 1, mass, dqsdm.first );
00040 
00041   AlgebraicMatrix dqminusdz( 3, dimension );
00042   dqminusdz.sub( 1, px, dqsdpx.second );
00043   dqminusdz.sub( 1, py, dqsdpy.second );
00044   dqminusdz.sub( 1, pz, dqsdpz.second );
00045   dqminusdz.sub( 1, theta, dqsdtheta.second );
00046   dqminusdz.sub( 1, phi, dqsdphi.second );
00047   dqminusdz.sub( 1, mass, dqsdm.second );
00048 
00049   return std::make_pair( dqplusdz, dqminusdz );
00050 }
00051 
00052 
00053 const std::pair< AlgebraicMatrix, AlgebraicMatrix >
00054 TwoBodyDecayDerivatives::selectedDerivatives( const TwoBodyDecay & tbd, const std::vector< bool > & selector ) const
00055 {
00056   return selectedDerivatives( tbd.decayParameters(), selector );
00057 }
00058 
00059 
00060 const std::pair< AlgebraicMatrix, AlgebraicMatrix >
00061 TwoBodyDecayDerivatives::selectedDerivatives( const TwoBodyDecayParameters & param, const std::vector< bool > & selector ) const
00062 {
00063   if ( selector.size() != dimension )
00064   {
00065     throw cms::Exception( "BadConfig" ) << "@SUB=TwoBodyDecayDerivatives::selectedDerivatives"
00066                                         << "selector has bad dimension (size=" << selector.size() << ").";
00067   }
00068 
00069   int nSelected = std::count( selector.begin(), selector.end(), true );
00070   int iSelected = 1;
00071 
00072   AlgebraicMatrix dqplusdz( 3, nSelected );
00073   AlgebraicMatrix dqminusdz( 3, nSelected );
00074   std::pair< AlgebraicMatrix, AlgebraicMatrix > dqsdzi;
00075 
00076   for ( unsigned int i = 1; i <= dimension; i++ )
00077   {
00078     if ( selector[i] )
00079     {
00080       dqsdzi = this->dqsdzi( param, DerivativeParameterName(i) );
00081       dqplusdz.sub( 1, iSelected, dqsdzi.first );
00082       dqminusdz.sub( 1, iSelected, dqsdzi.second );
00083       iSelected++;
00084     }
00085   }
00086 
00087   return std::make_pair( dqplusdz, dqminusdz );
00088 }
00089 
00090 
00091 const std::pair< AlgebraicMatrix, AlgebraicMatrix > TwoBodyDecayDerivatives::dqsdpx( const TwoBodyDecayParameters & param ) const
00092 {
00093   double px = param[TwoBodyDecayParameters::px];
00094   double py = param[TwoBodyDecayParameters::py];
00095   double pz = param[TwoBodyDecayParameters::pz];
00096   double theta = param[TwoBodyDecayParameters::theta];
00097   double phi = param[TwoBodyDecayParameters::phi];
00098 
00099   // compute transverse and absolute momentum
00100   double pT2 = px*px + py*py;
00101   double p2 = pT2 + pz*pz;
00102   double pT = sqrt( pT2 );
00103   double p = sqrt( p2 );
00104 
00105   double sphi = sin( phi );
00106   double cphi = cos( phi );
00107   double stheta = sin( theta );
00108   double ctheta = cos( theta );
00109 
00110   // some constants from kinematics
00111   double c1 = 0.5*thePrimaryMass/theSecondaryMass;
00112   double c2 = sqrt( c1*c1 - 1. );
00113   double c3 = 0.5*c2*ctheta/c1;
00114   double c4 = sqrt( p2 + thePrimaryMass*thePrimaryMass );
00115 
00116   // momentum of decay particle 1 in the primary's boosted frame
00117   AlgebraicMatrix pplus( 3, 1 );
00118   pplus[0][0] = theSecondaryMass*c2*stheta*cphi;
00119   pplus[1][0] = theSecondaryMass*c2*stheta*sphi;
00120   pplus[2][0] = 0.5*p + c3*c4;
00121 
00122   // momentum of decay particle 2 in the primary's boosted frame
00123   AlgebraicMatrix pminus( 3, 1 );
00124   pminus[0][0] = -pplus[0][0];
00125   pminus[1][0] = -pplus[1][0];
00126   pminus[2][0] = 0.5*p - c3*c4;
00127 
00128   // derivative of rotation matrix w.r.t. px
00129   AlgebraicMatrix dRotMatdpx( 3, 3 );
00130 
00131   dRotMatdpx[0][0] = pz/(pT*p)*(1.-px*px*(1./pT2+1./p2));
00132   dRotMatdpx[0][1] = px*py/(pT*pT2);
00133   dRotMatdpx[0][2] = (1.-px*px/p2)/p;
00134 
00135   dRotMatdpx[1][0] = -px*py*pz/(pT*p)*(1./pT2+1./p2);
00136   dRotMatdpx[1][1] = (1.-px*px/pT2)/pT;
00137   dRotMatdpx[1][2] = -px*py/(p*p2);
00138 
00139   dRotMatdpx[2][0] = -(1./pT-pT/p2)*px/p;
00140   dRotMatdpx[2][1] = 0.;
00141   dRotMatdpx[2][2] = -px*pz/(p*p2);
00142 
00143   // derivative of the momentum of particle 1 in the lab frame w.r.t. px
00144   double dpplusdpx = px*( 0.5/p + c3/c4 );
00145 
00146   AlgebraicMatrix dqplusdpx = dRotMatdpx*pplus;
00147   dqplusdpx[0][0] += px*dpplusdpx/p;
00148   dqplusdpx[1][0] += py*dpplusdpx/p;
00149   dqplusdpx[2][0] += pz*dpplusdpx/p;
00150 
00151   // derivative of the momentum of particle 2 in the lab frame w.r.t. px
00152   double dpminusdpx = px*( 0.5/p - c3/c4 );
00153 
00154   AlgebraicMatrix dqminusdpx = dRotMatdpx*pminus; 
00155   dqminusdpx[0][0] += px*dpminusdpx/p;
00156   dqminusdpx[1][0] += py*dpminusdpx/p;
00157   dqminusdpx[2][0] += pz*dpminusdpx/p;
00158 
00159   // return result
00160   return std::make_pair( dqplusdpx, dqminusdpx );
00161 }
00162 
00163 
00164 const std::pair< AlgebraicMatrix, AlgebraicMatrix > TwoBodyDecayDerivatives::dqsdpy( const TwoBodyDecayParameters & param ) const
00165 {
00166   double px = param[TwoBodyDecayParameters::px];
00167   double py = param[TwoBodyDecayParameters::py];
00168   double pz = param[TwoBodyDecayParameters::pz];
00169   double theta = param[TwoBodyDecayParameters::theta];
00170   double phi = param[TwoBodyDecayParameters::phi];
00171 
00172   // compute transverse and absolute momentum
00173   double pT2 = px*px + py*py;
00174   double p2 = pT2 + pz*pz;
00175   double pT = sqrt( pT2 );
00176   double p = sqrt( p2 );
00177 
00178   double sphi = sin( phi );
00179   double cphi = cos( phi );
00180   double stheta = sin( theta );
00181   double ctheta = cos( theta );
00182 
00183   // some constants from kinematics
00184   double c1 = 0.5*thePrimaryMass/theSecondaryMass;
00185   double c2 = sqrt( c1*c1 - 1. );
00186   double c3 = 0.5*c2*ctheta/c1;
00187   double c4 = sqrt( p2 + thePrimaryMass*thePrimaryMass );
00188 
00189   // momentum of decay particle 1 in the rest frame of the primary
00190   AlgebraicMatrix pplus( 3, 1 );
00191   pplus[0][0] = theSecondaryMass*c2*stheta*cphi;
00192   pplus[1][0] = theSecondaryMass*c2*stheta*sphi;
00193   pplus[2][0] = 0.5*p + c3*c4;
00194 
00195   // momentum of decay particle 2 in the rest frame of the primary
00196   AlgebraicMatrix pminus( 3, 1 );
00197   pminus[0][0] = -pplus[0][0];
00198   pminus[1][0] = -pplus[1][0];
00199   pminus[2][0] = 0.5*p - c3*c4;
00200 
00201   // derivative of rotation matrix w.r.t. py
00202   AlgebraicMatrix dRotMatdpy( 3, 3 );
00203 
00204   dRotMatdpy[0][0] = -px*py*pz/(pT*p)*(1./pT2+1./p2);
00205   dRotMatdpy[0][1] = (py*py/pT2-1.)/pT;
00206   dRotMatdpy[0][2] = -px*py/(p*p2);
00207 
00208   dRotMatdpy[1][0] = pz/(pT*p)*(1.-py*py*(1./pT2+1./p2));
00209   dRotMatdpy[1][1] = -px*py/(pT*pT2);
00210   dRotMatdpy[1][2] = (1.-py*py/p2)/p;
00211 
00212   dRotMatdpy[2][0] = -(1./pT-pT/p2)*py/p;
00213   dRotMatdpy[2][1] = 0.;
00214   dRotMatdpy[2][2] = -py*pz/(p*p2);
00215 
00216   // derivative of the momentum of particle 1 in the lab frame w.r.t. py
00217   double dpplusdpy = py*( 0.5/p + c3/c4 );
00218 
00219   AlgebraicMatrix dqplusdpy = dRotMatdpy*pplus;
00220   dqplusdpy[0][0] += px*dpplusdpy/p;
00221   dqplusdpy[1][0] += py*dpplusdpy/p;
00222   dqplusdpy[2][0] += pz*dpplusdpy/p;
00223 
00224   // derivative of the momentum of particle 2 in the lab frame w.r.t. py
00225   double dpminusdpy = py*( 0.5/p - c3/c4 );
00226 
00227   AlgebraicMatrix dqminusdpy = dRotMatdpy*pminus;
00228   dqminusdpy[0][0] += px*dpminusdpy/p;
00229   dqminusdpy[1][0] += py*dpminusdpy/p;
00230   dqminusdpy[2][0] += pz*dpminusdpy/p;
00231 
00232   // return result
00233   return std::make_pair( dqplusdpy, dqminusdpy );
00234 }
00235 
00236 
00237 const std::pair< AlgebraicMatrix, AlgebraicMatrix > TwoBodyDecayDerivatives::dqsdpz( const TwoBodyDecayParameters & param ) const
00238 {
00239   double px = param[TwoBodyDecayParameters::px];
00240   double py = param[TwoBodyDecayParameters::py];
00241   double pz = param[TwoBodyDecayParameters::pz];
00242   double theta = param[TwoBodyDecayParameters::theta];
00243   double phi = param[TwoBodyDecayParameters::phi];
00244 
00245   // compute transverse and absolute momentum
00246   double pT2 = px*px + py*py;
00247   double p2 = pT2 + pz*pz;
00248   double pT = sqrt( pT2 );
00249   double p = sqrt( p2 );
00250 
00251   double sphi = sin( phi );
00252   double cphi = cos( phi );
00253   double stheta = sin( theta );
00254   double ctheta = cos( theta );
00255 
00256   // some constants from kinematics
00257   double c1 = 0.5*thePrimaryMass/theSecondaryMass;
00258   double c2 = sqrt( c1*c1 - 1. );
00259   double c3 = 0.5*c2*ctheta/c1;
00260   double c4 = sqrt( p2 + thePrimaryMass*thePrimaryMass );
00261 
00262   // momentum of decay particle 1 in the rest frame of the primary
00263   AlgebraicMatrix pplus( 3, 1 );
00264   pplus[0][0] = theSecondaryMass*c2*stheta*cphi;
00265   pplus[1][0] = theSecondaryMass*c2*stheta*sphi;
00266   pplus[2][0] = 0.5*p + c3*c4;
00267 
00268   // momentum of decay particle 2 in the rest frame of the primary
00269   AlgebraicMatrix pminus( 3, 1 );
00270   pminus[0][0] = -pplus[0][0];
00271   pminus[1][0] = -pplus[1][0];
00272   pminus[2][0] = 0.5*p - c3*c4;
00273 
00274   // derivative of rotation matrix w.r.t. py
00275   AlgebraicMatrix dRotMatdpz( 3, 3 );
00276 
00277   dRotMatdpz[0][0] = px/(pT*p)*(1.-pz*pz/p2);
00278   dRotMatdpz[0][1] = 0.;
00279   dRotMatdpz[0][2] = -px*pz/(p*p2);
00280 
00281   dRotMatdpz[1][0] = py/(p*pT)*(1.-pz*pz/p2);
00282   dRotMatdpz[1][1] = 0.;
00283   dRotMatdpz[1][2] = -py*pz/(p*p2);
00284 
00285   dRotMatdpz[2][0] = pT*pz/(p*p2);
00286   dRotMatdpz[2][1] = 0.;
00287   dRotMatdpz[2][2] = (1.-pz*pz/p2)/p;
00288 
00289   // derivative of the momentum of particle 1 in the lab frame w.r.t. pz
00290   double dpplusdpz = pz*( 0.5/p + c3/c4 );
00291 
00292   AlgebraicMatrix dqplusdpz = dRotMatdpz*pplus;
00293   dqplusdpz[0][0] += px*dpplusdpz/p;
00294   dqplusdpz[1][0] += py*dpplusdpz/p;
00295   dqplusdpz[2][0] += pz*dpplusdpz/p;
00296 
00297   // derivative of the momentum of particle 2 in the lab frame w.r.t. pz
00298   double dpminusdpz = pz*( 0.5/p - c3/c4 );
00299 
00300   AlgebraicMatrix dqminusdpz = dRotMatdpz*pminus;
00301   dqminusdpz[0][0] += px*dpminusdpz/p;
00302   dqminusdpz[1][0] += py*dpminusdpz/p;
00303   dqminusdpz[2][0] += pz*dpminusdpz/p;
00304 
00305   // return result
00306   return std::make_pair( dqplusdpz, dqminusdpz );
00307 }
00308 
00309 
00310 const std::pair< AlgebraicMatrix, AlgebraicMatrix > TwoBodyDecayDerivatives::dqsdtheta( const TwoBodyDecayParameters & param ) const
00311 {
00312   double px = param[TwoBodyDecayParameters::px];
00313   double py = param[TwoBodyDecayParameters::py];
00314   double pz = param[TwoBodyDecayParameters::pz];
00315   double theta = param[TwoBodyDecayParameters::theta];
00316   double phi = param[TwoBodyDecayParameters::phi];
00317 
00318   // compute transverse and absolute momentum
00319   double pT2 = px*px + py*py;
00320   double p2 = pT2 + pz*pz;
00321 
00322   double sphi = sin( phi );
00323   double cphi = cos( phi );
00324   double stheta = sin( theta );
00325   double ctheta = cos( theta );
00326 
00327   // some constants from kinematics
00328   double c1 = 0.5*thePrimaryMass/theSecondaryMass;
00329   double c2 = sqrt( c1*c1 - 1. );
00330   double c3 = -0.5*c2*stheta/c1;
00331   double c4 = sqrt( p2 + thePrimaryMass*thePrimaryMass );
00332 
00333   // derivative of the momentum of particle 1 in the primary's rest frame w.r.t. theta
00334   AlgebraicMatrix dpplusdtheta( 3, 1 );
00335   dpplusdtheta[0][0] = theSecondaryMass*c2*ctheta*cphi;
00336   dpplusdtheta[1][0] = theSecondaryMass*c2*ctheta*sphi;
00337   dpplusdtheta[2][0] = c3*c4;
00338 
00339   // derivative of the momentum of particle 2 in the primary's rest frame w.r.t. theta
00340   AlgebraicMatrix dpminusdtheta( 3, 1 );
00341   dpminusdtheta[0][0] = -theSecondaryMass*c2*ctheta*cphi;
00342   dpminusdtheta[1][0] = -theSecondaryMass*c2*ctheta*sphi;
00343   dpminusdtheta[2][0] = -c3*c4;
00344 
00345   TwoBodyDecayModel decayModel;
00346   AlgebraicMatrix rotMat = decayModel.rotationMatrix( px, py, pz );
00347 
00348   AlgebraicMatrix dqplusdtheta = rotMat*dpplusdtheta;
00349   AlgebraicMatrix dqminusdtheta = rotMat*dpminusdtheta;
00350 
00351   return std::make_pair( dqplusdtheta, dqminusdtheta );
00352 }
00353 
00354 
00355 const std::pair< AlgebraicMatrix, AlgebraicMatrix > TwoBodyDecayDerivatives::dqsdphi( const TwoBodyDecayParameters & param ) const
00356 {
00357   double px = param[TwoBodyDecayParameters::px];
00358   double py = param[TwoBodyDecayParameters::py];
00359   double pz = param[TwoBodyDecayParameters::pz];
00360   double theta = param[TwoBodyDecayParameters::theta];
00361   double phi = param[TwoBodyDecayParameters::phi];
00362 
00363   double sphi = sin( phi );
00364   double cphi = cos( phi );
00365   double stheta = sin( theta );
00366 
00367   // some constants from kinematics
00368   double c1 = 0.5*thePrimaryMass/theSecondaryMass;
00369   double c2 = sqrt( c1*c1 - 1. );
00370 
00371   // derivative of the momentum of particle 1 in the primary's rest frame w.r.t. phi
00372   AlgebraicMatrix dpplusdphi( 3, 1 );
00373   dpplusdphi[0][0] = -theSecondaryMass*c2*stheta*sphi;
00374   dpplusdphi[1][0] = theSecondaryMass*c2*stheta*cphi;
00375   dpplusdphi[2][0] = 0.;
00376 
00377   // derivative of the momentum of particle 2 in the primary's rest frame w.r.t. phi
00378   AlgebraicMatrix dpminusdphi( 3, 1 );
00379   dpminusdphi[0][0] = theSecondaryMass*c2*stheta*sphi;
00380   dpminusdphi[1][0] = -theSecondaryMass*c2*stheta*cphi;
00381   dpminusdphi[2][0] = 0.;
00382 
00383   TwoBodyDecayModel decayModel;
00384   AlgebraicMatrix rotMat = decayModel.rotationMatrix( px, py, pz );
00385 
00386   AlgebraicMatrix dqplusdphi = rotMat*dpplusdphi;
00387   AlgebraicMatrix dqminusdphi = rotMat*dpminusdphi;
00388 
00389   return std::make_pair( dqplusdphi, dqminusdphi );
00390 }
00391 
00392 
00393 const std::pair< AlgebraicMatrix, AlgebraicMatrix > TwoBodyDecayDerivatives::dqsdm( const TwoBodyDecayParameters & param ) const
00394 {
00395   double px = param[TwoBodyDecayParameters::px];
00396   double py = param[TwoBodyDecayParameters::py];
00397   double pz = param[TwoBodyDecayParameters::pz];
00398   double theta = param[TwoBodyDecayParameters::theta];
00399   double phi = param[TwoBodyDecayParameters::phi];
00400 
00401   double pT2 = px*px + py*py;
00402   double p2 = pT2 + pz*pz;
00403 
00404   double sphi = sin( phi );
00405   double cphi = cos( phi );
00406   double ctheta = cos( theta );
00407   double stheta = sin( theta );
00408 
00409   // some constants from kinematics
00410   double c1 = 0.5*thePrimaryMass/theSecondaryMass;
00411   double c2 = 1./sqrt( c1*c1 - 1. );
00412   double m2 = thePrimaryMass*thePrimaryMass;
00413 
00414   // derivative of the momentum of particle 1 in the primary's rest frame w.r.t. the primary's mass
00415   AlgebraicMatrix dpplusdm( 3, 1 );
00416   dpplusdm[0][0] = c2*0.5*c1*stheta*cphi;
00417   dpplusdm[1][0] = c2*0.5*c1*stheta*sphi;
00418   dpplusdm[2][0] = c2*theSecondaryMass*( c1*c1 + p2/m2 )/sqrt( p2 + m2 )*ctheta;
00419 
00420   // derivative of the momentum of particle 2 in the primary's rest frame w.r.t. the primary's mass
00421   AlgebraicMatrix dpminusdm( 3, 1 );
00422   dpminusdm[0][0] = -dpplusdm[0][0];
00423   dpminusdm[1][0] = -dpplusdm[1][0];
00424   dpminusdm[2][0] = -dpplusdm[2][0];
00425 
00426   TwoBodyDecayModel decayModel;
00427   AlgebraicMatrix rotMat = decayModel.rotationMatrix( px, py, pz );
00428 
00429   AlgebraicMatrix dqplusdm = rotMat*dpplusdm;
00430   AlgebraicMatrix dqminusdm = rotMat*dpminusdm;
00431 
00432   return std::make_pair( dqplusdm, dqminusdm );
00433 }
00434 
00435 
00436 const std::pair< AlgebraicMatrix, AlgebraicMatrix >
00437 TwoBodyDecayDerivatives::dqsdzi( const TwoBodyDecayParameters & param, const DerivativeParameterName & i ) const
00438 {
00439   switch ( i )
00440   {
00441   case TwoBodyDecayDerivatives::px :
00442     return dqsdpx( param );
00443     break;
00444   case TwoBodyDecayDerivatives::py :
00445     return dqsdpy( param );
00446     break;
00447   case TwoBodyDecayDerivatives::pz :
00448     return dqsdpz( param );
00449     break;
00450   case TwoBodyDecayDerivatives::theta :
00451     return dqsdtheta( param );
00452     break;
00453   case TwoBodyDecayDerivatives::phi :
00454     return dqsdphi( param );
00455     break;
00456   case TwoBodyDecayDerivatives::mass :
00457     return dqsdm( param );
00458     break;
00459   default:
00460     throw cms::Exception( "BadConfig" ) << "@SUB=TwoBodyDecayDerivatives::dqsdzi"
00461                                         << "no decay parameter related to selector (" << i << ").";
00462   };
00463 
00464   return std::make_pair( AlgebraicMatrix( 3, 1, 0 ), AlgebraicMatrix( 3, 1, 0 ) );
00465 }
/data/refman/pasoursint/CMSSW_6_1_2_SLHC4_patch1/src/Alignment/TwoBodyDecay/src/TwoBodyDecayDerivatives.cc
