JetObjectCollectionMaker.cxx 43.9 KB
Newer Older
1
/*
2
   Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
3
*/
4

5
// $Id: JetObjectCollectionMaker.cxx 809674 2017-08-23 14:10:24Z iconnell $
6
7
#include "TopSystematicObjectMaker/JetObjectCollectionMaker.h"
#include "TopConfiguration/TopConfig.h"
8
#include "TopConfiguration/TreeFilter.h"
9
10
11
12
13
14
15
16
17
18
19
20
21
22
#include "TopEvent/EventTools.h"

#include "xAODJet/JetContainer.h"
#include "xAODJet/JetAuxContainer.h"
#include "xAODEventInfo/EventInfo.h"
#include "xAODCore/ShallowCopy.h"
#include "xAODBase/IParticleHelpers.h"
#include "xAODMissingET/MissingETContainer.h"
#include "PATInterfaces/SystematicsUtil.h"

#include "TopJetSubstructure/TopJetSubstructure.h"
#include "TopJetSubstructure/LargeJetTrimmer.h"
#include "TopJetSubstructure/SubjetMaker.h"

Tomas Dado's avatar
Tomas Dado committed
23
24
25
26
27
28
namespace top {
  JetObjectCollectionMaker::JetObjectCollectionMaker(const std::string& name) :
    asg::AsgTool(name),
    m_config(nullptr),
    m_doFull_JER(false),
    m_doFull_JER_Pseudodata(false),
29
    m_doOnly_JER_largeR(false),
Tomas Dado's avatar
Tomas Dado committed
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
    m_isMC(false),
    m_doMultipleJES(false),

    m_specifiedSystematics(),
    m_specifiedSystematicsLargeR(),
    m_specifiedSystematicsTrackJets(),
    m_recommendedSystematics(),
    m_nominalSystematicSet(),

    m_jetCalibrationTool("JetCalibrationTool"),
    m_jetCalibrationToolLargeR("JetCalibrationToolLargeR"),

    m_jetUncertaintiesTool("JetUncertaintiesTool"),
    m_jetUncertaintiesToolReducedNPScenario1("JetUncertaintiesToolReducedNPScenario1"),
    m_jetUncertaintiesToolReducedNPScenario2("JetUncertaintiesToolReducedNPScenario2"),
    m_jetUncertaintiesToolReducedNPScenario3("JetUncertaintiesToolReducedNPScenario3"),
    m_jetUncertaintiesToolReducedNPScenario4("JetUncertaintiesToolReducedNPScenario4"),
    m_jetUncertaintiesToolLargeR("JetUncertaintiesToolLargeR"),
48
    m_FFJetSmearingTool("FFJetSmearingTool"),
Tomas Dado's avatar
Tomas Dado committed
49
50

    m_jetUpdateJvtTool("JetUpdateJvtTool"),
51
    m_jetSelectfJvtTool("JetSelectfJvtTool"),
Tomas Dado's avatar
Tomas Dado committed
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71

    m_jetSubstructure(nullptr),

    m_systMap_AllNP(),
    m_systMap_ReducedNPScenario1(),
    m_systMap_ReducedNPScenario2(),
    m_systMap_ReducedNPScenario3(),
    m_systMap_ReducedNPScenario4(),
    m_systMap_LargeR() {
    declareProperty("config", m_config);

    declareProperty("JetCalibrationTool", m_jetCalibrationTool);
    declareProperty("JetCalibrationToolLargeR", m_jetCalibrationToolLargeR);

    declareProperty("JetUncertaintiesTool", m_jetUncertaintiesTool);
    declareProperty("JetUncertaintiesToolReducedNPScenario1", m_jetUncertaintiesToolReducedNPScenario1);
    declareProperty("JetUncertaintiesToolReducedNPScenario2", m_jetUncertaintiesToolReducedNPScenario2);
    declareProperty("JetUncertaintiesToolReducedNPScenario3", m_jetUncertaintiesToolReducedNPScenario3);
    declareProperty("JetUncertaintiesToolReducedNPScenario4", m_jetUncertaintiesToolReducedNPScenario4);
    declareProperty("JetUncertaintiesToolLargeR", m_jetUncertaintiesToolLargeR);
72
    declareProperty("FFJetSmearingTool", m_FFJetSmearingTool);
Tomas Dado's avatar
Tomas Dado committed
73
74
75
76

    declareProperty("JetUpdateJvtTool", m_jetUpdateJvtTool);

    declareProperty("TruthJetCollectionForHSTagging", m_truthJetCollForHS = "AntiKt4TruthJets");
77
78
  }

Tomas Dado's avatar
Tomas Dado committed
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
  StatusCode JetObjectCollectionMaker::initialize() {
    ATH_MSG_INFO(" top::JetObjectCollectionMaker initialize");

    ///-- Lets do the nominal systematics --///
    ///-- JetObjectCollectionMaker is a little different from the others --///
    m_specifiedSystematics.push_back(CP::SystematicSet());
    m_specifiedSystematicsLargeR.push_back(CP::SystematicSet());
    m_specifiedSystematicsTrackJets.push_back(CP::SystematicSet());
    m_recommendedSystematics.push_back(CP::SystematicSet());


    top::check(m_jetCalibrationTool.retrieve(), "Failed to retrieve JetCalibrationTool");
    if (m_config->useLargeRJets()) {
      top::check(m_jetCalibrationToolLargeR.retrieve(),
                 "Failed to retrieve JetCalibrationToolLargeR");
94
95
96
      if(m_config->largeRJESJMSConfig() != "UFOSDMass"){
	top::check(m_jetUncertaintiesToolLargeR.retrieve(),
		   "Failed to retrieve JetUncertaintiesToolLargeR");
97
98
        top::check(m_FFJetSmearingTool.retrieve(),
                   "Failed to retrieve FFJetSmearingTool");
99
      }
100
    }
Tomas Dado's avatar
Tomas Dado committed
101

102
103
104
105
    if (m_config->getDerivationStream() == "PHYS") {
      m_truthJetCollForHS = "AntiKt4TruthDressedWZJets";
    }

106
    ///-- Small-R JER (Pseudo-)Data Smearing Config --///
Tomas Dado's avatar
Tomas Dado committed
107
108
109
110
111
112
    if (m_config->jetJERSmearingModel() == "Full" || m_config->jetJERSmearingModel() == "All") m_doFull_JER = true;
    if (m_config->jetJERSmearingModel() == "Simple") m_doFull_JER = false;
    if (m_config->jetJERSmearingModel() == "Full_PseudoData" ||
        m_config->jetJERSmearingModel() == "All_PseudoData") m_doFull_JER_Pseudodata = true;
    else m_doFull_JER_Pseudodata = false;

113
114
115
116
117
118
119
120
    ///-- Large-R JER (Pseudo-)Data Smearing Config --///
    if ((!m_config->isMC() || m_config->largeRSysts_TreatMCasPseudodata()) 
        && (m_config->largeRJetUncertainties_NPModel()).find("_SimpleJER_") == std::string::npos) m_doOnly_JER_largeR = true;
    else m_doOnly_JER_largeR = false;
    bool skip_systs = false;
    if (m_config->largeRSysts_TreatMCasPseudodata() 
        && (m_config->largeRJetUncertainties_NPModel()).find("_SimpleJER_") != std::string::npos) skip_systs = true;
    
Tomas Dado's avatar
Tomas Dado committed
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142

    m_isMC = m_config->isMC();
    m_doMultipleJES = m_config->doMultipleJES();


    if (m_isMC || m_doFull_JER) {
      if (!m_doMultipleJES) {
        top::check(m_jetUncertaintiesTool.retrieve(), "Failed to retrieve JetUncertaintiesTool");
      } else {
        top::check(
          m_jetUncertaintiesToolReducedNPScenario1.retrieve(),
          "Failed to retrieve JetUncertaintiesToolReducedNPScenario1");
        top::check(
          m_jetUncertaintiesToolReducedNPScenario2.retrieve(),
          "Failed to retrieve JetUncertaintiesToolReducedNPScenario2");
        top::check(
          m_jetUncertaintiesToolReducedNPScenario3.retrieve(),
          "Failed to retrieve JetUncertaintiesToolReducedNPScenario3");
        top::check(
          m_jetUncertaintiesToolReducedNPScenario4.retrieve(),
          "Failed to retrieve JetUncertaintiesToolReducedNPScenario4");
      }
143
    }
144
145


Tomas Dado's avatar
Tomas Dado committed
146
    top::check(m_jetUpdateJvtTool.retrieve(), "Failed to retrieve JetUpdateJvtTool");
147
    //fJVT tool isn't setup unless requested
148
    if (m_config->doForwardJVTinMET() || m_config->getfJVTWP() != "None") {
149
150
      top::check(m_jetSelectfJvtTool.retrieve(), "Failed to retrieve JetSelectfJvtTool");
    }
Tomas Dado's avatar
Tomas Dado committed
151
152
153
154
155
156
    // Take this from the TopConfiguration
    // A blank vector will setup all systematics

    const std:: string& syststr = m_config->systematics();
    std::set<std::string> syst, systLargeR;

157

Tomas Dado's avatar
Tomas Dado committed
158
159
160
    if (!m_config->isSystNominal(syststr) && !m_config->isSystAll(syststr)) {
      bool ok = m_config->getSystematicsList(syststr, syst);
      bool okLargeR = m_config->getSystematicsList(syststr, systLargeR);
161
162

     
Tomas Dado's avatar
Tomas Dado committed
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
      if (!ok || !okLargeR) {
        ATH_MSG_ERROR(" top::JetObjectCollectionMaker could not determine systematic list");
        return StatusCode::FAILURE;
      }
      //here the idea is that if the user specifies AllXXX, we leave syst as an empty string, so that all recommended CP
      // systematics are then used
      if (m_config->contains(syst, "AllJets")) {
        syst.clear();
        systLargeR.clear();
      }
      if (m_config->contains(syst, "AllSmallRJets")) {
        syst.clear();
      }
      if (m_config->contains(systLargeR, "AllLargeRJets")) {
        systLargeR.clear();
      }
179
    }
Tomas Dado's avatar
Tomas Dado committed
180
181
182
183


    ///-- JES systematics --///
    if (m_isMC || m_doFull_JER) {
184
      std::string allNP(""),
Tomas Dado's avatar
Tomas Dado committed
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
      np1("SR_Scenario1_"), np2("SR_Scenario2_"), np3("SR_Scenario3_"), np4("SR_Scenario4_");

      bool onlyJER = ((!m_isMC) && m_doFull_JER) || (m_isMC && m_doFull_JER_Pseudodata);

      if (!m_doMultipleJES) {
        addSystematics(syst, m_jetUncertaintiesTool->recommendedSystematics(), m_systMap_AllNP, allNP, false, onlyJER);
      } else {
        addSystematics(syst,
                       m_jetUncertaintiesToolReducedNPScenario1->recommendedSystematics(), m_systMap_ReducedNPScenario1, np1, false,
                       onlyJER);
        addSystematics(syst,
                       m_jetUncertaintiesToolReducedNPScenario2->recommendedSystematics(), m_systMap_ReducedNPScenario2, np2, false,
                       onlyJER);
        addSystematics(syst,
                       m_jetUncertaintiesToolReducedNPScenario3->recommendedSystematics(), m_systMap_ReducedNPScenario3, np3, false,
                       onlyJER);
        addSystematics(syst,
                       m_jetUncertaintiesToolReducedNPScenario4->recommendedSystematics(), m_systMap_ReducedNPScenario4, np4, false,
                       onlyJER);
      }
205
    }
206
207


208
    ///-- Large-R JES/JER/JMS/JMR systematics --///
Tomas Dado's avatar
Tomas Dado committed
209
    CP::SystematicSet largeRsysts;
210
211
    if ((m_config->isMC() || m_doOnly_JER_largeR) && m_config->useLargeRJets()) {
      if (m_config->largeRJESJMSConfig() == "CombMass") { // Only CombMass is supported for large-R JES/JER/JMS/JMR systematics at the moment
Tomas Dado's avatar
Tomas Dado committed
212
        largeRsysts.insert(m_jetUncertaintiesToolLargeR->recommendedSystematics());
213
        largeRsysts.insert(m_FFJetSmearingTool->recommendedSystematics());
Tomas Dado's avatar
Tomas Dado committed
214
215
      } else {
        ATH_MSG_WARNING(
216
          "TA Mass & Calo Mass & TCCMass & UFO SD Mass are not supported for large-R jet uncertainties at the moment. Large-R jet systemtatics skipped!");
Tomas Dado's avatar
Tomas Dado committed
217
      }
218
    }
219

220
 
Tomas Dado's avatar
Tomas Dado committed
221
222
223
224
225
    ///-- Large R jet tagger scale factor uncertainties -- ///
    if (m_config->isMC() && m_config->useLargeRJets()) {
      for (const std::pair<std::string, std::string>& name : m_config->boostedTaggerSFnames()) {
        ToolHandle<ICPJetUncertaintiesTool> tmp_SF_uncert_tool("JetSFuncert_" + name.first);
        if (tmp_SF_uncert_tool.retrieve()) {
226
227
228
229
230
231
	  
	  m_tagSFUncorrelatedSystematics[name.first].clear();
	  CP::SystematicSet correlatedSys,uncorrelatedSys;
	  const CP::SystematicSet& recommendedSys = tmp_SF_uncert_tool->recommendedSystematics();
	  
	  for (const CP::SystematicVariation& sys : recommendedSys) {
232
233
234
	    // Splitting uncertainties into two categories
	    // correlated will get a special tree
	    // uncorrelated SFs will be stored in the nominal tree
235
236
	    bool res = ((sys.name().find("TopTag") == std::string::npos) &&
			(sys.name().find("WTag") == std::string::npos) &&
237
			(sys.name().find("ZTag") == std::string::npos) &&
238
			(sys.name().find("JetTag") == std::string::npos) &&
239
			(sys.name().find("bTag") == std::string::npos));   
240
	    res ? correlatedSys.insert(sys) : uncorrelatedSys.insert(sys);
241
242
243
	  }
	  
	  m_tagSFUncorrelatedSystematics[name.first] = CP::make_systematics_vector(uncorrelatedSys);
244
245
246
	  for(const CP::SystematicSet& sys : m_tagSFUncorrelatedSystematics[name.first]) {
	    if(sys.name()!="") m_tagSFSysNames[name.first].push_back(name.first + "_" + sys.name());
	  }
247
248
249

	  
          largeRsysts.insert(correlatedSys);
Tomas Dado's avatar
Tomas Dado committed
250
251
          m_tagSFuncertTool[name.first] = tmp_SF_uncert_tool;
        }
252
      }
253
254
      m_config->setBoostedTaggersSFSysNames(m_tagSFSysNames);
      
255
    }
256
    
Tomas Dado's avatar
Tomas Dado committed
257
258

    // add the merged set of systematics for large-R jets including the tagging SF systs
259
260
261
262
    if ((m_config->isMC() || m_doOnly_JER_largeR) && (m_config->useLargeRJets()) && (!skip_systs)) {
      std::string allNPlargeR("");     
      addSystematics(systLargeR, largeRsysts, m_systMap_LargeR, allNPlargeR, true, m_doOnly_JER_largeR);
    }
263
    
Tomas Dado's avatar
Tomas Dado committed
264
    ///-- Large R jet substructure --///
Tomas Dado's avatar
Tomas Dado committed
265
266
267
268
269
    if (m_config->jetSubstructureName() == "Trimmer") m_jetSubstructure.reset(new top::LargeJetTrimmer);

    if (m_config->jetSubstructureName() == "SubjetMaker") m_jetSubstructure.reset(new top::SubjetMaker);

    ///-- Large R jet truth labeling --///
270
    m_jetTruthLabelingTool = nullptr;
Tomas Dado's avatar
Tomas Dado committed
271
    if (m_config->isMC() && m_config->useLargeRJets()) {
272
      m_jetTruthLabelingTool = std::unique_ptr<JetTruthLabelingTool>(new JetTruthLabelingTool("JetTruthLabeling"));
273
274
      // For DAOD_PHYS we need to pass few more arguments as it uses TRUTH3
      if (m_config->getDerivationStream() == "PHYS") {
275
276
277
        top::check(m_jetTruthLabelingTool->setProperty("UseTRUTH3", true), "Failed to set UseTRUTH3 for m_jetTruthLabelingTool");
        top::check(m_jetTruthLabelingTool->setProperty("TruthBosonContainerName", "TruthBoson"), "Failed to set truth container name for m_jetTruthLabelingTool");
        top::check(m_jetTruthLabelingTool->setProperty("TruthTopQuarkContainerName", "TruthTop"), "Failed to set truth container name for m_jetTruthLabelingTool");
278
      }
279
      top::check(m_jetTruthLabelingTool->initialize(), "Failed to initialize m_jetTruthLabelingTool");
280
    }
281

282

Tomas Dado's avatar
Tomas Dado committed
283
284
285
286
287
288
    // set the systematics list
    m_config->systematicsJets(specifiedSystematics());
    m_config->systematicsLargeRJets(specifiedSystematicsLargeR());
    m_config->systematicsTrackJets(m_specifiedSystematicsTrackJets);

    ///-- DL1 Decoration --///
289
290
291
292
293
294
295
296
    for (const std::pair<std::string, std::string>& algo : m_config->bTagAlgo_selToolNames()) {
      m_btagSelToolsDL1Decor[algo.first] = algo.second;
      top::check(m_btagSelToolsDL1Decor[algo.first].retrieve(), "Failed to retrieve " + algo.first + " btagging selector for " + m_config->sgKeyJets() + ". This is required for b-tagging score decorations in EventSaver!");
      if (DLx.count(algo.first) == 0) {
        DLx.emplace(algo.first, SG::AuxElement::Decorator<float>("AnalysisTop_" + algo.first));
      }
    }

297
    if (m_config->useTrackJets()) {
298
299
300
301
302
303
304
      for (const std::pair<std::string, std::string>& algo : m_config->bTagAlgo_selToolNames_trkJet()) {
        m_btagSelToolsDL1Decor_trkJet[algo.first] = algo.second;
        top::check(m_btagSelToolsDL1Decor_trkJet[algo.first].retrieve(), "Failed to retrieve " + algo.first + " btagging selector for " + m_config->sgKeyTrackJets() + ". This is required for b-tagging score decorations in EventSaver!");
        if (DLx.count(algo.first) == 0) {
          DLx.emplace(algo.first, SG::AuxElement::Decorator<float>("AnalysisTop_" + algo.first));
        }
      }
305
306
    }

Tomas Dado's avatar
Tomas Dado committed
307
308
309
310
311
312
313
    // initialize boosted jet taggers -- we have to do it here instead pf TopObjectSelectionTools
    // because we have to apply tagger inbetween JES uncertainty tool and the tagging SF tool
    if (m_config->useLargeRJets()) {
      for (const std::pair<std::string, std::string>& name : m_config->boostedJetTaggers()) {
        std::string fullName = name.first + "_" + name.second;
        m_boostedJetTaggers[fullName] = ToolHandle<IJetSelectorTool>(fullName);
        top::check(m_boostedJetTaggers[fullName].retrieve(), "Failed to retrieve " + fullName);
314
315
      }
    }
Tomas Dado's avatar
Tomas Dado committed
316
317

    return StatusCode::SUCCESS;
318
319
  }

Tomas Dado's avatar
Tomas Dado committed
320
321
322
323
  StatusCode JetObjectCollectionMaker::executeJets(bool executeNominal) {
    bool isLargeR(false);

    return execute(isLargeR, executeNominal);
324
325
  }

Tomas Dado's avatar
Tomas Dado committed
326
327
328
329
  StatusCode JetObjectCollectionMaker::executeLargeRJets(bool executeNominal) {
    bool isLargeR(true);

    return execute(isLargeR, executeNominal);
330
331
  }

Tomas Dado's avatar
Tomas Dado committed
332
333
334
335
336
337
338
339
  StatusCode JetObjectCollectionMaker::execute(const bool isLargeR, bool executeNominal) {
    ///-- Run nominal first, if executing nominal
    if (executeNominal) {
      // decorating the HS jets with truth info on which are HS jets
      if (!isLargeR & m_isMC) {
        top::check(
          decorateHSJets(),
          "Failed to decorate jets with truth info of which are HS - this is needed for JVT scale-factors!");
340
341
342
        if (m_isMC && m_config->jetResponseMatchingDeltaR() > 0) {
          top::check(decorateMatchedTruth(), "Failed to decorate matched jet");
        }
Tomas Dado's avatar
Tomas Dado committed
343
344
345
346
347
348
349
350
351
352
353
      }

      // Decorate the DL1 variable
      top::check(decorateDL1(), "Failed to decorate jets with DL1 b-tagging discriminant");

      ///-- First calibrate the nominal jets, everything else comes from this, so let's only do it once not 3000 times
      // --///
      top::check(calibrate(isLargeR), "Failed to calibrate jets");

      ///-- Return after calibrating the nominal --///
      return StatusCode::SUCCESS;
354
    }
Tomas Dado's avatar
Tomas Dado committed
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374

    ///-- Systematics from here --///

    ///-- JES, JER regular atk4 for now --///
    if (!isLargeR) {
      ///-- JES --///
      if (m_isMC || m_doFull_JER) {
        if (!m_doMultipleJES) {
          top::check(applySystematic(m_jetUncertaintiesTool, m_systMap_AllNP), "Failed to apply JES");
        }
        if (m_doMultipleJES) {
          top::check(applySystematic(m_jetUncertaintiesToolReducedNPScenario1,
                                     m_systMap_ReducedNPScenario1), "Failed to apply JES");
          top::check(applySystematic(m_jetUncertaintiesToolReducedNPScenario2,
                                     m_systMap_ReducedNPScenario2), "Failed to apply JES");
          top::check(applySystematic(m_jetUncertaintiesToolReducedNPScenario3,
                                     m_systMap_ReducedNPScenario3), "Failed to apply JES");
          top::check(applySystematic(m_jetUncertaintiesToolReducedNPScenario4,
                                     m_systMap_ReducedNPScenario4), "Failed to apply JES");
        }
375
      }
Tomas Dado's avatar
Tomas Dado committed
376
    } else {
377
378
379
      // tag calibrated (nominal) jets -- the tagging information will be available
      // for systematically-shifted shallow copies as well
      top::check(tagNominalLargeRJets(), "Failed to tag large-R jets");
380
      if (m_isMC || m_doOnly_JER_largeR) {
381
	top::check(applyTaggingSFSystematic(),"Failed to apply large-R tagging SFs syst.");
Tomas Dado's avatar
Tomas Dado committed
382
383
        top::check(applySystematic(m_jetUncertaintiesToolLargeR, m_systMap_LargeR,
                                   true), "Failed to apply large-R syst.");
384
385
      }
    }
Tomas Dado's avatar
Tomas Dado committed
386
387

    return StatusCode::SUCCESS;
388
389
  }

Tomas Dado's avatar
Tomas Dado committed
390
391
392
  StatusCode JetObjectCollectionMaker::calibrate(const bool isLargeR) {
    ///-- Get base jets from xAOD --///
    std::string sgKey = isLargeR ? m_config->sgKeyLargeRJets() : m_config->sgKeyJets();
393

Tomas Dado's avatar
Tomas Dado committed
394
395
396
397
398
    const xAOD::JetContainer* xaod(nullptr);
    top::check(evtStore()->retrieve(xaod, sgKey), "Failed to retrieve Jets");

    ///-- Shallow copy of the xAOD --///
    std::pair< xAOD::JetContainer*, xAOD::ShallowAuxContainer* > shallow_xaod_copy = xAOD::shallowCopyContainer(*xaod);
399

Tomas Dado's avatar
Tomas Dado committed
400
401
402
403
    ///-- Loop over the xAOD Container --///
    for (auto jet : *(shallow_xaod_copy.first)) {
      if (!isLargeR) { ///-- small-R jets used in most analyses --///
        ///-- Calibrate jet --///
404

Tomas Dado's avatar
Tomas Dado committed
405
        top::check(m_jetCalibrationTool->applyCalibration(*jet), "Failed to applyCalibration");
406

Tomas Dado's avatar
Tomas Dado committed
407
408
409
        // only multiply by JSF and bJSF if one of them != 1.0 (used by top mass analysis)
        float JSF = m_config->JSF();
        float bJSF = m_config->bJSF();
410

Tomas Dado's avatar
Tomas Dado committed
411
412
413
414
415
        if (JSF != 1.0 || bJSF != 1.0) {
          int truthflav = -1;
          if (jet->isAvailable<int>("PartonTruthLabelID")) {
            jet->getAttribute("PartonTruthLabelID", truthflav);
          }
416

Tomas Dado's avatar
Tomas Dado committed
417
418
          xAOD::JetFourMom_t jet_p4 = jet->jetP4() * JSF;
          if (truthflav == 5) jet_p4 = jet_p4 * bJSF;
419

Tomas Dado's avatar
Tomas Dado committed
420
          jet->setJetP4(jet_p4);
421
        }
Tomas Dado's avatar
Tomas Dado committed
422
        // end application JSF/bJSF
423

Tomas Dado's avatar
Tomas Dado committed
424
425

        top::check(decorateBJets(*jet), "Failed to decorate if b-jet");
426
427
      }

Tomas Dado's avatar
Tomas Dado committed
428
429
430
      if (isLargeR && m_jetSubstructure.get() != nullptr) {
        m_jetSubstructure->correctJet(*jet);
      }
431

Tomas Dado's avatar
Tomas Dado committed
432
433
434
435
436
437
      if (isLargeR) {
        ///-- Calibrate jet --///

        float tau3 = jet->getAttribute<float>("Tau3_wta");
        float tau2 = jet->getAttribute<float>("Tau2_wta");
        float tau1 = jet->getAttribute<float>("Tau1_wta");
438
439
440
        float ECF1 = jet->getAttribute<float>("ECF1");
        float ECF2 = jet->getAttribute<float>("ECF2");
        float ECF3 = jet->getAttribute<float>("ECF3");
Tomas Dado's avatar
Tomas Dado committed
441
442
443
444
445

        jet->auxdecor<float>("Tau32_wta") = fabs(tau2) > 1.e-6 ? (tau3 / tau2) : -999;  // 999 to match
                                                                                        // JetSubStructureMomentTools/NSubjettinessRatiosTool
        jet->auxdecor<float>("Tau21_wta") = fabs(tau1) > 1.e-6 ? (tau2 / tau1) : -999;  // 999 to match
                                                                                        // JetSubStructureMomentTools/NSubjettinessRatiosTool
446
447
448
449
450
        // Same definition as in JetSubStructureMomentTools/Root/EnergyCorrelatorRatiosTool.cxx
        jet->auxdecor<float>("D2") = (ECF2 > 1e-8) ? (ECF3*ECF1*ECF1*ECF1) / (ECF2*ECF2*ECF2) : -999;
        jet->auxdecor<float>("C2") = (ECF2 > 1e-8) ? (ECF3*ECF1) / (ECF2*ECF2) : -999;
        jet->auxdecor<float>("E3") = (ECF1 > 1e-8) ? ECF3 / (ECF1*ECF1*ECF1) : -999;
                                                                                      
Tomas Dado's avatar
Tomas Dado committed
451
452
453
454
455
456
457

        top::check(m_jetCalibrationToolLargeR->applyCalibration(*jet), "Failed to applyCalibration");

        ///-- for TA mass or calo mass, the calibrated mass or pt needs special treatment --///
        const std::string calibChoice = m_config->largeRJESJMSConfig();
        if (m_config->isMC()) {
          ///-- Truth labeling required by the large-R jet uncertainties --///
458
          top::check(m_jetTruthLabelingTool->modifyJet(*jet), "Failed to do truth labeling for large-R jet");
Tomas Dado's avatar
Tomas Dado committed
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
          if (calibChoice == "TAMass") {
            xAOD::JetFourMom_t jet_calib_p4;
            jet->getAttribute<xAOD::JetFourMom_t>("JetJMSScaleMomentumTA", jet_calib_p4);
            jet->setJetP4(jet_calib_p4);
          }
        } else { //For data, there's only one config file so special method is required for TA mass and Calos mass
          if (calibChoice == "CaloMass") {
            xAOD::JetFourMom_t jetInsituP4_calo;
            jet->getAttribute<xAOD::JetFourMom_t>("JetInsituScaleMomentumCalo", jetInsituP4_calo);
            jet->setJetP4(jetInsituP4_calo);
          } else if (calibChoice == "TAMass") {
            xAOD::JetFourMom_t jetInsituP4_ta;
            jet->getAttribute<xAOD::JetFourMom_t>("JetInsituScaleMomentumTA", jetInsituP4_ta);
            jet->setJetP4(jetInsituP4_ta);
          }
        }
      }
476

Tomas Dado's avatar
Tomas Dado committed
477
478
479
480
      ///-- Update JVT --///
      if (!isLargeR) {
        jet->auxdecor<float>("AnalysisTop_JVT") = m_jetUpdateJvtTool->updateJvt(*jet);
      }
481
482
    }

Tomas Dado's avatar
Tomas Dado committed
483
484
485
    // Check if the derivation we are running on contains
    // MET_Track (once), if so apply the fJVT decoration
    // if not then don't
486
    if (!isLargeR && (m_config->doForwardJVTinMET() || m_config->getfJVTWP() != "None")) {
487
488
489
490
491
492
      static bool checked_track_MET = false;
      if (!checked_track_MET) {
	if (evtStore()->contains<xAOD::MissingETContainer>("MET_Track")) {
	  m_do_fjvt = true;
	} else {
	  ATH_MSG_ERROR(" Cannot retrieve MET_Track, fJVT values can't be calculated correctly!!"); 
493
	  return StatusCode::FAILURE; 
494
495
	}
	checked_track_MET = true;
496
      }
Tomas Dado's avatar
Tomas Dado committed
497
498
    }
    if (m_do_fjvt) {
499
      top::check(!m_jetSelectfJvtTool->modify(*shallow_xaod_copy.first),
Tomas Dado's avatar
Tomas Dado committed
500
                 "Failed to apply fJVT decoration");
501
502
503
    }

    if (!isLargeR) {
Tomas Dado's avatar
Tomas Dado committed
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
      ///-- Save calibrated jet to TStore --///
      ///-- set links to original objects- needed for MET calculation --///
      // NOTE, if we use one of the b-tagging re-trained collections, we need to load
      // the original uncalibrated jet container to which the b-tagging shallow-copy is pointing to
      const xAOD::JetContainer* xaod_original(nullptr);
      // for small-R jet collections, set links to uncalibrated *original* jets (not BTagging shallow-copy)
      if (m_config->sgKeyJets() != m_config->sgKeyJetsType()) {
        top::check(evtStore()->retrieve(xaod_original,
                                        m_config->sgKeyJetsType()),
                   "Failed to retrieve uncalibrated Jets for METMaker!");
      } else {
        xaod_original = xaod;
      }
      bool setLinks = xAOD::setOriginalObjectLink(*xaod_original, *shallow_xaod_copy.first);
      if (!setLinks) ATH_MSG_ERROR(" Cannot set original object links for jets, MET recalculation may struggle");
519
    }
Tomas Dado's avatar
Tomas Dado committed
520
521
522
523
524

    ///-- Save corrected xAOD Container to StoreGate / TStore --///
    std::string outputSGKey;
    if (!isLargeR) {
      outputSGKey = m_config->sgKeyJetsStandAlone(m_nominalSystematicSet.hash());
525
    }
Tomas Dado's avatar
Tomas Dado committed
526
527
    if (isLargeR) {
      outputSGKey = m_config->sgKeyLargeRJets(m_nominalSystematicSet.hash());
528
529
    }

Tomas Dado's avatar
Tomas Dado committed
530
    std::string outputSGKeyAux = outputSGKey + "Aux.";
531

Tomas Dado's avatar
Tomas Dado committed
532
533
534
535
    xAOD::TReturnCode save = evtStore()->tds()->record(shallow_xaod_copy.first, outputSGKey);
    xAOD::TReturnCode saveAux = evtStore()->tds()->record(shallow_xaod_copy.second, outputSGKeyAux);
    if (!save || !saveAux) {
      return StatusCode::FAILURE;
536
    }
537

Tomas Dado's avatar
Tomas Dado committed
538
    return StatusCode::SUCCESS;
539
540
  }

541
 
542
543
544
545
546
547
548
  StatusCode JetObjectCollectionMaker::applyTaggingSFSystematic() {
    
    ///-- Get calibrated jets --///
    const xAOD::JetContainer* ljets(nullptr);
    top::check(evtStore()->retrieve(ljets, m_config->sgKeyLargeRJets(
                                        m_nominalSystematicSet.hash())), "Failed to retrieve Jets");
    
549
550
    const size_t njets = ljets->size();
    
551
552
553
554
    const std::unordered_map<std::string,std::string>& sfNames = m_config->boostedTaggerSFnames();
   
    for(auto& it : m_tagSFuncertTool) {
      ToolHandle<ICPJetUncertaintiesTool>& tool = it.second;
555
      const std::string& fullName=it.first;   
556
557
558
      const SG::AuxElement::Accessor< char > accRange("passedRangeCheck_" + fullName);
      const std::string sfNameNominal =  sfNames.at(fullName);
      const SG::AuxElement::Accessor< float > accSF(sfNameNominal);
559
560
561
562
      // accessor to retrieve nominal efficiency decoration from BoostedJetTaggers
      std::string taggerName = sfNameNominal;
      taggerName.erase(taggerName.length()-3);
      const SG::AuxElement::Accessor<float> accEff(taggerName + "_efficiency");
563
      
564
      for(const CP::SystematicSet& sys : m_tagSFUncorrelatedSystematics[fullName]) {
565
566
567
568
569
570
571
572

        // shallow copy with nominal tagging SFs and efficiencies
        // JetUncertainties may apply relative variations, so we need to make sure we always apply on
        // shallow copy from nominal tagged jets
        std::pair< xAOD::JetContainer*, xAOD::ShallowAuxContainer* >
          shallow_xaod_copy = xAOD::shallowCopyContainer(*ljets);
        auto shallowJets = std::make_pair(std::unique_ptr<xAOD::JetContainer>{shallow_xaod_copy.first},
          std::unique_ptr<xAOD::ShallowAuxContainer>{shallow_xaod_copy.second});
573
	  
574
575
	top::check(tool->applySystematicVariation(sys), "Failed to applySystematicVariation");
	
576
577
578
	const std::string sfNameShifted = fullName + "_" + sys.name();
      
	for(size_t i=0;i<njets;i++) {
579
	  xAOD::Jet* shallowJet = shallowJets.first->at(i);
580
581
582
583
584
585
	  const xAOD::Jet* jet = ljets->at(i);
      
	  if(accRange.isAvailable(*shallowJet) && accRange(*shallowJet)) {
	    top::check(tool->applyCorrection(*shallowJet), "Failed to applyCorrection");
	    float sf = accSF.isAvailable(*shallowJet) ? accSF(*shallowJet) : -999.;
	    jet->auxdecor<float>(sfNameShifted.c_str()) = sf;
586
587
588
589
590
591
592

            // decorate efficiencies for inefficiency SF variations
            if (sys.name().find("TagEffUnc") != std::string::npos) {
              if (accEff.isAvailable(*shallowJet)) {
                jet->auxdecor<float>(fullName + "_" + sys.name() + "_efficiency") = accEff(*shallowJet);
              }
            }
593
594
595
596
597
598
599
600
	  }
	}
      }
    }
   
    return StatusCode::SUCCESS;
  }

Tomas Dado's avatar
Tomas Dado committed
601
602
603
604
605
606
  StatusCode JetObjectCollectionMaker::applySystematic(ToolHandle<ICPJetUncertaintiesTool>& tool,
                                                       const std::unordered_map<CP::SystematicSet,
                                                                                CP::SystematicSet>& map,
                                                       bool isLargeR) {
    ///-- Get calibrated jets --///
    const xAOD::JetContainer* xaod(nullptr);
607

Tomas Dado's avatar
Tomas Dado committed
608
609
610
611
612
613
614
    if (!isLargeR) {
      top::check(evtStore()->retrieve(xaod, m_config->sgKeyJetsStandAlone(
                                        m_nominalSystematicSet.hash())), "Failed to retrieve Jets");
    } else {
      top::check(evtStore()->retrieve(xaod, m_config->sgKeyLargeRJets(
                                        m_nominalSystematicSet.hash())), "Failed to retrieve Jets");
    }
615
616


617
    ///-- Loop over the systematics --///
Tomas Dado's avatar
Tomas Dado committed
618
619
620
    for (Itr syst = map.begin(); syst != map.end(); ++syst) {
      ///-- Don't do the nominal, we've already done that --///
      if ((*syst).second.hash() != m_nominalSystematicSet.hash()) {
621
622
623
624
625
626
627
628

        ///-- Grab systematic name, check if systematic is JMR-type --///
        bool isJMR = false;
        CP::SystematicSet::iterator itr = (*syst).second.begin();
        std::string systname = itr->name();
        if (systname.find("JMR") != std::string::npos) isJMR = true;
          

Tomas Dado's avatar
Tomas Dado committed
629
630
        ///-- Tell the tool which systematic to use --///
        ///-- Here we use the second, original CP::SystematicSet --///
631
632
633
        if (!isJMR) top::check(tool->applySystematicVariation((*syst).second),"Failed to applySystematicVariation");
        else top::check(m_FFJetSmearingTool->applySystematicVariation((*syst).second),"Failed to applySystematicVariation");
           
634

Tomas Dado's avatar
Tomas Dado committed
635
636
637
638
        if (isLargeR && m_config->isMC()) {
          // for boosted tagging SFs, apply syst variation for all initialized WPs
          for (std::pair<const std::string, ToolHandle<ICPJetUncertaintiesTool> >& tagSF : m_tagSFuncertTool) {
            top::check(tagSF.second->applySystematicVariation((*syst).second), "Failed to applySystematicVariation");
639
          }
640
        }
641

Tomas Dado's avatar
Tomas Dado committed
642
643
644
645
646
647
648
649
650
        ///-- Shallow copy of the xAOD --///
        std::pair< xAOD::JetContainer*, xAOD::ShallowAuxContainer* >
        shallow_xaod_copy = xAOD::shallowCopyContainer(*xaod);

        ///-- Loop over the xAOD Container --///
        for (auto jet : *(shallow_xaod_copy.first)) {
          if (isLargeR && m_config->isMC()) { //JES for large-R jets only exist for MC
            ///-- Only large R jets with the following properties can be calibrated.--///
            bool calibratable_jet = (std::fabs(jet->eta()) <= 2.0
651
                                     && jet->pt() > 175e3  //lower than 200e3 to allow studying on migration
652
653
                                     && jet->pt() < 6500e3
                                     && jet->m() < 700e3);
Tomas Dado's avatar
Tomas Dado committed
654
655
656
657
658
659
660
661
662
            std::string jetCalibrationNameLargeR = m_config->sgKeyLargeRJets();
            if (jetCalibrationNameLargeR.find("TrackCaloCluster") != std::string::npos) { //TCC jet
              calibratable_jet = (jet->m() / jet->pt() <= 1.
                                  && jet->m() / jet->pt() >= 0.
                                  && std::fabs(jet->eta()) <= 2.0
                                  && jet->pt() > 150e3
                                  && jet->pt() < 3000e3);
            }
            if (!calibratable_jet) continue;
663
664
	  

665
666
667
668
669
	    ///-- Apply large-R jet tagging SF uncertainties --///
	    for (std::pair<const std::string, ToolHandle<ICPJetUncertaintiesTool> >& tagSF : m_tagSFuncertTool) {
	      std::string fullName=tagSF.first;
	      std::replace(std::begin(fullName),std::end(fullName),':','_');
	      const SG::AuxElement::Accessor< char > acc("passedRangeCheck_" + fullName);
670

671
672
673
	      if(acc.isAvailable(*jet) && acc(*jet)) {
		top::check(tagSF.second->applyCorrection(*jet), "Failed to applyCorrection");
	      }
674

675
676
            }
	    
Tomas Dado's avatar
Tomas Dado committed
677
          }
678
	  
679

Tomas Dado's avatar
Tomas Dado committed
680
          ///-- Apply Corrrection --///
681
682
          if (!isJMR) top::check(tool->applyCorrection(*jet), "Failed to applyCorrection");
          else top::check(m_FFJetSmearingTool->applyCorrection(*jet), "Failed to applyCorrection");                      
683
684


Tomas Dado's avatar
Tomas Dado committed
685
686
          ///-- Update JVT --///
          if (!isLargeR) jet->auxdecor<float>("AnalysisTop_JVT") = m_jetUpdateJvtTool->updateJvt(*jet);
687
688
689
690
	  
          ///-- Decorate fJVT for systematics too --///
	  // Check if the derivation we are running on contains
	  // MET_Track (once) before applying the fJVT decoration
691
	  if (!isLargeR && (m_config->doForwardJVTinMET() || m_config->getfJVTWP() != "None")) {
692
693
694
695
696
697
	    static bool checked_track_MET = false;
	    if (!checked_track_MET) {
	      if (evtStore()->contains<xAOD::MissingETContainer>("MET_Track")) {
		m_do_fjvt = true;
	      } else {
		ATH_MSG_ERROR(" Cannot retrieve MET_Track, fJVT values can't be calculated correctly!!"); 
698
		return StatusCode::FAILURE;
699
700
701
702
703
704
705
706
707
	      }
	      checked_track_MET = true;
	    }
	  }
	  if (m_do_fjvt) {
	    top::check(!m_jetSelectfJvtTool->modify(*shallow_xaod_copy.first),
		       "Failed to apply fJVT decoration");
	  }

708
        }
709

Tomas Dado's avatar
Tomas Dado committed
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
        ///-- set links to original objects- needed for MET calculation --///
        bool setLinks = xAOD::setOriginalObjectLink(*xaod,
                                                    *shallow_xaod_copy.first);
        if (!setLinks) ATH_MSG_ERROR(" Cannot set original object links"
                                     " for jets, MET recalculation may struggle");

        ///-- Save corrected xAOD Container to StoreGate / TStore --///
        ///-- Here we use the first, AnalysisTop modified CP::SystematicSer --///
        ///-- This allows us to run multiple JES scenarios, which all have the same hash values --///
        std::string outputSGKey;
        if (isLargeR) {
          outputSGKey = m_config->sgKeyLargeRJets((*syst).first.hash());
        } else {
          outputSGKey = m_config->sgKeyJetsStandAlone((*syst).first.hash());
        }
        std::string outputSGKeyAux = outputSGKey + "Aux.";

        xAOD::TReturnCode save = evtStore()->tds()->record(shallow_xaod_copy.first, outputSGKey);
        xAOD::TReturnCode saveAux = evtStore()->tds()->record(shallow_xaod_copy.second, outputSGKeyAux);
        if (!save || !saveAux) {
          return StatusCode::FAILURE;
        }
732
733
      }
    }
Tomas Dado's avatar
Tomas Dado committed
734
735

    return StatusCode::SUCCESS;
736
737
  }

Tomas Dado's avatar
Tomas Dado committed
738
739
740
741
742
743
  StatusCode JetObjectCollectionMaker::executeTrackJets(bool executeNominal) {
    ///-- No calibrations or systematics yet --///
    ///-- Only run this on the nominal execution --///
    if (!executeNominal) return StatusCode::SUCCESS;

    top::check(decorateDL1(true), "Failed to decorate track jets with DL1 b-tagging discriminant");
744

745

Tomas Dado's avatar
Tomas Dado committed
746
    ///-- Just make a shallow copy to keep these in line with everything else --///
747

Tomas Dado's avatar
Tomas Dado committed
748
749
750
751
752
753
754
755
756
    const xAOD::JetContainer* xaod(nullptr);
    top::check(evtStore()->retrieve(xaod, m_config->sgKeyTrackJets()), "Failed to retrieve Jets");

    ///-- Shallow copy of the xAOD --///
    std::pair< xAOD::JetContainer*, xAOD::ShallowAuxContainer* > shallow_xaod_copy = xAOD::shallowCopyContainer(*xaod);

    ///-- set links to original objects --///
    bool setLinks = xAOD::setOriginalObjectLink(*xaod, *shallow_xaod_copy.first);
    if (!setLinks) ATH_MSG_ERROR(" Cannot set original object links for track jets");
757

Tomas Dado's avatar
Tomas Dado committed
758
759
760
    ///-- Save corrected xAOD Container to StoreGate / TStore --///
    std::string outputSGKey = m_config->sgKeyTrackJets(m_config->nominalHashValue());
    std::string outputSGKeyAux = outputSGKey + "Aux.";
761

Tomas Dado's avatar
Tomas Dado committed
762
763
764
765
766
767
768
    xAOD::TReturnCode save = evtStore()->tds()->record(shallow_xaod_copy.first, outputSGKey);
    xAOD::TReturnCode saveAux = evtStore()->tds()->record(shallow_xaod_copy.second, outputSGKeyAux);
    if (!save || !saveAux) {
      return StatusCode::FAILURE;
    }
    return StatusCode::SUCCESS;
  }
769

Tomas Dado's avatar
Tomas Dado committed
770
771
  StatusCode JetObjectCollectionMaker::printoutJets() {
    bool isLargeR(false);
772

Tomas Dado's avatar
Tomas Dado committed
773
774
    return this->printout(isLargeR);
  }
775

Tomas Dado's avatar
Tomas Dado committed
776
777
  StatusCode JetObjectCollectionMaker::printoutLargeRJets() {
    bool isLargeR(true);
778

Tomas Dado's avatar
Tomas Dado committed
779
    return this->printout(isLargeR);
780
781
  }

Tomas Dado's avatar
Tomas Dado committed
782
783
784
785
786
787
788
  StatusCode JetObjectCollectionMaker::printout(const bool isLargeR) {
    ///-- Loop over all systematics --///
    for (auto s : m_specifiedSystematics) {
      std::string sgKey = isLargeR ? m_config->sgKeyLargeRJets() : m_config->sgKeyJets();

      const xAOD::JetContainer* xaod(nullptr);
      top::check(evtStore()->retrieve(xaod, sgKey), "Failed to retrieve Jets");
789

Tomas Dado's avatar
Tomas Dado committed
790
791
792
793
      ATH_MSG_INFO(" Jets with sgKey = " << sgKey);
      for (auto x : *xaod) {
        ATH_MSG_INFO("   Jet pT , eta  = " << x->pt() << " , " << x->eta());
      }
794
    }
Tomas Dado's avatar
Tomas Dado committed
795
796

    return StatusCode::SUCCESS;
797
798
  }

799
800
801
802
803
  std::string JetObjectCollectionMaker::getLargeRModName(const std::string& NPModel) const {
    if (NPModel.find("CategoryReduction")!=std::string::npos) return "CategoryReduction_";
    return NPModel+"_"; 
  }

Tomas Dado's avatar
Tomas Dado committed
804
805
806
  void JetObjectCollectionMaker::addSystematics(const std::set<std::string>& specifiedSystematics,
                                                const CP::SystematicSet& recommendedSysts,
                                                std::unordered_map<CP::SystematicSet, CP::SystematicSet>& map,
807
                                                const std::string& modName, bool isLargeR, bool onlyJER) {
808

Tomas Dado's avatar
Tomas Dado committed
809
810
    ///-- Get the recommended systematics from the tool, in std::vector format --///
    const std::vector<CP::SystematicSet> systList = CP::make_systematics_vector(recommendedSysts);
811
        
Tomas Dado's avatar
Tomas Dado committed
812
813
814
    for (const CP::SystematicSet& s : systList) {
      if (s.size() == 1) {
        CP::SystematicSet::const_iterator ss = s.begin();
815
816

	if(!m_config->getTreeFilter()->filterTree(modName + ss->name())) continue; // Applying tree filter
Tomas Dado's avatar
Tomas Dado committed
817
818
819

        if (onlyJER && ss->name().find("JER") == std::string::npos) continue;

820
        CP::SystematicSet modSet(modName + ss->name());
Tomas Dado's avatar
Tomas Dado committed
821
822
823
824
825
826
827
828
829

        m_recommendedSystematics.push_back(modSet);
        if (!m_config->isSystNominal(m_config->systematics())) {
          if (specifiedSystematics.size() == 0) {
            if (!isLargeR) m_specifiedSystematics.push_back(modSet);
            else m_specifiedSystematicsLargeR.push_back(modSet);
            map.insert(std::make_pair(modSet, s));
          } else if (specifiedSystematics.size() > 0) {
            for (std::string i : specifiedSystematics) {
830
831
	      TreeFilter filter(i);
              if (!filter.filterTree(modSet.name())) {
Tomas Dado's avatar
Tomas Dado committed
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
                if (!isLargeR) m_specifiedSystematics.push_back(modSet);
                else m_specifiedSystematicsLargeR.push_back(modSet);
                map.insert(std::make_pair(modSet, s));
              }
            }
          } // User has specified a systematic
        } // Don't do anything if the user requests nominal only
      } // size() == 1 -- this means that there is a CP::SystematicVariation
    } // Loop over systList

    m_recommendedSystematics.sort();
    m_recommendedSystematics.unique();

    m_specifiedSystematics.sort();
    m_specifiedSystematics.unique();

    m_specifiedSystematicsLargeR.sort();
    m_specifiedSystematicsLargeR.unique();
850
  }
Tomas Dado's avatar
Tomas Dado committed
851
852
853
854
855
856
857
858
859
860
861

  StatusCode JetObjectCollectionMaker::tagNominalLargeRJets() {
    const xAOD::JetContainer* xaod_calibrated_jets(nullptr);

    top::check(evtStore()->retrieve(xaod_calibrated_jets, m_config->sgKeyLargeRJets(
                                      m_nominalSystematicSet.hash())),
               "Failed to retrieve nominal calibrated large-R jets");
    for (const xAOD::Jet* jet : *xaod_calibrated_jets) {
      top::check(tagLargeRJet(*jet), "Failed to tag large-R jet");
    }
    return StatusCode::SUCCESS;
862
  }
Tomas Dado's avatar
Tomas Dado committed
863
864
865

  StatusCode JetObjectCollectionMaker::tagLargeRJet(const xAOD::Jet& jet) {
    //decorate with boosted-tagging flags
866
867
868
869
870
871
872
873
874
875
876
877
    
    double lowMassCut{35000.};
    double highMassCut{1000000.};
    double lowPtCut{300000.};
    double highPtCut{3000000.};
    double etaCut{2.0};
    
    const double jet_pt=jet.pt();
    const double jet_m=jet.m();
    const double jet_eta=jet.eta();
    
    
Tomas Dado's avatar
Tomas Dado committed
878
    for (const std::pair<std::string, std::string>& name : m_config->boostedJetTaggers()) {
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
      const std::string fullName = name.first + "_" + name.second;

      // Checking kinematic range
      if(name.first=="SmoothedWZTagger") lowPtCut = 175000.;
      else lowPtCut=300000.;
      
      char passedRangeCheck = (jet_pt>lowPtCut &&
			       jet_pt<highPtCut && 
			       jet_m > lowMassCut &&
			       jet_m < highMassCut &&
			       std::abs(jet_eta)<etaCut);
      jet.auxdecor<char>("passedRangeCheck_" + fullName) = passedRangeCheck;
      if(!passedRangeCheck) continue;
      
      // Tagging jet
Tomas Dado's avatar
Tomas Dado committed
894
895
896
897
898
      const Root::TAccept& result = m_boostedJetTaggers[fullName]->tag(jet);
      // TAccept has bool operator overloaded, but let's be more explicit in the output to char
      jet.auxdecor<char>("isTagged_" + fullName) = (result ? 1 : 0);
      // for users to extract more detailed tagging result information in custom event saver
      jet.auxdecor<unsigned long>("TAccept_bitmap_" + fullName) = result.getCutResultBitSet().to_ulong();
899
    }
Tomas Dado's avatar
Tomas Dado committed
900
    return StatusCode::SUCCESS;
901
902
  }

Tomas Dado's avatar
Tomas Dado committed
903
904
905
906
907
908
909
910
911
912
  StatusCode JetObjectCollectionMaker::decorateBJets(xAOD::Jet& jet) {
    // initialise decorator and accessor
    static SG::AuxElement::Decorator<char> isbjet("IsBjet");
    static const std::string labelB = "PartonTruthLabelID";
    static SG::AuxElement::Accessor<int> truth_label(labelB);

    // Is b-jet if truth label == 5 and pT > 15 GeV
    isbjet(jet) = (jet.pt() > 15000. && truth_label(jet) == 5);

    return StatusCode::SUCCESS;
913
  }
914

Tomas Dado's avatar
Tomas Dado committed
915
916
917
918
919
920
921
922
923
924
  StatusCode JetObjectCollectionMaker::decorateHSJets() {
    // initialise decorator
    static SG::AuxElement::Decorator<char> isHS("AnalysisTop_isHS");

    // retrieve small-R jets collection
    const xAOD::JetContainer* jets(nullptr);

    top::check(evtStore()->retrieve(jets,
                                    m_config->sgKeyJets()),
               "Failed to retrieve small-R jet collection" + m_config->sgKeyJets());
925

Tomas Dado's avatar
Tomas Dado committed
926
927
928
929
930
931
932
933
934
935
936
    // retrieve truth jet collection
    const xAOD::JetContainer* truthJets = nullptr;
    top::check(evtStore()->retrieve(truthJets,
                                    m_truthJetCollForHS),
               "Unable to retrieve truth jet container " + m_truthJetCollForHS +
               " - this is needed to define HS jets for application of JVT");

    for (const auto& jet : *jets) {
      bool ishs = false;
      for (const auto& tjet : *truthJets) {
        if (tjet->p4().DeltaR(jet->p4()) < 0.3 && tjet->pt() > 10e3) ishs = true;
937
      }
Tomas Dado's avatar
Tomas Dado committed
938
      isHS(*jet) = ishs;
939
    }
Tomas Dado's avatar
Tomas Dado committed
940
941
942

    return StatusCode::SUCCESS;
  }
943
944
  
  StatusCode JetObjectCollectionMaker::decorateMatchedTruth() {
Tomas Dado's avatar
Tomas Dado committed
945
    static const SG::AuxElement::Decorator<float> matchedPt("AnalysisTop_MatchedTruthJetPt");
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
    // retrieve small-R jets collection
    const xAOD::JetContainer* jets(nullptr);

    top::check(evtStore()->retrieve(jets,
                                    m_config->sgKeyJets()),
               "Failed to retrieve small-R jet collection" + m_config->sgKeyJets());

    const xAOD::JetContainer* truthJets = nullptr;
    top::check(asg::AsgTool::evtStore()->retrieve(truthJets, m_config->sgKeyTruthJets()), "Failed to retrieve the truth jets");

    const xAOD::Jet* matchedTruthJet = nullptr;
    double deltaR(9999);
    
    for (const auto& jet : *jets) {
      // loop over truth jets
      for (const auto& iTruthJet : *truthJets) {
        TLorentzVector truthJetTLV;
        truthJetTLV.SetPtEtaPhiE(iTruthJet->pt(),iTruthJet->eta(),iTruthJet->phi(),iTruthJet->e());

        // do the matching
        if(!matchedTruthJet) {
          matchedTruthJet = iTruthJet;
        } else {
          const double newdR = jet->p4().DeltaR(iTruthJet->p4());
          if(newdR < deltaR) {
            deltaR = newdR;
            matchedTruthJet = iTruthJet;
          }
        }
      }
Tomas Dado's avatar
Tomas Dado committed
976
      if (deltaR > m_config->jetResponseMatchingDeltaR()) {
977
978
979
980
981
982
983
984
985
986
987
988
        matchedPt(*jet) = -9999;
        continue;
      }
      if (!matchedTruthJet) {
        matchedPt(*jet) = -9999;
        continue;
      }
      matchedPt(*jet) = matchedTruthJet->pt(); 
    }

    return StatusCode::SUCCESS;
  }
Tomas Dado's avatar
Tomas Dado committed
989
990
991
992
993
994
995
996
997
998
999
1000
1001

  StatusCode JetObjectCollectionMaker::decorateDL1(bool trackJets) {
    // retrieve small-R jets collection -- either calo or track jets
    const xAOD::JetContainer* jets(nullptr);

    if (trackJets) {
      top::check(evtStore()->retrieve(jets,
                                      m_config->sgKeyTrackJets()),
                 "Failed to retrieve track jet collection" + m_config->sgKeyTrackJets());
    } else {
      top::check(evtStore()->retrieve(jets,
                                      m_config->sgKeyJets()),
                 "Failed to retrieve small-R jet collection" + m_config->sgKeyJets());
1002
    }
Tomas Dado's avatar
Tomas Dado committed
1003
1004

    for (const auto& jet : *jets) {
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
      // loop over either calo or track jet btag selection tools to calculate the DL1x scores
      const std::unordered_map<std::string, ToolHandle<IBTaggingSelectionTool>>& m_btagDecorTools \
        = (trackJets ? m_btagSelToolsDL1Decor_trkJet : m_btagSelToolsDL1Decor);
      for (std::pair<std::string, ToolHandle<IBTaggingSelectionTool>> algo : m_btagDecorTools) {
        double DL1_weight = -999., dl1_pb = -10., dl1_pc = -10., dl1_pu = -10.;
        if (jet->btagging()->pb(algo.first, dl1_pb)
            && jet->btagging()->pc(algo.first, dl1_pc)
            && jet->btagging()->pu(algo.first, dl1_pu)) {
          if (!algo.second->getTaggerWeight(dl1_pb, dl1_pc, dl1_pu, DL1_weight)) {
            DL1_weight = -999.; // value for errors from retrieving DL1x weight
1015
          }
1016
1017
        } else {
          DL1_weight = -100.; // value for errors from nonexistence of probabilities
Tomas Dado's avatar
Tomas Dado committed
1018
        }
1019
        DLx.at(algo.first)(*jet) = DL1_weight;
Tomas Dado's avatar
Tomas Dado committed
1020
      }
1021
    }
Tomas Dado's avatar
Tomas Dado committed
1022
1023

    return StatusCode::SUCCESS;
1024
  }
1025
}  // namespace top