/[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.19 - (hide annotations) (download)
Thu Nov 29 16:11:53 2007 UTC (17 years ago) by mocchiut
Branch: MAIN
CVS Tags: v5r00
Changes since 1.18: +8 -8 lines 
More bugs related to plane shifting 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     INTEGER Nfitx,Nfity
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 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.13 IMPX = AX ! PAMELA REF
277 mocchiut 1.1 TANX = BX
278     ELSE
279     VARCFIT(2) = CHI2Y
280 mocchiut 1.13 IMPY = AY ! 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.13 Y(I) = (DISTX * TANX) + AX - XALIG
304 mocchiut 1.10 c CBAR(M,I) = Y(I)
305 mocchiut 1.1 BAR(M,I) = Y(I)
306 mocchiut 1.13 CBAR(M,I) = (Y(I) + XALIG)/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.13 YY(I) = (DISTY * TANY) + AY - YALIG
313 mocchiut 1.10 c CBAR(M,I) = YY(I)
314 mocchiut 1.1 BAR(M,I) = YY(I)
315 mocchiut 1.13 CBAR(M,I) = (YY(I) + YALIG)/10.
316 mocchiut 1.16 IF (I.EQ.NPLA) MY=ABS(Y(1)-Y(NPLA))/
317     & ABS(ZIN(1)-ZIN(NPLA))
318 mocchiut 1.1 C
319     ENDIF
320     CALL LASTRISCIA(BAR(M,I),IBAR(M,I))
321     cibar(M,I) = ibar(m,i)
322 mocchiut 1.8 IF (ibar(m,i).EQ.-1) THEN
323     CHTRACK = CHTRACK + 1
324     ELSE
325     IWPL(M) = IWPL(M) + 1
326     ENDIF
327 mocchiut 1.1 ENDDO
328     ENDIF
329     C
330     ENDDO
331     C
332     ENDIF
333 mocchiut 1.8 C
334     IF (TRIGTY.GE.2.AND.HZN.NE.0) THEN
335     IF (GOOD2.EQ.1) THEN
336     PRINT *,' CALORIMETER - WARNING F77: unknown request'
337     GOOD2 = 1
338     GOTO 50
339     ENDIF
340     IF ( NPCFIT(1).EQ.0.OR.NPCFIT(2).EQ.0 ) THEN
341     GOOD2 = 1
342     GOTO 50
343     ENDIF
344     ENDIF
345 mocchiut 1.1 C
346     6996 CONTINUE
347     C
348     DX0L = 0.
349     C
350     C IF THE TRACK IS OUTSIDE THE CALORIMETER GO OUT, IF NOT CALCULATE DX0L
351     C
352 mocchiut 1.7 IF (CHTRACK.EQ.44) THEN ! CHTRACK is the number of planes not touched by the track
353 mocchiut 1.1 GOOD2 = 0
354 mocchiut 1.3 c PRINT *,' CALORIMETER - WARNING F77: track outside calorimeter'
355 mocchiut 1.1 GOTO 50
356     ELSE
357     IF ( IWPL(1).LE.IWPL(2) ) THEN
358     IWPLU = IWPL(1)
359     ELSE
360     IWPLU = IWPL(2)
361     ENDIF
362     C
363     DX0L = IWPLU * SQRT((BAR(2,1)-(2.66*MY+BAR(2,1)))**2
364     & + (BAR(1,1)-(2.66*MX+BAR(1,1)))**2 + 2.66**2) /
365     & 3.6
366     C
367     ENDIF
368     C
369     C
370     C RIG IS RIGIDITY AS DETERMINED BY THE TRACKER
371     C OR by CALORIMETER IF IN SELFTRIGGER MODE
372     C
373     IF (GOOD2.EQ.1) THEN
374     GTR = 1
375     IF (TRKCHI2.LT.0) GTR = 2
376     IF ( AL_PP(GTR,5).NE.0. ) THEN
377     RIG = 1./(AL_PP(GTR,5))
378     ELSE
379     GOOD2 = 0
380     PRINT *,' CALORIMETER - WARNING F77: track with R = 0'
381     GOTO 50
382     ENDIF
383     ENDIF
384 mocchiut 1.7 IF (TRIGTY.GE.2.AND.HZN.EQ.0.AND.GOOD2.EQ.0) THEN
385 mocchiut 1.1 RIG = ELEN ! SELFTRIGGER RIGIDITY
386     IF ( RIG.EQ.0. ) THEN
387 mocchiut 1.8 GOOD2 = 1
388 mocchiut 1.1 PRINT *,' CALORIMETER - WARNING F77: ST track with R = 0'
389     GOTO 50
390     ENDIF
391     ENDIF
392     C
393 mocchiut 1.7 IF (GOOD2.EQ.0.AND.(TRIGTY.LT.2.OR.HZN.EQ.1)) THEN
394     RIG = RIGINPUT
395     ENDIF
396     C
397 mocchiut 1.1 RNSS = 0.
398     QTOTT = 0.
399     C
400 mocchiut 1.15 PPLANEMAX = 1.01*(LOG(ABS(RIG)/0.0081)-1.) / 0.74
401 mocchiut 1.1 C
402     IPLANE = INT(ANINT(PPLANEMAX)) + 5
403     C
404     IF (IPLANE.GT.NPLA) IPLANE=NPLA
405     IF (IPLANE.LT.1) IPLANE = 1
406 mocchiut 1.16 c print *,' calcolo...'
407 mocchiut 1.1 C
408     C CALCULATE QLOW AND NLOW
409     C
410     DO J = IPLANE,NPLA
411     DO I = 1,NCHA
412     IF (DEXY(1,J,I).GE.EMIN) THEN
413     NLOW = NLOW + 1
414     QLOW = QLOW + DEXY(1,J,I)
415     ENDIF
416     IF (DEXY(2,J,I).GE.EMIN) THEN
417     NLOW = NLOW + 1
418     QLOW = QLOW + DEXY(2,J,I)
419     ENDIF
420     ENDDO
421     ENDDO
422     C
423     C CALCULATE QCORE AND NCORE
424     C
425     C
426     C 8 STRIPS ARE 2.88 cm , A MOLIERE RADIUS IS ABOUT 0.7 cm .
427     C
428     DO J = 1,IPLANE
429     NNX = IBAR(1,J)
430 mocchiut 1.15 RNSS = 0. ! BACO!!
431     QTOTT = 0. ! BACO!!
432 mocchiut 1.1 IF (NNX.NE.-1) THEN
433     IF (NNX.LT.9) NNX = 9
434     IF (NNX.GT.88) NNX = 88
435     INFX = NNX - 8
436     ISUPX = NNX + 8
437     DO I = INFX,ISUPX
438     IF (DEXY(1,J,I).GE.EMIN) THEN
439     RNSS = RNSS + 1
440     QTOTT = QTOTT + DEXY(1,J,I)
441     ENDIF
442     ENDDO
443     ENDIF
444     C
445     NNY = IBAR(2,J)
446     IF (NNY.NE.-1) THEN
447     IF (NNY.LT.9) NNY = 9
448     IF (NNY.GT.88) NNY = 88
449     INFY = NNY - 8
450     ISUPY = NNY + 8
451     DO I = INFY,ISUPY
452     IF (DEXY(2,J,I).GE.EMIN) THEN
453     RNSS = RNSS + 1
454     QTOTT = QTOTT + DEXY(2,J,I)
455     ENDIF
456     ENDDO
457     ENDIF
458     NCORE = RNSS * FLOAT(J) + NCORE
459     QCORE = QTOTT * FLOAT(J) + QCORE
460     ENDDO
461     C
462     C CALCULATE NOINT
463     C
464     CALL NOINTER(NIN)
465     NOINT = FLOAT(NIN)
466     C
467     C
468     C QCYL = DETECTED ENERGY AND NCYL = NUMBER OF HIT STRIPS IN A CYLINDER oF
469     C RADIUS 8.5 STRIPS WITH AXIS DEFINED BY THE DIRECTION OF THE INCOMING
470     C PARTICLE .
471     C
472     C 8 STRIPS ARE 2.88 cm , A MOLIERE RADIUS IS ABOUT 0.7 cm .
473     C
474     DO J = 1,NPLA
475     NNX = IBAR(1,J)
476     IF (NNX.NE.-1) THEN
477     IF (NNX.LT.9) NNX = 9
478     IF (NNX.GT.88) NNX = 88
479     INFX = NNX - 8
480     ISUPX = NNX + 8
481     DO I = INFX,ISUPX
482     IF (DEXY(1,J,I).LT.EMIN) GO TO 710
483     NCYL = NCYL + 1
484     QCYL = QCYL + DEXY(1,J,I)
485     710 ENDDO
486     ENDIF
487     NNY = IBAR(2,J)
488     IF (NNY.NE.-1) THEN
489     IF (NNY.LT.9) NNY = 9
490     IF (NNY.GT.88) NNY = 88
491     INFY = NNY - 8
492     ISUPY = NNY + 8
493     DO I=INFY,ISUPY
494     IF (DEXY(2,J,I).LT.EMIN) GO TO 810
495     NCYL = NCYL + 1
496     QCYL = QCYL + DEXY(2,J,I)
497     810 ENDDO
498     ENDIF
499     C
500     C QTR = DETECTED ENERGY AND NTR = NUMBER OF HIT STRIPS IN A CYLINDER oF
501     C RADIUS 4.5 STRIPS WITH AXIS DEFINED BY THE DIRECTION OF THE INCOMING
502     C PARTICLE .
503     C
504     NNX = IBAR(1,J)
505     IF (NNX.NE.-1) THEN
506     IF (NNX.LT.5) NNX = 5
507     IF (NNX.GT.92) NNX = 92
508     INFX = NNX - 4
509     ISUPX = NNX + 4
510     DO I = INFX,ISUPX
511     IF (DEXY(1,J,I).GT.EMIN) THEN
512     NTR = NTR + 1
513     QTR = QTR + DEXY(1,J,I)
514     ENDIF
515     ENDDO
516     ENDIF
517     C
518     NNY = IBAR(2,J)
519     IF (NNY.NE.-1) THEN
520     IF (NNY.LT.5) NNY = 5
521     IF (NNY.GT.92) NNY = 92
522     INFY = NNY - 4
523     ISUPY = NNY + 4
524     DO I = INFY, ISUPY
525     IF (DEXY(2,J,I).GT.EMIN) THEN
526     NTR = NTR + 1
527     QTR = QTR + DEXY(2,J,I)
528     ENDIF
529     ENDDO
530     ENDIF
531     ENDDO
532     C
533     C CALCULATE QTRACK
534     C
535     CALL LATERALE(QTRACK,RQT)
536    
537     C
538     C CALCULATE NPRESH AND QPRESH
539     C
540     DO I = 1,4
541     NNX = IBAR(1,I)
542     IF (NNX.NE.-1) THEN
543     IF (NNX.LT.3) NNX = 3
544     IF (NNX.GT.94) NNX = 94
545     INFX = NNX - 2
546     ISUPX = NNX + 2
547     DO J = INFX,ISUPX
548     IF (DEXY(1,I,J).GE.EMIN) THEN
549     NPRESH = NPRESH + 1
550     QPRESH = QPRESH + DEXY(1,I,J)
551     ENDIF
552     ENDDO
553     ENDIF
554     C
555     NNY = IBAR(2,I)
556     IF (NNY.NE.-1) THEN
557     IF (NNY.LT.3) NNY = 3
558     IF (NNY.GT.94) NNY = 94
559     INFY = NNY - 2
560     ISUPY = NNY + 2
561     DO J = INFY,ISUPY
562     IF (DEXY(2,I,J).GE.EMIN) THEN
563     NPRESH = NPRESH + 1
564     QPRESH = QPRESH + DEXY(2,I,J)
565     ENDIF
566     ENDDO
567     ENDIF
568     ENDDO
569     C
570     C CALCULATE DXTRACK, DYTRACK, QTRACKX AND QTRACKY
571     C
572     ICONTROL5 = 0
573     CALL NSHOWER(ICONTROL5,DXTRACK,DYTRACK,QTRACKX,QTRACKY)
574     C
575     C CALCULATE QPRE AND NPRE
576     C
577     DO J = 1,3
578     NNX = IBAR(1,J)
579     IF (NNX.NE.-1) THEN
580     IF (NNX.LT.9) NNX = 9
581     IF (NNX.GT.88) NNX = 88
582     INFX = NNX - 8
583     ISUPX = NNX + 8
584     DO I = INFX,ISUPX
585     IF (DEXY(1,J,I).GE.EMIN) THEN
586     NPRE = NPRE + 1
587     QPRE = QPRE + DEXY(1,J,I)
588     ENDIF
589     ENDDO
590     ENDIF
591     C
592     NNY = IBAR(2,J)
593     IF (NNY.NE.-1) THEN
594     IF (NNY.LT.9) NNY = 9
595     IF (NNY.GT.88) NNY = 88
596     INFY = NNY - 8
597     ISUPY = NNY + 8
598     DO I=INFY,ISUPY
599     IF (DEXY(2,J,I).GE.EMIN) THEN
600     NPRE = NPRE + 1
601     QPRE = QPRE + DEXY(2,J,I)
602     ENDIF
603     ENDDO
604     ENDIF
605     ENDDO
606     C
607     C CALCULATE NLAST AND QLAST
608     C
609     DO J = NPLA-4,NPLA
610     NNX = IBAR(1,J)
611     IF (NNX.NE.-1) THEN
612 mocchiut 1.2 IF (NNX.LT.5) NNX = 5
613     IF (NNX.GT.92) NNX = 92
614     c IF (NNX.LT.9) NNX = 9
615     c IF (NNX.GT.88) NNX = 88
616     INFX = NNX - 4
617     ISUPX = NNX + 4
618     c INFX = NNX - 8
619     c ISUPX = NNX + 8
620 mocchiut 1.1 DO I = INFX,ISUPX
621     IF (DEXY(1,J,I).GE.EMIN) THEN
622     NLAST = NLAST + 1
623     QLAST = QLAST + DEXY(1,J,I)
624     ENDIF
625     ENDDO
626     ENDIF
627     C
628     NNY = IBAR(2,J)
629     IF (NNY.NE.-1) THEN
630 mocchiut 1.2 IF (NNY.LT.5) NNY = 5
631     IF (NNY.GT.92) NNY = 92
632     c IF (NNY.LT.9) NNY = 9
633     c IF (NNY.GT.88) NNY = 88
634     INFY = NNY - 4
635     ISUPY = NNY + 4
636     c INFY = NNY - 8
637     c ISUPY = NNY + 8
638 mocchiut 1.1 DO I=INFY,ISUPY
639     IF (DEXY(2,J,I).GE.EMIN) THEN
640     NLAST = NLAST + 1
641     QLAST = QLAST + DEXY(2,J,I)
642     ENDIF
643     ENDDO
644     ENDIF
645     ENDDO
646     C
647     C
648     C CALCULATE PLANETOT AND QMEAN
649     C
650     DO M = 1,2
651     RPIANO(M) = 0.
652     NTOT(M) = 0
653     ENDDO
654     NPIANI = 5
655     QMEAN = 0.
656     INDEX = 0
657 mocchiut 1.12 C
658     IF (TRIGTY.GE.2.AND.HZN.NE.0) THEN
659     EINF = 50.
660     ESUP = 15000.
661     CALL NUCLEI(RPIANO,NPIANI,QMEAN,NTOT,INDEX)
662     PLANETOT = RPIANO(1) + RPIANO(2)
663     ELSE
664     EINF = EMIN
665     ESUP = 15000.
666     CALL ELIO(RPIANO,NPIANI,QMEAN,NTOT,INDEX)
667     PLANETOT = RPIANO(1) + RPIANO(2)
668     ENDIF
669 mocchiut 1.1 C
670     50 CONTINUE
671     C
672 mocchiut 1.16 c print *,' esco'
673 mocchiut 1.1 RETURN
674     END
675    
676    
  ViewVC Help
Powered by ViewVC 1.1.23