subroutine fpp_track(icode) c c Written by Ed Brash - January 13, 1995 - packs data from raw c data structures and organizes the calls to the actual tracking c routines (trk4.f and associated routines). Fills most FPP c related histograms. c c Modified by Ed Brash - January-February, 1996 - accumulate c statistics for time offset and alignment calibration routines. c c implicit none c include 'fpp_local.h' c include 'parameter.h' include 'espace_type.h' include 'espace_user.h' include 'detector.h' include 'transport.h' include 'option.h' c integer*4 ntrack !# of tracks integer*4 itrk(5) !# of tracks found at each tracking interation integer*4 icode !signal to main code of tracking problem 0=o.k. integer*4 itrcode !tracking code integer*4 ntr(4) integer*4 ntrack_front,ntrack_rear,ntrack_tot integer*4 prev_straw real*8 track(40,25) !track information real*8 zplane,str,ddist,uvx,a,b,c,d,ee,ff,gg,hh,iii,ztr,yw real*8 jjj,dd2by2,dd2,dd,resres,fac real*8 residual,rbig real*8 t1,t2,t3,t4,t5,t6 real*8 a0_car,b1_car,b2_car,c1_car,c2_car,c3_car real*8 a0_CH2,b1_CH2,b2_CH2,c1_CH2,c2_CH2,c3_CH2 real*8 fg real*8 xfaold,xrold,yfaold,yrold,xfbold,yfbold real*8 tfold,trold,pfold,prold real*8 psif,psir,xhatr,yhatr,zhatr,xhatf,yhatf,zhatf real*8 xhat1f,yhat1f,zhat1f,xhat1r,yhat1r,zhat1r real*8 thetasc_car,phisc_car,psisc_car,thetafd,phifd real*8 theta_sc_CH2,phi_sc_CH2,theta_sc_car,phi_sc_car real*8 psia,psisc_CH2,thetasc_CH2,phisc_CH2 real*8 x_sc_car,y_sc_car,z_sc_car,psi_sc_car,r_sc_car real*8 x_sc_CH2,y_sc_CH2,z_sc_CH2,r_sc_CH2,psi_sc_CH2 c real*8 zone(25),xone(25),dxone(25),dxfone(25) integer*4 uone(25),lrfit(25) real*8 vat0,uat0 real*8 vchi, vslope ,uchi, uslope integer tr_num c real*8 xf3,xf4,yf3,yf4 real*8 xr3,xr4,yr3,yr4 real*8 dcone3,dcone4 c real*8 xf1,xf2,yf1,yf2 real*8 xr1,xr2,yr1,yr2 real*8 dcone1,dcone2 c real*8 usave(6),zusave(6),dusave(6) real*8 us,vs,zs1,vs1,vs3,us1,us2,us3,zs2 real*8 psi_tr integer*4 uplsave(6),vplsave(6),npl,npl_skip real*8 vsave(6),zvsave(6),dvsave(6) c real*8 ave1,ave2,zave1,zave2,npt1,npt2 real*8 ang_cr,int_cr real*8 chis(25) c real*8 xatz0f,yatz0f,xatz0r,yatz0r,xatz0x real*8 xtr(24),ytr(24),uvxtr(24),z(24) real*8 uvx_low,uvx_hi,strw_rad,uvx0 integer*4 straw_pred(24),nstraws,itrgood(24) real*8 tthetaf,tphif,tthetaa,tphia,thinitd,phiinitd real*8 tthetar,tphir real*8 theta_az_car,phi_az_car,theta_az_CH2,phi_az_CH2 real*8 xdiff_front,ydiff_front,xdiff_rear,ydiff_rear,xdiff_x real*8 theta_difff,phi_difff,theta_diffr,phi_diffr integer*4 jj,ii,kk,ntothits,nhinit,ngoodhit real*8 lambda1_CH2,lambda2_CH2,zclose_CH2,sclose_CH2 real*8 lambda1_car,lambda2_car,zclose_car,sclose_car integer*4 iwiregroup,igroup,itime real*8 letime,tetime,tdiff real*8 diff,rstraw,drdist integer*4 iss_front(8),imark,iss_rear_v3(8) integer*4 iss_rear_x3(8),iss_rear_u3(8),index integer*4 iss_rear_u4(8),iss_rear_v4(8),i data iss_front / 1,3,2,4,6,8,5,7 / data iss_rear_v3 / 1,3,2,4,6,8,5,7 / data iss_rear_u3 / 8,6,7,5,3,1,4,2 / data iss_rear_x3 / 1,3,2,4,6,8,5,7 / data iss_rear_u4 / 8,6,7,5,3,1,4,2 / data iss_rear_v4 / 8,6,7,5,3,1,4,2 / c if (debug_fpp) then c c c Load VDC tracking variables into local variables - write them out c for diagnostics if desired. c write(25,*)'ninterupts = ',nr_interrupts write(25,*)'FP_d coordinates' write(25,*)'vdc x in cm =',sptransport.l.particle.fp_d.x*100. write(25,*)'vdc y in cm =',sptransport.l.particle.fp_d.y*100. write(25,*)'vdc ph in rad =',sptransport.l.particle.fp_d.ph write(25,*)'vdc th in rad =',sptransport.l.particle.fp_d.th write(25,*)'FP_H coordinates' write(25,*)'vdc x in cm =',sptransport.l.particle.fp_h.x*100. write(25,*)'vdc y in cm =',sptransport.l.particle.fp_h.y*100. write(25,*)'vdc ph in rad =',sptransport.l.particle.fp_h.ph write(25,*)'vdc th in rad =',sptransport.l.particle.fp_h.th write(25,*)'FP_t coordinates' write(25,*)'vdc x in cm =',sptransport.l.particle.fp_t.x*100. write(25,*)'vdc y in cm =',sptransport.l.particle.fp_t.y*100. write(25,*)'vdc ph in rad =',sptransport.l.particle.fp_t.ph write(25,*)'vdc th in rad =',sptransport.l.particle.fp_t.th endif C c c IMPORTANT: changed to spec_e vdc track variables here for RCS, c using Left Spectrometer. djh May 02 c xinit = sptransport.l.particle.fp_t.x*100.0 yinit = sptransport.l.particle.fp_t.y*100.0 thinit = sptransport.l.particle.fp_t.ph thinit = datan(thinit) phiinit=datan(sptransport.l.particle.fp_t.th) c if (debug_fpp) then write(25,*) ' Calc vdc quantites' write(25,*) ' x = ',xinit,' y = ',yinit, > ' dydz = ',dtan(thinit),' dxdz = ',dtan(phiinit) endif c c Initialize Hit array c nevent = nevent + 1 do kk=1,5 do jj=1,35 do ii=1,24 hit(kk,jj,ii)=0.0 enddo enddo enddo c c Initialize Track array c do ii=1,40 do jj=1,25 track(ii,jj)=0.0 enddo enddo c c c Initialize nhit array c do ii=1,25 nhit(ii)=0 enddo c do ii=1,24 spdetector.l.fpp.plane(ii).nghits=0 enddo c c c Fill hit and nhit arrays with values from common blocks c ntothits=0 do jj=1,24 ngoodhit=0 nhit(jj)=spdetector.l.fpp.plane(jj).nhits nhinit=nhit(jj) prev_straw = -1 do ii=1,nhinit iwiregroup=spdetector.l.fpp.plane(jj).raw.wire(ii) letime=(spdetector.l.fpp.plane(jj).raw.ltdc(ii))/2.0 tetime=(spdetector.l.fpp.plane(jj).raw.ttdc(ii))/2.0 tdiff=letime-tetime spdetector.l.fpp.plane(jj).raw.letediff(ii) = tdiff imark=0 do kk=1,9 if(tdiff.gt.spdetector.l.fpp.plane(jj). > wiregroup(iwiregroup).cut(kk)) then imark=imark+1 endif enddo if(imark.gt.0.and.imark.le.8) then if(jj.ge.1.and.jj.le.12) diff=iss_front(imark) if(jj.ge.13.and.jj.le.14) diff=iss_rear_u3(imark) if(jj.ge.19.and.jj.le.21) diff=iss_rear_u4(imark) if(jj.ge.15.and.jj.le.16) diff=iss_rear_v3(imark) if(jj.ge.22.and.jj.le.24) diff=iss_rear_v4(imark) if(jj.ge.17.and.jj.le.18) diff=iss_rear_x3(imark) else if(imark.gt.8) then imark=0 if(jj.ge.1.and.jj.le.12) diff=iss_front(8) if(jj.ge.13.and.jj.le.14) diff=iss_rear_u3(8) if(jj.ge.19.and.jj.le.21) diff=iss_rear_u4(8) if(jj.ge.15.and.jj.le.16) diff=iss_rear_v3(8) if(jj.ge.22.and.jj.le.24) diff=iss_rear_v4(8) if(jj.ge.17.and.jj.le.18) diff=iss_rear_x3(8) else if(jj.ge.1.and.jj.le.12) diff=iss_front(1) if(jj.ge.13.and.jj.le.14) diff=iss_rear_u3(1) if(jj.ge.19.and.jj.le.21) diff=iss_rear_u4(1) if(jj.ge.15.and.jj.le.16) diff=iss_rear_v3(1) if(jj.ge.22.and.jj.le.24) diff=iss_rear_v4(1) if(jj.ge.17.and.jj.le.18) diff=iss_rear_x3(1) endif endif rstraw=8.*(iwiregroup-1)+diff istraw=rstraw spdetector.l.fpp.plane(jj).raw.straw(ii) = istraw c c Calculation of drift distance - changed to just be a call c of a subroutine. One can then put any calculation of the drift c distance that one likes in the fortran subroutine fpp_drift_dist c c call fpp_drift_dist(letime,toff(jj,iwiregroup), > jj,drdist) c if(letime.ge.300.and.letime.le.1400) then thits(jj,iwiregroup)=thits(jj,iwiregroup)+1.0 tsum(jj,iwiregroup)=tsum(jj,iwiregroup)+ > letime*1.0 endif if(tdiff.gt.15.0.and.tdiff.lt.110.) then do kk=1,8 if(icalib.ge.1.and.icalib.le.8) then if(jj.ge.1.and.jj.le.12) then if(icalib.eq.iss_front(kk))index=kk endif if(jj.ge.13.and.jj.le.14) then if(icalib.eq.iss_rear_u3(kk))index=kk endif if(jj.ge.19.and.jj.le.21) then if(icalib.eq.iss_rear_u4(kk))index=kk endif if(jj.ge.15.and.jj.le.16) then if(icalib.eq.iss_rear_v3(kk))index=kk endif if(jj.ge.22.and.jj.le.24) then if(icalib.eq.iss_rear_v4(kk))index=kk endif if(jj.ge.17.and.jj.le.18) then if(icalib.eq.iss_rear_x3(kk))index=kk endif else if(imark.ge.1.and.imark.le.8) then index=imark else index=1 endif endif enddo c c Try to decide if it was a reasonable hit or some cross talk - If c it is within 15 channels of the "standard" value, then include it. c If it is not then do not include it. 15 channels is probably ok c because of the ordering in the demuxing. We will never have cross c talk in the adjacent peak - it will always be about 2 peaks away. c if(tdiff.ge.spdetector.l.fpp.plane(jj). > wiregroup(iwiregroup).cut(index)-10.0.and. > tdiff.le.spdetector.l.fpp.plane(jj). > wiregroup(iwiregroup).cut(index)+20.0) then phits(jj,iwiregroup,index)= > phits(jj,iwiregroup,index)+1.0 psum(jj,iwiregroup,index)= > psum(jj,iwiregroup,index)+tdiff*1.0 endif endif itime=toff(jj,iwiregroup)-letime*1.0+40 if(itime.gt.0.and.itime.le.80) then thist(jj,iwiregroup,itime)= > thist(jj,iwiregroup,itime)+1 endif if((drdist.ge.0.0.and.drdist.le.0.522) > .or.pulser_fpp) then if(imark.eq.0) then nhit(jj)=nhit(jj)-1 else if(pulser_fpp) then ntothits=ntothits+1 ngoodhit=ngoodhit+1 spdetector.l.fpp.plane(jj).cor.ltdc(ngoodhit)= > toff(jj,iwiregroup)-letime*1.0 spdetector.l.fpp.plane(jj).cor. > letediff(ngoodhit)=tdiff spdetector.l.fpp.plane(jj).cor. > straw(ngoodhit)=istraw spdetector.l.fpp.plane(jj).cor. > drdist(ngoodhit)=drdist spdetector.l.fpp.plane(jj).cor. > wg(ngoodhit)=iwiregroup hit(1,ngoodhit,jj)=istraw*1.0 !straw number within plane hit(2,ngoodhit,jj)=drdist !drift distance in cm hit(5,ngoodhit,jj)=jj !write plane number else ntothits=ntothits+1 ngoodhit=ngoodhit+1 if ( ngoodhit .ge. 2 .and. > prev_straw .eq. istraw) then ngoodhit= ngoodhit - 1 nhit(jj)=nhit(jj)-1 endif ! if multi-hit in straw want last hit spdetector.l.fpp.plane(jj).cor.ltdc(ngoodhit)= > toff(jj,iwiregroup)-letime*1.0 spdetector.l.fpp.plane(jj).cor. > letediff(ngoodhit)=tdiff spdetector.l.fpp.plane(jj).cor. > straw(ngoodhit)=istraw prev_straw = istraw spdetector.l.fpp.plane(jj).cor. > drdist(ngoodhit)=drdist spdetector.l.fpp.plane(jj).cor. > wg(ngoodhit)=iwiregroup hit(1,ngoodhit,jj)=istraw*1.0 !straw number within plane hit(2,ngoodhit,jj)=drdist !drift distance in cm hit(5,ngoodhit,jj)=jj !plane number endif endif else nhit(jj)=nhit(jj)-1 endif enddo spdetector.l.fpp.plane(jj).nghits=ngoodhit enddo spdetector.l.fpp.nhits=ntothits c if (debug_fpp) then write(25,*) ' Raw hits ' write(25,*) 'plane #hit WG straw ttdc ltdc' do jj=1,24 do ii=1,spdetector.l.fpp.plane(jj).nhits write(25,'(3(3x,i3),2x,f5.1,2x,f8.2,2x,f8.2)') jj,ii, 1 spdetector.l.fpp.plane(jj).raw.wire(ii), 1 spdetector.l.fpp.plane(jj).raw.straw(ii), 1 spdetector.l.fpp.plane(jj).raw.ttdc(ii)/2.0, 1 spdetector.l.fpp.plane(jj).raw.ltdc(ii)/2.0 enddo enddo write(25,*) ' Good hits ' write(25,*) 'plane #hit WG straw letediff drdist' do jj=1,24 if (nhit(jj) .ne. spdetector.l.fpp.plane(jj).nghits) then write(25,*) 'nhits .ne. nghits' endif do ii=1,nhit(jj) write(25,'(3(3x,i3),2x,f5.1,2x,f7.1,2x,f8.2)') jj,ii, 1 spdetector.l.fpp.plane(jj).cor.wg(ii), 1 spdetector.l.fpp.plane(jj).cor.straw(ii), 1 spdetector.l.fpp.plane(jj).cor.letediff(ii), 1 spdetector.l.fpp.plane(jj).cor.drdist(ii) enddo enddo endif c c Now calculate the z and y (or x) position of the straw knowing what c plane it is and what straw it is within the plane. c EJB method if (.not. fpp_mkjtrack) then do jj=1,12 do ii=1,nhit(jj) hit(3,ii,jj)=uvx_zero(jj) + w_spacingf*(hit(1,ii,jj)-1) hit(4,ii,jj)=z_fpp(jj) enddo enddo do jj=13,24 do ii=1,nhit(jj) hit(3,ii,jj)=uvx_zero(jj) + w_spacingr*(hit(1,ii,jj)-1) hit(4,ii,jj)=z_fpp(jj) enddo enddo endif c if (fpp_mkjtrack) then c KW First find rough u,v positions for each wire plane using wire positions c only. This can be used to correct for chamber mis-alignments c do jj=1,24 do ii=1,nhit(jj) hit(3,ii,jj)=uvx_first(jj) + uvx_shift(jj) & + w_spacing(jj)*(hit(1,ii,jj)-1) hit(4,ii,jj)=zmid(jj)+zshift(jj) enddo enddo endif c c c Call tracking routine, sending it the nhit and hit arrays. c c Front U tracking - itrcode = 1 c Front V tracking - itrcode = 2 c Rear U tracking - itrcode = 3 c Rear V tracking - itrcode = 4 c Rear X tracking - itrcode = 5 c 100 format(1x,i3,1x,i3,4(1x,f10.5)) icode=0 if(pulser_fpp) then icode=1 do i=1,5 itrk(i)=0 enddo goto 1324 else if(ana_fpp_front) then itrcode = 1 itrk(itrcode)=0 call trk4(itrcode,nhit,hit,ntrack,track, & usave,zusave,dusave,uplsave) ntr(itrcode) = ntrack if (ntrack .ge. 1) itrk(itrcode)=1 if(ntrack.eq.0) then icode=1 spdetector.l.fpp.chisqfu = rbig else spdetector.l.fpp.chisqfu = track(12,1) endif itrcode = 2 itrk(itrcode)=0 call trk4(itrcode,nhit,hit,ntrack,track, & vsave,zvsave,dvsave,vplsave) ntr(itrcode) = ntrack if (ntrack .ge. 1) itrk(itrcode)=1 if(ntrack.eq.0) then icode=1 spdetector.l.fpp.chisqfv = rbig else spdetector.l.fpp.chisqfv = track(11,1) endif if (itrk(1) .eq. 1 .and. itrk(2) .eq. 1 & .and. fpp_mkjtrack) then if (debug_fpp) then write(25,*) ' Correct front u track ' endif ave1 = 0 ave2 = 0 zave1 = 0 zave2 = 0 npt1 = 0 npt2 = 0 do i=1,int(track(15,1)) npl = vplsave(i) if ( npl .le. 6) then ave1 = ave1 + vsave(i) zave1 = zave1 + zvsave(i) npt1 = npt1 + 1 else ave2 = ave2 + vsave(i) zave2 = zave2 + zvsave(i) npt2 = npt2 + 1 endif enddo ave1 = ave1/npt1 zave1 = zave1/npt1 ave2 = ave2/npt2 zave2 = zave2/npt2 if ( (zave2-zave1) .ne.0 )then ang_cr = (ave2-ave1)/(zave2-zave1) else ang_cr = 0 endif int_cr = ave1 - ang_cr*zave1 if (debug_fpp) then write(25,*) ' crude calc ' write(25,'(4f10.5)') & int_cr,ang_cr,track(3,1),track(4,1) write(25,*) ' vs us zs vs1 zs1 ', &' us2 zs2 us3 psi_tr ' endif do i=1,int(track(16,1)) npl = uplsave(i) c vs = track(3,1) + track(4,1)*zusave(i) vs = int_cr + ang_cr*zusave(i) us = usave(i) psi_tr = atan(tan(zurot(npl)) * cos(zvrot(npl))) zs = zshift(npl) vs1= cos(zvrot(npl))*vs + sin(zvrot(npl))*zs zs1= cos(zvrot(npl))*zs - sin(zvrot(npl))*vs us2= cos(psi_tr)*us - sin(psi_tr)*zs1 zs2= cos(psi_tr)*zs1 + sin(psi_tr)*us us3= cos(uvrot(npl))*us2 - sin(uvrot(npl))*vs1 if (debug_fpp) then write(25,'(9f9.5)') vs,us,zs,vs1,zs1,us2,zs2,us3,psi_tr endif zone(i) = zs2 + zmid(npl) xone(i) = us3 dxone(i) = dusave(i) uone(i) = 1 dxfone(i) = 0. enddo if (debug_fpp) then write(25,*) ' plane u du z ' do i=1,int(track(16,1)) write(25,*) uplsave(i),xone(i),dxone(i),zone(i) enddo endif tr_num = int(track(16,1)) call trackit(tr_num,zone,xone,dxone, > uone,dxfone,lrfit,uat0,uslope,uchi,chis) if (debug_fpp) write(25,'(i4,3f10.5)') & tr_num,uat0,uslope,uchi if (debug_fpp) then write(25,*) ' Correct front v track ' endif ave1 = 0 ave2 = 0 zave1 = 0 zave2 = 0 npt1 = 0 npt2 = 0 do i=1,int(track(16,1)) npl = uplsave(i) if ( npl .le. 6) then ave1 = ave1 + usave(i) zave1 = zave1 + zusave(i) npt1 = npt1 + 1 else ave2 = ave2 + usave(i) zave2 = zave2 + zusave(i) npt2 = npt2 + 1 endif enddo ave1 = ave1/npt1 zave1 = zave1/npt1 ave2 = ave2/npt2 zave2 = zave2/npt2 if ( (zave2-zave1) .ne.0 )then ang_cr = (ave2-ave1)/(zave2-zave1) else ang_cr = 0 endif int_cr = ave1 - ang_cr*zave1 if (debug_fpp) then write(25,*) ' crude calc ' write(25,'(4f10.5)') & int_cr,ang_cr,track(1,1),track(2,1) write(25,*) ' vs us zs vs1 zs1 ', &' us2 zs2 vs3 psi_tr ' endif do i=1,int(track(15,1)) npl = vplsave(i) c us = track(1,1) + track(2,1)*zvsave(i) us = int_cr + ang_cr*zvsave(i) vs = vsave(i) psi_tr = atan(tan(zurot(npl)) * cos(zvrot(npl))) zs = zshift(npl) vs1 = cos(zvrot(npl))*vs + sin(zvrot(npl))*zs zs1 = cos(zvrot(npl))*zs - sin(zvrot(npl))*vs us2 = cos(psi_tr)*us - sin(psi_tr)*zs1 zs2 = cos(psi_tr)*zs1 + sin(psi_tr)*us vs3 = sin(uvrot(npl))*us2 + cos(uvrot(npl))*vs1 if (debug_fpp) then write(25,'(9f9.5)') vs,us,zs,vs1,zs1,us2,zs2,vs3,psi_tr endif zone(i) = zs2 + zmid(npl) xone(i) = -sin(uvno(npl))*us2 + cos(uvno(npl))*vs3 dxone(i) = dvsave(i) uone(i) = 1 dxfone(i) = 0. enddo if (debug_fpp) then write(25,*) ' plane v dv z ' do i=1,int(track(15,1)) write(25,*) vplsave(i),xone(i),dxone(i),zone(i) enddo endif tr_num = int(track(15,1)) call trackit(tr_num,zone,xone,dxone, > uone,dxfone,lrfit,vat0,vslope,vchi,chis) if (debug_fpp) write(25,'(i4,3f10.5)') & tr_num,vat0,vslope,vchi track(1,1) = uat0 track(2,1) = uslope track(12,1) = uchi spdetector.l.fpp.chisqfu = track(12,1) track(3,1) = vat0 track(4,1) = vslope track(11,1) = vchi spdetector.l.fpp.chisqfv = track(11,1) endif else itrk(1)=0 itrk(2)=0 endif cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc if(ana_fpp_rear) then itrcode = 3 itrk(itrcode)=0 call trk4(itrcode,nhit,hit,ntrack,track, & usave,zusave,dusave,uplsave) ntr(itrcode) = ntrack if (ntrack .ge. 1) itrk(itrcode)=1 if(ntrack.eq.0) then icode=1 spdetector.l.fpp.chisqru = rbig else spdetector.l.fpp.chisqru = track(14,1) endif itrcode = 4 itrk(itrcode)=0 call trk4(itrcode,nhit,hit,ntrack,track, & vsave,zvsave,dvsave,vplsave) ntr(itrcode) = ntrack if (ntrack .ge. 1) itrk(itrcode)=1 if(ntrack.eq.0) then icode=1 spdetector.l.fpp.chisqrv = rbig else spdetector.l.fpp.chisqrv = track(13,1) endif if (itrk(3) .eq. 1 .and. itrk(4) .eq. 1 & .and. fpp_mkjtrack ) then if (debug_fpp) then write(25,*) ' Correct rear u track ' endif ave1 = 0 ave2 = 0 zave1 = 0 zave2 = 0 npt1 = 0 npt2 = 0 do i=1,int(track(17,1)) npl = vplsave(i) if ( npl .le. 16) then ave1 = ave1 + vsave(i) zave1 = zave1 + zvsave(i) npt1 = npt1 + 1 else ave2 = ave2 + vsave(i) zave2 = zave2 + zvsave(i) npt2 = npt2 + 1 endif enddo ave1 = ave1/npt1 zave1 = zave1/npt1 ave2 = ave2/npt2 zave2 = zave2/npt2 if ( (zave2-zave1) .ne.0 )then ang_cr = (ave2-ave1)/(zave2-zave1) else ang_cr = 0 endif int_cr = ave1 - ang_cr*zave1 if (debug_fpp) then write(25,*) ' crude calc ' write(25,'(4f10.5)') & ave1,zave1,ave2,zave2 write(25,'(4f10.5)') & int_cr,ang_cr,track(7,1),track(8,1) write(25,*) ' vs us zs vs1 zs1 ', &' us2 zs2 us3 psi_tr ' endif do i=1,int(track(18,1)) npl = uplsave(i) if ( npl .eq. 13 .or. npl .eq. 14) then npl_skip = 2 else npl_skip = 3 endif c vs = track(7,1) + track(8,1)*zusave(i) vs = int_cr + ang_cr*zusave(i) us = usave(i) psi_tr = atan(tan(zurot(npl)) * cos(zvrot(npl))) zs = zshift(npl) vs1= cos(zvrot(npl))*vs + sin(zvrot(npl))*zs zs1= cos(zvrot(npl))*zs - sin(zvrot(npl))*vs us2= cos(psi_tr)*us - sin(psi_tr)*zs1 zs2= cos(psi_tr)*zs1 + sin(psi_tr)*us us3= cos(uvrot(npl))*us2 - sin(uvrot(npl))*vs1 if (debug_fpp) then write(25,'(8f10.5)') vs,us,vs1,us1,us2,zs2,us3,zmid(npl) endif zone(i) = zs2 + zmid(npl) xone(i) = us3 dxone(i) = dusave(i) uone(i) = 1 dxfone(i) = 0. enddo if (debug_fpp) then write(25,*) ' plane u du z ' do i=1,int(track(18,1)) write(25,*) uplsave(i),xone(i),dxone(i),zone(i) enddo endif tr_num = int(track(18,1)) call trackit(tr_num,zone,xone,dxone, > uone,dxfone,lrfit,uat0,uslope,uchi,chis) if (debug_fpp) write(25,'(i4,3f10.5)') & tr_num,uat0,uslope,uchi if (debug_fpp) then write(25,*) ' Correct rear v track ' write(25,*) ' plane v dv z ' endif ave1 = 0 ave2 = 0 zave1 = 0 zave2 = 0 npt1 = 0 npt2 = 0 do i=1,int(track(18,1)) npl = uplsave(i) if ( npl .le. 16) then ave1 = ave1 + usave(i) zave1 = zave1 + zusave(i) npt1 = npt1 + 1 else ave2 = ave2 + usave(i) zave2 = zave2 + zusave(i) npt2 = npt2 + 1 endif enddo ave1 = ave1/npt1 zave1 = zave1/npt1 ave2 = ave2/npt2 zave2 = zave2/npt2 if ( (zave2-zave1) .ne.0 )then ang_cr = (ave2-ave1)/(zave2-zave1) else ang_cr = 0 endif int_cr = ave1 - ang_cr*zave1 if (debug_fpp) then write(25,*) ' crude calc ' write(25,'(4f10.5)') & ave1,zave1,ave2,zave2 write(25,'(4f10.5)') & int_cr,ang_cr,track(5,1),track(6,1) write(25,*) ' vs us zs vs1 zs1 ', &' us2 zs2 vs3 psi_tr ' endif do i=1,int(track(17,1)) npl = vplsave(i) if ( npl .eq. 15 .or. npl .eq. 16) then npl_skip = 2 else npl_skip = 3 endif c us = track(5,1) + track(6,1)*zvsave(i) us = int_cr + ang_cr*zvsave(i) vs = vsave(i) psi_tr = atan(tan(zurot(npl)) * cos(zvrot(npl))) zs = zshift(npl) vs1 = cos(zvrot(npl))*vs + sin(zvrot(npl))*zs zs1 = cos(zvrot(npl))*zs - sin(zvrot(npl))*vs us2 = cos(psi_tr)*us - sin(psi_tr)*zs1 zs2 = cos(psi_tr)*zs1 + sin(psi_tr)*us vs3 = sin(uvrot(npl))*us2 + cos(uvrot(npl))*vs1 if (debug_fpp) then write(25,'(9f9.5)') vs,us,zs,vs1,zs1,us2,zs2,vs3,zmid(npl) endif zone(i) = zs2 + zmid(npl) xone(i) = -sin(uvno(npl))*us2 + cos(uvno(npl))*vs3 dxone(i) = dvsave(i) uone(i) = 1 dxfone(i) = 0. enddo if (debug_fpp) then do i=1,int(track(17,1)) write(25,*) vplsave(i),xone(i),dxone(i),zone(i) enddo endif tr_num = int(track(17,1)) call trackit(tr_num,zone,xone,dxone, > uone,dxfone,lrfit,vat0,vslope,vchi,chis) if (debug_fpp) write(25,'(i4,3f10.5)') & tr_num,vat0,vslope,vchi track(7,1) = vat0 track(8,1) = vslope track(13,1) = vchi spdetector.l.fpp.chisqrv = track(13,1) track(5,1) = uat0 track(6,1) = uslope track(14,1) = uchi spdetector.l.fpp.chisqru = track(14,1) endif else itrk(3)=0 itrk(4)=0 itrk(5)=0 endif endif 1324 continue if (debug_fpp) then write(25,*) + '# Tracks for Front U Front V Rear U Rear V' write(25,*) + ' ------- ------ ------ ------' write(25,'(16x,i2,5x,i2,5x,i2,5x,i2)') (itrk (i),i=1,4) write(25,'(a,10x,4(g12.5))') 'chisq', + track(12,1),track(11,1),track(14,1),track(13,1) endif spdetector.l.fpp.ntrackfu = itrk(1) spdetector.l.fpp.ntrackfv = itrk(2) spdetector.l.fpp.ntrackru = itrk(3) spdetector.l.fpp.ntrackrv = itrk(4) ntrack_front = itrk(1) + itrk(2) ntrack_rear = itrk(3) + itrk(4) ntrack_tot = ntrack_front + ntrack_rear c spdetector.l.fpp.th_sc_CH2=rbig spdetector.l.fpp.ph_sc_CH2=rbig spdetector.l.fpp.th_az_CH2=rbig spdetector.l.fpp.ph_az_CH2=rbig spdetector.l.fpp.zclose_CH2=rbig spdetector.l.fpp.sclose_CH2=rbig spdetector.l.fpp.conetest_CH2=1000 spdetector.l.fpp.x_fa=rbig spdetector.l.fpp.y_fa=rbig spdetector.l.fpp.x_fb=rbig spdetector.l.fpp.y_fb=rbig spdetector.l.fpp.th_f=rbig spdetector.l.fpp.ph_f=rbig spdetector.l.fpp.x_diff_f=rbig spdetector.l.fpp.y_diff_f=rbig spdetector.l.fpp.th_diff_f=rbig spdetector.l.fpp.ph_diff_f=rbig c spdetector.l.fpp.conetest_car=1000 spdetector.l.fpp.x_fa=rbig spdetector.l.fpp.y_fa=rbig spdetector.l.fpp.x_fb=rbig spdetector.l.fpp.y_fb=rbig spdetector.l.fpp.x_r=rbig spdetector.l.fpp.y_r=rbig spdetector.l.fpp.x_x=rbig c spdetector.l.fpp.th_f=rbig spdetector.l.fpp.ph_f=rbig spdetector.l.fpp.th_r=rbig spdetector.l.fpp.ph_r=rbig c spdetector.l.fpp.x_diff_f=rbig spdetector.l.fpp.y_diff_f=rbig spdetector.l.fpp.x_diff_r=rbig spdetector.l.fpp.y_diff_r=rbig spdetector.l.fpp.x_diff_x=rbig spdetector.l.fpp.th_diff_f=rbig spdetector.l.fpp.ph_diff_f=rbig spdetector.l.fpp.th_diff_r=rbig spdetector.l.fpp.ph_diff_r=rbig c spdetector.l.fpp.th_sc_car=rbig spdetector.l.fpp.ph_sc_car=rbig c do kk=1,24 spdetector.l.fpp.plane(kk).itrgood=10000 spdetector.l.fpp.plane(kk).strwgood=10000 spdetector.l.fpp.plane(kk).strwbad=10000 do ii=1,nhit(kk) spdetector.l.fpp.plane(kk).cor.residual(ii)= & rbig enddo enddo c spdetector.l.fpp.th_az_car=rbig spdetector.l.fpp.ph_az_car=rbig spdetector.l.fpp.zclose_car=rbig spdetector.l.fpp.sclose_car=rbig spdetector.l.fpp.x_r=rbig spdetector.l.fpp.y_r=rbig spdetector.l.fpp.th_r=rbig spdetector.l.fpp.ph_r=rbig c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c FPP Reconstruction c c Modified for RCS FPP Setup (CH2 and Carbon Analyzers) c c djh, Feb 02 c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Calculate Track paramters for front and rear then align fg=180.00000/3.141592653589793 rbig=1.0d12 xatz0f=(track(1,1)-track(3,1))/1.414213562 yatz0f=(track(1,1)+track(3,1))/1.414213562 xatz0r=(track(5,1)-track(7,1))/1.414213562 yatz0r=(track(5,1)+track(7,1))/1.414213562 xatz0x=track(9,1) tthetaf=(track(2,1)+track(4,1))/1.414213562 tphif=(track(2,1)-track(4,1))/1.414213562 tthetar=(track(6,1)+track(8,1))/1.414213562 tphir=(track(6,1)-track(8,1))/1.414213562 thetafd=datan(tthetaf)*fg phifd=datan(tphif)*fg tthetaa=dtan(thinit) tphia=dtan(phiinit) thinitd=thinit*fg phiinitd=phiinit*fg if (.not. fpp_mkjtrack) then call fpp_align_corr (xatz0f,yatz0f,tthetaf,tphif, & xatz0r,yatz0r,tthetar,tphir) endif cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c CH2 Polarimeter - VDC/Front c c Variable_CH2 c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc if( ana_fpp_front .and. ntrack_front .eq. 2 ) then xdiff_front=xinit-xatz0f ydiff_front=yinit-yatz0f spdetector.l.fpp.x_fa=xatz0f spdetector.l.fpp.y_fa=yatz0f theta_difff=tthetaf-tthetaa phi_difff=tphif-tphia spdetector.l.fpp.th_f=fg*datan(tthetaf) spdetector.l.fpp.ph_f=fg*datan(tphif) spdetector.l.fpp.x_diff_f=xdiff_front spdetector.l.fpp.y_diff_f=ydiff_front spdetector.l.fpp.th_diff_f=fg*datan(theta_difff) spdetector.l.fpp.ph_diff_f=fg*datan(phi_difff) if( xinit .ne. 1.0d17 ) then if( ialign .eq. 10 ) then write(35,*) xatz0f,yatz0f,tthetaf,tphif,xinit,yinit,tthetaa, & tphia endif c Cartesian Scattering Angles psia=datan(tphia*dcos(datan(tthetaa))) psif=datan(tphif*dcos(datan(tthetaf))) xhatf=dsin(psif) yhatf=dcos(psif)*dsin(datan(tthetaf)) zhatf=dcos(psif)*dcos(datan(tthetaf)) xhat1f=xhatf*dcos(psia)- & yhatf*dsin(datan(tthetaa))*dsin(psia)- & zhatf*dcos(datan(tthetaa))*dsin(psia) yhat1f= & yhatf*dcos(datan(tthetaa))- & zhatf*dsin(datan(tthetaa)) zhat1f=xhatf*dsin(psia)+ & yhatf*dsin(datan(tthetaa))*dcos(psia)+ & zhatf*dcos(datan(tthetaa))*dcos(psia) psisc_CH2=dasin(xhat1f) if(dcos(psisc_CH2).ne.0.) then thetasc_CH2=dasin(yhat1f/dcos(psisc_CH2)) else write(*,*)'Zero problem with psisc_CH2' endif if(dcos(thetasc_CH2).ne.0.) then phisc_CH2=datan(dtan(psisc_CH2)/dcos(thetasc_CH2)) else write(*,*)'Zero problem with thetasc_CH2' endif theta_sc_CH2=thetasc_CH2 phi_sc_CH2=phisc_CH2 spdetector.l.fpp.th_sc_CH2=fg*theta_sc_CH2 spdetector.l.fpp.ph_sc_CH2=fg*phi_sc_CH2 c Zclose and Sclose a0_CH2 = (xatz0f-xinit)**2 + (yatz0f-yinit)**2 b1_CH2 = tphia**2 +tthetaa**2 + 1.0 b2_CH2 = tphif**2 +tthetaf**2 + 1.0 c1_CH2 = (xatz0f-xinit)*tphia + (yatz0f-yinit)*tthetaa c2_CH2 = (xatz0f-xinit)*tphif + (yatz0f-yinit)*tthetaf c3_CH2 = tphia*tphif + tthetaa*tthetaf + 1.0 lambda1_CH2 = (c1_CH2*b2_CH2 - c2_CH2*c3_CH2) / & (b1_CH2*b2_CH2 - c3_CH2*c3_CH2) lambda2_CH2 = (c1_CH2*c3_CH2 - c2_CH2*b1_CH2) / & (b1_CH2*b2_CH2 - c3_CH2*c3_CH2) zclose_CH2 = 0.5 * (lambda1_CH2 + lambda2_CH2) sclose_CH2 = sqrt(a0_CH2+lambda1_CH2**2*b1_CH2 +lambda2_CH2 & **2*b2_CH2+2.0*(-lambda1_CH2*c1_CH2+lambda2_CH2*c2_CH2 & -lambda1_CH2*lambda2_CH2*c3_CH2)) spdetector.l.fpp.zclose_CH2=zclose_CH2 spdetector.l.fpp.sclose_CH2=sclose_CH2 c Azimuthal Angles psi_sc_CH2=datan(dtan(phi_sc_CH2)*dcos(theta_sc_CH2)) x_sc_CH2=dsin(psi_sc_CH2) y_sc_CH2=dcos(psi_sc_CH2)*dsin(theta_sc_CH2) z_sc_CH2=dcos(psi_sc_CH2)*dcos(theta_sc_CH2) r_sc_CH2=sqrt(x_sc_CH2**2+y_sc_CH2**2) if(z_sc_CH2.ne.0..and.y_sc_CH2.ne.0.) then theta_az_CH2=datan(r_sc_CH2/z_sc_CH2) phi_az_CH2=datan(x_sc_CH2/y_sc_CH2) else write(*,*)'Zero problem with zsc or ysc' endif theta_az_CH2=theta_az_CH2*180./pi phi_az_CH2=phi_az_CH2*180./pi if(y_sc_CH2.lt.0.0) phi_az_CH2=phi_az_CH2+180.0 if(y_sc_CH2.ge.0.0.and.x_sc_CH2.lt.0.0) phi_az_CH2= & phi_az_CH2+360.0 spdetector.l.fpp.th_az_CH2=theta_az_CH2 spdetector.l.fpp.ph_az_CH2=phi_az_CH2 c Conetest for Front Chambers xf1=tphia*z_fpp(3)+xinit yf1=tthetaa*z_fpp(3)+yinit xr1=tphif*z_fpp(3)+xatz0f yr1=tthetaf*z_fpp(3)+yatz0f xf2=tphia*z_fpp(9)+xinit yf2=tthetaa*z_fpp(9)+yinit xr2=tphif*z_fpp(9)+xatz0f yr2=tthetaf*z_fpp(9)+yatz0f dcone1=sqrt((xf1-xr1)**2+(yf1-yr1)**2) dcone2=sqrt((xf2-xr2)**2+(yf2-yr2)**2) if ( (xf1+dcone1).lt.104.5 .and. (xf1+dcone1).gt.-104.5 $ .and.(xf2+dcone2).lt.104.5 .and. (xf2+dcone2).gt.-104.5 $ .and.(xf1-dcone1).lt.104.5 .and. (xf1-dcone1).gt.-104.5 $ .and.(xf2-dcone2).lt.104.5 .and. (xf2-dcone2).gt.-104.5 $ .and.(yf1+dcone1).lt.30. .and. (yf1+dcone1).gt.-30. $ .and.(yf1-dcone1).lt.30. .and. (yf1-dcone1).gt.-30. $ .and.(yf2+dcone2).lt.30. .and. (yf2+dcone2).gt.-30. $ .and.(yf2-dcone2).lt.30. .and. (yf2-dcone2).gt.-30. $ .and. xf1.lt.104.5 .and. xf1.gt.-104.5 .and. $ xf2.lt.104.5 .and. xf2.gt.-104.5 .and. $ yf1.lt.30. .and. yf1.gt.-30. .and. $ yf2.lt.30. .and. yf2.gt.-30. ) then spdetector.l.fpp.conetest_CH2 = 1 else spdetector.l.fpp.conetest_CH2 = -1 endif endif endif cccccccccccccccccccccccccccccccccccccccccccccccccccccc c Carbon Polarimeter - Front/Rear c c Variable_car c cccccccccccccccccccccccccccccccccccccccccccccccccccccc c if ( ntrack_rear .eq. 2 .and. ana_fpp_rear ) then xrold=xatz0r yrold=yatz0r trold=tthetar prold=tphir xatz0r=tphir*350.0+xatz0r yatz0r=tthetar*350.0+yatz0r spdetector.l.fpp.x_r=xatz0r spdetector.l.fpp.y_r=yatz0r spdetector.l.fpp.th_r=fg*datan(tthetar) spdetector.l.fpp.ph_r=fg*datan(tphir) if ( ntrack_front .eq. 2 .and. ana_fpp_front ) then xfaold=xatz0f yfaold=yatz0f xfbold=tphif*350.0+xatz0f yfbold=tthetaf*350.0+yatz0f tfold=tthetaf pfold=tphif xdiff_front=xinit-xatz0f ydiff_front=yinit-yatz0f spdetector.l.fpp.x_fa=xatz0f spdetector.l.fpp.y_fa=yatz0f theta_difff=tthetaf-tthetaa phi_difff=tphif-tphia theta_diffr=tthetar-tthetaf phi_diffr=tphir-tphif spdetector.l.fpp.th_f=fg*datan(tthetaf) spdetector.l.fpp.ph_f=fg*datan(tphif) xatz0f=tphif*350.0+xatz0f yatz0f=tthetaf*350.0+yatz0f xdiff_rear=xatz0f-xatz0r ydiff_rear=yatz0f-yatz0r xdiff_x=xatz0f-xatz0x spdetector.l.fpp.x_fb=xatz0f spdetector.l.fpp.y_fb=yatz0f spdetector.l.fpp.x_x=xatz0x spdetector.l.fpp.x_diff_f=xdiff_front spdetector.l.fpp.y_diff_f=ydiff_front spdetector.l.fpp.x_diff_r=xdiff_rear spdetector.l.fpp.y_diff_r=ydiff_rear spdetector.l.fpp.x_diff_x=xdiff_x spdetector.l.fpp.th_diff_f=fg*datan(theta_difff) spdetector.l.fpp.ph_diff_f=fg*datan(phi_difff) spdetector.l.fpp.th_diff_r=fg*datan(theta_diffr) spdetector.l.fpp.ph_diff_r=fg*datan(phi_diffr) c Cartesian Scattering Angles psif=datan(tphif*dcos(datan(tthetaf))) psir=datan(tphir*dcos(datan(tthetar))) xhatr=dsin(psir) yhatr=dcos(psir)*dsin(datan(tthetar)) zhatr=dcos(psir)*dcos(datan(tthetar)) xhat1r=xhatr*dcos(psif)- & yhatr*dsin(datan(tthetaf))*dsin(psif)- & zhatr*dcos(datan(tthetaf))*dsin(psif) yhat1r= & yhatr*dcos(datan(tthetaf))- & zhatr*dsin(datan(tthetaf)) zhat1r=xhatr*dsin(psif)+ & yhatr*dsin(datan(tthetaf))*dcos(psif)+ & zhatr*dcos(datan(tthetaf))*dcos(psif) psisc_car=dasin(xhat1r) if(dcos(psisc_car).ne.0.) then thetasc_car=dasin(yhat1r/dcos(psisc_car)) else write(*,*)'Zero problem with psisc_car' endif if(dcos(thetasc_car).ne.0.) then phisc_car=datan(dtan(psisc_car)/dcos(thetasc_car)) else write(*,*)'Zero problem with thetasc' endif theta_sc_car=thetasc_car phi_sc_car=phisc_car spdetector.l.fpp.th_sc_car=fg*theta_sc_car spdetector.l.fpp.ph_sc_car=fg*phi_sc_car xatz0f=xatz0f-350.0*tphif yatz0f=yatz0f-350.0*tthetaf xatz0r=xatz0r-350.0*tphir yatz0r=yatz0r-350.0*tthetar c ZClose and SClose a0_car = (xatz0r-xatz0f)**2 + (yatz0r-yatz0f)**2 b1_car = tphif**2 + tthetaf**2 + 1.0 b2_car = tphir**2 + tthetar**2 + 1.0 c1_car = (xatz0r-xatz0f)*tphif + (yatz0r-yatz0f)*tthetaf c2_car = (xatz0r-xatz0f)*tphir + (yatz0r-yatz0f)*tthetar c3_car = tphif*tphir + tthetaf*tthetar + 1.0 lambda1_car = (c1_car*b2_car - c2_car*c3_car) / & (b1_car*b2_car - c3_car*c3_car) lambda2_car = (c1_car*c3_car - c2_car*b1_car) / & (b1_car*b2_car - c3_car*c3_car) zclose_car = 0.5 * (lambda1_car + lambda2_car) sclose_car = sqrt(a0_car + lambda1_car**2 * b1_car + & lambda2_car**2 * b2_car + 2.0*(-lambda1_car*c1_car & + lambda2_car*c2_car - lambda1_car*lambda2_car*c3_car)) spdetector.l.fpp.zclose_car=zclose_car spdetector.l.fpp.sclose_car=sclose_car c Azimuthal Angles psi_sc_car=datan(dtan(phi_sc_car)*dcos(theta_sc_car)) x_sc_car=dsin(psi_sc_car) y_sc_car=dcos(psi_sc_car)*dsin(theta_sc_car) z_sc_car=dcos(psi_sc_car)*dcos(theta_sc_car) r_sc_car=sqrt(x_sc_car**2+y_sc_car**2) if(z_sc_car.ne.0..and.y_sc_car.ne.0.) then theta_az_car=datan(r_sc_car/z_sc_car) phi_az_car=datan(x_sc_car/y_sc_car) else write(*,*)'Zero problem with zsc_car or ysc_car' endif theta_az_car=theta_az_car*180./pi phi_az_car=phi_az_car*180./pi if(y_sc_car.lt.0.0) phi_az_car=phi_az_car+180.0 if(y_sc_car.ge.0.0.and.x_sc_car.lt.0.0) phi_az_car= & phi_az_car+360.0 spdetector.l.fpp.th_az_car=theta_az_car spdetector.l.fpp.ph_az_car=phi_az_car c Conetest for Rear Chambers xf3=tphif*z_fpp(15)+xatz0f yf3=tthetaf*z_fpp(15)+yatz0f xr3=tphir*z_fpp(15)+xatz0r yr3=tthetar*z_fpp(15)+yatz0r xf4=tphif*z_fpp(21)+xatz0f yf4=tthetaf*z_fpp(21)+yatz0f xr4=tphir*z_fpp(21)+xatz0r yr4=tthetar*z_fpp(21)+yatz0r if (fpp_mkjtrack) then xf3=tphif*(zmid(15)+zshift(15))+xatz0f yf3=tthetaf*(zmid(15)+zshift(15))+yatz0f xr3=tphir*(zmid(15)+zshift(15))+xatz0r yr3=tthetar*(zmid(15)+zshift(15))+yatz0r xf4=tphif*(zmid(21)+zshift(21))+xatz0f yf4=tthetaf*(zmid(21)+zshift(21))+yatz0f xr4=tphir*(zmid(21)+zshift(21))+xatz0r yr4=tthetar*(zmid(21)+zshift(21))+yatz0r endif dcone3=sqrt((xf3-xr3)**2+(yf3-yr3)**2) dcone4=sqrt((xf4-xr4)**2+(yf4-yr4)**2) if ( (xf3+dcone3).lt.134. .and. (xf3+dcone3).gt.-134. .and. $ (xf4+dcone4).lt.149. .and. (xf4+dcone4).gt.-149. $ .and.(xf3-dcone3).lt.134. .and. (xf3-dcone3).gt.-134. $ .and.(xf4-dcone4).lt.149. .and. (xf4-dcone4).gt.-149. $ .and.(yf3+dcone3).lt.62. .and. (yf3+dcone3).gt.-62. $ .and.(yf3-dcone3).lt.62. .and. (yf3-dcone3).gt.-62. $ .and.(yf4+dcone4).lt.72. .and. (yf4+dcone4).gt.-72. $ .and.(yf4-dcone4).lt.72. .and. (yf4-dcone4).gt.-72. .and. $ xf3.lt.134. .and. xf3.gt.-134. .and. $ xf4.lt.149. .and. xf4.gt.-149. .and. $ yf3.lt.62. .and. yf3.gt.-62. .and. $ yf4.lt.72. .and. yf4.gt.-72 ) then spdetector.l.fpp.conetest_car = 1 else spdetector.l.fpp.conetest_car = -1 endif cccccccccccccccccccccccccccccccccccccccccccc c end of reconstruction c cccccccccccccccccccccccccccccccccccccccccccc c Residual and Efficiency Calculations !!!!!!!!!!!! c c Note that all calculations are done with the raw angle and c position values before alignment corrections. c if (.not. fpp_mkjtrack) then do kk=1,24 if(kk.eq.1.or.kk.eq.2.or.kk.eq.3) igroup=1 if(kk.eq.4.or.kk.eq.5.or.kk.eq.6) igroup=2 if(kk.eq.7.or.kk.eq.8.or.kk.eq.9) igroup=1 if(kk.eq.10.or.kk.eq.11.or.kk.eq.12) igroup=2 if(kk.eq.13.or.kk.eq.14) igroup=2 if(kk.eq.15.or.kk.eq.16) igroup=1 if(kk.eq.17.or.kk.eq.18) igroup=3 if(kk.eq.19.or.kk.eq.20.or.kk.eq.21) igroup=2 if(kk.eq.22.or.kk.eq.23.or.kk.eq.24) igroup=1 z(kk)=z_fpp(kk) if(kk.le.12) then xtr(kk)=xfaold+z(kk)*pfold ytr(kk)=yfaold+z(kk)*tfold strw_rad=w_spacingf/2.0 if(igroup.eq.1) then uvxtr(kk)=(ytr(kk)-xtr(kk))/1.414213562 else uvxtr(kk)=(ytr(kk)+xtr(kk))/1.414213562 endif else xtr(kk)=xrold+(z(kk))*prold ytr(kk)=yrold+(z(kk))*trold strw_rad=w_spacingr/2.0 if(igroup.eq.1) then uvxtr(kk)=(ytr(kk)-xtr(kk))/1.414213562 else if(igroup.eq.2) then uvxtr(kk)=(ytr(kk)+xtr(kk))/1.414213562 else uvxtr(kk)=xtr(kk) endif endif straw_pred(kk)=0 nstraws=spdetector.l.fpp.geom.plane(kk).nwire uvx0=spdetector.l.fpp.geom.plane(kk).uvx0*100. do ii=1,nstraws uvx_low=uvx0+2.*(strw_rad)*(ii-1.)-0.522 uvx_hi=uvx0+2.*(strw_rad)*(ii-1.)+0.522 if(uvxtr(kk).ge.uvx_low.and.uvxtr(kk).le.uvx_hi) then straw_pred(kk)=ii endif enddo itrgood(kk)=0 do ii=1,nhit(kk) str=hit(1,ii,kk) if(str.eq.straw_pred(kk)) itrgood(kk)=1 enddo if(straw_pred(kk).eq.0) itrgood(kk)=2 spdetector.l.fpp.plane(kk).itrgood=itrgood(kk) spdetector.l.fpp.plane(kk).strwgood=straw_pred(kk) spdetector.l.fpp.plane(kk).strwbad=straw_pred(kk) do ii=1,nhit(kk) ddist=hit(2,ii,kk) drdist=ddist str=hit(1,ii,kk) zplane=hit(4,ii,kk) uvx=hit(3,ii,kk) if(igroup.eq.1) then fac=1.0 if(kk.le.12) then a=pfold b=sqrt(2.)*uvx+xfaold c=tfold d=yfaold else a=prold b=sqrt(2.)*uvx+xrold c=trold d=yrold endif ee=zplane endif if(igroup.eq.2) then fac=1.0 if(kk.le.12) then a=-1.0*pfold b=sqrt(2.)*uvx-xfaold c=tfold d=yfaold else a=-1.0*prold b=sqrt(2.)*uvx-xrold c=trold d=yrold endif ee=zplane endif if(igroup.eq.3) then fac=0.0 a=prold b=-1.0*uvx+xrold c=trold d=yrold ee=zplane endif c c if(kk.eq.27) then write(*,*)'****** Initial Residual Stats *******' write(*,*)ddist,str,zplane,uvx write(*,*)tphia,xinit,tthetaa,yinit write(*,*)xfaold,yfaold,xrold,yrold write(*,*)a,b,c,d,ee write(*,*)'*************************************' endif ff=-1.0*fac*b-1.0*d gg=(1.0+a**2+c**2)/2. hh=a*b+c*d-ee iii=-1.0*fac*a-1.0*c jjj=(b**2+d**2+ee**2)/2. c c 1st iteration c ztr=zplane yw=-1.0*(ff+iii*ztr)/(1.0+fac**2) if(kk.eq.27) write(*,*)'1st pass - ywire =',yw ztr=-0.5*(hh+iii*yw)/gg if(kk.eq.27) write(*,*)'1st pass - ztr =',ztr c c 2nd iteration c yw=-1.0*(ff+iii*ztr)/(1.0+fac**2) if(kk.eq.27) write(*,*)'2nd pass - ywire =',yw ztr=-0.5*(hh+iii*yw)/gg if(kk.eq.27) write(*,*)'2nd pass - ztr =',ztr c c Now calculate distance of closest approach between wire hit c and track. c t1=0.5*(1.0+fac**2)*yw**2 t2=ff*yw t3=gg*ztr**2 t4=hh*ztr t5=iii*yw*ztr t6=jjj dd2by2=t1+t2+t3+t4+t5+t6 dd2=dd2by2*2. if(kk.eq.27) write(*,*)'Distance squared =',dd2 if(dd2.ge.0.) then dd=sqrt(dd2) resres=ddist-dd else dd=0. resres=0.49 endif residual=resres if(kk.eq.27) write(*,*) & 'P D1 D2 Res=',kk,ddist,dd,resres spdetector.l.fpp.plane(kk).cor.residual(ii)= & residual enddo enddo endif c c c c write(*,*)theta_diffr,phi_diffr,theta_sc,phi_sc c write(*,*)zclose,sclose,residual c write(*,*)theta_sc,phi_sc c if(evdsp_fpp .and. xinit .ne. 1e+17 c & .and. ( spdetector.l.fpp.chisqfu .gt. 3. c & .or. spdetector.l.fpp.chisqfv .gt. 3. c & .or. spdetector.l.fpp.chisqrv .gt. 3. c & .or. spdetector.l.fpp.chisqru .gt. 3. ) & ) then write(30,*) nr_interrupts write(30,*) theta_az_car,phi_az_car,zclose_car,sclose_car write(30,*) (yinit-xinit)/1.4142,(yinit+xinit)/1.4142, & (thinit-phiinit)/1.4142,(phiinit+thinit)/1.4142 write(30,*) (ntr(i),i=1,4) write(30,*) (track(1,i),track(2,i),track(12,i),i=1,ntr(1)) write(30,*) (track(3,i),track(4,i),track(11,i),i=1,ntr(2)) write(30,*) (track(5,i),track(6,i),track(14,i),i=1,ntr(3)) write(30,*) (track(7,i),track(8,i),track(13,i),i=1,ntr(4)) write(30,*) (nhit(jj),jj=1,24) do jj=1,24 write(30,*) (hit(2,i,jj),hit(3,i,jj) & ,hit(4,i,jj),i=1,nhit(jj)) enddo endif if (debug_fpp) then if ( zclose_car .gt. 435 + .and. zclose_car .lt. 441 .and. theta_az_car .gt. 4.5 + .and. theta_az_car .lt. 4.9 ) then write(27,*) ' Raw hits ' write(27,*) 'plane #hit WG straw ttdc ltdc' do jj=1,24 do ii=1,spdetector.l.fpp.plane(jj).nhits write(27,'(3(3x,i3),2x,f5.1,2x,f8.2,2x,f8.2)') jj,ii, 1 spdetector.l.fpp.plane(jj).raw.wire(ii), 1 spdetector.l.fpp.plane(jj).raw.straw(ii), 1 spdetector.l.fpp.plane(jj).raw.ttdc(ii)/2.0, 1 spdetector.l.fpp.plane(jj).raw.ltdc(ii)/2.0 enddo enddo write(27,*) ' Good hits ' write(27,*) 'plane #hit WG straw letediff drdist' do jj=1,24 if (nhit(jj) .ne. spdetector.l.fpp.plane(jj).nghits) then write(27,*) 'nhits .ne. nghits' endif do ii=1,nhit(jj) write(27,'(3(3x,i3),2x,f5.1,2x,f7.1,2x,f8.2)') jj,ii, 1 spdetector.l.fpp.plane(jj).cor.wg(ii), 1 spdetector.l.fpp.plane(jj).cor.straw(ii), 1 spdetector.l.fpp.plane(jj).cor.letediff(ii), 1 spdetector.l.fpp.plane(jj).cor.drdist(ii) enddo enddo endif endif else c endif c endif c return end