Rivet analyses referenceLHCB_2017_I1509507J/$\psi$ production in jets at 13 TeVExperiment: LHCB (LHC) Inspire ID: 1509507 Status: VALIDATED NOHEPDATA Authors:
Beam energies: (6500.0, 6500.0) GeV Run details:
Measurement of the fragmentation function for the production of J/$\psi$ in jets at 13 TeV by LHCb. Source code: LHCB_2017_I1509507.cc 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
8
9namespace Rivet {
10
11
12 /// @brief J/psi production in jets at 13 TeV
13 class LHCB_2017_I1509507 : public Analysis {
14 public:
15
16 /// Constructor
17 RIVET_DEFAULT_ANALYSIS_CTOR(LHCB_2017_I1509507);
18
19
20 /// @name Analysis methods
21 /// @{
22
23 /// Book histograms and initialise projections before the run
24 void init() {
25 declare(UnstableParticles(), "UFS");
26 declare(FinalState(), "FS");
27 for(unsigned int ix=0;ix<2;++ix)
28 book(_h_frag[ix],1,1,ix+1);
29 }
30
31 /// Perform the per-event analysis
32 void analyze(const Event& event) {
33 // first see if we have any J/psi in the region
34 Particles Jpsi;
35 for(const Particle & p : apply<UnstableParticles>(event, "UFS").particles(Cuts::pid==443 &&
36 Cuts::etaIn(2.,4.5))) {
37 if(p.children().size()!=2) continue;
38 bool found=true;
39 for(const Particle & child :p.children()) {
40 double eta = child.eta();
41 if (child.abspid()!=13 || eta<2. || eta>4.5 || child.perp()<0.5 || child.p3().mod()<5.) {
42 found=false;
43 break;
44 }
45 }
46 if (found) Jpsi.push_back(p);
47 }
48 // no jpsi veto
49 if(Jpsi.empty()) vetoEvent;
50 // now get the final-state particles for the jets
51 const Particles & fsParticles = apply<FinalState>(event, "FS").particles();
52 vector<PseudoJet> particles;
53 particles.reserve(fsParticles.size());
54 // fs for fastjet omitting any J/psi decay products
55 for(const Particle & p : fsParticles) {
56 // skip muons and neutrinos (seems standard for LHCb)
57 if(p.abspid()==13 || p.abspid()==12 || p.abspid()==14 or p.abspid()==16) continue;
58 // skip anything coming from the decay of one of the jpsis
59 Particle parent = p;
60 while(!parent.parents().empty()) {
61 if(parent.pid()==443) break;
62 parent=parent.parents()[0];
63 }
64 bool match = parent.pid()==443;
65 if(match) {
66 match =false;
67 for(const Particle & psi : Jpsi) {
68 match = fuzzyEquals(parent.momentum(),psi.momentum());
69 if(match) break;
70 }
71 }
72 if(!match) {
73 PseudoJet j = p.pseudojet();
74 j.set_user_index(0);
75 particles.push_back(j);
76 }
77 }
78 // add the jpsis to the particles for fastjet
79 for(const Particle & p : Jpsi) {
80 PseudoJet j = p.pseudojet();
81 j.set_user_index(p.fromBottom()+1);
82 particles.push_back(j);
83 }
84 JetDefinition jet_def(fastjet::antikt_algorithm, 0.5);
85 fastjet::ClusterSequence clu = ClusterSequence(particles,jet_def);
86 vector<PseudoJet> jets = clu.inclusive_jets();
87 for(const PseudoJet & jet : jets) {
88 // pt and eta cut
89 if(jet.perp()<20. || jet.eta()<2. || jet.eta()>4.) continue;
90 // loop over constituents and find jpsi
91 for(const PseudoJet & sub : jet.constituents()) {
92 if(sub.user_index()==0) continue;
93 double z = sub.perp()/jet.perp();
94 _h_frag[sub.user_index()-1]->fill(z);
95 }
96 }
97 }
98
99
100 /// Normalise histograms etc., after the run
101 void finalize() {
102 for(unsigned int ix=0;ix<2;++ix)
103 normalize(_h_frag[ix],1.,false);
104 }
105
106 /// @}
107
108
109 /// @name Histograms
110 /// @{
111 Histo1DPtr _h_frag[2];
112 /// @}
113
114
115 };
116
117
118 RIVET_DECLARE_PLUGIN(LHCB_2017_I1509507);
119
120}
|