1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/UnstableParticles.hh"
  5
  6
  7namespace Rivet {
  8
  9
 10  /// @brief gamma gamma -> Lambda Lmabdabar, Sigma0 Sigmabar0
 11  class BELLE_2006_I727063 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2006_I727063);
 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(FinalState(), "FS");
 25      declare(UnstableParticles(Cuts::abspid==3122 ||	Cuts::abspid==3212), "UFS");
 26      // histograms
 27      for(unsigned int ix=0; ix<2; ++ix) {
 28        book(_sigma[ix],"TMP/c_"+toString(ix+1),refData(1,1,1+ix));
 29      }
 30    }
 31
 32    void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
 33      for (const Particle &child : p.children()) {
 34        if (child.children().empty()) {
 35          --nRes[child.pid()];
 36          --ncount;
 37        } else {
 38          findChildren(child,nRes,ncount);
 39        }
 40      }
 41    }
 42
 43    /// Perform the per-event analysis
 44    void analyze(const Event& event) {
 45      // check the final state
 46      const FinalState & fs = apply<FinalState>(event, "FS");
 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      const FinalState& ufs = apply<FinalState>(event, "UFS");
 54      Particle Lambda,Lambar;
 55      // loop over baryon
 56      bool matched = false;
 57      for (const Particle& p : ufs.particles()) {
 58        if(p.children().empty() || p.pid()<0) continue;
 59        map<long,int> nRes = nCount;
 60        int ncount = ntotal;
 61        findChildren(p,nRes,ncount);
 62      	// and antibaryons
 63      	for (const Particle& p2 : ufs.particles()) {
 64      	  if(p2.children().empty() || p2.pid()>0) continue;
 65      	  map<long,int> nRes2 = nRes;
 66      	  int ncount2 = ncount;
 67      	  findChildren(p2,nRes2,ncount2);
 68      	  if(ncount2!=0) continue;
 69      	  matched=2;
 70      	  for(auto const & val : nRes2) {
 71      	    if(val.second!=0) {
 72      	      matched = false;
 73      	      break;
 74      	    }
 75      	  }
 76      	  if(matched) {
 77      	    Lambda=p;
 78      	    Lambar=p2;
 79      	    break;
 80      	  }
 81      	}
 82      }
 83      if(!matched) vetoEvent;
 84      if(Lambda.pid()!=-Lambar.pid()) vetoEvent;
 85      double cTheta = abs(Lambda.momentum().z()/Lambda.momentum().p3().mod());
 86      if(cTheta>0.6) vetoEvent;
 87      if(Lambda.pid()==3122) _sigma[1]->fill(sqrtS());
 88      else                   _sigma[0]->fill(sqrtS());
 89    }
 90
 91
 92    /// Normalise histograms etc., after the run
 93    void finalize() {
 94      double fact = crossSection()/nanobarn/sumOfWeights();
 95      for(unsigned int ix=0;ix<2;++ix) {
 96        scale(_sigma[ix],fact);
 97        Estimate1DPtr tmp;
 98        book(tmp, 1, 1, 1+ix);
 99        barchart(_sigma[ix],tmp);
100      }
101    }
102
103    /// @}
104
105
106    /// @name Histograms
107    /// @{
108    Histo1DPtr _sigma[2];
109    /// @}
110
111
112  };
113
114
115  RIVET_DECLARE_PLUGIN(BELLE_2006_I727063);
116
117}