Rivet analyses referenceBESIII_2019_I1756876Cross section for $e^+e^-\to D^+D^-\pi^+\pi^-$ for $E_{\text{CMS}}=4.36\to4.60$ GeVExperiment: BESIII (BEPC) Inspire ID: 1756876 Status: VALIDATED Authors:
Beam energies: (2.2, 2.2); (2.2, 2.2); (2.2, 2.2); (2.2, 2.2); (2.3, 2.3); (2.3, 2.3); (2.3, 2.3) GeV Run details:
Cross section for $e^+e^-\to D^+D^-\pi^+\pi^-$ for $E_{\text{CMS}}=4.36\to4.60$ GeV. The cross sections for the subprocesses $e^+e^-\to D_1(2420)^\pm D^\mp$ and $e^+e^-\to \psi(3770) \pi^+\pi^-$ are measured. Source code: BESIII_2019_I1756876.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 Cross section for D+D- pi+ pi- (D1 D and psi(3770) pi+pi- subprocesses)
10 class BESIII_2019_I1756876 : public Analysis {
11 public:
12
13 /// Constructor
14 RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2019_I1756876);
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 histograms
26 book(_nD1D, "/TMP/nD1D");
27 book(_nPsi, "/TMP/nPsi");
28
29 }
30
31 void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
32 for(const Particle &child : p.children()) {
33 if(child.children().empty()) {
34 nRes[child.pid()]-=1;
35 --ncount;
36 }
37 else
38 findChildren(child,nRes,ncount);
39 }
40 }
41
42 /// Perform the per-event analysis
43 void analyze(const Event& event) {
44 const FinalState& fs = apply<FinalState>(event, "FS");
45
46 map<long,int> nCount;
47 int ntotal(0);
48 for (const Particle& p : fs.particles()) {
49 nCount[p.pid()] += 1;
50 ++ntotal;
51 }
52 const FinalState& ufs = apply<FinalState>(event, "UFS");
53 bool matched=false;
54 for(const Particle& p1 : ufs.particles(Cuts::abspid==411)) {
55 int sign = p1.pid()/411;
56 map<long,int> nRes = nCount;
57 int ncount = ntotal;
58 findChildren(p1,nRes,ncount);
59 matched=false;
60 for(const Particle& p2 : ufs.particles(Cuts::pid==-sign*411)) {
61 map<long,int> nRes2 = nRes;
62 int ncount2 = ncount;
63 findChildren(p2,nRes2,ncount2);
64 if(ncount2!=1) continue;
65 matched=true;
66 for(auto const & val : nRes2) {
67 if(abs(val.first)==211) {
68 if(val.second!=1) {
69 matched = false;
70 break;
71 }
72 }
73 else if(val.second!=0) {
74 matched = false;
75 break;
76 }
77 }
78 if(matched) break;
79 }
80 if(matched)
81 break;
82 }
83 if(!matched) vetoEvent;
84 // psi(3770) pi+pi-
85 for(const Particle& p1 : ufs.particles(Cuts::abspid==30443)) {
86 if(p1.children().empty()) continue;
87 map<long,int> nRes = nCount;
88 int ncount = ntotal;
89 findChildren(p1,nRes,ncount);
90 // psi(3770) pi+pi-
91 if(ncount!=2) continue;
92 bool matched = true;
93 for(auto const & val : nRes) {
94 if(abs(val.first)==211) {
95 if(val.second !=1) {
96 matched = false;
97 break;
98 }
99 }
100 else if(val.second!=0) {
101 matched = false;
102 break;
103 }
104 }
105 if(matched) {
106 _nPsi->fill();
107 break;
108 }
109 }
110 // D1 D
111 for(const Particle& p1 : ufs.particles(Cuts::abspid==10413)) {
112 if(p1.children().empty()) continue;
113 int sign = p1.pid()/10413;
114 map<long,int> nRes = nCount;
115 int ncount = ntotal;
116 findChildren(p1,nRes,ncount);
117 matched=false;
118 for(const Particle& p2 : ufs.particles(Cuts::pid==-sign*411)) {
119 map<long,int> nRes2 = nRes;
120 int ncount2 = ncount;
121 findChildren(p2,nRes2,ncount2);
122 if(ncount2!=1) continue;
123 matched=true;
124 for(auto const & val : nRes2) {
125 if(abs(val.first)==211) {
126 if(val.second!=1) {
127 matched = false;
128 break;
129 }
130 }
131 else if(val.second!=0) {
132 matched = false;
133 break;
134 }
135 }
136 if(matched) break;
137 }
138 if(matched) {
139 _nD1D->fill();
140 break;
141 }
142 }
143 }
144
145 /// Normalise histograms etc., after the run
146 void finalize() {
147 double fact = crossSection()/ sumOfWeights()/picobarn;
148 for(unsigned int iy=9;iy<11;++iy) {
149 double sigma,error;
150 if(iy==9) {
151 sigma = _nD1D->val()*fact;
152 error = _nD1D->err()*fact;
153 }
154 else if(iy==10) {
155 sigma = _nPsi->val()*fact;
156 error = _nPsi->err()*fact;
157 }
158 Scatter2D temphisto(refData(1, 1, iy));
159 Scatter2DPtr mult;
160 book(mult,1,1,iy);
161 for (size_t b = 0; b < temphisto.numPoints(); b++) {
162 const double x = temphisto.point(b).x();
163 pair<double,double> ex = temphisto.point(b).xErrs();
164 pair<double,double> ex2 = ex;
165 if(ex2.first ==0.) ex2. first=0.0001;
166 if(ex2.second==0.) ex2.second=0.0001;
167 if (inRange(sqrtS()/GeV, x-ex2.first, x+ex2.second)) {
168 mult ->addPoint(x, sigma, ex, make_pair(error,error));
169 }
170 else {
171 mult ->addPoint(x, 0., ex, make_pair(0.,.0));
172 }
173 }
174 }
175
176 }
177
178 //@}
179
180 /// @name Histograms
181 //@{
182 CounterPtr _nD1D, _nPsi;
183 //@}
184
185 };
186
187
188 // The hook for the plugin system
189 RIVET_DECLARE_PLUGIN(BESIII_2019_I1756876);
190
191
192}
|