1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5#include "Rivet/Projections/LeptonFinder.hh"
  6#include "Rivet/Projections/MissingMomentum.hh"
  7#include "Rivet/Projections/PromptFinalState.hh"
  8
  9namespace Rivet {
 10
 11
 12  /// @brief gamma gamma -> pi+pi-/K+ K-
 13  class TPC_1986_I228072 : public Analysis {
 14  public:
 15
 16    /// Constructor
 17    RIVET_DEFAULT_ANALYSIS_CTOR(TPC_1986_I228072);
 18
 19
 20    /// @name Analysis methods
 21    ///@{
 22
 23    /// Book histograms and initialise projections before the run
 24    void init() {
 25      // Final state
 26      declare(FinalState(),"FS");
 27
 28      // Check CMS energy in range
 29      if (sqrtS()<0.5*GeV || sqrtS()>3.5*GeV) throw Error("Invalid CMS energy for TPC_1986_I228072");
 30
 31      // Book histograms
 32      book(_cPi[0],"/TMP/nPi_06");
 33      book(_cPi[1],"/TMP/nPi_03");
 34      book(_cK    ,"/TMP/nK"    );
 35    }
 36
 37
 38    /// Perform the per-event analysis
 39    void analyze(const Event& event) {
 40      Particles part = apply<FinalState>(event,"FS").particles();
 41      if (part.size()!=2) vetoEvent;
 42      if (part[0].pid()!=-part[1].pid()) vetoEvent;
 43      double cTheta(0.);
 44      bool foundPi(false),foundK(false);
 45      for (const Particle& p : part) {
 46        if (p.pid()==PID::PIPLUS) {
 47          foundPi=true;
 48          cTheta = abs(p.momentum().z()/p.momentum().p3().mod());
 49        }
 50        else if (p.pid()==PID::KPLUS) {
 51          foundK=true;
 52          cTheta = abs(p.momentum().z()/p.momentum().p3().mod());
 53        }
 54      }
 55      if (!foundPi && !foundK) vetoEvent;
 56      if (cTheta>0.6) vetoEvent;
 57      if (foundPi) {
 58        _cPi[0]->fill();
 59        if (cTheta<0.3)
 60          _cPi[1]->fill();
 61      }
 62      else if (foundK)
 63        _cK->fill();
 64    }
 65
 66
 67    /// Normalise histograms etc., after the run
 68    void finalize() {
 69      double fact = crossSection()/nanobarn/sumOfWeights();
 70      for (unsigned int ih=1;ih<6;++ih) {
 71        if (ih==2||ih==3) continue;
 72        double sigma(0.),error(0.);
 73        CounterPtr count;
 74        if (ih==1) {
 75          count = _cPi[0];
 76        }
 77        else if (ih==4) {
 78          count = _cPi[1];
 79        }
 80        else if (ih==5) {
 81          count = _cK;
 82        }
 83        else
 84          assert(false);
 85        if (count->numEntries()==0) continue;
 86        sigma = count->val()*fact;
 87        error = count->err()*fact;
 88        Estimate1DPtr cross;
 89        book(cross, ih, 1, 1);
 90        for (auto& b : cross->bins()) {
 91          if (inRange(sqrtS(), b.xMin(), b.xMax())) {
 92            b.set(sigma, error);
 93          }
 94        }
 95      }
 96    }
 97
 98    ///@}
 99
100
101    /// @name Histograms
102    ///@{
103    CounterPtr _cPi[2],_cK;
104    ///@}
105
106
107  };
108
109
110  RIVET_DECLARE_PLUGIN(TPC_1986_I228072);
111
112}