| 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 |
#include "gpkey.inc" |
| 30 |
|
| 31 |
INTEGER IACT,NUMVOL(20),numsens |
| 32 |
REAL DELOSS, TRAPAR(7),xyzspa(3),VPOS(3),xyzspac(3) |
| 33 |
REAL BMAGNET(3),STRPOSL(3),STRPOSG(3) |
| 34 |
INTEGER ONCE |
| 35 |
DATA ONCE /0/ |
| 36 |
SAVE ONCE |
| 37 |
|
| 38 |
NSPEPLANE=INT((NUMVOL(2)-1)/6)+1 |
| 39 |
numsens=NUMVOL(2)-(NSPEPLANE-1)*6 |
| 40 |
|
| 41 |
|
| 42 |
VPOS(1)=TRAPAR(1)-STEP/2.*TRAPAR(4) |
| 43 |
VPOS(2)=TRAPAR(2)-STEP/2.*TRAPAR(5) |
| 44 |
VPOS(3)=TRAPAR(3)-STEP/2.*TRAPAR(6) |
| 45 |
|
| 46 |
|
| 47 |
delossmev=deloss*1000. |
| 48 |
|
| 49 |
|
| 50 |
call gmtod(VPOS,xyzspa,1) |
| 51 |
zup=TSPA(3)-xyzspa(3) |
| 52 |
zdown=TSPA(3)+xyzspa(3) |
| 53 |
|
| 54 |
|
| 55 |
nearstripx=nearstx(xyzspa(1),xyzspa(2)) |
| 56 |
if(nearstripx.ne.0) then |
| 57 |
|
| 58 |
|
| 59 |
dx=xyzspa(1)-xstrip(nearstripx) |
| 60 |
|
| 61 |
**************************************************************************** |
| 62 |
* |
| 63 |
* X-view strips collect holes, Y-view strips collect electrons. |
| 64 |
* Both charge carriers are shifted due to the magnetic field. |
| 65 |
* The shift for holes is significant, because it is |
| 66 |
* orthogonal to read-out strips. |
| 67 |
* A correction for this effect is introduced. |
| 68 |
* v is along -Z; B is along -Y --> shift is along -X |
| 69 |
* |
| 70 |
***************************************************************************** |
| 71 |
IF(FFIELD.NE.0) THEN |
| 72 |
CALL GUFLD(VPOS,BMAGNET) |
| 73 |
|
| 74 |
c |
| 75 |
c |
| 76 |
xshift=xyzspa(1)+zdown*hallmob*1.e-4*BMAGNET(2)/10. |
| 77 |
IF(NSPEPLANE.EQ.6) xshift=xyzspa(1)- |
| 78 |
+ zdown*hallmob*1.e-4*BMAGNET(2)/10. |
| 79 |
|
| 80 |
else |
| 81 |
xshift=xyzspa(1) |
| 82 |
endif |
| 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,numsens,j)= |
| 115 |
+ proxtanti(NSPEPLANE,numsens,j)+delossmev*qfract |
| 116 |
IF(GLOBSTRIPX(NSPEPLANE,numsens,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,numsens,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,numsens,j)= |
| 161 |
+ proytanti(NSPEPLANE,numsens,j)+delossmev*qfract |
| 162 |
|
| 163 |
IF(GLOBSTRIPY(NSPEPLANE,numsens,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,numsens,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 |
Parent Directory
| 
Revision Log