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/IdentifiedFinalState.hh"
  7#include "Rivet/Projections/VetoedFinalState.hh"
  8#include "Rivet/Projections/FastJets.hh"
  9#include "Rivet/Tools/Random.hh"
 10
 11namespace Rivet {
 12
 13
 14  class ATLAS_2012_I1126136 : public Analysis {
 15  public:
 16
 17    /// Constructor
 18    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2012_I1126136);
 19
 20
 21    /// @name Analysis methods
 22    /// @{
 23
 24    /// Book histograms and initialize projections before the run
 25    void init() {
 26
 27      // projection to find the electrons
 28      IdentifiedFinalState elecs(Cuts::abseta < 2.47 && Cuts::pT > 20*GeV);
 29      elecs.acceptIdPair(PID::ELECTRON);
 30      declare(elecs, "elecs");
 31
 32      // projection to find the muons
 33      IdentifiedFinalState muons(Cuts::abseta < 2.4 && Cuts::pT > 10*GeV);
 34      muons.acceptIdPair(PID::MUON);
 35      declare(muons, "muons");
 36
 37      // Jet finder
 38      VetoedFinalState vfs;
 39      vfs.addVetoPairId(PID::MUON);
 40      declare(FastJets(vfs, JetAlg::ANTIKT, 0.4), "AntiKtJets04");
 41
 42      // for pTmiss
 43      declare(VisibleFinalState(Cuts::abseta < 4.9),"vfs");
 44
 45      // Book histograms
 46      book(_count_SR_A     ,"count_SR_A"    , 1, 0., 1.);
 47      book(_count_SR_B     ,"count_SR_B"    , 1, 0., 1.);
 48
 49      book(_hist_mjjj1  ,"hist_mjjj1" , 30 , 0.   , 600.  );
 50      book(_hist_mjjj2  ,"hist_mjjj2" , 30 , 0.   , 600.  );
 51      book(_hist_ETmiss ,"hist_ETmiss", 20 , 100. , 600.  );
 52      book(_hist_mT2    ,"hist_mT2"   , 200, 0.   , 1000. );
 53
 54    }
 55
 56    /// Perform the per-event analysis
 57    void analyze(const Event& event) {
 58
 59      // pTmiss
 60      FourMomentum pTmiss;
 61      for (const Particle& p : apply<VisibleFinalState>(event, "vfs").particles() ) {
 62        pTmiss -= p.momentum();
 63      }
 64      double ETmiss = pTmiss.pT();
 65
 66      // require eTmiss > 150
 67      if (ETmiss < 150*GeV) vetoEvent;
 68
 69      // get the candiate jets
 70      Jets cand_jets = apply<FastJets>(event, "AntiKtJets04").jetsByPt(Cuts::pT > 20*GeV && Cuts::abseta < 4.5);
 71
 72      // find the electrons
 73      Particles cand_e;
 74      for( const Particle& e : apply<IdentifiedFinalState>(event, "elecs").particlesByPt()) {
 75        // remove any leptons within 0.4 of any candidate jets
 76        bool e_near_jet = false;
 77        for ( const Jet& jet : cand_jets ) {
 78          double dR = deltaR(e, jet);
 79          if (inRange(dR, 0.2, 0.4)) {
 80            e_near_jet = true;
 81            break;
 82          }
 83        }
 84        if ( e_near_jet ) continue;
 85        cand_e.push_back(e);
 86      }
 87
 88      // find the muons
 89      Particles cand_mu;
 90      for( const Particle& mu : apply<IdentifiedFinalState>(event, "muons").particlesByPt()) {
 91        // remove any leptons within 0.4 of any candidate jets
 92        bool mu_near_jet = false;
 93        for ( const Jet& jet : cand_jets ) {
 94          if ( deltaR(mu, jet) < 0.4 ) {
 95            mu_near_jet = true;
 96            break;
 97          }
 98        }
 99        if ( mu_near_jet ) continue;
100        cand_mu.push_back(mu);
101      }
102
103      // veto events with leptons
104      if( ! cand_e.empty() || ! cand_mu.empty() )
105        vetoEvent;
106
107      // discard jets that overlap with electrons
108      Jets recon_jets;
109      for ( const Jet& jet : cand_jets ) {
110        if (jet.abseta() > 2.8 || jet.pT() < 30*GeV) continue;
111        bool away_from_e = true;
112        for (const Particle& e : cand_e ) {
113          if ( deltaR(e, jet) < 0.2 ) {
114            away_from_e = false;
115            break;
116          }
117        }
118        if ( away_from_e ) recon_jets.push_back( jet );
119      }
120
121      // find b jets
122      Jets tight_bjets,loose_bjets;
123      for(const Jet& jet : recon_jets) {
124        /// @todo Should be abseta?
125        if (!jet.bTagged() && jet.eta()>2.5) continue;
126        double prob = rand01();
127        if (prob <= 0.60) tight_bjets.push_back(jet);
128        if (prob <= 0.75) loose_bjets.push_back(jet);
129      }
130
131      // require >=1 tight or >=2 loose b-jets
132      if (! ( !tight_bjets.empty() || loose_bjets.size()>=2) )
133        vetoEvent;
134
135      // must be at least 6 jets with pT>30
136      if (recon_jets.size()<6 ) vetoEvent;
137
138      // hardest > 130
139      if (recon_jets[0].perp() < 130. ) vetoEvent;
140
141      // three hardest jets must be separated from etmiss
142      for (unsigned int ix=0;ix<3;++ix) {
143        if (deltaPhi(recon_jets[ix].momentum(),pTmiss)<0.2*PI)
144          vetoEvent;
145      }
146
147      // remove events with tau like jets
148      for (unsigned int ix=3;ix<recon_jets.size();++ix) {
149        // skip jets seperated from eTmiss
150        if (deltaPhi(recon_jets[ix].momentum(),pTmiss)>=0.2*PI)
151          continue;
152        // check the number of tracks between 1 and 4
153        unsigned int ncharged=0;
154        for ( const Particle & particle : recon_jets[ix].particles()) {
155          if (PID::charge3(particle.pid())!=0) ++ncharged;
156        }
157        if (ncharged==0 || ncharged>4) continue;
158        // calculate transverse mass and reject if < 100
159        double mT = 2.*recon_jets[ix].perp()*ETmiss
160          -recon_jets[ix].px()*pTmiss.px()
161          -recon_jets[ix].py()*pTmiss.py();
162        if (mT<100.) vetoEvent;
163      }
164
165      // if 2 loose b-jets apply mT cut
166      if (loose_bjets.size()>=2) {
167        // find b-jet closest to eTmiss
168        double minR(1e30);
169        unsigned int ijet(0);
170        for(unsigned int ix=0;ix<loose_bjets.size();++ix) {
171          double dR = deltaR(loose_bjets[ix].momentum(),pTmiss);
172          if(dR<minR) {
173            minR=dR;
174            ijet = ix;
175          }
176        }
177        double  mT = 2.*loose_bjets[ijet].perp()*ETmiss
178          -loose_bjets[ijet].px()*pTmiss.px()
179          -loose_bjets[ijet].py()*pTmiss.py();
180        if(mT<170.) vetoEvent;
181      }
182
183      // 1 tight b-jet apply mT cut
184      if(tight_bjets.size()==1) {
185        for(unsigned int ix=0;ix<4;++ix) {
186          double mT = 2.*recon_jets[ix].perp()*ETmiss
187            -recon_jets[ix].px()*pTmiss.px()
188            -recon_jets[ix].py()*pTmiss.py();
189          if(mT<175.) vetoEvent;
190        }
191      }
192
193      // find the closest triplet of jets in (eta,phi)
194      unsigned int j1(0),j2(0),j3(0);
195      double minR2(1e30);
196      for(unsigned int i1=0;i1<recon_jets.size();++i1) {
197        for(unsigned int i2=i1+1;i2<recon_jets.size();++i2) {
198          for(unsigned int i3=i2+1;i3<recon_jets.size();++i3) {
199            double delR2 =
200              sqr(deltaR(recon_jets[i1].momentum(),recon_jets[i2].momentum())) +
201              sqr(deltaR(recon_jets[i1].momentum(),recon_jets[i3].momentum())) +
202              sqr(deltaR(recon_jets[i2].momentum(),recon_jets[i3].momentum()));
203            if(delR2<minR2) {
204              minR2=delR2;
205              j1=i1;
206              j2=i2;
207              j3=i3;
208            }
209          }
210        }
211      }
212      // 4-momentum and mass of first triplet
213      FourMomentum pjjj1 = recon_jets[j1].momentum() +
214        recon_jets[j2].momentum()+ recon_jets[j3].momentum();
215      double mjjj1 = pjjj1.mass();
216
217      // find the second triplet
218      unsigned int j4(0),j5(0),j6(0);
219      minR2=0.;
220      for(unsigned int i1=0;i1<recon_jets.size();++i1) {
221        if(i1==j1||i1==j2||i1==j3) continue;
222        for(unsigned int i2=i1+1;i2<recon_jets.size();++i2) {
223          if(i2==j1||i2==j2||i2==j3) continue;
224          for(unsigned int i3=i2+1;i3<recon_jets.size();++i3) {
225            if(i3==j1||i3==j2||i3==j3) continue;
226            double delR2 =
227              sqr(deltaR(recon_jets[i1].momentum(),recon_jets[i2].momentum())) +
228              sqr(deltaR(recon_jets[i1].momentum(),recon_jets[i3].momentum())) +
229              sqr(deltaR(recon_jets[i2].momentum(),recon_jets[i3].momentum()));
230            if(delR2<minR2) {
231              minR2=delR2;
232              j4=i1;
233              j5=i2;
234              j6=i3;
235            }
236          }
237        }
238      }
239
240      // 4-momentum and mass of first triplet
241      FourMomentum pjjj2 = recon_jets[j4].momentum() +
242        recon_jets[j5].momentum()+ recon_jets[j6].momentum();
243      double mjjj2 = pjjj2.mass();
244
245      _hist_mjjj1->fill(mjjj1);
246      _hist_mjjj2->fill(mjjj2);
247      // require triplets in 80<mjjj<270
248      if(mjjj1<80.||mjjj1>270.||mjjj2<80.||mjjj2>270.)
249        vetoEvent;
250
251      // counts in signal regions
252      _count_SR_A->fill(0.5);
253      if(ETmiss>260.) _count_SR_B->fill(0.5);
254
255      _hist_ETmiss->fill(ETmiss);
256      const double m_T2 = mT2(pjjj1, pjjj2, pTmiss, 0.0); // zero mass invisibles
257      _hist_mT2->fill(m_T2);
258    }
259    /// @}
260
261
262    void finalize() {
263
264      double norm = 4.7* crossSection()/sumOfWeights()/femtobarn;
265      scale(_count_SR_A ,     norm );
266      scale(_count_SR_B ,     norm );
267      scale(_hist_mjjj1 , 20.*norm );
268      scale(_hist_ETmiss, 50.*norm );
269      scale(_hist_mjjj2 , 20.*norm );
270      scale(_hist_mT2   ,     norm );
271
272    }
273
274  private:
275
276    /// @name Histograms
277    /// @{
278    Histo1DPtr _count_SR_A;
279    Histo1DPtr _count_SR_B;
280
281    Histo1DPtr _hist_mjjj1;
282    Histo1DPtr _hist_mjjj2;
283    Histo1DPtr _hist_ETmiss;
284    Histo1DPtr _hist_mT2;
285    /// @}
286
287  };
288
289  RIVET_DECLARE_PLUGIN(ATLAS_2012_I1126136);
290
291}