/[PAMELA software]/DarthVader/TrackerLevel2/src/F77/analysissubroutines.f
ViewVC logotype

Diff of /DarthVader/TrackerLevel2/src/F77/analysissubroutines.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.9 by pam-fi, Fri Oct 27 16:08:19 2006 UTC revision 1.17 by pam-fi, Thu Jan 11 10:20:58 2007 UTC
# Line 18  Line 18 
18        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
19        include 'common_mini_2.f'        include 'common_mini_2.f'
20        include 'calib.f'        include 'calib.f'
 c      include 'level1.f'  
21        include 'level2.f'        include 'level2.f'
22    
23  c      include 'momanhough_init.f'  c      include 'momanhough_init.f'
# Line 183  c      iflag=0 Line 182  c      iflag=0
182  c      include 'momanhough_init.f'  c      include 'momanhough_init.f'
183                
184        logical FIMAGE            !        logical FIMAGE            !
185          real*8 AL_GUESS(5)
186    
187  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
188  *     STEP 4   (ITERATED until any other physical track isn't found)  *     STEP 4   (ITERATED until any other physical track isn't found)
# Line 233  c$$$            endif Line 233  c$$$            endif
233  c$$$         enddo  c$$$         enddo
234  c$$$         if(ibest.eq.0)goto 880 !>> no good candidates  c$$$         if(ibest.eq.0)goto 880 !>> no good candidates
235    
236           rchi2best=1000000000.  *     -------------------------------------------------------
237    *     order track-candidates according to:
238    *     1st) decreasing n.points
239    *     2nd) increasing chi**2
240    *     -------------------------------------------------------
241             rchi2best=1000000000.
242           ndofbest=0             !(1)           ndofbest=0             !(1)
243           do i=1,ntracks           do i=1,ntracks
244             if(RCHI2_STORE(i).lt.rchi2best.and.             ndof=0               !(1)
245       $          RCHI2_STORE(i).gt.0)then             do ii=1,nplanes      !(1)
246               ndof=0             !(1)               ndof=ndof          !(1)
247               do ii=1,nplanes    !(1)       $            +int(xgood_store(ii,i)) !(1)
248                 ndof=ndof        !(1)       $            +int(ygood_store(ii,i)) !(1)
249       $              +int(xgood_store(ii,i)) !(1)             enddo                !(1)
250       $              +int(ygood_store(ii,i)) !(1)             if(ndof.gt.ndofbest)then !(1)
251               enddo              !(1)               ibest=i
252               if(ndof.ge.ndofbest)then !(1)               rchi2best=RCHI2_STORE(i)
253                 ndofbest=ndof      !(1)
254               elseif(ndof.eq.ndofbest)then !(1)
255                 if(RCHI2_STORE(i).lt.rchi2best.and.
256         $            RCHI2_STORE(i).gt.0)then
257                 ibest=i                 ibest=i
258                 rchi2best=RCHI2_STORE(i)                 rchi2best=RCHI2_STORE(i)
259                 ndofbest=ndof    !(1)                 ndofbest=ndof    !(1)
260               endif              !(1)               endif              !(1)
261             endif             endif
262           enddo           enddo
263    
264    c$$$         rchi2best=1000000000.
265    c$$$         ndofbest=0             !(1)
266    c$$$         do i=1,ntracks
267    c$$$           if(RCHI2_STORE(i).lt.rchi2best.and.
268    c$$$     $          RCHI2_STORE(i).gt.0)then
269    c$$$             ndof=0             !(1)
270    c$$$             do ii=1,nplanes    !(1)
271    c$$$               ndof=ndof        !(1)
272    c$$$     $              +int(xgood_store(ii,i)) !(1)
273    c$$$     $              +int(ygood_store(ii,i)) !(1)
274    c$$$             enddo              !(1)
275    c$$$             if(ndof.ge.ndofbest)then !(1)
276    c$$$               ibest=i
277    c$$$               rchi2best=RCHI2_STORE(i)
278    c$$$               ndofbest=ndof    !(1)
279    c$$$             endif              !(1)
280    c$$$           endif
281    c$$$         enddo
282    
283           if(ibest.eq.0)goto 880 !>> no good candidates           if(ibest.eq.0)goto 880 !>> no good candidates
284  *-------------------------------------------------------------------------------      *-------------------------------------------------------------------------------    
285  *     The best track candidate (ibest) is selected and a new fitting is performed.  *     The best track candidate (ibest) is selected and a new fitting is performed.
# Line 281  c$$$         if(ibest.eq.0)goto 880 !>> Line 310  c$$$         if(ibest.eq.0)goto 880 !>>
310  *     **********************************************************  *     **********************************************************
311  *     ************************** FIT *** FIT *** FIT *** FIT ***  *     ************************** FIT *** FIT *** FIT *** FIT ***
312  *     **********************************************************  *     **********************************************************
313             call guess()
314             do i=1,5
315                AL_GUESS(i)=AL(i)
316             enddo
317    c         print*,'## guess: ',al
318    
319           do i=1,5           do i=1,5
320              AL(i)=dble(AL_STORE(i,icand))                          AL(i)=dble(AL_STORE(i,icand))            
321           enddo           enddo
322            
323           IDCAND = icand         !fitted track-candidate           IDCAND = icand         !fitted track-candidate
324           ifail=0                !error flag in chi2 computation           ifail=0                !error flag in chi2 computation
325           jstep=0                !# minimization steps           jstep=0                !# minimization steps
326    
327           iprint=0           iprint=0
328           if(DEBUG)iprint=1  c         if(DEBUG)iprint=1
329             if(VERBOSE)iprint=1
330             if(DEBUG)iprint=2
331           call mini2(jstep,ifail,iprint)           call mini2(jstep,ifail,iprint)
332           if(ifail.ne.0) then           if(ifail.ne.0) then
333              if(DEBUG)then              if(VERBOSE)then
334                 print *,                 print *,
335       $              '*** MINIMIZATION FAILURE *** (mini2) '       $              '*** MINIMIZATION FAILURE *** (after refinement) '
336       $              ,iev       $              ,iev
337    
338    c$$$               print*,'guess:   ',(al_guess(i),i=1,5)
339    c$$$               print*,'previous: ',(al_store(i,icand),i=1,5)
340    c$$$               print*,'result:   ',(al(i),i=1,5)
341    c$$$               print*,'xgood ',xgood
342    c$$$               print*,'ygood ',ygood
343    c$$$               print*,'----------------------------------------------'
344              endif              endif
345              chi2=-chi2  c            chi2=-chi2
346           endif           endif
347                    
348           if(DEBUG)then           if(DEBUG)then
# Line 408  c     $        rchi2best.lt.15..and. Line 453  c     $        rchi2best.lt.15..and.
453        end        end
454    
455    
   
   
 c$$$************************************************************  
 c$$$  
 c$$$      subroutine readmipparam  
 c$$$              
 c$$$      include 'commontracker.f'  
 c$$$      include 'calib.f'  
 c$$$  
 c$$$      character*60 fname_param  
 c$$$ 201  format('trk-LADDER',i1,'-mip.dat')  
 c$$$      do ilad=1,nladders_view          
 c$$$         write(fname_param,201)ilad  
 c$$$         print *,'Opening file: ',fname_param  
 c$$$         open(10,  
 c$$$     $        FILE='./bin-aux/'//fname_param(1:LNBLNK(fname_param))  
 c$$$     $        ,STATUS='UNKNOWN'  
 c$$$     $        ,IOSTAT=iostat  
 c$$$     $        )  
 c$$$         if(iostat.ne.0)then  
 c$$$            print*,'READMIPPARAM: *** Error in opening file ***'  
 c$$$            return  
 c$$$         endif  
 c$$$         do iv=1,nviews  
 c$$$            read(10,*  
 c$$$     $           ,IOSTAT=iostat  
 c$$$     $           )pip,  
 c$$$     $            mip(int(pip),ilad)  
 c$$$c            print*,ilad,iv,pip,mip(int(pip),ilad)  
 c$$$         enddo  
 c$$$         close(10)  
 c$$$      enddo  
 c$$$  
 c$$$      return  
 c$$$      end  
 c$$$*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  
 c$$$      subroutine readchargeparam  
 c$$$        
 c$$$        
 c$$$      include 'commontracker.f'  
 c$$$      include 'calib.f'  
 c$$$  
 c$$$      character*60 fname_param  
 c$$$ 201  format('charge-l',i1,'.dat')  
 c$$$      do ilad=1,nladders_view          
 c$$$         write(fname_param,201)ilad  
 c$$$         print *,'Opening file: ',fname_param  
 c$$$         open(10,  
 c$$$     $        FILE='./bin-aux/'//fname_param(1:LNBLNK(fname_param))  
 c$$$     $        ,STATUS='UNKNOWN'  
 c$$$     $        ,IOSTAT=iostat  
 c$$$     $        )  
 c$$$         if(iostat.ne.0)then  
 c$$$            print*,'READCHARGEPARAM: *** Error in opening file ***'  
 c$$$            return  
 c$$$         endif  
 c$$$         do ip=1,nplanes  
 c$$$            read(10,*  
 c$$$     $           ,IOSTAT=iostat  
 c$$$     $           )pip,  
 c$$$     $            kch(ip,ilad),cch(ip,ilad),sch(ip,ilad)          
 c$$$c            print*,ilad,ip,pip,kch(ip,ilad),  
 c$$$c     $           cch(ip,ilad),sch(ip,ilad)  
 c$$$         enddo  
 c$$$         close(10)  
 c$$$      enddo  
 c$$$  
 c$$$      return  
 c$$$      end  
 c$$$*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  
 c$$$      subroutine readetaparam  
 c$$$*     -----------------------------------------  
 c$$$*     read eta2,3,4 calibration parameters  
 c$$$*     and fill variables:  
 c$$$*  
 c$$$*     eta2(netabin,nladders_view,nviews)  
 c$$$*     eta3(2*netabin,nladders_view,nviews)  
 c$$$*     eta4(2*netabin,nladders_view,nviews)  
 c$$$*  
 c$$$      include 'commontracker.f'  
 c$$$      include 'calib.f'  
 c$$$  
 c$$$      character*40 fname_binning  
 c$$$      character*40 fname_param  
 c$$$c      character*120 cmd1  
 c$$$c      character*120 cmd2  
 c$$$  
 c$$$  
 c$$$******retrieve ANGULAR BINNING info  
 c$$$      fname_binning='binning.dat'  
 c$$$      print *,'Opening file: ',fname_binning  
 c$$$      open(10,  
 c$$$     $     FILE='./bin-aux/'//fname_binning(1:LNBLNK(fname_binning))  
 c$$$     $     ,STATUS='UNKNOWN'  
 c$$$     $     ,IOSTAT=iostat  
 c$$$     $     )  
 c$$$      if(iostat.ne.0)then  
 c$$$         print*,'READETAPARAM: *** Error in opening file ***'  
 c$$$         return  
 c$$$      endif  
 c$$$      print*,'---- ANGULAR BINNING ----'  
 c$$$      print*,'Bin   -   angL   -   angR'  
 c$$$ 101  format(i2,'       ',f6.2,'     ',f6.2)  
 c$$$      do ibin=1,nangmax  
 c$$$         read(10,*  
 c$$$     $        ,IOSTAT=iostat  
 c$$$     $        )xnn,angL(ibin),angR(ibin)  
 c$$$         if(iostat.ne.0)goto 1000  
 c$$$         write(*,101)int(xnn),angL(ibin),angR(ibin)  
 c$$$      enddo          
 c$$$ 1000 nangbin=int(xnn)  
 c$$$      close(10)  
 c$$$      print*,'-------------------------'  
 c$$$        
 c$$$  
 c$$$  
 c$$$      do ieta=2,4               !loop on eta 2,3,4          
 c$$$******retrieve correction parameters  
 c$$$ 200     format(' Opening eta',i1,' files...')  
 c$$$         write(*,200)ieta  
 c$$$  
 c$$$ 201     format('eta',i1,'-bin',i1,'-l',i1,'.dat')  
 c$$$ 202     format('eta',i1,'-bin',i2,'-l',i1,'.dat')  
 c$$$         do iang=1,nangbin  
 c$$$            do ilad=1,nladders_view  
 c$$$               if(iang.lt.10)write(fname_param,201)ieta,iang,ilad  
 c$$$               if(iang.ge.10)write(fname_param,202)ieta,iang,ilad  
 c$$$c               print *,'Opening file: ',fname_param  
 c$$$               open(10,  
 c$$$     $             FILE='./bin-aux/'//fname_param(1:LNBLNK(fname_param))  
 c$$$     $              ,STATUS='UNKNOWN'  
 c$$$     $              ,IOSTAT=iostat  
 c$$$     $              )  
 c$$$               if(iostat.ne.0)then  
 c$$$                  print*,'READETAPARAM: *** Error in opening file ***'  
 c$$$                  return  
 c$$$               endif  
 c$$$               do ival=1,netavalmax  
 c$$$                  if(ieta.eq.2)read(10,*  
 c$$$     $                 ,IOSTAT=iostat  
 c$$$     $                 )  
 c$$$     $                 eta2(ival,iang),  
 c$$$     $                 (feta2(ival,iv,ilad,iang),iv=1,nviews)  
 c$$$                  if(ieta.eq.3)read(10,*  
 c$$$     $                 ,IOSTAT=iostat  
 c$$$     $                 )  
 c$$$     $                 eta3(ival,iang),  
 c$$$     $                 (feta3(ival,iv,ilad,iang),iv=1,nviews)  
 c$$$                  if(ieta.eq.4)read(10,*  
 c$$$     $                 ,IOSTAT=iostat  
 c$$$     $                 )  
 c$$$     $                 eta4(ival,iang),  
 c$$$     $                 (feta4(ival,iv,ilad,iang),iv=1,nviews)  
 c$$$                  if(iostat.ne.0)then  
 c$$$                     netaval=ival-1  
 c$$$c$$$                     if(eta2(1,iang).ne.-eta2(netaval,iang))  
 c$$$c$$$     $                    print*,'**** ERROR on parameters !!! ****'  
 c$$$                     goto 2000  
 c$$$                  endif  
 c$$$               enddo  
 c$$$ 2000          close(10)  
 c$$$*               print*,'... done'  
 c$$$            enddo  
 c$$$         enddo  
 c$$$  
 c$$$      enddo                     !end loop on eta 2,3,4  
 c$$$  
 c$$$  
 c$$$      return  
 c$$$      end  
 c$$$  
   
456                
457  ************************************************************  ************************************************************
458  ************************************************************  ************************************************************
# Line 660  c      include 'level1.f' Line 533  c      include 'level1.f'
533        include 'common_align.f'        include 'common_align.f'
534        include 'common_mech.f'        include 'common_mech.f'
535        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
536        include 'common_resxy.f'  c      include 'common_resxy.f'
537    
538  c      logical DEBUG  c      logical DEBUG
539  c      common/dbg/DEBUG  c      common/dbg/DEBUG
# Line 695  c      double precision xi_B,yi_B,zi_B Line 568  c      double precision xi_B,yi_B,zi_B
568        xPAM_B = 0.        xPAM_B = 0.
569        yPAM_B = 0.        yPAM_B = 0.
570        zPAM_B = 0.        zPAM_B = 0.
571    c      print*,'## xyz_PAM: ',icx,icy,sensor,PFAx,PFAy,angx,angy
572  *     -----------------  *     -----------------
573  *     CLUSTER X  *     CLUSTER X
574  *     -----------------  *     -----------------
# Line 814  c            resyPAM = ris_eta(icy,angy) Line 687  c            resyPAM = ris_eta(icy,angy)
687    
688        endif        endif
689    
690          c      print*,'## stripx,stripy ',stripx,stripy
691    
692  c===========================================================  c===========================================================
693  C     COUPLE  C     COUPLE
694  C===========================================================  C===========================================================
# Line 1016  c         print*,'A-(',xPAM_A,yPAM_A,') Line 890  c         print*,'A-(',xPAM_A,yPAM_A,')
890                            
891        endif        endif
892                    
893    
894    c      print*,'## xPAM,yPAM,zPAM       ',xPAM,yPAM,zPAM
895    c      print*,'## xPAM_A,yPAM_A,zPAM_A ',xPAM_A,yPAM_A,zPAM_A
896    c      print*,'## xPAM_B,yPAM_B,zPAM_B ',xPAM_B,yPAM_B,zPAM_B
897    
898   100  continue   100  continue
899        end        end
900    
# Line 1429  c      include 'common_analysis.f' Line 1308  c      include 'common_analysis.f'
1308  *************************************************************************  *************************************************************************
1309  *************************************************************************  *************************************************************************
1310  *************************************************************************  *************************************************************************
 c$$$      subroutine book_debug  
 c$$$  
 c$$$      include 'commontracker.f'  
 c$$$      include 'common_momanhough.f'  
 c$$$      include 'common_level2debug.f'  
 c$$$  
 c$$$      character*35 block1,block2,block3!,block4  
 c$$$     $     ,block5!,block6        
 c$$$        
 c$$$* * * * * * * * * * * * * * * * * * * * * * * *  
 c$$$*     HOUGH TRANSFORM PARAMETERS  
 c$$$        
 c$$$      call HBOOK2(1003  
 c$$$     $     ,'y vs tg thyz'  
 c$$$     $     ,300,-1.,1.         !x          
 c$$$     $     ,3000,-70.,70.,0.)   !y  
 c$$$  
 c$$$      call HBOOK1(1004  
 c$$$     $     ,'Dy'  
 c$$$     $     ,100,0.,2.,0.)    
 c$$$  
 c$$$      call HBOOK1(1005  
 c$$$     $     ,'D thyz'  
 c$$$     $     ,100,0.,.05,0.)    
 c$$$  
 c$$$  
 c$$$  
 c$$$*     DEBUG ntuple:  
 c$$$      call HBNT(ntp_level2+1,'LEVEL2',' ')  
 c$$$        
 c$$$      call HBNAME(ntp_level2+1,'EVENT',good2_nt,  
 c$$$     $     'GOOD2:L,NEV2:I')  
 c$$$  
 c$$$ 411  format('NDBLT:I::[0,',I5,']')  
 c$$$      write(block1,411) ndblt_max_nt  
 c$$$      call HBNAME(ntp_level2+1,'HOUGH YZ',ndblt_nt,  
 c$$$     $     block1//'  
 c$$$     $     ,ALFAYZ1(NDBLT):R  
 c$$$     $     ,ALFAYZ2(NDBLT):R  
 c$$$     $     ,DB_CLOUD(NDBLT):I  
 c$$$     $     ')    
 c$$$  
 c$$$ 412  format('NTRPT:I::[0,',I5,']')  
 c$$$      write(block2,412) ntrpt_max_nt  
 c$$$      call HBNAME(ntp_level2+1,'HOUGH XZ',NTRPT_nt,  
 c$$$     $     block2//'  
 c$$$     $     ,ALFAXZ1(NTRPT):R  
 c$$$     $     ,ALFAXZ2(NTRPT):R    
 c$$$     $     ,ALFAXZ3(NTRPT):R  
 c$$$     $     ,TR_CLOUD(NTRPT):I  
 c$$$     $     ')  
 c$$$        
 c$$$    
 c$$$ 413  format('NCLOUDS_YZ:I::[0,',I4,']')  
 c$$$c$$$ 414  format('DB_CLOUD(',I4,'):I')  
 c$$$      write(block3,413) ncloyz_max  
 c$$$c$$$      write(block4,414) ndblt_max_nt  
 c$$$      call HBNAME(ntp_level2+1,'CLOUD YZ',NCLOUDS_YZ,  
 c$$$     $     block3//'  
 c$$$     $     ,ALFAYZ1_AV(NCLOUDS_YZ):R  
 c$$$     $     ,ALFAYZ2_AV(NCLOUDS_YZ):R  
 c$$$     $     ,PTCLOUD_YZ(NCLOUDS_YZ):I'  
 c$$$c$$$     $     ,'//block4  
 c$$$     $     )  
 c$$$  
 c$$$ 415  format('NCLOUDS_XZ:I::[0,',I4,']')  
 c$$$c$$$ 416  format('TR_CLOUD(',I5,'):I')  
 c$$$      write(block5,415) ncloxz_max  
 c$$$c$$$      write(block6,416) ntrpt_max_nt  
 c$$$      call HBNAME(ntp_level2+1,'CLOUD XZ',NCLOUDS_XZ,  
 c$$$     $     block5//'  
 c$$$     $     ,ALFAXZ1_AV(NCLOUDS_XZ):R  
 c$$$     $     ,ALFAXZ2_AV(NCLOUDS_XZ):R  
 c$$$     $     ,ALFAXZ3_AV(NCLOUDS_XZ):R  
 c$$$     $     ,PTCLOUD_XZ(NCLOUDS_XZ):I'  
 c$$$c$$$     $     ,'//block6  
 c$$$     $     )  
 c$$$  
 c$$$        
 c$$$      return  
 c$$$      end  
 ***...***...***...***...***...***...***...***...***...***...***...***...***...***...***...***  
 *  
 *  
 *  
 *  
 *  
 *  
 *  
 *  
 *  
 ***...***...***...***...***...***...***...***...***...***...***...***...***...***...***...***  
 c$$$      subroutine book_level2  
 c$$$c*****************************************************  
 c$$$cccccc 11/9/2005 modified by david fedele  
 c$$$cccccc 07/10/2005 modified by elena vannuccini --> (2)  
 c$$$c*****************************************************  
 c$$$  
 c$$$      include 'commontracker.f'  
 c$$$      include 'common_momanhough.f'  
 c$$$      include 'level2.f'  
 c$$$  
 c$$$      character*35 block1,block2  
 c$$$  
 c$$$c      print*,'__________ booking LEVEL2 n-tuple __________'  
 c$$$  
 c$$$*     LEVEL1 ntuple:  
 c$$$      call HBNT(ntp_level2,'LEVEL2',' ')  
 c$$$        
 c$$$c*****************************************************  
 c$$$cccccc 11/9/2005 modified by david fedele  
 c$$$c      call HBNAME(ntp_level2,'EVENT',good2,'GOOD2:L,NEV2:I')  
 c$$$cccccc 06/10/2005 modified by elena vannuccini  
 c$$$c      call HBNAME(ntp_level2,'GENERAL',good2,'GOOD2:L,NEV2:I  
 c$$$c     $     ,WHIC_CALIB:I,SWCODE:I')  
 c$$$      call HBNAME(ntp_level2,'GENERAL',good2,'GOOD2:L,NEV2:I  
 c$$$     $     ,WHICH_CALIB:I,SWCODE:I,CRC(12):L')  
 c$$$c*********************************************************  
 c$$$    
 c$$$  
 c$$$c$$$# ifndef TEST2003  
 c$$$c$$$  
 c$$$c$$$      call HBNAME(ntp_level2,'CPU',pkt_type  
 c$$$c$$$     $     ,'PKT_TYPE:I::[0,50]  
 c$$$c$$$     $     ,PKT_NUM:I  
 c$$$c$$$     $     ,OBT:I'//  
 c$$$c$$$c********************************************************  
 c$$$c$$$cccccc 11/9/2005 modified by david fedele  
 c$$$c$$$c     $     ,WHICH_CALIB:I::[0,50]')  
 c$$$c$$$     $     ',CPU_CRC:L')  
 c$$$c$$$c********************************************************  
 c$$$c$$$  
 c$$$c$$$# endif  
 c$$$  
 c$$$ 417  format('NTRK:I::[0,',I4,']')  
 c$$$ 418  format(',IMAGE(NTRK):I::[0,',I4,']')  
 c$$$      write(block1,417)NTRKMAX  
 c$$$      write(block2,418)NTRKMAX  
 c$$$      call HBNAME(ntp_level2,'TRACKS',NTRK,  
 c$$$     $     block1//  
 c$$$     $     block2//'  
 c$$$     $     ,XM(6,NTRK):R  
 c$$$     $     ,YM(6,NTRK):R  
 c$$$     $     ,ZM(6,NTRK):R  
 c$$$     $     ,RESX(6,NTRK):R  
 c$$$     $     ,RESY(6,NTRK):R  
 c$$$     $     ,AL(5,NTRK):R  
 c$$$     $     ,COVAL(5,5,NTRK):R  
 c$$$     $     ,CHI2(NTRK):R  
 c$$$     $     ,XGOOD(6,NTRK):I::[0,1]  
 c$$$     $     ,YGOOD(6,NTRK):I::[0,1]  
 c$$$     $     ,XV(6,NTRK):R  
 c$$$     $     ,YV(6,NTRK):R  
 c$$$     $     ,ZV(6,NTRK):R  
 c$$$     $     ,AXV(6,NTRK):R  
 c$$$     $     ,AYV(6,NTRK):R'//  
 c$$$c*****************************************************  
 c$$$cccccc 11/9/2005 modified by david fedele  
 c$$$c     $     ,DEDXP(6,NTRK):R'//  
 c$$$c     $     ')  
 c$$$     $     ',DEDX_X(6,NTRK):R  
 c$$$     $     ,DEDX_Y(6,NTRK):R'//  
 c$$$c****************************************************  
 c$$$cccccc 06/10/2005 modified by elena vannuccini  
 c$$$c     $     ,CRC(12):L  
 c$$$     $     ',BdL(NTRK):R'  
 c$$$     $     )  
 c$$$c****************************************************  
 c$$$  
 c$$$    
 c$$$      call HBNAME(ntp_level2,'SINGLETX',nclsx,  
 c$$$c*****************************************************  
 c$$$cccccc 11/9/2005 modified by david fedele  
 c$$$c     $     'NCLSX(6):I,NCLSY(6):I')  
 c$$$     $     'NCLSX:I::[0,500],PLANEX(NCLSX):I  
 c$$$     $     ,XS(2,NCLSX):R,SGNLXS(NCLSX):R') !(2)  
 c$$$c    $     ,XS(NCLSX):R,SGNLXS(NCLSX):R')   !(2)  
 c$$$      call HBNAME(ntp_level2,'SINGLETY',nclsy,  
 c$$$     $     'NCLSY:I::[0,500],PLANEY(NCLSY):I  
 c$$$     $     ,YS(2,NCLSY):R,SGNLYS(NCLSY):R') !(2)  
 c$$$c    $     ,YS(NCLSY):R,SGNLYS(NCLSY):R')   !(2)  
 c$$$      return  
 c$$$      end  
   
 c$$$      subroutine fill_level2_clouds  
 c$$$c*****************************************************  
 c$$$c     29/11/2005 created by elena vannuccini  
 c$$$c*****************************************************  
 c$$$  
 c$$$*     -------------------------------------------------------  
 c$$$*     This routine fills the  variables related to the hough  
 c$$$*     transform, for the debig n-tuple  
 c$$$*     -------------------------------------------------------  
 c$$$  
 c$$$      include 'commontracker.f'  
 c$$$      include 'common_momanhough.f'  
 c$$$      include 'common_level2debug.f'  
 c$$$      include 'level2.f'  
 c$$$  
 c$$$      good2_nt=.true.!good2  
 c$$$c      nev2_nt=nev2  
 c$$$        
 c$$$      if(.false.  
 c$$$     $     .or.ntrpt.gt.ntrpt_max_nt  
 c$$$     $     .or.ndblt.gt.ndblt_max_nt  
 c$$$     $     .or.NCLOUDS_XZ.gt.ncloxz_max  
 c$$$     $     .or.NCLOUDS_yZ.gt.ncloyz_max  
 c$$$     $     )then  
 c$$$         good2_nt=.false.  
 c$$$         ntrpt_nt=0  
 c$$$         ndblt_nt=0  
 c$$$         NCLOUDS_XZ_nt=0  
 c$$$         NCLOUDS_YZ_nt=0  
 c$$$      else  
 c$$$         ndblt_nt=ndblt  
 c$$$         ntrpt_nt=ntrpt  
 c$$$         if(ndblt.ne.0)then  
 c$$$            do id=1,ndblt  
 c$$$               alfayz1_nt(id)=alfayz1(id) !Y0  
 c$$$               alfayz2_nt(id)=alfayz2(id) !tg theta-yz  
 c$$$c               db_cloud_nt(id)=db_cloud(id)  
 c$$$            enddo  
 c$$$         endif  
 c$$$         if(ndblt.ne.0)then  
 c$$$            do it=1,ntrpt  
 c$$$               alfaxz1_nt(it)=alfaxz1(it) !X0  
 c$$$               alfaxz2_nt(it)=alfaxz2(it) !tg theta-xz  
 c$$$               alfaxz3_nt(it)=alfaxz3(it) !1/r  
 c$$$c               tr_cloud_nt(it)=tr_cloud(it)  
 c$$$            enddo  
 c$$$         endif  
 c$$$         nclouds_yz_nt=nclouds_yz  
 c$$$         nclouds_xz_nt=nclouds_xz  
 c$$$         if(nclouds_yz.ne.0)then  
 c$$$            nnn=0  
 c$$$            do iyz=1,nclouds_yz  
 c$$$               ptcloud_yz_nt(iyz)=ptcloud_yz(iyz)  
 c$$$               alfayz1_av_nt(iyz)=alfayz1_av(iyz)  
 c$$$               alfayz2_av_nt(iyz)=alfayz2_av(iyz)  
 c$$$               nnn=nnn+ptcloud_yz(iyz)  
 c$$$            enddo  
 c$$$            do ipt=1,nnn  
 c$$$               db_cloud_nt(ipt)=db_cloud(ipt)  
 c$$$            enddo  
 c$$$c            print*,'#### ntupla #### ptcloud_yz '  
 c$$$c     $           ,(ptcloud_yz(i),i=1,nclouds_yz)  
 c$$$c            print*,'#### ntupla #### db_cloud '  
 c$$$c     $           ,(db_cloud(i),i=1,nnn)  
 c$$$         endif  
 c$$$         if(nclouds_xz.ne.0)then  
 c$$$            nnn=0  
 c$$$            do ixz=1,nclouds_xz  
 c$$$               ptcloud_xz_nt(ixz)=ptcloud_xz(ixz)  
 c$$$               alfaxz1_av_nt(ixz)=alfaxz1_av(ixz)  
 c$$$               alfaxz2_av_nt(ixz)=alfaxz2_av(ixz)  
 c$$$               alfaxz3_av_nt(ixz)=alfaxz3_av(ixz)  
 c$$$               nnn=nnn+ptcloud_xz(ixz)                
 c$$$            enddo  
 c$$$            do ipt=1,nnn  
 c$$$               tr_cloud_nt(ipt)=tr_cloud(ipt)  
 c$$$            enddo  
 c$$$c            print*,'#### ntupla #### ptcloud_xz '  
 c$$$c     $           ,(ptcloud_xz(i),i=1,nclouds_xz)  
 c$$$c            print*,'#### ntupla #### tr_cloud '  
 c$$$c     $           ,(tr_cloud(i),i=1,nnn)  
 c$$$         endif  
 c$$$      endif  
 c$$$      end  
1311                
1312    
1313  ***************************************************  ***************************************************
# Line 1754  c      include 'level1.f' Line 1365  c      include 'level1.f'
1365        enddo        enddo
1366  *     mask views with too many clusters  *     mask views with too many clusters
1367        do iv=1,nviews        do iv=1,nviews
1368           if( ncl_view(iv).gt. nclustermax)then           if( ncl_view(iv).gt. nclusterlimit)then
1369              mask_view(iv) = 1              mask_view(iv) = 1
1370              print*,' * WARNING * cl_to_couple: n.clusters > '              if(DEBUG)print*,' * WARNING * cl_to_couple: n.clusters > '
1371       $           ,nclustermax,' on view ', iv,' --> masked!'       $           ,nclusterlimit,' on view ', iv,' --> masked!'
1372           endif           endif
1373        enddo        enddo
1374    
# Line 1893  c                  cut = chcut * sch(npl Line 1504  c                  cut = chcut * sch(npl
1504                       goto 20    !charge not consistent                       goto 20    !charge not consistent
1505                    endif                    endif
1506                 endif                 endif
1507                  
1508  *     ------------------> COUPLE <------------------                 if(ncp_plane(nplx).gt.ncouplemax)then
 *     check to do not overflow vector dimentions  
 c$$$               if(ncp_plane(nplx).gt.ncouplemax)then  
 c$$$                  if(DEBUG)print*,  
 c$$$     $                    ' ** warning ** number of identified'//  
 c$$$     $                    ' couples on plane ',nplx,  
 c$$$     $                    ' exceeds vector dimention'//  
 c$$$     $                    ' ( ',ncouplemax,' )'  
 c$$$c     good2=.false.  
 c$$$c     goto 880   !fill ntp and go to next event  
 c$$$                  iflag=1  
 c$$$                  return  
 c$$$               endif  
                 
                if(ncp_plane(nplx).eq.ncouplemax)then  
1509                    if(verbose)print*,                    if(verbose)print*,
1510       $                 '** warning ** number of identified '//       $                 '** warning ** number of identified '//
1511       $                 'couples on plane ',nplx,       $                 'couples on plane ',nplx,
1512       $                 'exceeds vector dimention '       $                 'exceeds vector dimention '
1513       $                 ,'( ',ncouplemax,' ) NB - THIS SHOULD NOT HAPPEN'       $                 ,'( ',ncouplemax,' ) --> masked!'
1514  c     good2=.false.                    mask_view(nviewx(nplx)) = 2
1515  c     goto 880   !fill ntp and go to next event                                        mask_view(nviewy(nply)) = 2
1516                    iflag=1                    goto 10
                   return  
1517                 endif                 endif
1518                                
1519    *     ------------------> COUPLE <------------------
1520                 ncp_plane(nplx) = ncp_plane(nplx) + 1                 ncp_plane(nplx) = ncp_plane(nplx) + 1
1521                 clx(nplx,ncp_plane(nplx))=icx                 clx(nplx,ncp_plane(nplx))=icx
1522                 cly(nply,ncp_plane(nplx))=icy                 cly(nply,ncp_plane(nplx))=icy
1523                 cl_single(icx)=0                 cl_single(icx)=0
1524                 cl_single(icy)=0                 cl_single(icy)=0
             endif                                
1525  *     ----------------------------------------------  *     ----------------------------------------------
1526    
  20         continue  
          enddo                  !end loop on clusters(Y)  
           
  10      continue  
       enddo                     !end loop on clusters(X)  
         
         
       do icl=1,nclstr1  
          if(cl_single(icl).eq.1)then  
             ip=npl(VIEW(icl))  
             ncls(ip)=ncls(ip)+1  
             cls(ip,ncls(ip))=icl  
          endif  
       enddo  
         
         
       if(DEBUG)then  
          print*,'clusters  ',nclstr1  
          print*,'good    ',(cl_good(i),i=1,nclstr1)  
          print*,'singles ',(cl_single(i),i=1,nclstr1)  
          print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)  
       endif  
         
       do ip=1,6  
          ncp_tot=ncp_tot+ncp_plane(ip)  
       enddo  
 c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)  
         
       if(ncp_tot.gt.ncp_max)then  
          if(verbose)print*,  
      $           '** warning ** number of identified '//  
      $           'couples exceeds upper limit for Hough tr. '  
      $           ,'( ',ncp_max,' )'              
 c            good2=.false.  
 c     goto 880       !fill ntp and go to next event  
          iflag=1  
          return  
       endif  
         
       return  
       end  
         
 ***************************************************  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 **************************************************  
       subroutine cl_to_couples_nocharge(iflag)  
   
       include 'commontracker.f'  
       include 'level1.f'  
       include 'common_momanhough.f'  
 c      include 'momanhough_init.f'  
       include 'calib.f'  
 c      include 'level1.f'  
   
   
 *     output flag  
 *     --------------  
 *     0 = good event  
 *     1 = bad event  
 *     --------------  
       integer iflag  
   
       integer badseed,badcl  
   
 *     init variables  
       ncp_tot=0  
       do ip=1,nplanes  
          do ico=1,ncouplemax  
             clx(ip,ico)=0  
             cly(ip,ico)=0  
          enddo  
          ncp_plane(ip)=0  
          do icl=1,nclstrmax_level2  
             cls(ip,icl)=1  
          enddo  
          ncls(ip)=0  
       enddo  
       do icl=1,nclstrmax_level2  
          cl_single(icl)=1  
          cl_good(icl)=0  
       enddo  
         
 *     start association  
       ncouples=0  
       do icx=1,nclstr1          !loop on cluster (X)  
          if(mod(VIEW(icx),2).eq.1)goto 10  
           
 *     ----------------------------------------------------  
 *     cut on charge (X VIEW)  
          if(dedx(icx).lt.dedx_x_min)then  
             cl_single(icx)=0  
             goto 10  
          endif  
 *     cut BAD (X VIEW)              
          badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx)))  
          ifirst=INDSTART(icx)  
          if(icx.ne.nclstr1) then  
             ilast=INDSTART(icx+1)-1  
          else  
             ilast=TOTCLLENGTH  
          endif  
          badcl=badseed  
          do igood=-ngoodstr,ngoodstr  
             ibad=1  
             if((INDMAX(icx)+igood).gt.ifirst.and.  
      $           (INDMAX(icx)+igood).lt.ilast.and.  
      $           .true.)then  
                ibad=BAD(VIEW(icx),  
      $              nvk(MAXS(icx)+igood),  
      $              nst(MAXS(icx)+igood))  
             endif  
             badcl=badcl*ibad  
          enddo  
          if(badcl.eq.0)then     !<<<<<<<<<<<<<< BAD cut  
             cl_single(icx)=0    !<<<<<<<<<<<<<< BAD cut  
             goto 10             !<<<<<<<<<<<<<< BAD cut  
          endif                  !<<<<<<<<<<<<<< BAD cut  
 *     ----------------------------------------------------  
           
          cl_good(icx)=1  
          nplx=npl(VIEW(icx))  
          nldx=nld(MAXS(icx),VIEW(icx))  
           
          do icy=1,nclstr1       !loop on cluster (Y)  
             if(mod(VIEW(icy),2).eq.0)goto 20  
               
 *     ----------------------------------------------------  
 *     cut on charge (Y VIEW)  
             if(dedx(icy).lt.dedx_y_min)then  
                cl_single(icy)=0  
                goto 20  
             endif  
 *     cut BAD (Y VIEW)              
             badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy)))  
             ifirst=INDSTART(icy)  
             if(icy.ne.nclstr1) then  
                ilast=INDSTART(icy+1)-1  
             else  
                ilast=TOTCLLENGTH  
             endif  
             badcl=badseed  
             do igood=-ngoodstr,ngoodstr  
                ibad=1  
                if((INDMAX(icy)+igood).gt.ifirst.and.  
      $              (INDMAX(icy)+igood).lt.ilast.and.  
      $              .true.)  
      $              ibad=BAD(VIEW(icy),  
      $              nvk(MAXS(icy)+igood),  
      $              nst(MAXS(icy)+igood))  
                badcl=badcl*ibad  
             enddo  
             if(badcl.eq.0)then  !<<<<<<<<<<<<<< BAD cut  
                cl_single(icy)=0 !<<<<<<<<<<<<<< BAD cut  
                goto 20          !<<<<<<<<<<<<<< BAD cut  
             endif               !<<<<<<<<<<<<<< BAD cut  
 *     ----------------------------------------------------  
               
               
             cl_good(icy)=1                    
             nply=npl(VIEW(icy))  
             nldy=nld(MAXS(icy),VIEW(icy))  
               
 *     ----------------------------------------------  
 *     CONDITION TO FORM A COUPLE  
 *     ----------------------------------------------  
 *     geometrical consistency (same plane and ladder)  
             if(nply.eq.nplx.and.nldy.eq.nldx)then  
 *     charge correlation  
 *     ===========================================================  
 *     this version of the subroutine is used for the calibration  
 *     thus charge-correlation selection is obviously removed  
 *     ===========================================================  
 c$$$               ddd=(dedx(icy)  
 c$$$     $              -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))  
 c$$$               ddd=ddd/sqrt(kch(nplx,nldx)**2+1)  
 c$$$               cut=chcut*sch(nplx,nldx)  
 c$$$               if(abs(ddd).gt.cut)goto 20 !charge not consistent  
 *     ===========================================================  
                 
                 
 *     ------------------> COUPLE <------------------  
 *     check to do not overflow vector dimentions  
 c$$$               if(ncp_plane(nplx).gt.ncouplemax)then  
 c$$$                  if(DEBUG)print*,  
 c$$$     $                    ' ** warning ** number of identified'//  
 c$$$     $                    ' couples on plane ',nplx,  
 c$$$     $                    ' exceeds vector dimention'//  
 c$$$     $                    ' ( ',ncouplemax,' )'  
 c$$$c     good2=.false.  
 c$$$c     goto 880   !fill ntp and go to next event  
 c$$$                  iflag=1  
 c$$$                  return  
 c$$$               endif  
                 
                if(ncp_plane(nplx).eq.ncouplemax)then  
                   if(verbose)print*,  
      $                 '** warning ** number of identified '//  
      $                 'couples on plane ',nplx,  
      $                 'exceeds vector dimention '  
      $                 ,'( ',ncouplemax,' )'  
 c     good2=.false.  
 c     goto 880   !fill ntp and go to next event                      
                   iflag=1  
                   return  
                endif  
                 
                ncp_plane(nplx) = ncp_plane(nplx) + 1  
                clx(nplx,ncp_plane(nplx))=icx  
                cly(nply,ncp_plane(nplx))=icy  
                cl_single(icx)=0  
                cl_single(icy)=0  
1527              endif                                            endif                              
 *     ----------------------------------------------  
1528    
1529   20         continue   20         continue
1530           enddo                  !end loop on clusters(Y)           enddo                  !end loop on clusters(Y)
# Line 2171  c     goto 880   !fill ntp and go to nex Line 1550  c     goto 880   !fill ntp and go to nex
1550        endif        endif
1551                
1552        do ip=1,6        do ip=1,6
1553           ncp_tot=ncp_tot+ncp_plane(ip)           ncp_tot = ncp_tot + ncp_plane(ip)
1554        enddo        enddo
 c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)  
         
       if(ncp_tot.gt.ncp_max)then  
          if(verbose)print*,  
      $           '** warning ** number of identified '//  
      $           'couples exceeds upper limit for Hough tr. '  
      $           ,'( ',ncp_max,' )'              
 c            good2=.false.  
 c     goto 880       !fill ntp and go to next event  
          iflag=1  
          return  
       endif  
1555                
1556        return        return
1557        end        end
   
1558                
1559  ***************************************************  ***************************************************
1560  *                                                 *  *                                                 *
# Line 2200  c     goto 880       !fill ntp and go to Line 1566  c     goto 880       !fill ntp and go to
1566  **************************************************  **************************************************
1567    
1568        subroutine cp_to_doubtrip(iflag)        subroutine cp_to_doubtrip(iflag)
 c*****************************************************  
 c     02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
1569    
1570        include 'commontracker.f'        include 'commontracker.f'
1571        include 'level1.f'        include 'level1.f'
1572        include 'common_momanhough.f'        include 'common_momanhough.f'
 c      include 'momanhough_init.f'  
1573        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
1574        include 'common_mini_2.f'        include 'common_mini_2.f'
1575        include 'calib.f'        include 'calib.f'
 c      include 'level1.f'  
1576    
1577    
1578  *     output flag  *     output flag
# Line 2245  c      double precision xm3,ym3,zm3 Line 1606  c      double precision xm3,ym3,zm3
1606  *     -----------------------------  *     -----------------------------
1607    
1608    
1609    *     --------------------------------------------
1610    *     put a limit to the maximum number of couples
1611    *     per plane, in order to apply hough transform
1612    *     (couples recovered during track refinement)
1613    *     --------------------------------------------
1614          do ip=1,nplanes
1615             if(ncp_plane(ip).gt.ncouplelimit)then
1616                mask_view(nviewx(ip)) = 8
1617                mask_view(nviewy(ip)) = 8
1618             endif
1619          enddo
1620    
1621    
1622        ndblt=0                   !number of doublets        ndblt=0                   !number of doublets
1623        ntrpt=0                   !number of triplets        ntrpt=0                   !number of triplets
1624                
1625        do ip1=1,(nplanes-1)      !loop on planes  - COPPIA 1        do ip1=1,(nplanes-1)      !loop on planes  - COPPIA 1
1626           do is1=1,2             !loop on sensors - COPPIA 1           if(  mask_view(nviewx(ip1)).ne.0 .or.
1627                     $        mask_view(nviewy(ip1)).ne.0 )goto 10 !skip plane
1628             do is1=1,2             !loop on sensors - COPPIA 1            
1629              do icp1=1,ncp_plane(ip1) !loop on COPPIA 1              do icp1=1,ncp_plane(ip1) !loop on COPPIA 1
1630                 icx1=clx(ip1,icp1)                 icx1=clx(ip1,icp1)
1631                 icy1=cly(ip1,icp1)                 icy1=cly(ip1,icp1)
# Line 2261  c               call xyz_PAM(icx1,icy1,i Line 1635  c               call xyz_PAM(icx1,icy1,i
1635                 ym1=yPAM                 ym1=yPAM
1636                 zm1=zPAM                                   zm1=zPAM                  
1637  c     print*,'***',is1,xm1,ym1,zm1  c     print*,'***',is1,xm1,ym1,zm1
1638    
1639                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2
1640                      if(  mask_view(nviewx(ip2)).ne.0 .or.
1641         $                 mask_view(nviewy(ip2)).ne.0 )goto 20 !skip plane
1642                    do is2=1,2    !loop on sensors -ndblt COPPIA 2                    do is2=1,2    !loop on sensors -ndblt COPPIA 2
1643                                            
1644                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2
# Line 2286  c     $                       (icx2,icy2 Line 1663  c     $                       (icx2,icy2
1663       $                          ,'( ',ndblt_max,' )'       $                          ,'( ',ndblt_max,' )'
1664  c                           good2=.false.  c                           good2=.false.
1665  c                           goto 880 !fill ntp and go to next event  c                           goto 880 !fill ntp and go to next event
1666                               do iv=1,12
1667                                  mask_view(iv) = 3
1668                               enddo
1669                             iflag=1                             iflag=1
1670                             return                             return
1671                          endif                          endif
# Line 2319  c$$$ Line 1699  c$$$
1699  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1700    
1701    
1702                          if(ip2.eq.nplanes)goto 30 !no possible combination with 3 couples                          if(ip2.eq.nplanes)goto 31 !no possible combination with 3 couples
1703    
1704                          do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3                          do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3
1705                               if(  mask_view(nviewx(ip3)).ne.0 .or.
1706         $                          mask_view(nviewy(ip3)).ne.0 )goto 30 !skip plane
1707                             do is3=1,2 !loop on sensors - COPPIA 3                             do is3=1,2 !loop on sensors - COPPIA 3
1708                                                                
1709                                do icp3=1,ncp_plane(ip3) !loop on COPPIA 3                                do icp3=1,ncp_plane(ip3) !loop on COPPIA 3
# Line 2356  c     $                                 Line 1739  c     $                                
1739       $                    ,'( ',ntrpt_max,' )'       $                    ,'( ',ntrpt_max,' )'
1740  c                                    good2=.false.  c                                    good2=.false.
1741  c                                    goto 880 !fill ntp and go to next event  c                                    goto 880 !fill ntp and go to next event
1742                                        do iv=1,nviews
1743                                           mask_view(iv) = 4
1744                                        enddo
1745                                      iflag=1                                      iflag=1
1746                                      return                                      return
1747                                   endif                                   endif
# Line 2394  c     print*,alfaxz1(ntrpt),alfaxz2(ntrp Line 1780  c     print*,alfaxz1(ntrpt),alfaxz2(ntrp
1780                                endif                                endif
1781                             enddo !end loop on COPPIA 3                             enddo !end loop on COPPIA 3
1782                          enddo   !end loop on sensors - COPPIA 3                          enddo   !end loop on sensors - COPPIA 3
1783     30                     continue
1784                       enddo      !end loop on planes  - COPPIA 3                       enddo      !end loop on planes  - COPPIA 3
1785   30                  continue   31                  continue
1786                                            
1787   1                enddo         !end loop on COPPIA 2   1                enddo         !end loop on COPPIA 2
1788                 enddo            !end loop on sensors - COPPIA 2                 enddo            !end loop on sensors - COPPIA 2
1789     20            continue
1790              enddo               !end loop on planes  - COPPIA 2              enddo               !end loop on planes  - COPPIA 2
1791                            
1792           enddo                  !end loop on COPPIA1           enddo                  !end loop on COPPIA1
1793        enddo                     !end loop on sensors - COPPIA 1        enddo                     !end loop on sensors - COPPIA 1
1794     10   continue
1795        enddo                     !end loop on planes  - COPPIA 1        enddo                     !end loop on planes  - COPPIA 1
1796                
1797        if(DEBUG)then        if(DEBUG)then
# Line 2584  c         if(ncpused.lt.ncpyz_min)goto 2 Line 1973  c         if(ncpused.lt.ncpyz_min)goto 2
1973       $           ,'( ',ncloyz_max,' )'       $           ,'( ',ncloyz_max,' )'
1974  c               good2=.false.  c               good2=.false.
1975  c     goto 880         !fill ntp and go to next event  c     goto 880         !fill ntp and go to next event
1976                do iv=1,nviews
1977                   mask_view(iv) = 5
1978                enddo
1979              iflag=1              iflag=1
1980              return              return
1981           endif           endif
# Line 2803  c         if(ncpused.lt.ncpxz_min)goto 2 Line 2195  c         if(ncpused.lt.ncpxz_min)goto 2
2195       $           ,'( ',ncloxz_max,' )'       $           ,'( ',ncloxz_max,' )'
2196  c     good2=.false.  c     good2=.false.
2197  c     goto 880         !fill ntp and go to next event  c     goto 880         !fill ntp and go to next event
2198                do iv=1,nviews
2199                   mask_view(iv) = 6
2200                enddo
2201              iflag=1              iflag=1
2202              return              return
2203           endif           endif
# Line 2900  c      include 'momanhough_init.f' Line 2295  c      include 'momanhough_init.f'
2295  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2296  *     variables for track fitting  *     variables for track fitting
2297        double precision AL_INI(5)        double precision AL_INI(5)
2298        double precision tath  c      double precision tath
2299  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2300  c      real fitz(nplanes)        !z coordinates of the planes in cm  c      real fitz(nplanes)        !z coordinates of the planes in cm
2301    
# Line 2966  c      real fitz(nplanes)        !z coor Line 2361  c      real fitz(nplanes)        !z coor
2361                 nplused=nplused+ hit_plane(ip)                 nplused=nplused+ hit_plane(ip)
2362              enddo              enddo
2363                            
2364              if(nplused.lt.nplxz_min)goto 888 !next doublet  c            if(nplused.lt.nplxz_min)goto 888 !next doublet
2365                if(nplused.lt.nplyz_min)goto 888 !next doublet
2366              if(ncp_ok.lt.ncpok)goto 888 !next cloud              if(ncp_ok.lt.ncpok)goto 888 !next cloud
2367                            
2368              if(DEBUG)then              if(DEBUG)then
# Line 2993  c$$$  print*,'6 -- ',(cly(6,i),i=1,ncp_p Line 2389  c$$$  print*,'6 -- ',(cly(6,i),i=1,ncp_p
2389  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'
2390                            
2391  *     -------> INITIAL GUESS <-------  *     -------> INITIAL GUESS <-------
2392              AL_INI(1) = dreal(alfaxz1_av(ixz))  cccc       SBAGLIATO
2393              AL_INI(2) = dreal(alfayz1_av(iyz))  c$$$            AL_INI(1) = dreal(alfaxz1_av(ixz))
2394              AL_INI(4) = PIGR + datan(dreal(alfayz2_av(iyz))  c$$$            AL_INI(2) = dreal(alfayz1_av(iyz))
2395       $           /dreal(alfaxz2_av(ixz)))  c$$$            AL_INI(4) = PIGR + datan(dreal(alfayz2_av(iyz))
2396              tath      = -dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))  c$$$     $           /dreal(alfaxz2_av(ixz)))
2397              AL_INI(3) = tath/sqrt(1+tath**2)  c$$$            tath      = -dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))
2398              AL_INI(5) = (1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.  c$$$            AL_INI(3) = tath/sqrt(1+tath**2)
2399                c$$$            AL_INI(5) = (1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.
2400  c     print*,'*******',AL_INI(5)  cccc       GIUSTO (ma si sua guess())
2401              if(AL_INI(5).gt.defmax)goto 888 !next cloud  c$$$            AL_INI(1) = dreal(alfaxz1_av(ixz))
2402                c$$$            AL_INI(2) = dreal(alfayz1_av(iyz))
2403  c     print*,'alfaxz2, alfayz2 '  c$$$            tath      = -dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))
2404  c     $              ,alfaxz2_av(ixz),alfayz2_av(iyz)  c$$$            AL_INI(3) = tath/sqrt(1+tath**2)
2405                c$$$            IF(alfaxz2_av(ixz).NE.0)THEN
2406  *     -------> INITIAL GUESS <-------  c$$$            AL_INI(4) = PIGR + datan(dreal(alfayz2_av(iyz))
2407  c     print*,'AL_INI ',(al_ini(i),i=1,5)  c$$$     $           /dreal(alfaxz2_av(ixz)))
2408                c$$$            ELSE
2409    c$$$               AL_INI(4) = acos(-1.)/2
2410    c$$$               IF(alfayz2_av(iyz).LT.0)AL_INI(4) = AL_INI(4)+acos(-1.)
2411    c$$$            ENDIF
2412    c$$$            IF(alfaxz2_av(ixz).LT.0)AL_INI(4)= acos(-1.)+ AL_INI(4)
2413    c$$$            AL_INI(4) = -acos(-1.) + AL_INI(4) !from incidence direction to tracking rs
2414    c$$$            
2415    c$$$            AL_INI(5) = (1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.
2416    c$$$            
2417    c$$$            if(AL_INI(5).gt.defmax)goto 888 !next cloud
2418                            
2419              if(DEBUG)then              if(DEBUG)then
2420                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))
2421                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))
# Line 3077  c     $                                 Line 2483  c     $                                
2483  *     **********************************************************  *     **********************************************************
2484  *     ************************** FIT *** FIT *** FIT *** FIT ***  *     ************************** FIT *** FIT *** FIT *** FIT ***
2485  *     **********************************************************  *     **********************************************************
2486    cccc  scommentare se si usa al_ini della nuvola
2487    c$$$                              do i=1,5
2488    c$$$                                 AL(i)=AL_INI(i)
2489    c$$$                              enddo
2490                                  call guess()
2491                                do i=1,5                                do i=1,5
2492                                   AL(i)=AL_INI(i)                                   AL_INI(i)=AL(i)
2493                                enddo                                enddo
2494                                ifail=0 !error flag in chi^2 computation                                ifail=0 !error flag in chi^2 computation
2495                                jstep=0 !number of  minimization steps                                jstep=0 !number of  minimization steps
2496                                iprint=0                                iprint=0
2497                                if(DEBUG)iprint=1  c                              if(DEBUG)iprint=1
2498                                  if(DEBUG)iprint=2
2499                                call mini2(jstep,ifail,iprint)                                call mini2(jstep,ifail,iprint)
2500                                if(ifail.ne.0) then                                if(ifail.ne.0) then
2501                                   if(DEBUG)then                                   if(DEBUG)then
2502                                      print *,                                      print *,
2503       $                              '*** MINIMIZATION FAILURE *** '       $                              '*** MINIMIZATION FAILURE *** '
2504       $                              //'(mini2 in clouds_to_ctrack)'       $                              //'(clouds_to_ctrack)'
2505                                        print*,'initial guess: '
2506    
2507                                        print*,'AL_INI(1) = ',AL_INI(1)
2508                                        print*,'AL_INI(2) = ',AL_INI(2)
2509                                        print*,'AL_INI(3) = ',AL_INI(3)
2510                                        print*,'AL_INI(4) = ',AL_INI(4)
2511                                        print*,'AL_INI(5) = ',AL_INI(5)
2512                                   endif                                   endif
2513                                   chi2=-chi2  c                                 chi2=-chi2
2514                                endif                                endif
2515  *     **********************************************************  *     **********************************************************
2516  *     ************************** FIT *** FIT *** FIT *** FIT ***  *     ************************** FIT *** FIT *** FIT *** FIT ***
# Line 3110  c     $                                 Line 2529  c     $                                
2529       $                ,'( ',NTRACKSMAX,' )'       $                ,'( ',NTRACKSMAX,' )'
2530  c                                 good2=.false.  c                                 good2=.false.
2531  c                                 goto 880 !fill ntp and go to next event                      c                                 goto 880 !fill ntp and go to next event                    
2532                                     do iv=1,nviews
2533                                        mask_view(iv) = 7
2534                                     enddo
2535                                   iflag=1                                   iflag=1
2536                                   return                                   return
2537                                endif                                endif
# Line 3230  c      include 'level1.f' Line 2652  c      include 'level1.f'
2652        call track_init        call track_init
2653        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
2654    
2655    *     |||||||||||||||||||||||||||||||||||||||||||||||||
2656  *     -------------------------------------------------  *     -------------------------------------------------
2657  *     If the plane has been already included, it just  *     If the plane has been already included, it just
2658  *     computes again the coordinates of the x-y couple  *     computes again the coordinates of the x-y couple
2659  *     using improved PFAs  *     using improved PFAs
2660  *     -------------------------------------------------  *     -------------------------------------------------
2661    *     |||||||||||||||||||||||||||||||||||||||||||||||||
2662           if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and.           if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and.
2663       $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then       $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then
2664                            
# Line 3269  c            dedxtrk(nplanes-ip+1) = (de Line 2693  c            dedxtrk(nplanes-ip+1) = (de
2693              dedxtrk_x(nplanes-ip+1)=dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)              dedxtrk_x(nplanes-ip+1)=dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
2694              dedxtrk_y(nplanes-ip+1)=dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)              dedxtrk_y(nplanes-ip+1)=dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
2695                            
2696    *     |||||||||||||||||||||||||||||||||||||||||||||||||
2697  *     -------------------------------------------------  *     -------------------------------------------------
2698  *     If the plane has NOT  been already included,  *     If the plane has NOT  been already included,
2699  *     it tries to include a COUPLE or a single cluster  *     it tries to include a COUPLE or a single cluster
2700  *     -------------------------------------------------  *     -------------------------------------------------
2701    *     |||||||||||||||||||||||||||||||||||||||||||||||||
2702           else                             else                  
2703                                
2704              xgood(nplanes-ip+1)=0              xgood(nplanes-ip+1)=0
# Line 3327  c     $              'ETA2','ETA2', Line 2753  c     $              'ETA2','ETA2',
2753       $              AYV_STORE(nplanes-ip+1,ibest))       $              AYV_STORE(nplanes-ip+1,ibest))
2754                                
2755                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
2756                   distance = distance / RCHI2_STORE(ibest)!<<< MS
2757                 id=id_cp(ip,icp,ist)                 id=id_cp(ip,icp,ist)
2758                 if(DEBUG)print*,'( couple ',id                 if(DEBUG)print*,'( couple ',id
2759       $              ,' ) normalized distance ',distance       $              ,' ) normalized distance ',distance
# Line 3397  c     $              'ETA2','ETA2', Line 2824  c     $              'ETA2','ETA2',
2824       $              PFA,PFA,       $              PFA,PFA,
2825       $              AXV_STORE(nplanes-ip+1,ibest),0.)                     $              AXV_STORE(nplanes-ip+1,ibest),0.)              
2826                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
2827  c     if(DEBUG)print*,'normalized distance ',distance                 distance = distance / RCHI2_STORE(ibest)!<<< MS
2828                 if(DEBUG)print*,'( cl-X ',icx                 if(DEBUG)print*,'( cl-X ',icx
2829       $              ,' in cp ',id,' ) normalized distance ',distance       $              ,' in cp ',id,' ) normalized distance ',distance
2830                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
# Line 3425  c     $              'ETA2','ETA2', Line 2852  c     $              'ETA2','ETA2',
2852       $              PFA,PFA,       $              PFA,PFA,
2853       $              0.,AYV_STORE(nplanes-ip+1,ibest))       $              0.,AYV_STORE(nplanes-ip+1,ibest))
2854                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
2855                   distance = distance / RCHI2_STORE(ibest)!<<< MS
2856                 if(DEBUG)print*,'( cl-Y ',icy                 if(DEBUG)print*,'( cl-Y ',icy
2857       $              ,' in cp ',id,' ) normalized distance ',distance       $              ,' in cp ',id,' ) normalized distance ',distance
2858                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
# Line 3444  c                 dedxmm = dedx(icy)  !( Line 2872  c                 dedxmm = dedx(icy)  !(
2872                 endif                                   endif                  
2873  11882          continue  11882          continue
2874              enddo               !end loop on cluster inside couples              enddo               !end loop on cluster inside couples
2875  *----- single clusters -----------------------------------------------      *----- single clusters -----------------------------------------------  
2876    c            print*,'## ncls(',ip,') ',ncls(ip)
2877              do ic=1,ncls(ip)    !loop on single clusters              do ic=1,ncls(ip)    !loop on single clusters
2878                 icl=cls(ip,ic)                 icl=cls(ip,ic)
2879    c              print*,'## ic ',ic,' ist ',ist
2880  c               if(cl_used(icl).eq.1.or.     !if the cluster is already used  c               if(cl_used(icl).eq.1.or.     !if the cluster is already used
2881                 if(cl_used(icl).ne.0.or.     !if the cluster is already used !(3)                 if(cl_used(icl).ne.0.or.     !if the cluster is already used !(3)
2882       $              LADDER(icl).ne.nldt.or. !or the ladder number does not match       $              LADDER(icl).ne.nldt.or. !or the ladder number does not match
# Line 3464  c     $                 'ETA2','ETA2', Line 2894  c     $                 'ETA2','ETA2',
2894                 endif                 endif
2895    
2896                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
2897                   distance = distance / RCHI2_STORE(ibest)!<<< MS
2898                 if(DEBUG)print*,'( cl-s ',icl                 if(DEBUG)print*,'( cl-s ',icl
2899       $              ,' ) normalized distance ',distance       $              ,' ) normalized distance ',distance,'<',distmin,' ?'
2900                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
2901                      if(DEBUG)print*,'YES'
2902                    xmm_A = xPAM_A                    xmm_A = xPAM_A
2903                    ymm_A = yPAM_A                    ymm_A = yPAM_A
2904                    zmm_A = zPAM_A                    zmm_A = zPAM_A
# Line 3488  c                  dedxmm = dedx(icl)   Line 2920  c                  dedxmm = dedx(icl)  
2920                 endif                                   endif                  
2921  18882          continue  18882          continue
2922              enddo               !end loop on single clusters              enddo               !end loop on single clusters
2923    c            print*,'## distmin ', distmin,' clinc ',clinc
2924              if(distmin.le.clinc)then                                if(distmin.le.clinc)then                  
2925                                
2926                 CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<                     CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<    
2927  *              ----------------------------  *              ----------------------------
2928    c               print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
2929                 if(mod(VIEW(iclm),2).eq.0)then                 if(mod(VIEW(iclm),2).eq.0)then
2930                    XGOOD(nplanes-ip+1)=1.                    XGOOD(nplanes-ip+1)=1.
2931                    resx(nplanes-ip+1)=rxmm                    resx(nplanes-ip+1)=rxmm
2932                    if(DEBUG)print*,'%%%% included X-cl ',iclm                    if(DEBUG)print*,'%%%% included X-cl ',iclm
2933       $                 ,' ( norm.dist.= ',distmin,', cut ',clinc,' )'  c                  if(.true.)print*,'%%%% included X-cl ',iclm
2934         $                 ,'( chi^2, ',RCHI2_STORE(ibest)
2935         $                 ,', norm.dist.= ',distmin
2936         $                 ,', cut ',clinc,' )'
2937                 else                 else
2938                    YGOOD(nplanes-ip+1)=1.                    YGOOD(nplanes-ip+1)=1.
2939                    resy(nplanes-ip+1)=rymm                    resy(nplanes-ip+1)=rymm
2940                    if(DEBUG)print*,'%%%% included Y-cl ',iclm                    if(DEBUG)print*,'%%%% included Y-cl ',iclm
2941       $                 ,' ( norm.dist.= ',distmin,', cut ',clinc,' )'  c                  if(.true.)print*,'%%%% included Y-cl ',iclm
2942         $                 ,'( chi^2, ',RCHI2_STORE(ibest)
2943         $                 ,', norm.dist.= ', distmin
2944         $                 ,', cut ',clinc,' )'
2945                 endif                 endif
2946    c               print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
2947  *              ----------------------------  *              ----------------------------
2948                 xm_A(nplanes-ip+1) = xmm_A                 xm_A(nplanes-ip+1) = xmm_A
2949                 ym_A(nplanes-ip+1) = ymm_A                 ym_A(nplanes-ip+1) = ymm_A
# Line 3616  c               endif Line 3056  c               endif
3056    
3057    
3058    
 c$$$*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  
 c$$$      real function fbad_cog(ncog,ic)  
 c$$$  
 c$$$  
 c$$$      include 'commontracker.f'  
 c$$$      include 'level1.f'  
 c$$$      include 'calib.f'  
 c$$$  
 c$$$*     --> signal of the central strip  
 c$$$      sc = CLSIGNAL(INDMAX(ic)) !center  
 c$$$  
 c$$$*     signal of adjacent strips  
 c$$$*     --> left  
 c$$$      sl1 = 0                  !left 1  
 c$$$      if(  
 c$$$     $     (INDMAX(ic)-1).ge.INDSTART(ic)  
 c$$$     $     )  
 c$$$     $     sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))  
 c$$$  
 c$$$      sl2 = 0                  !left 2  
 c$$$      if(  
 c$$$     $     (INDMAX(ic)-2).ge.INDSTART(ic)  
 c$$$     $     )  
 c$$$     $     sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))  
 c$$$  
 c$$$*     --> right  
 c$$$      sr1 = 0                  !right 1  
 c$$$      if(  
 c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))  
 c$$$     $     .or.  
 c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)  
 c$$$     $     )  
 c$$$     $     sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))  
 c$$$  
 c$$$      sr2 = 0                  !right 2  
 c$$$      if(  
 c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))  
 c$$$     $     .or.  
 c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)  
 c$$$     $     )  
 c$$$     $     sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))  
 c$$$  
 c$$$  
 c$$$      if(mod(int(VIEW(ic)),2).eq.1)then !Y-view  
 c$$$         f  = 4.  
 c$$$         si = 8.4  
 c$$$      else                              !X-view  
 c$$$         f  = 6.  
 c$$$         si = 3.9  
 c$$$      endif  
 c$$$  
 c$$$      fbad_cog = 1.  
 c$$$      f0 = 1  
 c$$$      f1 = 1  
 c$$$      f2 = 1  
 c$$$      f3 = 1    
 c$$$      if(sl1.gt.sr1.and.sl1.gt.0.)then  
 c$$$          
 c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f  
 c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f1=f  
 c$$$c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f3=f  
 c$$$  
 c$$$         if(ncog.eq.2.and.sl1.ne.0)then  
 c$$$            fbad_cog = (f1**2*sc**2/sl1**2+f0**2)/(sc**2/sl1**2+1.)  
 c$$$         elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then  
 c$$$            fbad_cog = 1.  
 c$$$         elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then  
 c$$$            fbad_cog = 1.  
 c$$$         else  
 c$$$            fbad_cog = 1.  
 c$$$         endif  
 c$$$          
 c$$$      elseif(sl1.le.sr1.and.sr1.gt.0.)then  
 c$$$  
 c$$$  
 c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f  
 c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f1=f  
 c$$$c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f3=f  
 c$$$  
 c$$$         if(ncog.eq.2.and.sr1.ne.0)then  
 c$$$            fbad_cog = (f1**2*sc**2/sr1**2+f0**2)/(sc**2/sr1**2+1.)  
 c$$$         elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then  
 c$$$            fbad_cog = 1.  
 c$$$         elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then  
 c$$$            fbad_cog = 1.  
 c$$$         else  
 c$$$            fbad_cog = 1.  
 c$$$         endif  
 c$$$  
 c$$$      endif  
 c$$$  
 c$$$      fbad_cog = sqrt(fbad_cog)  
 c$$$  
 c$$$      return  
 c$$$      end  
 c$$$  
   
3059    
3060    
3061  *     ****************************************************  *     ****************************************************
# Line 3797  c      include 'level1.f' Line 3140  c      include 'level1.f'
3140           alfayz1_nt(idb)=0                 alfayz1_nt(idb)=0      
3141           alfayz2_nt(idb)=0                 alfayz2_nt(idb)=0      
3142        enddo        enddo
3143        do itr=1,ntrpl_max_nt        do itr=1,ntrpt_max_nt
3144           tr_cloud_nt(itr)=0           tr_cloud_nt(itr)=0
3145           alfaxz1_nt(itr)=0                 alfaxz1_nt(itr)=0      
3146           alfaxz2_nt(itr)=0                 alfaxz2_nt(itr)=0      
# Line 3826  c      include 'level1.f' Line 3169  c      include 'level1.f'
3169          alfayz1(idb)=0                    alfayz1(idb)=0          
3170          alfayz2(idb)=0                    alfayz2(idb)=0          
3171        enddo                            enddo                    
3172        do itr=1,ntrpl_max                do itr=1,ntrpt_max        
3173          tr_cloud(itr)=0                  tr_cloud(itr)=0        
3174          cpxz1(itr)=0                      cpxz1(itr)=0            
3175          cpxz2(itr)=0                      cpxz2(itr)=0            
Legend:
Removed from v.1.9  
changed lines
  Added in v.1.17
  ViewVC Help
Powered by ViewVC 1.1.23