File indexing completed on 2024-04-06 12:30:52
0001 #include "SimPPS/PPSPixelDigiProducer/interface/RPixChargeShare.h"
0002 #include "FWCore/ParameterSet/interface/FileInPath.h"
0003 #include <iostream>
0004 #include <fstream>
0005
0006 RPixChargeShare::RPixChargeShare(const edm::ParameterSet ¶ms, uint32_t det_id, const PPSPixelTopology &ppt)
0007 : det_id_(det_id) {
0008 verbosity_ = params.getParameter<int>("RPixVerbosity");
0009 signalCoupling_.clear();
0010 ChargeMapFile2E_[0] = params.getParameter<std::string>("ChargeMapFile2E");
0011 ChargeMapFile2E_[1] = params.getParameter<std::string>("ChargeMapFile2E_2X");
0012 ChargeMapFile2E_[2] = params.getParameter<std::string>("ChargeMapFile2E_2Y");
0013 ChargeMapFile2E_[3] = params.getParameter<std::string>("ChargeMapFile2E_2X2Y");
0014
0015 double coupling_constant_ = params.getParameter<double>("RPixCoupling");
0016
0017 signalCoupling_.push_back(coupling_constant_);
0018 signalCoupling_.push_back((1.0 - coupling_constant_) / 2);
0019
0020 no_of_pixels_ = ppt.getNoPixels();
0021
0022 double xMap, yMap;
0023 double chargeprobcollect;
0024 int xUpper[] = {75, 150, 75, 150};
0025 int yUpper[] = {50, 50, 100, 100};
0026 int ix, iy;
0027 for (int i = 0; i < 4; i++) {
0028 edm::FileInPath filename(ChargeMapFile2E_[i]);
0029 std::ifstream fChargeMap(filename.fullPath().c_str());
0030 if (fChargeMap.is_open()) {
0031 while (fChargeMap >> xMap >> yMap >> chargeprobcollect) {
0032 ix = int((xMap + xUpper[i]) / 5);
0033 iy = int((yMap + yUpper[i]) / 5);
0034 chargeMap2E_[i][ix][iy] = chargeprobcollect;
0035 }
0036 fChargeMap.close();
0037 } else
0038 throw cms::Exception("PPS RPixChargeShare") << "Charge map file not found";
0039 }
0040 }
0041
0042 std::map<unsigned short, double> RPixChargeShare::Share(const std::vector<RPixSignalPoint> &charge_map,
0043 const PPSPixelTopology &ppt) {
0044 std::map<unsigned short, double> thePixelChargeMap;
0045 if (verbosity_ > 1)
0046 edm::LogInfo("PPS") << "RPixChargeShare " << det_id_ << " : Clouds to be induced= " << charge_map.size();
0047
0048 double cH = 0;
0049
0050 for (std::vector<RPixSignalPoint>::const_iterator i = charge_map.begin(); i != charge_map.end(); ++i) {
0051 double hit_pos_x, hit_pos_y;
0052
0053 if (((*i).Position().x() + ppt.getSimXWidth() / 2.) < 0 ||
0054 ((*i).Position().x() + ppt.getSimXWidth() / 2.) > ppt.getSimXWidth()) {
0055 edm::LogInfo("PPS RPixChargeShare")
0056 << "**** Attention ((*i).Position().x()+simX_width_/2.)<0||((*i).Position().x()+simX_width_/2.)>simX_width ";
0057 edm::LogInfo("PPS RPixChargeShare") << "(*i).Position().x() = " << (*i).Position().x();
0058 continue;
0059 }
0060 if (((*i).Position().y() + ppt.getSimYWidth() / 2.) < 0 ||
0061 ((*i).Position().y() + ppt.getSimYWidth() / 2.) > ppt.getSimYWidth()) {
0062 edm::LogInfo("PPS RPixChargeShare")
0063 << "**** Attention ((*i).Position().y()+simY_width_/2.)<0||((*i).Position().y()+simY_width_/2.)>simY_width ";
0064 edm::LogInfo("PPS RPixChargeShare") << "(*i).Position().y() = " << (*i).Position().y();
0065 continue;
0066 }
0067
0068 PPSPixelTopology::PixelInfo relevant_pixels =
0069 ppt.getPixelsInvolved((*i).Position().x(), (*i).Position().y(), (*i).Sigma(), hit_pos_x, hit_pos_y);
0070 double effic = relevant_pixels.effFactor();
0071 if (effic < 1.e-4)
0072 continue;
0073
0074 unsigned short pixel_no = ppt.pixelIndex(relevant_pixels);
0075
0076 double charge_in_pixel = (*i).Charge() * effic;
0077
0078 cH += charge_in_pixel;
0079
0080 if (verbosity_ > 1)
0081 edm::LogInfo("PPS RPixChargeShare ") << "Efficiency in detector " << det_id_ << " and pixel no " << pixel_no
0082 << " : " << effic << " ch: " << charge_in_pixel << " CHtot: " << cH;
0083
0084 if (signalCoupling_[0] == 0.) {
0085 thePixelChargeMap[pixel_no] += charge_in_pixel;
0086 } else {
0087 int pixel_row = relevant_pixels.pixelRowNo();
0088 int pixel_col = relevant_pixels.pixelColNo();
0089 double pixel_lower_x = 0;
0090 double pixel_lower_y = 0;
0091 double pixel_upper_x = 0;
0092 double pixel_upper_y = 0;
0093 int psize = 0;
0094 ppt.pixelRange(pixel_row, pixel_col, pixel_lower_x, pixel_upper_x, pixel_lower_y, pixel_upper_y);
0095 double pixel_width_x = pixel_upper_x - pixel_lower_x;
0096 double pixel_width_y = pixel_upper_y - pixel_lower_y;
0097 double pixel_center_x = pixel_lower_x + (pixel_width_x) / 2.;
0098 double pixel_center_y = pixel_lower_y + (pixel_width_y) / 2.;
0099
0100
0101
0102
0103
0104
0105
0106
0107
0108
0109 if (pixel_row == 0) {
0110 pixel_width_x = ppt.getPitchSimX();
0111 pixel_center_x = pixel_lower_x + ppt.getPhysActiveEdgeDist() + ppt.getPitchSimX() / 2.;
0112 if (fabs((*i).Position().x() - pixel_center_x) > ppt.getPitchSimX() / 2.) {
0113 thePixelChargeMap[pixel_no] += charge_in_pixel;
0114 continue;
0115 }
0116 if (pixel_col == 0) {
0117 pixel_width_y = ppt.getPitchSimY();
0118 pixel_center_y = pixel_lower_y + ppt.getPhysActiveEdgeDist() + ppt.getPitchSimY() / 2.;
0119 if (fabs((*i).Position().y() - pixel_center_y) > ppt.getPitchSimY() / 2.) {
0120 thePixelChargeMap[pixel_no] += charge_in_pixel;
0121 continue;
0122 }
0123 }
0124 if (pixel_col == ppt.getNoPixelsSimY() - 1) {
0125 pixel_width_y = ppt.getPitchSimY();
0126 pixel_center_y = pixel_lower_y + ppt.getPitchSimY() / 2.;
0127 if (fabs((*i).Position().y() - pixel_center_y) > ppt.getPitchSimY() / 2.) {
0128 thePixelChargeMap[pixel_no] += charge_in_pixel;
0129 continue;
0130 }
0131 }
0132 }
0133 if (pixel_row == ppt.getNoPixelsSimX() - 1) {
0134 pixel_width_x = ppt.getPitchSimX();
0135 pixel_center_x = pixel_lower_x + ppt.getPitchSimX() / 2.;
0136 if (fabs((*i).Position().x() - pixel_center_x) > ppt.getPitchSimX() / 2.) {
0137 thePixelChargeMap[pixel_no] += charge_in_pixel;
0138 continue;
0139 }
0140 if (pixel_col == 0) {
0141 pixel_width_y = ppt.getPitchSimY();
0142 pixel_center_y = pixel_lower_y + ppt.getPhysActiveEdgeDist() + ppt.getPitchSimY() / 2.;
0143 if (fabs((*i).Position().y() - pixel_center_y) > ppt.getPitchSimY() / 2.) {
0144 thePixelChargeMap[pixel_no] += charge_in_pixel;
0145 continue;
0146 }
0147 }
0148 if (pixel_col == ppt.getNoPixelsSimY() - 1) {
0149 pixel_width_y = ppt.getPitchSimY();
0150 pixel_center_y = pixel_lower_y + ppt.getPitchSimY() / 2.;
0151 if (fabs((*i).Position().y() - pixel_center_y) > ppt.getPitchSimY() / 2.) {
0152 thePixelChargeMap[pixel_no] += charge_in_pixel;
0153 continue;
0154 }
0155 }
0156 }
0157 if (pixel_col == 0) {
0158 pixel_width_y = ppt.getPitchSimY();
0159 pixel_center_y = pixel_lower_y + ppt.getPhysActiveEdgeDist() + ppt.getPitchSimY() / 2.;
0160 if (fabs((*i).Position().y() - pixel_center_y) > ppt.getPitchSimY() / 2.) {
0161 thePixelChargeMap[pixel_no] += charge_in_pixel;
0162 continue;
0163 }
0164 }
0165 if (pixel_col == ppt.getNoPixelsSimY() - 1) {
0166 pixel_width_y = ppt.getPitchSimY();
0167 pixel_center_y = pixel_lower_y + ppt.getPitchSimY() / 2.;
0168 if (fabs((*i).Position().y() - pixel_center_y) > ppt.getPitchSimY() / 2.) {
0169 thePixelChargeMap[pixel_no] += charge_in_pixel;
0170 continue;
0171 }
0172 }
0173
0174
0175 int xbin = int((((*i).Position().y() - pixel_center_y) + pixel_width_y / 2.) * 1.e3 / 5.);
0176 int ybin = int((((*i).Position().x() - pixel_center_x) + pixel_width_x / 2.) * 1.e3 / 5.);
0177 if (pixel_width_x < 0.11 && pixel_width_y < 0.151) {
0178 psize = 0;
0179 if (xbin > xBinMax_[psize] || ybin > yBinMax_[psize]) {
0180 edm::LogError("PPS RPixChargeShare") << " Array index out of bounds 0";
0181 continue;
0182 }
0183 }
0184 if (pixel_width_x > 0.11 && pixel_width_y < 0.151) {
0185 psize = 2;
0186 if (xbin > xBinMax_[psize] || ybin > yBinMax_[psize]) {
0187 edm::LogError("PPS RPixChargeShare") << " Array index out of bounds 2";
0188 continue;
0189 }
0190 }
0191 if (pixel_width_x < 0.11 && pixel_width_y > 0.151) {
0192 psize = 1;
0193 if (xbin > xBinMax_[psize] || ybin > yBinMax_[psize]) {
0194 edm::LogError("PPS RPixChargeShare") << " Array index out of bounds 1";
0195 continue;
0196 }
0197 }
0198 if (pixel_width_x > 0.11 && pixel_width_y > 0.151) {
0199 psize = 3;
0200 if (xbin > xBinMax_[psize] || ybin > yBinMax_[psize]) {
0201 edm::LogError("PPS RPixChargeShare") << " Array index out of bounds 3";
0202 continue;
0203 }
0204 }
0205 if (xbin < 0 || ybin < 0) {
0206 edm::LogError("PPS RPixChargeShare") << " Negative array index xbin or ybin";
0207 continue;
0208 }
0209 double hit2neighbour[8];
0210 double collect_prob = chargeMap2E_[psize][xbin][ybin];
0211 thePixelChargeMap[pixel_no] += charge_in_pixel * collect_prob;
0212 unsigned short neighbour_no[8];
0213 unsigned short m = 0;
0214 double closer_neighbour = 0;
0215 unsigned short closer_no = 0;
0216
0217 for (int k = pixel_row - 1; k <= pixel_row + 1; k++) {
0218 for (int l = pixel_col - 1; l <= pixel_col + 1; l++) {
0219 if ((k < 0) || k > pxlRowSize_ - 1 || l < 0 || l > pxlColSize_ - 1)
0220 continue;
0221 if ((k == pixel_row) && (l == pixel_col))
0222 continue;
0223 double neighbour_pixel_lower_x = 0;
0224 double neighbour_pixel_lower_y = 0;
0225 double neighbour_pixel_upper_x = 0;
0226 double neighbour_pixel_upper_y = 0;
0227 double neighbour_pixel_center_x = 0;
0228 double neighbour_pixel_center_y = 0;
0229
0230 ppt.pixelRange(
0231 k, l, neighbour_pixel_lower_x, neighbour_pixel_upper_x, neighbour_pixel_lower_y, neighbour_pixel_upper_y);
0232 neighbour_pixel_center_x = neighbour_pixel_lower_x + (neighbour_pixel_upper_x - neighbour_pixel_lower_x) / 2.;
0233 neighbour_pixel_center_y = neighbour_pixel_lower_y + (neighbour_pixel_upper_y - neighbour_pixel_lower_y) / 2.;
0234 hit2neighbour[m] = sqrt(pow((*i).Position().x() - neighbour_pixel_center_x, 2.) +
0235 pow((*i).Position().y() - neighbour_pixel_center_y, -2.));
0236 neighbour_no[m] = l * pxlRowSize_ + k;
0237 if (hit2neighbour[m] > closer_neighbour) {
0238 closer_neighbour = hit2neighbour[m];
0239 closer_no = neighbour_no[m];
0240 }
0241 m++;
0242 }
0243 }
0244 double chargetransfereff = (1 - collect_prob) * signalCoupling_[0];
0245 thePixelChargeMap[closer_no] += charge_in_pixel * chargetransfereff;
0246 }
0247 }
0248
0249 return thePixelChargeMap;
0250 }