source: branches/start/MagicSoft/Simulation/Corsika/Mmcs/longft.f@ 18569

      SUBROUTINE LONGFT(FPARAM,CHI2)
 
C-----------------------------------------------------------------------
C  LONG(ITUDINAL) F(I)T
C
C  THIS ROUTINE PERFORMS A FIT TO THE LONGITUDINAL DISTRIBUTION OF AN
C  AIR SHOWER. DUE TO THE LARGE PARTICLE NUMBERS IN AN AIR SHOWER THE
C  STATISTICAL ERRORS ON THE PARTICLE NUMBER AT A GIVEN LEVEL ARE
C  MINUTE. THIS LEADS TO RATHER LARGE CHI**2/DOF FOR THE FITS EVEN IF
C  THE FITTED FUNCTION MATCHES THE POINTS BETTER THAN SAY 1%.
C  KEEP IN MIND THAT FITTING IS A DIFFICULT TASK AND THE RESULTS DO NOT
C  NECESSARILY REPRESENT THE ABOLUTE MINIMUM OR EVEN A GOOD
C  APPROXIMATION.
C
C  TRY A 6 PARAMETER FIT BASED ON J. BALL'S PROPOSED CURVE REPLACING HIS
C  CONSTANT WIDTH PARAMETER LAMBDA BY A POLYNOMIAL OF 3. DEGREE.
C   N(T) = NMAX * ((T-T0)/(TMAX-T0))**((TMAX-T)/(P1+P2*T+P3*T**2))
C   T    = DEPTH IN G/CM**2
C   T0   = STARTING DEPTH OF SHOWER
C   TMAX = DEPTH OF SHOWER MAXIMUM
C   NMAX = PARTICLE NUMBER AT TMAX
C   P1 .. P3 = PARAMETERS OF A POLYNOMIAL DESCRIBING THE WIDTH
C
C  THIS SUBROUTINE IS CALLED FROM MAIN
C-----------------------------------------------------------------------
 
      IMPLICIT NONE
*KEEP,CURVE.
      COMMON /CURVE/   CHAPAR,DEP,ERR,NSTP
      DOUBLE PRECISION CHAPAR(1100),DEP(1100),ERR(1100)
      INTEGER          NSTP
*KEEP,RUNPAR.
      COMMON /RUNPAR/  FIXHEI,THICK0,HILOECM,HILOELB,
     *                 STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
     *                 MONIOU,MDEBUG,NUCNUC,
     *                 CETAPE,
     *                 SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
     *                 N1STTR,MDBASE,
     *                 DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
     *                 FIX1I,FMUADD,FNKG,FPRINT,FDBASE
     *                ,GHEISH,GHESIG
      COMMON /RUNPAC/  DSN,HOST,USER
      DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
      REAL             STEPFC
      INTEGER          NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
     *                 SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
     *                 N1STTR,MDBASE
      INTEGER          CETAPE
      CHARACTER*79     DSN
      CHARACTER*20     HOST,USER
 
      LOGICAL          DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
     *                 FIX1I,FMUADD,FNKG,FPRINT,FDBASE
     *                ,GHEISH,GHESIG
*KEND.
 
      INTEGER NPAR
      PARAMETER (NPAR=6)
      DOUBLE PRECISION F(NPAR),FPARAM(NPAR),CHI2,CHISQ
      DOUBLE PRECISION P(NPAR+1,NPAR),Y(NPAR+1),EPS
      DOUBLE PRECISION T0,TMAX,NMAX,FAC
      INTEGER          I,J,JJ,K,ITER,IFLAG
      EXTERNAL         CHISQ
C-----------------------------------------------------------------------
 
      IF ( DEBUG ) WRITE(MDEBUG,*) 'LONGFT:'
 
C  FIND GOOD START VALUES FOR XMAX AND FMAX
      NMAX = 0.D0
      DO 2 I = 1,NSTP
        ERR(I) = MAX( 1.D0, SQRT(CHAPAR(I)) )
        IF ( CHAPAR(I) .GT. NMAX ) THEN
          NMAX = CHAPAR(I)
          TMAX = DEP(I)
        ENDIF
 2    CONTINUE
C  STARTVALUE FOR X0 IS ABOUT WHERE MORE THAN 1 PARTICLE SHOWS UP
      DO 3 I = 1,NSTP
        IF ( CHAPAR(I) .GT. 1.D0 ) GOTO 1
 3    CONTINUE
      I = NSTP
 1    CONTINUE
      T0 = DEP(I)
 
C-----------------------------------------------------------------------
C  FIT IS PERFORMED WITH THE ROUTINE AMOEBA FROM:
C      NUMERICAL RECIPES, W.H. PRESS ET AL.,
C      CAMBRIDGE UNIVERSITY PRESS, 1992  ISBN 0 521 43064 X
C  SEE THERE HOW IT HAS TO BE USED.
 
C  CREATE A SET OF NPAR+1 STARTING VERTICES
C  HERE IS THE FIRST ONE
      P(1,1) = NMAX
      P(1,2) = T0
      P(1,3) = TMAX
      P(1,4) = 200.D0
      P(1,5) = 1.D-1
      P(1,6) = 1.D-1
 
C  LOOP OVER THE FITTING ROUTINE (2 TIMES 5 FITS WITH VARYING PRECISION)
      DO 10 J = 1,2
        DO 9 JJ = 1,5
C  START WITH CRUDE PRECISION AND IMPROVE STEP BY STEP
C  AFTER FIVE STEPS ENLARGE AGAIN
          EPS = 10.D0**(-3.D0-JJ*0.5D0)
          FAC = 1.D0 + 2.D0**(2.1D0*(1.D0-JJ))
C  GO AS WELL IN DIFFERENT DIRECTIONS
          IF ( J .EQ. 2 ) FAC = 1.D0/FAC
 
C  GET OTHER NPAR STARTING VERTICES FROM THE STARTING POINT BY VARIATION
C  OF ONLY ONE OF THE COORDINATE VALUES
          DO 5 I = 2,NPAR+1
            DO 4 K = 1,NPAR
              P(I,K) = P(1,K)
 4          CONTINUE
            IF ( P(I,I-1) .EQ. 0.D0 ) THEN
              P(I,I-1) = 1.D0
            ELSE
              P(I,I-1) = P(I,I-1) * FAC
            ENDIF
 5        CONTINUE
          IF (DEBUG) WRITE(MDEBUG,*) 'LONGFT: TRIAL,FAC,EPS ',J,
     *                                   SNGL(FAC),SNGL(EPS)
 
C  CALCULATE FUNCTION VALUES AT THE START VERTICES
          DO 7 I = 1,NPAR+1
            DO 6 K = 1,NPAR
              F(K) = P(I,K)
 6          CONTINUE
            Y(I) = CHISQ(F)
 7        CONTINUE
C  PERFORM A FIT
          CALL AMOEBA(P,Y,NPAR+1,NPAR,NPAR,EPS,CHISQ,ITER,IFLAG)
          IF ( DEBUG ) THEN
            WRITE(MDEBUG,*) 'LONGFT: ITER/IFLAG=',ITER,IFLAG
            WRITE(MDEBUG,*) 'LONGFT: PARAMETERS=',1,(P(1,K),K=1,6)
            WRITE(MDEBUG,*) 'LONGFT: CHISQ     =',SNGL(Y(1))
          ENDIF
 
C  STORE VALUES AT FIRST TRIAL OR AT IMPROVED RESULT
          IF ( J .EQ. 1 .OR. Y(1) .LT. CHI2 ) THEN
            DO 8 I = 1,NPAR
              FPARAM(I) = P(1,I)
 8          CONTINUE
            CHI2 = Y(1)
          ENDIF
C  END OF LOOPS OVER THE FITTING ROUTINE
 9      CONTINUE
 10   CONTINUE
 
      RETURN
      END