File indexing completed on 2024-04-06 12:31:19
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/Top_Decaykin.h"
0036 #include "TopQuarkAnalysis/TopHitFit/interface/Lepjets_Event.h"
0037 #include "TopQuarkAnalysis/TopHitFit/interface/fourvec.h"
0038 #include <cmath>
0039 #include <algorithm>
0040 #include <ostream>
0041
0042 using std::abs;
0043 using std::ostream;
0044 using std::sqrt;
0045 using std::swap;
0046
0047 namespace hitfit {
0048
0049 namespace {
0050
0051
0052
0053
0054
0055
0056 Fourvec leptons(const Lepjets_Event& ev)
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066 {
0067 return (ev.sum(lepton_label) + ev.sum(electron_label) + ev.sum(muon_label));
0068 }
0069
0070 }
0071
0072 bool Top_Decaykin::solve_nu_tmass(const Lepjets_Event& ev, double tmass, double& nuz1, double& nuz2)
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082
0083
0084
0085
0086
0087
0088
0089
0090 {
0091 bool discrim_flag = true;
0092
0093 const Fourvec& vnu = ev.met();
0094 Fourvec cprime = leptons(ev) + ev.sum(lepb_label);
0095 double alpha1 = tmass * tmass - cprime.m2();
0096 double a = 2 * 4 * (cprime.z() * cprime.z() - cprime.e() * cprime.e());
0097 double alpha = alpha1 + 2 * (cprime.x() * vnu.x() + cprime.y() * vnu.y());
0098 double b = 4 * alpha * cprime.z();
0099 double c = alpha * alpha - 4 * cprime.e() * cprime.e() * vnu.vect().perp2();
0100 double d = b * b - 2 * a * c;
0101 if (d < 0) {
0102 discrim_flag = false;
0103 d = 0;
0104 }
0105
0106 double dd = sqrt(d);
0107 nuz1 = (-b + dd) / a;
0108 nuz2 = (-b - dd) / a;
0109 if (abs(nuz1) > abs(nuz2))
0110 swap(nuz1, nuz2);
0111
0112 return discrim_flag;
0113 }
0114
0115 bool Top_Decaykin::solve_nu_tmass(
0116 const Lepjets_Event& ev, double tmass, double& re_nuz1, double& im_nuz1, double& re_nuz2, double& im_nuz2)
0117
0118
0119
0120
0121
0122
0123
0124
0125
0126
0127
0128
0129
0130
0131
0132
0133
0134
0135
0136
0137
0138 {
0139 bool discrim_flag = true;
0140
0141 const Fourvec& vnu = ev.met();
0142 Fourvec cprime = leptons(ev) + ev.sum(lepb_label);
0143 double alpha1 = tmass * tmass - cprime.m2();
0144 double a = 2 * 4 * (cprime.z() * cprime.z() - cprime.e() * cprime.e());
0145
0146
0147 double alpha = alpha1 + 2 * (cprime.x() * vnu.x() + cprime.y() * vnu.y());
0148 double b = 4 * alpha * cprime.z();
0149 double c = alpha * alpha - 4 * cprime.e() * cprime.e() * vnu.vect().perp2();
0150 double d = b * b - 2 * a * c;
0151 if (d < 0) {
0152 discrim_flag = false;
0153 }
0154
0155 if (discrim_flag) {
0156 re_nuz1 = (-b + sqrt(d)) / a;
0157 im_nuz1 = 0;
0158 re_nuz2 = (-b - sqrt(d)) / a;
0159 im_nuz2 = 0;
0160 if (abs(re_nuz1) > abs(re_nuz2)) {
0161 swap(re_nuz1, re_nuz2);
0162 }
0163
0164 } else {
0165
0166
0167
0168 re_nuz1 = -b / a;
0169 im_nuz1 = -fabs(sqrt(-d) / a);
0170 re_nuz2 = -b / a;
0171 im_nuz2 = fabs(sqrt(-d) / a);
0172 }
0173
0174 return discrim_flag;
0175 }
0176
0177 bool Top_Decaykin::solve_nu(const Lepjets_Event& ev, double wmass, double& nuz1, double& nuz2)
0178
0179
0180
0181
0182
0183
0184
0185
0186
0187
0188
0189
0190
0191
0192
0193
0194
0195 {
0196 bool discrim_flag = true;
0197
0198 const Fourvec& vnu = ev.met();
0199 Fourvec vlep = leptons(ev);
0200
0201 double x = vlep.x() * vnu.x() + vlep.y() * vnu.y() + wmass * wmass / 2;
0202 double a = vlep.z() * vlep.z() - vlep.e() * vlep.e();
0203 double b = 2 * x * vlep.z();
0204 double c = x * x - vnu.perp2() * vlep.e() * vlep.e();
0205
0206 double d = b * b - 4 * a * c;
0207 if (d < 0) {
0208 d = 0;
0209 discrim_flag = false;
0210 }
0211
0212 nuz1 = (-b + sqrt(d)) / 2 / a;
0213 nuz2 = (-b - sqrt(d)) / 2 / a;
0214 if (abs(nuz1) > abs(nuz2))
0215 swap(nuz1, nuz2);
0216
0217 return discrim_flag;
0218 }
0219
0220 bool Top_Decaykin::solve_nu(
0221 const Lepjets_Event& ev, double wmass, double& re_nuz1, double& im_nuz1, double& re_nuz2, double& im_nuz2)
0222
0223
0224
0225
0226
0227
0228
0229
0230
0231
0232
0233
0234
0235
0236
0237
0238
0239
0240
0241
0242
0243 {
0244 bool discrim_flag = true;
0245
0246 const Fourvec& vnu = ev.met();
0247 Fourvec vlep = leptons(ev);
0248
0249 double x = vlep.x() * vnu.x() + vlep.y() * vnu.y() + wmass * wmass / 2;
0250 double a = vlep.z() * vlep.z() - vlep.e() * vlep.e();
0251 double b = 2 * x * vlep.z();
0252 double c = x * x - vnu.perp2() * vlep.e() * vlep.e();
0253
0254 double d = b * b - 4 * a * c;
0255 if (d < 0) {
0256 discrim_flag = false;
0257 }
0258
0259 if (discrim_flag) {
0260 re_nuz1 = (-b + sqrt(d)) / 2 / a;
0261 im_nuz1 = 0.0;
0262 re_nuz2 = (-b - sqrt(d)) / 2 / a;
0263 im_nuz2 = 0.0;
0264 if (fabs(re_nuz1) > fabs(re_nuz2)) {
0265 swap(re_nuz1, re_nuz2);
0266 }
0267
0268 } else {
0269
0270
0271
0272
0273 re_nuz1 = -b / 2 / a;
0274 im_nuz1 = -fabs(sqrt(-d) / a);
0275 re_nuz2 = -b / 2 / a;
0276 im_nuz2 = fabs(sqrt(-d) / a);
0277 }
0278
0279 return discrim_flag;
0280 }
0281
0282 Fourvec Top_Decaykin::hadw(const Lepjets_Event& ev)
0283
0284
0285
0286
0287
0288
0289
0290
0291
0292 {
0293 return (ev.sum(hadw1_label) + ev.sum(hadw2_label));
0294 }
0295
0296 Fourvec Top_Decaykin::hadw1(const Lepjets_Event& ev)
0297
0298
0299
0300
0301
0302
0303
0304
0305
0306 {
0307 return ev.sum(hadw1_label);
0308 }
0309
0310 Fourvec Top_Decaykin::hadw2(const Lepjets_Event& ev)
0311
0312
0313
0314
0315
0316
0317
0318
0319
0320
0321 {
0322 return ev.sum(hadw2_label);
0323 }
0324
0325 Fourvec Top_Decaykin::lepw(const Lepjets_Event& ev)
0326
0327
0328
0329
0330
0331
0332
0333
0334
0335 {
0336 return (leptons(ev) + ev.met());
0337 }
0338
0339 Fourvec Top_Decaykin::hadt(const Lepjets_Event& ev)
0340
0341
0342
0343
0344
0345
0346
0347
0348
0349 {
0350 return (ev.sum(hadb_label) + hadw(ev));
0351 }
0352
0353 Fourvec Top_Decaykin::lept(const Lepjets_Event& ev)
0354
0355
0356
0357
0358
0359
0360
0361
0362
0363 {
0364 return (ev.sum(lepb_label) + lepw(ev));
0365 }
0366
0367 ostream& Top_Decaykin::dump_ev(std::ostream& s, const Lepjets_Event& ev)
0368
0369
0370
0371
0372
0373
0374
0375
0376
0377
0378 {
0379 s << ev;
0380 Fourvec p;
0381
0382 p = lepw(ev);
0383 s << "lepw " << p << " " << p.m() << "\n";
0384 p = lept(ev);
0385 s << "lept " << p << " " << p.m() << "\n";
0386 p = hadw(ev);
0387 s << "hadw " << p << " " << p.m() << "\n";
0388 p = hadt(ev);
0389 s << "hadt " << p << " " << p.m() << "\n";
0390
0391 return s;
0392 }
0393
0394 double Top_Decaykin::cos_theta_star(const Fourvec& fermion, const Fourvec& W, const Fourvec& top)
0395
0396
0397
0398
0399
0400
0401
0402
0403
0404
0405 {
0406 if (W.isLightlike() || W.isSpacelike()) {
0407 return 100.0;
0408 }
0409
0410 CLHEP::HepBoost BoostWCM(W.findBoostToCM());
0411
0412 CLHEP::Hep3Vector boost_v3fermion = BoostWCM(fermion).vect();
0413 CLHEP::Hep3Vector boost_v3top = BoostWCM(top).vect();
0414
0415 double costhetastar = boost_v3fermion.cosTheta(-boost_v3top);
0416
0417 return costhetastar;
0418 }
0419
0420 double Top_Decaykin::cos_theta_star(const Lepjets_Event& ev)
0421
0422
0423
0424
0425
0426
0427
0428
0429 {
0430 return cos_theta_star(leptons(ev), lepw(ev), lept(ev));
0431 }
0432
0433 double Top_Decaykin::cos_theta_star_lept(const Lepjets_Event& ev)
0434
0435
0436
0437
0438
0439
0440
0441
0442 {
0443 return cos_theta_star(ev);
0444 }
0445
0446 double Top_Decaykin::cos_theta_star_hadt(const Lepjets_Event& ev)
0447
0448
0449
0450
0451
0452
0453
0454
0455
0456 {
0457 return fabs(cos_theta_star(hadw1(ev), hadw(ev), hadt(ev)));
0458 }
0459
0460 }