1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/ChargedFinalState.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief BELLE J/psi production
 10  class BELLE_2009_I811289 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2009_I811289);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22      declare("FS",FinalState());
 23      declare("CFS",FinalState());
 24      // book histos
 25      for (unsigned int ix=0; ix<5; ++ix) {
 26        book(_h_spect[ix], 2, 1, 1+ix);
 27      }
 28      for (unsigned int ix=0; ix<2; ++ix) {
 29        book(_h_sigma[ix], 1, 1, 1+ix);
 30        for (unsigned int iy=0; iy<3; ++iy) {
 31          book(_h_angle[ix][iy], 3, 1+ix, 1+iy);
 32        }
 33      }
 34    }
 35
 36    void plotHelicityAngle(const Particle& p, Histo1DPtr h) {
 37      if (p.children().size()!=2) return;
 38      if (p.children()[0].abspid()!=PID::MUON || p.children()[0].pid()!=-p.children()[1].pid()) return;
 39      Vector3 axis = p.p3().unit();
 40      const LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(p.mom().betaVec());
 41      for (const Particle& child : p.children()) {
 42        if (child.pid()!=PID::MUON) continue;
 43        h->fill(axis.dot(boost.transform(child.mom()).p3().unit()));
 44      }
 45    }
 46
 47    /// Perform the per-event analysis
 48    void analyze(const Event& event) {
 49      // 4 charged particle veto
 50      if (apply<FinalState>(event,"CFS").particles().size()<4) vetoEvent;
 51      // find the charmonium
 52      Particles D,Jpsi_direct,Jpsi_feed,onium;
 53      for (const Particle & p : apply<FinalState>(event,"FS").particles()) {
 54        Particle parent = p.parents()[0], Jpsi;
 55        while (!parent.parents().empty()) {
 56          // charm mesons
 57          if (parent.abspid()==411 || parent.abspid()==421 || parent.abspid()==431) {
 58            break;
 59          }
 60          // c cbar mesons
 61          if ((parent.abspid()%1000)/10 == 44 ) {
 62            if ((parent.parents()[0].abspid()%1000)/10!=44) break;
 63            if (parent.pid()==443) Jpsi=parent;
 64          }
 65          parent = parent.parents()[0];
 66        }
 67        // charm
 68        if (parent.abspid()==411 || parent.abspid()==421 || parent.abspid()==431) {
 69          bool found=false;
 70          for (auto & p : D) {
 71            // D already in list
 72            if (fuzzyEquals(p.mom(),parent.mom())) {
 73              found=true;
 74              break;
 75            }
 76          }
 77          if (!found) D.push_back(parent);
 78        }
 79        // c cbar
 80        else if ((parent.abspid()%1000)/10 == 44) {
 81          if (parent.pid()==443) {
 82            bool found=false;
 83            for (const auto& p : Jpsi_direct) {
 84              // psi already in list
 85              if (fuzzyEquals(p.mom(),parent.mom())) {
 86                found=true;
 87                break;
 88              }
 89            }
 90            if (!found) Jpsi_direct.push_back(parent);
 91          }
 92          else {
 93            bool found=false;
 94            for (const auto & p : onium) {
 95              // psi already in list
 96              if (fuzzyEquals(p.mom(),parent.mom())) {
 97                found=true;
 98                break;
 99              }
100            }
101            if (!found) onium.push_back(parent);
102          }
103        }
104        if (Jpsi.pid()==443) {
105          bool found=false;
106          for (const auto& p : Jpsi_feed) {
107            // psi already in list
108            if (fuzzyEquals(p.mom(),Jpsi.mom())) {
109              found=true;
110              break;
111            }
112          }
113          if (!found) Jpsi_feed.push_back(Jpsi);
114        }
115      }
116      if (Jpsi_feed.empty() && Jpsi_direct.empty()) vetoEvent;
117      // cross sections first
118      if (D.empty() && onium.empty() && Jpsi_direct.size()==1) {
119        _h_sigma[1]->fill("10.6"s);
120        for (const Particle & p : Jpsi_direct) {
121          double modp=p.p3().mod();
122          _h_spect[4]->fill(modp);
123          plotHelicityAngle(p,_h_angle[0][2]);
124          _h_angle[1][2]->fill(abs(p.p3().z()/modp));
125        }
126      }
127      else {
128        _h_sigma[0]->fill("10.6"s);
129        for (const Particle& p : Jpsi_direct) {
130          double modp=p.p3().mod();
131          _h_spect[3]->fill(modp);
132          plotHelicityAngle(p,_h_angle[0][1]);
133          _h_angle[1][1]->fill(abs(p.p3().z()/modp));
134        }
135        for (const Particle& p : Jpsi_feed) {
136          double modp=p.p3().mod();
137          _h_spect[3]->fill(modp);
138          plotHelicityAngle(p,_h_angle[0][1]);
139          _h_angle[1][1]->fill(abs(p.p3().z()/modp));
140        }
141        if (!D.empty()) {
142          for (const Particle& p : Jpsi_direct) {
143            _h_spect[1]->fill(p.p3().mod());
144          }
145          for (const Particle& p : Jpsi_feed) {
146            _h_spect[1]->fill(p.p3().mod());
147          }
148        }
149        else {
150          for (const Particle& p : Jpsi_direct) {
151            _h_spect[2]->fill(p.p3().mod());
152          }
153          for (const Particle& p : Jpsi_feed) {
154            _h_spect[2]->fill(p.p3().mod());
155          }
156        }
157      }
158      for (const Particle& p : Jpsi_direct) {
159        double modp=p.p3().mod();
160        _h_spect[0]->fill(modp);
161        plotHelicityAngle(p,_h_angle[0][0]);
162        _h_angle[1][0]->fill(abs(p.p3().z()/modp));
163      }
164      for (const Particle& p : Jpsi_feed) {
165        double modp=p.p3().mod();
166        _h_spect[0]->fill(modp);
167        plotHelicityAngle(p,_h_angle[0][0]);
168        _h_angle[1][0]->fill(abs(p.p3().z()/modp));
169      }
170    }
171
172
173    /// Normalise histograms etc., after the run
174    void finalize() {
175      scale(_h_spect, crossSection()/sumOfWeights()/femtobarn);
176      for (unsigned int ix=0; ix<2; ++ix) {
177        scale(_h_sigma[ix], crossSection()/sumOfWeights()/picobarn);
178        normalize(_h_angle[ix], 1.0, false);
179      }
180    }
181
182    /// @}
183
184
185    /// @name Histograms
186    /// @{
187    BinnedHistoPtr<string> _h_sigma[2];
188    Histo1DPtr _h_spect[5],_h_angle[2][3];
189    /// @}
190
191
192  };
193
194
195  RIVET_DECLARE_PLUGIN(BELLE_2009_I811289);
196
197}