Rivet analyses

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

namespace Rivet {


  /// ATLAS jet substructure measurement
  class ATLAS_2012_I1094564 : public Analysis {
  public:

    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2012_I1094564);


    // Returns constituents to make it easier to do the filtering
    PseudoJets splitjet(fastjet::PseudoJet jet, double& last_R, const FastJets& fj, bool& unclustered) const {

      // Build a new cluster sequence just using the constituents of this jet.
      fastjet::ClusterSequence cs(fj.clusterSeq()->constituents(jet), fastjet::JetDefinition(fastjet::cambridge_algorithm, M_PI/2.));

      // Get the jet back again
      vector<fastjet::PseudoJet> remadeJets = cs.inclusive_jets(0.);

      if ( remadeJets.size() != 1 ) return remadeJets;

      fastjet::PseudoJet remadeJet = remadeJets[0];
      fastjet::PseudoJet parent1, parent2;
      unclustered = false;

      while ( cs.has_parents(remadeJet, parent1, parent2) ) {
        if (parent1.squared_distance(parent2) < 0.09) break;
        if (parent1.m2() < parent2.m2()) {
          fastjet::PseudoJet tmp;
          tmp = parent1; parent1 = parent2; parent2 = tmp;
        }

        double ktdist = parent1.kt_distance(parent2);
        double rtycut2 = 0.3*0.3;
        if (parent1.m() < 0.67*remadeJet.m() && ktdist > rtycut2*remadeJet.m2()) {
          unclustered = true;
          break;
        } else {
          remadeJet = parent1;
        }
      }

      last_R = 0.5 * sqrt(parent1.squared_distance(parent2));
      return cs.constituents(remadeJet);
    }


    fastjet::PseudoJet filterjet(PseudoJets jets, double& stingy_R, const double def_R) const {
      if (stingy_R == 0.0) stingy_R = def_R;
      stingy_R = def_R < stingy_R ? def_R : stingy_R;
      fastjet::JetDefinition stingy_jet_def(fastjet::cambridge_algorithm, stingy_R);
      fastjet::ClusterSequence scs(jets, stingy_jet_def);
      vector<fastjet::PseudoJet> stingy_jets = sorted_by_pt(scs.inclusive_jets(0));
      fastjet::PseudoJet reconst_jet(0, 0, 0, 0);
      for (size_t isj = 0; isj < std::min((size_t) 3, stingy_jets.size()); ++isj) {
        reconst_jet += stingy_jets[isj];
      }
      return reconst_jet;
    }


    // These are custom functions for n-subjettiness.
    PseudoJets jetGetAxes(int n_jets, const PseudoJets& inputJets, double subR) const {
      // Sanity check
      if (inputJets.size() < (size_t) n_jets) {
        MSG_ERROR("Not enough input particles.");
        return inputJets;
      }

      // Get subjets, return
      fastjet::ClusterSequence sub_clust_seq(inputJets, fastjet::JetDefinition(fastjet::kt_algorithm, subR, fastjet::E_scheme, fastjet::Best));
      return sub_clust_seq.exclusive_jets(n_jets);
    }


    double jetTauValue(double beta, double jet_rad, const PseudoJets& particles, const PseudoJets& axes, double Rcut) const {
      double tauNum = 0.0;
      double tauDen = 0.0;

      if (particles.size() == 0) return 0.0;

      for (size_t i = 0; i < particles.size(); i++) {
        // find minimum distance (set R large to begin)
        double minR = 10000.0;
        for (size_t j = 0; j < axes.size(); j++) {
          double tempR = sqrt(particles[i].squared_distance(axes[j]));
          if (tempR < minR) minR = tempR;
        }
        if (minR > Rcut) minR = Rcut;
        // calculate nominator and denominator
        tauNum += particles[i].perp() * pow(minR,beta);
        tauDen += particles[i].perp() * pow(jet_rad,beta);
      }

      // return N-subjettiness (or 0 if denominator is 0)
      return safediv(tauNum, tauDen, 0);
    }


    void jetUpdateAxes(double beta, const PseudoJets& particles, PseudoJets& axes) const {
      vector<int> belongsto;
      for (size_t i = 0; i < particles.size(); i++) {
        // find minimum distance axis
        int assign = 0;
        double minR = 10000.0;
        for (size_t j = 0; j < axes.size(); j++) {
          double tempR = sqrt(particles[i].squared_distance(axes[j]));
          if (tempR < minR) {
            minR = tempR;
            assign = j;
          }
        }
        belongsto.push_back(assign);
      }

      // iterative step
      double deltaR2, distphi;
      vector<double> ynom, phinom, den;

      ynom.resize(axes.size());
      phinom.resize(axes.size());
      den.resize(axes.size());

      for (size_t i = 0; i < particles.size(); i++) {
        distphi = particles[i].phi() - axes[belongsto[i]].phi();
        deltaR2 = particles[i].squared_distance(axes[belongsto[i]]);
        if (deltaR2 == 0.) continue;
        if (abs(distphi) <= M_PI) phinom.at(belongsto[i]) += particles[i].perp() * particles[i].phi() * pow(deltaR2, (beta-2)/2);
        else if ( distphi > M_PI) phinom.at(belongsto[i]) += particles[i].perp() * (-2 * M_PI + particles[i].phi()) * pow(deltaR2, (beta-2)/2);
        else if ( distphi < M_PI) phinom.at(belongsto[i]) += particles[i].perp() * (+2 * M_PI + particles[i].phi()) * pow(deltaR2, (beta-2)/2);
        ynom.at(belongsto[i]) += particles[i].perp() * particles[i].rap() * pow(deltaR2, (beta-2)/2);
        den.at(belongsto[i])  += particles[i].perp() * pow(deltaR2, (beta-2)/2);
      }

      // reset to new axes
      for (size_t j = 0; j < axes.size(); j++) {
        if (den[j] == 0.) axes.at(j) = axes[j];
        else {
          double phi_new = fmod( 2*M_PI + (phinom[j] / den[j]), 2*M_PI );
          double pt_new  = axes[j].perp();
          double y_new   = ynom[j] / den[j];
          double px = pt_new * cos(phi_new);
          double py = pt_new * sin(phi_new);
          double pz = pt_new * sinh(y_new);
          axes.at(j).reset(px, py, pz, axes[j].perp()/2);
        }
      }
    }


    void init() {

      /// Projections:
      FinalState fs((Cuts::etaIn(-4.5, 4.5) && Cuts::pT >=  0.*GeV));
      declare(fs, "FS");
      declare(FastJets(fs, JetAlg::ANTIKT, 1.0), "AKT");
      declare(FastJets(fs, JetAlg::CAM, 1.2)   , "CA" );

      /// Histograms:
      book(_h_camass,   {200., 300., 400., 500., 600.},
                        {"d01-x01-y01", "d02-x01-y01", "d03-x01-y01", "d04-x01-y01"});
      book(_h_filtmass, {200., 300., 400., 500., 600.},
                        {"d05-x01-y01", "d06-x01-y01", "d07-x01-y01", "d08-x01-y01"});
      book(_h_ktmass,   {200., 300., 400., 500., 600.},
                        {"d09-x01-y01", "d10-x01-y01", "d11-x01-y01", "d12-x01-y01"});
      book(_h_ktd12,    {200., 300., 400., 500., 600.},
                        {"d13-x01-y01", "d14-x01-y01", "d15-x01-y01", "d16-x01-y01"});
      book(_h_ktd23,    {200., 300., 400., 500., 600.},
                        {"d17-x01-y01", "d18-x01-y01", "d19-x01-y01", "d20-x01-y01"});
      book(_h_cat21,    {200., 300., 400., 500., 600.},
                        {"d21-x01-y01", "d22-x01-y01", "d23-x01-y01", "d24-x01-y01"});
      book(_h_cat32,    {200., 300., 400., 500., 600.},
                        {"d25-x01-y01", "d26-x01-y01", "d27-x01-y01", "d28-x01-y01"});
      book(_h_ktt21,    {200., 300., 400., 500., 600.},
                        {"d29-x01-y01", "d30-x01-y01", "d31-x01-y01", "d32-x01-y01"});
      book(_h_ktt32,    {200., 300., 400., 500., 600.},
                        {"d33-x01-y01", "d34-x01-y01", "d35-x01-y01", "d36-x01-y01"});

    }


    /// Perform the per-event analysis
    void analyze(const Event& event) {

      using namespace fastjet;

      // Get anti-kt jets with p_T > 200 GeV, check abs(y) < 2, and fill mass histograms
      const FastJets& ktfj = apply<FastJets>(event, "AKT");
      PseudoJets ktjets = ktfj.pseudojetsByPt(200*GeV);
      for (const PseudoJet & ajet : ktjets) {
        if (abs(ajet.rap()) < 2) {
          _h_ktmass->fill(ajet.perp(), ajet.m());
        }
      }

      // Same as above but C/A jets
      const FastJets& cafj = apply<FastJets>(event, "CA");
      PseudoJets cajets = cafj.pseudojetsByPt(200*GeV);
      for (const PseudoJet & ajet : cajets) {
        if (abs(ajet.rap()) < 2) {
          _h_camass->fill(ajet.perp(), ajet.m());
        }
      }

      // Split and filter.
      // Only do this to C/A jets in this analysis.
      for (const PseudoJet & pjet : cajets) {
        if ( pjet.perp() > 600 || abs(pjet.rap()) > 2) continue;
        double dR = 0;
        bool unclustered = false;
        PseudoJets split_jets = splitjet(pjet, dR, cafj, unclustered);
        if ( (dR < 0.15) || (unclustered == false) ) continue;
        PseudoJet filt_jet = filterjet(split_jets, dR, 0.3);
        _h_filtmass->fill(filt_jet.perp(), filt_jet.m());
      }

      // Use the two last stages of clustering to get sqrt(d_12) and sqrt(d_23).
      // Only use anti-kt jets in this analysis.
      for (const PseudoJet & pjet : ktjets) {
        if (pjet.perp() > 600 || abs(pjet.rap()) > 2) continue;
        ClusterSequence subjet_cseq(ktfj.clusterSeq()->constituents(pjet), JetDefinition(kt_algorithm, M_PI/2.));
        double d_12 = subjet_cseq.exclusive_dmerge(1) * M_PI*M_PI/4.;
        double d_23 = subjet_cseq.exclusive_dmerge(2) * M_PI*M_PI/4.;
        _h_ktd12->fill(pjet.perp(), sqrt(d_12));
        _h_ktd23->fill(pjet.perp(), sqrt(d_23));
      }

      // N-subjettiness, use beta = 1 (no rationale given).
      // Uses the functions defined above.
      // C/A jets first, anti-kt after.
      double beta = 1.;
      //Rcut is used for particles that are very far from the closest axis. At 10
      //is has no impact on the outcome of the calculation
      double Rcut = 10.;
      for (const PseudoJet & pjet : cajets) {
        if (pjet.perp() > 600*GeV || fabs(pjet.rap()) > 2) continue;

        const PseudoJets constituents = cafj.clusterSeq()->constituents(pjet);
        if (constituents.size() < 3) continue;

        const PseudoJets axis1 = jetGetAxes(1, constituents, M_PI/2.0);
        const PseudoJets axis2 = jetGetAxes(2, constituents, M_PI/2.0);
        const PseudoJets axis3 = jetGetAxes(3, constituents, M_PI/2.0);

        const double radius = 1.2;
        const double tau1 = jetTauValue(beta, radius, constituents, axis1, Rcut);
        const double tau2 = jetTauValue(beta, radius, constituents, axis2, Rcut);
        const double tau3 = jetTauValue(beta, radius, constituents, axis3, Rcut);

        if (tau1 == 0 || tau2 == 0) continue;
        _h_cat21->fill(pjet.perp(), tau2/tau1);
        _h_cat32->fill(pjet.perp(), tau3/tau2);
      }

      for (const PseudoJet & pjet : ktjets) {
        if (pjet.perp() > 600*GeV || fabs(pjet.rap()) > 2) continue;

        const PseudoJets constituents = ktfj.clusterSeq()->constituents(pjet);
        if (constituents.size() < 3) continue;

        const PseudoJets axis1 = jetGetAxes(1, constituents, M_PI/2.0);
        const PseudoJets axis2 = jetGetAxes(2, constituents, M_PI/2.0);
        const PseudoJets axis3 = jetGetAxes(3, constituents, M_PI/2.0);

        const double radius = 1.0;
        const double tau1 = jetTauValue(beta, radius, constituents, axis1, Rcut);
        const double tau2 = jetTauValue(beta, radius, constituents, axis2, Rcut);
        const double tau3 = jetTauValue(beta, radius, constituents, axis3, Rcut);
        if (tau1 == 0 || tau2 == 0) continue;

        _h_ktt21->fill(pjet.perp(), tau2/tau1);
        _h_ktt32->fill(pjet.perp(), tau3/tau2);
      }
    }


    /// Normalise histograms etc., after the run
    void finalize() {
      normalize(_h_camass);
      normalize(_h_filtmass);
      normalize(_h_ktmass);
      normalize(_h_ktd12);
      normalize(_h_ktd23);
      normalize(_h_cat21);
      normalize(_h_cat32);
      normalize(_h_ktt21);
      normalize(_h_ktt32);
    }

  private:

    Histo1DGroupPtr _h_camass;
    Histo1DGroupPtr _h_filtmass;
    Histo1DGroupPtr _h_ktmass;
    Histo1DGroupPtr _h_ktd12;
    Histo1DGroupPtr _h_ktd23;
    Histo1DGroupPtr _h_cat21;
    Histo1DGroupPtr _h_cat32;
    Histo1DGroupPtr _h_ktt21;
    Histo1DGroupPtr _h_ktt32;

  };


  RIVET_DECLARE_PLUGIN(ATLAS_2012_I1094564);

}