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 D0 -> omega phi
 10  class BESIII_2022_I1900094 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2022_I1900094);
 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(UnstableParticles(Cuts::abspid==PID::D0), "UFS");
 24      for(unsigned int ix=0;ix<2;++ix)
 25	book(_h[ix],1,1,1+ix);
 26    }
 27
 28    void findDecayProducts(const Particle & mother, unsigned int & nstable,
 29			   Particles & pip , Particles & pim , Particles & pi0) {
 30      for(const Particle & p : mother.children()) {
 31        int id = p.pid();
 32        if ( id == PID::KPLUS || id == PID::KMINUS || id == PID::K0S || id == PID::K0L ) {
 33	  ++nstable;
 34	}
 35	else if (id == PID::PIPLUS) {
 36	  pip.push_back(p);
 37	  ++nstable;
 38	}
 39	else if (id == PID::PIMINUS) {
 40	  pim.push_back(p);
 41	  ++nstable;
 42	}
 43	else if (id == PID::PI0) {
 44	  pi0.push_back(p);
 45	  ++nstable;
 46	}
 47	else if ( !p.children().empty() ) {
 48	  findDecayProducts(p, nstable, pip, pim, pi0);
 49	}
 50	else
 51	  ++nstable;
 52      }
 53    }
 54
 55    /// Perform the per-event analysis
 56    void analyze(const Event& event) {
 57      for(const Particle & D0 : apply<UnstableParticles>(event, "UFS").particles()) {
 58	Particle omega,phi;
 59	if(D0.children()[0].pid()==PID::OMEGA && D0.children()[1].pid()==PID::PHI) {
 60	  omega = D0.children()[0];
 61	  phi   = D0.children()[1];
 62	}
 63	else if(D0.children()[1].pid()==PID::OMEGA && D0.children()[0].pid()==PID::PHI) {
 64	  omega = D0.children()[1];
 65	  phi   = D0.children()[0];
 66	}
 67	else
 68	  continue;
 69	LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(D0.momentum().betaVec());
 70	// first the phi
 71	if(phi.children().size()==2 &&
 72	   phi.children()[0].abspid()==PID::KPLUS &&
 73	   phi.children()[1].abspid()==PID::KPLUS ) {
 74	  // first boost all relevant momenta to D0 rest frame
 75	  FourMomentum pD0  = boost.transform(D0.momentum());
 76	  FourMomentum pphi = boost.transform(phi.momentum());
 77	  FourMomentum pK   = boost.transform(phi.children()[0].momentum());
 78	  LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(pphi.betaVec());
 79	  Vector3 axis1 = boost2.transform(pD0).p3().unit();
 80	  Vector3 axis2 = boost2.transform(pK ).p3().unit();
 81	  _h[1]->fill(abs(axis1.dot(axis2)));
 82	}
 83	// then the omega
 84	unsigned int nstable(0);
 85	Particles pip, pim, pi0;
 86	findDecayProducts(omega, nstable, pip, pim, pi0);
 87	if(nstable==3 && pip.size()==1 && pip.size()==1 && pi0.size()==1) {
 88	  FourMomentum pD0    = boost.transform(D0.momentum());
 89	  FourMomentum pomega = boost.transform(omega.momentum() );
 90	  FourMomentum ppip   = boost.transform(pip[0].momentum());
 91	  FourMomentum ppim   = boost.transform(pim[0].momentum());
 92	  LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(pomega.betaVec());
 93	  Vector3 axis = boost2.transform(pD0).p3().unit();
 94	  Vector3 pp   = boost2.transform(ppip).p3();
 95	  Vector3 pm   = boost2.transform(ppim).p3();
 96	  Vector3 axis2 = pp.cross(pm).unit();
 97	  _h[0]->fill(abs(axis.dot(axis2)));
 98	}
 99      }
100    }
101
102
103    /// Normalise histograms etc., after the run
104    void finalize() {
105      for(unsigned int ix=0;ix<2;++ix) normalize(_h[ix]);
106    }
107
108    /// @}
109
110
111    /// @name Histograms
112    /// @{
113    Histo1DPtr _h[2];
114    /// @}
115
116
117  };
118
119
120  RIVET_DECLARE_PLUGIN(BESIII_2022_I1900094);
121
122}