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- -> double charmonium
 10  class BELLE_2004_I653673 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2004_I653673);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22      // projections
 23      declare("FS",FinalState());
 24      declare("UFS",UnstableParticles(Cuts::pid==443 || Cuts::pid==100443 ||
 25                    Cuts::pid==441 || Cuts::pid==10441 || Cuts::pid==100441));
 26      // histograms
 27      for (unsigned int ix=0; ix<3; ++ix) {
 28        for (unsigned int iy=0;iy<3;++iy) {
 29          book(_p[ix][iy],"TMP/p_"+toString(ix+1)+"_"+toString(iy+1),
 30               refData<YODA::BinnedEstimate<string>>(3,1+ix,1+iy));
 31          if (ix==2) continue;
 32          book(_h_sigma[ix][iy], 1+ix, 1   , 1+2*iy);
 33          book(_h_angle[ix][iy], 4   , 1+ix, 1+iy  );
 34        }
 35      }
 36    }
 37
 38    void findChildren(const Particle& p,map<long,int>& nRes, int& ncount,
 39                      unsigned int & nCharged) {
 40      for (const Particle &child : p.children()) {
 41        if (child.children().empty()) {
 42          --nRes[child.pid()];
 43          --ncount;
 44          if (PID::isCharged(p.pid())) ++nCharged;
 45        }
 46        else {
 47          findChildren(child,nRes,ncount,nCharged);
 48        }
 49      }
 50    }
 51
 52    double helicityAngle(const Particle & p) const {
 53      if (p.children().size()!=2) return 10.;
 54      if (p.children()[0].abspid()!=PID::MUON || p.children()[0].pid()!=-p.children()[1].pid()) {
 55        return 10.;
 56      }
 57      Vector3 axis = p.p3().unit();
 58      const LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
 59      for (const Particle& child : p.children()) {
 60        if (child.pid()!=PID::MUON) continue;
 61        return axis.dot(boost.transform(child.momentum()).p3().unit());
 62      }
 63      return 10;
 64    }
 65
 66
 67    /// Perform the per-event analysis
 68    void analyze(const Event& event) {
 69      // final state particles
 70      const FinalState& fs = apply<FinalState>(event, "FS");
 71      map<long,int> nCount;
 72      int ntotal(0);
 73      for (const Particle& p : fs.particles()) {
 74      	nCount[p.pid()] += 1;
 75      	++ntotal;
 76      }
 77      // loop over J/psi and psi(2S)
 78      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 79      bool matched=false;
 80      for (const Particle& p : ufs.particles(Cuts::pid==443 || Cuts::pid==100443 )) {
 81      	if (p.children().empty()) continue;
 82      	map<long,int> nRes = nCount;
 83      	int ncount = ntotal;
 84        unsigned int nCharged=0;
 85      	findChildren(p,nRes,ncount,nCharged);
 86        // eta_c, chi_c0, eta_c(2S)
 87      	for (const Particle& p2 : ufs.particles(Cuts::pid==441 || Cuts::pid==10441 || Cuts::pid==100441)) {
 88          map<long,int> nResB = nRes;
 89          int ncountB = ncount;
 90          unsigned int nChargedB=0;
 91          findChildren(p2,nResB,ncountB,nChargedB);
 92          if (ncountB!=0) continue;
 93          matched = true;
 94          for (const auto& val : nResB) {
 95            if (val.second!=0) {
 96              matched = false;
 97              break;
 98            }
 99          }
100          if (matched) {
101            unsigned int ipsi = p.pid()/100000;
102            unsigned int ieta = p2.pid()/10000;
103            if (ieta>1) ieta=2;
104            // fill the cross sections
105            if ((ipsi==0 && nChargedB>2) || (ipsi==1 && nChargedB>0)) {
106              _h_sigma[ipsi][ieta]->fill(_ecms);
107            }
108            if (ipsi>0) break;
109            // angular dists for J/psi only
110            // production
111            const double cProd = p.p3().z()/p.p3().mod();
112            _h_angle[0][ieta]->fill(abs(cProd));
113            _p[0][ieta]->fill(_ecms,-1.25*(1.-3.*sqr(cProd)));
114            _p[2][ieta]->fill(_ecms,-1.25*(1.-3.*sqr(cProd)));
115            // helicity angle
116            const double cHel = helicityAngle(p);
117            if (cHel>1.) break;
118            _h_angle[1][ieta]->fill(abs(cHel));
119            _p[1][ieta]->fill(_ecms,-1.25*(1.-3.*sqr(cHel)));
120            _p[2][ieta]->fill(_ecms,-1.25*(1.-3.*sqr(cHel)));
121            break;
122          }
123        }
124        if (matched) break;
125      }
126    }
127
128
129    /// Normalise histograms etc., after the run
130    void finalize() {
131      for (unsigned int ix=0; ix<2; ++ix) {
132        scale(_h_sigma[ix], crossSection()/sumOfWeights()/femtobarn);
133        normalize(_h_angle[ix]);
134      }
135      // extract the alpha parameters
136      for (unsigned int ix=0; ix<3; ++ix) {
137        for (unsigned int iy=0; iy<3; ++iy) {
138          const double val = _p[ix][iy]->bin(1).mean(2);
139          const double err = _p[ix][iy]->bin(1).stdErr(2);
140          BinnedEstimatePtr<string> tmp;
141          book(tmp,3,1+ix,1+iy);
142          const double alpha = 3.*val/(1-val);
143          const double error = 3./sqr(1.-val)*err;
144          tmp->bin(1).set(alpha,error);
145        }
146      }
147    }
148
149    /// @}
150
151
152    /// @name Histograms
153    /// @{
154    BinnedHistoPtr<string> _h_sigma[2][3];
155    Histo1DPtr _h_angle[2][3];
156    BinnedProfilePtr<string> _p[3][3];
157    string _ecms="10.6";
158    /// @}
159
160
161  };
162
163
164  RIVET_DECLARE_PLUGIN(BELLE_2004_I653673);
165
166}