rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

ATLAS_2017_I1509919

Track-based underlying event at 13 TeV in ATLAS
Experiment: ATLAS (LHC)
Inspire ID: 1509919
Status: VALIDATED
Authors:
  • Roman Lysak
References: Beams: p+ p+
Beam energies: (6500.0, 6500.0) GeV
Run details:
  • $pp$ QCD interactions at 13 TeV. Diffractive events should be included.

We present charged-particle distributions sensitive to the underlying event, measured by the ATLAS detector in proton-proton collisions at a centre-of-mass energy of 13 TeV, in low-luminosity Large Hadron Collider fills corresponding to an integrated luminosity of 1.6 nb${}^{-1}$. The distributions were constructed using charged particles with absolute pseudorapidity less than 2.5 and with transverse momentum greater than 500 MeV, in events with at least one such charged particle with transverse momentum above 1 GeV. These distributions characterise the angular distribution of energy and particle flows with respect to the charged particle with highest transverse momentum, as a function of both that momentum and of charged-particle multiplicity. The results have been corrected for detector effects and are compared to the predictions of various Monte Carlo event generators, experimentally establishing the level of underlying-event activity at LHC Run 2 energies and providing inputs for the development of event generator modelling. The current models in use for UE modelling typically describe this data to 5% accuracy, compared with data uncertainties of less than 1%.

Source code: ATLAS_2017_I1509919.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/ChargedFinalState.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// Track-based underlying event at 13 TeV in ATLAS
  9  class ATLAS_2017_I1509919 : public Analysis {
 10  public:
 11
 12    // Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2017_I1509919);
 14
 15
 16    // Pre-run histogram and projection booking
 17    void init() {
 18
 19      declare(ChargedFinalState(Cuts::abseta < 2.5 && Cuts::pT > 500*MeV), "CFS500");
 20
 21      // Nch profiles vs. pT_lead
 22      book(_hist_nch[0], 22, 1, 1);
 23      book(_hist_nch[1], 23, 1, 1);
 24      book(_hist_nch[2], 21, 1, 1);
 25      book(_hist_nch[3],  3, 1, 1);
 26      book(_hist_nch[4],  2, 1, 1);
 27      book(_hist_nch[5],  4, 1, 1);
 28
 29      // pTsum profiles vs. pT_lead
 30      book(_hist_ptsum[0], 25, 1, 1);
 31      book(_hist_ptsum[1], 26, 1, 1);
 32      book(_hist_ptsum[2], 24, 1, 1);
 33      book(_hist_ptsum[3],  6, 1, 1);
 34      book(_hist_ptsum[4],  5, 1, 1);
 35      book(_hist_ptsum[5],  7, 1, 1);
 36
 37      // <pT> profiles vs pT_lead (not measured for trans diff)
 38      book(_hist_ptavg[0], 29, 1, 1);
 39      book(_hist_ptavg[1], 30, 1, 1);
 40      book(_hist_ptavg[2], 11, 1, 1);
 41      book(_hist_ptavg[3], 13, 1, 1);
 42      book(_hist_ptavg[4], 12, 1, 1);
 43
 44      // <pT> profiles vs. Nch (not measured for trans diff)
 45      book(_hist_dn_dpt[0], 27, 1, 1);
 46      book(_hist_dn_dpt[1], 28, 1, 1);
 47      book(_hist_dn_dpt[2],  8, 1, 1);
 48      book(_hist_dn_dpt[3], 10, 1, 1);
 49      book(_hist_dn_dpt[4],  9, 1, 1);
 50
 51      // Only measured for trans max/min
 52      book(_hist_dn_dpt2[3], 32, 1, 1);
 53      book(_hist_dn_dpt2[4], 31, 1, 1);
 54
 55      // Nch vs. Delta(phi) profiles
 56      book(_hist_N_vs_dPhi[0], 15, 1, 1);
 57      book(_hist_N_vs_dPhi[1], 16, 1, 1);
 58      book(_hist_N_vs_dPhi[2], 17, 1, 1);
 59
 60      // pT vs. Delta(phi) profiles
 61      book(_hist_pT_vs_dPhi[0], 18, 1, 1);
 62      book(_hist_pT_vs_dPhi[1], 19, 1, 1);
 63      book(_hist_pT_vs_dPhi[2], 20, 1, 1);
 64
 65      //ptLead histos only for 1 and 5 GeV cuts
 66      book(_hist_ptLead[0],  1, 1, 1);
 67      book(_hist_ptLead[1], 14, 1, 1);
 68
 69      for (size_t iC = 0; iC < NCUTS; ++iC) {
 70        book(_counters[iC], "Ctr_cut_" + toString(iC));
 71      }
 72
 73    }
 74
 75
 76    void analyze(const Event& event) {
 77
 78      // Get charged particles (tracks) with pT > 500 MeV
 79      const ChargedFinalState& charged500 = apply<ChargedFinalState>(event, "CFS500");
 80      const Particles& particlesAll = charged500.particlesByPt();
 81      MSG_DEBUG("Num tracks: " << particlesAll.size());
 82
 83      const Cut& pcut = ( (Cuts::abspid != PID::SIGMAMINUS) && (Cuts::abspid != PID::SIGMAPLUS) &&
 84                          (Cuts::abspid != PID::XIMINUS)    && (Cuts::abspid != PID::OMEGAMINUS) );
 85      const Particles& particles = charged500.particlesByPt(pcut);
 86      MSG_DEBUG("Num tracks without strange baryons: " <<  particles.size());
 87
 88      // Require at least one track in the event for pTlead histograms
 89      if (particles.empty()) vetoEvent;
 90      for (size_t iC = 0; iC < 2; ++iC) {
 91        if (particles[0].pT() < PTCUTS[iC]*GeV) continue;
 92        _counters[iC]->fill();
 93        _hist_ptLead[iC]->fill( particles[0].pT()/GeV);
 94      }
 95
 96      // Require at least one track in the event with pT >= 1 GeV for the rest
 97      if (particles[0].pT() < 1*GeV) vetoEvent;
 98
 99      // Identify leading track and its phi and pT
100      const Particle& p_lead = particles[0];
101      const double philead = p_lead.phi();
102      const double etalead = p_lead.eta();
103      const double pTlead  = p_lead.perp();
104      MSG_DEBUG("Leading track: pT = " << pTlead << ", eta = " << etalead << ", phi = " << philead);
105
106      // Iterate over all particles and count particles and scalar pTsum in three basic regions
107      vector<double> num(NREGIONS, 0), ptSum(NREGIONS, 0.0), avgpt(NREGIONS, 0.0);
108
109      // Temporary histos that bin Nch and pT in dPhi.
110      Histo1D hist_num_dphi(*_hist_N_vs_dPhi[0], "/hist_num_dphi");
111      Histo1D hist_pt_dphi(*_hist_pT_vs_dPhi[0], "/hist_pt_dphi");
112      hist_num_dphi.reset();
113      hist_pt_dphi .reset();
114
115      int    tmpnch[2]   = {0,0};
116      double tmpptsum[2] = {0,0};
117      for (const Particle& p : particles) {
118        const double pT   = p.pT()/GeV;
119        const double dPhi = deltaPhi(philead, p.phi()); // in range (0,pi)
120        const int    ir   = region_index(dPhi); // gives just toward/away/trans
121
122        // Toward/away/trans region: just count
123        num  [ir] += 1;
124        ptSum[ir] += pT;
125
126        // Determine which transverse side
127        if (ir == kTrans) {
128          const size_t iside = (mapAngleMPiToPi(p.phi() - philead) > 0) ? 0 : 1;
129          tmpnch  [iside] += 1;
130          tmpptsum[iside] += p.pT();
131        }
132
133        // Fill temp histos to bin Nch and pT in dPhi
134        if (p.genParticle() != p_lead.genParticle()) { // We don't want to fill all those zeros from the leading track...
135          hist_num_dphi.fill(dPhi/M_PI*180);
136          hist_pt_dphi .fill(dPhi/M_PI*180, pT/GeV);
137        }
138      }
139
140      // Construct max/min/diff regions
141      num[kTransMax ]   = std::max(tmpnch[0], tmpnch[1]);
142      num[kTransMin ]   = std::min(tmpnch[0], tmpnch[1]);
143      num[kTransDiff]   = num[kTransMax ] - num[kTransMin ];
144      ptSum[kTransMax ] = std::max(tmpptsum[0], tmpptsum[1]);
145      ptSum[kTransMin ] = std::min(tmpptsum[0], tmpptsum[1]);
146      ptSum[kTransDiff] = ptSum[kTransMax ] - ptSum[kTransMin ];
147      avgpt[kToward]    = (num[kToward] > 0 ) ? ptSum[kToward] / num[kToward] : 0. ;
148      avgpt[kAway]      = (num[kAway  ] > 0 ) ? ptSum[kAway]   / num[kAway]   : 0. ;
149      avgpt[kTrans]     = (num[kTrans ] > 0 ) ? ptSum[kTrans]  / num[kTrans]  : 0. ;
150      // Avg pt max/min regions determined according sumpt max/min
151      int sumptMaxRegID = (tmpptsum[0] >  tmpptsum[1]) ? 0 : 1 ;
152      int sumptMinRegID = (sumptMaxRegID == 0) ? 1 : 0;
153      avgpt[kTransMax ] = (tmpnch[sumptMaxRegID] > 0) ? tmpptsum[sumptMaxRegID] / tmpnch[sumptMaxRegID] : 0.;
154      avgpt[kTransMin ] = (tmpnch[sumptMinRegID] > 0) ? tmpptsum[sumptMinRegID] / tmpnch[sumptMinRegID] : 0.;
155      avgpt[kTransDiff] = ((tmpnch[sumptMaxRegID] > 0) && (tmpnch[sumptMinRegID] > 0)) ? avgpt[kTransMax ] - avgpt[kTransMin ] : 0.;
156
157
158      // Now fill underlying event histograms
159
160      // The densities are calculated by dividing the UE properties by dEta*dPhi
161      // -- each basic region has a dPhi of 2*PI/3 and dEta is two times 2.5
162      // min/max/diff regions are only half of that
163      const double dEtadPhi[NREGIONS] = { 2*2.5 * 2*PI/3.0, 2*2.5 * 2*PI/3.0, 2*2.5 * 2*PI/3.0,
164                                          2*2.5 *   PI/3.0, 2*2.5 *   PI/3.0, 2*2.5 *   PI/3.0 };
165      for (size_t iR = 0; iR < NREGIONS; ++iR) {
166
167        _hist_nch  [iR]->fill(pTlead/GeV, num[iR]   /dEtadPhi[iR]     );
168        _hist_ptsum[iR]->fill(pTlead/GeV, ptSum[iR] /GeV/dEtadPhi[iR] );
169
170        // <pT> profiles vs. pT_lead (first 3 are the same!)
171        switch (iR) {
172        case kToward    :
173        case kAway      :
174        case kTrans     :
175          if (num[iR] > 0) _hist_ptavg[iR]->fill(pTlead/GeV, avgpt[iR]/GeV);
176          break;
177        case kTransMax  :
178          if (tmpnch[sumptMaxRegID] > 0) _hist_ptavg[iR]->fill(pTlead/GeV, avgpt[iR]/GeV);
179          break;
180        case kTransMin  :
181          if (tmpnch[sumptMinRegID] > 0) _hist_ptavg[iR]->fill(pTlead/GeV, avgpt[iR]/GeV);
182          break;
183        case kTransDiff :
184          break;
185        default: //should not get here!!!
186          MSG_WARNING("Unknown region in <pT> profiles vs.pt lead switch!!! : " << iR);
187        }
188
189        // <pT> profiles vs. Nch (first 3 are the same!)
190        switch (iR) {
191        case kToward    :
192        case kAway      :
193        case kTrans     :
194          if (num[iR] > 0) _hist_dn_dpt[iR]->fill(num[iR] , avgpt[iR]/GeV);
195          break;
196        case kTransMax  :
197          if (tmpnch[sumptMaxRegID] > 0) {
198            _hist_dn_dpt [iR]->fill(num[kTrans]          , avgpt[iR]/GeV);
199            _hist_dn_dpt2[iR]->fill(tmpnch[sumptMaxRegID], avgpt[iR]/GeV);
200          }
201          break;
202        case kTransMin  :
203          if (tmpnch[sumptMinRegID] > 0) {
204            _hist_dn_dpt [iR]->fill(num[kTrans]          , avgpt[iR]/GeV);
205            _hist_dn_dpt2[iR]->fill(tmpnch[sumptMinRegID], avgpt[iR]/GeV);
206          }
207          break;
208        case kTransDiff :
209          break;
210        default : //should not get here!!!
211          MSG_INFO("unknown region in <pT> profiles vs. nch switch!!! : " <<  iR);
212        }
213
214      }
215
216
217      // Update the "proper" dphi profile histograms
218      // Note that we fill dN/dEtadPhi: dEta = 2*2.5, dPhi = 2*PI/nBins
219      const double dEtadPhi2 = (2*2.5 * 2) * M_PI/180.;
220      for (size_t i = 0; i < hist_num_dphi.numBins(); ++i) {
221
222        // First Nch
223        double mean = hist_num_dphi.bin(i).xMid() ;
224        double value = 0.;
225        if (hist_num_dphi.bin(i).numEntries() > 0) {
226          mean  = hist_num_dphi.bin(i).xMean() ;
227          value = hist_num_dphi.bin(i).area()/hist_num_dphi.bin(i).xWidth()/dEtadPhi2;
228        }
229        for (size_t iC = 0; iC < NCUTS; ++iC) {
230          if (pTlead >= PTCUTS[iC]*GeV) _hist_N_vs_dPhi[iC] ->fill(mean, value);
231        }
232
233        // Then pT
234        mean = hist_pt_dphi.bin(i).xMid();
235        value = 0.;
236        if (hist_pt_dphi.bin(i).numEntries() > 0) {
237          mean  = hist_pt_dphi.bin(i).xMean() ;
238          value = hist_pt_dphi.bin(i).area()/hist_pt_dphi.bin(i).xWidth()/dEtadPhi2;
239        }
240        for (size_t iC = 0; iC < NCUTS; ++iC) {
241          if (pTlead >= PTCUTS[iC]*GeV) _hist_pT_vs_dPhi[iC] ->fill(mean, value);
242        }
243      }
244
245    }
246
247
248    void finalize() {
249      for (size_t iC = 0; iC < NCUTS; ++iC) {
250        if (iC == 0 || iC == 1) scale(_hist_ptLead[iC], 1.0/_counters[iC]->sumW());
251      }
252    }
253
254
255  private:
256
257    enum regionID {
258      kToward = 0,
259      kAway,
260      kTrans,
261      kTransMax,
262      kTransMin,
263      kTransDiff,
264      NREGIONS
265    };
266
267    // Little helper function to identify basic Delta(phi) regions: toward/away/trans
268    int region_index(double dphi) {
269      assert(inRange(dphi, 0.0, PI, CLOSED, CLOSED));
270      if (dphi < PI/3.0)    return kToward;
271      if (dphi < 2*PI/3.0)  return kTrans;
272      return kAway;
273    }
274
275    const static size_t NCUTS = 3;
276    const vector<double> PTCUTS = { 1., 5., 10. };
277
278
279    /// @name Histograms
280    //@{
281
282    // Nch, sumpT, avgpT profiles vs. pTlead
283    Profile1DPtr _hist_nch   [NREGIONS]; //for regions: all 6 regions
284    Profile1DPtr _hist_ptsum [NREGIONS]; //for regions: all 6 regions
285    Profile1DPtr _hist_ptavg [NREGIONS]; //for regions: trans towards/away/all/min/max
286    // Nch, sumpT, avgpT profiles vs. Nch
287    Profile1DPtr _hist_dn_dpt [NREGIONS]; //regions: towards/away/ vs nch(region) & trans all/min/max vs nch(trans)
288    Profile1DPtr _hist_dn_dpt2[NREGIONS]; //regions: trans min/max vs. nch(region)
289
290    Profile1DPtr _hist_N_vs_dPhi [NCUTS];
291    Profile1DPtr _hist_pT_vs_dPhi[NCUTS];
292    Histo1DPtr   _hist_ptLead[NCUTS]; //for 1,5 GeV cuts only
293    CounterPtr   _counters[NCUTS];
294
295    //@}
296
297  };
298
299
300  RIVET_DECLARE_PLUGIN(ATLAS_2017_I1509919);
301
302}