13 #include "CLHEP/Geometry/Point3D.h"
14 #include "CLHEP/Units/SystemOfUnits.h"
17 #include "GaudiKernel/MsgStream.h"
18 #include "GaudiKernel/ThreadLocalContext.h"
50 static const double dRcut = 1.0e-7*
CLHEP::mm;
51 static const double dTcut = 1.0*
CLHEP::ns;
55 const EventContext& ctx = Gaudi::Hive::currentContext();
62 unsigned int endTruthID = 0;
63 unsigned int endId = 0;
64 unsigned int endHit = 0;
65 HepGeom::Point3D<double> lastEnd(0.0, 0.0, 0.0);
71 const int truthId = currentLink->
id();
78 if ( lastTruthId != truthId || lastIndex !=
index || (
idx - endTruthID > 65500) ) {
79 lastTruthId = truthId;
81 const unsigned short persIndex =
static_cast<unsigned short>(
index);
83 log << MSG::WARNING <<
"Attempting to persistify an eventIndex larger than max unsigned short!" <<
endmsg;
87 persCont->
m_truthID.push_back(
static_cast<unsigned int>(lastTruthId));
96 if ( (
int)trtHit->GetParticleEncoding() != lastId ||
idx - endId > 65500) {
98 lastId = trtHit->GetParticleEncoding();
99 persCont->
m_id.push_back(lastId);
101 persCont->
m_nId.push_back(
idx - endId);
106 const HepGeom::Point3D<double> hitStart(trtHit->GetPreStepX(), trtHit->GetPreStepY(), trtHit->GetPreStepZ());
108 const double meanTime = trtHit->GetGlobalTime();
109 const double dTLast = std::abs(meanTime - lastT);
110 const double dRLast = lastEnd.distance(hitStart);
115 if ( dRLast >= dRcut || dTLast >= dTcut ) {
127 const unsigned int strawId = trtHit->GetHitID();
128 persCont->
m_strawId1b.push_back( (
unsigned char)(strawId % 256) );
129 persCont->
m_strawId2b.push_back( (
unsigned short)(strawId / 256) );
130 if ( strawId>16777215 )
131 log << MSG::WARNING <<
"TRT_HitCollectionCnv: strawId > 2^24-1 cannot be persistified correctly! " <<
endmsg;
140 const double startR = sqrt( hitStart.x()*hitStart.x() + hitStart.y()*hitStart.y() );
141 unsigned short istartRflag;
154 const double startPhi = atan2( hitStart.y(), hitStart.x() );
155 persCont->
m_startPhi.push_back( (
unsigned char)( (startPhi+
M_PI)/(2.0*
M_PI)*256.0 ) );
171 istartZ = (istartZ << 1) | istartRflag;
172 persCont->
m_startZ.push_back( istartZ );
194 const HepGeom::Point3D<double> hitEnd(trtHit->GetPostStepX(), trtHit->GetPostStepY(), trtHit->GetPostStepZ());
195 const HepGeom::Point3D<double> hitLength = (hitEnd - hitStart);
214 double kinEne = trtHit->GetKineticEnergy() * 1.0e9;
215 double steplength = hitLength.distance() * 1.0e9;
216 if ( kinEne < 1.0 ) kinEne=1.0;
217 if ( steplength < 1.0 ) steplength=1.0;
218 if ( kinEne > 9.0e18 ) kinEne=9.0e18;
219 if ( steplength > 9.0e18 ) steplength=9.0e18;
220 const unsigned int kexponent = (
unsigned int)ceil(log10(kinEne)/0.30102999566398);
221 const unsigned int sexponent = (
unsigned int)ceil(log10(steplength)/0.30102999566398);
222 const unsigned int kmantissa = (
unsigned int)(kinEne/
pow(2.0,kexponent)*1024) - 512;
223 const unsigned int smantissa = (
unsigned int)(steplength/
pow(2.0,sexponent)*1024) - 512;
224 persCont->
m_kinEne.push_back( (kmantissa << 6) | kexponent );
225 persCont->
m_steplength.push_back( (smantissa << 6) | sexponent );
235 const double endR = sqrt( hitEnd.x()*hitEnd.x() + hitEnd.y()*hitEnd.y() );
236 unsigned short iendRflag;
242 persCont->
m_endR.push_back( (
unsigned char)(endR/(2.0*
CLHEP::mm)*256.0) );
250 const double endPhi = atan2( hitEnd.y(), hitEnd.x() );
251 persCont->
m_endPhi.push_back( (
unsigned char)( (endPhi+
M_PI)/(2.0*
M_PI)*256.0 ) );
259 unsigned short idZsign = (hitLength.z()>0.0) ? 1 : 0;
260 unsigned short imeanTime = ( meanTime < 75.0*
CLHEP::ns ) ? (
unsigned short)(meanTime/(75.0*
CLHEP::ns)*1024.0) : 1023;
261 if ( imeanTime == 1023 ) persCont->
m_meanTimeof.push_back( (
float)meanTime );
262 imeanTime = (imeanTime << 2) | (idZsign << 1) | iendRflag;
270 (
int)(abs(lastId)/100000) == 41 ||
271 (
int)(abs(lastId)/10000000) == 1
272 ) persCont->
m_hitEne.push_back( (
float)(trtHit->GetEnergyDeposit()) );
293 persCont->
m_nId.push_back(
idx - endId);
301 std::unique_ptr<TRTUncompressedHitCollection> trans(std::make_unique<TRTUncompressedHitCollection>(
"DefaultCollectionName",persObj->
m_nHits.size()));
303 return(trans.release());
310 const EventContext& ctx = Gaudi::Hive::currentContext();
315 unsigned int meanTimeofCount=0, startRCount=0, endRCount=0, hitEneCount=0;
316 unsigned int idxTruthID=0, idxId=0, endHit=0, endTruthID=0, endId=0;
322 for (
unsigned int i = 0;
i < persCont->
m_nHits.size();
i++ ) {
326 const unsigned int startHit = endHit;
334 const unsigned int strawId = i2*256+i1;
339 const unsigned int istartPhi = persCont->
m_startPhi[
i];
340 const double startPhi = -
M_PI + (istartPhi+0.5)*2.0*
M_PI/256.0;
345 const unsigned int istartZ = persCont->
m_startZ[
i] >> 1;
351 const unsigned int istartRflag = persCont->
m_startZ[
i] & 1;
357 if ( istartRflag == 1 ) {
361 const unsigned int istartR = persCont->
m_startR[startRCount++];
362 startR = (istartR+0.5)*2.0*
CLHEP::mm/256.0;
369 double startX = startR*
cos(startPhi);
370 double startY = startR*
sin(startPhi);
386 for (
unsigned int j = startHit; j < endHit; j++ ) {
388 if ( j >= endTruthID + persCont->
m_nTruthID[idxTruthID] ) endTruthID += persCont->
m_nTruthID[idxTruthID++];
389 if ( j >= endId + persCont->
m_nId[idxId] ) endId += persCont->
m_nId[idxId++];
394 const unsigned int imeanTime = persCont->
m_meanTime[j] >> 2;
395 double meanTime = (imeanTime+0.5)*75.0*
CLHEP::ns/1024.0;
396 if ( imeanTime == 1023 ) meanTime = (
double)persCont->
m_meanTimeof[meanTimeofCount++];
401 const unsigned int idZsign = (persCont->
m_meanTime[j] >> 1 ) & 1;
406 const unsigned int iendRflag = persCont->
m_meanTime[j] & 1;
411 const double hitEne = ( persCont->
m_id[idxId] == 22 ||
412 (
int)(abs(persCont->
m_id[idxId])/100000) == 41 ||
413 (
int)(abs(persCont->
m_id[idxId])/10000000) == 1
419 const unsigned int iendPhi = persCont->
m_endPhi[j];
420 double endPhi = -
M_PI + (iendPhi+0.5)*2.0*
M_PI/256.0;
426 if ( iendRflag==1 ) {
430 const unsigned int iendR = persCont->
m_endR[endRCount++];
438 double endX = endR*
cos(endPhi);
439 double endY = endR*
sin(endPhi);
451 const int kmantissa = persCont->
m_kinEne[j] >> 6;
453 const int kexponent = persCont->
m_kinEne[j] & 0x3F;
455 const double kinEne = (kmantissa+512.5)/1024 *
pow(2.0,kexponent) / 1.0e9;
456 double g4steplength = (smantissa+512.5)/1024 *
pow(2.0,sexponent) / 1.0e9;
457 if ( idZsign==0 ) g4steplength = -g4steplength;
462 double dX = endX-startX;
463 double dY = endY-startY;
465 double dXY2 = dX*dX+dY*dY;
466 double dL2 = g4steplength*g4steplength;
469 if (g4steplength<0.0) dZ=-dZ;
472 dX = dX * sqrt(dL2/dXY2);
473 dY = dY * sqrt(dL2/dXY2);
480 double endZ = startZ + dZ;
507 transCont->
Emplace( strawId, partLink, persCont->
m_id[idxId],
508 kinEne, hitEne, startX, startY, startZ,
509 endX, endY, endZ, meanTime );
515 startX = endXo; startY = endYo; startZ = endZ;