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
  9namespace Rivet {
 10
 11
 12  class ATLAS_2011_CONF_2011_098 : public Analysis {
 13  public:
 14
 15    /// Constructor
 16    ATLAS_2011_CONF_2011_098()
 17      : Analysis("ATLAS_2011_CONF_2011_098"),
 18        //debug variables
 19        threeJA(0), threeJB(0), threeJC(0), threeJD(0), bj(0), jets(0), zerolept(0), eTmisscut(0)
 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
 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      /// Jet finder
 41      declare(FastJets(FinalState(), JetAlg::ANTIKT, 0.4), "AntiKtJets04");
 42
 43
 44      // all tracks (to do deltaR with leptons)
 45      declare(ChargedFinalState(Cuts::abseta < 3.0),"cfs");
 46
 47      // for pTmiss
 48      declare(VisibleFinalState(Cuts::abseta < 4.9),"vfs");
 49
 50
 51      /// Book histograms
 52      book(_count_threeJA     ,"count_threeJA", 1, 0., 1.);
 53      book(_count_threeJB     ,"count_threeJB", 1, 0., 1.);
 54      book(_count_threeJC     ,"count_threeJC", 1, 0., 1.);
 55      book(_count_threeJD     ,"count_threeJD", 1, 0., 1.);
 56      book(_hist_meff_1bjet   ,"meff_1bjet", 32, 0., 1600.);
 57      book(_hist_eTmiss_1bjet ,"eTmiss_1bjet", 6, 0., 600.);
 58      book(_hist_pTj_1bjet    ,"pTjet_1bjet", 20, 0., 800.);
 59      book(_hist_meff_2bjet   ,"meff_2bjet", 32, 0., 1600.);
 60      book(_hist_eTmiss_2bjet ,"eTmiss_2bjet", 6, 0., 600.);
 61      book(_hist_pTj_2bjet    ,"pTjet_2bjet", 20, 0., 800.);
 62
 63
 64    }
 65
 66
 67    /// Perform the per-event analysis
 68    void analyze(const Event& event) {
 69
 70
 71      // Temp: calorimeter module failure with 10% acceptance loss;
 72      // region unknown ==> randomly choose 10% of events to be vetoed
 73
 74      if ( rand()/static_cast<double>(RAND_MAX) < 0.1 )
 75        vetoEvent;
 76
 77      Jets tmp_cand_jets = apply<FastJets>(event, "AntiKtJets04").jetsByPt(Cuts::pT>20*GeV && Cuts::abseta < 2.8);
 78
 79      Particles cand_e =
 80        apply<IdentifiedFinalState>(event, "elecs").particlesByPt();
 81      Particles cand_mu =
 82        apply<IdentifiedFinalState>(event, "muons").particlesByPt();
 83      Particles chg_tracks =
 84        apply<ChargedFinalState>(event, "cfs").particles();
 85
 86//cerr << "cand_e.size(): " << cand_e.size() << "   cand_mu.size(): " << cand_mu.size() << '\n';
 87
 88
 89      Jets cand_jets;
 90      for ( const Jet& jet : tmp_cand_jets ) {
 91        if ( fabs( jet.eta() ) >= 2.8 )
 92          cand_jets.push_back( jet );
 93        else {
 94          bool away_from_e = true;
 95          for ( const Particle & e : cand_e ) {
 96            if ( deltaR(e.momentum(),jet.momentum()) <= 0.2 ) {
 97              away_from_e = false;
 98              break;
 99            }
100          }
101          if ( away_from_e )
102            cand_jets.push_back( jet );
103        }
104      }
105
106      Particles cand_lept;
107
108      bool isolated_e;
109      for ( const Particle & e : cand_e ) {
110        isolated_e = true;
111        for ( const Jet& jet : cand_jets ) {
112          if ( deltaR(e.momentum(),jet.momentum()) < 0.4 )
113            isolated_e = false;
114        }
115        if ( isolated_e == true )
116          cand_lept.push_back( e );
117      }
118
119
120      bool isolated_mu;
121      for ( const Particle & mu : cand_mu ) {
122	isolated_mu = true;
123	for ( const Jet& jet : cand_jets ) {
124	  if ( deltaR(mu.momentum(),jet.momentum()) < 0.4 )
125	    isolated_mu = false;
126        }
127        if ( isolated_mu == true)
128          cand_lept.push_back( mu );
129      }
130
131
132      // pTmiss
133      Particles vfs_particles
134        = apply<VisibleFinalState>(event, "vfs").particles();
135      FourMomentum pTmiss;
136      for ( const Particle & p : vfs_particles ) {
137        pTmiss -= p.momentum();
138      }
139      double eTmiss = pTmiss.pT();
140
141
142      // bjets
143      Jets bjets,recon_jets;
144      for (const Jet& j : cand_jets) {
145	if(fabs( j.eta() ) <= 2.8) {
146	  recon_jets.push_back(j);
147	  if ( fabs( j.eta() ) <= 2.5 && j.perp()>50. &&
148	       j.bTagged() && rand()/static_cast<double>(RAND_MAX) < 0.5 )
149	    bjets.push_back(j);
150	}
151      }
152
153      if (bjets.empty())  {
154        MSG_DEBUG("No b-jet axes in acceptance");
155        vetoEvent;
156      }
157
158      ++bj;
159
160
161
162      // Jets event selection
163      if ( recon_jets.size() < 3 )
164        vetoEvent;
165      if ( recon_jets[0].pT() <= 130*GeV )
166        vetoEvent;
167      if ( recon_jets[1].pT() <= 50*GeV ||
168	   recon_jets[2].pT() <= 50*GeV )
169        vetoEvent;
170      ++jets;
171
172      // eTmiss cut
173      if ( eTmiss <= 130*GeV )
174        vetoEvent;
175
176      ++eTmisscut;
177
178      // 0-lepton requirement
179      if ( !cand_lept.empty() )
180        vetoEvent;
181      ++zerolept;
182
183      // m_eff cut
184      double m_eff = eTmiss
185        + recon_jets[0].pT()
186        + recon_jets[1].pT()
187        + recon_jets[2].pT();
188
189      if ( eTmiss / m_eff <= 0.25 )
190        vetoEvent;
191
192
193      // min_dPhi
194      double min_dPhi = 999.999;
195      for ( int i = 0; i < 3; ++i ) {
196        double dPhi = deltaPhi( pTmiss.phi(), recon_jets[i].phi() );
197        min_dPhi = min( min_dPhi, dPhi );
198      }
199
200      if ( min_dPhi <= 0.4 )
201        vetoEvent;
202
203
204
205    // ==================== FILL ====================
206
207
208      // 1 bjet
209      if ( bjets.size() >= 1 ) {
210
211        _hist_meff_1bjet->fill(m_eff);
212        _hist_eTmiss_1bjet->fill(eTmiss);
213        _hist_pTj_1bjet->fill(recon_jets[0].pT());
214
215        // 3JA region
216        if ( m_eff > 200*GeV ) {
217	  ++threeJA;
218        _count_threeJA->fill(0.5);
219        }
220
221        // 3JB region
222        if ( m_eff > 700*GeV ) {
223	  ++threeJB;
224        _count_threeJB->fill(0.5);
225        }
226      }
227
228      // 2 bjets
229      if ( bjets.size() >= 2 ) {
230
231        _hist_meff_2bjet->fill(m_eff);
232        _hist_eTmiss_2bjet->fill(eTmiss);
233        _hist_pTj_2bjet->fill(recon_jets[0].pT());
234
235        // 3JC region
236        if ( m_eff > 500*GeV ) {
237	  ++threeJC;
238          _count_threeJC->fill(0.5);
239        }
240
241        // 3JD region
242        if ( m_eff > 700*GeV ) {
243	  ++threeJD;
244          _count_threeJD->fill(0.5);
245        }
246      }
247
248
249
250
251    }
252
253    /// @}
254
255
256    void finalize() {
257      scale( _hist_meff_1bjet, 50. * 830. * crossSection()/picobarn/sumOfWeights() );
258      scale( _hist_eTmiss_1bjet, 100. * 830. * crossSection()/picobarn/sumOfWeights() );
259      scale( _hist_pTj_1bjet, 40. * 830. * crossSection()/picobarn/sumOfWeights() );
260      scale( _hist_meff_2bjet, 50. * 830. * crossSection()/picobarn/sumOfWeights() );
261      scale( _hist_eTmiss_2bjet, 100. * 830. * crossSection()/picobarn/sumOfWeights() );
262      scale( _hist_pTj_2bjet, 40. * 830. * crossSection()/picobarn/sumOfWeights() );
263
264// cerr<< '\n'<<'\n'
265// << "Saw "
266// << bj << " events aft bjets cut, "
267// << jets << " events aft jet cuts, "
268// << eTmisscut << " events aft eTmiss cut, "
269// << zerolept << " events after 0-lept cut. "
270// << '\n'
271// << threeJA << " 3JA events, "
272// << threeJB << " 3JB events, "
273// << threeJC << " 3JC events, "
274// << threeJD << " 3JD events. "
275// << '\n'
276// ;
277
278    }
279
280
281  private:
282
283    /// @name Histograms
284    /// @{
285    Histo1DPtr _count_threeJA;
286    Histo1DPtr _count_threeJB;
287    Histo1DPtr _count_threeJC;
288    Histo1DPtr _count_threeJD;
289    Histo1DPtr _hist_meff_1bjet;
290    Histo1DPtr _hist_eTmiss_1bjet;
291    Histo1DPtr _hist_pTj_1bjet;
292    Histo1DPtr _hist_meff_2bjet;
293    Histo1DPtr _hist_eTmiss_2bjet;
294    Histo1DPtr _hist_pTj_2bjet;
295
296    /// @}
297
298
299// debug variables
300int threeJA;
301int threeJB;
302int threeJC;
303int threeJD;
304int bj;
305int jets;
306int zerolept;
307int eTmisscut;
308
309  };
310
311
312
313  RIVET_DECLARE_PLUGIN(ATLAS_2011_CONF_2011_098);
314
315}