1// -*- C++ -*-
  2#include "Rivet/Analyses/MC_JETS_BASE.hh"
  3#include "Rivet/Projections/DileptonFinder.hh"
  4#include "Rivet/Projections/VetoedFinalState.hh"
  5#include "Rivet/Projections/IdentifiedFinalState.hh"
  6#include "Rivet/Projections/FastJets.hh"
  7
  8namespace Rivet {
  9
 10  
 11  /// @brief MC validation analysis for Higgs-pair events (stable Higgses)
 12  class MC_HHJETS : public MC_JETS_BASE {
 13  public:
 14
 15    /// Default constructor
 16    MC_HHJETS()
 17      : MC_JETS_BASE("MC_HHJETS", 4, "Jets")
 18    {    }
 19
 20
 21    /// @name Analysis methods
 22    /// @{
 23
 24    /// Book histograms
 25    void init() {
 26      IdentifiedFinalState ifs(Cuts::abseta < 10.0 && Cuts::pT > 0*GeV);
 27      ifs.acceptId(25);
 28      declare(ifs,"IFS");
 29
 30      VetoedFinalState vfs;
 31      vfs.addVetoPairId(25);
 32
 33      // set ptcut from input option
 34      const double jetptcut = getOption<double>("PTJMIN", 20.0);
 35      _jetptcut = jetptcut * GeV;
 36
 37      // set clustering radius from input option
 38      const double R = getOption<double>("R", 0.4);
 39
 40      // set clustering algorithm from input option
 41      JetAlg clusterAlgo;
 42      const string algoopt = getOption("ALGO", "ANTIKT");
 43      if ( algoopt == "KT" ) {
 44	clusterAlgo = JetAlg::KT;
 45      } else if ( algoopt == "CA" ) {
 46	clusterAlgo = JetAlg::CA;
 47      } else if ( algoopt == "ANTIKT" ) {
 48	clusterAlgo = JetAlg::ANTIKT;
 49      } else {
 50	MSG_WARNING("Unknown jet clustering algorithm option " + algoopt + ". Defaulting to anti-kT");
 51	clusterAlgo = JetAlg::ANTIKT;
 52      }
 53
 54      FastJets jetpro(vfs, clusterAlgo, R);
 55      declare(jetpro, "Jets");
 56
 57      book(_h_HH_mass ,"HH_mass", 250, 240, 4000.0);
 58      book(_h_HH_dR ,"HH_dR", 25, 0.5, 10.0);
 59      book(_h_HH_dPhi ,"HH_dPhi", 64, 0, 3.2);
 60      book(_h_HH_deta,"HH_deta", 50, -5, 5);
 61      book(_h_H_pT ,"H_pT", 50, 0, 2000.0);
 62      book(_h_HH_pT ,"HH_pT", 200, 0, 2000.0);
 63      book(_h_H_pT1 ,"H_pT1", 200, 0, 2000.0);
 64      book(_h_H_pT2 ,"H_pT2", 200, 0, 2000.0);
 65      book(_h_H_eta ,"H_eta", 50, -5.0, 5.0);
 66      book(_h_H_eta1 ,"H_eta1", 50, -5.0, 5.0);
 67      book(_h_H_eta2 ,"H_eta2", 50, -5.0, 5.0);
 68      book(_h_H_phi ,"H_phi", 25, 0.0, TWOPI);
 69      book(_h_H_jet1_deta ,"H_jet1_deta", 50, -5.0, 5.0);
 70      book(_h_H_jet1_dR ,"H_jet1_dR", 25, 0.5, 7.0);
 71
 72      MC_JETS_BASE::init();
 73    }
 74
 75
 76
 77    /// Do the analysis
 78    void analyze(const Event & e) {
 79
 80      const IdentifiedFinalState& ifs = apply<IdentifiedFinalState>(e, "IFS");
 81      Particles allp = ifs.particlesByPt();
 82      if (allp.empty()) vetoEvent;
 83
 84
 85      FourMomentum hmom = allp[0].momentum();
 86      if (allp.size() > 1) {
 87        FourMomentum hmom2(allp[1].momentum());
 88        _h_HH_dR->fill(deltaR(hmom, hmom2));
 89        _h_HH_dPhi->fill(deltaPhi(hmom, hmom2));
 90        _h_HH_deta->fill(hmom.eta()-hmom2.eta());
 91        _h_HH_pT->fill((hmom+hmom2).pT());
 92        _h_HH_mass->fill((hmom+hmom2).mass());
 93
 94        if (hmom.pT() > hmom2.pT()) {
 95          _h_H_pT1->fill(hmom.pT());
 96          _h_H_eta1->fill(hmom.eta());
 97          _h_H_pT2->fill(hmom2.pT());
 98          _h_H_eta2->fill(hmom2.eta());
 99        } else {
100          _h_H_pT1->fill(hmom2.pT());
101          _h_H_eta1->fill(hmom2.eta());
102          _h_H_pT2->fill(hmom.pT());
103          _h_H_eta2->fill(hmom.eta());
104        }
105      }
106      _h_H_pT->fill(hmom.pT());
107      _h_H_eta->fill(hmom.eta());
108      _h_H_phi->fill(hmom.azimuthalAngle());
109
110
111      // Get the jet candidates
112      Jets jets = apply<FastJets>(e, "Jets").jetsByPt(Cuts::pT > 20*GeV);
113      if (!jets.empty()) {
114        _h_H_jet1_deta->fill(deltaEta(hmom, jets[0]));
115        _h_H_jet1_dR->fill(deltaR(hmom, jets[0]));
116      }
117
118      MC_JETS_BASE::analyze(e);
119    }
120
121
122    /// Finalize
123    void finalize() {
124      scale(_h_HH_mass, crossSection()/picobarn/sumOfWeights());
125      scale(_h_HH_dR, crossSection()/picobarn/sumOfWeights());
126      scale(_h_HH_deta, crossSection()/picobarn/sumOfWeights());
127      scale(_h_HH_dPhi, crossSection()/picobarn/sumOfWeights());
128      scale(_h_H_pT, crossSection()/picobarn/sumOfWeights());
129      scale(_h_H_pT1, crossSection()/picobarn/sumOfWeights());
130      scale(_h_H_pT2, crossSection()/picobarn/sumOfWeights());
131      scale(_h_HH_pT, crossSection()/picobarn/sumOfWeights());
132      scale(_h_H_eta, crossSection()/picobarn/sumOfWeights());
133      scale(_h_H_eta1, crossSection()/picobarn/sumOfWeights());
134      scale(_h_H_eta2, crossSection()/picobarn/sumOfWeights());
135      scale(_h_H_phi, crossSection()/picobarn/sumOfWeights());
136      scale(_h_H_jet1_deta, crossSection()/picobarn/sumOfWeights());
137      scale(_h_H_jet1_dR, crossSection()/picobarn/sumOfWeights());
138
139      MC_JETS_BASE::finalize();
140    }
141
142    /// @}
143
144
145  private:
146
147    /// @name Histograms
148    /// @{
149    Histo1DPtr _h_HH_mass, _h_HH_pT, _h_HH_dR, _h_HH_deta, _h_HH_dPhi;
150    Histo1DPtr _h_H_pT, _h_H_pT1, _h_H_pT2, _h_H_eta, _h_H_eta1, _h_H_eta2, _h_H_phi;
151    Histo1DPtr _h_H_jet1_deta, _h_H_jet1_dR;
152    /// @}
153
154  };
155
156
157  RIVET_DECLARE_PLUGIN(MC_HHJETS);
158
159}