/[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.8 by pam-fi, Wed Oct 25 16:18:41 2006 UTC revision 1.47 by pam-ts, Wed Jun 4 07:57:04 2014 UTC
# Line 12  Line 12 
12        subroutine track_finding(iflag)        subroutine track_finding(iflag)
13    
14        include 'commontracker.f'        include 'commontracker.f'
15          include 'level1.f'
16        include 'common_momanhough.f'        include 'common_momanhough.f'
17        include 'common_mech.f'        include 'common_mech.f'
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'
       include 'level1.f'  
21        include 'level2.f'        include 'level2.f'
22    
23        include 'momanhough_init.f'  c      print*,'======================================================'
24    c$$$      do ic=1,NCLSTR1
25    c$$$         if(.false.
26    c$$$     $        .or.nsatstrips(ic).gt.0
27    c$$$c     $        .or.nbadstrips(0,ic).gt.0
28    c$$$c     $        .or.nbadstrips(4,ic).gt.0
29    c$$$c     $        .or.nbadstrips(3,ic).gt.0
30    c$$$     $        .or..false.)then
31    c$$$            print*,'--- cl-',ic,' ------------------------'
32    c$$$            istart = INDSTART(IC)
33    c$$$            istop  = TOTCLLENGTH
34    c$$$            if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
35    c$$$            print*,'ADC   ',(CLADC(i),i=istart,istop)
36    c$$$            print*,'s/n   ',(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
37    c$$$            print*,'sgnl  ',(CLSIGNAL(i),i=istart,istop)
38    c$$$            print*,'strip ',(i-INDMAX(ic),i=istart,istop)
39    c$$$            print*,'view ',VIEW(ic)
40    c$$$            print*,'maxs ',MAXS(ic)
41    c$$$            print*,'COG4 ',cog(4,ic)
42    c$$$            ff = fbad_cog(4,ic)
43    c$$$            print*,'fbad ',ff
44    c$$$            print*,(CLBAD(i),i=istart,istop)
45    c$$$            bb=nbadstrips(0,ic)
46    c$$$            print*,'#BAD (tot)',bb
47    c$$$            bb=nbadstrips(4,ic)
48    c$$$            print*,'#BAD (4)',bb
49    c$$$            bb=nbadstrips(3,ic)
50    c$$$            print*,'#BAD (3)',bb
51    c$$$            ss=nsatstrips(ic)
52    c$$$            print*,'#saturated ',ss
53    c$$$         endif
54    c$$$      enddo
55                
 c      logical DEBUG  
 c      common/dbg/DEBUG  
   
56  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
57  *     STEP 1  *     STEP 1
58  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
# Line 44  c      common/dbg/DEBUG Line 72  c      common/dbg/DEBUG
72  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
73  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
74    
 c      iflag=0  
75        call cl_to_couples(iflag)        call cl_to_couples(iflag)
76        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
77           goto 880               !fill ntp and go to next event                       goto 880               !go to next event
78        endif        endif
79                if(ncp_tot.eq.0)goto 880  !go to next event    
 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  
 *     selezione di tracce pulite per diagnostica  
 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  
 c$$$         if(DEBUG)then  
 c$$$            do ip=1,nplanes  
 c$$$               if(ncp_plane(ip).ne.1)good2=.false.  
 c$$$            enddo  
 c$$$c            if(good2.eq.0)goto 100!next event  
 c$$$c            if(good2.eq.0)goto 880!fill ntp and go to next event  
 c$$$         endif  
           
   
   
80  *-----------------------------------------------------  *-----------------------------------------------------
81  *-----------------------------------------------------  *-----------------------------------------------------
82  *     HOUGH TRASFORM  *     HOUGH TRASFORM
# Line 91  c$$$         endif Line 105  c$$$         endif
105  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
106  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
107    
108  c      iflag=0  
109        call cp_to_doubtrip(iflag)        call cp_to_doubtrip(iflag)
110        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
111           goto 880               !fill ntp and go to next event                       goto 880               !go to next event            
112        endif        endif
113          if(ndblt.eq.0.or.ntrpt.eq.0)goto 880 !go to next event    
114                
115                
116  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
# Line 123  c      iflag=0 Line 138  c      iflag=0
138  *     $     ,ptcloud_xz,tr_cloud,cpcloud_xz            *     $     ,ptcloud_xz,tr_cloud,cpcloud_xz          
139  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
140  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
141    *     count number of hit planes
142          planehit=0                
143          do np=1,nplanes          
144            if(ncp_plane(np).ne.0)then  
145              planehit=planehit+1  
146            endif                  
147          enddo                    
148          if(planehit.lt.3) goto 880 ! exit              
149          
150          nptxz_min=x_min_start              
151          nplxz_min=x_min_start              
152                
153          nptyz_min=y_min_start              
154          nplyz_min=y_min_start              
155    
156          cutdistyz=cutystart      
157          cutdistxz=cutxstart      
158    
159  c      iflag=0   878  continue
160        call doub_to_YZcloud(iflag)        call doub_to_YZcloud(iflag)
161        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
162           goto 880               !fill ntp and go to next event                       goto 880               !fill ntp and go to next event            
163          endif
164    *     ------------------------------------------------
165    *     first try to release the tolerance
166    *     ------------------------------------------------
167          if(nclouds_yz.eq.0.and.cutdistyz.lt.maxcuty)then
168             if(cutdistyz.lt.maxcuty/2)then
169                cutdistyz=cutdistyz+cutystep
170             else
171                cutdistyz=cutdistyz+(3*cutystep)
172             endif
173             if(DEBUG.EQ.1)print*,'try again - cutdistyz ',cutdistyz
174             goto 878                
175          endif    
176    *     ------------------------------------------------
177    *     hence reduce the minimum number of plane
178    *     ------------------------------------------------
179          if(nclouds_yz.eq.0.and.nplyz_min.gt.3)then
180             nplyz_min=nplyz_min-1
181             if(DEBUG.EQ.1)print*,'try again - nplyz_min ',nplyz_min
182             goto 878
183        endif        endif
184  c      iflag=0  
185          if(nclouds_yz.eq.0)goto 880 !go to next event    
186    
187    
188    ccc   if(planehit.eq.3) goto 881    
189          
190     879  continue  
191        call trip_to_XZcloud(iflag)        call trip_to_XZcloud(iflag)
192        if(iflag.eq.1)then        !bad event        if(iflag.eq.1)then        !bad event
193           goto 880               !fill ntp and go to next event                       goto 880               !fill ntp and go to next event            
194        endif        endif
195          *     ------------------------------------------------
196    *     first try to release the tolerance
197    *     ------------------------------------------------                          
198          if(nclouds_xz.eq.0.and.cutdistxz.lt.maxcutx)then
199            cutdistxz=cutdistxz+cutxstep
200             if(DEBUG.EQ.1)print*,'try again - cutdistxz ',cutdistxz
201            goto 879                
202          endif                    
203    *     ------------------------------------------------
204    *     hence reduce the minimum number of plane
205    *     ------------------------------------------------
206          if(nclouds_xz.eq.0.and.nplxz_min.gt.3)then
207             nplxz_min=nplxz_min-1
208             if(DEBUG.EQ.1)print*,'try again - nplxz_min ',nplxz_min
209             goto 879
210          endif
211    
212          if(nclouds_xz.eq.0)goto 880 !go to next event    
213    
214        
215    c$$$ 881  continue  
216    c$$$*     if there is at least three planes on the Y view decreases cuts on X view
217    c$$$      if(nclouds_xz.eq.0.and.nclouds_yz.gt.0.and.
218    c$$$     $     nplxz_min.ne.y_min_start)then
219    c$$$        nptxz_min=x_min_step    
220    c$$$        nplxz_min=x_min_start-x_min_step              
221    c$$$        goto 879                
222    c$$$      endif                    
223            
224   880  return   880  return
225        end        end
226    
# Line 144  c      iflag=0 Line 230  c      iflag=0
230        subroutine track_fitting(iflag)        subroutine track_fitting(iflag)
231    
232        include 'commontracker.f'        include 'commontracker.f'
233          include 'level1.f'
234        include 'common_momanhough.f'        include 'common_momanhough.f'
235        include 'common_mech.f'        include 'common_mech.f'
236        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
237        include 'common_mini_2.f'        include 'common_mini_2.f'
238        include 'calib.f'        include 'calib.f'
       include 'level1.f'  
239        include 'level2.f'        include 'level2.f'
240    
241        include 'momanhough_init.f'  c      include 'momanhough_init.f'
242                
 c      logical DEBUG  
 c      common/dbg/DEBUG  
   
243        logical FIMAGE            !        logical FIMAGE            !
244          real trackimage(NTRACKSMAX)
245          real*8 AL_GUESS(5)
246    
247  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
248  *     STEP 4   (ITERATED until any other physical track isn't found)  *     STEP 4   (ITERATED until any other physical track isn't found)
# Line 184  c      common/dbg/DEBUG Line 269  c      common/dbg/DEBUG
269  *  *
270  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
271  *-------------------------------------------------------------------------------  *-------------------------------------------------------------------------------
272           ntrk=0                 !counter of identified physical tracks  ccc         ntrk=0                 !counter of identified physical tracks
273    
274  11111    continue               !<<<<<<< come here when performing a new search  c11111    continue               !<<<<<<< come here when performing a new search
275             continue                  !<<<<<<< come here when performing a new search
276    
277             if(nclouds_xz.eq.0)goto 880 !go to next event    
278             if(nclouds_yz.eq.0)goto 880 !go to next event    
279    
280  c         iflag=0  c         iflag=0
281           call clouds_to_ctrack(iflag)           call clouds_to_ctrack(iflag)
282           if(iflag.eq.1)then     !no candidate tracks found           if(iflag.eq.1)then     !no candidate tracks found
283              goto 880            !fill ntp and go to next event                goto 880            !fill ntp and go to next event  
284           endif           endif
285             if(ntracks.eq.0)goto 880 !go to next event    
286    
287           FIMAGE=.false.         !processing best track (not track image)           FIMAGE=.false.         !processing best track (not track image)
288           ibest=0                !best track among candidates           ibest=0                !best track among candidates
289           iimage=0               !track image           iimage=0               !track image
290  *     ------------- select the best track -------------  *     ------------- select the best track -------------
291    c$$$         rchi2best=1000000000.
292    c$$$         do i=1,ntracks
293    c$$$            if(RCHI2_STORE(i).lt.rchi2best.and.
294    c$$$     $         RCHI2_STORE(i).gt.0)then
295    c$$$               ibest=i
296    c$$$               rchi2best=RCHI2_STORE(i)
297    c$$$            endif
298    c$$$         enddo
299    c$$$         if(ibest.eq.0)goto 880 !>> no good candidates
300    
301    *     -------------------------------------------------------
302    *     order track-candidates according to:
303    *     1st) decreasing n.points
304    *     2nd) increasing chi**2
305    *     -------------------------------------------------------
306           rchi2best=1000000000.           rchi2best=1000000000.
307             ndofbest=0            
308           do i=1,ntracks           do i=1,ntracks
309              if(RCHI2_STORE(i).lt.rchi2best.and.             ndof=0              
310       $         RCHI2_STORE(i).gt.0)then             do ii=1,nplanes    
311                 ndof=ndof        
312         $            +int(xgood_store(ii,i))
313         $            +int(ygood_store(ii,i))
314               enddo              
315               if(ndof.gt.ndofbest)then
316                 ibest=i
317                 rchi2best=RCHI2_STORE(i)
318                 ndofbest=ndof    
319               elseif(ndof.eq.ndofbest)then
320                 if(RCHI2_STORE(i).lt.rchi2best.and.
321         $            RCHI2_STORE(i).gt.0)then
322                 ibest=i                 ibest=i
323                 rchi2best=RCHI2_STORE(i)                 rchi2best=RCHI2_STORE(i)
324              endif                 ndofbest=ndof  
325                 endif            
326               endif
327           enddo           enddo
328    
329    
330           if(ibest.eq.0)goto 880 !>> no good candidates           if(ibest.eq.0)goto 880 !>> no good candidates
331  *-------------------------------------------------------------------------------      *-------------------------------------------------------------------------------    
332  *     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 224  c         iflag=0 Line 345  c         iflag=0
345              iimage=0              iimage=0
346           endif           endif
347           if(icand.eq.0)then           if(icand.eq.0)then
348              print*,'HAI FATTO UN CASINO!!!!!! icand = ',icand              if(VERBOSE.EQ.1)then
349       $           ,ibest,iimage                 print*,'HAI FATTO UN CASINO!!!!!! icand = ',icand
350         $              ,ibest,iimage
351                endif
352              return              return
353           endif           endif
354    
# Line 236  c         iflag=0 Line 359  c         iflag=0
359  *     **********************************************************  *     **********************************************************
360  *     ************************** FIT *** FIT *** FIT *** FIT ***  *     ************************** FIT *** FIT *** FIT *** FIT ***
361  *     **********************************************************  *     **********************************************************
362             call guess()
363             do i=1,5
364                AL_GUESS(i)=AL(i)
365             enddo
366    
367           do i=1,5           do i=1,5
368              AL(i)=dble(AL_STORE(i,icand))                          AL(i)=dble(AL_STORE(i,icand))            
369           enddo           enddo
370            
371           IDCAND = icand         !fitted track-candidate           IDCAND = icand         !fitted track-candidate
372           ifail=0                !error flag in chi2 computation           ifail=0                !error flag in chi2 computation
373           jstep=0                !# minimization steps           jstep=0                !# minimization steps
374    
375           iprint=0           iprint=0
376           if(DEBUG)iprint=1  c         if(DEBUG.EQ.1)iprint=1
377             if(VERBOSE.EQ.1)iprint=1
378             if(DEBUG.EQ.1)iprint=2
379           call mini2(jstep,ifail,iprint)           call mini2(jstep,ifail,iprint)
380           if(ifail.ne.0) then  cc         if(ifail.ne.0) then
381              if(DEBUG)then           if(ifail.ne.0.or.CHI2.ne.CHI2) then !new
382                if(CHI2.ne.CHI2)CHI2=-9999. !new
383                if(VERBOSE.EQ.1)then
384                 print *,                 print *,
385       $              '*** MINIMIZATION FAILURE *** (mini2) '       $              '*** MINIMIZATION FAILURE *** (after refinement) '
386       $              ,iev       $              ,iev
387              endif              endif
             chi2=-chi2  
388           endif           endif
389                    
390           if(DEBUG)then           if(DEBUG.EQ.1)then
391              print*,'----------------------------- improved track coord'              print*,'----------------------------- improved track coord'
392  22222       format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)  22222       format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)
393              do ip=1,6              do ip=1,6
# Line 265  c         iflag=0 Line 398  c         iflag=0
398           endif           endif
399    
400  c         rchi2=chi2/dble(ndof)  c         rchi2=chi2/dble(ndof)
401           if(DEBUG)then           if(DEBUG.EQ.1)then
402              print*,' '              print*,' '
403              print*,'****** SELECTED TRACK *************'              print*,'****** SELECTED TRACK *************'
404              print*,'#         R. chi2        RIG'              print*,'#         R. chi2        RIG'
# Line 281  c         rchi2=chi2/dble(ndof) Line 414  c         rchi2=chi2/dble(ndof)
414  *     ------------- search if the track has an IMAGE -------------  *     ------------- search if the track has an IMAGE -------------
415  *     ------------- (also this is stored )           -------------  *     ------------- (also this is stored )           -------------
416           if(FIMAGE)goto 122     !>>> jump! (this is already an image)           if(FIMAGE)goto 122     !>>> jump! (this is already an image)
417  *     now search for track-image, by comparing couples IDs  
418    *     -----------------------------------------------------
419    *     first check if the track is ambiguous
420    *     -----------------------------------------------------
421    *     (modified on august 2007 by ElenaV)
422             is1=0
423             do ip=1,NPLANES
424                if(SENSOR_STORE(ip,icand).ne.0)then
425                   is1=SENSOR_STORE(ip,icand)
426                   if(ip.eq.6)is1=3-is1 !last plane inverted
427                endif
428             enddo
429             if(is1.eq.0)then
430                if(WARNING.EQ.1)print*,'** WARNING ** sensor=0'
431                goto 122            !jump
432             endif
433             do ip=1,NPLANES
434                is2 = SENSOR_STORE(ip,icand) !sensor
435                if(ip.eq.6.and.is2.ne.0)is2=3-is2 !last plane inverted
436                if(
437         $           (is1.ne.is2.and.is2.ne.0)
438         $           )goto 122      !jump (this track cannot have an image)
439             enddo
440             if(DEBUG.eq.1)print*,' >>> ambiguous track! '
441    *     ---------------------------------------------------------------
442    *     take the candidate with the greatest number of matching couples
443    *     if more than one satisfies the requirement,
444    *     choose the one with more points and lower chi2
445    *     ---------------------------------------------------------------
446    *     count the number of matching couples
447             ncpmax = 0
448             iimage   = 0           !id of candidate with better matching
449           do i=1,ntracks           do i=1,ntracks
450              iimage=i              ncp=0
451              do ip=1,nplanes              do ip=1,nplanes
452                 if(     CP_STORE(nplanes-ip+1,icand).ne.                 if(CP_STORE(nplanes-ip+1,icand).ne.0)then
453       $              -1*CP_STORE(nplanes-ip+1,i) )iimage=0                    if(
454         $                 CP_STORE(nplanes-ip+1,i).ne.0
455         $                 .and.
456         $                 CP_STORE(nplanes-ip+1,icand).eq.
457         $                 -1*CP_STORE(nplanes-ip+1,i)
458         $                 )then
459                         ncp=ncp+1  !ok
460                      elseif(
461         $                    CP_STORE(nplanes-ip+1,i).ne.0
462         $                    .and.
463         $                    CP_STORE(nplanes-ip+1,icand).ne.
464         $                    -1*CP_STORE(nplanes-ip+1,i)
465         $                    )then
466                         ncp = 0
467                         goto 117   !it is not an image candidate
468                      else
469                        
470                      endif
471                   endif
472                enddo
473     117        continue
474                trackimage(i)=ncp   !number of matching couples
475                if(ncp>ncpmax)then
476                   ncpmax=ncp
477                   iimage=i
478                endif
479             enddo
480    *     check if there are more than one image candidates
481             ngood=0
482             do i=1,ntracks
483                if( ncpmax.ne.0.and.trackimage(i).eq.ncpmax )ngood=ngood+1
484             enddo
485             if(DEBUG.eq.1)print*,' n.image-candidates : ',ngood
486    *     if there are, choose the best one
487             if(ngood.gt.0)then
488    *     -------------------------------------------------------
489    *     order track-candidates according to:
490    *     1st) decreasing n.points
491    *     2nd) increasing chi**2
492    *     -------------------------------------------------------
493                rchi2best=1000000000.
494                ndofbest=0            
495                do i=1,ntracks
496                   if( trackimage(i).eq.ncpmax )then
497                      ndof=0              
498                      do ii=1,nplanes    
499                         ndof=ndof        
500         $                    +int(xgood_store(ii,i))
501         $                    +int(ygood_store(ii,i))
502                      enddo              
503                      if(ndof.gt.ndofbest)then
504                         iimage=i
505                         rchi2best=RCHI2_STORE(i)
506                         ndofbest=ndof    
507                      elseif(ndof.eq.ndofbest)then
508                         if(RCHI2_STORE(i).lt.rchi2best.and.
509         $                    RCHI2_STORE(i).gt.0)then
510                            iimage=i
511                            rchi2best=RCHI2_STORE(i)
512                            ndofbest=ndof  
513                         endif            
514                      endif
515                   endif
516              enddo              enddo
517              if(  iimage.ne.0.and.  
518  c     $           RCHI2_STORE(i).le.CHI2MAX.and.              if(DEBUG.EQ.1)then
519  c     $           RCHI2_STORE(i).gt.0.and.                 print*,'Track candidate ',iimage
      $           .true.)then  
                if(DEBUG)print*,'Track candidate ',iimage  
520       $              ,' >>> TRACK IMAGE >>> of'       $              ,' >>> TRACK IMAGE >>> of'
521       $              ,ibest       $              ,ibest
                goto 122         !image track found  
522              endif              endif
523           enddo              
524             endif
525    
526    
527   122     continue   122     continue
528    
 *     --- and store the results --------------------------------  
          ntrk = ntrk + 1                   !counter of found tracks  
          if(.not.FIMAGE  
      $        .and.iimage.eq.0) image(ntrk)= 0  
          if(.not.FIMAGE  
      $        .and.iimage.ne.0)image(ntrk)=ntrk+1 !this is the image of the next  
          if(FIMAGE)     image(ntrk)=ntrk-1 !this is the image of the previous  
          call fill_level2_tracks(ntrk)     !==> good2=.true.  
 c         print*,'++++++++++ iimage,fimage,ntrk,image '  
 c     $        ,iimage,fimage,ntrk,image(ntrk)  
529    
530    *     --- and store the results --------------------------------
531           if(ntrk.eq.NTRKMAX)then           if(ntrk.eq.NTRKMAX)then
532              if(verbose)              if(verbose.eq.1)
533       $           print*,       $           print*,
534       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
535       $           'tracks exceeds vector dimension '       $           'tracks exceeds vector dimension '
# Line 320  c     $        ,iimage,fimage,ntrk,image Line 537  c     $        ,iimage,fimage,ntrk,image
537  cc            good2=.false.  cc            good2=.false.
538              goto 880            !fill ntp and go to next event              goto 880            !fill ntp and go to next event
539           endif           endif
540    
541             ntrk = ntrk + 1                   !counter of found tracks
542             if(.not.FIMAGE
543         $        .and.iimage.eq.0) image(ntrk)= 0
544             if(.not.FIMAGE
545         $        .and.iimage.ne.0)image(ntrk)=ntrk+1 !this is the image of the next
546             if(FIMAGE)     image(ntrk)=ntrk-1 !this is the image of the previous
547             call fill_level2_tracks(ntrk)     !==> good2=.true.
548    
549    c$$$         if(ntrk.eq.NTRKMAX)then
550    c$$$            if(verbose.eq.1)
551    c$$$     $           print*,
552    c$$$     $           '** warning ** number of identified '//
553    c$$$     $           'tracks exceeds vector dimension '
554    c$$$     $           ,'( ',NTRKMAX,' )'
555    c$$$cc            good2=.false.
556    c$$$            goto 880            !fill ntp and go to next event
557    c$$$         endif
558           if(iimage.ne.0)then           if(iimage.ne.0)then
559              FIMAGE=.true.       !              FIMAGE=.true.       !
560              goto 1212           !>>> fit image-track              goto 1212           !>>> fit image-track
561           endif           endif
562    
 *     --- then remove selected clusters (ibest+iimage) from clouds ----  
          call clean_XYclouds(ibest,iflag)  
          if(iflag.eq.1)then     !bad event  
             goto 880            !fill ntp and go to next event              
          endif  
   
 *     **********************************************************  
 *     condition to start a new search  
 *     **********************************************************  
          ixznew=0  
          do ixz=1,nclouds_xz  
             if(ptcloud_xz(ixz).ge.nptxz_min)ixznew=1  
          enddo  
          iyznew=0  
          do iyz=1,nclouds_yz  
             if(ptcloud_yz(iyz).ge.nptyz_min)iyznew=1  
          enddo  
           
          if(ixznew.ne.0.and.  
      $      iyznew.ne.0.and.  
      $        rchi2best.le.CHI2MAX.and.  
 c     $        rchi2best.lt.15..and.  
      $        .true.)then  
             if(DEBUG)then  
                print*,'***** NEW SEARCH ****'  
             endif  
             goto 11111          !try new search  
               
          endif  
 *     **********************************************  
   
   
563    
564   880     return   880     return
565        end        end
566    
567    
   
   
 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$$$  
   
568                
569  ************************************************************  ************************************************************
570  ************************************************************  ************************************************************
# Line 558  c$$$ Line 589  c$$$
589  *     PFAy   - Position Finding Algorithm in y (COG2,ETA2,...)  *     PFAy   - Position Finding Algorithm in y (COG2,ETA2,...)
590  *     angx   - Projected angle in x  *     angx   - Projected angle in x
591  *     angy   - Projected angle in y  *     angy   - Projected angle in y
592    *     bfx    - x-component of magnetci field
593    *     bfy    - y-component of magnetic field
594  *  *
595  *     --------- COUPLES -------------------------------------------------------  *     --------- COUPLES -------------------------------------------------------
596  *     The couple defines a point in the space.  *     The couple defines a point in the space.
# Line 596  c$$$ Line 629  c$$$
629  *  *
630  *  *
631    
632        subroutine xyz_PAM(icx,icy,sensor,PFAx,PFAy,angx,angy)        subroutine xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
633    
 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*****************************************************  
634                
635        include 'commontracker.f'        include 'commontracker.f'
       include 'calib.f'  
636        include 'level1.f'        include 'level1.f'
637          include 'calib.f'
638        include 'common_align.f'        include 'common_align.f'
639        include 'common_mech.f'        include 'common_mech.f'
640        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
       include 'common_resxy.f'  
   
 c      logical DEBUG  
 c      common/dbg/DEBUG  
641    
642        integer icx,icy           !X-Y cluster ID        integer icx,icy           !X-Y cluster ID
643        integer sensor        integer sensor
644        integer viewx,viewy        integer viewx,viewy
645        character*4 PFAx,PFAy     !PFA to be used        character*4 PFAx,PFAy     !PFA to be used
646        real angx,angy            !X-Y angle        real ax,ay                !X-Y geometric angle
647          real angx,angy            !X-Y effective angle
648          real bfx,bfy              !X-Y b-field components
649    
650        real stripx,stripy        real stripx,stripy
651    
652          double precision xi,yi,zi
653          double precision xi_A,yi_A,zi_A
654          double precision xi_B,yi_B,zi_B
655        double precision xrt,yrt,zrt        double precision xrt,yrt,zrt
656        double precision xrt_A,yrt_A,zrt_A        double precision xrt_A,yrt_A,zrt_A
657        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  
658                
659    
660        parameter (ndivx=30)        parameter (ndivx=30)
661    
662    
663                
664        resxPAM = 0        resxPAM = 0
665        resyPAM = 0        resyPAM = 0
666    
667        xPAM = 0.        xPAM = 0.D0
668        yPAM = 0.        yPAM = 0.D0
669        zPAM = 0.        zPAM = 0.D0
670        xPAM_A = 0.        xPAM_A = 0.D0
671        yPAM_A = 0.        yPAM_A = 0.D0
672        zPAM_A = 0.        zPAM_A = 0.D0
673        xPAM_B = 0.        xPAM_B = 0.D0
674        yPAM_B = 0.        yPAM_B = 0.D0
675        zPAM_B = 0.        zPAM_B = 0.D0
676    
677          if(sensor.lt.1.or.sensor.gt.2)then
678             print*,'xyz_PAM   ***ERROR*** wrong input '
679             print*,'sensor ',sensor
680             icx=0
681             icy=0
682          endif
683    
684  *     -----------------  *     -----------------
685  *     CLUSTER X  *     CLUSTER X
686  *     -----------------  *     -----------------      
   
687        if(icx.ne.0)then        if(icx.ne.0)then
688           viewx = VIEW(icx)  
689           nldx = nld(MAXS(icx),VIEW(icx))           viewx   = VIEW(icx)
690           nplx = npl(VIEW(icx))           nldx    = nld(MAXS(icx),VIEW(icx))
691           resxPAM = RESXAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!!           nplx    = npl(VIEW(icx))
692            c         resxPAM = RESXAV
693           stripx = float(MAXS(icx))           stripx  = float(MAXS(icx))
694           if(PFAx.eq.'COG1')then  !(1)  
695              stripx = stripx      !(1)           if(
696              resxPAM = resxPAM    !(1)       $        viewx.lt.1.or.        
697           elseif(PFAx.eq.'COG2')then       $        viewx.gt.12.or.
698              stripx = stripx + cog(2,icx)                   $        nldx.lt.1.or.
699              resxPAM = resxPAM*fbad_cog(2,icx)       $        nldx.gt.3.or.
700           elseif(PFAx.eq.'ETA2')then       $        stripx.lt.1.or.
701  c            cog2 = cog(2,icx)       $        stripx.gt.3072.or.
702  c            etacorr = pfaeta2(cog2,viewx,nldx,angx)                   $        .false.)then
703  c            stripx = stripx + etacorr              print*,'xyz_PAM   ***ERROR*** wrong input '
704              stripx = stripx + pfaeta2(icx,angx)            !(3)              print*,'icx ',icx,'view ',viewx,'nld ',nldx,'strip ',stripx
705              resxPAM = risx_eta2(angx)                       !   (4)              icx = 0
706              if(DEBUG.and.fbad_cog(2,icx).ne.1)              goto 10
      $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)  
             resxPAM = resxPAM*fbad_cog(2,icx)  
          elseif(PFAx.eq.'ETA3')then                         !(3)  
             stripx = stripx + pfaeta3(icx,angx)            !(3)  
             resxPAM = risx_eta3(angx)                       !   (4)  
             if(DEBUG.and.fbad_cog(3,icx).ne.1)              !(3)  
      $           print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)!(3)  
             resxPAM = resxPAM*fbad_cog(3,icx)               !(3)  
          elseif(PFAx.eq.'ETA4')then                         !(3)  
             stripx = stripx + pfaeta4(icx,angx)            !(3)  
             resxPAM = risx_eta4(angx)                       !   (4)  
             if(DEBUG.and.fbad_cog(4,icx).ne.1)              !(3)  
      $           print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)!(3)  
             resxPAM = resxPAM*fbad_cog(4,icx)               !(3)  
          elseif(PFAx.eq.'ETA')then                          !(3)  
             stripx = stripx + pfaeta(icx,angx)             !(3)  
             resxPAM = ris_eta(icx,angx)                     !   (4)  
             if(DEBUG.and.fbad_cog(2,icx).ne.1)              !(3)  
      $           print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)!(3)  
 c            resxPAM = resxPAM*fbad_cog(2,icx)               !(3)TEMPORANEO  
             resxPAM = resxPAM*fbad_eta(icx,angx)             !(3)(4)  
          elseif(PFAx.eq.'COG')then           !(2)  
             stripx = stripx + cog(0,icx)     !(2)          
             resxPAM = risx_cog(angx)                        !   (4)  
             resxPAM = resxPAM*fbad_cog(0,icx)!(2)  
          else  
             print*,'*** Non valid p.f.a. (x) --> ',PFAx  
707           endif           endif
708    
709    *        --------------------------
710    *        magnetic-field corrections
711    *        --------------------------
712             stripx  = stripx +  fieldcorr(viewx,bfy)      
713             angx    = effectiveangle(ax,viewx,bfy)
714            
715             call applypfa(PFAx,icx,angx,corr,res)
716             stripx  = stripx + corr
717             resxPAM = res
718    
719     10   continue
720        endif        endif
721  c      if(icy.eq.0.and.icx.ne.0)      
 c     $     print*,PFAx,icx,angx,stripx,resxPAM,'***'  
         
722  *     -----------------  *     -----------------
723  *     CLUSTER Y  *     CLUSTER Y
724  *     -----------------  *     -----------------
725    
726        if(icy.ne.0)then        if(icy.ne.0)then
727    
728           viewy = VIEW(icy)           viewy = VIEW(icy)
729           nldy = nld(MAXS(icy),VIEW(icy))           nldy = nld(MAXS(icy),VIEW(icy))
730           nply = npl(VIEW(icy))           nply = npl(VIEW(icy))
731           resyPAM = RESYAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!!  c         resyPAM = RESYAV
732             stripy = float(MAXS(icy))
733    
734             if(
735         $        viewy.lt.1.or.        
736         $        viewy.gt.12.or.
737         $        nldy.lt.1.or.
738         $        nldy.gt.3.or.
739         $        stripy.lt.1.or.
740         $        stripy.gt.3072.or.
741         $        .false.)then
742                print*,'xyz_PAM   ***ERROR*** wrong input '
743                print*,'icy ',icy,'view ',viewy,'nld ',nldy,'strip ',stripy
744                icy = 0
745                goto 20
746             endif
747    
748           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
749              print*,'xyz_PAM   ***ERROR*** invalid cluster couple!!! '              if(DEBUG.EQ.1) then
750       $           ,icx,icy                 print*,'xyz_PAM   ***ERROR*** invalid cluster couple!!! '
751         $              ,icx,icy
752                endif
753              goto 100              goto 100
754           endif           endif
           
          stripy = float(MAXS(icy))  
          if(PFAy.eq.'COG1')then !(1)  
             stripy = stripy     !(1)  
             resyPAM = resyPAM   !(1)  
          elseif(PFAy.eq.'COG2')then  
             stripy = stripy + cog(2,icy)  
             resyPAM = resyPAM*fbad_cog(2,icy)  
          elseif(PFAy.eq.'ETA2')then  
 c            cog2 = cog(2,icy)  
 c            etacorr = pfaeta2(cog2,viewy,nldy,angy)  
 c            stripy = stripy + etacorr  
             stripy = stripy + pfaeta2(icy,angy)            !(3)  
             resyPAM = risy_eta2(angy)                       !   (4)  
             resyPAM = resyPAM*fbad_cog(2,icy)  
             if(DEBUG.and.fbad_cog(2,icy).ne.1)  
      $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)  
          elseif(PFAy.eq.'ETA3')then                         !(3)  
             stripy = stripy + pfaeta3(icy,angy)            !(3)  
             resyPAM = resyPAM*fbad_cog(3,icy)               !(3)  
             if(DEBUG.and.fbad_cog(3,icy).ne.1)              !(3)  
      $           print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)!(3)  
          elseif(PFAy.eq.'ETA4')then                         !(3)  
             stripy = stripy + pfaeta4(icy,angy)            !(3)  
             resyPAM = resyPAM*fbad_cog(4,icy)               !(3)  
             if(DEBUG.and.fbad_cog(4,icy).ne.1)              !(3)  
      $           print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)!(3)  
          elseif(PFAy.eq.'ETA')then                          !(3)  
             stripy = stripy + pfaeta(icy,angy)             !(3)  
             resyPAM = ris_eta(icy,angy)                     !   (4)  
 c            resyPAM = resyPAM*fbad_cog(2,icy)              !(3)TEMPORANEO  
             resyPAM = resyPAM*fbad_eta(icy,angy)            !   (4)  
             if(DEBUG.and.fbad_cog(2,icy).ne.1)              !(3)  
      $           print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)!(3)  
          elseif(PFAy.eq.'COG')then  
             stripy = stripy + cog(0,icy)              
             resyPAM = risy_cog(angy)                        !   (4)  
 c            resyPAM = ris_eta(icy,angy)                    !   (4)  
             resyPAM = resyPAM*fbad_cog(0,icy)  
          else  
             print*,'*** Non valid p.f.a. (x) --> ',PFAx  
          endif  
755    
756    *        --------------------------
757    *        magnetic-field corrections
758    *        --------------------------
759             stripy  = stripy +  fieldcorr(viewy,bfx)      
760             angy    = effectiveangle(ay,viewy,bfx)
761            
762             call applypfa(PFAy,icy,angy,corr,res)
763             stripy  = stripy + corr
764             resyPAM = res
765    
766     20   continue
767        endif        endif
768    
769          
770  c===========================================================  c===========================================================
771  C     COUPLE  C     COUPLE
772  C===========================================================  C===========================================================
# Line 778  c     (xi,yi,zi) = mechanical coordinate Line 777  c     (xi,yi,zi) = mechanical coordinate
777  c------------------------------------------------------------------------  c------------------------------------------------------------------------
778           if(((mod(int(stripx+0.5)-1,1024)+1).le.3)           if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
779       $        .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...
780              print*,'xyz_PAM (couple):',              if(DEBUG.EQ.1) then
781       $          ' WARNING: false X strip: strip ',stripx                 print*,'xyz_PAM (couple):',
782         $              ' WARNING: false X strip: strip ',stripx
783                endif
784           endif           endif
785           xi = acoordsi(stripx,viewx)           xi = dcoordsi(stripx,viewx)
786           yi = acoordsi(stripy,viewy)           yi = dcoordsi(stripy,viewy)
787           zi = 0.           zi = 0.D0
788                    
   
789  c------------------------------------------------------------------------  c------------------------------------------------------------------------
790  c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame  c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame
791  c------------------------------------------------------------------------  c------------------------------------------------------------------------
# Line 819  c--------------------------------------- Line 819  c---------------------------------------
819           yPAM = dcoord(yrt,viewy,nldy,sensor) / 1.d4           yPAM = dcoord(yrt,viewy,nldy,sensor) / 1.d4
820           zPAM = ( zrt + z_mech_sensor(nplx,nldx,sensor)*1000. ) / 1.d4           zPAM = ( zrt + z_mech_sensor(nplx,nldx,sensor)*1000. ) / 1.d4
821    
822           xPAM_A = 0.           xPAM_A = 0.D0
823           yPAM_A = 0.           yPAM_A = 0.D0
824           zPAM_A = 0.           zPAM_A = 0.D0
825    
826           xPAM_B = 0.           xPAM_B = 0.D0
827           yPAM_B = 0.           yPAM_B = 0.D0
828           zPAM_B = 0.           zPAM_B = 0.D0
829    
830        elseif(        elseif(
831       $        (icx.ne.0.and.icy.eq.0).or.       $        (icx.ne.0.and.icy.eq.0).or.
# Line 845  C======================================= Line 845  C=======================================
845              nldx = nldy              nldx = nldy
846              viewx = viewy + 1              viewx = viewy + 1
847    
848              yi   = acoordsi(stripy,viewy)              xi = dcoordsi(0.5*(nstrips+1),viewx) !sensor center
849                yi = dcoordsi(stripy,viewy)
850                zi = 0.D0
851    
852              xi_A = edgeY_d - SiDimX/2              xi_A = edgeY_d - SiDimX/2
853              yi_A = yi              yi_A = yi
# Line 855  C======================================= Line 857  C=======================================
857              yi_B = yi              yi_B = yi
858              zi_B = 0.              zi_B = 0.
859    
 c            print*,'Y-cl ',icy,stripy,' --> ',yi  
 c            print*,xi_A,' <--> ',xi_B  
860                            
861           elseif(icx.ne.0)then           elseif(icx.ne.0)then
862  c===========================================================  c===========================================================
# Line 867  C======================================= Line 867  C=======================================
867              nldy = nldx              nldy = nldx
868              viewy = viewx - 1              viewy = viewx - 1
869    
 c            print*,'X-singlet ',icx,nplx,nldx,viewx,stripx  
 c            if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips...  
870              if(((mod(int(stripx+0.5)-1,1024)+1).le.3)              if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
871       $           .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...
872                 print*,'xyz_PAM (X-singlet):',                 if(DEBUG.EQ.1) then
873       $             ' WARNING: false X strip: strip ',stripx                    print*,'xyz_PAM (X-singlet):',
874         $                 ' WARNING: false X strip: strip ',stripx
875                   endif
876              endif              endif
877              xi   = acoordsi(stripx,viewx)  
878                xi = dcoordsi(stripx,viewx)
879                yi = dcoordsi(0.5*(nstrips+1),viewy) !sensor center
880                zi = 0.D0
881    
882              xi_A = xi              xi_A = xi
883              yi_A = edgeX_d - SiDimY/2              yi_A = edgeX_d - SiDimY/2
# Line 890  c            if((stripx.le.3).or.(stripx Line 893  c            if((stripx.le.3).or.(stripx
893                 yi_B = yi                 yi_B = yi
894              endif              endif
895    
 c            print*,'X-cl ',icx,stripx,' --> ',xi  
 c            print*,yi_A,' <--> ',yi_B  
896    
897           else           else
898                if(DEBUG.EQ.1) then
899              print *,'routine xyz_PAM ---> not properly used !!!'                 print *,'routine xyz_PAM ---> not properly used !!!'
900              print *,'icx = ',icx                 print *,'icx = ',icx
901              print *,'icy = ',icy                 print *,'icy = ',icy
902                endif
903              goto 100              goto 100
904                            
905           endif           endif
# Line 936  c     N.B. I convert angles from microra Line 938  c     N.B. I convert angles from microra
938       $        + zi_B       $        + zi_B
939       $        + dz(nplx,nldx,sensor)       $        + dz(nplx,nldx,sensor)
940    
941    
942    
943             xrt = xi
944         $        - omega(nplx,nldx,sensor)*yi
945         $        + gamma(nplx,nldx,sensor)*zi
946         $        + dx(nplx,nldx,sensor)
947            
948             yrt = omega(nplx,nldx,sensor)*xi
949         $        + yi
950         $        - beta(nplx,nldx,sensor)*zi
951         $        + dy(nplx,nldx,sensor)
952            
953             zrt = -gamma(nplx,nldx,sensor)*xi
954         $        + beta(nplx,nldx,sensor)*yi
955         $        + zi
956         $        + dz(nplx,nldx,sensor)
957    
958    
959                    
960  c      xrt = xi  c      xrt = xi
961  c      yrt = yi  c      yrt = yi
# Line 946  c     (xPAM,yPAM,zPAM) = measured coordi Line 966  c     (xPAM,yPAM,zPAM) = measured coordi
966  c                        in PAMELA reference system  c                        in PAMELA reference system
967  c------------------------------------------------------------------------  c------------------------------------------------------------------------
968    
969           xPAM = 0.           xPAM = dcoord(xrt,viewx,nldx,sensor) / 1.d4
970           yPAM = 0.           yPAM = dcoord(yrt,viewy,nldy,sensor) / 1.d4
971           zPAM = 0.           zPAM = ( zrt + z_mech_sensor(nplx,nldx,sensor)*1000. ) / 1.d4
972    c$$$         xPAM = 0.D0
973    c$$$         yPAM = 0.D0
974    c$$$         zPAM = 0.D0
975    
976           xPAM_A = dcoord(xrt_A,viewx,nldx,sensor) / 1.d4           xPAM_A = dcoord(xrt_A,viewx,nldx,sensor) / 1.d4
977           yPAM_A = dcoord(yrt_A,viewy,nldy,sensor) / 1.d4           yPAM_A = dcoord(yrt_A,viewy,nldy,sensor) / 1.d4
# Line 959  c--------------------------------------- Line 982  c---------------------------------------
982           zPAM_B = ( zrt_B + z_mech_sensor(nplx,nldx,sensor)*1000.)/ 1.d4           zPAM_B = ( zrt_B + z_mech_sensor(nplx,nldx,sensor)*1000.)/ 1.d4
983                    
984    
 c         print*,'A-(',xPAM_A,yPAM_A,') B-(',xPAM_B,yPAM_B,')'  
985    
986        else        else
987                       if(DEBUG.EQ.1) then
988           print *,'routine xyz_PAM ---> not properly used !!!'              print *,'routine xyz_PAM ---> not properly used !!!'
989           print *,'icx = ',icx              print *,'icx = ',icx
990           print *,'icy = ',icy              print *,'icy = ',icy
991                         endif
992        endif        endif
993                    
994    
995    
996   100  continue   100  continue
997        end        end
998    
999    ************************************************************************
1000    *     Call xyz_PAM subroutine with default PFA and fill the mini2 common.
1001    *     (done to be called from c/c++)
1002    ************************************************************************
1003    
1004          subroutine xyzpam(ip,icx,icy,lad,sensor,ax,ay,bfx,bfy)
1005    
1006          include 'commontracker.f'
1007          include 'level1.f'
1008          include 'common_mini_2.f'
1009          include 'common_xyzPAM.f'
1010          include 'common_mech.f'
1011          include 'calib.f'
1012          
1013    *     flag to chose PFA
1014    c$$$      character*10 PFA
1015    c$$$      common/FINALPFA/PFA
1016    
1017          integer icx,icy           !X-Y cluster ID
1018          integer sensor
1019          character*4 PFAx,PFAy     !PFA to be used
1020          real ax,ay                !X-Y geometric angle
1021          real bfx,bfy              !X-Y b-field components
1022    
1023          ipx=0
1024          ipy=0      
1025          
1026    c$$$      PFAx = 'COG4'!PFA
1027    c$$$      PFAy = 'COG4'!PFA
1028    
1029    
1030          if(icx.gt.nclstr1.or.icy.gt.nclstr1)then
1031                print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1032         $           ,' do not exists (n.clusters=',nclstr1,')'
1033                icx = -1*icx
1034                icy = -1*icy
1035                return
1036            
1037          endif
1038          
1039          call idtoc(pfaid,PFAx)
1040          call idtoc(pfaid,PFAy)
1041    
1042          
1043          if(icx.ne.0.and.icy.ne.0)then
1044    
1045             ipx=npl(VIEW(icx))
1046             ipy=npl(VIEW(icy))
1047            
1048             if( (nplanes-ipx+1).ne.ip )then
1049                print*,'xyzpam: ***WARNING*** cluster icx=',icx
1050         $           ,' belongs to plane ',(nplanes-ipx+1)            
1051         $           ,' and not ',ip
1052                icx = -1*icx
1053                return
1054             endif
1055             if( (nplanes-ipy+1).ne.ip )then
1056                print*,'xyzpam: ***WARNING*** cluster icy=',icy
1057         $           ,' belongs to plane ',(nplanes-ipy+1)            
1058         $           ,' and not ',ip
1059                 icy = -1*icy
1060                return
1061             endif
1062    
1063             call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
1064    
1065             xgood(ip) = 1.
1066             ygood(ip) = 1.
1067             resx(ip)  = resxPAM
1068             resy(ip)  = resyPAM
1069    
1070             xm(ip) = xPAM
1071             ym(ip) = yPAM
1072             zm(ip) = zPAM
1073             xm_A(ip) = 0.D0
1074             ym_A(ip) = 0.D0
1075             zm_A(ip) = 0.D0
1076             xm_B(ip) = 0.D0
1077             ym_B(ip) = 0.D0
1078             zm_B(ip) = 0.D0
1079    
1080    c         zv(ip) = zPAM
1081    
1082          elseif(icx.eq.0.and.icy.ne.0)then
1083    
1084             ipy=npl(VIEW(icy))
1085             if( (nplanes-ipy+1).ne.ip )then
1086                print*,'xyzpam: ***WARNING*** cluster icy=',icy
1087         $           ,' belongs to plane ',(nplanes-ipy+1)            
1088         $           ,' and not ',ip
1089                icy = -1*icy
1090                return
1091             endif
1092    
1093             call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
1094            
1095             xgood(ip) = 0.
1096             ygood(ip) = 1.
1097             resx(ip)  = 1000.
1098             resy(ip)  = resyPAM
1099    
1100    c$$$         xm(ip) = -100.
1101    c$$$         ym(ip) = -100.
1102    c$$$         zm(ip) = (zPAM_A+zPAM_B)/2.
1103             xm(ip) = xPAM
1104             ym(ip) = yPAM
1105             zm(ip) = zPAM
1106             xm_A(ip) = xPAM_A
1107             ym_A(ip) = yPAM_A
1108             zm_A(ip) = zPAM_A
1109             xm_B(ip) = xPAM_B
1110             ym_B(ip) = yPAM_B
1111             zm_B(ip) = zPAM_B
1112    
1113    c         zv(ip) = (zPAM_A+zPAM_B)/2.
1114            
1115          elseif(icx.ne.0.and.icy.eq.0)then
1116    
1117             ipx=npl(VIEW(icx))
1118    
1119             if( (nplanes-ipx+1).ne.ip )then
1120                print*,'xyzpam: ***WARNING*** cluster icx=',icx
1121         $           ,' belongs to plane ',(nplanes-ipx+1)            
1122         $           ,' and not ',ip
1123                icx = -1*icx
1124                return
1125             endif
1126    
1127             call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
1128          
1129             xgood(ip) = 1.
1130             ygood(ip) = 0.
1131             resx(ip)  = resxPAM
1132             resy(ip)  = 1000.
1133    
1134    c$$$         xm(ip) = -100.
1135    c$$$         ym(ip) = -100.
1136    c$$$         zm(ip) = (zPAM_A+zPAM_B)/2.
1137             xm(ip) = xPAM
1138             ym(ip) = yPAM
1139             zm(ip) = zPAM
1140             xm_A(ip) = xPAM_A
1141             ym_A(ip) = yPAM_A
1142             zm_A(ip) = zPAM_A
1143             xm_B(ip) = xPAM_B
1144             ym_B(ip) = yPAM_B
1145             zm_B(ip) = zPAM_B
1146            
1147    c         zv(ip) = (zPAM_A+zPAM_B)/2.
1148    
1149          else
1150    
1151             il = 2
1152             if(lad.ne.0)il=lad
1153             is = 1
1154             if(sensor.ne.0)is=sensor
1155    
1156             xgood(ip) = 0.
1157             ygood(ip) = 0.
1158             resx(ip)  = 1000.
1159             resy(ip)  = 1000.
1160    
1161             xm(ip) = -100.
1162             ym(ip) = -100.          
1163             zm(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
1164             xm_A(ip) = 0.
1165             ym_A(ip) = 0.
1166             zm_A(ip) = 0.
1167             xm_B(ip) = 0.
1168             ym_B(ip) = 0.
1169             zm_B(ip) = 0.
1170    
1171    c         zv(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
1172    
1173          endif
1174    
1175          if(DEBUG.EQ.1)then
1176    22222    format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)
1177             write(*,22222)ip,zm(ip),xm(ip),ym(ip)
1178         $        ,xm_A(ip),ym_A(ip),xm_B(ip),ym_B(ip)
1179         $        ,xgood(ip),ygood(ip),resx(ip),resy(ip)
1180          endif
1181          end
1182    
1183  ********************************************************************************  ********************************************************************************
1184  ********************************************************************************  ********************************************************************************
# Line 993  c         print*,'A-(',xPAM_A,yPAM_A,') Line 1200  c         print*,'A-(',xPAM_A,yPAM_A,')
1200  *      *    
1201  ********************************************************************************  ********************************************************************************
1202    
1203        real function distance_to(XPP,YPP)        real function distance_to(rXPP,rYPP)
1204    
1205        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
1206    
# Line 1002  c         print*,'A-(',xPAM_A,yPAM_A,') Line 1209  c         print*,'A-(',xPAM_A,yPAM_A,')
1209  *     ( i.e. distance/resolution )  *     ( i.e. distance/resolution )
1210  *     -----------------------------------  *     -----------------------------------
1211    
1212          real rXPP,rYPP
1213          double precision XPP,YPP
1214        double precision distance,RE        double precision distance,RE
1215        double precision BETA,ALFA,xmi,ymi        double precision BETA,ALFA,xmi,ymi
1216    
1217          XPP=DBLE(rXPP)
1218          YPP=DBLE(rYPP)
1219    
1220  *     ----------------------  *     ----------------------
1221        if (        if (
1222       +     xPAM.eq.0.and.  c$$$     +     xPAM.eq.0.and.
1223       +     yPAM.eq.0.and.  c$$$     +     yPAM.eq.0.and.
1224       +     zPAM.eq.0.and.  c$$$     +     zPAM.eq.0.and.
1225       +     xPAM_A.ne.0.and.       +     xPAM_A.ne.0.and.
1226       +     yPAM_A.ne.0.and.       +     yPAM_A.ne.0.and.
1227       +     zPAM_A.ne.0.and.       +     zPAM_A.ne.0.and.
# Line 1047  c         print*,'A-(',xPAM_A,yPAM_A,') Line 1259  c         print*,'A-(',xPAM_A,yPAM_A,')
1259           endif                   endif        
1260    
1261           distance=           distance=
1262       $        ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2       $       ((xmi-XPP)**2+(ymi-YPP)**2)!QUIQUI
1263    cc     $        ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2
1264           distance=dsqrt(distance)                               distance=dsqrt(distance)                    
1265    
 c$$$         print*,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b  
 c$$$     $        ,' --- distance_to --- ',xpp,ypp  
 c$$$         print*,' resolution ',re  
1266    
1267                    
1268  *     ----------------------  *     ----------------------
1269        elseif(        elseif(
1270       +     xPAM.ne.0.and.  c$$$     +     xPAM.ne.0.and.
1271       +     yPAM.ne.0.and.  c$$$     +     yPAM.ne.0.and.
1272       +     zPAM.ne.0.and.  c$$$     +     zPAM.ne.0.and.
1273       +     xPAM_A.eq.0.and.       +     xPAM_A.eq.0.and.
1274       +     yPAM_A.eq.0.and.       +     yPAM_A.eq.0.and.
1275       +     zPAM_A.eq.0.and.       +     zPAM_A.eq.0.and.
# Line 1072  c$$$         print*,' resolution ',re Line 1282  c$$$         print*,' resolution ',re
1282  *     ----------------------  *     ----------------------
1283                    
1284           distance=           distance=
1285       $        ((xPAM-XPP)/resxPAM)**2       $       ((xPAM-XPP))**2 !QUIQUI
1286       $        +       $       +
1287       $        ((yPAM-YPP)/resyPAM)**2       $       ((yPAM-YPP))**2
1288    c$$$     $        ((xPAM-XPP)/resxPAM)**2
1289    c$$$     $        +
1290    c$$$     $        ((yPAM-YPP)/resyPAM)**2
1291           distance=dsqrt(distance)                               distance=dsqrt(distance)                    
1292    
 c$$$         print*,xPAM,yPAM,zPAM  
 c$$$     $        ,' --- distance_to --- ',xpp,ypp  
 c$$$         print*,' resolution ',resxPAM,resyPAM  
1293                    
1294        else        else
1295                    
          print*  
      $        ,' function distance_to ---> wrong usage!!!'  
          print*,' xPAM,yPAM,zPAM ',xPAM,yPAM,zPAM  
          print*,' xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b '  
      $        ,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b  
1296        endif          endif  
1297    
1298        distance_to = sngl(distance)        distance_to = sngl(distance)
# Line 1127  c$$$         print*,' resolution ',resxP Line 1332  c$$$         print*,' resolution ',resxP
1332        data c1/1.,0.,0.,1./        data c1/1.,0.,0.,1./
1333        data c2/1.,-1.,-1.,1./        data c2/1.,-1.,-1.,1./
1334        data c3/1.,1.,0.,0./        data c3/1.,1.,0.,0./
1335        real*8 yvvv,xvvv        double precision yvvv,xvvv
1336        double precision xi,yi,zi        double precision xi,yi,zi
1337        double precision xrt,yrt,zrt        double precision xrt,yrt,zrt
1338        real AA,BB        real AA,BB
1339        real yvv(4),xvv(4)        double precision yvv(4),xvv(4)
1340    
1341  *     tollerance to consider the track inside the sensitive area  *     tollerance to consider the track inside the sensitive area
1342        real ptoll        real ptoll
1343        data ptoll/150./          !um        data ptoll/150./          !um
1344    
1345        external nviewx,nviewy,acoordsi,dcoord        external nviewx,nviewy,dcoordsi,dcoord
1346    
1347        nplpt = nplPAM            !plane        nplpt = nplPAM            !plane
1348        viewx = nviewx(nplpt)        viewx = nviewx(nplpt)
# Line 1152  c$$$         print*,' resolution ',resxP Line 1357  c$$$         print*,' resolution ',resxP
1357  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1358  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame  c     (xi,yi,zi) = mechanical coordinates in the silicon sensor frame
1359  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1360                 if(((mod(int(stripx+0.5)-1,1024)+1).le.3)                 xi = dcoordsi(stripx,viewx)
1361       $              .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...                 yi = dcoordsi(stripy,viewy)
1362  c     if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips...                 zi = 0.D0
                   print*,'whichsensor: ',  
      $                ' WARNING: false X strip: strip ',stripx  
                endif  
                xi = acoordsi(stripx,viewx)  
                yi = acoordsi(stripy,viewy)  
                zi = 0.  
1363  c------------------------------------------------------------------------  c------------------------------------------------------------------------
1364  c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame  c     (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame
1365  c------------------------------------------------------------------------  c------------------------------------------------------------------------
# Line 1185  c--------------------------------------- Line 1384  c---------------------------------------
1384    
1385                 yvv(iv)=sngl(yvvv)                 yvv(iv)=sngl(yvvv)
1386                 xvv(iv)=sngl(xvvv)                 xvv(iv)=sngl(xvvv)
 c               print*,'LADDER ',il,' SENSOR ',is,' vertexes >> '  
 c     $              ,iv,xvv(iv),yvv(iv)  
1387              enddo               !end loop on sensor vertexes              enddo               !end loop on sensor vertexes
1388    
1389              dtot=0.              dtot=0.
# Line 1194  c     $              ,iv,xvv(iv),yvv(iv) Line 1391  c     $              ,iv,xvv(iv),yvv(iv)
1391                 iv1=iside                 iv1=iside
1392                 iv2=mod(iside,4)+1                 iv2=mod(iside,4)+1
1393  *     straight line passing trhough two consecutive vertexes  *     straight line passing trhough two consecutive vertexes
1394                 AA = (yvv(iv1)-yvv(iv2))/(xvv(iv1)-xvv(iv2))                 AA = REAL((yvv(iv1)-yvv(iv2))/(xvv(iv1)-xvv(iv2))) !EM GCC4.7
1395                 BB = yvv(iv1) - AA*xvv(iv1)                 BB = REAL(yvv(iv1) - AA*xvv(iv1)) !EM GCC4.7
1396  *     point along the straight line closer to the track  *     point along the straight line closer to the track
1397                 xoo = (xPAM+AA*yPAM-AA*BB)/(1+AA**2)                 xoo = (xPAM+AA*yPAM-AA*BB)/(1+AA**2)
1398                 yoo = AA*xoo + BB                 yoo = AA*xoo + BB
# Line 1207  c     $              ,iv,xvv(iv),yvv(iv) Line 1404  c     $              ,iv,xvv(iv),yvv(iv)
1404                 iv1=iside                 iv1=iside
1405                 iv2=mod(iside,4)+1                 iv2=mod(iside,4)+1
1406  *     straight line passing trhough two consecutive vertexes  *     straight line passing trhough two consecutive vertexes
1407                 AA = (xvv(iv1)-xvv(iv2))/(yvv(iv1)-yvv(iv2))                 AA = REAL((xvv(iv1)-xvv(iv2))/(yvv(iv1)-yvv(iv2))) !EM GCC4.7
1408                 BB = xvv(iv1) - AA*yvv(iv1)                 BB = REAL(xvv(iv1) - AA*yvv(iv1)) !EM GCC4.7
1409  *     point along the straight line closer to the track  *     point along the straight line closer to the track
1410                 yoo = (yPAM+AA*xPAM-AA*BB)/(1+AA**2)                 yoo = (yPAM+AA*xPAM-AA*BB)/(1+AA**2)
1411                 xoo = AA*yoo + BB                 xoo = AA*yoo + BB
# Line 1285  c     $              ,iv,xvv(iv),yvv(iv) Line 1482  c     $              ,iv,xvv(iv),yvv(iv)
1482  *     it returns the plane number  *     it returns the plane number
1483  *      *    
1484        include 'commontracker.f'        include 'commontracker.f'
1485          include 'level1.f'
1486  c      include 'common_analysis.f'  c      include 'common_analysis.f'
1487        include 'common_momanhough.f'        include 'common_momanhough.f'
1488                
# Line 1310  c      include 'common_analysis.f' Line 1508  c      include 'common_analysis.f'
1508        is_cp=0        is_cp=0
1509        if(id.lt.0)is_cp=1        if(id.lt.0)is_cp=1
1510        if(id.gt.0)is_cp=2        if(id.gt.0)is_cp=2
       if(id.eq.0)print*,'IS_CP ===> wrong couple id !!!'  
1511    
1512        return        return
1513        end        end
# Line 1322  c      include 'common_analysis.f' Line 1519  c      include 'common_analysis.f'
1519  *     it returns the id number ON THE PLANE  *     it returns the id number ON THE PLANE
1520  *      *    
1521        include 'commontracker.f'        include 'commontracker.f'
1522          include 'level1.f'
1523  c      include 'common_analysis.f'  c      include 'common_analysis.f'
1524        include 'common_momanhough.f'        include 'common_momanhough.f'
1525                
# Line 1350  c      include 'common_analysis.f' Line 1548  c      include 'common_analysis.f'
1548  *     positive if sensor =2  *     positive if sensor =2
1549  *  *
1550        include 'commontracker.f'        include 'commontracker.f'
1551          include 'level1.f'
1552  c      include 'calib.f'  c      include 'calib.f'
1553  c      include 'level1.f'  c      include 'level1.f'
1554  c      include 'common_analysis.f'  c      include 'common_analysis.f'
# Line 1379  c      include 'common_analysis.f' Line 1578  c      include 'common_analysis.f'
1578  *************************************************************************  *************************************************************************
1579  *************************************************************************  *************************************************************************
1580  *************************************************************************  *************************************************************************
 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  
1581                
1582    
1583  ***************************************************  ***************************************************
# Line 1661  c$$$      end Line 1592  c$$$      end
1592        subroutine cl_to_couples(iflag)        subroutine cl_to_couples(iflag)
1593    
1594        include 'commontracker.f'        include 'commontracker.f'
1595          include 'level1.f'
1596        include 'common_momanhough.f'        include 'common_momanhough.f'
1597        include 'momanhough_init.f'  c      include 'momanhough_init.f'
1598        include 'calib.f'        include 'calib.f'
1599        include 'level1.f'  c      include 'level1.f'
   
 c      logical DEBUG  
 c      common/dbg/DEBUG  
1600    
1601  *     output flag  *     output flag
1602  *     --------------  *     --------------
# Line 1678  c      common/dbg/DEBUG Line 1607  c      common/dbg/DEBUG
1607    
1608        integer badseed,badclx,badcly        integer badseed,badclx,badcly
1609    
1610          iflag = iflag
1611          if(DEBUG.EQ.1)print*,'cl_to_couples:'
1612    
1613    cc      if(RECOVER_SINGLETS.and..not.SECOND_SEARCH)goto 80
1614    
1615  *     init variables  *     init variables
       ncp_tot=0  
1616        do ip=1,nplanes        do ip=1,nplanes
1617           do ico=1,ncouplemax           do ico=1,ncouplemax
1618              clx(ip,ico)=0              clx(ip,ico)=0
# Line 1692  c      common/dbg/DEBUG Line 1625  c      common/dbg/DEBUG
1625           ncls(ip)=0           ncls(ip)=0
1626        enddo        enddo
1627        do icl=1,nclstrmax_level2        do icl=1,nclstrmax_level2
1628           cl_single(icl)=1           cl_single(icl) = 1     !all are single
1629           cl_good(icl)=0           cl_good(icl)   = 0     !all are bad
1630          enddo
1631          do iv=1,nviews
1632             ncl_view(iv)  = 0
1633             mask_view(iv) = 0      !all included
1634        enddo        enddo
1635                
1636    *     count number of cluster per view
1637          do icl=1,nclstr1
1638             ncl_view(VIEW(icl)) = ncl_view(VIEW(icl)) + 1        
1639          enddo
1640    *     mask views with too many clusters
1641          do iv=1,nviews
1642             if( ncl_view(iv).gt. nclusterlimit)then
1643    c            mask_view(iv) = 1
1644                mask_view(iv) = mask_view(iv) + 2**0
1645                if(DEBUG.EQ.1)
1646         $        print*,' * WARNING * cl_to_couple: n.clusters > '
1647         $        ,nclusterlimit,' on view ', iv,' --> masked!'
1648             endif
1649          enddo
1650    
1651    
1652  *     start association  *     start association
1653        ncouples=0        ncouples=0
1654        do icx=1,nclstr1          !loop on cluster (X)        do icx=1,nclstr1          !loop on cluster (X)
1655           if(mod(VIEW(icx),2).eq.1)goto 10           if(mod(VIEW(icx),2).eq.1)goto 10
1656                    
1657             if(cl_used(icx).ne.0)goto 10
1658    
1659  *     ----------------------------------------------------  *     ----------------------------------------------------
1660  *     cut on charge (X VIEW)  *     jump masked views (X VIEW)
1661  *     ----------------------------------------------------  *     ----------------------------------------------------
1662           if(dedx(icx).lt.dedx_x_min)then           if( mask_view(VIEW(icx)).ne.0 ) goto 10
             cl_single(icx)=0  
             goto 10  
          endif  
1663  *     ----------------------------------------------------  *     ----------------------------------------------------
1664  *     cut on multiplicity (X VIEW)  *     cut on charge (X VIEW)
1665  *     ----------------------------------------------------  *     ----------------------------------------------------
1666           if(mult(icx).ge.mult_x_max)then           if(sgnl(icx)/mip(VIEW(icx),LADDER(icx)).lt.dedx_x_min)then
1667              cl_single(icx)=0              cl_single(icx)=0
1668              goto 10              goto 10
1669           endif           endif
# Line 1752  c     endif Line 1704  c     endif
1704                    
1705           do icy=1,nclstr1       !loop on cluster (Y)           do icy=1,nclstr1       !loop on cluster (Y)
1706              if(mod(VIEW(icy),2).eq.0)goto 20              if(mod(VIEW(icy),2).eq.0)goto 20
1707    
1708                if(cl_used(icx).ne.0)goto 20
1709                            
1710  *     ----------------------------------------------------  *     ----------------------------------------------------
1711  *     cut on charge (Y VIEW)  *     jump masked views (Y VIEW)
1712  *     ----------------------------------------------------  *     ----------------------------------------------------
1713              if(dedx(icy).lt.dedx_y_min)then              if( mask_view(VIEW(icy)).ne.0 ) goto 20
1714                 cl_single(icy)=0  
                goto 20  
             endif  
1715  *     ----------------------------------------------------  *     ----------------------------------------------------
1716  *     cut on multiplicity (X VIEW)  *     cut on charge (Y VIEW)
1717  *     ----------------------------------------------------  *     ----------------------------------------------------
1718              if(mult(icy).ge.mult_y_max)then              if(sgnl(icy)/mip(VIEW(icy),LADDER(icy)) .lt.dedx_y_min)then
1719                 cl_single(icy)=0                 cl_single(icy)=0
1720                 goto 20                 goto 20
1721              endif              endif
# Line 1809  c     endif Line 1761  c     endif
1761  *     charge correlation  *     charge correlation
1762  *     (modified to be applied only below saturation... obviously)  *     (modified to be applied only below saturation... obviously)
1763    
1764  *     -------------------------------------------------------------                 if(  .not.(sgnl(icy).gt.chsaty.and.sgnl(icx).gt.chsatx)
 *     >>> eliminata (TEMPORANEAMENTE) la correlazione di carica <<<  
 *     -------------------------------------------------------------  
                if(  .not.(dedx(icy).gt.chsaty.and.dedx(icx).gt.chsatx)  
1765       $              .and.       $              .and.
1766       $              .not.(dedx(icy).lt.chmipy.and.dedx(icx).lt.chmipx)       $              .not.(sgnl(icy).lt.chmipy.and.sgnl(icx).lt.chmipx)
1767       $              .and.       $              .and.
1768       $              (badclx.eq.1.and.badcly.eq.1)       $              (badclx.eq.1.and.badcly.eq.1)
1769       $              .and.       $              .and.
1770       $              .true.)then       $              .true.)then
1771    
1772                    ddd=(dedx(icy)                    ddd=(sgnl(icy)
1773       $                 -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))       $                 -kch(nplx,nldx)*sgnl(icx)-cch(nplx,nldx))
1774                    ddd=ddd/sqrt(kch(nplx,nldx)**2+1)                    ddd=ddd/sqrt(kch(nplx,nldx)**2+1)
1775    
1776  c                  cut = chcut * sch(nplx,nldx)  c                  cut = chcut * sch(nplx,nldx)
1777    
1778                    sss=(kch(nplx,nldx)*dedx(icy)+dedx(icx)                    sss=(kch(nplx,nldx)*sgnl(icy)+sgnl(icx)
1779       $                 -kch(nplx,nldx)*cch(nplx,nldx))       $                 -kch(nplx,nldx)*cch(nplx,nldx))
1780                    sss=sss/sqrt(kch(nplx,nldx)**2+1)                    sss=sss/sqrt(kch(nplx,nldx)**2+1)
1781                    cut = chcut * (16 + sss/50.)                    cut = chcut * (16 + sss/50.)
# Line 1835  c                  cut = chcut * sch(npl Line 1784  c                  cut = chcut * sch(npl
1784                       goto 20    !charge not consistent                       goto 20    !charge not consistent
1785                    endif                    endif
1786                 endif                 endif
1787                  
1788  *     ------------------> COUPLE <------------------                 if(ncp_plane(nplx).gt.ncouplemax)then
1789  *     check to do not overflow vector dimentions                    if(verbose.eq.1)print*,
 c$$$               if(ncp_plane(nplx).gt.ncouplemax)then  
 c$$$                  if(DEBUG)print*,  
 c$$$     $                    ' ** warning ** number of identified'//  
 c$$$     $                    ' couples on plane ',nplx,  
 c$$$     $                    ' exceeds vector dimention'//  
 c$$$     $                    ' ( ',ncouplemax,' )'  
 c$$$c     good2=.false.  
 c$$$c     goto 880   !fill ntp and go to next event  
 c$$$                  iflag=1  
 c$$$                  return  
 c$$$               endif  
                 
                if(ncp_plane(nplx).eq.ncouplemax)then  
                   if(verbose)print*,  
1790       $                 '** warning ** number of identified '//       $                 '** warning ** number of identified '//
1791       $                 'couples on plane ',nplx,       $                 'couples on plane ',nplx,
1792       $                 'exceeds vector dimention '       $                 'exceeds vector dimention '
1793       $                 ,'( ',ncouplemax,' )'       $                 ,'( ',ncouplemax,' ) --> masked!'
1794  c     good2=.false.  c                  mask_view(nviewx(nplx)) = 2
1795  c     goto 880   !fill ntp and go to next event                      c                  mask_view(nviewy(nply)) = 2
1796                    iflag=1                   mask_view(nviewx(nplx))= mask_view(nviewx(nplx))+ 2**1
1797                    return                   mask_view(nviewy(nply))= mask_view(nviewy(nply))+ 2**1
1798                      goto 10
1799                 endif                 endif
1800                                
1801    *     ------------------> COUPLE <------------------
1802                 ncp_plane(nplx) = ncp_plane(nplx) + 1                 ncp_plane(nplx) = ncp_plane(nplx) + 1
1803                 clx(nplx,ncp_plane(nplx))=icx                 clx(nplx,ncp_plane(nplx))=icx
1804                 cly(nply,ncp_plane(nplx))=icy                 cly(nply,ncp_plane(nplx))=icy
1805                 cl_single(icx)=0                 cl_single(icx)=0
1806                 cl_single(icy)=0                 cl_single(icy)=0
             endif                                
1807  *     ----------------------------------------------  *     ----------------------------------------------
1808    
1809                endif                              
1810    
1811   20         continue   20         continue
1812           enddo                  !end loop on clusters(Y)           enddo                  !end loop on clusters(Y)
1813                    
1814   10      continue   10      continue
1815        enddo                     !end loop on clusters(X)        enddo                     !end loop on clusters(X)
1816                
         
1817        do icl=1,nclstr1        do icl=1,nclstr1
1818           if(cl_single(icl).eq.1)then           if(cl_single(icl).eq.1)then
1819              ip=npl(VIEW(icl))              ip=npl(VIEW(icl))
# Line 1884  c     goto 880   !fill ntp and go to nex Line 1821  c     goto 880   !fill ntp and go to nex
1821              cls(ip,ncls(ip))=icl              cls(ip,ncls(ip))=icl
1822           endif           endif
1823        enddo        enddo
1824    
1825    c 80   continue
1826          continue
1827                
1828                
1829        if(DEBUG)then        if(DEBUG.EQ.1)then
1830           print*,'clusters  ',nclstr1           print*,'clusters  ',nclstr1
1831           print*,'good    ',(cl_good(i),i=1,nclstr1)           print*,'good    ',(cl_good(i),i=1,nclstr1)
1832           print*,'singles ',(cl_single(i),i=1,nclstr1)           print*,'used    ',(cl_used(i),i=1,nclstr1)
1833             print*,'singlets',(cl_single(i),i=1,nclstr1)
1834           print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)           print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)
1835        endif        endif
         
       do ip=1,6  
          ncp_tot=ncp_tot+ncp_plane(ip)  
       enddo  
 c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)  
         
       if(ncp_tot.gt.ncp_max)then  
          if(verbose)print*,  
      $           '** warning ** number of identified '//  
      $           'couples exceeds upper limit for Hough tr. '  
      $           ,'( ',ncp_max,' )'              
 c            good2=.false.  
 c     goto 880       !fill ntp and go to next event  
          iflag=1  
          return  
       endif  
         
       return  
       end  
         
 ***************************************************  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 *                                                 *  
 **************************************************  
       subroutine cl_to_couples_nocharge(iflag)  
   
       include 'commontracker.f'  
       include 'common_momanhough.f'  
       include 'momanhough_init.f'  
       include 'calib.f'  
       include 'level1.f'  
   
 c      logical DEBUG  
 c      common/dbg/DEBUG  
1836    
1837  *     output flag    
1838  *     --------------        if(.not.RECOVER_SINGLETS)goto 81
 *     0 = good event  
 *     1 = bad event  
 *     --------------  
       integer iflag  
1839    
1840        integer badseed,badcl  *     ////////////////////////////////////////////////
1841    *     PATCH to recover events with less than 3 couples
1842    *     ////////////////////////////////////////////////    
1843    *     loop over singlet and create "fake" couples
1844    *     (with clx=0 or cly=0)
1845    *    
1846    
1847  *     init variables        if(DEBUG.EQ.1)
1848        ncp_tot=0       $     print*,'>>> Recover singlets '
1849        do ip=1,nplanes       $     ,'(creates fake couples) <<<'
1850           do ico=1,ncouplemax        do icl=1,nclstr1
1851              clx(ip,ico)=0           if(
1852              cly(ip,ico)=0       $        cl_single(icl).eq.1.and.
1853           enddo       $        cl_used(icl).eq.0.and.
1854           ncp_plane(ip)=0       $        .true.)then
          do icl=1,nclstrmax_level2  
             cls(ip,icl)=1  
          enddo  
          ncls(ip)=0  
       enddo  
       do icl=1,nclstrmax_level2  
          cl_single(icl)=1  
          cl_good(icl)=0  
       enddo  
         
 *     start association  
       ncouples=0  
       do icx=1,nclstr1          !loop on cluster (X)  
          if(mod(VIEW(icx),2).eq.1)goto 10  
           
1855  *     ----------------------------------------------------  *     ----------------------------------------------------
1856  *     cut on charge (X VIEW)  *     jump masked views
          if(dedx(icx).lt.dedx_x_min)then  
             cl_single(icx)=0  
             goto 10  
          endif  
 *     cut BAD (X VIEW)              
          badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx)))  
          ifirst=INDSTART(icx)  
          if(icx.ne.nclstr1) then  
             ilast=INDSTART(icx+1)-1  
          else  
             ilast=TOTCLLENGTH  
          endif  
          badcl=badseed  
          do igood=-ngoodstr,ngoodstr  
             ibad=1  
             if((INDMAX(icx)+igood).gt.ifirst.and.  
      $           (INDMAX(icx)+igood).lt.ilast.and.  
      $           .true.)then  
                ibad=BAD(VIEW(icx),  
      $              nvk(MAXS(icx)+igood),  
      $              nst(MAXS(icx)+igood))  
             endif  
             badcl=badcl*ibad  
          enddo  
          if(badcl.eq.0)then     !<<<<<<<<<<<<<< BAD cut  
             cl_single(icx)=0    !<<<<<<<<<<<<<< BAD cut  
             goto 10             !<<<<<<<<<<<<<< BAD cut  
          endif                  !<<<<<<<<<<<<<< BAD cut  
1857  *     ----------------------------------------------------  *     ----------------------------------------------------
1858                        if( mask_view(VIEW(icl)).ne.0 ) goto 21
1859           cl_good(icx)=1              if(mod(VIEW(icl),2).eq.0)then !=== X views
          nplx=npl(VIEW(icx))  
          nldx=nld(MAXS(icx),VIEW(icx))  
           
          do icy=1,nclstr1       !loop on cluster (Y)  
             if(mod(VIEW(icy),2).eq.0)goto 20  
               
1860  *     ----------------------------------------------------  *     ----------------------------------------------------
1861  *     cut on charge (Y VIEW)  *     cut on charge (X VIEW)
             if(dedx(icy).lt.dedx_y_min)then  
                cl_single(icy)=0  
                goto 20  
             endif  
 *     cut BAD (Y VIEW)              
             badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy)))  
             ifirst=INDSTART(icy)  
             if(icy.ne.nclstr1) then  
                ilast=INDSTART(icy+1)-1  
             else  
                ilast=TOTCLLENGTH  
             endif  
             badcl=badseed  
             do igood=-ngoodstr,ngoodstr  
                ibad=1  
                if((INDMAX(icy)+igood).gt.ifirst.and.  
      $              (INDMAX(icy)+igood).lt.ilast.and.  
      $              .true.)  
      $              ibad=BAD(VIEW(icy),  
      $              nvk(MAXS(icy)+igood),  
      $              nst(MAXS(icy)+igood))  
                badcl=badcl*ibad  
             enddo  
             if(badcl.eq.0)then  !<<<<<<<<<<<<<< BAD cut  
                cl_single(icy)=0 !<<<<<<<<<<<<<< BAD cut  
                goto 20          !<<<<<<<<<<<<<< BAD cut  
             endif               !<<<<<<<<<<<<<< BAD cut  
1862  *     ----------------------------------------------------  *     ----------------------------------------------------
1863                               if(sgnl(icl).lt.dedx_x_min) goto 21
               
             cl_good(icy)=1                    
             nply=npl(VIEW(icy))  
             nldy=nld(MAXS(icy),VIEW(icy))  
               
 *     ----------------------------------------------  
 *     CONDITION TO FORM A COUPLE  
 *     ----------------------------------------------  
 *     geometrical consistency (same plane and ladder)  
             if(nply.eq.nplx.and.nldy.eq.nldx)then  
 *     charge correlation  
 *     ===========================================================  
 *     this version of the subroutine is used for the calibration  
 *     thus charge-correlation selection is obviously removed  
 *     ===========================================================  
 c$$$               ddd=(dedx(icy)  
 c$$$     $              -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))  
 c$$$               ddd=ddd/sqrt(kch(nplx,nldx)**2+1)  
 c$$$               cut=chcut*sch(nplx,nldx)  
 c$$$               if(abs(ddd).gt.cut)goto 20 !charge not consistent  
 *     ===========================================================  
1864                                
1865                   nplx=npl(VIEW(icl))
1866    *     ------------------> (FAKE) COUPLE <-----------------
1867                   ncp_plane(nplx) = ncp_plane(nplx) + 1
1868                   clx(nplx,ncp_plane(nplx))=icl
1869                   cly(nplx,ncp_plane(nplx))=0
1870    c$$$  cl_single(icl)=0! I leave the cluster tagged as singlet!!!
1871    *     ----------------------------------------------------
1872                                
1873  *     ------------------> COUPLE <------------------              else                !=== Y views
1874  *     check to do not overflow vector dimentions  *     ----------------------------------------------------
1875  c$$$               if(ncp_plane(nplx).gt.ncouplemax)then  *     cut on charge (Y VIEW)
1876  c$$$                  if(DEBUG)print*,  *     ----------------------------------------------------
1877  c$$$     $                    ' ** warning ** number of identified'//                 if(sgnl(icl).lt.dedx_y_min) goto 21
 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  
1878                                
1879                 if(ncp_plane(nplx).eq.ncouplemax)then                 nply=npl(VIEW(icl))
1880                    if(verbose)print*,  *     ------------------> (FAKE) COUPLE <-----------------
1881       $                 '** warning ** number of identified '//                 ncp_plane(nply) = ncp_plane(nply) + 1
1882       $                 'couples on plane ',nplx,                 clx(nply,ncp_plane(nply))=0
1883       $                 'exceeds vector dimention '                 cly(nply,ncp_plane(nply))=icl
1884       $                 ,'( ',ncouplemax,' )'  c$$$  cl_single(icl)=0! I leave the cluster tagged as singlet!!!
1885  c     good2=.false.  *     ----------------------------------------------------
 c     goto 880   !fill ntp and go to next event                      
                   iflag=1  
                   return  
                endif  
1886                                
1887                 ncp_plane(nplx) = ncp_plane(nplx) + 1              endif
1888                 clx(nplx,ncp_plane(nplx))=icx           endif                  !end "single" condition
1889                 cly(nply,ncp_plane(nplx))=icy   21      continue
1890                 cl_single(icx)=0        enddo                     !end loop over clusters
                cl_single(icy)=0  
             endif                                
 *     ----------------------------------------------  
1891    
1892   20         continue        if(DEBUG.EQ.1)
1893           enddo                  !end loop on clusters(Y)       $     print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)
1894            
1895   10      continue  
1896        enddo                     !end loop on clusters(X)   81   continue
         
         
       do icl=1,nclstr1  
          if(cl_single(icl).eq.1)then  
             ip=npl(VIEW(icl))  
             ncls(ip)=ncls(ip)+1  
             cls(ip,ncls(ip))=icl  
          endif  
       enddo  
         
         
       if(DEBUG)then  
          print*,'clusters  ',nclstr1  
          print*,'good    ',(cl_good(i),i=1,nclstr1)  
          print*,'singles ',(cl_single(i),i=1,nclstr1)  
          print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)  
       endif  
1897                
1898        do ip=1,6        ncp_tot=0
1899           ncp_tot=ncp_tot+ncp_plane(ip)        do ip=1,NPLANES
1900             ncp_tot = ncp_tot + ncp_plane(ip)
1901        enddo        enddo
1902  c     if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)        if(DEBUG.EQ.1)
1903               $     print*,'n.couple tot:      ',ncp_tot
       if(ncp_tot.gt.ncp_max)then  
          if(verbose)print*,  
      $           '** warning ** number of identified '//  
      $           'couples exceeds upper limit for Hough tr. '  
      $           ,'( ',ncp_max,' )'              
 c            good2=.false.  
 c     goto 880       !fill ntp and go to next event  
          iflag=1  
          return  
       endif  
1904                
1905        return        return
1906        end        end
   
1907                
1908  ***************************************************  ***************************************************
1909  *                                                 *  *                                                 *
# Line 2143  c     goto 880       !fill ntp and go to Line 1915  c     goto 880       !fill ntp and go to
1915  **************************************************  **************************************************
1916    
1917        subroutine cp_to_doubtrip(iflag)        subroutine cp_to_doubtrip(iflag)
 c*****************************************************  
 c     02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
1918    
1919        include 'commontracker.f'        include 'commontracker.f'
1920          include 'level1.f'
1921        include 'common_momanhough.f'        include 'common_momanhough.f'
       include 'momanhough_init.f'  
1922        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
1923        include 'common_mini_2.f'        include 'common_mini_2.f'
1924        include 'calib.f'        include 'calib.f'
       include 'level1.f'  
1925    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
1926    
1927  *     output flag  *     output flag
1928  *     --------------  *     --------------
# Line 2188  c      double precision xm3,ym3,zm3 Line 1954  c      double precision xm3,ym3,zm3
1954        real xc,zc,radius        real xc,zc,radius
1955  *     -----------------------------  *     -----------------------------
1956    
1957          if(DEBUG.EQ.1)print*,'cp_to_doubtrip:'
1958    
1959    *     --------------------------------------------
1960    *     put a limit to the maximum number of couples
1961    *     per plane, in order to apply hough transform
1962    *     (couples recovered during track refinement)
1963    *     --------------------------------------------
1964          do ip=1,nplanes
1965             if(ncp_plane(ip).gt.ncouplelimit)then
1966                mask_view(nviewx(ip)) = mask_view(nviewx(ip)) + 2**7
1967                mask_view(nviewy(ip)) = mask_view(nviewy(ip)) + 2**7
1968             endif
1969          enddo
1970    
1971    
1972        ndblt=0                   !number of doublets        ndblt=0                   !number of doublets
1973        ntrpt=0                   !number of triplets        ntrpt=0                   !number of triplets
1974                
1975        do ip1=1,(nplanes-1)      !loop on planes  - COPPIA 1        do ip1=1,(nplanes-1)      !loop on planes  - COPPIA 1
1976           do is1=1,2             !loop on sensors - COPPIA 1  c$$$         print*,'(1) ip ',ip1
1977                c$$$     $        ,mask_view(nviewx(ip1))
1978    c$$$     $        ,mask_view(nviewy(ip1))        
1979             if(  mask_view(nviewx(ip1)).ne.0 .or.
1980         $        mask_view(nviewy(ip1)).ne.0 )goto 10 !skip plane
1981             do is1=1,2             !loop on sensors - COPPIA 1            
1982              do icp1=1,ncp_plane(ip1) !loop on COPPIA 1              do icp1=1,ncp_plane(ip1) !loop on COPPIA 1
1983                 icx1=clx(ip1,icp1)                 icx1=clx(ip1,icp1)
1984                 icy1=cly(ip1,icp1)                 icy1=cly(ip1,icp1)
1985                  
1986    c$$$               print*,'(1) ip ',ip1,' icp ',icp1
1987    
1988  c               call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)  c               call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)
1989                 call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)  c               call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)
1990                 xm1=xPAM                 call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.,0.,0.)
1991                 ym1=yPAM                 xm1=REAL(xPAM) !EM GCC4.7
1992                 zm1=zPAM                                   ym1=REAL(yPAM) !EM GCC4.7
1993  c     print*,'***',is1,xm1,ym1,zm1                 zm1=REAL(zPAM) !EM GCC4.7
1994    
1995                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2                 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2
1996                    do is2=1,2    !loop on sensors -ndblt COPPIA 2  c$$$                  print*,'(2) ip ',ip2
1997                        c$$$     $                 ,mask_view(nviewx(ip2))
1998    c$$$     $                 ,mask_view(nviewy(ip2))
1999                      if(  mask_view(nviewx(ip2)).ne.0 .or.
2000         $                 mask_view(nviewy(ip2)).ne.0 )goto 20 !skip plane
2001                      do is2=1,2    !loop on sensors -ndblt COPPIA 2                    
2002                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2                       do icp2=1,ncp_plane(ip2) !loop on COPPIA 2
2003                          icx2=clx(ip2,icp2)                          icx2=clx(ip2,icp2)
2004                          icy2=cly(ip2,icp2)                          icy2=cly(ip2,icp2)
2005    
2006    c$$$                        print*,'(2) ip ',ip2,' icp ',icp2
2007    
2008  c                        call xyz_PAM  c                        call xyz_PAM
2009  c     $                       (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1)  c     $                       (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1)
2010    c                        call xyz_PAM
2011    c     $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1)
2012                          call xyz_PAM                          call xyz_PAM
2013       $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1)       $                       (icx2,icy2,is2,PFAdef,PFAdef,0.,0.,0.,0.)
2014                          xm2=xPAM                          xm2=REAL(xPAM) !EM GCC4.7
2015                          ym2=yPAM                          ym2=REAL(yPAM) !EM GCC4.7
2016                          zm2=zPAM                          zm2=REAL(zPAM) !EM GCC4.7
2017                                                    
2018    *                       ---------------------------------------------------
2019    *                       both couples must have a y-cluster
2020    *                       (condition necessary when in RECOVER_SINGLETS mode)
2021    *                       ---------------------------------------------------
2022                            if(icy1.eq.0.or.icy2.eq.0)goto 111
2023    
2024                            if(cl_used(icy1).ne.0)goto 111
2025                            if(cl_used(icy2).ne.0)goto 111
2026    
2027                            
2028  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2029  *     track parameters on Y VIEW  *     track parameters on Y VIEW
2030  *     (2 couples needed)  *     (2 couples needed)
2031  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2032                          if(ndblt.eq.ndblt_max)then                          if(ndblt.eq.ndblt_max)then
2033                             if(verbose)print*,                             if(verbose.eq.1)print*,
2034       $                          '** warning ** number of identified '//       $                          '** warning ** number of identified '//
2035       $                          'doublets exceeds vector dimention '       $                          'doublets exceeds vector dimention '
2036       $                          ,'( ',ndblt_max,' )'       $                          ,'( ',ndblt_max,' )'
2037  c                           good2=.false.  c     good2=.false.
2038  c                           goto 880 !fill ntp and go to next event  c     goto 880 !fill ntp and go to next event
2039                               do iv=1,12
2040    c     mask_view(iv) = 3
2041                                  mask_view(iv) = mask_view(iv)+ 2**2
2042                               enddo
2043                             iflag=1                             iflag=1
2044                             return                             return
2045                          endif                          endif
2046                            
2047                            
2048    ccc                        print*,'<doublet> ',icp1,icp2
2049    
2050                          ndblt = ndblt + 1                          ndblt = ndblt + 1
2051  *     store doublet info  *     store doublet info
2052                          cpyz1(ndblt)=id_cp(ip1,icp1,is1)                          cpyz1(ndblt)=id_cp(ip1,icp1,is1)
# Line 2240  c                           goto 880 !fi Line 2054  c                           goto 880 !fi
2054  *     tg(th_yz)  *     tg(th_yz)
2055                          alfayz2(ndblt)=(ym1-ym2)/(zm1-zm2)                          alfayz2(ndblt)=(ym1-ym2)/(zm1-zm2)
2056  *     y0 (cm)  *     y0 (cm)
2057                          alfayz1(ndblt)=alfayz2(ndblt)*(zini-zm1)+ym1                      alfayz1(ndblt)=alfayz2(ndblt)*(REAL(zini)-zm1)+ym1! EM GCC4.7 zm1, ym1 and alfayz1/2 are REAL
2058                                                      
2059  ****  -----------------------------------------------****  ****  -----------------------------------------------****
2060  ****  reject non phisical couples                    ****  ****  reject non phisical couples                    ****
2061  ****  -----------------------------------------------****  ****  -----------------------------------------------****
2062                            if(SECOND_SEARCH)goto 111
2063                          if(                          if(
2064       $                       abs(alfayz2(ndblt)).gt.alfyz2_max       $                       abs(alfayz2(ndblt)).gt.alfyz2_max
2065       $                       .or.       $                       .or.
2066       $                       abs(alfayz1(ndblt)).gt.alfyz1_max       $                       abs(alfayz1(ndblt)).gt.alfyz1_max
2067       $                       )ndblt = ndblt-1       $                       )ndblt = ndblt-1
2068                                                    
2069  c$$$      if(iev.eq.33)then                          
2070  c$$$      print*,'********* ',ndblt,' -- ',icp1,icp2,is1,is2   111                    continue
 c$$$     $        ,' || ',icx1,icy1,icx2,icy2  
 c$$$     $        ,' || ',xm1,ym1,xm2,ym2  
 c$$$     $        ,' || ',alfayz2(ndblt),alfayz1(ndblt)  
 c$$$      endif  
 c$$$  
2071  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2072  *     track parameters on Y VIEW - end  *     track parameters on Y VIEW - end
2073  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2074    
2075    
2076                          if(ip2.eq.nplanes)goto 30 !no possible combination with 3 couples                          if(icx1.ne.0)then
2077                               if(cl_used(icx1).ne.0)goto 31
2078                            endif
2079                            if(icx2.ne.0)then
2080                               if(cl_used(icx2).ne.0)goto 31
2081                            endif
2082    
2083                            if(ip2.eq.nplanes)goto 31 !no possible combination with 3 couples
2084    
2085                          do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3                          do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3
2086    c$$$                           print*,'(3) ip ',ip3
2087    c$$$     $                          ,mask_view(nviewx(ip3))
2088    c$$$     $                          ,mask_view(nviewy(ip3))                          
2089                               if(  mask_view(nviewx(ip3)).ne.0 .or.
2090         $                          mask_view(nviewy(ip3)).ne.0 )goto 30 !skip plane
2091                             do is3=1,2 !loop on sensors - COPPIA 3                             do is3=1,2 !loop on sensors - COPPIA 3
2092                                                                
2093                                do icp3=1,ncp_plane(ip3) !loop on COPPIA 3                                do icp3=1,ncp_plane(ip3) !loop on COPPIA 3
2094                                   icx3=clx(ip3,icp3)                                   icx3=clx(ip3,icp3)
2095                                   icy3=cly(ip3,icp3)                                   icy3=cly(ip3,icp3)
2096    
2097    c$$$                                 print*,'(3) ip ',ip3,' icp ',icp3
2098    
2099    *     ---------------------------------------------------
2100    *     all three couples must have a x-cluster
2101    *     (condition necessary when in RECOVER_SINGLETS mode)
2102    *     ---------------------------------------------------
2103                                     if(
2104         $                                icx1.eq.0.or.
2105         $                                icx2.eq.0.or.
2106         $                                icx3.eq.0.or.
2107         $                                .false.)goto 29
2108                                    
2109                                     if(cl_used(icx1).ne.0)goto 29
2110                                     if(cl_used(icx2).ne.0)goto 29
2111                                     if(cl_used(icx3).ne.0)goto 29
2112    
2113  c                                 call xyz_PAM  c                                 call xyz_PAM
2114  c     $                               (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1)  c     $                               (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1)
2115    c                                 call xyz_PAM
2116    c     $                               (icx3,icy3,is3,PFAdef,PFAdef,0.,0.) !(1)
2117                                   call xyz_PAM                                   call xyz_PAM
2118       $                               (icx3,icy3,is3,PFAdef,PFAdef,0.,0.) !(1)       $                                (icx3,icy3,is3,PFAdef,PFAdef
2119                                   xm3=xPAM       $                                ,0.,0.,0.,0.)
2120                                   ym3=yPAM                                   xm3=REAL(xPAM) !EM GCC4.7
2121                                   zm3=zPAM                                   ym3=REAL(yPAM) !EM GCC4.7
2122                                     zm3=REAL(zPAM) !EM GCC4.7
2123    
2124    
2125  *     find the circle passing through the three points  *     find the circle passing through the three points
2126                                     iflag_t = DEBUG
2127                                   call tricircle(3,xp,zp,angp,resp,chi                                   call tricircle(3,xp,zp,angp,resp,chi
2128       $                                ,xc,zc,radius,iflag)       $                                ,xc,zc,radius,iflag_t)
 c     print*,xc,zc,radius  
2129  *     the circle must intersect the reference plane  *     the circle must intersect the reference plane
2130                                   if(  cc                                 if(iflag.ne.0)goto 29
2131  c     $                                 (xc.le.-1.*xclimit.or.                                   if(iflag_t.ne.0)then
2132  c     $                                 xc.ge.xclimit).and.  *     if tricircle fails, evaluate a straight line
2133       $                                radius**2.ge.(ZINI-zc)**2.and.                                      if(DEBUG.eq.1)
2134       $                                iflag.eq.0.and.       $                                   print*,'TRICIRCLE failure'
2135       $                                .true.)then       $                                   ,' >>> straight line'
2136                                                                        radius=0.
2137                                        xc=0.
2138                                        yc=0.
2139                                        
2140                                        SZZ=0.                  
2141                                        SZX=0.
2142                                        SSX=0.
2143                                        SZ=0.
2144                                        S1=0.
2145                                        X0=0.
2146                                        Ax=0.
2147                                        BX=0.
2148                                        DO I=1,3
2149                                           XX = XP(I)
2150                                           SZZ=SZZ+ZP(I)*ZP(I)
2151                                           SZX=SZX+ZP(I)*XX
2152                                           SSX=SSX+XX
2153                                           SZ=SZ+ZP(I)
2154                                           S1=S1+1.
2155                                        ENDDO
2156                                        DET=SZZ*S1-SZ*SZ
2157                                        AX=(SZX*S1-SZ*SSX)/DET
2158                                        BX=(SZZ*SSX-SZX*SZ)/DET
2159                                        X0  = AX*REAL(ZINI)+BX ! EM GCC4.7
2160                                        
2161                                     endif
2162    
2163                                     if(  .not.SECOND_SEARCH.and.
2164         $                                radius**2.lt.(ZINI-zc)**2)goto 29
2165                                                                      
2166  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2167  *     track parameters on X VIEW  *     track parameters on X VIEW
2168  *     (3 couples needed)  *     (3 couples needed)
2169  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2170                                   if(ntrpt.eq.ntrpt_max)then                                   if(ntrpt.eq.ntrpt_max)then
2171                                      if(verbose)print*,                                      if(verbose.eq.1)print*,
2172       $                     '** warning ** number of identified '//       $                                   '** warning **'//
2173       $                     'triplets exceeds vector dimention '       $                                   ' number of identified '//
2174       $                    ,'( ',ntrpt_max,' )'       $                                   'triplets exceeds'//
2175  c                                    good2=.false.       $                                   ' vector dimention '
2176  c                                    goto 880 !fill ntp and go to next event       $                                   ,'( ',ntrpt_max,' )'
2177    c     good2=.false.
2178    c     goto 880 !fill ntp and go to next event
2179                                        do iv=1,nviews
2180    c     mask_view(iv) = 4
2181                                           mask_view(iv) =
2182         $                                      mask_view(iv)+ 2**3
2183                                        enddo
2184                                      iflag=1                                      iflag=1
2185                                      return                                      return
2186                                   endif                                   endif
2187    
2188    ccc                                 print*,'<triplet> ',icp1,icp2,icp3
2189                                    
2190                                   ntrpt = ntrpt +1                                   ntrpt = ntrpt +1
2191  *     store triplet info  *     store triplet info
2192                                   cpxz1(ntrpt)=id_cp(ip1,icp1,is1)                                   cpxz1(ntrpt)=id_cp(ip1,icp1,is1)
2193                                   cpxz2(ntrpt)=id_cp(ip2,icp2,is2)                                   cpxz2(ntrpt)=id_cp(ip2,icp2,is2)
2194                                   cpxz3(ntrpt)=id_cp(ip3,icp3,is3)                                   cpxz3(ntrpt)=id_cp(ip3,icp3,is3)
2195                                                                    
2196                                   if(xc.lt.0)then                                   if(radius.ne.0.and.xc.lt.0)then
2197  *************POSITIVE DEFLECTION  *************POSITIVE DEFLECTION
2198                alfaxz1(ntrpt) = xc+sqrt(radius**2-(ZINI-zc)**2)             alfaxz1(ntrpt) = xc+sqrt(radius**2-(REAL(ZINI)-zc)**2) !EM GCC4.7
2199                alfaxz2(ntrpt) = (ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2)             alfaxz2(ntrpt) = (REAL(ZINI)-zc)/
2200                alfaxz3(ntrpt) = 1/radius       $          sqrt(radius**2-(REAL(ZINI)-zc)**2) !EM GCC4.7
2201                                   else             alfaxz3(ntrpt) = 1/radius
2202                                    else if(radius.ne.0.and.xc.ge.0)then
2203  *************NEGATIVE DEFLECTION  *************NEGATIVE DEFLECTION
2204                alfaxz1(ntrpt) = xc-sqrt(radius**2-(ZINI-zc)**2)             alfaxz1(ntrpt) = xc-sqrt(radius**2-(REAL(ZINI)-zc)**2)
2205                alfaxz2(ntrpt) = -(ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2)             alfaxz2(ntrpt) = -(REAL(ZINI)-zc)/
2206                alfaxz3(ntrpt) = -1/radius       $          sqrt(radius**2-(REAL(ZINI)-zc)**2) !EM GCC4.7
2207                                   endif             alfaxz3(ntrpt) = -1/radius
2208                                                                    else if(radius.eq.0)then
2209    *************straight fit
2210                 alfaxz1(ntrpt) = X0
2211                 alfaxz2(ntrpt) = AX
2212                 alfaxz3(ntrpt) = 0.
2213                                    endif
2214    
2215    c$$$                                print*,'alfaxz1 ', alfaxz1(ntrpt)
2216    c$$$                                print*,'alfaxz2 ', alfaxz2(ntrpt)
2217    c$$$                                print*,'alfaxz3 ', alfaxz3(ntrpt)
2218                                        
2219  ****  -----------------------------------------------****  ****  -----------------------------------------------****
2220  ****  reject non phisical triplets                   ****  ****  reject non phisical triplets                   ****
2221  ****  -----------------------------------------------****  ****  -----------------------------------------------****
2222                                   if(                                  if(SECOND_SEARCH)goto 29
2223       $                                abs(alfaxz2(ntrpt)).gt.alfxz2_max                                  if(
2224       $                                .or.       $                               abs(alfaxz2(ntrpt)).gt.
2225       $                                abs(alfaxz1(ntrpt)).gt.alfxz1_max       $                               alfxz2_max
2226       $                                )ntrpt = ntrpt-1       $                               .or.
2227                                         $                               abs(alfaxz1(ntrpt)).gt.
2228                                         $                               alfxz1_max
2229  c     print*,alfaxz1(ntrpt),alfaxz2(ntrpt),alfaxz3(ntrpt)       $                               )ntrpt = ntrpt-1
2230                                                                    
2231  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2232  *     track parameters on X VIEW - end  *     track parameters on X VIEW - end
2233  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -  *     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2234                                endif                                
2235     29                           continue
2236                             enddo !end loop on COPPIA 3                             enddo !end loop on COPPIA 3
2237                          enddo   !end loop on sensors - COPPIA 3                          enddo   !end loop on sensors - COPPIA 3
2238     30                     continue
2239                       enddo      !end loop on planes  - COPPIA 3                       enddo      !end loop on planes  - COPPIA 3
2240   30                  continue  
2241                         31                  continue                    
2242   1                enddo         !end loop on COPPIA 2  c 1                enddo         !end loop on COPPIA 2
2243                     enddo         !end loop on COPPIA 2
2244                 enddo            !end loop on sensors - COPPIA 2                 enddo            !end loop on sensors - COPPIA 2
2245     20            continue
2246              enddo               !end loop on planes  - COPPIA 2              enddo               !end loop on planes  - COPPIA 2
2247                
2248    c 11         continue
2249              continue
2250           enddo                  !end loop on COPPIA1           enddo                  !end loop on COPPIA1
2251        enddo                     !end loop on sensors - COPPIA 1        enddo                     !end loop on sensors - COPPIA 1
2252     10   continue
2253        enddo                     !end loop on planes  - COPPIA 1        enddo                     !end loop on planes  - COPPIA 1
2254                
2255        if(DEBUG)then        if(DEBUG.EQ.1)then
2256           print*,'--- doublets ',ndblt           print*,'--- doublets ',ndblt
2257           print*,'--- triplets ',ntrpt           print*,'--- triplets ',ntrpt
2258        endif        endif
# Line 2374  c     goto 880               !ntp fill Line 2277  c     goto 880               !ntp fill
2277        subroutine doub_to_YZcloud(iflag)        subroutine doub_to_YZcloud(iflag)
2278    
2279        include 'commontracker.f'        include 'commontracker.f'
2280          include 'level1.f'
2281        include 'common_momanhough.f'        include 'common_momanhough.f'
2282        include 'momanhough_init.f'  c      include 'momanhough_init.f'
2283    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
2284    
2285  *     output flag  *     output flag
2286  *     --------------  *     --------------
# Line 2397  c      common/dbg/DEBUG Line 2299  c      common/dbg/DEBUG
2299        integer cp_useds1(ncouplemaxtot) ! sensor 1        integer cp_useds1(ncouplemaxtot) ! sensor 1
2300        integer cp_useds2(ncouplemaxtot) ! sensor 2        integer cp_useds2(ncouplemaxtot) ! sensor 2
2301    
2302          if(DEBUG.EQ.1)print*,'doub_to_YZcloud:'
2303    
2304  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2305  *     classification of DOUBLETS  *     classification of DOUBLETS
# Line 2410  c      common/dbg/DEBUG Line 2313  c      common/dbg/DEBUG
2313        distance=0        distance=0
2314        nclouds_yz=0              !number of clouds        nclouds_yz=0              !number of clouds
2315        npt_tot=0        npt_tot=0
2316          nloop=0                  
2317     90   continue                  
2318        do idb1=1,ndblt           !loop (1) on DOUBLETS        do idb1=1,ndblt           !loop (1) on DOUBLETS
2319           if(db_used(idb1).eq.1)goto 2228 !db already included in a cloud           if(db_used(idb1).eq.1)goto 2228 !db already included in a cloud
2320                            
 c     print*,'--------------'  
 c     print*,'** ',idb1,' **'  
               
2321           do icp=1,ncp_tot           do icp=1,ncp_tot
2322              cp_useds1(icp)=0    !init              cp_useds1(icp)=0    !init
2323              cp_useds2(icp)=0    !init              cp_useds2(icp)=0    !init
# Line 2451  ccccc if(db_used(idbref).eq.1)goto 1188 Line 2353  ccccc if(db_used(idbref).eq.1)goto 1188
2353  *     doublet distance in parameter space  *     doublet distance in parameter space
2354                 distance=                 distance=
2355       $              ((alfayz1(idbref)-alfayz1(idb2))/Dalfayz1)**2       $              ((alfayz1(idbref)-alfayz1(idb2))/Dalfayz1)**2
2356       $              +((alfayz2(idbref)-alfayz2(idb2))/Dalfayz2)**2                     $              +((alfayz2(idbref)-alfayz2(idb2))/Dalfayz2)**2
2357                 distance = sqrt(distance)                 distance = sqrt(distance)
2358                                
 c$$$      if(iev.eq.33)then  
 c$$$      if(distance.lt.100)  
 c$$$     $ print*,'********* ',idb1,idbref,idb2,distance  
 c$$$      if(distance.lt.100)  
 c$$$     $ print*,'********* ',alfayz1(idbref),alfayz1(idb2)  
 c$$$     $                    ,alfayz2(idbref),alfayz2(idb2)  
 c$$$      endif  
2359                 if(distance.lt.cutdistyz)then                 if(distance.lt.cutdistyz)then
2360    
 c     print*,idb1,idb2,distance,' cloud ',nclouds_yz  
2361                    if(cpyz1(idb2).gt.0)cp_useds2(cpyz1(idb2))=1                    if(cpyz1(idb2).gt.0)cp_useds2(cpyz1(idb2))=1
2362                    if(cpyz1(idb2).lt.0)cp_useds1(-cpyz1(idb2))=1                    if(cpyz1(idb2).lt.0)cp_useds1(-cpyz1(idb2))=1
2363                    if(cpyz2(idb2).gt.0)cp_useds2(cpyz2(idb2))=1                    if(cpyz2(idb2).gt.0)cp_useds2(cpyz2(idb2))=1
# Line 2479  c     print*,idb1,idb2,distance,' cloud Line 2373  c     print*,idb1,idb2,distance,' cloud
2373    
2374                    temp1 = temp1 + alfayz1(idb2)                    temp1 = temp1 + alfayz1(idb2)
2375                    temp2 = temp2 + alfayz2(idb2)                    temp2 = temp2 + alfayz2(idb2)
 c     print*,'*   idbref,idb2 ',idbref,idb2  
2376                 endif                               endif              
2377                                
2378   1118          continue   1118          continue
2379              enddo               !end loop (2) on DOUBLETS              enddo               !end loop (2) on DOUBLETS
2380                            
2381   1188       continue  c 1188       continue
2382                continue
2383           enddo                  !end loop on... bo?           enddo                  !end loop on... bo?
2384                    
2385           nptloop=npv           nptloop=npv
# Line 2502  c     print*,'*   idbref,idb2 ',idbref,i Line 2396  c     print*,'*   idbref,idb2 ',idbref,i
2396           enddo           enddo
2397           ncpused=0           ncpused=0
2398           do icp=1,ncp_tot           do icp=1,ncp_tot
2399              if(cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0)then              if(
2400         $           (cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0).and.
2401         $           .true.)then
2402                 ncpused=ncpused+1                 ncpused=ncpused+1
2403                 ip=ip_cp(icp)                 ip=ip_cp(icp)
2404                 hit_plane(ip)=1                 hit_plane(ip)=1
# Line 2512  c     print*,'*   idbref,idb2 ',idbref,i Line 2408  c     print*,'*   idbref,idb2 ',idbref,i
2408           do ip=1,nplanes           do ip=1,nplanes
2409              nplused=nplused+ hit_plane(ip)              nplused=nplused+ hit_plane(ip)
2410           enddo           enddo
2411  c     print*,'>>>> ',ncpused,npt,nplused          
          if(ncpused.lt.ncpyz_min)goto 2228 !next doublet  
          if(npt.lt.nptyz_min)goto 2228 !next doublet  
2412           if(nplused.lt.nplyz_min)goto 2228 !next doublet           if(nplused.lt.nplyz_min)goto 2228 !next doublet
2413                    
2414  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
2415  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2416    
2417           if(nclouds_yz.ge.ncloyz_max)then           if(nclouds_yz.ge.ncloyz_max)then
2418              if(verbose)print*,              if(verbose.eq.1)print*,
2419       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
2420       $           'YZ clouds exceeds vector dimention '       $           'YZ clouds exceeds vector dimention '
2421       $           ,'( ',ncloyz_max,' )'       $           ,'( ',ncloyz_max,' )'
2422  c               good2=.false.  c               good2=.false.
2423  c     goto 880         !fill ntp and go to next event  c     goto 880         !fill ntp and go to next event
2424                do iv=1,nviews
2425    c               mask_view(iv) = 5
2426                   mask_view(iv) = mask_view(iv) + 2**4
2427                enddo
2428              iflag=1              iflag=1
2429              return              return
2430           endif           endif
# Line 2542  c     goto 880         !fill ntp and go Line 2440  c     goto 880         !fill ntp and go
2440  c     ptcloud_yz_nt(nclouds_yz)=npt  c     ptcloud_yz_nt(nclouds_yz)=npt
2441           do ipt=1,npt           do ipt=1,npt
2442              db_cloud(npt_tot+ipt) = db_all(ipt)              db_cloud(npt_tot+ipt) = db_all(ipt)
 c     print*,'>> ',ipt,db_all(ipt)  
2443           enddo             enddo  
2444           npt_tot=npt_tot+npt           npt_tot=npt_tot+npt
2445           if(DEBUG)then           if(DEBUG.EQ.1)then
             print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'  
2446              print*,'>>>> cloud ',nclouds_yz,' --- ',npt,' points'              print*,'>>>> cloud ',nclouds_yz,' --- ',npt,' points'
2447              print*,'- alfayz1 ',alfayz1_av(nclouds_yz)              print*,'- alfayz1  ',alfayz1_av(nclouds_yz)
2448              print*,'- alfayz2 ',alfayz2_av(nclouds_yz)              print*,'- alfayz2  ',alfayz2_av(nclouds_yz)
2449              print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot)              print*,'cp_useds1  ',(cp_useds1(icp),icp=1,ncp_tot)
2450              print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot)              print*,'cp_useds2  ',(cp_useds2(icp),icp=1,ncp_tot)
2451              print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)              print*,'cpcloud_yz '
2452  c$$$            print*,'nt-uple: ptcloud_yz(',nclouds_yz,') = '       $           ,(cpcloud_yz(nclouds_yz,icp),icp=1,ncp_tot)
2453  c$$$     $           ,ptcloud_yz(nclouds_yz)              print*,'hit_plane  ',(hit_plane(ip),ip=1,nplanes)
 c$$$            print*,'nt-uple: db_cloud(...) = '  
 c$$$     $           ,(db_cloud(iii),iii=npt_tot-npt+1,npt_tot)  
2454           endif           endif
2455  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2456  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
# Line 2564  c$$$     $           ,(db_cloud(iii),iii Line 2458  c$$$     $           ,(db_cloud(iii),iii
2458        enddo                     !end loop (1) on DOUBLETS        enddo                     !end loop (1) on DOUBLETS
2459                
2460                
2461        if(DEBUG)then        if(nloop.lt.nstepy)then      
2462           print*,'---------------------- '          cutdistyz = cutdistyz+cutystep
2463            nloop     = nloop+1          
2464            goto 90                
2465          endif                    
2466          
2467          if(DEBUG.EQ.1)then
2468           print*,'Y-Z total clouds ',nclouds_yz           print*,'Y-Z total clouds ',nclouds_yz
          print*,' '  
2469        endif        endif
2470                
2471                
# Line 2590  c$$$     $           ,(db_cloud(iii),iii Line 2488  c$$$     $           ,(db_cloud(iii),iii
2488        subroutine trip_to_XZcloud(iflag)        subroutine trip_to_XZcloud(iflag)
2489    
2490        include 'commontracker.f'        include 'commontracker.f'
2491          include 'level1.f'
2492        include 'common_momanhough.f'        include 'common_momanhough.f'
2493        include 'momanhough_init.f'  c      include 'momanhough_init.f'
2494    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
2495    
2496  *     output flag  *     output flag
2497  *     --------------  *     --------------
# Line 2614  c      common/dbg/DEBUG Line 2511  c      common/dbg/DEBUG
2511        integer cp_useds1(ncouplemaxtot) ! sensor 1        integer cp_useds1(ncouplemaxtot) ! sensor 1
2512        integer cp_useds2(ncouplemaxtot) ! sensor 2        integer cp_useds2(ncouplemaxtot) ! sensor 2
2513    
2514          if(DEBUG.EQ.1)print*,'trip_to_XZcloud:'
2515    
2516  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2517  *     classification of TRIPLETS  *     classification of TRIPLETS
2518  *     according to distance in parameter space  *     according to distance in parameter space
# Line 2625  c      common/dbg/DEBUG Line 2524  c      common/dbg/DEBUG
2524        distance=0        distance=0
2525        nclouds_xz=0              !number of clouds                nclouds_xz=0              !number of clouds        
2526        npt_tot=0                 !total number of selected triplets        npt_tot=0                 !total number of selected triplets
2527          nloop=0                  
2528     91   continue                  
2529        do itr1=1,ntrpt           !loop (1) on TRIPLETS        do itr1=1,ntrpt           !loop (1) on TRIPLETS
2530           if(tr_used(itr1).eq.1)goto 22288 !already included in a cloud           if(tr_used(itr1).eq.1)goto 22288 !already included in a cloud
 c     print*,'--------------'  
 c     print*,'** ',itr1,' **'  
2531                    
2532           do icp=1,ncp_tot           do icp=1,ncp_tot
2533              cp_useds1(icp)=0              cp_useds1(icp)=0
# Line 2662  c         tr_temp(1)=itr1 Line 2561  c         tr_temp(1)=itr1
2561              do itr2=1,ntrpt     !loop (2) on TRIPLETS              do itr2=1,ntrpt     !loop (2) on TRIPLETS
2562                 if(itr2.eq.itr1)goto 11188       !next triplet                 if(itr2.eq.itr1)goto 11188       !next triplet
2563                 if(tr_used(itr2).eq.1)goto 11188 !next triplet                               if(tr_used(itr2).eq.1)goto 11188 !next triplet              
2564    
2565    
2566  *     triplet distance in parameter space  *     triplet distance in parameter space
2567  *     solo i due parametri spaziali per il momemnto  *     solo i due parametri spaziali per il momemnto
2568                 distance=                 distance=
2569       $              ((alfaxz1(itrref)-alfaxz1(itr2))/Dalfaxz1)**2       $              ((alfaxz1(itrref)-alfaxz1(itr2))/Dalfaxz1)**2
2570       $              +((alfaxz2(itrref)-alfaxz2(itr2))/Dalfaxz2)**2                     $              +((alfaxz2(itrref)-alfaxz2(itr2))/Dalfaxz2)**2
2571                 distance = sqrt(distance)                 distance = sqrt(distance)
2572                  
2573                 if(distance.lt.cutdistxz)then  
2574  c     print*,idb1,idb2,distance,' cloud ',nclouds_yz  *     ------------------------------------------------------------------------
2575    *     FORCE INCLUSION OF TRIPLETS COMPOSED BY SAME COUPLES, IGNORING THE IMAGE
2576    *     ------------------------------------------------------------------------
2577    *     (added in august 2007)
2578                   istrimage=0
2579                   if(
2580         $              abs(cpxz1(itrref)).eq.abs(cpxz1(itr2)).and.
2581         $              abs(cpxz2(itrref)).eq.abs(cpxz2(itr2)).and.
2582         $              abs(cpxz3(itrref)).eq.abs(cpxz3(itr2)).and.
2583         $              .true.)istrimage=1
2584    
2585                   if(distance.lt.cutdistxz.or.istrimage.eq.1)then
2586                    if(cpxz1(itr2).gt.0)cp_useds2(cpxz1(itr2))=1                    if(cpxz1(itr2).gt.0)cp_useds2(cpxz1(itr2))=1
2587                    if(cpxz1(itr2).lt.0)cp_useds1(-cpxz1(itr2))=1                    if(cpxz1(itr2).lt.0)cp_useds1(-cpxz1(itr2))=1
2588                    if(cpxz2(itr2).gt.0)cp_useds2(cpxz2(itr2))=1                    if(cpxz2(itr2).gt.0)cp_useds2(cpxz2(itr2))=1
# Line 2689  c     print*,idb1,idb2,distance,' cloud Line 2601  c     print*,idb1,idb2,distance,' cloud
2601                    temp1 = temp1 + alfaxz1(itr2)                    temp1 = temp1 + alfaxz1(itr2)
2602                    temp2 = temp2 + alfaxz2(itr2)                    temp2 = temp2 + alfaxz2(itr2)
2603                    temp3 = temp3 + alfaxz3(itr2)                    temp3 = temp3 + alfaxz3(itr2)
 c     print*,'*   itrref,itr2 ',itrref,itr2,distance  
2604                 endif                               endif              
2605                                
2606  11188          continue  11188          continue
2607              enddo               !end loop (2) on TRIPLETS              enddo               !end loop (2) on TRIPLETS
2608                                                
2609  11888       continue  c11888       continue
2610                continue
2611           enddo                  !end loop on... bo?               enddo                  !end loop on... bo?    
2612                    
2613           nptloop=npv           nptloop=npv
# Line 2710  c     print*,'*   itrref,itr2 ',itrref,i Line 2622  c     print*,'*   itrref,itr2 ',itrref,i
2622  *     1bis)  *     1bis)
2623  *     2) it is not already stored  *     2) it is not already stored
2624  *     ------------------------------------------  *     ------------------------------------------
 c     print*,'check cp_used'  
2625           do ip=1,nplanes           do ip=1,nplanes
2626              hit_plane(ip)=0              hit_plane(ip)=0
2627           enddo           enddo
2628           ncpused=0           ncpused=0
2629           do icp=1,ncp_tot           do icp=1,ncp_tot
2630              if(cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0)then              if(
2631         $           (cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0).and.
2632         $           .true.)then
2633                 ncpused=ncpused+1                 ncpused=ncpused+1
2634                 ip=ip_cp(icp)                 ip=ip_cp(icp)
2635                 hit_plane(ip)=1                 hit_plane(ip)=1
# Line 2726  c     print*,'check cp_used' Line 2639  c     print*,'check cp_used'
2639           do ip=1,nplanes           do ip=1,nplanes
2640              nplused=nplused+ hit_plane(ip)              nplused=nplused+ hit_plane(ip)
2641           enddo           enddo
2642           if(ncpused.lt.ncpxz_min)goto 22288 !next triplet           if(nplused.lt.nplxz_min)goto 22288 !next triplet
          if(npt.lt.nptxz_min)goto 22288     !next triplet  
          if(nplused.lt.nplxz_min)goto 22288 !next doublet  
2643                    
2644  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
2645  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2646           if(nclouds_xz.ge.ncloxz_max)then           if(nclouds_xz.ge.ncloxz_max)then
2647              if(verbose)print*,              if(verbose.eq.1)print*,
2648       $           '** warning ** number of identified '//       $           '** warning ** number of identified '//
2649       $           'XZ clouds exceeds vector dimention '       $           'XZ clouds exceeds vector dimention '
2650       $           ,'( ',ncloxz_max,' )'       $           ,'( ',ncloxz_max,' )'
2651  c     good2=.false.  c     good2=.false.
2652  c     goto 880         !fill ntp and go to next event  c     goto 880         !fill ntp and go to next event
2653                do iv=1,nviews
2654    c               mask_view(iv) = 6
2655                   mask_view(iv) =  mask_view(iv) + 2**5
2656                enddo
2657              iflag=1              iflag=1
2658              return              return
2659           endif           endif
# Line 2756  c     goto 880         !fill ntp and go Line 2671  c     goto 880         !fill ntp and go
2671           enddo           enddo
2672           npt_tot=npt_tot+npt           npt_tot=npt_tot+npt
2673                    
2674           if(DEBUG)then           if(DEBUG.EQ.1)then
2675              print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'              print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points'
2676              print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points'                            print*,'- alfaxz1  ',alfaxz1_av(nclouds_xz)
2677              print*,'- alfaxz1 ',alfaxz1_av(nclouds_xz)              print*,'- alfaxz2  ',alfaxz2_av(nclouds_xz)
2678              print*,'- alfaxz2 ',alfaxz2_av(nclouds_xz)              print*,'- alfaxz3  ',alfaxz3_av(nclouds_xz)
2679              print*,'- alfaxz3 ',alfaxz3_av(nclouds_xz)              print*,'cp_useds1  ',(cp_useds1(icp),icp=1,ncp_tot)
2680              print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot)              print*,'cp_useds2  ',(cp_useds2(icp),icp=1,ncp_tot)
2681              print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot)              print*,'cpcloud_xz '
2682         $           ,(cpcloud_xz(nclouds_xz,icp),icp=1,ncp_tot)
2683              print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)              print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)
 c$$$            print*,'nt-uple: ptcloud_xz(',nclouds_xz,') = '  
 c$$$     $           ,ptcloud_xz(nclouds_xz)  
 c$$$            print*,'nt-uple: tr_cloud(...) = '  
 c$$$     $           ,(tr_cloud(iii),iii=npt_tot-npt+1,npt_tot)  
2684           endif           endif
2685  *     >>> NEW CLOUD <<<  *     >>> NEW CLOUD <<<
2686  *     ~~~~~~~~~~~~~~~~~  *     ~~~~~~~~~~~~~~~~~
2687  22288    continue  22288    continue
2688        enddo                     !end loop (1) on DOUBLETS        enddo                     !end loop (1) on DOUBLETS
2689          
2690        if(DEBUG)then         if(nloop.lt.nstepx)then      
2691           print*,'---------------------- '           cutdistxz=cutdistxz+cutxstep
2692             nloop=nloop+1          
2693             goto 91                
2694           endif                    
2695          
2696          if(DEBUG.EQ.1)then
2697           print*,'X-Z total clouds ',nclouds_xz           print*,'X-Z total clouds ',nclouds_xz
          print*,' '  
2698        endif        endif
2699                
2700                
# Line 2796  c$$$     $           ,(tr_cloud(iii),iii Line 2712  c$$$     $           ,(tr_cloud(iii),iii
2712  **************************************************  **************************************************
2713    
2714        subroutine clouds_to_ctrack(iflag)        subroutine clouds_to_ctrack(iflag)
 c*****************************************************  
 c     02/02/2006 modified by Elena Vannuccini --> (1)  
 c*****************************************************  
2715    
2716        include 'commontracker.f'        include 'commontracker.f'
2717          include 'level1.f'
2718        include 'common_momanhough.f'        include 'common_momanhough.f'
2719        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
2720        include 'common_mini_2.f'        include 'common_mini_2.f'
2721        include 'common_mech.f'        include 'common_mech.f'
       include 'momanhough_init.f'  
2722    
2723  c      logical DEBUG  
 c      common/dbg/DEBUG  
2724    
2725  *     output flag  *     output flag
2726  *     --------------  *     --------------
# Line 2824  c      common/dbg/DEBUG Line 2736  c      common/dbg/DEBUG
2736  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2737  *     list of matching couples in the combination  *     list of matching couples in the combination
2738  *     between a XZ and YZ cloud  *     between a XZ and YZ cloud
2739        integer cp_match(nplanes,ncouplemax)        integer cp_match(nplanes,2*ncouplemax)
2740        integer ncp_match(nplanes)        integer ncp_match(nplanes)
2741  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2742        integer hit_plane(nplanes)        integer hit_plane(nplanes)
2743  *     -----------------------------------------------------------  *     -----------------------------------------------------------
2744  *     variables for track fitting  *     variables for track fitting
2745        double precision AL_INI(5)        double precision AL_INI(5)
       double precision tath  
2746  *     -----------------------------------------------------------  *     -----------------------------------------------------------
 c      real fitz(nplanes)        !z coordinates of the planes in cm  
2747    
2748          if(DEBUG.EQ.1)print*,'clouds_to_ctrack:'
2749    
2750    
2751        ntracks=0                 !counter of track candidates        ntracks=0                 !counter of track candidates
2752                
2753        do iyz=1,nclouds_yz       !loop on YZ couds        do iyz=1,nclouds_yz       !loop on YZ clouds
2754           do ixz=1,nclouds_xz    !loop on XZ couds           do ixz=1,nclouds_xz    !loop on XZ clouds
2755                            
2756  *     --------------------------------------------------  *     --------------------------------------------------
2757  *     check of consistency of the clouds  *     check of consistency of the clouds
# Line 2849  c      real fitz(nplanes)        !z coor Line 2760  c      real fitz(nplanes)        !z coor
2760  *     of the two clouds  *     of the two clouds
2761  *     --------------------------------------------------  *     --------------------------------------------------
2762              do ip=1,nplanes              do ip=1,nplanes
2763                 hit_plane(ip)=0                 hit_plane(ip)=0  !n.matching couples (REAL couples, not SINGLETS)
2764                 ncp_match(ip)=0                 ncp_match(ip)=0  !n.matching couples per plane
2765                 do icpp=1,ncouplemax                 do icpp=1,ncouplemax
2766                    cp_match(ip,icpp)=0 !init couple list                    cp_match(ip,icpp)=0 !init couple list
2767                 enddo                 enddo
2768              enddo              enddo
2769              ncp_ok=0              ncp_ok=0            !count n.matching-couples
2770              do icp=1,ncp_tot    !loop on couples              ncpx_ok=0           !count n.matching-couples with x cluster
2771  *     get info on              ncpy_ok=0           !count n.matching-couples with y cluster
2772                 cpintersec(icp)=min(  
2773       $              cpcloud_yz(iyz,icp),  
2774       $              cpcloud_xz(ixz,icp))              do icp=1,ncp_tot    !loop over couples
2775                 if(  
2776       $    (cpcloud_yz(iyz,icp).eq.1.and.cpcloud_xz(ixz,icp).eq.2).or.                 if(.not.RECOVER_SINGLETS)then
2777       $    (cpcloud_yz(iyz,icp).eq.2.and.cpcloud_xz(ixz,icp).eq.1).or.  *     ------------------------------------------------------
2778       $              .false.)cpintersec(icp)=0  *     if NOT in RECOVER_SINGLETS mode, take the intersection
2779    *     between xz yz clouds
2780    *     ------------------------------------------------------
2781                      cpintersec(icp)=min(
2782         $                 cpcloud_yz(iyz,icp),
2783         $                 cpcloud_xz(ixz,icp))
2784    *     cpintersec is >0 if yz and xz clouds contain the same image of couple icp
2785    *     ------------------------------------------------------
2786    *     discard the couple if the sensor is in conflict
2787    *     ------------------------------------------------------
2788                      if(
2789         $       (cpcloud_yz(iyz,icp).eq.1.and.cpcloud_xz(ixz,icp).eq.2).or.
2790         $       (cpcloud_yz(iyz,icp).eq.2.and.cpcloud_xz(ixz,icp).eq.1).or.
2791         $                 .false.)cpintersec(icp)=0
2792                   else
2793    *     ------------------------------------------------------
2794    *     if RECOVER_SINGLETS take the union
2795    *     (otherwise the fake couples formed by singlets would be
2796    *     discarded)    
2797    *     ------------------------------------------------------
2798                      cpintersec(icp)=max(
2799         $                 cpcloud_yz(iyz,icp),
2800         $                 cpcloud_xz(ixz,icp))
2801    c$$$                  if(cpcloud_yz(iyz,icp).gt.0)
2802    c$$$     $                 cpintersec(icp)=cpcloud_yz(iyz,icp)
2803    *     cpintersec is >0 if either yz or xz cloud contains the couple icp
2804                   endif
2805    
2806    c$$$               print*,icp,ip_cp(icp),' -- ',cpintersec(icp)
2807    
2808                 if(cpintersec(icp).ne.0)then                 if(cpintersec(icp).ne.0)then
                   ncp_ok=ncp_ok+1    
2809                                        
2810                    ip=ip_cp(icp)                    ip=ip_cp(icp)
2811                    hit_plane(ip)=1                    hit_plane(ip)=1
2812    c$$$                  if(clx(ip,icp).gt.0.and.cly(ip,icp).gt.0)
2813    c$$$     $                 ncp_ok=ncp_ok+1  
2814    c$$$                  if(clx(ip,icp).gt.0.and.cly(ip,icp).eq.0)
2815    c$$$     $                 ncpx_ok=ncpx_ok+1
2816    c$$$                  if(clx(ip,icp).eq.0.and.cly(ip,icp).gt.0)
2817    c$$$     $                 ncpy_ok=ncpy_ok+1
2818    
2819                      if(  cpcloud_yz(iyz,icp).gt.0.and.
2820         $                 cpcloud_xz(ixz,icp).gt.0)
2821         $                 ncp_ok=ncp_ok+1  
2822                      if(  cpcloud_yz(iyz,icp).gt.0.and.
2823         $                 cpcloud_xz(ixz,icp).eq.0)
2824         $                 ncpy_ok=ncpy_ok+1  
2825                      if(  cpcloud_yz(iyz,icp).eq.0.and.
2826         $                 cpcloud_xz(ixz,icp).gt.0)
2827         $                 ncpx_ok=ncpx_ok+1  
2828    
2829                    if(cpintersec(icp).eq.1)then                    if(cpintersec(icp).eq.1)then
2830  *     1) only the couple image in sensor 1 matches  *     1) only the couple image in sensor 1 matches
2831                       id=-icp                       id=-icp
# Line 2896  c      real fitz(nplanes)        !z coor Line 2852  c      real fitz(nplanes)        !z coor
2852              do ip=1,nplanes              do ip=1,nplanes
2853                 nplused=nplused+ hit_plane(ip)                 nplused=nplused+ hit_plane(ip)
2854              enddo              enddo
2855                                        
2856              if(nplused.lt.nplxz_min)goto 888 !next doublet              if(nplused.lt.3)goto 888 !next combination
2857              if(ncp_ok.lt.ncpok)goto 888 !next cloud  ccc            if(nplused.lt.nplxz_min)goto 888 !next combination
2858                ccc            if(nplused.lt.nplyz_min)goto 888 !next combination
2859              if(DEBUG)then  *     -----------------------------------------------------------
2860                 print*,'Combination ',iyz,ixz  *     if in RECOVER_SINGLET mode, the two clouds must have
2861       $              ,' db ',ptcloud_yz(iyz)  *     at least ONE intersecting real couple
2862       $              ,' tr ',ptcloud_xz(ixz)  *     -----------------------------------------------------------
2863       $              ,'  -----> # matching couples ',ncp_ok              if(ncp_ok.lt.1)goto 888 !next combination
2864              endif  
2865  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'              if(DEBUG.EQ.1)then
2866  c$$$  print*,'Configurazione cluster XZ'                 print*,'////////////////////////////'
2867  c$$$  print*,'1 -- ',(clx(1,i),i=1,ncp_plane(1))                 print*,'Cloud combination (Y,X): ',iyz,ixz
2868  c$$$  print*,'2 -- ',(clx(2,i),i=1,ncp_plane(1))                 print*,' db ',ptcloud_yz(iyz)
2869  c$$$  print*,'3 -- ',(clx(3,i),i=1,ncp_plane(1))                 print*,' tr ',ptcloud_xz(ixz)
2870  c$$$  print*,'4 -- ',(clx(4,i),i=1,ncp_plane(1))                 print*,'  -----> # matching couples ',ncp_ok
2871  c$$$  print*,'5 -- ',(clx(5,i),i=1,ncp_plane(1))                 print*,'  -----> # fake couples (X)',ncpx_ok
2872  c$$$  print*,'6 -- ',(clx(6,i),i=1,ncp_plane(1))                 print*,'  -----> # fake couples (Y)',ncpy_ok
2873  c$$$  print*,'Configurazione cluster YZ'                 do icp=1,ncp_tot
2874  c$$$  print*,'1 -- ',(cly(1,i),i=1,ncp_plane(1))                    print*,'cp ',icp,' >'
2875  c$$$  print*,'2 -- ',(cly(2,i),i=1,ncp_plane(1))       $                 ,' x ',cpcloud_xz(ixz,icp)
2876  c$$$  print*,'3 -- ',(cly(3,i),i=1,ncp_plane(1))       $                 ,' y ',cpcloud_yz(iyz,icp)
2877  c$$$  print*,'4 -- ',(cly(4,i),i=1,ncp_plane(1))       $                 ,' ==> ',cpintersec(icp)
2878  c$$$  print*,'5 -- ',(cly(5,i),i=1,ncp_plane(1))                 enddo
2879  c$$$  print*,'6 -- ',(cly(6,i),i=1,ncp_plane(1))              endif
2880  c$$$  print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'                          
2881                            if(DEBUG.EQ.1)then
 *     -------> INITIAL GUESS <-------  
             AL_INI(1)=dreal(alfaxz1_av(ixz))  
             AL_INI(2)=dreal(alfayz1_av(iyz))  
             AL_INI(4)=datan(dreal(alfayz2_av(iyz))  
      $           /dreal(alfaxz2_av(ixz)))  
             tath=-dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))  
             AL_INI(3)=tath/sqrt(1+tath**2)  
             AL_INI(5)=(1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.  
               
 c     print*,'*******',AL_INI(5)  
             if(AL_INI(5).gt.defmax)goto 888 !next cloud  
               
 c     print*,'alfaxz2, alfayz2 '  
 c     $              ,alfaxz2_av(ixz),alfayz2_av(iyz)  
               
 *     -------> INITIAL GUESS <-------  
 c     print*,'AL_INI ',(al_ini(i),i=1,5)  
               
             if(DEBUG)then  
2882                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))                 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))
2883                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))                 print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))
2884                 print*,'3 >>> ',(cp_match(4,i),i=1,ncp_match(4))                 print*,'3 >>> ',(cp_match(4,i),i=1,ncp_match(4))
# Line 2974  c     print*,'AL_INI ',(al_ini(i),i=1,5) Line 2911  c     print*,'AL_INI ',(al_ini(i),i=1,5)
2911                                hit_plane(6)=icp6                                hit_plane(6)=icp6
2912                                if(ncp_match(6).eq.0)hit_plane(6)=0 !-icp6                                if(ncp_match(6).eq.0)hit_plane(6)=0 !-icp6
2913                                                                
2914                                                                if(DEBUG.eq.1)
2915         $                             print*,'combination: '
2916         $                             ,cp_match(6,icp1)
2917         $                             ,cp_match(5,icp2)
2918         $                             ,cp_match(4,icp3)
2919         $                             ,cp_match(3,icp4)
2920         $                             ,cp_match(2,icp5)
2921         $                             ,cp_match(1,icp6)
2922    
2923    
2924    *                             ---------------------------------------
2925    *                             check if this group of couples has been
2926    *                             already fitted    
2927    *                             ---------------------------------------
2928                                  do ica=1,ntracks
2929                                     isthesame=1
2930                                     do ip=1,NPLANES
2931                                        if(hit_plane(ip).ne.0)then
2932                                           if(  CP_STORE(nplanes-ip+1,ica)
2933         $                                      .ne.
2934         $                                      cp_match(ip,hit_plane(ip)) )
2935         $                                      isthesame=0
2936                                        else
2937                                           if(  CP_STORE(nplanes-ip+1,ica)
2938         $                                      .ne.
2939         $                                      0 )
2940         $                                      isthesame=0
2941                                        endif
2942                                     enddo
2943                                     if(isthesame.eq.1)then
2944                                        if(DEBUG.eq.1)
2945         $                                   print*,'(already fitted)'
2946                                        goto 666 !jump to next combination
2947                                     endif
2948                                  enddo
2949    
2950                                call track_init !init TRACK common                                call track_init !init TRACK common
2951    
2952                                do ip=1,nplanes !loop on planes                                do ip=1,nplanes !loop on planes (bottom to top)
2953                                   if(hit_plane(ip).ne.0)then                                   if(hit_plane(ip).ne.0)then
2954                                      id=cp_match(ip,hit_plane(ip))                                      id=cp_match(ip,hit_plane(ip))
2955                                      is=is_cp(id)                                      is=is_cp(id)
# Line 2991  c     print*,'AL_INI ',(al_ini(i),i=1,5) Line 2963  c     print*,'AL_INI ',(al_ini(i),i=1,5)
2963  *                                   *************************  *                                   *************************
2964  c                                    call xyz_PAM(icx,icy,is,  c                                    call xyz_PAM(icx,icy,is,
2965  c     $                                   'COG2','COG2',0.,0.)  c     $                                   'COG2','COG2',0.,0.)
2966    c                                    call xyz_PAM(icx,icy,is, !(1)
2967    c     $                                   PFAdef,PFAdef,0.,0.) !(1)
2968                                      call xyz_PAM(icx,icy,is, !(1)                                      call xyz_PAM(icx,icy,is, !(1)
2969       $                                   PFAdef,PFAdef,0.,0.) !(1)       $                                   PFAdef,PFAdef,0.,0.,0.,0.)
2970  *                                   *************************  *                                   *************************
2971  *                                   -----------------------------  *                                   -----------------------------
2972                                      xgood(nplanes-ip+1)=1.                                      if(icx.gt.0.and.icy.gt.0)then
2973                                      ygood(nplanes-ip+1)=1.                                         xgood(nplanes-ip+1)=1.
2974                                      xm(nplanes-ip+1)=xPAM                                         ygood(nplanes-ip+1)=1.
2975                                      ym(nplanes-ip+1)=yPAM                                         xm(nplanes-ip+1)=xPAM
2976                                      zm(nplanes-ip+1)=zPAM                                         ym(nplanes-ip+1)=yPAM
2977                                      resx(nplanes-ip+1)=resxPAM                                         zm(nplanes-ip+1)=zPAM
2978                                      resy(nplanes-ip+1)=resyPAM                                         resx(nplanes-ip+1)=resxPAM
2979                                           resy(nplanes-ip+1)=resyPAM
2980                                           if(DEBUG.EQ.1)print*,'(X,Y)'
2981         $                                      ,nplanes-ip+1,xPAM,yPAM
2982                                        else
2983                                           xm_A(nplanes-ip+1) = xPAM_A
2984                                           ym_A(nplanes-ip+1) = yPAM_A
2985                                           xm_B(nplanes-ip+1) = xPAM_B
2986                                           ym_B(nplanes-ip+1) = yPAM_B
2987                                           zm(nplanes-ip+1)
2988         $                                      = (zPAM_A+zPAM_B)/2.
2989                                           resx(nplanes-ip+1) = resxPAM
2990                                           resy(nplanes-ip+1) = resyPAM
2991                                           if(icx.eq.0.and.icy.gt.0)then
2992                                              xgood(nplanes-ip+1)=0.
2993                                              ygood(nplanes-ip+1)=1.
2994                                              resx(nplanes-ip+1) = 1000.
2995                                              if(DEBUG.EQ.1)print*,'(  Y)'
2996         $                                         ,nplanes-ip+1,xPAM,yPAM
2997                                           elseif(icx.gt.0.and.icy.eq.0)then
2998                                              xgood(nplanes-ip+1)=1.
2999                                              ygood(nplanes-ip+1)=0.
3000                                              if(DEBUG.EQ.1)print*,'(X  )'
3001         $                                         ,nplanes-ip+1,xPAM,yPAM
3002                                              resy(nplanes-ip+1) = 1000.
3003                                           else
3004                                              print*,'both icx=0 and icy=0'
3005         $                                         ,' ==> IMPOSSIBLE!!'
3006                                           endif
3007                                        endif
3008  *                                   -----------------------------  *                                   -----------------------------
3009                                   endif                                   endif
3010                                enddo !end loop on planes                                enddo !end loop on planes
3011  *     **********************************************************  *     **********************************************************
3012  *     ************************** FIT *** FIT *** FIT *** FIT ***  *     ************************** FIT *** FIT *** FIT *** FIT ***
3013  *     **********************************************************  *     **********************************************************
3014    cccc  scommentare se si usa al_ini della nuvola
3015    c$$$                              do i=1,5
3016    c$$$                                 AL(i)=AL_INI(i)
3017    c$$$                              enddo
3018                                  call guess()
3019                                do i=1,5                                do i=1,5
3020                                   AL(i)=AL_INI(i)                                   AL_INI(i)=AL(i)
3021                                enddo                                enddo
3022                                ifail=0 !error flag in chi^2 computation                                ifail=0 !error flag in chi^2 computation
3023                                jstep=0 !number of  minimization steps                                jstep=0 !number of  minimization steps
3024                                iprint=0                                iprint=0
3025                                if(DEBUG)iprint=1  c                              if(DEBUG.EQ.1)iprint=1
3026                                  if(DEBUG.EQ.1)iprint=2
3027                                call mini2(jstep,ifail,iprint)                                call mini2(jstep,ifail,iprint)
3028                                if(ifail.ne.0) then                                if(ifail.ne.0) then
3029                                   if(DEBUG)then                                   if(DEBUG.EQ.1)then
3030                                      print *,                                      print *,
3031       $                              '*** MINIMIZATION FAILURE *** '       $                              '*** MINIMIZATION FAILURE *** '
3032       $                              //'(mini2 in clouds_to_ctrack)'       $                              //'(clouds_to_ctrack)'
3033                                        print*,'initial guess: '
3034    
3035                                        print*,'AL_INI(1) = ',AL_INI(1)
3036                                        print*,'AL_INI(2) = ',AL_INI(2)
3037                                        print*,'AL_INI(3) = ',AL_INI(3)
3038                                        print*,'AL_INI(4) = ',AL_INI(4)
3039                                        print*,'AL_INI(5) = ',AL_INI(5)
3040                                   endif                                   endif
3041                                   chi2=-chi2  c                                 chi2=-chi2
3042                                endif                                endif
3043  *     **********************************************************  *     **********************************************************
3044  *     ************************** FIT *** FIT *** FIT *** FIT ***  *     ************************** FIT *** FIT *** FIT *** FIT ***
3045  *     **********************************************************  *     **********************************************************
3046    
3047                                if(chi2.le.0.)goto 666                                              if(chi2.le.0.)goto 666              
3048                                  if(chi2.ge.1.e08)goto 666 !OPTIMIZATION
3049                                  if(chi2.ne.chi2)goto 666  !OPTIMIZATION
3050    
3051  *     --------------------------  *     --------------------------
3052  *     STORE candidate TRACK INFO  *     STORE candidate TRACK INFO
3053  *     --------------------------  *     --------------------------
3054                                if(ntracks.eq.NTRACKSMAX)then                                if(ntracks.eq.NTRACKSMAX)then
3055                                                                    
3056                                   if(verbose)print*,                                   if(verbose.eq.1)print*,
3057       $                 '** warning ** number of candidate tracks '//       $                 '** warning ** number of candidate tracks '//
3058       $                 ' exceeds vector dimension '       $                 ' exceeds vector dimension '
3059       $                ,'( ',NTRACKSMAX,' )'       $                ,'( ',NTRACKSMAX,' )'
3060  c                                 good2=.false.  c                                 good2=.false.
3061  c                                 goto 880 !fill ntp and go to next event                      c                                 goto 880 !fill ntp and go to next event                    
3062                                     do iv=1,nviews
3063    c                                    mask_view(iv) = 7
3064                                        mask_view(iv) = mask_view(iv) + 2**6
3065                                     enddo
3066                                   iflag=1                                   iflag=1
3067                                   return                                   return
3068                                endif                                endif
3069                                                                
3070                                ntracks = ntracks + 1                                ntracks = ntracks + 1
3071                                                                
3072  c$$$                              ndof=0                                                                do ip=1,nplanes !top to bottom
3073                                do ip=1,nplanes  
 c$$$                                 ndof=ndof  
 c$$$     $                                +int(xgood(ip))  
 c$$$     $                                +int(ygood(ip))  
3074                                   XV_STORE(ip,ntracks)=sngl(xv(ip))                                   XV_STORE(ip,ntracks)=sngl(xv(ip))
3075                                   YV_STORE(ip,ntracks)=sngl(yv(ip))                                   YV_STORE(ip,ntracks)=sngl(yv(ip))
3076                                   ZV_STORE(ip,ntracks)=sngl(zv(ip))                                                                       ZV_STORE(ip,ntracks)=sngl(zv(ip))
3077                                   XM_STORE(ip,ntracks)=sngl(xm(ip))                                   XM_STORE(ip,ntracks)=sngl(xm(ip))
3078                                   YM_STORE(ip,ntracks)=sngl(ym(ip))                                   YM_STORE(ip,ntracks)=sngl(ym(ip))
3079                                   ZM_STORE(ip,ntracks)=sngl(zm(ip))                                   ZM_STORE(ip,ntracks)=sngl(zm(ip))
# Line 3067  c$$$     $                               Line 3086  c$$$     $                              
3086                                   AYV_STORE(ip,ntracks)=sngl(ayv(ip))                                   AYV_STORE(ip,ntracks)=sngl(ayv(ip))
3087                                   XGOOD_STORE(ip,ntracks)=sngl(xgood(ip))                                   XGOOD_STORE(ip,ntracks)=sngl(xgood(ip))
3088                                   YGOOD_STORE(ip,ntracks)=sngl(ygood(ip))                                   YGOOD_STORE(ip,ntracks)=sngl(ygood(ip))
3089    *                                NB! hit_plane is defined from bottom to top
3090                                   if(hit_plane(ip).ne.0)then                                   if(hit_plane(ip).ne.0)then
3091                                      CP_STORE(nplanes-ip+1,ntracks)=                                      CP_STORE(nplanes-ip+1,ntracks)=
3092       $                                   cp_match(ip,hit_plane(ip))       $                                   cp_match(ip,hit_plane(ip))
3093                                        SENSOR_STORE(nplanes-ip+1,ntracks)
3094         $                              = is_cp(cp_match(ip,hit_plane(ip)))
3095                                        
3096                                        icl=
3097         $                                   clx(ip,icp_cp(
3098         $                                   cp_match(ip,hit_plane(ip)
3099         $                                   )));
3100                                        if(icl.eq.0)
3101         $                                   icl=
3102         $                                   cly(ip,icp_cp(
3103         $                                   cp_match(ip,hit_plane(ip)
3104         $                                   )));
3105    
3106                                        LADDER_STORE(nplanes-ip+1,ntracks)
3107         $                                   = LADDER(icl);
3108                                   else                                   else
3109                                      CP_STORE(nplanes-ip+1,ntracks)=0                                      CP_STORE(nplanes-ip+1,ntracks)=0
3110                                        SENSOR_STORE(nplanes-ip+1,ntracks)=0
3111                                        LADDER_STORE(nplanes-ip+1,ntracks)=0
3112                                   endif                                   endif
3113                                   CLS_STORE(nplanes-ip+1,ntracks)=0                                   BX_STORE(ip,ntracks)=0!I dont need it now
3114                                     BY_STORE(ip,ntracks)=0!I dont need it now
3115                                     CLS_STORE(ip,ntracks)=0
3116                                   do i=1,5                                   do i=1,5
3117                                      AL_STORE(i,ntracks)=sngl(AL(i))                                      AL_STORE(i,ntracks)=sngl(AL(i))
3118                                   enddo                                   enddo
3119                                enddo                                enddo
3120                                                                
3121  c$$$  *                             Number of Degree Of Freedom                                RCHI2_STORE(ntracks)=REAL(chi2)
 c$$$  ndof=ndof-5                            
 c$$$  *                             reduced chi^2  
 c$$$  rchi2=chi2/dble(ndof)  
                               RCHI2_STORE(ntracks)=chi2  
3122                                                                
3123  *     --------------------------------  *     --------------------------------
3124  *     STORE candidate TRACK INFO - end  *     STORE candidate TRACK INFO - end
# Line 3103  c$$$  rchi2=chi2/dble(ndof) Line 3138  c$$$  rchi2=chi2/dble(ndof)
3138                
3139        if(ntracks.eq.0)then        if(ntracks.eq.0)then
3140           iflag=1           iflag=1
3141           return  cc         return
3142        endif        endif
3143                
3144        if(DEBUG)then        if(DEBUG.EQ.1)then
3145           print*,'****** TRACK CANDIDATES ***********'          print*,'****** TRACK CANDIDATES *****************'
3146           print*,'#         R. chi2        RIG'          print*,'#         R. chi2        RIG         ndof'
3147           do i=1,ntracks          do i=1,ntracks
3148              print*,i,' --- ',rchi2_store(i),' --- '            ndof=0                !(1)
3149       $           ,1./abs(AL_STORE(5,i))            do ii=1,nplanes       !(1)
3150           enddo              ndof=ndof           !(1)
3151           print*,'***********************************'       $           +int(xgood_store(ii,i)) !(1)
3152         $           +int(ygood_store(ii,i)) !(1)
3153              enddo                 !(1)
3154              print*,i,' --- ',rchi2_store(i),' --- '
3155         $         ,1./abs(AL_STORE(5,i)),' --- ',ndof
3156            enddo
3157            print*,'*****************************************'
3158        endif        endif
3159                
3160                
# Line 3132  c$$$  rchi2=chi2/dble(ndof) Line 3173  c$$$  rchi2=chi2/dble(ndof)
3173    
3174        subroutine refine_track(ibest)        subroutine refine_track(ibest)
3175    
 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******************************************************  
3176    
3177        include 'commontracker.f'        include 'commontracker.f'
3178          include 'level1.f'
3179        include 'common_momanhough.f'        include 'common_momanhough.f'
3180        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
3181        include 'common_mini_2.f'        include 'common_mini_2.f'
3182        include 'common_mech.f'        include 'common_mech.f'
       include 'momanhough_init.f'  
       include 'level1.f'  
3183        include 'calib.f'        include 'calib.f'
3184    
 c      logical DEBUG  
 c      common/dbg/DEBUG  
   
3185  *     flag to chose PFA  *     flag to chose PFA
3186        character*10 PFA        character*10 PFA
3187        common/FINALPFA/PFA        common/FINALPFA/PFA
3188    
3189          double precision xmm,rxmm,xmm_A,xmm_B !EM GCC4.7
3190          double precision ymm,rymm,ymm_A,ymm_B !EM GCC4.7
3191          double precision zmm,zmm_A,zmm_B !EM GCC4.7
3192          double precision clincnewc !EM GCC4.7
3193          double precision clincnew !EM GCC4.7
3194    
3195          real k(6)
3196          DATA k/1.099730,0.418900,0.220939,0.220907,0.418771,1.100674/
3197    
3198          real xp,yp,zp
3199          real xyzp(3),bxyz(3)
3200          equivalence (xp,xyzp(1)),(yp,xyzp(2)),(zp,xyzp(3))
3201    
3202          if(DEBUG.EQ.1)print*,'refine_track:'
3203  *     =================================================  *     =================================================
3204  *     new estimate of positions using ETA algorithm  *     new estimate of positions using ETA algorithm
3205  *                          and  *                          and
# Line 3162  c      common/dbg/DEBUG Line 3208  c      common/dbg/DEBUG
3208        call track_init        call track_init
3209        do ip=1,nplanes           !loop on planes        do ip=1,nplanes           !loop on planes
3210    
3211             if(DEBUG.EQ.1)print*,' ........... plane ',ip,' ........... '
3212    
3213             xP=XV_STORE(nplanes-ip+1,ibest)
3214             yP=YV_STORE(nplanes-ip+1,ibest)
3215             zP=ZV_STORE(nplanes-ip+1,ibest)
3216             call gufld(xyzp,bxyz)
3217             BX_STORE(nplanes-ip+1,ibest)=bxyz(1)
3218             BY_STORE(nplanes-ip+1,ibest)=bxyz(2)
3219    c$$$  bxyz(1)=0
3220    c$$$         bxyz(2)=0
3221    c$$$         bxyz(3)=0
3222    *     |||||||||||||||||||||||||||||||||||||||||||||||||
3223  *     -------------------------------------------------  *     -------------------------------------------------
3224  *     If the plane has been already included, it just  *     If the plane has been already included, it just
3225  *     computes again the coordinates of the x-y couple  *     computes again the coordinates of the x-y couple
3226  *     using improved PFAs  *     using improved PFAs
3227  *     -------------------------------------------------  *     -------------------------------------------------
3228           if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and.  *     |||||||||||||||||||||||||||||||||||||||||||||||||
3229    c$$$         if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and.
3230    c$$$     $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then
3231    c$$$            
3232    c$$$            id=CP_STORE(nplanes-ip+1,ibest)
3233    c$$$            
3234    c$$$            is=is_cp(id)
3235    c$$$            icp=icp_cp(id)
3236    c$$$            if(ip_cp(id).ne.ip)
3237    c$$$     $           print*,'OKKIO!!'
3238    c$$$     $           ,'id ',id,is,icp
3239    c$$$     $           ,ip_cp(id),ip
3240    c$$$            icx=clx(ip,icp)
3241    c$$$            icy=cly(ip,icp)
3242    c$$$c            call xyz_PAM(icx,icy,is,
3243    c$$$c     $           PFA,PFA,
3244    c$$$c     $           AXV_STORE(nplanes-ip+1,ibest),
3245    c$$$c     $           AYV_STORE(nplanes-ip+1,ibest))
3246    c$$$            call xyz_PAM(icx,icy,is,
3247    c$$$     $           PFA,PFA,
3248    c$$$     $           AXV_STORE(nplanes-ip+1,ibest),
3249    c$$$     $           AYV_STORE(nplanes-ip+1,ibest),
3250    c$$$     $           bxyz(1),
3251    c$$$     $           bxyz(2)
3252    c$$$     $           )
3253    c$$$
3254    c$$$            xm(nplanes-ip+1) = xPAM
3255    c$$$            ym(nplanes-ip+1) = yPAM
3256    c$$$            zm(nplanes-ip+1) = zPAM
3257    c$$$            xgood(nplanes-ip+1) = 1
3258    c$$$            ygood(nplanes-ip+1) = 1
3259    c$$$            resx(nplanes-ip+1) = resxPAM
3260    c$$$            resy(nplanes-ip+1) = resyPAM
3261    c$$$
3262    c$$$            dedxtrk_x(nplanes-ip+1)=sgnl(icx)/mip(VIEW(icx),LADDER(icx))
3263    c$$$            dedxtrk_y(nplanes-ip+1)=sgnl(icy)/mip(VIEW(icy),LADDER(icy))
3264             if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..or.
3265       $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then       $        YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then
3266                            
3267              id=CP_STORE(nplanes-ip+1,ibest)              id=CP_STORE(nplanes-ip+1,ibest)
# Line 3180  c      common/dbg/DEBUG Line 3274  c      common/dbg/DEBUG
3274       $           ,ip_cp(id),ip       $           ,ip_cp(id),ip
3275              icx=clx(ip,icp)              icx=clx(ip,icp)
3276              icy=cly(ip,icp)              icy=cly(ip,icp)
3277    c            call xyz_PAM(icx,icy,is,
3278    c     $           PFA,PFA,
3279    c     $           AXV_STORE(nplanes-ip+1,ibest),
3280    c     $           AYV_STORE(nplanes-ip+1,ibest))
3281              call xyz_PAM(icx,icy,is,              call xyz_PAM(icx,icy,is,
 c     $           'ETA2','ETA2',  
3282       $           PFA,PFA,       $           PFA,PFA,
3283       $           AXV_STORE(nplanes-ip+1,ibest),       $           AXV_STORE(nplanes-ip+1,ibest),
3284       $           AYV_STORE(nplanes-ip+1,ibest))       $           AYV_STORE(nplanes-ip+1,ibest),
3285  c$$$  call xyz_PAM(icx,icy,is,       $           bxyz(1),
3286  c$$$  $              'COG2','COG2',       $           bxyz(2)
3287  c$$$  $              0.,       $           )
3288  c$$$  $              0.)  
3289              xm(nplanes-ip+1) = xPAM              if(icx.gt.0.and.icy.gt.0)then
3290              ym(nplanes-ip+1) = yPAM                 xm(nplanes-ip+1) = xPAM
3291              zm(nplanes-ip+1) = zPAM                 ym(nplanes-ip+1) = yPAM
3292              xgood(nplanes-ip+1) = 1                 zm(nplanes-ip+1) = zPAM
3293              ygood(nplanes-ip+1) = 1                 xm_A(nplanes-ip+1) = 0.
3294              resx(nplanes-ip+1) = resxPAM                 ym_A(nplanes-ip+1) = 0.
3295              resy(nplanes-ip+1) = resyPAM                 xm_B(nplanes-ip+1) = 0.
3296                   ym_B(nplanes-ip+1) = 0.
3297  c            dedxtrk(nplanes-ip+1) = (dedx(icx)+dedx(icy))/2. !(1)                 xgood(nplanes-ip+1) = 1
3298              dedxtrk_x(nplanes-ip+1)=dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)                 ygood(nplanes-ip+1) = 1
3299              dedxtrk_y(nplanes-ip+1)=dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)                 resx(nplanes-ip+1) = resxPAM
3300                   resy(nplanes-ip+1) = resyPAM
3301                   dedxtrk_x(nplanes-ip+1)=
3302         $              sgnl(icx)/mip(VIEW(icx),LADDER(icx))
3303                   dedxtrk_y(nplanes-ip+1)=
3304         $              sgnl(icy)/mip(VIEW(icy),LADDER(icy))
3305                else
3306                   xm(nplanes-ip+1) = 0.
3307                   ym(nplanes-ip+1) = 0.
3308                   zm(nplanes-ip+1) = (zPAM_A+zPAM_B)/2.
3309                   xm_A(nplanes-ip+1) = xPAM_A
3310                   ym_A(nplanes-ip+1) = yPAM_A
3311                   xm_B(nplanes-ip+1) = xPAM_B
3312                   ym_B(nplanes-ip+1) = yPAM_B
3313                   xgood(nplanes-ip+1) = 0
3314                   ygood(nplanes-ip+1) = 0
3315                   resx(nplanes-ip+1) = 1000.!resxPAM
3316                   resy(nplanes-ip+1) = 1000.!resyPAM
3317                   dedxtrk_x(nplanes-ip+1)= 0
3318                   dedxtrk_y(nplanes-ip+1)= 0
3319                   if(icx.gt.0)then
3320                      xgood(nplanes-ip+1) = 1
3321                      resx(nplanes-ip+1) = resxPAM
3322                      dedxtrk_x(nplanes-ip+1)=
3323         $                 sgnl(icx)/mip(VIEW(icx),LADDER(icx))
3324                   elseif(icy.gt.0)then
3325                      ygood(nplanes-ip+1) = 1
3326                      resy(nplanes-ip+1) = resyPAM
3327                      dedxtrk_y(nplanes-ip+1)=
3328         $                 sgnl(icy)/mip(VIEW(icy),LADDER(icy))
3329                   endif
3330                endif
3331                            
3332    *     |||||||||||||||||||||||||||||||||||||||||||||||||
3333  *     -------------------------------------------------  *     -------------------------------------------------
3334  *     If the plane has NOT  been already included,  *     If the plane has NOT  been already included,
3335  *     it tries to include a COUPLE or a single cluster  *     it tries to include a COUPLE or a single cluster
3336  *     -------------------------------------------------  *     -------------------------------------------------
3337    *     |||||||||||||||||||||||||||||||||||||||||||||||||
3338           else                             else                  
3339                                
3340              xgood(nplanes-ip+1)=0              xgood(nplanes-ip+1)=0
3341              ygood(nplanes-ip+1)=0              ygood(nplanes-ip+1)=0
3342    
3343                CP_STORE(nplanes-ip+1,ibest)=0 !re-init
3344                CLS_STORE(nplanes-ip+1,ibest)=0
3345    
3346                                
3347  *     --------------------------------------------------------------  *     --------------------------------------------------------------
3348  *     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)  
3349              call whichsensor(ip,xP,yP,nldt,ist)              call whichsensor(ip,xP,yP,nldt,ist)
3350  *     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
3351              if(nldt.eq.0.or.ist.eq.0)goto 133              if(nldt.eq.0.or.ist.eq.0)goto 133
3352    
3353                SENSOR_STORE(nplanes-ip+1,IBEST)=ist
3354                LADDER_STORE(nplanes-ip+1,IBEST)=nldt
3355  *     --------------------------------------------------------------  *     --------------------------------------------------------------
3356    
3357              if(DEBUG)then              if(DEBUG.EQ.1)then
3358                 print*,                 print*,
3359       $              '------ Plane ',ip,' intersected on LADDER ',nldt       $              '------ Plane ',ip,' intersected on LADDER ',nldt
3360       $              ,' SENSOR ',ist       $              ,' SENSOR ',ist
# Line 3231  c            dedxtrk(nplanes-ip+1) = (de Line 3365  c            dedxtrk(nplanes-ip+1) = (de
3365  *     ===========================================  *     ===========================================
3366  *     STEP 1 >>>>>>>  try to include a new couple  *     STEP 1 >>>>>>>  try to include a new couple
3367  *     ===========================================  *     ===========================================
3368  c            if(DEBUG)print*,'>>>> try to include a new couple'              distmin=100000000.
             distmin=1000000.  
3369              xmm = 0.              xmm = 0.
3370              ymm = 0.              ymm = 0.
3371              zmm = 0.              zmm = 0.
# Line 3245  c            if(DEBUG)print*,'>>>> try t Line 3378  c            if(DEBUG)print*,'>>>> try t
3378              do icp=1,ncp_plane(ip) !loop on couples on plane icp              do icp=1,ncp_plane(ip) !loop on couples on plane icp
3379                 icx=clx(ip,icp)                 icx=clx(ip,icp)
3380                 icy=cly(ip,icp)                 icy=cly(ip,icp)
3381                   if(icx.eq.0.or.icy.eq.0)goto 1188!if fake couple, jump to next
3382                 if(LADDER(icx).ne.nldt.or. !If the ladder number does not match                 if(LADDER(icx).ne.nldt.or. !If the ladder number does not match
3383  c     $              cl_used(icx).eq.1.or. !or the X cluster is already used  c     $              cl_used(icx).eq.1.or. !or the X cluster is already used
3384  c     $              cl_used(icy).eq.1.or. !or the Y cluster is already used  c     $              cl_used(icy).eq.1.or. !or the Y cluster is already used
3385       $              cl_used(icx).ne.0.or. !or the X cluster is already used !(3)       $              cl_used(icx).ne.0.or. !or the X cluster is already used
3386       $              cl_used(icy).ne.0.or. !or the Y cluster is already used !(3)       $              cl_used(icy).ne.0.or. !or the Y cluster is already used
3387       $              .false.)goto 1188 !then jump to next couple.       $              .false.)goto 1188 !then jump to next couple.
3388  *            *          
3389                 call xyz_PAM(icx,icy,ist,                 call xyz_PAM(icx,icy,ist,
3390       $              PFA,PFA,       $              PFA,PFA,
 c     $              'ETA2','ETA2',  
3391       $              AXV_STORE(nplanes-ip+1,ibest),       $              AXV_STORE(nplanes-ip+1,ibest),
3392       $              AYV_STORE(nplanes-ip+1,ibest))       $              AYV_STORE(nplanes-ip+1,ibest),
3393         $              bxyz(1),
3394         $              bxyz(2)
3395         $              )
3396                                
3397                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3398    c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3399                 id=id_cp(ip,icp,ist)                 id=id_cp(ip,icp,ist)
3400                 if(DEBUG)print*,'( couple ',id                 if(DEBUG.EQ.1)
3401       $              ,' ) normalized distance ',distance       $              print*,'( couple ',id
3402         $              ,' ) distance ',distance
3403                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3404                    xmm = xPAM                    xmm = xPAM
3405                    ymm = yPAM                    ymm = yPAM
# Line 3270  c     $              'ETA2','ETA2', Line 3408  c     $              'ETA2','ETA2',
3408                    rymm = resyPAM                    rymm = resyPAM
3409                    distmin = distance                    distmin = distance
3410                    idm = id                                      idm = id                  
3411  c                 dedxmm = (dedx(icx)+dedx(icy))/2. !(1)                    dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
3412                    dedxmmx = dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)                    dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
3413                    dedxmmy = dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)                    clincnewc=10.*dsqrt(rymm**2+rxmm**2
3414         $            +DBLE(RCHI2_STORE(ibest)*k(ip)*(cov(1,1)+cov(2,2))))! EM GCC4.7
3415                 endif                 endif
3416   1188          continue   1188          continue
3417              enddo               !end loop on couples on plane icp              enddo               !end loop on couples on plane icp
3418              if(distmin.le.clinc)then                                if(distmin.le.clincnewc)then    
3419  *              -----------------------------------  *              -----------------------------------
3420                 xm(nplanes-ip+1) = xmm         !<<<                 xm(nplanes-ip+1) = xmm !<<<
3421                 ym(nplanes-ip+1) = ymm         !<<<                 ym(nplanes-ip+1) = ymm !<<<
3422                 zm(nplanes-ip+1) = zmm         !<<<                 zm(nplanes-ip+1) = zmm !<<<
3423                 xgood(nplanes-ip+1) = 1        !<<<                 xgood(nplanes-ip+1) = 1 !<<<
3424                 ygood(nplanes-ip+1) = 1        !<<<                 ygood(nplanes-ip+1) = 1 !<<<
3425                 resx(nplanes-ip+1)=rxmm        !<<<                 resx(nplanes-ip+1)=rxmm !<<<
3426                 resy(nplanes-ip+1)=rymm        !<<<                 resy(nplanes-ip+1)=rymm !<<<
3427  c              dedxtrk(nplanes-ip+1) = dedxmm !<<<  !(1)                 dedxtrk_x(nplanes-ip+1) = dedxmmx !<<<
3428                 dedxtrk_x(nplanes-ip+1) = dedxmmx    !(1)                 dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<
                dedxtrk_y(nplanes-ip+1) = dedxmmy    !(1)  
3429  *              -----------------------------------  *              -----------------------------------
3430                 CP_STORE(nplanes-ip+1,ibest)=idm                       CP_STORE(nplanes-ip+1,ibest)=idm      
3431                 if(DEBUG)print*,'%%%% included couple ',idm                 if(DEBUG.EQ.1)print*,'%%%% included couple ',idm
3432       $              ,' (norm.dist.= ',distmin,', cut ',clinc,' )'       $              ,' (dist.= ',distmin,', cut ',clincnewc,' )'
3433                 goto 133         !next plane                 goto 133         !next plane
3434              endif              endif
3435  *     ================================================  *     ================================================
3436  *     STEP 2 >>>>>>>  try to include a single cluster  *     STEP 2 >>>>>>>  try to include a single cluster
3437  *                     either from a couple or single  *                     either from a couple or single
3438  *     ================================================  *     ================================================
 c            if(DEBUG)print*,'>>>> try to include a new cluster'  
3439              distmin=1000000.              distmin=1000000.
3440              xmm_A = 0.          !---------------------------              xmm_A = 0.          !---------------------------
3441              ymm_A = 0.          ! init variables that              ymm_A = 0.          ! init variables that
# Line 3317  c            if(DEBUG)print*,'>>>> try t Line 3454  c            if(DEBUG)print*,'>>>> try t
3454              do icp=1,ncp_plane(ip) !loop on cluster inside couples              do icp=1,ncp_plane(ip) !loop on cluster inside couples
3455                 icx=clx(ip,icp)                 icx=clx(ip,icp)
3456                 icy=cly(ip,icp)                 icy=cly(ip,icp)
3457                   if(icx.eq.0.or.icy.eq.0)goto 11882!if fake couple, jump to next
3458                 id=id_cp(ip,icp,ist)                 id=id_cp(ip,icp,ist)
3459                 if(LADDER(icx).ne.nldt)goto 11882 !if the ladder number does not match                 if(LADDER(icx).ne.nldt)goto 11882 !if the ladder number does not match
3460  *                                                !jump to the next couple  *                                                !jump to the next couple
# Line 3324  c            if(DEBUG)print*,'>>>> try t Line 3462  c            if(DEBUG)print*,'>>>> try t
3462  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
3463                 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)
3464  *                                              !jump to the Y cluster  *                                              !jump to the Y cluster
3465    c               call xyz_PAM(icx,0,ist,
3466    c     $              PFA,PFA,
3467    c     $              AXV_STORE(nplanes-ip+1,ibest),0.)              
3468                 call xyz_PAM(icx,0,ist,                 call xyz_PAM(icx,0,ist,
 c     $              'ETA2','ETA2',  
3469       $              PFA,PFA,       $              PFA,PFA,
3470       $              AXV_STORE(nplanes-ip+1,ibest),0.)                     $              AXV_STORE(nplanes-ip+1,ibest),0.,
3471         $              bxyz(1),
3472         $              bxyz(2)
3473         $              )              
3474                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3475  c     if(DEBUG)print*,'normalized distance ',distance  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3476                 if(DEBUG)print*,'( cl-X ',icx                 if(DEBUG.EQ.1)
3477       $              ,' in cp ',id,' ) normalized distance ',distance       $              print*,'( cl-X ',icx
3478         $              ,' in cp ',id,' ) distance ',distance
3479                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3480                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3481                    ymm_A = yPAM_A                    ymm_A = yPAM_A
# Line 3343  c     if(DEBUG)print*,'normalized distan Line 3487  c     if(DEBUG)print*,'normalized distan
3487                    rymm = resyPAM                    rymm = resyPAM
3488                    distmin = distance                    distmin = distance
3489                    iclm = icx                    iclm = icx
3490  c                  dedxmm = dedx(icx) !(1)  c                  dedxmm = sgnl(icx) !(1)
3491                    dedxmmx = dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)                    dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
3492                    dedxmmy = 0.        !(1)                    dedxmmy = 0.        !(1)
3493                 endif                                   endif                  
3494  11881          continue  11881          continue
# Line 3352  c                  dedxmm = dedx(icx) !( Line 3496  c                  dedxmm = dedx(icx) !(
3496  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
3497                 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)
3498  *                                              !jump to the next couple  *                                              !jump to the next couple
3499    c               call xyz_PAM(0,icy,ist,
3500    c     $              PFA,PFA,
3501    c     $              0.,AYV_STORE(nplanes-ip+1,ibest))
3502                 call xyz_PAM(0,icy,ist,                 call xyz_PAM(0,icy,ist,
 c     $              'ETA2','ETA2',  
3503       $              PFA,PFA,       $              PFA,PFA,
3504       $              0.,AYV_STORE(nplanes-ip+1,ibest))       $              0.,AYV_STORE(nplanes-ip+1,ibest),
3505         $              bxyz(1),
3506         $              bxyz(2)
3507         $              )
3508                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3509                 if(DEBUG)print*,'( cl-Y ',icy  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3510       $              ,' in cp ',id,' ) normalized distance ',distance                 if(DEBUG.EQ.1)
3511         $              print*,'( cl-Y ',icy
3512         $              ,' in cp ',id,' ) distance ',distance
3513                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3514                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3515                    ymm_A = yPAM_A                    ymm_A = yPAM_A
# Line 3370  c     $              'ETA2','ETA2', Line 3521  c     $              'ETA2','ETA2',
3521                    rymm = resyPAM                    rymm = resyPAM
3522                    distmin = distance                    distmin = distance
3523                    iclm = icy                    iclm = icy
3524  c                 dedxmm = dedx(icy)  !(1)  c                 dedxmm = sgnl(icy)  !(1)
3525                    dedxmmx = 0.        !(1)                    dedxmmx = 0.        !(1)
3526                    dedxmmy = dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)                    dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
3527                 endif                                   endif                  
3528  11882          continue  11882          continue
3529              enddo               !end loop on cluster inside couples              enddo               !end loop on cluster inside couples
3530  *----- single clusters -----------------------------------------------      *----- single clusters -----------------------------------------------  
3531              do ic=1,ncls(ip)    !loop on single clusters              do ic=1,ncls(ip)    !loop on single clusters
3532                 icl=cls(ip,ic)                 icl=cls(ip,ic)
 c               if(cl_used(icl).eq.1.or.     !if the cluster is already used  
3533                 if(cl_used(icl).ne.0.or.     !if the cluster is already used !(3)                 if(cl_used(icl).ne.0.or.     !if the cluster is already used !(3)
3534       $              LADDER(icl).ne.nldt.or. !or the ladder number does not match       $              LADDER(icl).ne.nldt.or. !or the ladder number does not match
3535       $              .false.)goto 18882      !jump to the next singlet       $              .false.)goto 18882      !jump to the next singlet
3536                 if(mod(VIEW(icl),2).eq.0)then!<---- X view                 if(mod(VIEW(icl),2).eq.0)then!<---- X view
3537                    call xyz_PAM(icl,0,ist,                    call xyz_PAM(icl,0,ist,
 c     $                 'ETA2','ETA2',  
3538       $                 PFA,PFA,       $                 PFA,PFA,
3539       $                 AXV_STORE(nplanes-ip+1,ibest),0.)       $                 AXV_STORE(nplanes-ip+1,ibest),0.,
3540         $                 bxyz(1),
3541         $                 bxyz(2)
3542         $                 )
3543                 else                         !<---- Y view                 else                         !<---- Y view
3544                    call xyz_PAM(0,icl,ist,                    call xyz_PAM(0,icl,ist,
 c     $                 'ETA2','ETA2',  
3545       $                 PFA,PFA,       $                 PFA,PFA,
3546       $                 0.,AYV_STORE(nplanes-ip+1,ibest))       $                 0.,AYV_STORE(nplanes-ip+1,ibest),
3547         $                 bxyz(1),
3548         $                 bxyz(2)
3549         $                 )
3550                 endif                 endif
3551    
3552                 distance = distance_to(XP,YP)                 distance = distance_to(XP,YP)
3553                 if(DEBUG)print*,'( cl-s ',icl  c               distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3554       $              ,' ) normalized distance ',distance                 if(DEBUG.EQ.1)
3555         $              print*,'( cl-s ',icl
3556         $              ,' ) distance ',distance
3557                 if(distance.lt.distmin)then                 if(distance.lt.distmin)then
3558                    xmm_A = xPAM_A                    xmm_A = xPAM_A
3559                    ymm_A = yPAM_A                    ymm_A = yPAM_A
# Line 3409  c     $                 'ETA2','ETA2', Line 3565  c     $                 'ETA2','ETA2',
3565                    rymm = resyPAM                    rymm = resyPAM
3566                    distmin = distance                      distmin = distance  
3567                    iclm = icl                    iclm = icl
 c                  dedxmm = dedx(icl)                   !(1)  
3568                    if(mod(VIEW(icl),2).eq.0)then !<---- X view                    if(mod(VIEW(icl),2).eq.0)then !<---- X view
3569                       dedxmmx = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2)                       dedxmmx = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3570                       dedxmmy = 0.                       !(1)                       dedxmmy = 0.                  
3571                    else          !<---- Y view                    else          !<---- Y view
3572                       dedxmmx = 0.                       !(1)                       dedxmmx = 0.                  
3573                       dedxmmy = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2)                       dedxmmy = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3574                    endif                    endif
3575                 endif                                   endif                  
3576  18882          continue  18882          continue
3577              enddo               !end loop on single clusters              enddo               !end loop on single clusters
3578    
3579              if(distmin.le.clinc)then                                if(iclm.ne.0)then
                 
                CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<      
 *              ----------------------------  
3580                 if(mod(VIEW(iclm),2).eq.0)then                 if(mod(VIEW(iclm),2).eq.0)then
3581                    XGOOD(nplanes-ip+1)=1.                    clincnew=
3582                    resx(nplanes-ip+1)=rxmm       $            20.*     !EM GCC4.7
3583                    if(DEBUG)print*,'%%%% included X-cl ',iclm       $            dsqrt(rxmm**2 +
3584       $                 ,' ( norm.dist.= ',distmin,', cut ',clinc,' )'       $             DBLE(RCHI2_STORE(ibest)*k(ip))*cov(1,1))
3585                 else                 else if(mod(VIEW(iclm),2).ne.0)then
3586                    YGOOD(nplanes-ip+1)=1.                    clincnew=
3587                    resy(nplanes-ip+1)=rymm       $            10.* !EM GCC4.7
3588                    if(DEBUG)print*,'%%%% included Y-cl ',iclm       $            dsqrt(rymm**2 +
3589       $                 ,' ( norm.dist.= ',distmin,', cut ',clinc,' )'       $             DBLE(RCHI2_STORE(ibest)*k(ip))*cov(2,2))
3590                   endif
3591    
3592                   if(distmin.le.clincnew)then  
3593                      
3594                      CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<    
3595    *     ----------------------------
3596                      if(mod(VIEW(iclm),2).eq.0)then
3597                         XGOOD(nplanes-ip+1)=1.
3598                         resx(nplanes-ip+1)=rxmm
3599                         if(DEBUG.EQ.1)
3600         $                    print*,'%%%% included X-cl ',iclm
3601         $                    ,'( chi^2, ',RCHI2_STORE(ibest)
3602         $                    ,', dist.= ',distmin
3603         $                    ,', cut ',clincnew,' )'
3604                      else
3605                         YGOOD(nplanes-ip+1)=1.
3606                         resy(nplanes-ip+1)=rymm
3607                         if(DEBUG.EQ.1)
3608         $                    print*,'%%%% included Y-cl ',iclm
3609         $                    ,'( chi^2, ',RCHI2_STORE(ibest)
3610         $                    ,', dist.= ', distmin
3611         $                    ,', cut ',clincnew,' )'
3612                      endif
3613    *     ----------------------------
3614                      xm_A(nplanes-ip+1) = xmm_A
3615                      ym_A(nplanes-ip+1) = ymm_A
3616                      xm_B(nplanes-ip+1) = xmm_B
3617                      ym_B(nplanes-ip+1) = ymm_B
3618                      zm(nplanes-ip+1) = (zmm_A+zmm_B)/2.
3619                      dedxtrk_x(nplanes-ip+1) = dedxmmx !<<<
3620                      dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<
3621    *     ----------------------------
3622                 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)  
 *              ----------------------------  
3623              endif              endif
3624           endif           endif
3625   133     continue   133     continue
# Line 3456  c              dedxtrk(nplanes-ip+1) = d Line 3630  c              dedxtrk(nplanes-ip+1) = d
3630        return        return
3631        end        end
3632    
3633    
3634  ***************************************************  ***************************************************
3635  *                                                 *  *                                                 *
3636  *                                                 *  *                                                 *
# Line 3464  c              dedxtrk(nplanes-ip+1) = d Line 3639  c              dedxtrk(nplanes-ip+1) = d
3639  *                                                 *  *                                                 *
3640  *                                                 *  *                                                 *
3641  **************************************************  **************************************************
 cccccc 12/08/2006 modified by elena ---> (1)  
3642  *  *
       subroutine clean_XYclouds(ibest,iflag)  
   
       include 'commontracker.f'  
       include 'common_momanhough.f'  
       include 'momanhough_init.f'  
       include 'level2.f'        !(1)  
 c      include 'calib.f'  
 c      include 'level1.f'  
   
 c      logical DEBUG  
 c      common/dbg/DEBUG  
   
   
       do ip=1,nplanes           !loop on planes  
   
          id=CP_STORE(nplanes-ip+1,ibest)  
          icl=CLS_STORE(nplanes-ip+1,ibest)  
          if(id.ne.0.or.icl.ne.0)then                
             if(id.ne.0)then  
                iclx=clx(ip,icp_cp(id))  
                icly=cly(ip,icp_cp(id))  
 c               cl_used(iclx)=1  !tag used clusters  
 c               cl_used(icly)=1  !tag used clusters  
                cl_used(iclx)=ntrk  !tag used clusters !(1)  
                cl_used(icly)=ntrk  !tag used clusters !(1)  
             elseif(icl.ne.0)then  
 c               cl_used(icl)=1   !tag used clusters  
                cl_used(icl)=ntrk   !tag used clusters !1)  
             endif  
               
 c               if(DEBUG)then  
 c                  print*,ip,' <<< ',id  
 c               endif  
 *     -----------------------------  
 *     remove the couple from clouds  
 *     remove also vitual couples containing the  
 *     selected clusters  
 *     -----------------------------  
             do icp=1,ncp_plane(ip)  
                if(  
      $              clx(ip,icp).eq.iclx  
      $              .or.  
      $              clx(ip,icp).eq.icl  
      $              .or.  
      $              cly(ip,icp).eq.icly  
      $              .or.  
      $              cly(ip,icp).eq.icl  
      $              )then  
                   id=id_cp(ip,icp,1)  
                   if(DEBUG)then  
                      print*,ip,' <<< cp ',id  
      $                    ,' ( cl-x '  
      $                    ,clx(ip,icp)  
      $                    ,' cl-y '  
      $                    ,cly(ip,icp),' ) --> removed'  
                   endif  
 *     -----------------------------  
 *     remove the couple from clouds  
                   do iyz=1,nclouds_yz  
                      if(cpcloud_yz(iyz,abs(id)).ne.0)then  
                         ptcloud_yz(iyz)=ptcloud_yz(iyz)-1  
                         cpcloud_yz(iyz,abs(id))=0  
                      endif  
                   enddo  
                   do ixz=1,nclouds_xz  
                      if(cpcloud_xz(ixz,abs(id)).ne.0)then  
                         ptcloud_xz(ixz)=ptcloud_xz(ixz)-1  
                         cpcloud_xz(ixz,abs(id))=0  
                      endif  
                   enddo                      
 *     -----------------------------  
                endif  
             enddo  
               
          endif                
       enddo                     !end loop on planes  
         
       return  
       end  
   
   
   
   
 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$$$  
3643    
3644    
3645    
# Line 3652  c$$$ Line 3647  c$$$
3647    
3648        subroutine init_level2        subroutine init_level2
3649    
 c*****************************************************  
 c     07/10/2005 modified by elena vannuccini --> (1)  
 c*****************************************************  
   
3650        include 'commontracker.f'        include 'commontracker.f'
3651          include 'level1.f'
3652        include 'common_momanhough.f'        include 'common_momanhough.f'
3653        include 'level2.f'        include 'level2.f'
       include 'level1.f'  
3654    
3655    *     ---------------------------------
3656    *     variables initialized from level1
3657    *     ---------------------------------
3658        do i=1,nviews        do i=1,nviews
3659           good2(i)=good1(i)           good2(i)=good1(i)
3660             do j=1,nva1_view
3661                vkflag(i,j)=1
3662                if(cnnev(i,j).le.0)then
3663                   vkflag(i,j)=cnnev(i,j)
3664                endif
3665             enddo
3666        enddo        enddo
3667    *     ----------------
3668  c      good2 = 0!.false.  *     level2 variables
3669  c$$$      nev2 = nev1  *     ----------------
   
 c$$$# ifndef TEST2003  
 c$$$c*****************************************************  
 c$$$cccccc 11/9/2005 modified by david fedele  
 c$$$c      pkt_type = pkt_type1  
 c$$$c      pkt_num = pkt_num1  
 c$$$c      obt = obt1  
 c$$$c      which_calib = which_calib1  
 c$$$      swcode = 302  
 c$$$  
 c$$$      which_calib = which_calib1  
 c$$$      pkt_type = pkt_type1  
 c$$$      pkt_num = pkt_num1  
 c$$$      obt = obt1  
 c$$$      cpu_crc = cpu_crc1  
 c$$$      do iv=1,12  
 c$$$         crc(iv)=crc1(iv)  
 c$$$      enddo  
 c$$$# endif  
 c*****************************************************  
   
3670        NTRK = 0        NTRK = 0
3671        do it=1,NTRKMAX!NTRACKSMAX        do it=1,NTRKMAX
3672           IMAGE(IT)=0           IMAGE(IT)=0
3673           CHI2_nt(IT) = -100000.           CHI2_nt(IT) = -100000.
 c         BdL(IT) = 0.  
3674           do ip=1,nplanes           do ip=1,nplanes
3675              XM_nt(IP,IT) = 0              XM_nt(IP,IT) = 0
3676              YM_nt(IP,IT) = 0              YM_nt(IP,IT) = 0
3677              ZM_nt(IP,IT) = 0              ZM_nt(IP,IT) = 0
3678              RESX_nt(IP,IT) = 0              RESX_nt(IP,IT) = 0
3679              RESY_nt(IP,IT) = 0              RESY_nt(IP,IT) = 0
3680                TAILX_nt(IP,IT) = 0
3681                TAILY_nt(IP,IT) = 0
3682                XBAD(IP,IT) = 0
3683                YBAD(IP,IT) = 0
3684              XGOOD_nt(IP,IT) = 0              XGOOD_nt(IP,IT) = 0
3685              YGOOD_nt(IP,IT) = 0              YGOOD_nt(IP,IT) = 0
3686  c*****************************************************              LS(IP,IT) = 0
 cccccc 11/9/2005 modified by david fedele  
3687              DEDX_X(IP,IT) = 0              DEDX_X(IP,IT) = 0
3688              DEDX_Y(IP,IT) = 0              DEDX_Y(IP,IT) = 0
 c******************************************************  
 cccccc 17/8/2006 modified by elena  
3689              CLTRX(IP,IT) = 0              CLTRX(IP,IT) = 0
3690              CLTRY(IP,IT) = 0              CLTRY(IP,IT) = 0
3691                multmaxx(ip,it) = 0
3692                seedx(ip,it)    = 0
3693                xpu(ip,it)      = 0
3694                multmaxy(ip,it) = 0
3695                seedy(ip,it)    = 0
3696                ypu(ip,it)      = 0
3697           enddo           enddo
3698           do ipa=1,5           do ipa=1,5
3699              AL_nt(IPA,IT) = 0              AL_nt(IPA,IT) = 0
# Line 3717  cccccc 17/8/2006 modified by elena Line 3702  cccccc 17/8/2006 modified by elena
3702              enddo                                enddo                  
3703           enddo                             enddo                  
3704        enddo        enddo
         
         
 c*****************************************************  
 cccccc 11/9/2005 modified by david fedele  
3705        nclsx=0        nclsx=0
3706        nclsy=0              nclsy=0      
3707        do ip=1,NSINGMAX        do ip=1,NSINGMAX
3708          planex(ip)=0          planex(ip)=0
 c        xs(ip)=0  
3709          xs(1,ip)=0          xs(1,ip)=0
3710          xs(2,ip)=0          xs(2,ip)=0
3711          sgnlxs(ip)=0          sgnlxs(ip)=0
3712          planey(ip)=0          planey(ip)=0
 c        ys(ip)=0  
3713          ys(1,ip)=0          ys(1,ip)=0
3714          ys(2,ip)=0          ys(2,ip)=0
3715          sgnlys(ip)=0          sgnlys(ip)=0
3716            sxbad(ip)=0
3717            sybad(ip)=0
3718            multmaxsx(ip)=0
3719            multmaxsy(ip)=0
3720        enddo        enddo
 c*******************************************************  
3721        end        end
3722    
3723    
# Line 3750  c*************************************** Line 3732  c***************************************
3732  ************************************************************  ************************************************************
3733    
3734    
3735          subroutine init_hough
3736    
3737          include 'commontracker.f'
3738          include 'level1.f'
3739          include 'common_momanhough.f'
3740          include 'common_hough.f'
3741          include 'level2.f'
3742    
3743          ntrpt_nt=0
3744          ndblt_nt=0
3745          NCLOUDS_XZ_nt=0
3746          NCLOUDS_YZ_nt=0
3747          do idb=1,ndblt_max_nt
3748             db_cloud_nt(idb)=0
3749             alfayz1_nt(idb)=0      
3750             alfayz2_nt(idb)=0      
3751          enddo
3752          do itr=1,ntrpt_max_nt
3753             tr_cloud_nt(itr)=0
3754             alfaxz1_nt(itr)=0      
3755             alfaxz2_nt(itr)=0      
3756             alfaxz3_nt(itr)=0      
3757          enddo
3758          do idb=1,ncloyz_max      
3759            ptcloud_yz_nt(idb)=0    
3760            alfayz1_av_nt(idb)=0    
3761            alfayz2_av_nt(idb)=0    
3762          enddo                    
3763          do itr=1,ncloxz_max      
3764            ptcloud_xz_nt(itr)=0    
3765            alfaxz1_av_nt(itr)=0    
3766            alfaxz2_av_nt(itr)=0    
3767            alfaxz3_av_nt(itr)=0    
3768          enddo                    
3769    
3770          ntrpt=0                  
3771          ndblt=0                  
3772          NCLOUDS_XZ=0              
3773          NCLOUDS_YZ=0              
3774          do idb=1,ndblt_max        
3775            db_cloud(idb)=0        
3776            cpyz1(idb)=0            
3777            cpyz2(idb)=0            
3778            alfayz1(idb)=0          
3779            alfayz2(idb)=0          
3780          enddo                    
3781          do itr=1,ntrpt_max        
3782            tr_cloud(itr)=0        
3783            cpxz1(itr)=0            
3784            cpxz2(itr)=0            
3785            cpxz3(itr)=0            
3786            alfaxz1(itr)=0          
3787            alfaxz2(itr)=0          
3788            alfaxz3(itr)=0          
3789          enddo                    
3790          do idb=1,ncloyz_max      
3791            ptcloud_yz(idb)=0      
3792            alfayz1_av(idb)=0      
3793            alfayz2_av(idb)=0      
3794            do idbb=1,ncouplemaxtot
3795              cpcloud_yz(idb,idbb)=0
3796            enddo                  
3797          enddo                    
3798          do itr=1,ncloxz_max      
3799            ptcloud_xz(itr)=0      
3800            alfaxz1_av(itr)=0      
3801            alfaxz2_av(itr)=0      
3802            alfaxz3_av(itr)=0      
3803            do itrr=1,ncouplemaxtot
3804              cpcloud_xz(itr,itrr)=0
3805            enddo                  
3806          enddo                    
3807          end
3808    ************************************************************
3809    *
3810    *
3811    *
3812    *
3813    *
3814    *
3815    *
3816    ************************************************************
3817    
3818    
3819        subroutine fill_level2_tracks(ntr)        subroutine fill_level2_tracks(ntr)
3820    
3821  *     -------------------------------------------------------  *     -------------------------------------------------------
# Line 3761  c*************************************** Line 3827  c***************************************
3827            
3828        include 'commontracker.f'        include 'commontracker.f'
3829        include 'level1.f'        include 'level1.f'
3830          include 'common_momanhough.f'
3831        include 'level2.f'        include 'level2.f'
3832        include 'common_mini_2.f'        include 'common_mini_2.f'
3833        include 'common_momanhough.f'        include 'calib.f'
3834        real sinth,phi,pig        !(4)  
3835          character*10 PFA
3836          common/FINALPFA/PFA
3837    
3838          real sinth,phi,pig, npig ! EM GCC4.7
3839          integer ssensor,sladder
3840        pig=acos(-1.)        pig=acos(-1.)
3841    
3842  c      good2=1!.true.  
3843    
3844    *     -------------------------------------
3845        chi2_nt(ntr)        = sngl(chi2)        chi2_nt(ntr)        = sngl(chi2)
3846        nstep_nt(ntr)       = nstep        nstep_nt(ntr)       = nstep
3847    *     -------------------------------------
3848        phi   = al(4)             !(4)        phi   = REAL(al(4))
3849        sinth = al(3)             !(4)        sinth = REAL(al(3))
3850        if(sinth.lt.0)then        !(4)        if(sinth.lt.0)then      
3851           sinth = -sinth         !(4)           sinth = -sinth        
3852           phi = phi + pig        !(4)           phi = phi + pig      
3853        endif                     !(4)        endif                    
3854        npig = aint(phi/(2*pig))  !(4)        npig = aint(phi/(2*pig))
3855        phi = phi - npig*2*pig    !(4)        phi = phi - npig*2*pig  
3856        if(phi.lt.0)              !(4)        if(phi.lt.0)            
3857       $     phi = phi + 2*pig    !(4)       $     phi = phi + 2*pig  
3858        al(4) = phi               !(4)        al(4) = phi              
3859        al(3) = sinth             !(4)        al(3) = sinth            
 *****************************************************  
3860        do i=1,5        do i=1,5
3861           al_nt(i,ntr)     = sngl(al(i))           al_nt(i,ntr)     = sngl(al(i))
3862           do j=1,5           do j=1,5
3863              coval(i,j,ntr) = sngl(cov(i,j))              coval(i,j,ntr) = sngl(cov(i,j))
3864           enddo           enddo
 c     print*,al_nt(i,ntr)  
3865        enddo        enddo
3866          *     -------------------------------------      
3867        do ip=1,nplanes           ! loop on planes        do ip=1,nplanes           ! loop on planes
3868           xgood_nt(ip,ntr) = int(xgood(ip))           xgood_nt(ip,ntr) = int(xgood(ip))
3869           ygood_nt(ip,ntr) = int(ygood(ip))           ygood_nt(ip,ntr) = int(ygood(ip))
# Line 3800  c     print*,al_nt(i,ntr) Line 3872  c     print*,al_nt(i,ntr)
3872           zm_nt(ip,ntr)    = sngl(zm(ip))           zm_nt(ip,ntr)    = sngl(zm(ip))
3873           RESX_nt(IP,ntr)  = sngl(resx(ip))           RESX_nt(IP,ntr)  = sngl(resx(ip))
3874           RESY_nt(IP,ntr)  = sngl(resy(ip))           RESY_nt(IP,ntr)  = sngl(resy(ip))
3875             TAILX_nt(IP,ntr) = 0.
3876             TAILY_nt(IP,ntr) = 0.
3877           xv_nt(ip,ntr)    = sngl(xv(ip))           xv_nt(ip,ntr)    = sngl(xv(ip))
3878           yv_nt(ip,ntr)    = sngl(yv(ip))           yv_nt(ip,ntr)    = sngl(yv(ip))
3879           zv_nt(ip,ntr)    = sngl(zv(ip))           zv_nt(ip,ntr)    = sngl(zv(ip))
3880           axv_nt(ip,ntr)   = sngl(axv(ip))           axv_nt(ip,ntr)   = sngl(axv(ip))
3881           ayv_nt(ip,ntr)   = sngl(ayv(ip))           ayv_nt(ip,ntr)   = sngl(ayv(ip))  
3882  c        dedxp(ip,ntr)    = sngl(dedxtrk(ip))   !(1)  
3883           dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)) !(2)           factor = sqrt(
3884           dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)) !(2)         $        tan( acos(-1.) * sngl(axv(ip)) /180. )**2 +
3885         $        tan( acos(-1.) * sngl(ayv(ip)) /180. )**2 +
3886         $        1. )
3887    
3888             dedx_x(ip,ntr)   = sngl(dedxtrk_x(ip)/factor)
3889             dedx_y(ip,ntr)   = sngl(dedxtrk_y(ip)/factor)  
3890        
3891           id  = CP_STORE(ip,IDCAND)  
3892    ccc         print*,ip,'dedx >>> ',dedx_x(ip,ntr),dedx_y(ip,ntr)
3893    
3894             ax   = axv_nt(ip,ntr)
3895             ay   = ayv_nt(ip,ntr)
3896             bfx  = BX_STORE(ip,IDCAND)
3897             bfy  = BY_STORE(ip,IDCAND)
3898    c$$$         if(ip.eq.6) ax = -1. * axv_nt(ip,ntr)
3899    c$$$         if(ip.eq.6) bfy = -1. * BY_STORE(ip,IDCAND)
3900    c$$$         tgtemp   = tan(ax*acos(-1.)/180.) + pmuH_h*bfy*0.00001
3901    c$$$         angx     = 180.*atan(tgtemp)/acos(-1.)
3902    c$$$         tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001        
3903    c$$$         angy    = 180.*atan(tgtemp)/acos(-1.)
3904    
3905             angx = effectiveangle(ax,2*ip,bfy)
3906             angy = effectiveangle(ay,2*ip-1,bfx)
3907            
3908            
3909    
3910             id  = CP_STORE(ip,IDCAND) ! couple id
3911           icl = CLS_STORE(ip,IDCAND)           icl = CLS_STORE(ip,IDCAND)
3912           if(id.ne.0)then           ssensor = -1
3913             sladder = -1
3914             ssensor = SENSOR_STORE(ip,IDCAND)
3915             sladder = LADDER_STORE(ip,IDCAND)
3916             if(ip.eq.6.and.ssensor.ne.0)ssensor = 3 - ssensor !notazione paolo x align
3917             LS(IP,ntr)      = ssensor+10*sladder
3918    
3919    c         if(id.ne.0)then
3920    CCCCCC(10 novembre 2009) PATCH X NUCLEI
3921    C     non ho capito perche', ma durante il ritracciamento dei nuclei
3922    C     (quando una traccia viene trovata ma non e' buona) c'e' qualche variabile
3923    C     che non viene reinizializzata correttamente e i cluster esclusi
3924    C     dal fit risultano ancora inclusi...
3925    C
3926             cltrx(ip,ntr)   = 0
3927             cltry(ip,ntr)   = 0
3928    c$$$         if(
3929    c$$$     $        xgood_nt(ip,ntr).eq.1.and.ygood_nt(ip,ntr).eq.1
3930    c$$$     $        .and.
3931    c$$$     $        id.ne.0)then
3932             if(id.ne.0)then        !patch 30/12/09
3933    
3934    c           >>> is a couple
3935              cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))              cltrx(ip,ntr)   = clx(nplanes-ip+1,icp_cp(id))
3936              cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))              cltry(ip,ntr)   = cly(nplanes-ip+1,icp_cp(id))
3937  c            print*,ip,' ',cltrx(ip,ntr),cltry(ip,ntr)  
3938           elseif(icl.ne.0)then              if(clx(nplanes-ip+1,icp_cp(id)).gt.0)then
3939              if(mod(VIEW(icl),2).eq.0)cltrx(ip,ntr)=icl  
3940              if(mod(VIEW(icl),2).eq.1)cltry(ip,ntr)=icl                 cl_used(cltrx(ip,ntr)) = 1 !tag used clusters          
3941  c            print*,ip,' ',cltrx(ip,ntr),cltry(ip,ntr)  
3942                   xbad(ip,ntr)= nbadstrips(4,clx(nplanes-ip+1,icp_cp(id)))
3943    
3944                   if(nsatstrips(clx(nplanes-ip+1,icp_cp(id))).gt.0)
3945         $              dedx_x(ip,ntr)=-dedx_x(ip,ntr)
3946                  
3947                   multmaxx(ip,ntr) = maxs(cltrx(ip,ntr))
3948         $              +10000*mult(cltrx(ip,ntr))
3949                   seedx(ip,ntr)    = clsignal(indmax(cltrx(ip,ntr)))
3950         $              /clsigma(indmax(cltrx(ip,ntr)))
3951                   call applypfa(PFA,cltrx(ip,ntr),angx,corr,res)
3952                   xpu(ip,ntr)      = corr
3953    
3954                endif
3955                if(cly(nplanes-ip+1,icp_cp(id)).gt.0)then
3956    
3957                   cl_used(cltry(ip,ntr)) = 1 !tag used clusters          
3958    
3959                   ybad(ip,ntr)= nbadstrips(4,cly(nplanes-ip+1,icp_cp(id)))
3960    
3961                   if(nsatstrips(cly(nplanes-ip+1,icp_cp(id))).gt.0)
3962         $              dedx_y(ip,ntr)=-dedx_y(ip,ntr)
3963                  
3964                   multmaxy(ip,ntr) = maxs(cltry(ip,ntr))
3965         $              +10000*mult(cltry(ip,ntr))
3966                   seedy(ip,ntr)    = clsignal(indmax(cltry(ip,ntr)))
3967         $              /clsigma(indmax(cltry(ip,ntr)))
3968                   call applypfa(PFA,cltry(ip,ntr),angy,corr,res)
3969                   ypu(ip,ntr)      = corr
3970                endif
3971    
3972    c$$$         elseif(icl.ne.0)then
3973             endif                  !patch 30/12/09
3974             if(icl.ne.0)then       !patch 30/12/09
3975    
3976                cl_used(icl) = 1    !tag used clusters          
3977    
3978                if(mod(VIEW(icl),2).eq.0)then
3979                   cltrx(ip,ntr)=icl
3980                   xbad(ip,ntr) = nbadstrips(4,icl)
3981    
3982                   if(nsatstrips(icl).gt.0)dedx_x(ip,ntr)=-dedx_x(ip,ntr)
3983    
3984                   multmaxx(ip,ntr) = maxs(cltrx(ip,ntr))
3985         $                         +10000*mult(cltrx(ip,ntr))
3986                   seedx(ip,ntr)    = clsignal(indmax(cltrx(ip,ntr)))
3987         $           /clsigma(indmax(cltrx(ip,ntr)))
3988                   call applypfa(PFA,cltrx(ip,ntr),angx,corr,res)
3989                   xpu(ip,ntr)      = corr
3990    
3991                elseif(mod(VIEW(icl),2).eq.1)then
3992                   cltry(ip,ntr)=icl
3993                   ybad(ip,ntr) = nbadstrips(4,icl)
3994    
3995                   if(nsatstrips(icl).gt.0)dedx_y(ip,ntr)=-dedx_y(ip,ntr)
3996    
3997                   multmaxy(ip,ntr) = maxs(cltry(ip,ntr))
3998         $                         +10000*mult(cltry(ip,ntr))
3999                   seedy(ip,ntr)    = clsignal(indmax(cltry(ip,ntr)))
4000         $           /clsigma(indmax(cltry(ip,ntr)))
4001                   call applypfa(PFA,cltry(ip,ntr),angy,corr,res)
4002                   ypu(ip,ntr)      = corr
4003                  
4004                endif
4005    
4006           endif                     endif          
4007    
4008        enddo        enddo
 c      call CalcBdL(100,xxxx,IFAIL)  
 c      if(ifps(xxxx).eq.1)BdL(ntr) = xxxx  
 c$$$      print*,'xgood(ip,ntr) ',(xgood_nt(ip,ntr),ip=1,6)  
 c$$$      print*,'ygood(ip,ntr) ',(ygood_nt(ip,ntr),ip=1,6)  
 c$$$      print*,'dedx_x(ip,ntr) ',(dedx_x(ip,ntr),ip=1,6)  
 c$$$      print*,'dedx_y(ip,ntr) ',(dedx_y(ip,ntr),ip=1,6)  
4009    
4010          if(DEBUG.eq.1)then
4011             print*,'> STORING TRACK ',ntr
4012             print*,'clusters: '
4013             do ip=1,6
4014                print*,'> ',ip,' -- ',cltrx(ip,ntr),cltry(ip,ntr)
4015             enddo
4016             print*,'dedx: '
4017             do ip=1,6
4018                print*,'> ',ip,' -- ',dedx_x(ip,ntr),dedx_y(ip,ntr)
4019             enddo
4020          endif
4021    
4022        end        end
4023    
4024        subroutine fill_level2_siglets        subroutine fill_level2_siglets
 c*****************************************************  
 c     07/10/2005 created by elena vannuccini  
 c     31/01/2006 modified by elena vannuccini  
 *     to convert adc to mip  --> (2)  
 c*****************************************************  
4025    
4026  *     -------------------------------------------------------  *     -------------------------------------------------------
4027  *     This routine fills the  elements of the variables  *     This routine fills the  elements of the variables
# Line 3846  c*************************************** Line 4030  c***************************************
4030  *     -------------------------------------------------------  *     -------------------------------------------------------
4031    
4032        include 'commontracker.f'        include 'commontracker.f'
       include 'level1.f'  
       include 'level2.f'  
4033        include 'calib.f'        include 'calib.f'
4034          include 'level1.f'
4035        include 'common_momanhough.f'        include 'common_momanhough.f'
4036          include 'level2.f'
4037        include 'common_xyzPAM.f'        include 'common_xyzPAM.f'
4038    
4039  *     count #cluster per plane not associated to any track  *     count #cluster per plane not associated to any track
 c      good2=1!.true.  
4040        nclsx = 0        nclsx = 0
4041        nclsy = 0        nclsy = 0
4042    
4043          do iv = 1,nviews
4044    c         if( mask_view(iv).ne.0 )good2(iv) = 20+mask_view(iv)
4045             good2(iv) = good2(iv) + mask_view(iv)*2**8
4046          enddo
4047    
4048          if(DEBUG.eq.1)then
4049             print*,'> STORING SINGLETS '
4050          endif
4051    
4052        do icl=1,nclstr1        do icl=1,nclstr1
4053    
4054             ip=nplanes-npl(VIEW(icl))+1            
4055            
4056           if(cl_used(icl).eq.0)then !cluster not included in any track           if(cl_used(icl).eq.0)then !cluster not included in any track
4057              ip=nplanes-npl(VIEW(icl))+1              
4058              if(mod(VIEW(icl),2).eq.0)then !=== X views              if(mod(VIEW(icl),2).eq.0)then !=== X views
4059    
4060                 nclsx = nclsx + 1                 nclsx = nclsx + 1
4061                 planex(nclsx) = ip                 planex(nclsx) = ip
4062                 sgnlxs(nclsx) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)                 sgnlxs(nclsx) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
4063                   if(nsatstrips(icl).gt.0)sgnlxs(nclsx)=-sgnlxs(nclsx)
4064                 clsx(nclsx)   = icl                 clsx(nclsx)   = icl
4065                   sxbad(nclsx)  = nbadstrips(1,icl)
4066                   multmaxsx(nclsx) = maxs(icl)+10000*mult(icl)
4067                  
4068    
4069                 do is=1,2                 do is=1,2
4070  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)  c                  call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)
4071                    call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)  c                  call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)
4072                    xs(is,nclsx) = (xPAM_A+xPAM_B)/2                    call xyz_PAM(icl,0,is,PFAdef,'    ',0.,0.,0.,0.)
4073                      xs(is,nclsx) = REAL((xPAM_A+xPAM_B)/2.) ! EM GCC4.7
4074                 enddo                 enddo
 c$$$               print*,'nclsx         ',nclsx  
 c$$$               print*,'planex(nclsx) ',planex(nclsx)  
 c$$$               print*,'sgnlxs(nclsx) ',sgnlxs(nclsx)  
 c$$$               print*,'xs(1,nclsx)   ',xs(1,nclsx)  
 c$$$               print*,'xs(2,nclsx)   ',xs(2,nclsx)  
4075              else                          !=== Y views              else                          !=== Y views
4076                 nclsy = nclsy + 1                 nclsy = nclsy + 1
4077                 planey(nclsy) = ip                 planey(nclsy) = ip
4078                 sgnlys(nclsy) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)                 sgnlys(nclsy) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
4079                   if(nsatstrips(icl).gt.0)sgnlys(nclsy)=-sgnlys(nclsy)
4080                 clsy(nclsy)   = icl                 clsy(nclsy)   = icl
4081                   sybad(nclsy)  = nbadstrips(1,icl)
4082                   multmaxsy(nclsy) = maxs(icl)+10000*mult(icl)
4083    
4084    
4085                 do is=1,2                 do is=1,2
4086  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)  c                  call xyz_PAM(0,icl,is,' ','COG1',0.,0.)
4087                    call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)  c                  call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)
4088                    ys(is,nclsy) = (yPAM_A+yPAM_B)/2                    call xyz_PAM(0,icl,is,'    ',PFAdef,0.,0.,0.,0.)
4089                      ys(is,nclsy) = REAL((yPAM_A+yPAM_B)/2.) ! EM GCC4.7
4090                 enddo                 enddo
 c$$$               print*,'nclsy         ',nclsy  
 c$$$               print*,'planey(nclsy) ',planey(nclsy)  
 c$$$               print*,'sgnlys(nclsy) ',sgnlys(nclsy)  
 c$$$               print*,'ys(1,nclsy)   ',ys(1,nclsy)  
 c$$$               print*,'ys(2,nclsy)   ',ys(2,nclsy)  
4091              endif              endif
4092           endif           endif
 c      print*,icl,cl_used(icl),cl_good(icl),ip,VIEW(icl)!nclsx(ip),nclsy(ip)  
4093    
4094  ***** LO METTO QUI PERCHE` NON SO DOVE METTERLO  ***** LO METTO QUI PERCHE` NON SO DOVE METTERLO
4095           whichtrack(icl) = cl_used(icl)           whichtrack(icl) = cl_used(icl)
4096    *     --------------------------------------------------
4097    *     per non perdere la testa...
4098    *     whichtrack(icl) e` una variabile del common level1
4099    *     che serve solo per sapere quali cluster sono stati
4100    *     associati ad una traccia, e permettere di salvare
4101    *     solo questi nell'albero di uscita
4102    *     --------------------------------------------------
4103                    
4104        enddo        enddo
4105        end        end
4106    
4107    ***************************************************
4108    *                                                 *
4109    *                                                 *
4110    *                                                 *
4111    *                                                 *
4112    *                                                 *
4113    *                                                 *
4114    **************************************************
4115    
4116          subroutine fill_hough
4117    
4118    *     -------------------------------------------------------
4119    *     This routine fills the  variables related to the hough
4120    *     transform, for the debig n-tuple
4121    *     -------------------------------------------------------
4122    
4123          include 'commontracker.f'
4124          include 'level1.f'
4125          include 'common_momanhough.f'
4126          include 'common_hough.f'
4127          include 'level2.f'
4128    
4129          if(.false.
4130         $     .or.ntrpt.gt.ntrpt_max_nt
4131         $     .or.ndblt.gt.ndblt_max_nt
4132         $     .or.NCLOUDS_XZ.gt.ncloxz_max
4133         $     .or.NCLOUDS_yZ.gt.ncloyz_max
4134         $     )then
4135             ntrpt_nt=0
4136             ndblt_nt=0
4137             NCLOUDS_XZ_nt=0
4138             NCLOUDS_YZ_nt=0
4139          else
4140             ndblt_nt=ndblt
4141             ntrpt_nt=ntrpt
4142             if(ndblt.ne.0)then
4143                do id=1,ndblt
4144                   alfayz1_nt(id)=alfayz1(id) !Y0
4145                   alfayz2_nt(id)=alfayz2(id) !tg theta-yz
4146                enddo
4147             endif
4148             if(ndblt.ne.0)then
4149                do it=1,ntrpt
4150                   alfaxz1_nt(it)=alfaxz1(it) !X0
4151                   alfaxz2_nt(it)=alfaxz2(it) !tg theta-xz
4152                   alfaxz3_nt(it)=alfaxz3(it) !1/r
4153                enddo
4154             endif
4155             nclouds_yz_nt=nclouds_yz
4156             nclouds_xz_nt=nclouds_xz
4157             if(nclouds_yz.ne.0)then
4158                nnn=0
4159                do iyz=1,nclouds_yz
4160                   ptcloud_yz_nt(iyz)=ptcloud_yz(iyz)
4161                   alfayz1_av_nt(iyz)=alfayz1_av(iyz)
4162                   alfayz2_av_nt(iyz)=alfayz2_av(iyz)
4163                   nnn=nnn+ptcloud_yz(iyz)
4164                enddo
4165                do ipt=1,min(ndblt_max_nt,nnn)
4166                   db_cloud_nt(ipt)=db_cloud(ipt)
4167                 enddo
4168             endif
4169             if(nclouds_xz.ne.0)then
4170                nnn=0
4171                do ixz=1,nclouds_xz
4172                   ptcloud_xz_nt(ixz)=ptcloud_xz(ixz)
4173                   alfaxz1_av_nt(ixz)=alfaxz1_av(ixz)
4174                   alfaxz2_av_nt(ixz)=alfaxz2_av(ixz)
4175                   alfaxz3_av_nt(ixz)=alfaxz3_av(ixz)
4176                   nnn=nnn+ptcloud_xz(ixz)              
4177                enddo
4178                do ipt=1,min(ntrpt_max_nt,nnn)
4179                  tr_cloud_nt(ipt)=tr_cloud(ipt)
4180                 enddo
4181             endif
4182          endif
4183          end
4184          

Legend:
Removed from v.1.8  
changed lines
  Added in v.1.47

  ViewVC Help
Powered by ViewVC 1.1.23