diff --git a/Reconstruction/MET/METInterface/METInterface/IMETSignificance.h b/Reconstruction/MET/METInterface/METInterface/IMETSignificance.h index 6639ed71c081a89403451c8c650b51c6f5c0701c..cd098155d366a235ff57c2c5717c6068e45c8c3e 100644 --- a/Reconstruction/MET/METInterface/METInterface/IMETSignificance.h +++ b/Reconstruction/MET/METInterface/METInterface/IMETSignificance.h @@ -29,6 +29,14 @@ namespace met { PthardParam=1, TSTParam=2 }; + enum METSignificanceResoTerms{ ResoNone=0, + ResoJet =1, + ResoSoft=2, + ResoEle =3, + ResoMuo =4, + ResoPho =5, + ResoTau =6 + }; }// end met namespace class IMETSignificance : virtual public asg::IAsgTool { @@ -45,20 +53,25 @@ public: // Version with single soft term virtual StatusCode varianceMET(xAOD::MissingETContainer* metCont, std::string jetTermName, std::string softTermName, std::string totalMETName)=0; + // rotates the phi direction of the object resolutions & recomputes the MET significance virtual StatusCode RotateToPhi(float phi) = 0; + // subtracks the vector lambda from the MET & recomputes the met signficance in new MET - lambda direction + virtual StatusCode SetLambda(const float px, const float py, const bool GeV=true) = 0; + /////////////////////////////////////////////////////////////////// // Additional utility commands /////////////////////////////////////////////////////////////////// - virtual double GetMETOverSqrtSumET() const = 0 ; - virtual double GetMETOverSqrtHT () const = 0 ; - virtual double GetSignificance() const = 0 ; - virtual double GetSigDirectional() const = 0 ; - virtual double GetRho() const = 0 ; - virtual double GetVarL() const = 0 ; - virtual double GetVarT() const = 0 ; - + virtual double GetMETOverSqrtSumET() const = 0 ; + virtual double GetMETOverSqrtHT () const = 0 ; + virtual double GetSignificance() const = 0 ; + virtual double GetSigDirectional() const = 0 ; + virtual double GetRho() const = 0 ; + virtual double GetVarL() const = 0 ; + virtual double GetVarT() const = 0 ; + virtual double GetTermVarL(const int term) const = 0 ; + virtual double GetTermVarT(const int term) const = 0 ; }; #endif diff --git a/Reconstruction/MET/METUtilities/METUtilities/METMaker.h b/Reconstruction/MET/METUtilities/METUtilities/METMaker.h index 1362351f20a0d9f093ff6e04c2a5a09783b57a01..916a4102159a6787ee3c00e46c82d035576dafeb 100644 --- a/Reconstruction/MET/METUtilities/METUtilities/METMaker.h +++ b/Reconstruction/MET/METUtilities/METUtilities/METMaker.h @@ -175,6 +175,7 @@ namespace met { bool m_muEloss; bool m_orCaloTaggedMuon; + bool m_greedyPhotons; ToolHandle<InDet::IInDetTrackSelectionTool> m_trkseltool; /// Default constructor: diff --git a/Reconstruction/MET/METUtilities/METUtilities/METSignificance.h b/Reconstruction/MET/METUtilities/METUtilities/METSignificance.h index 9d74b92bba130b5708362ed06732fdfc6885d929..892fad523a5a967162923d29517d50a1b161a636 100644 --- a/Reconstruction/MET/METUtilities/METUtilities/METSignificance.h +++ b/Reconstruction/MET/METUtilities/METUtilities/METSignificance.h @@ -75,7 +75,11 @@ namespace met { StatusCode varianceMET(xAOD::MissingETContainer* metCont, std::string jetTermName, std::string softTermName, std::string totalMETName); + // rotates the phi direction of the object resolutions & recomputes the MET significance StatusCode RotateToPhi(float phi); + + // subtracks the vector lambda from the MET & recomputes the met signficance in new MET - lambda direction + StatusCode SetLambda(const float px, const float py, const bool GeV=true); /////////////////////////////////////////////////////////////////// // Const methods: @@ -87,6 +91,8 @@ namespace met { double GetRho() const { return m_rho; } double GetVarL() const { return m_VarL; } double GetVarT() const { return m_VarT; } + double GetTermVarL(const int term) const { if(m_term_VarL.find(term)!=m_term_VarL.end()) return m_term_VarL.find(term)->second; return -1.0e3; } + double GetTermVarT(const int term) const { if(m_term_VarT.find(term)!=m_term_VarT.end()) return m_term_VarT.find(term)->second; return -1.0e3; } /////////////////////////////////////////////////////////////////// // Non-const methods: @@ -120,7 +126,6 @@ namespace met { double Significance_LT(double Numerator, double var_parall, double var_perpen, double cov); - void InvertMatrix(double (&mat)[2][2], double (&m)[2][2]); void AddMatrix(double (&X)[2][2],double (&Y)[2][2], double (&mat_new)[2][2]); void RotateXY(const double (&mat)[2][2], double (&mat_new)[2][2], double phi); @@ -133,13 +138,18 @@ namespace met { double Bias_PtSoftParall(const double PtSoft_Parall); double Var_Ptsoft(const double PtSoft); + // Fill Reso map + void AddResoMap(const double varL, + const double varT, + const double CvTV, const int term); + // variables double m_GeV; TVector3 m_met_vect; TVector3 m_soft_vect; TVector3 m_pthard_vect; - + TVector3 m_lamda_vect; int m_softTermParam; double m_softTermReso; @@ -162,6 +172,10 @@ namespace met { double m_met_VarT; double m_met_CvLT; + std::map<int, double> m_term_VarL; + std::map<int, double> m_term_VarT; + std::map<int, double> m_term_CvLT; + double m_met; double m_metphi; double m_metsoft; diff --git a/Reconstruction/MET/METUtilities/Root/METMaker.cxx b/Reconstruction/MET/METUtilities/Root/METMaker.cxx index b9182381b6957ecc384a0e3a04336eb72eb1a9e4..4fbe51d59dd1c556fbba3284ae082c63d8463386 100644 --- a/Reconstruction/MET/METUtilities/Root/METMaker.cxx +++ b/Reconstruction/MET/METUtilities/Root/METMaker.cxx @@ -114,6 +114,7 @@ namespace met { declareProperty("DoMuonEloss", m_muEloss = false ); declareProperty("ORCaloTaggedMuons", m_orCaloTaggedMuon = true ); + declareProperty("GreedyPhotons", m_greedyPhotons = false ); declareProperty("UseGhostMuons", m_useGhostMuons = false ); declareProperty("DoRemoveMuonJets", m_doRemoveMuonJets = true ); @@ -313,6 +314,17 @@ namespace met { ATH_MSG_VERBOSE(obj->type() << " (" << orig <<") with pt " << obj->pt() << " is " << ( selected ? "non-" : "") << "overlapping"); + // reinstate greedy photon option for testing + if (selected && obj->type() == xAOD::Type::Photon && m_greedyPhotons) { + for (size_t i = 0; i < assocs.size(); i++) { + std::vector<size_t> ind = assocs[i]->overlapIndices(orig); + std::vector<const xAOD::IParticle*> allObjects = assocs[i]->objects(); + for (size_t indi = 0; indi < ind.size(); indi++) { + if (allObjects[ind[indi]] && allObjects[ind[indi]]->type()==xAOD::Type::Electron) assocs[i]->setObjSelectionFlag(allObjects[ind[indi]],true); + } + } + } + //Do special overlap removal for calo tagged muons if(m_orCaloTaggedMuon && !removeOverlap && orig->type()==xAOD::Type::Muon && static_cast<const xAOD::Muon*>(orig)->muonType()==xAOD::Muon::CaloTagged) { for (size_t i = 0; i < assocs.size(); i++) { diff --git a/Reconstruction/MET/METUtilities/Root/METSignificance.cxx b/Reconstruction/MET/METUtilities/Root/METSignificance.cxx index 66020429ea296e0af87d1bb69f46ab7b45f916d6..3d682672ed69f4d9389a6e3e6a0d34d3508c3213 100644 --- a/Reconstruction/MET/METUtilities/Root/METSignificance.cxx +++ b/Reconstruction/MET/METUtilities/Root/METSignificance.cxx @@ -111,7 +111,7 @@ namespace met { // Property declaration // declareProperty("SoftTermParam", m_softTermParam = met::Random ); - declareProperty("SoftTermReso", m_softTermReso = 10.0 ); + declareProperty("SoftTermReso", m_softTermReso = 10.0 ); declareProperty("TreatPUJets", m_treatPUJets = true ); declareProperty("DoPhiReso", m_doPhiReso = false ); declareProperty("ApplyBias", m_applyBias = false ); @@ -205,6 +205,7 @@ namespace met { m_VarT = 0.0; m_CvLT = 0.0; + int metTerm = 0; double particle_sum[2][2] = {{0.0,0.0}, {0.0,0.0}}; m_metphi = 0.0; //Angle for rotation of the cov matrix @@ -213,6 +214,10 @@ namespace met { m_metsoftphi = 0.0; m_sumet=-1.0; m_ht=0.0; + m_term_VarL.clear(); + m_term_VarT.clear(); + m_term_CvLT.clear(); + unsigned nIterSoft=0; double softSumET=0.0; @@ -260,6 +265,7 @@ namespace met { AddSoftTerm(met, m_met_vect, particle_sum); m_metsoft = met->met()/m_GeV; m_metsoftphi = met->phi(); + metTerm = 2; // this is actually filled in AddSoftTerm // done with the soft term. go to the next term. continue; } @@ -269,19 +275,23 @@ namespace met { float pt_reso=0.0, phi_reso=0.0; if(obj->type()==xAOD::Type::Muon){ ATH_CHECK(AddMuon(obj, pt_reso, phi_reso)); + metTerm=4; }else if(obj->type()==xAOD::Type::Jet){ // make sure the container name matches if(met->name()!=jetTermName) continue; AddJet(obj, pt_reso, phi_reso); - + metTerm=1; }else if(obj->type()==xAOD::Type::Electron){ AddElectron(obj, pt_reso, phi_reso); + metTerm=3; }else if(obj->type()==xAOD::Type::Photon){ AddPhoton(obj, pt_reso, phi_reso); + metTerm=5; }else if(obj->type()==xAOD::Type::Tau){ AddTau(obj, pt_reso, phi_reso); + metTerm=6; } // compute NEW @@ -294,12 +304,17 @@ namespace met { m_VarT+=particle_u_rot[1][1]; m_CvLT+=particle_u_rot[0][1]; + // Save the resolutions separated for each object type + AddResoMap(particle_u_rot[0][0], + particle_u_rot[1][1], + particle_u_rot[0][1], + metTerm); + RotateXY (particle_u, particle_u_rot, obj->p4().Phi()); // positive phi rotation AddMatrix(particle_sum, particle_u_rot, particle_sum); // END compute NEW - ATH_MSG_VERBOSE("Resolution: " << pt_reso << " phi reso: " << phi_reso ); /* @@ -402,6 +417,32 @@ namespace met { return StatusCode::SUCCESS; } + // subtracks the vector lambda from the MET & recomputes the met signficance in new MET - lambda direction + StatusCode METSignificance::SetLambda(const float px, const float py, const bool GeV){ + + // compute the new direction + double GeVConv = GeV ? 1.0 : m_GeV; + m_lamda_vect.SetXYZ(px/GeVConv, py/GeVConv, 0.0); + m_lamda_vect = (m_met_vect - m_lamda_vect); + const double met_m_lamda = m_lamda_vect.Pt(); + + // Rotation (components) + std::tie(m_VarL, m_VarT, m_CvLT) = CovMatrixRotation(m_met_VarL , m_met_VarT, m_met_CvLT, (m_lamda_vect.Phi()-m_metphi)); + + + if( m_VarL != 0 ){ + m_significance = Significance_LT(met_m_lamda,m_VarL,m_VarT,m_CvLT ); + m_rho = m_CvLT / sqrt( m_VarL * m_VarT ) ; + } + ATH_MSG_DEBUG(" Significance (squared) at new phi: " << m_significance + << " rho: " << GetRho() + << " MET: " << m_met + << " sigmaL: " << GetVarL() + << " sigmaT: " << GetVarT() ); + + return StatusCode::SUCCESS; + } + /////////////////////////////////////////////////////////////////// // Private methods /////////////////////////////////////////////////////////////////// @@ -522,6 +563,12 @@ namespace met { m_VarL+=particle_u_rot[0][0]; m_VarT+=particle_u_rot[1][1]; m_CvLT+=particle_u_rot[0][1]; + + // Save the resolutions separated for each object type + AddResoMap(particle_u_rot[0][0], + particle_u_rot[1][1], + particle_u_rot[0][1], + met::ResoSoft); RotateXY (particle_u, particle_u_rot,-1.0*soft->phi()); // negative phi rotation AddMatrix(particle_sum, particle_u_rot, particle_sum); @@ -550,6 +597,12 @@ namespace met { m_VarT+=particle_u_rot[1][1]; m_CvLT+=particle_u_rot[0][1]; + // Save the resolutions separated for each object type + AddResoMap(particle_u_rot[0][0], + particle_u_rot[1][1], + particle_u_rot[0][1], + met::ResoSoft); + RotateXY (particle_u, particle_u_rot,-1.0*m_pthard_vect.Phi()); // negative phi rotation AddMatrix(particle_sum, particle_u_rot, particle_sum); @@ -571,6 +624,12 @@ namespace met { m_VarT+=particle_u_rot[1][1]; m_CvLT+=particle_u_rot[0][1]; + // Save the resolutions separated for each object type + AddResoMap(particle_u_rot[0][0], + particle_u_rot[1][1], + particle_u_rot[0][1], + met::ResoSoft); + RotateXY (particle_u, particle_u_rot,-1.0*soft->phi()); // negative phi rotation AddMatrix(particle_sum, particle_u_rot, particle_sum); @@ -809,5 +868,20 @@ namespace met { return 0.0; } + void METSignificance::AddResoMap(const double varL, + const double varT, + const double CvLT, const int term){ + if(m_term_VarL.find(term)==m_term_VarL.end()){ + m_term_VarL[term] = 0.0; + m_term_VarT[term] = 0.0; + m_term_CvLT[term] = 0.0; + } + + m_term_VarL[term] += varL; + m_term_VarT[term] += varT; + m_term_CvLT[term] += CvLT; + + } + } //> end namespace met diff --git a/Reconstruction/MET/METUtilities/util/example_METSignificance.cxx b/Reconstruction/MET/METUtilities/util/example_METSignificance.cxx index 6d06be50a31d82cfb70bdbb6770747b810ac49e9..7931dc2d4ea0bc1f92ac0bbe4c8c7d989f374214 100644 --- a/Reconstruction/MET/METUtilities/util/example_METSignificance.cxx +++ b/Reconstruction/MET/METUtilities/util/example_METSignificance.cxx @@ -287,11 +287,31 @@ int main( int argc, char* argv[] ){std::cout << __PRETTY_FUNCTION__ << std::endl std::cout << "Soft term: " << static_cast<xAOD::MissingET*>(*(newMetContainer->find("PVSoftTrk")))->met() << " phi: " << static_cast<xAOD::MissingET*>(*(newMetContainer->find("PVSoftTrk")))->phi() << std::endl; + + // Print the METSignificance terms (e.g. jet, muon, ele, pho, etc) + // NOTE::: these are not currently rotated if the MET is rotated + std::cout << " jet VarL: " << metSignif->GetTermVarL(met::ResoJet) << " GeV VarT: " << metSignif->GetTermVarT(met::ResoJet) << " GeV" << std::endl; + std::cout << " muon VarL: " << metSignif->GetTermVarL(met::ResoMuo) << " GeV VarT: " << metSignif->GetTermVarT(met::ResoMuo) << " GeV" << std::endl; + std::cout << " electron VarL: " << metSignif->GetTermVarL(met::ResoEle) << " GeV VarT: " << metSignif->GetTermVarT(met::ResoEle) << " GeV" << std::endl; + std::cout << " photon VarL: " << metSignif->GetTermVarL(met::ResoPho) << " GeV VarT: " << metSignif->GetTermVarT(met::ResoPho) << " GeV" << std::endl; + std::cout << " tau VarL: " << metSignif->GetTermVarL(met::ResoTau) << " GeV VarT: " << metSignif->GetTermVarT(met::ResoTau) << " GeV" << std::endl; + std::cout << " Soft term VarL: " << metSignif->GetTermVarL(met::ResoSoft) << " GeV VarT: " << metSignif->GetTermVarT(met::ResoSoft) << " GeV" << std::endl; + std::cout << " other/bug VarL: " << metSignif->GetTermVarL(met::ResoNone) << " GeV VarT: " << metSignif->GetTermVarT(met::ResoNone) << " GeV" << std::endl; + } // extracting the MET significance std::cout << "MET significance: " << metSignif->GetSignificance() << std::endl; + if(debug){ + // Try a rotation to a new lambda parameter + std::cout << " Lambda Test Before: " << metSignif->GetSignificance() << " VarL: " << metSignif->GetVarL() << " VarT: " << metSignif->GetVarT() << std::endl; + metSignif->SetLambda(0.0, 0.0); + std::cout << " Lambda Test 0: " << metSignif->GetSignificance() << " VarL: " << metSignif->GetVarL() << " VarT: " << metSignif->GetVarT() << std::endl; + metSignif->SetLambda(10.0, 10.0); + std::cout << " Lambda Test 10: " << metSignif->GetSignificance() << " VarL: " << metSignif->GetVarL() << " VarT: " << metSignif->GetVarT() << std::endl; + } + ANA_CHECK(store->record( newMetContainer, "FinalMETContainer" )); ANA_CHECK(store->record( newMetAuxContainer, "FinalMETContainerAux."));