1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/UnstableParticles.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief Lambdab0 -> Lambda0 gamma
 10  class LHCB_2022_I1971920 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(LHCB_2022_I1971920);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22      // Initialise and register projections
 23      declare(UnstableParticles(Cuts::pid==5122), "UFS" );
 24      // histos
 25      book(_h_ctheta_gamma, "TMP/ctheta_gamma", 20, -1, 1);
 26    }
 27
 28
 29    /// Perform the per-event analysis
 30    void analyze(const Event& event) {
 31      // loop over Lambda_b baryons
 32      for (const Particle& lamB : apply<UnstableParticles>(event, "UFS").particles()) {
 33       	int sign = lamB.pid()/5122;
 34        // Lambda_b -> Lambda0 gamma
 35      	if (lamB.children().size()!=2) continue;
 36       	Particle lambda,gamma;
 37        if (lamB.children()[1].pid()==sign*3122 &&
 38           lamB.children()[0].pid()==22) {
 39          lambda = lamB.children()[1];
 40          gamma  = lamB.children()[0];
 41        }
 42        else if(lamB.children()[0].pid()==sign*3122 &&
 43          lamB.children()[1].pid()==22) {
 44          lambda = lamB.children()[0];
 45          gamma  = lamB.children()[1];
 46        }
 47        else continue;
 48        // Lambda0 -> p pi+
 49      	if (lambda.children().size()!=2) continue;
 50       	Particle proton, pion;
 51        if (lambda.children()[0].pid()== sign*2212 &&
 52           lambda.children()[1].pid()==-sign*211) {
 53          proton = lambda.children()[0];
 54          pion   = lambda.children()[1];
 55        }
 56        else if (lambda.children()[1].pid()== sign*2212 &&
 57          lambda.children()[0].pid()==-sign*211) {
 58          proton = lambda.children()[1];
 59          pion   = lambda.children()[0];
 60        }
 61        else  continue;
 62      	// first boost to the lamB rest frame
 63      	LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(lamB.momentum().betaVec());
 64      	FourMomentum plambda = boost1.transform(lambda.mom());
 65      	FourMomentum pproton = boost1.transform(proton.mom());
 66      	// to lambda rest frame
 67      	LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(plambda.betaVec());
 68      	Vector3 axis = plambda.p3().unit();
 69      	FourMomentum pp = boost2.transform(pproton);
 70      	// calculate angle
 71      	const double cTheta = pp.p3().unit().dot(axis);
 72        _h_ctheta_gamma->fill(cTheta);
 73      }
 74    }
 75
 76    pair<double,double> calcAlpha(Histo1DPtr hist) {
 77      if (hist->numEntries()==0.) return make_pair(0.,0.);
 78      double sum1(0.),sum2(0.);
 79      for (const auto& bin : hist->bins()) {
 80        double Oi = bin.sumW();
 81        if (Oi==0.) continue;
 82        double ai = 0.5*(bin.xMax()-bin.xMin());
 83        double bi = 0.5*ai*(bin.xMax()+bin.xMin());
 84        double Ei = bin.errW();
 85        sum1 += sqr(bi/Ei);
 86        sum2 += bi/sqr(Ei)*(Oi-ai);
 87      }
 88      return make_pair(sum2/sum1,sqrt(1./sum1));
 89    }
 90
 91    /// Normalise histograms etc., after the run
 92    void finalize() {
 93      normalize(_h_ctheta_gamma);
 94      pair<double,double> alpha = calcAlpha(_h_ctheta_gamma);
 95      const double aLam = 0.754;
 96      alpha.first  /= aLam;
 97      alpha.second /= aLam;
 98      Estimate0DPtr tmp;
 99      book(tmp,1,1,1);
100      tmp->set(alpha.first,alpha.second);
101    }
102
103    /// @}
104
105
106    /// @name Histograms
107    /// @{
108    Histo1DPtr _h_ctheta_gamma;
109    /// @}
110
111
112  };
113
114
115  RIVET_DECLARE_PLUGIN(LHCB_2022_I1971920);
116
117}