/[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.6 by pam-fi, Thu Sep 28 14:04:40 2006 UTC revision 1.23 by pam-fi, Fri Aug 31 14:56:52 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    c$$$      print*,debug,verbose,warning
31    c$$$      debug=1
32    c$$$      verbose=1
33    c$$$      warning=1
34    
35    *     //////////////////////////
36    *     initialize some parameters
37    *     //////////////////////////
38    
 c      call stripmask   !called later, after CN computation  
39        call init_level1        call init_level1
40    
41    c      debug=.true.
42    
43          if(debug.eq.1)print*,'-- check LEVEL0 status'
44    
45          ievco=-1
46          mismatch=0
47  c      good1 = good0  c      good1 = good0
48  c--------------------------------------------------  c--------------------------------------------------
49  c     check the LEVEL0 event status for missing  c     check the LEVEL0 event status for missing
# Line 41  c--------------------------------------- Line 55  c---------------------------------------
55  c           ------------------------  c           ------------------------
56  c           GOOD  c           GOOD
57  c           ------------------------  c           ------------------------
58              GOOD1(DSPnumber(iv))=0                  !OK              GOOD1(DSPnumber(iv))=0 !OK
59  c           ------------------------  c           ------------------------
60  c           CRC error  c           CRC error
61  c           ------------------------  c           ------------------------
62              if(crc(iv).eq.1) then              if(crc(iv).eq.1) then
63                 GOOD1(DSPnumber(iv)) = 2  c               GOOD1(DSPnumber(iv)) = 2
64                 goto 18 !next view  c               GOOD1(DSPnumber(iv)) = GOOD1(DSPnumber(iv)) + 2**1
65                   GOOD1(DSPnumber(iv)) = ior(GOOD1(DSPnumber(iv)),2**1)
66     102           format(' * WARNING * Event ',i7,' view',i3
67         $          ,' CRC error')
68                   if(debug.eq.1)write(*,102)eventn(1),DSPnumber(iv)
69    c               goto 18 !next view
70              endif              endif
71  c           ------------------------  c           ------------------------
72  c           online-software alarm  c           online-software alarm
# Line 62  c           ------------------------ Line 81  c           ------------------------
81       $           fc(iv).ne.0.or.       $           fc(iv).ne.0.or.
82       $           DATAlength(iv).eq.0.or.       $           DATAlength(iv).eq.0.or.
83       $           .false.)then       $           .false.)then
84                 GOOD1(DSPnumber(iv))=3  c               GOOD1(DSPnumber(iv))=3
85                 goto 18  c               GOOD1(DSPnumber(iv)) = GOOD1(DSPnumber(iv)) + 2**2
86                   GOOD1(DSPnumber(iv)) = ior(GOOD1(DSPnumber(iv)),2**2)
87     103           format(' * WARNING * Event ',i7,' view',i3
88         $          ,' software alarm')
89                   if(debug.eq.1)write(*,103)eventn(1),DSPnumber(iv)
90    c               goto 18
91              endif              endif
92  c           ------------------------  c           ------------------------
93  c           DSP-counter jump  c           DSP-counter jump
94  c           ------------------------  c           ------------------------
95              if(  c     commentato perche` non e` un controllo significativo nel caso in cui
96       $           eventn_old(iv).ne.0.and. !first event in this file  c     la subroutine venga chiamata per riprocessare l'evento
97       $           eventn(iv).ne.1.and.     !first event in run  c     sostituito con un check dei contatori dei vari dsp
98       $           good_old(DSPnumber(iv)).ne.0.and. !previous event corrupted  c$$$            if(
99       $           .true.)then  c$$$     $           eventn_old(iv).ne.0.and. !first event in this file
100    c$$$     $           eventn(iv).ne.1.and.     !first event in run
101                 if(eventn(iv).ne.(eventn_old(iv)+1))then  c$$$     $           good_old(DSPnumber(iv)).ne.0.and. !previous event corrupted
102                    GOOD1(DSPnumber(iv))=4  c$$$     $           .true.)then
103                    goto 18  c$$$
104    c$$$               if(eventn(iv).ne.(eventn_old(iv)+1))then
105    c$$$c                  GOOD1(DSPnumber(iv))=4
106    c$$$c                  GOOD1(DSPnumber(iv)) = GOOD1(DSPnumber(iv)) + 2**3
107    c$$$                  GOOD1(DSPnumber(iv)) = ior(GOOD1(DSPnumber(iv)),2**3)
108    c$$$ 104              format(' * WARNING * Event ',i7,' view',i3
109    c$$$     $          ,' counter jump ',i10,i10)
110    c$$$                  if(debug.eq.1)write(*,104)eventn(1),DSPnumber(iv)
111    c$$$     $                 ,eventn_old(iv),eventn(iv))
112    c$$$                  goto 18
113    c$$$               endif
114    c$$$
115    c$$$            endif
116    c           ------------------------
117    c 18         continue
118    c           ------------------------
119    c           DSP-counter
120    c           ------------------------
121                if( DSPnumber(iv).ne.0.and.GOOD1(DSPnumber(iv)).ne.1)then
122                   if(iv.ne.1.and.ievco.ne.-1)then
123                      if( eventn(iv).ne.ievco )then
124                         mismatch=1
125                      endif
126                 endif                 endif
127                   ievco = eventn(iv)
128              endif              endif
 c           ------------------------  
  18         continue  
129           endif           endif
130        enddo        enddo
131    
132    c      print*,'*** ',(eventn(iv),iv=1,12)
133          
134          if(mismatch.eq.1.and.debug.eq.1)
135         $     print*,' * WARNING * DSP counter mismatch: '
136         $     ,(eventn(iv),iv=1,12)
137    
138        ngood = 0        ngood = 0
139        do iv = 1,nviews        do iv = 1,nviews
140            
141             if(mismatch.eq.1.and.GOOD1(iv).ne.1)
142         $        GOOD1(iv)=ior(GOOD1(iv),2**3)
143    
144           eventn_old(iv) = eventn(iv)           eventn_old(iv) = eventn(iv)
145           good_old(iv)   = good1(iv)           good_old(iv)   = good1(iv)
146           ngood = ngood + good1(iv)           ngood = ngood + good1(iv)
147    
148        enddo        enddo
149  c      if(ngood.ne.0)print*,'* WARNING * LEVEL0 event status: '  c$$$      if(verbose.and.ngood.ne.0)print*,'* WARNING * Event ',eventn(1)
150  c     $     ,(good1(i),i=1,nviews)  c$$$     $     ,':LEVEL0 event status: '
151    c$$$     $     ,(good1(i),i=1,nviews)
152  c--------------------------------------------------  c--------------------------------------------------
153  c     read the variable DATATRACKER from LEVEL0  c     read the variable DATATRACKER from LEVEL0
154  c     and fill the variable ADC (invertin view 11)  c     and fill the variable ADC (invertin view 11)
155  c--------------------------------------------------  c--------------------------------------------------
156          
157          if(debug.eq.1)print*,'-- fill ADC vectors'
158    
159        call filladc(iflag)        call filladc(iflag)
160        if(iflag.ne.0)then        if(iflag.ne.0)then
 c        good1=0!<<<<<<<<<<<<<<<  
 c       if(DEBUG)print*,'event ',eventn(1),' >>>>>  decode ERROR'  
161           ierror = 220           ierror = 220
 c        goto 200  
 c         print*,'filladc error'  
162        endif        endif
163    
164  c--------------------------------------------------  c--------------------------------------------------
165  c     computes common noise for each VA1  c     computes common noise for each VA1
166  c     (excluding strips affected by signal,  c     (excluding strips with signal,
167  c     tagged with the flag CLSTR)  c     tagged with the flag CLSTR)
168  c--------------------------------------------------  c--------------------------------------------------
169          if(debug.eq.1)print*,'-- compute CN'
170    
171        do iv=1,nviews        do iv=1,nviews
172          do ik=1,nva1_view           ima=0
173            cn(iv,ik)  = 0           do ik=1,nva1_view
174            cnn(iv,ik) = -1              cn(iv,ik)    = 0
175            mask_vk_ev(iv,ik)=1              cnrms(iv,ik) = 0
176            iflag=0              cnn(iv,ik)   = -1
177            if(mask_vk(iv,ik).eq.1)call cncomp(iv,ik,iflag)              iflag = 0
178  c          if(iflag.ne.0)good1=0              mask_vk_ev(iv,ik) = 1
179            if(iflag.ne.0)then              call stripmask(iv,ik)      !compute mask(i,j,k), combining VA1-masks
180               mask_vk_ev(iv,ik)=0  *           --------------------------------------
181               ierror = 220  *           if chip is not masked ---> evaluate CN
182            endif  *           --------------------------------------
183          enddo              if( mask(iv,ik,1).eq.1 ) then !!!NBNB mask per la striscia 1 !!!!!!!!
184                   call cncomp(iv,ik,iflag)
185                   if(iflag.ne.0)then
186                      ima=ima+1
187                      mask_vk_ev(iv,ik)=0
188                      ierror = 220
189                   endif
190                   call stripmask(iv,ik) !compute mask(i,j,k), combining VA1-masks
191                endif
192             enddo
193     100     format(' * WARNING * Event ',i7,' view',i3,': VK MASK ',24i1)
194             if(ima.ne.0.and.verbose.eq.1)write(*,100)eventn(1),iv
195         $        ,(mask_vk_ev(iv,ik),ik=1,nva1_view)
196    c         if(ima.ne.0)write(*,100)eventn(1),iv
197    c     $        ,(mask_vk_ev(iv,ik),ik=1,nva1_view)  
198        enddo        enddo
 c      if(good1.eq.0)then  
 c         ierror = 220  
 c      endif  
199    
200        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
201    
202  c---------------------------------------------  c---------------------------------------------
203  c     loops on views, VA1 and strips,  c     loops on views, VA1 and strips,
204  c     and computes strips signals using  c     and computes strips signals using
205  c     badstrip, pedestals, and  c     badstrip, pedestals, and
206  c     sigma informations from histograms  c     sigma informations from histograms
207  c---------------------------------------------  c---------------------------------------------
       flag_shower = .false.  
208        ind=1                     !clsignal array index        ind=1                     !clsignal array index
209    
210          if(debug.eq.1)print*,'-- search clusters'
211        do iv=1,nviews            !loop on views        do iv=1,nviews            !loop on views
212          do is=1,nstrips_view    !loop on strips (1)          do is=1,nstrips_view    !loop on strips (1)
213            if(mod(iv,2).eq.1) then            if(mod(iv,2).eq.1) then
214  C===  > Y view  C===  > Y view
215    c             print*,iv,nvk(is),nst(is),adc(iv,nvk(is),nst(is))
216    c     $            ,cn(iv,nvk(is))
217    c     $            ,pedestal(iv,nvk(is),nst(is))
218              value(is)= -(DBLE(adc(iv,nvk(is),nst(is)))              value(is)= -(DBLE(adc(iv,nvk(is),nst(is)))
219       $           -cn(iv,nvk(is))-pedestal(iv,nvk(is),nst(is)))       $           -cn(iv,nvk(is))-pedestal(iv,nvk(is),nst(is)))
220       $           *mask(iv,nvk(is),nst(is))       $           *mask(iv,nvk(is),nst(is))
# Line 150  C===  > Y view Line 222  C===  > Y view
222       $           *mask(iv,nvk(is),nst(is))       $           *mask(iv,nvk(is),nst(is))
223              clinclcut(is)=incuty*sigma(iv,nvk(is),nst(is))              clinclcut(is)=incuty*sigma(iv,nvk(is),nst(is))
224       $           *mask(iv,nvk(is),nst(is))       $           *mask(iv,nvk(is),nst(is))
225  ccc            print*,"value(",is,")(reduction)= ",value(is)              sat(is)=0
226                if( adc(iv,nvk(is),nst(is)).lt.adc_saty )sat(is)=1
227            else                        else            
228  C===  > X view  C===  > X view
229              value(is)= (DBLE(adc(iv,nvk(is),nst(is)))              value(is)= (DBLE(adc(iv,nvk(is),nst(is)))
# Line 160  C===  > X view Line 233  C===  > X view
233       $           *mask(iv,nvk(is),nst(is))       $           *mask(iv,nvk(is),nst(is))
234              clinclcut(is)=incutx*sigma(iv,nvk(is),nst(is))              clinclcut(is)=incutx*sigma(iv,nvk(is),nst(is))
235       $           *mask(iv,nvk(is),nst(is))       $           *mask(iv,nvk(is),nst(is))
236                sat(is)=0
237                if( adc(iv,nvk(is),nst(is)).gt.adc_satx )sat(is)=1
238            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)  
239          enddo                   !end loop on strips (1)          enddo                   !end loop on strips (1)
240          call search_cluster(iv)          call search_cluster(iv)
241  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  
242          if(.not.flag_shower)then          if(.not.flag_shower)then
243             call save_cluster(iv)             call save_cluster(iv)
244               if(debug.eq.1)print*,'view ',iv,' #clusters ', nclstr_view
245          else          else
246             fshower(iv) = 1             fshower(iv) = 1
247             GOOD1(DSPn) = 11  c           GOOD1(DSPnumber(iv)) = 11 !AHAHAHAHA ORRORE!!
248    c           GOOD1(iv) = 11
249    c           GOOD1(iv) = GOOD1(iv) + 2**5
250               GOOD1(iv) = ior(GOOD1(iv),2**5)
251     101       format(' * WARNING * Event ',i7,' view',i3
252         $          ,' #clusters > ',i5,' --> MASKED')
253               if(verbose.eq.1)write(*,101)eventn(1),iv,nclstrmax_view
254          endif          endif
255        enddo                     ! end loop on views        enddo                     ! end loop on views
256        do iv=1,nviews        do iv=1,nviews
257          do ik=1,nva1_view          do ik=1,nva1_view
258            cnev(iv,ik)  = cn(iv,ik) !assigns computed CN to ntuple variables            cnev(iv,ik)    = cn(iv,ik) !assigns computed CN to ntuple variables
259            cnnev(iv,ik) = cnn(iv,ik) !assigns computed CN to ntuple variables            cnrmsev(iv,ik) = cnrms(iv,ik) !assigns computed CN to ntuple variables
260  ccc          print*,"cnev(",iv,",",ik,")(reduction)= ",cnev(iv,ik)            cnnev(iv,ik)   = cnn(iv,ik) !assigns computed CN to ntuple variables
261          enddo          enddo
262        enddo        enddo
263  C---------------------------------------------  C---------------------------------------------
# Line 200  C--------------------------------------- Line 270  C---------------------------------------
270        do iv = 1,nviews        do iv = 1,nviews
271           ngood = ngood + good1(iv)           ngood = ngood + good1(iv)
272        enddo        enddo
273        if(ngood.ne.0)print*,'* WARNING * LEVEL1 event status: '        if(verbose.eq.1.and.ngood.ne.0)
274         $     print*,'* WARNING * Event ',eventn(1)
275         $     ,':LEVEL1 event status: '
276       $     ,(good1(i),i=1,nviews)       $     ,(good1(i),i=1,nviews)
277  c------------------------------------------------------------------------  c------------------------------------------------------------------------
278  c  c
# Line 242  c      good1 = 0 Line 314  c      good1 = 0
314           indmax(ic) = 0           indmax(ic) = 0
315           maxs(ic) = 0           maxs(ic) = 0
316           mult(ic) = 0                     mult(ic) = 0          
317           dedx(ic) = 0           sgnl(ic) = 0
318           whichtrack(ic) = 0           whichtrack(ic) = 0     !assigned @ level2
319    
320        enddo        enddo
321        do id=1,maxlength         !???        do id=1,maxlength         !???
# Line 263  c        crc1(iv)=0 Line 335  c        crc1(iv)=0
335                
336        return        return
337        end        end
338    
339  *---***---***---***---***---***---***---***---***  *---***---***---***---***---***---***---***---***
340  *  *
341  *  *
# Line 283  c        crc1(iv)=0 Line 356  c        crc1(iv)=0
356    
357  c     local variables  c     local variables
358        integer rmax,lmax         !estremi del cluster        integer rmax,lmax         !estremi del cluster
359        integer rstop,lstop       !per decidere quali strip includere nel cluster        integer rstop,lstop      
360                                  ! oltre il seed        integer first,last
361        integer first,last,diff   !per includere le strip giuste... !???        integer fsat,lsat
   
       integer multtemp          !temporary multiplicity variable  
362    
363        external nst        external nst
364    
 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)  
   
365        iseed=-999                !cluster seed index initialization        iseed=-999                !cluster seed index initialization
366    
367          inext=-999                !index where to start new cluster search
368    
369          flag_shower = .false.
370        nclstr_view=0        nclstr_view=0
371    
372        do jl=1,nladders_view     !1..3 !loops on ladders        do jl=1,nladders_view              !1..3 !loops on ladders
373           first=1+nstrips_ladder*(jl-1) !1,1025,2049  
374           last=nstrips_ladder*jl !1024,2048,3072           first = 1+nstrips_ladder*(jl-1) !1,1025,2049
375  c     X views have 1018 strips instead of 1024           last  = nstrips_ladder*jl       !1024,2048,3072
376    
377    *        X views have 1018 strips instead of 1024
378           if(mod(iv,2).eq.0) then           if(mod(iv,2).eq.0) then
379              first=first+3              first=first+3
380              last=last-3              last=last-3
# Line 332  c     X views have 1018 strips instead o Line 382  c     X views have 1018 strips instead o
382    
383           do is=first,last       !loop on strips in each ladder           do is=first,last       !loop on strips in each ladder
384    
385              if(is.le.iseed+1) goto 220  c---------------------------------------------
386  *******************************************************  c     new-cluster search starts at index inext
387  *     Elena 08/2006  c---------------------------------------------
388  * 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  
 *******************************************************  
389    
390              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)  
391  c-----------------------------------------  c-----------------------------------------
392  c     possible SEED...  c     possible SEED...
393  c-----------------------------------------  c-----------------------------------------
394                 itemp=is                 itemp = is
395                   fsat = 0         ! first saturated strip
396                   lsat = 0         ! last saturated strip
397                 if(itemp.eq.last) goto 230 !estremo...                 if(itemp.eq.last) goto 230 !estremo...
398  ****************************************************  c              ------------------------                
399  *     modificato da Elena (08/2006) per salvare  c              search for first maximum
400  *     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  
401                 do while(                 do while(
402       $                   value(itemp).le.value(itemp+1)       $                   value(itemp).le.value(itemp+1)
403       $              .and.value(itemp+1).gt.clseedcut(itemp+1))       $              .and.value(itemp+1).gt.clseedcut(itemp+1))
404                    itemp=itemp+1                    itemp = itemp+1
405                    if(itemp.eq.last) goto 230 !stops if reaches last strip                    if(itemp.eq.last)   goto 230 !stops if reaches last strip
406                      if(sat(itemp).eq.1) goto 230 !stop if reaches a saturated strip
407                 enddo            ! of the ladder                 enddo            ! of the ladder
408   230           continue   230           continue
409  c-----------------------------------------  c              -----------------------------
410    c              check if strips are saturated
411    c              -----------------------------    
412                   if( sat(itemp).eq.1 )then
413                      fsat = itemp
414                      lsat = itemp
415                      if(itemp.eq.last) goto 231 !estremo...
416                      do while( sat(itemp+1).eq.1 )
417                         itemp = itemp+1
418                         lsat = itemp
419                         if(itemp.eq.last)   goto 231 !stops if reaches last strip
420                      enddo                  
421                   endif
422     231           continue
423    c---------------------------------------------------------------------------
424  c     fownd SEED!!!  c     fownd SEED!!!
425  c-----------------------------------------  c     (if there are saturated strips, the cluster is centered in the middle)
426                 iseed=itemp      c---------------------------------------------------------------------------
427  c----------------------------------------------------------                 if(fsat.eq.0.and.lsat.eq.0)then
428                      iseed = itemp ! <<< SEED
429                   else
430                      iseed = int((lsat+fsat)/2) ! <<< SEED
431    c$$$                  print*,'saturated strips ',fsat,lsat,iseed
432    c$$$                  print*,'--> ',(value(ii),ii=fsat,lsat)
433                   endif    
434    c---------------------------------------------------------------
435  c     after finding a cluster seed, checks also adjacent strips,  c     after finding a cluster seed, checks also adjacent strips,
436  C     and marks the ones exceeding clinclcut  C     and tags the ones exceeding clinclcut
437  c----------------------------------------------------------  c---------------------------------------------------------------
438                 ir=iseed         !indici destro                 ir=iseed         !indici destro
439                 il=iseed         ! e sinistro                 il=iseed         ! e sinistro
440                                
# Line 403  c--------------------------------------- Line 445  c---------------------------------------
445                 lstop=0          ! inclusion loop                 lstop=0          ! inclusion loop
446    
447                 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
448                    ir=ir+1       !position index for strips on right side of  
449                                  ! cluster seed  
450                    il=il-1       !and for left side                    ir=ir+1       !index for right side
451                      il=il-1       !index for left side
452  c------------------------------------------------------------------------  c------------------------------------------------------------------------
453  c     checks for last or first strip of the ladder  c     checks for last or first strip of the ladder
454  c------------------------------------------------------------------------  c------------------------------------------------------------------------
455                    if(ir.gt.last) then !when index goes beyond last strip                    if( ir.gt.last  ) rstop = 1                      
456                       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  
457                                        
458  c------------------------------------------------------------------------  c------------------------------------------------------------------------
459  c     check for clusters including more than nclstrp strips  c     add strips exceeding inclusion cut
460  c------------------------------------------------------------------------  c------------------------------------------------------------------------
461                    if((rmax-lmax+1).ge.nclstrp) then                    if( (rmax-lmax+1).ge.nclstrp )goto 210   !exits inclusion loop
462                       goto 210   !exits inclusion loop:  
463                                  ! lmax and rmax maintain last value                    if(rstop.eq.0) then !if right cluster border has not been reached
464                                  ! NB .ge.!???                       if(value(ir).gt.clinclcut(ir)) then
465                    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  
466                       else                       else
467                          rstop=1 !otherwise cluster ends on right and rstop                          rstop=1 !cluster right end
468                       endif      ! flag=1 signals it                       endif    
469                    endif                    endif
470                    if(lstop.eq.0) then  
471                      if( (rmax-lmax+1).ge.nclstrp )goto 210   !exits inclusion loop
472    
473                      if(lstop.eq.0) then !if left cluster border has not been reached
474                       if(value(il).gt.clinclcut(il)) then                       if(value(il).gt.clinclcut(il)) then
475                          lmax=il                          lmax=il !include a strip on the left
476                       else                       else
477                          lstop=1                          lstop=1 !cluster left end
478                       endif                       endif
479                    endif                    endif
480    
481    c                  if( (rmax-lmax+1).ge.nclstrp )goto 210   !exits inclusion loop
482    
483                 enddo            !ends strip inclusion loop                 enddo            !ends strip inclusion loop
484                   goto 211
485   210           continue         !jumps here if more than nclstrp have been included   210           continue         !jumps here if more than nclstrp have been included
486                        c               print*,'>>> nclstrp! '
487                 multtemp=rmax-lmax+1 !stores multiplicity in temp   211           continue
488                                  ! variable. NB rmax and lmax can change later in  c-----------------------------------------
489                                  ! order to include enough strips to calculate eta3  c     end of inclusion loop!
490                                  ! and eta4. so mult is not always equal to cllength  c-----------------------------------------
491  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  
492  c------------------------------------------------------------------------  c------------------------------------------------------------------------
493  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
494  c     strips, and the one adjacent to the greatest between this two strip, as the  c------------------------------------------------------------------------
495  c     fourth one. if the strips have the same value (!) the fourth one is chosen                 if(iseed.eq.lmax.and.lmax.ne.first)then
496  c     as the one having the greatest value between the second neighbors                    lmax = lmax-1
497  c------------------------------------------------------------------------                    if( (rmax-lmax+1).gt.nclstrp )rmax=rmax-1
498                 if(value(iseed+1).eq.value(iseed-1)) then                 endif
499                    if(value(iseed+2).ge.value(iseed-2)) then !??? qui cmq c'e'                 if(iseed.eq.rmax.and.rmax.ne.last )then
500                       diff=(iseed+2)-rmax                    rmax = rmax+1
501                       if(diff.gt.0) then                    if( (rmax-lmax+1).gt.nclstrp )lmax=lmax+1
                         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  
502                 endif                 endif
503   250           continue  c-------------------------------------------------------------------------------
504    c     adjust the cluster in order to have at least ANOTHER strip around the seed
505    c-------------------------------------------------------------------------------
506    c$$$               if(iseed-1.eq.lmax.and.lmax.ne.first)then
507    c$$$                  lmax = lmax-1
508    c$$$                  if( (rmax-lmax+1).gt.nclstrp )rmax=rmax-1
509    c$$$               endif
510    c$$$               if(iseed+1.eq.rmax.and.rmax.ne.last )then
511    c$$$                  rmax = rmax+1
512    c$$$                  if( (rmax-lmax+1).gt.nclstrp )lmax=lmax+1
513    c$$$               endif
514    c---------------------------------------------------
515    c     now we have 5 stored-strips around the maximum
516    c---------------------------------------------------
517    
518    c------------------------------------------------------------------------
519    c     adjust the cluster in order to store a minimum number of strips
520    c------------------------------------------------------------------------
521                   do while( (rmax-lmax+1).lt.nclstrpmin )
522    
523                      vl = -99999
524                      vr = -99999
525                      if(lmax-1.ge.first) vl = value(lmax-1)
526                      if(rmax+1.le.last ) vr = value(rmax+1)
527                      if(vr.ge.vl) then
528                         rmax = rmax+1
529                      else  
530                         lmax = lmax-1
531                      endif
532                      
533                   enddo
534    
535  c--------------------------------------------------------  c--------------------------------------------------------
536  c     fills cluster variables  c     store cluster info
537  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  
538                 nclstr_view = nclstr_view + 1 !cluster number                 nclstr_view = nclstr_view + 1 !cluster number
539  c               print*,'view ',iv,' -- search_cluster -- nclstr_view: '  
 c     $              ,nclstr_view  
540                 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:
541                    if(verbose) print*,'Event ',eventn(1),  c$$$                  if(verbose) print*,'Event ',eventn(1),
542       $                 ': more than ',nclstrmax_view  c$$$     $                 ': more than ',nclstrmax_view
543       $                 ,' clusters on view ',iv  c$$$     $                 ,' clusters on view ',iv
 c                  good1=0       ! event  
 c                  nclstr1=0  
 c                  totCLlength=0  
544                    flag_shower = .true.                    flag_shower = .true.
545                    goto 2000                    goto 2000
546                 endif                 endif
547    
548  c               view(nclstr1)   = iv !vista del cluster                 ladder_view(nclstr_view) = nld(iseed,iv)
549                 ladder_view(nclstr_view) = nld(iseed,iv) !ladder a cui appartiene il cluster seed                 maxs_view(nclstr_view)   = iseed
                maxs_view(nclstr_view)   = iseed !strip del cluster seed  
                mult_view(nclstr_view)   = multtemp !molteplicita'  
550                 rmax_view(nclstr_view)   = rmax                 rmax_view(nclstr_view)   = rmax
551                 lmax_view(nclstr_view)   = lmax                 lmax_view(nclstr_view)   = lmax
552    c               mult_view(nclstr_view)   = rmax-lmax+1
553                   mult_view(nclstr_view)   = 0
554                   do ii=lmax,rmax
555                      if(value(ii).gt.clinclcut(ii))  
556         $                 mult_view(nclstr_view) = mult_view(nclstr_view)+1
557                   enddo
558    
559    
560    c$$$               if(rmax-lmax+1.gt.25)
561    c$$$     $              print*,'view ',iv
562    c$$$     $              ,' cl ',nclstr_view,' mult ',rmax-lmax+1
563    c------------------------------------------------------------------------
564    c     search for a double peak inside the cluster                                                                                                            
565    c------------------------------------------------------------------------
566                   inext = rmax+1   !<< index where to start new-cluster search
567                  
568                   vmax = 0
569                   vmin = value(iseed)
570                   imax = iseed
571                   imin = iseed
572                   do iss = max(iseed+1,lsat+1),rmax
573                      if( value(iss).lt.vmin )then
574                         if( imax.ne.iseed )goto 221 !found dowble peek
575                         imin = iss
576                         vmin = value(iss)
577                      else
578                         delta = value(iss) - value(imin)
579                         cut = sqrt(clinclcut(iss)**2 + clinclcut(imin)**2)
580                         if(
581         $                    delta.gt.cut .and.
582         $                    value(iss).gt.clseedcut(iss).and.
583         $                    .true.)then
584                            if( value(iss).gt.vmax )then                        
585                               imax = iss
586                               vmax = value(iss)
587                            else
588                               goto 221 !found dowble peek
589                            endif
590                         endif
591                      endif
592                   enddo
593     221           continue
594                  
595                   if(imax.gt.iseed)then
596                      inext = imax    !<< index where to start new-cluster search
597    c$$$                  print*,'--- double peek ---'
598    c$$$                  print*,(value(ii),ii=lmax,rmax)
599    c$$$                  print*,'seed ',iseed,' imin ',imin,' imax ',imax
600                   endif
601  c--------------------------------------------------------  c--------------------------------------------------------
602  c  c
603  c--------------------------------------------------------  c--------------------------------------------------------
# Line 694  c        posizione del cluster seed nell Line 647  c        posizione del cluster seed nell
647                    
648           CLlength      = rmax_view(ic) - lmax_view(ic) + 1 !numero di strip salvate           CLlength      = rmax_view(ic) - lmax_view(ic) + 1 !numero di strip salvate
649           totCLlength   = totCLlength + CLlength           totCLlength   = totCLlength + CLlength
650           dedx(nclstr1) = 0           sgnl(nclstr1) = 0
651           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
652    
653              clsignal(ind) = value(j) ! clsignal array              clsignal(ind) = value(j) ! clsignal array
654    c$$$            print*,ind,clsignal(ind)
655              ivk=nvk(j)              ivk=nvk(j)
656              ist=nst(j)              ist=nst(j)
657    
# Line 710  c            clped(ind)   = pedestal(iv, Line 663  c            clped(ind)   = pedestal(iv,
663              ind=ind+1              ind=ind+1
664  c     if(value(j).gt.0)  c     if(value(j).gt.0)
665              if(value(j).gt.clinclcut(j))              if(value(j).gt.clinclcut(j))
666       $           dedx(nclstr1) = dedx(nclstr1) + value(j) !cluster charge       $           sgnl(nclstr1) = sgnl(nclstr1) + value(j) !cluster charge
667           enddo           enddo
668    
669  c         print*,'view ',iv,' -- save_cluster -- nclstr1: '  c$$$         print*,'view ',iv,' -- save_cluster -- nclstr1: '
670  c     $        ,nclstr1,maxs(nclstr1),mult(nclstr1),dedx(nclstr1)  c$$$     $        ,nclstr1,maxs(nclstr1),mult(nclstr1),sgnl(nclstr1)
671            c$$$         print*,'----------------------'
672    
673        enddo        enddo
674                
675        return        return
# Line 729  c     $        ,nclstr1,maxs(nclstr1),mu Line 683  c     $        ,nclstr1,maxs(nclstr1),mu
683  *---***---***---***---***---***---***---***---***  *---***---***---***---***---***---***---***---***
684    
685    
686        subroutine stripmask  c$$$      subroutine stripmask
687    c$$$
688    c$$$*     this routine set va1 and single-strip masks,
689    c$$$*     on the basis of the VA1 mask saved in the DB
690    c$$$*
691    c$$$*     mask(nviews,nva1_view,nstrips_va1) !strip mask
692    c$$$*     mask_vk(nviews,nva1_view)          !VA1 mask
693    c$$$*
694    c$$$      include 'commontracker.f'
695    c$$$      include 'level1.f'
696    c$$$      include 'common_reduction.f'
697    c$$$      include 'calib.f'
698    c$$$
699    c$$$*     init mask
700    c$$$      do iv=1,nviews
701    c$$$         do ivk=1,nva1_view
702    c$$$            do is=1,nstrips_va1
703    c$$$c               mask(iv,ivk,is) = mask_vk(iv,ivk)
704    c$$$               if( mask_vk(iv,ivk) .ne. -1)then
705    c$$$                  mask(iv,ivk,is) = 1
706    c$$$     $                 * mask_vk(iv,ivk)     !from DB
707    c$$$     $                 * mask_vk_ev(iv,ivk)  !from <SIG>
708    c$$$     $                 * mask_vk_run(iv,ivk) !from CN
709    c$$$               else
710    c$$$                  mask(iv,ivk,is) = -1
711    c$$$     $                 * mask_vk(iv,ivk)     !from DB
712    c$$$     $                 * mask_vk_ev(iv,ivk)  !from CN
713    c$$$               endif
714    c$$$            enddo
715    c$$$         enddo
716    c$$$      enddo
717    c$$$
718    c$$$
719    c$$$      return
720    c$$$      end
721    
722          subroutine stripmask(iv,ivk)
723    
724    *     -----------------------------------------------
725  *     this routine set va1 and single-strip masks,  *     this routine set va1 and single-strip masks,
726  *     on the basis of the VA1 mask saved in the DB  *     on the basis of the VA1 mask saved in the DB
727  *  *
728  *     mask(nviews,nva1_view,nstrips_va1) !strip mask  *     mask(nviews,nva1_view,nstrips_va1) !strip mask
729  *     mask_vk(nviews,nva1_view)          !VA1 mask  *     mask_vk(nviews,nva1_view)          !VA1 mask
730  *  *     -----------------------------------------------
731        include 'commontracker.f'        include 'commontracker.f'
732        include 'level1.f'        include 'level1.f'
733        include 'common_reduction.f'        include 'common_reduction.f'
734        include 'calib.f'        include 'calib.f'
735    
736  *     init mask  *     init mask
737        do iv=1,nviews        do is=1,nstrips_va1
738           do ivk=1,nva1_view  *        --------------------------------------------------------
739              do is=1,nstrips_va1  *        if VA1-mask from DB is 0 or 1, three masks are combined:
740  c               mask(iv,ivk,is) = mask_vk(iv,ivk)  *        - from DB (a-priori mask)
741                 mask(iv,ivk,is) = mask_vk(iv,ivk) * mask_vk_ev(iv,ivk)                *        - run-based (chip declared bad on the basis of <SIG>)
742              enddo  *        - event-based (failure in CN computation)
743           enddo  *        --------------------------------------------------------
744    c         print*,iv,ivk
745    c     $        ,mask_vk(iv,ivk),mask_vk_ev(iv,ivk),mask_vk_run(iv,ivk)
746             if( mask_vk(iv,ivk) .ne. -1)then            
747                mask(iv,ivk,is) = 1
748         $           * mask_vk(iv,ivk)     !from DB
749         $           * mask_vk_ev(iv,ivk)  !from <SIG>
750         $           * mask_vk_run(iv,ivk) !from CN
751    *        -----------------------------------------------------------
752    *        if VA1-mask from DB is -1 only event-based mask is applied
753    *        -----------------------------------------------------------
754             else
755                mask(iv,ivk,is) = -1
756         $           * mask_vk(iv,ivk)     !from DB
757         $           * mask_vk_ev(iv,ivk)  !from CN
758             endif
759        enddo        enddo
760          
761          
762        return        return
763        end        end
   

Legend:
Removed from v.1.6  
changed lines
  Added in v.1.23

  ViewVC Help
Powered by ViewVC 1.1.23