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