rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

BESIII_2016_I1384778

Azimuthal asymmetries in inclusive charged pion-pair production at $\sqrt{s}=3.65$ GeV
Experiment: BESIII (BEPC)
Inspire ID: 1384778
Status: VALIDATED
Authors:
  • Peter Richardson
References:
  • Phys.Rev.Lett. 116 (2016) no.4, 042001
Beams: e+ e-
Beam energies: (1.8, 1.8) GeV
Run details:
  • e+e- to hadrons

Measurement of azimuthal asymmetries in inclusive charged pion-pair production at $\sqrt{s}=3.65$ GeV by the BESII experiment

Source code: BESIII_2016_I1384778.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/Beam.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief Collins assymmetry
 10  class BESIII_2016_I1384778 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2016_I1384778);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22      declare(Beam(), "Beams");
 23      declare(FinalState(Cuts::abspid==PID::PIPLUS), "FS");
 24      // book the histograms
 25      _h_L.resize(6);
 26      _h_U.resize(6);
 27      _h_C.resize(6);
 28      for(size_t ix=0; _h_L.size(); ++ix) {
 29        const string pre = "/TMP/h_z1z2_"+to_string(ix+1);
 30        book(_h_L[ix], pre+"_L", 20, 0., M_PI);
 31        book(_h_U[ix], pre+"_U", 20, 0., M_PI);
 32        book(_h_C[ix], pre+"_C", 20, 0., M_PI);
 33      }
 34      book(_h_pT_L, {0., 0.2, 0.3, 0.45, 0.8, 1.4});
 35      book(_h_pT_U, {0., 0.2, 0.3, 0.45, 0.8, 1.4});
 36      book(_h_pT_C, {0., 0.2, 0.3, 0.45, 0.8, 1.4});
 37      for (size_t ix=1; ix < _h_pT_L->numBins()+1; ++ix) {
 38        const string pre = "/TMP/h_pT_"+to_string(ix);
 39        book(_h_pT_L->bin(ix), pre+"_L", 20, 0.0, M_PI);
 40        book(_h_pT_U->bin(ix), pre+"_U", 20, 0.0, M_PI);
 41        book(_h_pT_C->bin(ix), pre+"_C", 20, 0.0, M_PI);
 42      }
 43    }
 44
 45
 46    /// Perform the per-event analysis
 47    void analyze(const Event& event) {
 48      // get the axis, direction of incoming electron
 49      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 50      Vector3 axis;
 51      if (beams.first.pid()>0) {
 52        axis = beams.first .momentum().p3().unit();
 53      }
 54      else {
 55        axis = beams.second.momentum().p3().unit();
 56      }
 57      // loop over pions pair, using index to avoid double counting
 58      Particles pions = apply<FinalState>(event, "FS").particles();
 59      for (unsigned int i1=0; i1<pions.size(); ++i1) {
 60        const double x1=2.*pions[i1].momentum().t()/sqrtS();
 61        // cut on z1
 62        if (x1<0.2||x1>0.9) continue;
 63        // cos theta cut
 64        if (abs(cos(pions[i1].momentum().p3().polarAngle()))>0.93) continue;
 65        for(unsigned int i2=i1+1;i2<pions.size();++i2) {
 66          // cut on z2
 67          const double x2=2.*pions[i2].momentum().t()/sqrtS();
 68          if (x2<0.2||x2>0.9) continue;
 69          // cos theta cut
 70          if (abs(cos(pions[i2].momentum().p3().polarAngle()))>0.93) continue;
 71          // cut on opening angle (>120 degrees)
 72          if (pions[i1].momentum().p3().angle(pions[i2].momentum().p3())>2.*M_PI/3.) {
 73            continue;
 74          }
 75          Particle p1=pions[i1], p2=pions[i2];
 76          double z1(x1),z2(x2);
 77          // randomly order the particles
 78          if (rand()/static_cast<double>(RAND_MAX) < 0.5) {
 79            swap(p1,p2);
 80            swap(z1,z2);
 81          }
 82          // particle 2 defines the z axis
 83          Vector3 ez = p2.momentum().p3().unit();
 84          // beam and 2 define the plane (y is normal to plane)
 85          Vector3 ey = ez.cross(axis).unit();
 86          // x by cross product
 87          Vector3 ex = ey.cross(ez).unit();
 88          // phi
 89          double phi = ex.angle(p1.momentum().p3());
 90          // hists vs z1,z2
 91          unsigned int ibin=0;
 92          if (z1<=.3&&z2<=.3) {
 93            ibin=0;
 94          }
 95          else if (z1>0.5&&z2>0.5) {
 96            ibin=5;
 97          }
 98          else if (min(z1,z2)<=0.3) {
 99            if (max(z1,z2)>0.5) {
100              ibin=2;
101            }
102            else {
103              ibin=1;
104            }
105          }
106          else {
107            if (max(z1,z2)>0.5) {
108              ibin=4;
109            }
110            else {
111              ibin=3;
112            }
113          }
114          _h_C[ibin]->fill(phi);
115          if (p1.pid()==p2.pid()) {
116            _h_L[ibin]->fill(phi);
117          }
118          else {
119            _h_U[ibin]->fill(phi);
120          }
121          // hists vs pT
122          double pPar2 = sqr(ez.dot(p1.momentum().p3()));
123          double pPerp = sqrt(p1.momentum().p3().mod2()-pPar2);
124          _h_pT_C->fill(pPerp,phi);
125          if(p1.pid()==p2.pid()) {
126            _h_pT_L->fill(pPerp,phi);
127          }
128          else {
129            _h_pT_U->fill(pPerp,phi);
130          }
131        }
132      }
133    }
134
135    pair<double,double> calcAsymmetry(Estimate1DPtr hist) {
136      double sum1(0.),sum2(0.);
137      for (const auto& b : hist->bins()) {
138        double Oi = b.val();
139        if(Oi==0. || std::isnan(Oi) ) continue;
140        double ai = 1.;
141        double bi = 0.5*(sin(2.*b.xMax())-sin(2.*b.xMin()))/b.xWidth();
142        double Ei = b.errAvg();
143        sum1 += sqr(bi/Ei);
144        sum2 += bi/sqr(Ei)*(Oi-ai);
145      }
146      return make_pair(sum2/sum1, sqrt(1./sum1));
147    }
148
149
150    /// Normalise histograms etc., after the run
151    void finalize() {
152      // ratios
153      Estimate1DPtr _h_z_UL,_h_z_UC;
154      book(_h_z_UL,1,1,5);
155      book(_h_z_UC,1,1,6);
156      for (size_t ix=0; ix < _h_L.size(); ++ix) {
157        normalize(_h_L[ix]);
158        normalize(_h_U[ix]);
159        normalize(_h_C[ix]);
160        const string pre = "/TMP/R_z1z2_"+to_string(ix+1);
161        Estimate1DPtr R1;
162        book(R1, pre+"_UL");
163        divide(_h_U[ix], _h_L[ix], R1);
164        Estimate1DPtr R2;
165        book(R2, pre+"_UC");
166        divide(_h_U[ix], _h_C[ix], R2);
167        pair<double,double> asym1 = calcAsymmetry(R1);
168        _h_z_UL->bin(ix+1).set(asym1.first, asym1.second);
169        pair<double,double> asym2 = calcAsymmetry(R2);
170        _h_z_UC->bin(ix+1).set(asym2.first, asym2.second);
171      }
172      Estimate1DPtr _h_pT_UL,_h_pT_UC;
173      book(_h_pT_UL,2,1,4);
174      book(_h_pT_UC,2,1,5);
175      for (size_t ix=1; ix < _h_pT_L->numBins()+1; ++ix) {
176        normalize(_h_pT_L->bin(ix));
177        normalize(_h_pT_U->bin(ix));
178        normalize(_h_pT_C->bin(ix));
179        const string pre = "/TMP/R_pT_"+to_string(ix);
180        Estimate1DPtr R1;
181        book(R1, pre+"_UL");
182        divide(_h_pT_U->bin(ix), _h_pT_L->bin(ix), R1);
183        Estimate1DPtr R2;
184        book(R2, pre+"_UC");
185        divide(_h_pT_U->bin(ix), _h_pT_C->bin(ix), R2);
186        pair<double,double> asym1 = calcAsymmetry(R1);
187        _h_pT_UL->bin(ix).set(asym1.first, asym1.second);
188        pair<double,double> asym2 = calcAsymmetry(R2);
189        _h_pT_UC->bin(ix).set(asym2.first, asym2.second);
190      }
191    }
192    /// @}
193
194
195    /// @name Histograms
196    /// @{
197    vector<Histo1DPtr> _h_L,_h_U,_h_C;
198    Histo1DGroupPtr _h_pT_L,_h_pT_U,_h_pT_C;
199    /// @}
200
201
202  };
203
204
205  RIVET_DECLARE_PLUGIN(BESIII_2016_I1384778);
206
207}