1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief anti-deuteron spectrum in upslion decays
  9  class BABAR_2014_I1286317 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(BABAR_2014_I1286317);
 14
 15
 16    /// @name Analysis methods
 17    /// @{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21      // Initialise and register projections
 22      declare(UnstableParticles(), "UFS");
 23      // histograms
 24      book(_h_p[0],6,1,1);
 25      book(_h_p[1],6,2,1);
 26      book(_h_p[2],6,3,1);
 27      book(_h_p[3],6,4,1);
 28      book(_h_r[0],1,1,1);
 29      book(_h_r[1],2,1,1);
 30      book(_h_r[2],3,1,1);
 31      book(_h_r[3],5,1,1);
 32      book(_h_r[4],4,1,1);
 33      book(_w[0],"TMP/w_0");
 34      book(_w[1],"TMP/w_1");
 35      book(_w[2],"TMP/w_2");
 36      book(_w[3],"TMP/w_3");
 37    }
 38
 39    /// Recursively walk the decay tree to find decay products of @a p
 40    void findDecayProducts(Particle mother, Particles& deut) {
 41      for(const Particle & p: mother.children()) {
 42	if(p.pid() == _did) {
 43	  deut.push_back(p);
 44	}
 45	else if(!p.children().empty())
 46	  findDecayProducts(p, deut);
 47      }
 48    }
 49
 50    /// Perform the per-event analysis
 51    void analyze(const Event& event) {
 52      // find upsilon states
 53      UnstableParticles  ufs = apply<UnstableParticles>(event, "UFS");
 54      Particles ups = ufs.particles(Cuts::pid==553||Cuts::pid==100553||Cuts::pid==200553);
 55      // none, then continuum event
 56      if(ups.empty()) {
 57	Particles deut = ufs.particles(Cuts::pid==_did);
 58	_w[3]->fill();
 59	for(const Particle& p : deut) {
 60	  double mom = p.momentum().p3().mod();
 61	  _h_p[3]->fill(mom);
 62	  _h_r[3]->fill("10.58"s);
 63	  _h_r[4]->fill("10.58"s);
 64	}
 65      }
 66      // upsilon decays
 67      else {
 68	for(const Particle & Y : ups ) {
 69	  unsigned int ihist=2;
 70	  if(Y.pid()==100553) {
 71	    ihist=1;
 72	  }
 73	  else if(Y.pid()==200553) {
 74	    ihist=0;
 75	  }
 76	  _w  [ihist]->fill();
 77	  Particles deut;
 78	  findDecayProducts(Y, deut);
 79	  if(deut.empty()) continue;
 80	  LorentzTransform boost;
 81	  if (Y.p3().mod() > 1*MeV)
 82	    boost = LorentzTransform::mkFrameTransformFromBeta(Y.momentum().betaVec());
 83	  for(const Particle& p : deut) {
 84	    double mom = boost.transform(p.momentum()).p3().mod();
 85	    _h_r[ihist]->fill("10.58"s);
 86	    _h_p[ihist]->fill(mom);
 87	    _w  [ihist]->fill();
 88	  }
 89	}
 90      }
 91    }
 92
 93
 94    /// Normalise histograms etc., after the run
 95    void finalize() {
 96      // upsilon decays
 97      for(unsigned int ix=0;ix<3;++ix) {
 98	if(_w[ix]->effNumEntries()<=0.) continue;
 99	scale(_h_p[ix],1e6/ *_w[ix]);
100	scale(_h_r[ix],1./ *_w[ix]);
101	Histo1DPtr _h_p[4],_h_r[5];
102	CounterPtr _w[4];
103      }
104      // continuum
105      if(_w[3]->effNumEntries()>0.) {
106	scale(_h_p[3], crossSection()/sumOfWeights()/femtobarn);
107	scale(_h_r[4], crossSection()/sumOfWeights()/femtobarn);
108	scale(_h_r[3],1./ *_w[3]);
109      }
110    }
111
112    /// @}
113
114
115    /// @name Histograms
116    /// @{
117    Histo1DPtr _h_p[4];
118    BinnedHistoPtr<string> _h_r[5];
119    CounterPtr _w[4];
120
121    // deuteron id code
122    static const int _did = -1000010020;
123    /// @}
124
125
126  };
127
128
129  RIVET_DECLARE_PLUGIN(BABAR_2014_I1286317);
130
131}