/[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.5 by mocchiut, Thu Jun 29 12:56:46 2006 UTC revision 1.22 by mocchiut, Tue Aug 4 14:01:18 2009 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 21  C Line 21  C
21        REAL ax,bx,eax,ebx,chi2x        REAL ax,bx,eax,ebx,chi2x
22        REAL ay,by,eay,eby,chi2y        REAL ay,by,eay,eby,chi2y
23        REAL parzen3, TMISD        REAL parzen3, TMISD
24        INTEGER Nfitx,Nfity        INTEGER Nfitx,Nfity, MNPLA
25  C  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 FIRST TRACK IF WE HAVE A TRACK FROM TRACKER
96  C  C
97        T = 1        T = 1
98  C  C
99   10   CONTINUE  c 10   CONTINUE
100          CONTINUE
101  C  C
102        IF (GOOD2.EQ.1) THEN        IF (GOOD2.EQ.1) THEN
103  C      C    
104           CHTRACK = 0           CHTRACK = 0
105  C  C
106           CALL VZERO(IWPL,2)           CALL VZERO(IWPL,2)
107           CALL VZERO(BAR,2*NPLA)           CALL VZERO(BAR,2*NPLAV)
108           CALL VZERO(IBAR,2*NPLA)           CALL VZERO(IBAR,2*NPLAV)
109           CALL VZERO(TBAR,2*NPLA)           CALL VZERO(TBAR,2*NPLAV)
110           CALL VZERO(TIBAR,2*NPLA)           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           do m = 1, 5           do m = 1, 5
116              al_p(m) = al_pp(t,m)              al_p(m) = al_pp(t,m)
117    c            print *,' al_p(',m,') = ',al_p(m)
118           enddo           enddo
119           if (al_p(5).eq.0.) THEN           if (al_p(5).eq.0.) THEN
120         PRINT *,' CALORIMETER - WARNING F77: track with R = 0, discarded'         PRINT *,' CALORIMETER - WARNING F77: track with R = 0, discarded'
# Line 131  C Line 127  C
127                 YOUT(I) = 0.                 YOUT(I) = 0.
128                 IF (MOD(M,2).EQ.0) THEN                 IF (MOD(M,2).EQ.0) THEN
129                    DISTX = PIANO(I) + ZALIG                    DISTX = PIANO(I) + ZALIG
130    c                  print *,'T Y PLANE I= ',I,' Z = ',DISTX
131                 ELSE                 ELSE
132                    DISTX = PIANO(I) - 5.81 + ZALIG                    DISTX = PIANO(I) - 5.81 + ZALIG
133    c                  print *,'T X PLANE I= ',I,' Z = ',DISTX
134                 ENDIF                               ENDIF              
135                 ZIN(I) = distx / 10.                 ZIN(I) = distx / 10.
136    c               print *,' ZIN(',I,') = ',ZIN(I)
137                 TBAR(M,I) = 0.                 TBAR(M,I) = 0.
138                 TIBAR(M,I) = 0                 TIBAR(M,I) = 0
139              enddo              enddo
# Line 149  c               print *,' CALORIMETER - Line 148  c               print *,' CALORIMETER -
148              TY = TAN(ASIN(AL_P(3))) * SIN(AL_P(4))              TY = TAN(ASIN(AL_P(3))) * SIN(AL_P(4))
149              DO I = 1, NPLA              DO I = 1, NPLA
150                 NN = 0                 NN = 0
151                 IF (M.EQ.2) NN = 1  c               IF (M.EQ.2) NN = 1
152                 IF (MOD(I,2).EQ.NN) THEN                 IF (MOD(I,2).EQ.NN) THEN
153                    SHIFT = +0.5                    IF (REVERSE.EQ.0) THEN
154                         SHIFT = -0.5
155                      ELSE
156                         SHIFT = +0.5
157                      ENDIF
158                 ELSE                 ELSE
159                    SHIFT = -0.5                    IF (REVERSE.EQ.0) THEN
160                         SHIFT = +0.5
161                      ELSE
162                         SHIFT = -0.5
163                      ENDIF
164                 ENDIF                 ENDIF
165  C      C    
166  C     CHECK IF XOUT OR YOUT ARE NaN  C     CHECK IF XOUT OR YOUT ARE NaN
# Line 176  C     Line 183  C    
183                    Y(I) = CX                    Y(I) = CX
184                    BAR(M,I) = Y(I)                        BAR(M,I) = Y(I)    
185                    TBAR(M,I) = (Y(I) - XALIG)/10.                    TBAR(M,I) = (Y(I) - XALIG)/10.
186                    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))/
187         &                 ABS(ZIN(1)-ZIN(NPLA))
188                 ELSE                 ELSE
189                    YY(I) = CY                    YY(I) = CY
190                    BAR(M,I) = YY(I)                                      BAR(M,I) = YY(I)                  
191                    TBAR(M,I) = (-YALIG + YY(I))/10.                        TBAR(M,I) = (-YALIG + YY(I))/10.    
192                    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))/
193         &                 ABS(ZIN(1)-ZIN(NPLA))
194                 ENDIF                 ENDIF
195                 CALL LASTRISCIA(BAR(M,I),IBAR(M,I))                 CALL LASTRISCIA(BAR(M,I),IBAR(M,I))
196                 tibar(M,I) = ibar(m,i)                 tibar(M,I) = ibar(m,i)
# Line 235  C Line 244  C
244  C          C        
245  C     WE MUST PROCESS A SELFTRIGGER EVENT  C     WE MUST PROCESS A SELFTRIGGER EVENT
246  C  C
247        IF (TRIGTY.GE.2) THEN        IF (TRIGTY.GE.2.AND.HZN.EQ.0) THEN
248  C  C
249  C     CALL SELFTRIGGER SUBROUTINE  C     CALL SELFTRIGGER SUBROUTINE
250  C  C
251             CALL VZERO(IWPL,2)
252             CALL VZERO(VARCFIT,2)
253             CALL VZERO(NPCFIT,2)
254             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    C
263           CALL SELFTRIG()           CALL SELFTRIG()
264           ELEN = PARZEN3           ELEN = PARZEN3
265           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 269  C
269  C      C    
270           DO M = 1,2           DO M = 1,2
271  C  C
272    c            print *,' ax ',ax,' ay ',ay
273    c            print *,' bx ',bx,' by ',by
274              IF (NPCFIT(M).GE.2) THEN              IF (NPCFIT(M).GE.2) THEN
275                 IF (M.EQ.1) THEN                 IF (M.EQ.1) THEN
276                    VARCFIT(1) = CHI2X                    VARCFIT(1) = CHI2X
277                    IMPX = 10. * ( AX + 12.1 )                    IMPX = AX + BX * (ZALIG/10.) ! PAMELA REF
278                    TANX = BX                    TANX = BX
279                 ELSE                 ELSE
280                    VARCFIT(2) = CHI2Y                    VARCFIT(2) = CHI2Y
281                    IMPY = 10. * ( AY + 12.2 )                    IMPY = AY + BY * (ZALIG/10.) ! PAMELA REF
282                    TANY = BY                    TANY = BY
283                 ENDIF                 ENDIF
284  C  C
285                 DO I = 1,NPLA                     DO I = 1,NPLA    
286                    NN = 0                    NN = 0
287                    IF (M.EQ.2) NN = 1  c                  IF (M.EQ.2) NN = 1
288                    IF (MOD(I,2).EQ.NN) THEN                    IF (MOD(I,2).EQ.NN) THEN
289                       SHIFT = +0.5                       IF (REVERSE.EQ.0) THEN
290                            SHIFT = -0.5
291                         ELSE
292                            SHIFT = +0.5
293                         ENDIF
294                    ELSE                    ELSE
295                       SHIFT = -0.5                       IF (REVERSE.EQ.0) THEN
296                            SHIFT = +0.5
297                         ELSE
298                            SHIFT = -0.5
299                         ENDIF
300                    ENDIF                    ENDIF
301  C      C    
302                    IF (M.EQ.1) THEN                    IF (M.EQ.1) THEN
303                       DISTX = PIANO(I) - 5.81                       DISTX = PIANO(I) - 5.81
304                       Y(I) = DISTX * TANX + CX                       Y(I) = (DISTX * TANX) +  IMPX*10. !- XALIG
305                       BAR(M,I) = Y(I)  c                     CBAR(M,I) = Y(I)
306                       CBAR(M,I) = Y(I)                               BAR(M,I) =  Y(I) + XALIG
307                       IF (I.EQ.22) MX=ABS(Y(1)-Y(22))/ABS(ZIN(1)-ZIN(22))                       CBAR(M,I) = Y(I) / 10.
308                         IF (I.EQ.NPLA) MX=ABS(Y(1)-Y(NPLA))/
309         &                    ABS(ZIN(1)-ZIN(NPLA))
310  C      C    
311                    ELSE                    ELSE
312                       DISTY = PIANO(I)                                       DISTY = PIANO(I)                
313                       YY(I) = DISTY * TANY + CY                       YY(I) = (DISTY * TANY) + IMPY*10. !- YALIG
314                       BAR(M,I) = YY(I)  c                     print *,' I ',i,' YY ',YY(I)
315                       CBAR(M,I) = YY(I)  c                     CBAR(M,I) = YY(I)
316                       IF (I.EQ.22) MY=ABS(Y(1)-Y(22))/ABS(ZIN(1)-ZIN(22))                       BAR(M,I) =  YY(I) + YALIG
317                         CBAR(M,I) = YY(I) / 10.
318                         IF (I.EQ.NPLA) MY=ABS(Y(1)-Y(NPLA))/
319         &                    ABS(ZIN(1)-ZIN(NPLA))
320  C      C    
321                    ENDIF                    ENDIF
322                    CALL LASTRISCIA(BAR(M,I),IBAR(M,I))                    CALL LASTRISCIA(BAR(M,I),IBAR(M,I))
323                    cibar(M,I) = ibar(m,i)                    cibar(M,I) = ibar(m,i)
324                      IF (ibar(m,i).EQ.-1) THEN
325                         CHTRACK = CHTRACK + 1
326                      ELSE
327                         IWPL(M) = IWPL(M) + 1
328                      ENDIF
329                 ENDDO                             ENDDO            
330              ENDIF              ENDIF
331  C  C
332           ENDDO           ENDDO
333  C  C
334        ELSE        ENDIF
335           IF (GOOD2.EQ.0) THEN  C
336          IF (TRIGTY.GE.2.AND.HZN.NE.0) THEN
337             IF (GOOD2.EQ.1) THEN
338              PRINT *,' CALORIMETER - WARNING F77: unknown request'              PRINT *,' CALORIMETER - WARNING F77: unknown request'
339              GOOD2 = 1              GOOD2 = 1
340              GOTO 50              GOTO 50
341           ENDIF           ENDIF
342             IF ( NPCFIT(1).EQ.0.OR.NPCFIT(2).EQ.0 ) THEN
343                GOOD2 = 1
344                GOTO 50
345             ENDIF
346        ENDIF        ENDIF
347  C  C
348   6996 CONTINUE   6996 CONTINUE
# Line 304  C Line 351  C
351  C  C
352  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
353  C  C
354        IF (CHTRACK.EQ.44) THEN        IF (CHTRACK.EQ.44) THEN  ! CHTRACK is the number of planes not touched by the track
355           GOOD2 = 0           GOOD2 = 0
356  c         PRINT *,' CALORIMETER - WARNING F77: track outside calorimeter'  c         PRINT *,' CALORIMETER - WARNING F77: track outside calorimeter'
357           GOTO 50           GOTO 50
# Line 319  C Line 366  C
366       &      + (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) /
367       &      3.6         &      3.6  
368  C  C
 C         DX0L = X01PL * SQRT( (IWPL(1) * SQRT(1 + MX*MX))**2 +  
 C     &                        (IWPL(2) * SQRT(1 + MY*MY))**2 )/2.  
369        ENDIF        ENDIF
370  C  C
371  C  C
# Line 338  C Line 383  C
383              GOTO 50              GOTO 50
384           ENDIF           ENDIF
385        ENDIF        ENDIF
386        IF (TRIGTY.GE.2.AND.GOOD2.EQ.0) THEN        IF (TRIGTY.GE.2.AND.HZN.EQ.0.AND.GOOD2.EQ.0) THEN
387           RIG = ELEN ! SELFTRIGGER RIGIDITY           RIG = ELEN ! SELFTRIGGER RIGIDITY
388           IF ( RIG.EQ.0. ) THEN           IF ( RIG.EQ.0. ) THEN
389              GOOD2 = 0              GOOD2 = 1
390              PRINT *,' CALORIMETER - WARNING F77: ST track with R = 0'              PRINT *,' CALORIMETER - WARNING F77: ST track with R = 0'
391              GOTO 50              GOTO 50
392           ENDIF           ENDIF
393        ENDIF        ENDIF
394  C  C
395          IF (GOOD2.EQ.0.AND.(TRIGTY.LT.2.OR.HZN.EQ.1)) THEN
396             RIG = RIGINPUT
397          ENDIF
398    C
399        RNSS = 0.        RNSS = 0.
400        QTOTT = 0.        QTOTT = 0.
401  C  C
402        PPLANEMAX = 1.01*(LOG(ABS(RIG)/0.0081)-1.)        PPLANEMAX = 1.01*(LOG(ABS(RIG)/0.0081)-1.) / 0.74
403  C  C
404        IPLANE = INT(ANINT(PPLANEMAX)) + 5        IPLANE = INT(ANINT(PPLANEMAX)) + 5
405  C  C
406        IF (IPLANE.GT.NPLA) IPLANE=NPLA        IF (IPLANE.GT.NPLA) IPLANE=NPLA
407        IF (IPLANE.LT.1) IPLANE = 1        IF (IPLANE.LT.1) IPLANE = 1
408    c      print *,' calcolo...'
409  C  C
410  C     CALCULATE QLOW AND NLOW  C     CALCULATE QLOW AND NLOW
411  C  C
# Line 379  C     8 STRIPS ARE 2.88 cm , A MOLIERE R Line 429  C     8 STRIPS ARE 2.88 cm , A MOLIERE R
429  C      C    
430        DO J = 1,IPLANE        DO J = 1,IPLANE
431           NNX = IBAR(1,J)           NNX = IBAR(1,J)
432             RNSS = 0.               ! BACO!!
433             QTOTT = 0.               ! BACO!!
434           IF (NNX.NE.-1) THEN           IF (NNX.NE.-1) THEN
435              IF (NNX.LT.9) NNX = 9              IF (NNX.LT.9) NNX = 9
436              IF (NNX.GT.88) NNX = 88              IF (NNX.GT.88) NNX = 88
# Line 429  C     Line 481  C    
481              INFX = NNX - 8              INFX = NNX - 8
482              ISUPX = NNX + 8              ISUPX = NNX + 8
483              DO I = INFX,ISUPX              DO I = INFX,ISUPX
484                 IF (DEXY(1,J,I).LT.EMIN) GO TO 710                 IF (DEXY(1,J,I).GE.EMIN) THEN
485                 NCYL = NCYL + 1                    NCYL = NCYL + 1
486                 QCYL = QCYL + DEXY(1,J,I)                    QCYL = QCYL + DEXY(1,J,I)
487   710        ENDDO                 ENDIF
488                ENDDO
489           ENDIF           ENDIF
490           NNY = IBAR(2,J)           NNY = IBAR(2,J)
491           IF (NNY.NE.-1) THEN           IF (NNY.NE.-1) THEN
# Line 441  C     Line 494  C    
494              INFY = NNY - 8              INFY = NNY - 8
495              ISUPY = NNY + 8              ISUPY = NNY + 8
496              DO I=INFY,ISUPY              DO I=INFY,ISUPY
497                 IF (DEXY(2,J,I).LT.EMIN) GO TO 810                 IF (DEXY(2,J,I).GE.EMIN) THEN
498                 NCYL = NCYL + 1                    NCYL = NCYL + 1
499                 QCYL = QCYL + DEXY(2,J,I)                    QCYL = QCYL + DEXY(2,J,I)
500   810        ENDDO                 ENDIF
501                ENDDO
502           ENDIF           ENDIF
503  C      C    
504  C     QTR = DETECTED ENERGY AND NTR = NUMBER OF HIT STRIPS IN A CYLINDER oF  C     QTR = DETECTED ENERGY AND NTR = NUMBER OF HIT STRIPS IN A CYLINDER oF
# Line 556  C Line 610  C
610  C      C    
611  C     CALCULATE NLAST AND QLAST  C     CALCULATE NLAST AND QLAST
612  C  C
613        DO J = NPLA-4,NPLA        MNPLA = NPLA -4
614          IF ( MNPLA .LT. 1 ) MNPLA = 1
615          DO J = MNPLA,NPLA
616           NNX = IBAR(1,J)           NNX = IBAR(1,J)
617           IF (NNX.NE.-1) THEN           IF (NNX.NE.-1) THEN
618              IF (NNX.LT.5) NNX = 5              IF (NNX.LT.5) NNX = 5
# Line 594  c            ISUPY = NNY + 8 Line 650  c            ISUPY = NNY + 8
650           ENDIF           ENDIF
651        ENDDO        ENDDO
652  C  C
       EINF = EMIN  
       ESUP = 50.  
653  C  C
654  C     CALCULATE PLANETOT AND QMEAN  C     CALCULATE PLANETOT AND QMEAN
655  C  C
# Line 606  C Line 660  C
660        NPIANI = 5        NPIANI = 5
661        QMEAN = 0.        QMEAN = 0.
662        INDEX = 0        INDEX = 0
663        CALL ELIO(RPIANO,NPIANI,QMEAN,NTOT,INDEX)  C
664        PLANETOT = RPIANO(1) + RPIANO(2)          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  C  C
676   50   CONTINUE   50   CONTINUE
677  C  C
678    c      print *,' esco'
679        RETURN        RETURN
680        END        END
681    
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.22
  ViewVC Help
Powered by ViewVC 1.1.23