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