1#include "Rivet/Analysis.hh"
  2#include "Rivet/Projections/Beam.hh"
  3#include "Rivet/Projections/Thrust.hh"
  4#include "Rivet/Projections/DileptonFinder.hh"
  5#include "Rivet/Projections/FParameter.hh"
  6#include "Rivet/Projections/Spherocity.hh"
  7#include "Rivet/Projections/ChargedFinalState.hh"
  8#include "Rivet/Projections/VetoedFinalState.hh"
  9
 10namespace Rivet {
 11
 12
 13  /// @brief Event shapes in leptonic $Z$-events
 14  class ATLAS_2016_I1424838 : public Analysis {
 15  public:
 16
 17    /// Constructor
 18    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2016_I1424838);
 19
 20
 21    /// Book histograms and initialise projections before the run
 22    void init() {
 23
 24      // Charged particles inside acceptance region
 25      const ChargedFinalState cfs(Cuts::abseta < 2.5 && Cuts::pT > 500*MeV);
 26      declare(cfs, "CFS");
 27
 28      // Dilepton finders
 29      DileptonFinder zfinder(91.2*GeV, 0.1, Cuts::abseta < 2.4 && Cuts::pT > 20.0*GeV && Cuts::abspid == PID::ELECTRON, Cuts::massIn(66*GeV, 116*GeV));
 30      declare(zfinder, "DileptonFinder");
 31      DileptonFinder zfinder_mu(91.2*GeV, 0.1, Cuts::abseta < 2.4 && Cuts::pT > 20.0*GeV && Cuts::abspid == PID::MUON, Cuts::massIn(66*GeV, 116*GeV));
 32      declare(zfinder_mu, "DileptonFinderMu");
 33
 34      // This CFS only contains charged particles inside the acceptance excluding the leptons
 35      VetoedFinalState remfs(cfs);
 36      remfs.addVetoOnThisFinalState(zfinder);
 37      remfs.addVetoOnThisFinalState(zfinder_mu);
 38      declare(remfs, "REMFS");
 39
 40      const FParameter fparam(remfs);
 41      declare(fparam, "FParameter_");
 42
 43      const Spherocity sphero(remfs);
 44      declare(sphero, "Spherocity_");
 45
 46
 47      // Booking of ES histos
 48      for (size_t alg = 0; alg < 5; ++alg) {
 49        // Book the inclusive histograms
 50        size_t offset = alg * 6;
 51        book(_h_Elec_Ntrk[alg]         , offset+1,  1, 1);
 52        book(_h_Elec_SumPt[alg]        , offset+2,  1, 1);
 53        book(_h_Elec_Beamthrust[alg]   , offset+3,  1, 1);
 54        book(_h_Elec_Thrust[alg]       , offset+4,  1, 1);
 55        book(_h_Elec_FParam[alg]       , offset+5,  1, 1);
 56        book(_h_Elec_Spherocity[alg]   , offset+6,  1, 1);
 57        book(_h_Muon_Ntrk[alg]         , offset+31, 1, 1);
 58        book(_h_Muon_SumPt[alg]        , offset+32, 1, 1);
 59        book(_h_Muon_Beamthrust[alg]   , offset+33, 1, 1);
 60        book(_h_Muon_Thrust[alg]       , offset+34, 1, 1);
 61        book(_h_Muon_FParam[alg]       , offset+35, 1, 1);
 62        book(_h_Muon_Spherocity[alg]   , offset+36, 1, 1);
 63      }
 64    }
 65
 66
 67    /// Perform the per-event analysis
 68    void analyze(const Event& event) {
 69
 70      // Check for Z boson in event
 71      const DileptonFinder& zfinder    = apply<DileptonFinder>(event, "DileptonFinder");
 72      MSG_DEBUG("Num e+ e- pairs found = " << zfinder.bosons().size());
 73      const bool isElec = zfinder.bosons().size() == 1;
 74
 75      const DileptonFinder& zfinder_mu = apply<DileptonFinder>(event, "DileptonFinderMu");
 76      MSG_DEBUG("Num mu+ mu- pairs found = " << zfinder_mu.bosons().size());
 77      const bool isMuon = zfinder_mu.bosons().size() == 1;
 78
 79      // Only accept events with exactly two electrons or exactly two muons
 80      if (isElec && isMuon) vetoEvent;
 81      if (!(isElec || isMuon)) vetoEvent;
 82
 83      // This determines the Zpt phase-space
 84      double zpT = -1000;
 85      if (isElec) zpT = zfinder.bosons()[0].pT();
 86      if (isMuon) zpT = zfinder_mu.bosons()[0].pT();
 87
 88      unsigned int alg = 4; //< for > 25 GeV
 89      if (zpT < 6*GeV) alg = 1;
 90      else if (inRange(zpT/GeV, 6, 12)) alg = 2;
 91      else if (inRange(zpT/GeV, 12, 25)) alg = 3;
 92      assert(alg < 5);
 93      assert(alg > 0);
 94
 95      // All charged particles within |eta|<2.5 except the leptons from Z-decay
 96      const VetoedFinalState& remfs = apply<VetoedFinalState>(event, "REMFS");
 97      // sumPt and Beamthrust (the latter will only be filled if the min Nch criterion is met)
 98      // and Thrust preparation
 99      double sumPt = 0.0, beamThrust = 0.0;
100      vector<Vector3> momenta;
101      for (const Particle& p : remfs.particles()) {
102        const double pT = p.pT();
103        sumPt += pT;
104        beamThrust += pT*exp(-p.abseta());
105        const Vector3 mom = p.mom().pTvec();
106        momenta.push_back(mom);
107      }
108
109      // Fill inclusive histos
110      if (isElec) {
111        _h_Elec_Ntrk[alg]       ->fill(remfs.size());
112        _h_Elec_Ntrk[0]         ->fill(remfs.size());
113        _h_Elec_SumPt[alg]      ->fill(sumPt);
114        _h_Elec_SumPt[0]        ->fill(sumPt);
115      }
116      if (isMuon) {
117        _h_Muon_Ntrk[alg]       ->fill(remfs.size());
118        _h_Muon_Ntrk[0]         ->fill(remfs.size());
119        _h_Muon_SumPt[alg]      ->fill(sumPt);
120        _h_Muon_SumPt[0]        ->fill(sumPt);
121      }
122
123      // Skip event shape calculation if we don't match the minimum Nch criterion
124      if (remfs.size() >=2) {
125
126        // Eventshape calculations
127
128        // Calculate transverse Thrust using all charged FS particles except the lepton
129        // This is copied/inspired from the CMS_6000011_S8957746 analysis
130        if (momenta.size() == 2) {
131          // We need to use a ghost so that Thrust.calc() doesn't return 1.
132          momenta.push_back(Vector3(1e-10*MeV, 0., 0.));
133        }
134        Thrust thrustC;
135        thrustC.calc(momenta);
136
137        double thrust = thrustC.thrust();
138
139        // F-Parameter
140        const FParameter& fparam = apply<FParameter>(event, "FParameter_");
141        // Spherocity
142        const Spherocity& sphero = apply<Spherocity>(event, "Spherocity_");
143
144        // Histos differential in NMPI
145
146        // Fill inclusive histos
147        if (isElec) {
148          _h_Elec_Thrust[alg]     ->fill(thrust);
149          _h_Elec_Thrust[0]       ->fill(thrust);
150          _h_Elec_FParam[alg]     ->fill(fparam.F());
151          _h_Elec_FParam[0]       ->fill(fparam.F());
152          _h_Elec_Spherocity[alg] ->fill(sphero.spherocity());
153          _h_Elec_Spherocity[0]   ->fill(sphero.spherocity());
154          _h_Elec_Beamthrust[alg] ->fill(beamThrust/GeV);
155          _h_Elec_Beamthrust[0]   ->fill(beamThrust/GeV);
156        }
157        if (isMuon) {
158          _h_Muon_Thrust[alg]     ->fill(thrust);
159          _h_Muon_Thrust[0]       ->fill(thrust);
160          _h_Muon_FParam[alg]     ->fill(fparam.F());
161          _h_Muon_FParam[0]       ->fill(fparam.F());
162          _h_Muon_Spherocity[alg] ->fill(sphero.spherocity());
163          _h_Muon_Spherocity[0]   ->fill(sphero.spherocity());
164          _h_Muon_Beamthrust[alg] ->fill(beamThrust/GeV);
165          _h_Muon_Beamthrust[0]   ->fill(beamThrust/GeV);
166        }
167      }
168    }
169
170
171    /// Normalise histograms etc., after the run
172    void finalize() {
173      for (size_t alg = 0; alg < 5; ++alg) {
174        normalize(_h_Elec_Ntrk[alg],100.);
175        normalize(_h_Elec_SumPt[alg],100.);
176        normalize(_h_Elec_Beamthrust[alg],100.);
177        normalize(_h_Elec_Thrust[alg]);
178        normalize(_h_Elec_FParam[alg]);
179        normalize(_h_Elec_Spherocity[alg]);
180        normalize(_h_Muon_Ntrk[alg],100.);
181        normalize(_h_Muon_SumPt[alg],100.);
182        normalize(_h_Muon_Beamthrust[alg],100.);
183        normalize(_h_Muon_Thrust[alg]);
184        normalize(_h_Muon_FParam[alg]);
185        normalize(_h_Muon_Spherocity[alg]);
186      }
187    }
188
189
190  private:
191
192    Histo1DPtr _h_Elec_Ntrk[5];
193    Histo1DPtr _h_Elec_SumPt[5];
194    Histo1DPtr _h_Elec_Beamthrust[5];
195    Histo1DPtr _h_Elec_Thrust[5];
196    Histo1DPtr _h_Elec_FParam[5];
197    Histo1DPtr _h_Elec_Spherocity[5];
198
199    Histo1DPtr _h_Muon_Ntrk[5];
200    Histo1DPtr _h_Muon_SumPt[5];
201    Histo1DPtr _h_Muon_Beamthrust[5];
202    Histo1DPtr _h_Muon_Thrust[5];
203    Histo1DPtr _h_Muon_FParam[5];
204    Histo1DPtr _h_Muon_Spherocity[5];
205
206  };
207
208
209  RIVET_DECLARE_PLUGIN(ATLAS_2016_I1424838);
210
211
212}