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