/[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.21 - (show annotations) (download)
Fri Dec 19 13:39:34 2008 UTC (15 years, 11 months ago) by mocchiut
Branch: MAIN
CVS Tags: v6r01, v6r00
Changes since 1.20: +4 -2 lines
Bug when using small pre-sampler calorimeter (<=4 planes) fixed

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

  ViewVC Help
Powered by ViewVC 1.1.23