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      
 26      book(_weightSum_cont, "TMP/weightSum_cont");
 27      book(_weightSum_Ups1, "TMP/weightSum_Ups1");
 28      book(_weightSum_Ups2, "TMP/weightSum_Ups2");
 29      book(_charge_cont, "TMP/charge_cont");
 30      book(_charge_Ups1, "TMP/charge_Ups1");
 31      book(_charge_Ups2, "TMP/charge_Ups2");
 32
 33      if(isCompatibleWithSqrtS(9.5149,1e-2)) {
 34        book(_hist_T_cont ,4, 1, 1);
 35      }
 36      else if(isCompatibleWithSqrtS(9.9903,1e-2)) {
 37        book(_hist_T_cont ,4, 1, 2);
 38      }
 39      book(_hist_T_Ups1 ,4, 1, 3);
 40      book(_hist_T_Ups2 ,4, 1, 4);
 41    }
 42
 43    /// Recursively walk the decay tree to find the charged decay products of @a p
 44    void findDecayProducts(Particle mother, Particles& charged) {
 45      for(const Particle & p: mother.children()) {
 46        const int id = p.pid();
 47	if(!p.children().empty())
 48	  findDecayProducts(p, charged);
 49	else if(PID::isCharged(id))
 50	  charged.push_back(p);
 51      }
 52    }
 53
 54    // defn of thrust in paper used just the direction
 55    double thrustPrime(const LorentzTransform & boost, const Particles & particles) {
 56      vector<Vector3> vecs;
 57      for(const Particle & p : particles) {
 58	vecs.push_back(boost.transform(p.momentum()).p3().unit());
 59      }
 60      Thrust thrust;
 61      thrust.calc(vecs);
 62      return thrust.thrust();
 63    }
 64
 65    /// Perform the per-event analysis
 66    void analyze(const Event& event) {
 67      // Find the Upsilons among the unstables
 68      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 69      Particles upsilons = ufs.particles(Cuts::pid==553 or Cuts::pid==100553);
 70      if (upsilons.empty()) { 
 71        MSG_DEBUG("No Upsilons found => continuum event");
 72        _weightSum_cont->fill();
 73	Particles cfs = apply<ChargedFinalState>(event, "FS").particles();
 74	_charge_cont->fill(cfs.size());
 75	if(_hist_T_cont) {
 76          LorentzTransform boost;
 77	  _hist_T_cont->fill(thrustPrime(boost,cfs));
 78	}
 79      }
 80      // Upsilon(s) found
 81      else {
 82        for (const Particle& ups : upsilons) {
 83          const int parentId = ups.pid();
 84          Particles charged;
 85	  // boost to rest frame (if required)
 86          LorentzTransform boost;
 87          if (ups.p3().mod() > 1*MeV)
 88            boost = LorentzTransform::mkFrameTransformFromBeta(ups.momentum().betaVec());
 89          // Find the decay products we want
 90          findDecayProducts(ups, charged);
 91	  if(parentId==553) {
 92	    _weightSum_Ups1->fill();
 93	    _charge_Ups1->fill(charged.size());
 94	    _hist_T_Ups1->fill(thrustPrime(boost,charged));
 95	  }
 96	  else {
 97	    _weightSum_Ups2->fill();
 98	    _charge_Ups2->fill(charged.size());
 99	    _hist_T_Ups2->fill(thrustPrime(boost,charged));
100	  }
101	}
102      }
103
104    }
105
106
107    /// Normalise histograms etc., after the run
108    void finalize() {
109      // charged particle multiplicity
110      if(_weightSum_cont->val()>0. ) {
111	scale(_charge_cont,1./ *_weightSum_cont );
112	if(_hist_T_cont) scale(_hist_T_cont,1./ *_weightSum_cont );
113      }
114      if(_weightSum_Ups1->val()>0. ) {
115	scale(_charge_Ups1,1./ *_weightSum_Ups1 );
116	 scale(_hist_T_Ups1,1./ *_weightSum_Ups1 );
117      }
118      if(_weightSum_Ups2->val()>0. ) {
119	scale(_charge_Ups2,1./ *_weightSum_Ups2 );
120	scale(_hist_T_Ups2,1./ *_weightSum_Ups2 );
121      }
122      Scatter2D tempScat(refData(3, 1, 1));
123      Scatter2DPtr _mult;
124      book(_mult, 3, 1, 1);
125      for (size_t b = 0; b < tempScat.numPoints(); b++) {
126        const double x  = tempScat.point(b).x();
127        pair<double,double> ex = tempScat.point(b).xErrs();
128        pair<double,double> ex2 = ex;
129        if(ex2.first ==0.) ex2. first=0.02;
130        if(ex2.second==0.) ex2.second=0.02;
131	// Upsilon 1S
132	if(b==3) {
133	  if (_weightSum_Ups1->val()>0.) {
134	    _mult->addPoint(x, _charge_Ups1->val(), ex, make_pair(_charge_Ups1->err(),_charge_Ups1->err()));
135	  }
136	  else {
137	    _mult->addPoint(x, 0., ex, make_pair(0.,.0));
138	  }
139	}
140	// Upsilon 2S
141	else if(b==6) {
142	  if (_weightSum_Ups2->val()>0.) {
143	    _mult->addPoint(x, _charge_Ups2->val(), ex, make_pair(_charge_Ups2->err(),_charge_Ups2->err()));
144	  }
145	  else {
146	    _mult->addPoint(x, 0., ex, make_pair(0.,.0));
147	  }
148	}
149	else {
150	  if (inRange(sqrtS()/GeV, x-ex2.first, x+ex2.second) && _weightSum_cont->val()>0.) {
151	    _mult->addPoint(x, _charge_cont->val(), ex, make_pair(_charge_cont->err(),_charge_cont->err()));
152	  }
153	  else {
154	    _mult->addPoint(x, 0., ex, make_pair(0.,.0));
155	  }
156	}
157      }
158    }
159
160    //@}
161
162
163    /// @name Histograms
164    //@{
165    CounterPtr _weightSum_cont, _weightSum_Ups1, _weightSum_Ups2;
166    CounterPtr _charge_cont, _charge_Ups1, _charge_Ups2;
167    Histo1DPtr _hist_T_cont,_hist_T_Ups1,_hist_T_Ups2;
168    //@}
169
170
171  };
172
173
174  RIVET_DECLARE_PLUGIN(LENA_1981_I164397);
175
176}