1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5#include "Rivet/Projections/ChargedFinalState.hh"
  6#include "Rivet/Projections/Thrust.hh"
  7#include "Rivet/Projections/Sphericity.hh"
  8#include "Rivet/Projections/ParisiTensor.hh"
  9#include "Rivet/Projections/Hemispheres.hh"
 10#include "Rivet/Projections/Beam.hh"
 11
 12namespace Rivet {
 13
 14
 15  /// ALEPH jet rates and event shapes at LEP 1 and 2
 16  class ALEPH_2004_I636645 : public Analysis {
 17  public:
 18
 19    RIVET_DEFAULT_ANALYSIS_CTOR(ALEPH_2004_I636645);
 20
 21
 22    void init() {
 23      _initialisedJets    = true;
 24      _initialisedSpectra = true;
 25      // TODO: According to the paper they seem to discard neutral particles
 26      //       between 1 and 2 GeV. That correction is included in the systematic
 27      //       uncertainties and overly complicated to program, so we ignore it.
 28      const FinalState fs;
 29      declare(fs, "FS");
 30      FastJets durhamjets(fs, JetAlg::DURHAM, 0.7, JetMuons::ALL, JetInvisibles::ALL);
 31      declare(durhamjets, "DurhamJets");
 32
 33      const Thrust thrust(fs);
 34      declare(thrust, "Thrust");
 35      declare(Sphericity(fs), "Sphericity");
 36      declare(ParisiTensor(fs), "Parisi");
 37      declare(Hemispheres(thrust), "Hemispheres");
 38
 39      const ChargedFinalState cfs;
 40      declare(Beam(), "Beams");
 41      declare(cfs, "CFS");
 42
 43      // Histos
 44      // offset for the event shapes and jets
 45      int offset = 0;
 46      switch (int(sqrtS()/GeV + 0.5)) {
 47        case 91: offset = 0; break;
 48        case 133: offset = 1; break;
 49        case 161: offset = 2; break;
 50        case 172: offset = 3; break;
 51        case 183: offset = 4; break;
 52        case 189: offset = 5; break;
 53        case 200: offset = 6; break;
 54        case 206: offset = 7; break;
 55        default:
 56          _initialisedJets = false;
 57      }
 58      // event shapes
 59      if(_initialisedJets) {
 60        book(_h_thrust ,offset+54, 1, 1);
 61        book(_h_heavyjetmass ,offset+62, 1, 1);
 62        book(_h_totaljetbroadening ,offset+70, 1, 1);
 63        book(_h_widejetbroadening ,offset+78, 1, 1);
 64        book(_h_cparameter ,offset+86, 1, 1);
 65        book(_h_thrustmajor ,offset+94, 1, 1);
 66        book(_h_thrustminor ,offset+102, 1, 1);
 67        book(_h_jetmassdifference ,offset+110, 1, 1);
 68        book(_h_aplanarity ,offset+118, 1, 1);
 69        if ( offset != 0 )
 70          book(_h_planarity, offset+125, 1, 1);
 71        book(_h_oblateness ,offset+133, 1, 1);
 72        book(_h_sphericity ,offset+141, 1, 1);
 73
 74        // Durham n->m jet resolutions
 75        book(_h_y_Durham[0] ,offset+149, 1, 1);   // y12 d149 ... d156
 76        book(_h_y_Durham[1] ,offset+157, 1, 1);   // y23 d157 ... d164
 77        if (offset < 6) { // there is no y34, y45 and y56 for 200 gev
 78          book(_h_y_Durham[2] ,offset+165, 1, 1); // y34 d165 ... d172, but not 171
 79          book(_h_y_Durham[3] ,offset+173, 1, 1); // y45 d173 ... d179
 80          book(_h_y_Durham[4] ,offset+180, 1, 1); // y56 d180 ... d186
 81        }
 82        else if (offset == 6) {
 83          _h_y_Durham[2] = Histo1DPtr();
 84          _h_y_Durham[3] = Histo1DPtr();
 85          _h_y_Durham[4] = Histo1DPtr();
 86        }
 87        else if (offset == 7) {
 88          book(_h_y_Durham[2] ,172, 1, 1);
 89          book(_h_y_Durham[3] ,179, 1, 1);
 90          book(_h_y_Durham[4] ,186, 1, 1);
 91        }
 92
 93        // Durham n-jet fractions
 94        book(_h_R_Durham[0] ,offset+187, 1, 1); // R1 d187 ... d194
 95        book(_h_R_Durham[1] ,offset+195, 1, 1); // R2 d195 ... d202
 96        book(_h_R_Durham[2] ,offset+203, 1, 1); // R3 d203 ... d210
 97        book(_h_R_Durham[3] ,offset+211, 1, 1); // R4 d211 ... d218
 98        book(_h_R_Durham[4] ,offset+219, 1, 1); // R5 d219 ... d226
 99        book(_h_R_Durham[5] ,offset+227, 1, 1); // R>=6 d227 ... d234
100      }
101      // offset for the charged particle distributions
102      offset = 0;
103      switch (int(sqrtS() + 0.5)) {
104        case 133: offset = 0; break;
105        case 161: offset = 1; break;
106        case 172: offset = 2; break;
107        case 183: offset = 3; break;
108        case 189: offset = 4; break;
109        case 197: offset = 5; break;
110        case 200: offset = 6; break;
111        case 206: offset = 7; break;
112        default:
113          _initialisedSpectra = false;
114      }
115      if (_initialisedSpectra) {
116        book(_h_xp , 2+offset, 1, 1);
117        book(_h_xi ,11+offset, 1, 1);
118        book(_h_xe ,19+offset, 1, 1);
119        book(_h_pTin  ,27+offset, 1, 1);
120        if (offset == 7)
121          book(_h_pTout, 35, 1, 1);
122        book(_h_rapidityT ,36+offset, 1, 1);
123        book(_h_rapidityS ,44+offset, 1, 1);
124      }
125      book(_weightedTotalChargedPartNum, "_weightedTotalChargedPartNum");
126      if (!_initialisedSpectra && !_initialisedJets) {
127        MSG_WARNING("CoM energy of events sqrt(s) = " << sqrtS()/GeV
128                    << " doesn't match any available analysis energy .");
129      }
130
131      book(mult, 1, 1, 1);
132    }
133
134
135    void analyze(const Event& e) {
136
137      const Thrust& thrust = apply<Thrust>(e, "Thrust");
138      const Sphericity& sphericity = apply<Sphericity>(e, "Sphericity");
139
140      if(_initialisedJets) {
141        bool LEP1 = isCompatibleWithSqrtS(91.2*GeV,0.01);
142        // event shapes
143        double thr = LEP1 ? thrust.thrust() : 1.0 - thrust.thrust();
144        _h_thrust->fill(thr);
145        _h_thrustmajor->fill(thrust.thrustMajor());
146        if(LEP1)
147          _h_thrustminor->fill(log(thrust.thrustMinor()));
148        else
149          _h_thrustminor->fill(thrust.thrustMinor());
150        _h_oblateness->fill(thrust.oblateness());
151
152        const Hemispheres& hemi = apply<Hemispheres>(e, "Hemispheres");
153        _h_heavyjetmass->fill(hemi.scaledM2high());
154        _h_jetmassdifference->fill(hemi.scaledM2diff());
155        _h_totaljetbroadening->fill(hemi.Bsum());
156        _h_widejetbroadening->fill(hemi.Bmax());
157
158        const ParisiTensor& parisi = apply<ParisiTensor>(e, "Parisi");
159        _h_cparameter->fill(parisi.C());
160
161        _h_aplanarity->fill(sphericity.aplanarity());
162        if(_h_planarity)
163          _h_planarity->fill(sphericity.planarity());
164        _h_sphericity->fill(sphericity.sphericity());
165
166        // Jet rates
167        const FastJets& durjet = apply<FastJets>(e, "DurhamJets");
168        double log10e = log10(exp(1.));
169        if (durjet.clusterSeq()) {
170          double logynm1=0.;
171          double logyn;
172          for (size_t i=0; i<5; ++i) {
173            const double yn = durjet.clusterSeq()->exclusive_ymerge_max(i+1);
174            if (yn<=0.0) continue;
175            logyn = -log(yn);
176            if (_h_y_Durham[i]) {
177              _h_y_Durham[i]->fill(logyn);
178            }
179            if(!LEP1) logyn *= log10e;
180            for (const auto& b : _h_R_Durham[i]->bins()) {
181              const double xMin = b.xMin();
182              if (-xMin <= logynm1)  break;
183              if (-xMin<logyn)  _h_R_Durham[i]->fill(b.xMid(), b.xWidth());
184            }
185            logynm1 = logyn;
186          }
187          for (const auto& b : _h_R_Durham[5]->bins()) {
188            if (-b.xMin() <= logynm1)  break;
189            _h_R_Durham[5]->fill(b.xMid(), b.xWidth());
190          }
191        }
192        if( !_initialisedSpectra) {
193          const ChargedFinalState& cfs = apply<ChargedFinalState>(e, "CFS");
194          const size_t numParticles = cfs.particles().size();
195          _weightedTotalChargedPartNum->fill(numParticles);
196        }
197      }
198
199      // charged particle distributions
200      if(_initialisedSpectra) {
201        const ChargedFinalState& cfs = apply<ChargedFinalState>(e, "CFS");
202        const size_t numParticles = cfs.particles().size();
203        _weightedTotalChargedPartNum->fill(numParticles);
204        const ParticlePair& beams = apply<Beam>(e, "Beams").beams();
205        const double meanBeamMom = ( beams.first.p3().mod() +
206                                     beams.second.p3().mod() ) / 2.0;
207        for (const Particle& p : cfs.particles()) {
208          const double xp = p.p3().mod()/meanBeamMom;
209          _h_xp->fill(xp);
210          const double logxp = -std::log(xp);
211          _h_xi->fill(logxp);
212          const double xe = p.E()/meanBeamMom;
213          _h_xe->fill(xe);
214          const double pTinT  = dot(p.p3(), thrust.thrustMajorAxis());
215          const double pToutT = dot(p.p3(), thrust.thrustMinorAxis());
216          _h_pTin->fill(fabs(pTinT/GeV));
217          if(_h_pTout) _h_pTout->fill(fabs(pToutT/GeV));
218          const double momT = dot(thrust.thrustAxis()        ,p.p3());
219          const double rapidityT = 0.5 * std::log((p.E() + momT) /
220                                                  (p.E() - momT));
221          _h_rapidityT->fill(fabs(rapidityT));
222          const double momS = dot(sphericity.sphericityAxis(),p.p3());
223          const double rapidityS = 0.5 * std::log((p.E() + momS) /
224                                                  (p.E() - momS));
225          _h_rapidityS->fill(fabs(rapidityS));
226        }
227      }
228    }
229
230    void finalize() {
231      if(!_initialisedJets && !_initialisedSpectra) return;
232
233      if (_initialisedJets) {
234        normalize(_h_thrust);
235        normalize(_h_heavyjetmass);
236        normalize(_h_totaljetbroadening);
237        normalize(_h_widejetbroadening);
238        normalize(_h_cparameter);
239        normalize(_h_thrustmajor);
240        normalize(_h_thrustminor);
241        normalize(_h_jetmassdifference);
242        normalize(_h_aplanarity);
243        if(_h_planarity) normalize(_h_planarity);
244        normalize(_h_oblateness);
245        normalize(_h_sphericity);
246
247        for (size_t n=0; n<6; ++n) {
248          scale(_h_R_Durham[n], 1./sumOfWeights());
249        }
250
251        for (size_t n = 0; n < 5; ++n) {
252          if (_h_y_Durham[n]) {
253            scale(_h_y_Durham[n], 1.0/sumOfWeights());
254          }
255        }
256      }
257
258      const double avgNumParts = dbl(*_weightedTotalChargedPartNum) / sumOfWeights();
259
260
261      for (auto& b : mult->bins()) {
262        if (inRange(sqrtS()/GeV, b.xMin(), b.xMax())) {
263          b.set(avgNumParts, 0.);
264        }
265      }
266
267      if (_initialisedSpectra) {
268        normalize(_h_xp, avgNumParts);
269        normalize(_h_xi, avgNumParts);
270        normalize(_h_xe, avgNumParts);
271        normalize(_h_pTin , avgNumParts);
272        if (_h_pTout) normalize(_h_pTout, avgNumParts);
273        normalize(_h_rapidityT, avgNumParts);
274        normalize(_h_rapidityS, avgNumParts);
275      }
276    }
277
278
279  private:
280
281    bool _initialisedJets = false;
282    bool _initialisedSpectra = false;
283
284    Estimate1DPtr mult;
285    Histo1DPtr _h_xp;
286    Histo1DPtr _h_xi;
287    Histo1DPtr _h_xe;
288    Histo1DPtr _h_pTin;
289    Histo1DPtr _h_pTout;
290    Histo1DPtr _h_rapidityT;
291    Histo1DPtr _h_rapidityS;
292    Histo1DPtr _h_thrust;
293    Histo1DPtr _h_heavyjetmass;
294    Histo1DPtr _h_totaljetbroadening;
295    Histo1DPtr _h_widejetbroadening;
296    Histo1DPtr _h_cparameter;
297    Histo1DPtr _h_thrustmajor;
298    Histo1DPtr _h_thrustminor;
299    Histo1DPtr _h_jetmassdifference;
300    Histo1DPtr _h_aplanarity;
301    Histo1DPtr _h_planarity;
302    Histo1DPtr _h_oblateness;
303    Histo1DPtr _h_sphericity;
304
305    Histo1DPtr _h_R_Durham[6];
306    Histo1DPtr _h_y_Durham[5];
307
308    CounterPtr _weightedTotalChargedPartNum;
309
310  };
311
312
313
314  RIVET_DECLARE_ALIASED_PLUGIN(ALEPH_2004_I636645, ALEPH_2004_S5765862);
315
316}