/[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.10 by pam-fi, Thu Nov 2 11:19:47 2006 UTC revision 1.23 by pam-fi, Mon May 14 11:03:06 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'  
24    
25    c      print*,'======================================================'
26    c$$$      do ic=1,NCLSTR1
27    c$$$         if(.false.
28    c$$$     $        .or.nsatstrips(ic).gt.0
29    c$$$c     $        .or.nbadstrips(0,ic).gt.0
30    c$$$c     $        .or.nbadstrips(4,ic).gt.0
31    c$$$c     $        .or.nbadstrips(3,ic).gt.0
32    c$$$     $        .or..false.)then
33    c$$$            print*,'--- cl-',ic,' ------------------------'
34    c$$$            istart = INDSTART(IC)
35    c$$$            istop  = TOTCLLENGTH
36    c$$$            if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
37    c$$$            print*,'ADC   ',(CLADC(i),i=istart,istop)
38    c$$$            print*,'s/n   ',(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
39    c$$$            print*,'sgnl  ',(CLSIGNAL(i),i=istart,istop)
40    c$$$            print*,'strip ',(i-INDMAX(ic),i=istart,istop)
41    c$$$            print*,'view ',VIEW(ic)
42    c$$$            print*,'maxs ',MAXS(ic)
43    c$$$            print*,'COG4 ',cog(4,ic)
44    c$$$            ff = fbad_cog(4,ic)
45    c$$$            print*,'fbad ',ff
46    c$$$            print*,(CLBAD(i),i=istart,istop)
47    c$$$            bb=nbadstrips(0,ic)
48    c$$$            print*,'#BAD (tot)',bb
49    c$$$            bb=nbadstrips(4,ic)
50    c$$$            print*,'#BAD (4)',bb
51    c$$$            bb=nbadstrips(3,ic)
52    c$$$            print*,'#BAD (3)',bb
53    c$$$            ss=nsatstrips(ic)
54    c$$$            print*,'#saturated ',ss
55    c$$$         endif
56    c$$$      enddo
57                
58  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
59  *     STEP 1  *     STEP 1
# Line 42  c      include 'momanhough_init.f' Line 74  c      include 'momanhough_init.f'
74  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
75  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
76    
77  c      iflag=0  
78        call cl_to_couples(iflag)        call cl_to_couples(iflag)
79        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
80           goto 880               !go to next event           goto 880               !go to next event
# Line 76  c      iflag=0 Line 108  c      iflag=0
108  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
109  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
110    
111  c      iflag=0  
112        call cp_to_doubtrip(iflag)        call cp_to_doubtrip(iflag)
113        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
114           goto 880               !go to next event                       goto 880               !go to next event            
# Line 183  c      iflag=0 Line 215  c      iflag=0
215  c      include 'momanhough_init.f'  c      include 'momanhough_init.f'
216                
217        logical FIMAGE            !        logical FIMAGE            !
218          real*8 AL_GUESS(5)
219    
220  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
221  *     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 266  c$$$            endif
266  c$$$         enddo  c$$$         enddo
267  c$$$         if(ibest.eq.0)goto 880 !>> no good candidates  c$$$         if(ibest.eq.0)goto 880 !>> no good candidates
268    
269           rchi2best=1000000000.  *     -------------------------------------------------------
270           ndofbest=0             !(1)  *     order track-candidates according to:
271    *     1st) decreasing n.points
272    *     2nd) increasing chi**2
273    *     -------------------------------------------------------
274             rchi2best=1000000000.
275             ndofbest=0            
276           do i=1,ntracks           do i=1,ntracks
277             if(RCHI2_STORE(i).lt.rchi2best.and.             ndof=0              
278       $          RCHI2_STORE(i).gt.0)then             do ii=1,nplanes    
279               ndof=0             !(1)               ndof=ndof        
280               do ii=1,nplanes    !(1)       $            +int(xgood_store(ii,i))
281                 ndof=ndof        !(1)       $            +int(ygood_store(ii,i))
282       $              +int(xgood_store(ii,i)) !(1)             enddo              
283       $              +int(ygood_store(ii,i)) !(1)             if(ndof.gt.ndofbest)then
284               enddo              !(1)               ibest=i
285               if(ndof.ge.ndofbest)then !(1)               rchi2best=RCHI2_STORE(i)
286                 ndofbest=ndof    
287               elseif(ndof.eq.ndofbest)then
288                 if(RCHI2_STORE(i).lt.rchi2best.and.
289         $            RCHI2_STORE(i).gt.0)then
290                 ibest=i                 ibest=i
291                 rchi2best=RCHI2_STORE(i)                 rchi2best=RCHI2_STORE(i)
292                 ndofbest=ndof    !(1)                 ndofbest=ndof  
293               endif              !(1)               endif            
294             endif             endif
295           enddo           enddo
296    
297    c$$$         rchi2best=1000000000.
298    c$$$         ndofbest=0             !(1)
299    c$$$         do i=1,ntracks
300    c$$$           if(RCHI2_STORE(i).lt.rchi2best.and.
301    c$$$     $          RCHI2_STORE(i).gt.0)then
302    c$$$             ndof=0             !(1)
303    c$$$             do ii=1,nplanes    !(1)
304    c$$$               ndof=ndof        !(1)
305    c$$$     $              +int(xgood_store(ii,i)) !(1)
306    c$$$     $              +int(ygood_store(ii,i)) !(1)
307    c$$$             enddo              !(1)
308    c$$$             if(ndof.ge.ndofbest)then !(1)
309    c$$$               ibest=i
310    c$$$               rchi2best=RCHI2_STORE(i)
311    c$$$               ndofbest=ndof    !(1)
312    c$$$             endif              !(1)
313    c$$$           endif
314    c$$$         enddo
315    
316           if(ibest.eq.0)goto 880 !>> no good candidates           if(ibest.eq.0)goto 880 !>> no good candidates
317  *-------------------------------------------------------------------------------      *-------------------------------------------------------------------------------    
318  *     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 269  c$$$         if(ibest.eq.0)goto 880 !>> Line 331  c$$$         if(ibest.eq.0)goto 880 !>>
331              iimage=0              iimage=0
332           endif           endif
333           if(icand.eq.0)then           if(icand.eq.0)then
334              print*,'HAI FATTO UN CASINO!!!!!! icand = ',icand              if(VERBOSE)then
335       $           ,ibest,iimage                 print*,'HAI FATTO UN CASINO!!!!!! icand = ',icand
336         $              ,ibest,iimage
337                endif
338              return              return
339           endif           endif
340    
# Line 281  c$$$         if(ibest.eq.0)goto 880 !>> Line 345  c$$$         if(ibest.eq.0)goto 880 !>>
345  *     **********************************************************  *     **********************************************************
346  *     ************************** FIT *** FIT *** FIT *** FIT ***  *     ************************** FIT *** FIT *** FIT *** FIT ***
347  *     **********************************************************  *     **********************************************************
348             call guess()
349             do i=1,5
350                AL_GUESS(i)=AL(i)
351             enddo
352    c         print*,'## guess: ',al
353    
354           do i=1,5           do i=1,5
355              AL(i)=dble(AL_STORE(i,icand))                          AL(i)=dble(AL_STORE(i,icand))            
356           enddo           enddo
357            
358           IDCAND = icand         !fitted track-candidate           IDCAND = icand         !fitted track-candidate
359           ifail=0                !error flag in chi2 computation           ifail=0                !error flag in chi2 computation
360           jstep=0                !# minimization steps           jstep=0                !# minimization steps
361    
362           iprint=0           iprint=0
363           if(DEBUG)iprint=1  c         if(DEBUG)iprint=1
364             if(VERBOSE)iprint=1
365             if(DEBUG)iprint=2
366           call mini2(jstep,ifail,iprint)           call mini2(jstep,ifail,iprint)
367           if(ifail.ne.0) then           if(ifail.ne.0) then
368              if(DEBUG)then              if(VERBOSE)then
369                 print *,                 print *,
370       $              '*** MINIMIZATION FAILURE *** (mini2) '       $              '*** MINIMIZATION FAILURE *** (after refinement) '
371       $              ,iev       $              ,iev
372    
373    c$$$               print*,'guess:   ',(al_guess(i),i=1,5)
374    c$$$               print*,'previous: ',(al_store(i,icand),i=1,5)
375    c$$$               print*,'result:   ',(al(i),i=1,5)
376    c$$$               print*,'xgood ',xgood
377    c$$$               print*,'ygood ',ygood
378    c$$$               print*,'----------------------------------------------'
379              endif              endif
380              chi2=-chi2  c            chi2=-chi2
381           endif           endif
382                    
383           if(DEBUG)then           if(DEBUG)then
# Line 408  c     $        rchi2best.lt.15..and. Line 488  c     $        rchi2best.lt.15..and.
488        end        end
489    
490    
   
   
 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$$$  
   
491                
492  ************************************************************  ************************************************************
493  ************************************************************  ************************************************************
# Line 604  c$$$ Line 512  c$$$
512  *     PFAy   - Position Finding Algorithm in y (COG2,ETA2,...)  *     PFAy   - Position Finding Algorithm in y (COG2,ETA2,...)
513  *     angx   - Projected angle in x  *     angx   - Projected angle in x
514  *     angy   - Projected angle in y  *     angy   - Projected angle in y
515    *     bfx    - x-component of magnetci field
516    *     bfy    - y-component of magnetic field
517  *  *
518  *     --------- COUPLES -------------------------------------------------------  *     --------- COUPLES -------------------------------------------------------
519  *     The couple defines a point in the space.  *     The couple defines a point in the space.
# Line 642  c$$$ Line 552  c$$$
552  *  *
553  *  *
554    
555        subroutine xyz_PAM(icx,icy,sensor,PFAx,PFAy,angx,angy)        subroutine xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
556    
 c*****************************************************  
 c     07/10/2005 modified by elena vannuccini --> (1)  
 c     01/02/2006 modified by elena vannuccini --> (2)  
 c     02/02/2006 modified by Elena Vannuccini --> (3)  
 c                (implemented new p.f.a.)  
 c     03/02/2006 modified by Elena Vannuccini --> (4)  
 c                (implemented variable resolution)  
 c*****************************************************  
557                
558        include 'commontracker.f'        include 'commontracker.f'
559        include 'level1.f'        include 'level1.f'
560        include 'calib.f'        include 'calib.f'
 c      include 'level1.f'  
561        include 'common_align.f'        include 'common_align.f'
562        include 'common_mech.f'        include 'common_mech.f'
563        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
       include 'common_resxy.f'  
   
 c      logical DEBUG  
 c      common/dbg/DEBUG  
564    
565        integer icx,icy           !X-Y cluster ID        integer icx,icy           !X-Y cluster ID
566        integer sensor        integer sensor
567        integer viewx,viewy        integer viewx,viewy
568        character*4 PFAx,PFAy     !PFA to be used        character*4 PFAx,PFAy     !PFA to be used
569        real angx,angy            !X-Y angle        real ax,ay                !X-Y geometric angle
570          real angx,angy            !X-Y effective angle
571          real bfx,bfy              !X-Y b-field components
572    
573        real stripx,stripy        real stripx,stripy
574    
575        double precision xrt,yrt,zrt        double precision xrt,yrt,zrt
576        double precision xrt_A,yrt_A,zrt_A        double precision xrt_A,yrt_A,zrt_A
577        double precision xrt_B,yrt_B,zrt_B        double precision xrt_B,yrt_B,zrt_B
 c      double precision xi,yi,zi  
 c      double precision xi_A,yi_A,zi_A  
 c      double precision xi_B,yi_B,zi_B  
578                
579    
580        parameter (ndivx=30)        parameter (ndivx=30)
581    
582    
583    c$$$      print*,icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy
584                
585        resxPAM = 0        resxPAM = 0
586        resyPAM = 0        resyPAM = 0
# Line 695  c      double precision xi_B,yi_B,zi_B Line 594  c      double precision xi_B,yi_B,zi_B
594        xPAM_B = 0.        xPAM_B = 0.
595        yPAM_B = 0.        yPAM_B = 0.
596        zPAM_B = 0.        zPAM_B = 0.
597    c      print*,'## xyz_PAM: ',icx,icy,sensor,PFAx,PFAy,angx,angy
598    
599  *     -----------------  *     -----------------
600  *     CLUSTER X  *     CLUSTER X
601  *     -----------------  *     -----------------
602    
603        if(icx.ne.0)then        if(icx.ne.0)then
604           viewx = VIEW(icx)  
605           nldx = nld(MAXS(icx),VIEW(icx))           viewx   = VIEW(icx)
606           nplx = npl(VIEW(icx))           nldx    = nld(MAXS(icx),VIEW(icx))
607           resxPAM = RESXAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!!           nplx    = npl(VIEW(icx))
608                     resxPAM = RESXAV
609           stripx = float(MAXS(icx))           stripx  = float(MAXS(icx))
610           if(PFAx.eq.'COG1')then  !(1)  *        --------------------------
611              stripx = stripx      !(1)  *        magnetic-field corrections
612              resxPAM = resxPAM    !(1)  *        --------------------------
613             angtemp  = ax
614             bfytemp  = bfy
615    *        /////////////////////////////////
616    *        AAAAHHHHHHHH!!!!!!!!!!!!!!!!!!!!!
617    *        *grvzkkjsdgjhhhgngbn###>:(
618    *        /////////////////////////////////
619    c         if(nplx.eq.6) angtemp = -1. * ax
620    c         if(nplx.eq.6) bfytemp = -1. * bfy
621             if(viewx.eq.12) angtemp = -1. * ax
622             if(viewx.eq.12) bfytemp = -1. * bfy
623             tgtemp   = tan(angtemp*acos(-1.)/180.) + pmuH_h*bfytemp*0.00001
624             angx     = 180.*atan(tgtemp)/acos(-1.)
625             stripx   = stripx - 0.5*pmuH_h*bfytemp*0.00001*SiDimZ/pitchX
626    c$$$         print*,nplx,ax,bfy/10.
627    c$$$         print*,angx,0.5*pmuH_h*bfytemp*0.00001*SiDimZ/pitchX
628    c$$$         print*,'========================'
629    c$$$         if(bfy.ne.0.)print*,viewx,'-x- '
630    c$$$     $        ,bfy,-1*0.5*pmuH_h*bfytemp*0.00001*SiDimZ
631    *        --------------------------
632    
633    c$$$         print*,'--- x-cl ---'
634    c$$$         istart = INDSTART(ICX)
635    c$$$         istop  = TOTCLLENGTH
636    c$$$         if(icx.lt.NCLSTR1)istop=INDSTART(ICX+1)-1
637    c$$$         print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
638    c$$$         print*,(CLSIGNAL(i),i=istart,istop)
639    c$$$         print*,INDMAX(icx)
640    c$$$         print*,cog(4,icx)
641    c$$$         print*,fbad_cog(4,icx)
642            
643    
644             if(PFAx.eq.'COG1')then
645    
646                stripx  = stripx
647                resxPAM = 1e-4*pitchX/sqrt(12.)!!resxPAM
648    
649           elseif(PFAx.eq.'COG2')then           elseif(PFAx.eq.'COG2')then
650              stripx = stripx + cog(2,icx)              
651                stripx  = stripx + cog(2,icx)            
652                resxPAM = risx_cog(abs(angx))!TEMPORANEO              
653              resxPAM = resxPAM*fbad_cog(2,icx)              resxPAM = resxPAM*fbad_cog(2,icx)
654    
655             elseif(PFAx.eq.'COG3')then
656    
657                stripx  = stripx + cog(3,icx)            
658                resxPAM = risx_cog(abs(angx))!TEMPORANEO                      
659                resxPAM = resxPAM*fbad_cog(3,icx)
660    
661             elseif(PFAx.eq.'COG4')then
662    
663                stripx  = stripx + cog(4,icx)            
664                resxPAM = risx_cog(abs(angx))!TEMPORANEO                      
665                resxPAM = resxPAM*fbad_cog(4,icx)
666    
667           elseif(PFAx.eq.'ETA2')then           elseif(PFAx.eq.'ETA2')then
668  c            cog2 = cog(2,icx)  
669  c            etacorr = pfaeta2(cog2,viewx,nldx,angx)                          stripx  = stripx + pfaeta2(icx,angx)          
670  c            stripx = stripx + etacorr              resxPAM = risx_eta2(abs(angx))
671              stripx = stripx + pfaeta2(icx,angx)            !(3)              resxPAM = resxPAM*fbad_cog(2,icx)
             resxPAM = risx_eta2(angx)                       !   (4)  
672              if(DEBUG.and.fbad_cog(2,icx).ne.1)              if(DEBUG.and.fbad_cog(2,icx).ne.1)
673       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)
674              resxPAM = resxPAM*fbad_cog(2,icx)  
675           elseif(PFAx.eq.'ETA3')then                         !(3)           elseif(PFAx.eq.'ETA3')then                        
676              stripx = stripx + pfaeta3(icx,angx)            !(3)  
677              resxPAM = risx_eta3(angx)                       !   (4)              stripx  = stripx + pfaeta3(icx,angx)          
678              if(DEBUG.and.fbad_cog(3,icx).ne.1)              !(3)              resxPAM = risx_eta3(abs(angx))                      
679       $           print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)!(3)              resxPAM = resxPAM*fbad_cog(3,icx)              
680              resxPAM = resxPAM*fbad_cog(3,icx)               !(3)              if(DEBUG.and.fbad_cog(3,icx).ne.1)            
681           elseif(PFAx.eq.'ETA4')then                         !(3)       $           print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)
682              stripx = stripx + pfaeta4(icx,angx)            !(3)  
683              resxPAM = risx_eta4(angx)                       !   (4)           elseif(PFAx.eq.'ETA4')then                        
684              if(DEBUG.and.fbad_cog(4,icx).ne.1)              !(3)  
685       $           print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)!(3)              stripx  = stripx + pfaeta4(icx,angx)            
686              resxPAM = resxPAM*fbad_cog(4,icx)               !(3)              resxPAM = risx_eta4(abs(angx))                      
687           elseif(PFAx.eq.'ETA')then                          !(3)              resxPAM = resxPAM*fbad_cog(4,icx)              
688              stripx = stripx + pfaeta(icx,angx)             !(3)              if(DEBUG.and.fbad_cog(4,icx).ne.1)              
689              resxPAM = ris_eta(icx,angx)                     !   (4)       $           print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)
690              if(DEBUG.and.fbad_cog(2,icx).ne.1)              !(3)  
691       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)!(3)           elseif(PFAx.eq.'ETA')then  
692  c            resxPAM = resxPAM*fbad_cog(2,icx)               !(3)TEMPORANEO  
693              resxPAM = resxPAM*fbad_eta(icx,angx)             !(3)(4)              stripx  = stripx + pfaeta(icx,angx)            
694           elseif(PFAx.eq.'COG')then           !(2)              resxPAM = ris_eta(icx,angx)                    
695              stripx = stripx + cog(0,icx)     !(2)                      resxPAM = resxPAM*fbad_eta(icx,angx)            
696              resxPAM = risx_cog(angx)                        !   (4)              if(DEBUG.and.fbad_cog(2,icx).ne.1)              
697              resxPAM = resxPAM*fbad_cog(0,icx)!(2)       $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)
698    
699             elseif(PFAx.eq.'COG')then          
700    
701                stripx  = stripx + cog(0,icx)            
702                resxPAM = risx_cog(abs(angx))                    
703                resxPAM = resxPAM*fbad_cog(0,icx)
704    
705           else           else
706              print*,'*** Non valid p.f.a. (x) --> ',PFAx              if(DEBUG) print*,'*** Non valid p.f.a. (x) --> ',PFAx
707             endif
708    
709    
710    *     ======================================
711    *     temporary patch for saturated clusters
712    *     ======================================
713             if( nsatstrips(icx).gt.0 )then
714                stripx  = stripx + cog(4,icx)            
715                resxPAM = pitchX*1e-4/sqrt(12.)
716                cc=cog(4,icx)
717    c$$$            print*,icx,' *** ',cc
718    c$$$            print*,icx,' *** ',resxPAM
719           endif           endif
720    
721        endif        endif
 c      if(icy.eq.0.and.icx.ne.0)  
 c     $     print*,PFAx,icx,angx,stripx,resxPAM,'***'  
722                
723  *     -----------------  *     -----------------
724  *     CLUSTER Y  *     CLUSTER Y
725  *     -----------------  *     -----------------
726    
727        if(icy.ne.0)then        if(icy.ne.0)then
728    
729           viewy = VIEW(icy)           viewy = VIEW(icy)
730           nldy = nld(MAXS(icy),VIEW(icy))           nldy = nld(MAXS(icy),VIEW(icy))
731           nply = npl(VIEW(icy))           nply = npl(VIEW(icy))
732           resyPAM = RESYAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!!           resyPAM = RESYAV
733             stripy = float(MAXS(icy))
734    
735           if(icx.ne.0.and.(nply.ne.nplx.or.nldy.ne.nldx))then           if(icx.ne.0.and.(nply.ne.nplx.or.nldy.ne.nldx))then
736              print*,'xyz_PAM   ***ERROR*** invalid cluster couple!!! '              if(DEBUG) then
737       $           ,icx,icy                 print*,'xyz_PAM   ***ERROR*** invalid cluster couple!!! '
738         $              ,icx,icy
739                endif
740              goto 100              goto 100
741           endif           endif
742            *        --------------------------
743           stripy = float(MAXS(icy))  *        magnetic-field corrections
744           if(PFAy.eq.'COG1')then !(1)  *        --------------------------
745              stripy = stripy     !(1)           tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001        
746              resyPAM = resyPAM   !(1)           angy    = 180.*atan(tgtemp)/acos(-1.)
747             stripy = stripy + 0.5*pmuH_e*bfx*0.00001*SiDimZ/pitchY
748    c$$$         if(bfx.ne.0.)print*,viewy,'-y- '
749    c$$$     $        ,bfx,0.5*pmuH_e*bfx*0.00001*SiDimZ
750    *        --------------------------
751            
752    c$$$         print*,'--- y-cl ---'
753    c$$$         istart = INDSTART(ICY)
754    c$$$         istop  = TOTCLLENGTH
755    c$$$         if(icy.lt.NCLSTR1)istop=INDSTART(ICY+1)-1
756    c$$$         print*,(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
757    c$$$         print*,(CLSIGNAL(i),i=istart,istop)
758    c$$$         print*,INDMAX(icy)
759    c$$$         print*,cog(4,icy)
760    c$$$         print*,fbad_cog(4,icy)
761    
762             if(PFAy.eq.'COG1')then
763    
764                stripy  = stripy    
765                resyPAM = 1e-4*pitchY/sqrt(12.)!resyPAM  
766    
767           elseif(PFAy.eq.'COG2')then           elseif(PFAy.eq.'COG2')then
768              stripy = stripy + cog(2,icy)  
769                stripy  = stripy + cog(2,icy)
770                resyPAM = risy_cog(abs(angy))!TEMPORANEO
771              resyPAM = resyPAM*fbad_cog(2,icy)              resyPAM = resyPAM*fbad_cog(2,icy)
772    
773             elseif(PFAy.eq.'COG3')then
774    
775                stripy  = stripy + cog(3,icy)
776                resyPAM = risy_cog(abs(angy))!TEMPORANEO
777                resyPAM = resyPAM*fbad_cog(3,icy)
778    
779             elseif(PFAy.eq.'COG4')then
780    
781                stripy  = stripy + cog(4,icy)
782                resyPAM = risy_cog(abs(angy))!TEMPORANEO
783                resyPAM = resyPAM*fbad_cog(4,icy)
784    
785           elseif(PFAy.eq.'ETA2')then           elseif(PFAy.eq.'ETA2')then
786  c            cog2 = cog(2,icy)  
787  c            etacorr = pfaeta2(cog2,viewy,nldy,angy)              stripy  = stripy + pfaeta2(icy,angy)          
788  c            stripy = stripy + etacorr              resyPAM = risy_eta2(abs(angy))              
             stripy = stripy + pfaeta2(icy,angy)            !(3)  
             resyPAM = risy_eta2(angy)                       !   (4)  
789              resyPAM = resyPAM*fbad_cog(2,icy)              resyPAM = resyPAM*fbad_cog(2,icy)
790              if(DEBUG.and.fbad_cog(2,icy).ne.1)              if(DEBUG.and.fbad_cog(2,icy).ne.1)
791       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)
792           elseif(PFAy.eq.'ETA3')then                         !(3)  
793              stripy = stripy + pfaeta3(icy,angy)            !(3)           elseif(PFAy.eq.'ETA3')then                      
794              resyPAM = resyPAM*fbad_cog(3,icy)               !(3)  
795              if(DEBUG.and.fbad_cog(3,icy).ne.1)              !(3)              stripy  = stripy + pfaeta3(icy,angy)
796       $           print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)!(3)              resyPAM = resyPAM*fbad_cog(3,icy)  
797           elseif(PFAy.eq.'ETA4')then                         !(3)              if(DEBUG.and.fbad_cog(3,icy).ne.1)
798              stripy = stripy + pfaeta4(icy,angy)            !(3)       $           print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)
799              resyPAM = resyPAM*fbad_cog(4,icy)               !(3)  
800              if(DEBUG.and.fbad_cog(4,icy).ne.1)              !(3)           elseif(PFAy.eq.'ETA4')then  
801       $           print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)!(3)  
802           elseif(PFAy.eq.'ETA')then                          !(3)              stripy  = stripy + pfaeta4(icy,angy)
803              stripy = stripy + pfaeta(icy,angy)             !(3)              resyPAM = resyPAM*fbad_cog(4,icy)
804              resyPAM = ris_eta(icy,angy)                     !   (4)              if(DEBUG.and.fbad_cog(4,icy).ne.1)
805  c            resyPAM = resyPAM*fbad_cog(2,icy)              !(3)TEMPORANEO       $           print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)
806              resyPAM = resyPAM*fbad_eta(icy,angy)            !   (4)  
807              if(DEBUG.and.fbad_cog(2,icy).ne.1)              !(3)           elseif(PFAy.eq.'ETA')then
808       $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)!(3)  
809                stripy  = stripy + pfaeta(icy,angy)
810                resyPAM = ris_eta(icy,angy)  
811                resyPAM = resyPAM*fbad_eta(icy,angy)
812                if(DEBUG.and.fbad_cog(2,icy).ne.1)
813         $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)
814    
815           elseif(PFAy.eq.'COG')then           elseif(PFAy.eq.'COG')then
816              stripy = stripy + cog(0,icy)              
817              resyPAM = risy_cog(angy)                        !   (4)              stripy  = stripy + cog(0,icy)            
818  c            resyPAM = ris_eta(icy,angy)                    !   (4)              resyPAM = risy_cog(abs(angy))
819              resyPAM = resyPAM*fbad_cog(0,icy)              resyPAM = resyPAM*fbad_cog(0,icy)
820    
821           else           else
822              print*,'*** Non valid p.f.a. (x) --> ',PFAx              if(DEBUG) print*,'*** Non valid p.f.a. (x) --> ',PFAx
823           endif           endif
824    
825    
826    *     ======================================
827    *     temporary patch for saturated clusters
828    *     ======================================
829             if( nsatstrips(icy).gt.0 )then
830                stripy  = stripy + cog(4,icy)            
831                resyPAM = pitchY*1e-4/sqrt(12.)
832                cc=cog(4,icy)
833    c$$$            print*,icy,' *** ',cc
834    c$$$            print*,icy,' *** ',resyPAM
835             endif
836    
837    
838        endif        endif
839    
840          c$$$      print*,'## stripx,stripy ',stripx,stripy
841    
842  c===========================================================  c===========================================================
843  C     COUPLE  C     COUPLE
844  C===========================================================  C===========================================================
# Line 825  c     (xi,yi,zi) = mechanical coordinate Line 849  c     (xi,yi,zi) = mechanical coordinate
849  c------------------------------------------------------------------------  c------------------------------------------------------------------------
850           if(((mod(int(stripx+0.5)-1,1024)+1).le.3)           if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
851       $        .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...       $        .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
852              print*,'xyz_PAM (couple):',              if(DEBUG) then
853       $          ' WARNING: false X strip: strip ',stripx                 print*,'xyz_PAM (couple):',
854         $              ' WARNING: false X strip: strip ',stripx
855                endif
856           endif           endif
857           xi = acoordsi(stripx,viewx)           xi = acoordsi(stripx,viewx)
858           yi = acoordsi(stripy,viewy)           yi = acoordsi(stripy,viewy)
# Line 918  c            print*,'X-singlet ',icx,npl Line 944  c            print*,'X-singlet ',icx,npl
944  c            if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips...  c            if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips...
945              if(((mod(int(stripx+0.5)-1,1024)+1).le.3)              if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
946       $           .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...       $           .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
947                 print*,'xyz_PAM (X-singlet):',                 if(DEBUG) then
948       $             ' WARNING: false X strip: strip ',stripx                    print*,'xyz_PAM (X-singlet):',
949         $                 ' WARNING: false X strip: strip ',stripx
950                   endif
951              endif              endif
952              xi   = acoordsi(stripx,viewx)              xi   = acoordsi(stripx,viewx)
953    
# Line 941  c            print*,'X-cl ',icx,stripx,' Line 969  c            print*,'X-cl ',icx,stripx,'
969  c            print*,yi_A,' <--> ',yi_B  c            print*,yi_A,' <--> ',yi_B
970    
971           else           else
972                if(DEBUG) then
973              print *,'routine xyz_PAM ---> not properly used !!!'                 print *,'routine xyz_PAM ---> not properly used !!!'
974              print *,'icx = ',icx                 print *,'icx = ',icx
975              print *,'icy = ',icy                 print *,'icy = ',icy
976                endif
977              goto 100              goto 100
978                            
979           endif           endif
# Line 1009  c--------------------------------------- Line 1038  c---------------------------------------
1038  c         print*,'A-(',xPAM_A,yPAM_A,') B-(',xPAM_B,yPAM_B,')'  c         print*,'A-(',xPAM_A,yPAM_A,') B-(',xPAM_B,yPAM_B,')'
1039    
1040        else        else
1041                       if(DEBUG) then
1042           print *,'routine xyz_PAM ---> not properly used !!!'              print *,'routine xyz_PAM ---> not properly used !!!'
1043           print *,'icx = ',icx              print *,'icx = ',icx
1044           print *,'icy = ',icy              print *,'icy = ',icy
1045                         endif
1046        endif        endif
1047                    
1048    
1049    c      print*,'## xPAM,yPAM,zPAM       ',xPAM,yPAM,zPAM
1050    c      print*,'## xPAM_A,yPAM_A,zPAM_A ',xPAM_A,yPAM_A,zPAM_A
1051    c      print*,'## xPAM_B,yPAM_B,zPAM_B ',xPAM_B,yPAM_B,zPAM_B
1052    
1053   100  continue   100  continue
1054        end        end
1055    
1056    ************************************************************************
1057    *     Call xyz_PAM subroutine with default PFA and fill the mini2 common.
1058    *     (done to be called from c/c++)
1059    ************************************************************************
1060    
1061          subroutine xyzpam(ip,icx,icy,lad,sensor,ax,ay,bfx,bfy)
1062    
1063          include 'commontracker.f'
1064          include 'level1.f'
1065          include 'common_mini_2.f'
1066          include 'common_xyzPAM.f'
1067          include 'common_mech.f'
1068          include 'calib.f'
1069          
1070    *     flag to chose PFA
1071    c$$$      character*10 PFA
1072    c$$$      common/FINALPFA/PFA
1073    
1074          integer icx,icy           !X-Y cluster ID
1075          integer sensor
1076          character*4 PFAx,PFAy     !PFA to be used
1077          real ax,ay                !X-Y geometric angle
1078          real bfx,bfy              !X-Y b-field components
1079    
1080          ipx=0
1081          ipy=0      
1082          
1083    c$$$      PFAx = 'COG4'!PFA
1084    c$$$      PFAy = 'COG4'!PFA
1085          
1086          call idtoc(pfaid,PFAx)
1087          call idtoc(pfaid,PFAy)
1088    
1089          call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
1090    
1091    c$$$      print*,icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy
1092          
1093          if(icx.ne.0.and.icy.ne.0)then
1094    
1095             ipx=npl(VIEW(icx))
1096             ipy=npl(VIEW(icy))
1097             if( (nplanes-ipx+1).ne.ip.or.(nplanes-ipy+1).ne.ip )
1098         $        print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1099         $        ,' does not belong to the correct plane: ',ip,ipx,ipy
1100    
1101             xgood(ip) = 1.
1102             ygood(ip) = 1.
1103             resx(ip)  = resxPAM
1104             resy(ip)  = resyPAM
1105    
1106             xm(ip) = xPAM
1107             ym(ip) = yPAM
1108             zm(ip) = zPAM
1109             xm_A(ip) = 0.
1110             ym_A(ip) = 0.
1111             xm_B(ip) = 0.
1112             ym_B(ip) = 0.
1113    
1114    c         zv(ip) = zPAM
1115    
1116          elseif(icx.eq.0.and.icy.ne.0)then
1117    
1118             ipy=npl(VIEW(icy))
1119             if((nplanes-ipy+1).ne.ip)
1120         $        print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1121         $        ,' does not belong to the correct plane: ',ip,ipx,ipy
1122    
1123             xgood(ip) = 0.
1124             ygood(ip) = 1.
1125             resx(ip)  = 1000.
1126             resy(ip)  = resyPAM
1127    
1128             xm(ip) = -100.
1129             ym(ip) = -100.
1130             zm(ip) = (zPAM_A+zPAM_B)/2.
1131             xm_A(ip) = xPAM_A
1132             ym_A(ip) = yPAM_A
1133             xm_B(ip) = xPAM_B
1134             ym_B(ip) = yPAM_B
1135    
1136    c         zv(ip) = (zPAM_A+zPAM_B)/2.
1137            
1138          elseif(icx.ne.0.and.icy.eq.0)then
1139    
1140             ipx=npl(VIEW(icx))
1141             if((nplanes-ipx+1).ne.ip)
1142         $        print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1143         $        ,' does not belong to the correct plane: ',ip,ipx,ipy
1144    
1145             xgood(ip) = 1.
1146             ygood(ip) = 0.
1147             resx(ip)  = resxPAM
1148             resy(ip)  = 1000.
1149    
1150             xm(ip) = -100.
1151             ym(ip) = -100.
1152             zm(ip) = (zPAM_A+zPAM_B)/2.
1153             xm_A(ip) = xPAM_A
1154             ym_A(ip) = yPAM_A
1155             xm_B(ip) = xPAM_B
1156             ym_B(ip) = yPAM_B
1157            
1158    c         zv(ip) = (zPAM_A+zPAM_B)/2.
1159    
1160          else
1161    
1162             il = 2
1163             if(lad.ne.0)il=lad
1164             is = 1
1165             if(sensor.ne.0)is=sensor
1166    c         print*,nplanes-ip+1,il,is
1167    
1168             xgood(ip) = 0.
1169             ygood(ip) = 0.
1170             resx(ip)  = 1000.
1171             resy(ip)  = 1000.
1172    
1173             xm(ip) = -100.
1174             ym(ip) = -100.          
1175             zm(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
1176             xm_A(ip) = 0.
1177             ym_A(ip) = 0.
1178             xm_B(ip) = 0.
1179             ym_B(ip) = 0.
1180    
1181    c         zv(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
1182    
1183          endif
1184    
1185          if(DEBUG)then
1186    c         print*,'----------------------------- track coord'
1187    22222    format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)
1188             write(*,22222)ip,zm(ip),xm(ip),ym(ip)
1189         $        ,xm_A(ip),ym_A(ip),xm_B(ip),ym_B(ip)
1190         $        ,xgood(ip),ygood(ip),resx(ip),resy(ip)
1191    c$$$         print*,'-----------------------------'
1192          endif
1193          end
1194    
1195  ********************************************************************************  ********************************************************************************
1196  ********************************************************************************  ********************************************************************************
# Line 1094  c         print*,'A-(',xPAM_A,yPAM_A,') Line 1266  c         print*,'A-(',xPAM_A,yPAM_A,')
1266           endif                   endif        
1267    
1268           distance=           distance=
1269       $        ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2       $       ((xmi-XPP)**2+(ymi-YPP)**2)!QUIQUI
1270    cc     $        ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2
1271           distance=dsqrt(distance)                               distance=dsqrt(distance)                    
1272    
1273  c$$$         print*,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b  c$$$         print*,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b
# Line 1119  c$$$         print*,' resolution ',re Line 1292  c$$$         print*,' resolution ',re
1292  *     ----------------------  *     ----------------------
1293                    
1294           distance=           distance=
1295       $        ((xPAM-XPP)/resxPAM)**2       $       ((xPAM-XPP))**2 !QUIQUI
1296       $        +       $       +
1297       $        ((yPAM-YPP)/resyPAM)**2       $       ((yPAM-YPP))**2
1298    c$$$     $        ((xPAM-XPP)/resxPAM)**2
1299    c$$$     $        +
1300    c$$$     $        ((yPAM-YPP)/resyPAM)**2
1301           distance=dsqrt(distance)                               distance=dsqrt(distance)                    
1302    
1303  c$$$         print*,xPAM,yPAM,zPAM  c$$$         print*,xPAM,yPAM,zPAM
# Line 1130  c$$$         print*,' resolution ',resxP Line 1306  c$$$         print*,' resolution ',resxP
1306                    
1307        else        else
1308                    
1309           print*  c         print*
1310       $        ,' function distance_to ---> wrong usage!!!'  c     $        ,' function distance_to ---> wrong usage!!!'
1311           print*,' xPAM,yPAM,zPAM ',xPAM,yPAM,zPAM  c         print*,' xPAM,yPAM,zPAM ',xPAM,yPAM,zPAM
1312           print*,' xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b '  c         print*,' xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b '
1313       $        ,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b  c     $        ,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b
1314        endif          endif  
1315    
1316        distance_to = sngl(distance)        distance_to = sngl(distance)
# Line 1202  c--------------------------------------- Line 1378  c---------------------------------------
1378                 if(((mod(int(stripx+0.5)-1,1024)+1).le.3)                 if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
1379       $              .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...       $              .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
1380  c     if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips...  c     if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips...
1381                    print*,'whichsensor: ',  c                  print*,'whichsensor: ',
1382       $                ' WARNING: false X strip: strip ',stripx  c     $                ' WARNING: false X strip: strip ',stripx
1383                 endif                 endif
1384                 xi = acoordsi(stripx,viewx)                 xi = acoordsi(stripx,viewx)
1385                 yi = acoordsi(stripy,viewy)                 yi = acoordsi(stripy,viewy)
# Line 1358  c      include 'common_analysis.f' Line 1534  c      include 'common_analysis.f'
1534        is_cp=0        is_cp=0
1535        if(id.lt.0)is_cp=1        if(id.lt.0)is_cp=1
1536        if(id.gt.0)is_cp=2        if(id.gt.0)is_cp=2
1537        if(id.eq.0)print*,'IS_CP ===> wrong couple id !!!'  c      if(id.eq.0)print*,'IS_CP ===> wrong couple id !!!'
1538    
1539        return        return
1540        end        end
# Line 1429  c      include 'common_analysis.f' Line 1605  c      include 'common_analysis.f'
1605  *************************************************************************  *************************************************************************
1606  *************************************************************************  *************************************************************************
1607  *************************************************************************  *************************************************************************
 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  
1608                
1609    
1610  ***************************************************  ***************************************************
# Line 1754  c      include 'level1.f' Line 1662  c      include 'level1.f'
1662        enddo        enddo
1663  *     mask views with too many clusters  *     mask views with too many clusters
1664        do iv=1,nviews        do iv=1,nviews
1665           if( ncl_view(iv).gt. nclustermax)then           if( ncl_view(iv).gt. nclusterlimit)then
1666              mask_view(iv) = 1              mask_view(iv) = 1
1667              if(VERBOSE)print*,' * WARNING * cl_to_couple: n.clusters > '              if(DEBUG)print*,' * WARNING * cl_to_couple: n.clusters > '
1668       $           ,nclustermax,' on view ', iv,' --> masked!'       $           ,nclusterlimit,' on view ', iv,' --> masked!'
1669           endif           endif
1670        enddo        enddo
1671    
# Line 1774  c      include 'level1.f' Line 1682  c      include 'level1.f'
1682  *     ----------------------------------------------------  *     ----------------------------------------------------
1683  *     cut on charge (X VIEW)  *     cut on charge (X VIEW)
1684  *     ----------------------------------------------------  *     ----------------------------------------------------
1685           if(dedx(icx).lt.dedx_x_min)then           if(sgnl(icx).lt.dedx_x_min)then
1686              cl_single(icx)=0              cl_single(icx)=0
1687              goto 10              goto 10
1688           endif           endif
# Line 1824  c     endif Line 1732  c     endif
1732  *     ----------------------------------------------------  *     ----------------------------------------------------
1733  *     cut on charge (Y VIEW)  *     cut on charge (Y VIEW)
1734  *     ----------------------------------------------------  *     ----------------------------------------------------
1735              if(dedx(icy).lt.dedx_y_min)then              if(sgnl(icy).lt.dedx_y_min)then
1736                 cl_single(icy)=0                 cl_single(icy)=0
1737                 goto 20                 goto 20
1738              endif              endif
# Line 1870  c     endif Line 1778  c     endif
1778  *     charge correlation  *     charge correlation
1779  *     (modified to be applied only below saturation... obviously)  *     (modified to be applied only below saturation... obviously)
1780    
1781                 if(  .not.(dedx(icy).gt.chsaty.and.dedx(icx).gt.chsatx)                 if(  .not.(sgnl(icy).gt.chsaty.and.sgnl(icx).gt.chsatx)
1782       $              .and.       $              .and.
1783       $              .not.(dedx(icy).lt.chmipy.and.dedx(icx).lt.chmipx)       $              .not.(sgnl(icy).lt.chmipy.and.sgnl(icx).lt.chmipx)
1784       $              .and.       $              .and.
1785       $              (badclx.eq.1.and.badcly.eq.1)       $              (badclx.eq.1.and.badcly.eq.1)
1786       $              .and.       $              .and.
1787       $              .true.)then       $              .true.)then
1788    
1789                    ddd=(dedx(icy)                    ddd=(sgnl(icy)
1790       $                 -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))       $                 -kch(nplx,nldx)*sgnl(icx)-cch(nplx,nldx))
1791                    ddd=ddd/sqrt(kch(nplx,nldx)**2+1)                    ddd=ddd/sqrt(kch(nplx,nldx)**2+1)
1792    
1793  c                  cut = chcut * sch(nplx,nldx)  c                  cut = chcut * sch(nplx,nldx)
1794    
1795                    sss=(kch(nplx,nldx)*dedx(icy)+dedx(icx)                    sss=(kch(nplx,nldx)*sgnl(icy)+sgnl(icx)
1796       $                 -kch(nplx,nldx)*cch(nplx,nldx))       $                 -kch(nplx,nldx)*cch(nplx,nldx))
1797                    sss=sss/sqrt(kch(nplx,nldx)**2+1)                    sss=sss/sqrt(kch(nplx,nldx)**2+1)
1798                    cut = chcut * (16 + sss/50.)                    cut = chcut * (16 + sss/50.)
# Line 1893  c                  cut = chcut * sch(npl Line 1801  c                  cut = chcut * sch(npl
1801                       goto 20    !charge not consistent                       goto 20    !charge not consistent
1802                    endif                    endif
1803                 endif                 endif
1804                  
1805  *     ------------------> 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  
1806                    if(verbose)print*,                    if(verbose)print*,
1807       $                 '** warning ** number of identified '//       $                 '** warning ** number of identified '//
1808       $                 'couples on plane ',nplx,       $                 'couples on plane ',nplx,
1809       $                 'exceeds vector dimention '       $                 'exceeds vector dimention '
1810       $                 ,'( ',ncouplemax,' )'       $                 ,'( ',ncouplemax,' ) --> masked!'
 c     good2=.false.  
 c     goto 880   !fill ntp and go to next event  
1811                    mask_view(nviewx(nplx)) = 2                    mask_view(nviewx(nplx)) = 2
1812                    mask_view(nviewy(nply)) = 2                    mask_view(nviewy(nply)) = 2
1813  c                  iflag=1                    goto 10
 c                  return  
1814                 endif                 endif
1815                                
1816    *     ------------------> COUPLE <------------------
1817                 ncp_plane(nplx) = ncp_plane(nplx) + 1                 ncp_plane(nplx) = ncp_plane(nplx) + 1
1818                 clx(nplx,ncp_plane(nplx))=icx                 clx(nplx,ncp_plane(nplx))=icx
1819                 cly(nply,ncp_plane(nplx))=icy                 cly(nply,ncp_plane(nplx))=icy
1820                 cl_single(icx)=0                 cl_single(icx)=0
1821                 cl_single(icy)=0                 cl_single(icy)=0
             endif                                
1822  *     ----------------------------------------------  *     ----------------------------------------------
1823    
1824                endif                              
1825    
1826   20         continue   20         continue
1827           enddo                  !end loop on clusters(Y)           enddo                  !end loop on clusters(Y)
1828                    
# Line 1956  c                  return Line 1849  c                  return
1849        do ip=1,6        do ip=1,6
1850           ncp_tot = ncp_tot + ncp_plane(ip)           ncp_tot = ncp_tot + ncp_plane(ip)
1851        enddo        enddo
 c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)  
         
 c$$$      if(ncp_tot.gt.ncp_max)then  
 c$$$         if(verbose)print*,  
 c$$$     $           '** warning ** number of identified '//  
 c$$$     $           'couples exceeds upper limit for Hough tr. '  
 c$$$     $           ,'( ',ncp_max,' )'              
 c$$$         iflag=1  
 c$$$         return  
 c$$$      endif  
1852                
1853        return        return
1854        end        end
# Line 1978  c$$$      endif Line 1861  c$$$      endif
1861  *                                                 *  *                                                 *
1862  *                                                 *  *                                                 *
1863  **************************************************  **************************************************
 c$$$      subroutine cl_to_couples_nocharge(iflag)  
 c$$$  
 c$$$      include 'commontracker.f'  
 c$$$      include 'level1.f'  
 c$$$      include 'common_momanhough.f'  
 c$$$c      include 'momanhough_init.f'  
 c$$$      include 'calib.f'  
 c$$$c      include 'level1.f'  
 c$$$  
 c$$$  
 c$$$*     output flag  
 c$$$*     --------------  
 c$$$*     0 = good event  
 c$$$*     1 = bad event  
 c$$$*     --------------  
 c$$$      integer iflag  
 c$$$  
 c$$$      integer badseed,badcl  
 c$$$  
 c$$$*     init variables  
 c$$$      ncp_tot=0  
 c$$$      do ip=1,nplanes  
 c$$$         do ico=1,ncouplemax  
 c$$$            clx(ip,ico)=0  
 c$$$            cly(ip,ico)=0  
 c$$$         enddo  
 c$$$         ncp_plane(ip)=0  
 c$$$         do icl=1,nclstrmax_level2  
 c$$$            cls(ip,icl)=1  
 c$$$         enddo  
 c$$$         ncls(ip)=0  
 c$$$      enddo  
 c$$$      do icl=1,nclstrmax_level2  
 c$$$         cl_single(icl)=1  
 c$$$         cl_good(icl)=0  
 c$$$      enddo  
 c$$$        
 c$$$*     start association  
 c$$$      ncouples=0  
 c$$$      do icx=1,nclstr1          !loop on cluster (X)  
 c$$$         if(mod(VIEW(icx),2).eq.1)goto 10  
 c$$$          
 c$$$*     ----------------------------------------------------  
 c$$$*     cut on charge (X VIEW)  
 c$$$         if(dedx(icx).lt.dedx_x_min)then  
 c$$$            cl_single(icx)=0  
 c$$$            goto 10  
 c$$$         endif  
 c$$$*     cut BAD (X VIEW)              
 c$$$         badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx)))  
 c$$$         ifirst=INDSTART(icx)  
 c$$$         if(icx.ne.nclstr1) then  
 c$$$            ilast=INDSTART(icx+1)-1  
 c$$$         else  
 c$$$            ilast=TOTCLLENGTH  
 c$$$         endif  
 c$$$         badcl=badseed  
 c$$$         do igood=-ngoodstr,ngoodstr  
 c$$$            ibad=1  
 c$$$            if((INDMAX(icx)+igood).gt.ifirst.and.  
 c$$$     $           (INDMAX(icx)+igood).lt.ilast.and.  
 c$$$     $           .true.)then  
 c$$$               ibad=BAD(VIEW(icx),  
 c$$$     $              nvk(MAXS(icx)+igood),  
 c$$$     $              nst(MAXS(icx)+igood))  
 c$$$            endif  
 c$$$            badcl=badcl*ibad  
 c$$$         enddo  
 c$$$         if(badcl.eq.0)then     !<<<<<<<<<<<<<< BAD cut  
 c$$$            cl_single(icx)=0    !<<<<<<<<<<<<<< BAD cut  
 c$$$            goto 10             !<<<<<<<<<<<<<< BAD cut  
 c$$$         endif                  !<<<<<<<<<<<<<< BAD cut  
 c$$$*     ----------------------------------------------------  
 c$$$          
 c$$$         cl_good(icx)=1  
 c$$$         nplx=npl(VIEW(icx))  
 c$$$         nldx=nld(MAXS(icx),VIEW(icx))  
 c$$$          
 c$$$         do icy=1,nclstr1       !loop on cluster (Y)  
 c$$$            if(mod(VIEW(icy),2).eq.0)goto 20  
 c$$$              
 c$$$*     ----------------------------------------------------  
 c$$$*     cut on charge (Y VIEW)  
 c$$$            if(dedx(icy).lt.dedx_y_min)then  
 c$$$               cl_single(icy)=0  
 c$$$               goto 20  
 c$$$            endif  
 c$$$*     cut BAD (Y VIEW)              
 c$$$            badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy)))  
 c$$$            ifirst=INDSTART(icy)  
 c$$$            if(icy.ne.nclstr1) then  
 c$$$               ilast=INDSTART(icy+1)-1  
 c$$$            else  
 c$$$               ilast=TOTCLLENGTH  
 c$$$            endif  
 c$$$            badcl=badseed  
 c$$$            do igood=-ngoodstr,ngoodstr  
 c$$$               ibad=1  
 c$$$               if((INDMAX(icy)+igood).gt.ifirst.and.  
 c$$$     $              (INDMAX(icy)+igood).lt.ilast.and.  
 c$$$     $              .true.)  
 c$$$     $              ibad=BAD(VIEW(icy),  
 c$$$     $              nvk(MAXS(icy)+igood),  
 c$$$     $              nst(MAXS(icy)+igood))  
 c$$$               badcl=badcl*ibad  
 c$$$            enddo  
 c$$$            if(badcl.eq.0)then  !<<<<<<<<<<<<<< BAD cut  
 c$$$               cl_single(icy)=0 !<<<<<<<<<<<<<< BAD cut  
 c$$$               goto 20          !<<<<<<<<<<<<<< BAD cut  
 c$$$            endif               !<<<<<<<<<<<<<< BAD cut  
 c$$$*     ----------------------------------------------------  
 c$$$              
 c$$$              
 c$$$            cl_good(icy)=1                    
 c$$$            nply=npl(VIEW(icy))  
 c$$$            nldy=nld(MAXS(icy),VIEW(icy))  
 c$$$              
 c$$$*     ----------------------------------------------  
 c$$$*     CONDITION TO FORM A COUPLE  
 c$$$*     ----------------------------------------------  
 c$$$*     geometrical consistency (same plane and ladder)  
 c$$$            if(nply.eq.nplx.and.nldy.eq.nldx)then  
 c$$$*     charge correlation  
 c$$$*     ===========================================================  
 c$$$*     this version of the subroutine is used for the calibration  
 c$$$*     thus charge-correlation selection is obviously removed  
 c$$$*     ===========================================================  
 c$$$c$$$               ddd=(dedx(icy)  
 c$$$c$$$     $              -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))  
 c$$$c$$$               ddd=ddd/sqrt(kch(nplx,nldx)**2+1)  
 c$$$c$$$               cut=chcut*sch(nplx,nldx)  
 c$$$c$$$               if(abs(ddd).gt.cut)goto 20 !charge not consistent  
 c$$$*     ===========================================================  
 c$$$                
 c$$$                
 c$$$*     ------------------> COUPLE <------------------  
 c$$$*     check to do not overflow vector dimentions  
 c$$$c$$$               if(ncp_plane(nplx).gt.ncouplemax)then  
 c$$$c$$$                  if(DEBUG)print*,  
 c$$$c$$$     $                    ' ** warning ** number of identified'//  
 c$$$c$$$     $                    ' couples on plane ',nplx,  
 c$$$c$$$     $                    ' exceeds vector dimention'//  
 c$$$c$$$     $                    ' ( ',ncouplemax,' )'  
 c$$$c$$$c     good2=.false.  
 c$$$c$$$c     goto 880   !fill ntp and go to next event  
 c$$$c$$$                  iflag=1  
 c$$$c$$$                  return  
 c$$$c$$$               endif  
 c$$$                
 c$$$               if(ncp_plane(nplx).eq.ncouplemax)then  
 c$$$                  if(verbose)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  
 c$$$                
 c$$$               ncp_plane(nplx) = ncp_plane(nplx) + 1  
 c$$$               clx(nplx,ncp_plane(nplx))=icx  
 c$$$               cly(nply,ncp_plane(nplx))=icy  
 c$$$               cl_single(icx)=0  
 c$$$               cl_single(icy)=0  
 c$$$            endif                                
 c$$$*     ----------------------------------------------  
 c$$$  
 c$$$ 20         continue  
 c$$$         enddo                  !end loop on clusters(Y)  
 c$$$          
 c$$$ 10      continue  
 c$$$      enddo                     !end loop on clusters(X)  
 c$$$        
 c$$$        
 c$$$      do icl=1,nclstr1  
 c$$$         if(cl_single(icl).eq.1)then  
 c$$$            ip=npl(VIEW(icl))  
 c$$$            ncls(ip)=ncls(ip)+1  
 c$$$            cls(ip,ncls(ip))=icl  
 c$$$         endif  
 c$$$      enddo  
 c$$$        
 c$$$        
 c$$$      if(DEBUG)then  
 c$$$         print*,'clusters  ',nclstr1  
 c$$$         print*,'good    ',(cl_good(i),i=1,nclstr1)  
 c$$$         print*,'singles ',(cl_single(i),i=1,nclstr1)  
 c$$$         print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)  
 c$$$      endif  
 c$$$        
 c$$$      do ip=1,6  
 c$$$         ncp_tot=ncp_tot+ncp_plane(ip)  
 c$$$      enddo  
 c$$$c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)  
 c$$$        
 c$$$      if(ncp_tot.gt.ncp_max)then  
 c$$$         if(verbose)print*,  
 c$$$     $           '** warning ** number of identified '//  
 c$$$     $           'couples exceeds upper limit for Hough tr. '  
 c$$$     $           ,'( ',ncp_max,' )'              
 c$$$c            good2=.false.  
 c$$$c     goto 880       !fill ntp and go to next event  
 c$$$         iflag=1  
 c$$$         return  
 c$$$      endif  
 c$$$        
 c$$$      return  
 c$$$      end  
 c$$$  
         
 ***************************************************  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 **************************************************  
1864    
1865        subroutine cp_to_doubtrip(iflag)        subroutine cp_to_doubtrip(iflag)
 c*****************************************************  
 c     02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
1866    
1867        include 'commontracker.f'        include 'commontracker.f'
1868        include 'level1.f'        include 'level1.f'
1869        include 'common_momanhough.f'        include 'common_momanhough.f'
 c      include 'momanhough_init.f'  
1870        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
1871        include 'common_mini_2.f'        include 'common_mini_2.f'
1872        include 'calib.f'        include 'calib.f'
 c      include 'level1.f'  
1873    
1874    
1875  *     output flag  *     output flag
# Line 2245  c      double precision xm3,ym3,zm3 Line 1903  c      double precision xm3,ym3,zm3
1903  *     -----------------------------  *     -----------------------------
1904    
1905    
1906    *     --------------------------------------------
1907    *     put a limit to the maximum number of couples
1908    *     per plane, in order to apply hough transform
1909    *     (couples recovered during track refinement)
1910    *     --------------------------------------------
1911          do ip=1,nplanes
1912             if(ncp_plane(ip).gt.ncouplelimit)then
1913                mask_view(nviewx(ip)) = 8
1914                mask_view(nviewy(ip)) = 8
1915             endif
1916          enddo
1917    
1918    
1919        ndblt=0                   !number of doublets        ndblt=0                   !number of doublets
1920        ntrpt=0                   !number of triplets        ntrpt=0                   !number of triplets
1921                
1922        do ip1=1,(nplanes-1)      !loop on planes  - COPPIA 1        do ip1=1,(nplanes-1)      !loop on planes  - COPPIA 1
1923           do is1=1,2             !loop on sensors - COPPIA 1           if(  mask_view(nviewx(ip1)).ne.0 .or.
1924                     $        mask_view(nviewy(ip1)).ne.0 )goto 10 !skip plane
1925             do is1=1,2             !loop on sensors - COPPIA 1            
1926              do icp1=1,ncp_plane(ip1) !loop on COPPIA 1              do icp1=1,ncp_plane(ip1) !loop on COPPIA 1
1927                 icx1=clx(ip1,icp1)                 icx1=clx(ip1,icp1)
1928                 icy1=cly(ip1,icp1)                 icy1=cly(ip1,icp1)
1929  c               call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)  c               call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)
1930                 call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)  c               call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)
1931                   call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.,0.,0.)
1932                 xm1=xPAM                 xm1=xPAM
1933                 ym1=yPAM                 ym1=yPAM
1934                 zm1=zPAM                                   zm1=zPAM                  
1935  c     print*,'***',is1,xm1,ym1,zm1  c     print*,'***',is1,xm1,ym1,zm1
1936    
1937                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2
1938                      if(  mask_view(nviewx(ip2)).ne.0 .or.
1939         $                 mask_view(nviewy(ip2)).ne.0 )goto 20 !skip plane
1940                    do is2=1,2    !loop on sensors -ndblt COPPIA 2                    do is2=1,2    !loop on sensors -ndblt COPPIA 2
1941                                            
1942                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2
# Line 2269  c     print*,'***',is1,xm1,ym1,zm1 Line 1944  c     print*,'***',is1,xm1,ym1,zm1
1944                          icy2=cly(ip2,icp2)                          icy2=cly(ip2,icp2)
1945  c                        call xyz_PAM  c                        call xyz_PAM
1946  c     $                       (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1)  c     $                       (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1)
1947    c                        call xyz_PAM
1948    c     $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1)
1949                          call xyz_PAM                          call xyz_PAM
1950       $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1)       $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.,0.,0.)
1951                          xm2=xPAM                          xm2=xPAM
1952                          ym2=yPAM                          ym2=yPAM
1953                          zm2=zPAM                          zm2=zPAM
# Line 2322  c$$$ Line 1999  c$$$
1999  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2000    
2001    
2002                          if(ip2.eq.nplanes)goto 30 !no possible combination with 3 couples                          if(ip2.eq.nplanes)goto 31 !no possible combination with 3 couples
2003    
2004                          do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3                          do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3
2005                               if(  mask_view(nviewx(ip3)).ne.0 .or.
2006         $                          mask_view(nviewy(ip3)).ne.0 )goto 30 !skip plane
2007                             do is3=1,2 !loop on sensors - COPPIA 3                             do is3=1,2 !loop on sensors - COPPIA 3
2008                                                                
2009                                do icp3=1,ncp_plane(ip3) !loop on COPPIA 3                                do icp3=1,ncp_plane(ip3) !loop on COPPIA 3
# Line 2331  c$$$ Line 2011  c$$$
2011                                   icy3=cly(ip3,icp3)                                   icy3=cly(ip3,icp3)
2012  c                                 call xyz_PAM  c                                 call xyz_PAM
2013  c     $                               (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1)  c     $                               (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1)
2014    c                                 call xyz_PAM
2015    c     $                               (icx3,icy3,is3,PFAdef,PFAdef,0.,0.) !(1)
2016                                   call xyz_PAM                                   call xyz_PAM
2017       $                               (icx3,icy3,is3,PFAdef,PFAdef,0.,0.) !(1)       $                                (icx3,icy3,is3,PFAdef,PFAdef
2018         $                                ,0.,0.,0.,0.)
2019                                   xm3=xPAM                                   xm3=xPAM
2020                                   ym3=yPAM                                   ym3=yPAM
2021                                   zm3=zPAM                                   zm3=zPAM
# Line 2400  c     print*,alfaxz1(ntrpt),alfaxz2(ntrp Line 2083  c     print*,alfaxz1(ntrpt),alfaxz2(ntrp
2083                                endif                                endif
2084                             enddo !end loop on COPPIA 3                             enddo !end loop on COPPIA 3
2085                          enddo   !end loop on sensors - COPPIA 3                          enddo   !end loop on sensors - COPPIA 3
2086     30                     continue
2087                       enddo      !end loop on planes  - COPPIA 3                       enddo      !end loop on planes  - COPPIA 3
2088   30                  continue   31                  continue
2089                                            
2090   1                enddo         !end loop on COPPIA 2   1                enddo         !end loop on COPPIA 2
2091                 enddo            !end loop on sensors - COPPIA 2                 enddo            !end loop on sensors - COPPIA 2
2092     20            continue
2093              enddo               !end loop on planes  - COPPIA 2              enddo               !end loop on planes  - COPPIA 2
2094                            
2095           enddo                  !end loop on COPPIA1           enddo                  !end loop on COPPIA1
2096        enddo                     !end loop on sensors - COPPIA 1        enddo                     !end loop on sensors - COPPIA 1
2097     10   continue
2098        enddo                     !end loop on planes  - COPPIA 1        enddo                     !end loop on planes  - COPPIA 1
2099                
2100        if(DEBUG)then        if(DEBUG)then
# Line 2801  c     print*,'check cp_used' Line 2487  c     print*,'check cp_used'
2487           enddo           enddo
2488  c         if(ncpused.lt.ncpxz_min)goto 22288 !next triplet  c         if(ncpused.lt.ncpxz_min)goto 22288 !next triplet
2489           if(npt.lt.nptxz_min)goto 22288     !next triplet           if(npt.lt.nptxz_min)goto 22288     !next triplet
2490           if(nplused.lt.nplxz_min)goto 22288 !next doublet           if(nplused.lt.nplxz_min)goto 22288 !next triplet
2491                    
2492  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
2493  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
# Line 2878  c$$$     $           ,(tr_cloud(iii),iii Line 2564  c$$$     $           ,(tr_cloud(iii),iii
2564  **************************************************  **************************************************
2565    
2566        subroutine clouds_to_ctrack(iflag)        subroutine clouds_to_ctrack(iflag)
 c*****************************************************  
 c     02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
2567    
2568        include 'commontracker.f'        include 'commontracker.f'
2569        include 'level1.f'        include 'level1.f'
# Line 2888  c*************************************** Line 2571  c***************************************
2571        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
2572        include 'common_mini_2.f'        include 'common_mini_2.f'
2573        include 'common_mech.f'        include 'common_mech.f'
2574  c      include 'momanhough_init.f'  
2575    
2576    
2577  *     output flag  *     output flag
# Line 2912  c      include 'momanhough_init.f' Line 2595  c      include 'momanhough_init.f'
2595  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2596  *     variables for track fitting  *     variables for track fitting
2597        double precision AL_INI(5)        double precision AL_INI(5)
2598        double precision tath  c      double precision tath
2599  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2600  c      real fitz(nplanes)        !z coordinates of the planes in cm  c      real fitz(nplanes)        !z coordinates of the planes in cm
2601    
# Line 2978  c      real fitz(nplanes)        !z coor Line 2661  c      real fitz(nplanes)        !z coor
2661                 nplused=nplused+ hit_plane(ip)                 nplused=nplused+ hit_plane(ip)
2662              enddo              enddo
2663                            
2664              if(nplused.lt.nplxz_min)goto 888 !next doublet  c            if(nplused.lt.nplxz_min)goto 888 !next doublet
2665                if(nplused.lt.nplyz_min)goto 888 !next doublet
2666              if(ncp_ok.lt.ncpok)goto 888 !next cloud              if(ncp_ok.lt.ncpok)goto 888 !next cloud
2667                            
2668              if(DEBUG)then              if(DEBUG)then
# Line 3005  c$$$  print*,'6 -- ',(cly(6,i),i=1,ncp_p Line 2689  c$$$  print*,'6 -- ',(cly(6,i),i=1,ncp_p
2689  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'
2690                            
2691  *     -------> INITIAL GUESS <-------  *     -------> INITIAL GUESS <-------
2692              AL_INI(1) = dreal(alfaxz1_av(ixz))  cccc       SBAGLIATO
2693              AL_INI(2) = dreal(alfayz1_av(iyz))  c$$$            AL_INI(1) = dreal(alfaxz1_av(ixz))
2694              AL_INI(4) = PIGR + datan(dreal(alfayz2_av(iyz))  c$$$            AL_INI(2) = dreal(alfayz1_av(iyz))
2695       $           /dreal(alfaxz2_av(ixz)))  c$$$            AL_INI(4) = PIGR + datan(dreal(alfayz2_av(iyz))
2696              tath      = -dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))  c$$$     $           /dreal(alfaxz2_av(ixz)))
2697              AL_INI(3) = tath/sqrt(1+tath**2)  c$$$            tath      = -dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))
2698              AL_INI(5) = (1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.  c$$$            AL_INI(3) = tath/sqrt(1+tath**2)
2699                c$$$            AL_INI(5) = (1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.
2700  c     print*,'*******',AL_INI(5)  cccc       GIUSTO (ma si sua guess())
2701              if(AL_INI(5).gt.defmax)goto 888 !next cloud  c$$$            AL_INI(1) = dreal(alfaxz1_av(ixz))
2702                c$$$            AL_INI(2) = dreal(alfayz1_av(iyz))
2703  c     print*,'alfaxz2, alfayz2 '  c$$$            tath      = -dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))
2704  c     $              ,alfaxz2_av(ixz),alfayz2_av(iyz)  c$$$            AL_INI(3) = tath/sqrt(1+tath**2)
2705                c$$$            IF(alfaxz2_av(ixz).NE.0)THEN
2706  *     -------> INITIAL GUESS <-------  c$$$            AL_INI(4) = PIGR + datan(dreal(alfayz2_av(iyz))
2707  c     print*,'AL_INI ',(al_ini(i),i=1,5)  c$$$     $           /dreal(alfaxz2_av(ixz)))
2708                c$$$            ELSE
2709    c$$$               AL_INI(4) = acos(-1.)/2
2710    c$$$               IF(alfayz2_av(iyz).LT.0)AL_INI(4) = AL_INI(4)+acos(-1.)
2711    c$$$            ENDIF
2712    c$$$            IF(alfaxz2_av(ixz).LT.0)AL_INI(4)= acos(-1.)+ AL_INI(4)
2713    c$$$            AL_INI(4) = -acos(-1.) + AL_INI(4) !from incidence direction to tracking rs
2714    c$$$            
2715    c$$$            AL_INI(5) = (1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.
2716    c$$$            
2717    c$$$            if(AL_INI(5).gt.defmax)goto 888 !next cloud
2718                            
2719              if(DEBUG)then              if(DEBUG)then
2720                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))
2721                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))
# Line 3072  c     print*,'AL_INI ',(al_ini(i),i=1,5) Line 2766  c     print*,'AL_INI ',(al_ini(i),i=1,5)
2766  *                                   *************************  *                                   *************************
2767  c                                    call xyz_PAM(icx,icy,is,  c                                    call xyz_PAM(icx,icy,is,
2768  c     $                                   'COG2','COG2',0.,0.)  c     $                                   'COG2','COG2',0.,0.)
2769    c                                    call xyz_PAM(icx,icy,is, !(1)
2770    c     $                                   PFAdef,PFAdef,0.,0.) !(1)
2771                                      call xyz_PAM(icx,icy,is, !(1)                                      call xyz_PAM(icx,icy,is, !(1)
2772       $                                   PFAdef,PFAdef,0.,0.) !(1)       $                                   PFAdef,PFAdef,0.,0.,0.,0.)
2773  *                                   *************************  *                                   *************************
2774  *                                   -----------------------------  *                                   -----------------------------
2775                                      xgood(nplanes-ip+1)=1.                                      xgood(nplanes-ip+1)=1.
# Line 3089  c     $                                 Line 2785  c     $                                
2785  *     **********************************************************  *     **********************************************************
2786  *     ************************** FIT *** FIT *** FIT *** FIT ***  *     ************************** FIT *** FIT *** FIT *** FIT ***
2787  *     **********************************************************  *     **********************************************************
2788    cccc  scommentare se si usa al_ini della nuvola
2789    c$$$                              do i=1,5
2790    c$$$                                 AL(i)=AL_INI(i)
2791    c$$$                              enddo
2792                                  call guess()
2793                                do i=1,5                                do i=1,5
2794                                   AL(i)=AL_INI(i)                                   AL_INI(i)=AL(i)
2795                                enddo                                enddo
2796                                ifail=0 !error flag in chi^2 computation                                ifail=0 !error flag in chi^2 computation
2797                                jstep=0 !number of  minimization steps                                jstep=0 !number of  minimization steps
2798                                iprint=0                                iprint=0
2799                                if(DEBUG)iprint=1  c                              if(DEBUG)iprint=1
2800                                  if(DEBUG)iprint=2
2801                                call mini2(jstep,ifail,iprint)                                call mini2(jstep,ifail,iprint)
2802                                if(ifail.ne.0) then                                if(ifail.ne.0) then
2803                                   if(DEBUG)then                                   if(DEBUG)then
2804                                      print *,                                      print *,
2805       $                              '*** MINIMIZATION FAILURE *** '       $                              '*** MINIMIZATION FAILURE *** '
2806       $                              //'(mini2 in clouds_to_ctrack)'       $                              //'(clouds_to_ctrack)'
2807                                        print*,'initial guess: '
2808    
2809                                        print*,'AL_INI(1) = ',AL_INI(1)
2810                                        print*,'AL_INI(2) = ',AL_INI(2)
2811                                        print*,'AL_INI(3) = ',AL_INI(3)
2812                                        print*,'AL_INI(4) = ',AL_INI(4)
2813                                        print*,'AL_INI(5) = ',AL_INI(5)
2814                                   endif                                   endif
2815                                   chi2=-chi2  c                                 chi2=-chi2
2816                                endif                                endif
2817  *     **********************************************************  *     **********************************************************
2818  *     ************************** FIT *** FIT *** FIT *** FIT ***  *     ************************** FIT *** FIT *** FIT *** FIT ***
# Line 3131  c                                 goto 8 Line 2840  c                                 goto 8
2840                                                                
2841                                ntracks = ntracks + 1                                ntracks = ntracks + 1
2842                                                                
 c$$$                              ndof=0                                  
2843                                do ip=1,nplanes                                do ip=1,nplanes
2844  c$$$                                 ndof=ndof  
 c$$$     $                                +int(xgood(ip))  
 c$$$     $                                +int(ygood(ip))  
2845                                   XV_STORE(ip,ntracks)=sngl(xv(ip))                                   XV_STORE(ip,ntracks)=sngl(xv(ip))
2846                                   YV_STORE(ip,ntracks)=sngl(yv(ip))                                   YV_STORE(ip,ntracks)=sngl(yv(ip))
2847                                   ZV_STORE(ip,ntracks)=sngl(zv(ip))                                                                       ZV_STORE(ip,ntracks)=sngl(zv(ip))                                    
# Line 3154  c$$$     $                               Line 2860  c$$$     $                              
2860                                   if(hit_plane(ip).ne.0)then                                   if(hit_plane(ip).ne.0)then
2861                                      CP_STORE(nplanes-ip+1,ntracks)=                                      CP_STORE(nplanes-ip+1,ntracks)=
2862       $                                   cp_match(ip,hit_plane(ip))       $                                   cp_match(ip,hit_plane(ip))
2863                                        SENSOR_STORE(nplanes-ip+1,ntracks)
2864         $                              = is_cp(cp_match(ip,hit_plane(ip)))
2865                                        LADDER_STORE(nplanes-ip+1,ntracks)
2866         $                                   = LADDER(
2867         $                                   clx(ip,icp_cp(
2868         $                                   cp_match(ip,hit_plane(ip)
2869         $                                   ))));
2870                                   else                                   else
2871                                      CP_STORE(nplanes-ip+1,ntracks)=0                                      CP_STORE(nplanes-ip+1,ntracks)=0
2872                                        SENSOR_STORE(nplanes-ip+1,ntracks)=0
2873                                        LADDER_STORE(nplanes-ip+1,ntracks)=0
2874                                   endif                                   endif
2875                                     BX_STORE(nplanes-ip+1,ntracks)=0!I dont need it now
2876                                     BY_STORE(nplanes-ip+1,ntracks)=0!I dont need it now
2877                                   CLS_STORE(nplanes-ip+1,ntracks)=0                                   CLS_STORE(nplanes-ip+1,ntracks)=0
2878                                   do i=1,5                                   do i=1,5
2879                                      AL_STORE(i,ntracks)=sngl(AL(i))                                      AL_STORE(i,ntracks)=sngl(AL(i))
2880                                   enddo                                   enddo
2881                                enddo                                enddo
2882                                                                
 c$$$  *                             Number of Degree Of Freedom  
 c$$$  ndof=ndof-5                            
 c$$$  *                             reduced chi^2  
 c$$$  rchi2=chi2/dble(ndof)  
2883                                RCHI2_STORE(ntracks)=chi2                                RCHI2_STORE(ntracks)=chi2
2884                                                                
2885  *     --------------------------------  *     --------------------------------
# Line 3190  c$$$  rchi2=chi2/dble(ndof) Line 2903  c$$$  rchi2=chi2/dble(ndof)
2903           return           return
2904        endif        endif
2905                
2906    c$$$      if(DEBUG)then
2907    c$$$         print*,'****** TRACK CANDIDATES ***********'
2908    c$$$         print*,'#         R. chi2        RIG'
2909    c$$$         do i=1,ntracks
2910    c$$$            print*,i,' --- ',rchi2_store(i),' --- '
2911    c$$$     $           ,1./abs(AL_STORE(5,i))
2912    c$$$         enddo
2913    c$$$         print*,'***********************************'
2914    c$$$      endif
2915        if(DEBUG)then        if(DEBUG)then
2916           print*,'****** TRACK CANDIDATES ***********'          print*,'****** TRACK CANDIDATES *****************'
2917           print*,'#         R. chi2        RIG'          print*,'#         R. chi2        RIG         ndof'
2918           do i=1,ntracks          do i=1,ntracks
2919              print*,i,' --- ',rchi2_store(i),' --- '            ndof=0                !(1)
2920       $           ,1./abs(AL_STORE(5,i))            do ii=1,nplanes       !(1)
2921           enddo              ndof=ndof           !(1)
2922           print*,'***********************************'       $           +int(xgood_store(ii,i)) !(1)
2923         $           +int(ygood_store(ii,i)) !(1)
2924              enddo                 !(1)
2925              print*,i,' --- ',rchi2_store(i),' --- '
2926         $         ,1./abs(AL_STORE(5,i)),' --- ',ndof
2927            enddo
2928            print*,'*****************************************'
2929        endif        endif
2930                
2931                
# Line 3216  c$$$  rchi2=chi2/dble(ndof) Line 2944  c$$$  rchi2=chi2/dble(ndof)
2944    
2945        subroutine refine_track(ibest)        subroutine refine_track(ibest)
2946    
 c******************************************************  
 cccccc 06/10/2005 modified by elena vannuccini ---> (1)  
 cccccc 31/01/2006 modified by elena vannuccini ---> (2)  
 cccccc 12/08/2006 modified by elena vannucicni ---> (3)  
 c******************************************************  
2947    
2948        include 'commontracker.f'        include 'commontracker.f'
2949        include 'level1.f'        include 'level1.f'
# Line 3228  c*************************************** Line 2951  c***************************************
2951        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
2952        include 'common_mini_2.f'        include 'common_mini_2.f'
2953        include 'common_mech.f'        include 'common_mech.f'
 c      include 'momanhough_init.f'  
 c      include 'level1.f'  
2954        include 'calib.f'        include 'calib.f'
2955    
   
2956  *     flag to chose PFA  *     flag to chose PFA
2957        character*10 PFA        character*10 PFA
2958        common/FINALPFA/PFA        common/FINALPFA/PFA
2959    
2960          real k(6)
2961          DATA k/1.099730,0.418900,0.220939,0.220907,0.418771,1.100674/
2962    
2963          real xp,yp,zp
2964          real xyzp(3),bxyz(3)
2965          equivalence (xp,xyzp(1)),(yp,xyzp(2)),(zp,xyzp(3))
2966    
2967  *     =================================================  *     =================================================
2968  *     new estimate of positions using ETA algorithm  *     new estimate of positions using ETA algorithm
2969  *                          and  *                          and
# Line 3245  c      include 'level1.f' Line 2972  c      include 'level1.f'
2972        call track_init        call track_init
2973        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
2974    
2975             xP=XV_STORE(nplanes-ip+1,ibest)
2976             yP=YV_STORE(nplanes-ip+1,ibest)
2977             zP=ZV_STORE(nplanes-ip+1,ibest)
2978             call gufld(xyzp,bxyz)
2979             BX_STORE(nplanes-ip+1,ibest)=bxyz(1)
2980             BY_STORE(nplanes-ip+1,ibest)=bxyz(2)
2981    c$$$  bxyz(1)=0
2982    c$$$         bxyz(2)=0
2983    c$$$         bxyz(3)=0
2984    *     |||||||||||||||||||||||||||||||||||||||||||||||||
2985  *     -------------------------------------------------  *     -------------------------------------------------
2986  *     If the plane has been already included, it just  *     If the plane has been already included, it just
2987  *     computes again the coordinates of the x-y couple  *     computes again the coordinates of the x-y couple
2988  *     using improved PFAs  *     using improved PFAs
2989  *     -------------------------------------------------  *     -------------------------------------------------
2990    *     |||||||||||||||||||||||||||||||||||||||||||||||||
2991           if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and.           if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and.
2992       $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then       $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then
2993                            
# Line 3263  c      include 'level1.f' Line 3001  c      include 'level1.f'
3001       $           ,ip_cp(id),ip       $           ,ip_cp(id),ip
3002              icx=clx(ip,icp)              icx=clx(ip,icp)
3003              icy=cly(ip,icp)              icy=cly(ip,icp)
3004    c            call xyz_PAM(icx,icy,is,
3005    c     $           PFA,PFA,
3006    c     $           AXV_STORE(nplanes-ip+1,ibest),
3007    c     $           AYV_STORE(nplanes-ip+1,ibest))
3008              call xyz_PAM(icx,icy,is,              call xyz_PAM(icx,icy,is,
 c     $           'ETA2','ETA2',  
3009       $           PFA,PFA,       $           PFA,PFA,
3010       $           AXV_STORE(nplanes-ip+1,ibest),       $           AXV_STORE(nplanes-ip+1,ibest),
3011       $           AYV_STORE(nplanes-ip+1,ibest))       $           AYV_STORE(nplanes-ip+1,ibest),
3012  c$$$  call xyz_PAM(icx,icy,is,       $           bxyz(1),
3013  c$$$  $              'COG2','COG2',       $           bxyz(2)
3014  c$$$  $              0.,       $           )
3015  c$$$  $              0.)  
3016              xm(nplanes-ip+1) = xPAM              xm(nplanes-ip+1) = xPAM
3017              ym(nplanes-ip+1) = yPAM              ym(nplanes-ip+1) = yPAM
3018              zm(nplanes-ip+1) = zPAM              zm(nplanes-ip+1) = zPAM
# Line 3280  c$$$  $              0.) Line 3021  c$$$  $              0.)
3021              resx(nplanes-ip+1) = resxPAM              resx(nplanes-ip+1) = resxPAM
3022              resy(nplanes-ip+1) = resyPAM              resy(nplanes-ip+1) = resyPAM
3023    
3024  c            dedxtrk(nplanes-ip+1) = (dedx(icx)+dedx(icy))/2. !(1)              dedxtrk_x(nplanes-ip+1)=sgnl(icx)/mip(VIEW(icx),LADDER(icx))
3025              dedxtrk_x(nplanes-ip+1)=dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)              dedxtrk_y(nplanes-ip+1)=sgnl(icy)/mip(VIEW(icy),LADDER(icy))
             dedxtrk_y(nplanes-ip+1)=dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)  
3026                            
3027    *     |||||||||||||||||||||||||||||||||||||||||||||||||
3028  *     -------------------------------------------------  *     -------------------------------------------------
3029  *     If the plane has NOT  been already included,  *     If the plane has NOT  been already included,
3030  *     it tries to include a COUPLE or a single cluster  *     it tries to include a COUPLE or a single cluster
3031  *     -------------------------------------------------  *     -------------------------------------------------
3032    *     |||||||||||||||||||||||||||||||||||||||||||||||||
3033           else                             else                  
3034                                
3035              xgood(nplanes-ip+1)=0              xgood(nplanes-ip+1)=0
# Line 3295  c            dedxtrk(nplanes-ip+1) = (de Line 3037  c            dedxtrk(nplanes-ip+1) = (de
3037                                
3038  *     --------------------------------------------------------------  *     --------------------------------------------------------------
3039  *     determine which ladder and sensor are intersected by the track  *     determine which ladder and sensor are intersected by the track
             xP=XV_STORE(nplanes-ip+1,ibest)  
             yP=YV_STORE(nplanes-ip+1,ibest)  
             zP=ZV_STORE(nplanes-ip+1,ibest)  
3040              call whichsensor(ip,xP,yP,nldt,ist)              call whichsensor(ip,xP,yP,nldt,ist)
3041  *     if the track hit the plane in a dead area, go to the next plane  *     if the track hit the plane in a dead area, go to the next plane
3042              if(nldt.eq.0.or.ist.eq.0)goto 133              if(nldt.eq.0.or.ist.eq.0)goto 133
3043    
3044                SENSOR_STORE(nplanes-ip+1,IBEST)=ist
3045                LADDER_STORE(nplanes-ip+1,IBEST)=nldt
3046  *     --------------------------------------------------------------  *     --------------------------------------------------------------
3047    
3048              if(DEBUG)then              if(DEBUG)then
# Line 3337  c     $              cl_used(icy).eq.1.o Line 3079  c     $              cl_used(icy).eq.1.o
3079  *            *          
3080                 call xyz_PAM(icx,icy,ist,                 call xyz_PAM(icx,icy,ist,
3081       $              PFA,PFA,       $              PFA,PFA,
 c     $              'ETA2','ETA2',  
3082       $              AXV_STORE(nplanes-ip+1,ibest),       $              AXV_STORE(nplanes-ip+1,ibest),
3083       $              AYV_STORE(nplanes-ip+1,ibest))       $              AYV_STORE(nplanes-ip+1,ibest),
3084         $              bxyz(1),
3085         $              bxyz(2)
3086         $              )
3087                                
3088                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3089    c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3090                 id=id_cp(ip,icp,ist)                 id=id_cp(ip,icp,ist)
3091                 if(DEBUG)print*,'( couple ',id                 if(DEBUG)print*,'( couple ',id
3092       $              ,' ) normalized distance ',distance       $              ,' ) distance ',distance
3093                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3094                    xmm = xPAM                    xmm = xPAM
3095                    ymm = yPAM                    ymm = yPAM
# Line 3353  c     $              'ETA2','ETA2', Line 3098  c     $              'ETA2','ETA2',
3098                    rymm = resyPAM                    rymm = resyPAM
3099                    distmin = distance                    distmin = distance
3100                    idm = id                                      idm = id                  
3101  c                 dedxmm = (dedx(icx)+dedx(icy))/2. !(1)                    dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
3102                    dedxmmx = dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)                    dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
3103                    dedxmmy = dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)  c     QUIQUI --> non devo moltiplicare per clinc?!?!?!
3104                      clincnewc=10*sqrt(rymm**2+rxmm**2 !QUIQUI
3105         $                 +RCHI2_STORE(ibest)*k(ip)*(cov(1,1)+cov(2,2))) !QUIQUI
3106                 endif                 endif
3107   1188          continue   1188          continue
3108              enddo               !end loop on couples on plane icp              enddo               !end loop on couples on plane icp
3109              if(distmin.le.clinc)then                    c            if(distmin.le.clinc)then     !QUIQUI              
3110                if(distmin.le.clincnewc)then     !QUIQUI              
3111  *              -----------------------------------  *              -----------------------------------
3112                 xm(nplanes-ip+1) = xmm         !<<<                 xm(nplanes-ip+1) = xmm !<<<
3113                 ym(nplanes-ip+1) = ymm         !<<<                 ym(nplanes-ip+1) = ymm !<<<
3114                 zm(nplanes-ip+1) = zmm         !<<<                 zm(nplanes-ip+1) = zmm !<<<
3115                 xgood(nplanes-ip+1) = 1        !<<<                 xgood(nplanes-ip+1) = 1 !<<<
3116                 ygood(nplanes-ip+1) = 1        !<<<                 ygood(nplanes-ip+1) = 1 !<<<
3117                 resx(nplanes-ip+1)=rxmm        !<<<                 resx(nplanes-ip+1)=rxmm !<<<
3118                 resy(nplanes-ip+1)=rymm        !<<<                 resy(nplanes-ip+1)=rymm !<<<
3119  c              dedxtrk(nplanes-ip+1) = dedxmm !<<<  !(1)                 dedxtrk_x(nplanes-ip+1) = dedxmmx !<<<
3120                 dedxtrk_x(nplanes-ip+1) = dedxmmx    !(1)                 dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<
                dedxtrk_y(nplanes-ip+1) = dedxmmy    !(1)  
3121  *              -----------------------------------  *              -----------------------------------
3122                 CP_STORE(nplanes-ip+1,ibest)=idm                       CP_STORE(nplanes-ip+1,ibest)=idm      
3123                 if(DEBUG)print*,'%%%% included couple ',idm                 if(DEBUG)print*,'%%%% included couple ',idm
3124       $              ,' (norm.dist.= ',distmin,', cut ',clinc,' )'       $              ,' (dist.= ',distmin,', cut ',clinc,' )'
3125                 goto 133         !next plane                 goto 133         !next plane
3126              endif              endif
3127  *     ================================================  *     ================================================
# Line 3407  c            if(DEBUG)print*,'>>>> try t Line 3154  c            if(DEBUG)print*,'>>>> try t
3154  c               if(cl_used(icx).eq.1)goto 11881 !if the X cluster is already used  c               if(cl_used(icx).eq.1)goto 11881 !if the X cluster is already used
3155                 if(cl_used(icx).ne.0)goto 11881 !if the X cluster is already used  !(3)                 if(cl_used(icx).ne.0)goto 11881 !if the X cluster is already used  !(3)
3156  *                                              !jump to the Y cluster  *                                              !jump to the Y cluster
3157    c               call xyz_PAM(icx,0,ist,
3158    c     $              PFA,PFA,
3159    c     $              AXV_STORE(nplanes-ip+1,ibest),0.)              
3160                 call xyz_PAM(icx,0,ist,                 call xyz_PAM(icx,0,ist,
 c     $              'ETA2','ETA2',  
3161       $              PFA,PFA,       $              PFA,PFA,
3162       $              AXV_STORE(nplanes-ip+1,ibest),0.)                     $              AXV_STORE(nplanes-ip+1,ibest),0.,
3163         $              bxyz(1),
3164         $              bxyz(2)
3165         $              )              
3166                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3167  c     if(DEBUG)print*,'normalized distance ',distance  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3168                 if(DEBUG)print*,'( cl-X ',icx                 if(DEBUG)print*,'( cl-X ',icx
3169       $              ,' in cp ',id,' ) normalized distance ',distance       $              ,' in cp ',id,' ) distance ',distance
3170                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3171                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3172                    ymm_A = yPAM_A                    ymm_A = yPAM_A
# Line 3426  c     if(DEBUG)print*,'normalized distan Line 3178  c     if(DEBUG)print*,'normalized distan
3178                    rymm = resyPAM                    rymm = resyPAM
3179                    distmin = distance                    distmin = distance
3180                    iclm = icx                    iclm = icx
3181  c                  dedxmm = dedx(icx) !(1)  c                  dedxmm = sgnl(icx) !(1)
3182                    dedxmmx = dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)                    dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
3183                    dedxmmy = 0.        !(1)                    dedxmmy = 0.        !(1)
3184                 endif                                   endif                  
3185  11881          continue  11881          continue
# Line 3435  c                  dedxmm = dedx(icx) !( Line 3187  c                  dedxmm = dedx(icx) !(
3187  c               if(cl_used(icy).eq.1)goto 11882 !if the Y cluster is already used  c               if(cl_used(icy).eq.1)goto 11882 !if the Y cluster is already used
3188                 if(cl_used(icy).ne.0)goto 11882 !if the Y cluster is already used !(3)                 if(cl_used(icy).ne.0)goto 11882 !if the Y cluster is already used !(3)
3189  *                                              !jump to the next couple  *                                              !jump to the next couple
3190    c               call xyz_PAM(0,icy,ist,
3191    c     $              PFA,PFA,
3192    c     $              0.,AYV_STORE(nplanes-ip+1,ibest))
3193                 call xyz_PAM(0,icy,ist,                 call xyz_PAM(0,icy,ist,
 c     $              'ETA2','ETA2',  
3194       $              PFA,PFA,       $              PFA,PFA,
3195       $              0.,AYV_STORE(nplanes-ip+1,ibest))       $              0.,AYV_STORE(nplanes-ip+1,ibest),
3196         $              bxyz(1),
3197         $              bxyz(2)
3198         $              )
3199                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3200    c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3201                 if(DEBUG)print*,'( cl-Y ',icy                 if(DEBUG)print*,'( cl-Y ',icy
3202       $              ,' in cp ',id,' ) normalized distance ',distance       $              ,' in cp ',id,' ) distance ',distance
3203                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3204                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3205                    ymm_A = yPAM_A                    ymm_A = yPAM_A
# Line 3453  c     $              'ETA2','ETA2', Line 3211  c     $              'ETA2','ETA2',
3211                    rymm = resyPAM                    rymm = resyPAM
3212                    distmin = distance                    distmin = distance
3213                    iclm = icy                    iclm = icy
3214  c                 dedxmm = dedx(icy)  !(1)  c                 dedxmm = sgnl(icy)  !(1)
3215                    dedxmmx = 0.        !(1)                    dedxmmx = 0.        !(1)
3216                    dedxmmy = dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)                    dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
3217                 endif                                   endif                  
3218  11882          continue  11882          continue
3219              enddo               !end loop on cluster inside couples              enddo               !end loop on cluster inside couples
3220  *----- single clusters -----------------------------------------------      *----- single clusters -----------------------------------------------  
3221    c            print*,'## ncls(',ip,') ',ncls(ip)
3222              do ic=1,ncls(ip)    !loop on single clusters              do ic=1,ncls(ip)    !loop on single clusters
3223                 icl=cls(ip,ic)                 icl=cls(ip,ic)
3224  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
# Line 3468  c               if(cl_used(icl).eq.1.or. Line 3227  c               if(cl_used(icl).eq.1.or.
3227       $              .false.)goto 18882      !jump to the next singlet       $              .false.)goto 18882      !jump to the next singlet
3228                 if(mod(VIEW(icl),2).eq.0)then!<---- X view                 if(mod(VIEW(icl),2).eq.0)then!<---- X view
3229                    call xyz_PAM(icl,0,ist,                    call xyz_PAM(icl,0,ist,
 c     $                 'ETA2','ETA2',  
3230       $                 PFA,PFA,       $                 PFA,PFA,
3231       $                 AXV_STORE(nplanes-ip+1,ibest),0.)       $                 AXV_STORE(nplanes-ip+1,ibest),0.,
3232         $                 bxyz(1),
3233         $                 bxyz(2)
3234         $                 )
3235                 else                         !<---- Y view                 else                         !<---- Y view
3236                    call xyz_PAM(0,icl,ist,                    call xyz_PAM(0,icl,ist,
 c     $                 'ETA2','ETA2',  
3237       $                 PFA,PFA,       $                 PFA,PFA,
3238       $                 0.,AYV_STORE(nplanes-ip+1,ibest))       $                 0.,AYV_STORE(nplanes-ip+1,ibest),
3239         $                 bxyz(1),
3240         $                 bxyz(2)
3241         $                 )
3242                 endif                 endif
3243    
3244                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3245    c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3246                 if(DEBUG)print*,'( cl-s ',icl                 if(DEBUG)print*,'( cl-s ',icl
3247       $              ,' ) normalized distance ',distance       $              ,' ) distance ',distance,'<',distmin,' ?'
3248                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3249                      if(DEBUG)print*,'YES'
3250                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3251                    ymm_A = yPAM_A                    ymm_A = yPAM_A
3252                    zmm_A = zPAM_A                    zmm_A = zPAM_A
# Line 3492  c     $                 'ETA2','ETA2', Line 3257  c     $                 'ETA2','ETA2',
3257                    rymm = resyPAM                    rymm = resyPAM
3258                    distmin = distance                      distmin = distance  
3259                    iclm = icl                    iclm = icl
 c                  dedxmm = dedx(icl)                   !(1)  
3260                    if(mod(VIEW(icl),2).eq.0)then !<---- X view                    if(mod(VIEW(icl),2).eq.0)then !<---- X view
3261                       dedxmmx = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2)                       dedxmmx = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3262                       dedxmmy = 0.                       !(1)                       dedxmmy = 0.                  
3263                    else          !<---- Y view                    else          !<---- Y view
3264                       dedxmmx = 0.                       !(1)                       dedxmmx = 0.                  
3265                       dedxmmy = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2)                       dedxmmy = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3266                    endif                    endif
3267                 endif                                   endif                  
3268  18882          continue  18882          continue
3269              enddo               !end loop on single clusters              enddo               !end loop on single clusters
3270    c            print*,'## distmin ', distmin,' clinc ',clinc
3271    
3272              if(distmin.le.clinc)then                    c     QUIQUI------------
3273                  c     anche qui: non ci vuole clinc???
3274                 CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<                  if(iclm.ne.0)then
 *              ----------------------------  
3275                 if(mod(VIEW(iclm),2).eq.0)then                 if(mod(VIEW(iclm),2).eq.0)then
3276                    XGOOD(nplanes-ip+1)=1.                    clincnew=
3277                    resx(nplanes-ip+1)=rxmm       $                 20*
3278                    if(DEBUG)print*,'%%%% included X-cl ',iclm       $                 sqrt(rxmm**2+RCHI2_STORE(ibest)*k(ip)*cov(1,1))
3279       $                 ,' ( norm.dist.= ',distmin,', cut ',clinc,' )'                 else if(mod(VIEW(iclm),2).ne.0)then
3280                 else                    clincnew=
3281                    YGOOD(nplanes-ip+1)=1.       $                 10*
3282                    resy(nplanes-ip+1)=rymm       $                 sqrt(rymm**2+RCHI2_STORE(ibest)*k(ip)*cov(2,2))
3283                    if(DEBUG)print*,'%%%% included Y-cl ',iclm                 endif
3284       $                 ,' ( norm.dist.= ',distmin,', cut ',clinc,' )'  c     QUIQUI------------
3285                  
3286                   if(distmin.le.clincnew)then   !QUIQUI
3287    c     if(distmin.le.clinc)then          !QUIQUI          
3288                      
3289                      CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<    
3290    *     ----------------------------
3291    c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
3292                      if(mod(VIEW(iclm),2).eq.0)then
3293                         XGOOD(nplanes-ip+1)=1.
3294                         resx(nplanes-ip+1)=rxmm
3295                         if(DEBUG)print*,'%%%% included X-cl ',iclm
3296         $                    ,'( chi^2, ',RCHI2_STORE(ibest)
3297         $                    ,', dist.= ',distmin
3298         $                    ,', cut ',clinc,' )'
3299                      else
3300                         YGOOD(nplanes-ip+1)=1.
3301                         resy(nplanes-ip+1)=rymm
3302                         if(DEBUG)print*,'%%%% included Y-cl ',iclm
3303         $                    ,'( chi^2, ',RCHI2_STORE(ibest)
3304         $                    ,', dist.= ', distmin
3305         $                    ,', cut ',clinc,' )'
3306                      endif
3307    c     print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
3308    *     ----------------------------
3309                      xm_A(nplanes-ip+1) = xmm_A
3310                      ym_A(nplanes-ip+1) = ymm_A
3311                      xm_B(nplanes-ip+1) = xmm_B
3312                      ym_B(nplanes-ip+1) = ymm_B
3313                      zm(nplanes-ip+1) = (zmm_A+zmm_B)/2.
3314                      dedxtrk_x(nplanes-ip+1) = dedxmmx !<<<
3315                      dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<
3316    *     ----------------------------
3317                 endif                 endif
 *              ----------------------------  
                xm_A(nplanes-ip+1) = xmm_A  
                ym_A(nplanes-ip+1) = ymm_A  
                xm_B(nplanes-ip+1) = xmm_B  
                ym_B(nplanes-ip+1) = ymm_B  
                zm(nplanes-ip+1) = (zmm_A+zmm_B)/2.  
 c              dedxtrk(nplanes-ip+1) = dedxmm !<<<    !(1)  
                dedxtrk_x(nplanes-ip+1) = dedxmmx !<<< !(1)  
                dedxtrk_y(nplanes-ip+1) = dedxmmy !<<< !(1)  
 *              ----------------------------  
3318              endif              endif
3319           endif           endif
3320   133     continue   133     continue
# Line 3547  c              dedxtrk(nplanes-ip+1) = d Line 3333  c              dedxtrk(nplanes-ip+1) = d
3333  *                                                 *  *                                                 *
3334  *                                                 *  *                                                 *
3335  **************************************************  **************************************************
 cccccc 12/08/2006 modified by elena ---> (1)  
3336  *  *
3337        subroutine clean_XYclouds(ibest,iflag)        subroutine clean_XYclouds(ibest,iflag)
3338    
3339        include 'commontracker.f'        include 'commontracker.f'
3340        include 'level1.f'        include 'level1.f'
3341        include 'common_momanhough.f'        include 'common_momanhough.f'
3342  c      include 'momanhough_init.f'        include 'level2.f'      
       include 'level2.f'        !(1)  
 c      include 'calib.f'  
 c      include 'level1.f'  
   
   
3343    
3344        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
3345    
# Line 3569  c      include 'level1.f' Line 3349  c      include 'level1.f'
3349              if(id.ne.0)then              if(id.ne.0)then
3350                 iclx=clx(ip,icp_cp(id))                 iclx=clx(ip,icp_cp(id))
3351                 icly=cly(ip,icp_cp(id))                 icly=cly(ip,icp_cp(id))
3352  c               cl_used(iclx)=1  !tag used clusters                 cl_used(iclx)=ntrk  !tag used clusters
3353  c               cl_used(icly)=1  !tag used clusters                 cl_used(icly)=ntrk  !tag used clusters
                cl_used(iclx)=ntrk  !tag used clusters !(1)  
                cl_used(icly)=ntrk  !tag used clusters !(1)  
3354              elseif(icl.ne.0)then              elseif(icl.ne.0)then
3355  c               cl_used(icl)=1   !tag used clusters                 cl_used(icl)=ntrk   !tag used clusters
                cl_used(icl)=ntrk   !tag used clusters !1)  
3356              endif              endif
3357                            
 c               if(DEBUG)then  
 c                  print*,ip,' <<< ',id  
 c               endif  
3358  *     -----------------------------  *     -----------------------------
3359  *     remove the couple from clouds  *     remove the couple from clouds
3360  *     remove also vitual couples containing the  *     remove also vitual couples containing the
# Line 3631  c               endif Line 3405  c               endif
3405    
3406    
3407    
 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$$$  
   
3408    
3409    
3410  *     ****************************************************  *     ****************************************************
# Line 3738  c$$$ Line 3415  c$$$
3415        include 'level1.f'        include 'level1.f'
3416        include 'common_momanhough.f'        include 'common_momanhough.f'
3417        include 'level2.f'        include 'level2.f'
 c      include 'level1.f'  
3418    
3419    *     ---------------------------------
3420    *     variables initialized from level1
3421    *     ---------------------------------
3422        do i=1,nviews        do i=1,nviews
3423           good2(i)=good1(i)           good2(i)=good1(i)
3424             do j=1,nva1_view
3425                vkflag(i,j)=1
3426                if(cnnev(i,j).le.0)then
3427                   vkflag(i,j)=cnnev(i,j)
3428                endif
3429             enddo
3430        enddo        enddo
3431    *     ----------------
3432    *     level2 variables
3433    *     ----------------
3434        NTRK = 0        NTRK = 0
3435        do it=1,NTRKMAX        do it=1,NTRKMAX
3436           IMAGE(IT)=0           IMAGE(IT)=0
# Line 3755  c      include 'level1.f' Line 3441  c      include 'level1.f'
3441              ZM_nt(IP,IT) = 0              ZM_nt(IP,IT) = 0
3442              RESX_nt(IP,IT) = 0              RESX_nt(IP,IT) = 0
3443              RESY_nt(IP,IT) = 0              RESY_nt(IP,IT) = 0
3444                TAILX_nt(IP,IT) = 0
3445                TAILY_nt(IP,IT) = 0
3446                XBAD(IP,IT) = 0
3447                YBAD(IP,IT) = 0
3448              XGOOD_nt(IP,IT) = 0              XGOOD_nt(IP,IT) = 0
3449              YGOOD_nt(IP,IT) = 0              YGOOD_nt(IP,IT) = 0
3450                LS(IP,IT) = 0
3451              DEDX_X(IP,IT) = 0              DEDX_X(IP,IT) = 0
3452              DEDX_Y(IP,IT) = 0              DEDX_Y(IP,IT) = 0
3453              CLTRX(IP,IT) = 0              CLTRX(IP,IT) = 0
# Line 3812  c      include 'level1.f' Line 3503  c      include 'level1.f'
3503           alfayz1_nt(idb)=0                 alfayz1_nt(idb)=0      
3504           alfayz2_nt(idb)=0                 alfayz2_nt(idb)=0      
3505        enddo        enddo
3506        do itr=1,ntrpl_max_nt        do itr=1,ntrpt_max_nt
3507           tr_cloud_nt(itr)=0           tr_cloud_nt(itr)=0
3508           alfaxz1_nt(itr)=0                 alfaxz1_nt(itr)=0      
3509           alfaxz2_nt(itr)=0                 alfaxz2_nt(itr)=0      
# Line 3841  c      include 'level1.f' Line 3532  c      include 'level1.f'
3532          alfayz1(idb)=0                    alfayz1(idb)=0          
3533          alfayz2(idb)=0                    alfayz2(idb)=0          
3534        enddo                            enddo                    
3535        do itr=1,ntrpl_max                do itr=1,ntrpt_max        
3536          tr_cloud(itr)=0                  tr_cloud(itr)=0        
3537          cpxz1(itr)=0                      cpxz1(itr)=0            
3538          cpxz2(itr)=0                      cpxz2(itr)=0            
# Line 3889  c      include 'level1.f' Line 3580  c      include 'level1.f'
3580    
3581            
3582        include 'commontracker.f'        include 'commontracker.f'
 c      include 'level1.f'  
3583        include 'level1.f'        include 'level1.f'
3584        include 'common_momanhough.f'        include 'common_momanhough.f'
3585        include 'level2.f'        include 'level2.f'
3586        include 'common_mini_2.f'        include 'common_mini_2.f'
3587        real sinth,phi,pig              include 'calib.f'
3588    
3589          character*10 PFA
3590          common/FINALPFA/PFA
3591    
3592          real sinth,phi,pig
3593          integer ssensor,sladder
3594        pig=acos(-1.)        pig=acos(-1.)
3595    
3596    *     -------------------------------------
3597        chi2_nt(ntr)        = sngl(chi2)        chi2_nt(ntr)        = sngl(chi2)
3598        nstep_nt(ntr)       = nstep        nstep_nt(ntr)       = nstep
3599    *     -------------------------------------
3600        phi   = al(4)                  phi   = al(4)          
3601        sinth = al(3)                    sinth = al(3)            
3602        if(sinth.lt.0)then              if(sinth.lt.0)then      
# Line 3912  c      include 'level1.f' Line 3609  c      include 'level1.f'
3609       $     phi = phi + 2*pig         $     phi = phi + 2*pig  
3610        al(4) = phi                      al(4) = phi              
3611        al(3) = sinth                    al(3) = sinth            
   
3612        do i=1,5        do i=1,5
3613           al_nt(i,ntr)     = sngl(al(i))           al_nt(i,ntr)     = sngl(al(i))
3614           do j=1,5           do j=1,5
3615              coval(i,j,ntr) = sngl(cov(i,j))              coval(i,j,ntr) = sngl(cov(i,j))
3616           enddo           enddo
3617        enddo        enddo
3618          *     -------------------------------------      
3619        do ip=1,nplanes           ! loop on planes        do ip=1,nplanes           ! loop on planes
3620           xgood_nt(ip,ntr) = int(xgood(ip))           xgood_nt(ip,ntr) = int(xgood(ip))
3621           ygood_nt(ip,ntr) = int(ygood(ip))           ygood_nt(ip,ntr) = int(ygood(ip))
# Line 3928  c      include 'level1.f' Line 3624  c      include 'level1.f'
3624           zm_nt(ip,ntr)    = sngl(zm(ip))           zm_nt(ip,ntr)    = sngl(zm(ip))
3625           RESX_nt(IP,ntr)  = sngl(resx(ip))           RESX_nt(IP,ntr)  = sngl(resx(ip))
3626           RESY_nt(IP,ntr)  = sngl(resy(ip))           RESY_nt(IP,ntr)  = sngl(resy(ip))
3627             TAILX_nt(IP,ntr) = 0.
3628             TAILY_nt(IP,ntr) = 0.
3629           xv_nt(ip,ntr)    = sngl(xv(ip))           xv_nt(ip,ntr)    = sngl(xv(ip))
3630           yv_nt(ip,ntr)    = sngl(yv(ip))           yv_nt(ip,ntr)    = sngl(yv(ip))
3631           zv_nt(ip,ntr)    = sngl(zv(ip))           zv_nt(ip,ntr)    = sngl(zv(ip))
3632           axv_nt(ip,ntr)   = sngl(axv(ip))           axv_nt(ip,ntr)   = sngl(axv(ip))
3633           ayv_nt(ip,ntr)   = sngl(ayv(ip))           ayv_nt(ip,ntr)   = sngl(ayv(ip))  
3634           dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)) !(2)  c     l'avevo dimenticato!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3635           dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)) !(2)             factor = sqrt(
3636         $        tan( acos(-1.) * sngl(axv(ip)) /180. )**2 +
3637         $        tan( acos(-1.) * sngl(ayv(ip)) /180. )**2 +
3638         $        1. )
3639    c     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3640             dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)/factor)
3641             dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)/factor)  
3642        
3643           id  = CP_STORE(ip,IDCAND)           ax   = axv_nt(ip,ntr)
3644             ay   = ayv_nt(ip,ntr)
3645             bfx  = BX_STORE(ip,IDCAND)
3646             bfy  = BY_STORE(ip,IDCAND)
3647             if(ip.eq.6) ax = -1. * axv_nt(ip,ntr)
3648             if(ip.eq.6) bfy = -1. * BY_STORE(ip,IDCAND)
3649             tgtemp   = tan(ax*acos(-1.)/180.) + pmuH_h*bfy*0.00001
3650             angx     = 180.*atan(tgtemp)/acos(-1.)
3651             tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001        
3652             angy    = 180.*atan(tgtemp)/acos(-1.)
3653            
3654    c         print*,'* ',ip,bfx,bfy,angx,angy
3655    
3656             id  = CP_STORE(ip,IDCAND) ! couple id
3657           icl = CLS_STORE(ip,IDCAND)           icl = CLS_STORE(ip,IDCAND)
3658             ssensor = -1
3659             sladder = -1
3660             ssensor = SENSOR_STORE(ip,IDCAND)
3661             sladder = LADDER_STORE(ip,IDCAND)
3662             if(ip.eq.6.and.ssensor.ne.0)ssensor = 3 - ssensor !notazione paolo x align
3663             LS(IP,ntr)      = ssensor+10*sladder
3664    
3665           if(id.ne.0)then           if(id.ne.0)then
3666    c           >>> is a couple
3667              cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))              cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))
3668              cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))              cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))
3669  c            print*,ip,' ',cltrx(ip,ntr),cltry(ip,ntr)  
3670    c$$$            if(is_cp(id).ne.ssensor)
3671    c$$$     $           print*,'ERROR is sensor assignment (couple)'
3672    c$$$     $           ,is_cp(id),ssensor
3673    c$$$            if(LADDER(clx(nplanes-ip+1,icp_cp(id))).ne.sladder)
3674    c$$$     $           print*,'ERROR is ladder assignment (couple)'
3675    c$$$     $           ,LADDER(clx(nplanes-ip+1,icp_cp(id))),sladder
3676                
3677                nnnnx = npfastrips(clx(nplanes-ip+1,icp_cp(id)),PFA,angx)
3678                nnnny = npfastrips(cly(nplanes-ip+1,icp_cp(id)),PFA,angy)            
3679                xbad(ip,ntr)= nbadstrips(nnnnx,clx(nplanes-ip+1,icp_cp(id)))
3680                ybad(ip,ntr)= nbadstrips(nnnny,cly(nplanes-ip+1,icp_cp(id)))
3681    
3682                if(nsatstrips(clx(nplanes-ip+1,icp_cp(id))).gt.0)
3683         $           dedx_x(ip,ntr)=-dedx_x(ip,ntr)
3684                if(nsatstrips(cly(nplanes-ip+1,icp_cp(id))).gt.0)
3685         $           dedx_y(ip,ntr)=-dedx_y(ip,ntr)
3686    
3687           elseif(icl.ne.0)then           elseif(icl.ne.0)then
3688              if(mod(VIEW(icl),2).eq.0)cltrx(ip,ntr)=icl  c           >>> is a singlet
3689              if(mod(VIEW(icl),2).eq.1)cltry(ip,ntr)=icl  c$$$            if(LADDER(icl).ne.sladder)
3690  c            print*,ip,' ',cltrx(ip,ntr),cltry(ip,ntr)  c$$$     $           print*,'ERROR is ladder assignment (single)'
3691    c$$$     $           ,LADDER(icl),sladder
3692                if(mod(VIEW(icl),2).eq.0)then
3693                   cltrx(ip,ntr)=icl
3694                   nnnnn = npfastrips(icl,PFA,angx)
3695                   xbad(ip,ntr) = nbadstrips(nnnnn,icl)
3696                   if(nsatstrips(icl).gt.0)dedx_x(ip,ntr)=-dedx_x(ip,ntr)
3697                elseif(mod(VIEW(icl),2).eq.1)then
3698                   cltry(ip,ntr)=icl
3699                   nnnnn = npfastrips(icl,PFA,angy)
3700                   ybad(ip,ntr) = nbadstrips(nnnnn,icl)
3701                   if(nsatstrips(icl).gt.0)dedx_y(ip,ntr)=-dedx_y(ip,ntr)
3702                endif
3703           endif                     endif          
3704    
3705        enddo        enddo
3706    
3707    
3708    c$$$      print*,(xgood(i),i=1,6)
3709    c$$$      print*,(ygood(i),i=1,6)
3710    c$$$      print*,(ls(i,ntr),i=1,6)
3711    c$$$      print*,(dedx_x(i,ntr),i=1,6)
3712    c$$$      print*,(dedx_y(i,ntr),i=1,6)
3713    c$$$      print*,'-----------------------'
3714    
3715        end        end
3716    
3717        subroutine fill_level2_siglets        subroutine fill_level2_siglets
# Line 3962  c            print*,ip,' ',cltrx(ip,ntr) Line 3723  c            print*,ip,' ',cltrx(ip,ntr)
3723  *     -------------------------------------------------------  *     -------------------------------------------------------
3724    
3725        include 'commontracker.f'        include 'commontracker.f'
 c      include 'level1.f'  
3726        include 'calib.f'        include 'calib.f'
3727        include 'level1.f'        include 'level1.f'
3728        include 'common_momanhough.f'        include 'common_momanhough.f'
# Line 3970  c      include 'level1.f' Line 3730  c      include 'level1.f'
3730        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
3731    
3732  *     count #cluster per plane not associated to any track  *     count #cluster per plane not associated to any track
 c      good2=1!.true.  
3733        nclsx = 0        nclsx = 0
3734        nclsy = 0        nclsy = 0
3735    
# Line 3984  c      good2=1!.true. Line 3743  c      good2=1!.true.
3743              if(mod(VIEW(icl),2).eq.0)then !=== X views              if(mod(VIEW(icl),2).eq.0)then !=== X views
3744                 nclsx = nclsx + 1                 nclsx = nclsx + 1
3745                 planex(nclsx) = ip                 planex(nclsx) = ip
3746                 sgnlxs(nclsx) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)                 sgnlxs(nclsx) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3747                   if(nsatstrips(icl).gt.0)sgnlxs(nclsx)=-sgnlxs(nclsx)
3748                 clsx(nclsx)   = icl                 clsx(nclsx)   = icl
3749                 do is=1,2                 do is=1,2
3750  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)
3751                    call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)  c                  call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)
3752                      call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.,0.,0.)
3753                    xs(is,nclsx) = (xPAM_A+xPAM_B)/2                    xs(is,nclsx) = (xPAM_A+xPAM_B)/2
3754                 enddo                 enddo
3755  c$$$               print*,'nclsx         ',nclsx  c$$$               print*,'nclsx         ',nclsx
# Line 3999  c$$$               print*,'xs(2,nclsx)   Line 3760  c$$$               print*,'xs(2,nclsx)  
3760              else                          !=== Y views              else                          !=== Y views
3761                 nclsy = nclsy + 1                 nclsy = nclsy + 1
3762                 planey(nclsy) = ip                 planey(nclsy) = ip
3763                 sgnlys(nclsy) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)                 sgnlys(nclsy) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3764                   if(nsatstrips(icl).gt.0)sgnlys(nclsy)=-sgnlys(nclsy)
3765                 clsy(nclsy)   = icl                 clsy(nclsy)   = icl
3766                 do is=1,2                 do is=1,2
3767  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)
3768                    call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)  c                  call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)
3769                      call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.,0.,0.)
3770                    ys(is,nclsy) = (yPAM_A+yPAM_B)/2                    ys(is,nclsy) = (yPAM_A+yPAM_B)/2
3771                 enddo                 enddo
3772  c$$$               print*,'nclsy         ',nclsy  c$$$               print*,'nclsy         ',nclsy
# Line 4013  c$$$               print*,'ys(1,nclsy)   Line 3776  c$$$               print*,'ys(1,nclsy)  
3776  c$$$               print*,'ys(2,nclsy)   ',ys(2,nclsy)  c$$$               print*,'ys(2,nclsy)   ',ys(2,nclsy)
3777              endif              endif
3778           endif           endif
 c      print*,icl,cl_used(icl),cl_good(icl),ip,VIEW(icl)!nclsx(ip),nclsy(ip)  
3779    
3780  ***** LO METTO QUI PERCHE` NON SO DOVE METTERLO  ***** LO METTO QUI PERCHE` NON SO DOVE METTERLO
3781           whichtrack(icl) = cl_used(icl)           whichtrack(icl) = cl_used(icl)
Legend:
Removed from v.1.10  
changed lines
  Added in v.1.23
  ViewVC Help
Powered by ViewVC 1.1.23