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