Rivet analyses referenceBELLE_2017_I1590028$e^+e^-\to J/\psi D\bar{D}$ at $\sqrt{s}=10.58$GeVExperiment: BELLE (KEKB) Inspire ID: 1590028 Status: VALIDATED NOHEPDATA Authors:
Beam energies: (4.7, 4.7); (5.0, 5.0); (5.2, 5.2); (5.3, 5.3); (5.3, 5.3); (5.4, 5.4) GeV Run details:
Measurement of mass and angular distributions in $e^+e^-\to J/\psi D\bar{D}$ at $\sqrt{s}=10.58$GeV by BELLE. The cross section for $e^+e^-\to J/\psi X^*(3860)(\to D\bar{D}$) is also measured for production at the energies of the $\Upsilon(1\to5S)$ states and in the continuum at 10.52GeV. The status of the $X^*(3860)$ is not clear, although it is characterized as a $\chi_{c0}$ state we therefore take the PDG code to be 9010441, which can be changed using the PID option Source code: BELLE_2017_I1590028.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/Beam.hh"
4#include "Rivet/Projections/FinalState.hh"
5
6namespace Rivet {
7
8
9 /// @brief e+ e- > J/psi D Dbar
10 class BELLE_2017_I1590028 : public Analysis {
11 public:
12
13 /// Constructor
14 RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2017_I1590028);
15
16
17 /// @name Analysis methods
18 /// @{
19
20 /// Book histograms and initialise projections before the run
21 void init() {
22 // set the PDG code
23 _pid = getOption<int>("PID", 9010441);
24 // projections
25 declare(Beam(), "Beams");
26 declare(FinalState(), "FS");
27 // histograms
28 if(isCompatibleWithSqrtS(10.58,1e-4)) {
29 // histograms
30 for(unsigned int ix=0;ix<6;++ix) {
31 book(_h[ix],1,1,1+ix);
32 }
33 }
34 book(_sigma,2,1,1);
35 for (const string& en : _sigma.binning().edges<0>()) {
36 const double end = std::stod(en)*GeV;
37 if (isCompatibleWithSqrtS(end,1e-4)) {
38 _ecms = en;
39 break;
40 }
41 }
42 if (_ecms.empty()) MSG_ERROR("Beam energy incompatible with analysis.");
43 }
44
45
46 /// Perform the per-event analysis
47 void analyze(const Event& event) {
48 // get the axis, direction of incoming electron
49 const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
50 Vector3 axis;
51 FourMomentum pep;
52 if (beams.first.pid()>0) {
53 axis = beams.first .mom().p3().unit();
54 pep = beams.second.mom();
55 }
56 else {
57 axis = beams.second.mom().p3().unit();
58 pep = beams.first .mom();
59 }
60 Particles fs = apply<FinalState>(event, "FS").particles();
61 Particles DD,JPSI,other;
62 for (const Particle& p : fs) {
63 Particle parent=p;
64 while (!parent.parents().empty()) {
65 parent=parent.parents()[0];
66 if (parent.abspid()==PID::DPLUS ||
67 parent.abspid()==PID::D0||
68 parent.pid() ==PID::JPSI) break;
69 }
70 if (parent.abspid()!=PID::DPLUS &&
71 parent.abspid()!=PID::D0 &&
72 parent.pid() !=PID::JPSI) {
73 other.push_back(p);
74 continue;
75 }
76 bool found=false;
77 for (const auto& D : parent.pid()==PID::JPSI ? JPSI : DD) {
78 // D or J/psi already in list
79 if (fuzzyEquals(D.mom(),parent.mom())) {
80 found=true;
81 break;
82 }
83 }
84 if (!found) {
85 (parent.pid()==PID::JPSI ? JPSI : DD).push_back(parent);
86 }
87 }
88 // require J/psi D Dbar
89 if (JPSI.size()!=1 || DD.size()!=2 || other.size()!=0) vetoEvent;
90 if (DD[0].pid()!=-DD[1].pid()) vetoEvent;
91 FourMomentum pDD = DD[0].mom()+DD[1].mom();
92 // cross section
93 if(DD[0].parents()[0].pid()==_pid && DD[1].parents()[0].pid()==_pid) {
94 _sigma->fill(_ecms);
95 }
96 // histograms only at 10.58GeV
97 if(!_h[0]) return;
98 _h[0]->fill(pDD.mass());
99 const double cProd = axis.dot(JPSI[0].mom().p3().unit());
100 _h[1]->fill(cProd);
101 // finally the leptons from J/psi decay
102 if (JPSI[0].children().size()!=2) vetoEvent;
103 if (JPSI[0].children()[0].pid()!=-JPSI[0].children()[1].pid()) vetoEvent;
104 if (JPSI[0].children()[0].abspid()!=PID::EMINUS &&
105 JPSI[0].children()[0].abspid()!=PID::MUON) vetoEvent;
106 Particle lm = JPSI[0].children()[0];
107 Particle lp = JPSI[0].children()[1];
108 // variables in J/psi frame
109 Vector3 axis2 = -pDD.p3().unit();
110 LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(JPSI[0].mom().betaVec());
111 FourMomentum plm = boost1.transform(lm.mom());
112 Vector3 axis4 = plm.p3().unit();
113 double cPsi = axis2.dot(axis4);
114 _h[2]->fill(cPsi);
115 Vector3 axis3 = boost1.transform(pep).p3().unit();
116 Vector3 aep = axis3-axis3.dot(axis2)*axis2;
117 Vector3 alm = axis4-axis4.dot(axis2)*axis2;
118 double phil = atan2(-axis2.cross(aep).dot(alm),aep.dot(alm));
119 _h[4]->fill(phil);
120 // variables in DD frame
121 LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(pDD.betaVec());
122 FourMomentum pD = boost2.transform((DD[0].pid()>0 ? DD[0] : DD[1]).mom());
123 Vector3 axisD = pD.p3().unit();
124 axis2 *= -1;
125 double cX = axis2.dot(axisD);
126 _h[3]->fill(cX);
127 axis3 = boost2.transform(pep).p3().unit();
128 aep = axis3-axis3.dot(axis2)*axis2;
129 alm = axisD-axisD.dot(axis2)*axis2;
130 phil = atan2(-axis2.cross(aep).dot(alm),aep.dot(alm));
131 _h[5]->fill(phil);
132 }
133
134
135 /// Normalise histograms etc., after the run
136 void finalize() {
137 normalize(_h, 1.0, false);
138 scale(_sigma,crossSection()/ sumOfWeights() /femtobarn);
139 }
140
141 /// @}
142
143
144 /// @name Histograms
145 /// @{
146 int _pid;
147 Histo1DPtr _h[6];
148 BinnedHistoPtr<string> _sigma;
149 string _ecms;
150 /// @}
151
152 };
153
154
155 RIVET_DECLARE_PLUGIN(BELLE_2017_I1590028);
156
157}
|