1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/FinalState.hh"
 4#include "Rivet/Projections/PromptFinalState.hh"
 5
 6namespace Rivet {
 7
 8
 9  /// @brief A study of radiative muon pair events at Z0 energies
10  class DELPHI_1994_I375963 : public Analysis {
11  public:
12
13    /// Constructor
14    RIVET_DEFAULT_ANALYSIS_CTOR(DELPHI_1994_I375963);
15
16
17    /// @name Analysis methods
18    /// @{
19
20    /// Book histograms and initialise projections before the run
21    void init() {
22
23      // FinalState of prompt photons and muons
24      declare(PromptFinalState(Cuts::abspid == PID::MUON && Cuts::energy > 20*GeV), "muons");
25      declare(PromptFinalState(Cuts::abspid == PID::PHOTON && Cuts::energy > 2*GeV), "photons");
26
27      // Book histograms
28      book(_h["ph_energy"], 1, 1, 1);
29      book(_h["ph_angle"], 2, 1, 1);
30    }
31
32
33    /// Perform the per-event analysis
34    void analyze(const Event& event) {
35
36      Particles muons = apply<PromptFinalState>(event, "muons").particles();
37      Particles photons = apply<PromptFinalState>(event, "photons").particles();
38
39      iselect(muons, [](const Particle& muon) {
40        double theta = muon.theta()/M_PI * 180.;
41        return (theta > 20. && theta < 160.);
42      });
43
44      iselect(photons, [](const Particle& photon) {
45        double theta = photon.theta()/M_PI * 180.;
46        double phi = photon.phi()/M_PI * 180.;
47
48        bool femc = (theta > 10. && theta < 36.5) || (theta > 143.5 && theta < 170.);
49        bool hpc = theta > 43. && theta < 137. && abs(fmod(phi, 15.)) > 1.5 && abs(theta-90.) > 2.;
50
51        return (femc || hpc);
52      });
53
54      if (muons.size() != 2) vetoEvent;
55      if (photons.size() == 0) vetoEvent;
56
57      for (const Particle& photon : photons) {
58        double minAngle = 1000.;
59        for (const Particle& muon : muons) {
60          double angle = photon.angle(muon) / M_PI * 180.;
61          if (angle < minAngle) {
62            minAngle = angle;
63          }
64        }
65        if (minAngle > 5.) {
66          _h["ph_energy"]->fill(photon.energy()/GeV);
67          _h["ph_angle"]->fill(minAngle);
68        }
69      }
70
71    }
72
73
74    /// Normalise histograms etc., after the run
75    void finalize() {
76
77      scale({_h["ph_energy"], _h["ph_angle"]}, 1./sumW());
78
79    }
80
81    /// @}
82
83
84    /// @name Histograms
85    /// @{
86    map<string, Histo1DPtr> _h;
87    /// @}
88
89
90  };
91
92
93  RIVET_DECLARE_PLUGIN(DELPHI_1994_I375963);
94
95}