1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/Beam.hh"
 4#include "Rivet/Projections/FinalState.hh"
 5
 6namespace Rivet {
 7
 8
 9  /// @brief EEC at 29 GeV
10  class MARKII_1988_I250899 : public Analysis {
11  public:
12
13    /// Constructor
14    RIVET_DEFAULT_ANALYSIS_CTOR(MARKII_1988_I250899);
15
16
17    /// @name Analysis methods
18    ///@{
19
20    /// Book histograms and initialise projections before the run
21    void init() {
22      // Initialise and register projections
23      declare(FinalState(), "FS");
24      for(unsigned int ix=0;ix<3;++ix) {
25	book(_h_EEC [ix],2*ix+1,1,1);
26	book(_h_AEEC[ix],2*ix+2,1,1);
27      }
28      book(_weightSum, "TMP/weightSum");
29    }
30
31
32    /// Perform the per-event analysis
33    void analyze(const Event& event) {
34      const FinalState& fs = apply<FinalState>(event, "FS");
35      if ( fs.particles().size() < 2) vetoEvent;
36      _weightSum->fill();
37      // visible energy
38      double Evis = 0.0;
39      for (const Particle& p : fs.particles()) {
40        Evis += p.E();
41      }
42      double Evis2 = sqr(Evis);
43      // (A)EEC
44      // Need iterators since second loop starts at current outer loop iterator, i.e. no "foreach" here!
45      for (Particles::const_iterator p_i = fs.particles().begin(); p_i != fs.particles().end(); ++p_i) {
46        for (Particles::const_iterator p_j = p_i; p_j != fs.particles().end(); ++p_j) {
47          const Vector3 mom3_i = p_i->momentum().p3();
48          const Vector3 mom3_j = p_j->momentum().p3();
49          const double energy_i = p_i->momentum().E();
50          const double energy_j = p_j->momentum().E();
51          const double thetaij = 180.*mom3_i.unit().angle(mom3_j.unit())/M_PI;
52          double eec = (energy_i*energy_j) / Evis2;
53	  if(p_i != p_j) eec *= 2.;
54          for(unsigned int ix=0;ix<3;++ix) _h_EEC[ix]->fill(thetaij, eec);
55          if (thetaij <90.)
56            for(unsigned int ix=0;ix<3;++ix) _h_AEEC[ix]->fill( thetaij, -eec);
57          else
58            for(unsigned int ix=0;ix<3;++ix) _h_AEEC[ix]->fill( 180.-thetaij, eec);
59        }
60      }
61    }
62
63
64    /// Normalise histograms etc., after the run
65    void finalize() {
66      for(unsigned int ix=0;ix<3;++ix) {
67	// factor due angle in degress not radians
68	scale(_h_EEC [ix], 180.0/M_PI/ *_weightSum);
69	scale(_h_AEEC[ix], 180.0/M_PI/ *_weightSum);
70      }
71    }
72
73    ///@}
74
75
76    /// @name Histograms
77    ///@{
78    Histo1DPtr _h_EEC[3],_h_AEEC[3];
79    CounterPtr _weightSum;
80    ///@}
81
82
83  };
84
85
86  RIVET_DECLARE_PLUGIN(MARKII_1988_I250899);
87
88}