ATLAS Offline Software
SensorSim3DTool.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3 */
4 #include "SensorSim3DTool.h"
7 #include "InDetSimEvent/SiHit.h"
13 
14 #include "GaudiKernel/SystemOfUnits.h"
15 
16 #include "CLHEP/Random/RandFlat.h"
17 #include "CLHEP/Random/RandGaussZiggurat.h"
18 
19 #include "CLHEP/Units/PhysicalConstants.h"
21 
22 #include "TFile.h"
23 
24 #include <cmath>
25 #include <memory>
26 
27 using namespace InDetDD;
28 
29 //===============================================
30 // C O N S T R U C T O R
31 //===============================================
32 SensorSim3DTool::SensorSim3DTool(const std::string& type, const std::string& name, const IInterface* parent) :
34 }
35 
37 
38 //===============================================
39 // I N I T I A L I Z E
40 //===============================================
42  ATH_MSG_DEBUG("SensorSim3DTool::initialize()");
43 
45  ATH_CHECK(m_radDamageUtil.retrieve());
46 
49 
50  // read the template correction maps if template correction is used
51  if (m_radiationDamageSimulationType == RadiationDamageSimulationType::TEMPLATE_CORRECTION) {
52  ATH_MSG_INFO("Opening file: " << m_templateCorrectionRootFile << " for radiation damage correction");
53  std::unique_ptr<TFile> file(TFile::Open(PathResolverFindCalibFile(m_templateCorrectionRootFile).c_str(), "READ"));
54  if (!file) {
55  ATH_MSG_ERROR("Unable to read the ROOT file needed for radiation damage correction at: " << m_templateCorrectionRootFile);
56  return StatusCode::FAILURE;
57  }
58 
59  // we need separate corrections for the barrel and for the rings
60  constexpr std::size_t numberOfLayers = 2;
61  if (m_chargeCorrectionHistos.size() != numberOfLayers) {
62  ATH_MSG_ERROR("The size of the vector of charge correction histogram paths is " << m_chargeCorrectionHistos.size() << " instead of " << numberOfLayers);
63  return StatusCode::FAILURE;
64  }
65 
66  for (const std::string& i : m_chargeCorrectionHistos) {
67  ATH_MSG_INFO("Reading histogram: " << i << " for charge correction needed for Pixel radiation damage simulation");
68  std::unique_ptr<TH1> h(file->Get<TH1>(i.c_str()));
69  if (!h) {
70  ATH_MSG_ERROR("Cannot read histogram " << i << " needed for charge correction");
71  return StatusCode::FAILURE;
72  }
73 
74  h->SetDirectory(nullptr);
75 
76  m_chargeCorrection.emplace_back();
77  ATH_CHECK(m_chargeCorrection.back().setHisto1D(h.get()));
78  }
79 
80  file->Close();
81  }
82 
83  return StatusCode::SUCCESS;
84 }
85 
86 //===============================================
87 // F I N A L I Z E
88 //===============================================
90  ATH_MSG_DEBUG("SensorSim3DTool::finalize()");
91  return StatusCode::SUCCESS;
92 }
93 
94 //===============================================
95 // I N D U C E C H A R G E
96 //===============================================
98  SiChargedDiodeCollection& chargedDiodes,
99  const InDetDD::SiDetectorElement& Module,
100  const InDetDD::PixelModuleDesign& p_design,
101  std::vector< std::pair<double, double> >& trfHitRecord,
102  std::vector<double>& initialConditions,
103  CLHEP::HepRandomEngine* rndmEngine,
104  const EventContext &ctx) const {
105 
107  if (m_digitizeITk3Das3D) {
108  if (!p_design.is3D()) {
109  return StatusCode::SUCCESS;
110  }
111 
112  // for now skip pixel luminosity rings
113  if (Module.isPLR()) {
114  return StatusCode::SUCCESS;
115  }
116  } else {
117  return StatusCode::SUCCESS;
118  }
119  } else {
120  if (!Module.isBarrel()) {
121  return StatusCode::SUCCESS;
122  }
124  return StatusCode::SUCCESS;
125  }
126  if (p_design.numberOfCircuits() > 1) {
127  return StatusCode::SUCCESS;
128  }
129  }
130 
131  ATH_MSG_DEBUG("Applying SensorSim3D charge processor");
132  if (initialConditions.size() != 8) {
133  ATH_MSG_ERROR("Starting coordinates were not filled correctly in EnergyDepositionSvc.");
134  return StatusCode::FAILURE;
135  }
136 
137  double eta_0 = initialConditions[0];
138  double phi_0 = initialConditions[1];
139  double depth_0 = initialConditions[2];
140  double dEta = initialConditions[3];
141  double dPhi = initialConditions[4];
142  double dDepth = initialConditions[5];
143  double ncharges = initialConditions[6];
144  double iTotalLength = initialConditions[7];
145  ncharges = 50;
146 
147  ATH_MSG_VERBOSE("Applying 3D sensor simulation.");
148  double sensorThickness = Module.design().thickness();
149  const InDet::SiliconProperties& siProperties = m_siPropertiesTool->getSiProperties(Module.identifyHash(), ctx);
150  double eleholePairEnergy = siProperties.electronHolePairsPerEnergy();
151 
152 
153  // Charge Collection Probability Map bin size
154  const double x_bin_size = 0.001;
155  const double y_bin_size = 0.001;
156 
157  // determine which readout is used
158  // FEI4 : 50 X 250 microns
159  double pixel_size_x = Module.width() / p_design.rows();
160  double pixel_size_y = Module.length() / p_design.columns();
161  double module_size_x = Module.width();
162  double module_size_y = Module.length();
163 
164  const HepMcParticleLink particleLink = HepMcParticleLink::getRedirectedLink(phit->particleLink(), phit.eventId(), ctx); // This link should now correctly resolve to the TruthEvent McEventCollection in the main StoreGateSvc.
165  const double pHitTime = hitTime(phit);
166 
167  const int layerIndex = Module.isBarrel() ? 0 : 1;
168 
169  if (m_radiationDamageSimulationType == RadiationDamageSimulationType::RAMO_POTENTIAL) {
170  const bool doChunkCorrection = m_doChunkCorrection.value();
171  //**************************************//
172  //*** Now diffuse charges to surface *** //
173  //**************************************//
174  //Calculate trapping times based on fluence (already includes check for fluence=0)
176  const PixelRadiationDamageFluenceMapData *fluenceData = *fluenceDataHandle;
177 
178  auto [trappingTimeElectrons, trappingTimeHoles] = m_radDamageUtil->getTrappingTimes(fluenceData->getFluenceLayer3D(0)); //0 = IBL
179 
180  const PixelHistoConverter& ramoPotentialMap = fluenceData->getRamoPotentialMap3D(0);
181  const PixelHistoConverter& eFieldMap = fluenceData->getEFieldMap3D(0);
182  const PixelHistoConverter& xPositionMap_e = fluenceData->getXPositionMap3D_e(0);
183  const PixelHistoConverter& xPositionMap_h = fluenceData->getXPositionMap3D_h(0);
184  const PixelHistoConverter& yPositionMap_e = fluenceData->getYPositionMap3D_e(0);
185  const PixelHistoConverter& yPositionMap_h = fluenceData->getYPositionMap3D_h(0);
186  const PixelHistoConverter& timeMap_e = fluenceData->getTimeMap3D_e(0);
187  const PixelHistoConverter& timeMap_h = fluenceData->getTimeMap3D_h(0);
188  const PixelHistoConverter& avgChargeMap_e = fluenceData->getAvgChargeMap3D_e();
189  const PixelHistoConverter& avgChargeMap_h = fluenceData->getAvgChargeMap3D_h();
190 
191  //Parameters which will be smeared by a random value for each charge propagated,
192  //recomputed for every trfHitRecord but initialized only once
193  std::vector<double> DtElectron (ncharges, 0.);
194  std::vector<double> DtHole (ncharges, 0.);
195  std::vector<double> rdifElectron (ncharges, 0.);
196  std::vector<double> rdifHole (ncharges, 0.);
197 
198  for (const auto & iHitRecord : trfHitRecord) {
199  double eta_i = eta_0;
200  double phi_i = phi_0;
201  double depth_i = depth_0;
202 
203  if (iTotalLength) {
204  eta_i += 1.0 * iHitRecord.first / iTotalLength * dEta;
205  phi_i += 1.0 * iHitRecord.first / iTotalLength * dPhi;
206  depth_i += 1.0 * iHitRecord.first / iTotalLength * dDepth;
207  }
208 
209  const double energy_per_step = 1.0 * iHitRecord.second / 1.E+6 / ncharges; //in MeV
210  ATH_MSG_DEBUG("es_current: " << energy_per_step << " split between " << ncharges << " charges");
211 
212  double dist_electrode = 0.5 * sensorThickness - Module.design().readoutSide() * depth_i;
213  if (dist_electrode < 0) dist_electrode = 0;
214 
215  CLHEP::Hep3Vector chargepos;
216  chargepos.setX(phi_i);
217  chargepos.setY(eta_i);
218  chargepos.setZ(dist_electrode);
219 
220  bool coord = Module.isModuleFrame();
221 
223  "ismoduleframe " << coord << " -- startPosition (x,y,z) = " << chargepos.x() << ", " << chargepos.y() << ", " <<
224  chargepos.z());
225 
226  // -- change origin of coordinates to the left bottom of module
227  double x_new = chargepos.x() + 0.5*module_size_x;
228  double y_new = chargepos.y() + 0.5*module_size_y;
229 
230  // -- change from module frame to pixel frame
231  int nPixX = int(x_new / pixel_size_x);
232  int nPixY = int(y_new / pixel_size_y);
233  ATH_MSG_DEBUG(" -- nPixX = " << nPixX << " nPixY = " << nPixY);
234 
235  //position relative to the bottom left corner of the pixel
236  double x_pix = x_new - pixel_size_x * (nPixX);
237  double y_pix = y_new - pixel_size_y * (nPixY);
238  // -- change origin of coordinates to the center of the pixel
239  double x_pix_center = x_pix - pixel_size_x / 2;
240  double y_pix_center = y_pix - pixel_size_y / 2;
241  ATH_MSG_DEBUG(" -- current hit position w.r.t. pixel corner = " << x_pix << " " << y_pix);
242  ATH_MSG_DEBUG(" -- current hit position w.r.t. pixel center = " << x_pix_center << " " << y_pix_center);
243 
244  //only process hits which are not on the electrodes (E-field zero)
245  //all the maps have 250 as the x value, so need to invert x and y whenever reading maps
246  double efield = eFieldMap.getContent(eFieldMap.getBinX(1e3*y_pix),eFieldMap.getBinY(1e3*x_pix))*1.0E-7; //return efield in MV/mm (for mobility calculation)
247 
248  if (efield == 0) {
249  ATH_MSG_DEBUG("Skipping since efield = 0 for x_pix = " << x_pix << " y_pix = " << y_pix);
250  continue;
251  }
252 
253  const double mobilityElectron = getMobility(efield, false);
254  const double mobilityHole = getMobility(efield, true);
255  auto driftTimeElectron = getDriftTime(false, ncharges, rndmEngine, trappingTimeElectrons, trappingTimeHoles);
256  auto driftTimeHole = getDriftTime(true, ncharges, rndmEngine, trappingTimeElectrons, trappingTimeHoles);
257  //Need to determine how many elementary charges this charge chunk represents.
258  double chunk_size = energy_per_step * eleholePairEnergy; //number of electrons/holes
259  //set minimum limit to prevent dividing into smaller subcharges than one fundamental charge
260  if (chunk_size < 1) chunk_size = 1;
261  const double kappa = 1. / std::sqrt(chunk_size);
262 
263  const double timeToElectrodeElectron = timeMap_e.getContent(timeMap_e.getBinX(1e3*y_pix), timeMap_e.getBinY(1e3*x_pix));
264  const double timeToElectrodeHole = timeMap_h.getContent(timeMap_h.getBinX(1e3*y_pix), timeMap_h.getBinY(1e3*x_pix));
265 
266  /* Diffusion via the Einstein relation
267  D = mu * kB * T / q
268  D = (mu / mm^2/MV*ns) * (T/273 K) * 0.024 microns^2 / ns */
269  const auto prefactor_e = mobilityElectron*0.024*m_temperature / 273.;
270  const auto prefactor_h = mobilityHole*0.024*m_temperature / 273.;
271  //Apply diffusion. rdif is the max. diffusion
272  for(size_t i = 0; i < ncharges; ++i) {
273  DtElectron[i] = prefactor_e * std::min(driftTimeElectron[i], timeToElectrodeElectron);
274  DtHole[i] = prefactor_h * std::min(driftTimeHole[i], timeToElectrodeHole);
275  rdifElectron[i] = 1e-3*std::sqrt(DtElectron[i]); //in mm
276  rdifHole[i] = 1e-3*std::sqrt(DtHole[i]); //in mm
277  }
278 
279  const float average_chargeElectron = avgChargeMap_e.getContent(avgChargeMap_e.getBinY(1e3*y_pix), avgChargeMap_e.getBinX(1e3*x_pix));
280  const float average_chargeHole = avgChargeMap_h.getContent(avgChargeMap_h.getBinY(1e3*y_pix), avgChargeMap_h.getBinX(1e3*x_pix));
281 
282  //We're sticking to the "old" convention here where there was a loop over -1/0/1 for both directions
283  //(x/y) to the neighbouring pixels. We call them (p)lus, (z)ero, (m)inus in either i- or j- direction
284  //hence giving us plus-in-i or pi or zero-in-j or zj etc.
285  //i is in direction of x, whereas j is into the y direction, this code is ugly, but given that we don't have very few branching conditions
286  //it's pretty fast
287 
288  //This was: ramoPotentialMap.getBinY(1000*(x_pix + pixel_size_x * 3 - i * pixel_size_x)
289  const std::size_t ramo_init_bin_y_pi = ramoPotentialMap.getBinY(1000*(x_pix + pixel_size_x * 2));
290  const std::size_t ramo_init_bin_y_zi = ramoPotentialMap.getBinY(1000*(x_pix + pixel_size_x * 3));
291  const std::size_t ramo_init_bin_y_mi = ramoPotentialMap.getBinY(1000*(x_pix + pixel_size_x * 4));
292 
293  //This was: ramoPotentialMap.getBinX(1000*(y_pix + 0.5*pixel_size_y - j * pixel_size_y)
294  const std::size_t ramo_init_bin_x_pj = ramoPotentialMap.getBinX(1000*(y_pix - 0.5*pixel_size_y));
295  const std::size_t ramo_init_bin_x_zj = ramoPotentialMap.getBinX(1000*(y_pix + 0.5*pixel_size_y));
296  const std::size_t ramo_init_bin_x_mj = ramoPotentialMap.getBinX(1000*(y_pix + 1.5*pixel_size_y));
297 
298  //For some reason the x- and y-indices for some values are swapped, this is extremely confusing
299  float ramoInit_pipj = ramoPotentialMap.getContent(ramo_init_bin_x_pj, ramo_init_bin_y_pi);
300  float ramoInit_pizj = ramoPotentialMap.getContent(ramo_init_bin_x_zj, ramo_init_bin_y_pi);
301  float ramoInit_pimj = ramoPotentialMap.getContent(ramo_init_bin_x_mj, ramo_init_bin_y_pi);
302  float ramoInit_zipj = ramoPotentialMap.getContent(ramo_init_bin_x_pj, ramo_init_bin_y_zi);
303  float ramoInit_zizj = ramoPotentialMap.getContent(ramo_init_bin_x_zj, ramo_init_bin_y_zi);
304  float ramoInit_zimj = ramoPotentialMap.getContent(ramo_init_bin_x_mj, ramo_init_bin_y_zi);
305  float ramoInit_mipj = ramoPotentialMap.getContent(ramo_init_bin_x_pj, ramo_init_bin_y_mi);
306  float ramoInit_mizj = ramoPotentialMap.getContent(ramo_init_bin_x_zj, ramo_init_bin_y_mi);
307  float ramoInit_mimj = ramoPotentialMap.getContent(ramo_init_bin_x_mj, ramo_init_bin_y_mi);
308 
309  const auto hit_time = hitTime(phit);
310 
311  //Loop over charge-carrier pairs, we're looping over electrons and holes at the same time
312  for (int j = 0; j < ncharges; j++) {
313  // In case the charge moves into a neighboring pixel
314  int extraNPixXElectron = nPixX;
315  int extraNPixYElectron = nPixY;
316  int extraNPixXHole = nPixX;
317  int extraNPixYHole = nPixY;
318 
319  //Apply drift due to diffusion
320  std::array<double, 4> randomNumbers{};
321  CLHEP::RandGaussZiggurat::shootArray(rndmEngine, 4, randomNumbers.data());
322 
323  double xposDiffElectron = x_pix + rdifElectron[j] * randomNumbers[0];
324  double yposDiffElectron = y_pix + rdifElectron[j] * randomNumbers[1];
325  double xposDiffHole = x_pix + rdifHole[j] * randomNumbers[2];
326  double yposDiffHole = y_pix + rdifHole[j] * randomNumbers[3];
327 
328  // Account for drifting into another pixel
329  while (xposDiffElectron > pixel_size_x) {
330  extraNPixXElectron++; // increments or decrements pixel count in x
331  xposDiffElectron -= pixel_size_x; // moves xpos coordinate 1 pixel over in x
332  }
333  while (xposDiffElectron < 0) {
334  extraNPixXElectron--;
335  xposDiffElectron += pixel_size_x;
336  }
337  while (yposDiffElectron > pixel_size_y) {
338  extraNPixYElectron++; // increments or decrements pixel count in y
339  yposDiffElectron -= pixel_size_y; // moves xpos coordinate 1 pixel over in y
340  }
341  while (yposDiffElectron < 0) {
342  extraNPixYElectron--;
343  yposDiffElectron += pixel_size_y;
344  }
345 
346  //And drifting for for holes
347  while (xposDiffHole > pixel_size_x) {
348  extraNPixXHole++;
349  xposDiffHole -= pixel_size_x;
350  }
351  while (xposDiffHole < 0) {
352  extraNPixXHole--;
353  xposDiffHole += pixel_size_x;
354  }
355  while (yposDiffHole > pixel_size_y) {
356  extraNPixYHole++;
357  yposDiffHole -= pixel_size_y;
358  }
359  while (yposDiffHole < 0) {
360  extraNPixYHole--;
361  yposDiffHole += pixel_size_y;
362  }
363  auto xposFinalElectron = yPositionMap_e.getContent(
364  yPositionMap_e.getBinX(1e3*yposDiffElectron),
365  yPositionMap_e.getBinY(1e3*xposDiffElectron),
366  yPositionMap_e.getBinZ(std::min(driftTimeElectron[j],timeToElectrodeElectron))
367  ) * 1e-3; //[mm]
368  auto yposFinalElectron = xPositionMap_e.getContent(
369  xPositionMap_e.getBinX(1e3*yposDiffElectron),
370  xPositionMap_e.getBinY(1e3*xposDiffElectron),
371  xPositionMap_e.getBinZ(std::min(driftTimeElectron[j],timeToElectrodeElectron))
372  ) * 1e-3; //[mm]
373 
374  auto xposFinalHole = yPositionMap_h.getContent(
375  yPositionMap_h.getBinX(1e3*yposDiffHole),
376  yPositionMap_h.getBinY(1e3*xposDiffHole),
377  yPositionMap_h.getBinZ(std::min(driftTimeHole[j],timeToElectrodeHole))
378  ) * 1e-3; //[mm]
379  auto yposFinalHole = xPositionMap_h.getContent(
380  xPositionMap_h.getBinX(1e3*yposDiffHole),
381  xPositionMap_h.getBinY(1e3*xposDiffHole),
382  xPositionMap_h.getBinZ(std::min(driftTimeHole[j],timeToElectrodeHole))
383  ) * 1e-3; //[mm]
384 
385  // -- Calculate signal in current pixel and in the neighboring ones
386  //This was: ramoPotentialMap.getBinY(1000*(xposFinal + pixel_size_x * 3 - i * pixel_size_x)
387  const std::size_t ramo_final_bin_y_pi_electrons = ramoPotentialMap.getBinY(1000*(xposFinalElectron + pixel_size_x * 2));
388  const std::size_t ramo_final_bin_y_pi_holes = ramoPotentialMap.getBinY(1000*(xposFinalHole + pixel_size_x * 2));
389  const std::size_t ramo_final_bin_y_zi_electrons = ramoPotentialMap.getBinY(1000*(xposFinalElectron + pixel_size_x * 3));
390  const std::size_t ramo_final_bin_y_zi_holes = ramoPotentialMap.getBinY(1000*(xposFinalHole + pixel_size_x * 3));
391  const std::size_t ramo_final_bin_y_mi_electrons = ramoPotentialMap.getBinY(1000*(xposFinalElectron + pixel_size_x * 4));
392  const std::size_t ramo_final_bin_y_mi_holes = ramoPotentialMap.getBinY(1000*(xposFinalHole + pixel_size_x * 4));
393  //This was: ramoPotentialMap.getBinX(1000*(yposFinal + 0.5*pixel_size_y - j * pixel_size_y)
394  const std::size_t ramo_final_bin_x_pj_electrons = ramoPotentialMap.getBinX(1000*(yposFinalElectron - 0.5*pixel_size_y));
395  const std::size_t ramo_final_bin_x_pj_holes = ramoPotentialMap.getBinX(1000*(yposFinalHole - 0.5*pixel_size_y));
396  const std::size_t ramo_final_bin_x_zj_electrons = ramoPotentialMap.getBinX(1000*(yposFinalElectron + 0.5*pixel_size_y));
397  const std::size_t ramo_final_bin_x_zj_holes = ramoPotentialMap.getBinX(1000*(yposFinalHole + 0.5*pixel_size_y));
398  const std::size_t ramo_final_bin_x_mj_electrons = ramoPotentialMap.getBinX(1000*(yposFinalElectron + 1.5*pixel_size_y));
399  const std::size_t ramo_final_bin_x_mj_holes = ramoPotentialMap.getBinX(1000*(yposFinalHole + 1.5*pixel_size_y));
400 
401  float ramoFinal_pipj_electrons = ramoPotentialMap.getContent(ramo_final_bin_x_pj_electrons, ramo_final_bin_y_pi_electrons);
402  float ramoFinal_pizj_electrons = ramoPotentialMap.getContent(ramo_final_bin_x_zj_electrons, ramo_final_bin_y_pi_electrons);
403  float ramoFinal_pimj_electrons = ramoPotentialMap.getContent(ramo_final_bin_x_mj_electrons, ramo_final_bin_y_pi_electrons);
404  float ramoFinal_zipj_electrons = ramoPotentialMap.getContent(ramo_final_bin_x_pj_electrons, ramo_final_bin_y_zi_electrons);
405  float ramoFinal_zizj_electrons = ramoPotentialMap.getContent(ramo_final_bin_x_zj_electrons, ramo_final_bin_y_zi_electrons);
406  float ramoFinal_zimj_electrons = ramoPotentialMap.getContent(ramo_final_bin_x_mj_electrons, ramo_final_bin_y_zi_electrons);
407  float ramoFinal_mipj_electrons = ramoPotentialMap.getContent(ramo_final_bin_x_pj_electrons, ramo_final_bin_y_mi_electrons);
408  float ramoFinal_mizj_electrons = ramoPotentialMap.getContent(ramo_final_bin_x_zj_electrons, ramo_final_bin_y_mi_electrons);
409  float ramoFinal_mimj_electrons = ramoPotentialMap.getContent(ramo_final_bin_x_mj_electrons, ramo_final_bin_y_mi_electrons);
410 
411  float ramoFinal_pipj_holes = ramoPotentialMap.getContent(ramo_final_bin_x_pj_holes, ramo_final_bin_y_pi_holes);
412  float ramoFinal_pizj_holes = ramoPotentialMap.getContent(ramo_final_bin_x_zj_holes, ramo_final_bin_y_pi_holes);
413  float ramoFinal_pimj_holes = ramoPotentialMap.getContent(ramo_final_bin_x_mj_holes, ramo_final_bin_y_pi_holes);
414  float ramoFinal_zipj_holes = ramoPotentialMap.getContent(ramo_final_bin_x_pj_holes, ramo_final_bin_y_zi_holes);
415  float ramoFinal_zizj_holes = ramoPotentialMap.getContent(ramo_final_bin_x_zj_holes, ramo_final_bin_y_zi_holes);
416  float ramoFinal_zimj_holes = ramoPotentialMap.getContent(ramo_final_bin_x_mj_holes, ramo_final_bin_y_zi_holes);
417  float ramoFinal_mipj_holes = ramoPotentialMap.getContent(ramo_final_bin_x_pj_holes, ramo_final_bin_y_mi_holes);
418  float ramoFinal_mizj_holes = ramoPotentialMap.getContent(ramo_final_bin_x_zj_holes, ramo_final_bin_y_mi_holes);
419  float ramoFinal_mimj_holes = ramoPotentialMap.getContent(ramo_final_bin_x_mj_holes, ramo_final_bin_y_mi_holes);
420 
421  // Record deposit
422  double eHitRamo_pipj_electrons = energy_per_step * (ramoFinal_pipj_electrons - ramoInit_pipj);
423  double eHitRamo_pipj_holes = -1. * energy_per_step * (ramoFinal_pipj_holes - ramoInit_pipj);
424  double eHitRamo_pizj_electrons = energy_per_step * (ramoFinal_pizj_electrons - ramoInit_pizj);
425  double eHitRamo_pizj_holes = -1. * energy_per_step * (ramoFinal_pizj_holes - ramoInit_pizj);
426  double eHitRamo_pimj_electrons = energy_per_step * (ramoFinal_pimj_electrons - ramoInit_pimj);
427  double eHitRamo_pimj_holes = -1. * energy_per_step * (ramoFinal_pimj_holes - ramoInit_pimj);
428  double eHitRamo_zipj_electrons = energy_per_step * (ramoFinal_zipj_electrons - ramoInit_zipj);
429  double eHitRamo_zipj_holes = -1. * energy_per_step * (ramoFinal_zipj_holes - ramoInit_zipj);
430  double eHitRamo_zizj_electrons = energy_per_step * (ramoFinal_zizj_electrons - ramoInit_zizj);
431  double eHitRamo_zizj_holes = -1. * energy_per_step * (ramoFinal_zizj_holes - ramoInit_zizj);
432  double eHitRamo_zimj_electrons = energy_per_step * (ramoFinal_zimj_electrons - ramoInit_zimj);
433  double eHitRamo_zimj_holes = -1. * energy_per_step * (ramoFinal_zimj_holes - ramoInit_zimj);
434  double eHitRamo_mipj_electrons = energy_per_step * (ramoFinal_mipj_electrons - ramoInit_mipj);
435  double eHitRamo_mipj_holes = -1. * energy_per_step * (ramoFinal_mipj_holes - ramoInit_mipj);
436  double eHitRamo_mizj_electrons = energy_per_step * (ramoFinal_mizj_electrons - ramoInit_mizj);
437  double eHitRamo_mizj_holes = -1. * energy_per_step * (ramoFinal_mizj_holes - ramoInit_mizj);
438  double eHitRamo_mimj_electrons = energy_per_step * (ramoFinal_mimj_electrons - ramoInit_mimj);
439  double eHitRamo_mimj_holes = -1. * energy_per_step * (ramoFinal_mimj_holes - ramoInit_mimj);
440 
441  if(doChunkCorrection) {
442  eHitRamo_pipj_electrons = energy_per_step * average_chargeElectron + kappa * (eHitRamo_pipj_electrons - energy_per_step * average_chargeElectron);
443  eHitRamo_pipj_holes = energy_per_step * average_chargeHole + kappa * (eHitRamo_pipj_holes - energy_per_step * average_chargeHole);
444  eHitRamo_pizj_electrons = energy_per_step * average_chargeElectron + kappa * (eHitRamo_pizj_electrons - energy_per_step * average_chargeElectron);
445  eHitRamo_pizj_holes = energy_per_step * average_chargeHole + kappa * (eHitRamo_pizj_holes - energy_per_step * average_chargeHole);
446  eHitRamo_pimj_electrons = energy_per_step * average_chargeElectron + kappa * (eHitRamo_pimj_electrons - energy_per_step * average_chargeElectron);
447  eHitRamo_pimj_holes = energy_per_step * average_chargeHole + kappa * (eHitRamo_pimj_holes - energy_per_step * average_chargeHole);
448  eHitRamo_zipj_electrons = energy_per_step * average_chargeElectron + kappa * (eHitRamo_zipj_electrons - energy_per_step * average_chargeElectron);
449  eHitRamo_zipj_holes = energy_per_step * average_chargeHole + kappa * (eHitRamo_zipj_holes - energy_per_step * average_chargeHole);
450  eHitRamo_zizj_electrons = energy_per_step * average_chargeElectron + kappa * (eHitRamo_zizj_electrons - energy_per_step * average_chargeElectron);
451  eHitRamo_zizj_holes = energy_per_step * average_chargeHole + kappa * (eHitRamo_zizj_holes - energy_per_step * average_chargeHole);
452  eHitRamo_zimj_electrons = energy_per_step * average_chargeElectron + kappa * (eHitRamo_zimj_electrons - energy_per_step * average_chargeElectron);
453  eHitRamo_pimj_holes = energy_per_step * average_chargeHole + kappa * (eHitRamo_pimj_holes - energy_per_step * average_chargeHole);
454  eHitRamo_mipj_electrons = energy_per_step * average_chargeElectron + kappa * (eHitRamo_mipj_electrons - energy_per_step * average_chargeElectron);
455  eHitRamo_mipj_holes = energy_per_step * average_chargeHole + kappa * (eHitRamo_mipj_holes - energy_per_step * average_chargeHole);
456  eHitRamo_mizj_electrons = energy_per_step * average_chargeElectron + kappa * (eHitRamo_mizj_electrons - energy_per_step * average_chargeElectron);
457  eHitRamo_mizj_holes = energy_per_step * average_chargeHole + kappa * (eHitRamo_mizj_holes - energy_per_step * average_chargeHole);
458  eHitRamo_mimj_electrons = energy_per_step * average_chargeElectron + kappa * (eHitRamo_mimj_electrons - energy_per_step * average_chargeElectron);
459  eHitRamo_mimj_holes = energy_per_step * average_chargeHole + kappa * (eHitRamo_mimj_holes - energy_per_step * average_chargeHole);
460  }
461 
462  // -- pixel coordinates --> module coordinates
463  //This was: double x_mod = x_pix + pixel_size_x * extraNPixX - 0.5*module_size_x + i*pixel_size_x;
464  //This was: double y_mod = y_pix + pixel_size_y * extraNPixY - 0.5*module_size_y + j*pixel_size_y;
465  double x_mod_pi_electrons = x_pix + pixel_size_x * extraNPixXElectron - 0.5*module_size_x + pixel_size_x;
466  double y_mod_pj_electrons = y_pix + pixel_size_y * extraNPixYElectron - 0.5*module_size_y + pixel_size_y;
467  double x_mod_pi_holes = x_pix + pixel_size_x * extraNPixXHole - 0.5*module_size_x + pixel_size_x;
468  double y_mod_pj_holes = y_pix + pixel_size_y * extraNPixYHole - 0.5*module_size_y + pixel_size_y;
469  double x_mod_zi_electrons = x_pix + pixel_size_x * extraNPixXElectron - 0.5*module_size_x;
470  double y_mod_zj_electrons = y_pix + pixel_size_y * extraNPixYElectron - 0.5*module_size_y;
471  double x_mod_zi_holes = x_pix + pixel_size_x * extraNPixXHole - 0.5*module_size_x;
472  double y_mod_zj_holes = y_pix + pixel_size_y * extraNPixYHole - 0.5*module_size_y;
473  double x_mod_mi_electrons = x_pix + pixel_size_x * extraNPixXElectron - 0.5*module_size_x - pixel_size_x;
474  double y_mod_mj_electrons = y_pix + pixel_size_y * extraNPixYElectron - 0.5*module_size_y - pixel_size_y;
475  double x_mod_mi_holes = x_pix + pixel_size_x * extraNPixXHole - 0.5*module_size_x - pixel_size_x;
476  double y_mod_mj_holes = y_pix + pixel_size_y * extraNPixYHole - 0.5*module_size_y - pixel_size_y;
477 
478  //This was: const SiLocalPosition& chargePos = Module.hitLocalToLocal(y_mod, x_mod) (for whatever reason half the maps x/y are inverted...)
479  const SiLocalPosition& chargePos_pipj_electrons = Module.hitLocalToLocal(y_mod_pj_electrons, x_mod_pi_electrons);
480  const SiLocalPosition& chargePos_pizj_electrons = Module.hitLocalToLocal(y_mod_zj_electrons, x_mod_pi_electrons);
481  const SiLocalPosition& chargePos_pimj_electrons = Module.hitLocalToLocal(y_mod_mj_electrons, x_mod_pi_electrons);
482  const SiLocalPosition& chargePos_zipj_electrons = Module.hitLocalToLocal(y_mod_pj_electrons, x_mod_zi_electrons);
483  const SiLocalPosition& chargePos_zizj_electrons = Module.hitLocalToLocal(y_mod_zj_electrons, x_mod_zi_electrons);
484  const SiLocalPosition& chargePos_zimj_electrons = Module.hitLocalToLocal(y_mod_mj_electrons, x_mod_zi_electrons);
485  const SiLocalPosition& chargePos_mipj_electrons = Module.hitLocalToLocal(y_mod_pj_electrons, x_mod_mi_electrons);
486  const SiLocalPosition& chargePos_mizj_electrons = Module.hitLocalToLocal(y_mod_zj_electrons, x_mod_mi_electrons);
487  const SiLocalPosition& chargePos_mimj_electrons = Module.hitLocalToLocal(y_mod_mj_electrons, x_mod_mi_electrons);
488 
489  const SiLocalPosition& chargePos_pipj_holes = Module.hitLocalToLocal(y_mod_pj_holes, x_mod_pi_holes);
490  const SiLocalPosition& chargePos_pizj_holes = Module.hitLocalToLocal(y_mod_zj_holes, x_mod_pi_holes);
491  const SiLocalPosition& chargePos_pimj_holes = Module.hitLocalToLocal(y_mod_mj_holes, x_mod_pi_holes);
492  const SiLocalPosition& chargePos_zipj_holes = Module.hitLocalToLocal(y_mod_pj_holes, x_mod_zi_holes);
493  const SiLocalPosition& chargePos_zizj_holes = Module.hitLocalToLocal(y_mod_zj_holes, x_mod_zi_holes);
494  const SiLocalPosition& chargePos_zimj_holes = Module.hitLocalToLocal(y_mod_mj_holes, x_mod_zi_holes);
495  const SiLocalPosition& chargePos_mipj_holes = Module.hitLocalToLocal(y_mod_pj_holes, x_mod_mi_holes);
496  const SiLocalPosition& chargePos_mizj_holes = Module.hitLocalToLocal(y_mod_zj_holes, x_mod_mi_holes);
497  const SiLocalPosition& chargePos_mimj_holes = Module.hitLocalToLocal(y_mod_mj_holes, x_mod_mi_holes);
498 
499  //This was: const SiSurfaceCharge scharge(chargePos, SiCharge(eHitRamo*eleholePairEnergy, hit_time, SiCharge::track, particleLink));
500  const SiSurfaceCharge scharge_pipj_electrons(chargePos_pipj_electrons, SiCharge(eHitRamo_pipj_electrons*eleholePairEnergy, hit_time, SiCharge::track, particleLink));
501  const SiSurfaceCharge scharge_pizj_electrons(chargePos_pizj_electrons, SiCharge(eHitRamo_pizj_electrons*eleholePairEnergy, hit_time, SiCharge::track, particleLink));
502  const SiSurfaceCharge scharge_pimj_electrons(chargePos_pimj_electrons, SiCharge(eHitRamo_pimj_electrons*eleholePairEnergy, hit_time, SiCharge::track, particleLink));
503  const SiSurfaceCharge scharge_zipj_electrons(chargePos_zipj_electrons, SiCharge(eHitRamo_zipj_electrons*eleholePairEnergy, hit_time, SiCharge::track, particleLink));
504  const SiSurfaceCharge scharge_zizj_electrons(chargePos_zizj_electrons, SiCharge(eHitRamo_zizj_electrons*eleholePairEnergy, hit_time, SiCharge::track, particleLink));
505  const SiSurfaceCharge scharge_zimj_electrons(chargePos_zimj_electrons, SiCharge(eHitRamo_zimj_electrons*eleholePairEnergy, hit_time, SiCharge::track, particleLink));
506  const SiSurfaceCharge scharge_mipj_electrons(chargePos_mipj_electrons, SiCharge(eHitRamo_mipj_electrons*eleholePairEnergy, hit_time, SiCharge::track, particleLink));
507  const SiSurfaceCharge scharge_mizj_electrons(chargePos_mizj_electrons, SiCharge(eHitRamo_mizj_electrons*eleholePairEnergy, hit_time, SiCharge::track, particleLink));
508  const SiSurfaceCharge scharge_mimj_electrons(chargePos_mimj_electrons, SiCharge(eHitRamo_mimj_electrons*eleholePairEnergy, hit_time, SiCharge::track, particleLink));
509 
510  const SiSurfaceCharge scharge_pipj_holes(chargePos_pipj_holes, SiCharge(eHitRamo_pipj_holes*eleholePairEnergy, hit_time, SiCharge::track, particleLink));
511  const SiSurfaceCharge scharge_pizj_holes(chargePos_pizj_holes, SiCharge(eHitRamo_pizj_holes*eleholePairEnergy, hit_time, SiCharge::track, particleLink));
512  const SiSurfaceCharge scharge_pimj_holes(chargePos_pimj_holes, SiCharge(eHitRamo_pimj_holes*eleholePairEnergy, hit_time, SiCharge::track, particleLink));
513  const SiSurfaceCharge scharge_zipj_holes(chargePos_zipj_holes, SiCharge(eHitRamo_zipj_holes*eleholePairEnergy, hit_time, SiCharge::track, particleLink));
514  const SiSurfaceCharge scharge_zizj_holes(chargePos_zizj_holes, SiCharge(eHitRamo_zizj_holes*eleholePairEnergy, hit_time, SiCharge::track, particleLink));
515  const SiSurfaceCharge scharge_zimj_holes(chargePos_zimj_holes, SiCharge(eHitRamo_zimj_holes*eleholePairEnergy, hit_time, SiCharge::track, particleLink));
516  const SiSurfaceCharge scharge_mipj_holes(chargePos_mipj_holes, SiCharge(eHitRamo_mipj_holes*eleholePairEnergy, hit_time, SiCharge::track, particleLink));
517  const SiSurfaceCharge scharge_mizj_holes(chargePos_mizj_holes, SiCharge(eHitRamo_mizj_holes*eleholePairEnergy, hit_time, SiCharge::track, particleLink));
518  const SiSurfaceCharge scharge_mimj_holes(chargePos_mimj_holes, SiCharge(eHitRamo_mimj_holes*eleholePairEnergy, hit_time, SiCharge::track, particleLink));
519 
520  auto addCharge = [&Module, &chargedDiodes](SiSurfaceCharge const & scharge) {
521  const SiCellId& diode = Module.cellIdOfPosition(scharge.position());
522  if (diode.isValid()) {
523  const SiCharge& charge = scharge.charge();
524  chargedDiodes.add(diode, charge);
525  }
526  };
527 
528  addCharge(scharge_pipj_electrons);
529  addCharge(scharge_pizj_electrons);
530  addCharge(scharge_pimj_electrons);
531  addCharge(scharge_zipj_electrons);
532  addCharge(scharge_zizj_electrons);
533  addCharge(scharge_zimj_electrons);
534  addCharge(scharge_mipj_electrons);
535  addCharge(scharge_mizj_electrons);
536  addCharge(scharge_mimj_electrons);
537  addCharge(scharge_pipj_holes);
538  addCharge(scharge_pizj_holes);
539  addCharge(scharge_pimj_holes);
540  addCharge(scharge_zipj_holes);
541  addCharge(scharge_zizj_holes);
542  addCharge(scharge_zimj_holes);
543  addCharge(scharge_mipj_holes);
544  addCharge(scharge_mizj_holes);
545  addCharge(scharge_mimj_holes);
546  }
547  }
548  } else if (m_radiationDamageSimulationType == RadiationDamageSimulationType::TEMPLATE_CORRECTION) {
549  const PixelHistoConverter& chargeCorrectionHist = m_chargeCorrection[layerIndex];
550 
551  for (const auto& iHitRecord : trfHitRecord) {
552  double eta_i = eta_0;
553  double phi_i = phi_0;
554  double depth_i = depth_0;
555 
556  if (iTotalLength) {
557  eta_i += 1.0 * iHitRecord.first / iTotalLength * dEta;
558  phi_i += 1.0 * iHitRecord.first / iTotalLength * dPhi;
559  depth_i += 1.0 * iHitRecord.first / iTotalLength * dDepth;
560  }
561 
562  double es_current = 1.0 * iHitRecord.second / 1.E+6;
563 
564  double dist_electrode = 0.5 * sensorThickness - Module.design().readoutSide() * depth_i;
565  if (dist_electrode < 0) dist_electrode = 0;
566 
567  CLHEP::Hep3Vector chargepos;
568  chargepos.setX(phi_i);
569  chargepos.setY(eta_i);
570  chargepos.setZ(dist_electrode);
571 
572  bool coord = Module.isModuleFrame();
573 
575  "ismoduleframe " << coord << " -- startPosition (x,y,z) = " << chargepos.x() << ", " << chargepos.y() << ", " <<
576  chargepos.z());
577 
578  // -- change origin of coordinates to the left bottom of module
579  const double x_new = chargepos.x() + 0.5*module_size_x;
580  const double y_new = chargepos.y() + 0.5*module_size_y;
581 
582  // -- change from module frame to pixel frame
583  const int nPixX = int(x_new / pixel_size_x);
584  const int nPixY = int(y_new / pixel_size_y);
585  ATH_MSG_DEBUG(" -- nPixX = " << nPixX << " nPixY = " << nPixY);
586  const double x_pix = x_new - pixel_size_x * (nPixX);
587  const double y_pix = y_new - pixel_size_y * (nPixY);
588  // -- change origin of coordinates to the center of the pixel
589  const double x_pix_center = x_pix - pixel_size_x / 2;
590  const double y_pix_center = y_pix - pixel_size_y / 2;
591  ATH_MSG_DEBUG(" -- current hit position w.r.t. pixel center = " << x_pix_center << " " << y_pix_center);
592 
593  // radiation damage correction
594  const double distance = std::hypot(x_pix_center, y_pix_center) / Gaudi::Units::micrometer; // to get it in micrometers
595  const double chargeCorrection = chargeCorrectionHist.getContent(distance);
596 
597  // calculate the charge at the incident sensor
598  // no neighbours are considered as that is allready included in the template correction
599  const double ed = es_current * eleholePairEnergy * chargeCorrection;
600 
601  // -- pixel coordinates --> module coordinates
602  const double x_mod = x_pix_center + 0.5*pixel_size_x + pixel_size_x * nPixX - 0.5*module_size_x;
603  const double y_mod = y_pix_center + 0.5*pixel_size_y + pixel_size_y * nPixY - 0.5*module_size_y;
604  const SiLocalPosition& chargePos = Module.hitLocalToLocal(y_mod, x_mod);
605 
606  const SiSurfaceCharge scharge(chargePos,
607  SiCharge(ed, pHitTime, SiCharge::track, particleLink));
608  const SiCellId& diode = Module.cellIdOfPosition(scharge.position());
609  if (diode.isValid()) {
610  const SiCharge& charge = scharge.charge();
611  chargedDiodes.add(diode, charge);
612  }
613  }
614  } else {
615  //**************************************//
616  //*** Now diffuse charges to surface *** //
617  //**************************************//
618  for (const auto& iHitRecord : trfHitRecord) {
619  double eta_i = eta_0;
620  double phi_i = phi_0;
621  double depth_i = depth_0;
622 
623  if (iTotalLength) {
624  eta_i += 1.0 * iHitRecord.first / iTotalLength * dEta;
625  phi_i += 1.0 * iHitRecord.first / iTotalLength * dPhi;
626  depth_i += 1.0 * iHitRecord.first / iTotalLength * dDepth;
627  }
628 
629  double es_current = 1.0 * iHitRecord.second / 1.E+6;
630 
631  double dist_electrode = 0.5 * sensorThickness - Module.design().readoutSide() * depth_i;
632  if (dist_electrode < 0) dist_electrode = 0;
633 
634  CLHEP::Hep3Vector chargepos;
635  chargepos.setX(phi_i);
636  chargepos.setY(eta_i);
637  chargepos.setZ(dist_electrode);
638 
639  bool coord = Module.isModuleFrame();
640 
642  "ismoduleframe " << coord << " -- startPosition (x,y,z) = " << chargepos.x() << ", " << chargepos.y() << ", " <<
643  chargepos.z());
644 
645  // -- change origin of coordinates to the left bottom of module
646  double x_new = chargepos.x() + 0.5*module_size_x;
647  double y_new = chargepos.y() + 0.5*module_size_y;
648 
649  // -- change from module frame to pixel frame
650  int nPixX = int(x_new / pixel_size_x);
651  int nPixY = int(y_new / pixel_size_y);
652  ATH_MSG_DEBUG(" -- nPixX = " << nPixX << " nPixY = " << nPixY);
653  double x_pix = x_new - pixel_size_x * (nPixX);
654  double y_pix = y_new - pixel_size_y * (nPixY);
655  // -- change origin of coordinates to the center of the pixel
656  double x_pix_center = x_pix - pixel_size_x / 2;
657  double y_pix_center = y_pix - pixel_size_y / 2;
658  ATH_MSG_DEBUG(" -- current hit position w.r.t. pixel center = " << x_pix_center << " " << y_pix_center);
659 
660  double x_neighbor;
661  double y_neighbor;
662  CLHEP::Hep3Vector pos_neighbor;
663  // -- Calculate signal in current pixel and in the neighboring ones
664  // -- loop in the x-coordinate
665  for (int i = -1; i <= 1; i++) {
666  x_neighbor = x_pix_center - i * pixel_size_x;
667  // -- loop in the y-coordinate
668  for (int j = -1; j <= 1; j++) {
669  y_neighbor = y_pix_center - j * pixel_size_y;
670 
671  // -- check if the neighbor falls inside the charge collection prob map window
672  if ((std::abs(x_neighbor) < pixel_size_x) && (std::abs(y_neighbor) < pixel_size_y)) {
673  // -- change origin of coordinates to the bottom left of the charge
674  // collection prob map "window", i.e. shift of 1-pixel twd bottom left
675  double x_neighbor_map = x_neighbor + pixel_size_x;
676  double y_neighbor_map = y_neighbor + pixel_size_y;
677 
678  int x_bin_cc_map = x_neighbor_map / x_bin_size;
679  int y_bin_cc_map = y_neighbor_map / y_bin_size;
680 
681  // -- retrieve the charge collection probability from Svc
682  // -- swap x and y bins to match Map coord convention
683  double ccprob_neighbor = getProbMapEntry(InDetDD::PixelReadoutTechnology::FEI4, y_bin_cc_map, x_bin_cc_map);
684  if (ccprob_neighbor == -1.) return StatusCode::FAILURE;
685 
686  double ed = es_current * eleholePairEnergy * ccprob_neighbor;
687 
688  // -- pixel coordinates --> module coordinates
689  double x_mod = x_neighbor + 0.5*pixel_size_x + pixel_size_x * nPixX - 0.5*module_size_x;
690  double y_mod = y_neighbor + 0.5*pixel_size_y + pixel_size_y * nPixY - 0.5*module_size_y;
691  const SiLocalPosition& chargePos = Module.hitLocalToLocal(y_mod, x_mod);
692 
693  const SiSurfaceCharge scharge(chargePos,
694  SiCharge(ed, pHitTime, SiCharge::track, particleLink));
695  const SiCellId& diode = Module.cellIdOfPosition(scharge.position());
696  if (diode.isValid()) {
697  const SiCharge& charge = scharge.charge();
698  chargedDiodes.add(diode, charge);
699  }
700  }
701  }
702  }
703  }
704  }
705 
706  return StatusCode::SUCCESS;
707 }
708 
709 // read the Charge Collection Prob Map from text file
710 StatusCode SensorSim3DTool::readProbMap(const std::string& fileE) {
711  std::string line;
712  const std::string& fileName = fileE;
714  if (inputFile.empty()) {
715  ATH_MSG_ERROR("Could not open input file!!!!!");
716  return StatusCode::FAILURE;
717  }
718  ATH_MSG_DEBUG("opening file " << inputFile);
719  std::ifstream myfile(inputFile.c_str());
720  if (myfile.is_open()) {
721  ATH_MSG_DEBUG(" opened file!");
722  while (myfile.good()) {
723  while (std::getline(myfile, line)) {
724  std::istringstream sline(line);
725  int xpos, ypos;
726  double prob;
727  sline >> xpos >> ypos >> prob;
728  if (fileName.find("FEI4") != std::string::npos) {
729  m_probMapFEI4.insert(std::make_pair(std::make_pair(xpos, ypos), prob));
730  ATH_MSG_DEBUG("FEI4 inside xpos " << xpos << " ypos " << ypos << " prob " << prob);
731  } else if (fileName.find("FEI3") != std::string::npos) {
732  m_probMapFEI3.insert(std::make_pair(std::make_pair(xpos, ypos), prob));
733  ATH_MSG_DEBUG("FEI3 inside xpos " << xpos << " ypos " << ypos << " prob " << prob);
734  } else {
735  ATH_MSG_ERROR("Please check name of Charge Coll Prob Maps! (should contain FEI3 or FEI4) ");
736  return StatusCode::FAILURE;
737  }
738  }
739  }
740  myfile.close();
741  }
742  return StatusCode::SUCCESS;
743 }
744 
745 // -- Print out the Charge Collection Probability map (full map)
746 StatusCode SensorSim3DTool::printProbMap(const std::string& readout) const {
747  if (readout == "FEI4") {
748  for (const auto& it : m_probMapFEI4) {
750  "read full probMap FEI4 --- bin x " << it.first.first << " bin y " << it.first.second << " prob " <<
751  it.second);
752  }
753  } else if (readout == "FEI3") {
754  for (const auto& it : m_probMapFEI3) {
756  "read full probMap FEI3 --- bin x " << it.first.first << " bin y " << it.first.second << " prob " <<
757  it.second);
758  }
759  } else {
760  ATH_MSG_ERROR("Error in printout Charge Coll Prob Maps! (readout should contain FEI3 or FEI4 strings) ");
761  return StatusCode::FAILURE;
762  }
763  return StatusCode::SUCCESS;
764 }
765 
766 // -- Returns the Charge Collection Probability at a given point (bin_x,bin_y)
767 double SensorSim3DTool::getProbMapEntry(const InDetDD::PixelReadoutTechnology &readout, int binx, int biny) const {
768  std::pair<int, int> doublekey(binx, biny);
769  double echarge;
770  if (readout == InDetDD::PixelReadoutTechnology::FEI4) {
771  std::multimap<std::pair<int, int>, double>::const_iterator iter = m_probMapFEI4.find(doublekey);
772  echarge = iter->second;
773  } else if (readout == InDetDD::PixelReadoutTechnology::FEI3) {
774  std::multimap<std::pair<int, int>, double>::const_iterator iter = m_probMapFEI3.find(doublekey);
775  echarge = iter->second;
776  } else {
777  ATH_MSG_ERROR("No Map Entry available for the requested readout");
778  echarge = -1.;
779  }
780  return echarge;
781 }
782 
783 double SensorSim3DTool::getMobility(double electricField, bool isHoleBit) const {
784  //Not exactly the same as the getMobility function in RadDamageUtil, since we don't have a Hall effect for 3D sensors
785  // (B and E are parallel)
786  //Maybe good to do something about this in the future
787 
788  double vsat = 0;
789  double ecrit = 0;
790  double beta = 0;
791 
792  //These parameterizations come from C. Jacoboni et al., Solid-State Electronics 20 (1977) 77-89. (see also
793  // https://cds.cern.ch/record/684187/files/indet-2001-004.pdf).
794 
795  if (isHoleBit) {
796  vsat = 1.62 * std::pow(m_temperature, -0.52); // mm/ns
797  ecrit = 1.24E-7 * std::pow(m_temperature, 1.68); // MV/mm
798  beta = 0.46 * std::pow(m_temperature, 0.17);
799  } else {
800  vsat = 15.3 * std::pow(m_temperature, -0.87); // mm/ns
801  ecrit = 1.01E-7 * std::pow(m_temperature, 1.55); // MV/mm
802  beta = 2.57E-2 * std::pow(m_temperature, 0.66);
803  }
804 
805  double mobility = (vsat / ecrit) / std::pow(1 + std::pow((electricField / ecrit), beta), (1 / beta));
806  return mobility; // mm^2/(MV*ns)
807 }
808 
809 std::vector<double> SensorSim3DTool::getDriftTime(bool isHoleBit, size_t n,
810  CLHEP::HepRandomEngine* rndmEngine,
811  double trappingTimeElectrons,
812  double trappingTimeHoles) const
813 {
814  std::vector<double> rand (n, 0.);
815  std::vector<double> result (n, 0.);
816  CLHEP::RandFlat::shootArray(rndmEngine, n, rand.data(), 0., 1.);
817  for(size_t i = 0; i < n; i++) {
818  if (isHoleBit) {
819  result[i]= (-1.) * trappingTimeHoles * logf(rand[i]); // ns
820  } else {
821  result[i] = (-1.) * trappingTimeElectrons * logf(rand[i]); // ns
822  }
823  }
824  return result;
825 }
PixelID.h
This is an Identifier helper class for the Pixel subdetector. This class is a factory for creating co...
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
createLinkingScheme.iter
iter
Definition: createLinkingScheme.py:62
SiSurfaceCharge
Definition: SiSurfaceCharge.h:23
InDetDD::SolidStateDetectorElementBase::cellIdOfPosition
SiCellId cellIdOfPosition(const Amg::Vector2D &localPos) const
As in previous method but returns SiCellId.
Definition: SolidStateDetectorElementBase.cxx:224
InDetDD::PixelReadoutTechnology
PixelReadoutTechnology
Definition: PixelReadoutDefinitions.h:37
get_generator_info.result
result
Definition: get_generator_info.py:21
InDetDD::DetectorDesign::thickness
double thickness() const
Method which returns thickness of the silicon wafer.
Definition: DetectorDesign.h:271
SG::ReadCondHandle
Definition: ReadCondHandle.h:44
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
PixelHistoConverter::getBinY
std::size_t getBinY(Args &&...args) const
Definition: PixelHistoConverter.h:79
InDetDD::PixelModuleDesign
Definition: PixelModuleDesign.h:45
SensorSim3DTool::m_temperature
Gaudi::Property< double > m_temperature
Definition: SensorSim3DTool.h:97
InDetDD::PixelModuleDesign::columns
int columns() const
Number of cell columns per module:
Definition: PixelModuleDesign.h:331
SensorSim3DTool::m_probMapFEI3
std::multimap< std::pair< int, int >, double > m_probMapFEI3
Definition: SensorSim3DTool.h:70
SensorSimTool::m_radiationDamageSimulationType
Gaudi::Property< int > m_radiationDamageSimulationType
Definition: SensorSimTool.h:87
SiHit.h
SensorSim3DTool::finalize
virtual StatusCode finalize() override
Definition: SensorSim3DTool.cxx:89
SensorSim3DTool::readProbMap
StatusCode readProbMap(const std::string &)
Definition: SensorSim3DTool.cxx:710
SensorSim3DTool::initialize
virtual StatusCode initialize() override
Definition: SensorSim3DTool.cxx:41
min
constexpr double min()
Definition: ap_fixedTest.cxx:26
InDetDD::PixelModuleDesign::rows
int rows() const
Number of cell rows per module:
Definition: PixelModuleDesign.h:336
SiCharge::track
@ track
Definition: SiCharge.h:28
InDetDD::SiCellId::isValid
bool isValid() const
Test if its in a valid state.
Definition: SiCellId.h:136
skel.it
it
Definition: skel.GENtoEVGEN.py:407
SensorSim3DTool::m_cc_prob_file_fei3
Gaudi::Property< std::string > m_cc_prob_file_fei3
Definition: SensorSim3DTool.h:75
InDetDD::DetectorDesign::readoutSide
int readoutSide() const
ReadoutSide.
Definition: DetectorDesign.h:291
InDetDD::PixelReadoutTechnology::FEI3
@ FEI3
SiCharge
Definition: SiCharge.h:25
PixelRadiationDamageFluenceMapData::getTimeMap3D_e
const PixelHistoConverter & getTimeMap3D_e(int layer) const
Definition: PixelRadiationDamageFluenceMapData.cxx:65
TimedHitPtr< SiHit >
SensorSimTool::m_fluenceDataKey
SG::ReadCondHandleKey< PixelRadiationDamageFluenceMapData > m_fluenceDataKey
Definition: SensorSimTool.h:116
PixelRadiationDamageFluenceMapData::getEFieldMap3D
const PixelHistoConverter & getEFieldMap3D(int layer) const
Definition: PixelRadiationDamageFluenceMapData.cxx:60
SensorSimTool::m_chargeCorrectionHistos
Gaudi::Property< std::vector< std::string > > m_chargeCorrectionHistos
Definition: SensorSimTool.h:104
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
SensorSim3DTool::~SensorSim3DTool
virtual ~SensorSim3DTool()
covarianceTool.prob
prob
Definition: covarianceTool.py:678
SensorSim3DTool::m_probMapFEI4
std::multimap< std::pair< int, int >, double > m_probMapFEI4
Definition: SensorSim3DTool.h:69
SensorSim3DTool::getMobility
double getMobility(double electricField, bool isHoleBit) const
Definition: SensorSim3DTool.cxx:783
SensorSimTool::m_siPropertiesTool
ToolHandle< ISiPropertiesTool > m_siPropertiesTool
Definition: SensorSimTool.h:77
SensorSim3DTool::getDriftTime
std::vector< double > getDriftTime(bool isHoleBit, size_t number, CLHEP::HepRandomEngine *rndmEngine, double trappingTimeElectrons, double trappingTimeHoles) const
Definition: SensorSim3DTool.cxx:809
dq_defect_bulk_create_defects.line
line
Definition: dq_defect_bulk_create_defects.py:27
python.CaloAddPedShiftConfig.type
type
Definition: CaloAddPedShiftConfig.py:42
InDetDD::SolidStateDetectorElementBase::identifyHash
virtual IdentifierHash identifyHash() const override final
identifier hash (inline)
SiSurfaceCharge::position
const InDetDD::SiLocalPosition & position() const
Definition: SiSurfaceCharge.h:75
InDetDD::PixelReadoutTechnology::RD53
@ RD53
InDetDD::SiLocalPosition
Definition: SiLocalPosition.h:31
PixelHistoConverter::getBinZ
std::size_t getBinZ(Args &&...args) const
Definition: PixelHistoConverter.h:84
SensorSimTool::initialize
virtual StatusCode initialize()
Definition: SensorSimTool.h:54
PixelRadiationDamageFluenceMapData::getYPositionMap3D_e
const PixelHistoConverter & getYPositionMap3D_e(int layer) const
Definition: PixelRadiationDamageFluenceMapData.cxx:63
SiSurfaceCharge::charge
const SiCharge & charge() const
Definition: SiSurfaceCharge.h:80
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
CheckAppliedSFs.e3
e3
Definition: CheckAppliedSFs.py:264
LArG4FSStartPointFilter.rand
rand
Definition: LArG4FSStartPointFilter.py:80
SiChargedDiodeCollection
Definition: SiChargedDiodeCollection.h:109
PixelHistoConverter
Definition: PixelHistoConverter.h:25
InDetDD::SiDetectorElement::isModuleFrame
bool isModuleFrame() const
Check if the element and module frame are the same.
Definition: SiDetectorElement.cxx:212
CaloCondBlobAlgs_fillNoiseFromASCII.inputFile
string inputFile
Definition: CaloCondBlobAlgs_fillNoiseFromASCII.py:16
PixelRadiationDamageFluenceMapData
Definition: PixelRadiationDamageFluenceMapData.h:25
lumiFormat.i
int i
Definition: lumiFormat.py:85
h
beamspotman.n
n
Definition: beamspotman.py:729
python.SystemOfUnits.micrometer
float micrometer
Definition: SystemOfUnits.py:80
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
extractSporadic.h
list h
Definition: extractSporadic.py:96
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
TauGNNUtils::Variables::Track::dPhi
bool dPhi(const xAOD::TauJet &tau, const xAOD::TauTrack &track, double &out)
Definition: TauGNNUtils.cxx:549
InDetDD::PixelModuleDesign::numberOfCircuits
int numberOfCircuits() const
Total number of circuits:
Definition: PixelModuleDesign.h:306
file
TFile * file
Definition: tile_monitor.h:29
InDetDD::SolidStateDetectorElementBase::width
double width() const
Methods from design (inline)
SensorSim3DTool::SensorSim3DTool
SensorSim3DTool()
test_pyathena.parent
parent
Definition: test_pyathena.py:15
SensorSim3DTool::induceCharge
virtual StatusCode induceCharge(const TimedHitPtr< SiHit > &phit, SiChargedDiodeCollection &chargedDiodes, const InDetDD::SiDetectorElement &Module, const InDetDD::PixelModuleDesign &p_design, std::vector< std::pair< double, double > > &trfHitRecord, std::vector< double > &initialConditions, CLHEP::HepRandomEngine *rndmEngine, const EventContext &ctx) const override
Definition: SensorSim3DTool.cxx:97
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
SensorSim3DTool::printProbMap
StatusCode printProbMap(const std::string &) const
Definition: SensorSim3DTool.cxx:746
InDetDD::SolidStateDetectorElementBase::hitLocalToLocal
Amg::Vector2D hitLocalToLocal(double xEta, double xPhi) const
Simulation/Hit local frame to reconstruction local frame.
Definition: SolidStateDetectorElementBase.cxx:95
SensorSim3DTool::m_radDamageUtil
ToolHandle< RadDamageUtil > m_radDamageUtil
Definition: SensorSim3DTool.h:102
PixelRadiationDamageFluenceMapData::getFluenceLayer3D
double getFluenceLayer3D(int layer) const
Definition: PixelRadiationDamageFluenceMapData.cxx:46
PixelRadiationDamageFluenceMapData::getXPositionMap3D_h
const PixelHistoConverter & getXPositionMap3D_h(int layer) const
Definition: PixelRadiationDamageFluenceMapData.cxx:62
SensorSim3DTool::m_chargeCorrection
std::vector< PixelHistoConverter > m_chargeCorrection
Definition: SensorSim3DTool.h:72
PathResolver.h
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:240
charge
double charge(const T &p)
Definition: AtlasPID.h:986
InDetDD::SiDetectorElement
Definition: SiDetectorElement.h:109
InDetDD::PixelReadoutTechnology::FEI4
@ FEI4
InDetDD::SiDetectorElement::isBarrel
bool isBarrel() const
InDet::SiliconProperties
Definition: SiliconProperties.h:24
PixelRadiationDamageFluenceMapData::getTimeMap3D_h
const PixelHistoConverter & getTimeMap3D_h(int layer) const
Definition: PixelRadiationDamageFluenceMapData.cxx:66
PixelRadiationDamageFluenceMapData::getYPositionMap3D_h
const PixelHistoConverter & getYPositionMap3D_h(int layer) const
Definition: PixelRadiationDamageFluenceMapData.cxx:64
InDetDD::PixelModuleDesign::getReadoutTechnology
PixelReadoutTechnology getReadoutTechnology() const
Definition: PixelModuleDesign.h:391
SiDetectorElement.h
JetVoronoiDiagramHelpers::coord
double coord
Definition: JetVoronoiDiagramHelpers.h:45
SiChargedDiodeCollection::add
void add(const InDetDD::SiCellId &diode, const T &charge)
Definition: SiChargedDiodeCollection.h:299
PathResolverFindCalibFile
std::string PathResolverFindCalibFile(const std::string &logical_file_name)
Definition: PathResolver.cxx:283
InDetDD::SiCellId
Definition: SiCellId.h:29
SensorSimTool::m_digitizeITk3Das3D
Gaudi::Property< bool > m_digitizeITk3Das3D
Definition: SensorSimTool.h:122
PixelRadiationDamageFluenceMapData::getAvgChargeMap3D_h
const PixelHistoConverter & getAvgChargeMap3D_h() const
Definition: PixelRadiationDamageFluenceMapData.cxx:68
InDetDD::PixelModuleDesign::is3D
virtual bool is3D() const
Definition: PixelModuleDesign.h:386
python.CaloAddPedShiftConfig.int
int
Definition: CaloAddPedShiftConfig.py:45
TimedHitPtr::eventId
unsigned short eventId() const
the index of the component event in PileUpEventInfo.
Definition: TimedHitPtr.h:47
InDetDD
Message Stream Member.
Definition: FakeTrackBuilder.h:8
PixelModuleDesign.h
PixelRadiationDamageFluenceMapData::getAvgChargeMap3D_e
const PixelHistoConverter & getAvgChargeMap3D_e() const
Definition: PixelRadiationDamageFluenceMapData.cxx:67
SensorSim3DTool.h
InDet::SiliconProperties::electronHolePairsPerEnergy
double electronHolePairsPerEnergy() const
Definition: SiliconProperties.cxx:238
SensorSim3DTool::m_cc_prob_file_fei4
Gaudi::Property< std::string > m_cc_prob_file_fei4
Definition: SensorSim3DTool.h:81
SensorSimTool
Definition: SensorSimTool.h:38
SensorSim3DTool::m_doChunkCorrection
Gaudi::Property< bool > m_doChunkCorrection
Definition: SensorSim3DTool.h:92
jobOptions.fileName
fileName
Definition: jobOptions.SuperChic_ALP2.py:39
hitTime
float hitTime(const AFP_SIDSimHit &hit)
Definition: AFP_SIDSimHit.h:39
InDetDD::SolidStateDetectorElementBase::length
double length() const
Length in eta direction (z - barrel, r - endcap)
SensorSimTool::m_templateCorrectionRootFile
Gaudi::Property< std::string > m_templateCorrectionRootFile
Definition: SensorSimTool.h:92
MuonParameters::beta
@ beta
Definition: MuonParamDefs.h:144
PixelHistoConverter::getBinX
std::size_t getBinX(Args &&...args) const
Definition: PixelHistoConverter.h:74
InDetDD::SiDetectorElement::design
virtual const SiDetectorDesign & design() const override final
access to the local description (inline):
TauGNNUtils::Variables::Track::dEta
bool dEta(const xAOD::TauJet &tau, const xAOD::TauTrack &track, double &out)
Definition: TauGNNUtils.cxx:538
pow
constexpr int pow(int base, int exp) noexcept
Definition: ap_fixedTest.cxx:15
Amg::distance
float distance(const Amg::Vector3D &p1, const Amg::Vector3D &p2)
calculates the distance between two point in 3D space
Definition: GeoPrimitivesHelpers.h:54
PixelRadiationDamageFluenceMapData::getRamoPotentialMap3D
const PixelHistoConverter & getRamoPotentialMap3D(int layer) const
Definition: PixelRadiationDamageFluenceMapData.cxx:59
SensorSim3DTool::getProbMapEntry
double getProbMapEntry(const InDetDD::PixelReadoutTechnology &, int, int) const
Definition: SensorSim3DTool.cxx:767
InDetDD::SiDetectorElement::isPLR
bool isPLR() const
SiSurfaceCharge.h
PixelRadiationDamageFluenceMapData::getXPositionMap3D_e
const PixelHistoConverter & getXPositionMap3D_e(int layer) const
Definition: PixelRadiationDamageFluenceMapData.cxx:61
SiliconProperties.h
SiChargedDiodeCollection.h
PixelHistoConverter::getContent
float getContent(std::size_t x) const
Definition: PixelHistoConverter.h:34