1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/ChargedFinalState.hh"
  5#include "Rivet/Projections/Sphericity.hh"
  6#include "Rivet/Projections/Thrust.hh"
  7#include "Rivet/Projections/Hemispheres.hh"
  8
  9namespace Rivet {
 10
 11
 12  /// @brief Event shapes at 29 GeV
 13  class MARKII_1988_I246184 : public Analysis {
 14  public:
 15
 16    /// Constructor
 17    RIVET_DEFAULT_ANALYSIS_CTOR(MARKII_1988_I246184);
 18
 19
 20    /// @name Analysis methods
 21    ///@{
 22
 23    /// Book histograms and initialise projections before the run
 24    void init() {
 25      const FinalState fs;
 26      declare(fs, "FS");
 27      const ChargedFinalState cfs;
 28      declare(cfs, "CFS");
 29      Sphericity sphere(fs);
 30      declare(sphere, "Sphericity");
 31      declare(Thrust    (fs), "Thrust"    );
 32      declare(Hemispheres(sphere), "Hemispheres");
 33      // histograms
 34      unsigned int ioff=18;
 35      for(unsigned int ix=0;ix<3;++ix) {
 36	book(_histAplanarity [ix]  , 1+ioff*ix, 1, 1);
 37	book(_histQx         [ix]  , 2+ioff*ix, 1, 1);
 38	book(_histQ2Q1       [ix]  , 3+ioff*ix, 1, 1);
 39	book(_histSphericity [ix]  , 4+ioff*ix, 1, 1);
 40	book(_histThrust     [ix]  , 5+ioff*ix, 1, 1);
 41	book(_histMinor      [ix]  , 6+ioff*ix, 1, 1);
 42	book(_histOblateness [ix]  , 7+ioff*ix, 1, 1);
 43	book(_histMJetBroad  [ix]  , 8+ioff*ix, 1, 1);
 44	book(_histMJetSlim   [ix]  , 9+ioff*ix, 1, 1);
 45	book(_histMJetDiff   [ix]  ,10+ioff*ix, 1, 1);
 46	book(_histScaledMom  [ix]  ,15+ioff*ix, 1, 1);
 47	book(_histPt2S       [ix]  ,11+ioff*ix, 1, 1);
 48	book(_histPtS        [ix]  ,12+ioff*ix, 1, 1);
 49	book(_histPtSIn      [ix]  ,14+ioff*ix, 1, 1);
 50	book(_histPtSOut     [ix]  ,13+ioff*ix, 1, 1);
 51	book(_histRapidityS  [ix]  ,16+ioff*ix, 1, 1);
 52	book(_histTheta      [ix]  ,17+ioff*ix, 1, 1);
 53	book(_histETheta     [ix]  ,18+ioff*ix, 1, 1);
 54      }
 55    }
 56
 57
 58    /// Perform the per-event analysis
 59    void analyze(const Event& event) {
 60      // Sphericity related
 61      const Sphericity& sphericity = apply<Sphericity>(event, "Sphericity");
 62      for(unsigned int ix=0;ix<3;++ix) {
 63	_histSphericity[ix]->fill(sphericity.sphericity());
 64	_histAplanarity[ix]->fill(sphericity.aplanarity());
 65	_histQx        [ix]->fill((sphericity.lambda1()-sphericity.lambda2())/sqrt(3.));
 66	_histQ2Q1      [ix]->fill(sphericity.lambda2()-sphericity.lambda3());
 67      }
 68      // thrust related
 69      const Thrust& thrust = apply<Thrust>(event, "Thrust");
 70      for(unsigned int ix=0;ix<3;++ix) {
 71	_histThrust    [ix]->fill(thrust.thrust());
 72	_histMinor     [ix]->fill(thrust.thrustMinor());
 73	_histOblateness[ix]->fill(thrust.oblateness());
 74      }
 75      // hemisphere related
 76      const Hemispheres& hemi = apply<Hemispheres>(event, "Hemispheres");
 77      double mWide = hemi.scaledM2high(), mNarrow = hemi.scaledM2low();
 78      if(!hemi.massMatchesBroadening()) swap(mWide,mNarrow);
 79      for(unsigned int ix=0;ix<3;++ix) {
 80	_histMJetBroad[ix]->fill(mWide);
 81	_histMJetSlim [ix]->fill(mNarrow);
 82	_histMJetDiff [ix]->fill(hemi.scaledM2diff());
 83      }
 84      // dists w.r.t sphericity axis
 85      const FinalState& cfs = apply<FinalState>(event, "CFS");
 86      for (const Particle& p : cfs.particles()) {
 87        // Get momentum and energy of each particle.
 88        const Vector3 mom3 = p.p3();
 89        const double energy = p.E();
 90        // Scaled momenta.
 91        const double mom = mom3.mod();
 92        const double scaledMom = 2.*mom/sqrtS();
 93	// Get momenta components w.r.t. thrust and sphericity.
 94        const double pTinS = dot(mom3, sphericity.sphericityMajorAxis());
 95        const double pToutS = dot(mom3, sphericity.sphericityMinorAxis());
 96	double pT2 = sqr(pTinS)+sqr(pToutS);
 97	double pT  = sqrt(pT2);
 98        const double momS = dot(sphericity.sphericityAxis(), mom3);
 99        const double rapidityS = 0.5 * std::log((energy + momS) / (energy - momS));
100	// angle
101	double theta = sphericity.sphericityAxis().angle(mom3)/M_PI*180.;
102	if(theta>90.) theta=180.-theta;
103	// fill histos
104	for(unsigned int ix=0;ix<3;++ix) {
105	  _histScaledMom[ix]->fill(scaledMom);
106	  _histPt2S     [ix]->fill(fabs(pT2/GeV));
107	  _histPtS      [ix]->fill(fabs(pT/GeV));
108	  _histPtSIn    [ix]->fill(fabs(pTinS/GeV));
109	  _histPtSOut   [ix]->fill(fabs(pToutS/GeV));
110	  _histRapidityS[ix]->fill(fabs(rapidityS));
111	  _histTheta    [ix]->fill(theta);
112	}
113      }
114      // energy flow includes neutral w.r.t sphericity axis
115      const FinalState& fs = apply<FinalState>(event, "FS");
116      for (const Particle& p : fs.particles()) {
117        // Get momentum and energy of each particle.
118        const Vector3 mom3 = p.p3();
119        const double energy = p.E();
120	// angle
121	double theta = sphericity.sphericityAxis().angle(mom3)/M_PI*180.;
122	if(theta>90.) theta=180.-theta;
123	// fill histos
124	for(unsigned int ix=0;ix<3;++ix) {
125	  _histETheta   [ix]->fill(theta,energy);
126	}
127      }
128    }
129
130    /// Normalise histograms etc., after the run
131    void finalize() {
132      for(unsigned int ix=0;ix<3;++ix) {
133	scale(_histAplanarity [ix]  ,1./sumOfWeights());
134	scale(_histQx         [ix]  ,1./sumOfWeights());
135	scale(_histQ2Q1       [ix]  ,1./sumOfWeights());
136	scale(_histSphericity [ix]  ,1./sumOfWeights());
137	scale(_histThrust     [ix]  ,1./sumOfWeights());
138	scale(_histMinor      [ix]  ,1./sumOfWeights());
139	scale(_histOblateness [ix]  ,1./sumOfWeights());
140	scale(_histMJetBroad  [ix]  ,1./sumOfWeights());
141	scale(_histMJetSlim   [ix]  ,1./sumOfWeights());
142	scale(_histMJetDiff   [ix]  ,1./sumOfWeights());
143	scale(_histScaledMom  [ix]  ,1./sumOfWeights());
144	scale(_histPt2S       [ix]  ,1./sumOfWeights());
145	scale(_histPtS        [ix]  ,1./sumOfWeights());
146	scale(_histPtSIn      [ix]  ,1./sumOfWeights());
147	scale(_histPtSOut     [ix]  ,1./sumOfWeights());
148	scale(_histRapidityS  [ix]  ,1./sumOfWeights());
149	scale(_histTheta      [ix]  ,1./sumOfWeights());
150	scale(_histETheta     [ix]  ,1./sumOfWeights());
151      }
152    }
153
154    ///@}
155
156
157    /// @name Histograms
158    ///@{
159    Histo1DPtr _histAplanarity[3],_histQx[3],_histQ2Q1[3],_histSphericity[3];
160    Histo1DPtr _histThrust[3],_histMinor[3],_histOblateness[3];
161    Histo1DPtr _histMJetBroad[3],_histMJetSlim[3],_histMJetDiff[3];
162    Histo1DPtr _histScaledMom[3],_histPt2S[3],_histPtS[3],_histPtSIn[3],_histPtSOut[3],_histRapidityS[3];
163    Histo1DPtr _histTheta[3],_histETheta[3];
164    ///@}
165
166
167  };
168
169
170  RIVET_DECLARE_PLUGIN(MARKII_1988_I246184);
171
172}