/[PAMELA software]/tracker/ground/source/reduction/reduction.f
ViewVC logotype

Contents of /tracker/ground/source/reduction/reduction.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1.1.1 - (show annotations) (download) (vendor branch)
Wed Mar 8 15:00:39 2006 UTC (18 years, 8 months ago) by pam-fi
Branch: MAIN, trk-ground
CVS Tags: R3v02, HEAD
Changes since 1.1: +0 -0 lines
First CVS release of tracker ground software (R3v02) 

1 *************************************************************************
2 *
3 * Program reduction.f
4 *
5 * - reads readraw.f output files: LEVEL0 ntuple, and ped, sig and bad histograms
6 * - decodes raw data (DATATRACKER) using DSP ped, sig and bad values
7 * - looks for clusters information using ped, sig and bad values from
8 * DSP histograms
9 * - fills LEVEL1 ntuple
10 *
11 *
12 *
13 * 11/9/2005 modified by david fedele to include general variables
14 *
15 *************************************************************************
16
17
18 program reduction
19
20 include '../common/commontracker.f'
21 include '../common/level0.f'
22 include '../common/level1.f'
23 include '../common/common_reduction.f'
24 include '../common/calib.f'
25
26
27 c------------------------------------------------------------------------
28 c
29 c local variables
30 c
31 c------------------------------------------------------------------------
32
33
34 parameter(ncalibmax=50)
35
36 c local variables
37 character*74 data_file !data file name
38 character*74 data_dir !data file name
39 character*74 data_file_calib !calibration list file
40 character*40 file_calib(ncalibmax)
41 character*74 data_file_level0
42 character*74 data_file_level1
43 integer minevent !first event to be analysed
44 integer nevent !number of events to be analysed
45 character*24 processing_date
46
47 integer which_calib_last
48
49 parameter (lun_data_file=70) !data file id number
50 parameter (lun_data_level0=71) !data file id number
51 parameter (lun_data_level1=72) !data file id number
52 parameter(lun_calib_list=66) !calibration list file id
53
54
55 integer ind
56 real value(nstrips_view) !per trovare i cluster
57 real clseedcut(nstrips_view)
58 real clinclcut(nstrips_view)
59 common/searchcluster/value,clseedcut,clinclcut,ind
60
61 integer nshowers
62 logical flag_shower
63 common/shower/nshowers,flag_shower
64
65 COMMON/QUEST/IQUEST(100) !permette di ottenere ntuple funzionanti nonostante
66 ! il messaggio dei 64K di RZOUT...!???
67
68
69 c*****************************************************
70 cccccc 11/9/2005 modified by david fedele
71 swcode1=202
72 c****************************************************
73 c------------------------------------------------------------------------
74 c
75 c HBOOK initialization
76 c
77 c------------------------------------------------------------------------
78
79 call HLIMIT(NWPAWC)
80
81
82 c------------------------------------------------------------------------
83 c
84 c reads input informations (through < go_reduction)
85 c
86 c------------------------------------------------------------------------
87
88 call fdate(processing_date)
89 write(*,101)
90 $ processing_date
91 101 format(/
92 $ ,'*** *** *** *** *** *** *** *** *** *** *** *** ***',/
93 $ ,'* *',/
94 $ ,'* REDUCTION *',/
95 $ ,'* *',/
96 $ ,'*** *** *** *** *** *** *** *** *** *** *** *** ***',/
97 $ ,a24,/
98 $ )
99
100 111 format(a)
101 print*,'Data directory:'
102 read(*,111)data_dir
103 print*,data_dir
104 print*,'File identifier: (DATE_NUM)'
105 read(*,*)data_file
106 print*,data_file
107 c$$$ print*,
108 c$$$ $ data_dir(1:LNBLNK(data_dir))
109 c$$$ $ //'DW_'//data_file(1:LNBLNK(data_file))
110 minevent=1
111 print*,'Maximum number of events to be analized:' !20000
112 read(*,*) nevent
113 print*,nevent
114 print*,'---------------------------------------------------'
115
116
117 c-------------------------------------------------------------
118 c-------------------------------------------------------------
119 c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
120 c
121 c reads calibration file
122 c
123 c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
124 c-------------------------------------------------------------
125 c-------------------------------------------------------------
126
127 print*,' '
128 print*,'OPENING CALIBRATION-LIST FILE:'
129 501 format(a,'DW_',a,'_calib.txt')
130 write(data_file_calib,501)
131 $ data_dir(1:LNBLNK(data_dir))
132 $ ,data_file(1:LNBLNK(data_file))
133 print*,data_file_calib
134 open(lun_calib_list,
135 $ FILE=data_file_calib(1:LNBLNK(data_file_calib))
136 $ ,STATUS='UNKNOWN'
137 $ ,IOSTAT=iostat
138 $ )
139
140 113 format(i5,' ',a25)
141
142 do i=1,ncalibmax
143
144 read(lun_calib_list,113,IOSTAT=iostat) n_cal_list
145 $ ,file_calib(i)(1:LNBLNK(file_calib(i)))
146 if(iostat.ne.0)then
147 ncal=i-1
148 goto 2000
149 endif
150 print*,n_cal_list,' - '
151 $ ,file_calib(i)(1:LNBLNK(file_calib(i)))
152
153 enddo
154
155 2000 close(lun_calib_list)
156 which_calib_last=0
157
158
159 c--------------------------------------------------------------
160 c--------------------------------------------------------------
161 c
162 c opening LEVEL 0 ntuple
163 c
164 c--------------------------------------------------------------
165 c--------------------------------------------------------------
166
167 504 format(a,'DW_',a,'_level0.rz')
168 write(data_file_level0,504)
169 $ data_dir(1:LNBLNK(data_dir))
170 $ ,data_file(1:LNBLNK(data_file))
171 print*,''
172 print*,'OPENING LEVEL0 FILE:'
173 print*,data_file_level0
174 IQUEST(10)=65000 !permette di ottenere ntuple funzionanti nonostante
175 ! il messaggio dei 64K di RZOUT...!???
176 call HROPEN(lun_data_level0,
177 $ 'LEVEL0',data_file_level0,'QP',4096,istat) !opens rz
178 if(istat.ne.0) goto 19
179 call HRIN(0,9999,0)
180 call access_level0
181 c call HPRNTU(ntp_level0)
182 call HNOENT(ntp_level0,iemax0)
183 print*,' events',iemax0
184
185 c------------------------------------------------------------------------
186 c
187 c LEVEL 1 ntuple booking
188 c
189 c------------------------------------------------------------------------
190
191 503 format(a,'DW_',a,'_level1.rz')
192 write(data_file_level1,503)
193 $ data_dir(1:LNBLNK(data_dir))
194 $ ,data_file(1:LNBLNK(data_file))
195 print*,''
196 print*,'------------------------------------'
197 print*,' Creating LEVEL1 rz file'
198 print*,' ',data_file_level1
199 print*,'------------------------------------'
200
201 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
202 C largest RZ file: IQUEST(10) records x LREC words x 4 byte
203 C with the following settings: 65000 x 4096 x 4 = 1G
204 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
205 IQUEST(10)=65000
206 c permette di ottenere ntuple funzionanti nonostante
207 c il messaggio dei 64K di RZOUT
208 call HROPEN(lun_data_level1,
209 $ 'LEVEL1',data_file_level1,'QNP',4096,istat) !opens rz
210 if(istat.ne.0) goto 19
211 call book_level1
212
213 * -------------------------------------------------------
214 * STRIP MASK
215 * -------------------------------------------------------
216 call stripmask(data_file)
217
218
219 nshowers=0
220 nocalib=0
221 c counts number of events with more than nclstrmax clusters
222 maxevent=minevent+nevent-1
223 nevent = 0
224 nevent_good = 0
225
226 do iev = minevent,MIN(iemax0,maxevent) !loop on events
227
228 nevent = nevent + 1
229 if(mod(nevent,100).eq.0)print*,nevent
230
231 call init_level1
232
233 call HCDIR('//LEVEL0',' ')
234 call HGNT(ntp_level0,iev,ierr) !read event at LEVEL0
235 if(ierr.ne.0)then
236 print*,'HGNT ---> error reading event ',iev
237 good1 = .false.
238 nev1 = nev1 + 1
239 goto 200 !jump to fill n-tuple
240 endif
241 c if(mod(iev,100).eq.0)
242 c $ print*,' event number: ',iev
243 C---------------------------------------------------
244
245 c*****************************************************
246 cccccc 11/9/2005 modified by david fedele
247 ccC variables in blocks EVENT and CPU are anyway filled
248 C variables in blocks GENERAL and CPU are anyway filled
249 c********************************************************
250 C in order to mantain sincronization among
251 C events at different levels
252 C---------------------------------------------------
253 c*****************************************************
254 cccccc 11/9/2005 modified by david fedele
255 c good1=good0
256 c nev1=nev0
257 c pkt_type1=pkt_type
258 c pkt_num1=pkt_num
259 c obt1=obt
260 c which_calib1=which_calib
261 good1=good0
262 nev1=nev0
263 which_calib1=which_calib
264 pkt_type1=pkt_type
265 pkt_num1=pkt_num
266 obt1=obt
267 cpu_crc1=cpu_crc
268 do iv=1,12
269 crc1(iv)=crc(iv)
270 enddo
271 c*******************************************************
272
273 C---------------------------------------------------
274 C retrieve info about calib to be used
275 C---------------------------------------------------
276
277 if(which_calib.ne.which_calib_last.and.
278 $ which_calib.ne.0)then
279
280 data_file_calib=
281 $ data_dir(1:LNBLNK(data_dir))//
282 $ file_calib(which_calib)
283 $ (1:LNBLNK(file_calib(which_calib)))
284 c print*,data_file_calib
285 print*,''
286 print*,
287 $ '@ event ',nev0
288 $ ,' (CPU pkt N.',pkt_num,')'
289 print*,'--> ',data_file_calib
290 IQUEST(10)=65000
291 call HROPEN(lun_data_file,
292 $ 'CALIB',data_file_calib,'QP',4096,istat) !opens
293 if(istat.ne.0) goto 19
294 call HRIN(0,9999,0)
295 call fillpedsig
296 do iview=1,nviews
297 call HDELET(id_hi_bad+iview)
298 call HDELET(id_hi_ped+iview)
299 call HDELET(id_hi_sig+iview)
300 enddo
301 call HREND('CALIB')
302 close(lun_data_file)
303 which_calib_last=which_calib
304
305
306 elseif(which_calib.eq.0)then
307
308 nocalib=nocalib+1
309 good1=.false.
310
311 endif
312
313 if(.not.good1) goto 200 !jump to fill n-tuple
314
315 c--------------------------------------------------
316 c read the variable DATATRACKER from LEVEL0
317 c and fill the variable ADC (inverting view 11)
318 c--------------------------------------------------
319 call filladc(iflag)
320 if(iflag.ne.0)then
321 good1=.false.
322 print*,'event ',iev,' >>>>> decode ERROR'
323 goto 200
324 endif
325 ccc /////////////////////////////////////////////////// da qui
326 c--------------------------------------------------
327 c computes common noise for each VA1
328 c (excluding strips affected by signal,
329 c tagged with the flag CLSTR)
330 c--------------------------------------------------
331 do iv=1,nviews
332 do ik=1,nva1_view
333 cn(iv,ik)=0 !initializes cn variable
334 if(mask_vk(iv,ik).eq.1)call cncomp(iv,ik)
335 enddo
336 enddo
337 c---------------------------------------------
338 c loops on views, VA1 and strips,
339 c and computes strips signals using
340 c badstrip, pedestals, and
341 c sigma informations from histograms
342 c---------------------------------------------
343 flag_shower = .false.
344 ind=1 !clsignal array index
345 do iv=1,nviews !loop on views
346 do is=1,nstrips_view !loop on strips (1)
347 if(mod(iv,2).eq.1) then
348 C=== > Y view
349 value(is)= -(DBLE(adc(iv,nvk(is),nst(is)))
350 $ -cn(iv,nvk(is))-pedestal(iv,nvk(is),nst(is)))
351 $ *mask(iv,nvk(is),nst(is))
352 clseedcut(is)=clcuty*sigma(iv,nvk(is),nst(is))
353 $ *mask(iv,nvk(is),nst(is))
354 clinclcut(is)=incuty*sigma(iv,nvk(is),nst(is))
355 $ *mask(iv,nvk(is),nst(is))
356 else
357 C=== > X view
358 value(is)= (DBLE(adc(iv,nvk(is),nst(is)))
359 $ -cn(iv,nvk(is))-pedestal(iv,nvk(is),nst(is)))
360 $ *mask(iv,nvk(is),nst(is))
361 clseedcut(is)=clcutx*sigma(iv,nvk(is),nst(is))
362 $ *mask(iv,nvk(is),nst(is))
363 clinclcut(is)=incutx*sigma(iv,nvk(is),nst(is))
364 $ *mask(iv,nvk(is),nst(is))
365 endif
366 enddo !end loop on strips (1)
367 call search_cluster(iv)
368 if(flag_shower.eqv..true.)then
369 call init_level1
370 good1=.false.
371 nev1 = nev0
372 goto 200 !jump to next event
373 endif
374 enddo ! end loop on views
375 do iv=1,nviews
376 do ik=1,nva1_view
377 cnev(iv,ik)=cn(iv,ik) !assigns computed CN to ntuple variables
378 enddo
379 enddo
380 nevent_good = nevent_good + 1
381
382 200 continue
383 C---------------------------------------------
384 C come here if GOOD1=.false.
385 C or the event has too many clusters
386 C---------------------------------------------
387
388 call HCDIR('//LEVEL1',' ')
389 call HFNT(ntp_level1) !fills LEVEL1 ntuple
390
391 ccc /////////////////////////////////////////////////// a qui
392
393 enddo !end loop on events
394
395 c$$$ call HPRNT(ntp_level1)
396
397 write(*,108)
398 $ nevent
399 $ ,nevent_good
400 $ ,nshowers
401 $ ,nclstrmax
402 $ ,nocalib
403 108 format(/
404 $ ,'|------------------------------------------------|',/
405 $ ,'| Number of processed events |',/
406 $ ,'| Total ',i8,' |'/
407 $ ,'| Good ',i8,' |'/
408 $ ,'| (',i5,' events with > ',i5,' clusters) |'/
409 $ ,'| (',i5,' events without calibration) |'/
410 $ ,'|------------------------------------------------|',/
411 $ )
412
413 goto 9000
414
415 c------------------------------------------------------------------------
416 c
417 c no error exit
418 c
419 c------------------------------------------------------------------------
420
421 print*,' '
422 print*,'REDUCTION SUCCESSFULLY COMPLETED'
423 print*,' '
424 print*,' '
425
426 goto 9000 !happy ending
427
428
429 c------------------------------------------------------------------------
430 c
431 c data file opening error
432 c
433 c------------------------------------------------------------------------
434
435 19 continue
436
437 print*,' '
438 print*,'reduction: ERROR OPENING DATA FILE: ',data_file
439 print*,' '
440 print*,' '
441
442 goto 9000 !the end
443
444
445 c------------------------------------------------------------------------
446 c
447 c data file event error
448 c
449 c------------------------------------------------------------------------
450
451 20 continue
452
453 print*,' '
454 print*,'reduction: ERROR IN DATA FILE: ',data_file
455 print*,'NO EVENTS FOUND'
456 print*,' '
457 print*,' '
458
459 goto 9000 !the end
460
461
462 c------------------------------------------------------------------------
463 c
464 c data file eleboration level error
465 c
466 c------------------------------------------------------------------------
467
468 21 continue
469
470 print*,' '
471 print*,'reduction: ERROR IN DATA FILE: ',data_file
472 print*,'LEVEL0= ',level0
473 print*,' '
474 print*,' '
475
476 goto 9000 !the end
477
478
479 c------------------------------------------------------------------------
480 c
481 c level0 ntuple event reading error
482 c
483 c------------------------------------------------------------------------
484
485 22 continue
486
487 print*,' '
488 print*,'reduction: ERROR WHILE READING LEVEL0 NTUPLE, AT EVENT: ',
489 $ iev
490 print*,ierr
491 print*,' '
492 print*,' '
493
494 goto 9000 !the end
495
496
497 c------------------------------------------------------------------------
498 c
499 c closes files and exits
500 c
501 c------------------------------------------------------------------------
502
503 9000 continue
504
505
506 call HCDIR('//LEVEL1',' ')
507 call RZPURG(-1)
508 call HROUT(ntp_level1,ICYCLE,'T')
509 call HREND('level1')
510 close(lun_data_level1)
511
512 call HCDIR('//LEVEL0',' ')
513 call HREND('level0')
514 close(lun_data_level0)
515
516 c$$$ close(lun_data_runinfo)
517
518
519 stop
520 end
521 ***...***...***...***...***...***...***...***...***...***...***...***...***...***...***...***
522 *
523 *
524 *
525 *
526 *
527 *
528 *
529 *
530 *
531 ***...***...***...***...***...***...***...***...***...***...***...***...***...***...***...***
532 subroutine book_level1
533
534 include '../common/commontracker.f'
535 include '../common/level1.f'
536
537 character*35 block1,block2,block3,
538 $ block4,block5,block6,block7
539 $ ,block8
540
541 * LEVEL1 ntuple:
542 call HBNT(ntp_level1,'LEVEL1',' ')
543
544 c*****************************************************
545 cccccc 11/9/2005 modified by david fedele
546 c call HBNAME(ntp_level1,'EVENT',good1,'GOOD1:L,NEV1:I')
547 c call HBNAME(ntp_level1,'GENERAL',good1,'GOOD1:L,NEV1:I
548 c $ ,WHIC_CALIB1:I,SWCODE1:I')
549 cccccc 06/10/2005 modified by elena vannuccini
550 call HBNAME(ntp_level1,'GENERAL',good1,'GOOD1:L,NEV1:I
551 $ ,WHIC_CALIB1:I,SWCODE1:I,CRC1(12):L')
552 c*********************************************************
553 call HBNAME(ntp_level1,'CPU',pkt_type1
554 $ ,'PKT_TYPE1:I
555 $ ,PKT_NUM1:I
556 $ ,OBT1:I'//
557 c********************************************************
558 cccccc 11/9/2005 modified by david fedele
559 c $ ,WHICH_CALIB1:I')
560 $ ',CPU_CRC1:L')
561 c********************************************************
562 411 format('NCLSTR1:I::[0,',I4,']') !cluster number
563 412 format(',VIEW(NCLSTR1):I::[0,',I2,']') !view the cluster belongs to
564 413 format(',LADDER(NCLSTR1):I::[0,',I1,']') !ladder the cluster belongs to
565 414 format(',MAXS(NCLSTR1):I::[0,',I5,']') !cluster stip carrying the largest
566 ! signal value
567 415 format(',MULT(NCLSTR1):I::[0,',I2,']') !cluster multiplicity
568 416 format(',INDSTART(NCLSTR1):I::[0,',I5,']') !cluster starting point index in
569 ! clsignal array
570 417 format(',INDMAX(NCLSTR1):I::[0,',I5,']') !cluster maximum point index in
571 $ ! clsignal array
572 418 format(',TOTCLLENGTH:I::[0,',I5,']') !sum of all clusters length
573 write(block1,411) nclstrmax !maximum number of clusters per event
574 write(block2,412) nviews !number of views
575 write(block3,413) nladders_view !number of ladders per view
576 write(block4,414) nstrips_view !number of strips per view
577 write(block5,415) nclstrp !maximum number of strips in a cluster
578 write(block6,416) maxlength !maximum number of strip belonging to clusters
579 write(block7,417) maxlength ! for the whole event
580 write(block8,418) maxlength
581 call HBNAME(ntp_level1,'CLUSTER',nclstr1,block1//block2//block3
582 $ //block4//block5//',DEDX(NCLSTR1):r'//block6//block7
583 c*****************************************************
584 cccccc 11/9/2005 modified by david fedele
585 c $ //block8//',CLSIGNAL(TOTCLLENGTH):R')
586 c $ //block8//',CLSIGNAL(TOTCLLENGTH):R,CRC1(12):L')
587 cccccc 06/10/2005 modified by elena vannuccini
588 $ //block8//',CLSIGNAL(TOTCLLENGTH):R')
589 c*****************************************************
590 419 format('CNEV(',I4,',',I4,')') !common noise of the event for a certain view
591 write(block8,419) nviews,nva1_view ! and VA1
592 call HBNAME(ntp_level1,'CNOISE',cnev,block8)
593
594
595 return
596 end
597
598 c.............................................................
599
600 subroutine init_level1
601
602 include '../common/commontracker.f'
603 include '../common/level1.f'
604 include '../common/level0.f'
605
606 nclstr1=0
607 totCLlength=0
608 do ic=1,nclstrmax
609 view(ic)=0
610 ladder(ic)=0
611 indstart(ic)=0
612 indmax(ic)=0
613 maxs(ic)=0
614 mult(ic)=0
615 dedx(ic)=0
616 enddo
617 do id=1,maxlength !???
618 clsignal(id)=0.
619 enddo
620 do iv=1,nviews
621 c*****************************************************
622 cccccc 11/9/2005 modified by david fedele
623 crc1(iv)=.false.
624 c**************************************************
625 do ik=1,nva1_view
626 cnev(iv,ik)=0
627 enddo
628 enddo
629
630 return
631 end
632 *---***---***---***---***---***---***---***---***
633 *
634 *
635 *
636 *
637 *
638 *---***---***---***---***---***---***---***---***
639
640 subroutine search_cluster(iv)
641
642 include '../common/commontracker.f'
643 include '../common/level1.f'
644 include '../common/calib.f'
645
646 integer ind
647 real value(nstrips_view) !per trovare i cluster
648 real clseedcut(nstrips_view)
649 real clinclcut(nstrips_view)
650 common/searchcluster/value,clseedcut,clinclcut,ind
651
652 integer nshowers
653 logical flag_shower
654 common/shower/nshowers,flag_shower
655
656 c local variables
657 integer rmax,lmax !estremi del cluster
658 integer rstop,lstop !per decidere quali strip includere nel cluster
659 ! oltre il seed
660 integer first,last,diff !per includere le strip giuste... !???
661
662 integer multtemp !temporary multiplicity variable
663
664 integer CLlength !lunghezza in strip del cluster
665
666 c external nst
667
668 c------------------------------------------------------------------------
669 c looks for clusters on each view
670 C : CERCO STRIP SOPRA CLSEEDCUT, POI SCORRO A DX FINCHE'
671 c NON TROVO
672 C STRIP PIU' BASSA (in segnale/rumore)
673 C => L'ULTIMA DELLA SERIE CRESCENTE
674 C (LA PIU' ALTA) E' IL
675 C CLUSTER SEED. POI SCORRO A SX E DX INCLUDENDO TUTTE
676 C LE STRIP (FINO A 17 AL
677 C MAX) CHE SUPERANO CLINCLCUT.
678 C QUANDO CERCO IL CLUSTER SEED SUCCESSIVO SALTO LA STRIP
679 C ADIACENTE A DESTRA
680 C DELL'ULTIMO CLUSTER SEED (CHE SARA' NECESSARIAMENTE
681 C PIU' BASSA) E PRENDO
682 C COME SEED UNA STRIP SOLO SE IL SUO SEGNALE E'
683 C MAGGIORE DI QUELLO DELLA STRIP
684 C PRECEDENTE (PRATICAMENTE PER EVITARE CHE L'ULTIMA
685 C STRIP DI UN GRUPPO DI STRIP
686 C TUTTE SOPRA IL CLSEEDCUT VENGA AUTOMATICAMENTE PRESA
687 C COME SEED... DEVE ESSERE
688 C PRESA SOLO SE IL CLUSTER E' DOUBLE PEAKED...)
689 c------------------------------------------------------------------------
690 c 6 ottobre 2003
691 c Elena: CLSEEDCUT = 7 (old value 10)
692 c Elena: CLINCLCUT = 4 (old value 5)
693
694 iseed=-999 !cluster seed index initialization
695
696 do jl=1,nladders_view !1..3 !loops on ladders
697 first=1+nstrips_ladder*(jl-1) !1,1025,2049
698 last=nstrips_ladder*jl !1024,2048,3072
699 c X views have 1018 strips instead of 1024
700 if(mod(iv,2).eq.0) then
701 first=first+3
702 last=last-3
703 endif
704 do is=first,last !loop on strips in each ladder
705 if(is.le.iseed+1) goto 220
706 c-----------------------------------------
707 c after a cluster seed as been found,
708 c look for next one skipping one strip on the right
709 c (i.e. look for double peak cluster)
710 c-----------------------------------------
711 if(is.ne.first) then
712 if(value(is).le.value(is-1)) goto 220
713 endif
714 c-----------------------------------------
715 c skips cluster seed
716 c finding if strips values are descreasing (a strip
717 c can be a cluster seed only if previous strip value
718 c is lower)
719 c-----------------------------------------
720 if(value(is).gt.clseedcut(is)) then
721 c-----------------------------------------
722 c possible SEED...
723 c-----------------------------------------
724 itemp=is
725 if(itemp.eq.last) goto 230 !estremo...
726 c$$$ do while(value(itemp)
727 c$$$ $ /sigma(iv,nvk(itemp),nst(itemp))
728 c$$$ $ .le.value(itemp+1)
729 c$$$ $ /sigma(iv,nvk(itemp+1),nst(itemp+1))) !BIAS: aggiustare il caso uguale!???
730 do while(value(itemp)
731 $ .le.value(itemp+1)) !BIAS: aggiustare il caso uguale!???
732 itemp=itemp+1
733 if(itemp.eq.last) goto 230 !stops if reaches last strip
734 enddo ! of the ladder
735 230 continue
736 c-----------------------------------------
737 c fownd SEED!!!
738 c-----------------------------------------
739 iseed=itemp
740 c----------------------------------------------------------
741 c after finding a cluster seed, checks also adjacent strips,
742 C and marks the ones exceeding clinclcut
743 c----------------------------------------------------------
744 ir=iseed !indici destro
745 il=iseed ! e sinistro
746
747 rmax=ir !estremo destro del cluster
748 lmax=il ! e sinistro
749
750 rstop=0 !initialize flags used to exit from
751 lstop=0 ! inclusion loop
752
753 do while(lstop.eq.0.or.rstop.eq.0) !shifts left and right from
754 ir=ir+1 !position index for strips on right side of
755 ! cluster seed
756 il=il-1 !and for left side
757 c------------------------------------------------------------------------
758 c checks for last or first strip of the ladder
759 c------------------------------------------------------------------------
760 if(ir.gt.last) then !when index goes beyond last strip
761 rstop=1 ! of the ladder, change rstop flag in order
762 ! to "help" exiting from loop
763 endif
764
765 if(il.lt.first) then !idem when index goes beyond
766 lstop=1 ! first strip of the ladder
767 endif
768
769 c------------------------------------------------------------------------
770 c check for clusters including more than nclstrp strips
771 c------------------------------------------------------------------------
772 if((rmax-lmax+1).ge.nclstrp) then
773 goto 210 !exits inclusion loop:
774 ! lmax and rmax maintain last value
775 ! NB .ge.!???
776 endif
777 c------------------------------------------------------------------------
778 c marks strips exceeding inclusion cut
779 c------------------------------------------------------------------------
780 if(rstop.eq.0) then !if last strip of the ladder or last
781 ! over-cut strip has not been reached
782 if(value(ir).gt.clinclcut(ir)) then !puts in rmax the
783 rmax=ir ! last right over-cut strip
784 else
785 rstop=1 !otherwise cluster ends on right and rstop
786 endif ! flag=1 signals it
787 endif
788 if(lstop.eq.0) then
789 if(value(il).gt.clinclcut(il)) then
790 lmax=il
791 else
792 lstop=1
793 endif
794 endif
795
796 enddo !ends strip inclusion loop
797 210 continue !jumps here if more than nclstrp have been included
798
799 multtemp=rmax-lmax+1 !stores multiplicity in temp
800 ! variable. NB rmax and lmax can change later in
801 ! order to include enough strips to calculate eta3
802 ! and eta4. so mult is not always equal to cllength
803 c------------------------------------------------------------------------
804 c NB per essere sicuro di poter calcolare eta3 e eta4 devo includere
805 c sempre e comunque le 2 strip adiacenti al cluster seed e quella
806 c adiacente ulteriore dalla parte della piu' alta fra queste due
807 c (vedi oltre...)!???
808 c------------------------------------------------------------------------
809
810 c nel caso di estremi del ladder...!???
811
812 c ho meno di 4 strip nel cluster --> se sono sui bordi o quasi del ladder
813 c costruisco il cluster ad hoc e poi esco, se non sono sui bordi o quasi
814 c vado oltre (aggiungero' quindi strip a sx e dx in modo da poter calcolare
815 c eta3e4)
816 if((rmax-lmax+1).lt.4) then
817
818 if(iseed.eq.first) then !estremi...
819 rmax=iseed+2 !NB in questo modo puo' anche capitare di
820 lmax=iseed ! includere strip sotto taglio di inclusione
821 goto 250 ! che non serviranno per eta3e4!???
822 endif
823
824 if(iseed.eq.last) then !estremi...
825 rmax=iseed
826 lmax=iseed-2 !NB 2 e non 3, perche' altrimenti sarei in
827 goto 250 ! ((rmax-lmax+1).lt.4).eq.false. !???
828 endif !NMB questo e' l'unico caso di cllength=3!???
829
830 if(iseed.eq.first+1) then !quasi estremi...
831 rmax=iseed+2
832 lmax=iseed-1
833 goto 250
834 endif
835 if(iseed.eq.last-1) then
836 rmax=iseed+1
837 lmax=iseed-2
838 goto 250
839 endif
840 c se ho 4 o piu' strip --> se sono sui bordi esco, se sono sui quasi bordi
841 c includo la strip del bordo
842 else
843
844 if(iseed.eq.first) goto 250 !estremi... non includo altro
845 if(iseed.eq.last) goto 250
846 if(iseed.eq.first+1) then !quasi estremi... mi assicuro di
847 lmax=first ! avere le strip adiacenti al seed
848 if((rmax-lmax+1).gt.nclstrp) rmax=rmax-1 !NB effetto
849 goto 250 ! coperta: se la lunghezza del cluster era gia'
850 endif ! al limite (nclstrp), per poter aggiungere questa
851 ! strip a sinistra devo toglierne una a destra...!???
852 if(iseed.eq.last-1) then
853 rmax=last
854 if((rmax-lmax+1).gt.nclstrp) lmax=lmax+1
855 goto 250
856 endif
857 endif
858 c------------------------------------------------------------------------
859 c be sure to include in the cluster the cluster seed with its 2 adjacent
860 c strips, and the one adjacent to the greatest between this two strip, as the
861 c fourth one. if the strips have the same value (!) the fourth one is chosen
862 c as the one having the greatest value between the second neighbors
863 c------------------------------------------------------------------------
864 if(value(iseed+1).eq.value(iseed-1)) then
865 if(value(iseed+2).ge.value(iseed-2)) then !??? qui cmq c'e'
866 diff=(iseed+2)-rmax
867 if(diff.gt.0) then
868 rmax=rmax+diff
869 if((rmax-lmax+1).gt.nclstrp) then
870 lmax=rmax-nclstrp+1
871 endif
872 endif
873 diff=(iseed-1)-lmax
874 if(diff.lt.0) then
875 lmax=lmax+diff
876 if((rmax-lmax+1).gt.nclstrp) then
877 rmax=lmax+nclstrp-1
878 endif
879 endif
880 else
881 diff=(iseed-2)-lmax
882 if(diff.lt.0) then
883 lmax=lmax+diff
884 if((rmax-lmax+1).gt.nclstrp) then
885 rmax=lmax+nclstrp-1
886 endif
887 endif
888 diff=(iseed+1)-rmax
889 if(diff.gt.0) then
890 rmax=rmax+diff
891 if((rmax-lmax+1).gt.nclstrp) then
892 lmax=rmax-nclstrp+1
893 endif
894 endif
895
896 endif
897 elseif(value(iseed+1).gt.value(iseed-1)) then
898 c !??? sposto il limite del cluster a destra per includere sempre le strip
899 c necessarie al calcolo di eta-i
900 c se il cluster diventa troppo lungo lo accorcio a sinistra per avere non piu'
901 c di nclstrp (in questo caso sono sicuro di aver gia' incluso le strip
902 c necessarie al calcolo di eta-i a sinistra, quindi se voglio posso uscire)
903 diff=(iseed+2)-rmax
904 if(diff.gt.0) then
905 rmax=rmax+diff
906 if((rmax-lmax+1).gt.nclstrp) then
907 lmax=rmax-nclstrp+1
908 c goto 250
909 endif
910 endif
911 diff=(iseed-1)-lmax
912 if(diff.lt.0) then
913 lmax=lmax+diff
914 if((rmax-lmax+1).gt.nclstrp) then
915 rmax=lmax+nclstrp-1
916 c goto 250 !inutile!???
917 endif
918 endif
919 else
920 diff=(iseed-2)-lmax
921 if(diff.lt.0) then
922 lmax=lmax+diff
923 if((rmax-lmax+1).gt.nclstrp) then
924 rmax=lmax+nclstrp-1
925 c goto 250
926 endif
927 endif
928 diff=(iseed+1)-rmax
929 if(diff.gt.0) then
930 rmax=rmax+diff
931 if((rmax-lmax+1).gt.nclstrp) then
932 lmax=rmax-nclstrp+1
933 c goto 250 !inutile!???
934 endif
935 endif
936 endif
937 250 continue
938
939 c--------------------------------------------------------
940 c fills ntuple variables
941 c--------------------------------------------------------
942 nclstr1=nclstr1+1 !cluster number
943
944 if(nclstr1.gt.nclstrmax) then !too many clusters for the event:
945 nshowers=nshowers+1 ! skips variables filling and go to next
946 good1=.false. ! event
947 nclstr1=0
948 totCLlength=0
949 flag_shower = .true.
950 print*,'Event ',nev1,
951 $ ': more than ',nclstrmax,' clusters'
952 goto 2000
953 endif
954 view(nclstr1)=iv !vista del cluster
955 ladder(nclstr1)=nld(iseed,iv) !ladder a cui appartiene il cluster seed
956 maxs(nclstr1)=iseed !strip del cluster seed
957 mult(nclstr1)=multtemp !molteplicita'
958
959 indstart(nclstr1)=ind !posizione dell'inizio del cluster nell'
960 ! array clsignal
961 indmax(nclstr1)=indstart(nclstr1)+(iseed-lmax) !posizione del
962 ! cluster seed nell'array clsignal
963
964 CLlength=rmax-lmax+1 !numero di strip del cluster
965 totCLlength=totCLlength+CLlength
966 dedx(nclstr1)=0
967 do j=lmax,rmax !stores sequentially cluter strip values in
968 clsignal(ind)=value(j) ! clsignal array
969 ind=ind+1
970 c if(value(j).gt.0)
971 if(value(j).gt.clinclcut(j))
972 $ dedx(nclstr1)=dedx(nclstr1)+value(j) !cluster charge
973 enddo
974 c--------------------------------------------------------
975 c
976 c--------------------------------------------------------
977 endif !end possible seed conditio
978 220 continue !jumps here to skip strips left of last seed
979
980 enddo ! end loop on strips
981 enddo !end loop on ladders
982 2000 continue
983 return
984 end
985
986
987 *---***---***---***---***---***---***---***---***
988 *
989 *
990 *
991 *
992 *
993 *---***---***---***---***---***---***---***---***
994
995
996 subroutine stripmask(data_file)
997
998 * this routine set va1 and single-strip masks,
999 * on the basis of date and id downlink
1000 *
1001 * mask(nviews,nva1_view,nstrips_va1) !strip mask
1002 * mask_vk(nviews,nva1_view) !VA1 mask
1003 *
1004 include '../common/commontracker.f'
1005 include '../common/level1.f'
1006 include '../common/calib.f'
1007
1008 character*10 data_file
1009
1010 character*3 aid
1011 character*6 adate
1012 integer id
1013 integer date
1014
1015 * ----------------------
1016 * retrieve date and id
1017 aid=data_file(8:10)
1018 adate=data_file(1:6)
1019 READ (aid, '(I3)'), id
1020 READ (adate, '(I6)'), date
1021 * ----------------------
1022
1023 * init mask
1024 do iv=1,nviews
1025 do ivk=1,nva1_view
1026 mask_vk(iv,ivk)=1
1027 do is=1,nstrips_va1
1028 mask(iv,ivk,is)=1
1029 enddo
1030 enddo
1031 enddo
1032
1033 * ---------------------
1034 * VIEW 2 - VK 23-24
1035 * couple of vk damaged during integration
1036 if(date.ge.50208)then
1037 print*,'MASK: view 2 - vk 23/24'
1038 mask_vk(2,23)=0
1039 mask_vk(2,24)=0
1040 do is=1,nstrips_va1
1041 mask(2,23,is)=0
1042 mask(2,24,is)=0
1043 enddo
1044 endif
1045
1046 * ---------------------
1047 * VIEW 7 - VK 11-12
1048 if(date.ge.50209)then
1049 if(.not.(date.eq.50209.and.id.le.6)) then
1050 print*,'MASK: view 7 - vk 11/12'
1051 mask_vk(7,11)=0
1052 mask_vk(7,12)=0
1053 do is=1,nstrips_va1
1054 mask(7,11,is)=0
1055 mask(7,12,is)=0
1056 enddo
1057 endif
1058 endif
1059
1060 * ---------------------
1061 * VIEW 7 - VK 21-22
1062 if(date.ge.50316)then
1063 print*,'MASK: view 7 - vk 21/22'
1064 mask_vk(7,21)=0
1065 mask_vk(7,22)=0
1066 do is=1,nstrips_va1
1067 mask(7,21,is)=0
1068 mask(7,22,is)=0
1069 enddo
1070 endif
1071
1072 * ---------------------
1073 * VIEW 12 - VK 1-2-3-4
1074 if((date.eq.50317).and.(id.le.3))then
1075 print*,'MASK: view 12 - vk 1/2/3/4'
1076 mask_vk(12,1)=0
1077 mask_vk(12,2)=0
1078 mask_vk(12,3)=0
1079 mask_vk(12,4)=0
1080 do is=1,nstrips_va1
1081 mask(12,1,is)=0
1082 mask(12,2,is)=0
1083 mask(12,3,is)=0
1084 mask(12,4,is)=0
1085 enddo
1086 endif
1087
1088 * ---------------------
1089 * VIEW 7 - VK 5-6
1090 if(date.ge.50320)then
1091 if(.not.(date.eq.50320.and.id.le.3)) then
1092 print*,'MASK: view 7 - vk 5/6'
1093 mask_vk(7,5)=0
1094 mask_vk(7,6)=0
1095 do is=1,nstrips_va1
1096 mask(7,5,is)=0
1097 mask(7,6,is)=0
1098 enddo
1099 endif
1100 endif
1101
1102 * ---------------------
1103 * VIEW 7 - VK 13-14
1104 if(date.ge.50320)then
1105 if(.not.(date.eq.50320.and.id.le.5)) then
1106 print*,'MASK: view 7 - vk 13/14'
1107 mask_vk(7,13)=0
1108 mask_vk(7,14)=0
1109 do is=1,nstrips_va1
1110 mask(7,13,is)=0
1111 mask(7,14,is)=0
1112 enddo
1113 endif
1114 endif
1115
1116 *** SAMARA
1117 *** SAMARA
1118 *** SAMARA
1119 * it needs further checks...
1120
1121 * ---------------------
1122 * VIEW 7 - VK 9-10
1123 * VIEW 12 - VK 1-2-3-4
1124 if((date.eq.50516).and.(id.le.8))then
1125 print*,'MASK: view 7 - vk 9/10'
1126 print*,'MASK: view 12 - vk 1/2/3/4'
1127 mask_vk(7,9)=0
1128 mask_vk(7,10)=0
1129 mask_vk(12,1)=0
1130 mask_vk(12,2)=0
1131 mask_vk(12,3)=0
1132 mask_vk(12,4)=0
1133 do is=1,nstrips_va1
1134 mask(7,9,is)=0
1135 mask(7,10,is)=0
1136 mask(12,1,is)=0
1137 mask(12,2,is)=0
1138 mask(12,3,is)=0
1139 mask(12,4,is)=0
1140 enddo
1141 endif
1142
1143 * ---------------------
1144 * VIEW 7 - VK 9-10
1145 if(date.ge.50516)then
1146 if(.not.(date.eq.50516.and.id.le.8)) then
1147 print*,'MASK: view 7 - vk 9/10'
1148 mask_vk(7,9)=0
1149 mask_vk(7,10)=0
1150 do is=1,nstrips_va1
1151 mask(7,9,is)=0
1152 mask(7,10,is)=0
1153 enddo
1154 endif
1155 endif
1156
1157 * ---------------------
1158 * VIEW 12 - VK 7-8
1159 if(date.ge.50523)then
1160 if(.not.(date.eq.50523.and.id.le.3)) then
1161 print*,'MASK: view 12 - vk 7/8'
1162 mask_vk(12,7)=0
1163 mask_vk(12,8)=0
1164 do is=1,nstrips_va1
1165 mask(12,7,is)=0
1166 mask(12,8,is)=0
1167 enddo
1168 endif
1169 endif
1170
1171 return
1172 end

  ViewVC Help
Powered by ViewVC 1.1.23