rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

BABAR_2012_I1122034

$X(3915)$ production in $\gamma\gamma\to J/\psi \omega$
Experiment: BABAR (PEP-II)
Inspire ID: 1122034
Status: VALIDATED NOHEPDATA
Authors:
  • Peter Richardson
References:
  • Phys.Rev.D 86 (2012) 072002
Beams: e+ e-
Beam energies: (5.3, 5.3) GeV
Run details:
  • e+ e- > e+ e- gamma gamma with gamma gamma -> X(3915)

Measurement of the mass and angle distributions in $\gamma\gamma\to J/\psi \omega$.

Source code: BABAR_2012_I1122034.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/UnstableParticles.hh"
  5#include "Rivet/Projections/Beam.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// @brief gamma gamma -> X(3915) -> J/psi omega
 11  class BABAR_2012_I1122034 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(BABAR_2012_I1122034);
 16
 17
 18    /// @name Analysis methods
 19    /// @{
 20
 21    /// Book histograms and initialise projections before the run
 22    void init() {
 23      // Initialise and register projections
 24      declare(Beam(), "Beams");
 25      declare(FinalState(),"FS");
 26      declare(UnstableParticles(Cuts::pid==223 or Cuts::pid==443), "UFS");
 27      // histograms
 28      book(_h_mass,1,1,1);
 29      for (unsigned int ix=0;ix<4;++ix) {
 30        book(_h_angle1[ix],2,1,1+ix);
 31        if (ix<3) book(_h_angle2[ix],3,1,1+ix);
 32      }
 33    }
 34
 35    void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
 36      for (const Particle &child : p.children()) {
 37        if (child.children().empty()) {
 38          --nRes[child.pid()];
 39          --ncount;
 40        } else {
 41          findChildren(child,nRes,ncount);
 42        }
 43      }
 44    }
 45
 46    bool findScattered(Particle beam, double& q2) {
 47      bool found = false;
 48      Particle scat = beam;
 49      while (!scat.children().empty()) {
 50        found = false;
 51        for (const Particle & p : scat.children()) {
 52          if (p.pid()==scat.pid()) {
 53            scat=p;
 54            found=true;
 55            break;
 56          }
 57        }
 58        if (!found) break;
 59      }
 60      if (!found) return false;
 61      q2 = -(beam.mom() - scat.mom()).mass2();
 62      return true;
 63    }
 64
 65    void findChildren(const Particle & p, Particles & pim, Particles & pip,
 66		      Particles & pi0, unsigned int &ncount) {
 67      for (const Particle &child : p.children()) {
 68        if(child.pid()==PID::PIPLUS) {
 69          pip.push_back(child);
 70          ncount+=1;
 71        }
 72        else if(child.pid()==PID::PIMINUS) {
 73          pim.push_back(child);
 74          ncount+=1;
 75        }
 76        else if(child.pid()==PID::PI0) {
 77          pi0.push_back(child);
 78          ncount+=1;
 79        }
 80        else if(child.children().empty()) {
 81          ncount+=1;
 82        }
 83        else {
 84          findChildren(child,pim,pip,pi0,ncount);
 85        }
 86      }
 87    }
 88
 89    /// Perform the per-event analysis
 90    void analyze(const Event& event) {
 91      // find scattered leptons and calc Q2
 92      const Beam& beams = apply<Beam>(event, "Beams");
 93      double q12 = -1, q22 = -1;
 94      if (!findScattered(beams.beams().first,  q12)) vetoEvent;
 95      if (!findScattered(beams.beams().second, q22)) vetoEvent;
 96      // check the final state
 97      const FinalState& fs = apply<FinalState>(event, "FS");
 98      map<long,int> nCount;
 99      int ntotal(0);
100      for (const Particle& p : fs.particles()) {
101        nCount[p.pid()] += 1;
102        ++ntotal;
103      }
104      // find the J/psi
105      const FinalState& ufs = apply<FinalState>(event, "UFS");
106      Particle omega,psi;
107      bool found=false;
108      for (const Particle &p1 : ufs.particles(Cuts::pid==443)) {
109        if (p1.children().empty()) continue;
110        map<long,int> nRes = nCount;
111        int ncount = ntotal;
112        findChildren(p1,nRes,ntotal);
113        for (const Particle& p2 : ufs.particles(Cuts::pid==223)) {
114          if (p2.children().empty()) continue;
115          map<long,int> nRes2 = nRes;
116          int ncount2 = ncount;
117          findChildren(p2,nRes2,ncount2);
118          found = true;
119          for (const auto& val : nRes2) {
120            if (abs(val.first)==11) {
121              if (val.second!=1) {
122                found = false;
123                break;
124              }
125            }
126            else if(val.second!=0) {
127              found = false;
128              break;
129            }
130          }
131          if (found) {
132            psi   = p1;
133            omega = p2;
134            break;
135          }
136        }
137      }
138      if (!found) vetoEvent;
139      FourMomentum psum = omega.mom()+psi.mom();
140      if (psum.pT()>0.2) vetoEvent;
141      // mass distribution
142      _h_mass->fill(psum.mass());
143      // from now on we need specific decay modes of J/psi and omega
144      // first J/psi -> l+l-
145      if (psi.children().size()!=2) vetoEvent;
146      if (psi.children()[0].pid()!=-psi.children()[1].pid()) vetoEvent;
147      if (psi.children()[0].abspid()!=11 && psi.children()[0].abspid()!=13) vetoEvent;
148      Particle ep = psi.children()[0];
149      Particle em = psi.children()[1];
150      if (ep.pid()>0) swap(ep,em);
151      // omega decay
152      unsigned int ncount=0;
153      Particles pip,pim,pi0;
154      findChildren(omega,pim,pip,pi0,ncount);
155      if (ncount!=3 || !(pim.size()==1 && pip.size()==1 && pi0.size()==1)) vetoEvent;
156      // boost to gamma gamma frame
157      LorentzTransform boostCMS = LorentzTransform::mkFrameTransformFromBeta(psum.betaVec());
158      FourMomentum pPsi   = boostCMS.transform(psi   .mom());
159      FourMomentum pOmega = boostCMS.transform(omega .mom());
160      FourMomentum pLp    = boostCMS.transform(ep    .mom());
161      FourMomentum pPip   = boostCMS.transform(pip[0].mom());
162      FourMomentum pPim   = boostCMS.transform(pim[0].mom());
163      Vector3 axis(0.,0.,1.);
164      // lepton angle
165      double cosL = pLp.p3().unit().dot(axis);
166      _h_angle1[0]->fill(cosL);
167      LorentzTransform boostOmega = LorentzTransform::mkFrameTransformFromBeta(pOmega.betaVec());
168      pPip   = boostOmega.transform(pPip);
169      pPim   = boostOmega.transform(pPim);
170      // omega decay plane normal angle
171      Vector3 axisOmega = pPip.p3().cross(pPim.p3()).unit();
172      double cosN = axisOmega.dot(axis);
173      _h_angle1[1]->fill(cosN);
174      // angle lepton and omega
175      _h_angle1[2]->fill(pLp.p3().unit().dot(axisOmega));
176      // helicity angle
177      _h_angle1[3]->fill(pPsi.p3().unit().dot(psum.p3().unit()));
178      // now the new frame
179      Vector3 axisZ = pOmega.p3().unit();
180      Vector3 axisY = axisZ.cross(axisOmega);
181      Vector3 axisX = axisY.cross(axisZ);
182      // second set of angles
183      LorentzTransform boostPsi = LorentzTransform::mkFrameTransformFromBeta(pPsi.betaVec());
184      Vector3 axisL = boostPsi.transform(pLp).p3().unit();
185      _h_angle2[0]->fill(axisZ.dot(axisOmega));
186      _h_angle2[1]->fill(axisL.dot(pPsi.p3().unit()));
187      axisZ *=-1.;
188      axisX *=-1.;
189      Vector3 axisnp = axisL.cross(axisZ).unit();
190      double phiL = atan2(axisnp.dot(axisY),axisnp.dot(axisX));
191      double phiN = atan2(axisOmega.dot(axisY),axisOmega.dot(axisX));
192      _h_angle2[2]->fill(mapAngleMPiToPi(phiL-phiN)/M_PI*180.);
193    }
194
195
196    /// Normalise histograms etc., after the run
197    void finalize() {
198      normalize(_h_mass, 1.0, false);
199	    normalize(_h_angle1, 1.0, false);
200	    normalize(_h_angle2, 1.0, false);
201    }
202
203    /// @}
204
205
206    /// @name Histograms
207    /// @{
208    Histo1DPtr _h_mass,_h_angle1[4],_h_angle2[3];
209    /// @}
210
211
212  };
213
214
215  RIVET_DECLARE_PLUGIN(BABAR_2012_I1122034);
216
217}