1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief Xi* decay asymmetries
  9  class BABAR_2008_I781294 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(BABAR_2008_I781294);
 14
 15
 16    /// @name Analysis methods
 17    /// @{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21      // Initialise and register projections
 22      declare(UnstableParticles(), "UFS" );
 23      // book histos
 24      book(_h_ctheta,1,1,1);
 25      book(_h_P0    ,2,1,1);
 26      book(_h_P2    ,2,1,2);
 27      book(_wgtSum,"/TMP/WSUM");
 28    }
 29
 30    void findChildren(const Particle &p, int & sign,
 31		      unsigned int & nprod, Particles & Xi, Particles &pi,Particles &K) {
 32      for(const Particle & child : p.children()) {
 33	if(child.pid()==sign*3312) {
 34	  Xi.push_back(child);
 35	  ++nprod;
 36	}
 37	else if(child.pid()==211) {
 38	  pi.push_back(child);
 39	  ++nprod;
 40	}
 41	else if(child.pid()==321) {
 42	  K.push_back(child);
 43	  ++nprod;
 44	}
 45	else if(child.children().empty()) {
 46	  ++nprod;
 47	}
 48	else {
 49	  findChildren(child,sign,nprod,Xi,pi,K);
 50	}
 51      }
 52    }
 53
 54    /// Perform the per-event analysis
 55    void analyze(const Event& event) {
 56      // loop over Lambda_c
 57      for(const Particle& baryon : apply<UnstableParticles>(event, "UFS").particles(Cuts::abspid==4122)) {
 58	int sign = baryon.pid()/baryon.abspid();
 59	Particles Xi,pi,K;
 60	unsigned int nprod(0);
 61	findChildren(baryon,sign,nprod,Xi,pi,K);
 62	// check Lambda_c -> Xi- pi+ K+ decya mode
 63	if(nprod!=3||Xi.size()!=1||pi.size()!=1||K.size()!=1) continue;
 64	// first boost to the Lambda_c rest frame
 65	LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(baryon.momentum().betaVec());
 66	FourMomentum pbaryon1 = Xi[0].momentum()+pi[0].momentum();
 67	pbaryon1              = boost1.transform(pbaryon1);
 68	FourMomentum pbaryon2 = boost1.transform(Xi[0].momentum());
 69	// to Xi* frame
 70	LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(pbaryon1.betaVec());
 71	Vector3 axis = pbaryon1.p3().unit();
 72	FourMomentum pp = boost2.transform(pbaryon2);
 73	// calculate angle
 74	double cTheta = pp.p3().unit().dot(axis);
 75	double mass = pbaryon1.mass();
 76	if(mass>1.5 && mass< 1.65)
 77	  _h_ctheta->fill(cTheta);
 78	_h_P0->fill(mass);
 79	_h_P2->fill(mass,0.5*(3.*sqr(cTheta)-1.));
 80	_wgtSum->fill();
 81      }
 82    }
 83
 84
 85    /// Normalise histograms etc., after the run
 86    void finalize() {
 87      normalize(_h_ctheta);
 88      normalize(_h_P0,1.,false);
 89      normalize(_h_P2,1.,false);
 90    }
 91    /// @}
 92
 93
 94    /// @name Histograms
 95    /// @{
 96    Histo1DPtr _h_ctheta, _h_P0,_h_P2;
 97    CounterPtr _wgtSum;
 98    /// @}
 99
100
101  };
102
103
104  RIVET_DECLARE_PLUGIN(BABAR_2008_I781294);
105
106}