subroutine CHI_SQUARED(npar,iflag,t_index) c c This routine calculates a chi squared value for a set of data and a c function. In the process this routine returns the parameters thaty best c fit the data. c c Implementation of a different fitting package: Singular Value c Decomposition. Jan 01 1997 c --- Ketevi A. ASSAMAGAN c c Implementation of a robust fitting algorithm with Huber's rho or c Andrew's sine using the Multidimentional Simplex Method. c Sep 01 1997 c --- Ketevi A. Assamgan c c c READ THE NOTE IN THE FILE "ucorrect.f" BEFORE ATTEMPTING ANY FIT. c --- Ketevi A. Assamagan c c c INPUTS: c NPAR = number of currently variable parameters c PAR = array of constant and variable parameters c IFLAG= indicates what is to be calculated, c 1 -> initialization c 2 -> calculate a gradient, only used with set gradient command c 3 -> wrap everything up c FUTIL= function to be called c c OUTPUTS: c func_value = function value returned c gradient = array of first derivitives c c ?? TPW c implicit none include 'db.inc' include 'data.inc' include 'sieve.inc' include 'target.inc' include 'fit.inc' c real*8 sig(max_points) real*8 tmpfunc(max_points,MAX_NUM_PAR),ytmp(max_points) real*8 del(max_points),res(max_points) real*8 ptm(MAX_NUM_PAR),err(MAX_NUM_PAR),dchi,chisq_df,eprime,kine_fac real*8 ooptan,scattan,gasdev real*8 zdrift(max_points),pi real*8 chisq,chi,diffx,diffy,dist real*8 ddummy,ydummy,tdummy,pdummy real*8 hsyptar,hsxptar,scatt_ang,azim_ang real*4 r0,fictitious_radius,model,s1,s2,select integer*4 npar ! # of free parameters integer*4 i,k,l integer*4 iflag, t_index ! what to do integer*4 j,nx,ny,mfit,listp(MAX_NUM_PAR) ! index variables character*9 today ! today's date character*8 timenow ! current time c parameter (pi = 3.14159265358979323846) c write(*,*) ' in chisq_hms iflag = ',iflag open(unit=21,file='fit.log',status='unknown', 1 access='append') c first_type =.TRUE. if (iflag .eq. 1) then ! start life write(*,*) ' calling mapper_init' call mapper_init ! initialize mappings write (*,*) ' Initialization complete' endif if (iflag .eq. 2) then ! gradient calc write (*,*) ' Gradient type of calculation not supported' endif if ((npar.eq.0).or.(t_index.eq.1)) then write(*,*) 'No Fitting Requested ...' write(*,*) 'Mapping Data with Starting Matrix Elements !!!' write(*,*) ' ' call mapper_wrapup(npar) write (*,*) ' Wrapping up ... wait, ntotal = ',ntotal first_type=.FALSE. do i=1, ntotal if ( mod(i,1000) .eq. 0) write(*,*) ' at event =',i ni = i call mapper(pts(i)) enddo npass = npass + ipass return endif c if (icoll.eq.1.and.(t_index.eq.4.or.t_index.eq.2)) then type *, 'You Need Sieve Slit data to fit this quantity...' return endif c if (npass.eq.1) then open(unit=20,file='fit_parameters.out',status='unknown') else open(unit=20,file='fit_parameters.out',status='unknown', 1 access='append') write(20,*) ' ' endif **** added for pionct 5/2/05 open(unit=25,file='ben_oldvar.out',status='unknown') open(unit=26,file='ben_newvar.out',status='unknown') **** end pionct c open(unit=30,file='fit_matrix.db',status='unknown') c Write(*,*) 'Getting the Sample of Data to Fit...Patience ntotal =',ntotal c call date(today) call time(timenow) WRITE(6,'(1X,''Starting Date = '',A9, '' Starting Time = '',A8)') 1 today,timenow WRITE(21,'(1X,''Starting Date = '',A9, '' Starting Time = '',A8)') 1 today,timenow WRITE(21,'(1X,''Number of events considered = '',I8)') 1 ntotal c c Get the nominal errors (sigma) on each data point. Later on, these are smeared with c a Gaussian if least square svdfit is selected. --- Ketevi A. Assamagan c if (ifit.eq.0) then call uerrors(tdummy,pdummy,ydummy,ddummy) else if (t_index.eq.5) a0 = test_delta if (t_index.eq.4) a0 = test_phi1 if (t_index.eq.3) a0 = test_y_tar if (t_index.eq.2) a0 = test_theta1 endif c c Start looping over events. For each one, use the current matrix elements in c memory (this is either the starting elements or the elements generated in the c previous fit) to map focal plane quantities to target, then to slit and find (x,y) c of the event at the slit. Using the information in sieveslit.db, decide which slit c (defined by the slit numbers nx, ny) is closest to the event (x,y). --- Ketevi A. Assamagan c ndata = 0 k = 0 write(*,*) ' ntotal = ',ntotal do i = 1, ntotal ni = i ! loop over the data ndata = ndata + 1 ! go map the data k = k +1 call mapper (pts(i)) c c write(*,*) ' In chisq t_index = ',t_index if (t_index.eq.5) then c write(*,*) ' call func_fill ndata =',ndata, ' i = ',i call func_fill(pts(i)) c c Calculate delta-p from kinematics, subtract the excitation energy and compare c with the model or previous fit-reconstruted value. Select the sample of events c to fit based on this deviation. Ketevi A. Assamagan c c hsyptar = dble(phi(ndata)) c hsxptar = dble(theta(ndata)) c call scattering_angle(spec_ang,hsyptar,hsxptar,scatt_ang,azim_ang) c c kine_fac = 1 + ((2.0*ebeam-eloss)/target_mass)*((dsin(scatt_ang/2.0))**2) c eprime = (1/kine_fac)*(ebeam-eloss/2-eexc*(1+eexc/target_mass)) - eloss/2 c del(ndata)=(sngl(eprime)-deldata(i))/deldata(i) c tmp(ndata) = del(ndata) tmp(ndata) = deldata(i)/100. if (ifit.eq.0) then do j = 1, npar tmpfunc(ndata,j) = func(ndata,j) enddo sig(ndata) = gasdev(ddummy,t_index,i) c if (abs(del(ndata)-delta_t(ndata)).gt.sngl(test_delta)) c 1 ndata=ndata-1 c write(*,*) ' ndata = ',ndata,tmp(ndata),delta_t(ndata) else model = 0.0 do j=1,npar model = model + parm(j)*func(ndata,j) enddo res(k) = tmp(ndata) - model if (res(k).eq.0) k = k-1 endif endif c if (t_index.ge.2.and.t_index.le.4) then c theta(ndata) = theta_t(ndata) phi(ndata) = phi_t(ndata) call func_fill(pts(i)) ! this fills array func c if (t_index.eq.3) then c c find the ytarget of each event. c ytmp(ndata) = ydata(i) - Zzero(i)*phi(ndata) tmp(ndata) = ytmp(ndata) if (ifit.eq.0) then sig(ndata) = gasdev(ydummy,t_index,i) do j = 1, npar tmpfunc(ndata,j) = func(ndata,j) enddo else model = 0.0 do j=1,npar model = model + parm(j)*func(ndata,j) enddo res(k) = tmp(ndata) - model if (res(k).eq.0) k = k-1 endif c else c c if (t_index.eq.2) then tmp(ndata) = thdata(ndata) if (ifit.eq.0) then sig(ndata) = gasdev(tdummy,t_index,i) do j = 1, npar tmpfunc(ndata,j) = func(ndata,j) enddo else model = 0.0 do j=1,npar model = model + parm(j)*func(ndata,j) enddo res(k) = tmp(ndata) - model if (res(k).eq.0) k = k-1 endif endif c if (t_index.eq.4) then tmp(ndata) = phidata(ndata) if (ifit.eq.0) then sig(ndata) = gasdev(pdummy,t_index,i) do j = 1, npar tmpfunc(ndata,j) = func(ndata,j) enddo else model = 0.0 do j=1,npar model = model + parm(j)*func(ndata,j) enddo res(k) = tmp(ndata) - model if (res(k).eq.0) k = k-1 endif endif endif endif enddo !end looping over events c c Now, we are getting ready to do the fit: the sample data was collected above for the c target variable that we want fit. First, store the values of the fitting fuction foe c events. Second, identify which parameters we want to vary. Third, surrender control to c the fitting routine. Fourth, when the fit is done, calculate the chi-square. Fifth, with c the new matrix elements fromthe fit, remap the events to the target, then to the slits to c show the results of the fit (this is stored in the ntuple id 200 in pawc directory. When c all this is done, end the calculations and return control to paw so the user can study c the results of his/her work. ---Ketevi A. Assamagan c mfit = 0 do j=1,npar ptm(j) = parm(j) ! save initial parameters if (ind(j).ne.0) then mfit = mfit + 1 listp(mfit) = j ! array index of fitted parameters endif enddo c **** added for pionct 5/2/05 do i = 1, ntotal model = 0.0 do j=1,mfit model = model + parm(listp(j))*func(i,j) enddo write(25,*)' ',model ! write initial ytarg to file enddo close (25) **** end pionct if (ifit.eq.0) then c write(*,*)'number of events which pass the test',ndata write(*,*)'number of nonzero matixelements',npar write(21,*)'number of events which pass the test',ndata write(21,*)'number of nonzero matixelements',npar c c svdfit does the fit by singular value decomposition. You may want to study the c routine svdfit.f. c --- Ketevi A. Assamagan c call svdfit(tmpfunc,tmp,sig,ndata,npar,mfit,listp,err,dchi,tol) c c Calculate the chi square now that new parameters are generated. c --- Ketevi A. Assamagan c chisq = 0.0 do i=1,ndata model = 0.0 do j=1,npar model = model + parm(j)*func(i,j) enddo res(i) = tmp(i) - model chi = ( res(i)/sig(i) )**2 chisq = chisq + chi **** added for pionct 5/2/05 model = 0.0 do j=1,mfit *parm contains the new matrix elements fit by svd model = model + parm(listp(j))*func(i,j) enddo write(26,*)' ',model ! write final ytarg to file **** end pionct enddo close (26) chisq_df = chisq/(ndata-npar) c do j=1,npar write(20,22) ptm(j),parm(j) ! write initial and fitted parameters to file enddo 22 format(10x,2(4x,SPE15.7)) close (20) c WRITE(6,'(1X,''CHISQ/D.F. = '',E15.7, '' NPARMS = '',I3)') 1 chisq_df,npar WRITE(21,'(1X,''CHISQ/D.F. = '',E15.7, '' NPARMS = '',I3)') 1 chisq_df,npar write(*,*) 'Number of passes', npass write(21,*) 'Number of passes', npass else c write(*,*)'number of events with non zero residual',k write(*,*)'number of nonzero matixelements',npar write(21,*)'number of events with non zero residual',k write(21,*)'number of nonzero matixelements',npar c c if you have any other fitting package, insert the source code here. c note that the fitting package should return the array parm(npar) in c the common block defined in fit.inc. --- Ketevi A. Assamgan c c Robust estimate by downhill simplex method in multidimension using c Huber and Andrew's functions. --- Ketevi A. Assamagan c npr = npar type *, 'Searching for a robust measure of spread ...' l = k/2 if ( (k-2*l).eq.1 ) then l = (k+1)/2 scale = select(l,k,res,max_points) else l = k/2 s1 = select(l,k,res,max_points) l = k/2 + 1 s2 = select(l,k,res,max_points) scale = 0.5*(s1+s2) endif scale = scale/0.6745 type *, 'Robust estimate of scale = ', scale mstep = 0 do while (mstep.le.ipass) call robust_estimate c c Calculate the chi square now that new parameters are generated. c --- Ketevi A. Assamagan c chisq = 0.0 do i=1,ntotal model = 0.0 do j=1,npar model = model + parm(j)*func(i,j) enddo res(i) = tmp(i) - model chi = ( res(i)/scale )**2 chisq = chisq + chi enddo chisq_df = chisq/(ntotal-npar) do j=1,npar write(20,23) ptm(j),parm(j) ! write initial and fitted parameters to file enddo 23 format(10x,2(4x,SPE15.7)) WRITE(6,'(1X,''CHISQ/D.F. = '',E15.7, '' NPARMS = '',I3)') 1 chisq_df,npar WRITE(21,'(1X,''CHISQ/D.F. = '',E15.7, '' NPARMS = '',I3)') 1 chisq_df,npar enddo npass = mstep write(*,*) 'Number of passes', npass write(21,*) 'Number of passes', npass close (20) endif c c Goodness of fit. Not implemented at all. The defaut is just a return. The user c is free to any action in this routine ugoodness.f. c ---- Ketevi A. Assamagan c c call ugoodness(df,dchi,cl) c c Map events with the optimized matrix elements and produce a new ntuple id = 200 in c pawc directory with results reflecting the effects of the new parameters. c --- Ketevi A. Assamagan c call mapper_wrapup(npar) c call me_db_writefile(30,'set') close(30) c write(*,*) ' ' if (npass.eq.ipass) then write (*,*) ' Wrapping up ... wait' first_type=.FALSE. do i=1, ntotal ni = i call mapper(pts(i)) enddo c call date(today) call time(timenow) WRITE(6,'(1X,''Ending Date = '',A9, '' Ending Time = '',A8)') 1 today,timenow WRITE(21,'(1X,''Ending Date = '',A9, '' Ending Time = '',A8)') 1 today,timenow write(21,*) ' ******************************* ' close (21) endif end