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