1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/ChargedFinalState.hh"
  4
  5namespace Rivet {
  6
  7
  8  class ATLAS_2016_I1419652 : public Analysis {
  9  public:
 10
 11    /// Particle types included
 12    enum PartTypes {
 13      k_NoStrange,
 14      k_AllCharged,
 15      kNPartTypes
 16    };
 17
 18    /// Phase space regions
 19    enum RegionID {
 20      k_pt500_nch1_eta25,
 21      k_pt500_nch1_eta08,
 22      kNregions
 23    };
 24
 25    /// Nch cut for each region
 26    const static int nchCut[kNregions];
 27
 28
 29    /// Default constructor
 30    ATLAS_2016_I1419652() : Analysis("ATLAS_2016_I1419652") {}
 31
 32
 33    /// Initialization, called once before running
 34    void init() {
 35
 36      // Projections
 37      const ChargedFinalState cfs500_25((Cuts::etaIn(-2.5, 2.5) && Cuts::pT >=  500.0*MeV));
 38      declare(cfs500_25, "CFS500_25");
 39
 40      const ChargedFinalState cfs500_08((Cuts::etaIn(-0.8, 0.8) && Cuts::pT >=  500.0*MeV));
 41      declare(cfs500_08, "CFS500_08");
 42
 43      for (int iT = 0; iT < kNPartTypes; ++iT)  {
 44        for (int iR = 0; iR < kNregions; ++iR)  {
 45          size_t offset = 8 * iR + 4 * iT;
 46          book(_sumW[iT][iR], "_sumW" + to_str(iT) + to_str(iR));
 47          book(_hist_eta  [iT][iR], offset + 3, 1, 1);
 48          book(_hist_pt   [iT][iR], offset + 4, 1, 1);
 49          book(_hist_nch  [iT][iR], offset + 5, 1, 1);
 50          book(_hist_ptnch[iT][iR], offset + 6, 1, 1);
 51        }
 52      }
 53    }
 54
 55
 56    void analyze(const Event& event) {
 57      string fsName;
 58      for (int iR = 0; iR < kNregions; ++iR)  {
 59        switch (iR) {
 60        case k_pt500_nch1_eta25:  fsName = "CFS500_25"; break;
 61        case k_pt500_nch1_eta08:  fsName = "CFS500_08"; break;
 62        }
 63
 64        const ChargedFinalState& cfs = apply<ChargedFinalState>(event, fsName);
 65
 66        /// What's the benefit in separating this code which is only called from one place?!
 67        fillPtEtaNch(cfs, iR);
 68      }
 69    }
 70
 71
 72
 73    void finalize() {
 74      // Standard histograms
 75      for (int iT = 0; iT < kNPartTypes; ++iT)  {
 76        for (int iR = 0; iR < kNregions; ++iR)  {
 77
 78          double etaRangeSize = -999.0; //intentionally crazy
 79          switch (iR) {
 80            case k_pt500_nch1_eta25  : etaRangeSize = 5.0 ;  break;
 81            case k_pt500_nch1_eta08  : etaRangeSize = 1.6 ;  break;
 82            default: etaRangeSize = -999.0; break; //intentionally crazy
 83          }
 84
 85          if (_sumW[iT][iR]->val() > 0) {
 86            scale(_hist_nch[iT][iR], 1.0/ *_sumW[iT][iR]);
 87            scale(_hist_pt [iT][iR], 1.0/ dbl(*_sumW[iT][iR])/TWOPI/etaRangeSize);
 88            scale(_hist_eta[iT][iR], 1.0/ *_sumW[iT][iR]);
 89          } else {
 90            MSG_WARNING("Sum of weights is zero (!) in type/region: " << iT << " " << iR);
 91          }
 92        }
 93      }
 94    }
 95
 96
 97    /// Helper for collectively filling Nch, pT, eta, and pT vs. Nch histograms
 98    void fillPtEtaNch(const ChargedFinalState& cfs, int iRegion) {
 99
100      // Get number of particles
101      int nch[kNPartTypes];
102      int nch_noStrange = 0;
103      for (const Particle& p : cfs.particles()) {
104        PdgId pdg = p.abspid ();
105        if ( pdg == 3112 || // Sigma-
106             pdg == 3222 || // Sigma+
107             pdg == 3312 || // Xi-
108             pdg == 3334 )  // Omega-
109	        continue;
110	      nch_noStrange++;
111      }
112      nch[k_AllCharged] = cfs.size();
113      nch[k_NoStrange ] = nch_noStrange;
114
115      // Skip if event fails cut for all charged (noStrange will always be less)
116      if (nch[k_AllCharged] < nchCut[iRegion]) return;
117
118      // Fill event weight info
119      _sumW[k_AllCharged][iRegion]->fill();
120      if (nch[k_NoStrange ] >= nchCut[iRegion])	{
121        _sumW[k_NoStrange][iRegion]->fill();
122      }
123
124      // Fill nch
125      _hist_nch[k_AllCharged][iRegion]->fill(nch[k_AllCharged]);
126      if (nch[k_NoStrange ] >= nchCut[iRegion])	{
127        _hist_nch [k_NoStrange][iRegion]->fill(nch[k_NoStrange ]);
128      }
129
130      // Loop over particles, fill pT, eta and ptnch
131      for (const Particle& p : cfs.particles())  {
132        const double pt  = p.pT()/GeV;
133        const double eta = p.eta();
134        _hist_pt     [k_AllCharged][iRegion]->fill(pt , 1.0/pt);
135        _hist_eta    [k_AllCharged][iRegion]->fill(eta);
136        _hist_ptnch  [k_AllCharged][iRegion]->fill(nch[k_AllCharged], pt);
137
138        // Make sure nch cut is passed for nonStrange particles!
139        if (nch[k_NoStrange ] >= nchCut[iRegion])  {
140          PdgId pdg = p.abspid ();
141          if ( pdg == 3112 || // Sigma-
142               pdg == 3222 || // Sigma+
143               pdg == 3312 || // Xi-
144               pdg == 3334 )  // Omega-
145	          continue;
146          // Here we don't have strange particles anymore
147          _hist_pt   [k_NoStrange][iRegion]->fill(pt , 1.0/pt);
148          _hist_eta  [k_NoStrange][iRegion]->fill(eta);
149          _hist_ptnch[k_NoStrange][iRegion]->fill(nch[k_NoStrange], pt);
150        }
151      }
152    }
153
154
155  private:
156
157    CounterPtr _sumW[kNPartTypes][kNregions];
158
159    Histo1DPtr   _hist_nch  [kNPartTypes][kNregions];
160    Histo1DPtr   _hist_pt   [kNPartTypes][kNregions];
161    Histo1DPtr   _hist_eta  [kNPartTypes][kNregions];
162    Profile1DPtr _hist_ptnch[kNPartTypes][kNregions];
163
164  };
165
166
167  // Constants: pT & eta regions
168  const int ATLAS_2016_I1419652::nchCut[] = {1, 1};
169
170
171  RIVET_DECLARE_PLUGIN(ATLAS_2016_I1419652);
172
173}