1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/PromptFinalState.hh"
  5#include "Rivet/Projections/IdentifiedFinalState.hh"
  6#include "Rivet/Projections/LeptonFinder.hh"
  7#include "Rivet/Projections/VetoedFinalState.hh"
  8#include "Rivet/Projections/FastJets.hh"
  9
 10namespace Rivet {
 11
 12
 13  /// @brief Colour flow in hadronic top decay at 8 TeV
 14  class ATLAS_2015_I1376945 : public Analysis {
 15  public:
 16
 17    /// Constructor
 18    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2015_I1376945);
 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      PromptFinalState promptFs(fs);
 30      promptFs.acceptTauDecays(true);
 31      promptFs.acceptMuonDecays(false);
 32
 33      IdentifiedFinalState neutrino_fs(promptFs);
 34      neutrino_fs.acceptNeutrinos();
 35      declare(neutrino_fs, "NEUTRINO_FS");
 36
 37      IdentifiedFinalState Photon(fs);
 38      Photon.acceptIdPair(PID::PHOTON);
 39
 40      IdentifiedFinalState bare_muons_fs(promptFs);
 41      bare_muons_fs.acceptIdPair(PID::MUON);
 42
 43      IdentifiedFinalState bare_elecs_fs(promptFs);
 44      bare_elecs_fs.acceptIdPair(PID::ELECTRON);
 45
 46      Cut lep_cuts = (Cuts::abseta < 2.5) && (Cuts::pT > 1*MeV);
 47      LeptonFinder muons(bare_muons_fs, Photon, 0.1, lep_cuts);
 48      declare(muons, "MUONS");
 49
 50      LeptonFinder elecs(bare_elecs_fs, Photon, 0.1, lep_cuts);
 51      declare(elecs, "ELECS");
 52
 53      VetoedFinalState vfs;
 54      vfs.addVetoOnThisFinalState(muons);
 55      vfs.addVetoOnThisFinalState(elecs);
 56      vfs.addVetoOnThisFinalState(neutrino_fs);
 57
 58      FastJets fjets(vfs, JetAlg::ANTIKT, 0.4);
 59      fjets.useInvisibles();
 60      declare(fjets, "jets");
 61
 62      book(h_pull_all     ,4,1,1);
 63      book(h_pull_charged ,5,1,1);
 64    }
 65
 66
 67    /// Perform the per-event analysis
 68    void analyze(const Event& event) {
 69
 70      /**************
 71       *    JETS    *
 72       **************/
 73      const Jets& allJets = apply<FastJets>(event, "jets").jetsByPt(Cuts::pT > 25.0*GeV && Cuts::absrap < 2.5);
 74      const DressedLeptons& all_elecs = apply<LeptonFinder>(event, "ELECS").dressedLeptons();
 75      const DressedLeptons& all_muons = apply<LeptonFinder>(event, "MUONS").dressedLeptons();
 76      Jets goodJets;
 77      for (const Jet & j : allJets) {
 78        bool keep = true;
 79        for (const DressedLepton & el : all_elecs)  keep &= deltaR(j, el) >= 0.2;
 80        if (keep)  goodJets += j;
 81      }
 82      if ( goodJets.size() < 4 )  vetoEvent;
 83
 84      /****************
 85       *    LEPTONS   *
 86       ****************/
 87      DressedLeptons muons, vetoMuons;
 88      for (const DressedLepton & mu : all_muons) {
 89        bool keep = true;
 90        for (const Jet & j : goodJets)  keep &= deltaR(j, mu) >= 0.4;
 91        if (keep && mu.pt() > 15*GeV) {
 92          vetoMuons.push_back(mu);
 93          if (mu.pt() > 25*GeV)  muons.push_back(mu);
 94        }
 95      }
 96
 97      DressedLeptons elecs, vetoElecs;
 98      for (const DressedLepton & el : all_elecs) {
 99        bool keep = true;
100        for (const Jet & j : goodJets)  keep &= deltaR(j, el) >= 0.4;
101        if (keep && el.pt() > 15*GeV) {
102          vetoElecs.push_back(el);
103          if (el.pt() > 25*GeV)  elecs.push_back(el);
104        }
105      }
106
107      if (muons.empty() && elecs.empty())  vetoEvent;
108
109      bool muCandidate = !( muons.size() < 1 || vetoMuons.size() > 1 || vetoElecs.size() > 0 );
110      bool elCandidate = !( elecs.size() < 1 || vetoElecs.size() > 1 || vetoMuons.size() > 0 );
111
112      if (!elCandidate && !muCandidate)  vetoEvent;
113
114      /******************************
115       *    ELECTRON-MUON OVERLAP   *
116       ******************************/
117      for (const DressedLepton & electron : elecs) {
118        for (const DressedLepton & muon : muons) {
119          double d_theta = fabs(muon.theta() - electron.theta());
120          double d_phi = deltaPhi(muon.phi(), electron.phi());
121          if (d_theta < 0.005 && d_phi < 0.005)  vetoEvent;
122        }
123      }
124
125      /****************
126       *  NEUTRINOS   *
127       ****************/
128      const Particles& neutrinos = apply<IdentifiedFinalState>(event, "NEUTRINO_FS").particlesByPt();
129      FourMomentum metVector = FourMomentum(0.,0.,0.,0.);
130      for (const Particle& n : neutrinos) {
131        metVector += n.momentum();
132      }
133      double met = metVector.pt();
134      if (met <= 20*GeV)  vetoEvent;
135
136      if ( (_mT(muCandidate? muons[0] : elecs[0], metVector) + met) <= 60. )  vetoEvent;
137
138      /****************
139       *    B-JETS    *
140       ****************/
141      Jets bJets, wJets;
142      for(Jet & j : goodJets) {
143        bool b_tagged = false;
144        Particles bTags = j.bTags();
145        for ( Particle & b : bTags ) {
146          b_tagged |= b.pT() > 5*GeV;
147        }
148        if (b_tagged)  bJets += j;
149        if (!b_tagged && j.abseta() < 2.1)  wJets += j;
150      }
151
152      if ( bJets.size() < 2 || wJets.size() < 2 )  vetoEvent;
153
154      double pull_angle = fabs(CalculatePullAngle(wJets[0], wJets[1], 0));
155      h_pull_all->fill(pull_angle / Rivet::PI);
156
157      double pull_angle_charged = fabs(CalculatePullAngle(wJets[0], wJets[1], 1));
158      h_pull_charged->fill(pull_angle_charged / Rivet::PI);
159
160    }
161
162    Vector3 CalculatePull(Jet& jet, bool &isCharged) {
163      Vector3 pull(0.0, 0.0, 0.0);
164      double PT = jet.pT();
165      Particles& constituents = jet.particles();
166      Particles charged_constituents;
167      if (isCharged) {
168        for (Particle & p : constituents) {
169          if (p.charge3() != 0)  charged_constituents += p;
170        }
171        constituents = charged_constituents;
172      }
173      // calculate axis
174      FourMomentum axis;
175      for (Particle& p : constituents)  axis += p.momentum();
176      Vector3 J(axis.rap(), axis.phi(MINUSPI_PLUSPI), 0.0);
177      // calculate pull
178      for (Particle & p : constituents) {
179        Vector3 ri = Vector3(p.rap(), p.phi(MINUSPI_PLUSPI), 0.0) - J;
180        while (ri.y() >  Rivet::PI) ri.setY(ri.y() - Rivet::TWOPI);
181        while (ri.y() < -Rivet::PI) ri.setY(ri.y() + Rivet::TWOPI);
182        pull.setX(pull.x() + (ri.mod() * ri.x() * p.pT()) / PT);
183        pull.setY(pull.y() + (ri.mod() * ri.y() * p.pT()) / PT);
184      }
185      return pull;
186    }
187
188    double CalculatePullAngle(Jet& jet1, Jet& axisjet, bool isCharged) {
189      Vector3 pull_vector = CalculatePull(jet1, isCharged);
190      pull_vector = Vector3(1000.*pull_vector.x(), 1000.*pull_vector.y(), 0.);
191      double drap = axisjet.rap() - jet1.rap();
192      double dphi = axisjet.phi(MINUSPI_PLUSPI) - jet1.phi(MINUSPI_PLUSPI);
193      Vector3 j2_vector(drap, dphi, 0.0);
194      return mapAngleMPiToPi(deltaPhi(pull_vector, j2_vector));
195    }
196
197    double _mT(const FourMomentum &l, FourMomentum &nu) const {
198      return sqrt( 2 * l.pT() * nu.pT() * (1 - cos(deltaPhi(l, nu))) );
199    }
200
201    /// Normalise histograms etc., after the run
202    void finalize() {
203      normalize(h_pull_all);
204      normalize(h_pull_charged);
205    }
206
207    /// @}
208
209
210  private:
211
212    Histo1DPtr h_pull_all;
213    Histo1DPtr h_pull_charged;
214
215  };
216
217
218  RIVET_DECLARE_PLUGIN(ATLAS_2015_I1376945);
219
220}