@@ -77,6 +77,13 @@ class SiPixelChargeReweightingAlgorithm {
7777 const TrackerTopology* tTopo,
7878 CLHEP::HepRandomEngine* engine);
7979
80+ template <class AmplitudeType , typename SignalType>
81+ bool lateSignalReweight (const PixelGeomDetUnit* pixdet,
82+ std::vector<PixelDigi>& digis,
83+ PixelSimHitExtraInfoLite& loopTempSH,
84+ SignalType& theNewDigiSignal,
85+ const TrackerTopology* tTopo,
86+ CLHEP::HepRandomEngine* engine); // could be templated in future ?
8087private:
8188 // Internal typedef
8289 typedef GloballyPositioned<double > Frame;
@@ -841,4 +848,219 @@ bool SiPixelChargeReweightingAlgorithm::lateSignalReweight(const PixelGeomDetUni
841848 return true ;
842849}
843850
851+ template <class AmplitudeType , typename SignalType>
852+ bool SiPixelChargeReweightingAlgorithm::lateSignalReweight (const PixelGeomDetUnit* pixdet,
853+ std::vector<PixelDigi>& digis,
854+ PixelSimHitExtraInfoLite& loopTempSH,
855+ SignalType& theNewDigiSignal,
856+ const TrackerTopology* tTopo,
857+ CLHEP::HepRandomEngine* engine) {
858+ uint32_t detID = pixdet->geographicalId ().rawId ();
859+ const PixelTopology* topol = &pixdet->specificTopology ();
860+
861+ if (UseReweighting) {
862+ edm::LogError (" SiPixelChargeReweightingAlgorithm" ) << " ******************************** \n " ;
863+ edm::LogError (" SiPixelChargeReweightingAlgorithm" ) << " ******************************** \n " ;
864+ edm::LogError (" SiPixelChargeReweightingAlgorithm" ) << " ******************************** \n " ;
865+ edm::LogError (" SiPixelChargeReweightingAlgorithm" ) << " ***** INCONSISTENCY !!! ***** \n " ;
866+ edm::LogError (" SiPixelChargeReweightingAlgorithm" )
867+ << " applyLateReweighting_ and UseReweighting can not be true at the same time for PU ! \n " ;
868+ edm::LogError (" SiPixelChargeReweightingAlgorithm" ) << " ---> DO NOT APPLY CHARGE REWEIGHTING TWICE !!! \n " ;
869+ edm::LogError (" SiPixelChargeReweightingAlgorithm" ) << " ******************************** \n " ;
870+ edm::LogError (" SiPixelChargeReweightingAlgorithm" ) << " ******************************** \n " ;
871+ return false ;
872+ }
873+
874+ SignalType theDigiSignal;
875+
876+ int irow_min = topol->nrows ();
877+ int irow_max = 0 ;
878+ int icol_min = topol->ncolumns ();
879+ int icol_max = 0 ;
880+
881+ float chargeBefore = 0 ;
882+ float chargeAfter = 0 ;
883+
884+ // loop on digis
885+ std::vector<PixelDigi>::const_iterator loopDigi;
886+ for (loopDigi = digis.begin (); loopDigi != digis.end (); ++loopDigi) {
887+ unsigned int chan = loopDigi->channel ();
888+ if (loopTempSH.isInTheList (chan)) {
889+ float corresponding_charge = loopDigi->adc ();
890+ if constexpr (std::is_same_v<AmplitudeType, digitizerUtility::Ph2Amplitude>) {
891+ edm::LogError (" SiPixelChargeReweightingAlgorithm" )
892+ << " Phase-2 late charge reweighting not supported (not sure we need it at all)" ;
893+ return false ;
894+ } else {
895+ theDigiSignal[chan] += AmplitudeType (corresponding_charge, corresponding_charge);
896+ }
897+ chargeBefore += corresponding_charge;
898+ if (loopDigi->row () < irow_min)
899+ irow_min = loopDigi->row ();
900+ if (loopDigi->row () > irow_max)
901+ irow_max = loopDigi->row ();
902+ if (loopDigi->column () < icol_min)
903+ icol_min = loopDigi->column ();
904+ if (loopDigi->column () > icol_max)
905+ icol_max = loopDigi->column ();
906+ }
907+ }
908+ // end loop on digis
909+
910+ LocalVector direction = loopTempSH.exitPoint () - loopTempSH.entryPoint ();
911+ LocalPoint hitEntryPoint = loopTempSH.entryPoint ();
912+ float trajectoryScaleToPosition = std::abs (hitEntryPoint.z () / direction.z ());
913+ LocalPoint hitPosition = hitEntryPoint + trajectoryScaleToPosition * direction;
914+
915+ // addition as now the hit himself is not available
916+ // https://github.com/cms-sw/cmssw/blob/master/SimDataFormats/TrackingHit/interface/PSimHit.h#L52
917+ LocalPoint hitLocalPosition = hitEntryPoint + 0.5 * direction;
918+ MeasurementPoint hitPositionPixel = topol->measurementPosition (hitLocalPosition);
919+ std::pair<int , int > hitPixel =
920+ std::pair<int , int >(int (floor (hitPositionPixel.x ())), int (floor (hitPositionPixel.y ())));
921+
922+ MeasurementPoint originPixel = MeasurementPoint (hitPixel.first - THX + 0.5 , hitPixel.second - THY + 0.5 );
923+ LocalPoint origin = topol->localPosition (originPixel);
924+
925+ MeasurementPoint hitEntryPointPixel = topol->measurementPosition (loopTempSH.entryPoint ());
926+ MeasurementPoint hitExitPointPixel = topol->measurementPosition (loopTempSH.exitPoint ());
927+ std::pair<int , int > entryPixel =
928+ std::pair<int , int >(int (floor (hitEntryPointPixel.x ())), int (floor (hitEntryPointPixel.y ())));
929+ std::pair<int , int > exitPixel =
930+ std::pair<int , int >(int (floor (hitExitPointPixel.x ())), int (floor (hitExitPointPixel.y ())));
931+
932+ int hitcol_min, hitcol_max, hitrow_min, hitrow_max;
933+ if (entryPixel.first > exitPixel.first ) {
934+ hitrow_min = exitPixel.first ;
935+ hitrow_max = entryPixel.first ;
936+ } else {
937+ hitrow_min = entryPixel.first ;
938+ hitrow_max = exitPixel.first ;
939+ }
940+
941+ if (entryPixel.second > exitPixel.second ) {
942+ hitcol_min = exitPixel.second ;
943+ hitcol_max = entryPixel.second ;
944+ } else {
945+ hitcol_min = entryPixel.second ;
946+ hitcol_max = exitPixel.second ;
947+ }
948+
949+ if (!(irow_min <= hitrow_max && irow_max >= hitrow_min && icol_min <= hitcol_max && icol_max >= hitcol_min)) {
950+ // The clusters do not have an overlap, hence the hit is NOT reweighted
951+ return false ;
952+ }
953+
954+ track[0 ] = (hitPosition.x () - origin.x ()) * cmToMicrons;
955+ track[1 ] = (hitPosition.y () - origin.y ()) * cmToMicrons;
956+ track[2 ] = 0 .0f ; // Middle of sensor is origin for Z-axis
957+ track[3 ] = direction.x ();
958+ track[4 ] = direction.y ();
959+ track[5 ] = direction.z ();
960+
961+ array_2d pixrewgt (boost::extents[TXSIZE][TYSIZE]);
962+
963+ /*
964+ for (int row = 0; row < TXSIZE; ++row) {
965+ for (int col = 0; col < TYSIZE; ++col) {
966+ pixrewgt[row][col] = 0;
967+ }
968+ }
969+ */
970+
971+ for (int row = 0 ; row < TXSIZE; ++row) {
972+ xdouble[row] = topol->isItBigPixelInX (hitPixel.first + row - THX);
973+ }
974+
975+ for (int col = 0 ; col < TYSIZE; ++col) {
976+ ydouble[col] = topol->isItBigPixelInY (hitPixel.second + col - THY);
977+ }
978+
979+ for (int row = 0 ; row < TXSIZE; ++row) {
980+ for (int col = 0 ; col < TYSIZE; ++col) {
981+ // Fill charges into 21x13 Pixel Array with hitPixel in centre
982+ pixrewgt[row][col] =
983+ theDigiSignal[PixelDigi::pixelToChannel (hitPixel.first + row - THX, hitPixel.second + col - THY)];
984+ }
985+ }
986+
987+ if (PrintClusters) {
988+ LogDebug (" SiPixelChargeReweightingAlgorithm" ) << " Cluster before reweighting: " ;
989+ printCluster (pixrewgt);
990+ }
991+
992+ int ierr;
993+ // for unirradiated: 2nd argument is IDden
994+ // for irradiated: 2nd argument is IDnum
995+ int ID1 = dbobject_num->getTemplateID (detID);
996+ int ID0 = dbobject_den->getTemplateID (detID);
997+
998+ if (ID0 == ID1) {
999+ LogDebug (" SiPixelChargeReweightingAlgorithm" ) << " same template for num and den " ;
1000+ return false ;
1001+ }
1002+ ierr = PixelTempRewgt2D (ID0, ID1, pixrewgt);
1003+ if (ierr != 0 ) {
1004+ edm::LogError (" SiPixelChargeReweightingAlgorithm" ) << " Cluster Charge Reweighting did not work properly." ;
1005+ return false ;
1006+ }
1007+ if (PrintClusters) {
1008+ LogDebug (" SiPixelChargeReweightingAlgorithm" ) << " Cluster after reweighting: " ;
1009+ printCluster (pixrewgt);
1010+ }
1011+
1012+ for (int row = 0 ; row < TXSIZE; ++row) {
1013+ for (int col = 0 ; col < TYSIZE; ++col) {
1014+ float charge = 0 ;
1015+ charge = pixrewgt[row][col];
1016+ if ((hitPixel.first + row - THX) >= 0 && (hitPixel.first + row - THX) < topol->nrows () &&
1017+ (hitPixel.second + col - THY) >= 0 && (hitPixel.second + col - THY) < topol->ncolumns () && charge > 0 ) {
1018+ chargeAfter += charge;
1019+ if constexpr (std::is_same_v<AmplitudeType, digitizerUtility::Ph2Amplitude>) {
1020+ edm::LogError (" SiPixelChargeReweightingAlgorithm" )
1021+ << " Phase-2 late charge reweighting not supported (not sure we need it at all)" ;
1022+ return false ;
1023+ } else {
1024+ theNewDigiSignal[PixelDigi::pixelToChannel (hitPixel.first + row - THX, hitPixel.second + col - THY)] +=
1025+ AmplitudeType (charge, charge);
1026+ }
1027+ }
1028+ }
1029+ }
1030+
1031+ if (chargeBefore != 0 . && chargeAfter == 0 .) {
1032+ return false ;
1033+ }
1034+
1035+ if (PrintClusters) {
1036+ LogDebug (" SiPixelChargeReweightingAlgorithm" )
1037+ << " Charges (before->after): " << chargeBefore << " -> " << chargeAfter;
1038+ LogDebug (" SiPixelChargeReweightingAlgorithm" ) << " Charge loss: " << (1 - chargeAfter / chargeBefore) * 100 << " %" ;
1039+ }
1040+
1041+ // need to store the digi out of the 21x13 box.
1042+ for (auto & i : theDigiSignal) {
1043+ // check if in the 21x13 box
1044+ int chanDigi = i.first ;
1045+ std::pair<int , int > ip = PixelDigi::channelToPixel (chanDigi);
1046+ int row_digi = ip.first ;
1047+ int col_digi = ip.second ;
1048+ int i_transformed_row = row_digi - hitPixel.first + THX;
1049+ int i_transformed_col = col_digi - hitPixel.second + THY;
1050+ if (i_transformed_row < 0 || i_transformed_row > TXSIZE || i_transformed_col < 0 || i_transformed_col > TYSIZE) {
1051+ // not in the box
1052+ if (chanDigi >= 0 && i.second > 0 ) {
1053+ if constexpr (std::is_same_v<AmplitudeType, digitizerUtility::Ph2Amplitude>) {
1054+ edm::LogError (" SiPixelChargeReweightingAlgorithm" )
1055+ << " Phase-2 late charge reweighting not supported (not sure we need it at all)" ;
1056+ return false ;
1057+ } else {
1058+ theNewDigiSignal[chanDigi] += AmplitudeType (i.second , i.second );
1059+ }
1060+ }
1061+ }
1062+ }
1063+
1064+ return true ;
1065+ }
8441066#endif
0 commit comments