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	Scatter2D temphisto(refData(ih, 1, 1));
128        Scatter2DPtr     mult;
129        book(mult, ih, 1, 1);
130        for (size_t b = 0; b < temphisto.numPoints(); b++) {
131          const double x  = temphisto.point(b).x();
132          pair<double,double> ex = temphisto.point(b).xErrs();
133          pair<double,double> ex2 = ex;
134          if(ex2.first ==0.) ex2. first=0.0001;
135          if(ex2.second==0.) ex2.second=0.0001;
136          if (inRange(sqrtS()/GeV, x-ex2.first, x+ex2.second)) {
137            mult   ->addPoint(x, sigma, ex, make_pair(error,error));
138          }
139          else {
140            mult   ->addPoint(x, 0., ex, make_pair(0.,.0));
141          }
142        }
143      }
144    }
145
146    ///@}
147
148
149    /// @name Histograms
150    ///@{
151    CounterPtr _c_Ds,_c_DsStar;
152    ///@}
153
154
155  };
156
157
158  RIVET_DECLARE_PLUGIN(BESIII_2020_I1795949);
159
160}