1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief Add a short analysis description here
  9  class GAMMAGAMMA_1981_I158474 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(GAMMAGAMMA_1981_I158474);
 14
 15
 16    /// @name Analysis methods
 17    /// @{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21      // Initialise and register projections
 22      declare(FinalState(), "FS");
 23      book(_n3pi, "TMP/n3pi");
 24      book(_n4pi, "TMP/n4pi");
 25      book(_n5pi, "TMP/n5pi");
 26      book(_n6pi, "TMP/n6pi");
 27      book(_n35pi, "TMP/n35pi");
 28      book(_n46pi, "TMP/n46pi");
 29      book(_nC2, "TMP/nC2");
 30      book(_nC4, "TMP/nC4");
 31      book(_nmu, "TMP/nmu");
 32    }
 33
 34
 35    /// Perform the per-event analysis
 36    void analyze(const Event& event) {
 37      const FinalState& fs = apply<FinalState>(event, "FS");
 38
 39      map<long,int> nCount;
 40      int ntotal(0);
 41      for (const Particle& p : fs.particles()) {
 42	nCount[p.pid()] += 1;
 43	++ntotal;
 44      }
 45      // mu+mu- + photons
 46      if(nCount[-13]==1 and nCount[13]==1 &&
 47	 ntotal==2+nCount[22])
 48	_nmu->fill();
 49      else {
 50	if(ntotal==3 && nCount[211] == 1 && nCount[-211]==1 && nCount[111]==1 ) {
 51	  _n3pi->fill();
 52	}
 53	if(ntotal==4 && nCount[211] == 1 && nCount[-211]==1 && nCount[111]==2 ) {
 54	  _n4pi->fill();
 55	}
 56	if(ntotal==5 && nCount[211] == 2 && nCount[-211]==2 && nCount[111]==1 ) {
 57	  _n5pi->fill();
 58	}
 59	if(ntotal==6 && nCount[211] == 2 && nCount[-211]==2 && nCount[111]==2 ) {
 60	  _n6pi->fill();
 61	}
 62	if(nCount[211] == 1 && nCount[-211]==1 && ntotal == 2+nCount[111]) {
 63	  _nC2->fill();
 64	}
 65	if(nCount[211] == 2 && nCount[-211]==2 && ntotal == 4+nCount[111]) {
 66	  _nC4->fill();
 67	}
 68	if((nCount[211]+nCount[-211]+nCount[111])==ntotal ) {
 69	  if(ntotal==3 || ntotal ==5)
 70	    _n35pi->fill();
 71	  else if(ntotal==4 || ntotal==6)
 72	    _n46pi ->fill();
 73	}
 74      }
 75    }
 76
 77
 78    /// Normalise histograms etc., after the run
 79    void finalize() {
 80      double fact = crossSection()/ sumOfWeights() /nanobarn;
 81      for(unsigned int ix=1;ix<7;++ix) {
 82        double sigma = 0.0, error = 0.0;
 83        if(ix==1) {
 84          sigma = _n3pi->val()*fact;
 85          error = _n3pi->err()*fact;
 86        }
 87        else if(ix==2) {
 88          sigma = _n4pi->val()*fact;
 89          error = _n4pi->err()*fact;
 90        }
 91        else if(ix==3) {
 92          sigma = _n5pi->val()*fact;
 93          error = _n5pi->err()*fact;
 94        }
 95        else if(ix==4) {
 96          sigma = _n6pi->val()*fact;
 97          error = _n6pi->err()*fact;
 98        }
 99        else if(ix==5) {
100          sigma = _n35pi->val()*fact;
101          error = _n35pi->err()*fact;
102        }
103        else if(ix==6) {
104          sigma = _n46pi->val()*fact;
105          error = _n46pi->err()*fact;
106        }
107        Estimate1DPtr mult;
108        book(mult, 1, 1, ix);
109        for (auto& b : mult->bins()) {
110          if (inRange(sqrtS()/GeV, b.xMin(), b.xMax())) {
111            b.set(sigma, error);
112          }
113        }
114      }
115      for (unsigned int ix=1;ix<3;++ix) {
116        Estimate0D R = (ix==1? *_nC2 : *_nC4)/ *_nmu;
117        double sig_h = (ix ==1 ? _nC2 : _nC4)->val()*fact;
118        double err_h = (ix ==1 ? _nC2 : _nC4)->err()*fact;
119        double sig_m = _nmu->val()*fact;
120        double err_m = _nmu->err()*fact;
121        std::ostringstream title;
122        if(ix==1)
123          title << "sigma_2pi";
124        else
125          title << "sigma_4pi";
126        Estimate1DPtr hadrons;
127        book(hadrons, title.str());
128        Estimate1DPtr muons;
129        book(muons, "sigma_muons");
130        Estimate1DPtr mult;
131        book(mult, 2,1,ix);
132        for (auto& b : mult->bins()) {
133          if (inRange(sqrtS()/GeV, b.xMin(), b.xMax())) {
134            b.set(R.val(), R.errPos());
135            hadrons->bin(b.index()).set(sig_h, err_h);
136            if(ix==1) muons->bin(b.index()).set(sig_m, err_m);
137          }
138        }
139      }
140    }
141
142    /// @}
143
144
145    /// @name Histograms
146    /// @{
147    CounterPtr _n3pi,_n4pi,_n5pi,_n6pi,_n35pi,_n46pi,_nC2,_nC4,_nmu;
148    /// @}
149
150
151  };
152
153
154  RIVET_DECLARE_PLUGIN(GAMMAGAMMA_1981_I158474);
155
156
157}