Rivet analyses

Transverse cone and `Swiss cheese’ underlying event studies

Experiment: CDF (Tevatron Run 1)

Inspire ID: 647490

Status: VALIDATED

Authors: - Andy Buckley

References: - Phys. Rev. D70, 072002 (2004) - arXiv: hep-ex/0404004

Beams: p- p+

Beam energies: (315.0, 315.0); (900.0, 900.0)GeV

Run details: - QCD events at $\sqrt{s} = 630$ & 1800 GeV. Several $\pT^\text{min}$ cutoffs are probably required to fill the profile histograms, e.g. 0 (min bias), 30, 90, 150 GeV at 1800 GeV, and 0 (min bias), 20, 90, 150 GeV at 630 GeV.

.This analysis studies the underlying event via transverse cones of R = 0.7 at 90 degrees in relative to the leading (highest E) jet, at $\sqrt{s} = 630$ and 1800 GeV. This is similar to the 2001 CDF UE analysis, except that cones, rather than the whole central range are used. The transverse cones are categorised as “TransMIN” and “TransMAX” on an event-by-event basis, to give greater sensitivity to the UE component. “Swiss Cheese” distributions, where cones around the leading n jets are excluded from the distributions, are also included for n = 2, 3. This analysis is useful for constraining the energy evolution of the underlying event, since it performs the same analyses at two distinct CoM energies. WARNING! The pT plots are normalised to raw number of events. The min bias data have not been reproduced by MC, and are not recommended for tuning.’

Source code:CDF_2004_I647490.cc

// -*- C++ -*-
// "Acosta" underlying event analysis at CDF, inc. "Swiss Cheese"

#include "Rivet/Analysis.hh"
#include "Rivet/Jet.hh"
#include "Rivet/Projections/ChargedFinalState.hh"
#include "Rivet/Projections/FastJets.hh"
#include "Rivet/Projections/TriggerCDFRun0Run1.hh"

namespace Rivet {


  /// @brief CDF calo jet underlying event analysis at 630 and 1800 GeV
  ///
  /// CDF measurement of underlying event using calorimeter jet scales and
  /// alignment, particle flow activity in transverse cones, and the Swiss
  /// Cheese analysis method, where cones are excluded around the 2 and 3
  /// hardest jets.
  ///
  /// @author Andy Buckley
  class CDF_2004_I647490 : public Analysis {
  public:

    RIVET_DEFAULT_ANALYSIS_CTOR(CDF_2004_I647490);


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

    void init() {
      // Set up projections
      declare(TriggerCDFRun0Run1(), "Trigger");
      const FinalState calofs(Cuts::abseta < 1.2);
      declare(calofs, "CaloFS");
      declare(FastJets(calofs, JetAlg::CDFJETCLU, 0.7), "Jets");
      const ChargedFinalState trackfs(Cuts::abseta < 1.2 && Cuts::pT >= 0.4*GeV);
      declare(trackfs, "TrackFS");
      // Restrict tracks to |eta| < 0.7 for the min bias part.
      const ChargedFinalState mbfs(Cuts::abseta < 0.7 && Cuts::pT >= 0.4*GeV);
      declare(mbfs, "MBFS");
      // Restrict tracks to |eta| < 1 for the Swiss-Cheese part.
      const ChargedFinalState cheesefs(Cuts::abseta < 1.0 && Cuts::pT >= 0.4*GeV);
      declare(cheesefs, "CheeseFS");
      declare(FastJets(cheesefs, JetAlg::CDFJETCLU, 0.7), "CheeseJets");

      // Book histograms
      for (double eVal : allowedEnergies()) {
        const string en = toString(round(eVal));
        if (isCompatibleWithSqrtS(eVal))  _sqs = en;

        if (en == "1800"s) {
          book(_p[en+"pt90MaxAvg"] ,1, 1, 1);
          book(_p[en+"pt90MinAvg"] ,1, 1, 2);
          book(_p[en+"pt90Max"] ,2, 1, 1);
          book(_p[en+"pt90Min"] ,2, 1, 2);
          book(_p[en+"pt90Diff"] ,2, 1, 3);
          book(_p[en+"num90Max"] ,4, 1, 1);
          book(_p[en+"num90Min"] ,4, 1, 2);
          book(_p[en+"pTSum_2Jet"] ,7, 1, 1);
          book(_p[en+"pTSum_3Jet"] ,7, 1, 2);

          book(_h[en+"pt90DbnEt40"] ,3, 1, 1);
          book(_h[en+"pt90DbnEt80"] ,3, 1, 2);
          book(_h[en+"pt90DbnEt120"] ,3, 1, 3);
          book(_h[en+"pt90DbnEt160"] ,3, 1, 4);
          book(_h[en+"pt90DbnEt200"] ,3, 1, 5);
          book(_h[en+"numTracksDbnMB"] ,5, 1, 1);
          book(_h[en+"ptDbnMB"] ,6, 1, 1);
        }
        else {
          book(_p[en+"pt90Max"] ,8, 1, 1);
          book(_p[en+"pt90Min"] ,8, 1, 2);
          book(_p[en+"pt90Diff"] ,8, 1, 3);
          book(_p[en+"pTSum_2Jet"] ,9, 1, 1);
          book(_p[en+"pTSum_3Jet"] ,9, 1, 2);
          book(_h[en+"numTracksDbnMB"] ,10, 1, 1);
          book(_h[en+"ptDbnMB"] ,11, 1, 1);
        }
      }
      raiseBeamErrorIf(_sqs.empty());
    }


    /// Do the analysis
    void analyze(const Event& event) {
      // Trigger
      const bool trigger = apply<TriggerCDFRun0Run1>(event, "Trigger").minBiasDecision();
      if (!trigger) vetoEvent;

      {
        MSG_DEBUG("Running max/min analysis");
        Jets jets = apply<JetFinder>(event, "Jets").jets(cmpMomByE);
        if (!jets.empty()) {
          // Leading jet must be in central |eta| < 0.5 region
          const Jet leadingjet = jets.front();
          const double etaLead = leadingjet.eta();
          // Get Et of the leading jet: used to bin histograms
          const double ETlead = leadingjet.Et();
          MSG_DEBUG("Leading Et = " << ETlead/GeV << " GeV");
          if (fabs(etaLead) > 0.5 && ETlead < 15*GeV) {
            MSG_DEBUG("Leading jet eta = " << etaLead
                     << " not in |eta| < 0.5 & pT > 15 GeV");
          } else {
            // Multiplicity & pT distributions for sqrt(s) = 630 GeV, 1800 GeV
            const Particles tracks = apply<FinalState>(event, "TrackFS").particles();
            const ConesInfo cones = _calcTransCones(leadingjet.momentum(), tracks);
            if (_sqs == "630"s) {
              _p[_sqs+"pt90Max"]->fill(ETlead/GeV, cones.ptMax/GeV);
              _p[_sqs+"pt90Min"]->fill(ETlead/GeV, cones.ptMin/GeV);
              _p[_sqs+"pt90Diff"]->fill(ETlead/GeV, cones.ptDiff/GeV);
            }
            else {
              _p[_sqs+"num90Max"]->fill(ETlead/GeV, cones.numMax);
              _p[_sqs+"num90Min"]->fill(ETlead/GeV, cones.numMin);
              _p[_sqs+"pt90Max"]->fill(ETlead/GeV, cones.ptMax/GeV);
              _p[_sqs+"pt90Min"]->fill(ETlead/GeV, cones.ptMin/GeV);
              _p[_sqs+"pt90Diff"]->fill(ETlead/GeV, cones.ptDiff/GeV);
              _p[_sqs+"pt90MaxAvg"]->fill(ETlead/GeV, cones.ptMax/GeV); // /numMax
              _p[_sqs+"pt90MinAvg"]->fill(ETlead/GeV, cones.ptMin/GeV); // /numMin
              //
              const double ptTransTotal = cones.ptMax + cones.ptMin;
              if (inRange(ETlead/GeV, 40., 80.)) {
                _h[_sqs+"pt90DbnEt40"]->fill(ptTransTotal/GeV);
              } else if (inRange(ETlead/GeV, 80., 120.)) {
                _h[_sqs+"pt90DbnEt80"]->fill(ptTransTotal/GeV);
              } else if (inRange(ETlead/GeV, 120., 160.)) {
                _h[_sqs+"pt90DbnEt120"]->fill(ptTransTotal/GeV);
              } else if (inRange(ETlead/GeV, 160., 200.)) {
                _h[_sqs+"pt90DbnEt160"]->fill(ptTransTotal/GeV);
              } else if (inRange(ETlead/GeV, 200., 270.)) {
                _h[_sqs+"pt90DbnEt200"]->fill(ptTransTotal/GeV);
              }
            }

          }
        }
      }


      // Fill min bias total track multiplicity histos
      {
        MSG_DEBUG("Running min bias multiplicity analysis");
        const Particles mbtracks = apply<FinalState>(event, "MBFS").particles();
        _h[_sqs+"numTracksDbnMB"]->fill(mbtracks.size());
        // Run over all charged tracks
        for (const Particle& t : mbtracks) {
          FourMomentum trackMom = t.momentum();
          const double pt = trackMom.pT();
          // Plot total pT distribution for min bias
          _h[_sqs+"ptDbnMB"]->fill(pt/GeV);
        }
      }



      // Construct "Swiss Cheese" pT distributions, with pT contributions from
      // tracks within R = 0.7 of the 1st, 2nd (and 3rd) jets being ignored. A
      // different set of charged tracks, with |eta| < 1.0, is used here, and all
      // the removed jets must have Et > 5 GeV.
      {
        MSG_DEBUG("Running Swiss Cheese analysis");
        const Particles cheesetracks = apply<FinalState>(event, "CheeseFS").particles();
        Jets cheesejets = apply<JetFinder>(event, "Jets").jets(cmpMomByE);
        if (cheesejets.empty()) {
          MSG_DEBUG("No 'cheese' jets found in event");
          return;
        }
        if (cheesejets.size() > 1 &&
            fabs(cheesejets[0].eta()) <= 0.5 &&
            cheesejets[0].Et()/GeV > 5.0 &&
            cheesejets[1].Et()/GeV > 5.0) {

          const double cheeseETlead = cheesejets[0].Et();

          const double eta1 = cheesejets[0].eta();
          const double phi1 = cheesejets[0].phi();
          const double eta2 = cheesejets[1].eta();
          const double phi2 = cheesejets[1].phi();

          double ptSumSub2(0), ptSumSub3(0);
          for (const Particle& t : cheesetracks) {
            FourMomentum trackMom = t.mom();
            const double pt = trackMom.pT();

            // Subtracting 2 leading jets
            const double deltaR1 = deltaR(trackMom, eta1, phi1);
            const double deltaR2 = deltaR(trackMom, eta2, phi2);
            MSG_TRACE("Track vs jet(1): "
                     << "|(" << trackMom.eta() << ", " << trackMom.phi() << ") - "
                     << "|(" << eta1 << ", " << phi1 << ")| = " << deltaR1);
            MSG_TRACE("Track vs jet(2): "
                     << "|(" << trackMom.eta() << ", " << trackMom.phi() << ") - "
                     << "|(" << eta2 << ", " << phi2 << ")| = " << deltaR2);
            if (deltaR1 > 0.7 && deltaR2 > 0.7) {
              ptSumSub2 += pt;

              // Subtracting 3rd leading jet
              if (cheesejets.size() > 2 &&
                  cheesejets[2].Et()/GeV > 5.0) {
                const double eta3 = cheesejets[2].eta();
                const double phi3 = cheesejets[2].phi();
                const double deltaR3 = deltaR(trackMom, eta3, phi3);
                MSG_TRACE("Track vs jet(3): "
                         << "|(" << trackMom.eta() << ", " << trackMom.phi() << ") - "
                         << "|(" << eta3 << ", " << phi3 << ")| = " << deltaR3);
                if (deltaR3 > 0.7) {
                  ptSumSub3 += pt;
                }
              }
            }
          }

          // Swiss Cheese sub 2,3 jets distributions for sqrt(s) = 630 GeV, 1800 GeV
          if (_sqs == "630"s) {
            if (!isZero(ptSumSub2)) _p[_sqs+"pTSum_2Jet"]->fill(cheeseETlead/GeV, ptSumSub2/GeV);
            if (!isZero(ptSumSub3)) _p[_sqs+"pTSum_3Jet"]->fill(cheeseETlead/GeV, ptSumSub3/GeV);
          }
          else {
            if (!isZero(ptSumSub2)) _p[_sqs+"pTSum_2Jet"]->fill(cheeseETlead/GeV, ptSumSub2/GeV);
            if (!isZero(ptSumSub3)) _p[_sqs+"pTSum_3Jet"]->fill(cheeseETlead/GeV, ptSumSub3/GeV);
          }

        }
      }

    }


    void finalize() {
      /// @todo Take these normalisations from the data histo (it can't come from just the MC)

      for (double eVal : allowedEnergies()) {
        const string en = toString(round(eVal));
        if (en == "1800"s) {
          // Normalize to actual number of entries in pT dbn histos...
          normalize(_h[en+"pt90DbnEt40"],  1656.75); // norm OK
          normalize(_h[en+"pt90DbnEt80"],  4657.5); // norm OK
          normalize(_h[en+"pt90DbnEt120"], 5395.5); // norm OK
          normalize(_h[en+"pt90DbnEt160"], 7248.75); // norm OK
          normalize(_h[en+"pt90DbnEt200"], 2442.0); // norm OK
          normalize(_h[en+"numTracksDbnMB"], 309718.25); // norm OK
          normalize(_h[en+"ptDbnMB"], 33600.0); // norm OK
        }
        else {
          normalize(_h[en+"numTracksDbnMB"], 1101024.0); // norm OK
          normalize(_h[en+"ptDbnMB"], 105088.0); // norm OK
        }
      }
    }

    /// @}


  private:


    /// @name Cone machinery
    /// @{

    /// @cond CONEUE_DETAIL

    struct ConesInfo {
      ConesInfo() : numMax(0), numMin(0), ptMax(0), ptMin(0), ptDiff(0) {}
      unsigned int numMax, numMin;
      double ptMax, ptMin, ptDiff;
    };

    /// @endcond


    ConesInfo _calcTransCones(const double etaLead, const double phiLead,
                              const Particles& tracks) {
      const double phiTransPlus = mapAngle0To2Pi(phiLead + PI/2.0);
      const double phiTransMinus = mapAngle0To2Pi(phiLead - PI/2.0);
      MSG_DEBUG("phi_lead = " << phiLead
               << " -> trans = (" << phiTransPlus
               << ", " << phiTransMinus << ")");

      unsigned int numPlus(0), numMinus(0);
      double ptPlus(0), ptMinus(0);
      // Run over all charged tracks
      for (const Particle& t : tracks) {
        FourMomentum trackMom = t.momentum();
        const double pt = trackMom.pT();

        // Find if track mom is in either transverse cone
        if (deltaR(trackMom, etaLead, phiTransPlus) < 0.7) {
          ptPlus += pt;
          numPlus += 1;
        } else if (deltaR(trackMom, etaLead, phiTransMinus) < 0.7) {
          ptMinus += pt;
          numMinus += 1;
        }
      }

      ConesInfo rtn;
      // Assign N_{min,max} from N_{plus,minus}
      rtn.numMax = (ptPlus >= ptMinus) ? numPlus : numMinus;
      rtn.numMin = (ptPlus >= ptMinus) ? numMinus : numPlus;
      // Assign pT_{min,max} from pT_{plus,minus}
      rtn.ptMax = (ptPlus >= ptMinus) ? ptPlus : ptMinus;
      rtn.ptMin = (ptPlus >= ptMinus) ? ptMinus : ptPlus;
      rtn.ptDiff = fabs(rtn.ptMax - rtn.ptMin);

      MSG_DEBUG("Min cone has " << rtn.numMin << " tracks -> "
               << "pT_min = " << rtn.ptMin/GeV << " GeV");
      MSG_DEBUG("Max cone has " << rtn.numMax << " tracks -> "
               << "pT_max = " << rtn.ptMax/GeV << " GeV");

      return rtn;
    }


    ConesInfo _calcTransCones(const FourMomentum& leadvec,
                              const Particles& tracks) {
      const double etaLead = leadvec.eta();
      const double phiLead = leadvec.phi();
      return _calcTransCones(etaLead, phiLead, tracks);
    }

    /// @}


    /// @name Histogram collections
    /// @{

    /// Profile histograms
    map<string,Profile1DPtr> _p;

    /// Histogram
    map<string,Histo1DPtr> _h;

    string _sqs = "";

    /// @}

  };



  RIVET_DECLARE_ALIASED_PLUGIN(CDF_2004_I647490, CDF_2004_S5839831);

}

Aliases: - CDF_2004_S5839831