228{
229
230 SG::WriteHandle< PRD_MultiTruthCollection > sctPrdTruth(
m_sctPrdTruthKey, ctx);
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;
235 }
236 ATH_MSG_DEBUG(
"PRD_MultiTruthCollection " << sctPrdTruth.name() <<
" registered in StoreGate" );
237
238
241 rngWrapper->
setSeed( rngName, ctx );
242 CLHEP::HepRandomEngine *rndmEngine = rngWrapper->
getEngine(ctx);
243
244
245 SG::ReadCondHandle<InDetDD::SiDetectorElementCollection> sctDetEle(
m_SCTDetEleCollKey, ctx);
246 const InDetDD::SiDetectorElementCollection* elements(sctDetEle.retrieve());
247 if (elements==nullptr) {
249 return StatusCode::FAILURE;
250 }
251
256 {
258 std::vector<int> truthIdList;
259 std::vector<Identifier> detEl;
260 while (i != e)
261 {
262 const TimedHitPtr<SiHit>& currentSiHit(*i++);
263
264
266
267 const Identifier waferId =
m_sct_ID->wafer_id(currentSiHit->getBarrelEndcap(), currentSiHit->getLayerDisk(), currentSiHit->getPhiModule(), currentSiHit->getEtaModule(), currentSiHit->getSide());
268 const IdentifierHash waferHash =
m_sct_ID->wafer_hash(waferId);
269 const InDetDD::SiDetectorElement *hitSiDetElement = elements->getDetectorElement(waferHash);
270 if (!hitSiDetElement)
271 {
273 continue;
274 }
275
276
277 const InDetDD::SCT_ModuleSideDesign* design =
dynamic_cast<const InDetDD::SCT_ModuleSideDesign*
>(&hitSiDetElement->
design());
278 if (!design)
279 {
280 ATH_MSG_WARNING (
"SCT_DetailedSurfaceChargesGenerator::process can not get "<< design) ;
281 continue;
282 }
283
284 std::vector<HepMcParticleLink> hit_vector;
285
286
287 bool isRep = false;
288 const int truthID = (currentLink.
barcode() !=0 && currentLink.
id() == 0) ? 3 : currentLink.
id();
289 const Identifier detElId = hitSiDetElement->
identify();
290 for (int j : truthIdList)
291 {
292 for (auto & k : detEl)
293 {
294 if ((truthID > 0) && (truthID == j) && (detElId == k)) {isRep = true; break;}
295 }
296 if (isRep) { break; }
297 }
298 if (isRep) { continue; }
299 truthIdList.push_back(truthID);
300 detEl.push_back(detElId);
301
303
304 double shiftX,shiftY;
308 }
309 else{
312 }
313
314 HepGeom::Point3D<double> localStartPosition = hitSiDetElement->
hitLocalToLocal3D(currentSiHit->localStartPosition());
315 HepGeom::Point3D<double> localEndPosition = hitSiDetElement->
hitLocalToLocal3D(currentSiHit->localEndPosition());
316
317 Diffuse(localStartPosition,localEndPosition, shiftX * Gaudi::Units::micrometer, shiftY * Gaudi::Units::micrometer);
318
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();
325
326
327
330 const Amg::Vector3D localExit3D(localExitX,localExitY,localExitZ);
331
332 const double thickness = hitSiDetElement->
thickness();
333
334 const Trk::Surface *hitSurface = &hitSiDetElement->
surface();
335
336 const double distX = localExitX-localEntryX;
337 const double distY = localExitY-localEntryY;
338 const double distZ = localExitZ-localEntryZ;
339
341
342 double potentialClusterPath_Geom = localDirection.mag();
343 double potentialClusterPath_Step = 0.;
344 double potentialClusterPath_Used = 0.;
345
346
347 const double slopeYX = distY/distX;
348 const double slopeXZ = distX/distZ;
349 const double slopeZX = distZ/distX;
350
351 const int signX = distX > 0. ? 1 : -1 ;
352 const int signY = distY > 0. ? 1 : -1 ;
353 const int signZ = distZ > 0. ? 1 : -1 ;
354
355
358
361
362 const bool entryValid = entryCellId.
isValid();
363 const bool exitValid = exitCellId.
isValid();
364
365
366 const double par = -localEntryZ/(localExitZ-localEntryZ);
367 const double interX = localEntryX +
par*(localExitX-localEntryX);
368 const double interY = localEntryY +
par*(localExitY-localEntryY);
369
372
373
374 const IdentifierHash detElHash = hitSiDetElement->
identifyHash();
376 const int lorentzDirection = tanLorAng > 0. ? 1 : -1;
377 const bool useLorentzDrift = std::abs(tanLorAng) > 0.01;
378
380
381 const double xLoffset = -lorentzDirection*thickness*tanLorAng;
382
383
384
385 if (!entryValid && !exitValid)
386 {
387 continue;
388 }
389
390 std::vector<Identifier> potentialClusterRDOList;
391 std::map<Identifier, double> surfaceChargeWeights;
392
393 int phiIndexMinRaw = 1000;
394 int phiIndexMaxRaw = -1000;
395
396
397 const bool singleStrip = (entryCellId == exitCellId && entryValid);
398 if (singleStrip && !useLorentzDrift)
399 {
400
401
402 surfaceChargeWeights[intersectionId] = potentialClusterPath_Geom;
403 }
404 else
405 {
406
407 int strips = 0;
408
409 Identifier currentId = entryId;
410 InDetDD::SiCellId currentCellId = entryCellId;
412 Amg::Vector3D currentPosition3D(localEntryX,localEntryY,localEntryZ);
414
415
416
417
418
420
421
422
423 if (!entryValid)
424 {
425
426 const int numberOfDiodesPhi = design->
cells();
427
429 {
430
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;
435 }
436 else
437 {
438
439 const int phiExitIndex = exitCellId.
phiIndex() <= 2 ? 0 : numberOfDiodesPhi-1;
440
441 const InDetDD::SiCellId phiExitId(phiExitIndex);
443
444
446
447 localEntryX, localEntryY,
448 localExitX, localExitY,
449 slopeYX,
450 slopeZX,
451 phiExitCenter,
452 signX);
453 }
454
455
456 currentPosition3D += currentStep3D;
457
458 potentialClusterPath_Step += currentStep3D.mag();
459 ATH_MSG_VERBOSE(
"[ cluster - sct ] Strip entry shifted by " << currentStep3D <<
" ( " << potentialClusterPath_Step <<
" )");
460
461 const Amg::Vector2D positionInFirstValid(currentPosition3D.x()+0.01*signX,currentPosition3D.y()+0.01*signY);
462
466
467 }
468
469
470
471 double lplaneStepX = 0.;
472 double lplaneStepY = 0.;
476
477 Identifier lastId = currentId;
479 ATH_MSG_VERBOSE(
"[ cluster - sct ] CurrentPosition " << currentPosition3D );
480
481
482 for ( ; ; ++strips)
483 {
484
485
486 const int currentPhiIndex = currentCellId.
phiIndex();
487
488
489
490 if (lastStrip || currentPosition3D.z() > 0.5*thickness || strips > 4)
491 {
492 break;
493 }
494
495 if (!exitValid && !currentCellId.
isValid())
496 {
497 break;
498 }
499
500 phiIndexMinRaw = currentPhiIndex < phiIndexMinRaw ? currentPhiIndex : phiIndexMinRaw;
501 phiIndexMaxRaw = currentPhiIndex > phiIndexMaxRaw ? currentPhiIndex : phiIndexMaxRaw;
502
504
505 const double currentPitchX = hitSiDetElement->
phiPitch(currentCenterPosition);
506
507 std::vector<double> lorentzLinePositions;
508 const int trackDir = slopeXZ > 0 ? 1 : -1;
509 lorentzLinePositions.reserve(2);
510
511 lorentzLinePositions.push_back(currentCenterPosition.x() + trackDir*0.5*currentPitchX);
512 lorentzLinePositions.push_back(currentCenterPosition.x() - trackDir*0.5*currentPitchX);
513
514 if (xLoffset*xLoffset > currentPitchX*currentPitchX)
515 {
516 lorentzLinePositions.push_back(currentCenterPosition.x() + lorentzDirection*1.5*currentPitchX);
517 }
518
519 bool lplaneInterInCurrent(false);
520 double lorentzPlaneHalfX(0.);
522 std::vector<double>::iterator lorentzLinePosIter = lorentzLinePositions.begin();
523
524 int lplaneDirection = 100;
525
526 for ( ; lorentzLinePosIter != lorentzLinePositions.end(); ++lorentzLinePosIter)
527 {
528
529 Trk::LineIntersection2D intersectLorentzPlane(localEntryX,-0.5*signZ*thickness,localExitX,0.5*signZ*thickness,
530 (*lorentzLinePosIter)+xLoffset,-0.5*hitDepth*thickness,
531 (*lorentzLinePosIter),0.5*hitDepth*thickness);
532
533 const double formerPlaneStepZ = intersectLorentzPlane.interY-lplaneIntersection.z();
534 if (formerPlaneStepZ*formerPlaneStepZ < 10e-5)
535 {
536 lplaneInterInCurrent = false;
537 continue;
538 }
539
540 lplaneInterInCurrent = intersectLorentzPlane.interY > -0.5*thickness
541 && intersectLorentzPlane.interY < 0.5*thickness;
542
543
544
545 if (lplaneInterInCurrent)
546 {
547
548 lorentzPlaneHalfX = (*lorentzLinePosIter)+0.5*xLoffset;
549
550 lplaneStepX = intersectLorentzPlane.interX-localEntryX;
551 lplaneStepY = lplaneStepX*slopeYX;
552
553 lplaneCandidate =
Amg::Vector3D(intersectLorentzPlane.interX,
554 localEntryY+lplaneStepY,
555 intersectLorentzPlane.interY);
556
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.;
560
561 lplaneDirection = lplaneInterInCurrent ? 0 : lplaneDirection;
562 if (lplaneInterInCurrent) {break;}
563 }
564 }
565
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() );
571
572
573 if (lplaneDirection)
574 {
575
576 const double trackZatlplaneHalfX = (lorentzPlaneHalfX-localEntryX)*slopeXZ - 0.5*thickness;
577 lplaneDirection = trackZatlplaneHalfX < 0. ? -1 : 1;
578 }
579
580
581
582 currentStep3D =
lastStrip ? localExit3D-currentPosition3D :
584
585 currentPosition3D.x(), currentPosition3D.y(),
586 localExitX, localExitY,
587 slopeYX,
588 slopeZX,
589 currentCenterPosition,
590 signX);
591
592
593 currentPosition3D += currentStep3D;
594
596 currentPosition3D.y())))
597 {
599
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;
604
606 currentPosition3D += currentStep3D;
607 }
608
609 if (std::abs(currentPosition3D.z()) > 0.501*thickness)
610 {
611
613
614 currentPosition3D -= currentStep3D;
615 currentStep3D = localExit3D-currentPosition3D;
616 currentPosition3D = localExit3D;
617
618 ATH_MSG_VERBOSE(
"[ cluster - sct ] - current position set to local Exit position !");
619 }
620
621
622 const Amg::Vector2D currentInsidePosition(currentPosition3D.x()+0.01*signX,currentPosition3D.y()+0.01*signY);
625
627
628 potentialClusterPath_Step += currentStep3D.mag();
629 ATH_MSG_VERBOSE(
"[ cluster - sct ] CurrentPosition " << currentPosition3D
630 << " ( yielded by :" << currentStep3D << ")");
631
632 driftPostPosition3D = lplaneInterInCurrent ? lplaneIntersection : currentPosition3D;
633
635 {
636
637 const int currentDrifts = (
lastStrip && lplaneInterInCurrent) ? 2 : 1;
638 for (int idrift = 0; idrift < currentDrifts; ++idrift)
639 {
640
641 const Amg::Vector3D& currentDriftPrePosition = idrift ? driftPostPosition3D : driftPrePosition3D;
642 const Amg::Vector3D& currentDriftPostPosition = idrift ? localExit3D : driftPostPosition3D;
643
644 const Amg::Vector3D driftStepCenter = 0.5*(currentDriftPrePosition+currentDriftPostPosition);
645
646 const double driftZtoSurface = 0.5*hitDepth*thickness-driftStepCenter.z();
647
648 const double driftPositionAtSurfaceX = driftStepCenter.x() + tanLorAng*driftZtoSurface;
650 driftStepCenter.y());
653 {
654
655 if (surfaceChargeWeights.find(surfaceChargeId) == surfaceChargeWeights.end())
656 {
657 surfaceChargeWeights[surfaceChargeId] = (currentDriftPostPosition-currentDriftPrePosition).
mag();
658 }
659 else
660 {
661 surfaceChargeWeights[surfaceChargeId] += (currentDriftPostPosition-currentDriftPrePosition).
mag();
662 }
663 }
664
665 }
666 }
667 else
668 {
669 surfaceChargeWeights[lastId] = currentStep3D.mag();
670 }
671
672 driftPrePosition3D = driftPostPosition3D;
673 lastId = currentId;
674
675 }
676 }
677
678
679 double totalWeight = 0.;
681 std::map<Identifier,double>::iterator weightIter = surfaceChargeWeights.begin();
682 for ( ; weightIter != surfaceChargeWeights.end(); ++weightIter)
683 {
684
685 double chargeWeight = (weightIter)->second;
686 const Identifier chargeId = (weightIter)->first;
687
688
690 {
691
695 chargeWeight *= (1.+sPar);
696 }
697
698
700
701
703 potentialClusterPath_Used += chargeWeight;
704
706 totalWeight += takenWeight;
707 potentialClusterPosition += takenWeight * chargeCenterPosition;
708
709 potentialClusterRDOList.push_back(chargeId);
710 }
711
712
714 {
715 continue;
716 }
717
718
719
720
721
722
723
724
725 potentialClusterPosition *= 1./totalWeight;
726 Identifier potentialClusterId = hitSiDetElement->
identifierOfPosition(potentialClusterPosition);
727 if (!potentialClusterId.
is_valid()) {
continue;}
728
729 const IdentifierHash waferID =
m_sct_ID->wafer_hash(hitSiDetElement->
identify());
730
731
733 unsigned int countC(0);
734 ATH_MSG_INFO(
"Before cluster merging there were " << SCT_DetElClusterMap.size() <<
" clusters in the SCT_DetElClusterMap.");
735 for(SCT_detElement_RIO_map::iterator existingClusterIter = SCT_DetElClusterMap.begin(); existingClusterIter != SCT_DetElClusterMap.end();)
736 {
737 ++countC;
738 if (countC>100)
739 {
740 ATH_MSG_WARNING(
"Over 100 clusters checked for merging - bailing out!!");
741 break;
742 }
743 if (
m_sct_ID->wafer_hash(((existingClusterIter->second)->detectorElement())->identify()) != waferID) {
continue;}
744
745 SCT_detElement_RIO_map::iterator currentExistingClusterIter = existingClusterIter++;
746
747 const InDet::SCT_Cluster *existingCluster = (currentExistingClusterIter->second);
749 if(isNeighbour)
750 {
751
752 const std::vector<Identifier> &existingClusterRDOList = existingCluster->
rdoList();
753 potentialClusterRDOList.insert(potentialClusterRDOList.end(), existingClusterRDOList.begin(), existingClusterRDOList.end() );
755 potentialClusterPosition = (potentialClusterPosition + existingClusterPosition)/2;
757
758 SCT_DetElClusterMap.erase(currentExistingClusterIter);
759
760
761
762 std::pair<PRD_MultiTruthCollection::iterator,PRD_MultiTruthCollection::iterator> saved_hit = sctPrdTruth->equal_range(existingCluster->
identify());
763 for (PRD_MultiTruthCollection::iterator this_hit = saved_hit.first; this_hit != saved_hit.second; ++this_hit)
764 {
765 hit_vector.push_back(this_hit->second);
766 }
767
768 if (saved_hit.first != saved_hit.second) sctPrdTruth->erase(existingCluster->
identify());
769
770 delete existingCluster;
772
773 }
774 }
775 ATH_MSG_INFO(
"After cluster merging there were " << SCT_DetElClusterMap.size() <<
" clusters in the SCT_DetElClusterMap.");
776 }
777
779
780 std::unique_ptr<InDet::SCT_Cluster> potentialClusterUniq = nullptr;
781
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());
785 const InDet::SiWidth siWidth(
Amg::Vector2D(
int(potentialClusterRDOList.size()),1),
Amg::Vector2D(clusterWidth,stripLength) );
786
788 {
789 ATH_MSG_INFO(
"Using ClusterMakerTool to make cluster.");
790
791
793 {
795 }
796
798 {
799
800 const double newshift = 0.5*thickness*tanLorAng;
801 const double corr = ( shift - newshift );
803 }
804 bool not_valid = false;
805 for (auto & i : potentialClusterRDOList)
806 {
808 {
809 not_valid = true;
810 break;
811 }
812 }
813 if (not_valid)
814 {
815 continue;
816 }
817 potentialClusterUniq = std::make_unique<InDet::SCT_Cluster>(
819 potentialClusterId, potentialClusterPosition,
820 std::vector<Identifier>(potentialClusterRDOList), siWidth, hitSiDetElement,
822 }
823 else
824 {
826
828 const Amg::Vector2D lcorrectedPosition =
Amg::Vector2D(potentialClusterPosition.x()+appliedShift,potentialClusterPosition.y());
830 mat.setIdentity();
831 mat(Trk::locX,Trk::locX) = (clusterWidth*clusterWidth)/12.;
832 mat(Trk::locY,Trk::locY) = (stripLength*stripLength)/12.;
833
835 {
837 const double Sn2 = Sn*Sn;
838 const double Cs2 = 1.-Sn2;
839
840
846 }
847
848
849 potentialClusterUniq = std::make_unique<InDet::SCT_Cluster>(
850 potentialClusterId, lcorrectedPosition,
851 std::vector<Identifier>(potentialClusterRDOList), siWidth, hitSiDetElement,
853 }
854
855
857 SCT_DetElClusterMap.insert(SCT_detElement_RIO_map::value_type(
858 waferID, potentialClusterUniq.release()));
859
860
861 const InDet::SCT_Cluster* potentialCluster =
it->second;
862
863
866 sctPrdTruth->insert(std::make_pair(potentialCluster->
identify(), currentLink));
867 ATH_MSG_DEBUG(
"Truth map filled with cluster" << potentialCluster <<
" and link = " << currentLink);
868 }
869 }
870 else
871 {
872 ATH_MSG_DEBUG(
"Particle link NOT valid!! Truth map NOT filled with cluster" << potentialCluster <<
" and link = " << currentLink);
873 }
874
875
876 for(const HepMcParticleLink& p: hit_vector){
877 sctPrdTruth->insert(std::make_pair(potentialCluster->
identify(), p ));
878 }
879
880 hit_vector.clear();
881
882
883 }
884
885 (void)
m_sctClusterMap->insert(SCT_DetElClusterMap.begin(), SCT_DetElClusterMap.end());
886 }
887 return StatusCode::SUCCESS;
888}
Scalar mag() const
mag method
#define ATH_MSG_VERBOSE(x)
#define AmgSymMatrix(dim)
void setSeed(const std::string &algName, const EventContext &ctx)
Set the random seed using a string (e.g.
CLHEP::HepRandomEngine * getEngine(const EventContext &ctx) const
Retrieve the random engine corresponding to the provided EventContext.
int id() const
Return the id of the target particle.
bool isValid() const
Validity check.
static HepMcParticleLink getRedirectedLink(const HepMcParticleLink &particleLink, uint32_t eventIndex, const EventContext &ctx)
Return a HepMcParticleLink pointing at the same particle, but in a different GenEvent.
int barcode() const
Return the barcode of the target particle.
bool is_valid() const
Check if id is in a valid state.
virtual DetectorShape shape() const
Shape of element.
int cells() const
number of readout stips within module side:
virtual std::pair< SiLocalPosition, SiLocalPosition > endsOfStrip(const SiLocalPosition &position) const override=0
give the ends of strips
int phiIndex() const
Get phi index. Equivalent to strip().
bool isValid() const
Test if its in a valid state.
virtual SiCellId cellIdFromIdentifier(const Identifier &identifier) const override final
SiCellId from Identifier.
virtual const SiDetectorDesign & design() const override final
access to the local description (inline):
double phiPitch() const
Pitch (inline methods)
double sinStereoLocal(const Amg::Vector2D &localPos) const
Angle of strip in local frame with respect to the etaAxis.
double length() const
Length in eta direction (z - barrel, r - endcap)
HepGeom::Point3D< double > hitLocalToLocal3D(const HepGeom::Point3D< double > &hitPosition) const
Same as previuos method but 3D.
virtual IdentifierHash identifyHash() const override final
identifier hash (inline)
double hitDepthDirection() const
Directions of hit depth,phi,eta axes relative to reconstruction local position axes (LocalPosition).
Identifier identifierOfPosition(const Amg::Vector2D &localPos) const
Full identifier of the cell for a given position: assumes a raw local position (no Lorentz shift)
Amg::Vector2D rawLocalPositionOfCell(const SiCellId &cellId) const
Returns position (center) of cell.
Trk::Surface & surface()
Element Surface.
bool nextDetectorElement(const_iterator &b, const_iterator &e)
sets an iterator range with the hits of current detector element returns a bool when done
TimedVector::const_iterator const_iterator
const Amg::Vector2D & localPosition() const
return the local position reference
Identifier identify() const
return the identifier
const std::vector< Identifier > & rdoList() const
return the List of rdo identifiers (pointers)
virtual bool insideLoc2(const Amg::Vector2D &locpo, double tol2=0.) const =0
Extend the interface to for single inside Loc 1 / Loc2 tests.
virtual const SurfaceBounds & bounds() const =0
Surface Bounds method.
std::vector< std::string > intersection(std::vector< std::string > &v1, std::vector< std::string > &v2)
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > MatrixX
Dynamic Matrix - dynamic allocation.
Eigen::Matrix< double, 2, 1 > Vector2D
Eigen::Matrix< double, 3, 1 > Vector3D
bool ignoreTruthLink(const T &p, bool vetoPileUp)
Helper function for SDO creation in PileUpTools.