1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/UnstableParticles.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief gamma gamma -> p K+ pbar K-
 10  class BELLE_2016_I1444981 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2016_I1444981);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22      // Initialise and register projections
 23      declare(FinalState(), "FS");
 24      declare(UnstableParticles(Cuts::abspid==102134), "UFS");
 25      // counters
 26      for (unsigned int ix=0;ix<2;++ix) {
 27        book(_sigma[ix],"TMP/sigma_"+toString(ix+1),refData(ix+1, 1, 1));
 28      }
 29    }
 30
 31    void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
 32      for (const Particle &child : p.children()) {
 33        if (child.children().empty()) {
 34          nRes[child.pid()]-=1;
 35          --ncount;
 36        }
 37        else {
 38          findChildren(child,nRes,ncount);
 39        }
 40      }
 41    }
 42
 43    /// Perform the per-event analysis
 44    void analyze(const Event& event) {
 45      const FinalState& fs = apply<FinalState>(event, "FS");
 46      // find the final-state particles
 47      map<long,int> nCount;
 48      int ntotal(0);
 49      for (const Particle& p : fs.particles()) {
 50        nCount[p.pid()] += 1;
 51        ++ntotal;
 52      }
 53      // p pbar K+ K-
 54      if (ntotal==4 &&
 55          nCount[2212]==1 && nCount[-2212]==1 &&
 56          nCount[ 321]==1 && nCount[ -321]==1) _sigma[0]->fill(sqrtS());
 57      // intermediate Lambda(1520)
 58      for (const Particle& lam : apply<UnstableParticles>(event, "UFS").particles()) {
 59       	map<long,int> nRes=nCount;
 60       	int ncount = ntotal;
 61       	findChildren(lam,nRes,ncount);
 62        int sign = lam.pid()/lam.abspid();
 63        if (ncount!=2) continue;
 64        bool matched = true;
 65        for (const auto& val : nRes) {
 66          if (val.first==-sign*2212 || val.first==-sign*321) {
 67            if (val.second!=0) {
 68              matched = false;
 69              break;
 70            }
 71          }
 72          else if(val.second!=0) {
 73            matched = false;
 74            break;
 75          }
 76        }
 77        if (matched) {
 78          _sigma[1]->fill(sqrtS());
 79          break;
 80        }
 81      }
 82    }
 83
 84
 85    /// Normalise histograms etc., after the run
 86    void finalize() {
 87      const double fact = crossSection()/picobarn/sumOfWeights();
 88      // loop over tables in paper
 89      for(unsigned int ix=0;ix<2;++ix) {
 90        scale(_sigma[ix],fact);
 91        Estimate1DPtr tmp;
 92        book(tmp,ix+1,1,1);
 93        barchart(_sigma[ix],tmp);
 94      }
 95    }
 96
 97    /// @}
 98
 99
100    /// @name Histograms
101    /// @{
102    Histo1DPtr _sigma[2];
103    /// @}
104
105
106  };
107
108
109  RIVET_DECLARE_PLUGIN(BELLE_2016_I1444981);
110
111}