Skip to content
Snippets Groups Projects
Commit ff045d2f authored by Walter Lampl's avatar Walter Lampl
Browse files

Merge branch 'CandidateMatchHelpers_add_r_for_correct_eta' into 'master' 

CandidateMatchHelpers: Try to transform the cluster calo position vector wrt to (0,0,0) to a position vector wrt to track perigee position

See merge request atlas/athena!26703
parents 18ffb373 80d9a8e2
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 138 additions and 80 deletions
Reconstruction/egamma/egammaAlgs/src/egammaSelectedTrackCopy.cxx
+ 42
− 22
View file @ ff045d2f
......@@ -192,14 +192,17 @@ bool egammaSelectedTrackCopy::Select(const EventContext& ctx,
{
ATH_MSG_DEBUG("egammaSelectedTrackCopy::Select()" );
if (cluster == 0 || track == 0) return false;
if (cluster == nullptr || track == nullptr) return false;
const Trk::Perigee& candidatePerigee = track->perigeeParameters();
//Get Perigee Parameters
const double trkPhi = candidatePerigee.parameters()[Trk::phi];
const double trkEta = candidatePerigee.eta();
const double r_first=candidatePerigee.position().perp();
const double z_first=candidatePerigee.position().z();
const double z_perigee=candidatePerigee.position().z();
const double r_perigee=candidatePerigee.position().perp();
const Amg::Vector3D PerigeeXYZPosition(candidatePerigee.position().x(),
candidatePerigee.position().y(),
z_perigee);
//Get Cluster parameters
const double clusterEta=cluster->etaBE(2);
......@@ -215,42 +218,59 @@ bool egammaSelectedTrackCopy::Select(const EventContext& ctx,
return false;
}
//Calculate corrrected eta and Phi
const double etaclus_corrected = CandidateMatchHelpers::CorrectedEta(clusterEta,z_first,isEndCap);
const double phiRot = CandidateMatchHelpers::PhiROT(Et,trkEta, track->charge(),r_first ,isEndCap) ;
const double phiRotTrack = CandidateMatchHelpers::PhiROT(track->pt(),trkEta, track->charge(),r_first ,isEndCap) ;
//Calcualate deltaPhis
const double deltaPhiStd = P4Helpers::deltaPhi(cluster->phiBE(2), trkPhi);
const double trkPhiCorr = P4Helpers::deltaPhi(trkPhi, phiRot);
const double deltaPhi2 = P4Helpers::deltaPhi(cluster->phiBE(2), trkPhiCorr);
//Calculate the eta/phi of the cluster as would be seen from the perigee position of the Track
const Amg::Vector3D XYZClusterWrtTrackPerigee = CandidateMatchHelpers::approxXYZwrtPoint(*cluster,
PerigeeXYZPosition,
isEndCap);
//Calculate the possible rotation of the track
//Once assuming the cluster Et being the better estimate (e.g big brem)
const double phiRotRescaled = CandidateMatchHelpers::PhiROT(Et,trkEta,
track->charge(),
r_perigee,
isEndCap) ;
//And also assuming the track Pt being correct
const double phiRotTrack = CandidateMatchHelpers::PhiROT(track->pt(),
trkEta,
track->charge(),
r_perigee,
isEndCap) ;
//
const double clusterPhiCorrected=XYZClusterWrtTrackPerigee.phi();
//deltaPhi between the track and the cluster
const double deltaPhiStd = P4Helpers::deltaPhi(clusterPhiCorrected, trkPhi);
//deltaPhi between the track and the cluster accounting for rotation assuming cluster Et is a better estimator
const double trkPhiRescaled= P4Helpers::deltaPhi(trkPhi, phiRotRescaled);
const double deltaPhiRescaled = P4Helpers::deltaPhi(clusterPhiCorrected, trkPhiRescaled);
//deltaPhi between the track and the cluster accounting for rotation
const double trkPhiCorrTrack = P4Helpers::deltaPhi(trkPhi, phiRotTrack);
const double deltaPhi2Track = P4Helpers::deltaPhi(cluster->phiBE(2), trkPhiCorrTrack);
const double deltaPhiTrack = P4Helpers::deltaPhi(clusterPhiCorrected, trkPhiCorrTrack);
/*
* First we will see if it fails the quick match
* Then if it passed it will get 2 chances to be selected
* One if it matches from last measurement
* The second if it matched from Perigee rescales
*/
//Broad phi check
if ( (fabs(deltaPhi2) > 2*m_broadDeltaPhi)
&& (fabs(deltaPhi2Track) > 2.*m_broadDeltaPhi)
if ((fabs(deltaPhiRescaled) > 2*m_broadDeltaPhi)
&& (fabs(deltaPhiTrack) > 2.*m_broadDeltaPhi)
&& (fabs(deltaPhiStd) > 2*m_broadDeltaPhi)){
ATH_MSG_DEBUG("FAILS broad window phi match (track phi, phirotCluster , phiRotTrack ,cluster phi): ( "
<< trkPhi << ", " << phiRot<< ", "<<phiRotTrack<< ", " << cluster->phiBE(2) << ")" );
ATH_MSG_DEBUG("FAILS broad window phi match (track phi, phirotCluster , phiRotTrack , "
<<"cluster phi corrected, cluster phi): ( "
<< trkPhi << ", " << phiRotRescaled<< ", "<<phiRotTrack<< ", "
<<clusterPhiCorrected<< ", " <<cluster->phi() <<")" );
return false;
}
//if TRT we can stop here , we can not check much in eta really.
if(trkTRT){
return true;
}
const double clusterEtaCorrected=XYZClusterWrtTrackPerigee.eta();
//eta check
if (fabs(cluster->etaBE(2) - trkEta) > 2*m_broadDeltaEta &&
fabs( etaclus_corrected- trkEta) > 2.*m_broadDeltaEta){
//Broad eta check
if (fabs(cluster->eta() - trkEta) > 2*m_broadDeltaEta &&
fabs(clusterEtaCorrected- trkEta) > 2.*m_broadDeltaEta){
ATH_MSG_DEBUG("FAILS broad window eta match (track eta, cluster eta, cluster eta corrected): ( "
<< trkEta << ", " << cluster->etaBE(2) <<", "<<etaclus_corrected<<")" );
<< trkEta << ", " << cluster->etaBE(2) <<", "<<clusterEtaCorrected<<" )");
return false;
}
......
Reconstruction/egamma/egammaTools/src/EMTrackMatchBuilder.cxx
+ 44
− 22
View file @ ff045d2f
......@@ -380,43 +380,65 @@ EMTrackMatchBuilder::isCandidateMatch(const xAOD::CaloCluster* cluster,
//Decide whether to try the opposite direction (cosmics)
const double trkPhi = (!flip) ? candidatePerigee.parameters()[Trk::phi] : -candidatePerigee.parameters()[Trk::phi];
const double trkEta = (!flip) ? candidatePerigee.eta() : -candidatePerigee.eta();
const double r_first=candidatePerigee.position().perp();
const double z_first=candidatePerigee.position().z();
const double z_perigee=candidatePerigee.position().z();
const double r_perigee=candidatePerigee.position().perp();
const Amg::Vector3D PerigeeXYZPosition(candidatePerigee.position().x(),
candidatePerigee.position().y(),
z_perigee);
//Cluster variables
const double clusterEta=cluster->etaBE(2);
const double clusterEta=cluster->eta();
const bool isEndCap= !xAOD::EgammaHelpers::isBarrel(cluster);
const double Et= cluster->e()/cosh(trkEta);
const double clusterPhi=cluster->phiBE(2);
const double clusterPhi=cluster->phi();
//Avoid clusters with |eta| > 10 or Et less than 10 MeV
if(fabs(clusterEta)>10.0 || Et <10){
return false;
}
const double etaclus_corrected = CandidateMatchHelpers::CorrectedEta(clusterEta,z_first,isEndCap);
const double phiRot = CandidateMatchHelpers::PhiROT(Et,trkEta, track->charge(),r_first ,isEndCap) ;
const double phiRotTrack = CandidateMatchHelpers::PhiROT(track->pt(),trkEta, track->charge(),r_first ,isEndCap) ;
const double deltaPhiStd = P4Helpers::deltaPhi(clusterPhi, trkPhi);
const double trkPhiCorr = P4Helpers::deltaPhi(trkPhi, phiRot);
const double deltaPhi2 = P4Helpers::deltaPhi(clusterPhi, trkPhiCorr);
const double trkPhiCorrTrack = P4Helpers::deltaPhi(trkPhi, phiRotTrack);
const double deltaPhi2Track = P4Helpers::deltaPhi(clusterPhi, trkPhiCorrTrack);
//Calculate the eta/phi of the cluster as would be seen from the perigee position of the Track
const Amg::Vector3D XYZClusterWrtTrackPerigee = CandidateMatchHelpers::approxXYZwrtPoint(*cluster,
PerigeeXYZPosition,
isEndCap);
const double clusterEtaCorrected=XYZClusterWrtTrackPerigee.eta();
//check eta match . Both metrics need to fail in order to disgard the track
if ( (fabs(clusterEta - trkEta) > 2.*m_broadDeltaEta) &&
(fabs( etaclus_corrected- trkEta) > 2.*m_broadDeltaEta)){
(fabs(clusterEtaCorrected- trkEta) > 2.*m_broadDeltaEta)){
ATH_MSG_DEBUG(" Fails broad window eta match (track eta, cluster eta, cluster eta corrected): ( "
<< trkEta << ", " << clusterEta <<", "<<etaclus_corrected<<")" );
<< trkEta << ", " << clusterEta <<", "<<clusterEtaCorrected<<")" );
return false;
}
//Calculate the possible rotation of the track
//Once assuming the cluster Et being the better estimate (e.g big brem)
const double phiRotRescaled = CandidateMatchHelpers::PhiROT(Et,trkEta,
track->charge(),
r_perigee,
isEndCap) ;
//And also assuming the track Pt being correct
const double phiRotTrack = CandidateMatchHelpers::PhiROT(track->pt(),
trkEta,
track->charge(),
r_perigee,
isEndCap) ;
//
const double clusterPhiCorrected=XYZClusterWrtTrackPerigee.phi();
//deltaPhi between the track and the cluster
const double deltaPhiStd = P4Helpers::deltaPhi(clusterPhiCorrected, trkPhi);
//deltaPhi between the track and the cluster accounting for rotation assuming cluster Et is a better estimator
const double trkPhiRescaled= P4Helpers::deltaPhi(trkPhi, phiRotRescaled);
const double deltaPhiRescaled = P4Helpers::deltaPhi(clusterPhiCorrected, trkPhiRescaled);
//deltaPhi between the track and the cluster accounting for rotation
const double trkPhiCorrTrack = P4Helpers::deltaPhi(trkPhi, phiRotTrack);
const double deltaPhiTrack = P4Helpers::deltaPhi(clusterPhiCorrected, trkPhiCorrTrack);
//It has to fail all phi metrics in order to be disgarded
if ( (fabs(deltaPhi2) > 2.*m_broadDeltaPhi) &&
(fabs(deltaPhi2Track) > 2.*m_broadDeltaPhi) &&
if ( (fabs(deltaPhiRescaled) > 2.*m_broadDeltaPhi) &&
(fabs(deltaPhiTrack) > 2.*m_broadDeltaPhi) &&
(fabs(deltaPhiStd) > 2.*m_broadDeltaPhi) ){
ATH_MSG_DEBUG(" Fails broad window eta match (track eta, cluster eta, cluster eta corrected): ( "
<< trkEta << ", " << cluster->etaBE(2) <<", "<<etaclus_corrected<<")" );
ATH_MSG_DEBUG("FAILS broad window phi match (track phi, phirotCluster , phiRotTrack , "
<<"cluster phi corrected, cluster phi): ( "
<< trkPhi << ", " << phiRotRescaled<< ", "<<phiRotTrack<< ", "
<<clusterPhiCorrected<< ", " <<clusterPhi <<")" );
return false;
}
//if not false returned we end up here
......
Reconstruction/egamma/egammaUtils/CMakeLists.txt
+ 2
− 1
View file @ ff045d2f
......@@ -11,6 +11,7 @@ atlas_depends_on_subdirs( PUBLIC
Event/xAOD/xAODCaloEvent
Event/xAOD/xAODTracking
Event/xAOD/xAODEgamma
DetectorDescription/GeoPrimitives
PRIVATE
Event/FourMomUtils
PhysicsAnalysis/AnalysisCommon/AnalysisUtils
......@@ -25,7 +26,7 @@ atlas_add_library( egammaUtils
PUBLIC_HEADERS egammaUtils
INCLUDE_DIRS ${ROOT_INCLUDE_DIRS}
PRIVATE_INCLUDE_DIRS
LINK_LIBRARIES ${ROOT_LIBRARIES} AsgTools xAODCaloEvent xAODTracking xAODEgamma
LINK_LIBRARIES ${ROOT_LIBRARIES} AsgTools xAODCaloEvent xAODTracking xAODEgamma GeoPrimitives
PRIVATE_LINK_LIBRARIES FourMomUtils PathResolver AnalysisUtilsLib)
atlas_add_dictionary( egammaUtilsDict
......
Reconstruction/egamma/egammaUtils/Root/CandidateMatchHelpers.cxx
+ 33
− 29
View file @ ff045d2f
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
*/
#include "egammaUtils/CandidateMatchHelpers.h"
#include "xAODCaloEvent/CaloCluster.h"
#include <cmath>
double CandidateMatchHelpers::CorrectedEta(const double clusterEta,
const double z_first,
const bool isEndCap){
if (clusterEta == -999){
return clusterEta;
}
//===========================================================//
Amg::Vector3D CandidateMatchHelpers::approxXYZwrtPoint (const xAOD::CaloCluster& cluster,
const Amg::Vector3D& point,
const bool isEndCap){
return approxXYZwrtATLAS(cluster,isEndCap) - point;
}
Amg::Vector3D CandidateMatchHelpers::approxXYZwrtATLAS (const xAOD::CaloCluster& cluster,
const bool isEndCap){
//These are at the face of the calorimeter
const double Rcalo = 1500;
const double Zcalo = 3700;
double Rclus=Rcalo;
double Zclus=Zcalo;
const double RfaceCalo = 1500;
const double ZfaceCalo = 3700;
//Get the Rclus , Zclus given the cluster eta
double Rclus=RfaceCalo;
double Zclus=ZfaceCalo;
const double clusterEta=cluster.eta();
if(clusterEta!=0){
//eta !=0
double tanthetaclus=(1/cosh(clusterEta))/tanh(clusterEta);
/*
* tahn(eta) == cos(theta)
* 1/cosh(clusterEta) == sin(theta)
* tan =sin/cos
*/
const double tanthetaclus=(1/cosh(clusterEta))/tanh(clusterEta);
if (isEndCap) {
Rclus=fabs(Zcalo*(tanthetaclus));
Rclus=fabs(ZfaceCalo*(tanthetaclus));
//Negative Eta ---> negative Z
if(clusterEta<0){
Zclus = -Zclus;
......@@ -28,22 +38,16 @@ double CandidateMatchHelpers::CorrectedEta(const double clusterEta,
}
else{
if(tanthetaclus!=0){
Zclus=Rcalo/(tanthetaclus);
Zclus=RfaceCalo/(tanthetaclus);
}
}
}
else{
//eta 0
else{//when eta ==0
Zclus=0;
}
//correct for a possible shift in the track perigee wrt to (0,0)
double thetaclus_corrected=atan2(Rclus,Zclus-z_first);
double etaclus_corrected = 99999; //if theta =0 or M_PI the eta +-inf this happens when Rclus =0
if(Rclus!=0){
etaclus_corrected = -log(tan(thetaclus_corrected/2));
}
return etaclus_corrected;
const double clusterPhi=cluster.phi();
return Amg::Vector3D (Rclus*cos(clusterPhi),Rclus*sin(clusterPhi),Zclus);
}
double CandidateMatchHelpers::PhiROT(const double Et,
......@@ -58,12 +62,12 @@ double CandidateMatchHelpers::PhiROT(const double Et,
double phiRot = 0.0;
double ecCorr = 1.0;
if (isEndCap) {
double ecFactor = 1.0/(0.8*atan(15.0/34.0));
double sinTheta0 = 2.0*atan(exp(-fabs(Eta)));
const double ecFactor = 1.0/(0.8*atan(15.0/34.0));
const double sinTheta0 = 2.0*atan(exp(-fabs(Eta)));
ecCorr = sinTheta0*sqrt(sinTheta0)*ecFactor;
}
////
double Rscaled =(Rcalo-r_first)*(1./Rcalo);
const double Rscaled =(Rcalo-r_first)*(1./Rcalo);
phiRot = Rscaled*ecCorr*charge*430./(Et);
return phiRot;
}
......
Reconstruction/egamma/egammaUtils/egammaUtils/CandidateMatchHelpers.h
+ 17
− 6
View file @ ff045d2f
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
*/
#ifndef CANDIDATEMATCHHELPERS_H
#define CANDIDATEMATCHHELPERS_H
#include <cmath>
#include <iostream>
#include "GeoPrimitives/GeoPrimitives.h"
#include "xAODCaloEvent/CaloClusterFwd.h"
namespace CandidateMatchHelpers{
double CorrectedEta(const double clusterEta,const double z_first,const bool isEndCap);
double PhiROT(const double Et,const double Eta, const int charge, const double r_first ,const bool isEndCap);
/**Function to get the (x,y,z) of the cluster wrt to a point (x0,y0,z0)*/
Amg::Vector3D approxXYZwrtPoint (const xAOD::CaloCluster& cluster,
const Amg::Vector3D& point,
const bool isEndCap);
/**Function to get the (x,y,z) of the cluster in the ATLAS frame origin (0,0,0)*/
Amg::Vector3D approxXYZwrtATLAS (const xAOD::CaloCluster& clusterEta,
const bool isEndCap);
/** Function to calculate the approximate rotation in phi/bending of a track until it reaches the calo*/
double PhiROT(const double pt,const double eta, const int charge,
const double r_start,const bool isEndCap);
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment