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 DM2_1990_I297706 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(DM2_1990_I297706);
 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(_nProton,  "/TMP/nProton" );
 26      book(_nLambda,  "/TMP/nLambda" );
 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()]-=1;
 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      // total hadronic and muonic cross sections
 44      map<long,int> nCount;
 45      int ntotal(0);
 46      for (const Particle& p : fs.particles()) {
 47	nCount[p.pid()] += 1;
 48	++ntotal;
 49      }
 50      if(ntotal==2 && nCount[2212]==1 && nCount[-2212]==1)
 51	_nProton->fill();
 52      
 53      // find the Lambdas
 54      const FinalState& ufs = apply<UnstableParticles>(event, "UFS");
 55      for(unsigned int ix=0;ix<ufs.particles().size();++ix) {
 56	const Particle& p1 = ufs.particles()[ix];
 57	if(abs(p1.pid())!=3122) continue;
 58	bool matched = false;
 59	// check fs
 60	bool fs = true;
 61	for (const Particle & child : p1.children()) {
 62	  if(child.pid()==p1.pid()) {
 63	    fs = false;
 64	    break;
 65	  }
 66	}
 67	if(!fs) continue;
 68	// find the children
 69	map<long,int> nRes = nCount;
 70	int ncount = ntotal;
 71	findChildren(p1,nRes,ncount);
 72	for(unsigned int iy=ix+1;iy<ufs.particles().size();++iy) {
 73	  const Particle& p2 = ufs.particles()[iy];
 74	  if(abs(p2.pid())!=3122) continue;
 75	  // check fs
 76	  bool 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	  map<long,int> nRes2 = nRes;
 85	  int ncount2 = ncount;
 86	  findChildren(p2,nRes2,ncount2);
 87	  if(ncount2!=0) continue;
 88	  matched=true;
 89	  for(auto const & val : nRes2) {
 90	    if(val.second!=0) {
 91	      matched = false;
 92	      break;
 93	    }
 94	  }
 95	  if(matched) {
 96	    _nLambda->fill();
 97	    break;
 98	  }
 99	}
100	if(matched) break;
101      }
102    }
103
104
105    /// Normalise histograms etc., after the run
106    void finalize() {
107      double fact = crossSection()/ sumOfWeights() /picobarn;
108      for(unsigned int ix=1;ix<3;++ix) {
109	double sigma,error;
110	if(ix==1) {
111	  sigma = _nProton->val()*fact;
112	  error = _nProton->err()*fact;
113	}
114	else {
115	  sigma = _nLambda->val()*fact;
116	  error = _nLambda->err()*fact;
117	}
118	Scatter2D temphisto(refData(ix, 1, 1));
119	Scatter2DPtr mult;
120	book(mult, ix, 1, 1);
121	for (size_t b = 0; b < temphisto.numPoints(); b++) {
122	  const double x  = temphisto.point(b).x();
123	  pair<double,double> ex = temphisto.point(b).xErrs();
124	  pair<double,double> ex2 = ex;
125	  if(ex2.first ==0.) ex2. first=0.0001;
126	  if(ex2.second==0.) ex2.second=0.0001;
127	  if (inRange(sqrtS()/GeV, x-ex2.first, x+ex2.second)) {
128	    mult->addPoint(x, sigma, ex, make_pair(error,error));
129	  }
130	  else {
131	    mult->addPoint(x, 0., ex, make_pair(0.,.0));
132	  }
133	}
134      }
135    }
136    //@}
137
138    /// @name Histograms
139    //@{
140    CounterPtr _nProton,_nLambda;
141    //@}
142
143  };
144
145
146  // The hook for the plugin system
147  RIVET_DECLARE_PLUGIN(DM2_1990_I297706);
148
149
150}