rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

ATLAS_2012_I1094564

Jet mass and substructure of inclusive jets at 7 TeV
Experiment: ATLAS (LHC 7TeV)
Inspire ID: 1094564
Status: VALIDATED
Authors:
  • Karl Nordstrom
References: Beams: p+ p+
Beam energies: (3500.0, 3500.0) GeV
Run details:
  • pp QCD events at 7 TeV

In this analysis, the assumption that the internal substructure of jets generated by QCD radiation is well understood is tested on an inclusive sample of jets recorded with the ATLAS detector in 2010, which corresponds to 35 $pb^{-1}$ of $pp$ collisions delivered by the LHC at $\sqrt{s} = 7$ TeV. Jet invariant mass, $k_{t}$ splitting scales and $N$-subjettiness variables are presented for anti-$k_{t}$ $R = 1.0$ jets and Cambridge-Aachen $R = 1.2$ jets. Jet invariant-mass spectra for Cambridge-Aachen $R = 1.2$ jets after a splitting and filtering procedure are also presented.

Source code: ATLAS_2012_I1094564.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5#include "Rivet/Tools/BinnedHistogram.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// ATLAS jet substructure measurement
 11  class ATLAS_2012_I1094564 : public Analysis {
 12  public:
 13
 14    ATLAS_2012_I1094564()
 15      : Analysis("ATLAS_2012_I1094564")
 16    {}
 17
 18
 19    // Returns constituents to make it easier to do the filtering
 20    PseudoJets splitjet(fastjet::PseudoJet jet, double& last_R, const FastJets& fj, bool& unclustered) const {
 21
 22      // Build a new cluster sequence just using the constituents of this jet.
 23      fastjet::ClusterSequence cs(fj.clusterSeq()->constituents(jet), fastjet::JetDefinition(fastjet::cambridge_algorithm, M_PI/2.));
 24
 25      // Get the jet back again
 26      vector<fastjet::PseudoJet> remadeJets = cs.inclusive_jets(0.);
 27
 28      if ( remadeJets.size() != 1 ) return remadeJets;
 29
 30      fastjet::PseudoJet remadeJet = remadeJets[0];
 31      fastjet::PseudoJet parent1, parent2;
 32      unclustered = false;
 33
 34      while ( cs.has_parents(remadeJet, parent1, parent2) ) {
 35        if (parent1.squared_distance(parent2) < 0.09) break;
 36        if (parent1.m2() < parent2.m2()) {
 37          fastjet::PseudoJet tmp;
 38          tmp = parent1; parent1 = parent2; parent2 = tmp;
 39        }
 40
 41        double ktdist = parent1.kt_distance(parent2);
 42        double rtycut2 = 0.3*0.3;
 43        if (parent1.m() < 0.67*remadeJet.m() && ktdist > rtycut2*remadeJet.m2()) {
 44          unclustered = true;
 45          break;
 46        } else {
 47          remadeJet = parent1;
 48        }
 49      }
 50
 51      last_R = 0.5 * sqrt(parent1.squared_distance(parent2));
 52      return cs.constituents(remadeJet);
 53    }
 54
 55
 56    fastjet::PseudoJet filterjet(PseudoJets jets, double& stingy_R, const double def_R) const {
 57      if (stingy_R == 0.0) stingy_R = def_R;
 58      stingy_R = def_R < stingy_R ? def_R : stingy_R;
 59      fastjet::JetDefinition stingy_jet_def(fastjet::cambridge_algorithm, stingy_R);
 60      fastjet::ClusterSequence scs(jets, stingy_jet_def);
 61      vector<fastjet::PseudoJet> stingy_jets = sorted_by_pt(scs.inclusive_jets(0));
 62      fastjet::PseudoJet reconst_jet(0, 0, 0, 0);
 63      for (size_t isj = 0; isj < std::min((size_t) 3, stingy_jets.size()); ++isj) {
 64        reconst_jet += stingy_jets[isj];
 65      }
 66      return reconst_jet;
 67    }
 68
 69
 70    // These are custom functions for n-subjettiness.
 71    PseudoJets jetGetAxes(int n_jets, const PseudoJets& inputJets, double subR) const {
 72      // Sanity check
 73      if (inputJets.size() < (size_t) n_jets) {
 74        MSG_ERROR("Not enough input particles.");
 75        return inputJets;
 76      }
 77
 78      // Get subjets, return
 79      fastjet::ClusterSequence sub_clust_seq(inputJets, fastjet::JetDefinition(fastjet::kt_algorithm, subR, fastjet::E_scheme, fastjet::Best));
 80      return sub_clust_seq.exclusive_jets(n_jets);
 81    }
 82
 83
 84    double jetTauValue(double beta, double jet_rad, const PseudoJets& particles, const PseudoJets& axes, double Rcut) const {
 85      double tauNum = 0.0;
 86      double tauDen = 0.0;
 87
 88      if (particles.size() == 0) return 0.0;
 89
 90      for (size_t i = 0; i < particles.size(); i++) {
 91        // find minimum distance (set R large to begin)
 92        double minR = 10000.0;
 93        for (size_t j = 0; j < axes.size(); j++) {
 94          double tempR = sqrt(particles[i].squared_distance(axes[j]));
 95          if (tempR < minR) minR = tempR;
 96        }
 97        if (minR > Rcut) minR = Rcut;
 98        // calculate nominator and denominator
 99        tauNum += particles[i].perp() * pow(minR,beta);
100        tauDen += particles[i].perp() * pow(jet_rad,beta);
101      }
102
103      // return N-subjettiness (or 0 if denominator is 0)
104      return safediv(tauNum, tauDen, 0);
105    }
106
107
108    void jetUpdateAxes(double beta, const PseudoJets& particles, PseudoJets& axes) const {
109      vector<int> belongsto;
110      for (size_t i = 0; i < particles.size(); i++) {
111        // find minimum distance axis
112        int assign = 0;
113        double minR = 10000.0;
114        for (size_t j = 0; j < axes.size(); j++) {
115          double tempR = sqrt(particles[i].squared_distance(axes[j]));
116          if (tempR < minR) {
117            minR = tempR;
118            assign = j;
119          }
120        }
121        belongsto.push_back(assign);
122      }
123
124      // iterative step
125      double deltaR2, distphi;
126      vector<double> ynom, phinom, den;
127
128      ynom.resize(axes.size());
129      phinom.resize(axes.size());
130      den.resize(axes.size());
131
132      for (size_t i = 0; i < particles.size(); i++) {
133        distphi = particles[i].phi() - axes[belongsto[i]].phi();
134        deltaR2 = particles[i].squared_distance(axes[belongsto[i]]);
135        if (deltaR2 == 0.) continue;
136        if (abs(distphi) <= M_PI) phinom.at(belongsto[i]) += particles[i].perp() * particles[i].phi() * pow(deltaR2, (beta-2)/2);
137        else if ( distphi > M_PI) phinom.at(belongsto[i]) += particles[i].perp() * (-2 * M_PI + particles[i].phi()) * pow(deltaR2, (beta-2)/2);
138        else if ( distphi < M_PI) phinom.at(belongsto[i]) += particles[i].perp() * (+2 * M_PI + particles[i].phi()) * pow(deltaR2, (beta-2)/2);
139        ynom.at(belongsto[i]) += particles[i].perp() * particles[i].rap() * pow(deltaR2, (beta-2)/2);
140        den.at(belongsto[i])  += particles[i].perp() * pow(deltaR2, (beta-2)/2);
141      }
142
143      // reset to new axes
144      for (size_t j = 0; j < axes.size(); j++) {
145        if (den[j] == 0.) axes.at(j) = axes[j];
146        else {
147          double phi_new = fmod( 2*M_PI + (phinom[j] / den[j]), 2*M_PI );
148          double pt_new  = axes[j].perp();
149          double y_new   = ynom[j] / den[j];
150          double px = pt_new * cos(phi_new);
151          double py = pt_new * sin(phi_new);
152          double pz = pt_new * sinh(y_new);
153          axes.at(j).reset(px, py, pz, axes[j].perp()/2);
154        }
155      }
156    }
157
158
159    void init() {
160
161      /// Projections:
162      FinalState fs((Cuts::etaIn(-4.5, 4.5) && Cuts::pT >=  0.*GeV));
163      declare(fs, "FS");
164      declare(FastJets(fs, FastJets::ANTIKT, 1.0), "AKT");
165      declare(FastJets(fs, FastJets::CAM, 1.2)   , "CA" );
166
167      /// Histograms:
168      {Histo1DPtr tmp; _h_camass.add(200, 300, book(tmp, 1, 1, 1));}
169      {Histo1DPtr tmp; _h_camass.add(300, 400, book(tmp, 2, 1, 1));}
170      {Histo1DPtr tmp; _h_camass.add(400, 500, book(tmp, 3, 1, 1));}
171      {Histo1DPtr tmp; _h_camass.add(500, 600, book(tmp, 4, 1, 1));}
172
173      {Histo1DPtr tmp; _h_filtmass.add(200, 300, book(tmp, 5, 1, 1));}
174      {Histo1DPtr tmp; _h_filtmass.add(300, 400, book(tmp, 6, 1, 1));}
175      {Histo1DPtr tmp; _h_filtmass.add(400, 500, book(tmp, 7, 1, 1));}
176      {Histo1DPtr tmp; _h_filtmass.add(500, 600, book(tmp, 8, 1, 1));}
177
178      {Histo1DPtr tmp; _h_ktmass.add(200, 300, book(tmp,  9, 1, 1));}
179      {Histo1DPtr tmp; _h_ktmass.add(300, 400, book(tmp, 10, 1, 1));}
180      {Histo1DPtr tmp; _h_ktmass.add(400, 500, book(tmp, 11, 1, 1));}
181      {Histo1DPtr tmp; _h_ktmass.add(500, 600, book(tmp, 12, 1, 1));}
182
183      {Histo1DPtr tmp; _h_ktd12.add(200, 300, book(tmp, 13, 1, 1));}
184      {Histo1DPtr tmp; _h_ktd12.add(300, 400, book(tmp, 14, 1, 1));}
185      {Histo1DPtr tmp; _h_ktd12.add(400, 500, book(tmp, 15, 1, 1));}
186      {Histo1DPtr tmp; _h_ktd12.add(500, 600, book(tmp, 16, 1, 1));}
187
188      {Histo1DPtr tmp; _h_ktd23.add(200, 300, book(tmp, 17, 1 ,1));}
189      {Histo1DPtr tmp; _h_ktd23.add(300, 400, book(tmp, 18, 1 ,1));}
190      {Histo1DPtr tmp; _h_ktd23.add(400, 500, book(tmp, 19, 1 ,1));}
191      {Histo1DPtr tmp; _h_ktd23.add(500, 600, book(tmp, 20, 1 ,1));}
192
193      {Histo1DPtr tmp; _h_cat21.add(200, 300, book(tmp, 21, 1, 1));}
194      {Histo1DPtr tmp; _h_cat21.add(300, 400, book(tmp, 22, 1, 1));}
195      {Histo1DPtr tmp; _h_cat21.add(400, 500, book(tmp, 23, 1, 1));}
196      {Histo1DPtr tmp; _h_cat21.add(500, 600, book(tmp, 24, 1, 1));}
197
198      {Histo1DPtr tmp; _h_cat32.add(200, 300, book(tmp, 25, 1, 1));}
199      {Histo1DPtr tmp; _h_cat32.add(300, 400, book(tmp, 26, 1, 1));}
200      {Histo1DPtr tmp; _h_cat32.add(400, 500, book(tmp, 27, 1, 1));}
201      {Histo1DPtr tmp; _h_cat32.add(500, 600, book(tmp, 28, 1, 1));}
202
203      {Histo1DPtr tmp; _h_ktt21.add(200, 300, book(tmp, 29, 1, 1));}
204      {Histo1DPtr tmp; _h_ktt21.add(300, 400, book(tmp, 30, 1, 1));}
205      {Histo1DPtr tmp; _h_ktt21.add(400, 500, book(tmp, 31, 1, 1));}
206      {Histo1DPtr tmp; _h_ktt21.add(500, 600, book(tmp, 32, 1, 1));}
207
208      {Histo1DPtr tmp; _h_ktt32.add(200, 300, book(tmp, 33, 1, 1));}
209      {Histo1DPtr tmp; _h_ktt32.add(300, 400, book(tmp, 34, 1, 1));}
210      {Histo1DPtr tmp; _h_ktt32.add(400, 500, book(tmp, 35, 1, 1));}
211      {Histo1DPtr tmp; _h_ktt32.add(500, 600, book(tmp, 36, 1, 1));}
212    }
213
214
215    /// Perform the per-event analysis
216    void analyze(const Event& event) {
217      const double weight = 1.0;
218
219      using namespace fastjet;
220
221      // Get anti-kt jets with p_T > 200 GeV, check abs(y) < 2, and fill mass histograms
222      const FastJets& ktfj = apply<FastJets>(event, "AKT");
223      PseudoJets ktjets = ktfj.pseudoJetsByPt(200*GeV);
224      for (const PseudoJet & ajet : ktjets) {
225        if (abs(ajet.rap()) < 2) {
226          _h_ktmass.fill(ajet.perp(), ajet.m(), weight);
227        }
228      }
229
230      // Same as above but C/A jets
231      const FastJets& cafj = apply<FastJets>(event, "CA");
232      PseudoJets cajets = cafj.pseudoJetsByPt(200*GeV);
233      for (const PseudoJet & ajet : cajets) {
234        if (abs(ajet.rap()) < 2) {
235          _h_camass.fill(ajet.perp(), ajet.m(), weight);
236        }
237      }
238
239      // Split and filter.
240      // Only do this to C/A jets in this analysis.
241      for (const PseudoJet & pjet : cajets) {
242        if ( pjet.perp() > 600 || abs(pjet.rap()) > 2) continue;
243        double dR = 0;
244        bool unclustered = false;
245        PseudoJets split_jets = splitjet(pjet, dR, cafj, unclustered);
246        if ( (dR < 0.15) || (unclustered == false) ) continue;
247        PseudoJet filt_jet = filterjet(split_jets, dR, 0.3);
248        _h_filtmass.fill(filt_jet.perp(), filt_jet.m(), weight);
249      }
250
251      // Use the two last stages of clustering to get sqrt(d_12) and sqrt(d_23).
252      // Only use anti-kt jets in this analysis.
253      for (const PseudoJet & pjet : ktjets) {
254        if (pjet.perp() > 600 || abs(pjet.rap()) > 2) continue;
255        ClusterSequence subjet_cseq(ktfj.clusterSeq()->constituents(pjet), JetDefinition(kt_algorithm, M_PI/2.));
256        double d_12 = subjet_cseq.exclusive_dmerge(1) * M_PI*M_PI/4.;
257        double d_23 = subjet_cseq.exclusive_dmerge(2) * M_PI*M_PI/4.;
258        _h_ktd12.fill(pjet.perp(), sqrt(d_12), weight);
259        _h_ktd23.fill(pjet.perp(), sqrt(d_23), weight);
260      }
261
262      // N-subjettiness, use beta = 1 (no rationale given).
263      // Uses the functions defined above.
264      // C/A jets first, anti-kt after.
265      double beta = 1.;
266      //Rcut is used for particles that are very far from the closest axis. At 10
267      //is has no impact on the outcome of the calculation
268      double Rcut = 10.;
269      for (const PseudoJet & pjet : cajets) {
270        if (pjet.perp() > 600*GeV || fabs(pjet.rap()) > 2) continue;
271
272        const PseudoJets constituents = cafj.clusterSeq()->constituents(pjet);
273        if (constituents.size() < 3) continue;
274
275        const PseudoJets axis1 = jetGetAxes(1, constituents, M_PI/2.0);
276        const PseudoJets axis2 = jetGetAxes(2, constituents, M_PI/2.0);
277        const PseudoJets axis3 = jetGetAxes(3, constituents, M_PI/2.0);
278
279        const double radius = 1.2;
280        const double tau1 = jetTauValue(beta, radius, constituents, axis1, Rcut);
281        const double tau2 = jetTauValue(beta, radius, constituents, axis2, Rcut);
282        const double tau3 = jetTauValue(beta, radius, constituents, axis3, Rcut);
283
284        if (tau1 == 0 || tau2 == 0) continue;
285        _h_cat21.fill(pjet.perp(), tau2/tau1, weight);
286        _h_cat32.fill(pjet.perp(), tau3/tau2, weight);
287      }
288
289      for (const PseudoJet & pjet : ktjets) {
290        if (pjet.perp() > 600*GeV || fabs(pjet.rap()) > 2) continue;
291
292        const PseudoJets constituents = ktfj.clusterSeq()->constituents(pjet);
293        if (constituents.size() < 3) continue;
294
295        const PseudoJets axis1 = jetGetAxes(1, constituents, M_PI/2.0);
296        const PseudoJets axis2 = jetGetAxes(2, constituents, M_PI/2.0);
297        const PseudoJets axis3 = jetGetAxes(3, constituents, M_PI/2.0);
298
299        const double radius = 1.0;
300        const double tau1 = jetTauValue(beta, radius, constituents, axis1, Rcut);
301        const double tau2 = jetTauValue(beta, radius, constituents, axis2, Rcut);
302        const double tau3 = jetTauValue(beta, radius, constituents, axis3, Rcut);
303        if (tau1 == 0 || tau2 == 0) continue;
304
305        _h_ktt21.fill(pjet.perp(), tau2/tau1, weight);
306        _h_ktt32.fill(pjet.perp(), tau3/tau2, weight);
307      }
308    }
309
310
311    /// Normalise histograms etc., after the run
312    void finalize() {
313      for (Histo1DPtr h : _h_camass.histos()) normalize(h);
314      for (Histo1DPtr h : _h_filtmass.histos()) normalize(h);
315      for (Histo1DPtr h : _h_ktmass.histos()) normalize(h);
316      for (Histo1DPtr h : _h_ktd12.histos()) normalize(h);
317      for (Histo1DPtr h : _h_ktd23.histos()) normalize(h);
318      for (Histo1DPtr h : _h_cat21.histos()) normalize(h);
319      for (Histo1DPtr h : _h_cat32.histos()) normalize(h);
320      for (Histo1DPtr h : _h_ktt21.histos()) normalize(h);
321      for (Histo1DPtr h : _h_ktt32.histos()) normalize(h);
322    }
323
324  private:
325
326    BinnedHistogram _h_camass;
327    BinnedHistogram _h_filtmass;
328    BinnedHistogram _h_ktmass;
329    BinnedHistogram _h_ktd12;
330    BinnedHistogram _h_ktd23;
331    BinnedHistogram _h_cat21;
332    BinnedHistogram _h_cat32;
333    BinnedHistogram _h_ktt21;
334    BinnedHistogram _h_ktt32;
335
336  };
337
338
339  // The hook for the plugin system
340  RIVET_DECLARE_PLUGIN(ATLAS_2012_I1094564);
341
342}