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

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

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

revision 1.23 by pam-fi, Wed Mar 5 17:00:20 2008 UTC revision 1.24 by pam-fi, Sat Mar 22 08:32:50 2008 UTC
# Line 8  Line 8 
8  *  *
9  *     integer function npfastrips(ic,angle)  *     integer function npfastrips(ic,angle)
10  *  *
11    *     -----------------------------------------------------------------
12    *     p.f.a.
13    *     -----------------------------------------------------------------
14  *     real function pfaeta(ic,angle)  *     real function pfaeta(ic,angle)
15  *     real function pfaetal(ic,angle)  *     real function pfaetal(ic,angle)
16  *     real function pfaeta2(ic,angle)  *     real function pfaeta2(ic,angle)
# Line 15  Line 18 
18  *     real function pfaeta4(ic,angle)  *     real function pfaeta4(ic,angle)
19  *     real function cog(ncog,ic)  *     real function cog(ncog,ic)
20  *  *
21    *     -----------------------------------------------------------------
22    *     risoluzione spaziale media, stimata dalla simulazione (samuele)
23    *     -----------------------------------------------------------------
24    *     FUNCTION risxeta2(angle)
25    *     FUNCTION risxeta3(angle)
26    *     FUNCTION risxeta4(angle)
27    *     FUNCTION risyeta2(angle)
28    *     FUNCTION risy_cog(angle)
29    *     FUNCTION risx_cog(angle)
30    *     real function riseta(iview,angle)
31    *     -----------------------------------------------------------------
32    *     fattore moltiplicativo per tenere conto della dipendenza della
33    *     risoluzione dal rumore delle strip
34    *     -----------------------------------------------------------------
35  *     real function fbad_cog(ncog,ic)  *     real function fbad_cog(ncog,ic)
36  *     real function fbad_eta(ic,angle)  *     real function fbad_eta(ic,angle)
37  *  *
38  *     real function riseta(iview,angle)  *     -----------------------------------------------------------------
39  *     FUNCTION risxeta2(x)  *     NUOVO APPROCCIO PER LA STIMA DELLA RISOLUZIONE
40  *     FUNCTION risxeta3(x)  *     -----------------------------------------------------------------
41  *     FUNCTION risxeta4(x)  *     real function riscogtheor(ncog,ic)
42  *     FUNCTION risyeta2(x)  *     real function risetatheor(ncog,ic,angle)
 *     FUNCTION risy_cog(x)  
 *     FUNCTION risx_cog(x)  
43  *  *
44    *     -----------------------------------------------------------------
45    *     correzione landi
46    *     -----------------------------------------------------------------
47  *     real function pfacorr(ic,angle)  *     real function pfacorr(ic,angle)
48  *  *
49  *     real function effectiveangle(ang,iview,bbb)  *     real function effectiveangle(ang,iview,bbb)
# Line 552  c      if(mod(int(VIEW(ic)),2).eq.1)then Line 570  c      if(mod(int(VIEW(ic)),2).eq.1)then
570  *     resolution.  *     resolution.
571  *     It calls the function FBAD_COG(NCOG,IC),  *     It calls the function FBAD_COG(NCOG,IC),
572  *     accordingto the angle  *     accordingto the angle
573    *
574    *     >>> cosi` non e` corretto!!
575    *     >>> l'errore sulla coordinata eta si ottiene moltiplicando
576    *     >>> l'errore sulla coordinata cog per la derivata della
577    *     >>> distribuzione eta... pur sapendolo l'ho sempre ignorato...
578    *     >>> deve essere modificato!!!!
579    *
580  *-------------------------------------------------------  *-------------------------------------------------------
581    
582        include 'commontracker.f'        include 'commontracker.f'
# Line 589  c      if(mod(int(VIEW(ic)),2).eq.1)then Line 614  c      if(mod(int(VIEW(ic)),2).eq.1)then
614        end        end
615    
616  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
617        real function pfaeta2(ic,angle) !(1)        real function pfaeta2(ic,angle)
618  *--------------------------------------------------------------  *--------------------------------------------------------------
619  *     this function returns  *     this function returns
620  *  *
# Line 608  c      if(mod(int(VIEW(ic)),2).eq.1)then Line 633  c      if(mod(int(VIEW(ic)),2).eq.1)then
633        real cog2,angle        real cog2,angle
634        integer iview,lad        integer iview,lad
635    
636        iview = VIEW(ic)                    iview   = VIEW(ic)            
637        lad = nld(MAXS(ic),VIEW(ic))        lad     = nld(MAXS(ic),VIEW(ic))
638        cog2 = cog(2,ic)                  cog2    = cog(2,ic)          
639        pfaeta2=cog2        pfaeta2 = cog2
640    
641  *     ----------------  *     ----------------
642  *     find angular bin  *     find angular bin
# Line 1311  c            COG = 0. Line 1336  c            COG = 0.
1336        end        end
1337    
1338    
1339  c$$$*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1340  c$$$      real function fbad_cog0(ncog,ic)  
1341  c$$$*-------------------------------------------------------        real function riscogtheor(ncog,ic)
1342  c$$$*     this function returns a factor that takes into  *-------------------------------------------------------
1343  c$$$*     account deterioration of the spatial resolution  *
1344  c$$$*     in the case BAD strips are included in the cluster.  *     this function returns the expected resolution
1345  c$$$*     This factor should multiply the nominal spatial  *     obtained by propagating the strip noise
1346  c$$$*     resolution.  *     to the center-of-gravity coordinate
1347  c$$$*  *
1348  c$$$*     NB!!!  *     ncog = n.strip used in the coordinate evaluation
1349  c$$$*     (this is the old version. It consider only the two  *     (ncog=0 => all strips above threshold)
1350  c$$$*     strips with the greatest signal. The new one is  *
1351  c$$$*     fbad_cog(ncog,ic) )  *-------------------------------------------------------
1352  c$$$*      
1353  c$$$*-------------------------------------------------------        include 'commontracker.f'
1354  c$$$        include 'level1.f'
1355  c$$$      include 'commontracker.f'        include 'calib.f'
1356  c$$$      include 'level1.f'  
1357  c$$$      include 'calib.f'        if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
1358  c$$$           incut = incuty
1359  c$$$*     --> signal of the central strip           pitch = pitchY / 1.e4
1360  c$$$      sc = CLSIGNAL(INDMAX(ic)) !center        else                      !X-view
1361  c$$$           incut = incutx
1362  c$$$*     signal of adjacent strips           pitch = pitchX / 1.e4
1363  c$$$*     --> left        endif
1364  c$$$      sl1 = 0                  !left 1        
1365  c$$$      if(        func = 100000.
1366  c$$$     $     (INDMAX(ic)-1).ge.INDSTART(ic)        stot = 0.
1367  c$$$     $     )  
1368  c$$$     $     sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))        if (ncog.gt.0) then
1369  c$$$  
1370  c$$$      sl2 = 0                  !left 2  *     --> signal of the central strip
1371  c$$$      if(           sc = CLSIGNAL(INDMAX(ic)) !center
1372  c$$$     $     (INDMAX(ic)-2).ge.INDSTART(ic)           fsc = clsigma(INDMAX(ic))
1373  c$$$     $     )  *     --> signal of adjacent strips
1374  c$$$     $     sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))           sl1 = 0                !left 1
1375  c$$$           fsl1 = 1               !left 1
1376  c$$$*     --> right           if(
1377  c$$$      sr1 = 0                  !right 1       $        (INDMAX(ic)-1).ge.INDSTART(ic)
1378  c$$$      if(       $        )then
1379  c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))              sl1 = CLSIGNAL(INDMAX(ic)-1)
1380  c$$$     $     .or.              fsl1 = clsigma(INDMAX(ic)-1)
1381  c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)           endif
1382  c$$$     $     )  
1383  c$$$     $     sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))           sl2 = 0                !left 2
1384  c$$$           fsl2 = 1               !left 2
1385  c$$$      sr2 = 0                  !right 2           if(
1386  c$$$      if(       $        (INDMAX(ic)-2).ge.INDSTART(ic)
1387  c$$$     $     (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))       $        )then
1388  c$$$     $     .or.              sl2 = CLSIGNAL(INDMAX(ic)-2)
1389  c$$$     $     (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)              fsl2 = clsigma(INDMAX(ic)-2)
1390  c$$$     $     )           endif
1391  c$$$     $     sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))           sr1 = 0                !right 1
1392  c$$$           fsr1 = 1               !right 1
1393  c$$$           if(
1394  c$$$      if(mod(int(VIEW(ic)),2).eq.1)then !Y-view       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
1395  c$$$         f  = 4.       $        .or.
1396  c$$$         si = 8.4       $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
1397  c$$$      else                              !X-view       $        )then
1398  c$$$         f  = 6.              sr1 = CLSIGNAL(INDMAX(ic)+1)
1399  c$$$         si = 3.9              fsr1 = clsigma(INDMAX(ic)+1)
1400  c$$$      endif           endif    
1401  c$$$           sr2 = 0                !right 2
1402  c$$$      fbad_cog = 1.           fsr2 = 1               !right 2
1403  c$$$      f0 = 1           if(
1404  c$$$      f1 = 1       $        (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
1405  c$$$      f2 = 1       $        .or.
1406  c$$$      f3 = 1         $        (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
1407  c$$$      if(sl1.gt.sr1.and.sl1.gt.0.)then       $        )then
1408  c$$$                      sr2 = CLSIGNAL(INDMAX(ic)+2)
1409  c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f              fsr2 = clsigma(INDMAX(ic)+2)
1410  c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f1=f           endif
1411  c$$$c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f3=f  
1412  c$$$  
1413  c$$$         if(ncog.eq.2.and.sl1.ne.0)then  
1414  c$$$            fbad_cog = (f1**2*sc**2/sl1**2+f0**2)/(sc**2/sl1**2+1.)  ************************************************************
1415  c$$$         elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then  *     COG2-3-4 computation
1416  c$$$            fbad_cog = 1.  ************************************************************
1417  c$$$         elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then  
1418  c$$$            fbad_cog = 1.  c      print*,sl2,sl1,sc,sr1,sr2
1419  c$$$         else          
1420  c$$$            fbad_cog = 1.           vCOG = cog(ncog,ic)!0.
1421  c$$$         endif          
1422  c$$$                   if(ncog.eq.1)then
1423  c$$$      elseif(sl1.le.sr1.and.sr1.gt.0.)then              func = 1./12.
1424  c$$$              stot = 1.
1425  c$$$           elseif(ncog.eq.2)then
1426  c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic))  ).eq.0)f0=f              if(sl1.gt.sr1)then
1427  c$$$         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f1=f                 func = (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2)
1428  c$$$c         if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f3=f                 stot = sl1+sc
1429  c$$$              elseif(sl1.le.sr1)then
1430  c$$$         if(ncog.eq.2.and.sr1.ne.0)then                 func = (fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1431  c$$$            fbad_cog = (f1**2*sc**2/sr1**2+f0**2)/(sc**2/sr1**2+1.)                 stot = sc+sr1
1432  c$$$         elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then              endif
1433  c$$$            fbad_cog = 1.           elseif(ncog.eq.3)then
1434  c$$$         elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then              func =
1435  c$$$            fbad_cog = 1.       $           (fsl1*(-1-vCOG)**2+fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1436  c$$$         else              stot = sl1+sc+sr1
1437  c$$$            fbad_cog = 1.           elseif(ncog.eq.4)then
1438  c$$$         endif              if(sl2.gt.sr2)then
1439  c$$$                 func =
1440  c$$$      endif       $              (fsl2*(-2-vCOG)**2+fsl1*(-1-vCOG)**2
1441  c$$$       $              +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2)
1442  c$$$      fbad_cog0 = sqrt(fbad_cog)                 stot = sl2+sl1+sc+sr1
1443  c$$$              elseif(sl2.le.sr2)then
1444  c$$$      return                 func =
1445  c$$$      end       $              (fsl1*(-1-vCOG)**2
1446  c$$$       $              +fsc*(-vCOG)**2+fsr1*(1-vCOG)**2+fsr2*(2-vCOG)**2)
1447  c$$$                 stot = sl2+sl1+sc+sr1
1448  c$$$              endif
1449             else
1450                print*,'function riscogtheor(NCOG,IC) ==> NCOG=',NCOG
1451         $            ,' not implemented'
1452             endif
1453            
1454          elseif(ncog.eq.0)then
1455    *     =========================
1456    *     COG computation
1457    *     =========================
1458            
1459             vCOG = cog(0,ic)
1460    
1461             iv     = VIEW(ic)
1462             istart = INDSTART(IC)
1463             istop  = TOTCLLENGTH
1464             if(ic.lt.NCLSTR1)istop = INDSTART(IC+1)-1
1465    ccc         SGN = 0.
1466             SNU = 0.
1467    ccc         SDE = 0.
1468    
1469             do i=INDMAX(IC),istart,-1
1470                ipos = i-INDMAX(ic)
1471                cut  = incut*CLSIGMA(i)
1472                if(CLSIGNAL(i).gt.cut)then
1473                   fs = clsigma(i)
1474                   SNU  = SNU + fs*(ipos-vCOG)**2
1475                   stot = stot + CLSIGNAL(i)
1476                else
1477                   goto 10
1478                endif            
1479             enddo
1480     10      continue
1481             do i=INDMAX(IC)+1,istop
1482                ipos = i-INDMAX(ic)
1483                cut  = incut*CLSIGMA(i)
1484                if(CLSIGNAL(i).gt.cut)then
1485                   fs = clsigma(i)
1486                   SNU  = SNU + fs*(ipos-vCOG)**2
1487                   stot = stot + CLSIGNAL(i)
1488                else
1489                   goto 20
1490                endif            
1491             enddo
1492     20      continue
1493             if(SDE.ne.0)then
1494                FUNC=SNU
1495             else
1496                
1497             endif
1498    
1499          else
1500                      
1501             FUNC=0
1502             print*,'function FUNC(NCOG,IC) ==> WARNING!! NCOG=',NCOG
1503             print*,'                               (NCOG must be >= 0)'
1504            
1505    
1506          endif
1507    
1508    
1509          if(stot.gt.0..and.func.gt.0.)then
1510             func = sqrt(func)
1511             func = pitch * func/stot
1512          endif
1513    
1514          riscogtheor = func
1515    
1516          return
1517          end
1518    
1519    
1520    
1521    *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1522    
1523          real function risetatheor(ncog,ic,angle)
1524    *-------------------------------------------------------
1525    *
1526    *     this function returns the expected resolution
1527    *     obtained by propagating the strip noise
1528    *     to the coordinate evaluated with non-linear eta-function
1529    *
1530    *     ncog = n.strip used in the coordinate evaluation
1531    *     (ncog=0 => ncog=2,3,4 according to angle)
1532    *
1533    *-------------------------------------------------------
1534    
1535          include 'commontracker.f'
1536          include 'level1.f'
1537          include 'calib.f'
1538    
1539    
1540          func = 1.
1541    
1542          iview   = VIEW(ic)            
1543          lad     = nld(MAXS(ic),VIEW(ic))
1544    
1545    *     ------------------------------------------------
1546    *     number of strip to be used (in case of ncog = 0)
1547    *     ------------------------------------------------
1548    
1549          inoeta = 0
1550    
1551          if(mod(int(iview),2).eq.1)then !Y-view
1552    
1553             pitch = pitchY / 1.e4
1554    
1555             if(ncog.eq.0)then
1556                if( abs(angle).ge.e2fay.and.abs(angle).le.e2tay )then
1557                   ncog = 2
1558                elseif( abs(angle).ge.e3fay.and.abs(angle).le.e3tay )then
1559                   ncog = 3
1560                elseif( abs(angle).ge.e4fay.and.abs(angle).le.e4tay )then
1561                   ncog = 4
1562                else
1563                   ncog   = 4
1564                   inoeta = 1
1565                endif            
1566             endif
1567    
1568          else                      !X-view
1569    
1570             pitch = pitchX / 1.e4
1571    
1572             if(ncog.eq.0)then
1573                if( abs(angle).ge.e2fax.and.abs(angle).le.e2tax )then
1574                   ncog = 2
1575                elseif( abs(angle).ge.e3fax.and.abs(angle).le.e3tax )then
1576                   ncog = 3
1577                elseif( abs(angle).ge.e4fax.and.abs(angle).le.e4tax )then
1578                   ncog = 4
1579                else
1580                   ncog = 4
1581                   inoeta = 1
1582                endif            
1583             endif
1584    
1585          endif
1586          
1587          func = riscogtheor(ncog,ic)
1588    
1589          risetatheor = func
1590          
1591          if(inoeta.eq.1)return ! no eta correction is applied --> exit
1592          if(ncog.lt.1.or.ncog.gt.4)return
1593    
1594    *     ----------------
1595    *     find angular bin
1596    *     ----------------
1597    *     (in futuro possiamo pensare di interpolare anche sull'angolo)
1598          do iang=1,nangbin
1599             if(angL(iang).lt.angle.and.angR(iang).ge.angle)then
1600                iangle=iang
1601                goto 98
1602             endif
1603          enddo
1604          if(DEBUG.EQ.1)print*
1605         $     ,'risetatheor *** warning *** angle out of range: ',angle
1606          if(angle.le.angL(1))iang=1
1607          if(angle.ge.angR(nangbin))iang=nangbin
1608     98   continue                  !jump here if ok
1609    
1610    *     -------------
1611    *     within +/-0.5
1612    *     -------------
1613    
1614          vcog = cog(ncog,ic)          
1615    
1616          etamin = eta2(1,iang)
1617          etamax = eta2(netaval,iang)
1618    
1619          iaddmax=10
1620          iadd=0
1621     10   continue
1622          if(vcog.lt.etamin)then
1623             vcog = vcog + 1
1624             iadd = iadd + 1
1625             if(iadd>iaddmax)goto 111
1626             goto 10
1627          endif
1628     20   continue
1629          if(vcog.gt.etamax)then
1630             vcog = vcog - 1
1631             iadd = iadd - 1
1632             if(iadd<-1*iaddmax)goto 111
1633             goto 20
1634          endif
1635          goto 1111
1636     111  continue
1637          if(DEBUG.eq.1)
1638         $     print*,'risetatheor *** warning *** anomalous cluster'
1639          if(DEBUG.eq.1)
1640         $     print*,'--> COG(',ncog,') = ',vcog-iadd,' (set to zero)'
1641          vcog=0
1642     1111 continue
1643    
1644    *     ------------------------------------------------
1645    *     interpolation
1646    *     ------------------------------------------------
1647    
1648    
1649          ibin = netaval
1650          do i=2,netaval        
1651             if(ncog.eq.2)eta=eta2(i,iang)
1652             if(ncog.eq.3)eta=eta3(i,iang)
1653             if(ncog.eq.4)eta=eta4(i,iang)
1654             if(eta.ge.vcog)then
1655                ibin = i
1656                goto 99
1657             endif
1658          enddo
1659     99   continue
1660    
1661          if(ncog.eq.2)then
1662             x1 = eta2(ibin-1,iang)
1663             x2 = eta2(ibin,iang)
1664             y1 = feta2(ibin-1,iview,lad,iang)
1665             y2 = feta2(ibin,iview,lad,iang)
1666          elseif(ncog.eq.3)then
1667             x1 = eta3(ibin-1,iang)
1668             x2 = eta3(ibin,iang)
1669             y1 = feta3(ibin-1,iview,lad,iang)
1670             y2 = feta3(ibin,iview,lad,iang)
1671          elseif(ncog.eq.4)then
1672             x1 = eta4(ibin-1,iang)
1673             x2 = eta4(ibin,iang)
1674             y1 = feta4(ibin-1,iview,lad,iang)
1675             y2 = feta4(ibin,iview,lad,iang)
1676          endif
1677          
1678          func = func * (y2-y1)/(x2-x1)
1679    
1680          risetatheor = func
1681    
1682          return
1683          end
1684    
1685  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***  *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
1686    
# Line 1948  c$$$ Line 2208  c$$$
2208    
2209        pfacorr = fcorr(iview,lad,iang)        pfacorr = fcorr(iview,lad,iang)
2210    
2211        if(DEBUG.eq.1)print*,'CORR  (ic ',ic,' ang',angle,') -->',pfacorr        if(DEBUG.eq.1)print*,'LANDI (ic ',ic,' ang',angle,') -->',pfacorr
   
2212    
2213          
2214   100  return   100  return
2215        end        end
Legend:
Removed from v.1.23  
changed lines
  Added in v.1.24
  ViewVC Help
Powered by ViewVC 1.1.23