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 }
150
151 pair<double,double> calcAsymmetry(Scatter2DPtr hist,double fact=1.) {
152 double sum1(0.),sum2(0.);
153 for (auto point : hist->points() ) {
154 double Oi = point.y();
155 if(Oi==0. || std::isnan(Oi) ) continue;
156 double ai = 1.;
157 double bi = (sin(fact*point.xMax())-sin(fact*point.xMin()))/(point.xMax()-point.xMin())/fact;
158 double Ei = point.yErrAvg();
159 sum1 += sqr(bi/Ei);
160 sum2 += bi/sqr(Ei)*(Oi-ai);
161 }
162 if(sum1==0.) return make_pair(0.,0.);
163 return make_pair(sum2/sum1*1e4,sqrt(1./sum1)*1e4);
164 }
165
166 /// Normalise histograms etc., after the run
167 void finalize() {
168 for(unsigned int itype=0;itype<3;++itype) {
169 for(unsigned int icharge=0;icharge<3;++icharge) {
170 for(unsigned int ibin=0;ibin<16;++ibin) {
171 normalize(_h_thrust[itype][icharge][ibin]);
172 normalize(_h_hadron[itype][icharge][ibin]);
173 }
174 }
175 }
176 // construct ther ratios
177 // declare the histos for the distributions
178 string type [3] = {"pipi","Kpi","KK"};
179 string charge[3] = {"Like","Opposite","All"};
180 for(unsigned int itype=0;itype<3;++itype) {
181 Scatter3DPtr h3_thrust_UL;
182 book(h3_thrust_UL,2*itype+1,1,2,0);
183 Scatter3DPtr h3_thrust_UC;
184 book(h3_thrust_UC,2*itype+1,1,3,0);
185 Scatter3DPtr h3_hadron_UL;
186 book(h3_hadron_UL,2*itype+2,1,2,0);
187 Scatter3DPtr h3_hadron_UC;
188 book(h3_hadron_UC,2*itype+2,1,3,0);
189
190 unsigned int ihist=1;
191 Scatter2DPtr h2_thrust_UL;
192 book(h2_thrust_UL,7+2*itype,ihist,2);
193 Scatter2DPtr h2_thrust_UC;
194 book(h2_thrust_UC,7+2*itype,ihist,3);
195 Scatter2DPtr h2_hadron_UL;
196 book(h2_hadron_UL,8+2*itype,ihist,2);
197 Scatter2DPtr h2_hadron_UC;
198 book(h2_hadron_UC,8+2*itype,ihist,3);
199
200 Scatter3D temphisto1(refData<Scatter3D>(2*itype+1, 1, 2));
201 Scatter3D temphisto2(refData<Scatter3D>(2*itype+2, 1, 2));
202 for(unsigned int ibin=0;ibin<16;++ibin) {
203 const Point3D & p1 = temphisto1.points()[ibin];
204 const Point3D & p2 = temphisto2.points()[ibin];
205 if(ibin>0 && ibin%4==0) {
206 ++ihist;
207 book(h2_thrust_UL,7+2*itype,ihist,2);
208 book(h2_thrust_UC,7+2*itype,ihist,3);
209 book(h2_hadron_UL,8+2*itype,ihist,2);
210 book(h2_hadron_UC,8+2*itype,ihist,3);
211 }
212 // thrust direction
213 // opposite/like sign
214 std::ostringstream title1;
215 title1 << "/TMP/R_thrust_" << type[itype] << "_UL_" << ibin+1;
216 Scatter2DPtr htemp;
217 book(htemp,title1.str());
218 divide(_h_thrust[itype][1][ibin],
219 _h_thrust[itype][0][ibin],htemp);
220 pair<double,double> asym = calcAsymmetry(htemp);
221 h3_thrust_UL->addPoint(p1.x() ,p1.y() ,asym.first,
222 p1.xErrs(),p1.yErrs(),make_pair(asym.second,asym.second) );
223 h2_thrust_UL->addPoint(p1.y() ,asym.first,p1.yErrs(),make_pair(asym.second,asym.second) );
224 // opposite/all sign
225 std::ostringstream title2;
226 title2 << "/TMP/R_thrust_" << type[itype] << "_UC_" << ibin+1;
227 book(htemp,title2.str());
228 divide(_h_thrust[itype][1][ibin],
229 _h_thrust[itype][2][ibin],htemp);
230 asym = calcAsymmetry(htemp);
231 h3_thrust_UC->addPoint(p1.x() ,p1.y() ,asym.first,
232 p1.xErrs(),p1.yErrs(),make_pair(asym.second,asym.second) );
233 h2_thrust_UC->addPoint(p1.y() ,asym.first,p1.yErrs(),make_pair(asym.second,asym.second) );
234 // hadron dirn
235 // opposite/like sign
236 std::ostringstream title3;
237 title3 << "/TMP/R_hadron_" << type[itype] << "_UL_" << ibin+1;
238 book(htemp,title3.str());
239 divide(_h_hadron[itype][1][ibin],
240 _h_hadron[itype][0][ibin],htemp);
241 asym = calcAsymmetry(htemp,2.);
242 h3_hadron_UL->addPoint(p2.x() ,p2.y() ,asym.first,
243 p2.xErrs(),p2.yErrs(),make_pair(asym.second,asym.second) );
244 h2_hadron_UL->addPoint(p2.y() ,asym.first,p2.yErrs(),make_pair(asym.second,asym.second) );
245 // opposite/all sign
246 std::ostringstream title4;
247 title4 << "/TMP/R_hadron_" << type[itype] << "_UC_" << ibin+1;
248 book(htemp,title4.str());
249 divide(_h_hadron[itype][1][ibin],
250 _h_hadron[itype][2][ibin],htemp);
251 asym = calcAsymmetry(htemp,2.);
252 h3_hadron_UC->addPoint(p2.x() ,p2.y() ,asym.first,
253 p2.xErrs(),p2.yErrs(),make_pair(asym.second,asym.second) );
254 h2_hadron_UC->addPoint(p2.y() ,asym.first,p2.yErrs(),make_pair(asym.second,asym.second) );
255 }
256 }
257 }
258
259 //@}
260
261
262 /// @name Histograms
263 //@{
264 Histo1DPtr _h_thrust[3][3][16],_h_hadron[3][3][16];
265 //@}
266
267
268 };
269
270
271 RIVET_DECLARE_PLUGIN(BABAR_2015_I1377201);
272
273}
|