1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4#include "Rivet/Projections/ChargedFinalState.hh"
  5#include "Rivet/Projections/Thrust.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// @brief Thrust like variable at Upsilon(1S,2S)
 11  class LENA_1981_I164397 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(LENA_1981_I164397);
 16
 17
 18    /// @name Analysis methods
 19    /// @{
 20
 21    /// Book histograms and initialise projections before the run
 22    void init() {
 23      declare(UnstableParticles(), "UFS");
 24      declare(ChargedFinalState(), "FS");
 25      book(_mult,3,1,1);
 26      _ecms = "OTHER";
 27      if (inRange(sqrtS(),7.35,7.49)) {
 28        _ecms = "7.35 - 7.49"s;
 29      }
 30      else if (inRange(sqrtS(),8.629,9.142)) {
 31        _ecms = "8.629 - 9.142"s;
 32      }
 33      else if (inRange(sqrtS(),9.15,9.41)) {
 34        _ecms = "9.15 - 9.41"s;
 35      }
 36      else if(isCompatibleWithSqrtS(9.5149*GeV,1e-2)) {
 37        _ecms = "9.5149";
 38        book(_hist_T_cont ,4, 1, 1);
 39      }
 40      else if(isCompatibleWithSqrtS(9.9903*GeV,1e-2)) {
 41        _ecms = "9.9903";
 42        book(_hist_T_cont ,4, 1, 2);
 43      }
 44      book(_hist_T_Ups1 ,4, 1, 3);
 45      book(_hist_T_Ups2 ,4, 1, 4);
 46      book(_weightSum_cont, "TMP/weightSum_cont");
 47      book(_weightSum_Ups1, "TMP/weightSum_Ups1");
 48      book(_weightSum_Ups2, "TMP/weightSum_Ups2");
 49    }
 50
 51
 52    /// Recursively walk the decay tree to find the charged decay products of @a p
 53    void findDecayProducts(Particle mother, Particles& charged) {
 54      for(const Particle & p: mother.children()) {
 55        const int id = p.pid();
 56	if(!p.children().empty())
 57	  findDecayProducts(p, charged);
 58	else if(PID::isCharged(id))
 59	  charged.push_back(p);
 60      }
 61    }
 62
 63    // defn of thrust in paper used just the direction
 64    double thrustPrime(const LorentzTransform & boost, const Particles & particles) {
 65      vector<Vector3> vecs;
 66      for(const Particle & p : particles) {
 67	vecs.push_back(boost.transform(p.momentum()).p3().unit());
 68      }
 69      Thrust thrust;
 70      thrust.calc(vecs);
 71      return thrust.thrust();
 72    }
 73
 74    /// Perform the per-event analysis
 75    void analyze(const Event& event) {
 76      // Find the Upsilons among the unstables
 77      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 78      Particles upsilons = ufs.particles(Cuts::pid==553 or Cuts::pid==100553);
 79      if (upsilons.empty()) {
 80        MSG_DEBUG("No Upsilons found => continuum event");
 81        _weightSum_cont->fill();
 82	Particles cfs = apply<ChargedFinalState>(event, "FS").particles();
 83	_mult->fill(_ecms,cfs.size());
 84	if(_hist_T_cont) {
 85          LorentzTransform boost;
 86	  _hist_T_cont->fill(thrustPrime(boost,cfs));
 87	}
 88      }
 89      // Upsilon(s) found
 90      else {
 91        for (const Particle& ups : upsilons) {
 92          const int parentId = ups.pid();
 93          Particles charged;
 94	  // boost to rest frame (if required)
 95          LorentzTransform boost;
 96          if (ups.p3().mod() > 1*MeV)
 97            boost = LorentzTransform::mkFrameTransformFromBeta(ups.momentum().betaVec());
 98          // Find the decay products we want
 99          findDecayProducts(ups, charged);
100	  if(parentId==553) {
101	    _weightSum_Ups1->fill();
102	    _mult->fill("9.4624"s,charged.size());
103	    _hist_T_Ups1->fill(thrustPrime(boost,charged));
104	  }
105	  else {
106	    _weightSum_Ups2->fill();
107	    _mult->fill("10.0148"s,charged.size());
108	    _hist_T_Ups2->fill(thrustPrime(boost,charged));
109	  }
110	}
111      }
112
113    }
114
115
116    /// Normalise histograms etc., after the run
117    void finalize() {
118      // charged particle multiplicity
119      if ( _weightSum_cont->val()>0. && _hist_T_cont )
120        scale(_hist_T_cont,1./ *_weightSum_cont );
121      if(_weightSum_Ups1->val()>0. )
122         scale(_hist_T_Ups1,1./ *_weightSum_Ups1 );
123      if(_weightSum_Ups2->val()>0. )
124        scale(_hist_T_Ups2,1./ *_weightSum_Ups2 );
125    }
126
127    /// @}
128
129
130    /// @name Histograms
131    /// @{
132    BinnedProfilePtr<string> _mult;
133    CounterPtr _weightSum_cont, _weightSum_Ups1, _weightSum_Ups2;
134    Histo1DPtr _hist_T_cont,_hist_T_Ups1,_hist_T_Ups2;
135    string _ecms;
136    /// @}
137
138
139  };
140
141
142  RIVET_DECLARE_PLUGIN(LENA_1981_I164397);
143
144}