/[PAMELA software]/DarthVader/CalorimeterLevel2/src/calol2tr.for
ViewVC logotype

Contents of /DarthVader/CalorimeterLevel2/src/calol2tr.for

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.7 - (show annotations) (download)
Tue Nov 14 14:08:50 2006 UTC (18 years ago) by mocchiut
Branch: MAIN
CVS Tags: v2r01
Changes since 1.6: +29 -22 lines
Major calorimeter release, some news in DarthVader main

1 *****************************************************************************
2 INTEGER FUNCTION CALOL2TR()
3 c
4 IMPLICIT NONE
5 C
6 INCLUDE 'INTEST.TXT'
7 C
8 DOUBLE PRECISION al_p(5),
9 & xout(npla),yout(npla),zin(npla)
10 C
11 REAL PIANO(22), VARFIT(2)
12 REAL TX, TY, SHIFT
13 REAL BAR(2,NPLA), DISTY
14 REAL DISTX, Y(NPLA), YY(NPLA)
15 REAL RIG, PPLANEMAX, RMASS
16 REAL RNSS, QTOTT, RQT, MX, MY
17 REAL CHECK, ENER, CX, CY
18 REAL EINF, ESUP, RPIANO(2)
19 REAL hmemor(9000000), X01PL
20 C
21 REAL ax,bx,eax,ebx,chi2x
22 REAL ay,by,eay,eby,chi2y
23 REAL parzen3, TMISD
24 INTEGER Nfitx,Nfity
25 C
26 INTEGER INDEX, NTOT(2), NPIANI, GTR
27 INTEGER j, m, i, IWPL(2), timpx, timpy, T, nn
28 INTEGER IPLANE, NNX, NNY, INFX, INFY, ISUPX, ISUPY
29 INTEGER IBAR(2,NPLA), NPFIT(2), CHTRACK,IWPLU
30 INTEGER Iquest(100), ICONTROL5, nin, IFAIL
31 C
32 PARAMETER (X01PL=0.74)
33 C
34
35 C
36 COMMON / slftrig / tmisd,ax,bx,eax,ebx,chi2x,Nfitx,ay,by,eay,eby,
37 & chi2y,Nfity,parzen3
38 SAVE / slftrig /
39 C
40 COMMON / TAGLIOEN / EINF, ESUP, ENER(2)
41 SAVE / TAGLIOEN /
42 C
43 COMMON / SHIFT / SHIFT
44 SAVE / SHIFT /
45 C
46 COMMON / ANGOLO / BAR, IBAR
47 SAVE / ANGOLO /
48 C
49 COMMON / WHERE / CX, CY, PIANO
50 SAVE / WHERE /
51 C
52 COMMON / GENERAL / RIG, RMASS
53 SAVE / GENERAL /
54 C
55 COMMON / CH / CHECK
56 SAVE / CH /
57 C
58 COMMON / CALOFIT / VARFIT, NPFIT
59 SAVE / CALOFIT /
60 C
61 COMMON / pawcd / hmemor
62 save / pawcd /
63 C
64 Common / QUESTd / Iquest
65 save / questd /
66 C
67 C Begin !
68 C
69 CALOL2TR = 0;
70 NCORE = 0.
71 QCORE = 0.
72 NOINT = 0.
73 QCYL = 0.
74 NCYL = 0.
75 QLOW = 0.
76 NLOW = 0.
77 QTR = 0.
78 NTR = 0.
79 QLAST = 0.
80 QTRACK = 0.
81 QPRESH = 0.
82 NPRESH = 0.
83 QTRACKX = 0.
84 QTRACKY = 0.
85 DXTRACK = 0.
86 DYTRACK = 0.
87 QPRE = 0.
88 NPRE = 0.
89 NLAST = 0.
90 PLANETOT = 0.
91 QMEAN = 0.
92 SELFTRIGGER = 0
93 C
94 C BEGIN WITH THE FISRT TRACK IF WE HAVE A TRACK FROM TRACKER
95 C
96 T = 1
97 C
98 10 CONTINUE
99 C
100 IF (GOOD2.EQ.1) THEN
101 C
102 CHTRACK = 0
103 C
104 CALL VZERO(IWPL,2)
105 CALL VZERO(BAR,2*NPLA)
106 CALL VZERO(IBAR,2*NPLA)
107 CALL VZERO(TBAR,2*NPLA)
108 CALL VZERO(TIBAR,2*NPLA)
109 CALL VZERO(Y,NPLA)
110 CALL VZERO(YY,NPLA)
111 CALL VZERO(XOUT,NPLA)
112 CALL VZERO(YOUT,NPLA)
113 do m = 1, 5
114 al_p(m) = al_pp(t,m)
115 enddo
116 if (al_p(5).eq.0.) THEN
117 PRINT *,' CALORIMETER - WARNING F77: track with R = 0, discarded'
118 GOOD2 = 0
119 GOTO 969
120 ENDIF
121 DO M = 1,2
122 DO I = 1,NPLA
123 XOUT(I) = 0.
124 YOUT(I) = 0.
125 IF (MOD(M,2).EQ.0) THEN
126 DISTX = PIANO(I) + ZALIG
127 ELSE
128 DISTX = PIANO(I) - 5.81 + ZALIG
129 ENDIF
130 ZIN(I) = distx / 10.
131 TBAR(M,I) = 0.
132 TIBAR(M,I) = 0
133 enddo
134 IFAIL = 0
135 call DOTRACK(NPLA,ZIN,XOUT,YOUT,AL_P,IFAIL)
136 if(IFAIL.ne.0)then
137 GOOD2 = 0
138 c print *,' CALORIMETER - WARNING F77: tracking failed '
139 goto 969
140 endif
141 TX = TAN(ASIN(AL_P(3))) * COS(AL_P(4))
142 TY = TAN(ASIN(AL_P(3))) * SIN(AL_P(4))
143 DO I = 1, NPLA
144 NN = 0
145 C????? IF (M.EQ.2) NN = 1
146 IF (MOD(I,2).EQ.NN) THEN
147 SHIFT = +0.5
148 ELSE
149 SHIFT = -0.5
150 ENDIF
151 C
152 C CHECK IF XOUT OR YOUT ARE NaN
153 C
154 IF (XOUT(I).NE.XOUT(I).OR.YOUT(I).NE.YOUT(I)) THEN
155 c print *,
156 c & ' CALORIMETER - WARNING F77: tracking error (NaN values)'
157 GOOD2 = 0
158 GOTO 969
159 ENDIF
160 C
161 CX = XOUT(I)*10. + XALIG
162 CY = YOUT(I)*10. + YALIG
163 C
164 IF (I.EQ.1) THEN
165 TIMPX = CX
166 TIMPY = CY
167 ENDIF
168 IF (M.EQ.1) THEN
169 Y(I) = CX
170 BAR(M,I) = Y(I)
171 TBAR(M,I) = (Y(I) - XALIG)/10.
172 IF (I.EQ.22) MX=ABS(Y(1)-Y(22))/ABS(ZIN(1)-ZIN(22))
173 ELSE
174 YY(I) = CY
175 BAR(M,I) = YY(I)
176 TBAR(M,I) = (-YALIG + YY(I))/10.
177 IF (I.EQ.22) MY=ABS(Y(1)-Y(22))/ABS(ZIN(1)-ZIN(22))
178 ENDIF
179 CALL LASTRISCIA(BAR(M,I),IBAR(M,I))
180 tibar(M,I) = ibar(m,i)
181 IF (ibar(m,i).EQ.-1) THEN
182 CHTRACK = CHTRACK + 1
183 ELSE
184 IWPL(M) = IWPL(M) + 1
185 ENDIF
186 ENDDO
187 ENDDO
188 969 continue
189 cC
190 cC IF WE HAVE A GOOD CALORIMETER FIT DOES IT MATCH WITH TRACKER FIT?
191 cC
192 c IF (GOOD2.EQ.1.AND.NPFIT(2).GT.15.AND.VARFIT(2).LT.1000
193 c & .AND.TRKCHI2.EQ.1) THEN
194 c IF (ABS(TBAR(2,1)-CBAR(2,1))<40.) THEN
195 cC
196 cC GOOD, THE TWO TRACKS COINCIDE
197 cC
198 c IF (T.EQ.2) TRKCHI2 = 2
199 c GOTO 6996
200 c ELSE
201 cC
202 cC IT IS NOT A GOOD FIT BUT WE HAVE AN IMAGE AND IT IS THE FIRST TRACK
203 cC
204 c IF (T.EQ.1) THEN
205 c T = 2
206 c GOTO 10
207 c ENDIF
208 c IF (T.EQ.2) THEN
209 c TRKCHI2 = -1
210 c T = 1
211 c GOTO 10
212 c ENDIF
213 c ENDIF
214 c ENDIF
215 C
216 IF (GOOD2.EQ.0) THEN
217 c IF (T.EQ.1.AND.TRKCHI2.EQ.1) THEN
218 c GOOD2 = 1
219 c T = 2
220 c GOTO 10
221 c ENDIF
222 GOTO 50
223 ENDIF
224 C
225 GOTO 6996
226 C
227 ENDIF
228 C
229 C WE MUST PROCESS A SELFTRIGGER EVENT
230 C
231 IF (TRIGTY.GE.2.AND.HZN.EQ.0) THEN
232 C
233 C CALL SELFTRIGGER SUBROUTINE
234 C
235 CALL VZERO(VARCFIT,2)
236 CALL VZERO(NPCFIT,2)
237 CALL VZERO(TBAR,2*NPLA)
238 CALL VZERO(TIBAR,2*NPLA)
239 CALL VZERO(BAR,2*NPLA)
240 CALL VZERO(IBAR,2*NPLA)
241 CALL VZERO(Y,NPLA)
242 CALL VZERO(YY,NPLA)
243 CALL VZERO(XOUT,NPLA)
244 CALL VZERO(YOUT,NPLA)
245 C
246 CALL SELFTRIG()
247 ELEN = PARZEN3
248 SELEN = ABS(ELEN) * (11.98*1E-2 + 7.6 * EXP(-5736/ABS(ELEN)))
249 C
250 NPCFIT(1) = NFITX
251 NPCFIT(2) = NFITY
252 C
253 DO M = 1,2
254 C
255 IF (NPCFIT(M).GE.2) THEN
256 IF (M.EQ.1) THEN
257 VARCFIT(1) = CHI2X
258 IMPX = 10. * ( AX + 12.1 )
259 TANX = BX
260 ELSE
261 VARCFIT(2) = CHI2Y
262 IMPY = 10. * ( AY + 12.2 )
263 TANY = BY
264 ENDIF
265 C
266 DO I = 1,NPLA
267 NN = 0
268 C????? IF (M.EQ.2) NN = 1
269 IF (MOD(I,2).EQ.NN) THEN
270 SHIFT = +0.5
271 ELSE
272 SHIFT = -0.5
273 ENDIF
274 C
275 IF (M.EQ.1) THEN
276 DISTX = PIANO(I) - 5.81
277 Y(I) = DISTX * TANX + CX
278 BAR(M,I) = Y(I)
279 CBAR(M,I) = Y(I)
280 IF (I.EQ.22) MX=ABS(Y(1)-Y(22))/ABS(ZIN(1)-ZIN(22))
281 C
282 ELSE
283 DISTY = PIANO(I)
284 YY(I) = DISTY * TANY + CY
285 BAR(M,I) = YY(I)
286 CBAR(M,I) = YY(I)
287 IF (I.EQ.22) MY=ABS(Y(1)-Y(22))/ABS(ZIN(1)-ZIN(22))
288 C
289 ENDIF
290 CALL LASTRISCIA(BAR(M,I),IBAR(M,I))
291 cibar(M,I) = ibar(m,i)
292 ENDDO
293 ENDIF
294 C
295 ENDDO
296 C
297 ENDIF
298 C IF (TRIGTY.GE.2.AND.HZN.NE.0) THEN
299 C IF (GOOD2.EQ.0) THEN
300 C PRINT *,' CALORIMETER - WARNING F77: unknown request'
301 C GOOD2 = 1
302 C GOTO 50
303 C ENDIF
304 C ENDIF
305 C
306 6996 CONTINUE
307 C
308 DX0L = 0.
309 C
310 C IF THE TRACK IS OUTSIDE THE CALORIMETER GO OUT, IF NOT CALCULATE DX0L
311 C
312 IF (CHTRACK.EQ.44) THEN ! CHTRACK is the number of planes not touched by the track
313 GOOD2 = 0
314 c PRINT *,' CALORIMETER - WARNING F77: track outside calorimeter'
315 GOTO 50
316 ELSE
317 IF ( IWPL(1).LE.IWPL(2) ) THEN
318 IWPLU = IWPL(1)
319 ELSE
320 IWPLU = IWPL(2)
321 ENDIF
322 C
323 DX0L = IWPLU * SQRT((BAR(2,1)-(2.66*MY+BAR(2,1)))**2
324 & + (BAR(1,1)-(2.66*MX+BAR(1,1)))**2 + 2.66**2) /
325 & 3.6
326 C
327 ENDIF
328 C
329 C
330 C RIG IS RIGIDITY AS DETERMINED BY THE TRACKER
331 C OR by CALORIMETER IF IN SELFTRIGGER MODE
332 C
333 IF (GOOD2.EQ.1) THEN
334 GTR = 1
335 IF (TRKCHI2.LT.0) GTR = 2
336 IF ( AL_PP(GTR,5).NE.0. ) THEN
337 RIG = 1./(AL_PP(GTR,5))
338 ELSE
339 GOOD2 = 0
340 PRINT *,' CALORIMETER - WARNING F77: track with R = 0'
341 GOTO 50
342 ENDIF
343 ENDIF
344 IF (TRIGTY.GE.2.AND.HZN.EQ.0.AND.GOOD2.EQ.0) THEN
345 RIG = ELEN ! SELFTRIGGER RIGIDITY
346 IF ( RIG.EQ.0. ) THEN
347 GOOD2 = 0
348 PRINT *,' CALORIMETER - WARNING F77: ST track with R = 0'
349 GOTO 50
350 ENDIF
351 ENDIF
352 C
353 IF (GOOD2.EQ.0.AND.(TRIGTY.LT.2.OR.HZN.EQ.1)) THEN
354 RIG = RIGINPUT
355 ENDIF
356 C
357 RNSS = 0.
358 QTOTT = 0.
359 C
360 PPLANEMAX = 1.01*(LOG(ABS(RIG)/0.0081)-1.)
361 C
362 IPLANE = INT(ANINT(PPLANEMAX)) + 5
363 C
364 IF (IPLANE.GT.NPLA) IPLANE=NPLA
365 IF (IPLANE.LT.1) IPLANE = 1
366 C
367 C CALCULATE QLOW AND NLOW
368 C
369 DO J = IPLANE,NPLA
370 DO I = 1,NCHA
371 IF (DEXY(1,J,I).GE.EMIN) THEN
372 NLOW = NLOW + 1
373 QLOW = QLOW + DEXY(1,J,I)
374 ENDIF
375 IF (DEXY(2,J,I).GE.EMIN) THEN
376 NLOW = NLOW + 1
377 QLOW = QLOW + DEXY(2,J,I)
378 ENDIF
379 ENDDO
380 ENDDO
381 C
382 C CALCULATE QCORE AND NCORE
383 C
384 C
385 C 8 STRIPS ARE 2.88 cm , A MOLIERE RADIUS IS ABOUT 0.7 cm .
386 C
387 DO J = 1,IPLANE
388 NNX = IBAR(1,J)
389 IF (NNX.NE.-1) THEN
390 IF (NNX.LT.9) NNX = 9
391 IF (NNX.GT.88) NNX = 88
392 INFX = NNX - 8
393 ISUPX = NNX + 8
394 DO I = INFX,ISUPX
395 IF (DEXY(1,J,I).GE.EMIN) THEN
396 RNSS = RNSS + 1
397 QTOTT = QTOTT + DEXY(1,J,I)
398 ENDIF
399 ENDDO
400 ENDIF
401 C
402 NNY = IBAR(2,J)
403 IF (NNY.NE.-1) THEN
404 IF (NNY.LT.9) NNY = 9
405 IF (NNY.GT.88) NNY = 88
406 INFY = NNY - 8
407 ISUPY = NNY + 8
408 DO I = INFY,ISUPY
409 IF (DEXY(2,J,I).GE.EMIN) THEN
410 RNSS = RNSS + 1
411 QTOTT = QTOTT + DEXY(2,J,I)
412 ENDIF
413 ENDDO
414 ENDIF
415 NCORE = RNSS * FLOAT(J) + NCORE
416 QCORE = QTOTT * FLOAT(J) + QCORE
417 ENDDO
418 C
419 C CALCULATE NOINT
420 C
421 CALL NOINTER(NIN)
422 NOINT = FLOAT(NIN)
423 C
424 C
425 C QCYL = DETECTED ENERGY AND NCYL = NUMBER OF HIT STRIPS IN A CYLINDER oF
426 C RADIUS 8.5 STRIPS WITH AXIS DEFINED BY THE DIRECTION OF THE INCOMING
427 C PARTICLE .
428 C
429 C 8 STRIPS ARE 2.88 cm , A MOLIERE RADIUS IS ABOUT 0.7 cm .
430 C
431 DO J = 1,NPLA
432 NNX = IBAR(1,J)
433 IF (NNX.NE.-1) THEN
434 IF (NNX.LT.9) NNX = 9
435 IF (NNX.GT.88) NNX = 88
436 INFX = NNX - 8
437 ISUPX = NNX + 8
438 DO I = INFX,ISUPX
439 IF (DEXY(1,J,I).LT.EMIN) GO TO 710
440 NCYL = NCYL + 1
441 QCYL = QCYL + DEXY(1,J,I)
442 710 ENDDO
443 ENDIF
444 NNY = IBAR(2,J)
445 IF (NNY.NE.-1) THEN
446 IF (NNY.LT.9) NNY = 9
447 IF (NNY.GT.88) NNY = 88
448 INFY = NNY - 8
449 ISUPY = NNY + 8
450 DO I=INFY,ISUPY
451 IF (DEXY(2,J,I).LT.EMIN) GO TO 810
452 NCYL = NCYL + 1
453 QCYL = QCYL + DEXY(2,J,I)
454 810 ENDDO
455 ENDIF
456 C
457 C QTR = DETECTED ENERGY AND NTR = NUMBER OF HIT STRIPS IN A CYLINDER oF
458 C RADIUS 4.5 STRIPS WITH AXIS DEFINED BY THE DIRECTION OF THE INCOMING
459 C PARTICLE .
460 C
461 NNX = IBAR(1,J)
462 IF (NNX.NE.-1) THEN
463 IF (NNX.LT.5) NNX = 5
464 IF (NNX.GT.92) NNX = 92
465 INFX = NNX - 4
466 ISUPX = NNX + 4
467 DO I = INFX,ISUPX
468 IF (DEXY(1,J,I).GT.EMIN) THEN
469 NTR = NTR + 1
470 QTR = QTR + DEXY(1,J,I)
471 ENDIF
472 ENDDO
473 ENDIF
474 C
475 NNY = IBAR(2,J)
476 IF (NNY.NE.-1) THEN
477 IF (NNY.LT.5) NNY = 5
478 IF (NNY.GT.92) NNY = 92
479 INFY = NNY - 4
480 ISUPY = NNY + 4
481 DO I = INFY, ISUPY
482 IF (DEXY(2,J,I).GT.EMIN) THEN
483 NTR = NTR + 1
484 QTR = QTR + DEXY(2,J,I)
485 ENDIF
486 ENDDO
487 ENDIF
488 ENDDO
489 C
490 C CALCULATE QTRACK
491 C
492 CALL LATERALE(QTRACK,RQT)
493
494 C
495 C CALCULATE NPRESH AND QPRESH
496 C
497 DO I = 1,4
498 NNX = IBAR(1,I)
499 IF (NNX.NE.-1) THEN
500 IF (NNX.LT.3) NNX = 3
501 IF (NNX.GT.94) NNX = 94
502 INFX = NNX - 2
503 ISUPX = NNX + 2
504 DO J = INFX,ISUPX
505 IF (DEXY(1,I,J).GE.EMIN) THEN
506 NPRESH = NPRESH + 1
507 QPRESH = QPRESH + DEXY(1,I,J)
508 ENDIF
509 ENDDO
510 ENDIF
511 C
512 NNY = IBAR(2,I)
513 IF (NNY.NE.-1) THEN
514 IF (NNY.LT.3) NNY = 3
515 IF (NNY.GT.94) NNY = 94
516 INFY = NNY - 2
517 ISUPY = NNY + 2
518 DO J = INFY,ISUPY
519 IF (DEXY(2,I,J).GE.EMIN) THEN
520 NPRESH = NPRESH + 1
521 QPRESH = QPRESH + DEXY(2,I,J)
522 ENDIF
523 ENDDO
524 ENDIF
525 ENDDO
526 C
527 C CALCULATE DXTRACK, DYTRACK, QTRACKX AND QTRACKY
528 C
529 ICONTROL5 = 0
530 CALL NSHOWER(ICONTROL5,DXTRACK,DYTRACK,QTRACKX,QTRACKY)
531 C
532 C CALCULATE QPRE AND NPRE
533 C
534 DO J = 1,3
535 NNX = IBAR(1,J)
536 IF (NNX.NE.-1) THEN
537 IF (NNX.LT.9) NNX = 9
538 IF (NNX.GT.88) NNX = 88
539 INFX = NNX - 8
540 ISUPX = NNX + 8
541 DO I = INFX,ISUPX
542 IF (DEXY(1,J,I).GE.EMIN) THEN
543 NPRE = NPRE + 1
544 QPRE = QPRE + DEXY(1,J,I)
545 ENDIF
546 ENDDO
547 ENDIF
548 C
549 NNY = IBAR(2,J)
550 IF (NNY.NE.-1) THEN
551 IF (NNY.LT.9) NNY = 9
552 IF (NNY.GT.88) NNY = 88
553 INFY = NNY - 8
554 ISUPY = NNY + 8
555 DO I=INFY,ISUPY
556 IF (DEXY(2,J,I).GE.EMIN) THEN
557 NPRE = NPRE + 1
558 QPRE = QPRE + DEXY(2,J,I)
559 ENDIF
560 ENDDO
561 ENDIF
562 ENDDO
563 C
564 C CALCULATE NLAST AND QLAST
565 C
566 DO J = NPLA-4,NPLA
567 NNX = IBAR(1,J)
568 IF (NNX.NE.-1) THEN
569 IF (NNX.LT.5) NNX = 5
570 IF (NNX.GT.92) NNX = 92
571 c IF (NNX.LT.9) NNX = 9
572 c IF (NNX.GT.88) NNX = 88
573 INFX = NNX - 4
574 ISUPX = NNX + 4
575 c INFX = NNX - 8
576 c ISUPX = NNX + 8
577 DO I = INFX,ISUPX
578 IF (DEXY(1,J,I).GE.EMIN) THEN
579 NLAST = NLAST + 1
580 QLAST = QLAST + DEXY(1,J,I)
581 ENDIF
582 ENDDO
583 ENDIF
584 C
585 NNY = IBAR(2,J)
586 IF (NNY.NE.-1) THEN
587 IF (NNY.LT.5) NNY = 5
588 IF (NNY.GT.92) NNY = 92
589 c IF (NNY.LT.9) NNY = 9
590 c IF (NNY.GT.88) NNY = 88
591 INFY = NNY - 4
592 ISUPY = NNY + 4
593 c INFY = NNY - 8
594 c ISUPY = NNY + 8
595 DO I=INFY,ISUPY
596 IF (DEXY(2,J,I).GE.EMIN) THEN
597 NLAST = NLAST + 1
598 QLAST = QLAST + DEXY(2,J,I)
599 ENDIF
600 ENDDO
601 ENDIF
602 ENDDO
603 C
604 EINF = EMIN
605 ESUP = 50.
606 C
607 C CALCULATE PLANETOT AND QMEAN
608 C
609 DO M = 1,2
610 RPIANO(M) = 0.
611 NTOT(M) = 0
612 ENDDO
613 NPIANI = 5
614 QMEAN = 0.
615 INDEX = 0
616 CALL ELIO(RPIANO,NPIANI,QMEAN,NTOT,INDEX)
617 PLANETOT = RPIANO(1) + RPIANO(2)
618 C
619 50 CONTINUE
620 C
621 RETURN
622 END
623
624

  ViewVC Help
Powered by ViewVC 1.1.23