/[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.26 - (hide annotations) (download)
Thu May 24 16:45:48 2007 UTC (17 years, 6 months ago) by pam-fi
Branch: MAIN
Changes since 1.25: +130 -45 lines
several upgrades

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

  ViewVC Help
Powered by ViewVC 1.1.23