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(*)+ Ds*()-
 10  class BELLE_2011_I878228 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2011_I878228);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22
 23      // Initialise and register projections
 24      declare(FinalState(), "FS");
 25      declare(UnstableParticles(), "UFS");
 26
 27      // Book histograms
 28      book(  _c["DpDm"], "/TMP/sigma_DpDm"  , refData(1,1,1));
 29      book( _c["DpDmS"], "/TMP/sigma_DpDmS" , refData(1,1,2));
 30      book(_c["DpSDmS"], "/TMP/sigma_DpSDmS", refData(1,1,3));
 31      book(   _c["All"], "/TMP/sigma_All"   , refData(2,1,1));
 32      book(    _c["mu"], "/TMP/sigma_mu"    , refData(2,1,1));
 33    }
 34
 35    void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
 36      for (const Particle &child : p.children()) {
 37        if (child.children().empty()) {
 38          nRes[child.pid()]-=1;
 39          --ncount;
 40        }
 41        else {
 42          findChildren(child,nRes,ncount);
 43        }
 44      }
 45    }
 46
 47    /// Perform the per-event analysis
 48    void analyze(const Event& event) {
 49      const FinalState& fs = apply<FinalState>(event, "FS");
 50      // total hadronic and muonic cross sections
 51      map<long,int> nCount;
 52      int ntotal(0);
 53      for (const Particle& p : fs.particles()) {
 54      	nCount[p.pid()] += 1;
 55      	++ntotal;
 56      }
 57      // mu+mu- + photons
 58      if (nCount[-13]==1 and nCount[13]==1 && ntotal==2+nCount[22]) {
 59        _c["mu"]->fill(sqrtS());
 60        return;
 61      }
 62      // unstable charm analysis
 63
 64      Particles ds = apply<UnstableParticles>(event, "UFS").particles(Cuts::abspid==431 || Cuts::abspid==433);
 65      for (unsigned int ix=0; ix<ds.size(); ++ix) {
 66       	const Particle& p1 = ds[ix];
 67        int id1 = abs(p1.pid());
 68      	// check fs
 69      	bool fs = true;
 70      	for (const Particle & child : p1.children()) {
 71      	  if (child.pid()==p1.pid()) {
 72      	    fs = false;
 73      	    break;
 74      	  }
 75      	}
 76      	if (!fs) continue;
 77      	// find the children
 78      	map<long,int> nRes = nCount;
 79      	int ncount = ntotal;
 80      	findChildren(p1,nRes,ncount);
 81      	bool matched=false;
 82       	int sign = p1.pid()/id1;
 83        // loop over the other fs particles
 84        for (unsigned int iy=ix+1; iy<ds.size(); ++iy) {
 85          const Particle& p2 = ds[iy];
 86          fs = true;
 87          for (const Particle & child : p2.children()) {
 88            if (child.pid()==p2.pid()) {
 89              fs = false;
 90              break;
 91            }
 92          }
 93          if (!fs) continue;
 94          if (p2.pid()/abs(p2.pid())==sign) continue;
 95          int id2 = abs(p2.pid());
 96          if (!p2.parents().empty() && p2.parents()[0].pid()==p1.pid()) {
 97            continue;
 98          }
 99          map<long,int> nRes2 = nRes;
100          int ncount2 = ncount;
101          findChildren(p2,nRes2,ncount2);
102          if (ncount2!=0) continue;
103          matched=true;
104          for (const auto& val : nRes2) {
105            if (val.second!=0) {
106              matched = false;
107              break;
108            }
109          }
110          if (matched) {
111            if (id1==431 && id2==431) {
112              _c["DpDm"]->fill(sqrtS());
113            }
114            else if (id1==433 && id2==433) {
115              _c["DpSDmS"]->fill(sqrtS());
116            }
117            else if ((id1==431 && id2==433) || (id1==433 && id2==431)) {
118              _c["DpDmS"]->fill(sqrtS());
119            }
120            _c["All"]->fill(sqrtS());
121            break;
122          }
123        }
124        if (matched) break;
125      }
126    }
127
128
129    /// Normalise histograms etc., after the run
130    void finalize() {
131      scale(_c, crossSection()/ sumOfWeights()/nanobarn);
132      size_t idx = 0;
133      for (const string& s : { "DpDm"s, "DpDmS"s, "DpSDmS"s}) {
134        Estimate1DPtr mult;
135        book(mult, 1, 1, ++idx);
136        barchart(_c[s],mult);
137      }
138      // R
139      Estimate1DPtr mult;
140      book(mult, 2,1,1);
141      divide(_c["All"s],_c["mu"s],mult);
142    }
143
144    /// @}
145
146
147    /// @name Histograms
148    /// @{
149    map<string,Histo1DPtr> _c;
150    /// @}
151
152
153  };
154
155
156  RIVET_DECLARE_PLUGIN(BELLE_2011_I878228);
157
158}