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 ALEPH_1996_I402895 : public Analysis {
10  public:
11
12    /// Constructor
13    RIVET_DEFAULT_ANALYSIS_CTOR(ALEPH_1996_I402895);
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
29    }
30
31    void findDecayProducts(Particle p, Particles & lep, Particles & nu) {
32      for(const Particle & child : p.children()) {
33	if(PID::isHadron(child.pid())) continue;
34	if(child.abspid()==11 or child.abspid()==13)
35	  lep.push_back(child);
36	else if(child.abspid()==12 or child.abspid()==14)
37	  nu.push_back(child);
38	else if(child.abspid()!=15)
39	  findDecayProducts(child,lep,nu);
40      }
41    }
42
43    /// Perform the per-event analysis
44    void analyze(const Event& event) {
45      const FinalState& ufs = apply<UnstableParticles>(event, "UFS");
46      // loop over weakly decaying b-baryons
47      for (const Particle& p : ufs.particles(Cuts::abspid==5122 || Cuts::abspid==5132 ||
48					     Cuts::abspid==5232 || Cuts::abspid==5332)) {
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      }
55    }
56
57
58    /// Normalise histograms etc., after the run
59    void finalize() {
60      normalize(_h_El   );
61      normalize(_h_Ev   );
62      if(_h_El->effNumEntries()!=0. and _h_Ev->effNumEntries()!=0.) {
63	double Ev  = _h_Ev->xMean();
64	double El  = _h_El->xMean();
65	double dEv = _h_Ev->xStdErr();
66	double dEl = _h_El->xStdErr();
67	double ratio = El/Ev;
68	double dr    = (Ev*dEl-El*dEv)/sqr(Ev);
69	double rho = 0.091;
70	double P = (7. + rho*(30. - 40.*ratio) + 4.*(2.-3.*ratio)*ratio)/sqr(1.+2.*ratio);
71	double dP = (20.*(-1. + 4.*rho*(-2. + ratio) - 2.*ratio))/pow(1. + 2.*ratio,3)*dr;
72       	Scatter2DPtr h_pol;
73	book(h_pol,1,1,1);
74       	h_pol->addPoint(91.2, P, make_pair(0.5,0.5),make_pair(dP,dP) );
75      }
76    }
77
78    //@}
79
80
81    /// @name Histograms
82    //@{
83    Histo1DPtr _h_El, _h_Ev;
84    //@}
85
86
87  };
88
89
90  // The hook for the plugin system
91  RIVET_DECLARE_PLUGIN(ALEPH_1996_I402895);
92
93
94}