*************************************************************************
*     This version of momanhough performs tracker CALIBRATIONS:
*     - eta functions
*     - charge correlation
*     - ADC-to-MIP conversion
*     -------------------------------------------------------------------
*     It works as it follow:
*     - it reads cluster information (read a list of LEVEL1 ntuples, 
*       output of reduction.f)
*     - it performs track identification, select good tracks (6 hit planes) 
*       and perform first-order fit
*     - it calculate calibration variables
*     - it creates a LEVEL2 ntuple with selected events (in order to check)
*
*************************************************************************
      program momanhough

      include '../common/commontracker.f'
      include '../common/common_momanhough.f'
      include '../common/momanhough_init.f'
      include '../common/common_level2debug.f'
      include '../common/common_mech.f'
      include '../common/common_xyzPAM.f'
      include '../common/common_mini_2.f'
      include '../common/calib.f'
      include '../common/level1.f'
      include '../common/level2.f'
      include '../common/common_calibration.f'

*     flag to set debug mode
      logical DEBUG
      common/dbg/DEBUG
      logical DONTDO
      data DONTDO/.true./
c$$$*     flag to chose PFA
c$$$      character*10 PFA
c$$$      common/FINALPFA/PFA


c      parameter (inf=1.e8)      !just a huge number...
        
*     external functions
      external npl
      external acoordsi,coordsi,nld,coord,dcoord
c------------------------------------------------------------------------
c     
c     local variables
c     
c------------------------------------------------------------------------

      character*24 processing_date

      parameter (nfile_max=500)
      parameter (lun_data_level1=71) !data file id number
      parameter (lun_data_level2=72) !data file id number
      parameter (lun_data_calib=74)  !data file id number
c      parameter (lun_file_out=75)    !data file id number

      character*74 data_file    !data file name
      character*74 data_dir     !
      character*74 out_dir      !
      character*74 data_file_calib
      character*74 data_file_level1
      character*74 data_file_level2
      character*10 input_string(nfile_max)
      character*74 file_out
      character*74 fname_param
      character*7 rzdir



# ifndef TEST2003
      parameter(ncalibmax=50)
      character*40 file_calib(ncalibmax)
      parameter(lun_calib_list=66) !calibration list file id
      integer which_calib_last
# endif

      integer tottr
      integer minevent          !first event to be analysed
      logical FIMAGE            !


      COMMON/QUEST/IQUEST(100)  
c     !permette di ottenere ntuple funzionanti nonostante
c     ! il messaggio dei 64K di RZOUT...!???

# ifdef TEST2003
      parameter (nang=3)       !number of angular bins
      parameter (angstep=4.)
      parameter (angmin=-6.)
      print*,' '
      print*,'Set angular binning TEST 2003'
      print*,' '

      nangbin=nang
      do ia=1,nangbin
         angL(ia)=angmin+angstep*(ia-1)
         angR(ia)=angmin+angstep*ia
         print*,ia,' - ',angL(ia),angR(ia)
      enddo
      print*,' '


# else
c$$$      parameter (nang=21)       !number of angular bins
c$$$      parameter (angstep=2.)
c$$$      parameter (angmin=-21.)
c$$$      print*,' '
c$$$      print*,'Set angular binning Rome-2005'
c$$$      print*,' '
c$$$
c$$$      nangbin=nang
c$$$      do ia=1,nangbin
c$$$         angL(ia)=angmin+angstep*(ia-1)
c$$$         angR(ia)=angmin+angstep*ia
c$$$         print*,ia,' - ',angL(ia),angR(ia)
c$$$      enddo
c$$$      print*,' '

      parameter (nang=17)       !number of angular bins
*     NB angL and angR are created with length 21
      data angL
     $     /-21.,-17.,-13.,-11.,-9.,-7.,-5.,-3.,-1.
     $     ,1.,3.,5.,7.,9.,11.,13.,17.,0.,0.,0.,0./
      data angR
     $     /-17.,-13.,-11.,-9.,-7.,-5.,-3.,-1.
     $     ,1.,3.,5.,7.,9.,11.,13.,17.,21.,0.,0.,0.,0./
      print*,' '
      print*,'Set angular binning Rome-2005'
      print*,' '
      nangbin=nang
      do ia=1,nangbin
         print*,ia,' - ',angL(ia),angR(ia)
      enddo
      print*,' '

# endif

c------------------------------------------------------------------------
c     
c     HBOOK initialization
c     
c------------------------------------------------------------------------

      call HLIMIT(NWPAWC)


c------------------------------------------------------------------------
c     
c     reads input informations
c     
c------------------------------------------------------------------------
      call fdate(processing_date)
      write(*,101)
     $     processing_date
 101  format(/
     $      ,'*** *** *** *** *** *** *** *** *** *** *** *** ***',/
     $      ,'*                                                 *',/
     $      ,'* CALIBRATION                                     *',/     
     $      ,'*                                                 *',/
     $      ,'*** *** *** *** *** *** *** *** *** *** *** *** ***',/
     $     ,a24,/
     $     )

      print*,'Number of files to be analysed:'
      read(*,*)nfile
      print*,nfile
 111  format(a)
      print*,'Input data directory:'      
      read(*,111)data_dir
      print*,data_dir
      print*,'Output data directory:'      
      read(*,111)out_dir
      print*,out_dir
      print*,'List of file identifiers: (DATE_NUM)'
 300  format(A10)
      do i=1,nfile
         read(*,300)input_string(i)
         print*,input_string(i)
      enddo
      minevent=1
      print*,'Maximum number of events to be analized:' !20000
      read(*,*) ntotev
      print*,ntotev    
      print*,'DEBUG mode: (T, F)'
 11   format(l1)
      read(*,11)DEBUG
      print*,DEBUG
      print*,'---------------------------------------------------'


****** INITIALIZATIONS *************************************


*     2) read z coordinates of the planes
      print*,' '
      print*,'- read z coordinates of the planes'
      print*,' '
      call mech_sensor          !reads sensors centres coordinates
      do ip=1,nplanes
        fitz(ip)=z_mech_sensor(ip,1,1)*0.1 !cm
*     gets planes mechanical z positions
*     (in mm) and sets them in micrometers
      enddo                     

*     4) read magnetic field map      
      print*,' '
      print*,'- read magnetic field map'
      print*,' '
      call read_B

*     5) read allignment parameters
      print*,' '
      print*,'- read aligment parameters'
      print*,' '
      call readalignparam

************************************************************
************************************************************
************************************************************
************************************************************

      print*,' '
      print*,'---------------------------------------------------'
      print*,' '
      print*,'OPENING OUTPUT-HISTO FILES:'

 1010 format('LADDER',i1)
      do ilad=1,3

         write(rzdir,1010)ilad
         file_out=
     $        out_dir(1:LNBLNK(out_dir))//'trk-'
     $        //rzdir//'-histos.rz'
         print*,' -> '//file_out
         IQUEST(10)=65000
         call HROPEN(lun_file_out+ilad,rzdir,file_out,'PQN',4096,istat)
         if(istat.ne.0) goto 16 
         call HCDIR('//'//rzdir,' ')
         
         call HCDIR('//PAWC','')
         call HMDIR(rzdir,'S')

         call book_eta_histos
         call book_charge_histos
         call book_check_histos
         call HCDIR('//PAWC','')
         
      enddo

************************************************************
************************************************************
************************************************************
************************************************************
************************************************************
************************************************************
************************************************************
************************************************************
************************************************************
************************************************************
************************************************************
************************************************************

c------------------------------------------------------------------------
c     
c     LEVEL 2 ntuple booking
c     
c------------------------------------------------------------------------


 503  format(a,'DW__level2-calib.rz')
      write(data_file_level2,503)
     $     out_dir(1:LNBLNK(out_dir))
c     $     ,data_file(1:LNBLNK(data_file))
      
      print*,''
      print*,'------------------------------------'
      print*,'  Creating LEVEL2 rz file'
      print*,'  ',data_file_level2
      print*,'------------------------------------'
      

C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C     largest RZ file: IQUEST(10) records x LREC words x 4 byte
C     with the following settings: 65000 x 4096 x 4 = 1G
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
      IQUEST(10)=65000          
c     !permette di ottenere ntuple funzionanti nonostante
c     ! il messaggio dei 64K di RZOUT...!???
      call HROPEN(lun_data_level2,
     $     'LEVEL2',data_file_level2,'QNP',4096,istat) !opens rz  
      if(istat.ne.0) goto 19      
      call book_level2
         
      print*,' '
      print*,'---------------------------------------------------'
      print*,'           * * * START ANALYSIS * * *'
      print*,'---------------------------------------------------'
 
c------------------------------------------------------------------------
c     
c     loop on files
c     
c------------------------------------------------------------------------

      goodtr = 0
      tottr = 0
      do ifile = 1,nfile        !files loop

         print*,' '
         data_file = input_string(ifile)
         print*,'--> FILE ',ifile,': '//data_file(1:LNBLNK(data_file))

c------------------------------------------------------------------------
c     
c     fills bad variables with online DSP bad strips from calib histograms
c     for bad strip exclusion
c     
c------------------------------------------------------------------------
        
# ifdef TEST2003
c------------------------------------------------------------------------
c     
c     opens calib file
c     
c------------------------------------------------------------------------

 505     format(a,'output_',a,'_calib.rz')
         write(data_file_calib,505)
     $        data_dir(1:LNBLNK(data_dir))
     $        ,data_file(1:LNBLNK(data_file))
         
         print*,' '
         print*,'OPENING CALIB FILE...', data_file_calib
         
         IQUEST(10)=65000
         call HROPEN
     $        (lun_data_calib,'IN',data_file_calib,'QP',4096,istat)
         if(istat.ne.0) goto 17
         call HCDIR('//IN',' ')
         
         call HRIN(0,9999,0)    !puts histograms in memory
         
         print*,' '
         
         print*,' '
         print*,'READING PEDESTAL, SIGMA AND BADSTRIP HISTOGRAMS...'
         print*,' '
         
         call fillpedsig
         
         call HREND('IN')
         close (lun_data_calib)  
# else
      
         print*,' '
         print*,'OPENING CALIBRATION-LIST FILE:'
 501     format(a,'DW_',a,'_calib.txt')
         write(data_file_calib,501)
     $        data_dir(1:LNBLNK(data_dir))
     $        ,data_file(1:LNBLNK(data_file))
         print*,data_file_calib
         open(lun_calib_list,
     $        FILE=data_file_calib(1:LNBLNK(data_file_calib))
     $        ,STATUS='UNKNOWN'
     $        ,IOSTAT=iostat
     $        )
         
 113     format(i5,'    ',a25)
         
         do i=1,ncalibmax
            
            read(lun_calib_list,113,IOSTAT=iostat) n_cal_list
     $           ,file_calib(i)
c            file_calib(i)=file_calib(i)(1:LNBLNK(file_calib(i)))
            if(iostat.ne.0)then
               ncal=i-1
               goto 2000
            endif
            print*,n_cal_list,' - '
     $           ,file_calib(i)(1:LNBLNK(file_calib(i)))
            
         enddo
         
 2000    close(lun_calib_list)
         which_calib_last=0
         
# endif 
      
c------------------------------------------------------------------------
c     
c     opens level1 file
c     
c------------------------------------------------------------------------

# ifdef TEST2003
 504     format(a,'output_',a,'_level1.rz')
# else
 504     format(a,'DW_',a,'_level1.rz')
# endif
         write(data_file_level1,504)
     $      data_dir(1:LNBLNK(data_dir))
     $     ,data_file(1:LNBLNK(data_file))
         print*,''
         print*,'OPENING LEVEL1 FILE:'
         print*,data_file_level1
         IQUEST(10)=65000          
c     !permette di ottenere ntuple funzionanti nonostante
c     ! il messaggio dei 64K di RZOUT...!???
         call HROPEN(lun_data_level1,
     $        'LEVEL1',data_file_level1,'QP',4096,istat) !opens rz  
         if(istat.ne.0) goto 19      
         call HRIN(ntp_level1,9999,20)       
         call access_level1
c     call HPRNTU(ntp_level1+20)  
         call HNOENT(ntp_level1+20,iemax0)
         print*,'( ',iemax0,' events)'
         
         
************************************************************
************************************************************
************************************************************
*     
*     start track analysis
*
************************************************************
************************************************************
************************************************************
         nselect = 0
         maxevent=minevent+ntotev-1
         do iev = minevent,MIN(iemax0,maxevent) !loop on events
            
            call HCDIR('//LEVEL1',' ')
            call HGNT(ntp_level1+20,iev,ierr) !reads an event
            if(ierr.ne.0) goto 21

            tottr = tottr +1

*------------------------------------------------------
*        LEVEL2 N-TUPLE INITIALIZATIONS
            call init_level2
            if(.not.good1)then
*               goto 8800        !fill nt-uple and go to next event
               goto 100        !go to next event
            endif
*------------------------------------------------------


# ifndef TEST2003
      
            if(which_calib.ne.which_calib_last.and.
     $which_calib.ne.0)then
               
               data_file_calib=
     $              data_dir(1:LNBLNK(data_dir))//
     $              file_calib(which_calib)
     $              (1:LNBLNK(file_calib(which_calib)))
c     print*,data_file_calib
               print*,''
               print*,
     $              '@ event ',nev2
     $              ,'   (CPU pkt N.',pkt_num1,')'
               print*,'--> ',data_file_calib            
               IQUEST(10)=65000          
               call HROPEN(lun_data_calib,
     $              'CALIB',data_file_calib,'QP',4096,istat) !opens
               if(istat.ne.0) goto 19      
               call HRIN(0,9999,0)       
               call fillpedsig
               do iview=1,nviews
                  call HDELET(id_hi_bad+iview)
                  call HDELET(id_hi_ped+iview)
                  call HDELET(id_hi_sig+iview)
               enddo
               call HREND('CALIB')
               close(lun_data_calib)            
               which_calib_last=which_calib           
               
            elseif(which_calib.eq.0)then
               
               nocalib=nocalib+1
               good2=.false.
               goto 8800        !fill nt-uple and go to next event
               
            endif
            
# endif

      
*------------------------------------------------------
*     cut on maximum number of clusters
*------------------------------------------------------
*            if(nclstr1.gt.nclstrmax_level2)then
            if(nclstr1.gt.20)then
*               goto 8800        !fill nt-uple and go to next event
               goto 100        !go to next event
            endif
            
            do i=1,nclstr1
               cl_used(i)=0     !init mask of clusters associated to a track
            enddo
            
            if(DEBUG)then
               print*,'----------------------------------'
               print*,iev,'   ** ',nev2
            endif
            
*-------------------------------------------------------------------------------
*     STEP 1
*-------------------------------------------------------------------------------
*     X-Y cluster association
*     
*     Clusters are associated to form COUPLES
*     Clusters not associated in any couple are called SINGLETS
*
*     Track identification (Hough transform) and fitting is first done on couples. 
*     Hence singlets are possibly added to the track.
*    
*     Variables assigned by the routine "cl_to_couples" are those in the 
*     common blocks:
*     - common/clusters/cl_good
*     - common/couples/clx,cly,ncp_plane,ncp_tot,cp_useds1,cp_useds2
*     - common/singlets/ncls,cls,cl_single
*
*     NB - The routine version used for calibration is different from that used
*     for event processing. Different cluster selection is applied: 
*     - "ngoodstr" good strip around MAXS required (ngoodstr=2 => 5 good strips)
*     - no charge correlation required
*-------------------------------------------------------------------------------
*-------------------------------------------------------------------------------

            iflag=0
            call cl_to_couples_nocharge(iflag)
            if(iflag.eq.1)then  !bad event
*               goto 880         !fill ntp and go to next event            
               goto 100         !go to next event            
            endif
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*     PRE-SELEZIONE DI TRACCE PULITE 
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*     at least a couple per plane
            do ip=1,nplanes 
               if(ncp_plane(ip).eq.0)goto 100 !next event 
            enddo 
c     if(good2.eq.0)goto 880!fill ntp and go to next event 
            
*-----------------------------------------------------
*-----------------------------------------------------
*     HOUGH TRASFORM
*-----------------------------------------------------
*-----------------------------------------------------


*-------------------------------------------------------------------------------
*     STEP 2
*-------------------------------------------------------------------------------
*     
*     Association of couples to form
*     - DOUBLETS in YZ view
*     - TRIPLETS in XZ view
*    
*     Variables assigned by the routine "cp_to_doubtrip" are those in the 
*     common blocks:
*     - common/hough_param/
*     $     alfayz1,  !Y0
*     $     alfayz2,  !tg theta-yz
*     $     alfaxz1,  !X0
*     $     alfaxz2,  !tg theta-xz
*     $     alfaxz3   !1/r
*     - common/doublets/ndblt,cpyz1,cpyz2
*     - common/triplets/ntrpt,cpxz1,cpxz2,cpxz3
*-------------------------------------------------------------------------------
*-------------------------------------------------------------------------------

            iflag=0
            call cp_to_doubtrip(iflag)
            if(iflag.eq.1)then  !bad event
*               goto 880         !fill ntp and go to next event            
               goto 100         !go to next event            
            endif


*-------------------------------------------------------------------------------
*     STEP 3
*-------------------------------------------------------------------------------
*     
*     Classification of doublets and triplets to form CLOUDS,  
*     according to distance in parameter space. 
*     
*     cloud = cluster of points (doublets/triplets) in parameter space
*
*     
*    
*     Variables assigned by the routine "doub_to_YZcloud" are those in the 
*     common blocks:
*     - common/clouds_yz/                   
*     $     nclouds_yz                        
*     $     ,alfayz1_av,alfayz2_av           
*     $     ,ptcloud_yz,db_cloud,cpcloud_yz
*
*     Variables assigned by the routine "trip_to_XZcloud" are those in the 
*     common blocks:
*      common/clouds_xz/                   
*     $      nclouds_xz                        
*     $     ,alfaxz1_av,alfaxz2_av,alfaxz3_av 
*     $     ,ptcloud_xz,tr_cloud,cpcloud_xz           
*-------------------------------------------------------------------------------
*-------------------------------------------------------------------------------

            iflag=0
            call doub_to_YZcloud(iflag)
            if(iflag.eq.1)then  !bad event
*               goto 880         !fill ntp and go to next event            
               goto 100         !go to next event            
            endif
            iflag=0
            call trip_to_XZcloud(iflag)
            if(iflag.eq.1)then  !bad event
*               goto 880         !fill ntp and go to next event            
               goto 100         !go to next event            
            endif
            

*-------------------------------------------------------------------------------
*     STEP 4   (ITERATED until any other physical track isn't found)
*-------------------------------------------------------------------------------
*     
*     YZ and XZ clouds are combined in order to obtain the initial guess
*     of the candidate-track parameters.
*     A minimum number of matching couples between YZ and XZ clouds is required.
*
*     A TRACK CANDIDATE is defined by 
*     - the couples resulting from the INTERSECTION of the two clouds, and 
*     - the associated track parameters (evaluated by performing a zero-order  
*       track fitting) 
*
*     The NTRACKS candidate-track parameters are stored in common block:
*    
*     - common/track_candidates/NTRACKS,AL_STORE
*     $     ,XV_STORE,YV_STORE,ZV_STORE
*     $     ,XM_STORE,YM_STORE,ZM_STORE
*     $     ,RESX_STORE,RESY_STORE
*     $     ,AXV_STORE,AYV_STORE
*     $     ,XGOOD_STORE,YGOOD_STORE
*     $     ,CP_STORE,RCHI2_STORE
*
*-------------------------------------------------------------------------------
*-------------------------------------------------------------------------------
            ntrk=0              !counter of identified physical tracks
            
11111       continue            !<<<<<<< come here when performing a new search
            
            iflag=0
            call clouds_to_ctrack(iflag)
            if(iflag.eq.1)then  !no candidate tracks found
*               goto 880         !fill ntp and go to next event   
               goto 100         !go to next event            
            endif


*------------------------------------------------------
*     cut on maximum number of track candidates
*     (to have clean events)
*     NB: must be at least 2 to do not discard tracks
*     contained in a sensor column on Y view     
*
*     cut replaced with cut on the # cluster out of the track
*------------------------------------------------------
*            if(ntracks.gt.2)then
*               goto 8800        !fill nt-uple and go to next event
*               goto 100        !go to next event
*            endif
            

*            if(DEBUG)print*,'*GOOD!*'

            FIMAGE=.false.      !processing best track (not track image)
            ibest=0             !best track among candidates
            iimage=0            !track image
*     -----------------------------------------------------
*     --------------- select the best track ---------------
*     selection cuts:
*     - NX = NY = 6
*     -----------------------------------------------------
*
            rchi2best=100000.      
            do i=1,ntracks

               npx=0
               npy=0
               do ip=1,nplanes
                  npx=npx+XGOOD_STORE(ip,i)
                  npy=npy+YGOOD_STORE(ip,i)
               enddo
*               print*,
*     $         'Candidate ',i,' nx,ny ',npx,npy,' chi2 ',RCHI2_STORE(i)
               if(npx.eq.nplanes.and.npy.eq.nplanes)then
                  if(RCHI2_STORE(i).lt.rchi2best.and.
     $                 RCHI2_STORE(i).gt.0)then
                     ibest=i
                     rchi2best=RCHI2_STORE(i)
                  endif
               endif
            enddo
*     if(ibest.eq.0)goto 880 !>> no good candidates 
            if(ibest.eq.0)goto 100 !>> no good candidates 

*            print*,'BEST TRACK ',ibest,'  nx,ny',npx,npy

*            print*,'>>>> ',iev

*     -----------------------------------------------------
*     ------------- store the best track only -------------
*     -----------------------------------------------------
            ntr=1
            
            ntrk=ntr
            image(ntr)=0
            good2=.true.
            chi2_nt(ntr)        = RCHI2_STORE(ibest)
c     print*,chi2_nt(ntr)

*     ****************************************************
*     normalize the angles
*     ****************************************************
            pig=acos(-1.)
            phi   = AL_STORE(4,ibest)       
            sinth = AL_STORE(3,ibest)       
            if(sinth.lt.0)then  
               sinth = -sinth   
               phi = phi + pig  
            endif               
            npig = aint(phi/(2*pig)) 
            phi = phi - npig*2*pig 
            if(phi.lt.0)        
     $           phi = phi + 2*pig 
            AL_STORE(4,ibest) = phi         
            AL_STORE(3,ibest) = sinth       
*     ****************************************************
*     ... then assign the other variables
*     ****************************************************
            do i=1,5
               al_nt(i,ntr)     = AL_STORE(i,ibest)
c               print*,i,'al_nt(i,ntr)',al_nt(i,ntr)
               do j=1,5
                  coval(i,j,ntr) = 0 ! TEMPORANEO!!!!!!!!!!!!
               enddo
            enddo
            do ip=1,nplanes     ! loop on planes
               xgood_nt(ip,ntr) = XGOOD_STORE(ip,ibest)
               ygood_nt(ip,ntr) = YGOOD_STORE(ip,ibest)
               xm_nt(ip,ntr)    = XM_STORE(ip,ibest)
               ym_nt(ip,ntr)    = YM_STORE(ip,ibest)
               zm_nt(ip,ntr)    = ZM_STORE(ip,ibest)
               RESX_nt(IP,ntr)  = RESX_STORE(ip,ibest)
               RESY_nt(IP,ntr)  = RESY_STORE(ip,ibest)
               xv_nt(ip,ntr)    = XV_STORE(ip,ibest)
               yv_nt(ip,ntr)    = YV_STORE(ip,ibest)
               zv_nt(ip,ntr)    = ZV_STORE(ip,ibest)
               axv_nt(ip,ntr)   = AXV_STORE(ip,ibest)
               ayv_nt(ip,ntr)   = AYV_STORE(ip,ibest)
               
c               print*,ip,' xm_nt(ip,ntr)', xm_nt(ip,ntr)
c               print*,ip,' ym_nt(ip,ntr)', ym_nt(ip,ntr)
c               print*,ip,' zm_nt(ip,ntr)', zm_nt(ip,ntr)
c               print*,ip,' RESX_nt(ip,ntr)', RESX_nt(ip,ntr)
c               print*,ip,' RESY_nt(ip,ntr)', RESY_nt(ip,ntr)
c                print*,ip,' xv_nt(ip,ntr)',xv_nt(ip,ntr)
c                print*,ip,' yv_nt(ip,ntr)',yv_nt(ip,ntr)
c                print*,ip,' zv_nt(ip,ntr)',zv_nt(ip,ntr)
c                print*,ip,' axv_nt(ip,ntr)',axv_nt(ip,ntr)
c                print*,ip,' ayv_nt(ip,ntr)',ayv_nt(ip,ntr)

               id = CP_STORE(ip,ibest)
               icp = icp_cp(id)
               iclx = clx(nplanes-ip+1,icp)
               icly = cly(nplanes-ip+1,icp)
               which_cl_x(ip) = iclx
               which_cl_y(ip) = icly
               dedx_x(ip,ntr)   = dedx(iclx) 
               dedx_y(ip,ntr)   = dedx(icly)         
            enddo

c            call CalcBdL(100,xxxx,IFAIL)
c            if(ifps(xxxx).eq.1)BdL(ntr) = xxxx
            BdL(ntr) = 0

*     --- then remove selected clusters (ibest+iimage) from clouds ----
            call clean_XYclouds(ibest)
            if(iflag.eq.1)then  !bad event
*               goto 880         !fill ntp and go to next event            
               goto 100        !go to next event
            endif

*     >>>>>>>>>>>>>>>>>>> NT-UPLE filling <<<<<<<<<<<<<<<<<<<<
*       +   +   +   +   +   +   +   +   +   +   +   +   +   +   +   +   +
*      / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ /
*     +   +   +   +   +   +   +   +   +   +   +   +   +   +   +   +   +   

 880        continue
*     **********************************************************
*     stores info about clusters not associated with any track
*     **********************************************************
c            call fill_level2_siglets

            if(DEBUG)then
               
               print*,''
               print*,'DONE!'
               print*,''
               print*,'* summary *'
               print*,'tracks    ',ntrk
               print*,'cl used   ',(cl_used(i),i=1,nclstr1)
c$$$               print*,'cl unused (x-y)'
c$$$               do ip=1,nplanes
c$$$                  print*,ip,' << ',nclsx(ip),nclsy(ip)
c$$$               enddo
               print*,''
               print*,''
            endif
            
            
 8800       continue


*     -----------------------------------------------------
*     +++++++++++++++++++++++++++++++++++++++++++++++++++++
*     calibrations:
*     +++++++++++++++++++++++++++++++++++++++++++++++++++++
*     -----------------------------------------------------


*     ==================================
*     //////////////////////////////////
*     Track selection
*     \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
*     ==================================
*     - cut on number of planes
*      (fatto a monte... non mi ricordo perche')
c$$$            npx=0
c$$$            npy=0
c$$$            do ip=1,nplanes
c$$$               npx=npx+XGOOD_STORE(ip,i)
c$$$               npy=npy+YGOOD_STORE(ip,i)
c$$$            enddo
c$$$            if(.not.(npx.ge.5..and.npy.ge.5))goto 100
*     - chi^2 selection 
            rig = abs(1./al_nt(5,1))
c$$$            chicut = 5.         !100.
c$$$            chicut = chicut * ( 1.48 + 4.64/rig + 18.35/rig**2 )
c$$$            if(chi2_nt(1).gt.chicut)goto 100 !next event
            chicut = 20
            chicut = chicut *( 
c     $            3.9633*exp(-0.54485E-01*rig)
c     $           +0.99231
c     $           +10.835*exp(-0.67155*rig) )
     $            3.7032*exp(-0.41911E-01*rig)
     $           +0.85355
     $           +10.905*exp(-0.61893*rig) )
            if(chi2_nt(1).gt.chicut.or.chi2_nt(1).lt.0)goto 100 !next event
*     - cut on number of cluster out of the track (noise)
*      (no constraints on the last plane becouse of the 
*       back-scattering from calorimeter)
c$$$            do ip=1,5           !6
c$$$               nclstot=nclsx(ip)+nclsy(ip)
c$$$               if(nclstot.ge.1)goto 100 !next event
c$$$            enddo
            nclstot = 0
            do isx=1,nclsx
c               if(planex(isx).ne.6)nclstot =  nclstot +1
               nclstot =  nclstot +1
            enddo
            do isy=1,nclsy
c               if(planey(isy).ne.6)nclstot =  nclstot +1
               nclstot =  nclstot +1
            enddo
            if(nclstot.ne.0)print*,'iev ',iev,' ncls ',nclstot
c            if(nclstot.ge.1)goto 100 
*     ==================================
*     //////////////////////////////////
*     Track selection (end)
*     \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
*     ==================================

            call HCDIR('//LEVEL2',' ')
                        
            goodtr = goodtr +1
            nselect = nselect +1
c            print*,goodtr,' ****',iev
            call HFNT(ntp_level2) !fill LEVEL2 nt-uple         
            

            call HCDIR('//PAWC','')

*     -----------------------------------------------------
*     +++++++++++++++++++++++++++++++++++++++++++++++++++++
*     calibrations (sum over events)
*     +++++++++++++++++++++++++++++++++++++++++++++++++++++
*     -----------------------------------------------------

            if(goodtr.lt.ncorrmax)then
               do ip=1,nplanes  !1-->6
                  
                  iviewx = (nplanes - ip)*2+2 !12-->2
                  ilx=nld(maxs(which_cl_x(ip)),iviewx)
                  iviewy = (nplanes - ip)*2+1 !11-->1
                  ily=nld(maxs(which_cl_y(ip)),iviewy)
                  if(ilx.ne.ily)print*,'LADDERS DO NOT MATCH!!!'
                                    
                  chargexy(ip,ilx,goodtr,1)=dedx(which_cl_x(ip))
                  chargexy(ip,ily,goodtr,2)=dedx(which_cl_y(ip))
                  
               enddo
            endif
*     compute average amgle in each bin
            do i=1,nviews
               do ia=1,nangbin
                  if(navangbin(ia,i).ne.0)
     $                 avangbin(ia,i)=
     $                 avangbin(ia,i)/navangbin(ia,i)
               enddo
            enddo

*     apply a cut on the measured rigidity:
            if(rig.gt.rigmax)then
               call fill_eta_histos
               call fill_check_histos
            endif
            
*     -----------------------------------------------------
*     +++++++++++++++++++++++++++++++++++++++++++++++++++++
*     calibrations (end)
*     +++++++++++++++++++++++++++++++++++++++++++++++++++++
*     -----------------------------------------------------
 100        continue

         enddo                  !end loop on events             

         
         print*,''
         print *,'- Stored ',nselect,' entries in LEVEL2 nt-uple'
         call HCDIR('//LEVEL1',' ')
         call HREND('level1')
         close(lun_data_level1) 

         call HCDIR('//PAWC',' ')
         call HDELET(ntp_level1+20)

      enddo                     !end loop on files

      print*,''
      call HCDIR('//LEVEL2',' ')
      call RZPURG(-1)    
      call HROUT(ntp_level2,ICYCLE,'T')
      call HREND('level2')
      close(lun_data_level2)

      call HCDIR('//PAWC',' ')
      call HDELET(ntp_level2)


*     ===================================================
*     in the following the charge correlation parameters
*     and the eta correction functions are evaluated
*     and stored in txt files.
*     However I chosed to do it "off-line", by means of 
*     a kumac, so that I can group more files together.
*     ===================================================

*     -----------------------------------------------------
*     +++++++++++++++++++++++++++++++++++++++++++++++++++++
*     calibrations
*     +++++++++++++++++++++++++++++++++++++++++++++++++++++
*     -----------------------------------------------------

      do il=1,3
         write(rzdir,1010)il
**************************************************
**************************************************
*     CHARGE CORRELATION CALIBRATION
**************************************************
**************************************************
         if(DONTDO)goto 8881
         call charge_calibration(il)
 107     format(/
     $        ,/
     $        ,/
     $        ' >>> CHARGE-CORRELATION CALIBRATION <<<',/
     $        ,/
     $        'Save charge correlation parameters --->',/
     $        ,a,/
     $        ,/
     $        ,'Parameterization ------>  Sy = K * Sx + C',/
     $        ,'Plane        K       EK        C       EC',/
     $        ,'-----------------------------------------')
c     write(lun_file_log,107)file_out(1:LNBLNK(file_out))
         
         file_out=
     $        out_dir(1:LNBLNK(out_dir))//'trk-'
     $        //rzdir//'-charge.dat'
         write(*,107)file_out(1:LNBLNK(file_out))
         open(55,
     $        FILE=file_out(1:LNBLNK(file_out)),
c     $     STATUS='NEW',
     $        FORM='FORMATTED')
         
 108     format(i5,4f9.5)
         do ip=1,nplanes
c            write(lun_file_log,108) !log file
c     $           ip,syksx(ip),esyksx(ip),sycsx(ip),esycsx(ip)
            write(*,108)        !screen
     $           ip,syksx(ip),esyksx(ip),sycsx(ip),esycsx(ip)
            write(55,*)         !data
     $           ip,syksx(ip),esyksx(ip),sycsx(ip),esycsx(ip) 
         enddo
         close(55)
         
 8881    call fill_charge_histos(il)
**************************************************
**************************************************
*     ETA CALIBRATION
**************************************************
**************************************************
         file_out=
     $        out_dir(1:LNBLNK(out_dir))//'trk-'
     $        //rzdir
         
 109     format('Creating files:',/
     $        ,' ',a,'-eta2-binning.dat',/
     $        )
         write(*,109)file_out(1:LNBLNK(file_out))
         open(lun_file_eta,
     $        FILE=file_out(1:LNBLNK(file_out))//'-eta2-binning.dat'
     $        ,FORM='FORMATTED'
     $        )
         do ia=1,nangbin
            write(lun_file_eta,*)ia,angL(ia),angR(ia)
         enddo
         close(lun_file_eta)
      

         if(DONTDO)goto 8882
         do ia=1,nang

 201        format('-eta2-bin',i1,'.dat')
 202        format('-eta2-bin',i2,'.dat')
            if(ia.lt.10)write(fname_param,201)ia
            if(ia.ge.10)write(fname_param,202)ia
            
            open(lun_file_eta,
     $           FILE=file_out(1:LNBLNK(file_out))
     $           //fname_param(1:LNBLNK(fname_param))
     $           ,FORM='FORMATTED'
     $           )
            
            call eta_calibration(il,ia)
            
*     ********************** write **************************
 203        format(13f10.4)
            do ie=1,nhbin2+1
               write(lun_file_eta,203)
     $              eta2abs(ie),(eta2corr(ie,i),i=1,nviews)
            enddo   
            close(lun_file_eta)
*     ********************** write **************************
            
         enddo                  !end loop on angular bins
 8882    continue
         
      enddo                     !end loop on LADDERs

*     -----------------------------------------------------
*     +++++++++++++++++++++++++++++++++++++++++++++++++++++
*     calibrations
*     +++++++++++++++++++++++++++++++++++++++++++++++++++++
*     -----------------------------------------------------
 8888 continue


c------------------------------------------------------------------------
c     
c     no error exit
c     
c------------------------------------------------------------------------

      goto 9000                 !happy ending


c------------------------------------------------------------------------
c     
c     data file opening error
c     
c------------------------------------------------------------------------
 19   continue
      
      print*,' '
      print*,'ERROR OPENING DATA FILE: ',data_file
      print*,' '
      print*,' '
      
      goto 9000                 !the end

c------------------------------------------------------------------------
c     
c     level1 ntuple event reading error
c     
c------------------------------------------------------------------------
      
 21   continue
      
      print*,' '
      print*,'ERROR WHILE READING LEVEL1 NTUPLE, AT EVENT
     $     : ',iev
      print*,' '
      print*,' '
      
      goto 9000                 !the end

c------------------------------------------------------------------------
c     
c     calib file opening error
c     
c------------------------------------------------------------------------
      
 17   continue
      
      print*,' '
      print*,'preanalysis: ERROR OPENING INPUT FILE: ',data_file_calib
      print*,' '
      print*,' '
      
      goto 9000                 !the end

 16   continue
      
      print*,' '
      print*,'preanalysis: ERROR OPENING OUTPUT FILE: ',file_out
      print*,' '
      print*,' '
      
      goto 9000                 !the end

c------------------------------------------------------------------------
c     
c     closes files and exits
c     
c------------------------------------------------------------------------

 9000 continue
      
c$$$      if(DEBUG)call HPRNTU(ntp_level2+1)
c$$$      call HPRNTU(ntp_level2)
c$$$      print*,''
c$$$      call HCDIR('//LEVEL2',' ')
c$$$      call HROUT(ntp_level2,ICYCLE,'T')
c$$$      print *,'- Stored LEVEL2 nt-uple (',nev2,' entries )'
c$$$      if(DEBUG)call HROUT(ntp_level2+1,ICYCLE,'T')
c$$$      if(DEBUG)print *,'- Stored DEBUG nt-uple '
c$$$      call HREND('level2')
c$$$      close(lun_data_level2)

c$$$      call HCDIR('//LEVEL1',' ')
c$$$      call HREND('level1')
c$$$      close(lun_data_level1)

      write(*,105)
     $     tottr
     $     ,goodtr
 105  format(/
     $     ,/
     $     ,'|------------------------------------------------|',/
     $     ,'|  Number of processed events                    |',/
     $     ,'|  Total         ',i8,'                        |'/
     $     ,'|  Good tracks   ',i8,'                        |'/
     $     ,'|------------------------------------------------|',/
     $     )
      
      
      
      do ilad=1,3
         write(rzdir,1010)ilad
         
         call HCDIR('//PAWC/'//rzdir,' ')
         call HCDIR('//'//RZDIR,' ')
         
         call HROUT(id_hi_chi,ICYCLE,'T')
         call HROUT(id_hi_chi+100,ICYCLE,'T')
         call HROUT(id_hi_chi+200,ICYCLE,'T')
         call HROUT(id_hi_rig,ICYCLE,'T')
         call HROUT(id_hi_charge,ICYCLE,'T')
         do i=1,nviews
            call HROUT(id_hi_charge+i*100,ICYCLE,'T')
            call HROUT(id_hi_charge+i*100+1,ICYCLE,'T')
            call HROUT(id_hi_charge+i*100+2,ICYCLE,'T')
            call HROUT(80000+id_hi_charge+i*100+2,ICYCLE,'T')
            call HROUT(id_hi_mult+i*100,ICYCLE,'T')
            call HROUT(80000+id_hi_mult+i*100,ICYCLE,'T')
            do ia=1,nang
               call HROUT(id_hi_eta2+i*100+ia,ICYCLE,'T')
               call HROUT(id_hi_eta3+i*100+ia,ICYCLE,'T')
               call HROUT(id_hi_eta4+i*100+ia,ICYCLE,'T')
c               call HROUT(id_hi_charge+i*100+ia,ICYCLE,'T')
               call HROUT(id_hi_mult+i*100+ia,ICYCLE,'T')
               call HROUT(id_hi_mult+i*100+ia+50000,ICYCLE,'T')
            enddo            
c     call HROUT(id_hi_chi+100+ia,ICYCLE,'T')
c     call HROUT(id_hi_chi+200+ia,ICYCLE,'T')
         enddo
         do ip=1,nplanes
            call HROUT(id_hi_angle+100+ip,ICYCLE,'T')
            call HROUT(id_hi_angle+200+ip,ICYCLE,'T')
            call HROUT(id_hi_angle_stat+100+ip,ICYCLE,'T')
            call HROUT(id_hi_angle_stat+200+ip,ICYCLE,'T')
            call HROUT(id_hi_chargeco+ip*100+4,ICYCLE,'T')
            call HROUT(id_hi_chargeco+ip*100,ICYCLE,'T')
            do is=1,2
               call HROUT(id_hi_beam+ip*100+is+200000,ICYCLE,'T')
               call HROUT(id_hi_beam+ip*100+is+100000,ICYCLE,'T')
            enddo
            call HROUT(id_hi_bx+ip,ICYCLE,'T')
            call HROUT(id_hi_by+ip,ICYCLE,'T')
            call HROUT(id_hi_bz+ip,ICYCLE,'T')
         enddo
         
         call HREND(RZDIR)
         close(lun_file_out+ilad)
         
      enddo
      
     




      stop
      end


# include "momanhough-subroutines.F"


****** * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ******
******* * * * * * * * * * * * * * * * * * * * * * * * * * * * * *******
****** * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ******
******* * * * * * * * * * * * * * * * * * * * * * * * * * * * * *******
****** * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ******
******* * * * * * * * * * * * * * * * * * * * * * * * * * * * * *******
****** * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ******
******* * * * * * * * * * * * * * * * * * * * * * * * * * * * * *******
****** * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ******
******* * * * * * * * * * * * * * * * * * * * * * * * * * * * * *******
****** * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ******
******* * * * * * * * * * * * * * * * * * * * * * * * * * * * * *******

***************************************************************
*
*
*
*
*
***************************************************************
      subroutine book_eta_histos

      include '../common/commontracker.f'
      include '../common/calib.f'
      include '../common/common_calibration.f'  

      character*64 title        !histos title

 502  format('Eta-2 (view ',i2,', angle bin ',i2,' )')
 503  format('Eta-3 (view ',i2,', angle bin ',i2,' )')
 504  format('Eta-4 (view ',i2,', angle bin ',i2,' )')
      
      do i=1,nviews
         do ia=1,nangbin
            write(title,502)i,ia
            call HBOOK1(id_hi_eta2+i*100+ia
     $           ,title
c     $           ,nhbin2,0.,1.,0.)
     $           ,nhbin2,-0.5,0.5,0.)
            write(title,503)i,ia
            call HBOOK1(id_hi_eta3+i*100+ia
     $           ,title
c     $           ,nhbin3,0.,2.,0.)
     $           ,nhbin2,-0.5,0.5,0.)
            write(title,504)i,ia
            call HBOOK1(id_hi_eta4+i*100+ia
     $           ,title
c     $           ,nhbin4,0.,3.,0.)
c     $           ,nhbin3,-1.,1.,0.)
     $           ,nhbin3,-0.5,.5,0.)
         enddo            
      enddo


      return
      end
C............................................................
      subroutine fill_eta_histos

      include '../common/commontracker.f'
      include '../common/calib.f'
      include '../common/common_calibration.f'
      include '../common/level1.f'
      include '../common/level2.f'

      character*10 rzdir

*     flag to tag low-multilpicity cluster
*     (according to Samu definition)
      logical LMSX,LMSY,DRAY


 10   format('LADDER',i1)
*     start loop on planes to fill histos
      do ip=1,nplanes           !1-->6

*---------------------------------------------------------
*     main x-y cluster parameters
*---------------------------------------------------------
         xxx = xm_nt(ip,1)
         yyy = ym_nt(ip,1)
         call whichsensor(ip,xxx,yyy,iladder,isensor)
         if(iladder.eq.0.or.isensor.eq.0)then            
            goto 10102
         endif

         multx = mult(which_cl_x(ip))
         multy = mult(which_cl_y(ip))
         thetax = axv_nt(ip,1) 
         thetay = ayv_nt(ip,1) 
         LMSX=.false.
         LMSY=.false.
*     ===============================================
*     multiplicity selection of SAMUELE
c$$$         if(
c$$$     $        (thetax.le.15..and.multx.le.2).or.
c$$$     $        (thetax.gt.15..and.multx.le.3).or.
c$$$     $        .false.)LMSX=.true.
c$$$         if(
c$$$     $        (multx.le.2).or.
c$$$     $        .false.)LMSY=.true.
*     ===============================================
*     looser cut ... after LANDI warning
*     (I also require that both clx anc cly are low-multiplicity)
         if(
     $   (abs(thetax).le.10..and.multx.le.3).or.
     $  ((abs(thetax).gt.10..and.abs(thetax).le.15.).and.multx.le.4).or.
     $   (abs(thetax).gt.15..and.multx.le.5).or.
     $        .false.)LMSX=.true.
         if(
     $        (multy.le.3).or.
     $        .false.)LMSY=.true.


*     charge scaling factor (to vertical trajectory)
         pi=acos(-1.)
         fact=
     $        sqrt( 
     $        tan(pi*thetax/180.)**2 + 
     $        tan(pi*thetay/180.)**2 + 
     $        1. 
     $        )
         dedxscalx = dedx(which_cl_x(ip)) / fact
         dedxscaly = dedx(which_cl_y(ip)) / fact
         dedxscal = (dedxscalx+dedxscaly)/2
 
         DRAY=.false.
         if(dedxscal.gt.dedxmax)DRAY=.true.

*---------------------------------------------------------


         write(rzdir,10)iladder
         call HCDIR('//'//rzdir,'')
         call HCDIR('//PAWC/'//rzdir,'')

*---------------------------------------------------------
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
*                      X VIEWS
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
*---------------------------------------------------------

         if(  thetax.lt.angL(1).or.
     $        thetax.ge.angR(nangbin).or.
     $        .not.LMSX.or.
     $        .not.LMSY.or.
     $        DRAY.or.
     $        .false.)goto 10101 ! skip
         
         iview = (nplanes - ip)*2+2 !12-->2
         e2=cog(2,which_cl_x(ip))
         e3=cog(3,which_cl_x(ip))
         e4=cog(4,which_cl_x(ip))
         do ia=1,nangbin            
            if (thetax.lt.angR(ia))goto 111
         enddo
 111     ibin_ang = ia          !angular bin

         avangbin(ibin_ang,iview)=avangbin(ibin_ang,iview)+thetax
         navangbin(ibin_ang,iview)=navangbin(ibin_ang,iview)+1

*!!!!!!!!NB non sempre MAXS coincide con la strip massima !!!!
         if(e2.ge.0.5)e2=e2-1.
         if(e2.lt.-0.5)e2=e2+1.
         if(e3.ge.0.5)e3=e3-1.
         if(e3.lt.-0.5)e3=e3+1.
c         if(e4.ge.1.)e4=e4-2.
c         if(e4.lt.-1.)e4=e4+2.
         if(e4.ge.0.5)e4=e4-1.
         if(e4.lt.-0.5)e4=e4+1.
*!!!!!!!!NB

         if(e2.ne.0.)
     $        call hfill(id_hi_eta2+iview*100+ibin_ang,
     $        e2,0.,1.)
         if(e3.ne.0.)
     $        call hfill(id_hi_eta3+iview*100+ibin_ang,
     $        e3,0.,1.)
         if(e4.ne.0.)
     $        call hfill(id_hi_eta4+iview*100+ibin_ang,
     $        e4,0.,1.)

10101    continue
         
*---------------------------------------------------------
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
*                      Y VIEWS
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
*---------------------------------------------------------
         if(  thetay.lt.angL(1).or.
     $        thetay.ge.angR(nangbin).or.
     $        .not.LMSY.or.
     $        .not.LMSX.or.
     $        DRAY.or.
     $        .false.)goto 10102 !skip
         
         iview = (nplanes - ip)*2+1 !11-->1
         e2=cog(2,which_cl_y(ip))
         e3=cog(3,which_cl_y(ip))
         e4=cog(4,which_cl_y(ip))
         do ia=1,nangbin            
c            if (thetay.lt.(ang(ia)+angstep/2.))goto 112
c            if (thetay.lt.angbin(ia))goto 112
            if (thetay.lt.angR(ia))goto 112
         enddo
 112     ibin_ang = ia

         avangbin(ibin_ang,iview)=avangbin(ibin_ang,iview)+thetay
         navangbin(ibin_ang,iview)=navangbin(ibin_ang,iview)+1
c         print*,'****',avangbin(ibin_ang,iview),thetay


*!!!!!!!!NB non sempre MAXS coincide con la strip massima !!!!
         if(e2.gt.0.5)e2=e2-1.
         if(e2.lt.-0.5)e2=e2+1.
         if(e3.ge.0.5)e3=e3-1.
         if(e3.lt.-0.5)e3=e3+1.
c         if(e4.ge.1.)e4=e4-2.
c         if(e4.lt.-1.)e4=e4+2.
         if(e4.ge.0.5)e4=e4-1.
         if(e4.lt.-0.5)e4=e4+1.
*!!!!!!!!NB

         if(e2.ne.0.)
     $        call hfill(id_hi_eta2+iview*100+ibin_ang,
     $        e2,0.,1.)
         if(e3.ne.0.)
     $        call hfill(id_hi_eta3+iview*100+ibin_ang,
     $        e3,0.,1.)
         if(e4.ne.0.)
     $        call hfill(id_hi_eta4+iview*100+ibin_ang,
     $        e4,0.,1.)

10102    continue
      enddo
      call HCDIR('//PAWC','')
      return

      end

C............................................................
      subroutine book_charge_histos

      include '../common/commontracker.f'
      include '../common/calib.f'  
      include '../common/common_calibration.f'

      character*64 title        !histos titleinfo -f g77 M LEX

c 500  format('Charge correlation (ladder ',i1,' plane ',i2,')')
 500  format('Charge correlation ( plane ',i2,')')
 501  format('Charge residuals ( plane',i2,')')
      
      do ip=1,nplanes
         write(title,500)ip
         call HBOOK2(id_hi_chargeco+ip*100+4
     $           ,title
     $           ,200,0.,1000.
     $        ,200,0.,1000.,0.)
         write(title,501)ip
         call HBOOK1(id_hi_chargeco+ip*100
     $        ,title
     $        ,200,-100.,100.,0.)
      enddo



      return
      end
C............................................................
      subroutine fill_charge_histos(iladder)

      include '../common/commontracker.f'
      include '../common/calib.f'  
      include '../common/common_calibration.f'
      include '../common/level1.f'
      include '../common/level2.f'
*     start loop on planes

c      parameter (ncorrmax=20000)
c      real chargexy(nplanes,ncorrmax,2)
c      integer whichplane,goodtr
c      common/chargecorr/chargexy,whichplane,goodtr

      character*10 rzdir

 1010 format('LADDER',i1)
      write(rzdir,1010)iladder
      call HCDIR('//'//rzdir,'')
      call HCDIR('//PAWC/'//rzdir,'')
      
      do ip=1,nplanes           !1-->6
         
         do ie=1,min(goodtr,ncorrmax)
            
            sx = chargexy(ip,iladder,ie,1)
            sy = chargexy(ip,iladder,ie,2)
            
            if(sx.ne.0..and.sy.ne.0)then

               call HFILL(id_hi_chargeco+ip*100+4,sx,sy,1.) 
               d= sy - syksx(ip)*sx - sycsx(ip)
               d=d/sqrt(syksx(ip)**2+1)          
               call HFILL(id_hi_chargeco+ip*100,
     $              d,0.,1.)
            
            endif
            
         enddo
      enddo

      return
      end
C............................................................

      subroutine book_check_histos
*
*     some histos to check the sample tracks
*     and define quality cuts
*
      include '../common/commontracker.f'
      include '../common/calib.f'
      include '../common/common_calibration.f'  

      character*64 title        !histos title
      real xbins(nangbin+1)

      parameter (nrigbin=500)
      real xrigbin(nrigbin+1)

* * * * * * * * * * * * * * * * * * * * * * * *
*     CHI**2 + RIGIDITY


      rmin=0.2
      rmax=1000.
      step=abs(log10(rmin)-log10(rmax))/nrigbin
      xrigbin(1)=rmin
      do i=2,nrigbin+1
         xrigbin(i)=10**(log10(xrigbin(i-1))+step)
*         print*,'### ',i,xrigbin(i-1),xrigbin(i)
      enddo
      call HBOOKB(id_hi_rig
     $     ,'Rigidity'
     $     ,nrigbin,xrigbin,0.)


      call HBOOK1(id_hi_chi
     $     ,'Chi**2'
     $     ,1000,0.,100.,0.)
      call HBPROF(id_hi_chi+100
     $     ,'Chi**2 vs Rigidity'
     $     ,150,rmin,0.1*rmax,0.,100.,0.)
      call HBPROF(id_hi_chi+200
     $     ,'Chi**2 vs Angle x'
     $     ,50,0.,25.,0.,100.,0.)

      
*     * * * * * * * * * * * * * * * * * * * * * * *
*     ANGLE

      nbin=200.
      from=angL(1)!angbin(1)
      to=angR(nangbin)!angbin(nang+1)   
      do ip=1,nplanes
         call HBOOK1(id_hi_angle+100+ip
     $        ,'Angle (Y view)'
     $        ,nbin,from,to,0.)
         call HBOOK1(id_hi_angle+200+ip
     $        ,'Angle (X view)'
     $        ,nbin,from,to,0.)
         
c$$$         call HBOOK1(id_hi_angle_stat+100+ip
c$$$     $        ,'Angle (Y view)'
c$$$     $        ,nangbin,from,to,0.)
c$$$         call HBOOK1(id_hi_angle_stat+200+ip
c$$$     $        ,'Angle (X view)'
c$$$     $        ,nangbin,from,to,0.)

         xbins(1)=angL(1)
         do i=1,nangbin
            xbins(i+1)=angR(i) 
         enddo
         call HBOOKB(id_hi_angle_stat+100+ip
     $        ,'Angle (Y view)'
     $        ,nangbin,XBINS,0.)
         call HBOOKB(id_hi_angle_stat+200+ip
     $        ,'Angle (X view)'
     $        ,nangbin,XBINS,0.)
      enddo


*     * * * * * * * * * * * * * * * * * * * * * * *
*     MULTIPLICITY
 511  format('Multiplicity (view ',i2,', angle bin ',i2,' )')
 512  format('Multiplicity vs angle (view ',i2,' )')
      
      do i=1,nviews
         do ia=1,nangbin
            write(title,511)i,ia
            call HBOOK1(id_hi_mult+i*100+ia
     $           ,title
     $           ,10,0.5,10.5,0.)
         enddo            
         write(title,512)i
         call HBOOK2(id_hi_mult+i*100
     $        ,title
     $        ,50,0.,25.,10,0.5,10.5,0.)
         call HBOOK2(80000+id_hi_mult+i*100
     $        ,title
     $        ,50,0.,25.,10,0.5,10.5,0.)
      enddo
*     * * * * * * * * * * * * * * * * * * * * * * *
*     AVERAGE CLUSTER
 510  format('Average cluster (view ',i2,', angle bin ',i2,' )')

      do i=1,nviews
         do ia=1,nangbin
            write(title,510)i,ia
            call HBOOK1(id_hi_mult+i*100+ia+50000
     $           ,title
     $           ,int(2*nhside+1)
     $           ,-(float(nhside)+0.5),(float(nhside)+0.5),0.)
         enddo            
      enddo


* * * * * * * * * * * * * * * * * * * * * * * *
*     CHARGE
      
 5111 format('dE/dX scaled to 300 micron (view ',i2,' )')
 5112 format('dE/dX vs R (view ',i2,' )')
 5113 format('dE/dX vs multiplicity (view ',i2,' )')
      do i=1,nviews
         write(title,5111)i
         call HBOOK1(id_hi_charge+i*100
     $        ,title
     $        ,500,0.,1000.,0.)
         write(title,5112)i
         call HBPROF(id_hi_charge+i*100+1
     $        ,title
     $        ,1000,0.,100.,0.,1000.,0.)
         write(title,5113)i
         call HBOOK2(id_hi_charge+i*100+2
     $        ,title
     $        ,100,0.,1000.,10,0.5,10.5,0.)
         call HBOOK2(80000+id_hi_charge+i*100+2
     $        ,title
     $        ,100,0.,1000.,10,0.5,10.5,0.)
      enddo
      call HBOOK1(id_hi_charge
     $     ,'dE/dX scaled to 300 micron AVERAGE'
     $     ,500,0.,1000.,0.)

* * * * * * * * * * * * * * * * * * * * * * * *
*     EVENT SPATIAL DISTRIBUTION

51111 format('Strip Y (plane ',i2,', sensor ',i1,')')
51112 format('Strip X (plane ',i2,', sensor ',i1,')')
      do ip=1,nplanes
         do is=1,2
            write(title,51111)ip,is
            call HBOOK1(id_hi_beam+ip*100+is+100000
     $           ,title
     $           ,256,0.5,3072.5,0.)
            write(title,51112)ip,is
            call HBOOK1(id_hi_beam+ip*100+is+200000
     $           ,title
     $           ,256,0.5,3072.5,0.)
         enddo
      enddo

* * * * * * * * * * * * * * * * * * * * * * * *
*     MAGNETIC FIELD

 61   format('Bx (plane ',i2,')')
 62   format('By (plane ',i2,')')
 63   format('Bz (plane ',i2,')')
      do ip=1,nplanes
         write(title,61)ip
         call HBOOK1(id_hi_bx+ip
     $        ,title
     $        ,1000,-5.,5.,0.)
         write(title,62)ip
         call HBOOK1(id_hi_by+ip
     $        ,title
     $        ,1000,-5.,5.,0.)
         write(title,63)ip
         call HBOOK1(id_hi_bz+ip
     $        ,title
     $        ,1000,-5.,5.,0.)
      enddo



      return
      end
C............................................................
      subroutine fill_check_histos

      include '../common/commontracker.f'
      include '../common/calib.f'  
      include '../common/common_calibration.f'
      include '../common/common_mech.f'
      include '../common/level1.f'
      include '../common/level2.f'

      real v(3),b(3)
      character*10 rzdir

      real clusterx(2*nhside+1)
      real clustery(2*nhside+1)
      logical DRAY,LMSX,LMSY

*---------------------------------------------------------
*     track info parameters
*---------------------------------------------------------
      rig = abs(1./al_nt(5,1))
      chi = chi2_nt(1)



 1010 format('LADDER',i1)

*     i plot della rigidita` del chi2
*     non dipendono dal ladder...
      do il=1,3
         write(rzdir,1010)il
         call HCDIR('//'//rzdir,'')
         call HCDIR('//PAWC/'//rzdir,'')
         call HFILL(id_hi_chi,chi,0.,1.)
         call HFILL(id_hi_rig,rig,0.,1.)
         thetax = axv_nt(1,1) 
         thetay = ayv_nt(1,1) 
         call HFILL(id_hi_chi+100,chi,rig,1.)
         call HFILL(id_hi_chi+200,chi,thetax,1.)

      enddo
      
      dedxscalav   = 0
      npdedxscalav = 0
*     start loop on planes
      do ip=1,nplanes           !TOP-->BOTTOM
*---------------------------------------------------------
*     main x-y cluster parameters
*---------------------------------------------------------
         iviewx = nviewx(7-ip)
         iviewy = nviewy(7-ip)
         xxx = xm_nt(ip,1)
         yyy = ym_nt(ip,1)
         call whichsensor((7-ip),xxx,yyy,iladder,isensor)
         if(iladder.eq.0.or.isensor.eq.0)then            
            print*,'*** out *** ',ip,xxx,yyy,iladder,isensor
            goto 10101
         endif
         multx = mult(which_cl_x(ip))
         multy = mult(which_cl_y(ip))
         thetax = axv_nt(ip,1) 
         thetay = ayv_nt(ip,1) 
*     ===============================================
*     looser multiplicity cut ... after LANDI warning
*     (I also require that both clx anc cly are low-multiplicity)
         LMSX=.false.
         LMSY=.false.
         if(
     $   (abs(thetax).le.10..and.multx.le.3).or.
     $  ((abs(thetax).gt.10..and.abs(thetax).le.15.).and.multx.le.4).or.
     $   (abs(thetax).gt.15..and.multx.le.5).or.
     $        .false.)LMSX=.true.
         if(
     $        (multy.le.3).or.
     $        .false.)LMSY=.true.

*     charge scaling factor (to vertical trajectory)
         pi=acos(-1.)
         fact=
     $        sqrt( 
     $        tan(pi*thetax/180.)**2 + 
     $        tan(pi*thetay/180.)**2 + 
     $        1. 
     $        )
         dedxscalx = dedx(which_cl_x(ip)) / fact
         dedxscaly = dedx(which_cl_y(ip)) / fact         
         dedxscal = (dedxscalx+dedxscaly)/2     
         DRAY=.false.
         if(dedxscal.gt.dedxmax)DRAY=.true.
    
         if(.not.DRAY.and.LMSX.and.LMSX)then
            dedxscalav = dedxscalav + dedxscal
            npdedxscalav = npdedxscalav + 1
         endif



         ista  = indstart(which_cl_x(ip))
         ice   = indmax(which_cl_x(ip))
         ilast = totCLlength
         if(which_cl_x(ip).ne.nclstr1)
     $        ilast = indstart(which_cl_x(ip)+1)-1
         do i = -nhside,nhside
            ii = ice + i
            iii=i+nhside+1
            clusterx(iii) = 0
            if(ii.le.ilast.and.ii.ge.ista)
     $           clusterx(iii) = clsignal(ii) / dedx(which_cl_x(ip))
         enddo

         ista  = indstart(which_cl_y(ip))
         ice   = indmax(which_cl_y(ip))
         ilast = totCLlength
         if(which_cl_y(ip).ne.nclstr1)
     $        ilast = indstart(which_cl_y(ip)+1)-1
         do i = -nhside,nhside
            ii = ice + i
            iii=i+nhside+1
            clustery(iii) = 0
            if(ii.le.ilast.and.ii.ge.ista)
     $           clustery(iii) = clsignal(ii) / dedx(which_cl_y(ip))
         enddo

*---------------------------------------------------------

         write(rzdir,1010)iladder
         call HCDIR('//'//rzdir,'')
         call HCDIR('//PAWC/'//rzdir,'')


 
*---------------------------------------------------------
*     find the angular bins (ibin_angy ibin_angx)  
*     associated to the track intersection point with plane ip
*---------------------------------------------------------
         ia=0
         if(  thetax.lt.angL(1).or.
     $        thetax.ge.angR(nangbin))goto 111         
         do ia=1,nangbin            
            if (thetay.lt.angR(ia))goto 111
         enddo
 111     ibin_angy = ia

         ia=0
         if(  thetay.lt.angL(1).or.
     $        thetay.ge.angR(nangbin))goto 112
         do ia=1,nangbin            
            if (thetax.lt.angR(ia))goto 112
         enddo
 112     ibin_angx = ia



*---------------------------------------------------------
*     fill histos
*---------------------------------------------------------

************************************************ ANGLE      

         call HFILL(id_hi_angle+200+ip,thetax,0.,1.)
         call HFILL(id_hi_angle+100+ip,thetay,0.,1.)
      
         call HFILL(id_hi_angle_stat+200+ip,thetax,0.,1.)
         call HFILL(id_hi_angle_stat+100+ip,thetay,0.,1.)
         
************************************************ multiplicity & cluster
         if(ibin_angx.ne.0)then
            call HFILL(id_hi_mult+iviewx*100+ibin_angx,
     $           float(mult(which_cl_x(ip))),0.,1.)
            do i=-nhside,nhside
               ii = i + nhside +1
               xxxx=clusterx(ii)
               call HFILL(id_hi_mult+iviewx*100+ibin_angx+50000,
     $              float(i),0.,xxxx)
            enddo
         endif

         if(ibin_angy.ne.0)then
            call HFILL(id_hi_mult+iviewy*100+ibin_angy,
     $           float(mult(which_cl_y(ip))),0.,1.)
            do i=-nhside,nhside
               ii = i + nhside +1
               xxxx=clustery(ii)
               call HFILL(id_hi_mult+iviewy*100+ibin_angy+50000,
     $              float(i),0.,xxxx)
            enddo
         endif

************************************************ charge
         call HFILL(id_hi_charge+iviewx*100,
     $        dedxscalx,0.,1.)         
         call HFILL(id_hi_charge+iviewy*100,
     $        dedxscaly,0.,1.)
         
         call HFILL(id_hi_charge+iviewx*100+1,
     $        rig,dedxscalx,1.)         
         call HFILL(id_hi_charge+iviewy*100+1,
     $        rig,dedxscaly,1.)

*     --------------------------------------------
C     this are the condition applied to select the
C     particle sample for pfa calibration
*     --------------------------------------------
         if(.not.DRAY.and.LMSX.and.LMSY)then
            call HFILL(id_hi_charge+iviewx*100+2,
     $           dedxscalx,float(multx),1.)         
            call HFILL(id_hi_charge+iviewy*100+2,
     $           dedxscaly,float(multy),1.)
            call HFILL(id_hi_mult+iviewx*100,
     $           abs(thetax),float(multx),1.)
            call HFILL(id_hi_mult+iviewy*100,
     $           abs(thetay),float(multy),1.)
         endif
         call HFILL(80000+id_hi_charge+iviewx*100+2,
     $        dedxscalx,float(multx),1.)         
         call HFILL(80000+id_hi_charge+iviewy*100+2,
     $        dedxscaly,float(multy),1.)
         call HFILL(80000+id_hi_mult+iviewx*100,
     $        abs(thetax),float(multx),1.)
         call HFILL(80000+id_hi_mult+iviewy*100,
     $        abs(thetay),float(multy),1.)
         
************************************************ beam profle
         call hfill(id_hi_beam+ip*100+isensor+100000,
     $        float(MAXS(which_cl_y(ip))),0.,1.)
         call hfill(id_hi_beam+ip*100+isensor+200000,
     $        float(MAXS(which_cl_x(ip))),0.,1.)


************************************************ magnetic field

         v(1) = xv_nt(ip,1) 
         v(2) = yv_nt(ip,1) 
         v(3) = zv_nt(ip,1) 

         call gufld(v,b)

         bx=b(1)
         by=b(2)  
         bz=b(3)
         call hfill(id_hi_bx+ip,bx,0.,1.)
         call hfill(id_hi_by+ip,by,0.,1.)
         call hfill(id_hi_bz+ip,bz,0.,1.)

      enddo

      dedxscalav = dedxscalav/npdedxscalav
      call HFILL(id_hi_charge,
     $     dedxscalav,0.,1.)



10101 continue
      call HCDIR('//PAWC','')

      return
      end

***************************************************************************
*
*
*
*
*
*
*
*
*
***************************************************************************


      subroutine FCNcharge(npar,grad,fval,xval,iflag,futil)

      include '../common/commontracker.f'
      include '../common/calib.f'
      include '../common/common_calibration.f'

      double precision grad(*),xval(*),fval
      integer npar,iflag
      double precision KKK,CCC


c      parameter (ncorrmax=20000)
c      real chargexy(nplanes,ncorrmax,2)
c      integer whichplane,goodtr
c      common/chargecorr/chargexy,whichplane,goodtr
      
      KKK=xval(1)
      CCC=xval(2)


      fval = 0.


      sum2=0
      do i=1,min(goodtr,ncorrmax)
         sx=chargexy(whichplane,whichladder,i,1)
         sy=chargexy(whichplane,whichladder,i,2)
*         if(sx.gt.140..and.sy.gt.140)then
            d=(sy-KKK*sx-CCC)**2
            d=d/(KKK**2+1)
            d=d/(1.9*4.**2+1.2*8.**2) !errors
            sum2=sum2+d
c$$$  PRINT*,'xx ',i,d,' sx-sy  '
c$$$     $        ,chargexy(whichplane,whichladder,i,1)
c$$$     $        ,chargexy(whichplane,whichladder,i,2)
*         endif
      enddo


      fval = sum2



      return

      end




**************************************************
**************************************************
*     CHARGE CORRELATION CALIBRATION
**************************************************
**************************************************
      subroutine charge_calibration(iladder)

      include '../common/commontracker.f'
      include '../common/calib.f'  
      include '../common/common_calibration.f'
c      include '../common/common_mech.f'
c      include '../common/level1.f'
c      include '../common/level2.f'
      external FCNcharge

      character*12 chnam
      double precision val,error
      

      do ip=1,nplanes
         print*,' '
         print*,' '
         print*,' '
         print*,
     $        '               -----------------------------------'
         print*,
     $        '////////////// MINUIT for plane ',ip,' - ladder ',iladder
         print*,
     $        '               -----------------------------------'
         open(UNIT=40,FILE=
     $        'bin-aux/data-card.dat'
     $        ,STATUS='OLD')
         call mintio(40,6,7)
         whichplane = ip
         whichladder = iladder
         call minuit(FCNcharge,0)
c     call MNERRS(1,eplus,eminus,eparab,globcc)
         call MNPOUT(1,chnam,val,error,bnd1,bnd2,ivarbl)
         syksx(ip)=val
         esyksx(ip)=error
         call MNPOUT(2,chnam,val,error,bnd1,bnd2,ivarbl)
         sycsx(ip)=val
         esycsx(ip)=error
         close(40)
      enddo      
      
      return
      end

**************************************************
**************************************************
*     ETA CALIBRATION
**************************************************
**************************************************
      subroutine eta_calibration(iladder,ia)
      include '../common/commontracker.f'
      include '../common/calib.f'  
      include '../common/common_calibration.f'

*     --------------------------------------------
*     some local variables for eta correction
      real histo2(nhbin2)
c      real eta2corr(nhbin2+1,nviews) !eta2 correction
c      real eta2abs(nhbin2+1)         !eta2 abscissa
*     --------------------------------------------
      character*64 ctemp
      character*10 rzdir

 1010 format('LADDER',i1)
      write(rzdir,1010)iladder
      call HCDIR('//'//rzdir,'')
      call HCDIR('//PAWC/'//rzdir,'')

      do i=1,nviews
         id=id_hi_eta2+i*100+ia
*     eta2 histogram
         call HUNPAK(id,histo2,' ',0) !eta2 histogram
         call HGIVE(id,ctemp,nx,xmi,xma,ny,ymi,yma,nwt,loc)
         sum=HSUM(id)
c            print*,ctemp,nx,xmi,xma,ny,ymi,yma,nwt,loc
         if(i.eq.1)then
            do ie=0,nx
               eta2abs(ie+1)=xmi+ie*(xma-xmi)/nx
            enddo
         endif
         if(sum.eq.0)then
            do ie=0,nx
               eta2corr(ie+1,i)=eta2abs(ie+1)
            enddo            
            goto 3
         endif
*     ------------------------------------------------------
*     if the histogram contains less than NMIN4ETA2 entries,
*     the correction factor is not calculated
*     CHANGED!!!
*     it is calculated anyway....
*     ------------------------------------------------------
         nmin4eta2 = 200
         if(sum.lt.nmin4eta2)then
c$$$            do ie=0,nx
c$$$               eta2corr(ie+1,i)=eta2abs(ie+1)
c$$$            enddo            
            print*,'*** WARNING ***'
            print*,'ETA2 parameters: less than ',nmin4eta2
     $           ,' entries ----> ',rzdir,' view ',i,' ang. bin ',ia
c$$$            goto 3
         endif
*     integral of the eta2 histogram
         eta2corr(1,i)=0.
         do ie=1,nx
            eta2corr(ie+1,i)=eta2corr(ie,i)+histo2(ie)
         enddo
*     eta2 normalization
         angle=avangbin(ia,i)
         x0=eta2fix(i,angle)    !point of reference for correction function
         do ifix=1,nx+1
            if(eta2abs(ifix).ge.x0)goto 2
         enddo
 2       y0=eta2corr(ifix,i)/sum
c     shift=eta2shift(i,angle)
         shift = 0.
         do ie=0,nx
            eta2corr(ie+1,i)=
     $           x0-y0+eta2corr(ie+1,i)/sum-shift
         enddo     
 3       continue
      enddo                     !end loop on views       
      
      
      
      return
      end
*---***---***---***---***---***---***---***---***
*
*     functions to define the point where to "fix"
*     the eta correction
*     and the shift
*
*---***---***---***---***---***---***---***---***
      function eta2fix(iview,angle)

      if(mod(iview,2).eq.0)then !>>>>>>>> X

         eta2fix = -0.5

      else                      !>>>>>>>> Y

         eta2fix = -0.5

      endif


      return
      end


c$$$      function eta2shift(iview,angle)
c$$$
c$$$      parameter (hallm_e=1650.) !cm**2 / Vs
c$$$      parameter (hallm_h=310.)
c$$$
c$$$      include '../common/commontracker.f'
c$$$      include '../common/common_preanalysis.f'
c$$$      character*32 ctemp
c$$$      real hallm
c$$$
c$$$      pi=acos(-1.)
c$$$
c$$$      np=npl(iview)
c$$$      if(mod(iview,2).eq.0)then !>>>>>>>> X
c$$$         id=id_hi_by+np
c$$$         hallm=hallm_h*1e-4
c$$$      else 
c$$$         id=id_hi_bx+np
c$$$         hallm=hallm_e*1e-4
c$$$      endif
c$$$      bfield=HSTATI(id,1,' ',0)
c$$$      bfield=bfield*1e-1
c$$$      tgth_h=-hallm*bfield
c$$$
c$$$      tgth_eff=tan(pi*angle/180.)-tgth_h
c$$$      theta_eff=atan(tgth_eff)*180./pi
c$$$      if(iview.ne.12)theta_eff=-theta_eff     
c$$$      
c$$$      if(mod(iview,2).eq.0)then !>>>>>>>> X
c$$$         
c$$$         eta2shift = -1.549*sin(2*pi*theta_eff/13.59)         
c$$$         eta2shift = eta2shift/51.
c$$$
c$$$      else                      !>>>>>>>> Y
c$$$         
c$$$         eta2shift = -1.348*sin(2*pi*theta_eff/27.73)         
c$$$         eta2shift = eta2shift/67.
c$$$*     NBNBNBNBNB ... da risolvere il problema della simulazione!!!!
c$$$
c$$$      endif
c$$$
c$$$c      print*,iview,np,bfield,angle,theta_eff,'>>>>>>> ',eta2shift
c$$$
c$$$      return
c$$$      end
