Rivet analyses referenceBESIII_2019_I1725786$e^+e^-\to D\bar{D}\pi^+\pi^-$ for $\sqrt{s}=4.0854$ to $4.6$ GeVExperiment: BESIII (BEPC) Inspire ID: 1725786 Status: VALIDATED NOHEPDATA SINGLEWEIGHT Authors:
Beam energies: (2.0, 2.0); (2.1, 2.1); (2.1, 2.1); (2.1, 2.1); (2.1, 2.1); (2.1, 2.1); (2.1, 2.1); (2.2, 2.2); (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:
Analysis of the resonant contributions to $e^+e^-\to D\bar{D}\pi^+\pi^-$ for $\sqrt{s}=4.0854$ to $4.6$ GeV, where $D\bar{D}=D^+D^-, D^0\bar{D}^0$. The cross sections for the production of $\pi^+\pi^-\psi(3770)$ and $D_1(2420^0\bar{D}^0$+c.c. are measured as a function of energy and the mass distrinutions are implemented for three energies. Source code: BESIII_2019_I1725786.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- > D D pi pi
10 class BESIII_2019_I1725786 : public Analysis {
11 public:
12
13 /// Constructor
14 RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2019_I1725786);
15
16
17 /// @name Analysis methods
18 /// @{
19
20 /// Book histograms and initialise projections before the run
21 void init() {
22 // projections
23 declare(FinalState(), "FS");
24 // histograms
25 unsigned int iloc=0;
26 if (isCompatibleWithSqrtS(4.4156)) {
27 iloc=2;
28 for (unsigned int ix=0; ix<2; ++ix) {
29 book(_h_Dpi[ix],1,1,1+ix);
30 }
31 book (_h_DD,2,1,3);
32 }
33 else if (isCompatibleWithSqrtS(4.258)) {
34 book(_h_DD,2,1,1);
35 }
36 else if (isCompatibleWithSqrtS(4.3583)) {
37 book(_h_DD,2,1,2);
38 iloc=1;
39 }
40 else if (isCompatibleWithSqrtS(4.5995)) {
41 iloc=3;
42 }
43 if (iloc>0) {
44 for (unsigned int iy=0;iy<3;++iy) {
45 book(_h_mass[iy],3,iloc,1+iy);
46 }
47 }
48 iloc=0;
49 for (unsigned int ix=0;ix<5;++ix) {
50 if (ix==1) continue;
51 book(_sigma[ix], 4+ix,1,1);
52 if (ix!=0 && ix!=2 ) continue;
53 for (const string& en : _sigma[ix].binning().edges<0>()) {
54 const double end = std::stod(en)*GeV;
55 if (isCompatibleWithSqrtS(end)) {
56 _ecms[iloc] = en;
57 break;
58 }
59 }
60 iloc+=1;
61 }
62 if(_ecms[0].empty()) MSG_ERROR("Beam energy incompatible with analysis.");
63 }
64
65
66 /// Perform the per-event analysis
67 void analyze(const Event& event) {
68 Particles fs = apply<FinalState>(event, "FS").particles();
69 Particles DD,other;
70 for (const Particle & p : fs) {
71 Particle parent=p;
72 while(!parent.parents().empty()) {
73 parent=parent.parents()[0];
74 if(parent.abspid()==PID::D0 ||
75 parent.abspid()==PID::DPLUS) break;
76 }
77 if(parent.abspid()!=PID::D0 &&
78 parent.abspid()!=PID::DPLUS) {
79 other.push_back(p);
80 continue;
81 }
82 bool found=false;
83 for (auto& D : DD) {
84 // D already in list
85 if (fuzzyEquals(D.momentum(),parent.momentum())) {
86 found=true;
87 break;
88 }
89 }
90 if(!found) DD.push_back(parent);
91 }
92 // d dbar + 2 other particles
93 if(DD.size()!=2 || other.size()!=2) vetoEvent;
94 // other particles pi+ pi-
95 if(!(other[0].pid()==-other[1].pid() && other[0].abspid()==PID::PIPLUS)) vetoEvent;
96 // d dbar pair
97 if(DD[0].pid()!=-DD[1].pid()) vetoEvent;
98 if(other[0].pid()<0) swap(other[0],other[1]);
99 if(DD [0].pid()<0) swap(DD [0],DD [1]);
100 // fill the mass plots if needed
101 if (_h_Dpi[0] && DD[0].pid()==421) {
102 _h_Dpi[0]->fill((DD[0].momentum()+other[0].momentum()).mass());
103 _h_Dpi[1]->fill((DD[1].momentum()+other[1].momentum()).mass());
104 }
105 if (_h_DD) {
106 _h_DD->fill((DD[0].momentum()+DD[1].momentum()).mass());
107 }
108 const bool DstarP = (DD[0].pid()== 421 && DD[0].parents()[0].pid()== 413);
109 const bool DstarM = (DD[1].pid()==-421 && DD[1].parents()[0].pid()==-413);
110 const bool Dstar = DstarP || DstarM;
111 // D0 pi+ pi- mass
112 if (_h_mass[0] && DD[0].pid()==421&&!Dstar) {
113 _h_mass[0]->fill((DD[0].momentum()+other[0].momentum()+other[1].momentum()).mass());
114 }
115 // D*+ pi- mass
116 if (_h_mass[1] && DstarP) {
117 _h_mass[1]->fill((DD[0].momentum()+other[0].momentum()+other[1].momentum()).mass());
118 }
119 // D+ pi+ pi- mass
120 if (_h_mass[2] && DD[0].pid()==411) {
121 _h_mass[2]->fill((DD[0].momentum()+other[0].momentum()+other[1].momentum()).mass());
122 }
123 // now for the counters for the cross sections
124 if (DD[0].parents()[0].pid()==30443 && DD[1].parents()[0].pid()==30443) {
125 _sigma[0]->fill(_ecms[0]);
126 }
127 else if (DD[0].pid()==421 && (DD[0].parents()[0].pid()==10423 || DD[1].parents()[0].pid()==-10423)) {
128 _sigma[2]->fill(_ecms[1]);
129 }
130 else if (DD[0].pid()==421 && ((DstarP && DD[0].parents()[0].parents()[0].pid()== 10423) ||
131 (DstarM && DD[1].parents()[0].parents()[0].pid()==-10423))) {
132 _sigma[3]->fill(_ecms[1]);
133 }
134 else if (DD[0].pid()==411 && (DD[0].parents()[0].pid()==10413 || DD[1].parents()[0].pid()==-10413)) {
135 _sigma[4]->fill(_ecms[1]);
136 }
137 }
138
139
140 /// Normalise histograms etc., after the run
141 void finalize() {
142 // dists
143 if (_h_Dpi[0]) normalize(_h_Dpi, 1.0, false);
144 if (_h_DD) normalize(_h_DD, 1.0, false);
145 if (_h_mass[0]) {
146 for(unsigned int ix=0;ix<3;++ix)
147 normalize(_h_mass[ix], 1.0, false);
148 }
149 // cross sections
150 const vector<double> brs = { 0.93, 1., 1., 0.677, 1.};
151 const double fact = crossSection()/ sumOfWeights() /picobarn;;
152 for (unsigned int ii=0; ii<5; ++ii) {
153 if (ii==1) continue;
154 scale(_sigma[ii], fact/brs[ii]);
155 }
156 }
157
158 /// @}
159
160
161 /// @name Histograms
162 /// @{
163 Histo1DPtr _h_Dpi[2],_h_DD,_h_mass[3];
164 BinnedHistoPtr<string> _sigma[5];
165 string _ecms[2];
166 /// @}
167
168
169 };
170
171
172 RIVET_DECLARE_PLUGIN(BESIII_2019_I1725786);
173
174}
|