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