rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

ATLAS_2010_S8894728

Track-based underlying event at 900 GeV and 7 TeV in ATLAS
Experiment: ATLAS (LHC)
Inspire ID: 879407
Status: VALIDATED
Authors:
  • Andy Buckley
  • Holger Schulz
References: Beams: p+ p+
Beam energies: (450.0, 450.0); (3500.0, 3500.0) GeV
Run details:
  • pp QCD interactions at 900 GeV and 7 TeV. Diffractive events should be included, but only influence the lowest bins. Multiple kinematic cuts should not be required.

The underlying event measurements with the ATLAS detector at the LHC at the center of mass energies of 900 GeV and 7 TeV. The observables sensitive to the underlying event, i.e the charged particle density and charged pT sum, as well as their standard deviations and the average pT, are measured as functions of the leading track. A track pT cut of 500 MeV is applied for most observables, but the main profile plots are also shown for a lower track cut of 100 MeV, which includes much more of the soft cross-section. The angular distribution of the charged tracks with respect to the leading track is also studied, as are the correlation between mean transverse momentum and charged particle multiplicity, and the `plateau' height as a function of the leading track $|\eta|$.

Source code: ATLAS_2010_S8894728.cc
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/ChargedFinalState.hh"

namespace Rivet {


  class ATLAS_2010_S8894728 : public Analysis {
  public:

    ATLAS_2010_S8894728() : Analysis("ATLAS_2010_S8894728") {    }


    void init() {
      const ChargedFinalState cfs100(-2.5, 2.5, 100*MeV);
      declare(cfs100, "CFS100");
      const ChargedFinalState cfs500(-2.5, 2.5, 500*MeV);
      declare(cfs500, "CFS500");
      const ChargedFinalState cfslead(-2.5, 2.5, 1.0*GeV);
      declare(cfslead, "CFSlead");

      // Get an index for the beam energy
      int isqrts = -1;
      if (fuzzyEquals(sqrtS(), 900*GeV)) isqrts = 0;
      else if (fuzzyEquals(sqrtS(), 7*TeV)) isqrts = 1;
      assert(isqrts >= 0);

      // Nch profiles, 500 MeV track pT cut
      _hist_nch_transverse_500[0] = bookProfile1D(1+isqrts, 1, 1);
      _hist_nch_toward_500        = bookProfile1D(1+isqrts, 1, 2);
      _hist_nch_away_500          = bookProfile1D(1+isqrts, 1, 3);

      // pTsum profiles, 500 MeV track pT cut
      _hist_ptsum_transverse_500[0] = bookProfile1D(3+isqrts, 1, 1);
      _hist_ptsum_toward_500        = bookProfile1D(3+isqrts, 1, 2);
      _hist_ptsum_away_500          = bookProfile1D(3+isqrts, 1, 3);

      // Standard deviation profiles
      // First the higher moments of main profiles to calculate variance and error on variance...
      for (size_t i = 1; i < 4; ++i) {
        _hist_nch_transverse_500[i].reset(new Profile1D(refData(1+isqrts, 1, 1)));
        _hist_ptsum_transverse_500[i].reset(new Profile1D(refData(3+isqrts, 1, 1)));
      }
      // Then the data point sets into which the results will be inserted
      _dps_sdnch_transverse_500   = bookScatter2D(5+isqrts, 1, 1);
      _dps_sdptsum_transverse_500 = bookScatter2D(7+isqrts, 1, 1);

      // <pT> profiles, 500 MeV track pT cut
      _hist_ptavg_transverse_500 = bookProfile1D(9+isqrts, 1, 1);
      _hist_ptavg_toward_500     = bookProfile1D(9+isqrts, 1, 2);
      _hist_ptavg_away_500       = bookProfile1D(9+isqrts, 1, 3);

      // <pT> vs. Nch profiles, 500 MeV track pT cut
      _hist_dn_dpt_transverse_500 = bookProfile1D(11+isqrts, 1, 1);
      _hist_dn_dpt_toward_500     = bookProfile1D(11+isqrts, 1, 2);
      _hist_dn_dpt_away_500       = bookProfile1D(11+isqrts, 1, 3);

      // Nch vs. Delta(phi) profiles, 500 MeV track pT cut
      _hist_N_vs_dPhi_1_500 = bookProfile1D(13+isqrts, 1, 1);
      _hist_N_vs_dPhi_2_500 = bookProfile1D(13+isqrts, 1, 2);
      _hist_N_vs_dPhi_3_500 = bookProfile1D(13+isqrts, 1, 3);
      _hist_N_vs_dPhi_5_500 = bookProfile1D(13+isqrts, 1, 4);
      // pT vs. Delta(phi) profiles, 500 MeV track pT cut
      _hist_pT_vs_dPhi_1_500 = bookProfile1D(15+isqrts, 1, 1);
      _hist_pT_vs_dPhi_2_500 = bookProfile1D(15+isqrts, 1, 2);
      _hist_pT_vs_dPhi_3_500 = bookProfile1D(15+isqrts, 1, 3);
      _hist_pT_vs_dPhi_5_500 = bookProfile1D(15+isqrts, 1, 4);

      // Nch and pTsum profiles, 100 MeV track pT cut
      _hist_nch_transverse_100   = bookProfile1D(17+isqrts, 1, 1);
      _hist_nch_toward_100       = bookProfile1D(17+isqrts, 1, 2);
      _hist_nch_away_100         = bookProfile1D(17+isqrts, 1, 3);
      _hist_ptsum_transverse_100 = bookProfile1D(19+isqrts, 1, 1);
      _hist_ptsum_toward_100     = bookProfile1D(19+isqrts, 1, 2);
      _hist_ptsum_away_100       = bookProfile1D(19+isqrts, 1, 3);

      // Profiles vs. eta (7 TeV only)
      if (isqrts == 1) {
        _hist_nch_vs_eta_transverse_100   = bookProfile1D(21, 1, 1);
        _hist_ptsum_vs_eta_transverse_100 = bookProfile1D(22, 1, 1);
      }

    }


    // Little helper function to identify Delta(phi) regions
    inline int region_index(double dphi) {
      assert(inRange(dphi, 0.0, PI, CLOSED, CLOSED));
      if (dphi < PI/3.0) return 0;
      if (dphi < 2*PI/3.0) return 1;
      return 2;
    }


    void analyze(const Event& event) {
      const double weight = event.weight();

      // Require at least one track in the event with pT >= 1 GeV
      const ChargedFinalState& cfslead = apply<ChargedFinalState>(event, "CFSlead");
      if (cfslead.size() < 1) {
        vetoEvent;
      }

      // These are the charged particles (tracks) with pT > 500 MeV
      const ChargedFinalState& charged500 = apply<ChargedFinalState>(event, "CFS500");

      // Identify leading track and its phi and pT (this is the same for both the 100 MeV and 500 MeV track cuts)
      Particles particles500 = charged500.particlesByPt();
      Particle p_lead = particles500[0];
      const double philead = p_lead.phi();
      const double etalead = p_lead.eta();
      const double pTlead  = p_lead.perp();
      MSG_DEBUG("Leading track: pT = " << pTlead << ", eta = " << etalead << ", phi = " << philead);

      // Iterate over all > 500 MeV particles and count particles and scalar pTsum in the three regions
      vector<double> num500(3, 0), ptSum500(3, 0.0);
      // Temporary histos that bin Nch and pT in dPhi.
      // NB. Only one of each needed since binnings are the same for the energies and pT cuts
      Histo1D hist_num_dphi_500(refData(13,1,1));
      Histo1D hist_pt_dphi_500(refData(15,1,1));
      foreach (const Particle& p, particles500) {
        const double pT = p.pT();
        const double dPhi = deltaPhi(philead, p.phi());
        const int ir = region_index(dPhi);
        num500[ir] += 1;
        ptSum500[ir] += pT;

        // Fill temp histos to bin Nch and pT in dPhi
        if (p.genParticle() != p_lead.genParticle()) { // We don't want to fill all those zeros from the leading track...
          hist_num_dphi_500.fill(dPhi, 1);
          hist_pt_dphi_500.fill(dPhi, pT);
        }
      }


      // Iterate over charged particles again for profiles against Nch
      // This is necessary since the Nch are region-specific and so are only known after the first loop
      foreach (const Particle& p, particles500) {
        const double pT = p.pT();
        const double dPhi = deltaPhi(philead, p.phi());
        const int ir = region_index(dPhi);
        switch (ir) {
        case 0:
          _hist_dn_dpt_toward_500->fill(num500[0], pT, weight);
          break;
        case 1:
          _hist_dn_dpt_transverse_500->fill(num500[1], pT, weight);
          break;
        case 2:
          _hist_dn_dpt_away_500->fill(num500[2], pT, weight);
          break;
        default:
          assert(false && "How did we get here?");
        }
      }


      // Now fill underlying event histograms
      // The densities are calculated by dividing the UE properties by dEta*dPhi
      // -- each region has a dPhi of 2*PI/3 and dEta is two times 2.5
      const double dEtadPhi = (2*2.5 * 2*PI/3.0);
      // Transverse profiles need 4 orders of moments for stddev with errors
      for (int i = 0; i < 4; ++i) {
        _hist_nch_transverse_500[i]->fill(pTlead/GeV, intpow(num500[1]/dEtadPhi, i+1), weight);
        _hist_ptsum_transverse_500[i]->fill(pTlead/GeV, intpow(ptSum500[1]/GeV/dEtadPhi, i+1), weight);
      }
      // Toward and away profiles only need the first moment (the mean)
      _hist_nch_toward_500->fill(pTlead/GeV, num500[0]/dEtadPhi, weight);
      _hist_nch_away_500->fill(pTlead/GeV, num500[2]/dEtadPhi, weight);
      _hist_ptsum_toward_500->fill(pTlead/GeV, ptSum500[0]/GeV/dEtadPhi, weight);
      _hist_ptsum_away_500->fill(pTlead/GeV, ptSum500[2]/GeV/dEtadPhi, weight);
      // <pT> profiles
      //MSG_INFO("Trans pT1, pTsum, Nch, <pT>" << pTlead/GeV << ", " <<  ptSum500[1]/GeV << ", " << num500[1] << ", " << ptSum500[1]/GeV/num500[1]);
      if (num500[1] > 0) _hist_ptavg_transverse_500->fill(pTlead/GeV, ptSum500[1]/GeV/num500[1], weight);
      if (num500[0] > 0) _hist_ptavg_toward_500->fill(pTlead/GeV, ptSum500[0]/GeV/num500[0], weight);
      if (num500[2] > 0) _hist_ptavg_away_500->fill(pTlead/GeV, ptSum500[2]/GeV/num500[2], weight);


      // Update the "proper" dphi profile histograms
      // Note that we fill dN/dEtadPhi: dEta = 2*2.5, dPhi = 2*PI/nBins
      // The values tabulated in the note are for an (undefined) signed Delta(phi) rather than
      // |Delta(phi)| and so differ by a factor of 2: we have to actually norm for angular range = 2pi
      const size_t nbins = refData(13,1,1).numPoints();
      std::vector<double> ptcut;
      if (fuzzyEquals(sqrtS(), 900*GeV)) {
        ptcut += 1.0; ptcut += 1.5; ptcut += 2.0; ptcut += 2.5;
      }
      else if (fuzzyEquals(sqrtS(), 7*TeV)) {
        ptcut += 1.0; ptcut += 2.0; ptcut += 3.0; ptcut += 5.0;
      }
      assert(ptcut.size() == 4);
      for (size_t i = 0; i < nbins; ++i) {
        // First Nch
        double mean = hist_num_dphi_500.bin(i).xMid();
        double value = 0.;
        if (hist_num_dphi_500.bin(i).numEntries() > 0) {
          mean = hist_num_dphi_500.bin(i).xMean();
          value = hist_num_dphi_500.bin(i).area()/hist_num_dphi_500.bin(i).xWidth()/10.0;
        }
        if (pTlead/GeV >= ptcut[0]) _hist_N_vs_dPhi_1_500->fill(mean, value, weight);
        if (pTlead/GeV >= ptcut[1]) _hist_N_vs_dPhi_2_500->fill(mean, value, weight);
        if (pTlead/GeV >= ptcut[2]) _hist_N_vs_dPhi_3_500->fill(mean, value, weight);
        if (pTlead/GeV >= ptcut[3]) _hist_N_vs_dPhi_5_500->fill(mean, value, weight);

        // Then pT
        mean = hist_pt_dphi_500.bin(i).xMid();
        value = 0.;
        if (hist_pt_dphi_500.bin(i).numEntries() > 0) {
          mean = hist_pt_dphi_500.bin(i).xMean();
          value = hist_pt_dphi_500.bin(i).area()/hist_pt_dphi_500.bin(i).xWidth()/10.0;
        }
        if (pTlead/GeV >= ptcut[0]) _hist_pT_vs_dPhi_1_500->fill(mean, value, weight);
        if (pTlead/GeV >= ptcut[1]) _hist_pT_vs_dPhi_2_500->fill(mean, value, weight);
        if (pTlead/GeV >= ptcut[2]) _hist_pT_vs_dPhi_3_500->fill(mean, value, weight);
        if (pTlead/GeV >= ptcut[3]) _hist_pT_vs_dPhi_5_500->fill(mean, value, weight);
      }


      //////////////////////


      // These are the charged particles (tracks) with pT > 100 MeV
      const ChargedFinalState& charged100 = apply<ChargedFinalState>(event, "CFS100");

      // Iterate over all > 100 MeV particles and count particles and scalar pTsum in the three regions
      vector<double> num100(3, 0), ptSum100(3, 0.0);
      foreach (const Particle& p, charged100.particles()) {
        const double pT = p.pT();
        const double dPhi = deltaPhi(philead, p.phi());
        const int ir = region_index(dPhi);
        num100[ir] += 1;
        ptSum100[ir] += pT;
      }

      // Now fill the two sets of 100 MeV underlying event histograms
      _hist_nch_transverse_100->fill(pTlead/GeV, num100[1]/dEtadPhi, weight);
      _hist_nch_toward_100->fill(pTlead/GeV, num100[0]/dEtadPhi, weight);
      _hist_nch_away_100->fill(pTlead/GeV, num100[2]/dEtadPhi, weight);
      _hist_ptsum_transverse_100->fill(pTlead/GeV, ptSum100[1]/GeV/dEtadPhi, weight);
      _hist_ptsum_toward_100->fill(pTlead/GeV, ptSum100[0]/GeV/dEtadPhi, weight);
      _hist_ptsum_away_100->fill(pTlead/GeV, ptSum100[2]/GeV/dEtadPhi, weight);

      // And finally the Nch and pT vs eta_lead profiles (again from > 100 MeV tracks, and only at 7 TeV)
      if (fuzzyEquals(sqrtS(), 7*TeV) && pTlead > 5*GeV) {
        // MSG_INFO(sqrtS() << " " << pTlead << " " << ptSum100[1]/dEtadPhi << " " << num100[1]/dEtadPhi);
        _hist_nch_vs_eta_transverse_100->fill(etalead, num100[1]/dEtadPhi, weight);
        _hist_ptsum_vs_eta_transverse_100->fill(etalead, ptSum100[1]/GeV/dEtadPhi, weight);
      }

    }


    void finalize() {
      // Convert the various moments of the 500 MeV trans pT and Nch distributions to std devs with correct error
      _moments_to_stddev(_hist_nch_transverse_500, _dps_sdnch_transverse_500);
      _moments_to_stddev(_hist_ptsum_transverse_500, _dps_sdptsum_transverse_500);
    }


  private:


    inline void _moments_to_stddev(Profile1DPtr moment_profiles[], Scatter2DPtr target_dps) {
      for (size_t b = 0; b < moment_profiles[0]->numBins(); ++b) { // loop over points
        /// @todo Assuming unit weights here! Should use N_effective = sumW**2/sumW2?
        const double numentries = moment_profiles[0]->bin(b).numEntries();
        const double x = moment_profiles[0]->bin(b).xMid();
        const double ex = moment_profiles[0]->bin(b).xWidth()/2.;
        double var = 0.;
        double sd = 0.;
        if (numentries > 0) {
          var = moment_profiles[1]->bin(b).mean() - intpow(moment_profiles[0]->bin(b).mean(), 2);
          sd = fuzzyLessEquals(var,0.) ? 0 : sqrt(var); ///< Numerical safety check
        }
        if (sd == 0 || numentries < 3) {
          MSG_WARNING("Need at least 3 bin entries and a non-zero central value to calculate "
                      << "an error on standard deviation profiles (bin " << b << ")");
          target_dps->addPoint(x, sd, ex, 0);
          continue;
        }
        // c2(y) = m4(x) - 4 m3(x) m1(x) - m2(x)^2 + 8 m2(x) m1(x)^2 - 4 m1(x)^4
        const double var_on_var = moment_profiles[3]->bin(b).mean()
          - 4 * moment_profiles[2]->bin(b).mean() * moment_profiles[0]->bin(b).mean()
          - intpow(moment_profiles[1]->bin(b).mean(), 2)
          + 8 * moment_profiles[1]->bin(b).mean() * intpow(moment_profiles[0]->bin(b).mean(), 2)
          - 4 * intpow(moment_profiles[0]->bin(b).mean(), 4);
        const double stderr_on_var = sqrt(var_on_var/(numentries-2.0));
        const double stderr_on_sd = stderr_on_var / (2.0*sd);
        target_dps->addPoint(x, sd, ex, stderr_on_sd);
      }
    }


  private:

    Profile1DPtr _hist_nch_transverse_500[4];
    Profile1DPtr _hist_nch_toward_500;
    Profile1DPtr _hist_nch_away_500;

    Profile1DPtr _hist_ptsum_transverse_500[4];
    Profile1DPtr _hist_ptsum_toward_500;
    Profile1DPtr _hist_ptsum_away_500;

    Scatter2DPtr  _dps_sdnch_transverse_500;
    Scatter2DPtr  _dps_sdptsum_transverse_500;

    Profile1DPtr _hist_ptavg_transverse_500;
    Profile1DPtr _hist_ptavg_toward_500;
    Profile1DPtr _hist_ptavg_away_500;

    Profile1DPtr _hist_dn_dpt_transverse_500;
    Profile1DPtr _hist_dn_dpt_toward_500;
    Profile1DPtr _hist_dn_dpt_away_500;

    Profile1DPtr _hist_N_vs_dPhi_1_500;
    Profile1DPtr _hist_N_vs_dPhi_2_500;
    Profile1DPtr _hist_N_vs_dPhi_3_500;
    Profile1DPtr _hist_N_vs_dPhi_5_500;

    Profile1DPtr _hist_pT_vs_dPhi_1_500;
    Profile1DPtr _hist_pT_vs_dPhi_2_500;
    Profile1DPtr _hist_pT_vs_dPhi_3_500;
    Profile1DPtr _hist_pT_vs_dPhi_5_500;

    Profile1DPtr _hist_nch_transverse_100;
    Profile1DPtr _hist_nch_toward_100;
    Profile1DPtr _hist_nch_away_100;
    Profile1DPtr _hist_ptsum_transverse_100;
    Profile1DPtr _hist_ptsum_toward_100;
    Profile1DPtr _hist_ptsum_away_100;

    Profile1DPtr _hist_nch_vs_eta_transverse_100;
    Profile1DPtr _hist_ptsum_vs_eta_transverse_100;

  };



  // The hook for the plugin system
  DECLARE_RIVET_PLUGIN(ATLAS_2010_S8894728);

}