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