1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief Higgs differential cross-section combination between the ATLAS measurements in the yy and 4l channels
 10  ///
 11  /// Computes Higgs transverse momentum, rapidity, jet multiplicity and leading jet pT.
 12  ///
 13  /// @author Michaela Queitsch-Maitland <michaela.queitsch-maitland@cern.ch>
 14  /// @author Dag Gillberg <dag.gillberg@cern.ch>
 15  /// @author Florian Bernlochner <florian.bernlochner@cern.ch>
 16  /// @author Sarah Heim <sarah.heim@cern.ch>
 17  class ATLAS_2015_I1364361 : public Analysis {
 18  public:
 19
 20    /// Constructor
 21    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2015_I1364361);
 22
 23
 24    /// Book histograms and initialise projections before the run
 25    void init() {
 26
 27      // All final state particles
 28      const FinalState fs;
 29      declare(fs, "FS");
 30
 31      // Histograms with data bins
 32      book(_h_pTH_incl,   1, 1, 1);
 33      book(_h_yH_incl,    2, 1, 1);
 34      book(_h_Njets_incl, 3, 1, 1);
 35      book(_h_pTj1_incl,  4, 1, 1);
 36    }
 37
 38
 39    /// Perform the per-event analysis
 40    void analyze(const Event& event) {
 41
 42      // Get the final state particles ordered by pT
 43      const Particles& fs = apply<FinalState>(event, "FS").particlesByPt();
 44
 45      // Find a stable Higgs (mandatory)
 46      const auto higgsiter = std::find_if(fs.begin(), fs.end(), [](const Particle& p){ return p.pid() == PID::HIGGSBOSON; });
 47      if (higgsiter == fs.end()) vetoEvent;
 48      const Particle& higgs = *higgsiter;
 49
 50      // bool stable_higgs = false;
 51      // const Particle* higgs = 0;
 52      // for (const Particle& p : fs) {
 53      //   if (p.pid() == PID::HIGGSBOSON) {
 54      //     stable_higgs = true;
 55      //     higgs = &p;
 56      //     break;
 57      //   }
 58      // }
 59      // if (!stable_higgs) {
 60      //   MSG_WARNING("FATAL: No stable Higgs found in event record.\n");
 61      //   vetoEvent;
 62      // }
 63
 64
 65      // Loop over final state particles and fill various particle vectors
 66      Particles leptons, photons, jet_ptcls;
 67      for ( const Particle& ptcl : fs ) {
 68        // Do not include the Higgs in jet finding!
 69        if ( ptcl.pid() == PID::HIGGSBOSON ) continue;
 70        // Neutrinos not from hadronisation
 71        if ( ptcl.isNeutrino() && !ptcl.fromHadron() ) continue;
 72        // Electrons and muons not from hadronisation
 73        if ( ( ptcl.abspid() == PID::ELECTRON || ptcl.abspid() == PID::MUON ) && !ptcl.fromHadron() ) {
 74          leptons.push_back(ptcl);
 75          continue;
 76        }
 77        // Photons not from hadronisation
 78        if ( ptcl.abspid() == PID::PHOTON && !ptcl.fromHadron() ) {
 79          photons.push_back(ptcl);
 80          continue;
 81        }
 82        // Add particle to jet inputs
 83        jet_ptcls.push_back(ptcl);
 84      }
 85
 86      // Match FS photons to leptons within cone R=0.1
 87      // If they are not 'dressing' photons, add to jet particle vector
 88      for ( const Particle& ph : photons ) {
 89        bool fsr_photon = false;
 90        for ( const Particle& lep : leptons ) {
 91          if ( deltaR(ph.momentum(),lep.momentum()) < 0.1 ){
 92            fsr_photon=true;
 93            continue;
 94          }
 95        }
 96        if ( !fsr_photon ) jet_ptcls.push_back(ph);
 97      }
 98
 99      // Let's build the jets! By hand...
100      const PseudoJets pjs_in = mkPseudoJets(jet_ptcls);
101      const fastjet::JetDefinition jdef(fastjet::antikt_algorithm, 0.4);
102      const Jets alljets = mkJets(fastjet::ClusterSequence(pjs_in, jdef).inclusive_jets());
103      const Jets jets = sortByPt(select(alljets, Cuts::pT > 30*GeV && Cuts::absrap < 4.4));
104      // FastJets jet_pro(JetAlg::ANTIKT, 0.4);
105      // jet_pro.calc(jet_ptcls);
106      // Jets jets = jet_pro.jetsByPt(Cuts::pT > 30*GeV && Cuts::absrap < 4.4);
107
108      size_t njets = jets.size() > 3 ? 3 : jets.size();
109      _h_pTH_incl->fill(higgs.pT());
110      _h_yH_incl->fill(higgs.absrap());
111      const string edge = _h_Njets_incl->xEdges()[njets];
112      _h_Njets_incl->fill(edge);
113      _h_pTj1_incl->fill(jets.empty() ? 0 : jets[0].pT());
114    }
115
116
117    /// Normalise histograms etc., after the run
118    void finalize() {
119      const double xs = crossSectionPerEvent();
120      scale(_h_pTH_incl, xs);
121      scale(_h_yH_incl, xs);
122      scale(_h_Njets_incl, xs);
123      scale(_h_pTj1_incl, xs);
124    }
125
126
127  private:
128
129    Histo1DPtr _h_pTH_incl, _h_yH_incl, _h_pTj1_incl;
130    BinnedHistoPtr<string> _h_Njets_incl;
131
132  };
133
134
135  RIVET_DECLARE_PLUGIN(ATLAS_2015_I1364361);
136
137
138}