Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-03-17 11:26:09

0001 //
0002 //
0003 // File: src/Top_Fit.cc
0004 // Purpose: Handle jet permutations.
0005 // Created: Jul, 2000, sss, based on run 1 mass analysis code.
0006 //
0007 // XXX handle merging jets.
0008 // XXX btagging for ttH.
0009 //
0010 // CMSSW File      : src/Top_Fit.cc
0011 // Original Author : Scott Stuart Snyder <snyder@bnl.gov> for D0
0012 // Imported to CMSSW by Haryo Sumowidagdo <Suharyo.Sumowidagdo@cern.ch>
0013 //
0014 
0015 /**
0016     @file Top_Fit.cc
0017 
0018     @brief Handle and fit jet permutations of an event.  This is the
0019     primary interface between user's Lepjets_Event and HitFit kinematic
0020     fitting algorithm.  See the documentation for the header file
0021     Top_Fit.h for details.
0022 
0023     @author Scott Stuart Snyder <snyder@bnl.gov>
0024 
0025     @par Creation date:
0026     Jul 2000.
0027 
0028     @par Modification History:
0029     Apr 2009: Haryo Sumowidagdo <Suharyo.Sumowidagdo@cern.ch>:
0030     Imported to CMSSW.<br>
0031     Nov 2009: Haryo Sumowidagdo <Suharyo.Sumowidagdo@cern.ch>:
0032     Added doxygen tags for automatic generation of documentation.
0033 
0034     @par Terms of Usage:
0035     With consent for the original author (Scott Snyder).
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   // Argument handling.
0062   //
0063 
0064   Top_Fit_Args::Top_Fit_Args(const Defaults& defs)
0065       //
0066       // Purpose: Constructor.
0067       //
0068       // Inputs:
0069       //   defs -        The Defaults instance from which to initialize.
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   // Purpose: Return the print_event_flag parameter.
0084   //          See the header for documentation.
0085   //
0086   {
0087     return _print_event_flag;
0088   }
0089 
0090   bool Top_Fit_Args::do_higgs_flag() const
0091   //
0092   // Purpose: Return the do_higgs_flag parameter.
0093   //          See the header for documentation.
0094   //
0095   {
0096     return _do_higgs_flag;
0097   }
0098 
0099   double Top_Fit_Args::jet_mass_cut() const
0100   //
0101   // Purpose: Return the jet_mass_cut parameter.
0102   //          See the header for documentation.
0103   //
0104   {
0105     return _jet_mass_cut;
0106   }
0107 
0108   double Top_Fit_Args::mwhad_min_cut() const
0109   //
0110   // Purpose: Return the mwhad_min_cut parameter.
0111   //          See the header for documentation.
0112   //
0113   {
0114     return _mwhad_min_cut;
0115   }
0116 
0117   double Top_Fit_Args::mwhad_max_cut() const
0118   //
0119   // Purpose: Return the mwhad_max_cut parameter.
0120   //          See the header for documentation.
0121   //
0122   {
0123     return _mwhad_max_cut;
0124   }
0125 
0126   double Top_Fit_Args::mtdiff_max_cut() const
0127   //
0128   // Purpose: Return the mtdiff_max_cut parameter.
0129   //          See the header for documentation.
0130   //
0131   {
0132     return _mtdiff_max_cut;
0133   }
0134 
0135   int Top_Fit_Args::nkeep() const
0136   //
0137   // Purpose: Return the nkeep parameter.
0138   //          See the header for documentation.
0139   //
0140   {
0141     return _nkeep;
0142   }
0143 
0144   bool Top_Fit_Args::solve_nu_tmass() const
0145   //
0146   // Purpose: Return the solve_nu_tmass parameter
0147   //          See the header for documentation.
0148   //
0149   {
0150     return _solve_nu_tmass;
0151   }
0152 
0153   const Constrained_Top_Args& Top_Fit_Args::constrainer_args() const
0154   //
0155   // Purpose: Return the contained subobject parameters.
0156   //
0157   {
0158     return _args;
0159   }
0160 
0161   //*************************************************************************
0162   // Helper functions.
0163   //
0164 
0165   namespace {
0166 
0167     /**
0168     @brief Helper function: apply mass cuts to see if this
0169     event should be rejected before fitting.
0170 
0171     @param ev The event to test.
0172 
0173     @param args The parameter settings.
0174 
0175     @param mwhad  The hadronic  \f$ W- \f$ boson mass.
0176 
0177     @param umthad The mass of the hadronic top quark before fit.
0178 
0179     @param umtlep The mass of the leptonic top quark before fit.
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     // Purpose: Apply mass cuts to see if this event should be rejected
0185     //          without fitting.
0186     //
0187     // Inputs:
0188     //   ev -          The event to test.
0189     //   args -        Parameter setting.
0190     //   mwhad -       The hadronic W mass.
0191     //   umthad -      The hadronic top mass.
0192     //   umtlep -      The leptonic top mass.
0193     //
0194     // Returns:
0195     //   True if the event should be rejected.
0196     //
0197     {
0198       // Reject the event if any jet's mass is too large.
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       // Reject if if the hadronic W mass is outside the window.
0205       if (mwhad < args.mwhad_min_cut()) {
0206         return true;
0207       }
0208 
0209       // Reject if if the hadronic W mass is outside the window.
0210       if (mwhad > args.mwhad_max_cut()) {
0211         return true;
0212       }
0213 
0214       // And if the two top masses are too far apart.
0215       if (abs(umthad - umtlep) > args.mtdiff_max_cut()) {
0216         return true;
0217       }
0218 
0219       // It's ok.
0220       return false;
0221     }
0222 
0223     /**
0224     @brief Helper function: classify a jet permutation, to decide
0225     on what result lists it should be put.
0226 
0227     @param jet_types The vector representing a particular jet permutation,
0228     which is a vector of jet types.
0229 
0230     @param ev The original event being fit.
0231  */
0232     vector<int> classify_jetperm(const vector<int>& jet_types, const Lepjets_Event& ev)
0233     //
0234     // Purpose: Classify a jet permutation, to decide on what result
0235     //          lists it should be put.
0236     //
0237     // Inputs:
0238     //   jet_types -   Vector of jet types.
0239     //   ev -          The original event being fit.
0240     //
0241     // Returns:
0242     //   A list_flags vector, appropriate to pass to Fit_Results::push.
0243     //
0244     {
0245       // Start by assuming it's on all the lists.
0246       // We'll clear the flags if we see that it actually doesn't
0247       // belong.
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       // Loop over jets.
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];      // Current type of this jet.
0262           int t2 = ev.jet(i).type();  // `Correct' type of this jet.
0263 
0264           // Consider hadw1_label and hadw2_label the same.
0265           if (t1 == hadw2_label)
0266             t1 = hadw1_label;
0267           if (t2 == hadw2_label)
0268             t2 = hadw1_label;
0269 
0270           // If they're not the same, the permutation isn't correct.
0271           if (t1 != t2)
0272             out[noperm_list] = 0;
0273 
0274           // Test for a semicorrect permutation.
0275           // Here, all hadronic-side jets are considered equivalent.
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     @brief Helper function: update/overwrite the jet types in an event.
0302 
0303     @param jet_types The vector representing a particular jet permutation,
0304     which is a vector of jet types.
0305 
0306     @param ev Input: The event to update, output: the updated event.
0307  */
0308     void set_jet_types(const vector<int>& jet_types, Lepjets_Event& ev)
0309     //
0310     // Purpose: Update EV with a new set of jet types.
0311     //
0312     // Inputs:
0313     //   jet_types -   Vector of new jet types.
0314     //   ev -          The event to update.
0315     //
0316     // Outputs:
0317     //   ev -          The updated event.
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   }  // unnamed namespace
0334 
0335   //*************************************************************************
0336 
0337   Top_Fit::Top_Fit(const Top_Fit_Args& args, double lepw_mass, double hadw_mass, double top_mass)
0338       //
0339       // Purpose: Constructor.
0340       //
0341       // Inputs:
0342       //   args -        The parameter settings for this instance.
0343       //   lepw_mass -   The mass to which the leptonic W should be constrained,
0344       //                 or 0 to skip this constraint.
0345       //   hadw_mass -   The mass to which the hadronic W should be constrained,
0346       //                 or 0 to skip this constraint.
0347       //   top_mass -    The mass to which the top quarks should be constrained,
0348       //                 or 0 to skip this constraint.
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   // Purpose: Fit a single jet permutation.
0365   //
0366   // Inputs:
0367   //   ev -          The event to fit.
0368   //                 The object labels must have already been assigned.
0369   //   nuz -         Boolean flag to indicate which neutrino solution to be
0370   //                 used.
0371   //                 false = use smaller neutrino z solution
0372   //                 true  = use larger neutrino z solution
0373   //
0374   // Outputs:
0375   //   ev-           The event after the fit.
0376   //   umwhad -      Hadronic W mass before fitting.
0377   //   utmass -      Top mass before fitting, averaged from both sides.
0378   //   mt -          Top mass after fitting.
0379   //   sigmt -       Top mass uncertainty after fitting.
0380   //   pullx -       Vector of pull quantities for well-measured variables.
0381   //   pully -       Vector of pull quantities for poorly-measured variables.
0382   //
0383   // Returns:
0384   //   The fit chisq, or < 0 if the fit didn't converge.
0385   //
0386   // Adaptation note by Haryo Sumowidagdo:
0387   //   This function is rewritten in order to make its purpose reflects
0388   //   the function's name.  The function nows only fit one jet permutation
0389   //   with one neutrino solution only.
0390   //
0391   //
0392   {
0393     mt = 0;
0394     sigmt = 0;
0395 
0396     // Find the neutrino solutions by requiring either:
0397     // 1) that the leptonic top have the same mass as the hadronic top.
0398     // 2) that the mass of the lepton and neutrino is equal to the W mass
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     // Set up to use the selected neutrino solution
0411     if (!nuz) {
0412       ev.met().setZ(nuz1);
0413     } else {
0414       ev.met().setZ(nuz2);
0415     }
0416 
0417     // Note: We have set the neutrino Pz, but we haven't set the neutrino energy.
0418     // Remember that originally the neutrino energy was equal to
0419     // sqrt(nu_px*nu_px + nu_py*nu_py).  Calculating the invariant mass squared
0420     // for the neutrino will give negative mass squared.
0421     // Therefore we need to adjust (increase) the neutrino energy in order to
0422     // make its mass remain zero.
0423 
0424     adjust_e_for_mass(ev.met(), 0);
0425 
0426     // Find the unfit top mass as the average of the two sides.
0427     double umtlep = Top_Decaykin::lept(ev).m();
0428     utmass = (umthad + umtlep) / 2;
0429 
0430     // Trace, if requested.
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     // Maybe reject this event.
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     // Do the fit.
0443     double chisq = _constrainer.constrain(ev, mt, sigmt, pullx, pully);
0444 
0445     // Trace, if requested.
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     // Done!
0453     return chisq;
0454   }
0455 
0456   Fit_Results Top_Fit::fit(const Lepjets_Event& ev)
0457   //
0458   // Purpose: Fit all jet permutations for EV.
0459   //
0460   // Inputs:
0461   //   ev -          The event to fit.
0462   //
0463   // Returns:
0464   //   The results of the fit.
0465   //
0466   {
0467     // Make a new Fit_Results object.
0468     Fit_Results res(_args.nkeep(), n_lists);
0469 
0470     // Set up the vector of jet types.
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     // Must be in sorted order.
0484     stable_sort(jet_types.begin(), jet_types.end());
0485 
0486     do {
0487       // Loop over the two possible neutrino solution
0488       for (int nusol = 0; nusol != 2; nusol++) {
0489         // Set up the neutrino solution to be used
0490         bool nuz = bool(nusol);
0491 
0492         // Copy the event.
0493         Lepjets_Event fev = ev;
0494 
0495         // Install the new jet types.
0496         set_jet_types(jet_types, fev);
0497 
0498         // Figure out on what lists this permutation should go.
0499         vector<int> list_flags = classify_jetperm(jet_types, ev);
0500 
0501         // Set up the output variables for fit results.
0502         double umwhad, utmass, mt, sigmt;
0503         Column_Vector pullx;
0504         Column_Vector pully;
0505         double chisq;
0506 
0507         // Tracing.
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         // Do the fit.
0518         chisq = fit_one_perm(fev, nuz, umwhad, utmass, mt, sigmt, pullx, pully);
0519 
0520         // Print the result, if requested.
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         // Add it to the results.
0530         res.push(chisq, fev, pullx, pully, umwhad, utmass, mt, sigmt, list_flags);
0531 
0532       }  // end of for loop over the two neutrino solution
0533 
0534       // Step to the next permutation.
0535     } while (next_permutation(jet_types.begin(), jet_types.end()));
0536 
0537     return res;
0538   }
0539 
0540   /**
0541     @brief Output stream operator, print the content of this Top_Fit object
0542     to an output stream.
0543 
0544     @param s The output stream to which to write.
0545 
0546     @param fitter The instance of Top_Fit to be printed.
0547  */
0548   std::ostream& operator<<(std::ostream& s, const Top_Fit& fitter)
0549   //
0550   // Purpose: Print the object to S.
0551   //
0552   // Inputs:
0553   //   s -           The stream to which to write.
0554   //   fitter -      The object to write.
0555   //
0556   // Returns:
0557   //   The stream S.
0558   //
0559   {
0560     return s << fitter._constrainer;
0561   }
0562 
0563   const Top_Fit_Args& Top_Fit::args() const { return _args; }
0564 
0565 }  // namespace hitfit