RunHitFit.h

Go to the documentation of this file.

//
//

#ifndef HITFIT_RUNHITFIT_H
#define HITFIT_RUNHITFIT_H

#include <algorithm>

#include "TopQuarkAnalysis/TopHitFit/interface/Defaults_Text.h"
#include "TopQuarkAnalysis/TopHitFit/interface/LeptonTranslatorBase.h"
#include "TopQuarkAnalysis/TopHitFit/interface/JetTranslatorBase.h"
#include "TopQuarkAnalysis/TopHitFit/interface/METTranslatorBase.h"
#include "TopQuarkAnalysis/TopHitFit/interface/Fit_Result.h"
#include "TopQuarkAnalysis/TopHitFit/interface/Lepjets_Event.h"
#include "TopQuarkAnalysis/TopHitFit/interface/Top_Fit.h"

// Explanation about the MIN/MAX definitions:
//
// For a given number of jets, there is a corresponding number of
// permutations how to assign each jet in the event to the corresponding
// parton-level jet.
// The number of permutations up to 10 jets are given below for Tt and
// TtH events.
//
// NJet         Npermutation (Tt)       Npermutation (TtH)
// 4            24                      --
// 5            120                     --
// 6            360                     360
// 7            840                     2520
// 8            1680                    10080
// 9            3024                    30240
// 10           5040                    75600
//
// The formulas for the number of permutations for Tt and TtH events
// given n jets in the event are
//
//         n!
// Tt:  -------- ; n >= 4
//      (n - 4)!
//
//          n!
// TtH: ---------- ; n >= 6
//      (n - 6)!2!
//
// The current MAX settings are chosen for a maximum number of 8 jets
// Increasing this limit should be done with caution, as it will
// increase the number of permutations rapidly.
//

namespace hitfit {

  template <class AElectron, class AMuon, class AJet, class AMet>
  class RunHitFit {
  private:
    LeptonTranslatorBase<AElectron> _ElectronTranslator;

    LeptonTranslatorBase<AMuon> _MuonTranslator;

    JetTranslatorBase<AJet> _JetTranslator;

    METTranslatorBase<AMet> _METTranslator;

    Lepjets_Event _event;

    std::vector<AJet> _jets;

    bool _jetObjRes;

    Top_Fit _Top_Fit;

    std::vector<Lepjets_Event> _Unfitted_Events;

    std::vector<Fit_Result> _Fit_Results;

  public:
    RunHitFit(const LeptonTranslatorBase<AElectron>& el,
              const LeptonTranslatorBase<AMuon>& mu,
              const JetTranslatorBase<AJet>& jet,
              const METTranslatorBase<AMet>& met,
              const std::string default_file,
              double lepw_mass,
              double hadw_mass,
              double top_mass)
        : _ElectronTranslator(el),
          _MuonTranslator(mu),
          _JetTranslator(jet),
          _METTranslator(met),
          _event(0, 0),
          _jetObjRes(false),
          _Top_Fit(Top_Fit_Args(Defaults_Text(default_file)), lepw_mass, hadw_mass, top_mass) {}

    ~RunHitFit() {}

    void clear() {
      _event = Lepjets_Event(0, 0);
      _jets.clear();
      _jetObjRes = false;
      _Unfitted_Events.clear();
      _Fit_Results.clear();
    }

    void AddLepton(const AElectron& electron, bool useObjRes = false) {
      _event.add_lep(_ElectronTranslator(electron, electron_label, useObjRes));
      return;
    }

    void AddLepton(const AMuon& muon, bool useObjRes = false) {
      _event.add_lep(_MuonTranslator(muon, muon_label, useObjRes));
      return;
    }

    void AddJet(const AJet& jet, bool useObjRes = false) {
      // Only set flag when adding the first jet
      // the additional jets then WILL be treated in the
      // same way like the first jet.
      if (_jets.empty()) {
        _jetObjRes = useObjRes;
      }

      if (_jets.size() < MAX_HITFIT_JET) {
        _jets.push_back(jet);
      }
      return;
    }

    void SetMet(const AMet& met, bool useObjRes = false) {
      _event.met() = _METTranslator(met, useObjRes);
      _event.kt_res() = _METTranslator.KtResolution(met, useObjRes);
      return;
    }

    void SetKtResolution(const Resolution& res) {
      _event.kt_res() = res;
      return;
    }

    void SetMETResolution(const Resolution& res) {
      SetKtResolution(res);
      return;
    }

    const Top_Fit& GetTopFit() const { return _Top_Fit; }

    std::vector<Fit_Result>::size_type FitAllPermutation() {
      if (_jets.size() < MIN_HITFIT_JET) {
        // For ttbar lepton+jets, a minimum of MIN_HITFIT_JETS jets
        // is required
        return 0;
      }

      if (_jets.size() > MAX_HITFIT_JET) {
        // Restrict the maximum number of jets in the fit
        // to prevent loop overflow
        return 0;
      }

      _Unfitted_Events.clear();
      _Fit_Results.clear();

      // Prepare the array of jet types for permutation
      std::vector<int> jet_types(_jets.size(), unknown_label);
      jet_types[0] = lepb_label;
      jet_types[1] = hadb_label;
      jet_types[2] = hadw1_label;
      jet_types[3] = hadw1_label;

      if (_Top_Fit.args().do_higgs_flag() && _jets.size() >= MIN_HITFIT_TTH) {
        jet_types[4] = higgs_label;
        jet_types[5] = higgs_label;
      }

      std::stable_sort(jet_types.begin(), jet_types.end());

      do {
        // begin loop over all jet permutation
        for (int nusol = 0; nusol != 2; nusol++) {
          // loop over two neutrino solution
          bool nuz = bool(nusol);

          // Copy the event
          Lepjets_Event fev = _event;

          // Add jets into the event, with the assumed type
          // in accord with the permutation.
          // The translator _JetTranslator will correctly
          // return object of Lepjets_Event_Jet with
          // jet energy correction applied in accord with
          // the assumed jet type (b or light).
          for (size_t j = 0; j != _jets.size(); j++) {
            fev.add_jet(_JetTranslator(_jets[j], jet_types[j], _jetObjRes));
          }

          // Clone fev (intended to be fitted event)
          // to ufev (intended to be unfitted event)
          Lepjets_Event ufev = fev;

          // Set jet types.
          fev.set_jet_types(jet_types);
          ufev.set_jet_types(jet_types);

          // Store the unfitted event
          _Unfitted_Events.push_back(ufev);

          // Prepare the placeholder for various kinematic quantities
          double umwhad;
          double utmass;
          double mt;
          double sigmt;
          Column_Vector pullx;
          Column_Vector pully;

          // Do the fit
          double chisq = _Top_Fit.fit_one_perm(fev, nuz, umwhad, utmass, mt, sigmt, pullx, pully);
          // Store output of the fit
          _Fit_Results.push_back(Fit_Result(chisq, fev, pullx, pully, umwhad, utmass, mt, sigmt));

        }  // end loop over two neutrino solution

      } while (std::next_permutation(jet_types.begin(), jet_types.end()));
      // end loop over all jet permutations

      return _Fit_Results.size();
    }

    std::vector<Lepjets_Event> GetUnfittedEvent() { return _Unfitted_Events; }

    std::vector<Fit_Result> GetFitAllPermutation() { return _Fit_Results; }

    static const unsigned int MIN_HITFIT_JET = 4;

    static const unsigned int MIN_HITFIT_TTH = 6;

    static const unsigned int MAX_HITFIT_JET = 8;

    static const unsigned int MAX_HITFIT = 1680;

    static const unsigned int MAX_HITFIT_VAR = 32;
  };

}  // namespace hitfit

#endif  // #ifndef RUNHITFIT_H