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