1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/IdentifiedFinalState.hh"
  5#include "Rivet/Projections/DressedLeptons.hh"
  6
  7namespace Rivet {
  8
  9
 10  class ATLAS_2015_I1351916 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    ATLAS_2015_I1351916(const string name="ATLAS_2015_I1351916", size_t mode=0,
 15                        const string ref_data="ATLAS_2015_I1351916") : Analysis(name) {
 16      _mode = mode; // pick electron channel by default
 17      setRefDataName(ref_data);
 18    }
 19
 20
 21    /// @name Analysis methods
 22    //@{
 23
 24    /// Book histograms and initialise projections before the run
 25    void init() {
 26
 27      const FinalState fs;
 28
 29      IdentifiedFinalState bareleptons(fs);
 30      bareleptons.acceptIdPair(_mode? PID::MUON : PID::ELECTRON);
 31
 32      const Cut cuts = (_mode == 0) ? (Cuts::pT > 25*GeV && Cuts::abseta < 4.9) : (Cuts::pT > 20*GeV && Cuts::abseta < 2.47);
 33      DressedLeptons leptons(fs, bareleptons, 0.1, cuts, true);
 34      declare(leptons, "leptons");
 35
 36
 37      // Book dummy histograms for heterogeneous merging
 38      const Scatter2D& ref = refData(_mode? 4 : 2, 1, 2);
 39      book(_h["NCC_pos"], "_ncc_pos", ref);
 40      book(_h["NCC_neg"], "_ncc_neg", ref);
 41      book(_s["CC"], _mode ? 4 : 2, 1, 2, true);
 42
 43      if (_mode == 0) { // electron-channel only
 44        const Scatter2D& ref_cf = refData(3, 1, 2);
 45        book(_h["NCF_pos"], "_ncf_pos", ref_cf);
 46        book(_h["NCF_neg"], "_ncf_neg", ref_cf);
 47        book(_s["CF"], 3, 1, 2, true);
 48      }
 49    }
 50
 51
 52    /// Perform the per-event analysis
 53    void analyze(const Event& e) {
 54
 55      // Get and cut on dressed leptons
 56      const vector<DressedLepton>& leptons = apply<DressedLeptons>(e, "leptons").dressedLeptons();
 57      if (leptons.size() != 2) vetoEvent; // require exactly two leptons
 58      if (leptons[0].charge3() * leptons[1].charge3() > 0) vetoEvent; // require opposite charge
 59
 60      // Identify lepton vs antilepton
 61      const Particle& lpos = leptons[(leptons[0].charge3() > 0) ? 0 : 1];
 62      const Particle& lneg = leptons[(leptons[0].charge3() < 0) ? 0 : 1];
 63
 64      string label = "N";
 65      if (_mode == 1) {// electron channel
 66        label += "CC"; // only central-central for muons
 67      } else { // electron channel
 68        const double eta1 = lpos.abseta();
 69        const double eta2 = lneg.abseta();
 70        if ( (eta1 < 2.47 && inRange(eta2, 2.5, 4.9)) || (eta2 < 2.47 && inRange(eta1, 2.5, 4.9)) )
 71          label += "CF"; // central-forward
 72        else if (eta1 < 2.47 && eta2 < 2.47)
 73          label += "CC"; // central-central
 74        else vetoEvent; // ain't no forward-forward
 75      }
 76
 77      const double cosThetaStar = cosCollinsSoper(lneg, lpos);
 78      const double mll = (lpos.mom() + lneg.mom()).mass();
 79      label += cosThetaStar < 0.0?  "_neg" : "_pos";
 80      _h[label]->fill(mll/GeV);
 81    }
 82
 83
 84    /// Normalise histograms etc., after the run
 85    void finalize() {
 86      const double sf = crossSectionPerEvent() / picobarn;
 87      for (const auto& key_hist : _h) scale(key_hist.second, sf);
 88      divide(*_h["NCC_pos"] - *_h["NCC_neg"], *_h["NCC_pos"] + *_h["NCC_neg"], _s["CC"]);
 89      if (!_mode)  divide(*_h["NCF_pos"] - *_h["NCF_neg"], *_h["NCF_pos"] + *_h["NCF_neg"], _s["CF"]);
 90    }
 91
 92
 93    // Cosine of the decay angle in the Collins-Soper frame
 94    double cosCollinsSoper(const FourMomentum& l1, const FourMomentum& l2) {
 95      const FourMomentum ll = l1 + l2;
 96      const double nom  = (l1.E() + l1.pz()) * (l2.E() - l2.pz()) - (l1.E() - l1.pz()) * (l2.E() + l2.pz());
 97      const double denom = ll.mass() * sqrt( sqr(ll.mass()) + sqr(ll.pt()) );
 98      return sign(ll.pz()) * safediv(nom, denom); // protect against division by zero, you never know...
 99    }
100
101    //@}
102
103
104  protected:
105
106    /// Electron or muon mode = 0 or 1, for use by derived _EL, _MU analysis classes
107    size_t _mode;
108
109
110  private:
111
112    /// Histograms
113    map<string, Histo1DPtr> _h;
114    /// Asymmetries
115    map<string, Scatter2DPtr> _s;
116
117  };
118
119
120
121  class ATLAS_2015_I1351916_EL : public ATLAS_2015_I1351916 {
122  public:
123    ATLAS_2015_I1351916_EL() : ATLAS_2015_I1351916("ATLAS_2015_I1351916_EL", 0) { }
124  };
125
126
127  class ATLAS_2015_I1351916_MU : public ATLAS_2015_I1351916 {
128  public:
129    ATLAS_2015_I1351916_MU() : ATLAS_2015_I1351916("ATLAS_2015_I1351916_MU", 1) { }
130  };
131
132
133  RIVET_DECLARE_PLUGIN(ATLAS_2015_I1351916);
134  RIVET_DECLARE_PLUGIN(ATLAS_2015_I1351916_EL);
135  RIVET_DECLARE_PLUGIN(ATLAS_2015_I1351916_MU);
136
137}