/[PAMELA software]/calo/unpacking/nnplot_nt5.for
ViewVC logotype

Annotation of /calo/unpacking/nnplot_nt5.for

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1.1.1 - (hide annotations) (download) (vendor branch)
Mon Dec 5 16:23:21 2005 UTC (19 years ago) by mocchiut
Branch: MAIN, unpacking
CVS Tags: start, v1r00, HEAD
Changes since 1.1: +0 -0 lines
Imported sources

1 mocchiut 1.1 ***********************************************************************
2     *
3     * Subroutine READVME
4     * Routine to read a single event: use C subroutine readevnt.
5     * Use conditional compilation with cpp preprocessor, so this file must
6     * end in .F and cannot be included with INCLUDE directive (use #include)
7     * Input: IEV event # to read
8     * Output: fill ADC and DSP variables
9     *
10     * Rewritten from a previous routine from an HP machine
11     * S. Piccardi Nov. 1998
12     *
13     * $Id: readvme.F,v 1.1 1999/12/15 10:05:56 piccardi Exp $
14     *
15     *****************************************************************************
16     PROGRAM NNPLOT_NT5
17     c
18     IMPLICIT NONE
19     C
20     INCLUDE 'INTEST.TXT'
21     C
22     CHARACTER*6 file_name, file_name2
23     CHARACTER*40 file_name3
24     C
25     integer icount,icount2
26     c
27     integer RUNERROR, PEDERROR
28     C
29     C Normal variables definition
30     C
31     INTEGER FFD
32     C
33     INTEGER l, j, kk, m, mm, ll, il, ii, nn, ival, iev
34     INTEGER i, istat, icycle, ierr, icheck, icheck2, icheck3
35     INTEGER LVMIN
36    
37     INTEGER*2 length
38    
39     integer IP
40    
41     INTEGER ipre, ipre2, jc
42     INTEGER ic, k, jj
43     INTEGER status
44     INTEGER inf, isup, idet
45    
46     INTEGER LEVEL
47    
48     REAL ADCERR
49    
50     C$$$ INTEGER*4 lenght(4)
51    
52     INTEGER lundata, iosop
53    
54     INTEGER*2 st, st2
55    
56     CHARACTER*2 nr
57     c
58     integer*4 evnum, systime, len
59     integer nev
60     C
61     integer nset,nword(10)
62     C
63     INTEGER*2 STRINGA(200000)
64     C
65     REAL auto(11)
66    
67     integer lun(4)
68    
69     REAL SUM(6), SUM2(6), COUNT(6),
70     & BAS(4,11,6), DIS(6), DDIS(6), BASE(2)
71     REAL SSUM, SSUM2, CCOUNT, MME(2), SSIG
72    
73     REAL VAL
74    
75     REAL PED(4,11,96), GOOD(4,11,96)
76     COMMON / PEDE / PED, GOOD
77    
78     REAL MIP1(11,96), MIP2(11,96), MIP3(11,96), MIP4(11,96)
79     REAL SPOS1(3,11), SPOS2(3,11), SPOS3(3,11), SPOS4(3,11), DET
80    
81     REAL D1(11,96), D2(11,96), D3(11,96), D4(11,96)
82    
83     REAL DI, DS
84     REAL ME(6), SIG(6)
85    
86     REAL TG(2)
87     REAL SHIFT
88     COMMON /SHIFT/ SHIFT
89    
90     REAL BAR(2,NPLA)
91     INTEGER IBAR(2,NPLA)
92     COMMON/ANGOLO/BAR,IBAR
93    
94     REAL DISTX, DISTY, Y(NPLA), YY(NPLA)
95    
96     REAL CX, CY
97     COMMON/WHERE/CX,CY
98    
99     REAL RIG, PLANEMAX, RMASS
100     COMMON/GENERAL/RIG,RMASS
101    
102     INTEGER IPLANE, NNX, NNY, INFX, INFY, ISUPX, ISUPY
103     INTEGER NLK
104    
105     REAL RNSS, QTOTT, RQT
106    
107     real dedx1(11,96),dedx2(11,96),dedx3(11,96),dedx4(11,96)
108     real qpl(96), ene(96), qst(96)
109    
110     INTEGER nin
111    
112     REAL CHECK
113     COMMON / CH / CHECK
114    
115     REAL QMAXX, QMAXY
116     REAL VMAX
117    
118     real nstrip, qtot, ncore, qcore, impx, impy, tanx, tany, nint
119     REAL ncyl, qcyl, qtrack, qmax, cher, lead, nx22, qx22, qq(4)
120     INTEGER event
121     COMMON / Planes / nstrip, qtot, ncore, qcore, impx, impy,
122     & tanx, tany, nint, ncyl, qcyl, qtrack, qmax, cher, lead, nx22,
123     & qx22, qq, event
124    
125     INTEGER lrec
126     PARAMETER (lrec=6144)
127    
128     REAL hmemor(3000000)
129     COMMON /pawc/hmemor
130     CALL HLIMIT(3000000)
131    
132     C
133     C Begin !
134     C
135     RUNERROR = 0 ! error variable
136     PEDERROR = 0
137    
138     PRINT *,'Data file to read (only data and number included)? '
139     READ(*,905)file_name
140    
141     CALL PEDESTAL(PEDERROR)
142     IF (PEDERROR.EQ.1) THEN
143     PRINT *,'Problems with pedestal file'
144     GOTO 50
145     ENDIF
146     C
147     905 FORMAT(A6)
148    
149     lundata=33
150    
151     OPEN(UNIT=lundata,FILE='/wizard3/pamela/tb-0903/data/output_0309'
152     + //file_name//'.dat',STATUS='OLD',
153     + FORM='UNFORMATTED',ERR=50,IOSTAT=IOSOP)
154     C
155     C
156     C
157     PRINT *,'File to save? '
158     READ(*,904)file_name3
159    
160     904 FORMAT(A40)
161     C
162     C Histos creation
163     C
164     CALL HROPEN(59,'Event','/wizard3/pamela/tb-0903/paw/'//
165     + file_name3,'n',lrec,istat)
166     CALL HBNT(1,'Pamela Calo',' ')
167     CALL HBSET('BSIZE',lrec,ierr)
168    
169     *** /* Book ntuple variables */
170     CALL HBNAME(1,'calo',nstrip,'nstrip:R,qtot:R,
171     & ncore:R,qcore:R,impx:R,impy:R,tanx:R,tany:R,
172     & nint:R,ncyl:R,qcyl:R,qtrack:R,qmax:R,
173     & cher:R,lead:R,nx22:R,qx22:R,qq(4):R,event:I')
174     C
175     906 FORMAT(A2)
176     C
177     OPEN(71,FILE='/wizard3/pamela/tb-0903/source/mip1.dat',
178     & STATUS='OLD',FORM='FORMATTED',ERR=50,IOSTAT=IOSOP)
179     OPEN(72,FILE='/wizard3/pamela/tb-0903/source/mip2.dat',
180     & STATUS='OLD',FORM='FORMATTED',ERR=50,IOSTAT=IOSOP)
181     OPEN(73,FILE='/wizard3/pamela/tb-0903/source/nmip3.dat',
182     & STATUS='OLD',FORM='FORMATTED',ERR=50,IOSTAT=IOSOP)
183     OPEN(74,FILE='/wizard3/pamela/tb-0903/source/mip4.dat',
184     & STATUS='OLD',FORM='FORMATTED',ERR=50,IOSTAT=IOSOP)
185     C
186     OPEN(81,FILE='/wizard3/pamela/tb-0903/source/spos1.dat',
187     & STATUS='OLD',FORM='FORMATTED',ERR=50,IOSTAT=IOSOP)
188     OPEN(82,FILE='/wizard3/pamela/tb-0903/source/spos4.dat',
189     & STATUS='OLD',FORM='FORMATTED',ERR=50,IOSTAT=IOSOP)
190     OPEN(83,FILE='/wizard3/pamela/tb-0903/source/spos3.dat',
191     & STATUS='OLD',FORM='FORMATTED',ERR=50,IOSTAT=IOSOP)
192     OPEN(84,FILE='/wizard3/pamela/tb-0903/source/spos2.dat',
193     & STATUS='OLD',FORM='FORMATTED',ERR=50,IOSTAT=IOSOP)
194     C
195     DO I = 1,11
196     READ (71,102)(mip1(I,J),J=1,96)
197     READ (72,102)(mip2(I,J),J=1,96)
198     READ (73,102)(mip3(I,J),J=1,96)
199     READ (74,102)(mip4(I,J),J=1,96)
200     ENDDO
201     DO I = 1,3
202     READ (81,103)(SPOS1(I,J),J=1,11)
203     READ (82,103)(SPOS2(I,J),J=1,11)
204     READ (83,103)(SPOS3(I,J),J=1,11)
205     READ (84,103)(SPOS4(I,J),J=1,11)
206     print *,'I :',I,' J : 2 spos1 :',SPOS1(I,2)
207     ENDDO
208     CLOSE(21)
209     CLOSE(71)
210     CLOSE(72)
211     CLOSE(73)
212     CLOSE(74)
213     CLOSE(81)
214     CLOSE(82)
215     CLOSE(83)
216     CLOSE(84)
217    
218     102 FORMAT (4(2X,F9.3))
219     103 FORMAT (11(2X,F7.5))
220     C
221     C Take the file descriptor number (unix)
222     C
223     PRINT *,'Momentum:'
224     READ *,RIG
225     C
226     PRINT *,'Mass of the particle:'
227     READ *,RMASS
228     C
229     PRINT *,'How many events to read:'
230     READ *,nev
231     C
232     ffd=FNum(lundata) ! take the file descriptor number
233     C
234    
235     C Read event: use C readevent routine
236     C
237     icount = 0
238     icount2 = 0
239     event = 0
240     iev = 0.
241     CALL VZERO(BAS,66*4)
242     c
243     10 continue
244     C
245     event = event + 1
246     ival = 0
247     DO I = 1,11
248     DO J = 1,96
249     DEDX1(I,J) = 0.
250     DEDX2(I,J) = 0.
251     DEDX3(I,J) = 0.
252     DEDX4(I,J) = 0.
253     ENDDO
254     ENDDO
255     C
256     call readev(evnum,systime,len,stringa,RUNERROR,ffd) ! XXX
257    
258     if (runerror.eq.-1.or.runerror.eq.1) goto 50
259     C
260     c write (*,31) evnum
261     C write (*,32) systime
262     C write (*,33) len
263     C print *,'event number:',evnum
264     C print *,'systime:',systime
265     C print *,'event length:',len
266     31 format(1x,'event number :',1x,Z8)
267     32 format(1x,'systime :',1x,Z8)
268     33 format(1x,'event length :',1x,Z8)
269    
270     C$$$ write(6,*)nset,nword
271    
272     C$$ print *,'lenght 1 :',lenght(1)
273     C$$ print *,'lenght 2 :',lenght(2)
274     C$$ print *,'lenght 3 :',lenght(3)
275    
276     if (runerror.eq.-1.or.runerror.eq.1) goto 50
277    
278     ic = 0
279    
280     ic = ic + 4
281     length = ic
282     do kk = 1,4
283     C if (kk.eq.3) goto 45
284     ic = length + 1
285     C$$$ do k = ic-5,ic+5
286     k = ic
287     C write (*,41) stringa(k)
288     st2 = IAND(stringa(k),'FF00'x)
289     status = ISHFT(st2,-8)
290     C write (*,41) status
291     IF (KK.eq.1.and.status.ne.234) then
292     print *,'Problems with view:',KK
293     GOTO 45
294     endif
295     IF (KK.eq.2.and.status.ne.246) then
296     print *,'Problems with view:',KK
297     GOTO 45
298     endif
299     IF (KK.eq.3.and.status.ne.241) then
300     print *,'Problems with view:',KK
301     GOTO 45
302     endif
303     IF (KK.eq.4.and.status.ne.237) then
304     print *,'Problems with view:',KK
305     GOTO 45
306     endif
307    
308     C$$$ enddo
309     if (status.ne.-30460) then
310     c print *,'problems with DSP, event:',event
311     endif
312     C
313     41 format(1x,'status :',1x,Z4)
314     C
315     ic = ic + 1
316     lun(kk) = stringa(ic)
317     length = length + (stringa(ic) + 2)
318     C print *,'Lenght1 :',length
319     c write (*,41) stringa(ic)
320     C print *,'Lenght2 :',stringa(ic)
321    
322     C$$$ if (stringa(ic).ne.1064) then
323     C$$$ print *,'problems with DSP'
324     C$$$ endif
325     C
326     do i = 1,11
327     ic = ic + 1
328     auto(i) = stringa(ic)
329     enddo
330     C
331     do i = 1,11
332     do j = 1,96
333     ic = ic + 1
334     if (kk.eq.1) then
335     d1(i,j) = stringa(ic) -
336     & ped(kk,i,j)
337     if (i.eq.4.and.j.gt.64) then
338     c print *,'d1 :',d1(i,j),' ped :',ped(kk,i,j),
339     c & ' stringa :',stringa(ic),' j :',j
340     endif
341     endif
342     if (kk.eq.2) then
343     d2(i,j) = stringa(ic) -
344     & ped(kk,i,j)
345     endif
346     if (kk.eq.3) then
347     d3(i,j) = stringa(ic) -
348     & ped(kk,i,j)
349     endif
350     if (kk.eq.4) then
351     d4(i,j) = stringa(ic) -
352     & ped(kk,i,j)
353     C$ if (i.eq.2.and.j.gt.32.and.j.lt.49.and.event.eq.42)
354     C$ & then
355     C$ print *,'d4 :',d4(i,j),' ped :',ped(kk,i,j),
356     C$ & ' stringa :',stringa(ic),' j :',j
357     C$ endif
358     endif
359     enddo
360     enddo
361     C
362     DO I = 1,11
363     C
364     DO M = 1,6
365     INF = (M - 1) * 16 + 1
366     ISUP = M * 16
367     CHECK = 0.
368     C$ IF (M.EQ.3.AND.I.EQ.2.AND.KK.EQ.4.and.event.eq.42)
369     C$ & CHECK = 1.
370     IF (KK.EQ.1) CALL CALCBASE(KK,INF,ISUP,I,D1,BAS(KK,I,M))
371     c IF (KK.EQ.1.and.i.eq.4) PRINT *,'BASE 1:',BAS(KK,I,M),
372     c & ' m :',m
373     IF (KK.EQ.2) CALL CALCBASE(KK,INF,ISUP,I,D2,BAS(KK,I,M))
374     IF (KK.EQ.3) CALL CALCBASE(KK,INF,ISUP,I,D3,BAS(KK,I,M))
375     c IF (M.EQ.4.AND.I.EQ.4.AND.KK.EQ.4)
376     c & CHECK = 1.
377     IF (KK.EQ.4) CALL CALCBASE(KK,INF,ISUP,I,D4,BAS(KK,I,M))
378     ENDDO
379     C
380     C
381     do j = 1,96
382     IP = INT((J - 1) / 16) + 1
383     IF (KK.EQ.1) THEN
384     IF (BAS(1,I,IP).GT.0..AND.MIP1(I,J).GT.0.) THEN
385     DEDX1(I,J) = (D1(I,J) - BAS(1,I,IP)) / MIP1(I,J)
386     C if (i.eq.1) print *,'dedx :',dedx1(I,j),
387     C & ' base :',bas(1,i,ip),' mip :',mip1(i,j)
388     if (i.eq.4.and.j.gt.64) then
389     c print *,'dedx1 :',dedx1(i,j),' i :',i,' j :',j
390     endif
391     ELSE
392     DEDX1(I,J) = 0.
393     ENDIF
394     ENDIF
395     IF (KK.EQ.2) THEN
396     IF (BAS(2,I,IP).GT.0..AND.MIP2(I,J).GT.0.) THEN
397     DEDX2(I,J) = (D2(I,J) - BAS(2,I,IP)) / MIP2(I,J)
398     ELSE
399     DEDX2(I,J) = 0.
400     ENDIF
401     ENDIF
402     c IF (KK.EQ.3) THEN
403     c IF (BAS(3,I,IP).GT.0.) THEN
404     c DEDX3(I,J) = (D3(I,J) - BAS(3,I,IP)) / MIP3(I,J)
405     c ELSE
406     c DEDX3(I,J) = 0.
407     c ENDIF
408     c ENDIF
409     IF (KK.EQ.4) THEN
410     IF (BAS(4,I,IP).GT.0..AND.MIP4(I,J).GT.0.) THEN
411     DEDX4(I,J) = (D4(I,J) - BAS(4,I,IP)) / MIP4(I,J)
412     c IF (I.EQ.4.and.j.gt.48.and.j.lt.65)
413     c PRINT *,'DEDX4 :',DEDX4(I,J),' BASE :',
414     c & BAS(4,I,IP),' D4 :',D4(I,J),' J :',J
415     ELSE
416     DEDX4(I,J) = 0.
417     ENDIF
418     ENDIF
419     enddo
420     c
421     enddo
422     C
423     ic = ic + 1
424     C
425     45 continue
426     enddo
427     C
428     DO M = 1,4
429     DO K = 1,11
430     C
431     DO I = 1,96
432     IF (M.EQ.1.AND.MIP1(K,I).NE.0.) QPL(I) = DEDX1(K,I)
433     IF (M.EQ.2.AND.MIP2(K,I).NE.0.) QPL(I) = DEDX2(K,I)
434     IF (M.EQ.3.AND.MIP3(K,I).NE.0.) QPL(I) = DEDX3(K,I)
435     IF (M.EQ.4.AND.MIP4(K,I).NE.0.) QPL(I) = DEDX4(K,I)
436     ENE(I) = QPL(I)
437     QST(I) = QPL(I)
438     ENDDO
439     C
440     DO I = 1,96
441     IPRE = INT((I-1)/16) + 1
442     JC = MOD(I,32)
443     IF (QPL(I).GT.EMIN) THEN
444     C
445     INF = (IPRE - 1) * 16 + 1
446     ISUP = IPRE * 16
447     DO J = INF,ISUP
448     C$$$ ENE(J) = ENE(J) + QPL(I) * 0.00963
449     ENE(J) = ENE(J) + QPL(I) * 0.00478
450     ENDDO
451     C
452     IF (MOD(IPRE,2).EQ.0) IPRE2 = IPRE - 1
453     IF (MOD(IPRE,2).EQ.1) IPRE2 = IPRE + 1
454     INF = (IPRE2 - 1) * 16 + 1
455     ISUP = IPRE2 * 16
456     C$$$ DO J = INF,ISUP
457     C$$$ ENE(J) = ENE(J) + QPL(I) * 0.00485
458     C$$$ ENDDO
459     C
460     IF (JC.NE.1) ENE(I-1) = ENE(I-1) - QPL(I)*0.01581
461     IF (JC.NE.0) ENE(I+1) = ENE(I+1) - QPL(I)*0.01581
462     C
463     ENDIF
464     ENDDO
465     C
466     DO I = 1,96
467     IPRE = INT((I-1)/16) + 1
468     JC = MOD(I,32)
469     IF (ENE(I).GT.EMIN) THEN
470     C
471     INF = (IPRE - 1) * 16 + 1
472     ISUP = IPRE * 16
473     IDET = INT((I - 1) / 32) + 1
474     IF (M.EQ.1) DET = SPOS1(IDET,K)
475     IF (M.EQ.2) DET = SPOS2(IDET,K)
476     IF (M.EQ.3) DET = SPOS3(IDET,K)
477     IF (M.EQ.4) DET = SPOS4(IDET,K)
478     DO J = INF,ISUP
479     C$$$ QST(J) = QST(J) + ENE(I) * 0.00963
480     QST(J) = QST(J) + ENE(I) * (0.00478 +
481     & DET)
482     ENDDO
483     C
484     IF (MOD(IPRE,2).EQ.0) IPRE2 = IPRE - 1
485     IF (MOD(IPRE,2).EQ.1) IPRE2 = IPRE + 1
486     INF = (IPRE2 - 1) * 16 + 1
487     ISUP = IPRE2 * 16
488     DO J = INF,ISUP
489     QST(J) = QST(J) + ENE(I) * DET
490     ENDDO
491     C
492     IF (JC.NE.1) QST(I-1) = QST(I-1) - ENE(I)*0.01581
493     IF (JC.NE.0) QST(I+1) = QST(I+1) - ENE(I)*0.01581
494     C
495     ENDIF
496     ENDDO
497     C
498     DO I = 1,96
499     C IF (M.EQ.1.AND.MIP1(K,I).NE.0.) DEDX1(K,I) = QST(I)
500     c IF (M.EQ.2.AND.MIP2(K,I).NE.0.) DEDX2(K,I) = QST(I)
501     c IF (M.EQ.3.AND.MIP3(K,I).NE.0.) DEDX3(K,I) = QST(I)
502     c IF (M.EQ.4.AND.MIP4(K,I).NE.0.) DEDX4(K,I) = QST(I)
503     IF (M.EQ.4.AND.K.EQ.5) DEDX4(K,I) = 0.
504     ENDDO
505     C
506     ENDDO
507     C
508     ENDDO
509     C
510     CALL VZERO(DEXY,2*LENSEV)
511     CALL VZERO(QQ,4)
512     NSTRIP = 0.
513     QTOT = 0.
514     NX22 = 0.
515     QX22 = 0.
516     DO I = 1,11
517     DO J = 1,96
518     ICHECK3 = 1
519     IF (DEDX1(I,J).GT.EMIN)
520     & THEN
521     DEXY(2,2*I-1,97-J) = DEDX1(I,J)
522     NSTRIP = NSTRIP + 1.
523     QTOT = QTOT + DEDX1(I,J)
524     IF (I.NE.1) QQ(1) = QQ(1) + DEDX1(I,J)
525     c if (i.eq.1) print *,'qtot1 :',qtot,' dedx :',dedx1(i,j)
526     ENDIF
527     IF (DEDX2(I,J).GT.EMIN) THEN
528     DEXY(1,2*I-1,J) = DEDX2(I,J)
529     NSTRIP = NSTRIP + 1
530     QTOT = QTOT + DEDX2(I,J)
531     IF (I.EQ.11) THEN
532     NX22 = NX22 + 1.
533     QX22 = QX22 + DEDX2(I,J)
534     ENDIF
535     if (i.lt.11) QQ(2) = QQ(2) + DEDX2(I,J)
536     c print *,'qtot2 :',qtot
537     ENDIF
538     IF (DEDX3(I,J).GT.EMIN) THEN
539     C DEXY(2,2*I,J) = DEDX3(I,J)
540     C NSTRIP = NSTRIP + 1
541     C QTOT = QTOT + DEDX3(I,J)
542     C QQ(3) = QQ(3) + DEDX3(I,J)
543     c print *,'qtot3 :',qtot
544     ENDIF
545     IF (DEDX4(I,J).GT.EMIN) THEN
546     DEXY(1,2*I,J) = DEDX4(I,J)
547     NSTRIP = NSTRIP + 1
548     QTOT = QTOT + DEDX4(I,J)
549     IF (I.LT.11) QQ(4) = QQ(4) + DEDX4(I,J)
550     c if (i.eq.5) PRINT *,'DEDX4 :',DEDX4(I,J),' qtot :',
551     c & qtot,' J :',J,' I :',i
552     c print *,'qtot4 :',qtot
553     ENDIF
554     ENDDO
555     ENDDO
556     C
557     CALL CLUSTER
558     C
559     CALL DIRECTION(TG)
560     IMPX = CX - 241. / 2.
561     IMPY = CY - 241. / 2.
562     SHIFT = -0.5
563     CALL LASTRISCIA(CX,II)
564     IMPX = II
565     SHIFT = +0.5
566     CALL LASTRISCIA(CY,II)
567     IMPY = II
568     TANX = TG(1)
569     TANY = TG(2)
570     C
571     DO M = 1,2
572     DO I = 1,NPLA
573     C
574     NN = 0
575     IF (M.EQ.2) NN = 1
576     IF (MOD(I,2).EQ.NN) THEN
577     SHIFT = +0.5
578     ELSE
579     SHIFT = -0.5
580     ENDIF
581     C
582     DISTY = -PIANO * (I - 1.)
583     DISTX = -PIANO * (I - 1.)
584     C
585     IF (M.EQ.1) THEN
586     Y(I) = DISTX * TG(1) + CX
587     BAR(M,I) = Y(I)
588     C
589     ELSE
590     YY(I) = DISTY * TG(2) + CY
591     BAR(M,I) = YY(I)
592     C
593     ENDIF
594     CALL LASTRISCIA(BAR(M,I),IBAR(M,I))
595     ENDDO
596     ENDDO
597     C
598     NCORE = 0.
599     QCORE = 0.
600     PLANEMAX = 1.01*(LOG(ABS(RIG)/0.0081)-1.)
601     IPLANE = INT(ANINT(PLANEMAX)) + 5
602     IF (IPLANE.GT.NPLA) IPLANE=NPLA
603     DO J = 1,IPLANE
604     C
605     NNX = IBAR(1,J)
606     NNY = IBAR(2,J)
607     IF (NNX.LT.9) NNX = 9
608     IF (NNY.LT.9) NNY = 9
609     IF (NNX.GT.88) NNX = 88
610     IF (NNY.GT.88) NNY = 88
611     INFX = NNX - 8
612     INFY = NNY - 8
613     C
614     C 8 STRIPS ARE 2.88 cm , A MOLIERE RADIUS IS ABOUT 0.7 cm .
615     C
616     ISUPX = NNX + 8
617     ISUPY = NNY + 8
618     RNSS = 0.
619     QTOTT = 0.
620     DO I = INFX,ISUPX
621     IF (DEXY(1,J,I).GE.EMIN) THEN
622     RNSS = RNSS + 1
623     QTOTT = QTOTT + DEXY(1,J,I)
624     ENDIF
625     ENDDO
626     DO I = INFY,ISUPY
627     IF (DEXY(2,J,I).GE.EMIN) THEN
628     RNSS = RNSS + 1
629     QTOTT = QTOTT + DEXY(1,J,I)
630     ENDIF
631     ENDDO
632     NCORE = RNSS * FLOAT(J) + NCORE
633     QCORE = QTOTT * FLOAT(J) + QCORE
634     ENDDO
635     C
636     NINT = 0.
637     CALL NOINT(NIN) ! if NINT=1 not interacting particle
638     NINT = FLOAT(NIN)
639     C
640     C
641     QCYL = 0.
642     NCYL = 0.
643     C
644     C QCYL = DETECTED ENERGY AND NCYL = NUMBER OF HIT STRIPS IN A CYLINDER oF
645     C RADIUS 8.5 STRIPS WITH AXIS DEFINED BY THE DIRECTION OF THE INCOMING
646     C PARTICLE .
647     C
648     DO J = 1,NPLA
649     C
650     NNX = IBAR(1,J)
651     NNY = IBAR(2,J)
652     IF (NNX.LT.9) NNX = 9
653     IF (NNY.LT.9) NNY = 9
654     IF (NNX.GT.88) NNX = 88
655     IF (NNY.GT.88) NNY = 88
656     INFX = NNX - 8
657     INFY = NNY - 8
658     C
659     C 8 STRIPS ARE 2.88 cm , A MOLIERE RADIUS IS ABOUT 0.7 cm .
660     C
661     ISUPX = NNX + 8
662     ISUPY = NNY + 8
663     DO I = INFX,ISUPX
664     IF (DEXY(1,J,I).LT.EMIN) GO TO 710
665     NCYL = NCYL + 1
666     QCYL = QCYL + DEXY(1,J,I)
667     710 ENDDO
668     DO I=INFY,ISUPY
669     IF (DEXY(2,J,I).LT.EMIN) GO TO 810
670     NCYL = NCYL + 1
671     QCYL = QCYL + DEXY(1,J,I)
672     810 ENDDO
673     ENDDO
674     C
675     QTRACK = 0.
676     CALL LATERALE(QTRACK,RQT)
677     C
678     QMAX = 0.
679     DO M = 1,2
680     DO I = 1,NPLA
681     DO J = 1,NCHA
682     IF (DEXY(M,I,J).GT.QMAX) QMAX = DEXY(M,I,J)
683     ENDDO
684     ENDDO
685     ENDDO
686     C
687     if (ival.eq.1) iev = iev + 1
688     icount2 = icount2 + 1
689     if (icount2.ge.nev) goto 50
690    
691     c print *,'qtot :',qtot,' nstrip :',nstrip,' ncore :',ncore,
692     c & ' qcore :',qcore
693    
694     call hfnt(1)
695     C
696     goto 10
697    
698     50 continue
699    
700     print *,'Processed events :',icount2
701     print *,'Events with no calculated baseline :',iev
702    
703     *** /* Save histograms */
704    
705     CALL HROUT(0,icycle,' ')
706     CALL HREND('Event')
707     CLOSE (59)
708    
709     STOP
710     END
711    
712     C
713     C-------------------------------------------------------------------------
714     SUBROUTINE CALCBASE(K,IA,IB,IPL,DX,BASE)
715    
716     IMPLICIT NONE
717    
718     INTEGER IA, IB, IPL, K
719     REAL DX(11,96)
720     REAL S1, S2
721     REAL SIG
722     REAL C1
723     REAL V(16)
724     REAL RINF, RSUP
725     REAL BASE, BAS
726     REAL CHECK
727    
728     INTEGER J, I
729    
730     REAL PED(4,11,96),GOOD(4,11,96)
731    
732     COMMON / PEDE / PED, GOOD
733    
734     COMMON / CH / CHECK
735    
736     C$$$ print *,'ia :',ia,' ib :',ib,' di :',di,' ds :',ds
737    
738     S1 = 0.
739     S2 = 0.
740     C1 = 0.
741     DO J = IA,IB
742     IF (GOOD(K,IPL,J).EQ.0) THEN
743     S1 = S1 + DX(IPL,J)
744     S2 = S2 + DX(IPL,J)*DX(IPL,J)
745     C1 = C1 + 1.
746     ENDIF
747     ENDDO
748     C
749     BAS = 0.
750     SIG = 0.
751     IF (C1.GT.1) THEN
752     BAS = S1 / C1
753     SIG = SQRT(ABS((S2 - C1 * BAS*BAS) /(C1 - 1.)))
754     ENDIF
755     C
756     RINF = BAS - 3. * SIG
757     RSUP = BAS + 3. * SIG
758     C
759     S1 = 0.
760     S2 = 0.
761     C1 = 0.
762     DO J = IA,IB
763     IF (DX(IPL,J).GT.RINF.AND.DX(IPL,J).LT.RSUP
764     & .AND.GOOD(K,IPL,J).EQ.0) THEN
765     S1 = S1 + DX(IPL,J)
766     S2 = S2 + DX(IPL,J)*DX(IPL,J)
767     C1 = C1 + 1.
768     ENDIF
769     ENDDO
770     C
771     BAS = 0.
772     SIG = 0.
773     IF (C1.GT.1) THEN
774     BAS = S1 / C1
775     SIG = SQRT(ABS((S2 - C1 * BAS*BAS) /(C1 - 1.)))
776     ENDIF
777     C
778     IF (CHECK.EQ.1) PRINT *,'BAS :',BAS,' SIG :',SIG
779     IF (SIG.LT.30) THEN
780     BASE = BAS
781     RETURN
782     ENDIF
783     C
784     RINF = 1.E6
785     DO J = IA,IB
786     IF (DX(IPL,J).GT.0.AND.DX(IPL,J).LT.RINF.AND.
787     & GOOD(K,IPL,J).EQ.0) THEN
788     RINF = DX(IPL,J)
789     I = J
790     IF (CHECK.EQ.1) PRINT *,'RINF :',RINF,' I :',I
791     ENDIF
792     ENDDO
793     C
794     RINF = 1.E6
795     DO J = IA,IB
796     IF (DX(IPL,J).GT.0.AND.DX(IPL,J).LT.RINF.AND.
797     & GOOD(K,IPL,J).EQ.0.AND.J.NE.I) THEN
798     RINF = DX(IPL,J)
799     if (check.eq.1) print *,'rinf2 :',rinf
800     ENDIF
801     ENDDO
802     C
803     S1 = 0.
804     C1 = 0.
805     DO J = IA,IB
806     IF (DX(IPL,J).LT.(RINF+40.).AND.DX(IPL,J).GT.(RINF-40.)
807     & .AND.GOOD(K,IPL,J).EQ.0) THEN
808     S1 = S1 + DX(IPL,J)
809     C1 = C1 + 1.
810     if (check.eq.1) print *,'s1 :',s1,' c1 :',c1
811     ENDIF
812     ENDDO
813    
814     IF (C1.GT.1) BASE = S1 / C1
815     if (check.eq.1) print *,'BASE :',BASE
816     C PRINT *,'BASE :',BASE,' C1 :',C1
817    
818     RETURN
819     END
820    
821     C
822     C---------------------------------------------------------------------
823     SUBROUTINE LATERALE(RQT1,RQT2)
824     C---------------------------------------------------------------------
825     C RQT1 (IT WILL BE CALLED QTRACK IN THE N-TUPLE) IS THE SUM OF THE DETECTED
826     C ENERGY IN THE STRIP ALONG THE TRACK AND THE TWO CLOSEST STRIPS . FOR ALL THE
827     C LAYERS . RQT2 (IS NOT USED IN THE N-TUPLA) IS THE TOTAL ENERGY MINUS RQT1 .
828     C
829     INCLUDE 'INTEST.TXT'
830     REAL RQT1
831     INTEGER A,B
832     REAL BAR(2,NPLA)
833     REAL Q(0:NPLA)
834     INTEGER IBAR(2,NPLA)
835     COMMON/ANGOLO/BAR,IBAR
836    
837     RQT2=0.
838    
839     INPIA = 1
840     C
841     QQQ=0
842     MAX=0
843     Q(MAX)=0
844     C
845     DO I = INPIA,NPLA
846     A = IBAR(1,I)
847     B = IBAR(2,I)
848     IF (A.LE.2) A = 3
849     IF (B.LE.2) B = 3
850     IF (A.GE.(NCHA-1)) A = NCHA - 2
851     IF (B.GE.(NCHA-1)) B = NCHA - 2
852    
853     DO J = A-1,A+1
854     IF (DEXY(1,I,J).GE.EMIN) RQT1 = RQT1 + DEXY(1,I,J)
855     600 ENDDO
856     C
857     DO J = B-1,B+1
858     IF (DEXY(2,I,J).GE.EMIN) RQT1 = RQT1 + DEXY(2,I,J)
859     ENDDO
860     C
861     DO J=1,A-2
862     PXY = DEXY(1,I,J)
863     IF (PXY.GE.EMIN) RQT2 = RQT2 + PXY
864     650 ENDDO
865     C
866     DO J=A+2,NCHA
867     PXY = DEXY(1,I,J)
868     IF (PXY.GE.EMIN) RQT2 = RQT2 + PXY
869     700 ENDDO
870     C
871     DO J=1,B-2
872     PXY = DEXY(2,I,J)
873     IF (PXY.GE.EMIN) RQT2 = RQT2 + PXY
874     750 ENDDO
875     C
876     DO J=B+2,NCHA
877     PXY = DEXY(2,I,J)
878     IF (PXY.GE.EMIN) RQT2 = RQT2 + PXY
879     800 ENDDO
880     C
881     ENDDO
882     C
883     C
884     400 RETURN
885     END
886    
887     C------------------------------------------------
888     SUBROUTINE PEDESTAL(ERROR)
889     C------------------------------------------------
890    
891     IMPLICIT NONE
892     C
893     CHARACTER*6 file_name
894     C
895     integer icount,icount2
896     c
897     integer*4 evnum, systime, len
898     integer nev
899     integer RUNERROR, ERROR
900     C
901     C Normal variables definition
902     C
903     INTEGER FFD
904     C
905     integer event
906     C
907     INTEGER l, j, kk, m, mm, ll, ival, iev
908     INTEGER i, istat, icycle, ierr, icheck, icheck2, icheck3
909     INTEGER LVMIN
910    
911     integer nset,nword(10)
912     C
913     INTEGER*2 STRINGA(20000)
914     C
915     INTEGER ic, k, jj
916     INTEGER status
917     INTEGER inf, isup, idet
918    
919     INTEGER*2 length
920    
921     INTEGER lundata, iosop
922    
923     INTEGER*2 st, st2
924    
925     CHARACTER*2 nr
926     CHARACTER*1 alfa
927    
928     REAL auto(7)
929    
930     REAL RP,RI
931    
932     REAL SUM(6), SUM2(6), COUNT(6),
933     & BAS(4,11), DIS(6), DDIS(6)
934     REAL SSUM, SSUM2, CCOUNT, MME, SSIG
935    
936     REAL VAL
937    
938     REAL DI, DS
939     REAL ME(6), SIG(6)
940    
941     REAL PED(4,11,96), GOOD(4,11,96)
942     COMMON / PEDE / PED, GOOD
943    
944     C
945     C Begin !
946     C
947     RUNERROR=0 ! error variable
948    
949     PRINT *,'Pedestal file to read (only day and number included)? '
950     READ(*,905)file_name
951    
952     905 FORMAT(A6)
953    
954     lundata=44
955    
956     OPEN(UNIT=lundata,FILE='/wizard3/pamela/tb-0903/data/output_0309'
957     + //file_name//'.dat',STATUS='OLD',
958     + FORM='UNFORMATTED',ERR=50,IOSTAT=IOSOP)
959     C
960     ffd=FNum(lundata) ! take the file descriptor number
961     C
962    
963     C Read event: use C readevent routine
964     C
965    
966     10 continue
967     C
968     ival = 0
969     C
970     call readev(evnum,systime,len,stringa,RUNERROR,ffd)
971    
972     if (runerror.eq.-1.or.runerror.eq.1) goto 50
973    
974     ic = 0
975    
976     ic = ic + 1
977     length = ic
978     do kk = 1,4
979     ic = length + 1
980     k = ic
981     write (*,41) stringa(k)
982     st2 = IAND(stringa(k),'FF00'x)
983     status = ISHFT(st2,-8)
984     write (*,41) status
985     status = ISHFT(st2,-8)
986     IF (KK.eq.1.and.status.ne.170) then
987     print *,'Problems with view:',KK
988     ERROR = 1
989     GOTO 50
990     endif
991     IF (KK.eq.2.and.status.ne.182) then
992     print *,'Problems with view:',KK
993     ERROR = 1
994     GOTO 50
995     endif
996     IF (KK.eq.3.and.status.ne.177) then
997     print *,'Problems with view:',KK
998     ERROR = 1
999     GOTO 50
1000     endif
1001     IF (KK.eq.4.and.status.ne.173) then
1002     print *,'Problems with view:',KK
1003     ERROR = 1
1004     GOTO 50
1005     endif
1006     C
1007     41 format(1x,'status :',1x,Z4)
1008     C
1009     ic = ic + 1
1010     length = length + (stringa(ic) + 2)
1011     print *,'Lenght1 :',length
1012     print *,'Lenght2 :',stringa(ic)
1013    
1014     if (stringa(ic).ne.4629) then
1015     print *,'problems with view',KK
1016     ERROR = 1
1017     GOTO 50
1018     endif
1019     C
1020     do i = 1,11
1021     do j = 1,96
1022     ic = ic + 1
1023     ped(kk,i,j) = stringa(ic)
1024     good(kk,i,j) = stringa(ic+1)
1025     if (kk.eq.1.and.i.eq.4.and.j.eq.77)
1026     & print *,'ped :',stringa(ic),' g/b :',
1027     % stringa(ic+1)
1028     ic = ic + 1
1029     enddo
1030     enddo
1031     C
1032     enddo
1033     C
1034    
1035     50 continue
1036    
1037    
1038     CLOSE (44)
1039    
1040     RETURN
1041     END
1042    
1043    
1044    
1045    
1046    
1047    
1048    

  ViewVC Help
Powered by ViewVC 1.1.23