rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

ATLAS_2020_I1808726

Hadronic event shapes in multijet final states
Experiment: ATLAS (LHC)
Inspire ID: 1808726
Status: VALIDATED
Authors:
  • Javier Llorente
  • Deepak Kar
References: Beams: p+ p+
Beam energies: (6500.0, 6500.0) GeV
Run details:
  • Jet production at 13000 GeV.

A measurement of event-shape variables in proton-proton collisions at large momentum transfer is presented using data collected at $\sqrt{s}$=13 TeV with the ATLAS detector at the Large Hadron Collider. Six event-shape variables calculated using hadronic jets are studied in inclusive multijet events using data corresponding to an integrated luminosity of 139 fb$^{-1}$. Measurements are performed in bins of jet multiplicity and in different ranges of the scalar sum of the transverse momenta of the two leading jets, reaching scales beyond 2 TeV. These measurements are compared with predictions from Monte Carlo event generators containing leading-order or next-to-leading order matrix elements matched to parton showers simulated to leading-logarithm accuracy. At low jet multiplicities, shape discrepancies between the measurements and the Monte Carlo predictions are observed. At high jet multiplicities, the shapes are better described but discrepancies in the normalisation are observed.

Source code: ATLAS_2020_I1808726.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5#include "Rivet/Projections/Thrust.hh"
  6#include "Rivet/Projections/Sphericity.hh"
  7
  8
  9namespace Rivet {
 10
 11  /// @brief Multijet event shapes at 13 TeV
 12  class ATLAS_2020_I1808726 : public Analysis {
 13  public:
 14
 15    /// Constructor
 16    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2020_I1808726);
 17
 18    double xs1 = 0.0;
 19    double xs2 = 0.0;
 20    double xs3 = 0.0;
 21
 22    /// Initialization, called once before running
 23    void init() {
 24
 25      //Jet collection (excluding muons and neutrinos)
 26      const FinalState fs(Cuts::abseta < 4.5);
 27
 28      FastJets jets(fs, JetAlg::ANTIKT, 0.4, JetMuons::NONE, JetInvisibles::NONE);
 29      declare(jets, "Jets");
 30
 31      // Book histograms
 32      //Jet multiplicity
 33      book(_h["njet_h1"], 73, 1, 1);
 34      book(_h["njet_h2"], 74, 1, 1);
 35      book(_h["njet_h3"], 75, 1, 1);
 36
 37      //Transverse thrust
 38      book(_h["transThrust_j3_h1"], 1, 1, 1);
 39      book(_h["transThrust_j3_h2"], 5, 1, 1);
 40      book(_h["transThrust_j3_h3"], 9, 1, 1);
 41      book(_h["transThrust_j4_h1"], 2, 1, 1);
 42      book(_h["transThrust_j4_h2"], 6, 1, 1);
 43      book(_h["transThrust_j4_h3"], 10, 1, 1);
 44      book(_h["transThrust_j5_h1"], 3, 1, 1);
 45      book(_h["transThrust_j5_h2"], 7, 1, 1);
 46      book(_h["transThrust_j5_h3"], 11, 1, 1);
 47      book(_h["transThrust_j6_h1"], 4, 1, 1);
 48      book(_h["transThrust_j6_h2"], 8, 1, 1);
 49      book(_h["transThrust_j6_h3"], 12, 1, 1);
 50
 51      //Thrust minor
 52      book(_h["transMinor_j3_h1"], 13, 1, 1);
 53      book(_h["transMinor_j3_h2"], 17, 1, 1);
 54      book(_h["transMinor_j3_h3"], 21, 1, 1);
 55      book(_h["transMinor_j4_h1"], 14, 1, 1);
 56      book(_h["transMinor_j4_h2"], 18, 1, 1);
 57      book(_h["transMinor_j4_h3"], 22, 1, 1);
 58      book(_h["transMinor_j5_h1"], 15, 1, 1);
 59      book(_h["transMinor_j5_h2"], 19, 1, 1);
 60      book(_h["transMinor_j5_h3"], 23, 1, 1);
 61      book(_h["transMinor_j6_h1"], 16, 1, 1);
 62      book(_h["transMinor_j6_h2"], 20, 1, 1);
 63      book(_h["transMinor_j6_h3"], 24, 1, 1);
 64
 65      //Transverse sphericity
 66      book(_h["transSphericity_j3_h1"], 25, 1, 1);
 67      book(_h["transSphericity_j3_h2"], 29, 1, 1);
 68      book(_h["transSphericity_j3_h3"], 33, 1, 1);
 69      book(_h["transSphericity_j4_h1"], 26, 1, 1);
 70      book(_h["transSphericity_j4_h2"], 30, 1, 1);
 71      book(_h["transSphericity_j4_h3"], 34, 1, 1);
 72      book(_h["transSphericity_j5_h1"], 27, 1, 1);
 73      book(_h["transSphericity_j5_h2"], 31, 1, 1);
 74      book(_h["transSphericity_j5_h3"], 35, 1, 1);
 75      book(_h["transSphericity_j6_h1"], 28, 1, 1);
 76      book(_h["transSphericity_j6_h2"], 32, 1, 1);
 77      book(_h["transSphericity_j6_h3"], 36, 1, 1);
 78
 79      //Aplanarity
 80      book(_h["aplanarity_j3_h1"], 37, 1, 1);
 81      book(_h["aplanarity_j3_h2"], 41, 1, 1);
 82      book(_h["aplanarity_j3_h3"], 45, 1, 1);
 83      book(_h["aplanarity_j4_h1"], 38, 1, 1);
 84      book(_h["aplanarity_j4_h2"], 42, 1, 1);
 85      book(_h["aplanarity_j4_h3"], 46, 1, 1);
 86      book(_h["aplanarity_j5_h1"], 39, 1, 1);
 87      book(_h["aplanarity_j5_h2"], 43, 1, 1);
 88      book(_h["aplanarity_j5_h3"], 47, 1, 1);
 89      book(_h["aplanarity_j6_h1"], 40, 1, 1);
 90      book(_h["aplanarity_j6_h2"], 44, 1, 1);
 91      book(_h["aplanarity_j6_h3"], 48, 1, 1);
 92
 93      //C
 94      book(_h["C_j3_h1"], 49, 1, 1);
 95      book(_h["C_j3_h2"], 53, 1, 1);
 96      book(_h["C_j3_h3"], 57, 1, 1);
 97      book(_h["C_j4_h1"], 50, 1, 1);
 98      book(_h["C_j4_h2"], 54, 1, 1);
 99      book(_h["C_j4_h3"], 58, 1, 1);
100      book(_h["C_j5_h1"], 51, 1, 1);
101      book(_h["C_j5_h2"], 55, 1, 1);
102      book(_h["C_j5_h3"], 59, 1, 1);
103      book(_h["C_j6_h1"], 52, 1, 1);
104      book(_h["C_j6_h2"], 56, 1, 1);
105      book(_h["C_j6_h3"], 60, 1, 1);
106
107      //D
108      book(_h["D_j3_h1"], 61, 1, 1);
109      book(_h["D_j3_h2"], 65, 1, 1);
110      book(_h["D_j3_h3"], 69, 1, 1);
111      book(_h["D_j4_h1"], 62, 1, 1);
112      book(_h["D_j4_h2"], 66, 1, 1);
113      book(_h["D_j4_h3"], 70, 1, 1);
114      book(_h["D_j5_h1"], 63, 1, 1);
115      book(_h["D_j5_h2"], 67, 1, 1);
116      book(_h["D_j5_h3"], 71, 1, 1);
117      book(_h["D_j6_h1"], 64, 1, 1);
118      book(_h["D_j6_h2"], 68, 1, 1);
119      book(_h["D_j6_h3"], 72, 1, 1);
120    }
121
122    void analyze(const Event& event) {
123
124      const Jets& jets = apply<FastJets>(event, "Jets").jetsByPt(Cuts::pT > 7.0*GeV);
125
126      //Select jets passing kinematic cuts
127      std::vector<const Jet*> goodJets; goodJets.clear();
128      std::vector<Vector3> momenta2; momenta2.clear();
129      std::vector<Vector3> momenta3; momenta3.clear();
130
131      //foreach (const Jet& j, jets) {
132      for (const Jet& j : jets) {
133        if (j.abseta() < 2.4 && j.pt() > 100.0*GeV){
134	         goodJets.push_back(&j);
135
136             Vector3 jet2 = j.p3();
137             jet2.setZ(0.0);
138              momenta2.push_back(jet2);
139
140	          Vector3 jet3 = j.p3();
141	          momenta3.push_back(jet3);
142	       }
143        }
144
145      //Dijet event selection
146      if (goodJets.size() < 2) vetoEvent;
147      double ht2 = goodJets[0]->pt()+goodJets[1]->pt();
148      if (ht2 <= 1000.0*GeV) vetoEvent;
149
150      //Jet multiplicity
151      if (ht2 > 1000.0*GeV && ht2 < 1500.0*GeV) _h["njet_h1"]->fill(goodJets.size());
152      if (ht2 > 1500.0*GeV && ht2 < 2000.0*GeV) _h["njet_h2"]->fill(goodJets.size());
153      if (ht2 > 2000.0*GeV) _h["njet_h3"]->fill(goodJets.size());
154
155      //Thrust calculation
156      Thrust thrust;
157      thrust.calc(momenta2);
158      const double transThrust  = 1.0 - thrust.thrust();
159      const double transMinor = thrust.thrustMajor();
160
161      //Linearized sphericity calculation (2D)
162      double a11 = 0.0; double a22 = 0.0;
163      double a12 = 0.0;
164      double modSum2 = 0.0;
165
166      for (size_t k = 0; k < momenta2.size(); ++k) {
167        modSum2 += momenta2[k].mod();
168	      a11 += momenta2[k].x()*momenta2[k].x()/momenta2[k].mod();
169        a22 += momenta2[k].y()*momenta2[k].y()/momenta2[k].mod();
170        a12 += momenta2[k].x()*momenta2[k].y()/momenta2[k].mod();
171      }
172
173      double trc2 = (a11+a22)/modSum2;
174      double det2 = (a11*a22-a12*a12)/pow(modSum2,2);
175
176      double eigen21 = (trc2+sqrt(pow(trc2,2)-4*det2))/2;
177      double eigen22 = (trc2-sqrt(pow(trc2,2)-4*det2))/2;
178      double transSphericity = 2*eigen22/(eigen21+eigen22);
179
180      //Linearized sphericity calculation (3D)
181      double b11 = 0.0; double b12 = 0.0; double b13 = 0.0;
182      double b22 = 0.0; double b23 = 0.0;
183      double b33 = 0.0;
184      double modSum3 = 0.0;
185
186      for (size_t k = 0; k < momenta3.size(); ++k){
187        modSum3 += momenta3[k].mod();
188        b11 += momenta3[k].x()*momenta3[k].x()/momenta3[k].mod();
189        b22 += momenta3[k].y()*momenta3[k].y()/momenta3[k].mod();
190        b33 += momenta3[k].z()*momenta3[k].z()/momenta3[k].mod();
191        b12 += momenta3[k].x()*momenta3[k].y()/momenta3[k].mod();
192        b13 += momenta3[k].x()*momenta3[k].z()/momenta3[k].mod();
193        b23 += momenta3[k].y()*momenta3[k].z()/momenta3[k].mod();
194      }
195
196      Matrix3 sph3;
197      sph3.set(0,0, b11/modSum3); sph3.set(0,1, b12/modSum3); sph3.set(0,2, b13/modSum3);
198      sph3.set(1,0, b12/modSum3); sph3.set(1,1, b22/modSum3); sph3.set(1,2, b23/modSum3);
199      sph3.set(2,0, b13/modSum3); sph3.set(2,1, b23/modSum3); sph3.set(2,2, b33/modSum3);
200
201      double q = sph3.trace()/3.;
202      double p1 = sph3.get(0,1)*sph3.get(0,1) + sph3.get(0,2)*sph3.get(0,2) + sph3.get(1,2)*sph3.get(1,2);
203      double p2 = (sph3.get(0,0)-q)*(sph3.get(0,0)-q) + (sph3.get(1,1)-q)*(sph3.get(1,1)-q) + (sph3.get(2,2)-q)*(sph3.get(2,2)-q) + 2*p1;
204      double p = sqrt(p2/6.);
205
206      Matrix3 I3 = Matrix3::mkIdentity();
207      double r = ( 1./p * (sph3 - q*I3)).det()/2.;
208
209      double phi(0);
210      if (r <= -1) phi = M_PI / 3.;
211      else if (r >= 1) phi = 0;
212      else phi = acos(r) / 3.;
213
214      double eigen31 = q + 2 * p * cos(phi);
215      double eigen33 = q + 2 * p * cos(phi + (2*M_PI/3.));
216      double eigen32 = 3 * q - eigen31 - eigen33;
217
218      double aplanarity = (3./2)*eigen33;
219      double C = 3*(eigen31*eigen32 + eigen31*eigen33 + eigen32*eigen33);
220      double D = 27*eigen31*eigen32*eigen33;
221
222      //Fill event-shape histograms
223      if (ht2 > 1000.0*GeV && ht2 < 1500.0*GeV){
224
225      	if (goodJets.size() == 3){
226	       _h["transThrust_j3_h1"]->fill(transThrust); _h["transMinor_j3_h1"]->fill(transMinor);
227	       _h["transSphericity_j3_h1"]->fill(transSphericity); _h["aplanarity_j3_h1"]->fill(aplanarity);
228	       _h["C_j3_h1"]->fill(C); _h["D_j3_h1"]->fill(D);
229	      }
230
231	   if (goodJets.size() == 4){
232          _h["transThrust_j4_h1"]->fill(transThrust); _h["transMinor_j4_h1"]->fill(transMinor);
233          _h["transSphericity_j4_h1"]->fill(transSphericity); _h["aplanarity_j4_h1"]->fill(aplanarity);
234          _h["C_j4_h1"]->fill(C); _h["D_j4_h1"]->fill(D);
235	     }
236
237	   if (goodJets.size() == 5){
238          _h["transThrust_j5_h1"]->fill(transThrust); _h["transMinor_j5_h1"]->fill(transMinor);
239          _h["transSphericity_j5_h1"]->fill(transSphericity); _h["aplanarity_j5_h1"]->fill(aplanarity);
240          _h["C_j5_h1"]->fill(C); _h["D_j5_h1"]->fill(D);
241     	}
242
243	   if (goodJets.size() >= 6){
244          _h["transThrust_j6_h1"]->fill(transThrust); _h["transMinor_j6_h1"]->fill(transMinor);
245          _h["transSphericity_j6_h1"]->fill(transSphericity); _h["aplanarity_j6_h1"]->fill(aplanarity);
246          _h["C_j6_h1"]->fill(C); _h["D_j6_h1"]->fill(D);
247	    }
248      }
249
250
251      if (ht2 > 1500.0*GeV && ht2 < 2000.0*GeV){
252
253        if (goodJets.size() == 3){
254          _h["transThrust_j3_h2"]->fill(transThrust); _h["transMinor_j3_h2"]->fill(transMinor);
255          _h["transSphericity_j3_h2"]->fill(transSphericity); _h["aplanarity_j3_h2"]->fill(aplanarity);
256          _h["C_j3_h2"]->fill(C); _h["D_j3_h2"]->fill(D);
257        }
258
259        if (goodJets.size() == 4){
260          _h["transThrust_j4_h2"]->fill(transThrust); _h["transMinor_j4_h2"]->fill(transMinor);
261          _h["transSphericity_j4_h2"]->fill(transSphericity); _h["aplanarity_j4_h2"]->fill(aplanarity);
262          _h["C_j4_h2"]->fill(C); _h["D_j4_h2"]->fill(D);
263        }
264
265        if (goodJets.size() == 5){
266          _h["transThrust_j5_h2"]->fill(transThrust); _h["transMinor_j5_h2"]->fill(transMinor);
267          _h["transSphericity_j5_h2"]->fill(transSphericity); _h["aplanarity_j5_h2"]->fill(aplanarity);
268          _h["C_j5_h2"]->fill(C); _h["D_j5_h2"]->fill(D);
269        }
270
271        if (goodJets.size() >= 6){
272          _h["transThrust_j6_h2"]->fill(transThrust); _h["transMinor_j6_h2"]->fill(transMinor);
273          _h["transSphericity_j6_h2"]->fill(transSphericity); _h["aplanarity_j6_h2"]->fill(aplanarity);
274          _h["C_j6_h2"]->fill(C); _h["D_j6_h2"]->fill(D);
275        }
276      }
277
278      if (ht2 > 2000.0*GeV){
279
280        if (goodJets.size() == 3){
281          _h["transThrust_j3_h3"]->fill(transThrust); _h["transMinor_j3_h3"]->fill(transMinor);
282          _h["transSphericity_j3_h3"]->fill(transSphericity); _h["aplanarity_j3_h3"]->fill(aplanarity);
283          _h["C_j3_h3"]->fill(C); _h["D_j3_h3"]->fill(D);
284        }
285
286        if (goodJets.size() == 4){
287          _h["transThrust_j4_h3"]->fill(transThrust); _h["transMinor_j4_h3"]->fill(transMinor);
288          _h["transSphericity_j4_h3"]->fill(transSphericity); _h["aplanarity_j4_h3"]->fill(aplanarity);
289          _h["C_j4_h3"]->fill(C); _h["D_j4_h3"]->fill(D);
290        }
291
292        if (goodJets.size() == 5){
293          _h["transThrust_j5_h3"]->fill(transThrust); _h["transMinor_j5_h3"]->fill(transMinor);
294          _h["transSphericity_j5_h3"]->fill(transSphericity); _h["aplanarity_j5_h3"]->fill(aplanarity);
295          _h["C_j5_h3"]->fill(C); _h["D_j5_h3"]->fill(D);
296        }
297
298        if (goodJets.size() >= 6){
299          _h["transThrust_j6_h3"]->fill(transThrust); _h["transMinor_j6_h3"]->fill(transMinor);
300          _h["transSphericity_j6_h3"]->fill(transSphericity); _h["aplanarity_j6_h3"]->fill(aplanarity);
301          _h["C_j6_h3"]->fill(C); _h["D_j6_h3"]->fill(D);
302        }
303      }
304    }
305
306
307    void finalize() {
308
309      const double xs1 = _h["njet_h1"]->sumW();
310      const double xs2 = _h["njet_h2"]->sumW();
311      const double xs3 = _h["njet_h3"]->sumW();
312
313      for (auto& hist : _h) {
314        if (hist.first.find("njet_") != string::npos)  scale(hist.second, crossSectionPerEvent()/picobarn);
315        else if (hist.first.find("_h1") != string::npos) scale(hist.second, 1.0/xs1);
316        else if (hist.first.find("_h2") != string::npos) scale(hist.second, 1.0/xs2);
317        else if (hist.first.find("_h3") != string::npos) scale(hist.second, 1.0/xs3);
318      }
319    }
320
321  private:
322
323    //Jet multiplicity
324    map<string,Histo1DPtr> _h;
325
326  };
327
328  RIVET_DECLARE_PLUGIN(ATLAS_2020_I1808726);
329}