Rivet analyses

// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/UnstableParticles.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/Thrust.hh"

namespace Rivet {


  /// @brief DELPHI W decay analysis
  class DELPHI_2001_I526164 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(DELPHI_2001_I526164);


    /// @name Analysis methods
    ///@{

    /// Book histograms and initialise projections before the run
    void init() {
      // projections
      declare(UnstableParticles(),"UFS");
      declare(FinalState(),"FS");
      // Histogram booking
      for (size_t ih=1; ih<7; ++ih) {
        book(_e[ih], (ih==1)? 1 : 2, 1, (ih==1)? 1 : ih-1);
        if (ih<3)  book(_e[99+ih], 3, 1, 2*ih-1);
      }
      for (double eVal : allowedEnergies()) {
        const string en = toString(round(eVal));
        if (isCompatibleWithSqrtS(eVal))  _sqs = en;
        for (size_t ih=0; ih<7; ++ih) {
          book(_c[en+toString(ih)], "_n_qq_"s+toString(ih)+"_"s+en);
        }
        if (en == "183"s) {
          book(_q[en]["K0"],  16, 1, 1);
          book(_q[en]["Lam"], 16, 1, 3);
          book(_d[en+"p"],   8, 1, 1);
          book(_d[en+"xi"], 10, 1, 3);
          book(_d[en+"pT"], 12, 1, 3);
          for (size_t ih=0; ih<2; ++ih) {
            const string iw = toString(ih);
            book(_h[en+iw]["p_charged"], 7, 1, 2-ih);
            book(_h[en+iw]["p_chargedB"], "/TMP/h_7_1_"  +toString(2-ih), _d[en+"p"].binning().edges<0>());
            book(_h[en+iw]["xi_charged"], 10, 1, 2-ih);
            book(_h[en+iw]["xi_chargedB"], "/TMP/h_10_1_"+toString(2-ih), _d[en+"xi"].binning().edges<0>());
            book(_h[en+iw]["pT_charged"], 12, 1, 2-ih);
            book(_h[en+iw]["pT_chargedB"], "/TMP/h_12_1_"+toString(2-ih), _d[en+"pT"].binning().edges<0>());
          }
        }
        else if (en == "189"s) {
          book(_q[en]["K0"],  16, 1, 2);
          book(_q[en]["Lam"], 16, 1, 4);
          book(_d[en+"p"],   6, 1, 1);
          book(_d[en+"xi"],  9, 1, 3);
          book(_d[en+"pT"], 11, 1, 3);
          for (size_t ih=0; ih<2; ++ih) {
            const string iw = toString(ih);
            book(_h[en+iw]["p_charged"], 5, 1, 2-ih);
            book(_h[en+iw]["p_chargedB"], "/TMP/h_5_1_"  +toString(2-ih), _d[en+"p"].binning().edges<0>());
            book(_h[en+iw]["xi_charged"], 9, 1, 2-ih);
            book(_h[en+iw]["xi_chargedB"], "/TMP/h_9_1_" +toString(2-ih), _d[en+"xi"].binning().edges<0>());
            book(_h[en+iw]["pT_charged"], 11, 1, 2-ih);
            book(_h[en+iw]["pT_chargedB"], "/TMP/h_11_1_"+toString(2-ih), _d[en+"pT"].binning().edges<0>());
            for (size_t iy=0; iy<4; ++iy) {
              book(_h[en+iw]["xi_ident"+toString(iy)], 13+ih, 1, 1+iy);
              if (!ih) book(_h[en+iw]["xi_ident2"+toString(iy)], 15, 1, 1+iy);
              if (iy)  book(_e[10+2*iy-ih], 4, 1, 2*iy-ih);
            }
          }
        }
        for (size_t ix=0; ix<2; ++ix) {
          for (size_t iy=0; iy<5; ++iy) {
            const string suff = toString(ix)+toString(iy);
            book(_c[en+suff], "_n_WW_"s+suff+"_"s+en);
          }
        }
      }
      raiseBeamErrorIf(_sqs.empty());
      isWW = _sqs == "183"s || _sqs == "189"s;
    }

    /// Perform the per-event analysis
    void analyze(const Event& event) {
      const double mW=80.379*GeV;
      multimap<Particle,Particles,ParticleOrdering> Wbosons;
      // loop over final state particles to find W's
      const Particles& finalState = apply<FinalState>(event,"FS").particles();
      // loop over FS particles
      for (const Particle & p : finalState) {
        Particle parent=p;
        while (!parent.parents().empty()) {
          parent=parent.parents()[0];
          if (parent.abspid()==24 && parent.mass()>20.) break;
        }
        if (parent.abspid()!=24) continue;
        // find those which came from W
        bool found=false;
        for (auto& W : Wbosons) {
          // W already in list add particle to its decay products
          if (fuzzyEquals(W.first.mom(), parent.mom())) {
            W.second.push_back(p);
            found=true;
            break;
          }
        }
        if (!found) {
          // check W not child
          bool Wchild=false;
          for (const Particle& child : parent.children()) {
            if (child.abspid()==24) {
              Wchild=true;
              break;
            }
          }
          // add to list
          if (!Wchild) {
            Particles temp = {p};
            Wbosons.insert(make_pair(parent,temp));
          }
        }
      }
      // no W's => q qbar event
      if (Wbosons.empty()) {
        _c[_sqs+"0"s]->fill();
        UnstableParticles ufs=apply<UnstableParticles>(event,"UFS");
        for (const Particle & p : ufs.particles()) {
          if(p.abspid()==PID::PIPLUS) {
            _c[_sqs+"1"s]->fill();
            _c[_sqs+"2"s]->fill();
          }
          else if(p.abspid()==PID::KPLUS) {
            _c[_sqs+"1"s]->fill();
            _c[_sqs+"3"s]->fill();
          }
          else if(p.abspid()==PID::PROTON) {
            _c[_sqs+"1"s]->fill();
            _c[_sqs+"5"s]->fill();
          }
          else if(p.abspid()==130 || p.abspid()==310 ) {
            _c[_sqs+"4"s]->fill();
            if (isWW) _q[_sqs]["K0"]->fill(-log(2.*p.p3().mod()/sqrtS()));
          }
          else if(p.abspid()==PID::LAMBDA ) {
            _c[_sqs+"6"s]->fill();
            if (isWW) _q[_sqs]["Lam"]->fill(-log(2.*p.p3().mod()/sqrtS()));
          }
        }
      }
      else if (Wbosons.size()==2) {
        bool leptonic[2] = {false,false};
        unsigned int iboson=0;
        for (auto& W : Wbosons) {
          for (const Particle& child : W.first.children()) {
            // find lepton bosons
            if (child.abspid()==11 || child.abspid()==13) {
              leptonic[iboson]=true;
              break;
            }
            // veto events with W-> tau decays
            else if (child.abspid()==15) vetoEvent;
          }
          ++iboson;
        }
        // only fully hadronic or semi-leptonic events
        if (leptonic[0] && leptonic[1]) vetoEvent;
        // 0 SL 1 hadronic
        const string iw = (leptonic[0] || leptonic[1])? "0"s : "1"s;
        _c[_sqs+iw+"0"s]->fill();
        Particles particles;
        iboson=0;
        for (auto& W : Wbosons) {
          if (!leptonic[iboson]) {
            particles.insert(particles.end(),W.second.begin(),W.second.end());
          }
          ++iboson;
        }
        // calculate thrust
        Thrust thrust;
        thrust.calc(particles);
        // type of event
        // loop over particles and fill histos
        for (const Particle& p : particles) {
          if (!PID::isCharged(p.pid())) continue;
          // Get momentum of each particle.
          const Vector3 mom3 = p.p3();
          // Scaled momenta.
          const double mom = mom3.mod();
          const double xiW = -log(2.*mom/mW);
          const double xiL = -log(2.*mom/sqrtS());
          // Get momenta components w.r.t. thrust
          const double pTinT = dot(mom3, thrust.thrustMajorAxis());
          const double pToutT = dot(mom3, thrust.thrustMinorAxis());
          double pT = sqrt(sqr(pTinT)+sqr(pToutT));
          // fill charged particle hists
          if (isWW) {
            _h[_sqs+iw]["p_charged"s]->fill(mom);
            _h[_sqs+iw]["xi_charged"s]->fill(xiL);
            _h[_sqs+iw]["pT_charged"s]->fill(pT/GeV);
            _h[_sqs+iw]["p_chargedB"s]->fill(mom);
            _h[_sqs+iw]["xi_chargedB"s]->fill(xiL);
            _h[_sqs+iw]["pT_chargedB"s]->fill(pT/GeV);
          }
          // and identified particles
          _c[_sqs+iw+"1"s]->fill();
          if (p.abspid()==PID::PIPLUS)       _c[_sqs+iw+"2"s]->fill();
          else if (p.abspid()==PID::KPLUS)   _c[_sqs+iw+"3"s]->fill();
          else if (p.abspid()==PID::PROTON)  _c[_sqs+iw+"4"s]->fill();
          if (_sqs == "189"s) {
            _h[_sqs+iw]["xi_ident0"]->fill(xiW);
            if (p.abspid()==PID::PIPLUS)       _h[_sqs+iw]["xi_ident1"]->fill(xiW);
            else if (p.abspid()==PID::KPLUS)   _h[_sqs+iw]["xi_ident2"]->fill(xiW);
            else if (p.abspid()==PID::PROTON)  _h[_sqs+iw]["xi_ident3"]->fill(xiW);
          }
        }
        // boosted specta in W rest frame
        if (iw=="0"s && _sqs == "189"s) {
          iboson=0;
          for (auto& W : Wbosons) {
            ++iboson;
            if (leptonic[iboson-1]) continue;
            // boost to rest frame
            LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(W.first.momentum().betaVec());
            FourMomentum psum;
            for (const Particle& p : W.second) psum+=p.mom();
            // spectra
            for (const Particle& p : W.second) {
              if (!PID::isCharged(p.pid())) continue;
              // Get momentum of each particle.
              FourMomentum phad = boost.transform(p.mom());
              const Vector3 mom3 = phad.p3();
              // Scaled momenta.
              const double mom = mom3.mod();
              const double scaledMom = 2.*mom/mW;
              const double xi = -log(scaledMom);
              // /identified particle spectra
              _h[_sqs+iw]["xi_ident20"]->fill(xi);
              if (p.abspid()==PID::PIPLUS)       _h[_sqs+iw]["xi_ident21"]->fill(xi);
              else if (p.abspid()==PID::KPLUS)   _h[_sqs+iw]["xi_ident22"]->fill(xi);
              else if (p.abspid()==PID::PROTON)  _h[_sqs+iw]["xi_ident23"]->fill(xi);
            }
          }
        }
      }
    }


    /// Normalise histograms etc., after the run
    void finalize() {
      // qq
      scale(_h, crossSectionPerEvent());
      scale(_q, crossSectionPerEvent());
      scale(_c, crossSectionPerEvent());
      for (double eVal : allowedEnergies()) {
        const int edge = round(eVal);
        const string en = toString(edge);
        if (_c[en+"0"s]->effNumEntries()) {
          scale(_q[en],  1.0/ *_c[en+"0"s]);
          for (size_t ih=1; ih<7; ++ih) {
            if (!_c[en+toString(ih)]->effNumEntries())  continue;
            // hist for axis
            for (auto& b : _e[ih]->bins()) {
              if (edge == b.xEdge()) {
                b = *_c[en+toString(ih)] / *_c[en+"0"s];
              }
            }
          }
        }
        if (edge != 183 && edge != 189)  continue;
        for (size_t ih=0; ih<2; ++ih) {
          const string iw = toString(ih);
          if (!_c[en+iw+"0"s]->effNumEntries())  continue;
          const double sf = 1.0/ _c[en+iw+"0"s]->val();
          scale(_h[en+iw], sf);

          // charged multis
          Estimate0D ratio = *_c[en+toString(iw)+"1"s]/ *_c[en+toString(iw)+"0"s];
          _e[100+ih]->binAt(edge).set(ratio.val(), ratio.errPos());

          if (en != "189"s)  continue;
          for (size_t iy=1; iy<4; ++iy) {
            Estimate0D R = *_c[en+toString(iw)+toString(iy+1)] / *_c[en+toString(iw)+"0"s];
            for (auto& b : _e[10+2*iy-ih]->bins()) {
              b.set(R.val(), R.errPos());
            }
          }
        }

        // difference histos
        for (const string& obs : {"p"s, "xi"s, "pT"s}) {
          YODA::Estimate1D htemp = _h[en+"0"][obs+"_chargedB"]->mkEstimate();
          htemp.scale(2.0);
          const string hpath = _d[en+obs]->path();
          *_d[en+obs] = *_h[en+"1"][obs+"_chargedB"] - htemp;
          _d[en+obs]->setPath(hpath);
        }

      }
    }

    struct ParticleOrdering {
      bool operator()(const Particle& p1, const Particle& p2) const {
        return p1.pid() > p2.pid();
      }
    };

    ///@}


    /// @name Histograms
    ///@{
    // qqbar
    map<string,map<string,Histo1DPtr>> _h, _q;
    map<size_t,BinnedEstimatePtr<int>> _e;
    map<string,Estimate1DPtr> _d;
    map<string,CounterPtr> _c;
    string _sqs = "";
    bool isWW;
    ///@}


  };


  RIVET_DECLARE_PLUGIN(DELPHI_2001_I526164);

}