File indexing completed on 2023-10-25 10:05:00
0001 #include <cmath>
0002 #include <vector>
0003 #include <algorithm>
0004
0005 #include "CalibTracker/SiPixelESProducers/interface/SiPixelGainCalibrationOfflineSimService.h"
0006 #include "FWCore/Framework/interface/ConsumesCollector.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0010 #include "SimTracker/SiPhase2Digitizer/plugins/Pixel3DDigitizerAlgorithm.h"
0011
0012 using namespace sipixelobjects;
0013
0014 namespace {
0015
0016 constexpr double operator""_um(long double length) { return length * 1e-4; }
0017 constexpr double operator""_um_inv(long double length) { return length * 1e4; }
0018 }
0019
0020 Pixel3DDigitizerAlgorithm::Pixel3DDigitizerAlgorithm(const edm::ParameterSet& conf, edm::ConsumesCollector iC)
0021 : PixelDigitizerAlgorithm(conf, iC),
0022 np_column_radius_(
0023 (conf.getParameter<edm::ParameterSet>("Pixel3DDigitizerAlgorithm").getParameter<double>("NPColumnRadius")) *
0024 1.0_um),
0025 ohm_column_radius_(
0026 (conf.getParameter<edm::ParameterSet>("Pixel3DDigitizerAlgorithm").getParameter<double>("OhmicColumnRadius")) *
0027 1.0_um),
0028 np_column_gap_(
0029 (conf.getParameter<edm::ParameterSet>("Pixel3DDigitizerAlgorithm").getParameter<double>("NPColumnGap")) *
0030 1.0_um) {
0031
0032 pixelFlag_ = true;
0033
0034 edm::LogInfo("Pixel3DDigitizerAlgorithm")
0035 << "Algorithm constructed \n"
0036 << "Configuration parameters:\n"
0037 << "\n*** Threshold"
0038 << "\n Endcap = " << theThresholdInE_Endcap_ << " electrons"
0039 << "\n Barrel = " << theThresholdInE_Barrel_ << " electrons"
0040 << "\n*** Gain"
0041 << "\n Electrons per ADC:" << theElectronPerADC_ << "\n ADC Full Scale: " << theAdcFullScale_
0042 << "\n*** The delta cut-off is set to " << tMax_ << "\n*** Pixel-inefficiency: " << addPixelInefficiency_;
0043 }
0044
0045 Pixel3DDigitizerAlgorithm::~Pixel3DDigitizerAlgorithm() {}
0046
0047 const bool Pixel3DDigitizerAlgorithm::is_inside_n_column_(const LocalPoint& p, const float& sensor_thickness) const {
0048
0049 return (p.perp() <= np_column_radius_ && p.z() <= (sensor_thickness - np_column_gap_));
0050 }
0051
0052 const bool Pixel3DDigitizerAlgorithm::is_inside_ohmic_column_(const LocalPoint& p,
0053 const std::pair<float, float>& half_pitch) const {
0054
0055 return ((p - LocalVector(half_pitch.first, half_pitch.second, 0)).perp() <= ohm_column_radius_) ||
0056 ((p - LocalVector(-half_pitch.first, half_pitch.second, 0)).perp() <= ohm_column_radius_) ||
0057 ((p - LocalVector(half_pitch.first, -half_pitch.second, 0)).perp() <= ohm_column_radius_) ||
0058 ((p - LocalVector(-half_pitch.first, -half_pitch.second, 0)).perp() <= ohm_column_radius_);
0059 }
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070 std::vector<digitizerUtility::EnergyDepositUnit> Pixel3DDigitizerAlgorithm::diffusion(
0071 const LocalPoint& pos,
0072 const float& ncarriers,
0073 const std::function<LocalVector(float, float)>& u_drift,
0074 const std::pair<float, float> hpitches,
0075 const float& thickness) const {
0076
0077
0078
0079
0080 const float max_migration_radius = 0.4_um;
0081
0082 int displ_ind = -1;
0083 float pitch = 0.0;
0084
0085
0086
0087 if (hpitches.first - std::abs(pos.x()) < max_migration_radius) {
0088 displ_ind = 0;
0089 pitch = hpitches.first;
0090 } else if (hpitches.second - std::abs(pos.y()) < max_migration_radius) {
0091 displ_ind = 1;
0092 pitch = hpitches.second;
0093 } else {
0094
0095 return std::vector<digitizerUtility::EnergyDepositUnit>();
0096 }
0097
0098
0099
0100 std::vector<digitizerUtility::EnergyDepositUnit> migrated_charge;
0101
0102
0103 const float diffusion_step = 0.1_um;
0104
0105
0106 std::vector<float> pos_moving({pos.x(), pos.y(), pos.z()});
0107
0108 std::function<std::vector<float>(int)> do_step =
0109 [&pos_moving, &u_drift, diffusion_step](int i) -> std::vector<float> {
0110 auto dd = u_drift(pos_moving[0], pos_moving[1]);
0111 return std::vector<float>({i * diffusion_step * dd.x(), i * diffusion_step * dd.y(), i * diffusion_step * dd.z()});
0112 };
0113
0114 LogDebug("Pixel3DDigitizerAlgorithm::diffusion")
0115 << "\nMax. radius from the pixel edge to migrate charge: " << max_migration_radius * 1.0_um_inv << " [um]"
0116 << "\nMigration axis: " << displ_ind
0117 << "\n(super-)Charge distance to the pixel edge: " << (pitch - pos_moving[displ_ind]) * 1.0_um_inv << " [um]";
0118
0119
0120 const float N_SIGMA = 3.0;
0121
0122
0123
0124 float current_carriers = ncarriers;
0125 std::vector<float> newpos({pos_moving[0], pos_moving[1], pos_moving[2]});
0126 float distance_edge = 0.0_um;
0127
0128 const float sigma = 0.4_um;
0129 for (int i = 1;; ++i) {
0130 std::transform(pos_moving.begin(), pos_moving.end(), do_step(i).begin(), pos_moving.begin(), std::plus<float>());
0131 distance_edge = pitch - std::abs(pos_moving[displ_ind]);
0132
0133
0134 float migrated_e = current_carriers * 0.5 * (1.0 - std::erf(distance_edge / (sigma * std::sqrt(2.0))));
0135
0136 LogDebug("(super-)charge diffusion") << "step-" << i << ", Current carriers Ne= " << current_carriers << ","
0137 << "r=(" << pos_moving[0] * 1.0_um_inv << ", " << pos_moving[1] * 1.0_um_inv
0138 << ", " << pos_moving[2] * 1.0_um_inv << ") [um], "
0139 << "Migrated charge: " << migrated_e;
0140
0141
0142 current_carriers -= migrated_e;
0143
0144
0145 if (std::abs(distance_edge) >= max_migration_radius || current_carriers <= 0.5 * ncarriers) {
0146 break;
0147 }
0148
0149
0150
0151
0152
0153 std::vector<float> newpos(pos_moving);
0154
0155 newpos[displ_ind] += std::copysign(N_SIGMA * sigma, newpos[displ_ind]);
0156 migrated_charge.push_back(digitizerUtility::EnergyDepositUnit(migrated_e, newpos[0], newpos[1], newpos[2]));
0157 }
0158 return migrated_charge;
0159 }
0160
0161
0162
0163
0164
0165
0166
0167 std::vector<digitizerUtility::SignalPoint> Pixel3DDigitizerAlgorithm::drift(
0168 const PSimHit& hit,
0169 const Phase2TrackerGeomDetUnit* pixdet,
0170 const GlobalVector& bfield,
0171 const std::vector<digitizerUtility::EnergyDepositUnit>& ionization_points) const {
0172 return driftFor3DSensors(hit, pixdet, bfield, ionization_points, true);
0173 }
0174 std::vector<digitizerUtility::SignalPoint> Pixel3DDigitizerAlgorithm::driftFor3DSensors(
0175 const PSimHit& hit,
0176 const Phase2TrackerGeomDetUnit* pixdet,
0177 const GlobalVector& bfield,
0178 const std::vector<digitizerUtility::EnergyDepositUnit>& ionization_points,
0179 bool diffusion_activated) const {
0180
0181
0182
0183
0184
0185
0186
0187
0188 const auto pitch = pixdet->specificTopology().pitch();
0189 const auto half_pitch = std::make_pair<float, float>(pitch.first * 0.5, pitch.second * 0.5);
0190 const float thickness = pixdet->specificSurface().bounds().thickness();
0191 const int nrows = pixdet->specificTopology().nrows();
0192 const int ncolumns = pixdet->specificTopology().ncolumns();
0193 const float pix_rounding = 0.99;
0194
0195
0196 const float max_radial_distance =
0197 std::sqrt(half_pitch.first * half_pitch.first + half_pitch.second * half_pitch.second);
0198
0199
0200
0201
0202
0203
0204 LocalPoint center_proxy_cell(half_pitch.first, half_pitch.second, -0.5 * thickness);
0205
0206 LogDebug("Pixel3DDigitizerAlgorithm::drift")
0207 << "Pixel pitch:" << pitch.first * 1.0_um_inv << ", " << pitch.second * 1.0_um_inv << " [um]";
0208
0209
0210
0211 std::function<LocalVector(float, float)> drift_direction = [](float x, float y) -> LocalVector {
0212 const float theta = std::atan2(y, x);
0213 return LocalVector(-std::cos(theta), -std::sin(theta), 0.0);
0214 };
0215
0216 std::vector<digitizerUtility::SignalPoint> collection_points;
0217
0218 collection_points.reserve(ionization_points.size());
0219
0220
0221
0222 const float RAD_DAMAGE = 0.001;
0223
0224 for (const auto& super_charge : ionization_points) {
0225
0226 auto current_pixel = pixdet->specificTopology().pixel(LocalPoint(super_charge.x(), super_charge.y()));
0227
0228
0229
0230
0231
0232
0233 current_pixel.first = std::clamp(current_pixel.first, float(0.0), (nrows - 1) + pix_rounding);
0234 current_pixel.second = std::clamp(current_pixel.second, float(0.0), (ncolumns - 1) + pix_rounding);
0235
0236 const auto current_pixel_int = std::make_pair(std::floor(current_pixel.first), std::floor(current_pixel.second));
0237
0238
0239
0240 const auto relative_position_at_pc =
0241 std::make_pair((current_pixel.first - current_pixel_int.first) * pitch.first,
0242 (current_pixel.second - current_pixel_int.second) * pitch.second);
0243
0244
0245 LocalPoint position_at_pc(relative_position_at_pc.first - center_proxy_cell.x(),
0246 relative_position_at_pc.second - center_proxy_cell.y(),
0247 super_charge.z() - center_proxy_cell.z());
0248
0249 LogDebug("Pixel3DDigitizerAlgorithm::drift")
0250 << "(super-)Charge\nlocal position: (" << super_charge.x() * 1.0_um_inv << ", " << super_charge.y() * 1.0_um_inv
0251 << ", " << super_charge.z() * 1.0_um_inv << ") [um]"
0252 << "\nMeasurement Point (row,column) (" << current_pixel.first << ", " << current_pixel.second << ")"
0253 << "\nProxy pixel-cell frame (centered at left-back corner): (" << relative_position_at_pc.first * 1.0_um_inv
0254 << ", " << relative_position_at_pc.second * 1.0_um_inv << ") [um]"
0255 << "\nProxy pixel-cell frame (centered at n-column): (" << position_at_pc.x() * 1.0_um_inv << ", "
0256 << position_at_pc.y() * 1.0_um_inv << ") [um] "
0257 << "\nNe=" << super_charge.energy() << " electrons";
0258
0259
0260 if (is_inside_n_column_(position_at_pc, thickness) || is_inside_ohmic_column_(position_at_pc, half_pitch)) {
0261 LogDebug("Pixel3DDigitizerAlgorithm::drift") << "Remove charge, inside the n-column or p-column!!";
0262 continue;
0263 }
0264
0265 float nelectrons = super_charge.energy();
0266
0267 if (diffusion_activated) {
0268 auto migrated_charges = diffusion(position_at_pc, super_charge.energy(), drift_direction, half_pitch, thickness);
0269 for (auto& mc : migrated_charges) {
0270
0271 nelectrons -= mc.energy();
0272
0273
0274
0275 const float pixel_x = current_pixel_int.first + (mc.x() + center_proxy_cell.x()) / pitch.first;
0276 const float pixel_y = current_pixel_int.second + (mc.y() + center_proxy_cell.y()) / pitch.second;
0277 const auto lp = pixdet->specificTopology().localPosition(MeasurementPoint(pixel_x, pixel_y));
0278
0279
0280
0281
0282 mc.migrate_position(LocalPoint(lp.x(), lp.y(), mc.z() + center_proxy_cell.z()));
0283 }
0284 if (!migrated_charges.empty()) {
0285 LogDebug("Pixel3DDigitizerAlgorithm::drift") << "****************"
0286 << "MIGRATING (super-)charges"
0287 << "****************";
0288
0289 auto mig_colpoints = driftFor3DSensors(hit, pixdet, bfield, migrated_charges, false);
0290 collection_points.insert(std::end(collection_points), mig_colpoints.begin(), mig_colpoints.end());
0291 LogDebug("Pixel3DDigitizerAlgorithm::drift") << "*****************"
0292 << "DOME MIGRATION"
0293 << "****************";
0294 }
0295 }
0296
0297
0298
0299
0300 const float drift_distance = position_at_pc.perp() - np_column_radius_;
0301
0302
0303
0304 float energyOnCollector = nelectrons;
0305
0306
0307 if (pseudoRadDamage_ >= RAD_DAMAGE) {
0308 const float module_radius = pixdet->surface().position().perp();
0309 if (module_radius <= pseudoRadDamageRadius_) {
0310 const float kValue = pseudoRadDamage_ / (module_radius * module_radius);
0311 energyOnCollector = energyOnCollector * std::exp(-1.0 * kValue * drift_distance / max_radial_distance);
0312 }
0313 }
0314 LogDebug("Pixel3DDigitizerAlgorithm::drift")
0315 << "Drift distance = " << drift_distance * 1.0_um_inv << " [um], "
0316 << "Initial electrons = " << super_charge.energy()
0317 << " [electrons], Electrons after loss/diff= " << energyOnCollector << " [electrons] ";
0318
0319
0320 collection_points.push_back(digitizerUtility::SignalPoint(
0321 current_pixel_int.first, current_pixel_int.second, 0.0, 0.0, hit.tof(), energyOnCollector));
0322 }
0323
0324 return collection_points;
0325 }
0326
0327
0328
0329
0330
0331 void Pixel3DDigitizerAlgorithm::induce_signal(std::vector<PSimHit>::const_iterator inputBegin,
0332 const PSimHit& hit,
0333 const size_t hitIndex,
0334 const size_t firstHitIndex,
0335 const uint32_t tofBin,
0336 const Phase2TrackerGeomDetUnit* pixdet,
0337 const std::vector<digitizerUtility::SignalPoint>& collection_points) {
0338
0339
0340 const uint32_t detId = pixdet->geographicalId().rawId();
0341
0342 signal_map_type& the_signal = _signal[detId];
0343
0344
0345 std::function<int(int, int)> pixelToChannel =
0346 pixelFlag_ ? PixelDigi::pixelToChannel
0347 : static_cast<std::function<int(int, int)> >(Phase2TrackerDigi::pixelToChannel);
0348
0349
0350 for (const auto& pt : collection_points) {
0351
0352 const int channel = pixelToChannel(pt.position().x(), pt.position().y());
0353
0354 float corr_time = hit.tof() - pixdet->surface().toGlobal(hit.localPosition()).mag() * c_inv;
0355 if (makeDigiSimLinks_) {
0356 the_signal[channel] +=
0357 digitizerUtility::Ph2Amplitude(pt.amplitude(), &hit, pt.amplitude(), corr_time, hitIndex, tofBin);
0358 } else {
0359 the_signal[channel] += digitizerUtility::Ph2Amplitude(pt.amplitude(), nullptr, pt.amplitude());
0360 }
0361
0362 LogDebug("Pixel3DDigitizerAlgorithm")
0363 << " Induce charge at row,col:" << pt.position() << " N_electrons:" << pt.amplitude() << " [Channel:" << channel
0364 << "]\n [Accumulated signal in this channel:" << the_signal[channel].ampl() << "] "
0365 << " Global index linked PSimHit:" << hitIndex;
0366 }
0367 }