File indexing completed on 2023-10-25 09:33:09
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
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
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 ¶m) {
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
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
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
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
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 }