Back to home page

Project CMSSW displayed by LXR

 
 

    


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 &params, 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     // Used to avoid the abort due to hits out of detector
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   First/last rows/columns are bigger due to sensor edges. 0.2mm is added to the pixel pitch but 
0102   only 0.15mm is considered sensitive. 
0103 The pixel size is standard (100x150) but the actual one is +200 and the sensitive +150.
0104   Row 0: 100x150 -> 250x150 sensitive (lower x 50, higher x 300) but charge share only in 100x150 and
0105   x center in lower_x + active_edge + pitch_x/2.
0106   Row 159: 100x150 -> 250x150 sensitive (lower x 16300, higher x 16550) but charge share only in 100x150 and
0107   x center in lower_x + pitch_x/2.
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       // xbin and ybin are coordinates (um) inside the pixel as in the test beam, swapped wrt plane coordinates.
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) {  // pixel 100x150 um^2
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) {  // pixel 200x150 um^2
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) {  // pixel 100x300 um^2
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) {  // pixel 200x300 um^2
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       //      Considering the 8 neighbours to share charge
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           //      Check the hit approach to the neighbours
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 }