Rivet analyses referenceBESIII_2023_I2677290Cross section and form factors for $e^+e^-\to \Lambda_c^+\bar{\Lambda}_c^-$ below 4.95 GeV.Experiment: (BEPC) Inspire ID: 2677290 Status: VALIDATED NOHEPDATA Authors:
Beam energies: ANY Run details:
Cross section and form factors for $e^+e^-\to \Lambda_c^+\bar{\Lambda}_c^-$ below 4.95 GeV. Beam energy must be specified as analysis option "ENERGY" when rivet-merging samples. Source code: BESIII_2023_I2677290.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/Beam.hh"
4#include "Rivet/Projections/FinalState.hh"
5#include "Rivet/Projections/UnstableParticles.hh"
6
7namespace Rivet {
8
9
10 /// @brief e+ e- > Lambda_c+ Lambda_c-
11 class BESIII_2023_I2677290 : public Analysis {
12 public:
13
14 /// Constructor
15 RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2023_I2677290);
16
17
18 /// @name Analysis methods
19 /// @{
20
21 /// Book histograms and initialise projections before the run
22 void init() {
23 // Initialise and register projections
24 declare(Beam(), "Beams");
25 declare(UnstableParticles(Cuts::abspid==4122), "UFS");
26 declare(FinalState(), "FS");
27 // histograms
28 book(_h_cThetaL,"cThetaL",20,-1.,1.);
29 book(_wsum,"TMP/wsum");
30
31 vector<string> energies({"4.6119", "4.628", "4.6409", "4.6612", "4.6819", "4.6988",
32 "4.7397", "4.75", "4.7805", "4.8431", "4.918", "4.9509"});
33 for(const string& en : energies) {
34 double end = std::stod(en)*GeV;
35 if(isCompatibleWithSqrtS(end)) {
36 _ecms = en;
37 break;
38 }
39 }
40 if(_ecms.empty()) MSG_ERROR("Beam energy incompatible with analysis.");
41 }
42
43 void findChildren(const Particle& p, map<long,int>& nRes, int & ncount) {
44 for (const Particle& child : p.children()) {
45 if (child.children().empty()) {
46 nRes[child.pid()]-=1;
47 --ncount;
48 }
49 else {
50 findChildren(child,nRes,ncount);
51 }
52 }
53 }
54
55 /// Perform the per-event analysis
56 void analyze(const Event& event) {
57 // get the axis, direction of incoming electron
58 const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
59 Vector3 axis;
60 if (beams.first.pid()>0) {
61 axis = beams.first .mom().p3().unit();
62 }
63 else {
64 axis = beams.second.mom().p3().unit();
65 }
66 // types of final state particles
67 const FinalState& fs = apply<FinalState>(event, "FS");
68 map<long,int> nCount;
69 int ntotal(0);
70 for (const Particle& p : fs.particles()) {
71 nCount[p.pid()] += 1;
72 ++ntotal;
73 }
74 const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
75 Particle Lambda, LamBar;
76 bool matched(false);
77 for (const Particle& p : ufs.particles()) {
78 if (p.children().empty()) continue;
79 map<long,int> nRes=nCount;
80 int ncount = ntotal;
81 findChildren(p,nRes,ncount);
82 matched=false;
83 // check for antiparticle
84 for (const Particle& p2 : ufs.particles(Cuts::pid==-p.pid())) {
85 if (p2.children().empty()) continue;
86 map<long,int> nRes2=nRes;
87 int ncount2 = ncount;
88 findChildren(p2,nRes2,ncount2);
89 if (ncount2==0) {
90 matched = true;
91 for (const auto& val : nRes2) {
92 if (val.second!=0) {
93 matched = false;
94 break;
95 }
96 }
97 // found baryon and antibaryon
98 if (matched) {
99 if (p.pid()>0) {
100 Lambda = p;
101 LamBar = p2;
102 }
103 else {
104 Lambda = p2;
105 LamBar = p;
106 }
107 break;
108 }
109 }
110 }
111 if (matched) break;
112 }
113 if (!matched) vetoEvent;
114 const double cosL = axis.dot(Lambda.mom().p3().unit());
115 _wsum->fill();
116 _h_cThetaL->fill(cosL);
117 }
118
119 pair<double,pair<double,double> > calcAlpha0(Histo1DPtr hist) {
120 if (hist->numEntries()==0.) return make_pair(0.,make_pair(0.,0.));
121 const double d = 3./(pow(hist->xMax(),3)-pow(hist->xMin(),3));
122 const double c = 3.*(hist->xMax()-hist->xMin())/(pow(hist->xMax(),3)-pow(hist->xMin(),3));
123 double sum1(0.),sum2(0.),sum3(0.),sum4(0.),sum5(0.);
124 for (const auto& bin : hist->bins()) {
125 const double Oi = bin.sumW();
126 if (Oi==0.) continue;
127 const double a = d*(bin.xMax() - bin.xMin());
128 const double b = d/3.*(pow(bin.xMax(),3) - pow(bin.xMin(),3));
129 const double Ei = bin.errW();
130 sum1 += a*Oi/sqr(Ei);
131 sum2 += b*Oi/sqr(Ei);
132 sum3 += sqr(a)/sqr(Ei);
133 sum4 += sqr(b)/sqr(Ei);
134 sum5 += a*b/sqr(Ei);
135 }
136 // calculate alpha
137 const double alpha = (-c*sum1 + sqr(c)*sum2 + sum3 - c*sum5)/(sum1 - c*sum2 + c*sum4 - sum5);
138 // and error
139 const double cc = -pow((sum3 + sqr(c)*sum4 - 2*c*sum5),3);
140 const double bb = -2*sqr(sum3 + sqr(c)*sum4 - 2*c*sum5)*(sum1 - c*sum2 + c*sum4 - sum5);
141 const double aa = sqr(sum1 - c*sum2 + c*sum4 - sum5)*(-sum3 - sqr(c)*sum4 + sqr(sum1 - c*sum2 + c*sum4 - sum5) + 2*c*sum5);
142 double dis = sqr(bb)-4.*aa*cc;
143 if (dis>0.) {
144 dis = sqrt(dis);
145 return make_pair(alpha,make_pair(0.5*(-bb+dis)/aa,-0.5*(-bb-dis)/aa));
146 }
147 else {
148 return make_pair(alpha,make_pair(0.,0.));
149 }
150 }
151
152 /// Normalise histograms etc., after the run
153 void finalize() {
154 // storage of the values to fill histos
155 const double mn = 2.28646;
156 const double tau = 4.*sqr(mn/sqrtS());
157 const double beta = sqrt(1.-tau);
158 const double alpha = 7.2973525693e-3;
159 const double GeV2pb = 0.3893793721e9;
160 scale(_wsum, crossSection()/ sumOfWeights()/picobarn);
161 // prefactor and sigma
162 double sigma0 = 4.*M_PI*sqr(alpha/sqrtS())*beta*GeV2pb;
163 // calculate alpha0
164 pair<double,pair<double,double> > alpha0 = calcAlpha0(_h_cThetaL);
165 // Geff
166 pair<double,double> Geff = make_pair(1e2*sqrt(3.*_wsum->val()/(sigma0*(1 + 0.5*tau))),
167 1e2*1.5*_wsum->val()/(sigma0*(1 + 0.5*tau)));
168 // GM
169 double GMv = 1e2*sqrt(6.*((1+alpha0.first)*_wsum->val())/((3+alpha0.first)*sigma0));
170 double GMe1 = 1e2*sqrt((3*(sqr(_wsum->err())*sqr(1 + alpha0.first)*sqr(3 + alpha0.first) + 4*sqr(alpha0.second.first)*sqr(_wsum->val())))/
171 (2.*(1 + alpha0.first)*pow(3 + alpha0.first,3)*_wsum->val()*sigma0));
172 double GMe2 = 1e2*sqrt((3*(sqr(_wsum->err())*sqr(1 + alpha0.first)*sqr(3 + alpha0.first) + 4*sqr(alpha0.second.second)*sqr(_wsum->val())))/
173 (2.*(1 + alpha0.first)*pow(3 + alpha0.first,3)*_wsum->val()*sigma0));
174 pair<double,pair<double,double> > GM = make_pair(std::move(GMv), make_pair(std::move(GMe1), std::move(GMe2)));
175 // ratio
176 pair<double,pair<double,double>> R;
177 R.first = sqrt((1 - alpha0.first)/(tau + alpha0.first*tau));
178 R.second.first = R.first*alpha0.second.first /(1.-sqr(alpha0.first));
179 R.second.second = R.first*alpha0.second.second/(1.-sqr(alpha0.first));
180 for(unsigned int ix=1; ix<6; ++ix) {
181 double val;
182 pair<double,double> err;
183 if (ix==1) {
184 val = _wsum->val();
185 err = make_pair(_wsum->err(),_wsum->err());
186 }
187 else if (ix==2) {
188 val = Geff.first;
189 err = make_pair(Geff.second,Geff.second);
190 }
191 else if (ix==3) {
192 val = alpha0.first;
193 err = alpha0.second;
194 }
195 else if (ix==4) {
196 val = R.first;
197 err = R.second;
198 }
199 else {
200 val = GM.first;
201 err = GM.second;
202 }
203 BinnedEstimatePtr<string> tmp;
204 book(tmp, 1, 1, ix);
205 tmp->binAt(_ecms).set(val,err);
206 }
207 }
208
209 /// @}
210
211
212 /// @name Histograms
213 /// @{
214 Histo1DPtr _h_cThetaL;
215 CounterPtr _wsum;
216 string _ecms;
217 /// @}
218
219
220 };
221
222
223 RIVET_DECLARE_PLUGIN(BESIII_2023_I2677290);
224
225}
|