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