/[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.5 by bongi, Tue Sep 4 15:44:50 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    c------------------------------------------------------------------------
137    
138        function coordsi(istrip,view)        function coordsi(istrip,view)
139  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
# Line 325  c     NB mettere il 1024 nel commontrack Line 164  c     NB mettere il 1024 nel commontrack
164                
165        if(mod(view,2).eq.0) then !X view        if(mod(view,2).eq.0) then !X view
166    
167          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...
168            print*,'functions: WARNING: false X strip: strip ',is  c          print*,'functions: WARNING: false X strip: strip ',is
169          endif  c        endif
170    
171          is=is-3                 !4 =< is =< 1021 --> 1 =< is =< 1018          is=is-3                 !4 =< is =< 1021 --> 1 =< is =< 1018
172    
# Line 396  c     NB mettere il 1024 nel commontrack Line 235  c     NB mettere il 1024 nel commontrack
235                
236        if(mod(view,2).eq.0) then !X view        if(mod(view,2).eq.0) then !X view
237    
238          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...
239            print*,'functions: WARNING: false X strip: strip ',is  c          print*,'functions: WARNING: false X strip: strip ',is
240          endif  c        endif
241    
242          is=is-3                 !4 =< is =< 1021 --> 1 =< is =< 1018          is=is-3                 !4 =< is =< 1021 --> 1 =< is =< 1018
243    
# Line 436  ccccc coord1=(is-1)*p           !referre Line 275  ccccc coord1=(is-1)*p           !referre
275  c------------------------------------------------------------------------  c------------------------------------------------------------------------
276    
277    
278          double precision function dcoordsi(rstrip,view)
279    *
280    *     same as COORDSI, but accept a real value of strip!!!
281    *     and gives a double precision output
282    *
283    * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
284    *     it gives the strip coordinate in micrometers,
285    *     knowing the strip number (1..3072) and the view
286    *     number. the origin of the coordinate is on the
287    *     centre of the sensor the strip belongs to.
288    *     the axes directions are the same as in the PAMELA
289    *     reference frame (i.e.: the 11th view coordinate
290    *     direction has to be inverted here)
291    * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
292    
293    c      integer is,view,istrip
294    
295          integer view,is,istrip
296          real rstrip
297          double precision strip,stripladder,p
298          double precision edge,dim
299          double precision coord1
300    
301    
302          include 'commontracker.f'
303    
304    c     NB mettere il 1024 nel commontracker...!???
305    
306          strip=DBLE(rstrip)
307    
308          istrip = int(strip+0.5)   !istrip stores the closest integer to strip      
309    
310          is=istrip                 !it stores istrip number
311          is=mod(is-1,1024)+1       !it puts all clusters on a single ladder
312          
313          dcoordsi=0.
314          
315          if(mod(view,2).eq.0) then !X view
316    
317    c        if((is.le.3).or.(is.ge.1022)) then !X has 1018 strips...
318    c          print*,'functions: WARNING: false X strip: strip ',is
319    c        endif
320    
321            is=is-3                 !4 =< is =< 1021 --> 1 =< is =< 1018
322    
323            edge=edgeX
324            dim=SiDimX
325    
326          elseif(mod(view,2).eq.1) then !Y view
327    
328            edge=edgeY
329            dim=SiDimY
330    
331    c$$$        if(view.eq.11) then     !INVERSIONE!???
332    c$$$          is=1025-is
333    c$$$        endif
334    
335          endif
336    
337    
338          stripladder = DBLE(is)+(strip-DBLE(istrip))!cluster position relative to ladder
339          p=pitch(view)
340    
341    ccccc coord1=(is-1)*p           !referred to 1st sensor strip
342          coord1=(stripladder-1)*p  !referred to 1st sensor strip
343          coord1=coord1+edge        !referred to sensor edge
344          dcoordsi=coord1-dim/2     !referred to the centre of the sensor
345    
346          if(view.eq.11) then       !INVERSION: it puts y axis in the same direction for all views
347            dcoordsi=-dcoordsi
348          endif
349          
350          end
351    
352    
353    
354    c------------------------------------------------------------------------
355    
356    
357        function coord(coordsi,view,ladder,sen)        function coord(coordsi,view,ladder,sen)
358  *     it gives the coordinate in  *     it gives the coordinate in
359  *     micrometers, knowing the coordinate in the sensor  *     micrometers, knowing the coordinate in the sensor
# Line 520  c$$$                                ! CO Line 438  c$$$                                ! CO
438    
439        dcoord=0.        dcoord=0.
440                
441          if(
442         $     ladder.lt.1.or.
443         $     ladder.gt.nladders_view.or.
444         $     sen.lt.1.or.
445         $     sen.gt.2.or.
446         $     view.lt.1.or.
447         $     view.gt.nviews.or.
448         $     .false.)then
449             print*,'dcoord ---> wrong input: ladder ',ladder
450         $        ,' sensor ',sen
451         $        ,' view ',view
452             return
453          endif
454    
455        sx=ladder        sx=ladder
456        sy=sen        sy=sen
457        sz=npl(view)        sz=npl(view)
458    
459    
460        if(mod(view,2).eq.0) then !X view        if(mod(view,2).eq.0) then !X view
461    
462          trasl=x_mech_sensor(sz,sx,sy) !in mm          trasl=x_mech_sensor(sz,sx,sy) !in mm
# Line 544  c$$$        endif Line 477  c$$$        endif
477    
478    
479  c------------------------------------------------------------------------  c------------------------------------------------------------------------
480          integer function nsatstrips(ic)
481    *--------------------------------------------------------------
482    *     this function returns the number of saturated strips
483    *     inside a cluster
484    *--------------------------------------------------------------
485          include 'commontracker.f'
486          include 'level1.f'
487          include 'calib.f'
488    
489          
490    
491          integer nsat
492          nsat = 0
493          iv=VIEW(ic)              
494          if(mod(iv,2).eq.1)incut=incuty
495          if(mod(iv,2).eq.0)incut=incutx
496          istart = INDSTART(IC)
497          istop  = TOTCLLENGTH
498          if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
499          do i = INDMAX(IC),istart,-1
500    c         cut  = incut*CLSIGMA(i)
501    c         if(CLSIGNAL(i).ge.cut)then
502             if(  (mod(iv,2).eq.1.and.CLADC(i).lt.ADCsatx)
503         $        .or.
504         $        (mod(iv,2).eq.0.and.CLADC(i).gt.ADCsaty) )then
505                nsat = nsat +1
506             else
507                goto 10
508             endif
509          enddo
510     10   continue
511          do i = INDMAX(IC)+1,istop
512    c         cut  = incut*CLSIGMA(i)
513    c         if(CLSIGNAL(i).ge.cut)then
514             if(  (mod(iv,2).eq.1.and.CLADC(i).lt.ADCsatx)
515         $        .or.
516         $        (mod(iv,2).eq.0.and.CLADC(i).gt.ADCsaty) )then
517                nsat = nsat +1
518             else
519                goto 20
520             endif
521          enddo
522     20   continue
523                
524          nsatstrips = nsat
525          return
526          end
527    
528    c------------------------------------------------------------------------
529          integer function nbadstrips(ncog,ic)
530    *--------------------------------------------------------------
531    *     this function returns the number of BAD strips
532    *     inside a cluster:
533    *     - if NCOG=0, the number BAD strips inside the whole cluster
534    *     are given, according to the cluster multiplicity
535    *    
536    *     - if NCOG>0, the number BAD strips is evaluated using NCOG
537    *     strips, even if they have a negative signal (according to Landi)
538    *--------------------------------------------------------------
539          include 'commontracker.f'
540          include 'level1.f'
541          include 'calib.f'
542    
543          integer nbad
544          nbad = 0
545    
546          if (ncog.gt.0) then
547    
548    *     --> signal of the central strip
549             sc = CLSIGNAL(INDMAX(ic)) !center
550    *     signal of adjacent strips
551             sl1 = -100000                !left 1
552             if(
553         $        (INDMAX(ic)-1).ge.INDSTART(ic)
554         $        )
555         $        sl1 = CLSIGNAL(INDMAX(ic)-1)
556            
557             sl2 = -100000                !left 2
558             if(
559         $        (INDMAX(ic)-2).ge.INDSTART(ic)
560         $        )
561         $        sl2 = CLSIGNAL(INDMAX(ic)-2)
562            
563             sr1 = -100000                !right 1
564             if(
565         $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
566         $        .or.
567         $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
568         $        )
569         $        sr1 = CLSIGNAL(INDMAX(ic)+1)
570            
571             sr2 = -100000                !right 2
572             if(
573         $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
574         $        .or.
575         $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
576         $        )
577         $        sr2 = CLSIGNAL(INDMAX(ic)+2)
578            
579             if(ncog.ge.1)nbad = nbad+1-CLBAD(INDMAX(ic))
580            
581             if(ncog.ge.2)then
582                if(sl1.gt.sr1.and.(INDMAX(ic)-1).ge.1)then
583                   nbad=nbad+1-CLBAD(INDMAX(ic)-1)      
584                elseif(sl1.le.sr1.and.(INDMAX(ic)+1).le.NCLSTR1)then
585                   nbad=nbad+1-CLBAD(INDMAX(ic)+1)          
586                endif
587             endif
588            
589             if(ncog.ge.3)then
590                if(sl1.gt.sr1.and.(INDMAX(ic)+1).le.NCLSTR1)then
591                   nbad=nbad+1-CLBAD(INDMAX(ic)+1)
592                elseif(sl1.le.sr1.and.(INDMAX(ic)-1).ge.1)then
593    c               if(INDMAX(ic)-1.eq.0)
594    c     $              print*,' ======= ',sl2,sl1,sc,sr1,sr2
595                   nbad=nbad+1-CLBAD(INDMAX(ic)-1)
596                endif
597             endif
598            
599             if(ncog.ge.4)then
600                if(sl2.gt.sr2.and.(INDMAX(ic)-2).ge.1)then
601                   nbad=nbad+1-CLBAD(INDMAX(ic)-2)      
602                elseif(sl2.le.sr2.and.(INDMAX(ic)+2).le.NCLSTR1)then
603                   nbad=nbad+1-CLBAD(INDMAX(ic)+2)          
604                endif
605             endif
606            
607    c         if(ncog.ge.5)then
608    c            print*,'function CLBAD(NCOG,IC) ==> WARNING!! NCOG=',NCOG
609    c     $           ,' not implemented'
610    c         endif
611                
612          elseif(ncog.eq.0)then
613    *     =========================
614    *     COG computation
615    *     =========================
616    
617    
618            
619             iv=VIEW(ic)        
620             if(mod(iv,2).eq.1)incut=incuty
621             if(mod(iv,2).eq.0)incut=incutx
622            
623             istart = INDSTART(IC)
624             istop  = TOTCLLENGTH
625             if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1
626             nbad = 0
627    c$$$         do i=istart,istop
628    c$$$            cut  = incut*CLSIGMA(i)            
629    c$$$            if(CLSIGNAL(i).ge.cut)nbad = nbad +1 -CLBAD(i)
630    c$$$         enddo
631             do i = INDMAX(IC),istart,-1
632                cut  = incut*CLSIGMA(i)
633                if(CLSIGNAL(i).ge.cut)then
634                   nbad = nbad +1 -CLBAD(i)
635                else
636                   goto 10
637                endif
638             enddo
639     10      continue
640             do i = INDMAX(IC)+1,istop
641                cut  = incut*CLSIGMA(i)
642                if(CLSIGNAL(i).ge.cut)then
643                   nbad = nbad +1 -CLBAD(i)
644                else
645                   goto 20
646                endif
647             enddo
648     20      continue
649            
650          else
651            
652    c         print*,'function CLBAD(NCOG,IC) ==> WARNING!! NCOG=',NCOG
653    c     $        ,' not implemented'
654            
655    
656          endif
657    
658          nbadstrips = nbad
659    
660          return
661          end

Legend:
Removed from v.1.1.1.1  
changed lines
  Added in v.1.5

  ViewVC Help
Powered by ViewVC 1.1.23