MUON/reco_muon.F

   1 *  25/5/99
   2 **  Authors J.P. Cussonneau & P. Lautridou
   3
   4 ************** commentaires **********************
   5
   6 * bon muon = tous les hits de la trace (pas forcement 28) proviennent de muons
   7 * ghost = bonne trace dans laquelle certains hits proviennent délectrons, de muons ou sont ambigus
   8
   9 * ALPHATOP
  10 * alpha to P
  11
  12 * EFF
  13 * efficacite
  14
  15 * HHIT
  16 * h au vertex
  17
  18 * HTOP
  19 * h to P
  20
  21 * INDEXMAX
  22 * nombre de candidats a ordonner (a partir de l'ímpulsion)
  23
  24 * INDEXTAB
  25 * pour recuperer les candidats
  26
  27 * ISTAT
  28 * =1 si bon muon, =2 si ghost, =0 autrement
  29
  30 * ITCHECK
  31 * =1 si bonne trace, =0 autrement
  32
  33 * IT_LIST
  34 * permet a partir du numero de la trace de retrouver la numero du hit
  35
  36 * IT_NP
  37 * compte le nombre de plans touches par trace
  38
  39 * ITRACK
  40 * permet de retrouver le numero de la trace a partir du numero de hit
  41
  42 * ITTROUGH
  43 * pour une trace et une station donnee, dit si la trace est passee dans la chambre
  44
  45 * IVERTEX
  46 * =0 si point (0,0) du vertex impose, sinon coordonnes libres (pour le fit)
  47
  48 * JCAN
  49 * numero du hit pour une station et un candidat
  50
  51 * JCANTYP
  52 * nombre de hits par trace au niveau des stations 4 et 5 (3 ou 4)
  53
  54 * JJOUT
  55 * numero de hit associe a une chambre et a une trace
  56
  57 * NMUONALL
  58 * nombre de bons muons trouves
  59
  60 * NERR
  61 * = nombre de fois ou lón ná pas trouve la bonne trace dans la station
  62
  63 * NERRALL
  64 * nombre de cas ou pas de hit trouve par station
  65
  66 * NGHOSTALL
  67 * nombre de fantomes trouves
  68
  69 * NRES
  70 * nombre de resonances dans lácceptance
  71
  72 * NRESF
  73 * nombre de resonances trouvees
  74
  75 * NTRACFALL
  76 * nombre de traces totales trouvees
  77
  78 * NTRMUALL
  79 * nombre total de muons dans lácceptance
  80
  81
  82 ****************************************************************
  83       SUBROUTINE reconstmuon(IFIT,IDEBUGC,NEV,IDRES,IREADGEANT)
  84 ****************************************************************
  85       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
  86
  87       COMMON/DEBEVT/IDEBUG
  88
  89       IDEBUG=IDEBUGC
  90
  91 ** Read events
  92       CALL RECO_READEVT(NEV,IDRES,IREADGEANT)
  93
  94
  95 ** Trackfinding
  96       CALL RECO_TRACKF(IDRES,IREADGEANT)
  97
  98 ** Precision fit
  99       IF (IFIT.EQ.1) THEN
 100          CALL RECO_PRECISION
 101       ENDIF
 102
 103 ** Calculate
 104       CALL RECO_SUM
 105
 106       END
 107
 108 ********************************************
 109       BLOCK DATA
 110 ********************************************
 111       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
 112       PARAMETER(NBSTATION=5)
 113 * --
 114       COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
 115
 116       DATA DZ_PL/8.,8.,24.3,8.,8./
 117
 118 **      DATA DZ_PL/8.,8.,8.,8.,8./
 119 *      DATA DZ_PL/20.,20.,20.,20.,20./
 120 *      DATA DZ_PL/20.,20.,24.3,20.,20./
 121 *      DATA DZ_PL/20.,40.,40.,40.,40./
 122 *      DATA DZ_PL/20.,50.,50.,50.,50./
 123 *      DATA DZ_PL/20.,55.,55.,55.,55./
 124 *      DATA DZ_PL/20.,60.,60.,60.,60./
 125 *      DATA DZ_PL/20.,80.,80.,80.,80./
 126 *      DATA DZ_PL/20.,100.,100.,100.,100./
 127
 128       DATA ZPLANE/-511.,-686.0,-962.7,-1245.,-1445./ ! dstation=8cm
 129
 130 *      DATA ZPLANE/-507.5,-682.5,-967.5,-1241.5,-1441.5/ ! dstation=20cm CCC
 131 *      DATA ZPLANE/-505.,-680.,-965.,-1239.,-1439./ ! dstation=20cm
 132 **      DATA ZPLANE/-511.,-686.0,-971.,-1245.,-1445./ ! dstation=8cm
 133 *      DATA ZPLANE/-511.,-686.0,-962.7,-1245.,-1445./ ! dstation=8cm
 134 *      DATA ZPLANE/-505.,-680.,-962.7,-1239.,-1439./   ! dstation=20cm
 135 *      DATA ZPLANE/-505.,-670.,-954.,-1229.,-1429./ ! dstation=40cm
 136 *      DATA ZPLANE/-505.,-665.,-949.,-1224.,-1424./ ! dstation=50cm
 137 *      DATA ZPLANE/-505.,-663.,-947.,-1222.,-1422./ ! dstation=50cm
 138 *      DATA ZPLANE/-505.,-660.,-944.,-1219.,-1419./ ! dstation=60cm
 139 *      DATA ZPLANE/-505.,-650.,-935.,-1209.,-1409./ ! dstation=80cm
 140 *      DATA ZPLANE/-505.,-640.,-925.,-1199.,-1399./ ! dstation=100cm
 141
 142 **      DATA ZPLANE/-507.,-682.0,-950.7,-1241.,-1441./
 143
 144 **      DATA ZCOIL,ZMAGEND/-825.0,-1125./  ! Constant field 3 Tm
 145       DATA ZCOIL,ZMAGEND/-805.0,-1233./ ! CCC magn. field map M.B
 146 **      DATA ZCOIL,ZMAGEND/-795.1,-1242.9/ ! CCC magn. field map M.B
 147
 148 * --
 149       END
 150
 151 ********************************************
 152       SUBROUTINE cutpxz(spxzcut)
 153 ********************************************
 154       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
 155
 156       COMMON/ACUTPXZ/ACUTPXZ
 157
 158       ACUTPXZ = SPXZCUT
 159
 160       END
 161
 162 ********************************************
 163       SUBROUTINE sigmacut(ssigcut)
 164 ********************************************
 165       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
 166
 167       COMMON/TRACKFI/EFF,EFF1,EFF2,XPREC,YPREC,PHIPREC,ALAMPREC,
 168      &                 HCUT,LBKG,SIGCUT,ALPHATOP,HTOP
 169
 170 **      SIGCUT = SSIGCUT
 171       SIGCUT=5.  ! CCC
 172
 173       END
 174
 175 ********************************************
 176       SUBROUTINE xpreci(sxprec)
 177 ********************************************
 178       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
 179
 180       COMMON/TRACKFI/EFF,EFF1,EFF2,XPREC,YPREC,PHIPREC,ALAMPREC,
 181      &                 HCUT,LBKG,SIGCUT,ALPHATOP,HTOP
 182
 183 **      XPREC = SXPREC
 184       XPREC = 0.01 ! CCC
 185
 186       END
 187
 188 ********************************************
 189       SUBROUTINE ypreci(syprec)
 190 ********************************************
 191       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
 192
 193       COMMON/TRACKFI/EFF,EFF1,EFF2,XPREC,YPREC,PHIPREC,ALAMPREC,
 194      &                 HCUT,LBKG,SIGCUT,ALPHATOP,HTOP
 195
 196       YPREC = SYPREC
 197 **      YPREC = 0.2 ! CCC
 198
 199       END
 200
 201 *************************************************************************
 202       SUBROUTINE RECO_INIT(seff,sb0,sbl3)
 203 *************************************************************************
 204       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
 205
 206       CALL TRACKF_INIT(seff,sb0,sbl3)
 207
 208       CALL PREC_INIT
 209
 210       RETURN
 211       END
 212
 213 *************************************************************************
 214       SUBROUTINE TRACKF_INIT(seff,sb0,sbl3)
 215 *************************************************************************
 216       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
 217 **
 218       PARAMETER(NBSTATION=5,NTRMAX=500)
 219 **
 220       COMMON/REVENT/IEVBKGI,NBKGMAX,MAXUPSEV
 221 **
 222       COMMON/MAGNET/BL3,B0
 223 **
 224       COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
 225 **
 226       COMMON/FILED/FILERES,FILEBKG,FILEOUT,FILEMIN
 227 **
 228       COMMON/TRACKFI/EFF,EFF1,EFF2,XPREC,YPREC,PHIPREC,ALAMPREC,
 229      &                 HCUT,LBKG,SIGCUT,ALPHATOP,HTOP
 230 **
 231       COMMON/TRACKSUM/NRES(5),NRESF,NTRMUALL,NMUONALL,NGHOSTALL,
 232      &     NTRACKFALL,NERRALL(NBSTATION),IR
 233 **
 234       COMMON/ACUTPXZ/ACUTPXZ
 235 **
 236       COMMON/DEBEVT/IDEBUG
 237 **
 238       CALL HIST_CREATE
 239 *
 240       EFF = SEFF
 241       B0 = SB0
 242       BL3 = SBL3
 243
 244       PXZCUT = ACUTPXZ
 245       AMAGLEN = ZMAGEND-ZCOIL
 246       ZM = AMAGLEN/2.+ZCOIL
 247       ALPHATOP = 0.01*0.3*B0*ABS(AMAGLEN)
 248       HTOP = ALPHATOP*ZM
 249       HCUT = ABS(HTOP)/PXZCUT
 250
 251       print*,'TRACK_INIT hcut= ',hcut
 252       print*,'TRACK_INIT eff = ',eff
 253       print*,'TRACK_INIT b0 = ',b0
 254       print*,'TRACK_INIT bl3 = ',bl3
 255       print*,'TRACK_INIT sigmacut = ',sigcut
 256       print*,'TRACK_INIT cutpxz = ',pxzcut
 257       print*,'TRACK_INIT xprec = ',xprec
 258       print*,'TRACK_INIT yprec = ',yprec
 259
 260       EFF2 = EFF**2               ! PROBA. DEUX CHAMBRES TOUCHES
 261       EFF1 = EFF2+2.*EFF*(1.-EFF) ! PROBA. AU MOINS UNE CHAMBRE TOUCHE
 262 ** Used only for stations 4 & 5
 263       PHIPREC   = SQRT(2.)*XPREC/DZ_PL(5) ! PHI = (OZ , PROJ. DANS XOZ)
 264       ALAMPREC  = SQRT(2.)*YPREC/DZ_PL(5) ! LAM = (OM , PROJ. DANS XOZ)
 265
 266       DO I = 1,5
 267          NRES(I) = 0
 268       ENDDO
 269       NRESF = 0
 270       NTRMUALL = 0
 271       NMUONALL = 0
 272       NGHOSTALL = 0
 273       NTRACKFALL = 0
 274       DO I = 1,NBSTATION
 275          NERRALL(I) = 0
 276       ENDDO
 277       IR = 0
 278 *
 279       RETURN
 280       END
 281
 282 *************************************************************************
 283       SUBROUTINE PREC_INIT
 284 *************************************************************************
 285 *
 286 *
 287 *************************************************************************
 288
 289       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
 290 *
 291       PARAMETER(NPLANE=10,NBSTATION=5)
 292 *
 293       COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
 294 *
 295       COMMON/PARAM/ZPLANEP(NPLANE),THICK,XPREC,YPREC,B0,BL3,ZMAGS,
 296      &             ZMAGE,ZABS,XMAG,ZBP1,ZBP2,CONST
 297 *
 298       COMMON/PRECCUT/PCUT,PTCUT,CHI2CUT
 299 *
 300       COMMON/PRECSUM/NRESF1,NMUONALL1,NGHOSTALL1,NTRACKFALL1
 301 *
 302       COMMON/DEBEVT/IDEBUG
 303 *
 304
 305       DATA THICK/0.03D0/   ! X/X0=3%
 306 **      DATA THICK/0.02D0/   ! X/X0=2% chambre
 307
 308 **      DATA B0,BL3/10.,2.0/ ! Champ magnetique dans le dipole et dans L3 en kgauss
 309       DATA B0,BL3/7.,2.0/ ! Magnetic field in the dipole & L3 in kgauss
 310       DATA ZMAGS/805.0D0/,ZMAGE/1233.0D0/,ZABS/503.D0/ ! CCC not used when
 311                                                        ! magn. field map
 312 **      DATA ZMAGS/825.0D0/,ZMAGE/1125.0D0/,ZABS/503.D0/
 313       DATA XMAG/190.0D0/
 314
 315 **      DATA XPREC/0.0100D0/,YPREC/0.144337D0/ ! CCC
 316       DATA XPREC/0.0100D0/,YPREC/0.2D0/
 317
 318 * Input parameters
 319       CONST = 0.299792458D-3*B0*(ZMAGE-ZMAGS)
 320       J = 0
 321       DO I = 1,5
 322          J = J+1
 323          ZPLANEP(J) = -ZPLANE(I)
 324          J = J+1
 325          ZPLANEP(J) = -ZPLANE(I)+DZ_PL(I)
 326       ENDDO
 327
 328       PCUT = 3.        ! Coupure en PXZ muon (GeV/c)
 329       PTCUT = 0.       ! Coupure en Pt muon (GeV/c)
 330
 331       CHI2CUT = 1.E4  ! Coupure sur le CHI2 du fit
 332
 333       X01 = 18.8     ! C (cm)
 334       X02 = 10.397   ! Concrete (cm)
 335       X03 = 0.56     ! Plomb (cm)
 336       X04 = 47.26    ! Polyethylene (cm)
 337       X05 = 0.35     ! W (cm)
 338
 339 ** Calcul des parametres pour la correction de Branson de l'absorbeur
 340       ANBP = (315.**3-90.**3)/X01 +(467.**3-315.**3)/X02+
 341      &       (472.**3-467.**3)/X03+(477.**3-472.**3)/X04+
 342      &       (482.**3-477.**3)/X03+(487.**3-482.**3)/X04+
 343      &       (492.**3-487.**3)/X03+(497.**3-492.**3)/X04+
 344      &       (502.**3-497.**3)/X03
 345       ADBP = (315.**2-90.**2)/X01 +(467.**2-315.**2)/X02+
 346      &       (472.**2-467.**2)/X03+(477.**2-472.**2)/X04+
 347      &       (482.**2-477.**2)/X03+(487.**2-482.**2)/X04+
 348      &       (492.**2-487.**2)/X03+(497.**2-492.**2)/X04+
 349      &       (502.**2-497.**2)/X03
 350       ZBP1 = 2./3.*ANBP/ADBP
 351       ANBP = (315.**3-90.**3)/X01 +(467.**3-315.**3)/X02+
 352      &       (503.**3-467.**3)/X05
 353       ADBP = (315.**2-90.**2)/X01 +(467.**2-315.**2)/X02+
 354      &       (503.**2-467.**2)/X05
 355       ZBP2 = 2./3.*ANBP/ADBP
 356 *
 357       IF (IDEBUG.GE.1) THEN
 358          PRINT *,' PREC_INIT B0 (kgauss)',B0,' BL3 (kgauss)',BL3
 359          PRINT *,' PREC_INIT ZMAGE (cm)',ZMAGE,' ZMAGS (cm)',ZMAGS,
 360      &        ' XMAG (cm)',XMAG
 361          PRINT *,' PREC_INIT ZABS (cm)',ZABS,' ZBP1 (cm)',ZBP1,
 362      &        ' ZBP2 (cm)',ZBP2
 363          PRINT *,' PREC_INIT Radiation length absorber X01 (cm)',X01,
 364      &        ' X02 (cm)',X02
 365          PRINT *,' PREC_INIT X03(cm)',X03,' X04 (cm)',X04,' X05 (cm)',
 366      &           X05
 367          PRINT *,' PREC_INIT Radiation length chamber THICK (%)',
 368      &           THICK*100.
 369          PRINT *,' PREC_INIT XPREC (cm)',XPREC,' YPREC (cm)',YPREC
 370          PRINT *,' PREC_INIT Coupure en Pxz (GeV/c): ',PCUT
 371          PRINT *,' PREC_INIT Coupure en Pt (GeV/c): ',PTCUT
 372          PRINT *,' PREC_INIT Coupure en CHI2 : ',CHI2CUT
 373       ENDIF
 374
 375 *      PAUSE
 376
 377       NRESF1 = 0
 378       NMUONALL1 = 0
 379       NGHOSTALL1 = 0
 380       NTRACKFALL1 = 0
 381
 382 *  Magnetic Field Map GC
 383 **      OPEN (UNIT=40,FILE='data/field02.dat',
 384 **     &      STATUS='UNKNOWN')
 385
 386       CALL INITFIELD
 387
 388 **      CLOSE(40)
 389
 390       RETURN
 391       END
 392
 393 *************************************************************************
 394       SUBROUTINE RECO_READEVT(NEV,IDRES,IREADGEANT)
 395 *************************************************************************
 396 *
 397 *
 398 *************************************************************************
 399       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
 400
 401       PARAMETER(NTRMAX=500)
 402
 403       PARAMETER (NBSTATION=5,MAXIDG=20000,MAXHITTOT=20000,
 404      &           MAXHITCH=10000,MAXHIT=1000,NBCHAMBER=10)
 405
 406       COMMON/RHITG/ITYPG(MAXIDG),XTRG(MAXIDG),YTRG(MAXIDG),
 407      &             PTOTG(MAXIDG),IDG(MAXIDG),IZCH(MAXIDG),
 408      &             PVERT1G(MAXIDG),PVERT2G(MAXIDG),PVERT3G(MAXIDG),
 409      &     ZVERTG(MAXIDG),NHITTOT1,CX(MAXIDG),CY(MAXIDG),CZ(MAXIDG),
 410      &     XGEANT(MAXIDG),YGEANT(MAXIDG),CLSIZE1(MAXIDG),CLSIZE2(MAXIDG)
 411
 412       DIMENSION TYPG(MAXIDG),ZCH(MAXIDG)
 413
 414       REAL*4 R1,R2
 415       DATA R1,R2/0.,1./
 416
 417       IF (IREADGEANT.eq.1) THEN  ! GEANT hits
 418
 419          CALL TRACKF_READ_GEANT(ITYPG,XTRG,YTRG,PTOTG,IDG,IZCH,PVERT1G,
 420      &        PVERT2G,PVERT3G,ZVERTG,NHITTOT1,CX,CY,CZ,IEVR,NEV,
 421      &        XGEANT,YGEANT,CLSIZE1,CLSIZE2)
 422       ELSE
 423          CALL TRACKF_READ_SPOINT(ITYPG,XTRG,YTRG,PTOTG,IDG,IZCH,PVERT1G,
 424      &        PVERT2G,PVERT3G,ZVERTG,NHITTOT1,CX,CY,CZ,IEVR,NEV,
 425      &        XGEANT,YGEANT,CLSIZE1,CLSIZE2)
 426       ENDIF
 427
 428       do i=1,NHITTOT1
 429          TYPG(i)=ITYPG(i)
 430          call chfill(100,sngl(typg(i)),R1,R2)
 431          call chfill(101,sngl(ygeant(i)),R1,R2)
 432          call chfill(102,sngl(xgeant(i)),R1,R2)
 433          ZCH(i)=IZCH(i)
 434          call chfill(103,sngl(zch(i)),R1,R2)
 435          call chfill(104,sngl(ptotg(i)),R1,R2)
 436          call chfill(105,sngl(pvert2g(i)),R1,R2)
 437          call chfill(106,sngl(pvert1g(i)),R1,R2)
 438          call chfill(107,sngl(pvert3g(i)),R1,R2)
 439          call chfill(108,sngl(zvertg(i)),R1,R2)
 440          call chfill(109,sngl(ytrg(i)),R1,R2)
 441          call chfill(110,sngl(xtrg(i)),R1,R2)
 442          dx=xgeant(i)-xtrg(i)
 443          dy=ygeant(i)-ytrg(i)
 444          call chfill(111,sngl(dy),R1,R2)
 445          call chfill(112,sngl(dx),R1,R2)
 446          call chfill(119+int(zch(i)),sngl(dy),R1,R2)
 447          call chfill(129+int(zch(i)),sngl(dx),R1,R2)
 448       enddo
 449
 450       do i=1,NHITTOT1
 451          CALL CHFILL (999,SNGL(PTOTG(I)),R1,R2)
 452       enddo
 453
 454       CALL TRACKF_STAT(IDRES,IREADGEANT)
 455
 456       RETURN
 457       END
 458
 459 *************************************************************************
 460       SUBROUTINE TRACKF_STAT(IDRES,IREADGEANT)
 461 *************************************************************************
 462 *  Associate hits between two chambers inside a station
 463 *  Simulate spatial resolution and chamber efficiency
 464 *
 465 *************************************************************************
 466       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
 467 *
 468       PARAMETER (NBSTATION=5,MAXIDG=20000,MAXHITTOT=20000,
 469      &           MAXHITCH=10000,MAXHIT=1000,NBCHAMBER=10)
 470 *
 471       COMMON/TRACKFI/EFF,EFF1,EFF2,XPREC,YPREC,PHIPREC,ALAMPREC,
 472      &                 HCUT,LBKG,SIGCUT,ALPHATOP,HTOP
 473 *
 474       COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
 475 *
 476 * HITS GEANT initiaux par chambre
 477       COMMON/RHITG/ITYPG(MAXIDG),XTRG(MAXIDG),YTRG(MAXIDG),
 478      &     PTOTG(MAXIDG),IDG(MAXIDG),IZCH(MAXIDG),
 479      &     PVERT1G(MAXIDG),PVERT2G(MAXIDG),PVERT3G(MAXIDG),
 480      &     ZVERTG(MAXIDG),NHITTOT1,CX(MAXIDG),CY(MAXIDG),CZ(MAXIDG),
 481      &     XGEANT(MAXIDG),YGEANT(MAXIDG),CLSIZE1(MAXIDG),CLSIZE2(MAXIDG)
 482
 483 * HITS GEANT associes par station
 484       COMMON/RHIT/ITYP(MAXHITTOT),XTR(MAXHITTOT),YTR(MAXHITTOT),
 485      &            PTOT(MAXHITTOT),ID(MAXHITTOT),IZST(MAXHITTOT),
 486      &            PVERT1(MAXHITTOT),PVERT2(MAXHITTOT),PVERT3(MAXHITTOT),
 487      &            ZVERT(MAXHITTOT),NHITTOT
 488
 489 *
 490       COMMON/CHHIT/XM(NBSTATION,MAXHITCH),YM(NBSTATION,MAXHITCH),
 491      &             PHM(NBSTATION,MAXHITCH),ALM(NBSTATION,MAXHITCH),
 492      &             IZM(NBSTATION,MAXHITCH),
 493      &             IP(NBSTATION,MAXHITCH),JHIT(NBSTATION),
 494      &             XMR(NBSTATION,MAXHITCH,2),YMR(NBSTATION,MAXHITCH,2)
 495 *
 496       COMMON/DEBEVT/IDEBUG
 497 *
 498       DIMENSION RMIN(NBCHAMBER),RMAX1(NBCHAMBER)
 499       DIMENSION XMA(NBCHAMBER,MAXHITCH),YMA(NBCHAMBER,MAXHITCH),
 500      &          IMARK(NBCHAMBER,MAXHITCH)
 501 *
 502       DIMENSION IEFFI(MAXHITTOT)
 503       DIMENSION IH(NBCHAMBER,MAXHIT)
 504       DIMENSION NHIT(NBCHAMBER)
 505       DIMENSION DXMAX(NBSTATION),DYMAX(NBSTATION),VRES(2,5)
 506
 507       REAL*4 RNDM,RN,RN1,RN2,R1,R2
 508
 509 * Chambre 10 deg.
 510       DATA RMAX1/91.5,91.5,122.5,122.5,158.3,158.3,260.,260.,260.,260./
 511 * Zone de recherche entre deux plans d'une station
 512 **      DATA DXMAX/2.,1.5,2.5,3.,3./
 513       DATA DXMAX/1.5,1.5,3.,6.,6./  ! J psi 20cm
 514 **      DATA DXMAX/1.,1.,1.,1.,1./  ! Upsilon dz_ch = 8 cm
 515 **      DATA DYMAX/0.22,0.22,0.21,0.21,0.21/  ! CCC Upsilon dz_ch = 8cm
 516 **      DATA DXMAX/1.,1.,1.5,2.,2./  ! Upsilon dz_ch = 20 cm
 517       DATA DYMAX/0.22,0.22,0.22,0.22,0.22/  ! CCC Upsilon dz_ch = 20 cm
 518 **      DATA DXMAX/5.,5.,5.,5.,5./
 519 **      DATA DXMAX/10.,10.,10.,10.,10./
 520       DATA R1,R2/0.,1./
 521
 522       ICH = 0
 523       DO IZ=1,5
 524          ICH = ICH+1
 525          RMIN(ICH) =  ABS(ZPLANE(IZ)*TAN(2.*ACOS(-1.)/180))
 526          IF (IZ.GT.2) RMIN(ICH) = 30.
 527          ICH = ICH+1
 528          RMIN(ICH) =  ABS(ZPLANE(IZ)*TAN(2.*ACOS(-1.)/180))
 529          IF (IZ.GT.2) RMIN(ICH) = 30.
 530       ENDDO
 531
 532 *    Initialisations
 533       DO ICH = 1,10
 534          NHIT(ICH) = 0
 535       ENDDO
 536
 537 * 1 ere boucle Lecture des hits initiaux
 538
 539       print*,'TRACKF_STAT NHITTOT1=',NHITTOT1
 540
 541       IF (IREADGEANT.EQ.1) THEN
 542          DO I = 1,2
 543             DO J = 1,5
 544                VRES(I,J) = 0.
 545             ENDDO
 546          ENDDO
 547          IMU = 0
 548          DO I = 1,NHITTOT1
 549             ICH = IZCH(I)
 550             IF (ICH.EQ.9) THEN
 551                ISTAK = IDG(I)
 552                ISTAK = MOD(ISTAK,30000)
 553                ISTAK = MOD(ISTAK,10000)
 554
 555                IF ((ITYPG(I).EQ.5.OR.ITYPG(I).EQ.6).AND.
 556      &              ISTAK.EQ.IDRES) THEN ! upsilon
 557                   IMU = IMU+1
 558                   VRES(1,IMU) = XGEANT(I)
 559                   VRES(2,IMU) = YGEANT(I)
 560                ENDIF
 561             ENDIF
 562          ENDDO
 563       ENDIF
 564
 565
 566       DO I = 1,NHITTOT1 ! Boucle sur les hits GEANT de toutes les ch.
 567
 568 **        IF (ITYPG(I).NE.5.AND.ITYPG(I).NE.6) GOTO 1 ! CCC
 569
 570          ICH = IZCH(I)
 571
 572          IF (IREADGEANT.EQ.1) THEN ! GEANT hits
 573
 574             IF (ICH.EQ.9.OR.ICH.EQ.10) THEN
 575                DNUM = 999.
 576                DO IM = 1,IMU
 577                   DNU = SQRT((XGEANT(I)-VRES(1,IM))**2+
 578      &                 (YGEANT(I)-VRES(2,IM))**2)
 579                   IF (DNU.LT.DNUM) DNUM = DNU
 580                ENDDO
 581                IF (DNUM.GT.20.) GO TO 1
 582             ENDIF
 583
 584             CALL RANNOR(RN1,RN2)          ! CCCC
 585             X = XGEANT(I) + RN1 * XPREC
 586             Y = YGEANT(I) + RN2 * YPREC
 587 *     efficacite des chambres
 588             IEFFI(I) = 1
 589             RN = RNDM()
 590             IF (RN.GT.EFF) IEFFI(I) = 0
 591
 592
 593          ELSE  ! reconstructed hits
 594
 595             IEFFI(I) = 1
 596
 597             X = XTRG(I)
 598             Y = YTRG(I)
 599
 600          ENDIF
 601
 602          R = SQRT(X**2+Y**2)
 603          IF (R.LT.RMIN(ICH).OR.R.GT.RMAX1(ICH)) then
 604             if (ich.le.10.and.i.le.28) then
 605                print*,'* chambre ',ich,' * hit ',i
 606                print*,'ityp=',itypg(i)
 607                print*,'x=',X,' y=',Y
 608                print*,'R=',R,' RMIN=',RMIN(ICH),' RMAX1=',RMAX1(ICH)
 609            endif
 610             GO TO 1
 611          end if
 612
 613          NHIT(ICH) = NHIT(ICH)+1
 614          IH(ICH,NHIT(ICH)) = I
 615          XMA(ICH,NHIT(ICH)) = X
 616          YMA(ICH,NHIT(ICH)) = Y
 617          IMARK(ICH,NHIT(ICH)) = 0
 618
 619   1      CONTINUE
 620       ENDDO
 621
 622 * Association des hits entre chambres d'une station
 623       II = 0            ! nombre de hits GEANT par station
 624       DO ICH1 = 1,10,2
 625          IZ = INT(FLOAT(ICH1+1)/2.)
 626          JHIT(IZ) = 0
 627          ICH2 = ICH1+1
 628
 629          DO I1 = 1,NHIT(ICH1)
 630             II = II+1
 631             IFIND = 0
 632             I = IH(ICH1,I1)
 633
 634             ITYP(II) = ITYPG(I)
 635             XTR(II) = XTRG(I)
 636             YTR(II) = YTRG(I)
 637             PTOT(II) = PTOTG(I)
 638             ID(II) = IDG(I)
 639             IZST(II) = IZ
 640             PVERT1(II) = PVERT1G(I)
 641             PVERT2(II) = PVERT2G(I)
 642             PVERT3(II) = PVERT3G(I)
 643             ZVERT(II) = ZVERTG(I)
 644
 645             IF (IEFFI(I).EQ.1) THEN
 646                X1 = XMA(ICH1,I1)
 647                Y1 = YMA(ICH1,I1)
 648                ID1 = IDG(I)
 649                XEXT1 = (ZPLANE(IZ)-DZ_PL(IZ))/ZPLANE(IZ)*X1
 650                YEXT1 = (ZPLANE(IZ)-DZ_PL(IZ))/ZPLANE(IZ)*Y1
 651
 652                DO I2 = 1,NHIT(ICH2)
 653                   J = IH(ICH2,I2)
 654
 655                   IF (IEFFI(J).EQ.1) THEN
 656                      X2 = XMA(ICH2,I2)
 657                      Y2 = YMA(ICH2,I2)
 658                      ID2 = IDG(J)
 659                      DX = X2-XEXT1
 660                      DY = Y2-YEXT1
 661
 662                      IF (ID1.EQ.ID2.AND.
 663      &                    (ITYP(II).EQ.5.OR.ITYP(II).EQ.6)) THEN
 664                         CALL CHFILL(70+IZ,SNGL(DX),R1,R2)
 665                         CALL CHFILL(80+IZ,SNGL(DY),R1,R2)
 666                      ENDIF
 667                      DX = ABS(DX)
 668                      DY = ABS(DY)
 669
 670                      IF (DX.LT.DXMAX(IZ).AND.DY.LT.(SIGCUT*DYMAX(IZ)) ! CCC
 671      &                   ) THEN
 672                         IFIND = 1
 673                         IMARK(ICH2,I2) = 1
 674                         JHIT(IZ) =  JHIT(IZ)+1
 675                         XM(IZ,JHIT(IZ)) = X1
 676                         YM(IZ,JHIT(IZ)) = Y1
 677                         IZM(IZ,JHIT(IZ)) = 1
 678                         PHM(IZ,JHIT(IZ)) = -ATAN((X2-X1)/DZ_PL(IZ))
 679                         ALM(IZ,JHIT(IZ)) = ATAN((Y2-Y1)/DZ_PL(IZ)*
 680      &                                     COS(PHM(IZ,JHIT(IZ))))
 681                         XMR(IZ,JHIT(IZ),1) = X1
 682                         YMR(IZ,JHIT(IZ),1) = Y1
 683                         XMR(IZ,JHIT(IZ),2) = X2
 684                         YMR(IZ,JHIT(IZ),2) = Y2
 685
 686                         ISTAK = ID2
 687                         ISTAK = MOD(ISTAK,30000)
 688                         ISTAK = MOD(ISTAK,10000)
 689
 690                         IF ((ITYPG(J).EQ.5.OR.ITYPG(J).EQ.6).AND.
 691      &                       ISTAK.EQ.IDRES) THEN ! upsilon
 692
 693                            ITYP(II) = ITYPG(J)
 694                            XTR(II) = XTRG(J)
 695                            YTR(II) = YTRG(J)
 696                            PTOT(II) = PTOTG(J)
 697                            ID(II) = IDG(J)
 698                            PVERT1(II) = PVERT1G(J)
 699                            PVERT2(II) = PVERT2G(J)
 700                            PVERT3(II) = PVERT3G(J)
 701                            ZVERT(II) = ZVERTG(J)
 702
 703                         ENDIF
 704
 705                         IP(IZ,JHIT(IZ)) = II
 706
 707                      ENDIF
 708                   ENDIF
 709                ENDDO
 710
 711                IF (IFIND.EQ.0) THEN
 712                   JHIT(IZ) =  JHIT(IZ)+1
 713                   XM(IZ,JHIT(IZ)) = X1
 714                   YM(IZ,JHIT(IZ)) = Y1
 715                   IZM(IZ,JHIT(IZ)) = 1
 716                   IP(IZ,JHIT(IZ)) = II
 717                   PHM(IZ,JHIT(IZ)) = 10.
 718                   ALM(IZ,JHIT(IZ)) = 10.
 719                   XMR(IZ,JHIT(IZ),1) = X1
 720                   YMR(IZ,JHIT(IZ),1) = Y1
 721                   XMR(IZ,JHIT(IZ),2) = 0.
 722                   YMR(IZ,JHIT(IZ),2) = 0.
 723                ENDIF
 724                CALL CHFILL2(1000+IZ,SNGL(X1),SNGL(Y1),R2)
 725             ENDIF
 726          ENDDO
 727       ENDDO
 728
 729 * On conserve les HITS de la 2nde chambre des stations
 730
 731       DO ICH = 2,10,2
 732          IZ = INT(FLOAT(ICH+1)/2.)
 733          DO I = 1,NHIT(ICH)
 734             J = IH(ICH,I)
 735
 736             IF (IMARK(ICH,I).EQ.0) THEN
 737
 738                II = II+1
 739
 740                ITYP(II) = ITYPG(J)
 741                XTR(II) = XTRG(J)
 742                YTR(II) = YTRG(J)
 743                PTOT(II) = PTOTG(J)
 744                ID(II) = IDG(J)
 745                IZST(II) = IZ
 746                PVERT1(II) = PVERT1G(J)
 747                PVERT2(II) = PVERT2G(J)
 748                PVERT3(II) = PVERT3G(J)
 749                ZVERT(II) = ZVERTG(I)
 750
 751                IF (IEFFI(J).EQ.1) THEN
 752                   JHIT(IZ) =  JHIT(IZ)+1
 753                   XM(IZ,JHIT(IZ)) = XMA(ICH,I)
 754                   YM(IZ,JHIT(IZ)) = YMA(ICH,I)
 755                   IZM(IZ,JHIT(IZ)) = 2
 756                   PHM(IZ,JHIT(IZ)) = 10.
 757                   ALM(IZ,JHIT(IZ)) = 10.
 758                   IP(IZ,JHIT(IZ)) = II
 759                   XMR(IZ,JHIT(IZ),1) = 1000.
 760                   YMR(IZ,JHIT(IZ),1) = 1000.
 761                   XMR(IZ,JHIT(IZ),2) = XMA(ICH,I)
 762                   YMR(IZ,JHIT(IZ),2) = YMA(ICH,I)
 763                ENDIF
 764             ENDIF
 765          ENDDO
 766       ENDDO
 767
 768
 769       NHITTOT = II
 770 *
 771
 772       IF (IDEBUG.GE.2) THEN
 773          PRINT *,'TRACKF_STAT nb hits:',NHITTOT
 774       ENDIF
 775
 776 **      DO IZ = 1,NBSTATION
 777 **         PRINT *,' IZ=',IZ,' JHIT(IZ)=',JHIT(IZ)
 778 **         DO J = 1,JHIT(IZ)
 779 **            II = IP(IZ,J)
 780 **            print*,'II=',II
 781 **           PRINT *,' ID(II)=',ID(II)
 782 **         ENDDO
 783 **      ENDDO
 784
 785 **    DO IZ = 1,NBSTATION
 786 **       PRINT *,' IZ=',IZ,' JHIT(IZ)=',JHIT(IZ)
 787 **       DO J = 1,JHIT(IZ)
 788 **          PRINT *,' PHM(IZ,J)=',PHM(IZ,J),' ALM(IZ,J)=',ALM(IZ,J)
 789 **       ENDDO
 790 **     ENDDO
 791
 792       RETURN
 793       END
 794 *************************************************************************
 795       SUBROUTINE TRACKF_STAT_NEW(IDRES)
 796 *************************************************************************
 797 *  Associate hits between two chambers inside a station
 798 *  Simulate spatial resolution and chamber efficiency
 799 *
 800 *************************************************************************
 801       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
 802 *
 803       PARAMETER (NBSTATION=5,MAXIDG=20000,MAXHITTOT=20000,
 804      &           MAXHITCH=10000,MAXHIT=1000,NBCHAMBER=10)
 805 *
 806       COMMON/TRACKFI/EFF,EFF1,EFF2,XPREC,YPREC,PHIPREC,ALAMPREC,
 807      &                 HCUT,LBKG,SIGCUT,ALPHATOP,HTOP
 808 *
 809       COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
 810 *
 811 * HITS GEANT initiaux par chambre
 812       COMMON/RHITG/ITYPG(MAXIDG),XTRG(MAXIDG),YTRG(MAXIDG),
 813      &     PTOTG(MAXIDG),IDG(MAXIDG),IZCH(MAXIDG),
 814      &     PVERT1G(MAXIDG),PVERT2G(MAXIDG),PVERT3G(MAXIDG),
 815      &     ZVERTG(MAXIDG),NHITTOT1,CX(MAXIDG),CY(MAXIDG),CZ(MAXIDG),
 816      &     XGEANT(MAXIDG),YGEANT(MAXIDG),CLSIZE1(MAXIDG),CLSIZE2(MAXIDG)
 817
 818 * HITS GEANT associes par station
 819       COMMON/RHIT/ITYP(MAXHITTOT),XTR(MAXHITTOT),YTR(MAXHITTOT),
 820      &            PTOT(MAXHITTOT),ID(MAXHITTOT),IZST(MAXHITTOT),
 821      &            PVERT1(MAXHITTOT),PVERT2(MAXHITTOT),PVERT3(MAXHITTOT),
 822      &            ZVERT(MAXHITTOT),NHITTOT
 823
 824 *
 825       COMMON/CHHIT/XM(NBSTATION,MAXHITCH),YM(NBSTATION,MAXHITCH),
 826      &             PHM(NBSTATION,MAXHITCH),ALM(NBSTATION,MAXHITCH),
 827      &             IZM(NBSTATION,MAXHITCH),
 828      &             IP(NBSTATION,MAXHITCH),JHIT(NBSTATION),
 829      &             XMR(NBSTATION,MAXHITCH,2),YMR(NBSTATION,MAXHITCH,2)
 830 *
 831       COMMON/DEBEVT/IDEBUG
 832 *
 833       DIMENSION RMIN(NBCHAMBER),RMAX1(NBCHAMBER)
 834       DIMENSION XMA(NBCHAMBER,MAXHITCH),YMA(NBCHAMBER,MAXHITCH),
 835      &          IMARK(NBCHAMBER,MAXHITCH)
 836 *
 837       DIMENSION IEFFI(MAXHITTOT)
 838       DIMENSION IH(NBCHAMBER,MAXHIT)
 839       DIMENSION NHIT(NBCHAMBER)
 840       DIMENSION DXMAX(NBSTATION),DYMAX(NBSTATION),I2C(1000)
 841
 842       DIMENSION DIST(2),NMUON(2),NHITMUON(2,5),NMUONGOOD(2)
 843
 844 *      REAL*4 RNDM,RN,RN1,RN2,R1,R2
 845
 846 * Chambre 10 deg.
 847       DATA RMAX1/91.5,91.5,122.5,122.5,158.3,158.3,260.,260.,260.,260./
 848 * Zone de recherche entre deux plans d'une station
 849 **      DATA DXMAX/2.,1.5,2.5,3.,3./
 850       DATA DXMAX/1.5,1.5,3.,3.,3./
 851       DATA DYMAX/0.22,0.22,0.21,0.21,0.21/
 852       DATA R1,R2/0.,1./
 853
 854       ICH = 0
 855       DO IZ=1,5
 856          ICH = ICH+1
 857          RMIN(ICH) =  ABS(ZPLANE(IZ)*TAN(2.*ACOS(-1.)/180))
 858          IF (IZ.GT.2) RMIN(ICH) = 30.
 859          ICH = ICH+1
 860          RMIN(ICH) =  ABS(ZPLANE(IZ)*TAN(2.*ACOS(-1.)/180))
 861          IF (IZ.GT.2) RMIN(ICH) = 30.
 862       ENDDO
 863
 864 *    Initialisations
 865       DO ICH = 1,10
 866          NHIT(ICH) = 0
 867       ENDDO
 868
 869       DO NCH = 1,10
 870          DO J=1,2
 871             DIST(J)=999.
 872             NMUON(J)=0
 873          ENDDO
 874          DO I = 1,NHITTOT1
 875             IF (IZCH(I).EQ.NCH) THEN
 876                ISTAK = IDG(I)
 877                ISTAK = MOD(ISTAK,30000)
 878                ISTAK = MOD(ISTAK,10000)
 879                IF (ISTAK.EQ.IDRES.AND.IDG(I).EQ.50116) THEN
 880                   DISTMIN=(XTRG(I)-XGEANT(I))**2+(YTRG(I)-YGEANT(I))**2
 881                   IF (DISTMIN.LT.DIST(1)) THEN
 882                      DIST(1)=DISTMIN
 883                      NMUONGOOD(1)=I
 884                   ENDIF
 885                   NMUON(1)=NMUON(1)+1
 886                   NHITMUON(1,NMUON(1))=I
 887                ELSE
 888                 IF (ISTAK.EQ.IDRES.AND.IDG(I).EQ.70116) THEN
 889                    DISTMIN=(XTRG(I)-XGEANT(I))**2+(YTRG(I)-YGEANT(I))**2
 890                    IF (DISTMIN.LT.DIST(2)) THEN
 891                         DIST(2)=DISTMIN
 892                         NMUONGOOD(2)=I
 893                    ENDIF
 894                    NMUON(2)=NMUON(2)+1
 895                    NHITMUON(2,NMUON(2))=I
 896                   ENDIF
 897                ENDIF
 898             ENDIF
 899          ENDDO
 900          DO J=1,2
 901             IF (NMUON(J).GE.2) THEN
 902              print*,'j=',j,' nmuon=',nmuon(j)
 903              print*,'chambre',nch
 904              DO K=1,NMUON(J)
 905                IF (NHITMUON(J,K).NE.NMUONGOOD(J)) IDG(NHITMUON(J,K))=999 ! flag les mauvais hits MUONS
 906              ENDDO
 907             ENDIF
 908          ENDDO
 909       ENDDO
 910
 911
 912 * 1 ere boucle Lecture des hits initiaux
 913
 914       DO I = 1,NHITTOT1 ! Boucle sur les hits GEANT de toutes les ch.
 915
 916          ICH = IZCH(I)
 917
 918          X = XTRG(I)
 919          Y = YTRG(I)
 920
 921          R = SQRT(X**2+Y**2)
 922          IF (R.LT.RMIN(ICH).OR.R.GT.RMAX1(ICH)) then
 923             if (ich.le.10.and.i.le.28) then
 924                print*,'****** chambre ',ich,' ****** hit ',i
 925                print*,'ityp=',itypg(i)
 926                print*,'x=',XTRG(I),' y=',YTRG(I)
 927                print*,'R=',R,' RMIN=',RMIN(ICH),' RMAX1=',RMAX1(ICH)
 928            endif
 929             GO TO 1
 930          end if
 931
 932          IEFFI(I) = 1
 933
 934          NHIT(ICH) = NHIT(ICH)+1
 935          IH(ICH,NHIT(ICH)) = I
 936          XMA(ICH,NHIT(ICH)) = XTRG(I)
 937          YMA(ICH,NHIT(ICH)) = YTRG(I)
 938          IMARK(ICH,NHIT(ICH)) = 0
 939
 940          print*,' XTRG(I)=', XTRG(I),' YTRG(I)=', YTRG(I),' IDG(I)=',
 941      &        IDG(I),' ICH=',ICH
 942
 943   1      CONTINUE
 944       ENDDO
 945
 946
 947 * Association des hits entre chambres d'une station
 948       II = 0            ! nombre de hits GEANT par station
 949       DO ICH1 = 1,10,2
 950          IZ = INT(FLOAT(ICH1+1)/2.)
 951          JHIT(IZ) = 0
 952          ICH2 = ICH1+1
 953
 954          DO I1 = 1,NHIT(ICH1)
 955             II = II+1
 956             IFIND = 0
 957             I = IH(ICH1,I1)
 958
 959             ITYP(II) = ITYPG(I)
 960             XTR(II) = XTRG(I)
 961             YTR(II) = YTRG(I)
 962             PTOT(II) = PTOTG(I)
 963             ID(II) = IDG(I)
 964             IZST(II) = IZ
 965             PVERT1(II) = PVERT1G(I)
 966             PVERT2(II) = PVERT2G(I)
 967             PVERT3(II) = PVERT3G(I)
 968             ZVERT(II) = ZVERTG(I)
 969
 970             IF (IEFFI(I).EQ.1) THEN
 971                X1 = XMA(ICH1,I1)
 972                Y1 = YMA(ICH1,I1)
 973                ID1 = IDG(I)
 974                XEXT1 = (ZPLANE(IZ)-DZ_PL(IZ))/ZPLANE(IZ)*X1
 975                YEXT1 = (ZPLANE(IZ)-DZ_PL(IZ))/ZPLANE(IZ)*Y1
 976                KC = 0
 977                PRINT *,'***** DEBUT RECHERCHE',' ID1=',ID1,' ich1=',ICH1
 978                PRINT *,' XTR(II)=', XTR(II),' YTR(II)=', YTR(II)
 979                PRINT *,'  ITYP(II)=', ITYP(II)
 980                DO I2 = 1,NHIT(ICH2)
 981                   J = IH(ICH2,I2)
 982                   IF (IEFFI(J).EQ.1) THEN
 983                      X2 = XMA(ICH2,I2)
 984                      DX = X2-XEXT1
 985                      DX = ABS(DX)
 986                      IF (DX.LT.DXMAX(IZ)) THEN
 987                         KC = KC + 1
 988                         I2C(KC) = I2
 989                         ID2 = IDG(J)
 990                         print *,' DX=',DX,' KC=',KC,' ID2=',ID2
 991                      ENDIF
 992                   ENDIF
 993                ENDDO
 994                DYOLD = 999.
 995                I2FIND = 0
 996                DO IKC = 1,KC
 997                   I2 = I2C(IKC)
 998                   Y2 = YMA(ICH2,I2)
 999                   DY = Y2-YEXT1
1000                   DY = ABS(DY)
1001                   J = IH(ICH2,I2)
1002                   ID2 = IDG(J)
1003                   IF (DY.LT.DYOLD.AND.DY.LT.(SIGCUT*DYMAX(IZ))) THEN
1004                      DYOLD = DY
1005                      I2FIND = I2
1006                      PRINT *,' ID2=',ID2,' DY=',DY
1007                   ENDIF
1008                ENDDO
1009                IF (I2FIND.GT.0) THEN
1010                   I2 = I2FIND
1011                   J = IH(ICH2,I2)
1012                   ID2 = IDG(J)
1013                   IFIND = 1
1014                   IMARK(ICH2,I2) = 1
1015                   JHIT(IZ) =  JHIT(IZ)+1
1016                   X2 = XMA(ICH2,I2)
1017                   Y2 = YMA(ICH2,I2)
1018                   XM(IZ,JHIT(IZ)) = X1
1019                   YM(IZ,JHIT(IZ)) = Y1
1020                   IZM(IZ,JHIT(IZ)) = 1
1021                   PHM(IZ,JHIT(IZ)) = -ATAN((X2-X1)/DZ_PL(IZ))
1022                   ALM(IZ,JHIT(IZ)) = ATAN((Y2-Y1)/DZ_PL(IZ)*
1023      &                 COS(PHM(IZ,JHIT(IZ))))
1024                   XMR(IZ,JHIT(IZ),1) = X1
1025                   YMR(IZ,JHIT(IZ),1) = Y1
1026                   XMR(IZ,JHIT(IZ),2) = X2
1027                   YMR(IZ,JHIT(IZ),2) = Y2
1028                   IP(IZ,JHIT(IZ)) = II
1029
1030                   ISTAK = ID2
1031                   ISTAK = MOD(ISTAK,30000)
1032                   ISTAK = MOD(ISTAK,10000)
1033 * test
1034                   IF (ISTAK.EQ.IDRES.AND.ID1.NE.999) THEN
1035
1036                      ITYP(II) = ITYPG(J)
1037                      PTOT(II) = PTOTG(J)
1038                      XTR(II) = XTRG(I)
1039                      YTR(II) = YTRG(I)
1040                      ID(II) = IDG(J)
1041                      PVERT1(II) = PVERT1G(J)
1042                      PVERT2(II) = PVERT2G(J)
1043                      PVERT3(II) = PVERT3G(J)
1044                      ZVERT(II) = ZVERTG(J)
1045
1046                   ENDIF
1047                ENDIF
1048
1049
1050
1051                IF (IFIND.EQ.0) THEN
1052                   JHIT(IZ) =  JHIT(IZ)+1
1053                   XM(IZ,JHIT(IZ)) = X1
1054                   YM(IZ,JHIT(IZ)) = Y1
1055                   IZM(IZ,JHIT(IZ)) = 1
1056                   IP(IZ,JHIT(IZ)) = II
1057                   PHM(IZ,JHIT(IZ)) = 10.
1058                   ALM(IZ,JHIT(IZ)) = 10.
1059                   XMR(IZ,JHIT(IZ),1) = X1
1060                   YMR(IZ,JHIT(IZ),1) = Y1
1061                   XMR(IZ,JHIT(IZ),2) = 0.
1062                   YMR(IZ,JHIT(IZ),2) = 0.
1063                ENDIF
1064                CALL CHFILL2(1000+IZ,SNGL(X1),SNGL(Y1),R2)
1065             ENDIF
1066          ENDDO
1067       ENDDO
1068
1069 * On conserve les HITS de la 2nde chambre des stations
1070
1071       DO ICH = 2,10,2
1072          IZ = INT(FLOAT(ICH+1)/2.)
1073          DO I = 1,NHIT(ICH)
1074             J = IH(ICH,I)
1075
1076             IF (IMARK(ICH,I).EQ.0) THEN
1077 **               print *,' ich=',ich,' i=',i,' j=',j
1078
1079                II = II+1
1080
1081                ITYP(II) = ITYPG(J)
1082                XTR(II) = XTRG(J)
1083                YTR(II) = YTRG(J)
1084                PTOT(II) = PTOTG(J)
1085                ID(II) = IDG(J)
1086                IZST(II) = IZ
1087                PVERT1(II) = PVERT1G(J)
1088                PVERT2(II) = PVERT2G(J)
1089                PVERT3(II) = PVERT3G(J)
1090                ZVERT(II) = ZVERTG(I)
1091
1092                IF (IEFFI(J).EQ.1) THEN
1093                   JHIT(IZ) =  JHIT(IZ)+1
1094 **                  XM(IZ,JHIT(IZ)) = ZPLANE(IZ)/(ZPLANE(IZ)-DZ_PL)
1095 **     &                              *XMA(ICH,I)
1096 **                  YM(IZ,JHIT(IZ)) = ZPLANE(IZ)/(ZPLANE(IZ)-DZ_PL)
1097 **     &                              *YMA(ICH,I)
1098                   XM(IZ,JHIT(IZ)) = XMA(ICH,I)
1099                   YM(IZ,JHIT(IZ)) = YMA(ICH,I)
1100                   IZM(IZ,JHIT(IZ)) = 2
1101                   PHM(IZ,JHIT(IZ)) = 10.
1102                   ALM(IZ,JHIT(IZ)) = 10.
1103                   IP(IZ,JHIT(IZ)) = II
1104                   XMR(IZ,JHIT(IZ),1) = 1000.
1105                   YMR(IZ,JHIT(IZ),1) = 1000.
1106                   XMR(IZ,JHIT(IZ),2) = XMA(ICH,I)
1107                   YMR(IZ,JHIT(IZ),2) = YMA(ICH,I)
1108                ENDIF
1109             ENDIF
1110          ENDDO
1111       ENDDO
1112
1113
1114       NHITTOT = II
1115 *
1116       IF (IDEBUG.GE.2) THEN
1117          PRINT *,'TRACKF_MICRO nb hits:',NHITTOT
1118       ENDIF
1119       DO IZ = 1,NBSTATION
1120          PRINT *,' IZ=',IZ,' JHIT(IZ)=',JHIT(IZ)
1121          DO J = 1,JHIT(IZ)
1122             II = IP(IZ,J)
1123             PRINT *,' ID(II)=',ID(II)
1124             PRINT *,' XMR(IZ,J,1)=', XMR(IZ,J,1),
1125      &            ' YMR(IZ,J,1)=', YMR(IZ,J,1)
1126             PRINT *,' XMR(IZ,J,2)=', XMR(IZ,J,2),
1127      &            ' YMR(IZ,J,2)=', YMR(IZ,J,2)
1128         ENDDO
1129       ENDDO
1130
1131       RETURN
1132       END
1133
1134 *************************************************************************
1135       SUBROUTINE RECO_TRACKF(IDRES,IREADGEANT)
1136 *************************************************************************
1137 *
1138 *
1139 *************************************************************************
1140       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
1141 *
1142       PARAMETER (MAXIDG=20000)
1143 *
1144       PARAMETER (MAXHITCH=10000,MAXTRK=50000,MAXHITTOT=20000,
1145      &           NBSTATION=5,MAXCAN=1000,NBCHAMBER=10,NTRMAX=500)
1146 *
1147       COMMON/CHHIT/XM(NBSTATION,MAXHITCH),YM(NBSTATION,MAXHITCH),
1148      &             PHM(NBSTATION,MAXHITCH),ALM(NBSTATION,MAXHITCH),
1149      &             IZM(NBSTATION,MAXHITCH),
1150      &             IP(NBSTATION,MAXHITCH),JHIT(NBSTATION),
1151      &             XMR(NBSTATION,MAXHITCH,2),YMR(NBSTATION,MAXHITCH,2)
1152 *
1153       COMMON/RHIT/ITYP(MAXHITTOT),XTR(MAXHITTOT),YTR(MAXHITTOT),
1154      &            PTOT(MAXHITTOT),ID(MAXHITTOT),IZST(MAXHITTOT),
1155      &            PVERT1(MAXHITTOT),PVERT2(MAXHITTOT),PVERT3(MAXHITTOT),
1156      &            ZVERT(MAXHITTOT),NHITTOT
1157
1158 *
1159       COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
1160 *
1161       COMMON/VERIFGEANT/ITTROUGH(MAXTRK,NBSTATION),
1162      &                  IT_LIST(MAXTRK),IT_NP(MAXTRK),ITCHECK(MAXTRK),
1163      &                  ITRACK(MAXHITTOT)
1164 *
1165       COMMON/HCHHIT/HHIT(MAXHITCH),INDEXTAB(MAXHITCH),INDEXMAX
1166 *
1167       COMMON/PRECIS/EEXM(NBSTATION),EEYM(NBSTATION),EEPH(NBSTATION),
1168      &     EEAL(NBSTATION)
1169 *
1170       COMMON/MEASUR/XMES(NBSTATION),YMES(NBSTATION),IZMES(NBSTATION),
1171      &              PHMES(NBSTATION),ALMES(NBSTATION),MPOS(NBSTATION),
1172      &              MANG(NBSTATION)
1173 *
1174       COMMON/PLANE/XPL(NBSTATION,2),YPL(NBSTATION,2),
1175      &             PHPL(NBSTATION),ALPL(NBSTATION),CHI2PL
1176 *
1177       COMMON/CANDIDAT/JCAN(NBSTATION,MAXCAN),JCANTYP(MAXCAN),
1178      &          EEX(MAXCAN),EEY(MAXCAN),EEP(MAXCAN),EEA(MAXCAN)
1179 *
1180       COMMON /VERTEX/ERRV,IVERTEX
1181 *
1182       COMMON/TRACKSUM/NRES(5),NRESF,NTRMUALL,NMUONALL,NGHOSTALL,
1183      &     NTRACKFALL,NERRALL(NBSTATION),IR
1184 *
1185       COMMON/TRACKFOUT/IEVOUT,NTREVT,JJOUT(NBCHAMBER,NTRMAX),
1186      &               ISTAT(NTRMAX),PXZOUT(NTRMAX),TPHIOUT(NTRMAX),
1187      &               TALAMOUT(NTRMAX),XVERTOUT(NTRMAX),YVERTOUT(NTRMAX),
1188      &               CHI2OUT(NTRMAX),
1189      &               XMESOUT(NBCHAMBER,NTRMAX),YMESOUT(NBCHAMBER,NTRMAX)
1190      &              ,PXVOUT(NTRMAX),PYVOUT(NTRMAX),PZVOUT(NTRMAX)
1191 *
1192       COMMON/TRACKFI/EFF,EFF1,EFF2,XPREC,YPREC,PHIPREC,ALAMPREC,
1193      &                 HCUT,LBKG,SIGCUT,ALPHATOP,HTOP
1194 *
1195       COMMON/DEBEVT/IDEBUG
1196
1197       DIMENSION PEST(5),PSTEP(5),NERR(NBSTATION)
1198 *
1199       REAL*4 R2
1200       DATA R2/1./
1201
1202 ** GEANT informations
1203       DO I = 1,MAXTRK
1204          IT_LIST(I)= 0
1205       ENDDO
1206       DO I = 1,MAXHITTOT
1207          ITRACK(I) = 0
1208       ENDDO
1209       NTRTOT = 0
1210 *  BOUCLE SUR LES HITS
1211
1212       DO IH = 1,NHITTOT
1213          DO IT = 1,NTRTOT
1214             IF (ID(IH).EQ.IT_LIST(IT))  GOTO 4
1215          ENDDO
1216          NTRTOT = NTRTOT +1     ! NB de trace GEANT
1217          IT_LIST(NTRTOT) = ID(IH) ! ID DE LA TRACE NTRTOT
1218          IT_NP(NTRTOT) = 0      ! NOMBRE DE HITS PAR TRACE
1219          DO II=1,NBSTATION
1220             ITTROUGH(NTRTOT,II) = 0
1221          ENDDO
1222          IT = NTRTOT
1223  4       IT_NP(IT) = IT_NP(IT)+1
1224          ITTROUGH(IT,IZST(IH))=IH ! =IH si la trace IT touche le plan IZST
1225          ITRACK(IH) = IT
1226       ENDDO
1227       IF (IDEBUG.GE.2) THEN
1228          PRINT *,'RECO_TRACKF nb total de trace GEANT:',NTRTOT
1229       ENDIF
1230       NTRPART=0
1231       NTRMU = 0
1232
1233       DO IT = 1,NTRTOT
1234 *     print*,'(ITTROUGH(IT,1)=',ITTROUGH(IT,1)
1235 *     print*,'(ITTROUGH(IT,2)=',ITTROUGH(IT,2)
1236 *     print*,'(ITTROUGH(IT,3)=',ITTROUGH(IT,3)
1237 *     print*,'(ITTROUGH(IT,4)=',ITTROUGH(IT,4)
1238 *     print*,'(ITTROUGH(IT,5)=',ITTROUGH(IT,5)
1239          ITCHECK(IT) = 0
1240          IF ((ITTROUGH(IT,1)*ITTROUGH(IT,2)*ITTROUGH(IT,3)*
1241      &        ITTROUGH(IT,4)*ITTROUGH(IT,5)).NE.0)
1242      &        THEN
1243             NTRPART=NTRPART+1
1244             IH = ITTROUGH(IT,NBSTATION)
1245             IF (ITYP(IH).EQ.5.OR.ITYP(IH).EQ.6) THEN
1246                ISTAK = ID(IH)
1247                ISTAK = MOD(ISTAK,30000)
1248                ISTAK = MOD(ISTAK,10000)
1249 *     test
1250                pt=sqrt(pvert1(ih)**2+pvert2(ih)**2)
1251                thet=datan2(pt,pvert3(ih))*180/3.1416
1252                pp=sqrt(pt**2+pvert3(ih)**2)
1253
1254 *     print*,'istak=',istak
1255
1256                IF (ISTAK.EQ.IDRES) THEN ! psi or upsilon
1257                   print*,'resonance'
1258                   NTRMU = NTRMU+1
1259                   NTRMUALL = NTRMUALL+1
1260                   ITCHECK(IT) = 1
1261                ENDIF
1262             ENDIF
1263          ENDIF
1264       ENDDO
1265
1266       IF (IDEBUG.GE.2) THEN
1267          PRINT *,'RECO_TRACKF nb of part. GEANT crossing 5 st.:',
1268      &        NTRPART
1269          PRINT *,'RECO_TRACKF nb of muons/res. GEANT crossing 5 st.:',
1270      &        NTRMU
1271       ENDIF
1272
1273 **         CALL H_ACCEPTANCE(5)
1274 **         CALL H_ACCEPTANCE(4)
1275 **         PAUSE
1276
1277       NCAN = 0
1278
1279 *  Recherche 5 -> 4
1280
1281       CALL ORDONNE_HIT(5,HCUT)
1282
1283       DO IH = 1,INDEXMAX
1284          JJ = INDEXTAB(IH)
1285          X1 = XM(5,JJ)-(ZPLANE(5)-ZPLANE(4))*TAN(PHM(5,JJ))
1286          Y1 = YM(5,JJ)+(ZPLANE(5)-ZPLANE(4))*TAN(ALM(5,JJ))
1287      &                                            /COS(PHM(5,JJ))
1288          X2 = XM(5,JJ)-(ZPLANE(5)-ZPLANE(4)+DZ_PL(4))*TAN(PHM(5,JJ))
1289          Y2 = YM(5,JJ)+(ZPLANE(5)-ZPLANE(4)+DZ_PL(4))*TAN(ALM(5,JJ))
1290      &                                            /COS(PHM(5,JJ))
1291 *  Domaine de recherche dans la st. 4
1292          HCONST = 0.0136*SQRT(0.06)/HTOP  ! -> X/X0=0.03 % / chamber
1293          EPH2 = 2.0*PHIPREC**2  + (HCONST*HHIT(JJ))**2
1294          EAL2 = 2.0*ALAMPREC**2 + (HCONST*HHIT(JJ))**2
1295 ** 29.07
1296          EXM2 = 2.0*(XPREC/SQRT(2.))**2
1297      &          + (ZPLANE(5)-ZPLANE(4))**2/2.*EPH2
1298          EYM2 = 2.0*YPREC**2    + (ZPLANE(5)-ZPLANE(4))**2/2.*EAL2
1299          P1 = PTOT(IP(5,JJ))
1300          CALL CHFILL2(2500,SNGL(P1),SNGL(HHIT(JJ)),1.)
1301          CALL CHFILL2(2501,SNGL(P1),SNGL(HHIT(JJ)**2),1.)
1302          CALL CHFILL2(2502,SNGL(P1),SNGL(EPH2),1.)
1303          CALL CHFILL2(2503,SNGL(P1),SNGL(EAL2),1.)
1304          CALL CHFILL2(2504,SNGL(P1),SNGL(EXM2),1.)
1305          CALL CHFILL2(2505,SNGL(P1),SNGL(EYM2),1.)
1306
1307          CALL CHFILL2(2507,SNGL(P1),SNGL(SQRT(EPH2)),1.)
1308          CALL CHFILL2(2508,SNGL(P1),SNGL(SQRT(EAL2)),1.)
1309          CALL CHFILL2(2509,SNGL(P1),SNGL(SQRT(EXM2)),1.)
1310          CALL CHFILL2(2510,SNGL(P1),SNGL(SQRT(EYM2)),1.)
1311 ** end 29.07
1312          EXM12 = (PHIPREC**2+(HCONST*HHIT(JJ))**2)*
1313      &        (ZPLANE(5)-ZPLANE(4))**2 + XPREC**2
1314          EYM12 = (ALAMPREC**2+(HCONST*HHIT(JJ))**2)*
1315      &        (ZPLANE(5)-ZPLANE(4))**2 + YPREC**2
1316          EXM22 = (PHIPREC**2+(HCONST*HHIT(JJ))**2)*
1317      &        (ZPLANE(5)-ZPLANE(4)+DZ_PL(4))**2 + XPREC**2
1318          EYM22 = (ALAMPREC**2+(HCONST*HHIT(JJ))**2)*
1319      &        (ZPLANE(5)-ZPLANE(4)+DZ_PL(4))**2 + YPREC**2
1320
1321          EPH = SIGCUT*SQRT(EPH2)
1322          EAL = SIGCUT*SQRT(EAL2)
1323          EX1 = SIGCUT*SQRT(EXM12)
1324          EY1 = SIGCUT*SQRT(EYM12)
1325          EX2 = SIGCUT*SQRT(EXM22)
1326          EY2 = SIGCUT*SQRT(EYM22)
1327 ** 29.07
1328          EXM = SIGCUT*SQRT(EXM2)
1329          EYM = SIGCUT*SQRT(EYM2)
1330          CALL CHFILL2(2511,SNGL(P1),SNGL(EPH),1.)
1331          CALL CHFILL2(2512,SNGL(P1),SNGL(EAL),1.)
1332          CALL CHFILL2(2513,SNGL(P1),SNGL(EXM),1.)
1333          CALL CHFILL2(2514,SNGL(P1),SNGL(EYM),1.)
1334 ** end 29.07
1335 **         P2 = (HTOP/HHIT(JJ))**2
1336
1337 **         EPH=SIGCUT*SQRT(9.14D-7+1.2D-3/P2)
1338 **         EAL=SIGCUT*SQRT(1.84D-4)
1339 **         EX1=SIGCUT*SQRT(1.95D-2+6.37/P2)
1340 **         EY1=SIGCUT*SQRT(3.89+151./P2)
1341 **         EX2=EX1
1342 **         EY2=EY2
1343
1344 * renvoie le num de hit de 4 le plus pres dans le domaine de recherche
1345
1346          CALL DISTMIN4(X1,Y1,PHM(5,JJ),ALM(5,JJ),4,EX1,EY1,EPH,EAL,
1347      &        IFIND)
1348          P1 = PTOT(IP(5,JJ))
1349          CALL CHECK_HISTO4(11,4,IFIND,5,JJ,X1,Y1,PHM(5,JJ),ALM(5,JJ),P1,
1350      &                     EX1,EY1,EPH,EAL)
1351          IF (IFIND.GT.0) THEN
1352             CALL STOCK_CANDIDAT(5,JJ,4,IFIND,EX1,EY1,EPH,EAL,NCAN,1)
1353          ELSE
1354             CALL DISTMIN2(X1,Y1,X2,Y2,4,EX1,EY1,EX2,EY2,IFIND)
1355             CALL CHECK_HISTO2(0,4,IFIND,5,JJ,X1,Y1,X2,Y2,P1,EX1,EY1,
1356      &           EX2,EY2)
1357             IF (IFIND.GT.0) THEN
1358                CALL STOCK_CANDIDAT(5,JJ,4,IFIND,EX1,EY1,EPH,EAL,NCAN,2)
1359             ENDIF
1360          ENDIF
1361       ENDDO
1362
1363 *   Recherche 4 -> 5
1364
1365       CALL ORDONNE_HIT(4,HCUT)
1366
1367       DO IH = 1,INDEXMAX
1368          JJ = INDEXTAB(IH)
1369          X1 = XM(4,JJ)-(ZPLANE(4)-ZPLANE(5))*TAN(PHM(4,JJ))
1370          Y1 = YM(4,JJ)+(ZPLANE(4)-ZPLANE(5))*TAN(ALM(4,JJ))
1371      &                                            /COS(PHM(4,JJ))
1372          X2 = XM(4,JJ)-(ZPLANE(4)-ZPLANE(5)+DZ_PL(5))*TAN(PHM(4,JJ))
1373          Y2 = YM(4,JJ)+(ZPLANE(4)-ZPLANE(5)+DZ_PL(5))*TAN(ALM(4,JJ))
1374      &                                            /COS(PHM(4,JJ))
1375 *   Domaine de recherche dans la st. 5
1376          HCONST = 0.0136*SQRT(0.06)/HTOP  ! -> X/X0=0.03 / chamber
1377 **         EPH2 = 2.0*PHIPREC**2  + (HCONST*HHIT(JJ))**2
1378 **         EAL2 = 2.0*ALAMPREC**2 + (HCONST*HHIT(JJ))**2
1379 **         EX12 = 2.0*(XPREC/SQRT(2.))**2
1380 **     &          + (ZPLANE(5)-ZPLANE(4))**2/2.*EPH2
1381 **         EY12 = 2.0*YPREC**2 + (ZPLANE(5)-ZPLANE(4))**2/2.*EAL2
1382 **         EX22 = 2.0*(XPREC/SQRT(2.))**2
1383 **     &          + (ZPLANE(5)-ZPLANE(4)-DZ_PL(5))**2/2.*EPH2
1384 **         EY22 = 2.0*YPREC**2 + (ZPLANE(5)-ZPLANE(4)-DZ_PL(5))**2/2.
1385 **     &        *EAL2
1386          EPH2 = 2.0*PHIPREC**2  + (HCONST*HHIT(JJ))**2
1387          EAL2 = 2.0*ALAMPREC**2 + (HCONST*HHIT(JJ))**2
1388 ** 29.07
1389          EXM2 = 2.0*(XPREC/SQRT(2.))**2
1390      &          + (ZPLANE(5)-ZPLANE(4))**2/2.*EPH2
1391          EYM2 = 2.0*YPREC**2    + (ZPLANE(5)-ZPLANE(4))**2/2.*EAL2
1392          P1 = PTOT(IP(4,JJ))
1393          CALL CHFILL2(2400,SNGL(P1),SNGL(HHIT(JJ)),1.)
1394          CALL CHFILL2(2401,SNGL(P1),SNGL(HHIT(JJ)**2),1.)
1395          CALL CHFILL2(2402,SNGL(P1),SNGL(EPH2),1.)
1396          CALL CHFILL2(2403,SNGL(P1),SNGL(EAL2),1.)
1397          CALL CHFILL2(2404,SNGL(P1),SNGL(EXM2),1.)
1398          CALL CHFILL2(2405,SNGL(P1),SNGL(EYM2),1.)
1399
1400          CALL CHFILL2(2407,SNGL(P1),SNGL(SQRT(EPH2)),1.)
1401          CALL CHFILL2(2408,SNGL(P1),SNGL(SQRT(EAL2)),1.)
1402          CALL CHFILL2(2409,SNGL(P1),SNGL(SQRT(EXM2)),1.)
1403          CALL CHFILL2(2410,SNGL(P1),SNGL(SQRT(EYM2)),1.)
1404 ** end 29.07
1405
1406          EXM12 = (PHIPREC**2+(HCONST*HHIT(JJ))**2)*
1407      &        (ZPLANE(5)-ZPLANE(4))**2 + XPREC**2
1408          EYM12 = (ALAMPREC**2+(HCONST*HHIT(JJ))**2)*
1409      &        (ZPLANE(5)-ZPLANE(4))**2 + YPREC**2
1410          EXM22 = (PHIPREC**2+(HCONST*HHIT(JJ))**2)*
1411      &        (ZPLANE(5)-ZPLANE(4)-DZ_PL(5))**2 + XPREC**2
1412          EYM22 = (ALAMPREC**2+(HCONST*HHIT(JJ))**2)*
1413      &        (ZPLANE(5)-ZPLANE(4)-DZ_PL(5))**2 + YPREC**2
1414
1415          EPH = SIGCUT*SQRT(EPH2)
1416          EAL = SIGCUT*SQRT(EAL2)
1417          EX1 = SIGCUT*SQRT(EXM12)
1418          EY1 = SIGCUT*SQRT(EYM12)
1419          EX2 = SIGCUT*SQRT(EXM22)
1420          EY2 = SIGCUT*SQRT(EYM22)
1421 ** 29.07
1422          EXM = SIGCUT*SQRT(EXM2)
1423          EYM = SIGCUT*SQRT(EYM2)
1424          CALL CHFILL2(2411,SNGL(P1),SNGL(EPH),1.)
1425          CALL CHFILL2(2412,SNGL(P1),SNGL(EAL),1.)
1426          CALL CHFILL2(2413,SNGL(P1),SNGL(EXM),1.)
1427          CALL CHFILL2(2414,SNGL(P1),SNGL(EYM),1.)
1428 ** end 29.07
1429
1430 **         P2 = (HTOP/HHIT(JJ))**2
1431
1432 **         EPH=SIGCUT*SQRT(9.14D-7+1.2D-3/P2)
1433 **         EAL=SIGCUT*SQRT(1.84D-4)
1434 **         EX1=SIGCUT*SQRT(1.95D-2+6.37/P2)
1435 **         EY1=SIGCUT*SQRT(3.89+151./P2)
1436 **         EX2=EX1
1437 **         EY2=EY2
1438 *   Renvoie le num de hit de 5 le plus pres dans le domaine de recherche
1439          CALL DISTMIN2(X1,Y1,X2,Y2,5,EX1,EY1,EX2,EY2,IFIND)
1440          P1 = PTOT(IP(4,JJ))
1441          CALL CHECK_HISTO2(0,5,IFIND,4,JJ,X1,Y1,X2,Y2,P1,EX1,EY1,
1442      &        EX2,EY2)
1443          IF (IFIND.GT.0) THEN
1444             DO ICAN=1,NCAN
1445                IF (IFIND.EQ.JCAN(5,ICAN).AND.JJ.EQ.JCAN(4,ICAN)) GOTO 40
1446                IF (ABS(XM(5,JCAN(5,ICAN))-XM(5,IFIND)).LT.XPREC/10.)
1447      &            GO TO 40 ! elimine les doubles comptages de traces ccc
1448             ENDDO
1449             CALL STOCK_CANDIDAT(4,JJ,5,IFIND,EX1,EY1,EPH,EAL,NCAN,3)
1450          ENDIF
1451 40       CONTINUE
1452       ENDDO
1453
1454
1455       NMU45 = 0
1456       DO ICAN = 1,NCAN
1457          JJ4 = JCAN(4,ICAN)
1458          JJ5 = JCAN(5,ICAN)
1459          IF (JJ4.GT.0.AND.JJ5.GT.0) THEN
1460             ID4 = ID(IP(4,JJ4))
1461             ID5 = ID(IP(5,JJ5))
1462             IT = ITRACK(IP(5,JJ5))
1463             IF (ITCHECK(IT).EQ.1) THEN
1464                IF (ID4.EQ.ID5) THEN
1465                   NMU45 = NMU45 + 1
1466                ENDIF
1467             ENDIF
1468          ENDIF
1469       ENDDO
1470       IF (NMU45.GE.2) NRES(1) = NRES(1)+1
1471
1472       IF (IDEBUG.GE.2) THEN
1473          PRINT *,'RECO_TRACKF nb candidat recherche 4->5 et 5->4 :'
1474      &        ,NCAN
1475          PRINT *,'RECO_TRACKF nb of good muons 4->5 et 5->4 :'
1476      &        ,NMU45
1477
1478       ENDIF
1479
1480       NMU345 = 0
1481 *
1482 * -- Boucle sur les candidats (4,5) NCAN
1483 *
1484       DO I = 1,NBSTATION
1485          NERR(I) = 0
1486       ENDDO
1487       NMUONF = 0
1488       NGHOST = 0
1489       NTRACKF = 0
1490
1491       DO ICAN = 1,NCAN
1492          JJ1 = 0
1493          JJ2 = 0
1494          JJ3 = 0
1495          DO ICH = 1,NBSTATION
1496             MPOS(ICH) = 0
1497             MANG(ICH) = 0
1498          ENDDO
1499          MPOS(4) = 1
1500          MPOS(5) = 1
1501          MANG(4) = 1
1502          MANG(5) = 1
1503          IF (JCANTYP(ICAN).EQ.2) MANG(4)=0
1504          IF (JCANTYP(ICAN).EQ.3) MANG(5)=0
1505          EEXM(5) = EEX(ICAN)/SIGCUT
1506          EEYM(5) = EEY(ICAN)/SIGCUT
1507          EEPH(5) = EEP(ICAN)/SIGCUT
1508          EEAL(5) = EEA(ICAN)/SIGCUT
1509          EEXM(4) = EEX(ICAN)/SIGCUT
1510          EEYM(4) = EEY(ICAN)/SIGCUT
1511          EEPH(4) = EEP(ICAN)/SIGCUT
1512          EEAL(4) = EEA(ICAN)/SIGCUT
1513          JJ5 = JCAN(5,ICAN)
1514          JJ4 = JCAN(4,ICAN)
1515          P = PTOT(IP(5,JJ5))
1516          IF (IZM(4,JJ4).EQ.1) THEN
1517             ZPL4 = ZPLANE(4)
1518          ELSE
1519             ZPL4 = ZPLANE(4)-DZ_PL(4)
1520          ENDIF
1521          IF (IZM(5,JJ5).EQ.1) THEN
1522             ZPL5 = ZPLANE(5)
1523          ELSE
1524             ZPL5 = ZPLANE(5)-DZ_PL(5)
1525          ENDIF
1526          TPHEST = (XM(5,JJ5) - XM(4,JJ4))/(ZPL5-ZPL4)
1527          PHEST = ATAN(TPHEST)
1528          TALEST = -(YM(5,JJ5) - YM(4,JJ4))*COS(PHEST)
1529      &               /(ZPL5-ZPL4)
1530          PXZEST = -HTOP/(XM(4,JJ4) - ZPL4*TPHEST)
1531          PEST(1) = 1.0/PXZEST
1532          PEST(2) = TPHEST - ALPHATOP/PXZEST ! PHI emission au vertex
1533          PEST(3) = TALEST       ! tan(lambda)     !h=zm*ang deviation alpha
1534          PEST(4) = 0.0          !alpha=qbl/pxz
1535          PEST(5) = 0.0
1536          PSTEP(1) = 0.003       ! =d(1/p)=delta(p)/p**2
1537          PSTEP(2) = 0.001       ! 0.5 degre
1538          PSTEP(3) = 0.001       ! 0.5 degre
1539          PSTEP(4) = 0.0
1540          PSTEP(5) = 1.0
1541          XMES(4) = XM(4,JJ4)
1542          YMES(4) = YM(4,JJ4)
1543          IZMES(4) = IZM(4,JJ4)
1544          PHMES(4) = PHM(4,JJ4)
1545          ALMES(4) = ALM(4,JJ4)
1546          XMES(5) = XM(5,JJ5)
1547          YMES(5) = YM(5,JJ5)
1548          IZMES(5) = IZM(5,JJ5)
1549          PHMES(5) = PHM(5,JJ5)
1550          ALMES(5) = ALM(5,JJ5)
1551
1552          IVERTEX = 0 ! Vertex = (0,0,0)
1553
1554 * -- Fit 4,5,V pour la recherche dans 3
1555
1556          CALL TRACKF_FIT(IVERTEX,PEST,PSTEP,PXZINV45,TPHI45,TALAM45,
1557      &        XVERT,YVERT)
1558
1559 * -- Recherche dans la station 3
1560
1561          P2 = (1.0D0 + TALAM45**2)/(PXZINV45**2) ! P**2
1562 **         EPH=SIGCUT*SQRT(3.6D-6+0.011/P2)
1563 **         EAL=SIGCUT*SQRT(6.85D-4)
1564 **         EXM=SIGCUT*SQRT(0.034+446./P2)
1565 **         EYM=SIGCUT*SQRT(0.049+354./P2)
1566          EPH=SIGCUT*SQRT(4.1D-7+0.015/P2)
1567          EAL=SIGCUT*SQRT(1.04D-4)
1568          EXM=SIGCUT*SQRT(0.0+459./P2)
1569          EYM=SIGCUT*SQRT(0.042+345./P2)
1570          EEXM(3) = EXM/SIGCUT
1571          EEYM(3) = EYM/SIGCUT
1572          EEPH(3) = EPH/SIGCUT
1573          EEAL(3) = EAL/SIGCUT
1574 ** 29.07
1575          CALL CHF1(2301,SNGL(SQRT(P2)),1.)
1576          P1 = PTOT(IP(5,JJ5))
1577          CALL CHFILL2(2311,SNGL(P1),SNGL(EPH),1.)
1578          CALL CHFILL2(2312,SNGL(P1),SNGL(EAL),1.)
1579          CALL CHFILL2(2313,SNGL(P1),SNGL(EXM),1.)
1580          CALL CHFILL2(2314,SNGL(P1),SNGL(EYM),1.)
1581
1582          CALL CHFILL2(2315,SNGL(SQRT(P2)),SNGL(EPH),1.)
1583          CALL CHFILL2(2316,SNGL(SQRT(P2)),SNGL(EAL),1.)
1584          CALL CHFILL2(2317,SNGL(SQRT(P2)),SNGL(EXM),1.)
1585          CALL CHFILL2(2318,SNGL(SQRT(P2)),SNGL(EYM),1.)
1586
1587          CALL CHFILL2(2302,SNGL(SQRT(P2)),SNGL((EPH/SIGCUT)**2),1.)
1588          CALL CHFILL2(2303,SNGL(SQRT(P2)),SNGL((EAL/SIGCUT)**2),1.)
1589          CALL CHFILL2(2304,SNGL(SQRT(P2)),SNGL((EXM/SIGCUT)**2),1.)
1590          CALL CHFILL2(2305,SNGL(SQRT(P2)),SNGL((EYM/SIGCUT)**2),1.)
1591
1592          CALL CHFILL2(2307,SNGL(P1),SNGL((EPH/SIGCUT)**2),1.)
1593          CALL CHFILL2(2308,SNGL(P1),SNGL((EAL/SIGCUT)**2),1.)
1594          CALL CHFILL2(2309,SNGL(P1),SNGL((EXM/SIGCUT)**2),1.)
1595          CALL CHFILL2(2310,SNGL(P1),SNGL((EYM/SIGCUT)**2),1.)
1596 ** end 29.07
1597
1598 **         DO IEX = 4,5
1599 **            EEXM(IEX) = EEXM(3)
1600 **            EEYM(IEX) = EEYM(3)
1601 **            EEPH(IEX) = EEPH(3)
1602 **            EEAL(IEX) = EEAL(3)
1603 **         ENDDO
1604          X1 = XPL(3,1)
1605          Y1 = YPL(3,1)
1606          X2 = XPL(3,2)
1607          Y2 = YPL(3,2)
1608          PHI1 = PHPL(3)
1609          ALAM1 = ALPL(3)
1610
1611          CALL DISTMIN4(X1,Y1,PHI1,ALAM1,3,EXM,EYM,EPH,EAL,IPA3)
1612
1613          P1 = PTOT(IP(5,JJ5))
1614
1615          CALL CHECK_HISTO4(61,3,IPA3,5,JJ5,X1,Y1,PHI1,ALAM1,P1,
1616      &        EXM,EYM,EPH,EAL)
1617
1618          IF (IPA3.NE.0) THEN
1619             JJ3 = IPA3
1620             XMES(3) = XM(3,JJ3)
1621             YMES(3) = YM(3,JJ3)
1622             IZMES(3) = IZM(3,JJ3)
1623             PHMES(3) = PHM(3,JJ3)
1624             ALMES(3) = ALM(3,JJ3)
1625             MPOS(3) = 1
1626             MANG(3) = 1
1627             GO TO 124
1628          ELSE
1629             CALL DISTMIN2(X1,Y1,X2,Y2,3,EXM,EYM,0.D0,0.D0,IP3)
1630             CALL CHECK_HISTO2(0,3,IP3,5,JJ5,X1,Y1,X2,Y2,P1,EXM,EYM,
1631      &           0.D0,0.D0)
1632          ENDIF
1633          IF (IP3.NE.0) THEN
1634             JJ3 = IP3
1635             XMES(3) = XM(3,JJ3)
1636             YMES(3) = YM(3,JJ3)
1637             IZMES(3) = IZM(3,JJ3)
1638             MPOS(3) = 1
1639             MANG(3) = 0
1640          ELSE
1641             GO TO 123
1642          ENDIF
1643
1644
1645 *     -- Fit 3,4,5 pour la recherche dans 2
1646
1647  124     CONTINUE
1648
1649          IF (JJ5.GT.0.AND.JJ4.GT.0.AND.JJ3.GT.0) THEN
1650             ID4 = ID(IP(4,JJ4))
1651             ID5 = ID(IP(5,JJ5))
1652             ID3 = ID(IP(3,JJ3))
1653             IT = ITRACK(IP(5,JJ5))
1654             IF (ITCHECK(IT).EQ.1) THEN
1655                IF ((ID5.EQ.ID3).AND.(ID5.EQ.ID4)) THEN
1656                   NMU345 = NMU345 + 1
1657                ENDIF
1658             ENDIF
1659          ENDIF
1660
1661          IVERTEX = 0
1662
1663          PEST(1) = PXZINV45
1664          PEST(2) = TPHI45
1665          PEST(3) = TALAM45
1666          PEST(4) = 0.
1667          PEST(5) = 0.
1668          PSTEP(1) = 0.003
1669          PSTEP(2) = 0.001
1670          PSTEP(3) = 0.001
1671          PSTEP(4) = 1.
1672          PSTEP(5) = 1.
1673
1674          CALL TRACKF_FIT(IVERTEX,PEST,PSTEP,PXZINV345,TPHI345,
1675      &        TALAM345,XVERT,YVERT)
1676
1677 *     -- Recherche dans la st. 2
1678
1679          P2 = (1.0D0 + TALAM345**2)/(PXZINV345**2) ! P**2
1680 **         EPH=SIGCUT*SQRT(3.63D-6+4.8D-3/P2)
1681 **         EAL=SIGCUT*SQRT(6.49D-4)
1682 **         EXM=SIGCUT*SQRT(1.64D-2+821./P2)
1683 **         EYM=SIGCUT*SQRT(4.83D-2+866./P2)
1684          EPH=SIGCUT*SQRT(5.78D-7+5.97D-3/P2)
1685          EAL=SIGCUT*SQRT(1.D-4)
1686          EXM=SIGCUT*SQRT(0.+1453./P2)
1687          EYM=SIGCUT*SQRT(2.78D-2+504./P2)
1688          EEXM(2) = EXM/SIGCUT
1689          EEYM(2) = EYM/SIGCUT
1690          EEPH(2) = EPH/SIGCUT
1691          EEAL(2) = EAL/SIGCUT
1692 ** 29.07
1693          CALL CHF1(2201,SNGL(SQRT(P2)),1.)
1694          P1 = PTOT(IP(5,JJ5))
1695          CALL CHFILL2(2211,SNGL(P1),SNGL(EPH),1.)
1696          CALL CHFILL2(2212,SNGL(P1),SNGL(EAL),1.)
1697          CALL CHFILL2(2213,SNGL(P1),SNGL(EXM),1.)
1698          CALL CHFILL2(2214,SNGL(P1),SNGL(EYM),1.)
1699
1700          CALL CHFILL2(2215,SNGL(SQRT(P2)),SNGL(EPH),1.)
1701          CALL CHFILL2(2216,SNGL(SQRT(P2)),SNGL(EAL),1.)
1702          CALL CHFILL2(2217,SNGL(SQRT(P2)),SNGL(EXM),1.)
1703          CALL CHFILL2(2218,SNGL(SQRT(P2)),SNGL(EYM),1.)
1704
1705          CALL CHFILL2(2202,SNGL(SQRT(P2)),SNGL((EPH/SIGCUT)**2),1.)
1706          CALL CHFILL2(2203,SNGL(SQRT(P2)),SNGL((EAL/SIGCUT)**2),1.)
1707          CALL CHFILL2(2204,SNGL(SQRT(P2)),SNGL((EXM/SIGCUT)**2),1.)
1708          CALL CHFILL2(2205,SNGL(SQRT(P2)),SNGL((EYM/SIGCUT)**2),1.)
1709
1710          CALL CHFILL2(2207,SNGL(P1),SNGL((EPH/SIGCUT)**2),1.)
1711          CALL CHFILL2(2208,SNGL(P1),SNGL((EAL/SIGCUT)**2),1.)
1712          CALL CHFILL2(2209,SNGL(P1),SNGL((EXM/SIGCUT)**2),1.)
1713          CALL CHFILL2(2210,SNGL(P1),SNGL((EYM/SIGCUT)**2),1.)
1714 ** end 29.07
1715
1716 **         DO IEX = 3,5
1717 **            EEXM(IEX) = EEXM(2)
1718 **            EEYM(IEX) = EEYM(2)
1719 **            EEPH(IEX) = EEPH(2)
1720 **            EEAL(IEX) = EEAL(2)
1721 **         ENDDO
1722
1723          X1 = XPL(2,1)
1724          Y1 = YPL(2,1)
1725          PHI1 = PHPL(2)
1726          ALAM1 = ALPL(2)
1727          CALL DISTMIN4(X1,Y1,PHI1,ALAM1,2,EXM,EYM,EPH,EAL,IPA2)
1728          P1 = PTOT(IP(5,JJ5))
1729          CALL CHECK_HISTO4(21,2,IPA2,5,JJ5,X1,Y1,PHI1,ALAM1,P1,
1730      &        EXM,EYM,EPH,EAL)
1731 *     -- Recherche dans la st. 1
1732 **         EPH=SIGCUT*SQRT(3.54D-6+4.49D-3/P2)
1733 **         EAL=SIGCUT*SQRT(6.14D-4)
1734 **         EXM=SIGCUT*SQRT(6.43D-2+986./P2)
1735 **         EYM=SIGCUT*SQRT(4.66D-2+1444./P2)
1736          EPH=SIGCUT*SQRT(4.96D-7+5.87D-3/P2)
1737          EAL=SIGCUT*SQRT(1.D-4)
1738          EXM=SIGCUT*SQRT(6.98D-4+1467./P2)
1739          EYM=SIGCUT*SQRT(5.22D-2+1013./P2)
1740          EEXM(1) = EXM/SIGCUT
1741          EEYM(1) = EYM/SIGCUT
1742          EEPH(1) = EPH/SIGCUT
1743          EEAL(1) = EAL/SIGCUT
1744 ** 29.07
1745          CALL CHFILL2(2111,SNGL(P1),SNGL(EPH),1.)
1746          CALL CHFILL2(2112,SNGL(P1),SNGL(EAL),1.)
1747          CALL CHFILL2(2113,SNGL(P1),SNGL(EXM),1.)
1748          CALL CHFILL2(2114,SNGL(P1),SNGL(EYM),1.)
1749
1750          CALL CHFILL2(2115,SNGL(SQRT(P2)),SNGL(EPH),1.)
1751          CALL CHFILL2(2116,SNGL(SQRT(P2)),SNGL(EAL),1.)
1752          CALL CHFILL2(2117,SNGL(SQRT(P2)),SNGL(EXM),1.)
1753          CALL CHFILL2(2118,SNGL(SQRT(P2)),SNGL(EYM),1.)
1754
1755          CALL CHFILL2(2102,SNGL(SQRT(P2)),SNGL((EPH/SIGCUT)**2),1.)
1756          CALL CHFILL2(2103,SNGL(SQRT(P2)),SNGL((EAL/SIGCUT)**2),1.)
1757          CALL CHFILL2(2104,SNGL(SQRT(P2)),SNGL((EXM/SIGCUT)**2),1.)
1758          CALL CHFILL2(2105,SNGL(SQRT(P2)),SNGL((EYM/SIGCUT)**2),1.)
1759
1760          CALL CHFILL2(2107,SNGL(P1),SNGL((EPH/SIGCUT)**2),1.)
1761          CALL CHFILL2(2108,SNGL(P1),SNGL((EAL/SIGCUT)**2),1.)
1762          CALL CHFILL2(2109,SNGL(P1),SNGL((EXM/SIGCUT)**2),1.)
1763          CALL CHFILL2(2110,SNGL(P1),SNGL((EYM/SIGCUT)**2),1.)
1764 ** end 29.07
1765
1766
1767 **         DO IEX = 2,5
1768 **            EEXM(IEX) = EEXM(1)
1769 **            EEYM(IEX) = EEYM(1)
1770 **            EEPH(IEX) = EEPH(1)
1771 **            EEAL(IEX) = EEAL(1)
1772 **         ENDDO
1773
1774          X1 = XPL(1,1)
1775          Y1 = YPL(1,1)
1776          PHI1 = PHPL(1)
1777          ALAM1 = ALPL(1)
1778          CALL DISTMIN4(X1,Y1,PHI1,ALAM1,1,EXM,EYM,EPH,EAL,IPA1)
1779          CALL CHECK_HISTO4(31,1,IPA1,5,JJ5,X1,Y1,PHI1,ALAM1,P1,
1780      &        EXM,EYM,EPH,EAL)
1781 *     -- A partir de P+A dans la st. 2
1782          IPA2PA1 = 0
1783          IF (IPA2.GT.0) THEN
1784             PEST(1) = PXZINV345
1785             PEST(2) = TPHI345
1786             PEST(3) = TALAM345
1787             PEST(4) = 0.
1788             PEST(5) = 0.
1789             PSTEP(1) = 0.003
1790             PSTEP(2) = 0.001
1791             PSTEP(3) = 0.001
1792             PSTEP(4) = 1.
1793             PSTEP(5) = 1.
1794             XMES(2) = XM(2,IPA2)
1795             YMES(2) = YM(2,IPA2)
1796             IZMES(2) = IZM(2,IPA2)
1797             PHMES(2) = PHM(2,IPA2)
1798             ALMES(2) = ALM(2,IPA2)
1799             MPOS(2) = 1
1800             MANG(2) = 1
1801             IVERTEX = 1
1802 *     -- Fit V,2,3,4,5 pour la recherche dans 1
1803
1804             CALL TRACKF_FIT(IVERTEX,PEST,PSTEP,PXZINV,TPHI,TALAM,
1805      &           XVERT,YVERT)
1806
1807 *     !!! ATTENTION aux erreurs
1808
1809             P2 = (1.0D0 + TALAM**2)/(PXZINV**2)
1810 **            EPH=SIGCUT*SQRT(3.15D-6+9.21D-3/P2)
1811 **            EAL=SIGCUT*SQRT(5.94D-4)
1812 **            EXM=SIGCUT*SQRT(4.16D-2+182./P2)
1813 **            EYM=SIGCUT*SQRT(2.13)
1814             EPH=SIGCUT*SQRT(9.58D-7+8.93D-3/P2)
1815             EAL=SIGCUT*SQRT(1.D-4)
1816             EXM=SIGCUT*SQRT(1.92D-2+103.3/P2)
1817             EYM=SIGCUT*SQRT(6.3D-2+81.1/P2)
1818             EEXM(1) = EXM/SIGCUT
1819             EEYM(1) = EYM/SIGCUT
1820             EEPH(1) = EPH/SIGCUT
1821             EEAL(1) = EAL/SIGCUT
1822 ** 29.07
1823          CALL CHF1(2701,SNGL(SQRT(P2)),1.)
1824          CALL CHFILL2(2711,SNGL(P1),SNGL(EPH),1.)
1825          CALL CHFILL2(2712,SNGL(P1),SNGL(EAL),1.)
1826          CALL CHFILL2(2713,SNGL(P1),SNGL(EXM),1.)
1827          CALL CHFILL2(2714,SNGL(P1),SNGL(EYM),1.)
1828
1829          CALL CHFILL2(2715,SNGL(SQRT(P2)),SNGL(EPH),1.)
1830          CALL CHFILL2(2716,SNGL(SQRT(P2)),SNGL(EAL),1.)
1831          CALL CHFILL2(2717,SNGL(SQRT(P2)),SNGL(EXM),1.)
1832          CALL CHFILL2(2718,SNGL(SQRT(P2)),SNGL(EYM),1.)
1833
1834          CALL CHFILL2(2702,SNGL(SQRT(P2)),SNGL((EPH/SIGCUT)**2),1.)
1835          CALL CHFILL2(2703,SNGL(SQRT(P2)),SNGL((EAL/SIGCUT)**2),1.)
1836          CALL CHFILL2(2704,SNGL(SQRT(P2)),SNGL((EXM/SIGCUT)**2),1.)
1837          CALL CHFILL2(2705,SNGL(SQRT(P2)),SNGL((EYM/SIGCUT)**2),1.)
1838
1839
1840          CALL CHFILL2(2707,SNGL(P1),SNGL((EPH/SIGCUT)**2),1.)
1841          CALL CHFILL2(2708,SNGL(P1),SNGL((EAL/SIGCUT)**2),1.)
1842          CALL CHFILL2(2709,SNGL(P1),SNGL((EXM/SIGCUT)**2),1.)
1843          CALL CHFILL2(2710,SNGL(P1),SNGL((EYM/SIGCUT)**2),1.)
1844 ** end 29.07
1845
1846
1847 **            DO IEX = 2,5
1848 **               EEXM(IEX) = EEXM(1)
1849 **               EEYM(IEX) = EEYM(1)
1850 **               EEPH(IEX) = EEPH(1)
1851 **               EEAL(IEX) = EEAL(1)
1852 **            ENDDO
1853
1854             X1 = XPL(1,1)
1855             Y1 = YPL(1,1)
1856             X2 = XPL(1,2)
1857             Y2 = YPL(1,2)
1858             PHI1 = PHPL(1)
1859             ALAM1 = ALPL(1)
1860             CALL DISTMIN4(X1,Y1,PHI1,ALAM1,1,EXM,EYM,EPH,EAL,
1861      &           IPA2PA1)
1862             CALL CHECK_HISTO4(41,1,IPA2PA1,5,JJ5,X1,Y1,PHI1,ALAM1,
1863      &           P1,EXM,EYM,EPH,EAL)
1864             IF (IPA2PA1.GT.0) THEN
1865                JJ2 = IPA2
1866                JJ1 = IPA2PA1
1867                XMES(1) = XM(1,JJ1)
1868                YMES(1) = YM(1,JJ1)
1869                IZMES(1) = IZM(1,JJ1)
1870                PHMES(1) = PHM(1,JJ1)
1871                ALMES(1) = ALM(1,JJ1)
1872                MPOS(1) = 1
1873                MANG(1) = 1
1874                GOTO 123
1875             ELSE
1876                CALL DISTMIN2(X1,Y1,X2,Y2,1,EXM,EYM,0.D0,0.D0,IPA2P1)
1877                CALL CHECK_HISTO2(0,1,IPA2P1,5,JJ5,X1,Y1,X2,Y2,P1,
1878      &              EXM,EYM,0.D0,0.D0)
1879             ENDIF
1880          ENDIF
1881 *     --  A partir de P+A dans la st. 1
1882          IPA1PA2 = 0
1883          IF (IPA1.GT.0) THEN
1884             PEST(1) = PXZINV345
1885             PEST(2) = TPHI345
1886             PEST(3) = TALAM345
1887             PEST(4) = 0.
1888             PEST(5) = 0.
1889             PSTEP(1) = 0.003
1890             PSTEP(2) = 0.001
1891             PSTEP(3) = 0.001
1892             PSTEP(4) = 1.
1893             PSTEP(5) = 1.
1894             XMES(1) = XM(1,IPA1)
1895             YMES(1) = YM(1,IPA1)
1896             IZMES(1) = IZM(1,IPA1)
1897             PHMES(1) = PHM(1,IPA1)
1898             ALMES(1) = ALM(1,IPA1)
1899             MPOS(1) = 1
1900             MANG(1) = 1
1901             MPOS(2) = 0
1902             MANG(2) = 0
1903             IVERTEX = 1
1904 *     -- Fit V,1,3,4,5 pour la recherche dans 2
1905
1906             CALL TRACKF_FIT(IVERTEX,PEST,PSTEP,PXZINV,TPHI,TALAM,
1907      &           XVERT,YVERT)
1908
1909 *     !!! ATTENTION aux erreurs
1910
1911             P2 = (1.0D0 + TALAM**2)/(PXZINV**2)
1912 **            EPH=SIGCUT*SQRT(3.15D-6+9.21D-3/P2)
1913 **            EAL=SIGCUT*SQRT(5.94D-4)
1914 **            EXM=SIGCUT*SQRT(4.16D-2+182./P2)
1915 **            EYM=SIGCUT*SQRT(2.13)
1916             EPH=SIGCUT*SQRT(9.58D-7+8.93D-3/P2)
1917             EAL=SIGCUT*SQRT(1.D-4)
1918             EXM=SIGCUT*SQRT(4.0D-2+11.4/P2)
1919             EYM=SIGCUT*SQRT(5.5D-2+44.35/P2)
1920
1921             EEXM(2) = EXM/SIGCUT
1922             EEYM(2) = EYM/SIGCUT
1923             EEPH(2) = EPH/SIGCUT
1924             EEAL(2) = EAL/SIGCUT
1925 ** 29.07
1926          CALL CHF1(2801,SNGL(SQRT(P2)),1.)
1927          CALL CHFILL2(2811,SNGL(P1),SNGL(EPH),1.)
1928          CALL CHFILL2(2812,SNGL(P1),SNGL(EAL),1.)
1929          CALL CHFILL2(2813,SNGL(P1),SNGL(EXM),1.)
1930          CALL CHFILL2(2814,SNGL(P1),SNGL(EYM),1.)
1931
1932          CALL CHFILL2(2815,SNGL(SQRT(P2)),SNGL(EPH),1.)
1933          CALL CHFILL2(2816,SNGL(SQRT(P2)),SNGL(EAL),1.)
1934          CALL CHFILL2(2817,SNGL(SQRT(P2)),SNGL(EXM),1.)
1935          CALL CHFILL2(2818,SNGL(SQRT(P2)),SNGL(EYM),1.)
1936
1937          CALL CHFILL2(2802,SNGL(SQRT(P2)),SNGL((EPH/SIGCUT)**2),1.)
1938          CALL CHFILL2(2803,SNGL(SQRT(P2)),SNGL((EAL/SIGCUT)**2),1.)
1939          CALL CHFILL2(2804,SNGL(SQRT(P2)),SNGL((EXM/SIGCUT)**2),1.)
1940          CALL CHFILL2(2805,SNGL(SQRT(P2)),SNGL((EYM/SIGCUT)**2),1.)
1941
1942
1943          CALL CHFILL2(2807,SNGL(P1),SNGL((EPH/SIGCUT)**2),1.)
1944          CALL CHFILL2(2808,SNGL(P1),SNGL((EAL/SIGCUT)**2),1.)
1945          CALL CHFILL2(2809,SNGL(P1),SNGL((EXM/SIGCUT)**2),1.)
1946          CALL CHFILL2(2810,SNGL(P1),SNGL((EYM/SIGCUT)**2),1.)
1947
1948 ** end 29.07
1949
1950
1951
1952
1953 **            DO IEX = 1,5
1954 **               EEXM(IEX) = EEXM(2)
1955 **               EEYM(IEX) = EEYM(2)
1956 **               EEPH(IEX) = EEPH(2)
1957 **               EEAL(IEX) = EEAL(2)
1958 **            ENDDO
1959
1960             X1 = XPL(2,1)
1961             Y1 = YPL(2,1)
1962             X2 = XPL(2,2)
1963             Y2 = YPL(2,2)
1964             PHI1 = PHPL(2)
1965             ALAM1 = ALPL(2)
1966
1967             CALL DISTMIN4(X1,Y1,PHI1,ALAM1,2,EXM,EYM,EPH,EAL,
1968      &           IPA1PA2)
1969
1970             CALL CHECK_HISTO4(51,2,IPA1PA2,5,JJ5,X1,Y1,PHI1,ALAM1,
1971      &           P1,EXM,EYM,EPH,EAL)
1972             IF (IPA1PA2.GT.0) THEN
1973                JJ1 = IPA1
1974                JJ2 = IPA1PA2
1975                XMES(2) = XM(2,JJ2)
1976                YMES(2) = YM(2,JJ2)
1977                IZMES(2) = IZM(2,JJ2)
1978                PHMES(2) = PHM(2,JJ2)
1979                ALMES(2) = ALM(2,JJ2)
1980                MPOS(2) = 1
1981                MANG(2) = 1
1982                GOTO 123
1983             ELSE
1984                CALL DISTMIN2(X1,Y1,X2,Y2,2,EXM,EYM,0.D0,0.D0,IPA1P2)
1985                CALL CHECK_HISTO2(0,2,IPA1P2,5,JJ5,X1,Y1,X2,Y2,P1,
1986      &              EXM,EYM,0.D0,0.D0)
1987             ENDIF
1988          ENDIF
1989 *     -- Selection par type de candidat
1990          IF (IPA1PA2.EQ.0.OR.IPA2PA1.EQ.0) THEN
1991             IF (IPA2.GT.0.AND.IPA2P1.GT.0) THEN
1992                JJ2 = IPA2
1993                JJ1 = IPA2P1
1994                XMES(1) = XM(1,JJ1)
1995                YMES(1) = YM(1,JJ1)
1996                IZMES(1) = IZM(1,JJ1)
1997                MPOS(1) = 1
1998                MANG(1) = 0
1999                MPOS(2) = 1
2000                MANG(2) = 1
2001             ELSEIF(IPA1.GT.0.AND.IPA1P2.GT.0) THEN
2002                JJ1 = IPA1
2003                JJ2 = IPA1P2
2004                XMES(2) = XM(2,JJ2)
2005                YMES(2) = YM(2,JJ2)
2006                IZMES(2) = IZM(2,JJ2)
2007                MPOS(1) = 1
2008                MANG(1) = 1
2009                MPOS(2) = 1
2010                MANG(2) = 0
2011             ENDIF
2012          ENDIF
2013
2014  123     CONTINUE
2015 ***
2016          NTRACKFOLD = NTRACKF
2017          NMUONFOLD = NMUONF
2018          NGHOSTOLD = NGHOST
2019          CALL CHECK_MUON(JJ1,JJ2,JJ3,JJ4,JJ5,NTRACKF,NMUONF,
2020      &        NGHOST,NERR,PXZINV,TPHI,TALAM,XVERT,YVERT)
2021
2022          IF (NTRACKF.GT.NTRACKFOLD) THEN
2023             ISTAT(NTRACKF) = 0
2024             IF (NMUONF.GT.NMUONFOLD) ISTAT(NTRACKF) = 1 ! Good muon
2025             IF (NGHOST.GT.NGHOSTOLD) ISTAT(NTRACKF) = 2 ! ghost track
2026             PXZOUT(NTRACKF) = 1./PXZINV
2027             TPHIOUT(NTRACKF) = TPHI
2028             TALAMOUT(NTRACKF) = TALAM
2029             XVERTOUT(NTRACKF) = XVERT
2030             YVERTOUT(NTRACKF) = YVERT
2031             PXVOUT(NTRACKF) = PVERT1(IP(5,JJ5))
2032             PYVOUT(NTRACKF) = PVERT2(IP(5,JJ5))
2033             PZVOUT(NTRACKF) = PVERT3(IP(5,JJ5))
2034             CHI2OUT(NTRACKF) = CHI2PL
2035          ENDIF
2036 ***
2037       ENDDO                     ! end loop on candidats NCAN
2038
2039       IF (NMU345.GE.2) NRES(2) = NRES(2) + 1
2040
2041       IF (IDEBUG.GE.2) THEN
2042          PRINT *,'RECO_TRACKF nb of good muons 3 4 5 :'
2043      &        ,NMU345
2044
2045       ENDIF
2046
2047       IF (IDEBUG.GE.2) THEN
2048          PRINT *,'RECO_TRACKF nb of track /evt     :',NTRACKF
2049          PRINT *,'RECO_TRACKF nb of good muon /evt :',NMUONF
2050          PRINT *,'RECO_TRACKF nb of ghost /evt     :',NGHOST
2051          DO I = 1,4
2052             PRINT *,'RECO_TRACKF nb of error in st',I,'/evt:',NERR(I)
2053          ENDDO
2054       ENDIF
2055       IF (NTRMU.GE.2) NRES(5) = NRES(5)+1
2056       IF (NMUONF.GE.2) NRESF = NRESF+1
2057       NMUONALL = NMUONALL+NMUONF
2058       NGHOSTALL = NGHOSTALL+NGHOST
2059       NTRACKFALL = NTRACKFALL+NTRACKF
2060       DO I = 1,4
2061          NERRALL(I) = NERRALL(I)+NERR(I)
2062       ENDDO
2063
2064       CALL TRACKF_OUT(IR,NTRACKF)
2065 ***
2066       RETURN
2067       END
2068
2069 **********************************************************************
2070       SUBROUTINE TRACKF_OUT(IR,NTRACKF)
2071 **********************************************************************
2072 *
2073 *
2074 **********************************************************************
2075       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2076
2077       PARAMETER (MAXHITCH=10000,MAXTRK=50000,MAXHITTOT=20000,
2078      &           NBSTATION=5)
2079       PARAMETER (NBCHAMBER=10,NTRMAX=500)
2080
2081       COMMON/RECOUT/JJO(NBCHAMBER,NTRMAX),XMESO(NBCHAMBER,NTRMAX),
2082      &              YMESO(NBCHAMBER,NTRMAX)
2083
2084       COMMON/TRACKFOUT/IEVOUT,NTREVT,JJOUT(NBCHAMBER,NTRMAX),
2085      &               ISTAT(NTRMAX),PXZOUT(NTRMAX),TPHIOUT(NTRMAX),
2086      &               TALAMOUT(NTRMAX),XVERTOUT(NTRMAX),YVERTOUT(NTRMAX),
2087      &               CHI2OUT(NTRMAX),
2088      &               XMESOUT(NBCHAMBER,NTRMAX),YMESOUT(NBCHAMBER,NTRMAX)
2089      &              ,PXVOUT(NTRMAX),PYVOUT(NTRMAX),PZVOUT(NTRMAX)
2090 **
2091       IEVOUT = IR
2092       NTREVT = NTRACKF
2093       IF (NTREVT.GT.0) THEN
2094          DO ITR = 1,NTREVT
2095             ICH = 0
2096             DO IST = 1,NBSTATION
2097                DO ILOOP = 1,2
2098                   ICH = ICH + 1
2099                   JJOUT(ICH,ITR) = JJO(ICH,ITR)
2100                   IF (JJOUT(ICH,ITR).GT.0) THEN
2101                       XMESOUT(ICH,ITR) = XMESO(ICH,ITR)
2102                       YMESOUT(ICH,ITR) = YMESO(ICH,ITR)
2103                   ENDIF
2104                ENDDO
2105             ENDDO
2106          ENDDO
2107       ENDIF
2108
2109
2110       RETURN
2111       END
2112 **********************************************************************
2113       SUBROUTINE CHECK_MUON(JJ1,JJ2,JJ3,JJ4,JJ5,
2114      &                  NTRACKF,NMUONF,NGHOST,NERR,PXZINV,TPHI,TALAM,
2115      &                  XVERT,YVERT)
2116 **********************************************************************
2117 *   Check muon track candidate using GEANT informations
2118 *
2119 **********************************************************************
2120       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2121 * --
2122       PARAMETER (MAXHITCH=10000,MAXTRK=50000,MAXHITTOT=20000,
2123      &           NBSTATION=5)
2124       PARAMETER (NBCHAMBER=10,NTRMAX=500)
2125
2126       COMMON/CHHIT/XM(NBSTATION,MAXHITCH),YM(NBSTATION,MAXHITCH),
2127      &             PHM(NBSTATION,MAXHITCH),ALM(NBSTATION,MAXHITCH),
2128      &             IZM(NBSTATION,MAXHITCH),
2129      &             IP(NBSTATION,MAXHITCH),JHIT(NBSTATION),
2130      &             XMR(NBSTATION,MAXHITCH,2),YMR(NBSTATION,MAXHITCH,2)
2131
2132       COMMON/RHIT/ITYP(MAXHITTOT),XTR(MAXHITTOT),YTR(MAXHITTOT),
2133      &            PTOT(MAXHITTOT),ID(MAXHITTOT),IZST(MAXHITTOT),
2134      &            PVERT1(MAXHITTOT),PVERT2(MAXHITTOT),PVERT3(MAXHITTOT),
2135      &            ZVERT(MAXHITTOT),NHITTOT
2136
2137       COMMON/VERIFGEANT/ITTROUGH(MAXTRK,NBSTATION),
2138      &                  IT_LIST(MAXTRK),IT_NP(MAXTRK),ITCHECK(MAXTRK),
2139      &                  ITRACK(MAXHITTOT)
2140
2141       COMMON/MEASUR/XMES(NBSTATION),YMES(NBSTATION),IZMES(NBSTATION),
2142      &              PHMES(NBSTATION),ALMES(NBSTATION),MPOS(NBSTATION),
2143      &              MANG(NBSTATION)
2144
2145       COMMON/RECOUT/JJO(NBCHAMBER,NTRMAX),XMESO(NBCHAMBER,NTRMAX),
2146      &              YMESO(NBCHAMBER,NTRMAX)
2147
2148       COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
2149
2150       COMMON/MAGNET/BL3,B0
2151
2152       REAL*4 R1,R2
2153       DATA R1,R2/0.,1./
2154
2155       DIMENSION NERR(NBSTATION),JJK(NBSTATION)
2156
2157 *      print*,' *** appel de la subroutine check_muon ***'
2158
2159       IF (JJ1.GT.0.AND.JJ2.GT.0.AND.JJ3.GT.0.AND.JJ4.GT.0
2160      &     .AND.JJ5.GT.0) THEN
2161          IDFIND = ID(IP(5,JJ5))
2162          IT = ITRACK(IP(5,JJ5))
2163          JJK(1) = JJ1
2164          JJK(2) = JJ2
2165          JJK(3) = JJ3
2166          JJK(4) = JJ4
2167          JJK(5) = JJ5
2168          NTRACKF = NTRACKF+1
2169          DO I = 1,NBCHAMBER
2170             JJO(I,NTRACKF) = 0
2171          ENDDO
2172          ICH1 = -1
2173          DO IST = 1,5
2174             ICH1 = ICH1 + 2
2175             IF (XMR(IST,JJK(IST),1).LT.999.) THEN
2176                JJO(ICH1,NTRACKF) = JJK(IST)
2177                XMESO(ICH1,NTRACKF) = XMR(IST,JJK(IST),1)
2178                YMESO(ICH1,NTRACKF) = YMR(IST,JJK(IST),1)
2179                ICH2 = ICH1 + 1
2180                IF (MANG(IST).EQ.1) THEN
2181                   JJO(ICH2,NTRACKF) = JJK(IST)
2182                   XMESO(ICH2,NTRACKF) = XMR(IST,JJK(IST),2)
2183                   YMESO(ICH2,NTRACKF) = YMR(IST,JJK(IST),2)
2184                ENDIF
2185             ELSE
2186                ICH2 = ICH1+1
2187                JJO(ICH2,NTRACKF) = JJK(IST)
2188                XMESO(ICH2,NTRACKF) = XMR(IST,JJK(IST),2)
2189                YMESO(ICH2,NTRACKF) = YMR(IST,JJK(IST),2)
2190             ENDIF
2191          ENDDO
2192
2193          CALL CHFILL(700,SNGL(XVERT),R1,R2)
2194          CALL CHFILL(701,SNGL(YVERT),R1,R2)
2195
2196          IF (ITCHECK(IT).EQ.1) THEN
2197
2198
2199             IF (ID(IP(1,JJ1)).EQ.IDFIND .AND.
2200      &           ID(IP(2,JJ2)).EQ.IDFIND .AND.
2201      &           ID(IP(3,JJ3)).EQ.IDFIND .AND.
2202      &           ID(IP(4,JJ4)).EQ.IDFIND) THEN
2203
2204                NMUONF = NMUONF+1 ! Bon muon
2205
2206                PT = SQRT(PVERT1(IP(5,JJ5))**2+PVERT2(IP(5,JJ5))**2)
2207                PZ =  PVERT3(IP(5,JJ5))
2208                E = SQRT(PT**2+PZ**2+0.1056**2)
2209                YRAP = 0.5*DLOG((E+PZ)/(E-PZ))
2210                PHIM = DATAN2(PVERT2(IP(5,JJ5)),PVERT1(IP(5,JJ5)))
2211                CALL CHFILL(801,SNGL(YRAP),R1,R2)
2212                CALL CHFILL(901,SNGL(PT),R1,R2)
2213                CALL CHFILL(911,SNGL(PHIM),R1,R2)
2214             ELSE
2215
2216                NGHOST = NGHOST+1 ! ghost
2217
2218                PT = SQRT(PVERT1(IP(5,JJ5))**2+PVERT2(IP(5,JJ5))**2)
2219                PZ =  PVERT3(IP(5,JJ5))
2220                E = SQRT(PT**2+PZ**2+0.1056**2)
2221                YRAP = 0.5*DLOG((E+PZ)/(E-PZ))
2222                PHIM = DATAN2(PVERT2(IP(5,JJ5)),PVERT1(IP(5,JJ5)))
2223                CALL CHFILL(802,SNGL(YRAP),R1,R2)
2224                CALL CHFILL(902,SNGL(PT),R1,R2)
2225                CALL CHFILL(912,SNGL(PHIM),R1,R2)
2226             ENDIF
2227          ENDIF
2228       ENDIF
2229       IF (ITCHECK(ITRACK(IP(5,JJ5))).EQ.1) THEN
2230          IF (JJ3.EQ.0) NERR(3) = NERR(3)+1
2231          IF (JJ3.NE.0) THEN
2232             IF (JJ1.EQ.0) NERR(1) = NERR(1)+1
2233             IF (JJ2.EQ.0) NERR(2) = NERR(2)+1
2234          ENDIF
2235          PT = SQRT(PVERT1(IP(5,JJ5))**2+PVERT2(IP(5,JJ5))**2)
2236          PZ =  PVERT3(IP(5,JJ5))
2237          E = SQRT(PT**2+PZ**2+0.1056**2)
2238          YRAP = 0.5*DLOG((E+PZ)/(E-PZ))
2239          PHIM = DATAN2(PVERT2(IP(5,JJ5)),PVERT1(IP(5,JJ5)))
2240          CALL CHFILL(800,SNGL(YRAP),R1,R2)
2241          CALL CHFILL(900,SNGL(PT),R1,R2)
2242          CALL CHFILL(910,SNGL(PHIM),R1,R2)
2243       ENDIF
2244
2245       RETURN
2246       END
2247
2248
2249 *************************************************************************
2250       SUBROUTINE RECO_SUM
2251 *************************************************************************
2252 *
2253 *
2254 *************************************************************************
2255       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2256 *
2257       PARAMETER (MAXHITCH=10000,MAXTRK=50000,MAXHITTOT=20000,
2258      &           NBSTATION=5)
2259       PARAMETER (NBCHAMBER=10,NTRMAX=500)
2260 *
2261       COMMON/TRACKFOUT/IEVOUT,NTREVT,JJOUT(NBCHAMBER,NTRMAX),
2262      &               ISTAT(NTRMAX),PXZOUT(NTRMAX),TPHIOUT(NTRMAX),
2263      &               TALAMOUT(NTRMAX),XVERTOUT(NTRMAX),YVERTOUT(NTRMAX),
2264      &               CHI2OUT(NTRMAX),
2265      &               XMESOUT(NBCHAMBER,NTRMAX),YMESOUT(NBCHAMBER,NTRMAX)
2266      &              ,PXVOUT(NTRMAX),PYVOUT(NTRMAX),PZVOUT(NTRMAX)
2267 *
2268       COMMON/CHHIT/XM(NBSTATION,MAXHITCH),YM(NBSTATION,MAXHITCH),
2269      &             PHM(NBSTATION,MAXHITCH),ALM(NBSTATION,MAXHITCH),
2270      &             IZM(NBSTATION,MAXHITCH),
2271      &             IP(NBSTATION,MAXHITCH),JHIT(NBSTATION),
2272      &             XMR(NBSTATION,MAXHITCH,2),YMR(NBSTATION,MAXHITCH,2)
2273 *
2274       COMMON/RHIT/ITYP(MAXHITTOT),XTR(MAXHITTOT),YTR(MAXHITTOT),
2275      &            PTOT(MAXHITTOT),ID(MAXHITTOT),IZST(MAXHITTOT),
2276      &            PVERT1(MAXHITTOT),PVERT2(MAXHITTOT),PVERT3(MAXHITTOT),
2277      &            ZVERT(MAXHITTOT),NHITTOT
2278 *
2279       COMMON/RECOUT/JJO(NBCHAMBER,NTRMAX),XMESO(NBCHAMBER,NTRMAX),
2280      &              YMESO(NBCHAMBER,NTRMAX)
2281 *
2282       COMMON/TRACKSUM/NRES(5),NRESF,NTRMUALL,NMUONALL,NGHOSTALL,
2283      &     NTRACKFALL,NERRALL(NBSTATION),IR
2284 *
2285       COMMON/PRECSUM/NRESF1,NMUONALL1,NGHOSTALL1,NTRACKFALL1
2286 *
2287       REAL*4 PXR,PYR,PZR,ZVR,CHI2R,PXV,PYV,PZV
2288       COMMON/PAWCR4/IEVR,NTRACKR,ISTATR(NTRMAX),ISIGNR(NTRMAX),
2289      &              PXR(NTRMAX),PYR(NTRMAX),PZR(NTRMAX),ZVR(NTRMAX),
2290      &              CHI2R(NTRMAX),PXV(NTRMAX),PYV(NTRMAX),PZV(NTRMAX)
2291 *
2292       COMMON/DEBEVT/IDEBUG
2293 *
2294       NMUF = 0
2295       NGHF = 0
2296       DO ITR = 1,NTREVT
2297          NTRACKFALL1 =  NTRACKFALL1 + 1
2298          IF (ISTAT(ITR).EQ.1) THEN
2299             NMUF = NMUF + 1
2300             NMUONALL1 = NMUONALL1 + 1
2301          ELSEIF (ISTAT(ITR).EQ.2) THEN
2302             NGHF = NGHF + 1
2303             NGHOSTALL1 = NGHOSTALL1 + 1
2304          ENDIF
2305       ENDDO
2306       IF (NMUF.EQ.2) NRESF1 = NRESF1 + 1
2307 *
2308       NTRACKR = NTREVT
2309       DO ITR = 1,NTREVT
2310          ISTATR(ITR) = ISTAT(ITR)
2311          ISIGNR(ITR) = 1
2312          IF (PXZOUT(ITR).LT.0.) ISIGNR(ITR) = -1
2313          PXZ = ABS(PXZOUT(ITR))
2314          PHI = ATAN(TPHIOUT(ITR))
2315          ALAM = ATAN(TALAMOUT(ITR))
2316          PYR(ITR) = PXZ*SIN(PHI)
2317          PXR(ITR) = PXZ*TAN(ALAM)
2318          PZR(ITR) = PXZ*COS(PHI)
2319          ZVR(ITR) = 0.
2320          CHI2R(ITR) =  CHI2OUT(ITR)
2321          PXV(ITR) = PXVOUT(ITR)
2322          PYV(ITR) = PYVOUT(ITR)
2323          PZV(ITR) = PZVOUT(ITR)
2324       ENDDO
2325
2326       CALL CHFNT(IEVR,NTRACKR,ISTATR,ISIGNR,
2327      &     PXR,PYR,PZR,ZVR,CHI2R,PXV,PYV,PZV)
2328
2329       IF (IDEBUG.GE.2) THEN
2330          PRINT *,'RECO_SUM evt number :',IEVOUT
2331          PRINT *,'RECO_SUM nb of track /evt     :',NTREVT
2332          PRINT *,'RECO_SUM nb of good muon /evt :',NMUF
2333          PRINT *,'RECO_SUM nb of ghost /evt     :',NGHF
2334          IF (NTREVT.GT.0) THEN
2335             DO ITR = 1,NTREVT
2336                PRINT *,'RECO_SUM track number :',ITR
2337                PRINT *,'RECO_SUM PXZOUT :',PXZOUT(ITR)
2338                PRINT *,'RECO_SUM CHI2OUT :',CHI2OUT(ITR)
2339                PRINT *,' PX GEANT= ', PXV(ITR),' PX RECONS= ',PYR(ITR)
2340                PRINT *,' PY GEANT= ', PYV(ITR),' PY RECONS= ',PXR(ITR)
2341                PRINT *,' PZ GEANT= ', PZV(ITR),' PZ RECONS= ',PZR(ITR)
2342                PXZV = SQRT( PYV(ITR)**2+PZV(ITR)**2)
2343                PXZR = SQRT( PXR(ITR)**2+PZR(ITR)**2)
2344                PRINT *,' PXZ GEANT= ', PXZV,' PXZ RECONS= ',PXZR
2345                FIV=ATAN2(PYV(ITR),PZV(ITR))
2346                ALAMV=ATAN2(PXV(ITR),PXZV)
2347                FIR=ATAN2(PXR(ITR),PZR(ITR))
2348                ALAMR=ATAN2(PYR(ITR),PXZR)
2349 **               PRINT *,' PHI GEANT= ',FIV,' PXZ RECONS= ',FIR
2350 **               PRINT *,' ALAM GEANT= ',ALAMV,' ALAM RECONS= ',ALAMR
2351           ENDDO
2352          ENDIF
2353       ENDIF
2354
2355       RETURN
2356       END
2357
2358 *************************************************************************
2359       SUBROUTINE RECO_TERM
2360 *************************************************************************
2361 *
2362 *
2363 *************************************************************************
2364       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2365 *
2366       PARAMETER(NBSTATION=5)
2367 *
2368       COMMON/TRACKSUM/NRES(5),NRESF,NTRMUALL,NMUONALL,NGHOSTALL,
2369      &     NTRACKFALL,NERRALL(NBSTATION),IR
2370 *
2371       COMMON/PRECSUM/NRESF1,NMUONALL1,NGHOSTALL1,NTRACKFALL1
2372 *
2373       COMMON/DEBEVT/IDEBUG
2374 *
2375 *      CHARACTER*50 FILEBKG,FILERES,FILEOUT,FILEMIN
2376 *
2377       IF (IDEBUG.GE.1) THEN
2378          PRINT *,'    '
2379          PRINT *,'RECO_TERM ***** SUMMARY TRACK-FINDING *****'
2380          PRINT *,'RECO_TERM nb of resonances :',NRES(5)
2381          PRINT *,'RECO_TERM nb of resonances  45 :',NRES(1)
2382          PRINT *,'RECO_TERM nb of resonances 345 :',NRES(2)
2383          PRINT *,'RECO_TERM nb of resonances found :',NRESF
2384          PRINT *,'RECO_TERM nb of muon track       :',NTRMUALL
2385          PRINT *,'RECO_TERM nb of track found      :',NTRACKFALL
2386          PRINT *,'RECO_TERM nb of muon track found :',NMUONALL
2387          PRINT *,'RECO_TERM nb of ghost found      :',NGHOSTALL
2388          DO I = 1,4
2389             PRINT *,'RECO_TERM nb of error in st',I,':',NERRALL(I)
2390          ENDDO
2391
2392          PRINT *,'    '
2393          PRINT *,'RECO_TERM ***** SUMMARY PRECISION *****'
2394          PRINT *,'RECO_TERM nb of resonances found :',NRESF1
2395          PRINT *,'RECO_TERM nb of track found      :',NTRACKFALL1
2396          PRINT *,'RECO_TERM nb of muon track found :',NMUONALL1
2397          PRINT *,'RECO_TERM nb of ghost found      :',NGHOSTALL1
2398       ENDIF
2399 *
2400       CALL HIST_CLOSED
2401
2402       RETURN
2403       END
2404
2405 *************************************************************************
2406       SUBROUTINE RECO_PRECISION
2407 *************************************************************************
2408 *
2409 *
2410 *************************************************************************
2411       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2412 *
2413       PARAMETER (MAXHITCH=10000,MAXTRK=50000,MAXHITTOT=20000,
2414      &           NBSTATION=5,MAXCAN=1000,NTRMAX=500)
2415       PARAMETER (NPLANE=10)
2416 *
2417       COMMON/PARAM/ZPLANEP(NPLANE),THICK,XPREC,YPREC,B0,BL3,ZMAGS,
2418      &             ZMAGE,ZABS,XMAG,ZBP1,ZBP2,CONST
2419 *
2420       COMMON/CHHIT/XM(NBSTATION,MAXHITCH),YM(NBSTATION,MAXHITCH),
2421      &             PHM(NBSTATION,MAXHITCH),ALM(NBSTATION,MAXHITCH),
2422      &             IZM(NBSTATION,MAXHITCH),
2423      &             IP(NBSTATION,MAXHITCH),JHIT(NBSTATION),
2424      &             XMR(NBSTATION,MAXHITCH,2),YMR(NBSTATION,MAXHITCH,2)
2425 *
2426       COMMON/RHIT/ITYP(MAXHITTOT),XTR(MAXHITTOT),YTR(MAXHITTOT),
2427      &            PTOT(MAXHITTOT),ID(MAXHITTOT),IZST(MAXHITTOT),
2428      &            PVERT1(MAXHITTOT),PVERT2(MAXHITTOT),PVERT3(MAXHITTOT),
2429      &            ZVERT(MAXHITTOT),NHITTOT
2430 *
2431       COMMON/TRACKFOUT/IEVOUT,NTREVT,JJOUT(NPLANE,NTRMAX),
2432      &               ISTAT(NTRMAX),PXZOUT(NTRMAX),TPHIOUT(NTRMAX),
2433      &               TALAMOUT(NTRMAX),XVERTOUT(NTRMAX),YVERTOUT(NTRMAX),
2434      &               CHI2OUT(NTRMAX),
2435      &               XMESOUT(NPLANE,NTRMAX),YMESOUT(NPLANE,NTRMAX)
2436      &              ,PXVOUT(NTRMAX),PYVOUT(NTRMAX),PZVOUT(NTRMAX)
2437
2438 *
2439       COMMON/MEAS/LPLANE(NPLANE),XMP(NPLANE),YMP(NPLANE)
2440 *
2441       COMMON/FCNOUT/PXZEA,ALAMEA,PHIEA,XEA,YEA,NPLU,CHI2
2442 *
2443       COMMON/PRECCUT/PCUT,PTCUT,CHI2CUT
2444 *
2445       DIMENSION PARMU(MAXCAN,NPLANE,2),LPLANEMU(MAXCAN,NPLANE)
2446 *
2447       IF (NTREVT.EQ.0) RETURN
2448
2449       DO ITR = 1,NTREVT
2450          ICH = 0
2451          DO IST = 1,NBSTATION
2452             DO ILOOP = 1,2
2453                ICH = ICH + 1
2454 **               print *,' ich=',ich
2455                IF (JJOUT(ICH,ITR).GT.0) THEN
2456                   LPLANEMU(ITR,ICH) = 1
2457                   PARMU(ITR,ICH,1) = XMESOUT(ICH,ITR)
2458                   PARMU(ITR,ICH,2) = YMESOUT(ICH,ITR)
2459 **                  print *,' x,y ', PARMU(ITR,ICH,1),PARMU(ITR,ICH,2)
2460                ELSE
2461                   LPLANEMU(ITR,ICH) = 0
2462                ENDIF
2463             ENDDO
2464          ENDDO
2465       ENDDO
2466 *
2467       NTRACK = 0
2468       DO ICAN = 1,NTREVT
2469          DO ICH = 1,NPLANE
2470             LPLANE(ICH) = LPLANEMU(ICAN,ICH)
2471             XMP(ICH) = PARMU(ICAN,ICH,1)
2472             YMP(ICH) = PARMU(ICAN,ICH,2)
2473          ENDDO
2474
2475          IF (LPLANE(1).GT.0) THEN
2476             X1 = XMP(1)
2477             Y1 = YMP(1)
2478             IPL1 = 1
2479          ELSE
2480             X1 = XMP(2)
2481             Y1 = YMP(2)
2482             IPL1 = 2
2483          ENDIF
2484          IF (LPLANE(3).GT.0) THEN
2485             X2 = XMP(3)
2486             Y2 = YMP(3)
2487             IPL2 = 3
2488          ELSE
2489             X2 = XMP(4)
2490             Y2 = YMP(4)
2491             IPL2 = 4
2492          ENDIF
2493          IF (LPLANE(7).GT.0) THEN
2494             X3 = XMP(7)
2495             IPL3 = 7
2496          ELSE
2497             X3 = XMP(8)
2498             IPL3 = 8
2499          ENDIF
2500          IF (LPLANE(9).GT.0) THEN
2501             X4 = XMP(9)
2502             IPL4 = 9
2503          ELSE
2504             X4 = XMP(10)
2505             IPL4 = 10
2506          ENDIF
2507
2508          PHIAV = DATAN2((X2-X1),(ZPLANEP(IPL2)-ZPLANEP(IPL1)))
2509          PHIAP = DATAN2((X4-X3),(ZPLANEP(IPL4)-ZPLANEP(IPL3)))
2510          DPHI = (PHIAP-PHIAV)
2511          ASIGN = 1.
2512          IF (DPHI.LT.0.) ASIGN = -1. ! CCC
2513          PXZ = CONST/DABS(DPHI)
2514 ** Cuts PXZ
2515          IF (PXZ.LT.PCUT) GO TO 2
2516
2517          PXZINVI = ASIGN/PXZ ! CCC
2518 **         PXZINVI = 1./PXZOUT(ICAN) ! CCC
2519 **         PXZINVI = -1./49. ! CCC
2520          PHII = PHIAV
2521          ALAMI = DATAN2((Y2-Y1),DSQRT((X2-X1)**2
2522      &            +(ZPLANEP(IPL2)-ZPLANEP(IPL1))**2))
2523          XVR = X1
2524          YVR = Y1
2525 **         print *,' avant prec_fit pxzi phii alami x y',1./ PXZINVI,
2526 **     &         PHII, ALAMI ,XVR,YVR
2527 **         PRINT *,' X1 X2 X3 X4',X1,X2,X3,X4
2528 **         PRINT *,' Z1 Z2 Z3 Z4',ZPLANEP(IPL1),ZPLANEP(IPL2),
2529 **     &          ZPLANEP(IPL3),ZPLANEP(IPL4)
2530 **         PRINT *,' CONST= ',CONST
2531
2532          PRINT *,' RECO_PRECISION CHI2OUT avant fit ',CHI2OUT(ICAN)
2533
2534          IF (CHI2OUT(ICAN).GT.CHI2CUT) GO TO 2
2535
2536 ** Fit des traces apres l'absorbeur
2537          CALL PREC_FIT (PXZINVI,PHII,ALAMI,XVR,YVR,
2538      &          PXZINVF,PHIF,ALAMF,XVERTF,YVERTF,EPXZINV,EPHI,EALAM,
2539      &          EXVERT,EYVERT)
2540
2541 ** Correction de Branson
2542          CALL BRANSON(PXZEA,PHIEA,ALAMEA,XEA,YEA)
2543
2544          PXZ1 = DABS(PXZEA)
2545          PX1 = PXZ1*DSIN(PHIEA)
2546          PY1 = PXZ1*DTAN(ALAMEA)
2547          PT1 = DSQRT(PX1**2+PY1**2)
2548 ** Cuts PT
2549          IF (PT1.LT.PTCUT) GO TO 2
2550 ** Cuts CHI2
2551          IF ((CHI2/FLOAT(2*NPLU-5)).GT.CHI2CUT) GO TO 2
2552
2553          NTRACK = NTRACK + 1
2554          DO ICH = 1,NPLANE
2555             JJOUT(ICH,NTRACK) =  JJOUT(ICH,ICAN)
2556             XMESOUT(ICH,NTRACK) =  XMESOUT(ICH,ICAN)
2557             YMESOUT(ICH,NTRACK) =  YMESOUT(ICH,ICAN)
2558          ENDDO
2559          ISTAT(NTRACK) = ISTAT(ICAN)
2560          PXZOUT(NTRACK) = PXZEA
2561          TPHIOUT(NTRACK) = DTAN(PHIEA)
2562          TALAMOUT(NTRACK) = DTAN(ALAMEA)
2563          XVERTOUT(NTRACK) = XEA
2564          YVERTOUT(NTRACK) = YEA
2565          CHI2OUT(NTRACK) =  CHI2/FLOAT(2*NPLU-5)
2566
2567          PRINT *,' RECO_PRECISION CHI2OUT apres fit',CHI2OUT(NTRACK)
2568
2569 **         print *,' reco_precision pxz tphi talam xvert yvert chi2',
2570 **     &       PXZEA,PHIEA,ALAMEA,
2571 **     &       XEA,YEA,CHI2/FLOAT(2*NPLU-5)
2572   2      CONTINUE
2573       ENDDO
2574       NTREVT = NTRACK
2575
2576       RETURN
2577       END
2578
2579
2580 ************************************************************************
2581       SUBROUTINE ORDONNE_HIT(ICH,HCUT)
2582 **************************************************************************
2583 *
2584 *  Sort hits in station ICH according to the "impact parameter" HHIT
2585 *
2586 **************************************************************************
2587
2588       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2589
2590       PARAMETER (MAXHITCH=10000,MAXTRK=50000,MAXHITTOT=20000,
2591      &           NBSTATION=5)
2592
2593       COMMON/CHHIT/XM(NBSTATION,MAXHITCH),YM(NBSTATION,MAXHITCH),
2594      &             PHM(NBSTATION,MAXHITCH),ALM(NBSTATION,MAXHITCH),
2595      &             IZM(NBSTATION,MAXHITCH),
2596      &             IP(NBSTATION,MAXHITCH),JHIT(NBSTATION),
2597      &             XMR(NBSTATION,MAXHITCH,2),YMR(NBSTATION,MAXHITCH,2)
2598
2599       COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
2600
2601       COMMON/VERIFGEANT/ITTROUGH(MAXTRK,NBSTATION),
2602      &                  IT_LIST(MAXTRK),IT_NP(MAXTRK),ITCHECK(MAXTRK),
2603      &                  ITRACK(MAXHITTOT)
2604
2605       COMMON/HCHHIT/HHIT(MAXHITCH),INDEXTAB(MAXHITCH),INDEXMAX
2606 *
2607       COMMON/DEBEVT/IDEBUG
2608 *
2609       REAL*4 H4(MAXHITCH)
2610 * tri des impulsion par ordre decroissant
2611 * le tab INDEXTAB contient les j ordonnes
2612 * INDEXMAX est l indice max du tableau = NBHIT si pas de contrainte
2613
2614       JJ=0
2615       INDEXMAX = 0
2616 * boucle sur le nombre de hits candidats de la station
2617
2618       DO J=1,JHIT(ICH)
2619
2620           IF (PHM(ICH,J).LT.6.3) THEN        !2pi=6.3 radian
2621             JJ=JJ+1
2622 * calcul du h dans XOY a z=0
2623             HHIT(J)=ABS(XM(ICH,J)-ZPLANE(ICH)*PHM(ICH,J))
2624 * cut en Pxz
2625             IF (HHIT(J).LT.HCUT) THEN
2626                INDEXMAX=INDEXMAX+1
2627                INDEXTAB(INDEXMAX)=J
2628             ELSEIF(ITCHECK(ITRACK(IP(ICH,J))).EQ.1) THEN
2629                IF (IDEBUG.GE.2) THEN
2630                   PRINT *,'ORDONNE_HIT rejet muon/res in st.',ICH,
2631      &                 ' h=',HHIT(J)
2632                ENDIF
2633             ENDIF
2634          ENDIF
2635       ENDDO
2636
2637       NBHIT=JHIT(ICH)
2638
2639       DO I = 1,NBHIT
2640          H4(I) = SNGL(HHIT(I))
2641       ENDDO
2642
2643       CALL SORTZV(H4,INDEXTAB,INDEXMAX,-1,0,1)
2644
2645       DO I = 1,NBHIT
2646          HHIT(I) = DBLE(H4(I))
2647       ENDDO
2648 **      PRINT *,'ORDONNE st. numero',ICH
2649 **      PRINT *,'ORDONNE nb de hits initiaux dans st.',ICH,':',JHIT(ICH)
2650 **      PRINT *,'ORDONNE nb de hits avec mes. angulaire:',JJ
2651 **      PRINT *,'ORDONNE nb de hits avec mes. ang. et cut en Pxz:',INDEXMAX
2652       IF (IDEBUG.GE.2) THEN
2653          PRINT *,'ORDONNE_HIT nb de hits accepte dans st.',ICH,':',
2654      &           INDEXMAX
2655       ENDIF
2656
2657       RETURN
2658       END
2659 ***********************************************************************************
2660       SUBROUTINE DISTMIN4(X1,Y1,PHI1,ALAM1,ICH,EX,EY,EPHI,ELAM,IFIND)
2661 ***********************************************************************************
2662 *    Find the nearest hit in station ICH in the (X,Y,lambda,phi) phase space
2663 *
2664 ***********************************************************************************
2665
2666       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2667
2668       PARAMETER (MAXHITCH=10000,NBSTATION=5)
2669
2670       COMMON/CHHIT/XM(NBSTATION,MAXHITCH),YM(NBSTATION,MAXHITCH),
2671      &             PHM(NBSTATION,MAXHITCH),ALM(NBSTATION,MAXHITCH),
2672      &             IZM(NBSTATION,MAXHITCH),
2673      &             IP(NBSTATION,MAXHITCH),JHIT(NBSTATION),
2674      &             XMR(NBSTATION,MAXHITCH,2),YMR(NBSTATION,MAXHITCH,2)
2675
2676       IFIND = 0
2677       DISTMIN=4.
2678       DO I=1,JHIT(ICH)
2679          IF (PHM(ICH,I).LE.6.3) THEN ! angles measured
2680             IF (ABS(PHI1-PHM(ICH,I)) .LT. EPHI .AND.
2681      &         ABS(ALAM1-ALM(ICH,I)) .LT. ELAM .AND.
2682      &         ABS(X1-XM(ICH,I)) .LT. EX .AND.
2683      &         ABS(Y1-YM(ICH,I)) .LT. EY) THEN
2684                DIST = ((PHI1-PHM(ICH,I))/EPHI)**2 +
2685      &                ((ALAM1-ALM(ICH,I))/ELAM)**2 +
2686      &                ((X1-XM(ICH,I))/EX)**2 +
2687      &                ((Y1-YM(ICH,I))/EY)**2
2688                IF (DIST .LT. DISTMIN) THEN
2689                   DISTMIN = DIST
2690                   IFIND = I
2691                ENDIF
2692             ENDIF
2693          ENDIF
2694       ENDDO
2695
2696       RETURN
2697       END
2698 ***********************************************************************************
2699       SUBROUTINE DISTMIN2(X1,Y1,X2,Y2,ICH,EX1,EY1,EX2,EY2,IFIND)
2700 ***********************************************************************************
2701 *    Find the nearest hit in station ICH in the (X,Y) space
2702 *
2703 ***********************************************************************************
2704       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2705
2706       PARAMETER (MAXHITCH=10000,NBSTATION=5)
2707
2708       COMMON/CHHIT/XM(NBSTATION,MAXHITCH),YM(NBSTATION,MAXHITCH),
2709      &             PHM(NBSTATION,MAXHITCH),ALM(NBSTATION,MAXHITCH),
2710      &             IZM(NBSTATION,MAXHITCH),
2711      &             IP(NBSTATION,MAXHITCH),JHIT(NBSTATION),
2712      &             XMR(NBSTATION,MAXHITCH,2),YMR(NBSTATION,MAXHITCH,2)
2713
2714       IFIND = 0
2715       DISTMIN=2.
2716       DO I=1,JHIT(ICH)
2717          IF (IZM(ICH,I).EQ.1) THEN
2718             X = X1
2719             Y = Y1
2720          ELSE
2721             X = X2
2722             Y = Y2
2723          ENDIF
2724          EX = EX1
2725          EY = EY1
2726          IF (ICH.EQ.4.OR.ICH.EQ.5) THEN
2727             IF (IZM(ICH,I).EQ.1) THEN
2728                EX = EX1
2729                EY = EY1
2730             ELSE
2731                EX = EX2
2732                EY = EY2
2733             ENDIF
2734          ENDIF
2735          IF (ABS(X-XM(ICH,I)) .LT. EX .AND.
2736      &      ABS(Y-YM(ICH,I)) .LT. EY) THEN
2737             DIST = ((X-XM(ICH,I))/EX)**2 +
2738      &             ((Y-YM(ICH,I))/EY)**2
2739             IF (DIST .LT. DISTMIN) THEN
2740                DISTMIN = DIST
2741                IFIND = I
2742             ENDIF
2743          ENDIF
2744       ENDDO
2745
2746       RETURN
2747       END
2748 ********************************************************************************
2749       SUBROUTINE H_ACCEPTANCE(ICH)
2750 ********************************************************************************
2751 * Etude de l'acceptance des resonnances en  fonction du H
2752 * dans la station ICH
2753 *
2754 *  INPUT :    ICH
2755 *  OUTPUT :  Histo #1
2756 ********************************************************************************
2757
2758       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2759
2760       PARAMETER (MAXHITCH=10000,MAXTRK=50000,MAXHITTOT=20000,
2761      &           NBSTATION=5)
2762
2763       COMMON/CHHIT/XM(NBSTATION,MAXHITCH),YM(NBSTATION,MAXHITCH),
2764      &             PHM(NBSTATION,MAXHITCH),ALM(NBSTATION,MAXHITCH),
2765      &             IZM(NBSTATION,MAXHITCH),
2766      &             IP(NBSTATION,MAXHITCH),JHIT(NBSTATION),
2767      &             XMR(NBSTATION,MAXHITCH,2),YMR(NBSTATION,MAXHITCH,2)
2768
2769       COMMON/RHIT/ITYP(MAXHITTOT),XTR(MAXHITTOT),YTR(MAXHITTOT),
2770      &            PTOT(MAXHITTOT),ID(MAXHITTOT),IZST(MAXHITTOT),
2771      &            PVERT1(MAXHITTOT),PVERT2(MAXHITTOT),PVERT3(MAXHITTOT),
2772      &            ZVERT(MAXHITTOT),NHITTOT
2773
2774
2775       COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
2776
2777       COMMON/VERIFGEANT/ITTROUGH(MAXTRK,NBSTATION),
2778      &                  IT_LIST(MAXTRK),IT_NP(MAXTRK),ITCHECK(MAXTRK),
2779      &                  ITRACK(MAXHITTOT)
2780
2781       COMMON/HCHHIT/HHIT(MAXHITCH),INDEXTAB(MAXHITCH),INDEXMAX
2782
2783       REAL*4 R1,R2
2784       DATA R1,R2/0.,1./
2785
2786       NMUONI = 0
2787       DO J = 1,JHIT(ICH)
2788          IF (ITYP(IP(ICH,J)).EQ.5.OR.ITYP(IP(ICH,J)).EQ.6) THEN
2789             ISTAK = ID(IP(ICH,J))
2790             ISTAK = MOD(ISTAK,30000)
2791             ISTAK = MOD(ISTAK,10000)
2792             IF (ISTAK.EQ.0) THEN
2793 **              PRINT *,'ACCEPT. id du muon dans st.',ICH,':',ITYP(IP(ICH,J))
2794                NMUONI = NMUONI+1
2795             ENDIF
2796          ENDIF
2797       ENDDO
2798 *      PRINT *,'ACCEPT. nb de muons/res total dans st.',ICH,':',NMUONI
2799 *      pause
2800
2801       DO IH = 1,500
2802          HCUT = IH
2803 *   Sort hits in st. z
2804          CALL ORDONNE_HIT(ICH,HCUT)
2805          NMUON = 0
2806          DO IND = 1,INDEXMAX
2807             IIND = IP(ICH,INDEXTAB(IND))
2808             IDPART = ITYP(IIND)
2809             ISTAK = ID(IIND)
2810             ISTAK = MOD(ISTAK,30000)
2811             ISTAK = MOD(ISTAK,10000)
2812 **            PRINT *,' IDPART=',IDPART,' ISTAK=',ISTAK
2813             IF (IDPART.EQ.5.OR.IDPART.EQ.6.AND.ISTAK.EQ.0) THEN
2814                NMUON = NMUON+1
2815             ENDIF
2816          ENDDO
2817          IF (NMUON.EQ.2.AND.NMUONI.EQ.2) THEN
2818             CALL CHFILL(ICH*100,SNGL(HCUT),R1,R2)
2819          ENDIF
2820       ENDDO
2821
2822       RETURN
2823       END
2824
2825 ********************************************************************************
2826        SUBROUTINE FOLLOW(ZSTR,PEST,IFLAG,XPL,YPL,PHPL,ALPL)
2827 ********************************************************************************
2828 *   Calculate the particle trajectory in the spectrometer and
2829 *   (XPL,YPL,PHPL,ALPL)
2830 *   for the 5 stations.
2831 *
2832 ********************************************************************************
2833        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2834 *
2835        PARAMETER(NBSTATION=5)
2836 *
2837        DIMENSION XPL(NBSTATION,2),YPL(NBSTATION,2),PHPL(NBSTATION),
2838      &           ALPL(NBSTATION),PEST(NBSTATION)
2839
2840        COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
2841
2842        COMMON /MEASUR/XMES(NBSTATION),YMES(NBSTATION),IZMES(NBSTATION),
2843      &                PHMES(NBSTATION),ALMES(NBSTATION),MPOS(NBSTATION),
2844      &                MANG(NBSTATION)
2845        COMMON /MAGNET/BL3,B0
2846
2847        LOGICAL LFLAG, LFLAG1
2848
2849        XSTR = PEST(4)
2850        YSTR = PEST(5)
2851        PXZINV = PEST(1)
2852        TPHI = PEST(2)
2853        PHI = ATAN(TPHI)
2854        TALAM = PEST(3)
2855        PXZ = 1.0/PXZINV
2856        PY = ABS(PXZ)*TALAM
2857        PX = -ABS(PXZ)*SIN(PHI)
2858        PZ = -ABS(PXZ)*COS(PHI)
2859        PXY = SQRT(PX**2 + PY**2)
2860        FI=ATAN2(PY,PX)
2861        SINFI = SIN(FI)
2862        COSFI = COS(FI)
2863        TTHET = PZ/PXY
2864        RS = PXY*(100.0/(0.299792458*BL3))
2865        IF(PXZINV.LT.0.0) RS = -RS
2866 *       XC = XSTR + RS*SIN(FI)
2867 *       YC = YSTR - RS*COS(FI)
2868        PX0 = PX
2869        PY0 = PY
2870        LFLAG = .TRUE.
2871        LFLAG1 = .TRUE.
2872 *       PRINT *, XC,YC,RS,FI,TTHET,PXY,PZ
2873
2874        DO J = 1,5
2875
2876           IF (IFLAG.EQ.3 .OR. MPOS(J).EQ.1) THEN
2877              IF(ZPLANE(J) .GT. ZCOIL) THEN
2878 *                DFI =  (ZPLANE(J)-ZSTR)/(TTHET*RS)
2879 *                FIN = FI - DFI
2880 *                XPL(J,1) = XC - RS*SIN(FIN)
2881 *                YPL(J,1) = YC + RS*COS(FIN)
2882                  DFR =  (ZPLANE(J)-ZSTR)/TTHET
2883                  XPL(J,1) = XSTR + DFR*COSFI + 0.5D0*DFR*DFR*SINFI/RS
2884                  YPL(J,1) = YSTR + DFR*SINFI - 0.5D0*DFR*DFR*COSFI/RS
2885                  DFR2 =  (ZPLANE(J)-DZ_PL(J)-ZSTR)/TTHET
2886                  XPL(J,2) = XSTR + DFR2*COSFI + 0.5D0*DFR2*DFR2*SINFI/RS
2887                  YPL(J,2) = YSTR + DFR2*SINFI - 0.5D0*DFR2*DFR2*COSFI/RS
2888                 IF (IFLAG.EQ.3 .OR. MANG(J).EQ.1) THEN
2889 *                   PX=PXY*COS(FIN)
2890 *                   PY=PXY*SIN(FIN)
2891                    PX = PX0 + DFR * (PY0 - 0.5D0*PX0*DFR/RS) / RS
2892                    PY = PY0 - DFR * (PX0 + 0.5D0*PY0*DFR/RS) / RS
2893                    PHPL(J)=ATAN(PX/PZ)
2894                    ALPL(J)=ATAN(PY/SQRT(PX**2+PZ**2))
2895                 ENDIF
2896               ELSE
2897                 IF( LFLAG) THEN
2898 *                   DFI =  (ZCOIL-ZSTR)/(TTHET*RS)
2899 *                   FIN = FI - DFI
2900 *                   XCOIL = XC - RS*SIN(FIN)
2901 *                   YCOIL = YC + RS*COS(FIN)
2902                    DFR =  (ZCOIL-ZSTR)/TTHET
2903                    XCOIL = XSTR + DFR*COSFI + 0.5D0*DFR*DFR*SINFI/RS
2904                    YCOIL = YSTR + DFR*SINFI - 0.5D0*DFR*DFR*COSFI/RS
2905 *                   PX=PXY*COS(FIN)
2906 *                   PY=PXY*SIN(FIN)
2907                    PX = PX0 + DFR * (PY0 - 0.5D0*PX0*DFR/RS) / RS
2908                    PY = PY0 - DFR * (PX0 + 0.5D0*PY0*DFR/RS) / RS
2909                    PXZ = SQRT(PX**2 + PZ**2)
2910                    PHI=ATAN(PX/PZ)
2911                    TALAM = PY/PXZ
2912                    ALAM = ATAN(TALAM)
2913                    RD = PXZ*(100.0/(0.299792458*B0))
2914                    IF(PXZINV.LT.0.0) RD = -RD
2915                    ZC = ZCOIL - RD*SIN(PHI)
2916                    XC = XCOIL + RD*COS(PHI)
2917                    IF(ABS(ZMAGEND-ZC).GT.ABS(RD)) STOP 'FOLLOW'
2918                    LFLAG = .FALSE.
2919                 ENDIF
2920                 IF(ZPLANE(J) .GT. ZMAGEND) THEN
2921                   FIN = ASIN((ZPLANE(J) - ZC)/RD)
2922                   XPL(J,1)= XC - RD*COS(FIN)
2923                   YPL(J,1)= YCOIL - RD*(FIN - PHI)*TALAM
2924                   FIN2 = ASIN((ZPLANE(J)-DZ_PL(J) - ZC)/RD)
2925                   XPL(J,2)= XC - RD*COS(FIN2)
2926                   YPL(J,2)= YCOIL - RD*(FIN2 - PHI)*TALAM
2927                   PHPL(J)=FIN
2928                   ALPL(J)=ALAM
2929                 ELSE
2930                   IF (LFLAG1) THEN
2931                     FIN = ASIN((ZMAGEND - ZC)/RD)
2932                     XMAGEND = XC - RD*COS(FIN)
2933                     YMAGEND = YCOIL - RD*(FIN - PHI)*TALAM
2934                     TPHI = TAN(FIN)
2935                     CPHI = COS(FIN)
2936                     LFLAG1 = .FALSE.
2937                   ENDIF
2938                   XPL(J,1) = XMAGEND + (ZPLANE(J)-ZMAGEND)*TPHI
2939                   YPL(J,1) = YMAGEND - (ZPLANE(J)-ZMAGEND)*TALAM/CPHI
2940                   XPL(J,2) = XMAGEND + (ZPLANE(J)-DZ_PL(J)-ZMAGEND)*TPHI
2941                   YPL(J,2) = YMAGEND - (ZPLANE(J)-DZ_PL(J)-ZMAGEND)*
2942      &                      TALAM/CPHI
2943                   PHPL(J)=FIN
2944                   ALPL(J)=ALAM
2945                 ENDIF
2946              ENDIF
2947           ENDIF
2948        ENDDO
2949        RETURN
2950        END
2951 ********************************************************************************
2952        SUBROUTINE NEWFOLLOW(ZSTR,PEST,IFLAG,XPL,YPL,PHPL,ALPL)
2953 ********************************************************************************
2954 *   Calculate the particle trajectory in the spectrometer
2955 *   (XPL,YPL,PHPL,ALPL)
2956 *   for the 5 stations.
2957 *       Runge Kutta
2958 ********************************************************************************
2959        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2960 *
2961        PARAMETER(NBSTATION=5,NPLANE=10)
2962 *
2963        DIMENSION XPL(NBSTATION,2),YPL(NBSTATION,2),PHPL(NBSTATION),
2964      &           ALPL(NBSTATION),PEST(NBSTATION)
2965
2966        COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
2967
2968        COMMON/PARAM/ZPLANEP(NPLANE),THICK,XPREC,YPREC,B0,BL3,ZMAGS,
2969      &             ZMAGE,ZABS,XMAG,ZBP1,ZBP2,CONST
2970
2971        COMMON /MEASUR/XMES(NBSTATION),YMES(NBSTATION),IZMES(NBSTATION),
2972      &                PHMES(NBSTATION),ALMES(NBSTATION),MPOS(NBSTATION),
2973      &                MANG(NBSTATION)
2974
2975
2976        DIMENSION VECT(7),VOUT(7)
2977
2978        STEP = 6.                ! 1 cm
2979        NSTEPMAX = 5000
2980
2981        ASIGN = 1.
2982        IF (PEST(1).LT.0.) ASIGN = -1.
2983        TPHI = -1.*PEST(2)
2984        PHI = DATAN(TPHI)
2985        TALAM = PEST(3)
2986        ALAM = DATAN(TALAM)
2987        PXZ = DABS(1./PEST(1))
2988
2989        PX = PXZ*DSIN(PHI)
2990        PY = PXZ*DTAN(ALAM)
2991        PZ = PXZ*DCOS(PHI)
2992        PTOT = PXZ/DCOS(ALAM)
2993
2994        VECT(1) = PEST(4)
2995        VECT(2) = PEST(5)
2996        VECT(3) = 0.
2997        VECT(4) = PX/PTOT
2998        VECT(5) = PY/PTOT
2999        VECT(6) = PZ/PTOT
3000        VECT(7) = PTOT
3001
3002        Z = VECT(3)
3003        NSTEP = 0
3004 *
3005 ** Runge Kutta
3006 **      PRINT *,' AV GRKUTA ASIGN',ASIGN,' THET',THET
3007        ISTOLD = 0
3008        DO ICH = 1,NPLANE
3009
3010           IST = INT(FLOAT(ICH+1)/2.)
3011
3012
3013           DO WHILE (Z.GE.0..AND.Z.LT.ABS(ZPLANEP(ICH))
3014      &         .AND.NSTEP.LE.NSTEPMAX)
3015              NSTEP = NSTEP+1
3016 **          WRITE(6,*) NSTEP,(VECT(I),I=1,7)
3017              CALL RECO_GRKUTA (ASIGN,STEP,VECT,VOUT)
3018              DO I = 1,7
3019                 VECT(I) = VOUT(I)
3020              ENDDO
3021              Z = VECT(3)
3022           ENDDO
3023           IF (IST.NE.ISTOLD) THEN
3024              IPCH = 1
3025           ELSE
3026              IPCH = 2
3027           ENDIF
3028           XPL(IST,IPCH) = VECT(1)-(Z-ABS(ZPLANEP(ICH)))*VECT(4)/VECT(6)
3029           YPL(IST,IPCH) = VECT(2)-(Z-ABS(ZPLANEP(ICH)))*VECT(5)/VECT(6)
3030           IF (IPCH.EQ.2) THEN
3031              DX = XPL(IST,2)-XPL(IST,1)
3032              DY = YPL(IST,2)-YPL(IST,1)
3033              PHPL(IST) = -1.*DATAN2(DX,DZ_PL(IST))
3034              ALPL(IST) = DATAN2(DY,DSQRT(DX**2+DZ_PL(IST)**2))
3035           ENDIF
3036           ISTOLD = IST
3037        ENDDO
3038 **          print *,' vect= ',vect(1),vect(2),vect(3),vect(4),vect(5),
3039 **     &        vect(6),vect(7)
3040
3041
3042        RETURN
3043        END
3044 *******************************************************************************
3045        SUBROUTINE FCN(NPAR,GRAD,FVAL,PEST,IFLAG,FUTIL)
3046 *******************************************************************************
3047 *   Calculate FVAL=CHI2 the function minimized by minuit for a given track
3048 *
3049 *******************************************************************************
3050        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
3051
3052        PARAMETER(NBSTATION=5)
3053
3054 *       DIMENSION PEST(*),GRAD(*)
3055        DIMENSION PEST(5),GRAD(5)
3056        DIMENSION PEEST(NBSTATION)
3057
3058        COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
3059
3060
3061        COMMON/PRECIS/EEXM(NBSTATION),EEYM(NBSTATION),EEPH(NBSTATION),
3062      &     EEAL(NBSTATION)
3063
3064        COMMON /MEASUR/XMES(NBSTATION),YMES(NBSTATION),IZMES(NBSTATION),
3065      &               PHMES(NBSTATION),ALMES(NBSTATION), MPOS(NBSTATION),
3066      &               MANG(NBSTATION)
3067
3068        COMMON /PLANE/XPL(NBSTATION,2),YPL(NBSTATION,2),PHPL(NBSTATION),
3069      &               ALPL(NBSTATION),CHI2PL
3070
3071        COMMON/VERTEX/ERRV,IVERTEX
3072
3073        PEEST(1)=PEST(1)
3074        PEEST(2)=PEST(2)
3075        PEEST(3)=PEST(3)
3076        IF(IVERTEX.EQ.1) THEN
3077           PEEST(4)=PEST(4)      ! position du vertex
3078           PEEST(5)=PEST(5)
3079        ELSE
3080           PEEST(4)=0.0D0
3081           PEEST(5)=0.0D0
3082        ENDIF
3083
3084        CALL FOLLOW (0.0D0,PEEST,IFLAG,XPL,YPL,PHPL,ALPL) ! calcul des
3085        IF(IFLAG.EQ.1) THEN                       ! points d impacts dans les
3086           PRINT *,'FCN ',XPL(4,1),XMES(4)        ! plans
3087           PRINT *,'FCN ',YPL(4,1),YMES(4)
3088           PRINT *,'FCN ',XPL(5,1),XMES(5)
3089           PRINT *,'FCN ',YPL(5,1),YMES(5)
3090        ENDIF
3091
3092 *       IF (IVERTEX.EQ.1) THEN
3093 *          FVAL = (PEST(4)/ERRV)**2 + (PEST(5)/ERRV)**2
3094 *       ELSE
3095           FVAL = 0.0D0
3096 *       ENDIF
3097        NPLPL = 0
3098        DO J = 1,NBSTATION
3099           IF (MPOS(J).EQ.1) THEN
3100              NPLPL = NPLPL+1
3101              XPLC = XPL(J,IZMES(J))
3102              YPLC = YPL(J,IZMES(J))
3103              FF = (XMES(J) - XPLC)/EEXM(J)
3104              FVAL =FVAL + FF**2
3105              FF = (YMES(J) - YPLC)/EEYM(J)
3106              FVAL =FVAL + FF**2
3107           ENDIF
3108           IF (MANG(J).EQ.1) THEN
3109              NPLPL = NPLPL+1
3110              FF = (PHMES(J) - PHPL(J))/EEPH(J)
3111              FVAL =FVAL + FF**2
3112              FF = (ALMES(J) - ALPL(J))/EEAL(J)
3113              FVAL =FVAL + FF**2
3114           ENDIF
3115        ENDDO
3116 **       PRINT *,'ST 1',XPL(1,1),XMES(1),YPL(1,1),YMES(1)
3117 **       PRINT *,'ST 2',XPL(2,1),XMES(2),YPL(2,1),YMES(2)
3118 **       PRINT *,'ST 3',XPL(3,1),XMES(3),YPL(3,1),YMES(3)
3119 **       PRINT *,'ST 4',XPL(4,1),XMES(4),YPL(4,1),YMES(4)
3120 **       PRINT *,'ST 5',XPL(5,1),XMES(5),YPL(5,1),YMES(5)
3121        NPARAM = 3
3122        IF (IVERTEX.EQ.1) NPARAM = 5
3123        CHI2PL = FVAL/FLOAT(2*NPLPL-NPARAM)
3124
3125        RETURN
3126        END
3127 ********************************************************************************
3128       SUBROUTINE STOCK_CANDIDAT(ICH1,JHITCH1,ICH2,IFIND,EXM,EYM,EPH,EAL,
3129      &                          NCAN,ICODE)
3130 ********************************************************************************
3131 *   Fill common CANDIDAT with track candidates from the search in stations 4&5
3132 *
3133 ********************************************************************************
3134
3135       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
3136
3137       PARAMETER (MAXHITTOT=20000,MAXHITCH=10000,NBSTATION=5,MAXCAN=1000)
3138
3139       COMMON/CHHIT/XM(NBSTATION,MAXHITCH),YM(NBSTATION,MAXHITCH),
3140      &             PHM(NBSTATION,MAXHITCH),ALM(NBSTATION,MAXHITCH),
3141      &             IZM(NBSTATION,MAXHITCH),
3142      &             IP(NBSTATION,MAXHITCH),JHIT(NBSTATION),
3143      &             XMR(NBSTATION,MAXHITCH,2),YMR(NBSTATION,MAXHITCH,2)
3144
3145       COMMON/RHIT/ITYP(MAXHITTOT),XTR(MAXHITTOT),YTR(MAXHITTOT),
3146      &            PTOT(MAXHITTOT),ID(MAXHITTOT),IZST(MAXHITTOT),
3147      &            PVERT1(MAXHITTOT),PVERT2(MAXHITTOT),PVERT3(MAXHITTOT),
3148      &            ZVERT(MAXHITTOT),NHITTOT
3149
3150
3151       COMMON/CANDIDAT/JCAN(NBSTATION,MAXCAN),JCANTYP(MAXCAN),
3152      &          EEX(MAXCAN),EEY(MAXCAN),EEP(MAXCAN),EEA(MAXCAN)
3153
3154 **      PRINT *,'STOCK st. init.=',ICH1,'id. init.=',ID(IP(ICH1,JHITCH1))
3155 **      PRINT *,'STOCK st. finale=',ICH2,'id. final=',ID(IP(ICH2,IFIND))
3156 **      PRINT *,'STOCK ifind=',IFIND
3157 **      PRINT *,'STOCK icode=',ICODE
3158       NCAN = NCAN+1
3159       JCAN(ICH1,NCAN) = JHITCH1
3160       JCAN(ICH2,NCAN) = IFIND
3161       JCANTYP(NCAN) = ICODE
3162       EEX(NCAN) = EXM
3163       EEY(NCAN) = EYM
3164       EEP(NCAN) = EPH
3165       EEA(NCAN) = EAL
3166
3167       RETURN
3168       END
3169 *****************************************************************************
3170       SUBROUTINE CHECK_HISTO4(IDHIST,ICH2,IHIT2,ICH1,IHIT1,
3171      &                        X1,Y1,PHI1,ALAM1,P1,EXM,EYM,EPH,EAL)
3172 *****************************************************************************
3173 *   Check hit IHIT2 with GEANT informations from hit HIT1
3174 *
3175 *    INPUT : ICH2 : No st. de recherche
3176 *            IDCH1,X1,Y1,PHI1,ALAM1,P1 : Trace de reference
3177 *    OUTPUT : JOK : No hit dans ICH2 appartenant a la meme trace.
3178 *
3179 *****************************************************************************
3180       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
3181
3182       PARAMETER (MAXHITCH=10000,MAXHITTOT=20000,NBSTATION=5,
3183      &           MAXTRK=50000)
3184
3185       COMMON/CHHIT/XM(NBSTATION,MAXHITCH),YM(NBSTATION,MAXHITCH),
3186      &             PHM(NBSTATION,MAXHITCH),ALM(NBSTATION,MAXHITCH),
3187      &             IZM(NBSTATION,MAXHITCH),
3188      &             IP(NBSTATION,MAXHITCH),JHIT(NBSTATION),
3189      &             XMR(NBSTATION,MAXHITCH,2),YMR(NBSTATION,MAXHITCH,2)
3190
3191       COMMON/RHIT/ITYP(MAXHITTOT),XTR(MAXHITTOT),YTR(MAXHITTOT),
3192      &            PTOT(MAXHITTOT),ID(MAXHITTOT),IZST(MAXHITTOT),
3193      &            PVERT1(MAXHITTOT),PVERT2(MAXHITTOT),PVERT3(MAXHITTOT),
3194      &            ZVERT(MAXHITTOT),NHITTOT
3195
3196       COMMON/VERIFGEANT/ITTROUGH(MAXTRK,NBSTATION),
3197      &                  IT_LIST(MAXTRK),IT_NP(MAXTRK),ITCHECK(MAXTRK),
3198      &                  ITRACK(MAXHITTOT)
3199 *
3200       COMMON/DEBEVT/IDEBUG
3201
3202       REAL*4 R2
3203       DATA R2/1./
3204
3205       JOK = 0
3206
3207       DO I=1,JHIT(ICH2)
3208          IF (PHM(ICH2,I).LE.6.3) THEN
3209             IF (ID(IP(ICH1,IHIT1)).EQ.ID(IP(ICH2,I))) THEN
3210                JOK = I
3211                IF (IDHIST.GT.0) THEN
3212                   CALL CHFILL2(IDHIST,SNGL(P1),
3213      &                        SNGL((X1-XM(ICH2,I))**2),1.)
3214                   CALL CHFILL2(IDHIST+1,SNGL(P1),
3215      &                        SNGL((Y1-YM(ICH2,I))**2),1.)
3216                   CALL CHFILL2(IDHIST+2,SNGL(P1),
3217      &                     SNGL((PHI1-PHM(ICH2,I))**2),1.)
3218                   CALL CHFILL2(IDHIST+3,SNGL(P1),
3219      &                     SNGL((ALAM1-ALM(ICH2,I))**2),1.)
3220                   CALL CHF1(400+IDHIST,
3221      &                        SNGL(X1-XM(ICH2,I)),1.)
3222                   CALL CHF1(400+IDHIST+1,
3223      &                        SNGL(Y1-YM(ICH2,I)),1.)
3224                   CALL CHF1(400+IDHIST+2,
3225      &                     SNGL(PHI1-PHM(ICH2,I)),1.)
3226                   CALL CHF1(400+IDHIST+3,
3227      &                     SNGL(ALAM1-ALM(ICH2,I)),1.)
3228                   CALL CHF1(IDHIST+4,SNGL(P1),R2)
3229                ENDIF
3230             ENDIF
3231          ENDIF
3232       ENDDO
3233
3234       IF (JOK.GT.0) THEN
3235          IF (ITCHECK(ITRACK(IP(ICH1,IHIT1))).EQ.1) THEN
3236             IF (IDEBUG.GE.2) THEN
3237                IF (IHIT2.EQ.0) THEN
3238                   PRINT *,'CHECK4 histo nb:',IDHIST
3239                   PRINT *,'CHECK4 p de st.',ICH1,'=',P1
3240                   PRINT *,'CHECK4 track not found in st.',ICH2
3241                   PRINT *,'CHECK4 error X :',(XM(ICH2,JOK)-X1), EXM
3242                   PRINT *,'CHECK4 error Y :',(YM(ICH2,JOK)-Y1), EYM
3243                   PRINT *,'CHECK4 error PHI :',(PHM(ICH2,JOK)-PHI1),EPH
3244                   PRINT *,'CHECK4 error ALAM :',(ALM(ICH2,JOK)-ALAM1),
3245      &                     EAL
3246                ELSEIF (IHIT2.NE.JOK) THEN
3247                   PRINT *,'CHECK4 histo nb:',IDHIST
3248                   PRINT *,'CHECK4 p de st.',ICH1,'=',P1
3249                   PRINT *,'CHECK4 ghost in st.',ICH2
3250                   PRINT *,'CHECK4 id part. recherchee:',
3251      &                    ID(IP(ICH1,IHIT1))
3252                   PRINT *,'CHECK4 id ghost trouve    :',
3253      &                    ID(IP(ICH2,IHIT2))
3254                ENDIF
3255             ENDIF
3256          ENDIF
3257       ENDIF
3258
3259       RETURN
3260       END
3261 *****************************************************************************
3262       SUBROUTINE CHECK_HISTO2(IDHIST,ICH2,IHIT2,ICH1,IHIT1,
3263      &                        X1,Y1,X2,Y2,P1,EX1,EY1,EX2,EY2)
3264 *****************************************************************************
3265 *   Check hit IHIT2 with GEANT informations from hit HIT1
3266 *
3267 *    INPUT : IDHIST : No histo
3268 *            ICH2 : No st. de recherche
3269 *            IDCH1,X1,Y1,PHI1,ALAM1,P1 : Trace de reference
3270 *    OUTPUT : JOK : No hit dans ICH2 appartenant a la meme trace.
3271 *
3272 *****************************************************************************
3273       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
3274
3275       PARAMETER (MAXHITCH=10000,MAXHITTOT=20000,NBSTATION=5,
3276      &           MAXTRK=50000)
3277
3278       COMMON/CHHIT/XM(NBSTATION,MAXHITCH),YM(NBSTATION,MAXHITCH),
3279      &             PHM(NBSTATION,MAXHITCH),ALM(NBSTATION,MAXHITCH),
3280      &             IZM(NBSTATION,MAXHITCH),
3281      &             IP(NBSTATION,MAXHITCH),JHIT(NBSTATION),
3282      &             XMR(NBSTATION,MAXHITCH,2),YMR(NBSTATION,MAXHITCH,2)
3283
3284       COMMON/RHIT/ITYP(MAXHITTOT),XTR(MAXHITTOT),YTR(MAXHITTOT),
3285      &            PTOT(MAXHITTOT),ID(MAXHITTOT),IZST(MAXHITTOT),
3286      &            PVERT1(MAXHITTOT),PVERT2(MAXHITTOT),PVERT3(MAXHITTOT),
3287      &            ZVERT(MAXHITTOT),NHITTOT
3288
3289       COMMON/VERIFGEANT/ITTROUGH(MAXTRK,NBSTATION),
3290      &                  IT_LIST(MAXTRK),IT_NP(MAXTRK),ITCHECK(MAXTRK),
3291      &                  ITRACK(MAXHITTOT)
3292 *
3293       COMMON/DEBEVT/IDEBUG
3294
3295       REAL*4 R2
3296       DATA R2/1./
3297
3298       JOK = 0
3299
3300       DO I=1,JHIT(ICH2)
3301          IF (ID(IP(ICH1,IHIT1)).EQ.ID(IP(ICH2,I))) THEN
3302             JOK = I
3303             IF (IDHIST.GT.0) THEN
3304                IF (IZM(ICH2,I).EQ.1) THEN
3305                   X = X1
3306                   Y = Y1
3307                ELSE
3308                   X = X2
3309                   Y = Y2
3310                ENDIF
3311                CALL CHFILL2(200+IDHIST,SNGL(P1),
3312      &                       SNGL((X-XM(ICH2,I))**2),1.)
3313                CALL CHFILL2(200+IDHIST+1,SNGL(P1),
3314      &                       SNGL((Y-YM(ICH2,I))**2),1.)
3315                CALL CHF1(200+IDHIST+2,
3316      &                       SNGL(X-XM(ICH2,I)),1.)
3317                CALL CHF1(200+IDHIST+3,
3318      &                       SNGL(Y-YM(ICH2,I)),1.)
3319                CALL CHF1(200+IDHIST+4,SNGL(P1),R2)
3320             ENDIF
3321          ENDIF
3322       ENDDO
3323
3324       IF (JOK.GT.0) THEN
3325          IF (ITCHECK(ITRACK(IP(ICH1,IHIT1))).EQ.1) THEN
3326             EXM = EX1
3327             EYM = EY1
3328             IF (IZM(ICH2,JOK).EQ.1) THEN
3329                X = X1
3330                Y = Y1
3331                IF (ICH2.EQ.4.OR.ICH2.EQ.5) THEN
3332                   EXM = EX1
3333                   EYM = EY1
3334                ENDIF
3335             ELSE
3336                X = X2
3337                Y = Y2
3338                IF (ICH2.EQ.4.OR.ICH2.EQ.5) THEN
3339                   EXM = EX2
3340                   EYM = EY2
3341                ENDIF
3342             ENDIF
3343             IF (IDEBUG.GE.2) THEN
3344                IF (IHIT2.EQ.0) THEN
3345                   PRINT *,'CHECK2 histo nb:',IDHIST
3346                   PRINT *,'CHECK2 p de st.',ICH1,'=',P1
3347                   PRINT *,'CHECK2 track not found in st.',ICH2
3348                   PRINT *,'CHECK2 error X :',(XM(ICH2,JOK)-X), EXM
3349                   PRINT *,'CHECK2 error Y :',(YM(ICH2,JOK)-Y), EYM
3350                ELSEIF(IHIT2.NE.JOK) THEN
3351                   PRINT *,'CHECK2 histo nb:',IDHIST
3352                   PRINT *,'CHECK2 p de st.',ICH1,'=',P1
3353                   PRINT *,'CHECK2 ghost in st.',ICH2
3354                   PRINT *,'CHECK2 id part. recherchee:',
3355      &                    ID(IP(ICH1,IHIT1))
3356                   PRINT *,'CHECK2 id ghost trouve    :',
3357      &                    ID(IP(ICH2,IHIT2))
3358                ENDIF
3359             ENDIF
3360          ENDIF
3361       ENDIF
3362
3363       RETURN
3364       END
3365
3366
3367
3368 ************************************************************************
3369       DOUBLE PRECISION FUNCTION DEDX (P,THET,XEA,YEA)
3370 ************************************************************************
3371 *    DEDX est la nouvelle impulsion au vertex, corrigee de la perte
3372 *    d'energie dans l'absorbeur
3373 *
3374 ************************************************************************
3375
3376       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3377
3378       REA = DSQRT(XEA**2+YEA**2)
3379       IF (REA.GT.26.3611) THEN
3380 *   Plomb
3381          SPB = 5./DCOS(THET)
3382          P = P+SPB/1000.*11.35*(16.66D-3 * P+1.33)
3383 *   Polyethylene
3384          SPO = 5./DCOS(THET)
3385          P = P+SPO/1000.*0.935*(2.22D-3 * P+2.17)
3386 *   Plomb
3387          SPB = 5./DCOS(THET)
3388          P = P+SPB/1000.*11.35*(16.66D-3 * P+1.33)
3389 *   Polyethylene
3390          SPO = 5./DCOS(THET)
3391          P = P+SPO/1000.*0.935*(2.22D-3 * P+2.17)
3392 *   Plomb
3393          SPB = 5./DCOS(THET)
3394          P = P+SPB/1000.*11.35*(16.66D-3 * P+1.33)
3395 *   Polyethylene
3396          SPO = 5./DCOS(THET)
3397          P = P+SPO/1000.*0.935*(2.22D-3 * P+2.17)
3398 *   Plomb
3399          SPB = 5./DCOS(THET)
3400          P = P+SPB/1000.*11.35*(16.66D-3 * P+1.33)
3401
3402       ELSE
3403 *   Tungstene
3404          SW = (503.-467.)/DCOS(THET)
3405          P = P+SW/1000.*19.3*(16.66D-3 * P+1.33)
3406       ENDIF
3407
3408 *   Concrete
3409       SCONC = (467.-315.)/DCOS(THET)
3410       P = P+SCONC/1000.*2.5*(2.22D-3*P+2.17)
3411
3412 *   Carbone
3413       SC = (315.-90.)/DCOS(THET)
3414       P = P+SC/1000.*1.93*(2.22D-3*P+2.17) ! Carbone
3415
3416       DEDX = P
3417
3418       RETURN
3419
3420       END
3421 ************************************************************************
3422       SUBROUTINE BRANSON(PXZ,PHI,ALAM,XEA,YEA)
3423 ************************************************************************
3424 *
3425 *   Correction de Branson du multiple scattering dans l'absorbeur
3426 *
3427 ************************************************************************
3428
3429       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3430
3431       PARAMETER(NPLANE=10)
3432       COMMON/PARAM/ZPLANEP(NPLANE),THICK,XPREC,YPREC,B0,BL3,ZMAGS,
3433      &             ZMAGE,ZABS,XMAG,ZBP1,ZBP2,CONST
3434
3435       ASIGN = 1.
3436       IF (PXZ.LT.0.) ASIGN = -1.
3437       PXZ = DABS(PXZ)
3438       PX = PXZ*DSIN(PHI)
3439       PY = PXZ*DTAN(ALAM)
3440       PZ = PXZ*DCOS(PHI)
3441
3442       PTOT = PXZ/DCOS(ALAM)
3443
3444       ZEA = ZABS
3445
3446       REA = DSQRT(XEA**2+YEA**2)
3447       IF (REA.GT.26.3611) THEN
3448          ZBP = ZBP1
3449       ELSE            ! Abso. W  pour theta < 3 deg
3450          ZBP = ZBP2
3451       ENDIF
3452
3453       XBP = XEA-PX/PZ*(ZEA-ZBP)
3454       YBP = YEA-PY/PZ*(ZEA-ZBP)
3455       PZ = PTOT*ZBP/DSQRT(XBP**2+YBP**2+ZBP**2)
3456       PX = PZ*XBP/ZBP
3457       PY = PZ*YBP/ZBP
3458       PXZ = DSQRT(PX**2+PZ**2)
3459       PHI = DATAN2(PX,PZ)
3460       ALAM = DATAN2(PY,PXZ)
3461
3462 **      THET = DATAN2(REA,ZEA)
3463
3464       PT = DSQRT(PX**2+PY**2)
3465       THET = DATAN2(PT,PZ)
3466       PTOT =  DEDX(PTOT,THET,XEA,YEA)
3467
3468       PXZ = ASIGN*PTOT*DCOS(ALAM)
3469
3470       RETURN
3471       END
3472
3473 ***************************************************************
3474       SUBROUTINE DSINV(N,A,IDIM,IFAIL)
3475 ***************************************************************
3476
3477       DOUBLE PRECISION    A(IDIM,*),  ZERO,  ONE,  X, Y
3478
3479       REAL                PIVOTF
3480       CHARACTER*6         HNAME
3481
3482       DOUBLE PRECISION    S1, S31, S32, S33,  DOTF
3483
3484       PIVOTF(X)    =  SNGL(X)
3485       DOTF(X,Y,S1)  =  X * Y + S1
3486
3487       DATA      HNAME               /  'DSINV '  /
3488       DATA      ZERO, ONE           /  0.D0, 1.D0 /
3489
3490       IF(IDIM .LT. N  .OR.  N .LE. 0)  GOTO 900
3491 *
3492 * sfact.inc
3493 *
3494       IFAIL  =  0
3495       DO 144    J  =  1, N
3496          IF(PIVOTF(A(J,J)) .LE. 0.)  GOTO 150
3497          A(J,J)  =  ONE / A(J,J)
3498          IF(J .EQ. N)  GOTO 199
3499  140     JP1  =  J+1
3500          DO 143   L  =  JP1, N
3501             A(J,L)  =  A(J,J)*A(L,J)
3502             S1      =  -A(L,J+1)
3503             DO 141  I  =  1, J
3504                S1  =  DOTF(A(L,I),A(I,J+1),S1)
3505  141        CONTINUE
3506             A(L,J+1)  =  -S1
3507  143     CONTINUE
3508  144  CONTINUE
3509  150  IFAIL  =  -1
3510       RETURN
3511  199  CONTINUE
3512 *
3513 * sfinv.inc
3514 *
3515       IF(N .EQ. 1)  GOTO 399
3516       A(1,2)  =  -A(1,2)
3517       A(2,1)  =   A(1,2)*A(2,2)
3518       IF(N .EQ. 2)  GOTO 320
3519       DO 314    J  =  3, N
3520          JM2  =  J - 2
3521          DO 312 K  =  1, JM2
3522             S31  =  A(K,J)
3523             DO 311  I  =  K, JM2
3524                S31  =  DOTF(A(K,I+1),A(I+1,J),S31)
3525  311        CONTINUE
3526             A(K,J)  =  -S31
3527             A(J,K)  =  -S31*A(J,J)
3528  312     CONTINUE
3529          A(J-1,J)  =  -A(J-1,J)
3530          A(J,J-1)  =   A(J-1,J)*A(J,J)
3531  314  CONTINUE
3532  320  J  =  1
3533  323  S33  =  A(J,J)
3534       IF(J .EQ. N)  GOTO 325
3535       JP1  =  J + 1
3536       DO 324 I  =  JP1, N
3537          S33  =  DOTF(A(J,I),A(I,J),S33)
3538  324  CONTINUE
3539  325  A(J,J)  =  S33
3540       JM1  =  J
3541       J    =  JP1
3542       DO 328 K  =  1, JM1
3543          S32  =  ZERO
3544          DO 327  I  =  J, N
3545             S32  =  DOTF(A(K,I),A(I,J),S32)
3546  327     CONTINUE
3547          A(K,J)  =  S32
3548          A(J,K)  =  S32
3549  328  CONTINUE
3550       IF(J .LT. N)  GOTO 323
3551  399  CONTINUE
3552
3553       RETURN
3554  900  CALL TMPRNT(HNAME,N,IDIM,0)
3555       RETURN
3556       END
3557
3558 *******************************************************
3559       SUBROUTINE TMPRNT(NAME,N,IDIM,K)
3560 *******************************************************
3561
3562       CHARACTER*6         NAME
3563       LOGICAL             MFLAG,    RFLAG
3564
3565       IF(NAME(2:2) .EQ. 'S') THEN
3566          CALL KERMTR('F012.1',LGFILE,MFLAG,RFLAG)
3567       ELSE
3568          CALL KERMTR('F011.1',LGFILE,MFLAG,RFLAG)
3569       ENDIF
3570       IF(NAME(3:6) .EQ. 'FEQN') ASSIGN 1002 TO IFMT
3571       IF(NAME(3:6) .NE. 'FEQN') ASSIGN 1001 TO IFMT
3572       IF(MFLAG) THEN
3573          IF(LGFILE .EQ. 0) THEN
3574             IF(NAME(3:6) .EQ. 'FEQN') THEN
3575                WRITE(*,IFMT) NAME, N, IDIM, K
3576             ELSE
3577                WRITE(*,IFMT) NAME, N, IDIM
3578             ENDIF
3579          ELSE
3580             IF(NAME(3:6) .EQ. 'FEQN') THEN
3581                WRITE(LGFILE,IFMT) NAME, N, IDIM, K
3582             ELSE
3583                WRITE(LGFILE,IFMT) NAME, N, IDIM
3584             ENDIF
3585          ENDIF
3586       ENDIF
3587       IF(.NOT. RFLAG) CALL ABEND
3588       RETURN
3589  1001 FORMAT(7X, 31H PARAMETER ERROR IN SUBROUTINE , A6,
3590      +     27H ... (N.LT.1 OR IDIM.LT.N).,
3591      +             5X, 3HN =, I4, 5X, 6HIDIM =, I4, 1H. )
3592  1002 FORMAT(7X, 31H PARAMETER ERROR IN SUBROUTINE , A6,
3593      +     37H ... (N.LT.1 OR IDIM.LT.N OR K.LT.1).,
3594      +     5X, 3HN =, I4, 5X, 6HIDIM =, I4, 5X, 3HK =, I4,1H.)
3595       END
3596 *
3597 * $Id$
3598 *
3599 * $Log$
3600 * Revision 1.1  1999/07/30 10:53:20  fca
3601 * New version of MUON from M.Bondila & A.Morsch
3602 *
3603 * Revision 1.1.1.1  1996/02/15 17:48:35  mclareni
3604 * Kernlib
3605 *
3606 *
3607
3608 ***********************************************************
3609       SUBROUTINE KERSET(ERCODE,LGFILE,LIMITM,LIMITR)
3610 ***********************************************************
3611
3612       PARAMETER(KOUNTE  =  27)
3613       CHARACTER*6         ERCODE,   CODE(KOUNTE)
3614       LOGICAL             MFLAG,    RFLAG
3615       INTEGER             KNTM(KOUNTE),       KNTR(KOUNTE)
3616
3617       DATA      LOGF      /  0  /
3618       DATA      CODE(1), KNTM(1), KNTR(1)  / 'C204.1', 255, 255 /
3619       DATA      CODE(2), KNTM(2), KNTR(2)  / 'C204.2', 255, 255 /
3620       DATA      CODE(3), KNTM(3), KNTR(3)  / 'C204.3', 255, 255 /
3621       DATA      CODE(4), KNTM(4), KNTR(4)  / 'C205.1', 255, 255 /
3622       DATA      CODE(5), KNTM(5), KNTR(5)  / 'C205.2', 255, 255 /
3623       DATA      CODE(6), KNTM(6), KNTR(6)  / 'C305.1', 255, 255 /
3624       DATA      CODE(7), KNTM(7), KNTR(7)  / 'C308.1', 255, 255 /
3625       DATA      CODE(8), KNTM(8), KNTR(8)  / 'C312.1', 255, 255 /
3626       DATA      CODE(9), KNTM(9), KNTR(9)  / 'C313.1', 255, 255 /
3627       DATA      CODE(10),KNTM(10),KNTR(10) / 'C336.1', 255, 255 /
3628       DATA      CODE(11),KNTM(11),KNTR(11) / 'C337.1', 255, 255 /
3629       DATA      CODE(12),KNTM(12),KNTR(12) / 'C341.1', 255, 255 /
3630       DATA      CODE(13),KNTM(13),KNTR(13) / 'D103.1', 255, 255 /
3631       DATA      CODE(14),KNTM(14),KNTR(14) / 'D106.1', 255, 255 /
3632       DATA      CODE(15),KNTM(15),KNTR(15) / 'D209.1', 255, 255 /
3633       DATA      CODE(16),KNTM(16),KNTR(16) / 'D509.1', 255, 255 /
3634       DATA      CODE(17),KNTM(17),KNTR(17) / 'E100.1', 255, 255 /
3635       DATA      CODE(18),KNTM(18),KNTR(18) / 'E104.1', 255, 255 /
3636       DATA      CODE(19),KNTM(19),KNTR(19) / 'E105.1', 255, 255 /
3637       DATA      CODE(20),KNTM(20),KNTR(20) / 'E208.1', 255, 255 /
3638       DATA      CODE(21),KNTM(21),KNTR(21) / 'E208.2', 255, 255 /
3639       DATA      CODE(22),KNTM(22),KNTR(22) / 'F010.1', 255,   0 /
3640       DATA      CODE(23),KNTM(23),KNTR(23) / 'F011.1', 255,   0 /
3641       DATA      CODE(24),KNTM(24),KNTR(24) / 'F012.1', 255,   0 /
3642       DATA      CODE(25),KNTM(25),KNTR(25) / 'F406.1', 255,   0 /
3643       DATA      CODE(26),KNTM(26),KNTR(26) / 'G100.1', 255, 255 /
3644       DATA      CODE(27),KNTM(27),KNTR(27) / 'G100.2', 255, 255 /
3645
3646       LOGF  =  LGFILE
3647       L  =  0
3648       IF(ERCODE .NE. ' ')  THEN
3649          DO 10  L = 1, 6
3650             IF(ERCODE(1:L) .EQ. ERCODE)  GOTO 12
3651  10      CONTINUE
3652  12      CONTINUE
3653       ENDIF
3654       DO 14     I  =  1, KOUNTE
3655          IF(L .EQ. 0)  GOTO 13
3656          IF(CODE(I)(1:L) .NE. ERCODE(1:L))  GOTO 14
3657  13      IF(LIMITM.GE.0) KNTM(I)  =  LIMITM
3658          IF(LIMITR.GE.0) KNTR(I)  =  LIMITR
3659  14   CONTINUE
3660       RETURN
3661       ENTRY KERMTR(ERCODE,LOG,MFLAG,RFLAG)
3662       LOG  =  LOGF
3663       DO 20     I  =  1, KOUNTE
3664          IF(ERCODE .EQ. CODE(I))  GOTO 21
3665  20   CONTINUE
3666       WRITE(*,1000)  ERCODE
3667       CALL ABEND
3668       RETURN
3669  21   RFLAG  =  KNTR(I) .GE. 1
3670       IF(RFLAG  .AND.  (KNTR(I) .LT. 255))  KNTR(I)  =  KNTR(I) - 1
3671       MFLAG  =  KNTM(I) .GE. 1
3672       IF(MFLAG  .AND.  (KNTM(I) .LT. 255))  KNTM(I)  =  KNTM(I) - 1
3673       IF(.NOT. RFLAG)  THEN
3674          IF(LOGF .LT. 1)  THEN
3675             WRITE(*,1001)  CODE(I)
3676          ELSE
3677             WRITE(LOGF,1001)  CODE(I)
3678          ENDIF
3679       ENDIF
3680       IF(MFLAG .AND. RFLAG)  THEN
3681          IF(LOGF .LT. 1)  THEN
3682             WRITE(*,1002)  CODE(I)
3683          ELSE
3684             WRITE(LOGF,1002)  CODE(I)
3685          ENDIF
3686       ENDIF
3687       RETURN
3688  1000 FORMAT(' KERNLIB LIBRARY ERROR. ' /
3689      +     ' ERROR CODE ',A6,' NOT RECOGNIZED BY KERMTR',
3690      +     ' ERROR MONITOR. RUN ABORTED.')
3691  1001 FORMAT(/' ***** RUN TERMINATED BY CERN LIBRARY ERROR ',
3692      +     'CONDITION ',A6)
3693  1002 FORMAT(/' ***** CERN LIBRARY ERROR CONDITION ',A6)
3694       END
3695
3696 **************************
3697       subroutine abend
3698 **************************
3699
3700       stop 'abend!'
3701       end
3702
3703 ************************************************************************
3704       SUBROUTINE FCNFIT(NPAR, GRAD, FVAL, XVAL, IFLAG, FUTIL)
3705 ************************************************************************
3706 *    With magnetic Field Map GRKUTA
3707 *
3708 *    Calcule FVAL: la fonction minimisee par MINUIT
3709 *    With magnetic field map
3710 *
3711 ************************************************************************
3712       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3713
3714       PARAMETER(NPLANE=10)
3715
3716       COMMON/PARAM/ZPLANEP(NPLANE),THICK,XPREC,YPREC,B0,BL3,ZMAGS,
3717      &             ZMAGE,ZABS,XMAG,ZBP1,ZBP2,CONST
3718
3719       COMMON/MEAS/LPLANE(NPLANE),XMP(NPLANE),YMP(NPLANE)
3720
3721       COMMON/FCNOUT/PXZEA,ALAMEA,PHIEA,XEA,YEA,NPLU,CHI2
3722
3723       DIMENSION GRAD(*),XVAL(*),AMS(500),DISTAZ(500)
3724
3725       DIMENSION XP(NPLANE),YP(NPLANE),
3726      &          COV(NPLANE,NPLANE),AP(NPLANE),COVY(NPLANE,NPLANE)
3727       DIMENSION VECT(7),VOUT(7)
3728
3729       STEP = 2.   ! 1 cm
3730       NSTEPMAX = 5000
3731       PITODEG = 57.295
3732       XV = XVAL(4)
3733       YV = XVAL(5)
3734
3735       ASIGN = 1.
3736       IF (XVAL(1).LT.0.) ASIGN = -1.
3737       PHI  = XVAL(2)
3738       ALAM = XVAL(3)
3739       PXZ = DABS(1./XVAL(1))
3740
3741       PX = PXZ*DSIN(PHI)
3742       PY = PXZ*DTAN(ALAM)
3743       PZ = PXZ*DCOS(PHI)
3744       PTOT = PXZ/DCOS(ALAM)
3745
3746       A12 = 0.
3747       NTMAX = 0
3748
3749       ZEA = ZABS
3750       XEA = XV
3751       YEA = YV
3752       PXEA = PX
3753       PYEA = PY
3754       PHIEA = PHI
3755       PXZEA = ASIGN*PXZ
3756       ALAMEA = ALAM
3757
3758       VECT(1) = XV
3759       VECT(2) = YV
3760       VECT(3) = ZABS
3761       VECT(4) = PX/PTOT
3762       VECT(5) = PY/PTOT
3763       VECT(6) = PZ/PTOT
3764       VECT(7) = PTOT
3765
3766       R = SQRT(VECT(1)*VECT(1)+VECT(2)*VECT(2))
3767
3768       Z = VECT(3)
3769       NSTEP = 0
3770       IX = 0
3771       IY = 0
3772       IZ = 0
3773 **      PRINT *,' AV GRKUTA ASIGN',ASIGN,' THET',THET
3774       DO ICH = 1,NPLANE
3775          DO WHILE (Z.GE.ZABS.AND.Z.LT.ZPLANEP(ICH)
3776      &            .AND.NSTEP.LE.NSTEPMAX)
3777 **     &            .AND.(THETA*PITODEG).GT.2.
3778 **     &            .AND. (THETA*PITODEG).LT.9.)
3779             NSTEP = NSTEP+1
3780 **          WRITE(6,*) NSTEP,(VECT(I),I=1,7)
3781             CALL RECO_GRKUTA(ASIGN,STEP,VECT,VOUT)
3782             DO I = 1,7
3783                VECT(I) = VOUT(I)
3784             ENDDO
3785             Z = VECT(3)
3786             R = SQRT(VECT(1)*VECT(1)+VECT(2)*VECT(2))
3787          ENDDO
3788          IF (NSTEP.EQ.NSTEPMAX) RETURN
3789          XP(ICH) = VECT(1)-(Z-ZPLANEP(ICH))*VECT(4)/VECT(6)
3790          YP(ICH) = VECT(2)-(Z-ZPLANEP(ICH))*VECT(5)/VECT(6)
3791          AL   = THICK/ VECT(6)
3792          AP(ICH) = (0.0136D0/PTOT)*DSQRT(AL)*(1+0.038D0*DLOG(AL))
3793       ENDDO
3794 **    PRINT *,' AP GRKUTA ASIGN',ASIGN,' THET',THET
3795
3796
3797 ** Matrice de covariance
3798       I = 0
3799       DO II = 1,NPLANE
3800         IF (LPLANE(II).EQ.1) THEN
3801            I = I + 1
3802 *     I = II
3803            J = I - 1
3804            DO JJ = II, NPLANE
3805               IF (LPLANE(JJ).EQ.1) THEN
3806                  J = J + 1
3807 *     J = JJ
3808                  COV (I,J) = 0.0D0
3809                  COV (J,I) = A12
3810                  IF (I .EQ. J) THEN
3811                     COV(J,I) =COV(J,I) + XPREC**2
3812                  ENDIF
3813
3814 *     IF (I .EQ. 10 .AND. J .EQ. 10) PRINT *,'10 10   ',COV(J,I)
3815                  DO L = 1,NTMAX
3816                     COV(J,I) = COV(J,I)
3817      &                   +  (ZPLANEP(II) + DISTAZ(L))*(ZPLANEP(JJ) +
3818      &                   DISTAZ(L))*AMS(L)**2
3819                  ENDDO
3820                  DO K = 1, II-1
3821                     COV(J,I) = COV(J,I)
3822      &                   +  (ZPLANEP(II)-ZPLANEP(K))*
3823      &                   (ZPLANEP(JJ)-ZPLANEP(K))*AP(K)**2
3824 *     IF (I .EQ. 10 .AND. J .EQ. 10) PRINT *,'10 10   ',COV(J,I)
3825                  ENDDO
3826                  COVY(I,J) = 0.0D0
3827                  COVY(J,I) = COV(J,I)
3828                  IF (I .EQ. J) THEN
3829                     COVY(J,I) = COVY(J,I) - XPREC**2 + YPREC**2
3830                  ENDIF
3831               ENDIF
3832            ENDDO
3833         ENDIF
3834       ENDDO
3835
3836 *  Inversion des matrices de covariance
3837       NPLU = I
3838
3839       IFAIL = 0
3840       CALL DSINV(NPLU, COV, NPLANE, IFAIL)
3841 **      IF (JFAIL.NE.0 .AND. IFAIL .NE. 0) STOP 'ERROR'
3842       IF (IFAIL .NE. 0) STOP 'ERROR'
3843       IFAIL = 0
3844       CALL DSINV(NPLU, COVY, NPLANE, IFAIL)
3845 **      IF (JFAIL.NE.0 .AND. IFAIL .NE. 0) STOP 'ERROR'
3846       IF (IFAIL .NE. 0) STOP 'ERROR'
3847 *      PRINT *,' COVARIANCE MATRIX AFTER'
3848 *      DO I = 1, NPLANE
3849 *         PRINT *,(COV(J,I),J=1,NPLANE)
3850 *      ENDDO
3851
3852 ** Calcul de FVAL ou CHI2
3853       FVAL = 0.0D0
3854       I = 0
3855       DO II = 1,NPLANE
3856         IF (LPLANE(II).EQ.1) THEN
3857         I = I+1
3858 *        I = II
3859         J = 0
3860         DO JJ = 1,NPLANE
3861            IF (LPLANE(JJ).EQ.1) THEN
3862               J = J+1
3863 *             J = JJ
3864               FVAL = FVAL + COV(J,I)*(XMP(II)-XP(II))*(XMP(JJ)-XP(JJ))
3865               FVAL = FVAL + COVY(J,I)*(YMP(II)-YP(II))
3866      &                               *(YMP(JJ)-YP(JJ))
3867 **             IF (JJ.EQ.II) THEN
3868 **                 FVAL = FVAL + (XM(II)-XP(II))*(XM(JJ)-XP(JJ))/XPREC**2
3869 **                 FVAL = FVAL + (YM(II)-YP(II))
3870 **     &                               *(YM(JJ)-YP(JJ))/YPREC**2
3871 **             ENDIF
3872            ENDIF
3873         ENDDO
3874         ENDIF
3875       ENDDO
3876       CHI2 = FVAL
3877
3878 **      IF (CHI2.GT.1.E4) THEN
3879 **         PRINT *,'FCNFIT CHI2= ',CHI2
3880 **         FVAL = 0.
3881 **      ENDIF
3882
3883
3884  1000 FORMAT(I5,7F12.6)
3885
3886       RETURN
3887       END
3888
3889 ************************************************************************
3890       SUBROUTINE NEWFCNFIT(NPAR, GRAD, FVAL, XVAL, IFLAG, FUTIL)
3891 ************************************************************************
3892 *    With magnetic Field Map GRKUTA
3893 *      trackfinding
3894 *    Calcule FVAL: la fonction minimisee par MINUIT
3895 *    With magnetic field map
3896 *
3897 ************************************************************************
3898       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3899
3900       PARAMETER(NPLANE=10)
3901
3902       COMMON/PARAM/ZPLANEP(NPLANE),THICK,XPREC,YPREC,B0,BL3,ZMAGS,
3903      &             ZMAGE,ZABS,XMAG,ZBP1,ZBP2,CONST
3904
3905       COMMON/MEAS/LPLANE(NPLANE),XMP(NPLANE),YMP(NPLANE)
3906
3907       COMMON/FCNOUT/PXZEA,ALAMEA,PHIEA,XEA,YEA,NPLU,CHI2
3908
3909       DIMENSION GRAD(*),XVAL(*),AMS(500),DISTAZ(500)
3910
3911       DIMENSION XP(NPLANE),YP(NPLANE),
3912      &          COV(NPLANE,NPLANE),AP(NPLANE),COVY(NPLANE,NPLANE)
3913       DIMENSION VECT(7),VOUT(7)
3914
3915       STEP = 2.   ! 1 cm
3916       NSTEPMAX = 5000
3917       PITODEG = 57.295
3918       XV = XVAL(4)
3919       YV = XVAL(5)
3920
3921       ASIGN = 1.
3922       IF (XVAL(1).LT.0.) ASIGN = -1.
3923       PHI  = XVAL(2)
3924       ALAM = XVAL(3)
3925       PXZ = DABS(1./XVAL(1))
3926
3927       PX = PXZ*DSIN(PHI)
3928       PY = PXZ*DTAN(ALAM)
3929       PZ = PXZ*DCOS(PHI)
3930       PTOT = PXZ/DCOS(ALAM)
3931
3932       A12 = 0.
3933       NTMAX = 0
3934
3935       ZEA = ZABS
3936       XEA = XV
3937       YEA = YV
3938       PXEA = PX
3939       PYEA = PY
3940       PHIEA = PHI
3941       PXZEA = ASIGN*PXZ
3942       ALAMEA = ALAM
3943
3944       VECT(1) = XV
3945       VECT(2) = YV
3946       VECT(3) = ZABS
3947       VECT(4) = PX/PTOT
3948       VECT(5) = PY/PTOT
3949       VECT(6) = PZ/PTOT
3950       VECT(7) = PTOT
3951
3952       R = SQRT(VECT(1)*VECT(1)+VECT(2)*VECT(2))
3953
3954       Z = VECT(3)
3955       NSTEP = 0
3956       IX = 0
3957       IY = 0
3958       IZ = 0
3959 **      PRINT *,' AV GRKUTA ASIGN',ASIGN,' THET',THET
3960       DO ICH = 1,NPLANE
3961          DO WHILE (Z.GE.ZABS.AND.Z.LT.ZPLANEP(ICH)
3962      &            .AND.NSTEP.LE.NSTEPMAX)
3963 **     &            .AND.(THETA*PITODEG).GT.2.
3964 **     &            .AND. (THETA*PITODEG).LT.9.)
3965             NSTEP = NSTEP+1
3966 **          WRITE(6,*) NSTEP,(VECT(I),I=1,7)
3967             CALL RECO_GRKUTA (ASIGN,STEP,VECT,VOUT)
3968             DO I = 1,7
3969                VECT(I) = VOUT(I)
3970             ENDDO
3971             Z = VECT(3)
3972             R = SQRT(VECT(1)*VECT(1)+VECT(2)*VECT(2))
3973          ENDDO
3974          IF (NSTEP.EQ.NSTEPMAX) RETURN
3975          XP(ICH) = VECT(1)-(Z-ZPLANEP(ICH))*VECT(4)/VECT(6)
3976          YP(ICH) = VECT(2)-(Z-ZPLANEP(ICH))*VECT(5)/VECT(6)
3977          AL   = THICK/ VECT(6)
3978          AP(ICH) = (0.0136D0/PTOT)*DSQRT(AL)*(1+0.038D0*DLOG(AL))
3979       ENDDO
3980 **    PRINT *,' AP GRKUTA ASIGN',ASIGN,' THET',THET
3981
3982
3983 ** Matrice de covariance
3984       I = 0
3985       DO II = 1,NPLANE
3986         IF (LPLANE(II).EQ.1) THEN
3987            I = I + 1
3988 *     I = II
3989            J = I - 1
3990            DO JJ = II, NPLANE
3991               IF (LPLANE(JJ).EQ.1) THEN
3992                  J = J + 1
3993 *     J = JJ
3994                  COV (I,J) = 0.0D0
3995                  COV (J,I) = A12
3996                  IF (I .EQ. J) THEN
3997                     COV(J,I) =COV(J,I) + XPREC**2
3998                  ENDIF
3999
4000 *     IF (I .EQ. 10 .AND. J .EQ. 10) PRINT *,'10 10   ',COV(J,I)
4001                  DO L = 1,NTMAX
4002                     COV(J,I) = COV(J,I)
4003      &                   +  (ZPLANEP(II) + DISTAZ(L))*(ZPLANEP(JJ) +
4004      &                   DISTAZ(L))*AMS(L)**2
4005                  ENDDO
4006                  DO K = 1, II-1
4007                     COV(J,I) = COV(J,I)
4008      &                   +  (ZPLANEP(II)-ZPLANEP(K))*
4009      &                   (ZPLANEP(JJ)-ZPLANEP(K))*AP(K)**2
4010 *     IF (I .EQ. 10 .AND. J .EQ. 10) PRINT *,'10 10   ',COV(J,I)
4011                  ENDDO
4012                  COVY(I,J) = 0.0D0
4013                  COVY(J,I) = COV(J,I)
4014                  IF (I .EQ. J) THEN
4015                     COVY(J,I) = COVY(J,I) - XPREC**2 + YPREC**2
4016                  ENDIF
4017               ENDIF
4018            ENDDO
4019         ENDIF
4020       ENDDO
4021
4022 *  Inversion des matrices de covariance
4023       NPLU = I
4024
4025       IFAIL = 0
4026       CALL DSINV(NPLU, COV, NPLANE, IFAIL)
4027 **      IF (JFAIL.NE.0 .AND. IFAIL .NE. 0) STOP 'ERROR'
4028       IF (IFAIL .NE. 0) STOP 'ERROR'
4029       IFAIL = 0
4030       CALL DSINV(NPLU, COVY, NPLANE, IFAIL)
4031 **      IF (JFAIL.NE.0 .AND. IFAIL .NE. 0) STOP 'ERROR'
4032       IF (IFAIL .NE. 0) STOP 'ERROR'
4033 *      PRINT *,' COVARIANCE MATRIX AFTER'
4034 *      DO I = 1, NPLANE
4035 *         PRINT *,(COV(J,I),J=1,NPLANE)
4036 *      ENDDO
4037
4038 ** Calcul de FVAL ou CHI2
4039       FVAL = 0.0D0
4040       I = 0
4041       DO II = 1,NPLANE
4042         IF (LPLANE(II).EQ.1) THEN
4043         I = I+1
4044 *        I = II
4045         J = 0
4046         DO JJ = 1,NPLANE
4047            IF (LPLANE(JJ).EQ.1) THEN
4048               J = J+1
4049 *             J = JJ
4050               FVAL = FVAL + COV(J,I)*(XMP(II)-XP(II))*(XMP(JJ)-XP(JJ))
4051               FVAL = FVAL + COVY(J,I)*(YMP(II)-YP(II))
4052      &                               *(YMP(JJ)-YP(JJ))
4053 **             IF (JJ.EQ.II) THEN
4054 **                 FVAL = FVAL + (XM(II)-XP(II))*(XM(JJ)-XP(JJ))/XPREC**2
4055 **                 FVAL = FVAL + (YM(II)-YP(II))
4056 **     &                               *(YM(JJ)-YP(JJ))/YPREC**2
4057 **             ENDIF
4058            ENDIF
4059         ENDDO
4060         ENDIF
4061       ENDDO
4062       CHI2 = FVAL
4063
4064 **      IF (CHI2.GT.1.E4) THEN
4065 **         PRINT *,'FCNFIT CHI2= ',CHI2
4066 **         FVAL = 0.
4067 **      ENDIF
4068
4069
4070  1000 FORMAT(I5,7F12.6)
4071
4072       RETURN
4073       END
4074
4075 ***********************************************************************
4076       SUBROUTINE INITFIELDOLD
4077 *
4078 *   Galina
4079 ***********************************************************************
4080
4081       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4082
4083 **      IMPLICIT REAL*8(A-H,O-Z)
4084 **      REAL *4 BX,BY,BZ
4085       COMMON/DAT1/Z(81),X(81),Y(81,44),DX,DZ,LPX,LPY,LPZ
4086       COMMON/DAT2/BX(81,81,44),BY(81,81,44),BZ(81,81,44)
4087       COMMON/REG1/ZMAX,ZMIN,XMAX,XMIN
4088       COMMON/REG2/AY1,CY1,AY2,CY2
4089 **      REAL *4 BXP,BYP,BZP
4090       COMMON/SDAT1/ZP(51),RAD(10),FI(33),DZP,DFI,DR,YY0,LPPZ,NR,NFI
4091       COMMON/SDAT2/BXP(51,10,33),BYP(51,10,33),BZP(51,10,33)
4092       COMMON/SDAT4/B(2,2,32)
4093       COMMON/REG3/ZPMAX,ZPMIN,RMAX,RMIN
4094 cc      COMMON/CONST/PI2,EPS
4095       REWIND 40
4096  1000 FORMAT(5(1X,D15.7))
4097  2000 FORMAT(5(1X,I5))
4098       READ(40,2000) LPX,LPY,LPZ
4099       READ(40,1000) (Z(K),K=1,81)
4100       READ(40,1000) (X(K),K=1,81)
4101       READ(40,1000) DX,DY,DZ
4102       READ(40,1000) ZMAX,ZMIN,XMAX,XMIN
4103 c      write(*,*) 'zmin zmax',ZMIN,ZMAX
4104 c      write(*,*) 'xmin xmax',XMIN,XMAX
4105       READ(40,1000) AY1,CY1,AY2,CY2
4106 c      write(*,*) 'ay1,cy1,ay2,cy2', AY1,CY1,AY2,CY2
4107 cc      READ(40,1000) PI2,EPS
4108       READ(40,1000) (((BX(K,L,M),K=1,81),L=1,81),M=1,44)
4109       READ(40,1000) (((BY(K,L,M),K=1,81),L=1,81),M=1,44)
4110       READ(40,1000) (((BZ(K,L,M),K=1,81),L=1,81),M=1,44)
4111 **      RETURN
4112 **      END
4113 c Polar part
4114       READ(40,2000) LPPZ,NR,NFI
4115       READ(40,1000) (ZP(K),K=1,51)
4116       READ(40,1000) (RAD(K),K=1,10)
4117       READ(40,1000) (FI(L),L=1,33)
4118       READ(40,1000) DZP,DFI,DR
4119 c      write(*,*) 'dzp dfi dR',DZP,DFI,DR
4120       READ(40,1000) ZPMAX,ZPMIN,RMAX,RMIN
4121 c      write(*,*) 'zmin zmax',ZPMIN,ZPMAX
4122 c     write(*,*) 'Rmin Rmax',RMIN,RMAX
4123       READ(40,1000) (((BXP(K,L,M),K=1,51),L=1,10),M=1,33)
4124       READ(40,1000) (((BYP(K,L,M),K=1,51),L=1,10),M=1,33)
4125       READ(40,1000) (((BZP(K,L,M),K=1,51),L=1,10),M=1,33)
4126       READ(40,1000) (((B(K,L,M),K=1,2),L=1,2),M=1,32)
4127 **      RETURN
4128 **      END
4129
4130
4131
4132
4133       RETURN
4134       END
4135
4136 ***********************************************************************
4137       SUBROUTINE RECO_GUFLDOLD(X,F)
4138 C     ^^^^^^^^^^^^^^^^^^^^^^
4139 C   field map G. Chabratova
4140 C
4141 C  Field map 31/05/99
4142 ***********************************************************************
4143
4144
4145       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4146       COMMON/MAGERR/IMAGERR
4147
4148       DIMENSION X(7),F(3)
4149
4150       XT = X(2)
4151       X(2) = X(1)
4152       X(1) = XT
4153
4154       X0 = X(1)/100.
4155       Y0 = X(2)/100.
4156       Z0 = X(3)/100.
4157
4158       CALL FREG1(Z0,X0,Y0,FZ0,FX0,FY0,IND)
4159 **    PRINT 3000,Z0,X0,Y0,FZ0,FX0,FY0,IND
4160
4161       IF(IND.EQ.0) GOTO 1
4162       CALL FREG2(Z0,X0,Y0,FZ0,FX0,FY0,IND)
4163       IMAGERR = 0
4164       IF(IND.EQ.2) THEN
4165         IMAGERR = 1
4166 **        print 1000
4167 **        PRINT 3000,Z0,X0,Y0,FZ0,FX0,FY0,IND
4168       ENDIF
4169  1000 format(1x,'Attention!!! The point is out of range!!!')
4170
4171  3000 FORMAT(1X,'Z=',D13.7,1X,'X=',D13.7,1X,'Y=',D13.7,1X,
4172      & 'BZ=',D13.7,1X,'BX=',D13.7,1X,'BY=',D13.7,1X,'IND=',I3)
4173
4174   1   F(1) = FX0*10.
4175       F(2) = FY0*10.
4176       F(3) = FZ0*10.
4177
4178 **      X(1) = X0*100.
4179 **      X(2) = Y0*100.
4180 **      X(3) = Z0*100.
4181
4182
4183       FT = F(2)
4184       F(2) = F(1)
4185       F(1) = FT
4186
4187       XT = X(2)
4188       X(2) = X(1)
4189       X(1) = XT
4190
4191       RETURN
4192       END
4193
4194
4195 **************************************************
4196       SUBROUTINE FREG1(Z0,X0,Y0,FZ0,FX0,FY0,IND)
4197 **************************************************
4198       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4199 **      REAL *4 BX,BY,BZ
4200       COMMON/DAT1/Z(81),X(81),Y(81,44),DX,DZ,LPX,LPY,LPZ
4201       COMMON/DAT2/BX(81,81,44),BY(81,81,44),BZ(81,81,44)
4202       COMMON/REG1/ZMAX,ZMIN,XMAX,XMIN
4203       COMMON/REG2/AY1,CY1,AY2,CY2
4204       KC=1
4205       LC=1
4206       MC=1
4207       IND=0
4208       IF(Z0.LT.ZMIN.OR.Z0.GT.ZMAX)GO TO 100
4209       IF(X0.LT.XMIN.OR.X0.GT.XMAX)GO TO 100
4210       YY1=AY1*Z0+CY1
4211       YY2=AY2*Z0+CY2
4212  2000 FORMAT(1X,'YY1=',D15.7,1X,'YY2=',D15.7,1X,'Y0=',D15.7)
4213 c      PRINT 2000,YY1,YY2,Y0
4214       IF(Y0.LT.YY1)GO TO 100
4215       IF(Y0.GT.YY2)GO TO 100
4216       CALL FIZ(Z0,Z,DZ,KC,K0,Z1,Z2,81)
4217       CALL FIZ(X0,X,DX,LC,L0,X1,X2,81)
4218       DY=(YY2-YY1)/DFLOAT(LPY-1)
4219       YY=(Y0-YY1)
4220       M0=(YY/DY)
4221
4222       IF(Y0.GE.(YY1+DFLOAT(M0)*DY).AND.Y0.LE.(YY1+DFLOAT(M0+1)*DY))
4223      &GO TO 700
4224       M0=M0+1
4225  700  CONTINUE
4226       Y2=(Y0-(YY1+DFLOAT(M0)*DY))/DY
4227       Y1=1.-Y2
4228 **    write(*,*) 'm0 Y1 Y2',m0,Y1,Y2
4229 **    print *,' k0=',k0,' l0=',l0,' m0=',m0
4230 **    print *,' z1=',z1,' z2=',z2
4231       FX1=Z1*BX(K0,L0,M0)+Z2*BX(K0+1,L0,M0)
4232       FX2=Z2*BX(K0+1,L0+1,M0)+Z1*BX(K0,L0+1,M0)
4233       FFX1=X1*FX1+X2*FX2
4234       GX1=Z1*BX(K0,L0,M0+1)+Z2*BX(K0+1,L0,M0+1)
4235       GX2=Z2*BX(K0+1,L0+1,M0+1)+Z1*BX(K0,L0+1,M0+1)
4236       GGX1=X1*GX1+X2*GX2
4237       FX0=Y1*FFX1+Y2*GGX1
4238       FX1=Z1*BY(K0,L0,M0)+Z2*BY(K0+1,L0,M0)
4239       FX2=Z2*BY(K0+1,L0+1,M0)+Z1*BY(K0,L0+1,M0)
4240       FFX1=X1*FX1+X2*FX2
4241       GX1=Z1*BY(K0,L0,M0+1)+Z2*BY(K0+1,L0,M0+1)
4242       GX2=Z2*BY(K0+1,L0+1,M0+1)+Z1*BY(K0,L0+1,M0+1)
4243       GGX1=X1*GX1+X2*GX2
4244       FY0=Y1*FFX1+Y2*GGX1
4245       FX1=Z1*BZ(K0,L0,M0)+Z2*BZ(K0+1,L0,M0)
4246       FX2=Z2*BZ(K0+1,L0+1,M0)+Z1*BZ(K0,L0+1,M0)
4247       FFX1=X1*FX1+X2*FX2
4248       GX1=Z1*BZ(K0,L0,M0+1)+Z2*BZ(K0+1,L0,M0+1)
4249       GX2=Z2*BZ(K0+1,L0+1,M0+1)+Z1*BZ(K0,L0+1,M0+1)
4250       GGX1=X1*GX1+X2*GX2
4251       FZ0=Y1*FFX1+Y2*GGX1
4252       RETURN
4253  100  CONTINUE
4254       IND=1
4255  1000 format(1x,'Attention!!! The point is out of range!!!')
4256 C      print 1000
4257       RETURN
4258       END
4259
4260 *************************************************
4261       SUBROUTINE FIZ(Z0,Z,DEL,KI,K0,Z1,Z2,NDZ)
4262 *************************************************
4263       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4264 ** CC      DIMENSION Z(NDZ)
4265       DIMENSION Z(10000)
4266       DDEL=Z0-Z(KI)
4267       KDEL=INT(DDEL/DEL)
4268       KJ=KI+KDEL
4269       K0 = NDZ - 1 ! CCCC
4270 *      if (k0.gt.81) print*,'ndz=',ndz
4271       IF (KJ.GT.NDZ) THEN ! CCC
4272          K0 = NDZ-1
4273          GO TO 100
4274       ENDIF
4275       DO 1 K=KJ,NDZ-1 ! CCC
4276       IF(Z0.LT.Z(K)) THEN
4277         K0 = K
4278         GO TO 100
4279       ENDIF
4280  1    CONTINUE
4281  100  CONTINUE
4282 *      print *,'K0=',K0,' Z0',z0, Z(K0), Z(K0+1),z2
4283       if (k0.gt.81) print*,'k0=',k0
4284       Z2=(Z0-Z(K0))/(Z(K0+1)-Z(K0))
4285       Z1=1.-Z2
4286 **      write(*,*) 'ko z1 z2', K0,Z1,Z2,' ki=',ki,' kj=',kj,' K=',K
4287 **      write(*,*)' NDZ Z(K0) Z(K0+1)',NDZ,Z(K0), Z(K0+1)
4288       RETURN
4289       END
4290
4291 ***************************************************
4292       SUBROUTINE FREG2(Z0,X0,Y0,FZ0,FX0,FY0,IND)
4293 **************************************************
4294       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4295 **      REAL *4 BXP,BYP,BZP
4296       COMMON/SDAT1/ZP(51),RAD(10),FI(33),DZP,DFI,DR,YY0,LPPZ,NR,NFI
4297       COMMON/SDAT2/BXP(51,10,33),BYP(51,10,33),BZP(51,10,33)
4298       COMMON/SDAT4/B(2,2,32)
4299       COMMON/REG3/ZPMAX,ZPMIN,RMAX,RMIN
4300 cc      COMMON/CONST/PI2,EPS
4301       KP=32+1
4302       LP=1
4303       MP=1
4304       YY0=0.3
4305       EPS=0.1D-6
4306       PI2=0.6283185E+01
4307       R0=DSQRT((X0-YY0)**2+Y0**2)
4308 c       write (*,*)'ro=',R0
4309       IF(Z0.LT.ZPMIN.OR.Z0.GT.ZPMAX)GO TO 100
4310       IF(R0.LT.RMIN.OR.R0.GT.RMAX)GO TO 100
4311       IF(R0.LE.DR)GO TO 3000
4312       CALL FIZ(Z0,ZP,DZP,KP,K0,Z1,Z2,51)
4313       CALL FIZ(R0,RAD,DR,LP,L0,X1,X2,10)
4314 **      print *,' r0=',r0,' rad=',rad,' dr=',dr,' lp=',lp,' l0=',l0,
4315 **     &     ' x1=',x1,' x2=',x2
4316       FI0=DACOS((X0-YY0)/R0)
4317       IF(Y0.LT.0.D+0)FI0=PI2-FI0
4318       CALL FIZ(FI0,FI,DFI,MP,M0,Y1,Y2,32)
4319 **      print *,' Apres FIZ',' k0=',k0,' l0=',l0,' m0=',m0
4320       FX1=Z1*BXP(K0,L0,M0)+Z2*BXP(K0+1,L0,M0)
4321       FX2=Z2*BXP(K0+1,L0+1,M0)+Z1*BXP(K0,L0+1,M0)
4322       FFX1=X1*FX1+X2*FX2
4323       GX1=Z1*BXP(K0,L0,M0+1)+Z2*BXP(K0+1,L0,M0+1)
4324       GX2=Z2*BXP(K0+1,L0+1,M0+1)+Z1*BXP(K0,L0+1,M0+1)
4325       GGX1=X1*GX1+X2*GX2
4326       FX0=Y1*FFX1+Y2*GGX1
4327       FX1=Z1*BYP(K0,L0,M0)+Z2*BYP(K0+1,L0,M0)
4328       FX2=Z2*BYP(K0+1,L0+1,M0)+Z1*BYP(K0,L0+1,M0)
4329       FFX1=X1*FX1+X2*FX2
4330       GX1=Z1*BYP(K0,L0,M0+1)+Z2*BYP(K0+1,L0,M0+1)
4331       GX2=Z2*BYP(K0+1,L0+1,M0+1)+Z1*BYP(K0,L0+1,M0+1)
4332       GGX1=X1*GX1+X2*GX2
4333       FY0=Y1*FFX1+Y2*GGX1
4334       FX1=Z1*BZP(K0,L0,M0)+Z2*BZP(K0+1,L0,M0)
4335 ** CCC      FX2=Z2*BZ(K0+1,L0+1,M0)+Z1*BZ(K0,L0+1,M0)
4336       FX2=Z2*BZP(K0+1,L0+1,M0)+Z1*BZP(K0,L0+1,M0)
4337       FFX1=X1*FX1+X2*FX2
4338       GX1=Z1*BZP(K0,L0,M0+1)+Z2*BZP(K0+1,L0,M0+1)
4339       GX2=Z2*BZP(K0+1,L0+1,M0+1)+Z1*BZP(K0,L0+1,M0+1)
4340       GGX1=X1*GX1+X2*GX2
4341       FZ0=Y1*FFX1+Y2*GGX1
4342       IND=0
4343       RETURN
4344  100  CONTINUE
4345       IND=2
4346  1000 format(1x,'Attention!!! The point is out of range!!!')
4347 C      print 1000
4348       RETURN
4349  3000 CONTINUE
4350       IF(R0.LT.EPS)GO TO 4000
4351       CALL FIZ(Z0,ZP,DZP,KP,K0,Z1,Z2,51)
4352       XX=X0-YY0
4353       FI0=DACOS(XX/R0)
4354       IF(Y0.LT.0.D+0)FI0=PI2-FI0
4355       CALL FIZ(FI0,FI,DFI,MP,M0,Y1,Y2,32)
4356       ALF2=B(1,1,M0)*XX+B(1,2,M0)*Y0
4357       ALF3=B(2,1,M0)*XX+B(2,2,M0)*Y0
4358       ALF1=1.-ALF2-ALF3
4359       FX1=ALF1*BXP(K0,1,1)+ALF2*BXP(K0,1,M0)+ALF3*BXP(K0,1,M0+1)
4360       FX2=ALF1*BXP(K0+1,1,1)+ALF2*BXP(K0+1,1,M0)+ALF3*BXP(K0+1,1,M0+1)
4361       FX0=Z1*FX1+Z2*FX2
4362       FX1=ALF1*BYP(K0,1,1)+ALF2*BYP(K0,1,M0)+ALF3*BYP(K0,1,M0+1)
4363       FX2=ALF1*BYP(K0+1,1,1)+ALF2*BYP(K0+1,1,M0)+ALF3*BYP(K0+1,1,M0+1)
4364       FY0=Z1*FX1+Z2*FX2
4365       FX1=ALF1*BZP(K0,1,1)+ALF2*BZP(K0,1,M0)+ALF3*BZP(K0,1,M0+1)
4366       FX2=ALF1*BZP(K0+1,1,1)+ALF2*BZP(K0+1,1,M0)+ALF3*BZP(K0+1,1,M0+1)
4367       FZ0=Z1*FX1+Z2*FX2
4368 c      write(*,*) 'R<Dr:B(1,1,m0) B(1,2,m0) B(2,1,m0) B(2,2,m0)',
4369 c     +B(1,1,M0),B(1,2,M0),B(2,1,M0),B(2,2,M0)
4370 c      write(*,*)'BX(K0,1,1) BX(K0,1,M0)','BX(K0,1,M0+1) BX(K0+1,1,1)'
4371 c     + ,'BX(K0+1,1,M0) BX(K0+1,1,M0+1)', BX(K0,1,1),BX(K0,1,M0) ,
4372 c     + BX(K0,1,M0+1),BX(K0+1,1,1),BX(K0+1,1,M0),BX(K0+1,1,M0+1)
4373 c      write(*,*)'By(K0,1,1) By(K0,1,M0)','By(K0,1,M0+1) By(K0+1,1,1)'
4374 c     + ,'By(K0+1,1,M0) By(K0+1,1,M0+1)', By(K0,1,1),By(K0,1,M0) ,
4375 c     + By(K0,1,M0+1),By(K0+1,1,1),By(K0+1,1,M0),By(K0+1,1,M0+1)
4376 ccc      write (*,*)'Bz(K0,1,1) Bz(K0,1,M0) Bz(K0,1,M0+1) Bz(K0+1,1,1)
4377   77  FORMAT(5x,E15.7,2x,E15.7)
4378 c      PRINT 70
4379 **   70 FORMAT(22hBz(K0,1,1) Bz(K0,1,M0))
4380 cc      PRINT 77 , BzP(K0,1,1),BzP(K0,1,M0)
4381 cc      PRINT 71
4382 **   71 FORMAT(26hBz(K0,1,M0+1) Bz(K0+1,1,1))
4383 cc     PRINT 77, BzP(K0,1,M0+1),BzP(K0+1,1,1)
4384 cc     PRINT 72
4385 **   72 FORMAT(29hBz(K0+1,1,M0) Bz(K0+1,1,M0+1))
4386 cc      PRINT 77 ,BzP(K0+1,1,M0),BzP(K0+1,1,M0+1)
4387 cc   77 FORMAT(5x,D15.7,5x,D15.7)
4388
4389
4390       IND=0
4391       RETURN
4392  4000 CONTINUE
4393       CALL FIZ(Z0,ZP,DZP,KP,K0,Z1,Z2,51)
4394       FX0=Z1*BXP(K0,1,1)+Z2*BXP(K0+1,1,1)
4395       FY0=Z1*BYP(K0,1,1)+Z2*BYP(K0+1,1,1)
4396       FZ0=Z1*BZP(K0,1,1)+Z2*BZP(K0+1,1,1)
4397 c      write(*,*) ' R<eps: Bx(k) Bx(k+1) By(k) By(k+1) Bz(k) Bz(k+1)',
4398 c     + BXP(K0,1,1),BXP(K0+1,1,1),BYP(K0,1,1),BYP(K0+1,1,1),BZP(K0,1,1),
4399 c     + BZP(K0+1,1,1)
4400       IND=0
4401       RETURN
4402       END
4403
4404 ***********************************************************************
4405
4406       SUBROUTINE RECO_GRKUTA (CHARGE,STEP,VECT,VOUT)
4407 C.
4408 C.    ******************************************************************
4409 C.    *                                                                *
4410 C.    *  Runge-Kutta method for tracking a particle through a magnetic *
4411 C.    *  field. Uses Nystroem algorithm (See Handbook Nat. Bur. of     *
4412 C.    *  Standards, procedure 25.5.20)                                 *
4413 C.    *                                                                *
4414 C.    *  Input parameters                                              *
4415 C.    *       CHARGE    Particle charge                                *
4416 C.    *       STEP      Step size                                      *
4417 C.    *       VECT      Initial co-ords,direction cosines,momentum     *
4418 C.    *  Output parameters                                             *
4419 C.    *       VOUT      Output co-ords,direction cosines,momentum      *
4420 C.    *  User routine called                                           *
4421 C.    *       CALL GUFLD(X,F)                                          *
4422 C.    *                                                                *
4423 C.    *    ==>Called by : <USER>, GUSWIM                               *
4424 C.    *       Authors    R.Brun, M.Hansroul  *********                 *
4425 C.    *                  V.Perevoztchikov (CUT STEP implementation)    *
4426 C.    *                                                                *
4427 C.    *                                                                *
4428 C.    ******************************************************************
4429 C.
4430       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4431
4432 **      REAL CHARGE, STEP, VECT(*), VOUT(*), F(4)
4433 **      REAL XYZT(3), XYZ(3), X, Y, Z, XT, YT, ZT
4434       DIMENSION VECT(*), VOUT(*), F(3)
4435       DIMENSION XYZT(3), XYZ(3)
4436       DIMENSION SECXS(4),SECYS(4),SECZS(4),HXP(3)
4437       EQUIVALENCE (X,XYZ(1)),(Y,XYZ(2)),(Z,XYZ(3)),
4438      +            (XT,XYZT(1)),(YT,XYZT(2)),(ZT,XYZT(3))
4439 *
4440       PARAMETER (MAXIT = 1992, MAXCUT = 11)
4441       PARAMETER (EC=2.9979251D-4,DLT=1D-4,DLT32=DLT/32)
4442       PARAMETER (ZERO=0, ONE=1, TWO=2, THREE=3)
4443       PARAMETER (THIRD=ONE/THREE, HALF=ONE/TWO)
4444       PARAMETER (PISQUA=.986960440109D+01)
4445       PARAMETER      (IX=1,IY=2,IZ=3,IPX=4,IPY=5,IPZ=6)
4446 *.
4447 *.    ------------------------------------------------------------------
4448 *.
4449 *             This constant is for units CM,GEV/C and KGAUSS
4450 *
4451       ITER = 0
4452       NCUT = 0
4453       DO 10 J=1,7
4454          VOUT(J)=VECT(J)
4455    10 CONTINUE
4456       PINV   = EC * CHARGE / VECT(7)
4457       TL = 0.
4458       H      = STEP
4459 *
4460 *
4461    20 REST  = STEP-TL
4462       IF (ABS(H).GT.ABS(REST)) H = REST
4463       CALL RECO_GUFLD(VOUT,F)
4464 *
4465 *             Start of integration
4466 *
4467       X      = VOUT(1)
4468       Y      = VOUT(2)
4469       Z      = VOUT(3)
4470       A      = VOUT(4)
4471       B      = VOUT(5)
4472       C      = VOUT(6)
4473 *
4474       H2     = HALF * H
4475       H4     = HALF * H2
4476       PH     = PINV * H
4477       PH2    = HALF * PH
4478       SECXS(1) = (B * F(3) - C * F(2)) * PH2
4479       SECYS(1) = (C * F(1) - A * F(3)) * PH2
4480       SECZS(1) = (A * F(2) - B * F(1)) * PH2
4481       ANG2 = (SECXS(1)**2 + SECYS(1)**2 + SECZS(1)**2)
4482       IF (ANG2.GT.PISQUA) GO TO 40
4483       DXT    = H2 * A + H4 * SECXS(1)
4484       DYT    = H2 * B + H4 * SECYS(1)
4485       DZT    = H2 * C + H4 * SECZS(1)
4486       XT     = X + DXT
4487       YT     = Y + DYT
4488       ZT     = Z + DZT
4489 *
4490 *              Second intermediate point
4491 *
4492       EST = ABS(DXT)+ABS(DYT)+ABS(DZT)
4493       IF (EST.GT.H) GO TO 30
4494
4495       CALL RECO_GUFLD(XYZT,F)
4496       AT     = A + SECXS(1)
4497       BT     = B + SECYS(1)
4498       CT     = C + SECZS(1)
4499 *
4500       SECXS(2) = (BT * F(3) - CT * F(2)) * PH2
4501       SECYS(2) = (CT * F(1) - AT * F(3)) * PH2
4502       SECZS(2) = (AT * F(2) - BT * F(1)) * PH2
4503       AT     = A + SECXS(2)
4504       BT     = B + SECYS(2)
4505       CT     = C + SECZS(2)
4506       SECXS(3) = (BT * F(3) - CT * F(2)) * PH2
4507       SECYS(3) = (CT * F(1) - AT * F(3)) * PH2
4508       SECZS(3) = (AT * F(2) - BT * F(1)) * PH2
4509       DXT    = H * (A + SECXS(3))
4510       DYT    = H * (B + SECYS(3))
4511       DZT    = H * (C + SECZS(3))
4512       XT     = X + DXT
4513       YT     = Y + DYT
4514       ZT     = Z + DZT
4515       AT     = A + TWO*SECXS(3)
4516       BT     = B + TWO*SECYS(3)
4517       CT     = C + TWO*SECZS(3)
4518 *
4519       EST = ABS(DXT)+ABS(DYT)+ABS(DZT)
4520       IF (EST.GT.2.*ABS(H)) GO TO 30
4521
4522       CALL RECO_GUFLD(XYZT,F)
4523 *
4524       Z      = Z + (C + (SECZS(1) + SECZS(2) + SECZS(3)) * THIRD) * H
4525       Y      = Y + (B + (SECYS(1) + SECYS(2) + SECYS(3)) * THIRD) * H
4526       X      = X + (A + (SECXS(1) + SECXS(2) + SECXS(3)) * THIRD) * H
4527 *
4528       SECXS(4) = (BT*F(3) - CT*F(2))* PH2
4529       SECYS(4) = (CT*F(1) - AT*F(3))* PH2
4530       SECZS(4) = (AT*F(2) - BT*F(1))* PH2
4531       A      = A+(SECXS(1)+SECXS(4)+TWO * (SECXS(2)+SECXS(3))) * THIRD
4532       B      = B+(SECYS(1)+SECYS(4)+TWO * (SECYS(2)+SECYS(3))) * THIRD
4533       C      = C+(SECZS(1)+SECZS(4)+TWO * (SECZS(2)+SECZS(3))) * THIRD
4534 *
4535       EST    = ABS(SECXS(1)+SECXS(4) - (SECXS(2)+SECXS(3)))
4536      ++        ABS(SECYS(1)+SECYS(4) - (SECYS(2)+SECYS(3)))
4537      ++        ABS(SECZS(1)+SECZS(4) - (SECZS(2)+SECZS(3)))
4538 *
4539       IF (EST.GT.DLT .AND. ABS(H).GT.1.E-4) GO TO 30
4540       ITER = ITER + 1
4541       NCUT = 0
4542 *               If too many iterations, go to HELIX
4543       IF (ITER.GT.MAXIT) GO TO 40
4544 *
4545       TL = TL + H
4546       IF (EST.LT.(DLT32)) THEN
4547          H = H*TWO
4548       ENDIF
4549       CBA    = ONE/ SQRT(A*A + B*B + C*C)
4550       VOUT(1) = X
4551       VOUT(2) = Y
4552       VOUT(3) = Z
4553       VOUT(4) = CBA*A
4554       VOUT(5) = CBA*B
4555       VOUT(6) = CBA*C
4556       REST = STEP - TL
4557       IF (STEP.LT.0.) REST = -REST
4558       IF (REST .GT. 1.E-5*ABS(STEP)) GO TO 20
4559 *
4560       GO TO 999
4561 *
4562 **              CUT STEP
4563    30 NCUT = NCUT + 1
4564 *               If too many cuts , go to HELIX
4565       IF (NCUT.GT.MAXCUT)       GO TO 40
4566       H = H*HALF
4567       GO TO 20
4568 *
4569 **              ANGLE TOO BIG, USE HELIX
4570    40 F1  = F(1)
4571       F2  = F(2)
4572       F3  = F(3)
4573       F4  = SQRT(F1**2+F2**2+F3**2)
4574       RHO = -F4*PINV
4575       TET = RHO * STEP
4576       IF(TET.NE.0.) THEN
4577          HNORM = ONE/F4
4578          F1 = F1*HNORM
4579          F2 = F2*HNORM
4580          F3 = F3*HNORM
4581 *
4582          HXP(1) = F2*VECT(IPZ) - F3*VECT(IPY)
4583          HXP(2) = F3*VECT(IPX) - F1*VECT(IPZ)
4584          HXP(3) = F1*VECT(IPY) - F2*VECT(IPX)
4585
4586          HP = F1*VECT(IPX) + F2*VECT(IPY) + F3*VECT(IPZ)
4587 *
4588          RHO1 = ONE/RHO
4589          SINT = SIN(TET)
4590          COST = TWO*SIN(HALF*TET)**2
4591 *
4592          G1 = SINT*RHO1
4593          G2 = COST*RHO1
4594          G3 = (TET-SINT) * HP*RHO1
4595          G4 = -COST
4596          G5 = SINT
4597          G6 = COST * HP
4598
4599          VOUT(IX) = VECT(IX) + (G1*VECT(IPX) + G2*HXP(1) + G3*F1)
4600          VOUT(IY) = VECT(IY) + (G1*VECT(IPY) + G2*HXP(2) + G3*F2)
4601          VOUT(IZ) = VECT(IZ) + (G1*VECT(IPZ) + G2*HXP(3) + G3*F3)
4602
4603          VOUT(IPX) = VECT(IPX) + (G4*VECT(IPX) + G5*HXP(1) + G6*F1)
4604          VOUT(IPY) = VECT(IPY) + (G4*VECT(IPY) + G5*HXP(2) + G6*F2)
4605          VOUT(IPZ) = VECT(IPZ) + (G4*VECT(IPZ) + G5*HXP(3) + G6*F3)
4606 *
4607       ELSE
4608          VOUT(IX) = VECT(IX) + STEP*VECT(IPX)
4609          VOUT(IY) = VECT(IY) + STEP*VECT(IPY)
4610          VOUT(IZ) = VECT(IZ) + STEP*VECT(IPZ)
4611 *
4612       ENDIF
4613 *
4614   999 END
4615
4616
4617 *******************************************************************
4618       SUBROUTINE RECO_GUFLDOLD1(X,B)
4619 C
4620 C    CONSTANT FIELD
4621 C
4622 C *** ROUTINE DESCRIBING THE MAGNETIC FIELD IN THE ALICE SETUP ***
4623 C *** NVE 14-NOV-1990 CERN GENEVA ***
4624 C
4625 C CALLED BY : GUFLD
4626 C ORIGIN    : NICK VAN EIJNDHOVEN
4627 C
4628 C Input :
4629 C -------
4630 C X = (X,Y,Z) coordinates in cm
4631 C
4632 C Output :
4633 C --------
4634 C B    = Magnetic field components (BX,BY,BZ) in KG
4635 C
4636
4637       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4638
4639       PARAMETER(NBSTATION=5)
4640 C
4641       COMMON/ZDEFIN/ZPLANE(NBSTATION),ZCOIL,ZMAGEND,DZ_PL(NBSTATION)
4642
4643
4644 C
4645       DIMENSION X(3),B(3)
4646
4647       XT = X(2)
4648       X(2) = X(1)
4649       X(1) = XT
4650
4651       B(1) = 0.
4652       B(2) = 0.
4653       B(3) = 0.
4654
4655       IF (X(3).LT.(-1.*ZCOIL)) THEN
4656          B(3) = 2.
4657       ELSEIF ( X(3).LT.(-1.*ZMAGEND)) THEN
4658          B(1) = -10.
4659       ENDIF
4660
4661       BT = B(2)
4662       B(2) = B(1)
4663       B(1) = BT
4664
4665       XT = X(2)
4666       X(2) = X(1)
4667       X(1) = XT
4668 **      print *,' x =',X(1),X(2),X(3)
4669 **      print *,' B =',B(1),B(2),B(3)
4670
4671 C
4672
4673
4674  999  END
4675
4676 *******************************************************************
4677       SUBROUTINE RECO_GUFLD(X,B)
4678 C
4679 C    Field map Mariana
4680 C
4681 C *** ROUTINE DESCRIBING THE MAGNETIC FIELD IN THE ALICE SETUP ***
4682 C *** NVE 14-NOV-1990 CERN GENEVA ***
4683 C
4684 C CALLED BY : GUFLD
4685 C ORIGIN    : NICK VAN EIJNDHOVEN
4686 C
4687 C Input :
4688 C -------
4689 C X = (X,Y,Z) coordinates in cm
4690 C
4691 C Output :
4692 C --------
4693 C B    = Magnetic field components (BX,BY,BZ) in KG
4694 C
4695
4696       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4697 C
4698
4699 C --- Common containing magnetic field map data
4700       REAL DZ,DX,DY,UDX,UDY,UDZ
4701      $,XMBEG,YMBEG,ZMBEG,XMEND,YMEND,ZMEND
4702      $,BV
4703       INTEGER NX,NY,NZ
4704
4705       PARAMETER(MAXFLD=250000)
4706       COMMON /SCXMFD/ NX,NY,NZ,DZ,DX,DY,UDX,UDY,UDZ
4707      $,XMBEG,YMBEG,ZMBEG,XMEND,YMEND,ZMEND
4708      $,BV(MAXFLD)
4709 C
4710
4711 C
4712       DIMENSION X(3),B(3)
4713       DOUBLE PRECISION ONE, RATX, RATY, RATZ, HIX, HIY, HIZ
4714      $,  RATX1, RATY1, RATZ1
4715      $,  BHYHZ, BHYLZ, BLYHZ, BLYLZ, BHZ, BLZ
4716      $,  XL(3)
4717
4718       PARAMETER (ONE=1)
4719
4720       EXTERNAL BX,BY,BZ
4721
4722       XT = X(2)
4723       X(2) = X(1)
4724       X(1) = XT
4725
4726
4727       ISXFMAP = 3
4728
4729 **      BX(JX,JY,JZ)=BV(3*((JZ-1)*(NX*NY)+(JY-1)*NX+(JX-1))+1)
4730 **      BY(JX,JY,JZ)=BV(3*((JZ-1)*(NX*NY)+(JY-1)*NX+(JX-1))+2)
4731 **      BZ(JX,JY,JZ)=BV(3*((JZ-1)*(NX*NY)+(JY-1)*NX+(JX-1))+3)
4732
4733
4734
4735 *
4736 * --- Act accordingly to ISXFMAP
4737 *
4738       IF(ISXFMAP.EQ.1) THEN
4739          IF (ABS(X(3)) .LT. 700.
4740      +      .AND. X(1)**2+(X(2)+30.)**2 .LT. 560.**2 ) THEN
4741             B(1)=0.
4742             B(2)=0.
4743             B(3)=2.
4744          ELSE IF (X(3) .GE. 725. .AND. X(3) .LT. 1225.) THEN
4745             DZ=ABS(975.-X(3))/100.
4746             B(1)=(1.-0.1*DZ*DZ)*7.0
4747             B(2)=0.
4748             B(3)=0.
4749          ELSE
4750             B(1)=0.
4751             B(2)=0.
4752             B(3)=0.
4753           ENDIF
4754       ELSE IF(ISXFMAP.LE.3) THEN
4755          IF (-700.LT.X(3).AND.X(3).LT.ZMBEG
4756      +      .AND. X(1)**2+(X(2)+30.)**2 .LT. 560.**2 ) THEN
4757             B(1)=0.
4758             B(2)=0.
4759             B(3)=2.
4760          ELSE IF ((X(3) .GE. ZMBEG .AND. X(3) .LT. ZMEND) .AND.
4761      +            (XMBEG.LE.ABS(X(1)).AND.ABS(X(1)).LT.XMEND) .AND.
4762      +            (YMBEG.LE.ABS(X(2)).AND.ABS(X(2)).LT.YMEND)) THEN
4763
4764
4765
4766 C --- find the position in the grid ---
4767
4768            XL(1)=ABS(X(1))-XMBEG
4769            XL(2)=ABS(X(2))-YMBEG
4770            XL(3)=X(3)-ZMBEG
4771
4772 C --- Start with X
4773
4774            HIX=XL(1)*UDX
4775            RATX=HIX-AINT(HIX)
4776            IX=HIX+1
4777
4778            HIY=XL(2)*UDY
4779            RATY=HIY-AINT(HIY)
4780            IY=HIY+1
4781
4782            HIZ=XL(3)*UDZ
4783            RATZ=HIZ-AINT(HIZ)
4784            IZ=HIZ+1
4785
4786            IF(ISXFMAP.EQ.2) THEN
4787 * ... Simple interpolation
4788
4789               B(1) = BX(IX,IY,IZ)*(ONE-RATX) + BX(IX+1,IY+1,IZ+1)*RATX
4790               B(2) = BY(IX,IY,IZ)*(ONE-RATY) + BY(IX+1,IY+1,IZ+1)*RATY
4791               B(3) = BZ(IX,IY,IZ)*(ONE-RATZ) + BZ(IX+1,IY+1,IZ+1)*RATZ
4792            ELSE IF(ISXFMAP.EQ.3) THEN
4793 * ... more complicated interpolation
4794               RATX1=ONE-RATX
4795               RATY1=ONE-RATY
4796               RATZ1=ONE-RATZ
4797 **              print *,' bx by bz', BX(IX  ,IY+1,IZ+1),BY(IX  ,IY+1,IZ+1)
4798 **     &        , BZ(IX  ,IY+1,IZ+1)
4799               BHYHZ = BX(IX  ,IY+1,IZ+1)*RATX1+BX(IX+1,IY+1,IZ+1)*RATX
4800               BHYLZ = BX(IX  ,IY+1,IZ  )*RATX1+BX(IX+1,IY+1,IZ  )*RATX
4801               BLYHZ = BX(IX  ,IY  ,IZ+1)*RATX1+BX(IX+1,IY  ,IZ+1)*RATX
4802               BLYLZ = BX(IX  ,IY  ,IZ  )*RATX1+BX(IX+1,IY  ,IZ  )*RATX
4803               BHZ   = BLYHZ             *RATY1+BHYHZ             *RATY
4804               BLZ   = BLYLZ             *RATY1+BHYLZ             *RATY
4805               B(1)  = BLZ               *RATZ1+BHZ               *RATZ
4806 *
4807               BHYHZ = BY(IX  ,IY+1,IZ+1)*RATX1+BY(IX+1,IY+1,IZ+1)*RATX
4808               BHYLZ = BY(IX  ,IY+1,IZ  )*RATX1+BY(IX+1,IY+1,IZ  )*RATX
4809               BLYHZ = BY(IX  ,IY  ,IZ+1)*RATX1+BY(IX+1,IY  ,IZ+1)*RATX
4810               BLYLZ = BY(IX  ,IY  ,IZ  )*RATX1+BY(IX+1,IY  ,IZ  )*RATX
4811               BHZ   = BLYHZ             *RATY1+BHYHZ             *RATY
4812               BLZ   = BLYLZ             *RATY1+BHYLZ             *RATY
4813               B(2)  = BLZ               *RATZ1+BHZ               *RATZ
4814 *
4815               BHYHZ = BZ(IX  ,IY+1,IZ+1)*RATX1+BZ(IX+1,IY+1,IZ+1)*RATX
4816               BHYLZ = BZ(IX  ,IY+1,IZ  )*RATX1+BZ(IX+1,IY+1,IZ  )*RATX
4817               BLYHZ = BZ(IX  ,IY  ,IZ+1)*RATX1+BZ(IX+1,IY  ,IZ+1)*RATX
4818               BLYLZ = BZ(IX  ,IY  ,IZ  )*RATX1+BZ(IX+1,IY  ,IZ  )*RATX
4819               BHZ   = BLYHZ             *RATY1+BHYHZ             *RATY
4820               BLZ   = BLYLZ             *RATY1+BHYLZ             *RATY
4821               B(3)  = BLZ               *RATZ1+BHZ               *RATZ
4822 *
4823            ENDIF
4824 * ... use the dipole symmetry
4825
4826            IF (X(1)*X(2).LT.0) B(2)=-B(2)
4827            IF (X(1).LT.0) B(3)=-B(3)
4828
4829 *
4830
4831       ELSE
4832
4833            B(1)=0.
4834            B(2)=0.
4835            B(3)=0.
4836
4837       ENDIF ! z-coord for m.f.
4838
4839       ENDIF ! endif ISXFMAP
4840
4841       SXMAGN = 1.
4842       SXMGMX = 20.
4843  12   B(1)=B(1)*SXMAGN
4844       B(2)=B(2)*SXMAGN
4845       B(3)=B(3)*SXMAGN
4846       BTOT=SQRT(B(1)**2+B(2)**2+B(3)**2)
4847       IF(BTOT.GT.SXMGMX) THEN
4848          PRINT 10100, BTOT,SXMGMX
4849 10100    FORMAT(' *GUFLD* Field ',G10.4,' larger than max ',G10.4)
4850       ENDIF
4851
4852
4853       BT = B(2)
4854       B(2) = B(1)
4855       B(1) = BT
4856
4857       XT = X(2)
4858       X(2) = X(1)
4859       X(1) = XT
4860
4861 C
4862
4863       RETURN
4864       END
4865
4866
4867 *******************************************
4868       DOUBLE PRECISION FUNCTION BX(JX,JY,JZ)
4869
4870       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4871
4872 C --- Common containing magnetic field map data
4873       REAL DZ,DX,DY,UDX,UDY,UDZ
4874      $,XMBEG,YMBEG,ZMBEG,XMEND,YMEND,ZMEND
4875      $,BV
4876       INTEGER NX,NY,NZ
4877
4878       PARAMETER(MAXFLD=250000)
4879       COMMON /SCXMFD/ NX,NY,NZ,DZ,DX,DY,UDX,UDY,UDZ
4880      $,XMBEG,YMBEG,ZMBEG,XMEND,YMEND,ZMEND
4881      $,BV(MAXFLD)
4882
4883       BX=BV(3*((JZ-1)*(NX*NY)+(JY-1)*NX+(JX-1))+1)
4884
4885       END
4886
4887 *******************************************
4888        DOUBLE PRECISION FUNCTION BY(JX,JY,JZ)
4889
4890       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4891
4892 C --- Common containing magnetic field map data
4893       REAL DZ,DX,DY,UDX,UDY,UDZ
4894      $,XMBEG,YMBEG,ZMBEG,XMEND,YMEND,ZMEND
4895      $,BV
4896       INTEGER NX,NY,NZ
4897
4898       PARAMETER(MAXFLD=250000)
4899       COMMON /SCXMFD/ NX,NY,NZ,DZ,DX,DY,UDX,UDY,UDZ
4900      $,XMBEG,YMBEG,ZMBEG,XMEND,YMEND,ZMEND
4901      $,BV(MAXFLD)
4902
4903       BY=BV(3*((JZ-1)*(NX*NY)+(JY-1)*NX+(JX-1))+2)
4904
4905       END
4906
4907 *******************************************
4908       DOUBLE PRECISION FUNCTION BZ(JX,JY,JZ)
4909
4910       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4911
4912 C --- Common containing magnetic field map data
4913       REAL DZ,DX,DY,UDX,UDY,UDZ
4914      $,XMBEG,YMBEG,ZMBEG,XMEND,YMEND,ZMEND
4915      $,BV
4916       INTEGER NX,NY,NZ
4917
4918       PARAMETER(MAXFLD=250000)
4919       COMMON /SCXMFD/ NX,NY,NZ,DZ,DX,DY,UDX,UDY,UDZ
4920      $,XMBEG,YMBEG,ZMBEG,XMEND,YMEND,ZMEND
4921      $,BV(MAXFLD)
4922
4923       BZ=BV(3*((JZ-1)*(NX*NY)+(JY-1)*NX+(JX-1))+3)
4924
4925       END
4926
4927 ***********************************************************************
4928       SUBROUTINE INITFIELD
4929
4930 C
4931 C    Marianna
4932 C
4933 C *** INITIALISATION OF THE FIELD MAP ***
4934 C *** FCA 24-AUG-1998 CERN GENEVA ***
4935 C
4936 C CALLED BY : GALICE
4937 C ORIGIN    : NICK VAN EIJNDHOVEN
4938 C
4939       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4940
4941
4942 C --- Common containing magnetic field map data
4943       REAL DZ,DX,DY,UDX,UDY,UDZ
4944      $,XMBEG,YMBEG,ZMBEG,XMEND,YMEND,ZMEND
4945      $,BV
4946       INTEGER NX,NY,NZ
4947
4948       PARAMETER(MAXFLD=250000)
4949       COMMON /SCXMFD/ NX,NY,NZ,DZ,DX,DY,UDX,UDY,UDZ
4950      $,XMBEG,YMBEG,ZMBEG,XMEND,YMEND,ZMEND
4951      $,BV(MAXFLD)
4952
4953       ISXFMAP = 3
4954
4955       IF(ISXFMAP.EQ.1) THEN
4956 * ... constant field, nothing to do
4957       ELSE IF(ISXFMAP.LE.3) THEN
4958 * ... constant mesh field
4959          PRINT 10000, ISXFMAP
4960 10000    FORMAT(' *SXFMAP* Magnetic field map flag: ',I5
4961      $        ,'; Reading magnetic field map data ')
4962 *
4963          OPEN(77,FILE='data/field01.dat',
4964      $        FORM='FORMATTED',STATUS='OLD')
4965          READ(77,*) NX,NY,NZ,DX,DY,DZ,XMBEG,YMBEG,ZMBEG
4966          PRINT*,'NX,NY,NZ,DX,DY,DZ,XMBEG,YMBEG,ZMBEG',
4967      $   NX,NY,NZ,DX,DY,DZ,XMBEG,YMBEG,ZMBEG
4968          IF(3*NX*NY*NZ.GT.MAXFLD) THEN
4969             WRITE(6,10100) 3*NX*NY*NZ,MAXFLD
4970             STOP 'Increase MAXFLD'
4971          ENDIF
4972          UDX=1/DX
4973          UDY=1/DY
4974          UDZ=1/DZ
4975          XMEND=XMBEG+(NX-1)*DX
4976          YMEND=YMBEG+(NY-1)*DY
4977          ZMEND=ZMBEG+(NZ-1)*DZ
4978          DO IZ=1,NZ
4979             IPZ=3*(IZ-1)*(NX*NY)
4980             DO IY=1,NY
4981                IPY=IPZ+3*(IY-1)*NX
4982                DO IX=1,NX
4983                   IPX=IPY+3*(IX-1)
4984                   READ(77,*) BV(IPX+3),BV(IPX+2),BV(IPX+1)
4985                ENDDO
4986             ENDDO
4987          ENDDO
4988          CLOSE(77)
4989       ENDIF                     ! endif ISXFMAP
4990 *
4991 10100 FORMAT('*** SXFMAP: Need ',I7,' have ',I7)
4992       END
4993
4994 ****************************************
4995       SUBROUTINE RANNOR (A,B)
4996 ****************************************
4997 C
4998 C CERN PROGLIB# V100    RANNOR          .VERSION KERNFOR  4.18  880425
4999 C ORIG. 18/10/77
5000 C
5001       Y = RNDM()
5002       IF (Y.EQ.0.)  Y = RNDM()
5003       Z = RNDM()
5004
5005       X = 6.283185*Z
5006       R = SQRT (-2.0*LOG(Y))
5007       A = R*SIN (X)
5008       B = R*COS (X)
5009       RETURN
5010       END
5011 ************************************************************************
5012       SUBROUTINE OLDFCNFIT(NPAR,GRAD,FVAL,XVAL,IFLAG,FUTIL)
5013 ************************************************************************
5014 *    Calcule FVAL: la fonction minimisee par MINUIT
5015 *    with constant magnetic Field
5016 *
5017 ************************************************************************
5018
5019       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
5020
5021       PARAMETER(NPLANE=10)
5022
5023       COMMON/PARAM/ZPLANEP(NPLANE),THICK,XPREC,YPREC,B0,BL3,ZMAGS,
5024      &             ZMAGE,ZABS,XMAG,ZBP1,ZBP2,CONST
5025
5026       COMMON/MEAS/LPLANE(NPLANE),XMP(NPLANE),YMP(NPLANE)
5027
5028       COMMON/FCNOUT/PXZEA,ALAMEA,PHIEA,XEA,YEA,NPLU,CHI2
5029
5030       DIMENSION GRAD(*),XVAL(*),AMS(500),DISTAZ(500)
5031
5032       DIMENSION XP(NPLANE),YP(NPLANE),
5033      &          COV(NPLANE,NPLANE),AP(NPLANE),COVY(NPLANE,NPLANE)
5034       DIMENSION VECT(7),VOUT(7)
5035
5036       XV = XVAL(4)
5037       YV = XVAL(5)
5038
5039       ASIGN = 1.
5040       IF (XVAL(1).LT.0.) ASIGN = -1.
5041       PHI  = XVAL(2)
5042       ALAM = XVAL(3)
5043       PXZ = DABS(1./XVAL(1))
5044
5045       PX = PXZ*DSIN(PHI)
5046       PY = PXZ*DTAN(ALAM)
5047       PZ = PXZ*DCOS(PHI)
5048       PTOT = PXZ/DCOS(ALAM)
5049       TTHET = DSQRT(DTAN(ALAM)**2+DSIN(PHI)**2)/DCOS(PHI)
5050       THET = DATAN(TTHET)
5051       PT = DSQRT(PX**2+PY**2)
5052
5053       RL3 = ASIGN*PT / (0.299792458D-3 * BL3)
5054       ALPHA = DATAN2(PY,PX)
5055       XC = XV+RL3*DSIN(ALPHA)
5056       YC = YV-RL3*DCOS(ALPHA)
5057
5058       A12 = 0.
5059       NTMAX = 0
5060
5061       ZEA = ZABS
5062       XEA = XV
5063       YEA = YV
5064       PXEA = PX
5065       PYEA = PY
5066       PHIEA = PHI
5067       PXZEA = ASIGN*PXZ
5068       ALAMEA = ALAM
5069 * 1er plan
5070       ANGDEV = (ZPLANEP(1)-ZEA)*TTHET/RL3
5071       XP(1) = XC+RL3*DSIN(ANGDEV-ALPHA)
5072       YP(1) = YC+RL3*DCOS(ANGDEV-ALPHA)
5073       AL   = THICK/ DCOS(THET)
5074       AP(1)  = (0.0136D0/PTOT) * DSQRT(AL) * (1 + 0.038D0*DLOG(AL))
5075 * 2eme plan
5076       ANGDEV = (ZPLANEP(2)-ZEA)*TTHET/RL3
5077       XP(2) = XC+RL3*DSIN(ANGDEV-ALPHA)
5078       YP(2) = YC+RL3*DCOS(ANGDEV-ALPHA)
5079       AL   = THICK/ DCOS(THET)
5080       AP(2) = (0.0136D0/PTOT) * DSQRT(AL) * (1 + 0.038D0*DLOG(AL))
5081 * 3eme plan
5082       ANGDEV = (ZPLANEP(3)-ZEA)*TTHET/RL3
5083       XP(3) = XC+RL3*DSIN(ANGDEV-ALPHA)
5084       YP(3) = YC+RL3*DCOS(ANGDEV-ALPHA)
5085       AL   = THICK/ DCOS(THET)
5086       AP(3)  = (0.0136D0/PTOT) * DSQRT(AL) * (1 + 0.038D0*DLOG(AL))
5087 * 4eme plan
5088       ANGDEV = (ZPLANEP(4)-ZEA)*TTHET/RL3
5089       XP(4) = XC+RL3*DSIN(ANGDEV-ALPHA)
5090       YP(4) = YC+RL3*DCOS(ANGDEV-ALPHA)
5091       AL   = THICK/ DCOS(THET)
5092       AP(4) = (0.0136D0/PTOT) * DSQRT(AL) * (1 + 0.038D0*DLOG(AL))
5093 * Fin de L3
5094       ANGDEV = (700.D0-ZEA)*TTHET/RL3
5095       XPL3 = XC+RL3*DSIN(ANGDEV-ALPHA)
5096       YPL3 = YC+RL3*DCOS(ANGDEV-ALPHA)
5097       PX = PT*DCOS(ANGDEV-ALPHA)
5098       PY = -PT*DSIN(ANGDEV-ALPHA)
5099       PHIC = DATAN2(PY,PX)
5100       CX = DSIN(THET)*DCOS(PHIC)
5101       CY = DSIN(THET)*DSIN(PHIC)
5102       CZ = DCOS(THET)
5103 * Entree du dipole
5104       VECT(1) = XPL3+(ZMAGS-700.)*CX/CZ
5105       VECT(2) = YPL3+(ZMAGS-700.)*CY/CZ
5106       VECT(3) = ZMAGS
5107       VECT(4) = CX
5108       VECT(5) = CY
5109       VECT(6) = CZ
5110       VECT(7) = PTOT
5111
5112       PXZ  = PTOT*DSQRT(VECT(4)**2+VECT(6)**2)
5113       RDIP = ASIGN*PXZ / (0.299792458D-3 * B0)
5114       PHI1 = DATAN2(VECT(4),VECT(6))
5115       XC   = VECT(1)+RDIP*DCOS(PHI1)
5116       ZC   = VECT(3)-RDIP*DSIN(PHI1)
5117
5118       IF (DABS((ZMAGE-ZPLANEP(5))/RDIP).GE.1.D0) RETURN ! Particule boucle dans le champ
5119       XP(5) = XC-ASIGN*DSQRT(RDIP**2-(ZPLANEP(5)-ZC)**2)
5120       YP(5) = VECT(2)+(ZPLANEP(5)-ZMAGS)*VECT(5)/VECT(6)
5121       CX = (ZPLANEP(5)-ZC)/RDIP
5122       CY = VECT(5)
5123       CZ = DABS((XP(5)-XC)/RDIP)
5124       THET = DATAN2(DSQRT(CX**2+CY**2),CZ)
5125       AL   = THICK/ DCOS(THET)
5126       AP(5) = (0.0136D0/PTOT) * DSQRT(AL) * (1 + 0.038D0*DLOG(AL))
5127
5128       IF (DABS((ZMAGE-ZPLANEP(6))/RDIP).GE.1.D0) RETURN ! Particule boucle dans le champ
5129       XP(6) = XC-ASIGN*DSQRT(RDIP**2-(ZPLANEP(6)-ZC)**2)
5130       YP(6) = VECT(2)+(ZPLANEP(6)-ZMAGS)*VECT(5)/VECT(6)
5131       CX = (ZPLANEP(6)-ZC)/RDIP
5132       CY = VECT(5)
5133       CZ = DABS((XP(6)-XC)/RDIP)
5134       THET = DATAN2(DSQRT(CX**2+CY**2),CZ)
5135       AL   = THICK/ DCOS(THET)
5136       AP(6) = (0.0136D0/PTOT) * DSQRT(AL) * (1 + 0.038D0*DLOG(AL))
5137
5138
5139       IF (DABS((ZMAGE-ZC)/RDIP).GE.1.D0) RETURN ! Particule boucle dans le champ
5140       VOUT(1) = XC-ASIGN*DSQRT(RDIP**2-(ZMAGE-ZC)**2)
5141       VOUT(2) = VECT(2)+(ZMAGE-ZMAGS)*VECT(5)/VECT(6)
5142       VOUT(3) = ZMAGE
5143       VOUT(4) = (ZMAGE-ZC)/RDIP
5144       VOUT(5) = VECT(5)
5145       VOUT(6) = DABS((VOUT(1)-XC)/RDIP)
5146       VOUT(7) = PTOT
5147
5148       DO IV = 1,7
5149          VECT(IV) = VOUT(IV)
5150       ENDDO
5151
5152 * 7eme plan
5153       THET = DATAN2(DSQRT(VECT(4)**2+VECT(5)**2),VECT(6))
5154       XP(7) = VECT(1)+(ZPLANEP(7)-ZMAGE)*VECT(4)/VECT(6)
5155       YP(7) = VECT(2)+(ZPLANEP(7)-ZMAGE)*VECT(5)/VECT(6)
5156       AL   = THICK/ DCOS(THET)
5157       AP(7)  = (0.0136D0/PTOT) * DSQRT(AL) * (1 + 0.038D0*DLOG(AL))
5158 * 8eme plan
5159       XP(8) = XP(7)+(ZPLANEP(8)-ZPLANEP(7))*VECT(4)/VECT(6)
5160       YP(8) = YP(7)+(ZPLANEP(8)-ZPLANEP(7))*VECT(5)/VECT(6)
5161       AL   = THICK/ DCOS(THET)
5162       AP(8)  = (0.0136D0/PTOT) * DSQRT(AL) * (1 + 0.038D0*DLOG(AL))
5163
5164 * 9eme plan
5165       XP(9) = XP(8)+(ZPLANEP(9)-ZPLANEP(8))*VECT(4)/VECT(6)
5166       YP(9) = YP(8)+(ZPLANEP(9)-ZPLANEP(8))*VECT(5)/VECT(6)
5167       AL   = THICK/ DCOS(THET)
5168       AP(9)  = (0.0136D0/PTOT) * DSQRT(AL) * (1 + 0.038D0*DLOG(AL))
5169
5170 * 10eme plan
5171       XP(10) = XP(9)+(ZPLANEP(10)-ZPLANEP(9))*VECT(4)/VECT(6)
5172       YP(10) = YP(9)+(ZPLANEP(10)-ZPLANEP(9))*VECT(5)/VECT(6)
5173       AL   = THICK/ DCOS(THET)
5174       AP(10)  = (0.0136D0/PTOT) * DSQRT(AL) * (1 + 0.038D0*DLOG(AL))
5175
5176 ** Matrice de covariance
5177       I = 0
5178       DO II = 1,NPLANE
5179         IF (LPLANE(II).EQ.1) THEN
5180         I = I + 1
5181 *        I = II
5182         J = I - 1
5183         DO JJ = II, NPLANE
5184            IF (LPLANE(JJ).EQ.1) THEN
5185            J = J + 1
5186 *           J = JJ
5187            COV (I,J) = 0.0D0
5188            COV (J,I) = A12
5189            IF (I .EQ. J) THEN
5190                  COV(J,I) =COV(J,I) + XPREC**2
5191            ENDIF
5192
5193 *           IF (I .EQ. 10 .AND. J .EQ. 10) PRINT *,'10 10   ',COV(J,I)
5194            DO L = 1,NTMAX
5195               COV(J,I) = COV(J,I)
5196      &  +  (ZPLANEP(II) + DISTAZ(L))*(ZPLANEP(JJ) + DISTAZ(L))*AMS(L)**2
5197            ENDDO
5198            DO K = 1, II-1
5199               COV(J,I) = COV(J,I)
5200      &  + (ZPLANEP(II)-ZPLANEP(K))*(ZPLANEP(JJ)-ZPLANEP(K))*AP(K)**2
5201 *              IF (I .EQ. 10 .AND. J .EQ. 10) PRINT *,'10 10   ',COV(J,I)
5202            ENDDO
5203            COVY(I,J) = 0.0D0
5204            COVY(J,I) = COV(J,I)
5205            IF (I .EQ. J) THEN
5206                  COVY(J,I) = COVY(J,I) - XPREC**2 + YPREC**2
5207            ENDIF
5208         ENDIF
5209         ENDDO
5210         ENDIF
5211       ENDDO
5212
5213 *  Inversion des matrices de covariance
5214       NPLU = I
5215
5216       IFAIL = 0
5217       CALL DSINV(NPLU, COV, NPLANE, IFAIL)
5218 **      IF (JFAIL.NE.0 .AND. IFAIL .NE. 0) STOP 'ERROR'
5219       IF (IFAIL .NE. 0) STOP 'ERROR'
5220       IFAIL = 0
5221       CALL DSINV(NPLU, COVY, NPLANE, IFAIL)
5222 **      IF (JFAIL.NE.0 .AND. IFAIL .NE. 0) STOP 'ERROR'
5223       IF (IFAIL .NE. 0) STOP 'ERROR'
5224 *      PRINT *,' COVARIANCE MATRIX AFTER'
5225 *      DO I = 1, NPLANE
5226 *         PRINT *,(COV(J,I),J=1,NPLANE)
5227 *      ENDDO
5228
5229 ** Calcul de FVAL ou CHI2
5230       FVAL = 0.0D0
5231       I = 0
5232       DO II = 1,NPLANE
5233         IF (LPLANE(II).EQ.1) THEN
5234         I = I+1
5235 *        I = II
5236         J = 0
5237         DO JJ = 1,NPLANE
5238            IF (LPLANE(JJ).EQ.1) THEN
5239               J = J+1
5240 *             J = JJ
5241               FVAL = FVAL + COV(J,I)*(XMP(II)-XP(II))*(XMP(JJ)-XP(JJ))
5242               FVAL = FVAL + COVY(J,I)*(YMP(II)-YP(II))
5243      &                               *(YMP(JJ)-YP(JJ))
5244 **             IF (JJ.EQ.II) THEN
5245 **                 FVAL = FVAL + (XMP(II)-XP(II))*(XMP(JJ)-XP(JJ))/XPREC**2
5246 **                 FVAL = FVAL + (YMP(II)-YP(II))
5247 **     &                               *(YMP(JJ)-YP(JJ))/YPREC**2
5248 **             ENDIF
5249            ENDIF
5250         ENDDO
5251         ENDIF
5252       ENDDO
5253       CHI2 = FVAL
5254          print *,' fcnfit pxz tphi talam xvert yvert chi2',
5255      &       PXZEA,PHIEA,ALAMEA,
5256      &       XEA,YEA,CHI2/FLOAT(2*NPLU-5)
5257
5258  1000 FORMAT(I5,7F12.6)
5259
5260       RETURN
5261       END
5262
5263