Rivet analyses referenceBABAR_2018_I1667191Mass and angular distributions in the radiative decays Υ(1S)→π+π− and K+K−Experiment: BABAR (PEP-II) Inspire ID: 1667191 Status: VALIDATED NOHEPDATA Authors:
Beam energies: ANY Run details:
Mass and angular distributions in the radiative decays Υ(1S)→π+π− and K+K− Source code: BABAR_2018_I1667191.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/UnstableParticles.hh"
4#include "Rivet/Projections/DecayedParticles.hh"
5
6namespace Rivet {
7
8
9 /// @brief Upsilon(1S) -> gamma pi+pi- K+K-
10 class BABAR_2018_I1667191 : public Analysis {
11 public:
12
13 /// Constructor
14 RIVET_DEFAULT_ANALYSIS_CTOR(BABAR_2018_I1667191);
15
16
17 /// @name Analysis methods
18 /// @{
19
20 /// Book histograms and initialise projections before the run
21 void init() {
22 // projections
23 UnstableParticles ufs = UnstableParticles(Cuts::abspid==200553 or
24 Cuts::abspid==100553);
25 declare(ufs, "UFS");
26 DecayedParticles UPS(ufs);
27 UPS.addStable( 553);
28 declare(UPS, "UPS");
29 // histograms
30 for(unsigned int ix=0;ix<2;++ix)
31 book(_h_mass[ix],1+ix,1,1);
32 book(_h_pi,3,1,1);
33 for(unsigned int ix=0;ix<3;++ix)
34 for(unsigned int iy=0;iy<2;++iy)
35 book(_h_angle[ix][iy],4+ix,1,1+iy);
36 }
37
38
39 /// Recursively walk the decay tree to find decay products of @a p
40 void findDecayProducts(Particle mother, Particles& gamma,
41 Particles & pip, Particles & pim,
42 Particles & Kp , Particles & Km,unsigned int & nstable) {
43 for(const Particle & p: mother.children()) {
44 if (p.pid()== 211) pip.push_back(p);
45 else if(p.pid()==-211) pim.push_back(p);
46 else if(p.pid()== 321) Kp .push_back(p);
47 else if(p.pid()==-321) Km .push_back(p);
48 else if(p.pid()==22) gamma.push_back(p);
49 else if(p.children().empty())
50 nstable+=1;
51 else {
52 findDecayProducts(p, gamma,pip,pim,Kp,Km,nstable);
53
54 }
55 }
56 }
57
58 /// Perform the per-event analysis
59 void analyze(const Event& event) {
60 static const map<PdgId,unsigned int> & mode = { { 553,1},{ 211,1}, {-211,1}};
61 DecayedParticles UPS = apply<DecayedParticles>(event, "UPS");
62 // loop over particles
63 for(unsigned int ix=0;ix<UPS.decaying().size();++ix) {
64 // check pi+pi- upslion(1S) decay mode
65 if (!UPS.modeMatches(ix,3,mode)) continue;
66 const Particle & ups1 = UPS.decayProducts()[ix].at( 553)[0];
67 const Particle & pips = UPS.decayProducts()[ix].at( 211)[0];
68 const Particle & pims = UPS.decayProducts()[ix].at(-211)[0];
69 // boost to rest frame
70 LorentzTransform boost;
71 if (UPS.decaying()[ix].p3().mod() > 1*MeV)
72 boost = LorentzTransform::mkFrameTransformFromBeta(UPS.decaying()[ix].momentum().betaVec());
73 FourMomentum ppipi = boost.transform(pips.momentum()+pims.momentum());
74 Vector3 axis1 = ppipi.p3().unit();
75 LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(ppipi);
76 FourMomentum ppi = boost1.transform(boost.transform(pips.momentum()));
77 _h_pi->fill(abs(ppi.p3().unit().dot(axis1)));
78 unsigned int nstable=0;
79 Particles gamma,pip,pim,Kp,Km;
80 findDecayProducts(ups1, gamma,pip,pim,Kp,Km,nstable);
81 if(gamma.size()!=1 || nstable!=0) continue;
82 // gamma pi+pi-
83 if(Kp.empty()&&Km.empty()&&pip.size()==1&&pim.size()==1) {
84 ppipi = pip[0].momentum()+pim[0].momentum();
85 double mpipi = ppipi.mass();
86 _h_mass[0]->fill(mpipi);
87 axis1 *=-1;
88 LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(boost.transform(ups1.momentum()).betaVec());
89 FourMomentum pgamma = boost2.transform(boost.transform(gamma[0].momentum()));
90 Vector3 axis2 = pgamma.p3().unit();
91 double cGamma = axis2.dot(axis1);
92 ppipi = boost2.transform(boost1.transform(ppipi));
93 LorentzTransform boost3 = LorentzTransform::mkFrameTransformFromBeta(ppipi.betaVec());
94 Vector3 axis3 = boost3.transform(boost2.transform(boost1.transform(pip[0].momentum()))).p3().unit();
95 double cH = axis3.dot(axis2);
96 int iloc=-1;
97 if(mpipi>0.6&&mpipi<1.) iloc=0;
98 else if(mpipi>1.092&&mpipi<1.46) iloc=1;
99 if(iloc>=0) {
100 _h_angle[iloc][0]->fill(cGamma);
101 _h_angle[iloc][1]->fill(cH);
102 }
103 }
104 // gamma K+K-
105 else if (pip.empty()&&pim.empty()&&Kp.size()==1&&Km.size()==1) {
106 FourMomentum pKK = Kp[0].momentum()+Km[0].momentum();
107 double mKK = pKK.mass();
108 _h_mass[1]->fill(mKK);
109 axis1 *=-1;
110 LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(boost.transform(ups1.momentum()).betaVec());
111 FourMomentum pgamma = boost2.transform(boost.transform(gamma[0].momentum()));
112 Vector3 axis2 = pgamma.p3().unit();
113 double cGamma = axis2.dot(axis1);
114 pKK = boost2.transform(boost1.transform(pKK));
115 LorentzTransform boost3 = LorentzTransform::mkFrameTransformFromBeta(pKK.betaVec());
116 Vector3 axis3 = boost3.transform(boost2.transform(boost1.transform(Kp[0].momentum()))).p3().unit();
117 double cH = axis3.dot(axis2);
118 if(mKK>1.424 && mKK<1.62) {
119 _h_angle[2][0]->fill(cGamma);
120 _h_angle[2][1]->fill(cH);
121 }
122 }
123 }
124 }
125
126
127 /// Normalise histograms etc., after the run
128 void finalize() {
129 for(unsigned int ix=0;ix<2;++ix)
130 normalize(_h_mass[ix],1.,false);
131 normalize(_h_pi,1.,false);
132 for(unsigned int ix=0;ix<3;++ix)
133 for(unsigned int iy=0;iy<2;++iy)
134 normalize(_h_angle[ix][iy],1.,false);
135 }
136
137 /// @}
138
139
140 /// @name Histograms
141 /// @{
142 Histo1DPtr _h_mass[2];
143 Histo1DPtr _h_pi;
144 Histo1DPtr _h_angle[3][2];
145 /// @}
146
147 };
148
149
150 RIVET_DECLARE_PLUGIN(BABAR_2018_I1667191);
151
152}
|
