1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/LeptonFinder.hh"
  5
  6namespace Rivet {
  7
  8
  9  class ATLAS_2015_I1351916 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2015_I1351916);
 14
 15    /// @name Analysis methods
 16    /// @{
 17
 18    /// Book histograms and initialise projections before the run
 19    void init() {
 20
 21      _mode = 0;
 22      if ( getOption("LMODE") == "EL" ) _mode = 1;
 23      else if ( getOption("LMODE") == "MU" ) _mode = 2;
 24
 25      LeptonFinder elecs(0.1, Cuts::pT > 25*GeV && Cuts::abseta < 4.9 && Cuts::abspid == PID::ELECTRON);
 26      declare(elecs, "elecs");
 27
 28      LeptonFinder muons(0.1, Cuts::pT > 20*GeV && Cuts::abseta < 2.47 && Cuts::abspid == PID::MUON);
 29      declare(muons, "muons");
 30
 31      // Book dummy histograms for heterogeneous merging
 32      if (_mode == 0 || _mode == 1) {
 33        const Estimate1D& ref_el = refData(2, 1, 2);
 34        book(_h["el_NCC_pos"], "_el_ncc_pos", ref_el.xEdges());
 35        book(_h["el_NCC_neg"], "_el_ncc_neg", ref_el.xEdges());
 36        book(_s["el_CC"], 2, 1, 2);
 37
 38        // electron-channel only
 39        const Estimate1D& ref_cf = refData(3, 1, 2);
 40        book(_h["el_NCF_pos"], "_el_ncf_pos", ref_cf.xEdges());
 41        book(_h["el_NCF_neg"], "_el_ncf_neg", ref_cf.xEdges());
 42        book(_s["el_CF"], 3, 1, 2);
 43      }
 44      if (_mode == 0 || _mode == 2) {
 45        const Estimate1D& ref_mu = refData(4, 1, 2);
 46        book(_h["mu_NCC_pos"], "_mu_ncc_pos", ref_mu.xEdges());
 47        book(_h["mu_NCC_neg"], "_mu_ncc_neg", ref_mu.xEdges());
 48        book(_s["mu_CC"], 4, 1, 2);
 49      }
 50    }
 51
 52
 53    /// Perform the per-event analysis
 54    void analyze(const Event& event) {
 55
 56      // Get and cut on dressed leptons
 57      const DressedLeptons& elecs = apply<LeptonFinder>(event, "elecs").dressedLeptons();
 58      const DressedLeptons& muons = apply<LeptonFinder>(event, "muons").dressedLeptons();
 59
 60      if (_mode == 0 || _mode == 1)  fillHistos(elecs, true);
 61      if (_mode == 0 || _mode == 2)  fillHistos(muons, false);
 62    }
 63
 64    void fillHistos(const DressedLeptons& leptons, const bool doCF) {
 65
 66      if (leptons.size() != 2)  return; // require exactly two leptons
 67      if (leptons[0].charge3() * leptons[1].charge3() > 0) return; // require opposite charge
 68
 69      // Identify lepton vs antilepton
 70      const Particle& lpos = leptons[(leptons[0].charge3() > 0) ? 0 : 1];
 71      const Particle& lneg = leptons[(leptons[0].charge3() < 0) ? 0 : 1];
 72
 73      string label;
 74      if (doCF) { // electron channel
 75        label = "el_N";
 76        const double eta1 = lpos.abseta();
 77        const double eta2 = lneg.abseta();
 78        if ( (eta1 < 2.47 && inRange(eta2, 2.5, 4.9)) || (eta2 < 2.47 && inRange(eta1, 2.5, 4.9)) ) {
 79          label += "CF"; // central-forward
 80        }
 81        else if (eta1 < 2.47 && eta2 < 2.47) {
 82          label += "CC"; // central-central
 83        }
 84        else  return; // ain't no forward-forward
 85      }
 86      else { // muon channel
 87        label = "mu_NCC"; // only central-central for muons
 88     }
 89
 90      const double cosThetaStar = cosCollinsSoper(lneg, lpos);
 91      const double mll = (lpos.mom() + lneg.mom()).mass();
 92      label += cosThetaStar < 0.0?  "_neg" : "_pos";
 93      _h[label]->fill(mll/GeV);
 94    }
 95
 96
 97    /// Normalise histograms etc., after the run
 98    void finalize() {
 99      const double sf = crossSectionPerEvent() / picobarn;
100      scale(_h, sf);
101      if (_mode == 0 || _mode == 1) {
102        asymm(_h["el_NCC_pos"], _h["el_NCC_neg"], _s["el_CC"]);
103        asymm(_h["el_NCF_pos"], _h["el_NCF_neg"], _s["el_CF"]);
104      }
105      if (_mode == 0 || _mode == 2) {
106        asymm(_h["mu_NCC_pos"], _h["mu_NCC_neg"], _s["mu_CC"]);
107      }
108    }
109
110
111    // Cosine of the decay angle in the Collins-Soper frame
112    double cosCollinsSoper(const FourMomentum& l1, const FourMomentum& l2) {
113      const FourMomentum ll = l1 + l2;
114      const double nom  = (l1.E() + l1.pz()) * (l2.E() - l2.pz()) - (l1.E() - l1.pz()) * (l2.E() + l2.pz());
115      const double denom = ll.mass() * sqrt( sqr(ll.mass()) + sqr(ll.pt()) );
116      return sign(ll.pz()) * safediv(nom, denom); // protect against division by zero, you never know...
117    }
118
119    /// @}
120
121
122  private:
123
124    /// Histograms
125    map<string, Histo1DPtr> _h;
126    /// Asymmetries
127    map<string, Estimate1DPtr> _s;
128
129    // option to steer the lepton channel
130    size_t _mode;
131
132  };
133
134
135  RIVET_DECLARE_PLUGIN(ATLAS_2015_I1351916);
136
137}