Rivet analyses

// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/DISFinalState.hh"
#include "Rivet/Projections/DISKinematics.hh"
#include "Rivet/Projections/DISLepton.hh"
#include "Rivet/Projections/FastJets.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/UnstableParticles.hh"

namespace Rivet {


  /// @brief Measurement of inclusive production of D* mesons both with and without dijet production in DIS collisions (H1)
  class H1_2007_I736052 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(H1_2007_I736052);


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

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

      declare(DISLepton(), "Lepton");
      declare(DISKinematics(), "Kinematics");

      const FinalState fs;
      declare(fs, "FS");
      const UnstableParticles ufs;
      declare(ufs, "UFS");

      double jet_radius = 1.0;
      const DISFinalState DISfs(DISFrame::BREIT);
      declare(FastJets(DISfs, fastjet::JetAlgorithm::kt_algorithm, fastjet::RecombinationScheme::Et_scheme,
                       jet_radius),
              "jets_fs");

      Histo1DPtr tmp;
      book(_h["411"], 4, 1, 1);
      book(_h["511"], 5, 1, 1);
      book(_h["611"], 6, 1, 1);
      book(_h["711"], 7, 1, 1);
      book(_h["811"], 8, 1, 1);
      book(_h["911"], 9, 1, 1);
      book(_h_binned["Q2xbj"], {2., 4.22, 10., 17.8, 31.6, 100.},
           {"d10-x01-y01", "d11-x01-y01", "d12-x01-y01", "d13-x01-y01", "d14-x01-y01"});
      book(_h["1511"], 15, 1, 1);
      book(_h["1611"], 16, 1, 1);
      book(_h["1711"], 17, 1, 1);
      book(_h["1811"], 18, 1, 1);
      book(_h["1911"], 19, 1, 1);
      book(_h["2011"], 20, 1, 1);
      book(_h["2111"], 21, 1, 1);
      book(_h["2211"], 22, 1, 1);

      book(_h_binned["Q2phi"], {2., 10., 100.}, {"d23-x01-y01", "d24-x01-y01"});
      book(_h["2511"], 25, 1, 1);
      book(_h["2611"], 26, 1, 1);
      book(_h["2711"], 27, 1, 1);
      book(_h["2811"], 28, 1, 1);
      book(_h_binned["Q2xgam"], {2., 5., 10., 100.}, {"d29-x01-y01", "d29-x01-y02", "d29-x01-y03"});
      book(_h["3011"], 30, 1, 1);
      book(_h_binned["Q2xglue"], {2., 5., 10., 100.}, {"d31-x01-y01", "d31-x01-y02", "d31-x01-y03"});
    }


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

      const FinalState& fs = apply<FinalState>(event, "FS");
      const size_t numParticles = fs.particles().size();

      Jets jets_fs = apply<JetFinder>(event, "jets_fs").jetsByPt(); // Jets with cut on eta
      double jet_radius = 1.0;

      const UnstableParticles& ufs = apply<UnstableFinalState>(event, "UFS");

      if (numParticles < 2) {
        MSG_DEBUG("Failed leptonic event cut");
        vetoEvent;
      }

      Particles Dstar;
      for (const Particle& p : select(ufs.particles(),
                                      Cuts::pT > 1.5 * GeV and Cuts::pT < 15 * GeV and Cuts::abseta < 1.5
                                          and Cuts::abspid == 413)) {
        Dstar.push_back(p);
      }
      if (Dstar.size() == 0) { // Cut on Dstar
        MSG_DEBUG("Failed Dstar cut");
        vetoEvent;
      }

      const DISKinematics& dk = apply<DISKinematics>(event, "Kinematics");
      const DISLepton& dl = apply<DISLepton>(event, "Lepton");


      double Q2 = dk.Q2();
      double y = dk.y();

      if (y < 0.05 or y > 0.7 or Q2 < 2 or Q2 > 100) { // Cut on event kinematics
        MSG_DEBUG("Failed kinematics cut");
        vetoEvent;
      }

      // Extract the particles other than the lepton
      Particles particles;
      particles.reserve(fs.particles().size());
      ConstGenParticlePtr dislepGP = dl.out().genParticle();
      for (const Particle& p : fs.particles()) {
        ConstGenParticlePtr loopGP = p.genParticle();
        if (loopGP == dislepGP) continue;
        particles.push_back(p);
      }


      const LorentzTransform hcmboost = dk.boostHCM();        // Hadron cm system
      const LorentzTransform breitboost = dk.boostBreit();    //Breit system
      const LorentzTransform labboost = breitboost.inverse(); //Labsystem from Breit

      double xbj = dk.x();
      double W2 = dk.W2();
      double pT_cm(0);

      _h["411"]->fill(Q2);
      _h["511"]->fill(xbj);
      _h["611"]->fill(std::sqrt(W2));

      _h_binned["Q2xbj"]->fill(Q2, xbj);

      for (const Particle& p : Dstar) {
        _h["711"]->fill(p.pT());
        _h["811"]->fill(p.eta());
        _h["911"]->fill((p.E() - p.pz()) / (2 * y * dk.beamLepton().E()));

        const FourMomentum hcmMom = hcmboost.transform(p.momentum());
        if (pT_cm < hcmMom.pT()) pT_cm = hcmMom.pT();

        if (hcmMom.pT() > 2) { //Cut for 1711 and 1811 histograms
          _h["1711"]->fill(p.pT());
          _h["1811"]->fill(p.eta());
        }
      }

      if (pT_cm > 2) {
        _h["1511"]->fill(Q2);
        _h["1611"]->fill(xbj);
      }

      bool two_jets_with_cut(false);

      Jets jet_cut;

      for (const Jet& j : jets_fs) {
        double etalab = labboost.transform(j.momentum()).eta();
        if (j.momentum().Et() > 3 and etalab > -1 and etalab < 2.5) { // Cut on jet energy in Breit sys.
          jet_cut.push_back(j);
        }
      }

      if (jet_cut.size() >= 2 && jet_cut[0].momentum().Et() > 4) two_jets_with_cut = true;

      double delta_phi;
      FourMomentum jet1;
      FourMomentum jet2;

      Jets jet_DJ, jet_OJ;
      bool found_DJ(false), found_OJ(false);
      //cout << " check D and other jets " << endl;
      if (two_jets_with_cut) {
        jet1 = jet_cut[0].momentum(); // momentum of jet #1 in Breit sys.
        jet2 = jet_cut[1].momentum(); // momentum of jet #2 in Breit sys.

        delta_phi = deltaPhi(jet1, jet2);
        _h["1911"]->fill(Q2);
        _h["2011"]->fill(xbj);
        _h["2111"]->fill(jet_cut[0].momentum().Et() / GeV);
        _h["2211"]->fill(FourMomentum(jet1 + jet2).mass() / GeV); // Jets invariant mass in Breit sys.
        _h_binned["Q2phi"]->fill(Q2, delta_phi);
        for (const Jet& jet : jet_cut) {
          for (const Particle& p : Dstar) {
            if (deltaR(breitboost.transform(p.momentum()), jet.momentum()) < jet_radius) {
              jet_DJ.push_back(jet);
              found_DJ = true;
            }
            else {
              jet_OJ.push_back(jet);
              found_OJ = true;
            }
          }

          if (found_DJ && found_OJ) break;
        }
        if (jet_DJ.size() > 0 && jet_OJ.size() > 0) {
          double eta_DJ_breit = jet_DJ[0].momentum().eta();
          double eta_OJ_breit = jet_OJ[0].momentum().eta();

          _h["2511"]->fill(eta_DJ_breit);
          _h["2611"]->fill(eta_OJ_breit);
          _h["2711"]->fill(abs(eta_DJ_breit - eta_OJ_breit));

          double x_gamma, x_gluon;
          double E_star_p_z_star(0); // E() - pz() sum for x_gamma in \gamma p cm
          double E_star_p_z_had(0);

          // for jets in had cms: boost first back to lab and then to hcm
          Jet jet_DJ_hcm, jet_OJ_hcm;
          jet_DJ_hcm = hcmboost.transform(labboost(jet_DJ[0].momentum()));
          jet_OJ_hcm = hcmboost.transform(labboost(jet_OJ[0].momentum()));

          // Need to change sign: by default hcmboost has gamma* in +z dir,
          // and p in -z dir, but here we need: gamma* -in -z and proton in +z.
          // so - -> +: watch out for E-pz -> E+pz  and exp(eta_OJ_hcm) -> exp(-eta_OJ_hcm)
          // this was already noted in  H1_2007_I746380

          double Et_DJ_hcm = jet_DJ_hcm.Et();
          double eta_DJ_hcm = jet_DJ_hcm.eta();
          double Et_OJ_hcm = jet_OJ_hcm.Et();
          double eta_OJ_hcm = jet_OJ_hcm.eta();

          //observed fraction of the photon momentum carried by the parton involved in the hard subprocess

          E_star_p_z_star = jet_DJ_hcm.momentum().E() + jet_DJ_hcm.momentum().pz() + jet_OJ_hcm.momentum().E()
              + jet_OJ_hcm.momentum().pz();

          for (size_t ip1 = 0; ip1 < particles.size(); ++ip1) {
            const Particle& p = particles[ip1];
            E_star_p_z_had += (hcmboost.transform(p.momentum())).E()
                + (hcmboost.transform(p.momentum())).pz();
          }
          x_gamma = E_star_p_z_star / E_star_p_z_had;
          x_gluon = (Et_OJ_hcm * exp(-eta_OJ_hcm) + Et_DJ_hcm * exp(-eta_DJ_hcm))
              / (2. * hcmboost.transform(dk.beamHadron().momentum()).E());
          //observed fraction of the proton momentum carried by the gluon

          _h["2811"]->fill(x_gamma);
          _h_binned["Q2xgam"]->fill(Q2, x_gamma);

          _h["3011"]->fill(log10(x_gluon));
          _h_binned["Q2xglue"]->fill(Q2, log10(x_gluon));
        }
      }
    }


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

      const double norm = crossSection() / nanobarn / sumW();
      scale(_h, norm);
      scale(_h_binned, norm);
      divByGroupWidth(_h_binned);

      // new scaling needed, since x bins are in log10(x)
      for (auto& b : _h["3011"]->bins()) {
        const double scale_new = b.xWidth() / pow(10, b.xWidth());
        // jacobian d log10/dx = dlog10 / dlogx dlogx/dx = 2.3026 /x
        const double factor = pow(10, b.xMid()) / 2.3026;
        b.scaleW(scale_new / factor);
      }

      // new scaling needed, since x bins are in log10(x)
      for (auto& histo : _h_binned["Q2xglue"]->bins()) {
        for (auto& b : histo->bins()) {
          const double scale_new = b.xWidth() / pow(10, b.xWidth());
          // jacobian d log10/dx = dlog10 / dlogx dlogx/dx = 2.3026 /x
          const double factor = pow(10, b.xMid()) / 2.3026;
          b.scaleW(scale_new / factor);
        }
      }
    }

    ///@}


    /// @name Histograms
    ///@{
    map<string, Histo1DPtr> _h;
    map<string, Histo1DGroupPtr> _h_binned;
    ///@}
  };


  RIVET_DECLARE_PLUGIN(H1_2007_I736052);

}