1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/ChargedFinalState.hh"
  4#include "Rivet/Projections/NeutralFinalState.hh"
  5#include "Rivet/Projections/MergedFinalState.hh"
  6#include "Rivet/Projections/VetoedFinalState.hh"
  7#include "Rivet/Projections/FastJets.hh"
  8#include "fastjet/SISConePlugin.hh"
  9
 10namespace Rivet {
 11
 12
 13  /// @brief STAR underlying event
 14  ///
 15  /// @author Hendrik Hoeth
 16  class STAR_2009_UE_HELEN : public Analysis {
 17  public:
 18
 19    /// Constructor
 20    RIVET_DEFAULT_ANALYSIS_CTOR(STAR_2009_UE_HELEN);
 21
 22
 23    /// @name Analysis methods
 24    /// @{
 25
 26    void init() {
 27      // Charged final state, |eta|<1, pT>0.2GeV
 28      const Cut c = Cuts::abseta < 1.0 && Cuts::pT >= 0.2*GeV;
 29
 30      const ChargedFinalState cfs(c);
 31      declare(cfs, "CFS");
 32
 33      // Neutral final state, |eta|<1, ET>0.2GeV (needed for the jets)
 34      const NeutralFinalState nfs(c);
 35      declare(nfs, "NFS");
 36
 37      // STAR can't see neutrons and K^0_L
 38      VetoedFinalState vfs(nfs);
 39      vfs.vetoNeutrinos();
 40      vfs.addVetoPairId(PID::K0L);
 41      vfs.addVetoPairId(PID::NEUTRON);
 42      declare(vfs, "VFS");
 43
 44      // Jets are reconstructed from charged and neutral particles,
 45      // and the cuts are different (pT vs. ET), so we need to merge them.
 46      const MergedFinalState jfs(cfs, vfs);
 47      declare(jfs, "JFS");
 48
 49      // SISCone, R = 0.7, overlap_threshold = 0.75
 50      declare(FastJets(jfs, JetAlg::SISCONE, 0.7), "AllJets");
 51
 52      // Book histograms
 53      book(_hist_pmaxnchg, 1, 1, 1);
 54      book(_hist_pminnchg, 2, 1, 1);
 55      book(_hist_anchg,    3, 1, 1);
 56    }
 57
 58
 59    // Do the analysis
 60    void analyze(const Event& e) {
 61      const FinalState& cfs = apply<ChargedFinalState>(e, "CFS");
 62      if (cfs.particles().size() < 1) {
 63        MSG_DEBUG("Failed multiplicity cut");
 64        vetoEvent;
 65      }
 66
 67      const Jets& alljets = apply<FastJets>(e, "AllJets").jetsByPt();
 68      MSG_DEBUG("Total jet multiplicity = " << alljets.size());
 69
 70      // The jet acceptance region is |eta|<(1-R)=0.3  (with R = jet radius)
 71      // Jets also must have a neutral energy fraction of < 0.7
 72      Jets jets;
 73      for (const Jet& jet : alljets) {
 74        if (jet.neutralEnergy()/jet.totalEnergy() < 0.7 && jet.abseta() < 0.3) {
 75          jets.push_back(jet);
 76        }
 77      }
 78
 79      // This analysis requires a di-jet like event.
 80      // WARNING: There is more data in preparation, some of which
 81      //          does _not_ have this constraint!
 82      if (jets.size() != 2) {
 83        MSG_DEBUG("Failed jet multiplicity cut");
 84        vetoEvent;
 85      }
 86
 87      // The di-jet constraints in this analysis are:
 88      // - 2 and only 2 jets in the acceptance region
 89      // - delta(Phi) between the jets is > 150 degrees
 90      // - Pt_awayjet/Pt_towards_jet > 0.7
 91      if (deltaPhi(jets[0].phi(), jets[1].phi()) <= 5*PI/6 ||
 92          jets[1].pT()/jets[0].pT() <= 0.7)
 93      {
 94        MSG_DEBUG("Failed di-jet criteria");
 95        vetoEvent;
 96      }
 97
 98      // Now lets start ...
 99      const double jetphi = jets[0].phi();
100      const double jetpT  = jets[0].pT()/GeV;
101
102      size_t numTrans1(0), numTrans2(0), numAway(0);
103
104      // Calculate all the charged stuff
105      for (const Particle& p : cfs.particles()) {
106        const double dPhi = deltaPhi(p.phi(), jetphi);
107        const double pT = p.pT();
108        const double phi = p.phi();
109        double rotatedphi = phi - jetphi;
110        while (rotatedphi < 0) rotatedphi += 2*PI;
111
112        // @TODO: WARNING: The following lines are a hack to correct
113        //        for the STAR tracking efficiency. Once we have the
114        //        final numbers (corrected to hadron level), we need
115        //        to remove this!!!!
116        if (1.0*rand()/static_cast<double>(RAND_MAX) > 0.87834-exp(-1.48994-0.788432*pT)) {
117          continue;
118        }
119        // -------- end of efficiency hack -------
120
121        if (dPhi < PI/3.0) {
122          // toward
123        }
124        else if (dPhi < 2*PI/3.0) {
125          if (rotatedphi <= PI) {
126            ++numTrans1;
127          }
128          else {
129            ++numTrans2;
130          }
131        }
132        else {
133          ++numAway;
134        }
135      } // end charged particle loop
136
137      // Fill the histograms
138      _hist_pmaxnchg->fill(jetpT, double(numTrans1>numTrans2 ? numTrans1 : numTrans2)/(2*PI/3));
139      _hist_pminnchg->fill(jetpT, double(numTrans1<numTrans2 ? numTrans1 : numTrans2)/(2*PI/3));
140      _hist_anchg->fill(jetpT, (double)numAway/(PI*0.7*0.7)); // jet area = pi*R^2
141
142    }
143
144
145    void finalize() {
146      /// @todo Really nothing to do?
147    }
148
149    /// @}
150
151
152  private:
153
154    Profile1DPtr _hist_pmaxnchg;
155    Profile1DPtr _hist_pminnchg;
156    Profile1DPtr _hist_anchg;
157
158  };
159
160
161  RIVET_DECLARE_PLUGIN(STAR_2009_UE_HELEN);
162
163}