Rivet analyses referenceBABAR_2015_I1377201Azimuthal asymmetries in inclusive $\pi\pi$ $KK$ and $K\pi$ pairs at 10.58 GeVExperiment: BABAR (PEP-II) Inspire ID: 1377201 Status: VALIDATED Authors:
Beam energies: (5.3, 5.3) GeV Run details:
Measurement of azimuthal asymmetries in inclusive $\pi\pi$ $KK$ and $K\pi$ pair production at $\sqrt{s}=10.58$ GeV by the BABAR experiment Source code: BABAR_2015_I1377201.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/Thrust.hh"
4#include "Rivet/Projections/FinalState.hh"
5#include "Rivet/Projections/Beam.hh"
6
7namespace Rivet {
8
9
10 /// @brief azimuthal asymmetries in pipi Kpi and KK
11 class BABAR_2015_I1377201 : public Analysis {
12 public:
13
14 /// Constructor
15 RIVET_DEFAULT_ANALYSIS_CTOR(BABAR_2015_I1377201);
16
17
18 /// @name Analysis methods
19 /// @{
20
21 /// Book histograms and initialise projections before the run
22 void init() {
23 // projections
24 const FinalState fs;
25 declare(fs,"FS");
26 declare(Thrust(fs),"Thrust");
27 declare(Beam(), "Beams");
28 // declare the histos for the distributions
29 string type [3] = {"KK","Kpi","pipi"};
30 string charge[3] = {"Like","Opposite","All"};
31 unsigned int nbin=20;
32 for(unsigned int itype=0;itype<3;++itype) {
33 for(unsigned int icharge=0;icharge<3;++icharge) {
34 for(unsigned int ibin=0;ibin<16;++ibin) {
35 std::ostringstream title1;
36 title1 << "/TMP/h_thrust" << type[itype] << "_" << charge[icharge] << "_" << ibin+1;
37 book(_h_thrust[itype][icharge][ibin],title1.str(),nbin,0.,M_PI);
38 std::ostringstream title2;
39 title2 << "/TMP/h_hadron" << type[itype] << "_" << charge[icharge] << "_" << ibin+1;
40 book(_h_hadron[itype][icharge][ibin],title2.str(),nbin,0.,M_PI);
41 }
42 }
43 }
44 }
45
46 unsigned int iBin(double z) {
47 if (z<.2) return 0;
48 else if(z<.3) return 1;
49 else if(z<.5) return 2;
50 else return 3;
51 }
52
53 /// Perform the per-event analysis
54 void analyze(const Event& event) {
55 // get the axis, direction of incoming electron
56 const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
57 Vector3 axis1;
58 if(beams.first.pid()>0)
59 axis1 = beams.first .momentum().p3().unit();
60 else
61 axis1 = beams.second.momentum().p3().unit();
62 // apply thrust cuts T > 0.8 and | cos θ th | < 0.6
63 Thrust thrust = apply<Thrust>(event,"Thrust");
64 if(thrust.thrust()<=0.8) vetoEvent;
65 if(cos(thrust.thrustAxis().polarAngle())>=0.6) vetoEvent;
66 // construct x,y,z axes for thrust defn
67 ThreeVector t_z = thrust.thrustAxis();
68 ThreeVector t_x = (axis1-t_z.dot(axis1)*t_z).unit();
69 ThreeVector t_y = t_z.cross(t_x);
70 // loop over the particles
71 Particles charged = apply<FinalState>(event,"FS").particles(Cuts::abspid==PID::PIPLUS || Cuts::abspid==PID::KPLUS);
72 for(unsigned int ix=0;ix<charged.size();++ix) {
73 // z and angle cut
74 const double x1=2.*charged[ix].momentum().t()/sqrtS();
75 if(x1<0.16||x1>.9) continue;
76 if(abs(t_z.angle(charged[ix].momentum().p3()))>0.25*M_PI) continue;
77 double dot1 = t_z.dot(charged[ix].p3());
78 for(unsigned int iy=ix+1;iy<charged.size();++iy) {
79 const double x2=2.*charged[iy].momentum().t()/sqrtS();
80 // z and angle cut
81 if(x2<0.16||x2>.9) continue;
82 if(abs(t_z.angle(charged[ix].momentum().p3()))>0.25*M_PI) continue;
83 // different hemi
84 double dot2 = t_z.dot(charged[iy].p3());
85 if(dot1*dot2>0.) continue;
86 Particle p1=charged[ix], p2=charged[iy];
87 double z1(x1),z2(x2);
88 // randomly order the particles
89 if(rand()/static_cast<double>(RAND_MAX) < 0.5 ) {
90 swap(p1,p2);
91 swap(z1,z2);
92 }
93 // thrust def
94 double phi12 = atan2(p1.p3().dot(t_y),p1.p3().dot(t_x))+atan2(p2.p3().dot(t_y),p2.p3().dot(t_x));
95 if(phi12>M_PI) phi12 -= 2*M_PI;
96 if(phi12<-M_PI) phi12 += 2*M_PI;
97 if(phi12<0.) phi12 = -phi12;
98 // hadron defn
99 ThreeVector h_z = p2.p3().unit();
100 ThreeVector h_x = (axis1 - h_z.dot(axis1)*h_z).unit();
101 ThreeVector pt1 = p1.p3() - h_z.dot(p1.p3())*h_z;
102 double phi0 = pt1.angle(h_x);
103 if(phi0>M_PI) phi0 -= 2*M_PI;
104 if(phi0<-M_PI) phi0 += 2*M_PI;
105 int ibin = 4*iBin(z1)+iBin(z2);
106 // pi pi
107 if(p1.abspid()==PID::PIPLUS && p2.abspid()==PID::PIPLUS) {
108 if(p1.pid()==p2.pid()) {
109 _h_thrust[2][0][ibin]->fill(phi12);
110 _h_hadron[2][0][ibin]->fill(phi0);
111 }
112 else {
113 _h_thrust[2][1][ibin]->fill(phi12);
114 _h_hadron[2][1][ibin]->fill(phi0);
115 }
116 _h_thrust[2][2][ibin]->fill(phi12);
117 _h_hadron[2][2][ibin]->fill(phi0);
118 }
119 // K K
120 else if(p1.abspid()==PID::KPLUS && p2.abspid()==PID::KPLUS) {
121 if(p1.pid()==p2.pid()) {
122 _h_thrust[0][0][ibin]->fill(phi12);
123 _h_hadron[0][0][ibin]->fill(phi0);
124 }
125 else {
126 _h_thrust[0][1][ibin]->fill(phi12);
127 _h_hadron[0][1][ibin]->fill(phi0);
128 }
129 _h_thrust[0][2][ibin]->fill(phi12);
130 _h_hadron[0][2][ibin]->fill(phi0);
131 }
132 // K pi
133 else {
134 if(p1.pid()*p2.pid()>0) {
135 _h_thrust[1][0][ibin]->fill(phi12);
136 _h_hadron[1][0][ibin]->fill(phi0);
137 }
138 else {
139 _h_thrust[1][1][ibin]->fill(phi12);
140 _h_hadron[1][1][ibin]->fill(phi0);
141 }
142 _h_thrust[1][2][ibin]->fill(phi12);
143 _h_hadron[1][2][ibin]->fill(phi0);
144 }
145 }
146 }
147 }
148
149 pair<double,double> calcAsymmetry(Estimate1DPtr hist, double fact=1.) {
150 double sum1(0.),sum2(0.);
151 for (auto& b : hist->bins() ) {
152 double Oi = b.val();
153 if(Oi==0. || std::isnan(Oi) ) continue;
154 double ai = 1.;
155 double bi = (sin(fact*b.xMax())-sin(fact*b.xMin()))/(b.xWidth())/fact;
156 double Ei = b.errAvg();
157 sum1 += sqr(bi/Ei);
158 sum2 += bi/sqr(Ei)*(Oi-ai);
159 }
160 if(sum1==0.) return make_pair(0.,0.);
161 return make_pair(sum2/sum1*1e4,sqrt(1./sum1)*1e4);
162 }
163
164 /// Normalise histograms etc., after the run
165 void finalize() {
166 for(unsigned int itype=0;itype<3;++itype) {
167 for(unsigned int icharge=0;icharge<3;++icharge) {
168 for(unsigned int ibin=0;ibin<16;++ibin) {
169 normalize(_h_thrust[itype][icharge][ibin]);
170 normalize(_h_hadron[itype][icharge][ibin]);
171 }
172 }
173 }
174 // construct ther ratios
175 // declare the histos for the distributions
176 string type [3] = {"pipi","Kpi","KK"};
177 string charge[3] = {"Like","Opposite","All"};
178 for(unsigned int itype=0;itype<3;++itype) {
179 Estimate2DPtr h3_thrust_UL;
180 book(h3_thrust_UL,2*itype+1,1,2);
181 Estimate2DPtr h3_thrust_UC;
182 book(h3_thrust_UC,2*itype+1,1,3);
183 Estimate2DPtr h3_hadron_UL;
184 book(h3_hadron_UL,2*itype+2,1,2);
185 Estimate2DPtr h3_hadron_UC;
186 book(h3_hadron_UC,2*itype+2,1,3);
187
188 unsigned int ihist=1;
189 Estimate1DPtr h2_thrust_UL;
190 book(h2_thrust_UL,7+2*itype,ihist,2);
191 Estimate1DPtr h2_thrust_UC;
192 book(h2_thrust_UC,7+2*itype,ihist,3);
193 Estimate1DPtr h2_hadron_UL;
194 book(h2_hadron_UL,8+2*itype,ihist,2);
195 Estimate1DPtr h2_hadron_UC;
196 book(h2_hadron_UC,8+2*itype,ihist,3);
197
198 Estimate2D temphisto1(refData<Estimate2D>(2*itype+1, 1, 2));
199 Estimate2D temphisto2(refData<Estimate2D>(2*itype+2, 1, 2));
200 for (unsigned int ibin=0;ibin<16;++ibin) {
201 //const auto& b1 = temphisto1.bin(ibin);
202 //const auto& b2 = temphisto2.bin(ibin);
203 if (ibin>0 && ibin%4==0) {
204 ++ihist;
205 book(h2_thrust_UL,7+2*itype,ihist,2);
206 book(h2_thrust_UC,7+2*itype,ihist,3);
207 book(h2_hadron_UL,8+2*itype,ihist,2);
208 book(h2_hadron_UC,8+2*itype,ihist,3);
209 }
210 // thrust direction
211 // opposite/like sign
212 std::ostringstream title1;
213 title1 << "/TMP/R_thrust_" << type[itype] << "_UL_" << ibin+1;
214 Estimate1DPtr htemp;
215 book(htemp,title1.str());
216 divide(_h_thrust[itype][1][ibin],
217 _h_thrust[itype][0][ibin],htemp);
218 pair<double,double> asym = calcAsymmetry(htemp);
219 h3_thrust_UL->bin(ibin+1).set(asym.first, asym.second);
220 h2_thrust_UL->bin(ibin+1).set(asym.first, asym.second);
221 // opposite/all sign
222 std::ostringstream title2;
223 title2 << "/TMP/R_thrust_" << type[itype] << "_UC_" << ibin+1;
224 book(htemp,title2.str());
225 divide(_h_thrust[itype][1][ibin],
226 _h_thrust[itype][2][ibin],htemp);
227 asym = calcAsymmetry(htemp);
228 h3_thrust_UC->bin(ibin+1).set(asym.first, asym.second);
229 h2_thrust_UC->bin(ibin+1).set(asym.first, asym.second);
230 // hadron dirn
231 // opposite/like sign
232 std::ostringstream title3;
233 title3 << "/TMP/R_hadron_" << type[itype] << "_UL_" << ibin+1;
234 book(htemp,title3.str());
235 divide(_h_hadron[itype][1][ibin],
236 _h_hadron[itype][0][ibin],htemp);
237 asym = calcAsymmetry(htemp,2.);
238 h3_hadron_UL->bin(ibin+1).set(asym.first, asym.second);
239 h2_hadron_UL->bin(ibin+1).set(asym.first, asym.second);
240 // opposite/all sign
241 std::ostringstream title4;
242 title4 << "/TMP/R_hadron_" << type[itype] << "_UC_" << ibin+1;
243 book(htemp,title4.str());
244 divide(_h_hadron[itype][1][ibin],
245 _h_hadron[itype][2][ibin],htemp);
246 asym = calcAsymmetry(htemp,2.);
247 h3_hadron_UC->bin(ibin+1).set(asym.first, asym.second);
248 h2_hadron_UC->bin(ibin+1).set(asym.first, asym.second);
249 }
250 }
251 }
252
253 /// @}
254
255
256 /// @name Histograms
257 /// @{
258 Histo1DPtr _h_thrust[3][3][16],_h_hadron[3][3][16];
259 /// @}
260
261
262 };
263
264
265 RIVET_DECLARE_PLUGIN(BABAR_2015_I1377201);
266
267}
|