1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief e+ e- -> Ds+ Dbar0(*) K-
  9  class BESIII_2018_I1711382 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2018_I1711382);
 14
 15
 16    /// @name Analysis methods
 17    /// @{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21      // projections
 22      declare(FinalState(), "FS");
 23      // histograms
 24      for (unsigned int ix=0;ix<4;++ix) {
 25        book(_h_sigma[ix], 1, 1, 1+ix);
 26      }
 27      book(_h_cTheta, 2, 1, 1);
 28    }
 29
 30
 31    /// Perform the per-event analysis
 32    void analyze(const Event& event) {
 33      Particles fs = apply<FinalState>(event, "FS").particles();
 34      Particles DS,D0,other;
 35      // first find the kaon and Ds
 36      for (const Particle& p : fs) {
 37        Particle parent=p;
 38        while (!parent.parents().empty()) {
 39          parent=parent.parents()[0];
 40          if (parent.abspid()==PID::DSPLUS ||
 41              parent.abspid()==PID::D0) break;
 42        }
 43        if (parent.abspid()!=PID::DSPLUS &&
 44            parent.abspid()!=PID::D0 ) {
 45          other.push_back(p);
 46          continue;
 47        }
 48        bool found=false;
 49        for (auto& D : parent.abspid()==PID::DSPLUS ? DS : D0) {
 50          // D already in list
 51          if (fuzzyEquals(D.mom(),parent.mom())) {
 52            found=true;
 53            break;
 54          }
 55        }
 56        if(!found) {
 57          (parent.abspid()==PID::DSPLUS ? DS : D0).push_back(parent);
 58        }
 59      }
 60      // Ds and D0, one particle and one anti particle
 61      if (DS.size()!=1 || D0.size()!=1 || DS[0].pid()*D0[0].pid()>0) vetoEvent;
 62      int iK = DS[0].pid()>0 ? -321 : 321;
 63      if (other.size()==1 && other[0].pid()==iK) {
 64        _h_sigma[1]->fill("4.6"s);
 65        if (other[0].parents()[0].abspid()==435 &&
 66           D0   [0].parents()[0].abspid()==435) {
 67          _h_sigma[3]->fill("4.6"s);
 68          FourMomentum pDs2 = D0[0].momentum()+other[0].momentum();
 69          LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(pDs2.betaVec());
 70          Vector3 axis1 = pDs2.p3().unit();
 71          Vector3 axis2 = boost.transform(other[0].momentum()).p3().unit();
 72          _h_cTheta->fill(abs(axis1.dot(axis2)));
 73        }
 74      }
 75      else if (D0[0].parents()[0].abspid()==423) {
 76        bool Dstar = true;
 77        for (const Particle& p : other) {
 78          if (p.parents()[0].abspid()!=423 && p.pid()!=iK) {
 79            Dstar = false;
 80            break;
 81          }
 82        }
 83        if (Dstar) {
 84          _h_sigma[0]->fill("4.6"s);
 85          if (D0[0].parents()[0].parents()[0].abspid()==10433) {
 86            _h_sigma[2]->fill("4.6"s);
 87          }
 88        }
 89      }
 90    }
 91
 92    /// Normalise histograms etc., after the run
 93    void finalize() {
 94      scale(_h_sigma, crossSection()/ sumOfWeights() /picobarn);
 95      normalize(_h_cTheta);
 96    }
 97
 98    /// @}
 99
100
101    /// @name Histograms
102    /// @{
103    BinnedHistoPtr<string> _h_sigma[4];
104    Histo1DPtr _h_cTheta;
105    /// @}
106
107
108  };
109
110
111  RIVET_DECLARE_PLUGIN(BESIII_2018_I1711382);
112
113}