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