1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/UnstableParticles.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief baryon correlations
 10  class ARGUS_1990_I295621 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(ARGUS_1990_I295621);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22      declare(FinalState(Cuts::pid==2212), "FS");
 23      declare(UnstableParticles(), "UFS");
 24      // histos
 25      for (unsigned int ix=0;ix<2;++ix) {
 26        book(_c[ix],"TMP/c_"+toString(ix+1));
 27        book(_h_angle[ix], 9, 1, 1+ix);
 28      }
 29      for (unsigned int ix=0; ix<3; ++ix) {
 30        book(_n_cont[ix], 1+ix, 1, 1);
 31        book(_n_Ups [ix], 4+ix, 1, 1);
 32      }
 33    }
 34
 35    /// Recursively walk the decay tree to find decay products of @a p
 36    void findDecayProducts(const Particle& mother, Particles& protons, Particles& lam) {
 37      for (const Particle& p: mother.children()) {
 38        const int id = abs(p.pid());
 39        if (id == 2212) {
 40          protons.push_back(p);
 41        }
 42        else if (id == 3122) {
 43          lam.push_back(p);
 44        }
 45        if (!p.children().empty()) {
 46          findDecayProducts(p, protons, lam);
 47        }
 48      }
 49    }
 50
 51    void fillHistos(unsigned int imode, Particles& protons, Particles& lam, LorentzTransform& cms_boost) {
 52      Particles protons2;
 53      for (const Particle& p : protons) {
 54        const double modp = cms_boost.transform(p).p3().mod();
 55        if (modp<0.4 || modp>1.2) continue;
 56        if (imode==0) _n_cont[0]->fill("< 9.46"s);
 57        else          _n_Ups [0]->fill("9.460 + 10.023"s);
 58        protons2.push_back(p);
 59      }
 60      for (unsigned int ix=0; ix<protons2.size(); ++ix) {
 61        Vector3 axis1 = cms_boost(protons2[ix].mom()).p3().unit();
 62        if (protons2[ix].parents()[0].abspid()==3122) continue;
 63        for (unsigned int iy=ix+1;iy<protons2.size();++iy) {
 64          if (protons2[ix].pid()*protons2[iy].pid()<0) continue;
 65          if (protons2[iy].parents()[0].abspid()==3122) continue;
 66          Vector3 axis2 = cms_boost(protons2[iy].mom()).p3().unit();
 67          double cTheta = axis1.dot(axis2);
 68          if (imode==0) _n_cont[1]->fill("< 9.46"s);
 69          else         _n_Ups [1]->fill("9.460 + 10.023"s);
 70          _h_angle[imode]->fill(cTheta);
 71        }
 72      }
 73      for (unsigned int ix=0; ix<lam.size(); ++ix) {
 74        for (unsigned int iy=ix+1; iy<lam.size(); ++iy) {
 75          if (lam[ix].pid()*lam[iy].pid()<0) continue;
 76          if (imode==0) _n_cont[2]->fill("< 9.46"s);
 77          else          _n_Ups [2]->fill("9.460 + 10.023"s);
 78        }
 79      }
 80    }
 81
 82    /// Perform the per-event analysis
 83    void analyze(const Event& event) {
 84      // Find the upsilons
 85      // First in unstable final state
 86      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 87      Particles upsilons = ufs.particles(Cuts::pid==553 || Cuts::pid==100553);
 88      // continuum
 89      if (upsilons.empty()) {
 90        _c[0]->fill();
 91        Particles protons = apply<FinalState>(event,"FS").particles();
 92        Particles lam = ufs.particles(Cuts::abspid==3122);
 93        LorentzTransform cms_boost;
 94        fillHistos(0,protons,lam,cms_boost);
 95      }
 96      // found an upsilon
 97      else {
 98        for (const Particle& ups : upsilons) {
 99          _c[1]->fill();
100          Particles protons,lam;
101          findDecayProducts(ups, protons, lam);
102          LorentzTransform cms_boost;
103          if (ups.p3().mod() > 0.001) {
104            cms_boost = LorentzTransform::mkFrameTransformFromBeta(ups.mom().betaVec());
105          }
106          fillHistos(1,protons,lam,cms_boost);
107        }
108      }
109    }
110
111
112    /// Normalise histograms etc., after the run
113    void finalize() {
114      scale(_n_cont, 1.0/ *_c[0]);
115      scale(_n_Ups,  1.0/ *_c[1]);
116      scale(_h_angle[0], 1e4/ *_c[0]);
117      scale(_h_angle[1], 1e3/ *_c[1]);
118      BinnedEstimatePtr<string> tmp;
119      for (unsigned int ix=0;ix<2;++ix) {
120        book(tmp, 7+ix, 1, 1);
121        double val = _n_Ups[ix]->bin(1).sumW()/_n_cont[ix]->bin(1).sumW();
122        double err = val*sqrt( sqr(_n_Ups [ix]->bin(1).errW()/ _n_Ups[ix]->bin(1).sumW()) +
123                               sqr(_n_cont[ix]->bin(1).errW()/_n_cont[ix]->bin(1).sumW()));
124        tmp->bin(1).set(val,err);
125      }
126    }
127
128    /// @}
129
130
131    /// @name Histograms
132    /// @{
133    BinnedHistoPtr<string> _n_cont[3],_n_Ups[3];
134    Histo1DPtr _h_angle[2];
135    CounterPtr _c[2];
136
137    /// @}
138
139
140  };
141
142
143  RIVET_DECLARE_PLUGIN(ARGUS_1990_I295621);
144
145}