1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/ChargedFinalState.hh"
  4#include "Rivet/Projections/Thrust.hh"
  5#include "Rivet/Projections/Hemispheres.hh"
  6#include "Rivet/Projections/ParisiTensor.hh"
  7#include "Rivet/Projections/FastJets.hh"
  8
  9#define I_KNOW_THE_INITIAL_QUARKS_PROJECTION_IS_DODGY_BUT_NEED_TO_USE_IT
 10#include "Rivet/Projections/InitialQuarks.hh"
 11
 12namespace Rivet {
 13
 14
 15  /// @brief Event shapes at 197 GeV
 16  class L3_2008_I825820 : public Analysis {
 17  public:
 18
 19    /// Constructor
 20    RIVET_DEFAULT_ANALYSIS_CTOR(L3_2008_I825820);
 21
 22
 23    /// @name Analysis methods
 24    /// @{
 25
 26    /// Book histograms and initialise projections before the run
 27    void init() {
 28      // Projections to use
 29      const FinalState FS;
 30      declare(FS, "FS");
 31      const ChargedFinalState CFS;
 32      declare(CFS, "CFS");
 33      const Thrust thrust(FS);
 34      declare(thrust, "Thrust");
 35      declare(ParisiTensor(FS), "Parisi");
 36      declare(Hemispheres(thrust), "Hemispheres");
 37      declare(InitialQuarks(), "InitialQuarks");
 38      declare(FastJets(FS, JetAlg::JADE, 0.7), "Jets");
 39
 40      // histograms
 41      book(_h_T         ,1,1,1);
 42      book(_h_T_udsc    ,1,1,2);
 43      book(_h_T_bottom  ,1,1,3);
 44      book(_h_rho       ,2,1,1);
 45      book(_h_rho_udsc  ,2,1,2);
 46      book(_h_rho_bottom,2,1,3);
 47      book(_h_B_T       ,3,1,1);
 48      book(_h_B_T_udsc  ,3,1,2);
 49      book(_h_B_T_bottom,3,1,3);
 50      book(_h_B_W       ,4,1,1);
 51      book(_h_B_W_udsc  ,4,1,2);
 52      book(_h_B_W_bottom,4,1,3);
 53      book(_h_C         ,5,1,1);
 54      book(_h_C_udsc    ,5,1,2);
 55      book(_h_C_bottom  ,5,1,3);
 56      book(_h_y23       ,6,1,1);
 57      book(_h_y23_udsc  ,6,1,2);
 58      book(_h_y23_bottom,6,1,3);
 59      book(_sumW_udsc, "_sumW_udsc");
 60      book(_sumW_b, "_sumW_b");
 61    }
 62
 63
 64    /// Perform the per-event analysis
 65    void analyze(const Event& event) {
 66
 67      int flavour = 0;
 68      const InitialQuarks& iqf = apply<InitialQuarks>(event, "InitialQuarks");
 69      Particles quarks;
 70      if ( iqf.particles().size() == 2 ) {
 71        flavour = iqf.particles().front().abspid();
 72        quarks  = iqf.particles();
 73      } else {
 74        map<int, Particle> quarkmap;
 75        for (const Particle& p : iqf.particles()) {
 76          if (quarkmap.find(p.pid()) == quarkmap.end()) quarkmap[p.pid()] = p;
 77          else if (quarkmap[p.pid()].E() < p.E()) quarkmap[p.pid()] = p;
 78        }
 79        double max_energy = 0.;
 80        for (int i = 1; i <= 5; ++i) {
 81          double energy = 0.;
 82          if (quarkmap.find(i) != quarkmap.end())
 83            energy += quarkmap[ i].E();
 84          if (quarkmap.find(-i) != quarkmap.end())
 85            energy += quarkmap[-i].E();
 86          if (energy > max_energy)
 87            flavour = i;
 88        }
 89        if (quarkmap.find(flavour) != quarkmap.end())
 90          quarks.push_back(quarkmap[flavour]);
 91        if (quarkmap.find(-flavour) != quarkmap.end())
 92          quarks.push_back(quarkmap[-flavour]);
 93      }
 94
 95      // Flavour label
 96      int iflav = (flavour == PID::DQUARK || flavour == PID::UQUARK ||
 97                   flavour == PID::SQUARK || flavour == PID::CQUARK) ?
 98        1 : (flavour == PID::BQUARK) ? 5 : 0;
 99
100      // Update weight sums
101      if (iflav == 1) {
102        _sumW_udsc->fill();
103      } else if (iflav == 5) {
104        _sumW_b->fill();
105      }
106
107      // Thrust
108      const Thrust& thrust = apply<Thrust>(event, "Thrust");
109      if (iflav == 1) {
110        _h_T_udsc->fill(thrust.thrust());
111      } else if (iflav == 5) {
112        _h_T_bottom->fill(thrust.thrust());
113      }
114      _h_T->fill(thrust.thrust());
115
116
117      // The hemisphere variables
118      const Hemispheres& hemisphere = apply<Hemispheres>(event, "Hemispheres");
119      if (iflav == 1) {
120        _h_rho_udsc->fill(hemisphere.scaledM2high());
121        _h_B_T_udsc->fill(hemisphere.Bsum());
122        _h_B_W_udsc->fill(hemisphere.Bmax());
123      } else if (iflav == 5) {
124        _h_rho_bottom->fill(hemisphere.scaledM2high());
125        _h_B_T_bottom->fill(hemisphere.Bsum());
126        _h_B_W_bottom->fill(hemisphere.Bmax());
127      }
128      _h_rho->fill(hemisphere.scaledM2high());
129      _h_B_T->fill(hemisphere.Bsum());
130      _h_B_W->fill(hemisphere.Bmax());
131
132      const ParisiTensor& parisi = apply<ParisiTensor>(event, "Parisi");
133      if (iflav == 1) {
134        _h_C_udsc->fill(parisi.C());
135      } else if (iflav == 5) {
136        _h_C_bottom->fill(parisi.C());
137      }
138      _h_C->fill(parisi.C());
139      // y_23
140      const FastJets& durjet = apply<FastJets>(event, "Jets");
141      const double y23 = durjet.clusterSeq()->exclusive_ymerge_max(2);
142      if (iflav == 1) {
143        _h_y23_udsc->fill(y23);
144      } else if (iflav == 5) {
145        _h_y23_bottom->fill(y23);
146      }
147      _h_y23->fill(y23);
148    }
149
150    void finalize() {
151      scale(_h_T_udsc  , 1./ *_sumW_udsc);
152      scale(_h_T_bottom, 1./ *_sumW_b   );
153      normalize(_h_T);
154      scale(_h_rho_udsc  , 1./ *_sumW_udsc);
155      scale(_h_rho_bottom, 1./ *_sumW_b   );
156      normalize(_h_rho);
157      scale(_h_B_T_udsc  , 1./ *_sumW_udsc);
158      scale(_h_B_T_bottom, 1./ *_sumW_b   );
159      normalize(_h_B_T);
160      scale(_h_B_W_udsc  , 1./ *_sumW_udsc);
161      scale(_h_B_W_bottom, 1./ *_sumW_b   );
162      normalize(_h_B_W);
163      scale(_h_C_udsc  , 1./ *_sumW_udsc);
164      scale(_h_C_bottom, 1./ *_sumW_b   );
165      normalize(_h_C);
166      scale(_h_y23_udsc  , 1./ *_sumW_udsc);
167      scale(_h_y23_bottom, 1./ *_sumW_b   );
168      normalize(_h_y23);
169    }
170    /// @}
171
172
173    /// @name Histograms
174    /// @{
175    Histo1DPtr _h_T  ,_h_T_udsc  , _h_T_bottom;
176    Histo1DPtr _h_rho,_h_rho_udsc, _h_rho_bottom;
177    Histo1DPtr _h_B_T,_h_B_T_udsc, _h_B_T_bottom;
178    Histo1DPtr _h_B_W,_h_B_W_udsc, _h_B_W_bottom;
179    Histo1DPtr _h_C  ,_h_C_udsc  , _h_C_bottom;
180    Histo1DPtr _h_y23,_h_y23_udsc, _h_y23_bottom;
181    CounterPtr _sumW_udsc, _sumW_b;
182    /// @}
183
184
185  };
186
187
188  RIVET_DECLARE_PLUGIN(L3_2008_I825820);
189
190}