Rivet analyses referenceMARKII_1984_I195739$\gamma\gamma\to\pi^+\pi^-/K^+K^-$ for centre-of-mass energies between 1.6 and 2.5 GeVExperiment: MARKII (PEP) Inspire ID: 195739 Status: VALIDATED Authors:
Beam energies: ANY Run details:
Measurement of the differential cross section for $\gamma\gamma\to\pi^+\pi^-/K^+K^-$ for $1.6 \text{GeV} < W < 2.5 \text{GeV}$. Source code: MARKII_1984_I195739.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FinalState.hh"
4
5namespace Rivet {
6
7
8 /// @brief gamma gamma -> pi+pi-/K+ K-
9 class MARKII_1984_I195739 : public Analysis {
10 public:
11
12 /// Constructor
13 RIVET_DEFAULT_ANALYSIS_CTOR(MARKII_1984_I195739);
14
15
16 /// @name Analysis methods
17 ///@{
18
19 /// Book histograms and initialise projections before the run
20 void init() {
21 declare(FinalState(),"FS");
22 // check CMS energy in range
23 if(sqrtS()<1.6*GeV || sqrtS()>2.5*GeV)
24 throw Error("Invalid CMS energy for MARKII_1984_I195739");
25 book(_cPiK, "/TMP/nPiK_");
26 }
27
28
29 /// Perform the per-event analysis
30 void analyze(const Event& event) {
31 Particles part = applyProjection<FinalState>(event,"FS").particles();
32 if(part.size()!=2) vetoEvent;
33 if(part[0].pid()!=-part[1].pid()) vetoEvent;
34 double cTheta(0.);
35 bool foundPi(false),foundK(false);
36 for(const Particle & p : part) {
37 if(p.pid()==PID::PIPLUS) {
38 foundPi=true;
39 cTheta = abs(p.momentum().z()/p.momentum().p3().mod());
40 }
41 else if(p.pid()==PID::KPLUS) {
42 foundK=true;
43 cTheta = abs(p.momentum().z()/p.momentum().p3().mod());
44 }
45 }
46 if(!foundPi && !foundK) vetoEvent;
47 if(cTheta>0.3) vetoEvent;
48 _cPiK->fill();
49 }
50
51
52 /// Normalise histograms etc., after the run
53 void finalize() {
54 double fact = crossSection()/nanobarn/sumOfWeights();
55 double sigma = _cPiK->val()*fact;
56 double error = _cPiK->err()*fact;
57 Scatter2D temphisto(refData(1, 1, 1));
58 Scatter2DPtr cross;
59 book(cross, 1, 1, 1);
60 for (size_t b = 0; b < temphisto.numPoints(); b++) {
61 const double x = temphisto.point(b).x();
62 pair<double,double> ex = temphisto.point(b).xErrs();
63 pair<double,double> ex2 = ex;
64 if(ex2.first ==0.) ex2. first=0.0001;
65 if(ex2.second==0.) ex2.second=0.0001;
66 if (inRange(sqrtS(), x-ex2.first, x+ex2.second)) {
67 cross->addPoint(x, sigma, ex, make_pair(error,error));
68 }
69 else {
70 cross->addPoint(x, 0., ex, make_pair(0.,.0));
71 }
72 }
73 }
74
75 ///@}
76
77
78 /// @name Histograms
79 ///@{
80 CounterPtr _cPiK;
81 ///@}
82
83
84 };
85
86
87 RIVET_DECLARE_PLUGIN(MARKII_1984_I195739);
88
89}
|