JetCalibrationTool.cxx 24.4 KB
Newer Older
1
2
3
///////////////////////// -*- C++ -*- /////////////////////////////

/*
4
  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
*/

// JetCalibrationTool.cxx 
// Implementation file for class JetCalibrationTool
// Author: Joe Taenzer <joseph.taenzer@cern.ch>
/////////////////////////////////////////////////////////////////// 


// JetCalibTools includes
#include "JetCalibTools/JetCalibrationTool.h"
#include "PathResolver/PathResolver.h"

// Constructors
////////////////

JetCalibrationTool::JetCalibrationTool(const std::string& name)
21
  : JetCalibrationToolBase::JetCalibrationToolBase( name ),
22
23
24
    m_rhkEvtInfo("EventInfo"),
    m_rhkRhoKey(""),
    m_rhkPV("PrimaryVertices"),
25
    m_jetAlgo(""), m_config(""), m_calibSeq(""), m_calibAreaTag(""), m_originScale(""), m_devMode(false), m_isData(true), m_timeDependentCalib(false), m_rhoKey("auto"), m_dir(""), m_eInfoName(""), m_globalConfig(NULL),
26
    m_doJetArea(true), m_doResidual(true), m_doOrigin(true), m_doGSC(true), m_gscDepth("auto"),
27
    m_jetPileupCorr(NULL), m_etaJESCorr(NULL), m_globalSequentialCorr(NULL), m_insituDataCorr(NULL), m_jetMassCorr(NULL), m_jetSmearCorr(NULL)
28
29
30
{ 

  declareProperty( "JetCollection", m_jetAlgo = "AntiKt4LCTopo" );
31
  declareProperty( "RhoKey", m_rhoKey = "auto" );
32
33
34
  declareProperty( "ConfigFile", m_config = "" );
  declareProperty( "CalibSequence", m_calibSeq = "JetArea_Offset_AbsoluteEtaJES_Insitu" );
  declareProperty( "IsData", m_isData = true );
35
  declareProperty( "ConfigDir", m_dir = "JetCalibTools/CalibrationConfigs/" );
36
  declareProperty( "EventInfoName", m_eInfoName = "EventInfo");
37
38
39
  declareProperty( "DEVmode", m_devMode = false);
  declareProperty( "OriginScale", m_originScale = "JetOriginConstitScaleMomentum");
  declareProperty( "CalibArea", m_calibAreaTag = "00-04-82");
40
  declareProperty( "rhkRhoKey", m_rhkRhoKey);
41
  declareProperty( "PrimaryVerticesContainerName", m_rhkPV = "PrimaryVertices");
42
  declareProperty( "GSCDepth", m_gscDepth);
43
44
45
46
47
48
49
50
51
52
53
54
55

}


// Destructor
///////////////
JetCalibrationTool::~JetCalibrationTool() {

  if (m_globalConfig) delete m_globalConfig;
  if (m_jetPileupCorr) delete m_jetPileupCorr;
  if (m_etaJESCorr) delete m_etaJESCorr;
  if (m_globalSequentialCorr) delete m_globalSequentialCorr;
  if (m_insituDataCorr) delete m_insituDataCorr;
56
  if (m_jetMassCorr) delete m_jetMassCorr;
57
  if (m_jetSmearCorr) delete m_jetSmearCorr;
58
59
60
61
62
63
64
65
66
67
68

}


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

StatusCode JetCalibrationTool::initialize() {
  ATH_MSG_INFO ("Initializing " << name() << "...");
  ATH_CHECK( this->initializeTool( name() ) );
69
70
71

  // Initialise ReadHandle(s)
  ATH_CHECK( m_rhkEvtInfo.initialize() );
Teng Jian Khoo's avatar
Teng Jian Khoo committed
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
  if(m_rhkPV.empty()) {
    // No PV key: -- check if it is required
    if(m_doResidual) {
      // May require modification in case of residual that does not require NPV
      ATH_MSG_ERROR("Residual calibration requested but no primary vertex container specified!");
      return StatusCode::FAILURE;
    } else if(m_doGSC) {
      if(m_jetAlgo.find("PFlow")!=std::string::npos) {
	ATH_MSG_ERROR("GSC calibration for PFlow requested but no primary vertex container specified!");
	return StatusCode::FAILURE;
      } else if((m_gscDepth!="Tile0" && m_gscDepth!="EM3")) {
	ATH_MSG_ERROR("GSC calibration with tracks requested but no primary vertex container specified!");
	return StatusCode::FAILURE;
      }
    }
  } else {
    // Received a PV key, declare the data dependency
    ATH_CHECK( m_rhkPV.initialize() );
  }
91

92
93
94
95
96
97
98
99
100
101
102
103
104
105
  return StatusCode::SUCCESS;
}

StatusCode JetCalibrationTool::finalize() {
  ATH_MSG_INFO ("Finalizing " << name() << "...");
  return StatusCode::SUCCESS;
}


StatusCode JetCalibrationTool::initializeTool(const std::string& name) {

  TString jetAlgo = m_jetAlgo;
  TString config = m_config;
  TString calibSeq = m_calibSeq;
106
  std::string dir = m_dir;
107

108
109
  ATH_MSG_INFO("===================================");
  ATH_MSG_INFO("Initializing the xAOD Jet Calibration Tool for " << jetAlgo << "jets");
110
111
112
113
114
115
116

  //Make sure the necessary properties were set via the constructor or python configuration
  if ( jetAlgo.EqualTo("") || calibSeq.EqualTo("") ) {
    ATH_MSG_FATAL("JetCalibrationTool::initializeTool : At least one of your constructor arguments is not set. Did you use the copy constructor?");
    return StatusCode::FAILURE;
  }

117
  if ( config.EqualTo("") || !config ) { ATH_MSG_FATAL("No configuration file specified."); return StatusCode::FAILURE; } 
118
119
120
121
122
123
  // The calibration area tag is a property of the tool
  const std::string calibPath = "CalibArea-" + m_calibAreaTag + "/";
  if(m_devMode){
    ATH_MSG_WARNING("Dev Mode is ON!!!");
    ATH_MSG_WARNING("Dev Mode is NOT RECOMMENDED!!!");
    dir = "JetCalibTools/";
124
  }
125
  else{dir.insert(14,calibPath);} // Obtaining the path of the configuration file
126
  std::string configPath=dir+m_config; // Full path
127
  TString fn =  PathResolverFindCalibFile(configPath);
128

129
  ATH_MSG_INFO("Reading global JES settings from: " << m_config);
130
  ATH_MSG_INFO("resolved in: " << fn);
131
132
133
134
135
136
137
  
  m_globalConfig = new TEnv();
  //int status=m_globalConfig->ReadFile(FindFile(fn),EEnvLevel(0));
  int status=m_globalConfig->ReadFile(fn ,EEnvLevel(0));
  if (status!=0) { ATH_MSG_FATAL("Cannot read config file " << fn ); return StatusCode::FAILURE; }

  //Make sure that one of the standard jet collections is being used
138
139
140
141
142
143
  if ( calibSeq.Contains("JetArea") ) {
    if ( jetAlgo.Contains("PFlow") ) m_jetScale = PFLOW;
    else if ( jetAlgo.Contains("EM") ) m_jetScale = EM;
    else if ( jetAlgo.Contains("LC") ) m_jetScale = LC;
    else { ATH_MSG_FATAL("jetAlgo " << jetAlgo << " not recognized."); return StatusCode::FAILURE; }
  }
144
145
146
147

  //Set the default units to MeV, user can override by calling setUnitsGeV(true)
  setUnitsGeV(false);

148
  // Settings for R21/2.5.X
149
  m_originCorrectedClusters = m_globalConfig->GetValue("OriginCorrectedClusters",false);
150
  m_doSetDetectorEta = m_globalConfig->GetValue("SetDetectorEta",true);
151

152
153
154
155
  //Make sure the residual correction is turned on if requested, protect against applying it without the jet area subtraction                    
  if ( !calibSeq.Contains("JetArea") && !calibSeq.Contains("Residual") ) {
    m_doJetArea = false;
    m_doResidual = false;
156
  } else if ( calibSeq.Contains("JetArea") ) {
157
    if ( m_rhoKey.compare("auto") == 0 ) {
158
159
160
161
162
163
164
165
166
      if(!m_originCorrectedClusters){
        if ( m_jetScale == EM ) m_rhoKey = "Kt4EMTopoEventShape";
        else if ( m_jetScale == LC ) m_rhoKey = "Kt4LCTopoEventShape";
        else if ( m_jetScale == PFLOW ) m_rhoKey = "Kt4EMPFlowEventShape";
      } else{
        if ( m_jetScale == EM ) m_rhoKey = "Kt4EMTopoOriginEventShape";
        else if ( m_jetScale == LC ) m_rhoKey = "Kt4LCTopoOriginEventShape";
        else if ( m_jetScale == PFLOW ) m_rhoKey = "Kt4EMPFlowEventShape";
      }
167
    }
168
169
    ATH_CHECK( m_rhkRhoKey.assign(m_rhoKey)) ;  // set in `initializeTool'
    ATH_CHECK( m_rhkRhoKey.initialize() );
170
    if ( !calibSeq.Contains("Residual") ) m_doResidual = false;
171
  } else if ( !calibSeq.Contains("JetArea") && calibSeq.Contains("Residual") ) {
172
    ATH_MSG_FATAL("JetCalibrationTool::initializeTool : You are trying to initialize JetCalibTools with the jet area correction turned off and the residual offset correction turned on. This is inconsistent. Aborting.");
173
    return StatusCode::SUCCESS;
174
175
  }

176
177
  if ( !calibSeq.Contains("Origin") ) m_doOrigin = false;

178
179
  if ( !calibSeq.Contains("GSC") ) m_doGSC = false;

180
181
182
183
184
185
  //Protect against the in-situ calibration being requested when isData is false
  if ( calibSeq.Contains("Insitu") && !m_isData ) {
    ATH_MSG_FATAL("JetCalibrationTool::initializeTool : calibSeq string contains Insitu with isData set to false. Can't apply in-situ correction to MC!!");
    return StatusCode::FAILURE;
  }

186
187
  // Time-Dependent Insitu Calibration
  m_timeDependentCalib = m_globalConfig->GetValue("TimeDependentInsituCalibration",false);
188
  if(m_timeDependentCalib && calibSeq.Contains("Insitu")){ // Read Insitu Configs
189
190
191
192
193
194
195
196
197
198
199
200
201
202
    m_timeDependentInsituConfigs = JetCalibUtils::Vectorize( m_globalConfig->GetValue("InsituTimeDependentConfigs","") );
    if(m_timeDependentInsituConfigs.size()==0) ATH_MSG_ERROR("Please check there are at least two insitu configs");
    m_runBins = JetCalibUtils::VectorizeD( m_globalConfig->GetValue("InsituRunBins","") );
    if(m_runBins.size()!=m_timeDependentInsituConfigs.size()+1) ATH_MSG_ERROR("Please check the insitu run bins");
    for(unsigned int i=0;i<m_timeDependentInsituConfigs.size();++i){
      //InsituDataCorrection *insituTemp = NULL;
      //m_insituTimeDependentCorr.push_back(insituTemp);

      std::string configPath_insitu = dir+m_timeDependentInsituConfigs.at(i).Data(); // Full path
      TString fn_insitu =  PathResolverFindCalibFile(configPath_insitu);

      ATH_MSG_INFO("Reading time-dependent insitu settings from: " << m_timeDependentInsituConfigs.at(i));
      ATH_MSG_INFO("resolved in: " << fn_insitu);
  
203
204
      TEnv *globalConfig_insitu = new TEnv();
      int status = globalConfig_insitu->ReadFile(fn_insitu ,EEnvLevel(0));
205
      if (status!=0) { ATH_MSG_FATAL("Cannot read config file " << fn_insitu ); return StatusCode::FAILURE; }
206
      m_globalTimeDependentConfigs.push_back(globalConfig_insitu);
207
208
209
    }
  }

210
211
212
  //Loop over the request calib sequence
  //Initialize derived classes for applying the requested calibrations and add them to a vector
  std::vector<TString> vecCalibSeq = JetCalibUtils::Vectorize(calibSeq,"_");
213
  TString vecCalibSeqtmp;
214
  for ( unsigned int i=0; i<vecCalibSeq.size(); ++i) {
215
    if ( vecCalibSeq[i].EqualTo("Residual") || vecCalibSeq[i].EqualTo("Origin") || vecCalibSeq[i].EqualTo("DEV") ) continue;
216
    ATH_CHECK( getCalibClass(name,vecCalibSeq[i] ));
217
218
  }

219
  ATH_MSG_INFO("===================================");
220
221
222
223
224
225
226

  return StatusCode::SUCCESS;
}

//Method for initializing the requested calibration derived classes
StatusCode JetCalibrationTool::getCalibClass(const std::string&name, TString calibration) {
  TString jetAlgo = m_jetAlgo;
227
  const TString calibPath = "CalibArea-" + m_calibAreaTag + "/";
228
  std::string suffix = "";
229
  //ATH_MSG_INFO("Initializing sub tools.");
230
  if ( calibration.EqualTo("JetArea") ) {
231
    ATH_MSG_INFO("Initializing pileup correction.");
232
    suffix="_Pileup";
233
    if(m_devMode) suffix+="_DEV";
234
235
    m_jetPileupCorr = new JetPileupCorrection(name+suffix,m_globalConfig,jetAlgo,calibPath,m_doResidual,m_doOrigin,m_isData,m_devMode);
    ATH_CHECK( m_jetPileupCorr->setProperty("OriginScale",m_originScale.c_str()) );
236
    m_jetPileupCorr->msg().setLevel( this->msg().level() );
237
238
239
240
241
242
243
    if( m_jetPileupCorr->initializeTool(name+suffix).isFailure() ) { 
      ATH_MSG_FATAL("Couldn't initialize the pileup correction. Aborting"); 
      return StatusCode::FAILURE; 
    } else { 
      m_calibClasses.push_back(m_jetPileupCorr); 
      return StatusCode::SUCCESS; 
    }
244
  } else if ( calibration.EqualTo("EtaJES") || calibration.EqualTo("AbsoluteEtaJES") ) {
245
    ATH_MSG_INFO("Initializing JES correction.");
246
    suffix="_EtaJES";
247
    if(m_devMode) suffix+="_DEV";
248
    m_etaJESCorr = new EtaJESCorrection(name+suffix,m_globalConfig,jetAlgo,calibPath,false,m_devMode);
249
    m_etaJESCorr->msg().setLevel( this->msg().level() );
250
251
252
253
254
255
256
    if ( m_etaJESCorr->initializeTool(name+suffix).isFailure() ) {
      ATH_MSG_FATAL("Couldn't initialize the Monte Carlo JES correction. Aborting"); 
      return StatusCode::FAILURE; 
    } else { 
      m_calibClasses.push_back(m_etaJESCorr); 
      return StatusCode::SUCCESS; 
    }
257
  } else if ( calibration.EqualTo("EtaMassJES") ) {
258
    ATH_MSG_INFO("Initializing JES correction.");
259
    suffix="_EtaMassJES";
260
    if(m_devMode) suffix+="_DEV";
261
    m_etaJESCorr = new EtaJESCorrection(name+suffix,m_globalConfig,jetAlgo,calibPath,true,m_devMode);
262
    m_etaJESCorr->msg().setLevel( this->msg().level() );
263
264
265
266
267
268
269
    if ( m_etaJESCorr->initializeTool(name+suffix).isFailure() ) {
      ATH_MSG_FATAL("Couldn't initialize the Monte Carlo JES correction. Aborting");
      return StatusCode::FAILURE;
    } else {
      m_calibClasses.push_back(m_etaJESCorr);
      return StatusCode::SUCCESS;
    }
270
  } else if ( calibration.EqualTo("GSC") ) {
271
    ATH_MSG_INFO("Initializing GSC correction.");
272
    suffix="_GSC";
273
    if(m_devMode) suffix+="_DEV";
274
    m_globalSequentialCorr = new GlobalSequentialCorrection(name+suffix,m_globalConfig,jetAlgo,m_gscDepth,calibPath,m_devMode);
275
    m_globalSequentialCorr->msg().setLevel( this->msg().level() );
276
277
278
279
280
281
282
    if ( m_globalSequentialCorr->initializeTool(name+suffix).isFailure() ) {
      ATH_MSG_FATAL("Couldn't initialize the Global Sequential Calibration. Aborting"); 
      return StatusCode::FAILURE; 
    } else { 
      m_calibClasses.push_back(m_globalSequentialCorr); 
      return StatusCode::SUCCESS; 
    }
283
  } else if ( calibration.EqualTo("JMS") ) {
284
    ATH_MSG_INFO("Initializing JMS correction.");
285
    suffix="_JMS";
286
    if(m_devMode) suffix+="_DEV";
287
    m_jetMassCorr = new JMSCorrection(name+suffix,m_globalConfig,jetAlgo,calibPath,m_devMode);
288
    m_jetMassCorr->msg().setLevel( this->msg().level() );
289
290
291
292
293
294
295
    if ( m_jetMassCorr->initializeTool(name+suffix).isFailure() ) {
      ATH_MSG_FATAL("Couldn't initialize the JMS Calibration. Aborting");
      return StatusCode::FAILURE;
    } else {
      m_calibClasses.push_back(m_jetMassCorr);
      return StatusCode::SUCCESS;
    }
296
  } else if ( calibration.EqualTo("Insitu") ) {
297
298
299
300
    if(!m_timeDependentCalib){
      ATH_MSG_INFO("Initializing Insitu correction.");
      suffix="_Insitu";
      if(m_devMode) suffix+="_DEV";
301
      m_insituDataCorr = new InsituDataCorrection(name+suffix,m_globalConfig,jetAlgo,calibPath,m_devMode);
302
303
304
305
306
      m_insituDataCorr->msg().setLevel( this->msg().level() );
      if ( m_insituDataCorr->initializeTool(name+suffix).isFailure() ) {
        ATH_MSG_FATAL("Couldn't initialize the In-situ data correction. Aborting"); 
        return StatusCode::FAILURE; 
      } else { 
307
308
309
        m_calibClasses.push_back(m_insituDataCorr);
	m_relInsituPtMax.push_back(m_insituDataCorr->getRelHistoPtMax());
	m_absInsituPtMax.push_back(m_insituDataCorr->getAbsHistoPtMax());
310
311
312
313
314
315
316
        return StatusCode::SUCCESS; 
      }
    } else{
      ATH_MSG_INFO("Initializing Time-Dependent Insitu Corrections");
      for(unsigned int i=0;i<m_timeDependentInsituConfigs.size();++i){
        suffix="_Insitu"; suffix += "_"; suffix += std::to_string(i);
        if(m_devMode) suffix+="_DEV";
317
        InsituDataCorrection *insituTimeDependentCorr_Tmp = new InsituDataCorrection(name+suffix,m_globalTimeDependentConfigs.at(i),jetAlgo,calibPath,m_devMode);
318
319
        insituTimeDependentCorr_Tmp->msg().setLevel( this->msg().level() );
        if ( insituTimeDependentCorr_Tmp->initializeTool(name+suffix).isFailure() ) {
320
321
322
          ATH_MSG_FATAL("Couldn't initialize the In-situ data correction. Aborting"); 
          return StatusCode::FAILURE; 
        } else {     		
323
324
325
          m_insituTimeDependentCorr.push_back(insituTimeDependentCorr_Tmp);
	  m_relInsituPtMax.push_back(insituTimeDependentCorr_Tmp->getRelHistoPtMax());
	  m_absInsituPtMax.push_back(insituTimeDependentCorr_Tmp->getAbsHistoPtMax());
326
327
        }
      }
328
329
      return StatusCode::SUCCESS; 
    }
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
  } else if ( calibration.EqualTo("Smear") ) {
    if (m_isData)
    {
      ATH_MSG_FATAL("Asked for smearing of data, which is not supported.  Aborting.");
      return StatusCode::FAILURE;
    }

    ATH_MSG_INFO("Initializing jet smearing correction");
    suffix = "_Smear";
    if (m_devMode) suffix += "_DEV";
    
    m_jetSmearCorr = new JetSmearingCorrection(name+suffix,m_globalConfig,jetAlgo,calibPath,m_devMode);
    m_jetSmearCorr->msg().setLevel(this->msg().level());

    if (m_jetSmearCorr->initializeTool(name+suffix).isFailure())
    {
      ATH_MSG_FATAL("Couldn't initialize the jet smearing correction.  Aborting.");
      return StatusCode::FAILURE;
    }
    m_calibClasses.push_back(m_jetSmearCorr);
    return StatusCode::SUCCESS;
351
352
  }
  ATH_MSG_FATAL("Calibration string not recognized: " << calibration << ", aborting.");
353
  //ATH_MSG_INFO("Initializing of sub tools is complete.");
354
355
356
  return StatusCode::FAILURE;
}

357
StatusCode JetCalibrationTool::applyCalibration(xAOD::JetContainer& jets) const {
358
359
360
  //Grab necessary event info for pile up correction and store it in a JetEventInfo class object
  ATH_MSG_VERBOSE("Modifying jet collection.");
  JetEventInfo jetEventInfo;
361
  ATH_CHECK( initializeEvent(jetEventInfo) );
362
363
364
  xAOD::JetContainer::iterator jet_itr = jets.begin();
  xAOD::JetContainer::iterator jet_end = jets.end(); 
  for ( ; jet_itr != jet_end; ++jet_itr )
365
366
    ATH_CHECK( calibrateImpl(**jet_itr,jetEventInfo) );
 return StatusCode::SUCCESS;
367
368
369
370
371
372
}

// Private/Protected Methods
///////////////

StatusCode JetCalibrationTool::initializeEvent(JetEventInfo& jetEventInfo) const {
373
374
375
376
377
378
379

  // Check if the tool was initialized
  if( m_calibClasses.size() == 0 ){
    ATH_MSG_FATAL("   JetCalibrationTool::initializeEvent : The tool was not initialized.");
    return StatusCode::FAILURE;
  }

380
381
382
  // static accessor for PV index access
  static SG::AuxElement::ConstAccessor<int> PVIndexAccessor("PVIndex");
  
383
384
  ATH_MSG_VERBOSE("Initializing event.");

385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
  if( m_doJetArea ) {
    //Determine the rho value to use for the jet area subtraction
    //Should be determined using EventShape object, use hard coded values if EventShape doesn't exist
    double rho=0;
    const xAOD::EventShape * eventShape = 0;

    SG::ReadHandle<xAOD::EventShape> rhRhoKey(m_rhkRhoKey);

    if ( rhRhoKey.isValid() ) {
      ATH_MSG_VERBOSE("  Found event density container " << m_rhoKey);
      eventShape = rhRhoKey.cptr();
      if ( !rhRhoKey.isValid() ) {
	ATH_MSG_VERBOSE("  Event shape container not found.");
	ATH_MSG_FATAL("Could not retrieve xAOD::EventShape DataHandle.");
	return StatusCode::FAILURE;
      } else if ( !eventShape->getDensity( xAOD::EventShape::Density, rho ) ) {
	ATH_MSG_VERBOSE("  Event density not found in container.");
	ATH_MSG_FATAL("Could not retrieve xAOD::EventShape::Density from xAOD::EventShape.");
	return StatusCode::FAILURE;
      } else {
	ATH_MSG_VERBOSE("  Event density retrieved.");
      }
    } else if ( m_doJetArea && !rhRhoKey.isValid() ) {
      ATH_MSG_VERBOSE("  Rho container not found: " << m_rhoKey);
409
      ATH_MSG_FATAL("Could not retrieve xAOD::EventShape DataHandle.");
410
      return StatusCode::FAILURE;
411
    }
412
413
    jetEventInfo.setRho(rho);
    ATH_MSG_VERBOSE("  Rho = " << 0.001*rho << " GeV");
414
415
  }

416
  // Retrieve EventInfo object, which now has multiple uses
417
  if ( m_doResidual || m_doGSC ) {
418
419
    const xAOD::EventInfo * eventObj = 0;
    static unsigned int eventInfoWarnings = 0;
420
421
422
423
    SG::ReadHandle<xAOD::EventInfo> rhEvtInfo(m_rhkEvtInfo);
    if ( rhEvtInfo.isValid() ) {
      eventObj = rhEvtInfo.cptr();
    } else {
424
425
      ++eventInfoWarnings;
      if ( eventInfoWarnings < 20 )
426
        ATH_MSG_ERROR("   JetCalibrationTool::initializeEvent : Failed to retrieve event information.");
427
      jetEventInfo.setMu(0); //Hard coded value mu = 0 in case of failure (to prevent seg faults later).
428
      jetEventInfo.setPVIndex(0);
429
      return StatusCode::SUCCESS; //error is recoverable, so return SUCCESS
430
431
    }

432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
    // If we are applying the reisdual, then store mu
    if (m_doResidual)
      jetEventInfo.setMu( eventObj->averageInteractionsPerCrossing() );
    
    // If this is GSC, we need EventInfo to determine the PV to use
    // This is support for groups where PV0 is not the vertex of interest (H->gamgam, etc)
    if (m_doGSC)
    {
      // First retrieve the PVIndex if specified
      // Default is to not specify this, so no warning if it doesn't exist
      // However, if specified, it should be a sane value - fail if not
      if ( m_doGSC && PVIndexAccessor.isAvailable(*eventObj) )
        jetEventInfo.setPVIndex( PVIndexAccessor(*eventObj) );
      else
        jetEventInfo.setPVIndex(0);
      
    }

    // If PV index is not zero, we need to confirm it's a reasonable value
    // To do this, we need the primary vertices
    // However, other users of the GSC may not have the PV collection (in particular: trigger GSC in 2016)
    // So only retrieve vertices if needed for NPV (residual) or a non-zero PV index was specified (GSC)
    if (m_doResidual || (m_doGSC && jetEventInfo.PVIndex()))
    {
      //Retrieve VertexContainer object, use it to obtain NPV for the residual correction or check validity of GSC non-PV0 usage
      const xAOD::VertexContainer * vertices = 0;
      static unsigned int vertexContainerWarnings = 0;
459
460
461
462
463

      SG::ReadHandle<xAOD::VertexContainer> rhPV(m_rhkPV);
      if (rhPV.isValid()) {
        vertices = rhPV.cptr(); 
      } else {
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
        ++vertexContainerWarnings;
        if ( vertexContainerWarnings < 20 )
          ATH_MSG_ERROR("   JetCalibrationTool::initializeEvent : Failed to retrieve primary vertices.");
        jetEventInfo.setNPV(0); //Hard coded value NPV = 0 in case of failure (to prevent seg faults later).
        return StatusCode::SUCCESS; //error is recoverable, so return SUCCESS
      }

      // Calculate and set NPV if this is residual
      if (m_doResidual)
      {
        int eventNPV = 0;
        xAOD::VertexContainer::const_iterator vtx_itr = vertices->begin();
        xAOD::VertexContainer::const_iterator vtx_end = vertices->end(); 
        for ( ; vtx_itr != vtx_end; ++vtx_itr ) 
          if ( (*vtx_itr)->nTrackParticles() >= 2 ) ++eventNPV;
  
        jetEventInfo.setNPV(eventNPV);
      }
      
      // Validate value of non-standard PV index usage
      if (m_doGSC && jetEventInfo.PVIndex())
      {
        static unsigned int vertexIndexWarnings = 0;
        if (jetEventInfo.PVIndex() < 0 || static_cast<size_t>(jetEventInfo.PVIndex()) >= vertices->size())
        {
          ++vertexIndexWarnings;
          if (vertexIndexWarnings < 20)
            ATH_MSG_ERROR("   JetCalibrationTool::initializeEvent : PV index is out of bounds.");
          jetEventInfo.setPVIndex(0); // Hard coded value PVIndex = 0 in case of failure (to prevent seg faults later).
          return StatusCode::SUCCESS; // error is recoverable, so return SUCCESS
        }
      }
    }
497

498
499
500
501
502
503
504
505
506
507
508
    //Check if the input jets are coming from data or MC
    //if ( m_eventObj->eventType( xAOD::EventInfo::IS_SIMULATION ) ) {
    //     m_isData = false; // controls mu scaling in the pile up correction, no scaling for data
    //}

  }

    return StatusCode::SUCCESS;
}

StatusCode JetCalibrationTool::calibrateImpl(xAOD::Jet& jet, JetEventInfo& jetEventInfo) const {
509
510
511
512
513
514
515
516

  //Check for OriginCorrected and PileupCorrected attributes, assume they are false if not found
  int tmp = 0; //temporary int for checking getAttribute
  if ( !jet.getAttribute<int>("OriginCorrected",tmp) )
    jet.setAttribute<int>("OriginCorrected",false);
  if ( !jet.getAttribute<int>("PileupCorrected",tmp) )
    jet.setAttribute<int>("PileupCorrected",false);

517
  ATH_MSG_VERBOSE("Calibrating jet " << jet.index());
518
519
520
521
  if(m_doSetDetectorEta) {
    xAOD::JetFourMom_t jetconstitP4 = jet.getAttribute<xAOD::JetFourMom_t>("JetConstitScaleMomentum");
    jet.setAttribute<float>("DetectorEta",jetconstitP4.eta()); //saving constituent scale eta for later use
  }
522
523
  for (unsigned int i=0; i<m_calibClasses.size(); ++i) //Loop over requested calibations
    ATH_CHECK ( m_calibClasses[i]->calibrateImpl(jet,jetEventInfo) );
524
525
526
  TString CalibSeq = m_calibSeq;
  if(CalibSeq.Contains("Insitu") && m_timeDependentCalib){ // Insitu Time-Dependent Correction
    for(unsigned int i=0;i<m_timeDependentInsituConfigs.size();++i){
527
      // Retrieve EventInfo container
528
      const xAOD::EventInfo* eventInfo(nullptr);
529
530
531
532
533
      SG::ReadHandle<xAOD::EventInfo> rhEvtInfo(m_rhkEvtInfo);
      if ( rhEvtInfo.isValid() ) {
        eventInfo = rhEvtInfo.cptr();
      } else {
        ATH_MSG_ERROR("   JetCalibrationTool::calibrateImpl : Failed to retrieve EventInfo.");
534
535
536
537
538
539
      }
      // Run Number Dependent Correction
      double runNumber = eventInfo->runNumber();
      if(runNumber>m_runBins.at(i) && runNumber<=m_runBins.at(i+1)){ ATH_CHECK ( m_insituTimeDependentCorr.at(i)->calibrateImpl(jet,jetEventInfo) );}
    }
  }
540
541
  return StatusCode::SUCCESS; 
}
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563


StatusCode JetCalibrationTool::getNominalResolutionData(const xAOD::Jet& jet, double& resolution) const
{
    if (!m_jetSmearCorr)
    {
        ATH_MSG_ERROR("Cannot retrieve the nominal data resolution - smearing was not configured during initialization");
        return StatusCode::FAILURE;
    }
    return m_jetSmearCorr->getNominalResolutionData(jet,resolution);
}

StatusCode JetCalibrationTool::getNominalResolutionMC(const xAOD::Jet& jet, double& resolution) const
{
    if (!m_jetSmearCorr)
    {
        ATH_MSG_ERROR("Cannot retrieve the nominal MC resolution - smearing was not configured during initialization");
        return StatusCode::FAILURE;
    }
    return m_jetSmearCorr->getNominalResolutionMC(jet,resolution);
}