File indexing completed on 2023-03-17 11:26:09
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
0036
0037
0038
0039 #include "TopQuarkAnalysis/TopHitFit/interface/Top_Fit.h"
0040 #include "TopQuarkAnalysis/TopHitFit/interface/Lepjets_Event.h"
0041 #include "TopQuarkAnalysis/TopHitFit/interface/Top_Decaykin.h"
0042 #include "TopQuarkAnalysis/TopHitFit/interface/Defaults.h"
0043 #include "TopQuarkAnalysis/TopHitFit/interface/Fit_Results.h"
0044 #include "TopQuarkAnalysis/TopHitFit/interface/fourvec.h"
0045 #include <iostream>
0046 #include <algorithm>
0047 #include <cmath>
0048 #include <cassert>
0049
0050 using std::abs;
0051 using std::cout;
0052 using std::endl;
0053 using std::next_permutation;
0054 using std::ostream;
0055 using std::stable_sort;
0056 using std::vector;
0057
0058 namespace hitfit {
0059
0060
0061
0062
0063
0064 Top_Fit_Args::Top_Fit_Args(const Defaults& defs)
0065
0066
0067
0068
0069
0070
0071 : _print_event_flag(defs.get_bool("print_event_flag")),
0072 _do_higgs_flag(defs.get_bool("do_higgs_flag")),
0073 _jet_mass_cut(defs.get_float("jet_mass_cut")),
0074 _mwhad_min_cut(defs.get_float("mwhad_min_cut")),
0075 _mwhad_max_cut(defs.get_float("mwhad_max_cut")),
0076 _mtdiff_max_cut(defs.get_float("mtdiff_max_cut")),
0077 _nkeep(defs.get_int("nkeep")),
0078 _solve_nu_tmass(defs.get_bool("solve_nu_tmass")),
0079 _args(defs) {}
0080
0081 bool Top_Fit_Args::print_event_flag() const
0082
0083
0084
0085
0086 {
0087 return _print_event_flag;
0088 }
0089
0090 bool Top_Fit_Args::do_higgs_flag() const
0091
0092
0093
0094
0095 {
0096 return _do_higgs_flag;
0097 }
0098
0099 double Top_Fit_Args::jet_mass_cut() const
0100
0101
0102
0103
0104 {
0105 return _jet_mass_cut;
0106 }
0107
0108 double Top_Fit_Args::mwhad_min_cut() const
0109
0110
0111
0112
0113 {
0114 return _mwhad_min_cut;
0115 }
0116
0117 double Top_Fit_Args::mwhad_max_cut() const
0118
0119
0120
0121
0122 {
0123 return _mwhad_max_cut;
0124 }
0125
0126 double Top_Fit_Args::mtdiff_max_cut() const
0127
0128
0129
0130
0131 {
0132 return _mtdiff_max_cut;
0133 }
0134
0135 int Top_Fit_Args::nkeep() const
0136
0137
0138
0139
0140 {
0141 return _nkeep;
0142 }
0143
0144 bool Top_Fit_Args::solve_nu_tmass() const
0145
0146
0147
0148
0149 {
0150 return _solve_nu_tmass;
0151 }
0152
0153 const Constrained_Top_Args& Top_Fit_Args::constrainer_args() const
0154
0155
0156
0157 {
0158 return _args;
0159 }
0160
0161
0162
0163
0164
0165 namespace {
0166
0167
0168
0169
0170
0171
0172
0173
0174
0175
0176
0177
0178
0179
0180
0181 bool test_for_bad_masses(
0182 const Lepjets_Event& ev, const Top_Fit_Args& args, double mwhad, double umthad, double umtlep)
0183
0184
0185
0186
0187
0188
0189
0190
0191
0192
0193
0194
0195
0196
0197 {
0198
0199 if (ev.sum(lepb_label).m() > args.jet_mass_cut() || ev.sum(hadb_label).m() > args.jet_mass_cut() ||
0200 ev.sum(hadw1_label).m() > args.jet_mass_cut() || ev.sum(hadw2_label).m() > args.jet_mass_cut()) {
0201 return true;
0202 }
0203
0204
0205 if (mwhad < args.mwhad_min_cut()) {
0206 return true;
0207 }
0208
0209
0210 if (mwhad > args.mwhad_max_cut()) {
0211 return true;
0212 }
0213
0214
0215 if (abs(umthad - umtlep) > args.mtdiff_max_cut()) {
0216 return true;
0217 }
0218
0219
0220 return false;
0221 }
0222
0223
0224
0225
0226
0227
0228
0229
0230
0231
0232 vector<int> classify_jetperm(const vector<int>& jet_types, const Lepjets_Event& ev)
0233
0234
0235
0236
0237
0238
0239
0240
0241
0242
0243
0244 {
0245
0246
0247
0248 vector<int> out(n_lists);
0249 out[all_list] = 1;
0250 out[noperm_list] = 1;
0251 out[semicorrect_list] = 1;
0252 out[limited_isr_list] = 1;
0253 out[topfour_list] = 1;
0254 out[btag_list] = 1;
0255 out[htag_list] = 1;
0256
0257
0258 assert(jet_types.size() == ev.njets());
0259 for (vector<int>::size_type i = 0; i < jet_types.size(); i++) {
0260 {
0261 int t1 = jet_types[i];
0262 int t2 = ev.jet(i).type();
0263
0264
0265 if (t1 == hadw2_label)
0266 t1 = hadw1_label;
0267 if (t2 == hadw2_label)
0268 t2 = hadw1_label;
0269
0270
0271 if (t1 != t2)
0272 out[noperm_list] = 0;
0273
0274
0275
0276 if (t1 == hadw1_label)
0277 t1 = hadb_label;
0278 if (t2 == hadw1_label)
0279 t2 = hadb_label;
0280 if (t1 != t2)
0281 out[semicorrect_list] = 0;
0282 }
0283
0284 if (jet_types[i] == isr_label && i <= 2)
0285 out[limited_isr_list] = 0;
0286
0287 if ((jet_types[i] == isr_label && i <= 3) || (jet_types[i] != isr_label && i >= 4))
0288 out[topfour_list] = 0;
0289
0290 if ((ev.jet(i).svx_tag() || ev.jet(i).slt_tag()) && !(jet_types[i] == hadb_label || jet_types[i] == lepb_label))
0291 out[btag_list] = 0;
0292
0293 if ((ev.jet(i).svx_tag() || ev.jet(i).slt_tag()) &&
0294 !(jet_types[i] == hadb_label || jet_types[i] == lepb_label || jet_types[i] == higgs_label))
0295 out[htag_list] = 0;
0296 }
0297 return out;
0298 }
0299
0300
0301
0302
0303
0304
0305
0306
0307
0308 void set_jet_types(const vector<int>& jet_types, Lepjets_Event& ev)
0309
0310
0311
0312
0313
0314
0315
0316
0317
0318
0319 {
0320 assert(ev.njets() == jet_types.size());
0321 bool saw_hadw1 = false;
0322 for (vector<int>::size_type i = 0; i < ev.njets(); i++) {
0323 int t = jet_types[i];
0324 if (t == hadw1_label) {
0325 if (saw_hadw1)
0326 t = hadw2_label;
0327 saw_hadw1 = true;
0328 }
0329 ev.jet(i).type() = t;
0330 }
0331 }
0332
0333 }
0334
0335
0336
0337 Top_Fit::Top_Fit(const Top_Fit_Args& args, double lepw_mass, double hadw_mass, double top_mass)
0338
0339
0340
0341
0342
0343
0344
0345
0346
0347
0348
0349
0350 : _args(args),
0351 _constrainer(args.constrainer_args(), lepw_mass, hadw_mass, top_mass),
0352 _lepw_mass(lepw_mass),
0353 _hadw_mass(hadw_mass) {}
0354
0355 double Top_Fit::fit_one_perm(Lepjets_Event& ev,
0356 bool& nuz,
0357 double& umwhad,
0358 double& utmass,
0359 double& mt,
0360 double& sigmt,
0361 Column_Vector& pullx,
0362 Column_Vector& pully)
0363
0364
0365
0366
0367
0368
0369
0370
0371
0372
0373
0374
0375
0376
0377
0378
0379
0380
0381
0382
0383
0384
0385
0386
0387
0388
0389
0390
0391
0392 {
0393 mt = 0;
0394 sigmt = 0;
0395
0396
0397
0398
0399
0400 umwhad = Top_Decaykin::hadw(ev).m();
0401 double umthad = Top_Decaykin::hadt(ev).m();
0402 double nuz1, nuz2;
0403
0404 if (_args.solve_nu_tmass()) {
0405 Top_Decaykin::solve_nu_tmass(ev, umthad, nuz1, nuz2);
0406 } else {
0407 Top_Decaykin::solve_nu(ev, _lepw_mass, nuz1, nuz2);
0408 }
0409
0410
0411 if (!nuz) {
0412 ev.met().setZ(nuz1);
0413 } else {
0414 ev.met().setZ(nuz2);
0415 }
0416
0417
0418
0419
0420
0421
0422
0423
0424 adjust_e_for_mass(ev.met(), 0);
0425
0426
0427 double umtlep = Top_Decaykin::lept(ev).m();
0428 utmass = (umthad + umtlep) / 2;
0429
0430
0431 if (_args.print_event_flag()) {
0432 cout << "Top_Fit::fit_one_perm() : Before fit:\n";
0433 Top_Decaykin::dump_ev(cout, ev);
0434 }
0435
0436
0437 if (_hadw_mass > 0 && test_for_bad_masses(ev, _args, umwhad, umthad, umtlep)) {
0438 cout << "Top_Fit: bad mass comb.\n";
0439 return -999;
0440 }
0441
0442
0443 double chisq = _constrainer.constrain(ev, mt, sigmt, pullx, pully);
0444
0445
0446 if (_args.print_event_flag()) {
0447 cout << "Top_Fit::fit_one_perm() : After fit:\n";
0448 cout << "chisq: " << chisq << " mt: " << mt << " ";
0449 Top_Decaykin::dump_ev(cout, ev);
0450 }
0451
0452
0453 return chisq;
0454 }
0455
0456 Fit_Results Top_Fit::fit(const Lepjets_Event& ev)
0457
0458
0459
0460
0461
0462
0463
0464
0465
0466 {
0467
0468 Fit_Results res(_args.nkeep(), n_lists);
0469
0470
0471 vector<int> jet_types(ev.njets(), isr_label);
0472 assert(ev.njets() >= 4);
0473 jet_types[0] = lepb_label;
0474 jet_types[1] = hadb_label;
0475 jet_types[2] = hadw1_label;
0476 jet_types[3] = hadw1_label;
0477
0478 if (_args.do_higgs_flag() && ev.njets() >= 6) {
0479 jet_types[4] = higgs_label;
0480 jet_types[5] = higgs_label;
0481 }
0482
0483
0484 stable_sort(jet_types.begin(), jet_types.end());
0485
0486 do {
0487
0488 for (int nusol = 0; nusol != 2; nusol++) {
0489
0490 bool nuz = bool(nusol);
0491
0492
0493 Lepjets_Event fev = ev;
0494
0495
0496 set_jet_types(jet_types, fev);
0497
0498
0499 vector<int> list_flags = classify_jetperm(jet_types, ev);
0500
0501
0502 double umwhad, utmass, mt, sigmt;
0503 Column_Vector pullx;
0504 Column_Vector pully;
0505 double chisq;
0506
0507
0508 cout << "Top_Fit::fit(): Before fit: (";
0509 for (vector<int>::size_type i = 0; i < jet_types.size(); i++) {
0510 if (i)
0511 cout << " ";
0512 cout << jet_types[i];
0513 }
0514 cout << " nuz = " << nuz;
0515 cout << ") " << std::endl;
0516
0517
0518 chisq = fit_one_perm(fev, nuz, umwhad, utmass, mt, sigmt, pullx, pully);
0519
0520
0521 if (_args.print_event_flag()) {
0522 cout << "Top_Fit::fit(): After fit:\n";
0523 char buf[256];
0524 sprintf(
0525 buf, "chisq: %8.3f mt: %6.2f pm %5.2f %c\n", chisq, mt, sigmt, (list_flags[noperm_list] ? '*' : ' '));
0526 cout << buf;
0527 }
0528
0529
0530 res.push(chisq, fev, pullx, pully, umwhad, utmass, mt, sigmt, list_flags);
0531
0532 }
0533
0534
0535 } while (next_permutation(jet_types.begin(), jet_types.end()));
0536
0537 return res;
0538 }
0539
0540
0541
0542
0543
0544
0545
0546
0547
0548 std::ostream& operator<<(std::ostream& s, const Top_Fit& fitter)
0549
0550
0551
0552
0553
0554
0555
0556
0557
0558
0559 {
0560 return s << fitter._constrainer;
0561 }
0562
0563 const Top_Fit_Args& Top_Fit::args() const { return _args; }
0564
0565 }