| 1 |
nikolas |
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 |
|
|
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