1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5#include "Rivet/Projections/LeptonFinder.hh"
  6#include "Rivet/Projections/MissingMomentum.hh"
  7#include "Rivet/Projections/PromptFinalState.hh"
  8#include "Rivet/Projections/UnstableParticles.hh"
  9
 10namespace Rivet {
 11
 12
 13  /// @brief Charm-hadron production
 14  class BELLE_2005_I686014 : public Analysis {
 15  public:
 16
 17    /// Constructor
 18    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2005_I686014);
 19
 20
 21    /// @name Analysis methods
 22    /// @{
 23
 24    /// Book histograms and initialise projections before the run
 25    void init() {
 26      // projections
 27      declare(UnstableParticles(),"UFS");
 28
 29      // histos
 30      if (isCompatibleWithSqrtS(10.52*GeV)) _mode=1;
 31      else if (isCompatibleWithSqrtS(10.58*GeV,1e-4)) _mode=2;
 32      else MSG_ERROR("Beam energy not supported!");
 33
 34      for (unsigned int ix=0;ix<7;++ix) {
 35        if (_mode==1)
 36          book(_r[ix],2,1,ix+1);
 37        else
 38          book(_r[ix],1,1,ix+1);
 39        book(_h[ix],2+_mode,1,ix+1);
 40      }
 41      book(_c,"TMP/wgt");
 42    }
 43
 44    /// Perform the per-event analysis
 45    void analyze(const Event& event) {
 46      // unstable particles
 47      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 48      if (_mode==2 && ufs.particles(Cuts::pid==300553).size()!=1)
 49        vetoEvent;
 50      _c->fill();
 51      for (const Particle & p : ufs.particles(Cuts::abspid==411 ||
 52                                              Cuts::abspid==421 ||
 53                                              Cuts::abspid==431 ||
 54                                              Cuts::abspid==413 ||
 55                                              Cuts::abspid==423 ||
 56                                              Cuts::abspid==4122 )) {
 57        double pmax = sqrt(0.25*sqr(sqrtS())-sqr(p.mass()));
 58        double xp = p.momentum().p3().mod()/pmax;
 59        if (p.abspid()==421) {
 60          _r[0]->fill(0.5);
 61          _h[0]->fill(xp);
 62        }
 63        else if (p.abspid()==421) {
 64          _r[0]->fill(0.5);
 65          _h[0]->fill(xp);
 66        }
 67        else if (p.abspid()==411) {
 68          _r[1]->fill(0.5);
 69          _h[1]->fill(xp);
 70        }
 71        else if (p.abspid()==431) {
 72          _r[2]->fill(0.5);
 73          _h[2]->fill(xp);
 74        }
 75        else if (p.abspid()==4122) {
 76          _r[3]->fill(0.5);
 77          _h[3]->fill(xp);
 78        }
 79        else if (p.abspid()==413) {
 80          _r[4]->fill(0.5);
 81          _h[4]->fill(xp);
 82          _r[5]->fill(0.5);
 83          _h[5]->fill(xp);
 84        }
 85        else if (p.abspid()==423) {
 86          _r[6]->fill(0.5);
 87          _h[6]->fill(xp);
 88        }
 89      }
 90    }
 91
 92
 93    /// Normalise histograms etc., after the run
 94    void finalize() {
 95      if (_mode==1) {
 96        for (unsigned int ix=0;ix<7;++ix) {
 97          scale(_r[ix],crossSection()/picobarn/sumOfWeights());
 98          scale(_h[ix],crossSection()/nanobarn/sumOfWeights());
 99        }
100      }
101      else {
102        for (unsigned int ix=0;ix<7;++ix) {
103          scale(_r[ix], 0.5/ *_c);
104          scale(_h[ix],crossSection()/nanobarn/sumOfWeights());
105        }
106      }
107    }
108
109    /// @}
110
111
112    /// @name Histograms
113    /// @{
114    Histo1DPtr _h[7];
115    Histo1DPtr _r[7];
116    CounterPtr _c;
117    unsigned int _mode=0;
118    /// @}
119
120
121  };
122
123
124  RIVET_DECLARE_PLUGIN(BELLE_2005_I686014);
125
126}