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