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

Contents of /calo/unpacking/nnplot_nt5.for

Parent Directory Parent Directory | Revision Log Revision Log


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

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