program deteff

*     1. linearer Fit (A+xB) einer Detektoreichmessung              *
*     2. Normierung der Daten mit der verfeinerten Funktion         *
*     3. Speicherung im Format fuer deteff.dat                      *
*     27.4.99 keller, Aenderung 15.7.99,28.9.99                     *
*     Aenderung 3.11.99: schneidet Punkte die 0 sind weg (HRPT!)    * 
*               7.3.00,24.5.02,9.8.02                               *
*     16.01.03 zolliker (Fragt, ob untergrund flach sein soll)      *
*                                                                   *
*     > myfit -o deteff deteff.f                                    *

      implicit none

      external FIT_LIN_FUN

      real par(16),err(16),int(2000),tth(2000),par1,ti,step,tf,lim
      integer n,i,nret,l,count
      character*36 spec
      character*5 flag,flag2
      character*78 title
	real flat

      call fit_init

! Function title and parameter names
	call fit_userfun('STRAIGHT LINE', fit_lin_fun)	! function title, function
	call fit_userpar('Bg(0)')		! first parameter: background at zero
	call fit_userpar('dBg/dX')		! second parameter: slope

        call fit_dat_options('cal=0')   ! turns off calibration of raw data

      do i=1,2000
        int(i)=0.
      enddo


      call sys_getenv('dat_defspec',spec)
      if (spec.eq.' ') spec='DMC'
      call str_trim(spec,spec,l)            ! Laenge bestimmen (l)

  100 write(*,*)
      write(*,'(x,a)')'DMC = 1 / HRPT = 2'
      write(*,'(x,3a,$)')'Instrument (default: ',spec(1:l),'): '
      read(*,'(a)')flag                     !      ^ schreibt spec von
                                            !        Zeichen 1 bis l

      call str_upcase(flag,flag)           ! schreibt Inhalt von flag gross
      if (flag.eq.'1') flag='DMC'
      if (flag.eq.'2') flag='HRPT'
      if (flag.ne.' ') then
        if (flag.ne.'DMC' .and. flag.ne.'HRPT') goto 100
        spec=flag
      endif

      write(*,*)
      write(*,'(x,a,$)') 'Flat Background [y]'
      read(*,'(a)') flat
      call str_upcase(flat, flat)

      call sys_setenv('dat_defspec',spec)


C      WRITE (*,'(X,A,$)') 'Calibration data files (e.g.:  DMC/1999/
C     A4-34) :  '
C      READ (*,'(A)') FILES

C      CALL FIT_DAT(FILES)

C      call fit_dat(' ')
C      call fit_merge(0.02)
      call fit_dat_merge(' ',0.02)         ! ersetzt fit_dat und fit_merge
      call fit_get_array('X',tth,2000,nret)
      call fit_get_array('Y',int,2000,nret)

      write(*,*)
      write(*,'(x,a,$)')'Title of detector efficiency file: '
      read(*,'(a)')title

      if (title.eq.' ') then
        if (nret.eq.400) title='DMC detector efficiency file'
        if (nret.eq.1600) title='HRPT detector efficiency file'
      endif

      if (nret.ne.400) then
        if (nret.ne.1600) then
          write(*,*)
          write(*,*)
          write(*,*)
          write(*,*)
          write(*,*)'  **************** WARNING ***************'
          write(*,*)
          write(*,*)'  The number of data points does not match '
          write(*,*)'   the number of detectors of DMC or HRPT'
          write(*,*)
          write(*,*)
          write(*,'(a,$)')' continue? (y/n): '
          read(*,'(a)')flag2
          call str_upcase(flag2,flag2)   
          if (flag2.ne.'Y') then
            write(*,'(a)')' program aborted'
            goto 102
          endif
        endif
      endif

      TI=TTH(1)
      TF=TTH(NRET)
      STEP=(TTH(NRET)-TTH(1))/(NRET-1.)

      PAR1=0
      DO I=1,NRET
        PAR1=PAR1+INT(I)
      ENDDO
      PAR(1)=PAR1/NRET                          ! starting value for Bg(0)
      PAR(2)=0.                                 ! starting value for slope
      ERR(1)=1.
      ERR(2)=1.

      lim=par(1)/4.                             ! points below lim are excluded!

      call fit_exclude(0.,0.,0.,lim)

      CALL FIT_FUN(7,2,PAR,ERR)

      if (flat .eq. 'Y') call fit_fix(2)        ! fix slope

      CALL FIT_FIT(0)
      CALL FIT_FIT(0)
      CALL FIT_FIT(0)
      CALL FIT_GET_ARRAY('P',PAR,2,N)           ! reads fit parameters
      CALL FIT_GET_ARRAY('E',ERR,2,N)           ! reads parameter errors

      OPEN (1,FILE='deteff.dat',STATUS='UNKNOWN')

      WRITE (1,'(2A)')'#',TITLE
      if (nret.gt.999) write(1,'(I4)')nret
      if (nret.lt.1000) write(1,'(I3)')nret

      count=0                                   ! count: # of excluded points

      DO 101 I=1,NRET

        if (int(i).lt.lim) then
          write(1,'(F8.6)')0.0
          count=count+1
        else
          WRITE(1,'(F8.6)')INT(I)/(PAR(1)+((I-1)*STEP+TI)*PAR(2))
        endif

  101 CONTINUE 

      CLOSE(1)

      call str_trim(spec,spec,l)
      WRITE(*,*)
      WRITE(*,*)
      WRITE(*,*)
      WRITE(*,*)
      WRITE(*,*)'  New detector efficiency file:  deteff.dat'
      WRITE(*,'(3a)')'   Copy this file into the directory /data/
     Alnslib/data/',spec(1:l),'/<year>/'
      WRITE(*,*)
      WRITE(*,*)'  Normalization function: P(1)+2th*P(2)'
      WRITE(*,'(A,F12.5,A,F10.5,A,F10.5,A,F9.5,A)')'   P(1)='
     A,PAR(1),'(',ERR(1),'), P(2)=',PAR(2),'(',ERR(2),')'
      WRITE(*,*)
      if (count.gt.0) write (*,*)'  Number of excluded points: ',count
      WRITE(*,*)

  102 END


	real function fit_lin_fun(x,p,n,mode,cinfo)
!	-------------------------------------------

	implicit none

	real x		! x-value
	integer n	! number of parameters
	real p(n)	! parameters
	integer mode	! mode
	integer cinfo	! calculation information (see below)

	if (mode .eq. 0) then

! Define here your own function

	  fit_lin_fun=p(1)+x*p(2)

	elseif (mode .lt. 0) then

! Use this part to do some initialisations.
! (e.g. read files, write out comments on your user function)
! This section is called by FIT_FUN (command FUN)

	  type *
	  type *, 'to define your own user function leave FIT and type MYFIT'
	  type *, 'Example: STRAIGHT LINE'

	endif
	end