Rivet analyses referenceOPAL_2000_I502750Polarization of $\rho^\pm$ and $\omega^0$ mesons at LEP1Experiment: OPAL (LEP) Inspire ID: 502750 Status: VALIDATED Authors:
Beam energies: (45.6, 45.6) GeV Run details:
The measurement of the polarization of $\rho^\pm$ and $\omega^0$ mesons at LEP1 by the OPAL experiment. Source code: OPAL_2000_I502750.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/Beam.hh"
4#include "Rivet/Projections/ChargedFinalState.hh"
5#include "Rivet/Projections/UnstableParticles.hh"
6#include "Rivet/Tools/BinnedHistogram.hh"
7
8namespace Rivet {
9
10
11 /// @brief rho+/- and omega polarization
12 class OPAL_2000_I502750 : public Analysis {
13 public:
14
15 /// Constructor
16 RIVET_DEFAULT_ANALYSIS_CTOR(OPAL_2000_I502750);
17
18
19 /// @name Analysis methods
20 //@{
21
22 /// Book histograms and initialise projections before the run
23 void init() {
24
25 // Initialise and register projections
26 declare(Beam(), "Beams");
27 declare(ChargedFinalState(), "FS");
28 declare(UnstableParticles(), "UFS");
29
30 // Book histograms
31 {Histo1DPtr temp; _h_ctheta_rho .add(0.025,0.05,book(temp, "ctheta_rho_0",20,-1.,1.));}
32 {Histo1DPtr temp; _h_ctheta_rho .add(0.05 ,0.1 ,book(temp, "ctheta_rho_1",20,-1.,1.));}
33 {Histo1DPtr temp; _h_ctheta_rho .add(0.1 ,0.15,book(temp, "ctheta_rho_2",20,-1.,1.));}
34 {Histo1DPtr temp; _h_ctheta_rho .add(0.15 ,0.3 ,book(temp, "ctheta_rho_3",20,-1.,1.));}
35 {Histo1DPtr temp; _h_ctheta_rho .add(0.3 ,0.6 ,book(temp, "ctheta_rho_4",20,-1.,1.));}
36 {Histo1DPtr temp; _h_ctheta_omega.add(0.025,0.05,book(temp, "ctheta_omega_0",20,-1.,1.));}
37 {Histo1DPtr temp; _h_ctheta_omega.add(0.05 ,0.1 ,book(temp, "ctheta_omega_1",20,-1.,1.));}
38 {Histo1DPtr temp; _h_ctheta_omega.add(0.1 ,0.15,book(temp, "ctheta_omega_2",20,-1.,1.));}
39 {Histo1DPtr temp; _h_ctheta_omega.add(0.15 ,0.3 ,book(temp, "ctheta_omega_3",20,-1.,1.));}
40 {Histo1DPtr temp; _h_ctheta_omega.add(0.3 ,0.6 ,book(temp, "ctheta_omega_4",20,-1.,1.));}
41 book(_h_ctheta_omega_all, "ctheta_omega_all",20,-1.,1.);
42 }
43
44 pair<double,double> calcRho(Histo1DPtr hist) {
45 if(hist->numEntries()==0.) return make_pair(0.,0.);
46 double sum1(0.),sum2(0.);
47 for (auto bin : hist->bins() ) {
48 double Oi = bin.area();
49 if(Oi==0.) continue;
50 double ai = 0.25*(bin.xMax()*(3.-sqr(bin.xMax())) - bin.xMin()*(3.-sqr(bin.xMin())));
51 double bi = 0.75*(bin.xMin()*(1.-sqr(bin.xMin())) - bin.xMax()*(1.-sqr(bin.xMax())));
52 double Ei = bin.areaErr();
53 sum1 += sqr(bi/Ei);
54 sum2 += bi/sqr(Ei)*(Oi-ai);
55 }
56 return make_pair(sum2/sum1,sqrt(1./sum1));
57 }
58
59 bool findOmegaDecay(Particle omega,Particles & pi0, Particles & pip, Particles & pim) {
60 for(const Particle & child : omega.children()) {
61 if(child.pid()==211)
62 pip.push_back(child);
63 else if(child.pid()==-211)
64 pim.push_back(child);
65 else if(child.pid()==111)
66 pi0.push_back(child);
67 else if(!child.children().empty()) {
68 if(!findOmegaDecay(child,pi0,pip,pim)) return false;
69 }
70 else
71 return false;
72 }
73 return true;
74 }
75
76 /// Perform the per-event analysis
77 void analyze(const Event& event) {
78 // First, veto on leptonic events by requiring at least 4 charged FS particles
79 const FinalState& fs = apply<FinalState>(event, "FS");
80 const size_t numParticles = fs.particles().size();
81
82 // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
83 if (numParticles < 2) {
84 MSG_DEBUG("Failed leptonic event cut");
85 vetoEvent;
86 }
87 MSG_DEBUG("Passed leptonic event cut");
88 // Get beams and average beam momentum
89 const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
90 const double meanBeamMom = ( beams.first.p3().mod() +
91 beams.second.p3().mod() ) / 2.0;
92 MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
93 // loop over rho and omega mesons
94 const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
95 for (const Particle& p : ufs.particles(Cuts::abspid==213 || Cuts::abspid==223)) {
96 double xE = p.momentum().t()/meanBeamMom;
97 Vector3 e1z = p.momentum().p3().unit();
98 LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
99 if(p.abspid()==213) {
100 if(p.children().size()!=2) continue;
101 int sign = p.pid()/213;
102 Particle pion;
103 if(p.children()[0].pid()==sign*211 && p.children()[1].pid()==111) {
104 pion = p.children()[0];
105 }
106 else if(p.children()[1].pid()==sign*211 && p.children()[0].pid()==111) {
107 pion = p.children()[1];
108 }
109 else
110 continue;
111 Vector3 axis1 = boost.transform(pion.momentum()).p3().unit();
112 double ctheta = e1z.dot(axis1);
113 _h_ctheta_rho.fill(xE,ctheta);
114 }
115 else {
116 Particles pi0,pip,pim;
117 bool three_pi = findOmegaDecay(p,pi0,pip,pim);
118 if(!three_pi || pi0.size()!=1 || pip.size()!=1 || pim.size()!=1)
119 continue;
120 Vector3 v1 = boost.transform(pi0[0].momentum()).p3().unit();
121 Vector3 v2 = boost.transform(pip[0].momentum()).p3().unit();
122 Vector3 norm = v1.cross(v2).unit();
123 double ctheta = e1z.dot(norm);
124 _h_ctheta_omega.fill(xE,ctheta);
125 if(xE>0.025) _h_ctheta_omega_all->fill(ctheta);
126 }
127 }
128 }
129
130
131 /// Normalise histograms etc., after the run
132 void finalize() {
133 vector<double> x = {0.025,0.05,0.1,0.15,0.3,0.6};
134 Scatter2DPtr h_rho ;
135 book(h_rho, 1,1,1);
136 Scatter2DPtr h_omega;
137 book(h_omega, 2,1,1);
138 for(unsigned int ix=0;ix<_h_ctheta_rho.histos().size();++ix) {
139 // rho
140 normalize(_h_ctheta_rho.histos()[ix]);
141 pair<double,double> rho00 = calcRho(_h_ctheta_rho.histos()[ix]);
142 h_rho->addPoint(0.5*(x[ix]+x[ix+1]), rho00.first, make_pair(0.5*(x[ix+1]-x[ix]),0.5*(x[ix+1]-x[ix])),
143 make_pair(rho00.second,rho00.second) );
144 // omega
145 normalize(_h_ctheta_omega.histos()[ix]);
146 rho00 = calcRho(_h_ctheta_omega.histos()[ix]);
147 h_omega->addPoint(0.5*(x[ix]+x[ix+1]), rho00.first, make_pair(0.5*(x[ix+1]-x[ix]),0.5*(x[ix+1]-x[ix])),
148 make_pair(rho00.second,rho00.second) );
149 }
150 // omega over whole range
151 Scatter2DPtr h_omega_all;
152 book(h_omega_all,2,2,1);
153 normalize(_h_ctheta_omega_all);
154 pair<double,double> rho00 = calcRho(_h_ctheta_omega_all);
155 h_omega_all->addPoint(0.5125, rho00.first, make_pair(0.4875,0.4875),
156 make_pair(rho00.second,rho00.second) );
157 }
158
159 //@}
160
161
162 /// @name Histograms
163 //@{
164 BinnedHistogram _h_ctheta_rho,_h_ctheta_omega;
165 Histo1DPtr _h_ctheta_omega_all;
166 //@}
167
168
169 };
170
171
172 // The hook for the plugin system
173 RIVET_DECLARE_PLUGIN(OPAL_2000_I502750);
174
175
176}
|