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