/[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.32 - (hide annotations) (download)
Tue Sep 4 15:44:49 2007 UTC (17 years, 3 months ago) by bongi
Branch: MAIN
CVS Tags: v4r00
Changes since 1.31: +46 -38 lines
real to double precision change

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

  ViewVC Help
Powered by ViewVC 1.1.23