Rivet analyses referenceBABAR_2015_I1388182Charged particle asymmetries in $e^+e^-\to\mu^+\mu^-\gamma$ and $e^+e^-\to\pi^+\pi^-\gamma$Experiment: BABAR (PEP-II) Inspire ID: 1388182 Status: VALIDATED NOHEPDATA SINGLEWEIGHT Authors:
Beam energies: (3.5, 8.0); (5.3, 5.3) GeV Run details:
Measurement of Charged particle asymmetries in $e^+e^-\to\mu^+\mu^-\gamma$ and $e^+e^-\to\pi^+\pi^-\gamma$ which are senistive to interference between the initial- and final-state QED radiation Source code: BABAR_2015_I1388182.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FinalState.hh"
4#include "Rivet/Projections/Beam.hh"
5
6namespace Rivet {
7
8
9 /// @brief e+ e- > mu+ mu- gamma or pi+ pi- gamma
10 class BABAR_2015_I1388182 : public Analysis {
11 public:
12
13 /// Constructor
14 RIVET_DEFAULT_ANALYSIS_CTOR(BABAR_2015_I1388182);
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(Beam(), "Beams");
24 declare(FinalState(),"FS");
25 // histograms
26 Profile1DPtr tmp;
27 book(tmp,"TMP/pi0",refData(3,1,1));
28 _h_pipi.push_back(tmp);
29 for (unsigned int ix=0; ix<15; ++ix) {
30 book(tmp,3,1,1+ix);
31 _h_pipi.push_back(tmp);
32 if (ix>13) continue;
33 book(tmp,1,1,1+ix);
34 _h_mumu.push_back(tmp);
35 }
36 book(tmp,"TMP/pi16",refData(3,1,15));
37 _h_pipi.push_back(tmp);
38 book(tmp,"TMP/pi17",refData(3,1,15));
39 _h_pipi.push_back(tmp);
40 book(_h_nopsi,"TMP/nopsi",refData(1,1,7));
41 }
42
43
44 /// Perform the per-event analysis
45 void analyze(const Event& event) {
46 // get the axis, direction of incoming positron
47 const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
48 bool CMF = fuzzyEquals(-beams.first .mom().z()/beams.second.mom().z(),1);
49 Vector3 axis1 = beams.first .mom().p3().unit();
50 Vector3 axis2 = beams.second.mom().p3().unit();
51 if(beams.first.pid()<0) swap(axis1,axis2);
52 // find the final state final state particles
53 Particle mup,mum,pip,pim,gamma;
54 Particles fs = apply<FinalState>(event,"FS").particles();
55 // boost to CMF frame
56 LorentzTransform boost;
57 if (!CMF) {
58 boost = LorentzTransform::mkFrameTransformFromBeta((beams.first.mom()+
59 beams.second.mom()).betaVec());
60 }
61 else {
62 double E1=3.5,E2 = 0.25*sqr(sqrtS())/E1;
63 FourMomentum pnew(E1+E2,0,0,E1-E2);
64 boost = LorentzTransform::mkFrameTransformFromBeta(-pnew.betaVec());
65 }
66 FourMomentum pGamma;
67 for (const Particle& p : fs) {
68 FourMomentum pLab,pCMF;
69 if (CMF) {
70 pCMF = p.mom();
71 pLab = boost.transform(p.mom());
72 }
73 else {
74 pCMF = boost.transform(p.mom());
75 pLab = p.mom();
76 }
77 double theta = acos(pLab.p3().unit().dot(axis1));
78 if (p.isCharged()) {
79 if (theta<.4 || theta>2.45) continue;
80 if (pLab.p3().mod()<1.) continue;
81 if (p.pid()==PID::MUON && mum.pid()!=PID::MUON) {
82 mum = p;
83 }
84 else if (p.pid()==PID::ANTIMUON && mup.pid()!=PID::ANTIMUON) {
85 mup = p;
86 }
87 else if (p.pid()==PID::PIPLUS && pip.pid()!=PID::PIPLUS) {
88 pip = p;
89 }
90 else if (p.pid()==PID::PIMINUS && pim.pid()!=PID::PIMINUS) {
91 pim = p;
92 }
93 }
94 else if(p.pid()==PID::GAMMA) {
95 // angle cut on the photon
96 if (theta<.35 || theta>2.4) continue;
97 if (gamma.pid()!=PID::GAMMA) {
98 gamma = p;
99 pGamma = pCMF;
100 }
101 else {
102 if(pCMF.E()>pGamma.E()) {
103 gamma = p;
104 pGamma = pCMF;
105 }
106 }
107 }
108 else {
109 vetoEvent;
110 }
111 }
112 if (gamma.pid()!=PID::GAMMA) vetoEvent;
113 if (!( ((pip.pid()==PID::PIPLUS && pim.pid()==PID::PIMINUS ) ||
114 (mum.pid()==PID::MUON && mup.pid()==PID::ANTIMUON)))) {
115 vetoEvent;
116 }
117 if (pip.pid()==PID::PIPLUS && mum.pid()==PID::MUON) vetoEvent;
118 if (pGamma.E()<3.) vetoEvent;
119 Vector3 axisZ = pGamma.p3().unit();
120 Vector3 axisX = (axis2-axisZ.dot(axis2)*axisZ).unit();
121 Vector3 axisY = axisZ.cross(axisX);
122 FourMomentum pMinus,pPlus;
123 if (CMF) {
124 if (mum.pid()==PID::MUON) {
125 pMinus = mum.mom();
126 pPlus = mup.mom();
127 }
128 else {
129 pMinus = pim.mom();
130 pPlus = pip.mom();
131 }
132 }
133 else {
134 if (mum.pid()==PID::MUON) {
135 pMinus = boost.transform(mum.mom());
136 pPlus = boost.transform(mup.mom());
137 }
138 else {
139 pMinus = boost.transform(pim.mom());
140 pPlus = boost.transform(pip.mom());
141 }
142 }
143 double phiM = atan2(pMinus.p3().dot(axisY),pMinus.p3().dot(axisX));
144 if (phiM<0.) phiM+=2.*M_PI;
145 double phiP = atan2(pPlus .p3().dot(axisY),pPlus .p3().dot(axisX));
146 if (phiP<0.) phiP+=2.*M_PI;
147 double mass = (pMinus+pPlus).mass();
148 if (mum.pid()==PID::MUON) {
149 if (mass>0.2 && mass<7.) {
150 unsigned int imass = int(mass/.5);
151 if (phiM<M_PI) _h_mumu[imass]->fill(cos(phiM), 1.);
152 else _h_mumu[imass]->fill(cos(phiP),-1.);
153 if (imass==6 && mass>3.2) {
154 if (phiM<M_PI) _h_nopsi->fill(cos(phiM), 1.);
155 else _h_nopsi->fill(cos(phiP),-1.);
156 }
157 }
158 }
159 else if(pip.pid()==PID::PIPLUS) {
160 if (mass>0.2 && mass<2.0) {
161 unsigned int imass = int((mass-0.2)/.1);
162 if(phiM<M_PI) _h_pipi[imass]->fill(cos(phiM), 1.);
163 else _h_pipi[imass]->fill(cos(phiP),-1.);
164 }
165 }
166 }
167
168 pair<double,double> calcA0(Profile1DPtr hist) {
169 if(hist->numEntries()==0.) return make_pair(0.,0.);
170 double sum1(0.),sum2(0.);
171 for (const auto& bin : hist->bins() ) {
172 if (bin.numEntries()<2) continue;
173 double Oi = bin.mean(2);
174 double Ei = bin.stdErr(2);
175 if (Ei==0.) continue;
176 double xi = 0.5*(bin.xMin()+bin.xMax());
177 sum1 += sqr(xi/Ei);
178 sum2 += Oi*xi/sqr(Ei);
179 }
180 return make_pair(sum2/sum1,sqrt(1./sum1));
181 }
182
183 /// Normalise histograms etc., after the run
184 void finalize() {
185 // muon assymetries
186 Estimate1DPtr _A_mumu;
187 book(_A_mumu,2,1,1);
188 for (unsigned int ix=0;ix<14;++ix) {
189 pair<double,double> A0 = calcA0(ix!=6 ? _h_mumu[ix] : _h_nopsi);
190 _A_mumu->bin(ix+1).set(A0.first, A0.second);
191 }
192 // pion assymetries
193 Estimate1DPtr _A_pipi;
194 book(_A_pipi,4,1,1);
195 for (unsigned int ix=0;ix<_h_pipi.size();++ix) {
196 pair<double,double> A0 = calcA0(_h_pipi[ix]);
197 _A_pipi->bin(ix+1).set(A0.first, A0.second);
198 }
199 }
200
201 /// @}
202
203
204 /// @name Histograms
205 /// @{
206 vector<Profile1DPtr> _h_mumu, _h_pipi;
207 Profile1DPtr _h_nopsi;
208 /// @}
209
210
211 };
212
213
214 RIVET_DECLARE_PLUGIN(BABAR_2015_I1388182);
215
216}
|
