SUBROUTINE HDPM C----------------------------------------------------------------------- C H(ADRONIC) D(UAL) P(ARTON) M(ODEL) C C GENERATOR OF HADRONIC COLLISION INSPIRED BY DUAL PARTON MODEL C THIS SUBROUTINE IS CALLED FROM SDPM C----------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) *KEEP,CONST. COMMON /CONST/ PI,PI2,OB3,TB3,ENEPER DOUBLE PRECISION PI,PI2,OB3,TB3,ENEPER *KEEP,DPMFLG. COMMON /DPMFLG/ NFLAIN,NFLDIF,NFLPI0,NFLCHE,NFLPIF,NFRAGM INTEGER NFLAIN,NFLDIF,NFLPI0,NFLCHE,NFLPIF,NFRAGM *KEEP,ELADPM. COMMON /ELADPM/ ELMEAN,ELMEAA,IELDPM,IELDPA DOUBLE PRECISION ELMEAN(37),ELMEAA(37) INTEGER IELDPM(37,13),IELDPA(37,13) *KEEP,ELASTY. COMMON /ELASTY/ ELAST,IELIS,IELHM,IELNU,IELPI DOUBLE PRECISION ELAST INTEGER IELIS(20),IELHM(20),IELNU(20),IELPI(20) *KEEP,INDICE. COMMON /INDICE/ NNUCN,NKA0,NHYPN,NETA,NETAS,NPIZER, * NNC,NKC,NHC,NPC,NCH,NNN,NKN,NHN,NET,NPN INTEGER NNUCN(2:3),NKA0(2:3),NHYPN(2:3),NETA(2:3,1:4), * NETAS(2:3),NPIZER(2:3), * NNC,NKC,NHC,NPC,NCH,NNN,NKN,NHN,NET,NPN *KEEP,INTER. COMMON /INTER/ AVCH,AVCH3,DC0,DLOG,DMLOG,ECMDIF,ECMDPM,ELAB, * FNEUT,FNEUT2,GNU,PLAB,POSC2,POSC3,POSN2,POSN3, * RC3TO2,S,SEUGF,SEUGP,SLOG,SLOGSQ,SMLOG, * WIDC2,WIDC3,WIDN2,WIDN3,YCM,YY0,ZN, * IDIF,ITAR DOUBLE PRECISION AVCH,AVCH3,DC0,DLOG,DMLOG,ECMDIF,ECMDPM,ELAB, * FNEUT,FNEUT2,GNU,PLAB,POSC2,POSC3,POSN2,POSN3, * RC3TO2,S,SEUGF,SEUGP,SLOG,SLOGSQ,SMLOG, * WIDC2,WIDC3,WIDN2,WIDN3,YCM,YY0,ZN INTEGER IDIF,ITAR *KEEP,ISTA. COMMON /ISTA/ IFINET,IFINNU,IFINKA,IFINPI,IFINHY INTEGER IFINET,IFINNU,IFINKA,IFINPI,IFINHY *KEEP,LEPAR. COMMON /LEPAR/ LEPAR1,LEPAR2,LASTPI,NRESPC,NRESPN,NCPLUS INTEGER LEPAR1,LEPAR2,LASTPI,NRESPC,NRESPN,NCPLUS *KEEP,MULT. COMMON /MULT/ EKINL,MSMM,MULTMA,MULTOT DOUBLE PRECISION EKINL INTEGER MSMM,MULTMA(37,13),MULTOT(37,13) *KEEP,NEWPAR. COMMON /NEWPAR/ EA,PT2,PX,PY,TMAS,YR,ITYP, * IA1,IA2,IB1,IB2,IC1,IC2,ID1,ID2,IE1,IE2,IF1,IF2, * IG1,IG2,IH1,IH2,II1,II2,IJ1,NTOT DOUBLE PRECISION EA(3000),PT2(3000),PX(3000),PY(3000),TMAS(3000), * YR(3000) INTEGER ITYP(3000), * IA1,IA2,IB1,IB2,IC1,IC2,ID1,ID2,IE1,IE2,IF1,IF2, * IG1,IG2,IH1,IH2,II1,II2,IJ1,NTOT *KEEP,PAM. COMMON /PAM/ PAMA,SIGNUM DOUBLE PRECISION PAMA(6000),SIGNUM(6000) *KEEP,PARPAR. COMMON /PARPAR/ CURPAR,SECPAR,PRMPAR,OUTPAR,C, * E00,E00PN,PTOT0,PTOT0N,THICKH,ITYPE,LEVL DOUBLE PRECISION CURPAR(14),SECPAR(14),PRMPAR(14),OUTPAR(14), * C(50),E00,E00PN,PTOT0,PTOT0N,THICKH INTEGER ITYPE,LEVL *KEEP,PARPAE. DOUBLE PRECISION GAMMA,COSTHE,PHI,H,T,X,Y,CHI,BETA,GCM,ECM EQUIVALENCE (CURPAR(2),GAMMA), (CURPAR(3),COSTHE), * (CURPAR(4), PHI ), (CURPAR(5), H ), * (CURPAR(6), T ), (CURPAR(7), X ), * (CURPAR(8), Y ), (CURPAR(9), CHI ), * (CURPAR(10),BETA), (CURPAR(11),GCM ), * (CURPAR(12),ECM ) *KEEP,RANDPA. COMMON /RANDPA/ FAC,U1,U2,RD,NSEQ,ISEED,KNOR DOUBLE PRECISION FAC,U1,U2 REAL RD(3000) INTEGER ISEED(103,10),NSEQ LOGICAL KNOR *KEEP,RATIOS. COMMON /RATIOS/ RPI0R,RPIER,RPEKR,RPEKNR,PPICH,PPINCH,PPNKCH, * ISEL,NEUTOT,NTOTEM DOUBLE PRECISION RPI0R,RPIER,RPEKR,RPEKNR,PPICH,PPINCH,PPNKCH INTEGER ISEL,NEUTOT,NTOTEM *KEEP,RESON. COMMON /RESON/ RDRES,RESRAN,IRESPAR REAL RDRES(2),RESRAN(1000) INTEGER IRESPAR *KEEP,REST. COMMON /REST/ CONTNE,TAR,LT DOUBLE PRECISION CONTNE(3),TAR INTEGER LT *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. C----------------------------------------------------------------------- IF ( DEBUG ) WRITE(MDEBUG,444) (CURPAR(I),I=1,9) 444 FORMAT(' HDPM : CURPAR=',1P,9E10.3) C SET ANTI-LEADER TO PROTON OR NEUTRON; TARGET IS ALWAYS NUCLEON CALL RMMAR( RD,1,1 ) IF ( RD(1) .LE. CONTNE(LT) ) THEN ITAR = 13 ELSE ITAR = 14 ENDIF C CALCULATE LAB AND CM ENERGY IF ( ITYPE .NE. 1 ) THEN ELAB = PAMA(ITYPE) * GAMMA PLAB = ELAB * BETA S = PAMA(ITYPE)**2 + PAMA(ITAR)**2 + 2.D0*PAMA(ITAR)*ELAB ELSE C FOR GAMMA-INDUCED REACTION TAKE PI(0) AS LEADING PARTICLE ITYPE = 7 ELAB = GAMMA PLAB = ELAB S = PAMA(ITAR)**2 + 2.D0*PAMA(ITAR)*ELAB ENDIF ECMDPM = SQRT(S) IF ( DEBUG ) WRITE(MDEBUG,*) 'HDPM : ITYPE,ELAB,PLAB,S,ECMDPM=', * ITYPE,SNGL(ELAB),SNGL(PLAB),SNGL(S),SNGL(ECMDPM) C LN(S), LN(S)**2 AND RAPIDITY OF C. M. SYSTEM IN LAB SLOG = LOG(S) SLOGSQ = SLOG**2 SMLOG = LOG( 2.D0 * PAMA(ITAR) * ELAB ) ELABLG = LOG(ELAB) EPLUSP = ELAB + PLAB * YCM = 0.5D0 * LOG( (ELAB+PAMA(ITAR)+PLAB)/(ELAB+PAMA(ITAR)-PLAB) ) YCM = 0.5D0 * LOG( (EPLUSP**2 +PAMA(ITAR)*EPLUSP)/ * (PAMA(ITYPE)**2+PAMA(ITAR)*EPLUSP) ) IF ( DEBUG ) WRITE(MDEBUG,*)'HDPM : SLOG,SLOGSQ,YCM=', * SNGL(SLOG),SNGL(SLOGSQ),SNGL(YCM) C----------------------------------------------------------------------- C RETURN POINT IF CALCULATION OF PARTICLES GOES WRONG 1 CONTINUE IF ( ITYPE .NE. 7 ) THEN C CHOOSE NUMBER OF INTERACTIONS IN TARGET CALL TARINT ELSE C FOR GAMMA-INDUCED REACTIONS TAKE ALWAYS ONE COLLISION GNU = 1.D0 ENDIF C----------------------------------------------------------------------- C NO DIFFRACTION IF C OR THE NUMBER OF INTERACTIONS IN TARGET IS CHOSEN RANDOMLY C AND MORE THAN ONE INTERACTION TAKES PLACE C OR PRIMARY PARTICLE IS GAMMA (PI0) C NOW NFLDIF DECIDES WHETHER DIFFRACTIVE PROCESS POSSIBLE OR NOT IF ( ( NFLAIN.EQ.0 .AND. GNU.GT.1.D0 .AND. NFLDIF.EQ.0 ) * .OR. ( ITYPE .EQ. 7 ) ) THEN IDIF = 0 ELSE C SET DIFFRACTION FLAG IF RANDOM NUMBER < PROBABILITY CALL RMMAR( RD,1,1 ) C IDIF IS 0 : NO DIFFRACTION ; IDIF IS 1 : DIFFRACTION C DIFFRACTION RISES WITH ENERGY AND SATURATES AT 10000 GEV C ### DAS TUT ES ABER NICHT: ES IST KONSTANT 0.15 (SIEHE DPFUNC) !!!! IF ( RD(1) .GT. DPFUNC(ECMDPM) ) THEN IDIF = 0 ELSE IDIF = 1 ENDIF ENDIF C PRINTOUT FOR DEBUG IF ( DEBUG ) THEN WRITE(MDEBUG,*) ' DIFFRACTIVE INTERACTION (0/1) = ',IDIF ENDIF C SET COUNTER FOR REPEAT TO 0 NREPRD = 0 C GENERATION OF INTERACTION 1919 CONTINUE C FLAG TO CHECK NUMBER OF SECONDARIES; C IS CHANGED TO 1 IF SECONDARY MULTIPLICITY IS LOW ISEL = 0 C SET LEADING PARTICLE TO INCOMING PARTICLE AND ANTI-LEADER TO NUCLEON C (AS IT COMES FROM TARGET NUCLEUS) BOTH MAY BE CHANGED BY LEPACX LEPAR1 = ITYPE LEPAR2 = ITAR IF ( IDIF .EQ. 0 ) THEN C----------------------------------------------------------------------- C NON SINGLE DIFFRACTIVE PROCESS STARTS HERE CALL NSD C CHARGE EXCHANGE ENABLED? EXCHANGE LEADER AND ANTI-LEADER LASTPI = 0 NRESPC = 0 NRESPN = 0 NCPLUS = 0 IF ( NFLCHE .EQ. 0 ) THEN CALL LEPACX( ECMDPM,ELABLG,LEPAR1,1 ) CALL LEPACX( ECMDPM,ELABLG,LEPAR2,2 ) ENDIF 1921 CONTINUE CALL RNEGBI( NCH,AVCH,ECMDPM ) C NCH IS # OF ALL CHARGED PARTICLES INCLUDING EXCESS FROM TARGET IF ( NCH .LT. 1 ) THEN IF ( LEPAR1 .LT. 50 .OR. LEPAR2 .LT. 50 ) THEN NREPRD = NREPRD + 1 IF ( NREPRD .GT. 10 ) GOTO 1 GOTO 1921 ELSE C INTERACTION IS ONLY RESONANCE PRODUCTION ISEL = 1 ENDIF ENDIF C WIDTH PLATEAU FOR CLUSTERS AND FOR CALCULATION OF CENTR.RAP.DENSITY DELRAP = 0.6722D0 * (2.95D0 + 0.0302D0 * SLOG) C SET RSLOG FOR CALCULATION OF PARTICLE RATIOS RSLOG = SLOG C AVERAGE TRANSVERSE MOMENTUM CALL AVEPT( ECMDPM,SLOG ) ELSE C----------------------------------------------------------------------- C SINGLE DIFFRACTIVE PROCESS STARTS HERE 1920 CONTINUE CALL DIFRAC( NRETDF ) IF ( NRETDF .EQ. 1 ) GOTO 1 C CHARGE EXCHANGE ENABLED? EXCHANGE CHARGE OF DIFFRACTING PARTICLE LASTPI = 0 NRESPC = 0 NRESPN = 0 NCPLUS = 0 IF ( NFLCHE .EQ. 0 ) THEN IF ( YY0 .GT. 0.D0 ) THEN C PROJECTILE DIFFRACTION CALL LEPACX( ECMDIF,DMLOG,LEPAR1,1 ) ELSE C TARGET DIFFRACTION CALL LEPACX( ECMDIF,DMLOG,LEPAR2,2 ) ENDIF ENDIF C FLUCTUATION OF MULTIPLICITY ACCORDING TO NEG.BIN. DISTRIBUTION CALL RNEGBI( NCH,AVCH,ECMDIF ) C REPEAT CALCULATION AS SOMETHING WENT WRONG IF ( NCH .LT. 1 ) THEN IF ( (YY0 .GT. 0.D0 .AND. LEPAR1 .LT. 50) .OR. * (YY0 .LT. 0.D0 .AND. LEPAR2 .LT. 50) ) THEN NREPRD = NREPRD + 1 IF ( NREPRD .GT. 10 ) GOTO 1 GOTO 1920 ELSE C DIFFRACTIVE INTERACTION IS ONLY RESONANCE PRODUCTION ISEL = 1 ENDIF ENDIF C SET RSLOG FOR CALCULATION OF PARTICLE RATIOS RSLOG = DLOG C HERE WE USE ECMDPM, BECAUSE THE MOMENTUM TRANSFER IS DEPENDENT C ON THE ENERGY OF THE TOTAL SYSTEM AND NOT ON THE DIFFRACTING MASS CALL AVEPT( ECMDPM,SLOG ) ENDIF C----------------------------------------------------------------------- C NOW FOR NSD AND DIFFRACTIVE PROCESSES C IN CASE OF LOW MULTIPLICITY SET FLAG ISEL IF ( NCH .LE. 2 ) ISEL=1 C FNCH IS FLUCTUATING TOT.NUMBER OF CHARGED PARTICLES FOR ALL 3 STRINGS FNCH = DBLE(NCH) C RATIO ALL CHARGED PARTICLES WITH FLUCTUATION/WITHOUT FLUCTUATION XZ = FNCH / AVCH C FNCH3 IS FLUCTUATING NUMBER OF CHARGED PARTICLES FOR 3RD STRING FNCH3 = XZ * AVCH3 C FNCH2 IS FLUCTUATING NUMBER OF CHARGED PARTICLES 1ST AND 2ND STRING FNCH2 = FNCH - FNCH3 C RC3TO2 IS RATIO (CHARGED 3RD STRING)/(CHARGED 1ST AND 2ND STRING) IF ( FNCH2 .NE. 0.D0 ) THEN RC3TO2 = FNCH3 / FNCH2 ELSE RC3TO2 = 0.D0 ENDIF IF ( DEBUG ) WRITE(MDEBUG,*) ' FNCH,FNCH2,FNCH3,RC3TO2=', * SNGL(FNCH),SNGL(FNCH2),SNGL(FNCH3),SNGL(RC3TO2) C IS NUMBER OF NEUTRALS FLUCTUATING AS NUMBER OF CHARGED ? IF ( NFLPIF .EQ. 0 .OR. IDIF .EQ. 1 .OR. ECMDPM .LT. 60.D0 ) THEN C SET NUMBER OF GAMMAS ACCORDING TO NEG. BIN. VARIABLE XZ C AS NUMBER OF NEUTRALS FLUCTUATES AS CHARGED. SEUGF = SEUGP * XZ ZG = XZ ELSE C NFLPIF IS 1 MEANS: # OF PI(0) FLUCTUATES AS MEASURED AT COLLIDER IF ( ECMDPM .LT. 200.D0 ) THEN SEUGF = SEUGP * XZ * SEUGF = (0.0786D0*SLOG-0.010D0)*FNCH2 + (0.391D0*SLOG+0.305D0) ELSE C DETERMINE NEW NUMBER OF GAMMAS WITH FLUCTUATION AROUND SEUGP*XZ AGR = EXP(-XZ) DGR = SEUGP * XZ * (0.9823D0 - 0.3756D0 * AGR) SGS = DGR * (0.14718D0 + 2.53213D0 * AGR) 723 CONTINUE SEUGF = 0.88D0 * RANNOR(DGR,SGS) IF ( SEUGF .LT. 1.D0 ) GOTO 723 ENDIF C SET NEGATIVE BINOMIAL VARIABLE ZG FOR GAMMAS ZG = SEUGF / SEUGP ENDIF SEUGF = MAX( 1.D0, SEUGF ) IF ( DEBUG ) WRITE(MDEBUG,*) 'HDPM :XZ,ZG,SEUGF=', * SNGL(XZ),SNGL(ZG),SNGL(SEUGF) C----------------------------------------------------------------------- C RATIO ALL-NUCLEON/ALL-CHARGED C PARAMETRISATION FROM UA5, NUCL. PHYS. B291 (1987) 445 EQ.(2.4) RNUCCH = MAX( 0.D0, -0.008D0 + 0.00865D0 * RSLOG ) C NUMBER FOR DIRECT NEUTRON/ANTINEUTRON PRODUCTION 1ST AND 2ND STRING C MULTIPLY BY 0.5 BECAUSE RATIO RNUCCH GIVES (ALL-NUCL)/(ALL-CHARGED) C AND HERE ONLY THE NEUTRON-ANTINEUTRONS ARE COUNTED FNUCN = 0.5D0 * RNUCCH * FNCH2 C RATIO (ALL CHARGED SIGMAS)/(ALL CHARGED) IS 1/3 OF ALL STRANGE BARYON C PARAMETRISATION FORM UA5, NUCL. PHYS. B291 (1987) 445 EQ.(2.5) RHYPCH = MAX( 0.D0, (-0.007D0 + 0.0028D0 * RSLOG) * OB3 ) C NEUTRAL STRANGE BARYONS ARE DOUBLE OF CHARGED STRANGE BARYONS FHYPN = 2.D0 * RHYPCH * FNCH2 C CORRECT NUMBER OF GAMMAS FROM NEUTRAL HYPERON DECAY S0-->L+GAMMA SEUGFC = MAX( 0.D0, SEUGF - 0.5D0 * FHYPN ) C RATIO CHARGED-KAON/CHARGED PIONS C PARAMETRISATION FROM UA5, NUCL. PHYS. B291 (1987) 445 EQ.(2.7) RKPI = MAX (0.D0, 0.024D0 + 0.0062D0 * RSLOG ) C RKCH IS RATIO (CHARGED-KAON)/(ALL-CHARGED) DERIVED FROM RKPI; C THE FACTOR 0.5 IN FRONT OF RNUCCH IS BECAUSE ONLY HALF OF NUCLEONS C ARE P/PBAR. THE 1.17 IS AN APROXIMATE CONVERSION FACTOR FROM C (ALL-NUCL)/(ALL-CHARGED) TO (ALL-NUCL)/(CHARGED-PI), WHICH IS A BIT C ENERGY DEPENDENT (1.14 ...1.21) SEE GEICH-GIMBEL TABLE 7.1 RKCH = RKPI / (1.D0 + RKPI + (0.5D0*RNUCCH+RHYPCH) * 1.17D0) C K0/K0-BAR FOR 1ST AND 2ND STRING C NEUTRAL KAONS ARE PRODUCED WITH THE SAME RATE AS CHARGED KAONS FKA0 = RKCH * FNCH2 C RATIO ETA/PI(0) IS ASSUMED TO BE INDEPENDENT OF ENERGY = 0.19 C SEE: ANSORGE ET AL. (UA5-COLLABORATION) Z.PHYS.C43(1989)75 * RETPI0 = 0.19D0 C RATIO ETA/PI(0) IS ASSUMED TO BE DEPENDENT ON ENERGY C SEE: GEICH-GIMBEL,INT.J.MOD.PHYS.A4(1989)1527 TAB.7.1 C HECK'S FIT: RETPI0 IS 0.06 + 0.006*RSLOG + 0.0011*RSLOG**2 RETPI0 = 0.06D0 + 0.006D0 * RSLOG + 0.0011D0 * RSLOG**2 C AUXIL1 IS FRACTION OF PI(0)/(PI(0)+ETA) AUXIL1 = 1.D0 / (1.D0 + RETPI0) C NUMBER OF GAMMAS FROM PI(0) IS 2, FROM ETA IS 3.216 IN AVERAGE; C AUXIL2 IS NUMBER OF GAMMA-PRODUCING-PARTICLES: PI(0) AND ETA AUXIL2 = SEUGFC / ( AUXIL1 * 2.D0 + (1.D0 - AUXIL1) * 3.216D0 ) FETA = (1.D0 - AUXIL1) * AUXIL2 FPI0 = AUXIL1 * AUXIL2 C CORRECT FPI0 BY DECAYS OF STRANGE BARYONS; NEUTRAL: FHYPN*0.357 C CHARGED: 0.5*FNCH2*RHYPCH*0.5157; IT YIELDS FHYPN*(0.357+0.12893) FPI0 = MAX( 0.D0, FPI0 - FHYPN * 0.486D0 ) C SUMMED NEUTRAL PARTICLES FOR 1ST AND 2ND STRING FNEUT2 = FNUCN + FKA0 + FHYPN + FETA + FPI0 C NEUTRAL PARTICLES FROM 3RD STRING FNEUT3 = RC3TO2 * FNEUT2 C TOTAL NUMBER OF NEUTRALS FNEUT = FNEUT2 + FNEUT3 NEUTOT = NINT( FNEUT ) C CALCULATE TOTAL NUMBER OF PARTICLES TO BE CREATED NTOTEM = NCH + NEUTOT IF ( DEBUG ) WRITE(MDEBUG,*) * ' FNUCN,FKA0,FHYPN,FETA,FPI0,FNEUT2,FNEUT3,NTOTEM=', * SNGL(FNUCN),SNGL(FKA0),SNGL(FHYPN),SNGL(FETA),SNGL(FPI0), * SNGL(FNEUT2),SNGL(FNEUT3),NTOTEM C LIMIT OF SECONDARIES PRODUCED (GIVEN BY SIZE OF ARRAY : 3000) C LIMIT IS ARRAY SIZE - SIZE OF LARGEST TARGET NUCLEUS(40) IF ( NTOTEM .GE. 2956 ) THEN WRITE(MONIOU,*) 'HDPM : REJECT EVENT WITH ',NTOTEM, * ' SECONDARIES' GOTO 1 ENDIF C SPECIAL TREATMENT IF MULTIPLICITY IS TOO LOW IF ( NTOTEM .LE. 3 ) ISEL = 1 C FRACTION OF THE VARIOUS NEUTRAL PARTICLES (NN, K(0), L+S0 AS PAIRS) C NORMALIZE WITH THE SUM OF ALL NEUTRAL PARTICLES FNORML = 1.D0 / ( 0.5D0 * (FNUCN+FKA0+FHYPN) + FETA + FPI0 ) RNUCNR = FNUCN * FNORML * 0.5D0 RKA0R = FKA0 * FNORML * 0.5D0 RHYPNR = FHYPN * FNORML * 0.5D0 RETAR = FETA * FNORML RPI0R = FPI0 * FNORML C CUMULATED RATIOS (NN, K(0), LAMBDA+SIGMA0 AS PAIRS) RPIER = RPI0R + RETAR RPEKR = RPIER + RKA0R RPEKNR = RPEKR + RNUCNR C THEN THE REMAINDER (1-RPEKNR) MUST BE NEUTRAL HYPERON PAIRS IF ( DEBUG ) WRITE(MDEBUG,*) * ' RPI0R,RETAR,RKA0R,RNUCNR,RHYPNR,FNORML=', * SNGL(RPI0R),SNGL(RETAR),SNGL(RKA0R),SNGL(RNUCNR),SNGL(RHYPNR), * SNGL(FNORML) C PROBABILITY TO PRODUCE CHARGED PIONS IS PROBABILITY NOT TO PRODUCE C CHARGED KAONS OR PROTONS OR CHARGED HYPERONS, WHERE PROTON/ANTIPROTON C IS HALF OF (ALL-NUCL)/(ALL-CHARGED) AUXIL = RKCH + 0.5D0 * RNUCCH + RHYPCH AUXIL3 = 1.D0 - AUXIL C RENORMALIZATION AS P/P_BAR, K+-, AND HYPERONS ARE PRODUCED IN PAIRS C AUXIL2 IS INVERSE OF NORMALISATION AUXIL2 = 1.D0 / (1.D0 - 0.5D0 * AUXIL) C CUMULATED PROBABILITIES (PP, K+-, SIGMA+- AS PAIRS) PPICH = AUXIL3 * AUXIL2 PPINCH = PPICH + 0.25D0 * RNUCCH * AUXIL2 PPNKCH = PPINCH + 0.5D0 * RKCH * AUXIL2 C THEN THE REMAINDER (1-PPNKCH) MUST BE CHARGED HYPERON PAIRS IF ( DEBUG ) WRITE(MDEBUG,*) ' PPICH,PPINCH,PPNKCH=', * SNGL(PPICH),SNGL(PPINCH),SNGL(PPNKCH) C NOW SELECT HOW MANY PARTICLES OF EACH TYPE ARE PRODUCED CALL PARNUM( INUMFL ) IF ( INUMFL .NE. 0 ) GOTO 1919 C DEFINE PARTICLE NUMBERS WHERE SPECIAL RAPIDITY IS CALCULATED C FOR PARTICLES FROM TARGET (THIRD STRING) PPP = RC3TO2 / (1.D0+RC3TO2) C NUMBER OF PARTICLES IN PROTON ANTIPROTON PAIRS FROM TARGET ITA = NINT(PPP * 2.D0 * NNC) C NUMBER OF PARTICLES IN K+ K- PAIRS FROM TARGET ITB = NINT(PPP * 2.D0 * NKC) C NUMBER OF PARTICLES IN SIGMA+ SIGMA- PAIRS FROM TARGET ITC = NINT(PPP * 2.D0 * NHC) C NUMBER OF PI+ PI- FROM TARGET ITD = NINT(PPP * NPC ) C CALCULATE BOUNDARIES IA1 = 2 IA2 = IA1 + ITA IB1 = IA1 + 2 * NNC IB2 = IB1 + ITB IC1 = IB1 + 2 * NKC IC2 = IC1 + ITC ID1 = IC1 + 2 * NHC ID2 = ID1 + ITD IE1 = ID1 + NPC C NUMBER OF PARTICLES IN NEUTRON ANTINEUTRON PAIRS FROM TARGET IE2 = IE1 + 2 * NNUCN(3) IF1 = IE1 + 2 * NNN C NUMBER OF PARTICLES IN K0S K0L PAIRS FROM TARGET IF2 = IF1 + 2 * NKA0(3) IG1 = IF1 + 2 * NKN C NUMBER OF PARTICLES IN NEUTRAL HYPERON PAIRS FROM TARGET IG2 = IG1 + 2 * NHYPN(3) IH1 = IG1 + 2 * NHN C NUMBER OF ETA FROM TARGET IH2 = IH1 + NETAS(3) II1 = IH1 + NET C NUMBER OF PI(0) FROM TARGET II2 = II1 + NPIZER(3) IJ1 = II1 + NPN IF ( DEBUG ) THEN WRITE(MDEBUG,*) ' CHARGED FROM TARGET:',ITA,ITB,ITC,ITD WRITE(MDEBUG,*) ' NEUTRAL FROM TARGET:', * 2*NNUCN(3),2*NKA0(3),2*NHYPN(3),NETAS(3),NPIZER(3) WRITE(MDEBUG,*) ' NTOTEM,IJ1=',NTOTEM,IJ1 ENDIF C REDEFINE TOTAL NUMBER OF SECONDARY PARTICLES : NTOTEM C BY CHARGE EXCHANGE AND RESONANCE FORMATION THIS NUMBER MAY BE ALTERED NTOTEM = IJ1 - 2 C----------------------------------------------------------------------- C RATIO OF RAPIDITY DENSITY TO MEAN PSEUDORAPIDITY IN CENTER C PARAMETRISATION SEE CAPDEVIELLE, J.PHYS.G:NUCL.PHYS.15(1989)909,EQ.6 IF ( XZ .LT. 1.5D0 ) THEN RDS = (0.24396D0 + 0.70150424D0 * XZ)**2 ELSE RDS = (0.55685D0 + 0.48664753D0 * XZ)**2 ENDIF C CALCULATE NOW: DN/DY AT Y = 0; DC0 IS AVERAGE PSEUDORAPIDITY DENSITY C TRAP IS RATIO (RAPID.DENS.)/(PSEUDORAP.DENS.) IN CENTER OF RAPIDITY TRAP = 1.25D0 IF ( IDIF .EQ. 0 .AND. ECMDPM .GT. 19.4D0 ) * TRAP = MAX( 1.D0, 1.28852D0 - 0.0065D0 * SMLOG ) DCN2 = DC0 * RDS * TRAP IF ( DEBUG ) WRITE(MDEBUG,*) ' RDS,TRAP,DCN2=', * SNGL(RDS),SNGL(TRAP),SNGL(DCN2) C AMPLITUDE OF GAUSSIAN 1ST AND 2ND STRING ATG2 = FNCH2 / (5.0132566D0 * WIDC2) C NEW DEFINITION OF POSITION BASED ON SEMI INCLUSIVE DATA SQ2 = 2.D0 * ATG2 / DCN2 C FINAL POSITION OF GAUSSIAN; WIDTH WIDC2 IS UNCHANGED IF ( SQ2 .GT. 1.D0 ) POSC2 = WIDC2 * SQRT( 2.D0*LOG(SQ2) ) C DENSITY OF CHARGED IN EXCESS FROM TARGET IN CENTER OF RAPIDITY DCN3 = 0.5D0 * (GNU - 1.D0) * DCN2 IF (DEBUG) WRITE(MDEBUG,*) ' SQ2,POSC2,DCN3=', * SNGL(SQ2),SNGL(POSC2),SNGL(DCN3) IF ( DCN3 .GT. 0.D0 ) THEN C AMPLITUDE 3RD GAUSSIAN ATG3 = FNCH3 / (5.0132566D0 * WIDC3) C AMPLITUDE IS DIVIDED BY DENSITY FOR GETTING CENTER OF 3RD GAUSSIAN SQ3 = 2.D0 * ATG3 / DCN3 C CHECK IF ADDITIVE MULTIPLICITY IS TOO LOW IF ( SQ3 .GT. 1.D0 ) POSC3 = WIDC3 * SQRT( 2.D0*LOG(SQ3) ) IF (DEBUG) WRITE(MDEBUG,*)' SQ3,POSC3=',SNGL(SQ3),SNGL(POSC3) ENDIF C NFLPI0 .EQ. 0 MEANS TREAT PI(0) RAPIDITY ACCORDING TO COLLIDER DATA IF ( NFLPI0 .EQ. 0 ) THEN C RATIO OF RAPIDITY DENSITY TO MEAN PSEUDORAPIDITY AT CENTER WITH Z<1.5 IF ( ZG .LT. 1.5D0 ) THEN RDG = (0.24396D0 + 0.70150424D0 * ZG)**2 ELSE RDG = (0.55685D0 + 0.48664753D0 * ZG)**2 ENDIF C GAMMAS USE RATIO TRAG TO CALCULATE RATIO OF RAPIDITY TO C PSEUDO RAPIDITY DENSITY IN CENTER (TRAG = 1.1 * 0.5 ). C FACTOR 0.5 COMES FROM RATIO NEUTRAL/CHARGED, AS WE USE DC0, WHICH C IS AVERAGE PSEUDORAPIDITY DENSITY FOR CHARGED PIONS TRAG = 0.55D0 IF ( IDIF .EQ. 0 ) THEN IF ( ECMDPM .GT. 19.4D0 ) * TRAG = MAX( 0.4D0, 0.6658D0 - 0.01954D0 * SMLOG ) IF ( ECMDPM .LE. 50.D0 ) THEN DCG = DC0 * RDG * TRAG ELSEIF ( ECMDPM .LE. 200.D0 ) THEN DCG = DC0 * RDG * TRAG * (1.D0 + 0.18D0 * LOG(ECMDPM/50.D0)) ELSE DCG = DC0 * RDG * TRAG * 1.25D0 ENDIF ELSE DCG = DC0 * RDG * TRAG ENDIF C DEFINE WIDTH OF STRINGS FOR NEUTRAL PIONS AND ETAS WIDN2 = WIDC2 * MIN( 1.D0, 1.12275D0 - 0.0208D0 * RSLOG ) C NEW DEFINITION OF CENTER OF GAUSSIAN BASED ON SEMI INCLUSIVE DATA C USING AMPLITUDE OF THE GAUSSIAN FOR NEUTRALS AUXIL = 2.D0 / (5.0132566D0 * WIDN2 * DCG) C TOTAL MULTIPLICITY USED FOR 1ST AND 2ND STRING OF PI(0) AND ETA C IS GIVEN BY THEIR NUMBERS. ANALOGOUS FOR 3RD STRING SP2 = DBLE ( NPIZER(2)+NETAS(2)) * AUXIL C FINAL CENTER OF GAUSSIANS FOR PI(0) AND ETA (WIDC2 IS UNCHANGED) IF ( SP2 .GT. 1.D0 ) THEN POSN2 = WIDN2 * SQRT( 2.D0 * LOG(SP2) ) ELSE POSN2 = POSC2 ENDIF WIDN3 = WIDN2 SP3 = DBLE(NPIZER(3)+NETAS(3)) * AUXIL IF ( SP3 .GT. 1.D0 ) THEN POSN3 = WIDN3 * SQRT( 2.D0 * LOG(SP3) ) ELSE POSN3 = POSC3 ENDIF ELSE C NFLPI0 .EQ. 1 MEANS RAPIDITY OF PI(0) AND ETA SAME AS THAT OF CHARGED POSN2 = POSC2 WIDN2 = WIDC2 POSN3 = POSC3 WIDN3 = WIDC3 ENDIF IF ( DEBUG ) WRITE(MDEBUG,*) * ' ZG,RDG,DCG,SP2,SP3,POSN2,POSN3,WIDN2 =', * SNGL(ZG),SNGL(RDG),SNGL(DCG),SNGL(SP2),SNGL(SP3),SNGL(POSN2), * SNGL(POSN3),SNGL(WIDN2) C----------------------------------------------------------------------- NREPR1 = 0 C RETURN POINT. NUMBERS OF PARTICLES REMAIN UNCHANGED FOR NEXT TRY, C BUT INDIVIDUAL RAPIDITIES GET NEW VALUES. C START FROM BEGINNING IF NO MATCH AFTER 20 TRIES 30 CONTINUE NREPR1 = NREPR1 + 1 IF ( NREPR1 .GT. 20 ) THEN IF ( IDIF .EQ. 1 .AND. NREPRD .LE. 10 ) GOTO 1919 GOTO 1 ENDIF C FOR TOTAL NUMBER OF PARTICLES ADD 2 FOR LEADER AND ANTILEADER NTOT = NTOTEM + 2 C PRODUCTION OF INDIVIDUAL RAPIDITIES FOR ALL SECONDARY PARTICLES CALL PARRAP CC IF ( DEBUG ) THEN CC WRITE (MDEBUG,*) ' RAPIDITIES:' CC WRITE (MDEBUG,134) (I,YR(I), I=3,NTOT) C134 FORMAT(' ',1P, (1X, I4, 5X, G13.6 )) CC ENDIF C CALCULATION OF CENTRAL RAPIDITY WITHOUT (ANTI)LEADER C MULTIPLICITY IN CENTER OF RAPIDITY DISTRIBUTION IZN = 0.D0 IF ( IDIF .EQ. 0 ) THEN DO 111 I = 3,NTOT IF ( ABS(YR(I)) .LT. DELRAP ) IZN = IZN + 1 111 CONTINUE IF ( IZN .LT. 1 ) THEN IF ( ISEL .EQ. 0 ) GOTO 30 C IN CASE OF FEW PARTICLES, SET PARTICLE NUMBER IN PLATEAU TO 1 IZN = 1 ENDIF C CENTRAL RAPIDITY DENSITY FOR CHARGED PARTICLES IF ( NTOTEM .GE. 1 ) THEN ZNC = MAX( 1.1D0, DBLE(NCH)*IZN/(DBLE(NTOTEM)*2.D0*DELRAP) ) ELSE ZNC = 1.1D0 ENDIF ELSE C DIFFRACTION: SHIFT RAPIDITIES + TAKE CENT.RAP.DENS. FROM PARAMETRISAT DO 112 I = 3,NTOT YR(I) = YR(I) + YY0 112 CONTINUE ZNC = MAX( 1.1D0, DCN2 ) ENDIF C ZN ACCOUNTS FOR THE RISE OF PT WITH CENTRAL RAP.DENSITY. THE FORMULA C IS A FIT TO UA1 VALUES OF ARNISON ET AL, PHYS.LETT.B118(1982)167 C REGARD, THAT OUR ZN IS DEFINED DIFFERENT FROM LITERATURE N BY 1 C - - - - - - C MODIFICATION AFTER J.N. CAPDEVIELLE, (DEC.96) IF ( ECMDPM .LE. 500.D0 ) THEN ZN = MAX( 1.00001D0, 3.669D0 / ZNC**0.435D0 + 6.4D0 ) ELSE C TAKE INTO ACCOUNT THE RESULTS OF UA1/MIMI EXPERIMENT IF ( ZNC .GE. 3.D0 ) THEN ZN = 1.D0 /(ZNC*0.004111D0 + 0.145D0)+3.D0 ELSE C FOR ROCH < 3.00 (MIMI) (TO BE USED IN PTRAM) ZM = 0.0033D0 * (ZNC-1.56D0)**2 + 0.406D0 ZN = 2.64D0/ZM + 3.D0 ENDIF ENDIF C - - - - - - C NOW SET PARTICLE TYPE AND TRANSV. MOMENTA FOR NEW PARTICLES IN PPARAM C SET ALSO TRANSVERSE MASS FOR ALL PARTICLES (INCL. LEADER+ANTILEADER) CALL PPARAM IF ( IDIF .EQ. 0 ) THEN C NOW SET THE RAPIDITY OF THE ANTILEADER ACCORDING TO THE DISTRIBUTION C IN THE FEYNMAN X VARIABLE; SET THE RAPIDITY OF LEADER TO CONSUME C THE REMAINDER OF ENERGY CALL LEDENY( LEDEFL ) IF ( LEDEFL .NE. 0 ) THEN IF ( DEBUG ) WRITE(MDEBUG,*) ' LEDEFL=',LEDEFL GOTO 30 ENDIF C CALCULATE FOR SINGLE COLLISION SYSTEM C.M. ENERGY + RAPIDITY SHIFT IF ( GNU .LE. 1.D0 ) THEN JGNU = 0.D0 DYGNU = 0.D0 ECMJAD = ECMDPM ELSE C MULTIPLE COLLISION IN TARGET JGNU = NINT(GNU-1.D0) C ADD ADDITIONALLY INTERACTING C TARGET NUCLEONS TO GET CORRECT JADACH FILTERING C CHOSE RANDOMLY WHETHER PROTON OR NEUTRON CALL RMMAR( RD,JGNU,1 ) IPR = 0 INE = 0 TARMAS = PAMA(ITYP(2)) DO 114 I = 1,JGNU NTOT = NTOT + 1 IF ( RD(I) .LE. .5D0 ) THEN ITYP(NTOT) = 13 INE = INE + 1 ELSE ITYP(NTOT) = 14 IPR = IPR + 1 ENDIF TMAS(NTOT) = PAMA(ITYP(NTOT)) TARMAS = TARMAS + TMAS(NTOT) EA(NTOT) = TMAS(NTOT) PX(NTOT) = 0.D0 PY(NTOT) = 0.D0 PT2(NTOT) = 0.D0 114 CONTINUE C CALCULATE C.M. ENERGY + RAPIDITY SHIFT * YCMGNU = 0.5D0 * LOG( (ELAB+TARMAS+PLAB)/(ELAB+TARMAS-PLAB) ) YCMGNU = 0.5D0 * LOG( (EPLUSP**2 +TARMAS*EPLUSP)/ * (PAMA(ITYPE)**2+TARMAS*EPLUSP) ) DYGNU = YCM - YCMGNU C CALCULATE RAPIDITIES OF ADDITIONALLY INTERACTING TARGET NUCLEONS C IN THE CM SYSTEM OF NUCLEON-NUCLEON SYSTEM DO 115 I = NTOT-JGNU+1,NTOT YR(I) = - YCM 115 CONTINUE C SHIFT RAPIDITIES INTO CM SYSTEM OF GNU+1 MASSES DO 113 I = 1,NTOT YR(I) = YR(I) + DYGNU 113 CONTINUE C CENTER OF MASS ENERGY OF 1 PROJECTILE AND GNU TARGET NUCLEONS TO C BE USED IN THE JADACH FILTER. ECMJAD = SQRT( PAMA(ITYPE)**2 + TARMAS**2 + 2.D0*TARMAS*ELAB ) ENDIF ELSE C IN CASE OF DIFFRACTION SET THE RAPIDITY OF LEADER AND ANTILEADER C IN SUBROUTINE LEADDF DYGNU = 0.D0 ECMJAD = ECMDPM CALL LEADDF( IFLGLD ) IF ( IFLGLD .NE. 0 ) THEN IF ( DEBUG ) WRITE(MDEBUG,*) ' IFLGLD=',IFLGLD GOTO 30 ENDIF ENDIF C CORRECT THE RAPIDITIES TO CONSERVE LONGITUDINAL MOMENTA AND ENERGY C USING THE ALGORITHM OF JADACH (SIMPLIFIED VERSION BY R. ATTALLAH) CALL JADACH( ECMJAD,JADFLG ) IF ( JADFLG .NE. 0 ) THEN IF ( DEBUG ) WRITE(MDEBUG,*) ' JADFLG=', JADFLG IF ( JADFLG .GT. 0 ) GOTO 30 IF ( JADFLG .LT. 0 ) THEN NREPRD = NREPRD + 1 IF ( NREPRD .GT. 10 ) GOTO 1 GOTO 1919 ENDIF ENDIF C CALCULATE LAB ENERGIES OF SECONDARY PARTICLES FROM THE RAPIDITIES C INCLUDING THE ADDITIONAL TARGET NUCLEONS ETOT = 0.D0 DO 135 I = 1,NTOT YR(I) = YR(I) + YCM - DYGNU EA(I) = TMAS(I) * COSH( YR(I) ) ETOT = ETOT + EA(I) 135 CONTINUE IF ( DEBUG ) WRITE(MDEBUG,136) * (I,ITYP(I),PX(I),PY(I),YR(I),EA(I),I=1,NTOT) 136 FORMAT(' NO ITYP PX PY YR EA'/ * (' ',I4,I3,1X,1P,4G13.6) ) C----------------------------------------------------------------------- C LOOP OVER ALL SECONDARY PARTICLES AND THE LEADING PARTICLE C PUT THEM ON THE STACK DO 139 LK = 5,8 SECPAR(LK) = CURPAR(LK) 139 CONTINUE C PROCESS LOOP DO 140 J = 1,NTOT C REJECTION OF BACKWARD GOING PARTICLES IF ( YR(J) .LE. 0.D0 ) THEN IF ( DEBUG ) WRITE(MDEBUG,*)'HDPM : YR REJECT PARTICLE ',J GOTO 140 ENDIF C CALCULATE THE PROPERTIES OF ALL SECONDARIES C PARTICLE TYPE SECPAR(1) = ITYP(J) C CALCULATE GAMMA FACTOR SECPAR(2) = EA(J) / PAMA(ITYP(J)) C TOTAL MOMENTUM SQUARED PTM = EA(J)**2 - PAMA(ITYP(J))**2 IF ( PT2(J) .GT. PTM ) THEN IF ( DEBUG ) WRITE(MDEBUG,*)'HDPM : PT REJECT PARTICLE ',J GOTO 140 ENDIF C EMISSION ZENITH ANGLE AGAINST TRAJECTORY OF PROJECTILE IF ( PTM .EQ. 0.D0 ) THEN COSTET = 1.D0 ELSE COSTET = SQRT( 1.D0 - PT2(J) / PTM ) ENDIF C EMISSION AZIMUTH ANGLE IF ( PX(J) .NE. 0.D0 .OR. PY(J) .NE. 0.D0 ) THEN PHIJ = ATAN2( PY(J), PX(J) ) ELSE PHIJ = 0.D0 ENDIF CALL ADDANG( COSTHE,PHI, COSTET,PHIJ, SECPAR(3),SECPAR(4) ) IF ( SECPAR(3) .LT. C(29) ) THEN C OMIT UPWARD GOING PARTICLES IF (DEBUG) WRITE(MDEBUG,*)'HDPM : ANGLE REJECT PARTICLE ',J GOTO 140 ENDIF C PUT SECONDARY PARTICLES ON STACK, IF NOT GOING UPWARDS IF ( J .GT. 2 ) THEN CALL TSTACK ELSE C PUT LEADER OR ANTI-LEADER ON STACK, IF NOT GOING UPWARDS IF ( ITYP(J) .GT. 50 ) THEN C LEADER OR ANTI LEADER ARE RESONANCES AND DECAY IRESPAR = IRESPAR + 1 IF ( IRESPAR .GE. 1000 ) THEN WRITE(MONIOU,*) 'STACK OF RESDEC RANDOM NUMBERS FULL' IRESPAR = 999 ENDIF RESRAN(IRESPAR) = RDRES(J) C COUNTER FOR ENERGY-MULTIPLICITY MATRIX MSMM = MSMM + 1 ENDIF CALL TSTACK C CALCULATE ELASTICITY FROM ENERGY OF LEADER (REST OF RESONANCE DECAY) IF ( J. EQ. 1 ) ELASTI = SECPAR(2)*PAMA(NINT(SECPAR(1)))/ELAB ENDIF C COUNTERS FOR FIRST INTERACTION IF ( FIRSTI ) THEN IF ( SECPAR(1) .EQ. 7.D0 .OR. SECPAR(1) .EQ. 8.D0 * .OR. SECPAR(1) .EQ. 9.D0 ) THEN IFINPI = IFINPI + 1 ELSEIF ( SECPAR(1) .EQ. 13.D0 .OR. SECPAR(1) .EQ. 14.D0 * .OR. SECPAR(1) .EQ. 15.D0 .OR. SECPAR(1) .EQ. 25.D0 ) THEN IFINNU = IFINNU + 1 ELSEIF ( SECPAR(1) .EQ. 10.D0 .OR. SECPAR(1) .EQ. 11.D0 * .OR. SECPAR(1) .EQ. 12.D0 .OR. SECPAR(1) .EQ. 16.D0 ) THEN IFINKA = IFINKA + 1 ELSEIF ( SECPAR(1) .GE. 71.D0 .AND. SECPAR(1) .LE. 74.D0) THEN IFINET = IFINET + 1 ELSEIF ((SECPAR(1) .GE. 18.D0 .AND. SECPAR(1) .LE. 24.D0) * .OR. (SECPAR(1) .GE. 26.D0 .AND. SECPAR(1) .LE. 32.D0))THEN IFINHY = IFINHY + 1 ENDIF ENDIF 140 CONTINUE C COUNTER FOR ENERGY-MULTIPLICITY MATRIX MSMM = MSMM + NTOT - 2 C FILL ELASTICITY IN MATRICES MEL = MIN ( 1.D0+10.D0* MAX( 0.D0, ELASTI ) , 11.D0 ) MEN = MIN ( 4.D0+ 3.D0*LOG10(MAX( .1D0, EKINL )), 37.D0 ) IELDPM(MEN,MEL) = IELDPM(MEN,MEL) + 1 IELDPA(MEN,MEL) = IELDPA(MEN,MEL) + 1 IF ( ELASTI .LT. 1.D0 ) THEN ELMEAN(MEN) = ELMEAN(MEN) + ELASTI ELMEAA(MEN) = ELMEAA(MEN) + ELASTI ENDIF IF ( FIRSTI ) THEN ELAST = ELASTI FIRSTI = .FALSE. ENDIF IF ( DEBUG ) WRITE(MDEBUG,*) 'HDPM : ELAST=',SNGL(ELASTI), * SNGL(ETOT),SNGL(ELAB) RETURN END