ATLAS Offline Software
Egamma1_LArStrip_Fex_RowAware.cxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3  */
4 
5 /*
6  This algorithm outputs Calorimeter strip data in the region of
7  eFex RoIs.
8 */
9 
11 #include "../IO/LArStripNeighborhoodDumper.h"
12 
13 #include "CaloEvent/CaloCell.h"
14 
16 
17 #include <fstream>
18 #include <vector>
19 #include <algorithm>
20 #include <optional>
21 
22 namespace GlobalSim {
23 
24  std::optional<std::vector<std::size_t>> wrap5(std::size_t center) {
25  if (center > 63) {
26  return std::optional<std::vector<std::size_t>>{};
27  }
28 
29  if (center == 63) {
30  return std::make_optional(std::vector<std::size_t>({61ul, 62ul, 63ul, 0ul, 1ul}));
31  }
32 
33  if (center == 0) {
34  return std::make_optional(std::vector<std::size_t>({62ul, 63ul, 0ul, 1ul, 2ul}));
35  }
36 
37  return std::make_optional(std::vector<std::size_t>(
38  {center-2,
39  center-1,
40  center,
41  center+1,
42  center+2}));
43  }
44 
45  Egamma1_LArStrip_Fex_RowAware::Egamma1_LArStrip_Fex_RowAware(const std::string& name, ISvcLocator* pSvcLocator ) :
46  AthReentrantAlgorithm(name, pSvcLocator){
47  }
48 
49 
51  ATH_MSG_INFO ("Initializing " << name());
52 
53  CHECK(m_cellProducer.retrieve());
54  CHECK(m_roiAlgTool.retrieve());
55  CHECK(m_neighKey.initialize());
58 
59  return StatusCode::SUCCESS;
60  }
61 
62 
63  StatusCode Egamma1_LArStrip_Fex_RowAware::execute(const EventContext& ctx) const {
64  // Read in a CaloCell container. Ask producers to create
65  // vectors of CaloCells to be examined.
66 
67  ATH_MSG_DEBUG ("Executing");
68 
70  if(!eventInfo.isValid()) {
71  ATH_MSG_ERROR ("Error obtaining EventInfo object");
72  return StatusCode::FAILURE;
73  }
74 
75  std::vector<const CaloCell*> cells;
76  CHECK(m_cellProducer->cells(cells, ctx));
77  ATH_MSG_DEBUG(cells.size() <<"Cells read in");
78 
79  std::vector<const xAOD::eFexEMRoI*> rois;
80  CHECK(m_roiAlgTool->RoIs(rois, ctx));
81  ATH_MSG_DEBUG(rois.size() << " RoI(s) read in");
82 
83  // find cells in the neighborhood of RoIs.
84  // A neighborhood is a collection of CellData objects which
85  // contain cell eta, phi and Et.
86 
87  auto neighborhoods = std::make_unique<LArStripNeighborhoodContainer>();
88  auto phimax = std::make_unique<std::vector<int>>();
89 
90  CHECK(findNeighborhoods_RowAware(rois, cells, *neighborhoods, *phimax));
91 
94 
96  if (m_dump) {
97  CHECK(dumper.dump(name(), *eventInfo, *neighborhoods));
98  }
99 
100  if (m_dumpTerse) {
101  CHECK(dumper.dumpTerse(name(), *eventInfo, *neighborhoods));
102  }
103 
104  CHECK(h_write.record(std::move(neighborhoods)));
105  CHECK(h_phimax.record(std::move(phimax)));
106 
107  return StatusCode::SUCCESS;
108  }
109 
110  StatusCode
111  Egamma1_LArStrip_Fex_RowAware::findNeighborhoods_RowAware(const std::vector<const xAOD::eFexEMRoI*>& rois,
112  const std::vector<const CaloCell*>& cells,
113  LArStripNeighborhoodContainer& neighborhoods,
114  std::vector<int>& phimax) const{
115 
116  for (const auto& roi : rois) {
117  CHECK(findNeighborhood_RowAware(roi, cells, neighborhoods, phimax));
118  }
119 
120  return StatusCode::SUCCESS;
121  }
122 
123 
124  StatusCode
126  const std::vector<const CaloCell*>& cells,
127  LArStripNeighborhoodContainer& neighborhoods,
128  std::vector<int>& phimax) const {
129 
130  // this member function constructs an LArStripNeighborhood.
131  //
132  // A neighourhood is constructed from StripData objects constructed
133  // from CalCells (strips) in the vicinity of an EM RoI the following manner:
134  //
135  // - the cell in the vicinity of the RoI is identified. The neighboorhood
136  // strips are the strups centered on the maximum energy strip.
137  //
138  // In more detail:
139  // A rectangular eta-phi subset of all strips is made.
140  // This subset constains the cells needed for the maximum energy search
141  // and for any subsequent neigborhood strip selection.
142  //
143  // The eta window for the initial selection ismade with half-width of
144  // 0.05 + 8.5* DeltaEta where DeltaEta = 0.003125, the nominal strip width
145  // in eta. A cut-oof off eta = +-1.4 is applied in the selection.
146  //
147  // Phi selection consists of calculating a phi index for the the RoI,
148  // then requiring selected strips to be have a phi indix within +=1 of
149  // the RoI phi index.
150  //
151  // The strips used to identigy the max energy strip the are within +-1 of the
152  // as the phi index of the RoI, and eta that lies within 0.1 in eta to the RoI.
153  //
154  // The strips selected for the Neighborhood have eta within +- 8.5 DeltaEta
155  // of the max energy strip, and +- 1 of the max RoI strip phi Index.
156 
157  auto cells_near_roi = std::vector<const CaloCell*>();
158 
159  ATH_MSG_DEBUG("roi eta " << roi->eta() << " phi " << roi->phi());
160 
161  // lambda function to calculate the phi index of a eta-phi point
162  auto phi_ind = [](const auto& c) {
163  constexpr double dphi = std::numbers::pi/32;
164  std::size_t iphi = 32 + int(std::floor(c->phi()/dphi));
165  return iphi > 63 ?
166  std::optional<std::size_t>() : std::make_optional(iphi);
167  };
168 
169  const std::optional<std::size_t> roi_phi_index_opt = phi_ind(roi);
170  if (not roi_phi_index_opt.has_value()) {
171  return StatusCode::FAILURE;
172  }
173  const std::size_t roi_phi_index = *roi_phi_index_opt;
174 
175  // obtain adjacent phi indices
176  auto roi_phi_indices = *wrap5(roi_phi_index);
177 
178  // container for strips close to RoI in eta
179  auto close = std::deque<std::vector<const CaloCell*>>(5);
180  for (auto& v :close) {v.reserve(100);}
181 
182  //One ROI is 0.1 wide in eta, this is the half width
183  constexpr double half_deta_roi{0.05};
184  //One cell in the barrel is 0.003125, this is the half width of the 17 cell window
185  constexpr double half_deta_neigh{8.5*0.003125};
186  //To define the ±eta limit of our window, we need to sum the above
187  constexpr double half_deta_fid{half_deta_roi + half_deta_neigh};
188 
189 
190  double etalim_low = std::max(roi->eta()-half_deta_fid, -1.4);
191  double etalim_high = std::min(roi->eta()+half_deta_fid, 1.4);
192 
193 
194  // double loop. outer: all strips. inner: adjacent roi indices.
195  for (const auto& cell : cells) {
196  auto icell = *phi_ind(cell);
197  std::size_t pos{0};
198  for(const auto& iroi : roi_phi_indices) {
199  auto c_eta = cell->eta();
200  if (iroi == icell and c_eta >= etalim_low and c_eta < etalim_high) {
201  close[pos].push_back(cell);
202  break;
203  }
204  ++pos;
205  }
206  }
207 
208  // select the cells within a a tower width of the RoI. Then find the
209  // cell in this selction with the highest energy
210  auto roi_cells = std::deque<std::vector<const CaloCell*>>(5);
212 
213  for (std::size_t i{0ul}; i != close.size(); ++i) {
214  roi_cells[i].reserve(close[i].size());
215  std::copy_if(std::begin(close[i]),
216  std::end(close[i]),
217  std::back_inserter(roi_cells[i]),
218  [&roi](const auto& c) {
219  return std::abs(c->eta() - roi->eta()) < half_deta_roi;
220  });
221  //Work out where the max is in the central block.
222  if(i > 0 && i < 4){
223  auto it = std::max_element(std::begin(roi_cells[i]),
224  std::end(roi_cells[i]),
225  [](const auto& l,const auto& r) {
226  return l->e() < r->e();
227  });
228  ATH_MSG_DEBUG("max cell row "
229  << i << " position "
230  << std::distance(std::begin(roi_cells[i]), it)
231  << " :" << **it);
232  roi_max_it.push_back(it);
233  }
234  }
235 
236  auto max_row = std::max_element(std::begin(roi_max_it),
237  std::end(roi_max_it),
238  [](const auto& l,const auto& r) {
239  return (*l)->e() < (*r)->e();
240  });
241 
242  int max_row_pos = std::distance( roi_max_it.begin(), max_row );
243 
244  phimax.push_back(max_row_pos);
245 
246  ATH_MSG_DEBUG("max cell row "
247  << ' ' << std::distance( roi_max_it.begin(), max_row )+1);
248 
249  //If the maximum is not in the centre row, we need to re-seed
250  switch(max_row_pos)
251  {
252  case 0:
253  {
254  close.pop_back();
255  close.pop_back();
256  break;
257  }
258  case 1:
259  {
260  close.pop_front();
261  close.pop_back();
262  break;
263  }
264  case 2:
265  {
266  close.pop_front();
267  close.pop_front();
268  break;
269  }
270  }
271 
272  ATH_MSG_DEBUG("popped ");
273 
274  // set up Cell containers for the neighborhood. One container
275  // per adjacent RoI phi indices.
276  auto neigh_cells = std::vector<std::vector<const CaloCell*>>(3);
277 
278  const CaloCell* max_cell{*(*max_row)};
279  const auto max_cell_eta = max_cell->eta();
280 
281  ATH_MSG_DEBUG("Got the max cell");
282 
283  for (std::size_t iv{0ul}; iv != close.size(); ++iv) {
284  std::copy_if(std::begin(close[iv]),
285  std::end(close[iv]),
286  std::back_inserter(neigh_cells[iv]),
287  [&max_cell_eta, &half_deta_neigh](const auto& c){
288  return abs(c->eta()-max_cell_eta) < half_deta_neigh;
289  });
290  }
291 
292  ATH_MSG_DEBUG("Made our neighbourhood");
293 
294  auto max_neigh_cell_it = std::find(std::begin(neigh_cells[1]),
295  std::end(neigh_cells[1]),
296  max_cell);
297  if (max_neigh_cell_it == std::end(neigh_cells[1])){
298  ATH_MSG_ERROR("Lost the max cell");
299  return StatusCode::FAILURE;
300  }
301 
302  auto max_neigh_cell_pos{std::distance(std::begin(neigh_cells[1]),
303  max_neigh_cell_it)};
304 
305  ATH_MSG_DEBUG("Rediscovered our neighbourhood max");
306 
307  auto toStripData = [](const auto& fromCells){
308  auto stripdata = std::vector<StripData>();
309  stripdata.reserve(fromCells.size());
310  std::transform(std::begin(fromCells),
311  std::end(fromCells),
312  back_inserter(stripdata),
313  [](const auto& c) {
314  return StripData(c->eta(),
315  c->phi(),
316  c->e());});
317  return stripdata;
318  };
319 
320  auto low = toStripData(neigh_cells[0]);
321  auto center = toStripData(neigh_cells[1]);
322  auto high = toStripData(neigh_cells[2]);
323 
324  Coords roi_c{roi->eta(), roi->phi()};
325  Coords cell_c{max_cell->eta(), max_cell->phi()};
326 
327  ATH_MSG_DEBUG("Fill with strip data");
328 
329  neighborhoods.push_back(std::make_unique<LArStripNeighborhood>(low,
330  center,
331  high,
332  roi_c,
333  cell_c,
334  max_neigh_cell_pos));
335 
336  return StatusCode::SUCCESS;
337  }
338 }
xAOD::iterator
JetConstituentVector::iterator iterator
Definition: JetConstituentVector.cxx:68
beamspotman.r
def r
Definition: beamspotman.py:674
GlobalSim::Egamma1_LArStrip_Fex_RowAware::execute
virtual StatusCode execute(const EventContext &) const override
Definition: Egamma1_LArStrip_Fex_RowAware.cxx:63
RunTileCalibRec.cells
cells
Definition: RunTileCalibRec.py:281
GlobalSim::Egamma1_LArStrip_Fex_RowAware::Egamma1_LArStrip_Fex_RowAware
Egamma1_LArStrip_Fex_RowAware(const std::string &name, ISvcLocator *pSvcLocator)
Definition: Egamma1_LArStrip_Fex_RowAware.cxx:45
GlobalSim::StripData
Definition: Global/GlobalSimulation/src/IO/StripData.h:16
ReadCellNoiseFromCool.cell
cell
Definition: ReadCellNoiseFromCool.py:53
WriteCellNoiseToCool.icell
icell
Definition: WriteCellNoiseToCool.py:339
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
find
std::string find(const std::string &s)
return a remapped string
Definition: hcg.cxx:135
SG::ReadHandle
Definition: StoreGate/StoreGate/ReadHandle.h:67
max
constexpr double max()
Definition: ap_fixedTest.cxx:33
min
constexpr double min()
Definition: ap_fixedTest.cxx:26
GlobalSim::Egamma1_LArStrip_Fex_RowAware::m_cellProducer
ToolHandle< ICaloCellsProducer > m_cellProducer
Definition: Egamma1_LArStrip_Fex_RowAware.h:46
Egamma1_LArStrip_Fex_RowAware.h
PlotCalibFromCool.begin
begin
Definition: PlotCalibFromCool.py:94
skel.it
it
Definition: skel.GENtoEVGEN.py:407
GlobalSim::Egamma1_LArStrip_Fex_RowAware::initialize
virtual StatusCode initialize() override
Definition: Egamma1_LArStrip_Fex_RowAware.cxx:50
CaloCell.h
GlobalSim::Coords
std::pair< double, double > Coords
Definition: LArStripNeighborhood.h:26
UploadAMITag.l
list l
Definition: UploadAMITag.larcaf.py:157
xAOD::eFexEMRoI_v1
Class describing a LVL1 eFEX EM region of interest.
Definition: eFexEMRoI_v1.h:33
pi
#define pi
Definition: TileMuonFitter.cxx:65
mergePhysValFiles.end
end
Definition: DataQuality/DataQualityUtils/scripts/mergePhysValFiles.py:92
GlobalSim::Egamma1_LArStrip_Fex_RowAware::m_phimaxKey
SG::WriteHandleKey< std::vector< int > > m_phimaxKey
Definition: Egamma1_LArStrip_Fex_RowAware.h:79
AthReentrantAlgorithm
An algorithm that can be simultaneously executed in multiple threads.
Definition: AthReentrantAlgorithm.h:74
python.setupRTTAlg.size
int size
Definition: setupRTTAlg.py:39
GlobalSim
AlgTool to obtain a selection of eFex RoIs read in from the event store.
Definition: dump.h:8
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
lumiFormat.i
int i
Definition: lumiFormat.py:85
GlobalSim::Egamma1_LArStrip_Fex_RowAware::m_dump
Gaudi::Property< bool > m_dump
Definition: Egamma1_LArStrip_Fex_RowAware.h:59
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
GlobalSim::Egamma1_LArStrip_Fex_RowAware::m_neighKey
SG::WriteHandleKey< LArStripNeighborhoodContainer > m_neighKey
Definition: Egamma1_LArStrip_Fex_RowAware.h:72
Amg::transform
Amg::Vector3D transform(Amg::Vector3D &v, Amg::Transform3D &tr)
Transform a point from a Trasformation3D.
Definition: GeoPrimitivesHelpers.h:156
AnalysisUtils::copy_if
Out copy_if(In first, const In &last, Out res, const Pred &p)
Definition: IFilterUtils.h:30
GlobalSim::Egamma1_LArStrip_Fex_RowAware::m_eventInfoKey
SG::ReadHandleKey< xAOD::EventInfo > m_eventInfoKey
Definition: Egamma1_LArStrip_Fex_RowAware.h:39
CHECK
#define CHECK(...)
Evaluate an expression and check for errors.
Definition: Control/AthenaKernel/AthenaKernel/errorcheck.h:422
detail::ul
unsigned long ul
Definition: PrimitiveHelpers.h:46
SG::VarHandleKey::initialize
StatusCode initialize(bool used=true)
If this object is used as a property, then this should be called during the initialize phase.
Definition: AthToolSupport/AsgDataHandles/Root/VarHandleKey.cxx:103
DataVector
Derived DataVector<T>.
Definition: DataVector.h:794
SG::ReadHandle::isValid
virtual bool isValid() override final
Can the handle be successfully dereferenced?
GlobalSim::Egamma1_LArStrip_Fex_RowAware::findNeighborhood_RowAware
StatusCode findNeighborhood_RowAware(const xAOD::eFexEMRoI *, const std::vector< const CaloCell * > &, LArStripNeighborhoodContainer &, std::vector< int > &) const
Definition: Egamma1_LArStrip_Fex_RowAware.cxx:125
GlobalSim::wrap5
std::optional< std::vector< std::size_t > > wrap5(std::size_t center)
Definition: Egamma1_LArStrip_Fex_RowAware.cxx:24
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:240
DataVector::push_back
value_type push_back(value_type pElem)
Add an element to the end of the collection.
xAOD::eFexEMRoI_v1::eta
float eta() const
setter for the above
EventInfo.h
python.LumiBlobConversion.pos
pos
Definition: LumiBlobConversion.py:16
python.PyAthena.v
v
Definition: PyAthena.py:154
SG::WriteHandle
Definition: StoreGate/StoreGate/WriteHandle.h:73
GlobalSim::Egamma1_LArStrip_Fex_RowAware::m_roiAlgTool
ToolHandle< eFexRoIAlgTool > m_roiAlgTool
Definition: Egamma1_LArStrip_Fex_RowAware.h:54
python.CaloAddPedShiftConfig.int
int
Definition: CaloAddPedShiftConfig.py:45
SG::WriteHandle::record
StatusCode record(std::unique_ptr< T > data)
Record a const object to the store.
CaloCell
Data object for each calorimeter readout cell.
Definition: CaloCell.h:57
GlobalSim::LArStripNeighborhoodDumper
Definition: LArStripNeighborhoodDumper.h:14
TriggerTest.rois
rois
Definition: TriggerTest.py:23
GlobalSim::Egamma1_LArStrip_Fex_RowAware::m_dumpTerse
Gaudi::Property< bool > m_dumpTerse
Definition: Egamma1_LArStrip_Fex_RowAware.h:65
xAOD::eFexEMRoI_v1::phi
float phi() const
Seed supercell index within central tower (0 -> 3)
Amg::distance
float distance(const Amg::Vector3D &p1, const Amg::Vector3D &p2)
calculates the distance between two point in 3D space
Definition: GeoPrimitivesHelpers.h:54
python.compressB64.c
def c
Definition: compressB64.py:93
GlobalSim::Egamma1_LArStrip_Fex_RowAware::findNeighborhoods_RowAware
StatusCode findNeighborhoods_RowAware(const std::vector< const xAOD::eFexEMRoI * > &, const std::vector< const CaloCell * > &, LArStripNeighborhoodContainer &, std::vector< int > &) const
Definition: Egamma1_LArStrip_Fex_RowAware.cxx:111
CaloCell::eta
virtual double eta() const override final
get eta (through CaloDetDescrElement)
Definition: CaloCell.h:376