rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

ATLAS_2011_S8994773

Calo-based underlying event at 900 GeV and 7 TeV in ATLAS
Experiment: ATLAS (LHC)
Inspire ID: 891834
Status: VALIDATED
Authors:
  • Jinlong Zhang
  • Andy Buckley
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.

Underlying event measurements with the ATLAS detector at the LHC at center-of-mass energies of 900 GeV and 7 TeV, using calorimeter clusters rather than charged tracks.

Source code: ATLAS_2011_S8994773.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
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FinalState.hh"

namespace Rivet {


  /// @author Jinlong Zhang
  class ATLAS_2011_S8994773 : public Analysis {
  public:

    ATLAS_2011_S8994773()
      : Analysis("ATLAS_2011_S8994773") {    }


    void init() {
      const FinalState fs500(-2.5, 2.5, 500*MeV);
      declare(fs500, "FS500");
      const FinalState fslead(-2.5, 2.5, 1.0*GeV);
      declare(fslead, "FSlead");

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

      // N profiles, 500 MeV pT cut
      _hist_N_transverse_500 = bookProfile1D(1+isqrts, 1, 1);
      // pTsum profiles, 500 MeV pT cut
      _hist_ptsum_transverse_500 = bookProfile1D(3+isqrts, 1, 1);
      // N vs. Delta(phi) profiles, 500 MeV 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);
    }


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

      // Require at least one cluster in the event with pT >= 1 GeV
      const FinalState& fslead = apply<FinalState>(event, "FSlead");
      if (fslead.size() < 1) {
        vetoEvent;
      }

      // These are the particles  with pT > 500 MeV
      const FinalState& chargedNeutral500 = apply<FinalState>(event, "FS500");

      // Identify leading object and its phi and pT
      Particles particles500 = chargedNeutral500.particlesByPt();
      Particle p_lead = particles500[0];
      const double philead = p_lead.phi();
      const double etalead = p_lead.eta();
      const double pTlead  = p_lead.pT();
      MSG_DEBUG("Leading object: 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 N 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+isqrts,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 N 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);
        }
      }


      // 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);
      _hist_N_transverse_500->fill(pTlead/GeV,  num500[1]/dEtadPhi, weight);
      _hist_ptsum_transverse_500->fill(pTlead/GeV, ptSum500[1]/GeV/dEtadPhi, 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+isqrts,1,1).numPoints();
      for (size_t i = 0; i < nbins; ++i) {
        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 >= 1.0) _hist_N_vs_dPhi_1_500->fill(mean, value, weight);
        if (pTlead/GeV >= 2.0) _hist_N_vs_dPhi_2_500->fill(mean, value, weight);
        if (pTlead/GeV >= 3.0) _hist_N_vs_dPhi_3_500->fill(mean, value, weight);
      }

    }


    void finalize() {
    }


  private:

    // 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;
    }


  private:
    int isqrts;

    Profile1DPtr _hist_N_transverse_500;

    Profile1DPtr _hist_ptsum_transverse_500;

    Profile1DPtr _hist_N_vs_dPhi_1_500;
    Profile1DPtr _hist_N_vs_dPhi_2_500;
    Profile1DPtr _hist_N_vs_dPhi_3_500;

  };



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

}