SUBROUTINE GPUDIFFUSION(IACT,TRAPAR,NUMVOL,DELOSS,STEP,ITYPAR) ******************************************************************************** * * To perform diffusion of electron and holes bunch inside the silicon * detectors of the spectrometer * * Variables definition: * IN: * IACT, integer specifing the action to be taken. It is the INWVOL * variable in GCTRAK common * 0: Track is inside a volume * 1: Entering a new volume or is a new track * 2: Track is exiting current volume * TRAPAR, track parameter, is the VECT vector in GCTRAK common (x,y,z..) * NUMVOL, integr array of numbers identifying the DETECTOR (NUMBV di gustep) * * DELOSS, energy loss in the step (GeV) * ITYPAR, id particella della traccia(vhit(9) che sarà iparspe nell'entupla finale) * OUT: * * * Called by: GPUSPE * Author: Elena Taddei, 04/08/2005 , S. Bottai 30/01/06 * ************************************************************************************* #include "gpstripspe.inc" #include "gpgeo.inc" #include "gpgene.inc" INTEGER IACT,NUMVOL(20) REAL DELOSS, TRAPAR(7),xyzspa(3),VPOS(3),xyzspac(3) REAL BMAGNET(3),STRPOSL(3),STRPOSG(3) INTEGER ONCE DATA ONCE /0/ SAVE ONCE IF(NUMVOL(1).NE.0) THEN NSPEPLANE=NUMVOL(1) ELSE IF(NUMVOL(1).EQ.0) THEN NSPEPLANE=6 ENDIF VPOS(1)=TRAPAR(1)-STEP/2.*TRAPAR(4) VPOS(2)=TRAPAR(2)-STEP/2.*TRAPAR(5) VPOS(3)=TRAPAR(3)-STEP/2.*TRAPAR(6) delossmev=deloss*1000. call gmtod(VPOS,xyzspa,1) zup=TSPA(3)-xyzspa(3) zdown=TSPA(3)+xyzspa(3) nearstripx=nearstx(xyzspa(1),xyzspa(2)) if(nearstripx.ne.0) then dx=xyzspa(1)-xstrip(nearstripx) **************************************************************************** * * X-view strips collect holes, Y-view strips collect electrons. * Both charge carriers are shifted due to the magnetic field. * The shift for holes is significant, because it is * orthogonal to read-out strips. * A correction for this effect is introduced. * v is along -Z; B is along -Y --> shift is along -X * ***************************************************************************** CALL GUFLD(VPOS,BMAGNET) c c to be checked c xshift=xyzspa(1)-zdown*hallmob*1.e-4*BMAGNET(2)/10. IF(NSPEPLANE.EQ.6) xshift=xyzspa(1)+ + zdown*hallmob*1.e-4*BMAGNET(2)/10. * * Now widths of Gaussian functions can be calculated by means of * the routine sigmadiffus, that gives sigma in m --> *100 --> cm * sigxi=amax1(0.00014,100.*sigmadiffus(zdown)) !perchè min=1.4 um? * * Sharing of the charge on strips. * erf(x) from cernlib computes the (signed) integral of the gaussian * function from -x to x (sigma=sqrt(1./2.), x0=0). If you have gaussian * function with x0=a, sigma=b, area between -x and x is obtainable by the * following formula: * * A = erf((x-a)/(sqrt(2.)*b)) A>0 if x-a>0; A<0 if x-a<0 * * erfc(x) (ALWAYS > 0) computes the complementary function, i.e. * 2*integral between x and +infinity * --> 0.5*erfc(x)=area of the gaussian between x to +inf. * NSTRIPLOW=MIN(23,NEARSTRIPX) NSTRIPHIGH=MIN(15,NSTRIPX-NEARSTRIPX) DO J=(NEARSTRIPX-NSTRIPLOW+8),(NEARSTRIPX+NSTRIPHIGH-6) xqdivjm1=xstrip(j)-pitchx/2. xqdivj=xstrip(j)+pitchx/2. qfract=0.5*erfc((xqdivjm1-xshift)/(sqrt(2.)*sigxi)) + -0.5*erfc((xqdivj-xshift)/(sqrt(2.)*sigxi)) proxtanti(NSPEPLANE,numvol(2),j)= + proxtanti(NSPEPLANE,numvol(2),j)+delossmev*qfract IF(GLOBSTRIPX(NSPEPLANE,NUMVOL(2),J).EQ.0.) THEN STRPOSL(1)=XSTRIP(J) STRPOSL(2)=0. STRPOSL(3)=0. CALL GDTOM(STRPOSL,STRPOSG,1) GLOBSTRIPX(NSPEPLANE,NUMVOL(2),J)=STRPOSG(1) ENDIF enddo endif nearstripy=nearsty(xyzspa(1),xyzspa(2)) if(nearstripy.ne.0) then dy=xyzspa(2)-ystrip(nearstripy) sigyi=amax1(0.00023,100.*sigmadiffus(zup ) ) !perchè min=2.3 um? * * The standard deviation on the Y side is increased * according to a parabolic behaviour + a constant term near p-stop * py=pitchy if (abs(dy).lt.abs((py-psy2)/2.)) then sigyi=sigyi-psy1*(dy**2)+(py-psy2)*psy1*abs(dy) else sigyi=sigyi-psy1*(((py-psy2)/2.)**2) + +(py-psy2)*psy1*abs((py-psy2)/2.) endif NSTRIPLOW=MIN(7,NEARSTRIPY) NSTRIPHIGH=MIN(7,NSTRIPY-NEARSTRIPY) do j=(NEARSTRIPY-NSTRIPLOW+1),(NEARSTRIPY+NSTRIPHIGH) yqdivjm1=ystrip(j)-py/2. yqdivj=ystrip(j)+py/2. qfract=0.5*erfc((yqdivjm1-xyzspa(2))/(sqrt(2.)*sigyi)) + -0.5*erfc((yqdivj-xyzspa(2))/(sqrt(2.)*sigyi)) proytanti(NSPEPLANE,numvol(2),j)= + proytanti(NSPEPLANE,numvol(2),j)+delossmev*qfract IF(GLOBSTRIPY(NSPEPLANE,NUMVOL(2),J).EQ.0.) THEN STRPOSL(1)=0. STRPOSL(2)=YSTRIP(J) STRPOSL(3)=0. CALL GDTOM(STRPOSL,STRPOSG,1) GLOBSTRIPY(NSPEPLANE,NUMVOL(2),J)=STRPOSG(2) ENDIF enddo endif END * * //////////////////////////////////////////////////////////////////////////////////////// * real function sigmadiffus(zp) ********************************************************************* * Width of the Gaussian function due to diffusion spread is found. * x,y,z : where charge is generated (position in given in cm) * As output standard deviation (m) due to diffusion in silicon * Diffusion coefficients are proportional to mobility: D=kTm/q, * where m is mobility: this is true in the Internatinal System * of units, not in GCS. We compute this quantity in the * I.S. (renormalitation for m --> cm has been taken into account: * zpsi=zp/100. ! cm --> m * Efield=Efield*100. ! V/cm --> V/m --> 10^-4 ) * WARNING!! Sigma is independent on the carrier mobility m, * because hdiff = c*m but time = c/m. As a consequence, * sigma is independent on the dopant concentration. * E-h pairs created are mostly confined in a tube of about 1 um diameter. ********************************************************************** #include "gpstripspe.inc" zm=zp/100. ! cm --> m Evm=ebias*100. ! V/cm --> V/m vdepl=55. vappl=70. thick=3.e-4 * * timemu = collection time * mobility * timemu=abs(-(thick**2/(2.*vdepl))*log(1-(2*vdepl*zm)/ + ((vdepl+vappl)*thick))) sigmadiffus=sqrt((2.*boltis*temperature*timemu)/eis)+dsigma return end * //////////////////////////////////////////////////////////// real function xstrip(j) cv parameter......... #include "gpstripspe.inc" parameter (jlastx=2042) parameter (xlast=5.333/2.-0.07315) parameter (jfirstx=8) parameter (xfirst=0.07315-5.333/2.) px=pitchx py=pitchy if(j.lt.jfirstx.or.j.gt.jlastx) then write(6,*) 'error , stripx=',j,'not existing' xstrip=-1.e10 endif xstrip=(j-jfirstx)*px+xfirst end real function ystrip(j) cv parameter......... #include "gpstripspe.inc" parameter (jlasty=1024) parameter (ylast=7./2.-0.09855) parameter (jfirsty=1) parameter (yfirst=0.0985-7./2.) px=pitchx py=pitchy if(j.lt.jfirsty.or.j.gt.jlasty) then write(6,*) 'error , stripy=',j,'not existing' ystrip=-1.e10 endif ystrip=(j-jfirsty)*py+yfirst end function nearstx(x,y) cv parameter......... #include "gpstripspe.inc" parameter (jlastx=2042) parameter (xlast=5.333/2.-0.07315) parameter (jfirstx=8) parameter (xfirst=0.07315-5.333/2.) parameter (y1xstrip=0.1117-7./2.) parameter (y2xstrip=7./2.-0.09) px=pitchx py=pitchy if(x.lt.(xfirst-px/2.).or.x.gt.(xlast+px/2.)) then nearstx=0 return endif if(y.lt.y1xstrip.or.y.gt.y2xstrip) then nearstx=0 return endif nearstx=int((x-xfirst)/px)+jfirstx if( (x-xstrip(nearstx)).gt.(px/2.) ) nearstx=nearstx+1 end function nearsty(x,y) cv parameter......... #include "gpstripspe.inc" parameter (jlasty=1024) parameter (ylast=7./2.-0.09855) parameter (jfirsty=1) parameter (yfirst=0.0985-7./2.) parameter (x1ystrip=0.0894-5.333/2.) parameter (x2ystrip=5.333/2.-0.1221) px=pitchx py=pitchy if(y.lt.(yfirst-py/2.).or.y.gt.(ylast+py/2.)) then nearsty=0 return endif if(x.lt.x1ystrip.or.x.gt.x2ystrip) then nearsty=0 return endif nearsty=int((y-yfirst)/py)+jfirsty if( (y-ystrip(nearsty)).gt.(py/2.) ) nearsty=nearsty+1 end