/[PAMELA software]/DarthVader/TrackerLevel2/src/F77/reductionflight.f
ViewVC logotype

Diff of /DarthVader/TrackerLevel2/src/F77/reductionflight.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.4 by pam-fi, Fri Aug 4 08:18:06 2006 UTC revision 1.20 by pam-fi, Thu May 24 16:45:48 2007 UTC
# Line 18  Line 18 
18        include 'common_reduction.f'        include 'common_reduction.f'
19        include 'calib.f'        include 'calib.f'
20                
21          data eventn_old/nviews*0/
22    
23        integer ierror        integer ierror
24        ierror = 0        ierror = 0
25    
26  *     -------------------------------------------------------  c$$$      debug = .true.
27  *     STRIP MASK  c$$$      verbose = .true.
28  *     -------------------------------------------------------  c$$$      warning = .true.
29    
30    *     //////////////////////////
31    *     initialize some parameters
32    *     //////////////////////////
33    
 c      call stripmask   !called later, after CN computation  
34        call init_level1        call init_level1
35    
36        good1 = good0  c      debug=.true.
37    
38          if(debug)print*,'-- check LEVEL0 status'
39    
40          ievco=-1
41          mismatch=0
42    c      good1 = good0
43    c--------------------------------------------------
44    c     check the LEVEL0 event status for missing
45    c     sections or DSP alarms
46    c     ==> set the variable GOOD1(12)
47    c--------------------------------------------------
48          do iv=1,nviews
49             if(DSPnumber(iv).gt.0.and.DSPnumber(iv).le.12)then
50    c           ------------------------
51    c           GOOD
52    c           ------------------------
53                GOOD1(DSPnumber(iv))=0 !OK
54    c           ------------------------
55    c           CRC error
56    c           ------------------------
57                if(crc(iv).eq.1) then
58    c               GOOD1(DSPnumber(iv)) = 2
59    c               GOOD1(DSPnumber(iv)) = GOOD1(DSPnumber(iv)) + 2**1
60                   GOOD1(DSPnumber(iv)) = ior(GOOD1(DSPnumber(iv)),2**1)
61     102           format(' * WARNING * Event ',i7,' view',i3
62         $          ,' CRC error')
63                   if(debug)write(*,102)eventn(1),DSPnumber(iv)
64    c               goto 18 !next view
65                endif
66    c           ------------------------
67    c           online-software alarm
68    c           ------------------------
69                if(
70         $           fl1(iv).ne.0.or.
71         $           fl2(iv).ne.0.or.
72         $           fl3(iv).ne.0.or.
73         $           fl4(iv).ne.0.or.
74         $           fl5(iv).ne.0.or.
75         $           fl6(iv).ne.0.or.
76         $           fc(iv).ne.0.or.
77         $           DATAlength(iv).eq.0.or.
78         $           .false.)then
79    c               GOOD1(DSPnumber(iv))=3
80    c               GOOD1(DSPnumber(iv)) = GOOD1(DSPnumber(iv)) + 2**2
81                   GOOD1(DSPnumber(iv)) = ior(GOOD1(DSPnumber(iv)),2**2)
82     103           format(' * WARNING * Event ',i7,' view',i3
83         $          ,' software alarm')
84                   if(debug)write(*,103)eventn(1),DSPnumber(iv)
85    c               goto 18
86                endif
87    c           ------------------------
88    c           DSP-counter jump
89    c           ------------------------
90    c     commentato perche` non e` un controllo significativo nel caso in cui
91    c     la subroutine venga chiamata per riprocessare l'evento
92    c     sostituito con un check dei contatori dei vari dsp
93    c$$$            if(
94    c$$$     $           eventn_old(iv).ne.0.and. !first event in this file
95    c$$$     $           eventn(iv).ne.1.and.     !first event in run
96    c$$$     $           good_old(DSPnumber(iv)).ne.0.and. !previous event corrupted
97    c$$$     $           .true.)then
98    c$$$
99    c$$$               if(eventn(iv).ne.(eventn_old(iv)+1))then
100    c$$$c                  GOOD1(DSPnumber(iv))=4
101    c$$$c                  GOOD1(DSPnumber(iv)) = GOOD1(DSPnumber(iv)) + 2**3
102    c$$$                  GOOD1(DSPnumber(iv)) = ior(GOOD1(DSPnumber(iv)),2**3)
103    c$$$ 104              format(' * WARNING * Event ',i7,' view',i3
104    c$$$     $          ,' counter jump ',i10,i10)
105    c$$$                  if(debug)write(*,104)eventn(1),DSPnumber(iv)
106    c$$$     $                 ,eventn_old(iv),eventn(iv))
107    c$$$                  goto 18
108    c$$$               endif
109    c$$$
110    c$$$            endif
111    c           ------------------------
112    c 18         continue
113    c           ------------------------
114    c           DSP-counter
115    c           ------------------------
116                if( DSPnumber(iv).ne.0.and.GOOD1(DSPnumber(iv)).ne.1)then
117                   if(iv.ne.1.and.ievco.ne.-1)then
118                      if( eventn(iv).ne.ievco )then
119                         mismatch=1
120                      endif
121                   endif
122                   ievco = eventn(iv)
123                endif
124             endif
125          enddo
126    
127    c      print*,'*** ',(eventn(iv),iv=1,12)
128          
129          if(mismatch.eq.1.and.debug)
130         $     print*,' * WARNING * DSP counter mismatch: '
131         $     ,(eventn(iv),iv=1,12)
132    
133          ngood = 0
134          do iv = 1,nviews
135            
136             if(mismatch.eq.1.and.GOOD1(iv).ne.1)
137         $        GOOD1(iv)=ior(GOOD1(iv),2**3)
138    
139             eventn_old(iv) = eventn(iv)
140             good_old(iv)   = good1(iv)
141             ngood = ngood + good1(iv)
142    
143          enddo
144    c$$$      if(verbose.and.ngood.ne.0)print*,'* WARNING * Event ',eventn(1)
145    c$$$     $     ,':LEVEL0 event status: '
146    c$$$     $     ,(good1(i),i=1,nviews)
147  c--------------------------------------------------  c--------------------------------------------------
148  c     read the variable DATATRACKER from LEVEL0  c     read the variable DATATRACKER from LEVEL0
149  c     and fill the variable ADC (inverting view 11)  c     and fill the variable ADC (invertin view 11)
150  c--------------------------------------------------  c--------------------------------------------------
151          
152          if(debug)print*,'-- fill ADC vectors'
153    
154        call filladc(iflag)        call filladc(iflag)
155        if(iflag.ne.0)then        if(iflag.ne.0)then
156          good1=0           ierror = 220
 c       if(DEBUG)print*,'event ',eventn(1),' >>>>>  decode ERROR'  
             ierror = 220  
         goto 200  
157        endif        endif
158    
159  c--------------------------------------------------  c--------------------------------------------------
160  c     computes common noise for each VA1  c     computes common noise for each VA1
161  c     (excluding strips affected by signal,  c     (excluding strips with signal,
162  c     tagged with the flag CLSTR)  c     tagged with the flag CLSTR)
163  c--------------------------------------------------  c--------------------------------------------------
164          if(debug)print*,'-- compute CN'
165    
166        do iv=1,nviews        do iv=1,nviews
167          do ik=1,nva1_view           ima=0
168            cn(iv,ik)=0                     do ik=1,nva1_view
169            mask_vk_ev(iv,ik)=1              cn(iv,ik)    = 0
170            iflag=0              cnrms(iv,ik) = 0
171            if(mask_vk(iv,ik).eq.1)call cncomp(iv,ik,iflag)              cnn(iv,ik)   = -1
172  c          if(iflag.ne.0)good1=0              iflag = 0
173            if(iflag.ne.0)then              mask_vk_ev(iv,ik) = 1
174               mask_vk_ev(iv,ik)=0              call stripmask(iv,ik)      !compute mask(i,j,k), combining VA1-masks
175               ierror = 220  *           --------------------------------------
176            endif  *           if chip is not masked ---> evaluate CN
177          enddo  *           --------------------------------------
178                if( mask(iv,ik,1).eq.1 ) then !!!NBNB mask per la striscia 1 !!!!!!!!
179                   call cncomp(iv,ik,iflag)
180                   if(iflag.ne.0)then
181                      ima=ima+1
182                      mask_vk_ev(iv,ik)=0
183                      ierror = 220
184                   endif
185                   call stripmask(iv,ik) !compute mask(i,j,k), combining VA1-masks
186                endif
187             enddo
188     100     format(' * WARNING * Event ',i7,' view',i3,': VK MASK ',24i1)
189             if(ima.ne.0.and.verbose)write(*,100)eventn(1),iv
190         $        ,(mask_vk_ev(iv,ik),ik=1,nva1_view)
191    c         if(ima.ne.0)write(*,100)eventn(1),iv
192    c     $        ,(mask_vk_ev(iv,ik),ik=1,nva1_view)  
193        enddo        enddo
 c      if(good1.eq.0)then  
 c         ierror = 220  
 c      endif  
194    
195        call stripmask !compute mask(i,j,k)  cc      call stripmask !compute mask(i,j,k), combining mask_vk_ev and mask_vk
196    
197  c---------------------------------------------  c---------------------------------------------
198  c     loops on views, VA1 and strips,  c     loops on views, VA1 and strips,
199  c     and computes strips signals using  c     and computes strips signals using
200  c     badstrip, pedestals, and  c     badstrip, pedestals, and
201  c     sigma informations from histograms  c     sigma informations from histograms
202  c---------------------------------------------  c---------------------------------------------
       flag_shower = .false.  
203        ind=1                     !clsignal array index        ind=1                     !clsignal array index
204    
205          if(debug)print*,'-- search clusters'
206        do iv=1,nviews            !loop on views        do iv=1,nviews            !loop on views
207          do is=1,nstrips_view    !loop on strips (1)          do is=1,nstrips_view    !loop on strips (1)
208            if(mod(iv,2).eq.1) then            if(mod(iv,2).eq.1) then
209  C===  > Y view  C===  > Y view
210    c             print*,iv,nvk(is),nst(is),adc(iv,nvk(is),nst(is))
211    c     $            ,cn(iv,nvk(is))
212    c     $            ,pedestal(iv,nvk(is),nst(is))
213              value(is)= -(DBLE(adc(iv,nvk(is),nst(is)))              value(is)= -(DBLE(adc(iv,nvk(is),nst(is)))
214       $           -cn(iv,nvk(is))-pedestal(iv,nvk(is),nst(is)))       $           -cn(iv,nvk(is))-pedestal(iv,nvk(is),nst(is)))
215       $           *mask(iv,nvk(is),nst(is))       $           *mask(iv,nvk(is),nst(is))
# Line 83  C===  > Y view Line 217  C===  > Y view
217       $           *mask(iv,nvk(is),nst(is))       $           *mask(iv,nvk(is),nst(is))
218              clinclcut(is)=incuty*sigma(iv,nvk(is),nst(is))              clinclcut(is)=incuty*sigma(iv,nvk(is),nst(is))
219       $           *mask(iv,nvk(is),nst(is))       $           *mask(iv,nvk(is),nst(is))
220  ccc            print*,"value(",is,")(reduction)= ",value(is)              sat(is)=0
221                if( adc(iv,nvk(is),nst(is)).lt.adc_saty )sat(is)=1
222            else                        else            
223  C===  > X view  C===  > X view
224              value(is)= (DBLE(adc(iv,nvk(is),nst(is)))              value(is)= (DBLE(adc(iv,nvk(is),nst(is)))
# Line 93  C===  > X view Line 228  C===  > X view
228       $           *mask(iv,nvk(is),nst(is))       $           *mask(iv,nvk(is),nst(is))
229              clinclcut(is)=incutx*sigma(iv,nvk(is),nst(is))              clinclcut(is)=incutx*sigma(iv,nvk(is),nst(is))
230       $           *mask(iv,nvk(is),nst(is))       $           *mask(iv,nvk(is),nst(is))
231                sat(is)=0
232                if( adc(iv,nvk(is),nst(is)).gt.adc_satx )sat(is)=1
233            endif            endif
 c$$$          print*,iv,is,' --- ',adc(iv,nvk(is),nst(is)),cn(iv,nvk(is))  
 c$$$     $         ,pedestal(iv,nvk(is),nst(is)),value(is)  
 c$$$     $         ,sigma(iv,nvk(is),nst(is))  
 c          if(value(is).gt.clseedcut(is))  
 c     $         print*,iv,is,' --- (ADC_PED_CN) ',value(is),clseedcut(is)  
234          enddo                   !end loop on strips (1)          enddo                   !end loop on strips (1)
235          call search_cluster(iv)          call search_cluster(iv)
236          if(flag_shower.eqv..true.)then  
237            call init_level1                        if(.not.flag_shower)then
238            good1=0             call save_cluster(iv)
239            goto 200              !jump to next event             if(debug)print*,'view ',iv,' #clusters ', nclstr_view
240            else
241               fshower(iv) = 1
242    c           GOOD1(DSPnumber(iv)) = 11 !AHAHAHAHA ORRORE!!
243    c           GOOD1(iv) = 11
244    c           GOOD1(iv) = GOOD1(iv) + 2**5
245               GOOD1(iv) = ior(GOOD1(iv),2**5)
246     101       format(' * WARNING * Event ',i7,' view',i3
247         $          ,' #clusters > ',i5,' --> MASKED')
248               if(verbose)write(*,101)eventn(1),iv,nclstrmax_view
249          endif          endif
250        enddo                     ! end loop on views        enddo                     ! end loop on views
251        do iv=1,nviews        do iv=1,nviews
252          do ik=1,nva1_view          do ik=1,nva1_view
253            cnev(iv,ik)=cn(iv,ik) !assigns computed CN to ntuple variables            cnev(iv,ik)    = cn(iv,ik) !assigns computed CN to ntuple variables
254            cnevflag(iv,ik)=cnflag(iv,ik) !assigns computed CN to ntuple variables            cnrmsev(iv,ik) = cnrms(iv,ik) !assigns computed CN to ntuple variables
255  ccc          print*,"cnev(",iv,",",ik,")(reduction)= ",cnev(iv,ik)            cnnev(iv,ik)   = cnn(iv,ik) !assigns computed CN to ntuple variables
256          enddo          enddo
257        enddo        enddo
258  C---------------------------------------------  C---------------------------------------------
# Line 119  C     come here if GOOD1=0 Line 260  C     come here if GOOD1=0
260  C     or the event has too many clusters  C     or the event has too many clusters
261  C---------------------------------------------  C---------------------------------------------
262   200  continue   200  continue
263    
264          ngood = 0
265          do iv = 1,nviews
266             ngood = ngood + good1(iv)
267          enddo
268          if(verbose.and.ngood.ne.0)print*,'* WARNING * Event ',eventn(1)
269         $     ,':LEVEL1 event status: '
270         $     ,(good1(i),i=1,nviews)
271  c------------------------------------------------------------------------  c------------------------------------------------------------------------
272  c  c
273  c     closes files and exits  c     closes files and exits
# Line 146  c--------------------------------------- Line 295  c---------------------------------------
295        include 'level1.f'        include 'level1.f'
296        include 'level0.f'        include 'level0.f'
297    
298        good1=0  c      good1 = 0
299        nclstr1=0        do iv=1,12
300        totCLlength=0           good1(iv) = 1 !missing packet
301          enddo
302          nclstr1 = 0
303          totCLlength = 0
304        do ic=1,nclstrmax        do ic=1,nclstrmax
305           view(ic)=0           view(ic) = 0
306           ladder(ic)=0           ladder(ic) = 0
307           indstart(ic)=0           indstart(ic) = 0
308           indmax(ic)=0           indmax(ic) = 0
309           maxs(ic)=0           maxs(ic) = 0
310           mult(ic)=0                     mult(ic) = 0          
311           dedx(ic)=0           sgnl(ic) = 0
312             whichtrack(ic) = 0     !assigned @ level2
313    
314        enddo        enddo
315        do id=1,maxlength         !???        do id=1,maxlength         !???
316           clsignal(id)=0.           clsignal(id) = 0.
317             clsigma(id)  = 0.
318             cladc(id)    = 0.
319             clbad(id)    = 0.
320        enddo        enddo
321        do iv=1,nviews        do iv=1,nviews
322  c        crc1(iv)=0  c        crc1(iv)=0
323          do ik=1,nva1_view          do ik=1,nva1_view
324            cnev(iv,ik)=0            cnev(iv,ik) = 0
325              cnnev(iv,ik) = 0
326          enddo          enddo
327            fshower(iv) = 0
328        enddo        enddo
329                
330        return        return
331        end        end
332    
333  *---***---***---***---***---***---***---***---***  *---***---***---***---***---***---***---***---***
334  *  *
335  *  *
# Line 181  c        crc1(iv)=0 Line 341  c        crc1(iv)=0
341        subroutine search_cluster(iv)        subroutine search_cluster(iv)
342    
343        include 'commontracker.f'        include 'commontracker.f'
       include 'common_reduction.f'  
344        include 'level0.f'        include 'level0.f'
345        include 'level1.f'        include 'level1.f'
346        include 'calib.f'        include 'calib.f'
347    
348          include 'common_reduction.f'
349            
350    
351  c     local variables  c     local variables
352        integer rmax,lmax         !estremi del cluster        integer rmax,lmax         !estremi del cluster
353        integer rstop,lstop       !per decidere quali strip includere nel cluster        integer rstop,lstop      
354                                  ! oltre il seed        integer first,last
355        integer first,last,diff   !per includere le strip giuste... !???        integer fsat,lsat
356    
357        integer multtemp          !temporary multiplicity variable        external nst
358    
359        integer CLlength          !lunghezza in strip del cluster        iseed=-999                !cluster seed index initialization
360    
361        external nst        inext=-999                !index where to start new cluster search
362    
363  c------------------------------------------------------------------------        flag_shower = .false.
364  c     looks for clusters on each view        nclstr_view=0
 C     : CERCO STRIP SOPRA CLSEEDCUT, POI SCORRO A DX FINCHE'  
 c     NON TROVO  
 C     STRIP PIU' BASSA (in segnale/rumore)  
 C     => L'ULTIMA DELLA SERIE CRESCENTE  
 C     (LA PIU' ALTA) E' IL  
 C     CLUSTER SEED. POI SCORRO A SX E DX INCLUDENDO TUTTE  
 C     LE STRIP (FINO A 17 AL  
 C     MAX) CHE SUPERANO CLINCLCUT.  
 C     QUANDO CERCO IL CLUSTER SEED SUCCESSIVO SALTO LA STRIP  
 C     ADIACENTE A DESTRA  
 C     DELL'ULTIMO CLUSTER SEED (CHE SARA' NECESSARIAMENTE  
 C     PIU' BASSA) E PRENDO  
 C     COME SEED UNA STRIP SOLO SE IL SUO SEGNALE E'  
 C     MAGGIORE DI QUELLO DELLA STRIP  
 C     PRECEDENTE (PRATICAMENTE PER EVITARE CHE L'ULTIMA  
 C     STRIP DI UN GRUPPO DI STRIP  
 C     TUTTE SOPRA IL CLSEEDCUT VENGA AUTOMATICAMENTE PRESA  
 C     COME SEED... DEVE ESSERE  
 C     PRESA SOLO SE IL CLUSTER E' DOUBLE PEAKED...)  
 c------------------------------------------------------------------------  
 c     6 ottobre 2003  
 c     Elena: CLSEEDCUT = 7 (old value 10)  
 c     Elena: CLINCLCUT = 4 (old value 5)  
365    
366        iseed=-999                !cluster seed index initialization        do jl=1,nladders_view              !1..3 !loops on ladders
367    
368        do jl=1,nladders_view     !1..3 !loops on ladders           first = 1+nstrips_ladder*(jl-1) !1,1025,2049
369           first=1+nstrips_ladder*(jl-1) !1,1025,2049           last  = nstrips_ladder*jl       !1024,2048,3072
370           last=nstrips_ladder*jl !1024,2048,3072  
371  c     X views have 1018 strips instead of 1024  *        X views have 1018 strips instead of 1024
372           if(mod(iv,2).eq.0) then           if(mod(iv,2).eq.0) then
373              first=first+3              first=first+3
374              last=last-3              last=last-3
375           endif           endif
376    
377           do is=first,last       !loop on strips in each ladder           do is=first,last       !loop on strips in each ladder
378              if(is.le.iseed+1) goto 220  
379  c-----------------------------------------  c---------------------------------------------
380  c     after a cluster seed as been found,  c     new-cluster search starts at index inext
381  c     look for next one skipping one strip on the right  c---------------------------------------------
382  c     (i.e. look for double peak cluster)              if(is.lt.inext) goto 220 ! next strip
383  c-----------------------------------------  
             if(is.ne.first) then  
                if(value(is).le.value(is-1)) goto 220  
             endif  
 c-----------------------------------------  
 c     skips cluster seed  
 c     finding if strips values are descreasing (a strip  
 c     can be a cluster seed only if previous strip value  
 c     is lower)  
 c-----------------------------------------  
384              if(value(is).gt.clseedcut(is)) then              if(value(is).gt.clseedcut(is)) then
 ccc              print*,"value(",is,")=",value(is),  
 ccc     $             " .gt.clseedcut(",is,")=",clseedcut(is)  
385  c-----------------------------------------  c-----------------------------------------
386  c     possible SEED...  c     possible SEED...
387  c-----------------------------------------  c-----------------------------------------
388                 itemp=is                 itemp = is
389                   fsat = 0         ! first saturated strip
390                   lsat = 0         ! last saturated strip
391                 if(itemp.eq.last) goto 230 !estremo...                 if(itemp.eq.last) goto 230 !estremo...
392                 do while(value(itemp)  c              ------------------------                
393       $              /sigma(iv,nvk(itemp),nst(itemp))  c              search for first maximum
394       $              .le.value(itemp+1)  c              ------------------------
395       $              /sigma(iv,nvk(itemp+1),nst(itemp+1))) !BIAS: aggiustare il caso uguale!???                 do while(
396                    itemp=itemp+1       $                   value(itemp).le.value(itemp+1)
397                    if(itemp.eq.last) goto 230 !stops if reaches last strip       $              .and.value(itemp+1).gt.clseedcut(itemp+1))
398                      itemp = itemp+1
399                      if(itemp.eq.last)   goto 230 !stops if reaches last strip
400                      if(sat(itemp).eq.1) goto 230 !stop if reaches a saturated strip
401                 enddo            ! of the ladder                 enddo            ! of the ladder
402   230           continue   230           continue
403  c-----------------------------------------  c              -----------------------------
404    c              check if strips are saturated
405    c              -----------------------------    
406                   if( sat(itemp).eq.1 )then
407                      fsat = itemp
408                      lsat = itemp
409                      if(itemp.eq.last) goto 231 !estremo...
410                      do while( sat(itemp+1).eq.1 )
411                         itemp = itemp+1
412                         lsat = itemp
413                         if(itemp.eq.last)   goto 231 !stops if reaches last strip
414                      enddo                  
415                   endif
416     231           continue
417    c---------------------------------------------------------------------------
418  c     fownd SEED!!!  c     fownd SEED!!!
419  c-----------------------------------------  c     (if there are saturated strips, the cluster is centered in the middle)
420                 iseed=itemp      c---------------------------------------------------------------------------
421  c----------------------------------------------------------                 if(fsat.eq.0.and.lsat.eq.0)then
422                      iseed = itemp ! <<< SEED
423                   else
424                      iseed = int((lsat+fsat)/2) ! <<< SEED
425    c$$$                  print*,'saturated strips ',fsat,lsat,iseed
426    c$$$                  print*,'--> ',(value(ii),ii=fsat,lsat)
427                   endif    
428    c---------------------------------------------------------------
429  c     after finding a cluster seed, checks also adjacent strips,  c     after finding a cluster seed, checks also adjacent strips,
430  C     and marks the ones exceeding clinclcut  C     and tags the ones exceeding clinclcut
431  c----------------------------------------------------------  c---------------------------------------------------------------
432                 ir=iseed         !indici destro                 ir=iseed         !indici destro
433                 il=iseed         ! e sinistro                 il=iseed         ! e sinistro
434                                
# Line 286  c--------------------------------------- Line 439  c---------------------------------------
439                 lstop=0          ! inclusion loop                 lstop=0          ! inclusion loop
440    
441                 do while(lstop.eq.0.or.rstop.eq.0) !shifts left and right from                 do while(lstop.eq.0.or.rstop.eq.0) !shifts left and right from
442                    ir=ir+1       !position index for strips on right side of  
443                                  ! cluster seed  
444                    il=il-1       !and for left side                    ir=ir+1       !index for right side
445                      il=il-1       !index for left side
446  c------------------------------------------------------------------------  c------------------------------------------------------------------------
447  c     checks for last or first strip of the ladder  c     checks for last or first strip of the ladder
448  c------------------------------------------------------------------------  c------------------------------------------------------------------------
449                    if(ir.gt.last) then !when index goes beyond last strip                    if( ir.gt.last  ) rstop = 1                      
450                       rstop=1    ! of the ladder, change rstop flag in order                    if( il.lt.first ) lstop = 1
                                 ! to "help" exiting from loop  
                   endif  
                     
                   if(il.lt.first) then !idem when index goes beyond  
                      lstop=1    ! first strip of the ladder  
                   endif  
451                                        
452  c------------------------------------------------------------------------  c------------------------------------------------------------------------
453  c     check for clusters including more than nclstrp strips  c     add strips exceeding inclusion cut
454  c------------------------------------------------------------------------  c------------------------------------------------------------------------
455                    if((rmax-lmax+1).ge.nclstrp) then                    if( (rmax-lmax+1).ge.nclstrp )goto 210   !exits inclusion loop
456                       goto 210   !exits inclusion loop:  
457                                  ! lmax and rmax maintain last value                    if(rstop.eq.0) then !if right cluster morder has not been reached
458                                  ! NB .ge.!???                       if(value(ir).gt.clinclcut(ir)) then
459                    endif                          rmax=ir !include a strip on the right
 c------------------------------------------------------------------------  
 c     marks strips exceeding inclusion cut  
 c------------------------------------------------------------------------  
                   if(rstop.eq.0) then !if last strip of the ladder or last  
                                 ! over-cut strip has not been reached  
                      if(value(ir).gt.clinclcut(ir)) then !puts in rmax the  
                         rmax=ir ! last right over-cut strip  
460                       else                       else
461                          rstop=1 !otherwise cluster ends on right and rstop                          rstop=1 !cluster right end
462                       endif      ! flag=1 signals it                       endif    
463                    endif                    endif
464                    if(lstop.eq.0) then  
465                      if( (rmax-lmax+1).ge.nclstrp )goto 210   !exits inclusion loop
466    
467                      if(lstop.eq.0) then !if left cluster morder has not been reached
468                       if(value(il).gt.clinclcut(il)) then                       if(value(il).gt.clinclcut(il)) then
469                          lmax=il                          lmax=il !include a strip on the left
470                       else                       else
471                          lstop=1                          lstop=1 !cluster left end
472                       endif                       endif
473                    endif                    endif
474    
475                 enddo            !ends strip inclusion loop                 enddo            !ends strip inclusion loop
476                   goto 211
477   210           continue         !jumps here if more than nclstrp have been included   210           continue         !jumps here if more than nclstrp have been included
478                        c               print*,'>>> nclstrp! '
479                 multtemp=rmax-lmax+1 !stores multiplicity in temp   211           continue
480                                  ! variable. NB rmax and lmax can change later in  c-----------------------------------------
481                                  ! order to include enough strips to calculate eta3  c     end of inclusion loop!
482                                  ! and eta4. so mult is not always equal to cllength  c-----------------------------------------
483  c------------------------------------------------------------------------                
484  c     NB per essere sicuro di poter calcolare eta3 e eta4 devo includere  c------------------------------------------------------------------------
485  c     sempre e comunque le 2 strip adiacenti al cluster seed e quella  c     adjust the cluster in order to have at least a strip around the seed
486  c     adiacente ulteriore dalla parte della piu' alta fra queste due  c------------------------------------------------------------------------
487  c     (vedi oltre...)!???                 if(iseed.eq.lmax.and.lmax.ne.first)then
488  c------------------------------------------------------------------------                    lmax = lmax-1
489                      if( (rmax-lmax+1).gt.nclstrp )rmax=rmax-1
490  c     nel caso di estremi del ladder...!???                 endif
491                   if(iseed.eq.rmax.and.rmax.ne.last )then
492  c     ho meno di 4 strip nel cluster --> se sono sui bordi o quasi del ladder                    rmax = rmax+1
493  c     costruisco il cluster ad hoc e poi esco, se non sono sui bordi o quasi                    if( (rmax-lmax+1).gt.nclstrp )lmax=lmax+1
 c     vado oltre (aggiungero' quindi strip a sx e dx in modo da poter calcolare  
 c     eta3e4)  
                if((rmax-lmax+1).lt.4) then  
   
                   if(iseed.eq.first) then !estremi...  
                      rmax=iseed+2 !NB in questo modo puo' anche capitare di  
                      lmax=iseed ! includere strip sotto taglio di inclusione  
                      goto 250   ! che non serviranno per eta3e4!???  
                   endif  
                     
                   if(iseed.eq.last) then !estremi...  
                      rmax=iseed  
                      lmax=iseed-2 !NB 2 e non 3, perche' altrimenti sarei in  
                      goto 250   ! ((rmax-lmax+1).lt.4).eq.false. !???  
                   endif         !NMB questo e' l'unico caso di cllength=3!???  
                     
                   if(iseed.eq.first+1) then !quasi estremi...  
                      rmax=iseed+2  
                      lmax=iseed-1  
                      goto 250  
                   endif  
                   if(iseed.eq.last-1) then  
                      rmax=iseed+1  
                      lmax=iseed-2  
                      goto 250  
                   endif  
 c     se ho 4 o piu' strip --> se sono sui bordi esco, se sono sui quasi bordi  
 c     includo la strip del bordo  
                else  
                     
                   if(iseed.eq.first) goto 250 !estremi... non includo altro                    
                   if(iseed.eq.last) goto 250  
                   if(iseed.eq.first+1) then !quasi estremi... mi assicuro di  
                      lmax=first ! avere le strip adiacenti al seed  
                      if((rmax-lmax+1).gt.nclstrp) rmax=rmax-1 !NB effetto  
                      goto 250   ! coperta: se la lunghezza del cluster era gia'  
                   endif         ! al limite (nclstrp), per poter aggiungere questa  
                                 ! strip a sinistra devo toglierne una a destra...!???  
                   if(iseed.eq.last-1) then  
                      rmax=last  
                      if((rmax-lmax+1).gt.nclstrp) lmax=lmax+1  
                      goto 250  
                   endif                    
494                 endif                 endif
495    
496  c------------------------------------------------------------------------  c------------------------------------------------------------------------
497  c     be sure to include in the cluster the cluster seed with its 2 adjacent  c     adjust the cluster in order to store a minimum number of strips
498  c     strips, and the one adjacent to the greatest between this two strip, as the  c------------------------------------------------------------------------
499  c     fourth one. if the strips have the same value (!) the fourth one is chosen                 do while( (rmax-lmax+1).lt.nclstrpmin )
500  c     as the one having the greatest value between the second neighbors  
501  c------------------------------------------------------------------------                    vl = -99999
502                 if(value(iseed+1).eq.value(iseed-1)) then                    vr = -99999
503                    if(value(iseed+2).ge.value(iseed-2)) then !??? qui cmq c'e'                    if(lmax-1.ge.first) vl = value(lmax-1)
504                       diff=(iseed+2)-rmax                    if(rmax+1.le.last ) vr = value(rmax+1)
505                       if(diff.gt.0) then                    if(vr.ge.vl) then
506                          rmax=rmax+diff                       rmax = rmax+1
507                          if((rmax-lmax+1).gt.nclstrp) then                    else  
508                             lmax=rmax-nclstrp+1                       lmax = lmax-1
                         endif  
                      endif  
                      diff=(iseed-1)-lmax  
                      if(diff.lt.0) then  
                         lmax=lmax+diff  
                         if((rmax-lmax+1).gt.nclstrp) then  
                            rmax=lmax+nclstrp-1  
                         endif  
                      endif  
                   else  
                      diff=(iseed-2)-lmax  
                      if(diff.lt.0) then  
                         lmax=lmax+diff  
                         if((rmax-lmax+1).gt.nclstrp) then  
                            rmax=lmax+nclstrp-1  
                         endif  
                      endif  
                      diff=(iseed+1)-rmax  
                      if(diff.gt.0) then  
                         rmax=rmax+diff  
                         if((rmax-lmax+1).gt.nclstrp) then  
                            lmax=rmax-nclstrp+1  
                         endif  
                      endif  
                   endif  
                elseif(value(iseed+1).gt.value(iseed-1)) then  
 c     !??? sposto il limite del cluster a destra per includere sempre le strip  
 c     necessarie al calcolo di eta-i  
 c     se il cluster diventa  troppo lungo lo accorcio a sinistra per avere non piu'  
 c     di nclstrp (in questo caso sono sicuro di aver gia' incluso le strip  
 c     necessarie al calcolo di eta-i a sinistra, quindi se voglio posso uscire)  
                   diff=(iseed+2)-rmax  
                   if(diff.gt.0) then  
                      rmax=rmax+diff  
                      if((rmax-lmax+1).gt.nclstrp) then  
                         lmax=rmax-nclstrp+1  
 c     goto 250  
                      endif  
                   endif  
                   diff=(iseed-1)-lmax  
                   if(diff.lt.0) then  
                      lmax=lmax+diff  
                      if((rmax-lmax+1).gt.nclstrp) then  
                         rmax=lmax+nclstrp-1  
 c     goto 250 !inutile!???  
                      endif  
                   endif  
                else  
                   diff=(iseed-2)-lmax  
                   if(diff.lt.0) then  
                      lmax=lmax+diff  
                      if((rmax-lmax+1).gt.nclstrp) then  
                         rmax=lmax+nclstrp-1  
 c     goto 250  
                      endif  
                   endif  
                   diff=(iseed+1)-rmax  
                   if(diff.gt.0) then  
                      rmax=rmax+diff  
                      if((rmax-lmax+1).gt.nclstrp) then  
                         lmax=rmax-nclstrp+1  
 c     goto 250 !inutile!???  
                      endif  
509                    endif                    endif
510                 endif                    
511   250           continue                 enddo
512    
513  c--------------------------------------------------------  c--------------------------------------------------------
514  c     fills cluster variables  c     store cluster info
515  c--------------------------------------------------------  c--------------------------------------------------------
516                 nclstr1=nclstr1+1 !cluster number                 nclstr_view = nclstr_view + 1 !cluster number
517  ccc               print*,nclstr1,multtemp  
518                 if(nclstr1.gt.nclstrmax) then !too many clusters for the event:                 if(nclstr_view.gt.nclstrmax_view) then !too many clusters for the view:
519                    if(verbose)print*,'Event ',eventn(1),  c$$$                  if(verbose) print*,'Event ',eventn(1),
520       $                 ': more than ',nclstrmax,' clusters'  c$$$     $                 ': more than ',nclstrmax_view
521                    good1=0       ! event  c$$$     $                 ,' clusters on view ',iv
                   nclstr1=0  
                   totCLlength=0  
522                    flag_shower = .true.                    flag_shower = .true.
523                    goto 2000                    goto 2000
524                 endif                 endif
525                 view(nclstr1)=iv !vista del cluster  
526                 ladder(nclstr1)=nld(iseed,iv) !ladder a cui appartiene il cluster seed                 ladder_view(nclstr_view) = nld(iseed,iv)
527                 maxs(nclstr1)=iseed !strip del cluster seed                 maxs_view(nclstr_view)   = iseed
528                 mult(nclstr1)=multtemp !molteplicita'                 mult_view(nclstr_view)   = rmax-lmax+1
529                                 rmax_view(nclstr_view)   = rmax
530                 indstart(nclstr1)=ind !posizione dell'inizio del cluster nell'                 lmax_view(nclstr_view)   = lmax
531                                  ! array clsignal  
532                 indmax(nclstr1)=indstart(nclstr1)+(iseed-lmax) !posizione del  c$$$               if(rmax-lmax+1.gt.25)
533                                  ! cluster seed nell'array clsignal  c$$$     $              print*,'view ',iv
534    c$$$     $              ,' cl ',nclstr_view,' mult ',rmax-lmax+1
535    c------------------------------------------------------------------------
536    c     search for a double peak inside the cluster                                                                                                            
537    c------------------------------------------------------------------------
538                   inext = rmax+1   !<< index where to start new-cluster search
539                                
540                 CLlength=rmax-lmax+1 !numero di strip del cluster                 vmax = 0
541                 totCLlength=totCLlength+CLlength                 vmin = value(iseed)
542                 dedx(nclstr1)=0                 imax = iseed
543                 do j=lmax,rmax   !stores sequentially cluter strip values in                 imin = iseed
544                    clsignal(ind)=value(j) ! clsignal array                 do iss = max(iseed+1,lsat+1),rmax
545                    ind=ind+1                    if( value(iss).lt.vmin )then
546  c                  if(value(j).gt.0)                       if( imax.ne.iseed )goto 221 !found dowble peek
547                    if(value(j).gt.clinclcut(j))                       imin = iss
548       $                 dedx(nclstr1)=dedx(nclstr1)+value(j) !cluster charge                       vmin = value(iss)
549                      else
550                         delta = value(iss) - value(imin)
551                         cut = sqrt(clinclcut(iss)**2 + clinclcut(imin)**2)
552                         if(
553         $                    delta.gt.cut .and.
554         $                    value(iss).gt.clseedcut(iss).and.
555         $                    .true.)then
556                            if( value(iss).gt.vmax )then                        
557                               imax = iss
558                               vmax = value(iss)
559                            else
560                               goto 221 !found dowble peek
561                            endif
562                         endif
563                      endif
564                 enddo                 enddo
565     221           continue
566                  
567                   if(imax.gt.iseed)then
568                      inext = imax    !<< index where to start new-cluster search
569    c$$$                  print*,'--- double peek ---'
570    c$$$                  print*,(value(ii),ii=lmax,rmax)
571    c$$$                  print*,'seed ',iseed,' imin ',imin,' imax ',imax
572                   endif
573  c--------------------------------------------------------  c--------------------------------------------------------
574  c  c
575  c--------------------------------------------------------  c--------------------------------------------------------
# Line 525  c--------------------------------------- Line 591  c---------------------------------------
591  *  *
592  *---***---***---***---***---***---***---***---***  *---***---***---***---***---***---***---***---***
593    
594          subroutine save_cluster(iv)
595    *
596    *     (080/2006 Elena Vannuccini)
597    *     Save the clusters view by view
598    
599          include 'commontracker.f'
600          include 'level1.f'
601          include 'calib.f'
602          include 'common_reduction.f'
603    
604          integer CLlength          !lunghezza in strip del cluster
605    
606          do ic=1,nclstr_view
607    
608             nclstr1 = nclstr1+1
609             view(nclstr1)   = iv
610             ladder(nclstr1) = ladder_view(ic)
611             maxs(nclstr1)   = maxs_view(ic)
612             mult(nclstr1)   = mult_view(ic)
613                  
614    c        posizione dell'inizio del cluster nell' array clsignal
615             indstart(nclstr1) = ind
616    c        posizione del cluster seed nell'array clsignal
617             indmax(nclstr1)   = indstart(nclstr1)
618         $        +( maxs_view(ic) - lmax_view(ic) )
619            
620             CLlength      = rmax_view(ic) - lmax_view(ic) + 1 !numero di strip salvate
621             totCLlength   = totCLlength + CLlength
622             sgnl(nclstr1) = 0
623             do j=lmax_view(ic),rmax_view(ic)         !stores sequentially cluter strip values in
624    
625                clsignal(ind) = value(j) ! clsignal array
626    c$$$            print*,ind,clsignal(ind)
627                ivk=nvk(j)
628                ist=nst(j)
629    
630                clsigma(ind) = sigma(iv,ivk,ist)
631                cladc(ind)   = adc(iv,ivk,ist)
632                clbad(ind)   = bad(iv,ivk,ist)
633    c            clped(ind)   = pedestal(iv,ivk,ist)
634    
635                ind=ind+1
636    c     if(value(j).gt.0)
637                if(value(j).gt.clinclcut(j))
638         $           sgnl(nclstr1) = sgnl(nclstr1) + value(j) !cluster charge
639             enddo
640    
641    c$$$         print*,'view ',iv,' -- save_cluster -- nclstr1: '
642    c$$$     $        ,nclstr1,maxs(nclstr1),mult(nclstr1),sgnl(nclstr1)
643    c$$$         print*,'----------------------'
644    
645          enddo
646          
647          return
648          end
649    *---***---***---***---***---***---***---***---***
650    *
651    *
652    *
653    *
654    *
655    *---***---***---***---***---***---***---***---***
656    
657    
658        subroutine stripmask  c$$$      subroutine stripmask
659    c$$$
660    c$$$*     this routine set va1 and single-strip masks,
661    c$$$*     on the basis of the VA1 mask saved in the DB
662    c$$$*
663    c$$$*     mask(nviews,nva1_view,nstrips_va1) !strip mask
664    c$$$*     mask_vk(nviews,nva1_view)          !VA1 mask
665    c$$$*
666    c$$$      include 'commontracker.f'
667    c$$$      include 'level1.f'
668    c$$$      include 'common_reduction.f'
669    c$$$      include 'calib.f'
670    c$$$
671    c$$$*     init mask
672    c$$$      do iv=1,nviews
673    c$$$         do ivk=1,nva1_view
674    c$$$            do is=1,nstrips_va1
675    c$$$c               mask(iv,ivk,is) = mask_vk(iv,ivk)
676    c$$$               if( mask_vk(iv,ivk) .ne. -1)then
677    c$$$                  mask(iv,ivk,is) = 1
678    c$$$     $                 * mask_vk(iv,ivk)     !from DB
679    c$$$     $                 * mask_vk_ev(iv,ivk)  !from <SIG>
680    c$$$     $                 * mask_vk_run(iv,ivk) !from CN
681    c$$$               else
682    c$$$                  mask(iv,ivk,is) = -1
683    c$$$     $                 * mask_vk(iv,ivk)     !from DB
684    c$$$     $                 * mask_vk_ev(iv,ivk)  !from CN
685    c$$$               endif
686    c$$$            enddo
687    c$$$         enddo
688    c$$$      enddo
689    c$$$
690    c$$$
691    c$$$      return
692    c$$$      end
693    
694          subroutine stripmask(iv,ivk)
695    
696    *     -----------------------------------------------
697  *     this routine set va1 and single-strip masks,  *     this routine set va1 and single-strip masks,
698  *     on the basis of the VA1 mask saved in the DB  *     on the basis of the VA1 mask saved in the DB
699  *  *
700  *     mask(nviews,nva1_view,nstrips_va1) !strip mask  *     mask(nviews,nva1_view,nstrips_va1) !strip mask
701  *     mask_vk(nviews,nva1_view)          !VA1 mask  *     mask_vk(nviews,nva1_view)          !VA1 mask
702  *  *     -----------------------------------------------
703        include 'commontracker.f'        include 'commontracker.f'
704  c      include 'level1.f'        include 'level1.f'
705        include 'common_reduction.f'        include 'common_reduction.f'
706        include 'calib.f'        include 'calib.f'
707    
708  *     init mask  *     init mask
709        do iv=1,nviews        do is=1,nstrips_va1
710           do ivk=1,nva1_view  *        --------------------------------------------------------
711              do is=1,nstrips_va1  *        if VA1-mask from DB is 0 or 1, three masks are combined:
712  c               mask(iv,ivk,is) = mask_vk(iv,ivk)  *        - from DB (a-priori mask)
713                 mask(iv,ivk,is) = mask_vk(iv,ivk) * mask_vk_ev(iv,ivk)                *        - run-based (chip declared bad on the basis of <SIG>)
714              enddo  *        - event-based (failure in CN computation)
715           enddo  *        --------------------------------------------------------
716    c         print*,iv,ivk
717    c     $        ,mask_vk(iv,ivk),mask_vk_ev(iv,ivk),mask_vk_run(iv,ivk)
718             if( mask_vk(iv,ivk) .ne. -1)then            
719                mask(iv,ivk,is) = 1
720         $           * mask_vk(iv,ivk)     !from DB
721         $           * mask_vk_ev(iv,ivk)  !from <SIG>
722         $           * mask_vk_run(iv,ivk) !from CN
723    *        -----------------------------------------------------------
724    *        if VA1-mask from DB is -1 only event-based mask is applied
725    *        -----------------------------------------------------------
726             else
727                mask(iv,ivk,is) = -1
728         $           * mask_vk(iv,ivk)     !from DB
729         $           * mask_vk_ev(iv,ivk)  !from CN
730             endif
731        enddo        enddo
732          
733          
734        return        return
735        end        end
   

Legend:
Removed from v.1.4  
changed lines
  Added in v.1.20

  ViewVC Help
Powered by ViewVC 1.1.23