1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Tools/Cuts.hh"
  4#include "Rivet/Projections/UnstableParticles.hh"
  5
  6namespace Rivet {
  7
  8
  9  class PDG_UPSILON4S_HADRON_MULTIPLICITIES : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(PDG_UPSILON4S_HADRON_MULTIPLICITIES);
 14
 15
 16    /// @name Analysis methods
 17    /// @{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21
 22      // Initialise and register projections
 23      declare(UnstableParticles(), "UFS");
 24
 25      for(int ix : _ihistos)
 26	book(_histos[ix], ix, 1, 1);
 27      book(_wSum, "/TMP/SumWeights");
 28
 29    }
 30
 31
 32    /// Perform the per-event analysis
 33    void analyze(const Event& event) {
 34      for(const Particle& meson : apply<UnstableParticles>(event, "UFS").particles(Cuts::pid==300553)) {
 35	_wSum->fill();
 36	map<int,unsigned int> ncount;
 37	findDecayProducts(meson,ncount);
 38	_histos[29]->fill(10.579,double(ncount[411]+ncount[-411]));
 39	_histos[30]->fill(10.579,double(ncount[421]+ncount[-421]));
 40	_histos[31]->fill(10.579,double(ncount[413]+ncount[-413]));
 41	_histos[32]->fill(10.579,double(ncount[423]+ncount[-423]));
 42	_histos[33]->fill(10.579,double(ncount[431]+ncount[-431]));
 43	_histos[34]->fill(10.579,double(ncount[433]+ncount[-433]));
 44	_histos[48]->fill(10.579,double(ncount[443]));
 45	_histos[50]->fill(10.579,double(ncount[100443]));
 46	_histos[51]->fill(10.579,double(ncount[20443]));
 47	_histos[53]->fill(10.579,double(ncount[445]));
 48	_histos[60]->fill(10.579,double(ncount[321]+ncount[-321]));
 49	_histos[61]->fill(10.579,double(ncount[321]));
 50	_histos[62]->fill(10.579,double(ncount[-321]));
 51	_histos[63]->fill(10.579,double(ncount[130]+ncount[310]));
 52	_histos[64]->fill(10.579,double(ncount[323]+ncount[-323]));
 53	_histos[65]->fill(10.579,double(ncount[313]+ncount[-313]));
 54	_histos[87]->fill(10.579,double(ncount[211]+ncount[-211]));
 55	_histos[88]->fill(10.579,double(ncount[111]));
 56	_histos[89]->fill(10.579,double(ncount[221]));
 57	_histos[90]->fill(10.579,double(ncount[113]));
 58	_histos[92]->fill(10.579,double(ncount[333]));
 59	_histos[96]->fill(10.579,double(ncount[4122]+ncount[-4122]));
 60	_histos[104]->fill(10.579,double(ncount[-4222]));
 61	_histos[106]->fill(10.579,double(ncount[-4112]));
 62	_histos[110]->fill(10.579,double(ncount[2212]+ncount[-2212]));
 63	_histos[113]->fill(10.579,double(ncount[3122]+ncount[-3122]));
 64	_histos[116]->fill(10.579,double(ncount[3312]+ncount[-3312]));
 65      }
 66    }
 67
 68    void findDecayProducts(const Particle & mother,
 69			   map<int,unsigned int> & ncount) {
 70      for(const Particle & p : mother.children()) {
 71	int id = p.pid();
 72	if(p.children().empty()) {
 73	  ncount[id] += 1;
 74	}
 75	else {
 76	  // check particle is not a child or itself, eg copy or from photon radiation
 77	  bool isChild(false);
 78	  for(const Particle & p2 : p.children()) {
 79	    if(p2.pid()==id) {
 80	      isChild = true;
 81	      break;
 82	    }
 83	  }
 84	  if(!isChild) ncount[id] += 1;
 85	  findDecayProducts(p,ncount);
 86	}
 87      }
 88    }
 89
 90
 91    /// Normalise histograms etc., after the run
 92    void finalize() {
 93
 94      for(auto hist : _histos) {
 95	scale(hist.second, 1./_wSum->sumW());
 96      }
 97
 98    }
 99
100    /// @}
101
102
103    /// @name Histograms
104    /// @{
105    vector<int> _ihistos={29 ,30 ,31 ,32 ,33 ,34 ,48 ,50 ,51 ,53 ,60 ,61 ,62 ,63 ,
106			  64 ,65 ,87 ,88 ,89 ,90 ,92 ,96 ,104,106,110,113,116};
107    map<int,Histo1DPtr> _histos;
108    CounterPtr _wSum;
109    /// @}
110
111
112  };
113
114
115  RIVET_DECLARE_PLUGIN(PDG_UPSILON4S_HADRON_MULTIPLICITIES);
116
117}