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- > sigma+- sigmabar -+
 10  class BESIII_2020_I1814783 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2020_I1814783);
 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      book(_n_plus ,"/TMP/NPLUS" );
 26      book(_n_minus,"/TMP/NMINUS");
 27      if(isCompatibleWithSqrtS(2.396, 1E-2)) {
 28	book(_h_cTheta_A,3,1,1);
 29	book(_h_cTheta_B,3,1,2);
 30      }
 31    }
 32
 33    void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
 34      for(const Particle &child : p.children()) {
 35	if(child.children().empty()) {
 36	  nRes[child.pid()]-=1;
 37	  --ncount;
 38	}
 39	else
 40	  findChildren(child,nRes,ncount);
 41      }
 42    }
 43
 44    /// Perform the per-event analysis
 45    void analyze(const Event& event) {
 46      const FinalState& fs = apply<FinalState>(event, "FS");
 47      // total hadronic and muonic cross sections
 48      map<long,int> nCount;
 49      int ntotal(0);
 50      for (const Particle& p : fs.particles()) {
 51	nCount[p.pid()] += 1;
 52	++ntotal;
 53      }
 54      // find the Sigmas
 55      const FinalState& ufs = apply<UnstableParticles>(event, "UFS");
 56      for(unsigned int ix=0;ix<ufs.particles().size();++ix) {
 57      	const Particle& p1 = ufs.particles()[ix];
 58      	if(abs(p1.pid())!=3112&&abs(p1.pid())!=3222) continue;
 59      	bool matched = false;
 60       	// check fs
 61       	bool fs = true;
 62      	for(const Particle & child : p1.children()) {
 63      	  if(child.pid()==p1.pid()) {
 64      	    fs = false;
 65      	    break;
 66      	  }
 67      	}
 68      	if(!fs) continue;
 69	// find the children
 70	map<long,int> nRes = nCount;
 71	int ncount = ntotal;
 72	findChildren(p1,nRes,ncount);
 73	for(unsigned int iy=ix+1;iy<ufs.particles().size();++iy) {
 74       	  const Particle& p2 = ufs.particles()[iy];
 75       	  if(p2.pid() != -p1.pid()) continue;
 76       	  // check fs
 77      	  bool fs = true;
 78      	  for(const Particle & child : p2.children()) {
 79      	    if(child.pid()==p2.pid()) {
 80      	      fs = false;
 81      	      break;
 82      	    }
 83      	  }
 84      	  if(!fs) continue;
 85      	  map<long,int> nRes2 = nRes;
 86      	  int ncount2 = ncount;
 87      	  findChildren(p2,nRes2,ncount2);
 88      	  if(ncount2!=0) continue;
 89      	  matched=true;
 90      	  for(auto const & val : nRes2) {
 91      	    if(val.second!=0) {
 92      	      matched = false;
 93      	      break;
 94      	    }
 95      	  }
 96      	  if(matched) {
 97	    if(abs(p1.pid())==3222) {
 98	      _n_plus->fill();
 99	      if(_h_cTheta_A) {
100		double cTheta = p1.pid()>0 ?
101		  cos(p1.momentum().polarAngle()) :
102		  cos(p2.momentum().polarAngle());
103		_h_cTheta_A->fill(cTheta);
104		_h_cTheta_B->fill(cTheta);
105	      }
106	    }
107	    else if(abs(p1.pid())==3112)
108	      _n_minus->fill();
109	    break;
110	  }
111      	}
112	if(matched) break;
113      }
114      
115    }
116
117
118    /// Normalise histograms etc., after the run
119    void finalize() {
120      if(_h_cTheta_A) {
121      	normalize(_h_cTheta_A);
122      	normalize(_h_cTheta_B);
123      }
124      double fact = crossSection()/ sumOfWeights() /picobarn;
125      for(unsigned int iy=1;iy<3;++iy) {
126	double sigma,error;
127	if(iy==1) {
128	  sigma = _n_plus->val()*fact;
129	  error = _n_plus->err()*fact;
130	}
131	else {
132	  sigma = _n_minus->val()*fact;
133	  error = _n_minus->err()*fact;
134	}
135	for(unsigned int ix=1;ix<3;++ix) {
136	  Scatter2D temphisto(refData(ix, 1, iy));
137	  Scatter2DPtr  mult;
138	  book(mult, ix, 1, iy);
139	  for (size_t b = 0; b < temphisto.numPoints(); b++) {
140	    const double x  = temphisto.point(b).x();
141	    pair<double,double> ex = temphisto.point(b).xErrs();
142	    pair<double,double> ex2 = ex;
143	    if(ex2.first ==0.) ex2. first=0.0001;
144	    if(ex2.second==0.) ex2.second=0.0001;
145	    if (inRange(sqrtS()/GeV, x-ex2.first, x+ex2.second)) {
146	      mult->addPoint(x, sigma, ex, make_pair(error,error));
147	    }
148	    else {
149	      mult->addPoint(x, 0., ex, make_pair(0.,.0));
150	    }
151	  }
152	}
153      }
154    }
155    ///@}
156
157
158    /// @name Histograms
159    ///@{
160    CounterPtr _n_plus,_n_minus;
161    Histo1DPtr _h_cTheta_A,_h_cTheta_B;
162    ///@}
163
164
165  };
166
167
168  RIVET_DECLARE_PLUGIN(BESIII_2020_I1814783);
169
170}