1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/UnstableParticles.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief e+e- -> Ds(*) D_s1
 10  class BESIII_2020_I1795949 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2020_I1795949);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22      // Initialise and register projections
 23      declare(FinalState(), "FS");
 24      declare(UnstableParticles(), "UFS");
 25
 26      // Book histograms
 27      book(_c_Ds    ,"/TMP/c_Ds"    );
 28      book(_c_DsStar,"/TMP/c_DsStar");
 29    }
 30
 31    void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
 32      for(const Particle &child : p.children()) {
 33        if(child.children().empty()) {
 34          nRes[child.pid()]-=1;
 35          --ncount;
 36        }
 37        else
 38          findChildren(child,nRes,ncount);
 39      }
 40    }
 41
 42    /// Perform the per-event analysis
 43    void analyze(const Event& event) {
 44      const FinalState& fs = apply<FinalState>(event, "FS");
 45      // total analyse final state
 46      map<long,int> nCount;
 47      int ntotal(0);
 48      for (const Particle& p : fs.particles()) {
 49        nCount[p.pid()] += 1;
 50        ++ntotal;
 51      }
 52      // unstable charm analysis
 53      Particles ds = apply<UnstableParticles>(event, "UFS").particles(Cuts::abspid==431 or Cuts::abspid==433 or
 54                                                                      Cuts::abspid==20433);
 55      for(unsigned int ix=0;ix<ds.size();++ix) {
 56        const Particle& p1 = ds[ix];
 57        int id1 = abs(p1.pid());
 58        // check fs
 59        bool fs = true;
 60        for (const Particle & child : p1.children()) {
 61          if(child.pid()==p1.pid()) {
 62            fs = false;
 63            break;
 64          }
 65        }
 66        if(!fs) continue;
 67        // find the children
 68        map<long,int> nRes = nCount;
 69        int ncount = ntotal;
 70        findChildren(p1,nRes,ncount);
 71        bool matched=false;
 72        int sign = p1.pid()/id1;
 73        // loop over the other fs particles
 74        for(unsigned int iy=ix+1;iy<ds.size();++iy) {
 75          const Particle& p2 = ds[iy];
 76          fs = true;
 77          for (const Particle & child : p2.children()) {
 78            if(child.pid()==p2.pid()) {
 79              fs = false;
 80              break;
 81            }
 82          }
 83          if(!fs) continue;
 84          if(p2.pid()/abs(p2.pid())==sign) continue;
 85          int id2 = abs(p2.pid());
 86          if(!p2.parents().empty() && p2.parents()[0].pid()==p1.pid())
 87            continue;
 88          map<long,int> nRes2 = nRes;
 89          int ncount2 = ncount;
 90          findChildren(p2,nRes2,ncount2);
 91          if(ncount2!=0) continue;
 92          matched=true;
 93          for(auto const & val : nRes2) {
 94            if(val.second!=0) {
 95              matched = false;
 96              break;
 97            }
 98          }
 99          if(matched) {
100            if((id1==431 && id2==20433) || (id1==20433 && id2==421)) {
101              _c_Ds->fill();
102            }
103            else if((id1==433 && id2==20433) || (id1==20433 && id2==433)) {
104              _c_DsStar->fill();
105            }
106            break;
107          }
108        }
109        if(matched) break;
110      }
111    }
112
113
114    /// Normalise histograms etc., after the run
115    void finalize() {
116      double fact = crossSection()/ sumOfWeights()/picobarn;
117      for (unsigned int ih=1;ih<3;++ih) {
118        double sigma = 0.0, error = 0.0;
119        if(ih==1) {
120          sigma = _c_Ds->val()*fact;
121          error = _c_Ds->err()*fact;
122        }
123        else if(ih==2) {
124          sigma = _c_DsStar->val()*fact;
125          error = _c_DsStar->err()*fact;
126        }
127        Estimate1DPtr mult;
128        book(mult, ih, 1, 1);
129        for (auto& b : mult->bins()) {
130          if (inRange(sqrtS()/GeV, b.xMin(), b.xMax())) {
131            b.set(sigma, error);
132          }
133        }
134      }
135    }
136
137    /// @}
138
139
140    /// @name Histograms
141    /// @{
142    CounterPtr _c_Ds,_c_DsStar;
143    /// @}
144
145  };
146
147
148  RIVET_DECLARE_PLUGIN(BESIII_2020_I1795949);
149
150}