/[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.9 by pam-fi, Mon Oct 16 12:36:52 2006 UTC revision 1.21 by pam-fi, Tue Aug 7 13:56:29 2007 UTC
# Line 23  Line 23 
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    c      debug=.true.
37    
38          if(debug)print*,'-- check LEVEL0 status'
39    
40          ievco=-1
41          mismatch=0
42  c      good1 = good0  c      good1 = good0
43  c--------------------------------------------------  c--------------------------------------------------
44  c     check the LEVEL0 event status for missing  c     check the LEVEL0 event status for missing
# Line 46  c           ------------------------ Line 55  c           ------------------------
55  c           CRC error  c           CRC error
56  c           ------------------------  c           ------------------------
57              if(crc(iv).eq.1) then              if(crc(iv).eq.1) then
58                 GOOD1(DSPnumber(iv)) = 2  c               GOOD1(DSPnumber(iv)) = 2
59                 goto 18 !next view  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              endif
66  c           ------------------------  c           ------------------------
67  c           online-software alarm  c           online-software alarm
# Line 62  c           ------------------------ Line 76  c           ------------------------
76       $           fc(iv).ne.0.or.       $           fc(iv).ne.0.or.
77       $           DATAlength(iv).eq.0.or.       $           DATAlength(iv).eq.0.or.
78       $           .false.)then       $           .false.)then
79                 GOOD1(DSPnumber(iv))=3  c               GOOD1(DSPnumber(iv))=3
80                 goto 18  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              endif
87  c           ------------------------  c           ------------------------
88  c           DSP-counter jump  c           DSP-counter jump
89  c           ------------------------  c           ------------------------
90              if(  c     commentato perche` non e` un controllo significativo nel caso in cui
91       $           eventn_old(iv).ne.0.and. !first event in this file  c     la subroutine venga chiamata per riprocessare l'evento
92       $           eventn(iv).ne.1.and.     !first event in run  c     sostituito con un check dei contatori dei vari dsp
93       $           good_old(DSPnumber(iv)).ne.0.and. !previous event corrupted  c$$$            if(
94       $           .true.)then  c$$$     $           eventn_old(iv).ne.0.and. !first event in this file
95    c$$$     $           eventn(iv).ne.1.and.     !first event in run
96                 if(eventn(iv).ne.(eventn_old(iv)+1))then  c$$$     $           good_old(DSPnumber(iv)).ne.0.and. !previous event corrupted
97                    GOOD1(DSPnumber(iv))=4  c$$$     $           .true.)then
98                    goto 18  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                 endif
122                   ievco = eventn(iv)
123              endif              endif
 c           ------------------------  
  18         continue  
124           endif           endif
125        enddo        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        ngood = 0
134        do iv = 1,nviews        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)           eventn_old(iv) = eventn(iv)
140           good_old(iv)   = good1(iv)           good_old(iv)   = good1(iv)
141           ngood = ngood + good1(iv)           ngood = ngood + good1(iv)
142    
143        enddo        enddo
144  c      if(ngood.ne.0)print*,'* WARNING * LEVEL0 event status: '  c$$$      if(verbose.and.ngood.ne.0)print*,'* WARNING * Event ',eventn(1)
145  c     $     ,(good1(i),i=1,nviews)  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 (invertin 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
 c        good1=0!<<<<<<<<<<<<<<<  
 c       if(DEBUG)print*,'event ',eventn(1),' >>>>>  decode ERROR'  
156           ierror = 220           ierror = 220
 c        goto 200  
 c         print*,'filladc error'  
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           ima=0           ima=0
168           do ik=1,nva1_view           do ik=1,nva1_view
169              cn(iv,ik)  = 0              cn(iv,ik)    = 0
170              cnn(iv,ik) = -1              cnrms(iv,ik) = 0
171              mask_vk_ev(iv,ik)=1              cnn(iv,ik)   = -1
172              iflag=0              iflag = 0
173              if(mask_vk(iv,ik).eq.1)call cncomp(iv,ik,iflag)              mask_vk_ev(iv,ik) = 1
174  c     if(iflag.ne.0)good1=0              call stripmask(iv,ik)      !compute mask(i,j,k), combining VA1-masks
175              if(iflag.ne.0)then  *           --------------------------------------
176                 ima=ima+1  *           if chip is not masked ---> evaluate CN
177                 mask_vk_ev(iv,ik)=0  *           --------------------------------------
178                 ierror = 220              if( mask(iv,ik,1).eq.1 ) then !!!NBNB mask per la striscia 1 !!!!!!!!
179  c$$$               if(verbose)                 call cncomp(iv,ik,iflag)
180  c$$$     $              print*,' * WARNING * Event ',eventn(1)                 if(iflag.ne.0)then
181  c$$$     $              ,': masked vk ',ik,' on view',iv                    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              endif
187           enddo           enddo
188   100     format(' * WARNING * Event ',i7,' view',i3,': VK MASK ',24i1)   100     format(' * WARNING * Event ',i7,' view',i3,': VK MASK ',24i1)
 c         if(ima.ne.0.and.verbose)print*,' * WARNING * Event ',eventn(1)  
 c     $              ,' view',iv,': VK MASK '  
 c     $        ,(mask_vk_ev(iv,ik),ik=1,nva1_view)  
189           if(ima.ne.0.and.verbose)write(*,100)eventn(1),iv           if(ima.ne.0.and.verbose)write(*,100)eventn(1),iv
190       $        ,(mask_vk_ev(iv,ik),ik=1,nva1_view)       $        ,(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), combining mask_vk_ev and mask_vk  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 161  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 171  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  c$$$        if(flag_shower.eqv..true.)then  
 c$$$          call init_level1                
 c$$$          good1=0  
 c$$$          goto 200              !jump to next event  
 c$$$        endif  
 ccc  
 ccc    modified by Elena (08/2006)  
 ccc  
237          if(.not.flag_shower)then          if(.not.flag_shower)then
238             call save_cluster(iv)             call save_cluster(iv)
239               if(debug)print*,'view ',iv,' #clusters ', nclstr_view
240          else          else
241             fshower(iv) = 1             fshower(iv) = 1
242             GOOD1(DSPn) = 11  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            cnnev(iv,ik) = cnn(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 211  C--------------------------------------- Line 265  C---------------------------------------
265        do iv = 1,nviews        do iv = 1,nviews
266           ngood = ngood + good1(iv)           ngood = ngood + good1(iv)
267        enddo        enddo
268        if(ngood.ne.0)print*,'* WARNING * Event ',eventn(1)        if(verbose.and.ngood.ne.0)print*,'* WARNING * Event ',eventn(1)
269       $     ,':LEVEL1 event status: '       $     ,':LEVEL1 event status: '
270       $     ,(good1(i),i=1,nviews)       $     ,(good1(i),i=1,nviews)
271  c------------------------------------------------------------------------  c------------------------------------------------------------------------
# Line 254  c      good1 = 0 Line 308  c      good1 = 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           whichtrack(ic) = 0     !assigned @ level2
313    
314        enddo        enddo
315        do id=1,maxlength         !???        do id=1,maxlength         !???
# Line 275  c        crc1(iv)=0 Line 329  c        crc1(iv)=0
329                
330        return        return
331        end        end
332    
333  *---***---***---***---***---***---***---***---***  *---***---***---***---***---***---***---***---***
334  *  *
335  *  *
# Line 295  c        crc1(iv)=0 Line 350  c        crc1(iv)=0
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
   
       integer multtemp          !temporary multiplicity variable  
356    
357        external nst        external nst
358    
 c------------------------------------------------------------------------  
 c     looks for clusters on each view  
 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)  
   
359        iseed=-999                !cluster seed index initialization        iseed=-999                !cluster seed index initialization
360    
361          inext=-999                !index where to start new cluster search
362    
363          flag_shower = .false.
364        nclstr_view=0        nclstr_view=0
365    
366        do jl=1,nladders_view     !1..3 !loops on ladders        do jl=1,nladders_view              !1..3 !loops on ladders
367           first=1+nstrips_ladder*(jl-1) !1,1025,2049  
368           last=nstrips_ladder*jl !1024,2048,3072           first = 1+nstrips_ladder*(jl-1) !1,1025,2049
369  c     X views have 1018 strips instead of 1024           last  = nstrips_ladder*jl       !1024,2048,3072
370    
371    *        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
# Line 344  c     X views have 1018 strips instead o Line 376  c     X views have 1018 strips instead o
376    
377           do is=first,last       !loop on strips in each ladder           do is=first,last       !loop on strips in each ladder
378    
379              if(is.le.iseed+1) goto 220  c---------------------------------------------
380  *******************************************************  c     new-cluster search starts at index inext
381  *     Elena 08/2006  c---------------------------------------------
382  * QUESTA PARTE NON E` ADEGUATA per cluster con grossi rilasci di carica              if(is.lt.inext) goto 220 ! next strip
 * perche` salva molte volte lo stesso cluster  
 * (salvo il cluster rispetto al primo massimo e basta...)  
 *******************************************************  
 c$$$c-----------------------------------------  
 c$$$c     after a cluster seed as been found,  
 c$$$c     look for next one skipping one strip on the right  
 c$$$c     (i.e. look for double peak cluster)  
 c$$$c-----------------------------------------  
 c$$$            if(is.ne.first) then  
 c$$$               if(value(is).le.value(is-1)) goto 220  
 c$$$            endif  
 c$$$c-----------------------------------------  
 c$$$c     skips cluster seed  
 c$$$c     finding if strips values are descreasing (a strip  
 c$$$c     can be a cluster seed only if previous strip value  
 c$$$c     is lower)  
 c$$$c-----------------------------------------  
 *******************************************************  
 * LA RICERCA PARTE DALL'ULTIMA STRIP SALVATA (***TEMPORANEO****)  
 *******************************************************  
             if(is.le.iseed+rmax+1) goto 220  
 *******************************************************  
383    
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  ****************************************************  c              ------------------------                
393  *     modificato da Elena (08/2006) per salvare  c              search for first maximum
394  *     il cluster intorno al massimo assoluto  c              ------------------------
 ****************************************************  
 c$$$               do while(value(itemp)  
 c$$$     $              /sigma(iv,nvk(itemp),nst(itemp))  
 c$$$     $              .le.value(itemp+1)  
 c$$$     $              /sigma(iv,nvk(itemp+1),nst(itemp+1))) !BIAS: aggiustare il caso uguale!???  
 c$$$                  itemp=itemp+1  
 c$$$                  if(itemp.eq.last) goto 230 !stops if reaches last strip  
 c$$$               enddo            ! of the ladder  
395                 do while(                 do while(
396       $                   value(itemp).le.value(itemp+1)       $                   value(itemp).le.value(itemp+1)
397       $              .and.value(itemp+1).gt.clseedcut(itemp+1))       $              .and.value(itemp+1).gt.clseedcut(itemp+1))
398                    itemp=itemp+1                    itemp = itemp+1
399                    if(itemp.eq.last) goto 230 !stops if reaches last strip                    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 415  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
 c------------------------------------------------------------------------  
                   if((rmax-lmax+1).ge.nclstrp) then  
                      goto 210   !exits inclusion loop:  
                                 ! lmax and rmax maintain last value  
                                 ! NB .ge.!???  
                   endif  
 c------------------------------------------------------------------------  
 c     marks strips exceeding inclusion cut  
454  c------------------------------------------------------------------------  c------------------------------------------------------------------------
455                    if(rstop.eq.0) then !if last strip of the ladder or last                    if( (rmax-lmax+1).ge.nclstrp )goto 210   !exits inclusion loop
456                                  ! over-cut strip has not been reached  
457                       if(value(ir).gt.clinclcut(ir)) then !puts in rmax the                    if(rstop.eq.0) then !if right cluster border has not been reached
458                          rmax=ir ! last right over-cut strip                       if(value(ir).gt.clinclcut(ir)) then
459                            rmax=ir !include a strip on the right
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 border 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    c                  if( (rmax-lmax+1).ge.nclstrp )goto 210   !exits inclusion loop
476    
477                 enddo            !ends strip inclusion loop                 enddo            !ends strip inclusion loop
478                   goto 211
479   210           continue         !jumps here if more than nclstrp have been included   210           continue         !jumps here if more than nclstrp have been included
480                        c               print*,'>>> nclstrp! '
481                 multtemp=rmax-lmax+1 !stores multiplicity in temp   211           continue
482                                  ! variable. NB rmax and lmax can change later in  c-----------------------------------------
483                                  ! order to include enough strips to calculate eta3  c     end of inclusion loop!
484                                  ! and eta4. so mult is not always equal to cllength  c-----------------------------------------
485  c------------------------------------------------------------------------                
 c     NB per essere sicuro di poter calcolare eta3 e eta4 devo includere  
 c     sempre e comunque le 2 strip adiacenti al cluster seed e quella  
 c     adiacente ulteriore dalla parte della piu' alta fra queste due  
 c     (vedi oltre...)!???  
 c------------------------------------------------------------------------  
   
 c     nel caso di estremi del ladder...!???  
   
 c     ho meno di 4 strip nel cluster --> se sono sui bordi o quasi del ladder  
 c     costruisco il cluster ad hoc e poi esco, se non sono sui bordi o quasi  
 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                    
                endif  
486  c------------------------------------------------------------------------  c------------------------------------------------------------------------
487  c     be sure to include in the cluster the cluster seed with its 2 adjacent  c     adjust the cluster in order to have at least a strip around the seed
488  c     strips, and the one adjacent to the greatest between this two strip, as the  c------------------------------------------------------------------------
489  c     fourth one. if the strips have the same value (!) the fourth one is chosen                 if(iseed.eq.lmax.and.lmax.ne.first)then
490  c     as the one having the greatest value between the second neighbors                    lmax = lmax-1
491  c------------------------------------------------------------------------                    if( (rmax-lmax+1).gt.nclstrp )rmax=rmax-1
                if(value(iseed+1).eq.value(iseed-1)) then  
                   if(value(iseed+2).ge.value(iseed-2)) then !??? qui cmq c'e'  
                      diff=(iseed+2)-rmax  
                      if(diff.gt.0) then  
                         rmax=rmax+diff  
                         if((rmax-lmax+1).gt.nclstrp) then  
                            lmax=rmax-nclstrp+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  
                   endif  
492                 endif                 endif
493   250           continue                 if(iseed.eq.rmax.and.rmax.ne.last )then
494                      rmax = rmax+1
495                      if( (rmax-lmax+1).gt.nclstrp )lmax=lmax+1
496                   endif
497    c-------------------------------------------------------------------------------
498    c     adjust the cluster in order to have at least ANOTHER strip around the seed
499    c-------------------------------------------------------------------------------
500    c$$$               if(iseed-1.eq.lmax.and.lmax.ne.first)then
501    c$$$                  lmax = lmax-1
502    c$$$                  if( (rmax-lmax+1).gt.nclstrp )rmax=rmax-1
503    c$$$               endif
504    c$$$               if(iseed+1.eq.rmax.and.rmax.ne.last )then
505    c$$$                  rmax = rmax+1
506    c$$$                  if( (rmax-lmax+1).gt.nclstrp )lmax=lmax+1
507    c$$$               endif
508    c---------------------------------------------------
509    c     now we have 5 stored-strips around the maximum
510    c---------------------------------------------------
511    
512    c------------------------------------------------------------------------
513    c     adjust the cluster in order to store a minimum number of strips
514    c------------------------------------------------------------------------
515                   do while( (rmax-lmax+1).lt.nclstrpmin )
516    
517                      vl = -99999
518                      vr = -99999
519                      if(lmax-1.ge.first) vl = value(lmax-1)
520                      if(rmax+1.le.last ) vr = value(rmax+1)
521                      if(vr.ge.vl) then
522                         rmax = rmax+1
523                      else  
524                         lmax = lmax-1
525                      endif
526                      
527                   enddo
528    
529  c--------------------------------------------------------  c--------------------------------------------------------
530  c     fills cluster variables  c     store cluster info
531  c--------------------------------------------------------  c--------------------------------------------------------
 c$$$               nclstr1=nclstr1+1 !cluster number  
 c$$$ccc               print*,nclstr1,multtemp  
 c$$$               if(nclstr1.gt.nclstrmax) then !too many clusters for the event:  
 c$$$                  if(verbose)print*,'Event ',eventn(1),  
 c$$$     $                 ': more than ',nclstrmax,' clusters'  
 c$$$                  good1=0       ! event  
 c$$$                  nclstr1=0  
 c$$$                  totCLlength=0  
 c$$$                  flag_shower = .true.  
 c$$$                  goto 2000  
 c$$$               endif  
 c$$$               view(nclstr1)   = iv !vista del cluster  
 c$$$               ladder(nclstr1) = nld(iseed,iv) !ladder a cui appartiene il cluster seed  
 c$$$               maxs(nclstr1)   = iseed !strip del cluster seed  
 c$$$               mult(nclstr1)   = multtemp !molteplicita'  
 c$$$                
 c$$$               indstart(nclstr1) = ind !posizione dell'inizio del cluster nell'  
 c$$$c                                      ! array clsignal  
 c$$$               indmax(nclstr1)   = indstart(nclstr1)+(iseed-lmax) !posizione del  
 c$$$c                                      ! cluster seed nell'array clsignal  
 c$$$                
 c$$$               CLlength      = rmax-lmax+1 !numero di strip del cluster  
 c$$$               totCLlength   = totCLlength+CLlength  
 c$$$               dedx(nclstr1) = 0  
 c$$$               do j=lmax,rmax   !stores sequentially cluter strip values in  
 c$$$                  clsignal(ind) = value(j) ! clsignal array  
 c$$$                  ind=ind+1  
 c$$$c                  if(value(j).gt.0)  
 c$$$                  if(value(j).gt.clinclcut(j))  
 c$$$     $                 dedx(nclstr1) = dedx(nclstr1) + value(j) !cluster charge  
 c$$$               enddo  
 ccc  
 ccc            *** Modified by Elena (08/2006) ***  
 ccc  
532                 nclstr_view = nclstr_view + 1 !cluster number                 nclstr_view = nclstr_view + 1 !cluster number
533  c               print*,'view ',iv,' -- search_cluster -- nclstr_view: '  
 c     $              ,nclstr_view  
534                 if(nclstr_view.gt.nclstrmax_view) then !too many clusters for the view:                 if(nclstr_view.gt.nclstrmax_view) then !too many clusters for the view:
535                    if(verbose) print*,'Event ',eventn(1),  c$$$                  if(verbose) print*,'Event ',eventn(1),
536       $                 ': more than ',nclstrmax_view  c$$$     $                 ': more than ',nclstrmax_view
537       $                 ,' clusters on view ',iv  c$$$     $                 ,' clusters on view ',iv
 c                  good1=0       ! event  
 c                  nclstr1=0  
 c                  totCLlength=0  
538                    flag_shower = .true.                    flag_shower = .true.
539                    goto 2000                    goto 2000
540                 endif                 endif
541    
542  c               view(nclstr1)   = iv !vista del cluster                 ladder_view(nclstr_view) = nld(iseed,iv)
543                 ladder_view(nclstr_view) = nld(iseed,iv) !ladder a cui appartiene il cluster seed                 maxs_view(nclstr_view)   = iseed
544                 maxs_view(nclstr_view)   = iseed !strip del cluster seed                 mult_view(nclstr_view)   = rmax-lmax+1
                mult_view(nclstr_view)   = multtemp !molteplicita'  
545                 rmax_view(nclstr_view)   = rmax                 rmax_view(nclstr_view)   = rmax
546                 lmax_view(nclstr_view)   = lmax                 lmax_view(nclstr_view)   = lmax
547    
548    c$$$               if(rmax-lmax+1.gt.25)
549    c$$$     $              print*,'view ',iv
550    c$$$     $              ,' cl ',nclstr_view,' mult ',rmax-lmax+1
551    c------------------------------------------------------------------------
552    c     search for a double peak inside the cluster                                                                                                            
553    c------------------------------------------------------------------------
554                   inext = rmax+1   !<< index where to start new-cluster search
555                  
556                   vmax = 0
557                   vmin = value(iseed)
558                   imax = iseed
559                   imin = iseed
560                   do iss = max(iseed+1,lsat+1),rmax
561                      if( value(iss).lt.vmin )then
562                         if( imax.ne.iseed )goto 221 !found dowble peek
563                         imin = iss
564                         vmin = value(iss)
565                      else
566                         delta = value(iss) - value(imin)
567                         cut = sqrt(clinclcut(iss)**2 + clinclcut(imin)**2)
568                         if(
569         $                    delta.gt.cut .and.
570         $                    value(iss).gt.clseedcut(iss).and.
571         $                    .true.)then
572                            if( value(iss).gt.vmax )then                        
573                               imax = iss
574                               vmax = value(iss)
575                            else
576                               goto 221 !found dowble peek
577                            endif
578                         endif
579                      endif
580                   enddo
581     221           continue
582                  
583                   if(imax.gt.iseed)then
584                      inext = imax    !<< index where to start new-cluster search
585    c$$$                  print*,'--- double peek ---'
586    c$$$                  print*,(value(ii),ii=lmax,rmax)
587    c$$$                  print*,'seed ',iseed,' imin ',imin,' imax ',imax
588                   endif
589  c--------------------------------------------------------  c--------------------------------------------------------
590  c  c
591  c--------------------------------------------------------  c--------------------------------------------------------
# Line 706  c        posizione del cluster seed nell Line 635  c        posizione del cluster seed nell
635                    
636           CLlength      = rmax_view(ic) - lmax_view(ic) + 1 !numero di strip salvate           CLlength      = rmax_view(ic) - lmax_view(ic) + 1 !numero di strip salvate
637           totCLlength   = totCLlength + CLlength           totCLlength   = totCLlength + CLlength
638           dedx(nclstr1) = 0           sgnl(nclstr1) = 0
639           do j=lmax_view(ic),rmax_view(ic)         !stores sequentially cluter strip values in           do j=lmax_view(ic),rmax_view(ic)         !stores sequentially cluter strip values in
640    
641              clsignal(ind) = value(j) ! clsignal array              clsignal(ind) = value(j) ! clsignal array
642    c$$$            print*,ind,clsignal(ind)
643              ivk=nvk(j)              ivk=nvk(j)
644              ist=nst(j)              ist=nst(j)
645    
# Line 722  c            clped(ind)   = pedestal(iv, Line 651  c            clped(ind)   = pedestal(iv,
651              ind=ind+1              ind=ind+1
652  c     if(value(j).gt.0)  c     if(value(j).gt.0)
653              if(value(j).gt.clinclcut(j))              if(value(j).gt.clinclcut(j))
654       $           dedx(nclstr1) = dedx(nclstr1) + value(j) !cluster charge       $           sgnl(nclstr1) = sgnl(nclstr1) + value(j) !cluster charge
655           enddo           enddo
656    
657  c         print*,'view ',iv,' -- save_cluster -- nclstr1: '  c$$$         print*,'view ',iv,' -- save_cluster -- nclstr1: '
658  c     $        ,nclstr1,maxs(nclstr1),mult(nclstr1),dedx(nclstr1)  c$$$     $        ,nclstr1,maxs(nclstr1),mult(nclstr1),sgnl(nclstr1)
659            c$$$         print*,'----------------------'
660    
661        enddo        enddo
662                
663        return        return
# Line 741  c     $        ,nclstr1,maxs(nclstr1),mu Line 671  c     $        ,nclstr1,maxs(nclstr1),mu
671  *---***---***---***---***---***---***---***---***  *---***---***---***---***---***---***---***---***
672    
673    
674        subroutine stripmask  c$$$      subroutine stripmask
675    c$$$
676    c$$$*     this routine set va1 and single-strip masks,
677    c$$$*     on the basis of the VA1 mask saved in the DB
678    c$$$*
679    c$$$*     mask(nviews,nva1_view,nstrips_va1) !strip mask
680    c$$$*     mask_vk(nviews,nva1_view)          !VA1 mask
681    c$$$*
682    c$$$      include 'commontracker.f'
683    c$$$      include 'level1.f'
684    c$$$      include 'common_reduction.f'
685    c$$$      include 'calib.f'
686    c$$$
687    c$$$*     init mask
688    c$$$      do iv=1,nviews
689    c$$$         do ivk=1,nva1_view
690    c$$$            do is=1,nstrips_va1
691    c$$$c               mask(iv,ivk,is) = mask_vk(iv,ivk)
692    c$$$               if( mask_vk(iv,ivk) .ne. -1)then
693    c$$$                  mask(iv,ivk,is) = 1
694    c$$$     $                 * mask_vk(iv,ivk)     !from DB
695    c$$$     $                 * mask_vk_ev(iv,ivk)  !from <SIG>
696    c$$$     $                 * mask_vk_run(iv,ivk) !from CN
697    c$$$               else
698    c$$$                  mask(iv,ivk,is) = -1
699    c$$$     $                 * mask_vk(iv,ivk)     !from DB
700    c$$$     $                 * mask_vk_ev(iv,ivk)  !from CN
701    c$$$               endif
702    c$$$            enddo
703    c$$$         enddo
704    c$$$      enddo
705    c$$$
706    c$$$
707    c$$$      return
708    c$$$      end
709    
710          subroutine stripmask(iv,ivk)
711    
712    *     -----------------------------------------------
713  *     this routine set va1 and single-strip masks,  *     this routine set va1 and single-strip masks,
714  *     on the basis of the VA1 mask saved in the DB  *     on the basis of the VA1 mask saved in the DB
715  *  *
716  *     mask(nviews,nva1_view,nstrips_va1) !strip mask  *     mask(nviews,nva1_view,nstrips_va1) !strip mask
717  *     mask_vk(nviews,nva1_view)          !VA1 mask  *     mask_vk(nviews,nva1_view)          !VA1 mask
718  *  *     -----------------------------------------------
719        include 'commontracker.f'        include 'commontracker.f'
720        include 'level1.f'        include 'level1.f'
721        include 'common_reduction.f'        include 'common_reduction.f'
722        include 'calib.f'        include 'calib.f'
723    
724  *     init mask  *     init mask
725        do iv=1,nviews        do is=1,nstrips_va1
726           do ivk=1,nva1_view  *        --------------------------------------------------------
727              do is=1,nstrips_va1  *        if VA1-mask from DB is 0 or 1, three masks are combined:
728  c               mask(iv,ivk,is) = mask_vk(iv,ivk)  *        - from DB (a-priori mask)
729                 mask(iv,ivk,is) = mask_vk(iv,ivk) * mask_vk_ev(iv,ivk)                *        - run-based (chip declared bad on the basis of <SIG>)
730              enddo  *        - event-based (failure in CN computation)
731           enddo  *        --------------------------------------------------------
732    c         print*,iv,ivk
733    c     $        ,mask_vk(iv,ivk),mask_vk_ev(iv,ivk),mask_vk_run(iv,ivk)
734             if( mask_vk(iv,ivk) .ne. -1)then            
735                mask(iv,ivk,is) = 1
736         $           * mask_vk(iv,ivk)     !from DB
737         $           * mask_vk_ev(iv,ivk)  !from <SIG>
738         $           * mask_vk_run(iv,ivk) !from CN
739    *        -----------------------------------------------------------
740    *        if VA1-mask from DB is -1 only event-based mask is applied
741    *        -----------------------------------------------------------
742             else
743                mask(iv,ivk,is) = -1
744         $           * mask_vk(iv,ivk)     !from DB
745         $           * mask_vk_ev(iv,ivk)  !from CN
746             endif
747        enddo        enddo
748          
749          
750        return        return
751        end        end
   

Legend:
Removed from v.1.9  
changed lines
  Added in v.1.21

  ViewVC Help
Powered by ViewVC 1.1.23