I found typos in the attached fortran files used in the
optimizations of the optics matrix by espace which
needed to be corrected. The place
of the corrections in marked by 'mkj' and is commented
out with the replacement line(s) following.
Mark
c---------------------------------------------------------------------------
c Copyright (c) 1996 Southeastern Universities Research Association,
c Thomas Jefferson National Accelerator Facility
c
c This software was developed under a United States Government
c license described in the NOTICE file included as part of this
c distribution.
c
c Eddy A.J.M. Offermann, 12000 Jefferson Ave., Newport News, VA 23606
c Email: offerman@cebaf.gov Tel: (757) 249-7544 Fax: (757) 249-5800
c---------------------------------------------------------------------------
c
subroutine curfit_transp(idef,fitpar,unitnr,linear,xd,xt,xr,y,
> sigmay,yfit,dp_kin_ho,costh,sinth,tanth,
> pcentr,nlow,nhigh,a,
> dmin,dmax,thmin,thmax,phmin,phmax,
> ymin,ymax,sigmaa,nterms,itermax,chimin,
> conv,chisqr,nrfree,ifail)
c
c----------------------------------------------------------------------
c
c the subroutine curfit performs a least squares fit to the data, using
c a chi-square goodness-of-fit criterion based on gaussian statistics.
c the algorithm is based on the curfit routine of Bevington.
implicit none
include 'parameter.h'
include 'espace_type.h'
real*8 sigmaa(NPAR1,*),dp_kin_ho(*),costh(*),sinth(*),tanth(*),
> pcentr(*)
record /fit_param/ fitpar
record /matr_foc_tar/ a(*)
record /tar_event/ y(*),yfit(*)
record /foc_event/ xd(*),xt(*),xr(*)
record /vertex/ sigmay(*)
integer idef,nlow,nhigh,nterms,itermax,unitnr,nrfree,
> dmin,dmax,thmin,thmax,phmin,phmax,ymin,ymax,ifail
logical linear
real*8 fchisq_transp
integer pntr2
logical allocate_mem,deallocate_mem
record /matr_foc_tar/ a_int(NPAR1),b(NPAR1)
real*8 beta(NPAR1),alpha(NPAR1,NPAR1),wt1(*),delwt1(*),wt2(*),
> delwt2(*),deriv1(*),deriv2(*),array(NPAR,NPAR),
> chisq1,chisqr,denom,temp,chimin,conv,flamda,ajj,corr
integer nits,nlam,mpts,i,j,k,jj,mvar,indvc(NPAR1,PORDER),
> frvc(NPAR1*PORDER),id000,p1,p2
logical final,convrg,statusok
character*50 fname(3)
save fname
pointer (pderiv1,deriv1),(pderiv2,deriv2),
> (pwt1,wt1),(pwt2,wt2),(pdelwt1,delwt1),(pdelwt2,delwt2)
data fname(1) /' Dmnp*th^m*y^n*ph^p,kincor*ph'/
data fname(2) /' Ymnp*th^m*y^n*ph^p'/
data fname(3) /' Tmnp*th^m*y^n*ph^p,Pmnp*th^m*y^n*ph^p'/
c nterms = total number of transfer elements used
c mvar = total number of parameters fitted
c nits = total number of iterations
c flamda = projection parameter
c nlam = number of times flambda increased with no decrease
c in chisquare
ifail = 0
c Set up independent variable index array (indvc) based
c on vary codes a(i).vc. For each parameter, this array contains.
c The index of the free parameter to which it corresponds or
c is locked. It is used in computing the derivatives of the fitting
c function with respect to its free parameters. A parameter which
c is either fixed or locked to one that is fixed has an index of
c zero. Setup also an array (frvc) which gives the index "i" of
c a variable "a(i)" with varycode 1 if the indvc is known.
mvar = 0
do i = 1,nterms
a_int(i) = a(i)
do j = 1,PORDER
indvc(i,j) = 0
frvc(i*j) = 0
enddo
enddo
do i = 1,nterms
if (a(i).vc.ge.1) then
do j = 1,a(i).vc
mvar = mvar+1
indvc(i,j) = mvar
frvc(mvar) = i
enddo
endif
enddo
c Initialize memory for the arrays DERIV1 and DERIV2
statusok = allocate_mem(8*nhigh*mvar,pderiv1,deriv1)
> .and.
> allocate_mem(8*nhigh*mvar,pderiv2,deriv2)
> .and.
> allocate_mem(8*nhigh,pwt1,wt1)
> .and.
> allocate_mem(8*nhigh,pwt2,wt2)
> .and.
> allocate_mem(8*nhigh,pdelwt1,delwt1)
> .and.
> allocate_mem(8*nhigh,pdelwt2,delwt2)
if (.not.statusok) then
stop 'CURFIT_TRANSP: Error in memory allocation'
endif
c If the chisquare fit is going to be linear in the fitted
c parameters, make sure that FLAMDA is negligible.
if (linear) then
flamda = 1.0d-30
else
flamda = 1.0d-3
endif
final = .false.
nits = 0
nlam = 0
mpts = nhigh-nlow+1
nrfree = mpts-mvar
if (nrfree.le.0) then
ifail = 1
call derror('CURFIT_TRANSP: NFREE < 0 ?')
return
endif
id000 = 0
if (fitpar.option.eq.1) then
do i = 1,nterms
if (a(i).n.eq.'D000') id000 = i
enddo
if (id000.eq.0) then
ifail = 1
call derror('CURFIT_TRANSP: D000 is not being fitted ?')
return
endif
endif
c ***** compute initial estimate of fitting function *****
if (idef.eq.1) then
call functn_transp1(fitpar.mode,xd,xt,xr,dp_kin_ho,costh,sinth,
> tanth,pcentr,nlow,nhigh,
> a,nterms,dmin,dmax,thmin,thmax,ymin,ymax,
> phmin,phmax,fitpar.option,id000,yfit)
else if (idef.eq.2) then
call functn_transp2(fitpar.mode,xd,xt,dp_kin_ho,costh,sinth,
> tanth,pcentr,nlow,nhigh,
> a,nterms,dmin,dmax,thmin,thmax,ymin,ymax,
> phmin,phmax,fitpar.option,id000,yfit)
endif
chisq1 = fchisq_transp(y,sigmay,nlow,nhigh,nrfree,yfit,
> fitpar.option)
chisqr = chisq1
write(unitnr,'(/''fitfunction: '',a50)') fname(fitpar.option)
if (linear) then
write(unitnr,'(14x,''linear chisquare fit'')')
else
write(unitnr,'(14x,''non-linear chisquare fit'')')
endif
write(unitnr,'(/''initial reduced chisquare ='',1pe12.4)') chisqr
write(unitnr,'(/''nits reduced chisquare total chisquare'',7x,
> ''flambda'')')
c#######################################################################
c###### start main loop #######
c ***** compute alpha and beta matrices and the modified *****
c ***** curvature martrix "array". *****
100 if (flamda.eq.0.0) final = .true.
if (.not.final.and.nlam.eq.0) then
do i = 1,mvar
beta(i) = 0.0
do j = 1,i
alpha(i,j) = 0.0
enddo
enddo
c... Calculate the derivatives of the fitting function with respect
c to its parameters.
if (idef.eq.1) then
call fderiv_transp1(fitpar.mode,xd,xt,xr,dp_kin_ho,costh,
> sinth,tanth,pcentr,nlow,nhigh,a,mvar,
> dmin,dmax,thmin,thmax,ymin,ymax,
> phmin,phmax,indvc,fitpar.option,id000,
> deriv1,deriv2)
else if (idef.eq.2) then
call fderiv_transp2(fitpar.mode,xd,xt,dp_kin_ho,costh,
> sinth,tanth,pcentr,nlow,nhigh,a,mvar,
> dmin,dmax,thmin,thmax,ymin,ymax,
> phmin,phmax,indvc,fitpar.option,id000,
> deriv1,deriv2)
endif
c... Compute weights for the alpha and beta matrices.
c... gaussian statistics
do i = nlow,nhigh
if (fitpar.option.eq.1) then
wt1(i) = 1./sigmay(i).dp**2
delwt1(i) = (y(i).dp_kin-yfit(i).dp_kin)*wt1(i)
else if (fitpar.option.eq.2) then
wt1(i) = 1./sigmay(i).y**2
delwt1(i) = (y(i).y-yfit(i).y)*wt1(i)
else if (fitpar.option.eq.3) then
wt1(i) = 1./sigmay(i).th**2
delwt1(i) = (y(i).th-yfit(i).th)*wt1(i)
wt2(i) = 1./sigmay(i).ph**2
delwt2(i) = (y(i).ph-yfit(i).ph)*wt2(i)
endif
enddo
c... Compute beta
do j = 1,mvar
if (fitpar.option.eq.1.or.fitpar.option.eq.2) then
do i = nlow,nhigh
p1 = pntr2(i,j,nhigh)
beta(j) = beta(j)+delwt1(i)*deriv1(p1)
enddo
else if (fitpar.option.eq.3) then
if (frvc(j).ge.1.and.frvc(j).le.NMATR_FOC_TAR1) then
do i = nlow,nhigh
p2 = pntr2(i,j,nhigh)
beta(j) = beta(j)+delwt1(i)*deriv1(p2)+
> delwt2(i)*deriv2(p2)
enddo
else if (frvc(j).ge.thmin.and.frvc(j).le.thmax) then
do i = nlow,nhigh
p2 = pntr2(i,j,nhigh)
beta(j) = beta(j)+delwt1(i)*deriv1(p2)
enddo
else if (frvc(j).ge.phmin.and.frvc(j).le.phmax) then
do i = nlow,nhigh
p2 = pntr2(i,j,nhigh)
beta(j) = beta(j)+delwt2(i)*deriv2(p2)
enddo
endif
endif
c... Compute alpha
do k = 1,j
if (fitpar.option.eq.1.or.fitpar.option.eq.2) then
do i = nlow,nhigh
p1 = pntr2(i,j,nhigh)
p2 = pntr2(i,k,nhigh)
alpha(j,k) = alpha(j,k)+
> wt1(i)*deriv1(p1)*deriv1(p2)
enddo
else if (fitpar.option.eq.3) then
if (frvc(j).ge.1.and.frvc(j).le.NMATR_FOC_TAR1) then
if ((frvc(k).ge.1.and.frvc(k).le.NMATR_FOC_TAR1)
> .or.
> (frvc(k).ge.thmin.and.frvc(k).le.thmax)) then
do i = nlow,nhigh
p1 = pntr2(i,j,nhigh)
p2 = pntr2(i,k,nhigh)
alpha(j,k) = alpha(j,k)+
> wt1(i)*deriv1(p1)*deriv1(p2)
enddo
endif
if ((frvc(k).ge.1.and.frvc(k).le.NMATR_FOC_TAR1)
> .or.
c mkj
c > (frvc(k).ge.thmin.and.frvc(k).le.thmax)) then
c replace by
> (frvc(k).ge.phmin.and.frvc(k).le.phmax)) then
do i = nlow,nhigh
p1 = pntr2(i,j,nhigh)
p2 = pntr2(i,k,nhigh)
alpha(j,k) = alpha(j,k)+
> wt2(i)*deriv2(p1)*deriv2(p2)
enddo
endif
else if (frvc(j).ge.thmin.and.frvc(j).le.thmax) then
if ((frvc(k).ge.1.and.frvc(k).le.NMATR_FOC_TAR1)
> .or.
> (frvc(k).ge.thmin.and.frvc(k).le.thmax)) then
do i = nlow,nhigh
p1 = pntr2(i,j,nhigh)
p2 = pntr2(i,k,nhigh)
alpha(j,k) = alpha(j,k)+
> wt1(i)*deriv1(p1)*deriv1(p2)
enddo
endif
else if (frvc(j).ge.phmin.and.frvc(j).le.phmax) then
if ((frvc(k).ge.1.and.frvc(k).le.NMATR_FOC_TAR1)
> .or.
> (frvc(k).ge.phmin.and.frvc(k).le.phmax)) then
do i = nlow,nhigh
p1 = pntr2(i,j,nhigh)
p2 = pntr2(i,k,nhigh)
alpha(j,k) = alpha(j,k)+
> wt2(i)*deriv2(p1)*deriv2(p2)
enddo
endif
endif
endif
alpha(k,j) = alpha(j,k)
enddo
enddo
endif
c... Compute curvature matrix
do j = 1,mvar
ajj = alpha(j,j)
do k = 1,mvar
denom = abs(ajj*alpha(k,k))
array(j,k) = 0.0
if (denom.gt.0.) array(j,k) = alpha(j,k)/sqrt(denom)
enddo
array(j,j) = 1.0+flamda
enddo
call matinv(array,mvar,NPAR)
c-----------------------------------------------------------------------
c ***** compute new parameter estimates *****
if (.not.final) then
jj = 0
do j = 1,nterms
b(j) = a(j)
if (a(j).vc.ge.1) then
if (fitpar.mode.eq.0) then
jj = jj+1
ajj = alpha(jj,jj)
do k = 1,mvar
temp = alpha(k,k)*ajj
if (j.eq.id000.and.fitpar.option.eq.1) then
if (temp.gt.0.0) b(j).m(1) = b(j).m(1)+
> beta(k)*array(jj,k)/sqrt(temp)
else
if (temp.gt.0.0) b(j).v = b(j).v+
> beta(k)*array(jj,k)/sqrt(temp)
endif
enddo
else
do i = 1,a(j).vc
jj = jj+1
ajj = alpha(jj,jj)
do k = 1,mvar
temp = alpha(k,k)*ajj
if (temp.gt.0.0) b(j).m(i) = b(j).m(i)+
> beta(k)*array(jj,k)/sqrt(temp)
enddo
enddo
endif
endif
enddo
call protect_transp(b,a_int,nterms,fitpar.option)
c-----------------------------------------------------------------------
c ***** compute new estimate of fitting function and *****
c ***** test for convergence. *****
if (idef.eq.1) then
call functn_transp1(fitpar.mode,xd,xt,xr,dp_kin_ho,
> costh,sinth,tanth,pcentr,nlow,nhigh,b,
> nterms,dmin,dmax,thmin,thmax,ymin,ymax,
> phmin,phmax,fitpar.option,id000,yfit)
else if (idef.eq.2) then
call functn_transp2(fitpar.mode,xd,xt,dp_kin_ho,
> costh,sinth,tanth,pcentr,nlow,nhigh,b,
> nterms,dmin,dmax,thmin,thmax,ymin,ymax,
> phmin,phmax,fitpar.option,id000,yfit)
endif
chisq1 = fchisq_transp(y,sigmay,nlow,nhigh,nrfree,yfit,
> fitpar.option)
write(unitnr,'(i4,5x,1pe12.4,6x,e12.4,3x,e12.4)')
> nits,chisq1,nrfree*chisq1,flamda
c... Compare new chisquare with that of previous fit.
if ((chisqr-chisq1).ge.0.0) then
do j = 1,nterms
a(j) = b(j)
enddo
nits = nits+1
nlam = 0
if (flamda.ge.1.0e-4) flamda = flamda/10.0
else
flamda = flamda*10.0
nlam = nlam+1
endif
c... If convergence has been obtained with the present set of
c varycodes, stop the fit.
convrg = ((nits.ge.itermax) .or.
> (min(chisqr,chisq1).lt.chimin) .or.
> (nlam.ge.10).or.(flamda.gt.1.e6) .or.
> (((chisqr/chisq1-1.0).le.conv) .and.
> ((chisqr/chisq1-1.0).ge.0.0)))
convrg = convrg.or.linear
if (convrg) then
flamda = 0.0
endif
chisqr = min(chisqr,chisq1)
endif
if (.not.final) goto 100
c#####################################################################
c ***** compute proper error matrix and the uncertainties *****
c ***** in the parameters. *****
do j = 1,mvar
ajj = alpha(j,j)
do k = 1,mvar
denom = abs(ajj*alpha(k,k))
if (denom.le.0.) then
array(j,k) = 0.
if (j.eq.k) array(j,k) = 1.0
else
array(j,k) = array(j,k)/sqrt(denom)
endif
enddo
enddo
c... calculate the uncertainties in the parameters.
jj = 0
do j = 1,nterms
do i = 1,PORDER
if (a(j).vc.ge.i) then
jj = jj+1
sigmaa(j,i) = sqrt(abs(array(jj,jj)))
else
sigmaa(j,i) = 0.0
endif
enddo
enddo
c... check the correlation coefficients between the free parameters,
c a warning should be given if they are larger than 0.95.
do j = 1,mvar
do k = 1,j-1
corr = array(j,k)/sqrt(array(j,j)*array(k,k))
if (abs(corr).gt.0.95) then
write(unitnr,'(''WARNING: correlation between '',a,
> '' and '',a,'' > 0.95'',1pe10.3)')
> a(frvc(j)).n,a(frvc(k)).n,corr
endif
enddo
enddo
c... deallocate memory
statusok = deallocate_mem(pderiv1)
> .and.
> deallocate_mem(pderiv2)
> .and.
> deallocate_mem(pwt1)
> .and.
> deallocate_mem(pwt2)
> .and.
> deallocate_mem(pdelwt1)
> .and.
> deallocate_mem(pdelwt2)
return
end
integer function pntr2(i,j,n1)
c This function calculates for a two dimensional array of n1 rows
c the pointer to position where it is stored in a one dim array.
implicit none
integer i,j,n1
pntr2 = i+n1*(j-1)
return
end
c---------------------------------------------------------------------------
c Copyright (c) 1996 Southeastern Universities Research Association,
c Thomas Jefferson National Accelerator Facility
c
c This software was developed under a United States Government
c license described in the NOTICE file included as part of this
c distribution.
c
c Eddy A.J.M. Offermann, 12000 Jefferson Ave., Newport News, VA 23606
c Email: offerman@cebaf.gov Tel: (757) 249-7544 Fax: (757) 249-5800
c---------------------------------------------------------------------------
c
subroutine fderiv_transp1(fit_mode,xd,xt,xr,dp_kin_ho,costh,sinth,
> tanth,pcentr,nlow,nhigh,
> a,mvar,dmin,dmax,thmin,thmax,ymin,ymax,
> phmin,phmax,indvc,ifunc,id000,
> deriv1,deriv2)
c
c----------------------------------------------------------------------
c
implicit none
include 'parameter.h'
include 'espace_type.h'
include 'espace_user.h'
include 'optimize.h'
include 'experiment.h'
include 'kinema.h'
include 'fderiv_param.h'
record /foc_event/ xd(*),xt(*),xr(*)
record /matr_foc_tar/ a(*)
integer nlow,nhigh,dmin,dmax,thmin,thmax,ymin,ymax,phmin,phmax,
> indvc(NPAR1,PORDER),mvar,ifunc,id000,fit_mode
real*8 dp_kin_ho(*),deriv1(*),deriv2(*),costh(*),sinth(*),
> tanth(*),pcentr(*)
real*8 pow,calc_var_tg,foc_deriv
integer pntr2
integer i,j,k,p1,ifoc,io,focmin,focmax
real*8 const,cos_rho,cos_rho2,tan_rho,cos_rho_loc,cos_rho_loc2,
> tan_rho_loc,ph_tg,dp,dp_kin_shift,dx,dy,pmom
record /vertex/ fpow1(0:MAXPAR),fpow2(-1:MAXPAR),
> dfoc_coeff(3,PORDER)
c Reset all the derivatives
do j = 1,mvar
do i = nlow,nhigh
p1 = pntr2(i,j,nhigh)
deriv1(p1) = 0.0
deriv2(p1) = 0.0
enddo
enddo
focmin = 1
focmax = NMATR_FOC_TAR1
tan_rho = a(1).m(1)
cos_rho = 1./sqrt(1.0+tan_rho**2)
cos_rho2 = 1./(1.0+tan_rho**2)
do i = nlow,nhigh
if (a(1).vc.ge.1) then
xt(i).th = (xd(i).th+tan_rho)/(1-xd(i).th*tan_rho)
dx = xt(i).th*tan_rho
xt(i).x = xd(i).x*cos_rho*(1.+dx)
xr(i).x = xt(i).x
endif
if (fit_mode.eq.0) then
c It is extremely important that one calculates D000 event-by-event
c in FIT_MODE=0. This matrix element is indicated by ID000.
if (ifunc.eq.1)
> call calc_matrix(id000,id000,xr(i).x,a)
else
call calc_matrix(focmin,focmax,xr(i).x,a)
call calc_matrix(thmin,thmax,xr(i).x,a)
call calc_matrix(phmin,phmax,xr(i).x,a)
call calc_matrix(ymin,ymax,xr(i).x,a)
call calc_matrix(dmin,dmax,xr(i).x,a)
endif
tan_rho_loc = a(1).v
cos_rho_loc = 1./sqrt(1.0+tan_rho_loc**2)
cos_rho_loc2 = 1./(1.0+tan_rho_loc**2)
if (a(1).vc.ge.1.or.a(2).vc.ge.1.or.a(3).vc.ge.1) then
xt(i).ph = xd(i).ph/(1-xd(i).th*tan_rho)/cos_rho
dy = tan_rho*cos_rho*xt(i).ph*xd(i).x
xt(i).y = xd(i).y+dy
xr(i).th = (xd(i).th+a(1).v)/(1-xd(i).th*a(1).v)
xr(i).ph = (xd(i).ph-a(3).v)/(1-xd(i).th*a(1).v)/cos_rho_loc
xr(i).y = xt(i).y-a(2).v
endif
do j = -1,MAXPAR
fpow2(j).x = pow(xr(i).x,j)
fpow2(j).th = pow(xr(i).th,j)
fpow2(j).y = pow(xr(i).y,j)
fpow2(j).ph = pow(xr(i).ph,j)
if (j.ge.0) fpow1(j) = fpow2(j)
enddo
c derivatives of x,theta,y and phi at the focal plane wrt tf000,
c yf000 and pf000.
c The dependence of the rotated focal plane coordinates (XR..) on
c the tf000, yf000 and pf000 is as follows:
c xr(i).x : tf000(1) |
c xr(i).th: tf000 | The (1) indicates that it depends
c xr(i).y : tf000(1),p000,y000| on only the first coefficient.
c xr(i).th: tf000,pf000 |
c
c The array DFOC_COEFF(I,J) contains these derivatives, where
c I points to 1: tf000,
c 2: yf000,
c 3: pf000
c and J points to the order in xr(i).x dependence.
c d../dtf000
do io = 1,PORDER
dfoc_coeff(1,io).th =
> (1.0+xd(i).th**2)/(1.0-xd(i).th*tan_rho_loc)**2
if (io.eq.1) then
dfoc_coeff(1,1).x = xd(i).x*cos_rho*(xt(i).th*cos_rho2+
> dfoc_coeff(1,1).th*tan_rho-tan_rho*cos_rho2)
else
dfoc_coeff(1,io).x = 0.
endif
dfoc_coeff(1,io).ph =
> xr(i).ph*(cos_rho_loc2*tan_rho_loc+xd(i).th/
> (1.-xd(i).th*tan_rho_loc))
if (io.eq.1) then
dfoc_coeff(1,io).y = xd(i).x*cos_rho2*cos_rho*
> xt(i).ph*(1.+tan_rho*xt(i).th)
else
dfoc_coeff(1,io).y = 0.
endif
c d../dyf000
dfoc_coeff(2,io).th = 0.
dfoc_coeff(2,io).x = 0.
dfoc_coeff(2,io).ph = 0.
dfoc_coeff(2,io).y = -1.
c d../dpf000
dfoc_coeff(3,io).th = 0.
dfoc_coeff(3,io).x = 0.
dfoc_coeff(3,io).ph =
> -1./(1.-xd(i).th*tan_rho_loc)/cos_rho_loc
dfoc_coeff(3,io).y = 0.
enddo
dp = calc_var_tg(dmin,dmax,a,fpow1)
pmom = pcentr(i)*(1+dp)
ph_tg = calc_var_tg(phmin,phmax,a,fpow1)
if (particle_name.eq.'lepton') then
dp_kin_shift = pmom/ma*sinth(i)*ph_tg
else if (particle_name.eq.'hadron') then
dp_kin_shift = (tanth(i)+2*mom_e0**2*costh(i)*sinth(i)/
> ((mom_e0+mx)**2-(mom_e0*
> costh(i))**2))*ph_tg
else
dp_kin_shift = 0.0
endif
c /* >>>>> corrected momentum <<<<< */
if (ifunc.eq.1) then
do j = dmin,dmax
do ifoc = 1,3
do k = 1,a(ifoc).vc
if (indvc(ifoc,k).ne.0) then
p1 = pntr2(i,indvc(ifoc,k),nhigh)
deriv1(p1) = deriv1(p1)+(1+dp_kin_shift)*
> foc_deriv(j,k,a,fpow2,dfoc_coeff(ifoc,k),xr(i))
endif
enddo
enddo
do k = 1,a(j).vc
if (indvc(j,k).ne.0) then
p1 = pntr2(i,indvc(j,k),nhigh)
deriv1(p1) = (1+dp_kin_shift)*
> fpow2(a(j).pw(1)).th*
> fpow2(a(j).pw(2)).y*
> fpow2(a(j).pw(3)).ph*
> fpow2(k-1).x
endif
enddo
enddo
const = dp_kin_shift/ph_tg
do j = phmin,phmax
do ifoc = 1,3
do k = 1,a(ifoc).vc
if (indvc(ifoc,k).ne.0) then
p1 = pntr2(i,indvc(ifoc,k),nhigh)
deriv1(p1) = deriv1(p1)+const*(1+dp)*
> foc_deriv(j,k,a,fpow2,dfoc_coeff(ifoc,k),xr(i))
endif
enddo
enddo
do k = 1,a(j).vc
if (indvc(j,k).ne.0) then
p1 = pntr2(i,indvc(j,k),nhigh)
deriv1(p1) = const*(1+dp)*
> fpow2(a(j).pw(1)).th*
> fpow2(a(j).pw(2)).y*
> fpow2(a(j).pw(3)).ph*
> fpow2(k-1).x
endif
enddo
enddo
c /* >>>>> y target <<<<< */
else if (ifunc.eq.2) then
do j = ymin,ymax
do ifoc = 1,3
do k = 1,a(ifoc).vc
if (indvc(ifoc,k).ne.0) then
p1 = pntr2(i,indvc(ifoc,k),nhigh)
deriv1(p1) = deriv1(p1)+
> foc_deriv(j,k,a,fpow2,dfoc_coeff(ifoc,k),xr(i))
endif
enddo
enddo
do k = 1,a(j).vc
if (indvc(j,k).ne.0) then
p1 = pntr2(i,indvc(j,k),nhigh)
deriv1(p1) = fpow2(a(j).pw(1)).th*
> fpow2(a(j).pw(2)).y*
> fpow2(a(j).pw(3)).ph*
> fpow2(k-1).x
endif
enddo
enddo
c /* >>>>> theta and phi target <<<<< */
else if (ifunc.eq.3) then
do j = thmin,thmax
do ifoc = 1,3
do k = 1,a(ifoc).vc
if (indvc(ifoc,k).ne.0) then
p1 = pntr2(i,indvc(ifoc,k),nhigh)
deriv1(p1) = deriv1(p1)+
> foc_deriv(j,k,a,fpow2,dfoc_coeff(ifoc,k),xr(i))
endif
enddo
enddo
do k = 1,a(j).vc
if (indvc(j,k).ne.0) then
p1 = pntr2(i,indvc(j,k),nhigh)
deriv1(p1) = fpow2(a(j).pw(1)).th*
> fpow2(a(j).pw(2)).y*
> fpow2(a(j).pw(3)).ph*
> fpow2(k-1).x
endif
enddo
enddo
do j = phmin,phmax
do ifoc = 1,3
do k = 1,a(ifoc).vc
if (indvc(ifoc,k).ne.0) then
p1 = pntr2(i,indvc(ifoc,k),nhigh)
c mkj corrections : changed followling two lines to lines between ***
c deriv1(p1) = deriv1(p1)+
c > foc_deriv(j,k,a,fpow2,dfoc_coeff(ifoc,k),xr(i))
c ***
deriv2(p1) = deriv2(p1)+
> foc_deriv(j,k,a,fpow2,dfoc_coeff(ifoc,k),xr(i))
c ***
endif
enddo
enddo
do k = 1,a(j).vc
if (indvc(j,k).ne.0) then
p1 = pntr2(i,indvc(j,k),nhigh)
deriv2(p1) = fpow2(a(j).pw(1)).th*
> fpow2(a(j).pw(2)).y*
> fpow2(a(j).pw(3)).ph*
> fpow2(k-1).x
endif
enddo
enddo
endif
enddo
return
end
real*8 function foc_deriv(jcoef,k,coef,fpow,dpar,x)
implicit none
include 'parameter.h'
include 'espace_type.h'
include 'fderiv_param.h'
integer jcoef,k
record /matr_foc_tar/ coef(*)
record /vertex/ fpow(-1:MAXPAR),dpar
record /foc_event/ x
foc_deriv =
> fpow(coef(jcoef).pw(1)-1).th*
> fpow(coef(jcoef).pw(2)-1).y*
> fpow(coef(jcoef).pw(3)-1).ph*
> fpow(k-1).x*
> (coef(jcoef).v*coef(jcoef).pw(1)*x.y*x.ph*dpar.th+
> coef(jcoef).v*coef(jcoef).pw(2)*x.th*x.ph*dpar.y+
> coef(jcoef).v*coef(jcoef).pw(3)*x.th*x.y*dpar.ph+
> coef(jcoef).m(k)*(k-1)*fpow(k-2).x*x.th*x.y*x.ph*dpar.x)
return
end