1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief charm spectra at 10.55
  9  class CLEO_1988_I22954 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(CLEO_1988_I22954);
 14
 15
 16    /// @name Analysis methods
 17    /// @{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21      // projections
 22      declare(UnstableParticles(Cuts::abspid==411 ||
 23                                Cuts::abspid==421 ||
 24                                Cuts::abspid==413 ||
 25                                Cuts::abspid==423 ||
 26                                Cuts::abspid==431 ||
 27                                Cuts::abspid==4122), "UFS");
 28      // histos
 29      for (unsigned int ix=0; ix<4; ++ix) {
 30        book(_h_spect[ix],1+2*ix,1,1);
 31        book(_h_scale[ix],2+2*ix,1,1);
 32        book(_h_sigma[ix],14,1,1+ix);
 33        if (ix==3) continue;
 34        book(_h_spect[4+ix], 9+ix,1,1);
 35      }
 36      book(_h_scale[4],12,1,1);
 37      for (unsigned int ix=0;ix<2;++ix) {
 38	      for (unsigned int iy=0;iy<7;++iy) {
 39          book(_h_sigma_br[ix][iy],13,1+ix,1+iy);
 40        }
 41      }
 42    }
 43
 44
 45    /// Perform the per-event analysis
 46    void analyze(const Event& event) {
 47      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 48      for (const Particle& p : ufs.particles()) {
 49        double Emax = 0.5*sqrtS();
 50        double pmax = sqrt(sqr(Emax)-p.mass2());
 51        double xp = (p.E()+p.p3().mod())/(Emax+pmax);
 52        if (p.abspid()==413) {
 53          _h_spect[0]->fill(xp);
 54          _h_scale[0]->fill(xp);
 55          _h_spect[1]->fill(xp);
 56          _h_scale[1]->fill(xp);
 57          _h_sigma[2]->fill("10.55"s);
 58          for (unsigned int ix=0; ix<2; ++ix) {
 59            _h_sigma_br[ix][2]->fill(xp);
 60          }
 61          for (unsigned int ix=0; ix<2; ++ix) {
 62            _h_sigma_br[ix][3]->fill(xp);
 63          }
 64        }
 65        else if (p.abspid()==421) {
 66          _h_spect[2]->fill(xp);
 67          _h_scale[2]->fill(xp);
 68          _h_sigma[0]->fill("10.55"s);
 69          for (unsigned int ix=0; ix<2; ++ix) {
 70            _h_sigma_br[ix][0]->fill(xp);
 71          }
 72        }
 73        else if (p.abspid()==411) {
 74          _h_spect[3]->fill(xp);
 75          _h_scale[3]->fill(xp);
 76          _h_sigma[1]->fill("10.55"s);
 77          for (unsigned int ix=0; ix<2; ++ix) {
 78            _h_sigma_br[ix][1]->fill(xp);
 79          }
 80        }
 81        else if (p.abspid()==431) {
 82          _h_spect[4]->fill(xp);
 83          for (unsigned int ix=0;ix<2;++ix) {
 84            _h_sigma_br[ix][5]->fill(xp);
 85          }
 86        }
 87        else if (p.abspid()==4122) {
 88          _h_spect[5]->fill(xp);
 89          for (unsigned int ix=0;ix<2;++ix) {
 90            _h_sigma_br[ix][6]->fill(xp);
 91          }
 92        }
 93        else if (p.abspid()==423) {
 94          _h_spect[6]->fill(xp);
 95          _h_scale[4]->fill(xp);
 96          _h_sigma[3]->fill("10.55"s);
 97          for (unsigned int ix=0;ix<2;++ix) {
 98            _h_sigma_br[ix][4]->fill(xp);
 99          }
100        }
101      }
102    }
103
104
105    /// Normalise histograms etc., after the run
106    void finalize() {
107      const vector<double> br = {0.677*0.03947,0.677*0.0822,0.03947,0.0938,0.045,0.0628,0.03947};
108      const double fact1 = crossSection()/picobarn/sumOfWeights();
109      const double fact2 = crossSection()/nanobarn/sumOfWeights()*sqr(sqrtS());
110      for (unsigned int ix=0; ix<7; ++ix) {
111        scale(_h_spect[ix],fact1*br[ix]);
112        if (ix>4) continue;
113        scale(_h_scale[ix],fact2);
114        if (ix==4) continue;
115        scale(_h_sigma[ix],fact1*1e-3);
116      }
117      const double width[2]={1.,0.5};
118      for (unsigned int ix=0; ix<2; ++ix) {
119        scale(_h_sigma_br[ix][0],width[ix]*fact1*br[2]);
120        scale(_h_sigma_br[ix][1],width[ix]*fact1*br[3]);
121        scale(_h_sigma_br[ix][2],width[ix]*fact1*br[0]);
122        scale(_h_sigma_br[ix][3],width[ix]*fact1*br[1]);
123        scale(_h_sigma_br[ix][4],width[ix]*fact1*br[6]);
124        scale(_h_sigma_br[ix][5],width[ix]*fact1*br[4]);
125        scale(_h_sigma_br[ix][6],width[ix]*fact1*br[5]);
126      }
127    }
128
129    /// @}
130
131
132    /// @name Histograms
133    /// @{
134    Histo1DPtr _h_spect[7],_h_scale[5],_h_sigma_br[2][7];
135    BinnedHistoPtr<string> _h_sigma[4];
136    /// @}
137
138  };
139
140
141  RIVET_DECLARE_PLUGIN(CLEO_1988_I22954);
142
143}