Rivet analyses referenceBABAR_2014_I1254862Azimuthal asymmetries in inclusive charged $\pi\pi$ pair production at 10.58 GeVExperiment: BABAR (PEP-II) Inspire ID: 1254862 Status: VALIDATED NOHEPDATA Authors:
Beam energies: (5.3, 5.3) GeV Run details:
Measurement of azimuthal asymmetries in inclusive charged $\pi\pi$ pair production at $\sqrt{s}=10.58$ GeV by the BABAR experiment. Only the distributions in $z_{1,2}$ are currently implemented. Source code: BABAR_2014_I1254862.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#include "Rivet/Tools/Random.hh"
7
8namespace Rivet {
9
10
11 /// @brief azimuthal asymmetries in pipi
12 class BABAR_2014_I1254862 : public Analysis {
13 public:
14
15 /// Constructor
16 RIVET_DEFAULT_ANALYSIS_CTOR(BABAR_2014_I1254862);
17
18
19 /// @name Analysis methods
20 /// @{
21
22 /// Book histograms and initialise projections before the run
23 void init() {
24 // projections
25 const FinalState fs;
26 declare(fs,"FS");
27 declare(Thrust(fs),"Thrust");
28 declare(Beam(), "Beams");
29 // declare the histos for the distributions
30 string charge[3] = {"Like","Opposite","All"};
31 unsigned int nbin=20;
32 for (unsigned int icharge=0;icharge<3;++icharge) {
33 for (unsigned int ibin1=0;ibin1<6;++ibin1) {
34 for (unsigned int ibin2=0;ibin2<6;++ibin2) {
35 book(_h_thrust[icharge][ibin1][ibin2],
36 "TMP/h_thrust_"+charge[icharge]+"_" +toString(ibin1+1) + "_" + toString(ibin2+1),
37 nbin, 0., M_PI);
38 book(_h_hadron[icharge][ibin1][ibin2],
39 "TMP/h_hadron_"+charge[icharge]+"_" +toString(ibin1+1) + "_" + toString(ibin2+1),
40 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<.4) return 2;
50 else if(z<.5) return 3;
51 else if(z<.7) return 4;
52 else return 5;
53 }
54
55 /// Perform the per-event analysis
56 void analyze(const Event& event) {
57 // get the axis, direction of incoming electron
58 const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
59 Vector3 axis1;
60 if (beams.first.pid()>0) {
61 axis1 = beams.first .momentum().p3().unit();
62 }
63 else {
64 axis1 = beams.second.momentum().p3().unit();
65 }
66 // apply thrust cuts T > 0.8
67 Thrust thrust = apply<Thrust>(event,"Thrust");
68 if(thrust.thrust()<=0.8) vetoEvent;
69 // construct x,y,z axes for thrust defn
70 ThreeVector t_z = thrust.thrustAxis();
71 ThreeVector t_x = (axis1-t_z.dot(axis1)*t_z).unit();
72 ThreeVector t_y = t_z.cross(t_x);
73 // loop over the particles
74 Particles charged = apply<FinalState>(event,"FS").particles(Cuts::abspid==PID::PIPLUS);
75 for (unsigned int ix=0;ix<charged.size();++ix) {
76 // z and angle cut
77 const double x1=2.*charged[ix].momentum().t()/sqrtS();
78 if (x1<0.15||x1>.9) continue;
79 double dot1 = t_z.dot(charged[ix].p3().unit());
80 if (abs(dot1)<sqrt(0.5)) continue;
81 for (unsigned int iy=ix+1;iy<charged.size();++iy) {
82 const double x2=2.*charged[iy].momentum().t()/sqrtS();
83 // z and angle cut
84 if (x2<0.15||x2>.9) continue;
85 double dot2 = t_z.dot(charged[iy].p3().unit());
86 if (abs(dot2)<sqrt(.5) || dot1*dot2>0.) continue;
87 Particle p1=charged[ix], p2=charged[iy];
88 double z1(x1),z2(x2);
89 // randomly order the particles
90 if (rand01() < 0.5 ) {
91 swap(p1,p2);
92 swap(z1,z2);
93 }
94 // thrust def
95 double phi12 = atan2(p1.p3().dot(t_y),p1.p3().dot(t_x))+atan2(p2.p3().dot(t_y),p2.p3().dot(t_x));
96 if (phi12>M_PI) phi12 -= 2*M_PI;
97 if (phi12<-M_PI) phi12 += 2*M_PI;
98 if (phi12<0.) phi12 = -phi12;
99 // hadron defn
100 ThreeVector h_z = p2.p3().unit();
101 ThreeVector h_x = (axis1-h_z.dot(axis1)*h_z).unit();
102 ThreeVector pt1 = p1.p3()-h_z.dot(p1.p3())*h_z;
103 double phi0 = pt1.angle(h_x);
104 if (phi0>M_PI) phi0 -= 2*M_PI;
105 if (phi0<-M_PI) phi0 += 2*M_PI;
106 unsigned int ibin1=iBin(z1);
107 unsigned int ibin2=iBin(z2);
108 if (p1.pid()==p2.pid()) {
109 _h_thrust[0][ibin1][ibin2]->fill(phi12);
110 _h_hadron[0][ibin1][ibin2]->fill(phi0);
111 }
112 else {
113 _h_thrust[1][ibin1][ibin2]->fill(phi12);
114 _h_hadron[1][ibin1][ibin2]->fill(phi0);
115 }
116 _h_thrust[2][ibin1][ibin2]->fill(phi12);
117 _h_hadron[2][ibin1][ibin2]->fill(phi0);
118 }
119 }
120 }
121
122 pair<double,double> calcAsymmetry(Estimate1DPtr hist,double fact=1.) {
123 double sum1(0.),sum2(0.);
124 for (const auto& bin : hist->bins() ) {
125 double Oi = bin.val();
126 if (Oi==0. || std::isnan(Oi) ) continue;
127 double ai = 1.;
128 double bi = (sin(fact*bin.xMax())-sin(fact*bin.xMin()))/(bin.xMax()-bin.xMin())/fact;
129 double Ei = bin.errAvg();
130 sum1 += sqr(bi/Ei);
131 sum2 += bi/sqr(Ei)*(Oi-ai);
132 }
133 if (sum1==0.) return make_pair(0.,0.);
134 return make_pair(sum2/sum1*1e2,sqrt(1./sum1)*1e2);
135 }
136
137 /// Normalise histograms etc., after the run
138 void finalize() {
139 for (unsigned int ibin1=0;ibin1<6;++ibin1) {
140 Estimate1DPtr hthrustUL,hhadronUL;
141 book(hthrustUL,1,1+ibin1,1);
142 book(hhadronUL,1,1+ibin1,2);
143 Estimate1DPtr hthrustUC,hhadronUC;
144 book(hthrustUC,2,1+ibin1,1);
145 book(hhadronUC,2,1+ibin1,2);
146 for (unsigned int ibin2=0;ibin2<6;++ibin2) {
147 for (unsigned int icharge=0;icharge<3;++icharge) {
148 normalize(_h_thrust[icharge][ibin1][ibin2]);
149 normalize(_h_hadron[icharge][ibin1][ibin2]);
150 }
151 Estimate1DPtr htemp;
152 book(htemp,"TMP/R_thrust_UL_"+toString(ibin1)+"_"+toString(ibin2),_h_thrust[0][ibin1][ibin2]->xEdges());
153 // UL thrust
154 divide(_h_thrust[1][ibin1][ibin2],_h_thrust[0][ibin1][ibin2],htemp);
155 pair<double,double> asym = calcAsymmetry(htemp);
156 hthrustUL->bin(ibin2+1).set(asym.first,asym.second);
157 // UC thrust
158 book(htemp,"TMP/R_thrust_UC_"+toString(ibin1)+"_"+toString(ibin2),_h_thrust[1][ibin1][ibin2]->xEdges());
159 divide(_h_thrust[1][ibin1][ibin2],_h_thrust[2][ibin1][ibin2],htemp);
160 asym = calcAsymmetry(htemp);
161 hthrustUC->bin(ibin2+1).set(asym.first,asym.second);
162 // UL hadron
163 book(htemp,"TMP/R_hadron_UL_"+toString(ibin1)+"_"+toString(ibin2),_h_hadron[1][ibin1][ibin2]->xEdges());
164 divide(_h_hadron[1][ibin1][ibin2],_h_hadron[0][ibin1][ibin2],htemp);
165 asym = calcAsymmetry(htemp);
166 hhadronUL->bin(ibin2+1).set(asym.first,asym.second);
167 // UC hadron
168 book(htemp,"TMP/R_hadron_UC_"+toString(ibin1)+"_"+toString(ibin2),_h_hadron[1][ibin1][ibin2]->xEdges());
169 divide(_h_hadron[1][ibin1][ibin2],_h_hadron[2][ibin1][ibin2],htemp);
170 asym = calcAsymmetry(htemp);
171 hhadronUC->bin(ibin2+1).set(asym.first,asym.second);
172 }
173 }
174 }
175
176 /// @}
177
178
179 /// @name Histograms
180 /// @{
181 Histo1DPtr _h_thrust[3][6][6],_h_hadron[3][6][6];
182 /// @}
183
184
185 };
186
187
188 RIVET_DECLARE_PLUGIN(BABAR_2014_I1254862);
189
190}
|
