Rivet analyses

// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/Beam.hh"
#include "Rivet/Projections/ChargedFinalState.hh"
#include "Rivet/Projections/FastJets.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/Hemispheres.hh"
#include "Rivet/Projections/ParisiTensor.hh"
#include "Rivet/Projections/Sphericity.hh"
#include "Rivet/Projections/Thrust.hh"
#include "Rivet/Projections/UnstableParticles.hh"

namespace Rivet {


  /// @brief ALEPH QCD study with event shapes and identified particles
  ///
  /// @author Holger Schulz
  class ALEPH_1996_I428072 : public Analysis {
  public:

    RIVET_DEFAULT_ANALYSIS_CTOR(ALEPH_1996_I428072);


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

    void init() {
      // Set up projections
      declare(Beam(), "Beams");
      const ChargedFinalState cfs;
      declare(cfs, "FS");
      declare(UnstableParticles(), "UFS");
      declare(FastJets(cfs, JetAlg::DURHAM, 0.7), "DurhamJets");
      declare(Sphericity(cfs), "Sphericity");
      declare(ParisiTensor(cfs), "Parisi");
      const Thrust thrust(cfs);
      declare(thrust, "Thrust");
      declare(Hemispheres(thrust), "Hemispheres");

      // Book histograms
      book(_histSphericity, 1, 1, 1);
      book(_histAplanarity, 2, 1, 1);

      book(_hist1MinusT, 3, 1, 1);
      book(_histTMinor, 4, 1, 1);

      book(_histY3, 5, 1, 1);
      book(_histHeavyJetMass, 6, 1, 1);
      book(_histCParam, 7, 1, 1);
      book(_histOblateness, 8, 1, 1);

      book(_histScaledMom, 9, 1, 1);
      book(_histRapidityT, 10, 1, 1);

      book(_histPtSIn, 11, 1, 1);
      book(_histPtSOut, 12, 1, 1);

      book(_histLogScaledMom, 17, 1, 1);

      book(_histChMult, 18, 1, 1);

      book(_histMeanChMult, 19, 1, 1);
      book(_histMeanChMultRapt05, 20, 1, 1);
      book(_histMeanChMultRapt10, 21, 1, 1);
      book(_histMeanChMultRapt15, 22, 1, 1);
      book(_histMeanChMultRapt20, 23, 1, 1);


      // Particle spectra
      book(_histMultiPiPlus, 25, 1, 1);
      book(_histMultiKPlus, 26, 1, 1);
      book(_histMultiP, 27, 1, 1);
      book(_histMultiPhoton, 28, 1, 1);
      book(_histMultiPi0, 29, 1, 1);
      book(_histMultiEta, 30, 1, 1);
      book(_histMultiEtaPrime, 31, 1, 1);
      book(_histMultiK0, 32, 1, 1);
      book(_histMultiLambda0, 33, 1, 1);
      book(_histMultiXiMinus, 34, 1, 1);
      book(_histMultiSigma1385Plus, 35, 1, 1);
      book(_histMultiXi1530_0, 36, 1, 1);
      book(_histMultiRho, 37, 1, 1);
      book(_histMultiOmega782, 38, 1, 1);
      book(_histMultiKStar892_0, 39, 1, 1);
      book(_histMultiPhi, 40, 1, 1);

      book(_histMultiKStar892Plus, 43, 1, 1);

      // Mean multiplicities
      book(_histMeanMultiPi0, 44, 1, 2);
      book(_histMeanMultiEta, 44, 1, 3);
      book(_histMeanMultiEtaPrime, 44, 1, 4);
      book(_histMeanMultiK0, 44, 1, 5);
      book(_histMeanMultiRho, 44, 1, 6);
      book(_histMeanMultiOmega782, 44, 1, 7);
      book(_histMeanMultiPhi, 44, 1, 8);
      book(_histMeanMultiKStar892Plus, 44, 1, 9);
      book(_histMeanMultiKStar892_0, 44, 1, 10);
      book(_histMeanMultiLambda0, 44, 1, 11);
      book(_histMeanMultiSigma0, 44, 1, 12);
      book(_histMeanMultiXiMinus, 44, 1, 13);
      book(_histMeanMultiSigma1385Plus, 44, 1, 14);
      book(_histMeanMultiXi1530_0, 44, 1, 15);
      book(_histMeanMultiOmegaOmegaBar, 44, 1, 16);
      book(_weightedTotalPartNum, "/TMP/TotalPartNum");

      book(_weightedTotalPartNum, "/TMP/weightedTotalPartNum");
    }


    void analyze(const Event& e) {
      // First, veto on leptonic events by requiring at least 4 charged FS particles
      const FinalState& fs = apply<FinalState>(e, "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");

      _weightedTotalPartNum->fill(numParticles);

      // Get beams and average beam momentum
      const ParticlePair& beams = apply<Beam>(e, "Beams").beams();
      const double meanBeamMom = (beams.first.p3().mod() + beams.second.p3().mod()) / 2.0;
      MSG_DEBUG("Avg beam momentum = " << meanBeamMom);

      // Thrusts
      MSG_DEBUG("Calculating thrust");
      const Thrust& thrust = apply<Thrust>(e, "Thrust");
      _hist1MinusT->fill(1 - thrust.thrust());
      _histTMinor->fill(thrust.thrustMinor());
      _histOblateness->fill(thrust.oblateness());

      // Jets
      MSG_DEBUG("Calculating differential jet rate plots:");
      const FastJets& durjet = apply<FastJets>(e, "DurhamJets");
      if (durjet.clusterSeq()) {
        double y3 = durjet.clusterSeq()->exclusive_ymerge_max(2);
        if (y3 > 0.0) _histY3->fill(-1. * std::log(y3));
      }

      // Sphericities
      MSG_DEBUG("Calculating sphericity");
      const Sphericity& sphericity = apply<Sphericity>(e, "Sphericity");
      _histSphericity->fill(sphericity.sphericity());
      _histAplanarity->fill(sphericity.aplanarity());

      // C param
      MSG_DEBUG("Calculating Parisi params");
      const ParisiTensor& parisi = apply<ParisiTensor>(e, "Parisi");
      _histCParam->fill(parisi.C());

      // Hemispheres
      MSG_DEBUG("Calculating hemisphere variables");
      const Hemispheres& hemi = apply<Hemispheres>(e, "Hemispheres");
      _histHeavyJetMass->fill(hemi.scaledM2high());

      // Iterate over all the charged final state particles.
      double rapt05 = 0.;
      double rapt10 = 0.;
      double rapt15 = 0.;
      double rapt20 = 0.;
      MSG_DEBUG("About to iterate over charged FS particles");
      for (const Particle& p : fs.particles()) {
        // Get momentum and energy of each particle.
        const Vector3 mom3 = p.p3();
        const double energy = p.E();

        // Scaled momenta.
        const double mom = mom3.mod();
        const double scaledMom = mom / meanBeamMom;
        const double logInvScaledMom = -std::log(scaledMom);
        _histLogScaledMom->fill(logInvScaledMom);
        _histScaledMom->fill(scaledMom);

        // Get momenta components w.r.t. thrust and sphericity.
        const double momT = dot(thrust.thrustAxis(), mom3);
        const double pTinS = dot(mom3, sphericity.sphericityMajorAxis());
        const double pToutS = dot(mom3, sphericity.sphericityMinorAxis());
        _histPtSIn->fill(fabs(pTinS / GeV));
        _histPtSOut->fill(fabs(pToutS / GeV));

        // Calculate rapidities w.r.t. thrust.
        const double rapidityT = 0.5 * std::log((energy + momT) / (energy - momT));
        _histRapidityT->fill(fabs(rapidityT));
        if (std::fabs(rapidityT) <= 0.5) {
          rapt05 += 1.0;
        }
        if (std::fabs(rapidityT) <= 1.0) {
          rapt10 += 1.0;
        }
        if (std::fabs(rapidityT) <= 1.5) {
          rapt15 += 1.0;
        }
        if (std::fabs(rapidityT) <= 2.0) {
          rapt20 += 1.0;
        }
      }

      _histChMult->fill(numParticles);

      _histMeanChMultRapt05->fill(Ecms, rapt05);
      _histMeanChMultRapt10->fill(Ecms, rapt10);
      _histMeanChMultRapt15->fill(Ecms, rapt15);
      _histMeanChMultRapt20->fill(Ecms, rapt20);
      _histMeanChMult->fill(Ecms, numParticles);


      //// Final state of unstable particles to get particle spectra
      const UnstableParticles& ufs = apply<UnstableParticles>(e, "UFS");
      for (Particles::const_iterator p = ufs.particles().begin(); p != ufs.particles().end(); ++p) {
        const Vector3 mom3 = p->momentum().p3();
        int id = abs(p->pid());
        const double mom = mom3.mod();
        const double energy = p->momentum().E();
        const double scaledMom = mom / meanBeamMom;
        const double scaledEnergy = energy / meanBeamMom; // meanBeamMom is approximately beam energy
        switch (id) {
          case 22: _histMultiPhoton->fill(-1. * std::log(scaledMom)); break;
          case -321:
          case 321: _histMultiKPlus->fill(scaledMom); break;
          case 211:
          case -211: _histMultiPiPlus->fill(scaledMom); break;
          case 2212:
          case -2212: _histMultiP->fill(scaledMom); break;
          case 111:
            _histMultiPi0->fill(scaledMom);
            _histMeanMultiPi0->fill(Ecms);
            break;
          case 221:
            if (scaledMom >= 0.1) {
              _histMultiEta->fill(scaledEnergy);
              _histMeanMultiEta->fill(Ecms);
            }
            break;
          case 331:
            if (scaledMom >= 0.1) {
              _histMultiEtaPrime->fill(scaledEnergy);
              _histMeanMultiEtaPrime->fill(Ecms);
            }
            break;
          case 130: //klong
          case 310: //kshort
            _histMultiK0->fill(scaledMom);
            _histMeanMultiK0->fill(Ecms);
            break;
          case 113:
            _histMultiRho->fill(scaledMom);
            _histMeanMultiRho->fill(Ecms);
            break;
          case 223:
            _histMultiOmega782->fill(scaledMom);
            _histMeanMultiOmega782->fill(Ecms);
            break;
          case 333:
            _histMultiPhi->fill(scaledMom);
            _histMeanMultiPhi->fill(Ecms);
            break;
          case 313:
          case -313:
            _histMultiKStar892_0->fill(scaledMom);
            _histMeanMultiKStar892_0->fill(Ecms);
            break;
          case 323:
          case -323:
            _histMultiKStar892Plus->fill(scaledEnergy);
            _histMeanMultiKStar892Plus->fill(Ecms);
            break;
          case 3122:
          case -3122:
            _histMultiLambda0->fill(scaledMom);
            _histMeanMultiLambda0->fill(Ecms);
            break;
          case 3212:
          case -3212: _histMeanMultiSigma0->fill(Ecms); break;
          case 3312:
          case -3312:
            _histMultiXiMinus->fill(scaledEnergy);
            _histMeanMultiXiMinus->fill(Ecms);
            break;
          case 3114:
          case -3114:
          case 3224:
          case -3224:
            _histMultiSigma1385Plus->fill(scaledEnergy);
            _histMeanMultiSigma1385Plus->fill(Ecms);
            break;
          case 3324:
          case -3324:
            _histMultiXi1530_0->fill(scaledEnergy);
            _histMeanMultiXi1530_0->fill(Ecms);
            break;
          case 3334: _histMeanMultiOmegaOmegaBar->fill(Ecms); break;
        }
      }
    }


    /// Finalize
    void finalize() {
      // Normalize inclusive single particle distributions to the average number
      // of charged particles per event.
      const double avgNumParts = _weightedTotalPartNum->sumW() / sumOfWeights();

      normalize(_histPtSIn, avgNumParts);
      normalize(_histPtSOut, avgNumParts);

      normalize(_histRapidityT, avgNumParts);
      normalize(_histY3);

      normalize(_histLogScaledMom, avgNumParts);
      normalize(_histScaledMom, avgNumParts);

      // particle spectra
      scale(_histMultiPiPlus, 1. / sumOfWeights());
      scale(_histMultiKPlus, 1. / sumOfWeights());
      scale(_histMultiP, 1. / sumOfWeights());
      scale(_histMultiPhoton, 1. / sumOfWeights());
      scale(_histMultiPi0, 1. / sumOfWeights());
      scale(_histMultiEta, 1. / sumOfWeights());
      scale(_histMultiEtaPrime, 1. / sumOfWeights());
      scale(_histMultiK0, 1. / sumOfWeights());
      scale(_histMultiLambda0, 1. / sumOfWeights());
      scale(_histMultiXiMinus, 1. / sumOfWeights());
      scale(_histMultiSigma1385Plus, 1. / sumOfWeights());
      scale(_histMultiXi1530_0, 1. / sumOfWeights());
      scale(_histMultiRho, 1. / sumOfWeights());
      scale(_histMultiOmega782, 1. / sumOfWeights());
      scale(_histMultiKStar892_0, 1. / sumOfWeights());
      scale(_histMultiPhi, 1. / sumOfWeights());

      scale(_histMultiKStar892Plus, 1. / sumOfWeights());

      // event shape
      normalize(_hist1MinusT);
      normalize(_histTMinor);
      normalize(_histOblateness);

      normalize(_histSphericity);
      normalize(_histAplanarity);
      normalize(_histHeavyJetMass);
      normalize(_histCParam);


      // mean multiplicities
      scale(_histChMult, 1.0 / sumOfWeights()); // taking into account the binwidth of 2
      scale(_histMeanChMult, 1.0 / sumOfWeights());
      scale(_histMeanChMultRapt05, 1.0 / sumOfWeights());
      scale(_histMeanChMultRapt10, 1.0 / sumOfWeights());
      scale(_histMeanChMultRapt15, 1.0 / sumOfWeights());
      scale(_histMeanChMultRapt20, 1.0 / sumOfWeights());


      scale(_histMeanMultiPi0, 1.0 / sumOfWeights());
      scale(_histMeanMultiEta, 1.0 / sumOfWeights());
      scale(_histMeanMultiEtaPrime, 1.0 / sumOfWeights());
      scale(_histMeanMultiK0, 1.0 / sumOfWeights());
      scale(_histMeanMultiRho, 1.0 / sumOfWeights());
      scale(_histMeanMultiOmega782, 1.0 / sumOfWeights());
      scale(_histMeanMultiPhi, 1.0 / sumOfWeights());
      scale(_histMeanMultiKStar892Plus, 1.0 / sumOfWeights());
      scale(_histMeanMultiKStar892_0, 1.0 / sumOfWeights());
      scale(_histMeanMultiLambda0, 1.0 / sumOfWeights());
      scale(_histMeanMultiSigma0, 1.0 / sumOfWeights());
      scale(_histMeanMultiXiMinus, 1.0 / sumOfWeights());
      scale(_histMeanMultiSigma1385Plus, 1.0 / sumOfWeights());
      scale(_histMeanMultiXi1530_0, 1.0 / sumOfWeights());
      scale(_histMeanMultiOmegaOmegaBar, 1.0 / sumOfWeights());
    }

    /// @}


  private:

    /// Store the weighted sums of numbers of charged / charged+neutral
    /// particles - used to calculate average number of particles for the
    /// inclusive single particle distributions' normalisations.
    CounterPtr _weightedTotalPartNum;

    const string Ecms = "91.2";

    /// @name Histograms
    /// @{
    Histo1DPtr _histSphericity;
    Histo1DPtr _histAplanarity;

    Histo1DPtr _hist1MinusT;
    Histo1DPtr _histTMinor;

    Histo1DPtr _histY3;
    Histo1DPtr _histHeavyJetMass;
    Histo1DPtr _histCParam;
    Histo1DPtr _histOblateness;

    Histo1DPtr _histScaledMom;
    Histo1DPtr _histRapidityT;

    Histo1DPtr _histPtSIn;
    Histo1DPtr _histPtSOut;

    Histo1DPtr _histJetRate2Durham;
    Histo1DPtr _histJetRate3Durham;
    Histo1DPtr _histJetRate4Durham;
    Histo1DPtr _histJetRate5Durham;

    Histo1DPtr _histLogScaledMom;

    BinnedHistoPtr<int> _histChMult;

    Histo1DPtr _histMultiPiPlus;
    Histo1DPtr _histMultiKPlus;
    Histo1DPtr _histMultiP;
    Histo1DPtr _histMultiPhoton;
    Histo1DPtr _histMultiPi0;
    Histo1DPtr _histMultiEta;
    Histo1DPtr _histMultiEtaPrime;
    Histo1DPtr _histMultiK0;
    Histo1DPtr _histMultiLambda0;
    Histo1DPtr _histMultiXiMinus;
    Histo1DPtr _histMultiSigma1385Plus;
    Histo1DPtr _histMultiXi1530_0;
    Histo1DPtr _histMultiRho;
    Histo1DPtr _histMultiOmega782;
    Histo1DPtr _histMultiKStar892_0;
    Histo1DPtr _histMultiPhi;
    Histo1DPtr _histMultiKStar892Plus;

    // mean multiplicities
    BinnedHistoPtr<string> _histMeanChMult;
    BinnedHistoPtr<string> _histMeanChMultRapt05;
    BinnedHistoPtr<string> _histMeanChMultRapt10;
    BinnedHistoPtr<string> _histMeanChMultRapt15;
    BinnedHistoPtr<string> _histMeanChMultRapt20;

    BinnedHistoPtr<string> _histMeanMultiPi0;
    BinnedHistoPtr<string> _histMeanMultiEta;
    BinnedHistoPtr<string> _histMeanMultiEtaPrime;
    BinnedHistoPtr<string> _histMeanMultiK0;
    BinnedHistoPtr<string> _histMeanMultiRho;
    BinnedHistoPtr<string> _histMeanMultiOmega782;
    BinnedHistoPtr<string> _histMeanMultiPhi;
    BinnedHistoPtr<string> _histMeanMultiKStar892Plus;
    BinnedHistoPtr<string> _histMeanMultiKStar892_0;
    BinnedHistoPtr<string> _histMeanMultiLambda0;
    BinnedHistoPtr<string> _histMeanMultiSigma0;
    BinnedHistoPtr<string> _histMeanMultiXiMinus;
    BinnedHistoPtr<string> _histMeanMultiSigma1385Plus;
    BinnedHistoPtr<string> _histMeanMultiXi1530_0;
    BinnedHistoPtr<string> _histMeanMultiOmegaOmegaBar;
    /// @}
  };


  RIVET_DECLARE_ALIASED_PLUGIN(ALEPH_1996_I428072, ALEPH_1996_S3486095);

}