Rivet analyses referenceBESIII_2012_I1079921Mass and angular distributions in $J/\psi\to\gamma p\bar{p}$ near the $p\bar{p}$ thresholdExperiment: BESIII (BEPC) Inspire ID: 1079921 Status: VALIDATED NOHEPDATA Authors:
Beam energies: (1.6, 1.6); (1.8, 1.8) GeV Run details:
Mass and angular distributions in $J/\psi\to\gamma p\bar{p}$ near the $p\bar{p}$ threshold Source code: BESIII_2012_I1079921.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/Beam.hh"
4#include "Rivet/Projections/UnstableParticles.hh"
5#include "Rivet/Projections/DecayedParticles.hh"
6
7namespace Rivet {
8
9
10 /// @brief J/psi -> gamma p pbar
11 class BESIII_2012_I1079921 : public Analysis {
12 public:
13
14 /// Constructor
15 RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2012_I1079921);
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 UnstableParticles ufs = UnstableParticles(Cuts::abspid==443 or
25 Cuts::abspid==100443);
26 declare(ufs, "UFS");
27 DecayedParticles PSI(ufs);
28 declare(PSI, "PSI");
29 declare(Beam(), "Beams");
30 // book histos
31 if (isCompatibleWithSqrtS(3.1,0.001)) {
32 for(unsigned int ix=0;ix<4;++ix)
33 book(_h[ix],1,1,1+ix);
34 }
35 else if(isCompatibleWithSqrtS(3.686,0.001)) {
36 book(_h[0],2,1,1);
37 }
38 else {
39 cerr << "testing problem " << sqrtS() << "\n";
40 throw Error("Unexpected sqrtS ! Only 3.1 and 3.686 GeV are supported");
41 }
42 }
43
44 // angle cuts due regions of BES calorimeter
45 bool vetoPhoton(const double & cTheta) {
46 return cTheta>0.92 || (cTheta>0.8 && cTheta<0.86);
47 }
48
49 /// Perform the per-event analysis
50 void analyze(const Event& event) {
51 // get the axis, direction of incoming electron
52 const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
53 Vector3 axis;
54 if(beams.first.pid()>0)
55 axis = beams.first .momentum().p3().unit();
56 else
57 axis = beams.second.momentum().p3().unit();
58 // find the J/psi decays
59 static const map<PdgId,unsigned int> & mode = { { 2212,1}, { -2212,1},{ 22,1}};
60 DecayedParticles PSI = apply<DecayedParticles>(event, "PSI");
61 if( PSI.decaying().size()!=1) vetoEvent;
62 if(sqrtS()>3.2 && PSI.decaying()[0].pid()==443) vetoEvent;
63 if(!PSI.modeMatches(0,3,mode)) vetoEvent;
64 const Particle & pp = PSI.decayProducts()[0].at( 2212)[0];
65 const Particle & pbar = PSI.decayProducts()[0].at(-2212)[0];
66 const Particle & gam = PSI.decayProducts()[0].at( 22)[0];
67 double mass = (pp.momentum()+pbar.momentum()).mass()-pp.mass()-pbar.mass();
68 _h[0]->fill(mass);
69 if(!_h[1] || mass>0.05) return;
70 double cTheta = axis.dot(gam.p3().unit());
71 // photon angle
72 if(vetoPhoton(abs(cTheta))) vetoEvent;
73 _h[1]->fill(cTheta);
74 // remaining angles
75 LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(PSI.decaying()[0].momentum().betaVec());
76 FourMomentum pGamma = boost1.transform(gam.momentum());
77 FourMomentum pppbar = boost1.transform(pp.momentum()+pbar.momentum());
78 Vector3 e1z = pGamma.p3().unit();
79 Vector3 e1y = e1z.cross(axis).unit();
80 Vector3 e1x = e1y.cross(e1z).unit();
81 LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(pppbar.betaVec());
82 Vector3 axis2 = boost2.transform(boost1.transform(pp.momentum())).p3().unit();
83 _h[2]->fill(e1z.dot(axis2));
84 double phi = atan2(axis2.dot(e1y),axis2.dot(e1x))/M_PI*180.;
85 _h[3]->fill(phi);
86 }
87
88
89 /// Normalise histograms etc., after the run
90 void finalize() {
91 for(unsigned int ix=0;ix<4;++ix)
92 if(_h[ix]) normalize(_h[ix],1.,false);
93 }
94
95 /// @}
96
97
98 /// @name Histograms
99 /// @{
100 Histo1DPtr _h[4];
101 /// @}
102
103
104 };
105
106
107 RIVET_DECLARE_PLUGIN(BESIII_2012_I1079921);
108
109}
|