1 |
bottai |
3.1 |
SUBROUTINE GPUDIFFUSION(IACT,TRAPAR,NUMVOL,DELOSS,STEP,ITYPAR) |
2 |
|
|
******************************************************************************** |
3 |
|
|
* |
4 |
|
|
* To perform diffusion of electron and holes bunch inside the silicon |
5 |
|
|
* detectors of the spectrometer |
6 |
|
|
* |
7 |
|
|
* Variables definition: |
8 |
|
|
* IN: |
9 |
|
|
* IACT, integer specifing the action to be taken. It is the INWVOL |
10 |
|
|
* variable in GCTRAK common |
11 |
|
|
* 0: Track is inside a volume |
12 |
|
|
* 1: Entering a new volume or is a new track |
13 |
|
|
* 2: Track is exiting current volume |
14 |
|
|
* TRAPAR, track parameter, is the VECT vector in GCTRAK common (x,y,z..) |
15 |
|
|
* NUMVOL, integr array of numbers identifying the DETECTOR (NUMBV di gustep) |
16 |
|
|
* |
17 |
|
|
* DELOSS, energy loss in the step (GeV) |
18 |
|
|
* ITYPAR, id particella della traccia(vhit(9) che sarà iparspe nell'entupla finale) |
19 |
|
|
* OUT: |
20 |
|
|
* |
21 |
|
|
* |
22 |
|
|
* Called by: GPUSPE |
23 |
|
|
* Author: Elena Taddei, 04/08/2005 , S. Bottai 30/01/06 |
24 |
|
|
* |
25 |
|
|
************************************************************************************* |
26 |
|
|
#include "gpstripspe.inc" |
27 |
|
|
#include "gpgeo.inc" |
28 |
|
|
#include "gpgene.inc" |
29 |
|
|
|
30 |
|
|
INTEGER IACT,NUMVOL(20) |
31 |
|
|
REAL DELOSS, TRAPAR(7),xyzspa(3),VPOS(3),xyzspac(3) |
32 |
|
|
REAL BMAGNET(3),STRPOSL(3),STRPOSG(3) |
33 |
|
|
INTEGER ONCE |
34 |
|
|
DATA ONCE /0/ |
35 |
|
|
SAVE ONCE |
36 |
|
|
|
37 |
|
|
IF(NUMVOL(1).NE.0) THEN |
38 |
|
|
NSPEPLANE=NUMVOL(1) |
39 |
|
|
ELSE IF(NUMVOL(1).EQ.0) THEN |
40 |
|
|
NSPEPLANE=6 |
41 |
|
|
ENDIF |
42 |
|
|
|
43 |
|
|
|
44 |
|
|
VPOS(1)=TRAPAR(1)-STEP/2.*TRAPAR(4) |
45 |
|
|
VPOS(2)=TRAPAR(2)-STEP/2.*TRAPAR(5) |
46 |
|
|
VPOS(3)=TRAPAR(3)-STEP/2.*TRAPAR(6) |
47 |
|
|
|
48 |
|
|
|
49 |
|
|
delossmev=deloss*1000. |
50 |
|
|
|
51 |
|
|
|
52 |
|
|
call gmtod(VPOS,xyzspa,1) |
53 |
|
|
zup=TSPA(3)-xyzspa(3) |
54 |
|
|
zdown=TSPA(3)+xyzspa(3) |
55 |
|
|
|
56 |
|
|
nearstripx=nearstx(xyzspa(1),xyzspa(2)) |
57 |
|
|
if(nearstripx.ne.0) then |
58 |
|
|
|
59 |
|
|
|
60 |
|
|
dx=xyzspa(1)-xstrip(nearstripx) |
61 |
|
|
|
62 |
|
|
**************************************************************************** |
63 |
|
|
* |
64 |
|
|
* X-view strips collect holes, Y-view strips collect electrons. |
65 |
|
|
* Both charge carriers are shifted due to the magnetic field. |
66 |
|
|
* The shift for holes is significant, because it is |
67 |
|
|
* orthogonal to read-out strips. |
68 |
|
|
* A correction for this effect is introduced. |
69 |
|
|
* v is along -Z; B is along -Y --> shift is along -X |
70 |
|
|
* |
71 |
|
|
***************************************************************************** |
72 |
|
|
|
73 |
|
|
CALL GUFLD(VPOS,BMAGNET) |
74 |
|
|
|
75 |
|
|
c |
76 |
|
|
c to be checked |
77 |
|
|
c |
78 |
|
|
xshift=xyzspa(1)-zdown*hallmob*1.e-4*BMAGNET(2)/10. |
79 |
|
|
IF(NSPEPLANE.EQ.6) xshift=xyzspa(1)+ |
80 |
|
|
+ zdown*hallmob*1.e-4*BMAGNET(2)/10. |
81 |
|
|
|
82 |
|
|
|
83 |
|
|
|
84 |
|
|
* |
85 |
|
|
* Now widths of Gaussian functions can be calculated by means of |
86 |
|
|
* the routine sigmadiffus, that gives sigma in m --> *100 --> cm |
87 |
|
|
* |
88 |
|
|
sigxi=amax1(0.00014,100.*sigmadiffus(zdown)) !perchè min=1.4 um? |
89 |
|
|
|
90 |
|
|
* |
91 |
|
|
* Sharing of the charge on strips. |
92 |
|
|
* erf(x) from cernlib computes the (signed) integral of the gaussian |
93 |
|
|
* function from -x to x (sigma=sqrt(1./2.), x0=0). If you have gaussian |
94 |
|
|
* function with x0=a, sigma=b, area between -x and x is obtainable by the |
95 |
|
|
* following formula: |
96 |
|
|
* |
97 |
|
|
* A = erf((x-a)/(sqrt(2.)*b)) A>0 if x-a>0; A<0 if x-a<0 |
98 |
|
|
* |
99 |
|
|
* erfc(x) (ALWAYS > 0) computes the complementary function, i.e. |
100 |
|
|
* 2*integral between x and +infinity |
101 |
|
|
* --> 0.5*erfc(x)=area of the gaussian between x to +inf. |
102 |
|
|
* |
103 |
|
|
|
104 |
|
|
|
105 |
|
|
NSTRIPLOW=MIN(23,NEARSTRIPX) |
106 |
|
|
NSTRIPHIGH=MIN(15,NSTRIPX-NEARSTRIPX) |
107 |
|
|
|
108 |
|
|
DO J=(NEARSTRIPX-NSTRIPLOW+8),(NEARSTRIPX+NSTRIPHIGH-6) |
109 |
|
|
xqdivjm1=xstrip(j)-pitchx/2. |
110 |
|
|
xqdivj=xstrip(j)+pitchx/2. |
111 |
|
|
qfract=0.5*erfc((xqdivjm1-xshift)/(sqrt(2.)*sigxi)) |
112 |
|
|
+ -0.5*erfc((xqdivj-xshift)/(sqrt(2.)*sigxi)) |
113 |
|
|
|
114 |
|
|
proxtanti(NSPEPLANE,numvol(2),j)= |
115 |
|
|
+ proxtanti(NSPEPLANE,numvol(2),j)+delossmev*qfract |
116 |
|
|
IF(GLOBSTRIPX(NSPEPLANE,NUMVOL(2),J).EQ.0.) THEN |
117 |
|
|
STRPOSL(1)=XSTRIP(J) |
118 |
|
|
STRPOSL(2)=0. |
119 |
|
|
STRPOSL(3)=0. |
120 |
|
|
CALL GDTOM(STRPOSL,STRPOSG,1) |
121 |
|
|
GLOBSTRIPX(NSPEPLANE,NUMVOL(2),J)=STRPOSG(1) |
122 |
|
|
ENDIF |
123 |
|
|
|
124 |
|
|
enddo |
125 |
|
|
endif |
126 |
|
|
|
127 |
|
|
|
128 |
|
|
|
129 |
|
|
|
130 |
|
|
|
131 |
|
|
nearstripy=nearsty(xyzspa(1),xyzspa(2)) |
132 |
|
|
if(nearstripy.ne.0) then |
133 |
|
|
|
134 |
|
|
dy=xyzspa(2)-ystrip(nearstripy) |
135 |
|
|
|
136 |
|
|
|
137 |
|
|
sigyi=amax1(0.00023,100.*sigmadiffus(zup ) ) !perchè min=2.3 um? |
138 |
|
|
|
139 |
|
|
* |
140 |
|
|
* The standard deviation on the Y side is increased |
141 |
|
|
* according to a parabolic behaviour + a constant term near p-stop |
142 |
|
|
* |
143 |
|
|
py=pitchy |
144 |
|
|
if (abs(dy).lt.abs((py-psy2)/2.)) then |
145 |
|
|
sigyi=sigyi-psy1*(dy**2)+(py-psy2)*psy1*abs(dy) |
146 |
|
|
else |
147 |
|
|
sigyi=sigyi-psy1*(((py-psy2)/2.)**2) |
148 |
|
|
+ +(py-psy2)*psy1*abs((py-psy2)/2.) |
149 |
|
|
endif |
150 |
|
|
|
151 |
|
|
|
152 |
|
|
NSTRIPLOW=MIN(7,NEARSTRIPY) |
153 |
|
|
NSTRIPHIGH=MIN(7,NSTRIPY-NEARSTRIPY) |
154 |
|
|
do j=(NEARSTRIPY-NSTRIPLOW+1),(NEARSTRIPY+NSTRIPHIGH) |
155 |
|
|
yqdivjm1=ystrip(j)-py/2. |
156 |
|
|
yqdivj=ystrip(j)+py/2. |
157 |
|
|
qfract=0.5*erfc((yqdivjm1-xyzspa(2))/(sqrt(2.)*sigyi)) |
158 |
|
|
+ -0.5*erfc((yqdivj-xyzspa(2))/(sqrt(2.)*sigyi)) |
159 |
|
|
|
160 |
|
|
proytanti(NSPEPLANE,numvol(2),j)= |
161 |
|
|
+ proytanti(NSPEPLANE,numvol(2),j)+delossmev*qfract |
162 |
|
|
|
163 |
|
|
IF(GLOBSTRIPY(NSPEPLANE,NUMVOL(2),J).EQ.0.) THEN |
164 |
|
|
STRPOSL(1)=0. |
165 |
|
|
STRPOSL(2)=YSTRIP(J) |
166 |
|
|
STRPOSL(3)=0. |
167 |
|
|
CALL GDTOM(STRPOSL,STRPOSG,1) |
168 |
|
|
GLOBSTRIPY(NSPEPLANE,NUMVOL(2),J)=STRPOSG(2) |
169 |
|
|
ENDIF |
170 |
|
|
|
171 |
|
|
enddo |
172 |
|
|
|
173 |
|
|
endif |
174 |
|
|
|
175 |
|
|
|
176 |
|
|
|
177 |
|
|
|
178 |
|
|
|
179 |
|
|
END |
180 |
|
|
|
181 |
|
|
* |
182 |
|
|
* //////////////////////////////////////////////////////////////////////////////////////// |
183 |
|
|
* |
184 |
|
|
real function sigmadiffus(zp) |
185 |
|
|
********************************************************************* |
186 |
|
|
* Width of the Gaussian function due to diffusion spread is found. |
187 |
|
|
* x,y,z : where charge is generated (position in given in cm) |
188 |
|
|
* As output standard deviation (m) due to diffusion in silicon |
189 |
|
|
* Diffusion coefficients are proportional to mobility: D=kTm/q, |
190 |
|
|
* where m is mobility: this is true in the Internatinal System |
191 |
|
|
* of units, not in GCS. We compute this quantity in the |
192 |
|
|
* I.S. (renormalitation for m --> cm has been taken into account: |
193 |
|
|
* zpsi=zp/100. ! cm --> m |
194 |
|
|
* Efield=Efield*100. ! V/cm --> V/m --> 10^-4 ) |
195 |
|
|
* WARNING!! Sigma is independent on the carrier mobility m, |
196 |
|
|
* because hdiff = c*m but time = c/m. As a consequence, |
197 |
|
|
* sigma is independent on the dopant concentration. |
198 |
|
|
* E-h pairs created are mostly confined in a tube of about 1 um diameter. |
199 |
|
|
********************************************************************** |
200 |
|
|
#include "gpstripspe.inc" |
201 |
|
|
|
202 |
|
|
zm=zp/100. ! cm --> m |
203 |
|
|
Evm=ebias*100. ! V/cm --> V/m |
204 |
|
|
|
205 |
|
|
vdepl=55. |
206 |
|
|
vappl=70. |
207 |
|
|
thick=3.e-4 |
208 |
|
|
* |
209 |
|
|
* timemu = collection time * mobility |
210 |
|
|
* |
211 |
|
|
timemu=abs(-(thick**2/(2.*vdepl))*log(1-(2*vdepl*zm)/ |
212 |
|
|
+ ((vdepl+vappl)*thick))) |
213 |
|
|
|
214 |
|
|
sigmadiffus=sqrt((2.*boltis*temperature*timemu)/eis)+dsigma |
215 |
|
|
|
216 |
|
|
return |
217 |
|
|
end |
218 |
|
|
|
219 |
|
|
* //////////////////////////////////////////////////////////// |
220 |
|
|
|
221 |
|
|
|
222 |
|
|
|
223 |
|
|
real function xstrip(j) |
224 |
|
|
cv parameter......... |
225 |
|
|
#include "gpstripspe.inc" |
226 |
|
|
parameter (jlastx=2042) |
227 |
|
|
parameter (xlast=5.333/2.-0.07315) |
228 |
|
|
parameter (jfirstx=8) |
229 |
|
|
parameter (xfirst=0.07315-5.333/2.) |
230 |
|
|
|
231 |
|
|
px=pitchx |
232 |
|
|
py=pitchy |
233 |
|
|
if(j.lt.jfirstx.or.j.gt.jlastx) then |
234 |
|
|
write(6,*) 'error , stripx=',j,'not existing' |
235 |
|
|
xstrip=-1.e10 |
236 |
|
|
endif |
237 |
|
|
xstrip=(j-jfirstx)*px+xfirst |
238 |
|
|
|
239 |
|
|
end |
240 |
|
|
|
241 |
|
|
|
242 |
|
|
real function ystrip(j) |
243 |
|
|
cv parameter......... |
244 |
|
|
#include "gpstripspe.inc" |
245 |
|
|
parameter (jlasty=1024) |
246 |
|
|
parameter (ylast=7./2.-0.09855) |
247 |
|
|
parameter (jfirsty=1) |
248 |
|
|
parameter (yfirst=0.0985-7./2.) |
249 |
|
|
|
250 |
|
|
px=pitchx |
251 |
|
|
py=pitchy |
252 |
|
|
if(j.lt.jfirsty.or.j.gt.jlasty) then |
253 |
|
|
write(6,*) 'error , stripy=',j,'not existing' |
254 |
|
|
ystrip=-1.e10 |
255 |
|
|
endif |
256 |
|
|
ystrip=(j-jfirsty)*py+yfirst |
257 |
|
|
|
258 |
|
|
end |
259 |
|
|
|
260 |
|
|
|
261 |
|
|
|
262 |
|
|
function nearstx(x,y) |
263 |
|
|
cv parameter......... |
264 |
|
|
#include "gpstripspe.inc" |
265 |
|
|
parameter (jlastx=2042) |
266 |
|
|
parameter (xlast=5.333/2.-0.07315) |
267 |
|
|
parameter (jfirstx=8) |
268 |
|
|
parameter (xfirst=0.07315-5.333/2.) |
269 |
|
|
parameter (y1xstrip=0.1117-7./2.) |
270 |
|
|
parameter (y2xstrip=7./2.-0.09) |
271 |
|
|
|
272 |
|
|
px=pitchx |
273 |
|
|
py=pitchy |
274 |
|
|
if(x.lt.(xfirst-px/2.).or.x.gt.(xlast+px/2.)) then |
275 |
|
|
nearstx=0 |
276 |
|
|
return |
277 |
|
|
endif |
278 |
|
|
if(y.lt.y1xstrip.or.y.gt.y2xstrip) then |
279 |
|
|
nearstx=0 |
280 |
|
|
return |
281 |
|
|
endif |
282 |
|
|
|
283 |
|
|
nearstx=int((x-xfirst)/px)+jfirstx |
284 |
|
|
if( (x-xstrip(nearstx)).gt.(px/2.) ) nearstx=nearstx+1 |
285 |
|
|
|
286 |
|
|
|
287 |
|
|
end |
288 |
|
|
|
289 |
|
|
function nearsty(x,y) |
290 |
|
|
cv parameter......... |
291 |
|
|
#include "gpstripspe.inc" |
292 |
|
|
parameter (jlasty=1024) |
293 |
|
|
parameter (ylast=7./2.-0.09855) |
294 |
|
|
parameter (jfirsty=1) |
295 |
|
|
parameter (yfirst=0.0985-7./2.) |
296 |
|
|
|
297 |
|
|
parameter (x1ystrip=0.0894-5.333/2.) |
298 |
|
|
parameter (x2ystrip=5.333/2.-0.1221) |
299 |
|
|
|
300 |
|
|
px=pitchx |
301 |
|
|
py=pitchy |
302 |
|
|
if(y.lt.(yfirst-py/2.).or.y.gt.(ylast+py/2.)) then |
303 |
|
|
nearsty=0 |
304 |
|
|
return |
305 |
|
|
endif |
306 |
|
|
if(x.lt.x1ystrip.or.x.gt.x2ystrip) then |
307 |
|
|
nearsty=0 |
308 |
|
|
return |
309 |
|
|
endif |
310 |
|
|
|
311 |
|
|
nearsty=int((y-yfirst)/py)+jfirsty |
312 |
|
|
if( (y-ystrip(nearsty)).gt.(py/2.) ) nearsty=nearsty+1 |
313 |
|
|
|
314 |
|
|
|
315 |
|
|
|
316 |
|
|
end |