rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

BESIII_2023_I2636760

Decay assymetry in $\Sigma^+\to p\gamma$
Experiment: BESIII ()
Inspire ID: 2636760
Status: VALIDATED NOHEPDATA
Authors:
  • Your Name
References:
  • 2302.13568 [hep-ex]
Beams: e- e+
Beam energies: (1.6, 1.6) GeV
Run details:
  • e+e- > J/psi

Analysis of the angular distribution of the baryons, and decay products, produced in $e^+e^-\to J/\psi \to \Sigma^+\bar\Sigma^-$. Gives information about the decay and is useful for testing correlations in hadron decays.

Source code: BESIII_2023_I2636760.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 Sigma + -> p gamma
 11  class BESIII_2023_I2636760 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2023_I2636760);
 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(UnstableParticles(), "UFS");
 27      declare(FinalState(), "FS");
 28
 29      // Book histograms
 30      for (unsigned int ix=0; ix<2; ++ix) {
 31        book(_n[ix],"TMP/n_" + toString(ix+1));
 32        book(_t[ix],"TMP/t_" + toString(ix+1));
 33      }
 34      book(_n[2],"TMP/n_3");
 35      book(_t[2],"TMP/t_3");
 36    }
 37
 38
 39    void findChildren(const Particle& p,map<long,int>& nRes, int& ncount) {
 40      for (const Particle& child : p.children()) {
 41        if (child.children().empty()) {
 42          nRes[child.pid()]-=1;
 43          --ncount;
 44        }
 45        else {
 46          findChildren(child,nRes,ncount);
 47        }
 48      }
 49    }
 50
 51    /// Perform the per-event analysis
 52    void analyze(const Event& event) {
 53      // get the axis, direction of incoming electron
 54      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 55      Vector3 axis;
 56      if (beams.first.pid()>0) {
 57        axis = beams.first .mom().p3().unit();
 58      }
 59      else {
 60	      axis = beams.second.mom().p3().unit();
 61      }
 62      // types of final state particles
 63      const FinalState& fs = apply<FinalState>(event, "FS");
 64      map<long,int> nCount;
 65      int ntotal(0);
 66      for (const Particle& p :  fs.particles()) {
 67        nCount[p.pid()] += 1;
 68        ++ntotal;
 69      }
 70      // loop over Sigma+ baryons
 71      const UnstableParticles & ufs = apply<UnstableParticles>(event, "UFS");
 72      Particle Sigma, SigBar;
 73      bool matched(false);
 74      for (const Particle& p :  ufs.particles(Cuts::abspid==3222)) {
 75       	if (p.children().empty()) continue;
 76       	map<long,int> nRes=nCount;
 77       	int ncount = ntotal;
 78       	findChildren(p, nRes, ncount);
 79       	matched=false;
 80       	// check for antiparticle
 81      	for (const Particle& p2 :  ufs.particles(Cuts::pid==-p.pid())) {
 82      	  if(p2.children().empty()) continue;
 83      	  map<long,int> nRes2=nRes;
 84      	  int ncount2 = ncount;
 85      	  findChildren(p2, nRes2, ncount2);
 86      	  if (ncount2==0) {
 87      	    matched = true;
 88      	    for (const auto& val : nRes2) {
 89      	      if (val.second!=0) {
 90                matched = false;
 91                break;
 92      	      }
 93      	    }
 94      	    // fond baryon and antibaryon
 95      	    if (matched) {
 96              if (p.pid()>0) {
 97                Sigma = p;
 98                SigBar = p2;
 99              }
100              else {
101                Sigma = p2;
102                SigBar = p;
103              }
104       	      break;
105       	    }
106       	  }
107       	}
108      	if (matched) break;
109      }
110      if (!matched) vetoEvent;
111      if (Sigma .children().size()!=2) vetoEvent;
112      if (SigBar.children().size()!=2) vetoEvent;
113      // find proton
114      bool radiative[2]={false,false};
115      // identify Sigma decay
116      Particle baryon1;
117      if (Sigma.children()[0].pid()==PID::PROTON && Sigma.children()[1].pid()==PID::PI0) {
118        radiative[0]=false;
119        baryon1 = Sigma.children()[0];
120      }
121      else if (Sigma.children()[1].pid()==PID::PROTON && Sigma.children()[0].pid()==PID::PI0) {
122        radiative[0]=false;
123        baryon1 = Sigma.children()[1];
124      }
125      else if (Sigma.children()[0].pid()==PID::PROTON && Sigma.children()[1].pid()==PID::PHOTON) {
126        radiative[0]=true;
127        baryon1 = Sigma.children()[0];
128      }
129      else if (Sigma.children()[1].pid()==PID::PROTON && Sigma.children()[0].pid()==PID::PHOTON ) {
130        radiative[0]=true;
131        baryon1 = Sigma.children()[1];
132      }
133      else {
134        vetoEvent;
135      }
136      // antisigma decay
137      Particle baryon2;
138      if (SigBar.children()[0].pid()==PID::ANTIPROTON && SigBar.children()[1].pid()==PID::PI0 ) {
139        radiative[1]=false;
140        baryon2 = SigBar.children()[0];
141      }
142      else if (SigBar.children()[1].pid()==PID::ANTIPROTON && SigBar.children()[0].pid()==PID::PI0 ) {
143        radiative[1]=false;
144        baryon2 = SigBar.children()[1];
145      }
146      else if (SigBar.children()[0].pid()==PID::ANTIPROTON && SigBar.children()[1].pid()==PID::PHOTON ) {
147        radiative[1]=true;
148        baryon2 = SigBar.children()[0];
149      }
150      else if (SigBar.children()[1].pid()==PID::ANTIPROTON && SigBar.children()[0].pid()==PID::PHOTON ) {
151        radiative[1]=true;
152        baryon2 = SigBar.children()[1];
153      }
154      else {
155        vetoEvent;
156      }
157      if (radiative[0] == radiative[1]) vetoEvent;
158      // boost to the Sigma rest frame
159      LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(Sigma.mom().betaVec());
160      Vector3 e1z = Sigma.mom().p3().unit();
161      Vector3 e1y = e1z.cross(axis).unit();
162      Vector3 e1x = e1y.cross(e1z).unit();
163      Vector3 axis1 = boost1.transform(baryon1.mom()).p3().unit();
164      double n1x(e1x.dot(axis1)), n1z(e1z.dot(axis1));
165      // boost to the Sigma bar
166      LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(SigBar.mom().betaVec());
167      Vector3 axis2 = boost2.transform(baryon2.mom()).p3().unit();
168      double n2x(e1x.dot(axis2)),n2z(e1z.dot(axis2));
169      double cosL = axis.dot(Sigma.mom().p3().unit());
170      double sinL = sqrt(1.-sqr(cosL));
171      double T1 = sqr(sinL)*n1x*n2x+sqr(cosL)*n1z*n2z;
172      // sigma+ -> p gamma
173      if (radiative[0]) {
174        _n[0]->fill();
175        _n[2]->fill();
176        _t[0]->fill(T1);
177        _t[2]->fill(T1);
178      }
179      // sigmabar- -> pbar gamma
180      else {
181        _n[1]->fill();
182        _n[2]->fill();
183        _t[1]->fill(T1);
184        _t[2]->fill(T1);
185      }
186    }
187
188    /// Normalise histograms etc., after the run
189    void finalize() {
190      // values of constants
191      const double aPsi  =-0.508;
192      const double aPlus =-0.998;
193      const double factor = 45.*(3. +aPsi)/(11. + 5.*aPsi)/aPlus;
194      // alpha from the moments
195      for (unsigned int ix=0;ix<3;++ix) {
196        double value = _t[ix]->val()/_n[ix]->val();
197        double error = _t[ix]->err()/_n[ix]->val();
198        value *= factor;
199        error *= abs(factor);
200        if (ix==1) value *=-1.;
201        Estimate0DPtr alpha;
202        book(alpha, 1, 1, 1+ix);
203        alpha->set(value, error);
204      }
205
206    }
207
208    /// @}
209
210
211    /// @name Histograms
212    /// @{
213    CounterPtr _n[3];
214    Histo1DPtr _h_ctheta[3];
215    CounterPtr _t[3];
216    /// @}
217
218
219  };
220
221
222  RIVET_DECLARE_PLUGIN(BESIII_2023_I2636760);
223
224}