/data/refman/pasoursint/CMSSW_5_2_9/src/TopQuarkAnalysis/TopHitFit/src/Constrained_Top.cc

Go to the documentation of this file.

//
// $Id: Constrained_Top.cc,v 1.1 2011/05/26 09:46:59 mseidel Exp $
//
// File: src/Constrained_Top.cc
// Purpose: Do kinematic fitting for a ttbar -> ljets event.
// Created: Jul, 2000, sss, based on run 1 mass analysis code.
//
// CMSSW File      : src/Constrained_Top.cc
// Original Author : Scott Stuart Snyder <snyder@bnl.gov> for D0
// Imported to CMSSW by Haryo Sumowidagdo <Suharyo.Sumowidagdo@cern.ch>
//


#include "TopQuarkAnalysis/TopHitFit/interface/Constrained_Top.h"
#include "TopQuarkAnalysis/TopHitFit/interface/Fourvec_Event.h"
#include "TopQuarkAnalysis/TopHitFit/interface/Lepjets_Event.h"
#include "TopQuarkAnalysis/TopHitFit/interface/Defaults.h"
#include <ostream>
#include <cassert>
#include <stdio.h>


using std::ostream;


namespace hitfit {


//*************************************************************************
// Argument handling.
//


Constrained_Top_Args::Constrained_Top_Args (const Defaults& defs)
//
// Purpose: Constructor.
//
// Inputs:
//   defs -        The Defaults instance from which to initialize.
//
  : _bmass (defs.get_float ("bmass")),
    _fourvec_constrainer_args (defs),
    _equal_side(defs.get_bool("equal_side"))

{
}


double Constrained_Top_Args::bmass () const
//
// Purpose: Return the bmass parameter.
//          See the header for documentation.
//
{
  return _bmass;
}


const Fourvec_Constrainer_Args&
Constrained_Top_Args::fourvec_constrainer_args () const
//
// Purpose: Return the contained subobject parameters.
//
{
  return _fourvec_constrainer_args;
}


bool Constrained_Top_Args::equal_side () const
//
// Purpose: Return the equal_side parameter
//          See the header for documentation.
//
{
  return _equal_side;
}


//*************************************************************************


Constrained_Top::Constrained_Top (const Constrained_Top_Args& args,
                                  double lepw_mass,
                                  double hadw_mass,
                                  double top_mass)
//
// Purpose: Constructor.
//
// Inputs:
//   args -        The parameter settings for this instance.
//   lepw_mass -   The mass to which the leptonic W should be constrained,
//                 or 0 to skip this constraint.
//   hadw_mass -   The mass to which the hadronic W should be constrained,
//                 or 0 to skip this constraint.
//   top_mass -    The mass to which the top quarks should be constrained,
//                 or 0 to skip this constraint.
//
  : _args (args),
    _constrainer (args.fourvec_constrainer_args())
{
  char buf[256];
  if (lepw_mass > 0) {
    sprintf (buf, "(%d %d) = %f", nu_label, lepton_label, lepw_mass);
    _constrainer.add_constraint (buf);
  }

  if (hadw_mass > 0) {
    sprintf (buf, "(%d %d) = %f", hadw1_label, hadw2_label, hadw_mass);
    _constrainer.add_constraint (buf);
  }

  if (args.equal_side()) {
      sprintf (buf, "(%d %d %d) = (%d %d %d)",
               nu_label, lepton_label, lepb_label,
               hadw1_label, hadw2_label, hadb_label);
      _constrainer.add_constraint (buf);
  }

  if (top_mass > 0) {
    sprintf (buf, "(%d %d %d) = %f",
             hadw1_label, hadw2_label, hadb_label, top_mass);
    _constrainer.add_constraint (buf);
  }
  else {
    sprintf (buf, "(%d %d %d) = 0", hadw1_label, hadw2_label, hadb_label);
    _constrainer.mass_constraint (buf);
  }
}

namespace {


FE_Obj make_fe_obj (const Lepjets_Event_Lep& obj, double mass, int type)
//
// Purpose: Helper to create an object to put into the Fourvec_Event.
//
// Inputs:
//   obj -         The input object.
//   mass -        The mass to which it should be constrained.
//   type -        The type to assign it.
//
// Returns:
//   The constructed FE_Obj.
//
{
  return FE_Obj (obj.p(), mass, type,
                 obj.p_sigma(), obj.eta_sigma(), obj.phi_sigma(),
                 obj.res().p_res().inverse());
}


void do_import (const Lepjets_Event& ev, double bmass, Fourvec_Event& fe)
//
// Purpose: Convert from a Lepjets_Event to a Fourvec_Event.
//
// Inputs:
//   ev -          The input event.
//   bmass -       The mass to which b-jets should be fixed.
//
// Outputs:
//   fe -          The initialized Fourvec_Event.
//
{
  assert (ev.nleps() == 1);
  fe.add (make_fe_obj (ev.lep(0), 0, lepton_label));

  bool saw_lepb = false;
  bool saw_hadb = false;
  for (std::vector<Lepjets_Event_Jet>::size_type j=0; j < ev.njets(); j++) {
    if (ev.jet(j).type() == isr_label || ev.jet(j).type() == higgs_label)
      continue;
    double mass = 0;
    if (ev.jet(j).type() == lepb_label && !saw_lepb) {
      mass = bmass;
      saw_lepb = true;
    }
    else if (ev.jet(j).type() == hadb_label && !saw_hadb) {
      mass = bmass;
      saw_hadb = true;
    }
    fe.add (make_fe_obj (ev.jet(j), mass, ev.jet(j).type()));
  }

  fe.set_nu_p (ev.met());
  Fourvec kt = ev.kt ();
  fe.set_kt_error (ev.kt_res().sigma (kt.x()),
                   ev.kt_res().sigma (kt.y()),
                   0);
}


void do_export (const Fourvec_Event& fe, Lepjets_Event& ev)
//
// Purpose: Convert from a Fourvec_Event to a Lepjets_Event.
//
// Inputs:
//   fe -          The input event.
//   ev -          The original Lepjets_Event.
//
// Outputs:
//   ev -          The updated Lepjets_Event.
//
{
  ev.lep(0).p() = fe.obj(0).p;
  for (std::vector<Lepjets_Event_Jet>::size_type j=0, k=1; j < ev.njets(); j++) {
    if (ev.jet(j).type() == isr_label || ev.jet(j).type() == higgs_label)
      continue;
    ev.jet(j).p() = fe.obj(k++).p;
  }

  ev.met() = fe.nu();
}


} // unnamed namespace


double Constrained_Top::constrain (Lepjets_Event& ev,
                                   double& mt,
                                   double& sigmt,
                                   Column_Vector& pullx,
                                   Column_Vector& pully)
//
// Purpose: Do a constrained fit for EV.  Returns the top mass and
//          its error in MT and SIGMT, and the pull quantities in PULLX and
//          PULLY.  Returns the chisq; this will be < 0 if the fit failed
//          to converge.
//
// Inputs:
//   ev -          The event we're fitting.
//
// Outputs:
//   ev -          The fitted event.
//   mt -          Requested invariant mass.
//   sigmt -       Uncertainty on mt.
//   pullx -       Pull quantities for well-measured variables.
//   pully -       Pull quantities for poorly-measured variables.
//
// Returns:
//   The fit chisq, or < 0 if the fit didn't converge.
//
{
  Fourvec_Event fe;
  do_import (ev, _args.bmass (), fe);
  double chisq = _constrainer.constrain (fe, mt, sigmt, pullx, pully);
  do_export (fe, ev);

  return chisq;
}


std::ostream& operator<< (std::ostream& s, const Constrained_Top& ct)
//
// Purpose: Print the object to S.
//
// Inputs:
//   s -           The stream to which to write.
//   ct -          The object to write.
//
// Returns:
//   The stream S.
//
{
  return s << ct._constrainer;
}


} // namespace hitfit