/[PAMELA software]/gpamela/gpspe/gpudiffusion.F
ViewVC logotype

Annotation of /gpamela/gpspe/gpudiffusion.F

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3.1 - (hide annotations) (download)
Tue May 2 11:56:23 2006 UTC (18 years, 7 months ago) by bottai
Branch: MAIN
CVS Tags: v4r6, v4r7, v4r8, v4r9, v4r10, v4r11
subroutine for charge diffusion into the silicon sensors

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

  ViewVC Help
Powered by ViewVC 1.1.23