1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FastJets.hh"
  4#include "Rivet/Projections/FinalState.hh"
  5#include "Rivet/Projections/ChargedFinalState.hh"
  6#include "Rivet/Projections/Sphericity.hh"
  7#include "Rivet/Projections/Thrust.hh"
  8#include "Rivet/Projections/Hemispheres.hh"
  9#include "Rivet/Projections/ParisiTensor.hh"
 10
 11namespace Rivet {
 12
 13
 14  /// @brief Event shapes at MZ
 15  class L3_1992_I334954 : public Analysis {
 16  public:
 17
 18    /// Constructor
 19    RIVET_DEFAULT_ANALYSIS_CTOR(L3_1992_I334954);
 20
 21
 22    /// @name Analysis methods
 23    ///@{
 24
 25    /// Book histograms and initialise projections before the run
 26    void init() {
 27      // projections
 28      const FinalState fs;
 29      declare(fs, "FS");
 30      const ChargedFinalState cfs;
 31      declare(cfs, "CFS");
 32      declare(Sphericity(fs), "Sphericity");
 33      Thrust thrust(fs);
 34      declare(thrust, "Thrust"    );
 35      declare(Hemispheres(thrust), "Hemispheres");
 36      declare(ParisiTensor(fs), "Parisi");
 37      declare(FastJets(cfs, JetAlg::DURHAM, 0.7), "DurhamJets");
 38      declare(FastJets(cfs, JetAlg::JADE       ), "JadeJets"  );
 39      // histograms
 40      book(_histThrust    ,  1, 1, 1);
 41      book(_histMajor     ,  2, 1, 1);
 42      book(_histMinor     ,  3, 1, 1);
 43      book(_histOblateness,  4, 1, 1);
 44      book(_histJade      ,  6, 1, 1);
 45      book(_histDurham    ,  7, 1, 1);
 46      book(_histH3        ,  8, 1, 1);
 47      book(_histH4        ,  9, 1, 1);
 48      book(_histSphericity, 10, 1, 1);
 49      book(_histAplanarity, 11, 1, 1);
 50      book(_histC         , 12, 1, 1);
 51      book(_histD         , 13, 1, 1);
 52      book(_histMJetHeavy , 14, 1, 1);
 53      book(_histMJetLight , 15, 1, 1);
 54      book(_histMult      , 16, 1, 1);
 55    }
 56
 57
 58    /// Perform the per-event analysis
 59    void analyze(const Event& event) {
 60      // 5 charged particles
 61      const FinalState& cfs = apply<FinalState>(event, "CFS");
 62      if(cfs.particles().size()<5) vetoEvent;
 63      // charged particle mult
 64      const string multi_edge = std::to_string(cfs.particles().size()) + ".0";
 65      _histMult->fill(multi_edge);
 66      // Sphericity related
 67      const Sphericity& sphericity = apply<Sphericity>(event, "Sphericity");
 68      _histSphericity->fill(sphericity.sphericity());
 69      _histAplanarity->fill(sphericity.aplanarity());
 70      // thrust related
 71      const Thrust& thrust = apply<Thrust>(event, "Thrust");
 72      _histThrust    ->fill(thrust.thrust());
 73      _histMajor     ->fill(thrust.thrustMajor());
 74      _histMinor     ->fill(thrust.thrustMinor());
 75      _histOblateness->fill(thrust.oblateness());
 76      // hemisphere related
 77      const Hemispheres& hemi = apply<Hemispheres>(event, "Hemispheres");
 78      _histMJetHeavy->fill(hemi.scaledM2high());
 79      _histMJetLight->fill(hemi.scaledM2low());
 80      // C and D
 81      const ParisiTensor& parisi = apply<ParisiTensor>(event, "Parisi");
 82      _histC->fill(parisi.C());
 83      _histD->fill(parisi.D());
 84      // jet rates
 85      const FastJets& durjet = apply<FastJets>(event, "DurhamJets");
 86      double y23 = durjet.clusterSeq()->exclusive_ymerge_max(2);
 87      _histDurham->fill(y23);
 88      const FastJets& jadejet = apply<FastJets>(event, "JadeJets");
 89      y23 = jadejet.clusterSeq()->exclusive_ymerge_max(2);
 90      _histJade->fill(y23);
 91      // fox-wolfram moments
 92      double H3=0, H4=0;
 93      const FinalState&  fs = apply<FinalState>(event, "FS");
 94      for(const Particle & p1 : fs.particles()) {
 95	double modp1 = p1.p3().mod();
 96	for(const Particle & p2 : fs.particles()) {
 97	  double modp2 = p2.p3().mod();
 98	  double cTheta = p1.p3().dot(p2.p3())/modp1/modp2;
 99	  double pre = modp1*modp2/sqr(sqrtS());
100	  H3 +=   0.5*pre*cTheta*(5.*sqr(cTheta)-3.);
101	  H4 += 0.125*pre*((35.*sqr(cTheta)-30.)*sqr(cTheta)+3.);
102	}
103      }
104      _histH3->fill(H3);
105      _histH4->fill(H4);
106    }
107
108
109    /// Normalise histograms etc., after the run
110    void finalize() {
111      scale(_histThrust    , 1./sumOfWeights());
112      scale(_histMajor     , 1./sumOfWeights());
113      scale(_histMinor     , 1./sumOfWeights());
114      scale(_histOblateness, 1./sumOfWeights());
115      scale(_histJade      , 1./sumOfWeights());
116      scale(_histDurham    , 1./sumOfWeights());
117      scale(_histH3        , 1./sumOfWeights());
118      scale(_histH4        , 1./sumOfWeights());
119      scale(_histSphericity, 1./sumOfWeights());
120      scale(_histAplanarity, 1./sumOfWeights());
121      scale(_histC         , 1./sumOfWeights());
122      scale(_histD         , 1./sumOfWeights());
123      scale(_histMJetHeavy , 1./sumOfWeights());
124      scale(_histMJetLight , 1./sumOfWeights());
125      // percentage and bin width
126      scale(_histMult, 100./sumOfWeights());
127    }
128
129    ///@}
130
131
132    /// @name Histograms
133    ///@{
134    Histo1DPtr _histThrust, _histMajor, _histMinor, _histOblateness;
135    Histo1DPtr _histJade,_histDurham;
136    Histo1DPtr _histH3,_histH4;
137    Histo1DPtr _histSphericity, _histAplanarity;
138    Histo1DPtr _histC, _histD;
139    Histo1DPtr _histMJetHeavy, _histMJetLight;
140    BinnedHistoPtr<string> _histMult;
141    ///@}
142
143
144  };
145
146
147  RIVET_DECLARE_PLUGIN(L3_1992_I334954);
148
149}