1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief form factors for omega->pi and eta->gamma
  9  class A2_2017_I1486671 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(A2_2017_I1486671);
 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      book(_h_eta  , 1, 1, 1);
 25      book(_h_omega, 2, 1, 1);
 26      book(_weight_eta  ,"TMP/weight_eta");
 27      book(_weight_omega,"TMP/weight_omega");
 28    }
 29
 30    void findDecayProducts(const Particle & mother, unsigned int & nstable, unsigned int & npi,
 31                           unsigned int & nep, unsigned int & nem, unsigned int & ngamma,
 32			   FourMomentum & ptot) {
 33      for(const Particle & p : mother.children()) {
 34        int id = p.pid();
 35        if (id == PID::EMINUS ) {
 36	  ++nem;
 37          ++nstable;
 38	  ptot += p.momentum();
 39	}
 40        else if (id == PID::EPLUS) {
 41          ++nep;
 42          ++nstable;
 43	  ptot += p.momentum();
 44        }
 45        else if (id == PID::PI0) {
 46	  ++npi;
 47          ++nstable;
 48        }
 49	else if (id == PID::GAMMA && p.children().empty() ) {
 50	  ++ngamma;
 51	  ++nstable;
 52	}
 53        else if ( !p.children().empty() ) {
 54          findDecayProducts(p, nstable, npi,nep,nem,ngamma,ptot);
 55        }
 56        else
 57          ++nstable;
 58      }
 59    }
 60
 61    /// Perform the per-event analysis
 62    void analyze(const Event& event) {
 63      static double me     = 0.5109989461*MeV;
 64      static double momega = 782.65*MeV;
 65      static double meta   = 547.862 *MeV;
 66      static double mpi    = 134.9770*MeV;
 67
 68      // Loop over eta and omega mesons
 69      for( const Particle& p : apply<UnstableParticles>(event, "UFS").particles(Cuts::pid==223 or Cuts::pid==221)) {
 70       	unsigned nstable(0),npi(0),nep(0),nem(0),ngamma(0);
 71       	FourMomentum ptot;
 72      	findDecayProducts(p,nstable,npi,nep,nem,ngamma,ptot);
 73	if(p.pid()==221) {
 74	  if(nstable==3 && nem==1 && nem==1 && ngamma==1) {
 75	    double q = ptot.mass();
 76	    double beta = sqrt(1.-4*sqr(me/q));
 77	    double p = 1.-sqr(q/meta);
 78	    double fact = beta*MeV/q*(1.+2.*sqr(me/q))*pow(p,3);
 79	    _h_eta->fill(q/MeV,1./fact);
 80	  }
 81	  else if(nstable==2 && ngamma==2) {
 82	    _weight_eta->fill();
 83	  }
 84	}
 85	else {
 86	  if(nstable==3 && nem==1 && nem==1 && npi==1) {
 87	    double q = ptot.mass();
 88	    double beta = sqrt(1.-4*sqr(me/q));
 89	    double p = sqrt(sqr(1.+sqr(q)/(sqr(momega)-sqr(mpi)))-4.*sqr(momega*q/(sqr(momega)-sqr(mpi))));
 90	    double fact = beta*MeV/q*(1.+2.*sqr(me/q))*pow(p,3);
 91	    _h_omega->fill(q/MeV,1./fact);
 92	  }
 93	  else if(nstable==2 && ngamma ==1 && npi==1) {
 94	    _weight_omega->fill();
 95	  }
 96	}
 97      }
 98    }
 99
100
101    /// Normalise histograms etc., after the run
102    void finalize() {
103      static double alpha= 7.2973525664e-3;
104      scale(_h_eta  , 0.75*M_PI/alpha/ *_weight_eta  );
105      scale(_h_omega, 1.5 *M_PI/alpha/ *_weight_omega);
106    }
107
108    /// @}
109
110
111    /// @name Histograms
112    /// @{
113    Histo1DPtr _h_eta,_h_omega;
114    CounterPtr _weight_eta,_weight_omega;
115    /// @}
116
117
118  };
119
120
121  RIVET_DECLARE_PLUGIN(A2_2017_I1486671);
122
123
124}