1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief baryon rates and spectra
  9  class BELLE_2017_I1606201 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2017_I1606201);
 14
 15
 16    /// @name Analysis methods
 17    /// @{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21      // Initialise and register projections
 22      declare(UnstableParticles(), "UFS");
 23      // Book histograms
 24      for(unsigned int ix=1;ix<16;++ix) {
 25	book(_h[ix], ix, 1,  1);
 26	book(_r[ix], 16, 1, ix);
 27      }
 28    }
 29
 30
 31    /// Perform the per-event analysis
 32    void analyze(const Event& event) {
 33      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 34
 35      for (const Particle& p : ufs.particles()) {
 36	const Vector3 mom3 = p.p3();
 37	double pp = mom3.mod();
 38	double xp = pp/sqrt(0.25*sqr(sqrtS())-sqr(p.mass()));
 39	int id = abs(p.pid());
 40	if(id==3122) {
 41	  _h[ 1]->fill(xp);
 42	  _r[ 1]->fill(0.5);
 43	}
 44	else if(id==3212) {
 45	  _h[ 2]->fill(xp);
 46	  _r[ 2]->fill(0.5);
 47	}
 48	else if(id==3224) {
 49	  _h[ 3]->fill(xp);
 50	  _r[ 3]->fill(0.5);
 51	}
 52	else if(id==102134) {
 53	  _h[ 4]->fill(xp);
 54	  _r[ 4]->fill(0.5);
 55	}
 56	else if(id==3312) {
 57	  _h[ 5]->fill(xp);
 58	  _r[ 5]->fill(0.5);
 59	}
 60	else if(id==3334) {
 61	  _h[ 6]->fill(xp);
 62	  _r[ 6]->fill(0.5);
 63	}
 64	else if(id==3324) {
 65	  _h[ 7]->fill(xp);
 66	  _r[ 7]->fill(0.5);
 67	}
 68	else if(id==4122) {
 69	  _h[ 8]->fill(xp);
 70	  _r[ 8]->fill(0.5);
 71	}
 72	else if(id==102142) {
 73	  _h[ 9]->fill(xp);
 74	  _r[ 9]->fill(0.5);
 75	}
 76	else if(id==102144) {
 77	  _h[10]->fill(xp);
 78	  _r[10]->fill(0.5);
 79	}
 80	else if(id==4112) {
 81	  _h[11]->fill(xp);
 82	  _r[11]->fill(0.5);
 83	}
 84	else if(id==4114) {
 85	  _h[12]->fill(xp);
 86	  _r[12]->fill(0.5);
 87	}
 88	else if(id==4332) {
 89	  if(isDecay(p,{3334,-211})) {
 90	    _h[13]->fill(xp);
 91	    _r[13]->fill(0.5);
 92	  }
 93	}
 94	else if(id==4132) {
 95	  if(isDecay(p,{3312,211})) {
 96	    _h[14]->fill(xp);
 97	    _r[14]->fill(0.5);
 98	  }
 99	  else if(isDecay(p,{3334,321})) {
100	    _h[15]->fill(xp);
101	    _r[15]->fill(0.5);
102	  }
103	}
104      }
105    }
106
107    // Check for explicit decay into pdgids
108    bool isDecay(const Particle& mother, vector<int> ids) {
109      if(mother.pid()<0) {
110	for(unsigned int ix=0;ix<ids.size();++ix)
111	  ids[ix] *= -1;
112      }
113      // Trivial check to ignore any other decays but the one in question modulo photons
114      const Particles children = mother.children(Cuts::pid!=PID::PHOTON);
115      if (children.size()!=ids.size()) return false;
116      // Check for the explicit decay
117      return all(ids, [&](int i){return count(children, hasPID(i))==1;});
118    }
119
120    /// Normalise histograms etc., after the run
121    void finalize() {
122      // norm to cross section
123      for(unsigned int ix=1;ix<16;++ix) {
124	if( ix<=4 || (ix>=8 &&ix<=12) )
125	  scale(_h[ix], crossSection()/nanobarn/sumOfWeights());
126	else
127	  scale(_h[ix], crossSection()/picobarn/sumOfWeights());
128	scale(_r[ix], crossSection()/picobarn/sumOfWeights());
129      }
130    }
131
132    /// @}
133
134
135    /// @name Histograms
136    /// @{
137    Histo1DPtr _h[16],_r[16];
138    /// @}
139
140
141  };
142
143
144  RIVET_DECLARE_PLUGIN(BELLE_2017_I1606201);
145
146
147}