JetForwardJvtTool.cxx 8.65 KB
Newer Older
Mark Hodgkinson's avatar
Mark Hodgkinson committed
1
2
3
///////////////////////// -*- C++ -*- /////////////////////////////

/*
4
  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
Mark Hodgkinson's avatar
Mark Hodgkinson committed
5
6
7
8
9
10
11
12
13
*/

// JetForwardJvtTool.cxx
// Implementation file for class JetForwardJvtTool
// Author: Matt Klein<matthew.henry.klein@cern.ch>
///////////////////////////////////////////////////////////////////

// JetForwardJvtTool includes
#include "JetMomentTools/JetForwardJvtTool.h"
14
15
#include "AsgDataHandles/ReadDecorHandle.h"
#include "AsgDataHandles/WriteDecorHandle.h"
Mark Hodgkinson's avatar
Mark Hodgkinson committed
16

17
18
#include <TString.h>

Mark Hodgkinson's avatar
Mark Hodgkinson committed
19
20
21
22
23
24
25
26
27
// Jet EDM

  ///////////////////////////////////////////////////////////////////
  // Public methods:
  ///////////////////////////////////////////////////////////////////

  // Constructors
  ////////////////
  JetForwardJvtTool::JetForwardJvtTool(const std::string& name) :
28
    AsgTool(name) {
Mark Hodgkinson's avatar
Mark Hodgkinson committed
29
30
31
32
33
34
35
36
37
38
39
40
41
42
  }

  // Destructor
  ///////////////
  JetForwardJvtTool::~JetForwardJvtTool()
  {}

  // Athena algtool's Hooks
  ////////////////////////////
  StatusCode JetForwardJvtTool::initialize()
  {
    ATH_MSG_INFO ("Initializing " << name() << "...");
    if (m_tightOP) m_fjvtThresh = 0.4;
    else m_fjvtThresh = 0.5;
43

44
45
46
47
48
49
50
    ATH_CHECK(m_vertexContainerName.initialize());
    ATH_CHECK(m_trkMETName.initialize());

    if(m_jetContainerName.empty()) {
      ATH_MSG_ERROR("JetForwardJvtTool needs to have its input jet container configured!");
      return StatusCode::FAILURE;
    }
51
    if(!m_orKey.key().empty()) m_orKey = m_jetContainerName + "." + m_orKey.key();
52
53
54
55
56
57
58
59
    m_outKey = m_jetContainerName + "." + m_outKey.key();
    m_isHSKey = m_jetContainerName + "." + m_isHSKey.key();
    m_isPUKey = m_jetContainerName + "." + m_isPUKey.key();
    m_fjvtDecKey = m_jetContainerName + "." + m_fjvtDecKey.key();
    m_widthKey = m_jetContainerName + "." + m_widthKey.key();
    m_jvtMomentKey = m_jetContainerName + "." + m_jvtMomentKey.key();
    m_sumPtsKey = m_jetContainerName + "." + m_sumPtsKey.key();

60
    ATH_CHECK(m_orKey.initialize(!m_orKey.key().empty()));
61
62
63
64
65
66
67
    ATH_CHECK(m_outKey.initialize());
    ATH_CHECK(m_isHSKey.initialize());
    ATH_CHECK(m_isPUKey.initialize());
    ATH_CHECK(m_fjvtDecKey.initialize());
    ATH_CHECK(m_widthKey.initialize());
    ATH_CHECK(m_jvtMomentKey.initialize());
    ATH_CHECK(m_sumPtsKey.initialize());
68

Mark Hodgkinson's avatar
Mark Hodgkinson committed
69
70
71
    return StatusCode::SUCCESS;
  }

72
73
74
75
  StatusCode JetForwardJvtTool::decorate(const xAOD::JetContainer& jetCont) const {
    SG::WriteDecorHandle<xAOD::JetContainer, char> outHandle(m_outKey);
    SG::WriteDecorHandle<xAOD::JetContainer, float> fjvtDecHandle(m_fjvtDecKey);

Mark Hodgkinson's avatar
Mark Hodgkinson committed
76
    getPV();
77
    if (m_recalculateFjvt && jetCont.size() > 0) calculateVertexMomenta(&jetCont);
78
    for(const auto jetF : jetCont) {
79
      outHandle(*jetF) = 1;
80
      if(m_recalculateFjvt) fjvtDecHandle(*jetF) = 0;
Mark Hodgkinson's avatar
Mark Hodgkinson committed
81
82
      if (!forwardJet(jetF)) continue;
      double fjvt = getFJVT(jetF)/jetF->pt();
83
84
      if (fjvt>m_fjvtThresh) outHandle(*jetF) = 0;
      fjvtDecHandle(*jetF) = fjvt;
Mark Hodgkinson's avatar
Mark Hodgkinson committed
85
    }
86
    return StatusCode::SUCCESS;
Mark Hodgkinson's avatar
Mark Hodgkinson committed
87
88
89
  }

  float JetForwardJvtTool::getFJVT(const xAOD::Jet *jet) const {
90
91
92
93
94
    if(!m_recalculateFjvt){
      SG::WriteDecorHandle<xAOD::JetContainer, float> fjvtDecHandle(m_fjvtDecKey);
      return fjvtDecHandle(*jet);
    } 

Mark Hodgkinson's avatar
Mark Hodgkinson committed
95
96
97
98
99
100
101
102
103
104
105
106
107
    TVector2 fjet(-jet->pt()*cos(jet->phi()),-jet->pt()*sin(jet->phi()));
    double fjvt = 0;
    for (size_t pui = 0; pui < m_pileupMomenta.size(); pui++) {
      if (pui==m_pvind) continue;
      double projection = m_pileupMomenta[pui]*fjet/fjet.Mod();
      if (projection>fjvt) fjvt = projection;
    }
    //fjvt += getCombinedWidth(jet);
    return fjvt;
  }

  void JetForwardJvtTool::calculateVertexMomenta(const xAOD::JetContainer *jets) const {
    m_pileupMomenta.clear();
108

109
110
111
112
    SG::ReadHandle<xAOD::MissingETContainer> trkMetHandle(m_trkMETName);
    SG::ReadHandle<xAOD::VertexContainer> vertexContainerHandle(m_vertexContainerName);
    if( !trkMetHandle.isValid() ) {
      ATH_MSG_WARNING("xAOD::MissingETContainer " << m_trkMETName.key() << " is invalid");
113
      return;
Mark Hodgkinson's avatar
Mark Hodgkinson committed
114
    }
115
116
    if( !vertexContainerHandle.isValid() ) {
      ATH_MSG_WARNING("xAOD::VertexContainer " << m_vertexContainerName.key() << " is invalid");
117
      return;
Mark Hodgkinson's avatar
Mark Hodgkinson committed
118
    }
119

120
    for(const auto vx : *vertexContainerHandle) {
Mark Hodgkinson's avatar
Mark Hodgkinson committed
121
122
123
      if(vx->vertexType()!=xAOD::VxType::PriVtx && vx->vertexType()!=xAOD::VxType::PileUp) continue;
      TString vname = "PVTrack_vx";
      vname += vx->index();
124
125
126
      m_pileupMomenta.push_back(\
                                (vx->index()==m_pvind ? \
                                 0:\
127
                                 -(1./m_jetScaleFactor))*TVector2(0.5*(*trkMetHandle)[vname.Data()]->mpx(),0.5*(*trkMetHandle)[vname.Data()]->mpy()));
Mark Hodgkinson's avatar
Mark Hodgkinson committed
128
    }
129

130
    for (const auto jet : *jets) {
Mark Hodgkinson's avatar
Mark Hodgkinson committed
131
132
133
134
135
136
137
138
139
140
      if (!centralJet(jet)) continue;
      int jetvert = getJetVertex(jet);
      if (jetvert>=0) m_pileupMomenta[jetvert] += TVector2(0.5*jet->pt()*cos(jet->phi()),0.5*jet->pt()*sin(jet->phi())); 
    }
  }

  float JetForwardJvtTool::getCombinedWidth(const xAOD::Jet *jet) const {
    float Width = 0;
    float CWidth = 0;
    float ptsum = 0;
141
142
    SG::ReadDecorHandle<xAOD::JetContainer, float> widthHandle(m_widthKey);
    Width = widthHandle(*jet);
Mark Hodgkinson's avatar
Mark Hodgkinson committed
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
    xAOD::JetConstituentVector constvec = jet->getConstituents();
    for (xAOD::JetConstituentVector::iterator it = constvec.begin(); it != constvec.end(); it++) {
      const xAOD::CaloCluster *cl = static_cast<const xAOD::CaloCluster*>((*it)->rawConstituent());
      float secondR = cl->getMomentValue(xAOD::CaloCluster::MomentType::SECOND_R);
      float centermag = cl->getMomentValue(xAOD::CaloCluster::MomentType::CENTER_MAG);
      CWidth+=fabs(cl->pt()*atan(sqrt(secondR)/centermag)*cosh(cl->eta()));
      ptsum += cl->pt();
    }
    CWidth /= ptsum;
    return (CWidth + Width);
  }

  bool JetForwardJvtTool::forwardJet(const xAOD::Jet *jet) const {
    if (fabs(jet->eta())<m_etaThresh) return false;
    if (jet->pt()<m_forwardMinPt || jet->pt()>m_forwardMaxPt) return false;
    return true;
  }

  bool JetForwardJvtTool::centralJet(const xAOD::Jet *jet) const {
    if (fabs(jet->eta())>m_etaThresh) return false;
    if (jet->pt()<m_centerMinPt || (m_centerMaxPt>0 && jet->pt()>m_centerMaxPt)) return false;
164
165
166
167
    if(!m_orKey.key().empty()){
      SG::ReadDecorHandle<xAOD::JetContainer, char> orHandle(m_orKey);
      if(!orHandle(*jet)) return false;
    }
Mark Hodgkinson's avatar
Mark Hodgkinson committed
168
    float jvt = 0;
169
170
    SG::ReadDecorHandle<xAOD::JetContainer, float> jvtMomentHandle(m_jvtMomentKey);
    jvt = jvtMomentHandle(*jet);
Mark Hodgkinson's avatar
Mark Hodgkinson committed
171
172
173
174
175
176
177
    if (jvt>m_centerJvtThresh) return false;
    if (jet->pt()<m_maxStochPt && getDrpt(jet)<m_centerDrptThresh) return false;
    return true;
  }

  int JetForwardJvtTool::getJetVertex(const xAOD::Jet *jet) const {
    std::vector<float> sumpts;
178
179
    SG::ReadDecorHandle<xAOD::JetContainer, std::vector<float> > sumPtsHandle(m_sumPtsKey);
    sumpts = sumPtsHandle(*jet);
Mark Hodgkinson's avatar
Mark Hodgkinson committed
180
181
182
183
184
185
186
187
188
189
190
191
192
    double firstVal = 0;
    int bestMatch = -1;
    for (size_t i = 0; i < sumpts.size(); i++) {
      if (sumpts[i]>firstVal) {
        bestMatch = i;
        firstVal = sumpts[i];
      }
    }
    return bestMatch;
  }

  float JetForwardJvtTool::getDrpt(const xAOD::Jet *jet) const {
    std::vector<float> sumpts;
193
194
    SG::ReadDecorHandle<xAOD::JetContainer, std::vector<float> > sumPtsHandle(m_sumPtsKey);
    sumpts = sumPtsHandle(*jet);
Mark Hodgkinson's avatar
Mark Hodgkinson committed
195
196
197
198
199
200
201
202
203
204
    if (sumpts.size()<2) return 0;

    std::nth_element(sumpts.begin(),sumpts.begin()+sumpts.size()/2,sumpts.end(),std::greater<int>());
    double median = sumpts[sumpts.size()/2];
    std::nth_element(sumpts.begin(),sumpts.begin(),sumpts.end(),std::greater<int>());
    double max = sumpts[0];
    return (max-median)/jet->pt();
  }

  void JetForwardJvtTool::getPV() const {
205

206
    auto vertexContainer = SG::makeHandle (m_vertexContainerName);
207
208
209
210
211
212
    if (!vertexContainer.isValid()){
      ATH_MSG_WARNING("Invalid  xAOD::VertexContainer datahandle");
      return;
    }
    auto vxCont = vertexContainer.cptr();

Mark Hodgkinson's avatar
Mark Hodgkinson committed
213
    m_pvind = 0;
214
    if(vxCont->empty()) {
Mark Hodgkinson's avatar
Mark Hodgkinson committed
215
216
217
      ATH_MSG_WARNING("Event has no primary vertices!");
    } else {
      ATH_MSG_DEBUG("Successfully retrieved primary vertex container");
218
      for(const auto vx : *vxCont) {
Mark Hodgkinson's avatar
Mark Hodgkinson committed
219
220
221
222
223
224
225
        if(vx->vertexType()==xAOD::VxType::PriVtx)
          {m_pvind = vx->index(); break;}
      }
    }
  }

  StatusCode JetForwardJvtTool::tagTruth(const xAOD::JetContainer *jets,const xAOD::JetContainer *truthJets) {
226
227
    SG::WriteDecorHandle<xAOD::JetContainer, bool> isHSHandle(m_isHSKey);
    SG::WriteDecorHandle<xAOD::JetContainer, bool> isPUHandle(m_isPUKey);
228
    for(const auto jet : *jets) {
Mark Hodgkinson's avatar
Mark Hodgkinson committed
229
230
      bool ishs = false;
      bool ispu = true;
231
      for(const auto tjet : *truthJets) {
Mark Hodgkinson's avatar
Mark Hodgkinson committed
232
233
234
        if (tjet->p4().DeltaR(jet->p4())<0.3 && tjet->pt()>10e3) ishs = true;
        if (tjet->p4().DeltaR(jet->p4())<0.6) ispu = false;
      }
235
236
      isHSHandle(*jet)=ishs;
      isPUHandle(*jet)=ispu;
Mark Hodgkinson's avatar
Mark Hodgkinson committed
237
238
239
240
    }
    return StatusCode::SUCCESS;
  }