rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

OPAL_1997_I440103

Polarization of $\phi$, $D^{*+}$ and $B^*$ mesons at LEP1
Experiment: OPAL (LEP)
Inspire ID: 440103
Status: VALIDATED
Authors:
  • Peter Richardson
References:
  • Z.Phys. C74 (1997) 437-449
Beams: e- e+
Beam energies: ANY
Run details:
  • e+e- to hadrons

Mewasurement of the polarization of $\phi$, $D^{*+}$ and $B^*$ mesons at LEP1

Source code: OPAL_1997_I440103.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/Beam.hh"
  4#include "Rivet/Projections/FinalState.hh"
  5#include "Rivet/Projections/ChargedFinalState.hh"
  6#include "Rivet/Projections/UnstableParticles.hh"
  7#include "Rivet/Projections/Thrust.hh"
  8
  9#define I_KNOW_THE_INITIAL_QUARKS_PROJECTION_IS_DODGY_BUT_NEED_TO_USE_IT
 10#include "Rivet/Projections/InitialQuarks.hh"
 11
 12namespace Rivet {
 13
 14
 15  /// @brief Add a short analysis description here
 16  class OPAL_1997_I440103 : public Analysis {
 17  public:
 18
 19    /// Constructor
 20    RIVET_DEFAULT_ANALYSIS_CTOR(OPAL_1997_I440103);
 21
 22
 23    /// @name Analysis methods
 24    /// @{
 25
 26    /// Book histograms and initialise projections before the run
 27    void init() {
 28      // Initialise and register projections
 29      declare(Beam(), "Beams");
 30      declare(Thrust(FinalState()), "Thrust");
 31      declare(ChargedFinalState(), "FS");
 32      declare(InitialQuarks(), "IQF");
 33      declare(UnstableParticles(), "UFS" );
 34
 35      // Book histograms
 36      // B*
 37      book(_h_B , 8,1,1);
 38      book(_h_B2, "/TMP/c_theta_B", 20, -1.,1.);
 39      // phi
 40      book(_h_phi_ctheta , 5,1,1);
 41      book(_h_phi_ctheta2, "/TMP/c_theta_phi2", 20, -1.,1.   );
 42      book(_h_phi_ctheta3, "/TMP/c_theta_phi3", 20, -1.,1.   );
 43      book(_h_phi_ctheta4, "/TMP/c_theta_phi4", 20, -1.,1.   );
 44      book(_h_phi_alpha  , 5,1,2);
 45      book(_h_phi_alpha2 , "/TMP/alpha_phi2", 20, 0.,0.5*M_PI);
 46      book(_h_phi_alpha3 , "/TMP/alpha_phi3", 20, 0.,0.5*M_PI);
 47      book(_h_phi_alpha4 , "/TMP/alpha_phi4", 20, 0.,0.5*M_PI);
 48      book(_h_phi_beta   , 5,1,3);
 49      book(_h_phi_beta2  , "/TMP/beta_phi2", 20, 0.,0.5*M_PI );
 50      book(_h_phi_beta3  , "/TMP/beta_phi3", 20, 0.,0.5*M_PI );
 51      book(_h_phi_beta4  , "/TMP/beta_phi4", 20, 0.,0.5*M_PI );
 52      book(_c_phi_cos_plus , "/TMP/c_phi_cos_plus1");
 53      book(_c_phi_cos_neg  , "/TMP/c_phi_cos_neg1" );
 54      book(_c_phi_sin_plus , "/TMP/c_phi_sin_plus1");
 55      book(_c_phi_sin_neg  , "/TMP/c_phi_sin_neg1" );
 56      book(_c_phi_cos_plus2, "/TMP/c_phi_cos_plus2");
 57      book(_c_phi_cos_neg2 , "/TMP/c_phi_cos_neg2" );
 58      book(_c_phi_sin_plus2, "/TMP/c_phi_sin_plus2");
 59      book(_c_phi_sin_neg2 , "/TMP/c_phi_sin_neg2" );
 60      book(_c_phi_cos_plus3, "/TMP/c_phi_cos_plus3");
 61      book(_c_phi_cos_neg3 , "/TMP/c_phi_cos_neg3" );
 62      book(_c_phi_sin_plus3, "/TMP/c_phi_sin_plus3");
 63      book(_c_phi_sin_neg3 , "/TMP/c_phi_sin_neg3" );
 64      // D*
 65      book(_h_DS_ctheta , 6,1,1);
 66      book(_h_DS_ctheta2, "/TMP/c_theta_DS2", 20, -1.,1.   );
 67      book(_h_DS_alpha  , 7,1,1);
 68      book(_h_DS_alpha2 , "/TMP/alpha_DS2", 20, 0.,0.5*M_PI);
 69    }
 70
 71
 72    /// Perform the per-event analysis
 73    void analyze(const Event& event) {
 74      // First, veto on leptonic events by requiring at least 4 charged FS particles
 75      const FinalState& fs = apply<FinalState>(event, "FS");
 76      const size_t numParticles = fs.particles().size();
 77
 78      // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
 79      if (numParticles < 2) {
 80        MSG_DEBUG("Failed leptonic event cut");
 81        vetoEvent;
 82      }
 83      MSG_DEBUG("Passed leptonic event cut");
 84
 85      // Get beams and average beam momentum
 86      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 87      const double meanBeamMom = ( beams.first.p3().mod() +
 88                                   beams.second.p3().mod() ) / 2.0;
 89      MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
 90      Vector3 axis;
 91      if(beams.first.pid()>0)
 92        axis = beams.first.momentum().p3().unit();
 93      else
 94        axis = beams.second.momentum().p3().unit();
 95      // thrust, to define an axis
 96      const Thrust& thrust = apply<Thrust>(event, "Thrust");
 97
 98      int flavour = 0;
 99      const InitialQuarks& iqf = apply<InitialQuarks>(event, "IQF");
100
101      // If we only have two quarks (qqbar), just take the flavour.
102      // If we have more than two quarks, look for the highest energetic q-qbar pair.
103      /// @todo Yuck... does this *really* have to be quark-based?!?
104      if (iqf.particles().size() == 2) {
105        flavour = iqf.particles().front().abspid();
106      } else {
107        map<int, double> quarkmap;
108        for (const Particle& p : iqf.particles()) {
109          if (quarkmap[p.pid()] < p.E()) {
110            quarkmap[p.pid()] = p.E();
111          }
112        }
113        double maxenergy = 0.;
114        for (int i = 1; i <= 5; ++i) {
115          if (quarkmap[i]+quarkmap[-i] > maxenergy) {
116            flavour = i;
117          }
118        }
119      }
120
121      // loop over the particles
122      for (const Particle& p : apply<UnstableParticles>(event, "UFS").particles(Cuts::abspid==513 or Cuts::abspid==523 or
123                                                                                Cuts::pid==333 or Cuts::abspid==413)) {
124        int sign = p.pid()/p.abspid();
125        Particle decay;
126        if(p.children().size()!=2) continue;
127        // B*
128        if(p.abspid()==513 or p.abspid()==523) {
129          int mid = p.abspid()-2;
130          if(p.children()[0].pid()==sign*mid &&
131             p.children()[1].pid()==22) {
132            decay = p.children()[1];
133          }
134          else if(p.children()[1].pid()==sign*mid &&
135            p.children()[0].pid()==22) {
136            decay = p.children()[0];
137          }
138          else {
139            continue;
140          }
141        }
142        // phi
143        else if(p.pid()==333) {
144          // cut x_E > 0.7
145          double xE = p.momentum().E()/meanBeamMom;
146          if(xE<0.7) continue;
147          if(p.children()[0].pid()== 321 &&
148             p.children()[1].pid()==-321) {
149            decay = p.children()[0];
150          }
151          else if(p.children()[1].pid()== 321 &&
152            p.children()[0].pid()==-321) {
153            decay = p.children()[1];
154          }
155          else {
156            continue;
157          }
158        }
159        // D*
160        else if(p.abspid()==413) {
161          double xE = p.momentum().E()/meanBeamMom;
162          if(xE<0.5 || flavour!=4) continue;
163          if(p.children()[0].pid()==sign*421 &&
164             p.children()[1].pid()==sign*211) {
165            decay = p.children()[1];
166          }
167          else if(p.children()[1].pid()==sign*421 &&
168            p.children()[0].pid()==sign*211) {
169            decay = p.children()[0];
170          }
171          else {
172            continue;
173          }
174        }
175        LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
176        Vector3 e1z = p.p3().unit();
177        FourMomentum pp = boost.transform(decay.momentum());
178        Vector3 axis1 = boost.transform(decay.momentum()).p3().unit();
179        double ctheta = e1z.dot(axis1);
180        if(p.abspid()==513 or p.abspid()==523) {
181          _h_B ->fill(ctheta);
182          _h_B2->fill(ctheta);
183        }
184        // D*
185        else if(p.abspid()==413) {
186          // y and z axis
187          Vector3 e1y = e1z.cross(axis).unit();
188          Vector3 e1x = e1y.cross(e1z).unit();
189          // helicity beam axis, all phis
190          // cos theta_H
191          _h_DS_ctheta ->fill(ctheta);
192          _h_DS_ctheta2->fill(ctheta);
193          // alpha
194          double phi = atan2(e1y.dot(axis1),e1x.dot(axis1));
195          double alpha = abs(abs(phi)-0.5*M_PI);
196          _h_DS_alpha   ->fill(alpha);
197          _h_DS_alpha2  ->fill(alpha);
198        }
199        else if(p.pid()==333) {
200          // y and z axis
201          Vector3 e1y = e1z.cross(axis).unit();
202          Vector3 e1x = e1y.cross(e1z).unit();
203          // helicity beam axis, all phis
204          // cos theta_H
205          _h_phi_ctheta->fill(abs(ctheta));
206          _h_phi_ctheta2->fill(ctheta);
207          // alpha and beta
208          double phi = atan2(e1y.dot(axis1),e1x.dot(axis1));
209          double alpha = abs(abs(phi)-0.5*M_PI);
210          double beta =  abs(abs(phi+0.25*M_PI)-0.5*M_PI);
211          _h_phi_alpha   ->fill(alpha);
212          _h_phi_alpha2  ->fill(alpha);
213          _h_phi_beta    ->fill( beta);
214          _h_phi_beta2   ->fill( beta);
215          /// counters for asymmetries
216          double sin2H = 2.*ctheta*sqrt(1.-sqr(ctheta));
217          if(sin2H*cos(phi)>0.)
218            _c_phi_cos_plus->fill();
219          else
220            _c_phi_cos_neg->fill();
221          if(sin2H*sin(phi)>0.)
222            _c_phi_sin_plus->fill();
223          else
224            _c_phi_sin_neg->fill();
225          // whether or not is a primary hadron
226          Particle parent = p.parents()[0];
227          if(parent.children().size()==1 && parent.abspid()==p.abspid())
228            parent = parent.parents()[0];
229          bool primary = !PID::isHadron(parent.pid());
230          if (primary) {
231            // cos theta_H
232            _h_phi_ctheta3->fill(ctheta);
233            // alpha and beta
234            _h_phi_alpha3  ->fill(alpha);
235            _h_phi_beta3   ->fill( beta);
236            /// counters for asymmetries
237            if(sin2H*cos(phi)>0.)
238              _c_phi_cos_plus2->fill();
239            else
240              _c_phi_cos_neg2->fill();
241            if(sin2H*sin(phi)>0.)
242              _c_phi_sin_plus2->fill();
243            else
244              _c_phi_sin_neg2->fill();
245          }
246          // pT w.r.t thrust axis
247          double pT = sqrt(sqr(thrust.thrustMajorAxis().dot(p.momentum().p3()))+ sqr(thrust.thrustMinorAxis().dot(p.momentum().p3())));
248          // helicity-quark frame
249          if(pT>1.2) {
250            // cos theta H
251            _h_phi_ctheta4->fill(ctheta);
252            Vector3 axis2;
253            if(p.momentum().p3().dot(thrust.thrustAxis())>=0.) {
254              axis2 = thrust.thrustAxis();
255            }
256            else {
257              axis2 =-thrust.thrustAxis();
258            }
259            Vector3 e2y = e1z.cross(axis2).unit();
260            Vector3 e2x = e2y.cross(e1z).unit();
261            // alpha and beta
262            double phi = atan2(e2y.dot(axis1),e2x.dot(axis1));
263            double alpha = abs(abs(phi)-0.5*M_PI);
264            double beta =  abs(abs(phi+0.25*M_PI)-0.5*M_PI);
265            _h_phi_alpha4  ->fill(alpha);
266            _h_phi_beta4   ->fill( beta);
267            /// counters for asymmetries
268            double sin2H = 2.*ctheta*sqrt(1.-sqr(ctheta));
269            if(sin2H*cos(phi)>0.)
270              _c_phi_cos_plus3->fill();
271            else
272              _c_phi_cos_neg3->fill();
273            if(sin2H*sin(phi)>0.)
274              _c_phi_sin_plus3->fill();
275            else
276              _c_phi_sin_neg3->fill();
277
278          }
279        }
280      }
281    }
282
283    pair<double,double> calcRho(Histo1DPtr hist,unsigned int mode) {
284      if(hist->numEntries()==0.) return make_pair(0.,0.);
285      double sum1(0.),sum2(0.);
286      for (const auto& bin : hist->bins() ) {
287        double Oi = bin.sumW();
288        if (Oi==0.) continue;
289        double ai(0.),bi(0.);
290        if(mode==0) {
291          ai = 0.25*( -bin.xMin()*(3.-sqr(bin.xMin())) + bin.xMax()*(3.-sqr(bin.xMax())));
292          bi =-0.75*( -bin.xMin()*(1.-sqr(bin.xMin())) + bin.xMax()*(1.-sqr(bin.xMax())));
293        }
294        else if(mode==1) {
295          ai = 0.125*( -bin.xMin()*(3.+sqr(bin.xMin())) + bin.xMax()*(3.+sqr(bin.xMax())));
296          bi = 0.375*( -bin.xMin()*(1.-sqr(bin.xMin())) + bin.xMax()*(1.-sqr(bin.xMax())));
297        }
298        else if(mode==2) {
299          ai = -2.*(bin.xMin()-bin.xMax())/M_PI;
300          bi = -2.*(sin(2.*bin.xMin())-sin(2.*bin.xMax()))/M_PI;
301        }
302        double Ei = bin.errW();
303        sum1 += sqr(bi/Ei);
304        sum2 += bi/sqr(Ei)*(Oi-ai);
305      }
306      return make_pair(sum2/sum1,sqrt(1./sum1));
307    }
308
309    /// Normalise histograms etc., after the run
310    void finalize() {
311      // B*
312      normalize(_h_B,1.,false);
313      normalize(_h_B2);
314      pair<double,double> rho = calcRho(_h_B2,1);
315      Estimate1DPtr h_rhoB;
316      book(h_rhoB,4,1,1);
317      h_rhoB->bin(1).set(rho.first, rho.second);
318      // D*
319      normalize(_h_DS_ctheta );
320      normalize(_h_DS_ctheta2);
321      rho = calcRho(_h_DS_ctheta2,1);
322      Estimate1DPtr h_rhoD;
323      book(h_rhoD,3,1,1);
324      h_rhoD->bin(1).set(rho.first, rho.second);
325      normalize(_h_DS_alpha );
326      normalize(_h_DS_alpha2);
327      Estimate1DPtr h_reRho_D;
328      book(h_reRho_D,3,1,2);
329      rho = calcRho(_h_DS_alpha2,2);
330      h_reRho_D->bin(1).set(rho.first, rho.second);
331      // phi
332      // rho00
333      normalize(_h_phi_ctheta );
334      normalize(_h_phi_ctheta2);
335      normalize(_h_phi_ctheta3);
336      normalize(_h_phi_ctheta4);
337      Estimate1DPtr hrho_phi;
338      book(hrho_phi,1,1,1);
339      rho = calcRho(_h_phi_ctheta2,0);
340      hrho_phi->bin(1).set(rho.first, rho.second);
341      rho = calcRho(_h_phi_ctheta3,0);
342      hrho_phi->bin(2).set(rho.first, rho.second);
343      rho = calcRho(_h_phi_ctheta4,0);
344      hrho_phi->bin(3).set(rho.first, rho.second);
345      // Re rho
346      normalize(_h_phi_alpha );
347      normalize(_h_phi_alpha2);
348      normalize(_h_phi_alpha3);
349      normalize(_h_phi_alpha4);
350      Estimate1DPtr  hreRho_phi;
351      book(hreRho_phi,1,1,2);
352      rho = calcRho(_h_phi_alpha2,2);
353      hreRho_phi->bin(1).set(rho.first, rho.second);
354      rho = calcRho(_h_phi_alpha3,2);
355      hreRho_phi->bin(2).set(rho.first, rho.second);
356      rho = calcRho(_h_phi_alpha4,2);
357      hreRho_phi->bin(3).set(rho.first, rho.second);
358      // Im rho
359      normalize(_h_phi_beta );
360      normalize(_h_phi_beta2);
361      normalize(_h_phi_beta3);
362      normalize(_h_phi_beta4);
363      Estimate1DPtr himRho_phi;
364      book(himRho_phi,1,1,3);
365      rho = calcRho(_h_phi_beta2,2);
366      himRho_phi->bin(1).set(rho.first, rho.second);
367      rho = calcRho(_h_phi_beta3,2);
368      himRho_phi->bin(2).set(rho.first, rho.second);
369      rho = calcRho(_h_phi_beta4,2);
370      himRho_phi->bin(3).set(rho.first, rho.second);
371      // real diff
372      Estimate0D temp  = ((*_c_phi_cos_plus-*_c_phi_cos_neg)/(*_c_phi_cos_plus+*_c_phi_cos_neg));
373      Estimate0D temp2 = ((*_c_phi_cos_plus2-*_c_phi_cos_neg2)/(*_c_phi_cos_plus2+*_c_phi_cos_neg2));
374      Estimate0D temp3 = ((*_c_phi_cos_plus3-*_c_phi_cos_neg3)/(*_c_phi_cos_plus3+*_c_phi_cos_neg3));
375      Estimate1DPtr hreDiff_phi;
376      book(hreDiff_phi,1,1,4);
377      hreDiff_phi->bin(1) = temp;
378      hreDiff_phi->bin(2) = temp2;
379      hreDiff_phi->bin(3) = temp3;
380      // im diff
381      temp  = ((*_c_phi_sin_plus-*_c_phi_sin_neg)/(*_c_phi_sin_plus+*_c_phi_sin_neg));
382      temp2 = ((*_c_phi_sin_plus2-*_c_phi_sin_neg2)/(*_c_phi_sin_plus2+*_c_phi_sin_neg2));
383      temp3 = ((*_c_phi_sin_plus3-*_c_phi_sin_neg3)/(*_c_phi_sin_plus3+*_c_phi_sin_neg3));
384      Estimate1DPtr himDiff_phi;
385      book(himDiff_phi,1,1,5);
386      himDiff_phi->bin(1) = temp;
387      himDiff_phi->bin(2) = temp2;
388      himDiff_phi->bin(3) = temp3;
389    }
390
391    /// @}
392
393
394    /// @name Histograms
395    /// @{
396    Histo1DPtr _h_B,_h_B2;
397    Histo1DPtr _h_phi_ctheta, _h_phi_ctheta2, _h_phi_ctheta3, _h_phi_ctheta4;
398    Histo1DPtr _h_phi_alpha , _h_phi_alpha2 , _h_phi_alpha3 , _h_phi_alpha4 ;
399    Histo1DPtr _h_phi_beta  , _h_phi_beta2  , _h_phi_beta3  , _h_phi_beta4  ;
400    CounterPtr _c_phi_cos_plus, _c_phi_cos_neg, _c_phi_cos_plus2, _c_phi_cos_neg2, _c_phi_cos_plus3, _c_phi_cos_neg3;
401    CounterPtr _c_phi_sin_plus, _c_phi_sin_neg, _c_phi_sin_plus2, _c_phi_sin_neg2, _c_phi_sin_plus3, _c_phi_sin_neg3;
402    Histo1DPtr _h_DS_ctheta, _h_DS_ctheta2;
403    Histo1DPtr _h_DS_alpha , _h_DS_alpha2 ;
404    /// @}
405
406
407  };
408
409
410  RIVET_DECLARE_PLUGIN(OPAL_1997_I440103);
411
412
413}