File indexing completed on 2023-03-17 11:15:49
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 <iostream>
0036 #include <cmath>
0037 #include "PhysicsTools/Heppy/interface/mt2w_bisect.h"
0038
0039 using namespace std;
0040
0041 namespace heppy {
0042
0043 namespace mt2w_bisect {
0044
0045
0046
0047 const float mt2w::RELATIVE_PRECISION = 0.00001;
0048 const float mt2w::ABSOLUTE_PRECISION = 0.0;
0049 const float mt2w::MIN_MASS = 0.1;
0050 const float mt2w::ZERO_MASS = 0.0;
0051 const float mt2w::SCANSTEP = 0.1;
0052
0053 mt2w::mt2w(double upper_bound, double error_value, double scan_step) {
0054 solved = false;
0055 momenta_set = false;
0056 mt2w_b = 0.;
0057 this->upper_bound = upper_bound;
0058 this->error_value =
0059 error_value;
0060 this->scan_step = scan_step;
0061 }
0062
0063 double mt2w::get_mt2w() {
0064 if (!momenta_set) {
0065 cout << " Please set momenta first!" << endl;
0066 return error_value;
0067 }
0068
0069 if (!solved)
0070 mt2w_bisect();
0071 return mt2w_b;
0072 }
0073
0074 void mt2w::set_momenta(double *pl, double *pb1, double *pb2, double *pmiss) {
0075
0076
0077 set_momenta(
0078 pl[0], pl[1], pl[2], pl[3], pb1[0], pb1[1], pb1[2], pb1[3], pb2[0], pb2[1], pb2[2], pb2[3], pmiss[1], pmiss[2]);
0079 }
0080
0081 void mt2w::set_momenta(double El,
0082 double plx,
0083 double ply,
0084 double plz,
0085 double Eb1,
0086 double pb1x,
0087 double pb1y,
0088 double pb1z,
0089 double Eb2,
0090 double pb2x,
0091 double pb2y,
0092 double pb2z,
0093 double pmissx,
0094 double pmissy) {
0095 solved = false;
0096 momenta_set = true;
0097
0098 double msqtemp;
0099
0100
0101
0102 this->El = El;
0103 this->plx = plx;
0104 this->ply = ply;
0105 this->plz = plz;
0106
0107 Elsq = El * El;
0108
0109 msqtemp = El * El - plx * plx - ply * ply - plz * plz;
0110 if (msqtemp > 0.0) {
0111 mlsq = msqtemp;
0112 } else {
0113 mlsq = 0.0;
0114 }
0115 ml = sqrt(mlsq);
0116
0117
0118
0119 this->Eb1 = Eb1;
0120 this->pb1x = pb1x;
0121 this->pb1y = pb1y;
0122 this->pb1z = pb1z;
0123
0124 Eb1sq = Eb1 * Eb1;
0125
0126 msqtemp = Eb1 * Eb1 - pb1x * pb1x - pb1y * pb1y - pb1z * pb1z;
0127 if (msqtemp > 0.0) {
0128 mb1sq = msqtemp;
0129 } else {
0130 mb1sq = 0.0;
0131 }
0132 mb1 = sqrt(mb1sq);
0133
0134
0135
0136 this->Eb2 = Eb2;
0137 this->pb2x = pb2x;
0138 this->pb2y = pb2y;
0139 this->pb2z = pb2z;
0140
0141 Eb2sq = Eb2 * Eb2;
0142
0143 msqtemp = Eb2 * Eb2 - pb2x * pb2x - pb2y * pb2y - pb2z * pb2z;
0144 if (msqtemp > 0.0) {
0145 mb2sq = msqtemp;
0146 } else {
0147 mb2sq = 0.0;
0148 }
0149 mb2 = sqrt(mb2sq);
0150
0151
0152
0153 this->pmissx = pmissx;
0154 this->pmissy = pmissy;
0155
0156
0157
0158 mv = 0.0;
0159 mw = 80.4;
0160
0161
0162
0163 if (ABSOLUTE_PRECISION > 100. * RELATIVE_PRECISION)
0164 precision = ABSOLUTE_PRECISION;
0165 else
0166 precision = 100. * RELATIVE_PRECISION;
0167 }
0168
0169 void mt2w::mt2w_bisect() {
0170 solved = true;
0171 cout.precision(11);
0172
0173
0174 double mtop_high = upper_bound;
0175 double mtop_low;
0176
0177 if (mb1 >= mb2) {
0178 mtop_low = mw + mb1;
0179 } else {
0180 mtop_low = mw + mb2;
0181 }
0182
0183
0184
0185
0186
0187
0188 if (teco(mtop_high) == 0) {
0189 mtop_high = mtop_low;
0190 }
0191
0192
0193
0194 while (teco(mtop_high) == 0 && mtop_high < upper_bound + 2. * scan_step) {
0195 mtop_low = mtop_high;
0196 mtop_high = mtop_high + scan_step;
0197 }
0198
0199
0200 if (mtop_high > upper_bound) {
0201 mt2w_b = error_value;
0202 return;
0203 }
0204
0205
0206 while (mtop_high - mtop_low > precision) {
0207 double mtop_mid, teco_mid;
0208
0209 mtop_mid = (mtop_high + mtop_low) / 2.;
0210 teco_mid = teco(mtop_mid);
0211
0212 if (teco_mid == 0) {
0213 mtop_low = mtop_mid;
0214 } else {
0215 mtop_high = mtop_mid;
0216 }
0217 }
0218 mt2w_b = mtop_high;
0219 return;
0220 }
0221
0222
0223
0224 int mt2w::teco(double mtop) {
0225
0226
0227 if (mtop < mb1 + mw || mtop < mb2 + mw) {
0228 return 0;
0229 }
0230
0231
0232
0233 double ETb2sq = Eb2sq - pb2z * pb2z;
0234 double delta = (mtop * mtop - mw * mw - mb2sq) / (2. * ETb2sq);
0235
0236
0237
0238 double del1 = mw * mw - mv * mv - mlsq;
0239 double del2 = mtop * mtop - mw * mw - mb1sq - 2 * (El * Eb1 - plx * pb1x - ply * pb1y - plz * pb1z);
0240
0241
0242
0243 double aa = (El * pb1x - Eb1 * plx) / (Eb1 * plz - El * pb1z);
0244 double bb = (El * pb1y - Eb1 * ply) / (Eb1 * plz - El * pb1z);
0245 double cc = (El * del2 - Eb1 * del1) / (2. * Eb1 * plz - 2. * El * pb1z);
0246
0247
0248
0249
0250
0251
0252
0253
0254
0255 a1 = Eb1sq * (1. + aa * aa) - (pb1x + pb1z * aa) * (pb1x + pb1z * aa);
0256 b1 = Eb1sq * aa * bb - (pb1x + pb1z * aa) * (pb1y + pb1z * bb);
0257 c1 = Eb1sq * (1. + bb * bb) - (pb1y + pb1z * bb) * (pb1y + pb1z * bb);
0258 d1 = Eb1sq * aa * cc - (pb1x + pb1z * aa) * (pb1z * cc + del2 / 2.0);
0259 e1 = Eb1sq * bb * cc - (pb1y + pb1z * bb) * (pb1z * cc + del2 / 2.0);
0260 f1 = Eb1sq * (mv * mv + cc * cc) - (pb1z * cc + del2 / 2.0) * (pb1z * cc + del2 / 2.0);
0261
0262
0263
0264 double det1 = (a1 * (c1 * f1 - e1 * e1) - b1 * (b1 * f1 - d1 * e1) + d1 * (b1 * e1 - c1 * d1)) / (a1 + c1);
0265
0266 if (det1 > 0.0) {
0267 return 0;
0268 }
0269
0270
0271
0272 a2 = 1 - pb2x * pb2x / (ETb2sq);
0273 b2 = -pb2x * pb2y / (ETb2sq);
0274 c2 = 1 - pb2y * pb2y / (ETb2sq);
0275
0276
0277
0278 d2o = -delta * pb2x;
0279 e2o = -delta * pb2y;
0280 f2o = mw * mw - delta * delta * ETb2sq;
0281
0282 d2 = -d2o - a2 * pmissx - b2 * pmissy;
0283 e2 = -e2o - c2 * pmissy - b2 * pmissx;
0284 f2 = a2 * pmissx * pmissx + 2 * b2 * pmissx * pmissy + c2 * pmissy * pmissy + 2 * d2o * pmissx +
0285 2 * e2o * pmissy + f2o;
0286
0287
0288 double x0, h0, y0, r0;
0289 x0 = (c1 * d1 - b1 * e1) / (b1 * b1 - a1 * c1);
0290 h0 = (b1 * x0 + e1) * (b1 * x0 + e1) - c1 * (a1 * x0 * x0 + 2 * d1 * x0 + f1);
0291 if (h0 < 0.0) {
0292 return 0;
0293 }
0294 y0 = (-b1 * x0 - e1 + sqrt(h0)) / c1;
0295 r0 = a2 * x0 * x0 + 2 * b2 * x0 * y0 + c2 * y0 * y0 + 2 * d2 * x0 + 2 * e2 * y0 + f2;
0296 if (r0 < 0.0) {
0297 return 1;
0298 }
0299
0300
0301
0302 long double A4, A3, A2, A1, A0;
0303
0304 A4 = -4 * a2 * b1 * b2 * c1 + 4 * a1 * b2 * b2 * c1 + a2 * a2 * c1 * c1 + 4 * a2 * b1 * b1 * c2 -
0305 4 * a1 * b1 * b2 * c2 - 2 * a1 * a2 * c1 * c2 + a1 * a1 * c2 * c2;
0306
0307 A3 = (-4 * a2 * b2 * c1 * d1 + 8 * a2 * b1 * c2 * d1 - 4 * a1 * b2 * c2 * d1 - 4 * a2 * b1 * c1 * d2 +
0308 8 * a1 * b2 * c1 * d2 - 4 * a1 * b1 * c2 * d2 - 8 * a2 * b1 * b2 * e1 + 8 * a1 * b2 * b2 * e1 +
0309 4 * a2 * a2 * c1 * e1 - 4 * a1 * a2 * c2 * e1 + 8 * a2 * b1 * b1 * e2 - 8 * a1 * b1 * b2 * e2 -
0310 4 * a1 * a2 * c1 * e2 + 4 * a1 * a1 * c2 * e2) /
0311 Eb1;
0312
0313 A2 = (4 * a2 * c2 * d1 * d1 - 4 * a2 * c1 * d1 * d2 - 4 * a1 * c2 * d1 * d2 + 4 * a1 * c1 * d2 * d2 -
0314 8 * a2 * b2 * d1 * e1 - 8 * a2 * b1 * d2 * e1 + 16 * a1 * b2 * d2 * e1 + 4 * a2 * a2 * e1 * e1 +
0315 16 * a2 * b1 * d1 * e2 - 8 * a1 * b2 * d1 * e2 - 8 * a1 * b1 * d2 * e2 - 8 * a1 * a2 * e1 * e2 +
0316 4 * a1 * a1 * e2 * e2 - 4 * a2 * b1 * b2 * f1 + 4 * a1 * b2 * b2 * f1 + 2 * a2 * a2 * c1 * f1 -
0317 2 * a1 * a2 * c2 * f1 + 4 * a2 * b1 * b1 * f2 - 4 * a1 * b1 * b2 * f2 - 2 * a1 * a2 * c1 * f2 +
0318 2 * a1 * a1 * c2 * f2) /
0319 Eb1sq;
0320
0321 A1 = (-8 * a2 * d1 * d2 * e1 + 8 * a1 * d2 * d2 * e1 + 8 * a2 * d1 * d1 * e2 - 8 * a1 * d1 * d2 * e2 -
0322 4 * a2 * b2 * d1 * f1 - 4 * a2 * b1 * d2 * f1 + 8 * a1 * b2 * d2 * f1 + 4 * a2 * a2 * e1 * f1 -
0323 4 * a1 * a2 * e2 * f1 + 8 * a2 * b1 * d1 * f2 - 4 * a1 * b2 * d1 * f2 - 4 * a1 * b1 * d2 * f2 -
0324 4 * a1 * a2 * e1 * f2 + 4 * a1 * a1 * e2 * f2) /
0325 (Eb1sq * Eb1);
0326
0327 A0 = (-4 * a2 * d1 * d2 * f1 + 4 * a1 * d2 * d2 * f1 + a2 * a2 * f1 * f1 + 4 * a2 * d1 * d1 * f2 -
0328 4 * a1 * d1 * d2 * f2 - 2 * a1 * a2 * f1 * f2 + a1 * a1 * f2 * f2) /
0329 (Eb1sq * Eb1sq);
0330
0331 long double A3sq;
0332
0333
0334
0335
0336
0337
0338
0339 A3sq = A3 * A3;
0340
0341 long double B3, B2, B1, B0;
0342 B3 = 4 * A4;
0343 B2 = 3 * A3;
0344 B1 = 2 * A2;
0345 B0 = A1;
0346
0347 long double C2, C1, C0;
0348 C2 = -(A2 / 2 - 3 * A3sq / (16 * A4));
0349 C1 = -(3 * A1 / 4. - A2 * A3 / (8 * A4));
0350 C0 = -A0 + A1 * A3 / (16 * A4);
0351
0352 long double D1, D0;
0353 D1 = -B1 - (B3 * C1 * C1 / C2 - B3 * C0 - B2 * C1) / C2;
0354 D0 = -B0 - B3 * C0 * C1 / (C2 * C2) + B2 * C0 / C2;
0355
0356 long double E0;
0357 E0 = -C0 - C2 * D0 * D0 / (D1 * D1) + C1 * D0 / D1;
0358
0359 long double t1, t2, t3, t4, t5;
0360
0361 t1 = A4;
0362 t2 = A4;
0363 t3 = C2;
0364 t4 = D1;
0365 t5 = E0;
0366
0367
0368 int nsol;
0369 nsol = signchange_n(t1, t2, t3, t4, t5) - signchange_p(t1, t2, t3, t4, t5);
0370
0371
0372 if (nsol < 0)
0373 nsol = 0;
0374
0375 int out;
0376 if (nsol == 0) {
0377 out = 0;
0378 }
0379 else {
0380 out = 1;
0381 }
0382
0383 return out;
0384 }
0385
0386 inline int mt2w::signchange_n(long double t1, long double t2, long double t3, long double t4, long double t5) {
0387 int nsc;
0388 nsc = 0;
0389 if (t1 * t2 > 0)
0390 nsc++;
0391 if (t2 * t3 > 0)
0392 nsc++;
0393 if (t3 * t4 > 0)
0394 nsc++;
0395 if (t4 * t5 > 0)
0396 nsc++;
0397 return nsc;
0398 }
0399 inline int mt2w::signchange_p(long double t1, long double t2, long double t3, long double t4, long double t5) {
0400 int nsc;
0401 nsc = 0;
0402 if (t1 * t2 < 0)
0403 nsc++;
0404 if (t2 * t3 < 0)
0405 nsc++;
0406 if (t3 * t4 < 0)
0407 nsc++;
0408 if (t4 * t5 < 0)
0409 nsc++;
0410 return nsc;
0411 }
0412
0413 }
0414
0415 }