/[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.21 - (hide 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 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     10 CONTINUE
100     C
101     IF (GOOD2.EQ.1) THEN
102     C
103     CHTRACK = 0
104     C
105     CALL VZERO(IWPL,2)
106 mocchiut 1.16 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 mocchiut 1.1 do m = 1, 5
115     al_p(m) = al_pp(t,m)
116 mocchiut 1.16 c print *,' al_p(',m,') = ',al_p(m)
117 mocchiut 1.1 enddo
118     if (al_p(5).eq.0.) THEN
119 mocchiut 1.4 PRINT *,' CALORIMETER - WARNING F77: track with R = 0, discarded'
120 mocchiut 1.1 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 mocchiut 1.17 c print *,'T Y PLANE I= ',I,' Z = ',DISTX
130 mocchiut 1.1 ELSE
131     DISTX = PIANO(I) - 5.81 + ZALIG
132 mocchiut 1.17 c print *,'T X PLANE I= ',I,' Z = ',DISTX
133 mocchiut 1.1 ENDIF
134     ZIN(I) = distx / 10.
135 mocchiut 1.16 c print *,' ZIN(',I,') = ',ZIN(I)
136 mocchiut 1.1 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 mocchiut 1.3 c print *,' CALORIMETER - WARNING F77: tracking failed '
144 mocchiut 1.1 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 mocchiut 1.18 c IF (M.EQ.2) NN = 1
151 mocchiut 1.1 IF (MOD(I,2).EQ.NN) THEN
152 mocchiut 1.16 IF (REVERSE.EQ.0) THEN
153 mocchiut 1.19 SHIFT = -0.5
154     ELSE
155 mocchiut 1.16 SHIFT = +0.5
156     ENDIF
157 mocchiut 1.1 ELSE
158 mocchiut 1.16 IF (REVERSE.EQ.0) THEN
159 mocchiut 1.19 SHIFT = +0.5
160     ELSE
161 mocchiut 1.16 SHIFT = -0.5
162     ENDIF
163 mocchiut 1.1 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 mocchiut 1.3 c print *,
169     c & ' CALORIMETER - WARNING F77: tracking error (NaN values)'
170 mocchiut 1.1 GOOD2 = 0
171     GOTO 969
172     ENDIF
173     C
174     CX = XOUT(I)*10. + XALIG
175 mocchiut 1.5 CY = YOUT(I)*10. + YALIG
176 mocchiut 1.1 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 mocchiut 1.16 IF (I.EQ.NPLA) MX=ABS(Y(1)-Y(NPLA))/
186     & ABS(ZIN(1)-ZIN(NPLA))
187 mocchiut 1.1 ELSE
188     YY(I) = CY
189     BAR(M,I) = YY(I)
190 mocchiut 1.5 TBAR(M,I) = (-YALIG + YY(I))/10.
191 mocchiut 1.16 IF (I.EQ.NPLA) MY=ABS(Y(1)-Y(NPLA))/
192     & ABS(ZIN(1)-ZIN(NPLA))
193 mocchiut 1.1 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 mocchiut 1.7 IF (TRIGTY.GE.2.AND.HZN.EQ.0) THEN
247 mocchiut 1.1 C
248     C CALL SELFTRIGGER SUBROUTINE
249     C
250 mocchiut 1.8 CALL VZERO(IWPL,2)
251 mocchiut 1.7 CALL VZERO(VARCFIT,2)
252     CALL VZERO(NPCFIT,2)
253 mocchiut 1.16 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 mocchiut 1.7 C
262 mocchiut 1.1 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 mocchiut 1.13 c print *,' ax ',ax,' ay ',ay
272     c print *,' bx ',bx,' by ',by
273 mocchiut 1.1 IF (NPCFIT(M).GE.2) THEN
274     IF (M.EQ.1) THEN
275     VARCFIT(1) = CHI2X
276 mocchiut 1.20 IMPX = AX + BX * (ZALIG/10.) ! PAMELA REF
277 mocchiut 1.1 TANX = BX
278     ELSE
279     VARCFIT(2) = CHI2Y
280 mocchiut 1.20 IMPY = AY + BY * (ZALIG/10.) ! PAMELA REF
281 mocchiut 1.1 TANY = BY
282     ENDIF
283     C
284     DO I = 1,NPLA
285     NN = 0
286 mocchiut 1.18 c IF (M.EQ.2) NN = 1
287 mocchiut 1.1 IF (MOD(I,2).EQ.NN) THEN
288 mocchiut 1.16 IF (REVERSE.EQ.0) THEN
289 mocchiut 1.19 SHIFT = -0.5
290     ELSE
291 mocchiut 1.16 SHIFT = +0.5
292     ENDIF
293 mocchiut 1.1 ELSE
294 mocchiut 1.16 IF (REVERSE.EQ.0) THEN
295 mocchiut 1.19 SHIFT = +0.5
296     ELSE
297 mocchiut 1.16 SHIFT = -0.5
298     ENDIF
299 mocchiut 1.1 ENDIF
300     C
301     IF (M.EQ.1) THEN
302     DISTX = PIANO(I) - 5.81
303 mocchiut 1.20 Y(I) = (DISTX * TANX) + IMPX*10. !- XALIG
304 mocchiut 1.10 c CBAR(M,I) = Y(I)
305 mocchiut 1.20 BAR(M,I) = Y(I) + XALIG
306     CBAR(M,I) = Y(I) / 10.
307 mocchiut 1.16 IF (I.EQ.NPLA) MX=ABS(Y(1)-Y(NPLA))/
308     & ABS(ZIN(1)-ZIN(NPLA))
309 mocchiut 1.1 C
310     ELSE
311     DISTY = PIANO(I)
312 mocchiut 1.20 YY(I) = (DISTY * TANY) + IMPY*10. !- YALIG
313     c print *,' I ',i,' YY ',YY(I)
314 mocchiut 1.10 c CBAR(M,I) = YY(I)
315 mocchiut 1.20 BAR(M,I) = YY(I) + YALIG
316     CBAR(M,I) = YY(I) / 10.
317 mocchiut 1.16 IF (I.EQ.NPLA) MY=ABS(Y(1)-Y(NPLA))/
318     & ABS(ZIN(1)-ZIN(NPLA))
319 mocchiut 1.1 C
320     ENDIF
321     CALL LASTRISCIA(BAR(M,I),IBAR(M,I))
322     cibar(M,I) = ibar(m,i)
323 mocchiut 1.8 IF (ibar(m,i).EQ.-1) THEN
324     CHTRACK = CHTRACK + 1
325     ELSE
326     IWPL(M) = IWPL(M) + 1
327     ENDIF
328 mocchiut 1.1 ENDDO
329     ENDIF
330     C
331     ENDDO
332     C
333     ENDIF
334 mocchiut 1.8 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 mocchiut 1.1 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 mocchiut 1.7 IF (CHTRACK.EQ.44) THEN ! CHTRACK is the number of planes not touched by the track
354 mocchiut 1.1 GOOD2 = 0
355 mocchiut 1.3 c PRINT *,' CALORIMETER - WARNING F77: track outside calorimeter'
356 mocchiut 1.1 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 mocchiut 1.7 IF (TRIGTY.GE.2.AND.HZN.EQ.0.AND.GOOD2.EQ.0) THEN
386 mocchiut 1.1 RIG = ELEN ! SELFTRIGGER RIGIDITY
387     IF ( RIG.EQ.0. ) THEN
388 mocchiut 1.8 GOOD2 = 1
389 mocchiut 1.1 PRINT *,' CALORIMETER - WARNING F77: ST track with R = 0'
390     GOTO 50
391     ENDIF
392     ENDIF
393     C
394 mocchiut 1.7 IF (GOOD2.EQ.0.AND.(TRIGTY.LT.2.OR.HZN.EQ.1)) THEN
395     RIG = RIGINPUT
396     ENDIF
397     C
398 mocchiut 1.1 RNSS = 0.
399     QTOTT = 0.
400     C
401 mocchiut 1.15 PPLANEMAX = 1.01*(LOG(ABS(RIG)/0.0081)-1.) / 0.74
402 mocchiut 1.1 C
403     IPLANE = INT(ANINT(PPLANEMAX)) + 5
404     C
405     IF (IPLANE.GT.NPLA) IPLANE=NPLA
406     IF (IPLANE.LT.1) IPLANE = 1
407 mocchiut 1.16 c print *,' calcolo...'
408 mocchiut 1.1 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 mocchiut 1.15 RNSS = 0. ! BACO!!
432     QTOTT = 0. ! BACO!!
433 mocchiut 1.1 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 mocchiut 1.21 MNPLA = NPLA -4
611     IF ( MNPLA .LT. 1 ) MNPLA = 1
612     DO J = MNPLA,NPLA
613 mocchiut 1.1 NNX = IBAR(1,J)
614     IF (NNX.NE.-1) THEN
615 mocchiut 1.2 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 mocchiut 1.1 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 mocchiut 1.2 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 mocchiut 1.1 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 mocchiut 1.12 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 mocchiut 1.1 C
673     50 CONTINUE
674     C
675 mocchiut 1.16 c print *,' esco'
676 mocchiut 1.1 RETURN
677     END
678    
679    
  ViewVC Help
Powered by ViewVC 1.1.23