1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/Thrust.hh"
  4#include "Rivet/Projections/Sphericity.hh"
  5#include "Rivet/Projections/ChargedFinalState.hh"
  6
  7namespace Rivet {
  8
  9
 10  class TASSO_1990_I294755 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(TASSO_1990_I294755);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22      const ChargedFinalState cfs(Cuts::pT >=  0.1/GeV);
 23      declare(cfs, "CFS");
 24
 25      // Thrust and sphericity
 26      declare(Thrust(cfs), "Thrust");
 27      declare(Sphericity(cfs), "Sphericity");
 28
 29      // Histos
 30      for (double eVal : allowedEnergies()) {
 31
 32        const string en = toString(int(eVal/MeV));
 33        if (isCompatibleWithSqrtS(eVal, 1e-02))  _sqs = en;
 34
 35        int offset = 0;
 36        switch (int(eVal+0.5)) {
 37          case 14:
 38            offset = 0;
 39            break;
 40          case 22:
 41            offset = 1;
 42            break;
 43          case 35:
 44            offset = 2;
 45            break;
 46          case 44:
 47            offset = 3;
 48            break;
 49        }
 50
 51        book(_h[en+"xp0"],        2, 1, 1+offset);
 52        book(_h[en+"xp1"],        3, 1, 1+offset);
 53        book(_h[en+"xi"],         4, 1, 1+offset);
 54        book(_h[en+"pT"],         5, 1, 1+offset);
 55        book(_n[en+"sphericity"], 6, 1, 1+offset);
 56        book(_n[en+"aplanarity"], 7, 1, 1+offset);
 57        book(_n[en+"thrust"],     8, 1, 1+offset);
 58        book(_c[en], "/TMP/_sumWPassed_"+en);
 59      }
 60      if (_sqs == "" && !merging()) {
 61        throw BeamError("Invalid beam energy for " + name() + "\n");
 62      }
 63    }
 64
 65
 66    /// Perform the per-event analysis
 67    void analyze(const Event& event) {
 68      const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");
 69
 70      // TASSO hadronic event selection TODO: move this into a trigger definition
 71      // See page 2 in publication
 72      // Condition 1)  --- require at least 5 (4) 'good' tracks
 73      int nch = cfs.particles().size();
 74      if ( (int(sqrtS()/GeV) > 27 && nch < 5) || (int(sqrtS()/GeV) <= 27 && nch < 4 ) ) {
 75        MSG_DEBUG("Failed # good tracks cut: " << nch);
 76        vetoEvent;
 77      }
 78      // Condition 2) ---
 79      // Condition 5) --- scalar momentum (not pT!!!) sum >= 0.265*s
 80      double momsum = 0.0;
 81      for (const Particle& p : cfs.particles()) {
 82        const double mom = p.p3().mod();
 83        momsum += mom;
 84      }
 85      if (momsum <=0.265 * sqrtS()/GeV) {
 86        MSG_DEBUG("Failed pTsum cut: " << momsum << " < " << 0.265 * sqrtS()/GeV);
 87        vetoEvent;
 88      }
 89
 90      // Raise counter for events that pass trigger conditions
 91      _c[_sqs]->fill();
 92
 93      const Thrust& thrust = apply<Thrust>(event, "Thrust");
 94      //const Vector3 & thrustAxis = thrust.thrustAxis ();
 95      //double theta = thrustAxis.theta();
 96      //if ( fabs(cos(theta)) >= 0.8 ) {
 97        //MSG_DEBUG("Failed thrust angle cut: " << fabs(cos(theta)));
 98        //vetoEvent;
 99      //}
100
101      const Sphericity& sphericity = apply<Sphericity>(event, "Sphericity");
102
103      // Fill histograms in order of appearance in paper
104      for (const Particle& p : cfs.particles()) {
105        // Get momentum and energy of each particle.
106        const Vector3 mom3 = p.p3();
107        // Scaled momenta.
108        const double mom = mom3.mod();
109        const double scaledMom = 2.*mom/sqrtS();
110        const double pTin = dot(mom3, sphericity.sphericityMajorAxis());
111        const double pTout = dot(mom3, sphericity.sphericityMinorAxis());
112        const double pT = sqrt(sqr(pTin)+sqr(pTout));
113        _h[_sqs+"xp0"]->fill(scaledMom);
114        _h[_sqs+"xp1"]->fill(scaledMom);
115        _h[_sqs+"xi"]->fill(-log(scaledMom));
116        _h[_sqs+"pT"]->fill(pT);
117      }
118      // event shapes
119      _n[_sqs+"sphericity"]->fill(sphericity.sphericity());
120      _n[_sqs+"aplanarity"]->fill(sphericity.aplanarity());
121      _n[_sqs+"thrust"]->fill(thrust.thrust());
122    }
123
124
125    /// Normalise histograms etc., after the run
126    void finalize() {
127      for (auto& item : _h) {
128        const string en = item.first.substr(0,5);
129        if (_c[en]->sumW())  scale(item.second, 1./ *_c[en]);
130      }
131      normalize(_n);
132    }
133
134    /// @}
135
136
137  private:
138
139    /// @name Histograms
140    /// @{
141    map<string,Histo1DPtr> _h, _n;
142    map<string,CounterPtr> _c;
143
144    string _sqs = "";
145    /// @}
146
147  };
148
149
150
151  RIVET_DECLARE_ALIASED_PLUGIN(TASSO_1990_I294755, TASSO_1990_S2148048);
152
153}