1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/UnstableParticles.hh"
 4
 5namespace Rivet {
 6
 7
 8  /// @brief b-baryon polarization
 9  class OPAL_1998_I474012 : public Analysis {
10  public:
11
12    /// Constructor
13    RIVET_DEFAULT_ANALYSIS_CTOR(OPAL_1998_I474012);
14
15
16    /// @name Analysis methods
17    //@{
18
19    /// Book histograms and initialise projections before the run
20    void init() {
21
22      // Initialise and register projections
23      declare(UnstableParticles(), "UFS");
24
25      // Book histograms
26      book(_h_El   , "El",   45,0.,45.0);
27      book(_h_Ev   , "Ev",   45,0.,45.0);
28      book(_h_ratio, "ratio",20,0.,10.0);
29
30    }
31
32    void findDecayProducts(Particle p, Particles & lep, Particles & nu) {
33      for(const Particle & child : p.children()) {
34	if(PID::isHadron(child.pid())) continue;
35	if(child.abspid()==11 or child.abspid()==13)
36	  lep.push_back(child);
37	else if(child.abspid()==12 or child.abspid()==14)
38	  nu.push_back(child);
39	else if(child.abspid()!=15)
40	  findDecayProducts(child,lep,nu);
41      }
42    }
43
44    /// Perform the per-event analysis
45    void analyze(const Event& event) {
46      const FinalState& ufs = apply<UnstableParticles>(event, "UFS");
47      // loop over weakly decaying b-baryons
48      for (const Particle& p : ufs.particles(Cuts::abspid==5122)) {
49	Particles lep,nu;
50	findDecayProducts(p,lep,nu);
51	if(lep.size()!=1 || nu.size()!=1) continue;
52	_h_El   ->fill(lep[0].momentum().t());
53	_h_Ev   ->fill(nu [0].momentum().t());
54	_h_ratio->fill(nu [0].momentum().t()/lep[0].momentum().t());
55      }
56    }
57
58
59    /// Normalise histograms etc., after the run
60    void finalize() {
61      normalize(_h_El   );
62      normalize(_h_Ev   );
63      normalize(_h_ratio);
64      if(_h_El->effNumEntries()!=0. and _h_Ev->effNumEntries()!=0.) {
65	double Ev  = _h_Ev->xMean();
66	double El  = _h_El->xMean();
67	double dEv = _h_Ev->xStdErr();
68	double dEl = _h_El->xStdErr();
69	double ratio = Ev/El;
70	double dr    = (El*dEv-Ev*dEl)/sqr(El);
71	double rho = 0.091;
72	double  P = 7. - 20./(2. + ratio) + 10.*rho*ratio*(-4. + 3.*ratio)/sqr(2.+ratio);
73	double dP = (20.*(2. + ratio + rho*(-4. + 8.*ratio)))/pow(2.+ratio,3)*dr;
74       	Scatter2DPtr h_pol;
75	book(h_pol, 1,1,1);
76       	h_pol->addPoint(91.2, P, make_pair(0.5,0.5),make_pair(dP,dP) );
77      }
78    }
79
80    //@}
81
82
83    /// @name Histograms
84    //@{
85    Histo1DPtr _h_El, _h_Ev, _h_ratio;
86    //@}
87
88
89  };
90
91
92  // The hook for the plugin system
93  RIVET_DECLARE_PLUGIN(OPAL_1998_I474012);
94
95
96}