A new Monte Carlo code for full simulation of silicon ...mazziot/nima533_322.pdf · A new Monte...

22
Nuclear Instruments and Methods in Physics Research A 533 (2004) 322–343 A new Monte Carlo code for full simulation of silicon strip detectors M. Brigida, C. Favuzzi, P. Fusco, F. Gargano, N. Giglietto, F. Giordano, F. Loparco, B. Marangelli, M.N. Mazziotta , N. Mirizzi, S. Raino`, P. Spinelli Dipartimento di Fisica ‘‘Michelangelo Merlin’’dell’Universita` and INFN Sezione di Bari, via Amendola 173, I-70126 Bari, Italy Received 1 April 2004; received in revised form 10 May 2004; accepted 17 May 2004 Available online 3 August 2004 Abstract We have developed a full simulation code to evaluate the response of silicon strip detectors (SSDs) to ionizing particles. This simulation can be applied in the design stage of a SSD, when the detector parameters have to be chosen in order to optimize its performance. Our code allows to evaluate the electrical signals produced by ionizing particles crossing a SSD. All the physical processes leading to the generation of electron–hole pairs in silicon have been taken into account. Induced current signals on the readout strips are evaluated using the Shockley–Ramo theorem to the charge carriers propagating inside the detector volume. A simulation of the readout electronics has been implemented, to convert current signals into voltage signals. The predictions of our Monte Carlo simulation have been compared with experimental data taken using a 400 mm thick silicon strip detector with a 228 mm strip pitch exposed to a pion beam. r 2004 Elsevier B.V. All rights reserved. PACS: 07.05.Fb; 29.40.Gx; 29.40.Wk Keywords: Monte Carlo; Full simulation; Silicon strip detectors; Electron–hole pairs; Signal simulation 1. Introduction Silicon strip detectors (SSDs) were first intro- duced as high precision tracking devices for high energy physics experiments in the early 1980s. Since then, impressive technological progress has been made, and SSDs have been used both in accelerator and space experiments. SSDs allow to obtain space resolutions up to a few mm and, thanks to the reduction of their costs, large detector areas can be equipped with these devices. For these reasons, SSDs will be widely adopted in ARTICLE IN PRESS www.elsevier.com/locate/nima 0168-9002/$ - see front matter r 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.nima.2004.05.127 Corresponding author. Tel.: +39-080-5443163; fax: +39- 080-5442470. E-mail address: [email protected] (M.N. Mazziotta).

Transcript of A new Monte Carlo code for full simulation of silicon ...mazziot/nima533_322.pdf · A new Monte...

Page 1: A new Monte Carlo code for full simulation of silicon ...mazziot/nima533_322.pdf · A new Monte Carlo code for full simulation of silicon strip detectors M. Brigida, C. Favuzzi, P.

ARTICLE IN PRESS

0168-9002/$ - se

doi:10.1016/j.ni

�Correspondi

080-5442470.

E-mail addre

Nuclear Instruments and Methods in Physics Research A 533 (2004) 322–343www.elsevier.com/locate/nima

A new Monte Carlo code for full simulation of siliconstrip detectors

M. Brigida, C. Favuzzi, P. Fusco, F. Gargano, N. Giglietto, F. Giordano,F. Loparco, B. Marangelli, M.N. Mazziotta�, N. Mirizzi, S. Raino, P. Spinelli

Dipartimento di Fisica ‘‘Michelangelo Merlin’’dell’Universita and INFN Sezione di Bari, via Amendola 173, I-70126 Bari, Italy

Received 1 April 2004; received in revised form 10 May 2004; accepted 17 May 2004

Available online 3 August 2004

Abstract

We have developed a full simulation code to evaluate the response of silicon strip detectors (SSDs) to ionizing

particles. This simulation can be applied in the design stage of a SSD, when the detector parameters have to be chosen

in order to optimize its performance.

Our code allows to evaluate the electrical signals produced by ionizing particles crossing a SSD. All the physical

processes leading to the generation of electron–hole pairs in silicon have been taken into account. Induced current

signals on the readout strips are evaluated using the Shockley–Ramo theorem to the charge carriers propagating inside

the detector volume. A simulation of the readout electronics has been implemented, to convert current signals into

voltage signals.

The predictions of our Monte Carlo simulation have been compared with experimental data taken using a 400mm

thick silicon strip detector with a 228mm strip pitch exposed to a pion beam.

r 2004 Elsevier B.V. All rights reserved.

PACS: 07.05.Fb; 29.40.Gx; 29.40.Wk

Keywords: Monte Carlo; Full simulation; Silicon strip detectors; Electron–hole pairs; Signal simulation

1. Introduction

Silicon strip detectors (SSDs) were first intro-duced as high precision tracking devices for high

e front matter r 2004 Elsevier B.V. All rights reserve

ma.2004.05.127

ng author. Tel.: +39-080-5443163; fax: +39-

ss: [email protected] (M.N. Mazziotta).

energy physics experiments in the early 1980s.Since then, impressive technological progress hasbeen made, and SSDs have been used both inaccelerator and space experiments. SSDs allow toobtain space resolutions up to a few mm and,thanks to the reduction of their costs, largedetector areas can be equipped with these devices.For these reasons, SSDs will be widely adopted in

d.

Page 2: A new Monte Carlo code for full simulation of silicon ...mazziot/nima533_322.pdf · A new Monte Carlo code for full simulation of silicon strip detectors M. Brigida, C. Favuzzi, P.

ARTICLE IN PRESS

M. Brigida et al. / Nuclear Instruments and Methods in Physics Research A 533 (2004) 322–343 323

the next generation detectors, that will be mountedat the colliders (LHC, RHIC) or that will be flownon satellites.

During the design phase of a detector, there aresome parameters that need to be optimized inorder to ensure the performance required by theexperiment where the detector will be operated. Anaccurate detector model is therefore essential forunderstanding the behaviour of the device. In thelast years, many efforts in the SSD simulation havebeen done (see for instance Refs. [1,2]). We havewritten and developed a full Monte Carlo simula-tion code that simulates all the physical processestaking place in a SSD and allows to calculate theelectrical signals generated by ionizing particles.The details of this simulations will be presented inthis paper, and some results will be discussed.

The output of an electronic detector are thesignals induced on its readout electrodes by themotion of charge carriers (e.g. electron–ion pairsin gaseous detectors, or electron–hole pairs insemiconductor devices) in its sensitive volume.

Charge carriers are generated by the electro-magnetic interaction of ionizing particles inside thesensitive volume of the detector. The ionizationenergy loss D and the number of ionizing collisionsper unit path length, l, can be evaluated from thedifferential cross section sðEÞ, that is related to thekinetic properties of the incident particle and tothe atomic properties of the detector medium.

The energy loss is the result of a number i ofcollisions with energy transfers Ei to knock-onelectrons (d-rays). In a solid medium, e.g. silicon,collective excitations (‘‘plasmons’’) need to betaken into account. Each of these excitationsgenerates a number of charge carriers accordingboth to their energy and to the energy needed toproduce a pair. The fluctuation sn of the totalnumber of pairs n produced in each collision isrelated to the Fano factor F, according to theformula sn ¼

ffiffiffiffiffiffiFn

p[3]. For instance, in silicon the

average energy needed to produce an electron–hole (e–h) pair is about 3:6 eV and the Fano factoris F � 0:1.

In our Monte Carlo code a complete calculationof the production of e–h pairs is done to generatethe charge carriers produced by the absorptionboth of virtual (e.g. exchange photons of a charged

particle with the medium) and real photons (e.g.photons absorbed in the medium) (Sections 2 and3).

Holes generated in a silicon detector will drifttowards the p strips (grounded), electrons towardsthe n electrode (at positive voltage) under theaction of the electric field. During their drift,electrons and holes are diffused by multiplecollisions. Due to the motion of carriers, signalswill be induced on the electrodes. The time shapeof these signals can be evaluated by adding theindividual contributions from all the movingcarriers.

The current induced by a moving carrier on anelectrode can be evaluated by applying theShockley–Ramo theorem, and taking into accountthe drift velocity and the weighting field. In amulti-electrode geometry the weighting fieldsdescribe the geometrical coupling between themoving charge and the electrodes, and areresponsible, together with the diffusion of chargecarriers, of the charge sharing effect (Section 4.1).

The output voltage signals are finally evaluatedby feeding the current signals into the electronicchain, that is simulated as a preamplifier followedby a shaper (Section 4.2).

In our code, the detector noise as well as thenoise associated to the electronics, have been takeninto account. Only ‘‘white noise’’ was considered,and it has been simulated with a sequence of d-likepulses, that are transferred to the electronic chainand are converted into voltage signals (Section4.3).

The simulation code has been written both inFortran and C++ and calculates the ionizationcross-section using the algorithm described in Ref.[4].

2. Energy loss of a charged particle in silicon

A charged particle crossing a silicon layer leavesan ionization track along its trajectory. The mostprobable energy loss Dp is between 180 and300 eV=mm, increasing with the silicon thickness[4]. Since the energy required to produce an e–hpair in silicon is about 3:6 eV, these values

Page 3: A new Monte Carlo code for full simulation of silicon ...mazziot/nima533_322.pdf · A new Monte Carlo code for full simulation of silicon strip detectors M. Brigida, C. Favuzzi, P.

ARTICLE IN PRESS

M. Brigida et al. / Nuclear Instruments and Methods in Physics Research A 533 (2004) 322–343324

correspond to about 50 e–h pairs=mm and 80 e–hpairs=mm respectively.

A complete and accurate calculation of thedifferential cross section sðEÞ as a function ofthe energy loss E in silicon is given in Ref. [4]. Thecollision cross-section sðEÞ is shown in Fig. 1 fortwo bg values, namely for a minimum ionizingparticle (MIP, bg ¼ 5) and for a particle on the‘‘Fermi plateau’’ (bg=2000). The cross-sectionshows three peaks at energies of about 17, 150and 1840 eV, that are associated with the M, L andK shells of silicon.

The differential probability per unit ionizationenergy is given by

F ðEÞdE ¼sðEÞdER Emax

0 sðE0ÞdE0ð1Þ

where Emax is the maximum energy loss. Theprobability to have an energy loss in a singlecollision greater than E is given by

PðEÞ ¼

R Emax

EsðE0ÞdE0

R Emax

0 sðE0ÞdE0: ð2Þ

E (eV)

Na

E2 σ

(E)(

eV/c

m)

103

104

105

106

1 10 102 103 104 105

Fig. 1. Total differential cross section sðEÞ. The function

NaE2sðEÞ is plotted as function of the energy loss E for incident

particles with bg ¼ 5 (solid line) and bg ¼ 2000 (dashed line).

The quantity Na ¼ 4:9938 � 1022 cm�3 is the number of atoms

per unit volume.

Fig. 2 shows the differential probability, F ðEÞ, andthe probability to have an energy loss in a singlecollision greater than E, PðEÞ, for a MIP. Theenergy transfer in a single collision is mainly due tothe contribution of the M shell, while thecontribution from energies above the K shell isless than 1%.

The number l of ionizing collisions per unitpath length is given by

l ¼1R Emax

0 NasðE0ÞdE0ð3Þ

where Na ¼ 4:9938 � 1022 cm�3 is the density ofsilicon atoms. Fig. 3 shows the number ofcollisions per mm as a function of the bg factorof the charged particle. The number of collisionsper mm is about 3.8 for values of bg above the MIPvalue.

A Monte Carlo method is used to calculate thestraggling functions that describe the distributionsof energy loss for different silicon thicknesses. Theenergy loss Ei in the ith collision is determined byexctracting a random number from the differentialprobability distribution F ðEÞ. Individual energylosses are added up to give the total energy loss ofthe particle. The distance traveled between

E (eV)

F(E

), P

(E)

10-7

10-6

10-5

10-4

10-3

10-2

10-1

1

1 10 102 103 104 105

Fig. 2. Differential (solid line) and integral (dashed line) energy

loss probabilities as a function of the energy loss.

Page 4: A new Monte Carlo code for full simulation of silicon ...mazziot/nima533_322.pdf · A new Monte Carlo code for full simulation of silicon strip detectors M. Brigida, C. Favuzzi, P.

ARTICLE IN PRESSN

umbe

r of

col

lisio

n pe

r µm

3

3.5

4

4.5

5

5.5

6

1 10 102 103

βγ

Fig. 3. Number of collisions per mm as a function of bg.

Energy loss (keV)N

umbe

r of

eve

nts/

2.5

keV

0

25

50

75

100

125

150

175

200

225

0 200 400 600 800 1000

Fig. 4. Energy loss distribution for bg ¼ 5 electrons crossing a

400mm thick silicon detector.

p(MeV/c)

Mos

t pr

obab

le e

nerg

y lo

ss (

keV

)

100

105

110

115

120

125

130

135

140

1 10 102 103 104 105 106

Fig. 5. Most probable energy loss as a function of the

momentum for electrons (circles), muons (triangles), pions

(stars) and protons (squares), passing through a silicon detector

of 400mm thickness.

M. Brigida et al. / Nuclear Instruments and Methods in Physics Research A 533 (2004) 322–343 325

collisions is sorted according to an exponentialcollision probability, i.e.

PðxÞ ¼expð�x=lÞ

l: ð4Þ

The calculation is repeated for N particles and theenergy loss distribution is thus generated.

Fig. 4 shows the energy loss distributionobtained by simulating a sample of 4000 electronswith bg ¼ 5 crossing a silicon detector of 400 mmthickness. The most probable value is about107 keV, while the mean value is about 139 keVand the full width at half maximum (FWHM) isabout 43 keV. The energy loss distribution exhibitsa large tail up to about 1MeV.

Fig. 5 shows the most probable energy loss as afunction of the particle momentum for electrons,muons, pions and protons passing through a400mm thick silicon detector. The most probablevalue of the energy loss of particles on the Fermiplateau is about 7% greater than the MIP value.

Energy loss of a charged particle will causeionization in the silicon, both in the primarycollision and by the subsequent energy losses ofsecondary electrons. Usually, the energy lost by acharged particle traversing an absorber layer is

deposited along its trajectory. When a large energyloss occurs, some d-rays can have enough energyto escape from the detector volume. In this case,

Page 5: A new Monte Carlo code for full simulation of silicon ...mazziot/nima533_322.pdf · A new Monte Carlo code for full simulation of silicon strip detectors M. Brigida, C. Favuzzi, P.

ARTICLE IN PRESS

Initial electron kinetic energy (keV)R

ange

(µm

)

10-1

1

10

102

0 20 40 60 80 100 120 140 160 180 200

Fig. 6. Range of d-electrons as function of their initial kinetic

energy in silicon. Open circles: ranges evaluated without taking

into account the multiple scattering effect; black circles: ranges

evaluated taking into account the multiple scattering effect

(gaussian approximation). The solid line shows the best fit of

the data without multiple scattering, the dashed line shows the

best fit of the data with multiple scattering. For comparison, the

dotted line shows the parameterization reported in Ref. [5].

M. Brigida et al. / Nuclear Instruments and Methods in Physics Research A 533 (2004) 322–343326

the energy deposited in the absorber is less thanthe energy lost by the particle.

A d electron travels a distance in the absorberthat depends on its energy. The path length of a d-ray until it is stopped, the so called ‘‘practical

range’’, is randomized both by the multiplescattering in the medium and by the energy lostin the collisions that slow down the electron. Wehave calculated the range of d-rays in silicon bymeans of a Monte Carlo method, where thedistances traveled between collisions are evaluatedas discussed above. In each collision the electronkinetic energy is decreased by the energy loss andthe direction is randomized by taking into accountthe multiple scattering angle in the gaussianapproximation. In the next step the differentialcross-section sðEÞ and the number of ionizingcollisions per unit path length l are evaluatedstarting from the residual kinetic energy. Thecalculation is repeated for each step until theresidual energy falls under a threshold value thathas been set to 1 keV.

Fig. 6 shows the electron range calculated by theMonte Carlo method as a function of its initialkinetic energy in the range 10–100 keV. The MonteCarlo data can be fitted by the function [5]:

RðEÞ ¼ AE 1 �B

1 þ CE

� �: ð5Þ

The values of the parameters areA ¼ 2:9567mmkeV�1, B ¼ 0:9969 and C ¼

5:243 � 10�3 keV�1 when the multiple scatteringis not taken into account, while result to beA ¼ 1:9841mmkeV�1, B ¼ 0:9957 and C ¼

4:603 � 10�3 keV�1 when the multiple scatteringis taken into account. The parameters previouslyreported by [5] are A ¼ 2:305mmkeV�1, B ¼

0:9815 and C ¼ 3:123 � 10�3 keV�1.The practical range of a 10 keV electron in

silicon is about 1mm. It becomes comparable withthe typical detectors thicknesses of few a 100mmwhen the electron energy is larger than 100 keV.Since the probability of ionization energy losseslarger than 100 keV in a single collision is less than10�3 (see Fig. 2), we can expect that the number ofescaping d-rays will be negligible. However, the d-electron practical range becomes of the same orderof the typical SSD strip pitch (10–100mm) for

energies of few 10 keV. These d-electrons can losetheir energy in a region far away from the path ofthe incident particle, inducing current signals inseveral read-out strips.

As it will be discussed in Section 3.4, ionizationenergy losses are the starting point for thegeneration of charge carriers inside silicon. Acomplete description of the whole ionizationprocess allows a precise knowledge of the positionswhere charge carriers are generated. As it will beshown in Sections 4 and 5, this point is crucial forunderstanding the charge sharing effect occurringin SSDs.

3. Generation of electron–hole pairs in silicon

In this section we present our approach toevaluate the distribution of e–h pairs produced byionizing charged particles in silicon. The chargecarriers are produced both in the primary collision

Page 6: A new Monte Carlo code for full simulation of silicon ...mazziot/nima533_322.pdf · A new Monte Carlo code for full simulation of silicon strip detectors M. Brigida, C. Favuzzi, P.

ARTICLE IN PRESS

Valence band

Photon Energy (eV)P

roba

biit

iles

of

phot

oele

ctri

c ab

sorp

tion

0

0.2

0.4

0.6

0.8

1

1 10 102 103 104

Fig. 7. Relative probabilities of photoelectric absorption by the

silicon shells as a function of the photon energy. Solid line: M

shell absorption; dashed line: L2 and L3 shell absorptions;

dotted line: L1 shell absorption; dashed dotted line: K shell

absorption. The zero energy level corresponds to the top of the

valence energy band.

M. Brigida et al. / Nuclear Instruments and Methods in Physics Research A 533 (2004) 322–343 327

and by the subsequent energy losses of thesecondary e–h pairs. Usually, the energy loss in aprimary collision is larger than the energy W

needed to produce an e–h pair (� 3:6 eV). Thus,for each collision of the incident particle, severale–h pairs will be produced.

The primary collision occur with the exchangeof ‘‘virtual’’ photons between the incident chargedparticle and the medium. The absorption of thesephotons and the subsequent relaxation of theexcited atoms yield the photoelectrons, the Augerand Coster-Kronig ‘‘primary electrons’’ and thecorresponding ‘‘primary holes’’ in the valenceband. All primary electrons and holes lose partof their kinetic energy by phonon scattering andanother part by impact ionizations, producingsecondary e–h pairs [6].

3.1. Primary electron–hole pairs

When photons are absorbed by an inner shell x

(x ¼ K ;L1;L2;3), the photoelectron energy is de-termined by the photon energy, whereas theenergies of primary electrons and primary holesdepend only on the vacancies created in thevalence band by the photoionization process.The photoelectron energy, Epe, is given by Epe ¼

Eg � Ex � Egap where Egap ¼ 1:12 eV is the energyband gap at the temperature T ¼ 300 K. In case ofphotoabsorption by electrons in the valence band(M shell), the energies of the primary vacancies,and thus the energy of the photoelectrons, are notsharply defined, as the valence band has a finitewidth of 12 eV. Hence, the energies of electronsand holes originated in the valence band aredistributed according to a continuous energyspectrum.

Fig. 7 shows the relative probabilities of photo-electric absorption by the various silicon electronshells as a function of the photon energy. Thetransition rates are taken by Ref. [7]. For energiesabove the K edge (EK ¼ 1839 keV) there is aconstant 92% probability of absorption by the K

shell and an 8% probability of absorption by theL1 shell. In the energy region between the L1 and K

edges the absorption is distributed in an energy-dependent way between the L1 shell(EL1

¼ 148:7 eV), the L2 and the L3 shells (that

are assumed to be degenerate at EL2;3¼ 99:2 eV)

and the M shell (EM ¼ EV ¼ 12 eV). The M-shelledge stands as the bottom of valence band, thusthe zero energy level corresponds to the top of thevalence band.

The relaxation process following photoabsorp-tion in the atomic subshells produces electrons andvacancies in the K, L1 and L2;3 shells according tothe relaxation trees and probabilities shown inTables 1–3, that have been taken from Ref. [7]. Weignored the small probability of L2;3 shell fluores-cence and did not account for possible reabsorp-tion of K shell fluorescence photons inside thedetector volume. In other words, for the K and L2;3

shells relaxations, a 100% probability of Augeremission is assumed. The relaxation process endswhen the vacancies are transferred to the valenceband (M shell).

The primary hole generated from the vacanciesKLM and LLM is assigned a random energy Eh

extracted with a uniform probability in therange 0pEhpEV . The energy of the Auger and

Page 7: A new Monte Carlo code for full simulation of silicon ...mazziot/nima533_322.pdf · A new Monte Carlo code for full simulation of silicon strip detectors M. Brigida, C. Favuzzi, P.

ARTICLE IN PRESS

Table 1

Relaxation tree for the K shell

K-shell vacancy and photoelectron

Epe ¼ Eg � EK � Egap ðEK ¼ 1839 eVÞ

Auger emissions K-shell fluorescence

(95.6%) (4.4%)

Transition Probability Vacancies Probability Photon energy

chain (%) (%) and vacancy

KL1L1 19.2 2 L1 vacancies 59.3 1740 eV photon

L3 vacancy

KL1L2;3 38.9 L1, L2;3 vacancies 29.6 1740 eV photon

L2 vacancy

KL2;3L2;3 23.3 2 L2;3 vacancies 11.1 1836 eV photon

M vacancy

KL1M 7.5 L1, M vacancies

KL2;3M 10.4 L2;3, M vacancies

KMM 0.8 2 M vacancies

M. Brigida et al. / Nuclear Instruments and Methods in Physics Research A 533 (2004) 322–343328

Coster-Kroning primary electrons is defined asEe ¼ Ex � Ey � Eh � Egap where Ex � Ey is thedifference between the energies of the inner shells x

and y involved in the transition chain (e.g. EK �

EL for the KLM chain and EL1� EL2;3

for theLLM one).

For xMM Auger electrons, with x ¼ K ;L1;L2;3,the electron energy is randomly generated as [6]Ee ¼ Ex � ð1 þ rÞEV � Egap where r is randomlydistributed in the range �1prp1 with probabilitypðrÞ ¼ 1 � jrj, and Ex is the Auger energy, that isequal to the edge energy involved in the transitionchain, e.g. Ex ¼ EK and Ex ¼ EL for the transitionchains KMM and LMM, respectively. The energyof the first primary hole Eh1 corresponding to theAuger electron is then generated with a uniformprobability in the range 0pEh1pð1 þ rÞEV ; hence,the energy of the second primary hole is given byEh2 ¼ ð1 þ rÞEV � Eh1. Electron and the holeenergies are such that Ee þ Eh1 þ Eh2 ¼ Ex � Egap.

In the region between M and L2;3 edges, theprimary hole energy Eh is randomly chosen with auniform probability in the range 0pEhpEV andthe photoelectron is assigned an energyEpe ¼ Eg � Egap � Eh.

For the absorption of photons having energiesbetween Egap and EV , the primary hole energy

is extracted with a uniform probability in therange 0pEhpðEg � EgapÞ, and the photoelectronis thus assigned the residual energy Epe ¼

Eg � Egap � Eh.Finally, photons having energies less than Egap

are not absorbed: in other words, no latticedefects, mainly impurities, have been taken intoaccount.

3.2. Secondary electron–hole pairs

A primary electron (hole) with energy E largecompared to the semiconductor band gap caninteract both by ionization, i.e. creating secondarye–h pairs, and by phonon emission. The scattering-rate assumption from Ref. [8] is used to generatethe cascades of e–h pairs and phonons. Theprobability P0ðEÞ of the first interaction of acarrier of energy E being an ionization process(1 � P0ðEÞ is the probability that the first interac-tion is a phonon emission), is given by

P0ðEÞ ¼1

1 þ ðr0ðEÞ=rðEÞÞð6Þ

where rðEÞ and r0ðEÞ are the scattering rates byionization and by phonon emission respectively.

Page 8: A new Monte Carlo code for full simulation of silicon ...mazziot/nima533_322.pdf · A new Monte Carlo code for full simulation of silicon strip detectors M. Brigida, C. Favuzzi, P.

ARTICLE IN PRESS

M. Brigida et al. / Nuclear Instruments and Methods in Physics Research A 533 (2004) 322–343 329

The ratio between these rates is given by

r0ðEÞ

rðEÞ¼ A

105

2pðE � E0Þ

1=2

ðE � EgapÞ7=2

ð7Þ

where E0 ¼ 63 meV is the phonon energy at 300K,and A ¼ 5:2 eV3 is a constant [8]. The scatteringrate for the electron–phonon interaction due to thedeformation potential and to the polar-modeelectrostatic potential are not taken into account,since in silicon their contribution is generally small[8].

The generation of secondary e–h pairs is acascade process, that is simulated with a MonteCarlo method. Each particle with energy E isfollowed up from one step to the next. At the endof each step a random number R in the range0pRp1 is generated with a uniform probability;the scattering step will end with an ionizationprocess if RpP0ðEÞ, or with a phonon emission ifR4P0ðEÞ. If the scattering step ends with aphonon emission, the scattered particle is assignedthe energy E � E0. Otherwise, if the scattering stepends with an ionization, two new particles (theelectron and hole produced by ionization) arecreated, with kinetic energies E1 and E2. In thiscase, the scattered particle is assigned the energyE � ðE1 þ E2 þ EgapÞ. The energies E1 and E2 arerandomly generated according to the probabilitiesshown in Ref. [9] in the free-particle approxima-tion, i.e. E1 ¼ a1ðE � EgapÞ and E2 ¼

a2ðE � Egap � E1Þ, where a1 and a2 are randomnumbers in the range ½0; 1� extracted from theprobability density functions f ða1Þ / ð1 � a1Þ

2 ffiffiffiffiffia1

p

and f ða2Þ /ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffia2ð1 � a2Þ

p, respectively.

In the next step the particle is assigned itsresidual kinetic energy. The calculation is repeatedfor each step until the residual energy falls belowthe threshold energy Eth ¼ 3Egap=2, that is deter-mined [8] by the momentum conservation.1 Allsecondary electron–hole particles produced by theprimary particle are followed up with the sameprocedure until all particles have residual energiesbelow the threshold. These particles are the total

1No appreciable changes are observed if the threshold energy

is assumed to be Eth ¼ Egap, namely the value obtained by

neglecting momentum conservation.

charge carriers, nc, produced by the primaryparticle.

The total number of secondary e–h pairs is givenby np ¼ ðnc � 1Þ=2, where ‘‘1’’ stands for theprimary particle. The pair creation energy W ðEÞ

and the Fano factor F ðEÞ are given by

W ðEÞ ¼E

hnpðEÞið8Þ

F ðEÞ ¼hn2

pðEÞi � hnpðEÞi2

hnpðEÞið9Þ

where

hnpðEÞi ¼1

N

XN

i¼1

np;i ð10Þ

hn2pðEÞi ¼

1

N

XN

i¼1

n2p;i: ð11Þ

In the previous formulas N is the number ofsimulated events for each energy, and np;i is thetotal number of secondary e–h pairs in the ithevent.

Figs. 8–10 show the average number of second-ary e–h pairs, the mean energy needed to producea pair and the Fano factor as function of theprimary particle kinetic energy. The energy valuesare chosen in order to have 30 points in the ranges½Egap;EV � and ½EV ;EL2;3

�, and 10 points in each ofthe ranges ½EL2;3

;EL1�, ½EL1

;EK � and ½EK ; 10 keV�.The points are equally spaced on a logarithmicscale.

The mean energy needed to produce an e–h pairfor large primary energies (EX100 eV) approachesa constant value W1 � 3:65 eV. The value of W

increases with decreasing primary kinetic energiesfrom W1 to a maximum at E � 6 eV. A furtherdecrease of E shows that W reaches a minimum inthe region around 4 eV. In this region, the particleenergy becomes comparable with the pair creationenergy and np approaches 1 (see Fig. 8). The Fanofactor for the pair number distribution is almostconstant (F ’ 0:117): the fluctuations around thisconstant value in Fig. 10 are due to statisticalerrors associated to the number of simulatedevents.

Page 9: A new Monte Carlo code for full simulation of silicon ...mazziot/nima533_322.pdf · A new Monte Carlo code for full simulation of silicon strip detectors M. Brigida, C. Favuzzi, P.

ARTICLE IN PRESS

Electron (hole) energy (eV)

Fan

o fa

ctor

0

0.05

0.1

0.15

0.2

0.25

0.3

1 10 102 103 104

Fig. 10. Fano factor, F, as function of the primary electron

(hole) kinetic energy. For each energy bin 2000 events have

been simulated. The points are connected by a dashed line. The

dotted line shows the asymptotic value F1 ¼ 0:117 for large

primary kinetic energy.

Electron (hole) energy (eV)

Ave

rage

num

ber

of s

econ

dary

pai

rs

10-1

1

10

102

103

104

1 10 102 103 104

Fig. 8. Average number of secondary electron–hole pairs as

function of the primary electron (hole) kinetic energy. For each

energy bin 2000 events have been simulated. The dashed line

connecting the points is drawn to guide the eye.

Electron (hole) energy (eV)

Pai

r cr

eati

on e

nerg

y (e

V)

2

3

4

5

6

7

8

9

10

1 10 102 103 104

Fig. 9. Electron–hole pair creation energy, W, as a function of

the primary electron (hole) kinetic energy. For each energy bin

2000 events have been simulated. The points are connected by a

dashed line. The dotted line (W1 ¼ 3:645 eV) is the asymptotic

value for large primary kinetic energy.

M. Brigida et al. / Nuclear Instruments and Methods in Physics Research A 533 (2004) 322–343330

3.3. Electron–hole pairs produced by photons

In this section we show our Monte Carlo resultsfor the secondary e–h pairs produced by photons.The absorption of a photon and the followingrelaxation process yields the primary e–h pairs, asdiscussed in Section 3.1, that lose then their kineticenergy producing secondary e–h pairs as discussedin Section 3.2.

Figs. 11–13 show the average number ofsecondary e–h pairs and carriers, the mean energyto produce a pair and the Fano factor as afunction of the photon energy. The energy valueshave been chosen in order to have 30 points in theranges ½Egap;EV � and ½EV ;EL2;3

�, and 10 points ineach of the ranges ½EL2;3

;EL1�, ½EL1

;EK � and½EK ; 30 keV�. The points are equally spaced on alogarithmic scale. For each energy bin 2000 eventshave been simulated.

For each photon event, the number of second-ary e–h pairs is given by np ¼ ðnc � 2Þ=2, where‘‘�2’’ accounts for the primary electron and hole.The pair creation energy W and the Fano factor F

Page 10: A new Monte Carlo code for full simulation of silicon ...mazziot/nima533_322.pdf · A new Monte Carlo code for full simulation of silicon strip detectors M. Brigida, C. Favuzzi, P.

ARTICLE IN PRESS

Photon energy (eV)

Ave

rage

num

ber

of s

econ

dary

pai

rs (

carr

iers

)

10-1

1

10

102

103

104

1 10 102

103

104

Fig. 11. Average number of secondary e–h pairs, np (black

circles) and carriers, nc (open circles) as a function of the

photon energy. For each energy bin 2000 events have been

simulated. The points are connected by a dashed line. The inset

shows an expanded view of the behaviour of np in the low

energy region. For comparison, it is also shown the behaviour

of np assuming that the total energy is transferred to the

photoelectron (squares), or is equally split between the electron

and hole (stars).

Photon energy (eV)P

air

crea

tion

ene

rgy

(eV

)2

3

4

5

6

7

8

9

10

1 10 102

103

104

2

3

4

5

6

7

8

9

10

1 10 10

Fig. 12. Electron–hole pair creation energy, W, as function of

the photon energy. For each energy bin 2000 events have been

simulated. The points are connected by a dashed line. The

dotted line shows the asymptotic value W1 ¼ 3:645 eV. The

inset shows an expanded view of the low energy region. For

comparison, it is also shown the behaviour of W assuming that

the total energy is transferred to the photoelectron (squares), or

is equally split between the electron and hole (stars).

Photon energy (eV)

Fan

o fa

ctor

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

1 10 102

103

104

Fig. 13. Fano factor, F, as function of the photon energy. For

each energy bin 2000 events have been simulated. The points

are connected by a dashed line. The dotted line shows the

asymptotic value F1 ¼ 0:117 for large primary kinetic energy.

M. Brigida et al. / Nuclear Instruments and Methods in Physics Research A 533 (2004) 322–343 331

are calculated as described in Section 3.2. Ofcourse, the average number of carriers is twice theaverage number of pairs for large photon energies.

The behaviour of the e–h pair creation energy W

(Fig. 12) is similar to that shown in Fig. 9. Forenergies larger than the L1 edge, it approaches aconstant value W1 � 3:65 eV. The structure thatappears in the region of the L edges (around100 eV), is due to the many processes involved forthe photoelectric absorption, as shown in Fig. 7and Tables 2–3, and to the random energies of theholes produced in the valence band. The e–h pairproduction energy increases as the energy de-creases in the same way seen in Section 3.2, butwith a smoother behaviour due to the influence ofthe distribution of the photoelectrons and primaryholes. The average number of carriers nc shown inFig. 11 can be used to evaluate the chargecollection efficiency in a silicon device, aftertaking into account the linear photon absorption

Page 11: A new Monte Carlo code for full simulation of silicon ...mazziot/nima533_322.pdf · A new Monte Carlo code for full simulation of silicon strip detectors M. Brigida, C. Favuzzi, P.

ARTICLE IN PRESS

Table 2

Relaxation tree for L1 shell

L1-shell vacancy and photoelectron

Epe ¼ Eg � EL1� Egap ðEL1

¼ 148:7 eVÞ

Transition

chain

Probability

(%)

Vacancies Emission

L1MM 2.5 2 M vacancies Auger

L1L2;3M 97.5 L2;3, M vacancies Coster-Kroning

Table 3

Relaxation tree for L2;3 shell

L2;3-shell vacancy and photoelectron

Epe ¼ Eg � EL2;3� Egap (EL2;3

= 99.2 eV)

Transition

chain

Probability

(%)

Vacancies Emission

L2;3MM 99.9 2 M vacancies Auger

0.1 Fluorescence

M. Brigida et al. / Nuclear Instruments and Methods in Physics Research A 533 (2004) 322–343332

coefficient. The photoelectron and the Augerelectron range should be taken into account toevaluate the carrier production along the pathlength.

The results obtained with our simulation are inagreement with the data shown in Refs. [6,10]. Inour code, the band structure of silicon is notconsidered to determine the initial energy ofphotoelectron for low photon energies. In otherwords, we assume an uniform density of states inthe valence band instead of a more realisticdistribution [11]. This approximation does notappreciably affect the behaviour of the paircreation energy, except in the low photon energyregion, e.g. Egp10 eV. To check this simpleapproximation, we calculated the number ofsecondary e–h pairs in two extreme cases, i.e.when the total energy is transferred to thephotoelectron, and when for EphpEV the totalenergy is equally shared between the photoelectronand the hole or, for higher photon energies, themaximum energy EV is transferred to the hole andthe photolectron energy is Epe ¼ Eg � Egap � EV .

These results are shown in the insets of Figs. 11and 12, respectively.

3.4. Charge carriers produced by ionizing particles

In this section we present the simulation resultsfor the carriers produced by ionizing chargedparticles in silicon. The generation of chargecarriers can be divided in three steps:

(1)

The energy loss in each collision is evaluated asdiscussed in Section 2. The energy loss is due tothe interaction of the charged particle with thesilicon by the exchange of ‘‘virtual’’ photons.

(2)

The absorption of the virtual photons and thefollowing relaxation processes yield the pri-mary e–h pairs as discussed in Section 3.1. Thevirtual photons are absorbed close to theparticle path, because the ionization cross-section is related to an impact parameter� 1 (A.

(3)

All primary electrons and holes lose theirkinetic energy producing secondary e–h pairsas discussed in Section 3.2.

Since the practical range for electrons withEp10 keV is less than 1 mm (Fig. 6), carriersproduced by these electrons are located along theionizing particle path in the same points wherecollisions occurred. On the other hand, primaryelectrons with E410 keV (i.e. photoelectrons) areejected orthogonally with respect to the incidentparticle direction and are tracked taking intoaccount the practical range and the multiplescattering effect. These electrons, as well as thecharge carriers that they eventually generate, canescape from the path of the ionizing particle andsometimes can even escape from the detectorvolume.

Fig. 14 shows the distribution of the e–h pairsproduced by bg ¼ 5 electrons passing through asilicon detector of 400mm thickness. The numberof pairs is calculated by dividing by two thenumber of carriers. The most probable value in400mm silicon is np � 2:9 � 104, the average valueis np � 3:6 � 104 and the FWHM is about 1 � 104

pairs. The large tail of the energy loss distributionshown in Fig. 4 partially disappeared in the e–h

Page 12: A new Monte Carlo code for full simulation of silicon ...mazziot/nima533_322.pdf · A new Monte Carlo code for full simulation of silicon strip detectors M. Brigida, C. Favuzzi, P.

ARTICLE IN PRESS

Number of e-h pairs (×103 )

Num

ber

of e

vent

s/25

00

0

100

200

300

400

500

600

700

800

0 25 50 75 100 125 150 175 200

Fig. 14. Distribution of the number of pairs produced by bg ¼5 electrons passing through a silicon detector of 400mm

thickness. The number of simulated events is 4000.

Pairs along the particle direction, X(µm)

Pai

rs a

long

the

ort

hogo

nal p

arti

cle

dire

ctio

n, Y

(µm

)

-200

-150

-100

-50

0

50

100

150

200

-200 -100 0 100 200 300 400 500 600

Fig. 15. Space distribution of pairs produced by an electron

with bg ¼ 5 in a 400mm silicon layer. The dotted lines show the

detector boundaries. The particle direction is the X-axis.

M. Brigida et al. / Nuclear Instruments and Methods in Physics Research A 533 (2004) 322–343 333

pair distribution, due to the escape effect for theprimary electrons of higher energies.

Fig. 15 shows the positions of e–h pairs withrespect to the particle direction. There are many

pairs located far away from the particle trajectory,and sometime even outside the detector layer.

4. Signal simulation

In this section we describe our approach toevaluate the electrical output of a SSD. Thecurrents induced by moving carriers on the read-out electrodes are evaluated by means of theShockley–Ramo theorem [12–16] and the outputvoltage signals are evaluated by feeding the currentsignals into the electronic chain. Finally, the noiseis simulated by generating d-like pulses.

4.1. Generation of current signals

Silicon detectors can be described as p–njunction diodes, operated with a reverse biasvoltage that forms a sensitive region depleted ofmobile charge carriers, in order to allow thecollection of free charge carriers generated byradiation. SSDs are usually built with an asym-metric structure, consisting of highly doped p (pþ)strips implanted into a lightly doped n bulk, sothat the depletion region extends mainly into thelightly doped volume. A scheme of the typicalgeometry of a SSD is shown in Fig. 16, where d isthe detector thickness, p is the strip pitch, w is thestrip width and h is the depth of the pþ implant.

In order to simulate the motion of chargecarriers inside the detector, it is necessary toevaluate the electric potential and consequentlythe electric field inside the detector. This has beendone by dividing the detector volume in elemen-tary cells, according to the scheme shown inFig.16, and solving in a single cell the inhomoge-neous Maxwell equation:

~r � ~D ¼ r ð12Þ

with the following boundary conditions on theelectric potential:

fðx ¼ 0Þ ¼ f0 ð13Þ

fðx ¼ d;�w=2pypw=2Þ ¼ 0 ð14Þ

fðy ¼ �p=2Þ ¼ fðy ¼ þp=2Þ ð15Þ

Page 13: A new Monte Carlo code for full simulation of silicon ...mazziot/nima533_322.pdf · A new Monte Carlo code for full simulation of silicon strip detectors M. Brigida, C. Favuzzi, P.

ARTICLE IN PRESS

y

p+

p+

p+

p+

p+

n bulkn+

p w

h

d

Fig. 16. Schematic layout of a SSD: d is the detector thickness,

p is the strip pitch, w is the strip width and h is the depth of the

pþ implant. The origin of the local reference system is such that

the coordinates of the centers of the pþ strips are x ¼ d and

y ¼ kp, where k ¼ 0, �1, �2, . . . . The boundaries of

elementary detector cells have been represented with dashed

lines.

x (µm)y (µm)

V (

Vol

t)

050

100150

200250

300350

400

-100-75

-50-25

025

5075

100

0

20

40

60

80

100

Fig. 17. Electric potential in an elementary cell of a SSD. The

values of the geometrical parameters are: d ¼ 400mm,

p ¼ 228mm, w ¼ 60mm, h ¼ 5mm.

M. Brigida et al. / Nuclear Instruments and Methods in Physics Research A 533 (2004) 322–343334

The boundary condition (13) implies that the n

bulk is connected to a back plate kept at a voltagef0; condition (14) means that the pþ strips aregrounded (i.e. the input impedance of the pre-amplifiers is negligible). Condition (15) is requiredby the geometrical periodicity of the SSD alongthe y direction (Fig. 16). The electric field has beenevaluated on a two-dimensional discrete mesh byusing the MAXWELL package [18]. Fig. 17 showsthe electric potential fðx; yÞ evaluated for a SSDwith d ¼ 400 mm, p ¼ 228mm, w ¼ 60mm, h ¼

5mm and f0 ¼ 100V. A uniform doping has beenassumed for both the pþ strips and the n bulk.

After their production, under the action of theelectric field, holes and electrons will drift respec-tively towards the pþ strips and the n electrode.Due to the motion of carriers, current signals willbe induced on the electrodes. The time shape ofthese signals can be evaluated by adding theindividual contributions from all the movingcarriers.

The motion of a carrier (an electron or a hole)can be described by the equation

~v ¼ m~E ð16Þ

where ~v is the carrier velocity, ~E is the electric fieldand m is the mobility, that depends on the electricfield, i.e. m ¼ mðEÞ. For the mobilities of electronsand holes as a function of the electric field E

applied along a h1 1 1i direction, we assume thefollowing parameterization [17]:

m ¼vm=Ec

ð1 þ ðE=EcÞbÞ1=b

ð17Þ

where the values and temperature dependences ofthe parameters vm, Ec and b are shown in Table 4.

The first-order differential equation of motion(16) is solved numerically, using the Runge–Kuttamethod, in which the time step size dt is adjustedso that the overall space accuracy increases, whilekeeping the amount of CPU time used reasonablysmall. The value of the drift velocity ~v½~rðtÞ� at thebeginning of each time step is used to obtain thesize of the step according to the formula:

dt ¼�

j~v½~rðtÞ�jð18Þ

Page 14: A new Monte Carlo code for full simulation of silicon ...mazziot/nima533_322.pdf · A new Monte Carlo code for full simulation of silicon strip detectors M. Brigida, C. Favuzzi, P.

ARTICLE IN PRESS

Table 4

Parameters for the electric field and temperature dependence

(TðKÞ) of electron and hole mobilities in high-purity silicon

(TX250K)

Parameter Electrons Holes Units

vm 1:53 � 109 � T�0:87 1:62 � 108 � T�0:52 cm2=ðVsÞ

Ec 1:01 � T1:55 1:24 � T1:68 V/cm

b 2:57 � 10�2 � T0:66 0:46 � T0:17 ���

x (µm)y

(µm

)-150

-100

-50

0

50

100

150

0 50 100 150 200 250 300 350 400

Fig. 18. Drift line in a silicon cell. The geometrical parameters

are d ¼ 400mm, p ¼ 228mm, w ¼ 60mm, h ¼ 5mm. The bias

potential is f0 ¼ 100V.

M. Brigida et al. / Nuclear Instruments and Methods in Physics Research A 533 (2004) 322–343 335

where � is the absolute integration accuracy and~rðtÞ is the position at the time t. The drift lines ofcharge carriers in a silicon cell with d ¼ 400mm,p ¼ 228mm, w ¼ 60 mm, h ¼ 5mm and f0 ¼ 100Vare shown in Fig. 18.

During their drift, electrons and holes arediffused by multiple collisions, and their distribu-tion follows a Gaussian law:

dN

1ffiffiffiffiffiffiffiffiffiffiffi4pDt

p exp �r2

4Dt

� �dr ð19Þ

where dN=N is the fraction of particles in the lineelement dr at a distance r from the productionpoint and after a time t. The diffusion coefficient D

is related to the mobility by the Einstein law:

D ¼kT

em ð20Þ

where T is the absolute temperature, k is theBoltzmann constant and e is the electron charge. Ifthe mobility is assumed to be a constant, the rmsof the distribution is given by

s ¼ffiffiffiffiffiffiffiffi2Dt

pð21Þ

In our Monte Carlo code the diffusion effect istaken into account during the solution of theequation of motion. The step d~r ¼~rðt þ dtÞ �~rðtÞcan be decomposed in a sum of two terms:

d~r ¼ d~rE þ d~rD ð22Þ

where d~rE is the step due to the electric field, that isevaluated with the Runge–Kutta method, and d~rD

is the step due the random diffusion, that is

given by

d~rD ¼

cosj cos y � sinj � cosj sin y

sinj cos y cosj � sinj sin y

sin y 0 cos y

0BB@

1CCA

s1

s2

s3

0BB@

1CCA (23)

where si are three random values chosen from anormal distribution Nð0; sÞ, with s given by Eq.(21), where y and j are the angles of the directionof motion with respect to the electric field. Aftereach step, a check is done to verify whether thecarrier is still inside the detector volume: if thecarrier has abandoned the detector volume, it isstopped, otherwise a new step is performed.

The current induced at the time t by a movingcarrier on the kth electrode can be evaluated withthe Shockley–Ramo theorem, as

ikðtÞ ¼ �q0~v � ~Ekð~rðtÞÞ ð24Þ

Page 15: A new Monte Carlo code for full simulation of silicon ...mazziot/nima533_322.pdf · A new Monte Carlo code for full simulation of silicon strip detectors M. Brigida, C. Favuzzi, P.

ARTICLE IN PRESS

x (µm)

y (µm)

V (

Vol

t)

050100150200250300350400

-400-300

-200-100

0100

200300

400

0.2

0.4

0.6

0.8

1

Fig. 19. Weighting potential of a strip. The cell used for

computing the potential includes four neighbour strips. The

geometrical parameters are d ¼ 400mm, p ¼ 228mm,

w ¼ 60mm, h ¼ 5mm.

M. Brigida et al. / Nuclear Instruments and Methods in Physics Research A 533 (2004) 322–343336

where q0 is the charge of the carrier,~v is its velocityand ~Ek represents the weighting field associated tothe kth electrode. Eq. (24) holds when impedancefrom electrode to ground is negligible [15]. Theinduced currents are evaluated during the motionof the carriers, and the total currents are computedby adding the contributions from all carriers.

The weighting fields and potentials have beenevaluated by solving the same Maxwell equationfor the true electric field and potential, with r ¼ 0and the following boundary conditions:

fk ¼ 1V ð25Þ

fi ¼ 0; iak ð26Þ

i.e. imposing f ¼ 1V on the kth electrode and f ¼

0 on the other electrodes. The weighting field,~Ekð~rÞ describes the geometrical coupling between acharge at the position ~r and the kth electrode(Fig. 19).

4.2. Simulation of the front-end electronics

A typical detector front-end consists of apreamplifier, that provides a gain, followed by a

shaper, which tailors the overall frequency re-sponse for optimum noise performance, and limitsthe pulse length to accommodate the signal rate.

The analytical relation between input and out-put signals, X ðtÞ and Y ðtÞ, can be derived using theLaplace transform, as briefly outlined below. For alinear system

Y ðsÞ ¼ HðsÞ � X ðsÞ ð27Þ

where HðsÞ is the transfer function of the system,that can be expressed as the ratio between twopolynomials:

HðsÞ ¼AðsÞ

BðsÞ¼

Pi ais

iPj bjsj

: ð28Þ

In cases of our interest, both AðsÞ and BðsÞ aresecond degree polynomials, and HðsÞ is given by

HðsÞ ¼a0 þ a1s þ a2s2

b0 þ b1s þ b2s2ð29Þ

and consequently

ðb0 þ b1s þ b2s2ÞY ðsÞ ¼ ða0 þ a1s þ a2s2ÞX ðsÞ:

ð30Þ

Switching to the time domain, using inverseLaplace transform and assuming zero initialconditions, we get

b0Y ðtÞ þ b1dY ðtÞ

dtþ b2

d2Y ðtÞ

dt2

¼ a0X ðtÞ þ a1dX ðtÞ

dtþ a2

d2X ðtÞ

dt2: (31)

If t is allowed to assume the discrete valuestn ¼ nT s, where T s is the sampling time and n ¼

0; 1; 2; . . . ; the first two derivatives of a genericfunction f ðtÞ can be replaced by

f 0n ¼

df ðtnÞ

dt¼

f nþ1 � f n

T sð32Þ

f 00

n ¼d2f ðtnÞ

dt2¼

f nþ1 � 2f n þ f n�1

T2s

: ð33Þ

Hence, Eq. (31) becomes a recursive equation, thatcan be easily solved by starting from the initialvalues X�1 ¼ X 0 ¼ 0 and Y�1 ¼ Y 0 ¼ 0.

Considering now our front-end, each inducedcurrent signal ikðtÞ is fed into the electronic chain,that has been analytically simulated assuming the

Page 16: A new Monte Carlo code for full simulation of silicon ...mazziot/nima533_322.pdf · A new Monte Carlo code for full simulation of silicon strip detectors M. Brigida, C. Favuzzi, P.

ARTICLE IN PRESS

Zd

Zf

Vshap

Iin

C

R

GR

C

+-

Fig. 20. Equivalent electronic readout chain.

Zd

Vpreampl

Zf

ind inb ina

vna

inf

+-

+

-

Fig. 21. Noise equivalent electronic circuit.

M. Brigida et al. / Nuclear Instruments and Methods in Physics Research A 533 (2004) 322–343 337

scheme shown in Fig. 20, where a preamplifier isfollowed by a shaper. With reference to Fig. 20, Zd

is the detector impedance, Zf is the preamplifierfeedback impedance, C and R are the capacitanceand the resistance of the CR � RC shaper and G isthe gain of the shaper. The impedance Zd (Zf )consists of a capacitance Cd (Cf ) in parallel with aresistance Rd (Rf ).

The open loop voltage gain of the preamplifier,assuming a cascade configuration, is given by

A ¼ AðsÞ ¼ �gmZ0 ¼ �gm

R0

1 þ sR0C0ð34Þ

where gm is the input transistor transconductanceand Z0 is the amplifier load, that consists of acapacitance C0 in parallel with a resistance R0. Byusing the Miller theorem, we can set

Z0f ¼

Zf

1 þ A: ð35Þ

With this parameterization we get

VpðsÞ ¼ AðsÞI inðsÞZ0

fZd

Z0f þ Zd

ð36Þ

and consequently

HpðsÞ ¼VpðsÞ

I inðsÞ¼ AðsÞ

Z0fZd

Z0f þ Zd

: ð37Þ

Hence, for the preamplifier we have

HpðsÞ ¼a0

1 þ b1s þ b2s2ð38Þ

where a0, b1 and b2 are functions of all theelectrical parameters of the circuit.

In the same way, the shaper transfer functionhas been evaluated. For a CR � RC shaper

we have

HshðsÞ ¼V shðsÞ

VpðsÞ¼ G

sCR

ð1 þ sCRÞ2

ð39Þ

and consequently

HshðsÞ ¼a1s

1 þ b3s þ b4s2ð40Þ

where a1, b3 and b4 are functions of the electricalparameters RC of the shaper.

4.3. Noise simulation

The electronic noise is due both to the detectorand to the electronic front-end. Fig. 21 shows theequivalent electronic circuit used for the noiseanalysis. The main noise sources are the following:

the leakage current I l of the semiconductordetector, that can be represented by a currentnoise generator ind in parallel with the detector,and is a shot noise;

the noise due to the bias resistor, that can beassumed as a current generator inb in parallelwith the detector, and is a thermal noise;

the noise from the feedback resistance, that canbe described by a current generator inf and is athermal noise;

the electronic noise of the amplifier, that isdescribed by a combination of a voltage sourcevna and a current source ina at its input.

Shot noise and thermal noise have a whitefrequency distribution, i.e. the spectral powerdensities are constant: i2nd ¼ 2eI l, i2nb ¼ 4kT=Rb

and i2nf ¼ 4kT=Rf, where e is the electron charge,I l is the detector leakage current, k is theBoltzmann constant and T is the temperature.

Page 17: A new Monte Carlo code for full simulation of silicon ...mazziot/nima533_322.pdf · A new Monte Carlo code for full simulation of silicon strip detectors M. Brigida, C. Favuzzi, P.

ARTICLE IN PRESS

Time (µs)In

put

curr

ent

(µA

)

0

0.05

0.1

0.15

0.2

0.25

9.8 9.85 9.9 9.95 10 10.05 10.1 10.15 10.2 10.25 10.3

Fig. 22. Example of a current pulse. The noise current pulse is

superimposed to the signal pulse.

M. Brigida et al. / Nuclear Instruments and Methods in Physics Research A 533 (2004) 322–343338

The amplifier voltage noise consists of a whitecontribution, vnaw, and a ‘‘1=f ’’ contribution, e.g.v2naf ¼ v2

naw þ Af=f . Assuming a FET-input cas-cade, the amplifier current noise and ‘‘1=f ’’ noiseare typically negligible, due to the shaper transferfunction, while the typical value of v2

naw is2:7 kT=gm.

The rms noise voltage, vrms, at the shaper outputis vrms ¼

ffiffiffiffiffiv2

n

p. The value of v2

n is given by

v2n ¼

Z 1

0

ðði2nd þ i2nb þ i2nf ÞjG1ðoÞj2

þ v2nawjG2ðoÞj2Þ

do2p

(41)

where G1ðsÞ ¼ HpðsÞHshðsÞ and G2 ¼ H2ðsÞHshðsÞ

are the gain functions of each noise contribution.The function H2ðsÞ is given by

H2ðsÞ ¼ AðsÞZf þ Zd

Zf þ ð1 þ AðsÞÞZdð42Þ

Switching to the time domain, we get

v2n ¼

1

2

Z 1

�1

ðði2nd þ i2nb þ i2nf Þjg1ðtÞj2

þ v2nawjg2ðtÞj

2Þdt (43)

where g1ðtÞ and g2ðtÞ are the Fourier transform ofG1ðoÞ and G2ðoÞ respectively. Eq. (43) has beenobtained from Eq. (41) applying the Parsevalrelation.

In order to simulate the noise in the discrete-time domain, current and voltage pulses aregenerated according to the formula:

X ðtÞ ¼X

k

Akdðt � kT sÞ ð44Þ

where T s is the sampling time and the amplitudes

Ak are taken from a gaussian distribution with s ¼

en=ffiffiffiffiffiffiffiffi2T s

p(en corresponds to

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffii2nd þ i2nb þ i2nf

qandffiffiffiffiffiffiffiffiffi

v2naw

p). The noise pulses are processed by the

electronic chain, taking into account the correcttransfer function at preamplifier level. The result-ing voltage noise pulse is added to the signal pulseat the output of the preamplifier stage and thissignal is then transferred to shaper stage.

In our simulation we assumed the followingvalues for the parameters of the electronic front-end: Cd ¼ 10 pF (Cd takes into account the bulk

and the interstrip capacitances), Rd ¼ 200GO,Cf ¼ 0:6 pF, Rf ¼ 100MO, R0 ¼ 300 kO,C0 ¼ 1 pF, gm ¼ 0:7mS, I l ¼ 1 nA andRb ¼ 40MO, T s ¼ 1 ns. The shaper is a simpleCR � RC with a peaking time of 1:5 ms (C ¼ 10pFand R ¼ 150 kO) and gain G ¼ 20. The overall gainof this front-end circuit is of about 11:3mV=fC.

For example, let us assume that a current signalwith the shape shown in Fig. 22 has beengenerated by the detector. Figs. 22–24 show,respectively, the input current pulse, the pream-plifier voltage output pulse and the shaper voltageoutput pulse as a function of a time. The noisepulses are superimposed.

Fig. 25 shows the distribution of the overallnoise voltage amplitude at the shaper output. Therms of this distribution is about 0:53mV, corre-sponding to an ENC of about 300 e� rms: thisvalue is equal to the one evaluated with Eq. (41).

5. Monte Carlo results

In this section we present some results of oursimulation, concerning the charge sharing effect.Two different detector geometries have beeninvestigated. In the first case (large pitch

Page 18: A new Monte Carlo code for full simulation of silicon ...mazziot/nima533_322.pdf · A new Monte Carlo code for full simulation of silicon strip detectors M. Brigida, C. Favuzzi, P.

ARTICLE IN PRESS

Time (µs)

Shap

er v

olta

ge o

utpu

t (m

V)

-70

-60

-50

-40

-30

-20

-10

0

0 10 20 30 40 50

Fig. 24. Voltage pulse from the shaper. The noise voltage is

superimposed to the signal.

MeanRMS

0.8545E-02 0.5274

Vrms (mV)N

umbe

r of

ent

ries

0

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

-5 -4 -3 -2 -1 0 1 2 3 4 5

Fig. 25. Noise voltage amplitude distribution.

Time (µs)

Pre

ampl

ifie

r vo

ltag

e ou

tput

(m

V)

-10

-8

-6

-4

-2

0

10 12.5 15 17.5 20 22.5 25 27.5 30

Fig. 23. Voltage pulse from the preamplifier. The noise voltage

pulse is superimposed to the signal.

M. Brigida et al. / Nuclear Instruments and Methods in Physics Research A 533 (2004) 322–343 339

configuration) the values of parameters ared ¼ 400mm, p ¼ 228 mm, w ¼ 60mm, h ¼ 5mmand f0 ¼ 100V; in the second case (small pitch

configuration) the values are d ¼ 325mm,

p ¼ 25mm, w ¼ 12mm, h ¼ 5mm and f0 ¼ 100V.Small pitch configurations are generally chosen forvertex detectors to be mounted in collider experi-ments, where an extreme tracking precision isrequired. On the other hand, large pitch config-urations are typically adopted in silicon telescopesto be operated on satellite borne experiments,where a compromise between a good trackingprecision and a limited number of electronicchannels has to be reached.

The charge sharing occurs whenever a particlereleases its energy in a region including two ormore strips. The greater is the track length insidethe detector, the larger are the effects of chargesharing. These effects become particularly evidentfor tracks crossing the detector at large zenithangles, but they can be seen also for vertical tracksreleasing delta rays.

To understand how the charge sharing isdescribed by our code, we have simulated for boththe detector configurations described above, asample of 3GeV=c pions crossing the silicon layersat null zenith angle. The input track positions areuniformly distributed in the interval from y0 � p=2to y0 þ p=2, where y0 is the position of the centralstrip. In this case the charge sharing is not due tothe inclination of the tracks, but to the weightingfields and to the effects of diffusion.

Page 19: A new Monte Carlo code for full simulation of silicon ...mazziot/nima533_322.pdf · A new Monte Carlo code for full simulation of silicon strip detectors M. Brigida, C. Favuzzi, P.

ARTICLE IN PRESS

M. Brigida et al. / Nuclear Instruments and Methods in Physics Research A 533 (2004) 322–343340

For each event we have considered the chargesqj collected by the strips and we have studied thequantities:

f j ¼qjPk qk

ð45Þ

The charges qj have been evaluated by integratingthe current signals. Fig. 26 shows the averagefractions of charge collected by the various stripsfor both the configurations we have examined.From this figure, it is evident that the chargesharing effect increases when the strip pitchdecreases.

To estimate the contribution of diffusion to thecharge sharing, we have studied, for both theconfigurations, the motion of the charge carriers.If there were no diffusion, charge carriers shouldpropagate inside a silicon cell along the electricfield drift lines shown in Fig. 18. Hence, allthe holes generated along a vertical track withy0 � p=2oyoy0 þ p=2 should be collected by thecentral strip, and no holes should be collected by

Number of strips

Ave

rage

fra

ctio

n of

cha

rge

0

0.2

0.4

0.6

0.8

1

-6 -4 -2 0 2 4 6

Fig. 26. Average fractions of charge collected by the readout

strips for a sample of 3GeV=c pions, entering the detector at

null zenith angle. The input track positions are uniformly

distributed in the range from y0 � p=2 to y0 þ p=2, where y0 ¼

0 is the position of the central strip, and p is the strip pitch. The

results for both the large pitch (solid line) and small pitch

(dashed line) configurations are shown.

the adjacent strips. Due to diffusion, charge carrierscan deviate from the electric field drift lines, andsometimes can migrate from a cell to another. Forthis reason, not all the holes will be collected by thecentral strip, but part of them will be collected bythe adjacent strips, as shown in Fig. 27, where thefractions of holes collected by the various strips areshown for both the configurations. Comparing Fig.27 with Fig. 26, one can see that a part of the chargesharing effect is caused by diffusion. The remainingpart is due to the weighting fields, that areresponsible of the coupling between the movingcharges and the electrodes.

The charge sharing has been also studied bymeans of the Z function. For a particle crossing thedetector in a region between two strips, the Zfunction is defined as

Z ¼V left

V left þ V rightð46Þ

where V left (V right) are the values of the voltagesignals evaluated at the peaking time on the strips

Number of hit strips

Fra

ctio

n of

car

rier

eve

nts

0

0.2

0.4

0.6

0.8

1

-6 -4 -2 0 2 4 6

Fig. 27. Fraction of holes collected by the readout strips for a

sample of 3GeV=c pions crossing the detector at null zenith

angle. The input track positions are uniformly distributed in

the range from y0 � p=2 to y0 þ p=2, where y0 ¼ 0 is the

absolute position of the central strip, and p is the strip pitch.

The results for both the large pitch (solid line) and small pitch

(dashed line) configurations are shown.

Page 20: A new Monte Carlo code for full simulation of silicon ...mazziot/nima533_322.pdf · A new Monte Carlo code for full simulation of silicon strip detectors M. Brigida, C. Favuzzi, P.

ARTICLE IN PRESS

η

Num

ber

of e

vent

s

0

50

100

150

200

250

300

350

400

0 0.2 0.4 0.6 0.8 1

Fig. 28. Distribution of the Z function for a sample of 3GeV=c

pions crossing the detector at null zenith angle. The results for

the large (solid line) and small (dashed line) pitch configura-

tions are shown.

0

0.2

0.4

0.6

0.8

1

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

ηξ

Fig. 29. Dependence of Z from x for a sample of 3GeV=c pions

crossing the detector at null zenith angle. The results for the

large (solid line) and small (dashed line) pitch configurations

are shown.

Fra

ctio

n of

eve

nts

0

0.02

0.04

0.06

0.08

0.1

0.12

0.14

0.16

0.18

0.2

0 0.2 0.4 0.6 0.8 1

η

Fig. 30. Z distribution for a sample of vertical tracks crossing a

400mm thick SSD with 228mm strip pitch. Solid line:

experimental data; dashed line simulated data.

M. Brigida et al. / Nuclear Instruments and Methods in Physics Research A 533 (2004) 322–343 341

to the left (right) of the track. Fig. 28 shows the Zdistributions for the two configurations we haveexamined. Both the distributions are symmetricaround the value Z ¼ 0:5, as can be expected fromDefinition 46, and exhibit two peaks. The peaksare located at Z ’ 0 and 1 for the large pitch

configuration and at Z ’ 0:2 and 0:8 for the small

pitch configuration. These peaks correspond tolarger fractions of charge being collected by theright or by the left strip, respectively. In the large

pitch configuration the peaks are sharper respect tothe small pitch configuration because there is aweak coupling between adjacent strips.

Fig. 29 shows the Z value as function of thedistance of the track from the left strip, expressedin strip pitch units, i.e. x ¼ ðy � yÞ=p. While inthe large pitch configuration Z shows a stepfunction behaviour (Z ’ 1 for xo0:5 and Z ’ 0for x40:5), in the small pitch configurationit is a smoothly decreasing function of x, thatstarts from a value Z ’ 0:8 for x ¼ 0 and reachesthe value Z ’ 0:2 for x ¼ 1. This confirms that inthe small pitch configuration the coupling betweenthe strips is stronger than in the large pitch

configuration.

In order to check the reliability of our code, wehave compared its predictions with the data takenusing a 400mm thick SSD with 228mm strip pitch

Page 21: A new Monte Carlo code for full simulation of silicon ...mazziot/nima533_322.pdf · A new Monte Carlo code for full simulation of silicon strip detectors M. Brigida, C. Favuzzi, P.

ARTICLE IN PRESS

M. Brigida et al. / Nuclear Instruments and Methods in Physics Research A 533 (2004) 322–343342

that has been exposed to a 3GeV=c pion beam atthe CERN PS T9 facility. Fig. 30 shows acomparison of the measured Z distribution for asample of tracks orthogonal to the detector withthe predictions from our simulation. As can beseen from the figure, there is a quite goodagreement between the experimental data and theMonte Carlo prediction. Nevertheless, the peaksof the experimental Z distribution are weaklyshifted respect to those of the Monte Carlodistribution. This means that probably in oursimulation the coupling between adjacent stripshas been underestimated. This effect could beexplained by taking into account tracking errorsassociated to the experimental data.

6. Conclusions

We have developed a new Monte Carlo fullsimulation code that can be used for designingSSDs both for high energy and space physics. Ourcode includes all the physical processes takingplace in a SSD and allows to study its perfor-mance. It can be used for different detectorgeometries and front-end electronics. A completecalculation of the e–h pairs is performed togenerate the charge carries produced by theabsorption both of virtual, e.g. exchange photonsof charged particles with the silicon, and realphotons. The signal simulation includes a noisesimulation, that takes into account the detectorand the electronic noise contributions.

Our simulation code allows to evaluate the SSDresponse both to charged particles and photons. Acharge sharing analysis has been performed fortwo different SSD configurations, simulating asample of minimum ionizing particles. In a largepitch geometry the strips are weakly coupled, whilein a small pitch silicon detector there is a largecapacitive coupling between adjacent strips. Thisanalysis can be extended to study the SSDperformance in terms of efficiency andspace resolution both for charged particles andX-rays.

Silicon detectors are used in many X-rayapplications, such as X-ray astronomy. Our codeallows to evaluate the number npðEÞ of e–h pairs

created by photons of energy E and consequentlythe mean energy W ðEÞ required to produce an e–hpair. These results can be used to evaluate theefficiency of a silicon X-ray detector as a functionof the X-ray energy, after taking into account theX-ray attenuation coefficient mðEÞ, that can bederived from the complex dielectric function insilicon [4].

The results shown in this paper have beenobtained by simulating silicon detectors operatingat room temperature (T0 ¼ 300K). However, thephysical properties of silicon (i.e. the energy bandgap, the intrinsic carrier density and the diffusionconstants), and consequently the processes ofgeneration and motion of the charge carriers, arestrongly dependent on the temperature. Since ourcode allows to take into account the temperaturedependence of these processes, it can be useful forpredicting the performance of a SSD to beoperated at temperatures far away from T0, as ina space environment.

Acknowledgements

The authors are grateful to Prof. H. Bichsel,who gave us his ionization collision cross-sectionprogram in silicon.

References

[1] H.W. Kraner, et al., IEEE Trans. Nucl. Sci. NS-30 (1)

(1983).

[2] W. Dabrowski, P. Grybos, M. Idzik, Nucl. Instr. and

Meth. A 356 (1995) 241.

[3] G.F. Knoll, Radiation Detection and Measurement,

Wiley, New York, 1999, p. 115.

[4] H. Bichsel, Rev. Mod. Phys. 60 (1988) 663.

[5] W. Blum, L. Rolandi, Particle Detection with Drift

Chambers, Springer-Verlag, Berlin, 1994, p. 7.

[6] Scholze, H. Rabus, G. Ulm, J. Appl. Phys. 84 (1998) 2926.

[7] G.W. Fraser, et al., Nucl. Instr. and Meth. A 350 (1994)

368.

[8] R.C. Alig, S. Bloom, C.W. Struck, Phys. Rev. B 22 (1980)

5565.

[9] R.C. Alig, Phys. Rev. B 27 (1983) 968.

[10] A. Owens, et al., Nucl. Instr. and Meth. A 382 (1996) 503.

[11] C.S. Wang, B.M. Klein, Phys. Rev. B 24 (1981) 3393.

[12] W. Shockley, J. Appl. Phys. 9 (1938) 635.

Page 22: A new Monte Carlo code for full simulation of silicon ...mazziot/nima533_322.pdf · A new Monte Carlo code for full simulation of silicon strip detectors M. Brigida, C. Favuzzi, P.

ARTICLE IN PRESS

M. Brigida et al. / Nuclear Instruments and Methods in Physics Research A 533 (2004) 322–343 343

[13] S. Ramo, Currents induced by electron motion, in:

Proceeding of the I.R.E., September 1939, p. 584.

[14] G. Cavalleri, et al., Nucl. Instr. and Meth. 21 (1963) 177.

[15] E. Gatti, G. Padovini, V. Radeka, Nucl. Instr. and Meth.

193 (1982) 651.

[16] E. Gatti, P.F. Manfredi, Rivista del Nuovo Cimento 9

(1986) 1.

[17] C. Jacoboni, et al., Solid State Electron. 20 (1977) 77.

[18] Maxwell 2D, field simulator, v7.0, http://wwwce.web.

cern.ch/wwwce/ae/Maxwell/Maxwell.html.