diff --git a/Reconstruction/egamma/egammaTrackTools/python/egammaTrackToolsFactories.py b/Reconstruction/egamma/egammaTrackTools/python/egammaTrackToolsFactories.py index ca342bd3284ecc0d559d9b955e44c5ebb00d363b..97e1779bd76053acc5cb519c66a834ee1d93b493 100644 --- a/Reconstruction/egamma/egammaTrackTools/python/egammaTrackToolsFactories.py +++ b/Reconstruction/egamma/egammaTrackTools/python/egammaTrackToolsFactories.py @@ -15,16 +15,9 @@ ExtrapolateToCaloTool = ToolFactory( TrackToCaloConf.ExtrapolateToCaloTool, name = 'exToCalo', Extrapolator = egammaExtrapolator) -# Utility for conversion tool -from TrkVertexFitterUtils.TrkVertexFitterUtilsConf import Trk__NeutralParticleParameterCalculator -convUtils = ToolFactory( Trk__NeutralParticleParameterCalculator, - name = "convUtils", - LinearizedTrackFactory = None) - EMExtrapolationTools = ToolFactory( egammaTrackToolsConf.EMExtrapolationTools, ExtrapolateToCaloTool = ExtrapolateToCaloTool, - ConversionUtils = convUtils, doPrint = True) #--------------------------------------- diff --git a/Reconstruction/egamma/egammaTrackTools/src/EMExtrapolationTools.cxx b/Reconstruction/egamma/egammaTrackTools/src/EMExtrapolationTools.cxx index a87dd34db2d7c0080e6a4b34671886e48211fa23..e10879f07bebb81148d54d2f97fe0437421f6edb 100644 --- a/Reconstruction/egamma/egammaTrackTools/src/EMExtrapolationTools.cxx +++ b/Reconstruction/egamma/egammaTrackTools/src/EMExtrapolationTools.cxx @@ -10,16 +10,13 @@ #include "RecoToolInterfaces/IExtrapolateToCaloTool.h" #include "CaloDetDescr/CaloDepthTool.h" #include "InDetIdentifier/TRT_ID.h" -#include "TrkVertexFitterInterfaces/INeutralParticleParameterCalculator.h" #include "TrkNeutralParameters/NeutralParameters.h" -// #include "xAODTracking/TrackingPrimitives.h" #include "xAODCaloEvent/CaloCluster.h" #include "xAODTracking/TrackParticle.h" #include "xAODTracking/Vertex.h" - #include "TrkTrack/Track.h" -#include "TrkPseudoMeasurementOnTrack/TrkPseudoMeasurementOnTrack/PseudoMeasurementOnTrack.h" +#include "TrkPseudoMeasurementOnTrack/PseudoMeasurementOnTrack.h" #include "xAODEgamma/EgammaxAODHelpers.h" @@ -28,14 +25,10 @@ EMExtrapolationTools::EMExtrapolationTools(const std::string& type, const std::string& name, const IInterface* parent) : - AthAlgTool(type, name, parent), - m_convUtils("Trk::NeutralParticleParameterCalculator") -{ + AthAlgTool(type, name, parent){ declareProperty( "ExtrapolateToCaloTool", m_extrapolateToCalo); - // Name of the utility for conversion - declareProperty("ConversionUtils", m_convUtils, "Handle of the utility for conversion"); declareProperty( "BroadDeltaEta", m_broadDeltaEta = 0.05); declareProperty( "BroadDeltaPhi", m_broadDeltaPhi = 0.10); @@ -98,12 +91,6 @@ StatusCode EMExtrapolationTools::initialize() ATH_MSG_DEBUG("Could not get TRT_ID helper !"); } if (!m_trtId) m_guessTRTsection = true; - - if (m_convUtils.retrieve().isFailure()) { - ATH_MSG_ERROR("Could not find ConvUtils tool."); - return StatusCode::FAILURE; - } - else ATH_MSG_DEBUG("ConvUtils retrieved"); return StatusCode::SUCCESS; } @@ -114,8 +101,12 @@ StatusCode EMExtrapolationTools::finalize() return StatusCode::SUCCESS; } +// ======================= Electron/trackParticle Related Interfaces + ////////////////////////////////////////////////////////////////////////////// //* Extrapolate TrackParticle to electron CaloCluster +//* Returns also selection +//* i.e performs the cuts ////////////////////////////////////////////////////////////////////////////// bool EMExtrapolationTools::matchesAtCalo(const xAOD::CaloCluster* cluster, @@ -144,68 +135,11 @@ EMExtrapolationTools::matchesAtCalo(const xAOD::CaloCluster* cluster, extrapFrom) ); } -/////////////////////////////////////////////////////////////////////////////// -//* Extrapolate double-track conversion with the neutral perigee -/////////////////////////////////////////////////////////////////////////////// -bool -EMExtrapolationTools::matchesAtCalo(const xAOD::CaloCluster* cluster, - const Trk::NeutralParameters* phPerigee, - bool isTRT, - double& deltaEta, - double& deltaPhi) -{ - - double tmpDeltaEta = -999.; - double tmpDeltaPhi = -999.; - int trapped = 99; - - const Trk::NeutralParameters* result = 0; - double offset = 0.; - CaloCell_ID::CaloSample sample; - - - if ( xAOD::EgammaHelpers::isBarrel( cluster ) ) { - sample = CaloCell_ID::EMB2; - result = m_extrapolateToCalo->extrapolate((*phPerigee), sample, offset, Trk::photon, Trk::alongMomentum); - if (result) { - trapped = 0; - tmpDeltaEta = cluster->etaSample(sample) - result->position().eta(); - tmpDeltaPhi = m_phiHelper.diff(cluster->phiSample(sample),result->position().phi()); - } - }else { - sample = CaloCell_ID::EME2; - result = m_extrapolateToCalo->extrapolate((*phPerigee), sample, offset, Trk::photon, Trk::alongMomentum); - if (result) { - trapped = 0; - tmpDeltaEta = cluster->etaSample(sample) - result->position().eta(); - tmpDeltaPhi = m_phiHelper.diff(cluster->phiSample(sample),result->position().phi()); - } - } - if(result) delete result; - - deltaEta = tmpDeltaEta; - deltaPhi = tmpDeltaPhi; - - // Match in a broad eta / phi window. For TRT-only tracks we don't know if they are in the - // barrel or in the endcap. Apply the loosest cut (the code below would be the correct one) - // || (abs(isTRTB)==1 && fabs(deltaEta) > m_TRTbarrelDeltaEta) || - // || (abs(isTRTB)==2 && fabs(deltaEta) > m_TRTendcapDeltaEta) || - - if (trapped != 0 - || (!isTRT && (fabs(deltaEta) > m_broadDeltaEta) ) - || (isTRT && (fabs(deltaEta) > std::max(m_TRTbarrelDeltaEta, m_TRTendcapDeltaEta)) ) - || fabs(deltaPhi) > m_broadDeltaPhi) return false; - - if ((isTRT && fabs(deltaPhi) < std::max(m_narrowDeltaPhiTRTbarrel, m_narrowDeltaPhiTRTendcap)) || - (!isTRT && fabs(deltaEta)< m_narrowDeltaEta && fabs(deltaPhi) < m_narrowDeltaPhi) ) - return true; - - return false; -} - ////////////////////////////////////////////////////////////////////////////// //* Extrapolate TrackParticle to electron CaloCluster +//* Returns also selection +//* i.e performs the cuts ////////////////////////////////////////////////////////////////////////////// bool EMExtrapolationTools::matchesAtCalo(const xAOD::CaloCluster* cluster, @@ -241,149 +175,53 @@ EMExtrapolationTools::matchesAtCalo(const xAOD::CaloCluster* cluster, << " deltaPhi: " << deltaPhi.size()); return false; } - - // local variables - int trapped = 99; // Sample for testing match unsigned int iSampling = 2; - // result of extrapolation - const Trk::TrackParameters* result = 0; - // offset for extrapolation - double offset = 0.; - - - const Trk::TrackParameters* trkPar =0; - Trk::CurvilinearParameters temp; - if(isTRT || fromLastMeasurement == extrapFrom ) { - //For TRT it is always the last - int index(-1); - for(unsigned int i(0); i< trkPB->numberOfParameters() ; ++i ){ - if( xAOD::LastMeasurement == trkPB->parameterPosition(i) ){ - index = i; - break; - } - } - if(index!=-1){ - temp = trkPB->curvilinearParameters(index); - trkPar = &temp; - } else { - ATH_MSG_WARNING("Track Particle does not contain Last Measurement track parameters"); - } - } - else if (fromPerigee == extrapFrom) { - // From perigee - trkPar = &trkPB->perigeeParameters(); - } - else if (fromPerigeeRescaled == extrapFrom) { - // From rescaled perigee - trkPar = getRescaledPerigee(trkPB, cluster); - } - - if(!trkPar){ - ATH_MSG_WARNING("Track Particle does not have the requested track parameters -- extrapolation not possible"); - return false; - } - - // Should we flip the sign for deltaPhi? - bool flipSign = false; - if(trkPar->charge() > 0) flipSign = true; - - // In case of standalone TRT tracks get an idea where the last measurement is int isTRTB = 0; // barrel or endcap TRT if(isTRT){ + // Decide matching for TRT standalone + // In case of standalone TRT tracks get an idea where the last measurement is if (!m_trtId) { ATH_MSG_WARNING("Should have m_trtId defined for isTRT"); return false; } // Get last measurement for TRT check - const Trk::TrackParameters* trkParTRT = trkPar; + + Trk::CurvilinearParameters temp = getLastMeasurement(trkPB); + const Trk::TrackParameters* trkParTRT = &temp; + + if(!trkParTRT){ + ATH_MSG_ERROR("MatchAtCalo TRT: Cannot access track parameters"); + return StatusCode::FAILURE; + } + const Trk::Surface& sf = trkParTRT->associatedSurface(); const Identifier tid = sf.associatedDetectorElementIdentifier(); isTRTB = m_trtId->barrel_ec(tid); // First pass on TRT tracks, skip barrel tracks matching endcap clusters and vice-versa - if((isTRTB==2 && (cluster->eta()<=0.6 || cluster->eta()>=2.4) ) || + if((isTRTB==2 && (cluster->eta()<=0.6 || cluster->eta()>=2.4) ) || (isTRTB==-2 && (cluster->eta()>=-0.6 || cluster->eta()<=-2.4) ) || (isTRTB==1 && (cluster->eta()<=-0.1 || cluster->eta()>=1.3) ) || (isTRTB==-1 && (cluster->eta()>=0.1 || cluster->eta()<=-1.3) ) ) { return false; - } - } - - // Identifier barrel or endcap to identify correct calo sampling - int bec = 0; - if ( xAOD::EgammaHelpers::isBarrel( cluster ) ) { - // Barrel - bec = 0; - } else { - // Endcap - bec = 1; - } - - // Setup array for iteration over calo samplings - static CaloCell_ID::CaloSample samples[2][4] = - { { CaloCell_ID::PreSamplerB, - CaloCell_ID::EMB1, - CaloCell_ID::EMB2, - CaloCell_ID::EMB3 }, - { CaloCell_ID::PreSamplerE, - CaloCell_ID::EME1, - CaloCell_ID::EME2, - CaloCell_ID::EME3 } }; - - - static CaloSampling::CaloSample xAODsamples[2][4] = - { { CaloSampling::PreSamplerB, - CaloSampling::EMB1, - CaloSampling::EMB2, - CaloSampling::EMB3 }, - { CaloSampling::PreSamplerE, - CaloSampling::EME1, - CaloSampling::EME2, - CaloSampling::EME3 } }; - - if ( msgLvl (MSG::DEBUG) ) - { - ATH_MSG_DEBUG("Parameters for extrapolation:"); - trkPar->dump( msg() ); - } - - // Extrapolate to test sample - result = m_extrapolateToCalo->extrapolate ((*trkPar), samples[bec][iSampling], offset, - Trk::nonInteracting, direction); - // Save impacts for cleanup - std::vector<const Trk::TrackParameters*> impacts; - - if(result) { - // save deltaEta/deltaPhi - trapped = 0; - deltaEta[iSampling] = cluster->etaSample( xAODsamples[bec][iSampling] ) - result->position().eta(); - deltaPhi[iSampling] = m_phiHelper.diff(cluster->phiSample(xAODsamples[bec][iSampling]),result->position().phi()); - if(flipSign) deltaPhi[iSampling] = -deltaPhi[iSampling]; - impacts.push_back(result); + } } - - // For rescaled perigee, we must delete the trkPar which were - // created on the heap - if (fromPerigeeRescaled == extrapFrom) { - delete trkPar; + //Call getMatchAtCalo + if(getMatchAtCalo (cluster,trkPB,direction, doSample, eta,phi, deltaEta, deltaPhi,extrapFrom).isFailure()){ + ATH_MSG_WARNING("getMatchAtCalo failed"); + return false; } - // Selection in broad eta/phi window - if (trapped != 0 - || (isTRT && abs(isTRTB)==1 && fabs(deltaEta[iSampling]) > m_TRTbarrelDeltaEta) + if ((isTRT && abs(isTRTB)==1 && fabs(deltaEta[iSampling]) > m_TRTbarrelDeltaEta) || (isTRT && abs(isTRTB)==2 && fabs(deltaEta[iSampling]) > m_TRTendcapDeltaEta) || (!isTRT && (fabs(deltaEta[iSampling]) > m_broadDeltaEta)) || fabs(deltaPhi[iSampling]) > m_broadDeltaPhi) { - // cleanup - for (unsigned int i = 0; i < impacts.size(); i++) { - delete impacts[i]; - } return false; } - + // Selection in narrow eta/phi window if((isTRT && ((abs(isTRTB)==1 && deltaPhi[iSampling] < m_narrowDeltaPhiTRTbarrel && @@ -393,89 +231,13 @@ EMExtrapolationTools::matchesAtCalo(const xAOD::CaloCluster* cluster, || (!isTRT && fabs(deltaEta[iSampling]) < m_narrowDeltaEta && deltaPhi[iSampling] < m_narrowDeltaPhi && deltaPhi[iSampling] > -m_narrowDeltaPhiBrem) ) { - - // successful match - eta[iSampling] = result->position().eta(); - phi[iSampling] = result->position().phi(); - - // Now do other requested extrapolations - use previous - // extrapolation endpoint for efficiency - - // Save last impact for the next starting point for extrapolation - const Trk::TrackParameters* lastImpact = result; - // For the first extrapolation, we must check whether or not we - // must invert the extrapolation direction - bool firstExtrap = true; - for (unsigned int i = 0; i < doSample.size(); i++) { - if (i == iSampling || !doSample[i]) continue; // skip unrequested samples - Trk::PropDirection dir = direction; - if (firstExtrap) { - // invert if we go to an earlier sampling - if (i < iSampling) dir = (dir == Trk::alongMomentum) ? - Trk::oppositeMomentum : Trk::alongMomentum; - firstExtrap = false; - } - const Trk::TrackParameters* impact = - m_extrapolateToCalo->extrapolate ((*lastImpact), samples[bec][i], offset, - Trk::nonInteracting, dir); - if (impact) { - eta[i] = impact->position().eta(); - phi[i] = impact->position().phi(); - deltaEta[i] = cluster->etaSample( xAODsamples[bec][i]) - impact->position().eta(); - deltaPhi[i] = m_phiHelper.diff(cluster->phiSample( xAODsamples[bec][i]), impact->position().phi()); - if(flipSign) deltaPhi[i] = -deltaPhi[i]; - lastImpact = impact; - impacts.push_back(impact); - } - } - // cleanup - for (unsigned int i = 0; i < impacts.size(); i++) { - delete impacts[i]; - } return true; } - - // cleanup - for (unsigned int i = 0; i < impacts.size(); i++) { - delete impacts[i]; - } return false; } -// ================================================================= - /** Vertex-to-cluster match **/ -bool EMExtrapolationTools::matchesAtCalo(const xAOD::CaloCluster* cluster, - const xAOD::Vertex *vertex, - float etaAtCalo, - float phiAtCalo) const -{ - if (!cluster || !vertex) return false; - float deltaEta = fabs(etaAtCalo - cluster->etaBE(2)); - float deltaPhi = fabs(m_phiHelper.diff(cluster->phi(), phiAtCalo)); - - int TRTsection = 0; - if (xAOD::EgammaHelpers::numberOfSiTracks(vertex) == 0) - TRTsection = getTRTsection(vertex->trackParticle(0)); - - // First pass on TRT tracks, skip barrel tracks matching endcap clusters and vice-versa - if((TRTsection==2 && (cluster->eta()<=0.6 || cluster->eta()>=2.4) ) || - (TRTsection==-2 && (cluster->eta()>=-0.6 || cluster->eta()<=-2.4) ) || - (TRTsection==1 && (cluster->eta()<=-0.1 || cluster->eta()>=1.3) ) || - (TRTsection==-1 && (cluster->eta()>=0.1 || cluster->eta()<=-1.3) ) - ) { - return false; - } - - // The maximum deltaEta/deltaPhi for Si, TRT barrel, TRT endcap - std::vector<double> dEtaV{m_narrowDeltaEta, m_TRTbarrelDeltaEta, m_TRTendcapDeltaEta}; - std::vector<double> dPhiV{m_narrowDeltaPhi, m_narrowDeltaPhiTRTbarrel, m_narrowDeltaPhiTRTendcap}; - - return (deltaEta < dEtaV[abs(TRTsection)] && deltaPhi < dPhiV[abs(TRTsection)]); -} - - ///////////////////////////////////////////////////////////////////////////////////////////// -//* Extrapolate to calorimeter +//* Extrapolate to calorimeter for trackParticle ///////////////////////////////////////////////////////////////////////////////////////////// StatusCode EMExtrapolationTools::getMatchAtCalo (const xAOD::CaloCluster* cluster, @@ -512,26 +274,13 @@ EMExtrapolationTools::getMatchAtCalo (const xAOD::CaloCluster* cluster, << " deltaPhi: " << deltaPhi.size()); return StatusCode::FAILURE; } - - + const Trk::TrackParameters* trkPar = 0; Trk::CurvilinearParameters temp; if( fromLastMeasurement == extrapFrom ) { - //For TRT it is always the last - int index(-1); - for(unsigned int i(0); i< trkPB->numberOfParameters() ; ++i ){ - if( xAOD::LastMeasurement == trkPB->parameterPosition(i) ){ - index = i; - break; - } - } - if(index!=-1){ - temp = trkPB->curvilinearParameters(index); + temp= getLastMeasurement(trkPB); trkPar = &temp; - } else { - ATH_MSG_WARNING("Track Particle does not contain Last Measurement track parameters"); - } - } + } else if (fromPerigee == extrapFrom) { // From perigee trkPar = &trkPB->perigeeParameters(); @@ -625,128 +374,98 @@ EMExtrapolationTools::getMatchAtCalo (const xAOD::CaloCluster* cluster, } return StatusCode::SUCCESS; - } +// ======================= Photon/vertex Related Interfaces - /** get eta, phi at the four calorimeter layers given the - * Trk::TrackParameters. doSample indicates whether or not to - * extrapolate to each calo sample. Also return whether the - * extrapolation is in barrel or not. - */ -StatusCode -EMExtrapolationTools::getEtaPhiAtCalo (const xAOD::TrackParticle* trkPB, - Trk::PropDirection direction, - const std::vector<bool>& doSample, - std::vector<double>& eta, - std::vector<double>& phi, - std::vector<bool>& isBarrel, - unsigned int extrapFrom) const +/////////////////////////////////////////////////////////////////////////////// +//* Extrapolate double-track conversion with the neutral perigee +/////////////////////////////////////////////////////////////////////////////// +bool +EMExtrapolationTools::matchesAtCalo(const xAOD::CaloCluster* cluster, + const Trk::NeutralParameters* phPerigee, + bool isTRT, + double& deltaEta, + double& deltaPhi) { - // Extrapolate track to calo and return the extrapolated eta/phi - // We allow different ways to extrapolate: - // - // 1) from the last measurement point track parameters - // 2) from the perigee - // 3) from the perigee with the track momentum rescaled by the cluster energy - // + double tmpDeltaEta = -999.; + double tmpDeltaPhi = -999.; + int trapped = 99; - // Checks - // require input vectors to be of the same size: - unsigned int vsz = doSample.size(); - if ((vsz != eta.size()) || (vsz != phi.size())) { - ATH_MSG_ERROR("getMatchAtCalo: input vectors not the same size - doSample: " - << doSample.size() << " eta: " << eta.size() - << " phi: " << phi.size()); - return StatusCode::FAILURE; - } - - // Get track parameters from either last point or perigee - const Trk::TrackParameters* trkPar = 0; - Trk::CurvilinearParameters temp; + const Trk::NeutralParameters* result = 0; + double offset = 0.; + CaloCell_ID::CaloSample sample; - if (fromLastMeasurement == extrapFrom) { - // Use existing track parameters from last hit - int index(-1); - for(unsigned int i(0); i< trkPB->numberOfParameters() ; ++i ){ - if( xAOD::LastMeasurement == trkPB->parameterPosition(i) ){ - index = i; - break; - } + + if ( xAOD::EgammaHelpers::isBarrel( cluster ) ) { + sample = CaloCell_ID::EMB2; + result = m_extrapolateToCalo->extrapolate((*phPerigee), sample, offset, Trk::photon, Trk::alongMomentum); + if (result) { + trapped = 0; + tmpDeltaEta = cluster->etaSample(sample) - result->position().eta(); + tmpDeltaPhi = m_phiHelper.diff(cluster->phiSample(sample),result->position().phi()); } - if(index!=-1){ - temp = trkPB->curvilinearParameters(index); - trkPar = &temp; - } else { - ATH_MSG_WARNING("Track Particle does not contain Last Measurement track parameters"); + }else { + sample = CaloCell_ID::EME2; + result = m_extrapolateToCalo->extrapolate((*phPerigee), sample, offset, Trk::photon, Trk::alongMomentum); + if (result) { + trapped = 0; + tmpDeltaEta = cluster->etaSample(sample) - result->position().eta(); + tmpDeltaPhi = m_phiHelper.diff(cluster->phiSample(sample),result->position().phi()); } } - else if (fromPerigee == extrapFrom) { - trkPar = &trkPB->perigeeParameters(); - } - // return if track parameters not found! - if(!trkPar) { - ATH_MSG_ERROR("getMatchAtCalo: Cannot access track parameters"); - return StatusCode::FAILURE; - } + if(result) delete result; - // Now extrapolate to the samplings. - static CaloCell_ID::CaloSample samples[2][4] = - { { CaloCell_ID::PreSamplerB, - CaloCell_ID::EMB1, - CaloCell_ID::EMB2, - CaloCell_ID::EMB3 }, - { CaloCell_ID::PreSamplerE, - CaloCell_ID::EME1, - CaloCell_ID::EME2, - CaloCell_ID::EME3 } }; + deltaEta = tmpDeltaEta; + deltaPhi = tmpDeltaPhi; - double offset = 0.; - int bec = 0; - const Trk::TrackParameters* lastImpact = trkPar; - const Trk::TrackParameters* dummy(0); - std::vector<const Trk::TrackParameters*> impacts(4, dummy); - for (unsigned int i = 0; i < doSample.size(); i++) { - if (!doSample[i]) continue; // skip unrequested samples - - // Must choose between barrel and endcap - double trketa = lastImpact->position().eta(); - double distbar = m_calodepth->deta(samples[0][i], trketa); - double distec = m_calodepth->deta(samples[1][i], trketa); - if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << " getEtaPhiAtCalo> sample " << i << " : eta= " - << trketa << " dist to Barrel =" - << distbar << " to endcap =" << distec << endreq; - if (distbar < 0 ) bec = 0; - else if (distec < 0 ) bec = 1; - else if ( distbar < distec) bec = 0; - else bec = 1; - const Trk::TrackParameters* impact = - m_extrapolateToCalo->extrapolate ((*lastImpact), samples[bec][i], offset, - Trk::nonInteracting, direction); - if (impact) { - eta[i] = impact->position().eta(); - phi[i] = impact->position().phi(); - isBarrel[i] = (bec == 0); - lastImpact = impact; - impacts.push_back(impact); - } - } - // cleanup extrapolation generated objects - for (unsigned int i = 0; i < impacts.size(); i++) { - delete impacts[i]; + if (trapped != 0 + || (!isTRT && (fabs(deltaEta) > m_broadDeltaEta) ) + || (isTRT && (fabs(deltaEta) > std::max(m_TRTbarrelDeltaEta, m_TRTendcapDeltaEta)) ) + || fabs(deltaPhi) > m_broadDeltaPhi) return false; + + if ((isTRT && fabs(deltaPhi) < std::max(m_narrowDeltaPhiTRTbarrel, m_narrowDeltaPhiTRTendcap)) || + (!isTRT && fabs(deltaEta)< m_narrowDeltaEta && fabs(deltaPhi) < m_narrowDeltaPhi) ) + return true; + + return false; +} + + + +// ================================================================= + /** Vertex-to-cluster match **/ +bool EMExtrapolationTools::matchesAtCalo(const xAOD::CaloCluster* cluster, + const xAOD::Vertex *vertex, + float etaAtCalo, + float phiAtCalo) const +{ + if (!cluster || !vertex) return false; + float deltaEta = fabs(etaAtCalo - cluster->etaBE(2)); + float deltaPhi = fabs(m_phiHelper.diff(cluster->phi(), phiAtCalo)); + + int TRTsection = 0; + if (xAOD::EgammaHelpers::numberOfSiTracks(vertex) == 0) + TRTsection = getTRTsection(vertex->trackParticle(0)); + + // First pass on TRT tracks, skip barrel tracks matching endcap clusters and vice-versa + if((TRTsection==2 && (cluster->eta()<=0.6 || cluster->eta()>=2.4) ) || + (TRTsection==-2 && (cluster->eta()>=-0.6 || cluster->eta()<=-2.4) ) || + (TRTsection==1 && (cluster->eta()<=-0.1 || cluster->eta()>=1.3) ) || + (TRTsection==-1 && (cluster->eta()>=0.1 || cluster->eta()<=-1.3) ) + ) { + return false; } - // For rescaled perigee, we must delete the trkPar which were - // created on the heap - // Can't pass this condition, since this case wasn't considered above. - // Comment out to avoid coverity warning. - //if (fromPerigeeRescaled == extrapFrom) { - // delete trkPar; - //} + + // The maximum deltaEta/deltaPhi for Si, TRT barrel, TRT endcap + std::vector<double> dEtaV{m_narrowDeltaEta, m_TRTbarrelDeltaEta, m_TRTendcapDeltaEta}; + std::vector<double> dPhiV{m_narrowDeltaPhi, m_narrowDeltaPhiTRTbarrel, m_narrowDeltaPhiTRTendcap}; - return StatusCode::SUCCESS; + return (deltaEta < dEtaV[abs(TRTsection)] && deltaPhi < dPhiV[abs(TRTsection)]); } + // ================================================================= bool EMExtrapolationTools::getEtaPhiAtCalo (const xAOD::Vertex* vertex, float *etaAtCalo, @@ -780,100 +499,103 @@ bool EMExtrapolationTools::getEtaPhiAtCalo (const xAOD::Vertex* vertex, return true; } -// ================================================================= + Amg::Vector3D EMExtrapolationTools::getMomentumAtVertex(const xAOD::Vertex& vertex, bool reuse /* = true */) const { return m_extrapolateToCalo->getMomentumAtVertex(vertex, reuse); } -// ================================================================= +// ======================= HELPERS const Trk::TrackParameters* EMExtrapolationTools::getRescaledPerigee(const xAOD::TrackParticle* trkPB, const xAOD::CaloCluster* cluster) const { - const Trk::TrackParameters* oldPerigee = &trkPB->perigeeParameters(); - if (!oldPerigee) return 0; - Amg::Vector3D mom(oldPerigee->momentum().x(),oldPerigee->momentum().y(),oldPerigee->momentum().z()); - const double mag = mom.mag(); - if (mag == 0 || cluster->e() == 0) { - // should this be a warning or debug? For now I assume that it should not happen - ATH_MSG_WARNING("Trying to rescale perigee but cluster e = " << cluster->e() << " and perigee mag = " << mag); - return 0; + const Trk::TrackParameters* oldPerigee = &trkPB->perigeeParameters(); + if (!oldPerigee) return 0; + Amg::Vector3D mom(oldPerigee->momentum().x(),oldPerigee->momentum().y(),oldPerigee->momentum().z()); + const double mag = mom.mag(); + if (mag == 0 || cluster->e() == 0) { + // should this be a warning or debug? For now I assume that it should not happen + ATH_MSG_WARNING("Trying to rescale perigee but cluster e = " << cluster->e() << " and perigee mag = " << mag); + return 0; + } + + // rescale the momentum + mom = mom * Eigen::Scaling(cluster->e()/mom.mag()); + + // Create new perigee + double charge = oldPerigee->charge(); + //Construct the photon neutral parameters as neutral perigee + + Amg::Vector3D pos(oldPerigee->position().x(),oldPerigee->position().y(),oldPerigee->position().z()); + Trk::PerigeeSurface surface (pos); + AmgSymMatrix(5) *covariance = new AmgSymMatrix(5)(*(oldPerigee->covariance())); + const Trk::TrackParameters* result + = surface.createParameters<5,Trk::Charged>(oldPerigee->parameters()[0], + oldPerigee->parameters()[1], + mom.phi(), + mom.theta(), + charge/mom.mag(), + covariance); + return (result); +} + +Trk::CurvilinearParameters EMExtrapolationTools::getLastMeasurement(const xAOD::TrackParticle* trkPB) const{ + // Get last measurement for TRT check + Trk::CurvilinearParameters temp; + + //For TRT it is always the last + int index(-1); + for(unsigned int i(0); i< trkPB->numberOfParameters() ; ++i ){ + if( xAOD::LastMeasurement == trkPB->parameterPosition(i) ){ + index = i; + break; } - - // rescale the momentum - mom = mom * Eigen::Scaling(cluster->e()/mom.mag()); - - // Create new perigee - double charge = oldPerigee->charge(); - //Construct the photon neutral parameters as neutral perigee - - Amg::Vector3D pos(oldPerigee->position().x(),oldPerigee->position().y(),oldPerigee->position().z()); - Trk::PerigeeSurface surface (pos); - AmgSymMatrix(5) *covariance = new AmgSymMatrix(5)(*(oldPerigee->covariance())); - const Trk::TrackParameters* result - = surface.createParameters<5,Trk::Charged>(oldPerigee->parameters()[0], - oldPerigee->parameters()[1], - mom.phi(), - mom.theta(), - charge/mom.mag(), - covariance); - return (result); + } + if(index!=-1){ + temp= trkPB->curvilinearParameters(index); + } else { + ATH_MSG_WARNING("Track Particle does not contain Last Measurement track parameters"); + } + return temp; } + // ================================================================= -int EMExtrapolationTools::getTRTsection(const xAOD::TrackParticle* trkPB) const -{ - if (!trkPB) - { - ATH_MSG_DEBUG("Null pointer to TrackParticle"); - return 0; +int EMExtrapolationTools::getTRTsection(const xAOD::TrackParticle* trkPB) const{ + if (!trkPB){ + ATH_MSG_DEBUG("Null pointer to TrackParticle"); + return 0; } - if (!m_trtId || m_guessTRTsection) - { + if (!m_trtId || m_guessTRTsection){ ATH_MSG_DEBUG("Guessing TRT section based on eta: " << trkPB->eta()); return (trkPB->eta() > 0 ? 1 : -1) * (abs(trkPB->eta()) < 0.6 ? 1 : 2); } const Trk::TrackParameters* trkPar =0; - if (!m_useTrkTrackForTRT) - { - Trk::CurvilinearParameters temp; - int index(-1); - for(unsigned int i(0); i< trkPB->numberOfParameters() ; ++i ){ - if( xAOD::LastMeasurement == trkPB->parameterPosition(i) ){ - index = i; - break; - } - } - if(index!=-1){ - temp = trkPB->curvilinearParameters(index); - trkPar = &temp; - } else { - ATH_MSG_WARNING("Track Particle does not contain Last Measurement track parameters"); - } - + Trk::CurvilinearParameters temp; + + if (!m_useTrkTrackForTRT){ + temp= getLastMeasurement(trkPB); + trkPar = &temp; if(!trkPar){ ATH_MSG_WARNING("Track Particle does not have the requested track parameters"); return 0; } } - else if( trkPB->trackLink().isValid() && trkPB->track() != 0 ) - { + else if( trkPB->trackLink().isValid() && trkPB->track() != 0 ){ ATH_MSG_DEBUG("Will get TrackParameters from Trk::Track"); const DataVector<const Trk::TrackStateOnSurface> *trackStates = trkPB->track()->trackStateOnSurfaces(); - if (!trackStates) - { + if (!trackStates) { ATH_MSG_WARNING("NULL pointer to trackStateOnSurfaces"); return 0; } //Loop over the TrkStateOnSurfaces // search last valid TSOS first - for ( DataVector<const Trk::TrackStateOnSurface>::const_reverse_iterator rItTSoS = trackStates->rbegin(); rItTSoS != trackStates->rend(); ++rItTSoS) - { - if ( (*rItTSoS)->type(Trk::TrackStateOnSurface::Measurement) && (*rItTSoS)->trackParameters()!=0 && (*rItTSoS)->measurementOnTrack()!=0 && !dynamic_cast<const Trk::PseudoMeasurementOnTrack*>((*rItTSoS)->measurementOnTrack())) - { + for ( DataVector<const Trk::TrackStateOnSurface>::const_reverse_iterator rItTSoS = trackStates->rbegin(); rItTSoS != trackStates->rend(); ++rItTSoS){ + if ( (*rItTSoS)->type(Trk::TrackStateOnSurface::Measurement) && (*rItTSoS)->trackParameters()!=0 && (*rItTSoS)->measurementOnTrack()!=0 + && !dynamic_cast<const Trk::PseudoMeasurementOnTrack*>((*rItTSoS)->measurementOnTrack())){ trkPar = (*rItTSoS)->trackParameters(); break; } diff --git a/Reconstruction/egamma/egammaTrackTools/src/EMExtrapolationTools.h b/Reconstruction/egamma/egammaTrackTools/src/EMExtrapolationTools.h index b0fac538bb8f36e9d2500f17f719a8de09e966b1..e215eabdfb0594b6ebc7ed515e2308119f7f07a5 100644 --- a/Reconstruction/egamma/egammaTrackTools/src/EMExtrapolationTools.h +++ b/Reconstruction/egamma/egammaTrackTools/src/EMExtrapolationTools.h @@ -62,14 +62,7 @@ class EMExtrapolationTools : virtual public IEMExtrapolationTools, public AthAlg virtual bool matchesAtCalo(const xAOD::CaloCluster* cluster, const xAOD::TrackParticle* trkPB, bool isTRT, - unsigned int extrapFrom = fromLastMeasurement); - - /** test for cluster/extrapolated track match, from NeutralPerigee */ - virtual bool matchesAtCalo(const xAOD::CaloCluster* cluster, - const Trk::NeutralParameters* perigee, - bool isTRT, - double& deltaEta, - double& deltaPhi); + unsigned int extrapFrom = fromPerigee); /** test for cluster/extrapolated track match, from xAOD::TrackParticle, * returns true for good match, and @@ -84,14 +77,8 @@ class EMExtrapolationTools : virtual public IEMExtrapolationTools, public AthAlg std::vector<double>& phi, std::vector<double>& deltaEta, std::vector<double>& deltaPhi, - unsigned int extrapFrom = fromLastMeasurement) const; + unsigned int extrapFrom = fromPerigee) const; - /** test for vertex-to-cluster match given also the positions - * at the calorimeter from the vertex extrapolation **/ - bool matchesAtCalo(const xAOD::CaloCluster* cluster, - const xAOD::Vertex *vertex, - float etaAtCalo, - float phiAtCalo) const; /** get eta, phi, deltaEta, and deltaPhi at the four calorimeter * layers given the Trk::ParametersBase. doSample indicates @@ -105,22 +92,22 @@ class EMExtrapolationTools : virtual public IEMExtrapolationTools, public AthAlg std::vector<double>& phi, std::vector<double>& deltaEta, std::vector<double>& deltaPhi, - unsigned int extrapFrom = fromLastMeasurement) const; + unsigned int extrapFrom = fromPerigee) const; + /** test for vertex-to-cluster match given also the positions + * at the calorimeter from the vertex extrapolation **/ + bool matchesAtCalo(const xAOD::CaloCluster* cluster, + const xAOD::Vertex *vertex, + float etaAtCalo, + float phiAtCalo) const; - /** get eta, phi at the four calorimeter layers given the - * Trk::TrackParameters. doSample indicates whether or not to - * extrapolate to each calo sample. Also return whether the - * extrapolation is in barrel or not. - */ - virtual StatusCode getEtaPhiAtCalo (const xAOD::TrackParticle* trkPB, - Trk::PropDirection direction, - const std::vector<bool>& doSample, - std::vector<double>& eta, - std::vector<double>& phi, - std::vector<bool>& isBarrel, - unsigned int extrapFrom = fromLastMeasurement) const; + /** test for cluster/extrapolated track match, from NeutralPerigee */ + virtual bool matchesAtCalo(const xAOD::CaloCluster* cluster, + const Trk::NeutralParameters* perigee, + bool isTRT, + double& deltaEta, + double& deltaPhi); /** get eta, phi at EM2 given a vertex which is converted to NeutralParameters. Return false if the extrapolation fails **/ @@ -130,12 +117,14 @@ class EMExtrapolationTools : virtual public IEMExtrapolationTools, public AthAlg /** get sum of the momenta at the vertex (designed for conversions). Retrieve from auxdata if available and <reuse> is true **/ Amg::Vector3D getMomentumAtVertex(const xAOD::Vertex&, bool reuse = true) const; + + private: /** @brief Return +/- 1 (2) if track is in positive/negative TRT barrel (endcap) **/ int getTRTsection(const xAOD::TrackParticle* trkPB) const; - - private: - + + Trk::CurvilinearParameters getLastMeasurement(const xAOD::TrackParticle* trkPB) const; + const Trk::TrackParameters* getRescaledPerigee(const xAOD::TrackParticle* trkPB, const xAOD::CaloCluster* cluster) const; /** @brief TrackToCalo extrapolation tool. @@ -150,9 +139,6 @@ class EMExtrapolationTools : virtual public IEMExtrapolationTools, public AthAlg /** @brief */ CaloPhiRange m_phiHelper; - /** @brief Utility to calculate the converted photon track parameters and covariance matrix*/ - ToolHandle<Trk::INeutralParticleParameterCalculator> m_convUtils; - // Track-to-cluster match cuts double m_broadDeltaEta;