ATLAS Offline Software
Egamma1_LArStrip_Fex.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 
10 #include "Egamma1_LArStrip_Fex.h"
11 #include "CaloEvent/CaloCell.h"
12 
13 #include "../IO/LArStripNeighborhoodDumper.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>> wrap3(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>({62ul, 63ul, 0ul}));
31  }
32 
33  if (center == 0) {
34  return std::make_optional(std::vector<std::size_t>({63ul, 0ul, 1ul}));
35  }
36 
37  return std::make_optional(std::vector<std::size_t>(
38  {center-1,
39  center,
40  center+1}));
41  }
42 
43  Egamma1_LArStrip_Fex::Egamma1_LArStrip_Fex(const std::string& name, ISvcLocator* pSvcLocator ) :
44  AthReentrantAlgorithm(name, pSvcLocator){
45  }
46 
47 
49  ATH_MSG_INFO ("Initializing " << name());
50 
51  CHECK(m_cellProducer.retrieve());
52  CHECK(m_roiAlgTool.retrieve());
53  CHECK(m_neighKey.initialize());
55 
56 
57  return StatusCode::SUCCESS;
58  }
59 
60 
61  StatusCode Egamma1_LArStrip_Fex::execute(const EventContext& ctx) const {
62  // Read in a CaloCell container. Ask producers to create
63  // vectors of CaloCells to be examined.
64 
65  ATH_MSG_DEBUG ("Executing");
66 
68  if(!eventInfo.isValid()) {
69  ATH_MSG_ERROR ("Error obtaining EventInfo object");
70  return StatusCode::FAILURE;
71  }
72 
73  std::vector<const CaloCell*> cells;
74  CHECK(m_cellProducer->cells(cells, ctx));
75  ATH_MSG_DEBUG(cells.size() <<" cells read in");
76 
77  std::vector<const xAOD::eFexEMRoI*> rois;
78  CHECK(m_roiAlgTool->RoIs(rois, ctx));
79  ATH_MSG_DEBUG(rois.size() << " RoI(s) read in");
80 
81  // find cells in the neighborhood of RoIs.
82  // A neighborhood is a collection of CellData objects which
83  // contain cell eta, phi and Et.
84 
85  auto neighborhoods = std::make_unique<LArStripNeighborhoodContainer>();
86 
87  CHECK(findNeighborhoods(rois, cells, *neighborhoods));
88 
90 
92  if (m_dump) {
93  CHECK(dumper.dump(name(), *eventInfo, *neighborhoods));
94  }
95 
96  if (m_dumpTerse) {
97  CHECK(dumper.dumpTerse(name(), *eventInfo, *neighborhoods));
98  }
99 
100 
101  CHECK(h_write.record(std::move(neighborhoods)));
102 
103  ATH_MSG_INFO(evtStore()->dump());
104 
105  return StatusCode::SUCCESS;
106  }
107 
108  StatusCode
109  Egamma1_LArStrip_Fex::findNeighborhoods(const std::vector<const xAOD::eFexEMRoI*>& rois,
110  const std::vector<const CaloCell*>& cells,
111  LArStripNeighborhoodContainer& neighborhoods) const{
112 
113  for (const auto& roi : rois) {
114  CHECK(findNeighborhood(roi, cells, neighborhoods));
115  }
116 
117  return StatusCode::SUCCESS;
118  }
119 
120 
121  StatusCode
123  const std::vector<const CaloCell*>& cells,
124  LArStripNeighborhoodContainer& neighborhoods) const {
125 
126  // this member function constructs an LArStripNeighborhood.
127  //
128  // A neighourhood is constructed from StripData objects constructed
129  // from CalCells (strips) in the vicinity of an EM RoI the following manner:
130  //
131  // - the cell in the vicinity of the RoI is identified. The neighboorhood
132  // strips are the strups centered on the maximum energy strip.
133  //
134  // In more detail:
135  // A rectangular eta-phi subset of all strips is made.
136  // This subset constains the cells needed for the maximum energy search
137  // and for any subsequent neigborhood strip selection.
138  //
139  // The eta window for the initial selection ismade with half-width of
140  // 0.05 + 8.5* DeltaEta where DeltaEta = 0.0031, the nominal strip width
141  // in eta. A cut-oof off eta = +-1.4 is applied in the selection.
142  //
143  // Phi selection consists of calculating a phi index for the the RoI,
144  // then requiring selected strips to be have a phi indix within +=1 of
145  // the RoI phi index.
146  //
147  // The strips used to identigy the max energy strip have the same phi index
148  // as the RoI, and eta that lies witun 0.1 in eta to the RoI.
149  //
150  // The strips selected for the Neighborhood have eta within +- 8.5 DeltaEta
151  // of the max energy strip, and +- 1 of the RoI strip phi Index.
152 
153 
154 
155 
156  auto cells_near_roi = std::vector<const CaloCell*>();
157 
158  // lambda function to calculate the phi index of a eta-phi point
159  auto phi_ind = [](const auto& c) {
160  constexpr double dphi = std::numbers::pi/32;
161  std::size_t iphi = 32 + int(std::floor(c->phi()/dphi));
162  return iphi > 63 ?
163  std::optional<std::size_t>() : std::make_optional(iphi);
164  };
165 
166  // calculate roi phi index
167  auto roi_phi_index = *phi_ind(roi);
168 
169  // obtain adjacent phi indices
170  auto roi_phi_indices = *wrap3(roi_phi_index);
171  for (const auto& i : roi_phi_indices) {
172  ATH_MSG_DEBUG("roi index " << i);
173  }
174 
175  // container for strips close to RoI in eta
176  auto close = std::vector<std::vector<const CaloCell*>>(3);
177  for (auto& v :close) {v.reserve(100);}
178 
179 
180  constexpr double half_deta_roi{0.05};
181  constexpr double half_deta_neigh{8.5*0.0031};
182  constexpr double half_deta_fid{half_deta_roi + half_deta_neigh};
183 
184 
185  double etalim_low = std::max(roi->eta()-half_deta_fid, -1.4);
186  double etalim_high = std::min(roi->eta()+half_deta_fid, 1.4);
187 
188 
189  // double loop. outer: all strips. inner: adjacent roi indices.
190  for (const auto& cell : cells) {
191  auto icell = *phi_ind(cell);
192  std::size_t pos{0};
193  for(const auto& iroi : roi_phi_indices) {
194  auto c_eta = cell->eta();
195  if (iroi == icell and c_eta >= etalim_low and c_eta < etalim_high) {
196  close[pos].push_back(cell);
197  break;
198  }
199  ++pos;
200  }
201  }
202 
203 
204  for (const auto& v : close) {ATH_MSG_DEBUG("size close " << v.size());}
205 
206  // select the cells within a a tower width of the RoI. Then find the
207  // cell in this selction with the highest energy
208  auto roi_cells = std::vector<const CaloCell*> ();
209  roi_cells.reserve(close[1].size());
210  std::copy_if(std::begin(close[1]),
211  std::end(close[1]),
212  std::back_inserter(roi_cells),
213  [&roi](const auto& c) {
214  return std::abs(c->eta() - roi->eta()) < half_deta_roi;
215  });
216 
217 
218  auto it = std::max_element(std::begin(roi_cells),
219  std::end(roi_cells),
220  [](const auto& l,const auto& r) {
221  return l->e() < r->e();
222  });
223 
224  if (it == std::end(roi_cells)) {
225  ATH_MSG_ERROR("max_cell not found");
226  return StatusCode::FAILURE;
227  }
228 
229  ATH_MSG_DEBUG("max cell pos "
230  << std::distance(std::begin(roi_cells), it)
231  << ' ' << **it);
232 
233  // set up Cell containers for the neighborhood. One container
234  // per adjacent RoI phi indices.
235  auto neigh_cells = std::vector<std::vector<const CaloCell*>>(3);
236 
237  const CaloCell* max_cell{*it};
238  const auto max_cell_eta = max_cell->eta();
239 
240 
241  for (std::size_t iv{0ul}; iv != close.size(); ++iv) {
242  std::copy_if(std::begin(close[iv]),
243  std::end(close[iv]),
244  std::back_inserter(neigh_cells[iv]),
245  [&max_cell_eta, &half_deta_neigh](const auto& c){
246  return abs(c->eta()-max_cell_eta) < half_deta_neigh;
247  });
248  }
249 
250  auto max_neigh_cell_it = std::find(std::begin(neigh_cells[1]),
251  std::end(neigh_cells[1]),
252  max_cell);
253  if (max_neigh_cell_it == std::end(neigh_cells[1])){
254  ATH_MSG_ERROR("Lost the max cell");
255  return StatusCode::FAILURE;
256  }
257 
258  auto max_neigh_cell_pos{std::distance(std::begin(neigh_cells[1]),
259  max_neigh_cell_it)};
260 
261  auto toStripData = [](const auto& fromCells){
262  auto stripdata = std::vector<StripData>();
263  stripdata.reserve(fromCells.size());
264  std::transform(std::begin(fromCells),
265  std::end(fromCells),
266  back_inserter(stripdata),
267  [](const auto& c) {
268  return StripData(c->eta(),
269  c->phi(),
270  c->e());});
271  return stripdata;
272  };
273 
274  auto low = toStripData(neigh_cells[0]);
275  auto center = toStripData(neigh_cells[1]);
276  auto high = toStripData(neigh_cells[2]);
277 
278  Coords roi_c{roi->eta(), roi->phi()};
279  Coords cell_c{max_cell->eta(), max_cell->phi()};
280 
281 
282  neighborhoods.push_back(std::make_unique<LArStripNeighborhood>(low,
283  center,
284  high,
285  roi_c,
286  cell_c,
287  max_neigh_cell_pos));
288 
289  return StatusCode::SUCCESS;
290  }
291 
292 }
beamspotman.r
def r
Definition: beamspotman.py:674
GlobalSim::wrap3
std::optional< std::vector< std::size_t > > wrap3(std::size_t center)
Definition: Egamma1_LArStrip_Fex.cxx:24
RunTileCalibRec.cells
cells
Definition: RunTileCalibRec.py:281
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::m_cellProducer
ToolHandle< ICaloCellsProducer > m_cellProducer
Definition: Egamma1_LArStrip_Fex.h:45
PlotCalibFromCool.begin
begin
Definition: PlotCalibFromCool.py:94
skel.it
it
Definition: skel.GENtoEVGEN.py:407
GlobalSim::Egamma1_LArStrip_Fex::initialize
virtual StatusCode initialize() override
Definition: Egamma1_LArStrip_Fex.cxx:48
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
GlobalSim::Egamma1_LArStrip_Fex::findNeighborhood
StatusCode findNeighborhood(const xAOD::eFexEMRoI *, const std::vector< const CaloCell * > &, LArStripNeighborhoodContainer &) const
Definition: Egamma1_LArStrip_Fex.cxx:122
pi
#define pi
Definition: TileMuonFitter.cxx:65
mergePhysValFiles.end
end
Definition: DataQuality/DataQualityUtils/scripts/mergePhysValFiles.py:92
AthReentrantAlgorithm
An algorithm that can be simultaneously executed in multiple threads.
Definition: AthReentrantAlgorithm.h:74
GlobalSim::Egamma1_LArStrip_Fex::m_neighKey
SG::WriteHandleKey< LArStripNeighborhoodContainer > m_neighKey
Definition: Egamma1_LArStrip_Fex.h:72
python.setupRTTAlg.size
int size
Definition: setupRTTAlg.py:39
AthCommonDataStore< AthCommonMsg< Gaudi::Algorithm > >::evtStore
ServiceHandle< StoreGateSvc > & evtStore()
The standard StoreGateSvc (event store) Returns (kind of) a pointer to the StoreGateSvc.
Definition: AthCommonDataStore.h:85
GlobalSim
AlgTool to obtain a selection of eFex RoIs read in from the event store.
Definition: dump.h:8
GlobalSim::Egamma1_LArStrip_Fex::m_dump
Gaudi::Property< bool > m_dump
Definition: Egamma1_LArStrip_Fex.h:59
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
lumiFormat.i
int i
Definition: lumiFormat.py:85
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::m_eventInfoKey
SG::ReadHandleKey< xAOD::EventInfo > m_eventInfoKey
Definition: Egamma1_LArStrip_Fex.h:37
GlobalSim::Egamma1_LArStrip_Fex::Egamma1_LArStrip_Fex
Egamma1_LArStrip_Fex(const std::string &name, ISvcLocator *pSvcLocator)
Definition: Egamma1_LArStrip_Fex.cxx:43
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
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?
Egamma1_LArStrip_Fex.h
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
GlobalSim::Egamma1_LArStrip_Fex::m_dumpTerse
Gaudi::Property< bool > m_dumpTerse
Definition: Egamma1_LArStrip_Fex.h:65
EventInfo.h
python.LumiBlobConversion.pos
pos
Definition: LumiBlobConversion.py:16
GlobalSim::dump
void dump(const std::string &fn, const T &t)
python.PyAthena.v
v
Definition: PyAthena.py:154
SG::WriteHandle
Definition: StoreGate/StoreGate/WriteHandle.h:73
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::findNeighborhoods
StatusCode findNeighborhoods(const std::vector< const xAOD::eFexEMRoI * > &, const std::vector< const CaloCell * > &, LArStripNeighborhoodContainer &) const
Definition: Egamma1_LArStrip_Fex.cxx:109
xAOD::eFexEMRoI_v1::phi
float phi() const
Seed supercell index within central tower (0 -> 3)
GlobalSim::Egamma1_LArStrip_Fex::m_roiAlgTool
ToolHandle< eFexRoIAlgTool > m_roiAlgTool
Definition: Egamma1_LArStrip_Fex.h:54
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::execute
virtual StatusCode execute(const EventContext &) const override
Definition: Egamma1_LArStrip_Fex.cxx:61