1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/ChargedFinalState.hh"
  4#include "Rivet/Projections/FinalState.hh"
  5#include "Rivet/Projections/Thrust.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// @brief Single and di-hadron spectra
 11  class BELLE_2020_I1777678 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2020_I1777678);
 16
 17
 18    /// @name Analysis methods
 19    /// @{
 20
 21    /// Book histograms and initialise projections before the run
 22    void init() {
 23
 24      // Initialise and register projections
 25      declare(ChargedFinalState(Cuts::abspid==211 or Cuts::abspid==321 or Cuts::abspid==2212),"CFS");
 26      // projections
 27      FinalState fs;
 28      declare(fs,"FS");
 29      declare(Thrust(fs),"Thrust");
 30      // single particle hists
 31      vector<int> pdg={211,321,2212};
 32      for (size_t ix=0; ix<3; ++ix) {
 33        book(_s_all   [pdg[ix]], 1, ix+1, 1);
 34        book(_s_strong[pdg[ix]], 1, ix+1, 2);
 35      }
 36      // dihadron histograms
 37      const vector<double> edges{0.20,0.25,0.30,0.35,0.40,0.45,0.50,0.55,0.60,
 38                                 0.65,0.70,0.75,0.80,0.85,0.90,0.95,1.00};
 39      size_t i0=1;
 40      for (size_t defn=0; defn<3; ++defn) {
 41        for (size_t hemi=0; hemi<3; ++hemi) {
 42          for (size_t ip=0; ip<6; ++ip) {
 43            ++i0;
 44            size_t ymax=16;
 45            if(i0==7 || i0==19)      ymax=15;
 46            else if(i0>=8 && i0<=12) ymax=14;
 47            else if(i0==13)          ymax=13;
 48            else if(i0==26)          ymax=10;
 49            else if(i0==27||i0==30)  ymax= 9;
 50            else if(i0==28)          ymax= 7;
 51            else if(i0==29)          ymax= 8;
 52            else if(i0==31||i0==44)  ymax= 6;
 53            else if(i0==45||i0==48)  ymax= 5;
 54            else if(i0==46||i0==47)  ymax= 4;
 55            else if(i0==49        )  ymax= 3;
 56            book(_d_all[ip][defn][hemi], edges);
 57            book(_d_strong[ip][defn][hemi], edges);
 58            for (size_t iy=1; iy < _d_all[ip][defn][hemi]->numBins()+1; ++iy) {
 59              if (iy <= ymax) {
 60                book(_d_all[ip][defn][hemi]->bin(iy), i0, 1, iy);
 61                book(_d_strong[ip][defn][hemi]->bin(iy), i0, 2, iy);
 62              }
 63              else {
 64                _d_all[ip][defn][hemi]->maskBin(iy);
 65                _d_strong[ip][defn][hemi]->maskBin(iy);
 66              }
 67
 68            }
 69          }
 70        }
 71      }
 72    }
 73
 74    bool isWeak(const Particle & p) {
 75      bool weak = false;
 76      if(p.parents().empty()) return weak;
 77      Particle parent = p.parents()[0];
 78      while (!parent.parents().empty()) {
 79        if(parent.abspid()==411  || parent.abspid()==421  || parent.abspid()==431  ||
 80           parent.abspid()==4122 || parent.abspid()==4232 || parent.abspid()==4132 ||
 81           parent.abspid()==4332) {
 82          weak=true;
 83          break;
 84        }
 85        parent = parent.parents()[0];
 86      }
 87      return weak;
 88    }
 89
 90    void fillHistos(int ip,bool strong,bool same,bool opp,
 91                    const Particle & p1, const Particle & p2) {
 92      for (size_t def=0; def<3; ++def) {
 93        double z1 = 0., z2 = 0.;
 94        if(def==0) {
 95          z1 = 2.*p1.momentum().t()/sqrtS();
 96          z2 = 2.*p2.momentum().t()/sqrtS();
 97        }
 98        else if(def==1) {
 99          z1 = 2.*p1.momentum().t()/sqrtS();
100          z2 = (p1.momentum()*p2.momentum())/p1.momentum().t()/sqrtS();
101        }
102        else if(def==2) {
103          double p1p2 = p1.momentum()*p2.momentum();
104          double p1q = p1.momentum().t()*sqrtS();
105          double p2q = p2.momentum().t()*sqrtS();
106          z1 = (p1p2-p1.mass2()*p2.mass2()/p1p2)/(p2q-p2.mass2()*p1q/p1p2);
107          z2 = (p1p2-p1.mass2()*p2.mass2()/p1p2)/(p1q-p1.mass2()*p2q/p1p2);
108        }
109        _d_all[ip][def][0]->fill(z1,z2,0.5);
110        if (strong) _d_strong[ip][def][0]->fill(z1,z2,0.5);
111        if (same) {
112          _d_all[ip][def][1]->fill(z1,z2,0.5);
113          if (strong) _d_strong[ip][def][1]->fill(z1,z2,0.5);
114        }
115        if (opp) {
116          _d_all[ip][def][2]->fill(z1,z2,0.5);
117          if (strong) _d_strong[ip][def][2]->fill(z1,z2,0.5);
118        }
119      }
120    }
121
122    /// Perform the per-event analysis
123    void analyze(const Event& event) {
124      // apply projection
125      const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");
126      // fill single particle histos
127      for (const Particle& p : cfs.particles()) {
128        const double z = 2.*p.momentum().t()/sqrtS();
129        _s_all[p.abspid()]->fill(z);
130        if (!isWeak(p))  _s_strong[p.abspid()]->fill(z);
131      }
132      // get thrust
133      const Thrust thrust = apply<Thrust>(event,"Thrust");
134      ThreeVector axis = thrust.thrustAxis();
135      Particles piK = cfs.particles(Cuts::abspid==PID::KPLUS or Cuts::abspid==PID::PIPLUS);
136      for (size_t ix=0; ix<piK.size(); ++ix) {
137        double dot1 = axis.dot(piK[ix].momentum().p3());
138        bool weak1 = isWeak(piK[ix]);
139        for (size_t iy=0; iy<piK.size(); ++iy) {
140          if (ix==iy) continue;
141          double dot2 = axis.dot(piK[iy].momentum().p3());
142          bool weak2 = isWeak(piK[iy]);
143          bool strong = !weak1 && !weak2;
144          bool same = thrust.thrust()>0.8 && dot1*dot2>0.;
145          bool opp  = thrust.thrust()>0.8 && dot1*dot2<0.;
146          unsigned int ip=0;
147          if (piK[ix].pid()==PID::PIPLUS) {
148            if (piK[iy].pid()==PID::PIPLUS)       ip=1;
149            else if (piK[iy].pid()==PID::PIMINUS) ip=0;
150            else if (piK[iy].pid()==PID::KPLUS  ) ip=3;
151            else if (piK[iy].pid()==PID::KMINUS)  ip=2;
152          }
153          else if (piK[ix].pid()==PID::PIMINUS) {
154            if (piK[iy].pid()==PID::PIPLUS)       ip=0;
155            else if (piK[iy].pid()==PID::PIMINUS) ip=1;
156            else if (piK[iy].pid()==PID::KPLUS  ) ip=2;
157            else if (piK[iy].pid()==PID::KMINUS ) ip=3;
158          }
159          else if(piK[ix].pid()==PID::KPLUS) {
160            if(piK[iy].pid()==PID::PIPLUS)       ip=3;
161            else if(piK[iy].pid()==PID::PIMINUS) ip=2;
162            else if(piK[iy].pid()==PID::KPLUS)   ip=5;
163            else if(piK[iy].pid()==PID::KMINUS)  ip=4;
164          }
165          else if(piK[ix].pid()==PID::KMINUS) {
166            if(piK[iy].pid()==PID::PIPLUS)       ip=2;
167            else if(piK[iy].pid()==PID::PIMINUS) ip=3;
168            else if(piK[iy].pid()==PID::KPLUS)   ip=4;
169            else if(piK[iy].pid()==PID::KMINUS)  ip=5;
170          }
171          fillHistos(ip,strong,same,opp,piK[ix],piK[iy]);
172        }
173      }
174    }
175
176    /// Normalise histograms etc., after the run
177    void finalize() {
178
179      const double sf = crossSection()/femtobarn/sumOfWeights();
180      scale(_s_all, sf);
181      scale(_s_strong, sf);
182      for (size_t ix=0; ix<6; ++ix) {
183        for (size_t iy=0; iy<3; ++iy) {
184          scale(_d_all[ix][iy], sf);
185          divByGroupWidth(_d_all[ix][iy]);
186          scale(_d_strong[ix][iy], sf);
187          divByGroupWidth(_d_strong[ix][iy]);
188        }
189      }
190    }
191
192    /// @}
193
194
195    /// @name Histograms
196    /// @{
197    map<int,Histo1DPtr> _s_all,_s_strong;
198    Histo1DGroupPtr _d_all[6][3][3], _d_strong[6][3][3];
199    /// @}
200
201
202  };
203
204
205  RIVET_DECLARE_PLUGIN(BELLE_2020_I1777678);
206
207}