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 // // Book histograms
29 book(_c_hadrons , "/TMP/sigma_hadrons", refData<YODA::BinnedEstimate<string>>(1, 1, 1));
30 for(unsigned int ix=0;ix<2;++ix) {
31 book(_c_muons[ix], "/TMP/sigma_muons_"+toString(ix), refData<YODA::BinnedEstimate<string>>(1+2*ix, 1, 1));
32 book(_c_kaons[ix], "/TMP/sigma_kaons_"+toString(ix), refData<YODA::BinnedEstimate<string>>(1+2*ix, 1, 1));
33 for (const string& en : _c_muons[ix].binning().edges<0>()) {
34 // these points are really Upsilon(1S) decay
35 if(en == "9.45 - 9.466"s || en=="9.45 - 9.456"s) continue;
36 const size_t idx = en.find("-");
37 if(idx!=std::string::npos) {
38 const double emin = std::stod(en.substr(0,idx));
39 const double emax = std::stod(en.substr(idx+1,string::npos));
40 if(inRange(sqrtS()/GeV, emin, emax)) {
41 _ecms[ix] = en;
42 break;
43 }
44 }
45 else {
46 const double end = std::stod(en)*GeV;
47 if (isCompatibleWithSqrtS(end)) {
48 _ecms[ix] = en;
49 break;
50 }
51 }
52 }
53 }
54 if (isCompatibleWithSqrtS(9.4*GeV)) {
55 book(_h_spectrum1, 5, 1, 1);
56 }
57 else if (isCompatibleWithSqrtS(30.0*GeV, 1E-2)) {
58 book(_h_spectrum1, 4, 1, 1);
59 }
60 book(_h_spectrum2, 6, 1, 1);
61 book(_c_hadronsY,"TMP/nUps");
62 }
63
64 /// Recursively walk the decay tree to find decay products of @a p
65 void findDecayProducts(Particle mother, Particles &kaons, Particles& stable) {
66 for(const Particle & p: mother.children()) {
67 const int id = p.pid();
68 if(id==130 || id ==310) {
69 kaons.push_back(p);
70 }
71 if (id==111 or p.children().empty())
72 stable.push_back(p);
73 else
74 findDecayProducts(p, kaons, stable);
75 }
76 }
77
78
79 /// Perform the per-event analysis
80 void analyze(const Event& event) {
81 // Get beams and average beam momentum
82 const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
83 const double meanBeamMom = ( beams.first.p3().mod() +
84 beams.second.p3().mod() ) / 2.0;
85 MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
86 // Find the Upsilons among the unstables
87 const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
88 Particles upsilons = ufs.particles(Cuts::pid==553);
89 // Continuum
90 if (upsilons.empty()) {
91 MSG_DEBUG("No Upsilons found => continuum event");
92 // final state particles
93 const FinalState& fs = apply<FinalState>(event, "FS");
94 map<long,int> nCount;
95 int ntotal(0);
96 for (const Particle& p : fs.particles()) {
97 nCount[p.pid()] += 1;
98 ++ntotal;
99 }
100 // mu+mu- + photons
101 if(nCount[-13]==1 and nCount[13]==1 &&
102 ntotal==2+nCount[22]) {
103 for(unsigned int ix=0;ix<2;++ix) _c_muons[ix]->fill(_ecms[ix]);
104 }
105 // everything else
106 else
107 _c_hadrons->fill(_ecms[0]);
108 // unstable particles
109 for (const Particle& p : ufs.particles(Cuts::pid==130 or Cuts::pid==310)) {
110 if(_h_spectrum1) {
111 double xp = p.p3().mod()/meanBeamMom;
112 _h_spectrum1->fill(xp);
113 }
114 for(unsigned int ix=0;ix<2;++ix) _c_kaons[ix]->fill(_ecms[ix]);
115 }
116 }
117 else {
118 MSG_DEBUG("Upsilons found => resonance event");
119 for (const Particle& ups : upsilons) {
120 Particles kaons,stable;
121 // Find the decay products we want
122 findDecayProducts(ups, kaons, stable);
123 // boost to rest frame (if required)
124 LorentzTransform cms_boost;
125 if (ups.p3().mod() > 1*MeV)
126 cms_boost = LorentzTransform::mkFrameTransformFromBeta(ups.momentum().betaVec());
127 const double mass = ups.mass();
128
129 map<long,int> nCount;
130 int ntotal(0);
131 for (const Particle& p : stable) {
132 nCount[p.pid()] += 1;
133 ++ntotal;
134 }
135 for( const Particle & kaon : kaons) {
136 const FourMomentum p2 = cms_boost.transform(kaon.momentum());
137 const double xp = 2.*p2.p3().mod()/mass;
138 _h_spectrum2->fill(xp);
139 _c_kaons[0]->fill("9.45 - 9.466"s);
140 _c_kaons[1]->fill("9.45 - 9.456"s);
141 }
142 // mu+mu- + photons
143 if(nCount[-13]==1 and nCount[13]==1 &&
144 ntotal==2+nCount[22]) {
145 _c_muons[0]->fill("9.45 - 9.466"s);
146 _c_muons[1]->fill("9.45 - 9.456"s);
147 }
148 // everything else
149 else {
150 _c_hadrons->fill("9.45 - 9.466"s);
151 _c_hadronsY->fill();
152 }
153 }
154 }
155 }
156
157
158 /// Normalise histograms etc., after the run
159 void finalize() {
160 BinnedEstimatePtr<string> ratio;
161 book(ratio,1,1,1);
162 divide(_c_kaons[0],_c_muons[0],ratio);
163 book(ratio,2,1,1);
164 divide(_c_kaons[0],_c_hadrons,ratio);
165 book(ratio,3,1,1);
166 divide(_c_kaons[1],_c_muons[1],ratio);
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);
173 }
174 }
175
176 /// @}
177
178
179 /// @name Histograms
180 /// @{
181 BinnedHistoPtr<string> _c_hadrons,_c_muons[2],_c_kaons[2];
182 string _ecms[2];
183 Histo1DPtr _h_spectrum1, _h_spectrum2;
184 CounterPtr _c_hadronsY;
185 /// @}
186
187
188 };
189
190
191 RIVET_DECLARE_PLUGIN(PLUTO_1981_I165122);
192
193
194}
|