1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/Beam.hh"
  4#include "Rivet/Projections/FinalState.hh"
  5#include "Rivet/Projections/UnstableParticles.hh"
  6#include "Rivet/Projections/Sphericity.hh"
  7
  8namespace Rivet {
  9
 10
 11  /// @brief D*+/- polarization at 29 GeV
 12  class HRS_1987_I250823 : public Analysis {
 13  public:
 14
 15    /// Constructor
 16    RIVET_DEFAULT_ANALYSIS_CTOR(HRS_1987_I250823);
 17
 18
 19    /// @name Analysis methods
 20    /// @{
 21
 22    /// Book histograms and initialise projections before the run
 23    void init() {
 24      declare(Beam(), "Beams");
 25      declare(Sphericity(FinalState()), "Sphericity");
 26      declare(UnstableParticles(), "UFS" );
 27      for(unsigned int i=0;i<9;++i) {
 28       	unsigned int ix(0),iy(1);
 29       	if(i<3) {
 30       	  ix=1;
 31       	  iy=i+1;
 32       	}
 33       	else if(i==3) {
 34       	  ix=2;
 35       	}
 36       	else if(i<6) {
 37       	  ix=3;
 38       	  iy=i-3;
 39       	}
 40       	else {
 41       	  ix=i-2;
 42       	}
 43	book(_p_rho00[i],ix,iy,1);
 44	book(_p_rho11[i],ix,iy,2);
 45	book(_p_rho10[i],ix,iy,3);
 46      }
 47    }
 48
 49
 50    /// Perform the per-event analysis
 51    void analyze(const Event& event) {
 52      // Get beams and average beam momentum
 53      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 54      const double meanBeamMom = ( beams.first.p3().mod() +
 55                                   beams.second.p3().mod() ) / 2.0;
 56      MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
 57      Vector3 axis;
 58      if(beams.first.pid()>0)
 59	axis = beams.first .momentum().p3().unit();
 60      else
 61	axis = beams.second.momentum().p3().unit();
 62
 63      // sphericity, to define an axis
 64      const Sphericity& sphericity = apply<Sphericity>(event, "Sphericity");
 65
 66      // loop over the particles
 67      for(const Particle& p : apply<UnstableParticles>(event, "UFS").particles(Cuts::abspid==413)) {
 68	int sign = p.pid()/p.abspid();
 69	Particle decay;
 70	double xE = p.momentum().E()/meanBeamMom;
 71	if(p.children()[0].pid()==sign*421 &&
 72	   p.children()[1].pid()==sign*211) {
 73	  decay = p.children()[1];
 74	}
 75	else if(p.children()[1].pid()==sign*421 &&
 76		p.children()[0].pid()==sign*211) {
 77	  decay = p.children()[0];
 78	}
 79	else
 80	  continue;
 81	// axis and ctheta
 82	LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
 83	Vector3 e1z = p.p3().unit();
 84	FourMomentum pp = boost.transform(decay.momentum());
 85	Vector3 axis1 = boost.transform(decay.momentum()).p3().unit();
 86	double ctheta = e1z.dot(axis1);
 87	// y and z axis
 88	Vector3 e1y = e1z.cross(axis).unit();
 89	Vector3 e1x = e1y.cross(e1z).unit();
 90	double phi = atan2(e1y.dot(axis1),e1x.dot(axis1));
 91	double w1 = 0.5*(5.*sqr(ctheta)-1.);
 92	double w2 = -1.25*(1.-sqr(ctheta))*cos(2.*phi);
 93	double w3 = -1.25*sqrt(2)*ctheta*sqrt(1.-sqr(ctheta))*cos(phi);
 94	// fill the hists by x_E
 95	for(unsigned int ix=0;ix<6;++ix) {
 96	  _p_rho00[ix]->fill(xE,w1);
 97	  _p_rho11[ix]->fill(xE,w2);
 98	  _p_rho10[ix]->fill(xE,w3);
 99	}
100	_p_rho00[6]->fill(1.,w1);
101	_p_rho11[6]->fill(1.,w2);
102	_p_rho10[6]->fill(1.,w3);
103	// using jet axis
104	double pT = sqrt(sqr(sphericity.sphericityMajorAxis().dot(p.momentum().p3()))+
105			 sqr(sphericity.sphericityMinorAxis().dot(p.momentum().p3())));
106	Vector3 axis2;
107	if(p.momentum().p3().dot(sphericity.sphericityAxis())>=0.) {
108	  axis2 = sphericity.sphericityAxis();
109	}
110	else {
111	  axis2 =-sphericity.sphericityAxis();
112	}
113	Vector3 e2y = e1z.cross(axis2).unit();
114	Vector3 e2x = e2y.cross(e1z).unit();
115	// alpha and beta
116	phi = atan2(e2y.dot(axis1),e2x.dot(axis1));
117	w1 = 0.5*(5.*sqr(ctheta)-1.);
118	w2 = -1.25*(1.-sqr(ctheta))*cos(2.*phi);
119	w3 = -1.25*sqrt(2)*ctheta*sqrt(1.-sqr(ctheta))*cos(phi);
120	if(xE<0.4) continue;
121	_p_rho00[7]->fill(1.,w1);
122	_p_rho11[7]->fill(1.,w2);
123	_p_rho10[7]->fill(1.,w3);
124	if(pT<0.75) {
125	  _p_rho00[7]->fill(2.,w1);
126	  _p_rho11[7]->fill(2.,w2);
127	  _p_rho10[7]->fill(2.,w3);
128	}
129	else {
130	  _p_rho00[7]->fill(3.,w1);
131	  _p_rho11[7]->fill(3.,w2);
132	  _p_rho10[7]->fill(3.,w3);
133	}
134	if(sphericity.sphericity()>=0.1) {
135	  _p_rho00[8]->fill(1.,w1);
136	  _p_rho11[8]->fill(1.,w2);
137	  _p_rho10[8]->fill(1.,w3);
138	}
139	else {
140	  _p_rho00[8]->fill(2.,w1);
141	  _p_rho11[8]->fill(2.,w2);
142	  _p_rho10[8]->fill(2.,w3);
143	}
144      }
145    }
146
147
148    /// Normalise histograms etc., after the run
149    void finalize() {
150
151    }
152
153    /// @}
154
155
156    /// @name Histograms
157    /// @{
158    Profile1DPtr _p_rho00[9],_p_rho11[9],_p_rho10[9];
159    /// @}
160
161
162  };
163
164
165  RIVET_DECLARE_PLUGIN(HRS_1987_I250823);
166
167
168}