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)
88 common/dstation/idstation
96 CALL RECO_READEVT(NEV,IDRES,IREADGEANT)
100 CALL RECO_TRACKF(IDRES,IREADGEANT)
112 ********************************************
114 ********************************************
115 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
116 PARAMETER(NBSTATION=5)
118 COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
120 * ZPLANE = position de la premiere chambre de chaque station
121 * exemple (pour la premiere station) :
122 * zch (dans AliMUONv0) = 540 (position du centre des 2 chambres)
124 * ZPLANE = -530 si DZ_PL = 20cm
127 * DATA DZ_PL/8.,8.,8.,8.,8./
128 * DATA ZPLANE/-511.,-686.0,-971.,-1245.,-1445./
131 DATA DZ_PL/20.,20.,20.,20.,20./
132 DATA ZPLANE/-518.,-680.,-965.,-1239.,-1439./
135 * DATA ZCOIL,ZMAGEND/-825.0,-1125./ ! Constant field 3 Tm
136 DATA ZCOIL,ZMAGEND/-805.0,-1233./ ! CCC magn. field map M.B
137 ** DATA ZCOIL,ZMAGEND/-795.1,-1242.9/ ! CCC magn. field map M.B
142 ********************************************
143 SUBROUTINE cutpxz(spxzcut)
144 ********************************************
145 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
147 COMMON/ACUTPXZ/ACUTPXZ
153 ********************************************
154 SUBROUTINE sigmacut(ssigcut)
155 ********************************************
156 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
158 COMMON/TRACKFI/EFF,EFF1,EFF2,XPREC,YPREC,PHIPREC,ALAMPREC,
159 & HCUT,LBKG,SIGCUT,ALPHATOP,HTOP
165 ********************************************
166 SUBROUTINE xpreci(sxprec)
167 ********************************************
168 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
170 COMMON/TRACKFI/EFF,EFF1,EFF2,XPREC,YPREC,PHIPREC,ALAMPREC,
171 & HCUT,LBKG,SIGCUT,ALPHATOP,HTOP
177 ********************************************
178 SUBROUTINE ypreci(syprec)
179 ********************************************
180 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
182 COMMON/TRACKFI/EFF,EFF1,EFF2,XPREC,YPREC,PHIPREC,ALAMPREC,
183 & HCUT,LBKG,SIGCUT,ALPHATOP,HTOP
189 *************************************************************************
190 SUBROUTINE RECO_INIT(seff,sb0,sbl3)
191 *************************************************************************
192 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
194 CALL TRACKF_INIT(seff,sb0,sbl3)
201 *************************************************************************
202 SUBROUTINE TRACKF_INIT(seff,sb0,sbl3)
203 *************************************************************************
204 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
206 PARAMETER(NBSTATION=5,NTRMAX=500)
208 COMMON/REVENT/IEVBKGI,NBKGMAX,MAXUPSEV
212 COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
214 COMMON/FILED/FILERES,FILEBKG,FILEOUT,FILEMIN
216 COMMON/TRACKFI/EFF,EFF1,EFF2,XPREC,YPREC,PHIPREC,ALAMPREC,
217 & HCUT,LBKG,SIGCUT,ALPHATOP,HTOP
219 COMMON/TRACKSUM/NRES(5),NRESF,NTRMUALL,NMUONALL,NGHOSTALL,
220 & NTRACKFALL,NERRALL(NBSTATION),IR
222 COMMON/ACUTPXZ/ACUTPXZ
233 AMAGLEN = ZMAGEND-ZCOIL
234 ZM = AMAGLEN/2.+ZCOIL
235 ALPHATOP = 0.01*0.3*B0*ABS(AMAGLEN)
237 HCUT = ABS(HTOP)/PXZCUT
239 print*,'TRACK_INIT hcut= ',hcut
240 print*,'TRACK_INIT eff = ',eff
241 print*,'TRACK_INIT b0 = ',b0
242 print*,'TRACK_INIT bl3 = ',bl3
243 print*,'TRACK_INIT sigmacut = ',sigcut
244 print*,'TRACK_INIT cutpxz = ',pxzcut
245 print*,'TRACK_INIT xprec = ',xprec
246 print*,'TRACK_INIT yprec = ',yprec
248 EFF2 = EFF**2 ! PROBA. DEUX CHAMBRES TOUCHES
249 EFF1 = EFF2+2.*EFF*(1.-EFF) ! PROBA. AU MOINS UNE CHAMBRE TOUCHE
250 ** Used only for stations 4 & 5
251 PHIPREC = SQRT(2.)*XPREC/DZ_PL(5) ! PHI = (OZ , PROJ. DANS XOZ)
252 ALAMPREC = SQRT(2.)*YPREC/DZ_PL(5) ! LAM = (OM , PROJ. DANS XOZ)
270 *************************************************************************
272 *************************************************************************
275 *************************************************************************
277 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
279 PARAMETER(NPLANE=10,NBSTATION=5)
281 COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
283 COMMON/PARAM/ZPLANEP(NPLANE),THICK,XPREC,YPREC,B0,BL3,ZMAGS,
284 & ZMAGE,ZABS,XMAG,ZBP1,ZBP2,CONST
286 COMMON/PRECCUT/PCUT,PTCUT,CHI2CUT
288 COMMON/PRECSUM/NRESF1,NMUONALL1,NGHOSTALL1,NTRACKFALL1
293 DATA THICK/0.03D0/ ! X/X0=3%
294 ** DATA THICK/0.02D0/ ! X/X0=2% chambre
296 ** DATA B0,BL3/10.,2.0/ ! Champ magnetique dans le dipole et dans L3 en kgauss
297 DATA B0,BL3/7.,2.0/ ! Magnetic field in the dipole & L3 in kgauss
298 DATA ZMAGS/805.0D0/,ZMAGE/1233.0D0/,ZABS/503.D0/ ! CCC not used when
300 * DATA ZMAGS/825.0D0/,ZMAGE/1125.0D0/,ZABS/503.D0/
303 ** DATA XPREC/0.0100D0/,YPREC/0.144337D0/ ! CCC
304 DATA XPREC/0.0100D0/,YPREC/0.2D0/
307 CONST = 0.299792458D-3*B0*(ZMAGE-ZMAGS)
311 ZPLANEP(J) = -ZPLANE(I)
313 ZPLANEP(J) = -ZPLANE(I)+DZ_PL(I)
317 print*,'zplanep(',i,')=',ZPLANEP(I)
320 PCUT = 3. ! Coupure en PXZ muon (GeV/c)
321 PTCUT = 0. ! Coupure en Pt muon (GeV/c)
323 CHI2CUT = 1.E4 ! Coupure sur le CHI2 du fit
326 X02 = 10.397 ! Concrete (cm)
327 X03 = 0.56 ! Plomb (cm)
328 X04 = 47.26 ! Polyethylene (cm)
331 ** Calcul des parametres pour la correction de Branson de l'absorbeur
332 ANBP = (315.**3-90.**3)/X01 +(467.**3-315.**3)/X02+
333 & (472.**3-467.**3)/X03+(477.**3-472.**3)/X04+
334 & (482.**3-477.**3)/X03+(487.**3-482.**3)/X04+
335 & (492.**3-487.**3)/X03+(497.**3-492.**3)/X04+
336 & (502.**3-497.**3)/X03
337 ADBP = (315.**2-90.**2)/X01 +(467.**2-315.**2)/X02+
338 & (472.**2-467.**2)/X03+(477.**2-472.**2)/X04+
339 & (482.**2-477.**2)/X03+(487.**2-482.**2)/X04+
340 & (492.**2-487.**2)/X03+(497.**2-492.**2)/X04+
341 & (502.**2-497.**2)/X03
342 ZBP1 = 2./3.*ANBP/ADBP
343 ANBP = (315.**3-90.**3)/X01 +(467.**3-315.**3)/X02+
344 & (503.**3-467.**3)/X05
345 ADBP = (315.**2-90.**2)/X01 +(467.**2-315.**2)/X02+
346 & (503.**2-467.**2)/X05
347 ZBP2 = 2./3.*ANBP/ADBP
349 IF (IDEBUG.GE.1) THEN
350 PRINT *,' PREC_INIT B0 (kgauss)',B0,' BL3 (kgauss)',BL3
351 PRINT *,' PREC_INIT ZMAGE (cm)',ZMAGE,' ZMAGS (cm)',ZMAGS,
353 PRINT *,' PREC_INIT ZABS (cm)',ZABS,' ZBP1 (cm)',ZBP1,
355 PRINT *,' PREC_INIT Radiation length absorber X01 (cm)',X01,
357 PRINT *,' PREC_INIT X03(cm)',X03,' X04 (cm)',X04,' X05 (cm)',
359 PRINT *,' PREC_INIT Radiation length chamber THICK (%)',
361 PRINT *,' PREC_INIT XPREC (cm)',XPREC,' YPREC (cm)',YPREC
362 PRINT *,' PREC_INIT Coupure en Pxz (GeV/c): ',PCUT
363 PRINT *,' PREC_INIT Coupure en Pt (GeV/c): ',PTCUT
364 PRINT *,' PREC_INIT Coupure en CHI2 : ',CHI2CUT
374 * Magnetic Field Map GC
375 ** OPEN (UNIT=40,FILE='data/field02.dat',
376 ** & STATUS='UNKNOWN')
385 *************************************************************************
386 SUBROUTINE RECO_READEVT(NEV,IDRES,IREADGEANT)
387 *************************************************************************
390 *************************************************************************
391 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
393 PARAMETER(NTRMAX=500)
395 PARAMETER (NBSTATION=5,MAXIDG=20000,MAXHITTOT=20000,
396 & MAXHITCH=10000,MAXHIT=1000,NBCHAMBER=10)
398 COMMON/RHITG/ITYPG(MAXIDG),XTRG(MAXIDG),YTRG(MAXIDG),
399 & PTOTG(MAXIDG),IDG(MAXIDG),IZCH(MAXIDG),
400 & PVERT1G(MAXIDG),PVERT2G(MAXIDG),PVERT3G(MAXIDG),
401 & ZVERTG(MAXIDG),NHITTOT1,CX(MAXIDG),CY(MAXIDG),CZ(MAXIDG),
402 & XGEANT(MAXIDG),YGEANT(MAXIDG),CLSIZE1(MAXIDG),CLSIZE2(MAXIDG)
404 DIMENSION TYPG(MAXIDG),ZCH(MAXIDG)
409 IF (IREADGEANT.eq.1) THEN ! GEANT hits
411 CALL TRACKF_READ_GEANT(ITYPG,XTRG,YTRG,PTOTG,IDG,IZCH,PVERT1G,
412 & PVERT2G,PVERT3G,ZVERTG,NHITTOT1,CX,CY,CZ,IEVR,NEV,
413 & XGEANT,YGEANT,CLSIZE1,CLSIZE2)
414 ELSE ! reconstructed hits
415 CALL TRACKF_READ_SPOINT(ITYPG,XTRG,YTRG,PTOTG,IDG,IZCH,PVERT1G,
416 & PVERT2G,PVERT3G,ZVERTG,NHITTOT1,CX,CY,CZ,IEVR,NEV,
417 & XGEANT,YGEANT,CLSIZE1,CLSIZE2)
422 call chfill(100,sngl(typg(i)),R1,R2)
423 call chfill(101,sngl(ygeant(i)),R1,R2)
424 call chfill(102,sngl(xgeant(i)),R1,R2)
426 call chfill(103,sngl(zch(i)),R1,R2)
427 call chfill(104,sngl(ptotg(i)),R1,R2)
428 call chfill(105,sngl(pvert2g(i)),R1,R2)
429 call chfill(106,sngl(pvert1g(i)),R1,R2)
430 call chfill(107,sngl(pvert3g(i)),R1,R2)
431 call chfill(108,sngl(zvertg(i)),R1,R2)
432 call chfill(109,sngl(ytrg(i)),R1,R2)
433 call chfill(110,sngl(xtrg(i)),R1,R2)
437 CALL CHFILL (999,SNGL(PTOTG(I)),R1,R2)
440 CALL TRACKF_STAT(IDRES,IREADGEANT)
445 *************************************************************************
446 SUBROUTINE TRACKF_STAT(IDRES,IREADGEANT)
447 *************************************************************************
448 * Associate hits between two chambers inside a station
449 * Simulate spatial resolution and chamber efficiency
451 *************************************************************************
452 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
454 PARAMETER (NBSTATION=5,MAXIDG=20000,MAXHITTOT=20000,
455 & MAXHITCH=10000,MAXHIT=1000,NBCHAMBER=10)
457 COMMON/TRACKFI/EFF,EFF1,EFF2,XPREC,YPREC,PHIPREC,ALAMPREC,
458 & HCUT,LBKG,SIGCUT,ALPHATOP,HTOP
460 COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
462 * HITS GEANT initiaux par chambre
463 COMMON/RHITG/ITYPG(MAXIDG),XTRG(MAXIDG),YTRG(MAXIDG),
464 & PTOTG(MAXIDG),IDG(MAXIDG),IZCH(MAXIDG),
465 & PVERT1G(MAXIDG),PVERT2G(MAXIDG),PVERT3G(MAXIDG),
466 & ZVERTG(MAXIDG),NHITTOT1,CX(MAXIDG),CY(MAXIDG),CZ(MAXIDG),
467 & XGEANT(MAXIDG),YGEANT(MAXIDG),CLSIZE1(MAXIDG),CLSIZE2(MAXIDG)
469 * HITS GEANT associes par station
470 COMMON/RHIT/ITYP(MAXHITTOT),XTR(MAXHITTOT),YTR(MAXHITTOT),
471 & PTOT(MAXHITTOT),ID(MAXHITTOT),IZST(MAXHITTOT),
472 & PVERT1(MAXHITTOT),PVERT2(MAXHITTOT),PVERT3(MAXHITTOT),
473 & ZVERT(MAXHITTOT),NHITTOT
476 COMMON/CHHIT/XM(NBSTATION,MAXHITCH),YM(NBSTATION,MAXHITCH),
477 & PHM(NBSTATION,MAXHITCH),ALM(NBSTATION,MAXHITCH),
478 & IZM(NBSTATION,MAXHITCH),
479 & IP(NBSTATION,MAXHITCH),JHIT(NBSTATION),
480 & XMR(NBSTATION,MAXHITCH,2),YMR(NBSTATION,MAXHITCH,2)
483 common/dstation/idstation
485 DIMENSION RMIN(NBCHAMBER),RMAX1(NBCHAMBER)
486 DIMENSION XMA(NBCHAMBER,MAXHITCH),YMA(NBCHAMBER,MAXHITCH),
487 & IMARK(NBCHAMBER,MAXHITCH)
489 DIMENSION IEFFI(MAXHITTOT)
490 DIMENSION IH(NBCHAMBER,MAXHIT)
491 DIMENSION NHIT(NBCHAMBER)
492 DIMENSION DXMAX(NBSTATION),DYMAX(NBSTATION),VRES(2,5)
494 REAL*4 RNDM,RN,RN1,RN2,R1,R2
497 DATA RMAX1/91.5,91.5,122.5,122.5,158.3,158.3,260.,260.,260.,260./
498 * Zone de recherche entre deux plans d'une station
500 * if (idstation.eq.8) then
501 * DATA DXMAX/1.,1.,1.2,2.4,2.4/ ! dz_ch = 8 cm
502 * else if (idstation.eq.20) then
503 DATA DXMAX/1.5,1.5,3.,6.,6./ ! dz_ch = 20cm
506 DATA DYMAX/0.22,0.22,0.22,0.22,0.22/ ! CCC Upsilon dz_ch = 20 cm
514 RMIN(ICH) = ABS(ZPLANE(IZ)*TAN(2.*ACOS(-1.)/180))
515 IF (IZ.GT.2) RMIN(ICH) = 30.
517 RMIN(ICH) = ABS(ZPLANE(IZ)*TAN(2.*ACOS(-1.)/180))
518 IF (IZ.GT.2) RMIN(ICH) = 30.
526 * 1 ere boucle de lecture des hits initiaux
529 IF (IREADGEANT.EQ.1) THEN
538 IZ = INT(FLOAT(ICH+1)/2.)
541 IF (IMOD.NE.0.AND.IZ.LE.5) THEN
542 CALL CHFILL2(1000+IZ,SNGL(XGEANT(I)),SNGL(YGEANT(I)),R2)
546 ISTAK = MOD(ISTAK,30000)
547 ISTAK = MOD(ISTAK,10000)
549 IF ((ITYPG(I).EQ.5.OR.ITYPG(I).EQ.6).AND.
550 & ISTAK.EQ.IDRES) THEN ! upsilon
552 VRES(1,IMU) = XGEANT(I)
553 VRES(2,IMU) = YGEANT(I)
560 DO I = 1,NHITTOT1 ! Boucle sur les hits GEANT de toutes les ch.
562 ** IF (ITYPG(I).NE.5.AND.ITYPG(I).NE.6) GOTO 1 ! CCC
566 IF (ICH.GT.10) GO TO 1
568 IF (IREADGEANT.EQ.1) THEN ! GEANT hits
570 IF (ICH.EQ.9.OR.ICH.EQ.10) THEN
573 DNU = SQRT((XGEANT(I)-VRES(1,IM))**2+
574 & (YGEANT(I)-VRES(2,IM))**2)
575 IF (DNU.LT.DNUM) DNUM = DNU
577 IF (DNUM.GT.50.) GO TO 1 ! discard hits far from MUONS
580 CALL RANNOR(RN1,RN2) ! CCC
583 * X = XGEANT(I) + RN1 * XPREC
584 * Y = YGEANT(I) + RN2 * YPREC
585 * efficacite des chambres
588 IF (RN.GT.EFF) IEFFI(I) = 0
589 ** IF (ICH.EQ.9.OR.ICH.EQ.10) THEN
590 ** PRINT *,' HIT GEANT',' ICH=',ICH,' I=',I,' X =',X,' Y=',
591 ** & Y,' IDG=',IDG(I)
594 IF (ITYPG(I).EQ.5.OR.ITYPG(I).EQ.6) THEN
596 ISTAK = MOD(ISTAK,30000)
597 ISTAK = MOD(ISTAK,10000)
599 IF (ISTAK.EQ.IDRES) then
602 IZ = INT(FLOAT(IZCH(i)+1)/2.)
605 call chfill(ichx,sngl(dy),R1,R2)
606 call chfill(ichy,sngl(dx),R1,R2)
607 if (iz.eq.1) call chfill(116,sngl(dy),R1,R2)
608 if (iz.eq.2) call chfill(117,sngl(dy),R1,R2)
612 ELSE ! reconstructed hits
619 * étude des hits geant avec un seul fichier
620 * CALL RANNOR(RN1,RN2) ! CCCC
621 * X = XGEANT(I) + RN1 * XPREC
622 * Y = YGEANT(I) + RN2 * YPREC
624 IF (ITYPG(I).EQ.5.OR.ITYPG(I).EQ.6) THEN
626 ISTAK = MOD(ISTAK,30000)
627 ISTAK = MOD(ISTAK,10000)
629 IF (ISTAK.EQ.IDRES) then
632 IZ = INT(FLOAT(IZCH(i)+1)/2.)
635 call chfill(ichx,sngl(dy),R1,R2)
636 call chfill(ichy,sngl(dx),R1,R2)
637 if (iz.eq.1) call chfill(116,sngl(dy),R1,R2)
638 if (iz.eq.2) call chfill(117,sngl(dy),R1,R2)
645 ** IF (R.LT.RMIN(ICH).OR.R.GT.RMAX1(ICH)) then
646 ** if (ich.le.10) then
647 ** print*,'* chambre ',ich,' * hit ',i
648 ** print*,'ityp=',itypg(i),' x=',X,' y=',Y
649 ** print*,'R=',R,' RMIN=',RMIN(ICH),' RMAX1=',RMAX1(ICH)
654 NHIT(ICH) = NHIT(ICH)+1
655 IH(ICH,NHIT(ICH)) = I
656 XMA(ICH,NHIT(ICH)) = X
657 YMA(ICH,NHIT(ICH)) = Y
658 IMARK(ICH,NHIT(ICH)) = 0
663 * Association des hits entre chambres d'une station
664 II = 0 ! nombre de hits GEANT par station
665 DO ICH1 = 1,10,2 ! loop on chamber
666 IZ = INT(FLOAT(ICH1+1)/2.)
670 DO I1 = 1,NHIT(ICH1) ! loop on hits in 1st chamber
681 PVERT1(II) = PVERT1G(I)
682 PVERT2(II) = PVERT2G(I)
683 PVERT3(II) = PVERT3G(I)
684 ZVERT(II) = ZVERTG(I)
686 IF (IEFFI(I).EQ.1) THEN
690 XEXT1 = (ZPLANE(IZ)-DZ_PL(IZ))/ZPLANE(IZ)*X1
691 YEXT1 = (ZPLANE(IZ)-DZ_PL(IZ))/ZPLANE(IZ)*Y1
693 DO I2 = 1,NHIT(ICH2) ! loop on hits in 2nd chamber
696 IF (IEFFI(J).EQ.1) THEN
704 & (ITYP(II).EQ.5.OR.ITYP(II).EQ.6)) THEN
705 CALL CHFILL(70+IZ,SNGL(DX),R1,R2)
706 CALL CHFILL(80+IZ,SNGL(DY),R1,R2)
711 IF (DX.LT.DXMAX(IZ).AND.DY.LT.(SIGCUT*DYMAX(IZ)) ! CCC
715 JHIT(IZ) = JHIT(IZ)+1
719 PHM(IZ,JHIT(IZ)) = -ATAN((X2-X1)/DZ_PL(IZ))
720 ALM(IZ,JHIT(IZ)) = ATAN((Y2-Y1)/DZ_PL(IZ)*
721 & COS(PHM(IZ,JHIT(IZ))))
722 XMR(IZ,JHIT(IZ),1) = X1
723 YMR(IZ,JHIT(IZ),1) = Y1
724 XMR(IZ,JHIT(IZ),2) = X2
725 YMR(IZ,JHIT(IZ),2) = Y2
729 ISTAK = MOD(ISTAK,30000)
730 ISTAK = MOD(ISTAK,10000)
732 IF ((ITYPG(J).EQ.5.OR.ITYPG(J).EQ.6).AND.
733 & ISTAK.EQ.IDRES) THEN ! upsilon or J/psi
740 PVERT1(II) = PVERT1G(J)
741 PVERT2(II) = PVERT2G(J)
742 PVERT3(II) = PVERT3G(J)
743 ZVERT(II) = ZVERTG(J)
749 ENDDO ! loop on hits in 2nd chamber
751 IF (IFIND.EQ.0) THEN ! No possible association
752 JHIT(IZ) = JHIT(IZ)+1
757 PHM(IZ,JHIT(IZ)) = 10.
758 ALM(IZ,JHIT(IZ)) = 10.
759 XMR(IZ,JHIT(IZ),1) = X1
760 YMR(IZ,JHIT(IZ),1) = Y1
761 XMR(IZ,JHIT(IZ),2) = 0.
762 YMR(IZ,JHIT(IZ),2) = 0.
765 ENDDO ! end loop on hits in 1st chamber
766 ENDDO ! end loop on chamber
769 * On conserve les HITS (x,y) de la 2nde chambre des stations
771 DO ICH = 2,10,2 ! Loop on 2nd chambers
772 IZ = INT(FLOAT(ICH+1)/2.)
773 DO I = 1,NHIT(ICH) ! Loop on hits
776 IF (IMARK(ICH,I).EQ.0) THEN ! hit not already associated
786 PVERT1(II) = PVERT1G(J)
787 PVERT2(II) = PVERT2G(J)
788 PVERT3(II) = PVERT3G(J)
789 ZVERT(II) = ZVERTG(I)
791 IF (IEFFI(J).EQ.1) THEN
792 JHIT(IZ) = JHIT(IZ)+1
793 XM(IZ,JHIT(IZ)) = XMA(ICH,I)
794 YM(IZ,JHIT(IZ)) = YMA(ICH,I)
796 PHM(IZ,JHIT(IZ)) = 10.
797 ALM(IZ,JHIT(IZ)) = 10.
799 XMR(IZ,JHIT(IZ),1) = 1000.
800 YMR(IZ,JHIT(IZ),1) = 1000.
801 XMR(IZ,JHIT(IZ),2) = XMA(ICH,I)
802 YMR(IZ,JHIT(IZ),2) = YMA(ICH,I)
805 ENDDO ! End loop on hits
806 ENDDO ! End loop on 2nd chambers
809 NHITTOT = II ! total number of hits in stations
812 IF (IDEBUG.GE.2) THEN
813 PRINT *,'TRACKF_STAT nb hits:',NHITTOT
817 ** PRINT *,' IZ=',IZ,' JHIT(IZ)=',JHIT(IZ)
820 ** PRINT *,' HIT ASS.',' IZ=',IZ,' II=',II,' X =',
821 ** & XM(IZ,J),' Y=',YM(IZ,J),' ID=',ID(II)
822 ** PRINT *,' HIT ASS.',' IZ=',IZ,' II=',II,' X1 =',
823 ** & XMR(IZ,J,1),' Y1=',YMR(IZ,J,1),' ID=',ID(II)
824 ** PRINT *,' HIT ASS.',' IZ=',IZ,' II=',II,' X2 =',
825 ** & XMR(IZ,J,2),' Y2=',YMR(IZ,J,2),' ID=',ID(II)
829 ** DO IZ = 1,NBSTATION
830 ** PRINT *,' IZ=',IZ,' JHIT(IZ)=',JHIT(IZ)
832 ** PRINT *,' PHM(IZ,J)=',PHM(IZ,J),' ALM(IZ,J)=',ALM(IZ,J)
838 *************************************************************************
839 SUBROUTINE TRACKF_STAT_NEW(IDRES)
840 *************************************************************************
841 * Associate hits between two chambers inside a station
842 * Simulate spatial resolution and chamber efficiency
844 *************************************************************************
845 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
847 PARAMETER (NBSTATION=5,MAXIDG=20000,MAXHITTOT=20000,
848 & MAXHITCH=10000,MAXHIT=1000,NBCHAMBER=10)
850 COMMON/TRACKFI/EFF,EFF1,EFF2,XPREC,YPREC,PHIPREC,ALAMPREC,
851 & HCUT,LBKG,SIGCUT,ALPHATOP,HTOP
853 COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
855 * HITS GEANT initiaux par chambre
856 COMMON/RHITG/ITYPG(MAXIDG),XTRG(MAXIDG),YTRG(MAXIDG),
857 & PTOTG(MAXIDG),IDG(MAXIDG),IZCH(MAXIDG),
858 & PVERT1G(MAXIDG),PVERT2G(MAXIDG),PVERT3G(MAXIDG),
859 & ZVERTG(MAXIDG),NHITTOT1,CX(MAXIDG),CY(MAXIDG),CZ(MAXIDG),
860 & XGEANT(MAXIDG),YGEANT(MAXIDG),CLSIZE1(MAXIDG),CLSIZE2(MAXIDG)
862 * HITS GEANT associes par station
863 COMMON/RHIT/ITYP(MAXHITTOT),XTR(MAXHITTOT),YTR(MAXHITTOT),
864 & PTOT(MAXHITTOT),ID(MAXHITTOT),IZST(MAXHITTOT),
865 & PVERT1(MAXHITTOT),PVERT2(MAXHITTOT),PVERT3(MAXHITTOT),
866 & ZVERT(MAXHITTOT),NHITTOT
869 COMMON/CHHIT/XM(NBSTATION,MAXHITCH),YM(NBSTATION,MAXHITCH),
870 & PHM(NBSTATION,MAXHITCH),ALM(NBSTATION,MAXHITCH),
871 & IZM(NBSTATION,MAXHITCH),
872 & IP(NBSTATION,MAXHITCH),JHIT(NBSTATION),
873 & XMR(NBSTATION,MAXHITCH,2),YMR(NBSTATION,MAXHITCH,2)
877 DIMENSION RMIN(NBCHAMBER),RMAX1(NBCHAMBER)
878 DIMENSION XMA(NBCHAMBER,MAXHITCH),YMA(NBCHAMBER,MAXHITCH),
879 & IMARK(NBCHAMBER,MAXHITCH)
881 DIMENSION IEFFI(MAXHITTOT)
882 DIMENSION IH(NBCHAMBER,MAXHIT)
883 DIMENSION NHIT(NBCHAMBER)
884 DIMENSION DXMAX(NBSTATION),DYMAX(NBSTATION),I2C(1000)
886 DIMENSION DIST(2),NMUON(2),NHITMUON(2,5),NMUONGOOD(2)
888 REAL*4 RNDM,RN,RN1,RN2,R1,R2
891 DATA RMAX1/91.5,91.5,122.5,122.5,158.3,158.3,260.,260.,260.,260./
892 * Zone de recherche entre deux plans d'une station
893 ** DATA DXMAX/2.,1.5,2.5,3.,3./
894 DATA DXMAX/1.5,1.5,3.,3.,3./
895 DATA DYMAX/0.22,0.22,0.21,0.21,0.21/
901 RMIN(ICH) = ABS(ZPLANE(IZ)*TAN(2.*ACOS(-1.)/180))
902 IF (IZ.GT.2) RMIN(ICH) = 30.
904 RMIN(ICH) = ABS(ZPLANE(IZ)*TAN(2.*ACOS(-1.)/180))
905 IF (IZ.GT.2) RMIN(ICH) = 30.
919 IF (IZCH(I).EQ.NCH) THEN
921 ISTAK = MOD(ISTAK,30000)
922 ISTAK = MOD(ISTAK,10000)
923 IF (ISTAK.EQ.IDRES.AND.IDG(I).EQ.50116) THEN
924 DISTMIN=(XTRG(I)-XGEANT(I))**2+(YTRG(I)-YGEANT(I))**2
925 IF (DISTMIN.LT.DIST(1)) THEN
930 NHITMUON(1,NMUON(1))=I
932 IF (ISTAK.EQ.IDRES.AND.IDG(I).EQ.70116) THEN
933 DISTMIN=(XTRG(I)-XGEANT(I))**2+(YTRG(I)-YGEANT(I))**2
934 IF (DISTMIN.LT.DIST(2)) THEN
939 NHITMUON(2,NMUON(2))=I
945 IF (NMUON(J).GE.2) THEN
946 print*,'j=',j,' nmuon=',nmuon(j)
949 IF (NHITMUON(J,K).NE.NMUONGOOD(J)) IDG(NHITMUON(J,K))=999 ! flag les mauvais hits MUONS
956 * 1 ere boucle Lecture des hits initiaux
958 DO I = 1,NHITTOT1 ! Boucle sur les hits GEANT de toutes les ch.
966 IF (R.LT.RMIN(ICH).OR.R.GT.RMAX1(ICH)) then
967 if (ich.le.10.and.i.le.28) then
968 print*,'****** chambre ',ich,' ****** hit ',i
969 print*,'ityp=',itypg(i)
970 print*,'x=',XTRG(I),' y=',YTRG(I)
971 print*,'R=',R,' RMIN=',RMIN(ICH),' RMAX1=',RMAX1(ICH)
978 NHIT(ICH) = NHIT(ICH)+1
979 IH(ICH,NHIT(ICH)) = I
980 XMA(ICH,NHIT(ICH)) = XTRG(I)
981 YMA(ICH,NHIT(ICH)) = YTRG(I)
982 IMARK(ICH,NHIT(ICH)) = 0
984 print*,' XTRG(I)=', XTRG(I),' YTRG(I)=', YTRG(I),' IDG(I)=',
991 * Association des hits entre chambres d'une station
992 II = 0 ! nombre de hits GEANT par station
994 IZ = INT(FLOAT(ICH1+1)/2.)
1009 PVERT1(II) = PVERT1G(I)
1010 PVERT2(II) = PVERT2G(I)
1011 PVERT3(II) = PVERT3G(I)
1012 ZVERT(II) = ZVERTG(I)
1014 IF (IEFFI(I).EQ.1) THEN
1018 XEXT1 = (ZPLANE(IZ)-DZ_PL(IZ))/ZPLANE(IZ)*X1
1019 YEXT1 = (ZPLANE(IZ)-DZ_PL(IZ))/ZPLANE(IZ)*Y1
1021 PRINT *,'***** DEBUT RECHERCHE',' ID1=',ID1,' ich1=',ICH1
1022 PRINT *,' XTR(II)=', XTR(II),' YTR(II)=', YTR(II)
1023 PRINT *,' ITYP(II)=', ITYP(II)
1024 DO I2 = 1,NHIT(ICH2)
1026 IF (IEFFI(J).EQ.1) THEN
1030 IF (DX.LT.DXMAX(IZ)) THEN
1034 print *,' DX=',DX,' KC=',KC,' ID2=',ID2
1047 IF (DY.LT.DYOLD.AND.DY.LT.(SIGCUT*DYMAX(IZ))) THEN
1050 PRINT *,' ID2=',ID2,' DY=',DY
1053 IF (I2FIND.GT.0) THEN
1059 JHIT(IZ) = JHIT(IZ)+1
1062 XM(IZ,JHIT(IZ)) = X1
1063 YM(IZ,JHIT(IZ)) = Y1
1064 IZM(IZ,JHIT(IZ)) = 1
1065 PHM(IZ,JHIT(IZ)) = -ATAN((X2-X1)/DZ_PL(IZ))
1066 ALM(IZ,JHIT(IZ)) = ATAN((Y2-Y1)/DZ_PL(IZ)*
1067 & COS(PHM(IZ,JHIT(IZ))))
1068 XMR(IZ,JHIT(IZ),1) = X1
1069 YMR(IZ,JHIT(IZ),1) = Y1
1070 XMR(IZ,JHIT(IZ),2) = X2
1071 YMR(IZ,JHIT(IZ),2) = Y2
1072 IP(IZ,JHIT(IZ)) = II
1075 ISTAK = MOD(ISTAK,30000)
1076 ISTAK = MOD(ISTAK,10000)
1078 IF (ISTAK.EQ.IDRES.AND.ID1.NE.999) THEN
1085 PVERT1(II) = PVERT1G(J)
1086 PVERT2(II) = PVERT2G(J)
1087 PVERT3(II) = PVERT3G(J)
1088 ZVERT(II) = ZVERTG(J)
1095 IF (IFIND.EQ.0) THEN
1096 JHIT(IZ) = JHIT(IZ)+1
1097 XM(IZ,JHIT(IZ)) = X1
1098 YM(IZ,JHIT(IZ)) = Y1
1099 IZM(IZ,JHIT(IZ)) = 1
1100 IP(IZ,JHIT(IZ)) = II
1101 PHM(IZ,JHIT(IZ)) = 10.
1102 ALM(IZ,JHIT(IZ)) = 10.
1103 XMR(IZ,JHIT(IZ),1) = X1
1104 YMR(IZ,JHIT(IZ),1) = Y1
1105 XMR(IZ,JHIT(IZ),2) = 0.
1106 YMR(IZ,JHIT(IZ),2) = 0.
1108 CALL CHFILL2(1000+IZ,SNGL(X1),SNGL(Y1),R2)
1113 * On conserve les HITS de la 2nde chambre des stations
1116 IZ = INT(FLOAT(ICH+1)/2.)
1120 IF (IMARK(ICH,I).EQ.0) THEN
1121 ** print *,' ich=',ich,' i=',i,' j=',j
1131 PVERT1(II) = PVERT1G(J)
1132 PVERT2(II) = PVERT2G(J)
1133 PVERT3(II) = PVERT3G(J)
1134 ZVERT(II) = ZVERTG(I)
1136 IF (IEFFI(J).EQ.1) THEN
1137 JHIT(IZ) = JHIT(IZ)+1
1138 ** XM(IZ,JHIT(IZ)) = ZPLANE(IZ)/(ZPLANE(IZ)-DZ_PL)
1140 ** YM(IZ,JHIT(IZ)) = ZPLANE(IZ)/(ZPLANE(IZ)-DZ_PL)
1142 XM(IZ,JHIT(IZ)) = XMA(ICH,I)
1143 YM(IZ,JHIT(IZ)) = YMA(ICH,I)
1144 IZM(IZ,JHIT(IZ)) = 2
1145 PHM(IZ,JHIT(IZ)) = 10.
1146 ALM(IZ,JHIT(IZ)) = 10.
1147 IP(IZ,JHIT(IZ)) = II
1148 XMR(IZ,JHIT(IZ),1) = 1000.
1149 YMR(IZ,JHIT(IZ),1) = 1000.
1150 XMR(IZ,JHIT(IZ),2) = XMA(ICH,I)
1151 YMR(IZ,JHIT(IZ),2) = YMA(ICH,I)
1160 IF (IDEBUG.GE.2) THEN
1161 PRINT *,'TRACKF_MICRO nb hits:',NHITTOT
1164 PRINT *,' IZ=',IZ,' JHIT(IZ)=',JHIT(IZ)
1167 PRINT *,' ID(II)=',ID(II)
1168 PRINT *,' XMR(IZ,J,1)=', XMR(IZ,J,1),
1169 & ' YMR(IZ,J,1)=', YMR(IZ,J,1)
1170 PRINT *,' XMR(IZ,J,2)=', XMR(IZ,J,2),
1171 & ' YMR(IZ,J,2)=', YMR(IZ,J,2)
1178 *************************************************************************
1179 SUBROUTINE RECO_TRACKF(IDRES,IREADGEANT)
1180 *************************************************************************
1183 *************************************************************************
1184 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
1186 PARAMETER (MAXIDG=20000)
1188 PARAMETER (MAXHITCH=10000,MAXTRK=50000,MAXHITTOT=20000,
1189 & NBSTATION=5,MAXCAN=1000,NBCHAMBER=10,NTRMAX=500)
1191 COMMON/CHHIT/XM(NBSTATION,MAXHITCH),YM(NBSTATION,MAXHITCH),
1192 & PHM(NBSTATION,MAXHITCH),ALM(NBSTATION,MAXHITCH),
1193 & IZM(NBSTATION,MAXHITCH),
1194 & IP(NBSTATION,MAXHITCH),JHIT(NBSTATION),
1195 & XMR(NBSTATION,MAXHITCH,2),YMR(NBSTATION,MAXHITCH,2)
1197 COMMON/RHIT/ITYP(MAXHITTOT),XTR(MAXHITTOT),YTR(MAXHITTOT),
1198 & PTOT(MAXHITTOT),ID(MAXHITTOT),IZST(MAXHITTOT),
1199 & PVERT1(MAXHITTOT),PVERT2(MAXHITTOT),PVERT3(MAXHITTOT),
1200 & ZVERT(MAXHITTOT),NHITTOT
1203 COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
1205 COMMON/VERIFGEANT/ITTROUGH(MAXTRK,NBSTATION),
1206 & IT_LIST(MAXTRK),IT_NP(MAXTRK),ITCHECK(MAXTRK),
1209 COMMON/HCHHIT/HHIT(MAXHITCH),INDEXTAB(MAXHITCH),INDEXMAX
1211 COMMON/PRECIS/EEXM(NBSTATION),EEYM(NBSTATION),EEPH(NBSTATION),
1214 COMMON/MEASUR/XMES(NBSTATION),YMES(NBSTATION),IZMES(NBSTATION),
1215 & PHMES(NBSTATION),ALMES(NBSTATION),MPOS(NBSTATION),
1218 COMMON/PLANE/XPL(NBSTATION,2),YPL(NBSTATION,2),
1219 & PHPL(NBSTATION),ALPL(NBSTATION),CHI2PL
1221 COMMON/CANDIDAT/JCAN(NBSTATION,MAXCAN),JCANTYP(MAXCAN),
1222 & EEX(MAXCAN),EEY(MAXCAN),EEP(MAXCAN),EEA(MAXCAN)
1224 COMMON /VERTEX/ERRV,IVERTEX
1226 COMMON/TRACKSUM/NRES(5),NRESF,NTRMUALL,NMUONALL,NGHOSTALL,
1227 & NTRACKFALL,NERRALL(NBSTATION),IR
1229 COMMON/TRACKFOUT/IEVOUT,NTREVT,JJOUT(NBCHAMBER,NTRMAX),
1230 & ISTAT(NTRMAX),PXZOUT(NTRMAX),TPHIOUT(NTRMAX),
1231 & TALAMOUT(NTRMAX),XVERTOUT(NTRMAX),YVERTOUT(NTRMAX),
1233 & XMESOUT(NBCHAMBER,NTRMAX),YMESOUT(NBCHAMBER,NTRMAX)
1234 & ,PXVOUT(NTRMAX),PYVOUT(NTRMAX),PZVOUT(NTRMAX)
1236 COMMON/TRACKFI/EFF,EFF1,EFF2,XPREC,YPREC,PHIPREC,ALAMPREC,
1237 & HCUT,LBKG,SIGCUT,ALPHATOP,HTOP
1239 COMMON/DEBEVT/IDEBUG
1240 common/dstation/idstation
1242 DIMENSION PEST(5),PSTEP(5),NERR(NBSTATION),IFIND2(10),NMU45(2)
1243 DIMENSION NMU345(2),NMUONF(2)
1248 ** GEANT informations
1250 IT_LIST(I)= 0 ! ID of the GEANT tracks
1253 ITRACK(I) = 0 ! Track number associated to hit number I
1257 * BOUCLE SUR LES HITS
1260 IF (ID(IH).EQ.IT_LIST(IT)) THEN
1265 NTRTOT = NTRTOT +1 ! NB de trace GEANT
1266 IT_LIST(NTRTOT) = ID(IH) ! ID DE LA TRACE NTRTOT
1267 IT_NP(NTRTOT) = 0 ! NOMBRE DE HITS PAR TRACE
1269 ITTROUGH(NTRTOT,II) = 0
1272 4 IT_NP(IT1) = IT_NP(IT1)+1 ! Number of crossed stations per track
1273 ITTROUGH(IT1,IZST(IH))=IH ! =IH si la trace IT touche le plan IZST
1276 IF (IDEBUG.GE.2) THEN
1277 PRINT *,'RECO_TRACKF nb total de trace GEANT:',NTRTOT
1285 IF ((ITTROUGH(IT,1)*ITTROUGH(IT,2)*ITTROUGH(IT,3)*
1286 & ITTROUGH(IT,4)*ITTROUGH(IT,5)).NE.0)
1287 & THEN ! track crossing all stations
1289 IH = ITTROUGH(IT,NBSTATION)
1290 IF (ITYP(IH).EQ.5.OR.ITYP(IH).EQ.6) THEN
1292 ISTAK = MOD(ISTAK,30000)
1293 ISTAK = MOD(ISTAK,10000)
1295 pt=sqrt(pvert1(ih)**2+pvert2(ih)**2)
1296 thet=datan2(pt,pvert3(ih))*180/3.1416
1297 pp=sqrt(pt**2+pvert3(ih)**2)
1299 IF (ISTAK.EQ.IDRES) THEN ! psi or upsilon
1301 NTRMUALL = NTRMUALL+1
1308 IF (IDEBUG.GE.2) THEN
1309 PRINT *,'RECO_TRACKF nb of part. GEANT crossing 5 st.:',
1311 PRINT *,'RECO_TRACKF nb of muons/res. GEANT crossing 5 st.:',
1315 ** CALL H_ACCEPTANCE(5)
1316 ** CALL H_ACCEPTANCE(4)
1322 CALL ORDONNE_HIT(5,HCUT)
1326 X1 = XM(5,JJ)-(ZPLANE(5)-ZPLANE(4))*TAN(PHM(5,JJ))
1327 Y1 = YM(5,JJ)+(ZPLANE(5)-ZPLANE(4))*TAN(ALM(5,JJ))
1329 X2 = XM(5,JJ)-(ZPLANE(5)-ZPLANE(4)+DZ_PL(4))*TAN(PHM(5,JJ))
1330 Y2 = YM(5,JJ)+(ZPLANE(5)-ZPLANE(4)+DZ_PL(4))*TAN(ALM(5,JJ))
1333 * Domaine de recherche dans la st. 4
1334 HCONST = 0.0136*SQRT(0.06)/HTOP ! -> X/X0=0.03 % / chamber
1335 EPH2 = 2.0*PHIPREC**2 + (HCONST*HHIT(JJ))**2
1336 EAL2 = 2.0*ALAMPREC**2 + (HCONST*HHIT(JJ))**2
1337 EXM12 = (PHIPREC**2+(HCONST*HHIT(JJ))**2)*
1338 & (ZPLANE(5)-ZPLANE(4))**2 + XPREC**2
1339 EYM12 = (ALAMPREC**2+(HCONST*HHIT(JJ))**2)*
1340 & (ZPLANE(5)-ZPLANE(4))**2 + YPREC**2
1341 EXM22 = (PHIPREC**2+(HCONST*HHIT(JJ))**2)*
1342 & (ZPLANE(5)-ZPLANE(4)+DZ_PL(4))**2 + XPREC**2
1343 EYM22 = (ALAMPREC**2+(HCONST*HHIT(JJ))**2)*
1344 & (ZPLANE(5)-ZPLANE(4)+DZ_PL(4))**2 + YPREC**2
1346 EPH = SIGCUT*SQRT(EPH2)
1347 EAL = SIGCUT*SQRT(EAL2)
1348 EX1 = SIGCUT*SQRT(EXM12)
1349 EY1 = SIGCUT*SQRT(EYM12)
1350 EX2 = SIGCUT*SQRT(EXM22)
1351 EY2 = SIGCUT*SQRT(EYM22)
1353 ** P2 = (HTOP/HHIT(JJ))**2
1355 ** EPH=SIGCUT*SQRT(9.14D-7+1.2D-3/P2)
1356 ** EAL=SIGCUT*SQRT(1.84D-4)
1357 ** EX1=SIGCUT*SQRT(1.95D-2+6.37/P2)
1358 ** EY1=SIGCUT*SQRT(3.89+151./P2)
1362 * renvoie le num de hit de 4 le plus pres dans le domaine de recherche
1364 CALL DISTMIN4(X1,Y1,PHM(5,JJ),ALM(5,JJ),4,EX1,EY1,EPH,EAL,
1367 CALL CHECK_HISTO4(11,4,IFIND,5,JJ,X1,Y1,PHM(5,JJ),ALM(5,JJ),P1,
1369 IF (IFIND.GT.0) THEN
1370 CALL STOCK_CANDIDAT(5,JJ,4,IFIND,IFIND2,EX1,EY1,EPH,EAL,
1373 CALL DISTMIN2(X1,Y1,X2,Y2,4,EX1,EY1,EX2,EY2,IFIND,IFIND2)
1374 CALL CHECK_HISTO2(0,4,IFIND,5,JJ,X1,Y1,X2,Y2,P1,EX1,EY1,
1376 IF (IFIND.GT.0) THEN
1377 CALL STOCK_CANDIDAT(5,JJ,4,IFIND,IFIND2,EX1,EY1,EPH,EAL,
1385 CALL ORDONNE_HIT(4,HCUT)
1389 X1 = XM(4,JJ)-(ZPLANE(4)-ZPLANE(5))*TAN(PHM(4,JJ))
1390 Y1 = YM(4,JJ)+(ZPLANE(4)-ZPLANE(5))*TAN(ALM(4,JJ))
1392 X2 = XM(4,JJ)-(ZPLANE(4)-ZPLANE(5)+DZ_PL(5))*TAN(PHM(4,JJ))
1393 Y2 = YM(4,JJ)+(ZPLANE(4)-ZPLANE(5)+DZ_PL(5))*TAN(ALM(4,JJ))
1395 * Domaine de recherche dans la st. 5
1396 HCONST = 0.0136*SQRT(0.06)/HTOP ! -> X/X0=0.03 / chamber
1397 ** EPH2 = 2.0*PHIPREC**2 + (HCONST*HHIT(JJ))**2
1398 ** EAL2 = 2.0*ALAMPREC**2 + (HCONST*HHIT(JJ))**2
1399 ** EX12 = 2.0*(XPREC/SQRT(2.))**2
1400 ** & + (ZPLANE(5)-ZPLANE(4))**2/2.*EPH2
1401 ** EY12 = 2.0*YPREC**2 + (ZPLANE(5)-ZPLANE(4))**2/2.*EAL2
1402 ** EX22 = 2.0*(XPREC/SQRT(2.))**2
1403 ** & + (ZPLANE(5)-ZPLANE(4)-DZ_PL(5))**2/2.*EPH2
1404 ** EY22 = 2.0*YPREC**2 + (ZPLANE(5)-ZPLANE(4)-DZ_PL(5))**2/2.
1406 EPH2 = 2.0*PHIPREC**2 + (HCONST*HHIT(JJ))**2
1407 EAL2 = 2.0*ALAMPREC**2 + (HCONST*HHIT(JJ))**2
1408 EXM12 = (PHIPREC**2+(HCONST*HHIT(JJ))**2)*
1409 & (ZPLANE(5)-ZPLANE(4))**2 + XPREC**2
1410 EYM12 = (ALAMPREC**2+(HCONST*HHIT(JJ))**2)*
1411 & (ZPLANE(5)-ZPLANE(4))**2 + YPREC**2
1412 EXM22 = (PHIPREC**2+(HCONST*HHIT(JJ))**2)*
1413 & (ZPLANE(5)-ZPLANE(4)-DZ_PL(5))**2 + XPREC**2
1414 EYM22 = (ALAMPREC**2+(HCONST*HHIT(JJ))**2)*
1415 & (ZPLANE(5)-ZPLANE(4)-DZ_PL(5))**2 + YPREC**2
1417 EPH = SIGCUT*SQRT(EPH2)
1418 EAL = SIGCUT*SQRT(EAL2)
1419 EX1 = SIGCUT*SQRT(EXM12)
1420 EY1 = SIGCUT*SQRT(EYM12)
1421 EX2 = SIGCUT*SQRT(EXM22)
1422 EY2 = SIGCUT*SQRT(EYM22)
1425 ** P2 = (HTOP/HHIT(JJ))**2
1427 ** EPH=SIGCUT*SQRT(9.14D-7+1.2D-3/P2)
1428 ** EAL=SIGCUT*SQRT(1.84D-4)
1429 ** EX1=SIGCUT*SQRT(1.95D-2+6.37/P2)
1430 ** EY1=SIGCUT*SQRT(3.89+151./P2)
1433 * Renvoie le num de hit de 5 le plus pres dans le domaine de recherche
1434 CALL DISTMIN2(X1,Y1,X2,Y2,5,EX1,EY1,EX2,EY2,IFIND,IFIND2)
1436 CALL CHECK_HISTO2(0,5,IFIND,4,JJ,X1,Y1,X2,Y2,P1,EX1,EY1,
1438 IF (IFIND.GT.0) THEN
1440 IF (IFIND.EQ.JCAN(5,ICAN).AND.JJ.EQ.JCAN(4,ICAN)) GOTO 40
1441 ** IF (JJ.EQ.JCAN(4,ICAN).AND.
1442 ** & ABS(XM(5,JCAN(5,ICAN))-XM(5,IFIND)).LT.XPREC/10.)
1443 ** & GO TO 40 ! elimine les doubles comptages de traces ccc
1444 ** IF (IFIND.EQ.JCAN(5,ICAN).AND.
1445 ** & ABS(XM(4,JCAN(4,ICAN))-XM(4,JJ)).LT.XPREC/10.)
1446 ** & GO TO 40 ! elimine les doubles comptages de traces ccc
1448 DIST1 = SQRT(((XM(5,JCAN(5,ICAN))-XM(5,IFIND))
1449 & /(0.1*XPREC))**2+((YM(5,JCAN(5,ICAN))-YM(5,IFIND))
1451 DIST2 = SQRT(((XM(4,JCAN(4,ICAN))-XM(4,JJ))
1452 & /(0.1*XPREC))**2+((YM(4,JCAN(4,ICAN))-YM(4,JJ))
1454 IF (DIST1.LT.2..AND.DIST2.LT.2.)
1455 & GO TO 40 ! elimine les doubles comptages de traces ccc
1457 CALL STOCK_CANDIDAT(4,JJ,5,IFIND,IFIND2,EX1,EY1,EPH,EAL,NCAN
1469 IF (JJ4.GT.0.AND.JJ5.GT.0) THEN
1472 IT = ITRACK(IP(5,JJ5))
1473 IF (ITCHECK(IT).EQ.1) THEN
1474 IF (ID4.EQ.ID5) THEN
1475 IF (ITYP(IP(5,JJ5)).EQ.5) NMU45(1) = 1
1476 IF (ITYP(IP(5,JJ5)).EQ.6) NMU45(2) = 1
1481 IF (NMU45(1).GE.1.AND.NMU45(2).EQ.1) NRES(1) = NRES(1)+1
1483 IF (IDEBUG.GE.2) THEN
1484 PRINT *,'RECO_TRACKF nb candidat recherche 4->5 et 5->4 :'
1486 PRINT *,'RECO_TRACKF nb of good muons 4->5 et 5->4 :'
1487 & ,(NMU45(1)+NMU45(2))
1495 * -- Boucle sur les candidats (4,5) NCAN
1511 DO ICH = 1,NBSTATION
1519 IF (JCANTYP(ICAN).EQ.2) MANG(4)=0
1520 IF (JCANTYP(ICAN).EQ.3) MANG(5)=0
1521 EEXM(5) = EEX(ICAN)/SIGCUT
1522 EEYM(5) = EEY(ICAN)/SIGCUT
1523 EEPH(5) = EEP(ICAN)/SIGCUT
1524 EEAL(5) = EEA(ICAN)/SIGCUT
1525 EEXM(4) = EEX(ICAN)/SIGCUT
1526 EEYM(4) = EEY(ICAN)/SIGCUT
1527 EEPH(4) = EEP(ICAN)/SIGCUT
1528 EEAL(4) = EEA(ICAN)/SIGCUT
1532 IF (IZM(4,JJ4).EQ.1) THEN
1535 ZPL4 = ZPLANE(4)-DZ_PL(4)
1537 IF (IZM(5,JJ5).EQ.1) THEN
1540 ZPL5 = ZPLANE(5)-DZ_PL(5)
1542 TPHEST = (XM(5,JJ5) - XM(4,JJ4))/(ZPL5-ZPL4)
1543 PHEST = ATAN(TPHEST)
1544 TALEST = -(YM(5,JJ5) - YM(4,JJ4))*COS(PHEST)
1546 PXZEST = -HTOP/(XM(4,JJ4) - ZPL4*TPHEST)
1547 PEST(1) = 1.0/PXZEST
1548 PEST(2) = TPHEST - ALPHATOP/PXZEST ! PHI emission au vertex
1549 PEST(3) = TALEST ! tan(lambda) !h=zm*ang deviation alpha
1550 PEST(4) = 0.0 !alpha=qbl/pxz
1552 PSTEP(1) = 0.003 ! =d(1/p)=delta(p)/p**2
1553 PSTEP(2) = 0.001 ! 0.5 degre
1554 PSTEP(3) = 0.001 ! 0.5 degre
1559 IZMES(4) = IZM(4,JJ4)
1560 PHMES(4) = PHM(4,JJ4)
1561 ALMES(4) = ALM(4,JJ4)
1564 IZMES(5) = IZM(5,JJ5)
1565 PHMES(5) = PHM(5,JJ5)
1566 ALMES(5) = ALM(5,JJ5)
1568 IVERTEX = 0 ! Vertex = (0,0,0)
1570 * -- Fit 4,5,V pour la recherche dans 3
1572 CALL TRACKF_FIT(IVERTEX,PEST,PSTEP,PXZINV45,TPHI45,TALAM45,
1575 * -- Recherche dans la station 3
1577 P2 = (1.0D0 + TALAM45**2)/(PXZINV45**2) ! P**2
1579 if (idstation.eq.8) then ! DZ_CH = 8 CM
1580 EPH=SIGCUT*SQRT(3.6D-6+0.011/P2)
1581 EAL=SIGCUT*SQRT(6.85D-4)
1582 EXM=SIGCUT*SQRT(0.034+446./P2)
1583 EYM=SIGCUT*SQRT(0.049+354./P2)
1584 else if (idstation.eq.20) then ! DZ_CH = 20 CM
1585 EPH=SIGCUT*SQRT(4.1D-7+0.015/P2)
1586 EAL=SIGCUT*SQRT(1.04D-4)
1587 EXM=SIGCUT*SQRT(0.0+459./P2)
1588 EYM=SIGCUT*SQRT(0.042+345./P2)
1591 EEXM(3) = EXM/SIGCUT
1592 EEYM(3) = EYM/SIGCUT
1593 EEPH(3) = EPH/SIGCUT
1594 EEAL(3) = EAL/SIGCUT
1596 ** EEXM(IEX) = EEXM(3)
1597 ** EEYM(IEX) = EEYM(3)
1598 ** EEPH(IEX) = EEPH(3)
1599 ** EEAL(IEX) = EEAL(3)
1608 CALL DISTMIN4(X1,Y1,PHI1,ALAM1,3,EXM,EYM,EPH,EAL,IPA3,IFIND2)
1610 P1 = PTOT(IP(5,JJ5))
1612 CALL CHECK_HISTO4(61,3,IPA3,5,JJ5,X1,Y1,PHI1,ALAM1,P1,
1619 IZMES(3) = IZM(3,JJ3)
1620 PHMES(3) = PHM(3,JJ3)
1621 ALMES(3) = ALM(3,JJ3)
1626 CALL DISTMIN2(X1,Y1,X2,Y2,3,EXM,EYM,0.D0,0.D0,IP3,IFIND2)
1627 CALL CHECK_HISTO2(0,3,IP3,5,JJ5,X1,Y1,X2,Y2,P1,EXM,EYM,
1634 IZMES(3) = IZM(3,JJ3)
1642 * -- Fit 3,4,5 pour la recherche dans 2
1646 IF (JJ5.GT.0.AND.JJ4.GT.0.AND.JJ3.GT.0) THEN
1650 IT = ITRACK(IP(5,JJ5))
1651 IF (ITCHECK(IT).EQ.1) THEN
1652 IF ((ID5.EQ.ID3).AND.(ID5.EQ.ID4)) THEN
1653 IF (ITYP(IP(5,JJ5)).EQ.5) NMU345(1) = 1
1654 IF (ITYP(IP(5,JJ5)).EQ.6) NMU345(2) = 1
1672 CALL TRACKF_FIT(IVERTEX,PEST,PSTEP,PXZINV345,TPHI345,
1673 & TALAM345,XVERT,YVERT)
1675 * -- Recherche dans la st. 2
1677 P2 = (1.0D0 + TALAM345**2)/(PXZINV345**2) ! P**2
1679 if (idstation.eq.8) then ! DZ_CH = 8 CM
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 else if (idstation.eq.20) then ! DZ_CH = 20 CM
1685 EPH=SIGCUT*SQRT(5.78D-7+5.97D-3/P2)
1686 EAL=SIGCUT*SQRT(1.D-4)
1687 EXM=SIGCUT*SQRT(0.+1453./P2)
1688 EYM=SIGCUT*SQRT(2.78D-2+504./P2)
1691 EEXM(2) = EXM/SIGCUT
1692 EEYM(2) = EYM/SIGCUT
1693 EEPH(2) = EPH/SIGCUT
1694 EEAL(2) = EAL/SIGCUT
1696 ** EEXM(IEX) = EEXM(2)
1697 ** EEYM(IEX) = EEYM(2)
1698 ** EEPH(IEX) = EEPH(2)
1699 ** EEAL(IEX) = EEAL(2)
1706 CALL DISTMIN4(X1,Y1,PHI1,ALAM1,2,EXM,EYM,EPH,EAL,IPA2,IFIND2)
1707 P1 = PTOT(IP(5,JJ5))
1708 CALL CHECK_HISTO4(21,2,IPA2,5,JJ5,X1,Y1,PHI1,ALAM1,P1,
1710 * -- Recherche dans la st. 1
1712 if (idstation.eq.8) then ! DZ_CH = 8 CM
1713 EPH=SIGCUT*SQRT(3.54D-6+4.49D-3/P2)
1714 EAL=SIGCUT*SQRT(6.14D-4)
1715 EXM=SIGCUT*SQRT(6.43D-2+986./P2)
1716 EYM=SIGCUT*SQRT(4.66D-2+1444./P2)
1717 else if (idstation.eq.20) then ! DZ_CH = 20 CM
1718 EPH=SIGCUT*SQRT(4.96D-7+5.87D-3/P2)
1719 EAL=SIGCUT*SQRT(1.D-4)
1720 EXM=SIGCUT*SQRT(6.98D-4+1467./P2)
1721 EYM=SIGCUT*SQRT(5.22D-2+1013./P2)
1724 EEXM(1) = EXM/SIGCUT
1725 EEYM(1) = EYM/SIGCUT
1726 EEPH(1) = EPH/SIGCUT
1727 EEAL(1) = EAL/SIGCUT
1729 ** EEXM(IEX) = EEXM(1)
1730 ** EEYM(IEX) = EEYM(1)
1731 ** EEPH(IEX) = EEPH(1)
1732 ** EEAL(IEX) = EEAL(1)
1739 CALL DISTMIN4(X1,Y1,PHI1,ALAM1,1,EXM,EYM,EPH,EAL,IPA1,IFIND2)
1740 CALL CHECK_HISTO4(31,1,IPA1,5,JJ5,X1,Y1,PHI1,ALAM1,P1,
1742 * -- A partir de P+A dans la st. 2
1755 XMES(2) = XM(2,IPA2)
1756 YMES(2) = YM(2,IPA2)
1757 IZMES(2) = IZM(2,IPA2)
1758 PHMES(2) = PHM(2,IPA2)
1759 ALMES(2) = ALM(2,IPA2)
1763 * -- Fit V,2,3,4,5 pour la recherche dans 1
1765 CALL TRACKF_FIT(IVERTEX,PEST,PSTEP,PXZINV,TPHI,TALAM,
1768 * !!! ATTENTION aux erreurs
1770 P2 = (1.0D0 + TALAM**2)/(PXZINV**2)
1772 if (idstation.eq.8) then ! DZ_CH = 8 CM
1773 EPH=SIGCUT*SQRT(3.15D-6+9.21D-3/P2)
1774 EAL=SIGCUT*SQRT(5.94D-4)
1775 EXM=SIGCUT*SQRT(4.16D-2+182./P2)
1776 EYM=SIGCUT*SQRT(2.13)
1777 else if (idstation.eq.20) then ! DZ_CH = 20 CM
1778 EPH=SIGCUT*SQRT(9.58D-7+8.93D-3/P2)
1779 EAL=SIGCUT*SQRT(1.D-4)
1780 EXM=SIGCUT*SQRT(1.92D-2+103.3/P2)
1781 EYM=SIGCUT*SQRT(6.3D-2+81.1/P2)
1784 EEXM(1) = EXM/SIGCUT
1785 EEYM(1) = EYM/SIGCUT
1786 EEPH(1) = EPH/SIGCUT
1787 EEAL(1) = EAL/SIGCUT
1789 ** EEXM(IEX) = EEXM(1)
1790 ** EEYM(IEX) = EEYM(1)
1791 ** EEPH(IEX) = EEPH(1)
1792 ** EEAL(IEX) = EEAL(1)
1801 CALL DISTMIN4(X1,Y1,PHI1,ALAM1,1,EXM,EYM,EPH,EAL,
1803 CALL CHECK_HISTO4(41,1,IPA2PA1,5,JJ5,X1,Y1,PHI1,ALAM1,
1804 & P1,EXM,EYM,EPH,EAL)
1805 IF (IPA2PA1.GT.0) THEN
1812 IZMES(1) = IZM(1,JJ1)
1813 PHMES(1) = PHM(1,JJ1)
1814 ALMES(1) = ALM(1,JJ1)
1819 CALL DISTMIN2(X1,Y1,X2,Y2,1,EXM,EYM,0.D0,0.D0,IPA2P1,
1821 CALL CHECK_HISTO2(0,1,IPA2P1,5,JJ5,X1,Y1,X2,Y2,P1,
1822 & EXM,EYM,0.D0,0.D0)
1825 * -- A partir de P+A dans la st. 1
1838 XMES(1) = XM(1,IPA1)
1839 YMES(1) = YM(1,IPA1)
1840 IZMES(1) = IZM(1,IPA1)
1841 PHMES(1) = PHM(1,IPA1)
1842 ALMES(1) = ALM(1,IPA1)
1848 * -- Fit V,1,3,4,5 pour la recherche dans 2
1850 CALL TRACKF_FIT(IVERTEX,PEST,PSTEP,PXZINV,TPHI,TALAM,
1853 * !!! ATTENTION aux erreurs
1855 P2 = (1.0D0 + TALAM**2)/(PXZINV**2)
1857 if (idstation.eq.8) then ! DZ_CH = 8 CM
1858 EPH=SIGCUT*SQRT(3.15D-6+9.21D-3/P2)
1859 EAL=SIGCUT*SQRT(5.94D-4)
1860 EXM=SIGCUT*SQRT(4.16D-2+182./P2)
1861 EYM=SIGCUT*SQRT(2.13)
1862 else if (idstation.eq.20) then ! DZ_CH = 20 CM
1863 EPH=SIGCUT*SQRT(9.58D-7+8.93D-3/P2)
1864 EAL=SIGCUT*SQRT(1.D-4)
1865 EXM=SIGCUT*SQRT(4.0D-2+11.4/P2)
1866 EYM=SIGCUT*SQRT(5.5D-2+44.35/P2)
1869 EEXM(2) = EXM/SIGCUT
1870 EEYM(2) = EYM/SIGCUT
1871 EEPH(2) = EPH/SIGCUT
1872 EEAL(2) = EAL/SIGCUT
1874 ** EEXM(IEX) = EEXM(2)
1875 ** EEYM(IEX) = EEYM(2)
1876 ** EEPH(IEX) = EEPH(2)
1877 ** EEAL(IEX) = EEAL(2)
1887 CALL DISTMIN4(X1,Y1,PHI1,ALAM1,2,EXM,EYM,EPH,EAL,
1890 CALL CHECK_HISTO4(51,2,IPA1PA2,5,JJ5,X1,Y1,PHI1,ALAM1,
1891 & P1,EXM,EYM,EPH,EAL)
1892 IF (IPA1PA2.GT.0) THEN
1897 IZMES(2) = IZM(2,JJ2)
1898 PHMES(2) = PHM(2,JJ2)
1899 ALMES(2) = ALM(2,JJ2)
1904 CALL DISTMIN2(X1,Y1,X2,Y2,2,EXM,EYM,0.D0,0.D0,IPA1P2,
1906 CALL CHECK_HISTO2(0,2,IPA1P2,5,JJ5,X1,Y1,X2,Y2,P1,
1907 & EXM,EYM,0.D0,0.D0)
1910 * -- Selection par type de candidat
1911 IF (IPA1PA2.EQ.0.OR.IPA2PA1.EQ.0) THEN
1912 IF (IPA2.GT.0.AND.IPA2P1.GT.0) THEN
1917 IZMES(1) = IZM(1,JJ1)
1922 ELSEIF(IPA1.GT.0.AND.IPA1P2.GT.0) THEN
1927 IZMES(2) = IZM(2,JJ2)
1937 NTRACKFOLD = NTRACKF
1938 NMUONFOLD = NMUONF(1)+NMUONF(2)
1940 CALL CHECK_MUON(JJ1,JJ2,JJ3,JJ4,JJ5,NTRACKF,NMUONF,
1941 & NGHOST,NERR,PXZINV,TPHI,TALAM,XVERT,YVERT)
1943 IF (NTRACKF.GT.NTRACKFOLD) THEN
1945 IF ((NMUONF(1)+NMUONF(2)).GT.NMUONFOLD) ISTAT(NTRACKF) = 1 ! Good muon
1946 IF (NGHOST.GT.NGHOSTOLD) ISTAT(NTRACKF) = 2 ! ghost track
1947 PXZOUT(NTRACKF) = 1./PXZINV
1948 TPHIOUT(NTRACKF) = TPHI
1949 TALAMOUT(NTRACKF) = TALAM
1950 XVERTOUT(NTRACKF) = XVERT
1951 YVERTOUT(NTRACKF) = YVERT
1952 PXVOUT(NTRACKF) = PVERT1(IP(5,JJ5))
1953 PYVOUT(NTRACKF) = PVERT2(IP(5,JJ5))
1954 PZVOUT(NTRACKF) = PVERT3(IP(5,JJ5))
1955 CHI2OUT(NTRACKF) = CHI2PL
1958 ENDDO ! end loop on candidats NCAN
1960 125 IF (NMU345(1).EQ.1.AND.NMU345(2).EQ.1) NRES(2) = NRES(2) + 1
1962 IF (IDEBUG.GE.2) THEN
1963 PRINT *,'RECO_TRACKF nb of good muons 3 4 5 :'
1964 & ,(NMU345(1)+NMU345(2))
1968 IF (IDEBUG.GE.2) THEN
1969 PRINT *,'RECO_TRACKF nb of track/evt :',NTRACKF
1970 PRINT *,'RECO_TRACKF nb of good muon/evt :',NMUONF(1)+NMUONF(2)
1971 PRINT *,'RECO_TRACKF nb of ghost/evt :',NGHOST
1973 PRINT *,'RECO_TRACKF nb of error in st',I,'/evt:',NERR(I)
1976 IF (NTRMU.GE.2) NRES(5) = NRES(5)+1
1977 IF ((NMUONF(1)+NMUONF(2)).EQ.2) NRESF = NRESF+1
1978 NMUONALL = NMUONALL+NMUONF(1)+NMUONF(2)
1979 NGHOSTALL = NGHOSTALL+NGHOST
1980 NTRACKFALL = NTRACKFALL+NTRACKF
1982 NERRALL(I) = NERRALL(I)+NERR(I)
1985 CALL TRACKF_OUT(IR,NTRACKF)
1990 **********************************************************************
1991 SUBROUTINE TRACKF_OUT(IR,NTRACKF)
1992 **********************************************************************
1995 **********************************************************************
1996 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
1998 PARAMETER (MAXHITCH=10000,MAXTRK=50000,MAXHITTOT=20000,
2000 PARAMETER (NBCHAMBER=10,NTRMAX=500)
2002 COMMON/RECOUT/JJO(NBCHAMBER,NTRMAX),XMESO(NBCHAMBER,NTRMAX),
2003 & YMESO(NBCHAMBER,NTRMAX)
2005 COMMON/TRACKFOUT/IEVOUT,NTREVT,JJOUT(NBCHAMBER,NTRMAX),
2006 & ISTAT(NTRMAX),PXZOUT(NTRMAX),TPHIOUT(NTRMAX),
2007 & TALAMOUT(NTRMAX),XVERTOUT(NTRMAX),YVERTOUT(NTRMAX),
2009 & XMESOUT(NBCHAMBER,NTRMAX),YMESOUT(NBCHAMBER,NTRMAX)
2010 & ,PXVOUT(NTRMAX),PYVOUT(NTRMAX),PZVOUT(NTRMAX)
2014 IF (NTREVT.GT.0) THEN
2017 DO IST = 1,NBSTATION
2020 JJOUT(ICH,ITR) = JJO(ICH,ITR)
2021 IF (JJOUT(ICH,ITR).GT.0) THEN
2022 XMESOUT(ICH,ITR) = XMESO(ICH,ITR)
2023 YMESOUT(ICH,ITR) = YMESO(ICH,ITR)
2033 **********************************************************************
2034 SUBROUTINE CHECK_MUON(JJ1,JJ2,JJ3,JJ4,JJ5,
2035 & NTRACKF,NMUONF,NGHOST,NERR,PXZINV,TPHI,TALAM,
2037 **********************************************************************
2038 * Check muon track candidate using GEANT informations
2040 **********************************************************************
2041 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2043 PARAMETER (MAXHITCH=10000,MAXTRK=50000,MAXHITTOT=20000,
2045 PARAMETER (NBCHAMBER=10,NTRMAX=500)
2047 COMMON/CHHIT/XM(NBSTATION,MAXHITCH),YM(NBSTATION,MAXHITCH),
2048 & PHM(NBSTATION,MAXHITCH),ALM(NBSTATION,MAXHITCH),
2049 & IZM(NBSTATION,MAXHITCH),
2050 & IP(NBSTATION,MAXHITCH),JHIT(NBSTATION),
2051 & XMR(NBSTATION,MAXHITCH,2),YMR(NBSTATION,MAXHITCH,2)
2053 COMMON/RHIT/ITYP(MAXHITTOT),XTR(MAXHITTOT),YTR(MAXHITTOT),
2054 & PTOT(MAXHITTOT),ID(MAXHITTOT),IZST(MAXHITTOT),
2055 & PVERT1(MAXHITTOT),PVERT2(MAXHITTOT),PVERT3(MAXHITTOT),
2056 & ZVERT(MAXHITTOT),NHITTOT
2058 COMMON/VERIFGEANT/ITTROUGH(MAXTRK,NBSTATION),
2059 & IT_LIST(MAXTRK),IT_NP(MAXTRK),ITCHECK(MAXTRK),
2062 COMMON/MEASUR/XMES(NBSTATION),YMES(NBSTATION),IZMES(NBSTATION),
2063 & PHMES(NBSTATION),ALMES(NBSTATION),MPOS(NBSTATION),
2066 COMMON/RECOUT/JJO(NBCHAMBER,NTRMAX),XMESO(NBCHAMBER,NTRMAX),
2067 & YMESO(NBCHAMBER,NTRMAX)
2069 COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
2071 COMMON/MAGNET/BL3,B0
2076 DIMENSION NERR(NBSTATION),JJK(NBSTATION),NMUONF(2)
2078 * print*,' *** appel de la subroutine check_muon ***'
2080 IF (JJ1.GT.0.AND.JJ2.GT.0.AND.JJ3.GT.0.AND.JJ4.GT.0
2081 & .AND.JJ5.GT.0) THEN
2082 IDFIND = ID(IP(5,JJ5))
2083 IT = ITRACK(IP(5,JJ5))
2096 IF (XMR(IST,JJK(IST),1).LT.999.) THEN
2097 JJO(ICH1,NTRACKF) = JJK(IST)
2098 XMESO(ICH1,NTRACKF) = XMR(IST,JJK(IST),1)
2099 YMESO(ICH1,NTRACKF) = YMR(IST,JJK(IST),1)
2101 IF (MANG(IST).EQ.1) THEN
2102 JJO(ICH2,NTRACKF) = JJK(IST)
2103 XMESO(ICH2,NTRACKF) = XMR(IST,JJK(IST),2)
2104 YMESO(ICH2,NTRACKF) = YMR(IST,JJK(IST),2)
2108 JJO(ICH2,NTRACKF) = JJK(IST)
2109 XMESO(ICH2,NTRACKF) = XMR(IST,JJK(IST),2)
2110 YMESO(ICH2,NTRACKF) = YMR(IST,JJK(IST),2)
2114 CALL CHFILL(700,SNGL(XVERT),R1,R2)
2115 CALL CHFILL(701,SNGL(YVERT),R1,R2)
2117 IF (ITCHECK(IT).EQ.1) THEN
2120 IF (ID(IP(1,JJ1)).EQ.IDFIND .AND.
2121 & ID(IP(2,JJ2)).EQ.IDFIND .AND.
2122 & ID(IP(3,JJ3)).EQ.IDFIND .AND.
2123 & ID(IP(4,JJ4)).EQ.IDFIND) THEN
2125 IF (ITYP(IP(5,JJ5)).EQ.5) NMUONF(1) = 1 ! Bon muon
2126 IF (ITYP(IP(5,JJ5)).EQ.6) NMUONF(2) = 1 ! Bon muon
2128 PT = SQRT(PVERT1(IP(5,JJ5))**2+PVERT2(IP(5,JJ5))**2)
2129 PZ = PVERT3(IP(5,JJ5))
2130 E = SQRT(PT**2+PZ**2+0.1056**2)
2131 YRAP = 0.5*DLOG((E+PZ)/(E-PZ))
2132 PHIM = DATAN2(PVERT2(IP(5,JJ5)),PVERT1(IP(5,JJ5)))
2133 CALL CHFILL(801,SNGL(YRAP),R1,R2)
2134 CALL CHFILL(901,SNGL(PT),R1,R2)
2135 CALL CHFILL(911,SNGL(PHIM),R1,R2)
2138 NGHOST = NGHOST+1 ! ghost
2140 PT = SQRT(PVERT1(IP(5,JJ5))**2+PVERT2(IP(5,JJ5))**2)
2141 PZ = PVERT3(IP(5,JJ5))
2142 E = SQRT(PT**2+PZ**2+0.1056**2)
2143 YRAP = 0.5*DLOG((E+PZ)/(E-PZ))
2144 PHIM = DATAN2(PVERT2(IP(5,JJ5)),PVERT1(IP(5,JJ5)))
2145 CALL CHFILL(802,SNGL(YRAP),R1,R2)
2146 CALL CHFILL(902,SNGL(PT),R1,R2)
2147 CALL CHFILL(912,SNGL(PHIM),R1,R2)
2151 IF (ITCHECK(ITRACK(IP(5,JJ5))).EQ.1) THEN
2152 IF (JJ3.EQ.0) NERR(3) = NERR(3)+1
2155 IF (JJ1.EQ.0) NERR(1) = NERR(1)+1
2156 IF (JJ2.EQ.0) NERR(2) = NERR(2)+1
2157 IF (JJ1.EQ.0) print*,'hit not found stations 1', NERR(1)
2158 IF (JJ2.EQ.0) print*,'hit not found stations 2', NERR(2)
2161 PT = SQRT(PVERT1(IP(5,JJ5))**2+PVERT2(IP(5,JJ5))**2)
2162 PZ = PVERT3(IP(5,JJ5))
2163 E = SQRT(PT**2+PZ**2+0.1056**2)
2164 YRAP = 0.5*DLOG((E+PZ)/(E-PZ))
2165 PHIM = DATAN2(PVERT2(IP(5,JJ5)),PVERT1(IP(5,JJ5)))
2166 CALL CHFILL(800,SNGL(YRAP),R1,R2)
2167 CALL CHFILL(900,SNGL(PT),R1,R2)
2168 CALL CHFILL(910,SNGL(PHIM),R1,R2)
2175 *************************************************************************
2177 *************************************************************************
2180 *************************************************************************
2181 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2183 PARAMETER (MAXHITCH=10000,MAXTRK=50000,MAXHITTOT=20000,
2185 PARAMETER (NBCHAMBER=10,NTRMAX=500)
2187 COMMON/TRACKFOUT/IEVOUT,NTREVT,JJOUT(NBCHAMBER,NTRMAX),
2188 & ISTAT(NTRMAX),PXZOUT(NTRMAX),TPHIOUT(NTRMAX),
2189 & TALAMOUT(NTRMAX),XVERTOUT(NTRMAX),YVERTOUT(NTRMAX),
2191 & XMESOUT(NBCHAMBER,NTRMAX),YMESOUT(NBCHAMBER,NTRMAX)
2192 & ,PXVOUT(NTRMAX),PYVOUT(NTRMAX),PZVOUT(NTRMAX)
2194 COMMON/CHHIT/XM(NBSTATION,MAXHITCH),YM(NBSTATION,MAXHITCH),
2195 & PHM(NBSTATION,MAXHITCH),ALM(NBSTATION,MAXHITCH),
2196 & IZM(NBSTATION,MAXHITCH),
2197 & IP(NBSTATION,MAXHITCH),JHIT(NBSTATION),
2198 & XMR(NBSTATION,MAXHITCH,2),YMR(NBSTATION,MAXHITCH,2)
2200 COMMON/RHIT/ITYP(MAXHITTOT),XTR(MAXHITTOT),YTR(MAXHITTOT),
2201 & PTOT(MAXHITTOT),ID(MAXHITTOT),IZST(MAXHITTOT),
2202 & PVERT1(MAXHITTOT),PVERT2(MAXHITTOT),PVERT3(MAXHITTOT),
2203 & ZVERT(MAXHITTOT),NHITTOT
2205 COMMON/RECOUT/JJO(NBCHAMBER,NTRMAX),XMESO(NBCHAMBER,NTRMAX),
2206 & YMESO(NBCHAMBER,NTRMAX)
2208 COMMON/TRACKSUM/NRES(5),NRESF,NTRMUALL,NMUONALL,NGHOSTALL,
2209 & NTRACKFALL,NERRALL(NBSTATION),IR
2211 COMMON/PRECSUM/NRESF1,NMUONALL1,NGHOSTALL1,NTRACKFALL1
2213 REAL*4 PXR,PYR,PZR,ZVR,CHI2R,PXV,PYV,PZV
2214 COMMON/PAWCR4/IEVR,NTRACKR,ISTATR(NTRMAX),ISIGNR(NTRMAX),
2215 & PXR(NTRMAX),PYR(NTRMAX),PZR(NTRMAX),ZVR(NTRMAX),
2216 & CHI2R(NTRMAX),PXV(NTRMAX),PYV(NTRMAX),PZV(NTRMAX)
2218 COMMON/DEBEVT/IDEBUG
2220 DIMENSION ISEL(NTRMAX)
2223 CALL RECO_SELECT(ISEL)
2229 IF (ISEL(ITR).EQ.1) THEN
2230 NTRACKFALL1 = NTRACKFALL1 + 1
2231 IF (ISTAT(ITR).EQ.1) THEN
2233 NMUONALL1 = NMUONALL1 + 1
2234 ELSEIF (ISTAT(ITR).EQ.2) THEN
2236 NGHOSTALL1 = NGHOSTALL1 + 1
2241 IF (NMUF.GE.2) NRESF1 = NRESF1 + 1
2245 IF (ISEL(ITR).EQ.1) THEN
2246 NTRACKR = NTRACKR + 1
2247 ISTATR(NTRACKR) = ISTAT(ITR)
2249 IF (PXZOUT(ITR).LT.0.) ISIGNR(NTRACKR) = -1
2250 PXZ = ABS(PXZOUT(ITR))
2251 PHI = ATAN(TPHIOUT(ITR))
2252 ALAM = ATAN(TALAMOUT(ITR))
2253 PYR(NTRACKR) = PXZ*SIN(PHI)
2254 PXR(NTRACKR) = PXZ*TAN(ALAM)
2255 PZR(NTRACKR) = PXZ*COS(PHI)
2257 CHI2R(NTRACKR) = CHI2OUT(ITR)
2258 PXV(NTRACKR) = PXVOUT(ITR)
2259 PYV(NTRACKR) = PYVOUT(ITR)
2260 PZV(NTRACKR) = PZVOUT(ITR)
2264 CALL CHFNT(IEVR,NTRACKR,ISTATR,ISIGNR,
2265 & PXR,PYR,PZR,ZVR,CHI2R,PXV,PYV,PZV)
2267 IF (IDEBUG.GE.2) THEN
2268 PRINT *,'RECO_SUM evt number :',IEVOUT
2269 PRINT *,'RECO_SUM nb of track /evt :',NTRACKR
2270 PRINT *,'RECO_SUM nb of good muon /evt :',NMUF
2271 PRINT *,'RECO_SUM nb of ghost /evt :',NGHF
2272 IF (NTRACKR.GT.0) THEN
2274 PRINT *,'RECO_SUM track number :',ITR
2275 PRINT *,'RECO_SUM CHI2OUT :',CHI2R(ITR)
2276 PRINT *,' PX GEANT= ', PXV(ITR),' PX RECONS= ',PYR(ITR)
2277 PRINT *,' PY GEANT= ', PYV(ITR),' PY RECONS= ',PXR(ITR)
2278 PRINT *,' PZ GEANT= ', PZV(ITR),' PZ RECONS= ',PZR(ITR)
2279 PXZV = SQRT( PYV(ITR)**2+PZV(ITR)**2)
2280 PXZR = SQRT( PXR(ITR)**2+PZR(ITR)**2)
2281 PRINT *,' PXZ GEANT= ', PXZV,' PXZ RECONS= ',PXZR
2282 FIV=ATAN2(DBLE(PYV(ITR)),DBLE(PZV(ITR)))
2283 ALAMV=ATAN2(DBLE(PXV(ITR)),DBLE(PXZV))
2284 FIR=ATAN2(DBLE(PXR(ITR)),DBLE(PZR(ITR)))
2285 ALAMR=ATAN2(DBLE(PYR(ITR)),DBLE(PXZR))
2286 ** PRINT *,' PHI GEANT= ',FIV,' PXZ RECONS= ',FIR
2287 ** PRINT *,' ALAM GEANT= ',ALAMV,' ALAM RECONS= ',ALAMR
2295 *************************************************************************
2296 SUBROUTINE RECO_TERM
2297 *************************************************************************
2300 *************************************************************************
2301 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2303 PARAMETER(NBSTATION=5)
2305 COMMON/TRACKSUM/NRES(5),NRESF,NTRMUALL,NMUONALL,NGHOSTALL,
2306 & NTRACKFALL,NERRALL(NBSTATION),IR
2308 COMMON/PRECSUM/NRESF1,NMUONALL1,NGHOSTALL1,NTRACKFALL1
2310 COMMON/DEBEVT/IDEBUG
2312 CHARACTER*50 FILEBKG,FILERES,FILEOUT,FILEMIN
2314 IF (IDEBUG.GE.1) THEN
2316 PRINT *,'RECO_TERM ***** SUMMARY TRACK-FINDING *****'
2317 PRINT *,'RECO_TERM nb of resonances :',NRES(5)
2318 PRINT *,'RECO_TERM nb of resonances 45 :',NRES(1)
2319 PRINT *,'RECO_TERM nb of resonances 345 :',NRES(2)
2320 PRINT *,'RECO_TERM nb of resonances found :',NRESF
2321 PRINT *,'RECO_TERM nb of muon track :',NTRMUALL
2322 PRINT *,'RECO_TERM nb of track found :',NTRACKFALL
2323 PRINT *,'RECO_TERM nb of muon track found :',NMUONALL
2324 PRINT *,'RECO_TERM nb of ghost found :',NGHOSTALL
2326 PRINT *,'RECO_TERM nb of error in st',I,':',NERRALL(I)
2330 PRINT *,'RECO_TERM ***** SUMMARY PRECISION *****'
2331 PRINT *,'RECO_TERM nb of resonances found :',NRESF1
2332 PRINT *,'RECO_TERM nb of track found :',NTRACKFALL1
2333 PRINT *,'RECO_TERM nb of muon track found :',NMUONALL1
2334 PRINT *,'RECO_TERM nb of ghost found :',NGHOSTALL1
2342 *************************************************************************
2343 SUBROUTINE RECO_PRECISION
2344 *************************************************************************
2347 *************************************************************************
2348 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2350 PARAMETER (MAXHITCH=10000,MAXTRK=50000,MAXHITTOT=20000,
2351 & NBSTATION=5,MAXCAN=1000,NTRMAX=500)
2352 PARAMETER (NPLANE=10)
2354 COMMON/PARAM/ZPLANEP(NPLANE),THICK,XPREC,YPREC,B0,BL3,ZMAGS,
2355 & ZMAGE,ZABS,XMAG,ZBP1,ZBP2,CONST
2357 COMMON/CHHIT/XM(NBSTATION,MAXHITCH),YM(NBSTATION,MAXHITCH),
2358 & PHM(NBSTATION,MAXHITCH),ALM(NBSTATION,MAXHITCH),
2359 & IZM(NBSTATION,MAXHITCH),
2360 & IP(NBSTATION,MAXHITCH),JHIT(NBSTATION),
2361 & XMR(NBSTATION,MAXHITCH,2),YMR(NBSTATION,MAXHITCH,2)
2363 COMMON/RHIT/ITYP(MAXHITTOT),XTR(MAXHITTOT),YTR(MAXHITTOT),
2364 & PTOT(MAXHITTOT),ID(MAXHITTOT),IZST(MAXHITTOT),
2365 & PVERT1(MAXHITTOT),PVERT2(MAXHITTOT),PVERT3(MAXHITTOT),
2366 & ZVERT(MAXHITTOT),NHITTOT
2368 COMMON/TRACKFOUT/IEVOUT,NTREVT,JJOUT(NPLANE,NTRMAX),
2369 & ISTAT(NTRMAX),PXZOUT(NTRMAX),TPHIOUT(NTRMAX),
2370 & TALAMOUT(NTRMAX),XVERTOUT(NTRMAX),YVERTOUT(NTRMAX),
2372 & XMESOUT(NPLANE,NTRMAX),YMESOUT(NPLANE,NTRMAX)
2373 & ,PXVOUT(NTRMAX),PYVOUT(NTRMAX),PZVOUT(NTRMAX)
2376 COMMON/MEAS/LPLANE(NPLANE),XMP(NPLANE),YMP(NPLANE)
2378 COMMON/FCNOUT/PXZEA,ALAMEA,PHIEA,XEA,YEA,NPLU,CHI2
2380 COMMON/PRECCUT/PCUT,PTCUT,CHI2CUT
2382 DIMENSION PARMU(MAXCAN,NPLANE,2),LPLANEMU(MAXCAN,NPLANE)
2384 IF (NTREVT.EQ.0) RETURN
2388 DO IST = 1,NBSTATION
2391 ** print *,' ich=',ich
2392 IF (JJOUT(ICH,ITR).GT.0) THEN
2393 LPLANEMU(ITR,ICH) = 1
2394 PARMU(ITR,ICH,1) = XMESOUT(ICH,ITR)
2395 PARMU(ITR,ICH,2) = YMESOUT(ICH,ITR)
2396 ** print *,' x,y ', PARMU(ITR,ICH,1),PARMU(ITR,ICH,2)
2398 LPLANEMU(ITR,ICH) = 0
2407 LPLANE(ICH) = LPLANEMU(ICAN,ICH)
2408 XMP(ICH) = PARMU(ICAN,ICH,1)
2409 YMP(ICH) = PARMU(ICAN,ICH,2)
2412 IF (LPLANE(1).GT.0) THEN
2421 IF (LPLANE(3).GT.0) THEN
2430 IF (LPLANE(7).GT.0) THEN
2437 IF (LPLANE(9).GT.0) THEN
2445 PHIAV = DATAN2((X2-X1),(ZPLANEP(IPL2)-ZPLANEP(IPL1)))
2446 PHIAP = DATAN2((X4-X3),(ZPLANEP(IPL4)-ZPLANEP(IPL3)))
2447 DPHI = (PHIAP-PHIAV)
2449 IF (DPHI.LT.0.) ASIGN = -1. ! CCC
2450 PXZ = CONST/DABS(DPHI)
2452 IF (PXZ.LT.PCUT) GO TO 2
2454 PXZINVI = ASIGN/PXZ ! CCC
2455 ** PXZINVI = 1./PXZOUT(ICAN) ! CCC
2456 ** PXZINVI = -1./49. ! CCC
2458 ALAMI = DATAN2((Y2-Y1),DSQRT((X2-X1)**2
2459 & +(ZPLANEP(IPL2)-ZPLANEP(IPL1))**2))
2462 ** print *,' avant prec_fit pxzi phii alami x y',1./ PXZINVI,
2463 ** & PHII, ALAMI ,XVR,YVR
2464 ** PRINT *,' X1 X2 X3 X4',X1,X2,X3,X4
2465 ** PRINT *,' Z1 Z2 Z3 Z4',ZPLANEP(IPL1),ZPLANEP(IPL2),
2466 ** & ZPLANEP(IPL3),ZPLANEP(IPL4)
2467 ** PRINT *,' CONST= ',CONST
2470 IF (CHI2OUT(ICAN).GT.CHI2CUT) GO TO 2
2472 ** Fit des traces apres l'absorbeur
2473 CALL PREC_FIT (PXZINVI,PHII,ALAMI,XVR,YVR,
2474 & PXZINVF,PHIF,ALAMF,XVERTF,YVERTF,EPXZINV,EPHI,EALAM,
2477 ** Correction de Branson
2478 CALL BRANSON(PXZEA,PHIEA,ALAMEA,XEA,YEA)
2481 PX1 = PXZ1*DSIN(PHIEA)
2482 PY1 = PXZ1*DTAN(ALAMEA)
2483 PT1 = DSQRT(PX1**2+PY1**2)
2485 IF (PT1.LT.PTCUT) GO TO 2
2487 IF ((CHI2/FLOAT(2*NPLU-5)).GT.CHI2CUT) GO TO 2
2491 JJOUT(ICH,NTRACK) = JJOUT(ICH,ICAN)
2492 XMESOUT(ICH,NTRACK) = XMESOUT(ICH,ICAN)
2493 YMESOUT(ICH,NTRACK) = YMESOUT(ICH,ICAN)
2495 ISTAT(NTRACK) = ISTAT(ICAN)
2496 PXZOUT(NTRACK) = PXZEA
2497 TPHIOUT(NTRACK) = DTAN(PHIEA)
2498 TALAMOUT(NTRACK) = DTAN(ALAMEA)
2499 XVERTOUT(NTRACK) = XEA
2500 YVERTOUT(NTRACK) = YEA
2501 CHI2OUT(NTRACK) = CHI2/FLOAT(2*NPLU-5)
2504 ** print *,' reco_precision pxz tphi talam xvert yvert chi2',
2505 ** & PXZEA,PHIEA,ALAMEA,
2506 ** & XEA,YEA,CHI2/FLOAT(2*NPLU-5)
2515 ************************************************************************
2516 SUBROUTINE ORDONNE_HIT(ICH,HCUT)
2517 **************************************************************************
2519 * Sort hits in station ICH according to the "impact parameter" HHIT
2521 **************************************************************************
2523 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2525 PARAMETER (MAXHITCH=10000,MAXTRK=50000,MAXHITTOT=20000,
2528 COMMON/CHHIT/XM(NBSTATION,MAXHITCH),YM(NBSTATION,MAXHITCH),
2529 & PHM(NBSTATION,MAXHITCH),ALM(NBSTATION,MAXHITCH),
2530 & IZM(NBSTATION,MAXHITCH),
2531 & IP(NBSTATION,MAXHITCH),JHIT(NBSTATION),
2532 & XMR(NBSTATION,MAXHITCH,2),YMR(NBSTATION,MAXHITCH,2)
2534 COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
2536 COMMON/VERIFGEANT/ITTROUGH(MAXTRK,NBSTATION),
2537 & IT_LIST(MAXTRK),IT_NP(MAXTRK),ITCHECK(MAXTRK),
2540 COMMON/HCHHIT/HHIT(MAXHITCH),INDEXTAB(MAXHITCH),INDEXMAX
2542 COMMON/DEBEVT/IDEBUG
2545 * tri des impulsion par ordre decroissant
2546 * le tab INDEXTAB contient les j ordonnes
2547 * INDEXMAX est l indice max du tableau = NBHIT si pas de contrainte
2551 * boucle sur le nombre de hits candidats de la station
2555 IF (PHM(ICH,J).LT.6.3) THEN !2pi=6.3 radian
2557 * calcul du h dans XOY a z=0
2558 HHIT(J)=ABS(XM(ICH,J)-ZPLANE(ICH)*PHM(ICH,J))
2560 IF (HHIT(J).LT.HCUT) THEN
2562 INDEXTAB(INDEXMAX)=J
2563 ELSEIF(ITCHECK(ITRACK(IP(ICH,J))).EQ.1) THEN
2564 IF (IDEBUG.GE.2) THEN
2565 PRINT *,'ORDONNE_HIT rejet muon/res in st.',ICH,
2575 H4(I) = SNGL(HHIT(I))
2578 CALL SORTZV(H4,INDEXTAB,INDEXMAX,-1,0,1)
2581 HHIT(I) = DBLE(H4(I))
2583 ** PRINT *,'ORDONNE st. numero',ICH
2584 ** PRINT *,'ORDONNE nb de hits initiaux dans st.',ICH,':',JHIT(ICH)
2585 ** PRINT *,'ORDONNE nb de hits avec mes. angulaire:',JJ
2586 ** PRINT *,'ORDONNE nb de hits avec mes. ang. et cut en Pxz:',INDEXMAX
2587 IF (IDEBUG.GE.2) THEN
2588 PRINT *,'ORDONNE_HIT nb de hits accepte dans st.',ICH,':',
2594 ***********************************************************************************
2595 SUBROUTINE DISTMIN4(X1,Y1,PHI1,ALAM1,ICH,EX,EY,EPHI,ELAM,IFIND,
2597 ***********************************************************************************
2598 * Find the nearest hit in station ICH in the (X,Y,lambda,phi) phase space
2600 ***********************************************************************************
2602 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2604 PARAMETER (MAXHITCH=10000,NBSTATION=5)
2606 COMMON/CHHIT/XM(NBSTATION,MAXHITCH),YM(NBSTATION,MAXHITCH),
2607 & PHM(NBSTATION,MAXHITCH),ALM(NBSTATION,MAXHITCH),
2608 & IZM(NBSTATION,MAXHITCH),
2609 & IP(NBSTATION,MAXHITCH),JHIT(NBSTATION),
2610 & XMR(NBSTATION,MAXHITCH,2),YMR(NBSTATION,MAXHITCH,2)
2611 DIMENSION IFIND2(10)
2621 IF (PHM(ICH,I).LE.6.3) THEN ! vector measurement
2622 IF (ABS(PHI1-PHM(ICH,I)) .LT. EPHI .AND.
2623 & ABS(ALAM1-ALM(ICH,I)) .LT. ELAM .AND.
2624 & ABS(X1-XM(ICH,I)) .LT. EX .AND.
2625 & ABS(Y1-YM(ICH,I)) .LT. EY) THEN
2626 DIST = ((PHI1-PHM(ICH,I))/EPHI)**2 +
2627 & ((ALAM1-ALM(ICH,I))/ELAM)**2 +
2628 & ((X1-XM(ICH,I))/EX)**2 +
2629 & ((Y1-YM(ICH,I))/EY)**2
2631 IF (NF.LE.10) IFIND2(NF) = I
2632 IF (DIST .LT. DISTMIN) THEN
2642 ***********************************************************************************
2643 SUBROUTINE DISTMIN2(X1,Y1,X2,Y2,ICH,EX1,EY1,EX2,EY2,IFIND,IFIND2)
2644 ***********************************************************************************
2645 * Find the nearest hit in station ICH in the (X,Y) space
2647 ***********************************************************************************
2648 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2650 PARAMETER (MAXHITCH=10000,NBSTATION=5)
2652 COMMON/CHHIT/XM(NBSTATION,MAXHITCH),YM(NBSTATION,MAXHITCH),
2653 & PHM(NBSTATION,MAXHITCH),ALM(NBSTATION,MAXHITCH),
2654 & IZM(NBSTATION,MAXHITCH),
2655 & IP(NBSTATION,MAXHITCH),JHIT(NBSTATION),
2656 & XMR(NBSTATION,MAXHITCH,2),YMR(NBSTATION,MAXHITCH,2)
2657 DIMENSION IFIND2(10)
2667 IF (IZM(ICH,I).EQ.1) THEN ! 1st chamber
2676 IF (ICH.EQ.4.OR.ICH.EQ.5) THEN
2677 IF (IZM(ICH,I).EQ.1) THEN
2685 IF (ABS(X-XM(ICH,I)) .LT. EX .AND.
2686 & ABS(Y-YM(ICH,I)) .LT. EY) THEN
2687 DIST = ((X-XM(ICH,I))/EX)**2 +
2688 & ((Y-YM(ICH,I))/EY)**2
2690 IF (NF.LE.10) IFIND2(NF) = I
2691 IF (DIST .LT. DISTMIN) THEN
2700 ********************************************************************************
2701 SUBROUTINE H_ACCEPTANCE(ICH)
2702 ********************************************************************************
2703 * Etude de l'acceptance des resonnances en fonction du H
2704 * dans la station ICH
2708 ********************************************************************************
2710 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2712 PARAMETER (MAXHITCH=10000,MAXTRK=50000,MAXHITTOT=20000,
2715 COMMON/CHHIT/XM(NBSTATION,MAXHITCH),YM(NBSTATION,MAXHITCH),
2716 & PHM(NBSTATION,MAXHITCH),ALM(NBSTATION,MAXHITCH),
2717 & IZM(NBSTATION,MAXHITCH),
2718 & IP(NBSTATION,MAXHITCH),JHIT(NBSTATION),
2719 & XMR(NBSTATION,MAXHITCH,2),YMR(NBSTATION,MAXHITCH,2)
2721 COMMON/RHIT/ITYP(MAXHITTOT),XTR(MAXHITTOT),YTR(MAXHITTOT),
2722 & PTOT(MAXHITTOT),ID(MAXHITTOT),IZST(MAXHITTOT),
2723 & PVERT1(MAXHITTOT),PVERT2(MAXHITTOT),PVERT3(MAXHITTOT),
2724 & ZVERT(MAXHITTOT),NHITTOT
2727 COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
2729 COMMON/VERIFGEANT/ITTROUGH(MAXTRK,NBSTATION),
2730 & IT_LIST(MAXTRK),IT_NP(MAXTRK),ITCHECK(MAXTRK),
2733 COMMON/HCHHIT/HHIT(MAXHITCH),INDEXTAB(MAXHITCH),INDEXMAX
2740 IF (ITYP(IP(ICH,J)).EQ.5.OR.ITYP(IP(ICH,J)).EQ.6) THEN
2741 ISTAK = ID(IP(ICH,J))
2742 ISTAK = MOD(ISTAK,30000)
2743 ISTAK = MOD(ISTAK,10000)
2744 IF (ISTAK.EQ.0) THEN
2745 ** PRINT *,'ACCEPT. id du muon dans st.',ICH,':',ITYP(IP(ICH,J))
2750 * PRINT *,'ACCEPT. nb de muons/res total dans st.',ICH,':',NMUONI
2755 * Sort hits in st. z
2756 CALL ORDONNE_HIT(ICH,HCUT)
2759 IIND = IP(ICH,INDEXTAB(IND))
2762 ISTAK = MOD(ISTAK,30000)
2763 ISTAK = MOD(ISTAK,10000)
2764 ** PRINT *,' IDPART=',IDPART,' ISTAK=',ISTAK
2765 IF (IDPART.EQ.5.OR.IDPART.EQ.6.AND.ISTAK.EQ.0) THEN
2769 IF (NMUON.EQ.2.AND.NMUONI.EQ.2) THEN
2770 CALL CHFILL(ICH*100,SNGL(HCUT),R1,R2)
2777 ********************************************************************************
2778 SUBROUTINE OLDFOLLOW(ZSTR,PEST,IFLAG,XPL,YPL,PHPL,ALPL)
2779 ********************************************************************************
2780 * Calculate the particle trajectory in the spectrometer and
2781 * (XPL,YPL,PHPL,ALPL)
2782 * for the 5 stations.
2784 ********************************************************************************
2785 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2787 PARAMETER(NBSTATION=5)
2789 DIMENSION XPL(NBSTATION,2),YPL(NBSTATION,2),PHPL(NBSTATION),
2790 & ALPL(NBSTATION),PEST(NBSTATION)
2792 COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
2794 COMMON /MEASUR/XMES(NBSTATION),YMES(NBSTATION),IZMES(NBSTATION),
2795 & PHMES(NBSTATION),ALMES(NBSTATION),MPOS(NBSTATION),
2797 COMMON /MAGNET/BL3,B0
2799 LOGICAL LFLAG, LFLAG1
2809 PX = -ABS(PXZ)*SIN(PHI)
2810 PZ = -ABS(PXZ)*COS(PHI)
2811 PXY = SQRT(PX**2 + PY**2)
2812 FI=ATAN2(DBLE(PY),DBLE(PX))
2816 RS = PXY*(100.0/(0.299792458*BL3))
2817 IF(PXZINV.LT.0.0) RS = -RS
2818 * XC = XSTR + RS*SIN(FI)
2819 * YC = YSTR - RS*COS(FI)
2824 * PRINT *, XC,YC,RS,FI,TTHET,PXY,PZ
2828 IF (IFLAG.EQ.3 .OR. MPOS(J).EQ.1) THEN
2829 IF(ZPLANE(J) .GT. ZCOIL) THEN
2830 * DFI = (ZPLANE(J)-ZSTR)/(TTHET*RS)
2832 * XPL(J,1) = XC - RS*SIN(FIN)
2833 * YPL(J,1) = YC + RS*COS(FIN)
2834 DFR = (ZPLANE(J)-ZSTR)/TTHET
2835 XPL(J,1) = XSTR + DFR*COSFI + 0.5D0*DFR*DFR*SINFI/RS
2836 YPL(J,1) = YSTR + DFR*SINFI - 0.5D0*DFR*DFR*COSFI/RS
2837 DFR2 = (ZPLANE(J)-DZ_PL(J)-ZSTR)/TTHET
2838 XPL(J,2) = XSTR + DFR2*COSFI + 0.5D0*DFR2*DFR2*SINFI/RS
2839 YPL(J,2) = YSTR + DFR2*SINFI - 0.5D0*DFR2*DFR2*COSFI/RS
2840 IF (IFLAG.EQ.3 .OR. MANG(J).EQ.1) THEN
2843 PX = PX0 + DFR * (PY0 - 0.5D0*PX0*DFR/RS) / RS
2844 PY = PY0 - DFR * (PX0 + 0.5D0*PY0*DFR/RS) / RS
2846 ALPL(J)=ATAN(PY/SQRT(PX**2+PZ**2))
2850 * DFI = (ZCOIL-ZSTR)/(TTHET*RS)
2852 * XCOIL = XC - RS*SIN(FIN)
2853 * YCOIL = YC + RS*COS(FIN)
2854 DFR = (ZCOIL-ZSTR)/TTHET
2855 XCOIL = XSTR + DFR*COSFI + 0.5D0*DFR*DFR*SINFI/RS
2856 YCOIL = YSTR + DFR*SINFI - 0.5D0*DFR*DFR*COSFI/RS
2859 PX = PX0 + DFR * (PY0 - 0.5D0*PX0*DFR/RS) / RS
2860 PY = PY0 - DFR * (PX0 + 0.5D0*PY0*DFR/RS) / RS
2861 PXZ = SQRT(PX**2 + PZ**2)
2865 RD = PXZ*(100.0/(0.299792458*B0))
2866 IF(PXZINV.LT.0.0) RD = -RD
2867 ZC = ZCOIL - RD*SIN(PHI)
2868 XC = XCOIL + RD*COS(PHI)
2869 IF(ABS(ZMAGEND-ZC).GT.ABS(RD)) STOP 'FOLLOW'
2872 IF(ZPLANE(J) .GT. ZMAGEND) THEN
2873 FIN = ASIN((ZPLANE(J) - ZC)/RD)
2874 XPL(J,1)= XC - RD*COS(FIN)
2875 YPL(J,1)= YCOIL - RD*(FIN - PHI)*TALAM
2876 FIN2 = ASIN((ZPLANE(J)-DZ_PL(J) - ZC)/RD)
2877 XPL(J,2)= XC - RD*COS(FIN2)
2878 YPL(J,2)= YCOIL - RD*(FIN2 - PHI)*TALAM
2883 FIN = ASIN((ZMAGEND - ZC)/RD)
2884 XMAGEND = XC - RD*COS(FIN)
2885 YMAGEND = YCOIL - RD*(FIN - PHI)*TALAM
2890 XPL(J,1) = XMAGEND + (ZPLANE(J)-ZMAGEND)*TPHI
2891 YPL(J,1) = YMAGEND - (ZPLANE(J)-ZMAGEND)*TALAM/CPHI
2892 XPL(J,2) = XMAGEND + (ZPLANE(J)-DZ_PL(J)-ZMAGEND)*TPHI
2893 YPL(J,2) = YMAGEND - (ZPLANE(J)-DZ_PL(J)-ZMAGEND)*
2903 ********************************************************************************
2904 SUBROUTINE FOLLOW(ZSTR,PEST,IFLAG,XPL,YPL,PHPL,ALPL)
2905 ********************************************************************************
2906 * Calculate the particle trajectory in the spectrometer
2907 * (XPL,YPL,PHPL,ALPL)
2908 * for the 5 stations.
2910 ********************************************************************************
2911 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2913 PARAMETER(NBSTATION=5,NPLANE=10)
2915 DIMENSION XPL(NBSTATION,2),YPL(NBSTATION,2),PHPL(NBSTATION),
2916 & ALPL(NBSTATION),PEST(NBSTATION)
2918 COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
2920 COMMON/PARAM/ZPLANEP(NPLANE),THICK,XPREC,YPREC,B0,BL3,ZMAGS,
2921 & ZMAGE,ZABS,XMAG,ZBP1,ZBP2,CONST
2923 COMMON /MEASUR/XMES(NBSTATION),YMES(NBSTATION),IZMES(NBSTATION),
2924 & PHMES(NBSTATION),ALMES(NBSTATION),MPOS(NBSTATION),
2928 DIMENSION VECT(7),VOUT(7)
2934 IF (PEST(1).LT.0.) ASIGN = -1.
2939 PXZ = DABS(1./PEST(1))
2944 PTOT = PXZ/DCOS(ALAM)
2958 ** PRINT *,' AV GRKUTA ASIGN',ASIGN,' THET',THET
2962 IST = INT(FLOAT(ICH+1)/2.)
2965 DO WHILE (Z.GE.0..AND.Z.LT.ABS(ZPLANEP(ICH))
2966 & .AND.NSTEP.LE.NSTEPMAX)
2968 ** WRITE(6,*) NSTEP,(VECT(I),I=1,7)
2969 ** CALL RECO_GRKUTA (ASIGN,STEP,VECT,VOUT) ! CCC
2970 CALL RECO_GHELIX (ASIGN,STEP,VECT,VOUT)
2976 IF (IST.NE.ISTOLD) THEN
2981 XPL(IST,IPCH) = VECT(1)-(Z-ABS(ZPLANEP(ICH)))*VECT(4)/VECT(6)
2982 YPL(IST,IPCH) = VECT(2)-(Z-ABS(ZPLANEP(ICH)))*VECT(5)/VECT(6)
2984 DX = XPL(IST,2)-XPL(IST,1)
2985 DY = YPL(IST,2)-YPL(IST,1)
2986 PHPL(IST) = -1.*DATAN2(DX,DZ_PL(IST))
2987 ALPL(IST) = DATAN2(DY,DSQRT(DX**2+DZ_PL(IST)**2))
2991 ** print *,' vect= ',vect(1),vect(2),vect(3),vect(4),vect(5),
2992 ** & vect(6),vect(7)
2997 *******************************************************************************
2998 SUBROUTINE FCN(NPAR,GRAD,FVAL,PEST,IFLAG,FUTIL)
2999 *******************************************************************************
3000 * Calculate FVAL=CHI2 the function minimized by minuit for a given track
3002 *******************************************************************************
3003 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
3005 PARAMETER(NBSTATION=5)
3007 * DIMENSION PEST(*),GRAD(*)
3008 DIMENSION PEST(5),GRAD(5)
3009 DIMENSION PEEST(NBSTATION)
3011 COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
3014 COMMON/PRECIS/EEXM(NBSTATION),EEYM(NBSTATION),EEPH(NBSTATION),
3017 COMMON /MEASUR/XMES(NBSTATION),YMES(NBSTATION),IZMES(NBSTATION),
3018 & PHMES(NBSTATION),ALMES(NBSTATION), MPOS(NBSTATION),
3021 COMMON /PLANE/XPL(NBSTATION,2),YPL(NBSTATION,2),PHPL(NBSTATION),
3022 & ALPL(NBSTATION),CHI2PL
3024 COMMON/VERTEX/ERRV,IVERTEX
3026 DIMENSION XC(NBSTATION),YC(NBSTATION)
3033 IF(IVERTEX.EQ.1) THEN
3034 PEEST(4)=PEST(4) ! position du vertex
3041 ALAM = DATAN(PEST(3))
3042 PXZ = DABS(1./PEST(1))
3043 PTOT = PXZ/DCOS(ALAM)
3046 CALL FOLLOW (0.0D0,PEEST,IFLAG,XPL,YPL,PHPL,ALPL) ! calcul des
3047 IF(IFLAG.EQ.1) THEN ! points d impacts dans les
3048 PRINT *,'FCN ',XPL(4,1),XMES(4) ! plans
3049 PRINT *,'FCN ',YPL(4,1),YMES(4)
3050 PRINT *,'FCN ',XPL(5,1),XMES(5)
3051 PRINT *,'FCN ',YPL(5,1),YMES(5)
3059 IF (IVERTEX.EQ.1) THEN
3061 FVAL = RECOCHI2(MPOS,MANG,XMES,YMES,ALMES,PHMES,
3062 & XC,YC,ALPL,PHPL,PTOT,IZMES,NPLPL)
3070 IF (MPOS(J).EQ.1) THEN
3072 XPLC = XPL(J,IZMES(J))
3073 YPLC = YPL(J,IZMES(J))
3074 FF = (XMES(J) - XPLC)/EEXM(J)
3076 FF = (YMES(J) - YPLC)/EEYM(J)
3079 IF (MANG(J).EQ.1) THEN
3081 FF = (PHMES(J) - PHPL(J))/EEPH(J)
3083 FF = (ALMES(J) - ALPL(J))/EEAL(J)
3091 IF (IVERTEX.EQ.1) NPARAM = 5
3092 CHI2PL = FVAL/FLOAT(2*NPLPL-NPARAM)
3096 ********************************************************************************
3097 SUBROUTINE STOCK_CANDIDAT(ICH1,JHITCH1,ICH2,IFIND,IFIND2,EXM,EYM,
3098 & EPH,EAL,NCAN,ICODE)
3099 ********************************************************************************
3100 * Fill common CANDIDAT with track candidates from the search in stations 4&5
3102 ********************************************************************************
3104 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
3106 PARAMETER (MAXHITTOT=20000,MAXHITCH=10000,NBSTATION=5,MAXCAN=1000)
3108 COMMON/CHHIT/XM(NBSTATION,MAXHITCH),YM(NBSTATION,MAXHITCH),
3109 & PHM(NBSTATION,MAXHITCH),ALM(NBSTATION,MAXHITCH),
3110 & IZM(NBSTATION,MAXHITCH),
3111 & IP(NBSTATION,MAXHITCH),JHIT(NBSTATION),
3112 & XMR(NBSTATION,MAXHITCH,2),YMR(NBSTATION,MAXHITCH,2)
3114 COMMON/RHIT/ITYP(MAXHITTOT),XTR(MAXHITTOT),YTR(MAXHITTOT),
3115 & PTOT(MAXHITTOT),ID(MAXHITTOT),IZST(MAXHITTOT),
3116 & PVERT1(MAXHITTOT),PVERT2(MAXHITTOT),PVERT3(MAXHITTOT),
3117 & ZVERT(MAXHITTOT),NHITTOT
3120 COMMON/CANDIDAT/JCAN(NBSTATION,MAXCAN),JCANTYP(MAXCAN),
3121 & EEX(MAXCAN),EEY(MAXCAN),EEP(MAXCAN),EEA(MAXCAN)
3122 DIMENSION IFIND2(10)
3124 ** PRINT *,'STOCK st. init.=',ICH1,'id. init.=',ID(IP(ICH1,JHITCH1))
3125 ** PRINT *,'STOCK st. finale=',ICH2,'id. final=',ID(IP(ICH2,IFIND))
3126 ** PRINT *,'STOCK ifind=',IFIND
3127 ** PRINT *,'STOCK icode=',ICODE
3130 IF (IFIND2(I).GT.0) THEN
3132 JCAN(ICH1,NCAN) = JHITCH1
3133 JCAN(ICH2,NCAN) = IFIND2(I)
3134 JCANTYP(NCAN) = ICODE
3144 *******************************************************************************
3145 SUBROUTINE FCNOLD(NPAR,GRAD,FVAL,PEST,IFLAG,FUTIL)
3146 *******************************************************************************
3147 * Calculate FVAL=CHI2 the function minimized by minuit for a given track
3149 *******************************************************************************
3150 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
3152 PARAMETER(NBSTATION=5)
3154 * DIMENSION PEST(*),GRAD(*)
3155 DIMENSION PEST(5),GRAD(5)
3156 DIMENSION PEEST(NBSTATION)
3158 COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
3161 COMMON/PRECIS/EEXM(NBSTATION),EEYM(NBSTATION),EEPH(NBSTATION),
3164 COMMON /MEASUR/XMES(NBSTATION),YMES(NBSTATION),IZMES(NBSTATION),
3165 & PHMES(NBSTATION),ALMES(NBSTATION), MPOS(NBSTATION),
3168 COMMON /PLANE/XPL(NBSTATION,2),YPL(NBSTATION,2),PHPL(NBSTATION),
3169 & ALPL(NBSTATION),CHI2PL
3171 COMMON/VERTEX/ERRV,IVERTEX
3176 IF(IVERTEX.EQ.1) THEN
3177 PEEST(4)=PEST(4) ! position du vertex
3184 CALL FOLLOW (0.0D0,PEEST,IFLAG,XPL,YPL,PHPL,ALPL) ! calcul des
3185 IF(IFLAG.EQ.1) THEN ! points d impacts dans les
3186 PRINT *,'FCN ',XPL(4,1),XMES(4) ! plans
3187 PRINT *,'FCN ',YPL(4,1),YMES(4)
3188 PRINT *,'FCN ',XPL(5,1),XMES(5)
3189 PRINT *,'FCN ',YPL(5,1),YMES(5)
3192 * IF (IVERTEX.EQ.1) THEN
3193 * FVAL = (PEST(4)/ERRV)**2 + (PEST(5)/ERRV)**2
3199 IF (MPOS(J).EQ.1) THEN
3201 XPLC = XPL(J,IZMES(J))
3202 YPLC = YPL(J,IZMES(J))
3203 FF = (XMES(J) - XPLC)/EEXM(J)
3205 FF = (YMES(J) - YPLC)/EEYM(J)
3208 IF (MANG(J).EQ.1) THEN
3210 FF = (PHMES(J) - PHPL(J))/EEPH(J)
3212 FF = (ALMES(J) - ALPL(J))/EEAL(J)
3216 ** PRINT *,'ST 1',XPL(1,1),XMES(1),YPL(1,1),YMES(1)
3217 ** PRINT *,'ST 2',XPL(2,1),XMES(2),YPL(2,1),YMES(2)
3218 ** PRINT *,'ST 3',XPL(3,1),XMES(3),YPL(3,1),YMES(3)
3219 ** PRINT *,'ST 4',XPL(4,1),XMES(4),YPL(4,1),YMES(4)
3220 ** PRINT *,'ST 5',XPL(5,1),XMES(5),YPL(5,1),YMES(5)
3222 IF (IVERTEX.EQ.1) NPARAM = 5
3223 CHI2PL = FVAL/FLOAT(2*NPLPL-NPARAM)
3227 ****************************************************************************
3228 SUBROUTINE CHECK_HISTO4(IDHIST,ICH2,IHIT2,ICH1,IHIT1,
3229 & X1,Y1,PHI1,ALAM1,P1,EXM,EYM,EPH,EAL)
3230 *****************************************************************************
3231 * Check hit IHIT2 with GEANT informations from hit HIT1
3233 * INPUT : ICH2 : No st. de recherche
3234 * IDCH1,X1,Y1,PHI1,ALAM1,P1 : Trace de reference
3235 * OUTPUT : JOK : No hit dans ICH2 appartenant a la meme trace.
3237 *****************************************************************************
3238 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
3240 PARAMETER (MAXHITCH=10000,MAXHITTOT=20000,NBSTATION=5,
3243 COMMON/CHHIT/XM(NBSTATION,MAXHITCH),YM(NBSTATION,MAXHITCH),
3244 & PHM(NBSTATION,MAXHITCH),ALM(NBSTATION,MAXHITCH),
3245 & IZM(NBSTATION,MAXHITCH),
3246 & IP(NBSTATION,MAXHITCH),JHIT(NBSTATION),
3247 & XMR(NBSTATION,MAXHITCH,2),YMR(NBSTATION,MAXHITCH,2)
3249 COMMON/RHIT/ITYP(MAXHITTOT),XTR(MAXHITTOT),YTR(MAXHITTOT),
3250 & PTOT(MAXHITTOT),ID(MAXHITTOT),IZST(MAXHITTOT),
3251 & PVERT1(MAXHITTOT),PVERT2(MAXHITTOT),PVERT3(MAXHITTOT),
3252 & ZVERT(MAXHITTOT),NHITTOT
3254 COMMON/VERIFGEANT/ITTROUGH(MAXTRK,NBSTATION),
3255 & IT_LIST(MAXTRK),IT_NP(MAXTRK),ITCHECK(MAXTRK),
3258 COMMON/DEBEVT/IDEBUG
3266 IF (PHM(ICH2,I).LE.6.3) THEN ! vector measurement
3267 IF (ID(IP(ICH1,IHIT1)).EQ.ID(IP(ICH2,I))) THEN
3269 IF (IDHIST.GT.0) THEN
3270 * CALL CHF1(IDHIST,SNGL(P1),SNGL((X1-XM(ICH2,I))**2))
3271 * CALL CHF1(IDHIST+1,SNGL(P1),SNGL((Y1-YM(ICH2,I))**2))
3272 * CALL CHF1(IDHIST+2,SNGL(P1),
3273 * & SNGL((PHI1-PHM(ICH2,I))**2))
3274 * CALL CHF1(IDHIST+3,SNGL(P1),
3275 * & SNGL((ALAM1-ALM(ICH2,I))**2))
3276 * CALL CHF1(IDHIST+4,SNGL(P1),R2)
3283 IF (ITCHECK(ITRACK(IP(ICH1,IHIT1))).EQ.1) THEN
3284 IF (IDEBUG.GE.2) THEN
3285 IF (IHIT2.EQ.0) THEN
3286 PRINT *,'CHECK4 histo nb:',IDHIST
3287 PRINT *,'CHECK4 p de st.',ICH1,'=',P1
3288 PRINT *,'CHECK4 track not found in st.',ICH2
3289 PRINT *,'CHECK4 error X :',(XM(ICH2,JOK)-X1), EXM
3290 PRINT *,'CHECK4 error Y :',(YM(ICH2,JOK)-Y1), EYM
3291 PRINT *,'CHECK4 error PHI :',(PHM(ICH2,JOK)-PHI1),EPH
3292 PRINT *,'CHECK4 error ALAM :',(ALM(ICH2,JOK)-ALAM1),
3294 ELSEIF (IHIT2.NE.JOK) THEN
3295 PRINT *,'CHECK4 histo nb:',IDHIST
3296 PRINT *,'CHECK4 p de st.',ICH1,'=',P1
3297 PRINT *,'CHECK4 ghost in st.',ICH2
3298 PRINT *,'CHECK4 id part. recherchee:',
3299 & ID(IP(ICH1,IHIT1))
3300 PRINT *,'CHECK4 id ghost trouve :',
3301 & ID(IP(ICH2,IHIT2))
3302 PRINT *,'CHECK4 JOK=',JOK,' IHIT2=',IHIT2
3310 *****************************************************************************
3311 SUBROUTINE CHECK_HISTO2(IDHIST,ICH2,IHIT2,ICH1,IHIT1,
3312 & X1,Y1,X2,Y2,P1,EX1,EY1,EX2,EY2)
3313 *****************************************************************************
3314 * Check hit IHIT2 with GEANT informations from hit HIT1
3316 * INPUT : IDHIST : No histo
3317 * ICH2 : No st. de recherche
3318 * IDCH1,X1,Y1,PHI1,ALAM1,P1 : Trace de reference
3319 * OUTPUT : JOK : No hit dans ICH2 appartenant a la meme trace.
3321 *****************************************************************************
3322 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
3324 PARAMETER (MAXHITCH=10000,MAXHITTOT=20000,NBSTATION=5,
3327 COMMON/CHHIT/XM(NBSTATION,MAXHITCH),YM(NBSTATION,MAXHITCH),
3328 & PHM(NBSTATION,MAXHITCH),ALM(NBSTATION,MAXHITCH),
3329 & IZM(NBSTATION,MAXHITCH),
3330 & IP(NBSTATION,MAXHITCH),JHIT(NBSTATION),
3331 & XMR(NBSTATION,MAXHITCH,2),YMR(NBSTATION,MAXHITCH,2)
3333 COMMON/RHIT/ITYP(MAXHITTOT),XTR(MAXHITTOT),YTR(MAXHITTOT),
3334 & PTOT(MAXHITTOT),ID(MAXHITTOT),IZST(MAXHITTOT),
3335 & PVERT1(MAXHITTOT),PVERT2(MAXHITTOT),PVERT3(MAXHITTOT),
3336 & ZVERT(MAXHITTOT),NHITTOT
3338 COMMON/VERIFGEANT/ITTROUGH(MAXTRK,NBSTATION),
3339 & IT_LIST(MAXTRK),IT_NP(MAXTRK),ITCHECK(MAXTRK),
3342 COMMON/DEBEVT/IDEBUG
3350 IF (ID(IP(ICH1,IHIT1)).EQ.ID(IP(ICH2,I))) THEN
3352 IF (IDHIST.GT.0) THEN
3353 IF (IZM(ICH2,I).EQ.1) THEN ! 1st chamber
3360 CALL CHF1(IDHIST,SNGL(P1),SNGL((X-XM(ICH2,I))**2))
3361 CALL CHF1(IDHIST+1,SNGL(P1),SNGL((Y-YM(ICH2,I))**2))
3362 CALL CHF1(IDHIST+4,SNGL(P1),R2)
3368 IF (ITCHECK(ITRACK(IP(ICH1,IHIT1))).EQ.1) THEN
3371 IF (IZM(ICH2,JOK).EQ.1) THEN
3374 IF (ICH2.EQ.4.OR.ICH2.EQ.5) THEN
3381 IF (ICH2.EQ.4.OR.ICH2.EQ.5) THEN
3386 IF (IDEBUG.GE.2) THEN
3387 IF (IHIT2.EQ.0) THEN
3388 PRINT *,'CHECK2 histo nb:',IDHIST
3389 PRINT *,'CHECK2 p de st.',ICH1,'=',P1
3390 PRINT *,'CHECK2 track not found in st.',ICH2
3391 PRINT *,'CHECK2 error X :',(XM(ICH2,JOK)-X), EXM
3392 PRINT *,'CHECK2 error Y :',(YM(ICH2,JOK)-Y), EYM
3393 ELSEIF(IHIT2.NE.JOK) THEN
3394 PRINT *,'CHECK2 histo nb:',IDHIST
3395 PRINT *,'CHECK2 p de st.',ICH1,'=',P1
3396 PRINT *,'CHECK2 ghost in st.',ICH2
3397 PRINT *,'CHECK2 id part. recherchee:',
3398 & ID(IP(ICH1,IHIT1))
3399 PRINT *,'CHECK2 id ghost trouve :',
3400 & ID(IP(ICH2,IHIT2))
3401 PRINT *,'CHECK2 JOK=',JOK,' IHIT2=',IHIT2
3410 DOUBLE PRECISION FUNCTION DEDX(P,THET,XEA,YEA)
3411 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3412 REA = DSQRT(XEA**2+YEA**2)
3413 IF (REA.lT.26.3611) THEN
3414 if (p .lt. 15.) then
3415 DP=2.737+0.0494*p-0.001123*p*p
3420 if (p .lt. 15.) then
3421 DP = 2.1380+0.0351*p-0.000853*p*p
3423 DP = 2.407+0.00702*p
3431 ************************************************************************
3432 DOUBLE PRECISION FUNCTION DEDX_oldold(P,THET,XEA,YEA)
3433 ************************************************************************
3434 * DEDX est la nouvelle impulsion au vertex, corrigee de la perte
3435 * d'energie dans l'absorbeur
3437 ************************************************************************
3439 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3440 DIMENSION PPB(6), PW(6)
3441 * FIT RESULT FOR PB REGION (5TH ORDER POLY)
3442 * 1 p0 2.24358e+03 2.55765e+00 1.06982e-03 6.30474e-07
3443 * 2 p1 1.16393e+01 4.45081e-02 6.58627e-06 -1.15822e-04
3444 * 3 p2 -1.82314e-01 3.28429e-04 8.69340e-08 -1.27281e-02
3445 * 4 p3 1.60930e-03 2.23812e-06 7.67374e-10 -4.15573e-01
3446 * 5 p4 -6.96885e-06 1.35405e-08 3.32301e-12 1.22136e+02
3447 * 6 p5 1.16339e-08 5.91665e-11 1.06532e-14 3.33965e+04
3449 DATA PPB /2.24358d+03, 1.16393d+01, -1.82314d-01, 1.60930d-03
3450 + ,-6.96885d-06, 1.16339d-08/
3451 * FIT RESULT FOR W REGION (5TH ORDER POLY)
3452 * 1 p0 2.90155e+03 3.49066e+00 1.38357e-03 -2.79916e-05
3453 * 2 p1 1.57716e+01 6.09946e-02 9.03687e-06 -6.63098e-03
3454 * 3 p2 -2.48349e-01 4.50365e-04 1.18422e-07 -8.27199e-01
3455 * 4 p3 2.19908e-03 3.07148e-06 1.04860e-09 -1.03290e+02
3456 * 5 p4 -9.54046e-06 1.85908e-08 4.54924e-12 -1.51284e+04
3457 * 6 p5 1.59463e-08 8.11346e-11 1.46446e-14 -2.69491e+06
3458 DATA PW /2.90155d+03, 1.57716d+01, -2.48349d-01, 2.19908d-03
3459 + , -9.54046d-06, 1.59463d-08/
3461 REA = DSQRT(XEA**2+YEA**2)
3462 IF (REA.GT.26.3611) THEN
3463 DP=PPB(1)+PPB(2)*P+PPB(3)*P**2
3464 & +PPB(4)*P**3+PPB(5)*P**4+PPB(6)*P**5
3466 DP=PW(1)+PW(2)*P+PW(3)*P**2
3467 & +PW(4)*P**3+PW(5)*P**4+PW(6)*P**5
3469 P=P+DP/1000.D0/DCOS(THET)
3473 ************************************************************************
3474 DOUBLE PRECISION FUNCTION DEDX_OLD(P,THET,XEA,YEA)
3475 ************************************************************************
3476 * DEDX est la nouvelle impulsion au vertex, corrigee de la perte
3477 * d'energie dans l'absorbeur
3479 ************************************************************************
3481 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3483 REA = DSQRT(XEA**2+YEA**2)
3484 IF (REA.GT.26.3611) THEN
3487 P = P+SPB/1000.*11.35*(16.66D-3 * P+1.33)
3490 P = P+SPO/1000.*0.935*(2.22D-3 * P+2.17)
3493 P = P+SPB/1000.*11.35*(16.66D-3 * P+1.33)
3496 P = P+SPO/1000.*0.935*(2.22D-3 * P+2.17)
3499 P = P+SPB/1000.*11.35*(16.66D-3 * P+1.33)
3502 P = P+SPO/1000.*0.935*(2.22D-3 * P+2.17)
3505 P = P+SPB/1000.*11.35*(16.66D-3 * P+1.33)
3509 SW = (503.-467.)/DCOS(THET)
3510 P = P+SW/1000.*19.3*(16.66D-3 * P+1.33)
3514 SCONC = (467.-315.)/DCOS(THET)
3515 P = P+SCONC/1000.*2.5*(2.22D-3*P+2.17)
3518 SC = (315.-90.)/DCOS(THET)
3519 P = P+SC/1000.*1.93*(2.22D-3*P+2.17) ! Carbone
3528 ************************************************************************
3529 SUBROUTINE BRANSON(PXZ,PHI,ALAM,XEA,YEA)
3530 ************************************************************************
3532 * Correction de Branson du multiple scattering dans l'absorbeur
3534 ************************************************************************
3536 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3538 PARAMETER(NPLANE=10)
3539 COMMON/PARAM/ZPLANEP(NPLANE),THICK,XPREC,YPREC,B0,BL3,ZMAGS,
3540 & ZMAGE,ZABS,XMAG,ZBP1,ZBP2,CONST
3543 IF (PXZ.LT.0.) ASIGN = -1.
3549 PTOT = PXZ/DCOS(ALAM)
3553 REA = DSQRT(XEA**2+YEA**2)
3554 IF (REA.GT.26.3611) THEN
3556 ELSE ! Abso. W pour theta < 3 deg
3559 * ZBP = ZEA ! Andreas
3560 XBP = XEA-PX/PZ*(ZEA-ZBP)
3561 YBP = YEA-PY/PZ*(ZEA-ZBP)
3562 PZ = PTOT*ZBP/DSQRT(XBP**2+YBP**2+ZBP**2)
3565 PXZ = DSQRT(PX**2+PZ**2)
3567 ALAM = DATAN2(PY,PXZ)
3569 ** THET = DATAN2(REA,ZEA)
3571 PT = DSQRT(PX**2+PY**2)
3572 THET = DATAN2(PT,PZ)
3573 PTOT = DEDX(PTOT,THET,XEA,YEA)
3575 PXZ = ASIGN*PTOT*DCOS(ALAM)
3580 ***************************************************************
3581 SUBROUTINE DSINV(N,A,IDIM,IFAIL)
3582 ***************************************************************
3584 DOUBLE PRECISION A(IDIM,*), ZERO, ONE, X, Y
3589 DOUBLE PRECISION S1, S31, S32, S33, DOTF
3592 DOTF(X,Y,S1) = X * Y + S1
3594 DATA HNAME / 'DSINV ' /
3595 DATA ZERO, ONE / 0.D0, 1.D0 /
3597 IF(IDIM .LT. N .OR. N .LE. 0) GOTO 900
3603 IF(PIVOTF(A(J,J)) .LE. 0.) GOTO 150
3604 A(J,J) = ONE / A(J,J)
3605 IF(J .EQ. N) GOTO 199
3608 A(J,L) = A(J,J)*A(L,J)
3611 S1 = DOTF(A(L,I),A(I,J+1),S1)
3622 IF(N .EQ. 1) GOTO 399
3624 A(2,1) = A(1,2)*A(2,2)
3625 IF(N .EQ. 2) GOTO 320
3631 S31 = DOTF(A(K,I+1),A(I+1,J),S31)
3634 A(J,K) = -S31*A(J,J)
3636 A(J-1,J) = -A(J-1,J)
3637 A(J,J-1) = A(J-1,J)*A(J,J)
3641 IF(J .EQ. N) GOTO 325
3644 S33 = DOTF(A(J,I),A(I,J),S33)
3652 S32 = DOTF(A(K,I),A(I,J),S32)
3657 IF(J .LT. N) GOTO 323
3661 900 CALL TMPRNT(HNAME,N,IDIM,0)
3665 *******************************************************
3666 SUBROUTINE TMPRNT(NAME,N,IDIM,K)
3667 *******************************************************
3670 LOGICAL MFLAG, RFLAG
3672 IF(NAME(2:2) .EQ. 'S') THEN
3673 CALL KERMTR('F012.1',LGFILE,MFLAG,RFLAG)
3675 CALL KERMTR('F011.1',LGFILE,MFLAG,RFLAG)
3677 IF(NAME(3:6) .EQ. 'FEQN') ASSIGN 1002 TO IFMT
3678 IF(NAME(3:6) .NE. 'FEQN') ASSIGN 1001 TO IFMT
3680 IF(LGFILE .EQ. 0) THEN
3681 IF(NAME(3:6) .EQ. 'FEQN') THEN
3682 WRITE(*,IFMT) NAME, N, IDIM, K
3684 WRITE(*,IFMT) NAME, N, IDIM
3687 IF(NAME(3:6) .EQ. 'FEQN') THEN
3688 WRITE(LGFILE,IFMT) NAME, N, IDIM, K
3690 WRITE(LGFILE,IFMT) NAME, N, IDIM
3694 IF(.NOT. RFLAG) CALL ABEND
3696 1001 FORMAT(7X, 31H PARAMETER ERROR IN SUBROUTINE , A6,
3697 + 27H ... (N.LT.1 OR IDIM.LT.N).,
3698 + 5X, 3HN =, I4, 5X, 6HIDIM =, I4, 1H. )
3699 1002 FORMAT(7X, 31H PARAMETER ERROR IN SUBROUTINE , A6,
3700 + 37H ... (N.LT.1 OR IDIM.LT.N OR K.LT.1).,
3701 + 5X, 3HN =, I4, 5X, 6HIDIM =, I4, 5X, 3HK =, I4,1H.)
3707 * Revision 1.5 2000/06/15 07:58:49 morsch
3708 * Code from MUON-dev joined
3710 * Revision 1.4.4.2 2000/04/26 15:48:37 morsch
3711 * Some routines from obsolete algo.F are needed by reco_muon.F and have been
3714 * Revision 1.4.4.1 2000/01/12 16:00:55 morsch
3715 * New version of MUON code
3717 * Revision 1.1.1.1 1996/02/15 17:48:35 mclareni
3722 ***********************************************************
3723 SUBROUTINE KERSET(ERCODE,LGFILE,LIMITM,LIMITR)
3724 ***********************************************************
3726 PARAMETER(KOUNTE = 27)
3727 CHARACTER*6 ERCODE, CODE(KOUNTE)
3728 LOGICAL MFLAG, RFLAG
3729 INTEGER KNTM(KOUNTE), KNTR(KOUNTE)
3732 DATA CODE(1), KNTM(1), KNTR(1) / 'C204.1', 255, 255 /
3733 DATA CODE(2), KNTM(2), KNTR(2) / 'C204.2', 255, 255 /
3734 DATA CODE(3), KNTM(3), KNTR(3) / 'C204.3', 255, 255 /
3735 DATA CODE(4), KNTM(4), KNTR(4) / 'C205.1', 255, 255 /
3736 DATA CODE(5), KNTM(5), KNTR(5) / 'C205.2', 255, 255 /
3737 DATA CODE(6), KNTM(6), KNTR(6) / 'C305.1', 255, 255 /
3738 DATA CODE(7), KNTM(7), KNTR(7) / 'C308.1', 255, 255 /
3739 DATA CODE(8), KNTM(8), KNTR(8) / 'C312.1', 255, 255 /
3740 DATA CODE(9), KNTM(9), KNTR(9) / 'C313.1', 255, 255 /
3741 DATA CODE(10),KNTM(10),KNTR(10) / 'C336.1', 255, 255 /
3742 DATA CODE(11),KNTM(11),KNTR(11) / 'C337.1', 255, 255 /
3743 DATA CODE(12),KNTM(12),KNTR(12) / 'C341.1', 255, 255 /
3744 DATA CODE(13),KNTM(13),KNTR(13) / 'D103.1', 255, 255 /
3745 DATA CODE(14),KNTM(14),KNTR(14) / 'D106.1', 255, 255 /
3746 DATA CODE(15),KNTM(15),KNTR(15) / 'D209.1', 255, 255 /
3747 DATA CODE(16),KNTM(16),KNTR(16) / 'D509.1', 255, 255 /
3748 DATA CODE(17),KNTM(17),KNTR(17) / 'E100.1', 255, 255 /
3749 DATA CODE(18),KNTM(18),KNTR(18) / 'E104.1', 255, 255 /
3750 DATA CODE(19),KNTM(19),KNTR(19) / 'E105.1', 255, 255 /
3751 DATA CODE(20),KNTM(20),KNTR(20) / 'E208.1', 255, 255 /
3752 DATA CODE(21),KNTM(21),KNTR(21) / 'E208.2', 255, 255 /
3753 DATA CODE(22),KNTM(22),KNTR(22) / 'F010.1', 255, 0 /
3754 DATA CODE(23),KNTM(23),KNTR(23) / 'F011.1', 255, 0 /
3755 DATA CODE(24),KNTM(24),KNTR(24) / 'F012.1', 255, 0 /
3756 DATA CODE(25),KNTM(25),KNTR(25) / 'F406.1', 255, 0 /
3757 DATA CODE(26),KNTM(26),KNTR(26) / 'G100.1', 255, 255 /
3758 DATA CODE(27),KNTM(27),KNTR(27) / 'G100.2', 255, 255 /
3762 IF(ERCODE .NE. ' ') THEN
3764 IF(ERCODE(1:L) .EQ. ERCODE) GOTO 12
3769 IF(L .EQ. 0) GOTO 13
3770 IF(CODE(I)(1:L) .NE. ERCODE(1:L)) GOTO 14
3771 13 IF(LIMITM.GE.0) KNTM(I) = LIMITM
3772 IF(LIMITR.GE.0) KNTR(I) = LIMITR
3775 ENTRY KERMTR(ERCODE,LOG,MFLAG,RFLAG)
3778 IF(ERCODE .EQ. CODE(I)) GOTO 21
3780 WRITE(*,1000) ERCODE
3783 21 RFLAG = KNTR(I) .GE. 1
3784 IF(RFLAG .AND. (KNTR(I) .LT. 255)) KNTR(I) = KNTR(I) - 1
3785 MFLAG = KNTM(I) .GE. 1
3786 IF(MFLAG .AND. (KNTM(I) .LT. 255)) KNTM(I) = KNTM(I) - 1
3787 IF(.NOT. RFLAG) THEN
3788 IF(LOGF .LT. 1) THEN
3789 WRITE(*,1001) CODE(I)
3791 WRITE(LOGF,1001) CODE(I)
3794 IF(MFLAG .AND. RFLAG) THEN
3795 IF(LOGF .LT. 1) THEN
3796 WRITE(*,1002) CODE(I)
3798 WRITE(LOGF,1002) CODE(I)
3802 1000 FORMAT(' KERNLIB LIBRARY ERROR. ' /
3803 + ' ERROR CODE ',A6,' NOT RECOGNIZED BY KERMTR',
3804 + ' ERROR MONITOR. RUN ABORTED.')
3805 1001 FORMAT(/' ***** RUN TERMINATED BY CERN LIBRARY ERROR ',
3807 1002 FORMAT(/' ***** CERN LIBRARY ERROR CONDITION ',A6)
3810 **************************
3812 **************************
3817 ************************************************************************
3818 SUBROUTINE FCNFIT(NPAR, GRAD, FVAL, XVAL, IFLAG, FUTIL)
3819 ************************************************************************
3820 * With magnetic Field Map GRKUTA
3822 * Calcule FVAL: la fonction minimisee par MINUIT
3823 * With magnetic field map
3825 ************************************************************************
3826 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3828 PARAMETER(NPLANE=10)
3830 COMMON/PARAM/ZPLANEP(NPLANE),THICK,XPREC,YPREC,B0,BL3,ZMAGS,
3831 & ZMAGE,ZABS,XMAG,ZBP1,ZBP2,CONST
3833 COMMON/MEAS/LPLANE(NPLANE),XMP(NPLANE),YMP(NPLANE)
3835 COMMON/FCNOUT/PXZEA,ALAMEA,PHIEA,XEA,YEA,NPLU,CHI2
3838 DIMENSION GRAD(*),XVAL(*),AMS(500),DISTAZ(500)
3840 DIMENSION XP(NPLANE),YP(NPLANE),
3841 & COV(NPLANE,NPLANE),AP(NPLANE),COVY(NPLANE,NPLANE)
3842 DIMENSION VECT(7),VOUT(7)
3851 IF (XVAL(1).LT.0.) ASIGN = -1.
3854 PXZ = DABS(1./XVAL(1))
3859 PTOT = PXZ/DCOS(ALAM)
3881 R = SQRT(VECT(1)*VECT(1)+VECT(2)*VECT(2))
3888 ** PRINT *,' AV GRKUTA ASIGN',ASIGN,' THET',THET
3890 DO WHILE (Z.GE.ZABS.AND.Z.LT.ZPLANEP(ICH)
3891 & .AND.NSTEP.LE.NSTEPMAX)
3892 ** & .AND.(THETA*PITODEG).GT.2.
3893 ** & .AND. (THETA*PITODEG).LT.9.)
3895 ** WRITE(6,*) NSTEP,(VECT(I),I=1,7)
3896 ** CALL RECO_GRKUTA(ASIGN,STEP,VECT,VOUT) ! CCC
3897 CALL RECO_GHELIX(ASIGN,STEP,VECT,VOUT)
3902 R = SQRT(VECT(1)*VECT(1)+VECT(2)*VECT(2))
3904 IF (NSTEP.EQ.NSTEPMAX) RETURN
3905 XP(ICH) = VECT(1)-(Z-ZPLANEP(ICH))*VECT(4)/VECT(6)
3906 YP(ICH) = VECT(2)-(Z-ZPLANEP(ICH))*VECT(5)/VECT(6)
3908 AP(ICH) = (0.0136D0/PTOT)*DSQRT(AL)*(1+0.038D0*DLOG(AL))
3910 ** PRINT *,' AP GRKUTA ASIGN',ASIGN,' THET',THET
3913 ** Matrice de covariance
3916 IF (LPLANE(II).EQ.1) THEN
3921 IF (LPLANE(JJ).EQ.1) THEN
3927 COV(J,I) =COV(J,I) + XPREC**2
3930 * IF (I .EQ. 10 .AND. J .EQ. 10) PRINT *,'10 10 ',COV(J,I)
3932 ** COV(J,I) = COV(J,I)
3933 ** & + (ZPLANEP(II) + DISTAZ(L))*(ZPLANEP(JJ) +
3934 ** & DISTAZ(L))*AMS(L)**2
3938 & + (ZPLANEP(II)-ZPLANEP(K))*
3939 & (ZPLANEP(JJ)-ZPLANEP(K))*AP(K)**2
3940 * IF (I .EQ. 10 .AND. J .EQ. 10) PRINT *,'10 10 ',COV(J,I)
3943 COVY(J,I) = COV(J,I)
3945 COVY(J,I) = COVY(J,I) - XPREC**2 + YPREC**2
3952 * Inversion des matrices de covariance
3956 CALL DSINV(NPLU, COV, NPLANE, IFAIL)
3957 ** IF (JFAIL.NE.0 .AND. IFAIL .NE. 0) STOP 'ERROR'
3958 IF (IFAIL .NE. 0) STOP 'ERROR'
3960 CALL DSINV(NPLU, COVY, NPLANE, IFAIL)
3961 ** IF (JFAIL.NE.0 .AND. IFAIL .NE. 0) STOP 'ERROR'
3962 IF (IFAIL .NE. 0) STOP 'ERROR'
3963 * PRINT *,' COVARIANCE MATRIX AFTER'
3965 * PRINT *,(COV(J,I),J=1,NPLANE)
3968 ** Calcul de FVAL ou CHI2
3972 IF (LPLANE(II).EQ.1) THEN
3977 IF (LPLANE(JJ).EQ.1) THEN
3980 FVAL = FVAL + COV(J,I)*(XMP(II)-XP(II))*(XMP(JJ)-XP(JJ))
3981 FVAL = FVAL + COVY(J,I)*(YMP(II)-YP(II))
3983 ** IF (JJ.EQ.II) THEN
3984 ** FVAL = FVAL + (XM(II)-XP(II))*(XM(JJ)-XP(JJ))/XPREC**2
3985 ** FVAL = FVAL + (YM(II)-YP(II))
3986 ** & *(YM(JJ)-YP(JJ))/YPREC**2
3994 ** IF (CHI2.GT.1.E4) THEN
3995 ** PRINT *,'FCNFIT CHI2= ',CHI2
4000 1000 FORMAT(I5,7F12.6)
4005 ************************************************************************
4006 SUBROUTINE NEWFCNFIT(NPAR, GRAD, FVAL, XVAL, IFLAG, FUTIL)
4007 ************************************************************************
4008 * With magnetic Field Map GRKUTA
4010 * Calcule FVAL: la fonction minimisee par MINUIT
4011 * With magnetic field map
4013 ************************************************************************
4014 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
4016 PARAMETER(NPLANE=10)
4018 COMMON/PARAM/ZPLANEP(NPLANE),THICK,XPREC,YPREC,B0,BL3,ZMAGS,
4019 & ZMAGE,ZABS,XMAG,ZBP1,ZBP2,CONST
4021 COMMON/MEAS/LPLANE(NPLANE),XMP(NPLANE),YMP(NPLANE)
4023 COMMON/FCNOUT/PXZEA,ALAMEA,PHIEA,XEA,YEA,NPLU,CHI2
4025 DIMENSION GRAD(*),XVAL(*),AMS(500),DISTAZ(500)
4027 DIMENSION XP(NPLANE),YP(NPLANE),
4028 & COV(NPLANE,NPLANE),AP(NPLANE),COVY(NPLANE,NPLANE)
4029 DIMENSION VECT(7),VOUT(7)
4038 IF (XVAL(1).LT.0.) ASIGN = -1.
4041 PXZ = DABS(1./XVAL(1))
4046 PTOT = PXZ/DCOS(ALAM)
4068 R = SQRT(VECT(1)*VECT(1)+VECT(2)*VECT(2))
4075 ** PRINT *,' AV GRKUTA ASIGN',ASIGN,' THET',THET
4077 DO WHILE (Z.GE.ZABS.AND.Z.LT.ZPLANEP(ICH)
4078 & .AND.NSTEP.LE.NSTEPMAX)
4079 ** & .AND.(THETA*PITODEG).GT.2.
4080 ** & .AND. (THETA*PITODEG).LT.9.)
4082 ** WRITE(6,*) NSTEP,(VECT(I),I=1,7)
4083 CALL RECO_GRKUTA (ASIGN,STEP,VECT,VOUT)
4088 R = SQRT(VECT(1)*VECT(1)+VECT(2)*VECT(2))
4090 IF (NSTEP.EQ.NSTEPMAX) RETURN
4091 XP(ICH) = VECT(1)-(Z-ZPLANEP(ICH))*VECT(4)/VECT(6)
4092 YP(ICH) = VECT(2)-(Z-ZPLANEP(ICH))*VECT(5)/VECT(6)
4094 AP(ICH) = (0.0136D0/PTOT)*DSQRT(AL)*(1+0.038D0*DLOG(AL))
4096 ** PRINT *,' AP GRKUTA ASIGN',ASIGN,' THET',THET
4099 ** Matrice de covariance
4102 IF (LPLANE(II).EQ.1) THEN
4107 IF (LPLANE(JJ).EQ.1) THEN
4113 COV(J,I) =COV(J,I) + XPREC**2
4116 * IF (I .EQ. 10 .AND. J .EQ. 10) PRINT *,'10 10 ',COV(J,I)
4119 & + (ZPLANEP(II) + DISTAZ(L))*(ZPLANEP(JJ) +
4120 & DISTAZ(L))*AMS(L)**2
4124 & + (ZPLANEP(II)-ZPLANEP(K))*
4125 & (ZPLANEP(JJ)-ZPLANEP(K))*AP(K)**2
4126 * IF (I .EQ. 10 .AND. J .EQ. 10) PRINT *,'10 10 ',COV(J,I)
4129 COVY(J,I) = COV(J,I)
4131 COVY(J,I) = COVY(J,I) - XPREC**2 + YPREC**2
4138 * Inversion des matrices de covariance
4142 CALL DSINV(NPLU, COV, NPLANE, IFAIL)
4143 ** IF (JFAIL.NE.0 .AND. IFAIL .NE. 0) STOP 'ERROR'
4144 IF (IFAIL .NE. 0) STOP 'ERROR'
4146 CALL DSINV(NPLU, COVY, NPLANE, IFAIL)
4147 ** IF (JFAIL.NE.0 .AND. IFAIL .NE. 0) STOP 'ERROR'
4148 IF (IFAIL .NE. 0) STOP 'ERROR'
4149 * PRINT *,' COVARIANCE MATRIX AFTER'
4151 * PRINT *,(COV(J,I),J=1,NPLANE)
4154 ** Calcul de FVAL ou CHI2
4158 IF (LPLANE(II).EQ.1) THEN
4163 IF (LPLANE(JJ).EQ.1) THEN
4166 FVAL = FVAL + COV(J,I)*(XMP(II)-XP(II))*(XMP(JJ)-XP(JJ))
4167 FVAL = FVAL + COVY(J,I)*(YMP(II)-YP(II))
4169 ** IF (JJ.EQ.II) THEN
4170 ** FVAL = FVAL + (XM(II)-XP(II))*(XM(JJ)-XP(JJ))/XPREC**2
4171 ** FVAL = FVAL + (YM(II)-YP(II))
4172 ** & *(YM(JJ)-YP(JJ))/YPREC**2
4180 ** IF (CHI2.GT.1.E4) THEN
4181 ** PRINT *,'FCNFIT CHI2= ',CHI2
4186 1000 FORMAT(I5,7F12.6)
4191 ***********************************************************************
4192 SUBROUTINE INITFIELDOLD
4195 ***********************************************************************
4197 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4199 ** IMPLICIT REAL*8(A-H,O-Z)
4201 COMMON/DAT1/Z(81),X(81),Y(81,44),DX,DZ,LPX,LPY,LPZ
4202 COMMON/DAT2/BX(81,81,44),BY(81,81,44),BZ(81,81,44)
4203 COMMON/REG1/ZMAX,ZMIN,XMAX,XMIN
4204 COMMON/REG2/AY1,CY1,AY2,CY2
4205 ** REAL *4 BXP,BYP,BZP
4206 COMMON/SDAT1/ZP(51),RAD(10),FI(33),DZP,DFI,DR,YY0,LPPZ,NR,NFI
4207 COMMON/SDAT2/BXP(51,10,33),BYP(51,10,33),BZP(51,10,33)
4208 COMMON/SDAT4/B(2,2,32)
4209 COMMON/REG3/ZPMAX,ZPMIN,RMAX,RMIN
4210 cc COMMON/CONST/PI2,EPS
4212 1000 FORMAT(5(1X,D15.7))
4213 2000 FORMAT(5(1X,I5))
4214 READ(40,2000) LPX,LPY,LPZ
4215 READ(40,1000) (Z(K),K=1,81)
4216 READ(40,1000) (X(K),K=1,81)
4217 READ(40,1000) DX,DY,DZ
4218 READ(40,1000) ZMAX,ZMIN,XMAX,XMIN
4219 c write(*,*) 'zmin zmax',ZMIN,ZMAX
4220 c write(*,*) 'xmin xmax',XMIN,XMAX
4221 READ(40,1000) AY1,CY1,AY2,CY2
4222 c write(*,*) 'ay1,cy1,ay2,cy2', AY1,CY1,AY2,CY2
4223 cc READ(40,1000) PI2,EPS
4224 READ(40,1000) (((BX(K,L,M),K=1,81),L=1,81),M=1,44)
4225 READ(40,1000) (((BY(K,L,M),K=1,81),L=1,81),M=1,44)
4226 READ(40,1000) (((BZ(K,L,M),K=1,81),L=1,81),M=1,44)
4230 READ(40,2000) LPPZ,NR,NFI
4231 READ(40,1000) (ZP(K),K=1,51)
4232 READ(40,1000) (RAD(K),K=1,10)
4233 READ(40,1000) (FI(L),L=1,33)
4234 READ(40,1000) DZP,DFI,DR
4235 c write(*,*) 'dzp dfi dR',DZP,DFI,DR
4236 READ(40,1000) ZPMAX,ZPMIN,RMAX,RMIN
4237 c write(*,*) 'zmin zmax',ZPMIN,ZPMAX
4238 c write(*,*) 'Rmin Rmax',RMIN,RMAX
4239 READ(40,1000) (((BXP(K,L,M),K=1,51),L=1,10),M=1,33)
4240 READ(40,1000) (((BYP(K,L,M),K=1,51),L=1,10),M=1,33)
4241 READ(40,1000) (((BZP(K,L,M),K=1,51),L=1,10),M=1,33)
4242 READ(40,1000) (((B(K,L,M),K=1,2),L=1,2),M=1,32)
4252 ***********************************************************************
4253 SUBROUTINE RECO_GUFLDOLD(X,F)
4254 C ^^^^^^^^^^^^^^^^^^^^^^
4255 C field map G. Chabratova
4257 C Field map 31/05/99
4258 ***********************************************************************
4261 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4262 COMMON/MAGERR/IMAGERR
4274 CALL FREG1(Z0,X0,Y0,FZ0,FX0,FY0,IND)
4275 ** PRINT 3000,Z0,X0,Y0,FZ0,FX0,FY0,IND
4278 CALL FREG2(Z0,X0,Y0,FZ0,FX0,FY0,IND)
4283 ** PRINT 3000,Z0,X0,Y0,FZ0,FX0,FY0,IND
4285 1000 format(1x,'Attention!!! The point is out of range!!!')
4287 3000 FORMAT(1X,'Z=',D13.7,1X,'X=',D13.7,1X,'Y=',D13.7,1X,
4288 & 'BZ=',D13.7,1X,'BX=',D13.7,1X,'BY=',D13.7,1X,'IND=',I3)
4311 **************************************************
4312 SUBROUTINE FREG1(Z0,X0,Y0,FZ0,FX0,FY0,IND)
4313 **************************************************
4314 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4316 COMMON/DAT1/Z(81),X(81),Y(81,44),DX,DZ,LPX,LPY,LPZ
4317 COMMON/DAT2/BX(81,81,44),BY(81,81,44),BZ(81,81,44)
4318 COMMON/REG1/ZMAX,ZMIN,XMAX,XMIN
4319 COMMON/REG2/AY1,CY1,AY2,CY2
4324 IF(Z0.LT.ZMIN.OR.Z0.GT.ZMAX)GO TO 100
4325 IF(X0.LT.XMIN.OR.X0.GT.XMAX)GO TO 100
4328 2000 FORMAT(1X,'YY1=',D15.7,1X,'YY2=',D15.7,1X,'Y0=',D15.7)
4329 c PRINT 2000,YY1,YY2,Y0
4330 IF(Y0.LT.YY1)GO TO 100
4331 IF(Y0.GT.YY2)GO TO 100
4332 CALL FIZ(Z0,Z,DZ,KC,K0,Z1,Z2,81)
4333 CALL FIZ(X0,X,DX,LC,L0,X1,X2,81)
4334 DY=(YY2-YY1)/DFLOAT(LPY-1)
4338 IF(Y0.GE.(YY1+DFLOAT(M0)*DY).AND.Y0.LE.(YY1+DFLOAT(M0+1)*DY))
4342 Y2=(Y0-(YY1+DFLOAT(M0)*DY))/DY
4344 ** write(*,*) 'm0 Y1 Y2',m0,Y1,Y2
4345 ** print *,' k0=',k0,' l0=',l0,' m0=',m0
4346 ** print *,' z1=',z1,' z2=',z2
4347 FX1=Z1*BX(K0,L0,M0)+Z2*BX(K0+1,L0,M0)
4348 FX2=Z2*BX(K0+1,L0+1,M0)+Z1*BX(K0,L0+1,M0)
4350 GX1=Z1*BX(K0,L0,M0+1)+Z2*BX(K0+1,L0,M0+1)
4351 GX2=Z2*BX(K0+1,L0+1,M0+1)+Z1*BX(K0,L0+1,M0+1)
4354 FX1=Z1*BY(K0,L0,M0)+Z2*BY(K0+1,L0,M0)
4355 FX2=Z2*BY(K0+1,L0+1,M0)+Z1*BY(K0,L0+1,M0)
4357 GX1=Z1*BY(K0,L0,M0+1)+Z2*BY(K0+1,L0,M0+1)
4358 GX2=Z2*BY(K0+1,L0+1,M0+1)+Z1*BY(K0,L0+1,M0+1)
4361 FX1=Z1*BZ(K0,L0,M0)+Z2*BZ(K0+1,L0,M0)
4362 FX2=Z2*BZ(K0+1,L0+1,M0)+Z1*BZ(K0,L0+1,M0)
4364 GX1=Z1*BZ(K0,L0,M0+1)+Z2*BZ(K0+1,L0,M0+1)
4365 GX2=Z2*BZ(K0+1,L0+1,M0+1)+Z1*BZ(K0,L0+1,M0+1)
4371 1000 format(1x,'Attention!!! The point is out of range!!!')
4376 *************************************************
4377 SUBROUTINE FIZ(Z0,Z,DEL,KI,K0,Z1,Z2,NDZ)
4378 *************************************************
4379 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4380 ** CC DIMENSION Z(NDZ)
4386 * if (k0.gt.81) print*,'ndz=',ndz
4387 IF (KJ.GT.NDZ) THEN ! CCC
4391 DO 1 K=KJ,NDZ-1 ! CCC
4398 * print *,'K0=',K0,' Z0',z0, Z(K0), Z(K0+1),z2
4399 if (k0.gt.81) print*,'k0=',k0
4400 Z2=(Z0-Z(K0))/(Z(K0+1)-Z(K0))
4402 ** write(*,*) 'ko z1 z2', K0,Z1,Z2,' ki=',ki,' kj=',kj,' K=',K
4403 ** write(*,*)' NDZ Z(K0) Z(K0+1)',NDZ,Z(K0), Z(K0+1)
4407 ***************************************************
4408 SUBROUTINE FREG2(Z0,X0,Y0,FZ0,FX0,FY0,IND)
4409 **************************************************
4410 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4411 ** REAL *4 BXP,BYP,BZP
4412 COMMON/SDAT1/ZP(51),RAD(10),FI(33),DZP,DFI,DR,YY0,LPPZ,NR,NFI
4413 COMMON/SDAT2/BXP(51,10,33),BYP(51,10,33),BZP(51,10,33)
4414 COMMON/SDAT4/B(2,2,32)
4415 COMMON/REG3/ZPMAX,ZPMIN,RMAX,RMIN
4416 cc COMMON/CONST/PI2,EPS
4423 R0=DSQRT((X0-YY0)**2+Y0**2)
4424 c write (*,*)'ro=',R0
4425 IF(Z0.LT.ZPMIN.OR.Z0.GT.ZPMAX)GO TO 100
4426 IF(R0.LT.RMIN.OR.R0.GT.RMAX)GO TO 100
4427 IF(R0.LE.DR)GO TO 3000
4428 CALL FIZ(Z0,ZP,DZP,KP,K0,Z1,Z2,51)
4429 CALL FIZ(R0,RAD,DR,LP,L0,X1,X2,10)
4430 ** print *,' r0=',r0,' rad=',rad,' dr=',dr,' lp=',lp,' l0=',l0,
4431 ** & ' x1=',x1,' x2=',x2
4432 FI0=DACOS((X0-YY0)/R0)
4433 IF(Y0.LT.0.D+0)FI0=PI2-FI0
4434 CALL FIZ(FI0,FI,DFI,MP,M0,Y1,Y2,32)
4435 ** print *,' Apres FIZ',' k0=',k0,' l0=',l0,' m0=',m0
4436 FX1=Z1*BXP(K0,L0,M0)+Z2*BXP(K0+1,L0,M0)
4437 FX2=Z2*BXP(K0+1,L0+1,M0)+Z1*BXP(K0,L0+1,M0)
4439 GX1=Z1*BXP(K0,L0,M0+1)+Z2*BXP(K0+1,L0,M0+1)
4440 GX2=Z2*BXP(K0+1,L0+1,M0+1)+Z1*BXP(K0,L0+1,M0+1)
4443 FX1=Z1*BYP(K0,L0,M0)+Z2*BYP(K0+1,L0,M0)
4444 FX2=Z2*BYP(K0+1,L0+1,M0)+Z1*BYP(K0,L0+1,M0)
4446 GX1=Z1*BYP(K0,L0,M0+1)+Z2*BYP(K0+1,L0,M0+1)
4447 GX2=Z2*BYP(K0+1,L0+1,M0+1)+Z1*BYP(K0,L0+1,M0+1)
4450 FX1=Z1*BZP(K0,L0,M0)+Z2*BZP(K0+1,L0,M0)
4451 ** CCC FX2=Z2*BZ(K0+1,L0+1,M0)+Z1*BZ(K0,L0+1,M0)
4452 FX2=Z2*BZP(K0+1,L0+1,M0)+Z1*BZP(K0,L0+1,M0)
4454 GX1=Z1*BZP(K0,L0,M0+1)+Z2*BZP(K0+1,L0,M0+1)
4455 GX2=Z2*BZP(K0+1,L0+1,M0+1)+Z1*BZP(K0,L0+1,M0+1)
4462 1000 format(1x,'Attention!!! The point is out of range!!!')
4466 IF(R0.LT.EPS)GO TO 4000
4467 CALL FIZ(Z0,ZP,DZP,KP,K0,Z1,Z2,51)
4470 IF(Y0.LT.0.D+0)FI0=PI2-FI0
4471 CALL FIZ(FI0,FI,DFI,MP,M0,Y1,Y2,32)
4472 ALF2=B(1,1,M0)*XX+B(1,2,M0)*Y0
4473 ALF3=B(2,1,M0)*XX+B(2,2,M0)*Y0
4475 FX1=ALF1*BXP(K0,1,1)+ALF2*BXP(K0,1,M0)+ALF3*BXP(K0,1,M0+1)
4476 FX2=ALF1*BXP(K0+1,1,1)+ALF2*BXP(K0+1,1,M0)+ALF3*BXP(K0+1,1,M0+1)
4478 FX1=ALF1*BYP(K0,1,1)+ALF2*BYP(K0,1,M0)+ALF3*BYP(K0,1,M0+1)
4479 FX2=ALF1*BYP(K0+1,1,1)+ALF2*BYP(K0+1,1,M0)+ALF3*BYP(K0+1,1,M0+1)
4481 FX1=ALF1*BZP(K0,1,1)+ALF2*BZP(K0,1,M0)+ALF3*BZP(K0,1,M0+1)
4482 FX2=ALF1*BZP(K0+1,1,1)+ALF2*BZP(K0+1,1,M0)+ALF3*BZP(K0+1,1,M0+1)
4484 c write(*,*) 'R<Dr:B(1,1,m0) B(1,2,m0) B(2,1,m0) B(2,2,m0)',
4485 c +B(1,1,M0),B(1,2,M0),B(2,1,M0),B(2,2,M0)
4486 c write(*,*)'BX(K0,1,1) BX(K0,1,M0)','BX(K0,1,M0+1) BX(K0+1,1,1)'
4487 c + ,'BX(K0+1,1,M0) BX(K0+1,1,M0+1)', BX(K0,1,1),BX(K0,1,M0) ,
4488 c + BX(K0,1,M0+1),BX(K0+1,1,1),BX(K0+1,1,M0),BX(K0+1,1,M0+1)
4489 c write(*,*)'By(K0,1,1) By(K0,1,M0)','By(K0,1,M0+1) By(K0+1,1,1)'
4490 c + ,'By(K0+1,1,M0) By(K0+1,1,M0+1)', By(K0,1,1),By(K0,1,M0) ,
4491 c + By(K0,1,M0+1),By(K0+1,1,1),By(K0+1,1,M0),By(K0+1,1,M0+1)
4492 ccc write (*,*)'Bz(K0,1,1) Bz(K0,1,M0) Bz(K0,1,M0+1) Bz(K0+1,1,1)
4493 77 FORMAT(5x,E15.7,2x,E15.7)
4495 ** 70 FORMAT(22hBz(K0,1,1) Bz(K0,1,M0))
4496 cc PRINT 77 , BzP(K0,1,1),BzP(K0,1,M0)
4498 ** 71 FORMAT(26hBz(K0,1,M0+1) Bz(K0+1,1,1))
4499 cc PRINT 77, BzP(K0,1,M0+1),BzP(K0+1,1,1)
4501 ** 72 FORMAT(29hBz(K0+1,1,M0) Bz(K0+1,1,M0+1))
4502 cc PRINT 77 ,BzP(K0+1,1,M0),BzP(K0+1,1,M0+1)
4503 cc 77 FORMAT(5x,D15.7,5x,D15.7)
4509 CALL FIZ(Z0,ZP,DZP,KP,K0,Z1,Z2,51)
4510 FX0=Z1*BXP(K0,1,1)+Z2*BXP(K0+1,1,1)
4511 FY0=Z1*BYP(K0,1,1)+Z2*BYP(K0+1,1,1)
4512 FZ0=Z1*BZP(K0,1,1)+Z2*BZP(K0+1,1,1)
4513 c write(*,*) ' R<eps: Bx(k) Bx(k+1) By(k) By(k+1) Bz(k) Bz(k+1)',
4514 c + BXP(K0,1,1),BXP(K0+1,1,1),BYP(K0,1,1),BYP(K0+1,1,1),BZP(K0,1,1),
4520 ***********************************************************************
4522 SUBROUTINE RECO_GRKUTA (CHARGE,STEP,VECT,VOUT)
4524 C. ******************************************************************
4526 C. * Runge-Kutta method for tracking a particle through a magnetic *
4527 C. * field. Uses Nystroem algorithm (See Handbook Nat. Bur. of *
4528 C. * Standards, procedure 25.5.20) *
4530 C. * Input parameters *
4531 C. * CHARGE Particle charge *
4532 C. * STEP Step size *
4533 C. * VECT Initial co-ords,direction cosines,momentum *
4534 C. * Output parameters *
4535 C. * VOUT Output co-ords,direction cosines,momentum *
4536 C. * User routine called *
4537 C. * CALL GUFLD(X,F) *
4539 C. * ==>Called by : <USER>, GUSWIM *
4540 C. * Authors R.Brun, M.Hansroul ********* *
4541 C. * V.Perevoztchikov (CUT STEP implementation) *
4544 C. ******************************************************************
4546 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4548 ** REAL CHARGE, STEP, VECT(*), VOUT(*), F(4)
4549 ** REAL XYZT(3), XYZ(3), X, Y, Z, XT, YT, ZT
4550 DIMENSION VECT(*), VOUT(*), F(3)
4551 DIMENSION XYZT(3), XYZ(3)
4552 DIMENSION SECXS(4),SECYS(4),SECZS(4),HXP(3)
4553 EQUIVALENCE (X,XYZ(1)),(Y,XYZ(2)),(Z,XYZ(3)),
4554 + (XT,XYZT(1)),(YT,XYZT(2)),(ZT,XYZT(3))
4556 PARAMETER (MAXIT = 1992, MAXCUT = 11)
4557 PARAMETER (EC=2.9979251D-4,DLT=1D-4,DLT32=DLT/32)
4558 PARAMETER (ZERO=0, ONE=1, TWO=2, THREE=3)
4559 PARAMETER (THIRD=ONE/THREE, HALF=ONE/TWO)
4560 PARAMETER (PISQUA=.986960440109D+01)
4561 PARAMETER (IX=1,IY=2,IZ=3,IPX=4,IPY=5,IPZ=6)
4563 *. ------------------------------------------------------------------
4565 * This constant is for units CM,GEV/C and KGAUSS
4572 PINV = EC * CHARGE / VECT(7)
4578 IF (ABS(H).GT.ABS(REST)) H = REST
4579 CALL RECO_GUFLD(VOUT,F)
4581 * Start of integration
4594 SECXS(1) = (B * F(3) - C * F(2)) * PH2
4595 SECYS(1) = (C * F(1) - A * F(3)) * PH2
4596 SECZS(1) = (A * F(2) - B * F(1)) * PH2
4597 ANG2 = (SECXS(1)**2 + SECYS(1)**2 + SECZS(1)**2)
4598 IF (ANG2.GT.PISQUA) GO TO 40
4599 DXT = H2 * A + H4 * SECXS(1)
4600 DYT = H2 * B + H4 * SECYS(1)
4601 DZT = H2 * C + H4 * SECZS(1)
4606 * Second intermediate point
4608 EST = ABS(DXT)+ABS(DYT)+ABS(DZT)
4609 IF (EST.GT.H) GO TO 30
4611 CALL RECO_GUFLD(XYZT,F)
4616 SECXS(2) = (BT * F(3) - CT * F(2)) * PH2
4617 SECYS(2) = (CT * F(1) - AT * F(3)) * PH2
4618 SECZS(2) = (AT * F(2) - BT * F(1)) * PH2
4622 SECXS(3) = (BT * F(3) - CT * F(2)) * PH2
4623 SECYS(3) = (CT * F(1) - AT * F(3)) * PH2
4624 SECZS(3) = (AT * F(2) - BT * F(1)) * PH2
4625 DXT = H * (A + SECXS(3))
4626 DYT = H * (B + SECYS(3))
4627 DZT = H * (C + SECZS(3))
4631 AT = A + TWO*SECXS(3)
4632 BT = B + TWO*SECYS(3)
4633 CT = C + TWO*SECZS(3)
4635 EST = ABS(DXT)+ABS(DYT)+ABS(DZT)
4636 IF (EST.GT.2.*ABS(H)) GO TO 30
4638 CALL RECO_GUFLD(XYZT,F)
4640 Z = Z + (C + (SECZS(1) + SECZS(2) + SECZS(3)) * THIRD) * H
4641 Y = Y + (B + (SECYS(1) + SECYS(2) + SECYS(3)) * THIRD) * H
4642 X = X + (A + (SECXS(1) + SECXS(2) + SECXS(3)) * THIRD) * H
4644 SECXS(4) = (BT*F(3) - CT*F(2))* PH2
4645 SECYS(4) = (CT*F(1) - AT*F(3))* PH2
4646 SECZS(4) = (AT*F(2) - BT*F(1))* PH2
4647 A = A+(SECXS(1)+SECXS(4)+TWO * (SECXS(2)+SECXS(3))) * THIRD
4648 B = B+(SECYS(1)+SECYS(4)+TWO * (SECYS(2)+SECYS(3))) * THIRD
4649 C = C+(SECZS(1)+SECZS(4)+TWO * (SECZS(2)+SECZS(3))) * THIRD
4651 EST = ABS(SECXS(1)+SECXS(4) - (SECXS(2)+SECXS(3)))
4652 ++ ABS(SECYS(1)+SECYS(4) - (SECYS(2)+SECYS(3)))
4653 ++ ABS(SECZS(1)+SECZS(4) - (SECZS(2)+SECZS(3)))
4655 IF (EST.GT.DLT .AND. ABS(H).GT.1.E-4) GO TO 30
4658 * If too many iterations, go to HELIX
4659 IF (ITER.GT.MAXIT) GO TO 40
4662 IF (EST.LT.(DLT32)) THEN
4665 CBA = ONE/ SQRT(A*A + B*B + C*C)
4673 IF (STEP.LT.0.) REST = -REST
4674 IF (REST .GT. 1.E-5*ABS(STEP)) GO TO 20
4680 * If too many cuts , go to HELIX
4681 IF (NCUT.GT.MAXCUT) GO TO 40
4685 ** ANGLE TOO BIG, USE HELIX
4689 F4 = SQRT(F1**2+F2**2+F3**2)
4698 HXP(1) = F2*VECT(IPZ) - F3*VECT(IPY)
4699 HXP(2) = F3*VECT(IPX) - F1*VECT(IPZ)
4700 HXP(3) = F1*VECT(IPY) - F2*VECT(IPX)
4702 HP = F1*VECT(IPX) + F2*VECT(IPY) + F3*VECT(IPZ)
4706 COST = TWO*SIN(HALF*TET)**2
4710 G3 = (TET-SINT) * HP*RHO1
4715 VOUT(IX) = VECT(IX) + (G1*VECT(IPX) + G2*HXP(1) + G3*F1)
4716 VOUT(IY) = VECT(IY) + (G1*VECT(IPY) + G2*HXP(2) + G3*F2)
4717 VOUT(IZ) = VECT(IZ) + (G1*VECT(IPZ) + G2*HXP(3) + G3*F3)
4719 VOUT(IPX) = VECT(IPX) + (G4*VECT(IPX) + G5*HXP(1) + G6*F1)
4720 VOUT(IPY) = VECT(IPY) + (G4*VECT(IPY) + G5*HXP(2) + G6*F2)
4721 VOUT(IPZ) = VECT(IPZ) + (G4*VECT(IPZ) + G5*HXP(3) + G6*F3)
4724 VOUT(IX) = VECT(IX) + STEP*VECT(IPX)
4725 VOUT(IY) = VECT(IY) + STEP*VECT(IPY)
4726 VOUT(IZ) = VECT(IZ) + STEP*VECT(IPZ)
4733 *******************************************************************
4734 SUBROUTINE RECO_GUFLDOLD1(X,B)
4738 C *** ROUTINE DESCRIBING THE MAGNETIC FIELD IN THE ALICE SETUP ***
4739 C *** NVE 14-NOV-1990 CERN GENEVA ***
4742 C ORIGIN : NICK VAN EIJNDHOVEN
4746 C X = (X,Y,Z) coordinates in cm
4750 C B = Magnetic field components (BX,BY,BZ) in KG
4753 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4755 PARAMETER(NBSTATION=5)
4757 COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
4771 IF (X(3).LT.(-1.*ZCOIL)) THEN
4773 ELSEIF ( X(3).LT.(-1.*ZMAGEND)) THEN
4784 ** print *,' x =',X(1),X(2),X(3)
4785 ** print *,' B =',B(1),B(2),B(3)
4793 *******************************************************************
4794 SUBROUTINE RECO_GUFLD(X,B)
4798 C *** ROUTINE DESCRIBING THE MAGNETIC FIELD IN THE ALICE SETUP ***
4799 C *** NVE 14-NOV-1990 CERN GENEVA ***
4802 C ORIGIN : NICK VAN EIJNDHOVEN
4806 C X = (X,Y,Z) coordinates in cm
4810 C B = Magnetic field components (BX,BY,BZ) in KG
4813 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4816 C --- Common containing magnetic field map data
4817 REAL DZ,DX,DY,UDX,UDY,UDZ
4818 $,XMBEG,YMBEG,ZMBEG,XMEND,YMEND,ZMEND
4822 PARAMETER(MAXFLD=250000)
4823 COMMON /SCXMFD/ NX,NY,NZ,DZ,DX,DY,UDX,UDY,UDZ
4824 $,XMBEG,YMBEG,ZMBEG,XMEND,YMEND,ZMEND
4830 DOUBLE PRECISION ONE, RATX, RATY, RATZ, HIX, HIY, HIZ
4831 $, RATX1, RATY1, RATZ1
4832 $, BHYHZ, BHYLZ, BLYHZ, BLYLZ, BHZ, BLZ
4846 ** BX(JX,JY,JZ)=BV(3*((JZ-1)*(NX*NY)+(JY-1)*NX+(JX-1))+1)
4847 ** BY(JX,JY,JZ)=BV(3*((JZ-1)*(NX*NY)+(JY-1)*NX+(JX-1))+2)
4848 ** BZ(JX,JY,JZ)=BV(3*((JZ-1)*(NX*NY)+(JY-1)*NX+(JX-1))+3)
4853 * --- Act accordingly to ISXFMAP
4855 IF(ISXFMAP.EQ.1) THEN
4856 IF (ABS(X(3)) .LT. 700.
4857 + .AND. X(1)**2+(X(2)+30.)**2 .LT. 560.**2 ) THEN
4861 ELSE IF (X(3) .GE. 725. .AND. X(3) .LT. 1225.) THEN
4862 DZ=ABS(975.-X(3))/100.
4863 B(1)=(1.-0.1*DZ*DZ)*7.0
4871 ELSE IF(ISXFMAP.LE.3) THEN
4872 IF (-700.LT.X(3).AND.X(3).LT.ZMBEG
4873 + .AND. X(1)**2+(X(2)+30.)**2 .LT. 560.**2 ) THEN
4877 ELSE IF ((X(3) .GE. ZMBEG .AND. X(3) .LT. ZMEND) .AND.
4878 + (XMBEG.LE.ABS(X(1)).AND.ABS(X(1)).LT.XMEND) .AND.
4879 + (YMBEG.LE.ABS(X(2)).AND.ABS(X(2)).LT.YMEND)) THEN
4883 C --- find the position in the grid ---
4885 XL(1)=ABS(X(1))-XMBEG
4886 XL(2)=ABS(X(2))-YMBEG
4903 IF(ISXFMAP.EQ.2) THEN
4904 * ... Simple interpolation
4906 B(1) = BX(IX,IY,IZ)*(ONE-RATX) + BX(IX+1,IY+1,IZ+1)*RATX
4907 B(2) = BY(IX,IY,IZ)*(ONE-RATY) + BY(IX+1,IY+1,IZ+1)*RATY
4908 B(3) = BZ(IX,IY,IZ)*(ONE-RATZ) + BZ(IX+1,IY+1,IZ+1)*RATZ
4909 ELSE IF(ISXFMAP.EQ.3) THEN
4910 * ... more complicated interpolation
4914 ** print *,' bx by bz', BX(IX ,IY+1,IZ+1),BY(IX ,IY+1,IZ+1)
4915 ** & , BZ(IX ,IY+1,IZ+1)
4916 BHYHZ = BX(IX ,IY+1,IZ+1)*RATX1+BX(IX+1,IY+1,IZ+1)*RATX
4917 BHYLZ = BX(IX ,IY+1,IZ )*RATX1+BX(IX+1,IY+1,IZ )*RATX
4918 BLYHZ = BX(IX ,IY ,IZ+1)*RATX1+BX(IX+1,IY ,IZ+1)*RATX
4919 BLYLZ = BX(IX ,IY ,IZ )*RATX1+BX(IX+1,IY ,IZ )*RATX
4920 BHZ = BLYHZ *RATY1+BHYHZ *RATY
4921 BLZ = BLYLZ *RATY1+BHYLZ *RATY
4922 B(1) = BLZ *RATZ1+BHZ *RATZ
4924 BHYHZ = BY(IX ,IY+1,IZ+1)*RATX1+BY(IX+1,IY+1,IZ+1)*RATX
4925 BHYLZ = BY(IX ,IY+1,IZ )*RATX1+BY(IX+1,IY+1,IZ )*RATX
4926 BLYHZ = BY(IX ,IY ,IZ+1)*RATX1+BY(IX+1,IY ,IZ+1)*RATX
4927 BLYLZ = BY(IX ,IY ,IZ )*RATX1+BY(IX+1,IY ,IZ )*RATX
4928 BHZ = BLYHZ *RATY1+BHYHZ *RATY
4929 BLZ = BLYLZ *RATY1+BHYLZ *RATY
4930 B(2) = BLZ *RATZ1+BHZ *RATZ
4932 BHYHZ = BZ(IX ,IY+1,IZ+1)*RATX1+BZ(IX+1,IY+1,IZ+1)*RATX
4933 BHYLZ = BZ(IX ,IY+1,IZ )*RATX1+BZ(IX+1,IY+1,IZ )*RATX
4934 BLYHZ = BZ(IX ,IY ,IZ+1)*RATX1+BZ(IX+1,IY ,IZ+1)*RATX
4935 BLYLZ = BZ(IX ,IY ,IZ )*RATX1+BZ(IX+1,IY ,IZ )*RATX
4936 BHZ = BLYHZ *RATY1+BHYHZ *RATY
4937 BLZ = BLYLZ *RATY1+BHYLZ *RATY
4938 B(3) = BLZ *RATZ1+BHZ *RATZ
4941 * ... use the dipole symmetry
4943 IF (X(1)*X(2).LT.0) B(2)=-B(2)
4944 IF (X(1).LT.0) B(3)=-B(3)
4954 ENDIF ! z-coord for m.f.
4956 ENDIF ! endif ISXFMAP
4963 BTOT=SQRT(B(1)**2+B(2)**2+B(3)**2)
4964 IF(BTOT.GT.SXMGMX) THEN
4965 PRINT 10100, BTOT,SXMGMX
4966 10100 FORMAT(' *GUFLD* Field ',G10.4,' larger than max ',G10.4)
4984 *******************************************
4985 DOUBLE PRECISION FUNCTION BX(JX,JY,JZ)
4987 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4989 C --- Common containing magnetic field map data
4990 REAL DZ,DX,DY,UDX,UDY,UDZ
4991 $,XMBEG,YMBEG,ZMBEG,XMEND,YMEND,ZMEND
4995 PARAMETER(MAXFLD=250000)
4996 COMMON /SCXMFD/ NX,NY,NZ,DZ,DX,DY,UDX,UDY,UDZ
4997 $,XMBEG,YMBEG,ZMBEG,XMEND,YMEND,ZMEND
5000 BX=BV(3*((JZ-1)*(NX*NY)+(JY-1)*NX+(JX-1))+1)
5004 *******************************************
5005 DOUBLE PRECISION FUNCTION BY(JX,JY,JZ)
5007 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
5009 C --- Common containing magnetic field map data
5010 REAL DZ,DX,DY,UDX,UDY,UDZ
5011 $,XMBEG,YMBEG,ZMBEG,XMEND,YMEND,ZMEND
5015 PARAMETER(MAXFLD=250000)
5016 COMMON /SCXMFD/ NX,NY,NZ,DZ,DX,DY,UDX,UDY,UDZ
5017 $,XMBEG,YMBEG,ZMBEG,XMEND,YMEND,ZMEND
5020 BY=BV(3*((JZ-1)*(NX*NY)+(JY-1)*NX+(JX-1))+2)
5024 *******************************************
5025 DOUBLE PRECISION FUNCTION BZ(JX,JY,JZ)
5027 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
5029 C --- Common containing magnetic field map data
5030 REAL DZ,DX,DY,UDX,UDY,UDZ
5031 $,XMBEG,YMBEG,ZMBEG,XMEND,YMEND,ZMEND
5035 PARAMETER(MAXFLD=250000)
5036 COMMON /SCXMFD/ NX,NY,NZ,DZ,DX,DY,UDX,UDY,UDZ
5037 $,XMBEG,YMBEG,ZMBEG,XMEND,YMEND,ZMEND
5040 BZ=BV(3*((JZ-1)*(NX*NY)+(JY-1)*NX+(JX-1))+3)
5044 ***********************************************************************
5045 SUBROUTINE INITFIELD
5050 C *** INITIALISATION OF THE FIELD MAP ***
5051 C *** FCA 24-AUG-1998 CERN GENEVA ***
5053 C CALLED BY : GALICE
5054 C ORIGIN : NICK VAN EIJNDHOVEN
5056 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
5059 C --- Common containing magnetic field map data
5060 REAL DZ,DX,DY,UDX,UDY,UDZ
5061 $,XMBEG,YMBEG,ZMBEG,XMEND,YMEND,ZMEND
5065 PARAMETER(MAXFLD=250000)
5066 COMMON /SCXMFD/ NX,NY,NZ,DZ,DX,DY,UDX,UDY,UDZ
5067 $,XMBEG,YMBEG,ZMBEG,XMEND,YMEND,ZMEND
5072 IF(ISXFMAP.EQ.1) THEN
5073 * ... constant field, nothing to do
5074 ELSE IF(ISXFMAP.LE.3) THEN
5075 * ... constant mesh field
5076 PRINT 10000, ISXFMAP
5077 10000 FORMAT(' *SXFMAP* Magnetic field map flag: ',I5
5078 $ ,'; Reading magnetic field map data ')
5080 OPEN(77,FILE='/home/morsch/AliRoot/V3.02/data/field01.dat',
5081 $ FORM='FORMATTED',STATUS='OLD')
5082 READ(77,*) NX,NY,NZ,DX,DY,DZ,XMBEG,YMBEG,ZMBEG
5083 PRINT*,'NX,NY,NZ,DX,DY,DZ,XMBEG,YMBEG,ZMBEG',
5084 $ NX,NY,NZ,DX,DY,DZ,XMBEG,YMBEG,ZMBEG
5085 IF(3*NX*NY*NZ.GT.MAXFLD) THEN
5086 WRITE(6,10100) 3*NX*NY*NZ,MAXFLD
5087 STOP 'Increase MAXFLD'
5092 XMEND=XMBEG+(NX-1)*DX
5093 YMEND=YMBEG+(NY-1)*DY
5094 ZMEND=ZMBEG+(NZ-1)*DZ
5096 IPZ=3*(IZ-1)*(NX*NY)
5101 READ(77,*) BV(IPX+3),BV(IPX+2),BV(IPX+1)
5106 ENDIF ! endif ISXFMAP
5108 10100 FORMAT('*** SXFMAP: Need ',I7,' have ',I7)
5111 ****************************************
5112 SUBROUTINE RANNOR (A,B)
5113 ****************************************
5115 C CERN PROGLIB# V100 RANNOR .VERSION KERNFOR 4.18 880425
5119 IF (Y.EQ.0.) Y = RNDM()
5123 R = SQRT (-2.0*LOG(Y))
5128 ************************************************************************
5129 SUBROUTINE OLDFCNFIT(NPAR,GRAD,FVAL,XVAL,IFLAG,FUTIL)
5130 ************************************************************************
5131 * Calcule FVAL: la fonction minimisee par MINUIT
5132 * with constant magnetic Field
5134 ************************************************************************
5136 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
5138 PARAMETER(NPLANE=10)
5140 COMMON/PARAM/ZPLANEP(NPLANE),THICK,XPREC,YPREC,B0,BL3,ZMAGS,
5141 & ZMAGE,ZABS,XMAG,ZBP1,ZBP2,CONST
5143 COMMON/MEAS/LPLANE(NPLANE),XMP(NPLANE),YMP(NPLANE)
5145 COMMON/FCNOUT/PXZEA,ALAMEA,PHIEA,XEA,YEA,NPLU,CHI2
5147 DIMENSION GRAD(*),XVAL(*),AMS(500),DISTAZ(500)
5149 DIMENSION XP(NPLANE),YP(NPLANE),
5150 & COV(NPLANE,NPLANE),AP(NPLANE),COVY(NPLANE,NPLANE)
5151 DIMENSION VECT(7),VOUT(7)
5157 IF (XVAL(1).LT.0.) ASIGN = -1.
5160 PXZ = DABS(1./XVAL(1))
5165 PTOT = PXZ/DCOS(ALAM)
5166 TTHET = DSQRT(DTAN(ALAM)**2+DSIN(PHI)**2)/DCOS(PHI)
5168 PT = DSQRT(PX**2+PY**2)
5170 RL3 = ASIGN*PT / (0.299792458D-3 * BL3)
5171 ALPHA = DATAN2(PY,PX)
5172 XC = XV+RL3*DSIN(ALPHA)
5173 YC = YV-RL3*DCOS(ALPHA)
5187 ANGDEV = (ZPLANEP(1)-ZEA)*TTHET/RL3
5188 XP(1) = XC+RL3*DSIN(ANGDEV-ALPHA)
5189 YP(1) = YC+RL3*DCOS(ANGDEV-ALPHA)
5190 AL = THICK/ DCOS(THET)
5191 AP(1) = (0.0136D0/PTOT) * DSQRT(AL) * (1 + 0.038D0*DLOG(AL))
5193 ANGDEV = (ZPLANEP(2)-ZEA)*TTHET/RL3
5194 XP(2) = XC+RL3*DSIN(ANGDEV-ALPHA)
5195 YP(2) = YC+RL3*DCOS(ANGDEV-ALPHA)
5196 AL = THICK/ DCOS(THET)
5197 AP(2) = (0.0136D0/PTOT) * DSQRT(AL) * (1 + 0.038D0*DLOG(AL))
5199 ANGDEV = (ZPLANEP(3)-ZEA)*TTHET/RL3
5200 XP(3) = XC+RL3*DSIN(ANGDEV-ALPHA)
5201 YP(3) = YC+RL3*DCOS(ANGDEV-ALPHA)
5202 AL = THICK/ DCOS(THET)
5203 AP(3) = (0.0136D0/PTOT) * DSQRT(AL) * (1 + 0.038D0*DLOG(AL))
5205 ANGDEV = (ZPLANEP(4)-ZEA)*TTHET/RL3
5206 XP(4) = XC+RL3*DSIN(ANGDEV-ALPHA)
5207 YP(4) = YC+RL3*DCOS(ANGDEV-ALPHA)
5208 AL = THICK/ DCOS(THET)
5209 AP(4) = (0.0136D0/PTOT) * DSQRT(AL) * (1 + 0.038D0*DLOG(AL))
5211 ANGDEV = (700.D0-ZEA)*TTHET/RL3
5212 XPL3 = XC+RL3*DSIN(ANGDEV-ALPHA)
5213 YPL3 = YC+RL3*DCOS(ANGDEV-ALPHA)
5214 PX = PT*DCOS(ANGDEV-ALPHA)
5215 PY = -PT*DSIN(ANGDEV-ALPHA)
5216 PHIC = DATAN2(PY,PX)
5217 CX = DSIN(THET)*DCOS(PHIC)
5218 CY = DSIN(THET)*DSIN(PHIC)
5221 VECT(1) = XPL3+(ZMAGS-700.)*CX/CZ
5222 VECT(2) = YPL3+(ZMAGS-700.)*CY/CZ
5229 PXZ = PTOT*DSQRT(VECT(4)**2+VECT(6)**2)
5230 RDIP = ASIGN*PXZ / (0.299792458D-3 * B0)
5231 PHI1 = DATAN2(VECT(4),VECT(6))
5232 XC = VECT(1)+RDIP*DCOS(PHI1)
5233 ZC = VECT(3)-RDIP*DSIN(PHI1)
5235 IF (DABS((ZMAGE-ZPLANEP(5))/RDIP).GE.1.D0) RETURN ! Particule boucle dans le champ
5236 XP(5) = XC-ASIGN*DSQRT(RDIP**2-(ZPLANEP(5)-ZC)**2)
5237 YP(5) = VECT(2)+(ZPLANEP(5)-ZMAGS)*VECT(5)/VECT(6)
5238 CX = (ZPLANEP(5)-ZC)/RDIP
5240 CZ = DABS((XP(5)-XC)/RDIP)
5241 THET = DATAN2(DSQRT(CX**2+CY**2),CZ)
5242 AL = THICK/ DCOS(THET)
5243 AP(5) = (0.0136D0/PTOT) * DSQRT(AL) * (1 + 0.038D0*DLOG(AL))
5245 IF (DABS((ZMAGE-ZPLANEP(6))/RDIP).GE.1.D0) RETURN ! Particule boucle dans le champ
5246 XP(6) = XC-ASIGN*DSQRT(RDIP**2-(ZPLANEP(6)-ZC)**2)
5247 YP(6) = VECT(2)+(ZPLANEP(6)-ZMAGS)*VECT(5)/VECT(6)
5248 CX = (ZPLANEP(6)-ZC)/RDIP
5250 CZ = DABS((XP(6)-XC)/RDIP)
5251 THET = DATAN2(DSQRT(CX**2+CY**2),CZ)
5252 AL = THICK/ DCOS(THET)
5253 AP(6) = (0.0136D0/PTOT) * DSQRT(AL) * (1 + 0.038D0*DLOG(AL))
5256 IF (DABS((ZMAGE-ZC)/RDIP).GE.1.D0) RETURN ! Particule boucle dans le champ
5257 VOUT(1) = XC-ASIGN*DSQRT(RDIP**2-(ZMAGE-ZC)**2)
5258 VOUT(2) = VECT(2)+(ZMAGE-ZMAGS)*VECT(5)/VECT(6)
5260 VOUT(4) = (ZMAGE-ZC)/RDIP
5262 VOUT(6) = DABS((VOUT(1)-XC)/RDIP)
5270 THET = DATAN2(DSQRT(VECT(4)**2+VECT(5)**2),VECT(6))
5271 XP(7) = VECT(1)+(ZPLANEP(7)-ZMAGE)*VECT(4)/VECT(6)
5272 YP(7) = VECT(2)+(ZPLANEP(7)-ZMAGE)*VECT(5)/VECT(6)
5273 AL = THICK/ DCOS(THET)
5274 AP(7) = (0.0136D0/PTOT) * DSQRT(AL) * (1 + 0.038D0*DLOG(AL))
5276 XP(8) = XP(7)+(ZPLANEP(8)-ZPLANEP(7))*VECT(4)/VECT(6)
5277 YP(8) = YP(7)+(ZPLANEP(8)-ZPLANEP(7))*VECT(5)/VECT(6)
5278 AL = THICK/ DCOS(THET)
5279 AP(8) = (0.0136D0/PTOT) * DSQRT(AL) * (1 + 0.038D0*DLOG(AL))
5282 XP(9) = XP(8)+(ZPLANEP(9)-ZPLANEP(8))*VECT(4)/VECT(6)
5283 YP(9) = YP(8)+(ZPLANEP(9)-ZPLANEP(8))*VECT(5)/VECT(6)
5284 AL = THICK/ DCOS(THET)
5285 AP(9) = (0.0136D0/PTOT) * DSQRT(AL) * (1 + 0.038D0*DLOG(AL))
5288 XP(10) = XP(9)+(ZPLANEP(10)-ZPLANEP(9))*VECT(4)/VECT(6)
5289 YP(10) = YP(9)+(ZPLANEP(10)-ZPLANEP(9))*VECT(5)/VECT(6)
5290 AL = THICK/ DCOS(THET)
5291 AP(10) = (0.0136D0/PTOT) * DSQRT(AL) * (1 + 0.038D0*DLOG(AL))
5293 ** Matrice de covariance
5296 IF (LPLANE(II).EQ.1) THEN
5301 IF (LPLANE(JJ).EQ.1) THEN
5307 COV(J,I) =COV(J,I) + XPREC**2
5310 * IF (I .EQ. 10 .AND. J .EQ. 10) PRINT *,'10 10 ',COV(J,I)
5313 & + (ZPLANEP(II) + DISTAZ(L))*(ZPLANEP(JJ) + DISTAZ(L))*AMS(L)**2
5317 & + (ZPLANEP(II)-ZPLANEP(K))*(ZPLANEP(JJ)-ZPLANEP(K))*AP(K)**2
5318 * IF (I .EQ. 10 .AND. J .EQ. 10) PRINT *,'10 10 ',COV(J,I)
5321 COVY(J,I) = COV(J,I)
5323 COVY(J,I) = COVY(J,I) - XPREC**2 + YPREC**2
5330 * Inversion des matrices de covariance
5334 CALL DSINV(NPLU, COV, NPLANE, IFAIL)
5335 ** IF (JFAIL.NE.0 .AND. IFAIL .NE. 0) STOP 'ERROR'
5336 IF (IFAIL .NE. 0) STOP 'ERROR'
5338 CALL DSINV(NPLU, COVY, NPLANE, IFAIL)
5339 ** IF (JFAIL.NE.0 .AND. IFAIL .NE. 0) STOP 'ERROR'
5340 IF (IFAIL .NE. 0) STOP 'ERROR'
5341 * PRINT *,' COVARIANCE MATRIX AFTER'
5343 * PRINT *,(COV(J,I),J=1,NPLANE)
5346 ** Calcul de FVAL ou CHI2
5350 IF (LPLANE(II).EQ.1) THEN
5355 IF (LPLANE(JJ).EQ.1) THEN
5358 FVAL = FVAL + COV(J,I)*(XMP(II)-XP(II))*(XMP(JJ)-XP(JJ))
5359 FVAL = FVAL + COVY(J,I)*(YMP(II)-YP(II))
5361 ** IF (JJ.EQ.II) THEN
5362 ** FVAL = FVAL + (XMP(II)-XP(II))*(XMP(JJ)-XP(JJ))/XPREC**2
5363 ** FVAL = FVAL + (YMP(II)-YP(II))
5364 ** & *(YMP(JJ)-YP(JJ))/YPREC**2
5371 print *,' fcnfit pxz tphi talam xvert yvert chi2',
5372 & PXZEA,PHIEA,ALAMEA,
5373 & XEA,YEA,CHI2/FLOAT(2*NPLU-5)
5375 1000 FORMAT(I5,7F12.6)
5384 * Revision 1.5 2000/06/15 07:58:49 morsch
5385 * Code from MUON-dev joined
5387 * Revision 1.4.4.2 2000/04/26 15:48:37 morsch
5388 * Some routines from obsolete algo.F are needed by reco_muon.F and have been
5391 * Revision 1.4.4.1 2000/01/12 16:00:55 morsch
5392 * New version of MUON code
5394 * Revision 1.1.1.1 1995/10/24 10:21:41 cernlib
5398 *CMZ : 3.21/02 29/03/94 15.41.23 by S.Giani
5400 SUBROUTINE RECO_GHELIX (CHARGE, STEP, VECT, VOUT)
5402 C. ******************************************************************
5404 C. * Performs the tracking of one step in a magnetic field *
5405 C. * The trajectory is assumed to be a helix in a constant field *
5406 C. * taken at the mid point of the step. *
5409 C. * STEP =arc length of the step asked *
5410 C. * VECT =input vector (position,direction cos and momentum) *
5411 C. * CHARGE= electric charge of the particle *
5413 C. * VOUT = same as VECT after completion of the step *
5415 C. * ==>Called by : <USER>, GUSWIM *
5416 C. * Author M.Hansroul ********* *
5417 C. * Modified S.Egli, S.V.Levonian *
5418 C. * Modified V.Perevoztchikov
5420 C. ******************************************************************
5422 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
5424 DIMENSION VECT(7),VOUT(7)
5425 DIMENSION XYZ(3),H(4),HXP(3)
5426 PARAMETER (IX=1,IY=2,IZ=3,IPX=4,IPY=5,IPZ=6,IPP=7)
5427 PARAMETER (SIXTH = 1./6.)
5428 PARAMETER (EC=2.9979251E-4)
5430 C. ------------------------------------------------------------------
5432 C units are kgauss,centimeters,gev/c
5434 VOUT(IPP) = VECT(IPP)
5435 IF (CHARGE.EQ.0.) GO TO 10
5436 XYZ(1) = VECT(IX) + 0.5 * STEP * VECT(IPX)
5437 XYZ(2) = VECT(IY) + 0.5 * STEP * VECT(IPY)
5438 XYZ(3) = VECT(IZ) + 0.5 * STEP * VECT(IPZ)
5440 CALL RECO_GUFLD (XYZ, H)
5442 H2XY = H(1)**2 + H(2)**2
5443 H(4) = H(3)**2 + H2XY
5444 IF (H(4).LE.1.E-12) GO TO 10
5445 IF (H2XY.LE.1.E-12*H(4)) THEN
5446 CALL RECO_GHELX3 (CHARGE*H(3), STEP, VECT, VOUT)
5455 HXP(1) = H(2)*VECT(IPZ) - H(3)*VECT(IPY)
5456 HXP(2) = H(3)*VECT(IPX) - H(1)*VECT(IPZ)
5457 HXP(3) = H(1)*VECT(IPY) - H(2)*VECT(IPX)
5459 HP = H(1)*VECT(IPX) + H(2)*VECT(IPY) + H(3)*VECT(IPZ)
5461 RHO = -CHARGE*H(4)/VECT(IPP)
5463 IF (ABS(TET).GT.0.15) THEN
5466 TSINT = (TET-SINT)/TET
5467 COS1T = 2.*(SIN(0.5*TET))**2/TET
5469 TSINT = SIXTH*TET**2
5470 SINTT = (1. - TSINT)
5477 F3 = STEP * TSINT * HP
5480 F6 = TET * COS1T * HP
5482 VOUT(IX) = VECT(IX) + (F1*VECT(IPX) + F2*HXP(1) + F3*H(1))
5483 VOUT(IY) = VECT(IY) + (F1*VECT(IPY) + F2*HXP(2) + F3*H(2))
5484 VOUT(IZ) = VECT(IZ) + (F1*VECT(IPZ) + F2*HXP(3) + F3*H(3))
5486 VOUT(IPX) = VECT(IPX) + (F4*VECT(IPX) + F5*HXP(1) + F6*H(1))
5487 VOUT(IPY) = VECT(IPY) + (F4*VECT(IPY) + F5*HXP(2) + F6*H(2))
5488 VOUT(IPZ) = VECT(IPZ) + (F4*VECT(IPZ) + F5*HXP(3) + F6*H(3))
5494 VOUT(I) = VECT(I) + STEP * VECT(I+3)
5495 VOUT(I+3) = VECT(I+3)
5504 * Revision 1.5 2000/06/15 07:58:49 morsch
5505 * Code from MUON-dev joined
5507 * Revision 1.4.4.2 2000/04/26 15:48:37 morsch
5508 * Some routines from obsolete algo.F are needed by reco_muon.F and have been
5511 * Revision 1.4.4.1 2000/01/12 16:00:55 morsch
5512 * New version of MUON code
5514 * Revision 1.1.1.1 1995/10/24 10:21:41 cernlib
5519 *CMZ : 3.21/02 29/03/94 15.41.23 by S.Giani
5521 SUBROUTINE RECO_GHELX3 (FIELD, STEP, VECT, VOUT)
5523 C. ******************************************************************
5525 C. * Tracking routine in a constant field oriented *
5527 C. * Tracking is performed with a conventional *
5528 C. * helix step method *
5530 C. * ==>Called by : <USER>, GUSWIM *
5531 C. * Authors R.Brun, M.Hansroul ********* *
5532 C * Rewritten V.Perevoztchikov
5534 C. ******************************************************************
5536 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
5538 DIMENSION VECT(7),VOUT(7),HXP(3)
5539 PARAMETER (IX=1,IY=2,IZ=3,IPX=4,IPY=5,IPZ=6,IPP=7)
5540 PARAMETER (SIXTH = 1./6.)
5541 PARAMETER (EC=2.9979251E-4)
5543 C. ------------------------------------------------------------------
5545 C units are kgauss,centimeters,gev/c
5547 VOUT(IPP) = VECT(IPP)
5550 HXP(1) = - VECT(IPY)
5551 HXP(2) = + VECT(IPX)
5557 IF (ABS(TET).GT.0.15) THEN
5560 TSINT = (TET-SINT)/TET
5561 COS1T = 2.*(SIN(0.5*TET))**2/TET
5563 TSINT = SIXTH*TET**2
5564 SINTT = (1. - TSINT)
5571 F3 = STEP * TSINT * HP
5574 F6 = TET * COS1T * HP
5576 VOUT(IX) = VECT(IX) + (F1*VECT(IPX) + F2*HXP(1))
5577 VOUT(IY) = VECT(IY) + (F1*VECT(IPY) + F2*HXP(2))
5578 VOUT(IZ) = VECT(IZ) + (F1*VECT(IPZ) + F3)
5580 VOUT(IPX) = VECT(IPX) + (F4*VECT(IPX) + F5*HXP(1))
5581 VOUT(IPY) = VECT(IPY) + (F4*VECT(IPY) + F5*HXP(2))
5582 VOUT(IPZ) = VECT(IPZ) + (F4*VECT(IPZ) + F6)
5587 ************************************************************************
5588 DOUBLE PRECISION FUNCTION RECOCHI2 (MPOS,MANG,XM,YM,ALAMM,APHIM,
5589 & XC,YC,ALAMC,APHIC,PTOT,IZMES,NPLPL)
5591 C. ******************************************************************
5592 C. * Calculate chi2 taking into account MSC *
5594 C. ******************************************************************
5596 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
5597 PARAMETER(NBSTATION=5,NPLANE=10)
5599 COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
5601 COMMON/PARAM/ZPLANEP(NPLANE),THICK,XPREC,YPREC,B0,BL3,ZMAGS,
5602 & ZMAGE,ZABS,XMAG,ZBP1,ZBP2,CONST
5604 DIMENSION MPOS(NBSTATION),MANG(NBSTATION),
5605 & XM(NBSTATION),YM(NBSTATION),ALAMM(NBSTATION),APHIM(NBSTATION),
5606 & XC(NBSTATION),YC(NBSTATION),ALAMC(NBSTATION),APHIC(NBSTATION),
5607 & ALAMP(NBSTATION),APHIP(NBSTATION)
5609 DIMENSION AP(NPLANE),IZMES(NBSTATION),IPLANE(NPLANE)
5610 DIMENSION COVX(NPLANE,NPLANE),COVY(NPLANE,NPLANE)
5611 DIMENSION XM1(NPLANE),YM1(NPLANE),XC1(NPLANE),YC1(NPLANE)
5619 AL = THICK/(COS(ALAMC(IZ))*COS(APHIC(IZ)))
5620 AP(ICH) = (0.0136D0/PTOT)*DSQRT(AL)*(1.+0.038D0*DLOG(AL))
5621 IF (MPOS(IZ).EQ.1.AND.IZMES(IZ).EQ.1) THEN
5631 AP(ICH) = (0.0136D0/PTOT)*DSQRT(AL)*(1.+0.038D0*DLOG(AL))
5632 IF (MPOS(IZ).EQ.1.AND.IZMES(IZ).EQ.2) THEN
5637 IF (MANG(IZ).EQ.1) THEN
5639 XM1(ICH) = XM(IZ) - DZ_PL(IZ) * TAN(APHIM(IZ))
5640 YM1(ICH) = YM(IZ) + DZ_PL(IZ)/COS(APHIM(IZ))*TAN(ALAMM(IZ))
5642 XC1(ICH) = XC(IZ) - DZ_PL(IZ) * TAN(APHIC(IZ))
5643 YC1(ICH) = YC(IZ) + DZ_PL(IZ)/COS(APHIC(IZ))*TAN(ALAMC(IZ))
5644 IF (IPLANE(ICH).EQ.1) NPLPL = NPLPL+1
5649 ** Matrice de covariance X et Y
5653 IF (IPLANE(II).EQ.1) THEN
5658 IF (IPLANE(JJ).EQ.1) THEN
5664 COVX(J,I) =COVX(J,I) + XPREC**2
5668 COVX(J,I) = COVX(J,I)
5669 & + (-ZPLANEP(II)+ZPLANEP(K))*
5670 & (-ZPLANEP(JJ)+ZPLANEP(K))*AP(K)**2
5674 COVY(J,I) = COVX(J,I)
5676 COVY(J,I) = COVY(J,I) - XPREC**2 + YPREC**2
5686 * Inversion des matrices de covariance
5688 CALL DSINV(NPLUP, COVX, NPLANE, IFAIL)
5689 ** IF (JFAIL.NE.0 .AND. IFAIL .NE. 0) STOP 'ERROR'
5690 IF (IFAIL .NE. 0) STOP 'RECOCHI2 ERROR COVX'
5693 CALL DSINV(NPLUP, COVY, NPLANE, IFAIL)
5694 ** IF (JFAIL.NE.0 .AND. IFAIL .NE. 0) STOP 'ERROR'
5695 IF (IFAIL .NE. 0) STOP 'RECOCHI2 ERROR COVY'
5697 ** Calcul de FVAL ou CHI2
5702 IF (IPLANE(II).EQ.1) THEN
5705 ** print *,' II=',ii,' XM1,YM1 =',xm1(ii),ym1(ii),
5706 ** & ' XC,YC =',xc1(ii),yc1(ii)
5708 IF (IPLANE(JJ).EQ.1) THEN
5710 FVAL = FVAL + COVX(J,I)*(XM1(II)-XC1(II))
5711 & *(XM1(JJ)-XC1(JJ))
5712 FVAL = FVAL + COVY(J,I)*(YM1(II)-YC1(II))
5713 & *(YM1(JJ)-YC1(JJ))
5720 ** print *,' recochi2 =',recochi2
5726 ************************************************************************
5727 SUBROUTINE RECO_SELECT(ISEL)
5728 ************************************************************************ * ISEL(I) = 1 if track number I is OK, ISEL(I) = 0 otherwise
5729 ************************************************************************
5730 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
5732 PARAMETER (NPLANE=10,NBCHAMBER=10,NTRMAX=500)
5734 COMMON/DEBEVT/IDEBUG
5736 COMMON/PARAM/ZPLANEP(NPLANE),THICK,XPREC,YPREC,B0,BL3,ZMAGS,
5737 & ZMAGE,ZABS,XMAG,ZBP1,ZBP2,CONST
5739 COMMON/TRACKFOUT/IEVOUT,NTREVT,JJOUT(NBCHAMBER,NTRMAX),
5740 & ISTAT(NTRMAX),PXZOUT(NTRMAX),TPHIOUT(NTRMAX),
5741 & TALAMOUT(NTRMAX),XVERTOUT(NTRMAX),YVERTOUT(NTRMAX),
5743 & XMESOUT(NBCHAMBER,NTRMAX),YMESOUT(NBCHAMBER,NTRMAX)
5744 & ,PXVOUT(NTRMAX),PYVOUT(NTRMAX),PZVOUT(NTRMAX)
5746 DIMENSION ISEL(NTRMAX)
5756 IF (JJ1.EQ.0) ICH1 = 10
5757 X1 = XMESOUT(ICH1,I)
5758 Y1 = YMESOUT(ICH1,I)
5763 IF (JJ2.EQ.0) ICH2 = 10
5764 X2 = XMESOUT(ICH2,J)
5765 Y2 = YMESOUT(ICH2,J)
5766 DIST = SQRT(((X2-X1)/(10.*XPREC))**2
5767 & +((Y2-Y1)/(10.*YPREC))**2)
5768 IF (DIST.LT.2.) THEN
5771 IF (CHI22.LT.CHI21) THEN
5773 IF (IDEBUG.EQ.2) THEN
5774 PRINT *,' RECO_SELECT I,ISEL= ',I,ISEL(I)
5778 IF (IDEBUG.EQ.2) THEN
5779 PRINT *,' RECO_SELECT J,ISEL= ',J,ISEL(J)
5791 ***************************************************
5792 SUBROUTINE reconstmuon2(IFIT,IDEBUGC,NEV)
5793 ***************************************************
5794 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
5796 PARAMETER(NPLANE=10,MAXIDG=28,NTRMAX=500)
5798 REAL*4 PXR,PYR,PZR,ZVR,CHI2R,PXV,PYV,PZV
5800 COMMON/TRACKFI/EFF,EFF1,EFF2,XPREC,YPREC,PHIPREC,ALAMPREC,
5801 & HCUT,LBKG,SIGCUT,ALPHATOP,HTOP
5803 COMMON/PAWCR4/IEVR,NTRACKR,ISTATR(NTRMAX),ISIGNR(NTRMAX),
5804 & PXR(NTRMAX),PYR(NTRMAX),PZR(NTRMAX),ZVR(NTRMAX),
5805 & CHI2R(NTRMAX),PXV(NTRMAX),PYV(NTRMAX),PZV(NTRMAX)
5807 COMMON/MEAS/LPLANE(NPLANE),XMP(NPLANE),YMP(NPLANE)
5809 COMMON/FIT/NHITTOT1,izch(maxidg),xgeant(maxidg),
5812 COMMON/PRECCUT/PCUT,PTCUT,CHI2CUT
5817 CALL trackf_read_fit(IEVR,NEV,NHITTOT1,IZCH,XGEANT,YGEANT)
5818 print*,'nhittot1 ',nhittot1
5821 * print*,'x=',xgeant(i),' y=',ygeant(i),' ch=',izch(i)
5824 if (nhittot1.ne.20) goto 55
5827 do ntr=1,2 ! loop over tracks
5828 do nhit=nhit1,nhit2 ! loop over hits
5832 CALL RANNOR(RN1,RN2) ! CCC
5833 * xmp(ich)=xgeant(nhit)
5834 * ymp(ich)=ygeant(nhit)
5835 xmp(ich)=xgeant(nhit) + RN1 * XPREC
5836 ymp(ich)=ygeant(nhit) + RN2 * YPREC
5843 CALL CHFNT(IEVR,NTRACK,ISTATR,ISIGNR,
5844 & PXR,PYR,PZR,ZVR,CHI2R,PXV,PYV,PZV)
5850 *********************************
5851 subroutine fit_trace(ntr)
5852 *********************************
5854 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
5856 PARAMETER(NPLANE=10,NTRACK=2,NTRMAX=500)
5858 COMMON/FCNOUT/PXZEA,ALAMEA,PHIEA,XEA,YEA,NPLU,CHI2
5860 REAL*4 PXR,PYR,PZR,ZVR,CHI2R,PXV,PYV,PZV
5861 COMMON/PAWCR4/IEVR,NTRACKR,ISTATR(NTRMAX),ISIGNR(NTRMAX),
5862 & PXR(NTRMAX),PYR(NTRMAX),PZR(NTRMAX),ZVR(NTRMAX),
5863 & CHI2R(NTRMAX),PXV(NTRMAX),PYV(NTRMAX),PZV(NTRMAX)
5865 COMMON/MEAS/LPLANE(NPLANE),XMP(NPLANE),YMP(NPLANE)
5867 COMMON/PARAM/ZPLANEP(NPLANE),THICK,XPREC,YPREC,B0,BL3,ZMAGS,
5868 & ZMAGE,ZABS,XMAG,ZBP1,ZBP2,CONST
5870 COMMON/PRECCUT/PCUT,PTCUT,CHI2CUT
5872 dimension PXZOUT(NTRMAX),TPHIOUT(NTRMAX),TALAMOUT(NTRMAX),
5886 PHIAV = DATAN2((X2-X1),(ZPLANEP(IPL2)-ZPLANEP(IPL1)))
5887 PHIAP = DATAN2((X4-X3),(ZPLANEP(IPL4)-ZPLANEP(IPL3)))
5889 DPHI = (PHIAP-PHIAV)
5891 IF (DPHI.LT.0.) ASIGN = -1. ! CCC
5892 PXZ = CONST/DABS(DPHI)
5895 IF (PXZ.LT.PCUT) GO TO 66
5897 PXZINVI = ASIGN/PXZ ! CCC
5899 ALAMI = DATAN2((Y2-Y1),DSQRT((X2-X1)**2
5900 & +(ZPLANEP(IPL2)-ZPLANEP(IPL1))**2))
5904 * print *,' avant prec_fit pxzi phii alami x y',1./ PXZINVI,
5905 * & PHII, ALAMI ,XVR,YVR
5906 * PRINT *,' X1 X2 X3 X4',X1,X2,X3,X4
5907 * PRINT *,' Z1 Z2 Z3 Z4',ZPLANEP(IPL1),ZPLANEP(IPL2),
5908 * & ZPLANEP(IPL3),ZPLANEP(IPL4)
5909 * PRINT *,' CONST= ',CONST
5911 * * Fit des traces apres l'absorbeur
5912 CALL PREC_FIT (PXZINVI,PHII,ALAMI,XVR,YVR,
5913 & PXZINVF,PHIF,ALAMF,XVERTF,YVERTF,EPXZINV,EPHI,EALAM,
5916 * * Correction de Branson
5917 CALL BRANSON(PXZEA,PHIEA,ALAMEA,XEA,YEA)
5920 PX1 = PXZ1*DSIN(PHIEA)
5921 PY1 = PXZ1*DTAN(ALAMEA)
5922 PT1 = DSQRT(PX1**2+PY1**2)
5925 * print*,'ptcut=',ptcut
5926 * print*,'chi2=',CHI2/FLOAT(2*NPLU-5)
5927 * print*,'chi2cut=',CHI2CUT
5930 IF (PT1.LT.PTCUT) GO TO 66
5932 IF ((CHI2/FLOAT(2*NPLU-5)).GT.CHI2CUT) GO TO 66
5935 TPHIOUT(NTR) = DTAN(PHIEA)
5936 TALAMOUT(NTR) = DTAN(ALAMEA)
5937 CHI2OUT(NTR) = CHI2/FLOAT(2*NPLU-5)
5940 IF (PXZOUT(NTR).LT.0.) ISIGNR(NTR) = -1
5941 PXZ = ABS(PXZOUT(NTR))
5942 PHI = ATAN(TPHIOUT(NTR))
5943 ALAM = ATAN(TALAMOUT(NTR))
5944 PYR(NTR) = PXZ*SIN(PHI)
5945 PXR(NTR) = PXZ*TAN(ALAM)
5946 PZR(NTR) = PXZ*COS(PHI)
5953 SUBROUTINE SORTZV (A,INDEX,N1,MODE,NWAY,NSORT)
5955 C CERN PROGLIB# M101 SORTZV .VERSION KERNFOR 3.15 820113
5958 DIMENSION A(N1),INDEX(N1)
5963 IF (NSORT.NE.0) GO TO 2
5967 2 IF (N.EQ.1) RETURN
5969 10 CALL SORTTI (A,INDEX,N)
5972 20 CALL SORTTC(A,INDEX,N)
5975 30 CALL SORTTF (A,INDEX,N)
5977 40 IF (NWAY.EQ.0) GO TO 50
5987 SUBROUTINE SORTTI (A,INDEX,N1)
5990 DIMENSION A(N1),INDEX(N1)
6000 IF (AI.LE.A (I22)) GO TO 3
6006 INDEX (N) = INDEX (1)
6012 IF (I2.LE.N) I22= INDEX(I2)
6014 7 I222 = INDEX (I2+1)
6015 IF (A(I22)-A(I222)) 8,9,9
6018 9 IF (AI-A(I22)) 10,11,11
6028 * ========================================
6029 SUBROUTINE SORTTC (A,INDEX,N1)
6032 DIMENSION A(N1),INDEX(N1)
6042 IF(ICMPCH(AI,A(I22)))3,3,21
6048 INDEX (N) = INDEX (1)
6054 IF (I2.LE.N) I22= INDEX(I2)
6056 7 I222 = INDEX (I2+1)
6057 IF (ICMPCH(A(I22),A(I222))) 8,9,9
6060 9 IF (ICMPCH(AI,A(I22))) 10,11,11
6069 * ========================================
6070 FUNCTION ICMPCH(IC1,IC2)
6071 C FUNCTION TO COMPARE TWO 4 CHARACTER EBCDIC STRINGS - IC1,IC2
6072 C ICMPCH=-1 IF HEX VALUE OF IC1 IS LESS THAN IC2
6073 C ICMPCH=0 IF HEX VALUES OF IC1 AND IC2 ARE THE SAME
6074 C ICMPCH=+1 IF HEX VALUES OF IC1 IS GREATER THAN IC2
6077 IF(I1.GE.0.AND.I2.GE.0)GOTO 40
6083 40 IF(I1-I2)60,70,80
6092 SUBROUTINE SORTTF (A,INDEX,N1)
6094 DIMENSION A(N1),INDEX(N1)
6104 IF (AI.LE.A (I22)) GO TO 3
6110 INDEX (N) = INDEX (1)
6116 IF (I2.LE.N) I22= INDEX(I2)
6118 7 I222 = INDEX (I2+1)
6119 IF (A(I22)-A(I222)) 8,9,9
6122 9 IF (AI-A(I22)) 10,11,11