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