rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

PLUTO_1981_I165122

Ratio of the cross section for the production of $K^0_S$ to that for $\mu^+\mu^-$ between $3.6$ and $31.6$ GeV
Experiment: PLUTO (DORIS)
Inspire ID: 165122
Status: VALIDATED
Authors:
  • Peter Richardson
References:
  • Phys.Lett. B104 (1981) 79-83, 1981
Beams: e- e+
Beam energies: ANY
Run details:
  • e+ e- to hadrons and e+ e- to mu+ mu- (for normalization) Beam energy must be specified as analysis option "ENERGY" when rivet-merging samples.

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}