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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 1.22 - (hide annotations) (download)
Tue Aug 4 14:01:18 2009 UTC (15 years, 4 months ago) by mocchiut
Branch: MAIN
CVS Tags: v9r00, v9r01
Changes since 1.21: +12 -9 lines 
Changed to work with GCC 4.x (gfortran) + ROOT >= 5.24
1 mocchiut 1.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 mocchiut 1.16 & xout(nplav),yout(nplav),zin(nplav)
10 mocchiut 1.1 C
11 mocchiut 1.16 REAL PIANO(NPLAV), VARFIT(2)
12 mocchiut 1.1 REAL TX, TY, SHIFT
13 mocchiut 1.16 REAL BAR(2,NPLAV), DISTY
14     REAL DISTX, Y(NPLAV), YY(NPLAV)
15 mocchiut 1.1 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 mocchiut 1.21 INTEGER Nfitx,Nfity, MNPLA
25 mocchiut 1.1 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 mocchiut 1.16 INTEGER IBAR(2,NPLAV), NPFIT(2), CHTRACK,IWPLU
30 mocchiut 1.1 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 mocchiut 1.8 COMMON / CALOFIT / VARFIT, NPFIT, IWPL,CHTRACK
59 mocchiut 1.1 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 mocchiut 1.16 c print *,' sono qui'
70 mocchiut 1.1 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 mocchiut 1.8 C SELFTRIGGER = 0
94 mocchiut 1.1 C
95 mocchiut 1.18 C BEGIN WITH THE FIRST TRACK IF WE HAVE A TRACK FROM TRACKER
96 mocchiut 1.1 C
97     T = 1
98     C
99 mocchiut 1.22 c 10 CONTINUE
100     CONTINUE
101 mocchiut 1.1 C
102     IF (GOOD2.EQ.1) THEN
103     C
104     CHTRACK = 0
105     C
106     CALL VZERO(IWPL,2)
107 mocchiut 1.16 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 mocchiut 1.1 do m = 1, 5
116     al_p(m) = al_pp(t,m)
117 mocchiut 1.16 c print *,' al_p(',m,') = ',al_p(m)
118 mocchiut 1.1 enddo
119     if (al_p(5).eq.0.) THEN
120 mocchiut 1.4 PRINT *,' CALORIMETER - WARNING F77: track with R = 0, discarded'
121 mocchiut 1.1 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 mocchiut 1.17 c print *,'T Y PLANE I= ',I,' Z = ',DISTX
131 mocchiut 1.1 ELSE
132     DISTX = PIANO(I) - 5.81 + ZALIG
133 mocchiut 1.17 c print *,'T X PLANE I= ',I,' Z = ',DISTX
134 mocchiut 1.1 ENDIF
135     ZIN(I) = distx / 10.
136 mocchiut 1.16 c print *,' ZIN(',I,') = ',ZIN(I)
137 mocchiut 1.1 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 mocchiut 1.3 c print *,' CALORIMETER - WARNING F77: tracking failed '
145 mocchiut 1.1 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 mocchiut 1.18 c IF (M.EQ.2) NN = 1
152 mocchiut 1.1 IF (MOD(I,2).EQ.NN) THEN
153 mocchiut 1.16 IF (REVERSE.EQ.0) THEN
154 mocchiut 1.19 SHIFT = -0.5
155     ELSE
156 mocchiut 1.16 SHIFT = +0.5
157     ENDIF
158 mocchiut 1.1 ELSE
159 mocchiut 1.16 IF (REVERSE.EQ.0) THEN
160 mocchiut 1.19 SHIFT = +0.5
161     ELSE
162 mocchiut 1.16 SHIFT = -0.5
163     ENDIF
164 mocchiut 1.1 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 mocchiut 1.3 c print *,
170     c & ' CALORIMETER - WARNING F77: tracking error (NaN values)'
171 mocchiut 1.1 GOOD2 = 0
172     GOTO 969
173     ENDIF
174     C
175     CX = XOUT(I)*10. + XALIG
176 mocchiut 1.5 CY = YOUT(I)*10. + YALIG
177 mocchiut 1.1 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 mocchiut 1.16 IF (I.EQ.NPLA) MX=ABS(Y(1)-Y(NPLA))/
187     & ABS(ZIN(1)-ZIN(NPLA))
188 mocchiut 1.1 ELSE
189     YY(I) = CY
190     BAR(M,I) = YY(I)
191 mocchiut 1.5 TBAR(M,I) = (-YALIG + YY(I))/10.
192 mocchiut 1.16 IF (I.EQ.NPLA) MY=ABS(Y(1)-Y(NPLA))/
193     & ABS(ZIN(1)-ZIN(NPLA))
194 mocchiut 1.1 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 mocchiut 1.7 IF (TRIGTY.GE.2.AND.HZN.EQ.0) THEN
248 mocchiut 1.1 C
249     C CALL SELFTRIGGER SUBROUTINE
250     C
251 mocchiut 1.8 CALL VZERO(IWPL,2)
252 mocchiut 1.7 CALL VZERO(VARCFIT,2)
253     CALL VZERO(NPCFIT,2)
254 mocchiut 1.16 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 mocchiut 1.7 C
263 mocchiut 1.1 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 mocchiut 1.13 c print *,' ax ',ax,' ay ',ay
273     c print *,' bx ',bx,' by ',by
274 mocchiut 1.1 IF (NPCFIT(M).GE.2) THEN
275     IF (M.EQ.1) THEN
276     VARCFIT(1) = CHI2X
277 mocchiut 1.20 IMPX = AX + BX * (ZALIG/10.) ! PAMELA REF
278 mocchiut 1.1 TANX = BX
279     ELSE
280     VARCFIT(2) = CHI2Y
281 mocchiut 1.20 IMPY = AY + BY * (ZALIG/10.) ! PAMELA REF
282 mocchiut 1.1 TANY = BY
283     ENDIF
284     C
285     DO I = 1,NPLA
286     NN = 0
287 mocchiut 1.18 c IF (M.EQ.2) NN = 1
288 mocchiut 1.1 IF (MOD(I,2).EQ.NN) THEN
289 mocchiut 1.16 IF (REVERSE.EQ.0) THEN
290 mocchiut 1.19 SHIFT = -0.5
291     ELSE
292 mocchiut 1.16 SHIFT = +0.5
293     ENDIF
294 mocchiut 1.1 ELSE
295 mocchiut 1.16 IF (REVERSE.EQ.0) THEN
296 mocchiut 1.19 SHIFT = +0.5
297     ELSE
298 mocchiut 1.16 SHIFT = -0.5
299     ENDIF
300 mocchiut 1.1 ENDIF
301     C
302     IF (M.EQ.1) THEN
303     DISTX = PIANO(I) - 5.81
304 mocchiut 1.20 Y(I) = (DISTX * TANX) + IMPX*10. !- XALIG
305 mocchiut 1.10 c CBAR(M,I) = Y(I)
306 mocchiut 1.20 BAR(M,I) = Y(I) + XALIG
307     CBAR(M,I) = Y(I) / 10.
308 mocchiut 1.16 IF (I.EQ.NPLA) MX=ABS(Y(1)-Y(NPLA))/
309     & ABS(ZIN(1)-ZIN(NPLA))
310 mocchiut 1.1 C
311     ELSE
312     DISTY = PIANO(I)
313 mocchiut 1.20 YY(I) = (DISTY * TANY) + IMPY*10. !- YALIG
314     c print *,' I ',i,' YY ',YY(I)
315 mocchiut 1.10 c CBAR(M,I) = YY(I)
316 mocchiut 1.20 BAR(M,I) = YY(I) + YALIG
317     CBAR(M,I) = YY(I) / 10.
318 mocchiut 1.16 IF (I.EQ.NPLA) MY=ABS(Y(1)-Y(NPLA))/
319     & ABS(ZIN(1)-ZIN(NPLA))
320 mocchiut 1.1 C
321     ENDIF
322     CALL LASTRISCIA(BAR(M,I),IBAR(M,I))
323     cibar(M,I) = ibar(m,i)
324 mocchiut 1.8 IF (ibar(m,i).EQ.-1) THEN
325     CHTRACK = CHTRACK + 1
326     ELSE
327     IWPL(M) = IWPL(M) + 1
328     ENDIF
329 mocchiut 1.1 ENDDO
330     ENDIF
331     C
332     ENDDO
333     C
334     ENDIF
335 mocchiut 1.8 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 mocchiut 1.1 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 mocchiut 1.7 IF (CHTRACK.EQ.44) THEN ! CHTRACK is the number of planes not touched by the track
355 mocchiut 1.1 GOOD2 = 0
356 mocchiut 1.3 c PRINT *,' CALORIMETER - WARNING F77: track outside calorimeter'
357 mocchiut 1.1 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     PRINT *,' CALORIMETER - WARNING F77: track with R = 0'
383     GOTO 50
384     ENDIF
385     ENDIF
386 mocchiut 1.7 IF (TRIGTY.GE.2.AND.HZN.EQ.0.AND.GOOD2.EQ.0) THEN
387 mocchiut 1.1 RIG = ELEN ! SELFTRIGGER RIGIDITY
388     IF ( RIG.EQ.0. ) THEN
389 mocchiut 1.8 GOOD2 = 1
390 mocchiut 1.1 PRINT *,' CALORIMETER - WARNING F77: ST track with R = 0'
391     GOTO 50
392     ENDIF
393     ENDIF
394     C
395 mocchiut 1.7 IF (GOOD2.EQ.0.AND.(TRIGTY.LT.2.OR.HZN.EQ.1)) THEN
396     RIG = RIGINPUT
397     ENDIF
398     C
399 mocchiut 1.1 RNSS = 0.
400     QTOTT = 0.
401     C
402 mocchiut 1.15 PPLANEMAX = 1.01*(LOG(ABS(RIG)/0.0081)-1.) / 0.74
403 mocchiut 1.1 C
404     IPLANE = INT(ANINT(PPLANEMAX)) + 5
405     C
406     IF (IPLANE.GT.NPLA) IPLANE=NPLA
407     IF (IPLANE.LT.1) IPLANE = 1
408 mocchiut 1.16 c print *,' calcolo...'
409 mocchiut 1.1 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 mocchiut 1.15 RNSS = 0. ! BACO!!
433     QTOTT = 0. ! BACO!!
434 mocchiut 1.1 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 mocchiut 1.22 IF (DEXY(1,J,I).GE.EMIN) THEN
485     NCYL = NCYL + 1
486     QCYL = QCYL + DEXY(1,J,I)
487     ENDIF
488     ENDDO
489 mocchiut 1.1 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 mocchiut 1.22 IF (DEXY(2,J,I).GE.EMIN) THEN
498     NCYL = NCYL + 1
499     QCYL = QCYL + DEXY(2,J,I)
500     ENDIF
501     ENDDO
502 mocchiut 1.1 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 mocchiut 1.21 MNPLA = NPLA -4
614     IF ( MNPLA .LT. 1 ) MNPLA = 1
615     DO J = MNPLA,NPLA
616 mocchiut 1.1 NNX = IBAR(1,J)
617     IF (NNX.NE.-1) THEN
618 mocchiut 1.2 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 mocchiut 1.1 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 mocchiut 1.2 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 mocchiut 1.1 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 mocchiut 1.12 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 mocchiut 1.1 C
676     50 CONTINUE
677     C
678 mocchiut 1.16 c print *,' esco'
679 mocchiut 1.1 RETURN
680     END
681    
682    
  ViewVC Help
Powered by ViewVC 1.1.23