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    , 1, 1, 1);
 28      book(_c_DsStar, 2, 1, 1);
 29      for (const string& en : _c_Ds.binning().edges<0>()) {
 30        const double end = std::stod(en)*GeV;
 31        if (isCompatibleWithSqrtS(end)) {
 32          _ecms = en;
 33          break;
 34        }
 35      }
 36      if(_ecms.empty()) MSG_ERROR("Beam energy incompatible with analysis.");
 37    }
 38
 39    void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
 40      for(const Particle &child : p.children()) {
 41        if(child.children().empty()) {
 42          nRes[child.pid()]-=1;
 43          --ncount;
 44        }
 45        else
 46          findChildren(child,nRes,ncount);
 47      }
 48    }
 49
 50    /// Perform the per-event analysis
 51    void analyze(const Event& event) {
 52      const FinalState& fs = apply<FinalState>(event, "FS");
 53      // total analyse final state
 54      map<long,int> nCount;
 55      int ntotal(0);
 56      for (const Particle& p : fs.particles()) {
 57        nCount[p.pid()] += 1;
 58        ++ntotal;
 59      }
 60      // unstable charm analysis
 61      Particles ds = apply<UnstableParticles>(event, "UFS").particles(Cuts::abspid==431 or Cuts::abspid==433 or
 62                                                                      Cuts::abspid==20433);
 63      for(unsigned int ix=0;ix<ds.size();++ix) {
 64        const Particle& p1 = ds[ix];
 65        int id1 = abs(p1.pid());
 66        // check fs
 67        bool fs = true;
 68        for (const Particle & child : p1.children()) {
 69          if(child.pid()==p1.pid()) {
 70            fs = false;
 71            break;
 72          }
 73        }
 74        if(!fs) continue;
 75        // find the children
 76        map<long,int> nRes = nCount;
 77        int ncount = ntotal;
 78        findChildren(p1,nRes,ncount);
 79        bool matched=false;
 80        int sign = p1.pid()/id1;
 81        // loop over the other fs particles
 82        for(unsigned int iy=ix+1;iy<ds.size();++iy) {
 83          const Particle& p2 = ds[iy];
 84          fs = true;
 85          for (const Particle & child : p2.children()) {
 86            if(child.pid()==p2.pid()) {
 87              fs = false;
 88              break;
 89            }
 90          }
 91          if(!fs) continue;
 92          if(p2.pid()/abs(p2.pid())==sign) continue;
 93          int id2 = abs(p2.pid());
 94          if(!p2.parents().empty() && p2.parents()[0].pid()==p1.pid())
 95            continue;
 96          map<long,int> nRes2 = nRes;
 97          int ncount2 = ncount;
 98          findChildren(p2,nRes2,ncount2);
 99          if(ncount2!=0) continue;
100          matched=true;
101          for(auto const & val : nRes2) {
102            if(val.second!=0) {
103              matched = false;
104              break;
105            }
106          }
107          if(matched) {
108            if((id1==431 && id2==20433) || (id1==20433 && id2==421)) {
109              _c_Ds->fill(_ecms);
110            }
111            else if((id1==433 && id2==20433) || (id1==20433 && id2==433)) {
112              _c_DsStar->fill(_ecms);
113            }
114            break;
115          }
116        }
117        if(matched) break;
118      }
119    }
120
121
122    /// Normalise histograms etc., after the run
123    void finalize() {
124      double fact = crossSection()/ sumOfWeights()/picobarn;
125      scale(_c_Ds    ,fact);
126      scale(_c_DsStar,fact);
127    }
128
129    /// @}
130
131
132    /// @name Histograms
133    /// @{
134    BinnedHistoPtr<string> _c_Ds,_c_DsStar;
135    string _ecms;
136    /// @}
137
138  };
139
140
141  RIVET_DECLARE_PLUGIN(BESIII_2020_I1795949);
142
143}