1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief Xi_c0 -> Xi-pi+ asymmetry
  9  class CLEO_2000_I537236 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(CLEO_2000_I537236);
 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 histograms
 24      book(_h_ctheta, "ctheta"  , 20,-1,1);
 25    }
 26
 27
 28    /// Perform the per-event analysis
 29    void analyze(const Event& event) {
 30      // loop over Xi_c0 baryons
 31      for( const Particle& Xic : apply<UnstableParticles>(event, "UFS").particles(Cuts::abspid==4132)) {
 32	int sign = Xic.pid()/4132;
 33	if(Xic.children().size()!=2) continue;
 34	Particle baryon1,meson1;
 35	if(Xic.children()[0].pid()==sign*3312 &&
 36	   Xic.children()[1].pid()==sign*211) {
 37	  baryon1 = Xic.children()[0];
 38	  meson1  = Xic.children()[1];
 39	}
 40	else if(Xic.children()[1].pid()==sign*3312 &&
 41		Xic.children()[0].pid()==sign*211) {
 42	  baryon1 = Xic.children()[1];
 43	  meson1  = Xic.children()[0];
 44	}
 45	else
 46	  continue;
 47	Particle baryon2,meson2;
 48	if(baryon1.children()[0].pid()== sign*3122 &&
 49	   baryon1.children()[1].pid()==-sign*211) {
 50	  baryon2 = baryon1.children()[0];
 51	  meson2  = baryon1.children()[1];
 52	}
 53	else if(baryon1.children()[1].pid()== sign*3122 &&
 54		baryon1.children()[0].pid()==-sign*211) {
 55	  baryon2 = baryon1.children()[1];
 56	  meson2  = baryon1.children()[0];
 57	}
 58	else
 59	  continue;
 60	// first boost to the Xic rest frame
 61	LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(Xic.momentum().betaVec());
 62	FourMomentum pbaryon1 = boost1.transform(baryon1.momentum());
 63	FourMomentum pbaryon2 = boost1.transform(baryon2.momentum());
 64	// to lambda rest frame
 65	LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(pbaryon1.betaVec());
 66	Vector3 axis = pbaryon1.p3().unit();
 67	FourMomentum pp = boost2.transform(pbaryon2);
 68	// calculate angle
 69	double cTheta = pp.p3().unit().dot(axis);
 70	_h_ctheta->fill(cTheta,1.);
 71      }
 72    }
 73
 74    pair<double,double> calcAlpha(Histo1DPtr hist) {
 75      if(hist->numEntries()==0.) return make_pair(0.,0.);
 76      double sum1(0.),sum2(0.);
 77      for (const auto& bin : hist->bins() ) {
 78        double Oi = bin.sumW();
 79        if(Oi==0.) continue;
 80        double ai = 0.5*(bin.xMax()-bin.xMin());
 81        double bi = 0.5*ai*(bin.xMax()+bin.xMin());
 82        double Ei = bin.errW();
 83        sum1 += sqr(bi/Ei);
 84        sum2 += bi/sqr(Ei)*(Oi-ai);
 85      }
 86      return make_pair(sum2/sum1,sqrt(1./sum1));
 87    }
 88
 89    /// Normalise histograms etc., after the run
 90    void finalize() {
 91      normalize(_h_ctheta);
 92      Estimate0DPtr _h_alpha;
 93      book(_h_alpha,1,1,1);
 94      pair<double,double> alpha = calcAlpha(_h_ctheta);
 95      _h_alpha->set(alpha.first, alpha.second);
 96    }
 97
 98    /// @}
 99
100
101    /// @name Histograms
102    /// @{
103    Histo1DPtr _h_ctheta;
104    /// @}
105
106
107  };
108
109
110  RIVET_DECLARE_PLUGIN(CLEO_2000_I537236);
111
112
113}