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

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

Parent Directory Parent Directory | Revision Log Revision Log

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