Rivet analyses referencePLUTO_1981_I165122Ratio of the cross section for the production of $K^0_S$ to that for $\mu^+\mu^-$ between $3.6$ and $31.6$ GeVExperiment: PLUTO (DORIS) Inspire ID: 165122 Status: VALIDATED Authors:
Beam energies: ANY Run details:
Ratio of the cross section for the production of $K^0_S$ to that for $\mu^+\mu^-$ between $3.6$ and $31.6$ GeV. The average number of $K^0$ per hadronic event is also provided for some energies, together with the kaon spectrum at 9.4 and 31.6 GeV. N.B. The point at $9.45\to9.456$ GeV is from the $\Upsilon(1S)$ resonance. Beam energy must be specified as analysis option "ENERGY" when rivet-merging samples. Source code: PLUTO_1981_I165122.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 kaon production at low energies
11 class PLUTO_1981_I165122 : public Analysis {
12 public:
13
14 /// Constructor
15 RIVET_DEFAULT_ANALYSIS_CTOR(PLUTO_1981_I165122);
16
17
18 /// @name Analysis methods
19 /// @{
20
21 /// Book histograms and initialise projections before the run
22 void init() {
23
24 // Initialise and register projections
25 declare(Beam(), "Beams");
26 declare(FinalState(), "FS");
27 declare(UnstableParticles(), "UFS");
28
29 // Book histograms
30 book(_c_hadrons , "/TMP/sigma_hadrons");
31 book(_c_muons , "/TMP/sigma_muons");
32 book(_c_kaons , "/TMP/sigma_kaons");
33 book(_c_hadronsY, "/TMP/sigma_hadronsY");
34 book(_c_muonsY , "/TMP/sigma_muonsY");
35 book(_c_kaonsY , "/TMP/sigma_kaonsY");
36 if (isCompatibleWithSqrtS(9.4*GeV)) {
37 book(_h_spectrum1, 5, 1, 1);
38 }
39 else if (isCompatibleWithSqrtS(30.0*GeV, 1E-2)) {
40 book(_h_spectrum1, 4, 1, 1);
41 }
42 book(_h_spectrum2, 6, 1, 1);
43 }
44
45 /// Recursively walk the decay tree to find decay products of @a p
46 void findDecayProducts(Particle mother, Particles &kaons, Particles& stable) {
47 for(const Particle & p: mother.children()) {
48 const int id = p.pid();
49 if(id==130 || id ==310) {
50 kaons.push_back(p);
51 }
52 if (id==111 or p.children().empty())
53 stable.push_back(p);
54 else
55 findDecayProducts(p, kaons, stable);
56 }
57 }
58
59
60 /// Perform the per-event analysis
61 void analyze(const Event& event) {
62 // Get beams and average beam momentum
63 const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
64 const double meanBeamMom = ( beams.first.p3().mod() +
65 beams.second.p3().mod() ) / 2.0;
66 MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
67 // Find the Upsilons among the unstables
68 const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
69 Particles upsilons = ufs.particles(Cuts::pid==553);
70 // Continuum
71 if (upsilons.empty()) {
72 MSG_DEBUG("No Upsilons found => continuum event");
73 // final state particles
74 const FinalState& fs = apply<FinalState>(event, "FS");
75 map<long,int> nCount;
76 int ntotal(0);
77 for (const Particle& p : fs.particles()) {
78 nCount[p.pid()] += 1;
79 ++ntotal;
80 }
81 // mu+mu- + photons
82 if(nCount[-13]==1 and nCount[13]==1 &&
83 ntotal==2+nCount[22])
84 _c_muons->fill();
85 // everything else
86 else
87 _c_hadrons->fill();
88 // unstable particles
89 for (const Particle& p : ufs.particles(Cuts::pid==130 or Cuts::pid==310)) {
90 if(_h_spectrum1) {
91 double xp = p.p3().mod()/meanBeamMom;
92 _h_spectrum1->fill(xp);
93 }
94 _c_kaons->fill();
95 }
96 }
97 else {
98 MSG_DEBUG("Upsilons found => resonance event");
99 for (const Particle& ups : upsilons) {
100 Particles kaons,stable;
101 // Find the decay products we want
102 findDecayProducts(ups, kaons, stable);
103 // boost to rest frame (if required)
104 LorentzTransform cms_boost;
105 if (ups.p3().mod() > 1*MeV)
106 cms_boost = LorentzTransform::mkFrameTransformFromBeta(ups.momentum().betaVec());
107 const double mass = ups.mass();
108
109 map<long,int> nCount;
110 int ntotal(0);
111 for (const Particle& p : stable) {
112 nCount[p.pid()] += 1;
113 ++ntotal;
114 }
115 for( const Particle & kaon : kaons) {
116 const FourMomentum p2 = cms_boost.transform(kaon.momentum());
117 const double xp = 2.*p2.p3().mod()/mass;
118 _h_spectrum2->fill(xp);
119 _c_kaonsY->fill();
120 }
121 // mu+mu- + photons
122 if(nCount[-13]==1 and nCount[13]==1 &&
123 ntotal==2+nCount[22])
124 _c_muonsY->fill();
125 // everything else
126 else
127 _c_hadronsY->fill();
128 }
129 }
130 }
131
132
133 /// Normalise histograms etc., after the run
134 void finalize() {
135 // energy dependent
136 for(unsigned int ix=1;ix<4;++ix) {
137 CounterPtr denom = (ix==1 || ix==3 ) ? _c_muons : _c_hadrons;
138 Estimate0D R = *_c_kaons/ *denom;
139 const double rval = R.val();
140 const double rerr = R.errPos();
141 Estimate1DPtr mult;
142 book(mult,ix, 1, 1);
143 for (auto& b : mult->bins()) {
144 if( (ix==1 || ix==3) && _c_muonsY ->val()==0.) continue;
145 if( ix==2 && _c_hadrons->val()==0.) continue;
146 const double x = b.xMid();
147 if(x==9.458) {
148 if(_c_kaonsY->val()>0.) {
149 if (ix==1 ) {
150 b = *_c_kaonsY/ *_c_muonsY;
151 }
152 else if (ix==2) {
153 b = *_c_kaonsY/ *_c_hadronsY;
154 }
155 else if (ix==3) {
156 b = *_c_kaonsY/ *_c_muonsY;
157 }
158 }
159 }
160 else if (denom->numEntries()>0) {
161 if (inRange(sqrtS()/GeV, b.xMin(), b.xMax())) {
162 b.set(rval, rerr);
163 }
164 }
165 }
166 }
167 // normalize the spectra if required
168 if(_h_spectrum1) {
169 scale(_h_spectrum1, sqr(sqrtS())*crossSection()/microbarn/sumOfWeights());
170 }
171 if(_h_spectrum2) {
172 scale(_h_spectrum2, 1./_c_hadronsY->val());
173 }
174 }
175
176 /// @}
177
178
179 /// @name Histograms
180 /// @{
181 Histo1DPtr _h_spectrum1, _h_spectrum2;
182 CounterPtr _c_hadrons, _c_muons, _c_kaons;
183 CounterPtr _c_hadronsY, _c_muonsY, _c_kaonsY;
184 /// @}
185
186
187 };
188
189
190 RIVET_DECLARE_PLUGIN(PLUTO_1981_I165122);
191
192
193}
|