18 #include "Identifier/Identifier.h"
33 #include "CLHEP/Geometry/Point3D.h"
37 #define AUXDATA(OBJ, TYP, NAME) \
38 static const SG::AuxElement::Accessor<TYP> acc_##NAME (#NAME); acc_##NAME(*(OBJ))
41 unsigned int makeKey(
short phi,
char eta,
char layer) {
53 m_PixelHelper(nullptr),
54 m_useSiHitsGeometryMatching(true),
55 m_firstEventWarnings(true),
60 declareProperty(
"UseTruthInfo", m_useTruthInfo=
false);
61 declareProperty(
"UseSiHitsGeometryMatching", m_useSiHitsGeometryMatching=
true);
62 declareProperty(
"WriteSDOs", m_writeSDOs =
false);
63 declareProperty(
"WriteSiHits", m_writeSiHits =
false);
64 declareProperty(
"WriteNNinformation", m_writeNNinformation =
true);
65 declareProperty(
"WriteRDOinformation", m_writeRDOinformation =
true);
66 declareProperty(
"WriteExtendedPRDinformation", m_writeExtendedPRDinformation =
false);
69 declareProperty(
"SiClusterContainer", m_clustercontainer_key =
"PixelClusters");
70 declareProperty(
"MC_SDOs", m_SDOcontainer_key =
"PixelSDO_Map");
71 declareProperty(
"MC_Hits", m_sihitContainer_key =
"PixelHits");
72 declareProperty(
"PRD_MultiTruth", m_multiTruth_key =
"PRD_MultiTruthPixel");
74 declareProperty(
"OutputClusterContainer", m_write_xaod_key =
"PixelClusters");
77 declare(m_write_xaod_key);
78 declare(m_write_offsets);
122 return StatusCode::SUCCESS;
132 const EventContext& ctx = Gaudi::Hive::currentContext();
139 return StatusCode::FAILURE;
146 if (prdmtCollHandle.
isValid()) {
147 prdmtColl = &*prdmtCollHandle;
151 if (truthParticleLinksHandle.isValid()) {
152 truth_particle_links = truthParticleLinksHandle.cptr();
160 if (sdoCollectionHandle.
isValid()) {
161 sdoCollection = &*sdoCollectionHandle;
163 ATH_MSG_WARNING(
"SDO information requested, but SDO collection not available!");
168 bool foundSplitProbContainer =
false;
171 if (!splitProbContainer.
isValid()) {
174 foundSplitProbContainer =
true;
180 if (siHitCollectionHandle.
isValid()) {
181 for (
const SiHit& siHit: *siHitCollectionHandle) {
183 if (!siHit.isPixel())
continue;
186 siHit.getLayerDisk(),
187 siHit.getPhiModule(),
188 siHit.getEtaModule()));
191 siHits[wafer_hash].push_back(&siHit);
194 ATH_MSG_WARNING(
"SiHit information requested, but SiHit collection not available!");
201 if (!calibData_handle.
isValid()) {
204 calibData=calibData_handle.
cptr();
209 ATH_CHECK(xaod.
record(std::make_unique<xAOD::TrackMeasurementValidationContainer>(),
210 std::make_unique<xAOD::TrackMeasurementValidationAuxContainer>()));
215 unsigned int have_truth_link=0
u;
216 unsigned int missing_truth_particle=0
u;
217 unsigned int missing_parent_particle=0
u;
225 std::unordered_map< unsigned int , std::vector<unsigned int> > cluster_map;
229 (*offsets)[clusterCollection->identifyHash()] =
counter;
232 if( clusterCollection->empty() )
continue;
245 unsigned int cluster_idx = xaod->
size();
263 if(localCov.size() == 1){
266 }
else if(localCov.size() == 4){
275 std::vector< uint64_t > rdoIdentifierList;
276 rdoIdentifierList.reserve(prd->rdoList().size());
277 int rowmin=9999;
int rowmax=-9999;
278 int colmin=9999;
int colmax=-9999;
279 for(
const auto &hitIdentifier : prd->rdoList() ){
280 rdoIdentifierList.push_back( hitIdentifier.get_compact() );
284 if(rowmin >
row) rowmin =
row;
285 if(rowmax <
row) rowmax =
row;
286 if(colmin >
col) colmin =
col;
287 if(colmax <
col) colmax =
col;
297 AUXDATA(xprd,
int,phi_module) = the_phi ;
298 AUXDATA(xprd,
int,eta_module) = the_eta ;
302 cluster_map[ makeKey(the_phi, the_eta, the_layer)].push_back(cluster_idx);
307 AUXDATA(xprd,
int,nRDO) = (
int)prd->rdoList().size();
310 AUXDATA(xprd,
int,ToT) = prd->totalToT();
311 AUXDATA(xprd,
int,LVL1A) = prd->LVL1A();
314 AUXDATA(xprd,
char,gangedPixel) = (
char)prd->gangedPixel();
330 float deplVoltage = 0.0;
333 AUXDATA(xprd,
float,DepletionVoltage) = deplVoltage;
366 AUXDATA(xprd,
float,omegax) = prd->omegax();
367 AUXDATA(xprd,
float,omegay) = prd->omegay();
372 auto range{prdmtColl->equal_range(clusterId)};
373 if (truth_particle_links) {
374 std::vector<unsigned int> tp_indices;
377 if (a_truth_particle_link) {
379 if (truth_particle) {
381 tp_indices.push_back(
static_cast<int>(truth_particle->
index()));
384 ++missing_parent_particle;
389 ++missing_truth_particle;
393 AUXDATA(xprd,std::vector<unsigned int>, truth_index) = tp_indices;
403 std::vector< std::vector< int > > sdo_tracks;
425 for (
auto clusItr = xaod->
begin(); clusItr != xaod->
end(); ++clusItr ) {
426 AUXDATA(*clusItr,
char,broken) =
false;
436 for (
auto clusItr = xaod->
begin(); clusItr != xaod->
end(); ++clusItr)
443 for (
unsigned int cluster_idx : cluster_idx_list) {
444 auto pixelCluster2 = xaod->
at(cluster_idx);
445 if ( acc_layer(*pixelCluster2) !=
layer )
447 if ( acc_eta_module(*
pixelCluster) != acc_eta_module(*pixelCluster2) )
449 if ( acc_phi_module(*
pixelCluster) != acc_phi_module(*pixelCluster2) )
452 std::vector<int> barcodes2 = acc_sihit_barcode(*pixelCluster2);
455 if (
std::find(barcodes2.begin(), barcodes2.end(), bc ) == barcodes2.end())
continue;
458 acc_broken(*pixelCluster2) =
true;
468 return StatusCode::SUCCESS;
476 std::vector<int> sdo_word;
477 std::vector< std::vector< int > > sdo_depositsBarcode;
478 std::vector< std::vector< float > > sdo_depositsEnergy;
480 for(
const auto &hitIdentifier : prd->
rdoList() ){
481 auto pos = sdoCollection.find(hitIdentifier);
482 if(
pos == sdoCollection.end() )
continue;
483 sdo_word.push_back(
pos->second.word() ) ;
485 std::vector<float> sdoDepEnergy(
pos->second.getdeposits().size());
486 unsigned int nDepos{0};
487 for (
auto& deposit:
pos->second.getdeposits()) {
488 if (deposit.first) sdoDepBC[nDepos] =
HepMC::barcode(deposit.first);
490 sdoDepEnergy[nDepos] = deposit.second;
493 sdo_depositsBarcode.push_back( sdoDepBC );
494 sdo_depositsEnergy.push_back( sdoDepEnergy );
496 AUXDATA(xprd,std::vector<int>,sdo_words) = sdo_word;
497 AUXDATA(xprd,std::vector< std::vector<int> >,sdo_depositsBarcode) = sdo_depositsBarcode;
498 AUXDATA(xprd,std::vector< std::vector<float> >,sdo_depositsEnergy) = sdo_depositsEnergy;
500 return sdo_depositsBarcode;
507 const std::vector<SiHit> & matchingHits )
const
510 int numHits = matchingHits.size();
512 std::vector<float> sihit_energyDeposit(numHits,0);
513 std::vector<float> sihit_meanTime(numHits,0);
514 std::vector<int> sihit_barcode(numHits,0);
515 std::vector<int> sihit_pdgid(numHits,0);
517 std::vector<float> sihit_startPosX(numHits,0);
518 std::vector<float> sihit_startPosY(numHits,0);
519 std::vector<float> sihit_startPosZ(numHits,0);
521 std::vector<float> sihit_endPosX(numHits,0);
522 std::vector<float> sihit_endPosY(numHits,0);
523 std::vector<float> sihit_endPosZ(numHits,0);
528 for (
const auto& sihit : matchingHits ) {
529 sihit_energyDeposit[hitNumber] = sihit.energyLoss() ;
530 sihit_meanTime[hitNumber] = sihit.meanTime() ;
534 sihit_pdgid[hitNumber] = HMPL->pdg_id();
538 const HepGeom::Point3D<double>& startPos=sihit.localStartPosition();
541 sihit_startPosX[hitNumber] =
pos[0];
542 sihit_startPosY[hitNumber] =
pos[1];
543 sihit_startPosZ[hitNumber] = startPos.x();
546 const HepGeom::Point3D<double>& endPos=sihit.localEndPosition();
548 sihit_endPosX[hitNumber] =
pos[0];
549 sihit_endPosY[hitNumber] =
pos[1];
550 sihit_endPosZ[hitNumber] = endPos.x();
555 AUXDATA(xprd,std::vector<float>,sihit_energyDeposit) = sihit_energyDeposit;
556 AUXDATA(xprd,std::vector<float>,sihit_meanTime) = sihit_meanTime;
557 AUXDATA(xprd,std::vector<int>,sihit_barcode) = sihit_barcode;
558 AUXDATA(xprd,std::vector<int>,sihit_pdgid) = sihit_pdgid;
560 AUXDATA(xprd,std::vector<float>,sihit_startPosX) = sihit_startPosX;
561 AUXDATA(xprd,std::vector<float>,sihit_startPosY) = sihit_startPosY;
562 AUXDATA(xprd,std::vector<float>,sihit_startPosZ) = sihit_startPosZ;
564 AUXDATA(xprd,std::vector<float>,sihit_endPosX) = sihit_endPosX;
565 AUXDATA(xprd,std::vector<float>,sihit_endPosY) = sihit_endPosY;
566 AUXDATA(xprd,std::vector<float>,sihit_endPosZ) = sihit_endPosZ;
577 const std::vector<const SiHit*>* sihits,
578 std::vector< std::vector< int > > & trkBCs )
const
580 ATH_MSG_VERBOSE(
"Got " << sihits->size() <<
" SiHits to look through" );
581 std::vector<SiHit> matchingHits;
588 std::vector<const SiHit* > multiMatchingHits;
590 for (
const SiHit* siHit : *sihits) {
596 HepGeom::Point3D<double> averagePosition = siHit->localStartPosition() + siHit->localEndPosition();
597 averagePosition *= 0.5;
601 for(
const auto &hitIdentifier : prd->
rdoList() ){
607 multiMatchingHits.push_back(siHit);
615 for (
const auto& barcodeSDOColl : trkBCs ) {
616 if (
std::find(barcodeSDOColl.begin(),barcodeSDOColl.end(),bc) == barcodeSDOColl.end() )
continue;
617 multiMatchingHits.push_back(siHit);
625 ATH_MSG_DEBUG(
"Found " << multiMatchingHits.size() <<
" SiHit " );
626 for ( ; siHitIter != multiMatchingHits.end(); ++siHitIter) {
627 const SiHit* lowestXPos = *siHitIter;
628 const SiHit* highestXPos = *siHitIter;
632 std::vector<const SiHit* > ajoiningHits;
633 ajoiningHits.push_back( *siHitIter );
635 siHitIter2 = siHitIter+1;
636 while ( siHitIter2 != multiMatchingHits.end() ) {
645 if (std::abs((highestXPos->
localEndPosition().x()-(*siHitIter2)->localStartPosition().x()))<0.00005 &&
646 std::abs((highestXPos->
localEndPosition().y()-(*siHitIter2)->localStartPosition().y()))<0.00005 &&
647 std::abs((highestXPos->
localEndPosition().z()-(*siHitIter2)->localStartPosition().z()))<0.00005 )
649 highestXPos = *siHitIter2;
650 ajoiningHits.push_back( *siHitIter2 );
653 siHitIter2 = multiMatchingHits.erase( siHitIter2 );
654 }
else if (std::abs((lowestXPos->
localStartPosition().x()-(*siHitIter2)->localEndPosition().x()))<0.00005 &&
655 std::abs((lowestXPos->
localStartPosition().y()-(*siHitIter2)->localEndPosition().y()))<0.00005 &&
656 std::abs((lowestXPos->
localStartPosition().z()-(*siHitIter2)->localEndPosition().z()))<0.00005)
658 lowestXPos = *siHitIter2;
659 ajoiningHits.push_back( *siHitIter2 );
662 siHitIter2 = multiMatchingHits.erase( siHitIter2 );
668 if( ajoiningHits.size() == 0){
671 }
else if(ajoiningHits.size() == 1){
673 matchingHits.push_back( *ajoiningHits[0] );
677 ATH_MSG_DEBUG(
"Merging " << ajoiningHits.size() <<
" SiHits together." );
682 for(
const auto& siHit : ajoiningHits){
683 energyDep += siHit->energyLoss();
684 time += siHit->meanTime();
694 (*siHitIter)->getBarrelEndcap(),
695 (*siHitIter)->getLayerDisk(),
696 (*siHitIter)->getEtaModule(),
697 (*siHitIter)->getPhiModule(),
698 (*siHitIter)->getSide() );
699 ATH_MSG_DEBUG(
"Finished Merging " << ajoiningHits.size() <<
" SiHits together." );
716 const std::vector<Identifier>& rdos =
pixelCluster->rdoList();
718 const std::vector<float> &chList =
pixelCluster->chargeList();
719 const std::vector<int> &totList =
pixelCluster->totList();
723 std::vector<int> etaIndexList;
724 std::vector<int> phiIndexList;
725 std::vector<float> CTerm;
726 std::vector<float> ATerm;
727 std::vector<float> ETerm;
732 std::vector<Identifier>::const_iterator rdosBegin = rdos.begin();
733 std::vector<Identifier>::const_iterator rdosEnd = rdos.end();
735 ATH_MSG_VERBOSE(
" Putting together the n. " << rdos.size() <<
" rdos into a matrix.");
737 phiIndexList.reserve( rdos.size());
738 etaIndexList.reserve( rdos.size());
739 CTerm.reserve( rdos.size());
740 ATerm.reserve( rdos.size());
741 ETerm.reserve( rdos.size());
742 for (; rdosBegin!= rdosEnd; ++rdosBegin)
762 AUXDATA(xprd, std::vector<int>,rdo_phi_pixel_index) = phiIndexList;
763 AUXDATA(xprd, std::vector<int>,rdo_eta_pixel_index) = etaIndexList;
764 AUXDATA(xprd, std::vector<float>,rdo_charge) = chList;
765 AUXDATA(xprd, std::vector<int>,rdo_tot) = totList;
767 AUXDATA(xprd, std::vector<float>,rdo_Cterm) = CTerm;
768 AUXDATA(xprd, std::vector<float>,rdo_Aterm) = ATerm;
769 AUXDATA(xprd, std::vector<float>,rdo_Eterm) = ETerm;
776 const unsigned int sizeX,
const unsigned int sizeY )
const
789 ATH_MSG_WARNING(
"PixelModuleDesign was not retrieved in function 'addNNInformation'");
792 const std::vector<Identifier>& rdos =
pixelCluster->rdoList();
794 const std::vector<float>& chList =
pixelCluster->chargeList();
795 const std::vector<int>& totList =
pixelCluster->totList();
803 int phiPixelIndexMin, phiPixelIndexMax, etaPixelIndexMin, etaPixelIndexMax;
806 if (!cellIdWeightedPosition.
isValid())
811 int etaPixelIndexWeightedPosition=cellIdWeightedPosition.
etaIndex();
812 int phiPixelIndexWeightedPosition=cellIdWeightedPosition.
phiIndex();
815 ATH_MSG_DEBUG(
" weighted pos phiPixelIndex: " << phiPixelIndexWeightedPosition <<
" etaPixelIndex: " << etaPixelIndexWeightedPosition );
833 double localEtaPixelIndexWeightedPosition =
w.xEta();
834 double localPhiPixelIndexWeightedPosition =
w.xPhi();
836 int centralIndexX=(sizeX-1)/2;
837 int centralIndexY=(sizeY-1)/2;
843 if (abs(phiPixelIndexWeightedPosition-phiPixelIndexMin)>centralIndexX ||
844 abs(phiPixelIndexWeightedPosition-phiPixelIndexMax)>centralIndexX)
846 ATH_MSG_DEBUG(
" Cluster too large phiPixelIndexMin " << phiPixelIndexMin <<
" phiPixelIndexMax " << phiPixelIndexMax <<
" centralX " << centralIndexX);
850 if (abs(etaPixelIndexWeightedPosition-etaPixelIndexMin)>centralIndexY ||
851 abs(etaPixelIndexWeightedPosition-etaPixelIndexMax)>centralIndexY)
853 ATH_MSG_DEBUG(
" Cluster too large etaPixelIndexMin" << etaPixelIndexMin <<
" etaPixelIndexMax " << etaPixelIndexMax <<
" centralY " << centralIndexY);
857 std::vector< std::vector<float> > matrixOfToT (sizeX, std::vector<float>(sizeY,0) );
858 std::vector< std::vector<float> > matrixOfCharge(sizeX, std::vector<float>(sizeY,0));
859 std::vector<float> vectorOfPitchesY(sizeY,0.4);
863 std::vector<Identifier>::const_iterator rdosBegin = rdos.begin();
864 std::vector<Identifier>::const_iterator rdosEnd = rdos.end();
865 auto charge = chList.begin();
866 auto tot = totList.begin();
868 ATH_MSG_VERBOSE(
" Putting together the n. " << rdos.size() <<
" rdos into a matrix.");
870 for (; rdosBegin!= rdosEnd; ++rdosBegin)
876 if (
charge != chList.end()){
879 if (absphiPixelIndex <0 || absphiPixelIndex >= (
int)sizeX)
881 ATH_MSG_DEBUG(
" problem with index: " << absphiPixelIndex <<
" min: " << 0 <<
" max: " << sizeX);
885 if (absetaPixelIndex <0 || absetaPixelIndex >= (
int)sizeY)
887 ATH_MSG_DEBUG(
" problem with index: " << absetaPixelIndex <<
" min: " << 0 <<
" max: " << sizeY);
893 float pitchY = diodeParameters.
width().
xEta();
895 if ( (not totList.empty()) && tot != totList.end()) {
896 matrixOfToT[absphiPixelIndex][absetaPixelIndex] =*tot;
898 }
else matrixOfToT[absphiPixelIndex][absetaPixelIndex] = -1;
900 if ( (not chList.empty()) &&
charge != chList.end()){
901 matrixOfCharge[absphiPixelIndex][absetaPixelIndex]=*
charge;
903 }
else matrixOfCharge[absphiPixelIndex][absetaPixelIndex] = -1;
907 vectorOfPitchesY[absetaPixelIndex]=pitchY;
921 trackDir.normalize();
928 float trkphicomp = trackDir.dot(module_phiax);
929 float trketacomp = trackDir.dot(module_etaax);
930 float trknormcomp = trackDir.dot(module_normal);
931 double bowphi = atan2(trkphicomp,trknormcomp);
932 double boweta = atan2(trketacomp,trknormcomp);
934 if(bowphi > TMath::Pi()/2) bowphi -= TMath::Pi();
935 if(bowphi < -TMath::Pi()/2) bowphi += TMath::Pi();
936 int readoutside = design->readoutSide();
942 if (boweta>TMath::Pi()/2.) boweta-=TMath::Pi();
943 if (boweta<-TMath::Pi()/2.) boweta+=TMath::Pi();
947 ATH_MSG_VERBOSE(
" PhiPixelIndexWeightedPosition: " << phiPixelIndexWeightedPosition <<
" EtaPixelIndexWeightedPosition: " << etaPixelIndexWeightedPosition );
950 std::vector<float> vectorOfCharge(sizeX*sizeY,0);
951 std::vector<float> vectorOfToT(sizeX*sizeY,0);
953 for (
unsigned int u=0;
u<sizeX;
u++)
955 for (
unsigned int s=0;
s<sizeY;
s++)
957 vectorOfToT[
counter] = matrixOfToT[
u][
s];
958 vectorOfCharge[
counter] = matrixOfCharge[
u][
s];
967 AUXDATA(xprd,
int, NN_sizeX) = sizeX;
968 AUXDATA(xprd,
int, NN_sizeY) = sizeY;
971 AUXDATA(xprd,
float, NN_thetaBS) = boweta;
973 AUXDATA(xprd, std::vector<float>, NN_matrixOfToT) = vectorOfToT;
974 AUXDATA(xprd, std::vector<float>, NN_matrixOfCharge) = vectorOfCharge;
975 AUXDATA(xprd, std::vector<float>, NN_vectorOfPitchesY) = vectorOfPitchesY;
978 AUXDATA(xprd,
int, NN_etaPixelIndexWeightedPosition) = etaPixelIndexWeightedPosition;
979 AUXDATA(xprd,
int, NN_phiPixelIndexWeightedPosition) = phiPixelIndexWeightedPosition;
981 AUXDATA(xprd,
float, NN_localEtaPixelIndexWeightedPosition) = localEtaPixelIndexWeightedPosition;
982 AUXDATA(xprd,
float, NN_localPhiPixelIndexWeightedPosition) = localPhiPixelIndexWeightedPosition;
989 const std::vector<SiHit> & matchingHits )
const
1027 ATH_MSG_WARNING(
"PixelModuleDesign was not retrieved in function 'addNNTruthInfo'");
1032 unsigned hitNumber(0);
1033 for(
const auto& siHit : matchingHits ){
1035 HepGeom::Point3D<double> averagePosition = (siHit.localStartPosition() + siHit.localEndPosition()) * 0.5;
1037 ATH_MSG_VERBOSE(
"Truth Part X: " << averagePosition.y() <<
" shift " << shift <<
" Y: " << averagePosition.z() );
1040 float YposC = averagePosition.y()-shift;
1042 if (std::abs(YposC)>design->width()/2 &&
1043 std::abs(averagePosition.y())<design->width()/2)
1045 if (YposC>design->width()/2)
1047 YposC=design->width()/2-1
e-6;
1048 }
else if (YposC<-design->
width()/2)
1050 YposC=-design->width()/2+1
e-6;
1054 positionsX[hitNumber] = YposC;
1055 positionsY[hitNumber] = averagePosition.z();
1057 HepGeom::Point3D<double> deltaPosition = siHit.localEndPosition() - siHit.localStartPosition();
1059 pathlengthX[hitNumber] = deltaPosition.y();
1060 pathlengthY[hitNumber] = deltaPosition.z();
1061 pathlengthZ[hitNumber] = deltaPosition.x();
1066 InDetDD::SiCellId cellIdOfTruthPosition = design->cellIdOfPosition(siLocalTruthPosition);
1072 int truthEtaIndex = cellIdOfTruthPosition.
etaIndex();
1073 int truthPhiIndex = cellIdOfTruthPosition.
phiIndex();
1076 double pitchY = diodeParameters.
width().
xEta();
1077 double pitchX = diodeParameters.
width().
xPhi();
1094 double pixelCenterY = siLocalPositionCenter.
xEta();
1095 double pixelCenterX = siLocalPositionCenter.
xPhi();
1101 double truthIndexY = truthEtaIndex + (siLocalTruthPosition[
Trk::distEta] - pixelCenterY)/pitchY;
1102 double truthIndexX = truthPhiIndex + (siLocalTruthPosition[
Trk::distPhi] - pixelCenterX)/pitchX;
1105 positions_indexX[hitNumber] = truthIndexX - cellIdWeightedPosition.
phiIndex();
1106 positions_indexY[hitNumber] = truthIndexY - cellIdWeightedPosition.
etaIndex();
1108 HepGeom::Point3D<double> diffPositions = (siHit.localEndPosition() - siHit.localStartPosition());
1109 double bowphi = std::atan2( diffPositions.y(), diffPositions.x() );
1113 theta[hitNumber] = std::atan2(diffPositions.z() ,diffPositions.x());
1117 int readoutside = design->readoutSide();
1123 pdgid[hitNumber] =
particle->pdg_id();
1130 const auto& mother_of_particle=
vertex->particles_in().front();
1132 motherPdgid[hitNumber] = mother_of_particle->pdg_id();
1136 if(
vertex->particles_in_const_begin() !=
vertex->particles_in_const_end() ){
1138 motherPdgid[hitNumber] = (*
vertex->particles_in_const_begin())->pdg_id();
1143 chargeDep[hitNumber] = siHit.energyLoss() ;
1149 AUXDATA(xprd, std::vector<float>, NN_positionsX) = positionsX;
1150 AUXDATA(xprd, std::vector<float>, NN_positionsY) = positionsY;
1152 AUXDATA(xprd, std::vector<float>, NN_positions_indexX) = positions_indexX;
1153 AUXDATA(xprd, std::vector<float>, NN_positions_indexY) = positions_indexY;
1156 AUXDATA(xprd, std::vector<float>, NN_phi) =
phi;
1159 AUXDATA(xprd, std::vector<int>, NN_pdgid) = pdgid;
1160 AUXDATA(xprd, std::vector<float>, NN_energyDep) = chargeDep;
1161 AUXDATA(xprd, std::vector<float>, NN_trueP) = truep;
1163 AUXDATA(xprd, std::vector<int>, NN_motherBarcode) = motherBarcode;
1164 AUXDATA(xprd, std::vector<int>, NN_motherPdgid) = motherPdgid;
1167 AUXDATA(xprd, std::vector<float>, NN_pathlengthX) = pathlengthX;
1168 AUXDATA(xprd, std::vector<float>, NN_pathlengthY) = pathlengthY;
1169 AUXDATA(xprd, std::vector<float>, NN_pathlengthZ) = pathlengthZ;
1178 int *rphiPixelIndexMin,
1179 int *rphiPixelIndexMax,
1180 int *retaPixelIndexMin,
1181 int *retaPixelIndexMax )
const
1192 ATH_MSG_WARNING(
"PixelModuleDesign was not retrieved in function 'getCellIdWeightedPosition'");
1195 const std::vector<Identifier>& rdos =
pixelCluster->rdoList();
1198 const std::vector<float>& chList =
pixelCluster->chargeList();
1201 std::vector<Identifier>::const_iterator rdosBegin = rdos.begin();
1202 std::vector<Identifier>::const_iterator rdosEnd = rdos.end();
1204 auto charge = chList.begin();
1207 double sumOfCharge=0;
1209 int phiPixelIndexMin = 99999;
1210 int phiPixelIndexMax = -99999;
1211 int etaPixelIndexMin = 99999;
1212 int etaPixelIndexMax = -99999;
1214 for (; rdosBegin!= rdosEnd; ++rdosBegin, ++
charge)
1221 ATH_MSG_VERBOSE(
" Adding pixel phiPixelIndex: " << phiPixelIndex <<
" etaPixelIndex: " << etaPixelIndex <<
" charge: " << *
charge );
1239 sumOfWeightedPositions += (*charge)*siLocalPosition;
1240 sumOfCharge += (*charge);
1242 if (phiPixelIndex < phiPixelIndexMin)
1243 phiPixelIndexMin = phiPixelIndex;
1245 if (phiPixelIndex > phiPixelIndexMax)
1246 phiPixelIndexMax = phiPixelIndex;
1248 if (etaPixelIndex < etaPixelIndexMin)
1249 etaPixelIndexMin = etaPixelIndex;
1251 if (etaPixelIndex > etaPixelIndexMax)
1252 etaPixelIndexMax = etaPixelIndex;
1255 sumOfWeightedPositions /= sumOfCharge;
1257 ATH_MSG_VERBOSE (
"Wighted position: Row = " << sumOfWeightedPositions.
xRow() <<
", Col = " << sumOfWeightedPositions.
xColumn() );
1259 if(rphiPixelIndexMin) *rphiPixelIndexMin = phiPixelIndexMin;
1260 if(rphiPixelIndexMax) *rphiPixelIndexMax = phiPixelIndexMax;
1261 if(retaPixelIndexMin) *retaPixelIndexMin = etaPixelIndexMin;
1262 if(retaPixelIndexMax) *retaPixelIndexMax = etaPixelIndexMax;
1267 InDetDD::SiCellId cellIdWeightedPosition=design->cellIdOfPosition(sumOfWeightedPositions);
1270 return cellIdWeightedPosition;
1287 return StatusCode::SUCCESS;