1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/UnstableParticles.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief e+ e- > sigma+- sigmabar -+
 10  class BESIII_2020_I1814783 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2020_I1814783);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22      // Initialise and register projections
 23      declare(FinalState(), "FS");
 24      declare(UnstableParticles(), "UFS");
 25      book(_n_plus ,"/TMP/NPLUS" );
 26      book(_n_minus,"/TMP/NMINUS");
 27      if(isCompatibleWithSqrtS(2.396*GeV, 1E-2)) {
 28        book(_h_cTheta_A,3,1,1);
 29        book(_h_cTheta_B,3,1,2);
 30      }
 31    }
 32
 33    void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
 34      for(const Particle &child : p.children()) {
 35        if(child.children().empty()) {
 36          nRes[child.pid()]-=1;
 37          --ncount;
 38        }
 39        else
 40          findChildren(child,nRes,ncount);
 41      }
 42    }
 43
 44    /// Perform the per-event analysis
 45    void analyze(const Event& event) {
 46      const FinalState& fs = apply<FinalState>(event, "FS");
 47      // total hadronic and muonic cross sections
 48      map<long,int> nCount;
 49      int ntotal(0);
 50      for (const Particle& p : fs.particles()) {
 51        nCount[p.pid()] += 1;
 52        ++ntotal;
 53      }
 54      // find the Sigmas
 55      const FinalState& ufs = apply<UnstableParticles>(event, "UFS");
 56      for(unsigned int ix=0;ix<ufs.particles().size();++ix) {
 57        const Particle& p1 = ufs.particles()[ix];
 58        if(abs(p1.pid())!=3112&&abs(p1.pid())!=3222) continue;
 59        bool matched = false;
 60        // check fs
 61        bool fs = true;
 62        for(const Particle & child : p1.children()) {
 63          if(child.pid()==p1.pid()) {
 64            fs = false;
 65            break;
 66          }
 67        }
 68        if(!fs) continue;
 69        // find the children
 70        map<long,int> nRes = nCount;
 71        int ncount = ntotal;
 72        findChildren(p1,nRes,ncount);
 73        for(unsigned int iy=ix+1;iy<ufs.particles().size();++iy) {
 74          const Particle& p2 = ufs.particles()[iy];
 75          if(p2.pid() != -p1.pid()) continue;
 76          // check fs
 77          bool fs = true;
 78          for(const Particle & child : p2.children()) {
 79            if(child.pid()==p2.pid()) {
 80              fs = false;
 81              break;
 82            }
 83          }
 84          if(!fs) continue;
 85          map<long,int> nRes2 = nRes;
 86          int ncount2 = ncount;
 87          findChildren(p2,nRes2,ncount2);
 88          if(ncount2!=0) continue;
 89          matched=true;
 90          for(auto const & val : nRes2) {
 91            if(val.second!=0) {
 92              matched = false;
 93              break;
 94            }
 95          }
 96          if(matched) {
 97            if(abs(p1.pid())==3222) {
 98              _n_plus->fill();
 99              if(_h_cTheta_A) {
100                double cTheta = p1.pid()>0 ?
101                  cos(p1.momentum().polarAngle()) :
102                  cos(p2.momentum().polarAngle());
103                _h_cTheta_A->fill(cTheta);
104                _h_cTheta_B->fill(cTheta);
105              }
106            }
107            else if(abs(p1.pid())==3112)
108              _n_minus->fill();
109            break;
110          }
111        }
112        if(matched) break;
113      }
114
115    }
116
117
118    /// Normalise histograms etc., after the run
119    void finalize() {
120      if(_h_cTheta_A) {
121        normalize(_h_cTheta_A);
122        normalize(_h_cTheta_B);
123      }
124      double fact = crossSection()/ sumOfWeights() /picobarn;
125      for(unsigned int iy=1;iy<3;++iy) {
126        double sigma,error;
127        if(iy==1) {
128          sigma = _n_plus->val()*fact;
129          error = _n_plus->err()*fact;
130        }
131        else {
132          sigma = _n_minus->val()*fact;
133          error = _n_minus->err()*fact;
134        }
135        for(unsigned int ix=1;ix<3;++ix) {
136          Estimate1DPtr  mult;
137          book(mult, ix, 1, iy);
138          for (auto& b : mult->bins()) {
139            if (inRange(sqrtS()/GeV, b.xMin(), b.xMax())) {
140              b.set(sigma, error);
141            }
142          }
143        }
144      }
145    }
146    /// @}
147
148
149    /// @name Histograms
150    /// @{
151    CounterPtr _n_plus,_n_minus;
152    Histo1DPtr _h_cTheta_A,_h_cTheta_B;
153    /// @}
154
155  };
156
157
158  RIVET_DECLARE_PLUGIN(BESIII_2020_I1814783);
159
160}