PFRecoverSplitShowersTool.cxx 14 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
21
*/

#include "eflowRec/PFRecoverSplitShowersTool.h"
#include "eflowRec/eflowUtil.h"

#include "eflowRec/eflowRecCluster.h"
#include "eflowRec/eflowRecTrack.h"
#include "eflowRec/eflowTrackClusterLink.h"
#include "eflowRec/eflowCaloObject.h"
#include "eflowRec/eflowCaloObjectMaker.h"
#include "eflowRec/eflowTrackCaloPoints.h"
#include "eflowRec/eflowCellList.h"
#include "eflowRec/IEFlowCellEOverPTool.h"
#include "eflowRec/PFTrackClusterMatchingTool.h"
#include "eflowRec/eflowEEtaBinnedParameters.h"
#include "eflowRec/eflowLayerIntegrator.h"
#include "eflowRec/eflowRingSubtractionManager.h"
#include "eflowRec/eflowCellSubtractionFacilitator.h"
#include "eflowRec/eflowSubtractor.h"
22
#include "eflowRec/PFSubtractionStatusSetter.h"
23
#include "eflowRec/PFSubtractionEnergyRatioCalculator.h"
24
25
26
27
28
29
30

#include "CaloEvent/CaloClusterContainer.h"
#include "xAODCaloEvent/CaloClusterKineHelper.h"

using namespace eflowSubtract;

PFRecoverSplitShowersTool::PFRecoverSplitShowersTool(const std::string& type,const std::string& name,const IInterface* parent):
31
  base_class(type, name, parent),
32
  m_binnedParameters(std::make_unique<eflowEEtaBinnedParameters>())
33
34
35
36
37
38
39
40
{
}

PFRecoverSplitShowersTool::~PFRecoverSplitShowersTool() {}

StatusCode PFRecoverSplitShowersTool::initialize(){

  if (m_theEOverPTool.retrieve().isFailure()){
41
    ATH_MSG_WARNING("Cannot find eflowEOverPTool");
42
43
44
    return StatusCode::SUCCESS;
  }

45
  if (m_theEOverPTool->fillBinnedParameters(m_binnedParameters.get()).isFailure()){
46
    ATH_MSG_WARNING("Could not execute eflowCellEOverPTool");
47
48
49
50
51
52
    return StatusCode::SUCCESS;
  }

  return StatusCode::SUCCESS;
}

53
void PFRecoverSplitShowersTool::execute(eflowCaloObjectContainer* theEflowCaloObjectContainer, eflowRecTrackContainer* eflowRecTracks, eflowRecClusterContainer* eflowRecClusters) const {
54

55
  ATH_MSG_DEBUG("Executing");
56

57
58
  eflowData data;
  data.caloObjects = theEflowCaloObjectContainer;
59
60
  data.tracksToRecover.reserve(eflowRecTracks->size());
  data.clustersToConsider.reserve(eflowRecClusters->size());
61

62
63
  fillTracksToRecover(data);
  fillClustersToConsider(data);
64
65

  /* Add each track to its matching cluster's eflowCaloObject */
66
  int const nOriginalObj = matchAndCreateEflowCaloObj(data);
67
68

  /* Extrapolate tracks, match clusters, do shower simulation, run cell subtraction */
69
  performRecovery(nOriginalObj,data);
70
71
72
73
74
75
76
77

}

StatusCode PFRecoverSplitShowersTool::finalize(){
  return StatusCode::SUCCESS;

}

78
void PFRecoverSplitShowersTool::fillClustersToConsider(eflowData& data) const {
79

80
  data.clustersToConsider.clear();
81

82
  for (auto thisEflowCaloObject : *data.caloObjects){
83
84
85
86
87
88
89
90
91

    if (thisEflowCaloObject->nClusters() == 0) { continue; }

    for(unsigned i=0; i<thisEflowCaloObject->nClusters(); ++i){
        /* Skip empty clusters (subtraction remnants) */
        const CaloClusterCellLink* theCellLink = thisEflowCaloObject->efRecCluster(i)->getCluster()->getCellLinks();
        if (0 == (int)theCellLink->size()){ continue; }

        thisEflowCaloObject->efRecCluster(i)->clearTrackMatches();
92
        data.clustersToConsider.insert(thisEflowCaloObject->efRecCluster(i));
93
94
95
96
97
        thisEflowCaloObject->clearClusters();
    }
  }
}

98
void PFRecoverSplitShowersTool::fillTracksToRecover(eflowData& data) const {
99

100
101
  data.tracksToRecover.clear(); 
  for (auto thisEflowCaloObject : *data.caloObjects){
102
103
104
105

    /* Skip isolated tracks if flag set */
    if (!m_recoverIsolatedTracks && thisEflowCaloObject->nClusters() == 0) {
      unsigned int nTrk = thisEflowCaloObject->nTracks();
106
      // But make sure we get eflowObjects from them
107
      for (unsigned int iTrk = 0; iTrk < nTrk; ++iTrk) {
108
109
	      eflowRecTrack* thisEfRecTrack = thisEflowCaloObject->efRecTrack(iTrk);
      	if (!thisEfRecTrack->isSubtracted()) thisEfRecTrack->setSubtracted();
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
      }
      continue;
    }

    /* Add all tracks on the eflowCaloObject that haven't been subtracted yet*/
    std::vector<eflowRecTrack*> updatedTracks; updatedTracks.clear();
    unsigned int nTracks = thisEflowCaloObject->nTracks();

    /* For cluster-only CaloObjects */
    if(nTracks == 0) continue;

    /* For track-only and track-cluster CaloObjects */
    for (unsigned int iTrack = 0; iTrack < nTracks; ++iTrack){
      eflowRecTrack* thisEfRecTrack = thisEflowCaloObject->efRecTrack(iTrack);
      if (thisEfRecTrack->isSubtracted()){
        updatedTracks.push_back(thisEfRecTrack);
        continue;
      }
      thisEfRecTrack->clearClusterMatches();
129
      data.tracksToRecover.push_back(thisEfRecTrack);
130
131
132
133
134
135
136
137
138
139
    }
    thisEflowCaloObject->clearTracks();
    if (!updatedTracks.empty()) {
      thisEflowCaloObject->addTracks(updatedTracks);
      updatedTracks.clear();
    } else {
      thisEflowCaloObject->clearLinks();
    }
  }

140
  std::sort(data.tracksToRecover.begin(),data.tracksToRecover.end(),eflowRecTrack::SortDescendingPt());
141
142
143

}

144
unsigned int PFRecoverSplitShowersTool::matchAndCreateEflowCaloObj(eflowData& data) const {
scott snyder's avatar
scott snyder committed
145
146
147

  const EventContext& ctx = Gaudi::Hive::currentContext();

148
  /* Cache the original number of eflowCaloObjects */
149
  const unsigned int nCaloObj = data.caloObjects->size();
150

151
  /* loop tracks in data.tracksToRecover and do matching */
152
  for (auto *thisEfRecTrack : data.tracksToRecover){
153
 
154
155
156
157
    /* No point to do anything if bin does not exist */
    if (!thisEfRecTrack->hasBin()) {
      std::unique_ptr<eflowCaloObject> thisEflowCaloObject = std::make_unique<eflowCaloObject>();
      thisEflowCaloObject->addTrack(thisEfRecTrack);
158
      data.caloObjects->push_back(std::move(thisEflowCaloObject));
159
160
      continue;
    }
161
    if (msgLvl(MSG::DEBUG)){
162
      const xAOD::TrackParticle* track = thisEfRecTrack->getTrack();
163
      ATH_MSG_DEBUG("Recovering charged EFO with e,pt, eta and phi " << track->e() << ", " << track->pt() << ", " << track->eta() << " and " << track->phi());
164
    }
165
    // Get list of matched clusters in the dR<0.2 cone -- already identified
166
167
    const std::vector<eflowTrackClusterLink*>* matchedClusters_02 = thisEfRecTrack->getAlternativeClusterMatches("cone_02");
    if (!matchedClusters_02) { continue; }
168
169

    if (msgLvl(MSG::DEBUG)){
170
      for (auto *trkClusLink : *matchedClusters_02) {
171
	const eflowRecCluster* thisEFRecCluster = trkClusLink->getCluster();
Mark Hodgkinson's avatar
Mark Hodgkinson committed
172
	ATH_MSG_DEBUG("Have matched cluster with e, pt, eta, phi of " << thisEFRecCluster->getCluster()->e() << ", " <<  thisEFRecCluster->getCluster()->pt() << ", " << thisEFRecCluster->getCluster()->eta() << " and " << thisEFRecCluster->getCluster()->phi());
173
      }
174
175
    }

176
    if (matchedClusters_02->empty()) { continue; }
177
178

    /* Matched cluster: create TrackClusterLink and add it to both the track and the cluster (eflowCaloObject will be created later) */
179
    for (auto *trkClusLink : *matchedClusters_02){
180
181
      eflowRecCluster* thisEFRecCluster = trkClusLink->getCluster();
      // Look up whether this cluster is intended for recovery
182
      if( data.clustersToConsider.find(trkClusLink->getCluster()) == data.clustersToConsider.end() ) {continue;}
scott snyder's avatar
scott snyder committed
183
      eflowTrackClusterLink* trackClusterLink = eflowTrackClusterLink::getInstance(thisEfRecTrack,thisEFRecCluster, ctx);
184
      thisEfRecTrack->addClusterMatch(trackClusterLink);
185
      thisEFRecCluster->addTrackMatch(trackClusterLink);
186
187
188
189
190
    }
  }

  /* Create all eflowCaloObjects that contain eflowRecCluster */
  eflowCaloObjectMaker makeCaloObject;
191
192
193
  std::vector<eflowRecCluster*> v_clustersToConsider(data.clustersToConsider.begin(),data.clustersToConsider.end());
  std::sort(v_clustersToConsider.begin(),v_clustersToConsider.end(),eflowRecCluster::SortDescendingPt());
  unsigned int nCaloObjects = makeCaloObject.makeTrkCluCaloObjects(data.tracksToRecover, v_clustersToConsider,
194
								   data.caloObjects);
195
  ATH_MSG_DEBUG("PFRecoverSplitShowersTool created " << nCaloObjects << " CaloObjects");
196

197
198
199
200
201
202
203
  const double gaussianRadius = 0.032;
  const double gaussianRadiusError = 1.0e-3;
  const double maximumRadiusSigma = 3.0;
  
  // Should move to a common header or similar, as shared with PFCellLevelSubtractionTool
  eflowLayerIntegrator integrator(gaussianRadius, gaussianRadiusError, maximumRadiusSigma, m_isHLLHC);

204
  /* integrate cells; determine FLI; eoverp */
205
206
207
  for (unsigned int iCalo = nCaloObj; iCalo < data.caloObjects->size(); ++iCalo) {
    eflowCaloObject* thisEflowCaloObject = data.caloObjects->at(iCalo);
    thisEflowCaloObject->simulateShower(&integrator, m_binnedParameters.get(), m_useUpdated2015ChargedShowerSubtraction);
208
209
210
211
  }
  return nCaloObj;
}

212
void PFRecoverSplitShowersTool::performSubtraction(eflowCaloObject* thisEflowCaloObject) const {
213
214

  PFSubtractionStatusSetter pfSubtractionStatusSetter;
215
216
217
218
219
  PFSubtractionEnergyRatioCalculator pfSubtractionEnergyRatioCalculator;
  if( msgLevel( MSG::DEBUG ) ) {
    pfSubtractionStatusSetter.setLevel(MSG::DEBUG);
    pfSubtractionEnergyRatioCalculator.setLevel(MSG::DEBUG);
  }
220

221
  for (unsigned iTrack = 0; iTrack < thisEflowCaloObject->nTracks(); ++iTrack) {
222
223
224
    eflowRecTrack* thisEfRecTrack = thisEflowCaloObject->efRecTrack(iTrack);
    unsigned int trackIndex = thisEfRecTrack->getTrack()->index();

225
226
    ATH_MSG_DEBUG("About to recover track with e, pt, eta and phi of " << thisEfRecTrack->getTrack()->e() << ", " << thisEfRecTrack->getTrack()->pt() << ", " << thisEfRecTrack->getTrack()->eta() << " and "
    << thisEfRecTrack->getTrack()->eta());
227
228
229
230
    /* Get matched cluster via Links */
    std::vector<eflowRecCluster*> matchedClusters;
    matchedClusters.clear();
    std::vector<eflowTrackClusterLink*> links = thisEfRecTrack->getClusterMatches();
231
    for ( auto *thisEFlowTrackClusterLink : links) matchedClusters.push_back(thisEFlowTrackClusterLink->getCluster());
Teng Jian Khoo's avatar
Teng Jian Khoo committed
232
    std::sort(matchedClusters.begin(),matchedClusters.end(),eflowRecCluster::SortDescendingPt());
233

234
    if (msgLvl(MSG::DEBUG)){
235
      for (auto *thisClus : matchedClusters) ATH_MSG_DEBUG("Cluster " << thisClus->getCluster()->index() << " with e,pt, eta and phi of " << thisClus->getCluster()->e() << ", "<< thisClus->getCluster()->pt() << ", " << thisClus->getCluster()->eta() << " and " << thisClus->getCluster()->phi() << " will be subtracted");
236
    }
237
    /* Do subtraction */
238
    std::vector<std::pair<xAOD::CaloCluster*, bool> > clusterSubtractionList;
239
    clusterSubtractionList.reserve(matchedClusters.size());
240
    std::map<xAOD::CaloCluster*, double> clusterEnergyMap;
241
    for (auto *thisEFlowRecCluster : matchedClusters) {
242
      xAOD::CaloCluster* thisCluster = thisEFlowRecCluster->getCluster();
243
      clusterSubtractionList.emplace_back(thisCluster,false);
244
245
246
      clusterEnergyMap[thisCluster] = thisCluster->e();
    }

247
248
    ATH_MSG_DEBUG("Have filled clusterSubtractionList for this eflowCaloObject");

249
250
251
    if (getSumEnergy(clusterSubtractionList) - thisEfRecTrack->getEExpect() < m_subtractionSigmaCut
        * sqrt(thisEfRecTrack->getVarEExpect())) {
      /* Check if we can annihilate right away */
252
253
      if( msgLevel( MSG::DEBUG ) ) for (auto thisPair : clusterSubtractionList) ATH_MSG_DEBUG("Annihilating cluster with E and eta " << thisPair.first->e() << " and " << thisPair.first->eta());
    
254
      Subtractor::annihilateClusters(clusterSubtractionList);
255
256
      //Now we should mark all of these clusters as being subtracted
      //Now need to mark which clusters were modified in the subtraction procedure
257
258
259
      std::vector<float> clusterSubtractedEnergyRatios;
      pfSubtractionEnergyRatioCalculator.calculateSubtractedEnergyRatiosForAnnih(clusterSubtractionList,clusterEnergyMap,clusterSubtractedEnergyRatios);       
      pfSubtractionStatusSetter.markSubtractionStatus(clusterSubtractionList, clusterSubtractedEnergyRatios, *thisEflowCaloObject,trackIndex);	
260
261
    } else {
      /* Subtract the track from all matched clusters */
262
      const bool debugToggle = msgLvl(MSG::DEBUG);
263
264
      Subtractor::subtractTracksFromClusters(thisEfRecTrack, clusterSubtractionList, debugToggle);                

265
266
267
      ATH_MSG_DEBUG("Have performed subtraction for this eflowCaloObject");

      /* Annihilate the cluster(s) if the remnant is small (i.e. below k*sigma) */      
268
269
      if (getSumEnergy(clusterSubtractionList) < m_subtractionSigmaCut
          * sqrt(thisEfRecTrack->getVarEExpect())) {
270
        if( msgLevel( MSG::DEBUG ) ) for (auto thisPair : clusterSubtractionList) ATH_MSG_DEBUG("Annihilating remnant cluster with E and eta " << thisPair.first->e() << " and " << thisPair.first->eta());
271
        Subtractor::annihilateClusters(clusterSubtractionList);
272
	      //Now we should mark all of these clusters as being subtracted
273
274
275
        std::vector<float> clusterSubtractedEnergyRatios;
        pfSubtractionEnergyRatioCalculator.calculateSubtractedEnergyRatiosForAnnih(clusterSubtractionList,clusterEnergyMap,clusterSubtractedEnergyRatios); 	      
        pfSubtractionStatusSetter.markSubtractionStatus(clusterSubtractionList, clusterSubtractedEnergyRatios, *thisEflowCaloObject,trackIndex);	
276
      }
277
      else{
278
279
280
281
        std::vector<float> clusterSubtractedEnergyRatios;        
        pfSubtractionEnergyRatioCalculator.calculateSubtractedEnergyRatios(clusterSubtractionList,clusterEnergyMap,clusterSubtractedEnergyRatios);    
	      pfSubtractionStatusSetter.markSubtractionStatus(clusterSubtractionList, clusterSubtractedEnergyRatios, *thisEflowCaloObject,trackIndex);	
      } 
282

283
    }
284
285

    if (msgLvl(MSG::DEBUG)){
286
      for (auto *thisClus : matchedClusters) ATH_MSG_DEBUG("Cluster with e,pt, eta and phi of " << thisClus->getCluster()->e() << ", "<< thisClus->getCluster()->pt() << ", " << thisClus->getCluster()->eta() << " and " << thisClus->getCluster()->phi() << " has been subtracted");
287
288
    } 

289
    /* Flag tracks as subtracted */
290
    if (!thisEfRecTrack->isSubtracted()) thisEfRecTrack->setSubtracted();
291
292
293
  }
}

294
295
void PFRecoverSplitShowersTool::performRecovery(unsigned int const nOriginalObj, eflowData& data) const {
  unsigned int nEFCaloObs = data.caloObjects->size();
296
  for (unsigned int iCalo = nOriginalObj; iCalo < nEFCaloObs; ++iCalo) {
297
    eflowCaloObject* thisEflowCaloObject = data.caloObjects->at(iCalo);
298
    performSubtraction(thisEflowCaloObject);
299
300
301
302
  }

}

303
double PFRecoverSplitShowersTool::getSumEnergy(const std::vector<std::pair<xAOD::CaloCluster*, bool> >& clusters) const {
304
  double result = 0.0;
305
  for (auto thisPair : clusters) result += (thisPair.first)->e();
306
307
  return result;
}