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