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#include "Rivet/Projections/ParisiTensor.hh"
  9//#include "Rivet/Projections/FastJets.hh"
 10
 11namespace Rivet {
 12
 13
 14  /// @brief Event shapes at MZ
 15  class SLD_1995_I378545 : public Analysis {
 16  public:
 17
 18    /// Constructor
 19    RIVET_DEFAULT_ANALYSIS_CTOR(SLD_1995_I378545);
 20
 21
 22    /// @name Analysis methods
 23    ///@{
 24
 25    /// Book histograms and initialise projections before the run
 26    void init() {
 27      const FinalState fs;
 28      declare(fs, "FS");
 29      const ChargedFinalState cfs;
 30      declare(cfs, "CFS");
 31      Sphericity sphere(fs);
 32      declare(sphere, "Sphericity");
 33      declare(Thrust    (fs), "Thrust"    );
 34      declare(Hemispheres(sphere), "Hemispheres");
 35      declare(ParisiTensor(fs), "Parisi");
 36      // declare(FastJets(fs, JetAlg::DURHAM, 0.7), "DurhamJets");
 37      // declare(FastJets(fs, JetAlg::JADE  , 0.7), "JadeJets"  );
 38      // histograms
 39      book(_histThrust     , 2, 1, 1);
 40      book(_histMJetHigh   , 3, 1, 1);
 41      book(_histBT         , 4, 1, 1);
 42      book(_histBW         , 5, 1, 1);
 43      book(_histOblateness , 6, 1, 1);
 44      book(_histC          , 7, 1, 1);
 45      // book(_histy23JADE    , 8, 1, 1);
 46      // book(_histy23Durham  ,12, 1, 1);
 47      book(_histEEC        ,14, 1, 1);
 48      book(_histAEEC       ,15, 1, 1);
 49      book(_histJCEF       ,16, 1, 1);
 50    }
 51
 52
 53    /// Perform the per-event analysis
 54    void analyze(const Event& event) {
 55      const FinalState& fs = apply<FinalState>(event, "FS");
 56      if ( fs.particles().size() < 2) vetoEvent;
 57      // thrust
 58      const Thrust& thrust = apply<Thrust>(event, "Thrust");
 59      _histThrust    ->fill(1.-thrust.thrust());
 60      _histOblateness->fill(thrust.oblateness() );
 61      // visible energy
 62      double Evis = 0.0;
 63      for (const Particle& p : fs.particles()) {
 64        Evis += p.E();
 65      }
 66      double Evis2 = sqr(Evis);
 67      // (A)EEC
 68      // Need iterators since second loop starts at current outer loop iterator, i.e. no "foreach" here!
 69      for (Particles::const_iterator p_i = fs.particles().begin(); p_i != fs.particles().end(); ++p_i) {
 70        for (Particles::const_iterator p_j = p_i; p_j != fs.particles().end(); ++p_j) {
 71          const Vector3 mom3_i = p_i->momentum().p3();
 72          const Vector3 mom3_j = p_j->momentum().p3();
 73          const double energy_i = p_i->momentum().E();
 74          const double energy_j = p_j->momentum().E();
 75          const double thetaij = 180.*mom3_i.unit().angle(mom3_j.unit())/M_PI;
 76          double eec = (energy_i*energy_j) / Evis2;
 77          if (p_i != p_j)  eec *= 2.;
 78          _histEEC->fill(           thetaij, eec);
 79          if (thetaij <90.){
 80            _histAEEC->fill(thetaij, -eec);
 81          }
 82          else {
 83            _histAEEC->fill(180.-thetaij, eec);
 84          }
 85        }
 86      }
 87      // hemisphere related
 88      const Hemispheres& hemi = apply<Hemispheres>(event, "Hemispheres");
 89      _histMJetHigh->fill(hemi.scaledM2high());
 90      _histBW->fill(hemi.Bmax() );
 91      _histBT->fill(hemi.Bsum() );
 92      // C-parameter
 93      const ParisiTensor& parisi = apply<ParisiTensor>(event, "Parisi");
 94      _histC->fill(parisi.C());
 95      // jets
 96      // const FastJets&  durjet = apply<FastJets>(event, "DurhamJets");
 97      // const FastJets& jadejet = apply<FastJets>(event, "JadeJets");
 98      // if (durjet .clusterSeq()) _histy23Durham->fill( durjet.clusterSeq()->exclusive_ymerge_max(2));
 99      // if (jadejet.clusterSeq()) _histy23JADE  ->fill(jadejet.clusterSeq()->exclusive_ymerge_max(2));
100      // jet cone
101      Vector3 jetAxis=thrust.thrustAxis();
102      if(hemi.highMassDirection()) jetAxis *=-1.;
103      for (const Particle & p : fs.particles()) {
104        const double thetaij = 180.*jetAxis.angle(p.p3().unit())/M_PI;
105        double jcef = p.E()/ Evis;
106        _histJCEF->fill(thetaij,jcef);
107      }
108    }
109
110
111    /// Normalise histograms etc., after the run
112    void finalize() {
113      scale(_histThrust     , 1./sumOfWeights());
114      scale(_histMJetHigh   , 1./sumOfWeights());
115      scale(_histBT         , 1./sumOfWeights());
116      scale(_histBW         , 1./sumOfWeights());
117      scale(_histOblateness , 1./sumOfWeights());
118      scale(_histC          , 1./sumOfWeights());
119      // scale(_histy23JADE    , 1./sumOfWeights());
120      // scale(_histy23Durham  , 1./sumOfWeights());
121      scale(_histEEC        , 180./M_PI/sumOfWeights());
122      scale(_histAEEC       , 180./M_PI/sumOfWeights());
123      scale(_histJCEF       , 180./M_PI/sumOfWeights());
124    }
125
126    ///@}
127
128
129    /// @name Histograms
130    ///@{
131    Histo1DPtr  _histThrust, _histMJetHigh, _histBT, _histBW;
132    Histo1DPtr _histOblateness, _histC, _histy23JADE, _histy23Durham;
133    Histo1DPtr _histEEC, _histAEEC, _histJCEF;
134
135    ///@}
136
137
138  };
139
140
141  RIVET_DECLARE_PLUGIN(SLD_1995_I378545);
142
143}