Rivet Analyses Reference

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165

// Samantha Dooling DESY
// February 2012
//
// -*- C++ -*-
// =============================
//
// Ratio of the energy deposited in the pseudorapidity range
// -6.6 < eta < -5.2 for events with a charged particle jet
//
// =============================
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/ChargedFinalState.hh"
#include "Rivet/Projections/FastJets.hh"
#include "Rivet/Projections/Beam.hh"
#include "Rivet/Projections/VetoedFinalState.hh"

namespace Rivet {


  class CMS_2013_I1218372 : public Analysis {
  public:

  /// Constructor
  CMS_2013_I1218372()
    : Analysis("CMS_2013_I1218372")
    { }

    void init() {

      // gives the range of eta and min pT for the final state from which I get the jets
      FastJets jetpro (ChargedFinalState((Cuts::etaIn(-2.5, 2.5) && Cuts::pT >=  0.3*GeV)), FastJets::ANTIKT, 0.5);
      declare(jetpro, "Jets");

      // skip Neutrinos and Muons
      VetoedFinalState fsv(FinalState((Cuts::etaIn(-7.0, -4.0))));
      fsv.vetoNeutrinos();
      fsv.addVetoPairId(PID::MUON);
      declare(fsv, "fsv");

      FinalState a,b;
      a = b;

      // for the hadron level selection
      VetoedFinalState sfsv;
      sfsv.vetoNeutrinos();
      sfsv.addVetoPairId(PID::MUON);
      declare(sfsv, "sfsv");

      //counters
      book(passedSumOfWeights, "passedSumOfWeights");
      book(inclEflow, "inclEflow");

      // Temporary histograms to fill the energy flow for leading jet events.
      // Ratios are calculated in finalyze().
      int id = 0;
      if (isCompatibleWithSqrtS( 900)) id=1;
      if (isCompatibleWithSqrtS(2760)) id=2;
      if (isCompatibleWithSqrtS(7000)) id=3;
      book(_h_ratio, id, 1, 1);
      book(_tmp_jet , "TMP/eflow_jet"  ,refData(id, 1, 1));  // Leading jet energy flow in pt
      book(_tmp_njet, "TMP/number_jet" ,refData(id, 1, 1)); // Number of events in pt
    }


    /// Perform the per-event analysis
    void analyze(const Event& event) {

      // Skip if the event is empty
      const FinalState& fsv = apply<FinalState>(event, "fsv");
      if (fsv.empty()) vetoEvent;

      // ====================== Minimum Bias selection

      const FinalState& sfsv = apply<FinalState>(event, "sfsv");
      Particles parts = sfsv.particles(cmpMomByRap);
      if (parts.empty()) vetoEvent;

      // find dymax
      double dymax = 0;
      int gap_pos  = -1;
      for (size_t i = 0; i < parts.size()-1; ++i) {
        double dy = parts[i+1].rapidity() - parts[i].rapidity();
        if (dy > dymax) {
          dymax = dy;
          gap_pos = i;
        }
      }

      // calculate mx2 and my2
      FourMomentum xmom;
      for (int i=0; i<=gap_pos; ++i) {
        xmom += parts[i].momentum();
      }
      double mx2 = xmom.mass2();
      if (mx2<0) vetoEvent;

      FourMomentum ymom;
      for (size_t i=gap_pos+1; i<parts.size(); ++i) {
        ymom += parts[i].momentum();
      }
      double my2 = ymom.mass2();
      if (my2<0) vetoEvent;

      // calculate xix and xiy and xidd
      double xix  = mx2 / sqr(sqrtS());
      double xiy  = my2 / sqr(sqrtS());
      double xidd = mx2*my2 / sqr(sqrtS()*0.938*GeV);

      // combine the selection: xi cuts
      bool passedHadronCuts = false;
      if (isCompatibleWithSqrtS( 900) && (xix > 0.1  || xiy > 0.4 || xidd > 0.5)) passedHadronCuts = true;
      if (isCompatibleWithSqrtS(2760) && (xix > 0.07 || xiy > 0.2 || xidd > 0.5)) passedHadronCuts = true;
      if (isCompatibleWithSqrtS(7000) && (xix > 0.04 || xiy > 0.1 || xidd > 0.5)) passedHadronCuts = true;
      if (!passedHadronCuts) vetoEvent;

      //  ============================== MINIMUM BIAS EVENTS

      // loop over particles to calculate the energy
      passedSumOfWeights->fill();

      for (const Particle& p : fsv.particles()) {
        if (-5.2 > p.eta() && p.eta() > -6.6) inclEflow->fill(p.E()/GeV);
      }

      //  ============================== JET EVENTS

      const FastJets& jetpro = apply<FastJets>(event, "Jets");
      const Jets& jets = jetpro.jetsByPt(1.0*GeV);
      if (jets.size()<1) vetoEvent;

      if (fabs(jets[0].eta()) < 2.0) {
        _tmp_njet->fill(jets[0].pT()/GeV);

        // energy flow
        for (const Particle& p : fsv.particles()) {
          if (p.eta() > -6.6 && p.eta() < -5.2) {  // ask for the CASTOR region
            _tmp_jet->fill(jets[0].pT()/GeV, p.E()/GeV);
          }
        }
      }

    }// analysis

    void finalize() {
      scale(_tmp_jet, *passedSumOfWeights / *inclEflow);
      divide(_tmp_jet, _tmp_njet, _h_ratio);
    }

  private:
    // counters
    CounterPtr passedSumOfWeights;
    CounterPtr inclEflow;

    // histograms
    Scatter2DPtr _h_ratio;
    Histo1DPtr   _tmp_jet;
    Histo1DPtr   _tmp_njet;
  };


  // The hook for the plugin system
  RIVET_DECLARE_PLUGIN(CMS_2013_I1218372);

}