processSingleBTagEvent.cxx 8.16 KB
Newer Older
1
2
3
4
5
6
7
#include "processSingleBTagEvent.hh"

// the implementation is the same for every signature, we just have to
// recompile it a few times below. We do it here to avoid having the
// template headers interfere with anything that calls this function.
#include "SingleBTagTools.hh"
#include "SingleBTagConfig.hh"
8
#include "ShallowCopyBJets.hh"
9
10
11
12
13
14
15
16
#include "TruthTools.hh"

#include "xAODRootAccess/TEvent.h"
#include "xAODCore/ShallowCopy.h"
#include "xAODTruth/TruthParticleContainer.h"
#include "xAODJet/JetContainer.h"
#include "xAODTracking/TrackParticleContainer.h"
#include "xAODBTagging/BTaggingUtilities.h"
Dan Guest's avatar
Dan Guest committed
17
#include "xAODTracking/TrackMeasurementValidationContainer.h"
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

#include "PathResolver/PathResolver.h"

namespace {

  void check_rc(StatusCode code) {
    if (!code.isSuccess()) throw std::runtime_error("bad return code");
  }

  constexpr float operator "" _GeV(unsigned long long d) { return d*1e3; }
  constexpr float operator "" _GeV(long double d) { return d*1e3; }

  const xAOD::Vertex* primary(const xAOD::VertexContainer& vertices) {
    if (vertices.size() == 0) {
      throw std::runtime_error("no primary vertices");
    }
    for ( const xAOD::Vertex *vertex : vertices ) {
      if ( vertex->vertexType() == xAOD::VxType::PriVtx ) {
        return vertex;
      }
    }
    // if we find nothing else this should be the beam spot
    return vertices.front();
  }
}

// this is the templated code, the concrete instances are below
template <typename Event>
void processSingleBTagEventImpl(Event& event,
                                const SingleBTagConfig& jobcfg,
                                SingleBTagTools& tools) {
49
50
51
52
53
54
55
56
57
  const xAOD::EventInfo *event_info = nullptr;
  check_rc( event.retrieve(event_info, "EventInfo") );

  // jet event cleaning: do not process events with bad jets
  if (jobcfg.jetclean_option == JetCleanOption::event) {
    bool result = tools.acc.eventClean_looseBad(*event_info);
    // do not process event, if event-level jet cleaning flag is not passed
    if (!result) return;
  }
58
59
60
61
62
63

  const xAOD::JetContainer *raw_jets = nullptr;
  check_rc( event.retrieve(raw_jets, jobcfg.jet_collection) );

  // first loop: add decorations to jets. These have to be done
  // before calibration to be consistent with reconstruction
64
  auto [jets, aux] = shallowCopyBJets(*raw_jets);
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
  for (const xAOD::Jet* uncalib_jet: *jets) {

    // this is more important stuff
    const auto* btag = xAOD::BTaggingUtilities::getBTagging(*uncalib_jet);
    if (!btag) throw std::runtime_error("no btaggingLink");

    if (jobcfg.decorate.jet_aug) {
      tools.jet_augmenter.augment(*btag, *btag);
    }
    if (jobcfg.decorate.soft_muon) {
      tools.muon_augmenter.augment(*btag);
    }
    if (jobcfg.decorate.btag_jes){
      tools.jet_augmenter.augmentBtagJes(*btag, *btag);
    }
    for (const auto& dl2: tools.dl2s) {
      dl2.decorate(*btag);
    }
  }
  if (jobcfg.do_calibration) {
    check_rc(tools.calibration_tool.applyCalibration(*jets));
  }

  // sort jets by descending pt
  // we make a new container first to preserve the indices
  std::vector<const xAOD::Jet*> sorted_jets(jets->begin(), jets->end());
  std::sort(sorted_jets.begin(), sorted_jets.end(),
            [](const auto* j1, const auto* j2) {
              return j1->pt() > j2->pt();
            });

Dan Guest's avatar
Dan Guest committed
96
97
98
99
100
101
102
103
104
105
106
107
  // we need the pv to add the correct hits to jets
  const xAOD::VertexContainer* primary_vertices = nullptr;
  check_rc( event.retrieve(primary_vertices, jobcfg.vertex_collection));
  const xAOD::Vertex* pv = primary(*primary_vertices);
  if (jobcfg.decorate.hits) {
    const xAOD::TrackMeasurementValidationContainer* hits = nullptr;
    check_rc( event.retrieve(hits, "JetAssociatedPixelClusters"));
    for (const xAOD::Jet* jet: sorted_jets) {
      tools.hit_decorator.decorate(*jet, *hits, *pv);
    }
  }

108
109
110
111
112
113
114
  // decorate all tracks
  if ( jobcfg.decorate.track_ambiguity_solver_ranking ) {
    const xAOD::TrackParticleContainer* tracks = nullptr;
    check_rc( event.retrieve(tracks, "InDetTrackParticles"));
    tools.trkAmbiDecorator.decorateAll(*tracks);
  }

115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
  // second loop: select calibrated jets and write out to HDF5
  // we'll need truth particles and event info
  std::vector<const xAOD::TruthParticle*> truth_leptons;
  try {
    const xAOD::TruthParticleContainer *tpc = nullptr;
    check_rc( event.retrieve(tpc, "TruthParticles") );
    std::vector<const xAOD::TruthParticle*> truth_particles(
      tpc->begin(), tpc->end());
    truth_leptons = truth::getLeptonsFromWZ(truth_particles);
    tools.truth_counts.n_successful_truth_reads++;
  } catch (truth::TruthRecordError& err) {
    std::cerr << getEventInfoString(err, *event_info) << std::endl;
    tools.truth_counts.n_failed_truth_reads++;
    return;
  }

131
  // save some primary vertex information on eventinfo
132
  tools.dec.n_primary_vertices(*event_info) = primary_vertices->size();
Dan Guest's avatar
Dan Guest committed
133
  tools.dec.primary_vertex_detector_z(*event_info) = pv->z();
134
135
  tools.dec.primary_vertex_detector_z_uncertainty(*event_info) = (
    std::sqrt(pv->covariancePosition()(2,2)));
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
  const xAOD::TruthEventContainer* truthEventContainer = nullptr;
  const xAOD::TruthVertex* truth_PV = nullptr;
  if ( jobcfg.decorate.track_truth_info ) {
    check_rc( event.retrieve(truthEventContainer, "TruthEvents"));
    // truthEventContainer always has size == 1?
    truth_PV = truthEventContainer->at(0)->truthVertex(0);
  }

  std::vector<const xAOD::Jet*> jets_to_write;
  unsigned int rank = 0;
  for (const xAOD::Jet* calib_jet: sorted_jets) {
    if (truth::is_overlaping_lepton(*calib_jet, truth_leptons, 0.3)) {
      continue;
    }

    tools.dec.jet_rank(*calib_jet) = rank++;
    if (jobcfg.do_calibration){
      // if we're calibrating jets we need to check the JVT again
      float updated_jvt_value= tools.jvttool.updateJvt(*calib_jet);
      tools.dec.jvt(*calib_jet) = updated_jvt_value;
      bool fail_jvt = (
        calib_jet->pt() > 20_GeV &&
        calib_jet->pt() < 60_GeV &&
        std::abs(calib_jet->eta()) < 2.4 &&
        updated_jvt_value < jobcfg.jvt_cut );
      if (fail_jvt) continue;

163
164
      // only do jet-level jet cleaning if not doing event-level jet cleaning
      if (jobcfg.jetclean_option == JetCleanOption::jet && !tools.cleaning_tool.keep(*calib_jet)) continue;
165
166
167
168
169
170
171
172
173
    }
    else {
      tools.dec.jvt(*calib_jet) = NAN;
    }

    if (calib_jet->pt() < jobcfg.pt_cut || std::abs(calib_jet->eta()) > 2.5) {
      continue;
    }

174
    // write out tracks associated with jets and collect jets for output
175
176
177
178
179
    for (auto& tracktool: tools.tracks ) {
      const xAOD::Jet* uncalib_jet = raw_jets->at(calib_jet->index());
      auto tracks = tracktool.selector.get_tracks(*uncalib_jet);
      for (const auto& track: tracks) {
        if ( jobcfg.decorate.track_sv_info ) {
180
              tools.trkVertexDecorator.decorate(*track, *uncalib_jet, *pv);
181
182
183
184
185
        }
      }
      if ( jobcfg.decorate.track_truth_info ) {
            tools.trkTruthDecorator.decorateAll(tracks, truth_PV);
      }
Dan Guest's avatar
Dan Guest committed
186
187
188
189
190

      // TODO: remove the sort function's reliance on this decoration
      for (const auto& track: tracks) {
        tracktool.track_accessor.augment_with_ip(*track, *uncalib_jet);
      }
191
      sort(tracks.begin(), tracks.end(), tracktool.sort);
Dan Guest's avatar
Dan Guest committed
192

193
      tracktool.writer.write(tracks, *uncalib_jet);
Dan Guest's avatar
Dan Guest committed
194
      tracktool.n_tracks_decorator(*calib_jet) = tracks.size();
195
    }
Dan Guest's avatar
Dan Guest committed
196
197
    jets_to_write.push_back(calib_jet);

198
  }
199
  // write out jets
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
  if (!jets_to_write.empty()) {
    tools.jet_writer.write(jets_to_write, event_info);
  }

}

// Concrete versions of the templated function above. These are the
// ones that are exposed to be used in other files.
void processSingleBTagEvent(xAOD::TEvent& e,
                            const SingleBTagConfig& c,
                            SingleBTagTools& t) {
  processSingleBTagEventImpl(e, c, t);
}

// StoreGateSvc isn't defined in AnalysisBase...
#ifndef XAOD_STANDALONE
void processSingleBTagEvent(StoreGateSvc& e,
                            const SingleBTagConfig& c,
                            SingleBTagTools& t) {
  processSingleBTagEventImpl(e, c, t);
}
#else
// ... but SgTEvent is
void processSingleBTagEvent(asg::SgTEvent& e,
                            const SingleBTagConfig& c,
                            SingleBTagTools& t) {
  processSingleBTagEventImpl(e, c, t);
}
#endif  // XAOD_STANDALONE