13 #include "CLHEP/Geometry/Point3D.h"
15 #include "GaudiKernel/MsgStream.h"
77 static const double dRcut = 1.0e-7;
78 static const double dTcut = 1.0;
82 double stringFirstTheta = 0.0;
83 double stringFirstPhi = 0.0;
85 double persSumE = 0.0;
86 double transSumE = 0.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);
96 if ( siHit->particleLink().barcode() != lastBarcode || (
idx-endBC)==USHRT_MAX ) {
99 lastBarcode = siHit->particleLink().barcode();
103 persCont->
m_barcode.push_back(
static_cast<barcodeRep
>(lastBarcode));
106 persCont->
m_nBC.push_back(
idx - endBC);
111 if ( ( (
int)siHit->identify() != lastId ) || (
idx-endId)==USHRT_MAX) {
115 lastId = siHit->identify();
116 persCont->
m_id.push_back(lastId);
119 persCont->
m_nId.push_back(
idx - endId);
124 HepGeom::Point3D<double> st = siHit->localStartPosition();
125 HepGeom::Point3D<double>
en = siHit->localEndPosition();
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();
132 const double dRLast = sqrt(
dx *
dx +
dy *
dy + dz * dz);
133 const double dTLast = fabs(
t - lastT);
135 CLHEP::Hep3Vector direction(0.0, 0.0, 0.0);
138 bool startNewString = (dRLast >= dRcut || dTLast >= dTcut || (
idx - endHit) == USHRT_MAX);
140 if (!startNewString) {
144 direction = CLHEP::Hep3Vector(
en.x() - lastPersEnd.x(),
en.y() - lastPersEnd.y(),
en.z() - lastPersEnd.z() );
146 theta = direction.theta();
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);
160 startNewString =
true;
164 if (startNewString) {
168 direction = CLHEP::Hep3Vector(
en.x() - st.x(),
en.y() - st.y(),
en.z() - st.z() );
170 theta = direction.theta();
180 lastPersEnd = std::move(st);
182 stringFirstTheta =
theta;
183 stringFirstPhi =
phi;
191 lastTransEnd = std::move(
en);
192 transSumE += siHit->energyLoss();
200 double eneLoss = 0.0;
203 eneLoss = siHit->energyLoss();
224 CLHEP::Hep3Vector persDir(
length, 0.0, 0.0);
225 persDir.setTheta(
theta);
228 lastPersEnd = (CLHEP::Hep3Vector)lastPersEnd + persDir;
235 persCont->
m_nBC.push_back(
idx - endBC);
236 persCont->
m_nId.push_back(
idx - endId);
238 #ifdef ENABLE_SANITY_CHECKS
240 const unsigned int init(0);
241 const unsigned int transContSize = transCont->
size();
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;
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;
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;
256 std::unique_ptr<SiHitCollection> trans(std::make_unique<SiHitCollection>(
"DefaultCollectionName",persObj->
m_nHits.size()));
258 return(trans.release());
262 #ifdef ENABLE_SANITY_CHECKS
268 #ifdef ENABLE_SANITY_CHECKS
270 const unsigned int transContSize = persCont->
m_hitEne_2b.size();
271 const unsigned int init(0);
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;
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;
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;
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;
295 for (
unsigned int i = 0;
i < persCont->
m_nHits.size();
i++) {
299 const unsigned int start = endHit;
307 for (
unsigned int j =
start; j < endHit; j++) {
309 if (j >= endBC + persCont->
m_nBC[idxBC])
310 endBC += persCont->
m_nBC[idxBC++];
312 if (j >= endId + persCont->
m_nId[idxId])
313 endId += persCont->
m_nId[idxId++];
315 const double eneLoss_2b = persCont->
m_hitEne_2b[hitCount];
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++];
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;
324 const double meanTime =
t0;
325 const double theta = theta0 + dTheta;
328 CLHEP::Hep3Vector
r(
length, 0.0, 0.0);
332 HepGeom::Point3D<double> endThis( endLast +
r );
338 transCont->
Emplace( endLast, endThis, eneLoss, meanTime, partLink, persCont->
m_id[idxId]);
340 endLast = std::move(endThis);
343 if (j >
start) ++angleCount;