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

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

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

revision 1.1.1.1 by mocchiut, Fri May 19 13:15:55 2006 UTC revision 1.4 by pam-fi, Thu May 24 16:45:48 2007 UTC
# Line 132  c--------------------------------------- Line 132  c---------------------------------------
132    
133        end        end
134    
 c$$$        
 c$$$c$$$c------------------------------------------------------------------------  
 c$$$c$$$  
 c$$$c$$$c     NB: le coordinate in mech_pos.dat sono calcolate a partire da alcuni dati  
 c$$$c$$$c     contenuti in commontracker.f. forse si puo' evitare mech_pos.dat e mettere  
 c$$$c$$$c     tutto in commontracker.f  
 c$$$c$$$  
 c$$$c$$$  
 c$$$c$$$      subroutine mech_sensor    !it reads sensors coordinates (in PAMELA reference  
 c$$$c$$$                                ! frame) from a text file and it uses them to fill  
 c$$$c$$$                                ! x/y/z_mech_sensor variables, taking into account  
 c$$$c$$$                                ! last plane inversion  
 c$$$c$$$  
 c$$$c$$$      include './commontracker.f'  
 c$$$c$$$      include './common_tracks.f'  
 c$$$c$$$  
 c$$$c$$$      real xvec(nladders_view),yvec(2),zvec(nplanes)  
 c$$$c$$$  
 c$$$c$$$      integer id                !file identifier  
 c$$$c$$$      logical od                !.true. if the specified unit is connected to a file  
 c$$$c$$$  
 c$$$c$$$      do id=20,100,1            !opens the file using a free file id  
 c$$$c$$$        inquire (id, opened=od)  
 c$$$c$$$        if(.not.od) goto 666  
 c$$$c$$$      enddo        
 c$$$c$$$ 666  continue  
 c$$$c$$$  
 c$$$c$$$      open(id,FILE='../common/mech_pos.dat') !sensors centres coordinates in mm in  
 c$$$c$$$                                ! PAMELA reference frame:  
 c$$$c$$$                                ! the first plane is the one with lowest Z (the one  
 c$$$c$$$                                ! nearest the calorimeter)  
 c$$$c$$$                                ! the first ladder is the one with lowest X (the  
 c$$$c$$$                                ! one on which the first X strip is)  
 c$$$c$$$                                ! the first sensor is the one with lowest Y (the  
 c$$$c$$$                                ! one on which the first Y strip is) for planes  
 c$$$c$$$                                ! 2..6. for plane 1 the first sensor has higher Y  
 c$$$c$$$  
 c$$$c$$$      read(20,*) xvec  
 c$$$c$$$      read(20,*) yvec  
 c$$$c$$$      read(20,*) zvec  
 c$$$c$$$  
 c$$$c$$$      do i=1,nplanes  
 c$$$c$$$        do j=1,nladders_view  
 c$$$c$$$          do k=1,2  
 c$$$c$$$            x_mech_sensor(i,j,k)=xvec(j)  
 c$$$c$$$            y_mech_sensor(i,j,k)=yvec(k)  
 c$$$c$$$            z_mech_sensor(i,j,k)=zvec(i)  
 c$$$c$$$            if(i.eq.1) then     !y coordinates of first plane (11th view) are  
 c$$$c$$$              y_mech_sensor(i,j,k)=-yvec(k) ! exchanged due to last plane inversion  
 c$$$c$$$            endif                
 c$$$c$$$          enddo  
 c$$$c$$$        enddo  
 c$$$c$$$      enddo  
 c$$$c$$$  
 c$$$c$$$      close(id)  
 c$$$c$$$  
 c$$$c$$$  
 c$$$c$$$c$$$  ! *** INIZIO DEBUG ***  
 c$$$c$$$c$$$  do i=1,6  
 c$$$c$$$c$$$  do j=1,3  
 c$$$c$$$c$$$  do k=1,2  
 c$$$c$$$c$$$  c            print*,x_mech_sensor(1,j,k)  
 c$$$c$$$c$$$  print*,y_mech_sensor(i,j,k)  
 c$$$c$$$c$$$  c            print*,z_mech_sensor(i,j,k)  
 c$$$c$$$c$$$  enddo  
 c$$$c$$$c$$$  enddo  
 c$$$c$$$c$$$  print*,' '  
 c$$$c$$$c$$$  enddo  
 c$$$c$$$c$$$  ! *** FINE DEBUG ***  
 c$$$c$$$  
 c$$$c$$$        
 c$$$c$$$      return  
 c$$$c$$$      end  
   
   
 c$$$c------------------------------------------------------------------------  
 c$$$  
 c$$$c     NB: le coordinate in mech_pos.dat sono calcolate a partire da alcuni dati  
 c$$$c     contenuti in commontracker.f. forse si puo' evitare mech_pos.dat e mettere  
 c$$$c     tutto in commontracker.f  
 c$$$  
 c$$$  
 c$$$      subroutine mech_sensor      
 c$$$c     !it reads sensors coordinates (in PAMELA reference  
 c$$$c     ! frame) from a text file and it uses them to fill  
 c$$$c     ! x/y/z_mech_sensor variables, taking into account  
 c$$$c     ! last plane inversion  
 c$$$  
 c$$$      include './commontracker.f'  
 c$$$      include './common_tracks.f'  
 c$$$  
 c$$$      real xvec(nladders_view),yvec(2),zvec(nplanes)  
 c$$$  
 c$$$      integer id                !file identifier  
 c$$$      logical od                !.true. if the specified unit is connected to a file  
 c$$$  
 c$$$      do id=20,100,1            !opens the file using a free file id  
 c$$$        inquire (id, opened=od)  
 c$$$        if(.not.od) goto 666  
 c$$$      enddo        
 c$$$ 666  continue  
 c$$$  
 c$$$c      open(id,FILE='../common/mech_pos.dat') !sensors centres coordinates in mm in  
 c$$$c      open(id,FILE='source/common/mech_pos.dat')  
 c$$$c      call system('cp $TRK_GRND/source/common/mech_pos.dat .')  
 c$$$      print *,'Opening file: mech_pos.dat'  
 c$$$      open(id,FILE='./bin-aux/mech_pos.dat',IOSTAT=iostat)  
 c$$$c     !sensors centres coordinates in mm in  
 c$$$c     ! PAMELA reference frame:  
 c$$$c     ! the first plane is the one with lowest Z (the one  
 c$$$c     ! nearest the calorimeter)  
 c$$$c     ! the first ladder is the one with lowest X (the  
 c$$$c     ! one on which the first X strip is)  
 c$$$c     ! the first sensor is the one with lowest Y (the  
 c$$$c     ! one on which the first Y strip is) for planes  
 c$$$c     ! 2..6. for plane 1 the first sensor has higher Y  
 c$$$        
 c$$$      if(iostat.ne.0)then  
 c$$$         print*,'MECH_SENSOR: *** Error in opening file ***'  
 c$$$         return  
 c$$$      endif  
 c$$$  
 c$$$      read(id,*) xvec  
 c$$$      read(id,*) yvec  
 c$$$      read(id,*) zvec  
 c$$$  
 c$$$      do i=1,nplanes  
 c$$$        do j=1,nladders_view  
 c$$$          do k=1,2  
 c$$$            x_mech_sensor(i,j,k)=xvec(j)  
 c$$$            y_mech_sensor(i,j,k)=yvec(k)  
 c$$$            z_mech_sensor(i,j,k)=zvec(i)  
 c$$$            if(i.eq.1) then     !y coordinates of first plane (11th view) are  
 c$$$              y_mech_sensor(i,j,k)=-yvec(k) ! exchanged due to last plane inversion  
 c$$$            endif                
 c$$$          enddo  
 c$$$        enddo  
 c$$$      enddo  
 c$$$  
 c$$$      close(id)  
 c$$$c      call system('rm -f mech_pos.dat')  
 c$$$  
 c$$$c$$$                                ! *** INIZIO DEBUG ***  
 c$$$c$$$      do i=1,6  
 c$$$c$$$c        do j=1,3  
 c$$$c$$$          do k=1,2  
 c$$$c$$$            j=1  
 c$$$c$$$c            print*,x_mech_sensor(1,j,k)  
 c$$$c$$$            print*,y_mech_sensor(i,j,k)  
 c$$$c$$$c            print*,z_mech_sensor(i,j,k)  
 c$$$c$$$          enddo  
 c$$$c$$$c        enddo  
 c$$$c$$$        print*,' '  
 c$$$c$$$      enddo  
 c$$$c$$$                                ! *** FINE DEBUG ***  
 c$$$  
 c$$$        
 c$$$      return  
 c$$$      end  
 c$$$  
 c$$$  
 c$$$c------------------------------------------------------------------------  
135    
136    
137        function coordsi(istrip,view)        function coordsi(istrip,view)
# Line 325  c     NB mettere il 1024 nel commontrack Line 163  c     NB mettere il 1024 nel commontrack
163                
164        if(mod(view,2).eq.0) then !X view        if(mod(view,2).eq.0) then !X view
165    
166          if((is.le.3).or.(is.ge.1022)) then !X has 1018 strips...  c        if((is.le.3).or.(is.ge.1022)) then !X has 1018 strips...
167            print*,'functions: WARNING: false X strip: strip ',is  c          print*,'functions: WARNING: false X strip: strip ',is
168          endif  c        endif
169    
170          is=is-3                 !4 =< is =< 1021 --> 1 =< is =< 1018          is=is-3                 !4 =< is =< 1021 --> 1 =< is =< 1018
171    
# Line 396  c     NB mettere il 1024 nel commontrack Line 234  c     NB mettere il 1024 nel commontrack
234                
235        if(mod(view,2).eq.0) then !X view        if(mod(view,2).eq.0) then !X view
236    
237          if((is.le.3).or.(is.ge.1022)) then !X has 1018 strips...  c        if((is.le.3).or.(is.ge.1022)) then !X has 1018 strips...
238            print*,'functions: WARNING: false X strip: strip ',is  c          print*,'functions: WARNING: false X strip: strip ',is
239          endif  c        endif
240    
241          is=is-3                 !4 =< is =< 1021 --> 1 =< is =< 1018          is=is-3                 !4 =< is =< 1021 --> 1 =< is =< 1018
242    
# Line 520  c$$$                                ! CO Line 358  c$$$                                ! CO
358    
359        dcoord=0.        dcoord=0.
360                
361          if(
362         $     ladder.lt.1.or.
363         $     ladder.gt.nladders_view.or.
364         $     sen.lt.1.or.
365         $     sen.gt.2.or.
366         $     view.lt.1.or.
367         $     view.gt.nviews.or.
368         $     .false.)then
369             print*,'dcoord ---> wrong input: ladder ',ladder
370         $        ,' sensor ',sen
371         $        ,' view ',view
372             return
373          endif
374    
375        sx=ladder        sx=ladder
376        sy=sen        sy=sen
377        sz=npl(view)        sz=npl(view)
378    
379    
380        if(mod(view,2).eq.0) then !X view        if(mod(view,2).eq.0) then !X view
381    
382          trasl=x_mech_sensor(sz,sx,sy) !in mm          trasl=x_mech_sensor(sz,sx,sy) !in mm
# Line 544  c$$$        endif Line 397  c$$$        endif
397    
398    
399  c------------------------------------------------------------------------  c------------------------------------------------------------------------
400          integer function nsatstrips(ic)
401    *--------------------------------------------------------------
402    *     this function returns the number of saturated strips
403    *     inside a cluster
404    *--------------------------------------------------------------
405          include 'commontracker.f'
406          include 'level1.f'
407          include 'calib.f'
408    
409          
410    
411          integer nsat
412          nsat = 0
413          iv=VIEW(ic)              
414          if(mod(iv,2).eq.1)incut=incuty
415          if(mod(iv,2).eq.0)incut=incutx
416          istart = INDSTART(IC)
417          istop  = TOTCLLENGTH
418          if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
419          do i = INDMAX(IC),istart,-1
420    c         cut  = incut*CLSIGMA(i)
421    c         if(CLSIGNAL(i).ge.cut)then
422             if(  (mod(iv,2).eq.1.and.CLADC(i).lt.ADCsatx)
423         $        .or.
424         $        (mod(iv,2).eq.0.and.CLADC(i).gt.ADCsaty) )then
425                nsat = nsat +1
426             else
427                goto 10
428             endif
429          enddo
430     10   continue
431          do i = INDMAX(IC)+1,istop
432    c         cut  = incut*CLSIGMA(i)
433    c         if(CLSIGNAL(i).ge.cut)then
434             if(  (mod(iv,2).eq.1.and.CLADC(i).lt.ADCsatx)
435         $        .or.
436         $        (mod(iv,2).eq.0.and.CLADC(i).gt.ADCsaty) )then
437                nsat = nsat +1
438             else
439                goto 20
440             endif
441          enddo
442     20   continue
443                
444          nsatstrips = nsat
445          return
446          end
447    
448    c------------------------------------------------------------------------
449          integer function nbadstrips(ncog,ic)
450    *--------------------------------------------------------------
451    *     this function returns the number of BAD strips
452    *     inside a cluster:
453    *     - if NCOG=0, the number BAD strips inside the whole cluster
454    *     are given, according to the cluster multiplicity
455    *    
456    *     - if NCOG>0, the number BAD strips is evaluated using NCOG
457    *     strips, even if they have a negative signal (according to Landi)
458    *--------------------------------------------------------------
459          include 'commontracker.f'
460          include 'level1.f'
461          include 'calib.f'
462    
463          integer nbad
464          nbad = 0
465    
466          if (ncog.gt.0) then
467    
468    *     --> signal of the central strip
469             sc = CLSIGNAL(INDMAX(ic)) !center
470    *     signal of adjacent strips
471             sl1 = -100000                !left 1
472             if(
473         $        (INDMAX(ic)-1).ge.INDSTART(ic)
474         $        )
475         $        sl1 = CLSIGNAL(INDMAX(ic)-1)
476            
477             sl2 = -100000                !left 2
478             if(
479         $        (INDMAX(ic)-2).ge.INDSTART(ic)
480         $        )
481         $        sl2 = CLSIGNAL(INDMAX(ic)-2)
482            
483             sr1 = -100000                !right 1
484             if(
485         $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
486         $        .or.
487         $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
488         $        )
489         $        sr1 = CLSIGNAL(INDMAX(ic)+1)
490            
491             sr2 = -100000                !right 2
492             if(
493         $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
494         $        .or.
495         $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
496         $        )
497         $        sr2 = CLSIGNAL(INDMAX(ic)+2)
498            
499             if(ncog.ge.1)nbad = nbad+1-CLBAD(INDMAX(ic))
500            
501             if(ncog.ge.2)then
502                if(sl1.gt.sr1.and.(INDMAX(ic)-1).ge.1)then
503                   nbad=nbad+1-CLBAD(INDMAX(ic)-1)      
504                elseif(sl1.le.sr1.and.(INDMAX(ic)+1).le.NCLSTR1)then
505                   nbad=nbad+1-CLBAD(INDMAX(ic)+1)          
506                endif
507             endif
508            
509             if(ncog.ge.3)then
510                if(sl1.gt.sr1.and.(INDMAX(ic)+1).le.NCLSTR1)then
511                   nbad=nbad+1-CLBAD(INDMAX(ic)+1)
512                elseif(sl1.le.sr1.and.(INDMAX(ic)-1).ge.1)then
513    c               if(INDMAX(ic)-1.eq.0)
514    c     $              print*,' ======= ',sl2,sl1,sc,sr1,sr2
515                   nbad=nbad+1-CLBAD(INDMAX(ic)-1)
516                endif
517             endif
518            
519             if(ncog.ge.4)then
520                if(sl2.gt.sr2.and.(INDMAX(ic)-2).ge.1)then
521                   nbad=nbad+1-CLBAD(INDMAX(ic)-2)      
522                elseif(sl2.le.sr2.and.(INDMAX(ic)+2).le.NCLSTR1)then
523                   nbad=nbad+1-CLBAD(INDMAX(ic)+2)          
524                endif
525             endif
526            
527    c         if(ncog.ge.5)then
528    c            print*,'function CLBAD(NCOG,IC) ==> WARNING!! NCOG=',NCOG
529    c     $           ,' not implemented'
530    c         endif
531                
532          elseif(ncog.eq.0)then
533    *     =========================
534    *     COG computation
535    *     =========================
536    
537    
538            
539             iv=VIEW(ic)        
540             if(mod(iv,2).eq.1)incut=incuty
541             if(mod(iv,2).eq.0)incut=incutx
542            
543             istart = INDSTART(IC)
544             istop  = TOTCLLENGTH
545             if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
546             nbad = 0
547    c$$$         do i=istart,istop
548    c$$$            cut  = incut*CLSIGMA(i)            
549    c$$$            if(CLSIGNAL(i).ge.cut)nbad = nbad +1 -CLBAD(i)
550    c$$$         enddo
551             do i = INDMAX(IC),istart,-1
552                cut  = incut*CLSIGMA(i)
553                if(CLSIGNAL(i).ge.cut)then
554                   nbad = nbad +1 -CLBAD(i)
555                else
556                   goto 10
557                endif
558             enddo
559     10      continue
560             do i = INDMAX(IC)+1,istop
561                cut  = incut*CLSIGMA(i)
562                if(CLSIGNAL(i).ge.cut)then
563                   nbad = nbad +1 -CLBAD(i)
564                else
565                   goto 20
566                endif
567             enddo
568     20      continue
569            
570          else
571            
572    c         print*,'function CLBAD(NCOG,IC) ==> WARNING!! NCOG=',NCOG
573    c     $        ,' not implemented'
574            
575    
576          endif
577    
578          nbadstrips = nbad
579    
580          return
581          end
Legend:
Removed from v.1.1.1.1  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.23