TwoBodyDecay/src/TwoBodyDecayModel.cc

0001
0002 #include "Alignment/TwoBodyDecay/interface/TwoBodyDecayModel.h"
0003
0004 TwoBodyDecayModel::TwoBodyDecayModel(double mPrimary, double mSecondary)
0005     : thePrimaryMass(mPrimary), theSecondaryMass(mSecondary) {}
0006
0007 TwoBodyDecayModel::~TwoBodyDecayModel() {}
0008
0009 AlgebraicMatrix TwoBodyDecayModel::rotationMatrix(double px, double py, double pz) {
0010   // compute transverse and absolute momentum
0011   double pT2 = px * px + py * py;
0012   double p2 = pT2 + pz * pz;
0013   double pT = sqrt(pT2);
0014   double p = sqrt(p2);
0015
0016   AlgebraicMatrix rotMat(3, 3);
0017
0018   // compute rotation matrix
0019   rotMat[0][0] = px * pz / pT / p;
0020   rotMat[0][1] = -py / pT;
0021   rotMat[0][2] = px / p;
0022
0023   rotMat[1][0] = py * pz / pT / p;
0024   rotMat[1][1] = px / pT;
0025   rotMat[1][2] = py / p;
0026
0027   rotMat[2][0] = -pT / p;
0028   rotMat[2][1] = 0.;
0029   rotMat[2][2] = pz / p;
0030
0031   return rotMat;
0032 }
0033
0034 AlgebraicMatrix TwoBodyDecayModel::curvilinearToCartesianJacobian(double rho,
0035                                                                   double theta,
0036                                                                   double phi,
0037                                                                   double zMagField) {
0038   double q = ((rho < 0) ? -1. : 1.);
0039   double conv = q * zMagField;
0040
0041   double stheta = sin(theta);
0042   double ctheta = cos(theta);
0043   double sphi = sin(phi);
0044   double cphi = cos(phi);
0045
0046   AlgebraicMatrix curv2cart(3, 3);
0047
0048   curv2cart[0][0] = -rho * cphi;
0049   curv2cart[0][1] = -rho * sphi;
0050   curv2cart[0][2] = 0.;
0051
0052   curv2cart[1][0] = cphi * stheta * ctheta;
0053   curv2cart[1][1] = sphi * stheta * ctheta;
0054   curv2cart[1][2] = -stheta * stheta;
0055
0056   curv2cart[2][0] = -sphi;
0057   curv2cart[2][1] = cphi;
0058   curv2cart[2][2] = 0.;
0059
0060   curv2cart *= rho / conv;
0061
0062   return curv2cart;
0063 }
0064
0065 AlgebraicMatrix TwoBodyDecayModel::curvilinearToCartesianJacobian(const AlgebraicVector &curv, double zMagField) {
0066   return this->curvilinearToCartesianJacobian(curv[0], curv[1], curv[2], zMagField);
0067 }
0068
0069 AlgebraicVector TwoBodyDecayModel::convertCurvilinearToCartesian(const AlgebraicVector &curv, double zMagField) {
0070   double rt = fabs(zMagField / curv[0]);
0071
0072   AlgebraicVector cart(3);
0073   cart[0] = rt * cos(curv[2]);
0074   cart[1] = rt * sin(curv[2]);
0075   cart[2] = rt / tan(curv[1]);
0076
0077   return cart;
0078 }
0079
0080 const std::pair<AlgebraicVector, AlgebraicVector> TwoBodyDecayModel::cartesianSecondaryMomenta(
0081     const AlgebraicVector &param) {
0082   double px = param[TwoBodyDecayParameters::px];
0083   double py = param[TwoBodyDecayParameters::py];
0084   double pz = param[TwoBodyDecayParameters::pz];
0085   double theta = param[TwoBodyDecayParameters::theta];
0086   double phi = param[TwoBodyDecayParameters::phi];
0087
0088   // compute transverse and absolute momentum
0089   double pT2 = px * px + py * py;
0090   double p2 = pT2 + pz * pz;
0091   double p = sqrt(p2);
0092
0093   double sphi = sin(phi);
0094   double cphi = cos(phi);
0095   double stheta = sin(theta);
0096   double ctheta = cos(theta);
0097
0098   // some constants from kinematics
0099   double c1 = 0.5 * thePrimaryMass / theSecondaryMass;
0100   double c2 = sqrt(c1 * c1 - 1.);
0101   double c3 = 0.5 * c2 * ctheta / c1;
0102   double c4 = sqrt(p2 + thePrimaryMass * thePrimaryMass);
0103
0104   // momentum of decay particle 1 in the primary's boosted frame
0105   AlgebraicMatrix pplus(3, 1);
0106   pplus[0][0] = theSecondaryMass * c2 * stheta * cphi;
0107   pplus[1][0] = theSecondaryMass * c2 * stheta * sphi;
0108   pplus[2][0] = 0.5 * p + c3 * c4;
0109
0110   // momentum of decay particle 2 in the primary's boosted frame
0111   AlgebraicMatrix pminus(3, 1);
0112   pminus[0][0] = -pplus[0][0];
0113   pminus[1][0] = -pplus[1][0];
0114   pminus[2][0] = 0.5 * p - c3 * c4;
0115
0116   AlgebraicMatrix rotMat = rotationMatrix(px, py, pz);
0117
0118   return std::make_pair(rotMat * pplus, rotMat * pminus);
0119 }
0120
0121 const std::pair<AlgebraicVector, AlgebraicVector> TwoBodyDecayModel::cartesianSecondaryMomenta(const TwoBodyDecay &tbd) {
0122   return cartesianSecondaryMomenta(tbd.parameters());
0123 }
0124
0125 const std::pair<AlgebraicVector, AlgebraicVector> TwoBodyDecayModel::cartesianSecondaryMomenta(
0126     const TwoBodyDecayParameters &tbdparam) {
0127   return cartesianSecondaryMomenta(tbdparam.parameters());
0128 }