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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1.1.1 - (show 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 ************************************************************
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