rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

BELLE_2020_I1777678

Single and dihadron scaled momenta spectra at 10.58 GeV
Experiment: BELLE (KEKB)
Inspire ID: 1777678
Status: VALIDATED
Authors:
  • Peter Richardson
References:
  • Phys.Rev.D 101 (2020) 9, 092004
Beams: e+ e-
Beam energies: (5.3, 5.3) GeV
Run details:
  • e+ e- to hadrons

Measurement of single and di-hadron spectra by the BELLE collaboration at 10.58 GeV

Source code: BELLE_2020_I1777678.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/ChargedFinalState.hh"
  4#include "Rivet/Projections/FinalState.hh"
  5#include "Rivet/Projections/Thrust.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// @brief Single and di-hadron spectra
 11  class BELLE_2020_I1777678 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2020_I1777678);
 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(ChargedFinalState(Cuts::abspid==211 or Cuts::abspid==321 or Cuts::abspid==2212),"CFS");
 26      // projections
 27      FinalState fs;
 28      declare(fs,"FS");
 29      declare(Thrust(fs),"Thrust");
 30      // single particle hists
 31      vector<int> pdg={211,321,2212};
 32      for(unsigned int ix=0;ix<3;++ix) {
 33	book(_s_all   [pdg[ix]],1,ix+1,1);
 34	book(_s_strong[pdg[ix]],1,ix+1,2);
 35      }
 36      // dihadron histograms
 37      double bins[17]={0.20,0.25,0.30,0.35,0.40,0.45,0.50,0.55,0.60,
 38      		       0.65,0.70,0.75,0.80,0.85,0.90,0.95,1.00};
 39      unsigned int i0=1;
 40      for(unsigned int defn=0;defn<3;++defn) {
 41	for(unsigned int hemi=0;hemi<3;++hemi) {
 42	  for(unsigned int ip=0;ip<6;++ip) {
 43	    i0+=1;
 44	    unsigned int ymax=16;
 45	    if(i0==7 || i0==19)      ymax=15;
 46	    else if(i0>=8 && i0<=12) ymax=14;
 47	    else if(i0==13)          ymax=13;
 48	    else if(i0==26)          ymax=10;
 49	    else if(i0==27||i0==30)  ymax= 9;
 50	    else if(i0==28)          ymax= 7;
 51	    else if(i0==29)          ymax= 8;
 52	    else if(i0==31||i0==44)  ymax= 6;
 53	    else if(i0==45||i0==48)  ymax= 5;
 54	    else if(i0==46||i0==47)  ymax= 4;
 55	    else if(i0==49        )  ymax= 3;
 56	    for(unsigned int iy=0;iy<ymax;++iy) {
 57	      Histo1DPtr temp;
 58	      book(temp,i0,1,iy+1);
 59	      _d_all   [ip][defn][hemi].add(bins[iy],bins[iy+1],temp);
 60	      book(temp,i0,2,iy+1);
 61	      _d_strong[ip][defn][hemi].add(bins[iy],bins[iy+1],temp);
 62	    }
 63	  }
 64	}
 65      }
 66    }
 67
 68    bool isWeak(const Particle & p) {
 69      bool weak = false;
 70      if(p.parents().empty()) return weak;
 71      Particle parent = p.parents()[0];
 72      while (!parent.parents().empty()) {
 73	if(parent.abspid()==411  || parent.abspid()==421  || parent.abspid()==431  ||
 74	   parent.abspid()==4122 || parent.abspid()==4232 || parent.abspid()==4132 ||
 75	   parent.abspid()==4332) {
 76	  weak=true;
 77	  break;
 78	}
 79	parent = parent.parents()[0];
 80      }
 81      return weak;
 82    }
 83    
 84    void fillHistos(int ip,bool strong,bool same,bool opp,
 85		    const Particle & p1, const Particle & p2) {
 86      for(unsigned int def=0;def<3;++def) {
 87	double z1 = 0., z2 = 0.;
 88	if(def==0) {
 89	  z1 = 2.*p1.momentum().t()/sqrtS();
 90	  z2 = 2.*p2.momentum().t()/sqrtS();
 91	}
 92	else if(def==1) {
 93	  z1 = 2.*p1.momentum().t()/sqrtS();
 94	  z2 = (p1.momentum()*p2.momentum())/p1.momentum().t()/sqrtS();
 95	}
 96	else if(def==2) {
 97	  double p1p2 = p1.momentum()*p2.momentum();
 98	  double p1q = p1.momentum().t()*sqrtS();
 99	  double p2q = p2.momentum().t()*sqrtS();
100	  z1 = (p1p2-p1.mass2()*p2.mass2()/p1p2)/(p2q-p2.mass2()*p1q/p1p2);
101	  z2 = (p1p2-p1.mass2()*p2.mass2()/p1p2)/(p1q-p1.mass2()*p2q/p1p2);
102	}
103	_d_all[ip][def][0].fill(z1,z2,0.5);
104	if(strong) _d_strong[ip][def][0].fill(z1,z2,0.5);
105	if(same) {
106	  _d_all[ip][def][1].fill(z1,z2,0.5);
107	  if(strong) _d_strong[ip][def][1].fill(z1,z2,0.5);
108	}
109	if(opp) {
110	  _d_all[ip][def][2].fill(z1,z2,0.5);
111	  if(strong) _d_strong[ip][def][2].fill(z1,z2,0.5);
112	}
113      }
114    }
115
116    /// Perform the per-event analysis
117    void analyze(const Event& event) {
118      // apply projection
119      const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");
120      // fill single particle histos
121      for (const Particle& p : cfs.particles()) {
122	const double z = 2.*p.momentum().t()/sqrtS();
123	_s_all   [p.abspid()]->fill(z);
124	if(!isWeak(p)) _s_strong[p.abspid()]->fill(z);
125      }
126      // get thrust
127      const Thrust thrust = apply<Thrust>(event,"Thrust");
128      ThreeVector axis = thrust.thrustAxis();
129      Particles piK = cfs.particles(Cuts::abspid==PID::KPLUS or
130				    Cuts::abspid==PID::PIPLUS);
131      for(unsigned int ix=0;ix<piK.size();++ix) {
132	double dot1 = axis.dot(piK[ix].momentum().p3());
133	bool weak1 = isWeak(piK[ix]);
134	for(unsigned int iy=0;iy<piK.size();++iy) {
135	  if(ix==iy) continue;
136	  double dot2 = axis.dot(piK[iy].momentum().p3());
137	  bool weak2 = isWeak(piK[iy]);
138	  bool strong = !weak1 && !weak2;
139	  bool same = thrust.thrust()>0.8 && dot1*dot2>0.;
140	  bool opp  = thrust.thrust()>0.8 && dot1*dot2<0.;
141	  unsigned int ip=0;
142	  if(piK[ix].pid()==PID::PIPLUS) {
143	    if(piK[iy].pid()==PID::PIPLUS)       ip=1;
144	    else if(piK[iy].pid()==PID::PIMINUS) ip=0;
145	    else if(piK[iy].pid()==PID::KPLUS  ) ip=3;
146	    else if(piK[iy].pid()==PID::KMINUS)  ip=2;
147	  }
148	  else if(piK[ix].pid()==PID::PIMINUS) {
149	    if(piK[iy].pid()==PID::PIPLUS)       ip=0;
150	    else if(piK[iy].pid()==PID::PIMINUS) ip=1;
151	    else if(piK[iy].pid()==PID::KPLUS  ) ip=2;
152	    else if(piK[iy].pid()==PID::KMINUS ) ip=3;
153	  }
154	  else if(piK[ix].pid()==PID::KPLUS) {
155	    if(piK[iy].pid()==PID::PIPLUS)       ip=3;
156	    else if(piK[iy].pid()==PID::PIMINUS) ip=2;
157	    else if(piK[iy].pid()==PID::KPLUS)   ip=5;
158	    else if(piK[iy].pid()==PID::KMINUS)  ip=4;
159	  }
160	  else if(piK[ix].pid()==PID::KMINUS) {
161	    if(piK[iy].pid()==PID::PIPLUS)       ip=2;
162	    else if(piK[iy].pid()==PID::PIMINUS) ip=3;
163	    else if(piK[iy].pid()==PID::KPLUS)   ip=4;
164	    else if(piK[iy].pid()==PID::KMINUS)  ip=5;
165	  }
166	  fillHistos(ip,strong,same,opp,piK[ix],piK[iy]);
167	}
168      }
169    }
170
171    /// Normalise histograms etc., after the run
172    void finalize() {
173      for (const auto& kv : _s_all)
174	scale(kv.second,crossSection()/femtobarn/sumOfWeights());
175      for (const auto& kv : _s_strong)
176	scale(kv.second,crossSection()/femtobarn/sumOfWeights());
177      for(unsigned int ix=0;ix<6;++ix) {
178	for(unsigned int iy=0;iy<3;++iy) {
179	  for(unsigned int iz=0;iz<3;++iz) {
180	    _d_all   [ix][iy][iz].scale(crossSection()/femtobarn/sumOfWeights(),this);
181	    _d_strong[ix][iy][iz].scale(crossSection()/femtobarn/sumOfWeights(),this);
182	  }
183	}
184      }
185    }
186
187    ///@}
188
189
190    /// @name Histograms
191    ///@{
192    map<int,Histo1DPtr> _s_all,_s_strong;
193    BinnedHistogram _d_all[6][3][3],_d_strong[6][3][3];
194    ///@}
195
196
197  };
198
199
200  RIVET_DECLARE_PLUGIN(BELLE_2020_I1777678);
201
202}