/[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.23 - (show annotations) (download)
Wed May 23 12:28:33 2012 UTC (12 years, 6 months ago) by mocchiut
Branch: MAIN
Changes since 1.22: +3 -3 lines
Fortran printouts commented

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

  ViewVC Help
Powered by ViewVC 1.1.23