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