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