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