PROGRAM DIREZIO C C-O---------------------------------------------------------------------O C-O -O C-O PURPOSE AND METHODS : RICOSTRUIRE EVENTO PER EVENTO LA -O C-O DIREZIONE D'ARRIVO DELLA PARTICELLA (PUNTAMENTO) -O C-O -O C-O INPUTS : DATA FILE, OUTPUT FILE, CONTAINEMENT OR NOT -O C-O OUTPUTS : RZ FILE FOR PAW -O C-O CONTROLS: -O C-O -O C-O -O C-O CREATED FEB-2000 BY EMILIANO MOCCHIUTTI -O C-O -O C-O---------------------------------------------------------------------O C C PARAMETER (NUMEV=100000, PARAMETER (NVAR=1, TR=5.) !5. PARAMETER (NPLA=22,NCHA=96,LENSEV=NPLA*NCHA) PARAMETER (SOG=60.) PARAMETER (NWPAWC=3000000) PARAMETER (NWORD=100000) PARAMETER (NPAR=2, NPAR2=2, NPL=44) PARAMETER (PI=3.14159265358979324) PARAMETER (CALIB=0.0001059994) C INTEGER LREC PARAMETER (LREC=6144) C INTEGER EVCON, EVCOINC, EVCONMIS INTEGER TRIG INTEGER POSTRIX, POSTRIY INTEGER NPU, FLA, FLA2 INTEGER NDEP INTEGER EVENTO INTEGER NLK, NLK1, NLK2 INTEGER F INTEGER FINPAR1, FINPAR INTEGER IFIN2, NPLL INTEGER NSTR INTEGER NINF, NSUP INTEGER NLL INTEGER NPLOTX, NPLOTY, NPLOTX2, NPLOTY2, NNPLO C REAL TRASX, TRASY REAL POSIX, POSIY, POSIXM, POSIYM REAL QQX, QQY, QQZ, QQXM, QQYM, QQZM REAL VALORI(NWORD) REAL ZXY, PX, PY, THET, PH REAL X(NPLA), X2(NPLA), Y, Y2, EY, RN REAL XV(NPLA), YV, EYV(NPLA), PMI, EXV(NPLA) REAL EX REAL YN REAL ZZ REAL YM, EYM REAL THX, THY, THXM ,THYM REAL PPY, PPX REAL VENET, ENET, ERENET REAL ER1, ER2, ERENETM, ENETM REAL ENM, ERENM REAL TMIS, PMIS REAL ZETA, ZETAM, ZETA2, ZETA3 REAL PARZEN, PARZ, PARZE, ERZEN, ER3, ER4 REAL LIMLA, LIMLA2, LIMLA3, LIMLA4, LIMLA5 REAL FFLA REAL PMETX, PMETY REAL ENMEZX, ENMEZY REAL DELTX, DELTY REAL PO3S, FFLA2 REAL XEN(NPL), YEN(NPL), EYEN(NPL) REAL XEN2(NPL), YEN2(NPL), EYEN2(NPL) REAL POSMAX, PAR1 REAL QPLX, QPLY REAL X0ATT, FFLACO3, TOTTO2 REAL SCARTA, BANDI REAL VERT, ENCEN REAL ESTMA REAL NUMEV C DOUBLE PRECISION THETA, CPHI, SPHI, PHI DOUBLE PRECISION POSIXMD, POSIYMD, PMISD, TMISD, PMID DOUBLE PRECISION POSIXD, POSIYD DOUBLE PRECISION MX, QX, MY, QY, QXP, QYP DOUBLE PRECISION THETA2, POSIX2, POSIY2 DOUBLE PRECISION MMX, MMY, TTHX, TTHY C CHARACTER*6 CHOPT CHARACTER*20 FILENAME,NAME CHARACTER*80 CHTITLE CHARACTER*1 ALFA, BETA C EXTERNAL FUN EXTERNAL ENERGIF C COMMON/ENSTRIP/ZXY(2,LENSEV), YN(NCHA) C COMMON/VETFIT/YV(NPLA) C COMMON/DATI/THET, PH, TMIS, PMIS, THX, THY, THXM, THYM, & FLA2, FLA, QQX, QQY, QQZ, QQXM, QQYM, QQZM, & ZZ(2,NPLA), TRASX(NPLA), TRASY(NPLA), EVENTO, & Y(NPLA), EX(NPLA), Y2(NPLA), EY(NPLA), & VENET, PARZ(2,NPLA), PARZEN, FFLA, & ENMEZX(2,NPLA),ENMEZY(2,NPLA), & DELTX(NPLA), DELTY(NPLA), & FFLA2, POSMAX, NSTR(4,NPLA), PAR1, & QPLX(6,NPLA),QPLY(6,NPLA),VERT,ENCEN(12,NPLA), & NNPLO(2,NPLA),ESTMA(3) C C COMMON/DATI2/TOTTO2,SCARTA C COMMON/HCFITD/PA(NPAR),SIGPAR(NPAR),CHI2,PMIN(NPAR),PMAX(NPAR), & STEP(NPAR), & PAR(NPAR2),SIGPA(NPAR2),CHI,PMI2(NPAR2),PMA(NPAR2), & STE(NPAR2) COMMON/QUEST/IQUEST(100) COMMON/PAWC/HMEMOR(NWPAWC) C C ################################################################# C CALL HLIMIT(NWPAWC) CALL HDELET(0) C PRINT *,'FILENAME TO SAVE?' READ (*,5050)NAME PRINT *,'NUMERO DI EVENTI?' READ (*,6006),NUMEV IQUEST(10)=300000 CALL HROPEN(13,'SIMULATION', NAME,'N',LREC,ISTAT) IF (ISTAT.NE.0) GO TO 1000 C CALL HBNT(1,'DATI',' ') CALL HBSET('BSIZE',LREC,IERR) CALL HBNAME(1,'DATI',THET, & 'THET:R,PH:R,TMIS:R,PMIS:R,THX:R,THY:R,THXM:R,THYM:R,FLA2:I, & FLA:I,QQX:R,QQY:R,QQZ:R,QQXM:R,QQYM:R,QQZM:R,ZZ(2,22):R, & TRASX(22):R,TRASY(22):R,EVENTO:I,X(22):R,EX(22):R,Y(22):R, & EY(22):R,VENET:R,PARZ(2,22):R,PARZEN:R,FFLA:R,ENMEZX(2,22):R, & ENMEZY(2,22):R,DELTX(22):R,DELTY(22):R,FFLA2:R,POSMAX:R, & NSTR(4,22):I,PAR1:R,QPLX(6,22):R,QPLY(6,22):R,VERT:R, & ENCEN(12,22):R,NNPLO(2,22):I,ESTMA(3):R') C C CALL HBNT(2,'DATI2',' ') C CALL HBSET('BSIZE',LREC,IERR) C CALL HBNAME(2,'DATI2',TOTTO2,'TOTTO2:R,SCARTA:R') C PRINT *,'CONTAINED OR NOT ? (C/N)' C$ READ (*,6001),ALFA ALFA = 'C' C PRINT *,'PROFILO LONGITUDINALE? (Y/N)' C$ READ (*,6001),BETA BETA = 'Y' C IF (ALFA.EQ.'C'.OR.ALFA.EQ.'c') THEN PRINT *,'DA CHE PIANO VUOI CHE ENTRINO LE PARTICELLE?' C$ READ (*,6005),LIMLA2 LIMLA2 = 0. PRINT *,'PER QUANTI PIANI?' C$ READ (*,6005),LIMLA5 LIMLA5 = 8. PRINT *,'DA CHE PIANO VUOI CHE ESCANO?' C$ READ (*,6005),LIMLA4 LIMLA4 = 12. PRINT *,'LIMITAZIONE SULLA PROIEZIONE SUL PIANO 22: ' C$ READ (*,6006),LIMLA3 LIMLA3 = 999. LIMLA4 = 22. - LIMLA4 LIMLA = 12.1+LIMLA4/TAN((18.964-ABS(LIMLA2+LIMLA4)*.862)/24.2) PRINT *,'LIMLA CALCOLATO = ',LIMLA IF (LIMLA.GT.LIMLA3) LIMLA = LIMLA3 PRINT *,'LIMLA = ',LIMLA ELSE LIMLA = 10000. LIMLA2 = 0. LIMLA5 = 11. LIMLA4 = 11. LIMLA3 = 999. ENDIF C TRIG = 0 EVCONMIS = 0 EVCON = 0 EVCOINC = 0 EVENTO = 0 ENETM = 0. ERENETM =0. ER1 = 0. ER2 = 0. ER3 = 0. ER4 = 0. TOTTO2 = 0. SCARTA = 0. C CALL VZERO(YN,NCHA) CALL VZERO(X,NPLA) CALL VZERO(X2,NPLA) C C CONVERSIONE STRIP-CM C DO NN = 1, NCHA YN(NN) = (NN - 1.) * .244 +.218 !.2448 + .185 IF (NN.GT.32) YN(NN) = YN(NN) + .292 IF (NN.GT.64) YN(NN) = YN(NN) + .292 ENDDO C C CONVERSIONE XVIEW-CM, YVIEW-CM C DO NN = 1, NPLA RN = FLOAT(NN) X(NN) = (RN - 1.)*.862 + .169 X2(NN) = (RN - 1.)*.862 + .843 ENDDO C PRINT *,'DATA FILENAME : ' READ (*,5050)FILENAME C OPEN (UNIT=83,FILE=FILENAME,STATUS='OLD',FORM='UNFORMATTED') C DO JJ=1, NUMEV C READ (83,END=900),LENGHT CALL VZERO(ZXY,2*LENSEV) CALL VZERO(VALORI,NWORD) CALL VZERO(TRASX,22) CALL VZERO(TRASY,22) NEVENT = NEVENT + 1 CALL LEGGE(VALORI,LENGHT) C C VALORI IS A VECTOR OF VARIABLE LENGHT (LENGHT) WHICH CONTAINS ALL THE NEEDED C INFORMATION FROM THE SIMULATION C IF (LENGHT.LT.17) GO TO 449 C C PUT IN THE VECTOR ZXY ALL THE ENERGIES RELEASED BY THE C PARTICLE IN EVERY STRIP C DO I = 6,LENGHT-17,2 KK = VALORI(I) IF (KK.LT.5000) THEN ZXY(1,KK) = VALORI(I+1) C C SATURAZIONE DELLE STRIP A 1100MIP C IF (ZXY(1,KK).GT.1100.) ZXY(1,KK) = 1100. ELSE KS = KK - 5000 ZXY(2,KS) = VALORI(I+1) IF (ZXY(2,KS).GT.1100.) ZXY(2,KS) = 1100. ENDIF ENDDO EVENTO = EVENTO + 1 C C RICAVA GLI ANGOLI VERI DI INCIDENZA... C ZETA = ABS(VALORI(LENGHT)) THETA = DBLE(VALORI(LENGHT-10)) IF (THETA.EQ.0.) THEN CPHI = 0. SPHI = 1. ELSE CPHI = DBLE(VALORI(3))/DSIND(THETA) SPHI = DBLE(VALORI(4))/DSIND(THETA) ENDIF IF (CPHI.GT.1.) CPHI = 1. IF (CPHI.LT.-1.) CPHI = -1. IF (SPHI.GE.0.) PHI = DACOSD(CPHI) IF (SPHI.LT.0.) PHI =-DACOSD(CPHI)+DBLE(360.) C C SE LA PARTICELLA ARRIVA DALLA FACCIA X LATERALE GLI ANGOLI CHE CI DA` C LA SIMULAZIONE SONO RELATIVI AD UN SIS DI RIF RUOTATO E LI SI DEVONO C CONVERTIRE NEGLI ANGOLI GIUSTI C IF (ZETA.NE.18.964) THEN IF (PHI.LT.90.) THEN TTHY = PHI ELSE TTHY = PHI - DBLE(360.) ENDIF IF (PHI.LT.90.OR.PHI.GT.270.) THEN TTHX = DATAND(DABS(DTAND(THETA))/(DSQRT( & DBLE(1.)+DTAND(PHI)**2.)))-DBLE(90.) ENDIF IF (PHI.GT.90..AND.PHI.LT.270.) THEN TTHX = DATAND(-DABS(DTAND(THETA))/(DSQRT( & DBLE(1.)+DTAND(PHI)**2.)))-DBLE(90.) ENDIF MMX = DTAND(TTHX) MMY = DTAND(TTHY) THETA = DATAND(DSQRT((MMX*MMX) + (MMY*MMY))) C IF (MMX.EQ.0..AND.MMY.GT.0.) PHI = DBLE(90.) IF (MMX.EQ.0..AND.MMY.LT.0.) PHI = DBLE(270.) IF (MMY.EQ.0..AND.MMX.GE.0.) PHI = DBLE(0.) IF (MMY.EQ.0..AND.MMX.LT.0.) PHI = DBLE(180.) C IF (MMY.NE.0..AND.MMX.NE.0.) THEN PHI = DATAND(MMY/MMX) IF (MMY.LT.0..AND.MMX.GT.0.) PHI = PHI + DBLE(360.) IF (MMX.LT.0.) PHI = PHI + DBLE(180.) ENDIF CPHI = DCOSD(PHI) SPHI = DSIND(PHI) ENDIF C C ...E LA PROIEZIONE DELLA TRAIETTORIA SULL'ULTIMO PIANO C POSIXD=DBLE(VALORI(LENGHT-13))+DBLE(ZETA)* & DTAND(THETA)*CPHI POSIYD=DBLE(VALORI(LENGHT-12))+DBLE(ZETA)* & DTAND(THETA)*SPHI C C CONTAINED OR NOT? C IF (ALFA.EQ.'C'.OR.ALFA.EQ.'c') THEN IF (POSIXD.GT.LIMLA.OR.POSIXD.LT.-LIMLA.OR.POSIYD.GT. & LIMLA.OR.POSIYD.LT.-LIMLA) GO TO 449 ZETA2 = 18.964 - LIMLA2*.862 ZETA3 = 18.964 - (LIMLA2+LIMLA5)*.862 IF (ZETA.GT.ZETA2.OR.ZETA.LT.ZETA3) GO TO 449 ENDIF C C EVENTI REALMENTE CONTENUTI C IF (POSIXD.LE.LIMLA.AND.POSIXD.GE.-LIMLA.AND.POSIYD.LE. & LIMLA.AND.POSIYD.GE.-LIMLA) EVCON = EVCON + 1 C### c BANDI = 0. c QQZ = ZETA c PX = SNGL(POSIXD) c PY = SNGL(POSIYD) c THET = SNGL(THETA) c PH = SNGL(PHI) c QQX = SNGL(VALORI(LENGHT-13)) c IF ((QQZ.GE.12.964.AND.QQZ.LE.18.964.AND. c & (ABS(PX).GE.0.AND.ABS(PX).LT.(12.1+1.35*4)).AND. c & (ABS(PY).GE.0.AND.ABS(PY).LT.(12.1+1.35*4)).AND. c & (THET.GE.0.AND.THET.LT.63.5))) THEN cC c FFLACO3=(ABS(QQZ)+ABS((12.1-QQX)-18.964)/(TAND(THET)*COSD(PH))) cC c X0ATT = .76*(FFLACO3/(COSD(THET))) cC c IF (X0ATT.GT.10.) THEN c TOTTO2 = TOTTO2 + 1. c BANDI = 1. c ENDIF cC c ENDIF ccc IF (BANDI.NE.1.) GO TO 449 C### C TROVA IL PIANO SU CUI E` STATA RILASCIATA LA MASSIMA ENERGIA C (FINPAR) E REGISTRA LE ENERGIE RILASCIATE SU CIASCUN PIANO (PARZ) C ENET = 0. PARZE = 0. FINPAR = 1 FINPAR1 = 1 CALL VZERO(PARZ,2*NPLA) CALL VZERO(ESTMA,3) DO NN = 1, NPLA DO KK = 1, NCHA NLK = (KK - 1)*NPLA + NN IF (ZXY(1,NLK).GE.ESTMA(1).AND.ZXY(1,NLK).LT.1100.) THEN ESTMA(1) = ZXY(1,NLK) ESTMA(2) = NN ESTMA(3) = KK ENDIF IF (ZXY(2,NLK).GE.ESTMA(1).AND.ZXY(2,NLK).LT.1100.) THEN ESTMA(1) = ZXY(2,NLK) ESTMA(2) = NN ESTMA(3) = KK ENDIF PARZ(1,NN) = PARZ(1,NN) + ZXY(1,NLK) PARZ(2,NN) = PARZ(2,NN) + ZXY(2,NLK) ENET = ENET + ZXY(1,NLK) + ZXY(2,NLK) ENDDO IF (PARZ(1,NN).GE.PARZE) THEN PARZE = PARZ(1,NN) FINPAR = NN FINPAR1 = 1 ENDIF IF (PARZ(2,NN).GE.PARZE) THEN PARZE = PARZ(2,NN) FINPAR = NN FINPAR1 = 2 ENDIF ENDDO cc FINPAR = FINPAR + 3 cc IF (FINPAR.GT.NPLA) FINPAR = NPLA C$$$$ C C CERCA IL CLUSTER DI TRE STRIP CHE HA VISTO LA MAGGIORE ENERGIA C PER CIASCUN PIANO C PO3S = 0. DO I = 1, NPLA STRIMAX = 0. STRIMAY = 0. DO M = 1, 94 NLK = (M - 1) * NPLA + I NLK1 = M * NPLA + I NLK2 = (M + 1) * NPLA + I IF (ZXY(1,NLK).GT.TR.AND.ZXY(1,NLK1).GT.TR & .AND.ZXY(1,NLK2).GT.TR) THEN PPX = ZXY(1,NLK) + ZXY(1,NLK1) + ZXY(1,NLK2) ELSE PPX = 0. ENDIF IF (ZXY(2,NLK).GT.TR.AND.ZXY(2,NLK1).GT.TR & .AND.ZXY(2,NLK2).GT.TR) THEN PPY = ZXY(2,NLK) + ZXY(2,NLK1) + ZXY(2,NLK2) ELSE PPY = 0. ENDIF IF (PPX.GT.STRIMAX) THEN STRIMAX = PPX POSTRIX = M + 1 ENDIF IF (PPY.GT.STRIMAY) THEN STRIMAY = PPY POSTRIY = M + 1 ENDIF ENDDO C C REGISTRA LA POSIZIONE DELLA STRIP CENTRALE C IF (STRIMAX.GT.0.) THEN TRASX(I) = REAL(POSTRIX) + .122 ELSE TRASX(I) = -10. ENDIF IF (STRIMAY.GT.0.) THEN TRASY(I) = REAL(POSTRIY) + .122 ELSE TRASY(I) = -10. ENDIF C IF (STRIMAX.GT.PO3S) THEN PO3S = STRIMAX FFLA2 = FLOAT(I) ENDIF C ENDDO C C LA PRIMA VOLTA IL FIT VIENE ESEGUITO SULLE STRIP CENTRALI C DO F = 1, 2 cc F = 1 C C PREPARA I VETTORI PER IL SECONDO FIT; ADESSO LA STRIP "CENTRALE" C DIVENTA QUELLA DOVE PASSA LA TRAIETTORIA DELLA PRIMA RETTA C IF (F.EQ.2) THEN DO NN = 1, NPLA PY = X(NN)*SNGL(MX) + SNGL(QX) CALL TRAIE(PY,TRASX(NN)) C PY = X2(NN)*SNGL(MY) + SNGL(QY) CALL TRAIE(PY,TRASY(NN)) C C CORR. C IF (NN.GT.(FINPAR+1)) THEN IF (MX.LT.-0.58) & TRASX(NN) = TRASX(NN) - .4 IF (MX.GT.0.58) & TRASX(NN) = TRASX(NN) + .4 IF (MY.LT.-0.58) & TRASY(NN) = TRASY(NN) - .4 IF (MY.GT.0.58) & TRASY(NN) = TRASY(NN) + .4 ENDIF IF (NN.LT.(FINPAR-1)) THEN IF (MX.LT.-0.58) & TRASX(NN) = TRASX(NN) + .2 IF (MX.GT.0.58) & TRASX(NN) = TRASX(NN) - .2 IF (MY.LT.-0.58) & TRASY(NN) = TRASY(NN) + .2 IF (MY.GT.0.58) & TRASY(NN) = TRASY(NN) - .2 ENDIF C ENDDO ENDIF C C NUMERO DI STRIP SOPRA LA SOGLIA SOG ATTORNO (48 A SX E 48 A DX) C ALLA POSIZIONE DELLA TRAIETTORIA C IF (F.EQ.2) THEN CALL VZERO(NSTR,4*NPLA) CALL VZERO(ENCEN,12*NPLA) CALL VZERO(NNPLO,2*NPLA) VERT = -1. DO NN = 1, NPLA NPLOTX = 0 NPLOTY = 0 NPLOTX2 = 0 NPLOTY2 = 0 DO NJ = 1, NCHA NLK = (NJ - 1.)*NPLA + NN IF (ZXY(1,NLK).GT.0.) NSTR(3,NN) = NSTR(3,NN)+1 IF (ZXY(2,NLK).GT.0.) NSTR(4,NN) = NSTR(4,NN)+1 IF (ZXY(1,NLK).GT.0.) THEN NPLOTX = NPLOTX + 1 ELSE NPLOTX = 0 ENDIF IF (ZXY(2,NLK).GT.0.) THEN NPLOTY = NPLOTY + 1 ELSE NPLOTY = 0 ENDIF IF (NPLOTX.GT.NNPLO(1,NN)) NNPLO(1,NN) = NPLOTX IF (NPLOTY.GT.NNPLO(2,NN)) NNPLO(2,NN) = NPLOTY IF (NPLOTX.GE.2) NPLOTX2 = 1 IF (NPLOTY.GE.3) NPLOTY2 = 1 ENDDO IF (VERT.LT.0..AND.NPLOTX2.EQ.1.AND.NPLOTY2.EQ.1) & VERT = FLOAT(NN) IF (NINT(TRASX(NN)).GT.0.) THEN NINF = NINT(TRASX(NN)) - 48 NSUP = NINT(TRASX(NN)) + 48 IF (NINF.LT.1) NINF = 1 IF (NSUP.GT.96) NSUP = 96 DO NJ = NINF, NSUP NLK = (NJ - 1.)*NPLA + NN IF (ZXY(1,NLK).GT.SOG.AND.ZXY(1,NLK).GT.0.) & NSTR(1,NN) = NSTR(1,NN) + 1 C C 1 RM C IF (NJ.GE.NINT(TRASX(NN)-3.).AND. & NJ.LE.NINT(TRASX(NN)+3.)) THEN ENCEN(1,NN) = ENCEN(1,NN)+ZXY(1,NLK) IF (ZXY(1,NLK).GT.0.) THEN ENCEN(3,NN) = ENCEN(3,NN)+1. ENDIF ENDIF C C 3 RM C IF (NJ.GE.NINT(TRASX(NN)-12.).AND. & NJ.LE.NINT(TRASX(NN)+12.)) THEN ENCEN(5,NN) = ENCEN(5,NN)+ZXY(1,NLK) IF (ZXY(1,NLK).GT.0.) THEN ENCEN(7,NN) = ENCEN(7,NN)+1. ENDIF ENDIF C C 6 RM C IF (NJ.GE.NINT(TRASX(NN)-25.).AND. & NJ.LE.NINT(TRASX(NN)+25.)) THEN ENCEN(9,NN) = ENCEN(9,NN)+ZXY(1,NLK) IF (ZXY(1,NLK).GT.0.) THEN ENCEN(11,NN) = ENCEN(11,NN)+1. ENDIF ENDIF ENDDO ELSE NSTR(1,NN) = -1 C ENCEN(1,NN) = -1. DO NJ = 1, 12, 2 ENCEN(NJ,NN) = 0. ENDDO ENDIF IF (NINT(TRASY(NN)).GT.0.) THEN NINF = NINT(TRASY(NN)) - 48 NSUP = NINT(TRASY(NN)) + 48 IF (NINF.LT.1) NINF = 1 IF (NSUP.GT.96) NSUP = 96 DO NJ = NINF, NSUP NLK = (NJ - 1.)*NPLA + NN IF (ZXY(2,NLK).GT.SOG.AND.ZXY(2,NLK).GT.0.) & NSTR(2,NN) = NSTR(2,NN) + 1 C C 1 RM C IF (NJ.GE.NINT(TRASY(NN)-3.).AND. & NJ.LE.NINT(TRASY(NN)+3.)) THEN ENCEN(2,NN) = ENCEN(2,NN)+ZXY(2,NLK) IF (ZXY(2,NLK).GT.0.) THEN ENCEN(4,NN) = ENCEN(4,NN)+1. ENDIF ENDIF C C 3 RM C IF (NJ.GE.NINT(TRASY(NN)-12.).AND. & NJ.LE.NINT(TRASY(NN)+12.)) THEN ENCEN(6,NN) = ENCEN(6,NN)+ZXY(2,NLK) IF (ZXY(2,NLK).GT.0.) THEN ENCEN(8,NN) = ENCEN(8,NN)+1. ENDIF ENDIF C C 6 RM C IF (NJ.GE.NINT(TRASY(NN)-25.).AND. & NJ.LE.NINT(TRASY(NN)+25.)) THEN ENCEN(10,NN) = ENCEN(10,NN)+ZXY(2,NLK) IF (ZXY(2,NLK).GT.0.) THEN ENCEN(12,NN) = ENCEN(12,NN)+1. ENDIF ENDIF ENDDO ELSE NSTR(2,NN) = -1 C ENCEN(2,NN) = -1. DO NJ = 1, 12, 2 ENCEN(NJ+1,NN) = 0. ENDDO ENDIF ENDDO ENDIF C C CALCOLA IL BARICENTRO, PER CIASCUN PIANO, RISPETTO 7 STRIP C NEL PRIMO FIT, RISPETTO 3 STRIP NEL SECONDO C DO NN = 1, NPLA CALL BARIC(TRASX(NN),Y(NN),EX(NN),NN,F,1) C CALL BARIC(TRASY(NN),Y2(NN),EY(NN),NN,F,2) ENDDO C C --- PER LE VISTE X --- C METTE NEI VETTORI YV XV EXV I DATI DA FITTARE (ESCLUDE GLI ZERI) C CALL VZERO(XV,NPLA) CALL VZERO(YV,NPLA) CALL VZERO(EXV,NPLA) CALL VZERO(DELTX,NPLA) FLA2 = 0 DO KK = 1, NPLA !!finpar IF (Y(KK).GT.0.) THEN FLA2 = FLA2 + 1 XV(FLA2) = X(KK) YV(FLA2) = Y(KK) EXV(FLA2) = EX(KK) PY = X(KK)*TAN(THX) + QQX + 12.1 DELTX(KK) = PY - YV(FLA2) ENDIF ENDDO IF (FLA2.LT.2) THEN c IF (BANDI.EQ.1.) SCARTA = SCARTA + 1. GO TO 449 ENDIF C C RICHIAMA SOLO NEL PRIMO PASSAGGIO UN ANALISI TOPOLOGICA C DELL'EVENTO PER ESCLUDERE I PUNTI CHE CORRONO SUL BORDO C IF (F.EQ.1) THEN CALL BORDO(FLA2) IF (FLA2.LT.2) THEN c IF (BANDI.EQ.1.) SCARTA = SCARTA + 1. GO TO 449 ENDIF ENDIF C C ESEGUE IL FIT DELLA RETTA, LA PRIMA VOLTA A PARTIRE DA UNA RETTA C A PENDENZA 0, LA SECONDA DALLA RETTA DEL PRIMO FIT C CHOPT = 'QE' IF (F.EQ.1) THEN PA(1) = DBLE(12.1) PA(2) = DBLE(0.) ELSE PA(1) = QX PA(2) = MX ENDIF C CALL HFITV(FLA2,NPLA,NVAR,XV,YV,EXV,FUN,CHOPT,NPAR,PA,STEP, & PMIN,PMAX,SIGPAR,CHI2) C C MX COEFF. ANGOLARE, QX QUOTA C MX = PA(2) QX = PA(1) C C REGISTRA I PULL C IF (F.EQ.2) THEN DO I = 1, NPLA ZZ(1,I) = -100. ZZ(2,I) = -100. ENDDO C DO I = 1, FLA2 NDEP = NINT(1.+((XV(I)-.169)/.862)) YVM = XV(I)*SNGL(MX) + SNGL(QX) EYM = SQRT((XV(I)*SNGL(SIGPAR(2)))**2+(SNGL(SIGPAR(1)))**2) ZZ(1,NDEP) = (YV(I)-YVM)/SQRT(ABS((EXV(I))**2-EYM**2)) ENDDO ENDIF C C --- PER LE VISTE Y --- C METTE NEI VETTORI XV YV EYV I DATI DA FITTARE (ESCLUDE GLI ZERI) C CALL VZERO(XV,NPLA) CALL VZERO(YV,NPLA) CALL VZERO(EYV,NPLA) CALL VZERO(DELTY,NPLA) FLA = 0 DO KK = 1, NPLA !!finpar IF (Y2(KK).GT.0.) THEN FLA = FLA + 1 XV(FLA) = X2(KK) YV(FLA) = Y2(KK) EYV(FLA) = EY(KK) PY = X2(KK)*TAN(THY) + QQY + 12.1 DELTY(KK) = PY - YV(FLA) ENDIF ENDDO IF (FLA.LT.2) THEN c IF (BANDI.EQ.1.) SCARTA = SCARTA + 1. GO TO 449 ENDIF C C COME PER LE VISTE X C IF (F.EQ.1) THEN CALL BORDO(FLA) IF (FLA.LT.2) THEN c IF (BANDI.EQ.1.) SCARTA = SCARTA + 1. GO TO 449 ENDIF ENDIF C C ESEGUE IL FIT DELLA RETTA C CHOPT = 'QE' IF (F.EQ.1) THEN PA(1) = DBLE(12.1) PA(2) = DBLE(0.) ELSE PA(1) = QY PA(2) = MY ENDIF C CALL HFITV(FLA,NPLA,NVAR,XV,YV,EYV,FUN,CHOPT,NPAR,PA,STEP, & PMIN,PMAX,SIGPAR,CHI2) C C MY COEFF. ANGOLARE, QY QUOTA C MY = PA(2) QY = PA(1) C C REGISTRA I PULL C IF (F.EQ.2) THEN DO I = 1, FLA2 NDEP = NINT(1.+((XV(I)-.843)/.862)) YVM = XV(I)*SNGL(MY) + SNGL(QY) EYM = SQRT((XV(I)*SNGL(SIGPAR(2)))**2+(SNGL(SIGPAR(1)))**2) ZZ(2,NDEP) = (YV(I)-YVM)/SQRT(ABS((EYV(I))**2-EYM**2)) ENDDO ENDIF C C RIPETE UNA SECONDA VOLTA IL FIT C ENDDO C C CORREZIONE DOVUTA ALLA NON SIMMETRIA DELLO SCIAME C cc QX = QX - DBLE(.1) * DATAN(MX) !.038818 cc QY = QY - DBLE(.1) * DATAN(MY) !.243585 cc MX = DTAN(DATAN(MX)/DBLE(1.-.0064651)) !.00506855)) cc MY = DTAN(DATAN(MY)/DBLE(1.-.0064651)) !.00506855)) C C CALCOLA THETA (analoga e` DACOSD(1./SQRT(1.+(MX**2)+(MY**2))) ) C TMISD = DATAND(DSQRT((MX*MX) + (MY*MY))) C C CALCOLA PHI C IF (MX.EQ.0..AND.MY.GT.0.) PMISD = DBLE(90.) IF (MX.EQ.0..AND.MY.LT.0.) PMISD = DBLE(270.) IF (MY.EQ.0..AND.MX.GE.0.) PMISD = DBLE(0.) IF (MY.EQ.0..AND.MX.LT.0.) PMISD = DBLE(180.) C IF (MY.NE.0..AND.MX.NE.0.) THEN PMID = DATAND(MY/MX) IF (MY.LT.0..AND.MX.GT.0.) PMISD = PMID + DBLE(360.) IF (MX.LT.0.) PMISD = PMID + DBLE(180.) IF (MY.GT.0..AND.MX.GT.0.) PMISD = PMID ENDIF C C TROVA IL VERTICE DI INTERAZIONE SECONDO IL FIT C QXP = QX - DBLE(12.1) QYP = QY - DBLE(12.1) C ZETAM = 18.964 C IF (QXP.GT.12.1) THEN ZETAM = SNGL(DBLE(18.964)+(QXP-12.1)/MX) QXP = DBLE(12.1) ENDIF IF (QXP.LT.-12.1) THEN ZETAM = SNGL(DBLE(18.964)+(QXP+12.1)/MX) QXP = DBLE(-12.1) ENDIF IF (QYP.GT.12.1) THEN ZETAM = SNGL(DBLE(18.964)+(QYP-12.1)/MY) QYP = DBLE(12.1) ENDIF IF (QYP.LT.-12.1) THEN ZETAM = SNGL(DBLE(18.964)+(QYP+12.1)/MY) QYP = DBLE(-12.1) ENDIF C C C EVENTI CONTENUTI SECONDO IL FIT C POSIXMD = QXP + DBLE(ZETAM)*DTAND(TMISD)*DCOSD(PMISD) POSIYMD = QYP + DBLE(ZETAM)*DTAND(TMISD)*DSIND(PMISD) C IF (POSIXMD.LE.LIMLA.AND.POSIXMD.GE.-LIMLA.AND.POSIYMD.LE. & LIMLA.AND.POSIYMD.GE.-LIMLA) EVCONMIS = EVCONMIS + 1 C C EVENTI CONTENUTI SECONDO SIA SECONDO IL FIT CHE REALMENTE CONTENUTI C IF ((POSIXMD.LE.LIMLA.AND.POSIXMD.GE.-LIMLA.AND.POSIYMD.LE. & LIMLA.AND.POSIYMD.GE.-LIMLA).AND. & (POSIXD.LE.LIMLA.AND.POSIXD.GE.-LIMLA.AND.POSIYD.LE. & LIMLA.AND.POSIYD.GE.-LIMLA)) EVCOINC = EVCOINC + 1 C C CONVERTE IN RADIANTI E IN SINGOLA PRECISIONE PER L'ENTUPLA C C ANGOLI THETA E PHI C THET = SNGL(THETA*PI/DBLE(180.)) PH = SNGL(PHI*PI/DBLE(180.)) TMIS = SNGL(TMISD*PI/DBLE(180.)) PMIS = SNGL(PMISD*PI/DBLE(180.)) C C PROIEZIONE DEGLI ANGOLI SUI PIANI XZ E YZ C IF (PHI.LT.90..OR.PHI.GT.270.) THEN THX = SNGL(DATAN2(DABS(DTAND(THETA)),DSQRT( & DBLE(1.)+DTAND(PHI)**2.))) THY = SNGL(DATAN2(DTAND(PHI)*DABS(DTAND(THETA)), & DSQRT(DBLE(1.)+DTAND(PHI)**2.))) ENDIF IF (PHI.GT.90..AND.PHI.LT.270.) THEN THX = SNGL(DATAN2(-DABS(DTAND(THETA)),DSQRT( & DBLE(1.)+DTAND(PHI)**2.))) THY = SNGL(DATAN2(-DTAND(PHI)*DABS(DTAND(THETA)), & DSQRT(DBLE(1.)+DTAND(PHI)**2.))) ENDIF C C PROIEZIONE DEGLI ANGOLI MISURATI C THXM = SNGL(DATAN(MX)) THYM = SNGL(DATAN(MY)) C C LA PROIEZIONE SULL'ULTIMO PIANO C POSIXM = SNGL(POSIXMD) POSIYM = SNGL(POSIYMD) POSIX = SNGL(POSIXD) POSIY = SNGL(POSIYD) C C LE QUOTE REALI E MISURATE C QQX = SNGL(VALORI(LENGHT-13)) QQY = SNGL(VALORI(LENGHT-12)) QQZ = ZETA QQXM = SNGL(QXP) QQYM = SNGL(QYP) QQZM = ZETAM C C REGISTRA IN ENMEZX(Y) LE ENERGIE RILASCIATE SU TUTTE LE STRIP C A DX E A SX DELLA VERA TRAIETTORIA DELLA PARTICELLA PER I PIANI X(Y) C CALL VZERO(ENMEZX,2*22) CALL VZERO(ENMEZY,2*22) CALL VZERO(QPLX,6*22) CALL VZERO(QPLY,6*22) C DO LL = 1, NPLA PMETX = 0. PMETY = 0. PY = X(LL)*TAN(THX) + QQX + 12.1 CALL TRAIE(PY,PMETX) PY = X2(LL)*TAN(THY) + QQY + 12.1 CALL TRAIE(PY,PMETY) DO JG = 1, NCHA NLK = (JG - 1)*NPLA + LL IF (PMETX.LT.0.) GO TO 446 IF (JG.LT.NINT(PMETX)) & ENMEZX(1,LL) = ENMEZX(1,LL) + ZXY(1,NLK) IF (JG.EQ.NINT(PMETX)) THEN ENMEZX(1,LL) = ENMEZX(1,LL) + ZXY(1,NLK)/2. ENMEZX(2,LL) = ENMEZX(2,LL) + ZXY(1,NLK)/2. ENDIF IF (JG.GT.NINT(PMETX)) & ENMEZX(2,LL) = ENMEZX(2,LL) + ZXY(1,NLK) 446 CONTINUE IF (PMETY.LT.0.) GO TO 447 IF (JG.LT.NINT(PMETY)) & ENMEZY(1,LL) = ENMEZY(1,LL) + ZXY(2,NLK) IF (JG.EQ.NINT(PMETY)) THEN ENMEZY(1,LL) = ENMEZY(1,LL) + ZXY(2,NLK)/2. ENMEZY(2,LL) = ENMEZY(2,LL) + ZXY(2,NLK)/2. ENDIF IF (JG.GT.NINT(PMETY)) & ENMEZY(2,LL) = ENMEZY(2,LL) + ZXY(2,NLK) 447 CONTINUE ENDDO C DO RT = 1, 6 DO RE = 1, 16 NLL = (RE+(RT-1)*16 - 1)*NPLA + LL QPLX(RT,LL) = ZXY(1,NLL) + QPLX(RT,LL) QPLY(RT,LL) = ZXY(2,NLL) + QPLY(RT,LL) ENDDO ENDDO C ENDDO C C MISURA DELL'ENERGIA: SI INTEGRA L'ENERGIA RIVELATA FINO AL PIANO DI C MASSIMO C PARZEN = 0. DO LL = 1, FINPAR PARZEN = PARZEN + PARZ(1,LL) + PARZ(2,LL) ENDDO FFLA = FLOAT(FINPAR) IF (FINPAR1.EQ.1) THEN PARZEN = PARZEN - PARZ(2,FINPAR) FFLA = FFLA - 1. ENDIF C C CORREZIONE DELL'ENERGIA PER EVENTI USCITI DAL CALORIMETRO C cc IF ((FFLA/COS(TMIS)).LT.28..AND.(POSIXMD.LT.-12.1.OR.POSIXMD cc & .GT.12.1.OR.POSIYMD.LT.-12.1.OR.POSIYMD.GT.12.1)) THEN cc PARZEN = PARZEN + PARZEN*ENCORR(FFLA/COS(TMIS)) cc ENDIF C C ERZEN = ENERGER(PARZEN) ER3 = ER3 + PARZEN*(ERZEN**(-2)) ER4 = ER4 + ERZEN**(-2) C ERENET = ENERGER(ENET) VENET = VALORI(LENGHT-14) C ER1 = ER1 + ENET*(ERENET**(-2)) ER2 = ER2 + ERENET**(-2) C ENETM = ER1/ER2 C ERENETM = ER2**(-.5) C ENM = ENM + ENET ERENM = ERENM + ERENET CCC CALL VZERO(XEN,NPL) CALL VZERO(YEN,NPL) CALL VZERO(EYEN,NPL) CALL VZERO(XEN2,NPL) CALL VZERO(YEN2,NPL) CALL VZERO(EYEN2,NPL) IFIN2 = FINPAR !NINT(FFLA2) DO NN = 1, NPLA c NPIA = NN - 3 + IFIN2 NPIA = NN IF (NPIA.GT.22) NPIA = 22 IF (NPIA.LT.0) NPIA = 0 RN = FLOAT(NPIA) c c XEN(NN) = (RN - 1.)*0.75/COS(TMIS) ! (RN - 1.)*.862 + .169 c XEN(NN+22) = RN*0.75/COS(TMIS) !(RN - 1.)*.862 + .843 c YEN(NN) = PARZ(1,NPIA) YEN(NN+22) = PARZ(2,NPIA) EYEN(NN) = 5.*SQRT(YEN(NN)) EYEN(NN+22) = 5.*SQRT(YEN(NN+22)) ENDDO NPLL = 0 DO NN = 1, NPLA IF (YEN(NN).NE.0..AND.NN.EQ.1) THEN NPLL = NPLL + 1 XEN2(NPLL) = NPLL*.76/COS(TMIS) !XEN(NN) YEN2(NPLL) = YEN(NN) EYEN2(NPLL) = EYEN(NN) ENDIF IF (YEN(NN+22).NE.0..AND.NN.EQ.22) THEN NPLL = NPLL + 1 XEN2(NPLL) = NPLL*.76/COS(TMIS) !XEN(NN+22) YEN2(NPLL) = YEN(NN+22) EYEN2(NPLL) = EYEN(NN+22) ENDIF IF (NN.GT.1.AND.NN.LT.22.AND.YEN(NN+21). & NE.0..AND.YEN(NN).NE.0.) THEN NPLL = NPLL + 1 XEN2(NPLL) = NPLL*.76/COS(TMIS) !XEN(NN) YEN2(NPLL) = (YEN(NN) + YEN(NN+21))/2. EYEN2(NPLL) = (EYEN(NN) + EYEN(NN+21))/2. ENDIF IF (NN.GT.1.AND.NN.LT.22.AND.YEN(NN+21). & NE.0..AND.YEN(NN).EQ.0.) THEN NPLL = NPLL + 1 XEN2(NPLL) = NPLL*.76/COS(TMIS) !XEN(NN) YEN2(NPLL) = YEN(NN+21) EYEN2(NPLL) = EYEN(NN+21) ENDIF IF (NN.GT.1.AND.NN.LT.22.AND.YEN(NN+21). & EQ.0..AND.YEN(NN).NE.0.) THEN NPLL = NPLL + 1 XEN2(NPLL) = NPLL*.76/COS(TMIS) !XEN(NN) YEN2(NPLL) = YEN(NN) EYEN2(NPLL) = EYEN(NN) ENDIF ENDDO IF (NPLL.LT.13..OR.BETA.NE.'Y') THEN PAR(1) = -1. PAR(2) = -1. POSMAX = -1. PAR1 = -1. GO TO 448 ENDIF CHOPT = 'BEQ' PAR(1) = 1. PAR(2) = 1. PMI2(1) = 2. PMI2(2) = -0.5 !.4 PMA(1) = 20. !11. PMA(2) = 6. STE(1) = 0.01 STE(2) = 0.01 C CALL HFITV(NPLL,NPL,NVAR,XEN2,YEN2,EYEN2,ENERGIF,CHOPT, & NPAR2,PAR,STE,PMI2,PMA,SIGPA,CHI) CCC IF (PAR(2).NE.0.) THEN POSMAX = PAR(1)/PAR(2) PAR1 = PAR(1) ELSE PAR1 = PAR(1) POSMAX = 0. ENDIF 448 CONTINUE C C RIEMPIE LE ENTUPLE C CALL HFNT(1) C TRIG = TRIG + 1 C c IF ((THET-TMIS).LT.-.3) THEN c PRINT *,'FLA2 ',FLA2 c PRINT *,'FLA ',FLA c PRINT *,'QXP ',DBLE(VALORI(LENGHT-13)) c PRINT *,'QXPM ',QXP c PRINT *,'QYP ',DBLE(VALORI(LENGHT-12)) c PRINT *,'QYPM ',QYP c PRINT *,'THETA: ',THETA c PRINT *,'TMIS ',DATAND(DSQRT((MX**2) + (MY**2))) C PRINT *,'thet: ',thet c PRINT *,'TMIS: ',TMISD c PRINT *,'PHI: ',PHI c PRINT *,'PMIS ',PMISD c PRINT *,'TRIG ',TRIG c PRINT *,'***********' c ENDIF C 449 CONTINUE C C ESEGUE PER TUTTI GLI EVENTI C 450 ENDDO C CALL HFNT(2) 500 CONTINUE 900 CONTINUE C CLOSE(83) C C PRINT RESULTS C PRINT *,'FILE PROCESSATO: ',FILENAME PRINT *,'FILE DATI SALVATO: ',NAME PRINT *,'EVENTI TOTALI: ',EVENTO PRINT *,'EVENTI CONTENUTI REALEMENTE: ',EVCON PRINT *,'EVENTI ANALIZZATI: ',TRIG PRINT *,'EVENTI CONTENUTI MISURATI: ',EVCONMIS PRINT *,'EVENTI CONTENUTI COINCIDENTI: ',EVCOINC PRINT *,'EVENTI BUONI: ',TOTTO2 PRINT *,'EVENTI SCARTATI: ',SCARTA c PRINT *,'ENERGIA MISURATA PESATA: ',ENETM c PRINT *,'ERRORE EN. PESATA: ',ERENETM C PRINT *,'ENERGIA REALE ',VENET/CALIB C PRINT *,'ENERGIA MEDIA MISURATA: ',ENM/TRIG C PRINT *,'ERRORE MEDIO MISURATO: ',ERENM/TRIG c PRINT *,'EN. MIS. PARZIALE PESATA: ',ER3/ER4 c PRINT *,'ERRORE EN. PARZ. PESATA: ',ER4**(-.5) C 950 CONTINUE C CALL HROUT(0,ICYCLE,' ') CALL HREND('SIMULATION') 1000 CONTINUE C 5050 FORMAT(A40) 5060 FORMAT('Y',I2) 5070 FORMAT('NUEV') 5080 FORMAT(1X,I5) 6001 FORMAT(A1) 6005 FORMAT(F3.3) 6006 FORMAT(F4.3) C END C C--------------------------------------------------------------------- SUBROUTINE LEGGE(VALORI,LENGHT) C--------------------------------------------------------------------- C REAL VALORI(LENGHT) C READ(83),VALORI C RETURN END C C--------------------------------------------------------------------- SUBROUTINE BORDO(NB) C--------------------------------------------------------------------- C PARAMETER (STMA=23.616,STMI=.584) PARAMETER (NPLA=22) INTEGER CONF, CONF2, NB, STORE REAL YV C COMMON/VETFIT/YV(NPLA) C C ROUTINE PER SCARTARE DAL FIT I PUNTI CHE SCORRONO SUL BORDO DEL C CALORIMETRO E CHE NON SONO IL PROSEGUIMENTO DELLA VERA TRAIETTORIA C CONF = 0 CONF2 = 0 DO KK = 1, NB IF (YV(KK).LT.STMA.OR.YV(KK).GT.STMI) CONF = 1 IF ((YV(KK).GE.STMA.OR.YV(KK).LE.STMI).AND.CONF.EQ.1 & .AND.CONF2.EQ.0) STORE = KK IF (YV(KK).GE.STMA.OR.YV(KK).LE.STMI.AND.CONF.EQ.1) & CONF2 = CONF2 + 1 IF (CONF2.GE.2) THEN NB = STORE GO TO 7000 ENDIF ENDDO C 7000 CONTINUE RETURN END C C--------------------------------------------------------------------- SUBROUTINE TRAIE(YH,TRH) C--------------------------------------------------------------------- C REAL YH, TRH C IF (YH.GE.0..AND.YH.LE.24.2) THEN IF (YH.LE..218) TRH = 1. IF (YH.GT..218.AND.YH.LT.7.782) & TRH = 1. + ( YH - .218) / .244 IF (YH.LT.8.05.AND.YH.GE.7.782) & TRH = 32. IF (YH.LE.8.318.AND.YH.GT.8.05) & TRH = 33. IF (YH.GT.8.318.AND.YH.LT.15.882) & TRH = 1. + (YH - .51) / .244 IF (YH.LT.16.15.AND.YH.GE.15.882) & TRH = 64. IF (YH.LE.16.418.AND.YH.GT.16.15) & TRH = 65. IF (YH.GT.16.418.AND.YH.LT.23.982) & TRH = 1. + (YH - .802) / .244 IF (YH.GE.23.982) TRH = 96. ELSE TRH = -10. ENDIF RETURN END C C--------------------------------------------------------------------- SUBROUTINE BARIC(TRA,YY,EYY,NNF,FF,NX) C--------------------------------------------------------------------- C PARAMETER (STRIP1=7.,STRIP2=5.) !!7 5 PARAMETER (NPLA=22, NCHA=96) PARAMETER (LENSEV=NPLA*NCHA) C REAL TRA, YY, EYY INTEGER ST0, ST1, ST2, ST3 REAL UFF, ENTO REAL ZXY, YN C INTEGER NLK, NNF, FF, NX, NPOS C COMMON/ENSTRIP/ZXY(2,LENSEV),YN(NCHA) C ST0 = NINT(STRIP1) ST1 = NINT(STRIP1-STRIP2) ST2 = NINT((STRIP1-1.)/2. + 1.) ST3 = NINT(FLOAT(ST1)/2.) C IF (TRA.GT.0.) THEN UFF = 0. ENTO = 0. DO P = 1, (ST0-(FF-1)*ST1) IF (NINT(TRA).LT.(ST2-ST3*(FF-1))) THEN IF (P.GT.(ST0-ABS(NINT(TRA)-ST2-ST3*(FF-1)))) GO TO 7100 NPOS = P ELSE IF ((NINT(TRA)+P-ST2+ST3*(FF-1)).GT.96) GO TO 7100 NPOS = NINT(TRA)+P-ST2+ST3*(FF-1) ENDIF NLK = (NPOS-1)*NPLA + NNF UFF = UFF + YN(NPOS)*ZXY(NX,NLK) ENTO = ENTO + ZXY(NX,NLK) ENDDO 7100 CONTINUE IF (ENTO.GT.0.) THEN YY = UFF/ENTO c EYY = 1. EYY = YY/(ENTO**.79) c IF (NNF.GT.8.AND.NNF.LT.20) THEN cc EYY = YY/(ENTO**(.85 + 2.*EXP(-NNF*1.5))) c EYY = YY/(ENTO**(.71)) c ELSE c EYY = YY/(ENTO**(.89)) c ENDIF ELSE YY = -1. EYY = 1E5 ENDIF IF (NNF.LT.2.AND.YY.GT.0.) EYY = .244/SQRT(12.) ccc EYY = .244/SQRT(12.) ccc ELSE YY = -1. EYY = 1E5 ENDIF RETURN END C C--------------------------------------------------------------------------- REAL FUNCTION FUN(X) C--------------------------------------------------------------------------- C PARAMETER (NPAR=2, NPAR2=2) PARAMETER (NWPAWC=3000000) COMMON/HCFITD/PA(NPAR),SIGPAR(NPAR),CHI2,PMIN(NPAR),PMAX(NPAR), & STEP(NPAR), & PAR(NPAR2),SIGPA(NPAR2),CHI,PMI2(NPAR2),PMA(NPAR2), & STE(NPAR2) COMMON/PAWC/HMEMOR(NWPAWC) C FUN = PA(1) + X*PA(2) C END C C--------------------------------------------------------------------------- REAL FUNCTION ENERGER(X) C--------------------------------------------------------------------------- C PARAMETER (CALIB=0.0001059994) C ENERGER = .0205 * SQRT(7.4 * 6.76/(X * CALIB)) * X C END C C--------------------------------------------------------------------------- REAL FUNCTION ENCORR(X) C--------------------------------------------------------------------------- C REAL FFUN C FFUN = .5+.353*ATAN((X-21.972)*.374) c FFUN = .5466+.3831*ATAN((X-21.386)*.343) IF (X.LT.28.) THEN IF (FFUN.LT.0.) FFUN = 1.E-5 ENCORR = .45/FFUN ELSE ENCORR = 0. ENDIF C END C C--------------------------------------------------------------------------- REAL FUNCTION ENERGIF(X) C--------------------------------------------------------------------------- PARAMETER (NPAR=2, NPAR2=2) PARAMETER (NWPAWC=3000000) COMMON/HCFITD/PA(NPAR),SIGPAR(NPAR),CHI2,PMIN(NPAR),PMAX(NPAR), & STEP(NPAR), & PAR(NPAR2),SIGPA(NPAR2),CHI,PMI2(NPAR2),PMA(NPAR2), & STE(NPAR2) COMMON/PAWC/HMEMOR(NWPAWC) C C Rossi equation B C ENERGIF = (0.8*(X**(PAR(1)))+0.2)*EXP(-X*PAR(2)) C END