1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/DileptonFinder.hh"
  5#include "Rivet/Projections/Beam.hh"
  6#include "Rivet/Math/LorentzTrans.hh"
  7#include "Rivet/Math/Vector3.hh"
  8
  9namespace Rivet {
 10
 11
 12  /// @brief Drell-Yan dimuon production in proton-lead collisions at  8.16 TeV
 13  class CMS_2021_I1849180 : public Analysis {
 14  public:
 15
 16    /// Constructor
 17    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2021_I1849180);
 18
 19
 20    /// @name Analysis methods
 21    ///@{
 22
 23    /// Book histograms and initialise projections before the run
 24    void init() {
 25
 26      const ParticlePair& beam = beams();
 27      _pcom = beam.first.momentum() + beam.second.momentum();
 28      if (beam.first.mom().E() == beam.second.mom().E()) _y_shift = -0.465;
 29
 30      DileptonFinder zmumuFinder(91.2*GeV, 0.1, Cuts::abseta < 2.4 && Cuts::pT > 10*GeV &&
 31                                 Cuts::abspid == PID::MUON, Cuts::massIn(15*GeV, 600*GeV));
 32      declare(zmumuFinder, "ZmumuFinder");
 33
 34      DileptonFinder totzmumuFinder(91.2*GeV, 0.1, Cuts::abspid == PID::MUON, Cuts::mass < 1000*GeV);
 35      declare(totzmumuFinder, "TotzmumuFinder");
 36
 37      book(_h["111"], 1, 1, 1);
 38      book(_h["211"], 2, 1, 1);
 39      book(_h["311"], 3, 1, 1);
 40      book(_h["411"], 4, 1, 1);
 41      book(_h["511"], 5, 1, 1);
 42      book(_h["611"], 6, 1, 1);
 43      book(_h["711"], 7, 1, 1);
 44      book(_h["911"], 9, 1, 1);
 45      book(_h["1011"], 10, 1, 1);
 46      book(_h["1111"], 11, 1, 1);
 47      book(_h["1211"], 12, 1, 1);
 48      book(_h["1311"], 13, 1, 1);
 49      book(_h["1411"], 14, 1, 1);
 50      book(_h["1511"], 15, 1, 1);
 51    }
 52
 53    /// Perform the per-event analysis
 54    void analyze(const Event& event) {
 55      const DileptonFinder& ZmumuFinder = apply<DileptonFinder>(event, "ZmumuFinder");
 56      const DileptonFinder& TotzmumuFinder = apply<DileptonFinder>(event, "TotzmumuFinder");
 57      const Particles& zmumus = ZmumuFinder.bosons();
 58      const Particles& totzmumus = TotzmumuFinder.bosons();
 59
 60      const Vector3 betacom = _pcom.betaVec();
 61      const LorentzTransform comboost = LorentzTransform::mkFrameTransformFromBeta(betacom);
 62
 63      if (zmumus.size() == 1) {
 64        for (const Particle& zmumu : zmumus) {
 65          Particles leptons = sortBy(ZmumuFinder.constituents(), cmpMomByPt);
 66          if (leptons[0].pt() < 15 or leptons[1].pt() < 10) {
 67            continue;
 68          }
 69          const FourMomentum lminus = leptons[0].charge() < 0 ? leptons[0].momentum() : leptons[1].momentum();
 70          const FourMomentum lplus = leptons[0].charge() < 0 ? leptons[1].momentum() : leptons[0].momentum();
 71
 72          const double phi_acop = M_PI - deltaPhi(lminus, lplus);
 73          const double costhetastar = tanh((lminus.eta() - lplus.eta()) / 2.0);
 74          const double sin2thetastar = (costhetastar <= 1) ? 1.0 - sqr(costhetastar) : 0;
 75          const double phistar = tan(phi_acop / 2.0) * sqrt(sin2thetastar);
 76
 77          const FourMomentum Zcm = comboost.transform(zmumu.momentum());
 78
 79          const double Zmass = zmumu.mass() / GeV;
 80          const double Zpt = zmumu.momentum().pT() / GeV;
 81
 82          const double Zy = Zcm.rapidity() + _y_shift;
 83
 84          if (Zcm.rapidity() < 1.93 && Zcm.rapidity() > -2.87) {
 85            _h["111"]->fill(Zmass);
 86          }
 87          if (Zmass < 60 && Zmass > 15) {
 88            if (Zcm.rapidity() < 1.93 && Zcm.rapidity() > -2.87) {
 89              _h["411"]->fill(Zpt);
 90              _h["611"]->fill(phistar);
 91            }
 92            _h["211"]->fill(Zy);
 93          }
 94
 95          if (Zmass < 120 && Zmass > 60) {
 96            if (Zcm.rapidity() < 1.93 && Zcm.rapidity() > -2.87) {
 97              _h["511"]->fill(Zpt);
 98              _h["711"]->fill(phistar);
 99            }
100            _h["311"]->fill(Zy);
101          }
102        }
103      }
104
105      if (totzmumus.size() == 1) {
106        for (const Particle& totzmumu : totzmumus) {
107          Particles leptons = sortBy(TotzmumuFinder.constituents(), cmpMomByPt);
108
109          const FourMomentum lminus = leptons[0].charge() < 0 ? leptons[0].momentum() : leptons[1].momentum();
110          const FourMomentum lplus = leptons[0].charge() < 0 ? leptons[1].momentum() : leptons[0].momentum();
111
112          const double phi_acop = M_PI - deltaPhi(lminus, lplus);
113          const double costhetastar = tanh((lminus.eta() - lplus.eta()) / 2.0);
114          const double sin2thetastar = (costhetastar <= 1) ? 1.0 - sqr(costhetastar) : 0;
115          const double totphistar = tan(phi_acop / 2.0) * sqrt(sin2thetastar);
116
117          const FourMomentum totZcm = comboost.transform(totzmumu.momentum());
118
119          const double totZmass = totzmumu.mass() / GeV;
120          const double totZpt = totzmumu.momentum().pT() / GeV;
121          const double totZy = totZcm.rapidity() + _y_shift;
122
123          if (totZcm.rapidity() < 1.93 && totZcm.rapidity() > -2.87) {
124            _h["911"]->fill(totZmass);
125          }
126
127          if (totZmass < 60 && totZmass > 15) {
128            if (totZcm.rapidity() < 1.93 && totZcm.rapidity() > -2.87) {
129              _h["1211"]->fill(totZpt);
130              _h["1411"]->fill(totphistar);
131            }
132            _h["1011"]->fill(totZy);
133          }
134          if (totZmass < 120 && totZmass > 60) {
135            if (totZcm.rapidity() < 1.93 && totZcm.rapidity() > -2.87) {
136              _h["1311"]->fill(totZpt);
137              _h["1511"]->fill(totphistar);
138            }
139            _h["1111"]->fill(totZy);
140          }
141        }
142      }
143    }
144    /// Normalise histograms etc., after the run
145    void finalize() {
146      double norm = crossSection() * 208 / nanobarn / sumW();
147      scale(_h, norm);
148    }
149
150
151  private:
152
153    FourMomentum _pcom;
154
155    map<string, Histo1DPtr> _h;
156
157    /// Offset constant
158    double _y_shift = 0;
159
160  };
161
162
163  RIVET_DECLARE_PLUGIN(CMS_2021_I1849180);
164
165}