1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief chi_b(2S) decays
  9  class CLEOII_1992_I32611 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(CLEOII_1992_I32611);
 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      // averages
 25      book(_h_N_aver ,1,1,1);
 26      book(_h_R2_aver,1,1,2);
 27      // dists
 28      _h_N  = {Histo1DPtr(),Histo1DPtr(),Histo1DPtr()};
 29      _h_R2 = {Histo1DPtr(),Histo1DPtr(),Histo1DPtr()};
 30      book(_h_N [0],2,1,1);
 31      book(_h_N [1],2,1,2);
 32      book(_h_N [2],2,1,3);
 33      book(_h_R2[0],3,1,1);
 34      book(_h_R2[1],3,1,2);
 35      book(_h_R2[2],3,1,3);
 36    }
 37
 38    void findDecayProducts(Particle parent, Particles & children, unsigned int & nCharged) {
 39      for(const Particle & p: parent.children()) {
 40	if(p.children().empty()) {
 41	  if(isCharged(p)) ++nCharged;
 42	  children.push_back(p);
 43	}
 44	else
 45	  findDecayProducts(p,children,nCharged);
 46      }
 47    }
 48
 49    /// Perform the per-event analysis
 50    void analyze(const Event& event) {
 51      Particles chib = apply<UnstableParticles>(event,"UFS").particles(Cuts::pid==110551 or
 52								       Cuts::pid==120553 or
 53								       Cuts::pid==100555);
 54      for(const Particle & p : chib) {
 55	unsigned int iHist = (p.pid()%10)/2;
 56	unsigned int nCharged(0);
 57	Particles children;
 58	findDecayProducts(p,children,nCharged);
 59	// ncharged
 60	_h_N     [iHist]->fill(nCharged);
 61	_h_N_aver->fill(iHist,nCharged);
 62	// R_2
 63	LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
 64	vector<FourMomentum> mom;
 65	mom.reserve(children.size());
 66	for(const Particle & p2: children) {
 67	  mom.push_back(boost.transform(p2.momentum()));
 68	}
 69	// compute R2
 70	double H0(0.),H2(0.);
 71	for(const FourMomentum & p1:mom) {
 72	  double mod1 = p1.p3().mod();
 73	  Vector3 axis = p1.p3().unit();
 74	  for(const FourMomentum & p2:mom) {
 75	    double mod2 = p2.p3().mod();
 76	    double cTheta = axis.dot(p2.p3().unit());
 77	    H0 += mod1*mod2;
 78	    H2 += mod1*mod2*0.5*(3.*sqr(cTheta)-1.);
 79	  }
 80	}
 81	double R2=H2/H0;
 82	_h_R2     [iHist]->fill(R2);
 83	_h_R2_aver->fill(iHist,R2);
 84      }
 85    }
 86
 87
 88    /// Normalise histograms etc., after the run
 89    void finalize() {
 90      for(unsigned int ix=0;ix<3;++ix) {
 91	normalize( _h_N [ix]);
 92	normalize( _h_R2[ix]);
 93      }
 94    }
 95
 96    /// @}
 97
 98
 99    /// @name Histograms
100    /// @{
101    vector<Histo1DPtr> _h_N,_h_R2;
102    BinnedProfilePtr<int> _h_N_aver,_h_R2_aver;
103    /// @}
104
105
106  };
107
108
109  RIVET_DECLARE_PLUGIN(CLEOII_1992_I32611);
110
111}