/[PAMELA software]/tracker/ground/source/analysis/momanhough.F
ViewVC logotype

Annotation of /tracker/ground/source/analysis/momanhough.F

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.2 - (hide annotations) (download)
Mon Mar 20 19:43:33 2006 UTC (18 years, 8 months ago) by pam-fi
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +320 -299 lines
Some subroutines grouped in different files

1 pam-fi 1.1 *************************************************************************
2     * Program analysis.f
3     *
4     * - reads cluster information (LEVEL1, reduction.f output ntuple)
5     * - perform track identification and fit
6     * - create LEVEL2 ntuple
7     *
8     *************************************************************************
9     program momanhough
10    
11     include '../common/commontracker.f'
12     include '../common/common_momanhough.f'
13     c include '../common/common_level2debug.f'
14     include '../common/common_mech.f'
15     include '../common/common_xyzPAM.f'
16     include '../common/common_mini_2.f'
17     include '../common/calib.f'
18     include '../common/level1.f'
19     include '../common/level2.f'
20    
21 pam-fi 1.2 c include '../common/momanhough_init.f'
22 pam-fi 1.1
23     * flag to set debug mode
24     logical DEBUG
25     common/dbg/DEBUG
26     logical DBG_FILLED
27     data DBG_FILLED/.false./
28    
29     * flag to chose PFA
30     character*10 PFA
31     common/FINALPFA/PFA
32    
33    
34     c parameter (inf=1.e8) !just a huge number...
35    
36     * external functions
37     external npl
38     external acoordsi,coordsi,nld,coord,dcoord
39     c------------------------------------------------------------------------
40     c
41     c local variables
42     c
43     c------------------------------------------------------------------------
44    
45     character*24 processing_date
46    
47     parameter (lun_data_level1=71) !data file id number
48     parameter (lun_data_level2=72) !data file id number
49     parameter (lun_data_calib=74) !data file id number
50    
51     character*74 data_file !data file name
52     character*74 data_dir !data file name
53     character*74 data_file_calib
54     character*74 data_file_level1
55     character*74 data_file_level2
56     # ifndef TEST2003
57     parameter(ncalibmax=50)
58     character*40 file_calib(ncalibmax)
59     parameter(lun_calib_list=66) !calibration list file id
60     integer which_calib_last
61     # endif
62    
63     integer minevent !first event to be analysed
64 pam-fi 1.2 c logical FIMAGE !
65 pam-fi 1.1
66    
67     COMMON/QUEST/IQUEST(100)
68     c !permette di ottenere ntuple funzionanti nonostante
69     c ! il messaggio dei 64K di RZOUT...!???
70    
71     c*****************************************************
72     cccccc 11/9/2005 modified by david fedele
73     c swcode=202
74     cccccc 12/10/2005 modified by Elena Vannuccini
75     swcode=300
76     c****************************************************
77     c------------------------------------------------------------------------
78     c
79     c HBOOK initialization
80     c
81     c------------------------------------------------------------------------
82    
83     call HLIMIT(NWPAWC)
84    
85    
86     c------------------------------------------------------------------------
87     c
88     c reads input informations
89     c
90     c------------------------------------------------------------------------
91     call fdate(processing_date)
92     write(*,101)
93     $ processing_date
94     101 format(/
95     $ ,'*** *** *** *** *** *** *** *** *** *** *** *** ***',/
96     $ ,'* *',/
97     $ ,'* ANALYSIS *',/
98     $ ,'* *',/
99     $ ,'*** *** *** *** *** *** *** *** *** *** *** *** ***',/
100     $ ,a24,/
101     $ )
102    
103     111 format(a)
104     print*,'Data directory:'
105     read(*,111)data_dir
106     print*,data_dir
107     print*,'File identifier: (DATE_NUM)'
108     read(*,*)data_file
109     print*,data_file
110     minevent=1
111     print*,'Maximum number of events to be analized:' !20000
112     read(*,*) ntotev
113     print*,ntotev
114     print*,'Position-finding algorythm: (GOG2,ETA2)'
115     read(*,*)PFA
116     PFA=PFA(1:lnblnk(PFA))
117     print*,PFA
118     print*,'DEBUG mode: (T, F)'
119     11 format(l1)
120     read(*,11)DEBUG
121     print*,DEBUG
122     print*,'---------------------------------------------------'
123    
124    
125     ****** INITIALIZATIONS *************************************
126    
127     * 1) read charge-correlation parameters
128     print*,' '
129     print*,'- read charge-correlation parameters'
130     print*,' '
131     call readchargeparam
132    
133     * 1) read mip parameters
134     print*,' '
135     print*,'- read mip parameters'
136     print*,' '
137     call readmipparam
138    
139     * 2) read z coordinates of the planes
140     print*,' '
141     print*,'- read z coordinates of the planes'
142     print*,' '
143     call mech_sensor !reads sensors centres coordinates
144     do ip=1,nplanes
145     fitz(ip)=z_mech_sensor(ip,1,1)*0.1 !cm
146     * gets planes mechanical z positions
147     * (in mm) and sets them in micrometers
148     enddo
149    
150    
151     * 3) read eta PFA parameters
152     print*,' '
153     print*,'- read P.F.A. parameters'
154     print*,' '
155     call readetaparam
156    
157     print*,' '
158     print*,'- First guess P.F.A. >>>> ',PFAdef
159     print*,' '
160     print*,'- Final P.F.A. >>>> ',PFA
161     print*,' '
162    
163     * 4) read magnetic field map
164     print*,' '
165     print*,'- read magnetic field map'
166     print*,' '
167     c call read_B(5) !legge il nome da STI!!.. temporaneo
168     c call read_B_2maps(5) !legge il nome da STI!!.. temporaneo
169     call read_B
170    
171     * 5) read allignment parameters
172     print*,' '
173     print*,'- read aligment parameters'
174     print*,' '
175     call readalignparam
176    
177     ************************************************************
178    
179    
180    
181    
182    
183     c------------------------------------------------------------------------
184     c
185     c LEVEL 2 ntuple booking
186     c
187     c------------------------------------------------------------------------
188    
189    
190     503 format(a,'DW_',a,'_level2.rz')
191     write(data_file_level2,503)
192     $ data_dir(1:LNBLNK(data_dir))
193     $ ,data_file(1:LNBLNK(data_file))
194    
195     print*,''
196     print*,'------------------------------------'
197     print*,' Creating LEVEL2 rz file'
198     print*,' ',data_file_level2
199     print*,'------------------------------------'
200    
201    
202     C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
203     C largest RZ file: IQUEST(10) records x LREC words x 4 byte
204     C with the following settings: 65000 x 4096 x 4 = 1G
205     C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
206     IQUEST(10)=65000
207     c !permette di ottenere ntuple funzionanti nonostante
208     c ! il messaggio dei 64K di RZOUT...!???
209     call HROPEN(lun_data_level2,
210     $ 'LEVEL2',data_file_level2,'QNP',4096,istat) !opens rz
211     if(istat.ne.0) goto 19
212     call book_level2
213    
214     if(DEBUG)then
215     print*,'(creating DEBUG nt-uple and histos)'
216     call book_debug
217     endif
218    
219    
220    
221     c------------------------------------------------------------------------
222     c
223     c fills bad variables with online DSP bad strips from calib histograms
224     c for bad strip exclusion
225     c
226     c------------------------------------------------------------------------
227    
228     # ifdef TEST2003
229     c------------------------------------------------------------------------
230     c
231     c opens calib file
232     c
233     c------------------------------------------------------------------------
234    
235     505 format(a,'output_',a,'_calib.rz')
236     write(data_file_calib,505)
237     $ data_dir(1:LNBLNK(data_dir))
238     $ ,data_file(1:LNBLNK(data_file))
239    
240     print*,' '
241     print*,'OPENING CALIB FILE...', data_file_calib
242    
243     IQUEST(10)=65000
244     call HROPEN(lun_data_calib,'IN',data_file_calib,'QP',4096,istat)
245     if(istat.ne.0) goto 17
246     call HCDIR('//IN',' ')
247    
248     call HRIN(0,9999,0) !puts histograms in memory
249    
250     print*,' '
251    
252     print*,' '
253     print*,'READING PEDESTAL, SIGMA AND BADSTRIP HISTOGRAMS...'
254     print*,' '
255    
256     call fillpedsig
257    
258     call HREND('IN')
259     close (lun_data_calib)
260     # else
261    
262     print*,' '
263     print*,'OPENING CALIBRATION-LIST FILE:'
264     501 format(a,'DW_',a,'_calib.txt')
265     write(data_file_calib,501)
266     $ data_dir(1:LNBLNK(data_dir))
267     $ ,data_file(1:LNBLNK(data_file))
268     print*,data_file_calib
269     open(lun_calib_list,
270     $ FILE=data_file_calib(1:LNBLNK(data_file_calib))
271     $ ,STATUS='UNKNOWN'
272     $ ,IOSTAT=iostat
273     $ )
274    
275     113 format(i5,' ',a25)
276    
277     do i=1,ncalibmax
278    
279     read(lun_calib_list,113,IOSTAT=iostat) n_cal_list
280     $ ,file_calib(i)(1:LNBLNK(file_calib(i)))
281     if(iostat.ne.0)then
282     ncal=i-1
283     goto 2000
284     endif
285     print*,n_cal_list,' - '
286     $ ,file_calib(i)(1:LNBLNK(file_calib(i)))
287    
288     enddo
289    
290     2000 close(lun_calib_list)
291     which_calib_last=0
292    
293     # endif
294    
295     c------------------------------------------------------------------------
296     c
297     c opens level1 file
298     c
299     c------------------------------------------------------------------------
300    
301     # ifdef TEST2003
302     504 format(a,'output_',a,'_level1.rz')
303     # else
304     504 format(a,'DW_',a,'_level1.rz')
305     # endif
306     write(data_file_level1,504)
307     $ data_dir(1:LNBLNK(data_dir))
308     $ ,data_file(1:LNBLNK(data_file))
309     print*,''
310     print*,'OPENING LEVEL1 FILE:'
311     print*,data_file_level1
312     IQUEST(10)=65000
313     c !permette di ottenere ntuple funzionanti nonostante
314     c ! il messaggio dei 64K di RZOUT...!???
315     call HROPEN(lun_data_level1,
316     $ 'LEVEL1',data_file_level1,'QP',4096,istat) !opens rz
317     if(istat.ne.0) goto 19
318     call HRIN(ntp_level1,9999,20)
319     call access_level1
320     c call HPRNTU(ntp_level1+20)
321     call HNOENT(ntp_level1+20,iemax0)
322     print*,' events',iemax0
323    
324    
325     ************************************************************
326     ************************************************************
327     ************************************************************
328     *
329     * start track analysis
330     *
331     ************************************************************
332     ************************************************************
333     ************************************************************
334     maxevent=minevent+ntotev-1
335     do iev = minevent,MIN(iemax0,maxevent) !loop on events
336    
337     call HCDIR('//LEVEL1',' ')
338     call HGNT(ntp_level1+20,iev,ierr) !reads an event
339     if(ierr.ne.0) goto 21
340     *------------------------------------------------------
341     * LEVEL2 N-TUPLE INITIALIZATIONS
342     call init_level2
343     if(DEBUG)call init_level2_debug
344     if(DEBUG)DBG_FILLED=.false.
345    
346     if(.not.good1)then
347     goto 8800 !fill nt-uple and go to next event
348     endif
349     *------------------------------------------------------
350    
351    
352     # ifndef TEST2003
353    
354     if(which_calib.ne.which_calib_last.and.
355     $ which_calib.ne.0)then
356    
357     data_file_calib=
358     $ data_dir(1:LNBLNK(data_dir))//
359     $ file_calib(which_calib)
360     $ (1:LNBLNK(file_calib(which_calib)))
361     c print*,data_file_calib
362     print*,''
363     print*,
364     $ '@ event ',nev2
365     $ ,' (CPU pkt N.',pkt_num1,')'
366     print*,'--> ',data_file_calib
367     IQUEST(10)=65000
368     call HROPEN(lun_data_calib,
369     $ 'CALIB',data_file_calib,'QP',4096,istat) !opens
370     if(istat.ne.0) goto 19
371     call HRIN(0,9999,0)
372     call fillpedsig
373     do iview=1,nviews
374     call HDELET(id_hi_bad+iview)
375     call HDELET(id_hi_ped+iview)
376     call HDELET(id_hi_sig+iview)
377     enddo
378     call HREND('CALIB')
379     close(lun_data_calib)
380     which_calib_last=which_calib
381    
382     elseif(which_calib.eq.0)then
383    
384     nocalib=nocalib+1
385     good2=.false.
386     goto 8800 !fill nt-uple and go to next event
387    
388     endif
389    
390     # endif
391    
392    
393     *------------------------------------------------------
394     * cut on maximum number of clusters
395     *------------------------------------------------------
396     if(nclstr1.gt.nclstrmax_level2)then
397     goto 8800 !fill nt-uple and go to next event
398     endif
399    
400     do i=1,nclstr1
401     cl_used(i)=0 !init mask of clusters associated to a track
402     enddo
403    
404     if(DEBUG)then
405     print*,'----------------------------------'
406     print*,iev,' ** ',nev2
407     endif
408    
409 pam-fi 1.2 * ///////////////////////////////////////////////
410     * \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
411 pam-fi 1.1
412     iflag=0
413 pam-fi 1.2 call track_finding(iflag)
414 pam-fi 1.1 if(iflag.eq.1)then !bad event
415 pam-fi 1.2 goto 880 !fill ntp and go to next event
416     endif
417     c$$$*-------------------------------------------------------------------------------
418     c$$$* STEP 1
419     c$$$*-------------------------------------------------------------------------------
420     c$$$* X-Y cluster association
421     c$$$*
422     c$$$* Clusters are associated to form COUPLES
423     c$$$* Clusters not associated in any couple are called SINGLETS
424     c$$$*
425     c$$$* Track identification (Hough transform) and fitting is first done on couples.
426     c$$$* Hence singlets are possibly added to the track.
427     c$$$*
428     c$$$* Variables assigned by the routine "cl_to_couples" are those in the
429     c$$$* common blocks:
430     c$$$* - common/clusters/cl_good
431     c$$$* - common/couples/clx,cly,ncp_plane,ncp_tot,cp_useds1,cp_useds2
432     c$$$* - common/singlets/ncls,cls,cl_single
433     c$$$*-------------------------------------------------------------------------------
434     c$$$*-------------------------------------------------------------------------------
435     c$$$
436     c$$$ iflag=0
437     c$$$ call cl_to_couples(iflag)
438     c$$$ if(iflag.eq.1)then !bad event
439     c$$$ goto 880 !fill ntp and go to next event
440     c$$$ endif
441     c$$$
442     c$$$*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
443     c$$$* selezione di tracce pulite per diagnostica
444     c$$$*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
445     c$$$c$$$ if(DEBUG)then
446     c$$$c$$$ do ip=1,nplanes
447     c$$$c$$$ if(ncp_plane(ip).ne.1)good2=.false.
448     c$$$c$$$ enddo
449     c$$$c$$$c if(good2.eq.0)goto 100!next event
450     c$$$c$$$c if(good2.eq.0)goto 880!fill ntp and go to next event
451     c$$$c$$$ endif
452     c$$$
453     c$$$
454     c$$$
455     c$$$*-----------------------------------------------------
456     c$$$*-----------------------------------------------------
457     c$$$* HOUGH TRASFORM
458     c$$$*-----------------------------------------------------
459     c$$$*-----------------------------------------------------
460     c$$$
461     c$$$
462     c$$$*-------------------------------------------------------------------------------
463     c$$$* STEP 2
464     c$$$*-------------------------------------------------------------------------------
465     c$$$*
466     c$$$* Association of couples to form
467     c$$$* - DOUBLETS in YZ view
468     c$$$* - TRIPLETS in XZ view
469     c$$$*
470     c$$$* Variables assigned by the routine "cp_to_doubtrip" are those in the
471     c$$$* common blocks:
472     c$$$* - common/hough_param/
473     c$$$* $ alfayz1, !Y0
474     c$$$* $ alfayz2, !tg theta-yz
475     c$$$* $ alfaxz1, !X0
476     c$$$* $ alfaxz2, !tg theta-xz
477     c$$$* $ alfaxz3 !1/r
478     c$$$* - common/doublets/ndblt,cpyz1,cpyz2
479     c$$$* - common/triplets/ntrpt,cpxz1,cpxz2,cpxz3
480     c$$$*-------------------------------------------------------------------------------
481     c$$$*-------------------------------------------------------------------------------
482     c$$$
483     c$$$ iflag=0
484     c$$$ call cp_to_doubtrip(iflag)
485     c$$$ if(iflag.eq.1)then !bad event
486     c$$$ goto 880 !fill ntp and go to next event
487     c$$$ endif
488     c$$$
489     c$$$
490     c$$$*-------------------------------------------------------------------------------
491     c$$$* STEP 3
492     c$$$*-------------------------------------------------------------------------------
493     c$$$*
494     c$$$* Classification of doublets and triplets to form CLOUDS,
495     c$$$* according to distance in parameter space.
496     c$$$*
497     c$$$* cloud = cluster of points (doublets/triplets) in parameter space
498     c$$$*
499     c$$$*
500     c$$$*
501     c$$$* Variables assigned by the routine "doub_to_YZcloud" are those in the
502     c$$$* common blocks:
503     c$$$* - common/clouds_yz/
504     c$$$* $ nclouds_yz
505     c$$$* $ ,alfayz1_av,alfayz2_av
506     c$$$* $ ,ptcloud_yz,db_cloud,cpcloud_yz
507     c$$$*
508     c$$$* Variables assigned by the routine "trip_to_XZcloud" are those in the
509     c$$$* common blocks:
510     c$$$* common/clouds_xz/
511     c$$$* $ nclouds_xz xz2_av,alfaxz3_av
512     c$$$* $ ,ptcloud_xz,tr_cloud,cpcloud_xz
513     c$$$*-------------------------------------------------------------------------------
514     c$$$*-------------------------------------------------------------------------------
515     c$$$
516     c$$$ iflag=0
517     c$$$ call doub_to_YZcloud(iflag)
518     c$$$ if(iflag.eq.1)then !bad event
519     c$$$ goto 880 !fill ntp and go to next event
520     c$$$ endif
521     c$$$ iflag=0
522     c$$$ call trip_to_XZcloud(iflag)
523     c$$$ if(iflag.eq.1)then !bad event
524     c$$$ goto 880 !fill ntp and go to next event
525 pam-fi 1.1 c$$$ endif
526 pam-fi 1.2 c$$$
527 pam-fi 1.1 c*****************************************************
528     cccccc 01/12/2005 modified by elena
529     if(DEBUG)then
530     call fill_level2_clouds
531     call HCDIR('//LEVEL2',' ')
532     call HFNT(ntp_level2+1) !fill DEBUG nt-uple
533     DBG_FILLED=.true.
534     endif
535     c*****************************************************
536    
537    
538     iflag=0
539 pam-fi 1.2 call track_fitting(iflag)
540     if(iflag.eq.1)then !bad event
541     goto 880 !fill ntp and go to next event
542 pam-fi 1.1 endif
543    
544    
545 pam-fi 1.2 c$$$
546     c$$$*-------------------------------------------------------------------------------
547     c$$$* STEP 4 (ITERATED until any other physical track isn't found)
548     c$$$*-------------------------------------------------------------------------------
549     c$$$*
550     c$$$* YZ and XZ clouds are combined in order to obtain the initial guess
551     c$$$* of the candidate-track parameters.
552     c$$$* A minimum number of matching couples between YZ and XZ clouds is required.
553     c$$$*
554     c$$$* A TRACK CANDIDATE is defined by
555     c$$$* - the couples resulting from the INTERSECTION of the two clouds, and
556     c$$$* - the associated track parameters (evaluated by performing a zero-order
557     c$$$* track fitting)
558     c$$$*
559     c$$$* The NTRACKS candidate-track parameters are stored in common block:
560     c$$$*
561     c$$$* - common/track_candidates/NTRACKS,AL_STORE
562     c$$$* $ ,XV_STORE,YV_STORE,ZV_STORE
563     c$$$* $ ,XM_STORE,YM_STORE,ZM_STORE
564     c$$$* $ ,RESX_STORE,RESY_STORE
565     c$$$* $ ,AXV_STORE,AYV_STORE
566     c$$$* $ ,XGOOD_STORE,YGOOD_STORE
567     c$$$* $ ,CP_STORE,RCHI2_STORE
568     c$$$*
569     c$$$*-------------------------------------------------------------------------------
570     c$$$*-------------------------------------------------------------------------------
571     c$$$ ntrk=0 !counter of identified physical tracks
572     c$$$
573     c$$$11111 continue !<<<<<<< come here when performing a new search
574     c$$$
575     c$$$ iflag=0
576     c$$$ call clouds_to_ctrack(iflag)
577     c$$$ if(iflag.eq.1)then !no candidate tracks found
578     c$$$ goto 880 !fill ntp and go to next event
579     c$$$ endif
580     c$$$
581     c$$$ FIMAGE=.false. !processing best track (not track image)
582     c$$$ ibest=0 !best track among candidates
583     c$$$ iimage=0 !track image
584     c$$$* ------------- select the best track -------------
585     c$$$ rchi2best=1000000000.
586     c$$$ do i=1,ntracks
587     c$$$ if(RCHI2_STORE(i).lt.rchi2best.and.
588     c$$$ $ RCHI2_STORE(i).gt.0)then
589     c$$$ ibest=i
590     c$$$ rchi2best=RCHI2_STORE(i)
591     c$$$ endif
592     c$$$ enddo
593     c$$$ if(ibest.eq.0)goto 880 !>> no good candidates
594     c$$$*-------------------------------------------------------------------------------
595     c$$$* The best track candidate (ibest) is selected and a new fitting is performed.
596     c$$$* Previous to this, the track is refined by:
597     c$$$* - possibly adding new COUPLES or SINGLETS from the missing planes
598     c$$$* - evaluating the coordinates with improved PFAs
599     c$$$* ( angle-dependent ETA algorithms )
600     c$$$*-------------------------------------------------------------------------------
601     c$$$
602     c$$$ 1212 continue !<<<<< come here to fit track-image
603     c$$$
604     c$$$ if(.not.FIMAGE)then !processing best candidate
605     c$$$ icand=ibest
606     c$$$ else !processing image
607     c$$$ icand=iimage
608     c$$$ iimage=0
609     c$$$ endif
610     c$$$ if(icand.eq.0)then
611     c$$$ print*,'HAI FATTO UN CASINO!!!!!! icand = ',icand
612     c$$$ $ ,ibest,iimage
613     c$$$ return
614     c$$$ endif
615     c$$$
616     c$$$* *-*-*-*-*-*-*-*-*-*-*-*-*-*-*
617     c$$$ call refine_track(icand)
618     c$$$* *-*-*-*-*-*-*-*-*-*-*-*-*-*-*
619     c$$$
620     c$$$* **********************************************************
621     c$$$* ************************** FIT *** FIT *** FIT *** FIT ***
622     c$$$* **********************************************************
623     c$$$ do i=1,5
624     c$$$ AL(i)=dble(AL_STORE(i,icand))
625     c$$$ enddo
626     c$$$ ifail=0 !error flag in chi2 computation
627     c$$$ jstep=0 !# minimization steps
628     c$$$
629     c$$$ call mini_2(jstep,ifail)
630     c$$$ if(ifail.ne.0) then
631     c$$$ if(DEBUG)then
632     c$$$ print *,
633     c$$$ $ '*** MINIMIZATION FAILURE *** (mini_2) '
634     c$$$ $ ,iev
635     c$$$ endif
636     c$$$ chi2=-chi2
637     c$$$ endif
638     c$$$
639     c$$$ if(DEBUG)then
640     c$$$ print*,'----------------------------- improved track coord'
641     c$$$22222 format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)
642     c$$$ do ip=1,6
643     c$$$ write(*,22222)ip,zm(ip),xm(ip),ym(ip)
644     c$$$ $ ,xm_A(ip),ym_A(ip),xm_B(ip),ym_B(ip)
645     c$$$ $ ,xgood(ip),ygood(ip),resx(ip),resy(ip)
646     c$$$ enddo
647     c$$$ endif
648     c$$$
649     c$$$c rchi2=chi2/dble(ndof)
650     c$$$ if(DEBUG)then
651     c$$$ print*,' '
652     c$$$ print*,'****** SELECTED TRACK *************'
653     c$$$ print*,'# R. chi2 RIG'
654     c$$$ print*,' --- ',chi2,' --- '
655     c$$$ $ ,1./abs(AL(5))
656     c$$$ print*,'***********************************'
657     c$$$ endif
658     c$$$* **********************************************************
659     c$$$* ************************** FIT *** FIT *** FIT *** FIT ***
660     c$$$* **********************************************************
661     c$$$
662     c$$$
663     c$$$* ------------- search if the track has an IMAGE -------------
664     c$$$* ------------- (also this is stored ) -------------
665     c$$$ if(FIMAGE)goto 122 !>>> jump! (this is already an image)
666     c$$$* now search for track-image, by comparing couples IDs
667     c$$$ do i=1,ntracks
668     c$$$ iimage=i
669     c$$$ do ip=1,nplanes
670     c$$$ if( CP_STORE(nplanes-ip+1,icand).ne.
671     c$$$ $ -1*CP_STORE(nplanes-ip+1,i) )iimage=0
672     c$$$ enddo
673     c$$$ if( iimage.ne.0.and.
674     c$$$c $ RCHI2_STORE(i).le.CHI2MAX.and.
675     c$$$c $ RCHI2_STORE(i).gt.0.and.
676     c$$$ $ .true.)then
677     c$$$ if(DEBUG)print*,'Track candidate ',iimage
678     c$$$ $ ,' >>> TRACK IMAGE >>> of'
679     c$$$ $ ,ibest
680     c$$$ goto 122 !image track found
681     c$$$ endif
682     c$$$ enddo
683     c$$$ 122 continue
684     c$$$
685     c$$$* --- and store the results --------------------------------
686     c$$$ ntrk = ntrk + 1 !counter of found tracks
687     c$$$ if(.not.FIMAGE
688     c$$$ $ .and.iimage.eq.0) image(ntrk)= 0
689     c$$$ if(.not.FIMAGE
690     c$$$ $ .and.iimage.ne.0)image(ntrk)=ntrk+1 !this is the image of the next
691     c$$$ if(FIMAGE) image(ntrk)=ntrk-1 !this is the image of the previous
692     c$$$
693     c$$$ call fill_level2_tracks(ntrk) !==> good2=.true.
694     c$$$c print*,'++++++++++ iimage,fimage,ntrk,image '
695     c$$$c $ ,iimage,fimage,ntrk,image(ntrk)
696     c$$$
697     c$$$ if(ntrk.eq.NTRKMAX)then
698     c$$$ if(DEBUG)
699     c$$$ $ print*,
700     c$$$ $ '** warning ** number of identified '//
701     c$$$ $ 'tracks exceeds vector dimension '
702     c$$$ $ ,'( ',NTRKMAX,' )'
703     c$$$cc good2=.false.
704     c$$$ goto 880 !fill ntp and go to next event
705     c$$$ endif
706     c$$$ if(iimage.ne.0)then
707     c$$$ FIMAGE=.true. !
708     c$$$ goto 1212 !>>> fit image-track
709     c$$$ endif
710     c$$$
711     c$$$* --- then remove selected clusters (ibest+iimage) from clouds ----
712     c$$$ call clean_XYclouds(ibest)
713     c$$$ if(iflag.eq.1)then !bad event
714     c$$$ goto 880 !fill ntp and go to next event
715     c$$$ endif
716     c$$$
717     c$$$* **********************************************************
718     c$$$* condition to start a new search
719     c$$$* **********************************************************
720     c$$$ ixznew=0
721     c$$$ do ixz=1,nclouds_xz
722     c$$$ if(ptcloud_xz(ixz).ge.nptxz_min)ixznew=1
723     c$$$ enddo
724     c$$$ iyznew=0
725     c$$$ do iyz=1,nclouds_yz
726     c$$$ if(ptcloud_yz(iyz).ge.nptyz_min)iyznew=1
727     c$$$ enddo
728     c$$$
729     c$$$ if(ixznew.ne.0.and.
730     c$$$ $ iyznew.ne.0.and.
731     c$$$ $ rchi2best.le.CHI2MAX.and.
732     c$$$c $ rchi2best.lt.15..and.
733     c$$$ $ .true.)then
734     c$$$ if(DEBUG)then
735     c$$$ print*,'***** NEW SEARCH ****'
736     c$$$ endif
737     c$$$ goto 11111 !try new search
738     c$$$
739     c$$$ endif
740     c$$$* **********************************************
741     c$$$
742 pam-fi 1.1 * >>>>>>>>>>>>>>>>>>> NT-UPLE filling <<<<<<<<<<<<<<<<<<<<
743     * + + + + + + + + + + + + + + + + +
744     * / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ /
745     * + + + + + + + + + + + + + + + + +
746    
747     880 continue
748    
749     * **********************************************************
750     * stores info about clusters not associated with any track
751     * **********************************************************
752    
753     c*****************************************************
754     cccccc 27/9/2005 modified by david fedele
755     * count #cluster per plane not associated to any track
756     c$$$ do icl=1,nclstr1
757     c$$$ if(cl_used(icl).eq.0)then !cluster not included in any track
758     c$$$ ip=nplanes-npl(VIEW(icl))+1
759     c$$$ if(mod(VIEW(icl),2).eq.0)nclsx(ip)=nclsx(ip)+1
760     c$$$ if(mod(VIEW(icl),2).eq.1)nclsy(ip)=nclsy(ip)+1
761     c$$$ endif
762     c$$$c print*,icl,cl_used(icl),cl_good(icl),ip,VIEW(icl)!nclsx(ip),nclsy(ip)
763     c$$$ enddo
764     c**********************************************************************
765    
766     call fill_level2_siglets
767    
768     c print*,'****** ',iev,' - ',ntrk,nclsx,nclsy
769     c print*,'****** ',iev,' - ',ntrk,' --- ',(BdL(i),i=1,ntrk)
770    
771     if(DEBUG)then
772    
773     print*,''
774     print*,'DONE!'
775     print*,''
776     print*,'* summary *'
777     print*,'tracks ',ntrk
778     print*,'cl used ',(cl_used(i),i=1,nclstr1)
779     c*****************************************************
780     c$$$cccccc 27/9/2005 modified by david fedele
781     c$$$ print*,'cl unused (x-y)'
782     c$$$ do ip=1,nplanes
783     c$$$ print*,ip,' << ',nclsx(ip),nclsy(ip)
784     c$$$ enddo
785     c******************************************************
786     print*,''
787     print*,''
788     endif
789    
790     8800 continue
791    
792 pam-fi 1.2
793     * ///////////////////////////////////////////////
794     * \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
795    
796 pam-fi 1.1 call HCDIR('//LEVEL2',' ')
797     call HFNT(ntp_level2) !fill LEVEL2 nt-uple
798     if(.not.DBG_FILLED.and.DEBUG)
799     $ print*,'@@@@@@@@@@@@',ntp_level2+1,nev2_nt
800     if(.not.DBG_FILLED.and.DEBUG)call HFNT(ntp_level2+1)
801    
802     100 continue
803     enddo !end loop on events
804    
805     c------------------------------------------------------------------------
806     c
807     c no error exit
808     c
809     c------------------------------------------------------------------------
810    
811     c$$$ print*,' '
812     c$$$ print*,'REDUCTION SUCCESSFULLY COMPLETED'
813     c$$$ print*,' '
814     c$$$ print*,' '
815    
816     goto 9000 !happy ending
817    
818    
819     c------------------------------------------------------------------------
820     c
821     c data file opening error
822     c
823     c------------------------------------------------------------------------
824     19 continue
825    
826     print*,' '
827     print*,'ERROR OPENING DATA FILE: ',data_file
828     print*,' '
829     print*,' '
830    
831     goto 9000 !the end
832    
833     c------------------------------------------------------------------------
834     c
835     c level1 ntuple event reading error
836     c
837     c------------------------------------------------------------------------
838    
839     21 continue
840    
841     print*,' '
842     print*,'ERROR WHILE READING LEVEL1 NTUPLE, AT EVENT
843     $ : ',iev
844     print*,' '
845     print*,' '
846    
847     goto 9000 !the end
848    
849     c------------------------------------------------------------------------
850     c
851     c calib file opening error
852     c
853     c------------------------------------------------------------------------
854    
855     17 continue
856    
857     print*,' '
858     print*,'preanalysis: ERROR OPENING INPUT FILE: ',data_file_calib
859     print*,' '
860     print*,' '
861    
862     goto 9000 !the end
863    
864     c------------------------------------------------------------------------
865     c
866     c closes files and exits
867     c
868     c------------------------------------------------------------------------
869    
870     9000 continue
871    
872     if(DEBUG)call HPRNTU(ntp_level2+1)
873     call HPRNTU(ntp_level2)
874     print*,''
875     call HCDIR('//LEVEL2',' ')
876     call HROUT(ntp_level2,ICYCLE,'T')
877     print *,'- Stored LEVEL2 nt-uple (',nev2,' entries )'
878     if(DEBUG)call HROUT(ntp_level2+1,ICYCLE,'T')
879     if(DEBUG)print *,'- Stored DEBUG nt-uple '
880     call HREND('level2')
881     close(lun_data_level2)
882    
883     call HCDIR('//LEVEL1',' ')
884     call HREND('level1')
885     close(lun_data_level1)
886    
887     stop
888     end
889    
890    
891     ************************************************************
892    
893    
894     # include "momanhough-subroutines.F"
895    
896    

  ViewVC Help
Powered by ViewVC 1.1.23