1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/ChargedFinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5
  6namespace Rivet {
  7
  8
  9  class ATLAS_2010_CONF_2010_049 : public Analysis {
 10  public:
 11
 12    ATLAS_2010_CONF_2010_049()
 13      : Analysis("ATLAS_2010_CONF_2010_049")
 14    {    }
 15
 16
 17    void init() {
 18      ChargedFinalState cfs((Cuts::etaIn(-1.5, 1.5) && Cuts::pT >=  0.5*GeV));
 19      declare(cfs, "CFS");
 20
 21      FastJets jetsproj6(cfs, JetAlg::ANTIKT, 0.6);
 22      declare(jetsproj6, "Jets6");
 23
 24      FastJets jetsproj4(cfs, JetAlg::ANTIKT, 0.4);
 25      declare(jetsproj4, "Jets4");
 26
 27      /// @todo tmp YOs
 28      for (size_t i=0 ; i<2 ; i++) {
 29        book(_h_xsec[i]       ,1+i, 1, 1);
 30        book(_h_frag_04_06[i] ,3+i, 1, 1);
 31        book(_h_frag_06_10[i] ,3+i, 2, 1);
 32        book(_h_frag_10_15[i] ,3+i, 3, 1);
 33        book(_h_frag_15_24[i] ,3+i, 4, 1);
 34        book(_njets_04_06[i], "njets_04_06_"+to_string(i));
 35        book(_njets_06_10[i], "njets_06_10_"+to_string(i));
 36        book(_njets_10_15[i], "njets_10_15_"+to_string(i));
 37        book(_njets_15_24[i], "njets_15_24_"+to_string(i));
 38      }
 39    }
 40
 41
 42    void analyze(const Event& event) {
 43      const FastJets & jetsproj6 = apply<FastJets>(event, "Jets6");
 44      const FastJets & jetsproj4 = apply<FastJets>(event, "Jets4");
 45      Jets alljets[2];
 46      alljets[0] = jetsproj6.jetsByPt(Cuts::pT > 4*GeV);
 47      alljets[1] = jetsproj4.jetsByPt(Cuts::pT > 4*GeV);
 48
 49      for (size_t i=0 ; i<2 ; i++) {
 50        Jets jets;
 51
 52        // First we want to make sure that we only use jets within |eta|<0.57
 53        for (const Jet& jet : alljets[i]) {
 54          if (jet.abseta()<0.57) {
 55            jets.push_back(jet);
 56          }
 57        }
 58        for (const Jet& jet : jets) {
 59          const double pTjet = jet.pT();
 60          const double pjet = jet.p3().mod();
 61          _h_xsec[i]->fill(pTjet);
 62          if (pTjet > 24*GeV) continue;
 63          for (const Particle& p : jet.particles()) {
 64            double z = p.p3().mod()/pjet;
 65            if (z >= 1) z = 0.9999; // Make sure that z=1 doesn't go into overflow
 66            if (pTjet > 15*GeV) {
 67              _h_frag_15_24[i]->fill(z);
 68            }
 69            else if (pTjet > 10*GeV) {
 70              _h_frag_10_15[i]->fill(z);
 71            }
 72            else if (pTjet > 6*GeV) {
 73              _h_frag_06_10[i]->fill(z);
 74            }
 75            else {
 76              _h_frag_04_06[i]->fill(z);
 77            }
 78          }
 79          if (pTjet > 15*GeV) {
 80            _njets_15_24[i]->fill();
 81          }
 82          else if (pTjet > 10*GeV) {
 83            _njets_10_15[i]->fill();
 84          }
 85          else if (pTjet > 6*GeV) {
 86            _njets_06_10[i]->fill();
 87          }
 88          else {
 89            _njets_04_06[i]->fill();
 90          }
 91        }
 92      }
 93    }
 94
 95    void finalize() {
 96      for (size_t i=0 ; i<2 ; i++) {
 97        // deta = 2*0.57
 98        scale(_h_xsec[i], crossSection()/microbarn/sumOfWeights()/(2*0.57));
 99        scale(_h_frag_04_06[i], 1./_njets_04_06[i]->val());
100        scale(_h_frag_06_10[i], 1./_njets_06_10[i]->val());
101        scale(_h_frag_10_15[i], 1./_njets_10_15[i]->val());
102        scale(_h_frag_15_24[i], 1./_njets_15_24[i]->val());
103      }
104    }
105
106
107  private:
108
109    Histo1DPtr _h_xsec[2];
110    Histo1DPtr _h_frag_04_06[2];
111    Histo1DPtr _h_frag_06_10[2];
112    Histo1DPtr _h_frag_10_15[2];
113    Histo1DPtr _h_frag_15_24[2];
114    CounterPtr _njets_04_06[2];
115    CounterPtr _njets_06_10[2];
116    CounterPtr _njets_10_15[2];
117    CounterPtr _njets_15_24[2];
118
119  };
120
121
122  RIVET_DECLARE_PLUGIN(ATLAS_2010_CONF_2010_049);
123
124}