1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/FinalState.hh"
 4#include "Rivet/Projections/Beam.hh"
 5
 6namespace Rivet {
 7
 8
 9  /// @brief e+ e- -> e+e- / mu+mu-
10  class BESIII_2022_I2102287 : public Analysis {
11  public:
12
13    /// Constructor
14    RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2022_I2102287);
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(Beam(), "Beams");
24      declare(FinalState(), "FS");
25      // histograms
26      for (unsigned int ix=0; ix<2; ++ix) {
27	      book(_sigma[ix],1+ix,1,1);
28      }
29      // central beam energy
30      _eCent = getOption<string>("ECENT", std::to_string(sqrtS()/GeV));
31    }
32
33
34    /// Perform the per-event analysis
35    void analyze(const Event& event) {
36      static const double cos34 = cos(34./180.*M_PI);
37      // get the axis, direction of incoming positron
38      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
39      Vector3 axis1 = beams.first .momentum().p3().unit();
40      Vector3 axis2 = beams.second.momentum().p3().unit();
41      if (beams.first.pid()<0) swap(axis1,axis2);
42      // loop over FS particles
43      Particles fs = apply<FinalState>(event,"FS").particles();
44      Particles em,ep,mm,mp;
45      for (const Particle & p : fs) {
46        if (p.abspid()==PID::MUON) {
47          if (p.pid()>0) mm.push_back(p);
48          else           mp.push_back(p);
49        }
50        else if (p.abspid()==PID::ELECTRON) {
51          if (abs(axis1.dot(p.p3().unit()))>cos34) continue;
52          if (p.pid()>0) em.push_back(p);
53          else           ep.push_back(p);
54        }
55        else if (p.pid()!=PID::GAMMA) {
56          vetoEvent;
57        }
58      }
59      if (em.size()==1 && ep.size()==1 && mm.size()==0 && mp.size()==0) {
60        _sigma[0]->fill(_eCent);
61      }
62      else if(mm.size()==1 && mp.size()==1 && em.size()==0 && ep.size()==0) {
63        _sigma[1]->fill(_eCent);
64      }
65    }
66
67
68    /// Normalise histograms etc., after the run
69    void finalize() {
70      scale(_sigma, crossSection()/ sumOfWeights() /nanobarn);
71    }
72
73    /// @}
74
75
76    /// @name Histograms
77    /// @{
78    BinnedHistoPtr<string> _sigma[2];
79    string _eCent;
80    /// @}
81
82
83  };
84
85
86  RIVET_DECLARE_PLUGIN(BESIII_2022_I2102287);
87
88}