cd188f0d775c821d6098768de5637269dc632811
[u/mrichter/AliRoot.git] / MUON / reco_muon.F
1 *  25/5/99
2 **  Authors J.P. Cussonneau & P. Lautridou 
3
4 ************** commentaires **********************
5
6 * bon muon = tous les hits de la trace (pas forcement 28) proviennent de muons
7 * ghost = bonne trace dans laquelle certains hits proviennent délectrons, de muons ou sont ambigus
8  
9 * ALPHATOP
10 * alpha to P
11
12 * EFF
13 * efficacite
14
15 * HHIT 
16 * h au vertex
17
18 * HTOP
19 * h to P
20
21 * INDEXMAX
22 * nombre de candidats a ordonner (a partir de l'ímpulsion)
23
24 * INDEXTAB
25 * pour recuperer les candidats
26
27 * ISTAT
28 * =1 si bon muon, =2 si ghost, =0 autrement
29
30 * ITCHECK
31 * =1 si bonne trace, =0 autrement
32
33 * IT_LIST
34 * permet a partir du numero de la trace de retrouver la numero du hit
35
36 * IT_NP
37 * compte le nombre de plans touches par trace
38
39 * ITRACK
40 * permet de retrouver le numero de la trace a partir du numero de hit
41
42 * ITTROUGH
43 * pour une trace et une station donnee, dit si la trace est passee dans la chambre
44
45 * IVERTEX
46 * =0 si point (0,0) du vertex impose, sinon coordonnes libres (pour le fit)
47
48 * JCAN 
49 * numero du hit pour une station et un candidat
50
51 * JCANTYP
52 * nombre de hits par trace au niveau des stations 4 et 5 (3 ou 4)
53
54 * JJOUT
55 * numero de hit associe a une chambre et a une trace
56
57 * NMUONALL 
58 * nombre de bons muons trouves
59
60 * NERR
61 * = nombre de fois ou lón ná pas trouve la bonne trace dans la station
62
63 * NERRALL
64 * nombre de cas ou pas de hit trouve par station
65
66 * NGHOSTALL
67 * nombre de fantomes trouves
68
69 * NRES
70 * nombre de resonances dans lácceptance
71
72 * NRESF
73 * nombre de resonances trouvees
74
75 * NTRACFALL
76 * nombre de traces totales trouvees
77
78 * NTRMUALL
79 * nombre total de muons dans lácceptance
80
81
82 ****************************************************************
83       SUBROUTINE reconstmuon(IFIT,IDEBUGC,NEV,IDRES,IREADGEANT)
84 ****************************************************************
85       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
86
87       COMMON/DEBEVT/IDEBUG
88
89       IDEBUG=IDEBUGC
90
91 ** Read events          
92       CALL RECO_READEVT(NEV,IDRES,IREADGEANT)
93
94
95 ** Trackfinding         
96       CALL RECO_TRACKF(IDRES,IREADGEANT)
97
98 ** Precision fit  
99       IF (IFIT.EQ.1) THEN
100          CALL RECO_PRECISION
101       ENDIF
102
103 ** Calculate         
104       CALL RECO_SUM  
105
106       END
107
108 ********************************************
109       BLOCK DATA
110 ********************************************
111       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
112       PARAMETER(NBSTATION=5)
113 * --
114       COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
115
116       DATA DZ_PL/8.,8.,24.3,8.,8./     
117       
118 **      DATA DZ_PL/8.,8.,8.,8.,8./     
119 *      DATA DZ_PL/20.,20.,20.,20.,20./     
120 *      DATA DZ_PL/20.,20.,24.3,20.,20./  
121 *      DATA DZ_PL/20.,40.,40.,40.,40./   
122 *      DATA DZ_PL/20.,50.,50.,50.,50./  
123 *      DATA DZ_PL/20.,55.,55.,55.,55./  
124 *      DATA DZ_PL/20.,60.,60.,60.,60./ 
125 *      DATA DZ_PL/20.,80.,80.,80.,80./ 
126 *      DATA DZ_PL/20.,100.,100.,100.,100./ 
127
128       DATA ZPLANE/-511.,-686.0,-962.7,-1245.,-1445./ ! dstation=8cm
129
130 *      DATA ZPLANE/-507.5,-682.5,-967.5,-1241.5,-1441.5/ ! dstation=20cm CCC
131 *      DATA ZPLANE/-505.,-680.,-965.,-1239.,-1439./ ! dstation=20cm
132 **      DATA ZPLANE/-511.,-686.0,-971.,-1245.,-1445./ ! dstation=8cm
133 *      DATA ZPLANE/-511.,-686.0,-962.7,-1245.,-1445./ ! dstation=8cm
134 *      DATA ZPLANE/-505.,-680.,-962.7,-1239.,-1439./   ! dstation=20cm
135 *      DATA ZPLANE/-505.,-670.,-954.,-1229.,-1429./ ! dstation=40cm
136 *      DATA ZPLANE/-505.,-665.,-949.,-1224.,-1424./ ! dstation=50cm
137 *      DATA ZPLANE/-505.,-663.,-947.,-1222.,-1422./ ! dstation=50cm
138 *      DATA ZPLANE/-505.,-660.,-944.,-1219.,-1419./ ! dstation=60cm
139 *      DATA ZPLANE/-505.,-650.,-935.,-1209.,-1409./ ! dstation=80cm
140 *      DATA ZPLANE/-505.,-640.,-925.,-1199.,-1399./ ! dstation=100cm
141
142 **      DATA ZPLANE/-507.,-682.0,-950.7,-1241.,-1441./
143
144 **      DATA ZCOIL,ZMAGEND/-825.0,-1125./  ! Constant field 3 Tm
145       DATA ZCOIL,ZMAGEND/-805.0,-1233./ ! CCC magn. field map M.B 
146 **      DATA ZCOIL,ZMAGEND/-795.1,-1242.9/ ! CCC magn. field map M.B 
147
148 * --
149       END
150
151 ********************************************
152       SUBROUTINE cutpxz(spxzcut)
153 ********************************************
154       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
155
156       COMMON/ACUTPXZ/ACUTPXZ
157
158       ACUTPXZ = SPXZCUT
159
160       END
161
162 ********************************************
163       SUBROUTINE sigmacut(ssigcut)
164 ********************************************
165       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
166
167       COMMON/TRACKFI/EFF,EFF1,EFF2,XPREC,YPREC,PHIPREC,ALAMPREC,
168      &                 HCUT,LBKG,SIGCUT,ALPHATOP,HTOP
169
170 **      SIGCUT = SSIGCUT
171       SIGCUT=5.  ! CCC
172
173       END
174
175 ********************************************
176       SUBROUTINE xpreci(sxprec)
177 ********************************************
178       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
179
180       COMMON/TRACKFI/EFF,EFF1,EFF2,XPREC,YPREC,PHIPREC,ALAMPREC,
181      &                 HCUT,LBKG,SIGCUT,ALPHATOP,HTOP
182
183 **      XPREC = SXPREC
184       XPREC = 0.01 ! CCC
185  
186       END
187
188 ********************************************
189       SUBROUTINE ypreci(syprec)
190 ********************************************
191       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
192
193       COMMON/TRACKFI/EFF,EFF1,EFF2,XPREC,YPREC,PHIPREC,ALAMPREC,
194      &                 HCUT,LBKG,SIGCUT,ALPHATOP,HTOP
195
196       YPREC = SYPREC
197 **      YPREC = 0.2 ! CCC
198
199       END
200
201 *************************************************************************
202       SUBROUTINE RECO_INIT(seff,sb0,sbl3)
203 *************************************************************************
204       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
205
206       CALL TRACKF_INIT(seff,sb0,sbl3)
207
208       CALL PREC_INIT 
209
210       RETURN
211       END
212
213 *************************************************************************
214       SUBROUTINE TRACKF_INIT(seff,sb0,sbl3)
215 *************************************************************************
216       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
217 **
218       PARAMETER(NBSTATION=5,NTRMAX=500) 
219 **      
220       COMMON/REVENT/IEVBKGI,NBKGMAX,MAXUPSEV
221 **      
222       COMMON/MAGNET/BL3,B0
223 **      
224       COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
225 **
226       COMMON/FILED/FILERES,FILEBKG,FILEOUT,FILEMIN
227 **
228       COMMON/TRACKFI/EFF,EFF1,EFF2,XPREC,YPREC,PHIPREC,ALAMPREC,
229      &                 HCUT,LBKG,SIGCUT,ALPHATOP,HTOP
230 **
231       COMMON/TRACKSUM/NRES(5),NRESF,NTRMUALL,NMUONALL,NGHOSTALL,
232      &     NTRACKFALL,NERRALL(NBSTATION),IR               
233 **      
234       COMMON/ACUTPXZ/ACUTPXZ
235 **
236       COMMON/DEBEVT/IDEBUG
237 **
238       CALL HIST_CREATE
239 *
240       EFF = SEFF
241       B0 = SB0
242       BL3 = SBL3
243
244       PXZCUT = ACUTPXZ 
245       AMAGLEN = ZMAGEND-ZCOIL
246       ZM = AMAGLEN/2.+ZCOIL
247       ALPHATOP = 0.01*0.3*B0*ABS(AMAGLEN)
248       HTOP = ALPHATOP*ZM
249       HCUT = ABS(HTOP)/PXZCUT
250
251       print*,'TRACK_INIT hcut= ',hcut
252       print*,'TRACK_INIT eff = ',eff
253       print*,'TRACK_INIT b0 = ',b0
254       print*,'TRACK_INIT bl3 = ',bl3
255       print*,'TRACK_INIT sigmacut = ',sigcut
256       print*,'TRACK_INIT cutpxz = ',pxzcut
257       print*,'TRACK_INIT xprec = ',xprec
258       print*,'TRACK_INIT yprec = ',yprec
259
260       EFF2 = EFF**2               ! PROBA. DEUX CHAMBRES TOUCHES
261       EFF1 = EFF2+2.*EFF*(1.-EFF) ! PROBA. AU MOINS UNE CHAMBRE TOUCHE
262 ** Used only for stations 4 & 5
263       PHIPREC   = SQRT(2.)*XPREC/DZ_PL(5) ! PHI = (OZ , PROJ. DANS XOZ)
264       ALAMPREC  = SQRT(2.)*YPREC/DZ_PL(5) ! LAM = (OM , PROJ. DANS XOZ)
265
266       DO I = 1,5
267          NRES(I) = 0
268       ENDDO
269       NRESF = 0 
270       NTRMUALL = 0 
271       NMUONALL = 0
272       NGHOSTALL = 0
273       NTRACKFALL = 0
274       DO I = 1,NBSTATION
275          NERRALL(I) = 0
276       ENDDO   
277       IR = 0 
278 *
279       RETURN
280       END
281
282 *************************************************************************
283       SUBROUTINE PREC_INIT
284 *************************************************************************
285 *  
286 *
287 *************************************************************************
288
289       IMPLICIT DOUBLE PRECISION (A-H, O-Z)      
290 *            
291       PARAMETER(NPLANE=10,NBSTATION=5) 
292 *           
293       COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
294 *      
295       COMMON/PARAM/ZPLANEP(NPLANE),THICK,XPREC,YPREC,B0,BL3,ZMAGS,
296      &             ZMAGE,ZABS,XMAG,ZBP1,ZBP2,CONST
297 *           
298       COMMON/PRECCUT/PCUT,PTCUT,CHI2CUT
299 *
300       COMMON/PRECSUM/NRESF1,NMUONALL1,NGHOSTALL1,NTRACKFALL1 
301 *          
302       COMMON/DEBEVT/IDEBUG
303 *
304
305       DATA THICK/0.03D0/   ! X/X0=3% 
306 **      DATA THICK/0.02D0/   ! X/X0=2% chambre
307
308 **      DATA B0,BL3/10.,2.0/ ! Champ magnetique dans le dipole et dans L3 en kgauss
309       DATA B0,BL3/7.,2.0/ ! Magnetic field in the dipole & L3 in kgauss
310       DATA ZMAGS/805.0D0/,ZMAGE/1233.0D0/,ZABS/503.D0/ ! CCC not used when
311                                                        ! magn. field map 
312 **      DATA ZMAGS/825.0D0/,ZMAGE/1125.0D0/,ZABS/503.D0/
313       DATA XMAG/190.0D0/
314
315 **      DATA XPREC/0.0100D0/,YPREC/0.144337D0/ ! CCC
316       DATA XPREC/0.0100D0/,YPREC/0.2D0/
317
318 * Input parameters       
319       CONST = 0.299792458D-3*B0*(ZMAGE-ZMAGS)
320       J = 0
321       DO I = 1,5
322          J = J+1
323          ZPLANEP(J) = -ZPLANE(I)
324          J = J+1
325          ZPLANEP(J) = -ZPLANE(I)+DZ_PL(I)  
326       ENDDO 
327
328       PCUT = 3.        ! Coupure en PXZ muon (GeV/c)
329       PTCUT = 0.       ! Coupure en Pt muon (GeV/c)
330  
331       CHI2CUT = 1.E4  ! Coupure sur le CHI2 du fit
332
333       X01 = 18.8     ! C (cm)
334       X02 = 10.397   ! Concrete (cm)
335       X03 = 0.56     ! Plomb (cm)
336       X04 = 47.26    ! Polyethylene (cm)
337       X05 = 0.35     ! W (cm) 
338
339 ** Calcul des parametres pour la correction de Branson de l'absorbeur     
340       ANBP = (315.**3-90.**3)/X01 +(467.**3-315.**3)/X02+
341      &       (472.**3-467.**3)/X03+(477.**3-472.**3)/X04+
342      &       (482.**3-477.**3)/X03+(487.**3-482.**3)/X04+
343      &       (492.**3-487.**3)/X03+(497.**3-492.**3)/X04+
344      &       (502.**3-497.**3)/X03
345       ADBP = (315.**2-90.**2)/X01 +(467.**2-315.**2)/X02+
346      &       (472.**2-467.**2)/X03+(477.**2-472.**2)/X04+
347      &       (482.**2-477.**2)/X03+(487.**2-482.**2)/X04+
348      &       (492.**2-487.**2)/X03+(497.**2-492.**2)/X04+
349      &       (502.**2-497.**2)/X03
350       ZBP1 = 2./3.*ANBP/ADBP
351       ANBP = (315.**3-90.**3)/X01 +(467.**3-315.**3)/X02+
352      &       (503.**3-467.**3)/X05
353       ADBP = (315.**2-90.**2)/X01 +(467.**2-315.**2)/X02+
354      &       (503.**2-467.**2)/X05
355       ZBP2 = 2./3.*ANBP/ADBP
356 *      
357       IF (IDEBUG.GE.1) THEN
358          PRINT *,' PREC_INIT B0 (kgauss)',B0,' BL3 (kgauss)',BL3
359          PRINT *,' PREC_INIT ZMAGE (cm)',ZMAGE,' ZMAGS (cm)',ZMAGS,
360      &        ' XMAG (cm)',XMAG
361          PRINT *,' PREC_INIT ZABS (cm)',ZABS,' ZBP1 (cm)',ZBP1,
362      &        ' ZBP2 (cm)',ZBP2
363          PRINT *,' PREC_INIT Radiation length absorber X01 (cm)',X01,
364      &        ' X02 (cm)',X02
365          PRINT *,' PREC_INIT X03(cm)',X03,' X04 (cm)',X04,' X05 (cm)',
366      &           X05
367          PRINT *,' PREC_INIT Radiation length chamber THICK (%)',
368      &           THICK*100.
369          PRINT *,' PREC_INIT XPREC (cm)',XPREC,' YPREC (cm)',YPREC
370          PRINT *,' PREC_INIT Coupure en Pxz (GeV/c): ',PCUT
371          PRINT *,' PREC_INIT Coupure en Pt (GeV/c): ',PTCUT
372          PRINT *,' PREC_INIT Coupure en CHI2 : ',CHI2CUT
373       ENDIF
374
375 *      PAUSE
376
377       NRESF1 = 0 
378       NMUONALL1 = 0
379       NGHOSTALL1 = 0
380       NTRACKFALL1 = 0
381            
382 *  Magnetic Field Map GC
383 **      OPEN (UNIT=40,FILE='data/field02.dat',
384 **     &      STATUS='UNKNOWN')
385
386       CALL INITFIELD
387                   
388 **      CLOSE(40)
389
390       RETURN
391       END
392       
393 *************************************************************************
394       SUBROUTINE RECO_READEVT(NEV,IDRES,IREADGEANT)
395 *************************************************************************
396 *
397 *
398 *************************************************************************
399       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
400
401       PARAMETER(NTRMAX=500) 
402
403       PARAMETER (NBSTATION=5,MAXIDG=20000,MAXHITTOT=20000,
404      &           MAXHITCH=10000,MAXHIT=1000,NBCHAMBER=10)
405
406       COMMON/RHITG/ITYPG(MAXIDG),XTRG(MAXIDG),YTRG(MAXIDG),
407      &             PTOTG(MAXIDG),IDG(MAXIDG),IZCH(MAXIDG),
408      &             PVERT1G(MAXIDG),PVERT2G(MAXIDG),PVERT3G(MAXIDG),
409      &     ZVERTG(MAXIDG),NHITTOT1,CX(MAXIDG),CY(MAXIDG),CZ(MAXIDG),
410      &     XGEANT(MAXIDG),YGEANT(MAXIDG),CLSIZE1(MAXIDG),CLSIZE2(MAXIDG)
411
412       DIMENSION TYPG(MAXIDG),ZCH(MAXIDG)
413
414       REAL*4 R1,R2
415       DATA R1,R2/0.,1./
416
417       IF (IREADGEANT.eq.1) THEN  ! GEANT hits
418
419          CALL TRACKF_READ_GEANT(ITYPG,XTRG,YTRG,PTOTG,IDG,IZCH,PVERT1G,
420      &        PVERT2G,PVERT3G,ZVERTG,NHITTOT1,CX,CY,CZ,IEVR,NEV,
421      &        XGEANT,YGEANT,CLSIZE1,CLSIZE2)
422       ELSE
423          CALL TRACKF_READ_SPOINT(ITYPG,XTRG,YTRG,PTOTG,IDG,IZCH,PVERT1G,
424      &        PVERT2G,PVERT3G,ZVERTG,NHITTOT1,CX,CY,CZ,IEVR,NEV,
425      &        XGEANT,YGEANT,CLSIZE1,CLSIZE2)
426       ENDIF
427       
428       do i=1,NHITTOT1
429          TYPG(i)=ITYPG(i)
430          call chfill(100,sngl(typg(i)),R1,R2)
431          call chfill(101,sngl(ygeant(i)),R1,R2)
432          call chfill(102,sngl(xgeant(i)),R1,R2)
433          ZCH(i)=IZCH(i)
434          call chfill(103,sngl(zch(i)),R1,R2)
435          call chfill(104,sngl(ptotg(i)),R1,R2)
436          call chfill(105,sngl(pvert2g(i)),R1,R2)
437          call chfill(106,sngl(pvert1g(i)),R1,R2)
438          call chfill(107,sngl(pvert3g(i)),R1,R2)
439          call chfill(108,sngl(zvertg(i)),R1,R2)
440          call chfill(109,sngl(ytrg(i)),R1,R2)
441          call chfill(110,sngl(xtrg(i)),R1,R2)
442          dx=xgeant(i)-xtrg(i)
443          dy=ygeant(i)-ytrg(i)
444          call chfill(111,sngl(dy),R1,R2)
445          call chfill(112,sngl(dx),R1,R2)
446          call chfill(119+int(zch(i)),sngl(dy),R1,R2)
447          call chfill(129+int(zch(i)),sngl(dx),R1,R2)
448       enddo
449
450       do i=1,NHITTOT1
451          CALL CHFILL (999,SNGL(PTOTG(I)),R1,R2)
452       enddo
453
454       CALL TRACKF_STAT(IDRES,IREADGEANT)
455
456       RETURN
457       END
458
459 *************************************************************************
460       SUBROUTINE TRACKF_STAT(IDRES,IREADGEANT)
461 *************************************************************************
462 *  Associate hits between two chambers inside a station
463 *  Simulate spatial resolution and chamber efficiency
464 *
465 *************************************************************************
466       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
467 *
468       PARAMETER (NBSTATION=5,MAXIDG=20000,MAXHITTOT=20000,
469      &           MAXHITCH=10000,MAXHIT=1000,NBCHAMBER=10)
470 *      
471       COMMON/TRACKFI/EFF,EFF1,EFF2,XPREC,YPREC,PHIPREC,ALAMPREC,
472      &                 HCUT,LBKG,SIGCUT,ALPHATOP,HTOP
473 *
474       COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
475 *
476 * HITS GEANT initiaux par chambre
477       COMMON/RHITG/ITYPG(MAXIDG),XTRG(MAXIDG),YTRG(MAXIDG),
478      &     PTOTG(MAXIDG),IDG(MAXIDG),IZCH(MAXIDG),
479      &     PVERT1G(MAXIDG),PVERT2G(MAXIDG),PVERT3G(MAXIDG),
480      &     ZVERTG(MAXIDG),NHITTOT1,CX(MAXIDG),CY(MAXIDG),CZ(MAXIDG),
481      &     XGEANT(MAXIDG),YGEANT(MAXIDG),CLSIZE1(MAXIDG),CLSIZE2(MAXIDG)
482
483 * HITS GEANT associes par station
484       COMMON/RHIT/ITYP(MAXHITTOT),XTR(MAXHITTOT),YTR(MAXHITTOT),
485      &            PTOT(MAXHITTOT),ID(MAXHITTOT),IZST(MAXHITTOT),
486      &            PVERT1(MAXHITTOT),PVERT2(MAXHITTOT),PVERT3(MAXHITTOT),
487      &            ZVERT(MAXHITTOT),NHITTOT
488
489
490       COMMON/CHHIT/XM(NBSTATION,MAXHITCH),YM(NBSTATION,MAXHITCH),
491      &             PHM(NBSTATION,MAXHITCH),ALM(NBSTATION,MAXHITCH),
492      &             IZM(NBSTATION,MAXHITCH),
493      &             IP(NBSTATION,MAXHITCH),JHIT(NBSTATION),
494      &             XMR(NBSTATION,MAXHITCH,2),YMR(NBSTATION,MAXHITCH,2) 
495 *
496       COMMON/DEBEVT/IDEBUG
497 *
498       DIMENSION RMIN(NBCHAMBER),RMAX1(NBCHAMBER)
499       DIMENSION XMA(NBCHAMBER,MAXHITCH),YMA(NBCHAMBER,MAXHITCH),
500      &          IMARK(NBCHAMBER,MAXHITCH)
501
502       DIMENSION IEFFI(MAXHITTOT)
503       DIMENSION IH(NBCHAMBER,MAXHIT) 
504       DIMENSION NHIT(NBCHAMBER)
505       DIMENSION DXMAX(NBSTATION),DYMAX(NBSTATION),VRES(2,5)
506   
507       REAL*4 RNDM,RN,RN1,RN2,R1,R2
508
509 * Chambre 10 deg.      
510       DATA RMAX1/91.5,91.5,122.5,122.5,158.3,158.3,260.,260.,260.,260./
511 * Zone de recherche entre deux plans d'une station 
512 **      DATA DXMAX/2.,1.5,2.5,3.,3./
513       DATA DXMAX/1.5,1.5,3.,6.,6./  ! J psi 20cm
514 **      DATA DXMAX/1.,1.,1.,1.,1./  ! Upsilon dz_ch = 8 cm
515 **      DATA DYMAX/0.22,0.22,0.21,0.21,0.21/  ! CCC Upsilon dz_ch = 8cm
516 **      DATA DXMAX/1.,1.,1.5,2.,2./  ! Upsilon dz_ch = 20 cm
517       DATA DYMAX/0.22,0.22,0.22,0.22,0.22/  ! CCC Upsilon dz_ch = 20 cm
518 **      DATA DXMAX/5.,5.,5.,5.,5./ 
519 **      DATA DXMAX/10.,10.,10.,10.,10./ 
520       DATA R1,R2/0.,1./
521
522       ICH = 0
523       DO IZ=1,5
524          ICH = ICH+1
525          RMIN(ICH) =  ABS(ZPLANE(IZ)*TAN(2.*ACOS(-1.)/180))
526          IF (IZ.GT.2) RMIN(ICH) = 30.
527          ICH = ICH+1
528          RMIN(ICH) =  ABS(ZPLANE(IZ)*TAN(2.*ACOS(-1.)/180))
529          IF (IZ.GT.2) RMIN(ICH) = 30.
530       ENDDO   
531
532 *    Initialisations 
533       DO ICH = 1,10
534          NHIT(ICH) = 0
535       ENDDO
536          
537 * 1 ere boucle Lecture des hits initiaux 
538
539       print*,'TRACKF_STAT NHITTOT1=',NHITTOT1
540       
541       IF (IREADGEANT.EQ.1) THEN
542          DO I = 1,2
543             DO J = 1,5
544                VRES(I,J) = 0.
545             ENDDO
546          ENDDO    
547          IMU = 0 
548          DO I = 1,NHITTOT1
549             ICH = IZCH(I)
550             IF (ICH.EQ.9) THEN
551                ISTAK = IDG(I)
552                ISTAK = MOD(ISTAK,30000)
553                ISTAK = MOD(ISTAK,10000)
554
555                IF ((ITYPG(I).EQ.5.OR.ITYPG(I).EQ.6).AND.
556      &              ISTAK.EQ.IDRES) THEN ! upsilon
557                   IMU = IMU+1
558                   VRES(1,IMU) = XGEANT(I)
559                   VRES(2,IMU) = YGEANT(I)
560                ENDIF
561             ENDIF   
562          ENDDO
563       ENDIF   
564
565
566       DO I = 1,NHITTOT1 ! Boucle sur les hits GEANT de toutes les ch.
567       
568 **        IF (ITYPG(I).NE.5.AND.ITYPG(I).NE.6) GOTO 1 ! CCC
569
570          ICH = IZCH(I)
571          
572          IF (IREADGEANT.EQ.1) THEN ! GEANT hits
573
574             IF (ICH.EQ.9.OR.ICH.EQ.10) THEN
575                DNUM = 999.
576                DO IM = 1,IMU
577                   DNU = SQRT((XGEANT(I)-VRES(1,IM))**2+
578      &                 (YGEANT(I)-VRES(2,IM))**2)
579                   IF (DNU.LT.DNUM) DNUM = DNU
580                ENDDO
581                IF (DNUM.GT.20.) GO TO 1 
582             ENDIF
583
584             CALL RANNOR(RN1,RN2)          ! CCCC
585             X = XGEANT(I) + RN1 * XPREC
586             Y = YGEANT(I) + RN2 * YPREC
587 *     efficacite des chambres        
588             IEFFI(I) = 1
589             RN = RNDM()              
590             IF (RN.GT.EFF) IEFFI(I) = 0
591
592
593          ELSE  ! reconstructed hits 
594
595             IEFFI(I) = 1
596
597             X = XTRG(I)
598             Y = YTRG(I)
599
600          ENDIF
601
602          R = SQRT(X**2+Y**2)
603          IF (R.LT.RMIN(ICH).OR.R.GT.RMAX1(ICH)) then
604             if (ich.le.10.and.i.le.28) then
605                print*,'* chambre ',ich,' * hit ',i
606                print*,'ityp=',itypg(i)
607                print*,'x=',X,' y=',Y
608                print*,'R=',R,' RMIN=',RMIN(ICH),' RMAX1=',RMAX1(ICH)
609            endif
610             GO TO 1   
611          end if
612
613          NHIT(ICH) = NHIT(ICH)+1
614          IH(ICH,NHIT(ICH)) = I
615          XMA(ICH,NHIT(ICH)) = X
616          YMA(ICH,NHIT(ICH)) = Y
617          IMARK(ICH,NHIT(ICH)) = 0
618
619   1      CONTINUE
620       ENDDO
621       
622 * Association des hits entre chambres d'une station
623       II = 0            ! nombre de hits GEANT par station
624       DO ICH1 = 1,10,2
625          IZ = INT(FLOAT(ICH1+1)/2.)
626          JHIT(IZ) = 0
627          ICH2 = ICH1+1
628
629          DO I1 = 1,NHIT(ICH1)
630             II = II+1
631             IFIND = 0 
632             I = IH(ICH1,I1)
633
634             ITYP(II) = ITYPG(I) 
635             XTR(II) = XTRG(I) 
636             YTR(II) = YTRG(I)
637             PTOT(II) = PTOTG(I)
638             ID(II) = IDG(I)
639             IZST(II) = IZ
640             PVERT1(II) = PVERT1G(I)  
641             PVERT2(II) = PVERT2G(I)  
642             PVERT3(II) = PVERT3G(I)
643             ZVERT(II) = ZVERTG(I)
644  
645             IF (IEFFI(I).EQ.1) THEN
646                X1 = XMA(ICH1,I1)
647                Y1 = YMA(ICH1,I1)
648                ID1 = IDG(I)
649                XEXT1 = (ZPLANE(IZ)-DZ_PL(IZ))/ZPLANE(IZ)*X1
650                YEXT1 = (ZPLANE(IZ)-DZ_PL(IZ))/ZPLANE(IZ)*Y1
651     
652                DO I2 = 1,NHIT(ICH2)
653                   J = IH(ICH2,I2)
654
655                   IF (IEFFI(J).EQ.1) THEN
656                      X2 = XMA(ICH2,I2)
657                      Y2 = YMA(ICH2,I2)
658                      ID2 = IDG(J)
659                      DX = X2-XEXT1
660                      DY = Y2-YEXT1
661
662                      IF (ID1.EQ.ID2.AND.
663      &                    (ITYP(II).EQ.5.OR.ITYP(II).EQ.6)) THEN
664                         CALL CHFILL(70+IZ,SNGL(DX),R1,R2)
665                         CALL CHFILL(80+IZ,SNGL(DY),R1,R2)
666                      ENDIF   
667                      DX = ABS(DX) 
668                      DY = ABS(DY)
669
670                      IF (DX.LT.DXMAX(IZ).AND.DY.LT.(SIGCUT*DYMAX(IZ)) ! CCC
671      &                   ) THEN
672                         IFIND = 1
673                         IMARK(ICH2,I2) = 1
674                         JHIT(IZ) =  JHIT(IZ)+1
675                         XM(IZ,JHIT(IZ)) = X1
676                         YM(IZ,JHIT(IZ)) = Y1
677                         IZM(IZ,JHIT(IZ)) = 1
678                         PHM(IZ,JHIT(IZ)) = -ATAN((X2-X1)/DZ_PL(IZ))
679                         ALM(IZ,JHIT(IZ)) = ATAN((Y2-Y1)/DZ_PL(IZ)*
680      &                                     COS(PHM(IZ,JHIT(IZ)))) 
681                         XMR(IZ,JHIT(IZ),1) = X1
682                         YMR(IZ,JHIT(IZ),1) = Y1
683                         XMR(IZ,JHIT(IZ),2) = X2
684                         YMR(IZ,JHIT(IZ),2) = Y2
685                   
686                         ISTAK = ID2
687                         ISTAK = MOD(ISTAK,30000)
688                         ISTAK = MOD(ISTAK,10000)
689
690                         IF ((ITYPG(J).EQ.5.OR.ITYPG(J).EQ.6).AND.
691      &                       ISTAK.EQ.IDRES) THEN ! upsilon
692                           
693                            ITYP(II) = ITYPG(J) 
694                            XTR(II) = XTRG(J) 
695                            YTR(II) = YTRG(J)
696                            PTOT(II) = PTOTG(J)
697                            ID(II) = IDG(J)
698                            PVERT1(II) = PVERT1G(J)  
699                            PVERT2(II) = PVERT2G(J)  
700                            PVERT3(II) = PVERT3G(J)
701                            ZVERT(II) = ZVERTG(J)
702
703                         ENDIF   
704   
705                         IP(IZ,JHIT(IZ)) = II
706
707                      ENDIF
708                   ENDIF 
709                ENDDO
710
711                IF (IFIND.EQ.0) THEN
712                   JHIT(IZ) =  JHIT(IZ)+1
713                   XM(IZ,JHIT(IZ)) = X1
714                   YM(IZ,JHIT(IZ)) = Y1
715                   IZM(IZ,JHIT(IZ)) = 1
716                   IP(IZ,JHIT(IZ)) = II
717                   PHM(IZ,JHIT(IZ)) = 10.
718                   ALM(IZ,JHIT(IZ)) = 10.
719                   XMR(IZ,JHIT(IZ),1) = X1
720                   YMR(IZ,JHIT(IZ),1) = Y1
721                   XMR(IZ,JHIT(IZ),2) = 0.
722                   YMR(IZ,JHIT(IZ),2) = 0.
723                ENDIF
724                CALL CHFILL2(1000+IZ,SNGL(X1),SNGL(Y1),R2)
725             ENDIF
726          ENDDO
727       ENDDO  
728
729 * On conserve les HITS de la 2nde chambre des stations      
730
731       DO ICH = 2,10,2 
732          IZ = INT(FLOAT(ICH+1)/2.)
733          DO I = 1,NHIT(ICH)
734             J = IH(ICH,I)
735
736             IF (IMARK(ICH,I).EQ.0) THEN
737   
738                II = II+1
739             
740                ITYP(II) = ITYPG(J) 
741                XTR(II) = XTRG(J) 
742                YTR(II) = YTRG(J)
743                PTOT(II) = PTOTG(J)
744                ID(II) = IDG(J)
745                IZST(II) = IZ
746                PVERT1(II) = PVERT1G(J)  
747                PVERT2(II) = PVERT2G(J)  
748                PVERT3(II) = PVERT3G(J)
749                ZVERT(II) = ZVERTG(I)  
750                 
751                IF (IEFFI(J).EQ.1) THEN
752                   JHIT(IZ) =  JHIT(IZ)+1
753                   XM(IZ,JHIT(IZ)) = XMA(ICH,I)
754                   YM(IZ,JHIT(IZ)) = YMA(ICH,I) 
755                   IZM(IZ,JHIT(IZ)) = 2
756                   PHM(IZ,JHIT(IZ)) = 10.
757                   ALM(IZ,JHIT(IZ)) = 10.
758                   IP(IZ,JHIT(IZ)) = II  
759                   XMR(IZ,JHIT(IZ),1) = 1000.
760                   YMR(IZ,JHIT(IZ),1) = 1000.   
761                   XMR(IZ,JHIT(IZ),2) = XMA(ICH,I)
762                   YMR(IZ,JHIT(IZ),2) = YMA(ICH,I)                    
763                ENDIF 
764             ENDIF
765          ENDDO
766       ENDDO
767
768
769       NHITTOT = II
770 *      
771
772       IF (IDEBUG.GE.2) THEN
773          PRINT *,'TRACKF_STAT nb hits:',NHITTOT 
774       ENDIF     
775
776 **      DO IZ = 1,NBSTATION
777 **         PRINT *,' IZ=',IZ,' JHIT(IZ)=',JHIT(IZ)
778 **         DO J = 1,JHIT(IZ)
779 **            II = IP(IZ,J)
780 **            print*,'II=',II
781 **           PRINT *,' ID(II)=',ID(II)
782 **         ENDDO
783 **      ENDDO 
784
785 **    DO IZ = 1,NBSTATION
786 **       PRINT *,' IZ=',IZ,' JHIT(IZ)=',JHIT(IZ)
787 **       DO J = 1,JHIT(IZ)
788 **          PRINT *,' PHM(IZ,J)=',PHM(IZ,J),' ALM(IZ,J)=',ALM(IZ,J)
789 **       ENDDO
790 **     ENDDO 
791       
792       RETURN
793       END
794 *************************************************************************
795       SUBROUTINE TRACKF_STAT_NEW(IDRES)
796 *************************************************************************
797 *  Associate hits between two chambers inside a station
798 *  Simulate spatial resolution and chamber efficiency
799 *
800 *************************************************************************
801       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
802 *
803       PARAMETER (NBSTATION=5,MAXIDG=20000,MAXHITTOT=20000,
804      &           MAXHITCH=10000,MAXHIT=1000,NBCHAMBER=10)
805 *      
806       COMMON/TRACKFI/EFF,EFF1,EFF2,XPREC,YPREC,PHIPREC,ALAMPREC,
807      &                 HCUT,LBKG,SIGCUT,ALPHATOP,HTOP
808 *
809       COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
810 *
811 * HITS GEANT initiaux par chambre
812       COMMON/RHITG/ITYPG(MAXIDG),XTRG(MAXIDG),YTRG(MAXIDG),
813      &     PTOTG(MAXIDG),IDG(MAXIDG),IZCH(MAXIDG),
814      &     PVERT1G(MAXIDG),PVERT2G(MAXIDG),PVERT3G(MAXIDG),
815      &     ZVERTG(MAXIDG),NHITTOT1,CX(MAXIDG),CY(MAXIDG),CZ(MAXIDG),
816      &     XGEANT(MAXIDG),YGEANT(MAXIDG),CLSIZE1(MAXIDG),CLSIZE2(MAXIDG)
817
818 * HITS GEANT associes par station
819       COMMON/RHIT/ITYP(MAXHITTOT),XTR(MAXHITTOT),YTR(MAXHITTOT),
820      &            PTOT(MAXHITTOT),ID(MAXHITTOT),IZST(MAXHITTOT),
821      &            PVERT1(MAXHITTOT),PVERT2(MAXHITTOT),PVERT3(MAXHITTOT),
822      &            ZVERT(MAXHITTOT),NHITTOT
823
824
825       COMMON/CHHIT/XM(NBSTATION,MAXHITCH),YM(NBSTATION,MAXHITCH),
826      &             PHM(NBSTATION,MAXHITCH),ALM(NBSTATION,MAXHITCH),
827      &             IZM(NBSTATION,MAXHITCH),
828      &             IP(NBSTATION,MAXHITCH),JHIT(NBSTATION),
829      &             XMR(NBSTATION,MAXHITCH,2),YMR(NBSTATION,MAXHITCH,2) 
830 *
831       COMMON/DEBEVT/IDEBUG
832 *
833       DIMENSION RMIN(NBCHAMBER),RMAX1(NBCHAMBER)
834       DIMENSION XMA(NBCHAMBER,MAXHITCH),YMA(NBCHAMBER,MAXHITCH),
835      &          IMARK(NBCHAMBER,MAXHITCH)
836
837       DIMENSION IEFFI(MAXHITTOT)
838       DIMENSION IH(NBCHAMBER,MAXHIT) 
839       DIMENSION NHIT(NBCHAMBER)
840       DIMENSION DXMAX(NBSTATION),DYMAX(NBSTATION),I2C(1000)
841   
842       DIMENSION DIST(2),NMUON(2),NHITMUON(2,5),NMUONGOOD(2)
843
844       REAL*4 RNDM,RN,RN1,RN2,R1,R2
845
846 * Chambre 10 deg.      
847       DATA RMAX1/91.5,91.5,122.5,122.5,158.3,158.3,260.,260.,260.,260./
848 * Zone de recherche entre deux plans d'une station 
849 **      DATA DXMAX/2.,1.5,2.5,3.,3./
850       DATA DXMAX/1.5,1.5,3.,3.,3./
851       DATA DYMAX/0.22,0.22,0.21,0.21,0.21/
852       DATA R1,R2/0.,1./
853
854       ICH = 0
855       DO IZ=1,5
856          ICH = ICH+1
857          RMIN(ICH) =  ABS(ZPLANE(IZ)*TAN(2.*ACOS(-1.)/180))
858          IF (IZ.GT.2) RMIN(ICH) = 30.
859          ICH = ICH+1
860          RMIN(ICH) =  ABS(ZPLANE(IZ)*TAN(2.*ACOS(-1.)/180))
861          IF (IZ.GT.2) RMIN(ICH) = 30.
862       ENDDO   
863
864 *    Initialisations 
865       DO ICH = 1,10
866          NHIT(ICH) = 0
867       ENDDO
868          
869       DO NCH = 1,10
870          DO J=1,2
871             DIST(J)=999.
872             NMUON(J)=0
873          ENDDO
874          DO I = 1,NHITTOT1
875             IF (IZCH(I).EQ.NCH) THEN
876                ISTAK = IDG(I)
877                ISTAK = MOD(ISTAK,30000)
878                ISTAK = MOD(ISTAK,10000)
879                IF (ISTAK.EQ.IDRES.AND.IDG(I).EQ.50116) THEN
880                   DISTMIN=(XTRG(I)-XGEANT(I))**2+(YTRG(I)-YGEANT(I))**2
881                   IF (DISTMIN.LT.DIST(1)) THEN
882                      DIST(1)=DISTMIN
883                      NMUONGOOD(1)=I
884                   ENDIF
885                   NMUON(1)=NMUON(1)+1
886                   NHITMUON(1,NMUON(1))=I
887                ELSE 
888                 IF (ISTAK.EQ.IDRES.AND.IDG(I).EQ.70116) THEN
889                    DISTMIN=(XTRG(I)-XGEANT(I))**2+(YTRG(I)-YGEANT(I))**2
890                    IF (DISTMIN.LT.DIST(2)) THEN
891                         DIST(2)=DISTMIN
892                         NMUONGOOD(2)=I
893                    ENDIF
894                    NMUON(2)=NMUON(2)+1
895                    NHITMUON(2,NMUON(2))=I
896                   ENDIF
897                ENDIF
898             ENDIF
899          ENDDO
900          DO J=1,2
901             IF (NMUON(J).GE.2) THEN
902              print*,'j=',j,' nmuon=',nmuon(j)
903              print*,'chambre',nch
904              DO K=1,NMUON(J)
905                IF (NHITMUON(J,K).NE.NMUONGOOD(J)) IDG(NHITMUON(J,K))=999 ! flag les mauvais hits MUONS
906              ENDDO
907             ENDIF
908          ENDDO
909       ENDDO
910       
911
912 * 1 ere boucle Lecture des hits initiaux 
913
914       DO I = 1,NHITTOT1 ! Boucle sur les hits GEANT de toutes les ch.
915       
916          ICH = IZCH(I)
917          
918          X = XTRG(I)
919          Y = YTRG(I)
920
921          R = SQRT(X**2+Y**2)
922          IF (R.LT.RMIN(ICH).OR.R.GT.RMAX1(ICH)) then
923             if (ich.le.10.and.i.le.28) then
924                print*,'****** chambre ',ich,' ****** hit ',i
925                print*,'ityp=',itypg(i)
926                print*,'x=',XTRG(I),' y=',YTRG(I)
927                print*,'R=',R,' RMIN=',RMIN(ICH),' RMAX1=',RMAX1(ICH)
928            endif
929             GO TO 1   
930          end if
931
932          IEFFI(I) = 1
933
934          NHIT(ICH) = NHIT(ICH)+1
935          IH(ICH,NHIT(ICH)) = I
936          XMA(ICH,NHIT(ICH)) = XTRG(I)
937          YMA(ICH,NHIT(ICH)) = YTRG(I)
938          IMARK(ICH,NHIT(ICH)) = 0
939
940          print*,' XTRG(I)=', XTRG(I),' YTRG(I)=', YTRG(I),' IDG(I)=',
941      &        IDG(I),' ICH=',ICH
942          
943   1      CONTINUE
944       ENDDO
945       
946       
947 * Association des hits entre chambres d'une station
948       II = 0            ! nombre de hits GEANT par station
949       DO ICH1 = 1,10,2
950          IZ = INT(FLOAT(ICH1+1)/2.)
951          JHIT(IZ) = 0
952          ICH2 = ICH1+1
953
954          DO I1 = 1,NHIT(ICH1)
955             II = II+1
956             IFIND = 0 
957             I = IH(ICH1,I1)
958
959             ITYP(II) = ITYPG(I) 
960             XTR(II) = XTRG(I) 
961             YTR(II) = YTRG(I)
962             PTOT(II) = PTOTG(I)
963             ID(II) = IDG(I)
964             IZST(II) = IZ
965             PVERT1(II) = PVERT1G(I)  
966             PVERT2(II) = PVERT2G(I)  
967             PVERT3(II) = PVERT3G(I)
968             ZVERT(II) = ZVERTG(I)
969  
970             IF (IEFFI(I).EQ.1) THEN
971                X1 = XMA(ICH1,I1)
972                Y1 = YMA(ICH1,I1)
973                ID1 = IDG(I)
974                XEXT1 = (ZPLANE(IZ)-DZ_PL(IZ))/ZPLANE(IZ)*X1
975                YEXT1 = (ZPLANE(IZ)-DZ_PL(IZ))/ZPLANE(IZ)*Y1    
976                KC = 0
977                PRINT *,'***** DEBUT RECHERCHE',' ID1=',ID1,' ich1=',ICH1
978                PRINT *,' XTR(II)=', XTR(II),' YTR(II)=', YTR(II) 
979                PRINT *,'  ITYP(II)=', ITYP(II)
980                DO I2 = 1,NHIT(ICH2)
981                   J = IH(ICH2,I2)
982                   IF (IEFFI(J).EQ.1) THEN
983                      X2 = XMA(ICH2,I2)
984                      DX = X2-XEXT1
985                      DX = ABS(DX) 
986                      IF (DX.LT.DXMAX(IZ)) THEN
987                         KC = KC + 1
988                         I2C(KC) = I2
989                         ID2 = IDG(J)
990                         print *,' DX=',DX,' KC=',KC,' ID2=',ID2 
991                      ENDIF
992                   ENDIF
993                ENDDO   
994                DYOLD = 999.
995                I2FIND = 0
996                DO IKC = 1,KC
997                   I2 = I2C(IKC)
998                   Y2 = YMA(ICH2,I2)
999                   DY = Y2-YEXT1
1000                   DY = ABS(DY)
1001                   J = IH(ICH2,I2)
1002                   ID2 = IDG(J)
1003                   IF (DY.LT.DYOLD.AND.DY.LT.(SIGCUT*DYMAX(IZ))) THEN
1004                      DYOLD = DY
1005                      I2FIND = I2
1006                      PRINT *,' ID2=',ID2,' DY=',DY
1007                   ENDIF
1008                ENDDO   
1009                IF (I2FIND.GT.0) THEN
1010                   I2 = I2FIND
1011                   J = IH(ICH2,I2)
1012                   ID2 = IDG(J)
1013                   IFIND = 1
1014                   IMARK(ICH2,I2) = 1
1015                   JHIT(IZ) =  JHIT(IZ)+1
1016                   X2 = XMA(ICH2,I2)
1017                   Y2 = YMA(ICH2,I2)
1018                   XM(IZ,JHIT(IZ)) = X1
1019                   YM(IZ,JHIT(IZ)) = Y1
1020                   IZM(IZ,JHIT(IZ)) = 1
1021                   PHM(IZ,JHIT(IZ)) = -ATAN((X2-X1)/DZ_PL(IZ))
1022                   ALM(IZ,JHIT(IZ)) = ATAN((Y2-Y1)/DZ_PL(IZ)*
1023      &                 COS(PHM(IZ,JHIT(IZ)))) 
1024                   XMR(IZ,JHIT(IZ),1) = X1
1025                   YMR(IZ,JHIT(IZ),1) = Y1
1026                   XMR(IZ,JHIT(IZ),2) = X2
1027                   YMR(IZ,JHIT(IZ),2) = Y2
1028                   IP(IZ,JHIT(IZ)) = II
1029                   
1030                   ISTAK = ID2
1031                   ISTAK = MOD(ISTAK,30000)
1032                   ISTAK = MOD(ISTAK,10000)
1033 * test
1034                   IF (ISTAK.EQ.IDRES.AND.ID1.NE.999) THEN
1035                           
1036                      ITYP(II) = ITYPG(J) 
1037                      PTOT(II) = PTOTG(J)
1038                      XTR(II) = XTRG(I) 
1039                      YTR(II) = YTRG(I)
1040                      ID(II) = IDG(J)
1041                      PVERT1(II) = PVERT1G(J)  
1042                      PVERT2(II) = PVERT2G(J)  
1043                      PVERT3(II) = PVERT3G(J)
1044                      ZVERT(II) = ZVERTG(J)
1045
1046                   ENDIF  
1047                ENDIF   
1048   
1049                            
1050           
1051                IF (IFIND.EQ.0) THEN
1052                   JHIT(IZ) =  JHIT(IZ)+1
1053                   XM(IZ,JHIT(IZ)) = X1
1054                   YM(IZ,JHIT(IZ)) = Y1
1055                   IZM(IZ,JHIT(IZ)) = 1
1056                   IP(IZ,JHIT(IZ)) = II
1057                   PHM(IZ,JHIT(IZ)) = 10.
1058                   ALM(IZ,JHIT(IZ)) = 10.
1059                   XMR(IZ,JHIT(IZ),1) = X1
1060                   YMR(IZ,JHIT(IZ),1) = Y1
1061                   XMR(IZ,JHIT(IZ),2) = 0.
1062                   YMR(IZ,JHIT(IZ),2) = 0.
1063                ENDIF
1064                CALL CHFILL2(1000+IZ,SNGL(X1),SNGL(Y1),R2)
1065             ENDIF
1066          ENDDO
1067       ENDDO  
1068
1069 * On conserve les HITS de la 2nde chambre des stations      
1070
1071       DO ICH = 2,10,2 
1072          IZ = INT(FLOAT(ICH+1)/2.)
1073          DO I = 1,NHIT(ICH)
1074             J = IH(ICH,I)
1075
1076             IF (IMARK(ICH,I).EQ.0) THEN
1077 **               print *,' ich=',ich,' i=',i,' j=',j
1078   
1079                II = II+1
1080             
1081                ITYP(II) = ITYPG(J) 
1082                XTR(II) = XTRG(J) 
1083                YTR(II) = YTRG(J)
1084                PTOT(II) = PTOTG(J)
1085                ID(II) = IDG(J)
1086                IZST(II) = IZ
1087                PVERT1(II) = PVERT1G(J)  
1088                PVERT2(II) = PVERT2G(J)  
1089                PVERT3(II) = PVERT3G(J)
1090                ZVERT(II) = ZVERTG(I)  
1091                 
1092                IF (IEFFI(J).EQ.1) THEN
1093                   JHIT(IZ) =  JHIT(IZ)+1
1094 **                  XM(IZ,JHIT(IZ)) = ZPLANE(IZ)/(ZPLANE(IZ)-DZ_PL)
1095 **     &                              *XMA(ICH,I)
1096 **                  YM(IZ,JHIT(IZ)) = ZPLANE(IZ)/(ZPLANE(IZ)-DZ_PL)
1097 **     &                              *YMA(ICH,I) 
1098                   XM(IZ,JHIT(IZ)) = XMA(ICH,I)
1099                   YM(IZ,JHIT(IZ)) = YMA(ICH,I) 
1100                   IZM(IZ,JHIT(IZ)) = 2
1101                   PHM(IZ,JHIT(IZ)) = 10.
1102                   ALM(IZ,JHIT(IZ)) = 10.
1103                   IP(IZ,JHIT(IZ)) = II  
1104                   XMR(IZ,JHIT(IZ),1) = 1000.
1105                   YMR(IZ,JHIT(IZ),1) = 1000.   
1106                   XMR(IZ,JHIT(IZ),2) = XMA(ICH,I)
1107                   YMR(IZ,JHIT(IZ),2) = YMA(ICH,I)                    
1108                ENDIF 
1109             ENDIF
1110          ENDDO
1111       ENDDO
1112
1113
1114       NHITTOT = II
1115 *      
1116       IF (IDEBUG.GE.2) THEN
1117          PRINT *,'TRACKF_MICRO nb hits:',NHITTOT 
1118       ENDIF     
1119       DO IZ = 1,NBSTATION
1120          PRINT *,' IZ=',IZ,' JHIT(IZ)=',JHIT(IZ)
1121          DO J = 1,JHIT(IZ)
1122             II = IP(IZ,J)
1123             PRINT *,' ID(II)=',ID(II)
1124             PRINT *,' XMR(IZ,J,1)=', XMR(IZ,J,1),
1125      &            ' YMR(IZ,J,1)=', YMR(IZ,J,1)   
1126             PRINT *,' XMR(IZ,J,2)=', XMR(IZ,J,2),
1127      &            ' YMR(IZ,J,2)=', YMR(IZ,J,2)   
1128         ENDDO
1129       ENDDO 
1130       
1131       RETURN
1132       END
1133
1134 *************************************************************************
1135       SUBROUTINE RECO_TRACKF(IDRES,IREADGEANT)
1136 *************************************************************************
1137 *
1138 *
1139 *************************************************************************
1140       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
1141
1142       PARAMETER (MAXIDG=20000)
1143
1144       PARAMETER (MAXHITCH=10000,MAXTRK=50000,MAXHITTOT=20000,
1145      &           NBSTATION=5,MAXCAN=1000,NBCHAMBER=10,NTRMAX=500)
1146 *      
1147       COMMON/CHHIT/XM(NBSTATION,MAXHITCH),YM(NBSTATION,MAXHITCH),
1148      &             PHM(NBSTATION,MAXHITCH),ALM(NBSTATION,MAXHITCH),
1149      &             IZM(NBSTATION,MAXHITCH),
1150      &             IP(NBSTATION,MAXHITCH),JHIT(NBSTATION),
1151      &             XMR(NBSTATION,MAXHITCH,2),YMR(NBSTATION,MAXHITCH,2) 
1152 *      
1153       COMMON/RHIT/ITYP(MAXHITTOT),XTR(MAXHITTOT),YTR(MAXHITTOT),
1154      &            PTOT(MAXHITTOT),ID(MAXHITTOT),IZST(MAXHITTOT),
1155      &            PVERT1(MAXHITTOT),PVERT2(MAXHITTOT),PVERT3(MAXHITTOT),
1156      &            ZVERT(MAXHITTOT),NHITTOT
1157
1158
1159       COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
1160 *      
1161       COMMON/VERIFGEANT/ITTROUGH(MAXTRK,NBSTATION),
1162      &                  IT_LIST(MAXTRK),IT_NP(MAXTRK),ITCHECK(MAXTRK),
1163      &                  ITRACK(MAXHITTOT)
1164 *     
1165       COMMON/HCHHIT/HHIT(MAXHITCH),INDEXTAB(MAXHITCH),INDEXMAX
1166
1167       COMMON/PRECIS/EEXM(NBSTATION),EEYM(NBSTATION),EEPH(NBSTATION),
1168      &     EEAL(NBSTATION)
1169 *
1170       COMMON/MEASUR/XMES(NBSTATION),YMES(NBSTATION),IZMES(NBSTATION),
1171      &              PHMES(NBSTATION),ALMES(NBSTATION),MPOS(NBSTATION),
1172      &              MANG(NBSTATION)
1173 *
1174       COMMON/PLANE/XPL(NBSTATION,2),YPL(NBSTATION,2),
1175      &             PHPL(NBSTATION),ALPL(NBSTATION),CHI2PL
1176 *
1177       COMMON/CANDIDAT/JCAN(NBSTATION,MAXCAN),JCANTYP(MAXCAN),
1178      &          EEX(MAXCAN),EEY(MAXCAN),EEP(MAXCAN),EEA(MAXCAN)
1179 *     
1180       COMMON /VERTEX/ERRV,IVERTEX
1181 *
1182       COMMON/TRACKSUM/NRES(5),NRESF,NTRMUALL,NMUONALL,NGHOSTALL,
1183      &     NTRACKFALL,NERRALL(NBSTATION),IR              
1184 *            
1185       COMMON/TRACKFOUT/IEVOUT,NTREVT,JJOUT(NBCHAMBER,NTRMAX),
1186      &               ISTAT(NTRMAX),PXZOUT(NTRMAX),TPHIOUT(NTRMAX),
1187      &               TALAMOUT(NTRMAX),XVERTOUT(NTRMAX),YVERTOUT(NTRMAX),
1188      &               CHI2OUT(NTRMAX),
1189      &               XMESOUT(NBCHAMBER,NTRMAX),YMESOUT(NBCHAMBER,NTRMAX)  
1190      &              ,PXVOUT(NTRMAX),PYVOUT(NTRMAX),PZVOUT(NTRMAX)
1191 *
1192       COMMON/TRACKFI/EFF,EFF1,EFF2,XPREC,YPREC,PHIPREC,ALAMPREC,
1193      &                 HCUT,LBKG,SIGCUT,ALPHATOP,HTOP
1194 *
1195       COMMON/DEBEVT/IDEBUG
1196
1197       DIMENSION PEST(5),PSTEP(5),NERR(NBSTATION)
1198 *
1199       REAL*4 R2
1200       DATA R2/1./
1201
1202 ** GEANT informations   
1203       DO I = 1,MAXTRK
1204          IT_LIST(I)= 0
1205       ENDDO
1206       DO I = 1,MAXHITTOT
1207          ITRACK(I) = 0
1208       ENDDO
1209       NTRTOT = 0
1210 *  BOUCLE SUR LES HITS
1211       
1212       DO IH = 1,NHITTOT
1213          DO IT = 1,NTRTOT            
1214             IF (ID(IH).EQ.IT_LIST(IT))  GOTO 4
1215          ENDDO
1216          NTRTOT = NTRTOT +1     ! NB de trace GEANT
1217          IT_LIST(NTRTOT) = ID(IH) ! ID DE LA TRACE NTRTOT
1218          IT_NP(NTRTOT) = 0      ! NOMBRE DE HITS PAR TRACE
1219          DO II=1,NBSTATION
1220             ITTROUGH(NTRTOT,II) = 0
1221          ENDDO
1222          IT = NTRTOT
1223  4       IT_NP(IT) = IT_NP(IT)+1
1224          ITTROUGH(IT,IZST(IH))=IH ! =IH si la trace IT touche le plan IZST
1225          ITRACK(IH) = IT
1226       ENDDO
1227       IF (IDEBUG.GE.2) THEN
1228          PRINT *,'RECO_TRACKF nb total de trace GEANT:',NTRTOT
1229       ENDIF 
1230       NTRPART=0
1231       NTRMU = 0
1232       
1233       DO IT = 1,NTRTOT
1234 *     print*,'(ITTROUGH(IT,1)=',ITTROUGH(IT,1)
1235 *     print*,'(ITTROUGH(IT,2)=',ITTROUGH(IT,2)
1236 *     print*,'(ITTROUGH(IT,3)=',ITTROUGH(IT,3)
1237 *     print*,'(ITTROUGH(IT,4)=',ITTROUGH(IT,4)
1238 *     print*,'(ITTROUGH(IT,5)=',ITTROUGH(IT,5)
1239          ITCHECK(IT) = 0
1240          IF ((ITTROUGH(IT,1)*ITTROUGH(IT,2)*ITTROUGH(IT,3)*
1241      &        ITTROUGH(IT,4)*ITTROUGH(IT,5)).NE.0)
1242      &        THEN
1243             NTRPART=NTRPART+1
1244             IH = ITTROUGH(IT,NBSTATION)
1245             IF (ITYP(IH).EQ.5.OR.ITYP(IH).EQ.6) THEN
1246                ISTAK = ID(IH)
1247                ISTAK = MOD(ISTAK,30000)
1248                ISTAK = MOD(ISTAK,10000)
1249 *     test
1250                pt=sqrt(pvert1(ih)**2+pvert2(ih)**2)
1251                thet=datan2(pt,pvert3(ih))*180/3.1416
1252                pp=sqrt(pt**2+pvert3(ih)**2)
1253                
1254 *     print*,'istak=',istak
1255                
1256                IF (ISTAK.EQ.IDRES) THEN ! psi or upsilon
1257                   print*,'resonance'
1258                   NTRMU = NTRMU+1
1259                   NTRMUALL = NTRMUALL+1
1260                   ITCHECK(IT) = 1
1261                ENDIF
1262             ENDIF
1263          ENDIF
1264       ENDDO
1265
1266       IF (IDEBUG.GE.2) THEN
1267          PRINT *,'RECO_TRACKF nb of part. GEANT crossing 5 st.:',
1268      &        NTRPART
1269          PRINT *,'RECO_TRACKF nb of muons/res. GEANT crossing 5 st.:',
1270      &        NTRMU
1271       ENDIF
1272       
1273 **         CALL H_ACCEPTANCE(5)
1274 **         CALL H_ACCEPTANCE(4)
1275 **         PAUSE
1276          
1277       NCAN = 0 
1278
1279 *  Recherche 5 -> 4
1280
1281       CALL ORDONNE_HIT(5,HCUT)
1282
1283       DO IH = 1,INDEXMAX
1284          JJ = INDEXTAB(IH)
1285          X1 = XM(5,JJ)-(ZPLANE(5)-ZPLANE(4))*TAN(PHM(5,JJ))
1286          Y1 = YM(5,JJ)+(ZPLANE(5)-ZPLANE(4))*TAN(ALM(5,JJ))
1287      &                                            /COS(PHM(5,JJ))
1288          X2 = XM(5,JJ)-(ZPLANE(5)-ZPLANE(4)+DZ_PL(4))*TAN(PHM(5,JJ))
1289          Y2 = YM(5,JJ)+(ZPLANE(5)-ZPLANE(4)+DZ_PL(4))*TAN(ALM(5,JJ))
1290      &                                            /COS(PHM(5,JJ))
1291 *  Domaine de recherche dans la st. 4
1292          HCONST = 0.0136*SQRT(0.06)/HTOP  ! -> X/X0=0.03 % / chamber
1293          EPH2 = 2.0*PHIPREC**2  + (HCONST*HHIT(JJ))**2  
1294          EAL2 = 2.0*ALAMPREC**2 + (HCONST*HHIT(JJ))**2
1295 ** 29.07
1296          EXM2 = 2.0*(XPREC/SQRT(2.))**2    
1297      &          + (ZPLANE(5)-ZPLANE(4))**2/2.*EPH2
1298          EYM2 = 2.0*YPREC**2    + (ZPLANE(5)-ZPLANE(4))**2/2.*EAL2
1299          P1 = PTOT(IP(5,JJ))
1300          CALL CHFILL2(2500,SNGL(P1),SNGL(HHIT(JJ)),1.)
1301          CALL CHFILL2(2501,SNGL(P1),SNGL(HHIT(JJ)**2),1.)
1302          CALL CHFILL2(2502,SNGL(P1),SNGL(EPH2),1.)
1303          CALL CHFILL2(2503,SNGL(P1),SNGL(EAL2),1.)
1304          CALL CHFILL2(2504,SNGL(P1),SNGL(EXM2),1.)
1305          CALL CHFILL2(2505,SNGL(P1),SNGL(EYM2),1.)
1306
1307          CALL CHFILL2(2507,SNGL(P1),SNGL(SQRT(EPH2)),1.)
1308          CALL CHFILL2(2508,SNGL(P1),SNGL(SQRT(EAL2)),1.)
1309          CALL CHFILL2(2509,SNGL(P1),SNGL(SQRT(EXM2)),1.)
1310          CALL CHFILL2(2510,SNGL(P1),SNGL(SQRT(EYM2)),1.)
1311 ** end 29.07
1312          EXM12 = (PHIPREC**2+(HCONST*HHIT(JJ))**2)*
1313      &        (ZPLANE(5)-ZPLANE(4))**2 + XPREC**2
1314          EYM12 = (ALAMPREC**2+(HCONST*HHIT(JJ))**2)*
1315      &        (ZPLANE(5)-ZPLANE(4))**2 + YPREC**2  
1316          EXM22 = (PHIPREC**2+(HCONST*HHIT(JJ))**2)*
1317      &        (ZPLANE(5)-ZPLANE(4)+DZ_PL(4))**2 + XPREC**2
1318          EYM22 = (ALAMPREC**2+(HCONST*HHIT(JJ))**2)*
1319      &        (ZPLANE(5)-ZPLANE(4)+DZ_PL(4))**2 + YPREC**2  
1320  
1321          EPH = SIGCUT*SQRT(EPH2)
1322          EAL = SIGCUT*SQRT(EAL2)
1323          EX1 = SIGCUT*SQRT(EXM12)
1324          EY1 = SIGCUT*SQRT(EYM12)
1325          EX2 = SIGCUT*SQRT(EXM22)
1326          EY2 = SIGCUT*SQRT(EYM22)
1327 ** 29.07
1328          EXM = SIGCUT*SQRT(EXM2)
1329          EYM = SIGCUT*SQRT(EYM2)
1330          CALL CHFILL2(2511,SNGL(P1),SNGL(EPH),1.)
1331          CALL CHFILL2(2512,SNGL(P1),SNGL(EAL),1.)
1332          CALL CHFILL2(2513,SNGL(P1),SNGL(EXM),1.)
1333          CALL CHFILL2(2514,SNGL(P1),SNGL(EYM),1.)
1334 ** end 29.07        
1335 **         P2 = (HTOP/HHIT(JJ))**2
1336
1337 **         EPH=SIGCUT*SQRT(9.14D-7+1.2D-3/P2)
1338 **         EAL=SIGCUT*SQRT(1.84D-4)
1339 **         EX1=SIGCUT*SQRT(1.95D-2+6.37/P2)
1340 **         EY1=SIGCUT*SQRT(3.89+151./P2)
1341 **         EX2=EX1
1342 **         EY2=EY2
1343  
1344 * renvoie le num de hit de 4 le plus pres dans le domaine de recherche
1345
1346          CALL DISTMIN4(X1,Y1,PHM(5,JJ),ALM(5,JJ),4,EX1,EY1,EPH,EAL,
1347      &        IFIND)
1348          P1 = PTOT(IP(5,JJ))
1349          CALL CHECK_HISTO4(11,4,IFIND,5,JJ,X1,Y1,PHM(5,JJ),ALM(5,JJ),P1,
1350      &                     EX1,EY1,EPH,EAL)   
1351          IF (IFIND.GT.0) THEN 
1352             CALL STOCK_CANDIDAT(5,JJ,4,IFIND,EX1,EY1,EPH,EAL,NCAN,1)
1353          ELSE
1354             CALL DISTMIN2(X1,Y1,X2,Y2,4,EX1,EY1,EX2,EY2,IFIND)
1355             CALL CHECK_HISTO2(0,4,IFIND,5,JJ,X1,Y1,X2,Y2,P1,EX1,EY1,
1356      &           EX2,EY2)   
1357             IF (IFIND.GT.0) THEN
1358                CALL STOCK_CANDIDAT(5,JJ,4,IFIND,EX1,EY1,EPH,EAL,NCAN,2)
1359             ENDIF 
1360          ENDIF   
1361       ENDDO
1362          
1363 *   Recherche 4 -> 5
1364
1365       CALL ORDONNE_HIT(4,HCUT)
1366
1367       DO IH = 1,INDEXMAX
1368          JJ = INDEXTAB(IH)
1369          X1 = XM(4,JJ)-(ZPLANE(4)-ZPLANE(5))*TAN(PHM(4,JJ))
1370          Y1 = YM(4,JJ)+(ZPLANE(4)-ZPLANE(5))*TAN(ALM(4,JJ))
1371      &                                            /COS(PHM(4,JJ))
1372          X2 = XM(4,JJ)-(ZPLANE(4)-ZPLANE(5)+DZ_PL(5))*TAN(PHM(4,JJ))
1373          Y2 = YM(4,JJ)+(ZPLANE(4)-ZPLANE(5)+DZ_PL(5))*TAN(ALM(4,JJ))
1374      &                                            /COS(PHM(4,JJ))
1375 *   Domaine de recherche dans la st. 5
1376          HCONST = 0.0136*SQRT(0.06)/HTOP  ! -> X/X0=0.03 / chamber
1377 **         EPH2 = 2.0*PHIPREC**2  + (HCONST*HHIT(JJ))**2   
1378 **         EAL2 = 2.0*ALAMPREC**2 + (HCONST*HHIT(JJ))**2
1379 **         EX12 = 2.0*(XPREC/SQRT(2.))**2    
1380 **     &          + (ZPLANE(5)-ZPLANE(4))**2/2.*EPH2
1381 **         EY12 = 2.0*YPREC**2 + (ZPLANE(5)-ZPLANE(4))**2/2.*EAL2
1382 **         EX22 = 2.0*(XPREC/SQRT(2.))**2    
1383 **     &          + (ZPLANE(5)-ZPLANE(4)-DZ_PL(5))**2/2.*EPH2
1384 **         EY22 = 2.0*YPREC**2 + (ZPLANE(5)-ZPLANE(4)-DZ_PL(5))**2/2.
1385 **     &        *EAL2
1386          EPH2 = 2.0*PHIPREC**2  + (HCONST*HHIT(JJ))**2  
1387          EAL2 = 2.0*ALAMPREC**2 + (HCONST*HHIT(JJ))**2
1388 ** 29.07
1389          EXM2 = 2.0*(XPREC/SQRT(2.))**2    
1390      &          + (ZPLANE(5)-ZPLANE(4))**2/2.*EPH2
1391          EYM2 = 2.0*YPREC**2    + (ZPLANE(5)-ZPLANE(4))**2/2.*EAL2
1392          P1 = PTOT(IP(4,JJ))
1393          CALL CHFILL2(2400,SNGL(P1),SNGL(HHIT(JJ)),1.)
1394          CALL CHFILL2(2401,SNGL(P1),SNGL(HHIT(JJ)**2),1.)
1395          CALL CHFILL2(2402,SNGL(P1),SNGL(EPH2),1.)
1396          CALL CHFILL2(2403,SNGL(P1),SNGL(EAL2),1.)
1397          CALL CHFILL2(2404,SNGL(P1),SNGL(EXM2),1.)
1398          CALL CHFILL2(2405,SNGL(P1),SNGL(EYM2),1.)
1399
1400          CALL CHFILL2(2407,SNGL(P1),SNGL(SQRT(EPH2)),1.)
1401          CALL CHFILL2(2408,SNGL(P1),SNGL(SQRT(EAL2)),1.)
1402          CALL CHFILL2(2409,SNGL(P1),SNGL(SQRT(EXM2)),1.)
1403          CALL CHFILL2(2410,SNGL(P1),SNGL(SQRT(EYM2)),1.)
1404 ** end 29.07
1405
1406          EXM12 = (PHIPREC**2+(HCONST*HHIT(JJ))**2)*
1407      &        (ZPLANE(5)-ZPLANE(4))**2 + XPREC**2
1408          EYM12 = (ALAMPREC**2+(HCONST*HHIT(JJ))**2)*
1409      &        (ZPLANE(5)-ZPLANE(4))**2 + YPREC**2  
1410          EXM22 = (PHIPREC**2+(HCONST*HHIT(JJ))**2)*
1411      &        (ZPLANE(5)-ZPLANE(4)-DZ_PL(5))**2 + XPREC**2
1412          EYM22 = (ALAMPREC**2+(HCONST*HHIT(JJ))**2)*
1413      &        (ZPLANE(5)-ZPLANE(4)-DZ_PL(5))**2 + YPREC**2  
1414
1415          EPH = SIGCUT*SQRT(EPH2)
1416          EAL = SIGCUT*SQRT(EAL2)
1417          EX1 = SIGCUT*SQRT(EXM12)
1418          EY1 = SIGCUT*SQRT(EYM12)
1419          EX2 = SIGCUT*SQRT(EXM22)
1420          EY2 = SIGCUT*SQRT(EYM22)
1421 ** 29.07
1422          EXM = SIGCUT*SQRT(EXM2)
1423          EYM = SIGCUT*SQRT(EYM2)
1424          CALL CHFILL2(2411,SNGL(P1),SNGL(EPH),1.)
1425          CALL CHFILL2(2412,SNGL(P1),SNGL(EAL),1.)
1426          CALL CHFILL2(2413,SNGL(P1),SNGL(EXM),1.)
1427          CALL CHFILL2(2414,SNGL(P1),SNGL(EYM),1.)
1428 ** end 29.07
1429
1430 **         P2 = (HTOP/HHIT(JJ))**2
1431
1432 **         EPH=SIGCUT*SQRT(9.14D-7+1.2D-3/P2)
1433 **         EAL=SIGCUT*SQRT(1.84D-4)
1434 **         EX1=SIGCUT*SQRT(1.95D-2+6.37/P2)
1435 **         EY1=SIGCUT*SQRT(3.89+151./P2)
1436 **         EX2=EX1
1437 **         EY2=EY2
1438 *   Renvoie le num de hit de 5 le plus pres dans le domaine de recherche
1439          CALL DISTMIN2(X1,Y1,X2,Y2,5,EX1,EY1,EX2,EY2,IFIND)
1440          P1 = PTOT(IP(4,JJ))
1441          CALL CHECK_HISTO2(0,5,IFIND,4,JJ,X1,Y1,X2,Y2,P1,EX1,EY1,
1442      &        EX2,EY2)   
1443          IF (IFIND.GT.0) THEN
1444             DO ICAN=1,NCAN 
1445                IF (IFIND.EQ.JCAN(5,ICAN).AND.JJ.EQ.JCAN(4,ICAN)) GOTO 40 
1446                IF (ABS(XM(5,JCAN(5,ICAN))-XM(5,IFIND)).LT.XPREC/10.)
1447      &            GO TO 40 ! elimine les doubles comptages de traces ccc
1448             ENDDO   
1449             CALL STOCK_CANDIDAT(4,JJ,5,IFIND,EX1,EY1,EPH,EAL,NCAN,3)
1450          ENDIF
1451 40       CONTINUE           
1452       ENDDO
1453
1454
1455       NMU45 = 0
1456       DO ICAN = 1,NCAN
1457          JJ4 = JCAN(4,ICAN)  
1458          JJ5 = JCAN(5,ICAN)  
1459          IF (JJ4.GT.0.AND.JJ5.GT.0) THEN
1460             ID4 = ID(IP(4,JJ4))
1461             ID5 = ID(IP(5,JJ5))
1462             IT = ITRACK(IP(5,JJ5))
1463             IF (ITCHECK(IT).EQ.1) THEN
1464                IF (ID4.EQ.ID5) THEN
1465                   NMU45 = NMU45 + 1
1466                ENDIF   
1467             ENDIF
1468          ENDIF   
1469       ENDDO
1470       IF (NMU45.GE.2) NRES(1) = NRES(1)+1   
1471
1472       IF (IDEBUG.GE.2) THEN         
1473          PRINT *,'RECO_TRACKF nb candidat recherche 4->5 et 5->4 :'
1474      &        ,NCAN
1475          PRINT *,'RECO_TRACKF nb of good muons 4->5 et 5->4 :'
1476      &        ,NMU45
1477         
1478       ENDIF
1479
1480       NMU345 = 0
1481 *
1482 * -- Boucle sur les candidats (4,5) NCAN
1483 *         
1484       DO I = 1,NBSTATION
1485          NERR(I) = 0
1486       ENDDO
1487       NMUONF = 0
1488       NGHOST = 0
1489       NTRACKF = 0
1490      
1491       DO ICAN = 1,NCAN
1492          JJ1 = 0
1493          JJ2 = 0
1494          JJ3 = 0
1495          DO ICH = 1,NBSTATION
1496             MPOS(ICH) = 0
1497             MANG(ICH) = 0
1498          ENDDO
1499          MPOS(4) = 1
1500          MPOS(5) = 1
1501          MANG(4) = 1
1502          MANG(5) = 1
1503          IF (JCANTYP(ICAN).EQ.2) MANG(4)=0
1504          IF (JCANTYP(ICAN).EQ.3) MANG(5)=0
1505          EEXM(5) = EEX(ICAN)/SIGCUT
1506          EEYM(5) = EEY(ICAN)/SIGCUT
1507          EEPH(5) = EEP(ICAN)/SIGCUT
1508          EEAL(5) = EEA(ICAN)/SIGCUT
1509          EEXM(4) = EEX(ICAN)/SIGCUT
1510          EEYM(4) = EEY(ICAN)/SIGCUT
1511          EEPH(4) = EEP(ICAN)/SIGCUT
1512          EEAL(4) = EEA(ICAN)/SIGCUT
1513          JJ5 = JCAN(5,ICAN)
1514          JJ4 = JCAN(4,ICAN)
1515          P = PTOT(IP(5,JJ5))
1516          IF (IZM(4,JJ4).EQ.1) THEN
1517             ZPL4 = ZPLANE(4)
1518          ELSE
1519             ZPL4 = ZPLANE(4)-DZ_PL(4)
1520          ENDIF   
1521          IF (IZM(5,JJ5).EQ.1) THEN
1522             ZPL5 = ZPLANE(5)
1523          ELSE
1524             ZPL5 = ZPLANE(5)-DZ_PL(5)
1525          ENDIF   
1526          TPHEST = (XM(5,JJ5) - XM(4,JJ4))/(ZPL5-ZPL4)
1527          PHEST = ATAN(TPHEST)
1528          TALEST = -(YM(5,JJ5) - YM(4,JJ4))*COS(PHEST)
1529      &               /(ZPL5-ZPL4)
1530          PXZEST = -HTOP/(XM(4,JJ4) - ZPL4*TPHEST)
1531          PEST(1) = 1.0/PXZEST
1532          PEST(2) = TPHEST - ALPHATOP/PXZEST ! PHI emission au vertex
1533          PEST(3) = TALEST       ! tan(lambda)     !h=zm*ang deviation alpha
1534          PEST(4) = 0.0          !alpha=qbl/pxz
1535          PEST(5) = 0.0
1536          PSTEP(1) = 0.003       ! =d(1/p)=delta(p)/p**2
1537          PSTEP(2) = 0.001       ! 0.5 degre
1538          PSTEP(3) = 0.001       ! 0.5 degre
1539          PSTEP(4) = 0.0
1540          PSTEP(5) = 1.0
1541          XMES(4) = XM(4,JJ4)
1542          YMES(4) = YM(4,JJ4)
1543          IZMES(4) = IZM(4,JJ4)
1544          PHMES(4) = PHM(4,JJ4)
1545          ALMES(4) = ALM(4,JJ4)
1546          XMES(5) = XM(5,JJ5)
1547          YMES(5) = YM(5,JJ5)
1548          IZMES(5) = IZM(5,JJ5)
1549          PHMES(5) = PHM(5,JJ5)
1550          ALMES(5) = ALM(5,JJ5)
1551
1552          IVERTEX = 0 ! Vertex = (0,0,0)     
1553
1554 * -- Fit 4,5,V pour la recherche dans 3
1555
1556          CALL TRACKF_FIT(IVERTEX,PEST,PSTEP,PXZINV45,TPHI45,TALAM45,
1557      &        XVERT,YVERT)
1558          
1559 * -- Recherche dans la station 3
1560          
1561          P2 = (1.0D0 + TALAM45**2)/(PXZINV45**2) ! P**2
1562 **         EPH=SIGCUT*SQRT(3.6D-6+0.011/P2)
1563 **         EAL=SIGCUT*SQRT(6.85D-4)
1564 **         EXM=SIGCUT*SQRT(0.034+446./P2)
1565 **         EYM=SIGCUT*SQRT(0.049+354./P2)
1566          EPH=SIGCUT*SQRT(4.1D-7+0.015/P2)
1567          EAL=SIGCUT*SQRT(1.04D-4)
1568          EXM=SIGCUT*SQRT(0.0+459./P2)
1569          EYM=SIGCUT*SQRT(0.042+345./P2)
1570          EEXM(3) = EXM/SIGCUT
1571          EEYM(3) = EYM/SIGCUT
1572          EEPH(3) = EPH/SIGCUT
1573          EEAL(3) = EAL/SIGCUT
1574 ** 29.07
1575          CALL CHF1(2301,SNGL(SQRT(P2)),1.)
1576          P1 = PTOT(IP(5,JJ5))
1577          CALL CHFILL2(2311,SNGL(P1),SNGL(EPH),1.)
1578          CALL CHFILL2(2312,SNGL(P1),SNGL(EAL),1.)
1579          CALL CHFILL2(2313,SNGL(P1),SNGL(EXM),1.)
1580          CALL CHFILL2(2314,SNGL(P1),SNGL(EYM),1.)
1581
1582          CALL CHFILL2(2315,SNGL(SQRT(P2)),SNGL(EPH),1.)
1583          CALL CHFILL2(2316,SNGL(SQRT(P2)),SNGL(EAL),1.)
1584          CALL CHFILL2(2317,SNGL(SQRT(P2)),SNGL(EXM),1.)
1585          CALL CHFILL2(2318,SNGL(SQRT(P2)),SNGL(EYM),1.)
1586
1587          CALL CHFILL2(2302,SNGL(SQRT(P2)),SNGL((EPH/SIGCUT)**2),1.)
1588          CALL CHFILL2(2303,SNGL(SQRT(P2)),SNGL((EAL/SIGCUT)**2),1.)
1589          CALL CHFILL2(2304,SNGL(SQRT(P2)),SNGL((EXM/SIGCUT)**2),1.)
1590          CALL CHFILL2(2305,SNGL(SQRT(P2)),SNGL((EYM/SIGCUT)**2),1.)
1591
1592          CALL CHFILL2(2307,SNGL(P1),SNGL((EPH/SIGCUT)**2),1.)
1593          CALL CHFILL2(2308,SNGL(P1),SNGL((EAL/SIGCUT)**2),1.)
1594          CALL CHFILL2(2309,SNGL(P1),SNGL((EXM/SIGCUT)**2),1.)
1595          CALL CHFILL2(2310,SNGL(P1),SNGL((EYM/SIGCUT)**2),1.)
1596 ** end 29.07
1597
1598 **         DO IEX = 4,5
1599 **            EEXM(IEX) = EEXM(3) 
1600 **            EEYM(IEX) = EEYM(3)
1601 **            EEPH(IEX) = EEPH(3)
1602 **            EEAL(IEX) = EEAL(3)
1603 **         ENDDO   
1604          X1 = XPL(3,1)
1605          Y1 = YPL(3,1) 
1606          X2 = XPL(3,2)
1607          Y2 = YPL(3,2) 
1608          PHI1 = PHPL(3)
1609          ALAM1 = ALPL(3)
1610          
1611          CALL DISTMIN4(X1,Y1,PHI1,ALAM1,3,EXM,EYM,EPH,EAL,IPA3)
1612          
1613          P1 = PTOT(IP(5,JJ5))
1614          
1615          CALL CHECK_HISTO4(61,3,IPA3,5,JJ5,X1,Y1,PHI1,ALAM1,P1,
1616      &        EXM,EYM,EPH,EAL)
1617          
1618          IF (IPA3.NE.0) THEN
1619             JJ3 = IPA3
1620             XMES(3) = XM(3,JJ3)
1621             YMES(3) = YM(3,JJ3)
1622             IZMES(3) = IZM(3,JJ3)
1623             PHMES(3) = PHM(3,JJ3)
1624             ALMES(3) = ALM(3,JJ3)
1625             MPOS(3) = 1
1626             MANG(3) = 1
1627             GO TO 124
1628          ELSE
1629             CALL DISTMIN2(X1,Y1,X2,Y2,3,EXM,EYM,0.D0,0.D0,IP3)
1630             CALL CHECK_HISTO2(0,3,IP3,5,JJ5,X1,Y1,X2,Y2,P1,EXM,EYM,
1631      &           0.D0,0.D0)   
1632          ENDIF
1633          IF (IP3.NE.0) THEN
1634             JJ3 = IP3
1635             XMES(3) = XM(3,JJ3)
1636             YMES(3) = YM(3,JJ3)
1637             IZMES(3) = IZM(3,JJ3)
1638             MPOS(3) = 1
1639             MANG(3) = 0
1640          ELSE
1641             GO TO 123
1642          ENDIF     
1643          
1644          
1645 *     -- Fit 3,4,5 pour la recherche dans 2 
1646          
1647  124     CONTINUE
1648          
1649          IF (JJ5.GT.0.AND.JJ4.GT.0.AND.JJ3.GT.0) THEN           
1650             ID4 = ID(IP(4,JJ4))
1651             ID5 = ID(IP(5,JJ5))
1652             ID3 = ID(IP(3,JJ3))
1653             IT = ITRACK(IP(5,JJ5))
1654             IF (ITCHECK(IT).EQ.1) THEN
1655                IF ((ID5.EQ.ID3).AND.(ID5.EQ.ID4)) THEN
1656                   NMU345 = NMU345 + 1
1657                ENDIF   
1658             ENDIF
1659          ENDIF   
1660
1661          IVERTEX = 0    
1662
1663          PEST(1) = PXZINV45
1664          PEST(2) = TPHI45
1665          PEST(3) = TALAM45
1666          PEST(4) = 0.
1667          PEST(5) = 0.
1668          PSTEP(1) = 0.003
1669          PSTEP(2) = 0.001
1670          PSTEP(3) = 0.001
1671          PSTEP(4) = 1.
1672          PSTEP(5) = 1. 
1673          
1674          CALL TRACKF_FIT(IVERTEX,PEST,PSTEP,PXZINV345,TPHI345,
1675      &        TALAM345,XVERT,YVERT)
1676          
1677 *     -- Recherche dans la st. 2
1678          
1679          P2 = (1.0D0 + TALAM345**2)/(PXZINV345**2) ! P**2
1680 **         EPH=SIGCUT*SQRT(3.63D-6+4.8D-3/P2)
1681 **         EAL=SIGCUT*SQRT(6.49D-4)
1682 **         EXM=SIGCUT*SQRT(1.64D-2+821./P2)
1683 **         EYM=SIGCUT*SQRT(4.83D-2+866./P2)
1684          EPH=SIGCUT*SQRT(5.78D-7+5.97D-3/P2)
1685          EAL=SIGCUT*SQRT(1.D-4)
1686          EXM=SIGCUT*SQRT(0.+1453./P2)
1687          EYM=SIGCUT*SQRT(2.78D-2+504./P2)
1688          EEXM(2) = EXM/SIGCUT
1689          EEYM(2) = EYM/SIGCUT
1690          EEPH(2) = EPH/SIGCUT
1691          EEAL(2) = EAL/SIGCUT
1692 ** 29.07
1693          CALL CHF1(2201,SNGL(SQRT(P2)),1.)
1694          P1 = PTOT(IP(5,JJ5))
1695          CALL CHFILL2(2211,SNGL(P1),SNGL(EPH),1.)
1696          CALL CHFILL2(2212,SNGL(P1),SNGL(EAL),1.)
1697          CALL CHFILL2(2213,SNGL(P1),SNGL(EXM),1.)
1698          CALL CHFILL2(2214,SNGL(P1),SNGL(EYM),1.)
1699
1700          CALL CHFILL2(2215,SNGL(SQRT(P2)),SNGL(EPH),1.)
1701          CALL CHFILL2(2216,SNGL(SQRT(P2)),SNGL(EAL),1.)
1702          CALL CHFILL2(2217,SNGL(SQRT(P2)),SNGL(EXM),1.)
1703          CALL CHFILL2(2218,SNGL(SQRT(P2)),SNGL(EYM),1.)
1704
1705          CALL CHFILL2(2202,SNGL(SQRT(P2)),SNGL((EPH/SIGCUT)**2),1.)
1706          CALL CHFILL2(2203,SNGL(SQRT(P2)),SNGL((EAL/SIGCUT)**2),1.)
1707          CALL CHFILL2(2204,SNGL(SQRT(P2)),SNGL((EXM/SIGCUT)**2),1.)
1708          CALL CHFILL2(2205,SNGL(SQRT(P2)),SNGL((EYM/SIGCUT)**2),1.)
1709
1710          CALL CHFILL2(2207,SNGL(P1),SNGL((EPH/SIGCUT)**2),1.)
1711          CALL CHFILL2(2208,SNGL(P1),SNGL((EAL/SIGCUT)**2),1.)
1712          CALL CHFILL2(2209,SNGL(P1),SNGL((EXM/SIGCUT)**2),1.)
1713          CALL CHFILL2(2210,SNGL(P1),SNGL((EYM/SIGCUT)**2),1.)
1714 ** end 29.07
1715
1716 **         DO IEX = 3,5
1717 **            EEXM(IEX) = EEXM(2) 
1718 **            EEYM(IEX) = EEYM(2)
1719 **            EEPH(IEX) = EEPH(2)
1720 **            EEAL(IEX) = EEAL(2)
1721 **         ENDDO   
1722
1723          X1 = XPL(2,1)
1724          Y1 = YPL(2,1) 
1725          PHI1 = PHPL(2)
1726          ALAM1 = ALPL(2)
1727          CALL DISTMIN4(X1,Y1,PHI1,ALAM1,2,EXM,EYM,EPH,EAL,IPA2)
1728          P1 = PTOT(IP(5,JJ5))
1729          CALL CHECK_HISTO4(21,2,IPA2,5,JJ5,X1,Y1,PHI1,ALAM1,P1,
1730      &        EXM,EYM,EPH,EAL)   
1731 *     -- Recherche dans la st. 1
1732 **         EPH=SIGCUT*SQRT(3.54D-6+4.49D-3/P2)
1733 **         EAL=SIGCUT*SQRT(6.14D-4)
1734 **         EXM=SIGCUT*SQRT(6.43D-2+986./P2)  
1735 **         EYM=SIGCUT*SQRT(4.66D-2+1444./P2)  
1736          EPH=SIGCUT*SQRT(4.96D-7+5.87D-3/P2)
1737          EAL=SIGCUT*SQRT(1.D-4)
1738          EXM=SIGCUT*SQRT(6.98D-4+1467./P2)  
1739          EYM=SIGCUT*SQRT(5.22D-2+1013./P2)  
1740          EEXM(1) = EXM/SIGCUT
1741          EEYM(1) = EYM/SIGCUT
1742          EEPH(1) = EPH/SIGCUT
1743          EEAL(1) = EAL/SIGCUT
1744 ** 29.07
1745          CALL CHFILL2(2111,SNGL(P1),SNGL(EPH),1.)
1746          CALL CHFILL2(2112,SNGL(P1),SNGL(EAL),1.)
1747          CALL CHFILL2(2113,SNGL(P1),SNGL(EXM),1.)
1748          CALL CHFILL2(2114,SNGL(P1),SNGL(EYM),1.)
1749
1750          CALL CHFILL2(2115,SNGL(SQRT(P2)),SNGL(EPH),1.)
1751          CALL CHFILL2(2116,SNGL(SQRT(P2)),SNGL(EAL),1.)
1752          CALL CHFILL2(2117,SNGL(SQRT(P2)),SNGL(EXM),1.)
1753          CALL CHFILL2(2118,SNGL(SQRT(P2)),SNGL(EYM),1.)
1754
1755          CALL CHFILL2(2102,SNGL(SQRT(P2)),SNGL((EPH/SIGCUT)**2),1.)
1756          CALL CHFILL2(2103,SNGL(SQRT(P2)),SNGL((EAL/SIGCUT)**2),1.)
1757          CALL CHFILL2(2104,SNGL(SQRT(P2)),SNGL((EXM/SIGCUT)**2),1.)
1758          CALL CHFILL2(2105,SNGL(SQRT(P2)),SNGL((EYM/SIGCUT)**2),1.)
1759
1760          CALL CHFILL2(2107,SNGL(P1),SNGL((EPH/SIGCUT)**2),1.)
1761          CALL CHFILL2(2108,SNGL(P1),SNGL((EAL/SIGCUT)**2),1.)
1762          CALL CHFILL2(2109,SNGL(P1),SNGL((EXM/SIGCUT)**2),1.)
1763          CALL CHFILL2(2110,SNGL(P1),SNGL((EYM/SIGCUT)**2),1.)
1764 ** end 29.07
1765
1766
1767 **         DO IEX = 2,5
1768 **            EEXM(IEX) = EEXM(1) 
1769 **            EEYM(IEX) = EEYM(1)
1770 **            EEPH(IEX) = EEPH(1)
1771 **            EEAL(IEX) = EEAL(1)
1772 **         ENDDO   
1773
1774          X1 = XPL(1,1)
1775          Y1 = YPL(1,1) 
1776          PHI1 = PHPL(1)
1777          ALAM1 = ALPL(1)
1778          CALL DISTMIN4(X1,Y1,PHI1,ALAM1,1,EXM,EYM,EPH,EAL,IPA1)
1779          CALL CHECK_HISTO4(31,1,IPA1,5,JJ5,X1,Y1,PHI1,ALAM1,P1,
1780      &        EXM,EYM,EPH,EAL)   
1781 *     -- A partir de P+A dans la st. 2     
1782          IPA2PA1 = 0
1783          IF (IPA2.GT.0) THEN
1784             PEST(1) = PXZINV345
1785             PEST(2) = TPHI345
1786             PEST(3) = TALAM345
1787             PEST(4) = 0.
1788             PEST(5) = 0.
1789             PSTEP(1) = 0.003
1790             PSTEP(2) = 0.001
1791             PSTEP(3) = 0.001
1792             PSTEP(4) = 1.
1793             PSTEP(5) = 1. 
1794             XMES(2) = XM(2,IPA2)
1795             YMES(2) = YM(2,IPA2)
1796             IZMES(2) = IZM(2,IPA2)
1797             PHMES(2) = PHM(2,IPA2)
1798             ALMES(2) = ALM(2,IPA2)
1799             MPOS(2) = 1
1800             MANG(2) = 1
1801             IVERTEX = 1 
1802 *     -- Fit V,2,3,4,5 pour la recherche dans 1 
1803             
1804             CALL TRACKF_FIT(IVERTEX,PEST,PSTEP,PXZINV,TPHI,TALAM,
1805      &           XVERT,YVERT)
1806             
1807 *     !!! ATTENTION aux erreurs
1808             
1809             P2 = (1.0D0 + TALAM**2)/(PXZINV**2)
1810 **            EPH=SIGCUT*SQRT(3.15D-6+9.21D-3/P2)
1811 **            EAL=SIGCUT*SQRT(5.94D-4)
1812 **            EXM=SIGCUT*SQRT(4.16D-2+182./P2)
1813 **            EYM=SIGCUT*SQRT(2.13)
1814             EPH=SIGCUT*SQRT(9.58D-7+8.93D-3/P2)
1815             EAL=SIGCUT*SQRT(1.D-4)
1816             EXM=SIGCUT*SQRT(1.92D-2+103.3/P2)
1817             EYM=SIGCUT*SQRT(6.3D-2+81.1/P2)
1818             EEXM(1) = EXM/SIGCUT
1819             EEYM(1) = EYM/SIGCUT
1820             EEPH(1) = EPH/SIGCUT
1821             EEAL(1) = EAL/SIGCUT
1822 ** 29.07
1823          CALL CHF1(2701,SNGL(SQRT(P2)),1.)
1824          CALL CHFILL2(2711,SNGL(P1),SNGL(EPH),1.)
1825          CALL CHFILL2(2712,SNGL(P1),SNGL(EAL),1.)
1826          CALL CHFILL2(2713,SNGL(P1),SNGL(EXM),1.)
1827          CALL CHFILL2(2714,SNGL(P1),SNGL(EYM),1.)
1828
1829          CALL CHFILL2(2715,SNGL(SQRT(P2)),SNGL(EPH),1.)
1830          CALL CHFILL2(2716,SNGL(SQRT(P2)),SNGL(EAL),1.)
1831          CALL CHFILL2(2717,SNGL(SQRT(P2)),SNGL(EXM),1.)
1832          CALL CHFILL2(2718,SNGL(SQRT(P2)),SNGL(EYM),1.)
1833
1834          CALL CHFILL2(2702,SNGL(SQRT(P2)),SNGL((EPH/SIGCUT)**2),1.)
1835          CALL CHFILL2(2703,SNGL(SQRT(P2)),SNGL((EAL/SIGCUT)**2),1.)
1836          CALL CHFILL2(2704,SNGL(SQRT(P2)),SNGL((EXM/SIGCUT)**2),1.)
1837          CALL CHFILL2(2705,SNGL(SQRT(P2)),SNGL((EYM/SIGCUT)**2),1.)
1838
1839
1840          CALL CHFILL2(2707,SNGL(P1),SNGL((EPH/SIGCUT)**2),1.)
1841          CALL CHFILL2(2708,SNGL(P1),SNGL((EAL/SIGCUT)**2),1.)
1842          CALL CHFILL2(2709,SNGL(P1),SNGL((EXM/SIGCUT)**2),1.)
1843          CALL CHFILL2(2710,SNGL(P1),SNGL((EYM/SIGCUT)**2),1.)
1844 ** end 29.07
1845
1846
1847 **            DO IEX = 2,5
1848 **               EEXM(IEX) = EEXM(1) 
1849 **               EEYM(IEX) = EEYM(1)
1850 **               EEPH(IEX) = EEPH(1)
1851 **               EEAL(IEX) = EEAL(1)
1852 **            ENDDO   
1853
1854             X1 = XPL(1,1)
1855             Y1 = YPL(1,1) 
1856             X2 = XPL(1,2)
1857             Y2 = YPL(1,2) 
1858             PHI1 = PHPL(1)
1859             ALAM1 = ALPL(1)
1860             CALL DISTMIN4(X1,Y1,PHI1,ALAM1,1,EXM,EYM,EPH,EAL,
1861      &           IPA2PA1)
1862             CALL CHECK_HISTO4(41,1,IPA2PA1,5,JJ5,X1,Y1,PHI1,ALAM1,
1863      &           P1,EXM,EYM,EPH,EAL)
1864             IF (IPA2PA1.GT.0) THEN
1865                JJ2 = IPA2
1866                JJ1 = IPA2PA1
1867                XMES(1) = XM(1,JJ1)
1868                YMES(1) = YM(1,JJ1)
1869                IZMES(1) = IZM(1,JJ1)
1870                PHMES(1) = PHM(1,JJ1)
1871                ALMES(1) = ALM(1,JJ1)
1872                MPOS(1) = 1
1873                MANG(1) = 1
1874                GOTO 123
1875             ELSE
1876                CALL DISTMIN2(X1,Y1,X2,Y2,1,EXM,EYM,0.D0,0.D0,IPA2P1)
1877                CALL CHECK_HISTO2(0,1,IPA2P1,5,JJ5,X1,Y1,X2,Y2,P1,
1878      &              EXM,EYM,0.D0,0.D0)   
1879             ENDIF 
1880          ENDIF   
1881 *     --  A partir de P+A dans la st. 1     
1882          IPA1PA2 = 0
1883          IF (IPA1.GT.0) THEN
1884             PEST(1) = PXZINV345
1885             PEST(2) = TPHI345
1886             PEST(3) = TALAM345
1887             PEST(4) = 0.
1888             PEST(5) = 0.
1889             PSTEP(1) = 0.003
1890             PSTEP(2) = 0.001
1891             PSTEP(3) = 0.001
1892             PSTEP(4) = 1.
1893             PSTEP(5) = 1. 
1894             XMES(1) = XM(1,IPA1)
1895             YMES(1) = YM(1,IPA1)
1896             IZMES(1) = IZM(1,IPA1)
1897             PHMES(1) = PHM(1,IPA1)
1898             ALMES(1) = ALM(1,IPA1)
1899             MPOS(1) = 1
1900             MANG(1) = 1
1901             MPOS(2) = 0
1902             MANG(2) = 0
1903             IVERTEX = 1 
1904 *     -- Fit V,1,3,4,5 pour la recherche dans 2 
1905             
1906             CALL TRACKF_FIT(IVERTEX,PEST,PSTEP,PXZINV,TPHI,TALAM,
1907      &           XVERT,YVERT) 
1908             
1909 *     !!! ATTENTION aux erreurs
1910             
1911             P2 = (1.0D0 + TALAM**2)/(PXZINV**2)
1912 **            EPH=SIGCUT*SQRT(3.15D-6+9.21D-3/P2)
1913 **            EAL=SIGCUT*SQRT(5.94D-4)
1914 **            EXM=SIGCUT*SQRT(4.16D-2+182./P2)
1915 **            EYM=SIGCUT*SQRT(2.13)
1916             EPH=SIGCUT*SQRT(9.58D-7+8.93D-3/P2)
1917             EAL=SIGCUT*SQRT(1.D-4)
1918             EXM=SIGCUT*SQRT(4.0D-2+11.4/P2)
1919             EYM=SIGCUT*SQRT(5.5D-2+44.35/P2)
1920      
1921             EEXM(2) = EXM/SIGCUT
1922             EEYM(2) = EYM/SIGCUT
1923             EEPH(2) = EPH/SIGCUT
1924             EEAL(2) = EAL/SIGCUT
1925 ** 29.07
1926          CALL CHF1(2801,SNGL(SQRT(P2)),1.)
1927          CALL CHFILL2(2811,SNGL(P1),SNGL(EPH),1.)
1928          CALL CHFILL2(2812,SNGL(P1),SNGL(EAL),1.)
1929          CALL CHFILL2(2813,SNGL(P1),SNGL(EXM),1.)
1930          CALL CHFILL2(2814,SNGL(P1),SNGL(EYM),1.)
1931
1932          CALL CHFILL2(2815,SNGL(SQRT(P2)),SNGL(EPH),1.)
1933          CALL CHFILL2(2816,SNGL(SQRT(P2)),SNGL(EAL),1.)
1934          CALL CHFILL2(2817,SNGL(SQRT(P2)),SNGL(EXM),1.)
1935          CALL CHFILL2(2818,SNGL(SQRT(P2)),SNGL(EYM),1.)
1936
1937          CALL CHFILL2(2802,SNGL(SQRT(P2)),SNGL((EPH/SIGCUT)**2),1.)
1938          CALL CHFILL2(2803,SNGL(SQRT(P2)),SNGL((EAL/SIGCUT)**2),1.)
1939          CALL CHFILL2(2804,SNGL(SQRT(P2)),SNGL((EXM/SIGCUT)**2),1.)
1940          CALL CHFILL2(2805,SNGL(SQRT(P2)),SNGL((EYM/SIGCUT)**2),1.)
1941
1942
1943          CALL CHFILL2(2807,SNGL(P1),SNGL((EPH/SIGCUT)**2),1.)
1944          CALL CHFILL2(2808,SNGL(P1),SNGL((EAL/SIGCUT)**2),1.)
1945          CALL CHFILL2(2809,SNGL(P1),SNGL((EXM/SIGCUT)**2),1.)
1946          CALL CHFILL2(2810,SNGL(P1),SNGL((EYM/SIGCUT)**2),1.)
1947
1948 ** end 29.07
1949
1950
1951
1952
1953 **            DO IEX = 1,5
1954 **               EEXM(IEX) = EEXM(2) 
1955 **               EEYM(IEX) = EEYM(2)
1956 **               EEPH(IEX) = EEPH(2)
1957 **               EEAL(IEX) = EEAL(2)
1958 **            ENDDO   
1959
1960             X1 = XPL(2,1)
1961             Y1 = YPL(2,1) 
1962             X2 = XPL(2,2)
1963             Y2 = YPL(2,2) 
1964             PHI1 = PHPL(2)
1965             ALAM1 = ALPL(2)
1966             
1967             CALL DISTMIN4(X1,Y1,PHI1,ALAM1,2,EXM,EYM,EPH,EAL,
1968      &           IPA1PA2)
1969             
1970             CALL CHECK_HISTO4(51,2,IPA1PA2,5,JJ5,X1,Y1,PHI1,ALAM1,
1971      &           P1,EXM,EYM,EPH,EAL)
1972             IF (IPA1PA2.GT.0) THEN
1973                JJ1 = IPA1
1974                JJ2 = IPA1PA2
1975                XMES(2) = XM(2,JJ2)
1976                YMES(2) = YM(2,JJ2)
1977                IZMES(2) = IZM(2,JJ2)
1978                PHMES(2) = PHM(2,JJ2)
1979                ALMES(2) = ALM(2,JJ2)
1980                MPOS(2) = 1
1981                MANG(2) = 1
1982                GOTO 123
1983             ELSE
1984                CALL DISTMIN2(X1,Y1,X2,Y2,2,EXM,EYM,0.D0,0.D0,IPA1P2)
1985                CALL CHECK_HISTO2(0,2,IPA1P2,5,JJ5,X1,Y1,X2,Y2,P1,
1986      &              EXM,EYM,0.D0,0.D0) 
1987             ENDIF 
1988          ENDIF
1989 *     -- Selection par type de candidat            
1990          IF (IPA1PA2.EQ.0.OR.IPA2PA1.EQ.0) THEN
1991             IF (IPA2.GT.0.AND.IPA2P1.GT.0) THEN
1992                JJ2 = IPA2
1993                JJ1 = IPA2P1
1994                XMES(1) = XM(1,JJ1)
1995                YMES(1) = YM(1,JJ1)
1996                IZMES(1) = IZM(1,JJ1)
1997                MPOS(1) = 1
1998                MANG(1) = 0
1999                MPOS(2) = 1
2000                MANG(2) = 1
2001             ELSEIF(IPA1.GT.0.AND.IPA1P2.GT.0) THEN
2002                JJ1 = IPA1
2003                JJ2 = IPA1P2
2004                XMES(2) = XM(2,JJ2)
2005                YMES(2) = YM(2,JJ2)
2006                IZMES(2) = IZM(2,JJ2)
2007                MPOS(1) = 1
2008                MANG(1) = 1
2009                MPOS(2) = 1
2010                MANG(2) = 0
2011             ENDIF
2012          ENDIF
2013          
2014  123     CONTINUE 
2015 ***   
2016          NTRACKFOLD = NTRACKF
2017          NMUONFOLD = NMUONF
2018          NGHOSTOLD = NGHOST
2019          CALL CHECK_MUON(JJ1,JJ2,JJ3,JJ4,JJ5,NTRACKF,NMUONF,
2020      &        NGHOST,NERR,PXZINV,TPHI,TALAM,XVERT,YVERT)
2021          
2022          IF (NTRACKF.GT.NTRACKFOLD) THEN
2023             ISTAT(NTRACKF) = 0 
2024             IF (NMUONF.GT.NMUONFOLD) ISTAT(NTRACKF) = 1 ! Good muon
2025             IF (NGHOST.GT.NGHOSTOLD) ISTAT(NTRACKF) = 2 ! ghost track
2026             PXZOUT(NTRACKF) = 1./PXZINV
2027             TPHIOUT(NTRACKF) = TPHI  
2028             TALAMOUT(NTRACKF) = TALAM 
2029             XVERTOUT(NTRACKF) = XVERT 
2030             YVERTOUT(NTRACKF) = YVERT
2031             PXVOUT(NTRACKF) = PVERT1(IP(5,JJ5))
2032             PYVOUT(NTRACKF) = PVERT2(IP(5,JJ5))
2033             PZVOUT(NTRACKF) = PVERT3(IP(5,JJ5))
2034             CHI2OUT(NTRACKF) = CHI2PL
2035          ENDIF
2036 ***   
2037       ENDDO                     ! end loop on candidats NCAN
2038       
2039       IF (NMU345.GE.2) NRES(2) = NRES(2) + 1
2040
2041       IF (IDEBUG.GE.2) THEN         
2042          PRINT *,'RECO_TRACKF nb of good muons 3 4 5 :'
2043      &        ,NMU345
2044         
2045       ENDIF
2046      
2047       IF (IDEBUG.GE.2) THEN
2048          PRINT *,'RECO_TRACKF nb of track /evt     :',NTRACKF        
2049          PRINT *,'RECO_TRACKF nb of good muon /evt :',NMUONF         
2050          PRINT *,'RECO_TRACKF nb of ghost /evt     :',NGHOST         
2051          DO I = 1,4
2052             PRINT *,'RECO_TRACKF nb of error in st',I,'/evt:',NERR(I)
2053          ENDDO 
2054       ENDIF           
2055       IF (NTRMU.GE.2) NRES(5) = NRES(5)+1
2056       IF (NMUONF.GE.2) NRESF = NRESF+1
2057       NMUONALL = NMUONALL+NMUONF
2058       NGHOSTALL = NGHOSTALL+NGHOST
2059       NTRACKFALL = NTRACKFALL+NTRACKF
2060       DO I = 1,4
2061          NERRALL(I) = NERRALL(I)+NERR(I)
2062       ENDDO  
2063       
2064       CALL TRACKF_OUT(IR,NTRACKF) 
2065 ***   
2066       RETURN
2067       END
2068
2069 **********************************************************************         
2070       SUBROUTINE TRACKF_OUT(IR,NTRACKF)
2071 **********************************************************************         
2072 *   
2073 *  
2074 **********************************************************************
2075       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2076
2077       PARAMETER (MAXHITCH=10000,MAXTRK=50000,MAXHITTOT=20000,
2078      &           NBSTATION=5)
2079       PARAMETER (NBCHAMBER=10,NTRMAX=500)
2080   
2081       COMMON/RECOUT/JJO(NBCHAMBER,NTRMAX),XMESO(NBCHAMBER,NTRMAX),
2082      &              YMESO(NBCHAMBER,NTRMAX)
2083
2084       COMMON/TRACKFOUT/IEVOUT,NTREVT,JJOUT(NBCHAMBER,NTRMAX),
2085      &               ISTAT(NTRMAX),PXZOUT(NTRMAX),TPHIOUT(NTRMAX),
2086      &               TALAMOUT(NTRMAX),XVERTOUT(NTRMAX),YVERTOUT(NTRMAX),
2087      &               CHI2OUT(NTRMAX),
2088      &               XMESOUT(NBCHAMBER,NTRMAX),YMESOUT(NBCHAMBER,NTRMAX)
2089      &              ,PXVOUT(NTRMAX),PYVOUT(NTRMAX),PZVOUT(NTRMAX)
2090 **
2091       IEVOUT = IR  
2092       NTREVT = NTRACKF 
2093       IF (NTREVT.GT.0) THEN
2094          DO ITR = 1,NTREVT
2095             ICH = 0
2096             DO IST = 1,NBSTATION
2097                DO ILOOP = 1,2 
2098                   ICH = ICH + 1 
2099                   JJOUT(ICH,ITR) = JJO(ICH,ITR)
2100                   IF (JJOUT(ICH,ITR).GT.0) THEN
2101                       XMESOUT(ICH,ITR) = XMESO(ICH,ITR)
2102                       YMESOUT(ICH,ITR) = YMESO(ICH,ITR)
2103                   ENDIF
2104                ENDDO
2105             ENDDO
2106          ENDDO
2107       ENDIF
2108
2109
2110       RETURN
2111       END
2112 **********************************************************************          
2113       SUBROUTINE CHECK_MUON(JJ1,JJ2,JJ3,JJ4,JJ5,
2114      &                  NTRACKF,NMUONF,NGHOST,NERR,PXZINV,TPHI,TALAM,
2115      &                  XVERT,YVERT)  
2116 **********************************************************************          
2117 *   Check muon track candidate using GEANT informations
2118 *
2119 **********************************************************************
2120       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2121 * --
2122       PARAMETER (MAXHITCH=10000,MAXTRK=50000,MAXHITTOT=20000,
2123      &           NBSTATION=5)
2124       PARAMETER (NBCHAMBER=10,NTRMAX=500)
2125       
2126       COMMON/CHHIT/XM(NBSTATION,MAXHITCH),YM(NBSTATION,MAXHITCH),
2127      &             PHM(NBSTATION,MAXHITCH),ALM(NBSTATION,MAXHITCH),
2128      &             IZM(NBSTATION,MAXHITCH),
2129      &             IP(NBSTATION,MAXHITCH),JHIT(NBSTATION),
2130      &             XMR(NBSTATION,MAXHITCH,2),YMR(NBSTATION,MAXHITCH,2) 
2131       
2132       COMMON/RHIT/ITYP(MAXHITTOT),XTR(MAXHITTOT),YTR(MAXHITTOT),
2133      &            PTOT(MAXHITTOT),ID(MAXHITTOT),IZST(MAXHITTOT),
2134      &            PVERT1(MAXHITTOT),PVERT2(MAXHITTOT),PVERT3(MAXHITTOT),
2135      &            ZVERT(MAXHITTOT),NHITTOT
2136
2137       COMMON/VERIFGEANT/ITTROUGH(MAXTRK,NBSTATION),
2138      &                  IT_LIST(MAXTRK),IT_NP(MAXTRK),ITCHECK(MAXTRK),
2139      &                  ITRACK(MAXHITTOT)
2140      
2141       COMMON/MEASUR/XMES(NBSTATION),YMES(NBSTATION),IZMES(NBSTATION),
2142      &              PHMES(NBSTATION),ALMES(NBSTATION),MPOS(NBSTATION),
2143      &              MANG(NBSTATION)
2144
2145       COMMON/RECOUT/JJO(NBCHAMBER,NTRMAX),XMESO(NBCHAMBER,NTRMAX),
2146      &              YMESO(NBCHAMBER,NTRMAX)
2147
2148       COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)      
2149
2150       COMMON/MAGNET/BL3,B0
2151
2152       REAL*4 R1,R2
2153       DATA R1,R2/0.,1./
2154  
2155       DIMENSION NERR(NBSTATION),JJK(NBSTATION) 
2156       
2157 *      print*,' *** appel de la subroutine check_muon ***'
2158
2159       IF (JJ1.GT.0.AND.JJ2.GT.0.AND.JJ3.GT.0.AND.JJ4.GT.0
2160      &     .AND.JJ5.GT.0) THEN
2161          IDFIND = ID(IP(5,JJ5))
2162          IT = ITRACK(IP(5,JJ5))
2163          JJK(1) = JJ1
2164          JJK(2) = JJ2
2165          JJK(3) = JJ3
2166          JJK(4) = JJ4
2167          JJK(5) = JJ5
2168          NTRACKF = NTRACKF+1
2169          DO I = 1,NBCHAMBER
2170             JJO(I,NTRACKF) = 0
2171          ENDDO
2172          ICH1 = -1
2173          DO IST = 1,5
2174             ICH1 = ICH1 + 2  
2175             IF (XMR(IST,JJK(IST),1).LT.999.) THEN
2176                JJO(ICH1,NTRACKF) = JJK(IST)
2177                XMESO(ICH1,NTRACKF) = XMR(IST,JJK(IST),1)
2178                YMESO(ICH1,NTRACKF) = YMR(IST,JJK(IST),1)
2179                ICH2 = ICH1 + 1  
2180                IF (MANG(IST).EQ.1) THEN
2181                   JJO(ICH2,NTRACKF) = JJK(IST)
2182                   XMESO(ICH2,NTRACKF) = XMR(IST,JJK(IST),2)
2183                   YMESO(ICH2,NTRACKF) = YMR(IST,JJK(IST),2)
2184                ENDIF
2185             ELSE
2186                ICH2 = ICH1+1
2187                JJO(ICH2,NTRACKF) = JJK(IST)
2188                XMESO(ICH2,NTRACKF) = XMR(IST,JJK(IST),2)
2189                YMESO(ICH2,NTRACKF) = YMR(IST,JJK(IST),2)
2190             ENDIF
2191          ENDDO
2192          
2193          CALL CHFILL(700,SNGL(XVERT),R1,R2)
2194          CALL CHFILL(701,SNGL(YVERT),R1,R2)
2195          
2196          IF (ITCHECK(IT).EQ.1) THEN
2197             
2198             
2199             IF (ID(IP(1,JJ1)).EQ.IDFIND .AND.
2200      &           ID(IP(2,JJ2)).EQ.IDFIND .AND.
2201      &           ID(IP(3,JJ3)).EQ.IDFIND .AND.
2202      &           ID(IP(4,JJ4)).EQ.IDFIND) THEN
2203                
2204                NMUONF = NMUONF+1 ! Bon muon
2205                
2206                PT = SQRT(PVERT1(IP(5,JJ5))**2+PVERT2(IP(5,JJ5))**2)
2207                PZ =  PVERT3(IP(5,JJ5))
2208                E = SQRT(PT**2+PZ**2+0.1056**2)
2209                YRAP = 0.5*DLOG((E+PZ)/(E-PZ))
2210                PHIM = DATAN2(PVERT2(IP(5,JJ5)),PVERT1(IP(5,JJ5)))
2211                CALL CHFILL(801,SNGL(YRAP),R1,R2)
2212                CALL CHFILL(901,SNGL(PT),R1,R2)
2213                CALL CHFILL(911,SNGL(PHIM),R1,R2)
2214             ELSE
2215                
2216                NGHOST = NGHOST+1 ! ghost
2217                
2218                PT = SQRT(PVERT1(IP(5,JJ5))**2+PVERT2(IP(5,JJ5))**2)
2219                PZ =  PVERT3(IP(5,JJ5))
2220                E = SQRT(PT**2+PZ**2+0.1056**2)
2221                YRAP = 0.5*DLOG((E+PZ)/(E-PZ))
2222                PHIM = DATAN2(PVERT2(IP(5,JJ5)),PVERT1(IP(5,JJ5)))
2223                CALL CHFILL(802,SNGL(YRAP),R1,R2)
2224                CALL CHFILL(902,SNGL(PT),R1,R2)
2225                CALL CHFILL(912,SNGL(PHIM),R1,R2)
2226             ENDIF
2227          ENDIF    
2228       ENDIF
2229       IF (ITCHECK(ITRACK(IP(5,JJ5))).EQ.1) THEN
2230          IF (JJ3.EQ.0) NERR(3) = NERR(3)+1
2231          IF (JJ3.NE.0) THEN
2232             IF (JJ1.EQ.0) NERR(1) = NERR(1)+1
2233             IF (JJ2.EQ.0) NERR(2) = NERR(2)+1
2234          ENDIF   
2235          PT = SQRT(PVERT1(IP(5,JJ5))**2+PVERT2(IP(5,JJ5))**2)
2236          PZ =  PVERT3(IP(5,JJ5))
2237          E = SQRT(PT**2+PZ**2+0.1056**2)
2238          YRAP = 0.5*DLOG((E+PZ)/(E-PZ))
2239          PHIM = DATAN2(PVERT2(IP(5,JJ5)),PVERT1(IP(5,JJ5)))
2240          CALL CHFILL(800,SNGL(YRAP),R1,R2)
2241          CALL CHFILL(900,SNGL(PT),R1,R2)
2242          CALL CHFILL(910,SNGL(PHIM),R1,R2)
2243       ENDIF
2244       
2245       RETURN
2246       END
2247
2248
2249 *************************************************************************
2250       SUBROUTINE RECO_SUM
2251 *************************************************************************
2252 *
2253 *
2254 *************************************************************************
2255       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2256 *
2257       PARAMETER (MAXHITCH=10000,MAXTRK=50000,MAXHITTOT=20000,
2258      &           NBSTATION=5)
2259       PARAMETER (NBCHAMBER=10,NTRMAX=500)
2260 *
2261       COMMON/TRACKFOUT/IEVOUT,NTREVT,JJOUT(NBCHAMBER,NTRMAX),
2262      &               ISTAT(NTRMAX),PXZOUT(NTRMAX),TPHIOUT(NTRMAX),
2263      &               TALAMOUT(NTRMAX),XVERTOUT(NTRMAX),YVERTOUT(NTRMAX),
2264      &               CHI2OUT(NTRMAX), 
2265      &               XMESOUT(NBCHAMBER,NTRMAX),YMESOUT(NBCHAMBER,NTRMAX)
2266      &              ,PXVOUT(NTRMAX),PYVOUT(NTRMAX),PZVOUT(NTRMAX)
2267 *
2268       COMMON/CHHIT/XM(NBSTATION,MAXHITCH),YM(NBSTATION,MAXHITCH),
2269      &             PHM(NBSTATION,MAXHITCH),ALM(NBSTATION,MAXHITCH),
2270      &             IZM(NBSTATION,MAXHITCH),
2271      &             IP(NBSTATION,MAXHITCH),JHIT(NBSTATION),
2272      &             XMR(NBSTATION,MAXHITCH,2),YMR(NBSTATION,MAXHITCH,2) 
2273 *      
2274       COMMON/RHIT/ITYP(MAXHITTOT),XTR(MAXHITTOT),YTR(MAXHITTOT),
2275      &            PTOT(MAXHITTOT),ID(MAXHITTOT),IZST(MAXHITTOT),
2276      &            PVERT1(MAXHITTOT),PVERT2(MAXHITTOT),PVERT3(MAXHITTOT),
2277      &            ZVERT(MAXHITTOT),NHITTOT
2278 *  
2279       COMMON/RECOUT/JJO(NBCHAMBER,NTRMAX),XMESO(NBCHAMBER,NTRMAX),
2280      &              YMESO(NBCHAMBER,NTRMAX)
2281 *
2282       COMMON/TRACKSUM/NRES(5),NRESF,NTRMUALL,NMUONALL,NGHOSTALL,
2283      &     NTRACKFALL,NERRALL(NBSTATION),IR                
2284 *
2285       COMMON/PRECSUM/NRESF1,NMUONALL1,NGHOSTALL1,NTRACKFALL1          
2286 *
2287       REAL*4 PXR,PYR,PZR,ZVR,CHI2R,PXV,PYV,PZV
2288       COMMON/PAWCR4/IEVR,NTRACKR,ISTATR(NTRMAX),ISIGNR(NTRMAX),
2289      &              PXR(NTRMAX),PYR(NTRMAX),PZR(NTRMAX),ZVR(NTRMAX),
2290      &              CHI2R(NTRMAX),PXV(NTRMAX),PYV(NTRMAX),PZV(NTRMAX)
2291 *
2292       COMMON/DEBEVT/IDEBUG
2293 *  
2294       NMUF = 0
2295       NGHF = 0
2296       DO ITR = 1,NTREVT
2297          NTRACKFALL1 =  NTRACKFALL1 + 1
2298          IF (ISTAT(ITR).EQ.1) THEN
2299             NMUF = NMUF + 1
2300             NMUONALL1 = NMUONALL1 + 1
2301          ELSEIF (ISTAT(ITR).EQ.2) THEN 
2302             NGHF = NGHF + 1
2303             NGHOSTALL1 = NGHOSTALL1 + 1
2304          ENDIF
2305       ENDDO
2306       IF (NMUF.EQ.2) NRESF1 = NRESF1 + 1 
2307 *
2308       NTRACKR = NTREVT
2309       DO ITR = 1,NTREVT
2310          ISTATR(ITR) = ISTAT(ITR)
2311          ISIGNR(ITR) = 1 
2312          IF (PXZOUT(ITR).LT.0.) ISIGNR(ITR) = -1
2313          PXZ = ABS(PXZOUT(ITR))
2314          PHI = ATAN(TPHIOUT(ITR)) 
2315          ALAM = ATAN(TALAMOUT(ITR)) 
2316          PYR(ITR) = PXZ*SIN(PHI)
2317          PXR(ITR) = PXZ*TAN(ALAM)
2318          PZR(ITR) = PXZ*COS(PHI)
2319          ZVR(ITR) = 0.
2320          CHI2R(ITR) =  CHI2OUT(ITR)
2321          PXV(ITR) = PXVOUT(ITR)
2322          PYV(ITR) = PYVOUT(ITR) 
2323          PZV(ITR) = PZVOUT(ITR)         
2324       ENDDO     
2325
2326       CALL CHFNT(IEVR,NTRACKR,ISTATR,ISIGNR,
2327      &     PXR,PYR,PZR,ZVR,CHI2R,PXV,PYV,PZV)
2328
2329       IF (IDEBUG.GE.2) THEN 
2330          PRINT *,'RECO_SUM evt number :',IEVOUT
2331          PRINT *,'RECO_SUM nb of track /evt     :',NTREVT        
2332          PRINT *,'RECO_SUM nb of good muon /evt :',NMUF         
2333          PRINT *,'RECO_SUM nb of ghost /evt     :',NGHF         
2334          IF (NTREVT.GT.0) THEN
2335             DO ITR = 1,NTREVT
2336                PRINT *,'RECO_SUM track number :',ITR
2337                PRINT *,'RECO_SUM PXZOUT :',PXZOUT(ITR)
2338                PRINT *,'RECO_SUM CHI2OUT :',CHI2OUT(ITR)
2339                PRINT *,' PX GEANT= ', PXV(ITR),' PX RECONS= ',PYR(ITR)
2340                PRINT *,' PY GEANT= ', PYV(ITR),' PY RECONS= ',PXR(ITR)
2341                PRINT *,' PZ GEANT= ', PZV(ITR),' PZ RECONS= ',PZR(ITR)
2342                PXZV = SQRT( PYV(ITR)**2+PZV(ITR)**2)
2343                PXZR = SQRT( PXR(ITR)**2+PZR(ITR)**2) 
2344                PRINT *,' PXZ GEANT= ', PXZV,' PXZ RECONS= ',PXZR
2345                FIV=ATAN2(PYV(ITR),PZV(ITR))
2346                ALAMV=ATAN2(PXV(ITR),PXZV)
2347                FIR=ATAN2(PXR(ITR),PZR(ITR))
2348                ALAMR=ATAN2(PYR(ITR),PXZR)
2349 **               PRINT *,' PHI GEANT= ',FIV,' PXZ RECONS= ',FIR
2350 **               PRINT *,' ALAM GEANT= ',ALAMV,' ALAM RECONS= ',ALAMR
2351           ENDDO
2352          ENDIF
2353       ENDIF 
2354
2355       RETURN
2356       END
2357
2358 *************************************************************************
2359       SUBROUTINE RECO_TERM
2360 *************************************************************************
2361 *
2362 *
2363 *************************************************************************
2364       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2365 *
2366       PARAMETER(NBSTATION=5)
2367 *
2368       COMMON/TRACKSUM/NRES(5),NRESF,NTRMUALL,NMUONALL,NGHOSTALL,
2369      &     NTRACKFALL,NERRALL(NBSTATION),IR  
2370 *              
2371       COMMON/PRECSUM/NRESF1,NMUONALL1,NGHOSTALL1,NTRACKFALL1          
2372 *
2373       COMMON/DEBEVT/IDEBUG
2374 *
2375       CHARACTER*50 FILEBKG,FILERES,FILEOUT,FILEMIN
2376 *
2377       IF (IDEBUG.GE.1) THEN
2378          PRINT *,'    '
2379          PRINT *,'RECO_TERM ***** SUMMARY TRACK-FINDING *****'
2380          PRINT *,'RECO_TERM nb of resonances :',NRES(5)
2381          PRINT *,'RECO_TERM nb of resonances  45 :',NRES(1)
2382          PRINT *,'RECO_TERM nb of resonances 345 :',NRES(2)
2383          PRINT *,'RECO_TERM nb of resonances found :',NRESF
2384          PRINT *,'RECO_TERM nb of muon track       :',NTRMUALL
2385          PRINT *,'RECO_TERM nb of track found      :',NTRACKFALL
2386          PRINT *,'RECO_TERM nb of muon track found :',NMUONALL
2387          PRINT *,'RECO_TERM nb of ghost found      :',NGHOSTALL
2388          DO I = 1,4
2389             PRINT *,'RECO_TERM nb of error in st',I,':',NERRALL(I)
2390          ENDDO
2391
2392          PRINT *,'    '
2393          PRINT *,'RECO_TERM ***** SUMMARY PRECISION *****'
2394          PRINT *,'RECO_TERM nb of resonances found :',NRESF1
2395          PRINT *,'RECO_TERM nb of track found      :',NTRACKFALL1
2396          PRINT *,'RECO_TERM nb of muon track found :',NMUONALL1
2397          PRINT *,'RECO_TERM nb of ghost found      :',NGHOSTALL1
2398       ENDIF
2399
2400       CALL HIST_CLOSED
2401
2402       RETURN
2403       END
2404
2405 *************************************************************************
2406       SUBROUTINE RECO_PRECISION
2407 *************************************************************************
2408 *
2409 *
2410 *************************************************************************
2411       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2412 *
2413       PARAMETER (MAXHITCH=10000,MAXTRK=50000,MAXHITTOT=20000,
2414      &           NBSTATION=5,MAXCAN=1000,NTRMAX=500)
2415       PARAMETER (NPLANE=10) 
2416 *
2417       COMMON/PARAM/ZPLANEP(NPLANE),THICK,XPREC,YPREC,B0,BL3,ZMAGS,
2418      &             ZMAGE,ZABS,XMAG,ZBP1,ZBP2,CONST
2419 *      
2420       COMMON/CHHIT/XM(NBSTATION,MAXHITCH),YM(NBSTATION,MAXHITCH),
2421      &             PHM(NBSTATION,MAXHITCH),ALM(NBSTATION,MAXHITCH),
2422      &             IZM(NBSTATION,MAXHITCH),
2423      &             IP(NBSTATION,MAXHITCH),JHIT(NBSTATION),
2424      &             XMR(NBSTATION,MAXHITCH,2),YMR(NBSTATION,MAXHITCH,2) 
2425 *      
2426       COMMON/RHIT/ITYP(MAXHITTOT),XTR(MAXHITTOT),YTR(MAXHITTOT),
2427      &            PTOT(MAXHITTOT),ID(MAXHITTOT),IZST(MAXHITTOT),
2428      &            PVERT1(MAXHITTOT),PVERT2(MAXHITTOT),PVERT3(MAXHITTOT),
2429      &            ZVERT(MAXHITTOT),NHITTOT
2430 *
2431       COMMON/TRACKFOUT/IEVOUT,NTREVT,JJOUT(NPLANE,NTRMAX),
2432      &               ISTAT(NTRMAX),PXZOUT(NTRMAX),TPHIOUT(NTRMAX),
2433      &               TALAMOUT(NTRMAX),XVERTOUT(NTRMAX),YVERTOUT(NTRMAX),
2434      &               CHI2OUT(NTRMAX),
2435      &               XMESOUT(NPLANE,NTRMAX),YMESOUT(NPLANE,NTRMAX)
2436      &              ,PXVOUT(NTRMAX),PYVOUT(NTRMAX),PZVOUT(NTRMAX)
2437
2438 *
2439       COMMON/MEAS/LPLANE(NPLANE),XMP(NPLANE),YMP(NPLANE)
2440 *      
2441       COMMON/FCNOUT/PXZEA,ALAMEA,PHIEA,XEA,YEA,NPLU,CHI2
2442 *
2443       COMMON/PRECCUT/PCUT,PTCUT,CHI2CUT
2444 *
2445       DIMENSION PARMU(MAXCAN,NPLANE,2),LPLANEMU(MAXCAN,NPLANE)  
2446 *
2447       IF (NTREVT.EQ.0) RETURN
2448
2449       DO ITR = 1,NTREVT
2450          ICH = 0
2451          DO IST = 1,NBSTATION
2452             DO ILOOP = 1,2 
2453                ICH = ICH + 1 
2454 **               print *,' ich=',ich
2455                IF (JJOUT(ICH,ITR).GT.0) THEN
2456                   LPLANEMU(ITR,ICH) = 1
2457                   PARMU(ITR,ICH,1) = XMESOUT(ICH,ITR)
2458                   PARMU(ITR,ICH,2) = YMESOUT(ICH,ITR)
2459 **                  print *,' x,y ', PARMU(ITR,ICH,1),PARMU(ITR,ICH,2) 
2460                ELSE
2461                   LPLANEMU(ITR,ICH) = 0
2462                ENDIF
2463             ENDDO
2464          ENDDO
2465       ENDDO
2466 *
2467       NTRACK = 0
2468       DO ICAN = 1,NTREVT
2469          DO ICH = 1,NPLANE
2470             LPLANE(ICH) = LPLANEMU(ICAN,ICH)
2471             XMP(ICH) = PARMU(ICAN,ICH,1)
2472             YMP(ICH) = PARMU(ICAN,ICH,2)
2473          ENDDO 
2474
2475          IF (LPLANE(1).GT.0) THEN
2476             X1 = XMP(1)
2477             Y1 = YMP(1)
2478             IPL1 = 1
2479          ELSE
2480             X1 = XMP(2)
2481             Y1 = YMP(2)
2482             IPL1 = 2
2483          ENDIF
2484          IF (LPLANE(3).GT.0) THEN
2485             X2 = XMP(3)
2486             Y2 = YMP(3)
2487             IPL2 = 3
2488          ELSE
2489             X2 = XMP(4)
2490             Y2 = YMP(4)
2491             IPL2 = 4
2492          ENDIF
2493          IF (LPLANE(7).GT.0) THEN
2494             X3 = XMP(7)
2495             IPL3 = 7
2496          ELSE
2497             X3 = XMP(8)
2498             IPL3 = 8
2499          ENDIF
2500          IF (LPLANE(9).GT.0) THEN
2501             X4 = XMP(9)
2502             IPL4 = 9
2503          ELSE
2504             X4 = XMP(10)
2505             IPL4 = 10
2506          ENDIF
2507    
2508          PHIAV = DATAN2((X2-X1),(ZPLANEP(IPL2)-ZPLANEP(IPL1)))          
2509          PHIAP = DATAN2((X4-X3),(ZPLANEP(IPL4)-ZPLANEP(IPL3)))
2510          DPHI = (PHIAP-PHIAV)
2511          ASIGN = 1.
2512          IF (DPHI.LT.0.) ASIGN = -1. ! CCC
2513          PXZ = CONST/DABS(DPHI)
2514 ** Cuts PXZ           
2515          IF (PXZ.LT.PCUT) GO TO 2       
2516             
2517          PXZINVI = ASIGN/PXZ ! CCC
2518 **         PXZINVI = 1./PXZOUT(ICAN) ! CCC
2519 **         PXZINVI = -1./49. ! CCC
2520          PHII = PHIAV
2521          ALAMI = DATAN2((Y2-Y1),DSQRT((X2-X1)**2
2522      &            +(ZPLANEP(IPL2)-ZPLANEP(IPL1))**2))
2523          XVR = X1
2524          YVR = Y1
2525 **         print *,' avant prec_fit pxzi phii alami x y',1./ PXZINVI,
2526 **     &         PHII, ALAMI ,XVR,YVR             
2527 **         PRINT *,' X1 X2 X3 X4',X1,X2,X3,X4
2528 **         PRINT *,' Z1 Z2 Z3 Z4',ZPLANEP(IPL1),ZPLANEP(IPL2),
2529 **     &          ZPLANEP(IPL3),ZPLANEP(IPL4)
2530 **         PRINT *,' CONST= ',CONST 
2531
2532          PRINT *,' RECO_PRECISION CHI2OUT avant fit ',CHI2OUT(ICAN)
2533
2534          IF (CHI2OUT(ICAN).GT.CHI2CUT) GO TO 2
2535
2536 ** Fit des traces apres l'absorbeur           
2537          CALL PREC_FIT (PXZINVI,PHII,ALAMI,XVR,YVR,
2538      &          PXZINVF,PHIF,ALAMF,XVERTF,YVERTF,EPXZINV,EPHI,EALAM,
2539      &          EXVERT,EYVERT)
2540      
2541 ** Correction de Branson       
2542          CALL BRANSON(PXZEA,PHIEA,ALAMEA,XEA,YEA)
2543          
2544          PXZ1 = DABS(PXZEA)
2545          PX1 = PXZ1*DSIN(PHIEA)
2546          PY1 = PXZ1*DTAN(ALAMEA)
2547          PT1 = DSQRT(PX1**2+PY1**2)
2548 ** Cuts PT
2549          IF (PT1.LT.PTCUT) GO TO 2
2550 ** Cuts CHI2
2551          IF ((CHI2/FLOAT(2*NPLU-5)).GT.CHI2CUT) GO TO 2
2552
2553          NTRACK = NTRACK + 1
2554          DO ICH = 1,NPLANE           
2555             JJOUT(ICH,NTRACK) =  JJOUT(ICH,ICAN)
2556             XMESOUT(ICH,NTRACK) =  XMESOUT(ICH,ICAN)
2557             YMESOUT(ICH,NTRACK) =  YMESOUT(ICH,ICAN)
2558          ENDDO
2559          ISTAT(NTRACK) = ISTAT(ICAN)
2560          PXZOUT(NTRACK) = PXZEA
2561          TPHIOUT(NTRACK) = DTAN(PHIEA)
2562          TALAMOUT(NTRACK) = DTAN(ALAMEA) 
2563          XVERTOUT(NTRACK) = XEA        
2564          YVERTOUT(NTRACK) = YEA 
2565          CHI2OUT(NTRACK) =  CHI2/FLOAT(2*NPLU-5) 
2566
2567          PRINT *,' RECO_PRECISION CHI2OUT apres fit',CHI2OUT(NTRACK) 
2568
2569 **         print *,' reco_precision pxz tphi talam xvert yvert chi2',
2570 **     &       PXZEA,PHIEA,ALAMEA,
2571 **     &       XEA,YEA,CHI2/FLOAT(2*NPLU-5)  
2572   2      CONTINUE 
2573       ENDDO
2574       NTREVT = NTRACK 
2575
2576       RETURN
2577       END
2578
2579
2580 ************************************************************************
2581       SUBROUTINE ORDONNE_HIT(ICH,HCUT)
2582 **************************************************************************
2583 *    
2584 *  Sort hits in station ICH according to the "impact parameter" HHIT
2585 *
2586 **************************************************************************
2587  
2588       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2589  
2590       PARAMETER (MAXHITCH=10000,MAXTRK=50000,MAXHITTOT=20000,
2591      &           NBSTATION=5)
2592       
2593       COMMON/CHHIT/XM(NBSTATION,MAXHITCH),YM(NBSTATION,MAXHITCH),
2594      &             PHM(NBSTATION,MAXHITCH),ALM(NBSTATION,MAXHITCH),
2595      &             IZM(NBSTATION,MAXHITCH),
2596      &             IP(NBSTATION,MAXHITCH),JHIT(NBSTATION),
2597      &             XMR(NBSTATION,MAXHITCH,2),YMR(NBSTATION,MAXHITCH,2) 
2598      
2599       COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)      
2600       
2601       COMMON/VERIFGEANT/ITTROUGH(MAXTRK,NBSTATION),
2602      &                  IT_LIST(MAXTRK),IT_NP(MAXTRK),ITCHECK(MAXTRK),
2603      &                  ITRACK(MAXHITTOT)
2604      
2605       COMMON/HCHHIT/HHIT(MAXHITCH),INDEXTAB(MAXHITCH),INDEXMAX
2606 *
2607       COMMON/DEBEVT/IDEBUG
2608
2609       REAL*4 H4(MAXHITCH) 
2610 * tri des impulsion par ordre decroissant
2611 * le tab INDEXTAB contient les j ordonnes
2612 * INDEXMAX est l indice max du tableau = NBHIT si pas de contrainte
2613
2614       JJ=0
2615       INDEXMAX = 0
2616 * boucle sur le nombre de hits candidats de la station
2617
2618       DO J=1,JHIT(ICH)
2619
2620           IF (PHM(ICH,J).LT.6.3) THEN        !2pi=6.3 radian
2621             JJ=JJ+1
2622 * calcul du h dans XOY a z=0
2623             HHIT(J)=ABS(XM(ICH,J)-ZPLANE(ICH)*PHM(ICH,J))
2624 * cut en Pxz
2625             IF (HHIT(J).LT.HCUT) THEN
2626                INDEXMAX=INDEXMAX+1
2627                INDEXTAB(INDEXMAX)=J
2628             ELSEIF(ITCHECK(ITRACK(IP(ICH,J))).EQ.1) THEN
2629                IF (IDEBUG.GE.2) THEN
2630                   PRINT *,'ORDONNE_HIT rejet muon/res in st.',ICH,
2631      &                 ' h=',HHIT(J)
2632                ENDIF
2633             ENDIF
2634          ENDIF
2635       ENDDO
2636
2637       NBHIT=JHIT(ICH)
2638
2639       DO I = 1,NBHIT
2640          H4(I) = SNGL(HHIT(I))
2641       ENDDO   
2642
2643       CALL SORTZV(H4,INDEXTAB,INDEXMAX,-1,0,1)
2644
2645       DO I = 1,NBHIT
2646          HHIT(I) = DBLE(H4(I))
2647       ENDDO   
2648 **      PRINT *,'ORDONNE st. numero',ICH
2649 **      PRINT *,'ORDONNE nb de hits initiaux dans st.',ICH,':',JHIT(ICH)
2650 **      PRINT *,'ORDONNE nb de hits avec mes. angulaire:',JJ
2651 **      PRINT *,'ORDONNE nb de hits avec mes. ang. et cut en Pxz:',INDEXMAX
2652       IF (IDEBUG.GE.2) THEN 
2653          PRINT *,'ORDONNE_HIT nb de hits accepte dans st.',ICH,':',
2654      &           INDEXMAX
2655       ENDIF
2656
2657       RETURN
2658       END
2659 ***********************************************************************************
2660       SUBROUTINE DISTMIN4(X1,Y1,PHI1,ALAM1,ICH,EX,EY,EPHI,ELAM,IFIND)
2661 ***********************************************************************************
2662 *    Find the nearest hit in station ICH in the (X,Y,lambda,phi) phase space
2663 *
2664 ***********************************************************************************
2665  
2666       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2667       
2668       PARAMETER (MAXHITCH=10000,NBSTATION=5)
2669       
2670       COMMON/CHHIT/XM(NBSTATION,MAXHITCH),YM(NBSTATION,MAXHITCH),
2671      &             PHM(NBSTATION,MAXHITCH),ALM(NBSTATION,MAXHITCH),
2672      &             IZM(NBSTATION,MAXHITCH),
2673      &             IP(NBSTATION,MAXHITCH),JHIT(NBSTATION),
2674      &             XMR(NBSTATION,MAXHITCH,2),YMR(NBSTATION,MAXHITCH,2) 
2675  
2676       IFIND = 0
2677       DISTMIN=4.
2678       DO I=1,JHIT(ICH)
2679          IF (PHM(ICH,I).LE.6.3) THEN ! angles measured
2680             IF (ABS(PHI1-PHM(ICH,I)) .LT. EPHI .AND.
2681      &         ABS(ALAM1-ALM(ICH,I)) .LT. ELAM .AND.
2682      &         ABS(X1-XM(ICH,I)) .LT. EX .AND.
2683      &         ABS(Y1-YM(ICH,I)) .LT. EY) THEN
2684                DIST = ((PHI1-PHM(ICH,I))/EPHI)**2 +
2685      &                ((ALAM1-ALM(ICH,I))/ELAM)**2 +
2686      &                ((X1-XM(ICH,I))/EX)**2 +
2687      &                ((Y1-YM(ICH,I))/EY)**2
2688                IF (DIST .LT. DISTMIN) THEN
2689                   DISTMIN = DIST
2690                   IFIND = I
2691                ENDIF
2692             ENDIF
2693          ENDIF
2694       ENDDO
2695        
2696       RETURN
2697       END
2698 ***********************************************************************************
2699       SUBROUTINE DISTMIN2(X1,Y1,X2,Y2,ICH,EX1,EY1,EX2,EY2,IFIND)
2700 ***********************************************************************************
2701 *    Find the nearest hit in station ICH in the (X,Y) space
2702 *
2703 ***********************************************************************************
2704       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2705       
2706       PARAMETER (MAXHITCH=10000,NBSTATION=5)
2707       
2708       COMMON/CHHIT/XM(NBSTATION,MAXHITCH),YM(NBSTATION,MAXHITCH),
2709      &             PHM(NBSTATION,MAXHITCH),ALM(NBSTATION,MAXHITCH),
2710      &             IZM(NBSTATION,MAXHITCH),
2711      &             IP(NBSTATION,MAXHITCH),JHIT(NBSTATION),
2712      &             XMR(NBSTATION,MAXHITCH,2),YMR(NBSTATION,MAXHITCH,2) 
2713
2714       IFIND = 0
2715       DISTMIN=2.
2716       DO I=1,JHIT(ICH)
2717          IF (IZM(ICH,I).EQ.1) THEN
2718             X = X1
2719             Y = Y1
2720          ELSE
2721             X = X2
2722             Y = Y2
2723          ENDIF    
2724          EX = EX1
2725          EY = EY1
2726          IF (ICH.EQ.4.OR.ICH.EQ.5) THEN
2727             IF (IZM(ICH,I).EQ.1) THEN
2728                EX = EX1
2729                EY = EY1
2730             ELSE
2731                EX = EX2
2732                EY = EY2
2733             ENDIF    
2734          ENDIF    
2735          IF (ABS(X-XM(ICH,I)) .LT. EX .AND.
2736      &      ABS(Y-YM(ICH,I)) .LT. EY) THEN
2737             DIST = ((X-XM(ICH,I))/EX)**2 +
2738      &             ((Y-YM(ICH,I))/EY)**2
2739             IF (DIST .LT. DISTMIN) THEN
2740                DISTMIN = DIST
2741                IFIND = I
2742             ENDIF
2743          ENDIF
2744       ENDDO
2745        
2746       RETURN
2747       END
2748 ********************************************************************************
2749       SUBROUTINE H_ACCEPTANCE(ICH)
2750 ********************************************************************************
2751 * Etude de l'acceptance des resonnances en  fonction du H 
2752 * dans la station ICH
2753 *
2754 *  INPUT :    ICH
2755 *  OUTPUT :  Histo #1
2756 ********************************************************************************
2757
2758       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2759  
2760       PARAMETER (MAXHITCH=10000,MAXTRK=50000,MAXHITTOT=20000,
2761      &           NBSTATION=5)
2762       
2763       COMMON/CHHIT/XM(NBSTATION,MAXHITCH),YM(NBSTATION,MAXHITCH),
2764      &             PHM(NBSTATION,MAXHITCH),ALM(NBSTATION,MAXHITCH),
2765      &             IZM(NBSTATION,MAXHITCH),
2766      &             IP(NBSTATION,MAXHITCH),JHIT(NBSTATION),
2767      &             XMR(NBSTATION,MAXHITCH,2),YMR(NBSTATION,MAXHITCH,2) 
2768      
2769       COMMON/RHIT/ITYP(MAXHITTOT),XTR(MAXHITTOT),YTR(MAXHITTOT),
2770      &            PTOT(MAXHITTOT),ID(MAXHITTOT),IZST(MAXHITTOT),
2771      &            PVERT1(MAXHITTOT),PVERT2(MAXHITTOT),PVERT3(MAXHITTOT),
2772      &            ZVERT(MAXHITTOT),NHITTOT
2773
2774  
2775       COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
2776       
2777       COMMON/VERIFGEANT/ITTROUGH(MAXTRK,NBSTATION),
2778      &                  IT_LIST(MAXTRK),IT_NP(MAXTRK),ITCHECK(MAXTRK),
2779      &                  ITRACK(MAXHITTOT)
2780      
2781       COMMON/HCHHIT/HHIT(MAXHITCH),INDEXTAB(MAXHITCH),INDEXMAX
2782
2783       REAL*4 R1,R2
2784       DATA R1,R2/0.,1./
2785  
2786       NMUONI = 0
2787       DO J = 1,JHIT(ICH)
2788          IF (ITYP(IP(ICH,J)).EQ.5.OR.ITYP(IP(ICH,J)).EQ.6) THEN
2789             ISTAK = ID(IP(ICH,J))
2790             ISTAK = MOD(ISTAK,30000)
2791             ISTAK = MOD(ISTAK,10000)
2792             IF (ISTAK.EQ.0) THEN
2793 **              PRINT *,'ACCEPT. id du muon dans st.',ICH,':',ITYP(IP(ICH,J))
2794                NMUONI = NMUONI+1
2795             ENDIF
2796          ENDIF
2797       ENDDO
2798 *      PRINT *,'ACCEPT. nb de muons/res total dans st.',ICH,':',NMUONI
2799 *      pause
2800
2801       DO IH = 1,500
2802          HCUT = IH
2803 *   Sort hits in st. z
2804          CALL ORDONNE_HIT(ICH,HCUT)
2805          NMUON = 0
2806          DO IND = 1,INDEXMAX
2807             IIND = IP(ICH,INDEXTAB(IND))
2808             IDPART = ITYP(IIND)
2809             ISTAK = ID(IIND)
2810             ISTAK = MOD(ISTAK,30000)
2811             ISTAK = MOD(ISTAK,10000)
2812 **            PRINT *,' IDPART=',IDPART,' ISTAK=',ISTAK
2813             IF (IDPART.EQ.5.OR.IDPART.EQ.6.AND.ISTAK.EQ.0) THEN
2814                NMUON = NMUON+1
2815             ENDIF
2816          ENDDO
2817          IF (NMUON.EQ.2.AND.NMUONI.EQ.2) THEN
2818             CALL CHFILL(ICH*100,SNGL(HCUT),R1,R2)
2819          ENDIF
2820       ENDDO
2821  
2822       RETURN
2823       END
2824
2825 ********************************************************************************
2826        SUBROUTINE OLDFOLLOW(ZSTR,PEST,IFLAG,XPL,YPL,PHPL,ALPL)
2827 ********************************************************************************
2828 *   Calculate the particle trajectory in the spectrometer and
2829 *   (XPL,YPL,PHPL,ALPL)
2830 *   for the 5 stations.
2831 *       
2832 ********************************************************************************
2833        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2834 *
2835        PARAMETER(NBSTATION=5) 
2836 *
2837        DIMENSION XPL(NBSTATION,2),YPL(NBSTATION,2),PHPL(NBSTATION),
2838      &           ALPL(NBSTATION),PEST(NBSTATION)
2839
2840        COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
2841        
2842        COMMON /MEASUR/XMES(NBSTATION),YMES(NBSTATION),IZMES(NBSTATION),
2843      &                PHMES(NBSTATION),ALMES(NBSTATION),MPOS(NBSTATION),
2844      &                MANG(NBSTATION)
2845        COMMON /MAGNET/BL3,B0
2846
2847        LOGICAL LFLAG, LFLAG1
2848
2849        XSTR = PEST(4)
2850        YSTR = PEST(5)
2851        PXZINV = PEST(1)
2852        TPHI = PEST(2)
2853        PHI = ATAN(TPHI)
2854        TALAM = PEST(3)
2855        PXZ = 1.0/PXZINV
2856        PY = ABS(PXZ)*TALAM
2857        PX = -ABS(PXZ)*SIN(PHI)
2858        PZ = -ABS(PXZ)*COS(PHI)
2859        PXY = SQRT(PX**2 + PY**2)
2860        FI=ATAN2(PY,PX)
2861        SINFI = SIN(FI)
2862        COSFI = COS(FI)
2863        TTHET = PZ/PXY
2864        RS = PXY*(100.0/(0.299792458*BL3))
2865        IF(PXZINV.LT.0.0) RS = -RS
2866 *       XC = XSTR + RS*SIN(FI)
2867 *       YC = YSTR - RS*COS(FI)
2868        PX0 = PX
2869        PY0 = PY
2870        LFLAG = .TRUE.
2871        LFLAG1 = .TRUE.
2872 *       PRINT *, XC,YC,RS,FI,TTHET,PXY,PZ
2873
2874        DO J = 1,5
2875
2876           IF (IFLAG.EQ.3 .OR. MPOS(J).EQ.1) THEN
2877              IF(ZPLANE(J) .GT. ZCOIL) THEN
2878 *                DFI =  (ZPLANE(J)-ZSTR)/(TTHET*RS)
2879 *                FIN = FI - DFI
2880 *                XPL(J,1) = XC - RS*SIN(FIN)
2881 *                YPL(J,1) = YC + RS*COS(FIN)
2882                  DFR =  (ZPLANE(J)-ZSTR)/TTHET
2883                  XPL(J,1) = XSTR + DFR*COSFI + 0.5D0*DFR*DFR*SINFI/RS
2884                  YPL(J,1) = YSTR + DFR*SINFI - 0.5D0*DFR*DFR*COSFI/RS
2885                  DFR2 =  (ZPLANE(J)-DZ_PL(J)-ZSTR)/TTHET
2886                  XPL(J,2) = XSTR + DFR2*COSFI + 0.5D0*DFR2*DFR2*SINFI/RS
2887                  YPL(J,2) = YSTR + DFR2*SINFI - 0.5D0*DFR2*DFR2*COSFI/RS
2888                 IF (IFLAG.EQ.3 .OR. MANG(J).EQ.1) THEN
2889 *                   PX=PXY*COS(FIN)
2890 *                   PY=PXY*SIN(FIN)
2891                    PX = PX0 + DFR * (PY0 - 0.5D0*PX0*DFR/RS) / RS
2892                    PY = PY0 - DFR * (PX0 + 0.5D0*PY0*DFR/RS) / RS
2893                    PHPL(J)=ATAN(PX/PZ)
2894                    ALPL(J)=ATAN(PY/SQRT(PX**2+PZ**2))
2895                 ENDIF
2896               ELSE 
2897                 IF( LFLAG) THEN
2898 *                   DFI =  (ZCOIL-ZSTR)/(TTHET*RS)
2899 *                   FIN = FI - DFI
2900 *                   XCOIL = XC - RS*SIN(FIN)
2901 *                   YCOIL = YC + RS*COS(FIN)
2902                    DFR =  (ZCOIL-ZSTR)/TTHET
2903                    XCOIL = XSTR + DFR*COSFI + 0.5D0*DFR*DFR*SINFI/RS
2904                    YCOIL = YSTR + DFR*SINFI - 0.5D0*DFR*DFR*COSFI/RS
2905 *                   PX=PXY*COS(FIN)
2906 *                   PY=PXY*SIN(FIN)
2907                    PX = PX0 + DFR * (PY0 - 0.5D0*PX0*DFR/RS) / RS
2908                    PY = PY0 - DFR * (PX0 + 0.5D0*PY0*DFR/RS) / RS
2909                    PXZ = SQRT(PX**2 + PZ**2)
2910                    PHI=ATAN(PX/PZ)
2911                    TALAM = PY/PXZ
2912                    ALAM = ATAN(TALAM)
2913                    RD = PXZ*(100.0/(0.299792458*B0))
2914                    IF(PXZINV.LT.0.0) RD = -RD
2915                    ZC = ZCOIL - RD*SIN(PHI)
2916                    XC = XCOIL + RD*COS(PHI)
2917                    IF(ABS(ZMAGEND-ZC).GT.ABS(RD)) STOP 'FOLLOW'
2918                    LFLAG = .FALSE.
2919                 ENDIF
2920                 IF(ZPLANE(J) .GT. ZMAGEND) THEN 
2921                   FIN = ASIN((ZPLANE(J) - ZC)/RD)
2922                   XPL(J,1)= XC - RD*COS(FIN)
2923                   YPL(J,1)= YCOIL - RD*(FIN - PHI)*TALAM
2924                   FIN2 = ASIN((ZPLANE(J)-DZ_PL(J) - ZC)/RD)
2925                   XPL(J,2)= XC - RD*COS(FIN2)
2926                   YPL(J,2)= YCOIL - RD*(FIN2 - PHI)*TALAM
2927                   PHPL(J)=FIN
2928                   ALPL(J)=ALAM
2929                 ELSE
2930                   IF (LFLAG1) THEN
2931                     FIN = ASIN((ZMAGEND - ZC)/RD)
2932                     XMAGEND = XC - RD*COS(FIN)
2933                     YMAGEND = YCOIL - RD*(FIN - PHI)*TALAM
2934                     TPHI = TAN(FIN)
2935                     CPHI = COS(FIN)
2936                     LFLAG1 = .FALSE.
2937                   ENDIF
2938                   XPL(J,1) = XMAGEND + (ZPLANE(J)-ZMAGEND)*TPHI
2939                   YPL(J,1) = YMAGEND - (ZPLANE(J)-ZMAGEND)*TALAM/CPHI
2940                   XPL(J,2) = XMAGEND + (ZPLANE(J)-DZ_PL(J)-ZMAGEND)*TPHI
2941                   YPL(J,2) = YMAGEND - (ZPLANE(J)-DZ_PL(J)-ZMAGEND)*
2942      &                      TALAM/CPHI
2943                   PHPL(J)=FIN
2944                   ALPL(J)=ALAM
2945                 ENDIF
2946              ENDIF
2947           ENDIF
2948        ENDDO
2949        RETURN
2950        END
2951 ********************************************************************************
2952        SUBROUTINE FOLLOW(ZSTR,PEST,IFLAG,XPL,YPL,PHPL,ALPL)
2953 ********************************************************************************
2954 *   Calculate the particle trajectory in the spectrometer 
2955 *   (XPL,YPL,PHPL,ALPL)
2956 *   for the 5 stations.
2957 *       Runge Kutta
2958 ********************************************************************************
2959        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2960 *
2961        PARAMETER(NBSTATION=5,NPLANE=10) 
2962 *
2963        DIMENSION XPL(NBSTATION,2),YPL(NBSTATION,2),PHPL(NBSTATION),
2964      &           ALPL(NBSTATION),PEST(NBSTATION)
2965
2966        COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
2967
2968        COMMON/PARAM/ZPLANEP(NPLANE),THICK,XPREC,YPREC,B0,BL3,ZMAGS,
2969      &             ZMAGE,ZABS,XMAG,ZBP1,ZBP2,CONST
2970        
2971        COMMON /MEASUR/XMES(NBSTATION),YMES(NBSTATION),IZMES(NBSTATION),
2972      &                PHMES(NBSTATION),ALMES(NBSTATION),MPOS(NBSTATION),
2973      &                MANG(NBSTATION)
2974  
2975        
2976        DIMENSION VECT(7),VOUT(7)
2977        
2978        STEP = 6.                ! 1 cm
2979        NSTEPMAX = 5000
2980                 
2981        ASIGN = 1.
2982        IF (PEST(1).LT.0.) ASIGN = -1.
2983        TPHI = -1.*PEST(2)
2984        PHI = DATAN(TPHI)
2985        TALAM = PEST(3)
2986        ALAM = DATAN(TALAM)
2987        PXZ = DABS(1./PEST(1))
2988        
2989        PX = PXZ*DSIN(PHI)
2990        PY = PXZ*DTAN(ALAM)
2991        PZ = PXZ*DCOS(PHI)
2992        PTOT = PXZ/DCOS(ALAM)
2993        
2994        VECT(1) = PEST(4)
2995        VECT(2) = PEST(5)
2996        VECT(3) = 0.
2997        VECT(4) = PX/PTOT
2998        VECT(5) = PY/PTOT
2999        VECT(6) = PZ/PTOT
3000        VECT(7) = PTOT
3001               
3002        Z = VECT(3)
3003        NSTEP = 0
3004 *
3005 ** Runge Kutta  
3006 **      PRINT *,' AV GRKUTA ASIGN',ASIGN,' THET',THET
3007        ISTOLD = 0 
3008        DO ICH = 1,NPLANE
3009
3010           IST = INT(FLOAT(ICH+1)/2.)
3011
3012
3013           DO WHILE (Z.GE.0..AND.Z.LT.ABS(ZPLANEP(ICH))
3014      &         .AND.NSTEP.LE.NSTEPMAX)
3015              NSTEP = NSTEP+1 
3016 **          WRITE(6,*) NSTEP,(VECT(I),I=1,7)
3017              CALL RECO_GRKUTA (ASIGN,STEP,VECT,VOUT)
3018              DO I = 1,7
3019                 VECT(I) = VOUT(I)
3020              ENDDO   
3021              Z = VECT(3)
3022           ENDDO
3023           IF (IST.NE.ISTOLD) THEN
3024              IPCH = 1
3025           ELSE
3026              IPCH = 2
3027           ENDIF   
3028           XPL(IST,IPCH) = VECT(1)-(Z-ABS(ZPLANEP(ICH)))*VECT(4)/VECT(6)
3029           YPL(IST,IPCH) = VECT(2)-(Z-ABS(ZPLANEP(ICH)))*VECT(5)/VECT(6)
3030           IF (IPCH.EQ.2) THEN
3031              DX = XPL(IST,2)-XPL(IST,1)   
3032              DY = YPL(IST,2)-YPL(IST,1)
3033              PHPL(IST) = -1.*DATAN2(DX,DZ_PL(IST))
3034              ALPL(IST) = DATAN2(DY,DSQRT(DX**2+DZ_PL(IST)**2))
3035           ENDIF   
3036           ISTOLD = IST
3037        ENDDO   
3038 **          print *,' vect= ',vect(1),vect(2),vect(3),vect(4),vect(5),
3039 **     &        vect(6),vect(7)
3040
3041         
3042        RETURN
3043        END
3044 *******************************************************************************
3045        SUBROUTINE FCN(NPAR,GRAD,FVAL,PEST,IFLAG,FUTIL)
3046 *******************************************************************************
3047 *   Calculate FVAL=CHI2 the function minimized by minuit for a given track
3048 *
3049 *******************************************************************************
3050        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
3051
3052        PARAMETER(NBSTATION=5)
3053
3054 *       DIMENSION PEST(*),GRAD(*)
3055        DIMENSION PEST(5),GRAD(5)
3056        DIMENSION PEEST(NBSTATION)
3057
3058        COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
3059        
3060      
3061        COMMON/PRECIS/EEXM(NBSTATION),EEYM(NBSTATION),EEPH(NBSTATION),
3062      &     EEAL(NBSTATION)
3063
3064        COMMON /MEASUR/XMES(NBSTATION),YMES(NBSTATION),IZMES(NBSTATION),
3065      &               PHMES(NBSTATION),ALMES(NBSTATION), MPOS(NBSTATION),
3066      &               MANG(NBSTATION)
3067
3068        COMMON /PLANE/XPL(NBSTATION,2),YPL(NBSTATION,2),PHPL(NBSTATION),
3069      &               ALPL(NBSTATION),CHI2PL
3070
3071        COMMON/VERTEX/ERRV,IVERTEX
3072
3073        PEEST(1)=PEST(1)
3074        PEEST(2)=PEST(2)
3075        PEEST(3)=PEST(3)
3076        IF(IVERTEX.EQ.1) THEN
3077           PEEST(4)=PEST(4)      ! position du vertex
3078           PEEST(5)=PEST(5)
3079        ELSE
3080           PEEST(4)=0.0D0
3081           PEEST(5)=0.0D0
3082        ENDIF
3083
3084        CALL FOLLOW (0.0D0,PEEST,IFLAG,XPL,YPL,PHPL,ALPL) ! calcul des 
3085        IF(IFLAG.EQ.1) THEN                       ! points d impacts dans les
3086           PRINT *,'FCN ',XPL(4,1),XMES(4)        ! plans
3087           PRINT *,'FCN ',YPL(4,1),YMES(4)
3088           PRINT *,'FCN ',XPL(5,1),XMES(5)
3089           PRINT *,'FCN ',YPL(5,1),YMES(5)
3090        ENDIF
3091
3092 *       IF (IVERTEX.EQ.1) THEN
3093 *          FVAL = (PEST(4)/ERRV)**2 + (PEST(5)/ERRV)**2
3094 *       ELSE
3095           FVAL = 0.0D0
3096 *       ENDIF
3097        NPLPL = 0 
3098        DO J = 1,NBSTATION           
3099           IF (MPOS(J).EQ.1) THEN
3100              NPLPL = NPLPL+1
3101              XPLC = XPL(J,IZMES(J)) 
3102              YPLC = YPL(J,IZMES(J))
3103              FF = (XMES(J) - XPLC)/EEXM(J)
3104              FVAL =FVAL + FF**2
3105              FF = (YMES(J) - YPLC)/EEYM(J)
3106              FVAL =FVAL + FF**2
3107           ENDIF
3108           IF (MANG(J).EQ.1) THEN
3109              NPLPL = NPLPL+1
3110              FF = (PHMES(J) - PHPL(J))/EEPH(J)
3111              FVAL =FVAL + FF**2
3112              FF = (ALMES(J) - ALPL(J))/EEAL(J)
3113              FVAL =FVAL + FF**2
3114           ENDIF
3115        ENDDO
3116 **       PRINT *,'ST 1',XPL(1,1),XMES(1),YPL(1,1),YMES(1)          
3117 **       PRINT *,'ST 2',XPL(2,1),XMES(2),YPL(2,1),YMES(2)          
3118 **       PRINT *,'ST 3',XPL(3,1),XMES(3),YPL(3,1),YMES(3)          
3119 **       PRINT *,'ST 4',XPL(4,1),XMES(4),YPL(4,1),YMES(4)          
3120 **       PRINT *,'ST 5',XPL(5,1),XMES(5),YPL(5,1),YMES(5)          
3121        NPARAM = 3
3122        IF (IVERTEX.EQ.1) NPARAM = 5
3123        CHI2PL = FVAL/FLOAT(2*NPLPL-NPARAM)       
3124
3125        RETURN
3126        END
3127 ********************************************************************************
3128       SUBROUTINE STOCK_CANDIDAT(ICH1,JHITCH1,ICH2,IFIND,EXM,EYM,EPH,EAL,
3129      &                          NCAN,ICODE)
3130 ********************************************************************************
3131 *   Fill common CANDIDAT with track candidates from the search in stations 4&5
3132 *
3133 ********************************************************************************
3134       
3135       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
3136
3137       PARAMETER (MAXHITTOT=20000,MAXHITCH=10000,NBSTATION=5,MAXCAN=1000)
3138       
3139       COMMON/CHHIT/XM(NBSTATION,MAXHITCH),YM(NBSTATION,MAXHITCH),
3140      &             PHM(NBSTATION,MAXHITCH),ALM(NBSTATION,MAXHITCH),
3141      &             IZM(NBSTATION,MAXHITCH),
3142      &             IP(NBSTATION,MAXHITCH),JHIT(NBSTATION),
3143      &             XMR(NBSTATION,MAXHITCH,2),YMR(NBSTATION,MAXHITCH,2) 
3144       
3145       COMMON/RHIT/ITYP(MAXHITTOT),XTR(MAXHITTOT),YTR(MAXHITTOT),
3146      &            PTOT(MAXHITTOT),ID(MAXHITTOT),IZST(MAXHITTOT),
3147      &            PVERT1(MAXHITTOT),PVERT2(MAXHITTOT),PVERT3(MAXHITTOT),
3148      &            ZVERT(MAXHITTOT),NHITTOT
3149
3150             
3151       COMMON/CANDIDAT/JCAN(NBSTATION,MAXCAN),JCANTYP(MAXCAN),
3152      &          EEX(MAXCAN),EEY(MAXCAN),EEP(MAXCAN),EEA(MAXCAN)
3153       
3154 **      PRINT *,'STOCK st. init.=',ICH1,'id. init.=',ID(IP(ICH1,JHITCH1))
3155 **      PRINT *,'STOCK st. finale=',ICH2,'id. final=',ID(IP(ICH2,IFIND))
3156 **      PRINT *,'STOCK ifind=',IFIND 
3157 **      PRINT *,'STOCK icode=',ICODE 
3158       NCAN = NCAN+1
3159       JCAN(ICH1,NCAN) = JHITCH1
3160       JCAN(ICH2,NCAN) = IFIND
3161       JCANTYP(NCAN) = ICODE 
3162       EEX(NCAN) = EXM
3163       EEY(NCAN) = EYM
3164       EEP(NCAN) = EPH
3165       EEA(NCAN) = EAL
3166
3167       RETURN
3168       END
3169 *****************************************************************************
3170       SUBROUTINE CHECK_HISTO4(IDHIST,ICH2,IHIT2,ICH1,IHIT1,
3171      &                        X1,Y1,PHI1,ALAM1,P1,EXM,EYM,EPH,EAL)   
3172 *****************************************************************************
3173 *   Check hit IHIT2 with GEANT informations from hit HIT1
3174 *
3175 *    INPUT : ICH2 : No st. de recherche
3176 *            IDCH1,X1,Y1,PHI1,ALAM1,P1 : Trace de reference
3177 *    OUTPUT : JOK : No hit dans ICH2 appartenant a la meme trace.
3178
3179 *****************************************************************************
3180       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
3181       
3182       PARAMETER (MAXHITCH=10000,MAXHITTOT=20000,NBSTATION=5,
3183      &           MAXTRK=50000)
3184       
3185       COMMON/CHHIT/XM(NBSTATION,MAXHITCH),YM(NBSTATION,MAXHITCH),
3186      &             PHM(NBSTATION,MAXHITCH),ALM(NBSTATION,MAXHITCH),
3187      &             IZM(NBSTATION,MAXHITCH),
3188      &             IP(NBSTATION,MAXHITCH),JHIT(NBSTATION),
3189      &             XMR(NBSTATION,MAXHITCH,2),YMR(NBSTATION,MAXHITCH,2) 
3190
3191       COMMON/RHIT/ITYP(MAXHITTOT),XTR(MAXHITTOT),YTR(MAXHITTOT),
3192      &            PTOT(MAXHITTOT),ID(MAXHITTOT),IZST(MAXHITTOT),
3193      &            PVERT1(MAXHITTOT),PVERT2(MAXHITTOT),PVERT3(MAXHITTOT),
3194      &            ZVERT(MAXHITTOT),NHITTOT
3195
3196       COMMON/VERIFGEANT/ITTROUGH(MAXTRK,NBSTATION),
3197      &                  IT_LIST(MAXTRK),IT_NP(MAXTRK),ITCHECK(MAXTRK),
3198      &                  ITRACK(MAXHITTOT)
3199
3200       COMMON/DEBEVT/IDEBUG
3201
3202       REAL*4 R2
3203       DATA R2/1./
3204
3205       JOK = 0
3206        
3207       DO I=1,JHIT(ICH2)
3208          IF (PHM(ICH2,I).LE.6.3) THEN
3209             IF (ID(IP(ICH1,IHIT1)).EQ.ID(IP(ICH2,I))) THEN
3210                JOK = I 
3211                IF (IDHIST.GT.0) THEN
3212                   CALL CHFILL2(IDHIST,SNGL(P1),
3213      &                        SNGL((X1-XM(ICH2,I))**2),1.)
3214                   CALL CHFILL2(IDHIST+1,SNGL(P1),
3215      &                        SNGL((Y1-YM(ICH2,I))**2),1.)
3216                   CALL CHFILL2(IDHIST+2,SNGL(P1),
3217      &                     SNGL((PHI1-PHM(ICH2,I))**2),1.)
3218                   CALL CHFILL2(IDHIST+3,SNGL(P1),
3219      &                     SNGL((ALAM1-ALM(ICH2,I))**2),1.)
3220                   CALL CHF1(400+IDHIST,
3221      &                        SNGL(X1-XM(ICH2,I)),1.)
3222                   CALL CHF1(400+IDHIST+1,
3223      &                        SNGL(Y1-YM(ICH2,I)),1.)
3224                   CALL CHF1(400+IDHIST+2,
3225      &                     SNGL(PHI1-PHM(ICH2,I)),1.)
3226                   CALL CHF1(400+IDHIST+3,
3227      &                     SNGL(ALAM1-ALM(ICH2,I)),1.)
3228                   CALL CHF1(IDHIST+4,SNGL(P1),R2)
3229                ENDIF   
3230             ENDIF
3231          ENDIF
3232       ENDDO
3233       
3234       IF (JOK.GT.0) THEN
3235          IF (ITCHECK(ITRACK(IP(ICH1,IHIT1))).EQ.1) THEN
3236             IF (IDEBUG.GE.2) THEN 
3237                IF (IHIT2.EQ.0) THEN
3238                   PRINT *,'CHECK4 histo nb:',IDHIST 
3239                   PRINT *,'CHECK4 p de st.',ICH1,'=',P1
3240                   PRINT *,'CHECK4 track not found in st.',ICH2
3241                   PRINT *,'CHECK4 error X :',(XM(ICH2,JOK)-X1), EXM
3242                   PRINT *,'CHECK4 error Y :',(YM(ICH2,JOK)-Y1), EYM
3243                   PRINT *,'CHECK4 error PHI :',(PHM(ICH2,JOK)-PHI1),EPH
3244                   PRINT *,'CHECK4 error ALAM :',(ALM(ICH2,JOK)-ALAM1),
3245      &                     EAL
3246                ELSEIF (IHIT2.NE.JOK) THEN
3247                   PRINT *,'CHECK4 histo nb:',IDHIST 
3248                   PRINT *,'CHECK4 p de st.',ICH1,'=',P1
3249                   PRINT *,'CHECK4 ghost in st.',ICH2
3250                   PRINT *,'CHECK4 id part. recherchee:',
3251      &                    ID(IP(ICH1,IHIT1))
3252                   PRINT *,'CHECK4 id ghost trouve    :',
3253      &                    ID(IP(ICH2,IHIT2))  
3254                ENDIF
3255             ENDIF
3256          ENDIF
3257       ENDIF
3258                
3259       RETURN
3260       END     
3261 *****************************************************************************
3262       SUBROUTINE CHECK_HISTO2(IDHIST,ICH2,IHIT2,ICH1,IHIT1,
3263      &                        X1,Y1,X2,Y2,P1,EX1,EY1,EX2,EY2)   
3264 *****************************************************************************
3265 *   Check hit IHIT2 with GEANT informations from hit HIT1
3266 *
3267 *    INPUT : IDHIST : No histo
3268 *            ICH2 : No st. de recherche
3269 *            IDCH1,X1,Y1,PHI1,ALAM1,P1 : Trace de reference
3270 *    OUTPUT : JOK : No hit dans ICH2 appartenant a la meme trace.
3271
3272 *****************************************************************************
3273       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
3274       
3275       PARAMETER (MAXHITCH=10000,MAXHITTOT=20000,NBSTATION=5,
3276      &           MAXTRK=50000)
3277       
3278       COMMON/CHHIT/XM(NBSTATION,MAXHITCH),YM(NBSTATION,MAXHITCH),
3279      &             PHM(NBSTATION,MAXHITCH),ALM(NBSTATION,MAXHITCH),
3280      &             IZM(NBSTATION,MAXHITCH),
3281      &             IP(NBSTATION,MAXHITCH),JHIT(NBSTATION),
3282      &             XMR(NBSTATION,MAXHITCH,2),YMR(NBSTATION,MAXHITCH,2) 
3283
3284       COMMON/RHIT/ITYP(MAXHITTOT),XTR(MAXHITTOT),YTR(MAXHITTOT),
3285      &            PTOT(MAXHITTOT),ID(MAXHITTOT),IZST(MAXHITTOT),
3286      &            PVERT1(MAXHITTOT),PVERT2(MAXHITTOT),PVERT3(MAXHITTOT),
3287      &            ZVERT(MAXHITTOT),NHITTOT
3288      
3289       COMMON/VERIFGEANT/ITTROUGH(MAXTRK,NBSTATION),
3290      &                  IT_LIST(MAXTRK),IT_NP(MAXTRK),ITCHECK(MAXTRK),
3291      &                  ITRACK(MAXHITTOT)
3292
3293       COMMON/DEBEVT/IDEBUG
3294  
3295       REAL*4 R2
3296       DATA R2/1./
3297
3298       JOK = 0
3299        
3300       DO I=1,JHIT(ICH2)
3301          IF (ID(IP(ICH1,IHIT1)).EQ.ID(IP(ICH2,I))) THEN
3302             JOK = I 
3303             IF (IDHIST.GT.0) THEN
3304                IF (IZM(ICH2,I).EQ.1) THEN
3305                   X = X1
3306                   Y = Y1
3307                ELSE
3308                   X = X2
3309                   Y = Y2
3310                ENDIF 
3311                CALL CHFILL2(200+IDHIST,SNGL(P1),
3312      &                       SNGL((X-XM(ICH2,I))**2),1.)
3313                CALL CHFILL2(200+IDHIST+1,SNGL(P1),
3314      &                       SNGL((Y-YM(ICH2,I))**2),1.)
3315                CALL CHF1(200+IDHIST+2,
3316      &                       SNGL(X-XM(ICH2,I)),1.)
3317                CALL CHF1(200+IDHIST+3,
3318      &                       SNGL(Y-YM(ICH2,I)),1.)
3319                CALL CHF1(200+IDHIST+4,SNGL(P1),R2)
3320             ENDIF   
3321          ENDIF
3322       ENDDO
3323       
3324       IF (JOK.GT.0) THEN
3325          IF (ITCHECK(ITRACK(IP(ICH1,IHIT1))).EQ.1) THEN
3326             EXM = EX1
3327             EYM = EY1 
3328             IF (IZM(ICH2,JOK).EQ.1) THEN
3329                X = X1
3330                Y = Y1
3331                IF (ICH2.EQ.4.OR.ICH2.EQ.5) THEN
3332                   EXM = EX1
3333                   EYM = EY1
3334                ENDIF   
3335             ELSE
3336                X = X2
3337                Y = Y2
3338                IF (ICH2.EQ.4.OR.ICH2.EQ.5) THEN
3339                   EXM = EX2
3340                   EYM = EY2
3341                ENDIF   
3342             ENDIF      
3343             IF (IDEBUG.GE.2) THEN 
3344                IF (IHIT2.EQ.0) THEN
3345                   PRINT *,'CHECK2 histo nb:',IDHIST 
3346                   PRINT *,'CHECK2 p de st.',ICH1,'=',P1
3347                   PRINT *,'CHECK2 track not found in st.',ICH2
3348                   PRINT *,'CHECK2 error X :',(XM(ICH2,JOK)-X), EXM
3349                   PRINT *,'CHECK2 error Y :',(YM(ICH2,JOK)-Y), EYM
3350                ELSEIF(IHIT2.NE.JOK) THEN
3351                   PRINT *,'CHECK2 histo nb:',IDHIST 
3352                   PRINT *,'CHECK2 p de st.',ICH1,'=',P1
3353                   PRINT *,'CHECK2 ghost in st.',ICH2
3354                   PRINT *,'CHECK2 id part. recherchee:',
3355      &                    ID(IP(ICH1,IHIT1))
3356                   PRINT *,'CHECK2 id ghost trouve    :',
3357      &                    ID(IP(ICH2,IHIT2))  
3358                ENDIF
3359             ENDIF
3360          ENDIF
3361       ENDIF
3362                
3363       RETURN
3364       END
3365
3366
3367
3368 ************************************************************************        
3369       DOUBLE PRECISION FUNCTION DEDX (P,THET,XEA,YEA)
3370 ************************************************************************        
3371 *    DEDX est la nouvelle impulsion au vertex, corrigee de la perte    
3372 *    d'energie dans l'absorbeur
3373 *
3374 ************************************************************************        
3375       
3376       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3377
3378       REA = DSQRT(XEA**2+YEA**2)
3379       IF (REA.GT.26.3611) THEN
3380 *   Plomb      
3381          SPB = 5./DCOS(THET)
3382          P = P+SPB/1000.*11.35*(16.66D-3 * P+1.33) 
3383 *   Polyethylene      
3384          SPO = 5./DCOS(THET)
3385          P = P+SPO/1000.*0.935*(2.22D-3 * P+2.17) 
3386 *   Plomb      
3387          SPB = 5./DCOS(THET)
3388          P = P+SPB/1000.*11.35*(16.66D-3 * P+1.33) 
3389 *   Polyethylene      
3390          SPO = 5./DCOS(THET)
3391          P = P+SPO/1000.*0.935*(2.22D-3 * P+2.17) 
3392 *   Plomb      
3393          SPB = 5./DCOS(THET)
3394          P = P+SPB/1000.*11.35*(16.66D-3 * P+1.33) 
3395 *   Polyethylene      
3396          SPO = 5./DCOS(THET)
3397          P = P+SPO/1000.*0.935*(2.22D-3 * P+2.17) 
3398 *   Plomb      
3399          SPB = 5./DCOS(THET)
3400          P = P+SPB/1000.*11.35*(16.66D-3 * P+1.33) 
3401
3402       ELSE   
3403 *   Tungstene     
3404          SW = (503.-467.)/DCOS(THET) 
3405          P = P+SW/1000.*19.3*(16.66D-3 * P+1.33)
3406       ENDIF
3407
3408 *   Concrete      
3409       SCONC = (467.-315.)/DCOS(THET)
3410       P = P+SCONC/1000.*2.5*(2.22D-3*P+2.17)
3411       
3412 *   Carbone      
3413       SC = (315.-90.)/DCOS(THET)
3414       P = P+SC/1000.*1.93*(2.22D-3*P+2.17) ! Carbone
3415                
3416       DEDX = P  
3417       
3418       RETURN
3419  
3420       END
3421 ************************************************************************        
3422       SUBROUTINE BRANSON(PXZ,PHI,ALAM,XEA,YEA)
3423 ************************************************************************        
3424 *
3425 *   Correction de Branson du multiple scattering dans l'absorbeur
3426 *
3427 ************************************************************************
3428
3429       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3430       
3431       PARAMETER(NPLANE=10) 
3432       COMMON/PARAM/ZPLANEP(NPLANE),THICK,XPREC,YPREC,B0,BL3,ZMAGS,
3433      &             ZMAGE,ZABS,XMAG,ZBP1,ZBP2,CONST
3434      
3435       ASIGN = 1.      
3436       IF (PXZ.LT.0.) ASIGN = -1.       
3437       PXZ = DABS(PXZ)
3438       PX = PXZ*DSIN(PHI)
3439       PY = PXZ*DTAN(ALAM)
3440       PZ = PXZ*DCOS(PHI)
3441       
3442       PTOT = PXZ/DCOS(ALAM)
3443       
3444       ZEA = ZABS
3445       
3446       REA = DSQRT(XEA**2+YEA**2)
3447       IF (REA.GT.26.3611) THEN 
3448          ZBP = ZBP1
3449       ELSE            ! Abso. W  pour theta < 3 deg
3450          ZBP = ZBP2
3451       ENDIF
3452              
3453       XBP = XEA-PX/PZ*(ZEA-ZBP)
3454       YBP = YEA-PY/PZ*(ZEA-ZBP)
3455       PZ = PTOT*ZBP/DSQRT(XBP**2+YBP**2+ZBP**2)
3456       PX = PZ*XBP/ZBP
3457       PY = PZ*YBP/ZBP
3458       PXZ = DSQRT(PX**2+PZ**2)
3459       PHI = DATAN2(PX,PZ)
3460       ALAM = DATAN2(PY,PXZ)
3461
3462 **      THET = DATAN2(REA,ZEA)
3463
3464       PT = DSQRT(PX**2+PY**2)      
3465       THET = DATAN2(PT,PZ) 
3466       PTOT =  DEDX(PTOT,THET,XEA,YEA)
3467       
3468       PXZ = ASIGN*PTOT*DCOS(ALAM)
3469       
3470       RETURN
3471       END
3472
3473 ***************************************************************
3474       SUBROUTINE DSINV(N,A,IDIM,IFAIL)
3475 ***************************************************************
3476
3477       DOUBLE PRECISION    A(IDIM,*),  ZERO,  ONE,  X, Y
3478
3479       REAL                PIVOTF
3480       CHARACTER*6         HNAME
3481
3482       DOUBLE PRECISION    S1, S31, S32, S33,  DOTF
3483
3484       PIVOTF(X)    =  SNGL(X)
3485       DOTF(X,Y,S1)  =  X * Y + S1
3486
3487       DATA      HNAME               /  'DSINV '  /
3488       DATA      ZERO, ONE           /  0.D0, 1.D0 /
3489
3490       IF(IDIM .LT. N  .OR.  N .LE. 0)  GOTO 900
3491 *     
3492 * sfact.inc
3493 *
3494       IFAIL  =  0
3495       DO 144    J  =  1, N
3496          IF(PIVOTF(A(J,J)) .LE. 0.)  GOTO 150
3497          A(J,J)  =  ONE / A(J,J)
3498          IF(J .EQ. N)  GOTO 199
3499  140     JP1  =  J+1
3500          DO 143   L  =  JP1, N
3501             A(J,L)  =  A(J,J)*A(L,J)
3502             S1      =  -A(L,J+1)
3503             DO 141  I  =  1, J
3504                S1  =  DOTF(A(L,I),A(I,J+1),S1)
3505  141        CONTINUE
3506             A(L,J+1)  =  -S1
3507  143     CONTINUE
3508  144  CONTINUE
3509  150  IFAIL  =  -1
3510       RETURN
3511  199  CONTINUE
3512 *     
3513 * sfinv.inc
3514 *     
3515       IF(N .EQ. 1)  GOTO 399
3516       A(1,2)  =  -A(1,2)
3517       A(2,1)  =   A(1,2)*A(2,2)
3518       IF(N .EQ. 2)  GOTO 320
3519       DO 314    J  =  3, N
3520          JM2  =  J - 2
3521          DO 312 K  =  1, JM2
3522             S31  =  A(K,J)
3523             DO 311  I  =  K, JM2
3524                S31  =  DOTF(A(K,I+1),A(I+1,J),S31)
3525  311        CONTINUE
3526             A(K,J)  =  -S31
3527             A(J,K)  =  -S31*A(J,J)
3528  312     CONTINUE
3529          A(J-1,J)  =  -A(J-1,J)
3530          A(J,J-1)  =   A(J-1,J)*A(J,J)
3531  314  CONTINUE
3532  320  J  =  1
3533  323  S33  =  A(J,J)
3534       IF(J .EQ. N)  GOTO 325
3535       JP1  =  J + 1
3536       DO 324 I  =  JP1, N
3537          S33  =  DOTF(A(J,I),A(I,J),S33)
3538  324  CONTINUE
3539  325  A(J,J)  =  S33
3540       JM1  =  J
3541       J    =  JP1
3542       DO 328 K  =  1, JM1
3543          S32  =  ZERO
3544          DO 327  I  =  J, N
3545             S32  =  DOTF(A(K,I),A(I,J),S32)
3546  327     CONTINUE
3547          A(K,J)  =  S32
3548          A(J,K)  =  S32
3549  328  CONTINUE
3550       IF(J .LT. N)  GOTO 323
3551  399  CONTINUE
3552
3553       RETURN
3554  900  CALL TMPRNT(HNAME,N,IDIM,0)
3555       RETURN
3556       END
3557
3558 *******************************************************
3559       SUBROUTINE TMPRNT(NAME,N,IDIM,K)
3560 *******************************************************
3561
3562       CHARACTER*6         NAME
3563       LOGICAL             MFLAG,    RFLAG
3564
3565       IF(NAME(2:2) .EQ. 'S') THEN
3566          CALL KERMTR('F012.1',LGFILE,MFLAG,RFLAG)
3567       ELSE
3568          CALL KERMTR('F011.1',LGFILE,MFLAG,RFLAG)
3569       ENDIF
3570       IF(NAME(3:6) .EQ. 'FEQN') ASSIGN 1002 TO IFMT
3571       IF(NAME(3:6) .NE. 'FEQN') ASSIGN 1001 TO IFMT
3572       IF(MFLAG) THEN
3573          IF(LGFILE .EQ. 0) THEN
3574             IF(NAME(3:6) .EQ. 'FEQN') THEN
3575                WRITE(*,IFMT) NAME, N, IDIM, K
3576             ELSE
3577                WRITE(*,IFMT) NAME, N, IDIM
3578             ENDIF
3579          ELSE
3580             IF(NAME(3:6) .EQ. 'FEQN') THEN
3581                WRITE(LGFILE,IFMT) NAME, N, IDIM, K
3582             ELSE
3583                WRITE(LGFILE,IFMT) NAME, N, IDIM
3584             ENDIF
3585          ENDIF
3586       ENDIF
3587       IF(.NOT. RFLAG) CALL ABEND
3588       RETURN
3589  1001 FORMAT(7X, 31H PARAMETER ERROR IN SUBROUTINE , A6,
3590      +     27H ... (N.LT.1 OR IDIM.LT.N).,
3591      +             5X, 3HN =, I4, 5X, 6HIDIM =, I4, 1H. )
3592  1002 FORMAT(7X, 31H PARAMETER ERROR IN SUBROUTINE , A6,
3593      +     37H ... (N.LT.1 OR IDIM.LT.N OR K.LT.1).,
3594      +     5X, 3HN =, I4, 5X, 6HIDIM =, I4, 5X, 3HK =, I4,1H.)
3595       END
3596 *
3597 * $Id$
3598 *
3599 * $Log$
3600 * Revision 1.1.1.1  1996/02/15 17:48:35  mclareni
3601 * Kernlib
3602 *
3603 *
3604
3605 ***********************************************************
3606       SUBROUTINE KERSET(ERCODE,LGFILE,LIMITM,LIMITR)
3607 ***********************************************************
3608
3609       PARAMETER(KOUNTE  =  27)
3610       CHARACTER*6         ERCODE,   CODE(KOUNTE)
3611       LOGICAL             MFLAG,    RFLAG
3612       INTEGER             KNTM(KOUNTE),       KNTR(KOUNTE)
3613
3614       DATA      LOGF      /  0  /
3615       DATA      CODE(1), KNTM(1), KNTR(1)  / 'C204.1', 255, 255 /
3616       DATA      CODE(2), KNTM(2), KNTR(2)  / 'C204.2', 255, 255 /
3617       DATA      CODE(3), KNTM(3), KNTR(3)  / 'C204.3', 255, 255 /
3618       DATA      CODE(4), KNTM(4), KNTR(4)  / 'C205.1', 255, 255 /
3619       DATA      CODE(5), KNTM(5), KNTR(5)  / 'C205.2', 255, 255 /
3620       DATA      CODE(6), KNTM(6), KNTR(6)  / 'C305.1', 255, 255 /
3621       DATA      CODE(7), KNTM(7), KNTR(7)  / 'C308.1', 255, 255 /
3622       DATA      CODE(8), KNTM(8), KNTR(8)  / 'C312.1', 255, 255 /
3623       DATA      CODE(9), KNTM(9), KNTR(9)  / 'C313.1', 255, 255 /
3624       DATA      CODE(10),KNTM(10),KNTR(10) / 'C336.1', 255, 255 /
3625       DATA      CODE(11),KNTM(11),KNTR(11) / 'C337.1', 255, 255 /
3626       DATA      CODE(12),KNTM(12),KNTR(12) / 'C341.1', 255, 255 /
3627       DATA      CODE(13),KNTM(13),KNTR(13) / 'D103.1', 255, 255 /
3628       DATA      CODE(14),KNTM(14),KNTR(14) / 'D106.1', 255, 255 /
3629       DATA      CODE(15),KNTM(15),KNTR(15) / 'D209.1', 255, 255 /
3630       DATA      CODE(16),KNTM(16),KNTR(16) / 'D509.1', 255, 255 /
3631       DATA      CODE(17),KNTM(17),KNTR(17) / 'E100.1', 255, 255 /
3632       DATA      CODE(18),KNTM(18),KNTR(18) / 'E104.1', 255, 255 /
3633       DATA      CODE(19),KNTM(19),KNTR(19) / 'E105.1', 255, 255 /
3634       DATA      CODE(20),KNTM(20),KNTR(20) / 'E208.1', 255, 255 /
3635       DATA      CODE(21),KNTM(21),KNTR(21) / 'E208.2', 255, 255 /
3636       DATA      CODE(22),KNTM(22),KNTR(22) / 'F010.1', 255,   0 /
3637       DATA      CODE(23),KNTM(23),KNTR(23) / 'F011.1', 255,   0 /
3638       DATA      CODE(24),KNTM(24),KNTR(24) / 'F012.1', 255,   0 /
3639       DATA      CODE(25),KNTM(25),KNTR(25) / 'F406.1', 255,   0 /
3640       DATA      CODE(26),KNTM(26),KNTR(26) / 'G100.1', 255, 255 /
3641       DATA      CODE(27),KNTM(27),KNTR(27) / 'G100.2', 255, 255 /
3642
3643       LOGF  =  LGFILE
3644       L  =  0
3645       IF(ERCODE .NE. ' ')  THEN
3646          DO 10  L = 1, 6
3647             IF(ERCODE(1:L) .EQ. ERCODE)  GOTO 12
3648  10      CONTINUE
3649  12      CONTINUE
3650       ENDIF
3651       DO 14     I  =  1, KOUNTE
3652          IF(L .EQ. 0)  GOTO 13
3653          IF(CODE(I)(1:L) .NE. ERCODE(1:L))  GOTO 14
3654  13      IF(LIMITM.GE.0) KNTM(I)  =  LIMITM
3655          IF(LIMITR.GE.0) KNTR(I)  =  LIMITR
3656  14   CONTINUE
3657       RETURN
3658       ENTRY KERMTR(ERCODE,LOG,MFLAG,RFLAG)
3659       LOG  =  LOGF
3660       DO 20     I  =  1, KOUNTE
3661          IF(ERCODE .EQ. CODE(I))  GOTO 21
3662  20   CONTINUE
3663       WRITE(*,1000)  ERCODE
3664       CALL ABEND
3665       RETURN
3666  21   RFLAG  =  KNTR(I) .GE. 1
3667       IF(RFLAG  .AND.  (KNTR(I) .LT. 255))  KNTR(I)  =  KNTR(I) - 1
3668       MFLAG  =  KNTM(I) .GE. 1
3669       IF(MFLAG  .AND.  (KNTM(I) .LT. 255))  KNTM(I)  =  KNTM(I) - 1
3670       IF(.NOT. RFLAG)  THEN
3671          IF(LOGF .LT. 1)  THEN
3672             WRITE(*,1001)  CODE(I)
3673          ELSE
3674             WRITE(LOGF,1001)  CODE(I)
3675          ENDIF
3676       ENDIF
3677       IF(MFLAG .AND. RFLAG)  THEN
3678          IF(LOGF .LT. 1)  THEN
3679             WRITE(*,1002)  CODE(I)
3680          ELSE
3681             WRITE(LOGF,1002)  CODE(I)
3682          ENDIF
3683       ENDIF
3684       RETURN
3685  1000 FORMAT(' KERNLIB LIBRARY ERROR. ' /
3686      +     ' ERROR CODE ',A6,' NOT RECOGNIZED BY KERMTR',
3687      +     ' ERROR MONITOR. RUN ABORTED.')
3688  1001 FORMAT(/' ***** RUN TERMINATED BY CERN LIBRARY ERROR ',
3689      +     'CONDITION ',A6)
3690  1002 FORMAT(/' ***** CERN LIBRARY ERROR CONDITION ',A6)
3691       END
3692
3693 **************************
3694       subroutine abend
3695 **************************
3696
3697       stop 'abend!'
3698       end
3699
3700 ************************************************************************        
3701       SUBROUTINE FCNFIT(NPAR, GRAD, FVAL, XVAL, IFLAG, FUTIL)
3702 ************************************************************************        
3703 *    With magnetic Field Map GRKUTA
3704 *      
3705 *    Calcule FVAL: la fonction minimisee par MINUIT
3706 *    With magnetic field map
3707 *      
3708 ************************************************************************        
3709       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3710
3711       PARAMETER(NPLANE=10)
3712             
3713       COMMON/PARAM/ZPLANEP(NPLANE),THICK,XPREC,YPREC,B0,BL3,ZMAGS,
3714      &             ZMAGE,ZABS,XMAG,ZBP1,ZBP2,CONST
3715      
3716       COMMON/MEAS/LPLANE(NPLANE),XMP(NPLANE),YMP(NPLANE)
3717      
3718       COMMON/FCNOUT/PXZEA,ALAMEA,PHIEA,XEA,YEA,NPLU,CHI2
3719                  
3720       DIMENSION GRAD(*),XVAL(*),AMS(500),DISTAZ(500)
3721   
3722       DIMENSION XP(NPLANE),YP(NPLANE),
3723      &          COV(NPLANE,NPLANE),AP(NPLANE),COVY(NPLANE,NPLANE)
3724       DIMENSION VECT(7),VOUT(7) 
3725       
3726       STEP = 2.   ! 1 cm
3727       NSTEPMAX = 5000
3728       PITODEG = 57.295
3729       XV = XVAL(4)   
3730       YV = XVAL(5)   
3731               
3732       ASIGN = 1.
3733       IF (XVAL(1).LT.0.) ASIGN = -1.
3734       PHI  = XVAL(2)
3735       ALAM = XVAL(3)
3736       PXZ = DABS(1./XVAL(1))
3737            
3738       PX = PXZ*DSIN(PHI)
3739       PY = PXZ*DTAN(ALAM)
3740       PZ = PXZ*DCOS(PHI)
3741       PTOT = PXZ/DCOS(ALAM)
3742       
3743       A12 = 0.
3744       NTMAX = 0
3745       
3746       ZEA = ZABS
3747       XEA = XV
3748       YEA = YV
3749       PXEA = PX
3750       PYEA = PY
3751       PHIEA = PHI
3752       PXZEA = ASIGN*PXZ
3753       ALAMEA = ALAM
3754
3755       VECT(1) = XV
3756       VECT(2) = YV
3757       VECT(3) = ZABS
3758       VECT(4) = PX/PTOT
3759       VECT(5) = PY/PTOT
3760       VECT(6) = PZ/PTOT
3761       VECT(7) = PTOT
3762       
3763       R = SQRT(VECT(1)*VECT(1)+VECT(2)*VECT(2))
3764
3765       Z = VECT(3)
3766       NSTEP = 0
3767       IX = 0
3768       IY = 0 
3769       IZ = 0
3770 **      PRINT *,' AV GRKUTA ASIGN',ASIGN,' THET',THET
3771       DO ICH = 1,NPLANE 
3772          DO WHILE (Z.GE.ZABS.AND.Z.LT.ZPLANEP(ICH)
3773      &            .AND.NSTEP.LE.NSTEPMAX)
3774 **     &            .AND.(THETA*PITODEG).GT.2.
3775 **     &            .AND. (THETA*PITODEG).LT.9.) 
3776             NSTEP = NSTEP+1 
3777 **          WRITE(6,*) NSTEP,(VECT(I),I=1,7)
3778             CALL RECO_GRKUTA(ASIGN,STEP,VECT,VOUT)
3779             DO I = 1,7
3780                VECT(I) = VOUT(I)
3781             ENDDO   
3782             Z = VECT(3)
3783             R = SQRT(VECT(1)*VECT(1)+VECT(2)*VECT(2))
3784          ENDDO
3785          IF (NSTEP.EQ.NSTEPMAX) RETURN
3786          XP(ICH) = VECT(1)-(Z-ZPLANEP(ICH))*VECT(4)/VECT(6)
3787          YP(ICH) = VECT(2)-(Z-ZPLANEP(ICH))*VECT(5)/VECT(6)
3788          AL   = THICK/ VECT(6)
3789          AP(ICH) = (0.0136D0/PTOT)*DSQRT(AL)*(1+0.038D0*DLOG(AL))
3790       ENDDO
3791 **    PRINT *,' AP GRKUTA ASIGN',ASIGN,' THET',THET
3792
3793
3794 ** Matrice de covariance      
3795       I = 0
3796       DO II = 1,NPLANE
3797         IF (LPLANE(II).EQ.1) THEN
3798            I = I + 1
3799 *     I = II
3800            J = I - 1
3801            DO JJ = II, NPLANE
3802               IF (LPLANE(JJ).EQ.1) THEN
3803                  J = J + 1
3804 *     J = JJ
3805                  COV (I,J) = 0.0D0
3806                  COV (J,I) = A12
3807                  IF (I .EQ. J) THEN
3808                     COV(J,I) =COV(J,I) + XPREC**2
3809                  ENDIF      
3810                  
3811 *     IF (I .EQ. 10 .AND. J .EQ. 10) PRINT *,'10 10   ',COV(J,I)
3812                  DO L = 1,NTMAX
3813                     COV(J,I) = COV(J,I)
3814      &                   +  (ZPLANEP(II) + DISTAZ(L))*(ZPLANEP(JJ) + 
3815      &                   DISTAZ(L))*AMS(L)**2
3816                  ENDDO
3817                  DO K = 1, II-1
3818                     COV(J,I) = COV(J,I)
3819      &                   +  (ZPLANEP(II)-ZPLANEP(K))*
3820      &                   (ZPLANEP(JJ)-ZPLANEP(K))*AP(K)**2
3821 *     IF (I .EQ. 10 .AND. J .EQ. 10) PRINT *,'10 10   ',COV(J,I)
3822                  ENDDO
3823                  COVY(I,J) = 0.0D0
3824                  COVY(J,I) = COV(J,I)
3825                  IF (I .EQ. J) THEN
3826                     COVY(J,I) = COVY(J,I) - XPREC**2 + YPREC**2
3827                  ENDIF
3828               ENDIF   
3829            ENDDO
3830         ENDIF
3831       ENDDO
3832  
3833 *  Inversion des matrices de covariance
3834       NPLU = I
3835  
3836       IFAIL = 0
3837       CALL DSINV(NPLU, COV, NPLANE, IFAIL)
3838 **      IF (JFAIL.NE.0 .AND. IFAIL .NE. 0) STOP 'ERROR'
3839       IF (IFAIL .NE. 0) STOP 'ERROR'
3840       IFAIL = 0
3841       CALL DSINV(NPLU, COVY, NPLANE, IFAIL)
3842 **      IF (JFAIL.NE.0 .AND. IFAIL .NE. 0) STOP 'ERROR'
3843       IF (IFAIL .NE. 0) STOP 'ERROR'
3844 *      PRINT *,' COVARIANCE MATRIX AFTER'
3845 *      DO I = 1, NPLANE
3846 *         PRINT *,(COV(J,I),J=1,NPLANE)
3847 *      ENDDO
3848  
3849 ** Calcul de FVAL ou CHI2
3850       FVAL = 0.0D0
3851       I = 0
3852       DO II = 1,NPLANE
3853         IF (LPLANE(II).EQ.1) THEN
3854         I = I+1
3855 *        I = II
3856         J = 0
3857         DO JJ = 1,NPLANE
3858            IF (LPLANE(JJ).EQ.1) THEN
3859               J = J+1
3860 *             J = JJ
3861               FVAL = FVAL + COV(J,I)*(XMP(II)-XP(II))*(XMP(JJ)-XP(JJ))
3862               FVAL = FVAL + COVY(J,I)*(YMP(II)-YP(II))
3863      &                               *(YMP(JJ)-YP(JJ))
3864 **             IF (JJ.EQ.II) THEN
3865 **                 FVAL = FVAL + (XM(II)-XP(II))*(XM(JJ)-XP(JJ))/XPREC**2
3866 **                 FVAL = FVAL + (YM(II)-YP(II))
3867 **     &                               *(YM(JJ)-YP(JJ))/YPREC**2
3868 **             ENDIF
3869            ENDIF
3870         ENDDO
3871         ENDIF
3872       ENDDO
3873       CHI2 = FVAL
3874
3875 **      IF (CHI2.GT.1.E4) THEN
3876 **         PRINT *,'FCNFIT CHI2= ',CHI2
3877 **         FVAL = 0.
3878 **      ENDIF
3879
3880       
3881  1000 FORMAT(I5,7F12.6)
3882  
3883       RETURN
3884       END
3885
3886 ************************************************************************        
3887       SUBROUTINE NEWFCNFIT(NPAR, GRAD, FVAL, XVAL, IFLAG, FUTIL)
3888 ************************************************************************        
3889 *    With magnetic Field Map GRKUTA
3890 *      trackfinding
3891 *    Calcule FVAL: la fonction minimisee par MINUIT
3892 *    With magnetic field map
3893 *      
3894 ************************************************************************        
3895       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3896
3897       PARAMETER(NPLANE=10)
3898             
3899       COMMON/PARAM/ZPLANEP(NPLANE),THICK,XPREC,YPREC,B0,BL3,ZMAGS,
3900      &             ZMAGE,ZABS,XMAG,ZBP1,ZBP2,CONST
3901      
3902       COMMON/MEAS/LPLANE(NPLANE),XMP(NPLANE),YMP(NPLANE)
3903      
3904       COMMON/FCNOUT/PXZEA,ALAMEA,PHIEA,XEA,YEA,NPLU,CHI2
3905                  
3906       DIMENSION GRAD(*),XVAL(*),AMS(500),DISTAZ(500)
3907   
3908       DIMENSION XP(NPLANE),YP(NPLANE),
3909      &          COV(NPLANE,NPLANE),AP(NPLANE),COVY(NPLANE,NPLANE)
3910       DIMENSION VECT(7),VOUT(7) 
3911       
3912       STEP = 2.   ! 1 cm
3913       NSTEPMAX = 5000
3914       PITODEG = 57.295
3915       XV = XVAL(4)   
3916       YV = XVAL(5)   
3917               
3918       ASIGN = 1.
3919       IF (XVAL(1).LT.0.) ASIGN = -1.
3920       PHI  = XVAL(2)
3921       ALAM = XVAL(3)
3922       PXZ = DABS(1./XVAL(1))
3923            
3924       PX = PXZ*DSIN(PHI)
3925       PY = PXZ*DTAN(ALAM)
3926       PZ = PXZ*DCOS(PHI)
3927       PTOT = PXZ/DCOS(ALAM)
3928       
3929       A12 = 0.
3930       NTMAX = 0
3931       
3932       ZEA = ZABS
3933       XEA = XV
3934       YEA = YV
3935       PXEA = PX
3936       PYEA = PY
3937       PHIEA = PHI
3938       PXZEA = ASIGN*PXZ
3939       ALAMEA = ALAM
3940
3941       VECT(1) = XV
3942       VECT(2) = YV
3943       VECT(3) = ZABS
3944       VECT(4) = PX/PTOT
3945       VECT(5) = PY/PTOT
3946       VECT(6) = PZ/PTOT
3947       VECT(7) = PTOT
3948       
3949       R = SQRT(VECT(1)*VECT(1)+VECT(2)*VECT(2))
3950
3951       Z = VECT(3)
3952       NSTEP = 0
3953       IX = 0
3954       IY = 0 
3955       IZ = 0
3956 **      PRINT *,' AV GRKUTA ASIGN',ASIGN,' THET',THET
3957       DO ICH = 1,NPLANE 
3958          DO WHILE (Z.GE.ZABS.AND.Z.LT.ZPLANEP(ICH)
3959      &            .AND.NSTEP.LE.NSTEPMAX)
3960 **     &            .AND.(THETA*PITODEG).GT.2.
3961 **     &            .AND. (THETA*PITODEG).LT.9.) 
3962             NSTEP = NSTEP+1 
3963 **          WRITE(6,*) NSTEP,(VECT(I),I=1,7)
3964             CALL RECO_GRKUTA (ASIGN,STEP,VECT,VOUT)
3965             DO I = 1,7
3966                VECT(I) = VOUT(I)
3967             ENDDO   
3968             Z = VECT(3)
3969             R = SQRT(VECT(1)*VECT(1)+VECT(2)*VECT(2))
3970          ENDDO
3971          IF (NSTEP.EQ.NSTEPMAX) RETURN
3972          XP(ICH) = VECT(1)-(Z-ZPLANEP(ICH))*VECT(4)/VECT(6)
3973          YP(ICH) = VECT(2)-(Z-ZPLANEP(ICH))*VECT(5)/VECT(6)
3974          AL   = THICK/ VECT(6)
3975          AP(ICH) = (0.0136D0/PTOT)*DSQRT(AL)*(1+0.038D0*DLOG(AL))