*********************************************************************** * SUBROUTINE REFLAST * * * MODIFIED BY: DATE:09/05/93 17.33.11 * * * *********************************************************************** implicit none cSEQ,DRPARS,GUPSI,GUREF,GUMAT,GUFIL,GUCONV. #include "gurefp.inc" #include "guref.inc" #include "gumat.inc" #include "gupsi.inc" #include "gufil.inc" real tmatl(8,8),evecl(8,8),evall(8),work(100),d1s,d1p,d1e,dlo,dhi integer icycle,i1,i2,i3,i4,i5,inm,in,imax,ifo,irl,istat,it,ip character*24 file1,file2,file3 data file1/'file1'/ data file2/'file2'/ data file3/'file3'/ c output of reference tracks if(nref(1).eq.1)then open(unit=32,file=creff,form='formatted',status='new') c open(unit=32,file=file1,form='formatted',status='new') write(32,*)(pmcr(i1),i1=1,20) write(32,*)(iref0(i1),i1=1,5) write(32,*)x1pos write(32,*)y1pos write(32,*)x2pos write(32,*)y2pos write(32,*)x3pos write(32,*)y3pos write(32,*)x4pos write(32,*)y4pos close(unit=32) open(unit=33,file=csigf,form='formatted',status='new') c open(unit=33,file=file1,form='formatted',status='new') write(33,*)(pmcr(i1),i1=1,20) write(33,*)(iref0(i1),i1=1,5) write(33,*)sx1pos write(33,*)sy1pos write(33,*)sx2pos write(33,*)sy2pos write(33,*)sx3pos write(33,*)sy3pos write(33,*)sx4pos write(33,*)sy4pos close(unit=33) end if ************************************************************************ if(nref(1).eq.2.or.nref(1).eq.3)then do i4=1,nkik do i3=1,npad if(tmultp(i3,i4).ne.0)then d1p=1./float(tmultp(i3,i4)) else d1p=0. end if if(tmulte(i3,i4).ne.0)then d1e=1./float(tmulte(i3,i4)) else d1e=0. end if if(tmults(i3,i4).ne.0)then d1s=1./float(tmults(i3,i4)) else d1s=0. end if do i1=1,8 tcents(i1,i3,i4)=d1s*tcents(i1,i3,i4) tcentp(i1,i3,i4)=d1p*tcentp(i1,i3,i4) tcente(i1,i3,i4)=d1e*tcente(i1,i3,i4) do i2=1,8 tmats(i1,i2,i3,i4)=d1s*tmats(i1,i2,i3,i4) tmatp(i1,i2,i3,i4)=d1p*tmatp(i1,i2,i3,i4) tmate(i1,i2,i3,i4)=d1e*tmate(i1,i2,i3,i4) end do end do end do end do do i4=1,nkik do i3=1,npad do i1=1,8 do i2=1,8 tmats(i1,i2,i3,i4)=tmats(i1,i2,i3,i4)- 1 tcents(i1,i3,i4)*tcents(i2,i3,i4) tmatp(i1,i2,i3,i4)=tmatp(i1,i2,i3,i4)- 1 tcentp(i1,i3,i4)*tcentp(i2,i3,i4) tmate(i1,i2,i3,i4)=tmate(i1,i2,i3,i4)- 1 tcente(i1,i3,i4)*tcente(i2,i3,i4) end do end do end do end do ************************************************************************ do i4=1,nkik do i3=1,npad do i1=1,8 do i2=1,8 tmatl(i1,i2)=tmatp(i1,i2,i3,i4) end do end do inm=8 in=8 dlo=-100000. dhi=100000. imax=8 ifo=0 call eisrs3(inm,in,tmatl,dlo,dhi,imax,ifo,evall,evecl,irl,work) do i1=1,8 do i2=1,8 evecp(i1,i2,i3,i4)=evecl(i1,i2) end do evalp(i1,i3,i4)=evall(i1) end do end do end do ************************************************************************ do i4=1,nkik do i3=1,npad do i1=1,8 do i2=1,8 tmatl(i1,i2)=tmate(i1,i2,i3,i4) end do end do inm=8 in=8 dlo=-100000. dhi=100000. imax=8 ifo=0 call eisrs3(inm,in,tmatl,dlo,dhi,imax,ifo,evall,evecl,irl,work) do i1=1,8 do i2=1,8 evece(i1,i2,i3,i4)=evecl(i1,i2) end do evale(i1,i3,i4)=evall(i1) end do end do end do ************************************************************************ do i4=1,nkik do i3=1,npad do i1=1,8 do i2=1,8 tmatl(i1,i2)=tmats(i1,i2,i3,i4) end do end do inm=8 in=8 dlo=-100000. dhi=100000. imax=8 ifo=0 call eisrs3(inm,in,tmatl,dlo,dhi,imax,ifo,evall,evecl,irl,work) do i1=1,8 do i2=1,8 evecs(i1,i2,i3,i4)=evecl(i1,i2) end do evals(i1,i3,i4)=evall(i1) end do end do end do ************************************************************************ * begin of analysis * if it is a positron then use evecp instead of evecs etc * if it is a electron then use evece instead of evecs etc * know that track candidate is in pad it * do i1=1,8 * d1=0. * do i2=1,8 * d1=d1+evecs(i2,i1,it)*(xdet(i2)-tcents(i2,it)) * end do * psi(i1)=d1 * end do * * check whether psi(1) is within +- 3,4,5 psisigs(1,it) * check whether psi(2) is within +- 3,4,5 psisigs(2,it) * check whether psi(3) is within +- 3,4,5 psisigs(3,it) * do i4=1,nkik do it=1,npad do ip=1,tmultp(it,i4) do i1=1,8 d1p=0. d1s=0. do i2=1,8 d1p=d1p+evecp(i2,i1,it,i4)*(detp(i2,ip,it,i4) 1 -tcentp(i2,it,i4)) d1s=d1s+evecs(i2,i1,it,i4)*(detp(i2,ip,it,i4) 1 -tcents(i2,it,i4)) end do psisp(i1,it,i4)=psisp(i1,it,i4)+d1p psis2p(i1,it,i4)=psis2p(i1,it,i4)+d1p**2 psiss(i1,it,i4)=psiss(i1,it,i4)+d1s psis2s(i1,it,i4)=psis2s(i1,it,i4)+d1s**2 i3=10*(it-1)+i1 call hfill(i3,d1p,0.,1.) i3=10*(it-1)+i1+200 call hfill(i3,d1s,0.,1.) end do end do end do end do do i4=1,nkik do it=1,npad do ip=1,tmulte(it,i4) do i1=1,8 d1e=0. d1s=0. do i2=1,8 d1e=d1e+evece(i2,i1,it,i4)*(dete(i2,ip,it,i4) 1 -tcente(i2,it,i4)) d1s=d1s+evecs(i2,i1,it,i4)*(dete(i2,ip,it,i4) 1 -tcents(i2,it,i4)) end do psise(i1,it,i4)=psise(i1,it,i4)+d1e psis2e(i1,it,i4)=psis2e(i1,it,i4)+d1e**2 psiss(i1,it,i4)=psiss(i1,it,i4)+d1s psis2s(i1,it,i4)=psis2s(i1,it,i4)+d1s**2 i3=10*(it-1)+i1+100 call hfill(i3,d1e,0.,1.) i3=10*(it-1)+i1+200 call hfill(i3,d1s,0.,1.) end do end do end do end do do i4=1,nkik do i2=1,npad if(tmultp(i2,i4).ne.0)then d1p=1./float(tmultp(i2,i4)) else d1p=0. end if if(tmulte(i2,i4).ne.0)then d1e=1./float(tmulte(i2,i4)) else d1e=0. end if if(tmults(i2,i4).ne.0)then d1s=1./float(tmults(i2,i4)) else d1s=0. end if do i1=1,8 psisp(i1,i2,i4)=d1p*psisp(i1,i2,i4) psise(i1,i2,i4)=d1e*psise(i1,i2,i4) psiss(i1,i2,i4)=d1s*psiss(i1,i2,i4) psis2p(i1,i2,i4)=d1p*psis2p(i1,i2,i4) psis2e(i1,i2,i4)=d1e*psis2e(i1,i2,i4) psis2s(i1,i2,i4)=d1s*psis2s(i1,i2,i4) end do end do end do do i4=1,nkik do i2=1,npad do i1=1,8 psisigp(i1,i2,i4)=sqrt(psis2p(i1,i2,i4)-psisp(i1,i2,i4)**2) psisige(i1,i2,i4)=sqrt(psis2e(i1,i2,i4)-psise(i1,i2,i4)**2) psisigs(i1,i2,i4)=sqrt(psis2s(i1,i2,i4)-psiss(i1,i2,i4)**2) end do end do end do open(unit=32,file=cmatf,form='formatted',status='new') c open(unit=32,file=file2,form='formatted',status='new') write(32,*)(dkik0(i1),i1=1,nkik) do i4=1,nkik write(32,*)(((tmatp(i1,i2,i3,i4),i3=1,npad),i2=1,8),i1=1,8) write(32,*)(((evecp(i1,i2,i3,i4),i3=1,npad),i2=1,8),i1=1,8) write(32,*)((tcentp(i1,i2,i4),i2=1,npad),i1=1,8) write(32,*)((evalp(i1,i2,i4),i2=1,npad),i1=1,8) write(32,*)(tmultp(i1,i4),i1=1,npad) c write(32,*)((psisp(i1,i2,i4),i1=1,8),i2=1,npad) c write(32,*)((psis2p(i1,i2,i4),i1=1,8),i2=1,npad) c write(32,*)((psisigp(i1,i2,i4),i1=1,8),i2=1,npad) write(32,*)(((tmate(i1,i2,i3,i4),i3=1,npad),i2=1,8),i1=1,8) write(32,*)(((evece(i1,i2,i3,i4),i3=1,npad),i2=1,8),i1=1,8) write(32,*)((tcente(i1,i2,i4),i2=1,npad),i1=1,8) write(32,*)((evale(i1,i2,i4),i2=1,npad),i1=1,8) write(32,*)(tmulte(i1,i4),i1=1,npad) c write(32,*)((psise(i1,i2,i4),i1=1,8),i2=1,npad) c write(32,*)((psis2e(i1,i2,i4),i1=1,8),i2=1,npad) c write(32,*)((psisige(i1,i2,i4),i1=1,8),i2=1,npad) write(32,*)(((tmats(i1,i2,i3,i4),i3=1,npad),i2=1,8),i1=1,8) write(32,*)(((evecs(i1,i2,i3,i4),i3=1,npad),i2=1,8),i1=1,8) write(32,*)((tcents(i1,i2,i4),i2=1,npad),i1=1,8) write(32,*)((evals(i1,i2,i4),i2=1,npad),i1=1,8) write(32,*)(tmults(i1,i4),i1=1,npad) c write(32,*)((psiss(i1,i2,i4),i1=1,8),i2=1,npad) c write(32,*)((psis2s(i1,i2,i4),i1=1,8),i2=1,npad) c write(32,*)((psisigs(i1,i2,i4),i1=1,8),i2=1,npad) end do close(unit=32) call hropen(32,'hist',cmath,'N',1024,istat) c call hropen(32,'hist',file3,'N',1024,istat) c open(unit=32,file=crzhf,form='unformatted',recl=1024, c 1 access='direct',status='new',iostat=istat) c call hrfile(32,'hist','N') call hrout(951,icycle,' ') call hrout(952,icycle,' ') call hrout(953,icycle,' ') call hrout(954,icycle,' ') call hrout(955,icycle,' ') call hrout(956,icycle,' ') call hrout(957,icycle,' ') call hrout(958,icycle,' ') call hrout(961,icycle,' ') call hrout(962,icycle,' ') call hrout(963,icycle,' ') call hrout(964,icycle,' ') call hrout(965,icycle,' ') call hrout(966,icycle,' ') call hrout(967,icycle,' ') call hrout(968,icycle,' ') call hrout(971,icycle,' ') call hrout(972,icycle,' ') call hrout(973,icycle,' ') call hrout(974,icycle,' ') call hrout(975,icycle,' ') call hrout(976,icycle,' ') call hrout(977,icycle,' ') call hrout(978,icycle,' ') call hrout(979,icycle,' ') call hrout(980,icycle,' ') call hrout(981,icycle,' ') call hrout(982,icycle,' ') call hrout(983,icycle,' ') call hrout(984,icycle,' ') call hrout(985,icycle,' ') call hrout(986,icycle,' ') call hrout(987,icycle,' ') call hrout(988,icycle,' ') call hrout(989,icycle,' ') call hrout(990,icycle,' ') call hrend('hist') close(unit=32) end if ************************************************************************ RETURN END