Modelling a planar silicon sensor
In this example, we will simulate the signal in a 100 μm thick planar silicon sensor (shown schematically in the drawing below) due to the passage of a charged particle.
We will adopt a coordinate system where the back side of the sensor is at y = 0 and the front side (i. e. the strip or pixel side) is at y = 100 μm.
Transport properties
We start by creating a MediumSilicon object, which provides the transport parameters (e. g. drift velocity and diffusion coefficients) of electrons and holes as function of the electric field (and, in general, also the magnetic field, but we will assume that it is zero in this example).
MediumSilicon si;
// Set the temperature [K].
si.SetTemperature(293.15);
Unless the user overrides the default behaviour (by providing a table of velocities at different electric fields),
MediumSilicon calculates the drift velocities according to analytic parameterizations.
A description of the mobility models implemented in MediumSilicon can be found in the user guide.
In this example, we will use the default parameterizations, which correspond to the default models in Sentaurus Device.
The diffusion coefficients are calculated using the Einstein relation.
In order to plot the electron and hole drift velocities as function of the electric field, we can use the class ViewMedium.
ViewMedium mediumView(&si);
mediumView.PlotElectronVelocity('e');
const bool same = true;
mediumView.PlotHoleVelocity('e', same);
Electric field
For accurate calculations of the electric field in a segmented silicon sensor, one normally uses TCAD device simulation programs such as Synopsys Sentaurus Device or Silvaco Atlas.
In the present example, we will follow a simplified approach and approximate the electric field by that of an overdepleted pad sensor. In that case, the x and z components of the electric field vanish, and the y component varies linearly between
at the back side of the sensor (y = 0) and
at the front side of the sensor (y =d), where Vdep is the depletion voltage of the sensor and Vbias is the applied bias voltage.
In this example, we will use Vdep = -20 V and Vbias = -50 V.
In order to use this expression for the electric field in our simulation, we need to write a small function
// Sensor thickness [cm]
constexpr double d = 100.e-4;
// Bias voltage [V]
constexpr double vbias = -50.;
auto eLinear = [](const double /*x*/, const double y, const double /*z*/,
double& ex, double& ey, double& ez) {
// Depletion voltage [V]
constexpr double vdep = -20.;
ex = ez = 0.;
ey = (vbias - vdep) / d + 2 * y * vdep / (d * d);
};
and set up a ComponentUser object which delegates the calculation of the electric field to this function.
ComponentUser linearField;
linearField.SetElectricField(eLinear);
````
We restrict the active volume to a box 0 < *y* < *d*, -2*d* < *x*, *z* < 2*d* and set the medium in this box to silicon.
```cpp
linearField.SetArea(-2 * d, 0., -2 * d, 2 * d, d, 2 * d);
linearField.SetMedium(&si);
A pointer to this ComponentUser object is then passed to a Sensor object which acts as an interface to the transport classes discussed below.
Sensor sensor;
sensor.AddComponent(&linearField);
Weighting field
For signal simulations, we need to know not only the actual electric field in the sensor, but also the weighting field of the electrode for which we want to calculate the induced current (see this tutorial for more details).
In this example, we will use an analytic expression for the weighting field of a strip, as implemented in the class ComponentAnalyticField. We thus create a ComponentAnalyticField object, define the equipotential planes (y = 0 and y = d) and set the voltages at these planes to ground and V = Vbias. We will not use this class to calculate the "real" electric field in the sensor though, so the voltage settings don't actually matter for our purposes.
ComponentAnalyticField wField;
wField.SetMedium(&si);
wField.AddPlaneY(0, vbias, "back");
wField.AddPlaneY(d, 0, "front");
We now define a strip (55 μm width, centred at x = 0) on the front side of the sensor and aptly label it "strip".
constexpr double pitch = 55.e-4;
constexpr double halfpitch = 0.5 * pitch;
wField.AddStripOnPlaneY('z', d, -halfpitch, halfpitch, "strip");
Similarly we could have set up the weighting field of a pixel electrode.
wField.AddPixelOnPlaneY(d, -halfpitch, halfpitch, -halfpitch, halfpitch, "pixel");
ViewField.
ViewField* wfieldView = new ViewField();
wfieldView->SetComponent(&wField);
wfieldView->SetArea(-0.5 * d, 0, 0.5 * d, d);
wfieldView->PlotContourWeightingField("strip", "v");
If we want to plot the x or y component of the weighting field instead of the weighting potential, we just need to replace the second argument in PlotContourWeightingField by "ex" or "ey", respectively.
Finally, we need to instruct the Sensor to use the strip weighting field we just prepared for computing the induced signal
// Request signal calculation for the electrode named "strip",
// using the weighting field provided by the Component object wField.
sensor.AddElectrode(&wField, "strip");
and we need to set the granularity with which we want to record the signal (in our example: 1000 bins between 0 and 10 ns).
// Set the time bins.
const unsigned int nTimeBins = 1000;
const double tmin = 0.;
const double tmax = 10.;
const double tstep = (tmax - tmin) / nTimeBins;
sensor.SetTimeWindow(tmin, tstep, nTimeBins);
Primary ionization and charge carrier transport
We use Heed (see this tutorial) to simulate the electron/hole pairs produced by a 180 GeV / c charged pion traversing the sensor.
TrackHeed track(&sensor);
// Set the particle type and momentum [eV/c].
track.SetParticle("pion");
track.SetMomentum(180.e9);
For transporting the electrons and holes, we use the class AvalancheMC. When setting up the AvalancheMC object, we
need to set the step size used for the drift line calculation to a reasonable value. In this example, we use steps of 1 μm. This means that at each step, the electron/hole will be propagated by 1 μm in the direction of the drift velocity at the local field, followed by a random step based on the diffusion coefficient.
AvalancheMC drift(&sensor);
// Set the step size [cm].
drift.SetDistanceSteps(1.e-4);
We are now ready to run the simulation. In the snippet below, we simulate a perpendicularly incident charged particle track passing through the centre of the strip (x = 0), loop over the electron/hole pairs produced by the particle, and simulate a drift line for each electron and hole.
double x0 = 0., y0 = 0., z0 = 0., t0 = 0.;
double dx = 0., dy = 1., dz = 0.;
track.NewTrack(x0, y0, z0, t0, dx, dy, dz);
// Retrieve the clusters along the track.
for (const auto& cluster : track.GetClusters()) {
// Loop over the electrons in the cluster.
for (const auto& electron : cluster.electrons) {
// Simulate the electron and hole drift lines.
drift.DriftElectron(electron.x, electron.y, electron.z, electron.t);
drift.DriftHole(electron.x, electron.y, electron.z, electron.t);
}
}
To check whether the results are sensible, it can be instructive to visualize the drift lines using the class ViewDrift.
ViewDrift driftView;
driftView.SetArea(-0.5 * d, 0, -0.5 * d, 0.5 * d, d, 0.5 * d);
track.EnablePlotting(&driftView);
drift.EnablePlotting(&driftView);
With the plotting option switched on, AvalancheMC will pass the coordinates of all drift line points to a ViewDrift object. After having simulated all drift lines from a track, we can create a plot using
driftView.Plot2d();
An example of such a drift line plot is shown below (electron drift lines in orange, hole drift lines in red).
For reasons of clarity, only 5% of the drift lines are actually shown in this plot. Plotting the drift lines can slow down the execution time quite a bit, so it is advisable to switch it off when simulating a large number of tracks.
Retrieving the signal
After having simulated the charge carrier drift lines, we can plot the induced current using the class ViewSignal.
ViewSignal signalView(&sensor);
signalView.PlotSignal("strip");
It is possible to post-process the induced current pulse within Garfield++, e. g. by integrating it or by convoluting it with a transfer function that describes the response of the front-end electronics.
Often it can also be useful to save the signal to a file. An example for doing so is given in the snippet of code below.
std::ofstream outfile;
outfile.open("signal.txt", std::ios::out);
for (unsigned int i = 0; i < nTimeBins; ++i) {
const double t = (i + 0.5) * tstep;
const double f = sensor.GetSignal(label, i);
const double fe = sensor.GetElectronSignal(label, i);
const double fh = sensor.GetIonSignal(label, i);
outfile < < t < < " " < < f < < " " < < fe < < " " < < fh < < "\n";
}
Source files
A working program including the different steps discussed in this tutorial is available in Examples/Silicon. The same folder also includes an example for creating an animation of electron/hole drift and signal formation.
