719{
720
723 rngWrapper->
setSeed( rngName, ctx );
724 CLHEP::HepRandomEngine *rndmEngine = rngWrapper->
getEngine(ctx);
725
726 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: in SiSmearedDigizationTool::digitize() ---" );
727
728
729 const InDetDD::SiDetectorElementCollection* elementsPixel = nullptr;
731 SG::ReadCondHandle<InDetDD::SiDetectorElementCollection> pixelDetEle(
m_pixelDetEleCollKey, ctx);
732 elementsPixel = pixelDetEle.retrieve();
733 if (elementsPixel==nullptr) {
735 return StatusCode::FAILURE;
736 }
737 }
738
739 const InDetDD::SiDetectorElementCollection* elementsSCT = nullptr;
741 SG::ReadCondHandle<InDetDD::SiDetectorElementCollection> sctDetEle(
m_SCTDetEleCollKey, ctx);
742 elementsSCT = sctDetEle.retrieve();
743 if (elementsSCT==nullptr) {
745 return StatusCode::FAILURE;
746 }
747 }
748
750
753 } else {
755 }
756
758
759 while (i != e) {
761
762 const TimedHitPtr<SiHit>& hit(*i++);
763 int barrelEC = hit->getBarrelEndcap();
764 int layerDisk = hit->getLayerDisk();
766 int etaModule = hit->getEtaModule();
768
769 const InDetDD::SiDetectorElement* hitSiDetElement = nullptr;
770
772 Identifier wafer_id =
m_pixel_ID->wafer_id(barrelEC,layerDisk,phiModule,etaModule);
773 IdentifierHash wafer_hash =
m_pixel_ID->wafer_hash(wafer_id);
774 const InDetDD::SiDetectorElement* hitSiDetElement_temp = elementsPixel->
getDetectorElement(wafer_hash);
775 ATH_MSG_DEBUG(
"Pixel SiDetectorElement --> barrel_ec " << barrelEC <<
", layer_disk " << layerDisk <<
", phi_module " << phiModule <<
", eta_module " << etaModule );
776 hitSiDetElement = hitSiDetElement_temp;
777 } else {
778 side = hit->getSide();
779 Identifier idwafer =
m_sct_ID->wafer_id(barrelEC,layerDisk,phiModule,etaModule,side);
781 const InDetDD::SiDetectorElement* hitSiDetElement_temp = elementsSCT->
getDetectorElement(idhash);
782 ATH_MSG_DEBUG(
"SCT SiDetectorElement --> barrel_ec " << barrelEC <<
", layer_disk " << layerDisk <<
", phi_module " << phiModule <<
", eta_module " << etaModule <<
", side " << side);
783 hitSiDetElement = hitSiDetElement_temp;
784 }
785
786
787
788
789
790
791 if (not hitSiDetElement) {
792 ATH_MSG_FATAL(
"hitSiDetElement is null in SiSmearedDigitizationTool:"<<__LINE__);
793 throw std::runtime_error(std::string("hitSiDetElement is null in SiSmearedDigitizationTool::digitize() "));
794 }
795
796 if (
m_SmearPixel && !(hitSiDetElement->isPixel()))
continue;
797 if (!
m_SmearPixel && !(hitSiDetElement->isSCT()))
continue;
798
799 IdentifierHash waferID;
800
802 waferID =
m_pixel_ID->wafer_hash(hitSiDetElement->identify());
803 } else {
804 waferID =
m_sct_ID->wafer_hash(hitSiDetElement->identify());
805 }
806
807 HepGeom::Point3D<double> pix_localStartPosition = hit->localStartPosition();
808 HepGeom::Point3D<double> pix_localEndPosition = hit->localEndPosition();
809
810 pix_localStartPosition = hitSiDetElement->hitLocalToLocal3D(pix_localStartPosition);
811 pix_localEndPosition = hitSiDetElement->hitLocalToLocal3D(pix_localEndPosition);
812
813 double localEntryX = pix_localStartPosition.x();
814 double localEntryY = pix_localStartPosition.y();
815 double localEntryZ = pix_localStartPosition.z();
816 double localExitX = pix_localEndPosition.x();
817 double localExitY = pix_localEndPosition.y();
818 double localExitZ = pix_localEndPosition.z();
819
820 double thickness = 0.0;
821 thickness = hitSiDetElement->thickness();
822
823
825 HepGeom::Point3D<double> sct_localStartPosition = hit->localStartPosition();
826 HepGeom::Point3D<double> sct_localEndPosition = hit->localEndPosition();
827
828 sct_localStartPosition = hitSiDetElement->hitLocalToLocal3D(sct_localStartPosition);
829 sct_localEndPosition = hitSiDetElement->hitLocalToLocal3D(sct_localEndPosition);
830
831 localEntryX = sct_localStartPosition.x();
832 localEntryY = sct_localStartPosition.y();
833 localEntryZ = sct_localStartPosition.z();
834 localExitX = sct_localEndPosition.x();
835 localExitY = sct_localEndPosition.y();
836 localExitZ = sct_localEndPosition.z();
837 }
838
839 double distX = std::abs(std::abs(localExitX)-std::abs(localEntryX));
840 double distY = std::abs(std::abs(localExitY)-std::abs(localEntryY));
841
843 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: pixel start position --- " << localEntryX <<
", " << localEntryY <<
", " << localEntryZ );
844 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: pixel exit position --- " << localExitX <<
", " << localExitY <<
", " << localExitZ );
851 } else {
852 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: SCT start position --- " << localEntryX <<
", " << localEntryY <<
", " << localEntryZ );
853 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: SCT exit position --- " << localExitX <<
", " << localExitY <<
", " << localExitZ );
860 }
861
864
865
866 std::vector<Identifier> rdoList;
867
868 Amg::Vector3D localDirection(localExitX-localEntryX, localExitY-localEntryY, localExitZ-localEntryZ);
869
870 InDetDD::SiCellId entryCellId;
871 InDetDD::SiCellId exitCellId;
872
873
874 Identifier entryId = hitSiDetElement->identifierOfPosition(localEntry);
875 Identifier exitId = hitSiDetElement->identifierOfPosition(localExit);
876
877
878 entryCellId = hitSiDetElement->cellIdFromIdentifier(entryId);
879 exitCellId = hitSiDetElement->cellIdFromIdentifier(exitId);
880
881 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: entryId " << entryId <<
" --- exitId " << exitId );
882 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: entryCellId " << entryCellId <<
" --- exitCellId " << exitCellId );
883
884 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: surface " << hitSiDetElement->surface());
885
886
887 bool entryValid = entryCellId.
isValid();
888 bool exitValid = exitCellId.
isValid();
889
890 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: entryValid? " << entryValid <<
" --- exitValid? " << exitValid );
891
892 if (!entryValid && !exitValid) continue;
893
894
895 double interX = 0.5*(localEntryX+localExitX);
896 double interY = 0.5*(localEntryY+localExitY);
897
902
903
906
907 double newdistX = distX - (timesX*
m_pitch_X);
908 double newdistY = distY - (timesY*
m_pitch_Y);
909
910 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: times X --- " << timesX );
911 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: times Y --- " << timesY );
912 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: new dist X --- " << newdistX );
913 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: new dist Y --- " << newdistY );
914 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: thickness --- " << thickness );
915
916
917
918 double ProbY = 2*newdistY/(
m_pitch_Y+newdistY);
919 double ProbX = 2*newdistX/(
m_pitch_X+newdistX);
920
921 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: ProbX --- " << ProbX );
922 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: ProbY --- " << ProbY );
923
924
927
930
931 int elementX = timesX+1;
932 int elementY = timesY+1;
933
935 if (CLHEP::RandFlat::shoot(rndmEngine, 0.0, 1.0) < ProbY) {
937 elementY++;
938 } else
940
941 if (CLHEP::RandFlat::shoot(rndmEngine, 0.0, 1.0) < ProbX) {
943 elementX++;
944 } else
946 }
947
950
951
952 double temp_X = interX;
953 double temp_Y = interY;
954
956
958
959 Identifier intersectionId;
960 intersectionId = hitSiDetElement->identifierOfPosition(
intersection);
961
962 rdoList.push_back(intersectionId);
963 InDetDD::SiCellId currentCellId = hitSiDetElement->cellIdFromIdentifier(intersectionId);
964
965 if (!currentCellId.
isValid())
continue;
966
967 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: Intersection Id = " << intersectionId <<
" --- currentCellId = " << currentCellId );
968
970
973
975 ATH_MSG_WARNING(
"--- SiSmearedDigitizationTool: pitchX and/or pitchY are 0. Cluster length is forced to be 1. mm");
976
978
980
982 covariance.setIdentity();
983 covariance(Trk::locX,Trk::locX) = sigmaX*sigmaX;
984 covariance(Trk::locY,Trk::locY) = sigmaY*sigmaY;
985
986
987 pixelCluster = new InDet::PixelCluster(intersectionId,
989 std::vector<Identifier>(rdoList),
990 siWidth,
991 hitSiDetElement,
992 Amg::MatrixX(covariance));
993 m_pixelClusterMap->insert(std::pair<IdentifierHash, InDet::PixelCluster* >(waferID, pixelCluster));
994
997 return StatusCode::FAILURE;
998 }
999
1001
1002 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: pixelCluster --> " << *pixelCluster);
1003
1006
1009
1012
1017
1020
1022
1024 }
1025
1026 } else {
1027
1028
1029 InDet::SCT_Cluster * sctCluster = nullptr;
1030
1031
1032 const InDetDD::SCT_ModuleSideDesign* design_sct;
1033
1034 design_sct = dynamic_cast<const InDetDD::SCT_ModuleSideDesign*>(&hitSiDetElement->design());
1035
1036 if (!design_sct) {
1038 continue;
1039 }
1040
1041
1042 double clusterWidth = rdoList.size()*hitSiDetElement->phiPitch(
intersection);
1043 const std::pair<InDetDD::SiLocalPosition, InDetDD::SiLocalPosition> ends(design_sct->
endsOfStrip(
intersection));
1044 double stripLength = std::abs(ends.first.xEta()-ends.second.xEta());
1045
1046 InDet::SiWidth siWidth(
Amg::Vector2D(
int(rdoList.size()),1),
1048
1050
1052 mat.setIdentity();
1053 mat(Trk::locX,Trk::locX) = sigmaX*sigmaX;
1054 mat(Trk::locY,Trk::locY) = hitSiDetElement->
length()*hitSiDetElement->
length()/12.;
1055
1056
1057
1058
1059 InDetDD::
DetectorShape elShape = hitSiDetElement->design().shape();
1061 {
1062
1063 if(colRow.x() == 1) {
1064 mat(Trk::locX,Trk::locX) = pow(siWidth.phiR(),2)/12;
1065 }
1066 else if(colRow.x() == 2) {
1067 mat(Trk::locX,Trk::locX) = pow(0.27*siWidth.phiR(),2)/12;
1068 }
1069 else {
1070 mat(Trk::locX,Trk::locX) = pow(siWidth.phiR(),2)/12;
1071 }
1072
1074 double sn = hitSiDetElement->sinStereoLocal(
intersection);
1075 double sn2 = sn*sn;
1076 double cs2 = 1.-sn2;
1077 double w = hitSiDetElement->phiPitch(
intersection)/hitSiDetElement->phiPitch();
1083 }
1084
1085
1086 sctCluster = new InDet::SCT_Cluster(intersectionId,
1088 std::vector<Identifier>(rdoList),
1089 siWidth,
1090 hitSiDetElement,
1092
1093 m_sctClusterMap->insert(std::pair<IdentifierHash, InDet::SCT_Cluster* >(waferID, sctCluster));
1094
1097 return StatusCode::FAILURE;
1098 }
1099
1101
1102 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: SCT_Cluster --> " << *sctCluster);
1103
1104
1107
1108 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: BEFORE SMEARING LocalPosition --> X = " <<
m_x_SCT );
1110
1115
1117 }
1118 }
1119 }
1120 }
1121 return StatusCode::SUCCESS;
1122}
constexpr int pow(int base, int exp) noexcept
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.
virtual std::pair< SiLocalPosition, SiLocalPosition > endsOfStrip(const SiLocalPosition &position) const override=0
give the ends of strips
bool isValid() const
Test if its in a valid state.
const Amg::Vector3D & globalPosition() const
return global position reference
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
Eigen::Matrix< double, 3, 1 > Vector3D