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