1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/IdentifiedFinalState.hh"
 4#include "Rivet/Projections/ChargedFinalState.hh"
 5#include "Rivet/Projections/MissingMomentum.hh"
 6#include "Rivet/Projections/FinalState.hh"
 7
 8namespace Rivet {
 9
10
11  /// Muon charge asymmetry in W events at 7 TeV in ATLAS
12  class ATLAS_2011_I892704 : public Analysis {
13  public:
14
15    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2011_I892704);
16
17
18    void init() {
19      IdentifiedFinalState Muons(Cuts::abseta < 2.4 && Cuts::pT > 20*GeV);
20      Muons.acceptIdPair(PID::MUON);
21      declare(Muons, "muons");
22
23      ChargedFinalState CFS(Cuts::abseta < 2.8);
24      declare(CFS, "tracks");
25
26      MissingMomentum missmom(FinalState(Cuts::abseta < 5));
27      declare(missmom, "MissingMomentum");
28
29      /// @todo Will need to register TMP histograms for future histogramming
30      book(_tmp_h_plus, "TMP/plus", refData(1,1,1));
31      book(_tmp_h_minus, "TMP/minus", refData(1,1,1));
32      book(_h_asym, 1, 1, 1);
33    }
34
35
36    void analyze(const Event& event) {
37
38      const IdentifiedFinalState& muons = apply<IdentifiedFinalState>(event, "muons");
39      if (muons.size() < 1) vetoEvent;
40      const ChargedFinalState& tracks = apply<ChargedFinalState>(event, "tracks");
41
42      Particles selected_muons;
43      for (Particle muon : muons.particles()) {
44        FourMomentum testmom = muon.momentum();
45        double ptmu(testmom.pT()), ptsum(-ptmu), ratio(0.);
46        for (Particle track : tracks.particles()) {
47          const FourMomentum& trackmom = track.momentum();
48          if (deltaR(testmom, trackmom) < 0.4) {
49            ptsum += trackmom.pT();
50            ratio  = ptsum/ptmu;
51            if (ratio > 0.2) break;
52          }
53        }
54        if (ratio < 0.2) selected_muons.push_back(muon);
55      }
56      if (selected_muons.size() < 1) vetoEvent;
57
58      const FourMomentum muonmom = selected_muons[0].momentum();
59      const MissingMomentum& missmom = apply<MissingMomentum>(event, "MissingMomentum");
60      FourMomentum missvec = -missmom.visibleMomentum();
61      if (fabs(missvec.Et()) < 25*GeV) vetoEvent;
62
63      double MTW = sqrt( 2 * missvec.pT() * muonmom.pT() * (1 - cos( deltaPhi(missvec.phi(), muonmom.phi()) )) );
64      if (MTW < 40*GeV) vetoEvent;
65
66      Histo1DPtr & htmp = (selected_muons[0].pid() > 0) ? _tmp_h_minus : _tmp_h_plus;
67      htmp->fill(muonmom.eta());
68    }
69
70
71    /// Normalise histograms etc., after the run
72    void finalize() {
73      assert(*_tmp_h_plus == *_tmp_h_minus);
74      for (size_t i = 1; i < _tmp_h_plus->numBins()+1; ++i) {
75        const double num   = _tmp_h_plus->bin(i).sumW() - _tmp_h_minus->bin(i).sumW();
76        const double denom = _tmp_h_plus->bin(i).sumW() + _tmp_h_minus->bin(i).sumW();
77        const double relerr = _tmp_h_plus->bin(i).relErrW()  + _tmp_h_minus->bin(i).relErrW();
78        const double asym = (num != 0 && denom != 0) ? num / denom : 0;
79        const double asym_err = (num != 0 && denom != 0) ? asym*relerr : 0;
80        _h_asym->bin(i).set(asym, asym_err);
81      }
82    }
83
84
85  private:
86
87    Estimate1DPtr _h_asym;
88    /// @todo Will need to register TMP histograms for future histogramming
89    Histo1DPtr  _tmp_h_plus, _tmp_h_minus;
90
91  };
92
93
94  RIVET_DECLARE_ALIASED_PLUGIN(ATLAS_2011_I892704, ATLAS_2011_S9002537);
95
96}