egammaSelectedTrackCopy.cxx 12.3 KB
Newer Older
1
/*
2
   Copyright (C) 2002-2020  CERN for the benefit of the ATLAS collaboration
3
4
5
6
7
8
9
 */

/*
NAME:     egammaSelectedTrackCopy
PACKAGE:  offline/Reconstruction/egamma/egammaAlgs/egammaSelectedTrackCopy
AUTHORS:  Anastopoulos
CREATED:  25/06/2018
10

11
12
13
14
15
16
17
PURPOSE: Select track to be refitted later on with GSF
UPDATE : 25/06/2018
*/

#include "egammaSelectedTrackCopy.h"
//
#include "AthenaKernel/errorcheck.h"
18
19
#include "FourMomUtils/P4Helpers.h"
#include "egammaUtils/CandidateMatchHelpers.h"
20
21
22
#include "xAODCaloEvent/CaloCluster.h"
#include "xAODCaloEvent/CaloClusterContainer.h"
#include "xAODCore/ShallowCopy.h"
23
24
25
#include "xAODEgamma/EgammaxAODHelpers.h"
#include "xAODTracking/TrackParticle.h"
#include "xAODTracking/TrackParticleContainer.h"
26
#include "CaloDetDescr/CaloDetDescrManager.h"
27
#include "GaudiKernel/EventContext.h"
28
29
#include "StoreGate/ReadHandle.h"
#include "StoreGate/WriteHandle.h"
30
// std includes
31
32
33
#include <algorithm>
#include <cmath>
#include <memory>
34
#include <stdexcept>
35
#include <vector>
36

37
38
39
40
41
42
43
44
45
46
47
48
49
50
egammaSelectedTrackCopy::egammaSelectedTrackCopy(const std::string& name, ISvcLocator* pSvcLocator)
  : AthAlgorithm(name, pSvcLocator)
  , m_AllClusters{}
  , m_SelectedClusters{}
  , m_AllTracks{}
  , m_SelectedTracks{}
  , m_AllSiTracks{}
  , m_SelectedSiTracks{}
  , m_AllTRTTracks{}
  , m_SelectedTRTTracks{}
{}

StatusCode
egammaSelectedTrackCopy::initialize()
51
52
53
54
55
56
{

  ATH_CHECK(m_clusterContainerKey.initialize());
  ATH_CHECK(m_trackParticleContainerKey.initialize());
  ATH_CHECK(m_OutputTrkPartContainerKey.initialize());

57
58
59
  ATH_CHECK(m_pixelDetEleCollKey.initialize());
  ATH_CHECK(m_SCTDetEleCollKey.initialize());

60
  /* the extrapolation tool*/
61
62
  ATH_CHECK(m_extrapolationTool.retrieve());
  ATH_CHECK(m_extrapolationToolCommonCache.retrieve());
63

64
  ATH_CHECK(m_egammaCaloClusterSelector.retrieve());
65

66
  return StatusCode::SUCCESS;
67
}
68

69
70
71
StatusCode
egammaSelectedTrackCopy::egammaSelectedTrackCopy::finalize()
{
72

73
74
75
76
77
78
79
80
81
82
  ATH_MSG_INFO("--- egamma Selected Track Copy Statistics ---");
  ATH_MSG_INFO("--- All Clusters: " << m_AllClusters);
  ATH_MSG_INFO("--- Selected Clusters: " << m_SelectedClusters);
  ATH_MSG_INFO("--- All Tracks: " << m_AllTracks);
  ATH_MSG_INFO("--- Selected Tracks: " << m_SelectedTracks);
  ATH_MSG_INFO("--- All Si Tracks: " << m_AllSiTracks);
  ATH_MSG_INFO("--- Selected Si Tracks: " << m_SelectedSiTracks);
  ATH_MSG_INFO("--- All TRT Tracks: " << m_AllTRTTracks);
  ATH_MSG_INFO("--- Selected TRT Tracks: " << m_SelectedTRTTracks);
  ATH_MSG_INFO("----------------------------------------- ---");
83
84
85
  return StatusCode::SUCCESS;
}

86
87
StatusCode
egammaSelectedTrackCopy::execute()
88
89
{
  SG::ReadHandle<xAOD::CaloClusterContainer> clusterTES(m_clusterContainerKey);
90
91
  if (!clusterTES.isValid()) {
    ATH_MSG_FATAL("Failed to retrieve cluster container: " << m_clusterContainerKey.key());
92
93
94
    return StatusCode::FAILURE;
  }
  SG::ReadHandle<xAOD::TrackParticleContainer> trackTES(m_trackParticleContainerKey);
95
96
  if (!trackTES.isValid()) {
    ATH_MSG_FATAL("Failed to retrieve TrackParticle container: " << m_trackParticleContainerKey.key());
97
98
99
    return StatusCode::FAILURE;
  }

100
  SG::WriteHandle<ConstDataVector<xAOD::TrackParticleContainer>> outputTrkPartContainer(m_OutputTrkPartContainerKey);
101
  /*
102
103
   * Here it just needs to be a view copy ,
   * i.e the collection of selected trackParticles
104
105
   * we create does not really own its elements
   */
106
  auto viewCopy = std::make_unique<ConstDataVector<xAOD::TrackParticleContainer>>(SG::VIEW_ELEMENTS);
107

108
  // Local counters
109
110
  auto allClusters = m_AllClusters.buffer();
  auto selectedClusters = m_SelectedClusters.buffer();
111
112
113
114
115
116
  auto allTracks = m_AllTracks.buffer();
  auto selectedTracks = m_SelectedTracks.buffer();
  auto allSiTracks = m_AllSiTracks.buffer();
  auto selectedSiTracks = m_SelectedSiTracks.buffer();
  auto allTRTTracks = m_AllTRTTracks.buffer();
  auto selectedTRTTracks = m_SelectedTRTTracks.buffer();
117

118
119
  const CaloDetDescrManager* calodetdescrmgr = nullptr;
  ATH_CHECK( detStore()->retrieve(calodetdescrmgr,"CaloMgr") );
120
  // // lets first check which clusters to seed on;
121
122
123
  std::vector<const xAOD::CaloCluster*> passingClusters;
  for (const xAOD::CaloCluster* cluster : *clusterTES) {
    ++allClusters;
124
    if (m_egammaCaloClusterSelector->passSelection(cluster,*calodetdescrmgr)) {
125
      passingClusters.push_back(cluster);
126
      ++selectedClusters;
127
128
    }
  }
129

130
  // Extrapolation Cache
131
  IEMExtrapolationTools::Cache cache{};
132
  for (const xAOD::TrackParticle* track : *trackTES) {
133

134
    ++allTracks;
135
    bool isTRT = false;
136
    int nhits(0);
137
138
    uint8_t dummy(0);
    if (track->summaryValue(dummy, xAOD::numberOfPixelHits)) {
139
140
      nhits += dummy;
    }
141
142
    if (track->summaryValue(dummy, xAOD::numberOfSCTHits)) {
      nhits += dummy;
143
    }
144
    if (nhits < 4) {
145
146
      isTRT = true;
      ++allTRTTracks;
147
    } else {
148
149
150
      isTRT = false;
      ++allSiTracks;
    }
151

152
    for (const xAOD::CaloCluster* cluster : passingClusters) {
153
      /*
154
155
156
         check if it the track is selected due to this cluster.
         If not continue to next cluster
         */
157
      if (!Select(Gaudi::Hive::currentContext(), cluster, track, cache, isTRT)) {
158
159
        continue;
      }
160
      viewCopy->push_back(track);
161
      ++selectedTracks;
162
      if (isTRT) {
163
        ++selectedTRTTracks;
164
      } else {
165
166
167
        ++selectedSiTracks;
      }
      /*
168
       * The particular track got selected
169
170
171
172
       * due to a cluster (any one of them will do)
       * break here and move to the next track
       */
      break;
173
174
    } // Loop on clusters
  }   // Loop on tracks
175

176
  ATH_CHECK(outputTrkPartContainer.record(std::move(viewCopy)));
177
178
179
180

  return StatusCode::SUCCESS;
}

181
182
183
184
185
186
bool
egammaSelectedTrackCopy::Select(const EventContext& ctx,
                                const xAOD::CaloCluster* cluster,
                                const xAOD::TrackParticle* track,
                                IEMExtrapolationTools::Cache& cache,
                                bool trkTRT) const
187
188
{

189
190
191
192
  if (cluster == nullptr || track == nullptr) {
    return false;
  }

193
194
  const Trk::Perigee& candidatePerigee = track->perigeeParameters();

195
196
197
198
199
200
201
202
203
204
205
206
  // Get Perigee Parameters
  const double trkPhi = candidatePerigee.parameters()[Trk::phi];
  const double trkEta = candidatePerigee.eta();
  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);
  const bool isEndCap = !xAOD::EgammaHelpers::isBarrel(cluster);
  double Et = cluster->e() / cosh(trkEta);
  if (trkTRT) {
207
208
209
    Et = cluster->et();
  }
  // a few sanity checks
210
211
  if (fabs(clusterEta) > 10.0 || fabs(trkEta) > 10.0 || Et < 10.0) {
    ATH_MSG_DEBUG("FAILS sanity checks :  Track Eta : " << trkEta << ", Cluster Eta " << clusterEta);
212
213
214
    return false;
  }

215
216
217
218
219
220
221
222
  // 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);
223
  //
224
225
  const double clusterPhiCorrected = XYZClusterWrtTrackPerigee.phi();
  // deltaPhi between the track and the cluster
226
  const double deltaPhiStd = P4Helpers::deltaPhi(clusterPhiCorrected, trkPhi);
227
228
  // deltaPhi between the track and the cluster accounting for rotation assuming cluster Et is a better estimator
  const double trkPhiRescaled = P4Helpers::deltaPhi(trkPhi, phiRotRescaled);
229
  const double deltaPhiRescaled = P4Helpers::deltaPhi(clusterPhiCorrected, trkPhiRescaled);
230
  // deltaPhi between the track and the cluster accounting for rotation
231
  const double trkPhiCorrTrack = P4Helpers::deltaPhi(trkPhi, phiRotTrack);
232
  const double deltaPhiTrack = P4Helpers::deltaPhi(clusterPhiCorrected, trkPhiCorrTrack);
233
234
  /*
   * First we will see if it fails the quick match
235
236
237
   * 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
238
   */
239
240
241
  // Broad phi check
  if ((fabs(deltaPhiRescaled) > m_broadDeltaPhi) && (fabs(deltaPhiTrack) > m_broadDeltaPhi) &&
      (fabs(deltaPhiStd) > m_broadDeltaPhi)) {
242
    ATH_MSG_DEBUG("FAILS broad window phi match (track phi, phirotCluster , phiRotTrack , "
243
244
                  << "cluster phi corrected, cluster phi): ( " << trkPhi << ", " << phiRotRescaled << ", "
                  << phiRotTrack << ", " << clusterPhiCorrected << ", " << cluster->phi() << ")");
245
246
    return false;
  }
247
248
  // if TRT we can stop here , we can not check much in eta really.
  if (trkTRT) {
249
250
    return true;
  }
251
  const double clusterEtaCorrected = XYZClusterWrtTrackPerigee.eta();
252

253
254
255
256
257
  // Broad eta check
  if (fabs(cluster->eta() - trkEta) >  m_broadDeltaEta &&
      fabs(clusterEtaCorrected - trkEta) > m_broadDeltaEta) {
    ATH_MSG_DEBUG("FAILS broad window eta match (track eta, cluster eta, cluster eta corrected): ( "
                  << trkEta << ", " << cluster->etaBE(2) << ", " << clusterEtaCorrected << " )");
258
259
260
    return false;
  }

261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
  // Extrapolate from last measurement, since this is before brem fit last measurement is OK.
  std::array<double, 4> eta = { -999.0, -999.0, -999.0, -999.0 };
  std::array<double, 4> phi = { -999.0, -999.0, -999.0, -999.0 };
  std::array<double, 4> deltaEta = { -999.0, -999.0, -999.0, -999.0 };
  std::array<double, 4> deltaPhi = { -999.0, -999.0, -999.0, -999.0 };
  // We can try to use the common cache with Pflow/Tau if available or an in-algorithm cache
  if (m_extrapolationToolCommonCache
        ->getMatchAtCalo(ctx,
                         *cluster,
                         *track,
                         Trk::alongMomentum,
                         eta,
                         phi,
                         deltaEta,
                         deltaPhi,
                         IEMExtrapolationTools::fromLastMeasurement)
        .isFailure()) {
    if (m_extrapolationTool
          ->getMatchAtCalo(ctx,
                           *cluster,
                           *track,
                           Trk::alongMomentum,
                           eta,
                           phi,
                           deltaEta,
                           deltaPhi,
                           IEMExtrapolationTools::fromLastMeasurement,
                           &cache)
          .isFailure()) {
290
      return false;
291
    }
292
  }
293
  // Selection in narrow eta/phi window from last measurement
294
295
  if (fabs(deltaEta[2]) < m_narrowDeltaEta && deltaPhi[2] < m_narrowDeltaPhi && deltaPhi[2] > -m_narrowDeltaPhiBrem) {
    ATH_MSG_DEBUG("Match from Last measurement is successful :  " << deltaPhi[2]);
296
297
    return true;
  }
298
  /*
299
   * Cases where 
300
   * - it passes  deltaEta[2] from last measurement (rescaling should not affect the eta side)
301
302
   * - and we have a cluster with higher Et.
   * Rescale up the track to account for radiative loses and retry
303
   */
304
  if (fabs(deltaEta[2]) < m_narrowDeltaEta && cluster->et() > track->pt()) {
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
    // Extrapolate from Perigee Rescaled
    std::array<double, 4> etaRes = { -999.0, -999.0, -999.0, -999.0 };
    std::array<double, 4> phiRes = { -999.0, -999.0, -999.0, -999.0 };
    std::array<double, 4> deltaEtaRes = { -999.0, -999.0, -999.0, -999.0 };
    std::array<double, 4> deltaPhiRes = { -999.0, -999.0, -999.0, -999.0 };

    if (m_extrapolationTool
          ->getMatchAtCalo(ctx,
                           *cluster,
                           *track,
                           Trk::alongMomentum,
                           etaRes,
                           phiRes,
                           deltaEtaRes,
                           deltaPhiRes,
                           IEMExtrapolationTools::fromPerigeeRescaled)
          .isFailure()) {
322
323
      return false;
    }
324
325
326
327
    // Redo the check with rescale
    if (fabs(deltaEtaRes[2]) < m_narrowDeltaEta && deltaPhiRes[2] < m_narrowRescale &&
        deltaPhiRes[2] > -m_narrowRescaleBrem) {
      ATH_MSG_DEBUG("Rescale Match success " << deltaPhiRes[2]);
328
329
330
      return true;
    }
  }
331
  //default is fail
332
333
  return false;
}