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