D0_2001_S4674421.cc

Go to the documentation of this file.

// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/RivetYODA.hh"
#include "Rivet/Tools/Logging.hh"
#include "Rivet/Tools/ParticleIdUtils.hh"
#include "Rivet/Projections/LeadingParticlesFinalState.hh"
#include "Rivet/Projections/VetoedFinalState.hh"

namespace Rivet {


  /// @brief D0 Run I differential W/Z boson cross-section analysis
  /// @author Lars Sonnenschein
  class D0_2001_S4674421 : public Analysis {
  public:

    /// @name Constructors etc.
    //@{

    /// Constructor.
    D0_2001_S4674421()
      : Analysis("D0_2001_S4674421")
    {    }


    /// @name Analysis methods
    //@{

    void init() {
      // Final state projection
      FinalState fs(-5.0, 5.0); // corrected for detector acceptance
      addProjection(fs, "FS");

      // Z -> e- e+
      LeadingParticlesFinalState eeFS(FinalState(-5.0, 5.0, 0.)); //20.);
      eeFS.addParticleIdPair(ELECTRON);
      addProjection(eeFS, "eeFS");

      // W- -> e- nu_e~
      LeadingParticlesFinalState enuFS(FinalState(-5.0, 5.0, 0.)); //25.);
      enuFS.addParticleId(ELECTRON).addParticleId(NU_EBAR);
      addProjection(enuFS, "enuFS");

      // W+ -> e+ nu_e
      LeadingParticlesFinalState enubFS(FinalState(-5.0, 5.0, 0.)); //25.);
      enubFS.addParticleId(POSITRON).addParticleId(NU_E);
      addProjection(enubFS, "enubFS");

      // Remove neutrinos for isolation of final state particles
      VetoedFinalState vfs(fs);
      vfs.vetoNeutrinos();
      addProjection(vfs, "VFS");

      // Counters
      _eventsFilledW = 0.0;
      _eventsFilledZ = 0.0;

      // Histograms
      _h_dsigdpt_w = bookHisto1D(1, 1, 1);
      _h_dsigdpt_z = bookHisto1D(1, 1, 2);
      _h_dsigdpt_scaled_z = bookScatter2D(2, 1, 1);
    }



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

      const LeadingParticlesFinalState& eeFS = applyProjection<LeadingParticlesFinalState>(event, "eeFS");
      // Z boson analysis
      if (eeFS.particles().size() >= 2) {
        // If there is a Z candidate:
        // Fill Z pT distributions
        double deltaM2=1e30,mass2(0.);
        double pT=-1.;
        const ParticleVector& Zdaughters = eeFS.particles();
        for(unsigned int ix=0;ix<Zdaughters.size();++ix) {
          for(unsigned int iy=ix+1;iy<Zdaughters.size();++iy) {
            if(Zdaughters[ix].pdgId()!=-Zdaughters[iy].pdgId()) continue;
            const FourMomentum pmom = Zdaughters[ix].momentum() + Zdaughters[iy].momentum();
            double mz2 = pmom.mass2();
            double dm2 = abs(mz2-sqr(91.118*GeV));
            if(dm2<deltaM2) {
              pT = pmom.pT();
              deltaM2 = dm2;
              mass2 = mz2;
            }
          }
        }
        if (pT>0. && mass2 > 0. && inRange(sqrt(mass2)/GeV, 75.0, 105.0)) {
          _eventsFilledZ += weight;
          MSG_DEBUG("Z pmom.pT() = " << pT/GeV << " GeV");
          _h_dsigdpt_z->fill(pT/GeV, weight);
          // return if found a Z
          return;
        }
      }
      // There is no Z -> ee candidate... so this might be a W event
      const LeadingParticlesFinalState& enuFS = applyProjection<LeadingParticlesFinalState>(event, "enuFS");
      const LeadingParticlesFinalState& enubFS = applyProjection<LeadingParticlesFinalState>(event, "enubFS");

      double deltaM2=1e30;
      double pT=-1.;
      for(unsigned int iw=0;iw<2;++iw) {
        ParticleVector Wdaughters;
        Wdaughters = iw==0 ? enuFS.particles() : enubFS.particles();
        for(unsigned int ix=0;ix<Wdaughters.size();++ix) {
          for(unsigned int iy=ix+1;iy<Wdaughters.size();++iy) {
            if(Wdaughters[ix].pdgId()==Wdaughters[iy].pdgId())  continue;
            const FourMomentum pmom = Wdaughters[0].momentum() + Wdaughters[1].momentum();
            double dm2 = abs(pmom.mass2()-sqr(80.4*GeV));
            if(dm2<deltaM2) {
              pT = pmom.pT();
              deltaM2 = dm2;
            }
          }
        }
      }
      if(pT>0.) {
        _eventsFilledW += weight;
        _h_dsigdpt_w->fill(pT/GeV, weight);
      }
    }



    void finalize() {
      // Get cross-section per event (i.e. per unit weight) from generator
      const double xSecPerEvent = crossSectionPerEvent()/picobarn;

      // Correct W pT distribution to W cross-section
      const double xSecW = xSecPerEvent * _eventsFilledW;

      // Correct Z pT distribution to Z cross-section
      const double xSecZ = xSecPerEvent * _eventsFilledZ;

      // Get W and Z pT integrals
      const double wpt_integral = integral(_h_dsigdpt_w);
      const double zpt_integral = integral(_h_dsigdpt_z);

      // Divide and scale ratio histos
      if (xSecW == 0 || wpt_integral == 0 || xSecZ == 0 || zpt_integral == 0) {
        MSG_WARNING("Not filling ratio plot because input histos are empty");
      } else {
        // Scale factor converts event counts to cross-sections, and inverts the
        // branching ratios since only one decay channel has been analysed for each boson.
        // Oh, and we put MW/MZ in, like they do in the paper.
        const double MW_MZ = 0.8820; // Ratio M_W/M_Z
        const double BRZEE_BRWENU = 0.033632 / 0.1073; // Ratio of branching fractions
        const double scalefactor = (xSecW / wpt_integral) / (xSecZ / zpt_integral) * MW_MZ * BRZEE_BRWENU;
        for (size_t ibin=0; ibin<_h_dsigdpt_scaled_z->numPoints(); ibin++) {
          if (_h_dsigdpt_w->bin(ibin).area() == 0 || _h_dsigdpt_z->bin(ibin).area() == 0) {
            _h_dsigdpt_scaled_z->point(ibin) = Point2D(_h_dsigdpt_w->bin(ibin).midpoint(), 0.,
                                                       _h_dsigdpt_w->bin(ibin).width(), 0.);
          } else {
            double yval = scalefactor * _h_dsigdpt_w->bin(ibin).area() / _h_dsigdpt_z->bin(ibin).area();
            double dy2 = 0.;
            // binWidth(ibin) is needed because binHeight is actually sumofweights. It's AIDA. Don't ask.  :-((((
            dy2 += pow(_h_dsigdpt_w->bin(ibin).areaErr()/_h_dsigdpt_w->bin(ibin).height(),2);
            dy2 += pow(_h_dsigdpt_z->bin(ibin).areaErr()/_h_dsigdpt_z->bin(ibin).height(),2);
            double dy = scalefactor * _h_dsigdpt_w->bin(ibin).area()/_h_dsigdpt_z->bin(ibin).area() * sqrt(dy2);

            _h_dsigdpt_scaled_z->point(ibin) = Point2D(_h_dsigdpt_w->bin(ibin).midpoint(), yval,
                                                       _h_dsigdpt_w->bin(ibin).width(), dy);
          }
        }
      }

      // Normalize non-ratio histos
      normalize(_h_dsigdpt_w, xSecW);
      normalize(_h_dsigdpt_z, xSecZ);

    }


    //@}

  private:

    /// @name Event counters for cross section normalizations
    //@{
    double _eventsFilledW;
    double _eventsFilledZ;
    //@}

    //@{
    /// Histograms
    Histo1DPtr  _h_dsigdpt_w;
    Histo1DPtr  _h_dsigdpt_z;
    Scatter2DPtr _h_dsigdpt_scaled_z;
    //@}

  };



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

}