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 Add a short analysis description here
 10  class BESIII_2015_I1329785 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2015_I1329785);
 15
 16
 17    /// @name Analysis methods
 18    //@{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22      declare(FinalState(), "FS");
 23      declare(UnstableParticles(), "UFS");
 24      book(_nChi0, "TMP/chi0");
 25      book(_nChi1, "TMP/chi1");
 26      book(_nChi2, "TMP/chi2");
 27    }
 28
 29    void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
 30      for( const Particle &child : p.children()) {
 31	if(child.children().empty()) {
 32	  --nRes[child.pid()];
 33	  --ncount;
 34	}
 35	else
 36	  findChildren(child,nRes,ncount);
 37      }
 38    }
 39
 40    /// Perform the per-event analysis
 41    void analyze(const Event& event) {
 42      const FinalState& fs = apply<FinalState>(event, "FS");
 43      map<long,int> nCount;
 44      int ntotal(0);
 45      for (const Particle& p :  fs.particles()) {
 46	nCount[p.pid()] += 1;
 47	++ntotal;
 48      }
 49      const FinalState& ufs = apply<FinalState>(event, "UFS");
 50      for (const Particle& p :  ufs.particles(Cuts::pid==10441 or Cuts::pid==20443 or Cuts::pid==445)) {
 51	if(p.children().empty()) continue;
 52	map<long,int> nRes = nCount;
 53	int ncount = ntotal;
 54	findChildren(p,nRes,ncount);
 55	// chi gamma 
 56	if(ncount!=1) continue;
 57	bool matched = true;
 58	for(auto const & val : nRes) {
 59	  if(val.first==22) {
 60	    if(val.second !=1) {
 61	      matched = false;
 62	      break;
 63	    }
 64	  }
 65	  else if(val.second!=0) {
 66	    matched = false;
 67	    break;
 68	  }
 69	}
 70	if(matched) {
 71	  if(p.pid()==10441)
 72	    _nChi0->fill();
 73	  else if(p.pid()==20443)
 74	    _nChi1->fill();
 75	  else if(p.pid()==445)
 76	    _nChi2->fill();
 77	  break;
 78	}
 79      }
 80    }
 81
 82
 83    /// Normalise histograms etc., after the run
 84    void finalize() {
 85
 86      double fact =  crossSection()/ sumOfWeights() /picobarn;
 87      for(unsigned int ix=1;ix<4;++ix) {
 88	double sigma(0.),error(0.);
 89	if(ix==1) {
 90	  sigma = _nChi0->val()*fact;
 91	  error = _nChi0->err()*fact;
 92	}
 93	else if(ix==2) {
 94	  sigma = _nChi1->val()*fact;
 95	  error = _nChi1->err()*fact;
 96	}
 97	else if(ix==3) {
 98	  sigma = _nChi2->val()*fact;
 99	  error = _nChi2->err()*fact;
100	}
101	Scatter2D temphisto(refData(ix, 1, 8));
102	Scatter2DPtr  mult;
103	book(mult,ix, 1, 8);
104	for (size_t b = 0; b < temphisto.numPoints(); b++) {
105	  const double x  = temphisto.point(b).x();
106	  pair<double,double> ex = temphisto.point(b).xErrs();
107	  pair<double,double> ex2 = ex;
108	  if(ex2.first ==0.) ex2. first=0.0001;
109	  if(ex2.second==0.) ex2.second=0.0001;
110	  if (inRange(sqrtS()/GeV, x-ex2.first, x+ex2.second)) {
111	    mult->addPoint(x, sigma, ex, make_pair(error,error));
112	  }
113	  else {
114	    mult->addPoint(x, 0., ex, make_pair(0.,.0));
115	  }
116	}
117      }
118    }
119    //@}
120
121    /// @name Histograms
122    //@{
123    CounterPtr _nChi0,_nChi1,_nChi2;
124    //@}
125
126  };
127
128
129  // The hook for the plugin system
130  RIVET_DECLARE_PLUGIN(BESIII_2015_I1329785);
131
132
133}