1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/ChargedFinalState.hh"
  5#include "Rivet/Projections/TriggerCDFRun2.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// @brief CDF Run II min-bias cross-section
 11  ///
 12  /// @author Hendrik Hoeth
 13  ///
 14  /// Measurement of \f$ \langle p_T \rangle \f$ vs. \f$ n_\text{ch} \f$,
 15  /// the track \f$ p_T \f$ distribution, and the \f$ \sum E_T \f$ distribution.
 16  /// Particles are selected within |eta|<1 and with pT>0.4 GeV.
 17  /// There is no pT cut for the \f$ \sum E_T \f$ measurement.
 18  ///
 19  /// @par Run conditions
 20  ///
 21  /// @arg \f$ \sqrt{s} = \f$ 1960 GeV
 22  /// @arg Run with generic QCD events.
 23  /// @arg Set particles with c*tau > 10 mm stable
 24  class CDF_2009_I817466 : public Analysis {
 25  public:
 26
 27    RIVET_DEFAULT_ANALYSIS_CTOR(CDF_2009_I817466);
 28
 29
 30    /// @name Analysis methods
 31    /// @{
 32
 33    /// Book histograms and projections
 34    void init() {
 35      declare(TriggerCDFRun2(), "Trigger");
 36      declare(FinalState((Cuts::etaIn(-1.0, 1.0))), "EtFS");
 37      declare(ChargedFinalState((Cuts::etaIn(-1.0, 1.0) && Cuts::pT >=  0.4*GeV)), "CFS");
 38
 39      book(_hist_pt ,1, 1, 1);
 40      book(_hist_pt_vs_multiplicity ,2, 1, 1);
 41      book(_hist_sumEt ,3, 1, 1);
 42
 43      book(_sumWeightSelected,"_sumWeightSelected");
 44    }
 45
 46
 47    /// Do the analysis
 48    void analyze(const Event& evt) {
 49      // MinBias Trigger
 50      const bool trigger = apply<TriggerCDFRun2>(evt, "Trigger").minBiasDecision();
 51      if (!trigger) vetoEvent;
 52
 53      /// @todo The pT and sum(ET) distributions look slightly different from
 54      ///       Niccolo's Monte Carlo plots. Still waiting for his answer.
 55
 56      const ChargedFinalState& trackfs = apply<ChargedFinalState>(evt, "CFS");
 57      const size_t numParticles = trackfs.size();
 58      for (const Particle& p : trackfs.particles()) {
 59        const double pT = p.pT() / GeV;
 60        _hist_pt_vs_multiplicity->fill(numParticles, pT);
 61
 62        // The weight for entries in the pT distribution should be weight/(pT*dPhi*dy).
 63        //
 64        // - dPhi = 2*PI
 65        //
 66        // - dy depends on the pT: They calculate y assuming the particle has the
 67        //   pion mass and assuming that eta=1:
 68        //   dy = 2 * 1/2 * ln [(sqrt(m^2 + (a+1)*pT^2) + a*pT) / (sqrt(m^2 + (a+1)*pT^2) - a*pT)]
 69        //   with a = sinh(1).
 70        //
 71        // sinh(1) = 1.1752012
 72        // m(charged pion)^2 = (139.57 MeV)^2 = 0.019479785 GeV^2
 73        const double sinh1 = 1.1752012;
 74        const double apT  = sinh1 * pT;
 75        const double mPi = 139.57*MeV;
 76        const double root = sqrt(mPi*mPi + (1+sinh1)*pT*pT);
 77        const double dy = std::log((root+apT)/(root-apT));
 78        const double dphi = TWOPI;
 79        _hist_pt->fill(pT, 1.0/(pT*dphi*dy));
 80      }
 81
 82      // Calc sum(Et) from calo particles
 83      const FinalState& etfs = apply<FinalState>(evt, "EtFS");
 84      double sumEt = 0.0;
 85      for (const Particle& p : etfs.particles()) {
 86        sumEt += p.Et();
 87      }
 88      _hist_sumEt->fill(sumEt);
 89      _sumWeightSelected->fill();
 90    }
 91
 92
 93    /// Normalize histos
 94    void finalize() {
 95      scale(_hist_sumEt, crossSection()/millibarn/(4*M_PI*dbl(*_sumWeightSelected)));
 96      scale(_hist_pt, crossSection()/millibarn/dbl(*_sumWeightSelected));
 97      MSG_DEBUG("sumOfWeights()     = " << sumOfWeights());
 98      MSG_DEBUG("_sumWeightSelected = " << dbl(*_sumWeightSelected));
 99    }
100
101    /// @}
102
103
104  private:
105
106    CounterPtr _sumWeightSelected;
107    Profile1DPtr _hist_pt_vs_multiplicity;
108    Histo1DPtr _hist_pt;
109    Histo1DPtr _hist_sumEt;
110
111  };
112
113
114
115  RIVET_DECLARE_ALIASED_PLUGIN(CDF_2009_I817466, CDF_2009_S8233977);
116
117}