← Back to branch summary

~ubuntu-branches/ubuntu/saucy/rivet/saucy-proposed

« back to all changes in this revision

Viewing changes to src/Analyses/ATLAS_2012_I1095236.cc

  • Committer: Package Import Robot
  • Author(s): Lifeng Sun
  • Date: 2013-05-07 09:24:27 UTC
  • mfrom: (1.2.1) (2.1.2 experimental)
  • Revision ID: package-import@ubuntu.com-20130507092427-rpfijl4mn3y87ek7
Tags: 1.8.3-1
New upstream release.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/Analyses/ATLAS_2012_I1095236.cc
 
1
// -*- C++ -*-
 
2
#include "Rivet/Analysis.hh"
 
3
#include "Rivet/Tools/BinnedHistogram.hh"
 
4
#include "Rivet/RivetAIDA.hh"
 
5
#include "Rivet/Tools/Logging.hh"
 
6
#include "Rivet/Projections/FinalState.hh"
 
7
#include "Rivet/Projections/ChargedFinalState.hh"
 
8
#include "Rivet/Projections/VisibleFinalState.hh"
 
9
#include "Rivet/Projections/VetoedFinalState.hh"
 
10
#include "Rivet/Projections/IdentifiedFinalState.hh"
 
11
#include "Rivet/Projections/FastJets.hh"
 
12
#include "Rivet/Tools/ParticleIdUtils.hh"
 
13
 
 
14
namespace Rivet {
 
15
 
 
16
  /// @author Peter Richardson
 
17
  class ATLAS_2012_I1095236 : public Analysis {
 
18
  public:
 
19
 
 
20
    /// @name Constructors etc.
 
21
    //@{
 
22
 
 
23
    /// Constructor
 
24
    ATLAS_2012_I1095236()
 
25
      : Analysis("ATLAS_2012_I1095236")
 
26
    {    }
 
27
 
 
28
    //@}
 
29
 
 
30
 
 
31
  public:
 
32
 
 
33
    /// @name Analysis methods
 
34
    //@{
 
35
 
 
36
    /// Book histograms and initialise projections before the run
 
37
    void init() {
 
38
 
 
39
      // Projection to find the electrons
 
40
      std::vector<std::pair<double, double> > eta_e;
 
41
      eta_e.push_back(make_pair(-2.47,2.47));
 
42
      IdentifiedFinalState elecs(eta_e, 20.0*GeV);
 
43
      elecs.acceptIdPair(ELECTRON);
 
44
      addProjection(elecs, "elecs");
 
45
 
 
46
      // Projection to find the muons
 
47
      std::vector<std::pair<double, double> > eta_m;
 
48
      eta_m.push_back(make_pair(-2.4,2.4));
 
49
      IdentifiedFinalState muons(eta_m, 10.0*GeV);
 
50
      muons.acceptIdPair(MUON);
 
51
      addProjection(muons, "muons");
 
52
 
 
53
      // Jet finder
 
54
      VetoedFinalState vfs;
 
55
      vfs.addVetoPairId(MUON);
 
56
      addProjection(FastJets(vfs, FastJets::ANTIKT, 0.4), "AntiKtJets04");
 
57
 
 
58
      // All tracks (to do deltaR with leptons)
 
59
      addProjection(ChargedFinalState(-3.0,3.0),"cfs");
 
60
 
 
61
      // Used for pTmiss 
 
62
      addProjection(VisibleFinalState(-4.9,4.9),"vfs");
 
63
 
 
64
      // Book histograms
 
65
      _count_SR0_A1 = bookHistogram1D("count_SR0_A1", 1, 0., 1.);
 
66
      _count_SR0_B1 = bookHistogram1D("count_SR0_B1", 1, 0., 1.);
 
67
      _count_SR0_C1 = bookHistogram1D("count_SR0_C1", 1, 0., 1.);
 
68
      _count_SR0_A2 = bookHistogram1D("count_SR0_A2", 1, 0., 1.);
 
69
      _count_SR0_B2 = bookHistogram1D("count_SR0_B2", 1, 0., 1.);
 
70
      _count_SR0_C2 = bookHistogram1D("count_SR0_C2", 1, 0., 1.);
 
71
      _count_SR1_D  = bookHistogram1D("count_SR1_D" , 1, 0., 1.);
 
72
      _count_SR1_E  = bookHistogram1D("count_SR1_E" , 1, 0., 1.);
 
73
 
 
74
      _hist_meff_SR0_A1   = bookHistogram1D("hist_m_eff_SR0_A1", 14, 400., 1800.);
 
75
      _hist_meff_SR0_A2   = bookHistogram1D("hist_m_eff_SR0_A2", 14, 400., 1800.);
 
76
      _hist_meff_SR1_D_e  = bookHistogram1D("hist_meff_SR1_D_e" , 16, 600., 2200.);
 
77
      _hist_meff_SR1_D_mu = bookHistogram1D("hist_meff_SR1_D_mu", 16, 600., 2200.);
 
78
 
 
79
      _hist_met_SR0_A1    = bookHistogram1D("hist_met_SR0_A1", 14, 0., 700.);
 
80
      _hist_met_SR0_A2    = bookHistogram1D("hist_met_SR0_A2", 14, 0., 700.);
 
81
      _hist_met_SR0_D_e   = bookHistogram1D("hist_met_SR1_D_e" , 15, 0., 600.);
 
82
      _hist_met_SR0_D_mu  = bookHistogram1D("hist_met_SR1_D_mu", 15, 0., 600.);
 
83
 
 
84
    }
 
85
 
 
86
 
 
87
    /// Perform the per-event analysis
 
88
    void analyze(const Event& event) {
 
89
      const double weight = event.weight();
 
90
      
 
91
      Jets cand_jets;
 
92
      const Jets jets = applyProjection<FastJets>(event, "AntiKtJets04").jetsByPt(20.0*GeV);
 
93
      foreach (const Jet& jet, jets) {
 
94
        if ( fabs( jet.momentum().eta() ) < 2.8 ) {
 
95
          cand_jets.push_back(jet);
 
96
        }
 
97
      }
 
98
 
 
99
      const ParticleVector cand_e  = applyProjection<IdentifiedFinalState>(event, "elecs").particlesByPt();
 
100
 
 
101
      const ParticleVector cand_mu = applyProjection<IdentifiedFinalState>(event, "muons").particlesByPt();
 
102
      // Resolve jet-lepton overlap for jets with |eta| < 2.8
 
103
      Jets recon_jets;
 
104
      foreach ( const Jet& jet, cand_jets ) {
 
105
        if ( fabs( jet.momentum().eta() ) >= 2.8 ) continue;
 
106
        bool away_from_e = true;
 
107
        foreach ( const Particle & e, cand_e ) {
 
108
          if ( deltaR(e.momentum(),jet.momentum()) <= 0.2 ) {
 
109
            away_from_e = false;
 
110
            break;
 
111
          }
 
112
        }
 
113
        if ( away_from_e ) recon_jets.push_back( jet );
 
114
      }
 
115
      
 
116
      // get the loose leptons used to define the 0 lepton channel
 
117
      ParticleVector loose_e, loose_mu;
 
118
      foreach ( const Particle & e, cand_e ) {
 
119
        bool away = true;
 
120
        foreach ( const Jet& jet, recon_jets ) {
 
121
          if ( deltaR(e.momentum(),jet.momentum()) < 0.4 ) {
 
122
            away = false;
 
123
            break;
 
124
          }
 
125
        }
 
126
        if ( away ) loose_e.push_back( e );
 
127
      }
 
128
      foreach ( const Particle & mu, cand_mu ) {
 
129
        bool away = true;
 
130
        foreach ( const Jet& jet, recon_jets ) {
 
131
          if ( deltaR(mu.momentum(),jet.momentum()) < 0.4 ) {
 
132
            away = false;
 
133
            break;
 
134
          }
 
135
        }
 
136
        if ( away ) loose_mu.push_back( mu );
 
137
      }
 
138
      // tight leptons for the 1-lepton channel
 
139
      ParticleVector tight_mu;
 
140
      ParticleVector chg_tracks =
 
141
        applyProjection<ChargedFinalState>(event, "cfs").particles();
 
142
      foreach ( const Particle & mu, loose_mu) {
 
143
        if(mu.momentum().perp()<20.) continue;
 
144
        double pTinCone = -mu.momentum().pT();
 
145
        foreach ( const Particle & track, chg_tracks ) {
 
146
          if ( deltaR(mu.momentum(),track.momentum()) <= 0.2 )
 
147
            pTinCone += track.momentum().pT();
 
148
        }
 
149
        if ( pTinCone < 1.8*GeV )
 
150
          tight_mu.push_back(mu);
 
151
      }
 
152
      ParticleVector tight_e;
 
153
      foreach ( const Particle & e, loose_e ) {
 
154
        if(e.momentum().perp()<25.) continue;
 
155
        double pTinCone = -e.momentum().perp();
 
156
        foreach ( const Particle & track, chg_tracks ) {
 
157
          if ( deltaR(e.momentum(),track.momentum()) <= 0.2 )
 
158
            pTinCone += track.momentum().pT();
 
159
        }
 
160
        if (pTinCone/e.momentum().perp()<0.1) {
 
161
          tight_e.push_back(e);
 
162
        }
 
163
      }
 
164
 
 
165
      // pTmiss
 
166
      ParticleVector vfs_particles =
 
167
        applyProjection<VisibleFinalState>(event, "vfs").particles();
 
168
      FourMomentum pTmiss;
 
169
      foreach ( const Particle & p, vfs_particles ) {
 
170
        pTmiss -= p.momentum();
 
171
      }
 
172
      double eTmiss = pTmiss.pT();
 
173
 
 
174
      // get the number of b-tagged jets
 
175
      unsigned int ntagged=0;
 
176
      foreach (const Jet & jet, recon_jets ) {
 
177
        if(jet.momentum().perp()>50. && abs(jet.momentum().eta())<2.5 &&
 
178
           jet.containsBottom() && rand()/static_cast<double>(RAND_MAX)<=0.60)
 
179
          ++ntagged;
 
180
      }
 
181
 
 
182
      // ATLAS calo problem
 
183
      if(rand()/static_cast<double>(RAND_MAX)<=0.42) {
 
184
        foreach ( const Jet & jet, recon_jets ) {
 
185
          double eta = jet.momentum().rapidity();
 
186
          double phi = jet.momentum().azimuthalAngle(MINUSPI_PLUSPI);
 
187
          if(jet.momentum().perp()>50 && eta>-0.1&&eta<1.5&&phi>-0.9&&phi<-0.5)
 
188
            vetoEvent;
 
189
        }
 
190
      }
 
191
 
 
192
      // at least 1 b tag
 
193
      if(ntagged==0) vetoEvent;
 
194
 
 
195
      // minumum Et miss
 
196
      if(eTmiss<80.) vetoEvent;
 
197
 
 
198
      // at least 3 jets pT > 50
 
199
      if(recon_jets.size()<3 || recon_jets[2].momentum().perp()<50.)
 
200
        vetoEvent;
 
201
 
 
202
      // m_eff
 
203
      double m_eff =  eTmiss;
 
204
      for(unsigned int ix=0;ix<3;++ix)
 
205
        m_eff += recon_jets[ix].momentum().perp();
 
206
 
 
207
      // delta Phi
 
208
      double min_dPhi = 999.999;
 
209
      double pTmiss_phi = pTmiss.phi();
 
210
      for(unsigned int ix=0;ix<3;++ix) {
 
211
        min_dPhi = min( min_dPhi, deltaPhi( pTmiss_phi, recon_jets[ix].momentum().phi() ) );
 
212
      }
 
213
 
 
214
      // 0-lepton channels
 
215
      if(loose_e.empty() && loose_mu.empty() &&
 
216
         recon_jets[0].momentum().perp()>130.  && eTmiss>130. &&
 
217
         eTmiss/m_eff>0.25 && min_dPhi>0.4) {
 
218
        // jet charge cut
 
219
        bool jetCharge = true;
 
220
        for(unsigned int ix=0;ix<3;++ix) {
 
221
          if(fabs(recon_jets[ix].momentum().eta())>2.) continue;
 
222
          double trackpT=0;
 
223
          foreach(const Particle & p, recon_jets[ix].particles()) {
 
224
            if(PID::threeCharge(p.pdgId())==0) continue;
 
225
            trackpT += p.momentum().perp();
 
226
          }
 
227
          if(trackpT/recon_jets[ix].momentum().perp()<0.05)
 
228
            jetCharge = false;
 
229
        }
 
230
        if(jetCharge) {
 
231
          // SR0-A region
 
232
          if(m_eff>500.) {
 
233
            _count_SR0_A1->fill(0.5,weight);
 
234
            _hist_meff_SR0_A1->fill(m_eff,weight);
 
235
            _hist_met_SR0_A1 ->fill(eTmiss,weight);
 
236
            if(ntagged>=2) {
 
237
              _count_SR0_A2->fill(0.5,weight);
 
238
              _hist_meff_SR0_A2->fill(m_eff,weight);
 
239
              _hist_met_SR0_A2 ->fill(eTmiss,weight);
 
240
            }
 
241
          }
 
242
          // SR0-B
 
243
          if(m_eff>700.) {
 
244
            _count_SR0_B1->fill(0.5,weight);
 
245
            if(ntagged>=2) _count_SR0_B2->fill(0.5,weight);
 
246
          }
 
247
          // SR0-C
 
248
          if(m_eff>900.) {
 
249
            _count_SR0_C1->fill(0.5,weight);
 
250
            if(ntagged>=2) _count_SR0_C2->fill(0.5,weight);
 
251
          }
 
252
        }
 
253
      }
 
254
 
 
255
      // 1-lepton channels
 
256
      if(tight_e.size() + tight_mu.size() == 1 &&
 
257
         recon_jets.size()>=4 && recon_jets[3].momentum().perp()>50.&&
 
258
         recon_jets[0].momentum().perp()>60.) {
 
259
        Particle lepton = tight_e.empty() ? tight_mu[0] : tight_e[0];
 
260
        m_eff += lepton.momentum().perp() + recon_jets[3].momentum().perp();
 
261
        // transverse mass cut
 
262
        double mT = 2.*(lepton.momentum().perp()*eTmiss-
 
263
                        lepton.momentum().x()*pTmiss.x()-
 
264
                        lepton.momentum().y()*pTmiss.y());
 
265
        mT = sqrt(mT);
 
266
        if(mT>100.&&m_eff>700.) {
 
267
          // D region
 
268
          _count_SR1_D->fill(0.5,weight);
 
269
          if(abs(lepton.pdgId())==ELECTRON) {
 
270
            _hist_meff_SR1_D_e->fill(m_eff,weight);
 
271
            _hist_met_SR0_D_e->fill(eTmiss,weight);
 
272
          }
 
273
          else {
 
274
            _hist_meff_SR1_D_mu->fill(m_eff,weight);
 
275
            _hist_met_SR0_D_mu->fill(eTmiss,weight);
 
276
          }
 
277
          // E region
 
278
          if(eTmiss>200.) {
 
279
            _count_SR1_E->fill(0.5,weight);
 
280
          }
 
281
        }
 
282
      }
 
283
    }
 
284
 
 
285
 
 
286
    void finalize() {
 
287
 
 
288
      double norm = crossSection()/femtobarn*2.05/sumOfWeights();
 
289
      // these are number of events at 2.05fb^-1 per 100 GeV
 
290
      scale( _hist_meff_SR0_A1   , 100. * norm );
 
291
      scale( _hist_meff_SR0_A2   , 100. * norm );
 
292
      scale( _hist_meff_SR1_D_e  , 100. * norm );
 
293
      scale( _hist_meff_SR1_D_mu , 100. * norm );
 
294
      // these are number of events at 2.05fb^-1 per 50 GeV
 
295
      scale( _hist_met_SR0_A1, 50. * norm );
 
296
      scale( _hist_met_SR0_A2, 40. * norm );
 
297
      // these are number of events at 2.05fb^-1 per 40 GeV
 
298
      scale( _hist_met_SR0_D_e , 40. * norm );
 
299
      scale( _hist_met_SR0_D_mu, 40. * norm );
 
300
      // these are number of events at 2.05fb^-1
 
301
      scale(_count_SR0_A1,norm);
 
302
      scale(_count_SR0_B1,norm);
 
303
      scale(_count_SR0_C1,norm);
 
304
      scale(_count_SR0_A2,norm);
 
305
      scale(_count_SR0_B2,norm);
 
306
      scale(_count_SR0_C2,norm);
 
307
      scale(_count_SR1_D ,norm);
 
308
      scale(_count_SR1_E ,norm);
 
309
    }
 
310
 
 
311
    //@}
 
312
 
 
313
 
 
314
  private:
 
315
 
 
316
    AIDA::IHistogram1D* _count_SR0_A1;
 
317
    AIDA::IHistogram1D* _count_SR0_B1;
 
318
    AIDA::IHistogram1D* _count_SR0_C1;
 
319
    AIDA::IHistogram1D* _count_SR0_A2;
 
320
    AIDA::IHistogram1D* _count_SR0_B2;
 
321
    AIDA::IHistogram1D* _count_SR0_C2;
 
322
    AIDA::IHistogram1D* _count_SR1_D;
 
323
    AIDA::IHistogram1D* _count_SR1_E;
 
324
 
 
325
    AIDA::IHistogram1D* _hist_meff_SR0_A1;
 
326
    AIDA::IHistogram1D* _hist_meff_SR0_A2;
 
327
    AIDA::IHistogram1D* _hist_meff_SR1_D_e;
 
328
    AIDA::IHistogram1D* _hist_meff_SR1_D_mu;
 
329
    AIDA::IHistogram1D* _hist_met_SR0_A1;
 
330
    AIDA::IHistogram1D* _hist_met_SR0_A2;
 
331
    AIDA::IHistogram1D* _hist_met_SR0_D_e;
 
332
    AIDA::IHistogram1D* _hist_met_SR0_D_mu;
 
333
    
 
334
  };
 
335
 
 
336
 
 
337
  // This global object acts as a hook for the plugin system
 
338
  DECLARE_RIVET_PLUGIN(ATLAS_2012_I1095236);
 
339
 
 
340
}