1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/UnstableParticles.hh"
  5#include "fastjet/JetDefinition.hh"
  6#include "fastjet/ClusterSequence.hh"
  7
  8namespace Rivet {
  9
 10
 11  /// @brief J/psi in jets at 5.02 TeV
 12  class CMS_2022_I1870319 : public Analysis {
 13  public:
 14
 15    /// Constructor
 16    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2022_I1870319);
 17
 18
 19    /// @name Analysis methods
 20    /// @{
 21
 22    /// Book histograms and initialise projections before the run
 23    void init() {
 24      declare(UnstableParticles(), "UFS");
 25      declare(FinalState(), "FS");
 26      for (unsigned int ix=0; ix<2; ++ix) {
 27        book(_h_frag[ix],ix+1,1,1);
 28      }
 29    }
 30
 31
 32    /// Perform the per-event analysis
 33    void analyze(const Event& event) {
 34      // first see if we have any prompt J/psi in the region
 35      Particles Jpsi;
 36      for (const Particle& p :apply<UnstableParticles>(event, "UFS").particles(Cuts::pid==443 && Cuts::pT>6.5)) {
 37        if (!p.fromBottom()) Jpsi.push_back(p);
 38      }
 39      // no jpsi veto
 40      if (Jpsi.empty()) vetoEvent;
 41      // now get the final-state particles for the jets
 42      const Particles& fsParticles = apply<FinalState>(event, "FS").particles();
 43      vector<PseudoJet> particles;
 44      particles.reserve(fsParticles.size());
 45      // fs for fastjet omitting any J/psi decay products
 46      for (const Particle& p : fsParticles) {
 47        if (p.abspid()==12 || p.abspid()==14 or p.abspid()==16) continue;
 48        // skip anything coming from the decay of one of the jpsis
 49        Particle parent = p;
 50        while (!parent.parents().empty()) {
 51          if (parent.pid()==443) break;
 52          parent=parent.parents()[0];
 53        }
 54        bool match = parent.pid()==443;
 55        if (match) {
 56          match =false;
 57          for (const Particle& psi : Jpsi) {
 58            match = fuzzyEquals(parent.momentum(),psi.momentum());
 59            if (match) break;
 60          }
 61        }
 62        if (!match) {
 63          PseudoJet j = p.pseudojet();
 64          j.set_user_index(0);
 65          particles.push_back(j);
 66        }
 67      }
 68      // add the jpsis to the particles for fastjet
 69      for (const Particle& p : Jpsi) {
 70        PseudoJet j = p.pseudojet();
 71        j.set_user_index(1);
 72        particles.push_back(j);
 73      }
 74      JetDefinition jet_def(fastjet::antikt_algorithm, 0.3);
 75      fastjet::ClusterSequence clu = ClusterSequence(particles,jet_def);
 76      vector<PseudoJet> jets = clu.inclusive_jets();
 77      for (const PseudoJet& jet : jets) {
 78        // pt and eta cut
 79        if (jet.perp()<30. || jet.perp()>40. || abs(jet.eta())>2.) continue;
 80        // loop over constituents and find jpsi
 81        for (const PseudoJet& sub : jet.constituents()) {
 82          if (sub.user_index()==0) continue;
 83          const double z = sub.perp()/jet.perp();
 84          _h_frag[0]->fill(z);
 85          _h_frag[1]->fill(z);
 86        }
 87      }
 88    }
 89
 90
 91    /// Normalise histograms etc., after the run
 92    void finalize() {
 93      normalize(_h_frag[0]);
 94      scale(_h_frag[1], crossSection()/nanobarn/sumOfWeights());
 95    }
 96
 97    /// @}
 98
 99
100    /// @name Histograms
101    /// @{
102    Histo1DPtr _h_frag[2];
103    /// @}
104
105
106  };
107
108
109  RIVET_DECLARE_PLUGIN(CMS_2022_I1870319);
110
111}