1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/ChargedFinalState.hh"
  5#include "Rivet/Projections/VisibleFinalState.hh"
  6#include "Rivet/Projections/VetoedFinalState.hh"
  7#include "Rivet/Projections/IdentifiedFinalState.hh"
  8#include "Rivet/Projections/FastJets.hh"
  9#include "Rivet/Tools/RivetMT2.hh"
 10#include "Rivet/Tools/Random.hh"
 11
 12namespace Rivet {
 13
 14
 15  /// @author Peter Richardson
 16  class ATLAS_2012_I1112263 : public Analysis {
 17  public:
 18
 19    /// Constructor
 20    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2012_I1112263);
 21
 22
 23    /// @name Analysis methods
 24    /// @{
 25
 26    /// Book histograms and initialise projections before the run
 27    void init() {
 28
 29      // projection to find the electrons
 30      IdentifiedFinalState elecs(Cuts::abseta < 2.47 && Cuts::pT > 10*GeV);
 31      elecs.acceptIdPair(PID::ELECTRON);
 32      declare(elecs, "elecs");
 33
 34      // projection to find the muons
 35      IdentifiedFinalState muons(Cuts::abseta < 2.4 && Cuts::pT > 10*GeV);
 36      muons.acceptIdPair(PID::MUON);
 37      declare(muons, "muons");
 38
 39      // for pTmiss
 40      declare(VisibleFinalState(Cuts::abseta < 4.9),"vfs");
 41
 42      VetoedFinalState vfs;
 43      vfs.addVetoPairId(PID::MUON);
 44
 45      /// Jet finder
 46      declare(FastJets(vfs, JetAlg::ANTIKT, 0.4), "AntiKtJets04");
 47
 48      // all tracks (to do deltaR with leptons)
 49      declare(ChargedFinalState(Cuts::abseta < 3.0),"cfs");
 50
 51      // Book histograms
 52      {Histo1DPtr tmp; _hist_leptonpT_SR1.push_back(book(tmp,"hist_lepton_pT_1_SR1",11,0.,220.));}
 53      {Histo1DPtr tmp; _hist_leptonpT_SR1.push_back(book(tmp,"hist_lepton_pT_2_SR1", 7,0.,140.));}
 54      {Histo1DPtr tmp; _hist_leptonpT_SR1.push_back(book(tmp,"hist_lepton_pT_3_SR1", 8,0.,160.));}
 55      {Histo1DPtr tmp; _hist_leptonpT_SR2.push_back(book(tmp,"hist_lepton_pT_1_SR2",11,0.,220.));}
 56      {Histo1DPtr tmp; _hist_leptonpT_SR2.push_back(book(tmp,"hist_lepton_pT_2_SR2", 7,0.,140.));}
 57      {Histo1DPtr tmp; _hist_leptonpT_SR2.push_back(book(tmp,"hist_lepton_pT_3_SR2", 8,0.,160.));}
 58      book(_hist_etmiss_SR1_A ,"hist_etmiss_SR1_A",15,10.,310.);
 59      book(_hist_etmiss_SR1_B ,"hist_etmiss_SR1_B", 9,10.,190.);
 60      book(_hist_etmiss_SR2_A ,"hist_etmiss_SR2_A",15,10.,310.);
 61      book(_hist_etmiss_SR2_B ,"hist_etmiss_SR2_B", 9,10.,190.);
 62      book(_hist_mSFOS,"hist_mSFOF",9,0.,180.);
 63
 64      book(_count_SR1 ,"count_SR1", 1, 0., 1.);
 65      book(_count_SR2 ,"count_SR2", 1, 0., 1.);
 66    }
 67
 68
 69    /// Perform the per-event analysis
 70    void analyze(const Event& event) {
 71
 72      // Get the jet candidates
 73      Jets cand_jets = apply<FastJets>(event, "AntiKtJets04").jetsByPt(Cuts::pT > 20*GeV && Cuts::abseta < 2.8);
 74
 75      // Candidate muons
 76      Particles cand_mu;
 77      Particles chg_tracks =
 78        apply<ChargedFinalState>(event, "cfs").particles();
 79      for ( const Particle & mu : apply<IdentifiedFinalState>(event, "muons").particlesByPt() ) {
 80        double pTinCone = -mu.pT();
 81        for ( const Particle & track : chg_tracks ) {
 82          if ( deltaR(mu.momentum(),track.momentum()) <= 0.2 )
 83            pTinCone += track.pT();
 84        }
 85        if ( pTinCone < 1.8*GeV )
 86          cand_mu.push_back(mu);
 87      }
 88
 89      // Candidate electrons
 90      Particles cand_e;
 91      for ( const Particle & e : apply<IdentifiedFinalState>(event, "elecs").particlesByPt() ) {
 92        double eta = e.eta();
 93        // Remove electrons with pT<15 in old veto region
 94        // (NOT EXPLICIT IN THIS PAPER BUT IN SIMILAR 4 LEPTON PAPER and THIS DESCRPITION
 95        //  IS MUCH WORSE SO ASSUME THIS IS DONE)
 96        if ( fabs(eta)>1.37 && fabs(eta) < 1.52 && e.perp()< 15.*GeV)
 97          continue;
 98        double pTinCone = -e.perp();
 99        for ( const Particle & track : chg_tracks ) {
100          if ( deltaR(e.momentum(),track.momentum()) <= 0.2 )
101            pTinCone += track.pT();
102        }
103        if (pTinCone/e.perp()<0.1) {
104          cand_e.push_back(e);
105        }
106      }
107
108      // Resolve jet/lepton ambiguity
109      // (NOT EXPLICIT IN THIS PAPER BUT IN SIMILAR 4 LEPTON PAPER and THIS DESCRPITION
110      //  IS MUCH WORSE SO ASSUME THIS IS DONE)
111      Jets recon_jets;
112      for ( const Jet& jet : cand_jets ) {
113        bool away_from_e = true;
114        for ( const Particle & e : cand_e ) {
115          if ( deltaR(e.momentum(),jet.momentum()) <= 0.2 ) {
116            away_from_e = false;
117            break;
118          }
119        }
120        if ( away_from_e )
121          recon_jets.push_back( jet );
122      }
123
124      // Only keep electrons more than R=0.4 from jets
125      Particles recon_e;
126      for ( const Particle & e : cand_e ) {
127        bool away = true;
128        for ( const Jet& jet : recon_jets ) {
129          if ( deltaR(e.momentum(),jet.momentum()) < 0.4 ) {
130            away = false;
131            break;
132          }
133        }
134        // ... and 0.1 from any muons
135        if ( ! away ) {
136          for ( const Particle & mu : cand_e ) {
137            if ( deltaR(mu.momentum(),e.momentum()) < 0.1 ) {
138              away = false;
139              break;
140            }
141          }
142        }
143        if ( away )
144          recon_e.push_back( e );
145      }
146      // Only keep muons more than R=0.4 from jets
147      Particles recon_mu;
148      for ( const Particle & mu : cand_mu ) {
149        bool away = true;
150        for ( const Jet& jet : recon_jets ) {
151          if ( deltaR(mu.momentum(),jet.momentum()) < 0.4 ) {
152            away = false;
153            break;
154          }
155        }
156        // ... and 0.1 from any electrona
157        if ( ! away ) {
158          for ( const Particle & e : cand_e ) {
159            if ( deltaR(mu.momentum(),e.momentum()) < 0.1 ) {
160              away = false;
161              break;
162            }
163          }
164        }
165        if ( away )
166          recon_mu.push_back( mu );
167      }
168
169      // pTmiss
170      Particles vfs_particles =
171        apply<VisibleFinalState>(event, "vfs").particles();
172      FourMomentum pTmiss;
173      for ( const Particle & p : vfs_particles ) {
174        pTmiss -= p.momentum();
175      }
176      double eTmiss = pTmiss.pT();
177
178      // Now only use recon_jets, recon_mu, recon_e
179
180      // Reject events with wrong number of leptons
181      if ( recon_mu.size() + recon_e.size() != 3 ) {
182        MSG_DEBUG("To few charged leptons left after selection");
183        vetoEvent;
184      }
185
186      // ATLAS calo problem
187      if (rand01() <= 0.42) {
188        for ( const Particle & e : recon_e ) {
189          double eta = e.eta();
190          double phi = e.azimuthalAngle(MINUSPI_PLUSPI);
191          if (inRange(eta, -0.1, 1.5) && inRange(phi, -0.9, -0.5)) vetoEvent;
192        }
193        for ( const Jet & jet : recon_jets ) {
194          const double eta = jet.rapidity();
195          const double phi = jet.azimuthalAngle(MINUSPI_PLUSPI);
196          if (jet.perp() > 40*GeV && inRange(eta, -0.1, 1.5) && inRange(phi, -0.9, -0.5)) vetoEvent;
197        }
198      }
199
200      if ( !( !recon_e .empty() && recon_e[0] .perp() > 25*GeV) &&
201           !( !recon_mu.empty() && recon_mu[0].perp() > 20*GeV) ) {
202        MSG_DEBUG("Hardest lepton fails trigger");
203        vetoEvent;
204      }
205
206      // eTmiss cut
207      if (eTmiss < 50*GeV) vetoEvent;
208
209      // Check at least 1 SFOS pair
210      double mSFOS=1e30, mdiff=1e30*GeV;
211      size_t nSFOS=0;
212      for (size_t ix = 0; ix < recon_e.size(); ++ix) {
213        for (size_t iy = ix+1; iy < recon_e.size(); ++iy) {
214          if (recon_e[ix].pid()*recon_e[iy].pid() > 0) continue;
215          ++nSFOS;
216          double mtest = (recon_e[ix].momentum() + recon_e[iy].momentum()).mass();
217          // Veto is mass<20
218          if (mtest < 20*GeV) vetoEvent;
219          if (fabs(mtest - 90*GeV) < mdiff) {
220            mSFOS = mtest;
221            mdiff = fabs(mtest - 90*GeV);
222          }
223        }
224      }
225      for (size_t ix = 0; ix < recon_mu.size(); ++ix) {
226        for (size_t iy = ix+1; iy < recon_mu.size(); ++iy) {
227          if (recon_mu[ix].pid()*recon_mu[iy].pid() > 0) continue;
228          ++nSFOS;
229          double mtest = (recon_mu[ix].momentum() + recon_mu[iy].momentum()).mass();
230          // Veto is mass < 20*GeV
231          if (mtest < 20*GeV) vetoEvent;
232          if (fabs(mtest - 90*GeV) < mdiff) {
233            mSFOS = mtest;
234            mdiff = fabs(mtest - 90*GeV);
235          }
236        }
237      }
238      // Require at least 1 SFOS pair
239      if (nSFOS == 0) vetoEvent;
240      // b-jet veto in SR!
241      if (mdiff > 10*GeV) {
242        for (const Jet & jet : recon_jets ) {
243          if (jet.bTagged() && rand01() <= 0.60) vetoEvent;
244        }
245      }
246
247	  // Histogram filling
248
249      // Region SR1, Z depleted
250      if (mdiff > 10*GeV) {
251        _count_SR1->fill(0.5);
252        _hist_etmiss_SR1_A->fill(eTmiss);
253        _hist_etmiss_SR1_B->fill(eTmiss);
254        _hist_mSFOS->fill(mSFOS);
255      }
256      // Region SR2, Z enriched
257      else {
258        _count_SR2->fill(0.5);
259        _hist_etmiss_SR2_A->fill(eTmiss);
260        _hist_etmiss_SR2_B->fill(eTmiss);
261      }
262      // Make the control plots
263      // lepton pT
264      size_t ie=0, imu=0;
265      for (size_t ix = 0; ix < 3; ++ix) {
266        Histo1DPtr hist = (mdiff > 10*GeV) ? _hist_leptonpT_SR1[ix] :  _hist_leptonpT_SR2[ix];
267        double pTe  = (ie  < recon_e .size()) ? recon_e [ie ].perp() : -1*GeV;
268        double pTmu = (imu < recon_mu.size()) ? recon_mu[imu].perp() : -1*GeV;
269        if (pTe > pTmu) {
270          hist->fill(pTe);
271          ++ie;
272        } else {
273          hist->fill(pTmu);
274          ++imu;
275        }
276      }
277
278    }
279
280
281    void finalize() {
282      const double norm = crossSection()/femtobarn*2.06/sumOfWeights();
283      // These are number of events at 2.06fb^-1 per 20 GeV
284      for (size_t ix = 0; ix < 3; ++ix) {
285        scale(_hist_leptonpT_SR1[ix], norm*20.);
286        scale(_hist_leptonpT_SR2[ix], norm*20.);
287      }
288      scale(_hist_etmiss_SR1_A, norm*20.);
289      scale(_hist_etmiss_SR1_B, norm*20.);
290      scale(_hist_etmiss_SR2_A, norm*20.);
291      scale(_hist_etmiss_SR2_B, norm*20.);
292      scale(_hist_mSFOS, norm*20.);
293      // These are number of events at 2.06fb^-1
294      scale(_count_SR1, norm);
295      scale(_count_SR2, norm);
296    }
297
298    /// @}
299
300
301  private:
302
303  /// @name Histograms
304  /// @{
305  vector<Histo1DPtr> _hist_leptonpT_SR1, _hist_leptonpT_SR2;
306  Histo1DPtr _hist_etmiss_SR1_A, _hist_etmiss_SR1_B, _hist_etmiss_SR2_A, _hist_etmiss_SR2_B;
307  Histo1DPtr _hist_mSFOS, _count_SR1, _count_SR2;
308  /// @}
309
310  };
311
312  RIVET_DECLARE_PLUGIN(ATLAS_2012_I1112263);
313
314}