TopHitFit/src/Constrained_Top.cc

0001 //
0002 //
0003 // File: src/Constrained_Top.cc
0004 // Purpose: Do kinematic fitting for a ttbar -> ljets event.
0005 // Created: Jul, 2000, sss, based on run 1 mass analysis code.
0006 //
0007 // CMSSW File      : src/Constrained_Top.cc
0008 // Original Author : Scott Stuart Snyder <snyder@bnl.gov> for D0
0009 // Imported to CMSSW by Haryo Sumowidagdo <Suharyo.Sumowidagdo@cern.ch>
0010 //
0011
0012 /**
0013
0014     @file Constrained_Top.cc
0015
0016     @brief Do a constrained kinematic fit of a \f$t\bar{t}\to\ell +
0017     \rm{jets}\f$ event.  See the documentation for the header file
0018     Constrained_Top.h for details.
0019
0020     @par Creation date:
0021     July 2000.
0022
0023     @author
0024     Scott Stuart Snyder <snyder@bnl.gov>.
0025
0026     @par Modification History:
0027     Apr 2009: Haryo Sumowidagdo <Suharyo.Sumowidagdo@cern.ch>:
0028     Imported to CMSSW.<br>
0029     Oct 2009: Haryo Sumowidagdo <Suharyo.Sumowidagdo@cern.ch>:
0030     Added Doxygen tags for automatic generation of documentation.
0031
0032     @par Terms of Usage:
0033     With consent from the original author (Scott Snyder).
0034
0035  */
0036
0037 #include "TopQuarkAnalysis/TopHitFit/interface/Constrained_Top.h"
0038 #include "TopQuarkAnalysis/TopHitFit/interface/Fourvec_Event.h"
0039 #include "TopQuarkAnalysis/TopHitFit/interface/Lepjets_Event.h"
0040 #include "TopQuarkAnalysis/TopHitFit/interface/Defaults.h"
0041 #include <ostream>
0042 #include <cassert>
0043 #include <cstdio>
0044
0045 using std::ostream;
0046
0047 namespace hitfit {
0048
0049   //*************************************************************************
0050   // Argument handling.
0051   //
0052
0053   Constrained_Top_Args::Constrained_Top_Args(const Defaults& defs)
0054       //
0055       // Purpose: Constructor.
0056       //
0057       // Inputs:
0058       //   defs -        The Defaults instance from which to initialize.
0059       //
0060       : _bmass(defs.get_float("bmass")),
0061         _fourvec_constrainer_args(defs),
0062         _equal_side(defs.get_bool("equal_side"))
0063
0064   {}
0065
0066   double Constrained_Top_Args::bmass() const
0067   //
0068   // Purpose: Return the bmass parameter.
0069   //          See the header for documentation.
0070   //
0071   {
0072     return _bmass;
0073   }
0074
0075   const Fourvec_Constrainer_Args& Constrained_Top_Args::fourvec_constrainer_args() const
0076   //
0077   // Purpose: Return the contained subobject parameters.
0078   //
0079   {
0080     return _fourvec_constrainer_args;
0081   }
0082
0083   bool Constrained_Top_Args::equal_side() const
0084   //
0085   // Purpose: Return the equal_side parameter
0086   //          See the header for documentation.
0087   //
0088   {
0089     return _equal_side;
0090   }
0091
0092   //*************************************************************************
0093
0094   Constrained_Top::Constrained_Top(const Constrained_Top_Args& args, double lepw_mass, double hadw_mass, double top_mass)
0095       //
0096       // Purpose: Constructor.
0097       //
0098       // Inputs:
0099       //   args -        The parameter settings for this instance.
0100       //   lepw_mass -   The mass to which the leptonic W should be constrained,
0101       //                 or 0 to skip this constraint.
0102       //   hadw_mass -   The mass to which the hadronic W should be constrained,
0103       //                 or 0 to skip this constraint.
0104       //   top_mass -    The mass to which the top quarks should be constrained,
0105       //                 or 0 to skip this constraint.
0106       //
0107       : _args(args), _constrainer(args.fourvec_constrainer_args()) {
0108     char buf[256];
0109     if (lepw_mass > 0) {
0110       sprintf(buf, "(%d %d) = %f", nu_label, lepton_label, lepw_mass);
0111       _constrainer.add_constraint(buf);
0112     }
0113
0114     if (hadw_mass > 0) {
0115       sprintf(buf, "(%d %d) = %f", hadw1_label, hadw2_label, hadw_mass);
0116       _constrainer.add_constraint(buf);
0117     }
0118
0119     if (args.equal_side()) {
0120       sprintf(buf, "(%d %d %d) = (%d %d %d)", nu_label, lepton_label, lepb_label, hadw1_label, hadw2_label, hadb_label);
0121       _constrainer.add_constraint(buf);
0122     }
0123
0124     if (top_mass > 0) {
0125       sprintf(buf, "(%d %d %d) = %f", hadw1_label, hadw2_label, hadb_label, top_mass);
0126       _constrainer.add_constraint(buf);
0127     } else {
0128       sprintf(buf, "(%d %d %d) = 0", hadw1_label, hadw2_label, hadb_label);
0129       _constrainer.mass_constraint(buf);
0130     }
0131   }
0132
0133   namespace {
0134
0135     /**
0136
0137     @brief Helper function to create an object to put into Fourvec_Event.
0138
0139     @param obj The input object.
0140
0141     @param mass The mass to which it should be constrained.
0142
0143     @param type The type to assign to it.
0144
0145     @par Return:
0146     The constructed <i>FE_Obj</i>.
0147
0148  */
0149     FE_Obj make_fe_obj(const Lepjets_Event_Lep& obj, double mass, int type)
0150     //
0151     // Purpose: Helper to create an object to put into the Fourvec_Event.
0152     //
0153     // Inputs:
0154     //   obj -         The input object.
0155     //   mass -        The mass to which it should be constrained.
0156     //   type -        The type to assign it.
0157     //
0158     // Returns:
0159     //   The constructed FE_Obj.
0160     //
0161     {
0162       return FE_Obj(obj.p(), mass, type, obj.p_sigma(), obj.eta_sigma(), obj.phi_sigma(), obj.res().p_res().inverse());
0163     }
0164
0165     /**
0166
0167     @brief Convert from a Lepjets_Event to a Fourvec_Event.
0168
0169     @param ev The input event.
0170
0171     @param bmass The mass to which b-jets should be fixed.
0172
0173     @param fe The output Fourvec_Event.
0174
0175     @par Input:
0176     - Lepjets_Event <i>ev</i>.
0177     - double <i>bmass</i>.
0178
0179     @par Output:
0180     - Fourvec_Event <i>fe</i>.
0181
0182  */
0183     void do_import(const Lepjets_Event& ev, double bmass, Fourvec_Event& fe)
0184     //
0185     // Purpose: Convert from a Lepjets_Event to a Fourvec_Event.
0186     //
0187     // Inputs:
0188     //   ev -          The input event.
0189     //   bmass -       The mass to which b-jets should be fixed.
0190     //
0191     // Outputs:
0192     //   fe -          The initialized Fourvec_Event.
0193     //
0194     {
0195       assert(ev.nleps() == 1);
0196       fe.add(make_fe_obj(ev.lep(0), 0, lepton_label));
0197
0198       bool saw_lepb = false;
0199       bool saw_hadb = false;
0200       for (std::vector<Lepjets_Event_Jet>::size_type j = 0; j < ev.njets(); j++) {
0201         if (ev.jet(j).type() == isr_label || ev.jet(j).type() == higgs_label)
0202           continue;
0203         double mass = 0;
0204         if (ev.jet(j).type() == lepb_label && !saw_lepb) {
0205           mass = bmass;
0206           saw_lepb = true;
0207         } else if (ev.jet(j).type() == hadb_label && !saw_hadb) {
0208           mass = bmass;
0209           saw_hadb = true;
0210         }
0211         fe.add(make_fe_obj(ev.jet(j), mass, ev.jet(j).type()));
0212       }
0213
0214       fe.set_nu_p(ev.met());
0215       Fourvec kt = ev.kt();
0216       fe.set_kt_error(ev.kt_res().sigma(kt.x()), ev.kt_res().sigma(kt.y()), 0);
0217     }
0218
0219     /**
0220
0221     @brief Convert from a Fourvec_Event to a Lepjets_Event.
0222
0223     @param fe The input event.
0224
0225     @param ev The returned Lepjets_Event.
0226
0227     @par Input:
0228     - Fourvec_Event <i>fe</i>.
0229     - Lepjets_Event <i>ev</i>.
0230
0231     @par Output:
0232     - Lepjets_Event <i>ev</i>
0233
0234  */
0235     void do_export(const Fourvec_Event& fe, Lepjets_Event& ev)
0236     //
0237     // Purpose: Convert from a Fourvec_Event to a Lepjets_Event.
0238     //
0239     // Inputs:
0240     //   fe -          The input event.
0241     //   ev -          The original Lepjets_Event.
0242     //
0243     // Outputs:
0244     //   ev -          The updated Lepjets_Event.
0245     //
0246     {
0247       ev.lep(0).p() = fe.obj(0).p;
0248       for (std::vector<Lepjets_Event_Jet>::size_type j = 0, k = 1; j < ev.njets(); j++) {
0249         if (ev.jet(j).type() == isr_label || ev.jet(j).type() == higgs_label)
0250           continue;
0251         ev.jet(j).p() = fe.obj(k++).p;
0252       }
0253
0254       ev.met() = fe.nu();
0255     }
0256
0257   }  // unnamed namespace
0258
0259   double Constrained_Top::constrain(
0260       Lepjets_Event& ev, double& mt, double& sigmt, Column_Vector& pullx, Column_Vector& pully)
0261   //
0262   // Purpose: Do a constrained fit for EV.  Returns the top mass and
0263   //          its error in MT and SIGMT, and the pull quantities in PULLX and
0264   //          PULLY.  Returns the chisq; this will be < 0 if the fit failed
0265   //          to converge.
0266   //
0267   // Inputs:
0268   //   ev -          The event we're fitting.
0269   //
0270   // Outputs:
0271   //   ev -          The fitted event.
0272   //   mt -          Requested invariant mass.
0273   //   sigmt -       Uncertainty on mt.
0274   //   pullx -       Pull quantities for well-measured variables.
0275   //   pully -       Pull quantities for poorly-measured variables.
0276   //
0277   // Returns:
0278   //   The fit chisq, or < 0 if the fit didn't converge.
0279   //
0280   {
0281     Fourvec_Event fe;
0282     do_import(ev, _args.bmass(), fe);
0283     double chisq = _constrainer.constrain(fe, mt, sigmt, pullx, pully);
0284     do_export(fe, ev);
0285
0286     return chisq;
0287   }
0288
0289   /**
0290     @brief Output stream operator, print the content of this Constrained_Top
0291     object to an output stream.
0292
0293     @param s The output stream to which to write.
0294
0295     @param ct The instance of Constrained_Top to be printed.
0296
0297 */
0298   std::ostream& operator<<(std::ostream& s, const Constrained_Top& ct)
0299   //
0300   // Purpose: Print the object to S.
0301   //
0302   // Inputs:
0303   //   s -           The stream to which to write.
0304   //   ct -          The object to write.
0305   //
0306   // Returns:
0307   //   The stream S.
0308   //
0309   {
0310     return s << ct._constrainer;
0311   }
0312
0313 }  // namespace hitfit