      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 "gpdgeo.inc" 
#include "gpgene.inc"
#include "gpkey.inc"
#include "gpdkey.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

c        IF(NUMVOL(1).NE.0) THEN
	   NSPEPLANE=NUMVOL(1)
c	ELSE IF(NUMVOL(1).EQ.0) THEN
c           NSPEPLANE=6
c	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
*
*****************************************************************************
         

          IF(FFIELD.NE.0) THEN
          
	  CALL GUFIELD(VPOS,BMAGNET)
 
c          PRINT*,'VPOS:',VPOS(1),' ',VPOS(2),' ',VPOS(3),' ',BMAGNET(2)        
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.
	
          ELSE
             xshift=xyzspa(1)
          ENDIF
c          PRINT*,'NSPEPLENE: ',NSPEPLANE
c          PRINT*,'xyzspa: ', xyzspa(1),' ',xyzspa(2),' ',xyzspa(3)
c          PRINT*,'zdown:', zdown
c          PRINT*,'Bmagnet(2): ',BMAGNET(2)
c          PRINT*,'hallmob: ',hallmob
c          PRINT*,'shift(um):',(xshift-xyzspa(1))*10000,'pl:',NSPEPLANE
*
*        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))
cc        PRINT*,xyzspa(1),' ',xyzspa(2),' ', nearstripx, ' ', nearstripy

        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))
cc             PRINT*,'J ',j, 'qfract', qfract
             proytanti(NSPEPLANE,numvol(2),j)=
     +       proytanti(NSPEPLANE,numvol(2),j)+delossmev*qfract
cc             PRINT*,'proytanti',  proytanti(NSPEPLANE,numvol(2),j)
             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
