/[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.43 - (hide annotations) (download)
Mon Dec 14 10:51:40 2009 UTC (15 years ago) by pam-fi
Branch: MAIN
Changes since 1.42: +18 -8 lines 
bug fixed: check on maximum number of tracks
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.20 c print*,'======================================================'
24     c$$$ do ic=1,NCLSTR1
25     c$$$ if(.false.
26     c$$$ $ .or.nsatstrips(ic).gt.0
27     c$$$c $ .or.nbadstrips(0,ic).gt.0
28     c$$$c $ .or.nbadstrips(4,ic).gt.0
29     c$$$c $ .or.nbadstrips(3,ic).gt.0
30     c$$$ $ .or..false.)then
31     c$$$ print*,'--- cl-',ic,' ------------------------'
32     c$$$ istart = INDSTART(IC)
33     c$$$ istop = TOTCLLENGTH
34     c$$$ if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
35     c$$$ print*,'ADC ',(CLADC(i),i=istart,istop)
36     c$$$ print*,'s/n ',(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop)
37     c$$$ print*,'sgnl ',(CLSIGNAL(i),i=istart,istop)
38     c$$$ print*,'strip ',(i-INDMAX(ic),i=istart,istop)
39     c$$$ print*,'view ',VIEW(ic)
40     c$$$ print*,'maxs ',MAXS(ic)
41     c$$$ print*,'COG4 ',cog(4,ic)
42     c$$$ ff = fbad_cog(4,ic)
43     c$$$ print*,'fbad ',ff
44     c$$$ print*,(CLBAD(i),i=istart,istop)
45     c$$$ bb=nbadstrips(0,ic)
46     c$$$ print*,'#BAD (tot)',bb
47     c$$$ bb=nbadstrips(4,ic)
48     c$$$ print*,'#BAD (4)',bb
49     c$$$ bb=nbadstrips(3,ic)
50     c$$$ print*,'#BAD (3)',bb
51     c$$$ ss=nsatstrips(ic)
52     c$$$ print*,'#saturated ',ss
53     c$$$ endif
54     c$$$ enddo
55 mocchiut 1.1
56     *-------------------------------------------------------------------------------
57     * STEP 1
58     *-------------------------------------------------------------------------------
59     * X-Y cluster association
60     *
61     * Clusters are associated to form COUPLES
62     * Clusters not associated in any couple are called SINGLETS
63     *
64     * Track identification (Hough transform) and fitting is first done on couples.
65     * Hence singlets are possibly added to the track.
66     *
67     * Variables assigned by the routine "cl_to_couples" are those in the
68     * common blocks:
69     * - common/clusters/cl_good
70     * - common/couples/clx,cly,ncp_plane,ncp_tot,cp_useds1,cp_useds2
71     * - common/singlets/ncls,cls,cl_single
72     *-------------------------------------------------------------------------------
73     *-------------------------------------------------------------------------------
74    
75     call cl_to_couples(iflag)
76     if(iflag.eq.1)then !bad event
77 pam-fi 1.9 goto 880 !go to next event
78 mocchiut 1.1 endif
79 pam-fi 1.37 if(ncp_tot.eq.0)goto 880 !go to next event
80 mocchiut 1.1 *-----------------------------------------------------
81     *-----------------------------------------------------
82     * HOUGH TRASFORM
83     *-----------------------------------------------------
84     *-----------------------------------------------------
85    
86    
87     *-------------------------------------------------------------------------------
88     * STEP 2
89     *-------------------------------------------------------------------------------
90     *
91     * Association of couples to form
92     * - DOUBLETS in YZ view
93     * - TRIPLETS in XZ view
94     *
95     * Variables assigned by the routine "cp_to_doubtrip" are those in the
96     * common blocks:
97     * - common/hough_param/
98     * $ alfayz1, !Y0
99     * $ alfayz2, !tg theta-yz
100     * $ alfaxz1, !X0
101     * $ alfaxz2, !tg theta-xz
102     * $ alfaxz3 !1/r
103     * - common/doublets/ndblt,cpyz1,cpyz2
104     * - common/triplets/ntrpt,cpxz1,cpxz2,cpxz3
105     *-------------------------------------------------------------------------------
106     *-------------------------------------------------------------------------------
107    
108 pam-fi 1.20
109 mocchiut 1.1 call cp_to_doubtrip(iflag)
110     if(iflag.eq.1)then !bad event
111 pam-fi 1.9 goto 880 !go to next event
112 mocchiut 1.1 endif
113 pam-fi 1.37 if(ndblt.eq.0.or.ntrpt.eq.0)goto 880 !go to next event
114 mocchiut 1.1
115    
116     *-------------------------------------------------------------------------------
117     * STEP 3
118     *-------------------------------------------------------------------------------
119     *
120     * Classification of doublets and triplets to form CLOUDS,
121     * according to distance in parameter space.
122     *
123     * cloud = cluster of points (doublets/triplets) in parameter space
124     *
125     *
126     *
127     * Variables assigned by the routine "doub_to_YZcloud" are those in the
128     * common blocks:
129     * - common/clouds_yz/
130     * $ nclouds_yz
131     * $ ,alfayz1_av,alfayz2_av
132     * $ ,ptcloud_yz,db_cloud,cpcloud_yz
133     *
134     * Variables assigned by the routine "trip_to_XZcloud" are those in the
135     * common blocks:
136     * common/clouds_xz/
137     * $ nclouds_xz xz2_av,alfaxz3_av
138     * $ ,ptcloud_xz,tr_cloud,cpcloud_xz
139     *-------------------------------------------------------------------------------
140     *-------------------------------------------------------------------------------
141 pam-fi 1.9 * count number of hit planes
142     planehit=0
143     do np=1,nplanes
144     if(ncp_plane(np).ne.0)then
145     planehit=planehit+1
146     endif
147     enddo
148     if(planehit.lt.3) goto 880 ! exit
149    
150     nptxz_min=x_min_start
151     nplxz_min=x_min_start
152    
153     nptyz_min=y_min_start
154     nplyz_min=y_min_start
155    
156     cutdistyz=cutystart
157     cutdistxz=cutxstart
158 mocchiut 1.1
159 pam-fi 1.9 878 continue
160 mocchiut 1.1 call doub_to_YZcloud(iflag)
161     if(iflag.eq.1)then !bad event
162     goto 880 !fill ntp and go to next event
163 pam-fi 1.9 endif
164 pam-fi 1.37 * ------------------------------------------------
165     * first try to release the tolerance
166     * ------------------------------------------------
167 pam-fi 1.9 if(nclouds_yz.eq.0.and.cutdistyz.lt.maxcuty)then
168 pam-fi 1.37 if(cutdistyz.lt.maxcuty/2)then
169     cutdistyz=cutdistyz+cutystep
170     else
171     cutdistyz=cutdistyz+(3*cutystep)
172     endif
173     if(DEBUG.EQ.1)print*,'try again - cutdistyz ',cutdistyz
174     goto 878
175     endif
176     * ------------------------------------------------
177     * hence reduce the minimum number of plane
178     * ------------------------------------------------
179     if(nclouds_yz.eq.0.and.nplyz_min.gt.3)then
180     nplyz_min=nplyz_min-1
181     if(DEBUG.EQ.1)print*,'try again - nplyz_min ',nplyz_min
182     goto 878
183     endif
184    
185     if(nclouds_yz.eq.0)goto 880 !go to next event
186    
187 pam-fi 1.9
188 pam-fi 1.37 ccc if(planehit.eq.3) goto 881
189 pam-fi 1.9
190     879 continue
191 mocchiut 1.1 call trip_to_XZcloud(iflag)
192     if(iflag.eq.1)then !bad event
193     goto 880 !fill ntp and go to next event
194     endif
195 pam-fi 1.37 * ------------------------------------------------
196     * first try to release the tolerance
197     * ------------------------------------------------
198 pam-fi 1.9 if(nclouds_xz.eq.0.and.cutdistxz.lt.maxcutx)then
199     cutdistxz=cutdistxz+cutxstep
200 pam-fi 1.37 if(DEBUG.EQ.1)print*,'try again - cutdistxz ',cutdistxz
201 pam-fi 1.9 goto 879
202     endif
203 pam-fi 1.37 * ------------------------------------------------
204     * hence reduce the minimum number of plane
205     * ------------------------------------------------
206     if(nclouds_xz.eq.0.and.nplxz_min.gt.3)then
207     nplxz_min=nplxz_min-1
208     if(DEBUG.EQ.1)print*,'try again - nplxz_min ',nplxz_min
209     goto 879
210     endif
211    
212     if(nclouds_xz.eq.0)goto 880 !go to next event
213 pam-fi 1.9
214    
215 pam-fi 1.37 c$$$ 881 continue
216     c$$$* if there is at least three planes on the Y view decreases cuts on X view
217     c$$$ if(nclouds_xz.eq.0.and.nclouds_yz.gt.0.and.
218     c$$$ $ nplxz_min.ne.y_min_start)then
219     c$$$ nptxz_min=x_min_step
220     c$$$ nplxz_min=x_min_start-x_min_step
221     c$$$ goto 879
222     c$$$ endif
223 pam-fi 1.9
224 mocchiut 1.1 880 return
225     end
226    
227     ************************************************************
228    
229    
230     subroutine track_fitting(iflag)
231    
232     include 'commontracker.f'
233 pam-fi 1.9 include 'level1.f'
234 mocchiut 1.1 include 'common_momanhough.f'
235     include 'common_mech.f'
236     include 'common_xyzPAM.f'
237     include 'common_mini_2.f'
238     include 'calib.f'
239     include 'level2.f'
240    
241 pam-fi 1.9 c include 'momanhough_init.f'
242 mocchiut 1.1
243     logical FIMAGE !
244 pam-fi 1.28 real trackimage(NTRACKSMAX)
245 pam-fi 1.11 real*8 AL_GUESS(5)
246 mocchiut 1.1
247     *-------------------------------------------------------------------------------
248     * STEP 4 (ITERATED until any other physical track isn't found)
249     *-------------------------------------------------------------------------------
250     *
251     * YZ and XZ clouds are combined in order to obtain the initial guess
252     * of the candidate-track parameters.
253     * A minimum number of matching couples between YZ and XZ clouds is required.
254     *
255     * A TRACK CANDIDATE is defined by
256     * - the couples resulting from the INTERSECTION of the two clouds, and
257     * - the associated track parameters (evaluated by performing a zero-order
258     * track fitting)
259     *
260     * The NTRACKS candidate-track parameters are stored in common block:
261     *
262     * - common/track_candidates/NTRACKS,AL_STORE
263     * $ ,XV_STORE,YV_STORE,ZV_STORE
264     * $ ,XM_STORE,YM_STORE,ZM_STORE
265     * $ ,RESX_STORE,RESY_STORE
266     * $ ,AXV_STORE,AYV_STORE
267     * $ ,XGOOD_STORE,YGOOD_STORE
268     * $ ,CP_STORE,RCHI2_STORE
269     *
270     *-------------------------------------------------------------------------------
271     *-------------------------------------------------------------------------------
272 pam-fi 1.37 ccc ntrk=0 !counter of identified physical tracks
273 mocchiut 1.1
274 mocchiut 1.40 c11111 continue !<<<<<<< come here when performing a new search
275     continue !<<<<<<< come here when performing a new search
276 mocchiut 1.1
277 pam-fi 1.37 if(nclouds_xz.eq.0)goto 880 !go to next event
278     if(nclouds_yz.eq.0)goto 880 !go to next event
279    
280 mocchiut 1.1 c iflag=0
281     call clouds_to_ctrack(iflag)
282     if(iflag.eq.1)then !no candidate tracks found
283     goto 880 !fill ntp and go to next event
284     endif
285 pam-fi 1.37 if(ntracks.eq.0)goto 880 !go to next event
286 mocchiut 1.1
287     FIMAGE=.false. !processing best track (not track image)
288     ibest=0 !best track among candidates
289     iimage=0 !track image
290     * ------------- select the best track -------------
291 pam-fi 1.9 c$$$ rchi2best=1000000000.
292     c$$$ do i=1,ntracks
293     c$$$ if(RCHI2_STORE(i).lt.rchi2best.and.
294     c$$$ $ RCHI2_STORE(i).gt.0)then
295     c$$$ ibest=i
296     c$$$ rchi2best=RCHI2_STORE(i)
297     c$$$ endif
298     c$$$ enddo
299     c$$$ if(ibest.eq.0)goto 880 !>> no good candidates
300    
301 pam-fi 1.15 * -------------------------------------------------------
302     * order track-candidates according to:
303     * 1st) decreasing n.points
304     * 2nd) increasing chi**2
305     * -------------------------------------------------------
306 mocchiut 1.41 rchi2best=1000000000.
307 pam-fi 1.20 ndofbest=0
308 mocchiut 1.1 do i=1,ntracks
309 pam-fi 1.20 ndof=0
310     do ii=1,nplanes
311     ndof=ndof
312     $ +int(xgood_store(ii,i))
313     $ +int(ygood_store(ii,i))
314     enddo
315     if(ndof.gt.ndofbest)then
316 pam-fi 1.15 ibest=i
317     rchi2best=RCHI2_STORE(i)
318 pam-fi 1.20 ndofbest=ndof
319     elseif(ndof.eq.ndofbest)then
320 pam-fi 1.15 if(RCHI2_STORE(i).lt.rchi2best.and.
321     $ RCHI2_STORE(i).gt.0)then
322 mocchiut 1.1 ibest=i
323     rchi2best=RCHI2_STORE(i)
324 pam-fi 1.20 ndofbest=ndof
325     endif
326 pam-fi 1.9 endif
327 mocchiut 1.1 enddo
328 pam-fi 1.15
329    
330 mocchiut 1.1 if(ibest.eq.0)goto 880 !>> no good candidates
331     *-------------------------------------------------------------------------------
332     * The best track candidate (ibest) is selected and a new fitting is performed.
333     * Previous to this, the track is refined by:
334     * - possibly adding new COUPLES or SINGLETS from the missing planes
335     * - evaluating the coordinates with improved PFAs
336     * ( angle-dependent ETA algorithms )
337     *-------------------------------------------------------------------------------
338    
339     1212 continue !<<<<< come here to fit track-image
340    
341     if(.not.FIMAGE)then !processing best candidate
342     icand=ibest
343     else !processing image
344     icand=iimage
345     iimage=0
346     endif
347     if(icand.eq.0)then
348 pam-fi 1.28 if(VERBOSE.EQ.1)then
349 mocchiut 1.21 print*,'HAI FATTO UN CASINO!!!!!! icand = ',icand
350     $ ,ibest,iimage
351     endif
352 mocchiut 1.1 return
353     endif
354    
355     * *-*-*-*-*-*-*-*-*-*-*-*-*-*-*
356     call refine_track(icand)
357     * *-*-*-*-*-*-*-*-*-*-*-*-*-*-*
358    
359     * **********************************************************
360     * ************************** FIT *** FIT *** FIT *** FIT ***
361     * **********************************************************
362 pam-fi 1.11 call guess()
363     do i=1,5
364     AL_GUESS(i)=AL(i)
365     enddo
366    
367 mocchiut 1.1 do i=1,5
368 pam-fi 1.3 AL(i)=dble(AL_STORE(i,icand))
369 mocchiut 1.1 enddo
370 pam-fi 1.11
371 pam-fi 1.3 IDCAND = icand !fitted track-candidate
372 mocchiut 1.1 ifail=0 !error flag in chi2 computation
373     jstep=0 !# minimization steps
374 pam-fi 1.9
375 pam-fi 1.8 iprint=0
376 pam-fi 1.28 c if(DEBUG.EQ.1)iprint=1
377     if(VERBOSE.EQ.1)iprint=1
378     if(DEBUG.EQ.1)iprint=2
379 pam-fi 1.8 call mini2(jstep,ifail,iprint)
380 pam-fi 1.38 cc if(ifail.ne.0) then
381     if(ifail.ne.0.or.CHI2.ne.CHI2) then !new
382     if(CHI2.ne.CHI2)CHI2=-9999. !new
383 pam-fi 1.28 if(VERBOSE.EQ.1)then
384 mocchiut 1.1 print *,
385 pam-fi 1.11 $ '*** MINIMIZATION FAILURE *** (after refinement) '
386 mocchiut 1.1 $ ,iev
387     endif
388     endif
389    
390 pam-fi 1.28 if(DEBUG.EQ.1)then
391 mocchiut 1.1 print*,'----------------------------- improved track coord'
392     22222 format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)
393     do ip=1,6
394     write(*,22222)ip,zm(ip),xm(ip),ym(ip)
395     $ ,xm_A(ip),ym_A(ip),xm_B(ip),ym_B(ip)
396     $ ,xgood(ip),ygood(ip),resx(ip),resy(ip)
397     enddo
398     endif
399    
400     c rchi2=chi2/dble(ndof)
401 pam-fi 1.28 if(DEBUG.EQ.1)then
402 mocchiut 1.1 print*,' '
403     print*,'****** SELECTED TRACK *************'
404     print*,'# R. chi2 RIG'
405     print*,' --- ',chi2,' --- '
406     $ ,1./abs(AL(5))
407     print*,'***********************************'
408     endif
409     * **********************************************************
410     * ************************** FIT *** FIT *** FIT *** FIT ***
411     * **********************************************************
412    
413    
414     * ------------- search if the track has an IMAGE -------------
415     * ------------- (also this is stored ) -------------
416     if(FIMAGE)goto 122 !>>> jump! (this is already an image)
417 pam-fi 1.28
418     * -----------------------------------------------------
419     * first check if the track is ambiguous
420     * -----------------------------------------------------
421     * (modified on august 2007 by ElenaV)
422     is1=0
423     do ip=1,NPLANES
424     if(SENSOR_STORE(ip,icand).ne.0)then
425     is1=SENSOR_STORE(ip,icand)
426     if(ip.eq.6)is1=3-is1 !last plane inverted
427     endif
428     enddo
429     if(is1.eq.0)then
430     if(WARNING.EQ.1)print*,'** WARNING ** sensor=0'
431     goto 122 !jump
432     endif
433     do ip=1,NPLANES
434     is2 = SENSOR_STORE(ip,icand) !sensor
435     if(ip.eq.6.and.is2.ne.0)is2=3-is2 !last plane inverted
436     if(
437     $ (is1.ne.is2.and.is2.ne.0)
438     $ )goto 122 !jump (this track cannot have an image)
439     enddo
440     if(DEBUG.eq.1)print*,' >>> ambiguous track! '
441     * ---------------------------------------------------------------
442     * take the candidate with the greatest number of matching couples
443     * if more than one satisfies the requirement,
444     * choose the one with more points and lower chi2
445     * ---------------------------------------------------------------
446     * count the number of matching couples
447     ncpmax = 0
448     iimage = 0 !id of candidate with better matching
449 mocchiut 1.1 do i=1,ntracks
450 pam-fi 1.28 ncp=0
451 mocchiut 1.1 do ip=1,nplanes
452 pam-fi 1.28 if(CP_STORE(nplanes-ip+1,icand).ne.0)then
453     if(
454     $ CP_STORE(nplanes-ip+1,i).ne.0
455     $ .and.
456     $ CP_STORE(nplanes-ip+1,icand).eq.
457     $ -1*CP_STORE(nplanes-ip+1,i)
458     $ )then
459     ncp=ncp+1 !ok
460     elseif(
461     $ CP_STORE(nplanes-ip+1,i).ne.0
462     $ .and.
463     $ CP_STORE(nplanes-ip+1,icand).ne.
464     $ -1*CP_STORE(nplanes-ip+1,i)
465     $ )then
466     ncp = 0
467     goto 117 !it is not an image candidate
468     else
469    
470     endif
471     endif
472 mocchiut 1.1 enddo
473 pam-fi 1.28 117 continue
474     trackimage(i)=ncp !number of matching couples
475     if(ncp>ncpmax)then
476     ncpmax=ncp
477     iimage=i
478 mocchiut 1.1 endif
479     enddo
480 pam-fi 1.28 * check if there are more than one image candidates
481     ngood=0
482     do i=1,ntracks
483     if( ncpmax.ne.0.and.trackimage(i).eq.ncpmax )ngood=ngood+1
484     enddo
485 pam-fi 1.37 if(DEBUG.eq.1)print*,' n.image-candidates : ',ngood
486 pam-fi 1.28 * if there are, choose the best one
487 pam-fi 1.37 if(ngood.gt.0)then
488 pam-fi 1.28 * -------------------------------------------------------
489     * order track-candidates according to:
490     * 1st) decreasing n.points
491     * 2nd) increasing chi**2
492     * -------------------------------------------------------
493     rchi2best=1000000000.
494     ndofbest=0
495     do i=1,ntracks
496     if( trackimage(i).eq.ncpmax )then
497     ndof=0
498     do ii=1,nplanes
499     ndof=ndof
500     $ +int(xgood_store(ii,i))
501     $ +int(ygood_store(ii,i))
502     enddo
503     if(ndof.gt.ndofbest)then
504     iimage=i
505     rchi2best=RCHI2_STORE(i)
506     ndofbest=ndof
507     elseif(ndof.eq.ndofbest)then
508     if(RCHI2_STORE(i).lt.rchi2best.and.
509     $ RCHI2_STORE(i).gt.0)then
510     iimage=i
511     rchi2best=RCHI2_STORE(i)
512     ndofbest=ndof
513     endif
514     endif
515     endif
516     enddo
517 pam-fi 1.37
518     if(DEBUG.EQ.1)then
519     print*,'Track candidate ',iimage
520     $ ,' >>> TRACK IMAGE >>> of'
521     $ ,ibest
522     endif
523 pam-fi 1.28
524     endif
525    
526    
527 mocchiut 1.1 122 continue
528    
529 pam-fi 1.28
530 mocchiut 1.1 * --- and store the results --------------------------------
531     if(ntrk.eq.NTRKMAX)then
532 pam-fi 1.28 if(verbose.eq.1)
533 mocchiut 1.1 $ print*,
534     $ '** warning ** number of identified '//
535     $ 'tracks exceeds vector dimension '
536     $ ,'( ',NTRKMAX,' )'
537     cc good2=.false.
538     goto 880 !fill ntp and go to next event
539     endif
540 pam-fi 1.43
541     ntrk = ntrk + 1 !counter of found tracks
542     if(.not.FIMAGE
543     $ .and.iimage.eq.0) image(ntrk)= 0
544     if(.not.FIMAGE
545     $ .and.iimage.ne.0)image(ntrk)=ntrk+1 !this is the image of the next
546     if(FIMAGE) image(ntrk)=ntrk-1 !this is the image of the previous
547     call fill_level2_tracks(ntrk) !==> good2=.true.
548    
549     c$$$ if(ntrk.eq.NTRKMAX)then
550     c$$$ if(verbose.eq.1)
551     c$$$ $ print*,
552     c$$$ $ '** warning ** number of identified '//
553     c$$$ $ 'tracks exceeds vector dimension '
554     c$$$ $ ,'( ',NTRKMAX,' )'
555     c$$$cc good2=.false.
556     c$$$ goto 880 !fill ntp and go to next event
557     c$$$ endif
558 mocchiut 1.1 if(iimage.ne.0)then
559     FIMAGE=.true. !
560     goto 1212 !>>> fit image-track
561     endif
562    
563    
564     880 return
565     end
566    
567    
568    
569     ************************************************************
570     ************************************************************
571     ************************************************************
572     ************************************************************
573     *
574     * This routine provides the coordinates (in cm) in the PAMELA reference system:
575     * - of the point associated with a COUPLE ---> (xPAM,yPAM,zPAM)
576     * - of the extremes of the segment
577     * associated with a SINGLET ---------------> (xPAM_A,yPAM_A,zPAM_A)
578     * ---> (xPAM_B,yPAM_B,zPAM_B)
579     *
580     * It also assigns the spatial resolution to the evaluated coordinates,
581     * as a function (in principle) of the multiplicity, the angle, the PFA etc...
582     *
583     *
584     * To call the routine you must pass the arguments:
585     * icx - ID of cluster x
586     * icy - ID of cluster y
587     * sensor - sensor (1,2)
588     * PFAx - Position Finding Algorithm in x (COG2,ETA2,...)
589     * PFAy - Position Finding Algorithm in y (COG2,ETA2,...)
590     * angx - Projected angle in x
591     * angy - Projected angle in y
592 pam-fi 1.18 * bfx - x-component of magnetci field
593     * bfy - y-component of magnetic field
594 mocchiut 1.1 *
595     * --------- COUPLES -------------------------------------------------------
596     * The couple defines a point in the space.
597     * The coordinates of the point are evaluated as follows:
598     * 1 - the corrected coordinates relative to the sensor are evaluated
599     * according to the chosen PFA --> (xi,yi,0)
600     * 2 - coordinates are rotated and traslated, according to the aligmnet
601     * parameters, and expressed in the reference system of the mechanical
602     * sensor --> (xrt,yrt,zrt)
603     * 3 - coordinates are finally converted to the PAMELA reference system
604     * --> (xPAM,yPAM,zPAM)
605     *
606     * --------- SINGLETS -------------------------------------------------------
607     * Since a coordinate is missing, the singlet defines not a point
608     * in the space but a segment AB (parallel to the strips).
609     * In this case the routine returns the coordinates in the PAMELA reference
610     * system of the two extremes A and B of the segment:
611     * --> (xPAM_A,yPAM_A,zPAM_A)
612     * --> (xPAM_B,yPAM_B,zPAM_B)
613     *
614     * ==========================================================
615     *
616     * The output of the routine is stored in the commons:
617     *
618     * double precision xPAM,yPAM,zPAM
619     * common/coord_xyz_PAM/xPAM,yPAM,zPAM
620     *
621     * double precision xPAM_A,yPAM_A,zPAM_A
622     * double precision xPAM_B,yPAM_B,zPAM_B
623     * common/coord_AB_PAM/xPAM_A,yPAM_A,zPAM_A,xPAM_B,yPAM_B,zPAM_B
624     *
625     * double precision resxPAM,resyPAM
626     * common/resolution_PAM/resxPAM,resyPAM
627     *
628     * (in file common_xyzPAM.f)
629     *
630     *
631    
632 pam-fi 1.18 subroutine xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
633 mocchiut 1.1
634    
635     include 'commontracker.f'
636 pam-fi 1.9 include 'level1.f'
637 mocchiut 1.1 include 'calib.f'
638     include 'common_align.f'
639     include 'common_mech.f'
640     include 'common_xyzPAM.f'
641    
642     integer icx,icy !X-Y cluster ID
643     integer sensor
644     integer viewx,viewy
645     character*4 PFAx,PFAy !PFA to be used
646 pam-fi 1.18 real ax,ay !X-Y geometric angle
647     real angx,angy !X-Y effective angle
648     real bfx,bfy !X-Y b-field components
649 mocchiut 1.1
650     real stripx,stripy
651    
652 bongi 1.32 double precision xi,yi,zi
653     double precision xi_A,yi_A,zi_A
654     double precision xi_B,yi_B,zi_B
655 mocchiut 1.1 double precision xrt,yrt,zrt
656     double precision xrt_A,yrt_A,zrt_A
657     double precision xrt_B,yrt_B,zrt_B
658    
659    
660     parameter (ndivx=30)
661 pam-fi 1.23
662    
663 mocchiut 1.1
664     resxPAM = 0
665     resyPAM = 0
666    
667 bongi 1.32 xPAM = 0.D0
668     yPAM = 0.D0
669     zPAM = 0.D0
670     xPAM_A = 0.D0
671     yPAM_A = 0.D0
672     zPAM_A = 0.D0
673     xPAM_B = 0.D0
674     yPAM_B = 0.D0
675     zPAM_B = 0.D0
676 pam-fi 1.20
677 pam-fi 1.26 if(sensor.lt.1.or.sensor.gt.2)then
678     print*,'xyz_PAM ***ERROR*** wrong input '
679     print*,'sensor ',sensor
680     icx=0
681     icy=0
682     endif
683    
684 mocchiut 1.1 * -----------------
685     * CLUSTER X
686 pam-fi 1.26 * -----------------
687 mocchiut 1.1 if(icx.ne.0)then
688 pam-fi 1.18
689 pam-fi 1.20 viewx = VIEW(icx)
690     nldx = nld(MAXS(icx),VIEW(icx))
691     nplx = npl(VIEW(icx))
692 pam-fi 1.31 c resxPAM = RESXAV
693 pam-fi 1.20 stripx = float(MAXS(icx))
694 pam-fi 1.26
695     if(
696     $ viewx.lt.1.or.
697     $ viewx.gt.12.or.
698     $ nldx.lt.1.or.
699     $ nldx.gt.3.or.
700     $ stripx.lt.1.or.
701     $ stripx.gt.3072.or.
702     $ .false.)then
703     print*,'xyz_PAM ***ERROR*** wrong input '
704     print*,'icx ',icx,'view ',viewx,'nld ',nldx,'strip ',stripx
705     icx = 0
706     goto 10
707     endif
708    
709 pam-fi 1.18 * --------------------------
710     * magnetic-field corrections
711     * --------------------------
712 pam-fi 1.31 stripx = stripx + fieldcorr(viewx,bfy)
713     angx = effectiveangle(ax,viewx,bfy)
714 pam-fi 1.20
715 pam-fi 1.31 call applypfa(PFAx,icx,angx,corr,res)
716     stripx = stripx + corr
717     resxPAM = res
718 pam-fi 1.18
719 mocchiut 1.41 10 continue
720     endif
721 pam-fi 1.26
722 mocchiut 1.1 * -----------------
723     * CLUSTER Y
724     * -----------------
725    
726     if(icy.ne.0)then
727 pam-fi 1.18
728 mocchiut 1.1 viewy = VIEW(icy)
729     nldy = nld(MAXS(icy),VIEW(icy))
730     nply = npl(VIEW(icy))
731 pam-fi 1.31 c resyPAM = RESYAV
732 pam-fi 1.18 stripy = float(MAXS(icy))
733 mocchiut 1.1
734 pam-fi 1.26 if(
735     $ viewy.lt.1.or.
736     $ viewy.gt.12.or.
737     $ nldy.lt.1.or.
738     $ nldy.gt.3.or.
739     $ stripy.lt.1.or.
740     $ stripy.gt.3072.or.
741     $ .false.)then
742     print*,'xyz_PAM ***ERROR*** wrong input '
743     print*,'icy ',icy,'view ',viewy,'nld ',nldy,'strip ',stripy
744     icy = 0
745     goto 20
746     endif
747    
748 mocchiut 1.1 if(icx.ne.0.and.(nply.ne.nplx.or.nldy.ne.nldx))then
749 pam-fi 1.28 if(DEBUG.EQ.1) then
750 mocchiut 1.21 print*,'xyz_PAM ***ERROR*** invalid cluster couple!!! '
751     $ ,icx,icy
752     endif
753 mocchiut 1.1 goto 100
754     endif
755 pam-fi 1.31
756 pam-fi 1.18 * --------------------------
757     * magnetic-field corrections
758     * --------------------------
759 pam-fi 1.31 stripy = stripy + fieldcorr(viewy,bfx)
760     angy = effectiveangle(ay,viewy,bfx)
761 mocchiut 1.1
762 pam-fi 1.31 call applypfa(PFAy,icy,angy,corr,res)
763     stripy = stripy + corr
764     resyPAM = res
765 pam-fi 1.20
766 mocchiut 1.41 20 continue
767     endif
768 mocchiut 1.1
769 pam-fi 1.17
770 mocchiut 1.1 c===========================================================
771     C COUPLE
772     C===========================================================
773     if(icx.ne.0.and.icy.ne.0)then
774    
775     c------------------------------------------------------------------------
776     c (xi,yi,zi) = mechanical coordinates in the silicon sensor frame
777     c------------------------------------------------------------------------
778 pam-fi 1.4 if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
779     $ .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
780 pam-fi 1.28 if(DEBUG.EQ.1) then
781 mocchiut 1.21 print*,'xyz_PAM (couple):',
782     $ ' WARNING: false X strip: strip ',stripx
783     endif
784 pam-fi 1.4 endif
785 bongi 1.32 xi = dcoordsi(stripx,viewx)
786     yi = dcoordsi(stripy,viewy)
787     zi = 0.D0
788 mocchiut 1.1
789     c------------------------------------------------------------------------
790     c (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame
791     c------------------------------------------------------------------------
792     c N.B. I convert angles from microradiants to radiants
793    
794     xrt = xi
795     $ - omega(nplx,nldx,sensor)*yi
796     $ + gamma(nplx,nldx,sensor)*zi
797     $ + dx(nplx,nldx,sensor)
798    
799     yrt = omega(nplx,nldx,sensor)*xi
800     $ + yi
801     $ - beta(nplx,nldx,sensor)*zi
802     $ + dy(nplx,nldx,sensor)
803    
804     zrt = -gamma(nplx,nldx,sensor)*xi
805     $ + beta(nplx,nldx,sensor)*yi
806     $ + zi
807     $ + dz(nplx,nldx,sensor)
808    
809     c xrt = xi
810     c yrt = yi
811     c zrt = zi
812    
813     c------------------------------------------------------------------------
814     c (xPAM,yPAM,zPAM) = measured coordinates (in cm)
815     c in PAMELA reference system
816     c------------------------------------------------------------------------
817    
818     xPAM = dcoord(xrt,viewx,nldx,sensor) / 1.d4
819     yPAM = dcoord(yrt,viewy,nldy,sensor) / 1.d4
820     zPAM = ( zrt + z_mech_sensor(nplx,nldx,sensor)*1000. ) / 1.d4
821    
822 bongi 1.32 xPAM_A = 0.D0
823     yPAM_A = 0.D0
824     zPAM_A = 0.D0
825    
826     xPAM_B = 0.D0
827     yPAM_B = 0.D0
828     zPAM_B = 0.D0
829 mocchiut 1.1
830     elseif(
831     $ (icx.ne.0.and.icy.eq.0).or.
832     $ (icx.eq.0.and.icy.ne.0).or.
833     $ .false.
834     $ )then
835    
836     c------------------------------------------------------------------------
837     c (xi,yi,zi) = mechanical coordinates in the silicon sensor frame
838     c------------------------------------------------------------------------
839    
840     if(icy.ne.0)then
841     c===========================================================
842     C Y-SINGLET
843     C===========================================================
844     nplx = nply
845     nldx = nldy
846     viewx = viewy + 1
847    
848 pam-fi 1.37 xi = dcoordsi(0.5*(nstrips+1),viewx) !sensor center
849     yi = dcoordsi(stripy,viewy)
850     zi = 0.D0
851 mocchiut 1.1
852     xi_A = edgeY_d - SiDimX/2
853     yi_A = yi
854     zi_A = 0.
855    
856     xi_B = SiDimX/2 - edgeY_u
857     yi_B = yi
858     zi_B = 0.
859    
860    
861     elseif(icx.ne.0)then
862     c===========================================================
863     C X-SINGLET
864     C===========================================================
865    
866     nply = nplx
867     nldy = nldx
868     viewy = viewx - 1
869    
870 pam-fi 1.4 if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
871     $ .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
872 pam-fi 1.28 if(DEBUG.EQ.1) then
873 mocchiut 1.21 print*,'xyz_PAM (X-singlet):',
874     $ ' WARNING: false X strip: strip ',stripx
875     endif
876 pam-fi 1.4 endif
877 pam-fi 1.37
878     xi = dcoordsi(stripx,viewx)
879     yi = dcoordsi(0.5*(nstrips+1),viewy) !sensor center
880     zi = 0.D0
881 mocchiut 1.1
882     xi_A = xi
883     yi_A = edgeX_d - SiDimY/2
884     zi_A = 0.
885    
886     xi_B = xi
887     yi_B = SiDimY/2 - edgeX_u
888     zi_B = 0.
889    
890     if(viewy.eq.11)then
891     yi = yi_A
892     yi_A = yi_B
893     yi_B = yi
894     endif
895    
896    
897     else
898 pam-fi 1.28 if(DEBUG.EQ.1) then
899 mocchiut 1.21 print *,'routine xyz_PAM ---> not properly used !!!'
900     print *,'icx = ',icx
901     print *,'icy = ',icy
902     endif
903 mocchiut 1.1 goto 100
904    
905     endif
906     c------------------------------------------------------------------------
907     c (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame
908     c------------------------------------------------------------------------
909     c N.B. I convert angles from microradiants to radiants
910    
911     xrt_A = xi_A
912     $ - omega(nplx,nldx,sensor)*yi_A
913     $ + gamma(nplx,nldx,sensor)*zi_A
914     $ + dx(nplx,nldx,sensor)
915    
916     yrt_A = omega(nplx,nldx,sensor)*xi_A
917     $ + yi_A
918     $ - beta(nplx,nldx,sensor)*zi_A
919     $ + dy(nplx,nldx,sensor)
920    
921     zrt_A = -gamma(nplx,nldx,sensor)*xi_A
922     $ + beta(nplx,nldx,sensor)*yi_A
923     $ + zi_A
924     $ + dz(nplx,nldx,sensor)
925    
926     xrt_B = xi_B
927     $ - omega(nplx,nldx,sensor)*yi_B
928     $ + gamma(nplx,nldx,sensor)*zi_B
929     $ + dx(nplx,nldx,sensor)
930    
931     yrt_B = omega(nplx,nldx,sensor)*xi_B
932     $ + yi_B
933     $ - beta(nplx,nldx,sensor)*zi_B
934     $ + dy(nplx,nldx,sensor)
935    
936     zrt_B = -gamma(nplx,nldx,sensor)*xi_B
937     $ + beta(nplx,nldx,sensor)*yi_B
938     $ + zi_B
939     $ + dz(nplx,nldx,sensor)
940    
941 pam-fi 1.37
942    
943     xrt = xi
944     $ - omega(nplx,nldx,sensor)*yi
945     $ + gamma(nplx,nldx,sensor)*zi
946     $ + dx(nplx,nldx,sensor)
947    
948     yrt = omega(nplx,nldx,sensor)*xi
949     $ + yi
950     $ - beta(nplx,nldx,sensor)*zi
951     $ + dy(nplx,nldx,sensor)
952    
953     zrt = -gamma(nplx,nldx,sensor)*xi
954     $ + beta(nplx,nldx,sensor)*yi
955     $ + zi
956     $ + dz(nplx,nldx,sensor)
957    
958    
959 mocchiut 1.1
960     c xrt = xi
961     c yrt = yi
962     c zrt = zi
963    
964     c------------------------------------------------------------------------
965     c (xPAM,yPAM,zPAM) = measured coordinates (in cm)
966     c in PAMELA reference system
967     c------------------------------------------------------------------------
968    
969 pam-fi 1.37 xPAM = dcoord(xrt,viewx,nldx,sensor) / 1.d4
970     yPAM = dcoord(yrt,viewy,nldy,sensor) / 1.d4
971     zPAM = ( zrt + z_mech_sensor(nplx,nldx,sensor)*1000. ) / 1.d4
972     c$$$ xPAM = 0.D0
973     c$$$ yPAM = 0.D0
974     c$$$ zPAM = 0.D0
975 mocchiut 1.1
976     xPAM_A = dcoord(xrt_A,viewx,nldx,sensor) / 1.d4
977     yPAM_A = dcoord(yrt_A,viewy,nldy,sensor) / 1.d4
978     zPAM_A = ( zrt_A + z_mech_sensor(nplx,nldx,sensor)*1000.)/ 1.d4
979    
980     xPAM_B = dcoord(xrt_B,viewx,nldx,sensor) / 1.d4
981     yPAM_B = dcoord(yrt_B,viewy,nldy,sensor) / 1.d4
982     zPAM_B = ( zrt_B + z_mech_sensor(nplx,nldx,sensor)*1000.)/ 1.d4
983    
984    
985    
986     else
987 pam-fi 1.28 if(DEBUG.EQ.1) then
988 mocchiut 1.21 print *,'routine xyz_PAM ---> not properly used !!!'
989     print *,'icx = ',icx
990     print *,'icy = ',icy
991     endif
992 mocchiut 1.1 endif
993    
994 pam-fi 1.17
995    
996 mocchiut 1.1 100 continue
997     end
998    
999 pam-fi 1.23 ************************************************************************
1000     * Call xyz_PAM subroutine with default PFA and fill the mini2 common.
1001     * (done to be called from c/c++)
1002     ************************************************************************
1003    
1004     subroutine xyzpam(ip,icx,icy,lad,sensor,ax,ay,bfx,bfy)
1005    
1006     include 'commontracker.f'
1007     include 'level1.f'
1008     include 'common_mini_2.f'
1009     include 'common_xyzPAM.f'
1010     include 'common_mech.f'
1011     include 'calib.f'
1012    
1013     * flag to chose PFA
1014     c$$$ character*10 PFA
1015     c$$$ common/FINALPFA/PFA
1016    
1017     integer icx,icy !X-Y cluster ID
1018     integer sensor
1019     character*4 PFAx,PFAy !PFA to be used
1020     real ax,ay !X-Y geometric angle
1021     real bfx,bfy !X-Y b-field components
1022    
1023     ipx=0
1024     ipy=0
1025    
1026     c$$$ PFAx = 'COG4'!PFA
1027     c$$$ PFAy = 'COG4'!PFA
1028 pam-fi 1.26
1029 pam-fi 1.29
1030 pam-fi 1.26 if(icx.gt.nclstr1.or.icy.gt.nclstr1)then
1031     print*,'xyzpam: ***WARNING*** clusters ',icx,icy
1032 pam-fi 1.33 $ ,' do not exists (n.clusters=',nclstr1,')'
1033 pam-fi 1.26 icx = -1*icx
1034     icy = -1*icy
1035     return
1036    
1037     endif
1038 pam-fi 1.23
1039     call idtoc(pfaid,PFAx)
1040     call idtoc(pfaid,PFAy)
1041    
1042    
1043     if(icx.ne.0.and.icy.ne.0)then
1044    
1045     ipx=npl(VIEW(icx))
1046     ipy=npl(VIEW(icy))
1047 pam-fi 1.26
1048     if( (nplanes-ipx+1).ne.ip )then
1049     print*,'xyzpam: ***WARNING*** cluster ',icx
1050     $ ,' does not belong to plane: ',ip
1051     icx = -1*icx
1052     return
1053     endif
1054     if( (nplanes-ipy+1).ne.ip )then
1055     print*,'xyzpam: ***WARNING*** cluster ',icy
1056     $ ,' does not belong to plane: ',ip
1057     icy = -1*icy
1058     return
1059     endif
1060    
1061     call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
1062 pam-fi 1.23
1063     xgood(ip) = 1.
1064     ygood(ip) = 1.
1065     resx(ip) = resxPAM
1066     resy(ip) = resyPAM
1067    
1068     xm(ip) = xPAM
1069     ym(ip) = yPAM
1070     zm(ip) = zPAM
1071 bongi 1.32 xm_A(ip) = 0.D0
1072     ym_A(ip) = 0.D0
1073     xm_B(ip) = 0.D0
1074     ym_B(ip) = 0.D0
1075 pam-fi 1.23
1076     c zv(ip) = zPAM
1077    
1078     elseif(icx.eq.0.and.icy.ne.0)then
1079    
1080     ipy=npl(VIEW(icy))
1081 pam-fi 1.26 if( (nplanes-ipy+1).ne.ip )then
1082     print*,'xyzpam: ***WARNING*** cluster ',icy
1083     $ ,' does not belong to plane: ',ip
1084     icy = -1*icy
1085     return
1086     endif
1087 pam-fi 1.23
1088 pam-fi 1.26 call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
1089    
1090 pam-fi 1.23 xgood(ip) = 0.
1091     ygood(ip) = 1.
1092     resx(ip) = 1000.
1093     resy(ip) = resyPAM
1094    
1095 bongi 1.35 c$$$ xm(ip) = -100.
1096     c$$$ ym(ip) = -100.
1097     c$$$ zm(ip) = (zPAM_A+zPAM_B)/2.
1098 pam-fi 1.37 xm(ip) = xPAM
1099     ym(ip) = yPAM
1100     zm(ip) = zPAM
1101 pam-fi 1.23 xm_A(ip) = xPAM_A
1102     ym_A(ip) = yPAM_A
1103     xm_B(ip) = xPAM_B
1104     ym_B(ip) = yPAM_B
1105    
1106     c zv(ip) = (zPAM_A+zPAM_B)/2.
1107    
1108     elseif(icx.ne.0.and.icy.eq.0)then
1109    
1110     ipx=npl(VIEW(icx))
1111 pam-fi 1.26
1112     if( (nplanes-ipx+1).ne.ip )then
1113     print*,'xyzpam: ***WARNING*** cluster ',icx
1114     $ ,' does not belong to plane: ',ip
1115     icx = -1*icx
1116     return
1117     endif
1118 pam-fi 1.23
1119 pam-fi 1.26 call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
1120    
1121 pam-fi 1.23 xgood(ip) = 1.
1122     ygood(ip) = 0.
1123     resx(ip) = resxPAM
1124     resy(ip) = 1000.
1125    
1126 bongi 1.35 c$$$ xm(ip) = -100.
1127     c$$$ ym(ip) = -100.
1128     c$$$ zm(ip) = (zPAM_A+zPAM_B)/2.
1129 pam-fi 1.37 xm(ip) = xPAM
1130     ym(ip) = yPAM
1131     zm(ip) = zPAM
1132 pam-fi 1.23 xm_A(ip) = xPAM_A
1133     ym_A(ip) = yPAM_A
1134     xm_B(ip) = xPAM_B
1135     ym_B(ip) = yPAM_B
1136 pam-fi 1.37
1137 pam-fi 1.23 c zv(ip) = (zPAM_A+zPAM_B)/2.
1138    
1139     else
1140    
1141     il = 2
1142     if(lad.ne.0)il=lad
1143     is = 1
1144     if(sensor.ne.0)is=sensor
1145    
1146     xgood(ip) = 0.
1147     ygood(ip) = 0.
1148     resx(ip) = 1000.
1149     resy(ip) = 1000.
1150    
1151     xm(ip) = -100.
1152     ym(ip) = -100.
1153     zm(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
1154     xm_A(ip) = 0.
1155     ym_A(ip) = 0.
1156     xm_B(ip) = 0.
1157     ym_B(ip) = 0.
1158    
1159     c zv(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
1160    
1161     endif
1162    
1163 pam-fi 1.28 if(DEBUG.EQ.1)then
1164 pam-fi 1.23 22222 format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)
1165     write(*,22222)ip,zm(ip),xm(ip),ym(ip)
1166     $ ,xm_A(ip),ym_A(ip),xm_B(ip),ym_B(ip)
1167     $ ,xgood(ip),ygood(ip),resx(ip),resy(ip)
1168     endif
1169     end
1170 mocchiut 1.1
1171     ********************************************************************************
1172     ********************************************************************************
1173     ********************************************************************************
1174     *
1175     * The function distance_to(XP,YP) should be used after
1176     * a call to the xyz_PAM routine and it evaluate the
1177     * NORMALIZED distance (PROJECTED on the XY plane) between
1178     * the point (XP,YP), argument of the function,
1179     * and:
1180     *
1181     * - the point (xPAM,yPAM,zPAM), in the case of a COUPLE
1182     * or
1183     * - the segment (xPAM_A,yPAM_A,zPAM_A)-(xPAM_B,yPAM_B,zPAM_B),
1184     * in the case of a SINGLET.
1185     *
1186     * ( The routine xyz_PAM fills the common defined in "common_xyzPAM.f",
1187     * which stores the coordinates of the couple/singlet )
1188     *
1189     ********************************************************************************
1190    
1191 bongi 1.32 real function distance_to(rXPP,rYPP)
1192 mocchiut 1.1
1193     include 'common_xyzPAM.f'
1194    
1195     * -----------------------------------
1196     * it computes the normalized distance
1197     * ( i.e. distance/resolution )
1198     * -----------------------------------
1199    
1200 bongi 1.32 real rXPP,rYPP
1201     double precision XPP,YPP
1202 mocchiut 1.1 double precision distance,RE
1203     double precision BETA,ALFA,xmi,ymi
1204    
1205 bongi 1.32 XPP=DBLE(rXPP)
1206     YPP=DBLE(rYPP)
1207    
1208 mocchiut 1.1 * ----------------------
1209     if (
1210 pam-fi 1.37 c$$$ + xPAM.eq.0.and.
1211     c$$$ + yPAM.eq.0.and.
1212     c$$$ + zPAM.eq.0.and.
1213 mocchiut 1.1 + xPAM_A.ne.0.and.
1214     + yPAM_A.ne.0.and.
1215     + zPAM_A.ne.0.and.
1216     + xPAM_B.ne.0.and.
1217     + yPAM_B.ne.0.and.
1218     + zPAM_B.ne.0.and.
1219     + .true.)then
1220     * -----------------------
1221     * DISTANCE TO --- SINGLET
1222     * -----------------------
1223     if(abs(sngl(xPAM_B-xPAM_A)).lt.abs(sngl(yPAM_B-yPAM_A)))then
1224     * |||---------- X CLUSTER
1225    
1226     BETA = (xPAM_B-xPAM_A)/(yPAM_B-yPAM_A)
1227     ALFA = xPAM_A - BETA * yPAM_A
1228    
1229     ymi = ( YPP + BETA*XPP - BETA*ALFA )/(1+BETA**2)
1230     if(ymi.lt.dmin1(yPAM_A,yPAM_B))ymi=dmin1(yPAM_A,yPAM_B)
1231     if(ymi.gt.dmax1(yPAM_A,yPAM_B))ymi=dmax1(yPAM_A,yPAM_B)
1232     xmi = ALFA + BETA * ymi
1233     RE = resxPAM
1234    
1235     else
1236     * |||---------- Y CLUSTER
1237    
1238     BETA = (yPAM_B-yPAM_A)/(xPAM_B-xPAM_A)
1239     ALFA = yPAM_A - BETA * xPAM_A
1240    
1241     xmi = ( XPP + BETA*YPP - BETA*ALFA )/(1+BETA**2)
1242     if(xmi.lt.dmin1(xPAM_A,xPAM_B))xmi=dmin1(xPAM_A,xPAM_B)
1243     if(xmi.gt.dmax1(xPAM_A,xPAM_B))xmi=dmax1(xPAM_A,xPAM_B)
1244     ymi = ALFA + BETA * xmi
1245     RE = resyPAM
1246    
1247     endif
1248    
1249     distance=
1250 pam-fi 1.20 $ ((xmi-XPP)**2+(ymi-YPP)**2)!QUIQUI
1251     cc $ ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2
1252 mocchiut 1.1 distance=dsqrt(distance)
1253    
1254    
1255    
1256     * ----------------------
1257     elseif(
1258 pam-fi 1.37 c$$$ + xPAM.ne.0.and.
1259     c$$$ + yPAM.ne.0.and.
1260     c$$$ + zPAM.ne.0.and.
1261 mocchiut 1.1 + xPAM_A.eq.0.and.
1262     + yPAM_A.eq.0.and.
1263     + zPAM_A.eq.0.and.
1264     + xPAM_B.eq.0.and.
1265     + yPAM_B.eq.0.and.
1266     + zPAM_B.eq.0.and.
1267     + .true.)then
1268     * ----------------------
1269     * DISTANCE TO --- COUPLE
1270     * ----------------------
1271    
1272     distance=
1273 pam-fi 1.20 $ ((xPAM-XPP))**2 !QUIQUI
1274     $ +
1275     $ ((yPAM-YPP))**2
1276     c$$$ $ ((xPAM-XPP)/resxPAM)**2
1277     c$$$ $ +
1278     c$$$ $ ((yPAM-YPP)/resyPAM)**2
1279 mocchiut 1.1 distance=dsqrt(distance)
1280    
1281    
1282     else
1283    
1284     endif
1285    
1286     distance_to = sngl(distance)
1287    
1288     return
1289     end
1290    
1291     ********************************************************************************
1292     ********************************************************************************
1293     ********************************************************************************
1294     ********************************************************************************
1295    
1296     subroutine whichsensor(nplPAM,xPAM,yPAM,ladder,sensor)
1297     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
1298     * Given the plane (1-6 from BOTTOM to TOP!!) and the (xPAM,yPAM)
1299     * coordinates (in the PAMELA reference system), it returns
1300     * the ladder and the sensor which the point belongs to.
1301     *
1302     * The method to assign a point to a sensor consists in
1303     * - calculating the sum of the distances between the point
1304     * and the sensor edges
1305     * - requiring that it is less-equal than (SiDimX+SiDimY)
1306     *
1307     * NB -- SiDimX and SiDimY are not the dimentions of the SENSITIVE volume
1308     * but of the whole silicon sensor
1309     *
1310     * CONVENTION:
1311     * - sensor 1 is the one closest to the hybrid
1312     * - ladder 1 is the first to be read out (strips from 1 to 1024)
1313     *
1314     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
1315     include 'commontracker.f'
1316     include 'common_align.f'
1317    
1318     integer ladder,sensor,viewx,viewy
1319     real c1(4),c2(4),c3(4)
1320     data c1/1.,0.,0.,1./
1321     data c2/1.,-1.,-1.,1./
1322     data c3/1.,1.,0.,0./
1323 bongi 1.32 double precision yvvv,xvvv
1324 mocchiut 1.1 double precision xi,yi,zi
1325     double precision xrt,yrt,zrt
1326     real AA,BB
1327 bongi 1.32 double precision yvv(4),xvv(4)
1328 mocchiut 1.1
1329     * tollerance to consider the track inside the sensitive area
1330     real ptoll
1331     data ptoll/150./ !um
1332    
1333 bongi 1.32 external nviewx,nviewy,dcoordsi,dcoord
1334 mocchiut 1.1
1335     nplpt = nplPAM !plane
1336     viewx = nviewx(nplpt)
1337     viewy = nviewy(nplpt)
1338    
1339     do il=1,nladders_view
1340     do is=1,2
1341    
1342     do iv=1,4 !loop on sensor vertexes
1343     stripx = (il-c1(iv))*1024 + c1(iv) + c2(iv)*3
1344     stripy = (il-c3(iv))*1024 + c3(iv)
1345     c------------------------------------------------------------------------
1346     c (xi,yi,zi) = mechanical coordinates in the silicon sensor frame
1347     c------------------------------------------------------------------------
1348 bongi 1.32 xi = dcoordsi(stripx,viewx)
1349     yi = dcoordsi(stripy,viewy)
1350     zi = 0.D0
1351 mocchiut 1.1 c------------------------------------------------------------------------
1352     c (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame
1353     c------------------------------------------------------------------------
1354     c N.B. I convert angles from microradiants to radiants
1355     xrt = xi
1356     $ - omega(nplpt,il,is)*yi
1357     $ + gamma(nplpt,il,is)*zi
1358     $ + dx(nplpt,il,is)
1359     yrt = omega(nplpt,il,is)*xi
1360     $ + yi
1361     $ - beta(nplpt,il,is)*zi
1362     $ + dy(nplpt,il,is)
1363     zrt = -gamma(nplpt,il,is)*xi
1364     $ + beta(nplpt,il,is)*yi
1365     $ + zi
1366     $ + dz(nplpt,il,is)
1367     c------------------------------------------------------------------------
1368     c measured coordinates (in cm) in PAMELA reference system
1369     c------------------------------------------------------------------------
1370     yvvv = dcoord(yrt,viewy,il,is) / 1.d4
1371     xvvv = dcoord(xrt,viewx,il,is) / 1.d4
1372    
1373     yvv(iv)=sngl(yvvv)
1374     xvv(iv)=sngl(xvvv)
1375     enddo !end loop on sensor vertexes
1376    
1377     dtot=0.
1378     do iside=1,4,2 !loop on sensor edges X
1379     iv1=iside
1380     iv2=mod(iside,4)+1
1381     * straight line passing trhough two consecutive vertexes
1382     AA = (yvv(iv1)-yvv(iv2))/(xvv(iv1)-xvv(iv2))
1383     BB = yvv(iv1) - AA*xvv(iv1)
1384     * point along the straight line closer to the track
1385     xoo = (xPAM+AA*yPAM-AA*BB)/(1+AA**2)
1386     yoo = AA*xoo + BB
1387     * sum of the distances
1388     dtot = dtot +
1389     $ sqrt((xPAM-xoo)**2+(yPAM-yoo)**2)
1390     enddo !end loop on sensor edges
1391     do iside=2,4,2 !loop on sensor edges Y
1392     iv1=iside
1393     iv2=mod(iside,4)+1
1394     * straight line passing trhough two consecutive vertexes
1395     AA = (xvv(iv1)-xvv(iv2))/(yvv(iv1)-yvv(iv2))
1396     BB = xvv(iv1) - AA*yvv(iv1)
1397     * point along the straight line closer to the track
1398     yoo = (yPAM+AA*xPAM-AA*BB)/(1+AA**2)
1399     xoo = AA*yoo + BB
1400     * sum of the distances
1401     dtot = dtot +
1402     $ sqrt((xPAM-xoo)**2+(yPAM-yoo)**2)
1403     enddo !end loop on sensor edges
1404    
1405    
1406     * half-perimeter of sensitive area
1407     Perim =
1408     $ SiDimX - edgeX_l - edgeX_r
1409     $ +SiDimY - edgeY_l - edgeY_r
1410     Perim = (Perim + ptoll)/1.e4
1411     if(dtot.le.Perim)goto 100
1412    
1413    
1414     enddo
1415     enddo
1416    
1417     ladder = 0
1418     sensor = 0
1419     goto 200
1420    
1421     100 continue
1422     ladder = il
1423     sensor = is
1424    
1425    
1426     200 return
1427     end
1428    
1429    
1430    
1431     *************************************************************************
1432    
1433     subroutine reverse(v,n,temp) !invert the order of the components of v(n) vector
1434    
1435     implicit double precision (A-H,O-Z)
1436    
1437     dimension v(*)
1438     dimension temp(*)
1439     integer i,n
1440    
1441     do i=1,n
1442     temp(i)=v(n+1-i)
1443     enddo
1444    
1445     do i=1,n
1446     v(i)=temp(i)
1447     enddo
1448    
1449     return
1450     end
1451    
1452     *************************************************************************
1453     *************************************************************************
1454     *************************************************************************
1455     *************************************************************************
1456     *************************************************************************
1457     *************************************************************************
1458     *************************************************************************
1459     *************************************************************************
1460     *************************************************************************
1461     *************************************************************************
1462     *************************************************************************
1463     *************************************************************************
1464     *************************************************************************
1465     *************************************************************************
1466     *************************************************************************
1467     integer function ip_cp(id)
1468     *
1469     * given the couple id,
1470     * it returns the plane number
1471     *
1472     include 'commontracker.f'
1473 pam-fi 1.9 include 'level1.f'
1474 mocchiut 1.1 c include 'common_analysis.f'
1475     include 'common_momanhough.f'
1476    
1477     ip_cp=0
1478     ncpp=0
1479     do ip=1,nplanes
1480     ncpp=ncpp+ncp_plane(ip)
1481     if(ncpp.ge.abs(id))then
1482     ip_cp=ip
1483     goto 100
1484     endif
1485     enddo
1486     100 continue
1487     return
1488     end
1489    
1490    
1491     integer function is_cp(id)
1492     *
1493     * given the couple id,
1494     * it returns the sensor number
1495     *
1496     is_cp=0
1497     if(id.lt.0)is_cp=1
1498     if(id.gt.0)is_cp=2
1499    
1500     return
1501     end
1502    
1503    
1504     integer function icp_cp(id)
1505     *
1506     * given the couple id,
1507     * it returns the id number ON THE PLANE
1508     *
1509     include 'commontracker.f'
1510 pam-fi 1.9 include 'level1.f'
1511 mocchiut 1.1 c include 'common_analysis.f'
1512     include 'common_momanhough.f'
1513    
1514     icp_cp=0
1515    
1516     ncpp=0
1517     do ip=1,nplanes
1518     ncppold=ncpp
1519     ncpp=ncpp+ncp_plane(ip)
1520     if(ncpp.ge.abs(id))then
1521     icp_cp=abs(id)-ncppold
1522     goto 100
1523     endif
1524     enddo
1525     100 continue
1526     return
1527     end
1528    
1529    
1530    
1531     integer function id_cp(ip,icp,is)
1532     *
1533     * given a plane, a couple and the sensor
1534     * it returns the absolute couple id
1535     * negative if sensor =1
1536     * positive if sensor =2
1537     *
1538     include 'commontracker.f'
1539 pam-fi 1.9 include 'level1.f'
1540 mocchiut 1.1 c include 'calib.f'
1541     c include 'level1.f'
1542     c include 'common_analysis.f'
1543     include 'common_momanhough.f'
1544    
1545     id_cp=0
1546    
1547     if(ip.gt.1)then
1548     do i=1,ip-1
1549     id_cp = id_cp + ncp_plane(i)
1550     enddo
1551     endif
1552    
1553     id_cp = id_cp + icp
1554    
1555     if(is.eq.1) id_cp = -id_cp
1556    
1557     return
1558     end
1559    
1560    
1561    
1562    
1563     *************************************************************************
1564     *************************************************************************
1565     *************************************************************************
1566     *************************************************************************
1567     *************************************************************************
1568     *************************************************************************
1569    
1570    
1571     ***************************************************
1572     * *
1573     * *
1574     * *
1575     * *
1576     * *
1577     * *
1578     **************************************************
1579    
1580     subroutine cl_to_couples(iflag)
1581    
1582     include 'commontracker.f'
1583 pam-fi 1.9 include 'level1.f'
1584 mocchiut 1.1 include 'common_momanhough.f'
1585 pam-fi 1.9 c include 'momanhough_init.f'
1586 mocchiut 1.1 include 'calib.f'
1587 pam-fi 1.9 c include 'level1.f'
1588 mocchiut 1.1
1589     * output flag
1590     * --------------
1591     * 0 = good event
1592     * 1 = bad event
1593     * --------------
1594     integer iflag
1595    
1596 pam-fi 1.6 integer badseed,badclx,badcly
1597 mocchiut 1.40
1598     iflag = iflag
1599 pam-fi 1.28 if(DEBUG.EQ.1)print*,'cl_to_couples:'
1600 mocchiut 1.1
1601 pam-fi 1.37 cc if(RECOVER_SINGLETS.and..not.SECOND_SEARCH)goto 80
1602    
1603 mocchiut 1.1 * init variables
1604     do ip=1,nplanes
1605     do ico=1,ncouplemax
1606     clx(ip,ico)=0
1607     cly(ip,ico)=0
1608     enddo
1609     ncp_plane(ip)=0
1610     do icl=1,nclstrmax_level2
1611     cls(ip,icl)=1
1612     enddo
1613     ncls(ip)=0
1614     enddo
1615     do icl=1,nclstrmax_level2
1616 pam-fi 1.37 cl_single(icl) = 1 !all are single
1617     cl_good(icl) = 0 !all are bad
1618 pam-fi 1.9 enddo
1619     do iv=1,nviews
1620     ncl_view(iv) = 0
1621     mask_view(iv) = 0 !all included
1622 mocchiut 1.1 enddo
1623    
1624 pam-fi 1.9 * count number of cluster per view
1625     do icl=1,nclstr1
1626     ncl_view(VIEW(icl)) = ncl_view(VIEW(icl)) + 1
1627     enddo
1628     * mask views with too many clusters
1629     do iv=1,nviews
1630 pam-fi 1.14 if( ncl_view(iv).gt. nclusterlimit)then
1631 pam-fi 1.26 c mask_view(iv) = 1
1632     mask_view(iv) = mask_view(iv) + 2**0
1633 pam-fi 1.28 if(DEBUG.EQ.1)
1634     $ print*,' * WARNING * cl_to_couple: n.clusters > '
1635     $ ,nclusterlimit,' on view ', iv,' --> masked!'
1636 pam-fi 1.9 endif
1637     enddo
1638    
1639    
1640 mocchiut 1.1 * start association
1641     ncouples=0
1642     do icx=1,nclstr1 !loop on cluster (X)
1643     if(mod(VIEW(icx),2).eq.1)goto 10
1644    
1645 pam-fi 1.37 if(cl_used(icx).ne.0)goto 10
1646    
1647 mocchiut 1.1 * ----------------------------------------------------
1648 pam-fi 1.9 * jump masked views (X VIEW)
1649     * ----------------------------------------------------
1650     if( mask_view(VIEW(icx)).ne.0 ) goto 10
1651     * ----------------------------------------------------
1652 mocchiut 1.1 * cut on charge (X VIEW)
1653 pam-fi 1.4 * ----------------------------------------------------
1654 pam-fi 1.37 if(sgnl(icx)/mip(VIEW(icx),LADDER(icx)).lt.dedx_x_min)then
1655 mocchiut 1.1 cl_single(icx)=0
1656     goto 10
1657     endif
1658 pam-fi 1.4 * ----------------------------------------------------
1659 mocchiut 1.1 * cut BAD (X VIEW)
1660 pam-fi 1.4 * ----------------------------------------------------
1661 mocchiut 1.1 badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx)))
1662     ifirst=INDSTART(icx)
1663     if(icx.ne.nclstr1) then
1664     ilast=INDSTART(icx+1)-1
1665     else
1666     ilast=TOTCLLENGTH
1667     endif
1668 pam-fi 1.6 badclx=badseed
1669 mocchiut 1.1 do igood=-ngoodstr,ngoodstr
1670     ibad=1
1671     if((INDMAX(icx)+igood).gt.ifirst.and.
1672     $ (INDMAX(icx)+igood).lt.ilast.and.
1673     $ .true.)then
1674     ibad=BAD(VIEW(icx),
1675     $ nvk(MAXS(icx)+igood),
1676     $ nst(MAXS(icx)+igood))
1677     endif
1678 pam-fi 1.6 badclx=badclx*ibad
1679 mocchiut 1.1 enddo
1680 pam-fi 1.4 * ----------------------------------------------------
1681     * >>> eliminato il taglio sulle BAD <<<
1682     * ----------------------------------------------------
1683 mocchiut 1.1 c if(badcl.eq.0)then
1684     c cl_single(icx)=0
1685     c goto 10
1686     c endif
1687     * ----------------------------------------------------
1688    
1689     cl_good(icx)=1
1690     nplx=npl(VIEW(icx))
1691     nldx=nld(MAXS(icx),VIEW(icx))
1692    
1693     do icy=1,nclstr1 !loop on cluster (Y)
1694     if(mod(VIEW(icy),2).eq.0)goto 20
1695 pam-fi 1.37
1696     if(cl_used(icx).ne.0)goto 20
1697 mocchiut 1.1
1698     * ----------------------------------------------------
1699 pam-fi 1.9 * jump masked views (Y VIEW)
1700     * ----------------------------------------------------
1701     if( mask_view(VIEW(icy)).ne.0 ) goto 20
1702    
1703     * ----------------------------------------------------
1704 mocchiut 1.1 * cut on charge (Y VIEW)
1705 pam-fi 1.4 * ----------------------------------------------------
1706 pam-fi 1.37 if(sgnl(icy)/mip(VIEW(icy),LADDER(icy)) .lt.dedx_y_min)then
1707 mocchiut 1.1 cl_single(icy)=0
1708     goto 20
1709     endif
1710 pam-fi 1.4 * ----------------------------------------------------
1711 mocchiut 1.1 * cut BAD (Y VIEW)
1712 pam-fi 1.4 * ----------------------------------------------------
1713 mocchiut 1.1 badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy)))
1714     ifirst=INDSTART(icy)
1715     if(icy.ne.nclstr1) then
1716     ilast=INDSTART(icy+1)-1
1717     else
1718     ilast=TOTCLLENGTH
1719     endif
1720 pam-fi 1.6 badcly=badseed
1721 mocchiut 1.1 do igood=-ngoodstr,ngoodstr
1722     ibad=1
1723     if((INDMAX(icy)+igood).gt.ifirst.and.
1724     $ (INDMAX(icy)+igood).lt.ilast.and.
1725     $ .true.)
1726     $ ibad=BAD(VIEW(icy),
1727     $ nvk(MAXS(icy)+igood),
1728     $ nst(MAXS(icy)+igood))
1729 pam-fi 1.6 badcly=badcly*ibad
1730 mocchiut 1.1 enddo
1731 pam-fi 1.4 * ----------------------------------------------------
1732     * >>> eliminato il taglio sulle BAD <<<
1733     * ----------------------------------------------------
1734 mocchiut 1.1 c if(badcl.eq.0)then
1735     c cl_single(icy)=0
1736     c goto 20
1737     c endif
1738     * ----------------------------------------------------
1739    
1740     cl_good(icy)=1
1741     nply=npl(VIEW(icy))
1742     nldy=nld(MAXS(icy),VIEW(icy))
1743    
1744     * ----------------------------------------------
1745     * CONDITION TO FORM A COUPLE
1746     * ----------------------------------------------
1747     * geometrical consistency (same plane and ladder)
1748     if(nply.eq.nplx.and.nldy.eq.nldx)then
1749     * charge correlation
1750 pam-fi 1.4 * (modified to be applied only below saturation... obviously)
1751    
1752 pam-fi 1.18 if( .not.(sgnl(icy).gt.chsaty.and.sgnl(icx).gt.chsatx)
1753 pam-fi 1.6 $ .and.
1754 pam-fi 1.18 $ .not.(sgnl(icy).lt.chmipy.and.sgnl(icx).lt.chmipx)
1755 pam-fi 1.6 $ .and.
1756     $ (badclx.eq.1.and.badcly.eq.1)
1757     $ .and.
1758     $ .true.)then
1759    
1760 pam-fi 1.18 ddd=(sgnl(icy)
1761     $ -kch(nplx,nldx)*sgnl(icx)-cch(nplx,nldx))
1762 pam-fi 1.6 ddd=ddd/sqrt(kch(nplx,nldx)**2+1)
1763    
1764     c cut = chcut * sch(nplx,nldx)
1765    
1766 pam-fi 1.18 sss=(kch(nplx,nldx)*sgnl(icy)+sgnl(icx)
1767 pam-fi 1.6 $ -kch(nplx,nldx)*cch(nplx,nldx))
1768     sss=sss/sqrt(kch(nplx,nldx)**2+1)
1769     cut = chcut * (16 + sss/50.)
1770    
1771     if(abs(ddd).gt.cut)then
1772     goto 20 !charge not consistent
1773     endif
1774     endif
1775 pam-fi 1.14
1776 pam-fi 1.15 if(ncp_plane(nplx).gt.ncouplemax)then
1777 pam-fi 1.28 if(verbose.eq.1)print*,
1778 pam-fi 1.5 $ '** warning ** number of identified '//
1779     $ 'couples on plane ',nplx,
1780     $ 'exceeds vector dimention '
1781 pam-fi 1.14 $ ,'( ',ncouplemax,' ) --> masked!'
1782 pam-fi 1.26 c mask_view(nviewx(nplx)) = 2
1783     c mask_view(nviewy(nply)) = 2
1784 mocchiut 1.40 mask_view(nviewx(nplx))= mask_view(nviewx(nplx))+ 2**1
1785     mask_view(nviewy(nply))= mask_view(nviewy(nply))+ 2**1
1786 pam-fi 1.15 goto 10
1787 pam-fi 1.5 endif
1788 pam-fi 1.15
1789     * ------------------> COUPLE <------------------
1790     ncp_plane(nplx) = ncp_plane(nplx) + 1
1791     clx(nplx,ncp_plane(nplx))=icx
1792     cly(nply,ncp_plane(nplx))=icy
1793     cl_single(icx)=0
1794     cl_single(icy)=0
1795 pam-fi 1.14 * ----------------------------------------------
1796    
1797 mocchiut 1.1 endif
1798    
1799     20 continue
1800     enddo !end loop on clusters(Y)
1801    
1802     10 continue
1803     enddo !end loop on clusters(X)
1804    
1805     do icl=1,nclstr1
1806     if(cl_single(icl).eq.1)then
1807     ip=npl(VIEW(icl))
1808     ncls(ip)=ncls(ip)+1
1809     cls(ip,ncls(ip))=icl
1810     endif
1811     enddo
1812 pam-fi 1.37
1813 mocchiut 1.40 c 80 continue
1814     continue
1815 mocchiut 1.1
1816    
1817 pam-fi 1.28 if(DEBUG.EQ.1)then
1818 mocchiut 1.1 print*,'clusters ',nclstr1
1819     print*,'good ',(cl_good(i),i=1,nclstr1)
1820 pam-fi 1.37 print*,'used ',(cl_used(i),i=1,nclstr1)
1821 pam-fi 1.28 print*,'singlets',(cl_single(i),i=1,nclstr1)
1822 mocchiut 1.1 print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)
1823     endif
1824 pam-fi 1.37
1825    
1826     if(.not.RECOVER_SINGLETS)goto 81
1827    
1828     * ////////////////////////////////////////////////
1829     * PATCH to recover events with less than 3 couples
1830     * ////////////////////////////////////////////////
1831     * loop over singlet and create "fake" couples
1832     * (with clx=0 or cly=0)
1833     *
1834    
1835     if(DEBUG.EQ.1)
1836     $ print*,'>>> Recover singlets '
1837     $ ,'(creates fake couples) <<<'
1838     do icl=1,nclstr1
1839     if(
1840     $ cl_single(icl).eq.1.and.
1841     $ cl_used(icl).eq.0.and.
1842     $ .true.)then
1843     * ----------------------------------------------------
1844     * jump masked views
1845     * ----------------------------------------------------
1846     if( mask_view(VIEW(icl)).ne.0 ) goto 21
1847     if(mod(VIEW(icl),2).eq.0)then !=== X views
1848     * ----------------------------------------------------
1849     * cut on charge (X VIEW)
1850     * ----------------------------------------------------
1851     if(sgnl(icl).lt.dedx_x_min) goto 21
1852    
1853     nplx=npl(VIEW(icl))
1854     * ------------------> (FAKE) COUPLE <-----------------
1855     ncp_plane(nplx) = ncp_plane(nplx) + 1
1856     clx(nplx,ncp_plane(nplx))=icl
1857     cly(nplx,ncp_plane(nplx))=0
1858     c$$$ cl_single(icl)=0! I leave the cluster tagged as singlet!!!
1859     * ----------------------------------------------------
1860    
1861     else !=== Y views
1862     * ----------------------------------------------------
1863     * cut on charge (Y VIEW)
1864     * ----------------------------------------------------
1865     if(sgnl(icl).lt.dedx_y_min) goto 21
1866    
1867     nply=npl(VIEW(icl))
1868     * ------------------> (FAKE) COUPLE <-----------------
1869     ncp_plane(nply) = ncp_plane(nply) + 1
1870     clx(nply,ncp_plane(nply))=0
1871     cly(nply,ncp_plane(nply))=icl
1872     c$$$ cl_single(icl)=0! I leave the cluster tagged as singlet!!!
1873     * ----------------------------------------------------
1874    
1875     endif
1876     endif !end "single" condition
1877     21 continue
1878     enddo !end loop over clusters
1879    
1880     if(DEBUG.EQ.1)
1881     $ print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)
1882    
1883    
1884     81 continue
1885 mocchiut 1.1
1886 pam-fi 1.37 ncp_tot=0
1887     do ip=1,NPLANES
1888 pam-fi 1.10 ncp_tot = ncp_tot + ncp_plane(ip)
1889 mocchiut 1.1 enddo
1890 pam-fi 1.37 if(DEBUG.EQ.1)
1891     $ print*,'n.couple tot: ',ncp_tot
1892 mocchiut 1.1
1893     return
1894     end
1895    
1896     ***************************************************
1897     * *
1898     * *
1899     * *
1900     * *
1901     * *
1902     * *
1903     **************************************************
1904    
1905     subroutine cp_to_doubtrip(iflag)
1906    
1907     include 'commontracker.f'
1908 pam-fi 1.9 include 'level1.f'
1909 mocchiut 1.1 include 'common_momanhough.f'
1910     include 'common_xyzPAM.f'
1911     include 'common_mini_2.f'
1912     include 'calib.f'
1913    
1914    
1915     * output flag
1916     * --------------
1917     * 0 = good event
1918     * 1 = bad event
1919     * --------------
1920     integer iflag
1921    
1922    
1923     * -----------------------------
1924     * DOUBLETS/TRIPLETS coordinates
1925     c double precision xm1,ym1,zm1
1926     c double precision xm2,ym2,zm2
1927     c double precision xm3,ym3,zm3
1928    
1929     real xm1,ym1,zm1
1930     real xm2,ym2,zm2
1931     real xm3,ym3,zm3
1932     * -----------------------------
1933     * variable needed for tricircle:
1934     real xp(3),zp(3)!TRIPLETS coordinates, to find a circle
1935     EQUIVALENCE (xm1,xp(1))
1936     EQUIVALENCE (xm2,xp(2))
1937     EQUIVALENCE (xm3,xp(3))
1938     EQUIVALENCE (zm1,zp(1))
1939     EQUIVALENCE (zm2,zp(2))
1940     EQUIVALENCE (zm3,zp(3))
1941     real angp(3),resp(3),chi
1942     real xc,zc,radius
1943     * -----------------------------
1944    
1945 pam-fi 1.28 if(DEBUG.EQ.1)print*,'cp_to_doubtrip:'
1946 mocchiut 1.1
1947 pam-fi 1.14 * --------------------------------------------
1948     * put a limit to the maximum number of couples
1949     * per plane, in order to apply hough transform
1950     * (couples recovered during track refinement)
1951     * --------------------------------------------
1952     do ip=1,nplanes
1953     if(ncp_plane(ip).gt.ncouplelimit)then
1954 pam-fi 1.26 mask_view(nviewx(ip)) = mask_view(nviewx(ip)) + 2**7
1955     mask_view(nviewy(ip)) = mask_view(nviewy(ip)) + 2**7
1956 pam-fi 1.14 endif
1957     enddo
1958    
1959 mocchiut 1.1
1960     ndblt=0 !number of doublets
1961     ntrpt=0 !number of triplets
1962    
1963     do ip1=1,(nplanes-1) !loop on planes - COPPIA 1
1964 pam-fi 1.37 c$$$ print*,'(1) ip ',ip1
1965     c$$$ $ ,mask_view(nviewx(ip1))
1966     c$$$ $ ,mask_view(nviewy(ip1))
1967 pam-fi 1.14 if( mask_view(nviewx(ip1)).ne.0 .or.
1968     $ mask_view(nviewy(ip1)).ne.0 )goto 10 !skip plane
1969     do is1=1,2 !loop on sensors - COPPIA 1
1970 mocchiut 1.1 do icp1=1,ncp_plane(ip1) !loop on COPPIA 1
1971     icx1=clx(ip1,icp1)
1972     icy1=cly(ip1,icp1)
1973 pam-fi 1.37
1974     c$$$ print*,'(1) ip ',ip1,' icp ',icp1
1975    
1976 mocchiut 1.1 c call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)
1977 pam-fi 1.18 c call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)
1978     call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.,0.,0.)
1979 mocchiut 1.1 xm1=xPAM
1980     ym1=yPAM
1981     zm1=zPAM
1982 pam-fi 1.14
1983 mocchiut 1.1 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2
1984 pam-fi 1.37 c$$$ print*,'(2) ip ',ip2
1985     c$$$ $ ,mask_view(nviewx(ip2))
1986     c$$$ $ ,mask_view(nviewy(ip2))
1987 pam-fi 1.14 if( mask_view(nviewx(ip2)).ne.0 .or.
1988     $ mask_view(nviewy(ip2)).ne.0 )goto 20 !skip plane
1989 pam-fi 1.28 do is2=1,2 !loop on sensors -ndblt COPPIA 2
1990 mocchiut 1.1 do icp2=1,ncp_plane(ip2) !loop on COPPIA 2
1991     icx2=clx(ip2,icp2)
1992     icy2=cly(ip2,icp2)
1993 pam-fi 1.37
1994     c$$$ print*,'(2) ip ',ip2,' icp ',icp2
1995    
1996 mocchiut 1.1 c call xyz_PAM
1997     c $ (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1)
1998 pam-fi 1.18 c call xyz_PAM
1999     c $ (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1)
2000 mocchiut 1.1 call xyz_PAM
2001 pam-fi 1.18 $ (icx2,icy2,is2,PFAdef,PFAdef,0.,0.,0.,0.)
2002 mocchiut 1.1 xm2=xPAM
2003     ym2=yPAM
2004     zm2=zPAM
2005 pam-fi 1.37
2006     * ---------------------------------------------------
2007     * both couples must have a y-cluster
2008     * (condition necessary when in RECOVER_SINGLETS mode)
2009     * ---------------------------------------------------
2010     if(icy1.eq.0.or.icy2.eq.0)goto 111
2011    
2012     if(cl_used(icy1).ne.0)goto 111
2013     if(cl_used(icy2).ne.0)goto 111
2014    
2015    
2016 mocchiut 1.1 * - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2017     * track parameters on Y VIEW
2018     * (2 couples needed)
2019     * - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2020     if(ndblt.eq.ndblt_max)then
2021 pam-fi 1.28 if(verbose.eq.1)print*,
2022 mocchiut 1.1 $ '** warning ** number of identified '//
2023     $ 'doublets exceeds vector dimention '
2024     $ ,'( ',ndblt_max,' )'
2025 pam-fi 1.37 c good2=.false.
2026     c goto 880 !fill ntp and go to next event
2027 pam-fi 1.10 do iv=1,12
2028 pam-fi 1.37 c mask_view(iv) = 3
2029 pam-fi 1.26 mask_view(iv) = mask_view(iv)+ 2**2
2030 pam-fi 1.10 enddo
2031 mocchiut 1.1 iflag=1
2032     return
2033     endif
2034 pam-fi 1.37
2035    
2036     ccc print*,'<doublet> ',icp1,icp2
2037    
2038 mocchiut 1.1 ndblt = ndblt + 1
2039     * store doublet info
2040     cpyz1(ndblt)=id_cp(ip1,icp1,is1)
2041     cpyz2(ndblt)=id_cp(ip2,icp2,is2)
2042     * tg(th_yz)
2043     alfayz2(ndblt)=(ym1-ym2)/(zm1-zm2)
2044     * y0 (cm)
2045     alfayz1(ndblt)=alfayz2(ndblt)*(zini-zm1)+ym1
2046 pam-fi 1.37
2047 mocchiut 1.1 **** -----------------------------------------------****
2048     **** reject non phisical couples ****
2049     **** -----------------------------------------------****
2050 pam-fi 1.37 if(SECOND_SEARCH)goto 111
2051 mocchiut 1.1 if(
2052     $ abs(alfayz2(ndblt)).gt.alfyz2_max
2053     $ .or.
2054     $ abs(alfayz1(ndblt)).gt.alfyz1_max
2055     $ )ndblt = ndblt-1
2056    
2057 pam-fi 1.37
2058     111 continue
2059 mocchiut 1.1 * - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2060     * track parameters on Y VIEW - end
2061     * - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2062    
2063    
2064 pam-fi 1.37 if(icx1.ne.0)then
2065     if(cl_used(icx1).ne.0)goto 31
2066     endif
2067     if(icx2.ne.0)then
2068     if(cl_used(icx2).ne.0)goto 31
2069     endif
2070    
2071 pam-fi 1.14 if(ip2.eq.nplanes)goto 31 !no possible combination with 3 couples
2072    
2073 mocchiut 1.1 do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3
2074 pam-fi 1.37 c$$$ print*,'(3) ip ',ip3
2075     c$$$ $ ,mask_view(nviewx(ip3))
2076     c$$$ $ ,mask_view(nviewy(ip3))
2077 pam-fi 1.14 if( mask_view(nviewx(ip3)).ne.0 .or.
2078     $ mask_view(nviewy(ip3)).ne.0 )goto 30 !skip plane
2079 mocchiut 1.1 do is3=1,2 !loop on sensors - COPPIA 3
2080    
2081     do icp3=1,ncp_plane(ip3) !loop on COPPIA 3
2082     icx3=clx(ip3,icp3)
2083     icy3=cly(ip3,icp3)
2084 pam-fi 1.37
2085     c$$$ print*,'(3) ip ',ip3,' icp ',icp3
2086    
2087     * ---------------------------------------------------
2088     * all three couples must have a x-cluster
2089     * (condition necessary when in RECOVER_SINGLETS mode)
2090     * ---------------------------------------------------
2091     if(
2092     $ icx1.eq.0.or.
2093     $ icx2.eq.0.or.
2094     $ icx3.eq.0.or.
2095     $ .false.)goto 29
2096    
2097     if(cl_used(icx1).ne.0)goto 29
2098     if(cl_used(icx2).ne.0)goto 29
2099     if(cl_used(icx3).ne.0)goto 29
2100    
2101 mocchiut 1.1 c call xyz_PAM
2102     c $ (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1)
2103 pam-fi 1.18 c call xyz_PAM
2104     c $ (icx3,icy3,is3,PFAdef,PFAdef,0.,0.) !(1)
2105 mocchiut 1.1 call xyz_PAM
2106 pam-fi 1.18 $ (icx3,icy3,is3,PFAdef,PFAdef
2107     $ ,0.,0.,0.,0.)
2108 mocchiut 1.1 xm3=xPAM
2109     ym3=yPAM
2110     zm3=zPAM
2111 pam-fi 1.37
2112    
2113 mocchiut 1.1 * find the circle passing through the three points
2114 pam-fi 1.37 iflag_t = DEBUG
2115 mocchiut 1.1 call tricircle(3,xp,zp,angp,resp,chi
2116 pam-fi 1.37 $ ,xc,zc,radius,iflag_t)
2117 mocchiut 1.1 * the circle must intersect the reference plane
2118 pam-fi 1.37 cc if(iflag.ne.0)goto 29
2119     if(iflag_t.ne.0)then
2120     * if tricircle fails, evaluate a straight line
2121     if(DEBUG.eq.1)
2122     $ print*,'TRICIRCLE failure'
2123     $ ,' >>> straight line'
2124     radius=0.
2125     xc=0.
2126     yc=0.
2127    
2128     SZZ=0.
2129     SZX=0.
2130     SSX=0.
2131     SZ=0.
2132     S1=0.
2133     X0=0.
2134     Ax=0.
2135     BX=0.
2136     DO I=1,3
2137     XX = XP(I)
2138     SZZ=SZZ+ZP(I)*ZP(I)
2139     SZX=SZX+ZP(I)*XX
2140     SSX=SSX+XX
2141     SZ=SZ+ZP(I)
2142 mocchiut 1.40 S1=S1+1.
2143 pam-fi 1.37 ENDDO
2144     DET=SZZ*S1-SZ*SZ
2145     AX=(SZX*S1-SZ*SSX)/DET
2146     BX=(SZZ*SSX-SZX*SZ)/DET
2147     X0 = AX*ZINI+BX
2148    
2149     endif
2150    
2151     if( .not.SECOND_SEARCH.and.
2152     $ radius**2.lt.(ZINI-zc)**2)goto 29
2153    
2154 mocchiut 1.1 * - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2155     * track parameters on X VIEW
2156     * (3 couples needed)
2157     * - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2158     if(ntrpt.eq.ntrpt_max)then
2159 pam-fi 1.28 if(verbose.eq.1)print*,
2160 pam-fi 1.37 $ '** warning **'//
2161     $ ' number of identified '//
2162     $ 'triplets exceeds'//
2163     $ ' vector dimention '
2164     $ ,'( ',ntrpt_max,' )'
2165     c good2=.false.
2166     c goto 880 !fill ntp and go to next event
2167 pam-fi 1.10 do iv=1,nviews
2168 pam-fi 1.37 c mask_view(iv) = 4
2169     mask_view(iv) =
2170     $ mask_view(iv)+ 2**3
2171 pam-fi 1.10 enddo
2172 mocchiut 1.1 iflag=1
2173     return
2174     endif
2175 pam-fi 1.37
2176     ccc print*,'<triplet> ',icp1,icp2,icp3
2177    
2178 mocchiut 1.1 ntrpt = ntrpt +1
2179     * store triplet info
2180     cpxz1(ntrpt)=id_cp(ip1,icp1,is1)
2181     cpxz2(ntrpt)=id_cp(ip2,icp2,is2)
2182     cpxz3(ntrpt)=id_cp(ip3,icp3,is3)
2183    
2184 pam-fi 1.37 if(radius.ne.0.and.xc.lt.0)then
2185 mocchiut 1.1 *************POSITIVE DEFLECTION
2186 pam-fi 1.37 alfaxz1(ntrpt) = xc+sqrt(radius**2-(ZINI-zc)**2)
2187     alfaxz2(ntrpt) = (ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2)
2188     alfaxz3(ntrpt) = 1/radius
2189     else if(radius.ne.0.and.xc.ge.0)then
2190 mocchiut 1.1 *************NEGATIVE DEFLECTION
2191 pam-fi 1.37 alfaxz1(ntrpt) = xc-sqrt(radius**2-(ZINI-zc)**2)
2192     alfaxz2(ntrpt) = -(ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2)
2193     alfaxz3(ntrpt) = -1/radius
2194     else if(radius.eq.0)then
2195     *************straight fit
2196     alfaxz1(ntrpt) = X0
2197     alfaxz2(ntrpt) = AX
2198     alfaxz3(ntrpt) = 0.
2199     endif
2200    
2201     c$$$ print*,'alfaxz1 ', alfaxz1(ntrpt)
2202     c$$$ print*,'alfaxz2 ', alfaxz2(ntrpt)
2203     c$$$ print*,'alfaxz3 ', alfaxz3(ntrpt)
2204    
2205 mocchiut 1.1 **** -----------------------------------------------****
2206     **** reject non phisical triplets ****
2207     **** -----------------------------------------------****
2208 pam-fi 1.37 if(SECOND_SEARCH)goto 29
2209     if(
2210     $ abs(alfaxz2(ntrpt)).gt.
2211     $ alfxz2_max
2212     $ .or.
2213     $ abs(alfaxz1(ntrpt)).gt.
2214     $ alfxz1_max
2215     $ )ntrpt = ntrpt-1
2216    
2217 mocchiut 1.1 * - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2218     * track parameters on X VIEW - end
2219     * - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2220 pam-fi 1.37
2221     29 continue
2222 mocchiut 1.1 enddo !end loop on COPPIA 3
2223     enddo !end loop on sensors - COPPIA 3
2224 pam-fi 1.14 30 continue
2225 mocchiut 1.1 enddo !end loop on planes - COPPIA 3
2226 pam-fi 1.37
2227     31 continue
2228 mocchiut 1.40 c 1 enddo !end loop on COPPIA 2
2229     enddo !end loop on COPPIA 2
2230 mocchiut 1.1 enddo !end loop on sensors - COPPIA 2
2231 pam-fi 1.14 20 continue
2232 mocchiut 1.1 enddo !end loop on planes - COPPIA 2
2233 pam-fi 1.37
2234 mocchiut 1.40 c 11 continue
2235     continue
2236 mocchiut 1.1 enddo !end loop on COPPIA1
2237     enddo !end loop on sensors - COPPIA 1
2238 pam-fi 1.14 10 continue
2239 mocchiut 1.1 enddo !end loop on planes - COPPIA 1
2240    
2241 pam-fi 1.28 if(DEBUG.EQ.1)then
2242 mocchiut 1.1 print*,'--- doublets ',ndblt
2243     print*,'--- triplets ',ntrpt
2244     endif
2245    
2246     c goto 880 !ntp fill
2247    
2248    
2249     return
2250     end
2251    
2252    
2253    
2254     ***************************************************
2255     * *
2256     * *
2257     * *
2258     * *
2259     * *
2260     * *
2261     **************************************************
2262    
2263     subroutine doub_to_YZcloud(iflag)
2264    
2265     include 'commontracker.f'
2266 pam-fi 1.9 include 'level1.f'
2267 mocchiut 1.1 include 'common_momanhough.f'
2268 pam-fi 1.9 c include 'momanhough_init.f'
2269 mocchiut 1.1
2270    
2271     * output flag
2272     * --------------
2273     * 0 = good event
2274     * 1 = bad event
2275     * --------------
2276     integer iflag
2277    
2278     integer db_used(ndblt_max)
2279     integer db_temp(ndblt_max)
2280     integer db_all(ndblt_max) !stores db ID in each cloud
2281    
2282     integer hit_plane(nplanes)
2283    
2284     * mask for used couples
2285     integer cp_useds1(ncouplemaxtot) ! sensor 1
2286     integer cp_useds2(ncouplemaxtot) ! sensor 2
2287    
2288 pam-fi 1.28 if(DEBUG.EQ.1)print*,'doub_to_YZcloud:'
2289 mocchiut 1.1
2290     *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2291     * classification of DOUBLETS
2292     * according to distance in parameter space
2293     * (cloud = group of points (doublets) in parameter space)
2294     *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2295     do idb=1,ndblt
2296     db_used(idb)=0
2297     enddo
2298    
2299     distance=0
2300     nclouds_yz=0 !number of clouds
2301     npt_tot=0
2302 pam-fi 1.9 nloop=0
2303     90 continue
2304 mocchiut 1.1 do idb1=1,ndblt !loop (1) on DOUBLETS
2305     if(db_used(idb1).eq.1)goto 2228 !db already included in a cloud
2306    
2307     do icp=1,ncp_tot
2308     cp_useds1(icp)=0 !init
2309     cp_useds2(icp)=0 !init
2310     enddo
2311     do idb=1,ndblt
2312     db_all(idb)=0
2313     enddo
2314     if(cpyz1(idb1).gt.0)cp_useds2(cpyz1(idb1))=1
2315     if(cpyz1(idb1).lt.0)cp_useds1(-cpyz1(idb1))=1
2316     if(cpyz2(idb1).gt.0)cp_useds2(cpyz2(idb1))=1
2317     if(cpyz2(idb1).lt.0)cp_useds1(-cpyz2(idb1))=1
2318     temp1 = alfayz1(idb1)
2319     temp2 = alfayz2(idb1)
2320     npt=1 !counter of points in the cloud
2321    
2322     db_all(npt) = idb1
2323    
2324     nptloop=1
2325     db_temp(1)=idb1
2326    
2327     88 continue
2328    
2329     npv=0 !# new points inlcuded
2330     do iloop=1,nptloop
2331     idbref=db_temp(iloop) !local point of reference
2332     ccccc if(db_used(idbref).eq.1)goto 1188 !next
2333    
2334     do idb2=1,ndblt !loop (2) on DOUBLETS
2335     if(idb2.eq.idbref)goto 1118 !next doublet
2336     if(db_used(idb2).eq.1)goto 1118
2337    
2338    
2339     * doublet distance in parameter space
2340     distance=
2341     $ ((alfayz1(idbref)-alfayz1(idb2))/Dalfayz1)**2
2342 mocchiut 1.40 $ +((alfayz2(idbref)-alfayz2(idb2))/Dalfayz2)**2
2343 mocchiut 1.1 distance = sqrt(distance)
2344    
2345     if(distance.lt.cutdistyz)then
2346    
2347     if(cpyz1(idb2).gt.0)cp_useds2(cpyz1(idb2))=1
2348     if(cpyz1(idb2).lt.0)cp_useds1(-cpyz1(idb2))=1
2349     if(cpyz2(idb2).gt.0)cp_useds2(cpyz2(idb2))=1
2350     if(cpyz2(idb2).lt.0)cp_useds1(-cpyz2(idb2))=1
2351     npt = npt + 1 !counter of points in the cloud
2352    
2353     npv = npv +1
2354     db_temp(npv) = idb2
2355     db_used(idbref) = 1
2356     db_used(idb2) = 1
2357    
2358     db_all(npt) = idb2
2359    
2360     temp1 = temp1 + alfayz1(idb2)
2361     temp2 = temp2 + alfayz2(idb2)
2362     endif
2363    
2364     1118 continue
2365     enddo !end loop (2) on DOUBLETS
2366    
2367 mocchiut 1.40 c 1188 continue
2368     continue
2369 mocchiut 1.1 enddo !end loop on... bo?
2370    
2371     nptloop=npv
2372     if(nptloop.ne.0)goto 88
2373    
2374     * ------------------------------------------
2375     * stores the cloud only if
2376     * 1) it includes a minimum number of REAL couples
2377     * 1bis) it inlcudes a minimum number of doublets
2378     * 2) it is not already stored
2379     * ------------------------------------------
2380     do ip=1,nplanes
2381     hit_plane(ip)=0
2382     enddo
2383     ncpused=0
2384     do icp=1,ncp_tot
2385 pam-fi 1.37 if(
2386     $ (cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0).and.
2387     $ .true.)then
2388 mocchiut 1.1 ncpused=ncpused+1
2389     ip=ip_cp(icp)
2390     hit_plane(ip)=1
2391     endif
2392     enddo
2393     nplused=0
2394     do ip=1,nplanes
2395     nplused=nplused+ hit_plane(ip)
2396     enddo
2397 pam-fi 1.37
2398 mocchiut 1.1 if(nplused.lt.nplyz_min)goto 2228 !next doublet
2399    
2400     * ~~~~~~~~~~~~~~~~~
2401     * >>> NEW CLOUD <<<
2402    
2403     if(nclouds_yz.ge.ncloyz_max)then
2404 pam-fi 1.28 if(verbose.eq.1)print*,
2405 mocchiut 1.1 $ '** warning ** number of identified '//
2406     $ 'YZ clouds exceeds vector dimention '
2407     $ ,'( ',ncloyz_max,' )'
2408     c good2=.false.
2409     c goto 880 !fill ntp and go to next event
2410 pam-fi 1.10 do iv=1,nviews
2411 pam-fi 1.26 c mask_view(iv) = 5
2412     mask_view(iv) = mask_view(iv) + 2**4
2413 pam-fi 1.10 enddo
2414 mocchiut 1.1 iflag=1
2415     return
2416     endif
2417    
2418     nclouds_yz = nclouds_yz + 1 !increase counter
2419     alfayz1_av(nclouds_yz) = temp1/npt !store average parameter
2420     alfayz2_av(nclouds_yz) = temp2/npt ! "
2421     do icp=1,ncp_tot
2422     cpcloud_yz(nclouds_yz,icp)=
2423     $ cp_useds1(icp)+2*cp_useds2(icp) !store cp info
2424     enddo
2425     ptcloud_yz(nclouds_yz)=npt
2426     c ptcloud_yz_nt(nclouds_yz)=npt
2427     do ipt=1,npt
2428     db_cloud(npt_tot+ipt) = db_all(ipt)
2429     enddo
2430     npt_tot=npt_tot+npt
2431 pam-fi 1.28 if(DEBUG.EQ.1)then
2432 mocchiut 1.1 print*,'>>>> cloud ',nclouds_yz,' --- ',npt,' points'
2433 pam-fi 1.28 print*,'- alfayz1 ',alfayz1_av(nclouds_yz)
2434     print*,'- alfayz2 ',alfayz2_av(nclouds_yz)
2435     print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot)
2436     print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot)
2437     print*,'cpcloud_yz '
2438     $ ,(cpcloud_yz(nclouds_yz,icp),icp=1,ncp_tot)
2439     print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)
2440 mocchiut 1.1 endif
2441     * >>> NEW CLOUD <<<
2442     * ~~~~~~~~~~~~~~~~~
2443     2228 continue
2444     enddo !end loop (1) on DOUBLETS
2445    
2446    
2447 pam-fi 1.9 if(nloop.lt.nstepy)then
2448     cutdistyz = cutdistyz+cutystep
2449     nloop = nloop+1
2450     goto 90
2451     endif
2452    
2453 pam-fi 1.28 if(DEBUG.EQ.1)then
2454 mocchiut 1.1 print*,'Y-Z total clouds ',nclouds_yz
2455     endif
2456    
2457    
2458     return
2459     end
2460    
2461    
2462    
2463    
2464    
2465     ***************************************************
2466     * *
2467     * *
2468     * *
2469     * *
2470     * *
2471     * *
2472     **************************************************
2473    
2474     subroutine trip_to_XZcloud(iflag)
2475    
2476     include 'commontracker.f'
2477 pam-fi 1.9 include 'level1.f'
2478 mocchiut 1.1 include 'common_momanhough.f'
2479 pam-fi 1.9 c include 'momanhough_init.f'
2480 mocchiut 1.1
2481    
2482     * output flag
2483     * --------------
2484     * 0 = good event
2485     * 1 = bad event
2486     * --------------
2487     integer iflag
2488    
2489     integer tr_used(ntrpt_max)
2490     integer tr_temp(ntrpt_max)
2491     integer tr_incl(ntrpt_max)
2492     integer tr_all(ntrpt_max) !stores tr ID in each cloud
2493    
2494     integer hit_plane(nplanes)
2495    
2496     * mask for used couples
2497     integer cp_useds1(ncouplemaxtot) ! sensor 1
2498     integer cp_useds2(ncouplemaxtot) ! sensor 2
2499    
2500 pam-fi 1.28 if(DEBUG.EQ.1)print*,'trip_to_XZcloud:'
2501    
2502 mocchiut 1.1 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2503     * classification of TRIPLETS
2504     * according to distance in parameter space
2505     *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2506     do itr=1,ntrpt
2507     tr_used(itr)=0
2508     enddo
2509    
2510     distance=0
2511     nclouds_xz=0 !number of clouds
2512     npt_tot=0 !total number of selected triplets
2513 pam-fi 1.9 nloop=0
2514     91 continue
2515 mocchiut 1.1 do itr1=1,ntrpt !loop (1) on TRIPLETS
2516     if(tr_used(itr1).eq.1)goto 22288 !already included in a cloud
2517    
2518     do icp=1,ncp_tot
2519     cp_useds1(icp)=0
2520     cp_useds2(icp)=0
2521     enddo
2522     do itr=1,ntrpt
2523     tr_all(itr)=0 !list of included triplets
2524     enddo
2525     if(cpxz1(itr1).gt.0)cp_useds2(cpxz1(itr1))=1
2526     if(cpxz1(itr1).lt.0)cp_useds1(-cpxz1(itr1))=1
2527     if(cpxz2(itr1).gt.0)cp_useds2(cpxz2(itr1))=1
2528     if(cpxz2(itr1).lt.0)cp_useds1(-cpxz2(itr1))=1
2529     if(cpxz3(itr1).gt.0)cp_useds2(cpxz3(itr1))=1
2530     if(cpxz3(itr1).lt.0)cp_useds1(-cpxz3(itr1))=1
2531     temp1 = alfaxz1(itr1)
2532     temp2 = alfaxz2(itr1)
2533     temp3 = alfaxz3(itr1)
2534     npt=1 !counter of points in the cloud
2535    
2536     tr_all(npt) = itr1
2537    
2538     nptloop=1
2539     c tr_temp(1)=itr1
2540     tr_incl(1)=itr1
2541    
2542     8881 continue
2543    
2544     npv=0 !# new points inlcuded
2545     do iloop=1,nptloop
2546     itrref=tr_incl(iloop) !local point of reference
2547     do itr2=1,ntrpt !loop (2) on TRIPLETS
2548     if(itr2.eq.itr1)goto 11188 !next triplet
2549     if(tr_used(itr2).eq.1)goto 11188 !next triplet
2550 pam-fi 1.37
2551    
2552 mocchiut 1.1 * triplet distance in parameter space
2553     * solo i due parametri spaziali per il momemnto
2554     distance=
2555     $ ((alfaxz1(itrref)-alfaxz1(itr2))/Dalfaxz1)**2
2556 mocchiut 1.40 $ +((alfaxz2(itrref)-alfaxz2(itr2))/Dalfaxz2)**2
2557 mocchiut 1.1 distance = sqrt(distance)
2558 pam-fi 1.37
2559    
2560 pam-fi 1.28 * ------------------------------------------------------------------------
2561     * FORCE INCLUSION OF TRIPLETS COMPOSED BY SAME COUPLES, IGNORING THE IMAGE
2562     * ------------------------------------------------------------------------
2563     * (added in august 2007)
2564     istrimage=0
2565     if(
2566     $ abs(cpxz1(itrref)).eq.abs(cpxz1(itr2)).and.
2567     $ abs(cpxz2(itrref)).eq.abs(cpxz2(itr2)).and.
2568     $ abs(cpxz3(itrref)).eq.abs(cpxz3(itr2)).and.
2569     $ .true.)istrimage=1
2570    
2571     if(distance.lt.cutdistxz.or.istrimage.eq.1)then
2572 mocchiut 1.1 if(cpxz1(itr2).gt.0)cp_useds2(cpxz1(itr2))=1
2573     if(cpxz1(itr2).lt.0)cp_useds1(-cpxz1(itr2))=1
2574     if(cpxz2(itr2).gt.0)cp_useds2(cpxz2(itr2))=1
2575     if(cpxz2(itr2).lt.0)cp_useds1(-cpxz2(itr2))=1
2576     if(cpxz3(itr2).gt.0)cp_useds2(cpxz3(itr2))=1
2577     if(cpxz3(itr2).lt.0)cp_useds1(-cpxz3(itr2))=1
2578     npt = npt + 1 !counter of points in the cloud
2579    
2580     npv = npv +1
2581     tr_temp(npv) = itr2
2582     tr_used(itrref) = 1
2583     tr_used(itr2) = 1
2584    
2585     tr_all(npt) = itr2
2586    
2587     temp1 = temp1 + alfaxz1(itr2)
2588     temp2 = temp2 + alfaxz2(itr2)
2589     temp3 = temp3 + alfaxz3(itr2)
2590     endif
2591    
2592     11188 continue
2593     enddo !end loop (2) on TRIPLETS
2594    
2595 mocchiut 1.40 c11888 continue
2596     continue
2597 mocchiut 1.1 enddo !end loop on... bo?
2598    
2599     nptloop=npv
2600     do i=1,npv
2601     tr_incl(i)=tr_temp(i)
2602     enddo
2603     if(nptloop.ne.0)goto 8881
2604    
2605     * ------------------------------------------
2606     * stores the cloud only if
2607     * 1) it includes a minimum number of REAL couples
2608     * 1bis)
2609     * 2) it is not already stored
2610     * ------------------------------------------
2611     do ip=1,nplanes
2612     hit_plane(ip)=0
2613     enddo
2614     ncpused=0
2615     do icp=1,ncp_tot
2616 pam-fi 1.37 if(
2617     $ (cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0).and.
2618     $ .true.)then
2619 mocchiut 1.1 ncpused=ncpused+1
2620     ip=ip_cp(icp)
2621     hit_plane(ip)=1
2622     endif
2623     enddo
2624     nplused=0
2625     do ip=1,nplanes
2626     nplused=nplused+ hit_plane(ip)
2627     enddo
2628 pam-fi 1.20 if(nplused.lt.nplxz_min)goto 22288 !next triplet
2629 mocchiut 1.1
2630     * ~~~~~~~~~~~~~~~~~
2631     * >>> NEW CLOUD <<<
2632     if(nclouds_xz.ge.ncloxz_max)then
2633 pam-fi 1.28 if(verbose.eq.1)print*,
2634 mocchiut 1.1 $ '** warning ** number of identified '//
2635     $ 'XZ clouds exceeds vector dimention '
2636     $ ,'( ',ncloxz_max,' )'
2637     c good2=.false.
2638     c goto 880 !fill ntp and go to next event
2639 pam-fi 1.10 do iv=1,nviews
2640 pam-fi 1.26 c mask_view(iv) = 6
2641     mask_view(iv) = mask_view(iv) + 2**5
2642 pam-fi 1.10 enddo
2643 mocchiut 1.1 iflag=1
2644     return
2645     endif
2646     nclouds_xz = nclouds_xz + 1 !increase counter
2647     alfaxz1_av(nclouds_xz) = temp1/npt !store average parameter
2648     alfaxz2_av(nclouds_xz) = temp2/npt ! "
2649     alfaxz3_av(nclouds_xz) = temp3/npt ! "
2650     do icp=1,ncp_tot
2651     cpcloud_xz(nclouds_xz,icp)=
2652     $ cp_useds1(icp)+2*cp_useds2(icp) !store cp info
2653     enddo
2654     ptcloud_xz(nclouds_xz)=npt
2655     do ipt=1,npt
2656     tr_cloud(npt_tot+ipt) = tr_all(ipt)
2657     enddo
2658     npt_tot=npt_tot+npt
2659    
2660 pam-fi 1.28 if(DEBUG.EQ.1)then
2661 mocchiut 1.40 print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points'
2662 pam-fi 1.28 print*,'- alfaxz1 ',alfaxz1_av(nclouds_xz)
2663     print*,'- alfaxz2 ',alfaxz2_av(nclouds_xz)
2664     print*,'- alfaxz3 ',alfaxz3_av(nclouds_xz)
2665     print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot)
2666     print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot)
2667     print*,'cpcloud_xz '
2668     $ ,(cpcloud_xz(nclouds_xz,icp),icp=1,ncp_tot)
2669 mocchiut 1.1 print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)
2670     endif
2671     * >>> NEW CLOUD <<<
2672     * ~~~~~~~~~~~~~~~~~
2673     22288 continue
2674     enddo !end loop (1) on DOUBLETS
2675 pam-fi 1.9
2676     if(nloop.lt.nstepx)then
2677     cutdistxz=cutdistxz+cutxstep
2678     nloop=nloop+1
2679     goto 91
2680     endif
2681    
2682 pam-fi 1.28 if(DEBUG.EQ.1)then
2683 mocchiut 1.1 print*,'X-Z total clouds ',nclouds_xz
2684     endif
2685    
2686    
2687     return
2688     end
2689    
2690    
2691     ***************************************************
2692     * *
2693     * *
2694     * *
2695     * *
2696     * *
2697     * *
2698     **************************************************
2699    
2700     subroutine clouds_to_ctrack(iflag)
2701    
2702     include 'commontracker.f'
2703 pam-fi 1.9 include 'level1.f'
2704 mocchiut 1.1 include 'common_momanhough.f'
2705     include 'common_xyzPAM.f'
2706     include 'common_mini_2.f'
2707     include 'common_mech.f'
2708 pam-fi 1.20
2709 mocchiut 1.1
2710    
2711     * output flag
2712     * --------------
2713     * 0 = good event
2714     * 1 = bad event
2715     * --------------
2716     integer iflag
2717    
2718     * -----------------------------------------------------------
2719     * mask to store (locally) the couples included
2720     * in the intersection bewteen a XZ and YZ cloud
2721     integer cpintersec(ncouplemaxtot)
2722     * -----------------------------------------------------------
2723     * list of matching couples in the combination
2724     * between a XZ and YZ cloud
2725 pam-fi 1.9 integer cp_match(nplanes,2*ncouplemax)
2726 mocchiut 1.1 integer ncp_match(nplanes)
2727     * -----------------------------------------------------------
2728     integer hit_plane(nplanes)
2729     * -----------------------------------------------------------
2730     * variables for track fitting
2731     double precision AL_INI(5)
2732     * -----------------------------------------------------------
2733    
2734 pam-fi 1.28 if(DEBUG.EQ.1)print*,'clouds_to_ctrack:'
2735 mocchiut 1.1
2736    
2737     ntracks=0 !counter of track candidates
2738    
2739 pam-fi 1.37 do iyz=1,nclouds_yz !loop on YZ clouds
2740     do ixz=1,nclouds_xz !loop on XZ clouds
2741 mocchiut 1.1
2742     * --------------------------------------------------
2743     * check of consistency of the clouds
2744     * ---> required a minimum number of matching couples
2745     * the track fit will be performed on the INTERSECTION
2746     * of the two clouds
2747     * --------------------------------------------------
2748     do ip=1,nplanes
2749 pam-fi 1.37 hit_plane(ip)=0 !n.matching couples (REAL couples, not SINGLETS)
2750     ncp_match(ip)=0 !n.matching couples per plane
2751 mocchiut 1.1 do icpp=1,ncouplemax
2752     cp_match(ip,icpp)=0 !init couple list
2753     enddo
2754     enddo
2755 pam-fi 1.37 ncp_ok=0 !count n.matching-couples
2756     ncpx_ok=0 !count n.matching-couples with x cluster
2757     ncpy_ok=0 !count n.matching-couples with y cluster
2758    
2759    
2760 pam-fi 1.28 do icp=1,ncp_tot !loop over couples
2761 pam-fi 1.37
2762     if(.not.RECOVER_SINGLETS)then
2763     * ------------------------------------------------------
2764     * if NOT in RECOVER_SINGLETS mode, take the intersection
2765     * between xz yz clouds
2766     * ------------------------------------------------------
2767     cpintersec(icp)=min(
2768     $ cpcloud_yz(iyz,icp),
2769     $ cpcloud_xz(ixz,icp))
2770 pam-fi 1.28 * cpintersec is >0 if yz and xz clouds contain the same image of couple icp
2771 pam-fi 1.37 * ------------------------------------------------------
2772     * discard the couple if the sensor is in conflict
2773     * ------------------------------------------------------
2774     if(
2775     $ (cpcloud_yz(iyz,icp).eq.1.and.cpcloud_xz(ixz,icp).eq.2).or.
2776     $ (cpcloud_yz(iyz,icp).eq.2.and.cpcloud_xz(ixz,icp).eq.1).or.
2777     $ .false.)cpintersec(icp)=0
2778     else
2779     * ------------------------------------------------------
2780     * if RECOVER_SINGLETS take the union
2781     * (otherwise the fake couples formed by singlets would be
2782     * discarded)
2783     * ------------------------------------------------------
2784     cpintersec(icp)=max(
2785     $ cpcloud_yz(iyz,icp),
2786     $ cpcloud_xz(ixz,icp))
2787     c$$$ if(cpcloud_yz(iyz,icp).gt.0)
2788     c$$$ $ cpintersec(icp)=cpcloud_yz(iyz,icp)
2789     * cpintersec is >0 if either yz or xz cloud contains the couple icp
2790     endif
2791    
2792     c$$$ print*,icp,ip_cp(icp),' -- ',cpintersec(icp)
2793    
2794 mocchiut 1.1 if(cpintersec(icp).ne.0)then
2795    
2796     ip=ip_cp(icp)
2797     hit_plane(ip)=1
2798 pam-fi 1.37 c$$$ if(clx(ip,icp).gt.0.and.cly(ip,icp).gt.0)
2799     c$$$ $ ncp_ok=ncp_ok+1
2800     c$$$ if(clx(ip,icp).gt.0.and.cly(ip,icp).eq.0)
2801     c$$$ $ ncpx_ok=ncpx_ok+1
2802     c$$$ if(clx(ip,icp).eq.0.and.cly(ip,icp).gt.0)
2803     c$$$ $ ncpy_ok=ncpy_ok+1
2804    
2805     if( cpcloud_yz(iyz,icp).gt.0.and.
2806     $ cpcloud_xz(ixz,icp).gt.0)
2807     $ ncp_ok=ncp_ok+1
2808     if( cpcloud_yz(iyz,icp).gt.0.and.
2809     $ cpcloud_xz(ixz,icp).eq.0)
2810     $ ncpy_ok=ncpy_ok+1
2811     if( cpcloud_yz(iyz,icp).eq.0.and.
2812     $ cpcloud_xz(ixz,icp).gt.0)
2813     $ ncpx_ok=ncpx_ok+1
2814    
2815 mocchiut 1.1 if(cpintersec(icp).eq.1)then
2816     * 1) only the couple image in sensor 1 matches
2817     id=-icp
2818     ncp_match(ip)=ncp_match(ip)+1
2819     cp_match(ip,ncp_match(ip))=id
2820     elseif(cpintersec(icp).eq.2)then
2821     * 2) only the couple image in sensor 2 matches
2822     id=icp
2823     ncp_match(ip)=ncp_match(ip)+1
2824     cp_match(ip,ncp_match(ip))=id
2825     else
2826     * 3) both couple images match
2827     id=icp
2828     do is=1,2
2829     id=-id
2830     ncp_match(ip)=ncp_match(ip)+1
2831     cp_match(ip,ncp_match(ip))=id
2832     enddo
2833     endif
2834     endif !end matching condition
2835     enddo !end loop on couples
2836    
2837     nplused=0
2838     do ip=1,nplanes
2839     nplused=nplused+ hit_plane(ip)
2840     enddo
2841 pam-fi 1.37
2842     if(nplused.lt.3)goto 888 !next combination
2843     ccc if(nplused.lt.nplxz_min)goto 888 !next combination
2844     ccc if(nplused.lt.nplyz_min)goto 888 !next combination
2845     * -----------------------------------------------------------
2846     * if in RECOVER_SINGLET mode, the two clouds must have
2847     * at least ONE intersecting real couple
2848     * -----------------------------------------------------------
2849     if(ncp_ok.lt.1)goto 888 !next combination
2850    
2851 pam-fi 1.28 if(DEBUG.EQ.1)then
2852 pam-fi 1.37 print*,'////////////////////////////'
2853     print*,'Cloud combination (Y,X): ',iyz,ixz
2854     print*,' db ',ptcloud_yz(iyz)
2855     print*,' tr ',ptcloud_xz(ixz)
2856     print*,' -----> # matching couples ',ncp_ok
2857     print*,' -----> # fake couples (X)',ncpx_ok
2858     print*,' -----> # fake couples (Y)',ncpy_ok
2859     do icp=1,ncp_tot
2860     print*,'cp ',icp,' >'
2861     $ ,' x ',cpcloud_xz(ixz,icp)
2862     $ ,' y ',cpcloud_yz(iyz,icp)
2863     $ ,' ==> ',cpintersec(icp)
2864     enddo
2865     endif
2866 pam-fi 1.11
2867 pam-fi 1.28 if(DEBUG.EQ.1)then
2868 mocchiut 1.1 print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6))
2869     print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5))
2870     print*,'3 >>> ',(cp_match(4,i),i=1,ncp_match(4))
2871     print*,'4 >>> ',(cp_match(3,i),i=1,ncp_match(3))
2872     print*,'5 >>> ',(cp_match(2,i),i=1,ncp_match(2))
2873     print*,'6 >>> ',(cp_match(1,i),i=1,ncp_match(1))
2874     endif
2875    
2876     do icp1=1,max(1,ncp_match(1))
2877     hit_plane(1)=icp1
2878     if(ncp_match(1).eq.0)hit_plane(1)=0 !-icp1
2879    
2880     do icp2=1,max(1,ncp_match(2))
2881     hit_plane(2)=icp2
2882     if(ncp_match(2).eq.0)hit_plane(2)=0 !-icp2
2883    
2884     do icp3=1,max(1,ncp_match(3))
2885     hit_plane(3)=icp3
2886     if(ncp_match(3).eq.0)hit_plane(3)=0 !-icp3
2887    
2888     do icp4=1,max(1,ncp_match(4))
2889     hit_plane(4)=icp4
2890     if(ncp_match(4).eq.0)hit_plane(4)=0 !-icp4
2891    
2892     do icp5=1,max(1,ncp_match(5))
2893     hit_plane(5)=icp5
2894     if(ncp_match(5).eq.0)hit_plane(5)=0 !-icp5
2895    
2896     do icp6=1,max(1,ncp_match(6))
2897     hit_plane(6)=icp6
2898     if(ncp_match(6).eq.0)hit_plane(6)=0 !-icp6
2899    
2900 pam-fi 1.37 if(DEBUG.eq.1)
2901     $ print*,'combination: '
2902     $ ,cp_match(6,icp1)
2903     $ ,cp_match(5,icp2)
2904     $ ,cp_match(4,icp3)
2905     $ ,cp_match(3,icp4)
2906     $ ,cp_match(2,icp5)
2907     $ ,cp_match(1,icp6)
2908    
2909    
2910 pam-fi 1.28 * ---------------------------------------
2911     * check if this group of couples has been
2912     * already fitted
2913     * ---------------------------------------
2914     do ica=1,ntracks
2915     isthesame=1
2916     do ip=1,NPLANES
2917     if(hit_plane(ip).ne.0)then
2918     if( CP_STORE(nplanes-ip+1,ica)
2919     $ .ne.
2920     $ cp_match(ip,hit_plane(ip)) )
2921     $ isthesame=0
2922     else
2923     if( CP_STORE(nplanes-ip+1,ica)
2924     $ .ne.
2925     $ 0 )
2926     $ isthesame=0
2927     endif
2928     enddo
2929     if(isthesame.eq.1)then
2930     if(DEBUG.eq.1)
2931     $ print*,'(already fitted)'
2932     goto 666 !jump to next combination
2933     endif
2934     enddo
2935    
2936 mocchiut 1.1 call track_init !init TRACK common
2937    
2938 pam-fi 1.28 do ip=1,nplanes !loop on planes (bottom to top)
2939 mocchiut 1.1 if(hit_plane(ip).ne.0)then
2940     id=cp_match(ip,hit_plane(ip))
2941     is=is_cp(id)
2942     icp=icp_cp(id)
2943     if(ip_cp(id).ne.ip)
2944     $ print*,'OKKIO!!'
2945     $ ,'id ',id,is,icp
2946     $ ,ip_cp(id),ip
2947     icx=clx(ip,icp)
2948     icy=cly(ip,icp)
2949     * *************************
2950     c call xyz_PAM(icx,icy,is,
2951     c $ 'COG2','COG2',0.,0.)
2952 pam-fi 1.18 c call xyz_PAM(icx,icy,is, !(1)
2953     c $ PFAdef,PFAdef,0.,0.) !(1)
2954 mocchiut 1.1 call xyz_PAM(icx,icy,is, !(1)
2955 pam-fi 1.18 $ PFAdef,PFAdef,0.,0.,0.,0.)
2956 mocchiut 1.1 * *************************
2957     * -----------------------------
2958 pam-fi 1.37 if(icx.gt.0.and.icy.gt.0)then
2959     xgood(nplanes-ip+1)=1.
2960     ygood(nplanes-ip+1)=1.
2961     xm(nplanes-ip+1)=xPAM
2962     ym(nplanes-ip+1)=yPAM
2963     zm(nplanes-ip+1)=zPAM
2964 mocchiut 1.40 resx(nplanes-ip+1)=resxPAM
2965 pam-fi 1.37 resy(nplanes-ip+1)=resyPAM
2966     if(DEBUG.EQ.1)print*,'(X,Y)'
2967     $ ,nplanes-ip+1,xPAM,yPAM
2968     else
2969     xm_A(nplanes-ip+1) = xPAM_A
2970     ym_A(nplanes-ip+1) = yPAM_A
2971     xm_B(nplanes-ip+1) = xPAM_B
2972     ym_B(nplanes-ip+1) = yPAM_B
2973     zm(nplanes-ip+1)
2974     $ = (zPAM_A+zPAM_B)/2.
2975 mocchiut 1.40 resx(nplanes-ip+1) = resxPAM
2976 pam-fi 1.37 resy(nplanes-ip+1) = resyPAM
2977     if(icx.eq.0.and.icy.gt.0)then
2978     xgood(nplanes-ip+1)=0.
2979     ygood(nplanes-ip+1)=1.
2980     resx(nplanes-ip+1) = 1000.
2981     if(DEBUG.EQ.1)print*,'( Y)'
2982     $ ,nplanes-ip+1,xPAM,yPAM
2983     elseif(icx.gt.0.and.icy.eq.0)then
2984     xgood(nplanes-ip+1)=1.
2985     ygood(nplanes-ip+1)=0.
2986     if(DEBUG.EQ.1)print*,'(X )'
2987     $ ,nplanes-ip+1,xPAM,yPAM
2988     resy(nplanes-ip+1) = 1000.
2989     else
2990     print*,'both icx=0 and icy=0'
2991     $ ,' ==> IMPOSSIBLE!!'
2992     endif
2993     endif
2994 mocchiut 1.1 * -----------------------------
2995     endif
2996     enddo !end loop on planes
2997     * **********************************************************
2998     * ************************** FIT *** FIT *** FIT *** FIT ***
2999     * **********************************************************
3000 pam-fi 1.11 cccc scommentare se si usa al_ini della nuvola
3001     c$$$ do i=1,5
3002     c$$$ AL(i)=AL_INI(i)
3003     c$$$ enddo
3004     call guess()
3005 mocchiut 1.1 do i=1,5
3006 pam-fi 1.11 AL_INI(i)=AL(i)
3007 mocchiut 1.1 enddo
3008     ifail=0 !error flag in chi^2 computation
3009     jstep=0 !number of minimization steps
3010 pam-fi 1.8 iprint=0
3011 pam-fi 1.28 c if(DEBUG.EQ.1)iprint=1
3012     if(DEBUG.EQ.1)iprint=2
3013 pam-fi 1.8 call mini2(jstep,ifail,iprint)
3014 mocchiut 1.1 if(ifail.ne.0) then
3015 pam-fi 1.28 if(DEBUG.EQ.1)then
3016 mocchiut 1.1 print *,
3017     $ '*** MINIMIZATION FAILURE *** '
3018 pam-fi 1.11 $ //'(clouds_to_ctrack)'
3019     print*,'initial guess: '
3020    
3021     print*,'AL_INI(1) = ',AL_INI(1)
3022     print*,'AL_INI(2) = ',AL_INI(2)
3023     print*,'AL_INI(3) = ',AL_INI(3)
3024     print*,'AL_INI(4) = ',AL_INI(4)
3025     print*,'AL_INI(5) = ',AL_INI(5)
3026 mocchiut 1.1 endif
3027 pam-fi 1.11 c chi2=-chi2
3028 mocchiut 1.1 endif
3029     * **********************************************************
3030     * ************************** FIT *** FIT *** FIT *** FIT ***
3031     * **********************************************************
3032    
3033     if(chi2.le.0.)goto 666
3034 pam-fi 1.38 if(chi2.ge.1.e08)goto 666 !OPTIMIZATION
3035     if(chi2.ne.chi2)goto 666 !OPTIMIZATION
3036 mocchiut 1.1
3037     * --------------------------
3038     * STORE candidate TRACK INFO
3039     * --------------------------
3040     if(ntracks.eq.NTRACKSMAX)then
3041    
3042 pam-fi 1.28 if(verbose.eq.1)print*,
3043 mocchiut 1.1 $ '** warning ** number of candidate tracks '//
3044     $ ' exceeds vector dimension '
3045     $ ,'( ',NTRACKSMAX,' )'
3046     c good2=.false.
3047     c goto 880 !fill ntp and go to next event
3048 pam-fi 1.10 do iv=1,nviews
3049 pam-fi 1.26 c mask_view(iv) = 7
3050     mask_view(iv) = mask_view(iv) + 2**6
3051 pam-fi 1.10 enddo
3052 mocchiut 1.1 iflag=1
3053     return
3054     endif
3055    
3056     ntracks = ntracks + 1
3057    
3058 pam-fi 1.28 do ip=1,nplanes !top to bottom
3059 pam-fi 1.20
3060 mocchiut 1.1 XV_STORE(ip,ntracks)=sngl(xv(ip))
3061     YV_STORE(ip,ntracks)=sngl(yv(ip))
3062 mocchiut 1.40 ZV_STORE(ip,ntracks)=sngl(zv(ip))
3063 mocchiut 1.1 XM_STORE(ip,ntracks)=sngl(xm(ip))
3064     YM_STORE(ip,ntracks)=sngl(ym(ip))
3065     ZM_STORE(ip,ntracks)=sngl(zm(ip))
3066     RESX_STORE(ip,ntracks)=sngl(resx(ip))
3067     RESY_STORE(ip,ntracks)=sngl(resy(ip))
3068     XV_STORE(ip,ntracks)=sngl(xv(ip))
3069     YV_STORE(ip,ntracks)=sngl(yv(ip))
3070     ZV_STORE(ip,ntracks)=sngl(zv(ip))
3071     AXV_STORE(ip,ntracks)=sngl(axv(ip))
3072     AYV_STORE(ip,ntracks)=sngl(ayv(ip))
3073     XGOOD_STORE(ip,ntracks)=sngl(xgood(ip))
3074     YGOOD_STORE(ip,ntracks)=sngl(ygood(ip))
3075 pam-fi 1.28 * NB! hit_plane is defined from bottom to top
3076 mocchiut 1.1 if(hit_plane(ip).ne.0)then
3077     CP_STORE(nplanes-ip+1,ntracks)=
3078     $ cp_match(ip,hit_plane(ip))
3079 pam-fi 1.20 SENSOR_STORE(nplanes-ip+1,ntracks)
3080     $ = is_cp(cp_match(ip,hit_plane(ip)))
3081 pam-fi 1.37
3082     icl=
3083 pam-fi 1.20 $ clx(ip,icp_cp(
3084     $ cp_match(ip,hit_plane(ip)
3085 pam-fi 1.37 $ )));
3086     if(icl.eq.0)
3087     $ icl=
3088     $ cly(ip,icp_cp(
3089     $ cp_match(ip,hit_plane(ip)
3090     $ )));
3091    
3092     LADDER_STORE(nplanes-ip+1,ntracks)
3093     $ = LADDER(icl);
3094 mocchiut 1.1 else
3095     CP_STORE(nplanes-ip+1,ntracks)=0
3096 pam-fi 1.20 SENSOR_STORE(nplanes-ip+1,ntracks)=0
3097     LADDER_STORE(nplanes-ip+1,ntracks)=0
3098 mocchiut 1.1 endif
3099 pam-fi 1.28 BX_STORE(ip,ntracks)=0!I dont need it now
3100     BY_STORE(ip,ntracks)=0!I dont need it now
3101     CLS_STORE(ip,ntracks)=0
3102 mocchiut 1.1 do i=1,5
3103     AL_STORE(i,ntracks)=sngl(AL(i))
3104     enddo
3105     enddo
3106    
3107     RCHI2_STORE(ntracks)=chi2
3108    
3109     * --------------------------------
3110     * STORE candidate TRACK INFO - end
3111     * --------------------------------
3112    
3113     666 continue
3114     enddo !end loop on cp in plane 6
3115     enddo !end loop on cp in plane 5
3116     enddo !end loop on cp in plane 4
3117     enddo !end loop on cp in plane 3
3118     enddo !end loop on cp in plane 2
3119     enddo !end loop on cp in plane 1
3120    
3121     888 continue
3122     enddo !end loop on XZ couds
3123     enddo !end loop on YZ couds
3124    
3125     if(ntracks.eq.0)then
3126     iflag=1
3127 pam-fi 1.37 cc return
3128 mocchiut 1.1 endif
3129    
3130 pam-fi 1.28 if(DEBUG.EQ.1)then
3131 pam-fi 1.20 print*,'****** TRACK CANDIDATES *****************'
3132     print*,'# R. chi2 RIG ndof'
3133     do i=1,ntracks
3134     ndof=0 !(1)
3135     do ii=1,nplanes !(1)
3136     ndof=ndof !(1)
3137     $ +int(xgood_store(ii,i)) !(1)
3138     $ +int(ygood_store(ii,i)) !(1)
3139     enddo !(1)
3140     print*,i,' --- ',rchi2_store(i),' --- '
3141     $ ,1./abs(AL_STORE(5,i)),' --- ',ndof
3142     enddo
3143     print*,'*****************************************'
3144 mocchiut 1.1 endif
3145    
3146    
3147     return
3148     end
3149    
3150    
3151     ***************************************************
3152     * *
3153     * *
3154     * *
3155     * *
3156     * *
3157     * *
3158     **************************************************
3159    
3160     subroutine refine_track(ibest)
3161    
3162    
3163     include 'commontracker.f'
3164 pam-fi 1.9 include 'level1.f'
3165 mocchiut 1.1 include 'common_momanhough.f'
3166     include 'common_xyzPAM.f'
3167     include 'common_mini_2.f'
3168     include 'common_mech.f'
3169     include 'calib.f'
3170    
3171     * flag to chose PFA
3172     character*10 PFA
3173     common/FINALPFA/PFA
3174    
3175 pam-fi 1.20 real k(6)
3176     DATA k/1.099730,0.418900,0.220939,0.220907,0.418771,1.100674/
3177    
3178 pam-fi 1.18 real xp,yp,zp
3179     real xyzp(3),bxyz(3)
3180     equivalence (xp,xyzp(1)),(yp,xyzp(2)),(zp,xyzp(3))
3181    
3182 pam-fi 1.28 if(DEBUG.EQ.1)print*,'refine_track:'
3183 mocchiut 1.1 * =================================================
3184     * new estimate of positions using ETA algorithm
3185     * and
3186     * search for new couples and single clusters to add
3187     * =================================================
3188     call track_init
3189     do ip=1,nplanes !loop on planes
3190    
3191 pam-fi 1.34 if(DEBUG.EQ.1)print*,' ........... plane ',ip,' ........... '
3192    
3193 pam-fi 1.18 xP=XV_STORE(nplanes-ip+1,ibest)
3194     yP=YV_STORE(nplanes-ip+1,ibest)
3195     zP=ZV_STORE(nplanes-ip+1,ibest)
3196     call gufld(xyzp,bxyz)
3197 pam-fi 1.20 BX_STORE(nplanes-ip+1,ibest)=bxyz(1)
3198     BY_STORE(nplanes-ip+1,ibest)=bxyz(2)
3199     c$$$ bxyz(1)=0
3200 pam-fi 1.18 c$$$ bxyz(2)=0
3201     c$$$ bxyz(3)=0
3202 pam-fi 1.11 * |||||||||||||||||||||||||||||||||||||||||||||||||
3203 mocchiut 1.1 * -------------------------------------------------
3204     * If the plane has been already included, it just
3205     * computes again the coordinates of the x-y couple
3206     * using improved PFAs
3207     * -------------------------------------------------
3208 pam-fi 1.11 * |||||||||||||||||||||||||||||||||||||||||||||||||
3209 pam-fi 1.37 c$$$ if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and.
3210     c$$$ $ YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then
3211     c$$$
3212     c$$$ id=CP_STORE(nplanes-ip+1,ibest)
3213     c$$$
3214     c$$$ is=is_cp(id)
3215     c$$$ icp=icp_cp(id)
3216     c$$$ if(ip_cp(id).ne.ip)
3217     c$$$ $ print*,'OKKIO!!'
3218     c$$$ $ ,'id ',id,is,icp
3219     c$$$ $ ,ip_cp(id),ip
3220     c$$$ icx=clx(ip,icp)
3221     c$$$ icy=cly(ip,icp)
3222     c$$$c call xyz_PAM(icx,icy,is,
3223     c$$$c $ PFA,PFA,
3224     c$$$c $ AXV_STORE(nplanes-ip+1,ibest),
3225     c$$$c $ AYV_STORE(nplanes-ip+1,ibest))
3226     c$$$ call xyz_PAM(icx,icy,is,
3227     c$$$ $ PFA,PFA,
3228     c$$$ $ AXV_STORE(nplanes-ip+1,ibest),
3229     c$$$ $ AYV_STORE(nplanes-ip+1,ibest),
3230     c$$$ $ bxyz(1),
3231     c$$$ $ bxyz(2)
3232     c$$$ $ )
3233     c$$$
3234     c$$$ xm(nplanes-ip+1) = xPAM
3235     c$$$ ym(nplanes-ip+1) = yPAM
3236     c$$$ zm(nplanes-ip+1) = zPAM
3237     c$$$ xgood(nplanes-ip+1) = 1
3238     c$$$ ygood(nplanes-ip+1) = 1
3239     c$$$ resx(nplanes-ip+1) = resxPAM
3240     c$$$ resy(nplanes-ip+1) = resyPAM
3241     c$$$
3242     c$$$ dedxtrk_x(nplanes-ip+1)=sgnl(icx)/mip(VIEW(icx),LADDER(icx))
3243     c$$$ dedxtrk_y(nplanes-ip+1)=sgnl(icy)/mip(VIEW(icy),LADDER(icy))
3244     if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..or.
3245 mocchiut 1.1 $ YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then
3246    
3247     id=CP_STORE(nplanes-ip+1,ibest)
3248    
3249     is=is_cp(id)
3250     icp=icp_cp(id)
3251     if(ip_cp(id).ne.ip)
3252     $ print*,'OKKIO!!'
3253     $ ,'id ',id,is,icp
3254     $ ,ip_cp(id),ip
3255     icx=clx(ip,icp)
3256     icy=cly(ip,icp)
3257 pam-fi 1.18 c call xyz_PAM(icx,icy,is,
3258     c $ PFA,PFA,
3259     c $ AXV_STORE(nplanes-ip+1,ibest),
3260     c $ AYV_STORE(nplanes-ip+1,ibest))
3261 mocchiut 1.1 call xyz_PAM(icx,icy,is,
3262     $ PFA,PFA,
3263     $ AXV_STORE(nplanes-ip+1,ibest),
3264 pam-fi 1.18 $ AYV_STORE(nplanes-ip+1,ibest),
3265     $ bxyz(1),
3266     $ bxyz(2)
3267     $ )
3268 pam-fi 1.20
3269 pam-fi 1.37 if(icx.gt.0.and.icy.gt.0)then
3270     xm(nplanes-ip+1) = xPAM
3271     ym(nplanes-ip+1) = yPAM
3272     zm(nplanes-ip+1) = zPAM
3273     xm_A(nplanes-ip+1) = 0.
3274     ym_A(nplanes-ip+1) = 0.
3275     xm_B(nplanes-ip+1) = 0.
3276     ym_B(nplanes-ip+1) = 0.
3277     xgood(nplanes-ip+1) = 1
3278     ygood(nplanes-ip+1) = 1
3279     resx(nplanes-ip+1) = resxPAM
3280     resy(nplanes-ip+1) = resyPAM
3281     dedxtrk_x(nplanes-ip+1)=
3282     $ sgnl(icx)/mip(VIEW(icx),LADDER(icx))
3283     dedxtrk_y(nplanes-ip+1)=
3284     $ sgnl(icy)/mip(VIEW(icy),LADDER(icy))
3285     else
3286     xm(nplanes-ip+1) = 0.
3287     ym(nplanes-ip+1) = 0.
3288     zm(nplanes-ip+1) = (zPAM_A+zPAM_B)/2.
3289     xm_A(nplanes-ip+1) = xPAM_A
3290     ym_A(nplanes-ip+1) = yPAM_A
3291     xm_B(nplanes-ip+1) = xPAM_B
3292     ym_B(nplanes-ip+1) = yPAM_B
3293     xgood(nplanes-ip+1) = 0
3294     ygood(nplanes-ip+1) = 0
3295     resx(nplanes-ip+1) = 1000.!resxPAM
3296     resy(nplanes-ip+1) = 1000.!resyPAM
3297     dedxtrk_x(nplanes-ip+1)= 0
3298     dedxtrk_y(nplanes-ip+1)= 0
3299     if(icx.gt.0)then
3300     xgood(nplanes-ip+1) = 1
3301     resx(nplanes-ip+1) = resxPAM
3302     dedxtrk_x(nplanes-ip+1)=
3303     $ sgnl(icx)/mip(VIEW(icx),LADDER(icx))
3304     elseif(icy.gt.0)then
3305     ygood(nplanes-ip+1) = 1
3306     resy(nplanes-ip+1) = resyPAM
3307     dedxtrk_y(nplanes-ip+1)=
3308     $ sgnl(icy)/mip(VIEW(icy),LADDER(icy))
3309     endif
3310     endif
3311 mocchiut 1.1
3312 pam-fi 1.11 * |||||||||||||||||||||||||||||||||||||||||||||||||
3313 mocchiut 1.1 * -------------------------------------------------
3314     * If the plane has NOT been already included,
3315     * it tries to include a COUPLE or a single cluster
3316     * -------------------------------------------------
3317 pam-fi 1.11 * |||||||||||||||||||||||||||||||||||||||||||||||||
3318 mocchiut 1.1 else
3319    
3320     xgood(nplanes-ip+1)=0
3321     ygood(nplanes-ip+1)=0
3322 pam-fi 1.37
3323     CP_STORE(nplanes-ip+1,ibest)=0 !re-init
3324     CLS_STORE(nplanes-ip+1,ibest)=0
3325    
3326 mocchiut 1.1
3327     * --------------------------------------------------------------
3328     * determine which ladder and sensor are intersected by the track
3329     call whichsensor(ip,xP,yP,nldt,ist)
3330     * if the track hit the plane in a dead area, go to the next plane
3331     if(nldt.eq.0.or.ist.eq.0)goto 133
3332 pam-fi 1.20
3333     SENSOR_STORE(nplanes-ip+1,IBEST)=ist
3334     LADDER_STORE(nplanes-ip+1,IBEST)=nldt
3335 mocchiut 1.1 * --------------------------------------------------------------
3336    
3337 pam-fi 1.28 if(DEBUG.EQ.1)then
3338 mocchiut 1.1 print*,
3339     $ '------ Plane ',ip,' intersected on LADDER ',nldt
3340     $ ,' SENSOR ',ist
3341     print*,
3342     $ '------ coord: ',XP,YP
3343     endif
3344    
3345     * ===========================================
3346     * STEP 1 >>>>>>> try to include a new couple
3347     * ===========================================
3348 pam-fi 1.37 distmin=100000000.
3349 mocchiut 1.1 xmm = 0.
3350     ymm = 0.
3351     zmm = 0.
3352     rxmm = 0.
3353     rymm = 0.
3354     dedxmmx = 0. !(1)
3355     dedxmmy = 0. !(1)
3356     idm = 0 !ID of the closer couple
3357     distance=0.
3358     do icp=1,ncp_plane(ip) !loop on couples on plane icp
3359     icx=clx(ip,icp)
3360     icy=cly(ip,icp)
3361 pam-fi 1.37 if(icx.eq.0.or.icy.eq.0)goto 1188!if fake couple, jump to next
3362 mocchiut 1.1 if(LADDER(icx).ne.nldt.or. !If the ladder number does not match
3363 pam-fi 1.3 c $ cl_used(icx).eq.1.or. !or the X cluster is already used
3364     c $ cl_used(icy).eq.1.or. !or the Y cluster is already used
3365 pam-fi 1.33 $ cl_used(icx).ne.0.or. !or the X cluster is already used
3366     $ cl_used(icy).ne.0.or. !or the Y cluster is already used
3367 mocchiut 1.1 $ .false.)goto 1188 !then jump to next couple.
3368     *
3369     call xyz_PAM(icx,icy,ist,
3370     $ PFA,PFA,
3371     $ AXV_STORE(nplanes-ip+1,ibest),
3372 pam-fi 1.18 $ AYV_STORE(nplanes-ip+1,ibest),
3373     $ bxyz(1),
3374     $ bxyz(2)
3375     $ )
3376 mocchiut 1.1
3377     distance = distance_to(XP,YP)
3378 pam-fi 1.20 c distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3379 mocchiut 1.1 id=id_cp(ip,icp,ist)
3380 pam-fi 1.37 if(DEBUG.EQ.1)
3381     $ print*,'( couple ',id
3382 pam-fi 1.20 $ ,' ) distance ',distance
3383 mocchiut 1.1 if(distance.lt.distmin)then
3384     xmm = xPAM
3385     ymm = yPAM
3386     zmm = zPAM
3387     rxmm = resxPAM
3388     rymm = resyPAM
3389     distmin = distance
3390     idm = id
3391 pam-fi 1.18 dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
3392     dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
3393 pam-fi 1.37 clincnewc=10*sqrt(rymm**2+rxmm**2
3394     $ +RCHI2_STORE(ibest)*k(ip)*(cov(1,1)+cov(2,2)))
3395 mocchiut 1.1 endif
3396     1188 continue
3397     enddo !end loop on couples on plane icp
3398 pam-fi 1.37 if(distmin.le.clincnewc)then
3399 mocchiut 1.1 * -----------------------------------
3400 pam-fi 1.20 xm(nplanes-ip+1) = xmm !<<<
3401     ym(nplanes-ip+1) = ymm !<<<
3402     zm(nplanes-ip+1) = zmm !<<<
3403     xgood(nplanes-ip+1) = 1 !<<<
3404     ygood(nplanes-ip+1) = 1 !<<<
3405     resx(nplanes-ip+1)=rxmm !<<<
3406     resy(nplanes-ip+1)=rymm !<<<
3407     dedxtrk_x(nplanes-ip+1) = dedxmmx !<<<
3408     dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<
3409 mocchiut 1.1 * -----------------------------------
3410     CP_STORE(nplanes-ip+1,ibest)=idm
3411 pam-fi 1.28 if(DEBUG.EQ.1)print*,'%%%% included couple ',idm
3412 pam-fi 1.37 $ ,' (dist.= ',distmin,', cut ',clincnewc,' )'
3413 mocchiut 1.1 goto 133 !next plane
3414     endif
3415     * ================================================
3416     * STEP 2 >>>>>>> try to include a single cluster
3417     * either from a couple or single
3418     * ================================================
3419     distmin=1000000.
3420     xmm_A = 0. !---------------------------
3421     ymm_A = 0. ! init variables that
3422     zmm_A = 0. ! define the SINGLET
3423     xmm_B = 0. !
3424     ymm_B = 0. !
3425     zmm_B = 0. !
3426     rxmm = 0. !
3427     rymm = 0. !
3428     dedxmmx = 0. !(1)
3429     dedxmmy = 0. !(1)
3430     iclm=0 !---------------------------
3431     distance=0.
3432    
3433     *----- clusters inside couples -------------------------------------
3434     do icp=1,ncp_plane(ip) !loop on cluster inside couples
3435     icx=clx(ip,icp)
3436     icy=cly(ip,icp)
3437 pam-fi 1.37 if(icx.eq.0.or.icy.eq.0)goto 11882!if fake couple, jump to next
3438 mocchiut 1.1 id=id_cp(ip,icp,ist)
3439     if(LADDER(icx).ne.nldt)goto 11882 !if the ladder number does not match
3440     * !jump to the next couple
3441     *----- try cluster x -----------------------------------------------
3442 pam-fi 1.3 c if(cl_used(icx).eq.1)goto 11881 !if the X cluster is already used
3443     if(cl_used(icx).ne.0)goto 11881 !if the X cluster is already used !(3)
3444 mocchiut 1.1 * !jump to the Y cluster
3445 pam-fi 1.18 c call xyz_PAM(icx,0,ist,
3446     c $ PFA,PFA,
3447     c $ AXV_STORE(nplanes-ip+1,ibest),0.)
3448 mocchiut 1.1 call xyz_PAM(icx,0,ist,
3449     $ PFA,PFA,
3450 pam-fi 1.18 $ AXV_STORE(nplanes-ip+1,ibest),0.,
3451     $ bxyz(1),
3452     $ bxyz(2)
3453     $ )
3454 mocchiut 1.1 distance = distance_to(XP,YP)
3455 pam-fi 1.20 c distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3456 pam-fi 1.37 if(DEBUG.EQ.1)
3457     $ print*,'( cl-X ',icx
3458 pam-fi 1.20 $ ,' in cp ',id,' ) distance ',distance
3459 mocchiut 1.1 if(distance.lt.distmin)then
3460     xmm_A = xPAM_A
3461     ymm_A = yPAM_A
3462     zmm_A = zPAM_A
3463     xmm_B = xPAM_B
3464     ymm_B = yPAM_B
3465     zmm_B = zPAM_B
3466     rxmm = resxPAM
3467     rymm = resyPAM
3468     distmin = distance
3469     iclm = icx
3470 pam-fi 1.18 c dedxmm = sgnl(icx) !(1)
3471     dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
3472 mocchiut 1.1 dedxmmy = 0. !(1)
3473     endif
3474     11881 continue
3475     *----- try cluster y -----------------------------------------------
3476 pam-fi 1.3 c if(cl_used(icy).eq.1)goto 11882 !if the Y cluster is already used
3477     if(cl_used(icy).ne.0)goto 11882 !if the Y cluster is already used !(3)
3478 mocchiut 1.1 * !jump to the next couple
3479 pam-fi 1.18 c call xyz_PAM(0,icy,ist,
3480     c $ PFA,PFA,
3481     c $ 0.,AYV_STORE(nplanes-ip+1,ibest))
3482 mocchiut 1.1 call xyz_PAM(0,icy,ist,
3483     $ PFA,PFA,
3484 pam-fi 1.18 $ 0.,AYV_STORE(nplanes-ip+1,ibest),
3485     $ bxyz(1),
3486     $ bxyz(2)
3487     $ )
3488 mocchiut 1.1 distance = distance_to(XP,YP)
3489 pam-fi 1.20 c distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3490 pam-fi 1.37 if(DEBUG.EQ.1)
3491     $ print*,'( cl-Y ',icy
3492 pam-fi 1.20 $ ,' in cp ',id,' ) distance ',distance
3493 mocchiut 1.1 if(distance.lt.distmin)then
3494     xmm_A = xPAM_A
3495     ymm_A = yPAM_A
3496     zmm_A = zPAM_A
3497     xmm_B = xPAM_B
3498     ymm_B = yPAM_B
3499     zmm_B = zPAM_B
3500     rxmm = resxPAM
3501     rymm = resyPAM
3502     distmin = distance
3503     iclm = icy
3504 pam-fi 1.18 c dedxmm = sgnl(icy) !(1)
3505 mocchiut 1.1 dedxmmx = 0. !(1)
3506 pam-fi 1.18 dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
3507 mocchiut 1.1 endif
3508     11882 continue
3509     enddo !end loop on cluster inside couples
3510 pam-fi 1.17 *----- single clusters -----------------------------------------------
3511 mocchiut 1.1 do ic=1,ncls(ip) !loop on single clusters
3512     icl=cls(ip,ic)
3513 pam-fi 1.3 if(cl_used(icl).ne.0.or. !if the cluster is already used !(3)
3514 mocchiut 1.1 $ LADDER(icl).ne.nldt.or. !or the ladder number does not match
3515     $ .false.)goto 18882 !jump to the next singlet
3516     if(mod(VIEW(icl),2).eq.0)then!<---- X view
3517     call xyz_PAM(icl,0,ist,
3518     $ PFA,PFA,
3519 pam-fi 1.18 $ AXV_STORE(nplanes-ip+1,ibest),0.,
3520     $ bxyz(1),
3521     $ bxyz(2)
3522     $ )
3523 mocchiut 1.1 else !<---- Y view
3524     call xyz_PAM(0,icl,ist,
3525     $ PFA,PFA,
3526 pam-fi 1.18 $ 0.,AYV_STORE(nplanes-ip+1,ibest),
3527     $ bxyz(1),
3528     $ bxyz(2)
3529     $ )
3530 mocchiut 1.1 endif
3531    
3532     distance = distance_to(XP,YP)
3533 pam-fi 1.20 c distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
3534 pam-fi 1.37 if(DEBUG.EQ.1)
3535     $ print*,'( cl-s ',icl
3536 pam-fi 1.28 $ ,' ) distance ',distance
3537 mocchiut 1.1 if(distance.lt.distmin)then
3538     xmm_A = xPAM_A
3539     ymm_A = yPAM_A
3540     zmm_A = zPAM_A
3541     xmm_B = xPAM_B
3542     ymm_B = yPAM_B
3543     zmm_B = zPAM_B
3544     rxmm = resxPAM
3545     rymm = resyPAM
3546     distmin = distance
3547     iclm = icl
3548     if(mod(VIEW(icl),2).eq.0)then !<---- X view
3549 pam-fi 1.20 dedxmmx = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3550     dedxmmy = 0.
3551 mocchiut 1.1 else !<---- Y view
3552 pam-fi 1.20 dedxmmx = 0.
3553     dedxmmy = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
3554 mocchiut 1.1 endif
3555     endif
3556     18882 continue
3557     enddo !end loop on single clusters
3558 pam-fi 1.20
3559     if(iclm.ne.0)then
3560     if(mod(VIEW(iclm),2).eq.0)then
3561     clincnew=
3562     $ 20*
3563     $ sqrt(rxmm**2+RCHI2_STORE(ibest)*k(ip)*cov(1,1))
3564     else if(mod(VIEW(iclm),2).ne.0)then
3565     clincnew=
3566     $ 10*
3567     $ sqrt(rymm**2+RCHI2_STORE(ibest)*k(ip)*cov(2,2))
3568     endif
3569 pam-fi 1.37
3570     if(distmin.le.clincnew)then
3571 pam-fi 1.20
3572     CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<
3573     * ----------------------------
3574     if(mod(VIEW(iclm),2).eq.0)then
3575     XGOOD(nplanes-ip+1)=1.
3576     resx(nplanes-ip+1)=rxmm
3577 pam-fi 1.37 if(DEBUG.EQ.1)
3578     $ print*,'%%%% included X-cl ',iclm
3579 pam-fi 1.20 $ ,'( chi^2, ',RCHI2_STORE(ibest)
3580     $ ,', dist.= ',distmin
3581 pam-fi 1.37 $ ,', cut ',clincnew,' )'
3582 pam-fi 1.20 else
3583     YGOOD(nplanes-ip+1)=1.
3584     resy(nplanes-ip+1)=rymm
3585 pam-fi 1.37 if(DEBUG.EQ.1)
3586     $ print*,'%%%% included Y-cl ',iclm
3587 pam-fi 1.20 $ ,'( chi^2, ',RCHI2_STORE(ibest)
3588     $ ,', dist.= ', distmin
3589 pam-fi 1.37 $ ,', cut ',clincnew,' )'
3590 pam-fi 1.20 endif
3591     * ----------------------------
3592     xm_A(nplanes-ip+1) = xmm_A
3593     ym_A(nplanes-ip+1) = ymm_A
3594     xm_B(nplanes-ip+1) = xmm_B
3595     ym_B(nplanes-ip+1) = ymm_B
3596     zm(nplanes-ip+1) = (zmm_A+zmm_B)/2.
3597     dedxtrk_x(nplanes-ip+1) = dedxmmx !<<<
3598     dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<
3599     * ----------------------------
3600 mocchiut 1.1 endif
3601     endif
3602     endif
3603     133 continue
3604     enddo !end loop on planes
3605    
3606    
3607    
3608     return
3609     end
3610    
3611 pam-fi 1.37
3612 mocchiut 1.1 ***************************************************
3613     * *
3614     * *
3615     * *
3616     * *
3617     * *
3618     * *
3619     **************************************************
3620 pam-fi 1.3 *
3621 mocchiut 1.1
3622    
3623    
3624     * ****************************************************
3625    
3626     subroutine init_level2
3627    
3628     include 'commontracker.f'
3629 pam-fi 1.9 include 'level1.f'
3630 mocchiut 1.1 include 'common_momanhough.f'
3631     include 'level2.f'
3632    
3633 pam-fi 1.20 * ---------------------------------
3634     * variables initialized from level1
3635     * ---------------------------------
3636 pam-fi 1.4 do i=1,nviews
3637     good2(i)=good1(i)
3638 pam-fi 1.20 do j=1,nva1_view
3639     vkflag(i,j)=1
3640     if(cnnev(i,j).le.0)then
3641     vkflag(i,j)=cnnev(i,j)
3642     endif
3643     enddo
3644 pam-fi 1.4 enddo
3645 pam-fi 1.20 * ----------------
3646     * level2 variables
3647     * ----------------
3648 mocchiut 1.1 NTRK = 0
3649 pam-fi 1.9 do it=1,NTRKMAX
3650 mocchiut 1.1 IMAGE(IT)=0
3651     CHI2_nt(IT) = -100000.
3652     do ip=1,nplanes
3653     XM_nt(IP,IT) = 0
3654     YM_nt(IP,IT) = 0
3655     ZM_nt(IP,IT) = 0
3656     RESX_nt(IP,IT) = 0
3657     RESY_nt(IP,IT) = 0
3658 pam-fi 1.20 TAILX_nt(IP,IT) = 0
3659     TAILY_nt(IP,IT) = 0
3660     XBAD(IP,IT) = 0
3661     YBAD(IP,IT) = 0
3662 mocchiut 1.1 XGOOD_nt(IP,IT) = 0
3663     YGOOD_nt(IP,IT) = 0
3664 pam-fi 1.20 LS(IP,IT) = 0
3665 mocchiut 1.1 DEDX_X(IP,IT) = 0
3666     DEDX_Y(IP,IT) = 0
3667 pam-fi 1.3 CLTRX(IP,IT) = 0
3668     CLTRY(IP,IT) = 0
3669 pam-fi 1.31 multmaxx(ip,it) = 0
3670     seedx(ip,it) = 0
3671     xpu(ip,it) = 0
3672     multmaxy(ip,it) = 0
3673     seedy(ip,it) = 0
3674     ypu(ip,it) = 0
3675 mocchiut 1.1 enddo
3676     do ipa=1,5
3677     AL_nt(IPA,IT) = 0
3678     do ipaa=1,5
3679     coval(ipa,ipaa,IT)=0
3680     enddo
3681     enddo
3682     enddo
3683     nclsx=0
3684     nclsy=0
3685     do ip=1,NSINGMAX
3686     planex(ip)=0
3687     xs(1,ip)=0
3688     xs(2,ip)=0
3689     sgnlxs(ip)=0
3690     planey(ip)=0
3691     ys(1,ip)=0
3692     ys(2,ip)=0
3693     sgnlys(ip)=0
3694 pam-fi 1.34 sxbad(ip)=0
3695     sybad(ip)=0
3696     multmaxsx(ip)=0
3697     multmaxsy(ip)=0
3698 mocchiut 1.1 enddo
3699     end
3700    
3701    
3702     ************************************************************
3703     *
3704     *
3705     *
3706     *
3707     *
3708     *
3709     *
3710     ************************************************************
3711    
3712    
3713 pam-fi 1.9 subroutine init_hough
3714    
3715     include 'commontracker.f'
3716     include 'level1.f'
3717     include 'common_momanhough.f'
3718     include 'common_hough.f'
3719     include 'level2.f'
3720    
3721     ntrpt_nt=0
3722     ndblt_nt=0
3723     NCLOUDS_XZ_nt=0
3724     NCLOUDS_YZ_nt=0
3725     do idb=1,ndblt_max_nt
3726     db_cloud_nt(idb)=0
3727     alfayz1_nt(idb)=0
3728     alfayz2_nt(idb)=0
3729     enddo
3730 pam-fi 1.14 do itr=1,ntrpt_max_nt
3731 pam-fi 1.9 tr_cloud_nt(itr)=0
3732     alfaxz1_nt(itr)=0
3733     alfaxz2_nt(itr)=0
3734     alfaxz3_nt(itr)=0
3735     enddo
3736     do idb=1,ncloyz_max
3737     ptcloud_yz_nt(idb)=0
3738     alfayz1_av_nt(idb)=0
3739     alfayz2_av_nt(idb)=0
3740     enddo
3741     do itr=1,ncloxz_max
3742     ptcloud_xz_nt(itr)=0
3743     alfaxz1_av_nt(itr)=0
3744     alfaxz2_av_nt(itr)=0
3745     alfaxz3_av_nt(itr)=0
3746     enddo
3747    
3748     ntrpt=0
3749     ndblt=0
3750     NCLOUDS_XZ=0
3751     NCLOUDS_YZ=0
3752     do idb=1,ndblt_max
3753     db_cloud(idb)=0
3754     cpyz1(idb)=0
3755     cpyz2(idb)=0
3756     alfayz1(idb)=0
3757     alfayz2(idb)=0
3758     enddo
3759 pam-fi 1.14 do itr=1,ntrpt_max
3760 pam-fi 1.9 tr_cloud(itr)=0
3761     cpxz1(itr)=0
3762     cpxz2(itr)=0
3763     cpxz3(itr)=0
3764     alfaxz1(itr)=0
3765     alfaxz2(itr)=0
3766     alfaxz3(itr)=0
3767     enddo
3768     do idb=1,ncloyz_max
3769     ptcloud_yz(idb)=0
3770     alfayz1_av(idb)=0
3771     alfayz2_av(idb)=0
3772     do idbb=1,ncouplemaxtot
3773     cpcloud_yz(idb,idbb)=0
3774     enddo
3775     enddo
3776     do itr=1,ncloxz_max
3777     ptcloud_xz(itr)=0
3778     alfaxz1_av(itr)=0
3779     alfaxz2_av(itr)=0
3780     alfaxz3_av(itr)=0
3781     do itrr=1,ncouplemaxtot
3782     cpcloud_xz(itr,itrr)=0
3783     enddo
3784     enddo
3785     end
3786     ************************************************************
3787     *
3788     *
3789     *
3790     *
3791     *
3792     *
3793     *
3794     ************************************************************
3795    
3796    
3797 mocchiut 1.1 subroutine fill_level2_tracks(ntr)
3798    
3799     * -------------------------------------------------------
3800     * This routine fills the ntr-th element of the variables
3801     * inside the level2_tracks common, which correspond
3802     * to the ntr-th track info.
3803     * -------------------------------------------------------
3804    
3805    
3806     include 'commontracker.f'
3807 pam-fi 1.3 include 'level1.f'
3808 pam-fi 1.9 include 'common_momanhough.f'
3809 mocchiut 1.1 include 'level2.f'
3810     include 'common_mini_2.f'
3811 pam-fi 1.20 include 'calib.f'
3812    
3813     character*10 PFA
3814     common/FINALPFA/PFA
3815    
3816     real sinth,phi,pig
3817     integer ssensor,sladder
3818 mocchiut 1.1 pig=acos(-1.)
3819    
3820 pam-fi 1.37
3821    
3822 pam-fi 1.20 * -------------------------------------
3823 mocchiut 1.1 chi2_nt(ntr) = sngl(chi2)
3824 pam-fi 1.8 nstep_nt(ntr) = nstep
3825 pam-fi 1.20 * -------------------------------------
3826 pam-fi 1.9 phi = al(4)
3827     sinth = al(3)
3828     if(sinth.lt.0)then
3829     sinth = -sinth
3830     phi = phi + pig
3831     endif
3832     npig = aint(phi/(2*pig))
3833     phi = phi - npig*2*pig
3834     if(phi.lt.0)
3835     $ phi = phi + 2*pig
3836     al(4) = phi
3837     al(3) = sinth
3838 mocchiut 1.1 do i=1,5
3839     al_nt(i,ntr) = sngl(al(i))
3840     do j=1,5
3841     coval(i,j,ntr) = sngl(cov(i,j))
3842     enddo
3843     enddo
3844 pam-fi 1.20 * -------------------------------------
3845 mocchiut 1.1 do ip=1,nplanes ! loop on planes
3846     xgood_nt(ip,ntr) = int(xgood(ip))
3847     ygood_nt(ip,ntr) = int(ygood(ip))
3848     xm_nt(ip,ntr) = sngl(xm(ip))
3849     ym_nt(ip,ntr) = sngl(ym(ip))
3850     zm_nt(ip,ntr) = sngl(zm(ip))
3851     RESX_nt(IP,ntr) = sngl(resx(ip))
3852     RESY_nt(IP,ntr) = sngl(resy(ip))
3853 pam-fi 1.20 TAILX_nt(IP,ntr) = 0.
3854     TAILY_nt(IP,ntr) = 0.
3855 mocchiut 1.1 xv_nt(ip,ntr) = sngl(xv(ip))
3856     yv_nt(ip,ntr) = sngl(yv(ip))
3857     zv_nt(ip,ntr) = sngl(zv(ip))
3858     axv_nt(ip,ntr) = sngl(axv(ip))
3859 pam-fi 1.18 ayv_nt(ip,ntr) = sngl(ayv(ip))
3860 pam-fi 1.28
3861 pam-fi 1.18 factor = sqrt(
3862 pam-fi 1.19 $ tan( acos(-1.) * sngl(axv(ip)) /180. )**2 +
3863     $ tan( acos(-1.) * sngl(ayv(ip)) /180. )**2 +
3864 pam-fi 1.18 $ 1. )
3865 pam-fi 1.28
3866 pam-fi 1.18 dedx_x(ip,ntr) = sngl(dedxtrk_x(ip)/factor)
3867     dedx_y(ip,ntr) = sngl(dedxtrk_y(ip)/factor)
3868 pam-fi 1.3
3869 pam-fi 1.37
3870     ccc print*,ip,'dedx >>> ',dedx_x(ip,ntr),dedx_y(ip,ntr)
3871    
3872 pam-fi 1.20 ax = axv_nt(ip,ntr)
3873     ay = ayv_nt(ip,ntr)
3874     bfx = BX_STORE(ip,IDCAND)
3875     bfy = BY_STORE(ip,IDCAND)
3876 pam-fi 1.31 c$$$ if(ip.eq.6) ax = -1. * axv_nt(ip,ntr)
3877     c$$$ if(ip.eq.6) bfy = -1. * BY_STORE(ip,IDCAND)
3878     c$$$ tgtemp = tan(ax*acos(-1.)/180.) + pmuH_h*bfy*0.00001
3879     c$$$ angx = 180.*atan(tgtemp)/acos(-1.)
3880     c$$$ tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001
3881     c$$$ angy = 180.*atan(tgtemp)/acos(-1.)
3882    
3883     angx = effectiveangle(ax,2*ip,bfy)
3884     angy = effectiveangle(ay,2*ip-1,bfx)
3885    
3886 pam-fi 1.20
3887    
3888     id = CP_STORE(ip,IDCAND) ! couple id
3889 pam-fi 1.3 icl = CLS_STORE(ip,IDCAND)
3890 pam-fi 1.20 ssensor = -1
3891     sladder = -1
3892     ssensor = SENSOR_STORE(ip,IDCAND)
3893     sladder = LADDER_STORE(ip,IDCAND)
3894     if(ip.eq.6.and.ssensor.ne.0)ssensor = 3 - ssensor !notazione paolo x align
3895     LS(IP,ntr) = ssensor+10*sladder
3896    
3897 pam-fi 1.42 c if(id.ne.0)then
3898     CCCCCC(10 novembre 2009) PATCH X NUCLEI
3899     C non ho capito perche', ma durante il ritracciamento dei nuclei
3900     C (quando una traccia viene trovata ma non e' buona) c'e' qualche variabile
3901     C che non viene reinizializzata correttamente e i cluster esclusi
3902     C dal fit risultano ancora inclusi...
3903     C
3904     cltrx(ip,ntr) = 0
3905     cltry(ip,ntr) = 0
3906     if(
3907     $ xgood_nt(ip,ntr).eq.1.and.ygood_nt(ip,ntr).eq.1
3908     $ .and.
3909     $ id.ne.0)then
3910    
3911 pam-fi 1.20 c >>> is a couple
3912 pam-fi 1.3 cltrx(ip,ntr) = clx(nplanes-ip+1,icp_cp(id))
3913     cltry(ip,ntr) = cly(nplanes-ip+1,icp_cp(id))
3914 pam-fi 1.28
3915 pam-fi 1.37 if(clx(nplanes-ip+1,icp_cp(id)).gt.0)then
3916 pam-fi 1.28
3917 pam-fi 1.37 cl_used(cltrx(ip,ntr)) = 1 !tag used clusters
3918 pam-fi 1.20
3919 mocchiut 1.40 xbad(ip,ntr)= nbadstrips(4,clx(nplanes-ip+1,icp_cp(id)))
3920 pam-fi 1.37
3921     if(nsatstrips(clx(nplanes-ip+1,icp_cp(id))).gt.0)
3922     $ dedx_x(ip,ntr)=-dedx_x(ip,ntr)
3923    
3924     multmaxx(ip,ntr) = maxs(cltrx(ip,ntr))
3925     $ +10000*mult(cltrx(ip,ntr))
3926     seedx(ip,ntr) = clsignal(indmax(cltrx(ip,ntr)))
3927     $ /clsigma(indmax(cltrx(ip,ntr)))
3928     call applypfa(PFA,cltrx(ip,ntr),angx,corr,res)
3929     xpu(ip,ntr) = corr
3930    
3931     endif
3932     if(cly(nplanes-ip+1,icp_cp(id)).gt.0)then
3933 pam-fi 1.20
3934 pam-fi 1.37 cl_used(cltry(ip,ntr)) = 1 !tag used clusters
3935 pam-fi 1.20
3936 pam-fi 1.37 ybad(ip,ntr)= nbadstrips(4,cly(nplanes-ip+1,icp_cp(id)))
3937 pam-fi 1.31
3938 pam-fi 1.37 if(nsatstrips(cly(nplanes-ip+1,icp_cp(id))).gt.0)
3939     $ dedx_y(ip,ntr)=-dedx_y(ip,ntr)
3940    
3941     multmaxy(ip,ntr) = maxs(cltry(ip,ntr))
3942     $ +10000*mult(cltry(ip,ntr))
3943     seedy(ip,ntr) = clsignal(indmax(cltry(ip,ntr)))
3944     $ /clsigma(indmax(cltry(ip,ntr)))
3945     call applypfa(PFA,cltry(ip,ntr),angy,corr,res)
3946     ypu(ip,ntr) = corr
3947     endif
3948 pam-fi 1.31
3949 pam-fi 1.3 elseif(icl.ne.0)then
3950 pam-fi 1.26
3951 pam-fi 1.28 cl_used(icl) = 1 !tag used clusters
3952    
3953 pam-fi 1.20 if(mod(VIEW(icl),2).eq.0)then
3954     cltrx(ip,ntr)=icl
3955 pam-fi 1.26 xbad(ip,ntr) = nbadstrips(4,icl)
3956    
3957 pam-fi 1.20 if(nsatstrips(icl).gt.0)dedx_x(ip,ntr)=-dedx_x(ip,ntr)
3958 pam-fi 1.31
3959     multmaxx(ip,ntr) = maxs(cltrx(ip,ntr))
3960     $ +10000*mult(cltrx(ip,ntr))
3961     seedx(ip,ntr) = clsignal(indmax(cltrx(ip,ntr)))
3962     $ /clsigma(indmax(cltrx(ip,ntr)))
3963     call applypfa(PFA,cltrx(ip,ntr),angx,corr,res)
3964     xpu(ip,ntr) = corr
3965    
3966 pam-fi 1.20 elseif(mod(VIEW(icl),2).eq.1)then
3967     cltry(ip,ntr)=icl
3968 pam-fi 1.26 ybad(ip,ntr) = nbadstrips(4,icl)
3969 pam-fi 1.31
3970 pam-fi 1.20 if(nsatstrips(icl).gt.0)dedx_y(ip,ntr)=-dedx_y(ip,ntr)
3971 pam-fi 1.31
3972     multmaxy(ip,ntr) = maxs(cltry(ip,ntr))
3973     $ +10000*mult(cltry(ip,ntr))
3974     seedy(ip,ntr) = clsignal(indmax(cltry(ip,ntr)))
3975     $ /clsigma(indmax(cltry(ip,ntr)))
3976     call applypfa(PFA,cltry(ip,ntr),angy,corr,res)
3977     ypu(ip,ntr) = corr
3978    
3979 pam-fi 1.20 endif
3980 pam-fi 1.26
3981 pam-fi 1.3 endif
3982    
3983 mocchiut 1.1 enddo
3984    
3985 pam-fi 1.28 if(DEBUG.eq.1)then
3986     print*,'> STORING TRACK ',ntr
3987     print*,'clusters: '
3988     do ip=1,6
3989     print*,'> ',ip,' -- ',cltrx(ip,ntr),cltry(ip,ntr)
3990     enddo
3991 pam-fi 1.37 print*,'dedx: '
3992     do ip=1,6
3993     print*,'> ',ip,' -- ',dedx_x(ip,ntr),dedx_y(ip,ntr)
3994     enddo
3995 pam-fi 1.28 endif
3996 mocchiut 1.1
3997     end
3998    
3999     subroutine fill_level2_siglets
4000    
4001     * -------------------------------------------------------
4002     * This routine fills the elements of the variables
4003     * inside the level2_singletsx and level2_singletsy commons,
4004     * which store info on clusters outside the tracks
4005     * -------------------------------------------------------
4006    
4007     include 'commontracker.f'
4008 pam-fi 1.9 include 'calib.f'
4009 mocchiut 1.1 include 'level1.f'
4010 pam-fi 1.9 include 'common_momanhough.f'
4011 mocchiut 1.1 include 'level2.f'
4012     include 'common_xyzPAM.f'
4013    
4014     * count #cluster per plane not associated to any track
4015     nclsx = 0
4016     nclsy = 0
4017    
4018 pam-fi 1.9 do iv = 1,nviews
4019 pam-fi 1.26 c if( mask_view(iv).ne.0 )good2(iv) = 20+mask_view(iv)
4020     good2(iv) = good2(iv) + mask_view(iv)*2**8
4021 pam-fi 1.9 enddo
4022    
4023 pam-fi 1.28 if(DEBUG.eq.1)then
4024     print*,'> STORING SINGLETS '
4025     endif
4026    
4027 mocchiut 1.1 do icl=1,nclstr1
4028 pam-fi 1.28
4029     ip=nplanes-npl(VIEW(icl))+1
4030    
4031 mocchiut 1.1 if(cl_used(icl).eq.0)then !cluster not included in any track
4032 pam-fi 1.34
4033 mocchiut 1.1 if(mod(VIEW(icl),2).eq.0)then !=== X views
4034 pam-fi 1.34
4035 mocchiut 1.1 nclsx = nclsx + 1
4036     planex(nclsx) = ip
4037 pam-fi 1.20 sgnlxs(nclsx) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
4038     if(nsatstrips(icl).gt.0)sgnlxs(nclsx)=-sgnlxs(nclsx)
4039 pam-fi 1.3 clsx(nclsx) = icl
4040 pam-fi 1.34 sxbad(nclsx) = nbadstrips(1,icl)
4041     multmaxsx(nclsx) = maxs(icl)+10000*mult(icl)
4042    
4043    
4044 mocchiut 1.1 do is=1,2
4045     c call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)
4046 pam-fi 1.18 c call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)
4047     call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.,0.,0.)
4048 mocchiut 1.1 xs(is,nclsx) = (xPAM_A+xPAM_B)/2
4049     enddo
4050     else !=== Y views
4051     nclsy = nclsy + 1
4052     planey(nclsy) = ip
4053 pam-fi 1.20 sgnlys(nclsy) = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
4054     if(nsatstrips(icl).gt.0)sgnlys(nclsy)=-sgnlys(nclsy)
4055 pam-fi 1.3 clsy(nclsy) = icl
4056 pam-fi 1.34 sybad(nclsy) = nbadstrips(1,icl)
4057     multmaxsy(nclsy) = maxs(icl)+10000*mult(icl)
4058    
4059    
4060 mocchiut 1.1 do is=1,2
4061     c call xyz_PAM(0,icl,is,' ','COG1',0.,0.)
4062 pam-fi 1.18 c call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)
4063     call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.,0.,0.)
4064 mocchiut 1.1 ys(is,nclsy) = (yPAM_A+yPAM_B)/2
4065     enddo
4066     endif
4067     endif
4068 pam-fi 1.3
4069     ***** LO METTO QUI PERCHE` NON SO DOVE METTERLO
4070     whichtrack(icl) = cl_used(icl)
4071 pam-fi 1.28 * --------------------------------------------------
4072     * per non perdere la testa...
4073     * whichtrack(icl) e` una variabile del common level1
4074     * che serve solo per sapere quali cluster sono stati
4075     * associati ad una traccia, e permettere di salvare
4076     * solo questi nell'albero di uscita
4077     * --------------------------------------------------
4078 pam-fi 1.33
4079 mocchiut 1.1 enddo
4080     end
4081    
4082 pam-fi 1.9 ***************************************************
4083     * *
4084     * *
4085     * *
4086     * *
4087     * *
4088     * *
4089     **************************************************
4090 mocchiut 1.1
4091 pam-fi 1.9 subroutine fill_hough
4092 mocchiut 1.1
4093 pam-fi 1.9 * -------------------------------------------------------
4094     * This routine fills the variables related to the hough
4095     * transform, for the debig n-tuple
4096     * -------------------------------------------------------
4097 mocchiut 1.1
4098 pam-fi 1.9 include 'commontracker.f'
4099     include 'level1.f'
4100     include 'common_momanhough.f'
4101     include 'common_hough.f'
4102     include 'level2.f'
4103 mocchiut 1.1
4104 pam-fi 1.9 if(.false.
4105     $ .or.ntrpt.gt.ntrpt_max_nt
4106     $ .or.ndblt.gt.ndblt_max_nt
4107     $ .or.NCLOUDS_XZ.gt.ncloxz_max
4108     $ .or.NCLOUDS_yZ.gt.ncloyz_max
4109     $ )then
4110     ntrpt_nt=0
4111     ndblt_nt=0
4112     NCLOUDS_XZ_nt=0
4113     NCLOUDS_YZ_nt=0
4114     else
4115     ndblt_nt=ndblt
4116     ntrpt_nt=ntrpt
4117     if(ndblt.ne.0)then
4118     do id=1,ndblt
4119     alfayz1_nt(id)=alfayz1(id) !Y0
4120     alfayz2_nt(id)=alfayz2(id) !tg theta-yz
4121     enddo
4122     endif
4123     if(ndblt.ne.0)then
4124     do it=1,ntrpt
4125     alfaxz1_nt(it)=alfaxz1(it) !X0
4126     alfaxz2_nt(it)=alfaxz2(it) !tg theta-xz
4127     alfaxz3_nt(it)=alfaxz3(it) !1/r
4128     enddo
4129     endif
4130     nclouds_yz_nt=nclouds_yz
4131     nclouds_xz_nt=nclouds_xz
4132     if(nclouds_yz.ne.0)then
4133     nnn=0
4134     do iyz=1,nclouds_yz
4135     ptcloud_yz_nt(iyz)=ptcloud_yz(iyz)
4136     alfayz1_av_nt(iyz)=alfayz1_av(iyz)
4137     alfayz2_av_nt(iyz)=alfayz2_av(iyz)
4138     nnn=nnn+ptcloud_yz(iyz)
4139     enddo
4140 pam-fi 1.39 do ipt=1,min(ndblt_max_nt,nnn)
4141 pam-fi 1.9 db_cloud_nt(ipt)=db_cloud(ipt)
4142     enddo
4143     endif
4144     if(nclouds_xz.ne.0)then
4145     nnn=0
4146     do ixz=1,nclouds_xz
4147     ptcloud_xz_nt(ixz)=ptcloud_xz(ixz)
4148     alfaxz1_av_nt(ixz)=alfaxz1_av(ixz)
4149     alfaxz2_av_nt(ixz)=alfaxz2_av(ixz)
4150     alfaxz3_av_nt(ixz)=alfaxz3_av(ixz)
4151     nnn=nnn+ptcloud_xz(ixz)
4152     enddo
4153 pam-fi 1.39 do ipt=1,min(ntrpt_max_nt,nnn)
4154 pam-fi 1.9 tr_cloud_nt(ipt)=tr_cloud(ipt)
4155     enddo
4156     endif
4157     endif
4158     end
4159    
  ViewVC Help
Powered by ViewVC 1.1.23