rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

CLEOC_2008_I777917

Measurement of charm final-states, the total hadronic cross section and $R$ for energies between 3.92 and 4.26 GeV
Experiment: CLEOC (CESR)
Inspire ID: 777917
Status: VALIDATED
Authors:
  • Peter Richardson
References:
  • Phys.Rev. D80 (2009) 072001, 2009
Beams: e- e+
Beam energies: ANY
Run details:
  • e+ e- to hadrons and e+ e- to mu+ mu- (for normalization)

Measurement of charm final-states, the total hadronic cross section and $R$ for energies between 3.92 and 4.26 GeV The muonic cross section is also outputted to the yoda file so that ratio $R$ can be recalcuated if runs are combined.

Source code: CLEOC_2008_I777917.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/UnstableParticles.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief Add a short analysis description here
 10  class CLEOC_2008_I777917 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(CLEOC_2008_I777917);
 15
 16
 17    /// @name Analysis methods
 18    //@{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22
 23      // Initialise and register projections
 24      declare(FinalState(), "FS");
 25      declare(UnstableParticles(), "UFS");
 26      book(_c_hadrons, "/TMP/sigma_hadrons");
 27      book(_c_muons, "/TMP/sigma_muons");
 28      book(_c_D0D0, "/TMP/sigma_D0D0");
 29      book(_c_DpDm, "/TMP/sigma_DpDm");
 30      book(_c_DsDs, "/TMP/sigma_DsDs");
 31      book(_c_D0D0S, "/TMP/sigma_D0D0S");
 32      book(_c_DpDmS, "/TMP/sigma_DpDmS");
 33      book(_c_DsDsS, "/TMP/sigma_DsDsS");
 34      book(_c_D0SD0S, "/TMP/sigma_D0SD0S");
 35      book(_c_DpSDmS, "/TMP/sigma_DpSDmS");
 36      book(_c_DsSDsS, "/TMP/sigma_DsSDsS");
 37      book(_c_DD, "/TMP/sigma_DD");
 38      book(_c_DDX, "/TMP/sigma_DDX");
 39      book(_c_DSDpi, "/TMP/sigma_DSDpi");
 40      book(_c_DSDSpi, "/TMP/sigma_DSDSpi");
 41    }
 42
 43    void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
 44      for (const Particle &child : p.children()) {
 45	if(child.children().empty()) {
 46	  nRes[child.pid()]-=1;
 47	  --ncount;
 48	}
 49	else
 50	  findChildren(child,nRes,ncount);
 51      }
 52    }
 53
 54
 55    /// Perform the per-event analysis
 56    void analyze(const Event& event) {
 57      const FinalState& fs = apply<FinalState>(event, "FS");
 58      // total hadronic and muonic cross sections
 59      map<long,int> nCount;
 60      int ntotal(0);
 61      for (const Particle& p : fs.particles()) {
 62	nCount[p.pid()] += 1;
 63	++ntotal;
 64      }
 65      // mu+mu- + photons
 66      if(nCount[-13]==1 and nCount[13]==1 &&
 67	 ntotal==2+nCount[22])
 68	_c_muons->fill();
 69      // everything else
 70      else
 71	_c_hadrons->fill();
 72      // identified final state with D mesons
 73      const FinalState& ufs = apply<UnstableParticles>(event, "UFS");
 74      for(unsigned int ix=0;ix<ufs.particles().size();++ix) {
 75	const Particle& p1 = ufs.particles()[ix];
 76	int id1 = abs(p1.pid());
 77	if(id1 != 411 && id1 != 413 && id1 != 421 && id1 != 423 &&
 78	   id1 != 431 && id1 != 433)
 79	  continue;
 80	// check fs
 81	bool fs = true;
 82	for (const Particle & child : p1.children()) {
 83	  if(child.pid()==p1.pid()) {
 84	    fs = false;
 85	    break;
 86	  }
 87	}
 88	if(!fs) continue;
 89	// find the children
 90	map<long,int> nRes = nCount;
 91	int ncount = ntotal;
 92	findChildren(p1,nRes,ncount);
 93	bool matched=false;
 94	int sign = p1.pid()/id1;
 95	// loop over the other fs particles
 96	for(unsigned int iy=ix+1;iy<ufs.particles().size();++iy) {
 97	  const Particle& p2 = ufs.particles()[iy];
 98	  fs = true;
 99	  for (const Particle & child : p2.children()) {
100	    if(child.pid()==p2.pid()) {
101	      fs = false;
102	      break;
103	    }
104	  }
105	  if(!fs) continue;
106	  if(p2.pid()/abs(p2.pid())==sign) continue;
107	  int id2 = abs(p2.pid());
108	  if(id2 != 411 && id2 != 413 && id2 != 421 && id2 != 423 &&
109	     id2 != 431 && id2 != 433)
110	    continue;
111	  if(!p2.parents().empty() && p2.parents()[0].pid()==p1.pid())
112	    continue;
113	  if((id1==411 || id1==421 || id1==431) && (id2==411 || id2==421 || id2==431 ))
114	    _c_DDX->fill();
115	  map<long,int> nRes2 = nRes;
116	  int ncount2 = ncount;
117	  findChildren(p2,nRes2,ncount2);
118	  if(ncount2==0) {
119	    matched=true;
120	    for(auto const & val : nRes2) {
121	      if(val.second!=0) {
122		matched = false;
123		break;
124	      }
125	    }
126	    if(matched) {
127	      if(id1==411 && id2==411) {
128		_c_DpDm->fill();
129		_c_DD  ->fill();
130	      }
131	      else if(id1==421&& id2==421) {
132		_c_D0D0->fill();
133		_c_DD  ->fill();
134	      }
135	      else if(id1==431&& id2==431) {
136		_c_DsDs->fill();
137	      }
138	      else if(id1==413 && id2==413) {
139		_c_DpSDmS->fill();
140	      }
141	      else if(id1==423&& id2==423) {
142		_c_D0SD0S->fill();
143	      }
144	      else if(id1==433&& id2==433) {
145		_c_DsSDsS->fill();
146	      }
147	      else if((id1==421 && id2==423) ||
148		      (id1==423 && id2==421)) {
149		_c_D0D0S->fill();
150	      }
151	      else if((id1==411 && id2==413) ||
152		      (id1==413 && id2==411)) {
153		_c_DpDmS->fill();
154	      }
155	      else if((id1==431 && id2==433) ||
156		      (id1==433 && id2==431)) {
157		_c_DsDsS->fill();
158	      }
159	    }
160	  }
161	  else if(ncount2==1) {
162	    int ipi=0;
163	    if(nRes2[111]==1 && nRes2[211]==0 && nRes[-211]==0 )
164	      ipi = 111;
165	    else if(nRes2[111]==0 && nRes2[211]==1 && nRes[-211]==0 )
166	      ipi = 211;
167	    else if(nRes2[111]==0 && nRes2[211]==0 && nRes[-211]==1 )
168	      ipi =-211;
169	    if(ipi==0) continue;
170	    matched=true;
171	    for(auto const & val : nRes2) {
172	      if(val.first==ipi)
173		continue;
174	      else if(val.second!=0) {
175		matched = false;
176		break;
177	      }
178	    }
179	    if(matched) {
180	      bool Ddecay = false;
181	      Particle mother = p1;
182	      while (!mother.parents().empty()) {
183		mother = mother.parents()[0];
184		if(PID::isCharmMeson(mother.pid()) && mother.pid()!=p1.pid()) {
185		  Ddecay = true;
186		  break;
187		}
188	      }
189	      mother = p2;
190	      while (!mother.parents().empty()) {
191		mother = mother.parents()[0];
192		if(PID::isCharmMeson(mother.pid()) && mother.pid()!=p1.pid()) {
193		  Ddecay = true;
194		  break;
195		}
196	      }
197	      if(Ddecay) continue;
198	      if((id1==413 || id1==423 ) &&
199		 (id2==413 || id2==423 )) {
200		_c_DSDSpi->fill();
201	      }
202	      else if((id1==411 || id1==421 ) &&
203		      (id2==413 || id2==423 )) {
204		_c_DSDpi->fill();
205	      }
206	      else if((id1==413 || id1==423 ) &&
207		      (id2==411 || id2==421 )) {
208		_c_DSDpi->fill();
209	      }
210	    }
211	  }
212	}
213      }
214    }
215
216    /// Normalise histograms etc., after the run
217    void finalize() {
218      // R
219      Scatter1D R = *_c_hadrons/ *_c_muons;
220      double              rval = R.point(0).x();
221      pair<double,double> rerr = R.point(0).xErrs();
222      double fact = crossSection()/ sumOfWeights() /nanobarn;
223      double sig_h = _c_hadrons->val()*fact;
224      double err_h = _c_hadrons->err()*fact;
225      double sig_c = _c_DDX->val()*fact;
226      double err_c = _c_DDX->err()*fact;
227      double sig_m = _c_muons  ->val()*fact;
228      double err_m = _c_muons  ->err()*fact;
229      Scatter2D temphisto(refData(6, 1, 2));
230      Scatter2DPtr charm;
231      book(charm, 6,1,1);
232      Scatter2DPtr hadrons;
233      book(hadrons, "sigma_hadrons");
234      Scatter2DPtr muons;
235      book(muons, "sigma_muons"  );
236      Scatter2DPtr mult;
237      book(mult, 6,1,2);
238      for (size_t b = 0; b < temphisto.numPoints(); b++) {
239	const double x  = temphisto.point(b).x();
240	pair<double,double> ex = temphisto.point(b).xErrs();
241	pair<double,double> ex2 = ex;
242	if(ex2.first ==0.) ex2. first=0.0001;
243	if(ex2.second==0.) ex2.second=0.0001;
244	if (inRange(sqrtS()/MeV, x-ex2.first, x+ex2.second)) {
245	  mult   ->addPoint(x, rval, ex, rerr);
246	  hadrons->addPoint(x, sig_h, ex, make_pair(err_h,err_h));
247	  charm  ->addPoint(x, sig_c, ex, make_pair(err_c,err_c));
248	  muons  ->addPoint(x, sig_m, ex, make_pair(err_m,err_m));
249	}
250	else {
251	  mult   ->addPoint(x, 0., ex, make_pair(0.,.0));
252	  hadrons->addPoint(x, 0., ex, make_pair(0.,.0));
253	  charm  ->addPoint(x, 0., ex, make_pair(0.,.0));
254	  muons  ->addPoint(x, 0., ex, make_pair(0.,.0));
255	}
256      }
257      for(unsigned int ix=1;ix<6;++ix) {
258	unsigned int imax(0);
259	if     (ix<=3) imax = 4;
260	else           imax = 3;
261	for(unsigned int iy=1;iy<imax;++iy) {
262	  double sigma(0),error(0);
263	  if(ix==1) {
264	    if(iy==1) {
265	      sigma = _c_D0D0->val()/picobarn;
266	      error = _c_D0D0->err()/picobarn;
267	    }
268	    else if(iy==2) {
269	      sigma = _c_D0D0S->val()/picobarn;
270	      error = _c_D0D0S->err()/picobarn;
271	    }
272	    else if(iy==3) {
273	      sigma = _c_D0SD0S->val()/picobarn;
274	      error = _c_D0SD0S->err()/picobarn;
275	    }
276	  }
277	  else if(ix==2) {
278	    if(iy==1) {
279	      sigma = _c_DpDm->val()/picobarn;
280	      error = _c_DpDm->err()/picobarn;
281	    }
282	    else if(iy==2) {
283	      sigma = _c_DpDmS->val()/picobarn;
284	      error = _c_DpDmS->err()/picobarn;
285	    }
286	    else if(iy==3) {
287	      sigma = _c_DpSDmS->val()/picobarn;
288	      error = _c_DpSDmS->err()/picobarn;
289	    }
290	  }
291	  else if(ix==3) {
292	    if(iy==1) {
293	      sigma = _c_DsDs->val()/picobarn;
294	      error = _c_DsDs->err()/picobarn;
295	    }
296	    else if(iy==2) {
297	      sigma = _c_DsDsS->val()/picobarn;
298	      error = _c_DsDsS->err()/picobarn;
299	    }
300	    else if(iy==3) {
301	      sigma = _c_DsSDsS->val()/picobarn;
302	      error = _c_DsSDsS->err()/picobarn;
303	    }
304	  }
305	  else if(ix==4) {
306	    if(iy==1) {
307	      sigma = _c_DSDpi->val()/picobarn;
308	      error = _c_DSDpi->err()/picobarn;
309	    }
310	    else if(iy==2) {
311	      sigma = _c_DSDSpi->val()/picobarn;
312	      error = _c_DSDSpi->err()/picobarn;
313	    }
314	  }
315	  else if(ix==5) {
316	    if(iy==1) {
317	      sigma = _c_DD->val()/nanobarn;
318	      error = _c_DD->err()/nanobarn;
319	    }
320	    else if(iy==2) {
321	      sigma = _c_DDX->val()/nanobarn;
322	      error = _c_DDX->err()/nanobarn;
323	    }
324	  }
325	  sigma *= crossSection()/ sumOfWeights();
326	  error *= crossSection()/ sumOfWeights();
327	  Scatter2D temphisto(refData(ix, 1, iy));
328	  Scatter2DPtr mult;
329	  book(mult, ix,1,iy);
330	  for (size_t b = 0; b < temphisto.numPoints(); b++) {
331	    const double x  = temphisto.point(b).x();
332	    pair<double,double> ex = temphisto.point(b).xErrs();
333	    pair<double,double> ex2 = ex;
334	    if(ex2.first ==0.) ex2. first=0.0001;
335	    if(ex2.second==0.) ex2.second=0.0001;
336	    if (inRange(sqrtS()/MeV, x-ex2.first, x+ex2.second)) {
337	      mult   ->addPoint(x, sigma, ex, make_pair(error,error));
338	    }
339	    else {
340	      mult   ->addPoint(x, 0., ex, make_pair(0.,.0));
341	    }
342	  }
343	}
344      }
345    }
346
347    //@}
348
349
350    /// @name Histograms
351    //@{
352    CounterPtr _c_D0D0, _c_DpDm,_c_DsDs;
353    CounterPtr _c_D0D0S, _c_DpDmS,_c_DsDsS;
354    CounterPtr _c_D0SD0S, _c_DpSDmS,_c_DsSDsS;
355    CounterPtr _c_DD, _c_DDX;
356    CounterPtr _c_DSDpi, _c_DSDSpi;
357    CounterPtr _c_hadrons, _c_muons;
358    //@}
359
360
361  };
362
363
364  // The hook for the plugin system
365  RIVET_DECLARE_PLUGIN(CLEOC_2008_I777917);
366
367
368}