RPixChargeShare.cc

Go to the documentation of this file.

#include "SimPPS/PPSPixelDigiProducer/interface/RPixChargeShare.h"
#include "FWCore/ParameterSet/interface/FileInPath.h"
#include <iostream>
#include <fstream>
 
RPixChargeShare::RPixChargeShare(const edm::ParameterSet &params, uint32_t det_id, const PPSPixelTopology &ppt)
    : det_id_(det_id) {
  verbosity_ = params.getParameter<int>("RPixVerbosity");
  signalCoupling_.clear();
  ChargeMapFile2E_[0] = params.getParameter<std::string>("ChargeMapFile2E");
  ChargeMapFile2E_[1] = params.getParameter<std::string>("ChargeMapFile2E_2X");
  ChargeMapFile2E_[2] = params.getParameter<std::string>("ChargeMapFile2E_2Y");
  ChargeMapFile2E_[3] = params.getParameter<std::string>("ChargeMapFile2E_2X2Y");
 
  double coupling_constant_ = params.getParameter<double>("RPixCoupling");
 
  signalCoupling_.push_back(coupling_constant_);
  signalCoupling_.push_back((1.0 - coupling_constant_) / 2);
 
  no_of_pixels_ = ppt.getNoPixels();
 
  double xMap, yMap;
  double chargeprobcollect;
  int xUpper[] = {75, 150, 75, 150};
  int yUpper[] = {50, 50, 100, 100};
  int ix, iy;
  for (int i = 0; i < 4; i++) {
    edm::FileInPath filename(ChargeMapFile2E_[i]);
    std::ifstream fChargeMap(filename.fullPath().c_str());
    if (fChargeMap.is_open()) {
      while (fChargeMap >> xMap >> yMap >> chargeprobcollect) {
        ix = int((xMap + xUpper[i]) / 5);
        iy = int((yMap + yUpper[i]) / 5);
        chargeMap2E_[i][ix][iy] = chargeprobcollect;
      }
      fChargeMap.close();
    } else
      throw cms::Exception("PPS RPixChargeShare") << "Charge map file not found";
  }
}
 
std::map<unsigned short, double> RPixChargeShare::Share(const std::vector<RPixSignalPoint> &charge_map,
                                                        const PPSPixelTopology &ppt) {
  std::map<unsigned short, double> thePixelChargeMap;
  if (verbosity_ > 1)
    edm::LogInfo("PPS") << "RPixChargeShare " << det_id_ << " : Clouds to be induced= " << charge_map.size();
 
  double cH = 0;
 
  for (std::vector<RPixSignalPoint>::const_iterator i = charge_map.begin(); i != charge_map.end(); ++i) {
    double hit_pos_x, hit_pos_y;
    // Used to avoid the abort due to hits out of detector
    if (((*i).Position().x() + ppt.getSimXWidth() / 2.) < 0 ||
        ((*i).Position().x() + ppt.getSimXWidth() / 2.) > ppt.getSimXWidth()) {
      edm::LogInfo("PPS RPixChargeShare")
          << "**** Attention ((*i).Position().x()+simX_width_/2.)<0||((*i).Position().x()+simX_width_/2.)>simX_width  ";
      edm::LogInfo("PPS RPixChargeShare") << "(*i).Position().x() = " << (*i).Position().x();
      continue;
    }
    if (((*i).Position().y() + ppt.getSimYWidth() / 2.) < 0 ||
        ((*i).Position().y() + ppt.getSimYWidth() / 2.) > ppt.getSimYWidth()) {
      edm::LogInfo("PPS RPixChargeShare")
          << "**** Attention ((*i).Position().y()+simY_width_/2.)<0||((*i).Position().y()+simY_width_/2.)>simY_width  ";
      edm::LogInfo("PPS RPixChargeShare") << "(*i).Position().y() = " << (*i).Position().y();
      continue;
    }
 
    PPSPixelTopology::PixelInfo relevant_pixels =
        ppt.getPixelsInvolved((*i).Position().x(), (*i).Position().y(), (*i).Sigma(), hit_pos_x, hit_pos_y);
    double effic = relevant_pixels.effFactor();
    if (effic < 1.e-4)
      continue;
 
    unsigned short pixel_no = ppt.pixelIndex(relevant_pixels);
 
    double charge_in_pixel = (*i).Charge() * effic;
 
    cH += charge_in_pixel;
 
    if (verbosity_ > 1)
      edm::LogInfo("PPS RPixChargeShare ") << "Efficiency in detector " << det_id_ << " and pixel no " << pixel_no
                                           << "  : " << effic << "  ch: " << charge_in_pixel << "   CHtot: " << cH;
 
    if (signalCoupling_[0] == 0.) {
      thePixelChargeMap[pixel_no] += charge_in_pixel;
    } else {
      int pixel_row = relevant_pixels.pixelRowNo();
      int pixel_col = relevant_pixels.pixelColNo();
      double pixel_lower_x = 0;
      double pixel_lower_y = 0;
      double pixel_upper_x = 0;
      double pixel_upper_y = 0;
      int psize = 0;
      ppt.pixelRange(pixel_row, pixel_col, pixel_lower_x, pixel_upper_x, pixel_lower_y, pixel_upper_y);
      double pixel_width_x = pixel_upper_x - pixel_lower_x;
      double pixel_width_y = pixel_upper_y - pixel_lower_y;
      double pixel_center_x = pixel_lower_x + (pixel_width_x) / 2.;
      double pixel_center_y = pixel_lower_y + (pixel_width_y) / 2.;
      // xbin and ybin are coordinates (um) inside the pixel as in the test beam, swapped wrt plane coordinates.
      int xbin = int((((*i).Position().y() - pixel_center_y) + pixel_width_y / 2.) * 1.e3 / 5.);
      int ybin = int((((*i).Position().x() - pixel_center_x) + pixel_width_x / 2.) * 1.e3 / 5.);
      if (pixel_width_x < 0.11 && pixel_width_y < 0.151) {  // pixel 100x150 um^2
        psize = 0;
        if (xbin > xBinMax_[psize] || ybin > yBinMax_[psize]) {
          edm::LogError("PPS RPixChargeShare") << " Array index out of bounds";
          continue;
        }
      }
      if (pixel_width_x > 0.11 && pixel_width_y < 0.151) {  // pixel 200x150 um^2
        psize = 2;
        if (xbin > xBinMax_[psize] || ybin > yBinMax_[psize]) {
          edm::LogError("PPS RPixChargeShare") << " Array index out of bounds";
          continue;
        }
      }
      if (pixel_width_x < 0.11 && pixel_width_y > 0.151) {  // pixel 100x300 um^2
        psize = 1;
        if (xbin > xBinMax_[psize] || ybin > yBinMax_[psize]) {
          edm::LogError("PPS RPixChargeShare") << " Array index out of bounds";
          continue;
        }
      }
      if (pixel_width_x > 0.11 && pixel_width_y > 0.151) {  // pixel 200x300 um^2
        psize = 3;
        if (xbin > xBinMax_[psize] || ybin > yBinMax_[psize]) {
          edm::LogError("PPS RPixChargeShare") << " Array index out of bounds";
          continue;
        }
      }
      if (xbin < 0 || ybin < 0) {
        edm::LogError("PPS RPixChargeShare") << " Negative array index xbin or ybin";
        continue;
      }
      double hit2neighbour[8];
      double collect_prob = chargeMap2E_[psize][xbin][ybin];
      thePixelChargeMap[pixel_no] += charge_in_pixel * collect_prob;
      unsigned short neighbour_no[8];
      unsigned short m = 0;
      double closer_neighbour = 0;
      unsigned short closer_no = 0;
      //      Considering the 8 neighbours to share charge
      for (int k = pixel_row - 1; k <= pixel_row + 1; k++) {
        for (int l = pixel_col - 1; l <= pixel_col + 1; l++) {
          if ((k < 0) || k > pxlRowSize_ - 1 || l < 0 || l > pxlColSize_ - 1)
            continue;
          if ((k == pixel_row) && (l == pixel_col))
            continue;
          double neighbour_pixel_lower_x = 0;
          double neighbour_pixel_lower_y = 0;
          double neighbour_pixel_upper_x = 0;
          double neighbour_pixel_upper_y = 0;
          double neighbour_pixel_center_x = 0;
          double neighbour_pixel_center_y = 0;
          //      Check the hit approach to the neighbours
          ppt.pixelRange(
              k, l, neighbour_pixel_lower_x, neighbour_pixel_upper_x, neighbour_pixel_lower_y, neighbour_pixel_upper_y);
          neighbour_pixel_center_x = neighbour_pixel_lower_x + (neighbour_pixel_upper_x - neighbour_pixel_lower_x) / 2.;
          neighbour_pixel_center_y = neighbour_pixel_lower_y + (neighbour_pixel_upper_y - neighbour_pixel_lower_y) / 2.;
          hit2neighbour[m] = sqrt(pow((*i).Position().x() - neighbour_pixel_center_x, 2.) +
                                  pow((*i).Position().y() - neighbour_pixel_center_y, -2.));
          neighbour_no[m] = l * pxlRowSize_ + k;
          if (hit2neighbour[m] > closer_neighbour) {
            closer_neighbour = hit2neighbour[m];
            closer_no = neighbour_no[m];
          }
          m++;
        }
      }
      double chargetransfereff = (1 - collect_prob) * signalCoupling_[0];
      thePixelChargeMap[closer_no] += charge_in_pixel * chargetransfereff;
    }
  }
 
  return thePixelChargeMap;
}