      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

