/[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.17 by pam-fi, Thu Mar 15 12:17:10 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    
 c      call stripmask   !called later, after CN computation  
30        call init_level1        call init_level1
31    
32          if(debug)print*,'-- check LEVEL0 status'
33    
34  c      good1 = good0  c      good1 = good0
35  c--------------------------------------------------  c--------------------------------------------------
36  c     check the LEVEL0 event status for missing  c     check the LEVEL0 event status for missing
# Line 41  c--------------------------------------- Line 42  c---------------------------------------
42  c           ------------------------  c           ------------------------
43  c           GOOD  c           GOOD
44  c           ------------------------  c           ------------------------
45              GOOD1(DSPnumber(iv))=0                  !OK              GOOD1(DSPnumber(iv))=0 !OK
46  c           ------------------------  c           ------------------------
47  c           CRC error  c           CRC error
48  c           ------------------------  c           ------------------------
# Line 91  c           ------------------------ Line 92  c           ------------------------
92           good_old(iv)   = good1(iv)           good_old(iv)   = good1(iv)
93           ngood = ngood + good1(iv)           ngood = ngood + good1(iv)
94        enddo        enddo
95  c      if(ngood.ne.0)print*,'* WARNING * LEVEL0 event status: '        if(debug.and.ngood.ne.0)print*,'* WARNING * LEVEL0 event status: '
96  c     $     ,(good1(i),i=1,nviews)       $     ,(good1(i),i=1,nviews)
97  c--------------------------------------------------  c--------------------------------------------------
98  c     read the variable DATATRACKER from LEVEL0  c     read the variable DATATRACKER from LEVEL0
99  c     and fill the variable ADC (invertin view 11)  c     and fill the variable ADC (invertin view 11)
100  c--------------------------------------------------  c--------------------------------------------------
101          
102          if(debug)print*,'-- fill ADC vectors'
103    
104        call filladc(iflag)        call filladc(iflag)
105        if(iflag.ne.0)then        if(iflag.ne.0)then
 c        good1=0!<<<<<<<<<<<<<<<  
 c       if(DEBUG)print*,'event ',eventn(1),' >>>>>  decode ERROR'  
106           ierror = 220           ierror = 220
 c        goto 200  
 c         print*,'filladc error'  
107        endif        endif
108    
109  c--------------------------------------------------  c--------------------------------------------------
110  c     computes common noise for each VA1  c     computes common noise for each VA1
111  c     (excluding strips affected by signal,  c     (excluding strips with signal,
112  c     tagged with the flag CLSTR)  c     tagged with the flag CLSTR)
113  c--------------------------------------------------  c--------------------------------------------------
114          if(debug)print*,'-- compute CN'
115    
116        do iv=1,nviews        do iv=1,nviews
117          do ik=1,nva1_view           ima=0
118            cn(iv,ik)  = 0           do ik=1,nva1_view
119            cnn(iv,ik) = -1              cn(iv,ik)  = 0
120            mask_vk_ev(iv,ik)=1              cnrms(iv,ik)  = 0
121            iflag=0              cnn(iv,ik) = -1
122            if(mask_vk(iv,ik).eq.1)call cncomp(iv,ik,iflag)              iflag=0
123  c          if(iflag.ne.0)good1=0              mask_vk_ev(iv,ik)=1
124            if(iflag.ne.0)then              call stripmask(iv,ik)      !compute mask(i,j,k), combining VA1-masks
125               mask_vk_ev(iv,ik)=0  c     NBNBNBNBNB mask per la striscia 1 !!!!!!!!
126               ierror = 220              if(mask(iv,ik,1).eq.1)call cncomp(iv,ik,iflag)
127            endif              if(iflag.ne.0)then
128          enddo                 ima=ima+1
129                   mask_vk_ev(iv,ik)=0
130                   ierror = 220
131                endif
132                call stripmask(iv,ik)      !compute mask(i,j,k), combining VA1-masks
133                
134             enddo
135     100     format(' * WARNING * Event ',i7,' view',i3,': VK MASK ',24i1)
136             if(ima.ne.0.and.debug)write(*,100)eventn(1),iv
137         $        ,(mask_vk_ev(iv,ik),ik=1,nva1_view)
138        enddo        enddo
 c      if(good1.eq.0)then  
 c         ierror = 220  
 c      endif  
139    
140        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
141    
142  c---------------------------------------------  c---------------------------------------------
143  c     loops on views, VA1 and strips,  c     loops on views, VA1 and strips,
144  c     and computes strips signals using  c     and computes strips signals using
145  c     badstrip, pedestals, and  c     badstrip, pedestals, and
146  c     sigma informations from histograms  c     sigma informations from histograms
147  c---------------------------------------------  c---------------------------------------------
       flag_shower = .false.  
148        ind=1                     !clsignal array index        ind=1                     !clsignal array index
149    
150          if(debug)print*,'-- search clusters'
151        do iv=1,nviews            !loop on views        do iv=1,nviews            !loop on views
152          do is=1,nstrips_view    !loop on strips (1)          do is=1,nstrips_view    !loop on strips (1)
153            if(mod(iv,2).eq.1) then            if(mod(iv,2).eq.1) then
# Line 150  C===  > Y view Line 159  C===  > Y view
159       $           *mask(iv,nvk(is),nst(is))       $           *mask(iv,nvk(is),nst(is))
160              clinclcut(is)=incuty*sigma(iv,nvk(is),nst(is))              clinclcut(is)=incuty*sigma(iv,nvk(is),nst(is))
161       $           *mask(iv,nvk(is),nst(is))       $           *mask(iv,nvk(is),nst(is))
162  ccc            print*,"value(",is,")(reduction)= ",value(is)              sat(is)=0
163                if( adc(iv,nvk(is),nst(is)).lt.adc_saty )sat(is)=1
164            else                        else            
165  C===  > X view  C===  > X view
166              value(is)= (DBLE(adc(iv,nvk(is),nst(is)))              value(is)= (DBLE(adc(iv,nvk(is),nst(is)))
# Line 160  C===  > X view Line 170  C===  > X view
170       $           *mask(iv,nvk(is),nst(is))       $           *mask(iv,nvk(is),nst(is))
171              clinclcut(is)=incutx*sigma(iv,nvk(is),nst(is))              clinclcut(is)=incutx*sigma(iv,nvk(is),nst(is))
172       $           *mask(iv,nvk(is),nst(is))       $           *mask(iv,nvk(is),nst(is))
173                sat(is)=0
174                if( adc(iv,nvk(is),nst(is)).gt.adc_satx )sat(is)=1
175            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)  
176          enddo                   !end loop on strips (1)          enddo                   !end loop on strips (1)
177          call search_cluster(iv)          call search_cluster(iv)
178  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  
179          if(.not.flag_shower)then          if(.not.flag_shower)then
180             call save_cluster(iv)             call save_cluster(iv)
181               if(debug)print*,'view ',iv,' #clusters ', nclstr_view
182          else          else
183             fshower(iv) = 1             fshower(iv) = 1
184             GOOD1(DSPn) = 11  c           GOOD1(DSPnumber(iv)) = 11 !AHAHAHAHA ORRORE!!
185               GOOD1(iv) = 11
186     101       format(' * WARNING * Event ',i7,' view',i3
187         $          ,' #clusters > ',i5,' --> MASKED')
188               if(debug)write(*,101)eventn(1),iv,nclstrmax_view
189          endif          endif
190        enddo                     ! end loop on views        enddo                     ! end loop on views
191        do iv=1,nviews        do iv=1,nviews
192          do ik=1,nva1_view          do ik=1,nva1_view
193            cnev(iv,ik)  = cn(iv,ik) !assigns computed CN to ntuple variables            cnev(iv,ik)    = cn(iv,ik) !assigns computed CN to ntuple variables
194            cnnev(iv,ik) = cnn(iv,ik) !assigns computed CN to ntuple variables            cnrmsev(iv,ik) = cnrms(iv,ik) !assigns computed CN to ntuple variables
195  ccc          print*,"cnev(",iv,",",ik,")(reduction)= ",cnev(iv,ik)            cnnev(iv,ik)   = cnn(iv,ik) !assigns computed CN to ntuple variables
196          enddo          enddo
197        enddo        enddo
198  C---------------------------------------------  C---------------------------------------------
# Line 200  C--------------------------------------- Line 205  C---------------------------------------
205        do iv = 1,nviews        do iv = 1,nviews
206           ngood = ngood + good1(iv)           ngood = ngood + good1(iv)
207        enddo        enddo
208        if(ngood.ne.0)print*,'* WARNING * LEVEL1 event status: '        if(debug.and.ngood.ne.0)print*,'* WARNING * Event ',eventn(1)
209         $     ,':LEVEL1 event status: '
210       $     ,(good1(i),i=1,nviews)       $     ,(good1(i),i=1,nviews)
211  c------------------------------------------------------------------------  c------------------------------------------------------------------------
212  c  c
# Line 242  c      good1 = 0 Line 248  c      good1 = 0
248           indmax(ic) = 0           indmax(ic) = 0
249           maxs(ic) = 0           maxs(ic) = 0
250           mult(ic) = 0                     mult(ic) = 0          
251           dedx(ic) = 0           sgnl(ic) = 0
252           whichtrack(ic) = 0           whichtrack(ic) = 0     !assigned @ level2
253    
254        enddo        enddo
255        do id=1,maxlength         !???        do id=1,maxlength         !???
# Line 263  c        crc1(iv)=0 Line 269  c        crc1(iv)=0
269                
270        return        return
271        end        end
272    
273  *---***---***---***---***---***---***---***---***  *---***---***---***---***---***---***---***---***
274  *  *
275  *  *
# Line 283  c        crc1(iv)=0 Line 290  c        crc1(iv)=0
290    
291  c     local variables  c     local variables
292        integer rmax,lmax         !estremi del cluster        integer rmax,lmax         !estremi del cluster
293        integer rstop,lstop       !per decidere quali strip includere nel cluster        integer rstop,lstop      
294                                  ! oltre il seed        integer first,last
295        integer first,last,diff   !per includere le strip giuste... !???        integer fsat,lsat
   
       integer multtemp          !temporary multiplicity variable  
296    
297        external nst        external nst
298    
 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)  
   
299        iseed=-999                !cluster seed index initialization        iseed=-999                !cluster seed index initialization
300    
301          inext=-999                !index where to start new cluster search
302    
303          flag_shower = .false.
304        nclstr_view=0        nclstr_view=0
305    
306        do jl=1,nladders_view     !1..3 !loops on ladders        do jl=1,nladders_view              !1..3 !loops on ladders
307           first=1+nstrips_ladder*(jl-1) !1,1025,2049  
308           last=nstrips_ladder*jl !1024,2048,3072           first = 1+nstrips_ladder*(jl-1) !1,1025,2049
309  c     X views have 1018 strips instead of 1024           last  = nstrips_ladder*jl       !1024,2048,3072
310    
311    *        X views have 1018 strips instead of 1024
312           if(mod(iv,2).eq.0) then           if(mod(iv,2).eq.0) then
313              first=first+3              first=first+3
314              last=last-3              last=last-3
# Line 332  c     X views have 1018 strips instead o Line 316  c     X views have 1018 strips instead o
316    
317           do is=first,last       !loop on strips in each ladder           do is=first,last       !loop on strips in each ladder
318    
319              if(is.le.iseed+1) goto 220  c---------------------------------------------
320  *******************************************************  c     new-cluster search starts at index inext
321  *     Elena 08/2006  c---------------------------------------------
322  * 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  
 *******************************************************  
323    
324              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)  
325  c-----------------------------------------  c-----------------------------------------
326  c     possible SEED...  c     possible SEED...
327  c-----------------------------------------  c-----------------------------------------
328                 itemp=is                 itemp = is
329                   fsat = 0         ! first saturated strip
330                   lsat = 0         ! last saturated strip
331                 if(itemp.eq.last) goto 230 !estremo...                 if(itemp.eq.last) goto 230 !estremo...
332  ****************************************************  c              ------------------------                
333  *     modificato da Elena (08/2006) per salvare  c              search for first maximum
334  *     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  
335                 do while(                 do while(
336       $                   value(itemp).le.value(itemp+1)       $                   value(itemp).le.value(itemp+1)
337       $              .and.value(itemp+1).gt.clseedcut(itemp+1))       $              .and.value(itemp+1).gt.clseedcut(itemp+1))
338                    itemp=itemp+1                    itemp = itemp+1
339                    if(itemp.eq.last) goto 230 !stops if reaches last strip                    if(itemp.eq.last)   goto 230 !stops if reaches last strip
340                      if(sat(itemp).eq.1) goto 230 !stop if reaches a saturated strip
341                 enddo            ! of the ladder                 enddo            ! of the ladder
342   230           continue   230           continue
343  c-----------------------------------------  c              -----------------------------
344    c              check if strips are saturated
345    c              -----------------------------    
346                   if( sat(itemp).eq.1 )then
347                      fsat = itemp
348                      lsat = itemp
349                      if(itemp.eq.last) goto 231 !estremo...
350                      do while( sat(itemp+1).eq.1 )
351                         itemp = itemp+1
352                         lsat = itemp
353                         if(itemp.eq.last)   goto 231 !stops if reaches last strip
354                      enddo                  
355                   endif
356     231           continue
357    c---------------------------------------------------------------------------
358  c     fownd SEED!!!  c     fownd SEED!!!
359  c-----------------------------------------  c     (if there are saturated strips, the cluster is centered in the middle)
360                 iseed=itemp      c---------------------------------------------------------------------------
361  c----------------------------------------------------------                 if(fsat.eq.0.and.lsat.eq.0)then
362                      iseed = itemp ! <<< SEED
363                   else
364                      iseed = int((lsat+fsat)/2) ! <<< SEED
365    c$$$                  print*,'saturated strips ',fsat,lsat,iseed
366    c$$$                  print*,'--> ',(value(ii),ii=fsat,lsat)
367                   endif    
368    c---------------------------------------------------------------
369  c     after finding a cluster seed, checks also adjacent strips,  c     after finding a cluster seed, checks also adjacent strips,
370  C     and marks the ones exceeding clinclcut  C     and tags the ones exceeding clinclcut
371  c----------------------------------------------------------  c---------------------------------------------------------------
372                 ir=iseed         !indici destro                 ir=iseed         !indici destro
373                 il=iseed         ! e sinistro                 il=iseed         ! e sinistro
374                                
# Line 403  c--------------------------------------- Line 379  c---------------------------------------
379                 lstop=0          ! inclusion loop                 lstop=0          ! inclusion loop
380    
381                 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
382                    ir=ir+1       !position index for strips on right side of  
383                                  ! cluster seed  
384                    il=il-1       !and for left side                    ir=ir+1       !index for right side
385                      il=il-1       !index for left side
386  c------------------------------------------------------------------------  c------------------------------------------------------------------------
387  c     checks for last or first strip of the ladder  c     checks for last or first strip of the ladder
388  c------------------------------------------------------------------------  c------------------------------------------------------------------------
389                    if(ir.gt.last) then !when index goes beyond last strip                    if( ir.gt.last  ) rstop = 1                      
390                       rstop=1    ! of the ladder, change rstop flag in order                    if( il.lt.first ) lstop = 1
                                 ! to "help" exiting from loop  
                   endif  
391                                        
                   if(il.lt.first) then !idem when index goes beyond  
                      lstop=1    ! first strip of the ladder  
                   endif  
                     
 c------------------------------------------------------------------------  
 c     check for clusters including more than nclstrp strips  
392  c------------------------------------------------------------------------  c------------------------------------------------------------------------
393                    if((rmax-lmax+1).ge.nclstrp) then  c     add strips exceeding inclusion cut
                      goto 210   !exits inclusion loop:  
                                 ! lmax and rmax maintain last value  
                                 ! NB .ge.!???  
                   endif  
 c------------------------------------------------------------------------  
 c     marks strips exceeding inclusion cut  
394  c------------------------------------------------------------------------  c------------------------------------------------------------------------
395                    if(rstop.eq.0) then !if last strip of the ladder or last                    if( (rmax-lmax+1).ge.nclstrp )goto 210   !exits inclusion loop
396                                  ! over-cut strip has not been reached  
397                       if(value(ir).gt.clinclcut(ir)) then !puts in rmax the                    if(rstop.eq.0) then !if right cluster morder has not been reached
398                          rmax=ir ! last right over-cut strip                       if(value(ir).gt.clinclcut(ir)) then
399                            rmax=ir !include a strip on the right
400                       else                       else
401                          rstop=1 !otherwise cluster ends on right and rstop                          rstop=1 !cluster right end
402                       endif      ! flag=1 signals it                       endif    
403                    endif                    endif
404                    if(lstop.eq.0) then  
405                      if( (rmax-lmax+1).ge.nclstrp )goto 210   !exits inclusion loop
406    
407                      if(lstop.eq.0) then !if left cluster morder has not been reached
408                       if(value(il).gt.clinclcut(il)) then                       if(value(il).gt.clinclcut(il)) then
409                          lmax=il                          lmax=il !include a strip on the left
410                       else                       else
411                          lstop=1                          lstop=1 !cluster left end
412                       endif                       endif
413                    endif                    endif
414    
415                 enddo            !ends strip inclusion loop                 enddo            !ends strip inclusion loop
416                   goto 211
417   210           continue         !jumps here if more than nclstrp have been included   210           continue         !jumps here if more than nclstrp have been included
418                        c               print*,'>>> nclstrp! '
419                 multtemp=rmax-lmax+1 !stores multiplicity in temp   211           continue
420                                  ! variable. NB rmax and lmax can change later in  c-----------------------------------------
421                                  ! order to include enough strips to calculate eta3  c     end of inclusion loop!
422                                  ! and eta4. so mult is not always equal to cllength  c-----------------------------------------
423  c------------------------------------------------------------------------                
424  c     NB per essere sicuro di poter calcolare eta3 e eta4 devo includere  c------------------------------------------------------------------------
425  c     sempre e comunque le 2 strip adiacenti al cluster seed e quella  c     adjust the cluster in order to have at least a strip around the seed
426  c     adiacente ulteriore dalla parte della piu' alta fra queste due  c------------------------------------------------------------------------
427  c     (vedi oltre...)!???                 if(iseed.eq.lmax.and.lmax.ne.first)then
428  c------------------------------------------------------------------------                    lmax = lmax-1
429                      if( (rmax-lmax+1).gt.nclstrp )rmax=rmax-1
430  c     nel caso di estremi del ladder...!???                 endif
431                   if(iseed.eq.rmax.and.rmax.ne.last )then
432  c     ho meno di 4 strip nel cluster --> se sono sui bordi o quasi del ladder                    rmax = rmax+1
433  c     costruisco il cluster ad hoc e poi esco, se non sono sui bordi o quasi                    if( (rmax-lmax+1).gt.nclstrp )lmax=lmax+1
 c     vado oltre (aggiungero' quindi strip a sx e dx in modo da poter calcolare  
 c     eta3e4)  
                if((rmax-lmax+1).lt.4) then  
   
                   if(iseed.eq.first) then !estremi...  
                      rmax=iseed+2 !NB in questo modo puo' anche capitare di  
                      lmax=iseed ! includere strip sotto taglio di inclusione  
                      goto 250   ! che non serviranno per eta3e4!???  
                   endif  
                     
                   if(iseed.eq.last) then !estremi...  
                      rmax=iseed  
                      lmax=iseed-2 !NB 2 e non 3, perche' altrimenti sarei in  
                      goto 250   ! ((rmax-lmax+1).lt.4).eq.false. !???  
                   endif         !NMB questo e' l'unico caso di cllength=3!???  
                     
                   if(iseed.eq.first+1) then !quasi estremi...  
                      rmax=iseed+2  
                      lmax=iseed-1  
                      goto 250  
                   endif  
                   if(iseed.eq.last-1) then  
                      rmax=iseed+1  
                      lmax=iseed-2  
                      goto 250  
                   endif  
 c     se ho 4 o piu' strip --> se sono sui bordi esco, se sono sui quasi bordi  
 c     includo la strip del bordo  
                else  
                     
                   if(iseed.eq.first) goto 250 !estremi... non includo altro                    
                   if(iseed.eq.last) goto 250  
                   if(iseed.eq.first+1) then !quasi estremi... mi assicuro di  
                      lmax=first ! avere le strip adiacenti al seed  
                      if((rmax-lmax+1).gt.nclstrp) rmax=rmax-1 !NB effetto  
                      goto 250   ! coperta: se la lunghezza del cluster era gia'  
                   endif         ! al limite (nclstrp), per poter aggiungere questa  
                                 ! strip a sinistra devo toglierne una a destra...!???  
                   if(iseed.eq.last-1) then  
                      rmax=last  
                      if((rmax-lmax+1).gt.nclstrp) lmax=lmax+1  
                      goto 250  
                   endif                    
434                 endif                 endif
435    
436  c------------------------------------------------------------------------  c------------------------------------------------------------------------
437  c     be sure to include in the cluster the cluster seed with its 2 adjacent  c     adjust the cluster in order to store a minimum number of strips
438  c     strips, and the one adjacent to the greatest between this two strip, as the  c------------------------------------------------------------------------
439  c     fourth one. if the strips have the same value (!) the fourth one is chosen                 do while( (rmax-lmax+1).lt.nclstrpmin )
440  c     as the one having the greatest value between the second neighbors  
441  c------------------------------------------------------------------------                    vl = -99999
442                 if(value(iseed+1).eq.value(iseed-1)) then                    vr = -99999
443                    if(value(iseed+2).ge.value(iseed-2)) then !??? qui cmq c'e'                    if(lmax-1.ge.first) vl = value(lmax-1)
444                       diff=(iseed+2)-rmax                    if(rmax+1.le.last ) vr = value(rmax+1)
445                       if(diff.gt.0) then                    if(vr.ge.vl) then
446                          rmax=rmax+diff                       rmax = rmax+1
447                          if((rmax-lmax+1).gt.nclstrp) then                    else  
448                             lmax=rmax-nclstrp+1                       lmax = lmax-1
                         endif  
                      endif  
                      diff=(iseed-1)-lmax  
                      if(diff.lt.0) then  
                         lmax=lmax+diff  
                         if((rmax-lmax+1).gt.nclstrp) then  
                            rmax=lmax+nclstrp-1  
                         endif  
                      endif  
                   else  
                      diff=(iseed-2)-lmax  
                      if(diff.lt.0) then  
                         lmax=lmax+diff  
                         if((rmax-lmax+1).gt.nclstrp) then  
                            rmax=lmax+nclstrp-1  
                         endif  
                      endif  
                      diff=(iseed+1)-rmax  
                      if(diff.gt.0) then  
                         rmax=rmax+diff  
                         if((rmax-lmax+1).gt.nclstrp) then  
                            lmax=rmax-nclstrp+1  
                         endif  
                      endif  
                   endif  
                elseif(value(iseed+1).gt.value(iseed-1)) then  
 c     !??? sposto il limite del cluster a destra per includere sempre le strip  
 c     necessarie al calcolo di eta-i  
 c     se il cluster diventa  troppo lungo lo accorcio a sinistra per avere non piu'  
 c     di nclstrp (in questo caso sono sicuro di aver gia' incluso le strip  
 c     necessarie al calcolo di eta-i a sinistra, quindi se voglio posso uscire)  
                   diff=(iseed+2)-rmax  
                   if(diff.gt.0) then  
                      rmax=rmax+diff  
                      if((rmax-lmax+1).gt.nclstrp) then  
                         lmax=rmax-nclstrp+1  
 c     goto 250  
                      endif  
                   endif  
                   diff=(iseed-1)-lmax  
                   if(diff.lt.0) then  
                      lmax=lmax+diff  
                      if((rmax-lmax+1).gt.nclstrp) then  
                         rmax=lmax+nclstrp-1  
 c     goto 250 !inutile!???  
                      endif  
                   endif  
                else  
                   diff=(iseed-2)-lmax  
                   if(diff.lt.0) then  
                      lmax=lmax+diff  
                      if((rmax-lmax+1).gt.nclstrp) then  
                         rmax=lmax+nclstrp-1  
 c     goto 250  
                      endif  
                   endif  
                   diff=(iseed+1)-rmax  
                   if(diff.gt.0) then  
                      rmax=rmax+diff  
                      if((rmax-lmax+1).gt.nclstrp) then  
                         lmax=rmax-nclstrp+1  
 c     goto 250 !inutile!???  
                      endif  
449                    endif                    endif
450                 endif                    
451   250           continue                 enddo
452    
453  c--------------------------------------------------------  c--------------------------------------------------------
454  c     fills cluster variables  c     store cluster info
455  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  
456                 nclstr_view = nclstr_view + 1 !cluster number                 nclstr_view = nclstr_view + 1 !cluster number
457  c               print*,'view ',iv,' -- search_cluster -- nclstr_view: '  
 c     $              ,nclstr_view  
458                 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:
459                    if(verbose) print*,'Event ',eventn(1),  c$$$                  if(verbose) print*,'Event ',eventn(1),
460       $                 ': more than ',nclstrmax_view  c$$$     $                 ': more than ',nclstrmax_view
461       $                 ,' clusters on view ',iv  c$$$     $                 ,' clusters on view ',iv
 c                  good1=0       ! event  
 c                  nclstr1=0  
 c                  totCLlength=0  
462                    flag_shower = .true.                    flag_shower = .true.
463                    goto 2000                    goto 2000
464                 endif                 endif
465    
466  c               view(nclstr1)   = iv !vista del cluster                 ladder_view(nclstr_view) = nld(iseed,iv)
467                 ladder_view(nclstr_view) = nld(iseed,iv) !ladder a cui appartiene il cluster seed                 maxs_view(nclstr_view)   = iseed
468                 maxs_view(nclstr_view)   = iseed !strip del cluster seed                 mult_view(nclstr_view)   = rmax-lmax+1
                mult_view(nclstr_view)   = multtemp !molteplicita'  
469                 rmax_view(nclstr_view)   = rmax                 rmax_view(nclstr_view)   = rmax
470                 lmax_view(nclstr_view)   = lmax                 lmax_view(nclstr_view)   = lmax
471    
472    c$$$               if(rmax-lmax+1.gt.25)
473    c$$$     $              print*,'view ',iv
474    c$$$     $              ,' cl ',nclstr_view,' mult ',rmax-lmax+1
475    c------------------------------------------------------------------------
476    c     search for a double peak inside the cluster                                                                                                            
477    c------------------------------------------------------------------------
478                   inext = rmax+1   !<< index where to start new-cluster search
479                  
480                   vmax = 0
481                   vmin = value(iseed)
482                   imax = iseed
483                   imin = iseed
484                   do iss = max(iseed+1,lsat+1),rmax
485                      if( value(iss).lt.vmin )then
486                         if( imax.ne.iseed )goto 221 !found dowble peek
487                         imin = iss
488                         vmin = value(iss)
489                      else
490                         delta = value(iss) - value(imin)
491                         cut = sqrt(clinclcut(iss)**2 + clinclcut(imin)**2)
492                         if(
493         $                    delta.gt.cut .and.
494         $                    value(iss).gt.clseedcut(iss).and.
495         $                    .true.)then
496                            if( value(iss).gt.vmax )then                        
497                               imax = iss
498                               vmax = value(iss)
499                            else
500                               goto 221 !found dowble peek
501                            endif
502                         endif
503                      endif
504                   enddo
505     221           continue
506                  
507                   if(imax.gt.iseed)then
508                      inext = imax    !<< index where to start new-cluster search
509    c$$$                  print*,'--- double peek ---'
510    c$$$                  print*,(value(ii),ii=lmax,rmax)
511    c$$$                  print*,'seed ',iseed,' imin ',imin,' imax ',imax
512                   endif
513  c--------------------------------------------------------  c--------------------------------------------------------
514  c  c
515  c--------------------------------------------------------  c--------------------------------------------------------
# Line 694  c        posizione del cluster seed nell Line 559  c        posizione del cluster seed nell
559                    
560           CLlength      = rmax_view(ic) - lmax_view(ic) + 1 !numero di strip salvate           CLlength      = rmax_view(ic) - lmax_view(ic) + 1 !numero di strip salvate
561           totCLlength   = totCLlength + CLlength           totCLlength   = totCLlength + CLlength
562           dedx(nclstr1) = 0           sgnl(nclstr1) = 0
563           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
564    
565              clsignal(ind) = value(j) ! clsignal array              clsignal(ind) = value(j) ! clsignal array
# Line 710  c            clped(ind)   = pedestal(iv, Line 575  c            clped(ind)   = pedestal(iv,
575              ind=ind+1              ind=ind+1
576  c     if(value(j).gt.0)  c     if(value(j).gt.0)
577              if(value(j).gt.clinclcut(j))              if(value(j).gt.clinclcut(j))
578       $           dedx(nclstr1) = dedx(nclstr1) + value(j) !cluster charge       $           sgnl(nclstr1) = sgnl(nclstr1) + value(j) !cluster charge
579           enddo           enddo
580    
581  c         print*,'view ',iv,' -- save_cluster -- nclstr1: '  c         print*,'view ',iv,' -- save_cluster -- nclstr1: '
582  c     $        ,nclstr1,maxs(nclstr1),mult(nclstr1),dedx(nclstr1)  c     $        ,nclstr1,maxs(nclstr1),mult(nclstr1),sgnl(nclstr1)
583                    
584        enddo        enddo
585                
# Line 729  c     $        ,nclstr1,maxs(nclstr1),mu Line 594  c     $        ,nclstr1,maxs(nclstr1),mu
594  *---***---***---***---***---***---***---***---***  *---***---***---***---***---***---***---***---***
595    
596    
597        subroutine stripmask  c$$$      subroutine stripmask
598    c$$$
599    c$$$*     this routine set va1 and single-strip masks,
600    c$$$*     on the basis of the VA1 mask saved in the DB
601    c$$$*
602    c$$$*     mask(nviews,nva1_view,nstrips_va1) !strip mask
603    c$$$*     mask_vk(nviews,nva1_view)          !VA1 mask
604    c$$$*
605    c$$$      include 'commontracker.f'
606    c$$$      include 'level1.f'
607    c$$$      include 'common_reduction.f'
608    c$$$      include 'calib.f'
609    c$$$
610    c$$$*     init mask
611    c$$$      do iv=1,nviews
612    c$$$         do ivk=1,nva1_view
613    c$$$            do is=1,nstrips_va1
614    c$$$c               mask(iv,ivk,is) = mask_vk(iv,ivk)
615    c$$$               if( mask_vk(iv,ivk) .ne. -1)then
616    c$$$                  mask(iv,ivk,is) = 1
617    c$$$     $                 * mask_vk(iv,ivk)     !from DB
618    c$$$     $                 * mask_vk_ev(iv,ivk)  !from <SIG>
619    c$$$     $                 * mask_vk_run(iv,ivk) !from CN
620    c$$$               else
621    c$$$                  mask(iv,ivk,is) = -1
622    c$$$     $                 * mask_vk(iv,ivk)     !from DB
623    c$$$     $                 * mask_vk_ev(iv,ivk)  !from CN
624    c$$$               endif
625    c$$$            enddo
626    c$$$         enddo
627    c$$$      enddo
628    c$$$
629    c$$$
630    c$$$      return
631    c$$$      end
632    
633          subroutine stripmask(iv,ivk)
634    
635  *     this routine set va1 and single-strip masks,  *     this routine set va1 and single-strip masks,
636  *     on the basis of the VA1 mask saved in the DB  *     on the basis of the VA1 mask saved in the DB
# Line 743  c     $        ,nclstr1,maxs(nclstr1),mu Line 644  c     $        ,nclstr1,maxs(nclstr1),mu
644        include 'calib.f'        include 'calib.f'
645    
646  *     init mask  *     init mask
647        do iv=1,nviews        do is=1,nstrips_va1
648           do ivk=1,nva1_view           if( mask_vk(iv,ivk) .ne. -1)then            
649              do is=1,nstrips_va1              mask(iv,ivk,is) = 1
650  c               mask(iv,ivk,is) = mask_vk(iv,ivk)       $           * mask_vk(iv,ivk) !from DB
651                 mask(iv,ivk,is) = mask_vk(iv,ivk) * mask_vk_ev(iv,ivk)                     $           * mask_vk_ev(iv,ivk) !from <SIG>
652              enddo       $           * mask_vk_run(iv,ivk) !from CN
653           enddo           else
654                mask(iv,ivk,is) = -1
655         $           * mask_vk(iv,ivk) !from DB
656         $           * mask_vk_ev(iv,ivk) !from CN
657             endif
658        enddo        enddo
659          
660          
661        return        return
662        end        end
   

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

  ViewVC Help
Powered by ViewVC 1.1.23