RunHitFit.h

Go to the documentation of this file.

//
//     $Id: RunHitFit.h,v 1.1 2011/05/26 09:46:53 mseidel Exp $
//

#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