2 ** Authors J.P. Cussonneau & P. Lautridou
4 ************** commentaires **********************
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
22 * nombre de candidats a ordonner (a partir de l'ímpulsion)
25 * pour recuperer les candidats
28 * =1 si bon muon, =2 si ghost, =0 autrement
31 * =1 si bonne trace, =0 autrement
34 * permet a partir du numero de la trace de retrouver la numero du hit
37 * compte le nombre de plans touches par trace
40 * permet de retrouver le numero de la trace a partir du numero de hit
43 * pour une trace et une station donnee, dit si la trace est passee dans la chambre
46 * =0 si point (0,0) du vertex impose, sinon coordonnes libres (pour le fit)
49 * numero du hit pour une station et un candidat
52 * nombre de hits par trace au niveau des stations 4 et 5 (3 ou 4)
55 * numero de hit associe a une chambre et a une trace
58 * nombre de bons muons trouves
61 * = nombre de fois ou lón ná pas trouve la bonne trace dans la station
64 * nombre de cas ou pas de hit trouve par station
67 * nombre de fantomes trouves
70 * nombre de resonances dans lácceptance
73 * nombre de resonances trouvees
76 * nombre de traces totales trouvees
79 * nombre total de muons dans lácceptance
82 ****************************************************************
83 SUBROUTINE reconstmuon(IFIT,IDEBUGC,NEV,IDRES,IREADGEANT)
84 ****************************************************************
85 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
92 CALL RECO_READEVT(NEV,IDRES,IREADGEANT)
96 CALL RECO_TRACKF(IDRES,IREADGEANT)
108 ********************************************
110 ********************************************
111 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
112 PARAMETER(NBSTATION=5)
114 COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
116 DATA DZ_PL/8.,8.,24.3,8.,8./
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./
128 DATA ZPLANE/-511.,-686.0,-962.7,-1245.,-1445./ ! dstation=8cm
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
142 ** DATA ZPLANE/-507.,-682.0,-950.7,-1241.,-1441./
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
151 ********************************************
152 SUBROUTINE cutpxz(spxzcut)
153 ********************************************
154 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
156 COMMON/ACUTPXZ/ACUTPXZ
162 ********************************************
163 SUBROUTINE sigmacut(ssigcut)
164 ********************************************
165 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
167 COMMON/TRACKFI/EFF,EFF1,EFF2,XPREC,YPREC,PHIPREC,ALAMPREC,
168 & HCUT,LBKG,SIGCUT,ALPHATOP,HTOP
175 ********************************************
176 SUBROUTINE xpreci(sxprec)
177 ********************************************
178 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
180 COMMON/TRACKFI/EFF,EFF1,EFF2,XPREC,YPREC,PHIPREC,ALAMPREC,
181 & HCUT,LBKG,SIGCUT,ALPHATOP,HTOP
188 ********************************************
189 SUBROUTINE ypreci(syprec)
190 ********************************************
191 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
193 COMMON/TRACKFI/EFF,EFF1,EFF2,XPREC,YPREC,PHIPREC,ALAMPREC,
194 & HCUT,LBKG,SIGCUT,ALPHATOP,HTOP
201 *************************************************************************
202 SUBROUTINE RECO_INIT(seff,sb0,sbl3)
203 *************************************************************************
204 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
206 CALL TRACKF_INIT(seff,sb0,sbl3)
213 *************************************************************************
214 SUBROUTINE TRACKF_INIT(seff,sb0,sbl3)
215 *************************************************************************
216 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
218 PARAMETER(NBSTATION=5,NTRMAX=500)
220 COMMON/REVENT/IEVBKGI,NBKGMAX,MAXUPSEV
224 COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
226 COMMON/FILED/FILERES,FILEBKG,FILEOUT,FILEMIN
228 COMMON/TRACKFI/EFF,EFF1,EFF2,XPREC,YPREC,PHIPREC,ALAMPREC,
229 & HCUT,LBKG,SIGCUT,ALPHATOP,HTOP
231 COMMON/TRACKSUM/NRES(5),NRESF,NTRMUALL,NMUONALL,NGHOSTALL,
232 & NTRACKFALL,NERRALL(NBSTATION),IR
234 COMMON/ACUTPXZ/ACUTPXZ
245 AMAGLEN = ZMAGEND-ZCOIL
246 ZM = AMAGLEN/2.+ZCOIL
247 ALPHATOP = 0.01*0.3*B0*ABS(AMAGLEN)
249 HCUT = ABS(HTOP)/PXZCUT
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
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)
282 *************************************************************************
284 *************************************************************************
287 *************************************************************************
289 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
291 PARAMETER(NPLANE=10,NBSTATION=5)
293 COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
295 COMMON/PARAM/ZPLANEP(NPLANE),THICK,XPREC,YPREC,B0,BL3,ZMAGS,
296 & ZMAGE,ZABS,XMAG,ZBP1,ZBP2,CONST
298 COMMON/PRECCUT/PCUT,PTCUT,CHI2CUT
300 COMMON/PRECSUM/NRESF1,NMUONALL1,NGHOSTALL1,NTRACKFALL1
305 DATA THICK/0.03D0/ ! X/X0=3%
306 ** DATA THICK/0.02D0/ ! X/X0=2% chambre
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
312 ** DATA ZMAGS/825.0D0/,ZMAGE/1125.0D0/,ZABS/503.D0/
315 ** DATA XPREC/0.0100D0/,YPREC/0.144337D0/ ! CCC
316 DATA XPREC/0.0100D0/,YPREC/0.2D0/
319 CONST = 0.299792458D-3*B0*(ZMAGE-ZMAGS)
323 ZPLANEP(J) = -ZPLANE(I)
325 ZPLANEP(J) = -ZPLANE(I)+DZ_PL(I)
328 PCUT = 3. ! Coupure en PXZ muon (GeV/c)
329 PTCUT = 0. ! Coupure en Pt muon (GeV/c)
331 CHI2CUT = 1.E4 ! Coupure sur le CHI2 du fit
334 X02 = 10.397 ! Concrete (cm)
335 X03 = 0.56 ! Plomb (cm)
336 X04 = 47.26 ! Polyethylene (cm)
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
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,
361 PRINT *,' PREC_INIT ZABS (cm)',ZABS,' ZBP1 (cm)',ZBP1,
363 PRINT *,' PREC_INIT Radiation length absorber X01 (cm)',X01,
365 PRINT *,' PREC_INIT X03(cm)',X03,' X04 (cm)',X04,' X05 (cm)',
367 PRINT *,' PREC_INIT Radiation length chamber THICK (%)',
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
382 * Magnetic Field Map GC
383 ** OPEN (UNIT=40,FILE='data/field02.dat',
384 ** & STATUS='UNKNOWN')
393 *************************************************************************
394 SUBROUTINE RECO_READEVT(NEV,IDRES,IREADGEANT)
395 *************************************************************************
398 *************************************************************************
399 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
401 PARAMETER(NTRMAX=500)
403 PARAMETER (NBSTATION=5,MAXIDG=20000,MAXHITTOT=20000,
404 & MAXHITCH=10000,MAXHIT=1000,NBCHAMBER=10)
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)
412 DIMENSION TYPG(MAXIDG),ZCH(MAXIDG)
417 IF (IREADGEANT.eq.1) THEN ! GEANT hits
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)
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)
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)
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)
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)
451 CALL CHFILL (999,SNGL(PTOTG(I)),R1,R2)
454 CALL TRACKF_STAT(IDRES,IREADGEANT)
459 *************************************************************************
460 SUBROUTINE TRACKF_STAT(IDRES,IREADGEANT)
461 *************************************************************************
462 * Associate hits between two chambers inside a station
463 * Simulate spatial resolution and chamber efficiency
465 *************************************************************************
466 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
468 PARAMETER (NBSTATION=5,MAXIDG=20000,MAXHITTOT=20000,
469 & MAXHITCH=10000,MAXHIT=1000,NBCHAMBER=10)
471 COMMON/TRACKFI/EFF,EFF1,EFF2,XPREC,YPREC,PHIPREC,ALAMPREC,
472 & HCUT,LBKG,SIGCUT,ALPHATOP,HTOP
474 COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
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)
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
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)
498 DIMENSION RMIN(NBCHAMBER),RMAX1(NBCHAMBER)
499 DIMENSION XMA(NBCHAMBER,MAXHITCH),YMA(NBCHAMBER,MAXHITCH),
500 & IMARK(NBCHAMBER,MAXHITCH)
502 DIMENSION IEFFI(MAXHITTOT)
503 DIMENSION IH(NBCHAMBER,MAXHIT)
504 DIMENSION NHIT(NBCHAMBER)
505 DIMENSION DXMAX(NBSTATION),DYMAX(NBSTATION),VRES(2,5)
507 REAL*4 RNDM,RN,RN1,RN2,R1,R2
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./
525 RMIN(ICH) = ABS(ZPLANE(IZ)*TAN(2.*ACOS(-1.)/180))
526 IF (IZ.GT.2) RMIN(ICH) = 30.
528 RMIN(ICH) = ABS(ZPLANE(IZ)*TAN(2.*ACOS(-1.)/180))
529 IF (IZ.GT.2) RMIN(ICH) = 30.
537 * 1 ere boucle Lecture des hits initiaux
539 print*,'TRACKF_STAT NHITTOT1=',NHITTOT1
541 IF (IREADGEANT.EQ.1) THEN
552 ISTAK = MOD(ISTAK,30000)
553 ISTAK = MOD(ISTAK,10000)
555 IF ((ITYPG(I).EQ.5.OR.ITYPG(I).EQ.6).AND.
556 & ISTAK.EQ.IDRES) THEN ! upsilon
558 VRES(1,IMU) = XGEANT(I)
559 VRES(2,IMU) = YGEANT(I)
566 DO I = 1,NHITTOT1 ! Boucle sur les hits GEANT de toutes les ch.
568 ** IF (ITYPG(I).NE.5.AND.ITYPG(I).NE.6) GOTO 1 ! CCC
572 IF (IREADGEANT.EQ.1) THEN ! GEANT hits
574 IF (ICH.EQ.9.OR.ICH.EQ.10) THEN
577 DNU = SQRT((XGEANT(I)-VRES(1,IM))**2+
578 & (YGEANT(I)-VRES(2,IM))**2)
579 IF (DNU.LT.DNUM) DNUM = DNU
581 IF (DNUM.GT.20.) GO TO 1
584 CALL RANNOR(RN1,RN2) ! CCCC
585 X = XGEANT(I) + RN1 * XPREC
586 Y = YGEANT(I) + RN2 * YPREC
587 * efficacite des chambres
590 IF (RN.GT.EFF) IEFFI(I) = 0
593 ELSE ! reconstructed hits
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)
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
622 * Association des hits entre chambres d'une station
623 II = 0 ! nombre de hits GEANT par station
625 IZ = INT(FLOAT(ICH1+1)/2.)
640 PVERT1(II) = PVERT1G(I)
641 PVERT2(II) = PVERT2G(I)
642 PVERT3(II) = PVERT3G(I)
643 ZVERT(II) = ZVERTG(I)
645 IF (IEFFI(I).EQ.1) THEN
649 XEXT1 = (ZPLANE(IZ)-DZ_PL(IZ))/ZPLANE(IZ)*X1
650 YEXT1 = (ZPLANE(IZ)-DZ_PL(IZ))/ZPLANE(IZ)*Y1
655 IF (IEFFI(J).EQ.1) THEN
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)
670 IF (DX.LT.DXMAX(IZ).AND.DY.LT.(SIGCUT*DYMAX(IZ)) ! CCC
674 JHIT(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
687 ISTAK = MOD(ISTAK,30000)
688 ISTAK = MOD(ISTAK,10000)
690 IF ((ITYPG(J).EQ.5.OR.ITYPG(J).EQ.6).AND.
691 & ISTAK.EQ.IDRES) THEN ! upsilon
698 PVERT1(II) = PVERT1G(J)
699 PVERT2(II) = PVERT2G(J)
700 PVERT3(II) = PVERT3G(J)
701 ZVERT(II) = ZVERTG(J)
712 JHIT(IZ) = JHIT(IZ)+1
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.
724 CALL CHFILL2(1000+IZ,SNGL(X1),SNGL(Y1),R2)
729 * On conserve les HITS de la 2nde chambre des stations
732 IZ = INT(FLOAT(ICH+1)/2.)
736 IF (IMARK(ICH,I).EQ.0) THEN
746 PVERT1(II) = PVERT1G(J)
747 PVERT2(II) = PVERT2G(J)
748 PVERT3(II) = PVERT3G(J)
749 ZVERT(II) = ZVERTG(I)
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)
756 PHM(IZ,JHIT(IZ)) = 10.
757 ALM(IZ,JHIT(IZ)) = 10.
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)
772 IF (IDEBUG.GE.2) THEN
773 PRINT *,'TRACKF_STAT nb hits:',NHITTOT
776 ** DO IZ = 1,NBSTATION
777 ** PRINT *,' IZ=',IZ,' JHIT(IZ)=',JHIT(IZ)
781 ** PRINT *,' ID(II)=',ID(II)
785 ** DO IZ = 1,NBSTATION
786 ** PRINT *,' IZ=',IZ,' JHIT(IZ)=',JHIT(IZ)
788 ** PRINT *,' PHM(IZ,J)=',PHM(IZ,J),' ALM(IZ,J)=',ALM(IZ,J)
794 *************************************************************************
795 SUBROUTINE TRACKF_STAT_NEW(IDRES)
796 *************************************************************************
797 * Associate hits between two chambers inside a station
798 * Simulate spatial resolution and chamber efficiency
800 *************************************************************************
801 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
803 PARAMETER (NBSTATION=5,MAXIDG=20000,MAXHITTOT=20000,
804 & MAXHITCH=10000,MAXHIT=1000,NBCHAMBER=10)
806 COMMON/TRACKFI/EFF,EFF1,EFF2,XPREC,YPREC,PHIPREC,ALAMPREC,
807 & HCUT,LBKG,SIGCUT,ALPHATOP,HTOP
809 COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
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)
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
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)
833 DIMENSION RMIN(NBCHAMBER),RMAX1(NBCHAMBER)
834 DIMENSION XMA(NBCHAMBER,MAXHITCH),YMA(NBCHAMBER,MAXHITCH),
835 & IMARK(NBCHAMBER,MAXHITCH)
837 DIMENSION IEFFI(MAXHITTOT)
838 DIMENSION IH(NBCHAMBER,MAXHIT)
839 DIMENSION NHIT(NBCHAMBER)
840 DIMENSION DXMAX(NBSTATION),DYMAX(NBSTATION),I2C(1000)
842 DIMENSION DIST(2),NMUON(2),NHITMUON(2,5),NMUONGOOD(2)
844 REAL*4 RNDM,RN,RN1,RN2,R1,R2
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/
857 RMIN(ICH) = ABS(ZPLANE(IZ)*TAN(2.*ACOS(-1.)/180))
858 IF (IZ.GT.2) RMIN(ICH) = 30.
860 RMIN(ICH) = ABS(ZPLANE(IZ)*TAN(2.*ACOS(-1.)/180))
861 IF (IZ.GT.2) RMIN(ICH) = 30.
875 IF (IZCH(I).EQ.NCH) THEN
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
886 NHITMUON(1,NMUON(1))=I
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
895 NHITMUON(2,NMUON(2))=I
901 IF (NMUON(J).GE.2) THEN
902 print*,'j=',j,' nmuon=',nmuon(j)
905 IF (NHITMUON(J,K).NE.NMUONGOOD(J)) IDG(NHITMUON(J,K))=999 ! flag les mauvais hits MUONS
912 * 1 ere boucle Lecture des hits initiaux
914 DO I = 1,NHITTOT1 ! Boucle sur les hits GEANT de toutes les ch.
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)
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
940 print*,' XTRG(I)=', XTRG(I),' YTRG(I)=', YTRG(I),' IDG(I)=',
947 * Association des hits entre chambres d'une station
948 II = 0 ! nombre de hits GEANT par station
950 IZ = INT(FLOAT(ICH1+1)/2.)
965 PVERT1(II) = PVERT1G(I)
966 PVERT2(II) = PVERT2G(I)
967 PVERT3(II) = PVERT3G(I)
968 ZVERT(II) = ZVERTG(I)
970 IF (IEFFI(I).EQ.1) THEN
974 XEXT1 = (ZPLANE(IZ)-DZ_PL(IZ))/ZPLANE(IZ)*X1
975 YEXT1 = (ZPLANE(IZ)-DZ_PL(IZ))/ZPLANE(IZ)*Y1
977 PRINT *,'***** DEBUT RECHERCHE',' ID1=',ID1,' ich1=',ICH1
978 PRINT *,' XTR(II)=', XTR(II),' YTR(II)=', YTR(II)
979 PRINT *,' ITYP(II)=', ITYP(II)
982 IF (IEFFI(J).EQ.1) THEN
986 IF (DX.LT.DXMAX(IZ)) THEN
990 print *,' DX=',DX,' KC=',KC,' ID2=',ID2
1003 IF (DY.LT.DYOLD.AND.DY.LT.(SIGCUT*DYMAX(IZ))) THEN
1006 PRINT *,' ID2=',ID2,' DY=',DY
1009 IF (I2FIND.GT.0) THEN
1015 JHIT(IZ) = JHIT(IZ)+1
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
1031 ISTAK = MOD(ISTAK,30000)
1032 ISTAK = MOD(ISTAK,10000)
1034 IF (ISTAK.EQ.IDRES.AND.ID1.NE.999) THEN
1041 PVERT1(II) = PVERT1G(J)
1042 PVERT2(II) = PVERT2G(J)
1043 PVERT3(II) = PVERT3G(J)
1044 ZVERT(II) = ZVERTG(J)
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.
1064 CALL CHFILL2(1000+IZ,SNGL(X1),SNGL(Y1),R2)
1069 * On conserve les HITS de la 2nde chambre des stations
1072 IZ = INT(FLOAT(ICH+1)/2.)
1076 IF (IMARK(ICH,I).EQ.0) THEN
1077 ** print *,' ich=',ich,' i=',i,' j=',j
1087 PVERT1(II) = PVERT1G(J)
1088 PVERT2(II) = PVERT2G(J)
1089 PVERT3(II) = PVERT3G(J)
1090 ZVERT(II) = ZVERTG(I)
1092 IF (IEFFI(J).EQ.1) THEN
1093 JHIT(IZ) = JHIT(IZ)+1
1094 ** XM(IZ,JHIT(IZ)) = ZPLANE(IZ)/(ZPLANE(IZ)-DZ_PL)
1096 ** YM(IZ,JHIT(IZ)) = ZPLANE(IZ)/(ZPLANE(IZ)-DZ_PL)
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)
1116 IF (IDEBUG.GE.2) THEN
1117 PRINT *,'TRACKF_MICRO nb hits:',NHITTOT
1120 PRINT *,' IZ=',IZ,' JHIT(IZ)=',JHIT(IZ)
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)
1134 *************************************************************************
1135 SUBROUTINE RECO_TRACKF(IDRES,IREADGEANT)
1136 *************************************************************************
1139 *************************************************************************
1140 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
1142 PARAMETER (MAXIDG=20000)
1144 PARAMETER (MAXHITCH=10000,MAXTRK=50000,MAXHITTOT=20000,
1145 & NBSTATION=5,MAXCAN=1000,NBCHAMBER=10,NTRMAX=500)
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)
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
1159 COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
1161 COMMON/VERIFGEANT/ITTROUGH(MAXTRK,NBSTATION),
1162 & IT_LIST(MAXTRK),IT_NP(MAXTRK),ITCHECK(MAXTRK),
1165 COMMON/HCHHIT/HHIT(MAXHITCH),INDEXTAB(MAXHITCH),INDEXMAX
1167 COMMON/PRECIS/EEXM(NBSTATION),EEYM(NBSTATION),EEPH(NBSTATION),
1170 COMMON/MEASUR/XMES(NBSTATION),YMES(NBSTATION),IZMES(NBSTATION),
1171 & PHMES(NBSTATION),ALMES(NBSTATION),MPOS(NBSTATION),
1174 COMMON/PLANE/XPL(NBSTATION,2),YPL(NBSTATION,2),
1175 & PHPL(NBSTATION),ALPL(NBSTATION),CHI2PL
1177 COMMON/CANDIDAT/JCAN(NBSTATION,MAXCAN),JCANTYP(MAXCAN),
1178 & EEX(MAXCAN),EEY(MAXCAN),EEP(MAXCAN),EEA(MAXCAN)
1180 COMMON /VERTEX/ERRV,IVERTEX
1182 COMMON/TRACKSUM/NRES(5),NRESF,NTRMUALL,NMUONALL,NGHOSTALL,
1183 & NTRACKFALL,NERRALL(NBSTATION),IR
1185 COMMON/TRACKFOUT/IEVOUT,NTREVT,JJOUT(NBCHAMBER,NTRMAX),
1186 & ISTAT(NTRMAX),PXZOUT(NTRMAX),TPHIOUT(NTRMAX),
1187 & TALAMOUT(NTRMAX),XVERTOUT(NTRMAX),YVERTOUT(NTRMAX),
1189 & XMESOUT(NBCHAMBER,NTRMAX),YMESOUT(NBCHAMBER,NTRMAX)
1190 & ,PXVOUT(NTRMAX),PYVOUT(NTRMAX),PZVOUT(NTRMAX)
1192 COMMON/TRACKFI/EFF,EFF1,EFF2,XPREC,YPREC,PHIPREC,ALAMPREC,
1193 & HCUT,LBKG,SIGCUT,ALPHATOP,HTOP
1195 COMMON/DEBEVT/IDEBUG
1197 DIMENSION PEST(5),PSTEP(5),NERR(NBSTATION)
1202 ** GEANT informations
1210 * BOUCLE SUR LES HITS
1214 IF (ID(IH).EQ.IT_LIST(IT)) GOTO 4
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
1220 ITTROUGH(NTRTOT,II) = 0
1223 4 IT_NP(IT) = IT_NP(IT)+1
1224 ITTROUGH(IT,IZST(IH))=IH ! =IH si la trace IT touche le plan IZST
1227 IF (IDEBUG.GE.2) THEN
1228 PRINT *,'RECO_TRACKF nb total de trace GEANT:',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)
1240 IF ((ITTROUGH(IT,1)*ITTROUGH(IT,2)*ITTROUGH(IT,3)*
1241 & ITTROUGH(IT,4)*ITTROUGH(IT,5)).NE.0)
1244 IH = ITTROUGH(IT,NBSTATION)
1245 IF (ITYP(IH).EQ.5.OR.ITYP(IH).EQ.6) THEN
1247 ISTAK = MOD(ISTAK,30000)
1248 ISTAK = MOD(ISTAK,10000)
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)
1254 * print*,'istak=',istak
1256 IF (ISTAK.EQ.IDRES) THEN ! psi or upsilon
1259 NTRMUALL = NTRMUALL+1
1266 IF (IDEBUG.GE.2) THEN
1267 PRINT *,'RECO_TRACKF nb of part. GEANT crossing 5 st.:',
1269 PRINT *,'RECO_TRACKF nb of muons/res. GEANT crossing 5 st.:',
1273 ** CALL H_ACCEPTANCE(5)
1274 ** CALL H_ACCEPTANCE(4)
1281 CALL ORDONNE_HIT(5,HCUT)
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))
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))
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
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
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.)
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.)
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
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)
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.)
1335 ** P2 = (HTOP/HHIT(JJ))**2
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)
1344 * renvoie le num de hit de 4 le plus pres dans le domaine de recherche
1346 CALL DISTMIN4(X1,Y1,PHM(5,JJ),ALM(5,JJ),4,EX1,EY1,EPH,EAL,
1349 CALL CHECK_HISTO4(11,4,IFIND,5,JJ,X1,Y1,PHM(5,JJ),ALM(5,JJ),P1,
1351 IF (IFIND.GT.0) THEN
1352 CALL STOCK_CANDIDAT(5,JJ,4,IFIND,EX1,EY1,EPH,EAL,NCAN,1)
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,
1357 IF (IFIND.GT.0) THEN
1358 CALL STOCK_CANDIDAT(5,JJ,4,IFIND,EX1,EY1,EPH,EAL,NCAN,2)
1365 CALL ORDONNE_HIT(4,HCUT)
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))
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))
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.
1386 EPH2 = 2.0*PHIPREC**2 + (HCONST*HHIT(JJ))**2
1387 EAL2 = 2.0*ALAMPREC**2 + (HCONST*HHIT(JJ))**2
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
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.)
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.)
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
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)
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.)
1430 ** P2 = (HTOP/HHIT(JJ))**2
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)
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)
1441 CALL CHECK_HISTO2(0,5,IFIND,4,JJ,X1,Y1,X2,Y2,P1,EX1,EY1,
1443 IF (IFIND.GT.0) THEN
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
1449 CALL STOCK_CANDIDAT(4,JJ,5,IFIND,EX1,EY1,EPH,EAL,NCAN,3)
1459 IF (JJ4.GT.0.AND.JJ5.GT.0) THEN
1462 IT = ITRACK(IP(5,JJ5))
1463 IF (ITCHECK(IT).EQ.1) THEN
1464 IF (ID4.EQ.ID5) THEN
1470 IF (NMU45.GE.2) NRES(1) = NRES(1)+1
1472 IF (IDEBUG.GE.2) THEN
1473 PRINT *,'RECO_TRACKF nb candidat recherche 4->5 et 5->4 :'
1475 PRINT *,'RECO_TRACKF nb of good muons 4->5 et 5->4 :'
1482 * -- Boucle sur les candidats (4,5) NCAN
1495 DO ICH = 1,NBSTATION
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
1516 IF (IZM(4,JJ4).EQ.1) THEN
1519 ZPL4 = ZPLANE(4)-DZ_PL(4)
1521 IF (IZM(5,JJ5).EQ.1) THEN
1524 ZPL5 = ZPLANE(5)-DZ_PL(5)
1526 TPHEST = (XM(5,JJ5) - XM(4,JJ4))/(ZPL5-ZPL4)
1527 PHEST = ATAN(TPHEST)
1528 TALEST = -(YM(5,JJ5) - YM(4,JJ4))*COS(PHEST)
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
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
1543 IZMES(4) = IZM(4,JJ4)
1544 PHMES(4) = PHM(4,JJ4)
1545 ALMES(4) = ALM(4,JJ4)
1548 IZMES(5) = IZM(5,JJ5)
1549 PHMES(5) = PHM(5,JJ5)
1550 ALMES(5) = ALM(5,JJ5)
1552 IVERTEX = 0 ! Vertex = (0,0,0)
1554 * -- Fit 4,5,V pour la recherche dans 3
1556 CALL TRACKF_FIT(IVERTEX,PEST,PSTEP,PXZINV45,TPHI45,TALAM45,
1559 * -- Recherche dans la station 3
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
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.)
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.)
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.)
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.)
1599 ** EEXM(IEX) = EEXM(3)
1600 ** EEYM(IEX) = EEYM(3)
1601 ** EEPH(IEX) = EEPH(3)
1602 ** EEAL(IEX) = EEAL(3)
1611 CALL DISTMIN4(X1,Y1,PHI1,ALAM1,3,EXM,EYM,EPH,EAL,IPA3)
1613 P1 = PTOT(IP(5,JJ5))
1615 CALL CHECK_HISTO4(61,3,IPA3,5,JJ5,X1,Y1,PHI1,ALAM1,P1,
1622 IZMES(3) = IZM(3,JJ3)
1623 PHMES(3) = PHM(3,JJ3)
1624 ALMES(3) = ALM(3,JJ3)
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,
1637 IZMES(3) = IZM(3,JJ3)
1645 * -- Fit 3,4,5 pour la recherche dans 2
1649 IF (JJ5.GT.0.AND.JJ4.GT.0.AND.JJ3.GT.0) THEN
1653 IT = ITRACK(IP(5,JJ5))
1654 IF (ITCHECK(IT).EQ.1) THEN
1655 IF ((ID5.EQ.ID3).AND.(ID5.EQ.ID4)) THEN
1674 CALL TRACKF_FIT(IVERTEX,PEST,PSTEP,PXZINV345,TPHI345,
1675 & TALAM345,XVERT,YVERT)
1677 * -- Recherche dans la st. 2
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
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.)
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.)
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.)
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.)
1717 ** EEXM(IEX) = EEXM(2)
1718 ** EEYM(IEX) = EEYM(2)
1719 ** EEPH(IEX) = EEPH(2)
1720 ** EEAL(IEX) = EEAL(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,
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
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.)
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.)
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.)
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.)
1768 ** EEXM(IEX) = EEXM(1)
1769 ** EEYM(IEX) = EEYM(1)
1770 ** EEPH(IEX) = EEPH(1)
1771 ** EEAL(IEX) = EEAL(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,
1781 * -- A partir de P+A dans la st. 2
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)
1802 * -- Fit V,2,3,4,5 pour la recherche dans 1
1804 CALL TRACKF_FIT(IVERTEX,PEST,PSTEP,PXZINV,TPHI,TALAM,
1807 * !!! ATTENTION aux erreurs
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
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.)
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.)
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.)
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.)
1848 ** EEXM(IEX) = EEXM(1)
1849 ** EEYM(IEX) = EEYM(1)
1850 ** EEPH(IEX) = EEPH(1)
1851 ** EEAL(IEX) = EEAL(1)
1860 CALL DISTMIN4(X1,Y1,PHI1,ALAM1,1,EXM,EYM,EPH,EAL,
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
1869 IZMES(1) = IZM(1,JJ1)
1870 PHMES(1) = PHM(1,JJ1)
1871 ALMES(1) = ALM(1,JJ1)
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)
1881 * -- A partir de P+A dans la st. 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)
1904 * -- Fit V,1,3,4,5 pour la recherche dans 2
1906 CALL TRACKF_FIT(IVERTEX,PEST,PSTEP,PXZINV,TPHI,TALAM,
1909 * !!! ATTENTION aux erreurs
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)
1921 EEXM(2) = EXM/SIGCUT
1922 EEYM(2) = EYM/SIGCUT
1923 EEPH(2) = EPH/SIGCUT
1924 EEAL(2) = EAL/SIGCUT
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.)
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.)
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.)
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.)
1954 ** EEXM(IEX) = EEXM(2)
1955 ** EEYM(IEX) = EEYM(2)
1956 ** EEPH(IEX) = EEPH(2)
1957 ** EEAL(IEX) = EEAL(2)
1967 CALL DISTMIN4(X1,Y1,PHI1,ALAM1,2,EXM,EYM,EPH,EAL,
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
1977 IZMES(2) = IZM(2,JJ2)
1978 PHMES(2) = PHM(2,JJ2)
1979 ALMES(2) = ALM(2,JJ2)
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)
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
1996 IZMES(1) = IZM(1,JJ1)
2001 ELSEIF(IPA1.GT.0.AND.IPA1P2.GT.0) THEN
2006 IZMES(2) = IZM(2,JJ2)
2016 NTRACKFOLD = NTRACKF
2019 CALL CHECK_MUON(JJ1,JJ2,JJ3,JJ4,JJ5,NTRACKF,NMUONF,
2020 & NGHOST,NERR,PXZINV,TPHI,TALAM,XVERT,YVERT)
2022 IF (NTRACKF.GT.NTRACKFOLD) THEN
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
2037 ENDDO ! end loop on candidats NCAN
2039 IF (NMU345.GE.2) NRES(2) = NRES(2) + 1
2041 IF (IDEBUG.GE.2) THEN
2042 PRINT *,'RECO_TRACKF nb of good muons 3 4 5 :'
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
2052 PRINT *,'RECO_TRACKF nb of error in st',I,'/evt:',NERR(I)
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
2061 NERRALL(I) = NERRALL(I)+NERR(I)
2064 CALL TRACKF_OUT(IR,NTRACKF)
2069 **********************************************************************
2070 SUBROUTINE TRACKF_OUT(IR,NTRACKF)
2071 **********************************************************************
2074 **********************************************************************
2075 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2077 PARAMETER (MAXHITCH=10000,MAXTRK=50000,MAXHITTOT=20000,
2079 PARAMETER (NBCHAMBER=10,NTRMAX=500)
2081 COMMON/RECOUT/JJO(NBCHAMBER,NTRMAX),XMESO(NBCHAMBER,NTRMAX),
2082 & YMESO(NBCHAMBER,NTRMAX)
2084 COMMON/TRACKFOUT/IEVOUT,NTREVT,JJOUT(NBCHAMBER,NTRMAX),
2085 & ISTAT(NTRMAX),PXZOUT(NTRMAX),TPHIOUT(NTRMAX),
2086 & TALAMOUT(NTRMAX),XVERTOUT(NTRMAX),YVERTOUT(NTRMAX),
2088 & XMESOUT(NBCHAMBER,NTRMAX),YMESOUT(NBCHAMBER,NTRMAX)
2089 & ,PXVOUT(NTRMAX),PYVOUT(NTRMAX),PZVOUT(NTRMAX)
2093 IF (NTREVT.GT.0) THEN
2096 DO IST = 1,NBSTATION
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)
2112 **********************************************************************
2113 SUBROUTINE CHECK_MUON(JJ1,JJ2,JJ3,JJ4,JJ5,
2114 & NTRACKF,NMUONF,NGHOST,NERR,PXZINV,TPHI,TALAM,
2116 **********************************************************************
2117 * Check muon track candidate using GEANT informations
2119 **********************************************************************
2120 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2122 PARAMETER (MAXHITCH=10000,MAXTRK=50000,MAXHITTOT=20000,
2124 PARAMETER (NBCHAMBER=10,NTRMAX=500)
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)
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
2137 COMMON/VERIFGEANT/ITTROUGH(MAXTRK,NBSTATION),
2138 & IT_LIST(MAXTRK),IT_NP(MAXTRK),ITCHECK(MAXTRK),
2141 COMMON/MEASUR/XMES(NBSTATION),YMES(NBSTATION),IZMES(NBSTATION),
2142 & PHMES(NBSTATION),ALMES(NBSTATION),MPOS(NBSTATION),
2145 COMMON/RECOUT/JJO(NBCHAMBER,NTRMAX),XMESO(NBCHAMBER,NTRMAX),
2146 & YMESO(NBCHAMBER,NTRMAX)
2148 COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
2150 COMMON/MAGNET/BL3,B0
2155 DIMENSION NERR(NBSTATION),JJK(NBSTATION)
2157 * print*,' *** appel de la subroutine check_muon ***'
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))
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)
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)
2187 JJO(ICH2,NTRACKF) = JJK(IST)
2188 XMESO(ICH2,NTRACKF) = XMR(IST,JJK(IST),2)
2189 YMESO(ICH2,NTRACKF) = YMR(IST,JJK(IST),2)
2193 CALL CHFILL(700,SNGL(XVERT),R1,R2)
2194 CALL CHFILL(701,SNGL(YVERT),R1,R2)
2196 IF (ITCHECK(IT).EQ.1) THEN
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
2204 NMUONF = NMUONF+1 ! Bon muon
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)
2216 NGHOST = NGHOST+1 ! ghost
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)
2229 IF (ITCHECK(ITRACK(IP(5,JJ5))).EQ.1) THEN
2230 IF (JJ3.EQ.0) NERR(3) = NERR(3)+1
2232 IF (JJ1.EQ.0) NERR(1) = NERR(1)+1
2233 IF (JJ2.EQ.0) NERR(2) = NERR(2)+1
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)
2249 *************************************************************************
2251 *************************************************************************
2254 *************************************************************************
2255 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2257 PARAMETER (MAXHITCH=10000,MAXTRK=50000,MAXHITTOT=20000,
2259 PARAMETER (NBCHAMBER=10,NTRMAX=500)
2261 COMMON/TRACKFOUT/IEVOUT,NTREVT,JJOUT(NBCHAMBER,NTRMAX),
2262 & ISTAT(NTRMAX),PXZOUT(NTRMAX),TPHIOUT(NTRMAX),
2263 & TALAMOUT(NTRMAX),XVERTOUT(NTRMAX),YVERTOUT(NTRMAX),
2265 & XMESOUT(NBCHAMBER,NTRMAX),YMESOUT(NBCHAMBER,NTRMAX)
2266 & ,PXVOUT(NTRMAX),PYVOUT(NTRMAX),PZVOUT(NTRMAX)
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)
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
2279 COMMON/RECOUT/JJO(NBCHAMBER,NTRMAX),XMESO(NBCHAMBER,NTRMAX),
2280 & YMESO(NBCHAMBER,NTRMAX)
2282 COMMON/TRACKSUM/NRES(5),NRESF,NTRMUALL,NMUONALL,NGHOSTALL,
2283 & NTRACKFALL,NERRALL(NBSTATION),IR
2285 COMMON/PRECSUM/NRESF1,NMUONALL1,NGHOSTALL1,NTRACKFALL1
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)
2292 COMMON/DEBEVT/IDEBUG
2297 NTRACKFALL1 = NTRACKFALL1 + 1
2298 IF (ISTAT(ITR).EQ.1) THEN
2300 NMUONALL1 = NMUONALL1 + 1
2301 ELSEIF (ISTAT(ITR).EQ.2) THEN
2303 NGHOSTALL1 = NGHOSTALL1 + 1
2306 IF (NMUF.EQ.2) NRESF1 = NRESF1 + 1
2310 ISTATR(ITR) = ISTAT(ITR)
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)
2320 CHI2R(ITR) = CHI2OUT(ITR)
2321 PXV(ITR) = PXVOUT(ITR)
2322 PYV(ITR) = PYVOUT(ITR)
2323 PZV(ITR) = PZVOUT(ITR)
2326 CALL CHFNT(IEVR,NTRACKR,ISTATR,ISIGNR,
2327 & PXR,PYR,PZR,ZVR,CHI2R,PXV,PYV,PZV)
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
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
2358 *************************************************************************
2359 SUBROUTINE RECO_TERM
2360 *************************************************************************
2363 *************************************************************************
2364 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2366 PARAMETER(NBSTATION=5)
2368 COMMON/TRACKSUM/NRES(5),NRESF,NTRMUALL,NMUONALL,NGHOSTALL,
2369 & NTRACKFALL,NERRALL(NBSTATION),IR
2371 COMMON/PRECSUM/NRESF1,NMUONALL1,NGHOSTALL1,NTRACKFALL1
2373 COMMON/DEBEVT/IDEBUG
2375 CHARACTER*50 FILEBKG,FILERES,FILEOUT,FILEMIN
2377 IF (IDEBUG.GE.1) THEN
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
2389 PRINT *,'RECO_TERM nb of error in st',I,':',NERRALL(I)
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
2405 *************************************************************************
2406 SUBROUTINE RECO_PRECISION
2407 *************************************************************************
2410 *************************************************************************
2411 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2413 PARAMETER (MAXHITCH=10000,MAXTRK=50000,MAXHITTOT=20000,
2414 & NBSTATION=5,MAXCAN=1000,NTRMAX=500)
2415 PARAMETER (NPLANE=10)
2417 COMMON/PARAM/ZPLANEP(NPLANE),THICK,XPREC,YPREC,B0,BL3,ZMAGS,
2418 & ZMAGE,ZABS,XMAG,ZBP1,ZBP2,CONST
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)
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
2431 COMMON/TRACKFOUT/IEVOUT,NTREVT,JJOUT(NPLANE,NTRMAX),
2432 & ISTAT(NTRMAX),PXZOUT(NTRMAX),TPHIOUT(NTRMAX),
2433 & TALAMOUT(NTRMAX),XVERTOUT(NTRMAX),YVERTOUT(NTRMAX),
2435 & XMESOUT(NPLANE,NTRMAX),YMESOUT(NPLANE,NTRMAX)
2436 & ,PXVOUT(NTRMAX),PYVOUT(NTRMAX),PZVOUT(NTRMAX)
2439 COMMON/MEAS/LPLANE(NPLANE),XMP(NPLANE),YMP(NPLANE)
2441 COMMON/FCNOUT/PXZEA,ALAMEA,PHIEA,XEA,YEA,NPLU,CHI2
2443 COMMON/PRECCUT/PCUT,PTCUT,CHI2CUT
2445 DIMENSION PARMU(MAXCAN,NPLANE,2),LPLANEMU(MAXCAN,NPLANE)
2447 IF (NTREVT.EQ.0) RETURN
2451 DO IST = 1,NBSTATION
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)
2461 LPLANEMU(ITR,ICH) = 0
2470 LPLANE(ICH) = LPLANEMU(ICAN,ICH)
2471 XMP(ICH) = PARMU(ICAN,ICH,1)
2472 YMP(ICH) = PARMU(ICAN,ICH,2)
2475 IF (LPLANE(1).GT.0) THEN
2484 IF (LPLANE(3).GT.0) THEN
2493 IF (LPLANE(7).GT.0) THEN
2500 IF (LPLANE(9).GT.0) THEN
2508 PHIAV = DATAN2((X2-X1),(ZPLANEP(IPL2)-ZPLANEP(IPL1)))
2509 PHIAP = DATAN2((X4-X3),(ZPLANEP(IPL4)-ZPLANEP(IPL3)))
2510 DPHI = (PHIAP-PHIAV)
2512 IF (DPHI.LT.0.) ASIGN = -1. ! CCC
2513 PXZ = CONST/DABS(DPHI)
2515 IF (PXZ.LT.PCUT) GO TO 2
2517 PXZINVI = ASIGN/PXZ ! CCC
2518 ** PXZINVI = 1./PXZOUT(ICAN) ! CCC
2519 ** PXZINVI = -1./49. ! CCC
2521 ALAMI = DATAN2((Y2-Y1),DSQRT((X2-X1)**2
2522 & +(ZPLANEP(IPL2)-ZPLANEP(IPL1))**2))
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
2532 PRINT *,' RECO_PRECISION CHI2OUT avant fit ',CHI2OUT(ICAN)
2534 IF (CHI2OUT(ICAN).GT.CHI2CUT) GO TO 2
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,
2541 ** Correction de Branson
2542 CALL BRANSON(PXZEA,PHIEA,ALAMEA,XEA,YEA)
2545 PX1 = PXZ1*DSIN(PHIEA)
2546 PY1 = PXZ1*DTAN(ALAMEA)
2547 PT1 = DSQRT(PX1**2+PY1**2)
2549 IF (PT1.LT.PTCUT) GO TO 2
2551 IF ((CHI2/FLOAT(2*NPLU-5)).GT.CHI2CUT) GO TO 2
2555 JJOUT(ICH,NTRACK) = JJOUT(ICH,ICAN)
2556 XMESOUT(ICH,NTRACK) = XMESOUT(ICH,ICAN)
2557 YMESOUT(ICH,NTRACK) = YMESOUT(ICH,ICAN)
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)
2567 PRINT *,' RECO_PRECISION CHI2OUT apres fit',CHI2OUT(NTRACK)
2569 ** print *,' reco_precision pxz tphi talam xvert yvert chi2',
2570 ** & PXZEA,PHIEA,ALAMEA,
2571 ** & XEA,YEA,CHI2/FLOAT(2*NPLU-5)
2580 ************************************************************************
2581 SUBROUTINE ORDONNE_HIT(ICH,HCUT)
2582 **************************************************************************
2584 * Sort hits in station ICH according to the "impact parameter" HHIT
2586 **************************************************************************
2588 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2590 PARAMETER (MAXHITCH=10000,MAXTRK=50000,MAXHITTOT=20000,
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)
2599 COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
2601 COMMON/VERIFGEANT/ITTROUGH(MAXTRK,NBSTATION),
2602 & IT_LIST(MAXTRK),IT_NP(MAXTRK),ITCHECK(MAXTRK),
2605 COMMON/HCHHIT/HHIT(MAXHITCH),INDEXTAB(MAXHITCH),INDEXMAX
2607 COMMON/DEBEVT/IDEBUG
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
2616 * boucle sur le nombre de hits candidats de la station
2620 IF (PHM(ICH,J).LT.6.3) THEN !2pi=6.3 radian
2622 * calcul du h dans XOY a z=0
2623 HHIT(J)=ABS(XM(ICH,J)-ZPLANE(ICH)*PHM(ICH,J))
2625 IF (HHIT(J).LT.HCUT) THEN
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,
2640 H4(I) = SNGL(HHIT(I))
2643 CALL SORTZV(H4,INDEXTAB,INDEXMAX,-1,0,1)
2646 HHIT(I) = DBLE(H4(I))
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,':',
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
2664 ***********************************************************************************
2666 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2668 PARAMETER (MAXHITCH=10000,NBSTATION=5)
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)
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
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
2703 ***********************************************************************************
2704 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2706 PARAMETER (MAXHITCH=10000,NBSTATION=5)
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)
2717 IF (IZM(ICH,I).EQ.1) THEN
2726 IF (ICH.EQ.4.OR.ICH.EQ.5) THEN
2727 IF (IZM(ICH,I).EQ.1) THEN
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
2748 ********************************************************************************
2749 SUBROUTINE H_ACCEPTANCE(ICH)
2750 ********************************************************************************
2751 * Etude de l'acceptance des resonnances en fonction du H
2752 * dans la station ICH
2756 ********************************************************************************
2758 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2760 PARAMETER (MAXHITCH=10000,MAXTRK=50000,MAXHITTOT=20000,
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)
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
2775 COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
2777 COMMON/VERIFGEANT/ITTROUGH(MAXTRK,NBSTATION),
2778 & IT_LIST(MAXTRK),IT_NP(MAXTRK),ITCHECK(MAXTRK),
2781 COMMON/HCHHIT/HHIT(MAXHITCH),INDEXTAB(MAXHITCH),INDEXMAX
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))
2798 * PRINT *,'ACCEPT. nb de muons/res total dans st.',ICH,':',NMUONI
2803 * Sort hits in st. z
2804 CALL ORDONNE_HIT(ICH,HCUT)
2807 IIND = IP(ICH,INDEXTAB(IND))
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
2817 IF (NMUON.EQ.2.AND.NMUONI.EQ.2) THEN
2818 CALL CHFILL(ICH*100,SNGL(HCUT),R1,R2)
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.
2832 ********************************************************************************
2833 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2835 PARAMETER(NBSTATION=5)
2837 DIMENSION XPL(NBSTATION,2),YPL(NBSTATION,2),PHPL(NBSTATION),
2838 & ALPL(NBSTATION),PEST(NBSTATION)
2840 COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
2842 COMMON /MEASUR/XMES(NBSTATION),YMES(NBSTATION),IZMES(NBSTATION),
2843 & PHMES(NBSTATION),ALMES(NBSTATION),MPOS(NBSTATION),
2845 COMMON /MAGNET/BL3,B0
2847 LOGICAL LFLAG, LFLAG1
2857 PX = -ABS(PXZ)*SIN(PHI)
2858 PZ = -ABS(PXZ)*COS(PHI)
2859 PXY = SQRT(PX**2 + PY**2)
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)
2872 * PRINT *, XC,YC,RS,FI,TTHET,PXY,PZ
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)
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
2891 PX = PX0 + DFR * (PY0 - 0.5D0*PX0*DFR/RS) / RS
2892 PY = PY0 - DFR * (PX0 + 0.5D0*PY0*DFR/RS) / RS
2894 ALPL(J)=ATAN(PY/SQRT(PX**2+PZ**2))
2898 * DFI = (ZCOIL-ZSTR)/(TTHET*RS)
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
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)
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'
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
2931 FIN = ASIN((ZMAGEND - ZC)/RD)
2932 XMAGEND = XC - RD*COS(FIN)
2933 YMAGEND = YCOIL - RD*(FIN - PHI)*TALAM
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)*
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.
2958 ********************************************************************************
2959 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2961 PARAMETER(NBSTATION=5,NPLANE=10)
2963 DIMENSION XPL(NBSTATION,2),YPL(NBSTATION,2),PHPL(NBSTATION),
2964 & ALPL(NBSTATION),PEST(NBSTATION)
2966 COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
2968 COMMON/PARAM/ZPLANEP(NPLANE),THICK,XPREC,YPREC,B0,BL3,ZMAGS,
2969 & ZMAGE,ZABS,XMAG,ZBP1,ZBP2,CONST
2971 COMMON /MEASUR/XMES(NBSTATION),YMES(NBSTATION),IZMES(NBSTATION),
2972 & PHMES(NBSTATION),ALMES(NBSTATION),MPOS(NBSTATION),
2976 DIMENSION VECT(7),VOUT(7)
2982 IF (PEST(1).LT.0.) ASIGN = -1.
2987 PXZ = DABS(1./PEST(1))
2992 PTOT = PXZ/DCOS(ALAM)
3006 ** PRINT *,' AV GRKUTA ASIGN',ASIGN,' THET',THET
3010 IST = INT(FLOAT(ICH+1)/2.)
3013 DO WHILE (Z.GE.0..AND.Z.LT.ABS(ZPLANEP(ICH))
3014 & .AND.NSTEP.LE.NSTEPMAX)
3016 ** WRITE(6,*) NSTEP,(VECT(I),I=1,7)
3017 CALL RECO_GRKUTA (ASIGN,STEP,VECT,VOUT)
3023 IF (IST.NE.ISTOLD) THEN
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)
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))
3038 ** print *,' vect= ',vect(1),vect(2),vect(3),vect(4),vect(5),
3039 ** & vect(6),vect(7)
3044 *******************************************************************************
3045 SUBROUTINE FCN(NPAR,GRAD,FVAL,PEST,IFLAG,FUTIL)
3046 *******************************************************************************
3047 * Calculate FVAL=CHI2 the function minimized by minuit for a given track
3049 *******************************************************************************
3050 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
3052 PARAMETER(NBSTATION=5)
3054 * DIMENSION PEST(*),GRAD(*)
3055 DIMENSION PEST(5),GRAD(5)
3056 DIMENSION PEEST(NBSTATION)
3058 COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
3061 COMMON/PRECIS/EEXM(NBSTATION),EEYM(NBSTATION),EEPH(NBSTATION),
3064 COMMON /MEASUR/XMES(NBSTATION),YMES(NBSTATION),IZMES(NBSTATION),
3065 & PHMES(NBSTATION),ALMES(NBSTATION), MPOS(NBSTATION),
3068 COMMON /PLANE/XPL(NBSTATION,2),YPL(NBSTATION,2),PHPL(NBSTATION),
3069 & ALPL(NBSTATION),CHI2PL
3071 COMMON/VERTEX/ERRV,IVERTEX
3076 IF(IVERTEX.EQ.1) THEN
3077 PEEST(4)=PEST(4) ! position du vertex
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)
3092 * IF (IVERTEX.EQ.1) THEN
3093 * FVAL = (PEST(4)/ERRV)**2 + (PEST(5)/ERRV)**2
3099 IF (MPOS(J).EQ.1) THEN
3101 XPLC = XPL(J,IZMES(J))
3102 YPLC = YPL(J,IZMES(J))
3103 FF = (XMES(J) - XPLC)/EEXM(J)
3105 FF = (YMES(J) - YPLC)/EEYM(J)
3108 IF (MANG(J).EQ.1) THEN
3110 FF = (PHMES(J) - PHPL(J))/EEPH(J)
3112 FF = (ALMES(J) - ALPL(J))/EEAL(J)
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)
3122 IF (IVERTEX.EQ.1) NPARAM = 5
3123 CHI2PL = FVAL/FLOAT(2*NPLPL-NPARAM)
3127 ********************************************************************************
3128 SUBROUTINE STOCK_CANDIDAT(ICH1,JHITCH1,ICH2,IFIND,EXM,EYM,EPH,EAL,
3130 ********************************************************************************
3131 * Fill common CANDIDAT with track candidates from the search in stations 4&5
3133 ********************************************************************************
3135 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
3137 PARAMETER (MAXHITTOT=20000,MAXHITCH=10000,NBSTATION=5,MAXCAN=1000)
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)
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
3151 COMMON/CANDIDAT/JCAN(NBSTATION,MAXCAN),JCANTYP(MAXCAN),
3152 & EEX(MAXCAN),EEY(MAXCAN),EEP(MAXCAN),EEA(MAXCAN)
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
3159 JCAN(ICH1,NCAN) = JHITCH1
3160 JCAN(ICH2,NCAN) = IFIND
3161 JCANTYP(NCAN) = ICODE
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
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.
3179 *****************************************************************************
3180 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
3182 PARAMETER (MAXHITCH=10000,MAXHITTOT=20000,NBSTATION=5,
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)
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
3196 COMMON/VERIFGEANT/ITTROUGH(MAXTRK,NBSTATION),
3197 & IT_LIST(MAXTRK),IT_NP(MAXTRK),ITCHECK(MAXTRK),
3200 COMMON/DEBEVT/IDEBUG
3208 IF (PHM(ICH2,I).LE.6.3) THEN
3209 IF (ID(IP(ICH1,IHIT1)).EQ.ID(IP(ICH2,I))) THEN
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)
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),
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))
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
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.
3272 *****************************************************************************
3273 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
3275 PARAMETER (MAXHITCH=10000,MAXHITTOT=20000,NBSTATION=5,
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)
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
3289 COMMON/VERIFGEANT/ITTROUGH(MAXTRK,NBSTATION),
3290 & IT_LIST(MAXTRK),IT_NP(MAXTRK),ITCHECK(MAXTRK),
3293 COMMON/DEBEVT/IDEBUG
3301 IF (ID(IP(ICH1,IHIT1)).EQ.ID(IP(ICH2,I))) THEN
3303 IF (IDHIST.GT.0) THEN
3304 IF (IZM(ICH2,I).EQ.1) THEN
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)
3325 IF (ITCHECK(ITRACK(IP(ICH1,IHIT1))).EQ.1) THEN
3328 IF (IZM(ICH2,JOK).EQ.1) THEN
3331 IF (ICH2.EQ.4.OR.ICH2.EQ.5) THEN
3338 IF (ICH2.EQ.4.OR.ICH2.EQ.5) THEN
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))
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
3374 ************************************************************************
3376 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3378 REA = DSQRT(XEA**2+YEA**2)
3379 IF (REA.GT.26.3611) THEN
3382 P = P+SPB/1000.*11.35*(16.66D-3 * P+1.33)
3385 P = P+SPO/1000.*0.935*(2.22D-3 * P+2.17)
3388 P = P+SPB/1000.*11.35*(16.66D-3 * P+1.33)
3391 P = P+SPO/1000.*0.935*(2.22D-3 * P+2.17)
3394 P = P+SPB/1000.*11.35*(16.66D-3 * P+1.33)
3397 P = P+SPO/1000.*0.935*(2.22D-3 * P+2.17)
3400 P = P+SPB/1000.*11.35*(16.66D-3 * P+1.33)
3404 SW = (503.-467.)/DCOS(THET)
3405 P = P+SW/1000.*19.3*(16.66D-3 * P+1.33)
3409 SCONC = (467.-315.)/DCOS(THET)
3410 P = P+SCONC/1000.*2.5*(2.22D-3*P+2.17)
3413 SC = (315.-90.)/DCOS(THET)
3414 P = P+SC/1000.*1.93*(2.22D-3*P+2.17) ! Carbone
3421 ************************************************************************
3422 SUBROUTINE BRANSON(PXZ,PHI,ALAM,XEA,YEA)
3423 ************************************************************************
3425 * Correction de Branson du multiple scattering dans l'absorbeur
3427 ************************************************************************
3429 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3431 PARAMETER(NPLANE=10)
3432 COMMON/PARAM/ZPLANEP(NPLANE),THICK,XPREC,YPREC,B0,BL3,ZMAGS,
3433 & ZMAGE,ZABS,XMAG,ZBP1,ZBP2,CONST
3436 IF (PXZ.LT.0.) ASIGN = -1.
3442 PTOT = PXZ/DCOS(ALAM)
3446 REA = DSQRT(XEA**2+YEA**2)
3447 IF (REA.GT.26.3611) THEN
3449 ELSE ! Abso. W pour theta < 3 deg
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)
3458 PXZ = DSQRT(PX**2+PZ**2)
3460 ALAM = DATAN2(PY,PXZ)
3462 ** THET = DATAN2(REA,ZEA)
3464 PT = DSQRT(PX**2+PY**2)
3465 THET = DATAN2(PT,PZ)
3466 PTOT = DEDX(PTOT,THET,XEA,YEA)
3468 PXZ = ASIGN*PTOT*DCOS(ALAM)
3473 ***************************************************************
3474 SUBROUTINE DSINV(N,A,IDIM,IFAIL)
3475 ***************************************************************
3477 DOUBLE PRECISION A(IDIM,*), ZERO, ONE, X, Y
3482 DOUBLE PRECISION S1, S31, S32, S33, DOTF
3485 DOTF(X,Y,S1) = X * Y + S1
3487 DATA HNAME / 'DSINV ' /
3488 DATA ZERO, ONE / 0.D0, 1.D0 /
3490 IF(IDIM .LT. N .OR. N .LE. 0) GOTO 900
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
3501 A(J,L) = A(J,J)*A(L,J)
3504 S1 = DOTF(A(L,I),A(I,J+1),S1)
3515 IF(N .EQ. 1) GOTO 399
3517 A(2,1) = A(1,2)*A(2,2)
3518 IF(N .EQ. 2) GOTO 320
3524 S31 = DOTF(A(K,I+1),A(I+1,J),S31)
3527 A(J,K) = -S31*A(J,J)
3529 A(J-1,J) = -A(J-1,J)
3530 A(J,J-1) = A(J-1,J)*A(J,J)
3534 IF(J .EQ. N) GOTO 325
3537 S33 = DOTF(A(J,I),A(I,J),S33)
3545 S32 = DOTF(A(K,I),A(I,J),S32)
3550 IF(J .LT. N) GOTO 323
3554 900 CALL TMPRNT(HNAME,N,IDIM,0)
3558 *******************************************************
3559 SUBROUTINE TMPRNT(NAME,N,IDIM,K)
3560 *******************************************************
3563 LOGICAL MFLAG, RFLAG
3565 IF(NAME(2:2) .EQ. 'S') THEN
3566 CALL KERMTR('F012.1',LGFILE,MFLAG,RFLAG)
3568 CALL KERMTR('F011.1',LGFILE,MFLAG,RFLAG)
3570 IF(NAME(3:6) .EQ. 'FEQN') ASSIGN 1002 TO IFMT
3571 IF(NAME(3:6) .NE. 'FEQN') ASSIGN 1001 TO IFMT
3573 IF(LGFILE .EQ. 0) THEN
3574 IF(NAME(3:6) .EQ. 'FEQN') THEN
3575 WRITE(*,IFMT) NAME, N, IDIM, K
3577 WRITE(*,IFMT) NAME, N, IDIM
3580 IF(NAME(3:6) .EQ. 'FEQN') THEN
3581 WRITE(LGFILE,IFMT) NAME, N, IDIM, K
3583 WRITE(LGFILE,IFMT) NAME, N, IDIM
3587 IF(.NOT. RFLAG) CALL ABEND
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.)
3600 * Revision 1.1.1.1 1996/02/15 17:48:35 mclareni
3605 ***********************************************************
3606 SUBROUTINE KERSET(ERCODE,LGFILE,LIMITM,LIMITR)
3607 ***********************************************************
3609 PARAMETER(KOUNTE = 27)
3610 CHARACTER*6 ERCODE, CODE(KOUNTE)
3611 LOGICAL MFLAG, RFLAG
3612 INTEGER KNTM(KOUNTE), KNTR(KOUNTE)
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 /
3645 IF(ERCODE .NE. ' ') THEN
3647 IF(ERCODE(1:L) .EQ. ERCODE) GOTO 12
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
3658 ENTRY KERMTR(ERCODE,LOG,MFLAG,RFLAG)
3661 IF(ERCODE .EQ. CODE(I)) GOTO 21
3663 WRITE(*,1000) ERCODE
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)
3674 WRITE(LOGF,1001) CODE(I)
3677 IF(MFLAG .AND. RFLAG) THEN
3678 IF(LOGF .LT. 1) THEN
3679 WRITE(*,1002) CODE(I)
3681 WRITE(LOGF,1002) CODE(I)
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 ',
3690 1002 FORMAT(/' ***** CERN LIBRARY ERROR CONDITION ',A6)
3693 **************************
3695 **************************
3700 ************************************************************************
3701 SUBROUTINE FCNFIT(NPAR, GRAD, FVAL, XVAL, IFLAG, FUTIL)
3702 ************************************************************************
3703 * With magnetic Field Map GRKUTA
3705 * Calcule FVAL: la fonction minimisee par MINUIT
3706 * With magnetic field map
3708 ************************************************************************
3709 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3711 PARAMETER(NPLANE=10)
3713 COMMON/PARAM/ZPLANEP(NPLANE),THICK,XPREC,YPREC,B0,BL3,ZMAGS,
3714 & ZMAGE,ZABS,XMAG,ZBP1,ZBP2,CONST
3716 COMMON/MEAS/LPLANE(NPLANE),XMP(NPLANE),YMP(NPLANE)
3718 COMMON/FCNOUT/PXZEA,ALAMEA,PHIEA,XEA,YEA,NPLU,CHI2
3720 DIMENSION GRAD(*),XVAL(*),AMS(500),DISTAZ(500)
3722 DIMENSION XP(NPLANE),YP(NPLANE),
3723 & COV(NPLANE,NPLANE),AP(NPLANE),COVY(NPLANE,NPLANE)
3724 DIMENSION VECT(7),VOUT(7)
3733 IF (XVAL(1).LT.0.) ASIGN = -1.
3736 PXZ = DABS(1./XVAL(1))
3741 PTOT = PXZ/DCOS(ALAM)
3763 R = SQRT(VECT(1)*VECT(1)+VECT(2)*VECT(2))
3770 ** PRINT *,' AV GRKUTA ASIGN',ASIGN,' THET',THET
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.)
3777 ** WRITE(6,*) NSTEP,(VECT(I),I=1,7)
3778 CALL RECO_GRKUTA(ASIGN,STEP,VECT,VOUT)
3783 R = SQRT(VECT(1)*VECT(1)+VECT(2)*VECT(2))
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)
3789 AP(ICH) = (0.0136D0/PTOT)*DSQRT(AL)*(1+0.038D0*DLOG(AL))
3791 ** PRINT *,' AP GRKUTA ASIGN',ASIGN,' THET',THET
3794 ** Matrice de covariance
3797 IF (LPLANE(II).EQ.1) THEN
3802 IF (LPLANE(JJ).EQ.1) THEN
3808 COV(J,I) =COV(J,I) + XPREC**2
3811 * IF (I .EQ. 10 .AND. J .EQ. 10) PRINT *,'10 10 ',COV(J,I)
3814 & + (ZPLANEP(II) + DISTAZ(L))*(ZPLANEP(JJ) +
3815 & DISTAZ(L))*AMS(L)**2
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)
3824 COVY(J,I) = COV(J,I)
3826 COVY(J,I) = COVY(J,I) - XPREC**2 + YPREC**2
3833 * Inversion des matrices de covariance
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'
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'
3846 * PRINT *,(COV(J,I),J=1,NPLANE)
3849 ** Calcul de FVAL ou CHI2
3853 IF (LPLANE(II).EQ.1) THEN
3858 IF (LPLANE(JJ).EQ.1) THEN
3861 FVAL = FVAL + COV(J,I)*(XMP(II)-XP(II))*(XMP(JJ)-XP(JJ))
3862 FVAL = FVAL + COVY(J,I)*(YMP(II)-YP(II))
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
3875 ** IF (CHI2.GT.1.E4) THEN
3876 ** PRINT *,'FCNFIT CHI2= ',CHI2
3881 1000 FORMAT(I5,7F12.6)
3886 ************************************************************************
3887 SUBROUTINE NEWFCNFIT(NPAR, GRAD, FVAL, XVAL, IFLAG, FUTIL)
3888 ************************************************************************
3889 * With magnetic Field Map GRKUTA
3891 * Calcule FVAL: la fonction minimisee par MINUIT
3892 * With magnetic field map
3894 ************************************************************************
3895 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3897 PARAMETER(NPLANE=10)
3899 COMMON/PARAM/ZPLANEP(NPLANE),THICK,XPREC,YPREC,B0,BL3,ZMAGS,
3900 & ZMAGE,ZABS,XMAG,ZBP1,ZBP2,CONST
3902 COMMON/MEAS/LPLANE(NPLANE),XMP(NPLANE),YMP(NPLANE)
3904 COMMON/FCNOUT/PXZEA,ALAMEA,PHIEA,XEA,YEA,NPLU,CHI2
3906 DIMENSION GRAD(*),XVAL(*),AMS(500),DISTAZ(500)
3908 DIMENSION XP(NPLANE),YP(NPLANE),
3909 & COV(NPLANE,NPLANE),AP(NPLANE),COVY(NPLANE,NPLANE)
3910 DIMENSION VECT(7),VOUT(7)
3919 IF (XVAL(1).LT.0.) ASIGN = -1.
3922 PXZ = DABS(1./XVAL(1))
3927 PTOT = PXZ/DCOS(ALAM)
3949 R = SQRT(VECT(1)*VECT(1)+VECT(2)*VECT(2))
3956 ** PRINT *,' AV GRKUTA ASIGN',ASIGN,' THET',THET
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.)
3963 ** WRITE(6,*) NSTEP,(VECT(I),I=1,7)
3964 CALL RECO_GRKUTA (ASIGN,STEP,VECT,VOUT)
3969 R = SQRT(VECT(1)*VECT(1)+VECT(2)*VECT(2))
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)
3975 AP(ICH) = (0.0136D0/PTOT)*DSQRT(AL)*(1+0.038D0*DLOG(AL))
3977 ** PRINT *,' AP GRKUTA ASIGN',ASIGN,' THET',THET
3980 ** Matrice de covariance
3983 IF (LPLANE(II).EQ.1) THEN
3988 IF (LPLANE(JJ).EQ.1) THEN
3994 COV(J,I) =COV(J,I) + XPREC**2
3997 * IF (I .EQ. 10 .AND. J .EQ. 10) PRINT *,'10 10 ',COV(J,I)
4000 & + (ZPLANEP(II) + DISTAZ(L))*(ZPLANEP(JJ) +
4001 & DISTAZ(L))*AMS(L)**2
4005 & + (ZPLANEP(II)-ZPLANEP(K))*
4006 & (ZPLANEP(JJ)-ZPLANEP(K))*AP(K)**2
4007 * IF (I .EQ. 10 .AND. J .EQ. 10) PRINT *,'10 10 ',COV(J,I)
4010 COVY(J,I) = COV(J,I)
4012 COVY(J,I) = COVY(J,I) - XPREC**2 + YPREC**2
4019 * Inversion des matrices de covariance
4023 CALL DSINV(NPLU, COV, NPLANE, IFAIL)
4024 ** IF (JFAIL.NE.0 .AND. IFAIL .NE. 0) STOP 'ERROR'
4025 IF (IFAIL .NE. 0) STOP 'ERROR'
4027 CALL DSINV(NPLU, COVY, NPLANE, IFAIL)
4028 ** IF (JFAIL.NE.0 .AND. IFAIL .NE. 0) STOP 'ERROR'
4029 IF (IFAIL .NE. 0) STOP 'ERROR'
4030 * PRINT *,' COVARIANCE MATRIX AFTER'
4032 * PRINT *,(COV(J,I),J=1,NPLANE)
4035 ** Calcul de FVAL ou CHI2
4039 IF (LPLANE(II).EQ.1) THEN
4044 IF (LPLANE(JJ).EQ.1) THEN
4047 FVAL = FVAL + COV(J,I)*(XMP(II)-XP(II))*(XMP(JJ)-XP(JJ))
4048 FVAL = FVAL + COVY(J,I)*(YMP(II)-YP(II))
4050 ** IF (JJ.EQ.II) THEN
4051 ** FVAL = FVAL + (XM(II)-XP(II))*(XM(JJ)-XP(JJ))/XPREC**2
4052 ** FVAL = FVAL + (YM(II)-YP(II))
4053 ** & *(YM(JJ)-YP(JJ))/YPREC**2
4061 ** IF (CHI2.GT.1.E4) THEN
4062 ** PRINT *,'FCNFIT CHI2= ',CHI2
4067 1000 FORMAT(I5,7F12.6)
4072 ***********************************************************************
4073 SUBROUTINE INITFIELDOLD
4076 ***********************************************************************
4078 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4080 ** IMPLICIT REAL*8(A-H,O-Z)
4082 COMMON/DAT1/Z(81),X(81),Y(81,44),DX,DZ,LPX,LPY,LPZ
4083 COMMON/DAT2/BX(81,81,44),BY(81,81,44),BZ(81,81,44)
4084 COMMON/REG1/ZMAX,ZMIN,XMAX,XMIN
4085 COMMON/REG2/AY1,CY1,AY2,CY2
4086 ** REAL *4 BXP,BYP,BZP
4087 COMMON/SDAT1/ZP(51),RAD(10),FI(33),DZP,DFI,DR,YY0,LPPZ,NR,NFI
4088 COMMON/SDAT2/BXP(51,10,33),BYP(51,10,33),BZP(51,10,33)
4089 COMMON/SDAT4/B(2,2,32)
4090 COMMON/REG3/ZPMAX,ZPMIN,RMAX,RMIN
4091 cc COMMON/CONST/PI2,EPS
4093 1000 FORMAT(5(1X,D15.7))
4094 2000 FORMAT(5(1X,I5))
4095 READ(40,2000) LPX,LPY,LPZ
4096 READ(40,1000) (Z(K),K=1,81)
4097 READ(40,1000) (X(K),K=1,81)
4098 READ(40,1000) DX,DY,DZ
4099 READ(40,1000) ZMAX,ZMIN,XMAX,XMIN
4100 c write(*,*) 'zmin zmax',ZMIN,ZMAX
4101 c write(*,*) 'xmin xmax',XMIN,XMAX
4102 READ(40,1000) AY1,CY1,AY2,CY2
4103 c write(*,*) 'ay1,cy1,ay2,cy2', AY1,CY1,AY2,CY2
4104 cc READ(40,1000) PI2,EPS
4105 READ(40,1000) (((BX(K,L,M),K=1,81),L=1,81),M=1,44)
4106 READ(40,1000) (((BY(K,L,M),K=1,81),L=1,81),M=1,44)
4107 READ(40,1000) (((BZ(K,L,M),K=1,81),L=1,81),M=1,44)
4111 READ(40,2000) LPPZ,NR,NFI
4112 READ(40,1000) (ZP(K),K=1,51)
4113 READ(40,1000) (RAD(K),K=1,10)
4114 READ(40,1000) (FI(L),L=1,33)
4115 READ(40,1000) DZP,DFI,DR
4116 c write(*,*) 'dzp dfi dR',DZP,DFI,DR
4117 READ(40,1000) ZPMAX,ZPMIN,RMAX,RMIN
4118 c write(*,*) 'zmin zmax',ZPMIN,ZPMAX
4119 c write(*,*) 'Rmin Rmax',RMIN,RMAX
4120 READ(40,1000) (((BXP(K,L,M),K=1,51),L=1,10),M=1,33)
4121 READ(40,1000) (((BYP(K,L,M),K=1,51),L=1,10),M=1,33)
4122 READ(40,1000) (((BZP(K,L,M),K=1,51),L=1,10),M=1,33)
4123 READ(40,1000) (((B(K,L,M),K=1,2),L=1,2),M=1,32)
4133 ***********************************************************************
4134 SUBROUTINE RECO_GUFLDOLD(X,F)
4135 C ^^^^^^^^^^^^^^^^^^^^^^
4136 C field map G. Chabratova
4138 C Field map 31/05/99
4139 ***********************************************************************
4142 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4143 COMMON/MAGERR/IMAGERR
4155 CALL FREG1(Z0,X0,Y0,FZ0,FX0,FY0,IND)
4156 ** PRINT 3000,Z0,X0,Y0,FZ0,FX0,FY0,IND
4159 CALL FREG2(Z0,X0,Y0,FZ0,FX0,FY0,IND)
4164 ** PRINT 3000,Z0,X0,Y0,FZ0,FX0,FY0,IND
4166 1000 format(1x,'Attention!!! The point is out of range!!!')
4168 3000 FORMAT(1X,'Z=',D13.7,1X,'X=',D13.7,1X,'Y=',D13.7,1X,
4169 & 'BZ=',D13.7,1X,'BX=',D13.7,1X,'BY=',D13.7,1X,'IND=',I3)
4192 **************************************************
4193 SUBROUTINE FREG1(Z0,X0,Y0,FZ0,FX0,FY0,IND)
4194 **************************************************
4195 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4197 COMMON/DAT1/Z(81),X(81),Y(81,44),DX,DZ,LPX,LPY,LPZ
4198 COMMON/DAT2/BX(81,81,44),BY(81,81,44),BZ(81,81,44)
4199 COMMON/REG1/ZMAX,ZMIN,XMAX,XMIN
4200 COMMON/REG2/AY1,CY1,AY2,CY2
4205 IF(Z0.LT.ZMIN.OR.Z0.GT.ZMAX)GO TO 100
4206 IF(X0.LT.XMIN.OR.X0.GT.XMAX)GO TO 100
4209 2000 FORMAT(1X,'YY1=',D15.7,1X,'YY2=',D15.7,1X,'Y0=',D15.7)
4210 c PRINT 2000,YY1,YY2,Y0
4211 IF(Y0.LT.YY1)GO TO 100
4212 IF(Y0.GT.YY2)GO TO 100
4213 CALL FIZ(Z0,Z,DZ,KC,K0,Z1,Z2,81)
4214 CALL FIZ(X0,X,DX,LC,L0,X1,X2,81)
4215 DY=(YY2-YY1)/DFLOAT(LPY-1)
4219 IF(Y0.GE.(YY1+DFLOAT(M0)*DY).AND.Y0.LE.(YY1+DFLOAT(M0+1)*DY))
4223 Y2=(Y0-(YY1+DFLOAT(M0)*DY))/DY
4225 ** write(*,*) 'm0 Y1 Y2',m0,Y1,Y2
4226 ** print *,' k0=',k0,' l0=',l0,' m0=',m0
4227 ** print *,' z1=',z1,' z2=',z2
4228 FX1=Z1*BX(K0,L0,M0)+Z2*BX(K0+1,L0,M0)
4229 FX2=Z2*BX(K0+1,L0+1,M0)+Z1*BX(K0,L0+1,M0)
4231 GX1=Z1*BX(K0,L0,M0+1)+Z2*BX(K0+1,L0,M0+1)
4232 GX2=Z2*BX(K0+1,L0+1,M0+1)+Z1*BX(K0,L0+1,M0+1)
4235 FX1=Z1*BY(K0,L0,M0)+Z2*BY(K0+1,L0,M0)
4236 FX2=Z2*BY(K0+1,L0+1,M0)+Z1*BY(K0,L0+1,M0)
4238 GX1=Z1*BY(K0,L0,M0+1)+Z2*BY(K0+1,L0,M0+1)
4239 GX2=Z2*BY(K0+1,L0+1,M0+1)+Z1*BY(K0,L0+1,M0+1)
4242 FX1=Z1*BZ(K0,L0,M0)+Z2*BZ(K0+1,L0,M0)
4243 FX2=Z2*BZ(K0+1,L0+1,M0)+Z1*BZ(K0,L0+1,M0)
4245 GX1=Z1*BZ(K0,L0,M0+1)+Z2*BZ(K0+1,L0,M0+1)
4246 GX2=Z2*BZ(K0+1,L0+1,M0+1)+Z1*BZ(K0,L0+1,M0+1)
4252 1000 format(1x,'Attention!!! The point is out of range!!!')
4257 *************************************************
4258 SUBROUTINE FIZ(Z0,Z,DEL,KI,K0,Z1,Z2,NDZ)
4259 *************************************************
4260 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4261 ** CC DIMENSION Z(NDZ)
4267 * if (k0.gt.81) print*,'ndz=',ndz
4268 IF (KJ.GT.NDZ) THEN ! CCC
4272 DO 1 K=KJ,NDZ-1 ! CCC
4279 * print *,'K0=',K0,' Z0',z0, Z(K0), Z(K0+1),z2
4280 if (k0.gt.81) print*,'k0=',k0
4281 Z2=(Z0-Z(K0))/(Z(K0+1)-Z(K0))
4283 ** write(*,*) 'ko z1 z2', K0,Z1,Z2,' ki=',ki,' kj=',kj,' K=',K
4284 ** write(*,*)' NDZ Z(K0) Z(K0+1)',NDZ,Z(K0), Z(K0+1)
4288 ***************************************************
4289 SUBROUTINE FREG2(Z0,X0,Y0,FZ0,FX0,FY0,IND)
4290 **************************************************
4291 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4292 ** REAL *4 BXP,BYP,BZP
4293 COMMON/SDAT1/ZP(51),RAD(10),FI(33),DZP,DFI,DR,YY0,LPPZ,NR,NFI
4294 COMMON/SDAT2/BXP(51,10,33),BYP(51,10,33),BZP(51,10,33)
4295 COMMON/SDAT4/B(2,2,32)
4296 COMMON/REG3/ZPMAX,ZPMIN,RMAX,RMIN
4297 cc COMMON/CONST/PI2,EPS
4304 R0=DSQRT((X0-YY0)**2+Y0**2)
4305 c write (*,*)'ro=',R0
4306 IF(Z0.LT.ZPMIN.OR.Z0.GT.ZPMAX)GO TO 100
4307 IF(R0.LT.RMIN.OR.R0.GT.RMAX)GO TO 100
4308 IF(R0.LE.DR)GO TO 3000
4309 CALL FIZ(Z0,ZP,DZP,KP,K0,Z1,Z2,51)
4310 CALL FIZ(R0,RAD,DR,LP,L0,X1,X2,10)
4311 ** print *,' r0=',r0,' rad=',rad,' dr=',dr,' lp=',lp,' l0=',l0,
4312 ** & ' x1=',x1,' x2=',x2
4313 FI0=DACOS((X0-YY0)/R0)
4314 IF(Y0.LT.0.D+0)FI0=PI2-FI0
4315 CALL FIZ(FI0,FI,DFI,MP,M0,Y1,Y2,32)
4316 ** print *,' Apres FIZ',' k0=',k0,' l0=',l0,' m0=',m0
4317 FX1=Z1*BXP(K0,L0,M0)+Z2*BXP(K0+1,L0,M0)
4318 FX2=Z2*BXP(K0+1,L0+1,M0)+Z1*BXP(K0,L0+1,M0)
4320 GX1=Z1*BXP(K0,L0,M0+1)+Z2*BXP(K0+1,L0,M0+1)
4321 GX2=Z2*BXP(K0+1,L0+1,M0+1)+Z1*BXP(K0,L0+1,M0+1)
4324 FX1=Z1*BYP(K0,L0,M0)+Z2*BYP(K0+1,L0,M0)
4325 FX2=Z2*BYP(K0+1,L0+1,M0)+Z1*BYP(K0,L0+1,M0)
4327 GX1=Z1*BYP(K0,L0,M0+1)+Z2*BYP(K0+1,L0,M0+1)
4328 GX2=Z2*BYP(K0+1,L0+1,M0+1)+Z1*BYP(K0,L0+1,M0+1)
4331 FX1=Z1*BZP(K0,L0,M0)+Z2*BZP(K0+1,L0,M0)
4332 ** CCC FX2=Z2*BZ(K0+1,L0+1,M0)+Z1*BZ(K0,L0+1,M0)
4333 FX2=Z2*BZP(K0+1,L0+1,M0)+Z1*BZP(K0,L0+1,M0)
4335 GX1=Z1*BZP(K0,L0,M0+1)+Z2*BZP(K0+1,L0,M0+1)
4336 GX2=Z2*BZP(K0+1,L0+1,M0+1)+Z1*BZP(K0,L0+1,M0+1)
4343 1000 format(1x,'Attention!!! The point is out of range!!!')
4347 IF(R0.LT.EPS)GO TO 4000
4348 CALL FIZ(Z0,ZP,DZP,KP,K0,Z1,Z2,51)
4351 IF(Y0.LT.0.D+0)FI0=PI2-FI0
4352 CALL FIZ(FI0,FI,DFI,MP,M0,Y1,Y2,32)
4353 ALF2=B(1,1,M0)*XX+B(1,2,M0)*Y0
4354 ALF3=B(2,1,M0)*XX+B(2,2,M0)*Y0
4356 FX1=ALF1*BXP(K0,1,1)+ALF2*BXP(K0,1,M0)+ALF3*BXP(K0,1,M0+1)
4357 FX2=ALF1*BXP(K0+1,1,1)+ALF2*BXP(K0+1,1,M0)+ALF3*BXP(K0+1,1,M0+1)
4359 FX1=ALF1*BYP(K0,1,1)+ALF2*BYP(K0,1,M0)+ALF3*BYP(K0,1,M0+1)
4360 FX2=ALF1*BYP(K0+1,1,1)+ALF2*BYP(K0+1,1,M0)+ALF3*BYP(K0+1,1,M0+1)
4362 FX1=ALF1*BZP(K0,1,1)+ALF2*BZP(K0,1,M0)+ALF3*BZP(K0,1,M0+1)
4363 FX2=ALF1*BZP(K0+1,1,1)+ALF2*BZP(K0+1,1,M0)+ALF3*BZP(K0+1,1,M0+1)
4365 c write(*,*) 'R<Dr:B(1,1,m0) B(1,2,m0) B(2,1,m0) B(2,2,m0)',
4366 c +B(1,1,M0),B(1,2,M0),B(2,1,M0),B(2,2,M0)
4367 c write(*,*)'BX(K0,1,1) BX(K0,1,M0)','BX(K0,1,M0+1) BX(K0+1,1,1)'
4368 c + ,'BX(K0+1,1,M0) BX(K0+1,1,M0+1)', BX(K0,1,1),BX(K0,1,M0) ,
4369 c + BX(K0,1,M0+1),BX(K0+1,1,1),BX(K0+1,1,M0),BX(K0+1,1,M0+1)
4370 c write(*,*)'By(K0,1,1) By(K0,1,M0)','By(K0,1,M0+1) By(K0+1,1,1)'
4371 c + ,'By(K0+1,1,M0) By(K0+1,1,M0+1)', By(K0,1,1),By(K0,1,M0) ,
4372 c + By(K0,1,M0+1),By(K0+1,1,1),By(K0+1,1,M0),By(K0+1,1,M0+1)
4373 ccc write (*,*)'Bz(K0,1,1) Bz(K0,1,M0) Bz(K0,1,M0+1) Bz(K0+1,1,1)
4374 77 FORMAT(5x,E15.7,2x,E15.7)
4376 ** 70 FORMAT(22hBz(K0,1,1) Bz(K0,1,M0))
4377 cc PRINT 77 , BzP(K0,1,1),BzP(K0,1,M0)
4379 ** 71 FORMAT(26hBz(K0,1,M0+1) Bz(K0+1,1,1))
4380 cc PRINT 77, BzP(K0,1,M0+1),BzP(K0+1,1,1)
4382 ** 72 FORMAT(29hBz(K0+1,1,M0) Bz(K0+1,1,M0+1))
4383 cc PRINT 77 ,BzP(K0+1,1,M0),BzP(K0+1,1,M0+1)
4384 cc 77 FORMAT(5x,D15.7,5x,D15.7)
4390 CALL FIZ(Z0,ZP,DZP,KP,K0,Z1,Z2,51)
4391 FX0=Z1*BXP(K0,1,1)+Z2*BXP(K0+1,1,1)
4392 FY0=Z1*BYP(K0,1,1)+Z2*BYP(K0+1,1,1)
4393 FZ0=Z1*BZP(K0,1,1)+Z2*BZP(K0+1,1,1)
4394 c write(*,*) ' R<eps: Bx(k) Bx(k+1) By(k) By(k+1) Bz(k) Bz(k+1)',
4395 c + BXP(K0,1,1),BXP(K0+1,1,1),BYP(K0,1,1),BYP(K0+1,1,1),BZP(K0,1,1),
4401 ***********************************************************************
4403 SUBROUTINE RECO_GRKUTA (CHARGE,STEP,VECT,VOUT)
4405 C. ******************************************************************
4407 C. * Runge-Kutta method for tracking a particle through a magnetic *
4408 C. * field. Uses Nystroem algorithm (See Handbook Nat. Bur. of *
4409 C. * Standards, procedure 25.5.20) *
4411 C. * Input parameters *
4412 C. * CHARGE Particle charge *
4413 C. * STEP Step size *
4414 C. * VECT Initial co-ords,direction cosines,momentum *
4415 C. * Output parameters *
4416 C. * VOUT Output co-ords,direction cosines,momentum *
4417 C. * User routine called *
4418 C. * CALL GUFLD(X,F) *
4420 C. * ==>Called by : <USER>, GUSWIM *
4421 C. * Authors R.Brun, M.Hansroul ********* *
4422 C. * V.Perevoztchikov (CUT STEP implementation) *
4425 C. ******************************************************************
4427 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4429 ** REAL CHARGE, STEP, VECT(*), VOUT(*), F(4)
4430 ** REAL XYZT(3), XYZ(3), X, Y, Z, XT, YT, ZT
4431 DIMENSION VECT(*), VOUT(*), F(3)
4432 DIMENSION XYZT(3), XYZ(3)
4433 DIMENSION SECXS(4),SECYS(4),SECZS(4),HXP(3)
4434 EQUIVALENCE (X,XYZ(1)),(Y,XYZ(2)),(Z,XYZ(3)),
4435 + (XT,XYZT(1)),(YT,XYZT(2)),(ZT,XYZT(3))
4437 PARAMETER (MAXIT = 1992, MAXCUT = 11)
4438 PARAMETER (EC=2.9979251D-4,DLT=1D-4,DLT32=DLT/32)
4439 PARAMETER (ZERO=0, ONE=1, TWO=2, THREE=3)
4440 PARAMETER (THIRD=ONE/THREE, HALF=ONE/TWO)
4441 PARAMETER (PISQUA=.986960440109D+01)
4442 PARAMETER (IX=1,IY=2,IZ=3,IPX=4,IPY=5,IPZ=6)
4444 *. ------------------------------------------------------------------
4446 * This constant is for units CM,GEV/C and KGAUSS
4453 PINV = EC * CHARGE / VECT(7)
4459 IF (ABS(H).GT.ABS(REST)) H = REST
4460 CALL RECO_GUFLD(VOUT,F)
4462 * Start of integration
4475 SECXS(1) = (B * F(3) - C * F(2)) * PH2
4476 SECYS(1) = (C * F(1) - A * F(3)) * PH2
4477 SECZS(1) = (A * F(2) - B * F(1)) * PH2
4478 ANG2 = (SECXS(1)**2 + SECYS(1)**2 + SECZS(1)**2)
4479 IF (ANG2.GT.PISQUA) GO TO 40
4480 DXT = H2 * A + H4 * SECXS(1)
4481 DYT = H2 * B + H4 * SECYS(1)
4482 DZT = H2 * C + H4 * SECZS(1)
4487 * Second intermediate point
4489 EST = ABS(DXT)+ABS(DYT)+ABS(DZT)
4490 IF (EST.GT.H) GO TO 30
4492 CALL RECO_GUFLD(XYZT,F)
4497 SECXS(2) = (BT * F(3) - CT * F(2)) * PH2
4498 SECYS(2) = (CT * F(1) - AT * F(3)) * PH2
4499 SECZS(2) = (AT * F(2) - BT * F(1)) * PH2
4503 SECXS(3) = (BT * F(3) - CT * F(2)) * PH2
4504 SECYS(3) = (CT * F(1) - AT * F(3)) * PH2
4505 SECZS(3) = (AT * F(2) - BT * F(1)) * PH2
4506 DXT = H * (A + SECXS(3))
4507 DYT = H * (B + SECYS(3))
4508 DZT = H * (C + SECZS(3))
4512 AT = A + TWO*SECXS(3)
4513 BT = B + TWO*SECYS(3)
4514 CT = C + TWO*SECZS(3)
4516 EST = ABS(DXT)+ABS(DYT)+ABS(DZT)
4517 IF (EST.GT.2.*ABS(H)) GO TO 30
4519 CALL RECO_GUFLD(XYZT,F)
4521 Z = Z + (C + (SECZS(1) + SECZS(2) + SECZS(3)) * THIRD) * H
4522 Y = Y + (B + (SECYS(1) + SECYS(2) + SECYS(3)) * THIRD) * H
4523 X = X + (A + (SECXS(1) + SECXS(2) + SECXS(3)) * THIRD) * H
4525 SECXS(4) = (BT*F(3) - CT*F(2))* PH2
4526 SECYS(4) = (CT*F(1) - AT*F(3))* PH2
4527 SECZS(4) = (AT*F(2) - BT*F(1))* PH2
4528 A = A+(SECXS(1)+SECXS(4)+TWO * (SECXS(2)+SECXS(3))) * THIRD
4529 B = B+(SECYS(1)+SECYS(4)+TWO * (SECYS(2)+SECYS(3))) * THIRD
4530 C = C+(SECZS(1)+SECZS(4)+TWO * (SECZS(2)+SECZS(3))) * THIRD
4532 EST = ABS(SECXS(1)+SECXS(4) - (SECXS(2)+SECXS(3)))
4533 ++ ABS(SECYS(1)+SECYS(4) - (SECYS(2)+SECYS(3)))
4534 ++ ABS(SECZS(1)+SECZS(4) - (SECZS(2)+SECZS(3)))
4536 IF (EST.GT.DLT .AND. ABS(H).GT.1.E-4) GO TO 30
4539 * If too many iterations, go to HELIX
4540 IF (ITER.GT.MAXIT) GO TO 40
4543 IF (EST.LT.(DLT32)) THEN
4546 CBA = ONE/ SQRT(A*A + B*B + C*C)
4554 IF (STEP.LT.0.) REST = -REST
4555 IF (REST .GT. 1.E-5*ABS(STEP)) GO TO 20
4561 * If too many cuts , go to HELIX
4562 IF (NCUT.GT.MAXCUT) GO TO 40
4566 ** ANGLE TOO BIG, USE HELIX
4570 F4 = SQRT(F1**2+F2**2+F3**2)
4579 HXP(1) = F2*VECT(IPZ) - F3*VECT(IPY)
4580 HXP(2) = F3*VECT(IPX) - F1*VECT(IPZ)
4581 HXP(3) = F1*VECT(IPY) - F2*VECT(IPX)
4583 HP = F1*VECT(IPX) + F2*VECT(IPY) + F3*VECT(IPZ)
4587 COST = TWO*SIN(HALF*TET)**2
4591 G3 = (TET-SINT) * HP*RHO1
4596 VOUT(IX) = VECT(IX) + (G1*VECT(IPX) + G2*HXP(1) + G3*F1)
4597 VOUT(IY) = VECT(IY) + (G1*VECT(IPY) + G2*HXP(2) + G3*F2)
4598 VOUT(IZ) = VECT(IZ) + (G1*VECT(IPZ) + G2*HXP(3) + G3*F3)
4600 VOUT(IPX) = VECT(IPX) + (G4*VECT(IPX) + G5*HXP(1) + G6*F1)
4601 VOUT(IPY) = VECT(IPY) + (G4*VECT(IPY) + G5*HXP(2) + G6*F2)
4602 VOUT(IPZ) = VECT(IPZ) + (G4*VECT(IPZ) + G5*HXP(3) + G6*F3)
4605 VOUT(IX) = VECT(IX) + STEP*VECT(IPX)
4606 VOUT(IY) = VECT(IY) + STEP*VECT(IPY)
4607 VOUT(IZ) = VECT(IZ) + STEP*VECT(IPZ)
4614 *******************************************************************
4615 SUBROUTINE RECO_GUFLDOLD1(X,B)
4619 C *** ROUTINE DESCRIBING THE MAGNETIC FIELD IN THE ALICE SETUP ***
4620 C *** NVE 14-NOV-1990 CERN GENEVA ***
4623 C ORIGIN : NICK VAN EIJNDHOVEN
4627 C X = (X,Y,Z) coordinates in cm
4631 C B = Magnetic field components (BX,BY,BZ) in KG
4634 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4636 PARAMETER(NBSTATION=5)
4638 COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
4652 IF (X(3).LT.(-1.*ZCOIL)) THEN
4654 ELSEIF ( X(3).LT.(-1.*ZMAGEND)) THEN
4665 ** print *,' x =',X(1),X(2),X(3)
4666 ** print *,' B =',B(1),B(2),B(3)
4673 *******************************************************************
4674 SUBROUTINE RECO_GUFLD(X,B)
4678 C *** ROUTINE DESCRIBING THE MAGNETIC FIELD IN THE ALICE SETUP ***
4679 C *** NVE 14-NOV-1990 CERN GENEVA ***
4682 C ORIGIN : NICK VAN EIJNDHOVEN
4686 C X = (X,Y,Z) coordinates in cm
4690 C B = Magnetic field components (BX,BY,BZ) in KG
4693 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4696 C --- Common containing magnetic field map data
4697 REAL DZ,DX,DY,UDX,UDY,UDZ
4698 $,XMBEG,YMBEG,ZMBEG,XMEND,YMEND,ZMEND
4702 PARAMETER(MAXFLD=250000)
4703 COMMON /SCXMFD/ NX,NY,NZ,DZ,DX,DY,UDX,UDY,UDZ
4704 $,XMBEG,YMBEG,ZMBEG,XMEND,YMEND,ZMEND
4710 DOUBLE PRECISION ONE, RATX, RATY, RATZ, HIX, HIY, HIZ
4711 $, RATX1, RATY1, RATZ1
4712 $, BHYHZ, BHYLZ, BLYHZ, BLYLZ, BHZ, BLZ
4726 ** BX(JX,JY,JZ)=BV(3*((JZ-1)*(NX*NY)+(JY-1)*NX+(JX-1))+1)
4727 ** BY(JX,JY,JZ)=BV(3*((JZ-1)*(NX*NY)+(JY-1)*NX+(JX-1))+2)
4728 ** BZ(JX,JY,JZ)=BV(3*((JZ-1)*(NX*NY)+(JY-1)*NX+(JX-1))+3)
4733 * --- Act accordingly to ISXFMAP
4735 IF(ISXFMAP.EQ.1) THEN
4736 IF (ABS(X(3)) .LT. 700.
4737 + .AND. X(1)**2+(X(2)+30.)**2 .LT. 560.**2 ) THEN
4741 ELSE IF (X(3) .GE. 725. .AND. X(3) .LT. 1225.) THEN
4742 DZ=ABS(975.-X(3))/100.
4743 B(1)=(1.-0.1*DZ*DZ)*7.0
4751 ELSE IF(ISXFMAP.LE.3) THEN
4752 IF (-700.LT.X(3).AND.X(3).LT.ZMBEG
4753 + .AND. X(1)**2+(X(2)+30.)**2 .LT. 560.**2 ) THEN
4757 ELSE IF ((X(3) .GE. ZMBEG .AND. X(3) .LT. ZMEND) .AND.
4758 + (XMBEG.LE.ABS(X(1)).AND.ABS(X(1)).LT.XMEND) .AND.
4759 + (YMBEG.LE.ABS(X(2)).AND.ABS(X(2)).LT.YMEND)) THEN
4763 C --- find the position in the grid ---
4765 XL(1)=ABS(X(1))-XMBEG
4766 XL(2)=ABS(X(2))-YMBEG
4783 IF(ISXFMAP.EQ.2) THEN
4784 * ... Simple interpolation
4786 B(1) = BX(IX,IY,IZ)*(ONE-RATX) + BX(IX+1,IY+1,IZ+1)*RATX
4787 B(2) = BY(IX,IY,IZ)*(ONE-RATY) + BY(IX+1,IY+1,IZ+1)*RATY
4788 B(3) = BZ(IX,IY,IZ)*(ONE-RATZ) + BZ(IX+1,IY+1,IZ+1)*RATZ
4789 ELSE IF(ISXFMAP.EQ.3) THEN
4790 * ... more complicated interpolation
4794 ** print *,' bx by bz', BX(IX ,IY+1,IZ+1),BY(IX ,IY+1,IZ+1)
4795 ** & , BZ(IX ,IY+1,IZ+1)
4796 BHYHZ = BX(IX ,IY+1,IZ+1)*RATX1+BX(IX+1,IY+1,IZ+1)*RATX
4797 BHYLZ = BX(IX ,IY+1,IZ )*RATX1+BX(IX+1,IY+1,IZ )*RATX
4798 BLYHZ = BX(IX ,IY ,IZ+1)*RATX1+BX(IX+1,IY ,IZ+1)*RATX
4799 BLYLZ = BX(IX ,IY ,IZ )*RATX1+BX(IX+1,IY ,IZ )*RATX
4800 BHZ = BLYHZ *RATY1+BHYHZ *RATY
4801 BLZ = BLYLZ *RATY1+BHYLZ *RATY
4802 B(1) = BLZ *RATZ1+BHZ *RATZ
4804 BHYHZ = BY(IX ,IY+1,IZ+1)*RATX1+BY(IX+1,IY+1,IZ+1)*RATX
4805 BHYLZ = BY(IX ,IY+1,IZ )*RATX1+BY(IX+1,IY+1,IZ )*RATX
4806 BLYHZ = BY(IX ,IY ,IZ+1)*RATX1+BY(IX+1,IY ,IZ+1)*RATX
4807 BLYLZ = BY(IX ,IY ,IZ )*RATX1+BY(IX+1,IY ,IZ )*RATX
4808 BHZ = BLYHZ *RATY1+BHYHZ *RATY
4809 BLZ = BLYLZ *RATY1+BHYLZ *RATY
4810 B(2) = BLZ *RATZ1+BHZ *RATZ
4812 BHYHZ = BZ(IX ,IY+1,IZ+1)*RATX1+BZ(IX+1,IY+1,IZ+1)*RATX
4813 BHYLZ = BZ(IX ,IY+1,IZ )*RATX1+BZ(IX+1,IY+1,IZ )*RATX
4814 BLYHZ = BZ(IX ,IY ,IZ+1)*RATX1+BZ(IX+1,IY ,IZ+1)*RATX
4815 BLYLZ = BZ(IX ,IY ,IZ )*RATX1+BZ(IX+1,IY ,IZ )*RATX
4816 BHZ = BLYHZ *RATY1+BHYHZ *RATY
4817 BLZ = BLYLZ *RATY1+BHYLZ *RATY
4818 B(3) = BLZ *RATZ1+BHZ *RATZ
4821 * ... use the dipole symmetry
4823 IF (X(1)*X(2).LT.0) B(2)=-B(2)
4824 IF (X(1).LT.0) B(3)=-B(3)
4834 ENDIF ! z-coord for m.f.
4836 ENDIF ! endif ISXFMAP
4843 BTOT=SQRT(B(1)**2+B(2)**2+B(3)**2)
4844 IF(BTOT.GT.SXMGMX) THEN
4845 PRINT 10100, BTOT,SXMGMX
4846 10100 FORMAT(' *GUFLD* Field ',G10.4,' larger than max ',G10.4)
4864 *******************************************
4865 DOUBLE PRECISION FUNCTION BX(JX,JY,JZ)
4867 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4869 C --- Common containing magnetic field map data
4870 REAL DZ,DX,DY,UDX,UDY,UDZ
4871 $,XMBEG,YMBEG,ZMBEG,XMEND,YMEND,ZMEND
4875 PARAMETER(MAXFLD=250000)
4876 COMMON /SCXMFD/ NX,NY,NZ,DZ,DX,DY,UDX,UDY,UDZ
4877 $,XMBEG,YMBEG,ZMBEG,XMEND,YMEND,ZMEND
4880 BX=BV(3*((JZ-1)*(NX*NY)+(JY-1)*NX+(JX-1))+1)
4884 *******************************************
4885 DOUBLE PRECISION FUNCTION BY(JX,JY,JZ)
4887 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4889 C --- Common containing magnetic field map data
4890 REAL DZ,DX,DY,UDX,UDY,UDZ
4891 $,XMBEG,YMBEG,ZMBEG,XMEND,YMEND,ZMEND
4895 PARAMETER(MAXFLD=250000)
4896 COMMON /SCXMFD/ NX,NY,NZ,DZ,DX,DY,UDX,UDY,UDZ
4897 $,XMBEG,YMBEG,ZMBEG,XMEND,YMEND,ZMEND
4900 BY=BV(3*((JZ-1)*(NX*NY)+(JY-1)*NX+(JX-1))+2)
4904 *******************************************
4905 DOUBLE PRECISION FUNCTION BZ(JX,JY,JZ)
4907 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4909 C --- Common containing magnetic field map data
4910 REAL DZ,DX,DY,UDX,UDY,UDZ
4911 $,XMBEG,YMBEG,ZMBEG,XMEND,YMEND,ZMEND
4915 PARAMETER(MAXFLD=250000)
4916 COMMON /SCXMFD/ NX,NY,NZ,DZ,DX,DY,UDX,UDY,UDZ
4917 $,XMBEG,YMBEG,ZMBEG,XMEND,YMEND,ZMEND
4920 BZ=BV(3*((JZ-1)*(NX*NY)+(JY-1)*NX+(JX-1))+3)
4924 ***********************************************************************
4925 SUBROUTINE INITFIELD
4930 C *** INITIALISATION OF THE FIELD MAP ***
4931 C *** FCA 24-AUG-1998 CERN GENEVA ***
4933 C CALLED BY : GALICE
4934 C ORIGIN : NICK VAN EIJNDHOVEN
4936 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4939 C --- Common containing magnetic field map data
4940 REAL DZ,DX,DY,UDX,UDY,UDZ
4941 $,XMBEG,YMBEG,ZMBEG,XMEND,YMEND,ZMEND
4945 PARAMETER(MAXFLD=250000)
4946 COMMON /SCXMFD/ NX,NY,NZ,DZ,DX,DY,UDX,UDY,UDZ
4947 $,XMBEG,YMBEG,ZMBEG,XMEND,YMEND,ZMEND
4952 IF(ISXFMAP.EQ.1) THEN
4953 * ... constant field, nothing to do
4954 ELSE IF(ISXFMAP.LE.3) THEN
4955 * ... constant mesh field
4956 PRINT 10000, ISXFMAP
4957 10000 FORMAT(' *SXFMAP* Magnetic field map flag: ',I5
4958 $ ,'; Reading magnetic field map data ')
4960 OPEN(77,FILE='data/field01.dat',
4961 $ FORM='FORMATTED',STATUS='OLD')
4962 READ(77,*) NX,NY,NZ,DX,DY,DZ,XMBEG,YMBEG,ZMBEG
4963 PRINT*,'NX,NY,NZ,DX,DY,DZ,XMBEG,YMBEG,ZMBEG',
4964 $ NX,NY,NZ,DX,DY,DZ,XMBEG,YMBEG,ZMBEG
4965 IF(3*NX*NY*NZ.GT.MAXFLD) THEN
4966 WRITE(6,10100) 3*NX*NY*NZ,MAXFLD
4967 STOP 'Increase MAXFLD'
4972 XMEND=XMBEG+(NX-1)*DX
4973 YMEND=YMBEG+(NY-1)*DY
4974 ZMEND=ZMBEG+(NZ-1)*DZ
4976 IPZ=3*(IZ-1)*(NX*NY)
4981 READ(77,*) BV(IPX+3),BV(IPX+2),BV(IPX+1)
4986 ENDIF ! endif ISXFMAP
4988 10100 FORMAT('*** SXFMAP: Need ',I7,' have ',I7)
4991 ****************************************
4992 SUBROUTINE RANNOR (A,B)
4993 ****************************************
4995 C CERN PROGLIB# V100 RANNOR .VERSION KERNFOR 4.18 880425
4999 IF (Y.EQ.0.) Y = RNDM()
5003 R = SQRT (-2.0*LOG(Y))
5008 ************************************************************************
5009 SUBROUTINE OLDFCNFIT(NPAR,GRAD,FVAL,XVAL,IFLAG,FUTIL)
5010 ************************************************************************
5011 * Calcule FVAL: la fonction minimisee par MINUIT
5012 * with constant magnetic Field
5014 ************************************************************************
5016 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
5018 PARAMETER(NPLANE=10)
5020 COMMON/PARAM/ZPLANEP(NPLANE),THICK,XPREC,YPREC,B0,BL3,ZMAGS,
5021 & ZMAGE,ZABS,XMAG,ZBP1,ZBP2,CONST
5023 COMMON/MEAS/LPLANE(NPLANE),XMP(NPLANE),YMP(NPLANE)
5025 COMMON/FCNOUT/PXZEA,ALAMEA,PHIEA,XEA,YEA,NPLU,CHI2
5027 DIMENSION GRAD(*),XVAL(*),AMS(500),DISTAZ(500)
5029 DIMENSION XP(NPLANE),YP(NPLANE),
5030 & COV(NPLANE,NPLANE),AP(NPLANE),COVY(NPLANE,NPLANE)
5031 DIMENSION VECT(7),VOUT(7)
5037 IF (XVAL(1).LT.0.) ASIGN = -1.
5040 PXZ = DABS(1./XVAL(1))
5045 PTOT = PXZ/DCOS(ALAM)
5046 TTHET = DSQRT(DTAN(ALAM)**2+DSIN(PHI)**2)/DCOS(PHI)
5048 PT = DSQRT(PX**2+PY**2)
5050 RL3 = ASIGN*PT / (0.299792458D-3 * BL3)
5051 ALPHA = DATAN2(PY,PX)
5052 XC = XV+RL3*DSIN(ALPHA)
5053 YC = YV-RL3*DCOS(ALPHA)
5067 ANGDEV = (ZPLANEP(1)-ZEA)*TTHET/RL3
5068 XP(1) = XC+RL3*DSIN(ANGDEV-ALPHA)
5069 YP(1) = YC+RL3*DCOS(ANGDEV-ALPHA)
5070 AL = THICK/ DCOS(THET)
5071 AP(1) = (0.0136D0/PTOT) * DSQRT(AL) * (1 + 0.038D0*DLOG(AL))
5073 ANGDEV = (ZPLANEP(2)-ZEA)*TTHET/RL3
5074 XP(2) = XC+RL3*DSIN(ANGDEV-ALPHA)
5075 YP(2) = YC+RL3*DCOS(ANGDEV-ALPHA)
5076 AL = THICK/ DCOS(THET)
5077 AP(2) = (0.0136D0/PTOT) * DSQRT(AL) * (1 + 0.038D0*DLOG(AL))
5079 ANGDEV = (ZPLANEP(3)-ZEA)*TTHET/RL3
5080 XP(3) = XC+RL3*DSIN(ANGDEV-ALPHA)
5081 YP(3) = YC+RL3*DCOS(ANGDEV-ALPHA)
5082 AL = THICK/ DCOS(THET)
5083 AP(3) = (0.0136D0/PTOT) * DSQRT(AL) * (1 + 0.038D0*DLOG(AL))
5085 ANGDEV = (ZPLANEP(4)-ZEA)*TTHET/RL3
5086 XP(4) = XC+RL3*DSIN(ANGDEV-ALPHA)
5087 YP(4) = YC+RL3*DCOS(ANGDEV-ALPHA)
5088 AL = THICK/ DCOS(THET)
5089 AP(4) = (0.0136D0/PTOT) * DSQRT(AL) * (1 + 0.038D0*DLOG(AL))
5091 ANGDEV = (700.D0-ZEA)*TTHET/RL3
5092 XPL3 = XC+RL3*DSIN(ANGDEV-ALPHA)
5093 YPL3 = YC+RL3*DCOS(ANGDEV-ALPHA)
5094 PX = PT*DCOS(ANGDEV-ALPHA)
5095 PY = -PT*DSIN(ANGDEV-ALPHA)
5096 PHIC = DATAN2(PY,PX)
5097 CX = DSIN(THET)*DCOS(PHIC)
5098 CY = DSIN(THET)*DSIN(PHIC)
5101 VECT(1) = XPL3+(ZMAGS-700.)*CX/CZ
5102 VECT(2) = YPL3+(ZMAGS-700.)*CY/CZ
5109 PXZ = PTOT*DSQRT(VECT(4)**2+VECT(6)**2)
5110 RDIP = ASIGN*PXZ / (0.299792458D-3 * B0)
5111 PHI1 = DATAN2(VECT(4),VECT(6))
5112 XC = VECT(1)+RDIP*DCOS(PHI1)
5113 ZC = VECT(3)-RDIP*DSIN(PHI1)
5115 IF (DABS((ZMAGE-ZPLANEP(5))/RDIP).GE.1.D0) RETURN ! Particule boucle dans le champ
5116 XP(5) = XC-ASIGN*DSQRT(RDIP**2-(ZPLANEP(5)-ZC)**2)
5117 YP(5) = VECT(2)+(ZPLANEP(5)-ZMAGS)*VECT(5)/VECT(6)
5118 CX = (ZPLANEP(5)-ZC)/RDIP
5120 CZ = DABS((XP(5)-XC)/RDIP)
5121 THET = DATAN2(DSQRT(CX**2+CY**2),CZ)
5122 AL = THICK/ DCOS(THET)
5123 AP(5) = (0.0136D0/PTOT) * DSQRT(AL) * (1 + 0.038D0*DLOG(AL))
5125 IF (DABS((ZMAGE-ZPLANEP(6))/RDIP).GE.1.D0) RETURN ! Particule boucle dans le champ
5126 XP(6) = XC-ASIGN*DSQRT(RDIP**2-(ZPLANEP(6)-ZC)**2)
5127 YP(6) = VECT(2)+(ZPLANEP(6)-ZMAGS)*VECT(5)/VECT(6)
5128 CX = (ZPLANEP(6)-ZC)/RDIP
5130 CZ = DABS((XP(6)-XC)/RDIP)
5131 THET = DATAN2(DSQRT(CX**2+CY**2),CZ)
5132 AL = THICK/ DCOS(THET)
5133 AP(6) = (0.0136D0/PTOT) * DSQRT(AL) * (1 + 0.038D0*DLOG(AL))
5136 IF (DABS((ZMAGE-ZC)/RDIP).GE.1.D0) RETURN ! Particule boucle dans le champ
5137 VOUT(1) = XC-ASIGN*DSQRT(RDIP**2-(ZMAGE-ZC)**2)
5138 VOUT(2) = VECT(2)+(ZMAGE-ZMAGS)*VECT(5)/VECT(6)
5140 VOUT(4) = (ZMAGE-ZC)/RDIP
5142 VOUT(6) = DABS((VOUT(1)-XC)/RDIP)
5150 THET = DATAN2(DSQRT(VECT(4)**2+VECT(5)**2),VECT(6))
5151 XP(7) = VECT(1)+(ZPLANEP(7)-ZMAGE)*VECT(4)/VECT(6)
5152 YP(7) = VECT(2)+(ZPLANEP(7)-ZMAGE)*VECT(5)/VECT(6)
5153 AL = THICK/ DCOS(THET)
5154 AP(7) = (0.0136D0/PTOT) * DSQRT(AL) * (1 + 0.038D0*DLOG(AL))
5156 XP(8) = XP(7)+(ZPLANEP(8)-ZPLANEP(7))*VECT(4)/VECT(6)
5157 YP(8) = YP(7)+(ZPLANEP(8)-ZPLANEP(7))*VECT(5)/VECT(6)
5158 AL = THICK/ DCOS(THET)
5159 AP(8) = (0.0136D0/PTOT) * DSQRT(AL) * (1 + 0.038D0*DLOG(AL))
5162 XP(9) = XP(8)+(ZPLANEP(9)-ZPLANEP(8))*VECT(4)/VECT(6)
5163 YP(9) = YP(8)+(ZPLANEP(9)-ZPLANEP(8))*VECT(5)/VECT(6)
5164 AL = THICK/ DCOS(THET)
5165 AP(9) = (0.0136D0/PTOT) * DSQRT(AL) * (1 + 0.038D0*DLOG(AL))
5168 XP(10) = XP(9)+(ZPLANEP(10)-ZPLANEP(9))*VECT(4)/VECT(6)
5169 YP(10) = YP(9)+(ZPLANEP(10)-ZPLANEP(9))*VECT(5)/VECT(6)
5170 AL = THICK/ DCOS(THET)
5171 AP(10) = (0.0136D0/PTOT) * DSQRT(AL) * (1 + 0.038D0*DLOG(AL))
5173 ** Matrice de covariance
5176 IF (LPLANE(II).EQ.1) THEN
5181 IF (LPLANE(JJ).EQ.1) THEN
5187 COV(J,I) =COV(J,I) + XPREC**2
5190 * IF (I .EQ. 10 .AND. J .EQ. 10) PRINT *,'10 10 ',COV(J,I)
5193 & + (ZPLANEP(II) + DISTAZ(L))*(ZPLANEP(JJ) + DISTAZ(L))*AMS(L)**2
5197 & + (ZPLANEP(II)-ZPLANEP(K))*(ZPLANEP(JJ)-ZPLANEP(K))*AP(K)**2
5198 * IF (I .EQ. 10 .AND. J .EQ. 10) PRINT *,'10 10 ',COV(J,I)
5201 COVY(J,I) = COV(J,I)
5203 COVY(J,I) = COVY(J,I) - XPREC**2 + YPREC**2
5210 * Inversion des matrices de covariance
5214 CALL DSINV(NPLU, COV, NPLANE, IFAIL)
5215 ** IF (JFAIL.NE.0 .AND. IFAIL .NE. 0) STOP 'ERROR'
5216 IF (IFAIL .NE. 0) STOP 'ERROR'
5218 CALL DSINV(NPLU, COVY, NPLANE, IFAIL)
5219 ** IF (JFAIL.NE.0 .AND. IFAIL .NE. 0) STOP 'ERROR'
5220 IF (IFAIL .NE. 0) STOP 'ERROR'
5221 * PRINT *,' COVARIANCE MATRIX AFTER'
5223 * PRINT *,(COV(J,I),J=1,NPLANE)
5226 ** Calcul de FVAL ou CHI2
5230 IF (LPLANE(II).EQ.1) THEN
5235 IF (LPLANE(JJ).EQ.1) THEN
5238 FVAL = FVAL + COV(J,I)*(XMP(II)-XP(II))*(XMP(JJ)-XP(JJ))
5239 FVAL = FVAL + COVY(J,I)*(YMP(II)-YP(II))
5241 ** IF (JJ.EQ.II) THEN
5242 ** FVAL = FVAL + (XMP(II)-XP(II))*(XMP(JJ)-XP(JJ))/XPREC**2
5243 ** FVAL = FVAL + (YMP(II)-YP(II))
5244 ** & *(YMP(JJ)-YP(JJ))/YPREC**2
5251 print *,' fcnfit pxz tphi talam xvert yvert chi2',
5252 & PXZEA,PHIEA,ALAMEA,
5253 & XEA,YEA,CHI2/FLOAT(2*NPLU-5)
5255 1000 FORMAT(I5,7F12.6)