Rivet analyses

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


namespace {

  unsigned int factorial(unsigned int n) {
    if(n<=1) return 0;
    unsigned int out=1;
    for (unsigned int i=n;i>0;--i) out *= i;
    return out;
  }

  double rapidityAxis(Rivet::Vector4 momentum, Rivet::Vector3 axis) {
    // assume axis is unit vector
    double plong = momentum.vector3().dot(axis);
    return 0.5*log((momentum.t()+plong)/(momentum.t()-plong));
  }

  double cosThetaStar(const Rivet::Particle & p1, const Rivet::Particle & p2, Rivet::Vector3 axis) {
    Rivet::FourMomentum psum = p1.momentum()+p2.momentum();
    Rivet::LorentzTransform trans(Rivet::LorentzTransform::mkObjTransformFromBeta(-psum.betaVec()));
    Rivet::FourMomentum m1 = trans.transform(p1.momentum());
    Rivet::FourVector a4;
    a4.setX(axis.x());
    a4.setY(axis.y());
    a4.setZ(axis.z());
    a4.setT(1.);
    a4 = trans.transform(a4);
    return fabs(m1.vector3().dot(a4.vector3())/m1.vector3().mod()/a4.vector3().mod());
  }

}

namespace Rivet {


  /// @brief lambda anti-lambda correlations
  class OPAL_2000_I474010 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(OPAL_2000_I474010);


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

    /// Book histograms and initialise projections before the run
    void init() {
      // projections
      const FinalState fs;
      declare(fs, "FS");
      FastJets durhamjets(fs, JetAlg::DURHAM, 0.7);
      durhamjets.useInvisibles(JetInvisibles::ALL);
      declare(durhamjets, "DurhamJets");
      const Thrust thrust(fs);
      declare(thrust, "Thrust");
      declare(UnstableParticles(), "UFS");
      // histograms
      book(_histLLBar    ,  1, 1, 1);
      book(_histLL       ,  1, 1, 2);
      book(_histLLBarCorr,  1, 1, 3);
      book(_histLLBar_2jet    ,  2, 1, 1);
      book(_histLL_2jet       ,  2, 1, 2);
      book(_histLLBarCorr_2jet,  2, 1, 3);
      book(_histLLBar_3jet    ,  3, 1, 1);
      book(_histLL_3jet       ,  3, 1, 2);
      book(_histLLBarCorr_3jet,  3, 1, 3);
      book(_histdcosTheta,  4, 1, 1);
      book(_histdy       ,  5, 1, 1);
      // weights
      book(_weight2Jet,"/TMP/W2Jet");
      book(_weight3Jet,"/TMP/W3Jet");

    }


    /// Perform the per-event analysis
    void analyze(const Event& event) {
      // First, veto on leptonic events by requiring at least 4 charged FS particles
      const FinalState& fs = apply<FinalState>(event, "FS");
      const size_t numParticles = fs.particles().size();

      // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
      if (numParticles < 2) {
        MSG_DEBUG("Failed leptonic event cut");
        vetoEvent;
      }
      MSG_DEBUG("Passed leptonic event cut");

      const Thrust& thrust = apply<Thrust>(event, "Thrust");
      Vector3 axis = thrust.thrustAxis();

      // Final state of unstable particles to get particle spectra
      const UnstableParticles & ufs = apply<UnstableParticles>(event, "UFS");

      Particles lambda,lambdaBar;
      for (const Particle& p : ufs.particles()) {
        const int id = p.pid();
        if(id==3122) lambda.push_back(p);
        else if(id==-3122) lambdaBar.push_back(p);
      }

      // get the number of jets
      const FastJets& durjet = apply<FastJets>(event, "DurhamJets");

      unsigned int njet = durjet.clusterSeq()->n_exclusive_jets_ycut(0.005);
      if(njet==2)  _weight2Jet->fill();
      else if(njet==3)  _weight3Jet->fill();

      // need a pair of baryons
      if(lambda.size()+lambdaBar.size()<2)  vetoEvent;

      // multiplicities
      int nsame = (factorial(lambda.size())+factorial(lambdaBar.size()))/2;
      int ndiff = lambda.size()*lambdaBar.size();
      _histLLBar    ->fill(Ecm, ndiff);
      _histLL       ->fill(Ecm, nsame);
      _histLLBarCorr->fill(Ecm, (ndiff-nsame));

      // multiplicities for different numbers of jets
      if(njet==2) {
        _histLLBar_2jet    ->fill(Ecm, ndiff);
        _histLL_2jet       ->fill(Ecm, nsame);
        _histLLBarCorr_2jet->fill(Ecm, (ndiff-nsame));
      }
      else if(njet==3) {
        _histLLBar_3jet    ->fill(Ecm, ndiff);
        _histLL_3jet       ->fill(Ecm, nsame);
        _histLLBarCorr_3jet->fill(Ecm, (ndiff-nsame));
      }
      // uncorrelated lambda
      for (unsigned int ix=0;ix<lambda.size();++ix) {
        double y1 = rapidityAxis(lambda[ix].momentum(),axis);
        for (unsigned int iy=ix+1;iy<lambda.size();++iy) {
          double y2 = rapidityAxis(lambda[iy].momentum(),axis);
          if(njet==2) _histdy->fill(fabs(y1-y2),-1.0);
          double ctheta = cosThetaStar(lambda[ix],lambda[iy],axis);
          _histdcosTheta->fill(ctheta,-1.0);
        }
      }
      // uncorrelated lambdabar
      for (unsigned int ix=0;ix<lambdaBar.size();++ix) {
        double y1 = rapidityAxis(lambdaBar[ix].momentum(),axis);
        for(unsigned int iy=ix+1;iy<lambdaBar.size();++iy) {
          double y2 = rapidityAxis(lambdaBar[iy].momentum(),axis);
          if(njet==2) _histdy->fill(fabs(y1-y2),-1.0);
          double ctheta = cosThetaStar(lambdaBar[ix],lambdaBar[iy],axis);
          _histdcosTheta->fill(ctheta,-1.0);
        }
      }
      // all lambda lambdabar
      for (unsigned int ix=0;ix<lambda.size();++ix) {
        double y1 = rapidityAxis(lambda[ix].momentum(),axis);
        for(unsigned int iy=0;iy<lambdaBar.size();++iy) {
          double y2 = rapidityAxis(lambdaBar[iy].momentum(),axis);
          if(njet==2) _histdy->fill(fabs(y1-y2));
          double ctheta = cosThetaStar(lambda[ix],lambdaBar[iy],axis);
          _histdcosTheta->fill(ctheta);
        }
      }
    }


    /// Normalise histograms etc., after the run
    void finalize() {
      scale(_histLLBar,1./sumOfWeights());
      scale(_histLL,1./sumOfWeights());
      scale(_histLLBarCorr,1./sumOfWeights());
      scale(_histLLBar_2jet,1./ *_weight2Jet);
      scale(_histLL_2jet,1./ *_weight2Jet);
      scale(_histLLBarCorr_2jet,1./ *_weight2Jet);
      scale(_histLLBar_3jet,1./ *_weight3Jet);
      scale(_histLL_3jet,1./ *_weight3Jet);
      scale(_histLLBarCorr_3jet,1./ *_weight3Jet);
      normalize(_histdcosTheta, 1.);
      normalize(_histdy       , 1.);
    }

    /// @}


    /// @name Histograms
    /// @{

    CounterPtr _weight2Jet;
    CounterPtr _weight3Jet;

    BinnedHistoPtr<string> _histLLBar, _histLL, _histLLBarCorr;
    BinnedHistoPtr<string> _histLLBar_2jet, _histLL_2jet, _histLLBarCorr_2jet;
    BinnedHistoPtr<string> _histLLBar_3jet, _histLL_3jet, _histLLBarCorr_3jet;

    Histo1DPtr _histdcosTheta;
    Histo1DPtr _histdy;

    const string Ecm = "91.2";
    /// @}


  };


  RIVET_DECLARE_PLUGIN(OPAL_2000_I474010);


}