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