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