1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/ChargedFinalState.hh"
  5#include "Rivet/Projections/MissingMomentum.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// UA1 minbias track multiplicities, \f$ p_\perp \f$ and \f$ E_\perp \f$
 11  class UA1_1990_S2044935 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(UA1_1990_S2044935);
 16
 17
 18    /// @name Analysis methods
 19    //@{
 20
 21    /// Book projections and histograms
 22    void init() {
 23      declare(ChargedFinalState((Cuts::etaIn(-5.5, 5.5))), "TriggerFS");
 24      declare(ChargedFinalState((Cuts::etaIn(-2.5, 2.5))), "TrackFS");
 25      const FinalState trkcalofs((Cuts::etaIn(-2.5, 2.5)));
 26      declare(MissingMomentum(trkcalofs), "MET25");
 27      const FinalState calofs((Cuts::etaIn(-6.0, 6.0)));
 28      declare(MissingMomentum(calofs), "MET60");
 29
 30      if (isCompatibleWithSqrtS(63)) {
 31        book(_hist_Pt ,8,1,1);
 32      } else if (isCompatibleWithSqrtS(200)) {
 33        book(_hist_Nch ,1,1,1);
 34        book(_hist_Esigd3p ,2,1,1);
 35        book(_hist_Pt ,6,1,1);
 36        book(_hist_Et ,9,1,1);
 37        book(_hist_Etavg ,12,1,1);
 38      } else if (isCompatibleWithSqrtS(500)) {
 39        book(_hist_Nch ,1,1,2);
 40        book(_hist_Esigd3p ,2,1,2);
 41        book(_hist_Et ,10,1,1);
 42        book(_hist_Etavg ,12,1,2);
 43      } else if (isCompatibleWithSqrtS(900)) {
 44        book(_hist_Nch ,1,1,3);
 45        book(_hist_Esigd3p ,2,1,3);
 46        book(_hist_Pt ,7,1,1);
 47        book(_hist_Et ,11,1,1);
 48        book(_hist_Etavg ,12,1,3);
 49        book(_hist_Esigd3p08 ,3,1,1);
 50        book(_hist_Esigd3p40 ,4,1,1);
 51        book(_hist_Esigd3p80 ,5,1,1);
 52      }
 53      book(_sumwTrig, "TMP/sumwTrig");
 54      // book(_sumwTrig08, "TMP/sumwTrig08");
 55      // book(_sumwTrig40, "TMP/sumwTrig40");
 56      // book(_sumwTrig80, "TMP/sumwTrig80");
 57
 58    }
 59
 60
 61    void analyze(const Event& event) {
 62      // Trigger
 63      const FinalState& trigfs = apply<FinalState>(event, "TriggerFS");
 64      unsigned int n_minus(0), n_plus(0);
 65      for (const Particle& p : trigfs.particles()) {
 66        const double eta = p.eta();
 67        if (inRange(eta, -5.5, -1.5)) n_minus++;
 68        else if (inRange(eta, 1.5, 5.5)) n_plus++;
 69      }
 70      MSG_DEBUG("Trigger -: " << n_minus << ", Trigger +: " << n_plus);
 71      if (n_plus == 0 || n_minus == 0) vetoEvent;
 72      _sumwTrig->fill();
 73
 74      // Use good central detector tracks
 75      const FinalState& cfs = apply<FinalState>(event, "TrackFS");
 76      const double Et25 = apply<MissingMomentum>(event, "MET25").scalarEt();
 77      const double Et60 = apply<MissingMomentum>(event, "MET60").scalarEt();
 78      const unsigned int nch = cfs.size();
 79
 80      // Event level histos
 81      if (!isCompatibleWithSqrtS(63)) {
 82        _hist_Nch->fill(nch);
 83        _hist_Et->fill(Et60/GeV);
 84        _hist_Etavg->fill(nch, Et25/GeV);
 85      }
 86
 87      // Particle/track level histos
 88      const double deta = 2 * 5.0;
 89      const double dphi = TWOPI;
 90      const double dnch_deta = nch/deta;
 91      for (const Particle& p : cfs.particles()) {
 92        const double pt = p.pT();
 93        const double scaled_weight = 1.0/(deta*dphi*pt/GeV);
 94        if (!isCompatibleWithSqrtS(500)) {
 95          _hist_Pt->fill(nch, pt/GeV);
 96        }
 97        if (!isCompatibleWithSqrtS(63)) {
 98          _hist_Esigd3p->fill(pt/GeV, scaled_weight);
 99        }
100        // Also fill for specific dn/deta ranges at 900 GeV
101        if (isCompatibleWithSqrtS(900, 1E-3)) {
102          if (inRange(dnch_deta, 0.8, 4.0)) {
103            //_sumwTrig08 ->fill();
104            _hist_Esigd3p08->fill(pt/GeV, scaled_weight);
105          } else if (inRange(dnch_deta, 4.0, 8.0)) {
106            //_sumwTrig40 ->fill();
107            _hist_Esigd3p40->fill(pt/GeV, scaled_weight);
108          } else {
109            //MSG_WARNING(dnch_deta);
110            if (dnch_deta > 8.0) {
111              //_sumwTrig80 ->fill();
112              _hist_Esigd3p80->fill(pt/GeV, scaled_weight);
113            }
114          }
115        }
116      }
117
118    }
119
120
121    void finalize() {
122      if (_sumwTrig->val() <= 0) {
123        MSG_WARNING("No events passed the trigger!");
124        return;
125      }
126      const double xsec = crossSectionPerEvent();
127      if (!isCompatibleWithSqrtS(63)) {
128        scale(_hist_Nch, 2*xsec/millibarn); ///< Factor of 2 for Nch bin widths?
129        scale(_hist_Esigd3p, xsec/millibarn);
130        scale(_hist_Et, xsec/millibarn);
131      }
132      if (isCompatibleWithSqrtS(900)) {
133        // NB. Ref data is normalised to a fixed value not reproducible from MC.
134        const double scale08 =  (_hist_Esigd3p08->bin(0).area() > 0) ?
135          0.933e5/_hist_Esigd3p08->bin(0).height() : 0;
136        scale(_hist_Esigd3p08, scale08);
137        const double scale40 = (_hist_Esigd3p40->bin(0).area() > 0) ?
138          1.369e5/_hist_Esigd3p40->bin(0).height() : 0;
139        scale(_hist_Esigd3p40, scale40);
140        const double scale80 = (_hist_Esigd3p80->bin(0).area() > 0) ?
141          1.657e5/_hist_Esigd3p80->bin(0).height() : 0;
142        scale(_hist_Esigd3p80, scale80);
143      }
144    }
145
146    //@}
147
148
149  private:
150
151    /// @name Weight counters
152    //@{
153    CounterPtr _sumwTrig; //, _sumwTrig08, _sumwTrig40, _sumwTrig80;
154    //@}
155
156    /// @name Histogram collections
157    //@{
158    Histo1DPtr _hist_Nch;
159    Histo1DPtr _hist_Esigd3p;
160    Histo1DPtr _hist_Esigd3p08;
161    Histo1DPtr _hist_Esigd3p40;
162    Histo1DPtr _hist_Esigd3p80;
163    Profile1DPtr _hist_Pt;
164    Profile1DPtr _hist_Etavg;
165    Histo1DPtr _hist_Et;
166    //@}
167
168  };
169
170
171
172  RIVET_DECLARE_ALIASED_PLUGIN(UA1_1990_S2044935, UA1_1990_I280412);
173
174}