1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief chi_cJ -> phi phi
  9  class BESIII_2023_I2627838 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2023_I2627838);
 14
 15
 16    /// @name Analysis methods
 17    /// @{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21      declare(UnstableParticles(Cuts::abspid==100443), "UFS");
 22      // counters
 23      for (unsigned int ix=0; ix<3; ++ix) {
 24        book(_n[ix],"TMP/n_"+toString(ix+1));
 25        for (unsigned int iy=0; iy<ix+1; ++iy)
 26          book(_m[ix][iy],"TMP/m_"+toString(ix+1)+"_"+toString(iy+1));
 27      }
 28    }
 29
 30
 31    /// Perform the per-event analysis
 32    void analyze(const Event& event) {
 33      for (const Particle& psi : apply<UnstableParticles>(event, "UFS").particles()) {
 34        if (psi.children().size()!=2) continue;
 35        // find chi_c gamma decay and type of chi_c
 36        Particle chi;
 37        if (psi.children()[0].pid()==PID::GAMMA) {
 38          if (psi.children()[1].pid()==10441 ||
 39              psi.children()[1].pid()==20443 ||
 40              psi.children()[1].pid()==445 ) {
 41            chi = psi.children()[1];
 42          }
 43        }
 44        else if (psi.children()[1].pid()==PID::GAMMA) {
 45          if (psi.children()[0].pid()==10441 ||
 46              psi.children()[0].pid()==20443 ||
 47              psi.children()[0].pid()==445) {
 48	          chi = psi.children()[0];
 49          }
 50        }
 51        else {
 52          continue;
 53        }
 54        unsigned int iloc=0;
 55        if      (chi.pid()==10441) iloc=0;
 56        else if (chi.pid()==20443) iloc=1;
 57        else if (chi.pid()==  445) iloc=2;
 58        else {
 59          continue;
 60        }
 61        // require chi_c -> phi phi
 62        if (chi.children().size()!=2) continue;
 63        if (chi.children()[0].pid()!=333 ||
 64            chi.children()[0].pid()!=333) continue;
 65        bool found = true;
 66        Particle Km[2],Kp[2];
 67        for (unsigned int ix=0; ix<2; ++ix) {
 68          // required K+K- decay
 69          if (chi.children()[ix].children().size()!=2) {
 70            found = false;
 71            break;
 72          }
 73          if (chi.children()[ix].children()[0].pid()!=-chi.children()[ix].children()[1].pid() ||
 74              chi.children()[ix].children()[0].abspid()!=321) {
 75            found = false;
 76            break;
 77          }
 78          if (chi.children()[ix].children()[0].pid()>0) {
 79            Kp[ix] = chi.children()[ix].children()[0];
 80            Km[ix] = chi.children()[ix].children()[1];
 81          }
 82          else {
 83            Kp[ix] = chi.children()[ix].children()[1];
 84            Km[ix] = chi.children()[ix].children()[0];
 85          }
 86        }
 87        if (!found) continue;
 88        // fill count
 89        _n[iloc]->fill();
 90        // boost to psi(2S) frame
 91        LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(psi.mom().betaVec());
 92        LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(chi.mom().betaVec());
 93        Vector3 axis1 = boost1.transform(chi.mom()).p3().unit();
 94        double cTheta[3]={0.,0.,0.};
 95        for (unsigned int ix=0; ix<2; ++ix) {
 96          FourMomentum pPhi = boost2.transform(boost1.transform(chi.children()[ix].mom()));
 97          Vector3 axis2 = pPhi.p3().unit();
 98          if (ix==0) cTheta[0] = axis1.dot(axis2);
 99          LorentzTransform boost3 = LorentzTransform::mkFrameTransformFromBeta(pPhi.betaVec());
100          FourMomentum pK = boost3.transform(boost2.transform(boost1.transform(Kp[ix].mom())));
101          cTheta[ix+1] = axis2.dot(pK.p3().unit());
102        }
103        if (iloc==0) {
104          _m[0][0]->fill(0.25*(3. - 5.*sqr(cTheta[1])));
105        }
106        else if (iloc==1) {
107          _m[1][0]->fill(0.625*(1. + sqr(cTheta[1]) - 4.*sqr(cTheta[0])));
108          _m[1][1]->fill(     -(3.                  -10.*sqr(cTheta[0])));
109        }
110        else {
111          _m[2][0]->fill(-0.125*(25.*sqr(cTheta[1])*sqr(cTheta[2]) - 10*sqr(cTheta[1]) - 10*sqr(cTheta[2]) + 3.));
112          _m[2][1]->fill( 0.05 *(8 - 140.*sqr(cTheta[1]) + 325.*sqr(cTheta[0])*sqr(cTheta[2])
113                         + 250.*sqr(cTheta[0])*sqr(cTheta[1])*sqr(cTheta[2])));
114          _m[2][2]->fill( 0.75 *(1. + 5.*sqr(cTheta[1]) - 25.* sqr(cTheta[0])*sqr(cTheta[2])));
115        }
116      }
117    }
118
119
120    /// Normalise histograms etc., after the run
121    void finalize() {
122      for (unsigned int ix=0; ix<3; ++ix) {
123        if (_n[ix]->numEntries()==0)  continue;
124        scale(_m[ix], 1.0/ *_n[ix]);
125        if (ix==0) {
126          double x = _m[0][0]->val()/(1.-2.*_m[0][0]->val());
127          pair<double,double> dx = make_pair(x-(-_m[0][0]->err() + _m[0][0]->val())/(1 + 2*_m[0][0]->err() - 2*_m[0][0]->val()),
128                                            (_m[0][0]->err() + _m[0][0]->val())/(1 - 2*_m[0][0]->err() - 2*_m[0][0]->val())-x);
129          double rx = sqrt(abs(x));
130          dx.first  *= 0.5/rx;
131          dx.second *= 0.5/rx;
132          if (x<0.) rx*=-1.;
133          Estimate0DPtr h_x;
134          book(h_x, 1, 1, 1);
135          h_x->set(rx, dx);
136        }
137        else if (ix==1) {
138          double u1 = -4.*_m[1][0]->val()/(-1. + 4.*_m[1][0]->val() + _m[1][1]->val());
139          double u2 = -4.*_m[1][1]->val()/(-1. + 4.*_m[1][0]->val() + _m[1][1]->val());
140          double O1 = _m[1][0]->val(), DO1 = _m[1][0]->err();
141          double O2 = _m[1][1]->val(), DO2 = _m[1][1]->err();
142          double root1 = sqrt(sqr(-1 + 4*O1 + O2)*(sqr(DO2)*sqr(1 - 4*O1) + 16*sqr(DO1)*sqr(O2)));
143          pair<double,double> du1 = make_pair(-4*DO1*DO2/(4*DO1*DO2*(-1 + 4*O1 + O2) - root1),
144                                               4*DO1*DO2/(4*DO1*DO2*(-1 + 4*O1 + O2) + root1));
145          double ru1 = sqrt(abs(u1));
146          du1.first  *=0.5/ru1;
147          du1.second *=0.5/ru1;
148          if (u1<0.) ru1*=-1.;
149          double root2 = sqrt((sqr(DO2)*sqr(O1) + sqr(DO1)*sqr(-1 + O2))*sqr(-1 + 4*O1 + O2));
150          pair<double,double> du2 = make_pair(-4*DO1*DO2/(DO1*DO2*(-1 + 4*O1 + O2) - root2),
151                                               4*DO1*DO2/(DO1*DO2*(-1 + 4*O1 + O2) + root2));
152          double ru2 = sqrt(abs(u2));
153          du2.first  *=0.5/ru2;
154          du2.second *=0.5/ru2;
155          if (u2<0.) ru2*=-1.;
156          Estimate0DPtr h_u1;
157          book(h_u1, 1, 2, 1);
158          h_u1->set(ru1, du1);
159          Estimate0DPtr h_u2;
160          book(h_u2, 1, 2, 2);
161          h_u2->set(ru2, du2);
162        }
163        else {
164          double O1 = _m[2][0]->val(), DO1 = _m[2][0]->err();
165          double O2 = _m[2][1]->val(), DO2 = _m[2][1]->err();
166          double O4 = _m[2][2]->val(), DO4 = _m[2][2]->err();
167          double w1 = -O1/(-1 + 4*O1 + 2*O2 + 2*O4);
168          double w2 = -O2/(-1 + 4*O1 + 2*O2 + 2*O4);
169          double w4 = -O4/(-1 + 4*O1 + 2*O2 + 2*O4);
170          double root1 = sqrt((16*sqr(DO1)*sqr(DO4)*sqr(O2) + sqr(DO2)*(sqr(DO4) *
171                         sqr(1 - 4*O1) + 16*sqr(DO1)*sqr(O4)))/sqr(-1 + 4*O1 + 2*O2 + 2*O4));
172          pair<double,double> dw1 = make_pair( DO1*DO2*DO4*(4*DO1*DO2*DO4 + (-1 + 4*O1 + 2*O2 + 2*O4)*root1) /
173                                               ((-1 + 4*O1 + 2*O2 + 2*O4)*(16*sqr(DO1)*sqr(DO4)*sqr(O2) +
174                                               sqr(DO2)*(sqr(DO4)*(-16*sqr(DO1) + sqr(1 - 4*O1)) + 16*sqr(DO1)*sqr(O4)))),
175                                               DO1*DO2*DO4*(-4*DO1*DO2*DO4 + (-1 + 4*O1 + 2*O2 + 2*O4)*root1) /
176                                               ((-1 + 4*O1 + 2*O2 + 2*O4)*(16*sqr(DO1)*sqr(DO4)*sqr(O2) +
177                                               sqr(DO2)*(sqr(DO4)*(-16*sqr(DO1) +
178                                               sqr(1 - 4*O1)) + 16*sqr(DO1)*sqr(O4)))));
179          double rw1 = sqrt(abs(w1));
180          dw1.first  *=0.5/rw1;
181          dw1.second *=0.5/rw1;
182          if (w1<0.) rw1*=-1.;
183          double root2 = sqrt((4*sqr(DO2)*sqr(DO4)*sqr(O1) + sqr(DO1)*(sqr(DO4) *
184                         sqr(1 - 2*O2) + 4*sqr(DO2)*sqr(O4)))/sqr(-1 + 4*O1 + 2*O2 + 2*O4));
185          pair<double,double> dw2 = make_pair(DO1*DO2*DO4*(2*DO1*DO2*DO4 + (-1 + 4*O1 + 2*O2 + 2*O4)*root2) /
186                                              ((-1 + 4*O1 + 2*O2 + 2*O4)*(4*sqr(DO2)*sqr(DO4)*sqr(O1) +
187                                              sqr(DO1)*(sqr(DO4)*(-4*sqr(DO2) + sqr(1 - 2*O2)) + 4*sqr(DO2)*sqr(O4)))),
188                                              DO1*DO2*DO4*(-2*DO1*DO2*DO4 + (-1 + 4*O1 + 2*O2 + 2*O4)*root2) /
189                                              ((-1 + 4*O1 + 2*O2 + 2*O4)*(4*sqr(DO2)*sqr(DO4)*sqr(O1) +
190                                              sqr(DO1)*(sqr(DO4)*(-4*sqr(DO2) +
191                                              sqr(1 - 2*O2)) + 4*sqr(DO2)*sqr(O4)))));
192          double rw2 = sqrt(abs(w2));
193          dw2.first  *=0.5/rw2;
194          dw2.second *=0.5/rw2;
195          if (w2<0.) rw2*=-1.;
196          double root3 = sqrt((4*sqr(DO2)*sqr(DO4)*sqr(O1) + sqr(DO1)*(4*sqr(DO4)*sqr(O2) +
197                         sqr(DO2)*sqr(1 - 2*O4)))/sqr(-1 + 4*O1 + 2*O2 + 2*O4));
198          pair<double,double> dw4 = make_pair(DO1*DO2*DO4*(2*DO1*DO2*DO4 + (-1 + 4*O1 + 2*O2 + 2*O4)*root3) /
199                                              ((4*sqr(DO2)*sqr(DO4)*sqr(O1) + sqr(DO1)*(4*sqr(DO4)*sqr(O2) +
200                                              sqr(DO2)*(-4*sqr(DO4) + sqr(1 - 2*O4)))) * (-1 + 4*O1 + 2*O2 + 2*O4)),
201                                              DO1*DO2*DO4*(-2*DO1*DO2*DO4 + (-1 + 4*O1 + 2*O2 + 2*O4)*root3) /
202                                              ((4*sqr(DO2)*sqr(DO4)*sqr(O1) + sqr(DO1)*(4*sqr(DO4)*sqr(O2) +
203                                              sqr(DO2)*(-4*sqr(DO4) + sqr(1 - 2*O4)))) * (-1 + 4*O1 + 2*O2 + 2*O4)));
204          double rw4 = sqrt(abs(w4));
205          dw4.first  *=0.5/rw4;
206          dw4.second *=0.5/rw4;
207          if (w4<0.) rw4*=-1.;
208          Estimate0DPtr h_w1;
209          book(h_w1, 1, 3, 1);
210          h_w1->set(rw1, dw1);
211          Estimate0DPtr h_w2;
212          book(h_w2, 1, 3, 2);
213          h_w2->set(rw2, dw2);
214          Estimate0DPtr h_w4;
215          book(h_w4, 1, 3, 3);
216          h_w4->set(rw4, dw4);
217        }
218      }
219    }
220
221    /// @}
222
223
224    /// @name Histograms
225    /// @{
226    CounterPtr _n[3],_m[3][3];
227    /// @}
228
229
230  };
231
232
233  RIVET_DECLARE_PLUGIN(BESIII_2023_I2627838);
234
235}