716{
717
720 rngWrapper->
setSeed( rngName, ctx );
721 CLHEP::HepRandomEngine *rndmEngine = rngWrapper->
getEngine(ctx);
722
723 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: in SiSmearedDigizationTool::digitize() ---" );
724
725
726 const InDetDD::SiDetectorElementCollection* elementsPixel = nullptr;
728 SG::ReadCondHandle<InDetDD::SiDetectorElementCollection> pixelDetEle(
m_pixelDetEleCollKey, ctx);
729 elementsPixel = pixelDetEle.retrieve();
730 if (elementsPixel==nullptr) {
732 return StatusCode::FAILURE;
733 }
734 }
735
736 const InDetDD::SiDetectorElementCollection* elementsSCT = nullptr;
738 SG::ReadCondHandle<InDetDD::SiDetectorElementCollection> sctDetEle(
m_SCTDetEleCollKey, ctx);
739 elementsSCT = sctDetEle.retrieve();
740 if (elementsSCT==nullptr) {
742 return StatusCode::FAILURE;
743 }
744 }
745
747
750 } else {
752 }
753
755
756 while (i != e) {
758
759 const TimedHitPtr<SiHit>& hit(*i++);
760 int barrelEC = hit->getBarrelEndcap();
761 int layerDisk = hit->getLayerDisk();
763 int etaModule = hit->getEtaModule();
765
766 const InDetDD::SiDetectorElement* hitSiDetElement = nullptr;
767
769 Identifier wafer_id =
m_pixel_ID->wafer_id(barrelEC,layerDisk,phiModule,etaModule);
770 IdentifierHash wafer_hash =
m_pixel_ID->wafer_hash(wafer_id);
771 const InDetDD::SiDetectorElement* hitSiDetElement_temp = elementsPixel->
getDetectorElement(wafer_hash);
772 ATH_MSG_DEBUG(
"Pixel SiDetectorElement --> barrel_ec " << barrelEC <<
", layer_disk " << layerDisk <<
", phi_module " << phiModule <<
", eta_module " << etaModule );
773 hitSiDetElement = hitSiDetElement_temp;
774 } else {
775 side = hit->getSide();
776 Identifier idwafer =
m_sct_ID->wafer_id(barrelEC,layerDisk,phiModule,etaModule,side);
778 const InDetDD::SiDetectorElement* hitSiDetElement_temp = elementsSCT->
getDetectorElement(idhash);
779 ATH_MSG_DEBUG(
"SCT SiDetectorElement --> barrel_ec " << barrelEC <<
", layer_disk " << layerDisk <<
", phi_module " << phiModule <<
", eta_module " << etaModule <<
", side " << side);
780 hitSiDetElement = hitSiDetElement_temp;
781 }
782
783
784
785
786
787
788 if (not hitSiDetElement) {
789 ATH_MSG_FATAL(
"hitSiDetElement is null in SiSmearedDigitizationTool:"<<__LINE__);
790 throw std::runtime_error(std::string("hitSiDetElement is null in SiSmearedDigitizationTool::digitize() "));
791 }
792
793 if (
m_SmearPixel && !(hitSiDetElement->isPixel()))
continue;
794 if (!
m_SmearPixel && !(hitSiDetElement->isSCT()))
continue;
795
796 IdentifierHash waferID;
797
799 waferID =
m_pixel_ID->wafer_hash(hitSiDetElement->identify());
800 } else {
801 waferID =
m_sct_ID->wafer_hash(hitSiDetElement->identify());
802 }
803
804 HepGeom::Point3D<double> pix_localStartPosition = hit->localStartPosition();
805 HepGeom::Point3D<double> pix_localEndPosition = hit->localEndPosition();
806
807 pix_localStartPosition = hitSiDetElement->hitLocalToLocal3D(pix_localStartPosition);
808 pix_localEndPosition = hitSiDetElement->hitLocalToLocal3D(pix_localEndPosition);
809
810 double localEntryX = pix_localStartPosition.x();
811 double localEntryY = pix_localStartPosition.y();
812 double localEntryZ = pix_localStartPosition.z();
813 double localExitX = pix_localEndPosition.x();
814 double localExitY = pix_localEndPosition.y();
815 double localExitZ = pix_localEndPosition.z();
816
817 double thickness = 0.0;
818 thickness = hitSiDetElement->thickness();
819
820
822 HepGeom::Point3D<double> sct_localStartPosition = hit->localStartPosition();
823 HepGeom::Point3D<double> sct_localEndPosition = hit->localEndPosition();
824
825 sct_localStartPosition = hitSiDetElement->hitLocalToLocal3D(sct_localStartPosition);
826 sct_localEndPosition = hitSiDetElement->hitLocalToLocal3D(sct_localEndPosition);
827
828 localEntryX = sct_localStartPosition.x();
829 localEntryY = sct_localStartPosition.y();
830 localEntryZ = sct_localStartPosition.z();
831 localExitX = sct_localEndPosition.x();
832 localExitY = sct_localEndPosition.y();
833 localExitZ = sct_localEndPosition.z();
834 }
835
836 double distX = std::abs(std::abs(localExitX)-std::abs(localEntryX));
837 double distY = std::abs(std::abs(localExitY)-std::abs(localEntryY));
838
840 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: pixel start position --- " << localEntryX <<
", " << localEntryY <<
", " << localEntryZ );
841 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: pixel exit position --- " << localExitX <<
", " << localExitY <<
", " << localExitZ );
848 } else {
849 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: SCT start position --- " << localEntryX <<
", " << localEntryY <<
", " << localEntryZ );
850 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: SCT exit position --- " << localExitX <<
", " << localExitY <<
", " << localExitZ );
857 }
858
861
862
863 std::vector<Identifier> rdoList;
864
865 Amg::Vector3D localDirection(localExitX-localEntryX, localExitY-localEntryY, localExitZ-localEntryZ);
866
867 InDetDD::SiCellId entryCellId;
868 InDetDD::SiCellId exitCellId;
869
870
871 Identifier entryId = hitSiDetElement->identifierOfPosition(localEntry);
872 Identifier exitId = hitSiDetElement->identifierOfPosition(localExit);
873
874
875 entryCellId = hitSiDetElement->cellIdFromIdentifier(entryId);
876 exitCellId = hitSiDetElement->cellIdFromIdentifier(exitId);
877
878 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: entryId " << entryId <<
" --- exitId " << exitId );
879 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: entryCellId " << entryCellId <<
" --- exitCellId " << exitCellId );
880
881 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: surface " << hitSiDetElement->surface());
882
883
884 bool entryValid = entryCellId.
isValid();
885 bool exitValid = exitCellId.
isValid();
886
887 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: entryValid? " << entryValid <<
" --- exitValid? " << exitValid );
888
889 if (!entryValid && !exitValid) continue;
890
891
892 double interX = 0.5*(localEntryX+localExitX);
893 double interY = 0.5*(localEntryY+localExitY);
894
899
900
903
904 double newdistX = distX - (timesX*
m_pitch_X);
905 double newdistY = distY - (timesY*
m_pitch_Y);
906
907 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: times X --- " << timesX );
908 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: times Y --- " << timesY );
909 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: new dist X --- " << newdistX );
910 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: new dist Y --- " << newdistY );
911 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: thickness --- " << thickness );
912
913
914
915 double ProbY = 2*newdistY/(
m_pitch_Y+newdistY);
916 double ProbX = 2*newdistX/(
m_pitch_X+newdistX);
917
918 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: ProbX --- " << ProbX );
919 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: ProbY --- " << ProbY );
920
921
924
927
928 int elementX = timesX+1;
929 int elementY = timesY+1;
930
932 if (CLHEP::RandFlat::shoot(rndmEngine, 0.0, 1.0) < ProbY) {
934 elementY++;
935 } else
937
938 if (CLHEP::RandFlat::shoot(rndmEngine, 0.0, 1.0) < ProbX) {
940 elementX++;
941 } else
943 }
944
947
948
949 double temp_X = interX;
950 double temp_Y = interY;
951
953
955
956 Identifier intersectionId;
957 intersectionId = hitSiDetElement->identifierOfPosition(
intersection);
958
959 rdoList.push_back(intersectionId);
960 InDetDD::SiCellId currentCellId = hitSiDetElement->cellIdFromIdentifier(intersectionId);
961
962 if (!currentCellId.
isValid())
continue;
963
964 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: Intersection Id = " << intersectionId <<
" --- currentCellId = " << currentCellId );
965
967
970
972 ATH_MSG_WARNING(
"--- SiSmearedDigitizationTool: pitchX and/or pitchY are 0. Cluster length is forced to be 1. mm");
973
975
977
979 covariance.setIdentity();
980 covariance(Trk::locX,Trk::locX) = sigmaX*sigmaX;
981 covariance(Trk::locY,Trk::locY) = sigmaY*sigmaY;
982
983
984 pixelCluster = new InDet::PixelCluster(intersectionId,
986 std::vector<Identifier>(rdoList),
987 siWidth,
988 hitSiDetElement,
989 Amg::MatrixX(covariance));
990 m_pixelClusterMap->insert(std::pair<IdentifierHash, InDet::PixelCluster* >(waferID, pixelCluster));
991
994 return StatusCode::FAILURE;
995 }
996
998
999 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: pixelCluster --> " << *pixelCluster);
1000
1003
1006
1009
1014
1017
1019
1021 }
1022
1023 } else {
1024
1025
1026 InDet::SCT_Cluster * sctCluster = nullptr;
1027
1028
1029 const InDetDD::SCT_ModuleSideDesign* design_sct;
1030
1031 design_sct = dynamic_cast<const InDetDD::SCT_ModuleSideDesign*>(&hitSiDetElement->design());
1032
1033 if (!design_sct) {
1035 continue;
1036 }
1037
1038
1039 double clusterWidth = rdoList.size()*hitSiDetElement->phiPitch(
intersection);
1040 const std::pair<InDetDD::SiLocalPosition, InDetDD::SiLocalPosition> ends(design_sct->
endsOfStrip(
intersection));
1041 double stripLength = std::abs(ends.first.xEta()-ends.second.xEta());
1042
1043 InDet::SiWidth siWidth(
Amg::Vector2D(
int(rdoList.size()),1),
1045
1047
1049 mat.setIdentity();
1050 mat(Trk::locX,Trk::locX) = sigmaX*sigmaX;
1051 mat(Trk::locY,Trk::locY) = hitSiDetElement->
length()*hitSiDetElement->
length()/12.;
1052
1053
1054
1055
1056 InDetDD::
DetectorShape elShape = hitSiDetElement->design().shape();
1058 {
1059
1060 if(colRow.x() == 1) {
1061 mat(Trk::locX,Trk::locX) = pow(siWidth.phiR(),2)/12;
1062 }
1063 else if(colRow.x() == 2) {
1064 mat(Trk::locX,Trk::locX) = pow(0.27*siWidth.phiR(),2)/12;
1065 }
1066 else {
1067 mat(Trk::locX,Trk::locX) = pow(siWidth.phiR(),2)/12;
1068 }
1069
1071 double sn = hitSiDetElement->sinStereoLocal(
intersection);
1072 double sn2 = sn*sn;
1073 double cs2 = 1.-sn2;
1074 double w = hitSiDetElement->phiPitch(
intersection)/hitSiDetElement->phiPitch();
1080 }
1081
1082
1083 sctCluster = new InDet::SCT_Cluster(intersectionId,
1085 std::vector<Identifier>(rdoList),
1086 siWidth,
1087 hitSiDetElement,
1089
1090 m_sctClusterMap->insert(std::pair<IdentifierHash, InDet::SCT_Cluster* >(waferID, sctCluster));
1091
1094 return StatusCode::FAILURE;
1095 }
1096
1098
1099 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: SCT_Cluster --> " << *sctCluster);
1100
1101
1104
1105 ATH_MSG_DEBUG(
"--- SiSmearedDigitizationTool: BEFORE SMEARING LocalPosition --> X = " <<
m_x_SCT );
1107
1112
1114 }
1115 }
1116 }
1117 }
1118 return StatusCode::SUCCESS;
1119}
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