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    pair<double,double> calcAsymmetry(Estimate1DPtr hist, double fact=1.) {
150      double sum1(0.),sum2(0.);
151      for (auto& b : hist->bins() ) {
152        double Oi = b.val();
153        if(Oi==0. || std::isnan(Oi) ) continue;
154        double ai = 1.;
155        double bi = (sin(fact*b.xMax())-sin(fact*b.xMin()))/(b.xWidth())/fact;
156        double Ei = b.errAvg();
157        sum1 += sqr(bi/Ei);
158        sum2 += bi/sqr(Ei)*(Oi-ai);
159      }
160      if(sum1==0.) return make_pair(0.,0.);
161      return make_pair(sum2/sum1*1e4,sqrt(1./sum1)*1e4);
162    }
163
164    /// Normalise histograms etc., after the run
165    void finalize() {
166      for(unsigned int itype=0;itype<3;++itype) {
167        for(unsigned int icharge=0;icharge<3;++icharge) {
168          for(unsigned int ibin=0;ibin<16;++ibin) {
169            normalize(_h_thrust[itype][icharge][ibin]);
170            normalize(_h_hadron[itype][icharge][ibin]);
171          }
172        }
173      }
174      // construct ther ratios
175      // declare the histos for the distributions
176      string type  [3] = {"pipi","Kpi","KK"};
177      string charge[3] = {"Like","Opposite","All"};
178      for(unsigned int itype=0;itype<3;++itype) {
179        Estimate2DPtr h3_thrust_UL;
180        book(h3_thrust_UL,2*itype+1,1,2);
181        Estimate2DPtr h3_thrust_UC;
182        book(h3_thrust_UC,2*itype+1,1,3);
183        Estimate2DPtr h3_hadron_UL;
184        book(h3_hadron_UL,2*itype+2,1,2);
185        Estimate2DPtr h3_hadron_UC;
186        book(h3_hadron_UC,2*itype+2,1,3);
187
188        unsigned int ihist=1;
189        Estimate1DPtr h2_thrust_UL;
190        book(h2_thrust_UL,7+2*itype,ihist,2);
191        Estimate1DPtr h2_thrust_UC;
192        book(h2_thrust_UC,7+2*itype,ihist,3);
193        Estimate1DPtr h2_hadron_UL;
194        book(h2_hadron_UL,8+2*itype,ihist,2);
195        Estimate1DPtr h2_hadron_UC;
196        book(h2_hadron_UC,8+2*itype,ihist,3);
197
198        Estimate2D temphisto1(refData<Estimate2D>(2*itype+1, 1, 2));
199        Estimate2D temphisto2(refData<Estimate2D>(2*itype+2, 1, 2));
200        for (unsigned int ibin=0;ibin<16;++ibin) {
201          //const auto& b1 = temphisto1.bin(ibin);
202          //const auto& b2 = temphisto2.bin(ibin);
203          if (ibin>0 && ibin%4==0) {
204            ++ihist;
205            book(h2_thrust_UL,7+2*itype,ihist,2);
206            book(h2_thrust_UC,7+2*itype,ihist,3);
207            book(h2_hadron_UL,8+2*itype,ihist,2);
208            book(h2_hadron_UC,8+2*itype,ihist,3);
209          }
210          // thrust direction
211          // opposite/like sign
212          std::ostringstream title1;
213          title1 << "/TMP/R_thrust_" << type[itype] << "_UL_" << ibin+1;
214          Estimate1DPtr htemp;
215          book(htemp,title1.str());
216          divide(_h_thrust[itype][1][ibin],
217          _h_thrust[itype][0][ibin],htemp);
218          pair<double,double> asym = calcAsymmetry(htemp);
219          h3_thrust_UL->bin(ibin+1).set(asym.first, asym.second);
220          h2_thrust_UL->bin(ibin+1).set(asym.first, asym.second);
221          // opposite/all sign
222          std::ostringstream title2;
223          title2 << "/TMP/R_thrust_" << type[itype] << "_UC_" << ibin+1;
224          book(htemp,title2.str());
225          divide(_h_thrust[itype][1][ibin],
226          _h_thrust[itype][2][ibin],htemp);
227          asym = calcAsymmetry(htemp);
228          h3_thrust_UC->bin(ibin+1).set(asym.first, asym.second);
229          h2_thrust_UC->bin(ibin+1).set(asym.first, asym.second);
230          // hadron dirn
231          // opposite/like sign
232          std::ostringstream title3;
233          title3 << "/TMP/R_hadron_" << type[itype] << "_UL_" << ibin+1;
234          book(htemp,title3.str());
235          divide(_h_hadron[itype][1][ibin],
236          _h_hadron[itype][0][ibin],htemp);
237          asym = calcAsymmetry(htemp,2.);
238          h3_hadron_UL->bin(ibin+1).set(asym.first, asym.second);
239          h2_hadron_UL->bin(ibin+1).set(asym.first, asym.second);
240          // opposite/all sign
241          std::ostringstream title4;
242          title4 << "/TMP/R_hadron_" << type[itype] << "_UC_" << ibin+1;
243          book(htemp,title4.str());
244          divide(_h_hadron[itype][1][ibin],
245          _h_hadron[itype][2][ibin],htemp);
246          asym = calcAsymmetry(htemp,2.);
247          h3_hadron_UC->bin(ibin+1).set(asym.first, asym.second);
248          h2_hadron_UC->bin(ibin+1).set(asym.first, asym.second);
249        }
250      }
251    }
252
253    /// @}
254
255
256    /// @name Histograms
257    /// @{
258    Histo1DPtr _h_thrust[3][3][16],_h_hadron[3][3][16];
259    /// @}
260
261
262  };
263
264
265  RIVET_DECLARE_PLUGIN(BABAR_2015_I1377201);
266
267}