/[PAMELA software]/PamVMC_update/trk/src/f77/gpudiffusion.F
ViewVC logotype

Annotation of /PamVMC_update/trk/src/f77/gpudiffusion.F

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (hide annotations) (download)
Tue Oct 15 15:51:25 2013 UTC (11 years, 1 month ago) by formato
Branch point for: MAIN, rel
Initial revision

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

  ViewVC Help
Powered by ViewVC 1.1.23