Rivet analyses referenceBESIII_2015_I1391798$e^+e^-\to (D\bar{D}^*)^\pm\pi^\mp$ for $\sqrt{s}=4.23$ and 4.26 GeVExperiment: BESIII (BEPC) Inspire ID: 1391798 Status: VALIDATED NOHEPDATA SINGLEWEIGHT Authors:
Beam energies: (2.1, 2.1); (2.1, 2.1) GeV Run details:
Measurement of mass distributions for $e^+e^-\to (D\bar{D}^*)^\pm\pi^\mp$ for $\sqrt{s}=4.23$ and 4.26 GeV by BES. The cross section for $e^+e^-\to Z_c(3885)^\pm\pi^\mp\to (D\bar{D}^*)^\pm\pi^\mp$ is also measured. As there is no PDG code for the $Z_c(3885)^+ we take it to be 9044213, although this can be changed using the PID option. Source code: BESIII_2015_I1391798.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FinalState.hh"
4#include "Rivet/Projections/UnstableParticles.hh"
5#include "Rivet/Projections/Beam.hh"
6
7namespace Rivet {
8
9
10 /// @brief e+e-> Z+- pi-+
11 class BESIII_2015_I1391798 : public Analysis {
12 public:
13
14 /// Constructor
15 RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2015_I1391798);
16
17
18 /// @name Analysis methods
19 /// @{
20
21 /// Book histograms and initialise projections before the run
22 void init() {
23 // set the PDG code
24 _pid = getOption<double>("PID", 9044213);
25 // projections
26 declare(Beam(), "Beams");
27 declare(FinalState(), "FS");
28 // histograms
29 unsigned int iEnergy=0;
30 if (isCompatibleWithSqrtS(4.23)) {
31 _ecms="4.23";
32 iEnergy=1;
33 }
34 else if (isCompatibleWithSqrtS(4.26)) {
35 _ecms="4.26";
36 iEnergy=2;
37 }
38 else {
39 MSG_ERROR("Beam energy incompatible with analysis.");
40 }
41 // histos
42 for (unsigned int ix=0;ix<2;++ix) {
43 book(_h[ix],2,ix+1,iEnergy);
44 book(_h[2+ix],3,1,1+ix);
45 book(_sigma[ix],1,1,1+ix);
46 }
47 }
48
49
50 /// Perform the per-event analysis
51 void analyze(const Event& event) {
52 // get the axis, direction of incoming electron
53 const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
54 Vector3 axis;
55 if (beams.first.pid()>0) {
56 axis = beams.first .momentum().p3().unit();
57 }
58 else {
59 axis = beams.second.momentum().p3().unit();
60 }
61 Particles fs = apply<FinalState>(event, "FS").particles();
62 Particles DD,other;
63 for (const Particle& p : fs) {
64 Particle parent=p;
65 while (!parent.parents().empty()) {
66 parent=parent.parents()[0];
67 if (parent.abspid()==411 || parent.abspid()==413 ||
68 parent.abspid()==421 || parent.abspid()==423) break;
69 }
70 if ((parent.abspid()==411 || parent.abspid()==421)
71 && !parent.parents().empty()) {
72 Particle Dstar = parent.parents()[0];
73 if (Dstar.abspid()==413 || Dstar.abspid()==423) {
74 parent=Dstar;
75 }
76 }
77 if (parent.abspid()==411 || parent.abspid()==413 ||
78 parent.abspid()==421 || parent.abspid()==423) {
79 bool found=false;
80 for (const auto& D : DD) {
81 // D already in list
82 if (fuzzyEquals(D.mom(), parent.mom())) {
83 found=true;
84 break;
85 }
86 }
87 if(!found) DD.push_back(parent);
88 }
89 else {
90 other.push_back(p);
91 }
92 }
93 // D Dbar + charged pion
94 if (DD.size()!=2 || other.size()!=1) vetoEvent;
95 if (DD[0].pid()*DD[1].pid()>0) vetoEvent;
96 if (other[0].abspid()!=211) vetoEvent;
97 if (DD[0].abspid()%10!=3) swap(DD[0],DD[1]);
98 if (DD[0].abspid()%10!=3 || DD[1].abspid()%10!=1) vetoEvent;
99 const double mass = (DD[0].momentum()+DD[1].momentum()).mass();
100 const double cTheta = abs(axis.dot(other[0].momentum().p3().unit()));
101 unsigned int iloc=0;
102 // D0 D*- pi+ +cc
103 if ((DD[0].pid()==-413 && DD[1].pid()== 421 && other[0].pid()== 211) ||
104 (DD[0].pid()== 413 && DD[1].pid()==-421 && other[0].pid()==-211)) {
105 iloc=0;
106 }
107 // D- D*0 pi+ +cc
108 else if((DD[0].pid()== 423 && DD[1].pid()==-411 && other[0].pid()== 211) ||
109 (DD[0].pid()==-423 && DD[1].pid()== 411 && other[0].pid()==-211)) {
110 iloc=1;
111 }
112 // otherwise veto event
113 else {
114 vetoEvent;
115 }
116 _h[iloc]->fill(mass);
117 // parent Z+/-
118 if (DD[0].parents()[0].abspid()==_pid && DD[1].parents()[0].abspid()==_pid &&
119 fuzzyEquals(DD[0].parents()[0].momentum(),DD[1].parents()[0].momentum()) ) {
120 _sigma[0]->fill(_ecms);
121 if (_ecms=="4.26") _sigma[1]->fill(_ecms);
122 _h[2+iloc]->fill(cTheta);
123 }
124 }
125
126 /// Normalise histograms etc., after the run
127 void finalize() {
128 // distributions
129 normalize(_h, 1.0, false);
130 // cross section
131 for(unsigned int ix=0;ix<2;++ix)
132 scale(_sigma[ix], crossSection()/ sumOfWeights() /picobarn);
133 }
134
135 /// @}
136
137
138 /// @name Histograms
139 /// @{
140 BinnedHistoPtr<string> _sigma[2];
141 Histo1DPtr _h[4];
142 int _pid;
143 string _ecms;
144 /// @}
145
146
147 };
148
149
150 RIVET_DECLARE_PLUGIN(BESIII_2015_I1391798);
151
152}
|
