Rivet analyses

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

namespace Rivet {


  /// @brief Forward Z boson production in pp collisions at 7 TeV
  class LHCB_2015_I1373300 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(LHCB_2015_I1373300);


    /// @name Analysis methods
    ///@{

    /// Book histograms and initialise projections before the run
    void init() {

      // Initialise and register projections
      Cut cuts = (Cuts::pT >= 20.0*GeV && Cuts::absetaIn(2.0, 4.5) && Cuts::abspid == PID::MUON);
      DileptonFinder zFinder(91.2*GeV, 0.0, cuts, Cuts::massIn(60*GeV, 120*GeV));
      DileptonFinder zFinder_Fsr(91.2*GeV, 0.1, cuts, Cuts::massIn(60*GeV, 120*GeV));
      LeptonFinder wFinder(cuts, 0.0);
      LeptonFinder wFinder_Fsr(cuts, 0.1);
      declare(zFinder, "DileptonFinder");
      declare(zFinder_Fsr, "DileptonFinder_Fsr");
      declare(wFinder, "LeptonFinder");
      declare(wFinder_Fsr, "LeptonFinder_Fsr");

      // Book histograms
      book(_h_zProdRap, 1, 1, 1);
      book(_h_zProdPt, 2, 1, 1);
      book(_h_zProdPhi, 3, 1, 1);
      book(_h_wPlusProdEta, 4, 1, 1);
      book(_h_wMinusProdEta, 4, 1, 3);
      book(_h_zProdRap_Fsr, 1, 1, 3);
      book(_h_zProdPt_Fsr, 2, 1, 3);
      book(_h_zProdPhi_Fsr, 3, 1, 3);
      book(_h_wPlusProdEta_Fsr, 4, 1, 5);
      book(_h_wMinusProdEta_Fsr, 4, 1, 6);
      book(_h_wRatioEta, 5, 1, 1);
      book(_h_wChargeAsymEta, 6, 1, 1);
      book(_h_wRatioEta_Fsr, 5, 1, 2);
      book(_h_wChargeAsymEta_Fsr, 6, 1, 2);
    }


    /// Perform the per-event analysis
    void analyze(const Event& event) {

      // Searching for Z bosons
      const DileptonFinder& z = apply<DileptonFinder>(event, "DileptonFinder");
      if (!z.empty() && z.bosons().size() == 1) {
        const FourMomentum zmom = z.bosons()[0].momentum();
        const Particle& muPlus = z.leptons()[0];
        const Particle& muMinus = z.leptons()[1];
        const double phi = tan(0.5*(M_PI - deltaPhi(muPlus, muMinus))) / cosh(0.5*deltaEta(muPlus, muMinus));
        _h_zProdRap->fill(zmom.absrap());
        _h_zProdPt->fill(zmom.pt()/GeV);
        _h_zProdPhi->fill(phi);
      }

      // Searching for Z bosons + FSR photon
      const DileptonFinder& z_Fsr = apply<DileptonFinder>(event, "DileptonFinder_Fsr");
      if (!z_Fsr.empty() && z_Fsr.bosons().size() == 1) {
        const FourMomentum zmom = z_Fsr.bosons()[0].momentum();
        const Particle& muPlus = z_Fsr.leptons()[0];
        const Particle& muMinus = z_Fsr.leptons()[1];
        const double phi = tan((M_PI - deltaPhi(muPlus, muMinus))/2) / cosh(deltaEta(muPlus, muMinus)/2);
        _h_zProdRap_Fsr->fill(zmom.absrapidity());
        _h_zProdPt_Fsr->fill(zmom.pt()/GeV);
        _h_zProdPhi_Fsr->fill(phi);
      }

      // Searching for W bosons
      const DressedLeptons w = apply<LeptonFinder>(event, "LeptonFinder").dressedLeptons();
      if (w.size() > 0) {
        const Particle& mu = w[0];
        (mu.charge3() > 0 ? _h_wPlusProdEta : _h_wMinusProdEta)->fill(mu.abseta());
      }

      // Searching for W bosons + FSR photon
      const DressedLeptons w_Fsr = apply<LeptonFinder>(event, "LeptonFinder_Fsr").dressedLeptons();
      if (w_Fsr.size() > 0) {
        const Particle& mu = w_Fsr[0];
        (mu.charge3() > 0 ? _h_wPlusProdEta_Fsr : _h_wMinusProdEta_Fsr)->fill(mu.abseta());
      }
    }


    /// Normalise histograms etc., after the run
    void finalize() {

      // Scaling of the histograms and correction on FSR
      const double xs = crossSection()/picobarn;
      const double sf = xs/sumOfWeights()/2.0;

      scale(_h_zProdRap, sf);
      scale(_h_zProdRap_Fsr, sf);
      scale(_h_zProdPt, sf);
      scale(_h_zProdPt_Fsr, sf);
      scale(_h_zProdPhi, sf);
      scale(_h_zProdPhi_Fsr, sf);
      scale(_h_wPlusProdEta, sf);
      scale(_h_wPlusProdEta_Fsr, sf);
      scale(_h_wMinusProdEta, sf);
      scale(_h_wMinusProdEta_Fsr, sf);

      //Calculation of ratio and charge asymmetry
      divide(_h_wPlusProdEta, _h_wMinusProdEta, _h_wRatioEta);
      divide(_h_wPlusProdEta_Fsr, _h_wMinusProdEta_Fsr, _h_wRatioEta_Fsr);
      asymm(_h_wPlusProdEta, _h_wMinusProdEta, _h_wChargeAsymEta);
      asymm(_h_wPlusProdEta_Fsr, _h_wMinusProdEta_Fsr, _h_wChargeAsymEta_Fsr);
      scale(_h_wChargeAsymEta, 100);
      scale(_h_wChargeAsymEta_Fsr, 100);
    }

    ///@}


    /// @name Histograms
    ///@{
    Histo1DPtr _h_zProdRap;
    Histo1DPtr _h_zProdPt;
    Histo1DPtr _h_zProdPhi;
    Histo1DPtr _h_wPlusProdEta;
    Histo1DPtr _h_wMinusProdEta;
    Histo1DPtr _h_zProdRap_Fsr;
    Histo1DPtr _h_zProdPt_Fsr;
    Histo1DPtr _h_zProdPhi_Fsr;
    Histo1DPtr _h_wPlusProdEta_Fsr;
    Histo1DPtr _h_wMinusProdEta_Fsr;
    Estimate1DPtr _h_wRatioEta;
    Estimate1DPtr _h_wChargeAsymEta;
    Estimate1DPtr _h_wRatioEta_Fsr;
    Estimate1DPtr _h_wChargeAsymEta_Fsr;
    ///@}


  };


  RIVET_DECLARE_PLUGIN(LHCB_2015_I1373300);

}