1#include "Rivet/Analysis.hh"
  2#include "Rivet/Projections/FinalState.hh"
  3#include "Rivet/Projections/ChargedFinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// g -> bb at 13 TeV
 10  class ATLAS_2018_I1711114: public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2018_I1711114);
 15
 16
 17    /// Book cuts and projections
 18    void init() {
 19      // All final state particles
 20      FinalState fs(Cuts::abseta<5.0);
 21
 22      ChargedFinalState cfs(Cuts::pT > 0.5*GeV && Cuts::abseta < 2.5);
 23      FastJets smallR_jets(cfs, JetAlg::ANTIKT, 0.2, JetMuons::NONE, JetInvisibles::NONE);
 24      declare(smallR_jets, "track_jets");
 25
 26      FastJets largeR_jets(fs, JetAlg::ANTIKT, 1.0, JetMuons::NONE, JetInvisibles::NONE);
 27      declare(largeR_jets, "largeR_jets");
 28
 29      book(_h_R, 1,1,1);
 30      book(_h_phi, 2,1,1);
 31      book(_h_z, 3,1,1);
 32      book(_h_rho, 4,1,1);
 33    }
 34
 35
 36    void analyze(const Event& event) {
 37
 38      const PseudoJets& myJets = apply<FastJets>(event, "largeR_jets").pseudojetsByPt(450*GeV);
 39      if (myJets.empty()) vetoEvent;
 40
 41      fastjet::Filter trimmer(fastjet::JetDefinition(fastjet::kt_algorithm, 0.2), fastjet::SelectorPtFractionMin(0.05));
 42      Jets myTrimmedJets;
 43      for (const Jet jet : myJets)  myTrimmedJets += Jet(trimmer(jet));
 44      std::sort(myTrimmedJets.begin(), myTrimmedJets.end(), cmpMomByPt);
 45      if (myTrimmedJets[0].pT() < 450*GeV) vetoEvent;
 46
 47      const Jets& myJets_charged = apply<FastJets>(event, "track_jets").jetsByPt(Cuts::pT > 10*GeV);
 48      PseudoJets pjs;
 49      for (auto tj : myTrimmedJets[0].pseudojet().constituents()) {
 50        fastjet::PseudoJet pj = tj;
 51        pj.set_user_index(-1); // dummmy
 52        pjs.push_back(pj);
 53      }
 54      for (size_t i = 0; i < myJets_charged.size(); ++i) {
 55        fastjet::PseudoJet pj = myJets_charged[i];
 56        pj *= 1e-20; // ghostify momentum
 57        pj.set_user_index(i);
 58        pjs.push_back(pj);
 59      }
 60      fastjet::ClusterSequence tagged_seq(pjs, fastjet::JetDefinition(fastjet::antikt_algorithm, 1.0));
 61      PseudoJets GAfatjet = fastjet::sorted_by_pt(tagged_seq.inclusive_jets(50.0));
 62      Jets associated_trackjets;
 63      for (auto pj : GAfatjet[0].constituents()) {
 64        if (pj.user_index() >= 0)  associated_trackjets += myJets_charged[pj.user_index()];
 65      }
 66      if (associated_trackjets.size() < 2) vetoEvent;
 67      std::sort(associated_trackjets.begin(), associated_trackjets.end(), cmpMomByPt);
 68
 69      size_t nbtags = 0;
 70      for (unsigned int ij = 0; ij < 2; ++ij) {
 71        if (associated_trackjets[ij].bTagged(Cuts::pT > 5*GeV))  ++nbtags;
 72      }
 73      if (nbtags != 2) vetoEvent;
 74
 75      // Now, time to compute the observables!
 76      const Vector3 b1v = associated_trackjets[0].p3();
 77      const Vector3 b2v = associated_trackjets[1].p3();
 78      const Vector3 plane1 = b1v.cross(b2v);
 79      const Vector3 gv = b1v + b2v;
 80      const Vector3 beam = Vector3(0,0,1);
 81      const Vector3 plane2 = beam.cross(gv);
 82
 83      const double fTz = associated_trackjets[1].pT() / (associated_trackjets[0].pT()+associated_trackjets[1].pT());
 84      const double fTR = deltaR(associated_trackjets[0], associated_trackjets[1], RapScheme::YRAP); //N.B. this uses y and not eta
 85      const FourMomentum dijet = associated_trackjets[0].mom() + associated_trackjets[1].mom();
 86      const double fTrho = log(dijet.mass() / dijet.pT());
 87      const double fTphi = (plane1).angle(plane2)/M_PI;
 88
 89      _h_R->fill(  fTR);
 90      _h_phi->fill(fTphi);
 91      _h_z->fill(  fTz);
 92      _h_rho->fill(fTrho);
 93    }
 94
 95
 96    /// Scale histos
 97    void finalize() {
 98      normalize(_h_R  ,  1.0, false);
 99      normalize(_h_phi,  1.0, false);
100      normalize(_h_z  ,  1.0, false);
101      normalize(_h_rho,  1.0, false);
102    }
103
104
105  private:
106
107    /// Histograms
108    Histo1DPtr _h_z, _h_R, _h_rho, _h_phi;
109
110  };
111
112
113  RIVET_DECLARE_PLUGIN(ATLAS_2018_I1711114);
114
115}