1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/LeadingParticlesFinalState.hh"
  4#include "Rivet/Projections/VetoedFinalState.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief D0 Run I differential W/Z boson cross-section analysis
 10  ///
 11  /// @author Lars Sonnenschein
 12  /// @author Andy Buckley
 13  class D0_2001_I559624 : public Analysis {
 14  public:
 15
 16    RIVET_DEFAULT_ANALYSIS_CTOR(D0_2001_I559624);
 17
 18
 19    /// @name Analysis methods
 20    /// @{
 21
 22    void init() {
 23      // Final state projection
 24      FinalState fs((Cuts::etaIn(-5.0, 5.0))); // corrected for detector acceptance
 25      declare(fs, "FS");
 26
 27      // Z -> e- e+
 28      LeadingParticlesFinalState eeFS(FinalState((Cuts::etaIn(-5.0, 5.0)))); //20.);
 29      eeFS.addParticleIdPair(PID::ELECTRON);
 30      declare(eeFS, "eeFS");
 31
 32      // W- -> e- nu_e~
 33      LeadingParticlesFinalState enuFS(FinalState((Cuts::etaIn(-5.0, 5.0)))); //25.);
 34      enuFS.addParticleId(PID::ELECTRON).addParticleId(PID::NU_EBAR);
 35      declare(enuFS, "enuFS");
 36
 37      // W+ -> e+ nu_e
 38      LeadingParticlesFinalState enubFS(FinalState((Cuts::etaIn(-5.0, 5.0)))); //25.);
 39      enubFS.addParticleId(PID::POSITRON).addParticleId(PID::NU_E);
 40      declare(enubFS, "enubFS");
 41
 42      // Remove neutrinos for isolation of final state particles
 43      VetoedFinalState vfs(fs);
 44      vfs.vetoNeutrinos();
 45      declare(vfs, "VFS");
 46
 47      // Counters
 48      book(_eventsFilledW,"eventsFilledW");
 49      book(_eventsFilledZ,"eventsFilledZ");
 50
 51      // Histograms
 52      book(_h_dsigdpt_w ,1, 1, 1);
 53      book(_h_dsigdpt_z ,1, 1, 2);
 54      book(_h_dsigdpt_scaled_z, 2, 1, 1);
 55    }
 56
 57
 58    void analyze(const Event& event) {
 59      const LeadingParticlesFinalState& eeFS = apply<LeadingParticlesFinalState>(event, "eeFS");
 60      // Z boson analysis
 61      if (eeFS.particles().size() >= 2) {
 62        // If there is a Z candidate:
 63        // Fill Z pT distributions
 64        double deltaM2=1e30,mass2(0.);
 65        double pT=-1.;
 66        const Particles& Zdaughters = eeFS.particles();
 67        for (size_t ix = 0; ix < Zdaughters.size(); ++ix) {
 68          for (size_t iy = ix+1; iy < Zdaughters.size(); ++iy) {
 69            if (Zdaughters[ix].pid()!=-Zdaughters[iy].pid()) continue;
 70            const FourMomentum pmom = Zdaughters[ix].momentum() + Zdaughters[iy].momentum();
 71            double mz2 = pmom.mass2();
 72            double dm2 = fabs(mz2 - sqr(91.118*GeV));
 73            if (dm2 < deltaM2) {
 74              pT = pmom.pT();
 75              deltaM2 = dm2;
 76              mass2 = mz2;
 77            }
 78          }
 79        }
 80        if (pT > 0. && mass2 > 0. && inRange(sqrt(mass2)/GeV, 75.0, 105.0)) {
 81          _eventsFilledZ->fill();
 82          MSG_DEBUG("Z pmom.pT() = " << pT/GeV << " GeV");
 83          _h_dsigdpt_z->fill(pT/GeV);
 84          // return if found a Z
 85          return;
 86        }
 87      }
 88      // There is no Z -> ee candidate... so this might be a W event
 89      const LeadingParticlesFinalState& enuFS = apply<LeadingParticlesFinalState>(event, "enuFS");
 90      const LeadingParticlesFinalState& enubFS = apply<LeadingParticlesFinalState>(event, "enubFS");
 91
 92      double deltaM2=1e30;
 93      double pT=-1.;
 94      for (size_t iw = 0; iw < 2; ++iw) {
 95        Particles Wdaughters;
 96        Wdaughters = (iw == 0) ? enuFS.particles() : enubFS.particles();
 97        for (size_t ix = 0; ix < Wdaughters.size(); ++ix) {
 98          for (size_t iy = ix+1; iy < Wdaughters.size(); ++iy) {
 99            if (Wdaughters[ix].pid() == Wdaughters[iy].pid())  continue;
100            const FourMomentum pmom = Wdaughters[0].momentum() + Wdaughters[1].momentum();
101            double dm2 = abs(pmom.mass2() - sqr(80.4*GeV));
102            if (dm2 < deltaM2) {
103              pT = pmom.pT();
104              deltaM2 = dm2;
105            }
106          }
107        }
108      }
109      if (pT > 0.) {
110        _eventsFilledW->fill();
111        _h_dsigdpt_w->fill(pT/GeV);
112      }
113    }
114
115
116    void finalize() {
117      // Get cross-section per event (i.e. per unit weight) from generator
118      const double xSecPerEvent = crossSectionPerEvent()/picobarn;
119
120      // Correct W pT distribution to W cross-section
121      const double xSecW = xSecPerEvent * dbl(*_eventsFilledW);
122
123      // Correct Z pT distribution to Z cross-section
124      const double xSecZ = xSecPerEvent * dbl(*_eventsFilledZ);
125
126      // Get W and Z pT integrals
127      const double wpt_integral = _h_dsigdpt_w->integral();
128      const double zpt_integral = _h_dsigdpt_z->integral();
129
130      // Divide and scale ratio histos
131      if (xSecW == 0 || wpt_integral == 0 || xSecZ == 0 || zpt_integral == 0) {
132        MSG_WARNING("Not filling ratio plot because input histos are empty");
133      } else {
134        // Scale factor converts event counts to cross-sections, and inverts the
135        // branching ratios since only one decay channel has been analysed for each boson.
136        // Oh, and we put MW/MZ in, like they do in the paper.
137        const double MW_MZ = 0.8820; // Ratio M_W/M_Z
138        const double BRZEE_BRWENU = 0.033632 / 0.1073; // Ratio of branching fractions
139        const double scalefactor = (xSecW / wpt_integral) / (xSecZ / zpt_integral) * MW_MZ * BRZEE_BRWENU;
140        for (size_t ibin = 1; ibin < _h_dsigdpt_w->numBins()+1; ibin++) {
141          double yval(0), yerr(0);
142          if (_h_dsigdpt_w->bin(ibin).sumW() != 0 && _h_dsigdpt_z->bin(ibin).sumW() != 0) {
143            yval = scalefactor * _h_dsigdpt_w->bin(ibin).sumW() / _h_dsigdpt_z->bin(ibin).sumW();
144            yerr = yval * sqrt( sqr(_h_dsigdpt_w->bin(ibin).relErrW()) + sqr(_h_dsigdpt_z->bin(ibin).errW()) );
145          }
146          _h_dsigdpt_scaled_z->bin(ibin).set(yval, yerr);
147        }
148      }
149
150      // Normalize non-ratio histos
151      normalize(_h_dsigdpt_w, xSecW);
152      normalize(_h_dsigdpt_z, xSecZ);
153    }
154
155    /// @}
156
157
158  private:
159
160    /// @name Event counters for cross section normalizations
161    /// @{
162    CounterPtr _eventsFilledW;
163    CounterPtr _eventsFilledZ;
164    /// @}
165
166    /// @{
167    /// Histograms
168    Histo1DPtr  _h_dsigdpt_w;
169    Histo1DPtr  _h_dsigdpt_z;
170    Estimate1DPtr _h_dsigdpt_scaled_z;
171    /// @}
172
173  };
174
175
176
177  RIVET_DECLARE_ALIASED_PLUGIN(D0_2001_I559624, D0_2001_S4674421);
178
179}