ATLAS Offline Software
Loading...
Searching...
No Matches
METSignificance.cxx
Go to the documentation of this file.
1
2/*
3 Copyright (C) 2002-2026 CERN for the benefit of the ATLAS collaboration
4*/
5// METSignificance.cxx
6// Implementation file for class METSignificance
7// Author: P.Francavilla<francav@cern.ch>
8// Author: D.Schaefer<schae@cern.ch>
10
11// METUtilities includes
13
14// MET EDM
18
19// Jet EDM
21
22// Other xAOD EDM
28
32
33// Needed for xAOD::get_eta_calo() function
34
36
37#ifdef XAOD_STANDALONE
39#endif
40
41namespace met {
42
43 using iplink_t = ElementLink<xAOD::IParticleContainer>;
44
45 static const SG::ConstAccessor<float> acc_varX("varX");
46 static const SG::ConstAccessor<float> acc_varY("varY");
47 static const SG::ConstAccessor<float> acc_covXY("covXY");
48 static const SG::ConstAccessor<float> acc_jvt("Jvt");
49 static const SG::ConstAccessor<float> acc_fjvt("fJvt");
50 static const SG::ConstAccessor<float> acc_fjvt_der("DFCommonJets_fJvt");
52 const static MissingETBase::Types::bitmask_t invisSource = 0x100000; // doesn't overlap with any other
53
56 m_GeV(1.0e3),
57 m_jetOK(true),
58 m_muonOK(true),
59 m_egammaOK(true),
60 m_tauOK(true),
62 m_jerForEMu(false),
63 m_jetPtThr(-1.0),
64 m_jetEtaThr(-1.0),
65 m_significance(0.0),
66 m_rho(0.0),
67 m_VarL(0.0),
68 m_VarT(0.0),
69 m_CvLT(0.0),
70 m_met_VarL(0.0),
71 m_met_VarT(0.0),
72 m_met_CvLT(0.0),
73 m_met(0.0),
74 m_metx(0.0),
75 m_mety(0.0),
76 m_metphi(0.0),
77 m_metsoft(0.0),
78 m_metsoftphi(0.0),
79 m_ht(0.0),
80 m_sumet(0.0),
81 m_file(nullptr),
82 m_phi_reso_pt20(nullptr),
83 m_phi_reso_pt50(nullptr),
84 m_phi_reso_pt100(nullptr)
85 {
86 declareProperty("SoftTermParam", m_softTermParam = met::Random );
87 declareProperty("SoftTermReso", m_softTermReso = 8.5 );
88 declareProperty("TreatPUJets", m_treatPUJets = true );
89 declareProperty("DoPhiReso", m_doPhiReso = false );
90 declareProperty("ApplyBias", m_applyBias = false );
91 declareProperty("DoJerForEMu", m_jerForEMu = false ); // run jet resolution for all electrons and muons
92 declareProperty("ScalarBias", m_scalarBias = 0.0 );
93 declareProperty("JetPtThr", m_jetPtThr = -1.0 );
94 declareProperty("JetEtaThr", m_jetEtaThr = -1.0 );
95 declareProperty("ConfigPrefix", m_configPrefix = "METUtilities/data17_13TeV/metsig_Aug15/");
96 declareProperty("ConfigJetPhiResoFile", m_configJetPhiResoFile = "jet_unc.root" );
97 declareProperty("JetResoAux", m_JetResoAux = "" ); // relative pT resolution in addition to normal JES
98 declareProperty("EMuResoAux", m_EMuResoAux = "" ); // aux string sets a bool for the leptons to run the jet resolation
99 declareProperty("JetCollection", m_JetCollection = "AntiKt4EMPFlow" );
100
101 // properties to delete eventually
102 declareProperty("IsDataJet", m_isDataJet = false );
103 declareProperty("IsDataMuon", m_isDataMuon = false );
104
105 m_file = nullptr;
106 }
107
109
111
112 ATH_MSG_INFO ("Initializing " << name() << "...");
113 ATH_MSG_INFO("Set up JER tools");
114 if(m_JetCollection == "AntiKt4EMTopoJets"){
115 ATH_MSG_WARNING(" tool wasn't updated for EMTopo jets so far and is not supported.");
116 }
117 // Phi resolution
119 m_file = TFile::Open(configpath.c_str());
120 if(m_file){
121 m_phi_reso_pt20 = static_cast<TH2F *>(m_file->Get("phi_reso_pt20"));
122 if(!m_phi_reso_pt20) ATH_MSG_ERROR("PU Jet Uncertainty Histogram not valid");
123 m_phi_reso_pt50 = static_cast<TH2F *>(m_file->Get("phi_reso_pt50"));
124 if(!m_phi_reso_pt50) ATH_MSG_ERROR("PU Jet Uncertainty Histogram not valid");
125 m_phi_reso_pt100 = static_cast<TH2F *>(m_file->Get("phi_reso_pt100"));
126 if(!m_phi_reso_pt100) ATH_MSG_ERROR("PU Jet Uncertainty Histogram not valid");
127 }
128 else{
129 ATH_MSG_ERROR("PU Jet Uncertainty TFile is not valid: " << configpath);
130 return StatusCode::FAILURE;
131 }
132
133 if (m_jetCalibTool.empty()){
134 if(m_jetCalibConfig.empty() || m_jetCalibSeq.empty())
135 m_jetOK = false;
136 else{
137 asg::AsgToolConfig toolConfig ("JetCalibrationTool/MetSigAutoConf_JetCalibTool_"+m_JetCollection);
138 ATH_CHECK( toolConfig.setProperty("JetCollection",m_JetCollection) );
139 ATH_CHECK( toolConfig.setProperty("ConfigFile",m_jetCalibConfig) );
140 ATH_CHECK( toolConfig.setProperty("CalibSequence",m_jetCalibSeq) );
141 if(!m_jetCalibArea.empty()) ATH_CHECK( toolConfig.setProperty("CalibArea",m_jetCalibArea) );
142 ATH_CHECK( toolConfig.setProperty("IsData",false) ); // configure for MC due to technical reasons. Both data and MC smearing are available with this setting.
144 }
145 }
146 if(m_jetOK) ATH_CHECK( m_jetCalibTool.retrieve() );
147 else ATH_MSG_WARNING("No jet calibration tool or config provided for MET Significance");
148
150 if(m_muonCalibMode == -1)
151 m_muonOK = false;
152 else{
153 asg::AsgToolConfig toolConfig ("CP::MuonCalibTool/METSigAutoConf_MuonCalibrationAndSmearingTool");
154 ATH_CHECK(toolConfig.setProperty("calibMode", m_muonCalibMode));
156 }
157 }
159 else ATH_MSG_WARNING("No muon calibration tool or config provided for MET Significance");
160
161 if (m_egammaCalibTool.empty()){
162 if(m_egESModel.empty() || m_egDecorrModel.empty())
163 m_egammaOK = false;
164 else{
165 asg::AsgToolConfig toolConfig ("CP::EgammaCalibrationAndSmearingTool/METSigAutoConf_EgammaCalibrationAndSmearingTool");
166 ATH_CHECK(toolConfig.setProperty("ESModel", m_egESModel));
167 ATH_CHECK(toolConfig.setProperty("decorrelationModel", m_egDecorrModel));
168 ATH_CHECK(toolConfig.setProperty("useFastSim", m_egUseFastsim ? 1 : 0));
169 // only resolution() is used, which does not need the MVA BDTs
170 ATH_CHECK(toolConfig.setProperty("useMVACalibration", 0));
172 }
173 }
174 if(m_egammaOK) ATH_CHECK( m_egammaCalibTool.retrieve() );
175 else ATH_MSG_WARNING("No egamma calibration tool or config provided for MET Significance");
176
177 if (m_tauResolutionTool.empty()){
178 if(m_tauTESConfig.empty())
179 m_tauOK = false;
180 else{
181 asg::AsgToolConfig toolConfig ("TauAnalysisTools::TauResolutionTool/METSigAutoConf_TauPerfTool");
182 ATH_CHECK( toolConfig.setProperty("WeightFileName", m_tauTESConfig) );
184 }
185 }
186 if(m_tauOK) ATH_CHECK( m_tauResolutionTool.retrieve() );
187 else ATH_MSG_WARNING("No tau calibration tool or config provided for MET Significance");
188
189
190 return StatusCode::SUCCESS;
191 }
192
194
195 ATH_MSG_INFO ("Finalizing " << name() << "...");
196 delete m_phi_reso_pt20;
197 delete m_phi_reso_pt50;
198 delete m_phi_reso_pt100;
199
200 return StatusCode::SUCCESS;
201 }
202
203 // **** Rebuild generic MET term ****
204 StatusCode METSignificance::varianceMET(xAOD::MissingETContainer* metCont, float avgmu, const std::string& jetTermName, const std::string& softTermName, const std::string& totalMETName){
205
206 // reset variables
207 m_VarL = 0.0;
208 m_VarT = 0.0;
209 m_CvLT = 0.0;
210
211 int metTerm = 0;
212 double particle_sum[2][2] = {{0.0,0.0}, {0.0,0.0}};
213 m_metphi = 0.0; //Angle for rotation of the cov matrix
214 m_met = -1.0; // Numerator
215 m_metsoft = 0.0;
216 m_metsoftphi = 0.0;
217 m_sumet=-1.0;
218 m_ht=0.0;
219 m_term_VarL.clear();
220 m_term_VarT.clear();
221 m_term_CvLT.clear();
222
223 unsigned nIterSoft=0;
224 double softSumET=0.0;
225
226 // first fill the total MET
227 if(metCont->find(totalMETName)!=metCont->end()){
228 const auto &tot_met = static_cast<xAOD::MissingET*>(*(metCont->find(totalMETName)));
229 if(!MissingETBase::Source::isTotalTerm(tot_met->source())){
230 ATH_MSG_ERROR("NOT the total MET with name:" <<totalMETName);
231 return StatusCode::SUCCESS;
232 }
233 m_met = tot_met->met()/m_GeV;
234 m_metx = tot_met->mpx()/m_GeV;
235 m_mety = tot_met->mpy()/m_GeV;
236 m_metphi = tot_met->phi();
237 m_sumet = tot_met->sumet()/m_GeV;
238 m_ht = m_sumet;
239 ATH_MSG_VERBOSE("total MET: " << m_met << " phi: " << m_metphi << " name: " << tot_met->name());
240 }
241 else{
242 ATH_MSG_ERROR("Could not find the total MET with name:" <<totalMETName);
243 return StatusCode::SUCCESS;
244 }
245 m_met_vect.SetXYZ(m_metx, m_mety, 0);
246
247 // Fill the remaining terms
248 for(const auto met : *metCont) {
249
250 // skip the invisible and total MET
252 ATH_MSG_VERBOSE("Total: " << met->name() << " val: " << met->met());
253 continue;
254 }
255 if(met->source()==invisSource) continue;
256
257 // Soft term collection
259
261 ATH_MSG_VERBOSE("Soft Name: " << met->name());
262 // make sure the container name matches
263 if(met->name()!=softTermName || nIterSoft>0){
264 if(nIterSoft>0) ATH_MSG_ERROR("Found multiple soft terms with the name:" <<softTermName << ". Your MET configuration is wrong!!!");
265 continue;
266 }
267 ++nIterSoft;
268 softSumET=(met->sumet()/m_GeV);
269
270 AddSoftTerm(met, m_met_vect, particle_sum);
271 m_metsoft = met->met()/m_GeV;
272 m_metsoftphi = met->phi();
273 metTerm = 2; // this is actually filled in AddSoftTerm
274 // done with the soft term. go to the next term.
275 continue;
276 }
277 ATH_MSG_VERBOSE("Add MET term " << met->name() );
278 for(const auto& el : acc_constitObjLinks(*met)) {
279 const xAOD::IParticle* obj(*el);
280 float pt_reso=0.0, phi_reso=0.0;
281 if(!obj){
282 ATH_MSG_ERROR("Particle pointer is not valid. This will likely result in a crash " << obj);
283 return StatusCode::FAILURE;
284 }
285 ATH_MSG_VERBOSE("pT: " << obj->pt() << " type: " << obj->type() << " truth: " << (obj->type()==xAOD::Type::TruthParticle));
286 if(obj->type()==xAOD::Type::Muon || (obj->type()==xAOD::Type::TruthParticle && std::abs(static_cast<const xAOD::TruthParticle*>(obj)->pdgId())==13)){
287 ATH_CHECK(AddMuon(obj, pt_reso, phi_reso, avgmu));
288 metTerm=4;
289 }
290 else if(obj->type()==xAOD::Type::Jet){
291 // make sure the container name matches
292 if(met->name()!=jetTermName) continue;
293 ATH_CHECK(AddJet(obj, pt_reso, phi_reso, avgmu));
294 metTerm=1;
295 }
296 else if(obj->type()==xAOD::Type::Electron || (obj->type()==xAOD::Type::TruthParticle && std::abs(static_cast<const xAOD::TruthParticle*>(obj)->pdgId())==11)){
297 ATH_CHECK(AddElectron(obj, pt_reso, phi_reso, avgmu));
298 metTerm=3;
299 }
300 else if(obj->type()==xAOD::Type::Photon || (obj->type()==xAOD::Type::TruthParticle && std::abs(static_cast<const xAOD::TruthParticle*>(obj)->pdgId())==22)){
301 ATH_CHECK(AddPhoton(obj, pt_reso, phi_reso));
302 metTerm=5;
303 }
304 else if(obj->type()==xAOD::Type::Tau || (obj->type()==xAOD::Type::TruthParticle && std::abs(static_cast<const xAOD::TruthParticle*>(obj)->pdgId())==15)){
305 ATH_CHECK(AddTau(obj, pt_reso, phi_reso));
306 metTerm=6;
307 }
308
309 // compute NEW
310 double particle_u[2][2] = {{pt_reso*pt_reso*obj->pt()*obj->pt()/m_GeV/m_GeV,0.0},
311 {0.0,phi_reso*phi_reso/m_GeV/m_GeV}};
312 double particle_u_rot[2][2] = {{pt_reso*pt_reso*obj->pt()*obj->pt()/m_GeV/m_GeV,0.0},
313 {0.0,phi_reso*phi_reso/m_GeV/m_GeV}};
314 RotateXY(particle_u, particle_u_rot,m_met_vect.DeltaPhi(obj->p4().Vect()));
315 m_VarL+=particle_u_rot[0][0];
316 m_VarT+=particle_u_rot[1][1];
317 m_CvLT+=particle_u_rot[0][1];
318
319 // Save the resolutions separated for each object type
320 AddResoMap(particle_u_rot[0][0],
321 particle_u_rot[1][1],
322 particle_u_rot[0][1],
323 metTerm);
324
325 RotateXY (particle_u, particle_u_rot, obj->p4().Phi()); // positive phi rotation
326 AddMatrix(particle_sum, particle_u_rot, particle_sum);
327 // END compute NEW
328
329 ATH_MSG_VERBOSE("Resolution: " << pt_reso << " phi reso: " << phi_reso );
330 }
331 }
332
333 // setting the MET directed variables for later phi rotations if requested
337
338 if( m_VarL != 0 ){
339
340 if(m_applyBias){
341 TVector3 met_vect = m_met_vect;
342 TVector3 soft_vect = m_soft_vect;
343
344 // should be done to reset the phi as well...
346 Double_t Bias_TST = BiasPtSoftdir(m_metsoft);
347 Double_t MEx = m_met * std::cos(m_metphi) - Bias_TST * std::cos(m_metsoftphi);
348 Double_t MEy = m_met * std::sin(m_metphi) - Bias_TST * std::sin(m_metsoftphi);
349 met_vect.SetXYZ(MEx,MEy,0.0);
350 }
352 m_soft_vect.SetPtEtaPhi(m_metsoft, 0.0, m_metsoftphi);
354 Double_t PtSoftparaPH = m_pthard_vect.Mag()>0.0 ? (m_soft_vect.Dot(m_pthard_vect))/m_pthard_vect.Mag() : 0.0;
355 Double_t Bias_pthard = Bias_PtSoftParall(PtSoftparaPH);
356 Double_t MEx = m_met * std::cos(m_metphi) - Bias_pthard * std::cos(m_metsoftphi);
357 Double_t MEy = m_met * std::sin(m_metphi) - Bias_pthard * std::sin(m_metsoftphi);
358 met_vect.SetXYZ(MEx,MEy,0.0);
359 }
360 // Rotate & compute
361 ATH_CHECK(RotateToPhi(met_vect.Phi()));
363 m_rho = m_CvLT / std::sqrt( m_VarL * m_VarT ) ;
364 }
365 else{
366 // standard calculation
368 m_rho = m_CvLT / std::sqrt( m_VarL * m_VarT ) ;
369 }
370 m_ht-=softSumET;
371 ATH_MSG_VERBOSE(" Significance (squared): " << m_significance << " rho: " << GetRho()
372 << " MET: " << m_met << " phi: " << m_metphi << " SUMET: " << m_sumet << " HT: " << m_ht << " sigmaL: " << GetVarL()
373 << " sigmaT: " << GetVarT() << " MET/sqrt(SumEt): " << GetMETOverSqrtSumET()
374 << " MET/sqrt(HT): " << GetMETOverSqrtHT()
375 << " sqrt(signif): " << GetSignificance()
376 << " sqrt(signifDirectional): " << GetSigDirectional());
377 }
378 else
379 ATH_MSG_DEBUG("Var_L is 0");
380
381 return StatusCode::SUCCESS;
382 }
383
385
386 // Rotation (components)
388
389 if( m_VarL != 0 ){
391 m_rho = m_CvLT / std::sqrt( m_VarL * m_VarT ) ;
392 }
393 ATH_MSG_DEBUG(" Significance (squared) at new phi: " << m_significance
394 << " rho: " << GetRho()
395 << " MET: " << m_met
396 << " sigmaL: " << GetVarL()
397 << " sigmaT: " << GetVarT() );
398
399 return StatusCode::SUCCESS;
400 }
401
402 StatusCode METSignificance::SetLambda(const float px, const float py, const bool GeV){
403
404 // compute the new direction
405 double GeVConv = GeV ? 1.0 : m_GeV;
406 m_lamda_vect.SetXYZ(px/GeVConv, py/GeVConv, 0.0);
408 const double met_m_lamda = m_lamda_vect.Pt();
409
410 // Rotation (components)
412
413 if( m_VarL != 0 ){
415 m_rho = m_CvLT / std::sqrt( m_VarL * m_VarT ) ;
416 }
417 ATH_MSG_DEBUG(" Significance (squared) at new phi: " << m_significance
418 << " rho: " << GetRho()
419 << " MET: " << m_met
420 << " sigmaL: " << GetVarL()
421 << " sigmaT: " << GetVarT() );
422
423 return StatusCode::SUCCESS;
424 }
425
426 // Muon propagation of resolution
427 StatusCode METSignificance::AddMuon(const xAOD::IParticle* obj, float &pt_reso, float &phi_reso, float avgmu){
428
429 int dettype = 0;
430 bool DoEMuReso = false;
431 ATH_MSG_VERBOSE("Particle type: " << obj->type());
432
433 if(obj->type()==xAOD::Type::TruthParticle){
434 pt_reso =0.01;
435 if(obj->pt()>0.5e6) pt_reso=0.03;
436 if(obj->pt()>1.0e6) pt_reso=0.1;// this is just a rough estimate for the time being until the interface can handle truth muons
437 }
438 else{
439 if(!m_muonOK){
440 ATH_MSG_ERROR("MET Significance received a muon but was not configured for muons!");
441 return StatusCode::FAILURE;
442 }
443 const xAOD::Muon* muon(static_cast<const xAOD::Muon*>(obj));
444 if(muon->muonType()==xAOD::Muon::MuonType::Combined){
445 dettype=3;//CB
446 }
447 else if(muon->muonType()==xAOD::Muon::MuonType::MuonStandAlone){
448 dettype=1;//MS
449 }
450 else {//Segment, Calo, Silicon
451 dettype=2;//ID
452 }
453
454 pt_reso=m_muonCalibrationAndSmearingTool->expectedResolution(dettype,*muon,!m_isDataMuon);
455 if(m_doPhiReso) phi_reso = muon->pt()*0.001;
456 // run the jet resolution for muons. for validation region extrapolation
457 if(!m_EMuResoAux.empty()){
459 DoEMuReso = acc_EMReso.isAvailable(*muon) ? acc_EMReso(*muon) : false;
460 }
461 ATH_MSG_VERBOSE("muon: " << pt_reso << " dettype: " << dettype << " " << muon->pt() << " " << muon->p4().Eta() << " " << muon->p4().Phi());
462 }// end reco setup
463
464 // Common setup
465 if(m_doPhiReso) phi_reso = obj->pt()*0.001;
466 ATH_MSG_VERBOSE("muon: " << pt_reso << " dettype: " << dettype << " " << obj->pt() << " " << obj->p4().Eta() << " " << obj->p4().Phi());
467
468 if(m_jerForEMu || DoEMuReso){
469 bool treatPUJets = m_treatPUJets;
470 m_treatPUJets=false; //turn off pileup jet treatement for this electron
471 ATH_CHECK(AddJet(obj, pt_reso, phi_reso, avgmu));
472 m_treatPUJets = treatPUJets; // reset value
473 }
474
475 return StatusCode::SUCCESS;
476 }
477
478 // Electron propagation of resolution
479 StatusCode METSignificance::AddElectron(const xAOD::IParticle* obj, float &pt_reso, float &phi_reso, float avgmu){
480
481 if(!m_egammaOK){
482 ATH_MSG_ERROR("MET Significance received an electron but was not configured for egamma!");
483 return StatusCode::FAILURE;
484 }
485
486 bool DoEMuReso = false;
487 if(obj->type()==xAOD::Type::TruthParticle){
488 pt_reso=m_egammaCalibTool->resolution(obj->e(),obj->eta(),obj->eta(),PATCore::ParticleType::Electron);
489 if(m_doPhiReso) phi_reso = obj->pt()*0.004;
490 }
491 else{
492 const xAOD::Electron* ele(static_cast<const xAOD::Electron*>(obj));
493 const auto cl_etaCalo = xAOD::get_eta_calo(*(ele->caloCluster()), ele->author());
494 pt_reso=m_egammaCalibTool->resolution(ele->e(),ele->caloCluster()->eta(),cl_etaCalo,PATCore::ParticleType::Electron);
495 if(m_doPhiReso) phi_reso = ele->pt()*0.004;
496 ATH_MSG_VERBOSE("el: " << pt_reso << " " << ele->pt() << " " << ele->p4().Eta() << " " << ele->p4().Phi());
497
498 // run the jet resolution for muons. for validation region extrapolation
499 if(!m_EMuResoAux.empty()){
501 DoEMuReso = acc_EMReso.isAvailable(*ele) ? acc_EMReso(*ele) : false;
502 }
503 }
504
505 if(m_jerForEMu || DoEMuReso){
506 bool treatPUJets = m_treatPUJets;
507 m_treatPUJets=false; //turn off pileup jet treatement for this electron
508 ATH_CHECK(AddJet(obj, pt_reso, phi_reso, avgmu));
509 m_treatPUJets = treatPUJets; // reset value
510 }
511 return StatusCode::SUCCESS;
512 }
513
514 // Photon propagation of resolution
515 StatusCode METSignificance::AddPhoton(const xAOD::IParticle* obj, float &pt_reso, float &phi_reso){
516
517 if(!m_egammaOK){
518 ATH_MSG_ERROR("MET Significance received a photon but was not configured for egamma!");
519 return StatusCode::FAILURE;
520 }
521
522 if(obj->type()==xAOD::Type::TruthParticle){
523 pt_reso=m_egammaCalibTool->resolution(obj->e(),obj->eta(),obj->eta(),PATCore::ParticleType::Electron); // leaving as an electron for the truth implementation rather than declaring a reco photon
524 if(m_doPhiReso) phi_reso = obj->pt()*0.004;
525 }
526 else{
527 const xAOD::Egamma* pho(static_cast<const xAOD::Egamma*>(obj));
528 pt_reso=m_egammaCalibTool->getResolution(*pho);
529 if(m_doPhiReso) phi_reso = pho->pt()*0.004;
530 ATH_MSG_VERBOSE("pho: " << pt_reso << " " << pho->pt() << " " << pho->p4().Eta() << " " << pho->p4().Phi());
531 }
532 return StatusCode::SUCCESS;
533 }
534
535 // Jet propagation of resolution. returns the relative pT and phi resolution.
536 StatusCode METSignificance::AddJet(const xAOD::IParticle* obj, float &pt_reso, float &phi_reso, float &avgmu){
537
538 const xAOD::Jet* jet(static_cast<const xAOD::Jet*>(obj));
539 double pt_reso_dbl_data=0.0, pt_reso_dbl_mc=0.0, pt_reso_dbl_max=0.0;
540
541 // setting limits on jets if requested
542 if(m_jetPtThr>0.0 && m_jetPtThr>jet->pt()) return StatusCode::SUCCESS;
543 if(m_jetEtaThr>0.0 && m_jetEtaThr<std::abs(jet->eta())) return StatusCode::SUCCESS;
544
545 if(!m_jetOK){
546 ATH_MSG_ERROR("MET Significance received a jet but was not configured for jets!");
547 return StatusCode::FAILURE;
548 }
549
550 ATH_CHECK(m_jetCalibTool->getNominalResolutionData(*jet, pt_reso_dbl_data));
551 ATH_CHECK(m_jetCalibTool->getNominalResolutionMC(*jet, pt_reso_dbl_mc));
552 pt_reso_dbl_max = std::max(pt_reso_dbl_data,pt_reso_dbl_mc);
553 pt_reso = pt_reso_dbl_max;
554
555 ATH_MSG_VERBOSE("jet: " << pt_reso << " jetpT: " << jet->pt() << " " << jet->p4().Eta() << " " << jet->p4().Phi());
556
557 // Add extra uncertainty for PU jets based on JVT
558 if(m_treatPUJets){
559 double jet_pu_unc = 0.;
560 if(acc_fjvt.isAvailable(*jet))
561 jet_pu_unc = GetPUProb(jet->eta(), jet->phi(),jet->pt()/m_GeV, acc_jvt(*jet), acc_fjvt(*jet), avgmu);
562 else if(acc_fjvt_der.isAvailable(*jet))
563 jet_pu_unc = GetPUProb(jet->eta(), jet->phi(),jet->pt()/m_GeV, acc_jvt(*jet), acc_fjvt_der(*jet), avgmu);
564 else{
565 ATH_MSG_ERROR("No fJVT decoration available - must have treat pileup jets set to off or provide fJVT!");
566 return StatusCode::FAILURE;
567 }
568 pt_reso = std::sqrt(jet_pu_unc*jet_pu_unc + pt_reso*pt_reso);
569 ATH_MSG_VERBOSE("jet_pu_unc: " << jet_pu_unc);
570 }
571
572 // Use the phi resolution of the jets
573 // needs to be finished
574 if(m_doPhiReso){
575 double jet_phi_unc = std::abs(GetPhiUnc(jet->eta(), jet->phi(),jet->pt()/m_GeV));
576 phi_reso = jet->pt()*jet_phi_unc;
577 }
578
579 // Add user defined additional resolutions. For example, b-tagged jets
580 if(!m_JetResoAux.empty()){
582 if(acc_extra.isAvailable(*jet)){
583 float extra_relative_pt_reso = acc_extra(*jet);
584 pt_reso = std::sqrt(pt_reso*pt_reso + extra_relative_pt_reso*extra_relative_pt_reso);
585 }
586 }
587
588 return StatusCode::SUCCESS;
589 }
590
591 // Tau propagation of resolution
592 StatusCode METSignificance::AddTau(const xAOD::IParticle* obj, float &pt_reso, float &phi_reso){
593
594 // tau objects
595 if(obj->type()==xAOD::Type::TruthParticle){
596 pt_reso= 0.1;
597 if(m_doPhiReso) phi_reso = obj->pt()*0.01;
598 }
599 else{
600 if(!m_tauOK){
601 ATH_MSG_ERROR("MET Significance received a tau but was not configured for taus!");
602 return StatusCode::FAILURE;
603 }
604 const xAOD::TauJet* tau(static_cast<const xAOD::TauJet*>(obj));
605 if (auto *combp4 = dynamic_cast<TauAnalysisTools::TauResolutionTool*>(&*m_tauResolutionTool)) {
606 pt_reso = combp4->getMvaEnergyResolution(*tau);
607 }
608
609 if(m_doPhiReso) phi_reso = tau->pt()*0.01;
610 ATH_MSG_VERBOSE("tau: " << pt_reso << " " << tau->pt() << " " << tau->p4().Eta() << " " << tau->p4().Phi() << " phi reso: " << phi_reso);
611 }
612 return StatusCode::SUCCESS;
613 }
614
615 //
616 // Soft term propagation of resolution
617 //
618 void METSignificance::AddSoftTerm(const xAOD::MissingET* soft, const TVector3 &met_vect, double (&particle_sum)[2][2]){
619
621
622 ATH_MSG_VERBOSE("Resolution Soft term set to 10GeV");
623
624 m_soft_vect.SetXYZ(soft->mpx()/m_GeV, soft->mpy()/m_GeV, 0);
625
626 double particle_u[2][2] = {{m_softTermReso*m_softTermReso,0.0},
628 double particle_u_rot[2][2] = {{m_softTermReso*m_softTermReso,0.0},
630
631 RotateXY(particle_u, particle_u_rot,met_vect.DeltaPhi(m_soft_vect));
632 m_VarL+=particle_u_rot[0][0];
633 m_VarT+=particle_u_rot[1][1];
634 m_CvLT+=particle_u_rot[0][1];
635
636 // Save the resolutions separated for each object type
637 AddResoMap(particle_u_rot[0][0],
638 particle_u_rot[1][1],
639 particle_u_rot[0][1],
641
642 RotateXY (particle_u, particle_u_rot,-1.0*soft->phi()); // negative phi rotation
643 AddMatrix(particle_sum, particle_u_rot, particle_sum);
644
645 ATH_MSG_VERBOSE("SOFT " << soft->name() <<" - pt_reso: " << m_softTermReso << " soft: " << soft->met() << " phi: " << soft->phi()
646 << " Var_L: " << particle_u_rot[0][0] << " Var_T: " << particle_u_rot[1][1]
647 << " " << particle_u_rot[0][1]);
648 }
650
651 ATH_MSG_VERBOSE("Resolution Soft term parameterized in pthard direction");
652
653 m_soft_vect.SetXYZ(soft->mpx()/m_GeV, soft->mpy()/m_GeV, 0);
654
655 m_pthard_vect = m_soft_vect - met_vect;
656
657 double varTST = Var_Ptsoft(soft->met()/m_GeV);
658
659 double particle_u[2][2] = {{varTST,0.0},
660 {0.0,varTST}};
661 double particle_u_rot[2][2] = {{varTST,0.0},
662 {0.0,varTST}};
663
664 RotateXY(particle_u, particle_u_rot,met_vect.DeltaPhi(m_pthard_vect));
665 m_VarL+=particle_u_rot[0][0];
666 m_VarT+=particle_u_rot[1][1];
667 m_CvLT+=particle_u_rot[0][1];
668
669 // Save the resolutions separated for each object type
670 AddResoMap(particle_u_rot[0][0],
671 particle_u_rot[1][1],
672 particle_u_rot[0][1],
674
675 RotateXY (particle_u, particle_u_rot,-1.0*m_pthard_vect.Phi()); // negative phi rotation
676 AddMatrix(particle_sum, particle_u_rot, particle_sum);
677
678 }
680
681 ATH_MSG_VERBOSE("Resolution Soft term parameterized in TST");
682
683 m_soft_vect.SetXYZ(soft->mpx()/m_GeV, soft->mpy()/m_GeV, 0);
684
685 double varTST = VarparPtSoftdir(soft->met()/m_GeV, soft->sumet()/m_GeV);
686
687 double particle_u[2][2] = {{varTST,0.0},
688 {0.0,varTST}};
689 double particle_u_rot[2][2] = {{varTST,0.0},
690 {0.0,varTST}};
691
692 RotateXY(particle_u, particle_u_rot,met_vect.DeltaPhi(m_soft_vect));
693 m_VarL+=particle_u_rot[0][0];
694 m_VarT+=particle_u_rot[1][1];
695 m_CvLT+=particle_u_rot[0][1];
696
697 // Save the resolutions separated for each object type
698 AddResoMap(particle_u_rot[0][0],
699 particle_u_rot[1][1],
700 particle_u_rot[0][1],
702
703 RotateXY (particle_u, particle_u_rot,-1.0*soft->phi()); // negative phi rotation
704 AddMatrix(particle_sum, particle_u_rot, particle_sum);
705
706 }
707 else{
708 ATH_MSG_ERROR("Soft term parameterization is NOT defined for:" << m_softTermParam);
709 }
710
711 }
712
713 double METSignificance::GetPUProb(double jet_eta, double /*jet_phi*/,
714 double jet_pt, double jet_jvt,
715 double jet_fjvt,
716 float avgmu) {
717
718 double unc=0.0;
719
720 // Coefficients from Doug Schaefer <schae@cern.ch> and the MET subgroup
721 if(m_JetCollection == "AntiKt4EMTopoJets"){
722 if(std::abs(jet_eta)<2.4){
723 if(jet_pt<30){
724 if(jet_jvt<0.11) unc = 1;
725 else if(jet_jvt<0.25) unc = 0.0730 + 0.0024 * avgmu + 0.00001 * avgmu * avgmu;
726 else if(jet_jvt<0.85) unc = 0.0995 + 0.0031 * avgmu + 0.00005 * avgmu * avgmu;
727 else if(jet_jvt<0.95) unc = 0.0311 + 0.0025 * avgmu + 0.00005 * avgmu * avgmu;
728 else unc = 0.0308 -0.0010 * avgmu + 0.00006 * avgmu * avgmu ;
729 }else if(jet_pt<40){
730 if(jet_jvt<0.11) unc = 1.;
731 else if(jet_jvt<0.25) unc = 1.;
732 else if(jet_jvt<0.85) unc = -0.0188 + 0.0039 * avgmu + 0.00002 * avgmu * avgmu;
733 else if(jet_jvt<0.95) unc = 0.0252 -0.0009 * avgmu + 0.00006 * avgmu * avgmu ;
734 else unc = 0.0085 -0.0003 * avgmu + 0.00002 * avgmu * avgmu ;
735 }else if(jet_pt<50){
736 if(jet_jvt<0.11) unc = 1;
737 else if(jet_jvt<0.25) unc = 0.0345 -0.0006 * avgmu + 0.00004 * avgmu * avgmu ;
738 else if(jet_jvt<0.85) unc = 0.1078 -0.0051 * avgmu + 0.00011 * avgmu * avgmu ;
739 else if(jet_jvt<0.95) unc = -0.0026 + 0.0005 * avgmu + 0.00002 * avgmu * avgmu;
740 else unc = 0.0090 -0.0004 * avgmu + 0.00001 * avgmu * avgmu ;
741 }else if(jet_pt<60){
742 if(jet_jvt<0.11) unc = 1;
743 else if(jet_jvt<0.25) unc = -0.0321 + 0.0030 * avgmu -0.00002 * avgmu * avgmu;
744 else if(jet_jvt<0.85) unc = 0.0260 -0.0007 * avgmu + 0.00003 * avgmu * avgmu ;
745 else unc = -0.0040 + 0.0003 * avgmu;
746 }else if(jet_pt<100){
747 unc = 0.9492 -2.0757 * jet_jvt + 1.13328 * jet_jvt * jet_jvt;
748 }else if(jet_pt<150){
749 unc = 0.7888 -1.8372 * jet_jvt + 1.05539 * jet_jvt * jet_jvt;
750 }
751 }else if(std::abs(jet_eta)<2.6){
752 if(jet_pt<30){
753 if(jet_jvt<0.11) unc = 0.2633 + 0.0091 * avgmu + -0.00009 * avgmu * avgmu;
754 else if(jet_jvt<0.25) unc = 0.1841 + 0.0144 * avgmu + -0.00008 * avgmu * avgmu;
755 else if(jet_jvt<0.85) unc = 0.1401 + 0.0048 * avgmu + 0.00006 * avgmu * avgmu ;
756 else if(jet_jvt<0.95) unc = -0.0118 + 0.0076 * avgmu + 0.00003 * avgmu * avgmu;
757 else unc = 0.0534 + -0.0011 * avgmu + 0.00010 * avgmu * avgmu;
758 }else if(jet_pt<40){
759 if(jet_jvt<0.11) unc = 0.1497 + 0.0133 * avgmu + -0.00015 * avgmu * avgmu ;
760 else if(jet_jvt<0.25) unc = -0.2260 + 0.0276 * avgmu + -0.00021 * avgmu * avgmu ;
761 else if(jet_jvt<0.85) unc = 0.2743 + -0.0093 * avgmu + 0.00022 * avgmu * avgmu ;
762 else if(jet_jvt<0.95) unc = 0.0604 + 0.0006 * avgmu + 0.00006 * avgmu * avgmu ;
763 else unc = 0.0478 + -0.0009 * avgmu + 0.00004 * avgmu * avgmu ;
764 }else if(jet_pt<50){
765 if(jet_jvt<0.11) unc = -0.2187 + 0.0317 * avgmu + -0.00037 * avgmu * avgmu ;
766 else if(jet_jvt<0.25) unc = 0.0964 + 0.0053 * avgmu + 0.00002 * avgmu * avgmu ;
767 else if(jet_jvt<0.85) unc = 1.1730 + -0.0624 * avgmu + 0.00088 * avgmu * avgmu ;
768 else if(jet_jvt<0.95) unc = -0.2011 + 0.0151 * avgmu + -0.00018 * avgmu * avgmu ;
769 else unc = 0.0145 + -0.0003 * avgmu + 0.00002 * avgmu * avgmu ;
770 }else if(jet_pt<60){
771 if(jet_jvt<0.11) unc = 0.0051 + 0.0113 * avgmu + -0.00008 * avgmu * avgmu ;
772 else if(jet_jvt<0.25) unc = -0.1024 + 0.0109 * avgmu + -0.00006 * avgmu * avgmu ;
773 else if(jet_jvt<0.85) unc = 1.2491 + -0.0501 * avgmu + 0.00052 * avgmu * avgmu ;
774 else unc = 0.0267 + -0.0014 * avgmu + 0.00003 * avgmu * avgmu ;
775 }else if(jet_pt<100){
776 unc = 0.8951 -2.4995 * jet_jvt + 1.63229 * jet_jvt * jet_jvt;
777 }else if(jet_pt<150){
778 unc = 0.9998 -1.7319 * jet_jvt + 0.72680 * jet_jvt * jet_jvt;
779 }
780 }else if(std::abs(jet_eta)<2.7){
781 if(jet_pt<30){
782 if(jet_jvt<0.11) unc = 0.3001 + 0.0054 * avgmu -0.00004 * avgmu * avgmu ;
783 else if(jet_jvt<0.25) unc = 0.0663 + 0.0198 * avgmu -0.00013 * avgmu * avgmu ;
784 else if(jet_jvt<0.85) unc = -0.0842 + 0.0163 * avgmu -0.00008 * avgmu * avgmu ;
785 else if(jet_jvt<0.95) unc = -0.0219 + 0.0080 * avgmu + 0.00003 * avgmu * avgmu;
786 else unc = 0.0461 -0.0003 * avgmu + 0.00012 * avgmu * avgmu ;
787 }else if(jet_pt<40){
788 if(jet_jvt<0.11) unc = 0.1885 + 0.0083 * avgmu -0.00006 * avgmu * avgmu ;
789 else if(jet_jvt<0.25) unc = -0.0286 + 0.0150 * avgmu -0.00007 * avgmu * avgmu;
790 else if(jet_jvt<0.85) unc = 0.0152 + 0.0028 * avgmu + 0.00005 * avgmu * avgmu;
791 else if(jet_jvt<0.95) unc = 0.1815 -0.0076 * avgmu + 0.00018 * avgmu * avgmu ;
792 else unc = 0.0192 -0.0003 * avgmu + 0.00007 * avgmu * avgmu ;
793 }else if(jet_pt<50){
794 if(jet_jvt<0.11) unc = 0.1257 + 0.0074 * avgmu -0.00004 * avgmu * avgmu ;
795 else if(jet_jvt<0.25) unc = -0.0276 + 0.0080 * avgmu + 0.00000 * avgmu * avgmu;
796 else if(jet_jvt<0.85) unc = 0.1403 -0.0051 * avgmu + 0.00009 * avgmu * avgmu ;
797 else if(jet_jvt<0.95) unc = 0.2078 -0.0101 * avgmu + 0.00017 * avgmu * avgmu ;
798 else unc = 0.2597 -0.0132 * avgmu + 0.00020 * avgmu * avgmu ;
799 }else if(jet_pt<60){
800 if(jet_jvt<0.11) unc = 0.1111 + 0.0045 * avgmu -0.00000 * avgmu * avgmu ;
801 else if(jet_jvt<0.25) unc = 0.0975 -0.0011 * avgmu + 0.00008 * avgmu * avgmu ;
802 else if(jet_jvt<0.85) unc = 0.0920 -0.0053 * avgmu + 0.00013 * avgmu * avgmu ;
803 else unc = -0.0071 + 0.0016 * avgmu -0.00001 * avgmu * avgmu;
804 }else if(jet_pt<100){
805 unc = 0.4660 -1.2116 * jet_jvt + 0.78807 * jet_jvt * jet_jvt;
806 }else if(jet_pt<150){
807 unc = 0.2254 -0.5476 * jet_jvt + 0.32617 * jet_jvt * jet_jvt;
808 }
809 }// end eta 2.7
810 else{//forward jets
811 float fjvt = jet_fjvt>0.6 ? 0.6 : jet_fjvt; // the pileup more or less plateaus at 0.6
812 if(jet_pt<30) unc = 0.5106 + 1.2566 * fjvt -1.15060 * fjvt * fjvt;
813 else if(jet_pt<40) unc = 0.2972 + 1.9418 * fjvt -1.82694 * fjvt * fjvt;
814 else if(jet_pt<50) unc = 0.1543 + 1.9864 * fjvt -1.48429 * fjvt * fjvt;
815 else if(jet_pt<60) unc = 0.1050 + 1.3196 * fjvt + 0.03554 * fjvt * fjvt;
816 else if(jet_pt<120) unc = 0.0400 + 0.5653 * fjvt + 1.96323 * fjvt * fjvt;
817 // max of 0.9 seems reasonable
818 if(jet_fjvt>0.6) unc = 0.9;
819 }
820 // end emtopo
821 // Coefficients from Badr-eddine Ngair <badr-eddine.ngair@cern.ch>
822 // Pile-UP (resolution) estimatetd using Z->ee Events using exclusive jet_pt binning using RUN3 samples
823 }else{//p-flow inputs
824 if(std::abs(jet_eta)<2.4){
825 if(jet_pt<30){
826 if(jet_jvt<0.11) unc = 0.524466 + 0.00750057 * avgmu -4.73422e-05 * avgmu * avgmu ;
827 else if(jet_jvt<0.25) unc = 4.17584e-01 + 1.00112e-02 * avgmu -7.43546e-05 * avgmu * avgmu ;
828 else if(jet_jvt<0.85) unc = 2.12625e-01 + 1.03484e-02 * avgmu -5.68063e-05 * avgmu * avgmu;
829 else if(jet_jvt<0.95) unc = 1.08396e-01 + 1.04273e-02 * avgmu -5.00299e-05 * avgmu * avgmu;
830 else unc = 1.26304e-03 +2.10385e-04 * avgmu + 1.10086e-06 * avgmu * avgmu ;
831 }else if(jet_pt<40){
832 if(jet_jvt<0.11) unc = 3.78090e-01 + 8.83535e-03 * avgmu -5.38873e-05* avgmu * avgmu;
833 else if(jet_jvt<0.25) unc = 2.67244e-01+ 9.81193e-03* avgmu -5.87765e-05 * avgmu * avgmu;
834 else if(jet_jvt<0.85) unc = 9.60892e-02 +8.11069e-03 * avgmu -3.73101e-05 * avgmu * avgmu ;
835 else if(jet_jvt<0.95) unc = 5.16235e-02 +6.27371e-03 * avgmu -1.95433e-05 * avgmu * avgmu ;
836 else unc = -2.74714e-03 +2.45273e-04* avgmu -1.44731e-06* avgmu * avgmu ;
837 }else if(jet_pt<50){
838 if(jet_jvt<0.11) unc = 2.44953e-01 + 1.23246e-02 * avgmu -8.48696e-05 * avgmu * avgmu;
839 else if(jet_jvt<0.25) unc = 1.78141e-01 + 1.09451e-02 * avgmu -6.90796e-05 * avgmu * avgmu;
840 else if(jet_jvt<0.85) unc = 9.60998e-02+6.21945e-03* avgmu -2.76203e-05* avgmu * avgmu ;
841 else if(jet_jvt<0.95) unc = 5.79210e-02+ 4.49780e-03 * avgmu + -1.15125e-05 * avgmu * avgmu ;
842 else unc = -2.96644e-03 +2.27707e-04 * avgmu -1.86712e-06 * avgmu * avgmu ;
843 }else if(jet_pt<60){
844 if(jet_jvt<0.11) unc = 1.97017e-01 + 1.34089e-02 * avgmu -9.18923e-05 * avgmu * avgmu ;
845 else if(jet_jvt<0.25) unc = 1.27602e-01 + 1.12287e-02 * avgmu -6.62192e-05 * avgmu * avgmu ;
846 else if(jet_jvt<0.85) unc = 6.94905e-02 + 6.27784e-03 * avgmu -3.07298e-05 * avgmu * avgmu;
847 else if(jet_jvt<0.95) unc = 3.58417e-02 +4.62268e-03 * avgmu -2.12417e-05 * avgmu * avgmu ;
848 else unc = 1.35616e-03 + 5.46723e-06* avgmu + 1.92327e-07 * avgmu * avgmu;
849 }else if(jet_pt<100){
850 unc = 6.19009e-01 -8.96042e-01 * jet_jvt + 2.89066e-01 * jet_jvt * jet_jvt;
851 }else if(jet_pt<150){
852 unc = 6.18350e-01 -8.97327e-01 * jet_jvt + 2.90998e-01 * jet_jvt * jet_jvt;
853 }
854 }else if(std::abs(jet_eta)<2.6){
855 if(jet_pt<30){
856 if(jet_jvt<0.11) unc = 5.06496e-01 + 8.21123e-03 * avgmu -5.17501e-05 * avgmu * avgmu;
857 else if(jet_jvt<0.25) unc = 4.26616e-01 + 9.25936e-03 * avgmu -5.68847e-05 * avgmu * avgmu;
858 else if(jet_jvt<0.85) unc = 2.03333e-01 + 1.11951e-02 * avgmu-6.09233e-05 * avgmu * avgmu ;
859 else if(jet_jvt<0.95) unc = 1.03167e-01 + 1.13444e-02 * avgmu -5.43274e-05 * avgmu * avgmu;
860 else unc = 1.51480e-03 + 2.08394e-04 * avgmu + 1.39579e-06 * avgmu * avgmu;
861 }else if(jet_pt<40){
862 if(jet_jvt<0.11) unc = 3.40612e-01 + 9.94199e-03 * avgmu + -5.93760e-05* avgmu * avgmu ;
863 else if(jet_jvt<0.25) unc = 2.43360e-01 + 1.05579e-02* avgmu + -6.05403e-05* avgmu * avgmu ;
864 else if(jet_jvt<0.85) unc = 8.34364e-02 + 8.76364e-03 * avgmu -3.64035e-05 * avgmu * avgmu ;
865 else if(jet_jvt<0.95) unc = 4.40362e-02 + 6.92580e-03 * avgmu -1.79853e-05 * avgmu * avgmu ;
866 else unc = -2.68670e-03 + 2.50861e-04 * avgmu + -1.46410e-06* avgmu * avgmu ;
867 }else if(jet_pt<50){
868 if(jet_jvt<0.11) unc = 2.36561e-01 + 1.14078e-02 * avgmu +-7.10025e-05 * avgmu * avgmu ;
869 else if(jet_jvt<0.25) unc = 1.86653e-01 +9.61140e-03 * avgmu +-5.15356e-05 * avgmu * avgmu ;
870 else if(jet_jvt<0.85) unc = 9.37026e-02 +5.93028e-03 * avgmu +-2.02571e-05 * avgmu * avgmu ;
871 else if(jet_jvt<0.95) unc = 5.79210e-02+ 4.49780e-03 * avgmu + -1.15125e-05 * avgmu * avgmu ;
872 else unc = -3.02487e-03 + 2.31337e-04* avgmu + -1.85225e-06 * avgmu * avgmu ;
873 }else if(jet_pt<60){
874 if(jet_jvt<0.11) unc =1.75215e-01+ 1.21805e-02 * avgmu + -7.48846e-05 * avgmu * avgmu ;
875 else if(jet_jvt<0.25) unc = 1.26276e-01+ 9.80117e-03 * avgmu + -4.99913e-05 * avgmu * avgmu ;
876 else if(jet_jvt<0.85) unc = 7.91422e-02 + 5.26009e-03 * avgmu +-1.87388e-05 * avgmu * avgmu ;
877 else if(jet_jvt<0.95) unc = 4.39136e-02+ 4.09435e-03* avgmu +-1.35926e-05 * avgmu * avgmu ;
878 else unc = 1.21410e-03 + 1.14188e-05 * avgmu + 1.53654e-07 * avgmu * avgmu ;
879 }else if( jet_pt<100){
880 unc = 6.44179e-01 -9.20194e-01* jet_jvt + 2.89686e-01 * jet_jvt * jet_jvt;
881 }else if(jet_pt<150){
882 unc = 6.43423e-01 -9.21407e-01 * jet_jvt + 2.91648e-01 * jet_jvt * jet_jvt;
883 }
884 }else if(std::abs(jet_eta)<2.7){
885 if(jet_pt<30){
886 if(jet_jvt<0.11) unc = 4.76243e-01 + 9.22046e-03 * avgmu -5.88765e-05 * avgmu * avgmu;
887 else if(jet_jvt<0.25) unc = 4.07406e-01+ 1.01167e-02 * avgmu -6.30429e-05 * avgmu * avgmu;
888 else if(jet_jvt<0.85) unc = 2.01324e-01+ 1.20631e-02 * avgmu -6.75582e-05 * avgmu * avgmu;
889 else if(jet_jvt<0.95) unc = 1.03815e-01 + 1.24007e-02 * avgmu -6.26892e-05 * avgmu * avgmu;
890 else unc = 1.63714e-03 +2.00682e-04 * avgmu + 1.53621e-06 * avgmu * avgmu;
891 }else if(jet_pt<40){
892 if(jet_jvt<0.11) unc = 2.89505e-01 + 1.11643e-02 * avgmu -6.45475e-05* avgmu * avgmu;
893 else if(jet_jvt<0.25) unc = 2.15968e-01 + 1.14451e-02* avgmu -6.32545e-05 * avgmu * avgmu;
894 else if(jet_jvt<0.85) unc = 7.92319e-02 +9.66239e-03* avgmu -3.81872e-05 * avgmu * avgmu ;
895 else if(jet_jvt<0.95) unc = 4.25501e-02 +7.90022e-03 * avgmu -1.93561e-05 * avgmu * avgmu;
896 else unc = -2.68089e-03 +2.50117e-04 * avgmu -1.43591e-06 * avgmu * avgmu ;
897 }else if(jet_pt<50){
898 if(jet_jvt<0.11) unc = 1.66062e-01+ 1.21029e-02 * avgmu -7.00743e-05 * avgmu * avgmu;
899 else if(jet_jvt<0.25) unc = 1.36874e-01+ 1.03543e-02 * avgmu -5.13482e-05 * avgmu * avgmu;
900 else if(jet_jvt<0.85) unc = 7.24015e-02 +6.72611e-03 * avgmu -1.98316e-05 * avgmu * avgmu;
901 else if(jet_jvt<0.95) unc = 4.84508e-02 +5.10548e-03* avgmu -8.65067e-06 * avgmu * avgmu;
902 else unc = -3.02710e-03 +2.31422e-04* avgmu -1.84795e-06 * avgmu * avgmu;
903 }else if(jet_pt<60){
904 if(jet_jvt<0.11) unc =1.33626e-01 + 1.02813e-02 * avgmu -4.87698e-05 * avgmu * avgmu;
905 else if(jet_jvt<0.25) unc = 1.06724e-01 + 8.45131e-03 * avgmu -3.26833e-05 * avgmu * avgmu;
906 else if(jet_jvt<0.85) unc = 7.47468e-02 +4.85387e-03 * avgmu -1.01430e-05* avgmu * avgmu;
907 else if(jet_jvt<0.95) unc = 4.57521e-02 + 3.86782e-03 * avgmu -6.38948e-06 * avgmu * avgmu;
908 else unc = 1.20495e-03 +1.18941e-05* avgmu +1.47846e-07 * avgmu * avgmu;
909 }else if( jet_pt<100){
910 unc = 6.59079e-01 -9.29754e-01 * jet_jvt + 2.83653e-01 * jet_jvt * jet_jvt;
911 }else if(jet_pt<150){
912 unc = 6.58295e-01 -9.31032e-01 * jet_jvt + 2.85724e-01 * jet_jvt * jet_jvt;
913 }
914 }// end eta 2.7
915 else{//forward jets
916 float fjvt = jet_fjvt>0.6 ? 0.6 : jet_fjvt; // the pileup more or less plateaus at 0.6
917 if(jet_pt<30) unc = 0.605329 + 0.625734 * fjvt -0.42484 * fjvt * fjvt;
918 else if(jet_pt<40) unc = 0.409696 + 1.00173 * fjvt -0.609179 * fjvt * fjvt;
919 else if(jet_pt<50) unc = 0.173755 + 1.48847 * fjvt -0.803771 * fjvt * fjvt;
920 else if(jet_pt<60) unc = 0.0140303 + 1.79909 * fjvt -0.889274 * fjvt * fjvt;
921 else if(jet_pt<120) unc = -0.0828333 + 1.81167 * fjvt -0.716881 * fjvt * fjvt;
922 // max of 0.9 seems reasonable
923 if(jet_fjvt>0.6) unc = 0.9;
924 }
925 }// end pflow
926
927 unc = std::min(unc, 1.0);
928 unc = std::max(unc, 0.0);
929
930 return unc;
931 }
932
933 double METSignificance::GetPhiUnc(double jet_eta, double jet_phi,double jet_pt){
934
935 unsigned int xbin = getEtaBin(jet_eta);
936 unsigned int ybin = jet_phi>0.0 ? int(jet_phi/0.4)+9 : int(jet_phi/0.4)+8;
937
938 // Stored as bin content = Mean, error = RMS, we want to use the RMS.
940 ATH_MSG_ERROR("Jet Phi Resolution histograms are invalid.");
941 return 0.0;
942 }
943
944 // Collect the phi resolution
945 if(jet_pt<50.0)
946 return m_phi_reso_pt20->GetBinError(xbin, ybin);
947 else if(jet_pt<100.0)
948 return m_phi_reso_pt50->GetBinError(xbin, ybin);
949 return m_phi_reso_pt100->GetBinError(xbin, ybin);
950 }
951
952 unsigned int METSignificance::getEtaBin(double jet_eta){
953 // For the phi uncertainty lookup
954 if(-4.5<jet_eta && -3.8>=jet_eta) return 1;
955 else if(-3.8<jet_eta && -3.5>=jet_eta) return 2;
956 else if(-3.5<jet_eta && -3.0>=jet_eta) return 3;
957 else if(-3.0<jet_eta && -2.7>=jet_eta) return 4;
958 else if(-2.7<jet_eta && -2.4>=jet_eta) return 5;
959 else if(-2.4<jet_eta && -1.5>=jet_eta) return 6;
960 else if(-1.5<jet_eta && -0.5>=jet_eta) return 7;
961 else if(-0.5<jet_eta && 0.0>=jet_eta) return 8;
962 else if(0.0<jet_eta && 0.5>=jet_eta) return 9;
963 else if(0.5<jet_eta && 1.5>=jet_eta) return 10;
964 else if(1.5<jet_eta && 2.4>=jet_eta) return 11;
965 else if(2.4<jet_eta && 2.7>=jet_eta) return 12;
966 else if(2.7<jet_eta && 3.0>=jet_eta) return 13;
967 else if(3.0<jet_eta && 3.5>=jet_eta) return 14;
968 else if(3.5<jet_eta && 3.8>=jet_eta) return 15;
969 else if(3.8<jet_eta ) return 16;
970 return 0;
971 }
972
973 std::tuple<double,double,double> METSignificance::CovMatrixRotation(double var_x, double var_y, double cv_xy, double Phi){
974 // Covariance matrix parallel and transverse to the Phi direction
975 Double_t V11 = std::pow(std::cos(Phi),2)*var_x + 2*std::sin(Phi)*std::cos(Phi)*cv_xy + std::pow(std::sin(Phi),2)*var_y;
976 Double_t V22 = std::pow(std::sin(Phi),2)*var_x - 2*std::sin(Phi)*std::cos(Phi)*cv_xy + std::pow(std::cos(Phi),2)*var_y;
977 Double_t V12 = std::pow(std::cos(Phi),2)*cv_xy -std::sin(Phi)*std::cos(Phi)*var_x + std::sin(Phi)*std::cos(Phi)*var_y - std::pow(std::sin(Phi),2)*cv_xy; // rho is equal to one for just one jet
978 return std::make_tuple( V11, V22, V12);
979 }
980
981 double METSignificance::Significance_LT(double Numerator, double var_parall, double var_perpen, double cov){
982
983 Double_t rho = cov / std::sqrt( var_parall * var_perpen ) ;
984 Double_t Significance = 0;
985 if (std::abs( rho ) >= 0.9 ){ //Cov Max not invertible -> Significance diverges
986 ATH_MSG_VERBOSE("rho is large: " << rho);
987 Significance = std::pow( Numerator - m_scalarBias , 2 ) / ( var_parall ) ;
988 }
989 else
990 Significance = std::pow( Numerator - m_scalarBias , 2 ) / ( var_parall * ( 1 - std::pow(rho,2) ) ) ;
991
992 if( std::abs(Significance) >= 10e+15)
993 ATH_MSG_WARNING("warning -->"<< Significance);
994
995 return Significance;
996 }
997
998 void METSignificance::InvertMatrix(double (&mat)[2][2], double (&m)[2][2]){
999
1000 // determinant
1001 double det = mat[0][0]*mat[1][1]-mat[0][1]*mat[1][0];
1002
1003 m[0][0]=0.0;
1004 m[0][1]=0.0;
1005 m[1][0]=0.0;
1006 m[1][1]=0.0;
1007
1008 if(det==0.0) return;
1009
1010 m[0][0]= 1.0/det*(mat[1][1]);
1011 m[1][0]=-1.0/det*(mat[1][0]);
1012 m[0][1]=-1.0/det*(mat[0][1]);
1013 m[1][1]= 1.0/det*(mat[0][0]);
1014 }
1015
1016 void METSignificance::AddMatrix(double (&X)[2][2],double (&Y)[2][2], double (&mat_new)[2][2]){
1017 mat_new[0][0]=X[0][0]+Y[0][0];
1018 mat_new[0][1]=X[0][1]+Y[0][1];
1019 mat_new[1][0]=X[1][0]+Y[1][0];
1020 mat_new[1][1]=X[1][1]+Y[1][1];
1021 }
1022
1023 void METSignificance::RotateXY(const double (&mat)[2][2], double (&mat_new)[2][2], double phi){
1024
1025 double c = std::cos(phi);
1026 double s = std::sin(phi);
1027 double cc = c*c;
1028 double ss = s*s;
1029 double cs = c*s;
1030
1031 double V11 = mat[0][0]*cc + mat[1][1]*ss - cs*(mat[1][0] + mat[0][1]);
1032 double V12 = mat[0][1]*cc - mat[1][0]*ss + cs*(mat[0][0] - mat[1][1]);
1033 double V21 = mat[1][0]*cc - mat[0][1]*ss + cs*(mat[0][0] - mat[1][1]);
1034 double V22 = mat[0][0]*ss + mat[1][1]*cc + cs*(mat[1][0] + mat[0][1]);
1035
1036 mat_new[0][0]=V11;
1037 mat_new[0][1]=V12;
1038 mat_new[1][0]=V21;
1039 mat_new[1][1]=V22;
1040 }
1041
1043 // Coefficients from Doug Schaefer <schae@cern.ch> and the MET subgroup
1044 double METSignificance::BiasPtSoftdir(const double PtSoft){
1045 if (PtSoft<60.) return (0.145)+(-0.45)*PtSoft;
1046 else return (0.145)+(-0.45)*(60.);
1047 }
1048
1049 // variation in ptsoft direction
1050 // Coefficients from Doug Schaefer <schae@cern.ch> and the MET subgroup
1051 double METSignificance::VarparPtSoftdir(const double PtSoft, const double SoftSumet){
1052 if (SoftSumet<25){
1053 if (PtSoft<50.) return 41.9+3.8*PtSoft+0.1*std::pow(PtSoft,2)-12.7+ 1.39*SoftSumet-0.03*std::pow(SoftSumet,2);
1054 else return 41.9+3.8*50.+0.1*std::pow(50.,2)-12.7+ 1.39*SoftSumet-0.03*std::pow(SoftSumet,2);
1055 }
1056 else{
1057 if (PtSoft<50.) return 41.9+3.8*PtSoft+0.1*std::pow(PtSoft,2);
1058 else return (40.5614)+(4.10965)*50.+(0.0955044)*std::pow(50.,2);
1059 }
1060 }
1061
1062 // Coefficients from Doug Schaefer <schae@cern.ch> and the MET subgroup
1063 double METSignificance::Var_Ptsoft(const double PtSoft){
1064 if (PtSoft<45.) return 40. + 2*PtSoft + 0.1*std::pow(PtSoft,2);
1065 else return 40. + 2*45 + 0.1*std::pow(45,2);
1066 }
1067
1068 // Coefficients from Doug Schaefer <schae@cern.ch> and the MET subgroup
1069 double METSignificance::Bias_PtSoftParall(const double PtSoft_Parall)
1070 {
1071 if (-60.<=PtSoft_Parall && PtSoft_Parall<0.) return -8. -0.4*PtSoft_Parall;
1072 if (-60.>PtSoft_Parall) return -8. -0.4 * (-60.);
1073 if( PtSoft_Parall>=0. && PtSoft_Parall<60.) return -8. -PtSoft_Parall;
1074 if(PtSoft_Parall>60.) return -8. -60.;
1075 return 0.0;
1076 }
1077
1078 void METSignificance::AddResoMap(const double varL, const double varT, const double CvLT, const int term){
1079
1080 m_term_VarL[term] += varL;
1081 m_term_VarT[term] += varT;
1082 m_term_CvLT[term] += CvLT;
1083
1084 }
1085
1086} //> end namespace met
#define ATH_CHECK
Evaluate an expression and check for errors.
#define ATH_MSG_DEBUG(x,...)
#define ATH_MSG_ERROR(x,...)
#define ATH_MSG_WARNING(x,...)
#define ATH_MSG_VERBOSE(x,...)
#define ATH_MSG_INFO(x,...)
Defines enum to access jet attribute and associated particles/objects.
static Double_t ss
std::string PathResolverFindCalibFile(const std::string &logical_file_name)
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T, V, H > &t)
const_iterator end() const noexcept
Return a const_iterator pointing past the end of the collection.
Helper class to provide constant type-safe access to aux data.
bool isAvailable(const ELT &e) const
Test to see if this variable exists in the store.
StatusCode setProperty(const std::string &name, const T &value)
set the given property
an object that can create a AsgTool
::StatusCode makePrivateTool(ToolHandle< T > &toolHandle) const
make a private tool with the given configuration
AsgTool(const std::string &name)
Constructor specifying the tool instance's name.
Definition AsgTool.cxx:58
virtual ~METSignificance()
Destructor:
double GetRho() const
void AddSoftTerm(const xAOD::MissingET *soft, const TVector3 &met_vect, double(&particle_sum)[2][2])
std::tuple< double, double, double > CovMatrixRotation(double var_x, double var_y, double cv_xy, double Phi)
double Var_Ptsoft(const double PtSoft)
ToolHandle< TauAnalysisTools::ITauResolutionTool > m_tauResolutionTool
double Significance_LT(double Numerator, double var_parall, double var_perpen, double cov)
StatusCode RotateToPhi(float phi)
double GetMETOverSqrtSumET() const
Gaudi::Property< std::string > m_tauTESConfig
Gaudi::Property< std::string > m_egESModel
double BiasPtSoftdir(const double PtSoft)
Parameterization with PtSoft Direction //.
void AddMatrix(double(&X)[2][2], double(&Y)[2][2], double(&mat_new)[2][2])
Gaudi::Property< bool > m_egUseFastsim
StatusCode initialize()
Dummy implementation of the initialisation function.
Gaudi::Property< std::string > m_jetCalibSeq
std::map< int, double > m_term_VarT
ToolHandle< IJetCalibrationTool > m_jetCalibTool
double GetSigDirectional() const
double GetSignificance() const
double GetVarL() const
static unsigned int getEtaBin(double jet_eta)
std::map< int, double > m_term_VarL
void AddResoMap(const double varL, const double varT, const double CvTV, const int term)
Gaudi::Property< std::string > m_egDecorrModel
Gaudi::Property< int > m_muonCalibMode
void RotateXY(const double(&mat)[2][2], double(&mat_new)[2][2], double phi)
StatusCode AddPhoton(const xAOD::IParticle *obj, float &pt_reso, float &phi_reso)
std::string m_configJetPhiResoFile
double GetPhiUnc(double jet_eta, double jet_phi, double jet_pt)
METSignificance()
Default constructor:
StatusCode SetLambda(const float px, const float py, const bool GeV=true)
StatusCode AddElectron(const xAOD::IParticle *obj, float &pt_reso, float &phi_reso, float avgmu)
double Bias_PtSoftParall(const double PtSoft_Parall)
Gaudi::Property< std::string > m_jetCalibArea
ToolHandle< CP::IEgammaCalibrationAndSmearingTool > m_egammaCalibTool
StatusCode AddTau(const xAOD::IParticle *obj, float &pt_reso, float &phi_reso)
StatusCode varianceMET(xAOD::MissingETContainer *metCont, float avgmu, const std::string &jetTermName, const std::string &softTermName, const std::string &totalMETName)
std::map< int, double > m_term_CvLT
Gaudi::Property< std::string > m_jetCalibConfig
void InvertMatrix(double(&mat)[2][2], double(&m)[2][2])
double VarparPtSoftdir(const double PtSoft, const double SoftSumet)
double GetVarT() const
double GetMETOverSqrtHT() const
StatusCode AddMuon(const xAOD::IParticle *obj, float &pt_reso, float &phi_reso, float avgmu)
ToolHandle< CP::IMuonCalibrationAndSmearingTool > m_muonCalibrationAndSmearingTool
double GetPUProb(double jet_eta, double jet_phi, double jet_pt, double jet_jvt, double jet_fjvt, float avgmu)
StatusCode AddJet(const xAOD::IParticle *obj, float &pt_reso, float &phi_reso, float &avgmu)
virtual double eta() const
The pseudorapidity ( ) of the particle.
virtual double pt() const override final
The transverse momentum ( ) of the particle.
Definition Egamma_v1.cxx:66
virtual double e() const override
The total energy of the particle.
Definition Egamma_v1.cxx:86
uint16_t author(uint16_t bitmask=EgammaParameters::AuthorALL) const
Get author.
virtual FourMom_t p4() const override final
The full 4-momentum of the particle as a TLoretzVector.
Definition Egamma_v1.cxx:94
const xAOD::CaloCluster * caloCluster(size_t index=0) const
Pointer to the xAOD::CaloCluster/s that define the electron candidate.
Class providing the definition of the 4-vector interface.
const_iterator find(const std::string &name) const
Find non-modifiable MET object by name.
float sumet() const
Returns.
float met() const
Returns .
float phi() const
Returns .
const std::string & name() const
Identifier getters.
float mpx() const
Returns .
float mpy() const
Returns .
uint64_t bitmask_t
Type for status word bit mask.
ElementLink< xAOD::IParticleContainer > iplink_t
static const SG::ConstAccessor< float > acc_covXY("covXY")
static const SG::ConstAccessor< float > acc_varY("varY")
static const SG::ConstAccessor< float > acc_varX("varX")
static const MissingETBase::Types::bitmask_t invisSource
Definition METHelpers.h:38
static const SG::ConstAccessor< float > acc_fjvt_der("DFCommonJets_fJvt")
static const SG::ConstAccessor< float > acc_jvt("Jvt")
static const SG::ConstAccessor< std::vector< iplink_t > > acc_constitObjLinks("ConstitObjectLinks")
static const SG::ConstAccessor< float > acc_fjvt("fJvt")
@ Jet
The object is a jet.
Definition ObjectType.h:40
@ Photon
The object is a photon.
Definition ObjectType.h:47
@ Muon
The object is a muon.
Definition ObjectType.h:48
@ TruthParticle
The object is a truth particle.
Definition ObjectType.h:67
@ Electron
The object is an electron.
Definition ObjectType.h:46
@ Tau
The object is a tau (jet).
Definition ObjectType.h:49
Jet_v1 Jet
Definition of the current "jet version".
MissingET_v1 MissingET
Version control by type defintion.
Definition MissingET.h:15
setRawEt setRawPhi int
Egamma_v1 Egamma
Definition of the current "egamma version".
Definition Egamma.h:17
TauJet_v3 TauJet
Definition of the current "tau version".
Definition TauJet.h:17
MissingETContainer_v1 MissingETContainer
float get_eta_calo(const xAOD::CaloCluster &cluster, int author, bool do_throw=false)
TruthParticle_v1 TruthParticle
Typedef to implementation.
Muon_v1 Muon
Reference the current persistent version:
Definition Muon.h:13
Electron_v1 Electron
Definition of the current "egamma version".
Definition Electron.h:17
static bool isSoftTerm(Types::bitmask_t bits, Region reg=Region::FullAcceptance)
@ Track
Indicator for MET contribution from reconstructed charged particle tracks.
static bool hasPattern(E bits, F mask)
Generic check for given pattern.
static bool isTotalTerm(Types::bitmask_t bits, Region reg=Region::FullAcceptance)