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/FastJets.hh"
  8#include "Rivet/Projections/VetoedFinalState.hh"
  9
 10namespace Rivet {
 11
 12
 13  /// @author Peter Richardson
 14  class ATLAS_2012_CONF_2012_109 : public Analysis {
 15  public:
 16
 17    /// Constructor
 18    ATLAS_2012_CONF_2012_109()
 19      : Analysis("ATLAS_2012_CONF_2012_109")
 20    {    }
 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 > 20*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      // Jet finder
 40      VetoedFinalState vfs;
 41      vfs.addVetoPairId(PID::MUON);
 42      declare(FastJets(vfs, JetAlg::ANTIKT, 0.4), "AntiKtJets04");
 43
 44      // All tracks (to do deltaR with leptons)
 45      declare(ChargedFinalState(Cuts::abseta < 3.0), "cfs");
 46
 47      // Used for pTmiss (N.B. the real 'vfs' extends beyond 4.5 to |eta| = 4.9)
 48      declare(VisibleFinalState(Cuts::abseta < 4.5), "vfs");
 49
 50      // Book histograms
 51      book(_count_A_tight   ,"count_A_tight"   , 1, 0., 1.);
 52      book(_count_A_medium  ,"count_A_medium"  , 1, 0., 1.);
 53      book(_count_A_loose   ,"count_A_loose"   , 1, 0., 1.);
 54      book(_count_B_tight   ,"count_B_tight"   , 1, 0., 1.);
 55      book(_count_B_medium  ,"count_B_medium"  , 1, 0., 1.);
 56      book(_count_C_tight   ,"count_C_tight"   , 1, 0., 1.);
 57      book(_count_C_medium  ,"count_C_medium"  , 1, 0., 1.);
 58      book(_count_C_loose   ,"count_C_loose"   , 1, 0., 1.);
 59      book(_count_D_tight   ,"count_D_tight"   , 1, 0., 1.);
 60      book(_count_E_tight   ,"count_E_tight"   , 1, 0., 1.);
 61      book(_count_E_medium  ,"count_E_medium"  , 1, 0., 1.);
 62      book(_count_E_loose   ,"count_E_loose"   , 1, 0., 1.);
 63
 64      book(_hist_meff_A_medium ,"meff_A_medium" , 40, 0., 4000.);
 65      book(_hist_meff_A_tight  ,"meff_A_tight"  , 40, 0., 4000.);
 66      book(_hist_meff_B_medium ,"meff_B_medium" , 40, 0., 4000.);
 67      book(_hist_meff_B_tight  ,"meff_B_tight"  , 40, 0., 4000.);
 68      book(_hist_meff_C_medium ,"meff_C_medium" , 40, 0., 4000.);
 69      book(_hist_meff_C_tight  ,"meff_C_tight"  , 40, 0., 4000.);
 70      book(_hist_meff_D        ,"meff_D"        , 40, 0., 4000.);
 71      book(_hist_meff_E_loose  ,"meff_E_loose"  , 40, 0., 4000.);
 72      book(_hist_meff_E_medium ,"meff_E_medium" , 40, 0., 4000.);
 73      book(_hist_meff_E_tight  ,"meff_E_tight"  , 40, 0., 4000.);
 74
 75    }
 76
 77
 78    /// Perform the per-event analysis
 79    void analyze(const Event& event) {
 80
 81      Jets cand_jets = apply<FastJets>(event, "AntiKtJets04").jetsByPt(Cuts::pT > 20*GeV && Cuts::abseta < 4.9);
 82
 83      const Particles cand_e  = apply<IdentifiedFinalState>(event, "elecs").particlesByPt();
 84
 85      // Muon isolation not mentioned in hep-exp 1109.6572 but assumed to still be applicable
 86      Particles cand_mu;
 87      const Particles chg_tracks = apply<ChargedFinalState>(event, "cfs").particles();
 88      const Particles muons = apply<IdentifiedFinalState>(event, "muons").particlesByPt();
 89      for (const Particle& mu : muons) {
 90        double pTinCone = -mu.pT();
 91        for (const Particle& track : chg_tracks) {
 92          if ( deltaR(mu.momentum(),track.momentum()) <= 0.2 ) {
 93            pTinCone += track.pT();
 94          }
 95        }
 96        if ( pTinCone < 1.8*GeV ) cand_mu.push_back(mu);
 97      }
 98
 99      // Resolve jet-lepton overlap for jets with |eta| < 2.8
100      Jets recon_jets;
101      for ( const Jet& jet : cand_jets ) {
102        if ( fabs( jet.eta() ) >= 2.8 ) continue;
103        bool away_from_e = true;
104        for ( const Particle & e : cand_e ) {
105          if ( deltaR(e.momentum(),jet.momentum()) <= 0.2 ) {
106            away_from_e = false;
107            break;
108          }
109        }
110        if ( away_from_e ) recon_jets.push_back( jet );
111      }
112
113      Particles recon_e, recon_mu;
114
115      for ( const Particle & e : cand_e ) {
116        bool away = true;
117        for ( const Jet& jet : recon_jets ) {
118          if ( deltaR(e.momentum(),jet.momentum()) < 0.4 ) {
119            away = false;
120            break;
121          }
122        }
123        if ( away ) recon_e.push_back( e );
124      }
125
126      for ( const Particle & mu : cand_mu ) {
127        bool away = true;
128        for ( const Jet& jet : recon_jets ) {
129          if ( deltaR(mu.momentum(),jet.momentum()) < 0.4 ) {
130            away = false;
131            break;
132          }
133        }
134        if ( away ) recon_mu.push_back( mu );
135      }
136
137      // pTmiss
138      // Based on all candidate electrons, muons and jets, plus everything else with |eta| < 4.5
139      // i.e. everything in our projection "vfs" plus the jets with |eta| > 4.5
140      Particles vfs_particles = apply<VisibleFinalState>(event, "vfs").particles();
141      FourMomentum pTmiss;
142      for ( const Particle & p : vfs_particles ) {
143        pTmiss -= p.momentum();
144      }
145      for ( const Jet& jet : cand_jets ) {
146        if ( fabs( jet.eta() ) > 4.5 ) pTmiss -= jet.momentum();
147      }
148      double eTmiss = pTmiss.pT();
149
150      // no electron pT> 20 or muons pT>10
151      if ( !recon_mu.empty() || !recon_e.empty() ) {
152        MSG_DEBUG("Charged leptons left after selection");
153        vetoEvent;
154      }
155
156      if ( eTmiss <= 160 * GeV ) {
157        MSG_DEBUG("Not enough eTmiss: " << eTmiss << " < 130");
158        vetoEvent;
159      }
160
161      // check the hardest two jets
162      if ( recon_jets.size()<2 ||
163           recon_jets[0].pT() <= 130.0 * GeV ||
164           recon_jets[0].pT() <=  60.0 * GeV ) {
165        MSG_DEBUG("No hard leading jet in " << recon_jets.size() << " jets");
166        vetoEvent;
167      }
168
169      // check the charged and EM fractions of the hard jets to avoid photons
170      for (unsigned int ix = 0; ix < 2; ++ix) {
171        // jets over 100 GeV
172        if (recon_jets[ix].pT() < 100*GeV ||
173            recon_jets[ix].eta() > 2.) continue; ///< @todo Should be |eta|?
174        double fch(0.), fem(0.), eTotal(0.);
175        for(const Particle & part : recon_jets[ix].particles()) {
176          long id = part.abspid();
177          if(PID::charge3(id)!=0)
178            fch += part.E();
179          if (id == PID::PHOTON || id == PID::ELECTRON || id == PID::PI0)
180            fem += part.E();
181        }
182        fch /= eTotal;
183        fem /= eTotal;
184        // remove events with hard photon
185        if (fch < 0.02 || (fch < 0.05 && fem > 0.09)) vetoEvent;
186      }
187
188      // ==================== observables ====================
189
190      int Njets = 0;
191      double min_dPhi_All = 999.999; ///< @todo Use std::numeric_limits!
192      double min_dPhi_2   = 999.999; ///< @todo Use std::numeric_limits!
193      double min_dPhi_3   = 999.999; ///< @todo Use std::numeric_limits!
194      double pTmiss_phi = pTmiss.phi();
195
196      for ( const Jet& jet : recon_jets ) {
197        if ( jet.pT() < 40*GeV ) continue;
198        double dPhi = deltaPhi( pTmiss_phi, jet.phi());
199        if ( Njets < 2 ) min_dPhi_2 = min( min_dPhi_2, dPhi );
200        if ( Njets < 3 ) min_dPhi_3 = min( min_dPhi_3, dPhi );
201        min_dPhi_All = min( min_dPhi_All, dPhi );
202        ++Njets;
203      }
204
205      // inclusive meff
206      double m_eff_inc =  eTmiss;
207      for ( const Jet& jet : recon_jets ) {
208        double perp =  jet.pT();
209        if(perp>40.) m_eff_inc += perp;
210      }
211
212      // region A
213      double m_eff_Nj = eTmiss + recon_jets[0].pT() + recon_jets[1].pT();
214      if( min_dPhi_2 > 0.4 && eTmiss/m_eff_Nj > 0.3 ) {
215        _hist_meff_A_tight ->fill(m_eff_inc);
216        if(eTmiss/m_eff_Nj > 0.4)
217          _hist_meff_A_medium->fill(m_eff_inc);
218        if(m_eff_inc>1900.)
219          _count_A_tight ->fill(0.5);
220        if(m_eff_inc>1300. && eTmiss/m_eff_Nj > 0.4)
221          _count_A_medium->fill(0.5);
222        if(m_eff_inc>1300. && eTmiss/m_eff_Nj > 0.4)
223          _count_A_loose ->fill(0.5);
224      }
225
226      // for rest of regions 3 jets pT> 60 needed
227      if(recon_jets.size()<3 || recon_jets[2].perp()<60.)
228        vetoEvent;
229
230      // region B
231      m_eff_Nj +=  recon_jets[2].perp();
232      if( min_dPhi_3 > 0.4 && eTmiss/m_eff_Nj > 0.25 ) {
233        _hist_meff_B_tight->fill(m_eff_inc);
234        if(eTmiss/m_eff_Nj > 0.3)
235          _hist_meff_B_medium->fill(m_eff_inc);
236        if(m_eff_inc>1900.)
237          _count_B_tight ->fill(0.5);
238        if(m_eff_inc>1300. && eTmiss/m_eff_Nj > 0.3)
239          _count_B_medium->fill(0.5);
240      }
241
242      // for rest of regions 4 jets pT> 60 needed
243      if(recon_jets.size()<4 || recon_jets[3].perp()<60.)
244        vetoEvent;
245
246      // region C
247      m_eff_Nj +=  recon_jets[3].perp();
248      if( min_dPhi_3 > 0.4 && min_dPhi_All > 0.2 && eTmiss/m_eff_Nj > 0.25 ) {
249        _hist_meff_C_tight->fill(m_eff_inc);
250        if( eTmiss/m_eff_Nj > 0.3 )
251          _hist_meff_C_medium->fill(m_eff_inc);
252        if(m_eff_inc>1900.)
253          _count_C_tight ->fill(0.5);
254        if(m_eff_inc>1300. && eTmiss/m_eff_Nj > 0.3)
255          _count_C_medium->fill(0.5);
256        if(m_eff_inc>1000. && eTmiss/m_eff_Nj > 0.3)
257          _count_C_loose ->fill(0.5);
258      }
259
260      // for rest of regions 5 jets pT> 40 needed
261      if(recon_jets.size()<5 || recon_jets[4].perp()<40.)
262        vetoEvent;
263
264      // region D
265      m_eff_Nj +=  recon_jets[4].perp();
266      if( min_dPhi_3 > 0.4 && min_dPhi_All > 0.2 && eTmiss/m_eff_Nj > 0.15 ) {
267        _hist_meff_D->fill(m_eff_inc);
268        if(m_eff_inc>1700.) _count_D_tight ->fill(0.5);
269      }
270
271      // for rest of regions 6 jets pT> 40 needed
272      if(recon_jets.size()<6 || recon_jets[5].perp()<40.)
273        vetoEvent;
274
275      // region E
276      m_eff_Nj +=  recon_jets[5].perp();
277      if( min_dPhi_3 > 0.4 && min_dPhi_All > 0.2 && eTmiss/m_eff_Nj > 0.15 ) {
278        _hist_meff_E_tight->fill(m_eff_inc);
279        if( eTmiss/m_eff_Nj > 0.25 )
280          _hist_meff_E_medium->fill(m_eff_inc);
281        if( eTmiss/m_eff_Nj > 0.3 )
282          _hist_meff_E_loose->fill(m_eff_inc);
283        if(m_eff_inc>1400.) _count_E_tight ->fill(0.5);
284        if(m_eff_inc>1300.&& eTmiss/m_eff_Nj > 0.25 )
285          _count_E_medium->fill(0.5);
286        if(m_eff_inc>1000.&& eTmiss/m_eff_Nj > 0.3  )
287          _count_E_loose ->fill(0.5);
288      }
289    }
290
291
292    void finalize() {
293
294      double norm = crossSection()/femtobarn*5.8/sumOfWeights();
295      // these are number of events at 5.8fb^-1 per 100 GeV
296      scale( _hist_meff_A_medium , 100. * norm );
297      scale( _hist_meff_A_tight  , 100. * norm );
298      scale( _hist_meff_B_medium , 100. * norm );
299      scale( _hist_meff_B_tight  , 100. * norm );
300      scale( _hist_meff_C_medium , 100. * norm );
301      scale( _hist_meff_C_tight  , 100. * norm );
302      scale( _hist_meff_D    , 100. * norm );
303      scale( _hist_meff_E_loose  , 100. * norm );
304      scale( _hist_meff_E_medium , 100. * norm );
305      scale( _hist_meff_E_tight  , 100. * norm );
306      // these are number of events at 5.8fb^-1
307      scale(_count_A_tight  ,norm);
308      scale(_count_A_medium ,norm);
309      scale(_count_A_loose  ,norm);
310      scale(_count_B_tight  ,norm);
311      scale(_count_B_medium ,norm);
312      scale(_count_C_tight  ,norm);
313      scale(_count_C_medium ,norm);
314      scale(_count_C_loose  ,norm);
315      scale(_count_D_tight  ,norm);
316      scale(_count_E_tight  ,norm);
317      scale(_count_E_medium ,norm);
318      scale(_count_E_loose  ,norm);
319    }
320
321    /// @}
322
323  private:
324
325    Histo1DPtr _count_A_tight;
326    Histo1DPtr _count_A_medium;
327    Histo1DPtr _count_A_loose;
328    Histo1DPtr _count_B_tight;
329    Histo1DPtr _count_B_medium;
330    Histo1DPtr _count_C_tight;
331    Histo1DPtr _count_C_medium;
332    Histo1DPtr _count_C_loose;
333    Histo1DPtr _count_D_tight;
334    Histo1DPtr _count_E_tight;
335    Histo1DPtr _count_E_medium;
336    Histo1DPtr _count_E_loose;
337
338    Histo1DPtr _hist_meff_A_medium;
339    Histo1DPtr _hist_meff_A_tight;
340    Histo1DPtr _hist_meff_B_medium;
341    Histo1DPtr _hist_meff_B_tight;
342    Histo1DPtr _hist_meff_C_medium;
343    Histo1DPtr _hist_meff_C_tight;
344    Histo1DPtr _hist_meff_D;
345    Histo1DPtr _hist_meff_E_loose;
346    Histo1DPtr _hist_meff_E_medium;
347    Histo1DPtr _hist_meff_E_tight;
348
349  };
350
351
352  RIVET_DECLARE_PLUGIN(ATLAS_2012_CONF_2012_109);
353
354}