mcdet

Fortran source

0001 c     Simulation par la methode de Monte Carlo de la reponse d'un
0002 c     scintillateur au rayonnement incident en fonction de sa
0003 c     taille et de la position et de l'extension de la zone eclairee.
0004 c
0005       real hnu0,alpha0,ctheta,stheta,cphi,sphi,hnu
0006       real tableau(10000000,8),x,y,z,r1,xmax,ymax,zmax,mhu
0007       real salpha0,calpha0,mhuph,mhuco,omegax,omegay,omegaz
0008       real dze,dzet,hnu1,cosdev,sindev,qui,omegax1,omegay1,omegaz1
0009       real mhuph1,mhuco1,mhu1,dze1,x1,y1,z1,ctheta1,stheta1
0010       real cphi1,sphi1,mec2,prof,nbtotcoll,nbmoycoll,distri(2,100)
0011       real nbtotlt,nbtotzm
0012       real nbmoylt,nbmoyzm
0013       real a,x0,y0,xr,yr,atom,xmhuph,xmhuco
0014 c
0015       integer n,i,nphot,nbdet,w,reset0,reset,nbcoll,incidence
0016       integer nv,nbabs,choixdistri,iv,nbdiff,type
0017 c
0018       common/etiq1/reset
0019 c
0020 c
0021 c        hnu0 = energie des photons incidents
0022 c        stheta= cos(theta) 
0023 c        cphi= cos(azimuth)
0024 c        hnu= energie du photon
0025 c        x,y,z=position 
0026 c        xr=position du centre du diaphragme / axe Ox
0027 c        yr=position du centre du diaphragme / axe Oy
0028 c        r1 = rayon du collimateur interieur (cm)
0029 c        xmax = demi hauteur de la paroi
0030 c        ymax = demi profondeur de la paroi
0031 c        zmax = epaisseur de la paroi (valeur negative !)
0032 c        n = nombre de photon envoyes
0033 c        nphot = numero du photon traite
0034 c        nbdet=nombre de photons detectes
0035 c        w = survie
0036 c        reset0 et reset = pour l'initialisation de ran2
0037 c        nbcoll = nombre de collisions
0038 c        mhu = coefficient d'absorption de la paroi
0039 c        choixdistri = mode sortie complete ou distributions seules
0040 c        incidence=1:normale au detecteur et 0 isotrope
0041 c
0042       data mec2,pi,rad /511.0,3.1415926535,360/ 
0043 c
0044       open(3,FILE='dmcdet.data',STATUS='UNKNOWN')
0045 
0046 c
0047 c     ************* entree donnees et ecriture ****************
0048 c          
0049       read(3,*) incidence
0050       read(3,*) type
0051       read(3,*) r1
0052       read(3,*) xr
0053       read(3,*) yr
0054       read(3,*) xmax
0055       read(3,*) ymax
0056       read(3,*) zmax
0057       read(3,*) hnu0
0058       read(3,*) n
0059       read(3,*) reset0
0060       read(3,*) choixdistri
0061       read(3,*) nv
0062 c
0063       close(3)
0064 c    
0065 c     --------------------CALCUL ATTENUATION LINEIQUE-------------------
0066 c
0067       if (type.eq.1) then
0068          call xBGO(hnu0,mhuph,mhuco)
0069       elseif (type.eq.2) then
0070          call xNaI(hnu0,mhuph,mhuco)
0071       elseif (type.eq.3) then
0072          call xCsI(hnu0,mhuph,mhuco)
0073       elseif (type.eq.4) then
0074          call xgerma(hnu0,mhuph,mhuco)
0075       elseif (type.eq.5) then
0076          call xCdTe(hnu0,mhuph,mhuco)
0077       elseif (type.eq.6) then
0078          call xLYSO(hnu0,mhuph,mhuco)
0079       endif
0080       mhu=mhuph+mhuco
0081 c
0082 c     ----------------------ECRITURE DES PARAMETRES---------------------
0083 c
0084       open(2,FILE='rmcdet.data',STATUS='UNKNOWN')
0085 c
0086       write(2,*) '----------------------------------------------------'
0087       write(2,*) 'Code MCdet : Interaction radiation-detector :'
0088       write(2,*) '----------------------------------------------------'
0089       write(2,*) ' ' 
0090       write(2,*) 'Photon direction (1:perp et 0: isotropic)--------->',
0091      *incidence
0092       write(2,*) '(1):BGO (2):NaI (3):CsI (4):Ge (5):CdTe (6):LYSO-->',
0093      *type
0094       write(2,1) 'Collimateur radius N1 (detector side) (cm)-------->',
0095      *r1
0096       write(2,1) 'Impact location on the detector along Ox (cm)----->',
0097      *xr
0098       write(2,1) 'Impact location on the detector along Oy (cm)----->',
0099      *yr
0100       write(2,1) 'Half-height of the detector (cm)------------------>',
0101      *xmax
0102       write(2,1) 'Half-width of the detector (cm)------------------->',
0103      *ymax
0104       write(2,1) 'Thickness of the detector (cm)-------------------->',
0105      *zmax
0106       write(2,1) 'Energy of the incident photon (keV)--------------->',
0107      *hnu0
0108       write(2,1) 'Lineic attenuation  (cm-1)------------------------>',
0109      *mhu
0110       write(2,*) 'Number of incoming photons ----------------------->',
0111      *n
0112       write(2,*) 'Initialisation random variable generator (>0)----->',
0113      *reset0
0114       write(2,*) 'Number of histogram steps in energy (<100)------->',
0115      *nv
0116       write(2,*) ' ' 
0117       write(2,*) '----------------------------------------------------'
0118       write(2,*) ' ' 
0119 c
0120 c     ------------MISE AU FORMAT DES DONNEES ET INITIALISATION--------
0121 c
0122 c
0123       reset=-1*reset0
0124       nbdet=0
0125       nbdiff=0
0126       nbabs=0
0127       nphot=1
0128 c
0129   30  nbcoll=1
0130       prof=0.0
0131       dzet=0.0
0132       if (incidence.eq.1) then
0133          alpha0=0.0
0134          cphi=1.0
0135          sphi=0.0
0136       elseif (incidence.eq.0) then
0137          alpha0=acos(ran2(reset))
0138          phi=2*pi*ran2(reset)
0139          cphi=cos(phi)
0140          sphi=sin(phi)
0141       endif
0142       salpha0=sin(alpha0*2*pi/rad)
0143       calpha0=cos(alpha0*2*pi/rad)
0144       ctheta=-1.0*calpha0
0145       stheta=-1.0*salpha0
0146       omegax=stheta*cphi
0147       omegay=stheta*sphi
0148       omegaz=ctheta
0149       call posi(r1,a)
0150       qui=tirazi(reset)
0151       x0=a*cos(qui)+xr
0152       y0=a*sin(qui)+yr
0153       x=x0
0154       y=y0
0155       z=0.0
0156       hnu=hnu0
0157       w=1
0158 c
0159 c     -----------------CALCUL DU PARCOURS DANS LA MATIERE--------------
0160 c
0161       if (type.eq.1) then
0162          call xBGO(hnu,mhuph,mhuco)
0163       elseif (type.eq.2) then
0164          call xNaI(hnu,mhuph,mhuco)
0165       elseif (type.eq.3) then
0166          call xCsI(hnu,mhuph,mhuco)
0167       elseif (type.eq.4) then
0168          call xgerma(hnu,mhuph,mhuco)
0169       elseif (type.eq.5) then
0170          call xCdTe(hnu,mhuph,mhuco) 
0171       elseif (type.eq.6) then
0172          call xLYSO(hnu,mhuph,mhuco)      
0173       endif
0174       mhu=mhuph+mhuco
0175       call tirpos(mhu,dze)
0176       x=x+omegax*dze
0177       y=y+omegay*dze
0178       z=z+omegaz*dze
0179       dzet=dze
0180       prof=z
0181       if ((abs(x).gt.xmax).or.(abs(y).gt.ymax)) then
0182          if (nphot.ge.n) then
0183             goto 50
0184          else 
0185             nphot=nphot+1
0186             goto 30
0187          endif
0188       elseif (z.lt.(-1.0*zmax)) then
0189          if (nphot.ge.n) then
0190             goto 50
0191          else 
0192             nphot=nphot+1
0193             goto 30
0194          endif
0195       endif
0196 c
0197   20  continue
0198       atom=ran2(reset)
0199       if (type.eq.1) then
0200          if ((atom.ge.0.0).and.(atom.le.(4.0/19.0))) then
0201             call xbismu(hnu,xmhuph,xmhuco)
0202          endif
0203          if ((atom.gt.(4.0/19.0)).and.(atom.le.(7.0/19.0))) then
0204             call xgerma(hnu,xmhuph,xmhuco)
0205          endif
0206          if ((atom.gt.(7.0/19.0)).and.(atom.le.1.0)) then
0207             call xoxy(hnu,xmhuph,xmhuco)
0208          endif
0209       elseif (type.eq.2) then
0210          if ((atom.ge.0.0).and.(atom.le.0.5)) then
0211             call xNa(hnu,xmhuph,xmhuco)
0212          endif
0213          if ((atom.gt.0.5).and.(atom.le.1.0)) then
0214             call xI(hnu,xmhuph,xmhuco)
0215          endif
0216       elseif (type.eq.3) then
0217          if ((atom.ge.0.0).and.(atom.le.0.5)) then
0218             call xCs(hnu,xmhuph,xmhuco)
0219          endif
0220          if ((atom.gt.0.5).and.(atom.le.1.0)) then
0221             call xI(hnu,xmhuph,xmhuco)
0222          endif
0223       elseif (type.eq.4) then
0224          call xgerma(hnu,xmhuph,xmhuco)
0225       elseif (type.eq.5) then
0226          if ((atom.ge.0.0).and.(atom.le.0.5)) then
0227             call xCd(hnu,xmhuph,xmhuco)
0228          endif
0229          if ((atom.gt.0.5).and.(atom.le.1.0)) then
0230             call xTe(hnu,xmhuph,xmhuco)
0231          endif
0232       elseif (type.eq.6) then
0233          if ((atom.ge.0.0).and.(atom.le.(1.8/8.0))) then
0234             call xlu(hnu,xmhuph,xmhuco)
0235          endif
0236          if ((atom.gt.(1.8/8.0)).and.(atom.le.(2.0/8.0))) then
0237             call xyt(hnu,xmhuph,xmhuco)
0238          endif
0239          if ((atom.gt.(2.0/8.0)).and.(atom.le.(3.0/8.0))) then
0240             call xsi(hnu,xmhuph,xmhuco)
0241          endif
0242          if ((atom.gt.(3.0/8.0)).and.(atom.le.1.0)) then
0243             call xoxy(hnu,xmhuph,xmhuco)
0244          endif
0245       endif
0246       call tirsur(w,xmhuph,xmhuco)
0247       if (w.eq.0)  then
0248         nbabs=nbabs+1
0249         nbdet=nbdet+1
0250         ibdet=nbdet
0251         tableau(ibdet,1)=hnu0
0252         tableau(ibdet,3)=float(nbcoll)
0253         tableau(ibdet,4)=dzet
0254         tableau(ibdet,5)=(-1.0)*prof
0255         tableau(ibdet,6)=x0
0256         tableau(ibdet,7)=hnu0+gasdev(type,hnu0)
0257         tableau(ibdet,8)=y0
0258         if (nphot.ge.n) then
0259            goto 50
0260         else 
0261            nphot=nphot+1
0262            goto 30
0263         endif
0264       else
0265          nbcoll=nbcoll+1
0266          hnu1=hnu*tircom(hnu)
0267          cosdev=1.0-(1.0/hnu1-1.0/hnu)*mec2
0268          sindev=sqrt(1.0-cosdev**2)
0269          qui=tirazi(reset)
0270          omegax1=sindev*cos(qui)*ctheta*cphi-sindev*sin(qui)*sphi
0271      *   +cosdev*stheta*cphi
0272          omegay1=sindev*cos(qui)*ctheta*sphi+sindev*sin(qui)*cphi
0273      *   +cosdev*stheta*sphi
0274          omegaz1=-sindev*cos(qui)*stheta+cosdev*ctheta
0275          if (hnu1.lt.3.99) hnu1=3.99
0276          if (type.eq.1) then
0277             call xBGO(hnu1,mhuph1,mhuco1)
0278          elseif (type.eq.2) then
0279             call xNaI(hnu1,mhuph1,mhuco1)
0280          elseif (type.eq.3) then
0281             call xCsI(hnu1,mhuph1,mhuco1)
0282          elseif (type.eq.4) then
0283             call xgerma(hnu1,mhuph1,mhuco1)
0284          elseif (type.eq.5) then
0285             call xCdTe(hnu1,mhuph1,mhuco1)
0286          elseif (type.eq.6) then
0287            call xLYSO(hnu1,mhuph1,mhuco1)
0288          endif
0289          mhu1=mhuph1+mhuco1
0290          call tirpos(mhu1,dze1)
0291          x1=x+omegax1*dze1
0292          y1=y+omegay1*dze1
0293          z1=z+omegaz1*dze1
0294          ctheta1=omegaz1
0295          stheta1=sqrt(1.0-omegaz1**2)
0296          if (stheta1.eq.0.0) then
0297            cphi1=1.0
0298            sphi1=0.0
0299          else
0300            cphi1=omegax1/stheta1
0301            sphi1=omegay1/stheta1
0302          endif
0303          if ((abs(x1).gt.xmax).or.(abs(y1).gt.ymax)) then
0304             nbdiff=nbdiff+1
0305             nbdet=nbdet+1
0306             ibdet=nbdet
0307             tableau(ibdet,1)=hnu0-hnu1
0308             tableau(ibdet,3)=float(nbcoll)
0309             tableau(ibdet,4)=dzet
0310             tableau(ibdet,5)=(-1.0)*prof
0311             tableau(ibdet,6)=x0
0312             tableau(ibdet,7)=hnu0+gasdev(type,hnu0-hnu1)-hnu1
0313             tableau(ibdet,8)=y0
0314             if (nphot.ge.n) then
0315                goto 50
0316             else
0317                nphot=nphot+1
0318                goto 30
0319            endif
0320          elseif (z1.lt.(-1.0*zmax)) then
0321            nbdiff=nbdiff+1
0322            nbdet=nbdet+1
0323            ibdet=nbdet
0324            tableau(ibdet,1)=hnu0-hnu1
0325            tableau(ibdet,3)=float(nbcoll)
0326            tableau(ibdet,4)=dzet
0327            tableau(ibdet,5)=(-1.0)*prof
0328            tableau(ibdet,6)=x0
0329            tableau(ibdet,7)=hnu0+gasdev(type,hnu0-hnu1)-hnu1
0330            tableau(ibdet,8)=y0
0331            if (nphot.ge.n) then
0332               goto 50
0333            else
0334               nphot=nphot+1
0335               goto 30
0336            endif
0337         elseif (z1.ge.0.0) then
0338            nbdiff=nbdiff+1
0339            nbdet=nbdet+1
0340            ibdet=nbdet
0341            tableau(ibdet,1)=hnu0-hnu1
0342            tableau(ibdet,3)=float(nbcoll)
0343            tableau(ibdet,4)=dzet
0344            tableau(ibdet,5)=(-1.0)*prof
0345            tableau(ibdet,6)=x0
0346            tableau(ibdet,7)=hnu0+gasdev(type,hnu0-hnu1)-hnu1
0347            tableau(ibdet,8)=y0
0348            if (nphot.ge.n) then
0349               goto 50
0350            else
0351               nphot=nphot+1
0352            goto 30
0353            endif
0354         else
0355            dzet=dzet+dze1
0356            if (z1.lt.prof) prof=z1
0357            omegax=omegax1
0358            omegay=omegay1
0359            omegaz=omegaz1
0360            ctheta=ctheta1
0361            stheta=stheta1
0362            cphi=cphi1
0363            sphi=sphi1
0364            hnu=hnu1
0365            x=x1
0366            y=y1
0367            z=z1
0368            mhuph=mhuph1
0369            mhuco=mhuco1
0370         endif
0371         if (nphot.ge.n) then
0372           goto 50
0373         else
0374           goto 20
0375         endif
0376       endif
0377    50 continue
0378       nbtotcoll=0
0379       nbtotlt=0
0380       nbtotzm=0
0381       do 80 i=1,nbdet
0382          nbtotcoll=nbtotcoll+tableau(i,3)
0383          nbtotlt=nbtotlt+tableau(i,4)
0384          nbtotzm=nbtotzm+tableau(i,5)
0385    80 continue
0386 c
0387       nbmoycoll=nbtotcoll/nbdet
0388       nbmoylt=nbtotlt/nbdet
0389       nbmoyzm=nbtotzm/nbdet
0390 c
0391 c     ----------------------AFFICHAGE DES RESULTATS--------------------
0392 c
0393       write(2,110) 'Number of incoming photons   = ',nphot
0394       write(2,110) 'Number of detected photons   = ',nbdet
0395       write(2,110) 'Number of absorbed photons   = ',nbabs
0396       write(2,110) 'Number of scattered photons  = ',nbdiff
0397       write(2,120) 'Mean number of collisions per detected photon  =',
0398      *nbmoycoll
0399       write(2,120) 'Mean free path per detected photon (cm)        =',
0400      *nbmoylt
0401       write(2,120) 'Mean depth reached per detected photon (cm)    =',
0402      *nbmoyzm
0403       write(2,*) '-----------------------------------------------------'
0404 c
0405       if (choixdistri.eq.1) then
0406       write(2,*) 'i,energie(i),nbcollision(i),Ltotal(i),Zmax(i),a(i),E(i
0407      *)'
0408       write(2,*) '-----------------------------------------------------'
0409       do 60 i=1,nbdet
0410    60 write(2,100) i,tableau(i,1),tableau(i,3),tableau(i,4),tableau(i,5)
0411      *,tableau(i,6),tableau(i,7)
0412       write(2,*) '-----------------------------------------------------'
0413       endif
0414 c
0415       if (nbdet.ne.0) then
0416          iv=1
0417       write(2,*) '-----------------------------------------------------'
0418          write(2,*) 'Energy distribution of absorbed photons : k,f(k)'
0419       write(2,*) '-----------------------------------------------------'
0420          call histo(nv,iv,nbdet,tableau,distri)
0421          do 71 i=1,nv
0422             write(2,150) distri(1,i),distri(2,i)
0423    71    continue
0424          iv=7
0425       write(2,*) '-----------------------------------------------------'
0426          write(2,*) 'Scintillation energy distribution : E,f(E)'
0427       write(2,*) '-----------------------------------------------------'
0428          call histo(nv,iv,nbdet,tableau,distri)
0429          do 78 i=1,nv
0430             write(2,150) distri(1,i),distri(2,i)
0431    78    continue
0432          iv=6
0433       write(2,*) '-----------------------------------------------------'
0434          write(2,*) 'Impact distribution : x0,f(x0)'
0435       write(2,*) '-----------------------------------------------------'
0436          call histo(nv,iv,nbdet,tableau,distri)
0437          do 73 i=1,nv
0438             write(2,150) distri(1,i),distri(2,i)
0439    73    continue
0440          iv=8
0441       write(2,*) '-----------------------------------------------------'
0442          write(2,*) 'Impact distribution : y0,f(y0)'
0443       write(2,*) '-----------------------------------------------------'
0444          call histo(nv,iv,nbdet,tableau,distri)
0445          do 79 i=1,nv
0446             write(2,150) distri(1,i),distri(2,i)
0447    79    continue
0448          iv=3
0449       write(2,*) '-----------------------------------------------------'
0450          write(2,*) 'Number of collisions distribution : nc,f(nc)'
0451       write(2,*) '-----------------------------------------------------'
0452          call histo(nv,iv,nbdet,tableau,distri)
0453          do 75 i=1,nv
0454             write(2,150) distri(1,i),distri(2,i)
0455    75    continue
0456          iv=4
0457       write(2,*) '-----------------------------------------------------'
0458          write(2,*) 'Mean path distribution : lt,f(lt)'
0459       write(2,*) '-----------------------------------------------------'
0460          call histo(nv,iv,nbdet,tableau,distri)
0461          do 76 i=1,nv
0462             write(2,150) distri(1,i),distri(2,i)
0463    76    continue
0464          iv=5
0465       write(2,*) '-----------------------------------------------------'
0466          write(2,*) 'Depth distribution max : zm,f(zm)'
0467       write(2,*) '-----------------------------------------------------'
0468          call histo(nv,iv,nbdet,tableau,distri)
0469          do 77 i=1,nv
0470             write(2,150) distri(1,i),distri(2,i)
0471    77    continue
0472          write(2,*) '--------------------------------------------------'
0473          write(2,*) '--------------------------------------------------'
0474       else
0475          write(2,*) 'Ni photon detected'
0476       endif
0477       close(2)
0478 c
0479 c     -----------------------------FORMATS-----------------------------
0480 c
0481     1 format(a51,f10.2)
0482   100 format(i6,1x,f5.1,1x,f6.2,1x,f6.2,1x,f6.2,1x,f6.2,1x,f6.2,1x,f6.
0483      *2,1x,f6.2,1x,f6.2)
0484   110 format(a30,i7)
0485   120 format(a48,f6.2)
0486   121 format(a47,e10.3)
0487   150 format(e10.3,e10.3)
0488   151 format(a43,f6.2)
0489   152 format(e10.3,e10.3,e10.3,e10.3)
0490   753 format(a51,i6)
0491   752 format(a51,f8.3)
0492   751 format(a51,e8.3)
0493   750 format(a59)
0494       end
0495 
0496 
0497 
0498       subroutine histo(nv,iv,nt,tab,distri)
0499       real vmin,vmax,vmin1,vmax1,tab(10000000,8)
0500       real dv,distri(2,100),cn,ntt
0501       integer nv,i,iv,ivm
0502 c
0503 c         Cette routine determine le spectre
0504 c
0505       do 40 i=1,nt
0506          tab(i,2)=0.0
0507   40  continue
0508       do 41 i=1,nv
0509          distri(1,i)=0.0
0510          distri(2,i)=0.0
0511   41  continue
0512 c
0513       if ((iv.eq.1).or.(iv.eq.7)) then
0514          ivm=7
0515       else
0516          ivm=iv
0517       endif
0518       vmin=+1.0e+30
0519       vmax=-1.0e+30
0520       do 42 i=1,nt
0521         if (tab(i,ivm).lt.vmin) vmin=tab(i,ivm)
0522         if (tab(i,ivm).gt.vmax) vmax=tab(i,ivm)
0523   42  continue
0524       if ((iv.eq.1).or.(iv.eq.7)) vmax=vmax+(vmax-vmin)/nv
0525 c
0526       dv=(vmax-vmin)/nv
0527 c      cn=1.0/(dv*nt)
0528       ntt=0.0
0529       do 20 j=0,nv-1
0530          distri(2,j+1)=0.0
0531          vmin1=vmin+j*dv
0532          vmax1=vmin+(j+1)*dv
0533          if (iv.ne.3) then
0534             distri(1,j+1)=(vmin1+vmax1)/2
0535          else
0536             distri(1,j+1)=vmin1
0537          endif
0538          do 10 i=1,nt
0539            if (tab(i,2).eq.0.0) then
0540              if ((tab(i,iv).ge.vmin1).and.(tab(i,iv).lt.vmax1)) then
0541 c                distri(2,j+1)=distri(2,j+1)+1.0*cn
0542                 distri(2,j+1)=distri(2,j+1)+1.0
0543                 tab(i,2)=1.0
0544                 ntt=ntt+1.0
0545                 goto 10
0546               else
0547                 goto 10
0548               endif
0549            else
0550               goto 10
0551            endif
0552    10    continue
0553    20 continue 
0554 c
0555 c      do 30 i=1,nv
0556 c   30 distri(2,i)=distri(2,i)*nt/ntt
0557 c      return
0558       end
0559 
0560 
0561 
0562       function ran2(idum)
0563       integer m,ia,ic,ir(97),iy,idum,j,iff
0564       real rm 
0565 c
0566 c     Variable aleatoire uniforme comprise entre 0 et 1
0567 c     ATTENTION : periodicite m=714025 -> ne pas tirer plus de m fois
0568 c
0569       m=714025
0570       ia=1366
0571       ic=150889
0572       iff=0
0573       idum=-idum
0574       rm=1./m
0575       if (idum.lt.0.or.iff.eq.0) then
0576          iff=1
0577          idum=mod(ic-idum,m)
0578          do 10 j=1,97
0579              idum=mod(ia*idum+ic,m)
0580              ir(j)=idum
0581    10        continue
0582          idum=mod(ia*idum+ic,m)
0583          iy=idum
0584       endif
0585       j=1+(97*iy)/m
0586       if(j.gt.97.or.j.lt.1) write(2,*) 'erreur random'
0587       iy=ir(j)
0588       ran2=iy*rm
0589       idum=mod(ia*idum+ic,m)
0590       ir(j)=idum
0591       return
0592       end
0593 
0594 
0595 
0596 
0597 
0598 
0599       function tirazi(idum)
0600       integer m,ia,ic,ir(97),iy,idum,j,iff
0601       real*4 rm,pi
0602       data pi/3.1415926535/
0603 c
0604 c     Variable aleatoire uniforme comprise entre 0 et 2*pi
0605 c     ATTENTION : periodicite m=714025 -> ne pas tirer plus de m fois
0606 c     (tirazi=2*pi*ran2 -> ran2 : variable aleatoire uniforme (0,1)
0607 c
0608       m=714025
0609       ia=1366
0610       ic=150889
0611       iff=0
0612       idum=-idum
0613       rm=1./m
0614       if (idum.lt.0.or.iff.eq.0) then
0615          iff=1
0616          idum=mod(ic-idum,m)
0617          do 10 j=1,97
0618              idum=mod(ia*idum+ic,m)
0619              ir(j)=idum
0620    10        continue
0621          idum=mod(ia*idum+ic,m)
0622          iy=idum
0623       endif
0624       j=1+(97*iy)/m
0625       if(j.gt.97.or.j.lt.1) write(2,*) 'erreur random'
0626       iy=ir(j)
0627       tirazi=iy*rm*2.0*pi
0628       idum=mod(ia*idum+ic,m)
0629       ir(j)=idum
0630       return
0631       end
0632 
0633 
0634 
0635 
0636 
0637       subroutine tirpos(mhut,var)
0638       real var,mhut,aleat
0639       integer reset
0640 c
0641       common/etiq1/reset
0642 c
0643       aleat=ran2(reset)
0644       if (aleat.ne.0.0) then
0645          var=-log(aleat)/mhut
0646       else
0647          write(2,*) 'Erreur tirpos'
0648 c         pause
0649       endif
0650       return
0651       end
0652 
0653 
0654 
0655       subroutine ener(hnu10,hnu0,hnu)
0656       real hnu10,hnu0,hnu,aleat
0657       integer reset
0658 c
0659       common/etiq1/reset
0660 c
0661       aleat=ran2(reset)
0662       if (aleat.ne.0.0) then
0663          hnu=hnu10-hnu0*log(1.0-aleat)
0664       else
0665          write(2,*) 'Erreur tirener'
0666 c         pause
0667       endif
0668       return
0669       end
0670 
0671 
0672       subroutine posi(a,x)
0673       real a,x,aleat
0674       integer reset
0675 c
0676       common/etiq1/reset
0677 c
0678       aleat=ran2(reset)
0679       x=a*sqrt(aleat)
0680       return
0681       end
0682 
0683 
0684       subroutine tirsur(survie,mhuph,mhuco)
0685       real mhuco,mhuph,r,va
0686       integer reset,survie
0687 c
0688       common/etiq1/reset
0689 c
0690       r=mhuco/(mhuco+mhuph)
0691 c
0692       va=ran2(reset)
0693 c      write(1,*) va
0694       if (va.le.r) then
0695         survie=1
0696       else
0697         survie=0
0698       endif
0699       return
0700       end
0701 
0702 
0703       function tircom(hnu0)
0704       real q1,q2,q3,lambda,r,hnu0
0705       real test1,test2a,test2b
0706       integer reset
0707 c
0708       common/etiq1/reset
0709 c
0710       data mec2/511.0/
0711 c
0712 c     Variable aleatoire associee a l'effet Compton
0713 c
0714    10 continue
0715       q1=ran2(reset)
0716       q2=ran2(reset)    
0717       q3=ran2(reset)
0718       lambda=mec2/hnu0
0719       test1=(1.0+2.0/lambda)/(9.0+2.0/lambda)
0720       if (q1.le.test1) then
0721           r=1.0+2.0*q2/lambda
0722           test2a=4.0*(1.0/r-1.0/r**2)
0723           if (q3.le.test2a) then
0724              tircom=1.0/r
0725           else 
0726              goto 10
0727           endif
0728       else 
0729           r=(1.0+2.0/lambda)/(1.0+2.0*q2/lambda)
0730           test2b=0.5*((lambda-r*lambda+1.0)**2+1.0/r)
0731           if (q3.le.test2b) then
0732              tircom=1.0/r
0733           else 
0734              goto10
0735           endif
0736       endif
0737       return
0738       end
0739 
0740       subroutine xLYSO(Eph,mhuph,mhuco)
0741       real Eph,mhuph,mhuco
0742       real mhuphlu,mhuphyt,mhuphsi,mhuphoxy
0743       real mhucolu,mhucoyt,mhucosi,mhucooxy
0744       data dens,denslu,densyt/7.1,9.842,4.405/
0745       data denssi,densoxy/2.33,1.43e-3/
0746       data plu,pyt,psi,poxy/0.71,0.0403,0.0637,0.1815/
0747 c      if ((Eph.ge.3.99).and.(Eph.le.1000.0)) then
0748       if (Eph.ge.3.99) then
0749         call xlu(Eph,mhuphlu,mhucolu)
0750         call xyt(Eph,mhuphyt,mhucoyt)
0751         call xsi(Eph,mhuphsi,mhucosi)
0752         call xoxy(Eph,mhuphoxy,mhucooxy)
0753         mhuph=mhuphlu*plu*dens/denslu
0754         mhuph=mhuph+mhuphyt*pyt*dens/densyt
0755         mhuph=mhuph+mhuphsi*psi*dens/denssi
0756         mhuph=mhuph+mhuphoxy*poxy*dens/densoxy
0757         mhuco=mhucolu*plu*dens/denslu
0758         mhuco=mhuco+mhucoyt*pyt*dens/densyt
0759         mhuco=mhuco+mhucosi*psi*dens/denssi
0760         mhuco=mhuco+mhucooxy*poxy*dens/densoxy
0761       else
0762         write(2,*) 'LYSO (Lu2(1-x)Y2xSiO5,x = 0.1): k<3.99 ou k>1000'
0763       endif
0764       return
0765       end      
0766 
0767       subroutine xlu(Eph,mhup,mhuc)
0768       real Eph
0769       real mhu0,mhup,mhuc
0770       data dens,pa,avo,barn/9.842,174.99,6.023e+23,1.0e-24/
0771       data ap10,ap11,ap12,ap13/12.6387,1.55476,-0.881094,0.0602036/
0772       data ap20,ap21,ap22,ap23/17.2638,-2.37189,-0.0495994,0.0/
0773       data ap30,ap31,ap32,ap33/16.2289,-2.67128,0.0,0.0/
0774       data ap40,ap41,ap42,ap43/13.9813,-2.40841,0.0,0.0/
0775       data ac0,ac1,ac2,ac3/0.197176,0.150264,-0.192474,0.00385751/
0776       mhu0=dens*avo*barn/pa
0777       if (Eph.lt.63.314) then
0778         mhup=exp(ap10+ap11*log(Eph)+ap12*log(Eph)**2+ap13*log(Eph)**3)
0779       elseif (Eph.lt.10.8700) then
0780         mhup=exp(ap20+ap21*log(Eph)+ap22*log(Eph)**2+ap23*log(Eph)**3)
0781       elseif (Eph.lt.2.4920) then
0782         mhup=exp(ap30+ap31*log(Eph)+ap32*log(Eph)**2+ap33*log(Eph)**3)
0783       else
0784         mhup=exp(ap40+ap41*log(Eph)+ap42*log(Eph)**2+ap43*log(Eph)**3)
0785       endif
0786       mhuc=exp(ac0+ac1*log(Eph)+ac2*log(Eph)**2+ac3*log(Eph)**3)
0787       mhup=mhup*mhu0
0788       mhuc=mhuc*mhu0
0789       return
0790       end
0791 
0792       subroutine xyt(Eph,mhup,mhuc)
0793       real Eph
0794       real mhu0,mhup,mhuc
0795       data dens,pa,avo,barn/4.405,88.9050,6.023e+23,1.0e-24/
0796       data ap10,ap11,ap12,ap13/13.4674,0.191023,-0.686616,0.0497356/
0797       data ap20,ap21,ap22,ap23/15.1822,-2.38946,-0.0881174,0.0/
0798       data ap30,ap31,ap32,ap33/13.2775,-2.43174,0.0,0.0/
0799       data ac0,ac1,ac2,ac3/0.0629057,  1.41577, -0.199713,  0.00533312/
0800       mhu0=dens*avo*barn/pa
0801       if (Eph.lt.17.08) then
0802         mhup=exp(ap10+ap11*log(Eph)+ap12*log(Eph)**2+ap13*log(Eph)**3)
0803       elseif (Eph.lt.2.373) then
0804         mhup=exp(ap20+ap21*log(Eph)+ap22*log(Eph)**2+ap23*log(Eph)**3)
0805       else
0806         mhup=exp(ap30+ap31*log(Eph)+ap32*log(Eph)**2+ap33*log(Eph)**3)
0807       endif
0808       mhuc=exp(ac0+ac1*log(Eph)+ac2*log(Eph)**2+ac3*log(Eph)**3)
0809       mhup=mhup*mhu0
0810       mhuc=mhuc*mhu0
0811       return
0812       end
0813 
0814       subroutine xsi(Eph,mhup,mhuc)
0815       real Eph
0816       real mhu0,mhup,mhuc
0817       data dens,pa,avo,barn/2.33,28.086,6.023e+23,1.0e-24/
0818       data ap10,ap11,ap12,ap13/13.2682,-1.98174,-0.31695,0.0273928/
0819       data ac0,ac1,ac2,ac3/-0.414971,1.34868,-0.222315,0.00841959/
0820       mhu0=dens*avo*barn/pa
0821       mhup=exp(ap10+ap11*log(Eph)+ap12*log(Eph)**2+ap13*log(Eph)**3)
0822       mhuc=exp(ac0+ac1*log(Eph)+ac2*log(Eph)**2+ac3*log(Eph)**3)
0823       mhup=mhup*mhu0
0824       mhuc=mhuc*mhu0
0825       return
0826       end
0827 
0828       subroutine xCdTe(Eph,mhuph,mhuco)
0829       real Eph,mhuph,mhuco
0830       real mhuphCd,mhuphTe
0831       real mhucoCd,mhucoTe
0832       data dens,densCd,densTe/5.856,8.65,6.24/
0833       data pCd,pTe/0.468,0.532/
0834 c      if ((Eph.ge.3.99).and.(Eph.le.1000.0)) then
0835       if (Eph.ge.3.99) then
0836         call xCd(Eph,mhuphCd,mhucoCd)
0837         call xTe(Eph,mhuphTe,mhucoTe)
0838         mhuph=mhuphCd*pCd*dens/densCd
0839         mhuph=mhuph+mhuphTe*pTe*dens/densTe
0840         mhuco=mhucoCd*pCd*dens/densCd
0841         mhuco=mhuco+mhucoTe*pTe*dens/densTe
0842       else
0843         write(2,*) 'CdTe: k<3.99 ou k>1000'
0844       endif
0845       return
0846       end      
0847       
0848       
0849       subroutine xCd(Eph,mhup,mhuc)
0850       real Eph
0851       real mhu0,mhup,mhuc
0852       data dens,pa,avo,barn/8.65,112.41,6.023e+23,1.0e-24/
0853       data ap10,ap11,ap12,ap13/15.9668,-2.38363,-0.0801104,0.0/
0854       data ap20,ap21,ap22,ap23/11.5254,1.07714,-0.831424,0.0579120/
0855       data ac0,ac1,ac2,ac3/-0.0516701,1.57426,-0.227646,0.0070565/
0856       mhu0=dens*avo*barn/pa
0857       if (Eph.lt.26.711) then
0858         mhup=exp(ap10+ap11*log(Eph)+ap12*log(Eph)**2+ap13*log(Eph)**3)
0859       else
0860         mhup=exp(ap20+ap21*log(Eph)+ap22*log(Eph)**2+ap23*log(Eph)**3)
0861       endif
0862       mhuc=exp(ac0+ac1*log(Eph)+ac2*log(Eph)**2+ac3*log(Eph)**3)
0863       mhup=mhup*mhu0
0864       mhuc=mhuc*mhu0
0865       return
0866       end
0867       
0868       subroutine xTe(Eph,mhup,mhuc)
0869       real Eph
0870       real mhu0,mhup,mhuc
0871       data dens,pa,avo,barn/6.24,127.6,6.023e+23,1.0e-24/
0872       data ap10,ap11,ap12,ap13/16.1087,-2.27876,-0.0929405,0.0/
0873       data ap20,ap21,ap22,ap23/11.6656,1.71052,-0.948281,0.0653213/
0874       data ac0,ac1,ac2,ac3/-0.0407579,1.64267,-0.247897,0.00880567/
0875       mhu0=dens*avo*barn/pa
0876       if (Eph.lt.31.813) then
0877         mhup=exp(ap10+ap11*log(Eph)+ap12*log(Eph)**2+ap13*log(Eph)**3)
0878       else
0879         mhup=exp(ap20+ap21*log(Eph)+ap22*log(Eph)**2+ap23*log(Eph)**3)
0880       endif
0881       mhuc=exp(ac0+ac1*log(Eph)+ac2*log(Eph)**2+ac3*log(Eph)**3)
0882       mhup=mhup*mhu0
0883       mhuc=mhuc*mhu0
0884       return
0885       end
0886 
0887       subroutine xBGO(Eph,mhuph,mhuco)
0888       real Eph,mhuph,mhuco
0889       real mhuphbismu,mhuphgerma,mhuphoxy
0890       real mhucobismu,mhucogerma,mhucooxy
0891       data dens,densbismu,densgerma,densoxy/7.13,9.8,5.32,1.43e-3/
0892       data pbismu,pgerma,poxy/0.671,0.175,0.154/
0893 c     if ((Eph.ge.3.99).and.(Eph.le.1000.0)) then
0894       if (Eph.ge.3.99) then
0895         call xbismu(Eph,mhuphbismu,mhucobismu)
0896         call xgerma(Eph,mhuphgerma,mhucogerma)
0897         call xoxy(Eph,mhuphoxy,mhucooxy)
0898         mhuph=mhuphbismu*pbismu*dens/densbismu
0899         mhuph=mhuph+mhuphgerma*pgerma*dens/densgerma
0900         mhuph=mhuph+mhuphoxy*poxy*dens/densoxy
0901         mhuco=mhucobismu*pbismu*dens/densbismu
0902         mhuco=mhuco+mhucogerma*pgerma*dens/densgerma
0903         mhuco=mhuco+mhucooxy*poxy*dens/densoxy
0904       else
0905         write(2,*) 'BGO: k<3.99 ou k>1000'
0906       endif
0907       return
0908       end
0909 
0910 
0911       subroutine xoxy(Eph,mhup,mhuc)
0912       real Eph
0913       real mhu0,mhup,mhuc
0914       data dens,pa,avo,barn/1.43e-3,16.00,6.023e+23,1.0e-24/
0915       data ap0,ap1,ap2,ap3/11.7130,-2.57229,-0.205893,0.0199244/
0916       data ac0,ac1,ac2,ac3/-1.73679,2.17686,-0.449050,0.0264733/
0917       mhu0=dens*avo*barn/pa
0918       mhup=exp(ap0+ap1*log(Eph)+ap2*log(Eph)**2+ap3*log(Eph)**3)
0919       mhuc=exp(ac0+ac1*log(Eph)+ac2*log(Eph)**2+ac3*log(Eph)**3)
0920       mhup=mhup*mhu0
0921       mhuc=mhuc*mhu0
0922       return
0923       end
0924 
0925 
0926       subroutine xbismu(Eph,mhup,mhuc)
0927       real Eph
0928       real mhu0,mhup,mhuc
0929       data dens,pa,avo,barn/9.8,209.00,6.023e+23,1.0e-24/
0930       data ap10,ap11,ap12,ap13/16.7078,-2.58648,0.0,0.0/
0931       data ap20,ap21,ap22,ap23/17.5348,-2.23353,-0.0596161,0.0/
0932       data ap30,ap31,ap32,ap33/9.44293,3.44965,-1.19886,0.0783484/
0933       data ac0,ac1,ac2,ac3/0.189860,1.56125,-0.200932,0.00436768/
0934       mhu0=dens*avo*barn/pa
0935       if (Eph.lt.16.39) then
0936         mhup=exp(ap10+ap11*log(Eph)+ap12*log(Eph)**2+ap13*log(Eph)**3)
0937       elseif ((Eph.ge.16.39).and.(Eph.lt.90.53)) then
0938         mhup=exp(ap20+ap21*log(Eph)+ap22*log(Eph)**2+ap23*log(Eph)**3)
0939       else
0940         mhup=exp(ap30+ap31*log(Eph)+ap32*log(Eph)**2+ap33*log(Eph)**3)
0941       endif
0942       mhuc=exp(ac0+ac1*log(Eph)+ac2*log(Eph)**2+ac3*log(Eph)**3)
0943       mhup=mhup*mhu0
0944       mhuc=mhuc*mhu0
0945       return
0946       end
0947 
0948 
0949 
0950       subroutine xgerma(Eph,mhup,mhuc)
0951       real Eph
0952       real mhu0,mhup,mhuc
0953       data dens,pa,avo,barn/5.32,72.6,6.023e+23,1.0e-24/
0954       data ap10,ap11,ap12,ap13/14.6813,-2.59285,-0.0208355,0.0/
0955       data ap20,ap21,ap22,ap23/13.9288,-0.479613,-0.572897,0.0431277/
0956       data ac0,ac1,ac2,ac3/ -0.334383,1.60237,-0.245555,0.00871239/
0957       mhu0=dens*avo*barn/pa
0958       if (Eph.lt.11.10) then
0959         mhup=exp(ap10+ap11*log(Eph)+ap12*log(Eph)**2+ap13*log(Eph)**3)
0960       else
0961         mhup=exp(ap20+ap21*log(Eph)+ap22*log(Eph)**2+ap23*log(Eph)**3)
0962       endif
0963       mhuc=exp(ac0+ac1*log(Eph)+ac2*log(Eph)**2+ac3*log(Eph)**3)
0964       mhup=mhup*mhu0
0965       mhuc=mhuc*mhu0
0966       return
0967       end
0968 
0969       function gasdev(type,k)
0970       real sigma,k
0971       integer reset,type
0972 c
0973       common/etiq1/reset
0974 c
0975       if (type.eq.1) then
0976          sigma=4.5338*k**0.48/(2.0*sqrt(2.0*log(2.0)))
0977       elseif (type.eq.2) then
0978          sigma=(2.315*sqrt(k)+0.01*k)/(2.0*sqrt(2.0*log(2.0)))
0979       elseif (type.eq.3) then
0980          sigma=13.12*k**0.2051/(2.0*sqrt(2.0*log(2.0)))
0981       elseif (type.eq.4) then
0982          sigma=0
0983       elseif (type.eq.5) then
0984          sigma = 0
0985       elseif (type.eq.6) then
0986          sigma = 0
0987       endif
0988 c
0989     1 v1=2.0*ran2(reset)-1
0990       v2=2.0*ran2(reset)-1
0991       r=v1**2+v2**2
0992       if (r.ge.1.0) goto 1
0993       fac=sqrt(-2*log(r)/r)
0994       gasdev=v2*fac*sigma
0995       return
0996       end
0997 
0998 
0999       subroutine xNaI(Eph,mhuph,mhuco)
1000       real Eph,mhuph,mhuco
1001       real mhucoNa,mhuphNa,mhucoI,mhuphI
1002       data dens,densNa,densI/3.67,0.97,4.94/
1003       data pNa,pI/0.1534,0.8466/
1004 c      if ((Eph.ge.5.198).and.(Eph.le.1000.0)) then
1005       if (Eph.ge.5.198) then
1006          call xNa(Eph,mhuphNa,mhucoNa)
1007          call xI(Eph,mhuphI,mhucoI)
1008          mhuph=mhuphNa*pNa*dens/densNa
1009          mhuph=mhuph+mhuphI*pI*dens/densI
1010          mhuco=mhucoNa*pNa*dens/densNa
1011          mhuco=mhuco+mhucoI*pI*dens/densI
1012       else
1013         write(2,*) 'NaI: k<5.198 ou k>1000'
1014       endif
1015       return
1016       end
1017 
1018       subroutine xCsI(Eph,mhuph,mhuco)
1019       real Eph,mhuph,mhuco
1020       real mhuphCs,mhucoCs,mhuphI,mhucoI
1021       data dens,densCs,densI/4.51,1.873,4.94/
1022       data pCs,pI/0.5115,0.4885/
1023 c      if ((Eph.ge.5.713).and.(Eph.le.1000.0)) then
1024       if (Eph.ge.5.713) then
1025          call xCs(Eph,mhuphCs,mhucoCs)
1026          call xI(Eph,mhuphI,mhucoI)
1027          mhuph=mhuphCs*pCs*dens/densCs
1028          mhuph=mhuph+mhuphI*pI*dens/densI
1029          mhuco=mhucoCs*pCs*dens/densCs
1030          mhuco=mhuco+mhucoI*pI*dens/densI
1031       else
1032         write(2,*) 'CsI: k<5.713 ou k>1000'
1033       endif
1034       return
1035       end
1036 
1037 
1038 
1039       subroutine xCs(Eph,mhup,mhuc)
1040       real Eph
1041       real mhu0,mhup,mhuc
1042       data dens,pa,avo,barn/1.873,132.91,6.023e+23,1.0e-24/
1043       data ap10,ap11,ap12,ap13/16.5418,-2.46363,-0.0542849,0.0/
1044       data ap20,ap21,ap22,ap23/11.3757,1.94161,-0.983232,0.0671986/
1045       data ac0,ac1,ac2,ac3/0.184861,1.50030,-0.213333,0.00624264/
1046       mhu0=dens*avo*barn/pa
1047       if (Eph.lt.35.985) then
1048          mhup=exp(ap10+ap11*log(Eph)+ap12*log(Eph)**2+ap13*log(Eph)**3)
1049       else
1050          mhup=exp(ap20+ap21*log(Eph)+ap22*log(Eph)**2+ap23*log(Eph)**3)
1051       endif
1052       mhuc=exp(ac0+ac1*log(Eph)+ac2*log(Eph)**2+ac3*log(Eph)**3)
1053       mhup=mhu0*mhup
1054       mhuc=mhu0*mhuc
1055       return
1056       end
1057 
1058 
1059       subroutine xNa(Eph,mhup,mhuc)
1060       real Eph
1061       real mhu0,mhup,mhuc
1062       data dens,pa,avo,barn/0.97,22.997,6.023e+23,1.0e-24/
1063       data ap10,ap11,ap12,ap13/12.6777,-2.24521,-0.274873,0.02507/
1064       data ac0,ac1,ac2,ac3/-0.967717,1.61794,-0.287191,0.0131526/
1065       mhu0=dens*avo*barn/pa
1066       mhup=exp(ap10+ap11*log(Eph)+ap12*log(Eph)**2+ap13*log(Eph)**3)
1067       mhuc=exp(ac0+ac1*log(Eph)+ac2*log(Eph)**2+ac3*log(Eph)**3)
1068       mhup=mhu0*mhup
1069       mhuc=mhu0*mhuc
1070       return
1071       end
1072 
1073 
1074 
1075       subroutine xI(Eph,mhup,mhuc)
1076       real Eph
1077       real mhu0,mhup,mhuc
1078       data dens,pa,avo,barn/4.94,126.91,6.023e+23,1.0e-24/
1079       data ap10,ap11,ap12,ap13/16.4086,-2.48214,-0.0507179,0.0/
1080       data ap20,ap21,ap22,ap23/12.1075,1.43635,-0.882038,0.0603575/
1081       data ac0,ac1,ac2,ac3/-0.0404420,1.65596,-0.251067,0.00904874/
1082       mhu0=dens*avo*barn/pa
1083       if (Eph.lt.33.169) then
1084          mhup=exp(ap10+ap11*log(Eph)+ap12*log(Eph)**2+ap13*log(Eph)**3)
1085       else
1086          mhup=exp(ap20+ap21*log(Eph)+ap22*log(Eph)**2+ap23*log(Eph)**3)
1087       endif
1088       mhuc=exp(ac0+ac1*log(Eph)+ac2*log(Eph)**2+ac3*log(Eph)**3)
1089       mhuc=mhu0*mhuc
1090       mhup=mhu0*mhup
1091       return
1092       end