File indexing completed on 2023-03-17 11:26:08
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035 #include "TopQuarkAnalysis/TopHitFit/interface/fourvec.h"
0036 #include <cmath>
0037
0038 using std::atan;
0039 using std::atan2;
0040 using std::cos;
0041 using std::exp;
0042 using std::fabs;
0043 using std::log;
0044 using std::sin;
0045 using std::sqrt;
0046 using std::tan;
0047
0048 namespace {
0049
0050 double cal_th(double theta, double z)
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062 {
0063 const double R_CC = 91.6;
0064 const double Z_EC = 178.9;
0065 const double BIGG = 1e8;
0066
0067 double tanth;
0068 if (fabs(cos(theta)) < 1 / BIGG)
0069 tanth = BIGG * sin(theta);
0070 else
0071 tanth = tan(theta);
0072
0073 double z_cc = R_CC / tanth + z;
0074
0075 if (fabs(z_cc) < Z_EC)
0076 theta = atan2(R_CC, z_cc);
0077
0078 else {
0079 double zz = Z_EC;
0080 if (tanth < 0)
0081 zz = -zz;
0082 double r_ec = (zz - z) * tanth;
0083 theta = atan2(r_ec, zz);
0084 }
0085
0086 if (theta < 0)
0087 theta += 2 * M_PI;
0088 return theta;
0089 }
0090
0091 }
0092
0093 namespace hitfit {
0094
0095 void adjust_p_for_mass(Fourvec& v, double mass)
0096
0097
0098
0099
0100
0101
0102
0103
0104
0105
0106
0107 {
0108 CLHEP::Hep3Vector vect = v.vect();
0109 double old_p2 = vect.mag2();
0110 if (old_p2 == 0)
0111 return;
0112 double new_p2 = v.e() * v.e() - mass * mass;
0113 if (new_p2 < 0)
0114 new_p2 = 0;
0115 vect *= sqrt(new_p2 / old_p2);
0116 v.setVect(vect);
0117 }
0118
0119 void adjust_e_for_mass(Fourvec& v, double mass)
0120
0121
0122
0123
0124
0125
0126
0127
0128
0129
0130
0131 {
0132 v.setE(sqrt(v.vect().mag2() + mass * mass));
0133 }
0134
0135 void rottheta(Fourvec& v, double theta)
0136
0137
0138
0139
0140
0141
0142
0143
0144
0145
0146 {
0147 double s = sin(theta), c = cos(theta);
0148 double old_pt = v.perp();
0149 double new_pt = old_pt * c - v.z() * s;
0150 v.setZ(old_pt * s + v.z() * c);
0151
0152 v.setX(v.x() * new_pt / old_pt);
0153 v.setY(v.y() * new_pt / old_pt);
0154 }
0155
0156 void roteta(Fourvec& v, double eta)
0157
0158
0159
0160
0161
0162
0163
0164
0165
0166
0167
0168 {
0169 double theta1 = v.theta();
0170 double eta1 = theta_to_eta(theta1);
0171 double eta2 = eta1 + eta;
0172 double theta2 = eta_to_theta(eta2);
0173
0174 rottheta(v, theta1 - theta2);
0175 }
0176
0177 double eta_to_theta(double eta)
0178
0179
0180
0181
0182
0183
0184
0185
0186
0187 {
0188 return 2 * atan(exp(-eta));
0189 }
0190
0191 double theta_to_eta(double theta)
0192
0193
0194
0195
0196
0197
0198
0199
0200
0201 {
0202 return -log(tan(theta / 2));
0203 }
0204
0205 double deteta(const Fourvec& v, double zvert)
0206
0207
0208
0209
0210
0211
0212
0213
0214
0215
0216 {
0217 return theta_to_eta(cal_th(v.theta(), zvert));
0218 }
0219
0220 double phidiff(double phi)
0221
0222
0223
0224
0225
0226
0227
0228
0229
0230 {
0231 while (phi < -M_PI)
0232 phi += 2 * M_PI;
0233 while (phi > M_PI)
0234 phi -= 2 * M_PI;
0235 return phi;
0236 }
0237
0238 double delta_r(const Fourvec& a, const Fourvec& b)
0239
0240
0241
0242
0243
0244
0245
0246
0247
0248
0249 {
0250 double deta = a.pseudoRapidity() - b.pseudoRapidity();
0251 double dphi = phidiff(a.phi() - b.phi());
0252 return sqrt(deta * deta + dphi * dphi);
0253 }
0254
0255 }