1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/UnstableParticles.hh"
  5#include "Rivet/Projections/Beam.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// @brief gamma gamma -> omega J/psi
 11  class BELLE_2010_I841003 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2010_I841003);
 16
 17
 18    /// @name Analysis methods
 19    /// @{
 20
 21    /// Book histograms and initialise projections before the run
 22    void init() {
 23      // Initialise and register projections
 24      declare(Beam(), "Beams");
 25      declare(FinalState(),"FS");
 26      declare(UnstableParticles(Cuts::pid==223 || Cuts::pid==443), "UFS");
 27      // histograms
 28      book(_h,1,1,1);
 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()];
 35          --ncount;
 36        } else {
 37          findChildren(child,nRes,ncount);
 38        }
 39      }
 40    }
 41    bool findScattered(Particle beam, double& q2) {
 42      bool found = false;
 43      Particle scat = beam;
 44      while (!scat.children().empty()) {
 45        found = false;
 46        for (const Particle & p : scat.children()) {
 47          if (p.pid()==scat.pid()) {
 48            scat=p;
 49            found=true;
 50            break;
 51          }
 52        }
 53        if (!found) break;
 54      }
 55      if (!found) return false;
 56      q2 = -(beam.momentum() - scat.momentum()).mass2();
 57      return true;
 58    }
 59
 60    /// Perform the per-event analysis
 61    void analyze(const Event& event) {
 62      // find scattered leptons and calc Q2
 63      const Beam& beams = apply<Beam>(event, "Beams");
 64      double q12 = -1, q22 = -1;
 65      if (!findScattered(beams.beams().first,  q12)) vetoEvent;
 66      if (!findScattered(beams.beams().second, q22)) vetoEvent;
 67      // check the final state
 68      const FinalState & fs = apply<FinalState>(event, "FS");
 69      map<long,int> nCount;
 70      int ntotal(0);
 71      for (const Particle& p : fs.particles()) {
 72        nCount[p.pid()] += 1;
 73        ++ntotal;
 74      }
 75      // find the meson
 76      const FinalState& ufs = apply<FinalState>(event, "UFS");
 77      for (const Particle& p1 : ufs.particles(Cuts::pid==223)) {
 78        if(p1.children().empty()) continue;
 79        bool matched=false;
 80        map<long,int> nRes = nCount;
 81        int ncount = ntotal;
 82        findChildren(p1,nRes,ncount);
 83        for (const Particle& p2 : ufs.particles(Cuts::pid==443)) {
 84          map<long,int> nRes2 = nRes;
 85          int ncount2 = ncount;
 86          findChildren(p2,nRes2,ncount2);
 87          matched = true;
 88          for(const auto& val : nRes2) {
 89            if (abs(val.first)==11) {
 90              if (val.second!=1) {
 91                matched = false;
 92                break;
 93              }
 94            }
 95            else if(val.second!=0) {
 96              matched = false;
 97              break;
 98            }
 99          }
100          if (matched) {
101            _h->fill((p1.momentum()+p2.momentum()).mass());
102            break;
103          }
104        }
105        if (matched) break;
106      }
107    }
108
109
110    /// Normalise histograms etc., after the run
111    void finalize() {
112      normalize(_h, 1.0, false);
113    }
114
115    /// @}
116
117
118    /// @name Histograms
119    /// @{
120    Histo1DPtr _h;
121    /// @}
122
123
124  };
125
126
127  RIVET_DECLARE_PLUGIN(BELLE_2010_I841003);
128
129}