rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

MC_VH2BB

MC unboosted VH2bb validation plots
Experiment: (LHC)
Status: VALIDATED
Authors:
  • Ben Smart
  • Andy Buckley
No references listed
Beams: * *
Beam energies: ANY
Run details:
  • $VH$ with $H \to b \bar{b}$ and the vector boson decaying to electron or muon channels.

Various plots for characterising the process $V H \to b\bar{b}$

Source code: MC_VH2BB.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/ZFinder.hh"
  5#include "Rivet/Projections/WFinder.hh"
  6#include "Rivet/Projections/UnstableParticles.hh"
  7#include "Rivet/Projections/FastJets.hh"
  8#include "Rivet/Math/LorentzTrans.hh"
  9
 10namespace Rivet {
 11
 12
 13  class MC_VH2BB : public Analysis {
 14  public:
 15
 16    /// @name Constructors etc.
 17    //@{
 18
 19    /// Constructor
 20    MC_VH2BB()
 21      : Analysis("MC_VH2BB")
 22    {    }
 23
 24    //@}
 25
 26
 27    /// @name Analysis methods
 28    //@{
 29
 30    /// Book histograms and initialise projections before the run
 31    void init() {
 32
 33      FinalState fs;
 34      Cut cut = Cuts::abseta < 3.5 && Cuts::pT > 25*GeV;
 35      ZFinder zeefinder(fs, cut, PID::ELECTRON, 65*GeV, 115*GeV, 0.2);
 36      declare(zeefinder, "ZeeFinder");
 37      ZFinder zmmfinder(fs, cut, PID::MUON, 65*GeV, 115*GeV, 0.2);
 38      declare(zmmfinder, "ZmmFinder");
 39      WFinder wefinder(fs, cut, PID::ELECTRON, 60*GeV, 100*GeV, 25*GeV, 0.2);
 40      declare(wefinder, "WeFinder");
 41      WFinder wmfinder(fs, cut, PID::MUON, 60*GeV, 100*GeV, 25*GeV, 0.2);
 42      declare(wmfinder, "WmFinder");
 43
 44      declare(fs, "FinalState");
 45      declare(FastJets(fs, FastJets::ANTIKT, 0.4), "AntiKT04");
 46      declare(FastJets(fs, FastJets::ANTIKT, 0.5), "AntiKT05");
 47      declare(FastJets(fs, FastJets::ANTIKT, 0.6), "AntiKT06");
 48
 49      /// Book histograms
 50      book(_h_jet_bb_Delta_eta ,"jet_bb_Delta_eta", 50, 0, 4);
 51      book(_h_jet_bb_Delta_phi ,"jet_bb_Delta_phi", 50, 0, 1);
 52      book(_h_jet_bb_Delta_pT ,"jet_bb_Delta_pT", 50,0, 500);
 53      book(_h_jet_bb_Delta_R ,"jet_bb_Delta_R", 50, 0, 5);
 54      book(_h_jet_b_jet_eta ,"jet_b_jet_eta", 50, -4, 4);
 55      book(_h_jet_b_jet_multiplicity ,"jet_b_jet_multiplicity", 11, -0.5, 10.5);
 56      book(_h_jet_b_jet_phi ,"jet_b_jet_phi", 50, 0, 1);
 57      book(_h_jet_b_jet_pT ,"jet_b_jet_pT", 50, 0, 500);
 58      book(_h_jet_H_eta_using_bb ,"jet_H_eta_using_bb", 50, -4, 4);
 59      book(_h_jet_H_mass_using_bb ,"jet_H_mass_using_bb", 50, 50, 200);
 60      book(_h_jet_H_phi_using_bb ,"jet_H_phi_using_bb", 50, 0, 1);
 61      book(_h_jet_H_pT_using_bb ,"jet_H_pT_using_bb", 50, 0, 500);
 62      book(_h_jet_eta ,"jet_eta", 50, -4, 4);
 63      book(_h_jet_multiplicity ,"jet_multiplicity", 11, -0.5, 10.5);
 64      book(_h_jet_phi ,"jet_phi", 50, 0, 1);
 65      book(_h_jet_pT ,"jet_pT", 50, 0, 500);
 66      book(_h_jet_VBbb_Delta_eta ,"jet_VBbb_Delta_eta", 50, 0, 4);
 67      book(_h_jet_VBbb_Delta_phi ,"jet_VBbb_Delta_phi", 50, 0, 1);
 68      book(_h_jet_VBbb_Delta_pT ,"jet_VBbb_Delta_pT", 50, 0, 500);
 69      book(_h_jet_VBbb_Delta_R ,"jet_VBbb_Delta_R", 50, 0, 8);
 70
 71      book(_h_VB_eta ,"VB_eta", 50, -4, 4);
 72      book(_h_VB_mass ,"VB_mass", 50, 60, 110);
 73      book(_h_Z_multiplicity ,"Z_multiplicity", 11, -0.5, 10.5);
 74      book(_h_W_multiplicity ,"W_multiplicity", 11, -0.5, 10.5);
 75      book(_h_VB_phi ,"VB_phi", 50, 0, 1);
 76      book(_h_VB_pT ,"VB_pT", 50, 0, 500);
 77
 78      book(_h_jet_bVB_angle_Hframe ,"jet_bVB_angle_Hframe", 50, 0, 1);
 79      book(_h_jet_bb_angle_Hframe ,"jet_bb_angle_Hframe", 50, 0, 1);
 80      book(_h_jet_bVB_cosangle_Hframe ,"jet_bVB_cosangle_Hframe", 50, -1, 1);
 81      book(_h_jet_bb_cosangle_Hframe ,"jet_bb_cosangle_Hframe", 50, -1, 1);
 82    }
 83
 84
 85    /// Perform the per-event analysis
 86    void analyze(const Event& event) {
 87      const double JETPTCUT = 30*GeV;
 88
 89      const ZFinder& zeefinder = apply<ZFinder>(event, "ZeeFinder");
 90      const ZFinder& zmmfinder = apply<ZFinder>(event, "ZmmFinder");
 91      const WFinder& wefinder = apply<WFinder>(event, "WeFinder");
 92      const WFinder& wmfinder = apply<WFinder>(event, "WmFinder");
 93      const Particles vectorBosons = zeefinder.bosons() + zmmfinder.bosons() + wefinder.bosons() + wmfinder.bosons();
 94      _h_Z_multiplicity->fill(zeefinder.bosons().size() + zmmfinder.bosons().size());
 95      _h_W_multiplicity->fill(wefinder.bosons().size() + wmfinder.bosons().size());
 96
 97      const Jets jets = apply<FastJets>(event, "AntiKT04").jetsByPt(JETPTCUT);
 98      _h_jet_multiplicity->fill(jets.size());
 99
100      // Identify the b-jets
101      Jets bjets;
102      for (const Jet& jet : jets) {
103        const double jetEta = jet.eta();
104        const double jetPhi = jet.phi();
105        const double jetPt = jet.pT();
106        _h_jet_eta->fill(jetEta);
107        _h_jet_phi->fill(jetPhi/2/M_PI);
108        _h_jet_pT->fill(jetPt/GeV);
109
110        if (jet.bTagged() && jet.pT() > JETPTCUT) {
111          bjets.push_back(jet);
112          _h_jet_b_jet_eta->fill(jetEta);
113          _h_jet_b_jet_phi->fill(jetPhi/2/M_PI);
114          _h_jet_b_jet_pT->fill(jetPt);
115        }
116      }
117      _h_jet_b_jet_multiplicity->fill(bjets.size());
118
119      // Plot vector boson properties
120      for (const Particle& v : vectorBosons) {
121        _h_VB_phi->fill(v.phi()/2/M_PI);
122        _h_VB_pT->fill(v.pT());
123        _h_VB_eta->fill(v.eta());
124        _h_VB_mass->fill(v.mass());
125      }
126
127      // rest of analysis requires at least 1 b jets
128      if(bjets.empty()) vetoEvent;
129
130      // Construct Higgs candidates from pairs of b-jets
131      for (size_t i = 0; i < bjets.size()-1; ++i) {
132        for (size_t j = i+1; j < bjets.size(); ++j) {
133          const Jet& jet1 = bjets[i];
134          const Jet& jet2 = bjets[j];
135
136          const double deltaEtaJJ = fabs(jet1.eta() - jet2.eta());
137          const double deltaPhiJJ = deltaPhi(jet1.momentum(), jet2.momentum());
138          const double deltaRJJ = deltaR(jet1.momentum(), jet2.momentum());
139          const double deltaPtJJ = fabs(jet1.pT() - jet2.pT());
140          _h_jet_bb_Delta_eta->fill(deltaEtaJJ);
141          _h_jet_bb_Delta_phi->fill(deltaPhiJJ/M_PI);
142          _h_jet_bb_Delta_pT->fill(deltaPtJJ);
143          _h_jet_bb_Delta_R->fill(deltaRJJ);
144
145          const FourMomentum phiggs = jet1.momentum() + jet2.momentum();
146          _h_jet_H_eta_using_bb->fill(phiggs.eta());
147          _h_jet_H_mass_using_bb->fill(phiggs.mass());
148          _h_jet_H_phi_using_bb->fill(phiggs.phi()/2/M_PI);
149          _h_jet_H_pT_using_bb->fill(phiggs.pT());
150
151          for (const Particle& v : vectorBosons) {
152            const double deltaEtaVH = fabs(phiggs.eta() - v.eta());
153            const double deltaPhiVH = deltaPhi(phiggs, v.momentum());
154            const double deltaRVH = deltaR(phiggs, v.momentum());
155            const double deltaPtVH = fabs(phiggs.pT() - v.pT());
156            _h_jet_VBbb_Delta_eta->fill(deltaEtaVH);
157            _h_jet_VBbb_Delta_phi->fill(deltaPhiVH/M_PI);
158            _h_jet_VBbb_Delta_pT->fill(deltaPtVH);
159            _h_jet_VBbb_Delta_R->fill(deltaRVH);
160
161            // Calculate boost angles
162            const vector<double> angles = boostAngles(jet1.momentum(), jet2.momentum(), v.momentum());
163            _h_jet_bVB_angle_Hframe->fill(angles[0]/M_PI);
164            _h_jet_bb_angle_Hframe->fill(angles[1]/M_PI);
165            _h_jet_bVB_cosangle_Hframe->fill(cos(angles[0]));
166            _h_jet_bb_cosangle_Hframe->fill(cos(angles[1]));
167          }
168
169        }
170      }
171    }
172
173
174    /// Normalise histograms etc., after the run
175    void finalize() {
176      scale(_h_jet_bb_Delta_eta, crossSection()/sumOfWeights());
177      scale(_h_jet_bb_Delta_phi, crossSection()/sumOfWeights());
178      scale(_h_jet_bb_Delta_pT, crossSection()/sumOfWeights());
179      scale(_h_jet_bb_Delta_R, crossSection()/sumOfWeights());
180      scale(_h_jet_b_jet_eta, crossSection()/sumOfWeights());
181      scale(_h_jet_b_jet_multiplicity, crossSection()/sumOfWeights());
182      scale(_h_jet_b_jet_phi, crossSection()/sumOfWeights());
183      scale(_h_jet_b_jet_pT, crossSection()/sumOfWeights());
184      scale(_h_jet_H_eta_using_bb, crossSection()/sumOfWeights());
185      scale(_h_jet_H_mass_using_bb, crossSection()/sumOfWeights());
186      scale(_h_jet_H_phi_using_bb, crossSection()/sumOfWeights());
187      scale(_h_jet_H_pT_using_bb, crossSection()/sumOfWeights());
188      scale(_h_jet_eta, crossSection()/sumOfWeights());
189      scale(_h_jet_multiplicity, crossSection()/sumOfWeights());
190      scale(_h_jet_phi, crossSection()/sumOfWeights());
191      scale(_h_jet_pT, crossSection()/sumOfWeights());
192      scale(_h_jet_VBbb_Delta_eta, crossSection()/sumOfWeights());
193      scale(_h_jet_VBbb_Delta_phi, crossSection()/sumOfWeights());
194      scale(_h_jet_VBbb_Delta_pT, crossSection()/sumOfWeights());
195      scale(_h_jet_VBbb_Delta_R, crossSection()/sumOfWeights());
196
197      scale(_h_VB_eta, crossSection()/sumOfWeights());
198      scale(_h_VB_mass, crossSection()/sumOfWeights());
199      scale(_h_Z_multiplicity, crossSection()/sumOfWeights());
200      scale(_h_W_multiplicity, crossSection()/sumOfWeights());
201      scale(_h_VB_phi, crossSection()/sumOfWeights());
202      scale(_h_VB_pT, crossSection()/sumOfWeights());
203
204      scale(_h_jet_bVB_angle_Hframe, crossSection()/sumOfWeights());
205      scale(_h_jet_bb_angle_Hframe, crossSection()/sumOfWeights());
206      scale(_h_jet_bVB_cosangle_Hframe, crossSection()/sumOfWeights());
207      scale(_h_jet_bb_cosangle_Hframe, crossSection()/sumOfWeights());
208    }
209
210
211    /// This should take in the four-momenta of two b's (jets/hadrons) and a vector boson, for the process VB*->VBH with H->bb
212    /// It should return the smallest angle between the virtual vector boson and one of the b's, in the rest frame of the Higgs boson.
213    /// It should also return (as the second element of the vector) the angle between the b's, in the rest frame of the Higgs boson.
214    vector<double> boostAngles(const FourMomentum& b1, const FourMomentum& b2, const FourMomentum& vb) {
215      const FourMomentum higgsMomentum = b1 + b2;
216      const FourMomentum virtualVBMomentum = higgsMomentum + vb;
217      const LorentzTransform lt = LorentzTransform::mkFrameTransformFromBeta(higgsMomentum.betaVec());
218
219      const FourMomentum virtualVBMomentumBOOSTED = lt.transform(virtualVBMomentum);
220      const FourMomentum b1BOOSTED = lt.transform(b1);
221      const FourMomentum b2BOOSTED = lt.transform(b2);
222
223      const double angle1 = b1BOOSTED.angle(virtualVBMomentumBOOSTED);
224      const double angle2 = b2BOOSTED.angle(virtualVBMomentumBOOSTED);
225
226      const double anglebb = mapAngle0ToPi(b1BOOSTED.angle(b2BOOSTED));
227
228      vector<double> rtn;
229      rtn.push_back(angle1 < angle2 ? angle1 : angle2);
230      rtn.push_back(anglebb);
231      return rtn;
232    }
233
234    //@}
235
236
237  private:
238
239    /// @name Histograms
240    //@{
241
242    Histo1DPtr _h_Z_multiplicity, _h_W_multiplicity;
243    Histo1DPtr _h_jet_bb_Delta_eta, _h_jet_bb_Delta_phi, _h_jet_bb_Delta_pT, _h_jet_bb_Delta_R;
244    Histo1DPtr _h_jet_b_jet_eta, _h_jet_b_jet_multiplicity, _h_jet_b_jet_phi, _h_jet_b_jet_pT;
245    Histo1DPtr _h_jet_H_eta_using_bb, _h_jet_H_mass_using_bb, _h_jet_H_phi_using_bb, _h_jet_H_pT_using_bb;
246    Histo1DPtr _h_jet_eta, _h_jet_multiplicity, _h_jet_phi, _h_jet_pT;
247    Histo1DPtr _h_jet_VBbb_Delta_eta, _h_jet_VBbb_Delta_phi, _h_jet_VBbb_Delta_pT, _h_jet_VBbb_Delta_R;
248    Histo1DPtr _h_VB_eta, _h_VB_mass, _h_VB_phi, _h_VB_pT;
249    Histo1DPtr _h_jet_bVB_angle_Hframe, _h_jet_bb_angle_Hframe, _h_jet_bVB_cosangle_Hframe, _h_jet_bb_cosangle_Hframe;
250    //Histo1DPtr _h_jet_cuts_bb_deltaR_v_HpT;
251
252    //@}
253
254  };
255
256
257  // This global object acts as a hook for the plugin system
258  RIVET_DECLARE_PLUGIN(MC_VH2BB);
259
260}