744{
745
748 rngWrapper->
setSeed( rngName, ctx );
749 CLHEP::HepRandomEngine *rndmEngine = rngWrapper->
getEngine(ctx);
750
751 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: in SiSmearedDigizationTool::digitize() ---" );
752
753
754 const InDetDD::SiDetectorElementCollection* elementsPixel = nullptr;
756 SG::ReadCondHandle<InDetDD::SiDetectorElementCollection> pixelDetEle(
m_pixelDetEleCollKey, ctx);
757 elementsPixel = pixelDetEle.retrieve();
758 if (elementsPixel==nullptr) {
760 return StatusCode::FAILURE;
761 }
762 }
763
764 const InDetDD::SiDetectorElementCollection* elementsSCT = nullptr;
766 SG::ReadCondHandle<InDetDD::SiDetectorElementCollection> sctDetEle(
m_SCTDetEleCollKey, ctx);
767 elementsSCT = sctDetEle.retrieve();
768 if (elementsSCT==nullptr) {
770 return StatusCode::FAILURE;
771 }
772 }
773
775
778 } else {
780 }
781
783
784 while (i != e) {
786
787 const TimedHitPtr<SiHit>&
hit(*i++);
788 int barrelEC =
hit->getBarrelEndcap();
789 int layerDisk =
hit->getLayerDisk();
791 int etaModule =
hit->getEtaModule();
793
794 const InDetDD::SiDetectorElement* hitSiDetElement = nullptr;
795
797 Identifier wafer_id =
m_pixel_ID->wafer_id(barrelEC,layerDisk,phiModule,etaModule);
798 IdentifierHash wafer_hash =
m_pixel_ID->wafer_hash(wafer_id);
799 const InDetDD::SiDetectorElement* hitSiDetElement_temp = elementsPixel->
getDetectorElement(wafer_hash);
800 ATH_MSG_DEBUG(
"Pixel SiDetectorElement --> barrel_ec " << barrelEC <<
", layer_disk " << layerDisk <<
", phi_module " << phiModule <<
", eta_module " << etaModule );
801 hitSiDetElement = hitSiDetElement_temp;
802 } else {
804 Identifier idwafer =
m_sct_ID->wafer_id(barrelEC,layerDisk,phiModule,etaModule,side);
806 const InDetDD::SiDetectorElement* hitSiDetElement_temp = elementsSCT->
getDetectorElement(idhash);
807 ATH_MSG_DEBUG(
"SCT SiDetectorElement --> barrel_ec " << barrelEC <<
", layer_disk " << layerDisk <<
", phi_module " << phiModule <<
", eta_module " << etaModule <<
", side " << side);
808 hitSiDetElement = hitSiDetElement_temp;
809 }
810
811
812
813
814
815
816 if (not hitSiDetElement) {
817 ATH_MSG_FATAL(
"hitSiDetElement is null in SiSmearedDigitizationTool:"<<__LINE__);
818 throw std::runtime_error(std::string("hitSiDetElement is null in SiSmearedDigitizationTool::digitize() "));
819 }
820
821 if (
m_SmearPixel && !(hitSiDetElement->isPixel()))
continue;
822 if (!
m_SmearPixel && !(hitSiDetElement->isSCT()))
continue;
823
824 IdentifierHash waferID;
825
827 waferID =
m_pixel_ID->wafer_hash(hitSiDetElement->identify());
828 } else {
829 waferID =
m_sct_ID->wafer_hash(hitSiDetElement->identify());
830 }
831
832 HepGeom::Point3D<double> pix_localStartPosition =
hit->localStartPosition();
833 HepGeom::Point3D<double> pix_localEndPosition =
hit->localEndPosition();
834
835 pix_localStartPosition = hitSiDetElement->hitLocalToLocal3D(pix_localStartPosition);
836 pix_localEndPosition = hitSiDetElement->hitLocalToLocal3D(pix_localEndPosition);
837
838 double localEntryX = pix_localStartPosition.x();
839 double localEntryY = pix_localStartPosition.y();
840 double localEntryZ = pix_localStartPosition.z();
841 double localExitX = pix_localEndPosition.x();
842 double localExitY = pix_localEndPosition.y();
843 double localExitZ = pix_localEndPosition.z();
844
845 double thickness = 0.0;
846 thickness = hitSiDetElement->thickness();
847
848
850 HepGeom::Point3D<double> sct_localStartPosition =
hit->localStartPosition();
851 HepGeom::Point3D<double> sct_localEndPosition =
hit->localEndPosition();
852
853 sct_localStartPosition = hitSiDetElement->hitLocalToLocal3D(sct_localStartPosition);
854 sct_localEndPosition = hitSiDetElement->hitLocalToLocal3D(sct_localEndPosition);
855
856 localEntryX = sct_localStartPosition.x();
857 localEntryY = sct_localStartPosition.y();
858 localEntryZ = sct_localStartPosition.z();
859 localExitX = sct_localEndPosition.x();
860 localExitY = sct_localEndPosition.y();
861 localExitZ = sct_localEndPosition.z();
862 }
863
864 double distX = std::abs(std::abs(localExitX)-std::abs(localEntryX));
865 double distY = std::abs(std::abs(localExitY)-std::abs(localEntryY));
866
868 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: pixel start position --- " << localEntryX <<
", " << localEntryY <<
", " << localEntryZ );
869 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: pixel exit position --- " << localExitX <<
", " << localExitY <<
", " << localExitZ );
876 } else {
877 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: SCT start position --- " << localEntryX <<
", " << localEntryY <<
", " << localEntryZ );
878 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: SCT exit position --- " << localExitX <<
", " << localExitY <<
", " << localExitZ );
885 }
886
889
890
891 std::vector<Identifier> rdoList;
892
893 Amg::Vector3D localDirection(localExitX-localEntryX, localExitY-localEntryY, localExitZ-localEntryZ);
894
895 InDetDD::SiCellId entryCellId;
896 InDetDD::SiCellId exitCellId;
897
898
899 Identifier entryId = hitSiDetElement->identifierOfPosition(localEntry);
900 Identifier exitId = hitSiDetElement->identifierOfPosition(localExit);
901
902
903 entryCellId = hitSiDetElement->cellIdFromIdentifier(entryId);
904 exitCellId = hitSiDetElement->cellIdFromIdentifier(exitId);
905
906 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: entryId " << entryId <<
" --- exitId " << exitId );
907 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: entryCellId " << entryCellId <<
" --- exitCellId " << exitCellId );
908
909 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: surface " << hitSiDetElement->surface());
910
911
912 bool entryValid = entryCellId.
isValid();
913 bool exitValid = exitCellId.
isValid();
914
915 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: entryValid? " << entryValid <<
" --- exitValid? " << exitValid );
916
917 if (!entryValid && !exitValid) continue;
918
919
920 double interX = 0.5*(localEntryX+localExitX);
921 double interY = 0.5*(localEntryY+localExitY);
922
927
928
931
932 double newdistX = distX - (timesX*
m_pitch_X);
933 double newdistY = distY - (timesY*
m_pitch_Y);
934
935 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: times X --- " << timesX );
936 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: times Y --- " << timesY );
937 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: new dist X --- " << newdistX );
938 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: new dist Y --- " << newdistY );
939 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: thickness --- " << thickness );
940
941
942
943 double ProbY = 2*newdistY/(
m_pitch_Y+newdistY);
944 double ProbX = 2*newdistX/(
m_pitch_X+newdistX);
945
946 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: ProbX --- " << ProbX );
947 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: ProbY --- " << ProbY );
948
949
952
955
956 int elementX = timesX+1;
957 int elementY = timesY+1;
958
960 if (CLHEP::RandFlat::shoot(rndmEngine, 0.0, 1.0) < ProbY) {
962 elementY++;
963 } else
965
966 if (CLHEP::RandFlat::shoot(rndmEngine, 0.0, 1.0) < ProbX) {
968 elementX++;
969 } else
971 }
972
975
976
977 double temp_X = interX;
978 double temp_Y = interY;
979
981
983
984 Identifier intersectionId;
985 intersectionId = hitSiDetElement->identifierOfPosition(
intersection);
986
987 rdoList.push_back(intersectionId);
988 InDetDD::SiCellId currentCellId = hitSiDetElement->cellIdFromIdentifier(intersectionId);
989
990 if (!currentCellId.
isValid())
continue;
991
992 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: Intersection Id = " << intersectionId <<
" --- currentCellId = " << currentCellId );
993
995
998
1000 ATH_MSG_WARNING(
"--- SiSmearedDigitizationTool: pitchX and/or pitchY are 0. Cluster length is forced to be 1. mm");
1001
1003
1005
1007 covariance.setIdentity();
1008 covariance(Trk::locX,Trk::locX) = sigmaX*sigmaX;
1009 covariance(Trk::locY,Trk::locY) = sigmaY*sigmaY;
1010
1011
1012 pixelCluster = new InDet::PixelCluster(intersectionId,
1014 std::move(rdoList),
1015 siWidth,
1016 hitSiDetElement,
1017 Amg::MatrixX(covariance));
1018 m_pixelClusterMap->insert(std::pair<IdentifierHash, InDet::PixelCluster* >(waferID, pixelCluster));
1019
1022 return StatusCode::FAILURE;
1023 }
1024
1026
1027 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: pixelCluster --> " << *pixelCluster);
1028
1031
1034
1037
1042
1045
1047
1049 }
1050
1051 } else {
1052
1053
1054 InDet::SCT_Cluster * sctCluster = nullptr;
1055
1056
1057 const InDetDD::SCT_ModuleSideDesign* design_sct;
1058
1059 design_sct = dynamic_cast<const InDetDD::SCT_ModuleSideDesign*>(&hitSiDetElement->design());
1060
1061 if (!design_sct) {
1063 continue;
1064 }
1065
1066
1067 double clusterWidth = rdoList.size()*hitSiDetElement->phiPitch(
intersection);
1068 const std::pair<InDetDD::SiLocalPosition, InDetDD::SiLocalPosition> ends(design_sct->
endsOfStrip(
intersection));
1069 double stripLength = std::abs(ends.first.xEta()-ends.second.xEta());
1070
1071 InDet::SiWidth siWidth(
Amg::Vector2D(
int(rdoList.size()),1),
1073
1075
1077 mat.setIdentity();
1078 mat(Trk::locX,Trk::locX) = sigmaX*sigmaX;
1079 mat(Trk::locY,Trk::locY) = hitSiDetElement->
length()*hitSiDetElement->
length()/12.;
1080
1081
1082
1083
1084 InDetDD::
DetectorShape elShape = hitSiDetElement->design().shape();
1086 {
1087
1088 if(colRow.x() == 1) {
1089 mat(Trk::locX,Trk::locX) = pow(siWidth.phiR(),2)/12;
1090 }
1091 else if(colRow.x() == 2) {
1092 mat(Trk::locX,Trk::locX) = pow(0.27*siWidth.phiR(),2)/12;
1093 }
1094 else {
1095 mat(Trk::locX,Trk::locX) = pow(siWidth.phiR(),2)/12;
1096 }
1097
1099 double sn = hitSiDetElement->sinStereoLocal(
intersection);
1100 double sn2 = sn*sn;
1101 double cs2 = 1.-sn2;
1102 double w = hitSiDetElement->phiPitch(
intersection)/hitSiDetElement->phiPitch();
1108 }
1109
1110
1111 sctCluster = new InDet::SCT_Cluster(intersectionId,
1113 std::move(rdoList),
1114 siWidth,
1115 hitSiDetElement,
1117
1118 m_sctClusterMap->insert(std::pair<IdentifierHash, InDet::SCT_Cluster* >(waferID, sctCluster));
1119
1122 return StatusCode::FAILURE;
1123 }
1124
1126
1127 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: SCT_Cluster --> " << *sctCluster);
1128
1129
1132
1133 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: BEFORE SMEARING LocalPosition --> X = " <<
m_x_SCT );
1135
1140
1142 }
1143 }
1144 }
1145 }
1146 return StatusCode::SUCCESS;
1147}
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