TauTrackFinder.cxx 24.4 KB
Newer Older
1
/*
2
  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
3
4
*/

5
#ifndef XAOD_ANALYSIS
6
7
8
9
#include "TrkToolInterfaces/ITrackSelectorTool.h"
#include "TrkParametersIdentificationHelpers/TrackParametersIdHelper.h"

#include "xAODTau/TauJet.h"
10
#include "xAODTau/TauTrackContainer.h"
11
12
13
14
15

#include "TauTrackFinder.h"
#include "tauRecTools/TrackSort.h"


16
17
TauTrackFinder::TauTrackFinder(const std::string& name) :
  TauRecToolBase(name) {
18
19
  m_EMSamplings = {CaloSampling::EME1, CaloSampling::EMB1};
  m_HadSamplings = {CaloSampling::TileBar1, CaloSampling::HEC1, CaloSampling::TileExt1};
20
21
22
23
24
25
26
27
}

TauTrackFinder::~TauTrackFinder() {
}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
StatusCode TauTrackFinder::initialize() {

28
29
30
31
32
33
34
  // retrieve tools
  ATH_CHECK( m_trackSelectorTool_tau.retrieve() );
  ATH_CHECK( m_trackToVertexTool.retrieve() );
  ATH_CHECK( m_caloExtensionTool.retrieve() );
  ATH_CHECK( m_trackToVertexIPEstimator.retrieve() );

  // initialize ReadHandleKey
35
  ATH_CHECK( m_trackPartInputContainer.initialize() );
36
37
38
39
40
  // use CaloExtensionTool when key is empty 
  ATH_CHECK( m_ParticleCacheKey.initialize(SG::AllowEmpty) );
  // allow empty for LRT
  ATH_CHECK( m_largeD0TracksInputContainer.initialize(SG::AllowEmpty) );

41
42
43
44
45
46
  if(m_useGhostTracks) {
    if(inTrigger()) {
      ATH_MSG_ERROR ("Ghost matching is not a valid tau-track association scheme for trigger, use cone association. Aborting.");
      return StatusCode::FAILURE;
    }
    ATH_MSG_INFO ("Using ghost matching for tau-track association" << m_ghostTrackDR );
47
48
    // allow empty for trigger 
    ATH_CHECK( m_jetContainer.initialize(SG::AllowEmpty) );
49
50
  }

51
52
53
  if (inTrigger()) {
    ATH_CHECK(m_beamSpotKey.initialize());
  }
54

55
  return StatusCode::SUCCESS;
56
57
58
}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
59
StatusCode TauTrackFinder::executeTrackFinder(xAOD::TauJet& pTau, xAOD::TauTrackContainer& tauTrackCon) const {
60
61
62
63
  
  std::vector<const xAOD::TrackParticle*> tauTracks;
  std::vector<const xAOD::TrackParticle*> wideTracks;
  std::vector<const xAOD::TrackParticle*> otherTracks;
64

65
  //Retrieve standard track container
66
67
  const xAOD::TrackParticleContainer* trackParticleCont = nullptr; 
  
68
69
70
71
  SG::ReadHandle<xAOD::TrackParticleContainer> trackPartInHandle( m_trackPartInputContainer );
  if (!trackPartInHandle.isValid()) {
    ATH_MSG_ERROR ("Could not retrieve HiveDataObj with key " << trackPartInHandle.key());
    return StatusCode::FAILURE;
72
  }
73
  trackParticleCont = trackPartInHandle.cptr();
74

75
  //Retrieve LRT container
76
77
78
79
80
81
82
83
  const xAOD::TrackParticleContainer* largeD0TracksParticleCont = nullptr; 
  
  if (! m_largeD0TracksInputContainer.empty()) { 
    SG::ReadHandle<xAOD::TrackParticleContainer> trackPartInHandle( m_largeD0TracksInputContainer );
    if (!trackPartInHandle.isValid()) {
      ATH_MSG_VERBOSE ("Could not retrieve HiveDataObj with key " << trackPartInHandle.key());
      ATH_MSG_VERBOSE ("LRT container " << trackPartInHandle.key()<<" is not being used for tau tracks");
    }
84
85
86
    else { 
      largeD0TracksParticleCont = trackPartInHandle.cptr();
    }  
87
88
  }

89
90
91
92
93
94
95
96
97
98
99
  // retrieve the seed jet container when using ghost-matching
  const xAOD::JetContainer* jetContainer = nullptr; 
  if (! m_jetContainer.empty()) {
    SG::ReadHandle<xAOD::JetContainer> jetContHandle( m_jetContainer );
    if (!jetContHandle.isValid()) {
      ATH_MSG_ERROR ("Could not retrieve HiveDataObj with key " << jetContHandle.key());
      return StatusCode::FAILURE;
    }
    jetContainer = jetContHandle.cptr();
  }

100
  // get the primary vertex
101
102
  const xAOD::Vertex* pVertex = nullptr;
  if (pTau.vertexLink().isValid()) pVertex = pTau.vertex();
103
104
105

  // retrieve tracks wrt a vertex                                                                                                                              
  // as a vertex is used: tau origin / PV / beamspot / 0,0,0 (in this order, depending on availability)                                                        
106

107
  getTauTracksFromPV(pTau, *trackParticleCont, pVertex, m_useGhostTracks, jetContainer, tauTracks, wideTracks, otherTracks);
108

109
110
111
  bool foundLRTCont = bool (largeD0TracksParticleCont != nullptr);
  // additional LRT with vertex association added to tracks
  if (foundLRTCont){
112
    // for now, use cone association for LRTs, not ghost association
113
    getTauTracksFromPV(pTau, *largeD0TracksParticleCont, pVertex, false, nullptr, tauTracks, wideTracks, otherTracks);
114
  }
115

116
117
118
119
120
121
122
123
124
125
126
  // remove core and wide tracks outside a maximal delta z0 wrt lead core track                                                                                
  if (m_applyZ0cut) {
    this->removeOffsideTracksWrtLeadTrk(tauTracks, wideTracks, otherTracks, pVertex, m_z0maxDelta);
  }

  if(m_removeDuplicateCoreTracks){
    bool alreadyUsed = false;
    for (std::vector<const xAOD::TrackParticle*>::iterator track_it = tauTracks.begin(); track_it != tauTracks.end() ;)
      {
	alreadyUsed = false;
	//loop over all up-to-now core tracks	
127
	for( const xAOD::TauTrack* tau_trk : tauTrackCon ) {
128
129
	  //originally it was coreTrack&passTrkSelector
	  if(! tau_trk->flagWithMask( (1<<xAOD::TauJetParameters::TauTrackFlag::coreTrack) | (1<<xAOD::TauJetParameters::TauTrackFlag::passTrkSelector))) continue; 
130
131
132
133
	  if( (*track_it) == tau_trk->track()) alreadyUsed = true;
	}
	//if this track has already been used by another tau, don't associate it to this new one                                                               
	if(alreadyUsed) ATH_MSG_INFO( "Found Already Used track new, now removing: " << *track_it );
134
	if (alreadyUsed) track_it = tauTracks.erase(track_it);
135
136
137
	else ++track_it;
      }
  }
138

139
  // associate track to tau candidate and calculate charge                                                                                                    
Bertrand Martin's avatar
Bertrand Martin committed
140
  float charge = 0.;  
141
142
  for (unsigned int i = 0; i < tauTracks.size(); ++i) {
    const xAOD::TrackParticle* trackParticle = tauTracks.at(i);
143

144
    ATH_MSG_VERBOSE(name() << " adding core track nr: " << i
145
		    << " eta " << trackParticle->eta()
146
		    << " phi " << trackParticle->phi());
147

148
    charge += trackParticle->charge();
149

150
    xAOD::TauTrack* track = new xAOD::TauTrack();
151
    tauTrackCon.push_back(track);
152

153
    ElementLink<xAOD::TrackParticleContainer> linkToTrackParticle;
154
155
    if (foundLRTCont && isLargeD0Track(trackParticle)){//Check LRT track and link to container
      linkToTrackParticle.toContainedElement(*largeD0TracksParticleCont, trackParticle);
156
      track->setFlag(xAOD::TauJetParameters::TauTrackFlag::LargeRadiusTrack, true);
157
158
    }
    else {
159
      linkToTrackParticle.toContainedElement(*trackParticleCont, trackParticle);
160
      track->setFlag(xAOD::TauJetParameters::TauTrackFlag::LargeRadiusTrack, false);
161
    }
162
    track->addTrackLink(linkToTrackParticle);
163

164
165
166
167
168
169
    track->setP4(trackParticle->pt(), trackParticle->eta(), trackParticle->phi(), trackParticle->m());
    track->setFlag(xAOD::TauJetParameters::TauTrackFlag::coreTrack, true);
    track->setFlag(xAOD::TauJetParameters::TauTrackFlag::passTrkSelector, true);
    // in case TrackClassifier is not run, still get sensible results                                                                                        
    track->setFlag(xAOD::TauJetParameters::TauTrackFlag::classifiedCharged, true);
    track->setFlag(xAOD::TauJetParameters::TauTrackFlag::unclassified, true);
170
    
171
    ElementLink<xAOD::TauTrackContainer> linkToTauTrack;
172
    linkToTauTrack.toContainedElement(tauTrackCon, track);
173
174
    pTau.addTauTrackLink(linkToTauTrack);

175
    ATH_MSG_VERBOSE(name() << " added core track nr: " << i
176
		    << " eta " << pTau.track(i)->eta()
177
		    << " phi " << pTau.track(i)->phi());
178
179
180
181
182
183
184
  }
  // set the charge, which is defined by the core tau tracks only                                                                                              
  pTau.setCharge(charge);

  for (unsigned int i = 0; i < wideTracks.size(); ++i) {
    const xAOD::TrackParticle* trackParticle = wideTracks.at(i);

185
    ATH_MSG_VERBOSE(name() << " adding wide track nr: " << i
186
		    << " eta " << trackParticle->eta()
187
		    << " phi " << trackParticle->phi());
188
189

    xAOD::TauTrack* track = new xAOD::TauTrack();
190
    tauTrackCon.push_back(track);
191
192

    ElementLink<xAOD::TrackParticleContainer> linkToTrackParticle;
193
194
    if (foundLRTCont && isLargeD0Track(trackParticle)){//Check LRT track and link to container
      linkToTrackParticle.toContainedElement(*largeD0TracksParticleCont, trackParticle);
195
      track->setFlag(xAOD::TauJetParameters::TauTrackFlag::LargeRadiusTrack, true);
196
197
    }
    else {
198
      linkToTrackParticle.toContainedElement(*trackParticleCont, trackParticle);
199
      track->setFlag(xAOD::TauJetParameters::TauTrackFlag::LargeRadiusTrack, false);
200
    }
201
202
203
204
205
206
207
208
209
210
211
    track->addTrackLink(linkToTrackParticle);

    track->setP4(trackParticle->pt(), trackParticle->eta(), trackParticle->phi(), trackParticle->m());
    track->setFlag(xAOD::TauJetParameters::TauTrackFlag::wideTrack, true);
    track->setFlag(xAOD::TauJetParameters::TauTrackFlag::passTrkSelector, true);
    // in case TrackClassifier is not run, still get sensible results                                                                                        
    track->setFlag(xAOD::TauJetParameters::TauTrackFlag::classifiedIsolation, true); // for sake of trigger, reset in TauTrackClassifier
    track->setFlag(xAOD::TauJetParameters::TauTrackFlag::modifiedIsolationTrack, true); // for sake of trigger, reset in TauTrackClassifier
    track->setFlag(xAOD::TauJetParameters::TauTrackFlag::unclassified, true);

    ElementLink<xAOD::TauTrackContainer> linkToTauTrack;
212
    linkToTauTrack.toContainedElement(tauTrackCon, track);
213
214
215
216
217
218
219
220
221
222
223
224
    pTau.addTauTrackLink(linkToTauTrack);
  }

  //These are set again in TauTrackClassifier                                                                                                                  
  pTau.setDetail(xAOD::TauJetParameters::nChargedTracks, (int) pTau.nTracks());
  pTau.setDetail(xAOD::TauJetParameters::nIsolatedTracks, (int) pTau.nTracks(xAOD::TauJetParameters::classifiedIsolation));

  for (unsigned int i = 0; i < otherTracks.size(); ++i) {
    const xAOD::TrackParticle* trackParticle = otherTracks.at(i);

    ATH_MSG_VERBOSE(name()     << " adding other track nr: " << i
		    << " eta " << trackParticle->eta()
225
		    << " phi " << trackParticle->phi());
226
227

    xAOD::TauTrack* track = new xAOD::TauTrack();
228
    tauTrackCon.push_back(track);
229
230

    ElementLink<xAOD::TrackParticleContainer> linkToTrackParticle;
231
232
    if (foundLRTCont && isLargeD0Track(trackParticle)){//Check LRT track and link to container
      linkToTrackParticle.toContainedElement(*largeD0TracksParticleCont, trackParticle);
233
      track->setFlag(xAOD::TauJetParameters::TauTrackFlag::LargeRadiusTrack, true); 
234
235
    }
    else {
236
      linkToTrackParticle.toContainedElement(*trackParticleCont, trackParticle);
237
      track->setFlag(xAOD::TauJetParameters::TauTrackFlag::LargeRadiusTrack, false);
238
    }
239
240
241
242
243
244
245
246
247
    track->addTrackLink(linkToTrackParticle);

    track->setP4(trackParticle->pt(), trackParticle->eta(), trackParticle->phi(), trackParticle->m());
    float dR = track->p4().DeltaR(pTau.p4());
    if(dR<=0.2) track->setFlag(xAOD::TauJetParameters::TauTrackFlag::coreTrack, true);
    else track->setFlag(xAOD::TauJetParameters::TauTrackFlag::wideTrack, true);
    track->setFlag(xAOD::TauJetParameters::TauTrackFlag::unclassified, true);

    ElementLink<xAOD::TauTrackContainer> linkToTauTrack;
248
    linkToTauTrack.toContainedElement(tauTrackCon, track);
249
250
251
    pTau.addTauTrackLink(linkToTauTrack);
  }

252
253
254
255
  pTau.setDetail(xAOD::TauJetParameters::nLargeRadiusTracks, (int) pTau.nTracks(xAOD::TauJetParameters::LargeRadiusTrack));
  // keep track of total number of associated tracks, in case of tau track thinning
  pTau.setDetail(xAOD::TauJetParameters::nAllTracks, (int) pTau.nAllTracks());

256
257
258
  ATH_MSG_DEBUG("numTrack: " << "/" << pTau.nTracks());
  ATH_MSG_DEBUG("charge: " << "/" << pTau.charge());

259
260
261
  // impact parameter variables w.r.t. tau vertex 
  const xAOD::Vertex* vxcand = nullptr;

262
263
264
265
266
  xAOD::Vertex vxbkp;
  vxbkp.makePrivateStore();

  if (pTau.vertexLink().isValid() && pTau.vertex()->vertexType() != xAOD::VxType::NoVtx) {
    vxcand = pTau.vertex();
267
268
  }
  else if (inTrigger()) { // online: use vertex with x-y coordinates from the beamspot and the z from the leading track
269
    vxbkp.setX(0); vxbkp.setY(0); vxbkp.setZ(0);
270
271
272

    SG::ReadCondHandle<InDet::BeamSpotData> beamSpotHandle { m_beamSpotKey };
    if(beamSpotHandle.isValid()) {
273
      vxbkp.setPosition(beamSpotHandle->beamPos());
274
      const auto& cov = beamSpotHandle->beamVtx().covariancePosition();
275
276
277
278
279
      vxbkp.setCovariancePosition(cov);

      if(!tauTracks.empty()) {
         vxbkp.setZ(tauTracks.at(0)->z0());
      }
280
281
282
283
    }
    else {
      ATH_MSG_DEBUG("No Beamspot object in tau candidate");
    }
284
    vxcand = & vxbkp;
285
  }
286

287

288
289
290
291
292
  // this could be replaced with TauTrack::setDetail
  static const SG::AuxElement::Accessor<float> dec_d0TJVA("d0TJVA");
  static const SG::AuxElement::Accessor<float> dec_z0sinthetaTJVA("z0sinthetaTJVA");
  static const SG::AuxElement::Accessor<float> dec_d0SigTJVA("d0SigTJVA");
  static const SG::AuxElement::Accessor<float> dec_z0sinthetaSigTJVA("z0sinthetaSigTJVA");
293

294
  for(auto *track : pTau.allTracks()) {      
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
    dec_d0TJVA(*track) = track->track()->d0();
    dec_z0sinthetaTJVA(*track) = track->z0sinThetaTJVA(pTau);
    dec_d0SigTJVA(*track) = -999.;
    dec_z0sinthetaSigTJVA(*track) = -999.;

    // in the trigger, z0sintheta and corresponding significance are meaningless if we use the beamspot
    if(vxcand) {
      std::unique_ptr<const Trk::ImpactParametersAndSigma> myIPandSigma 
	= std::unique_ptr<const Trk::ImpactParametersAndSigma>(m_trackToVertexIPEstimator->estimate(track->track(), vxcand));
      
      if(myIPandSigma) {
	dec_d0TJVA(*track) = myIPandSigma->IPd0;
	dec_z0sinthetaTJVA(*track) = myIPandSigma->IPz0SinTheta;
	dec_d0SigTJVA(*track) = (myIPandSigma->sigmad0 != 0.) ? (float)( myIPandSigma->IPd0 / myIPandSigma->sigmad0 ) : -999.;
	dec_z0sinthetaSigTJVA(*track) = (myIPandSigma->sigmaz0SinTheta != 0.) ? (float)( myIPandSigma->IPz0SinTheta / myIPandSigma->sigmaz0SinTheta ) : -999.;
      }
    }
  }

314
315
316
317
  // extrapolate core tracks to calorimeter surface
  // store information only in ExtraDetailsContainer
  if(!m_bypassExtrapolator)
    {
318
      StatusCode sc = extrapolateToCaloSurface(pTau);
319
320
321
      if (sc.isFailure() && !sc.isRecoverable()) {
	ATH_MSG_ERROR("couldn't extrapolate tracks to calo surface");
	return StatusCode::FAILURE;
322
      }
323
324
325
    }
  
  return StatusCode::SUCCESS;
326
327
328
329
330
}


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
TauTrackFinder::TauTrackType TauTrackFinder::tauTrackType( const xAOD::TauJet& pTau,
331
332
							   const xAOD::TrackParticle& trackParticle,
							   const xAOD::Vertex* primaryVertex) const
333
{
334
  double dR = pTau.p4().DeltaR(trackParticle.p4());
335

336
  if (dR > m_maxJetDr_wide) return NotTauTrack;
337

338
339
340
  bool goodTrack = true;
  if(!m_bypassSelector)
    goodTrack = m_trackSelectorTool_tau->decision(trackParticle, primaryVertex);
341
    
342
343
344
345
346
347
348
  if (goodTrack) {
    if (dR > m_maxJetDr_tau)
      return TauTrackWide;
    else
      return TauTrackCore;
  } else
    return TauTrackOther;
349
350
351
352
}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
void TauTrackFinder::getTauTracksFromPV( const xAOD::TauJet& pTau,
353
354
					 const xAOD::TrackParticleContainer& trackParticleCont,
					 const xAOD::Vertex* primaryVertex,
355
					 const bool& useGhostTracks,
356
					 const xAOD::JetContainer* jetContainer,
357
358
359
					 std::vector<const xAOD::TrackParticle*> &tauTracks,
					 std::vector<const xAOD::TrackParticle*> &wideTracks,
					 std::vector<const xAOD::TrackParticle*> &otherTracks) const
360
{
361
362
363
364
365
366
367
368
369
370
371
  std::vector<const xAOD::TrackParticle*> ghostTracks;
  if(useGhostTracks) {
    const xAOD::Jet* seedJet = pTau.jet();
    if(seedJet) {
      if(!seedJet->getAssociatedObjects("GhostTrack", ghostTracks)) {
	ATH_MSG_WARNING("Could not retrieve GhostTrack from seed jet.");
      }
    }
  }
  
  for (const xAOD::TrackParticle *trackParticle : trackParticleCont) {
372
    TauTrackType type = tauTrackType(pTau, *trackParticle, primaryVertex);
373
374
375
376
377
378
    if(type == NotTauTrack) continue;

    if(useGhostTracks) {
      // require that tracks are ghost-matched with the seed jet at large dR(tau,track), to avoid using tracks from another tau
      double dR = pTau.p4().DeltaR(trackParticle->p4());
      if (dR > m_ghostTrackDR) {
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
	if (std::find(ghostTracks.begin(), ghostTracks.end(), trackParticle) == ghostTracks.end()) {
	  // check whether the jet closest to the track is the current seed jet
	  // if so, recover the track even if not ghost-matched, to improve tau-track association efficiency at low pt (esp. for 3p)
	  double dRmin = 999.;
	  bool isSeedClosest = false;
	  for (const xAOD::Jet* jet : *jetContainer) {
	    xAOD::JetFourMom_t jetP4 = jet->jetP4(xAOD::JetScale::JetConstitScaleMomentum);
	    TLorentzVector jetLV;
	    jetLV.SetPtEtaPhiM(jetP4.Pt(), jetP4.Eta(), jetP4.Phi(), jetP4.M());
	    double dRjet = trackParticle->p4().DeltaR(jetLV);
	    if(dRjet < dRmin) {
	      dRmin = dRjet;
	      isSeedClosest = (jet == pTau.jet());
	    }
	  }
	  if(!isSeedClosest) continue;
395
396
397
398
	}
      }
    }
    
399
400
401
402
403
404
405
406
407
408
    if (type == TauTrackCore)
      tauTracks.push_back(trackParticle);
    else if (type == TauTrackWide)
      wideTracks.push_back(trackParticle);
    else if (type == TauTrackOther)
      otherTracks.push_back(trackParticle);
  }
  std::sort(tauTracks.begin(), tauTracks.end(), TrackSort());
  std::sort(wideTracks.begin(), wideTracks.end(), TrackSort());
  std::sort(otherTracks.begin(), otherTracks.end(), TrackSort());
409
410
411
412
}


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
413
StatusCode TauTrackFinder::extrapolateToCaloSurface(xAOD::TauJet& pTau) const {
414

415
416
417
418
419
420
421
  Trk::TrackParametersIdHelper parsIdHelper;

  int trackIndex = -1;
  const Trk::CaloExtension * caloExtension = nullptr;
  std::unique_ptr<Trk::CaloExtension> uniqueExtension;
  for( xAOD::TauTrack* tauTrack : pTau.allTracks() ) {
    const xAOD::TrackParticle *orgTrack = tauTrack->track();
Shaun Roe's avatar
Shaun Roe committed
422
    if( !orgTrack ) continue;
423
424
    trackIndex = orgTrack->index();

Shaun Roe's avatar
Shaun Roe committed
425
    
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444

    // set default values
    float etaEM = -10.0;
    float phiEM = -10.0;
    float etaHad = -10.0;
    float phiHad = -10.0;

    // get the extrapolation into the calo
    ATH_MSG_DEBUG( "Try extrapolation of track with pt = " << orgTrack->pt() 
		   << ", eta " << orgTrack->eta() 
		   << ", phi" << orgTrack->phi() );

    if(!m_ParticleCacheKey.key().empty()){
      /*get the CaloExtension object*/
      ATH_MSG_VERBOSE("Using the CaloExtensionBuilder Cache");
      SG::ReadHandle<CaloExtensionCollection>  particleCache {m_ParticleCacheKey};
      caloExtension = (*particleCache)[trackIndex];
      ATH_MSG_VERBOSE("Getting element " << trackIndex << " from the particleCache");
      if( not caloExtension ){
445
446
447
448
449
        ATH_MSG_VERBOSE("Cache does not contain a calo extension -> "
                        "Calculating with the a CaloExtensionTool");
        uniqueExtension = m_caloExtensionTool->caloExtension(
          Gaudi::Hive::currentContext(), *orgTrack);
        caloExtension = uniqueExtension.get();
450
      }
451
452
    }
    else {
453
454
      /* If CaloExtensionBuilder is unavailable, use the calo extension tool */
      ATH_MSG_VERBOSE("Using the CaloExtensionTool");
455
456
      uniqueExtension = m_caloExtensionTool->caloExtension(
        Gaudi::Hive::currentContext(), *orgTrack);
457
458
459
460
461
462
463
464
      caloExtension = uniqueExtension.get();
    }

    if (!caloExtension)
      { 
	ATH_MSG_DEBUG("Track extrapolation failed");
      }
    else {
465
      const std::vector<Trk::CurvilinearParameters>& clParametersVector = caloExtension->caloLayerIntersections();
466
467
468
469
470
471
472
      if (clParametersVector.empty()) {
	ATH_MSG_DEBUG("Track extrapolation failed");
      }

      ATH_MSG_DEBUG("Scanning samplings");
      bool validECal = false;
      bool validHCal = false;
473
474
      for( const Trk::CurvilinearParameters& cur : clParametersVector ){
	ATH_MSG_DEBUG("Sampling " << parsIdHelper.caloSample(cur.cIdentifier()) );
475
                
476
	// only use entry layer
477
	if( not parsIdHelper.isEntryToVolume(cur.cIdentifier()) ) continue;
478

479
	CaloSampling::CaloSample sample = parsIdHelper.caloSample(cur.cIdentifier());
480
481
482
483
484

	// ECal
	if( not validECal and m_EMSamplings.count(sample))
	  {
	    validECal = true;
485
486
	    etaEM = cur.position().eta();
	    phiEM = cur.position().phi();
487
488
489
490
491
492
493
	    ATH_MSG_DEBUG("Extrapolated to ECal layer " << sample);
	  }

	// HCal
	if( not validHCal and m_HadSamplings.count(sample))
	  {
	    validHCal = true;
494
495
	    etaHad = cur.position().eta();
	    phiHad = cur.position().phi();
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
	    ATH_MSG_DEBUG("Extrapolated to HCal layer " << sample);
	  }
	if( validECal and validHCal ) break;
      }
      // EM failure warn if within acceptance 
      if( not validECal and std::abs(orgTrack->pt()) < 2.48 ){
	ATH_MSG_DEBUG("Failed extrapolation to ECal");
      }
      // Had failure warn if enough pt to reach HCal
      if( not validHCal and orgTrack->pt() > 2000. ){
	ATH_MSG_DEBUG("Failed extrapolation to HCal");
      }

      ATH_MSG_DEBUG( "Extrapolated track with eta=" << orgTrack->eta()
		     << " phi="<<orgTrack->phi()
		     << " to ECal eta=" << etaEM 
		     << " phi="<< phiEM
		     << " HCal eta=" << etaHad 
		     << " phi="<< phiHad
		     );
516
    }
517
518
519
520
521
    tauTrack->setDetail(xAOD::TauJetParameters::CaloSamplingEtaEM, etaEM);
    tauTrack->setDetail(xAOD::TauJetParameters::CaloSamplingPhiEM, phiEM);
    tauTrack->setDetail(xAOD::TauJetParameters::CaloSamplingEtaHad, etaHad);
    tauTrack->setDetail(xAOD::TauJetParameters::CaloSamplingPhiHad, phiHad);
  }
522

523
  return StatusCode::SUCCESS;
524
525
526
527
528

}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
void TauTrackFinder::removeOffsideTracksWrtLeadTrk(std::vector<const xAOD::TrackParticle*> &tauTracks,
529
530
531
532
						   std::vector<const xAOD::TrackParticle*> &wideTracks,
						   std::vector<const xAOD::TrackParticle*> &otherTracks,
						   const xAOD::Vertex* tauOrigin,
						   double maxDeltaZ0) const
533
{
534
  float MAX=1e5;
535

536
  // need at least one core track to have a leading trk to compare with
537
  if (tauTracks.empty()) return;
538

539
540
541
  // get lead trk parameters
  const xAOD::TrackParticle *leadTrack = tauTracks.at(0);
  float z0_leadTrk = getZ0(leadTrack, tauOrigin);
542

543
  if (z0_leadTrk > MAX-1) return; // bad lead trk -> do nothing
544

545
  ATH_MSG_VERBOSE("before z0 cut: #coreTracks=" << tauTracks.size() << ", #wideTracks=" << wideTracks.size() << ", #otherTracks=" << otherTracks.size());
546

547
  std::vector<const xAOD::TrackParticle*>::iterator itr;
548
    
549
550
551
552
553
554
555
556
557
558
559
  // skip leading track, because it is the reference
  itr = tauTracks.begin()+1;
  while (itr!=tauTracks.end()) {
    float z0 = getZ0(*itr, tauOrigin);
    float deltaZ0=z0 - z0_leadTrk;
    ATH_MSG_VERBOSE("core Trks: deltaZ0= " << deltaZ0);

    if ( std::abs(deltaZ0) < maxDeltaZ0 ) {++itr;}
    else {
      if (m_storeInOtherTrks) otherTracks.push_back(*itr);
      itr = tauTracks.erase(itr); //remove from core track collection
560
    }
561
562
563
564
565
566
567
568
  }

  // check wide tracks
  itr = wideTracks.begin();
  while (itr!=wideTracks.end()) {
    float z0 = getZ0(*itr, tauOrigin);
    float deltaZ0=z0 - z0_leadTrk;
    ATH_MSG_VERBOSE("wide Trks: deltaZ0= " << deltaZ0);
569

570
571
572
573
    if ( std::abs(deltaZ0) < maxDeltaZ0 ) { ++itr; }
    else {
      if (m_storeInOtherTrks) otherTracks.push_back(*itr);
      itr = wideTracks.erase(itr); //remove from wide track collection
574
    }
575
  }
576

577
  ATH_MSG_VERBOSE("after z0 cut: #coreTracks=" << tauTracks.size() << ", #wideTracks=" << wideTracks.size() << ", #otherTracks=" << otherTracks.size());
578

579
580
581
582
  // sort again
  std::sort(tauTracks.begin(), tauTracks.end(), TrackSort());
  std::sort(wideTracks.begin(), wideTracks.end(), TrackSort());
  std::sort(otherTracks.begin(), otherTracks.end(), TrackSort());
583
584
585
}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
586
float TauTrackFinder::getZ0(const xAOD::TrackParticle* track, const xAOD::Vertex* vertex) const
587
{
588
  float MAX=1e5;
589

590
  if (!track) return MAX;
591

592
593
594
  const Trk::Perigee* perigee = nullptr;
  if (vertex) perigee = m_trackToVertexTool->perigeeAtVertex(*track, vertex->position());
  else        perigee = m_trackToVertexTool->perigeeAtVertex(*track); //will use beamspot or 0,0,0 instead
595

596
597
598
599
  if (!perigee) {
    ATH_MSG_WARNING("Bad track; can't find perigee at vertex.");
    return MAX;
  }
600

601
  float z0 = perigee->parameters()[Trk::z0];
602

603
  delete perigee; //cleanup necessary to prevent mem leak
604

605
  return z0;
606
}
607
608
609
610
611
612
613
614
615

bool TauTrackFinder::isLargeD0Track(const xAOD::TrackParticle* track) const
{
  const std::bitset<xAOD::NumberOfTrackRecoInfo> patternReco = track->patternRecoInfo();
  if (patternReco.test(xAOD::TrackPatternRecoInfo::SiSpacePointsSeedMaker_LargeD0)) {
    ATH_MSG_DEBUG("LargeD0Track found");
    return true;
  }

616
  return false;
617
}
618
#endif