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