rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

BABAR_2015_I1377201

Azimuthal asymmetries in inclusive $\pi\pi$ $KK$ and $K\pi$ pairs at 10.58 GeV
Experiment: BABAR (PEP-II)
Inspire ID: 1377201
Status: VALIDATED
Authors:
  • Peter Richardson
References:
  • Phys.Rev. D92 (2015) no.11, 111101
Beams: e+ e-
Beam energies: (5.3, 5.3) GeV
Run details:
  • e+e- to hadrons

Measurement of azimuthal asymmetries in inclusive $\pi\pi$ $KK$ and $K\pi$ pair production at $\sqrt{s}=10.58$ GeV by the BABAR experiment

Source code: BABAR_2015_I1377201.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/Thrust.hh"
  4#include "Rivet/Projections/FinalState.hh"
  5#include "Rivet/Projections/Beam.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// @brief azimuthal asymmetries in pipi Kpi and KK
 11  class BABAR_2015_I1377201 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(BABAR_2015_I1377201);
 16
 17
 18    /// @name Analysis methods
 19    //@{
 20
 21    /// Book histograms and initialise projections before the run
 22    void init() {
 23      // projections
 24      const FinalState fs;
 25      declare(fs,"FS");
 26      declare(Thrust(fs),"Thrust");
 27      declare(Beam(), "Beams");
 28      // declare the histos for the distributions
 29      string type  [3] = {"KK","Kpi","pipi"};
 30      string charge[3] = {"Like","Opposite","All"};
 31      unsigned int nbin=20;
 32      for(unsigned int itype=0;itype<3;++itype) {
 33	for(unsigned int icharge=0;icharge<3;++icharge) {
 34	  for(unsigned int ibin=0;ibin<16;++ibin) {
 35	    std::ostringstream title1;
 36	    title1 << "/TMP/h_thrust" << type[itype] << "_" << charge[icharge] << "_" << ibin+1;
 37	    book(_h_thrust[itype][icharge][ibin],title1.str(),nbin,0.,M_PI);
 38	    std::ostringstream title2;
 39	    title2 << "/TMP/h_hadron" << type[itype] << "_" << charge[icharge] << "_" << ibin+1;
 40	    book(_h_hadron[itype][icharge][ibin],title2.str(),nbin,0.,M_PI);
 41	  }
 42	}
 43      }
 44    }
 45
 46    unsigned int iBin(double z) {
 47      if     (z<.2) return 0;
 48      else if(z<.3) return 1;
 49      else if(z<.5) return 2;
 50      else          return 3;
 51    }
 52    
 53    /// Perform the per-event analysis
 54    void analyze(const Event& event) {
 55      // get the axis, direction of incoming electron
 56      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 57      Vector3 axis1;
 58      if(beams.first.pid()>0)
 59	axis1 = beams.first .momentum().p3().unit();
 60      else
 61	axis1 = beams.second.momentum().p3().unit();
 62      // apply thrust cuts  T > 0.8  and | cos θ th | < 0.6
 63      Thrust thrust = apply<Thrust>(event,"Thrust");
 64      if(thrust.thrust()<=0.8) vetoEvent;
 65      if(cos(thrust.thrustAxis().polarAngle())>=0.6) vetoEvent;
 66      // construct x,y,z axes for thrust defn
 67      ThreeVector t_z = thrust.thrustAxis();
 68      ThreeVector t_x = (axis1-t_z.dot(axis1)*t_z).unit();
 69      ThreeVector t_y = t_z.cross(t_x);
 70      // loop over the particles
 71      Particles charged = apply<FinalState>(event,"FS").particles(Cuts::abspid==PID::PIPLUS || Cuts::abspid==PID::KPLUS);
 72      for(unsigned int ix=0;ix<charged.size();++ix) {
 73	// z and angle cut
 74	const double x1=2.*charged[ix].momentum().t()/sqrtS();
 75	if(x1<0.16||x1>.9) continue;
 76	if(abs(t_z.angle(charged[ix].momentum().p3()))>0.25*M_PI) continue;
 77	double dot1 = t_z.dot(charged[ix].p3());
 78	for(unsigned int iy=ix+1;iy<charged.size();++iy) {
 79	  const double x2=2.*charged[iy].momentum().t()/sqrtS();
 80	  // z and angle cut
 81	  if(x2<0.16||x2>.9) continue;
 82	  if(abs(t_z.angle(charged[ix].momentum().p3()))>0.25*M_PI) continue;
 83	  // different hemi
 84	  double dot2 = t_z.dot(charged[iy].p3());
 85	  if(dot1*dot2>0.) continue;
 86	  Particle p1=charged[ix], p2=charged[iy];
 87	  double z1(x1),z2(x2);
 88	  // randomly order the particles
 89	  if(rand()/static_cast<double>(RAND_MAX) < 0.5 ) {
 90	    swap(p1,p2);
 91	    swap(z1,z2);
 92	  }
 93	  // thrust def
 94	  double phi12 = atan2(p1.p3().dot(t_y),p1.p3().dot(t_x))+atan2(p2.p3().dot(t_y),p2.p3().dot(t_x));
 95	  if(phi12>M_PI)  phi12 -= 2*M_PI;
 96	  if(phi12<-M_PI) phi12 += 2*M_PI;
 97	  if(phi12<0.) phi12 = -phi12;
 98	  // hadron defn
 99	  ThreeVector h_z = p2.p3().unit();
100	  ThreeVector h_x = (axis1-h_z.dot(axis1)*h_z).unit();
101	  ThreeVector pt1 = p1.p3()-h_z.dot(p1.p3())*h_z;
102	  double phi0 = pt1.angle(h_x);
103	  if(phi0>M_PI)  phi0 -= 2*M_PI;
104	  if(phi0<-M_PI) phi0 += 2*M_PI;
105	  int ibin = 4*iBin(z1)+iBin(z2);
106	  // pi pi
107	  if(p1.abspid()==PID::PIPLUS && p2.abspid()==PID::PIPLUS) {
108	    if(p1.pid()==p2.pid()) {
109	      _h_thrust[2][0][ibin]->fill(phi12);
110	      _h_hadron[2][0][ibin]->fill(phi0);
111	    }
112	    else {
113	      _h_thrust[2][1][ibin]->fill(phi12);
114	      _h_hadron[2][1][ibin]->fill(phi0);
115	    }
116	    _h_thrust[2][2][ibin]->fill(phi12);
117	    _h_hadron[2][2][ibin]->fill(phi0);
118	  }
119	  // K K
120	  else if(p1.abspid()==PID::KPLUS && p2.abspid()==PID::KPLUS) {
121	    if(p1.pid()==p2.pid()) {
122	      _h_thrust[0][0][ibin]->fill(phi12);
123	      _h_hadron[0][0][ibin]->fill(phi0);
124	    }
125	    else {
126	      _h_thrust[0][1][ibin]->fill(phi12);
127	      _h_hadron[0][1][ibin]->fill(phi0);
128	    }
129	    _h_thrust[0][2][ibin]->fill(phi12);
130	    _h_hadron[0][2][ibin]->fill(phi0);
131	  }
132	  // K pi
133	  else {
134	    if(p1.pid()*p2.pid()>0) {
135	      _h_thrust[1][0][ibin]->fill(phi12);
136	      _h_hadron[1][0][ibin]->fill(phi0);
137	    }
138	    else {
139	      _h_thrust[1][1][ibin]->fill(phi12);
140	      _h_hadron[1][1][ibin]->fill(phi0);
141	    }
142	    _h_thrust[1][2][ibin]->fill(phi12);
143	    _h_hadron[1][2][ibin]->fill(phi0);
144	  }
145
146	  
147	}
148      }
149    }
150    
151    pair<double,double> calcAsymmetry(Scatter2DPtr hist,double fact=1.) {
152      double sum1(0.),sum2(0.);
153      for (auto point : hist->points() ) {
154	double Oi = point.y();
155	if(Oi==0. || std::isnan(Oi) ) continue;
156	double ai = 1.;
157	double bi = (sin(fact*point.xMax())-sin(fact*point.xMin()))/(point.xMax()-point.xMin())/fact;
158	double Ei = point.yErrAvg();
159	sum1 += sqr(bi/Ei);
160	sum2 += bi/sqr(Ei)*(Oi-ai);
161      }
162      if(sum1==0.) return make_pair(0.,0.);
163      return make_pair(sum2/sum1*1e4,sqrt(1./sum1)*1e4);
164    }
165
166    /// Normalise histograms etc., after the run
167    void finalize() {
168      for(unsigned int itype=0;itype<3;++itype) {
169	for(unsigned int icharge=0;icharge<3;++icharge) {
170	  for(unsigned int ibin=0;ibin<16;++ibin) {
171	    normalize(_h_thrust[itype][icharge][ibin]);
172	    normalize(_h_hadron[itype][icharge][ibin]);
173	  }
174	}
175      }
176      // construct ther ratios
177      // declare the histos for the distributions
178      string type  [3] = {"pipi","Kpi","KK"};
179      string charge[3] = {"Like","Opposite","All"};
180      for(unsigned int itype=0;itype<3;++itype) {
181	Scatter3DPtr h3_thrust_UL;
182	book(h3_thrust_UL,2*itype+1,1,2,0);
183	Scatter3DPtr h3_thrust_UC;
184	book(h3_thrust_UC,2*itype+1,1,3,0);
185	Scatter3DPtr h3_hadron_UL;
186	book(h3_hadron_UL,2*itype+2,1,2,0);
187	Scatter3DPtr h3_hadron_UC;
188	book(h3_hadron_UC,2*itype+2,1,3,0);
189	
190	unsigned int ihist=1;
191	Scatter2DPtr h2_thrust_UL;
192	book(h2_thrust_UL,7+2*itype,ihist,2);
193	Scatter2DPtr h2_thrust_UC;
194	book(h2_thrust_UC,7+2*itype,ihist,3);
195	Scatter2DPtr h2_hadron_UL;
196	book(h2_hadron_UL,8+2*itype,ihist,2);
197	Scatter2DPtr h2_hadron_UC;
198	book(h2_hadron_UC,8+2*itype,ihist,3);
199	
200	Scatter3D temphisto1(refData<Scatter3D>(2*itype+1, 1, 2));
201	Scatter3D temphisto2(refData<Scatter3D>(2*itype+2, 1, 2));
202	for(unsigned int ibin=0;ibin<16;++ibin) {
203	  const Point3D & p1 = temphisto1.points()[ibin];
204	  const Point3D & p2 = temphisto2.points()[ibin];
205	  if(ibin>0 && ibin%4==0) {
206	    ++ihist;
207	    book(h2_thrust_UL,7+2*itype,ihist,2);
208	    book(h2_thrust_UC,7+2*itype,ihist,3);
209	    book(h2_hadron_UL,8+2*itype,ihist,2);
210	    book(h2_hadron_UC,8+2*itype,ihist,3);
211	  }
212	  // thrust direction
213	  // opposite/like sign
214	  std::ostringstream title1;
215	  title1 << "/TMP/R_thrust_" << type[itype] << "_UL_" << ibin+1;
216	  Scatter2DPtr htemp;
217	  book(htemp,title1.str());
218	  divide(_h_thrust[itype][1][ibin],
219		 _h_thrust[itype][0][ibin],htemp);
220	  pair<double,double> asym = calcAsymmetry(htemp);
221	  h3_thrust_UL->addPoint(p1.x()    ,p1.y()    ,asym.first,
222				 p1.xErrs(),p1.yErrs(),make_pair(asym.second,asym.second) );
223	  h2_thrust_UL->addPoint(p1.y()    ,asym.first,p1.yErrs(),make_pair(asym.second,asym.second) );
224	  // opposite/all sign
225	  std::ostringstream title2;
226	  title2 << "/TMP/R_thrust_" << type[itype] << "_UC_" << ibin+1;
227	  book(htemp,title2.str());
228	  divide(_h_thrust[itype][1][ibin],
229		 _h_thrust[itype][2][ibin],htemp);
230	  asym = calcAsymmetry(htemp);
231	  h3_thrust_UC->addPoint(p1.x()    ,p1.y()    ,asym.first,
232				 p1.xErrs(),p1.yErrs(),make_pair(asym.second,asym.second) );
233	  h2_thrust_UC->addPoint(p1.y()    ,asym.first,p1.yErrs(),make_pair(asym.second,asym.second) );
234	  // hadron dirn
235	  // opposite/like sign
236	  std::ostringstream title3;
237	  title3 << "/TMP/R_hadron_" << type[itype] << "_UL_" << ibin+1;
238	  book(htemp,title3.str());
239	  divide(_h_hadron[itype][1][ibin],
240		 _h_hadron[itype][0][ibin],htemp);
241	  asym = calcAsymmetry(htemp,2.);
242	  h3_hadron_UL->addPoint(p2.x()    ,p2.y()    ,asym.first,
243				 p2.xErrs(),p2.yErrs(),make_pair(asym.second,asym.second) );
244	  h2_hadron_UL->addPoint(p2.y()    ,asym.first,p2.yErrs(),make_pair(asym.second,asym.second) );
245	  // opposite/all sign
246	  std::ostringstream title4;
247	  title4 << "/TMP/R_hadron_" << type[itype] << "_UC_" << ibin+1;
248	  book(htemp,title4.str());
249	  divide(_h_hadron[itype][1][ibin],
250		 _h_hadron[itype][2][ibin],htemp);
251	  asym = calcAsymmetry(htemp,2.);
252	  h3_hadron_UC->addPoint(p2.x()    ,p2.y()    ,asym.first,
253				 p2.xErrs(),p2.yErrs(),make_pair(asym.second,asym.second) );
254	  h2_hadron_UC->addPoint(p2.y()    ,asym.first,p2.yErrs(),make_pair(asym.second,asym.second) );
255	}
256      }
257    }
258
259    //@}
260
261
262    /// @name Histograms
263    //@{
264    Histo1DPtr _h_thrust[3][3][16],_h_hadron[3][3][16];
265    //@}
266
267
268  };
269
270
271  RIVET_DECLARE_PLUGIN(BABAR_2015_I1377201);
272
273}