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