TwoBodyDecay/src/TwoBodyDecayDerivatives.cc

0001
0002 #include "Alignment/TwoBodyDecay/interface/TwoBodyDecayDerivatives.h"
0003 #include "Alignment/TwoBodyDecay/interface/TwoBodyDecayModel.h"
0004
0005 #include "FWCore/Utilities/interface/Exception.h"
0006
0007 #include <algorithm>
0008
0009 TwoBodyDecayDerivatives::TwoBodyDecayDerivatives(double mPrimary, double mSecondary)
0010     : thePrimaryMass(mPrimary), theSecondaryMass(mSecondary) {}
0011
0012 TwoBodyDecayDerivatives::~TwoBodyDecayDerivatives() {}
0013
0014 const std::pair<AlgebraicMatrix, AlgebraicMatrix> TwoBodyDecayDerivatives::derivatives(const TwoBodyDecay &tbd) const {
0015   return derivatives(tbd.decayParameters());
0016 }
0017
0018 const std::pair<AlgebraicMatrix, AlgebraicMatrix> TwoBodyDecayDerivatives::derivatives(
0019     const TwoBodyDecayParameters &param) const {
0020   // get the derivatives with respect to all parameters
0021   std::pair<AlgebraicMatrix, AlgebraicMatrix> dqsdpx = this->dqsdpx(param);
0022   std::pair<AlgebraicMatrix, AlgebraicMatrix> dqsdpy = this->dqsdpy(param);
0023   std::pair<AlgebraicMatrix, AlgebraicMatrix> dqsdpz = this->dqsdpz(param);
0024   std::pair<AlgebraicMatrix, AlgebraicMatrix> dqsdtheta = this->dqsdtheta(param);
0025   std::pair<AlgebraicMatrix, AlgebraicMatrix> dqsdphi = this->dqsdphi(param);
0026   std::pair<AlgebraicMatrix, AlgebraicMatrix> dqsdm = this->dqsdm(param);
0027
0028   AlgebraicMatrix dqplusdz(3, dimension);
0029   dqplusdz.sub(1, px, dqsdpx.first);
0030   dqplusdz.sub(1, py, dqsdpy.first);
0031   dqplusdz.sub(1, pz, dqsdpz.first);
0032   dqplusdz.sub(1, theta, dqsdtheta.first);
0033   dqplusdz.sub(1, phi, dqsdphi.first);
0034   dqplusdz.sub(1, mass, dqsdm.first);
0035
0036   AlgebraicMatrix dqminusdz(3, dimension);
0037   dqminusdz.sub(1, px, dqsdpx.second);
0038   dqminusdz.sub(1, py, dqsdpy.second);
0039   dqminusdz.sub(1, pz, dqsdpz.second);
0040   dqminusdz.sub(1, theta, dqsdtheta.second);
0041   dqminusdz.sub(1, phi, dqsdphi.second);
0042   dqminusdz.sub(1, mass, dqsdm.second);
0043
0044   return std::make_pair(dqplusdz, dqminusdz);
0045 }
0046
0047 const std::pair<AlgebraicMatrix, AlgebraicMatrix> TwoBodyDecayDerivatives::selectedDerivatives(
0048     const TwoBodyDecay &tbd, const std::vector<bool> &selector) const {
0049   return selectedDerivatives(tbd.decayParameters(), selector);
0050 }
0051
0052 const std::pair<AlgebraicMatrix, AlgebraicMatrix> TwoBodyDecayDerivatives::selectedDerivatives(
0053     const TwoBodyDecayParameters &param, const std::vector<bool> &selector) const {
0054   if (selector.size() != dimension) {
0055     throw cms::Exception("BadConfig") << "@SUB=TwoBodyDecayDerivatives::selectedDerivatives"
0056                                       << "selector has bad dimension (size=" << selector.size() << ").";
0057   }
0058
0059   int nSelected = std::count(selector.begin(), selector.end(), true);
0060   int iSelected = 1;
0061
0062   AlgebraicMatrix dqplusdz(3, nSelected);
0063   AlgebraicMatrix dqminusdz(3, nSelected);
0064   std::pair<AlgebraicMatrix, AlgebraicMatrix> dqsdzi;
0065
0066   for (unsigned int i = 1; i <= dimension; i++) {
0067     if (selector[i]) {
0068       dqsdzi = this->dqsdzi(param, DerivativeParameterName(i));
0069       dqplusdz.sub(1, iSelected, dqsdzi.first);
0070       dqminusdz.sub(1, iSelected, dqsdzi.second);
0071       iSelected++;
0072     }
0073   }
0074
0075   return std::make_pair(dqplusdz, dqminusdz);
0076 }
0077
0078 const std::pair<AlgebraicMatrix, AlgebraicMatrix> TwoBodyDecayDerivatives::dqsdpx(
0079     const TwoBodyDecayParameters &param) const {
0080   double px = param[TwoBodyDecayParameters::px];
0081   double py = param[TwoBodyDecayParameters::py];
0082   double pz = param[TwoBodyDecayParameters::pz];
0083   double theta = param[TwoBodyDecayParameters::theta];
0084   double phi = param[TwoBodyDecayParameters::phi];
0085
0086   // compute transverse and absolute momentum
0087   double pT2 = px * px + py * py;
0088   double p2 = pT2 + pz * pz;
0089   double pT = sqrt(pT2);
0090   double p = sqrt(p2);
0091
0092   double sphi = sin(phi);
0093   double cphi = cos(phi);
0094   double stheta = sin(theta);
0095   double ctheta = cos(theta);
0096
0097   // some constants from kinematics
0098   double c1 = 0.5 * thePrimaryMass / theSecondaryMass;
0099   double c2 = sqrt(c1 * c1 - 1.);
0100   double c3 = 0.5 * c2 * ctheta / c1;
0101   double c4 = sqrt(p2 + thePrimaryMass * thePrimaryMass);
0102
0103   // momentum of decay particle 1 in the primary's boosted frame
0104   AlgebraicMatrix pplus(3, 1);
0105   pplus[0][0] = theSecondaryMass * c2 * stheta * cphi;
0106   pplus[1][0] = theSecondaryMass * c2 * stheta * sphi;
0107   pplus[2][0] = 0.5 * p + c3 * c4;
0108
0109   // momentum of decay particle 2 in the primary's boosted frame
0110   AlgebraicMatrix pminus(3, 1);
0111   pminus[0][0] = -pplus[0][0];
0112   pminus[1][0] = -pplus[1][0];
0113   pminus[2][0] = 0.5 * p - c3 * c4;
0114
0115   // derivative of rotation matrix w.r.t. px
0116   AlgebraicMatrix dRotMatdpx(3, 3);
0117
0118   dRotMatdpx[0][0] = pz / (pT * p) * (1. - px * px * (1. / pT2 + 1. / p2));
0119   dRotMatdpx[0][1] = px * py / (pT * pT2);
0120   dRotMatdpx[0][2] = (1. - px * px / p2) / p;
0121
0122   dRotMatdpx[1][0] = -px * py * pz / (pT * p) * (1. / pT2 + 1. / p2);
0123   dRotMatdpx[1][1] = (1. - px * px / pT2) / pT;
0124   dRotMatdpx[1][2] = -px * py / (p * p2);
0125
0126   dRotMatdpx[2][0] = -(1. / pT - pT / p2) * px / p;
0127   dRotMatdpx[2][1] = 0.;
0128   dRotMatdpx[2][2] = -px * pz / (p * p2);
0129
0130   // derivative of the momentum of particle 1 in the lab frame w.r.t. px
0131   double dpplusdpx = px * (0.5 / p + c3 / c4);
0132
0133   AlgebraicMatrix dqplusdpx = dRotMatdpx * pplus;
0134   dqplusdpx[0][0] += px * dpplusdpx / p;
0135   dqplusdpx[1][0] += py * dpplusdpx / p;
0136   dqplusdpx[2][0] += pz * dpplusdpx / p;
0137
0138   // derivative of the momentum of particle 2 in the lab frame w.r.t. px
0139   double dpminusdpx = px * (0.5 / p - c3 / c4);
0140
0141   AlgebraicMatrix dqminusdpx = dRotMatdpx * pminus;
0142   dqminusdpx[0][0] += px * dpminusdpx / p;
0143   dqminusdpx[1][0] += py * dpminusdpx / p;
0144   dqminusdpx[2][0] += pz * dpminusdpx / p;
0145
0146   // return result
0147   return std::make_pair(dqplusdpx, dqminusdpx);
0148 }
0149
0150 const std::pair<AlgebraicMatrix, AlgebraicMatrix> TwoBodyDecayDerivatives::dqsdpy(
0151     const TwoBodyDecayParameters &param) const {
0152   double px = param[TwoBodyDecayParameters::px];
0153   double py = param[TwoBodyDecayParameters::py];
0154   double pz = param[TwoBodyDecayParameters::pz];
0155   double theta = param[TwoBodyDecayParameters::theta];
0156   double phi = param[TwoBodyDecayParameters::phi];
0157
0158   // compute transverse and absolute momentum
0159   double pT2 = px * px + py * py;
0160   double p2 = pT2 + pz * pz;
0161   double pT = sqrt(pT2);
0162   double p = sqrt(p2);
0163
0164   double sphi = sin(phi);
0165   double cphi = cos(phi);
0166   double stheta = sin(theta);
0167   double ctheta = cos(theta);
0168
0169   // some constants from kinematics
0170   double c1 = 0.5 * thePrimaryMass / theSecondaryMass;
0171   double c2 = sqrt(c1 * c1 - 1.);
0172   double c3 = 0.5 * c2 * ctheta / c1;
0173   double c4 = sqrt(p2 + thePrimaryMass * thePrimaryMass);
0174
0175   // momentum of decay particle 1 in the rest frame of the primary
0176   AlgebraicMatrix pplus(3, 1);
0177   pplus[0][0] = theSecondaryMass * c2 * stheta * cphi;
0178   pplus[1][0] = theSecondaryMass * c2 * stheta * sphi;
0179   pplus[2][0] = 0.5 * p + c3 * c4;
0180
0181   // momentum of decay particle 2 in the rest frame of the primary
0182   AlgebraicMatrix pminus(3, 1);
0183   pminus[0][0] = -pplus[0][0];
0184   pminus[1][0] = -pplus[1][0];
0185   pminus[2][0] = 0.5 * p - c3 * c4;
0186
0187   // derivative of rotation matrix w.r.t. py
0188   AlgebraicMatrix dRotMatdpy(3, 3);
0189
0190   dRotMatdpy[0][0] = -px * py * pz / (pT * p) * (1. / pT2 + 1. / p2);
0191   dRotMatdpy[0][1] = (py * py / pT2 - 1.) / pT;
0192   dRotMatdpy[0][2] = -px * py / (p * p2);
0193
0194   dRotMatdpy[1][0] = pz / (pT * p) * (1. - py * py * (1. / pT2 + 1. / p2));
0195   dRotMatdpy[1][1] = -px * py / (pT * pT2);
0196   dRotMatdpy[1][2] = (1. - py * py / p2) / p;
0197
0198   dRotMatdpy[2][0] = -(1. / pT - pT / p2) * py / p;
0199   dRotMatdpy[2][1] = 0.;
0200   dRotMatdpy[2][2] = -py * pz / (p * p2);
0201
0202   // derivative of the momentum of particle 1 in the lab frame w.r.t. py
0203   double dpplusdpy = py * (0.5 / p + c3 / c4);
0204
0205   AlgebraicMatrix dqplusdpy = dRotMatdpy * pplus;
0206   dqplusdpy[0][0] += px * dpplusdpy / p;
0207   dqplusdpy[1][0] += py * dpplusdpy / p;
0208   dqplusdpy[2][0] += pz * dpplusdpy / p;
0209
0210   // derivative of the momentum of particle 2 in the lab frame w.r.t. py
0211   double dpminusdpy = py * (0.5 / p - c3 / c4);
0212
0213   AlgebraicMatrix dqminusdpy = dRotMatdpy * pminus;
0214   dqminusdpy[0][0] += px * dpminusdpy / p;
0215   dqminusdpy[1][0] += py * dpminusdpy / p;
0216   dqminusdpy[2][0] += pz * dpminusdpy / p;
0217
0218   // return result
0219   return std::make_pair(dqplusdpy, dqminusdpy);
0220 }
0221
0222 const std::pair<AlgebraicMatrix, AlgebraicMatrix> TwoBodyDecayDerivatives::dqsdpz(
0223     const TwoBodyDecayParameters &param) const {
0224   double px = param[TwoBodyDecayParameters::px];
0225   double py = param[TwoBodyDecayParameters::py];
0226   double pz = param[TwoBodyDecayParameters::pz];
0227   double theta = param[TwoBodyDecayParameters::theta];
0228   double phi = param[TwoBodyDecayParameters::phi];
0229
0230   // compute transverse and absolute momentum
0231   double pT2 = px * px + py * py;
0232   double p2 = pT2 + pz * pz;
0233   double pT = sqrt(pT2);
0234   double p = sqrt(p2);
0235
0236   double sphi = sin(phi);
0237   double cphi = cos(phi);
0238   double stheta = sin(theta);
0239   double ctheta = cos(theta);
0240
0241   // some constants from kinematics
0242   double c1 = 0.5 * thePrimaryMass / theSecondaryMass;
0243   double c2 = sqrt(c1 * c1 - 1.);
0244   double c3 = 0.5 * c2 * ctheta / c1;
0245   double c4 = sqrt(p2 + thePrimaryMass * thePrimaryMass);
0246
0247   // momentum of decay particle 1 in the rest frame of the primary
0248   AlgebraicMatrix pplus(3, 1);
0249   pplus[0][0] = theSecondaryMass * c2 * stheta * cphi;
0250   pplus[1][0] = theSecondaryMass * c2 * stheta * sphi;
0251   pplus[2][0] = 0.5 * p + c3 * c4;
0252
0253   // momentum of decay particle 2 in the rest frame of the primary
0254   AlgebraicMatrix pminus(3, 1);
0255   pminus[0][0] = -pplus[0][0];
0256   pminus[1][0] = -pplus[1][0];
0257   pminus[2][0] = 0.5 * p - c3 * c4;
0258
0259   // derivative of rotation matrix w.r.t. py
0260   AlgebraicMatrix dRotMatdpz(3, 3);
0261
0262   dRotMatdpz[0][0] = px / (pT * p) * (1. - pz * pz / p2);
0263   dRotMatdpz[0][1] = 0.;
0264   dRotMatdpz[0][2] = -px * pz / (p * p2);
0265
0266   dRotMatdpz[1][0] = py / (p * pT) * (1. - pz * pz / p2);
0267   dRotMatdpz[1][1] = 0.;
0268   dRotMatdpz[1][2] = -py * pz / (p * p2);
0269
0270   dRotMatdpz[2][0] = pT * pz / (p * p2);
0271   dRotMatdpz[2][1] = 0.;
0272   dRotMatdpz[2][2] = (1. - pz * pz / p2) / p;
0273
0274   // derivative of the momentum of particle 1 in the lab frame w.r.t. pz
0275   double dpplusdpz = pz * (0.5 / p + c3 / c4);
0276
0277   AlgebraicMatrix dqplusdpz = dRotMatdpz * pplus;
0278   dqplusdpz[0][0] += px * dpplusdpz / p;
0279   dqplusdpz[1][0] += py * dpplusdpz / p;
0280   dqplusdpz[2][0] += pz * dpplusdpz / p;
0281
0282   // derivative of the momentum of particle 2 in the lab frame w.r.t. pz
0283   double dpminusdpz = pz * (0.5 / p - c3 / c4);
0284
0285   AlgebraicMatrix dqminusdpz = dRotMatdpz * pminus;
0286   dqminusdpz[0][0] += px * dpminusdpz / p;
0287   dqminusdpz[1][0] += py * dpminusdpz / p;
0288   dqminusdpz[2][0] += pz * dpminusdpz / p;
0289
0290   // return result
0291   return std::make_pair(dqplusdpz, dqminusdpz);
0292 }
0293
0294 const std::pair<AlgebraicMatrix, AlgebraicMatrix> TwoBodyDecayDerivatives::dqsdtheta(
0295     const TwoBodyDecayParameters &param) const {
0296   double px = param[TwoBodyDecayParameters::px];
0297   double py = param[TwoBodyDecayParameters::py];
0298   double pz = param[TwoBodyDecayParameters::pz];
0299   double theta = param[TwoBodyDecayParameters::theta];
0300   double phi = param[TwoBodyDecayParameters::phi];
0301
0302   // compute transverse and absolute momentum
0303   double pT2 = px * px + py * py;
0304   double p2 = pT2 + pz * pz;
0305
0306   double sphi = sin(phi);
0307   double cphi = cos(phi);
0308   double stheta = sin(theta);
0309   double ctheta = cos(theta);
0310
0311   // some constants from kinematics
0312   double c1 = 0.5 * thePrimaryMass / theSecondaryMass;
0313   double c2 = sqrt(c1 * c1 - 1.);
0314   double c3 = -0.5 * c2 * stheta / c1;
0315   double c4 = sqrt(p2 + thePrimaryMass * thePrimaryMass);
0316
0317   // derivative of the momentum of particle 1 in the primary's rest frame w.r.t.
0318   // theta
0319   AlgebraicMatrix dpplusdtheta(3, 1);
0320   dpplusdtheta[0][0] = theSecondaryMass * c2 * ctheta * cphi;
0321   dpplusdtheta[1][0] = theSecondaryMass * c2 * ctheta * sphi;
0322   dpplusdtheta[2][0] = c3 * c4;
0323
0324   // derivative of the momentum of particle 2 in the primary's rest frame w.r.t.
0325   // theta
0326   AlgebraicMatrix dpminusdtheta(3, 1);
0327   dpminusdtheta[0][0] = -theSecondaryMass * c2 * ctheta * cphi;
0328   dpminusdtheta[1][0] = -theSecondaryMass * c2 * ctheta * sphi;
0329   dpminusdtheta[2][0] = -c3 * c4;
0330
0331   TwoBodyDecayModel decayModel;
0332   AlgebraicMatrix rotMat = decayModel.rotationMatrix(px, py, pz);
0333
0334   AlgebraicMatrix dqplusdtheta = rotMat * dpplusdtheta;
0335   AlgebraicMatrix dqminusdtheta = rotMat * dpminusdtheta;
0336
0337   return std::make_pair(dqplusdtheta, dqminusdtheta);
0338 }
0339
0340 const std::pair<AlgebraicMatrix, AlgebraicMatrix> TwoBodyDecayDerivatives::dqsdphi(
0341     const TwoBodyDecayParameters &param) const {
0342   double px = param[TwoBodyDecayParameters::px];
0343   double py = param[TwoBodyDecayParameters::py];
0344   double pz = param[TwoBodyDecayParameters::pz];
0345   double theta = param[TwoBodyDecayParameters::theta];
0346   double phi = param[TwoBodyDecayParameters::phi];
0347
0348   double sphi = sin(phi);
0349   double cphi = cos(phi);
0350   double stheta = sin(theta);
0351
0352   // some constants from kinematics
0353   double c1 = 0.5 * thePrimaryMass / theSecondaryMass;
0354   double c2 = sqrt(c1 * c1 - 1.);
0355
0356   // derivative of the momentum of particle 1 in the primary's rest frame w.r.t.
0357   // phi
0358   AlgebraicMatrix dpplusdphi(3, 1);
0359   dpplusdphi[0][0] = -theSecondaryMass * c2 * stheta * sphi;
0360   dpplusdphi[1][0] = theSecondaryMass * c2 * stheta * cphi;
0361   dpplusdphi[2][0] = 0.;
0362
0363   // derivative of the momentum of particle 2 in the primary's rest frame w.r.t.
0364   // phi
0365   AlgebraicMatrix dpminusdphi(3, 1);
0366   dpminusdphi[0][0] = theSecondaryMass * c2 * stheta * sphi;
0367   dpminusdphi[1][0] = -theSecondaryMass * c2 * stheta * cphi;
0368   dpminusdphi[2][0] = 0.;
0369
0370   TwoBodyDecayModel decayModel;
0371   AlgebraicMatrix rotMat = decayModel.rotationMatrix(px, py, pz);
0372
0373   AlgebraicMatrix dqplusdphi = rotMat * dpplusdphi;
0374   AlgebraicMatrix dqminusdphi = rotMat * dpminusdphi;
0375
0376   return std::make_pair(dqplusdphi, dqminusdphi);
0377 }
0378
0379 const std::pair<AlgebraicMatrix, AlgebraicMatrix> TwoBodyDecayDerivatives::dqsdm(
0380     const TwoBodyDecayParameters &param) const {
0381   double px = param[TwoBodyDecayParameters::px];
0382   double py = param[TwoBodyDecayParameters::py];
0383   double pz = param[TwoBodyDecayParameters::pz];
0384   double theta = param[TwoBodyDecayParameters::theta];
0385   double phi = param[TwoBodyDecayParameters::phi];
0386
0387   double pT2 = px * px + py * py;
0388   double p2 = pT2 + pz * pz;
0389
0390   double sphi = sin(phi);
0391   double cphi = cos(phi);
0392   double ctheta = cos(theta);
0393   double stheta = sin(theta);
0394
0395   // some constants from kinematics
0396   double c1 = 0.5 * thePrimaryMass / theSecondaryMass;
0397   double c2 = 1. / sqrt(c1 * c1 - 1.);
0398   double m2 = thePrimaryMass * thePrimaryMass;
0399
0400   // derivative of the momentum of particle 1 in the primary's rest frame w.r.t.
0401   // the primary's mass
0402   AlgebraicMatrix dpplusdm(3, 1);
0403   dpplusdm[0][0] = c2 * 0.5 * c1 * stheta * cphi;
0404   dpplusdm[1][0] = c2 * 0.5 * c1 * stheta * sphi;
0405   dpplusdm[2][0] = c2 * theSecondaryMass * (c1 * c1 + p2 / m2) / sqrt(p2 + m2) * ctheta;
0406
0407   // derivative of the momentum of particle 2 in the primary's rest frame w.r.t.
0408   // the primary's mass
0409   AlgebraicMatrix dpminusdm(3, 1);
0410   dpminusdm[0][0] = -dpplusdm[0][0];
0411   dpminusdm[1][0] = -dpplusdm[1][0];
0412   dpminusdm[2][0] = -dpplusdm[2][0];
0413
0414   TwoBodyDecayModel decayModel;
0415   AlgebraicMatrix rotMat = decayModel.rotationMatrix(px, py, pz);
0416
0417   AlgebraicMatrix dqplusdm = rotMat * dpplusdm;
0418   AlgebraicMatrix dqminusdm = rotMat * dpminusdm;
0419
0420   return std::make_pair(dqplusdm, dqminusdm);
0421 }
0422
0423 const std::pair<AlgebraicMatrix, AlgebraicMatrix> TwoBodyDecayDerivatives::dqsdzi(
0424     const TwoBodyDecayParameters &param, const DerivativeParameterName &i) const {
0425   switch (i) {
0426     case TwoBodyDecayDerivatives::px:
0427       return dqsdpx(param);
0428       break;
0429     case TwoBodyDecayDerivatives::py:
0430       return dqsdpy(param);
0431       break;
0432     case TwoBodyDecayDerivatives::pz:
0433       return dqsdpz(param);
0434       break;
0435     case TwoBodyDecayDerivatives::theta:
0436       return dqsdtheta(param);
0437       break;
0438     case TwoBodyDecayDerivatives::phi:
0439       return dqsdphi(param);
0440       break;
0441     case TwoBodyDecayDerivatives::mass:
0442       return dqsdm(param);
0443       break;
0444     default:
0445       throw cms::Exception("BadConfig") << "@SUB=TwoBodyDecayDerivatives::dqsdzi"
0446                                         << "no decay parameter related to selector (" << i << ").";
0447   };
0448
0449   return std::make_pair(AlgebraicMatrix(3, 1, 0), AlgebraicMatrix(3, 1, 0));
0450 }