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