1// Samantha Dooling DESY
  2// February 2012
  3//
  4// -*- C++ -*-
  5// =============================
  6//
  7// Ratio of the energy deposited in the pseudorapidity range
  8// -6.6 < eta < -5.2 for events with a charged particle jet
  9//
 10// =============================
 11#include "Rivet/Analysis.hh"
 12#include "Rivet/Projections/FinalState.hh"
 13#include "Rivet/Projections/ChargedFinalState.hh"
 14#include "Rivet/Projections/FastJets.hh"
 15#include "Rivet/Projections/Beam.hh"
 16#include "Rivet/Projections/VetoedFinalState.hh"
 17
 18namespace Rivet {
 19
 20
 21  class CMS_2013_I1218372 : public Analysis {
 22  public:
 23
 24  /// Constructor
 25  CMS_2013_I1218372()
 26    : Analysis("CMS_2013_I1218372")
 27    { }
 28
 29    void init() {
 30
 31      // gives the range of eta and min pT for the final state from which I get the jets
 32      FastJets jetpro (ChargedFinalState((Cuts::etaIn(-2.5, 2.5) && Cuts::pT >=  0.3*GeV)), FastJets::ANTIKT, 0.5);
 33      declare(jetpro, "Jets");
 34
 35      // skip Neutrinos and Muons
 36      VetoedFinalState fsv(FinalState((Cuts::etaIn(-7.0, -4.0))));
 37      fsv.vetoNeutrinos();
 38      fsv.addVetoPairId(PID::MUON);
 39      declare(fsv, "fsv");
 40
 41      FinalState a,b;
 42      a = b;
 43
 44      // for the hadron level selection
 45      VetoedFinalState sfsv;
 46      sfsv.vetoNeutrinos();
 47      sfsv.addVetoPairId(PID::MUON);
 48      declare(sfsv, "sfsv");
 49
 50      //counters
 51      book(passedSumOfWeights, "passedSumOfWeights");
 52      book(inclEflow, "inclEflow");
 53
 54      // Temporary histograms to fill the energy flow for leading jet events.
 55      // Ratios are calculated in finalyze().
 56      int id = 0;
 57      if (isCompatibleWithSqrtS( 900)) id=1;
 58      if (isCompatibleWithSqrtS(2760)) id=2;
 59      if (isCompatibleWithSqrtS(7000)) id=3;
 60      book(_h_ratio, id, 1, 1);
 61      book(_tmp_jet , "TMP/eflow_jet"  ,refData(id, 1, 1));  // Leading jet energy flow in pt
 62      book(_tmp_njet, "TMP/number_jet" ,refData(id, 1, 1)); // Number of events in pt
 63    }
 64
 65
 66    /// Perform the per-event analysis
 67    void analyze(const Event& event) {
 68
 69      // Skip if the event is empty
 70      const FinalState& fsv = apply<FinalState>(event, "fsv");
 71      if (fsv.empty()) vetoEvent;
 72
 73      // ====================== Minimum Bias selection
 74
 75      const FinalState& sfsv = apply<FinalState>(event, "sfsv");
 76      Particles parts = sfsv.particles(cmpMomByRap);
 77      if (parts.empty()) vetoEvent;
 78
 79      // find dymax
 80      double dymax = 0;
 81      int gap_pos  = -1;
 82      for (size_t i = 0; i < parts.size()-1; ++i) {
 83        double dy = parts[i+1].rapidity() - parts[i].rapidity();
 84        if (dy > dymax) {
 85          dymax = dy;
 86          gap_pos = i;
 87        }
 88      }
 89
 90      // calculate mx2 and my2
 91      FourMomentum xmom;
 92      for (int i=0; i<=gap_pos; ++i) {
 93        xmom += parts[i].momentum();
 94      }
 95      double mx2 = xmom.mass2();
 96      if (mx2<0) vetoEvent;
 97
 98      FourMomentum ymom;
 99      for (size_t i=gap_pos+1; i<parts.size(); ++i) {
100        ymom += parts[i].momentum();
101      }
102      double my2 = ymom.mass2();
103      if (my2<0) vetoEvent;
104
105      // calculate xix and xiy and xidd
106      double xix  = mx2 / sqr(sqrtS());
107      double xiy  = my2 / sqr(sqrtS());
108      double xidd = mx2*my2 / sqr(sqrtS()*0.938*GeV);
109
110      // combine the selection: xi cuts
111      bool passedHadronCuts = false;
112      if (isCompatibleWithSqrtS( 900) && (xix > 0.1  || xiy > 0.4 || xidd > 0.5)) passedHadronCuts = true;
113      if (isCompatibleWithSqrtS(2760) && (xix > 0.07 || xiy > 0.2 || xidd > 0.5)) passedHadronCuts = true;
114      if (isCompatibleWithSqrtS(7000) && (xix > 0.04 || xiy > 0.1 || xidd > 0.5)) passedHadronCuts = true;
115      if (!passedHadronCuts) vetoEvent;
116
117      //  ============================== MINIMUM BIAS EVENTS
118
119      // loop over particles to calculate the energy
120      passedSumOfWeights->fill();
121
122      for (const Particle& p : fsv.particles()) {
123        if (-5.2 > p.eta() && p.eta() > -6.6) inclEflow->fill(p.E()/GeV);
124      }
125
126      //  ============================== JET EVENTS
127
128      const FastJets& jetpro = apply<FastJets>(event, "Jets");
129      const Jets& jets = jetpro.jetsByPt(1.0*GeV);
130      if (jets.size()<1) vetoEvent;
131
132      if (fabs(jets[0].eta()) < 2.0) {
133        _tmp_njet->fill(jets[0].pT()/GeV);
134
135        // energy flow
136        for (const Particle& p : fsv.particles()) {
137          if (p.eta() > -6.6 && p.eta() < -5.2) {  // ask for the CASTOR region
138            _tmp_jet->fill(jets[0].pT()/GeV, p.E()/GeV);
139          }
140        }
141      }
142
143    }// analysis
144
145    void finalize() {
146      scale(_tmp_jet, *passedSumOfWeights / *inclEflow);
147      divide(_tmp_jet, _tmp_njet, _h_ratio);
148    }
149
150  private:
151    // counters
152    CounterPtr passedSumOfWeights;
153    CounterPtr inclEflow;
154
155    // histograms
156    Scatter2DPtr _h_ratio;
157    Histo1DPtr   _tmp_jet;
158    Histo1DPtr   _tmp_njet;
159  };
160
161
162  // The hook for the plugin system
163  RIVET_DECLARE_PLUGIN(CMS_2013_I1218372);
164
165}