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