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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1 by mocchiut, Fri May 19 13:15:50 2006 UTC revision 1.16 by mocchiut, Fri Jul 20 08:24:53 2007 UTC
# Line 6  C Line 6  C
6        INCLUDE 'INTEST.TXT'        INCLUDE 'INTEST.TXT'
7  C  C
8        DOUBLE PRECISION al_p(5),        DOUBLE PRECISION al_p(5),
9       &     xout(npla),yout(npla),zin(npla)       &     xout(nplav),yout(nplav),zin(nplav)
10  C  C
11        REAL PIANO(22), VARFIT(2)        REAL PIANO(NPLAV), VARFIT(2)
12        REAL TX, TY, SHIFT        REAL TX, TY, SHIFT
13        REAL BAR(2,NPLA), DISTY        REAL BAR(2,NPLAV), DISTY
14        REAL DISTX, Y(NPLA), YY(NPLA)        REAL DISTX, Y(NPLAV), YY(NPLAV)
15        REAL RIG, PPLANEMAX, RMASS        REAL RIG, PPLANEMAX, RMASS
16        REAL RNSS, QTOTT, RQT, MX, MY        REAL RNSS, QTOTT, RQT, MX, MY
17        REAL CHECK, ENER, CX, CY        REAL CHECK, ENER, CX, CY
# Line 26  C Line 26  C
26        INTEGER INDEX, NTOT(2), NPIANI, GTR        INTEGER INDEX, NTOT(2), NPIANI, GTR
27        INTEGER j, m, i, IWPL(2), timpx, timpy, T, nn        INTEGER j, m, i, IWPL(2), timpx, timpy, T, nn
28        INTEGER IPLANE, NNX, NNY, INFX, INFY, ISUPX, ISUPY        INTEGER IPLANE, NNX, NNY, INFX, INFY, ISUPX, ISUPY
29        INTEGER IBAR(2,NPLA), NPFIT(2), CHTRACK,IWPLU        INTEGER IBAR(2,NPLAV), NPFIT(2), CHTRACK,IWPLU
30        INTEGER Iquest(100), ICONTROL5, nin, IFAIL        INTEGER Iquest(100), ICONTROL5, nin, IFAIL
31  C  C
32        PARAMETER (X01PL=0.74)        PARAMETER (X01PL=0.74)
# Line 55  C Line 55  C
55        COMMON / CH / CHECK        COMMON / CH / CHECK
56        SAVE / CH /        SAVE / CH /
57  C  C
58        COMMON / CALOFIT / VARFIT, NPFIT        COMMON / CALOFIT / VARFIT, NPFIT, IWPL,CHTRACK
59        SAVE / CALOFIT /        SAVE / CALOFIT /
60  C  C
61        COMMON / pawcd / hmemor        COMMON / pawcd / hmemor
# Line 66  C Line 66  C
66  C  C
67  C Begin !  C Begin !
68  C  C
69    c      print *,' sono qui'
70        CALOL2TR = 0;        CALOL2TR = 0;
71        NCORE = 0.        NCORE = 0.
72        QCORE = 0.        QCORE = 0.
# Line 89  C Line 90  C
90        NLAST = 0.        NLAST = 0.
91        PLANETOT = 0.        PLANETOT = 0.
92        QMEAN = 0.        QMEAN = 0.
93        SELFTRIGGER = 0  C      SELFTRIGGER = 0
       CALL VZERO(VARCFIT,2)  
       CALL VZERO(NPCFIT,2)  
       CALL VZERO(TBAR,2*NPLA)  
       CALL VZERO(TIBAR,2*NPLA)  
       CALL VZERO(BAR,2*NPLA)  
       CALL VZERO(IBAR,2*NPLA)  
       CALL VZERO(IBAR,2*NPLA)  
       CALL VZERO(Y,NPLA)  
       CALL VZERO(YY,NPLA)  
       CALL VZERO(XOUT,NPLA)  
       CALL VZERO(YOUT,NPLA)  
94  C  C
95  C     BEGIN WITH THE FISRT TRACK IF WE HAVE A TRACK FROM TRACKER  C     BEGIN WITH THE FISRT TRACK IF WE HAVE A TRACK FROM TRACKER
96  C  C
# Line 113  C     Line 103  C    
103           CHTRACK = 0           CHTRACK = 0
104  C  C
105           CALL VZERO(IWPL,2)           CALL VZERO(IWPL,2)
106           CALL VZERO(BAR,2*NPLA)           CALL VZERO(BAR,2*NPLAV)
107           CALL VZERO(IBAR,2*NPLA)           CALL VZERO(IBAR,2*NPLAV)
108           CALL VZERO(TBAR,2*NPLA)           CALL VZERO(TBAR,2*NPLAV)
109           CALL VZERO(TIBAR,2*NPLA)           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           do m = 1, 5           do m = 1, 5
115              al_p(m) = al_pp(t,m)              al_p(m) = al_pp(t,m)
116    c            print *,' al_p(',m,') = ',al_p(m)
117           enddo           enddo
118           if (al_p(5).eq.0.) THEN           if (al_p(5).eq.0.) THEN
119         PRINT *,' CALORIMETER - WARNING F77: track with R = 0, discarded'         PRINT *,' CALORIMETER - WARNING F77: track with R = 0, discarded'
# Line 135  C Line 130  C
130                    DISTX = PIANO(I) - 5.81 + ZALIG                    DISTX = PIANO(I) - 5.81 + ZALIG
131                 ENDIF                               ENDIF              
132                 ZIN(I) = distx / 10.                 ZIN(I) = distx / 10.
133    c               print *,' ZIN(',I,') = ',ZIN(I)
134                 TBAR(M,I) = 0.                 TBAR(M,I) = 0.
135                 TIBAR(M,I) = 0                 TIBAR(M,I) = 0
136              enddo              enddo
# Line 142  C Line 138  C
138              call DOTRACK(NPLA,ZIN,XOUT,YOUT,AL_P,IFAIL)              call DOTRACK(NPLA,ZIN,XOUT,YOUT,AL_P,IFAIL)
139              if(IFAIL.ne.0)then              if(IFAIL.ne.0)then
140                 GOOD2 = 0                 GOOD2 = 0
141                 print *,' CALORIMETER - WARNING F77: tracking failed '  c               print *,' CALORIMETER - WARNING F77: tracking failed '
142                 goto 969                 goto 969
143              endif              endif
144              TX = TAN(ASIN(AL_P(3))) * COS(AL_P(4))              TX = TAN(ASIN(AL_P(3))) * COS(AL_P(4))
# Line 151  C Line 147  C
147                 NN = 0                 NN = 0
148                 IF (M.EQ.2) NN = 1                 IF (M.EQ.2) NN = 1
149                 IF (MOD(I,2).EQ.NN) THEN                 IF (MOD(I,2).EQ.NN) THEN
150                    SHIFT = +0.5                    IF (REVERSE.EQ.0) THEN
151                         SHIFT = +0.5
152                      ELSE
153                         SHIFT = -0.5
154                      ENDIF
155                 ELSE                 ELSE
156                    SHIFT = -0.5                    IF (REVERSE.EQ.0) THEN
157                         SHIFT = -0.5
158                      ELSE
159                         SHIFT = +0.5
160                      ENDIF
161                 ENDIF                 ENDIF
162  C      C    
163  C     CHECK IF XOUT OR YOUT ARE NaN  C     CHECK IF XOUT OR YOUT ARE NaN
164  C      C    
165                 IF (XOUT(I).NE.XOUT(I).OR.YOUT(I).NE.YOUT(I)) THEN                 IF (XOUT(I).NE.XOUT(I).OR.YOUT(I).NE.YOUT(I)) THEN
166                    print *,  c                  print *,
167       &         ' CALORIMETER - WARNING F77: tracking error (NaN values)'                    c     &         ' CALORIMETER - WARNING F77: tracking error (NaN values)'                  
168                    GOOD2 = 0                    GOOD2 = 0
169                    GOTO 969                    GOTO 969
170                 ENDIF                 ENDIF
171  C  C
172                 CX = XOUT(I)*10. + XALIG                 CX = XOUT(I)*10. + XALIG
173                 CY = -YOUT(I)*10. + YALIG                 CY = YOUT(I)*10. + YALIG
174  C      C    
175                 IF (I.EQ.1) THEN                 IF (I.EQ.1) THEN
176                    TIMPX = CX                    TIMPX = CX
# Line 176  C     Line 180  C    
180                    Y(I) = CX                    Y(I) = CX
181                    BAR(M,I) = Y(I)                        BAR(M,I) = Y(I)    
182                    TBAR(M,I) = (Y(I) - XALIG)/10.                    TBAR(M,I) = (Y(I) - XALIG)/10.
183                    IF (I.EQ.22) MX=ABS(Y(1)-Y(22))/ABS(ZIN(1)-ZIN(22))                    IF (I.EQ.NPLA) MX=ABS(Y(1)-Y(NPLA))/
184         &                 ABS(ZIN(1)-ZIN(NPLA))
185                 ELSE                 ELSE
186                    YY(I) = CY                    YY(I) = CY
187                    BAR(M,I) = YY(I)                                      BAR(M,I) = YY(I)                  
188                    TBAR(M,I) = (YALIG - YY(I))/10.                        TBAR(M,I) = (-YALIG + YY(I))/10.    
189                    IF (I.EQ.22) MY=ABS(Y(1)-Y(22))/ABS(ZIN(1)-ZIN(22))                    IF (I.EQ.NPLA) MY=ABS(Y(1)-Y(NPLA))/
190         &                 ABS(ZIN(1)-ZIN(NPLA))
191                 ENDIF                 ENDIF
192                 CALL LASTRISCIA(BAR(M,I),IBAR(M,I))                 CALL LASTRISCIA(BAR(M,I),IBAR(M,I))
193                 tibar(M,I) = ibar(m,i)                 tibar(M,I) = ibar(m,i)
# Line 235  C Line 241  C
241  C          C        
242  C     WE MUST PROCESS A SELFTRIGGER EVENT  C     WE MUST PROCESS A SELFTRIGGER EVENT
243  C  C
244        IF (TRIGTY.GE.2) THEN        IF (TRIGTY.GE.2.AND.HZN.EQ.0) THEN
245  C  C
246  C     CALL SELFTRIGGER SUBROUTINE  C     CALL SELFTRIGGER SUBROUTINE
247  C  C
248             CALL VZERO(IWPL,2)
249             CALL VZERO(VARCFIT,2)
250             CALL VZERO(NPCFIT,2)
251             CALL VZERO(TBAR,2*NPLAV)
252             CALL VZERO(TIBAR,2*NPLAV)
253             CALL VZERO(BAR,2*NPLAV)
254             CALL VZERO(IBAR,2*NPLAV)
255             CALL VZERO(Y,NPLAV)
256             CALL VZERO(YY,NPLAV)
257             CALL VZERO(XOUT,NPLAV)
258             CALL VZERO(YOUT,NPLAV)
259    C
260           CALL SELFTRIG()           CALL SELFTRIG()
261           ELEN = PARZEN3           ELEN = PARZEN3
262           SELEN = ABS(ELEN) * (11.98*1E-2 + 7.6 * EXP(-5736/ABS(ELEN)))           SELEN = ABS(ELEN) * (11.98*1E-2 + 7.6 * EXP(-5736/ABS(ELEN)))
# Line 248  C Line 266  C
266  C      C    
267           DO M = 1,2           DO M = 1,2
268  C  C
269    c            print *,' ax ',ax,' ay ',ay
270    c            print *,' bx ',bx,' by ',by
271              IF (NPCFIT(M).GE.2) THEN              IF (NPCFIT(M).GE.2) THEN
272                 IF (M.EQ.1) THEN                 IF (M.EQ.1) THEN
273                    VARCFIT(1) = CHI2X                    VARCFIT(1) = CHI2X
274                    IMPX = 10. * ( AX + 12.1 )                    IMPX = AX ! PAMELA REF
275                    TANX = BX                    TANX = BX
276                 ELSE                 ELSE
277                    VARCFIT(2) = CHI2Y                    VARCFIT(2) = CHI2Y
278                    IMPY = 10. * ( AY + 12.2 )                    IMPY = AY ! PAMELA REF
279                    TANY = BY                    TANY = BY
280                 ENDIF                 ENDIF
281  C  C
# Line 263  C Line 283  C
283                    NN = 0                    NN = 0
284                    IF (M.EQ.2) NN = 1                    IF (M.EQ.2) NN = 1
285                    IF (MOD(I,2).EQ.NN) THEN                    IF (MOD(I,2).EQ.NN) THEN
286                       SHIFT = +0.5                       IF (REVERSE.EQ.0) THEN
287                            SHIFT = +0.5
288                         ELSE
289                            SHIFT = -0.5
290                         ENDIF
291                    ELSE                    ELSE
292                       SHIFT = -0.5                       IF (REVERSE.EQ.0) THEN
293                            SHIFT = -0.5
294                         ELSE
295                            SHIFT = +0.5
296                         ENDIF
297                    ENDIF                    ENDIF
298  C      C    
299                    IF (M.EQ.1) THEN                    IF (M.EQ.1) THEN
300                       DISTX = PIANO(I) - 5.81                       DISTX = PIANO(I) - 5.81
301                       Y(I) = DISTX * TANX + CX                       Y(I) = (DISTX * TANX) +  AX - XALIG
302    c                     CBAR(M,I) = Y(I)
303                       BAR(M,I) = Y(I)                       BAR(M,I) = Y(I)
304                       CBAR(M,I) = Y(I)                               CBAR(M,I) = (Y(I) + XALIG)/10.
305                       IF (I.EQ.22) MX=ABS(Y(1)-Y(22))/ABS(ZIN(1)-ZIN(22))                       IF (I.EQ.NPLA) MX=ABS(Y(1)-Y(NPLA))/
306         &                    ABS(ZIN(1)-ZIN(NPLA))
307  C      C    
308                    ELSE                    ELSE
309                       DISTY = PIANO(I)                                       DISTY = PIANO(I)                
310                       YY(I) = DISTY * TANY + CY                       YY(I) = (DISTY * TANY) + AY - YALIG
311    c                     CBAR(M,I) = YY(I)
312                       BAR(M,I) = YY(I)                       BAR(M,I) = YY(I)
313                       CBAR(M,I) = YY(I)                       CBAR(M,I) = (YY(I) + YALIG)/10.
314                       IF (I.EQ.22) MY=ABS(Y(1)-Y(22))/ABS(ZIN(1)-ZIN(22))                       IF (I.EQ.NPLA) MY=ABS(Y(1)-Y(NPLA))/
315         &                    ABS(ZIN(1)-ZIN(NPLA))
316  C      C    
317                    ENDIF                    ENDIF
318                    CALL LASTRISCIA(BAR(M,I),IBAR(M,I))                    CALL LASTRISCIA(BAR(M,I),IBAR(M,I))
319                    cibar(M,I) = ibar(m,i)                    cibar(M,I) = ibar(m,i)
320                      IF (ibar(m,i).EQ.-1) THEN
321                         CHTRACK = CHTRACK + 1
322                      ELSE
323                         IWPL(M) = IWPL(M) + 1
324                      ENDIF
325                 ENDDO                             ENDDO            
326              ENDIF              ENDIF
327  C  C
328           ENDDO           ENDDO
329  C  C
330        ELSE        ENDIF
331           IF (GOOD2.EQ.0) THEN  C
332          IF (TRIGTY.GE.2.AND.HZN.NE.0) THEN
333             IF (GOOD2.EQ.1) THEN
334              PRINT *,' CALORIMETER - WARNING F77: unknown request'              PRINT *,' CALORIMETER - WARNING F77: unknown request'
335              GOOD2 = 1              GOOD2 = 1
336              GOTO 50              GOTO 50
337           ENDIF           ENDIF
338             IF ( NPCFIT(1).EQ.0.OR.NPCFIT(2).EQ.0 ) THEN
339                GOOD2 = 1
340                GOTO 50
341             ENDIF
342        ENDIF        ENDIF
343  C  C
344   6996 CONTINUE   6996 CONTINUE
# Line 304  C Line 347  C
347  C  C
348  C IF THE TRACK IS OUTSIDE THE CALORIMETER GO OUT, IF NOT CALCULATE DX0L  C IF THE TRACK IS OUTSIDE THE CALORIMETER GO OUT, IF NOT CALCULATE DX0L
349  C  C
350        IF (CHTRACK.EQ.44) THEN        IF (CHTRACK.EQ.44) THEN  ! CHTRACK is the number of planes not touched by the track
351           GOOD2 = 0           GOOD2 = 0
352           PRINT *,' CALORIMETER - WARNING F77: track outside calorimeter'  c         PRINT *,' CALORIMETER - WARNING F77: track outside calorimeter'
353           GOTO 50           GOTO 50
354        ELSE        ELSE
355           IF ( IWPL(1).LE.IWPL(2) ) THEN           IF ( IWPL(1).LE.IWPL(2) ) THEN
# Line 319  C Line 362  C
362       &      + (BAR(1,1)-(2.66*MX+BAR(1,1)))**2 + 2.66**2) /       &      + (BAR(1,1)-(2.66*MX+BAR(1,1)))**2 + 2.66**2) /
363       &      3.6         &      3.6  
364  C  C
 C         DX0L = X01PL * SQRT( (IWPL(1) * SQRT(1 + MX*MX))**2 +  
 C     &                        (IWPL(2) * SQRT(1 + MY*MY))**2 )/2.  
365        ENDIF        ENDIF
366  C  C
367  C  C
# Line 338  C Line 379  C
379              GOTO 50              GOTO 50
380           ENDIF           ENDIF
381        ENDIF        ENDIF
382        IF (TRIGTY.GE.2.AND.GOOD2.EQ.0) THEN        IF (TRIGTY.GE.2.AND.HZN.EQ.0.AND.GOOD2.EQ.0) THEN
383           RIG = ELEN ! SELFTRIGGER RIGIDITY           RIG = ELEN ! SELFTRIGGER RIGIDITY
384           IF ( RIG.EQ.0. ) THEN           IF ( RIG.EQ.0. ) THEN
385              GOOD2 = 0              GOOD2 = 1
386              PRINT *,' CALORIMETER - WARNING F77: ST track with R = 0'              PRINT *,' CALORIMETER - WARNING F77: ST track with R = 0'
387              GOTO 50              GOTO 50
388           ENDIF           ENDIF
389        ENDIF        ENDIF
390  C  C
391          IF (GOOD2.EQ.0.AND.(TRIGTY.LT.2.OR.HZN.EQ.1)) THEN
392             RIG = RIGINPUT
393          ENDIF
394    C
395        RNSS = 0.        RNSS = 0.
396        QTOTT = 0.        QTOTT = 0.
397  C  C
398        PPLANEMAX = 1.01*(LOG(ABS(RIG)/0.0081)-1.)        PPLANEMAX = 1.01*(LOG(ABS(RIG)/0.0081)-1.) / 0.74
399  C  C
400        IPLANE = INT(ANINT(PPLANEMAX)) + 5        IPLANE = INT(ANINT(PPLANEMAX)) + 5
401  C  C
402        IF (IPLANE.GT.NPLA) IPLANE=NPLA        IF (IPLANE.GT.NPLA) IPLANE=NPLA
403        IF (IPLANE.LT.1) IPLANE = 1        IF (IPLANE.LT.1) IPLANE = 1
404    c      print *,' calcolo...'
405  C  C
406  C     CALCULATE QLOW AND NLOW  C     CALCULATE QLOW AND NLOW
407  C  C
# Line 379  C     8 STRIPS ARE 2.88 cm , A MOLIERE R Line 425  C     8 STRIPS ARE 2.88 cm , A MOLIERE R
425  C      C    
426        DO J = 1,IPLANE        DO J = 1,IPLANE
427           NNX = IBAR(1,J)           NNX = IBAR(1,J)
428             RNSS = 0.               ! BACO!!
429             QTOTT = 0.               ! BACO!!
430           IF (NNX.NE.-1) THEN           IF (NNX.NE.-1) THEN
431              IF (NNX.LT.9) NNX = 9              IF (NNX.LT.9) NNX = 9
432              IF (NNX.GT.88) NNX = 88              IF (NNX.GT.88) NNX = 88
# Line 559  C Line 607  C
607        DO J = NPLA-4,NPLA        DO J = NPLA-4,NPLA
608           NNX = IBAR(1,J)           NNX = IBAR(1,J)
609           IF (NNX.NE.-1) THEN           IF (NNX.NE.-1) THEN
610              IF (NNX.LT.9) NNX = 9              IF (NNX.LT.5) NNX = 5
611              IF (NNX.GT.88) NNX = 88              IF (NNX.GT.92) NNX = 92
612              INFX = NNX - 8  c            IF (NNX.LT.9) NNX = 9
613              ISUPX = NNX + 8  c            IF (NNX.GT.88) NNX = 88
614                INFX = NNX - 4
615                ISUPX = NNX + 4
616    c            INFX = NNX - 8
617    c            ISUPX = NNX + 8
618              DO I = INFX,ISUPX              DO I = INFX,ISUPX
619                 IF (DEXY(1,J,I).GE.EMIN) THEN                 IF (DEXY(1,J,I).GE.EMIN) THEN
620                    NLAST = NLAST + 1                    NLAST = NLAST + 1
# Line 573  C Line 625  C
625  C  C
626           NNY = IBAR(2,J)           NNY = IBAR(2,J)
627           IF (NNY.NE.-1) THEN           IF (NNY.NE.-1) THEN
628              IF (NNY.LT.9) NNY = 9              IF (NNY.LT.5) NNY = 5
629              IF (NNY.GT.88) NNY = 88              IF (NNY.GT.92) NNY = 92
630              INFY = NNY - 8  c            IF (NNY.LT.9) NNY = 9
631              ISUPY = NNY + 8  c            IF (NNY.GT.88) NNY = 88
632                INFY = NNY - 4
633                ISUPY = NNY + 4
634    c            INFY = NNY - 8
635    c            ISUPY = NNY + 8
636              DO I=INFY,ISUPY              DO I=INFY,ISUPY
637                 IF (DEXY(2,J,I).GE.EMIN) THEN                 IF (DEXY(2,J,I).GE.EMIN) THEN
638                    NLAST = NLAST + 1                    NLAST = NLAST + 1
# Line 586  C Line 642  C
642           ENDIF           ENDIF
643        ENDDO        ENDDO
644  C  C
       EINF = EMIN  
       ESUP = 50.  
645  C  C
646  C     CALCULATE PLANETOT AND QMEAN  C     CALCULATE PLANETOT AND QMEAN
647  C  C
# Line 598  C Line 652  C
652        NPIANI = 5        NPIANI = 5
653        QMEAN = 0.        QMEAN = 0.
654        INDEX = 0        INDEX = 0
655        CALL ELIO(RPIANO,NPIANI,QMEAN,NTOT,INDEX)  C
656        PLANETOT = RPIANO(1) + RPIANO(2)          IF (TRIGTY.GE.2.AND.HZN.NE.0) THEN
657             EINF = 50.
658             ESUP = 15000.
659             CALL NUCLEI(RPIANO,NPIANI,QMEAN,NTOT,INDEX)
660             PLANETOT = RPIANO(1) + RPIANO(2)  
661          ELSE
662             EINF = EMIN
663             ESUP = 15000.
664             CALL ELIO(RPIANO,NPIANI,QMEAN,NTOT,INDEX)
665             PLANETOT = RPIANO(1) + RPIANO(2)  
666          ENDIF
667  C  C
668   50   CONTINUE   50   CONTINUE
669  C  C
670    c      print *,' esco'
671        RETURN        RETURN
672        END        END
673    
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.16
  ViewVC Help
Powered by ViewVC 1.1.23