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)) {
37 book(_h_spectrum1, 5, 1, 1);
38 }
39 else if (isCompatibleWithSqrtS(30.0, 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 Scatter1D R = *_c_kaons/ *denom;
139 double rval = R.point(0).x();
140 pair<double,double> rerr = R.point(0).xErrs();
141 Scatter2D temphisto(refData(ix, 1, 1));
142 Scatter2DPtr mult;
143 book(mult,ix, 1, 1);
144 for (size_t b = 0; b < temphisto.numPoints(); b++) {
145 if( (ix==1 || ix==3) && _c_muonsY ->val()==0.) continue;
146 if( ix==2 && _c_hadrons->val()==0.) continue;
147 const double x = temphisto.point(b).x();
148 pair<double,double> ex = temphisto.point(b).xErrs();
149 if(x==9.458) {
150 if(_c_kaonsY->val()>0.) {
151 Scatter1D R2;
152 if(ix==1 ) {
153 R2 = *_c_kaonsY/ *_c_muonsY;
154 }
155 else if(ix==2) {
156 R2 = *_c_kaonsY/ *_c_hadronsY;
157 }
158 else if(ix==3) {
159 R2 = *_c_kaonsY/ *_c_muonsY;
160 }
161 mult ->addPoint(x, R2.point(0).x(), ex, R2.point(0).xErrs());
162 }
163 else {
164 mult ->addPoint(x, 0., ex, make_pair(0.,.0));
165 }
166 }
167 else if (denom->numEntries()>0) {
168 pair<double,double> ex2 = ex;
169 if(ex2.first ==0.) ex2. first=0.0001;
170 if(ex2.second==0.) ex2.second=0.0001;
171
172 if (inRange(sqrtS()/GeV, x-ex2.first, x+ex2.second)) {
173 mult ->addPoint(x, rval, ex, rerr);
174 }
175 else {
176 mult ->addPoint(x, 0., ex, make_pair(0.,.0));
177 }
178 }
179 }
180 }
181 // normalize the spectra if required
182 if(_h_spectrum1) {
183 scale(_h_spectrum1, sqr(sqrtS())*crossSection()/microbarn/sumOfWeights());
184 }
185 if(_h_spectrum2) {
186 scale(_h_spectrum2, 1./_c_hadronsY->val());
187 }
188 }
189
190 //@}
191
192
193 /// @name Histograms
194 //@{
195 Histo1DPtr _h_spectrum1, _h_spectrum2;
196 CounterPtr _c_hadrons, _c_muons, _c_kaons;
197 CounterPtr _c_hadronsY, _c_muonsY, _c_kaonsY;
198 //@}
199
200
201 };
202
203
204 // The hook for the plugin system
205 RIVET_DECLARE_PLUGIN(PLUTO_1981_I165122);
206
207
208}
|