ATLAS Offline Software
SiHitCollectionCnv_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 #include "InDetSimEvent/SiHit.h"
9 
10 #include <cmath>
11 
12 //CLHEP
13 #include "CLHEP/Geometry/Point3D.h"
14 // Gaudi
15 #include "GaudiKernel/MsgStream.h"
16 // Athena
18 #include "StoreGate/StoreGateSvc.h"
19 
20 // * * * stolen from eflowRec * * * //
21 inline double phicorr(double a)
22 {
23  if (a <= -M_PI)
24  {
25  return a+(2*M_PI*floor(-(a-M_PI)/(2*M_PI)));
26  }
27  else if (a > M_PI)
28  {
29  return a-(2*M_PI*floor((a+M_PI)/(2*M_PI)));
30  }
31  else
32  {
33  return a;
34  }
35 }
36 
37 // * * * stolen from eflowRec * * * //
38 inline double cycle(double a, double b)
39 {
40  double del = b-a;
41  if (del > M_PI)
42  {
43  return a+2.0*M_PI;
44  }
45  else if (del < -M_PI)
46  {
47  return a-2.0*M_PI;
48  }
49  else
50  {
51  return a;
52  }
53 }
54 
55 
56 const double SiHitCollectionCnv_p2::m_persEneUnit = 1.0e-5;
57 const double SiHitCollectionCnv_p2::m_persLenUnit = 1.0e-5;
58 const double SiHitCollectionCnv_p2::m_persAngUnit = 1.0e-5;
59 const double SiHitCollectionCnv_p2::m_2bHalfMaximum = pow(2.0, 15.0);
60 const int SiHitCollectionCnv_p2::m_2bMaximum = (unsigned short)(-1);
61 
62 
63 void SiHitCollectionCnv_p2::transToPers(const SiHitCollection* transCont, SiHitCollection_p2* persCont, MsgStream &/*log*/)
64 {
65  // Finds hits belonging to a "string" (in which the end point of one hit is the same as the start point of the next) and
66  // persistifies the end point of each hit plus the start point of the first hit in each string.
67  //
68  // Further compression is achieved by optimising the storage of the position vectors:- start (x,y,z) and (theta,phi) of
69  // first hit are stored as floats, (delta_theta,delta_phi) relative to the fisrst hit are stored as 2 byte numbers and
70  // used to specify the hit direction. All hit lengths are stored as 2 byte numbers.
71  //
72  // Additional savings are achieved by storing the energy loss for each hit as a 2 byte number and only storing the mean
73  // time of the first hit per string.
74  //
75  // See http://indico.cern.ch/getFile.py/access?contribId=11&resId=2&materialId=slides&confId=30893 for more info.
76 
77  static const double dRcut = 1.0e-7;
78  static const double dTcut = 1.0;
79 
80  int lastBarcode = -1;
81  int lastId = -1;
82  double stringFirstTheta = 0.0;
83  double stringFirstPhi = 0.0;
84  double lastT = 0.0;
85  double persSumE = 0.0;
86  double transSumE = 0.0;
87  unsigned int idx = 0;
88  unsigned int endBC = 0;
89  unsigned int endId = 0;
90  unsigned int endHit = 0;
91  HepGeom::Point3D<double> lastTransEnd(0.0, 0.0, 0.0);
92  HepGeom::Point3D<double> lastPersEnd(0.0, 0.0, 0.0);
93 
94  for (SiHitCollection::const_iterator it = transCont->begin(); it != transCont->end(); ++it) {
96  if ( siHit->particleLink().barcode() != lastBarcode || (idx-endBC)==USHRT_MAX ) {
97 
98  // store barcode once for set of consecutive hits with same barcode
99  lastBarcode = siHit->particleLink().barcode();
100  //m_barcode has type std::vector<unsigned long>, but lastBarcode could be -1
101  //make the conversion explicit here with a static_cast
102  using barcodeRep = decltype(persCont->m_barcode)::value_type;
103  persCont->m_barcode.push_back(static_cast<barcodeRep>(lastBarcode));
104 
105  if (idx > 0) {
106  persCont->m_nBC.push_back(idx - endBC);
107  endBC = idx;
108  }
109  }
110 
111  if ( ( (int)siHit->identify() != lastId ) || (idx-endId)==USHRT_MAX) {
112 
113  // store id once for set of consecutive hits with same barcode
114 
115  lastId = siHit->identify();
116  persCont->m_id.push_back(lastId);
117 
118  if (idx > 0) {
119  persCont->m_nId.push_back(idx - endId);
120  endId = idx;
121  }
122  }
123 
124  HepGeom::Point3D<double> st = siHit->localStartPosition();
125  HepGeom::Point3D<double> en = siHit->localEndPosition();
126 
127  const double dx = st.x() - lastTransEnd.x();
128  const double dy = st.y() - lastTransEnd.y();
129  const double dz = st.z() - lastTransEnd.z();
130  const double t = siHit->meanTime();
131 
132  const double dRLast = sqrt(dx * dx + dy * dy + dz * dz); // dR between end of previous hit and start of current one
133  const double dTLast = fabs(t - lastT);
134 
135  CLHEP::Hep3Vector direction(0.0, 0.0, 0.0);
136  double theta = 0.0;
137  double phi = 0.0;
138  bool startNewString = (dRLast >= dRcut || dTLast >= dTcut || (idx - endHit) == USHRT_MAX);
139 
140  if (!startNewString) {
141 
142  // hit is part of existing string
143 
144  direction = CLHEP::Hep3Vector( en.x() - lastPersEnd.x(), en.y() - lastPersEnd.y(), en.z() - lastPersEnd.z() );
145 
146  theta = direction.theta();
147  phi = phicorr( direction.phi() );
148 
149  const int dTheta_2b = (int)( (theta - stringFirstTheta) / m_persAngUnit + m_2bHalfMaximum + 0.5 );
150  const int dPhi_2b = (int)( (cycle(phi, stringFirstPhi) - stringFirstPhi) / m_persAngUnit + m_2bHalfMaximum + 0.5 );
151 
152  if ( dTheta_2b < m_2bMaximum && dTheta_2b >= 0 && dPhi_2b < m_2bMaximum && dPhi_2b >= 0) {
153  persCont->m_dTheta.push_back(dTheta_2b);
154  persCont->m_dPhi.push_back(dPhi_2b);
155  theta = stringFirstTheta + ( (double)dTheta_2b - m_2bHalfMaximum ) * m_persAngUnit;
156  phi = stringFirstPhi + ( (double)dPhi_2b - m_2bHalfMaximum ) * m_persAngUnit;
157  phi = phicorr(phi);
158  }
159  else {
160  startNewString = true;
161  }
162  }
163 
164  if (startNewString) {
165 
166  // begin new hit string
167 
168  direction = CLHEP::Hep3Vector( en.x() - st.x(), en.y() - st.y(), en.z() - st.z() );
169 
170  theta = direction.theta();
171  phi = phicorr( direction.phi() );
172 
173  persCont->m_hit1_meanTime.push_back(t);
174  persCont->m_hit1_x0.push_back(st.x());
175  persCont->m_hit1_y0.push_back(st.y());
176  persCont->m_hit1_z0.push_back(st.z());
177  persCont->m_hit1_theta.push_back(theta);
178  persCont->m_hit1_phi.push_back(phi);
179 
180  lastPersEnd = std::move(st);
181 
182  stringFirstTheta = theta;
183  stringFirstPhi = phi;
184 
185  if (idx > 0) {
186  persCont->m_nHits.push_back(idx - endHit);
187  endHit = idx;
188  }
189  }
190 
191  lastTransEnd = std::move(en);
192  transSumE += siHit->energyLoss();
193 
194  const int eneLoss_2b = (int)((transSumE - persSumE) / m_persEneUnit + 0.5); // calculated to allow recovery sum over
195  // whole hit string to chosen precision
196 
197  const int hitLength_2b = (int)(direction.mag() / m_persLenUnit + 0.5); // calculated to give the correct position to
198  // the chosen precision, NOT the length of the
199  // hit (small difference in practice).
200  double eneLoss = 0.0;
201 
202  if (eneLoss_2b >= m_2bMaximum) {
203  eneLoss = siHit->energyLoss();
204  persCont->m_hitEne_2b.push_back(m_2bMaximum);
205  persCont->m_hitEne_4b.push_back(eneLoss);
206  }
207  else {
208  eneLoss = eneLoss_2b * m_persEneUnit;
209  persCont->m_hitEne_2b.push_back(eneLoss_2b);
210  }
211 
212  double length = 0.0;
213 
214  if (hitLength_2b >= m_2bMaximum) {
215  length = direction.mag();
216  persCont->m_hitLength_2b.push_back(m_2bMaximum);
217  persCont->m_hitLength_4b.push_back(direction.mag());
218  }
219  else {
220  length = hitLength_2b * m_persLenUnit;
221  persCont->m_hitLength_2b.push_back(hitLength_2b);
222  }
223 
224  CLHEP::Hep3Vector persDir(length, 0.0, 0.0);
225  persDir.setTheta(theta);
226  persDir.setPhi(phi);
227 
228  lastPersEnd = (CLHEP::Hep3Vector)lastPersEnd + persDir;
229  persSumE += eneLoss;
230  lastT = t;
231 
232  ++idx;
233  }
234 
235  persCont->m_nBC.push_back(idx - endBC);
236  persCont->m_nId.push_back(idx - endId);
237  persCont->m_nHits.push_back(idx - endHit);
238 #ifdef ENABLE_SANITY_CHECKS
239  // Sanity check
240  const unsigned int init(0);
241  const unsigned int transContSize = transCont->size();
242  if (std::accumulate(persCont->m_nBC.begin(), persCont->m_nBC.end(), init)!=transContSize) {
243  log << MSG::ERROR << "transToPers: sum of entries of persCont->m_nBC (" << std::accumulate(persCont->m_nBC.begin(), persCont->m_nBC.end(), init) << ") does not match transient container size = " << transContSize << endmsg;
244  }
245  if (std::accumulate(persCont->m_nId.begin(), persCont->m_nId.end(), init)!=transContSize) {
246  log << MSG::ERROR << "transToPers: sum of entries of persCont->m_nId (" << std::accumulate(persCont->m_nId.begin(), persCont->m_nId.end(), init) << ") does not match transient container size = " << transContSize << endmsg;
247  }
248  if (std::accumulate(persCont->m_nHits.begin(), persCont->m_nHits.end(), init)!=transContSize) {
249  log << MSG::ERROR << "transToPers: sum of entries of persCont->m_nHits (" << std::accumulate(persCont->m_nHits.begin(), persCont->m_nHits.end(), init) << ") does not match transient container size = " << transContSize << endmsg;
250  }
251 #endif
252 }
253 
254 
256  std::unique_ptr<SiHitCollection> trans(std::make_unique<SiHitCollection>("DefaultCollectionName",persObj->m_nHits.size()));
257  persToTrans(persObj, trans.get(), log);
258  return(trans.release());
259 }
260 
261 
262 #ifdef ENABLE_SANITY_CHECKS
263 void SiHitCollectionCnv_p2::persToTrans(const SiHitCollection_p2* persCont, SiHitCollection* transCont, MsgStream &log)
264 #else
265 void SiHitCollectionCnv_p2::persToTrans(const SiHitCollection_p2* persCont, SiHitCollection* transCont, MsgStream &/*log*/)
266 #endif
267 {
268 #ifdef ENABLE_SANITY_CHECKS
269  // Sanity check
270  const unsigned int transContSize = persCont->m_hitEne_2b.size(); // this vector has one entry per transient SiHit, so its size can be used as a proxy for the transient Container size
271  const unsigned int init(0);
272  if (std::accumulate(persCont->m_nBC.begin(), persCont->m_nBC.end(), init)!=transContSize) {
273  log << MSG::ERROR << "persToTrans: sum of entries of persCont->m_nBC (" << std::accumulate(persCont->m_nBC.begin(), persCont->m_nBC.end(), init) << ") does not match transient container size = " << transContSize << endmsg;
274  }
275  if (std::accumulate(persCont->m_nId.begin(), persCont->m_nId.end(), init)!=transContSize) {
276  log << MSG::ERROR << "persToTrans: sum of entries of persCont->m_nId (" << std::accumulate(persCont->m_nId.begin(), persCont->m_nId.end(), init) << ") does not match transient container size = " << transContSize << endmsg;
277  }
278  if (std::accumulate(persCont->m_nHits.begin(), persCont->m_nHits.end(), init)!=transContSize) {
279  log << MSG::ERROR << "persToTrans: sum of entries of persCont->m_nHits (" << std::accumulate(persCont->m_nHits.begin(), persCont->m_nHits.end(), init) << ") does not match transient container size = " << transContSize << endmsg;
280  }
281 #endif
282 
283  unsigned int hitCount = 0;
284  unsigned int angleCount = 0;
285  unsigned int idxBC = 0;
286  unsigned int idxId = 0;
287  unsigned int idxEne4b = 0;
288  unsigned int idxLen4b = 0;
289  unsigned int endHit = 0;
290  unsigned int endBC = 0;
291  unsigned int endId = 0;
292  // Assume that all Hits should be linked to the hard-scatter GenEvent
294  const int event_number = HepMcParticleLink::getEventNumberAtPosition (0, sg);
295  for (unsigned int i = 0; i < persCont->m_nHits.size(); i++) {
296 
297  if (persCont->m_nHits[i]) {
298 
299  const unsigned int start = endHit;
300  endHit += persCont->m_nHits[i];
301 
302  const double t0 = persCont->m_hit1_meanTime[i];
303  const double theta0 = persCont->m_hit1_theta[i];
304  const double phi0 = persCont->m_hit1_phi[i];
305  HepGeom::Point3D<double> endLast(persCont->m_hit1_x0[i], persCont->m_hit1_y0[i], persCont->m_hit1_z0[i]);
306 
307  for (unsigned int j = start; j < endHit; j++) {
308 
309  if (j >= endBC + persCont->m_nBC[idxBC])
310  endBC += persCont->m_nBC[idxBC++];
311 
312  if (j >= endId + persCont->m_nId[idxId])
313  endId += persCont->m_nId[idxId++];
314 
315  const double eneLoss_2b = persCont->m_hitEne_2b[hitCount];
316  const double hitLength_2b = persCont->m_hitLength_2b[hitCount];
317 
318  const double eneLoss = (eneLoss_2b < m_2bMaximum) ? eneLoss_2b * m_persEneUnit : persCont->m_hitEne_4b[idxEne4b++];
319  const double length = (hitLength_2b < m_2bMaximum) ? hitLength_2b * m_persLenUnit : persCont->m_hitLength_4b[idxLen4b++];
320 
321  const double dTheta = (j > start) ? ((double)persCont->m_dTheta[angleCount] - m_2bHalfMaximum) * m_persAngUnit : 0.0;
322  const double dPhi = (j > start) ? ((double)persCont->m_dPhi[angleCount] - m_2bHalfMaximum) * m_persAngUnit : 0.0;
323 
324  const double meanTime = t0;
325  const double theta = theta0 + dTheta;
326  const double phi = phicorr(phi0 + dPhi);
327 
328  CLHEP::Hep3Vector r(length, 0.0, 0.0);
329  r.setTheta(theta);
330  r.setPhi(phi);
331 
332  HepGeom::Point3D<double> endThis( endLast + r );
333 
334  HepMcParticleLink partLink( persCont->m_barcode[idxBC], event_number, HepMcParticleLink::IS_EVENTNUM, HepMcParticleLink::IS_BARCODE, sg);
335  if ( HepMC::BarcodeBased::is_truth_suppressed_pileup(static_cast<int>(persCont->m_barcode[idxBC])) ) {
337  }
338  transCont->Emplace( endLast, endThis, eneLoss, meanTime, partLink, persCont->m_id[idxId]);
339 
340  endLast = std::move(endThis);
341 
342  ++hitCount;
343  if (j > start) ++angleCount;
344  }
345  }
346  }
347 }
beamspotman.r
def r
Definition: beamspotman.py:674
SiHitCollection_p2::m_hit1_phi
std::vector< float > m_hit1_phi
Definition: SiHitCollection_p2.h:29
phi
Scalar phi() const
phi method
Definition: AmgMatrixBasePlugin.h:67
SiHitCollectionCnv_p2::m_2bMaximum
static const int m_2bMaximum
Definition: SiHitCollectionCnv_p2.h:37
SiHitCollection_p2::m_hitLength_4b
std::vector< float > m_hitLength_4b
Definition: SiHitCollection_p2.h:40
SiHit.h
xAOD::short
short
Definition: Vertex_v1.cxx:165
SiHitCollection_p2::m_hitLength_2b
std::vector< unsigned short > m_hitLength_2b
Definition: SiHitCollection_p2.h:33
SiHitCollection_p2
Definition: SiHitCollection_p2.h:18
SiHitCollection_p2::m_hit1_theta
std::vector< float > m_hit1_theta
Definition: SiHitCollection_p2.h:28
InDetAccessor::phi0
@ phi0
Definition: InDetAccessor.h:33
mergePhysValFiles.start
start
Definition: DataQuality/DataQualityUtils/scripts/mergePhysValFiles.py:13
cycle
double cycle(double a, double b)
Definition: SiHitCollectionCnv_p2.cxx:38
theta
Scalar theta() const
theta method
Definition: AmgMatrixBasePlugin.h:75
ALFA_EventTPCnv_Dict::t0
std::vector< ALFA_RawData_p1 > t0
Definition: ALFA_EventTPCnvDict.h:42
AtlasHitsVector< SiHit >
skel.it
it
Definition: skel.GENtoEVGEN.py:407
M_PI
#define M_PI
Definition: ActiveFraction.h:11
SiHitCollection_p2::m_hit1_y0
std::vector< float > m_hit1_y0
Definition: SiHitCollection_p2.h:26
SG::CurrentEventStore::store
static IProxyDict * store()
Fetch the current store.
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
AtlasHitsVector::begin
const_iterator begin() const
Definition: AtlasHitsVector.h:130
IProxyDict
A proxy dictionary.
Definition: AthenaKernel/AthenaKernel/IProxyDict.h:47
AtlasHitsVector< SiHit >::const_iterator
CONT::const_iterator const_iterator
Definition: AtlasHitsVector.h:42
SiHitCollection_p2::m_nBC
std::vector< unsigned short > m_nBC
Definition: SiHitCollection_p2.h:43
SiHitCollection_p2::m_id
std::vector< unsigned long > m_id
Definition: SiHitCollection_p2.h:45
AtlasHitsVector::Emplace
void Emplace(Args &&... args)
Definition: AtlasHitsVector.h:80
SiHitCollectionCnv_p2::m_persAngUnit
static const double m_persAngUnit
Definition: SiHitCollectionCnv_p2.h:35
SiHitCollection_p2::m_hitEne_4b
std::vector< float > m_hitEne_4b
Definition: SiHitCollection_p2.h:38
lumiFormat.i
int i
Definition: lumiFormat.py:85
endmsg
#define endmsg
Definition: AnalysisConfig_Ntuple.cxx:63
SiHitCollection_p2::m_nHits
std::vector< unsigned short > m_nHits
Definition: SiHitCollection_p2.h:30
TauGNNUtils::Variables::Track::dPhi
bool dPhi(const xAOD::TauJet &tau, const xAOD::TauTrack &track, double &out)
Definition: TauGNNUtils.cxx:549
SiHitCollection_p2::m_hit1_x0
std::vector< float > m_hit1_x0
Definition: SiHitCollection_p2.h:25
SiHitCollectionCnv_p2::persToTrans
virtual void persToTrans(const SiHitCollection_p2 *persCont, SiHitCollection *transCont, MsgStream &log)
Definition: SiHitCollectionCnv_p2.cxx:265
SiHitCollection_p2::m_barcode
std::vector< unsigned long > m_barcode
Definition: SiHitCollection_p2.h:42
xAOD::double
double
Definition: CompositeParticle_v1.cxx:159
SiHitCollectionCnv_p2.h
SiHitCollection_p2::m_hit1_meanTime
std::vector< float > m_hit1_meanTime
Definition: SiHitCollection_p2.h:24
MagicNumbers.h
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:76
PlotCalibFromCool.en
en
Definition: PlotCalibFromCool.py:399
HepMC::BarcodeBased::is_truth_suppressed_pileup
bool is_truth_suppressed_pileup(const T &p)
Method to establish if a particle (or barcode) corresponds to truth-suppressed pile-up.
Definition: MagicNumbers.h:196
SiHitCollectionCnv_p2::m_persEneUnit
static const double m_persEneUnit
Definition: SiHitCollectionCnv_p2.h:33
python.PyKernel.init
def init(v_theApp, v_rootStream=None)
Definition: PyKernel.py:45
SiHitCollection_p2::m_dTheta
std::vector< unsigned short > m_dTheta
Definition: SiHitCollection_p2.h:35
makeTRTBarrelCans.dy
tuple dy
Definition: makeTRTBarrelCans.py:21
SiHitCollection_p2.h
SiHitCollectionCnv_p2::transToPers
virtual void transToPers(const SiHitCollection *transCont, SiHitCollection_p2 *persCont, MsgStream &log)
Definition: SiHitCollectionCnv_p2.cxx:63
runIDAlign.accumulate
accumulate
Update flags based on parser line args.
Definition: runIDAlign.py:60
phicorr
double phicorr(double a)
Definition: SiHitCollectionCnv_p2.cxx:21
a
TList * a
Definition: liststreamerinfos.cxx:10
python.CaloAddPedShiftConfig.int
int
Definition: CaloAddPedShiftConfig.py:45
AtlasHitsVector::end
const_iterator end() const
Definition: AtlasHitsVector.h:133
makeTRTBarrelCans.dx
tuple dx
Definition: makeTRTBarrelCans.py:20
python.CaloCondTools.log
log
Definition: CaloCondTools.py:20
LArNewCalib_DelayDump_OFC_Cali.idx
idx
Definition: LArNewCalib_DelayDump_OFC_Cali.py:69
AtlasHitsVector::size
size_type size() const
Definition: AtlasHitsVector.h:142
EBC_PU_SUPPRESSED
@ EBC_PU_SUPPRESSED
Definition: MagicNumbers.h:30
SiHitCollection_p2::m_hit1_z0
std::vector< float > m_hit1_z0
Definition: SiHitCollection_p2.h:27
SiHitCollectionCnv_p2::createTransient
virtual SiHitCollection * createTransient(const SiHitCollection_p2 *persObj, MsgStream &log)
Definition: SiHitCollectionCnv_p2.cxx:255
SiHitCollection_p2::m_dPhi
std::vector< unsigned short > m_dPhi
Definition: SiHitCollection_p2.h:36
SiHitCollectionCnv_p2::m_persLenUnit
static const double m_persLenUnit
Definition: SiHitCollectionCnv_p2.h:34
value_type
Definition: EDM_MasterSearch.h:11
SiHitCollection_p2::m_nId
std::vector< unsigned short > m_nId
Definition: SiHitCollection_p2.h:46
pow
constexpr int pow(int base, int exp) noexcept
Definition: ap_fixedTest.cxx:15
StoreGateSvc.h
SiHitCollectionCnv_p2::m_2bHalfMaximum
static const double m_2bHalfMaximum
Definition: SiHitCollectionCnv_p2.h:36
length
double length(const pvec &v)
Definition: FPGATrackSimLLPDoubletHoughTransformTool.cxx:26
SiHitCollection_p2::m_hitEne_2b
std::vector< unsigned short > m_hitEne_2b
Definition: SiHitCollection_p2.h:32
SiHitCollection.h