Rivet Analyses Reference

  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

// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FinalState.hh"

namespace Rivet {


  /// @brief Calo-based underlying event at 900 GeV and 7 TeV in ATLAS
  ///
  /// @author Jinlong Zhang
  class ATLAS_2011_S8994773 : public Analysis {
  public:

    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2011_S8994773);


    void init() {
      const FinalState fs500((Cuts::etaIn(-2.5, 2.5) && Cuts::pT >=  500*MeV));
      declare(fs500, "FS500");
      const FinalState fslead((Cuts::etaIn(-2.5, 2.5) && Cuts::pT >=  1.0*GeV));
      declare(fslead, "FSlead");

      // Get an index for the beam energy
      isqrts = -1;
      if (isCompatibleWithSqrtS(900)) isqrts = 0;
      else if (isCompatibleWithSqrtS( 7000)) isqrts = 1;
      assert(isqrts >= 0);

      // N profiles, 500 MeV pT cut
      book(_hist_N_transverse_500 ,1+isqrts, 1, 1);
      // pTsum profiles, 500 MeV pT cut
      book(_hist_ptsum_transverse_500 ,3+isqrts, 1, 1);
      // N vs. Delta(phi) profiles, 500 MeV pT cut
      book(_hist_N_vs_dPhi_1_500 ,13+isqrts, 1, 1);
      book(_hist_N_vs_dPhi_2_500 ,13+isqrts, 1, 2);
      book(_hist_N_vs_dPhi_3_500 ,13+isqrts, 1, 3);
    }


    void analyze(const Event& event) {
      // 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));
      for (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);
      _hist_ptsum_transverse_500->fill(pTlead/GeV, ptSum500[1]/GeV/dEtadPhi);

      // 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);
        if (pTlead/GeV >= 2.0) _hist_N_vs_dPhi_2_500->fill(mean, value);
        if (pTlead/GeV >= 3.0) _hist_N_vs_dPhi_3_500->fill(mean, value);
      }

    }


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


    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;

  };



  RIVET_DECLARE_ALIASED_PLUGIN(ATLAS_2011_S8994773, ATLAS_2011_I891834);

}