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#include "Rivet/Projections/Sphericity.hh"
  7
  8namespace Rivet {
  9
 10
 11  /// @brief Event shapes at Upsilon(1S)
 12  class ARGUS_1986_I227324 : public Analysis {
 13  public:
 14
 15    /// Constructor
 16    RIVET_DEFAULT_ANALYSIS_CTOR(ARGUS_1986_I227324);
 17
 18
 19    /// @name Analysis methods
 20    /// @{
 21
 22    /// Book histograms and initialise projections before the run
 23    void init() {
 24      // projections
 25      declare(UnstableParticles(), "UFS");
 26      declare(ChargedFinalState(), "CFS");
 27      const FinalState fs;
 28      declare(Thrust(fs)    ,"Thrust");
 29      declare(Sphericity(fs),"Sphericity");
 30      // histograms
 31      if(isCompatibleWithSqrtS(9.98*GeV,1e-2)) {
 32        book(_h_T_cont ,2, 1, 2);
 33        book(_h_S_cont ,1, 1, 2);
 34      }
 35      book(_h_T_Ups ,2, 1, 1);
 36      book(_h_S_Ups ,1, 1, 1);
 37    }
 38
 39    /// Recursively walk the decay tree to find the stable decay products of @a p
 40    void findDecayProducts(Particle mother, Particles& charged, Particles & neutral) {
 41      for(const Particle & p: mother.children()) {
 42	if(!p.children().empty())
 43	  findDecayProducts(p, charged, neutral);
 44	else {
 45	  if(isCharged(p))
 46	    charged.push_back(p);
 47	  else
 48	    neutral.push_back(p);
 49	}
 50      }
 51    }
 52
 53
 54    /// Perform the per-event analysis
 55    void analyze(const Event& event) {
 56      // Find the Upsilons among the unstables
 57      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 58      Particles upsilons = ufs.particles(Cuts::pid==553 or Cuts::pid==100553);
 59      if (upsilons.empty() && _h_T_cont) {
 60	Particles charged = apply<ChargedFinalState>(event, "CFS").particles();
 61	// at least 6 charged particles
 62	if(charged.size()<6) vetoEvent;
 63	// cut on high momentum particles
 64	unsigned int nHigh(0);
 65	for(const Particle & p : charged) {
 66	  if(p.momentum().p3().mod()>2.5) ++nHigh;
 67	}
 68	if(nHigh>1) vetoEvent;
 69        MSG_DEBUG("No Upsilons found => continuum event");
 70	Thrust thrust = apply<Thrust>(event, "Thrust");
 71	_h_T_cont->fill(thrust.thrust());
 72	Sphericity sphericity = apply<Sphericity>(event, "Sphericity");
 73	_h_S_cont->fill(sphericity.sphericity());
 74      }
 75      else {
 76        for (const Particle& ups : upsilons) {
 77          LorentzTransform boost;
 78          if (ups.p3().mod() > 1*MeV)
 79            boost = LorentzTransform::mkFrameTransformFromBeta(ups.momentum().betaVec());
 80          // Find the decay products we want
 81          Particles charged,neutral;
 82	  // 6 charged particles
 83          findDecayProducts(ups, charged, neutral);
 84	  if(charged.size()<6) continue;
 85	  // at most 1 |p|>2.5
 86	  vector<FourMomentum> mom;
 87	  mom.reserve(neutral.size()+charged.size());
 88	  unsigned int nHigh(0);
 89	  for(const Particle & p : charged) {
 90	    mom.push_back(boost.transform(p.momentum()));
 91	    if(mom.back().p3().mod()>2.5) ++nHigh;
 92	  }
 93	  if(nHigh>1) continue;
 94	  for(const Particle & p : neutral) {
 95	    mom.push_back(boost.transform(p.momentum()));
 96	  }
 97	  Thrust thrust;
 98	  thrust.calc(mom);
 99	  _h_T_Ups->fill(thrust.thrust());
100	  Sphericity sphericity;
101	  sphericity.calc(mom);
102	  _h_S_Ups->fill(sphericity.sphericity());
103	}
104      }
105    }
106
107
108    /// Normalise histograms etc., after the run
109    void finalize() {
110      if(_h_T_cont) {
111	normalize(_h_T_cont);
112	normalize(_h_S_cont);
113      }
114      if(_h_T_Ups->numEntries()!=0.) {
115	normalize(_h_T_Ups);
116	normalize(_h_S_Ups);
117      }
118    }
119
120    /// @}
121
122
123    /// @name Histograms
124    /// @{
125    Histo1DPtr _h_T_Ups,_h_T_cont;
126    Histo1DPtr _h_S_Ups,_h_S_cont;
127    /// @}
128
129
130  };
131
132
133  RIVET_DECLARE_PLUGIN(ARGUS_1986_I227324);
134
135}