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_I1190891 : public Analysis {
 16  public:
 17
 18    /// Constructor
 19    ATLAS_2012_I1190891()
 20      : Analysis("ATLAS_2012_I1190891")
 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, FastJets::ANTIKT, 0.4), "AntiKtJets04");
 48
 49      // all tracks (to do deltaR with leptons)
 50      declare(ChargedFinalState((Cuts::etaIn(-3.0,3.0))),"cfs");
 51
 52      // Book histograms
 53      book(_hist_etmiss ,"hist_etmiss",10,0.,500.);
 54      book(_hist_meff   ,"hist_m_eff",7,0.,1050.);
 55      book(_count_SR1 ,"count_SR1", 1, 0., 1.);
 56      book(_count_SR2 ,"count_SR2", 1, 0., 1.);
 57    }
 58
 59
 60    /// Perform the per-event analysis
 61    void analyze(const Event& event) {
 62      const double weight = 1.0;
 63      // get the jet candidates
 64      Jets cand_jets;
 65      for (const Jet& jet :
 66               apply<FastJets>(event, "AntiKtJets04").jetsByPt(20.0*GeV) ) {
 67        if ( fabs( jet.eta() ) < 2.5 ) {
 68          cand_jets.push_back(jet);
 69        }
 70      }
 71
 72      // candidate muons
 73      Particles cand_mu;
 74      Particles chg_tracks =
 75        apply<ChargedFinalState>(event, "cfs").particles();
 76      for ( const Particle & mu :
 77                apply<IdentifiedFinalState>(event, "muons").particlesByPt() ) {
 78        double pTinCone = -mu.pT();
 79        for ( const Particle & track : chg_tracks ) {
 80          if ( deltaR(mu.momentum(),track.momentum()) <= 0.2 )
 81            pTinCone += track.pT();
 82        }
 83        if ( pTinCone < 1.8*GeV )
 84          cand_mu.push_back(mu);
 85      }
 86
 87      // candidate electrons
 88      Particles cand_e;
 89      for ( const Particle & e :
 90                apply<IdentifiedFinalState>(event, "elecs").particlesByPt() ) {
 91        double pTinCone = -e.perp();
 92        for ( const Particle & track : chg_tracks ) {
 93          if ( deltaR(e.momentum(),track.momentum()) <= 0.2 )
 94            pTinCone += track.pT();
 95        }
 96        if (pTinCone/e.perp()<0.1) {
 97          cand_e.push_back(e);
 98        }
 99      }
100
101      // resolve jet/lepton ambiguity
102      Jets recon_jets;
103      for ( const Jet& jet : cand_jets ) {
104        bool away_from_e = true;
105        for ( const Particle & e : cand_e ) {
106          if ( deltaR(e.momentum(),jet.momentum()) <= 0.2 ) {
107            away_from_e = false;
108            break;
109          }
110        }
111        if ( away_from_e )
112          recon_jets.push_back( jet );
113      }
114
115      // only keep electrons more than R=0.4 from jets
116      Particles cand2_e;
117      for(unsigned int ie=0;ie<cand_e.size();++ie) {
118        const Particle & e = cand_e[ie];
119        // at least 0.4 from any jets
120        bool away = true;
121        for ( const Jet& jet : recon_jets ) {
122          if ( deltaR(e.momentum(),jet.momentum()) < 0.4 ) {
123            away = false;
124            break;
125          }
126        }
127        // and 0.1 from any muons
128        if ( away ) {
129          for ( const Particle & mu : cand_mu ) {
130            if ( deltaR(mu.momentum(),e.momentum()) < 0.1 ) {
131              away = false;
132              break;
133            }
134          }
135        }
136        // and 0.1 from electrons ( keep higher energy)
137        for(unsigned int ie2=0;ie2<cand_e.size();++ie2) {
138          if(ie==ie2) continue;
139          if ( deltaR(e.momentum(),cand_e[ie2].momentum()) < 0.1 &&
140               e.E() < cand_e[ie2].E() ) {
141            away = false;
142            break;
143          }
144        }
145        // if isolated keep it
146        if ( away )
147          cand2_e.push_back( e );
148      }
149
150      // remove e+e- pairs with mass < 20.
151      Particles recon_e;
152      for(unsigned int ie=0;ie<cand2_e.size();++ie) {
153        bool pass = true;
154        for(unsigned int ie2=0;ie2<cand2_e.size();++ie2) {
155          if(cand2_e[ie].pid()*cand2_e[ie2].pid()>0) continue;
156          double mtest = (cand2_e[ie].momentum()+cand2_e[ie2].momentum()).mass();
157          if(mtest<=20.) {
158            pass = false;
159            break;
160          }
161        }
162        if(pass) recon_e.push_back(cand2_e[ie]);
163      }
164
165      // only keep muons more than R=0.4 from jets
166      Particles cand2_mu;
167      for(unsigned int imu=0;imu<cand_mu.size();++imu) {
168        const Particle & mu = cand_mu[imu];
169        bool away = true;
170        // at least 0.4 from any jets
171        for ( const Jet& jet : recon_jets ) {
172          if ( deltaR(mu.momentum(),jet.momentum()) < 0.4 ) {
173            away = false;
174            break;
175          }
176        }
177        // and 0.1 from any electrona
178        if ( away ) {
179          for ( const Particle & e : cand_e ) {
180            if ( deltaR(mu.momentum(),e.momentum()) < 0.1 ) {
181              away = false;
182              break;
183            }
184          }
185        }
186        if ( away )
187          cand2_mu.push_back( mu );
188      }
189
190      // remove mu+mu- pairs with mass < 20.
191      Particles recon_mu;
192      for(unsigned int imu=0;imu<cand2_mu.size();++imu) {
193        bool pass = true;
194        for(unsigned int imu2=0;imu2<cand2_mu.size();++imu2) {
195          if(cand2_mu[imu].pid()*cand2_mu[imu2].pid()>0) continue;
196          double mtest = (cand2_mu[imu].momentum()+cand2_mu[imu2].momentum()).mass();
197          if(mtest<=20.) {
198            pass = false;
199            break;
200          }
201        }
202        if(pass) recon_mu.push_back(cand2_mu[imu]);
203      }
204
205      // pTmiss
206      Particles vfs_particles =
207        apply<VisibleFinalState>(event, "vfs").particles();
208      FourMomentum pTmiss;
209      for ( const Particle & p : vfs_particles ) {
210        pTmiss -= p.momentum();
211      }
212      double eTmiss = pTmiss.pT();
213
214      // now only use recon_jets, recon_mu, recon_e
215
216      // reject events with less than 4 electrons and muons
217      if ( recon_mu.size() + recon_e.size() < 4 ) {
218        MSG_DEBUG("To few charged leptons left after selection");
219        vetoEvent;
220      }
221
222      // check if passes single lepton trigger
223      bool passSingle =
224        ( !recon_e .empty() && recon_e[0] .perp()>25. )||
225        ( !recon_mu.empty() && recon_mu[0].perp()>20.);
226
227      // or two lepton trigger
228      bool passDouble =
229        ( recon_mu.size()>=2 && recon_mu[1].perp()>12.) ||
230        ( recon_e .size()>=2 && recon_e [1].perp()>17.) ||
231        ( !recon_e.empty() && !recon_mu.empty() &&
232          recon_e[0].perp()>15. &&  recon_mu[0].perp()>10.);
233
234      // must pass a trigger
235      if( !passSingle && !passDouble ) {
236        MSG_DEBUG("Hardest lepton fails trigger");
237        vetoEvent;
238      }
239
240      // calculate meff
241      double meff = eTmiss;
242      for ( const Particle & e  : recon_e  )
243        meff += e.perp();
244      for ( const Particle & mu : recon_mu )
245        meff += mu.perp();
246      for ( const Jet & jet : recon_jets ) {
247        double pT = jet.perp();
248        if(pT>40.) meff += pT;
249      }
250
251      // mass of SFOS pairs closest to the Z mass
252      for(unsigned int ix=0;ix<recon_e.size();++ix) {
253        for(unsigned int iy=ix+1;iy<recon_e.size();++iy) {
254          if(recon_e[ix].pid()*recon_e[iy].pid()>0) continue;
255          double mtest = (recon_e[ix].momentum()+recon_e[iy].momentum()).mass();
256          if(mtest>81.2 && mtest<101.2) vetoEvent;
257        }
258      }
259      for(unsigned int ix=0;ix<recon_mu.size();++ix) {
260        for(unsigned int iy=ix+1;iy<recon_mu.size();++iy) {
261          if(recon_mu[ix].pid()*recon_mu[iy].pid()>0) continue;
262          double mtest = (recon_mu[ix].momentum()+recon_mu[iy].momentum()).mass();
263          if(mtest>81.2 && mtest<101.2) vetoEvent;
264        }
265      }
266
267      // make the control plots
268      _hist_etmiss ->fill(eTmiss,weight);
269      _hist_meff   ->fill(meff  ,weight);
270      // finally the counts
271      if(eTmiss>50.) _count_SR1->fill(0.5,weight);
272      if(meff >300.) _count_SR2->fill(0.5,weight);
273    }
274
275    //@}
276
277    void finalize() {
278      double norm = crossSection()/femtobarn*4.7/sumOfWeights();
279      scale(_hist_etmiss,norm* 50.);
280      scale(_hist_meff  ,norm*150.);
281      scale(_count_SR1,norm);
282      scale(_count_SR2,norm);
283    }
284
285  private:
286
287    /// @name Histograms
288    //@{
289    Histo1DPtr _hist_etmiss;
290    Histo1DPtr _hist_meff;
291    Histo1DPtr _count_SR1;
292    Histo1DPtr _count_SR2;
293    //@}
294
295  };
296
297  // The hook for the plugin system
298  RIVET_DECLARE_PLUGIN(ATLAS_2012_I1190891);
299
300}