PFCellLevelSubtractionTool.cxx 20.8 KB
Newer Older
1
/*
scott snyder's avatar
scott snyder committed
2
  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
*/

#include "eflowRec/PFCellLevelSubtractionTool.h"

#include "eflowRec/eflowRecCluster.h"
#include "eflowRec/eflowRecTrack.h"
#include "eflowRec/eflowTrackClusterLink.h"
#include "eflowRec/eflowDatabase.h"
#include "eflowRec/eflowDepthCalculator.h"
#include "eflowRec/eflowTrackCaloPoints.h"
#include "eflowRec/IEFlowCellEOverPTool.h"
#include "eflowRec/eflowEEtaBinnedParameters.h"
#include "eflowRec/eflowLayerIntegrator.h"
#include "eflowRec/PFTrackClusterMatchingTool.h"
#include "eflowRec/eflowRingSubtractionManager.h"
#include "eflowRec/eflowCellSubtractionFacilitator.h"
#include "eflowRec/eflowSubtractor.h"
#include "eflowRec/eflowRingThicknesses.h"
21
#include "eflowRec/PFSubtractionStatusSetter.h"
22
#include "eflowRec/PFSubtractionEnergyRatioCalculator.h"
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

#include "CaloEvent/CaloCluster.h"

#include "xAODCaloEvent/CaloCluster.h"
#include "xAODCaloEvent/CaloClusterKineHelper.h"

#include "xAODPFlow/PFO.h"

#include "eflowRec/eflowCaloObjectMaker.h"
#include "GaudiKernel/MsgStream.h"

#include "GaudiKernel/ListItem.h"

#include <algorithm>
#include <iostream>
#include <cmath>
#include <vector>
40
#include <map>
41
42
43
44
45


using namespace eflowSubtract;

PFCellLevelSubtractionTool::PFCellLevelSubtractionTool(const std::string& type,const std::string& name,const IInterface* parent) :
46
  base_class( type, name, parent),
47
  m_binnedParameters(nullptr)
48
{
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
}

PFCellLevelSubtractionTool::~PFCellLevelSubtractionTool() {
}

StatusCode PFCellLevelSubtractionTool::initialize(){
  StatusCode sc;

  if (m_matchingTool.retrieve().isFailure()) {
    msg(MSG::WARNING) << "Cannot find PFTrackClusterMatchingTool" << endmsg;
  }

  if (m_matchingToolForPull_02.retrieve().isFailure()) {
    msg(MSG::WARNING) << "Cannot find PFTrackClusterMatchingTool_2" << endmsg;
  }

  m_binnedParameters = std::make_unique<eflowEEtaBinnedParameters>();
66

67
  sc = m_theEOverPTool->fillBinnedParameters(m_binnedParameters.get());
68

69
70
  m_runInGoldenMode = ((m_goldenModeString == "golden1") || (m_goldenModeString == "golden2"));

71
72
73
74
75
  if (sc.isFailure()) {
    msg(MSG::WARNING) << "Could not execute eflowCellEOverPTool " << endmsg;
    return StatusCode::SUCCESS;
  }

76
77
78
79
80
81
  m_trkpos.reset( dynamic_cast<PFMatch::TrackEtaPhiInFixedLayersProvider*>(PFMatch::TrackPositionFactory::Get("EM2EtaPhi").release()) );
  if(!m_trkpos.get()) {
    ATH_MSG_ERROR("Failed to get TrackPositionProvider for cluster preselection!");
    return StatusCode::FAILURE;
  }

82
83
84
85
  return StatusCode::SUCCESS;

}

86
void PFCellLevelSubtractionTool::execute(eflowCaloObjectContainer* theEflowCaloObjectContainer, eflowRecTrackContainer* recTrackContainer, eflowRecClusterContainer* recClusterContainer) const {
87
88
89

  ATH_MSG_VERBOSE("Executing PFCellLevelSubtractionTool");

90
91
92
93
  eflowData data;
  data.caloObjects = theEflowCaloObjectContainer;
  data.tracks = recTrackContainer;
  data.clusters = recClusterContainer;
94
  data.clusterLUT.fill(*recClusterContainer);
95

96
  /* Add each track to its best matching cluster's eflowCaloObject */
97
  matchAndCreateEflowCaloObj(m_nMatchesInCellLevelSubtraction, data);
98
99

  if (msgLvl(MSG::DEBUG)) printAllClusters(*recClusterContainer);
100

101
  /* Check e/p mode - only perform subtraction if not in this mode */
102
  if (!m_calcEOverP) {performSubtraction(data);}
103
104

  ATH_MSG_DEBUG("Have executed performSubtraction");
105

106
  /* Check e/p mode - only perform radial profiles calculations if in this mode */
107
  if (m_calcEOverP) {calculateRadialEnergyProfiles(data);}
108

109
  ATH_MSG_DEBUG("Have executed calculateRadialEnergyProfiles");
110

111
112
}

113
114
unsigned int PFCellLevelSubtractionTool::matchAndCreateEflowCaloObj(unsigned int n, eflowData& data) const {
  unsigned int nMatches(0);
115

scott snyder's avatar
scott snyder committed
116
117
  const EventContext& ctx = Gaudi::Hive::currentContext();

118
  /* Create eflowTrackClusterLink after matching */
119
  unsigned int nTrack = data.tracks->size();
120
  for (unsigned int iTrack=0; iTrack<nTrack; ++iTrack) {
121
    eflowRecTrack *thisEfRecTrack = static_cast<eflowRecTrack*>(data.tracks->at(iTrack));
122

123
124
125
126
127
128
129
130
    // Preselect clusters in a local cone
    eflowRecMatchTrack matchTrack(thisEfRecTrack);
    PFMatch::EtaPhi etaphi = m_trkpos->getPosition(&matchTrack);
    float tracketa = etaphi.getEta();
    float trackphi = etaphi.getPhiD();
    float dRpresel = 0.4; // Some safety margin, but still << full detector
    std::vector<eflowRecCluster*> nearbyClusters = data.clusterLUT.clustersInCone(tracketa,trackphi,dRpresel);

131
     /* Add cluster matches needed for pull calculation*/
132
133
134
135
136
    std::vector<std::pair<eflowRecCluster*,float> > bestClusters_02 = m_matchingToolForPull_02->doMatches(thisEfRecTrack, data.clusters, -1);

    for (auto& matchpair : bestClusters_02) {
      eflowRecCluster* theCluster = matchpair.first;
      float distancesq = matchpair.second;
scott snyder's avatar
scott snyder committed
137
      eflowTrackClusterLink* trackClusterLink = eflowTrackClusterLink::getInstance(thisEfRecTrack, theCluster, ctx);
138
139
140
141
142
      if(distancesq<0.15*0.15) {
	// Narrower cone is a subset of the selected clusters
	// Distance returned is deltaR^2
	thisEfRecTrack->addAlternativeClusterMatch(trackClusterLink,"cone_015");
      }
143
      thisEfRecTrack->addAlternativeClusterMatch(trackClusterLink,"cone_02");
144
145
    }

146
    std::vector<std::pair<eflowRecCluster*,float> > bestClusters = m_matchingTool->doMatches(thisEfRecTrack, data.clusters, n);
147
148
149
150
    if (bestClusters.empty()) { continue; }

    /* Matched cluster: create TrackClusterLink and add it to both the track and the cluster (eflowCaloObject will be created later) */
    nMatches++;
151
152
    for (auto& matchpair : bestClusters) {
      eflowRecCluster* theCluster = matchpair.first;
scott snyder's avatar
scott snyder committed
153
      eflowTrackClusterLink* trackClusterLink = eflowTrackClusterLink::getInstance(thisEfRecTrack, theCluster, ctx);
154
155
      thisEfRecTrack->addClusterMatch(trackClusterLink);
      theCluster->addTrackMatch(trackClusterLink);
156
    }
157

158
159
160
161
  }

  /* Create 3 types eflowCaloObjects: track-only, cluster-only, track-cluster-link */
  eflowCaloObjectMaker makeCaloObject;
162
  unsigned int nCaloObjects = makeCaloObject.makeTrkCluCaloObjects(data.tracks, data.clusters, data.caloObjects);
163
  ATH_MSG_DEBUG("PFCellLevelSubtractionTool created total " << nCaloObjects << " CaloObjects.");
164

165
166
167
168
169
170
171
  // Should move to a common header or similar, as shared with PFRecoverSplitShowersTool
  const double gaussianRadius = 0.032;
  const double gaussianRadiusError = 1.0e-3;
  const double maximumRadiusSigma = 3.0;

  eflowLayerIntegrator integrator(gaussianRadius, gaussianRadiusError, maximumRadiusSigma, m_isHLLHC);

172
  /* integrate cells; determine FLI; eoverp */
173
174
  for (unsigned int iCalo=0; iCalo<data.caloObjects->size(); ++iCalo) {
    eflowCaloObject *thisEflowCaloObject = static_cast<eflowCaloObject*>(data.caloObjects->at(iCalo));
175

176
    thisEflowCaloObject->simulateShower(&integrator, m_binnedParameters.get(), m_useUpdated2015ChargedShowerSubtraction);
177
178
179
180
181
  }

  return nMatches;
}

182
void PFCellLevelSubtractionTool::calculateRadialEnergyProfiles(eflowData& data) const {
183

184
  ATH_MSG_DEBUG("Accessed radial energy profile function");
185

186
  unsigned int nEFCaloObs = data.caloObjects->size();
187

188
  for (unsigned int iEFCalOb = 0; iEFCalOb < nEFCaloObs; ++iEFCalOb) {
189

190
    eflowCaloObject* thisEflowCaloObject = data.caloObjects->at(iEFCalOb);
191

192
    unsigned int nClusters = thisEflowCaloObject->nClusters();
193

194
195
196
    if (nClusters < 1) {
    continue;
    }
197

198
    const std::vector<std::pair<eflowTrackClusterLink*,float> >& matchedTrackList = thisEflowCaloObject->efRecLink();
199

200
    unsigned int nTrackMatches = thisEflowCaloObject->nTracks();
201

202
    for (unsigned int iTrack = 0; iTrack < nTrackMatches; ++iTrack) {
203

204
      eflowRecTrack* efRecTrack = (matchedTrackList[iTrack].first)->getTrack();
205

206
207
208
      std::vector<eflowRecCluster*> matchedClusters;
      matchedClusters.clear();
      std::vector<eflowTrackClusterLink*> links = efRecTrack->getClusterMatches();
209
      for (auto *thisEFlowTrackClusterLink : links) matchedClusters.push_back(thisEFlowTrackClusterLink->getCluster());
210

211
      std::vector<std::pair<xAOD::CaloCluster*, bool> > clusterSubtractionList;
212
213
      clusterSubtractionList.reserve(matchedClusters.size());

214
      for (auto *thisEFlowRecCluster : matchedClusters) clusterSubtractionList.emplace_back(thisEFlowRecCluster->getCluster(),false);
215
216
217
218
219
220

      eflowCellList calorimeterCellList;
      Subtractor::makeOrderedCellList(efRecTrack->getTrackCaloPoints(),clusterSubtractionList,calorimeterCellList);
      std::vector<int> layerToStoreVector;
      std::vector<float> radiusToStoreVector;
      std::vector<float> avgEdensityToStoreVector;
221

222
      for (int i=0; i < eflowCalo::nRegions ;i++){
223

224
	eflowCaloENUM layer = (eflowCaloENUM)i;
225
	ATH_MSG_DEBUG("layer is "<<layer);
226
	double ringThickness = eflowRingThicknesses::ringThickness((eflowCaloENUM)i);
227
	ATH_MSG_DEBUG("ring thickness is "<<ringThickness);
228

229
	double eta_extr = calorimeterCellList.etaFF(layer);
230
	ATH_MSG_DEBUG("extrapolated eta ["<<layer<<"] is "<<eta_extr);
231
	double phi_extr = calorimeterCellList.phiFF(layer);
232
	ATH_MSG_DEBUG("extrapolated phi ["<<layer<<"] is "<<phi_extr);
233

234
235
236
	if (eta_extr == -999.0){
	  continue;
	}
237

238
239
	int indexofRing = 0;
	int layerToStore = -999;
240

241
242
	double radiusToStore = 0;
	double totalEnergyPerCell = 0;
243

244
245
	double energyDensityPerCell = -666;
	double totalEnergyPerRing = 0;
246

247
248
	double cellVolume = 0;
	int totalCells = 0;
249

250
	/* 100 is chosen as a number higher than the number of cells found in a normal list */
251
	bool breakloop = false;
252
	for (unsigned int n=1; n<100; n++){
253

254
255
256
257
258
	  CellIt beginRing = calorimeterCellList.getLowerBound((eflowCaloENUM)i, ringThickness*(n-1));

	  if(beginRing == calorimeterCellList.end()){
	    break;
	  }
259

260
	  indexofRing += 1;
261

262
263
264
	  int totalCellsinRing = 0;
	  double energyDensityPerRing = 0;
	  double averageEnergyDensityPerRing = 0;
265

266
267
	  std::vector<std::pair<const CaloCell*,int> > tempVector = (*beginRing).second;

268
269
	  for (auto thisPair : tempVector){
	    const CaloDetDescrElement* DDE = (thisPair.first)->caloDDE();
270
	    CaloCell_ID::CaloSample sampling = DDE->getSampling();
271
272
273
274
275

	    if(eflowCalo::translateSampl(sampling)!=(eflowCaloENUM)i){
	      breakloop = true;
	      break;
	    }
276

277
278
	    ATH_MSG_DEBUG(" cell eta and phi are " << (thisPair.first)->eta() << " and " << (thisPair.first)->phi() << " with index " << thisPair.second << " and sampling of " << sampling);
	    ATH_MSG_DEBUG(" cell energy is " << (thisPair.first)->energy());
279

280
281
	    totalCells += 1;
	    totalCellsinRing += 1;
282

283
284
	    totalEnergyPerRing += (thisPair.first)->energy();
	    totalEnergyPerCell = (thisPair.first)->energy();
285
286
	    ATH_MSG_DEBUG(" Total E per Cell is " << totalEnergyPerCell);
	    ATH_MSG_DEBUG(" Total E per Ring is " << totalEnergyPerRing);
287

288
	    cellVolume = DDE->volume();
289
	    ATH_MSG_DEBUG(" cell volume is " << cellVolume/1000.);
290

291
	    energyDensityPerCell = totalEnergyPerCell/(cellVolume/1000.);
292
293
	    ATH_MSG_DEBUG(" E density per Cell is " << energyDensityPerCell);
	    ATH_MSG_DEBUG(" Initial added E density per Cell is " << energyDensityPerRing);
294
	    energyDensityPerRing += energyDensityPerCell;
295
	    ATH_MSG_DEBUG(" Final added E density per Cell is " << energyDensityPerRing);
296
297
	    averageEnergyDensityPerRing = energyDensityPerRing/((totalCellsinRing)*(efRecTrack->getTrack()->e()/1000.));
	  }
298

299
	  ATH_MSG_DEBUG(" track eta is " << efRecTrack->getTrack()->eta());
300
301
	  ATH_MSG_DEBUG(" track E is " << efRecTrack->getTrack()->e()/1000.);
	  ATH_MSG_DEBUG(" Average E density per Ring is " << averageEnergyDensityPerRing);
302

303
304
305
	  if (averageEnergyDensityPerRing != 0){
	    avgEdensityToStoreVector.push_back(averageEnergyDensityPerRing);
	    layerToStore = (eflowCaloENUM)i;
306
	    ATH_MSG_DEBUG("layerToStore is "<< layerToStore);
307
308
	    layerToStoreVector.push_back(layerToStore);
	    radiusToStore = (indexofRing)*ringThickness;
309
	    ATH_MSG_DEBUG("radiusToStore is "<< radiusToStore);
310
311
	    radiusToStoreVector.push_back(radiusToStore);
	  }
312
	  else {ATH_MSG_DEBUG("averageEnergyDensityPerRing = 0");}
313
	  if (breakloop) break;
314
315
	}//loop on 100 cells
      }//loop on calo regions
316

317
318
319
      efRecTrack->setLayerCellOrderVector(layerToStoreVector);
      efRecTrack->setRadiusCellOrderVector(radiusToStoreVector);
      efRecTrack->setAvgEDensityCellOrderVector(avgEdensityToStoreVector);
320

321
322
323
      layerToStoreVector.clear();
      radiusToStoreVector.clear();
      avgEdensityToStoreVector.clear();
324
    }//loop on track matches
325
326
327
  }//loop on eflowCaloObjects
}

328
void PFCellLevelSubtractionTool::performSubtraction(eflowData& data) const {
329

330
  ATH_MSG_DEBUG("In performSubtraction");
331
332

  PFSubtractionStatusSetter pfSubtractionStatusSetter;
333
334
335
336
337
  PFSubtractionEnergyRatioCalculator pfSubtractionEnergyRatioCalculator;
  if( msgLevel( MSG::DEBUG ) ) {
    pfSubtractionStatusSetter.setLevel(MSG::DEBUG);
    pfSubtractionEnergyRatioCalculator.setLevel(MSG::DEBUG);
  }
338

339
  unsigned int nEFCaloObs = data.caloObjects->size();
340
  for (unsigned int iEFCalOb = 0; iEFCalOb < nEFCaloObs; ++iEFCalOb) {
341
    eflowCaloObject* thisEflowCaloObject = data.caloObjects->at(iEFCalOb);
342

343
    unsigned int nClusters = thisEflowCaloObject->nClusters();
344

345
    ATH_MSG_DEBUG("Have got an eflowCaloObject with " << nClusters << " clusters");
346

347
348
349
    if (nClusters < 1) {
      continue;
    }
350

351
352
353
354
355
    /* Get matched cluster via Links */

    double expectedEnergy = thisEflowCaloObject->getExpectedEnergy();
    double clusterEnergy = thisEflowCaloObject->getClusterEnergy();
    double expectedSigma = sqrt(thisEflowCaloObject->getExpectedVariance());
356

357
    ATH_MSG_DEBUG("Have got reference values for eflowCaloObject");
358

359
360
    /* Check e/p */
    if (isEOverPFail(expectedEnergy, expectedSigma, clusterEnergy, m_consistencySigmaCut,
361
362
                     m_runInGoldenMode)
        || (m_runInGoldenMode && thisEflowCaloObject->nTracks() > 1)) {
363
364
      continue;
    }
365

366
    ATH_MSG_DEBUG("eflowCaloObject passed eOverP checks");
367

368
369
    /* Do subtraction */

370
371
    unsigned int nTrackMatches = thisEflowCaloObject->nTracks();
    if (nTrackMatches == 0 || m_runInGoldenMode) {
372
373
374
      continue;
    }

375
    const std::vector<std::pair<eflowTrackClusterLink*, float> >& matchedTrackList = thisEflowCaloObject->efRecLink();
376

377
378
379
380
381
382
383
384
385
    if( msgLevel( MSG::DEBUG ) ){ 
      for (unsigned int iTrack = 0; iTrack < nTrackMatches; ++iTrack) {
	      const xAOD::TrackParticle* thisTrack = (matchedTrackList[iTrack].first)->getTrack()->getTrack();
        ATH_MSG_DEBUG("eflowCaloObject has track with E, pt and eta " << thisTrack->e() << ", " << thisTrack->pt() << " and " << thisTrack->eta());
      }
    }

    ATH_MSG_DEBUG("About to perform subtraction for this eflowCaloObject");

386
    /* Check if we can annihilate right away */
387
    if (canAnnihilated(expectedEnergy, expectedSigma, clusterEnergy)) {
388

389
      std::vector<std::pair<xAOD::CaloCluster*, bool> > clusterList;
390
      std::map<xAOD::CaloCluster*, double> clusterEnergyMap;
391
      unsigned nCluster = thisEflowCaloObject->nClusters();
392
      for (unsigned iCluster = 0; iCluster < nCluster; ++iCluster) clusterList.emplace_back(thisEflowCaloObject->efRecCluster(iCluster)->getCluster(),false);
393

394
395
396
      ATH_MSG_DEBUG("We are going to annihilate. ExpectedEnergy, expectedSigma and clusterEnergy are " << expectedEnergy << ", " << expectedSigma << " and " << clusterEnergy);
      if( msgLevel( MSG::DEBUG ) ) for (auto thisPair : clusterList) ATH_MSG_DEBUG("Annihilating cluster with E and eta " << thisPair.first->e() << " and " << thisPair.first->eta());

397
      Subtractor::annihilateClusters(clusterList);
398

399
400
      //Now we should mark all of these clusters as being subtracted      
      pfSubtractionStatusSetter.markAllTracksAnnihStatus(*thisEflowCaloObject);
401

402
    } else {
403

404
      /* Subtract the track from all matched clusters */
405

406
407
      ATH_MSG_DEBUG("Subtract tracks one by one");

408
      for (unsigned int iTrack = 0; iTrack < nTrackMatches; ++iTrack) {
409
	      eflowRecTrack* efRecTrack = (matchedTrackList[iTrack].first)->getTrack();
410
        unsigned int trackIndex = efRecTrack->getTrack()->index();
411

412
	      ATH_MSG_DEBUG("Have got eflowRecTrack number " << iTrack << " for this eflowCaloObject");
413

414
415
416
417
	      /* Can't subtract without e/p */
	      if (!efRecTrack->hasBin()) {
	        continue;
	      }
418

419
	      if (efRecTrack->isInDenseEnvironment()) continue;
420

421
	      ATH_MSG_DEBUG("Have bin and am not in dense environment for this eflowCaloObject");
422

423
424
425
	      std::vector<eflowRecCluster*> matchedClusters;
	      matchedClusters.clear();
	      std::vector<eflowTrackClusterLink*> links = efRecTrack->getClusterMatches();
426
	      for (auto *thisEFlowTrackClusterLink : links) matchedClusters.push_back(thisEFlowTrackClusterLink->getCluster());
427

428
	      ATH_MSG_DEBUG("Have filled matchedClusters list for this eflowCaloObject");
429

430
431
	      std::vector<std::pair<xAOD::CaloCluster*, bool> > clusterSubtractionList;
        std::map<xAOD::CaloCluster*, double> clusterEnergyMap;
432
	      for (auto *thisEFlowRecCluster : matchedClusters) {
433
          xAOD::CaloCluster* thisCluster = thisEFlowRecCluster->getCluster();
434
          clusterSubtractionList.emplace_back(thisCluster,false);
435
436
437
438
          clusterEnergyMap[thisCluster] = thisCluster->e();
        }

	      ATH_MSG_DEBUG("Have filled clusterSubtractionList for this eflowCaloObject");
439

440
441
442
        bool debugToggle = false;
        if( msgLevel( MSG::DEBUG ) ) debugToggle = true;
	      Subtractor::subtractTracksFromClusters(efRecTrack, clusterSubtractionList,debugToggle);
443

444
	      ATH_MSG_DEBUG("Have performed subtraction for this eflowCaloObject");
445

446
447
	      /* Annihilate the cluster(s) if the remnant is small (i.e. below k*sigma) */
	      if (canAnnihilated(0, expectedSigma, clusterEnergy)) {
448
          if( msgLevel( MSG::DEBUG ) ) for (auto thisPair : clusterSubtractionList) ATH_MSG_DEBUG("Annihilating cluster with E and eta " << thisPair.first->e() << " and " << thisPair.first->eta());
449
450
	        Subtractor::annihilateClusters(clusterSubtractionList);
	        //Now we should mark all of these clusters as being subtracted
451
452
453
          std::vector<float> clusterSubtractedEnergyRatios;
          pfSubtractionEnergyRatioCalculator.calculateSubtractedEnergyRatiosForAnnih(clusterSubtractionList,clusterEnergyMap,clusterSubtractedEnergyRatios);           
          pfSubtractionStatusSetter.markSubtractionStatus(clusterSubtractionList, clusterSubtractedEnergyRatios, *thisEflowCaloObject, trackIndex);
454
455
456
	      }
        else{
          std::vector<float> clusterSubtractedEnergyRatios;
457
458
          pfSubtractionEnergyRatioCalculator.calculateSubtractedEnergyRatios(clusterSubtractionList,clusterEnergyMap,clusterSubtractedEnergyRatios);          
	        pfSubtractionStatusSetter.markSubtractionStatus(clusterSubtractionList, clusterSubtractedEnergyRatios, *thisEflowCaloObject, trackIndex);
459
460
        }

461
	      ATH_MSG_DEBUG("Have checked if can annihilate clusters for this eflowCaloOject");
462

463
      }
464
    }
465

466
467
468
469
470
471
472
473
474
475
476
477
478
    /* Flag tracks as subtracted */
    for (unsigned int iTrack = 0; iTrack < thisEflowCaloObject->nTracks(); ++iTrack) {
      if (!thisEflowCaloObject->efRecTrack(iTrack)->isInDenseEnvironment()) thisEflowCaloObject->efRecTrack(iTrack)->setSubtracted();
    }
  }
}

StatusCode PFCellLevelSubtractionTool::finalize(){

  return StatusCode::SUCCESS;
}


479
bool PFCellLevelSubtractionTool::isEOverPFail(double expectedEnergy, double sigma, double clusterEnergy, bool consistencySigmaCut, bool useGoldenMode) const {
480
481
482
483
484
485
486
487
488
  if ((expectedEnergy == 0) && (clusterEnergy > 0)) {
    return false;
  }

  bool result = useGoldenMode ? fabs(clusterEnergy - expectedEnergy) > consistencySigmaCut*sigma
                              : clusterEnergy < expectedEnergy - consistencySigmaCut*sigma;
  return result;
}

489
bool PFCellLevelSubtractionTool::canAnnihilated(double expectedEnergy, double sigma, double clusterEnergy) const {
490
491
492
   return  clusterEnergy - expectedEnergy < m_subtractionSigmaCut * sigma;
}

493
std::string PFCellLevelSubtractionTool::printTrack(const xAOD::TrackParticle* track) const {
494
495
496
497
498
  std::stringstream result;
  result << " track with E, eta and phi "<< track->e() << ", " << track->eta() << " and " << track->phi();
  return result.str();
}

499
std::string PFCellLevelSubtractionTool::printCluster(const xAOD::CaloCluster* cluster) const {
500
501
502
503
504
  std::stringstream result;
  result << " cluster with E, eta and phi of " << cluster->e() << ", " << cluster->eta() << " and " << cluster->phi();
  return result.str();
}

505
void PFCellLevelSubtractionTool::printAllClusters(const eflowRecClusterContainer& recClusterContainer) const {
506
507
508
509
510
511
512
513
  unsigned int nClusters = recClusterContainer.size();
  for (unsigned int iCluster = 0; iCluster < nClusters; ++iCluster) {
    /* Create the eflowCaloObject and put it in the container */
    const eflowRecCluster* thisEFRecCluster = recClusterContainer.at(iCluster);
    if (thisEFRecCluster->getTrackMatches().empty()) {
      ATH_MSG_DEBUG("Isolated" << printCluster(thisEFRecCluster->getCluster()));
    } else {
      ATH_MSG_DEBUG("Matched" << printCluster(thisEFRecCluster->getCluster()));
514
      unsigned int nTrackMatches = thisEFRecCluster->getNTracks();
515
      std::vector<eflowTrackClusterLink*> theTrackLinks = thisEFRecCluster->getTrackMatches();
516
      for (unsigned int iTrackMatch = 0; iTrackMatch < nTrackMatches; ++iTrackMatch) {
517
518
519
520
521
	ATH_MSG_DEBUG("Matched" << printTrack(theTrackLinks[iTrackMatch]->getTrack()->getTrack()));
      }
    }
  }
}