ATLAS Offline Software
TileEopFilterAlg.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 /*
6  * File: TileEopFilterAlg.cxx
7  * Author: Carlos.Solans@cern.ch
8  */
9 
11 #include "TileEopFilterAlg.h"
12 #include "CaloGeoHelpers/CaloSampling.h"
14 
15 using namespace std;
16 
17 //=======================================
18 TileEopFilterAlg::TileEopFilterAlg( const std::string& name, ISvcLocator* pSvcLocator ):
19  AthAlgorithm( name, pSvcLocator ), m_trackInCalo("TrackExtrapolator"){
20 //=======================================
21 
22  declareProperty("InputTracks", m_inputTracks = "InDetTrackParticles");
23  declareProperty("OutputTracks", m_outputTracks = "SelectedTracks");
24  declareProperty("InputClusters", m_inputClusters = "CaloTopoClusters");
25  declareProperty("OutputClusters", m_outputClusters = "SelectedClusters");
26  declareProperty("InputCells", m_inputCells = "AllCalo");
27  declareProperty("OutputCells", m_outputCells = "SelectedCells");
28  declareProperty("TrackClusters", m_trackClusters = "TrackClusters");
29  declareProperty("TrackCells", m_trackCells = "TrackCells");
30  declareProperty("ClusterCells", m_clusterCells = "ClusterCells");
31  declareProperty("TrackMomCut", m_trackMom = 1000.);
32  declareProperty("TrackEtaCut", m_trackEta = 2.4);
33  declareProperty("TrackIsoCut", m_trackIso = 0.8);
34  declareProperty("TrackPtrelCut", m_trackPtrel = 0.15);
35  declareProperty("TrackClusterR", m_trackClusterR = 0.2);
36  declareProperty("TrackCellR", m_trackCellR = 1.0);
37  declareProperty("DumpLArCells", m_dumpLarCells = false);
38  declareProperty("TrackTools", m_trackInCalo);
39 }
40 
41 //=======================================
43 //=======================================
44 
45  CHECK(m_trackInCalo.retrieve());
46  //CHECK(m_trackIsoTool.retrieve());
47 
48  return StatusCode::SUCCESS;
49 }
50 
51 //=======================================
53 //=======================================
54 
55  //Get the input tracks
56  const TRACKCONTAINER* inputTracks = 0;
58 
59  //Allocate the output tracks container
60  TRACKCONTAINER * outputTracks = new TRACKCONTAINER;
61  TRACKAUXCONTAINER* outputAuxTracks = new TRACKAUXCONTAINER;
62  outputTracks->setStore( outputAuxTracks );
63 
64 
65  //Select tracks
66  TRACKCONTAINER::const_iterator trackItr = inputTracks->begin();
68  for(; trackItr != trackEnd; ++trackItr ){
69  const TRACK* track = *trackItr;
70  if(!track){ ATH_MSG_INFO("Not a valid track"); continue; }
71  //Cut 1. Momentum
72  if(track->p4().P() >= m_trackMom){
73  //Cut 2. Eta
74  if(fabs(track->eta()) <= m_trackEta){
75  /*
76  //Cut 3. Isolation: Only one track in that cone
77  //Extrapolate the current track and all other tracks and check they don't overlap in the range
78  int tracksInCone=1;
79  TRACKCONTAINER::const_iterator trackItr2 = inputTracks->begin();
80  TRACKCONTAINER::const_iterator trackEnd2 = inputTracks->end();
81  vector< vector<double> > etrack = m_trackInCalo->getXYZEtaPhiPerLayer(track);
82  for(; trackItr2 != trackEnd2; ++trackItr2 ){
83  if(trackItr==trackItr2) continue;
84  const TRACK* track2 = *trackItr2;
85  if(track2->p() < m_trackMom) continue;
86  vector< vector<double> > etrack2 = m_trackInCalo->getXYZEtaPhiPerLayer(track2);
87  for(int i=0;i<11;i++){
88  double dif_eta = etrack2[i][3] - etrack[i][3];
89  double dif_phi = etrack2[i][4] - etrack[i][4];
90  if(dif_phi<0) dif_phi=-dif_phi;
91  if(dif_phi>M_PI){dif_phi=2*M_PI-dif_phi;}
92  double deltaR = sqrt(dif_eta*dif_eta+dif_phi*dif_phi);
93  if(deltaR<m_trackIso) tracksInCone++;
94  }
95  }
96  if(tracksInCone==1){
97  TRACK* tmpTrack = new TRACK();
98  tmpTrack->makePrivateStore( track );
99  outputTracks->push_back(tmpTrack);
100  }
101  */
102 
103  //Cut 3. Isolation: require ptcone{IsoCut}/track_pt < ptrelCut
104  double ptcone=0.;
105  TRACKCONTAINER::const_iterator trackItr2 = inputTracks->begin();
106  TRACKCONTAINER::const_iterator trackEnd2 = inputTracks->end();
107  for(; trackItr2 != trackEnd2; ++trackItr2 ){
108  if(trackItr==trackItr2) continue;
109  const TRACK* track2 = *trackItr2;
110  double dif_eta = track->eta() - track2->eta();
111  double dif_phi = track->phi() - track2->phi();
112  if(dif_phi<0) dif_phi=-dif_phi;
113  if(dif_phi>M_PI){dif_phi=2*M_PI-dif_phi;}
114  double deltaR = sqrt(dif_eta*dif_eta+dif_phi*dif_phi);
115  if(deltaR<m_trackIso) ptcone += track2->pt();
116  }
117  if(ptcone/track->pt() < m_trackPtrel ){
118  TRACK* tmpTrack = new TRACK();
119  tmpTrack->makePrivateStore( track );
120  outputTracks->push_back(tmpTrack);
121  }
122  }
123  }
124  }
125 
126  ATH_MSG_DEBUG("Number of selected tracks: " << outputTracks->size());
127 
128  //Get input clusters
129  const CLUSTERCONTAINER* inputClusters = 0;
130  CHECK( evtStore()->retrieve( inputClusters, m_inputClusters ) );
131 
132  //Allocate output clusters container
135  outputClusters->setStore( outputAuxClusters );
136 
137 
138  //Allocate output association between tracks and clusters
139  ASSOCCONTAINER* trackClusters = new ASSOCCONTAINER_CONSTRUCTOR(outputTracks->size());
140  CHECK( evtStore()->record( trackClusters, m_trackClusters ) );
141 
142  //Select clusters that match a track
143  //For each cluster loop over all selected tracks
144  //extrapolate the track to each calorimeter layer
145  //compute distance of cluster center in given layer to track pointer in same layer
146  //select cluster if distance is smaller than threshold
147  CLUSTERCONTAINER::const_iterator clusterItr = inputClusters->begin();
148  CLUSTERCONTAINER::const_iterator clusterEnd = inputClusters->end();
149  for(;clusterItr != clusterEnd; ++clusterItr){
150  const CLUSTER* cluster = *clusterItr;
151  TRACKCONTAINER::const_iterator trackItr = outputTracks->begin();
152  TRACKCONTAINER::const_iterator trackEnd = outputTracks->end();
153  ASSOCCONTAINER::iterator assocItr = trackClusters->begin();
154  for( ; trackItr != trackEnd; ++trackItr ){
155  const TRACK* track = *trackItr;
156  vector< vector<double> > etrack = m_trackInCalo->getXYZEtaPhiPerLayer(track);
157  double deltaR=999;
158  for(int cell_sid=0;cell_sid<CaloSampling::Unknown;cell_sid++){
159  int lay=-1;
160  if (cell_sid==0 || cell_sid==4 ){lay=0;}
161  else if(cell_sid==1 || cell_sid==5 ){lay=1;}
162  else if(cell_sid==2 || cell_sid==6 ){lay=2;}
163  else if(cell_sid==3 || cell_sid==7 ){lay=3;}
164  else if(cell_sid==12 || cell_sid==18){lay=4;}
165  else if(cell_sid==13 || cell_sid==19){lay=5;}
166  else if(cell_sid==15 || cell_sid==17){lay=5;}
167  else if(cell_sid==14 || cell_sid==20){lay=6;}
168  else if(cell_sid==16) {lay=6;}
169  else if(cell_sid==8 ){lay=7;}
170  else if(cell_sid==9 ){lay=8;}
171  else if(cell_sid==10 ){lay=9;}
172  else if(cell_sid==11 ){lay=10;}
173  if(lay==-1) continue;
174  double dif_eta = cluster->etaSample((CaloSampling::CaloSample)cell_sid) - etrack[lay][3];
175  double dif_phi = cluster->phiSample((CaloSampling::CaloSample)cell_sid) - etrack[lay][4];
176  if(dif_phi<0) dif_phi=-dif_phi;
177  if(dif_phi>M_PI){dif_phi=2*M_PI-dif_phi;}
178  double tmp = sqrt(dif_eta*dif_eta+dif_phi*dif_phi);
179  deltaR = (tmp<deltaR ? tmp : deltaR);
180  }
181  if(deltaR<=m_trackClusterR){
182  unsigned int i=0;
183  while(i!=outputClusters->size()){if(outputClusters->at(i)==cluster){break;}i++;}
184  if(i==outputClusters->size()){
185  CLUSTER* tmpCluster = new CLUSTER();
186  tmpCluster->makePrivateStore( cluster );
187  outputClusters->push_back(tmpCluster);
188  }
189  //Fill association
190  assocItr->push_back(i); //FIXME
191  }
192  ++assocItr;
193  }
194  }
195 
196  ATH_MSG_DEBUG("Number of selected clusters: " << outputClusters->size());
197 
198  //Get input cells
199  const CELLCONTAINER* inputCells = 0;
200  CHECK( evtStore()->retrieve( inputCells, m_inputCells ) );
201 
202  //Allocate output cells container
204  CHECK( evtStore()->record( outputCells, m_outputCells ) );
205 
206  //Allocate output association between clusters and cells
207  ASSOCCONTAINER* clusterCells = new ASSOCCONTAINER_CONSTRUCTOR(outputClusters->size());
208  CHECK( evtStore()->record( clusterCells, m_clusterCells ) );
209 
210  //Select cells associated to clusters
211  clusterItr = outputClusters->begin();
212  clusterEnd = outputClusters->end();
213  ASSOCCONTAINER::iterator assocItr = clusterCells->begin();
214  //Loop over clusters
215  for(;clusterItr != clusterEnd; ++clusterItr){
216  const CLUSTER* cluster = *clusterItr;
217  if(cluster->getCellLinks() != 0){
218  CLUSTER::const_cell_iterator cellItr = cluster->cell_begin();
219  CLUSTER::const_cell_iterator cellEnd = cluster->cell_end();
220  for(;cellItr != cellEnd; ++cellItr){
221  const CELL* cell = *cellItr;
222  unsigned int i=0;
223  while(i!=outputCells->size()){if(outputCells->at(i)==cell){break;}i++;}
224  if(i==outputCells->size()){
225  outputCells->push_back(cell);
226  }
227  //Fill association
228  assocItr->push_back(i);//FIXME
229  }
230  }
231  ++assocItr;
232  }
233 
234  //Allocate output association between tracks and cells
235  ASSOCCONTAINER* trackCells = new ASSOCCONTAINER_CONSTRUCTOR(outputTracks->size());
236  CHECK( evtStore()->record( trackCells, m_trackCells ) );
237 
238  //Also select cells associated to tracks
239  //cell matches track in a cone of DeltaR and is not in ouputCells already
240  trackItr = outputTracks->begin();
241  trackEnd = outputTracks->end();
242  assocItr = trackCells->begin();
243  //Loop over tracks
244  for( ; trackItr != trackEnd; ++trackItr ){
245  const TRACK* track = *trackItr;
246  vector< vector<double> > etrack = m_trackInCalo->getXYZEtaPhiPerLayer(track);
247  //loop cells
248  CELLCONTAINER::const_iterator cellItr = inputCells->begin();
249  CELLCONTAINER::const_iterator cellEnd = inputCells->end();
250  for(;cellItr != cellEnd; ++cellItr){
251  const CELL* cell = *cellItr;
252  CaloSampling::CaloSample cell_sid = (CaloSampling::CaloSample)cell->caloDDE()->getSampling();
253  //http://acode-browser.usatlas.bnl.gov/lxr/source/atlas/Calorimeter/CaloGeoHelpers/CaloGeoHelpers/CaloSampling.def
254  int lay=-1;
255  if (cell_sid==0 || cell_sid==4 ){lay=0;}
256  else if(cell_sid==1 || cell_sid==5 ){lay=1;}
257  else if(cell_sid==2 || cell_sid==6 ){lay=2;}
258  else if(cell_sid==3 || cell_sid==7 ){lay=3;}
259  else if(cell_sid==12 || cell_sid==18){lay=4;}
260  else if(cell_sid==13 || cell_sid==19){lay=5;}
261  else if(cell_sid==15 || cell_sid==17){lay=5;}
262  else if(cell_sid==14 || cell_sid==20){lay=6;}
263  else if(cell_sid==16) {lay=6;}
264  else if(cell_sid==8 ){lay=7;}
265  else if(cell_sid==9 ){lay=8;}
266  else if(cell_sid==10 ){lay=9;}
267  else if(cell_sid==11 ){lay=10;}
268  if(lay==-1) continue;
269  double dif_eta = cell->eta() - etrack[lay][3];
270  double dif_phi = cell->phi() - etrack[lay][4];
271  if(dif_phi<0) dif_phi=-dif_phi;
272  if(dif_phi>M_PI){dif_phi=2*M_PI-dif_phi;}
273  double deltaR = sqrt(dif_eta*dif_eta+dif_phi*dif_phi);
274  if(deltaR<=m_trackCellR){
275  unsigned int i=0;
276  while(i!=outputCells->size()){if(outputCells->at(i)==cell){break;}i++;}
277  if(i==outputCells->size()){
278  outputCells->push_back(cell);
279  }
280  assocItr->push_back(i);//FIXME
281  }
282  }
283  ++assocItr;
284  }
285 
286  ATH_MSG_DEBUG("Number of selected cells: " << outputCells->size());
287 
288  CHECK( evtStore()->record(outputClusters, m_outputClusters ) );
289  CHECK( evtStore()->record(outputAuxClusters, m_outputClusters+"Aux.") );
290  CHECK( evtStore()->record(outputTracks, m_outputTracks) );
291  CHECK( evtStore()->record(outputAuxTracks, m_outputTracks+"Aux.") );
292 
293 
294  return StatusCode::SUCCESS;
295 
296 }
297 
python.PyKernel.retrieve
def retrieve(aClass, aKey=None)
Definition: PyKernel.py:110
xAOD::iterator
JetConstituentVector::iterator iterator
Definition: JetConstituentVector.cxx:68
xAOD::TrackParticle_v1::pt
virtual double pt() const override final
The transverse momentum ( ) of the particle.
Definition: TrackParticle_v1.cxx:73
GetLCDefs::Unknown
@ Unknown
Definition: GetLCDefs.h:21
ConstDataVector::at
ElementProxy at(size_type n)
Access an element, as an lvalue.
xAOD::CaloCluster_v1::cell_begin
const_cell_iterator cell_begin() const
Iterator of the underlying CaloClusterCellLink (const version)
Definition: CaloCluster_v1.h:812
DataModel_detail::const_iterator
Const iterator class for DataVector/DataList.
Definition: DVLIterator.h:82
xAOD::TrackParticleAuxContainer_v5
Temporary container used until we have I/O for AuxStoreInternal.
Definition: TrackParticleAuxContainer_v5.h:35
CLUSTER
xAOD::CaloCluster CLUSTER
Definition: D3PDMaker/TileD3PDMaker/src/ITrackTools.h:89
ReadCellNoiseFromCool.cell
cell
Definition: ReadCellNoiseFromCool.py:53
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
SG::VIEW_ELEMENTS
@ VIEW_ELEMENTS
this data object is a view, it does not own its elmts
Definition: OwnershipPolicy.h:18
TileEopFilterAlg::m_dumpLarCells
bool m_dumpLarCells
Definition: TileEopFilterAlg.h:51
ConstDataVector.h
DataVector adapter that acts like it holds const pointers.
TileEopFilterAlg::m_trackCells
std::string m_trackCells
Definition: TileEopFilterAlg.h:43
AthCommonDataStore< AthCommonMsg< Algorithm > >::declareProperty
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T > &t)
Definition: AthCommonDataStore.h:145
xAOD::TrackParticle_v1::eta
virtual double eta() const override final
The pseudorapidity ( ) of the particle.
Definition: TrackParticle_v1.cxx:77
TileEopFilterAlg::m_inputTracks
std::string m_inputTracks
Definition: TileEopFilterAlg.h:36
TRACKCONTAINER
xAOD::TrackParticleContainer TRACKCONTAINER
Definition: D3PDMaker/TileD3PDMaker/src/ITrackTools.h:87
M_PI
#define M_PI
Definition: ActiveFraction.h:11
ASSOCCONTAINER_CONSTRUCTOR
#define ASSOCCONTAINER_CONSTRUCTOR(size)
Definition: TileAssocFillerTool.h:23
TileEopFilterAlg::m_trackIso
float m_trackIso
Definition: TileEopFilterAlg.h:47
TileEopFilterAlg::m_trackClusters
std::string m_trackClusters
Definition: TileEopFilterAlg.h:42
xAOD::CaloCluster_v1::etaSample
float etaSample(const CaloSample sampling) const
Retrieve barycenter in a given sample.
Definition: CaloCluster_v1.cxx:532
GeoPrimitives.h
AthCommonDataStore< AthCommonMsg< Algorithm > >::evtStore
ServiceHandle< StoreGateSvc > & evtStore()
The standard StoreGateSvc (event store) Returns (kind of) a pointer to the StoreGateSvc.
Definition: AthCommonDataStore.h:85
TileEopFilterAlg::m_trackInCalo
ToolHandle< ITrackTools > m_trackInCalo
Definition: TileEopFilterAlg.h:53
xAOD::CaloClusterContainer
CaloClusterContainer_v1 CaloClusterContainer
Define the latest version of the calorimeter cluster container.
Definition: Event/xAOD/xAODCaloEvent/xAODCaloEvent/CaloClusterContainer.h:17
xAOD::CaloCluster_v1
Description of a calorimeter cluster.
Definition: CaloCluster_v1.h:59
TileEopFilterAlg::m_clusterCells
std::string m_clusterCells
Definition: TileEopFilterAlg.h:44
python.TrigInDetConfig.inputTracks
inputTracks
Definition: TrigInDetConfig.py:168
lumiFormat.i
int i
Definition: lumiFormat.py:92
CaloSampling::CaloSample
CaloSample
Definition: Calorimeter/CaloGeoHelpers/CaloGeoHelpers/CaloSampling.h:22
TileEopFilterAlg::execute
virtual StatusCode execute()
Definition: TileEopFilterAlg.cxx:52
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
TileEopFilterAlg::m_inputClusters
std::string m_inputClusters
Definition: TileEopFilterAlg.h:38
TileEopFilterAlg::m_trackCellR
float m_trackCellR
Definition: TileEopFilterAlg.h:50
CHECK
#define CHECK(...)
Evaluate an expression and check for errors.
Definition: Control/AthenaKernel/AthenaKernel/errorcheck.h:422
xAOD::Iso::ptcone
@ ptcone
Track isolation.
Definition: IsolationFlavour.h:22
TileEopFilterAlg.h
xAOD::CaloCluster_v1::phiSample
float phiSample(const CaloSample sampling) const
Retrieve barycenter in a given sample.
Definition: CaloCluster_v1.cxx:547
DeMoUpdate.tmp
string tmp
Definition: DeMoUpdate.py:1167
DataVector< xAOD::TrackParticle_v1 >
xAOD::CaloCluster_v1::getCellLinks
const CaloClusterCellLink * getCellLinks() const
Get a pointer to the CaloClusterCellLink object (const version)
Definition: CaloCluster_v1.cxx:905
AthAlgorithm
Definition: AthAlgorithm.h:47
TRACKAUXCONTAINER
xAOD::TrackParticleAuxContainer TRACKAUXCONTAINER
Definition: D3PDMaker/TileD3PDMaker/src/ITrackTools.h:88
TileEopFilterAlg::m_outputCells
std::string m_outputCells
Definition: TileEopFilterAlg.h:41
xAOD::CaloClusterAuxContainer_v2
Auxiliary container for calorimeter cluster containers.
Definition: CaloClusterAuxContainer_v2.h:30
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:195
xAOD::CaloClusterAuxContainer
CaloClusterAuxContainer_v2 CaloClusterAuxContainer
Define the latest version of the calorimeter cluster auxiliary container.
Definition: CaloClusterAuxContainer.h:16
TileEopFilterAlg::m_trackClusterR
float m_trackClusterR
Definition: TileEopFilterAlg.h:49
ConstDataVector::push_back
value_type push_back(value_type pElem)
Add an element to the end of the collection.
DataVector::push_back
value_type push_back(value_type pElem)
Add an element to the end of the collection.
SG::AuxElement::makePrivateStore
void makePrivateStore()
Create a new (empty) private store for this object.
Definition: AuxElement.cxx:172
CaloCellContainer
Container class for CaloCell.
Definition: CaloCellContainer.h:55
DataVector::end
const_iterator end() const noexcept
Return a const_iterator pointing past the end of the collection.
TileEopFilterAlg::initialize
virtual StatusCode initialize()
Definition: TileEopFilterAlg.cxx:42
TileEopFilterAlg::m_inputCells
std::string m_inputCells
Definition: TileEopFilterAlg.h:40
CaloCell
Data object for each calorimeter readout cell.
Definition: CaloCell.h:57
TileEopFilterAlg::m_outputTracks
std::string m_outputTracks
Definition: TileEopFilterAlg.h:37
ConstDataVector
DataVector adapter that acts like it holds const pointers.
Definition: ConstDataVector.h:76
TileEopFilterAlg::m_trackEta
float m_trackEta
Definition: TileEopFilterAlg.h:46
xAOD::CaloCluster_v1::cell_end
const_cell_iterator cell_end() const
Definition: CaloCluster_v1.h:813
TileEopFilterAlg::TileEopFilterAlg
TileEopFilterAlg(const std::string &name, ISvcLocator *pSvcLocator)
Definition: TileEopFilterAlg.cxx:18
TRACK
xAOD::TrackParticle TRACK
Definition: D3PDMaker/TileD3PDMaker/src/ITrackTools.h:86
TileEopFilterAlg::m_trackPtrel
float m_trackPtrel
Definition: TileEopFilterAlg.h:48
xAOD::track
@ track
Definition: TrackingPrimitives.h:512
xAOD::TrackParticle_v1
Class describing a TrackParticle.
Definition: TrackParticle_v1.h:43
DataVector::at
const T * at(size_type n) const
Access an element, as an rvalue.
TileEopFilterAlg::m_trackMom
float m_trackMom
Definition: TileEopFilterAlg.h:45
TileEopFilterAlg::m_outputClusters
std::string m_outputClusters
Definition: TileEopFilterAlg.h:39
makeComparison.deltaR
float deltaR
Definition: makeComparison.py:36
DataVector::size
size_type size() const noexcept
Returns the number of elements in the collection.
DataVector::begin
const_iterator begin() const noexcept
Return a const_iterator pointing at the beginning of the collection.
xAOD::TrackParticle_v1::phi
virtual double phi() const override final
The azimuthal angle ( ) of the particle (has range to .)
ASSOCCONTAINER
std::vector< std::vector< int > > ASSOCCONTAINER
Definition: TileAssocFillerTool.h:22