231 sctPrdTruth = std::make_unique< PRD_MultiTruthCollection >();
232 if ( !sctPrdTruth.isValid() ) {
233 ATH_MSG_FATAL(
"Could not record collection " << sctPrdTruth.name() );
234 return StatusCode::FAILURE;
236 ATH_MSG_DEBUG(
"PRD_MultiTruthCollection " << sctPrdTruth.name() <<
" registered in StoreGate" );
241 rngWrapper->
setSeed( rngName, ctx );
242 CLHEP::HepRandomEngine *rndmEngine = rngWrapper->
getEngine(ctx);
247 if (elements==
nullptr) {
249 return StatusCode::FAILURE;
258 std::vector<int> truthIdList;
259 std::vector<Identifier> detEl;
267 const Identifier waferId =
m_sct_ID->
wafer_id(currentSiHit->getBarrelEndcap(), currentSiHit->getLayerDisk(), currentSiHit->getPhiModule(), currentSiHit->getEtaModule(), currentSiHit->getSide());
270 if (!hitSiDetElement)
280 ATH_MSG_WARNING (
"SCT_DetailedSurfaceChargesGenerator::process can not get "<< design) ;
284 std::vector<HepMcParticleLink> hit_vector;
288 const int truthID = (currentLink.
barcode() !=0 && currentLink.
id() == 0) ? 3 : currentLink.
id();
290 for (
int j : truthIdList)
292 for (
auto &
k : detEl)
294 if ((truthID > 0) && (truthID == j) && (detElId ==
k)) {isRep =
true;
break;}
296 if (isRep) {
break; }
298 if (isRep) {
continue; }
299 truthIdList.push_back(truthID);
300 detEl.push_back(detElId);
304 double shiftX,shiftY;
314 HepGeom::Point3D<double> localStartPosition = hitSiDetElement->
hitLocalToLocal3D(currentSiHit->localStartPosition());
315 HepGeom::Point3D<double> localEndPosition = hitSiDetElement->
hitLocalToLocal3D(currentSiHit->localEndPosition());
319 const double localEntryX = localStartPosition.x();
320 const double localEntryY = localStartPosition.y();
321 const double localEntryZ = localStartPosition.z();
322 const double localExitX = localEndPosition.x();
323 const double localExitY = localEndPosition.y();
324 const double localExitZ = localEndPosition.z();
330 const Amg::Vector3D localExit3D(localExitX,localExitY,localExitZ);
332 const double thickness = hitSiDetElement->
thickness();
336 const double distX = localExitX-localEntryX;
337 const double distY = localExitY-localEntryY;
338 const double distZ = localExitZ-localEntryZ;
342 double potentialClusterPath_Geom = localDirection.mag();
343 double potentialClusterPath_Step = 0.;
344 double potentialClusterPath_Used = 0.;
347 const double slopeYX = distY/distX;
348 const double slopeXZ = distX/distZ;
349 const double slopeZX = distZ/distX;
351 const int signX = distX > 0. ? 1 : -1 ;
352 const int signY = distY > 0. ? 1 : -1 ;
353 const int signZ = distZ > 0. ? 1 : -1 ;
362 const bool entryValid = entryCellId.
isValid();
363 const bool exitValid = exitCellId.
isValid();
366 const double par = -localEntryZ/(localExitZ-localEntryZ);
367 const double interX = localEntryX +
par*(localExitX-localEntryX);
368 const double interY = localEntryY +
par*(localExitY-localEntryY);
376 const int lorentzDirection = tanLorAng > 0. ? 1 : -1;
377 const bool useLorentzDrift = std::abs(tanLorAng) > 0.01;
381 const double xLoffset = -lorentzDirection*thickness*tanLorAng;
385 if (!entryValid && !exitValid)
390 std::vector<Identifier> potentialClusterRDOList;
391 std::map<Identifier, double> surfaceChargeWeights;
393 int phiIndexMinRaw = 1000;
394 int phiIndexMaxRaw = -1000;
397 const bool singleStrip = (entryCellId == exitCellId && entryValid);
398 if (singleStrip && !useLorentzDrift)
402 surfaceChargeWeights[intersectionId] = potentialClusterPath_Geom;
412 Amg::Vector3D currentPosition3D(localEntryX,localEntryY,localEntryZ);
426 const int numberOfDiodesPhi = design->
cells();
431 const double stepInY = signY*(std::abs(localEntryY)-0.5*hitSiDetElement->
length());
432 const double stepInX = stepInY*distX/distY;
433 const double stepInZ = stepInY*distZ/distY;
439 const int phiExitIndex = exitCellId.
phiIndex() <= 2 ? 0 : numberOfDiodesPhi-1;
447 localEntryX, localEntryY,
448 localExitX, localExitY,
456 currentPosition3D += currentStep3D;
458 potentialClusterPath_Step += currentStep3D.mag();
459 ATH_MSG_VERBOSE(
"[ cluster - sct ] Strip entry shifted by " << currentStep3D <<
" ( " << potentialClusterPath_Step <<
" )");
461 const Amg::Vector2D positionInFirstValid(currentPosition3D.x()+0.01*signX,currentPosition3D.y()+0.01*signY);
471 double lplaneStepX = 0.;
472 double lplaneStepY = 0.;
479 ATH_MSG_VERBOSE(
"[ cluster - sct ] CurrentPosition " << currentPosition3D );
486 const int currentPhiIndex = currentCellId.
phiIndex();
490 if (
lastStrip || currentPosition3D.z() > 0.5*thickness || strips > 4)
495 if (!exitValid && !currentCellId.
isValid())
500 phiIndexMinRaw = currentPhiIndex < phiIndexMinRaw ? currentPhiIndex : phiIndexMinRaw;
501 phiIndexMaxRaw = currentPhiIndex > phiIndexMaxRaw ? currentPhiIndex : phiIndexMaxRaw;
505 const double currentPitchX = hitSiDetElement->
phiPitch(currentCenterPosition);
507 std::vector<double> lorentzLinePositions;
508 const int trackDir = slopeXZ > 0 ? 1 : -1;
509 lorentzLinePositions.reserve(2);
511 lorentzLinePositions.push_back(currentCenterPosition.x() + trackDir*0.5*currentPitchX);
512 lorentzLinePositions.push_back(currentCenterPosition.x() - trackDir*0.5*currentPitchX);
514 if (xLoffset*xLoffset > currentPitchX*currentPitchX)
516 lorentzLinePositions.push_back(currentCenterPosition.x() + lorentzDirection*1.5*currentPitchX);
519 bool lplaneInterInCurrent(
false);
520 double lorentzPlaneHalfX(0.);
524 int lplaneDirection = 100;
526 for ( ; lorentzLinePosIter != lorentzLinePositions.end(); ++lorentzLinePosIter)
530 (*lorentzLinePosIter)+xLoffset,-0.5*hitDepth*thickness,
531 (*lorentzLinePosIter),0.5*hitDepth*thickness);
533 const double formerPlaneStepZ = intersectLorentzPlane.interY-lplaneIntersection.z();
534 if (formerPlaneStepZ*formerPlaneStepZ < 10
e-5)
536 lplaneInterInCurrent =
false;
540 lplaneInterInCurrent = intersectLorentzPlane.interY > -0.5*thickness
541 && intersectLorentzPlane.interY < 0.5*thickness;
545 if (lplaneInterInCurrent)
548 lorentzPlaneHalfX = (*lorentzLinePosIter)+0.5*xLoffset;
550 lplaneStepX = intersectLorentzPlane.interX-localEntryX;
551 lplaneStepY = lplaneStepX*slopeYX;
553 lplaneCandidate =
Amg::Vector3D(intersectLorentzPlane.interX,
554 localEntryY+lplaneStepY,
555 intersectLorentzPlane.interY);
557 const double distToNextLineY = 0.5*hitSiDetElement->
length()-lplaneCandidate.y();
558 const double distToPrevLineY = -0.5*hitSiDetElement->
length()-lplaneCandidate.y();
559 lplaneInterInCurrent = (distToNextLineY*distToPrevLineY) < 0.;
561 lplaneDirection = lplaneInterInCurrent ? 0 : lplaneDirection;
562 if (lplaneInterInCurrent) {
break;}
566 if (lplaneInterInCurrent) {lplaneIntersection = lplaneCandidate;}
567 ATH_MSG_VERBOSE(
"[ cluster - pix ] Lorentz plane intersection x/y/z = "
568 << lplaneCandidate.x() <<
", "
569 << lplaneCandidate.y() <<
", "
570 << lplaneCandidate.z() );
576 const double trackZatlplaneHalfX = (lorentzPlaneHalfX-localEntryX)*slopeXZ - 0.5*thickness;
577 lplaneDirection = trackZatlplaneHalfX < 0. ? -1 : 1;
582 currentStep3D =
lastStrip ? localExit3D-currentPosition3D :
585 currentPosition3D.x(), currentPosition3D.y(),
586 localExitX, localExitY,
589 currentCenterPosition,
593 currentPosition3D += currentStep3D;
596 currentPosition3D.y())))
600 currentPosition3D -= currentStep3D;
601 const double stepInY = signY*0.5*hitSiDetElement->
length()-currentPosition3D.y();
602 const double stepInX = distX/distY*stepInY;
603 const double stepInZ = slopeZX*stepInX;
606 currentPosition3D += currentStep3D;
609 if (std::abs(currentPosition3D.z()) > 0.501*thickness)
614 currentPosition3D -= currentStep3D;
615 currentStep3D = localExit3D-currentPosition3D;
616 currentPosition3D = localExit3D;
618 ATH_MSG_VERBOSE(
"[ cluster - sct ] - current position set to local Exit position !");
622 const Amg::Vector2D currentInsidePosition(currentPosition3D.x()+0.01*signX,currentPosition3D.y()+0.01*signY);
628 potentialClusterPath_Step += currentStep3D.mag();
629 ATH_MSG_VERBOSE(
"[ cluster - sct ] CurrentPosition " << currentPosition3D
630 <<
" ( yielded by :" << currentStep3D <<
")");
632 driftPostPosition3D = lplaneInterInCurrent ? lplaneIntersection : currentPosition3D;
637 const int currentDrifts = (
lastStrip && lplaneInterInCurrent) ? 2 : 1;
638 for (
int idrift = 0; idrift < currentDrifts; ++idrift)
641 const Amg::Vector3D& currentDriftPrePosition = idrift ? driftPostPosition3D : driftPrePosition3D;
642 const Amg::Vector3D& currentDriftPostPosition = idrift ? localExit3D : driftPostPosition3D;
644 const Amg::Vector3D driftStepCenter = 0.5*(currentDriftPrePosition+currentDriftPostPosition);
646 const double driftZtoSurface = 0.5*hitDepth*thickness-driftStepCenter.z();
648 const double driftPositionAtSurfaceX = driftStepCenter.x() + tanLorAng*driftZtoSurface;
650 driftStepCenter.y());
655 if (surfaceChargeWeights.find(surfaceChargeId) == surfaceChargeWeights.end())
657 surfaceChargeWeights[surfaceChargeId] = (currentDriftPostPosition-currentDriftPrePosition).
mag();
661 surfaceChargeWeights[surfaceChargeId] += (currentDriftPostPosition-currentDriftPrePosition).
mag();
669 surfaceChargeWeights[lastId] = currentStep3D.mag();
672 driftPrePosition3D = driftPostPosition3D;
679 double totalWeight = 0.;
682 for ( ; weightIter != surfaceChargeWeights.end(); ++weightIter)
685 double chargeWeight = (weightIter)->
second;
695 chargeWeight *= (1.+sPar);
703 potentialClusterPath_Used += chargeWeight;
706 totalWeight += takenWeight;
707 potentialClusterPosition += takenWeight * chargeCenterPosition;
709 potentialClusterRDOList.push_back(chargeId);
725 potentialClusterPosition *= 1./totalWeight;
727 if (!potentialClusterId.
is_valid()) {
continue;}
733 unsigned int countC(0);
734 ATH_MSG_INFO(
"Before cluster merging there were " << SCT_DetElClusterMap.size() <<
" clusters in the SCT_DetElClusterMap.");
740 ATH_MSG_WARNING(
"Over 100 clusters checked for merging - bailing out!!");
743 if (
m_sct_ID->
wafer_hash(((existingClusterIter->second)->detectorElement())->identify()) != waferID) {
continue;}
752 const std::vector<Identifier> &existingClusterRDOList = existingCluster->
rdoList();
753 potentialClusterRDOList.insert(potentialClusterRDOList.end(), existingClusterRDOList.begin(), existingClusterRDOList.end() );
755 potentialClusterPosition = (potentialClusterPosition + existingClusterPosition)/2;
758 SCT_DetElClusterMap.erase(currentExistingClusterIter);
762 std::pair<PRD_MultiTruthCollection::iterator,PRD_MultiTruthCollection::iterator> saved_hit = sctPrdTruth->equal_range(existingCluster->
identify());
765 hit_vector.push_back(this_hit->second);
768 if (saved_hit.first != saved_hit.second) sctPrdTruth->erase(existingCluster->
identify());
770 delete existingCluster;
775 ATH_MSG_INFO(
"After cluster merging there were " << SCT_DetElClusterMap.size() <<
" clusters in the SCT_DetElClusterMap.");
780 std::unique_ptr<InDet::SCT_Cluster> potentialClusterUniq =
nullptr;
782 const double clusterWidth(potentialClusterRDOList.size()*hitSiDetElement->
phiPitch(potentialClusterPosition));
783 const std::pair<InDetDD::SiLocalPosition, InDetDD::SiLocalPosition> ends(design->
endsOfStrip(potentialClusterPosition));
784 const double stripLength = std::abs(ends.first.xEta()-ends.second.xEta());
789 ATH_MSG_INFO(
"Using ClusterMakerTool to make cluster.");
800 const double newshift = 0.5*thickness*tanLorAng;
801 const double corr = ( shift - newshift );
804 bool not_valid =
false;
805 for (
auto &
i : potentialClusterRDOList)
817 potentialClusterUniq = std::make_unique<InDet::SCT_Cluster>(
819 potentialClusterId, potentialClusterPosition,
820 std::vector<Identifier>(potentialClusterRDOList), siWidth, hitSiDetElement,
828 const Amg::Vector2D lcorrectedPosition =
Amg::Vector2D(potentialClusterPosition.x()+appliedShift,potentialClusterPosition.y());
837 const double Sn2 = Sn*Sn;
838 const double Cs2 = 1.-Sn2;
849 potentialClusterUniq = std::make_unique<InDet::SCT_Cluster>(
850 potentialClusterId, lcorrectedPosition,
851 std::vector<Identifier>(potentialClusterRDOList), siWidth, hitSiDetElement,
857 SCT_DetElClusterMap.insert(SCT_detElement_RIO_map::value_type(
858 waferID, potentialClusterUniq.release()));
866 sctPrdTruth->insert(std::make_pair(potentialCluster->
identify(), currentLink));
867 ATH_MSG_DEBUG(
"Truth map filled with cluster" << potentialCluster <<
" and link = " << currentLink);
872 ATH_MSG_DEBUG(
"Particle link NOT valid!! Truth map NOT filled with cluster" << potentialCluster <<
" and link = " << currentLink);
877 sctPrdTruth->insert(std::make_pair(potentialCluster->
identify(),
p ));
885 (void)
m_sctClusterMap->insert(SCT_DetElClusterMap.begin(), SCT_DetElClusterMap.end());
887 return StatusCode::SUCCESS;