Rivet analyses referenceSND_2020_I1809286Cross section for $e^+e^-\to\pi^+\pi^-\pi^0$ between 1.15 and 2 GeVExperiment: SND (VEPP-2000) Inspire ID: 1809286 Status: VALIDATED Authors:
Beam energies: ANY Run details:
Cross section for $e^+e^-\to\pi^+\pi^-\pi^0$ between 1.15 and 2 GeV measured by SND. The $\omega\pi^0$, $\rho\pi$ and $\rho(1450)\pi$ subprocesses are also extracted. Source code: SND_2020_I1809286.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FinalState.hh"
4#include "Rivet/Projections/UnstableParticles.hh"
5
6namespace Rivet {
7
8
9 /// @brief e+e- > pi+pi-pi0
10 class SND_2020_I1809286 : public Analysis {
11 public:
12
13 /// Constructor
14 RIVET_DEFAULT_ANALYSIS_CTOR(SND_2020_I1809286);
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 declare(UnstableParticles(), "UFS");
25 book(_c_total, "/TMP/total");
26 book(_c_omega, "/TMP/omega");
27 book(_c_rho , "/TMP/rho" );
28 book(_c_rhop , "/TMP/rhop" );
29 if(inRange(sqrtS()/GeV,1.42,1.48)) {
30 book(_h_x,4,1,1);
31 book(_h_m,4,1,3);
32 }
33 else if(inRange(sqrtS()/GeV,1.65,1.68)) {
34 book(_h_x,4,1,2);
35 book(_h_m,4,1,4);
36 }
37 }
38
39 void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
40 for( const Particle &child : p.children()) {
41 if(child.children().empty()) {
42 --nRes[child.pid()];
43 --ncount;
44 }
45 else
46 findChildren(child,nRes,ncount);
47 }
48 }
49
50 /// Perform the per-event analysis
51 void analyze(const Event& event) {
52
53 const FinalState& fs = apply<FinalState>(event, "FS");
54
55 map<long,int> nCount;
56 int ntotal(0);
57 Particle pip,pim;
58 for (const Particle& p : fs.particles()) {
59 nCount[p.pid()] += 1;
60 ++ntotal;
61 if(p.pid() == 211) pip=p;
62 else if(p.pid()==-211) pim=p;
63 }
64 if(ntotal!=3) vetoEvent;
65 if(nCount[-211]==1&&nCount[211]==1&&nCount[111]==1)
66 _c_total->fill();
67 else
68 vetoEvent;
69 if(_h_x) {
70 _h_x->fill(pip.momentum().p()/sqrtS());
71 _h_x->fill(pim.momentum().p()/sqrtS());
72 _h_m->fill((pip.momentum()+pim.momentum()).mass()/MeV);
73 }
74 const FinalState& ufs = apply<FinalState>(event, "UFS");
75 for (const Particle& p : ufs.particles(Cuts::abspid==223 or
76 Cuts::abspid==113 or Cuts::abspid==213 or
77 Cuts::abspid==100113 or Cuts::abspid==100213)) {
78 if(p.children().empty()) continue;
79 map<long,int> nRes = nCount;
80 int ncount = ntotal;
81 findChildren(p,nRes,ncount);
82 if(ncount!=1) continue;
83 int idOther = 211;
84 if(p.pid()==223 || p.pid()==113 || p.pid()==100113)
85 idOther = 111;
86 else if(p.pid()==213 || p.pid()==100213)
87 idOther = -211;
88 bool matched=true;
89 for(auto const & val : nRes) {
90 if(val.first==idOther ) {
91 if(val.second !=1) {
92 matched = false;
93 break;
94 }
95 }
96 else if(val.second!=0) {
97 matched = false;
98 break;
99 }
100 }
101 if(!matched) continue;
102 if(matched) {
103 if(p.pid()==223)
104 _c_omega->fill();
105 else if(p.pid()==213 || p.pid()==113)
106 _c_rho->fill();
107 else
108 _c_rhop->fill();
109 break;
110 }
111 }
112 }
113
114
115 /// Normalise histograms etc., after the run
116 void finalize() {
117 if(_h_x) {
118 normalize(_h_x,1.,false);
119 normalize(_h_m,1.,false);
120 }
121 double fact = crossSection()/nanobarn/sumOfWeights();
122 for(unsigned int ix=1;ix<4;++ix) {
123 unsigned int ymax = ix!=3 ? 2 : 4;
124 for(unsigned int iy=1;iy<ymax;++iy) {
125 double sigma(0.),error(0.);
126 if(ix<3) {
127 sigma = _c_total->val()*fact;
128 error = _c_total->err()*fact;
129 }
130 else if(ix==3) {
131 if(iy==1) {
132 sigma = _c_rho->val()*fact;
133 error = _c_rho->err()*fact;
134 }
135 else if(iy==2) {
136 sigma = _c_rhop->val()*fact;
137 error = _c_rhop->err()*fact;
138 }
139 else {
140 sigma = _c_omega->val()*fact;
141 error = _c_omega->err()*fact;
142 }
143 }
144 Scatter2D temphisto(refData(ix, 1, iy));
145 Scatter2DPtr mult;
146 book(mult, ix, 1, iy);
147 for (size_t b = 0; b < temphisto.numPoints(); b++) {
148 const double x = temphisto.point(b).x();
149 pair<double,double> ex = temphisto.point(b).xErrs();
150 pair<double,double> ex2 = ex;
151 if(ex2.first ==0.) ex2. first=0.0001;
152 if(ex2.second==0.) ex2.second=0.0001;
153 if (inRange(sqrtS()/GeV, x-ex2.first, x+ex2.second)) {
154 mult->addPoint(x, sigma, ex, make_pair(error,error));
155 }
156 else {
157 mult->addPoint(x, 0., ex, make_pair(0.,.0));
158 }
159 }
160 }
161 }
162 }
163
164 ///@}
165
166
167 /// @name Histograms
168 ///@{
169 CounterPtr _c_total,_c_omega,_c_rho,_c_rhop;
170 Histo1DPtr _h_x,_h_m;
171 ///@}
172
173
174 };
175
176
177 RIVET_DECLARE_PLUGIN(SND_2020_I1809286);
178
179}
|