ATLAS Offline Software
CaloCellLinkContainerCnv_p2.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 /* @file CaloCellLinkContainerCnv_p2.cxx
6  * @author Ilija Vukotic <ivukotic@cern.ch>, scott snyder <snyder@bnl.gov>
7  * @date May, 2007
8  * @brief T/P conversions for CaloCellLinkContainerCnv_p2.
9  */
10 
11 
12 #include "AthLinks/ElementLink.h"
18 
19 #include "CxxUtils/fpcompare.h"
20 
21 #include <vector>
22 #include <algorithm>
23 #include <cmath>
24 
25 using namespace std;
26 
27 
28 namespace {
29 
30 
34 struct cell_t
35 {
36  cell_t (int the_cell = 0,
37  float the_weight = 0)
38  : cell (the_cell), weight (the_weight) {}
39  bool operator< (const cell_t& other) const { return cell < other.cell; }
40  bool operator== (const cell_t& other) const { return cell == other.cell; }
41 
43  int cell;
44 
46  float weight;
47 };
48 
49 
51 const float NWEIGHT_SCALE = 1000;
52 
54 const int INDEX_BITS = 18;
55 
57 const unsigned int INDEX_MASK = (1<<INDEX_BITS) - 1;
58 
60 const int ISEQ_BITS = 32 - INDEX_BITS;
61 
63 const unsigned int ISEQ_MAX = (1 << ISEQ_BITS) - 1;
64 
65 
66 } // anonymous namespace
67 
68 
75 void
77  CaloCellLinkContainer* trans,
78  const std::string& /*key*/,
79  MsgStream& /*log*/) const
80 {
81  trans->clear();
82  trans->resize(pers->m_nClusters);
83 
84  if (pers->m_nClusters > pers->m_vISizes.size() ||
85  pers->m_nClusters > pers->m_vWSizes.size())
86  {
87  REPORT_MESSAGE_WITH_CONTEXT (MSG::WARNING,
88  "CaloCellLinkContainerCnv_p2")
89  << "Inconsistent sequence array lengths: "
90  << pers->m_nClusters
91  << " " << pers->m_vISizes.size()
92  << " " << pers->m_vWSizes.size() << "\n";
93  return;
94  }
95 
96  std::vector<unsigned>::const_iterator cIterI = pers->m_linkI.begin();
97  std::vector<unsigned>::const_iterator cIterIe = pers->m_linkI.end();
98  std::vector<float>::const_iterator cIterW = pers->m_linkW.begin();
99  std::vector<float>::const_iterator cIterWe = pers->m_linkW.end();
100 
101  for (unsigned int iCluster=0; iCluster<pers->m_nClusters; iCluster++) {
102 
103  std::vector<float>::const_iterator cIterWb = cIterW;
104 
105  bool bad = false;
106  CaloCellLink* lnk=new CaloCellLink();
107  trans->at(iCluster)=lnk; // puts an empty link == cluster
108 
109  float weight = 0;
110  unsigned int nw = 0;
111 
112  // Get number of cells and pre-allocate vector.
113  unsigned int ncells = 0;
114  for(unsigned iCell=0; iCell<pers->m_vISizes[iCluster]; iCell++) {
115  if (cIterI >= cIterIe) {
116  REPORT_MESSAGE_WITH_CONTEXT (MSG::WARNING,
117  "CaloCellLinkContainerCnv_p2")
118  << "Cell index array overrun";
119  bad = true;
120  goto endcells;
121  }
122  ncells += (*(cIterI+iCell)) >> INDEX_BITS;
123  }
124  if (ncells > 1000) ncells = 1000; // Don't go crazy.
125 
126  // Unpacking indices and making links
127  for(unsigned iCell=0; iCell<pers->m_vISizes[iCluster]; iCell++) {
128  unsigned int compositeI = *cIterI++;
129  unsigned int realI = compositeI & INDEX_MASK;
130  unsigned int lSEQ = compositeI>>INDEX_BITS;
131  for(unsigned seq=0;seq<lSEQ;seq++) {
132  if (nw == 0) {
133  if (cIterW >= cIterWe) {
134  REPORT_MESSAGE_WITH_CONTEXT (MSG::WARNING,
135  "CaloCellLinkContainerCnv_p2")
136  << "Cell weight array overrun";
137  bad = true;
138  goto endcells;
139  }
140  float compositeW = *cIterW++;
141  nw = static_cast<unsigned int> (compositeW*(1./NWEIGHT_SCALE));
142  weight = compositeW - (nw*NWEIGHT_SCALE);
143  }
144 
145  ElementLink<CaloCellContainer> el_link(pers->m_contName, realI+seq);
146  lnk->insertElement (el_link, weight, ncells);
147  --nw;
148  }
149  }//End loop over cells in cluster
150 
151  if (nw != 0) {
152  REPORT_MESSAGE_WITH_CONTEXT (MSG::WARNING,
153  "CaloCellLinkContainerCnv_p2")
154  << "Weight/index arrays out of sync";
155  bad = true;
156  }
157 
158  if (cIterW - cIterWb != static_cast<int> (pers->m_vWSizes[iCluster])) {
159  REPORT_MESSAGE_WITH_CONTEXT (MSG::WARNING,
160  "CaloCellLinkContainerCnv_p2")
161  << "Wrong number of weight sequences";
162  bad = true;
163  }
164  endcells:
165 
166  if (lnk->size() > 0 && pers->m_contName.empty()) {
167  REPORT_MESSAGE_WITH_CONTEXT (MSG::WARNING,
168  "CaloCellLinkContainerCnv_p2")
169  << "Missing cell container name";
170  bad = true;
171  }
172 
173  if (bad) {
174  REPORT_MESSAGE_WITH_CONTEXT (MSG::WARNING,
175  "CaloCellLinkContainerCnv_p2")
176  << "Discarding cluster cells";
177  lnk->removeAll();
178  }
179 
180  } // End loop over clusters
181 }
182 
183 
190 void
193  const std::string& /*key*/,
194  MsgStream& /*log*/) const
195 {
196  unsigned int nclus = trans->size();
197  pers->m_nClusters = nclus; // number of clusters in container
198 
199  // Gives it size of trans container ->
200  // will store number of indices per cluster.
201  pers->m_vISizes.resize(nclus);
202  pers->m_vWSizes.resize(nclus);
203 
204  // Guess.
205  pers->m_linkI.reserve (5*nclus);
206  pers->m_linkW.reserve (5*nclus);
207 
208  pers->m_contName.clear();
209  SG::sgkey_t first_key = 0;
210 
211 
212  const SG::ThinningCache* thinningCache = SG::getThinningCache();
213 
214  // Declare this outside the loop to save on memory allocations.
215  std::vector<cell_t> cells;
216  for (unsigned int iCluster = 0; iCluster < nclus; ++iCluster) {
217 
218  // This is actually a Navigable - contains links to the cells.
219  const CaloCellLink& lnk = *(*trans)[iCluster];
220 
221  cells.clear();
222  cells.reserve (lnk.size()+1);
223  CaloCellLink::cell_iterator it_cell = lnk.begin();
224  CaloCellLink::cell_iterator it_cell_e = lnk.end();
225 
226  for (;it_cell!=it_cell_e;it_cell++) { // Loop over cells in cluster.
228  it_cell.getInternalIterator();
229 
230  // Put stuff into array for later sorting.
231  float weight = citer->second;
232  // Trying to store weights > NWEIGHT_SCALE will lead to corrupted data.
233  if (fabs (weight) > NWEIGHT_SCALE) {
234  REPORT_MESSAGE_WITH_CONTEXT (MSG::WARNING,
235  "CaloCellLinkContainerCnv_p2")
236  << "Crazy cluster cell weight " << weight
237  << " reset to 1";
238  weight = 1;
239  }
240 
241  size_t index = citer->first.index();
242 
243  // Handle thinning.
244  if (thinningCache){
245  const SG::ThinningDecisionBase* dec =
246  thinningCache->thinning (citer->first.key());
247  if (dec) {
248  size_t persIdx = dec->index (index);
249 
250  if( SG::ThinningDecisionBase::RemovedIdx == persIdx ) {
251  // Cell index has been thinned away; skip all cells.
252  cells.clear();
253  break;
254 
255  }else
256  {
257  index = persIdx;
258  }
259  }
260  }
261 
262  SG::sgkey_t key = citer->first.key();
263  citer->first.source()->tryELRemap (key, index, key, index);
264 
265  if (pers->m_contName.empty()) {
266  pers->m_contName = *citer->first.source()->keyToString (key);
267  first_key = key;
268  }
269  else if (first_key != key) {
270  REPORT_MESSAGE_WITH_CONTEXT (MSG::ERROR,
271  "CaloCellLinkContainerCnv_p2")
272  << "Multiple container names seen; "
273  << "cannot be persistified correctly: "
274  << pers->m_contName
275  << " != " << *citer->first.source()->keyToString (key)
276  << "; cell skipped";
277  continue;
278  }
279 
280  if ((index & ~INDEX_MASK) != 0) {
281  REPORT_MESSAGE_WITH_CONTEXT (MSG::WARNING,
282  "CaloCellLinkContainerCnv_p2")
283  << "Bad cell index " << index
284  << " had high bits masked off.";
285  index &= INDEX_MASK;
286  }
287  cells.push_back (cell_t (index, weight));
288  } // End loop over cells
289 
290  // Sort links according to cell indices.
291  std::sort (cells.begin(), cells.end());
292 
293  unsigned int& nlseq = pers->m_vISizes[iCluster];
294  unsigned int& nwseq = pers->m_vWSizes[iCluster];
295  nlseq = 0;
296  nwseq = 0;
297  int last_cell = -9999;
298  float last_weight = -9999;
299  unsigned int nl_in_seq = 0;
300  unsigned int nw_in_seq = 0;
301 
302  // Sentinel.
303  cells.push_back (cell_t (-99999, -99999));
304 
305  size_t ncells = cells.size();
306  for (size_t i = 0; i < ncells; ++i) {
307  const cell_t& cell = cells[i];
308 
309  if (last_cell >= 0 && static_cast<int>(last_cell + nl_in_seq) == cell.cell &&
310  nl_in_seq < ISEQ_MAX)
311  {
312  ++nl_in_seq;
313  }
314  else {
315  if (nl_in_seq) {
316  pers->m_linkI.push_back (last_cell | (nl_in_seq << INDEX_BITS));
317  ++nlseq;
318  }
319  last_cell = cell.cell;
320  nl_in_seq = 1;
321  }
322 
323  // Don't let the weight sequence length get too large,
324  // or we'll lose precision on the weight. A limit of 255
325  // allows the weight to be represented within ~0.02.
326  //
327  // We also need to be careful comparing weights.
328  // Because we lose precision due to adding the sequence length,
329  // it's possible to have a case where last_weight != cell.weight,
330  // but they're equal after packing. That means that if we
331  // then copy the data, the persistent data won't be exactly
332  // the same. Try to take this into account in the comparison.
333  float w1 = (nw_in_seq+1)*NWEIGHT_SCALE + last_weight;
334  float w2 = (nw_in_seq+1)*NWEIGHT_SCALE + cell.weight;
336  if (equal (nextafterf (w1, w2), w2) && nw_in_seq < 255)
337  {
338  ++nw_in_seq;
339  }
340  else {
341  if (nw_in_seq) {
342  pers->m_linkW.push_back (nw_in_seq*NWEIGHT_SCALE + last_weight);
343  ++nwseq;
344  }
345  last_weight = cell.weight;
346  nw_in_seq = 1;
347  }
348  }
349  } // End loop over clusters
350 }
test_athena_ntuple_filter.seq
seq
filter configuration ## -> we use the special sequence 'AthMasterSeq' which is run before any other a...
Definition: test_athena_ntuple_filter.py:18
bad
@ bad
Definition: SUSYToolsTester.cxx:95
operator<
bool operator<(const DataVector< T > &a, const DataVector< T > &b)
Vector ordering relation.
RunTileCalibRec.cells
cells
Definition: RunTileCalibRec.py:281
operator==
bool operator==(const DataVector< T > &a, const DataVector< T > &b)
Vector equality comparison.
ReadCellNoiseFromCool.cell
cell
Definition: ReadCellNoiseFromCool.py:53
Navigable::removeAll
bool removeAll()
Definition: Navigable.h:237
index
Definition: index.py:1
SG::ThinningDecisionBase
Hold thinning decisions for one container.
Definition: ThinningDecisionBase.h:39
Navigable::insertElement
void insertElement(const CONT *objectContainer, const constituent_type *constituentObject, const RPAR &objectParameter=RPAR(), size_t sizeHint=0)
Navigable::begin
virtual object_iter begin() const
ThinningDecisionBase.h
Hold thinning decisions for one container.
Navigable::end
virtual object_iter end() const
CxxUtils::fpcompare::equal
bool equal(double a, double b)
Compare two FP numbers, working around x87 precision issues.
Definition: fpcompare.h:114
dqt_zlumi_pandas.weight
int weight
Definition: dqt_zlumi_pandas.py:190
CaloCellLinkContainerCnv_p2::persToTransWithKey
virtual void persToTransWithKey(const CaloCellLinkContainer_p2 *pers, CaloCellLinkContainer *trans, const std::string &key, MsgStream &log) const override
Convert from persistent to transient object.
Definition: CaloCellLinkContainerCnv_p2.cxx:76
lumiFormat.i
int i
Definition: lumiFormat.py:85
fpcompare.h
Workaround x86 precision issues for FP inequality comparisons.
CaloCellLinkContainerCnv_p2::transToPersWithKey
virtual void transToPersWithKey(const CaloCellLinkContainer *trans, CaloCellLinkContainer_p2 *pers, const std::string &key, MsgStream &log) const override
Convert from transient to persistent object.
Definition: CaloCellLinkContainerCnv_p2.cxx:191
getThinningCache.h
Helpers to retrieve the current thinning cache from the event context.
SG::ThinningDecisionBase::RemovedIdx
static const std::size_t RemovedIdx
Flag used to show that an index has been thinned away.
Definition: ThinningDecisionBase.h:42
SG::getThinningCache
const SG::ThinningCache * getThinningCache(const EventContext &ctx)
Retrieve the current thinning cache from the event context.
Definition: getThinningCache.cxx:28
DataVector::clear
void clear()
Erase all the elements in the collection.
REPORT_MESSAGE_WITH_CONTEXT
#define REPORT_MESSAGE_WITH_CONTEXT(LVL, CONTEXT_NAME)
Report a message, with an explicitly specified context name.
Definition: Control/AthenaKernel/AthenaKernel/errorcheck.h:345
DataVector::resize
void resize(size_type sz)
Resizes the collection to the specified number of elements.
CaloCellLinkContainerCnv_p2.h
errorcheck.h
Helpers for checking error return status codes and reporting errors.
SG::sgkey_t
uint32_t sgkey_t
Type used for hashed StoreGate key+CLID pairs.
Definition: CxxUtils/CxxUtils/sgkey_t.h:32
Navigable
Navigable template generalization to handle navigation.
Definition: Navigable.h:93
SG::ThinningDecisionBase::index
size_t index(size_t ndxOrig) const
Return the index corresponding to ndxOrig after thinning.
SG::ThinningCache::thinning
const ThinningDecisionBase * thinning(const std::string &key) const
Return thinning for key.
Definition: ThinningCache.cxx:36
InDetDD::other
@ other
Definition: InDetDD_Defs.h:16
Navigable::size
virtual unsigned int size() const
std::sort
void sort(typename std::reverse_iterator< DataModel_detail::iterator< DVL > > beg, typename std::reverse_iterator< DataModel_detail::iterator< DVL > > end, const Compare &comp)
Specialization of sort for DataVector/List.
Definition: DVL_algorithms.h:623
if
if(febId1==febId2)
Definition: LArRodBlockPhysicsV0.cxx:567
DataVector::at
const T * at(size_type n) const
Access an element, as an rvalue.
SG::ThinningCache
Cache thinning decisions for converters.
Definition: ThinningCache.h:48
DataVector::size
size_type size() const noexcept
Returns the number of elements in the collection.
mapkey::key
key
Definition: TElectronEfficiencyCorrectionTool.cxx:37
ThinningCache.h