TauTrackFinder.cxx 24.5 KB
Newer Older
1
/*
2
  Copyright (C) 2002-2022 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
  ATH_CHECK( m_beamSpotKey.initialize(inTrigger()) );
52

53
  return StatusCode::SUCCESS;
54
55
56
}

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

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

73
  //Retrieve LRT container
74
75
76
77
78
79
80
81
  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");
    }
82
83
84
    else { 
      largeD0TracksParticleCont = trackPartInHandle.cptr();
    }  
85
86
  }

87
88
89
90
91
92
93
94
95
96
97
  // 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();
  }

98
  // get the primary vertex
99
  const xAOD::Vertex* pVertex = pTau.vertex();
100
101
102

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

104
  getTauTracksFromPV(pTau, *trackParticleCont, pVertex, m_useGhostTracks, jetContainer, tauTracks, wideTracks, otherTracks);
105

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

113
114
115
116
117
118
119
120
121
122
123
  // 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	
124
	for( const xAOD::TauTrack* tau_trk : tauTrackCon ) {
125
126
	  //originally it was coreTrack&passTrkSelector
	  if(! tau_trk->flagWithMask( (1<<xAOD::TauJetParameters::TauTrackFlag::coreTrack) | (1<<xAOD::TauJetParameters::TauTrackFlag::passTrkSelector))) continue; 
127
128
129
130
	  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 );
131
	if (alreadyUsed) track_it = tauTracks.erase(track_it);
132
133
134
	else ++track_it;
      }
  }
135

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

141
    ATH_MSG_VERBOSE(name() << " adding core track nr: " << i
142
		    << " eta " << trackParticle->eta()
143
		    << " phi " << trackParticle->phi());
144

145
    charge += trackParticle->charge();
146

147
    xAOD::TauTrack* track = new xAOD::TauTrack();
148
    tauTrackCon.push_back(track);
149

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

161
162
163
164
165
166
    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);
167
    
168
    ElementLink<xAOD::TauTrackContainer> linkToTauTrack;
169
    linkToTauTrack.toContainedElement(tauTrackCon, track);
170
171
    pTau.addTauTrackLink(linkToTauTrack);

172
    ATH_MSG_VERBOSE(name() << " added core track nr: " << i
173
		    << " eta " << pTau.track(i)->eta()
174
		    << " phi " << pTau.track(i)->phi());
175
176
177
178
179
180
181
  }
  // 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);

182
    ATH_MSG_VERBOSE(name() << " adding wide track nr: " << i
183
		    << " eta " << trackParticle->eta()
184
		    << " phi " << trackParticle->phi());
185
186

    xAOD::TauTrack* track = new xAOD::TauTrack();
187
    tauTrackCon.push_back(track);
188
189

    ElementLink<xAOD::TrackParticleContainer> linkToTrackParticle;
190
191
    if (foundLRTCont && isLargeD0Track(trackParticle)){//Check LRT track and link to container
      linkToTrackParticle.toContainedElement(*largeD0TracksParticleCont, trackParticle);
192
      track->setFlag(xAOD::TauJetParameters::TauTrackFlag::LargeRadiusTrack, true);
193
194
    }
    else {
195
      linkToTrackParticle.toContainedElement(*trackParticleCont, trackParticle);
196
      track->setFlag(xAOD::TauJetParameters::TauTrackFlag::LargeRadiusTrack, false);
197
    }
198
199
200
201
202
203
204
205
206
207
208
    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;
209
    linkToTauTrack.toContainedElement(tauTrackCon, track);
210
211
212
213
214
215
216
217
218
219
220
221
    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()
222
		    << " phi " << trackParticle->phi());
223
224

    xAOD::TauTrack* track = new xAOD::TauTrack();
225
    tauTrackCon.push_back(track);
226
227

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

249
250
251
252
  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());

253
254
255
  ATH_MSG_DEBUG("numTrack: " << "/" << pTau.nTracks());
  ATH_MSG_DEBUG("charge: " << "/" << pTau.charge());

256
257
258
  // impact parameter variables w.r.t. tau vertex 
  const xAOD::Vertex* vxcand = nullptr;

259
260
261
  xAOD::Vertex vxbkp;
  vxbkp.makePrivateStore();

262
263
  // FIXME: don't we want to use the beamspot in the offline reconstruction too, when there is no reconstructed primary vertex?
  if (pTau.vertex()!=nullptr && pTau.vertex()->vertexType() != xAOD::VxType::NoVtx) {
264
    vxcand = pTau.vertex();
265
266
  }
  else if (inTrigger()) { // online: use vertex with x-y coordinates from the beamspot and the z from the leading track
267
    vxbkp.setX(0); vxbkp.setY(0); vxbkp.setZ(0);
268
269
270

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

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

285

286
287
288
289
290
  // 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");
291

292
  for(auto *track : pTau.allTracks()) {      
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
    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.;
      }
    }
  }

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


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

334
  if (dR > m_maxJetDr_wide) return NotTauTrack;
335

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

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
void TauTrackFinder::getTauTracksFromPV( const xAOD::TauJet& pTau,
351
352
					 const xAOD::TrackParticleContainer& trackParticleCont,
					 const xAOD::Vertex* primaryVertex,
353
					 const bool& useGhostTracks,
354
					 const xAOD::JetContainer* jetContainer,
355
356
357
					 std::vector<const xAOD::TrackParticle*> &tauTracks,
					 std::vector<const xAOD::TrackParticle*> &wideTracks,
					 std::vector<const xAOD::TrackParticle*> &otherTracks) const
358
{
359
360
361
362
363
364
365
366
367
368
369
  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) {
370
    TauTrackType type = tauTrackType(pTau, *trackParticle, primaryVertex);
371
372
373
374
375
376
    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) {
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
	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;
393
394
395
396
	}
      }
    }
    
397
398
399
400
401
402
403
404
405
406
    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());
407
408
409
410
}


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
411
StatusCode TauTrackFinder::extrapolateToCaloSurface(xAOD::TauJet& pTau) const {
412

413
414
415
416
417
418
419
  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
420
    if( !orgTrack ) continue;
421
    trackIndex = orgTrack->index();   
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440

    // 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 ){
441
442
443
444
445
        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();
446
      }
447
448
    }
    else {
449
450
      /* If CaloExtensionBuilder is unavailable, use the calo extension tool */
      ATH_MSG_VERBOSE("Using the CaloExtensionTool");
451
452
      uniqueExtension = m_caloExtensionTool->caloExtension(
        Gaudi::Hive::currentContext(), *orgTrack);
453
454
455
456
457
458
459
460
      caloExtension = uniqueExtension.get();
    }

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

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

475
	CaloSampling::CaloSample sample = parsIdHelper.caloSample(cur.cIdentifier());
476
477
478
479
480

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

	// HCal
	if( not validHCal and m_HadSamplings.count(sample))
	  {
	    validHCal = true;
490
491
	    etaHad = cur.position().eta();
	    phiHad = cur.position().phi();
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
	    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
		     );
512
    }
513
514
515
516
517
    tauTrack->setDetail(xAOD::TauJetParameters::CaloSamplingEtaEM, etaEM);
    tauTrack->setDetail(xAOD::TauJetParameters::CaloSamplingPhiEM, phiEM);
    tauTrack->setDetail(xAOD::TauJetParameters::CaloSamplingEtaHad, etaHad);
    tauTrack->setDetail(xAOD::TauJetParameters::CaloSamplingPhiHad, phiHad);
  }
518

519
  return StatusCode::SUCCESS;
520
521
522
523
524

}

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

532
  // need at least one core track to have a leading trk to compare with
533
  if (tauTracks.empty()) return;
534

535
536
537
  // get lead trk parameters
  const xAOD::TrackParticle *leadTrack = tauTracks.at(0);
  float z0_leadTrk = getZ0(leadTrack, tauOrigin);
538

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

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

543
  std::vector<const xAOD::TrackParticle*>::iterator itr;
544
    
545
546
547
548
549
550
551
552
553
554
555
  // 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
556
    }
557
558
559
560
561
562
563
564
  }

  // 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);
565

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

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

575
576
577
578
  // sort again
  std::sort(tauTracks.begin(), tauTracks.end(), TrackSort());
  std::sort(wideTracks.begin(), wideTracks.end(), TrackSort());
  std::sort(otherTracks.begin(), otherTracks.end(), TrackSort());
579
580
581
}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
582
float TauTrackFinder::getZ0(const xAOD::TrackParticle* track, const xAOD::Vertex* vertex) const
583
{
584
  float MAX=1e5;
585

586
  if (!track) return MAX;
587

588
589
590
  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
591

592
593
594
595
  if (!perigee) {
    ATH_MSG_WARNING("Bad track; can't find perigee at vertex.");
    return MAX;
  }
596

597
  float z0 = perigee->parameters()[Trk::z0];
598

599
  delete perigee; //cleanup necessary to prevent mem leak
600

601
  return z0;
602
}
603
604
605
606
607
608
609
610
611

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;
  }

612
  return false;
613
}
614
#endif