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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.15 - (hide annotations) (download)
Tue Nov 21 17:13:31 2006 UTC (18 years ago) by pam-fi
Branch: MAIN
Changes since 1.14: +50 -19 lines
some refinements

1 mocchiut 1.1 ************************************************************
2     * The following subroutines
3     * - track_finding >> hough transform
4     * - track_fitting >> bob golden fitting
5     * all the procedures to create LEVEL2 data, starting from LEVEL1 data.
6     *
7     *
8     *
9     * (This subroutine and all the dependent subroutines
10     * will be included in the flight software)
11     ************************************************************
12     subroutine track_finding(iflag)
13    
14     include 'commontracker.f'
15 pam-fi 1.9 include 'level1.f'
16 mocchiut 1.1 include 'common_momanhough.f'
17     include 'common_mech.f'
18     include 'common_xyzPAM.f'
19     include 'common_mini_2.f'
20     include 'calib.f'
21     include 'level2.f'
22    
23 pam-fi 1.9 c include 'momanhough_init.f'
24 mocchiut 1.1
25     *-------------------------------------------------------------------------------
26     * STEP 1
27     *-------------------------------------------------------------------------------
28     * X-Y cluster association
29     *
30     * Clusters are associated to form COUPLES
31     * Clusters not associated in any couple are called SINGLETS
32     *
33     * Track identification (Hough transform) and fitting is first done on couples.
34     * Hence singlets are possibly added to the track.
35     *
36     * Variables assigned by the routine "cl_to_couples" are those in the
37     * common blocks:
38     * - common/clusters/cl_good
39     * - common/couples/clx,cly,ncp_plane,ncp_tot,cp_useds1,cp_useds2
40     * - common/singlets/ncls,cls,cl_single
41     *-------------------------------------------------------------------------------
42     *-------------------------------------------------------------------------------
43    
44     c iflag=0
45     call cl_to_couples(iflag)
46     if(iflag.eq.1)then !bad event
47 pam-fi 1.9 goto 880 !go to next event
48 mocchiut 1.1 endif
49    
50     *-----------------------------------------------------
51     *-----------------------------------------------------
52     * HOUGH TRASFORM
53     *-----------------------------------------------------
54     *-----------------------------------------------------
55    
56    
57     *-------------------------------------------------------------------------------
58     * STEP 2
59     *-------------------------------------------------------------------------------
60     *
61     * Association of couples to form
62     * - DOUBLETS in YZ view
63     * - TRIPLETS in XZ view
64     *
65     * Variables assigned by the routine "cp_to_doubtrip" are those in the
66     * common blocks:
67     * - common/hough_param/
68     * $ alfayz1, !Y0
69     * $ alfayz2, !tg theta-yz
70     * $ alfaxz1, !X0
71     * $ alfaxz2, !tg theta-xz
72     * $ alfaxz3 !1/r
73     * - common/doublets/ndblt,cpyz1,cpyz2
74     * - common/triplets/ntrpt,cpxz1,cpxz2,cpxz3
75     *-------------------------------------------------------------------------------
76     *-------------------------------------------------------------------------------
77    
78     c iflag=0
79     call cp_to_doubtrip(iflag)
80     if(iflag.eq.1)then !bad event
81 pam-fi 1.9 goto 880 !go to next event
82 mocchiut 1.1 endif
83    
84    
85     *-------------------------------------------------------------------------------
86     * STEP 3
87     *-------------------------------------------------------------------------------
88     *
89     * Classification of doublets and triplets to form CLOUDS,
90     * according to distance in parameter space.
91     *
92     * cloud = cluster of points (doublets/triplets) in parameter space
93     *
94     *
95     *
96     * Variables assigned by the routine "doub_to_YZcloud" are those in the
97     * common blocks:
98     * - common/clouds_yz/
99     * $ nclouds_yz
100     * $ ,alfayz1_av,alfayz2_av
101     * $ ,ptcloud_yz,db_cloud,cpcloud_yz
102     *
103     * Variables assigned by the routine "trip_to_XZcloud" are those in the
104     * common blocks:
105     * common/clouds_xz/
106     * $ nclouds_xz xz2_av,alfaxz3_av
107     * $ ,ptcloud_xz,tr_cloud,cpcloud_xz
108     *-------------------------------------------------------------------------------
109     *-------------------------------------------------------------------------------
110 pam-fi 1.9 * count number of hit planes
111     planehit=0
112     do np=1,nplanes
113     if(ncp_plane(np).ne.0)then
114     planehit=planehit+1
115     endif
116     enddo
117     if(planehit.lt.3) goto 880 ! exit
118    
119     nptxz_min=x_min_start
120     nplxz_min=x_min_start
121    
122     nptyz_min=y_min_start
123     nplyz_min=y_min_start
124    
125     cutdistyz=cutystart
126     cutdistxz=cutxstart
127 mocchiut 1.1
128 pam-fi 1.9 878 continue
129 mocchiut 1.1 call doub_to_YZcloud(iflag)
130     if(iflag.eq.1)then !bad event
131     goto 880 !fill ntp and go to next event
132 pam-fi 1.9 endif
133     if(nclouds_yz.eq.0.and.cutdistyz.lt.maxcuty)then
134     if(cutdistyz.lt.maxcuty/2)then
135     cutdistyz=cutdistyz+cutystep
136     else
137     cutdistyz=cutdistyz+(3*cutystep)
138     endif
139     goto 878
140     endif
141    
142     if(planehit.eq.3) goto 881
143    
144     879 continue
145 mocchiut 1.1 call trip_to_XZcloud(iflag)
146     if(iflag.eq.1)then !bad event
147     goto 880 !fill ntp and go to next event
148     endif
149 pam-fi 1.9
150     if(nclouds_xz.eq.0.and.cutdistxz.lt.maxcutx)then
151     cutdistxz=cutdistxz+cutxstep
152     goto 879
153     endif
154    
155    
156     881 continue
157     * if there is at least three planes on the Y view decreases cuts on X view
158     if(nclouds_xz.eq.0.and.nclouds_yz.gt.0.and.
159     $ nplxz_min.ne.y_min_start)then
160     nptxz_min=x_min_step
161     nplxz_min=x_min_start-x_min_step
162     goto 879
163     endif
164    
165 mocchiut 1.1 880 return
166     end
167    
168     ************************************************************
169    
170    
171     subroutine track_fitting(iflag)
172    
173     include 'commontracker.f'
174 pam-fi 1.9 include 'level1.f'
175 mocchiut 1.1 include 'common_momanhough.f'
176     include 'common_mech.f'
177     include 'common_xyzPAM.f'
178     include 'common_mini_2.f'
179     include 'calib.f'
180     include 'level2.f'
181    
182 pam-fi 1.9 c include 'momanhough_init.f'
183 mocchiut 1.1
184     logical FIMAGE !
185 pam-fi 1.11 real*8 AL_GUESS(5)
186 mocchiut 1.1
187     *-------------------------------------------------------------------------------
188     * STEP 4 (ITERATED until any other physical track isn't found)
189     *-------------------------------------------------------------------------------
190     *
191     * YZ and XZ clouds are combined in order to obtain the initial guess
192     * of the candidate-track parameters.
193     * A minimum number of matching couples between YZ and XZ clouds is required.
194     *
195     * A TRACK CANDIDATE is defined by
196     * - the couples resulting from the INTERSECTION of the two clouds, and
197     * - the associated track parameters (evaluated by performing a zero-order
198     * track fitting)
199     *
200     * The NTRACKS candidate-track parameters are stored in common block:
201     *
202     * - common/track_candidates/NTRACKS,AL_STORE
203     * $ ,XV_STORE,YV_STORE,ZV_STORE
204     * $ ,XM_STORE,YM_STORE,ZM_STORE
205     * $ ,RESX_STORE,RESY_STORE
206     * $ ,AXV_STORE,AYV_STORE
207     * $ ,XGOOD_STORE,YGOOD_STORE
208     * $ ,CP_STORE,RCHI2_STORE
209     *
210     *-------------------------------------------------------------------------------
211     *-------------------------------------------------------------------------------
212     ntrk=0 !counter of identified physical tracks
213    
214     11111 continue !<<<<<<< come here when performing a new search
215    
216     c iflag=0
217     call clouds_to_ctrack(iflag)
218     if(iflag.eq.1)then !no candidate tracks found
219     goto 880 !fill ntp and go to next event
220     endif
221    
222     FIMAGE=.false. !processing best track (not track image)
223     ibest=0 !best track among candidates
224     iimage=0 !track image
225     * ------------- select the best track -------------
226 pam-fi 1.9 c$$$ rchi2best=1000000000.
227     c$$$ do i=1,ntracks
228     c$$$ if(RCHI2_STORE(i).lt.rchi2best.and.
229     c$$$ $ RCHI2_STORE(i).gt.0)then
230     c$$$ ibest=i
231     c$$$ rchi2best=RCHI2_STORE(i)
232     c$$$ endif
233     c$$$ enddo
234     c$$$ if(ibest.eq.0)goto 880 !>> no good candidates
235    
236 pam-fi 1.15 * -------------------------------------------------------
237     * order track-candidates according to:
238     * 1st) decreasing n.points
239     * 2nd) increasing chi**2
240     * -------------------------------------------------------
241     rchi2best=1000000000.
242 pam-fi 1.9 ndofbest=0 !(1)
243 mocchiut 1.1 do i=1,ntracks
244 pam-fi 1.15 ndof=0 !(1)
245     do ii=1,nplanes !(1)
246     ndof=ndof !(1)
247     $ +int(xgood_store(ii,i)) !(1)
248     $ +int(ygood_store(ii,i)) !(1)
249     enddo !(1)
250     if(ndof.gt.ndofbest)then !(1)
251     ibest=i
252     rchi2best=RCHI2_STORE(i)
253     ndofbest=ndof !(1)
254     elseif(ndof.eq.ndofbest)then !(1)
255     if(RCHI2_STORE(i).lt.rchi2best.and.
256     $ RCHI2_STORE(i).gt.0)then
257 mocchiut 1.1 ibest=i
258     rchi2best=RCHI2_STORE(i)
259 pam-fi 1.9 ndofbest=ndof !(1)
260     endif !(1)
261     endif
262 mocchiut 1.1 enddo
263 pam-fi 1.15
264     c$$$ rchi2best=1000000000.
265     c$$$ ndofbest=0 !(1)
266     c$$$ do i=1,ntracks
267     c$$$ if(RCHI2_STORE(i).lt.rchi2best.and.
268     c$$$ $ RCHI2_STORE(i).gt.0)then
269     c$$$ ndof=0 !(1)
270     c$$$ do ii=1,nplanes !(1)
271     c$$$ ndof=ndof !(1)
272     c$$$ $ +int(xgood_store(ii,i)) !(1)
273     c$$$ $ +int(ygood_store(ii,i)) !(1)
274     c$$$ enddo !(1)
275     c$$$ if(ndof.ge.ndofbest)then !(1)
276     c$$$ ibest=i
277     c$$$ rchi2best=RCHI2_STORE(i)
278     c$$$ ndofbest=ndof !(1)
279     c$$$ endif !(1)
280     c$$$ endif
281     c$$$ enddo
282    
283 mocchiut 1.1 if(ibest.eq.0)goto 880 !>> no good candidates
284     *-------------------------------------------------------------------------------
285     * The best track candidate (ibest) is selected and a new fitting is performed.
286     * Previous to this, the track is refined by:
287     * - possibly adding new COUPLES or SINGLETS from the missing planes
288     * - evaluating the coordinates with improved PFAs
289     * ( angle-dependent ETA algorithms )
290     *-------------------------------------------------------------------------------
291    
292     1212 continue !<<<<< come here to fit track-image
293    
294     if(.not.FIMAGE)then !processing best candidate
295     icand=ibest
296     else !processing image
297     icand=iimage
298     iimage=0
299     endif
300     if(icand.eq.0)then
301     print*,'HAI FATTO UN CASINO!!!!!! icand = ',icand
302     $ ,ibest,iimage
303     return
304     endif
305    
306     * *-*-*-*-*-*-*-*-*-*-*-*-*-*-*
307     call refine_track(icand)
308     * *-*-*-*-*-*-*-*-*-*-*-*-*-*-*
309    
310     * **********************************************************
311     * ************************** FIT *** FIT *** FIT *** FIT ***
312     * **********************************************************
313 pam-fi 1.11 call guess()
314     do i=1,5
315     AL_GUESS(i)=AL(i)
316     enddo
317    
318 mocchiut 1.1 do i=1,5
319 pam-fi 1.3 AL(i)=dble(AL_STORE(i,icand))
320 mocchiut 1.1 enddo
321 pam-fi 1.11
322 pam-fi 1.3 IDCAND = icand !fitted track-candidate
323 mocchiut 1.1 ifail=0 !error flag in chi2 computation
324     jstep=0 !# minimization steps
325 pam-fi 1.9
326 pam-fi 1.8 iprint=0
327 pam-fi 1.11 c if(DEBUG)iprint=1
328 pam-fi 1.13 if(VERBOSE)iprint=1
329 pam-fi 1.8 call mini2(jstep,ifail,iprint)
330 mocchiut 1.1 if(ifail.ne.0) then
331 pam-fi 1.11 if(.true.)then
332 mocchiut 1.1 print *,
333 pam-fi 1.11 $ '*** MINIMIZATION FAILURE *** (after refinement) '
334 mocchiut 1.1 $ ,iev
335 pam-fi 1.11
336 pam-fi 1.12 c$$$ print*,'guess: ',(al_guess(i),i=1,5)
337     c$$$ print*,'previous: ',(al_store(i,icand),i=1,5)
338     c$$$ print*,'result: ',(al(i),i=1,5)
339     c$$$ print*,'xgood ',xgood
340     c$$$ print*,'ygood ',ygood
341     c$$$ print*,'----------------------------------------------'
342 mocchiut 1.1 endif
343 pam-fi 1.11 c chi2=-chi2
344 mocchiut 1.1 endif
345    
346     if(DEBUG)then
347     print*,'----------------------------- improved track coord'
348     22222 format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)
349     do ip=1,6
350     write(*,22222)ip,zm(ip),xm(ip),ym(ip)
351     $ ,xm_A(ip),ym_A(ip),xm_B(ip),ym_B(ip)
352     $ ,xgood(ip),ygood(ip),resx(ip),resy(ip)
353     enddo
354     endif
355    
356     c rchi2=chi2/dble(ndof)
357     if(DEBUG)then
358     print*,' '
359     print*,'****** SELECTED TRACK *************'
360     print*,'# R. chi2 RIG'
361     print*,' --- ',chi2,' --- '
362     $ ,1./abs(AL(5))
363     print*,'***********************************'
364     endif
365     * **********************************************************
366     * ************************** FIT *** FIT *** FIT *** FIT ***
367     * **********************************************************
368    
369    
370     * ------------- search if the track has an IMAGE -------------
371     * ------------- (also this is stored ) -------------
372     if(FIMAGE)goto 122 !>>> jump! (this is already an image)
373     * now search for track-image, by comparing couples IDs
374     do i=1,ntracks
375     iimage=i
376     do ip=1,nplanes
377     if( CP_STORE(nplanes-ip+1,icand).ne.
378     $ -1*CP_STORE(nplanes-ip+1,i) )iimage=0
379     enddo
380     if( iimage.ne.0.and.
381     c $ RCHI2_STORE(i).le.CHI2MAX.and.
382     c $ RCHI2_STORE(i).gt.0.and.
383     $ .true.)then
384     if(DEBUG)print*,'Track candidate ',iimage
385     $ ,' >>> TRACK IMAGE >>> of'
386     $ ,ibest
387     goto 122 !image track found
388     endif
389     enddo
390     122 continue
391    
392     * --- and store the results --------------------------------
393     ntrk = ntrk + 1 !counter of found tracks
394     if(.not.FIMAGE
395     $ .and.iimage.eq.0) image(ntrk)= 0
396     if(.not.FIMAGE
397     $ .and.iimage.ne.0)image(ntrk)=ntrk+1 !this is the image of the next
398     if(FIMAGE) image(ntrk)=ntrk-1 !this is the image of the previous
399     call fill_level2_tracks(ntrk) !==> good2=.true.
400     c print*,'++++++++++ iimage,fimage,ntrk,image '
401     c $ ,iimage,fimage,ntrk,image(ntrk)
402    
403     if(ntrk.eq.NTRKMAX)then
404 pam-fi 1.5 if(verbose)
405 mocchiut 1.1 $ print*,
406     $ '** warning ** number of identified '//
407     $ 'tracks exceeds vector dimension '
408     $ ,'( ',NTRKMAX,' )'
409     cc good2=.false.
410     goto 880 !fill ntp and go to next event
411     endif
412     if(iimage.ne.0)then
413     FIMAGE=.true. !
414     goto 1212 !>>> fit image-track
415     endif
416    
417     * --- then remove selected clusters (ibest+iimage) from clouds ----
418     call clean_XYclouds(ibest,iflag)
419     if(iflag.eq.1)then !bad event
420     goto 880 !fill ntp and go to next event
421     endif
422    
423     * **********************************************************
424     * condition to start a new search
425     * **********************************************************
426     ixznew=0
427     do ixz=1,nclouds_xz
428     if(ptcloud_xz(ixz).ge.nptxz_min)ixznew=1
429     enddo
430     iyznew=0
431     do iyz=1,nclouds_yz
432     if(ptcloud_yz(iyz).ge.nptyz_min)iyznew=1
433     enddo
434    
435     if(ixznew.ne.0.and.
436     $ iyznew.ne.0.and.
437     $ rchi2best.le.CHI2MAX.and.
438     c $ rchi2best.lt.15..and.
439     $ .true.)then
440     if(DEBUG)then
441     print*,'***** NEW SEARCH ****'
442     endif
443     goto 11111 !try new search
444    
445     endif
446     * **********************************************
447    
448    
449    
450     880 return
451     end
452    
453    
454    
455     ************************************************************
456     ************************************************************
457     ************************************************************
458     ************************************************************
459     *
460     * This routine provides the coordinates (in cm) in the PAMELA reference system:
461     * - of the point associated with a COUPLE ---> (xPAM,yPAM,zPAM)
462     * - of the extremes of the segment
463     * associated with a SINGLET ---------------> (xPAM_A,yPAM_A,zPAM_A)
464     * ---> (xPAM_B,yPAM_B,zPAM_B)
465     *
466     * It also assigns the spatial resolution to the evaluated coordinates,
467     * as a function (in principle) of the multiplicity, the angle, the PFA etc...
468     *
469     *
470     * To call the routine you must pass the arguments:
471     * icx - ID of cluster x
472     * icy - ID of cluster y
473     * sensor - sensor (1,2)
474     * PFAx - Position Finding Algorithm in x (COG2,ETA2,...)
475     * PFAy - Position Finding Algorithm in y (COG2,ETA2,...)
476     * angx - Projected angle in x
477     * angy - Projected angle in y
478     *
479     * --------- COUPLES -------------------------------------------------------
480     * The couple defines a point in the space.
481     * The coordinates of the point are evaluated as follows:
482     * 1 - the corrected coordinates relative to the sensor are evaluated
483     * according to the chosen PFA --> (xi,yi,0)
484     * 2 - coordinates are rotated and traslated, according to the aligmnet
485     * parameters, and expressed in the reference system of the mechanical
486     * sensor --> (xrt,yrt,zrt)
487     * 3 - coordinates are finally converted to the PAMELA reference system
488     * --> (xPAM,yPAM,zPAM)
489     *
490     * --------- SINGLETS -------------------------------------------------------
491     * Since a coordinate is missing, the singlet defines not a point
492     * in the space but a segment AB (parallel to the strips).
493     * In this case the routine returns the coordinates in the PAMELA reference
494     * system of the two extremes A and B of the segment:
495     * --> (xPAM_A,yPAM_A,zPAM_A)
496     * --> (xPAM_B,yPAM_B,zPAM_B)
497     *
498     * ==========================================================
499     *
500     * The output of the routine is stored in the commons:
501     *
502     * double precision xPAM,yPAM,zPAM
503     * common/coord_xyz_PAM/xPAM,yPAM,zPAM
504     *
505     * double precision xPAM_A,yPAM_A,zPAM_A
506     * double precision xPAM_B,yPAM_B,zPAM_B
507     * common/coord_AB_PAM/xPAM_A,yPAM_A,zPAM_A,xPAM_B,yPAM_B,zPAM_B
508     *
509     * double precision resxPAM,resyPAM
510     * common/resolution_PAM/resxPAM,resyPAM
511     *
512     * (in file common_xyzPAM.f)
513     *
514     *
515    
516     subroutine xyz_PAM(icx,icy,sensor,PFAx,PFAy,angx,angy)
517    
518     c*****************************************************
519     c 07/10/2005 modified by elena vannuccini --> (1)
520     c 01/02/2006 modified by elena vannuccini --> (2)
521     c 02/02/2006 modified by Elena Vannuccini --> (3)
522     c (implemented new p.f.a.)
523     c 03/02/2006 modified by Elena Vannuccini --> (4)
524     c (implemented variable resolution)
525     c*****************************************************
526    
527     include 'commontracker.f'
528 pam-fi 1.9 include 'level1.f'
529 mocchiut 1.1 include 'calib.f'
530 pam-fi 1.9 c include 'level1.f'
531 mocchiut 1.1 include 'common_align.f'
532     include 'common_mech.f'
533     include 'common_xyzPAM.f'
534 pam-fi 1.11 c include 'common_resxy.f'
535 mocchiut 1.1
536 pam-fi 1.2 c logical DEBUG
537     c common/dbg/DEBUG
538 mocchiut 1.1
539     integer icx,icy !X-Y cluster ID
540     integer sensor
541     integer viewx,viewy
542     character*4 PFAx,PFAy !PFA to be used
543     real angx,angy !X-Y angle
544    
545     real stripx,stripy
546    
547     double precision xrt,yrt,zrt
548     double precision xrt_A,yrt_A,zrt_A
549     double precision xrt_B,yrt_B,zrt_B
550     c double precision xi,yi,zi
551     c double precision xi_A,yi_A,zi_A
552     c double precision xi_B,yi_B,zi_B
553    
554    
555     parameter (ndivx=30)
556    
557     resxPAM = 0
558     resyPAM = 0
559    
560     xPAM = 0.
561     yPAM = 0.
562     zPAM = 0.
563     xPAM_A = 0.
564     yPAM_A = 0.
565     zPAM_A = 0.
566     xPAM_B = 0.
567     yPAM_B = 0.
568     zPAM_B = 0.
569    
570     * -----------------
571     * CLUSTER X
572     * -----------------
573    
574     if(icx.ne.0)then
575     viewx = VIEW(icx)
576     nldx = nld(MAXS(icx),VIEW(icx))
577     nplx = npl(VIEW(icx))
578     resxPAM = RESXAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!!
579    
580     stripx = float(MAXS(icx))
581     if(PFAx.eq.'COG1')then !(1)
582     stripx = stripx !(1)
583     resxPAM = resxPAM !(1)
584     elseif(PFAx.eq.'COG2')then
585     stripx = stripx + cog(2,icx)
586     resxPAM = resxPAM*fbad_cog(2,icx)
587     elseif(PFAx.eq.'ETA2')then
588     c cog2 = cog(2,icx)
589 pam-fi 1.7 c etacorr = pfaeta2(cog2,viewx,nldx,angx)
590 mocchiut 1.1 c stripx = stripx + etacorr
591 pam-fi 1.7 stripx = stripx + pfaeta2(icx,angx) !(3)
592 mocchiut 1.1 resxPAM = risx_eta2(angx) ! (4)
593     if(DEBUG.and.fbad_cog(2,icx).ne.1)
594     $ print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)
595     resxPAM = resxPAM*fbad_cog(2,icx)
596     elseif(PFAx.eq.'ETA3')then !(3)
597 pam-fi 1.7 stripx = stripx + pfaeta3(icx,angx) !(3)
598 mocchiut 1.1 resxPAM = risx_eta3(angx) ! (4)
599     if(DEBUG.and.fbad_cog(3,icx).ne.1) !(3)
600     $ print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)!(3)
601     resxPAM = resxPAM*fbad_cog(3,icx) !(3)
602     elseif(PFAx.eq.'ETA4')then !(3)
603 pam-fi 1.7 stripx = stripx + pfaeta4(icx,angx) !(3)
604 mocchiut 1.1 resxPAM = risx_eta4(angx) ! (4)
605     if(DEBUG.and.fbad_cog(4,icx).ne.1) !(3)
606     $ print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)!(3)
607     resxPAM = resxPAM*fbad_cog(4,icx) !(3)
608     elseif(PFAx.eq.'ETA')then !(3)
609 pam-fi 1.7 stripx = stripx + pfaeta(icx,angx) !(3)
610 mocchiut 1.1 resxPAM = ris_eta(icx,angx) ! (4)
611     if(DEBUG.and.fbad_cog(2,icx).ne.1) !(3)
612     $ print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)!(3)
613     c resxPAM = resxPAM*fbad_cog(2,icx) !(3)TEMPORANEO
614     resxPAM = resxPAM*fbad_eta(icx,angx) !(3)(4)
615     elseif(PFAx.eq.'COG')then !(2)
616     stripx = stripx + cog(0,icx) !(2)
617     resxPAM = risx_cog(angx) ! (4)
618     resxPAM = resxPAM*fbad_cog(0,icx)!(2)
619     else
620     print*,'*** Non valid p.f.a. (x) --> ',PFAx
621     endif
622    
623     endif
624 pam-fi 1.4 c if(icy.eq.0.and.icx.ne.0)
625     c $ print*,PFAx,icx,angx,stripx,resxPAM,'***'
626 mocchiut 1.1
627     * -----------------
628     * CLUSTER Y
629     * -----------------
630    
631     if(icy.ne.0)then
632     viewy = VIEW(icy)
633     nldy = nld(MAXS(icy),VIEW(icy))
634     nply = npl(VIEW(icy))
635     resyPAM = RESYAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!!
636    
637    
638     if(icx.ne.0.and.(nply.ne.nplx.or.nldy.ne.nldx))then
639     print*,'xyz_PAM ***ERROR*** invalid cluster couple!!! '
640     $ ,icx,icy
641     goto 100
642     endif
643    
644     stripy = float(MAXS(icy))
645     if(PFAy.eq.'COG1')then !(1)
646     stripy = stripy !(1)
647     resyPAM = resyPAM !(1)
648     elseif(PFAy.eq.'COG2')then
649     stripy = stripy + cog(2,icy)
650     resyPAM = resyPAM*fbad_cog(2,icy)
651     elseif(PFAy.eq.'ETA2')then
652     c cog2 = cog(2,icy)
653 pam-fi 1.7 c etacorr = pfaeta2(cog2,viewy,nldy,angy)
654 mocchiut 1.1 c stripy = stripy + etacorr
655 pam-fi 1.7 stripy = stripy + pfaeta2(icy,angy) !(3)
656 mocchiut 1.1 resyPAM = risy_eta2(angy) ! (4)
657     resyPAM = resyPAM*fbad_cog(2,icy)
658     if(DEBUG.and.fbad_cog(2,icy).ne.1)
659     $ print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)
660     elseif(PFAy.eq.'ETA3')then !(3)
661 pam-fi 1.7 stripy = stripy + pfaeta3(icy,angy) !(3)
662 mocchiut 1.1 resyPAM = resyPAM*fbad_cog(3,icy) !(3)
663     if(DEBUG.and.fbad_cog(3,icy).ne.1) !(3)
664     $ print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)!(3)
665     elseif(PFAy.eq.'ETA4')then !(3)
666 pam-fi 1.7 stripy = stripy + pfaeta4(icy,angy) !(3)
667 mocchiut 1.1 resyPAM = resyPAM*fbad_cog(4,icy) !(3)
668     if(DEBUG.and.fbad_cog(4,icy).ne.1) !(3)
669     $ print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)!(3)
670     elseif(PFAy.eq.'ETA')then !(3)
671 pam-fi 1.7 stripy = stripy + pfaeta(icy,angy) !(3)
672 mocchiut 1.1 resyPAM = ris_eta(icy,angy) ! (4)
673     c resyPAM = resyPAM*fbad_cog(2,icy) !(3)TEMPORANEO
674     resyPAM = resyPAM*fbad_eta(icy,angy) ! (4)
675     if(DEBUG.and.fbad_cog(2,icy).ne.1) !(3)
676     $ print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)!(3)
677     elseif(PFAy.eq.'COG')then
678     stripy = stripy + cog(0,icy)
679     resyPAM = risy_cog(angy) ! (4)
680     c resyPAM = ris_eta(icy,angy) ! (4)
681     resyPAM = resyPAM*fbad_cog(0,icy)
682     else
683     print*,'*** Non valid p.f.a. (x) --> ',PFAx
684     endif
685    
686     endif
687    
688    
689     c===========================================================
690     C COUPLE
691     C===========================================================
692     if(icx.ne.0.and.icy.ne.0)then
693    
694     c------------------------------------------------------------------------
695     c (xi,yi,zi) = mechanical coordinates in the silicon sensor frame
696     c------------------------------------------------------------------------
697 pam-fi 1.4 if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
698     $ .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
699     print*,'xyz_PAM (couple):',
700     $ ' WARNING: false X strip: strip ',stripx
701     endif
702 mocchiut 1.1 xi = acoordsi(stripx,viewx)
703     yi = acoordsi(stripy,viewy)
704     zi = 0.
705    
706    
707     c------------------------------------------------------------------------
708     c (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame
709     c------------------------------------------------------------------------
710     c N.B. I convert angles from microradiants to radiants
711    
712     xrt = xi
713     $ - omega(nplx,nldx,sensor)*yi
714     $ + gamma(nplx,nldx,sensor)*zi
715     $ + dx(nplx,nldx,sensor)
716    
717     yrt = omega(nplx,nldx,sensor)*xi
718     $ + yi
719     $ - beta(nplx,nldx,sensor)*zi
720     $ + dy(nplx,nldx,sensor)
721    
722     zrt = -gamma(nplx,nldx,sensor)*xi
723     $ + beta(nplx,nldx,sensor)*yi
724     $ + zi
725     $ + dz(nplx,nldx,sensor)
726    
727     c xrt = xi
728     c yrt = yi
729     c zrt = zi
730    
731     c------------------------------------------------------------------------
732     c (xPAM,yPAM,zPAM) = measured coordinates (in cm)
733     c in PAMELA reference system
734     c------------------------------------------------------------------------
735    
736     xPAM = dcoord(xrt,viewx,nldx,sensor) / 1.d4
737     yPAM = dcoord(yrt,viewy,nldy,sensor) / 1.d4
738     zPAM = ( zrt + z_mech_sensor(nplx,nldx,sensor)*1000. ) / 1.d4
739    
740     xPAM_A = 0.
741     yPAM_A = 0.
742     zPAM_A = 0.
743    
744     xPAM_B = 0.
745     yPAM_B = 0.
746     zPAM_B = 0.
747    
748     elseif(
749     $ (icx.ne.0.and.icy.eq.0).or.
750     $ (icx.eq.0.and.icy.ne.0).or.
751     $ .false.
752     $ )then
753    
754     c------------------------------------------------------------------------
755     c (xi,yi,zi) = mechanical coordinates in the silicon sensor frame
756     c------------------------------------------------------------------------
757    
758     if(icy.ne.0)then
759     c===========================================================
760     C Y-SINGLET
761     C===========================================================
762     nplx = nply
763     nldx = nldy
764     viewx = viewy + 1
765    
766     yi = acoordsi(stripy,viewy)
767    
768     xi_A = edgeY_d - SiDimX/2
769     yi_A = yi
770     zi_A = 0.
771    
772     xi_B = SiDimX/2 - edgeY_u
773     yi_B = yi
774     zi_B = 0.
775    
776     c print*,'Y-cl ',icy,stripy,' --> ',yi
777     c print*,xi_A,' <--> ',xi_B
778    
779     elseif(icx.ne.0)then
780     c===========================================================
781     C X-SINGLET
782     C===========================================================
783    
784     nply = nplx
785     nldy = nldx
786     viewy = viewx - 1
787    
788 pam-fi 1.4 c print*,'X-singlet ',icx,nplx,nldx,viewx,stripx
789     c if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips...
790     if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
791     $ .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
792     print*,'xyz_PAM (X-singlet):',
793     $ ' WARNING: false X strip: strip ',stripx
794     endif
795 mocchiut 1.1 xi = acoordsi(stripx,viewx)
796    
797     xi_A = xi
798     yi_A = edgeX_d - SiDimY/2
799     zi_A = 0.
800    
801     xi_B = xi
802     yi_B = SiDimY/2 - edgeX_u
803     zi_B = 0.
804    
805     if(viewy.eq.11)then
806     yi = yi_A
807     yi_A = yi_B
808     yi_B = yi
809     endif
810    
811     c print*,'X-cl ',icx,stripx,' --> ',xi
812     c print*,yi_A,' <--> ',yi_B
813    
814     else
815    
816     print *,'routine xyz_PAM ---> not properly used !!!'
817     print *,'icx = ',icx
818     print *,'icy = ',icy
819     goto 100
820    
821     endif
822     c------------------------------------------------------------------------
823     c (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame
824     c------------------------------------------------------------------------
825     c N.B. I convert angles from microradiants to radiants
826    
827     xrt_A = xi_A
828     $ - omega(nplx,nldx,sensor)*yi_A
829     $ + gamma(nplx,nldx,sensor)*zi_A
830     $ + dx(nplx,nldx,sensor)
831    
832     yrt_A = omega(nplx,nldx,sensor)*xi_A
833     $ + yi_A
834     $ - beta(nplx,nldx,sensor)*zi_A
835     $ + dy(nplx,nldx,sensor)
836    
837     zrt_A = -gamma(nplx,nldx,sensor)*xi_A
838     $ + beta(nplx,nldx,sensor)*yi_A
839     $ + zi_A
840     $ + dz(nplx,nldx,sensor)
841    
842     xrt_B = xi_B
843     $ - omega(nplx,nldx,sensor)*yi_B
844     $ + gamma(nplx,nldx,sensor)*zi_B
845     $ + dx(nplx,nldx,sensor)
846    
847     yrt_B = omega(nplx,nldx,sensor)*xi_B
848     $ + yi_B
849     $ - beta(nplx,nldx,sensor)*zi_B
850     $ + dy(nplx,nldx,sensor)
851    
852     zrt_B = -gamma(nplx,nldx,sensor)*xi_B
853     $ + beta(nplx,nldx,sensor)*yi_B
854     $ + zi_B
855     $ + dz(nplx,nldx,sensor)
856    
857    
858     c xrt = xi
859     c yrt = yi
860     c zrt = zi
861    
862     c------------------------------------------------------------------------
863     c (xPAM,yPAM,zPAM) = measured coordinates (in cm)
864     c in PAMELA reference system
865     c------------------------------------------------------------------------
866    
867     xPAM = 0.
868     yPAM = 0.
869     zPAM = 0.
870    
871     xPAM_A = dcoord(xrt_A,viewx,nldx,sensor) / 1.d4
872     yPAM_A = dcoord(yrt_A,viewy,nldy,sensor) / 1.d4
873     zPAM_A = ( zrt_A + z_mech_sensor(nplx,nldx,sensor)*1000.)/ 1.d4
874    
875     xPAM_B = dcoord(xrt_B,viewx,nldx,sensor) / 1.d4
876     yPAM_B = dcoord(yrt_B,viewy,nldy,sensor) / 1.d4
877     zPAM_B = ( zrt_B + z_mech_sensor(nplx,nldx,sensor)*1000.)/ 1.d4
878    
879    
880     c print*,'A-(',xPAM_A,yPAM_A,') B-(',xPAM_B,yPAM_B,')'
881    
882     else
883    
884     print *,'routine xyz_PAM ---> not properly used !!!'
885     print *,'icx = ',icx
886     print *,'icy = ',icy
887    
888     endif
889    
890     100 continue
891     end
892    
893    
894     ********************************************************************************
895     ********************************************************************************
896     ********************************************************************************
897     *
898     * The function distance_to(XP,YP) should be used after
899     * a call to the xyz_PAM routine and it evaluate the
900     * NORMALIZED distance (PROJECTED on the XY plane) between
901     * the point (XP,YP), argument of the function,
902     * and:
903     *
904     * - the point (xPAM,yPAM,zPAM), in the case of a COUPLE
905     * or
906     * - the segment (xPAM_A,yPAM_A,zPAM_A)-(xPAM_B,yPAM_B,zPAM_B),
907     * in the case of a SINGLET.
908     *
909     * ( The routine xyz_PAM fills the common defined in "common_xyzPAM.f",
910     * which stores the coordinates of the couple/singlet )
911     *
912     ********************************************************************************
913    
914     real function distance_to(XPP,YPP)
915    
916     include 'common_xyzPAM.f'
917    
918     * -----------------------------------
919     * it computes the normalized distance
920     * ( i.e. distance/resolution )
921     * -----------------------------------
922    
923     double precision distance,RE
924     double precision BETA,ALFA,xmi,ymi
925    
926     * ----------------------
927     if (
928     + xPAM.eq.0.and.
929     + yPAM.eq.0.and.
930     + zPAM.eq.0.and.
931     + xPAM_A.ne.0.and.
932     + yPAM_A.ne.0.and.
933     + zPAM_A.ne.0.and.
934     + xPAM_B.ne.0.and.
935     + yPAM_B.ne.0.and.
936     + zPAM_B.ne.0.and.
937     + .true.)then
938     * -----------------------
939     * DISTANCE TO --- SINGLET
940     * -----------------------
941     if(abs(sngl(xPAM_B-xPAM_A)).lt.abs(sngl(yPAM_B-yPAM_A)))then
942     * |||---------- X CLUSTER
943    
944     BETA = (xPAM_B-xPAM_A)/(yPAM_B-yPAM_A)
945     ALFA = xPAM_A - BETA * yPAM_A
946    
947     ymi = ( YPP + BETA*XPP - BETA*ALFA )/(1+BETA**2)
948     if(ymi.lt.dmin1(yPAM_A,yPAM_B))ymi=dmin1(yPAM_A,yPAM_B)
949     if(ymi.gt.dmax1(yPAM_A,yPAM_B))ymi=dmax1(yPAM_A,yPAM_B)
950     xmi = ALFA + BETA * ymi
951     RE = resxPAM
952    
953     else
954     * |||---------- Y CLUSTER
955    
956     BETA = (yPAM_B-yPAM_A)/(xPAM_B-xPAM_A)
957     ALFA = yPAM_A - BETA * xPAM_A
958    
959     xmi = ( XPP + BETA*YPP - BETA*ALFA )/(1+BETA**2)
960     if(xmi.lt.dmin1(xPAM_A,xPAM_B))xmi=dmin1(xPAM_A,xPAM_B)
961     if(xmi.gt.dmax1(xPAM_A,xPAM_B))xmi=dmax1(xPAM_A,xPAM_B)
962     ymi = ALFA + BETA * xmi
963     RE = resyPAM
964    
965     endif
966    
967     distance=
968     $ ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2
969     distance=dsqrt(distance)
970    
971     c$$$ print*,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b
972     c$$$ $ ,' --- distance_to --- ',xpp,ypp
973     c$$$ print*,' resolution ',re
974    
975    
976     * ----------------------
977     elseif(
978     + xPAM.ne.0.and.
979     + yPAM.ne.0.and.
980     + zPAM.ne.0.and.
981     + xPAM_A.eq.0.and.
982     + yPAM_A.eq.0.and.
983     + zPAM_A.eq.0.and.
984     + xPAM_B.eq.0.and.
985     + yPAM_B.eq.0.and.
986     + zPAM_B.eq.0.and.
987     + .true.)then
988     * ----------------------
989     * DISTANCE TO --- COUPLE
990     * ----------------------
991    
992     distance=
993     $ ((xPAM-XPP)/resxPAM)**2
994     $ +
995     $ ((yPAM-YPP)/resyPAM)**2
996     distance=dsqrt(distance)
997    
998     c$$$ print*,xPAM,yPAM,zPAM
999     c$$$ $ ,' --- distance_to --- ',xpp,ypp
1000     c$$$ print*,' resolution ',resxPAM,resyPAM
1001    
1002     else
1003    
1004     print*
1005     $ ,' function distance_to ---> wrong usage!!!'
1006     print*,' xPAM,yPAM,zPAM ',xPAM,yPAM,zPAM
1007     print*,' xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b '
1008     $ ,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b
1009     endif
1010    
1011     distance_to = sngl(distance)
1012    
1013     return
1014     end
1015    
1016     ********************************************************************************
1017     ********************************************************************************
1018     ********************************************************************************
1019     ********************************************************************************
1020    
1021     subroutine whichsensor(nplPAM,xPAM,yPAM,ladder,sensor)
1022     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
1023     * Given the plane (1-6 from BOTTOM to TOP!!) and the (xPAM,yPAM)
1024     * coordinates (in the PAMELA reference system), it returns
1025     * the ladder and the sensor which the point belongs to.
1026     *
1027     * The method to assign a point to a sensor consists in
1028     * - calculating the sum of the distances between the point
1029     * and the sensor edges
1030     * - requiring that it is less-equal than (SiDimX+SiDimY)
1031     *
1032     * NB -- SiDimX and SiDimY are not the dimentions of the SENSITIVE volume
1033     * but of the whole silicon sensor
1034     *
1035     * CONVENTION:
1036     * - sensor 1 is the one closest to the hybrid
1037     * - ladder 1 is the first to be read out (strips from 1 to 1024)
1038     *
1039     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
1040     include 'commontracker.f'
1041     include 'common_align.f'
1042    
1043     integer ladder,sensor,viewx,viewy
1044     real c1(4),c2(4),c3(4)
1045     data c1/1.,0.,0.,1./
1046     data c2/1.,-1.,-1.,1./
1047     data c3/1.,1.,0.,0./
1048     real*8 yvvv,xvvv
1049     double precision xi,yi,zi
1050     double precision xrt,yrt,zrt
1051     real AA,BB
1052     real yvv(4),xvv(4)
1053    
1054     * tollerance to consider the track inside the sensitive area
1055     real ptoll
1056     data ptoll/150./ !um
1057    
1058     external nviewx,nviewy,acoordsi,dcoord
1059    
1060     nplpt = nplPAM !plane
1061     viewx = nviewx(nplpt)
1062     viewy = nviewy(nplpt)
1063    
1064     do il=1,nladders_view
1065     do is=1,2
1066    
1067     do iv=1,4 !loop on sensor vertexes
1068     stripx = (il-c1(iv))*1024 + c1(iv) + c2(iv)*3
1069     stripy = (il-c3(iv))*1024 + c3(iv)
1070     c------------------------------------------------------------------------
1071     c (xi,yi,zi) = mechanical coordinates in the silicon sensor frame
1072     c------------------------------------------------------------------------
1073 pam-fi 1.4 if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
1074     $ .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
1075     c if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips...
1076     print*,'whichsensor: ',
1077     $ ' WARNING: false X strip: strip ',stripx
1078     endif
1079 mocchiut 1.1 xi = acoordsi(stripx,viewx)
1080     yi = acoordsi(stripy,viewy)
1081     zi = 0.
1082     c------------------------------------------------------------------------
1083     c (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame
1084     c------------------------------------------------------------------------
1085     c N.B. I convert angles from microradiants to radiants
1086     xrt = xi
1087     $ - omega(nplpt,il,is)*yi
1088     $ + gamma(nplpt,il,is)*zi
1089     $ + dx(nplpt,il,is)
1090     yrt = omega(nplpt,il,is)*xi
1091     $ + yi
1092     $ - beta(nplpt,il,is)*zi
1093     $ + dy(nplpt,il,is)
1094     zrt = -gamma(nplpt,il,is)*xi
1095     $ + beta(nplpt,il,is)*yi
1096     $ + zi
1097     $ + dz(nplpt,il,is)
1098     c------------------------------------------------------------------------
1099     c measured coordinates (in cm) in PAMELA reference system
1100     c------------------------------------------------------------------------
1101     yvvv = dcoord(yrt,viewy,il,is) / 1.d4
1102     xvvv = dcoord(xrt,viewx,il,is) / 1.d4
1103    
1104     yvv(iv)=sngl(yvvv)
1105     xvv(iv)=sngl(xvvv)
1106     c print*,'LADDER ',il,' SENSOR ',is,' vertexes >> '
1107     c $ ,iv,xvv(iv),yvv(iv)
1108     enddo !end loop on sensor vertexes
1109    
1110     dtot=0.
1111     do iside=1,4,2 !loop on sensor edges X
1112     iv1=iside
1113     iv2=mod(iside,4)+1
1114     * straight line passing trhough two consecutive vertexes
1115     AA = (yvv(iv1)-yvv(iv2))/(xvv(iv1)-xvv(iv2))
1116     BB = yvv(iv1) - AA*xvv(iv1)
1117     * point along the straight line closer to the track
1118     xoo = (xPAM+AA*yPAM-AA*BB)/(1+AA**2)
1119     yoo = AA*xoo + BB
1120     * sum of the distances
1121     dtot = dtot +
1122     $ sqrt((xPAM-xoo)**2+(yPAM-yoo)**2)
1123     enddo !end loop on sensor edges
1124     do iside=2,4,2 !loop on sensor edges Y
1125     iv1=iside
1126     iv2=mod(iside,4)+1
1127     * straight line passing trhough two consecutive vertexes
1128     AA = (xvv(iv1)-xvv(iv2))/(yvv(iv1)-yvv(iv2))
1129     BB = xvv(iv1) - AA*yvv(iv1)
1130     * point along the straight line closer to the track
1131     yoo = (yPAM+AA*xPAM-AA*BB)/(1+AA**2)
1132     xoo = AA*yoo + BB
1133     * sum of the distances
1134     dtot = dtot +
1135     $ sqrt((xPAM-xoo)**2+(yPAM-yoo)**2)
1136     enddo !end loop on sensor edges
1137    
1138    
1139     * half-perimeter of sensitive area
1140     Perim =
1141     $ SiDimX - edgeX_l - edgeX_r
1142     $ +SiDimY - edgeY_l - edgeY_r
1143     Perim = (Perim + ptoll)/1.e4
1144     if(dtot.le.Perim)goto 100
1145    
1146    
1147     enddo
1148     enddo
1149    
1150     ladder = 0
1151     sensor = 0
1152     goto 200
1153    
1154     100 continue
1155     ladder = il
1156     sensor = is
1157    
1158    
1159     200 return
1160     end
1161    
1162    
1163    
1164     *************************************************************************
1165    
1166     subroutine reverse(v,n,temp) !invert the order of the components of v(n) vector
1167    
1168     implicit double precision (A-H,O-Z)
1169    
1170     dimension v(*)
1171     dimension temp(*)
1172     integer i,n
1173    
1174     do i=1,n
1175     temp(i)=v(n+1-i)
1176     enddo
1177    
1178     do i=1,n
1179     v(i)=temp(i)
1180     enddo
1181    
1182     return
1183     end
1184    
1185     *************************************************************************
1186     *************************************************************************
1187     *************************************************************************
1188     *************************************************************************
1189     *************************************************************************
1190     *************************************************************************
1191     *************************************************************************
1192     *************************************************************************
1193     *************************************************************************
1194     *************************************************************************
1195     *************************************************************************
1196     *************************************************************************
1197     *************************************************************************
1198     *************************************************************************
1199     *************************************************************************
1200     integer function ip_cp(id)
1201     *
1202     * given the couple id,
1203     * it returns the plane number
1204     *
1205     include 'commontracker.f'
1206 pam-fi 1.9 include 'level1.f'
1207 mocchiut 1.1 c include 'common_analysis.f'
1208     include 'common_momanhough.f'
1209    
1210     ip_cp=0
1211     ncpp=0
1212     do ip=1,nplanes
1213     ncpp=ncpp+ncp_plane(ip)
1214     if(ncpp.ge.abs(id))then
1215     ip_cp=ip
1216     goto 100
1217     endif
1218     enddo
1219     100 continue
1220     return
1221     end
1222    
1223    
1224     integer function is_cp(id)
1225     *
1226     * given the couple id,
1227     * it returns the sensor number
1228     *
1229     is_cp=0
1230     if(id.lt.0)is_cp=1
1231     if(id.gt.0)is_cp=2
1232     if(id.eq.0)print*,'IS_CP ===> wrong couple id !!!'
1233    
1234     return
1235     end
1236    
1237    
1238     integer function icp_cp(id)
1239     *
1240     * given the couple id,
1241     * it returns the id number ON THE PLANE
1242     *
1243     include 'commontracker.f'
1244 pam-fi 1.9 include 'level1.f'
1245 mocchiut 1.1 c include 'common_analysis.f'
1246     include 'common_momanhough.f'
1247    
1248     icp_cp=0
1249    
1250     ncpp=0
1251     do ip=1,nplanes
1252     ncppold=ncpp
1253     ncpp=ncpp+ncp_plane(ip)
1254     if(ncpp.ge.abs(id))then
1255     icp_cp=abs(id)-ncppold
1256     goto 100
1257     endif
1258     enddo
1259     100 continue
1260     return
1261     end
1262    
1263    
1264    
1265     integer function id_cp(ip,icp,is)
1266     *
1267     * given a plane, a couple and the sensor
1268     * it returns the absolute couple id
1269     * negative if sensor =1
1270     * positive if sensor =2
1271     *
1272     include 'commontracker.f'
1273 pam-fi 1.9 include 'level1.f'
1274 mocchiut 1.1 c include 'calib.f'
1275     c include 'level1.f'
1276     c include 'common_analysis.f'
1277     include 'common_momanhough.f'
1278    
1279     id_cp=0
1280    
1281     if(ip.gt.1)then
1282     do i=1,ip-1
1283     id_cp = id_cp + ncp_plane(i)
1284     enddo
1285     endif
1286    
1287     id_cp = id_cp + icp
1288    
1289     if(is.eq.1) id_cp = -id_cp
1290    
1291     return
1292     end
1293    
1294    
1295    
1296    
1297     *************************************************************************
1298     *************************************************************************
1299     *************************************************************************
1300     *************************************************************************
1301     *************************************************************************
1302     *************************************************************************
1303    
1304    
1305     ***************************************************
1306     * *
1307     * *
1308     * *
1309     * *
1310     * *
1311     * *
1312     **************************************************
1313    
1314     subroutine cl_to_couples(iflag)
1315    
1316     include 'commontracker.f'
1317 pam-fi 1.9 include 'level1.f'
1318 mocchiut 1.1 include 'common_momanhough.f'
1319 pam-fi 1.9 c include 'momanhough_init.f'
1320 mocchiut 1.1 include 'calib.f'
1321 pam-fi 1.9 c include 'level1.f'
1322 mocchiut 1.1
1323     * output flag
1324     * --------------
1325     * 0 = good event
1326     * 1 = bad event
1327     * --------------
1328     integer iflag
1329    
1330 pam-fi 1.6 integer badseed,badclx,badcly
1331 mocchiut 1.1
1332     * init variables
1333     ncp_tot=0
1334     do ip=1,nplanes
1335     do ico=1,ncouplemax
1336     clx(ip,ico)=0
1337     cly(ip,ico)=0
1338     enddo
1339     ncp_plane(ip)=0
1340     do icl=1,nclstrmax_level2
1341     cls(ip,icl)=1
1342     enddo
1343     ncls(ip)=0
1344     enddo
1345     do icl=1,nclstrmax_level2
1346 pam-fi 1.9 cl_single(icl) = 1
1347     cl_good(icl) = 0
1348     enddo
1349     do iv=1,nviews
1350     ncl_view(iv) = 0
1351     mask_view(iv) = 0 !all included
1352 mocchiut 1.1 enddo
1353    
1354 pam-fi 1.9 * count number of cluster per view
1355     do icl=1,nclstr1
1356     ncl_view(VIEW(icl)) = ncl_view(VIEW(icl)) + 1
1357     enddo
1358     * mask views with too many clusters
1359     do iv=1,nviews
1360 pam-fi 1.14 if( ncl_view(iv).gt. nclusterlimit)then
1361 pam-fi 1.9 mask_view(iv) = 1
1362 pam-fi 1.11 if(DEBUG)print*,' * WARNING * cl_to_couple: n.clusters > '
1363 pam-fi 1.14 $ ,nclusterlimit,' on view ', iv,' --> masked!'
1364 pam-fi 1.9 endif
1365     enddo
1366    
1367    
1368 mocchiut 1.1 * start association
1369     ncouples=0
1370     do icx=1,nclstr1 !loop on cluster (X)
1371     if(mod(VIEW(icx),2).eq.1)goto 10
1372    
1373     * ----------------------------------------------------
1374 pam-fi 1.9 * jump masked views (X VIEW)
1375     * ----------------------------------------------------
1376     if( mask_view(VIEW(icx)).ne.0 ) goto 10
1377     * ----------------------------------------------------
1378 mocchiut 1.1 * cut on charge (X VIEW)
1379 pam-fi 1.4 * ----------------------------------------------------
1380 mocchiut 1.1 if(dedx(icx).lt.dedx_x_min)then
1381     cl_single(icx)=0
1382     goto 10
1383     endif
1384 pam-fi 1.4 * ----------------------------------------------------
1385 mocchiut 1.1 * cut BAD (X VIEW)
1386 pam-fi 1.4 * ----------------------------------------------------
1387 mocchiut 1.1 badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx)))
1388     ifirst=INDSTART(icx)
1389     if(icx.ne.nclstr1) then
1390     ilast=INDSTART(icx+1)-1
1391     else
1392     ilast=TOTCLLENGTH
1393     endif
1394 pam-fi 1.6 badclx=badseed
1395 mocchiut 1.1 do igood=-ngoodstr,ngoodstr
1396     ibad=1
1397     if((INDMAX(icx)+igood).gt.ifirst.and.
1398     $ (INDMAX(icx)+igood).lt.ilast.and.
1399     $ .true.)then
1400     ibad=BAD(VIEW(icx),
1401     $ nvk(MAXS(icx)+igood),
1402     $ nst(MAXS(icx)+igood))
1403     endif
1404 pam-fi 1.6 badclx=badclx*ibad
1405 mocchiut 1.1 enddo
1406 pam-fi 1.4 * ----------------------------------------------------
1407     * >>> eliminato il taglio sulle BAD <<<
1408     * ----------------------------------------------------
1409 mocchiut 1.1 c if(badcl.eq.0)then
1410     c cl_single(icx)=0
1411     c goto 10
1412     c endif
1413     * ----------------------------------------------------
1414    
1415     cl_good(icx)=1
1416     nplx=npl(VIEW(icx))
1417     nldx=nld(MAXS(icx),VIEW(icx))
1418    
1419     do icy=1,nclstr1 !loop on cluster (Y)
1420     if(mod(VIEW(icy),2).eq.0)goto 20
1421    
1422     * ----------------------------------------------------
1423 pam-fi 1.9 * jump masked views (Y VIEW)
1424     * ----------------------------------------------------
1425     if( mask_view(VIEW(icy)).ne.0 ) goto 20
1426    
1427     * ----------------------------------------------------
1428 mocchiut 1.1 * cut on charge (Y VIEW)
1429 pam-fi 1.4 * ----------------------------------------------------
1430 mocchiut 1.1 if(dedx(icy).lt.dedx_y_min)then
1431     cl_single(icy)=0
1432     goto 20
1433     endif
1434 pam-fi 1.4 * ----------------------------------------------------
1435 mocchiut 1.1 * cut BAD (Y VIEW)
1436 pam-fi 1.4 * ----------------------------------------------------
1437 mocchiut 1.1 badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy)))
1438     ifirst=INDSTART(icy)
1439     if(icy.ne.nclstr1) then
1440     ilast=INDSTART(icy+1)-1
1441     else
1442     ilast=TOTCLLENGTH
1443     endif
1444 pam-fi 1.6 badcly=badseed
1445 mocchiut 1.1 do igood=-ngoodstr,ngoodstr
1446     ibad=1
1447     if((INDMAX(icy)+igood).gt.ifirst.and.
1448     $ (INDMAX(icy)+igood).lt.ilast.and.
1449     $ .true.)
1450     $ ibad=BAD(VIEW(icy),
1451     $ nvk(MAXS(icy)+igood),
1452     $ nst(MAXS(icy)+igood))
1453 pam-fi 1.6 badcly=badcly*ibad
1454 mocchiut 1.1 enddo
1455 pam-fi 1.4 * ----------------------------------------------------
1456     * >>> eliminato il taglio sulle BAD <<<
1457     * ----------------------------------------------------
1458 mocchiut 1.1 c if(badcl.eq.0)then
1459     c cl_single(icy)=0
1460     c goto 20
1461     c endif
1462     * ----------------------------------------------------
1463    
1464     cl_good(icy)=1
1465     nply=npl(VIEW(icy))
1466     nldy=nld(MAXS(icy),VIEW(icy))
1467    
1468     * ----------------------------------------------
1469     * CONDITION TO FORM A COUPLE
1470     * ----------------------------------------------
1471     * geometrical consistency (same plane and ladder)
1472     if(nply.eq.nplx.and.nldy.eq.nldx)then
1473     * charge correlation
1474 pam-fi 1.4 * (modified to be applied only below saturation... obviously)
1475    
1476 pam-fi 1.6 if( .not.(dedx(icy).gt.chsaty.and.dedx(icx).gt.chsatx)
1477     $ .and.
1478     $ .not.(dedx(icy).lt.chmipy.and.dedx(icx).lt.chmipx)
1479     $ .and.
1480     $ (badclx.eq.1.and.badcly.eq.1)
1481     $ .and.
1482     $ .true.)then
1483    
1484     ddd=(dedx(icy)
1485     $ -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))
1486     ddd=ddd/sqrt(kch(nplx,nldx)**2+1)
1487    
1488     c cut = chcut * sch(nplx,nldx)
1489    
1490     sss=(kch(nplx,nldx)*dedx(icy)+dedx(icx)
1491     $ -kch(nplx,nldx)*cch(nplx,nldx))
1492     sss=sss/sqrt(kch(nplx,nldx)**2+1)
1493     cut = chcut * (16 + sss/50.)
1494    
1495     if(abs(ddd).gt.cut)then
1496     goto 20 !charge not consistent
1497     endif
1498     endif
1499 pam-fi 1.14
1500 pam-fi 1.15 if(ncp_plane(nplx).gt.ncouplemax)then
1501 pam-fi 1.5 if(verbose)print*,
1502     $ '** warning ** number of identified '//
1503     $ 'couples on plane ',nplx,
1504     $ 'exceeds vector dimention '
1505 pam-fi 1.14 $ ,'( ',ncouplemax,' ) --> masked!'
1506 pam-fi 1.10 mask_view(nviewx(nplx)) = 2
1507     mask_view(nviewy(nply)) = 2
1508 pam-fi 1.15 goto 10
1509 pam-fi 1.5 endif
1510 pam-fi 1.15
1511     * ------------------> COUPLE <------------------
1512     ncp_plane(nplx) = ncp_plane(nplx) + 1
1513     clx(nplx,ncp_plane(nplx))=icx
1514     cly(nply,ncp_plane(nplx))=icy
1515     cl_single(icx)=0
1516     cl_single(icy)=0
1517 pam-fi 1.14 * ----------------------------------------------
1518    
1519 mocchiut 1.1 endif
1520    
1521     20 continue
1522     enddo !end loop on clusters(Y)
1523    
1524     10 continue
1525     enddo !end loop on clusters(X)
1526    
1527    
1528     do icl=1,nclstr1
1529     if(cl_single(icl).eq.1)then
1530     ip=npl(VIEW(icl))
1531     ncls(ip)=ncls(ip)+1
1532     cls(ip,ncls(ip))=icl
1533     endif
1534     enddo
1535    
1536    
1537     if(DEBUG)then
1538     print*,'clusters ',nclstr1
1539     print*,'good ',(cl_good(i),i=1,nclstr1)
1540     print*,'singles ',(cl_single(i),i=1,nclstr1)
1541     print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)
1542     endif
1543    
1544     do ip=1,6
1545 pam-fi 1.10 ncp_tot = ncp_tot + ncp_plane(ip)
1546 mocchiut 1.1 enddo
1547    
1548     return
1549     end
1550    
1551     ***************************************************
1552     * *
1553     * *
1554     * *
1555     * *
1556     * *
1557     * *
1558     **************************************************
1559    
1560     subroutine cp_to_doubtrip(iflag)
1561    
1562     include 'commontracker.f'
1563 pam-fi 1.9 include 'level1.f'
1564 mocchiut 1.1 include 'common_momanhough.f'
1565     include 'common_xyzPAM.f'
1566     include 'common_mini_2.f'
1567     include 'calib.f'
1568    
1569    
1570     * output flag
1571     * --------------
1572     * 0 = good event
1573     * 1 = bad event
1574     * --------------
1575     integer iflag
1576    
1577    
1578     * -----------------------------
1579     * DOUBLETS/TRIPLETS coordinates
1580     c double precision xm1,ym1,zm1
1581     c double precision xm2,ym2,zm2
1582     c double precision xm3,ym3,zm3
1583    
1584     real xm1,ym1,zm1
1585     real xm2,ym2,zm2
1586     real xm3,ym3,zm3
1587     * -----------------------------
1588     * variable needed for tricircle:
1589     real xp(3),zp(3)!TRIPLETS coordinates, to find a circle
1590     EQUIVALENCE (xm1,xp(1))
1591     EQUIVALENCE (xm2,xp(2))
1592     EQUIVALENCE (xm3,xp(3))
1593     EQUIVALENCE (zm1,zp(1))
1594     EQUIVALENCE (zm2,zp(2))
1595     EQUIVALENCE (zm3,zp(3))
1596     real angp(3),resp(3),chi
1597     real xc,zc,radius
1598     * -----------------------------
1599    
1600    
1601 pam-fi 1.14 * --------------------------------------------
1602     * put a limit to the maximum number of couples
1603     * per plane, in order to apply hough transform
1604     * (couples recovered during track refinement)
1605     * --------------------------------------------
1606     do ip=1,nplanes
1607     if(ncp_plane(ip).gt.ncouplelimit)then
1608     mask_view(nviewx(ip)) = 8
1609     mask_view(nviewy(ip)) = 8
1610     endif
1611     enddo
1612    
1613 mocchiut 1.1
1614     ndblt=0 !number of doublets
1615     ntrpt=0 !number of triplets
1616    
1617     do ip1=1,(nplanes-1) !loop on planes - COPPIA 1
1618 pam-fi 1.14 if( mask_view(nviewx(ip1)).ne.0 .or.
1619     $ mask_view(nviewy(ip1)).ne.0 )goto 10 !skip plane
1620     do is1=1,2 !loop on sensors - COPPIA 1
1621 mocchiut 1.1 do icp1=1,ncp_plane(ip1) !loop on COPPIA 1
1622     icx1=clx(ip1,icp1)
1623     icy1=cly(ip1,icp1)
1624     c call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)
1625     call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)
1626     xm1=xPAM
1627     ym1=yPAM
1628     zm1=zPAM
1629     c print*,'***',is1,xm1,ym1,zm1
1630 pam-fi 1.14
1631 mocchiut 1.1 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2
1632 pam-fi 1.14 if( mask_view(nviewx(ip2)).ne.0 .or.
1633     $ mask_view(nviewy(ip2)).ne.0 )goto 20 !skip plane
1634 mocchiut 1.1 do is2=1,2 !loop on sensors -ndblt COPPIA 2
1635    
1636     do icp2=1,ncp_plane(ip2) !loop on COPPIA 2
1637     icx2=clx(ip2,icp2)
1638     icy2=cly(ip2,icp2)
1639     c call xyz_PAM
1640     c $ (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1)
1641     call xyz_PAM
1642     $ (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1)
1643     xm2=xPAM
1644     ym2=yPAM
1645     zm2=zPAM
1646    
1647     * - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1648     * track parameters on Y VIEW
1649     * (2 couples needed)
1650     * - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1651     if(ndblt.eq.ndblt_max)then
1652 pam-fi 1.5 if(verbose)print*,
1653 mocchiut 1.1 $ '** warning ** number of identified '//
1654     $ 'doublets exceeds vector dimention '
1655     $ ,'( ',ndblt_max,' )'
1656     c good2=.false.
1657     c goto 880 !fill ntp and go to next event
1658 pam-fi 1.10 do iv=1,12
1659     mask_view(iv) = 3
1660     enddo
1661 mocchiut 1.1 iflag=1
1662     return
1663     endif
1664     ndblt = ndblt + 1
1665     * store doublet info
1666     cpyz1(ndblt)=id_cp(ip1,icp1,is1)
1667     cpyz2(ndblt)=id_cp(ip2,icp2,is2)
1668     * tg(th_yz)
1669     alfayz2(ndblt)=(ym1-ym2)/(zm1-zm2)
1670     * y0 (cm)
1671     alfayz1(ndblt)=alfayz2(ndblt)*(zini-zm1)+ym1
1672    
1673     **** -----------------------------------------------****
1674     **** reject non phisical couples ****
1675     **** -----------------------------------------------****
1676     if(
1677     $ abs(alfayz2(ndblt)).gt.alfyz2_max
1678     $ .or.
1679     $ abs(alfayz1(ndblt)).gt.alfyz1_max
1680     $ )ndblt = ndblt-1
1681    
1682     c$$$ if(iev.eq.33)then
1683     c$$$ print*,'********* ',ndblt,' -- ',icp1,icp2,is1,is2
1684     c$$$ $ ,' || ',icx1,icy1,icx2,icy2
1685     c$$$ $ ,' || ',xm1,ym1,xm2,ym2
1686     c$$$ $ ,' || ',alfayz2(ndblt),alfayz1(ndblt)
1687     c$$$ endif
1688     c$$$
1689     * - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1690     * track parameters on Y VIEW - end
1691     * - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1692    
1693    
1694 pam-fi 1.14 if(ip2.eq.nplanes)goto 31 !no possible combination with 3 couples
1695    
1696 mocchiut 1.1 do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3
1697 pam-fi 1.14 if( mask_view(nviewx(ip3)).ne.0 .or.
1698     $ mask_view(nviewy(ip3)).ne.0 )goto 30 !skip plane
1699 mocchiut 1.1 do is3=1,2 !loop on sensors - COPPIA 3
1700    
1701     do icp3=1,ncp_plane(ip3) !loop on COPPIA 3
1702     icx3=clx(ip3,icp3)
1703     icy3=cly(ip3,icp3)
1704     c call xyz_PAM
1705     c $ (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1)
1706     call xyz_PAM
1707     $ (icx3,icy3,is3,PFAdef,PFAdef,0.,0.) !(1)
1708     xm3=xPAM
1709     ym3=yPAM
1710     zm3=zPAM
1711     * find the circle passing through the three points
1712     call tricircle(3,xp,zp,angp,resp,chi
1713     $ ,xc,zc,radius,iflag)
1714     c print*,xc,zc,radius
1715     * the circle must intersect the reference plane
1716     if(
1717     c $ (xc.le.-1.*xclimit.or.
1718     c $ xc.ge.xclimit).and.
1719     $ radius**2.ge.(ZINI-zc)**2.and.
1720     $ iflag.eq.0.and.
1721     $ .true.)then
1722    
1723     * - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1724     * track parameters on X VIEW
1725     * (3 couples needed)
1726     * - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1727     if(ntrpt.eq.ntrpt_max)then
1728 pam-fi 1.5 if(verbose)print*,
1729 mocchiut 1.1 $ '** warning ** number of identified '//
1730     $ 'triplets exceeds vector dimention '
1731     $ ,'( ',ntrpt_max,' )'
1732     c good2=.false.
1733     c goto 880 !fill ntp and go to next event
1734 pam-fi 1.10 do iv=1,nviews
1735     mask_view(iv) = 4
1736     enddo
1737 mocchiut 1.1 iflag=1
1738     return
1739     endif
1740     ntrpt = ntrpt +1
1741     * store triplet info
1742     cpxz1(ntrpt)=id_cp(ip1,icp1,is1)
1743     cpxz2(ntrpt)=id_cp(ip2,icp2,is2)
1744     cpxz3(ntrpt)=id_cp(ip3,icp3,is3)
1745    
1746     if(xc.lt.0)then
1747     *************POSITIVE DEFLECTION
1748     alfaxz1(ntrpt) = xc+sqrt(radius**2-(ZINI-zc)**2)
1749     alfaxz2(ntrpt) = (ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2)
1750     alfaxz3(ntrpt) = 1/radius
1751     else
1752     *************NEGATIVE DEFLECTION
1753     alfaxz1(ntrpt) = xc-sqrt(radius**2-(ZINI-zc)**2)
1754     alfaxz2(ntrpt) = -(ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2)
1755     alfaxz3(ntrpt) = -1/radius
1756     endif
1757    
1758     **** -----------------------------------------------****
1759     **** reject non phisical triplets ****
1760     **** -----------------------------------------------****
1761     if(
1762     $ abs(alfaxz2(ntrpt)).gt.alfxz2_max
1763     $ .or.
1764     $ abs(alfaxz1(ntrpt)).gt.alfxz1_max
1765     $ )ntrpt = ntrpt-1
1766    
1767    
1768     c print*,alfaxz1(ntrpt),alfaxz2(ntrpt),alfaxz3(ntrpt)
1769     * - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1770     * track parameters on X VIEW - end
1771     * - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1772     endif
1773     enddo !end loop on COPPIA 3
1774     enddo !end loop on sensors - COPPIA 3
1775 pam-fi 1.14 30 continue
1776 mocchiut 1.1 enddo !end loop on planes - COPPIA 3
1777 pam-fi 1.14 31 continue
1778 mocchiut 1.1
1779     1 enddo !end loop on COPPIA 2
1780     enddo !end loop on sensors - COPPIA 2
1781 pam-fi 1.14 20 continue
1782 mocchiut 1.1 enddo !end loop on planes - COPPIA 2
1783    
1784     enddo !end loop on COPPIA1
1785     enddo !end loop on sensors - COPPIA 1
1786 pam-fi 1.14 10 continue
1787 mocchiut 1.1 enddo !end loop on planes - COPPIA 1
1788    
1789     if(DEBUG)then
1790     print*,'--- doublets ',ndblt
1791     print*,'--- triplets ',ntrpt
1792     endif
1793    
1794     c goto 880 !ntp fill
1795    
1796    
1797     return
1798     end
1799    
1800    
1801    
1802     ***************************************************
1803     * *
1804     * *
1805     * *
1806     * *
1807     * *
1808     * *
1809     **************************************************
1810    
1811     subroutine doub_to_YZcloud(iflag)
1812    
1813     include 'commontracker.f'
1814 pam-fi 1.9 include 'level1.f'
1815 mocchiut 1.1 include 'common_momanhough.f'
1816 pam-fi 1.9 c include 'momanhough_init.f'
1817 mocchiut 1.1
1818    
1819     * output flag
1820     * --------------
1821     * 0 = good event
1822     * 1 = bad event
1823     * --------------
1824     integer iflag
1825    
1826     integer db_used(ndblt_max)
1827     integer db_temp(ndblt_max)
1828     integer db_all(ndblt_max) !stores db ID in each cloud
1829    
1830     integer hit_plane(nplanes)
1831    
1832     * mask for used couples
1833     integer cp_useds1(ncouplemaxtot) ! sensor 1
1834     integer cp_useds2(ncouplemaxtot) ! sensor 2
1835    
1836    
1837     *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1838     * classification of DOUBLETS
1839     * according to distance in parameter space
1840     * (cloud = group of points (doublets) in parameter space)
1841     *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1842     do idb=1,ndblt
1843     db_used(idb)=0
1844     enddo
1845    
1846     distance=0
1847     nclouds_yz=0 !number of clouds
1848     npt_tot=0
1849 pam-fi 1.9 nloop=0
1850     90 continue
1851 mocchiut 1.1 do idb1=1,ndblt !loop (1) on DOUBLETS
1852     if(db_used(idb1).eq.1)goto 2228 !db already included in a cloud
1853    
1854     c print*,'--------------'
1855     c print*,'** ',idb1,' **'
1856    
1857     do icp=1,ncp_tot
1858     cp_useds1(icp)=0 !init
1859     cp_useds2(icp)=0 !init
1860     enddo
1861     do idb=1,ndblt
1862     db_all(idb)=0
1863     enddo
1864     if(cpyz1(idb1).gt.0)cp_useds2(cpyz1(idb1))=1
1865     if(cpyz1(idb1).lt.0)cp_useds1(-cpyz1(idb1))=1
1866     if(cpyz2(idb1).gt.0)cp_useds2(cpyz2(idb1))=1
1867     if(cpyz2(idb1).lt.0)cp_useds1(-cpyz2(idb1))=1
1868     temp1 = alfayz1(idb1)
1869     temp2 = alfayz2(idb1)
1870     npt=1 !counter of points in the cloud
1871    
1872     db_all(npt) = idb1
1873    
1874     nptloop=1
1875     db_temp(1)=idb1
1876    
1877     88 continue
1878    
1879     npv=0 !# new points inlcuded
1880     do iloop=1,nptloop
1881     idbref=db_temp(iloop) !local point of reference
1882     ccccc if(db_used(idbref).eq.1)goto 1188 !next
1883    
1884     do idb2=1,ndblt !loop (2) on DOUBLETS
1885     if(idb2.eq.idbref)goto 1118 !next doublet
1886     if(db_used(idb2).eq.1)goto 1118
1887    
1888    
1889     * doublet distance in parameter space
1890     distance=
1891     $ ((alfayz1(idbref)-alfayz1(idb2))/Dalfayz1)**2
1892     $ +((alfayz2(idbref)-alfayz2(idb2))/Dalfayz2)**2
1893     distance = sqrt(distance)
1894    
1895     c$$$ if(iev.eq.33)then
1896     c$$$ if(distance.lt.100)
1897     c$$$ $ print*,'********* ',idb1,idbref,idb2,distance
1898     c$$$ if(distance.lt.100)
1899     c$$$ $ print*,'********* ',alfayz1(idbref),alfayz1(idb2)
1900     c$$$ $ ,alfayz2(idbref),alfayz2(idb2)
1901     c$$$ endif
1902     if(distance.lt.cutdistyz)then
1903    
1904     c print*,idb1,idb2,distance,' cloud ',nclouds_yz
1905     if(cpyz1(idb2).gt.0)cp_useds2(cpyz1(idb2))=1
1906     if(cpyz1(idb2).lt.0)cp_useds1(-cpyz1(idb2))=1
1907     if(cpyz2(idb2).gt.0)cp_useds2(cpyz2(idb2))=1
1908     if(cpyz2(idb2).lt.0)cp_useds1(-cpyz2(idb2))=1
1909     npt = npt + 1 !counter of points in the cloud
1910    
1911     npv = npv +1
1912     db_temp(npv) = idb2
1913     db_used(idbref) = 1
1914     db_used(idb2) = 1
1915    
1916     db_all(npt) = idb2
1917    
1918     temp1 = temp1 + alfayz1(idb2)
1919     temp2 = temp2 + alfayz2(idb2)
1920     c print*,'* idbref,idb2 ',idbref,idb2
1921     endif
1922    
1923     1118 continue
1924     enddo !end loop (2) on DOUBLETS
1925    
1926     1188 continue
1927     enddo !end loop on... bo?
1928    
1929     nptloop=npv
1930     if(nptloop.ne.0)goto 88
1931    
1932     * ------------------------------------------
1933     * stores the cloud only if
1934     * 1) it includes a minimum number of REAL couples
1935     * 1bis) it inlcudes a minimum number of doublets
1936     * 2) it is not already stored
1937     * ------------------------------------------
1938     do ip=1,nplanes
1939     hit_plane(ip)=0
1940     enddo
1941     ncpused=0
1942     do icp=1,ncp_tot
1943     if(cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0)then
1944     ncpused=ncpused+1
1945     ip=ip_cp(icp)
1946     hit_plane(ip)=1
1947     endif
1948     enddo
1949     nplused=0
1950     do ip=1,nplanes
1951     nplused=nplused+ hit_plane(ip)
1952     enddo
1953     c print*,'>>>> ',ncpused,npt,nplused
1954 pam-fi 1.9 c if(ncpused.lt.ncpyz_min)goto 2228 !next doublet
1955 mocchiut 1.1 if(npt.lt.nptyz_min)goto 2228 !next doublet
1956     if(nplused.lt.nplyz_min)goto 2228 !next doublet
1957    
1958     * ~~~~~~~~~~~~~~~~~
1959     * >>> NEW CLOUD <<<
1960    
1961     if(nclouds_yz.ge.ncloyz_max)then
1962 pam-fi 1.5 if(verbose)print*,
1963 mocchiut 1.1 $ '** warning ** number of identified '//
1964     $ 'YZ clouds exceeds vector dimention '
1965     $ ,'( ',ncloyz_max,' )'
1966     c good2=.false.
1967     c goto 880 !fill ntp and go to next event
1968 pam-fi 1.10 do iv=1,nviews
1969     mask_view(iv) = 5
1970     enddo
1971 mocchiut 1.1 iflag=1
1972     return
1973     endif
1974    
1975     nclouds_yz = nclouds_yz + 1 !increase counter
1976     alfayz1_av(nclouds_yz) = temp1/npt !store average parameter
1977     alfayz2_av(nclouds_yz) = temp2/npt ! "
1978     do icp=1,ncp_tot
1979     cpcloud_yz(nclouds_yz,icp)=
1980     $ cp_useds1(icp)+2*cp_useds2(icp) !store cp info
1981     enddo
1982     ptcloud_yz(nclouds_yz)=npt
1983     c ptcloud_yz_nt(nclouds_yz)=npt
1984     do ipt=1,npt
1985     db_cloud(npt_tot+ipt) = db_all(ipt)
1986     c print*,'>> ',ipt,db_all(ipt)
1987     enddo
1988     npt_tot=npt_tot+npt
1989     if(DEBUG)then
1990     print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'
1991     print*,'>>>> cloud ',nclouds_yz,' --- ',npt,' points'
1992     print*,'- alfayz1 ',alfayz1_av(nclouds_yz)
1993     print*,'- alfayz2 ',alfayz2_av(nclouds_yz)
1994     print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot)
1995     print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot)
1996     print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)
1997     c$$$ print*,'nt-uple: ptcloud_yz(',nclouds_yz,') = '
1998     c$$$ $ ,ptcloud_yz(nclouds_yz)
1999     c$$$ print*,'nt-uple: db_cloud(...) = '
2000     c$$$ $ ,(db_cloud(iii),iii=npt_tot-npt+1,npt_tot)
2001     endif
2002     * >>> NEW CLOUD <<<
2003     * ~~~~~~~~~~~~~~~~~
2004     2228 continue
2005     enddo !end loop (1) on DOUBLETS
2006    
2007    
2008 pam-fi 1.9 if(nloop.lt.nstepy)then
2009     cutdistyz = cutdistyz+cutystep
2010     nloop = nloop+1
2011     goto 90
2012     endif
2013    
2014 mocchiut 1.1 if(DEBUG)then
2015     print*,'---------------------- '
2016     print*,'Y-Z total clouds ',nclouds_yz
2017     print*,' '
2018     endif
2019    
2020    
2021     return
2022     end
2023    
2024    
2025    
2026    
2027    
2028     ***************************************************
2029     * *
2030     * *
2031     * *
2032     * *
2033     * *
2034     * *
2035     **************************************************
2036    
2037     subroutine trip_to_XZcloud(iflag)
2038    
2039     include 'commontracker.f'
2040 pam-fi 1.9 include 'level1.f'
2041 mocchiut 1.1 include 'common_momanhough.f'
2042 pam-fi 1.9 c include 'momanhough_init.f'
2043 mocchiut 1.1
2044    
2045     * output flag
2046     * --------------
2047     * 0 = good event
2048     * 1 = bad event
2049     * --------------
2050     integer iflag
2051    
2052     integer tr_used(ntrpt_max)
2053     integer tr_temp(ntrpt_max)
2054     integer tr_incl(ntrpt_max)
2055     integer tr_all(ntrpt_max) !stores tr ID in each cloud
2056    
2057     integer hit_plane(nplanes)
2058    
2059     * mask for used couples
2060     integer cp_useds1(ncouplemaxtot) ! sensor 1
2061     integer cp_useds2(ncouplemaxtot) ! sensor 2
2062    
2063     *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2064     * classification of TRIPLETS
2065     * according to distance in parameter space
2066     *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2067     do itr=1,ntrpt
2068     tr_used(itr)=0
2069     enddo
2070    
2071     distance=0
2072     nclouds_xz=0 !number of clouds
2073     npt_tot=0 !total number of selected triplets
2074 pam-fi 1.9 nloop=0
2075     91 continue
2076 mocchiut 1.1 do itr1=1,ntrpt !loop (1) on TRIPLETS
2077     if(tr_used(itr1).eq.1)goto 22288 !already included in a cloud
2078     c print*,'--------------'
2079     c print*,'** ',itr1,' **'
2080    
2081     do icp=1,ncp_tot
2082     cp_useds1(icp)=0
2083     cp_useds2(icp)=0
2084     enddo
2085     do itr=1,ntrpt
2086     tr_all(itr)=0 !list of included triplets
2087     enddo
2088     if(cpxz1(itr1).gt.0)cp_useds2(cpxz1(itr1))=1
2089     if(cpxz1(itr1).lt.0)cp_useds1(-cpxz1(itr1))=1
2090     if(cpxz2(itr1).gt.0)cp_useds2(cpxz2(itr1))=1
2091     if(cpxz2(itr1).lt.0)cp_useds1(-cpxz2(itr1))=1
2092     if(cpxz3(itr1).gt.0)cp_useds2(cpxz3(itr1))=1
2093     if(cpxz3(itr1).lt.0)cp_useds1(-cpxz3(itr1))=1
2094     temp1 = alfaxz1(itr1)
2095     temp2 = alfaxz2(itr1)
2096     temp3 = alfaxz3(itr1)
2097     npt=1 !counter of points in the cloud
2098    
2099     tr_all(npt) = itr1
2100    
2101     nptloop=1
2102     c tr_temp(1)=itr1
2103     tr_incl(1)=itr1
2104    
2105     8881 continue
2106    
2107     npv=0 !# new points inlcuded
2108     do iloop=1,nptloop
2109     itrref=tr_incl(iloop) !local point of reference
2110     do itr2=1,ntrpt !loop (2) on TRIPLETS
2111     if(itr2.eq.itr1)goto 11188 !next triplet
2112     if(tr_used(itr2).eq.1)goto 11188 !next triplet
2113     * triplet distance in parameter space
2114     * solo i due parametri spaziali per il momemnto
2115     distance=
2116     $ ((alfaxz1(itrref)-alfaxz1(itr2))/Dalfaxz1)**2
2117     $ +((alfaxz2(itrref)-alfaxz2(itr2))/Dalfaxz2)**2
2118     distance = sqrt(distance)
2119    
2120     if(distance.lt.cutdistxz)then
2121     c print*,idb1,idb2,distance,' cloud ',nclouds_yz
2122     if(cpxz1(itr2).gt.0)cp_useds2(cpxz1(itr2))=1
2123     if(cpxz1(itr2).lt.0)cp_useds1(-cpxz1(itr2))=1
2124     if(cpxz2(itr2).gt.0)cp_useds2(cpxz2(itr2))=1
2125     if(cpxz2(itr2).lt.0)cp_useds1(-cpxz2(itr2))=1
2126     if(cpxz3(itr2).gt.0)cp_useds2(cpxz3(itr2))=1
2127     if(cpxz3(itr2).lt.0)cp_useds1(-cpxz3(itr2))=1
2128     npt = npt + 1 !counter of points in the cloud
2129    
2130     npv = npv +1
2131     tr_temp(npv) = itr2
2132     tr_used(itrref) = 1
2133     tr_used(itr2) = 1
2134    
2135     tr_all(npt) = itr2
2136    
2137     temp1 = temp1 + alfaxz1(itr2)
2138     temp2 = temp2 + alfaxz2(itr2)
2139     temp3 = temp3 + alfaxz3(itr2)
2140     c print*,'* itrref,itr2 ',itrref,itr2,distance
2141     endif
2142    
2143     11188 continue
2144     enddo !end loop (2) on TRIPLETS
2145    
2146     11888 continue
2147     enddo !end loop on... bo?
2148    
2149     nptloop=npv
2150     do i=1,npv
2151     tr_incl(i)=tr_temp(i)
2152     enddo
2153     if(nptloop.ne.0)goto 8881
2154    
2155     * ------------------------------------------
2156     * stores the cloud only if
2157     * 1) it includes a minimum number of REAL couples
2158     * 1bis)
2159     * 2) it is not already stored
2160     * ------------------------------------------
2161     c print*,'check cp_used'
2162     do ip=1,nplanes
2163     hit_plane(ip)=0
2164     enddo
2165     ncpused=0
2166     do icp=1,ncp_tot
2167     if(cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0)then
2168     ncpused=ncpused+1
2169     ip=ip_cp(icp)
2170     hit_plane(ip)=1
2171     endif
2172     enddo
2173     nplused=0
2174     do ip=1,nplanes
2175     nplused=nplused+ hit_plane(ip)
2176     enddo
2177 pam-fi 1.9 c if(ncpused.lt.ncpxz_min)goto 22288 !next triplet
2178 mocchiut 1.1 if(npt.lt.nptxz_min)goto 22288 !next triplet
2179     if(nplused.lt.nplxz_min)goto 22288 !next doublet
2180    
2181     * ~~~~~~~~~~~~~~~~~
2182     * >>> NEW CLOUD <<<
2183     if(nclouds_xz.ge.ncloxz_max)then
2184 pam-fi 1.5 if(verbose)print*,
2185 mocchiut 1.1 $ '** warning ** number of identified '//
2186     $ 'XZ clouds exceeds vector dimention '
2187     $ ,'( ',ncloxz_max,' )'
2188     c good2=.false.
2189     c goto 880 !fill ntp and go to next event
2190 pam-fi 1.10 do iv=1,nviews
2191     mask_view(iv) = 6
2192     enddo
2193 mocchiut 1.1 iflag=1
2194     return
2195     endif
2196     nclouds_xz = nclouds_xz + 1 !increase counter
2197     alfaxz1_av(nclouds_xz) = temp1/npt !store average parameter
2198     alfaxz2_av(nclouds_xz) = temp2/npt ! "
2199     alfaxz3_av(nclouds_xz) = temp3/npt ! "
2200     do icp=1,ncp_tot
2201     cpcloud_xz(nclouds_xz,icp)=
2202     $ cp_useds1(icp)+2*cp_useds2(icp) !store cp info
2203     enddo
2204     ptcloud_xz(nclouds_xz)=npt
2205     do ipt=1,npt
2206     tr_cloud(npt_tot+ipt) = tr_all(ipt)
2207     enddo
2208     npt_tot=npt_tot+npt
2209    
2210     if(DEBUG)then
2211     print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'
2212     print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points'
2213     print*,'- alfaxz1 ',alfaxz1_av(nclouds_xz)
2214     print*,'- alfaxz2 ',alfaxz2_av(nclouds_xz)
2215     print*,'- alfaxz3 ',alfaxz3_av(nclouds_xz)
2216     print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot)
2217     print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot)
2218     print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)
2219     c$$$ print*,'nt-uple: ptcloud_xz(',nclouds_xz,') = '
2220     c$$$ $ ,ptcloud_xz(nclouds_xz)
2221     c$$$ print*,'nt-uple: tr_cloud(...) = '
2222     c$$$ $ ,(tr_cloud(iii),iii=npt_tot-npt+1,npt_tot)
2223     endif
2224     * >>> NEW CLOUD <<<
2225     * ~~~~~~~~~~~~~~~~~
2226     22288 continue
2227     enddo !end loop (1) on DOUBLETS
2228 pam-fi 1.9
2229     if(nloop.lt.nstepx)then
2230     cutdistxz=cutdistxz+cutxstep
2231     nloop=nloop+1
2232     goto 91
2233     endif
2234    
2235 mocchiut 1.1 if(DEBUG)then
2236     print*,'---------------------- '
2237     print*,'X-Z total clouds ',nclouds_xz
2238     print*,' '
2239     endif
2240    
2241    
2242     return
2243     end
2244    
2245    
2246     ***************************************************
2247     * *
2248     * *
2249     * *
2250     * *
2251     * *
2252     * *
2253     **************************************************
2254    
2255     subroutine clouds_to_ctrack(iflag)
2256     c*****************************************************
2257     c 02/02/2006 modified by Elena Vannuccini --> (1)
2258     c*****************************************************
2259    
2260     include 'commontracker.f'
2261 pam-fi 1.9 include 'level1.f'
2262 mocchiut 1.1 include 'common_momanhough.f'
2263     include 'common_xyzPAM.f'
2264     include 'common_mini_2.f'
2265     include 'common_mech.f'
2266 pam-fi 1.9 c include 'momanhough_init.f'
2267 mocchiut 1.1
2268    
2269     * output flag
2270     * --------------
2271     * 0 = good event
2272     * 1 = bad event
2273     * --------------
2274     integer iflag
2275    
2276     * -----------------------------------------------------------
2277     * mask to store (locally) the couples included
2278     * in the intersection bewteen a XZ and YZ cloud
2279     integer cpintersec(ncouplemaxtot)
2280     * -----------------------------------------------------------
2281     * list of matching couples in the combination
2282     * between a XZ and YZ cloud
2283 pam-fi 1.9 integer cp_match(nplanes,2*ncouplemax)
2284 mocchiut 1.1 integer ncp_match(nplanes)
2285     * -----------------------------------------------------------
2286     integer hit_plane(nplanes)
2287     * -----------------------------------------------------------
2288     * variables for track fitting
2289     double precision AL_INI(5)
2290 pam-fi 1.11 c double precision tath
2291 mocchiut 1.1 * -----------------------------------------------------------
2292     c real fitz(nplanes) !z coordinates of the planes in cm
2293    
2294    
2295    
2296     ntracks=0 !counter of track candidates
2297    
2298     do iyz=1,nclouds_yz !loop on YZ couds
2299     do ixz=1,nclouds_xz !loop on XZ couds
2300    
2301     * --------------------------------------------------
2302     * check of consistency of the clouds
2303     * ---> required a minimum number of matching couples
2304     * the track fit will be performed on the INTERSECTION
2305     * of the two clouds
2306     * --------------------------------------------------
2307     do ip=1,nplanes
2308     hit_plane(ip)=0
2309     ncp_match(ip)=0
2310     do icpp=1,ncouplemax
2311     cp_match(ip,icpp)=0 !init couple list
2312     enddo
2313     enddo
2314     ncp_ok=0
2315     do icp=1,ncp_tot !loop on couples
2316     * get info on
2317     cpintersec(icp)=min(
2318     $ cpcloud_yz(iyz,icp),
2319     $ cpcloud_xz(ixz,icp))
2320     if(
2321     $ (cpcloud_yz(iyz,icp).eq.1.and.cpcloud_xz(ixz,icp).eq.2).or.
2322     $ (cpcloud_yz(iyz,icp).eq.2.and.cpcloud_xz(ixz,icp).eq.1).or.
2323     $ .false.)cpintersec(icp)=0
2324     if(cpintersec(icp).ne.0)then
2325     ncp_ok=ncp_ok+1
2326    
2327     ip=ip_cp(icp)
2328     hit_plane(ip)=1
2329     if(cpintersec(icp).eq.1)then
2330     * 1) only the couple image in sensor 1 matches
2331     id=-icp
2332     ncp_match(ip)=ncp_match(ip)+1
2333     cp_match(ip,ncp_match(ip))=id
2334     elseif(cpintersec(icp).eq.2)then
2335     * 2) only the couple image in sensor 2 matches
2336     id=icp
2337     ncp_match(ip)=ncp_match(ip)+1
2338     cp_match(ip,ncp_match(ip))=id
2339     else
2340     * 3) both couple images match
2341     id=icp
2342     do is=1,2
2343     id=-id
2344     ncp_match(ip)=ncp_match(ip)+1
2345     cp_match(ip,ncp_match(ip))=id
2346     enddo
2347     endif
2348     endif !end matching condition
2349     enddo !end loop on couples
2350    
2351     nplused=0
2352     do ip=1,nplanes
2353     nplused=nplused+ hit_plane(ip)
2354     enddo
2355    
2356 pam-fi 1.15 c if(nplused.lt.nplxz_min)goto 888 !next doublet
2357     if(nplused.lt.nplyz_min)goto 888 !next doublet
2358 mocchiut 1.1 if(ncp_ok.lt.ncpok)goto 888 !next cloud
2359    
2360     if(DEBUG)then
2361     print*,'Combination ',iyz,ixz
2362     $ ,' db ',ptcloud_yz(iyz)
2363     $ ,' tr ',ptcloud_xz(ixz)
2364     $ ,' -----> # matching couples ',ncp_ok
2365     endif
2366     c$$$ print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'
2367     c$$$ print*,'Configurazione cluster XZ'
2368     c$$$ print*,'1 -- ',(clx(1,i),i=1,ncp_plane(1))
2369     c$$$ print*,'2 -- ',(clx(2,i),i=1,ncp_plane(1))
2370     c$$$ print*,'3 -- ',(clx(3,i),i=1,ncp_plane(1))
2371     c$$$ print*,'4 -- ',(clx(4,i),i=1,ncp_plane(1))
2372     c$$$ print*,'5 -- ',(clx(5,i),i=1,ncp_plane(1))
2373     c$$$ print*,'6 -- ',(clx(6,i),i=1,ncp_plane(1))
2374     c$$$ print*,'Configurazione cluster YZ'
2375     c$$$ print*,'1 -- ',(cly(1,i),i=1,ncp_plane(1))
2376     c$$$ print*,'2 -- ',(cly(2,i),i=1,ncp_plane(1))
2377     c$$$ print*,'3 -- ',(cly(3,i),i=1,ncp_plane(1))
2378     c$$$ print*,'4 -- ',(cly(4,i),i=1,ncp_plane(1))
2379     c$$$ print*,'5 -- ',(cly(5,i),i=1,ncp_plane(1))
2380     c$$$ print*,'6 -- ',(cly(6,i),i=1,ncp_plane(1))
2381     c$$$ print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'
2382    
2383     * -------> INITIAL GUESS <-------
2384 pam-fi 1.11 cccc SBAGLIATO
2385     c$$$ AL_INI(1) = dreal(alfaxz1_av(ixz))
2386     c$$$ AL_INI(2) = dreal(alfayz1_av(iyz))
2387     c$$$ AL_INI(4) = PIGR + datan(dreal(alfayz2_av(iyz))
2388     c$$$ $ /dreal(alfaxz2_av(ixz)))
2389     c$$$ tath = -dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))
2390     c$$$ AL_INI(3) = tath/sqrt(1+tath**2)
2391     c$$$ AL_INI(5) = (1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.
2392     cccc GIUSTO (ma si sua guess())
2393     c$$$ AL_INI(1) = dreal(alfaxz1_av(ixz))
2394     c$$$ AL_INI(2) = dreal(alfayz1_av(iyz))
2395     c$$$ tath = -dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))
2396     c$$$ AL_INI(3) = tath/sqrt(1+tath**2)
2397     c$$$ IF(alfaxz2_av(ixz).NE.0)THEN
2398     c$$$ AL_INI(4) = PIGR + datan(dreal(alfayz2_av(iyz))
2399     c$$$ $ /dreal(alfaxz2_av(ixz)))
2400     c$$$ ELSE
2401     c$$$ AL_INI(4) = acos(-1.)/2
2402     c$$$ IF(alfayz2_av(iyz).LT.0)AL_INI(4) = AL_INI(4)+acos(-1.)
2403     c$$$ ENDIF
2404     c$$$ IF(alfaxz2_av(ixz).LT.0)AL_INI(4)= acos(-1.)+ AL_INI(4)
2405     c$$$ AL_INI(4) = -acos(-1.) + AL_INI(4) !from incidence direction to tracking rs
2406     c$$$
2407     c$$$ AL_INI(5) = (1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.
2408     c$$$
2409     c$$$ if(AL_INI(5).gt.defmax)goto 888 !next cloud
2410    
2411 mocchiut 1.1 if(DEBUG)then
2412     print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))
2413     print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))
2414     print*,'3 >>> ',(cp_match(4,i),i=1,ncp_match(4))
2415     print*,'4 >>> ',(cp_match(3,i),i=1,ncp_match(3))
2416     print*,'5 >>> ',(cp_match(2,i),i=1,ncp_match(2))
2417     print*,'6 >>> ',(cp_match(1,i),i=1,ncp_match(1))
2418     endif
2419    
2420     do icp1=1,max(1,ncp_match(1))
2421     hit_plane(1)=icp1
2422     if(ncp_match(1).eq.0)hit_plane(1)=0 !-icp1
2423    
2424     do icp2=1,max(1,ncp_match(2))
2425     hit_plane(2)=icp2
2426     if(ncp_match(2).eq.0)hit_plane(2)=0 !-icp2
2427    
2428     do icp3=1,max(1,ncp_match(3))
2429     hit_plane(3)=icp3
2430     if(ncp_match(3).eq.0)hit_plane(3)=0 !-icp3
2431    
2432     do icp4=1,max(1,ncp_match(4))
2433     hit_plane(4)=icp4
2434     if(ncp_match(4).eq.0)hit_plane(4)=0 !-icp4
2435    
2436     do icp5=1,max(1,ncp_match(5))
2437     hit_plane(5)=icp5
2438     if(ncp_match(5).eq.0)hit_plane(5)=0 !-icp5
2439    
2440     do icp6=1,max(1,ncp_match(6))
2441     hit_plane(6)=icp6
2442     if(ncp_match(6).eq.0)hit_plane(6)=0 !-icp6
2443    
2444    
2445     call track_init !init TRACK common
2446    
2447     do ip=1,nplanes !loop on planes
2448     if(hit_plane(ip).ne.0)then
2449     id=cp_match(ip,hit_plane(ip))
2450     is=is_cp(id)
2451     icp=icp_cp(id)
2452     if(ip_cp(id).ne.ip)
2453     $ print*,'OKKIO!!'
2454     $ ,'id ',id,is,icp
2455     $ ,ip_cp(id),ip
2456     icx=clx(ip,icp)
2457     icy=cly(ip,icp)
2458     * *************************
2459     c call xyz_PAM(icx,icy,is,
2460     c $ 'COG2','COG2',0.,0.)
2461     call xyz_PAM(icx,icy,is, !(1)
2462     $ PFAdef,PFAdef,0.,0.) !(1)
2463     * *************************
2464     * -----------------------------
2465     xgood(nplanes-ip+1)=1.
2466     ygood(nplanes-ip+1)=1.
2467     xm(nplanes-ip+1)=xPAM
2468     ym(nplanes-ip+1)=yPAM
2469     zm(nplanes-ip+1)=zPAM
2470     resx(nplanes-ip+1)=resxPAM
2471     resy(nplanes-ip+1)=resyPAM
2472     * -----------------------------
2473     endif
2474     enddo !end loop on planes
2475     * **********************************************************
2476     * ************************** FIT *** FIT *** FIT *** FIT ***
2477     * **********************************************************
2478 pam-fi 1.11 cccc scommentare se si usa al_ini della nuvola
2479     c$$$ do i=1,5
2480     c$$$ AL(i)=AL_INI(i)
2481     c$$$ enddo
2482     call guess()
2483 mocchiut 1.1 do i=1,5
2484 pam-fi 1.11 AL_INI(i)=AL(i)
2485 mocchiut 1.1 enddo
2486     ifail=0 !error flag in chi^2 computation
2487     jstep=0 !number of minimization steps
2488 pam-fi 1.8 iprint=0
2489 pam-fi 1.11 c if(DEBUG)iprint=1
2490 pam-fi 1.8 if(DEBUG)iprint=1
2491     call mini2(jstep,ifail,iprint)
2492 mocchiut 1.1 if(ifail.ne.0) then
2493     if(DEBUG)then
2494     print *,
2495     $ '*** MINIMIZATION FAILURE *** '
2496 pam-fi 1.11 $ //'(clouds_to_ctrack)'
2497     print*,'initial guess: '
2498    
2499     print*,'AL_INI(1) = ',AL_INI(1)
2500     print*,'AL_INI(2) = ',AL_INI(2)
2501     print*,'AL_INI(3) = ',AL_INI(3)
2502     print*,'AL_INI(4) = ',AL_INI(4)
2503     print*,'AL_INI(5) = ',AL_INI(5)
2504 mocchiut 1.1 endif
2505 pam-fi 1.11 c chi2=-chi2
2506 mocchiut 1.1 endif
2507     * **********************************************************
2508     * ************************** FIT *** FIT *** FIT *** FIT ***
2509     * **********************************************************
2510    
2511     if(chi2.le.0.)goto 666
2512    
2513     * --------------------------
2514     * STORE candidate TRACK INFO
2515     * --------------------------
2516     if(ntracks.eq.NTRACKSMAX)then
2517    
2518 pam-fi 1.5 if(verbose)print*,
2519 mocchiut 1.1 $ '** warning ** number of candidate tracks '//
2520     $ ' exceeds vector dimension '
2521     $ ,'( ',NTRACKSMAX,' )'
2522     c good2=.false.
2523     c goto 880 !fill ntp and go to next event
2524 pam-fi 1.10 do iv=1,nviews
2525     mask_view(iv) = 7
2526     enddo
2527 mocchiut 1.1 iflag=1
2528     return
2529     endif
2530    
2531     ntracks = ntracks + 1
2532    
2533     c$$$ ndof=0
2534     do ip=1,nplanes
2535     c$$$ ndof=ndof
2536     c$$$ $ +int(xgood(ip))
2537     c$$$ $ +int(ygood(ip))
2538     XV_STORE(ip,ntracks)=sngl(xv(ip))
2539     YV_STORE(ip,ntracks)=sngl(yv(ip))
2540     ZV_STORE(ip,ntracks)=sngl(zv(ip))
2541     XM_STORE(ip,ntracks)=sngl(xm(ip))
2542     YM_STORE(ip,ntracks)=sngl(ym(ip))
2543     ZM_STORE(ip,ntracks)=sngl(zm(ip))
2544     RESX_STORE(ip,ntracks)=sngl(resx(ip))
2545     RESY_STORE(ip,ntracks)=sngl(resy(ip))
2546     XV_STORE(ip,ntracks)=sngl(xv(ip))
2547     YV_STORE(ip,ntracks)=sngl(yv(ip))
2548     ZV_STORE(ip,ntracks)=sngl(zv(ip))
2549     AXV_STORE(ip,ntracks)=sngl(axv(ip))
2550     AYV_STORE(ip,ntracks)=sngl(ayv(ip))
2551     XGOOD_STORE(ip,ntracks)=sngl(xgood(ip))
2552     YGOOD_STORE(ip,ntracks)=sngl(ygood(ip))
2553     if(hit_plane(ip).ne.0)then
2554     CP_STORE(nplanes-ip+1,ntracks)=
2555     $ cp_match(ip,hit_plane(ip))
2556     else
2557     CP_STORE(nplanes-ip+1,ntracks)=0
2558     endif
2559     CLS_STORE(nplanes-ip+1,ntracks)=0
2560     do i=1,5
2561     AL_STORE(i,ntracks)=sngl(AL(i))
2562     enddo
2563     enddo
2564    
2565     c$$$ * Number of Degree Of Freedom
2566     c$$$ ndof=ndof-5
2567     c$$$ * reduced chi^2
2568     c$$$ rchi2=chi2/dble(ndof)
2569     RCHI2_STORE(ntracks)=chi2
2570    
2571     * --------------------------------
2572     * STORE candidate TRACK INFO - end
2573     * --------------------------------
2574    
2575     666 continue
2576     enddo !end loop on cp in plane 6
2577     enddo !end loop on cp in plane 5
2578     enddo !end loop on cp in plane 4
2579     enddo !end loop on cp in plane 3
2580     enddo !end loop on cp in plane 2
2581     enddo !end loop on cp in plane 1
2582    
2583     888 continue
2584     enddo !end loop on XZ couds
2585     enddo !end loop on YZ couds
2586    
2587     if(ntracks.eq.0)then
2588     iflag=1
2589     return
2590     endif
2591    
2592     if(DEBUG)then
2593     print*,'****** TRACK CANDIDATES ***********'
2594     print*,'# R. chi2 RIG'
2595     do i=1,ntracks
2596     print*,i,' --- ',rchi2_store(i),' --- '
2597     $ ,1./abs(AL_STORE(5,i))
2598     enddo
2599     print*,'***********************************'
2600     endif
2601    
2602    
2603     return
2604     end
2605    
2606    
2607     ***************************************************
2608     * *
2609     * *
2610     * *
2611     * *
2612     * *
2613     * *
2614     **************************************************
2615    
2616     subroutine refine_track(ibest)
2617    
2618     c******************************************************
2619     cccccc 06/10/2005 modified by elena vannuccini ---> (1)
2620     cccccc 31/01/2006 modified by elena vannuccini ---> (2)
2621 pam-fi 1.3 cccccc 12/08/2006 modified by elena vannucicni ---> (3)
2622 mocchiut 1.1 c******************************************************
2623    
2624     include 'commontracker.f'
2625 pam-fi 1.9 include 'level1.f'
2626 mocchiut 1.1 include 'common_momanhough.f'
2627     include 'common_xyzPAM.f'
2628     include 'common_mini_2.f'
2629     include 'common_mech.f'
2630 pam-fi 1.9 c include 'momanhough_init.f'
2631     c include 'level1.f'
2632 mocchiut 1.1 include 'calib.f'
2633    
2634    
2635     * flag to chose PFA
2636     character*10 PFA
2637     common/FINALPFA/PFA
2638    
2639     * =================================================
2640     * new estimate of positions using ETA algorithm
2641     * and
2642     * search for new couples and single clusters to add
2643     * =================================================
2644     call track_init
2645     do ip=1,nplanes !loop on planes
2646    
2647 pam-fi 1.11 * |||||||||||||||||||||||||||||||||||||||||||||||||
2648 mocchiut 1.1 * -------------------------------------------------
2649     * If the plane has been already included, it just
2650     * computes again the coordinates of the x-y couple
2651     * using improved PFAs
2652     * -------------------------------------------------
2653 pam-fi 1.11 * |||||||||||||||||||||||||||||||||||||||||||||||||
2654 mocchiut 1.1 if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and.
2655     $ YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then
2656    
2657     id=CP_STORE(nplanes-ip+1,ibest)
2658    
2659     is=is_cp(id)
2660     icp=icp_cp(id)
2661     if(ip_cp(id).ne.ip)
2662     $ print*,'OKKIO!!'
2663     $ ,'id ',id,is,icp
2664     $ ,ip_cp(id),ip
2665     icx=clx(ip,icp)
2666     icy=cly(ip,icp)
2667     call xyz_PAM(icx,icy,is,
2668     c $ 'ETA2','ETA2',
2669     $ PFA,PFA,
2670     $ AXV_STORE(nplanes-ip+1,ibest),
2671     $ AYV_STORE(nplanes-ip+1,ibest))
2672     c$$$ call xyz_PAM(icx,icy,is,
2673     c$$$ $ 'COG2','COG2',
2674     c$$$ $ 0.,
2675     c$$$ $ 0.)
2676     xm(nplanes-ip+1) = xPAM
2677     ym(nplanes-ip+1) = yPAM
2678     zm(nplanes-ip+1) = zPAM
2679     xgood(nplanes-ip+1) = 1
2680     ygood(nplanes-ip+1) = 1
2681     resx(nplanes-ip+1) = resxPAM
2682     resy(nplanes-ip+1) = resyPAM
2683    
2684     c dedxtrk(nplanes-ip+1) = (dedx(icx)+dedx(icy))/2. !(1)
2685     dedxtrk_x(nplanes-ip+1)=dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
2686     dedxtrk_y(nplanes-ip+1)=dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
2687    
2688 pam-fi 1.11 * |||||||||||||||||||||||||||||||||||||||||||||||||
2689 mocchiut 1.1 * -------------------------------------------------
2690     * If the plane has NOT been already included,
2691     * it tries to include a COUPLE or a single cluster
2692     * -------------------------------------------------
2693 pam-fi 1.11 * |||||||||||||||||||||||||||||||||||||||||||||||||
2694 mocchiut 1.1 else
2695    
2696     xgood(nplanes-ip+1)=0
2697     ygood(nplanes-ip+1)=0
2698    
2699     * --------------------------------------------------------------
2700     * determine which ladder and sensor are intersected by the track
2701     xP=XV_STORE(nplanes-ip+1,ibest)
2702     yP=YV_STORE(nplanes-ip+1,ibest)
2703     zP=ZV_STORE(nplanes-ip+1,ibest)
2704     call whichsensor(ip,xP,yP,nldt,ist)
2705     * if the track hit the plane in a dead area, go to the next plane
2706     if(nldt.eq.0.or.ist.eq.0)goto 133
2707     * --------------------------------------------------------------
2708    
2709     if(DEBUG)then
2710     print*,
2711     $ '------ Plane ',ip,' intersected on LADDER ',nldt
2712     $ ,' SENSOR ',ist
2713     print*,
2714     $ '------ coord: ',XP,YP
2715     endif
2716    
2717     * ===========================================
2718     * STEP 1 >>>>>>> try to include a new couple
2719     * ===========================================
2720     c if(DEBUG)print*,'>>>> try to include a new couple'
2721     distmin=1000000.
2722     xmm = 0.
2723     ymm = 0.
2724     zmm = 0.
2725     rxmm = 0.
2726     rymm = 0.
2727     dedxmmx = 0. !(1)
2728     dedxmmy = 0. !(1)
2729     idm = 0 !ID of the closer couple
2730     distance=0.
2731     do icp=1,ncp_plane(ip) !loop on couples on plane icp
2732     icx=clx(ip,icp)
2733     icy=cly(ip,icp)
2734     if(LADDER(icx).ne.nldt.or. !If the ladder number does not match
2735 pam-fi 1.3 c $ cl_used(icx).eq.1.or. !or the X cluster is already used
2736     c $ cl_used(icy).eq.1.or. !or the Y cluster is already used
2737     $ cl_used(icx).ne.0.or. !or the X cluster is already used !(3)
2738     $ cl_used(icy).ne.0.or. !or the Y cluster is already used !(3)
2739 mocchiut 1.1 $ .false.)goto 1188 !then jump to next couple.
2740     *
2741     call xyz_PAM(icx,icy,ist,
2742     $ PFA,PFA,
2743     c $ 'ETA2','ETA2',
2744     $ AXV_STORE(nplanes-ip+1,ibest),
2745     $ AYV_STORE(nplanes-ip+1,ibest))
2746    
2747     distance = distance_to(XP,YP)
2748 pam-fi 1.11 distance = distance / RCHI2_STORE(ibest)!<<< MS
2749 mocchiut 1.1 id=id_cp(ip,icp,ist)
2750     if(DEBUG)print*,'( couple ',id
2751     $ ,' ) normalized distance ',distance
2752     if(distance.lt.distmin)then
2753     xmm = xPAM
2754     ymm = yPAM
2755     zmm = zPAM
2756     rxmm = resxPAM
2757     rymm = resyPAM
2758     distmin = distance
2759     idm = id
2760     c dedxmm = (dedx(icx)+dedx(icy))/2. !(1)
2761     dedxmmx = dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
2762     dedxmmy = dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
2763     endif
2764     1188 continue
2765     enddo !end loop on couples on plane icp
2766     if(distmin.le.clinc)then
2767     * -----------------------------------
2768     xm(nplanes-ip+1) = xmm !<<<
2769     ym(nplanes-ip+1) = ymm !<<<
2770     zm(nplanes-ip+1) = zmm !<<<
2771     xgood(nplanes-ip+1) = 1 !<<<
2772     ygood(nplanes-ip+1) = 1 !<<<
2773     resx(nplanes-ip+1)=rxmm !<<<
2774     resy(nplanes-ip+1)=rymm !<<<
2775     c dedxtrk(nplanes-ip+1) = dedxmm !<<< !(1)
2776     dedxtrk_x(nplanes-ip+1) = dedxmmx !(1)
2777     dedxtrk_y(nplanes-ip+1) = dedxmmy !(1)
2778     * -----------------------------------
2779     CP_STORE(nplanes-ip+1,ibest)=idm
2780     if(DEBUG)print*,'%%%% included couple ',idm
2781     $ ,' (norm.dist.= ',distmin,', cut ',clinc,' )'
2782     goto 133 !next plane
2783     endif
2784     * ================================================
2785     * STEP 2 >>>>>>> try to include a single cluster
2786     * either from a couple or single
2787     * ================================================
2788     c if(DEBUG)print*,'>>>> try to include a new cluster'
2789     distmin=1000000.
2790     xmm_A = 0. !---------------------------
2791     ymm_A = 0. ! init variables that
2792     zmm_A = 0. ! define the SINGLET
2793     xmm_B = 0. !
2794     ymm_B = 0. !
2795     zmm_B = 0. !
2796     rxmm = 0. !
2797     rymm = 0. !
2798     dedxmmx = 0. !(1)
2799     dedxmmy = 0. !(1)
2800     iclm=0 !---------------------------
2801     distance=0.
2802    
2803     *----- clusters inside couples -------------------------------------
2804     do icp=1,ncp_plane(ip) !loop on cluster inside couples
2805     icx=clx(ip,icp)
2806     icy=cly(ip,icp)
2807     id=id_cp(ip,icp,ist)
2808     if(LADDER(icx).ne.nldt)goto 11882 !if the ladder number does not match
2809     * !jump to the next couple
2810     *----- try cluster x -----------------------------------------------
2811 pam-fi 1.3 c if(cl_used(icx).eq.1)goto 11881 !if the X cluster is already used
2812     if(cl_used(icx).ne.0)goto 11881 !if the X cluster is already used !(3)
2813 mocchiut 1.1 * !jump to the Y cluster
2814     call xyz_PAM(icx,0,ist,
2815     c $ 'ETA2','ETA2',
2816     $ PFA,PFA,
2817     $ AXV_STORE(nplanes-ip+1,ibest),0.)
2818     distance = distance_to(XP,YP)
2819 pam-fi 1.11 distance = distance / RCHI2_STORE(ibest)!<<< MS
2820 mocchiut 1.1 if(DEBUG)print*,'( cl-X ',icx
2821     $ ,' in cp ',id,' ) normalized distance ',distance
2822     if(distance.lt.distmin)then
2823     xmm_A = xPAM_A
2824     ymm_A = yPAM_A
2825     zmm_A = zPAM_A
2826     xmm_B = xPAM_B
2827     ymm_B = yPAM_B
2828     zmm_B = zPAM_B
2829     rxmm = resxPAM
2830     rymm = resyPAM
2831     distmin = distance
2832     iclm = icx
2833     c dedxmm = dedx(icx) !(1)
2834     dedxmmx = dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
2835     dedxmmy = 0. !(1)
2836     endif
2837     11881 continue
2838     *----- try cluster y -----------------------------------------------
2839 pam-fi 1.3 c if(cl_used(icy).eq.1)goto 11882 !if the Y cluster is already used
2840     if(cl_used(icy).ne.0)goto 11882 !if the Y cluster is already used !(3)
2841 mocchiut 1.1 * !jump to the next couple
2842     call xyz_PAM(0,icy,ist,
2843     c $ 'ETA2','ETA2',
2844     $ PFA,PFA,
2845     $ 0.,AYV_STORE(nplanes-ip+1,ibest))
2846     distance = distance_to(XP,YP)
2847 pam-fi 1.11 distance = distance / RCHI2_STORE(ibest)!<<< MS
2848 mocchiut 1.1 if(DEBUG)print*,'( cl-Y ',icy
2849     $ ,' in cp ',id,' ) normalized distance ',distance
2850     if(distance.lt.distmin)then
2851     xmm_A = xPAM_A
2852     ymm_A = yPAM_A
2853     zmm_A = zPAM_A
2854     xmm_B = xPAM_B
2855     ymm_B = yPAM_B
2856     zmm_B = zPAM_B
2857     rxmm = resxPAM
2858     rymm = resyPAM
2859     distmin = distance
2860     iclm = icy
2861     c dedxmm = dedx(icy) !(1)
2862     dedxmmx = 0. !(1)
2863     dedxmmy = dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
2864     endif
2865     11882 continue
2866     enddo !end loop on cluster inside couples
2867     *----- single clusters -----------------------------------------------
2868     do ic=1,ncls(ip) !loop on single clusters
2869     icl=cls(ip,ic)
2870 pam-fi 1.3 c if(cl_used(icl).eq.1.or. !if the cluster is already used
2871     if(cl_used(icl).ne.0.or. !if the cluster is already used !(3)
2872 mocchiut 1.1 $ LADDER(icl).ne.nldt.or. !or the ladder number does not match
2873     $ .false.)goto 18882 !jump to the next singlet
2874     if(mod(VIEW(icl),2).eq.0)then!<---- X view
2875     call xyz_PAM(icl,0,ist,
2876     c $ 'ETA2','ETA2',
2877     $ PFA,PFA,
2878     $ AXV_STORE(nplanes-ip+1,ibest),0.)
2879     else !<---- Y view
2880     call xyz_PAM(0,icl,ist,
2881     c $ 'ETA2','ETA2',
2882     $ PFA,PFA,
2883     $ 0.,AYV_STORE(nplanes-ip+1,ibest))
2884     endif
2885    
2886     distance = distance_to(XP,YP)
2887 pam-fi 1.11 distance = distance / RCHI2_STORE(ibest)!<<< MS
2888 mocchiut 1.1 if(DEBUG)print*,'( cl-s ',icl
2889     $ ,' ) normalized distance ',distance
2890     if(distance.lt.distmin)then
2891     xmm_A = xPAM_A
2892     ymm_A = yPAM_A
2893     zmm_A = zPAM_A
2894     xmm_B = xPAM_B
2895     ymm_B = yPAM_B
2896     zmm_B = zPAM_B
2897     rxmm = resxPAM
2898     rymm = resyPAM
2899     distmin = distance
2900     iclm = icl
2901     c dedxmm = dedx(icl) !(1)
2902     if(mod(VIEW(icl),2).eq.0)then !<---- X view
2903     dedxmmx = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2)
2904     dedxmmy = 0. !(1)
2905     else !<---- Y view
2906     dedxmmx = 0. !(1)
2907     dedxmmy = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2)
2908     endif
2909     endif
2910     18882 continue
2911     enddo !end loop on single clusters
2912    
2913     if(distmin.le.clinc)then
2914    
2915     CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<
2916     * ----------------------------
2917 pam-fi 1.11 c print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
2918 mocchiut 1.1 if(mod(VIEW(iclm),2).eq.0)then
2919     XGOOD(nplanes-ip+1)=1.
2920     resx(nplanes-ip+1)=rxmm
2921     if(DEBUG)print*,'%%%% included X-cl ',iclm
2922 pam-fi 1.11 c if(.true.)print*,'%%%% included X-cl ',iclm
2923     $ ,'( chi^2, ',RCHI2_STORE(ibest)
2924     $ ,', norm.dist.= ',distmin
2925     $ ,', cut ',clinc,' )'
2926 mocchiut 1.1 else
2927     YGOOD(nplanes-ip+1)=1.
2928     resy(nplanes-ip+1)=rymm
2929     if(DEBUG)print*,'%%%% included Y-cl ',iclm
2930 pam-fi 1.11 c if(.true.)print*,'%%%% included Y-cl ',iclm
2931     $ ,'( chi^2, ',RCHI2_STORE(ibest)
2932     $ ,', norm.dist.= ', distmin
2933     $ ,', cut ',clinc,' )'
2934 mocchiut 1.1 endif
2935 pam-fi 1.11 c print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
2936 mocchiut 1.1 * ----------------------------
2937     xm_A(nplanes-ip+1) = xmm_A
2938     ym_A(nplanes-ip+1) = ymm_A
2939     xm_B(nplanes-ip+1) = xmm_B
2940     ym_B(nplanes-ip+1) = ymm_B
2941     zm(nplanes-ip+1) = (zmm_A+zmm_B)/2.
2942     c dedxtrk(nplanes-ip+1) = dedxmm !<<< !(1)
2943     dedxtrk_x(nplanes-ip+1) = dedxmmx !<<< !(1)
2944     dedxtrk_y(nplanes-ip+1) = dedxmmy !<<< !(1)
2945     * ----------------------------
2946     endif
2947     endif
2948     133 continue
2949     enddo !end loop on planes
2950    
2951    
2952    
2953     return
2954     end
2955    
2956     ***************************************************
2957     * *
2958     * *
2959     * *
2960     * *
2961     * *
2962     * *
2963     **************************************************
2964 pam-fi 1.3 cccccc 12/08/2006 modified by elena ---> (1)
2965     *
2966 mocchiut 1.1 subroutine clean_XYclouds(ibest,iflag)
2967    
2968     include 'commontracker.f'
2969 pam-fi 1.9 include 'level1.f'
2970 mocchiut 1.1 include 'common_momanhough.f'
2971 pam-fi 1.9 c include 'momanhough_init.f'
2972 pam-fi 1.3 include 'level2.f' !(1)
2973 mocchiut 1.1 c include 'calib.f'
2974     c include 'level1.f'
2975    
2976    
2977    
2978     do ip=1,nplanes !loop on planes
2979    
2980     id=CP_STORE(nplanes-ip+1,ibest)
2981     icl=CLS_STORE(nplanes-ip+1,ibest)
2982     if(id.ne.0.or.icl.ne.0)then
2983     if(id.ne.0)then
2984     iclx=clx(ip,icp_cp(id))
2985     icly=cly(ip,icp_cp(id))
2986 pam-fi 1.3 c cl_used(iclx)=1 !tag used clusters
2987     c cl_used(icly)=1 !tag used clusters
2988     cl_used(iclx)=ntrk !tag used clusters !(1)
2989     cl_used(icly)=ntrk !tag used clusters !(1)
2990 mocchiut 1.1 elseif(icl.ne.0)then
2991 pam-fi 1.3 c cl_used(icl)=1 !tag used clusters
2992     cl_used(icl)=ntrk !tag used clusters !1)
2993 mocchiut 1.1 endif
2994    
2995     c if(DEBUG)then
2996     c print*,ip,' <<< ',id
2997     c endif
2998     * -----------------------------
2999     * remove the couple from clouds
3000     * remove also vitual couples containing the
3001     * selected clusters
3002     * -----------------------------
3003     do icp=1,ncp_plane(ip)
3004     if(
3005     $ clx(ip,icp).eq.iclx
3006     $ .or.
3007     $ clx(ip,icp).eq.icl
3008     $ .or.
3009     $ cly(ip,icp).eq.icly
3010     $ .or.
3011     $ cly(ip,icp).eq.icl
3012     $ )then
3013     id=id_cp(ip,icp,1)
3014     if(DEBUG)then
3015     print*,ip,' <<< cp ',id
3016     $ ,' ( cl-x '
3017     $ ,clx(ip,icp)
3018     $ ,' cl-y '
3019     $ ,cly(ip,icp),' ) --> removed'
3020     endif
3021     * -----------------------------
3022     * remove the couple from clouds
3023     do iyz=1,nclouds_yz
3024     if(cpcloud_yz(iyz,abs(id)).ne.0)then
3025     ptcloud_yz(iyz)=ptcloud_yz(iyz)-1
3026     cpcloud_yz(iyz,abs(id))=0
3027     endif
3028     enddo
3029     do ixz=1,nclouds_xz
3030     if(cpcloud_xz(ixz,abs(id)).ne.0)then
3031     ptcloud_xz(ixz)=ptcloud_xz(ixz)-1
3032     cpcloud_xz(ixz,abs(id))=0
3033     endif
3034     enddo
3035     * -----------------------------
3036     endif
3037     enddo
3038    
3039     endif
3040     enddo !end loop on planes
3041    
3042     return
3043     end
3044    
3045    
3046    
3047    
3048    
3049    
3050     * ****************************************************
3051    
3052     subroutine init_level2
3053    
3054     include 'commontracker.f'
3055 pam-fi 1.9 include 'level1.f'
3056 mocchiut 1.1 include 'common_momanhough.f'
3057     include 'level2.f'
3058 pam-fi 1.9 c include 'level1.f'
3059 mocchiut 1.1
3060 pam-fi 1.4 do i=1,nviews
3061     good2(i)=good1(i)
3062     enddo
3063 mocchiut 1.1
3064    
3065     NTRK = 0
3066 pam-fi 1.9 do it=1,NTRKMAX
3067 mocchiut 1.1 IMAGE(IT)=0
3068     CHI2_nt(IT) = -100000.
3069     do ip=1,nplanes
3070     XM_nt(IP,IT) = 0
3071     YM_nt(IP,IT) = 0
3072     ZM_nt(IP,IT) = 0
3073     RESX_nt(IP,IT) = 0
3074     RESY_nt(IP,IT) = 0
3075     XGOOD_nt(IP,IT) = 0
3076     YGOOD_nt(IP,IT) = 0
3077     DEDX_X(IP,IT) = 0
3078     DEDX_Y(IP,IT) = 0
3079 pam-fi 1.3 CLTRX(IP,IT) = 0
3080     CLTRY(IP,IT) = 0
3081 mocchiut 1.1 enddo
3082     do ipa=1,5
3083     AL_nt(IPA,IT) = 0
3084     do ipaa=1,5
3085     coval(ipa,ipaa,IT)=0
3086     enddo
3087     enddo
3088     enddo
3089     nclsx=0
3090     nclsy=0
3091     do ip=1,NSINGMAX
3092     planex(ip)=0
3093     xs(1,ip)=0
3094     xs(2,ip)=0
3095     sgnlxs(ip)=0
3096     planey(ip)=0
3097     ys(1,ip)=0
3098     ys(2,ip)=0
3099     sgnlys(ip)=0
3100     enddo
3101     end
3102    
3103    
3104     ************************************************************
3105     *
3106     *
3107     *
3108     *
3109     *
3110     *
3111     *
3112     ************************************************************
3113    
3114    
3115 pam-fi 1.9 subroutine init_hough
3116    
3117     include 'commontracker.f'
3118     include 'level1.f'
3119     include 'common_momanhough.f'
3120     include 'common_hough.f'
3121     include 'level2.f'
3122    
3123     ntrpt_nt=0
3124     ndblt_nt=0
3125     NCLOUDS_XZ_nt=0
3126     NCLOUDS_YZ_nt=0
3127     do idb=1,ndblt_max_nt
3128     db_cloud_nt(idb)=0
3129     alfayz1_nt(idb)=0
3130     alfayz2_nt(idb)=0
3131     enddo
3132 pam-fi 1.14 do itr=1,ntrpt_max_nt
3133 pam-fi 1.9 tr_cloud_nt(itr)=0
3134     alfaxz1_nt(itr)=0
3135     alfaxz2_nt(itr)=0
3136     alfaxz3_nt(itr)=0
3137     enddo
3138     do idb=1,ncloyz_max
3139     ptcloud_yz_nt(idb)=0
3140     alfayz1_av_nt(idb)=0
3141     alfayz2_av_nt(idb)=0
3142     enddo
3143     do itr=1,ncloxz_max
3144     ptcloud_xz_nt(itr)=0
3145     alfaxz1_av_nt(itr)=0
3146     alfaxz2_av_nt(itr)=0
3147     alfaxz3_av_nt(itr)=0
3148     enddo
3149    
3150     ntrpt=0
3151     ndblt=0
3152     NCLOUDS_XZ=0
3153     NCLOUDS_YZ=0
3154     do idb=1,ndblt_max
3155     db_cloud(idb)=0
3156     cpyz1(idb)=0
3157     cpyz2(idb)=0
3158     alfayz1(idb)=0
3159     alfayz2(idb)=0
3160     enddo
3161 pam-fi 1.14 do itr=1,ntrpt_max
3162 pam-fi 1.9 tr_cloud(itr)=0
3163     cpxz1(itr)=0
3164     cpxz2(itr)=0
3165     cpxz3(itr)=0
3166     alfaxz1(itr)=0
3167     alfaxz2(itr)=0
3168     alfaxz3(itr)=0
3169     enddo
3170     do idb=1,ncloyz_max
3171     ptcloud_yz(idb)=0
3172     alfayz1_av(idb)=0
3173     alfayz2_av(idb)=0
3174     do idbb=1,ncouplemaxtot
3175     cpcloud_yz(idb,idbb)=0
3176     enddo
3177     enddo
3178     do itr=1,ncloxz_max
3179     ptcloud_xz(itr)=0
3180     alfaxz1_av(itr)=0
3181     alfaxz2_av(itr)=0
3182     alfaxz3_av(itr)=0
3183     do itrr=1,ncouplemaxtot
3184     cpcloud_xz(itr,itrr)=0
3185     enddo
3186     enddo
3187     end
3188     ************************************************************
3189     *
3190     *
3191     *
3192     *
3193     *
3194     *
3195     *
3196     ************************************************************
3197    
3198    
3199 mocchiut 1.1 subroutine fill_level2_tracks(ntr)
3200    
3201     * -------------------------------------------------------
3202     * This routine fills the ntr-th element of the variables
3203     * inside the level2_tracks common, which correspond
3204     * to the ntr-th track info.
3205     * -------------------------------------------------------
3206    
3207    
3208     include 'commontracker.f'
3209 pam-fi 1.9 c include 'level1.f'
3210 pam-fi 1.3 include 'level1.f'
3211 pam-fi 1.9 include 'common_momanhough.f'
3212 mocchiut 1.1 include 'level2.f'
3213     include 'common_mini_2.f'
3214 pam-fi 1.9 real sinth,phi,pig
3215 mocchiut 1.1 pig=acos(-1.)
3216    
3217     chi2_nt(ntr) = sngl(chi2)
3218 pam-fi 1.8 nstep_nt(ntr) = nstep
3219 mocchiut 1.1
3220 pam-fi 1.9 phi = al(4)
3221     sinth = al(3)
3222     if(sinth.lt.0)then
3223     sinth = -sinth
3224     phi = phi + pig
3225     endif
3226     npig = aint(phi/(2*pig))
3227     phi = phi - npig*2*pig
3228     if(phi.lt.0)
3229     $ phi = phi + 2*pig
3230     al(4) = phi
3231     al(3) = sinth
3232    
3233 mocchiut 1.1 do i=1,5
3234     al_nt(i,ntr) = sngl(al(i))
3235     do j=1,5
3236     coval(i,j,ntr) = sngl(cov(i,j))
3237     enddo
3238     enddo
3239    
3240     do ip=1,nplanes ! loop on planes
3241     xgood_nt(ip,ntr) = int(xgood(ip))
3242     ygood_nt(ip,ntr) = int(ygood(ip))
3243     xm_nt(ip,ntr) = sngl(xm(ip))
3244     ym_nt(ip,ntr) = sngl(ym(ip))
3245     zm_nt(ip,ntr) = sngl(zm(ip))
3246     RESX_nt(IP,ntr) = sngl(resx(ip))
3247     RESY_nt(IP,ntr) = sngl(resy(ip))
3248     xv_nt(ip,ntr) = sngl(xv(ip))
3249     yv_nt(ip,ntr) = sngl(yv(ip))
3250     zv_nt(ip,ntr) = sngl(zv(ip))
3251     axv_nt(ip,ntr) = sngl(axv(ip))
3252     ayv_nt(ip,ntr) = sngl(ayv(ip))
3253     dedx_x(ip,ntr) = sngl(dedxtrk_x(ip)) !(2)
3254 pam-fi 1.3 dedx_y(ip,ntr) = sngl(dedxtrk_y(ip)) !(2)
3255    
3256     id = CP_STORE(ip,IDCAND)
3257     icl = CLS_STORE(ip,IDCAND)
3258     if(id.ne.0)then
3259     cltrx(ip,ntr) = clx(nplanes-ip+1,icp_cp(id))
3260     cltry(ip,ntr) = cly(nplanes-ip+1,icp_cp(id))
3261     c print*,ip,' ',cltrx(ip,ntr),cltry(ip,ntr)
3262     elseif(icl.ne.0)then
3263     if(mod(VIEW(icl),2).eq.0)cltrx(ip,ntr)=icl
3264     if(mod(VIEW(icl),2).eq.1)cltry(ip,ntr)=icl
3265     c print*,ip,' ',cltrx(ip,ntr),cltry(ip,ntr)
3266     endif
3267    
3268 mocchiut 1.1 enddo
3269    
3270    
3271     end
3272    
3273     subroutine fill_level2_siglets
3274    
3275     * -------------------------------------------------------
3276     * This routine fills the elements of the variables
3277     * inside the level2_singletsx and level2_singletsy commons,
3278     * which store info on clusters outside the tracks
3279     * -------------------------------------------------------
3280    
3281     include 'commontracker.f'
3282 pam-fi 1.9 c include 'level1.f'
3283     include 'calib.f'
3284 mocchiut 1.1 include 'level1.f'
3285 pam-fi 1.9 include 'common_momanhough.f'
3286 mocchiut 1.1 include 'level2.f'
3287     include 'common_xyzPAM.f'
3288    
3289     * count #cluster per plane not associated to any track
3290 pam-fi 1.4 c good2=1!.true.
3291 mocchiut 1.1 nclsx = 0
3292     nclsy = 0
3293    
3294 pam-fi 1.9 do iv = 1,nviews
3295     if( mask_view(iv).ne.0 )good2(iv) = 20+mask_view(iv)
3296     enddo
3297    
3298 mocchiut 1.1 do icl=1,nclstr1
3299     if(cl_used(icl).eq.0)then !cluster not included in any track
3300     ip=nplanes-npl(VIEW(icl))+1
3301     if(mod(VIEW(icl),2).eq.0)then !=== X views
3302     nclsx = nclsx + 1
3303     planex(nclsx) = ip
3304     sgnlxs(nclsx) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)
3305 pam-fi 1.3 clsx(nclsx) = icl
3306 mocchiut 1.1 do is=1,2
3307     c call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)
3308     call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)
3309     xs(is,nclsx) = (xPAM_A+xPAM_B)/2
3310     enddo
3311     c$$$ print*,'nclsx ',nclsx
3312     c$$$ print*,'planex(nclsx) ',planex(nclsx)
3313     c$$$ print*,'sgnlxs(nclsx) ',sgnlxs(nclsx)
3314     c$$$ print*,'xs(1,nclsx) ',xs(1,nclsx)
3315     c$$$ print*,'xs(2,nclsx) ',xs(2,nclsx)
3316     else !=== Y views
3317     nclsy = nclsy + 1
3318     planey(nclsy) = ip
3319     sgnlys(nclsy) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)
3320 pam-fi 1.3 clsy(nclsy) = icl
3321 mocchiut 1.1 do is=1,2
3322     c call xyz_PAM(0,icl,is,' ','COG1',0.,0.)
3323     call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)
3324     ys(is,nclsy) = (yPAM_A+yPAM_B)/2
3325     enddo
3326     c$$$ print*,'nclsy ',nclsy
3327     c$$$ print*,'planey(nclsy) ',planey(nclsy)
3328     c$$$ print*,'sgnlys(nclsy) ',sgnlys(nclsy)
3329     c$$$ print*,'ys(1,nclsy) ',ys(1,nclsy)
3330     c$$$ print*,'ys(2,nclsy) ',ys(2,nclsy)
3331     endif
3332     endif
3333     c print*,icl,cl_used(icl),cl_good(icl),ip,VIEW(icl)!nclsx(ip),nclsy(ip)
3334 pam-fi 1.3
3335     ***** LO METTO QUI PERCHE` NON SO DOVE METTERLO
3336     whichtrack(icl) = cl_used(icl)
3337    
3338 mocchiut 1.1 enddo
3339     end
3340    
3341 pam-fi 1.9 ***************************************************
3342     * *
3343     * *
3344     * *
3345     * *
3346     * *
3347     * *
3348     **************************************************
3349 mocchiut 1.1
3350 pam-fi 1.9 subroutine fill_hough
3351 mocchiut 1.1
3352 pam-fi 1.9 * -------------------------------------------------------
3353     * This routine fills the variables related to the hough
3354     * transform, for the debig n-tuple
3355     * -------------------------------------------------------
3356 mocchiut 1.1
3357 pam-fi 1.9 include 'commontracker.f'
3358     include 'level1.f'
3359     include 'common_momanhough.f'
3360     include 'common_hough.f'
3361     include 'level2.f'
3362 mocchiut 1.1
3363 pam-fi 1.9 if(.false.
3364     $ .or.ntrpt.gt.ntrpt_max_nt
3365     $ .or.ndblt.gt.ndblt_max_nt
3366     $ .or.NCLOUDS_XZ.gt.ncloxz_max
3367     $ .or.NCLOUDS_yZ.gt.ncloyz_max
3368     $ )then
3369     ntrpt_nt=0
3370     ndblt_nt=0
3371     NCLOUDS_XZ_nt=0
3372     NCLOUDS_YZ_nt=0
3373     else
3374     ndblt_nt=ndblt
3375     ntrpt_nt=ntrpt
3376     if(ndblt.ne.0)then
3377     do id=1,ndblt
3378     alfayz1_nt(id)=alfayz1(id) !Y0
3379     alfayz2_nt(id)=alfayz2(id) !tg theta-yz
3380     enddo
3381     endif
3382     if(ndblt.ne.0)then
3383     do it=1,ntrpt
3384     alfaxz1_nt(it)=alfaxz1(it) !X0
3385     alfaxz2_nt(it)=alfaxz2(it) !tg theta-xz
3386     alfaxz3_nt(it)=alfaxz3(it) !1/r
3387     enddo
3388     endif
3389     nclouds_yz_nt=nclouds_yz
3390     nclouds_xz_nt=nclouds_xz
3391     if(nclouds_yz.ne.0)then
3392     nnn=0
3393     do iyz=1,nclouds_yz
3394     ptcloud_yz_nt(iyz)=ptcloud_yz(iyz)
3395     alfayz1_av_nt(iyz)=alfayz1_av(iyz)
3396     alfayz2_av_nt(iyz)=alfayz2_av(iyz)
3397     nnn=nnn+ptcloud_yz(iyz)
3398     enddo
3399     do ipt=1,nnn
3400     db_cloud_nt(ipt)=db_cloud(ipt)
3401     enddo
3402     endif
3403     if(nclouds_xz.ne.0)then
3404     nnn=0
3405     do ixz=1,nclouds_xz
3406     ptcloud_xz_nt(ixz)=ptcloud_xz(ixz)
3407     alfaxz1_av_nt(ixz)=alfaxz1_av(ixz)
3408     alfaxz2_av_nt(ixz)=alfaxz2_av(ixz)
3409     alfaxz3_av_nt(ixz)=alfaxz3_av(ixz)
3410     nnn=nnn+ptcloud_xz(ixz)
3411     enddo
3412     do ipt=1,nnn
3413     tr_cloud_nt(ipt)=tr_cloud(ipt)
3414     enddo
3415     endif
3416     endif
3417     end
3418    

  ViewVC Help
Powered by ViewVC 1.1.23