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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1.1.1 - (hide annotations) (download) (vendor branch)
Wed Mar 8 15:00:39 2006 UTC (18 years, 9 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 pam-fi 1.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